王 童

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

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

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

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

本文針對傳統(tǒng)四節(jié)點膜單元無法進行彎曲的缺陷，基于單元內(nèi)力和彎矩的平衡原理，提出一種可進行彎曲修正的新型四節(jié)點膜單元模型。將其應(yīng)用于客車側(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］：

式中：ε為單元應(yīng)變；σ為單元Cauchy應(yīng)力；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)：

應(yīng)用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é)點法線等距的原理，計算彎曲后單元上下表面應(yīng)力。圖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，同時分別計算其兩個表面的應(yīng)力。規(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下對應(yīng)節(jié)點的坐標。

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

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

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

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

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

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

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

需要注意的是式（15）中，σx，σy，τxy為由板彎曲變形引起的三個應(yīng)力分量，式（13）給出的是由面內(nèi)變形和彎曲變形共同引起的應(yīng)力。故需進一步將應(yīng)力分量分解為面內(nèi)變形應(yīng)力σp和彎曲變形應(yīng)力σb兩部分。根據(jù)薄板彎曲理論，彎曲應(yīng)力應(yīng)該滿足σ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é)點法向力。應(yīng)用最終的單元節(jié)點內(nèi)力Fint對客車側(cè)翻一步碰撞算法進行改進，可有效修正彎曲變形的影響，同時提高算法的計算效率。

3 應(yīng)用實例

本文以某12 m 公路客車車身結(jié)構(gòu)的典型車身段作為分析對象，進行側(cè)翻碰撞試驗；將經(jīng)過彎曲修正的四節(jié)點膜單元應(yīng)用于側(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)抬升，車身段試驗?zāi)Ｐ碗S著試驗臺的升高旋轉(zhuǎn)，當(dāng)重心位置越過垂直中心線達到側(cè)翻臨界狀態(tài)后開始自由側(cè)翻下落，車身段骨架撞向地面，試驗結(jié)束。

圖4 所示為圖3 所示待試驗典型車身段有限元模型，共離散四邊形單元259 976 個，節(jié)點258 368個。材料性能參數(shù)采用Holloman冪次硬化法則強化系數(shù)k=710 MPa，應(yīng)變硬化指數(shù)n為0.2，初始應(yīng)變ε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è)翻試驗制備及測量過程的偶然性影響導(dǎo)致，但誤差均在工程允許范圍內(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%，兩者計算精度基本相當(dāng)；改進算法與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

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

4 結(jié)論

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

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

本文改進的客車側(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)碰撞模擬中的應(yīng)用進行探索，為更復(fù)雜的客車及乘用車正碰撞快速算法研究奠定堅實的前期基礎(chǔ)。