陳軼嵩，陳子童，王 童，劉永濤，羅 耿

（長安大學(xué)汽車學(xué)院，西安710064）

前言

我國綜合國力不斷提升，人均可支配收入逐年攀升，汽車保有量持續(xù)增長，逐漸成為家庭生活必需品，極大促進(jìn)了汽車工業(yè)發(fā)展［1］。但與此同時城市交通狀況也越發(fā)擁堵，大氣環(huán)境污染進(jìn)一步惡化，要求我們大力發(fā)展公共交通［2］。隨著公路客運(yùn)蓬勃發(fā)展，城市公交車及客運(yùn)大巴是人民短途出行優(yōu)選的交通工具［3］；但一些惡劣客車交通事故也時有發(fā)生。由于客車載客量較大，一次事故動輒造成大量人員傷亡，帶來巨大損失，其中客車側(cè)翻死亡率最高，危害極大。如何提高客車側(cè)翻安全性，保護(hù)駕駛員及乘員生命安全，已成為世界各國研究人員亟待解決的迫切任務(wù)［4］。

客車側(cè)翻一步碰撞算法是作者參考板料沖壓一步成型算法思想［5］，提出的一種用于客車側(cè)翻碰撞模擬分析的新算法。對于板料沖壓一步成型算法，Batoz等［6］在1989年開始提出“逆方法”，對沖壓件厚度及應(yīng)力分布進(jìn)行了預(yù)測；柳玉起等［7］對應(yīng)用一步逆有限元法模擬沖壓成形過程中的褶皺、回彈及失穩(wěn)問題進(jìn)行了研究；陳偉等［8］又在金屬板料沖壓一步逆成形算法的基礎(chǔ)上，初步探討了更具工程應(yīng)用前景的“一步正”板料成形模擬算法，快速進(jìn)行金屬板料沖壓成形性能評價，并收到了良好效果。借鑒該算法核心思想，客車側(cè)翻一步碰撞算法同樣基于非線性全量理論，利用了側(cè)翻碰撞過程中能量轉(zhuǎn)換關(guān)系，與現(xiàn)有LS?DYNA等增量法軟件相比，在略犧牲一些計算精度情況下，大幅提升了模擬效率。

現(xiàn)有初始算法是對基于全量理論的側(cè)翻一步碰撞理論體系的初步探討，計算精度和效率仍有一定的改進(jìn)空間。為進(jìn)一步提高算法計算效率且保證精度，使算法能更好預(yù)測客車結(jié)構(gòu)側(cè)翻安全性，為后續(xù)對客車車身結(jié)構(gòu)進(jìn)行基于側(cè)翻安全的拓?fù)鋬?yōu)化提供支撐，本文中在客車側(cè)翻一步碰撞初始算法基礎(chǔ)上，分析算法運(yùn)行過程中的關(guān)鍵環(huán)節(jié)，從計算量較大的積分和迭代部分入手，對初始算法積分和迭代方法進(jìn)行探索和改進(jìn)。首先，對算法單元積分方法，在單元厚度方向（下簡稱厚向）采用不同積分點數(shù)量并引入面內(nèi)減縮積分和局部減縮積分，對厚向及面內(nèi)積分點模式進(jìn)行合理匹配，尋找算法能量和失衡力平衡迭代計算過程中最優(yōu)積分方案。其次，對廣義失衡力平衡迭代，采用更優(yōu)的力與力矩分離迭代法，確定合理權(quán)重比，對力與力矩進(jìn)行加權(quán)，提高迭代速度，完成迭代優(yōu)化。改進(jìn)后算法可在車身設(shè)計初期，對結(jié)構(gòu)側(cè)翻碰撞安全性進(jìn)行快速評價，縮短產(chǎn)品開發(fā)周期，同時，可為后續(xù)針對客車側(cè)翻安全性的靈敏度分析、參數(shù)優(yōu)化和拓?fù)鋬?yōu)化算法研究提供必要的支撐條件。

1 側(cè)翻一步碰撞算法基本理論

