一種彎曲修正膜單元及對側(cè)翻一步碰撞算法的改進

王 童

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

近年來公共交通的發(fā)展呈現(xiàn)迅猛發(fā)展態(tài)勢，客車作為一種公共出行工具，逐步占據(jù)社會主流地位。同時，與客車相關(guān)的重大交通事故數(shù)量也逐步呈現(xiàn)上升趨勢，給人類的生命財產(chǎn)帶來巨大損失。側(cè)翻是客車最嚴重的事故類型之一，往往造成群死群傷的惡劣后果［1］。目前，國內(nèi)外許多學者已經(jīng)投入到客車側(cè)翻安全性的探索與改進研究中。比如，Keith等［2］發(fā)現(xiàn)高強鋼頂梁或很便宜的纖維增強復合塑料(FRP）玻璃可以阻止盒狀現(xiàn)象出現(xiàn)，并預防客車側(cè)翻的發(fā)生.

客車側(cè)翻一步碰撞算法，是筆者參考板料沖壓一步成型算法思想［3-5］，提出的一種用于客車側(cè)翻碰撞模擬的新算法。對于板料沖壓一步成型算法，王鵬等［6］已經(jīng)做過了一些研究。在他們的研究中，大部分金屬成形過程都通過應用基于塑性變形理論的一步有限元方法進行模擬，它可以快速預測坯料最終構(gòu)型及應力應變分布。借鑒該算法的核心思想，基于非線性全量理論，客車側(cè)翻一步碰撞算法與現(xiàn)有的LS-DYNA 等增量法軟件相比，在略微犧牲一點計算精度的條件下，計算效率可以得到大幅度提升，在車身設(shè)計初期，對客車的側(cè)翻安全性能進行快速評價。

算法中單元模型的選擇需兼顧精度和計算效率，需在兩者間尋找一個平衡點［7］。側(cè)翻過程結(jié)構(gòu)的變形特性主要為彎曲變形，所有膜單元及較簡單的三節(jié)點三角形殼單元均無法滿足要求［8-10］。由于具有低階、簡單、位移連續(xù)等優(yōu)點，初始側(cè)翻一步碰撞算法選擇了基于Mindlin 四邊形板單元和四節(jié)點膜單元組合而成的空間殼單元對結(jié)構(gòu)進行離散。而應用四節(jié)點膜單元進行模擬，結(jié)構(gòu)的總自由度會比殼單元減少一半，其計算時間約為原時間的60%左右，對算法模擬更有優(yōu)勢。同時，膜單元彎曲效應的模擬問題也必須解決。

本文針對傳統(tǒng)四節(jié)點膜單元無法進行彎曲的缺陷，基于單元內(nèi)力和彎矩的平衡原理，提出一種可進行彎曲修正的新型四節(jié)點膜單元模型。將其應用于客車側(cè)翻一步碰撞算法，簡便有效地考慮了客車側(cè)翻過程中的單元節(jié)點彎曲內(nèi)力。因其計算量遠小于傳統(tǒng)殼單元，故在有效保證計算精度的同時，計算效率可以得到顯著提升。

1 客車側(cè)翻一步碰撞算法基本理論[11]

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

將碰撞開始狀態(tài)的結(jié)構(gòu)作為原始構(gòu)形X0。此時車體未發(fā)生變形，結(jié)構(gòu)動能Ed由式（1）計算：

式中：M為車體質(zhì)量；Δh為車體重心下降高度；J為車體繞假定轉(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）計算：

式中：x0為結(jié)構(gòu)變形最大構(gòu)型。

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

式中：ε為單元應變；σ為單元Cauchy應力；Ve為單元體積；N為單元數(shù)。

判斷此刻的結(jié)構(gòu)形變能W與車體動能Ed是否滿足式（5）的能量關(guān)系假定：

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

由于車體結(jié)構(gòu)在空間內(nèi)變形過程無外力作用，滿足能量轉(zhuǎn)換關(guān)系的初始解U，節(jié)點失衡力R(U)已處于不平衡狀態(tài)：

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

對于式（7）中切線剛度矩陣的求解，由于需要迭代求解大型線性方程組，算法計算效率大大降低。故在初始算法中，課題組采用改進型擬共軛梯度法代替該切線剛度矩陣的計算，同時引入攝動原理，在無需計算切線剛度矩陣的情況下，對廣義失衡力采用一維帶寬存儲，通過梯度因子快速尋找廣義失衡力下降方向，使廣義失衡力“2”范數(shù)達到極小值，結(jié)構(gòu)碰撞能量達到平衡，獲得結(jié)構(gòu)的最終變形。

2 彎曲修正的四節(jié)點膜單元模型

本文提出的可進行彎曲修正的新型四節(jié)點膜單元模型如圖 1 所示。圖中i，j，k，l為單元節(jié)點在整體坐標系下的坐標；iu，ju，ku，lu為單元上表面 4 個節(jié)點在整體坐標系下的坐標；id，jd，kd，ld為單元下表面 4 個節(jié)點在整體坐標系下的坐標；ni，nj，nk，nl為單元各節(jié)點法線。與傳統(tǒng)殼單元模型不同，這種彎曲修正的四節(jié)點膜單元節(jié)點內(nèi)力由兩部分構(gòu)成：一部分是由面內(nèi)變形引起的節(jié)點切向力，另一部分是由彎曲變形引起的節(jié)點彎矩和法向力，如式（8）所示：

式中：Fint為節(jié)點的總內(nèi)力；Fm為因面內(nèi)變形而引起的節(jié)點內(nèi)力；Fb為因彎曲變形而引起的節(jié)點彎曲內(nèi)力。Fm和Fb的向量形式由式（9）表示：

面內(nèi)變形引起的節(jié)點內(nèi)力Fm的計算與普通膜單元相同。為計算單元節(jié)點彎曲內(nèi)力Fb，首先利用四節(jié)點膜單元沿節(jié)點法線等距的原理，計算彎曲后單元上下表面應力。圖1所示為最大變形構(gòu)形x0下的某一單元A(i，j，k，l)，單元厚度為t。該單元節(jié)點在 整 體 坐 標 系 下 的 坐 標 分 別 為i(xi，yi，zi)，j(xj，yj，zj)，k(xk，yk，zk)，l(xl，yl，zl)，各節(jié)點法線分別為ni(li，m，iqi)，nj(lj，m，jqj)，nk(lk，m k，qk)，nl(ll，m，lql)，節(jié)點法線由節(jié)點相鄰單元加權(quán)平均求得。沿著單元A的4個節(jié)點法線方向及負方向平行等距移動半個板厚可得其上、下兩個表面Au和Ad，同時分別計算其兩個表面的應力。規(guī)定上表面Au四 個 節(jié) 點 的 坐 標 分 別 為i(xui，yui，zui)，j(xuj，yuj，zuj)，k(xuk，yuk，zuk)，l(xul，yul，zul)，下 表 面Ad四個節(jié)點的坐標分別為i(xdi，ydi，zdi)，j(xdj，ydj，zdj)，k(xdk，ydk，zdk)，l(xdl，ydl，zdl)。單元上、下表面節(jié)點處的坐標可由式（10）計算：

圖1 可進行彎曲修正的四節(jié)點膜單元Fig.1 Four-node membrane element considering bending modification

這樣，可得到單元上、下表面節(jié)點處的位移。以上表面單元Au為例，整體坐標系下各節(jié)點的位移可由式（11）計算：

式中：uur表示節(jié)點ru(r=i，j，k，l)在整體坐標系中的位移向量；xur表示節(jié)點ru(r=i，j，k，l)在最大變形構(gòu)形x0下的整體坐標；Xur表示節(jié)點ru(r=i，j，k，l)在原始構(gòu)形X0下對應節(jié)點的坐標。

根據(jù)大變形幾何關(guān)系可以求得Green 變形張量下的單元主伸長λur(r=1，2，3)，在主應變空間的對數(shù)應變εur(r=1，2，3)可由式（12）計算：

考慮材料各向異性的塑性本構(gòu)關(guān)系，在比例加載條件下，上表面單元Au任一點的應力值σu可由式（13）計算［12］：

這樣，即可獲得上表面單元Au任一點的應力值σu，同理，可獲得下表面單元Ad任一點的應力值σd。

根據(jù)薄板彎曲問題的有限元理論，薄板的廣義內(nèi)力可由式（14）計算：

式中：Mx為垂直于x軸截面單位長度彎矩；My為垂直于y軸截面單位長度彎矩；Mxy為垂直于x（y）軸截面單位長度扭矩，Mxy=Myx。

由于沿z軸方向的應力成線性分布，薄板內(nèi)任意一點的應力可由式（15）計算：

對于圖1 中的模型，z=t/2。則虛擬板單元某一截面單位長度的彎矩和扭矩可由式（16）計算：

需要注意的是式（15）中，σx，σy，τxy為由板彎曲變形引起的三個應力分量，式（13）給出的是由面內(nèi)變形和彎曲變形共同引起的應力。故需進一步將應力分量分解為面內(nèi)變形應力σp和彎曲變形應力σb兩部分。根據(jù)薄板彎曲理論，彎曲應力應該滿足σub=-σdb。令Δ。

因單元內(nèi)部廣義內(nèi)力相互平衡，故需在單元邊界對式（16）進行積分。由彎曲變形引起的合彎矩可由式（17）計算：

式中：s為單元的邊界。

考慮到碰撞結(jié)束后單元處于平衡狀態(tài)，如圖2所示，圖中Fri，F(xiàn)rj，F(xiàn)rk，F(xiàn)rl，r=x，y，z為單元4個節(jié)點所受到的法向力，其作用在單元上的合外力為零。由單元平衡條件可得式（18）：

簡化式（18）可得到式（19）：

由于式（19）中包含F(xiàn)zi，F(xiàn)zj，F(xiàn)zk，F(xiàn)zl四個未知參數(shù)，但只有三個方程，故求解該方程組仍需要一個補充條件。本文提出一種基于單元節(jié)點法向力方差的最小極值法，通過對單元4個節(jié)點法向力求方差，對所構(gòu)成的函數(shù)求極小值。結(jié)合式（19）求解單元4個節(jié)點的法向力，所得各節(jié)點法向力的離散程度最小，計算結(jié)果最趨于穩(wěn)定。補充條件可由式（20）表示：

將式（19）代入式（20），可得到一個以Fzk為變量的一元二次方程f(Fzk)。對該方程求極小值，得到滿足各節(jié)點法向力方差最小的Fzk值，從而得到各節(jié)點法向力。應用最終的單元節(jié)點內(nèi)力Fint對客車側(cè)翻一步碰撞算法進行改進，可有效修正彎曲變形的影響，同時提高算法的計算效率。

3 應用實例

本文以某12 m 公路客車車身結(jié)構(gòu)的典型車身段作為分析對象，進行側(cè)翻碰撞試驗；將經(jīng)過彎曲修正的四節(jié)點膜單元應用于側(cè)翻一步碰撞算法進行模擬，并與初始算法及LS-DYNA 結(jié)果和側(cè)翻試驗結(jié)果對比，檢驗所提單元模型的實際工程效果。

圖2 單元局部坐標系及節(jié)點內(nèi)力示意圖Fig.2 Schematic diagram of element local coordinate system and nodal internal forces

圖3 為所選待側(cè)翻試驗典型車身段模型。根據(jù)標準GB 17578—2013［13］及ECE R66法規(guī)［14］，按照車身段實際重量配重，保證重心位置與實車相同，車身段由一個簡單的翻滾支架固定，最小離地間隙與實車保持一致。

圖3 典型車身段模型Fig.3 Typical bus body section model

將該車身段模型置于圖3所示側(cè)翻試驗臺固定軸邊緣位置，車身段處于自由靜止狀態(tài)；側(cè)翻試驗臺繞著固定軸以0.1(。)·S-1的準靜態(tài)速度緩慢旋轉(zhuǎn)抬升，車身段試驗模型隨著試驗臺的升高旋轉(zhuǎn)，當重心位置越過垂直中心線達到側(cè)翻臨界狀態(tài)后開始自由側(cè)翻下落，車身段骨架撞向地面，試驗結(jié)束。

圖4 所示為圖3 所示待試驗典型車身段有限元模型，共離散四邊形單元259 976 個，節(jié)點258 368個。材料性能參數(shù)采用Holloman冪次硬化法則強化系數(shù)k=710 MPa，應變硬化指數(shù)n為0.2，初始應變ε0為0，厚向異性系數(shù)r為1.3。

圖5 所示為經(jīng)過改進的客車側(cè)翻一步碰撞算法（簡稱“改進算法”）模擬的結(jié)構(gòu)最終變形云圖，圖6為初始算法結(jié)果，圖7 所示為LS-DYNA 仿真結(jié)果，圖8為側(cè)翻試驗結(jié)果。

據(jù)圖5～圖8 對比發(fā)現(xiàn)，4 種方式結(jié)構(gòu)最終變形趨勢吻合，改進算法計算結(jié)果與初始算法、LSDYNA及側(cè)翻碰撞結(jié)果趨勢基本一致。

為進一步對結(jié)構(gòu)變形進行定量分析，將側(cè)翻試驗提供的若干測點作為樣本點，在車身段封閉環(huán)①和②兩側(cè)立柱上各選取11個測點進行數(shù)據(jù)采集，得到兩側(cè)立柱的變形量情況，并作為與改進算法、初始算法和LS-DYNA 仿真分析結(jié)構(gòu)變形的對比數(shù)據(jù)，測點位置如圖9所示。參考側(cè)翻試驗中測點選取方式，有限元模型測點的選取位置如圖10 所示，兩點之間的距離為100 mm。

圖4 典型車身段有限元模型Fig.4 Finite element analysis model of typical bus body section

圖5 改進算法模擬的結(jié)構(gòu)最終變形云圖（單位：m）Fig.5 Final deformation contour plot of improved algorithm(unit:m)

圖6 初始算法模擬的結(jié)構(gòu)最終變形云圖（單位：m）Fig.6 Final deformation contour plot of original algorithm(unit:m)

圖7 LS-DYNA仿真的結(jié)構(gòu)最終變形云圖Fig.7 Final deformation contour plot of LS-DYNA

圖8 側(cè)翻試驗結(jié)構(gòu)最終變形Fig.8 Final deformation of structure of rollover test

圖9 側(cè)翻試驗測點選取位置Fig.9 Gauging point location for selection of roll over test

改進算法、初始算法、LS-DYNA 仿真及側(cè)翻試驗的封閉環(huán)①和②兩側(cè)立柱的變形量數(shù)據(jù)如表1所示。由表1可知，部分數(shù)據(jù)存在微小偏差，是由于受實際側(cè)翻試驗制備及測量過程的偶然性影響導致，但誤差均在工程允許范圍內(nèi)。

圖10 有限元模型測點選取位置Fig.10 Gauging point location for selection of finite element model

表1 各封閉環(huán)兩側(cè)立柱變形量統(tǒng)計Tab.1 Data of deformations of both sides of each pillar

圖11 和圖12 所示為典型車身段的封閉環(huán)①和②兩側(cè)立柱變形量對比柱狀圖。

圖11 封閉環(huán)①兩側(cè)立柱變形量對比Fig.11 Comparison of deformations of closed-loop ①

圖12 封閉環(huán)②兩側(cè)立柱變形量對比Fig.12 Comparison of deformations of closed-loop ②

從圖12可看出，由于受到試驗制備及實際測量過程的誤差影響，側(cè)翻試驗數(shù)據(jù)在穩(wěn)定性方面波動相對比較明顯，但對改進算法的有效性檢驗，具有一定工程參考價值。通過對比分析圖11 和圖12，4 種結(jié)果獲得的各測點數(shù)據(jù)走勢基本一致，驗證了結(jié)構(gòu)最終變形趨勢吻合的結(jié)論；改進算法與初始算法間的誤差平均為2.2%，兩者計算精度基本相當；改進算法與LS-DYNA 仿真分析之間誤差平均約為4.1%，計算精度基本滿足實際需要；改進算法與側(cè)翻試驗結(jié)果的誤差平均約為11.3%，小于有限元工程計算誤差的經(jīng)驗值范圍15%，計算精度在實際可接受范圍內(nèi)。

接著對比計算效率。算法模擬采用的硬件平臺為 DELL 工作站，配置 8 核 Intel i7 處理器、32G 內(nèi)存及8T固態(tài)硬盤。表2為改進算法與初始算法和LSDYNA仿真的模擬效率對比結(jié)果，由表2可知，改進算法模擬時間約為初始算法的3/5，改進算法的模擬時間約為LS-DYNA仿真的1/15。

表2 模擬效率對比Tab.2 Comparison of simulation efficiency

綜合對比計算精度和計算效率，改進算法雖精度與初始算法基本相當，但計算效率得到了大幅提升，檢驗了本文所提新型單元在實際工程應用中的有效性，可在客車結(jié)構(gòu)設(shè)計初期，在基本保證計算精度的前提下，更快速地獲得側(cè)翻碰撞結(jié)構(gòu)的最終變形，對結(jié)構(gòu)側(cè)翻碰撞安全性進行快速評價。

4 結(jié)論

本文針對傳統(tǒng)四節(jié)點膜單元無法進行彎曲的缺陷，基于單元內(nèi)力和彎矩的平衡原理，提出一種可進行彎曲修正的新型四節(jié)點膜單元模型。將其應用于客車側(cè)翻一步碰撞快速模擬算法，簡便有效地考慮了客車側(cè)翻過程中由彎曲變形引起的節(jié)點內(nèi)力。通過對典型車身段進行側(cè)翻碰撞模擬，并與初始算法、LS-DYNA及側(cè)翻試驗結(jié)果對比，4種方法結(jié)構(gòu)最終變形趨勢比較吻合，改進算法在略微犧牲一點計算精度的情況下，計算效率得到了顯著提高。

本文所提出可進行彎曲修正的四節(jié)點膜單元模型，對一般膜單元無法計算單元法向力的缺點進行了修正。整個改進過程依舊采用膜單元理論，節(jié)點的整體自由度并未增加，整體剛度矩陣帶寬也未改變。既保持了原始膜單元計算效率高的優(yōu)點，又獲得了可觀的計算精度，具有一定的實際工程應用價值。

本文改進的客車側(cè)翻一步碰撞算法，在基本保證計算精度的同時，計算效率進一步得到了大幅提高，可在結(jié)構(gòu)設(shè)計初期對客車側(cè)翻安全性能進行快速評價，縮短產(chǎn)品開發(fā)周期，降低研發(fā)成本；同時，為后續(xù)開展基于客車側(cè)翻碰撞安全性能的結(jié)構(gòu)尺寸優(yōu)化、參數(shù)優(yōu)化、拓撲優(yōu)化及靈敏度分析等算法提供必要支撐條件；并可對非線性隱式理論在結(jié)構(gòu)碰撞模擬中的應用進行探索，為更復雜的客車及乘用車正碰撞快速算法研究奠定堅實的前期基礎(chǔ)。