積分迭代模式優化對側翻一步碰撞算法的改進*

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

（長安大學汽車學院，西安710064）

前言

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

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

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

1 側翻一步碰撞算法基本理論

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

將碰撞開始狀態結構作為初始構形{X0}。此時車體未發生變形，結構動能Ed為

式中：M為車體質量；Δh為車體質心下降高度；J為車體繞側翻假定轉軸的轉動慣量；ω為車體角速度。

碰撞開始狀態結構各節點速度{v0}由式（2）計算：

式中：ri為各節點到側翻假定轉軸距離；n為節點數。

在結構最大變形狀態，車體結構產生明顯變形。側翻一步碰撞算法中，結構最大變形不確定須假定一個最大變形構形{x0}。各節點位移{U0}由式（3）計算：

車體動能Ed在碰撞中主要轉換為結構形變能，結構形變能W由式（4）計算［13］：

式中：{ε}為單元塑性應變；{σ}為單元塑性應力；Ve為單元體積；N為單元數。

結構形變能W與車體動能Ed理應相等，即

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

對滿足能量轉換關系初始解{U}，節點失衡力{R(U)}已處于不平衡狀態：

應用Newton?Raphson法，解決節點失衡力不平衡問題。對初始解{U}按式（7）迭代求解，使式（6）達到平衡，得到結構最終變形。

2 單元積分點模式優化

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

2.1 單元面內積分方案

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

式中：r、s為積分點坐標；h為積分系數；i=1，j=1。

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

為實現局部減縮積分，將單元［B］矩陣和應力向量{σ}進行分解，令：

其中

單元內力為

令：

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

2.2 單元厚向積分方案

2.2.1 厚向五點高斯積分

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

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

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

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

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

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

接著，對應用厚向五點3個積分方案分析計算時間與LS?DYNA時間進行對比，結果如表2所示。

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

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

2.2.2 厚向七點高斯積分

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

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

厚向七點積分與圖1中所示的單元厚向五點積分示意圖類似，只是積分點數量和位置有所不同。厚向七點積分中積分點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，對應的積分系數分別為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。

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

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

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

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

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

2.2.3 厚向九點高斯積分

當單元厚向采用九點積分時，隨著單元厚向積分點數進一步增加，其積分方案匹配方式繼續增加，共設計5種單元積分方案。單元厚向九點積分各積分層面內減縮積分與局部減縮積分匹配方案設計如表5所示。

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

厚向九點積分與圖1中所示的單元厚向五點積分示意圖類似，只是積分點數量達到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，對應的積分系數分別 為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。

同樣，將厚向九點各積分方案應用于客車側翻一步碰撞快速算法進行改進，獲得封閉環兩側立柱各測點變形量和計算時間。統計各方案中11個測點變形量數據，將厚向九點各積分方案計算結果與LS?DYNA結果進行對比，封閉環兩側立柱各測點變形量對比柱狀圖如圖4所示。

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

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

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

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

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

3 廣義失衡力迭代模式優化

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

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

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

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

節點廣義失衡力計算完成后，將每個節點中的失衡力分離出來單獨計算失衡力“2”范數：

式中n為車體結構中單元節點總數。

同樣，將每個節點中失衡力矩分離出來單獨計算力矩“2”范數：

取一組權系數將上述力和力矩“2”范數進行整合，得到加權調整后的廣義失衡力“2”范數，即

式中α、β分別為力和力矩的權系數，其和為2。

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

圖5 計算時間-權重比曲線

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

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

此時，車體結構處于穩定狀態，該狀態下車體結構變形即為客車側翻碰撞后的結構最終變形。

4 應用實例

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

4.1 模型建立

車身段側翻試驗參考歐洲ECE R66法規設計，試驗臺高800 mm，試驗模型配重后質量2 409 kg，質心高度1.5 m，典型車身段試驗模型如圖6所示。

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

對該車身段試驗模型進行簡化并假設，建立車身段有限元模型如圖7所示。模型共離散270 171個單元，269 894個節點，骨架鋼材選用Q345鋼，材料彈性模量E=2.06×1011Pa，泊松比μ=0.3，密度ρ=7800 kg/m3，屈服強度σs=345 MPa。

圖7 車身段有限元模型

車身段側翻試驗模型以及相應的有限元模型中所包含的A、B兩個封閉環，參考企業側翻試驗過程中測點選取方式，在封閉環兩側立柱共選取11個測點，左側立柱內表面選取6個測點，間隔100 mm，右側立柱內表面選取5個測點，間隔100 mm，各測點具體位置見圖8。

圖8 封閉環模型測點選取方式示意圖

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

圖9 封閉環模型LS?DYNA分析測點變形量柱狀圖

4.2 結果對比

初始算法結構變形情況如圖10所示，改進算法結構變形情況如圖11所示，LS?DYNA分析結構變形情況如圖12所示。

圖10 初始算法結構變形情況

圖12 LS?DYNA結構變形情況

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

圖11 改進算法結構變形情況

側翻試驗測點變形量數據來自企業實車開發過程中車身段側翻試驗，數據由企業技術人員協助測得。側翻試驗、LS?DYNA分析、初始算法和改進算法所得封閉環A、B兩側立柱變形量數據統計如表6所示。

表6 典型車身段測點變形量統計 mm

封閉環A兩側立柱各測點變形量對比見圖13，封閉環B兩側立柱各測點變形量對比如圖14所示。

圖13 封閉環A各測點變形量對比

圖14 封閉環B各測點變形量對比

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

接著，對計算效率進行對比，結果如表7所示。

表7 計算效率對比

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

5 結論

為進一步提高側翻一步碰撞算法計算效率且保證計算精度，為后續對客車車身結構進行基于側翻安全的拓撲優化提供支撐，在客車側翻一步碰撞初始算法基礎上，分析算法運行過程中的關鍵環節，從計算量較大的積分及迭代部分入手，對初始算法積分及迭代方法進行探索和改進。經過計算分析，改進算法在基本保證計算精度的同時，計算效率得到顯著提升。本文中所提單元積分點模式及廣義失衡力迭代模式改進方案，具有良好應用效果。

但是，初始算法本身單元模型基于Mindlin殼單元，受單元模型理論限制，算法計算精度不能大幅提升，后續仍需對算法單元模型進行改進研究，進一步提升算法計算精度和效率，以便在車身設計初期，對結構側翻碰撞安全性進行更加快速準確評價，進一步縮短產品開發周期，為后續針對客車側翻安全性的靈敏度分析、參數優化和拓撲優化算法研究奠定基礎。