客車側(cè)翻一步碰撞算法，基于非線性全量理論和比例加載假定，依據(jù)ECE R66法規(guī)［9-10］，忽略中間狀態(tài)和構(gòu)形變化［11］，只考慮結(jié)構(gòu)碰撞開始和最大變形兩個狀態(tài)。根據(jù)車體側(cè)翻碰撞過程運(yùn)動變形特點和能量轉(zhuǎn)換關(guān)系，得到滿足變形條件初始解，采用Newton?Raphson法迭代求解［12］，快速獲得結(jié)構(gòu)最終變形。

將碰撞開始狀態(tài)結(jié)構(gòu)作為初始構(gòu)形{X0}。此時車體未發(fā)生變形，結(jié)構(gòu)動能Ed為

式中：M為車體質(zhì)量；Δh為車體質(zhì)心下降高度；J為車體繞側(cè)翻假定轉(zhuǎn)軸的轉(zhuǎn)動慣量；ω為車體角速度。

碰撞開始狀態(tài)結(jié)構(gòu)各節(jié)點速度{v0}由式（2）計算：

式中：ri為各節(jié)點到側(cè)翻假定轉(zhuǎn)軸距離；n為節(jié)點數(shù)。

在結(jié)構(gòu)最大變形狀態(tài)，車體結(jié)構(gòu)產(chǎn)生明顯變形。側(cè)翻一步碰撞算法中，結(jié)構(gòu)最大變形不確定須假定一個最大變形構(gòu)形{x0}。各節(jié)點位移{U0}由式（3）計算：

車體動能Ed在碰撞中主要轉(zhuǎn)換為結(jié)構(gòu)形變能，結(jié)構(gòu)形變能W由式（4）計算［13］：

式中：{ε}為單元塑性應(yīng)變；{σ}為單元塑性應(yīng)力；Ve為單元體積；N為單元數(shù)。

結(jié)構(gòu)形變能W與車體動能Ed理應(yīng)相等，即

若式（5）不滿足，則對節(jié)點位移{U0}進(jìn)行修正，按照式（4）重新計算結(jié)構(gòu)形變能。將滿足式（5）能量關(guān)系假定的節(jié)點位移{U}作為Newton?Raphson迭代初始解。

對滿足能量轉(zhuǎn)換關(guān)系初始解{U}，節(jié)點失衡力{R(U)}已處于不平衡狀態(tài)：

應(yīng)用Newton?Raphson法，解決節(jié)點失衡力不平衡問題。對初始解{U}按式（7）迭代求解，使式（6）達(dá)到平衡，得到結(jié)構(gòu)最終變形。

2 單元積分點模式優(yōu)化

側(cè)翻一步碰撞初始算法采用Mindlin殼單元厚向五點且面內(nèi)全積分方案，為在基本保證計算精度基礎(chǔ)上進(jìn)一步提高計算效率，本文在積分方案上做進(jìn)一步研究。單元面內(nèi)，在四點全積分基礎(chǔ)上，引入減縮積分和局部減縮積分方案；單元厚向，引入五點、七點和九點3種高斯積分方案；并在厚向各積分層采用減縮積分與局部減縮積分的不同匹配方案。經(jīng)多次研究對比發(fā)現(xiàn)，薄殼單元中性層積分點數(shù)量對模擬結(jié)果影響較大，且中性層局部減縮積分方案具有更高的計算精度。因此，為保證一定精度，本文中所有積分方案均采用在中性層面內(nèi)應(yīng)用局部減縮積分方式進(jìn)行匹配，同時為保證計算效率提升效果，每個積分方案中其余積分層最多引入兩層局部減縮積分。

2.1 單元面內(nèi)積分方案

面內(nèi)減縮積分即單點高斯積分，計算效率較高，可避免使用全積分導(dǎo)致剛度過大及剪切鎖死問題，但數(shù)學(xué)上看，計算精度相對較低，其積分點位置和積分系數(shù)為

式中：r、s為積分點坐標(biāo)；h為積分系數(shù)；i=1，j=1。

本文中所提局部減縮積分是一種復(fù)合式積分方法，即利用單點積分計算內(nèi)力中切應(yīng)力相關(guān)分量，利用四點積分計算正應(yīng)力相關(guān)分量，既可克服單點積分引起的沙漏模式，又可避免四點積分引起的剪切鎖死。

為實現(xiàn)局部減縮積分，將單元［B］矩陣和應(yīng)力向量{σ}進(jìn)行分解，令：

其中

單元內(nèi)力為

令：

式中：{fintσ}為只與正應(yīng)力σx和σy的相關(guān)項，數(shù)值求解時，為提高計算精度，沿r、s方向均取2點積分；{fintτ}為只與切應(yīng)力τxy的相關(guān)項，為避免剪切鎖死，求解時采用單點減縮積分，即

2.2 單元厚向積分方案

2.2.1 厚向五點高斯積分

當(dāng)厚向采用五點高斯積分時，各單元厚向積分點位置如圖1所示。1點到5點的高斯積分點位置因子分別為0.906 179 85、0.538 469 31、0、-0.538 469 31、-0.906 179 85，對應(yīng)的積分系數(shù)分別 為0.236 926 89、0.478 628 67、0.568 888 89、0.478 628 67、0.236 926 89。

圖1 單元厚向五點積分示意圖

同時，考慮到積分方案對中性層對稱性設(shè)計，單元厚向五點積分各積分層面內(nèi)減縮積分與局部減縮積分匹配方案設(shè)計如表1所示。

表1 厚向五點高斯積分方案

以某客車車身段模型作為分析對象，將厚向五點各積分方案應(yīng)用于客車側(cè)翻一步碰撞快速算法進(jìn)行計算，獲得車身段兩側(cè)立柱各測點變形量和模擬時間，并與LS?DYNA計算結(jié)果對比，車身段封閉環(huán)兩側(cè)立柱各測點變形量對比的柱狀圖如圖2所示。

圖2 厚向五點高斯積分與LS?DYNA測點變形量對比

接著，對應(yīng)用厚向五點3個積分方案分析計算時間與LS?DYNA時間進(jìn)行對比，結(jié)果如表2所示。

表2 厚向五點高斯積分與LS?DYNA分析效率對比

本文中選取11個測點變形量作為評價指標(biāo)，對算法計算精度進(jìn)行考量。通過對比分析發(fā)現(xiàn)，計算精度上3種積分方案計算結(jié)果十分接近，與LS?DYNA計算結(jié)果平均相對誤差分別為8.895%、8.996%和9.102%，方案三計算精度略低。計算效率上，3個積分方案模擬時間比LS?DYNA都有明顯縮短，且方案三計算效率比前兩個方案略有提高。

2.2.2 厚向七點高斯積分

當(dāng)單元厚向采用七點積分時，隨著單元厚向積分點數(shù)的增加，其積分方案的匹配方式也有所增加，共設(shè)計4種積分方案。單元厚向七點積分各積分層面內(nèi)減縮積分與局部減縮積分的匹配方案設(shè)計如表3所示。

表3 厚向七點高斯積分方案

厚向七點積分與圖1中所示的單元厚向五點積分示意圖類似，只是積分點數(shù)量和位置有所不同。厚向七點積分中積分點1到積分點7的高斯積分點位置因子分別為0.949 107 91、0.741 531 19、0.405 845 15、0、-0.405 845 15、-0.741 531 19、-0.949 107 91，對應(yīng)的積分系數(shù)分別為0.101 228 54、0.279 705 39、0.381 830 05、0.417 959 18、0.381 830 05、0.279 705 39、0.101 228 54。

同樣，將厚向七點各積分方案應(yīng)用于客車側(cè)翻一步碰撞快速算法進(jìn)行改進(jìn)，獲得封閉環(huán)兩側(cè)立柱各測點變形量及計算時間，對比結(jié)果如圖3所示。

圖3 厚向七點高斯積分與LS?DYNA測點變形量對比

接著，將厚向七點4個積分方案計算時間與LS?DYNA分析所用時間進(jìn)行對比，確定各方案計算效率優(yōu)劣，結(jié)果如表4所示。

表4 厚向七點高斯積分與LS?DYNA分析效率對比

同樣，以11個測點變形量平均相對誤差作為每個積分方案算法計算精度的評價標(biāo)準(zhǔn)。通過誤差計算和對比發(fā)現(xiàn)，在計算精度上厚向七點4種積分方案計算結(jié)果與LS?DYNA平均相對誤差分別為3.278%、3.675%、3.667%、4.087%，方案一計算精度最高，方案二和方案三計算精度相差不多，方案四計算精度比前3個方案有所降低。在計算效率上，4個方案模擬時間相差不多，比LS?DYNA仍有明顯縮短。與厚向五點積分相比，計算精度提高較大，計算效率略微降低。

2.2.3 厚向九點高斯積分

當(dāng)單元厚向采用九點積分時，隨著單元厚向積分點數(shù)進(jìn)一步增加，其積分方案匹配方式繼續(xù)增加，共設(shè)計5種單元積分方案。單元厚向九點積分各積分層面內(nèi)減縮積分與局部減縮積分匹配方案設(shè)計如表5所示。

表5 厚向九點高斯積分方案

厚向九點積分與圖1中所示的單元厚向五點積分示意圖類似，只是積分點數(shù)量達(dá)到9個。厚向九點積分中積分點1到積分點9的高斯積分點位置因子分別為0.968 160 24、0.836 031 11、0.613 371 43、0.324 253 42、0、-0.324 253 42、-0.613 371 43、-0.836 031 11、-0.968 160 24，對應(yīng)的積分系數(shù)分別 為0.081 274 39、0.180 648 16、0.260 610 69、0.312 347 08、0.330 239 36、0.312 347 08、0.260 610 69、0.180 648 16、0.081 274 39。

同樣，將厚向九點各積分方案應(yīng)用于客車側(cè)翻一步碰撞快速算法進(jìn)行改進(jìn)，獲得封閉環(huán)兩側(cè)立柱各測點變形量和計算時間。統(tǒng)計各方案中11個測點變形量數(shù)據(jù)，將厚向九點各積分方案計算結(jié)果與LS?DYNA結(jié)果進(jìn)行對比，封閉環(huán)兩側(cè)立柱各測點變形量對比柱狀圖如圖4所示。

圖4 厚向九點高斯積分與LS?DYNA測點變形量對比

接著，將厚向九點5個積分方案計算時間與LS?DYNA進(jìn)行對比，確定各方案計算效率優(yōu)劣，結(jié)果如表6所示。

表6 厚向九點高斯積分與LS?DYNA分析效率對比

同樣，以11個測點變形量平均相對誤差作為每個積分方案算法計算精度評價標(biāo)準(zhǔn)。通過誤差計算和對比分析發(fā)現(xiàn)，在計算精度上厚向九點5種方案計算結(jié)果與LS?DYNA平均相對誤差分別為2.868%、2.927%、3.032%、3.038%、3.210%，方案一計算精度稍好，方案五計算精度比前4個方案略低。在計算效率上，方案五計算效率比前4個方案略有提高，比LS?DYNA仍有明顯縮短。對比五點、七點積分來說效率下降較為明顯，且比厚向七點積分精度提高并不明顯。

綜合對比以上結(jié)果可知，隨著厚向積分點個數(shù)增加，計算精度逐步提高，但計算效率相應(yīng)下降。從上述分析可看出，厚向七點積分方案一具有較高計算精度，滿足算法精度要求，同時由于減縮積分和局部減縮積分引入，使其具有較高計算效率。綜合考慮精度與效率，本文最終確定厚向七點積分中方案一作為初始算法改進(jìn)積分方案。

3 廣義失衡力迭代模式優(yōu)化

客車側(cè)翻碰撞變形復(fù)雜，單元除有面內(nèi)拉壓引起的面內(nèi)變形，還有彎曲變形。側(cè)翻一步碰撞初始算法迭代過程采用力和力矩融合迭代方式，沒有引入比例系數(shù)進(jìn)行調(diào)整，該迭代方式存在不合理性，忽略了兩點：第一，力和力矩量綱不同，數(shù)值差別也較大，直接融合在一起進(jìn)行計算失衡力“2”范數(shù)，會影響廣義失衡力迭代收斂效率；第二，力和力矩對廣義失衡力平衡迭代收斂效率貢獻(xiàn)度也有所不同，不考慮權(quán)重直接進(jìn)行迭代，同樣會影響計算效率。故本文將力和力矩進(jìn)行分離，分別計算各自失衡力“2”范數(shù)，取不同權(quán)重比對其進(jìn)行加權(quán)調(diào)整，對側(cè)翻一步碰撞算法進(jìn)行改進(jìn)。

首先，將計算得到的Mindlin組合殼單元節(jié)點廣義內(nèi)力，由式（19）表示：

接著，計算各節(jié)點廣義失衡力向量R(ui)，如式（20）所示：

其中每個節(jié)點共包含3個失衡力Rix、Riy、Riz和3個失衡力矩Rixx、Riyy、Rizz，即

節(jié)點廣義失衡力計算完成后，將每個節(jié)點中的失衡力分離出來單獨(dú)計算失衡力“2”范數(shù)：

式中n為車體結(jié)構(gòu)中單元節(jié)點總數(shù)。

同樣，將每個節(jié)點中失衡力矩分離出來單獨(dú)計算力矩“2”范數(shù)：

取一組權(quán)系數(shù)將上述力和力矩“2”范數(shù)進(jìn)行整合，得到加權(quán)調(diào)整后的廣義失衡力“2”范數(shù)，即

式中α、β分別為力和力矩的權(quán)系數(shù)，其和為2。

利用算得的結(jié)構(gòu)變形情況初步判斷計算結(jié)果能否滿足工程計算的合理要求，統(tǒng)計不同權(quán)重比下算法的計算時間，繪制計算時間-權(quán)重比重線，如圖5所示。通過選取曲線中計算時間較短的區(qū)段進(jìn)行多項式擬合，并求得擬合曲線的極值確定合理權(quán)重比，以確定改進(jìn)算法中力和力矩的合理權(quán)系數(shù)值。

圖5 計算時間-權(quán)重比曲線

算法計算時間曲線開始時有一個較短的下降段，隨后則逐漸上升，計算時間最短的區(qū)間段在權(quán)重比為0.11~1之間。通過對比計算，可確定計算時間最短的點所對應(yīng)的權(quán)重比為0.479 2，故最后得到合理的權(quán)重值α、β分別為0.684和1.352［1］。

選取合適松弛因子，利用式（7）中迭代過程對廣義失衡力進(jìn)行迭代，在每步迭代中，均需對比式（20）中所示廣義失衡力“2”范數(shù)的大小，直到達(dá)到最小值時結(jié)束迭代，即

此時，車體結(jié)構(gòu)處于穩(wěn)定狀態(tài)，該狀態(tài)下車體結(jié)構(gòu)變形即為客車側(cè)翻碰撞后的結(jié)構(gòu)最終變形。

4 應(yīng)用實例

為檢驗單元積分點模式和廣義失衡力迭代模式優(yōu)化改進(jìn)后側(cè)翻一步碰撞改進(jìn)算法的有效性，本文中選取某12 m公路客車車身結(jié)構(gòu)典型車身段作為研究對象。利用改進(jìn)算法對該車身段進(jìn)行側(cè)翻碰撞分析，將計算結(jié)果與初始算法、LS?DYNA和側(cè)翻試驗結(jié)果進(jìn)行精度及效率對比。

4.1 模型建立

車身段側(cè)翻試驗參考?xì)W洲ECE R66法規(guī)設(shè)計，試驗臺高800 mm，試驗?zāi)Ｐ团渲睾筚|(zhì)量2 409 kg，質(zhì)心高度1.5 m，典型車身段試驗?zāi)Ｐ腿鐖D6所示。

圖6 典型車身段模型示意圖

對該車身段試驗?zāi)Ｐ瓦M(jìn)行簡化并假設(shè)，建立車身段有限元模型如圖7所示。模型共離散270 171個單元，269 894個節(jié)點，骨架鋼材選用Q345鋼，材料彈性模量E=2.06×1011Pa，泊松比μ=0.3，密度ρ=7800 kg/m3，屈服強(qiáng)度σs=345 MPa。

圖7 車身段有限元模型

車身段側(cè)翻試驗?zāi)Ｐ鸵约跋鄳?yīng)的有限元模型中所包含的A、B兩個封閉環(huán)，參考企業(yè)側(cè)翻試驗過程中測點選取方式，在封閉環(huán)兩側(cè)立柱共選取11個測點，左側(cè)立柱內(nèi)表面選取6個測點，間隔100 mm，右側(cè)立柱內(nèi)表面選取5個測點，間隔100 mm，各測點具體位置見圖8。

圖8 封閉環(huán)模型測點選取方式示意圖

將有限元模型導(dǎo)入LS?DYNA求解計算，兩側(cè)立柱各測點變形量柱狀圖如圖9所示。

圖9 封閉環(huán)模型LS?DYNA分析測點變形量柱狀圖

4.2 結(jié)果對比

初始算法結(jié)構(gòu)變形情況如圖10所示，改進(jìn)算法結(jié)構(gòu)變形情況如圖11所示，LS?DYNA分析結(jié)構(gòu)變形情況如圖12所示。

圖10 初始算法結(jié)構(gòu)變形情況

圖12 LS?DYNA結(jié)構(gòu)變形情況

由圖10和圖11可看出，兩種算法結(jié)構(gòu)變形趨勢一致性良好，改進(jìn)算法具有基本工程計算合理性。為對改進(jìn)算法計算精度進(jìn)行更準(zhǔn)確分析，需對各測點變形量進(jìn)行定量對比。

圖11 改進(jìn)算法結(jié)構(gòu)變形情況

側(cè)翻試驗測點變形量數(shù)據(jù)來自企業(yè)實車開發(fā)過程中車身段側(cè)翻試驗，數(shù)據(jù)由企業(yè)技術(shù)人員協(xié)助測得。側(cè)翻試驗、LS?DYNA分析、初始算法和改進(jìn)算法所得封閉環(huán)A、B兩側(cè)立柱變形量數(shù)據(jù)統(tǒng)計如表6所示。

表6 典型車身段測點變形量統(tǒng)計 mm

封閉環(huán)A兩側(cè)立柱各測點變形量對比見圖13，封閉環(huán)B兩側(cè)立柱各測點變形量對比如圖14所示。

圖13 封閉環(huán)A各測點變形量對比

圖14 封閉環(huán)B各測點變形量對比

通過對比分析發(fā)現(xiàn)，改進(jìn)算法與LS?DYNA平均誤差為2.291%，初始算法與LS?DYNA平均誤差為3.37%；改進(jìn)算法與側(cè)翻試驗平均誤差為10.349%，初始算法與側(cè)翻試驗的平均誤差為10.765%?？梢哉J(rèn)為，改進(jìn)算法在計算精度上有所提升，且小于工程有限元仿真分析誤差容許值15%，基本滿足算法所需計算精度。

接著，對計算效率進(jìn)行對比，結(jié)果如表7所示。

表7 計算效率對比

可以看出，改進(jìn)算法計算效率在基本保證精度的情況下明顯提高，檢驗了本文所提改進(jìn)方法在實際應(yīng)用中的有效性。

5 結(jié)論

為進(jìn)一步提高側(cè)翻一步碰撞算法計算效率且保證計算精度，為后續(xù)對客車車身結(jié)構(gòu)進(jìn)行基于側(cè)翻安全的拓?fù)鋬?yōu)化提供支撐，在客車側(cè)翻一步碰撞初始算法基礎(chǔ)上，分析算法運(yùn)行過程中的關(guān)鍵環(huán)節(jié)，從計算量較大的積分及迭代部分入手，對初始算法積分及迭代方法進(jìn)行探索和改進(jìn)。經(jīng)過計算分析，改進(jìn)算法在基本保證計算精度的同時，計算效率得到顯著提升。本文中所提單元積分點模式及廣義失衡力迭代模式改進(jìn)方案，具有良好應(yīng)用效果。

但是，初始算法本身單元模型基于Mindlin殼單元，受單元模型理論限制，算法計算精度不能大幅提升，后續(xù)仍需對算法單元模型進(jìn)行改進(jìn)研究，進(jìn)一步提升算法計算精度和效率，以便在車身設(shè)計初期，對結(jié)構(gòu)側(cè)翻碰撞安全性進(jìn)行更加快速準(zhǔn)確評價，進(jìn)一步縮短產(chǎn)品開發(fā)周期，為后續(xù)針對客車側(cè)翻安全性的靈敏度分析、參數(shù)優(yōu)化和拓?fù)鋬?yōu)化算法研究奠定基礎(chǔ)。