新型DKT組合薄殼單元及對側(cè)翻一步碰撞算法的改進(jìn)*

王 童，陳軼嵩，劉永濤

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

前言

汽車保有量急劇上漲，人們出行生活越來越便利［1］。同時(shí)，惡劣交通事故數(shù)量呈現(xiàn)逐步上升趨勢，給乘客生命財(cái)產(chǎn)帶來重大損失［2］。其中，客車側(cè)翻常造成群死群傷的嚴(yán)重后果，是最嚴(yán)重的交通事故類型之一［3］。如何最大限度預(yù)防側(cè)翻事故的發(fā)生，加強(qiáng)客車結(jié)構(gòu)側(cè)翻碰撞安全性研究，不斷提高客車側(cè)翻碰撞安全性能，已成為全社會共同關(guān)注并努力解決的重要問題之一［4］。

客車側(cè)翻一步碰撞算法，是參考板料沖壓一步成型算法思想，提出的一種用于客車側(cè)翻碰撞模擬分析的新算法［5］。借鑒該算法核心思想，側(cè)翻一步碰撞算法同樣基于塑性全量理論，利用側(cè)翻碰撞過程中能量的轉(zhuǎn)換關(guān)系。與現(xiàn)有LS?DYNA等增量法軟件相比，在略微犧牲一點(diǎn)計(jì)算精度情況下，大幅提升了計(jì)算效率。算法可在車身設(shè)計(jì)初期，對結(jié)構(gòu)側(cè)翻碰撞安全性進(jìn)行快速評價(jià)，縮短產(chǎn)品開發(fā)周期，同時(shí)，可為后續(xù)針對客車側(cè)翻安全性的靈敏度分析、參數(shù)優(yōu)化和拓?fù)鋬?yōu)化算法研究提供必要支撐條件［6］。

在有限元數(shù)值模擬過程中，不同單元模型選擇會對算法計(jì)算結(jié)果和效率產(chǎn)生較大影響［7］。因側(cè)翻一步碰撞算法基于塑性全量理論，車體骨架由若干大幾何尺寸薄壁桿件構(gòu)成，側(cè)翻碰撞過程主要受力特性為彎曲，故所有膜單元及比較簡單的三角形薄殼單元均無法滿足計(jì)算要求［8］。作者在客車側(cè)翻一步碰撞初始算法中采用了四節(jié)點(diǎn)Mindlin組合殼單元作為分析單元類型。為進(jìn)一步提高算法計(jì)算精度和效率，對四節(jié)點(diǎn)薄殼單元進(jìn)行翹曲修正改進(jìn)，基于改進(jìn)后的四節(jié)點(diǎn)薄殼單元和塑性全量理論，將DKT（discreted kirchhoff theory）四節(jié)點(diǎn)薄板單元和傳統(tǒng)四節(jié)點(diǎn)等參膜單元組合，構(gòu)造一種新型四節(jié)點(diǎn)DKT組合薄殼單元，替代初始算法四節(jié)點(diǎn)Mindlin組合殼單元，對側(cè)翻一步碰撞算法改進(jìn)。預(yù)計(jì)改進(jìn)后算法在模擬效率基本不變的同時(shí)，計(jì)算精度得到一定程度的提高。

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

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

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

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

碰撞開始狀態(tài)結(jié)構(gòu)各節(jié)點(diǎn)速度{v0}為

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

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

車體結(jié)構(gòu)動能Ed在碰撞中主要轉(zhuǎn)換為結(jié)構(gòu)形變能，結(jié)構(gòu)形變能W為［14］

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

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

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

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

式中：{Fext(Ui)}為廣義失衡外力；{Fint(Ui)}為廣義失衡內(nèi)力。

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

2 四節(jié)點(diǎn)薄殼單元翹曲修正改進(jìn)

客車側(cè)翻碰撞結(jié)構(gòu)變形過程中，部分結(jié)構(gòu)可能存在變形過大問題，導(dǎo)致較明顯單元翹曲［15］，影響算法計(jì)算精度。對所有四節(jié)點(diǎn)薄殼單元，首先考慮單元翹曲修正問題，對其進(jìn)行改進(jìn)，接著建立本文所需四節(jié)點(diǎn)DKT組合薄殼單元模型，對側(cè)翻一步碰撞初始算法進(jìn)行改進(jìn)。

2.1 考慮翹曲修正的四節(jié)點(diǎn)薄殼單元坐標(biāo)轉(zhuǎn)換關(guān)系

提取四節(jié)點(diǎn)薄殼單元中性層。在整體坐標(biāo)系下，產(chǎn)生翹曲的四節(jié)點(diǎn)單元中性層為i'j'k'l'，各節(jié)點(diǎn)坐 標(biāo) 為i'(Xi'，Yi'，Zi')、j'(Xj'，Yj'，Zj')、k'(Xk'，Yk'，Zk')、l'(Xl'，Yl'，Zl')。將翹曲后單元中性層，以對角線i'k'作為軸線，分割為兩個三角形Δi'l'k'和Δi'j'k'，展開到過k'點(diǎn)、以n→為法向量的平面內(nèi)，可以得到不產(chǎn)生翹曲的四節(jié)點(diǎn)單元i'j″k'l″，如圖1所示。

圖1 四節(jié)點(diǎn)薄殼單元中性層翹曲展開前后示意圖

整體坐標(biāo)系內(nèi)，翹曲單元4個節(jié)點(diǎn)坐標(biāo)均已知。對于三角形 Δi'j″k'，邊長Li'j'=Li'j″、Lk'j'=Lk'j″、Lk'l'=Lk'l″、Li'l'=Li'l″及軸線Li'k'長度均可計(jì)算。

2.2 基于形心位置不變假設(shè)的最終展開構(gòu)形

翹曲四節(jié)點(diǎn)薄殼單元在展開前后形心位置保持不變。在整體坐標(biāo)系內(nèi)，將展開后的四節(jié)點(diǎn)單元形心平移至展開前初始位置。由于單元翹曲時(shí)形心位置，在單元展開后已發(fā)生變化，需將其平移至展開前位置。以展開后四節(jié)點(diǎn)單元i'j″k'l″形心為基點(diǎn)，因展開后四節(jié)點(diǎn)單元形心需積分，計(jì)算不便，故將單元各節(jié)點(diǎn)坐標(biāo)平均值近似作為形心O″，如圖2所示。

圖2 整體坐標(biāo)系下展開后單元形心平移示意圖

因翹曲單元i'j'k'l'形心無法準(zhǔn)確確定，故本文以翹曲單元各節(jié)點(diǎn)坐標(biāo)平均值為近似形心O'，將展開后四節(jié)點(diǎn)單元形心O″平移到翹曲展開前點(diǎn)O'，得到最終翹曲修正單元構(gòu)形。翹曲展開后單元形心O″坐標(biāo)為

翹曲單元展開前形心O'坐標(biāo)為

可得平移形心后四節(jié)點(diǎn)單元ijkl各節(jié)點(diǎn)坐標(biāo)為

式中ξ分別取i、j、k、l，?分別取i'、j″、k'、l″。

經(jīng)過上述計(jì)算過程，可對已發(fā)生翹曲四節(jié)點(diǎn)薄殼單元進(jìn)行修正，在整體坐標(biāo)系內(nèi)獲得翹曲修正改進(jìn)后單元（中性層）各節(jié)點(diǎn)坐標(biāo)。

3 改進(jìn)的四節(jié)點(diǎn)DKT組合薄殼單元構(gòu)造及對側(cè)翻一步碰撞初始算法的改進(jìn)

對經(jīng)過翹曲修正的四節(jié)點(diǎn)薄殼單元，為便于計(jì)算單元應(yīng)力應(yīng)變和節(jié)點(diǎn)內(nèi)力，須重新將整體坐標(biāo)系下翹曲修正改進(jìn)后的單元各節(jié)點(diǎn)坐標(biāo)轉(zhuǎn)換至局部坐標(biāo)下，進(jìn)行算法后續(xù)計(jì)算過程［16］。

3.1 四節(jié)點(diǎn)DKT薄板單元

離散Kirchhoff理論基于以下基本假定：（1）單元中面的垂直法線，變形前后仍與中面垂直，且長度不變；（2）單元厚向應(yīng)變εz=0，厚向應(yīng)力σz=0，切向應(yīng)力γxz=0、γyz=0。

在上述基本假定基礎(chǔ)上，建立基于離散Kirchhoff理論的四節(jié)點(diǎn)薄板單元。單元局部坐標(biāo)系內(nèi)，模型約定僅單元4個角節(jié)點(diǎn)和4條邊的4個邊中節(jié)點(diǎn)保持Kirchhoff直法線假設(shè)，且每個角節(jié)點(diǎn)有自由 度(wξ，θξx，θξy)(ξ=i，j，k，l)，邊 中 節(jié) 點(diǎn) 有 自 由 度(θ?x，θ?y)(?=i'，j'，k'，l')，如圖3所示。

圖3 離散Kirchhoff理論節(jié)點(diǎn)自由度假設(shè)

四節(jié)點(diǎn)DKT薄板單元中性層內(nèi)任意一點(diǎn)轉(zhuǎn)動自由度θx和θy可由單元4個角節(jié)點(diǎn)進(jìn)行插值，即

式中 ：{}為單 元各節(jié) 點(diǎn)自由 度{wi,θix,θiy,wj,θjx,θjy,wk,θkx,θky,wl,θlx,θly}T，且有[Hx]=[Hx1Hx2...Hx12]，[Hy]=[Hy1Hy2...Hy12]。

四節(jié)點(diǎn)DKT薄板單元Bb矩陣為

其中

3.2 四節(jié)點(diǎn)等參膜單元

對側(cè)翻碰撞過程中單元拉伸效應(yīng)，本文中選擇雙線性四節(jié)點(diǎn)等參膜單元模擬，各節(jié)點(diǎn)自由度向量如圖4所示。

圖4 四節(jié)點(diǎn)等參膜單元示意圖

經(jīng)等參變換后單元形函數(shù)為

式中ri和si(i=1，2，3，4)為經(jīng)過等參變換后各節(jié)點(diǎn)的自然坐標(biāo)，其值為±1。

四節(jié)點(diǎn)等參膜單元B m矩陣為

3.3 基于塑性全量理論的節(jié)點(diǎn)廣義內(nèi)力計(jì)算

將上述3.1節(jié)部分所獲四節(jié)點(diǎn)DKT薄板單元和3.2節(jié)部分所獲四節(jié)點(diǎn)等參膜單元B矩陣進(jìn)行組合，得到所需四節(jié)點(diǎn)DKT組合薄殼單元。將其應(yīng)用于基于塑性全量理論的側(cè)翻一步碰撞算法，計(jì)算客車側(cè)翻碰撞過程產(chǎn)生的單元變形及各節(jié)點(diǎn)廣義內(nèi)力。

3.3.1 塑性應(yīng)力計(jì)算

左Cauchy?Green變形矩陣逆陣[G]-1為

其中

對空間四節(jié)點(diǎn)DKT組合薄殼單元內(nèi)部任意一點(diǎn)i，由于受到彎曲變形影響，該點(diǎn)所在橫向上任意一點(diǎn)位移不同。令薄殼厚度為t，為計(jì)算單元應(yīng)變、應(yīng)力及節(jié)點(diǎn)廣義內(nèi)力，經(jīng)對多組方案對比分析，本文在單元厚向選擇5點(diǎn)高斯積分［17］，各積分層均采用全積分方式，對板彎曲和膜拉伸兩種變形同時(shí)進(jìn)行考慮，如圖5所示。

圖5 薄殼單元的向量關(guān)系示意圖

四節(jié)點(diǎn)等參膜單元內(nèi)部任意一點(diǎn)位移，可由位移插值函數(shù)計(jì)算：

四節(jié)點(diǎn)DKT薄板單元內(nèi)部任意點(diǎn)轉(zhuǎn)動自由度θx和θy可由單元4個角節(jié)點(diǎn)進(jìn)行插值，由式（12）計(jì)算。在圖5中，L1、L2、L3、L4、L5為5點(diǎn)高斯積分求積位置因 子 ，數(shù) 值 分 別 為0.906 179 85、0.538 469 31、0、-0.906 179 85、-0.538 469 31。

各積分層單元局部坐標(biāo)系下各個方向?qū)?shù)應(yīng)變?yōu)椋?8］

單元各積分層面內(nèi)或彎曲變形產(chǎn)生的塑性應(yīng)力σs1、σs2、σs3、σs4、σs5可由式（21）計(jì)算［19］，其中σs3為膜應(yīng)力，σs1、σs2、σs4、σs5為板應(yīng)力。

3.3.2 節(jié)點(diǎn)廣義內(nèi)力計(jì)算

由3.3.1節(jié)部分所獲各積分層塑性應(yīng)力，各積分層廣義內(nèi)力為

單元各節(jié)點(diǎn)廣義內(nèi)力可對各積分層廣義內(nèi)力組合得到，即

應(yīng)用坐標(biāo)轉(zhuǎn)換矩陣[T]-1可將局部坐標(biāo)系下單元各節(jié)點(diǎn)內(nèi)力及內(nèi)力矩轉(zhuǎn)換至空間坐標(biāo)系下，即

其中

車體結(jié)構(gòu)各節(jié)點(diǎn)廣義內(nèi)力向量{F}int可由空間內(nèi)各單元節(jié)點(diǎn)廣義內(nèi)力累加，即

接著可對初始解各節(jié)點(diǎn)廣義失衡力進(jìn)行Newton?Raphson迭代，得到算法所需結(jié)構(gòu)最終變形。

4 應(yīng)用實(shí)例

為檢驗(yàn)本文所提新型DKT組合薄殼單元對客車側(cè)翻一步碰撞初始算法改進(jìn)的有效性，以某企業(yè)實(shí)際開發(fā)的一款12 m公路客車典型車身段作為分析對象，通過對車身段模型適當(dāng)簡化，將側(cè)翻一步碰撞改進(jìn)算法與初始算法、LS?DYNA和實(shí)車試驗(yàn)結(jié)果進(jìn)行對比，檢驗(yàn)所提單元模型實(shí)際效果。

車身段試驗(yàn)?zāi)Ｐ腿鐖D6所示，試驗(yàn)?zāi)Ｐ突緟?shù)如下：總質(zhì)量2 409 kg，質(zhì)心離地高度1.5 m，X方向距離前端面為0.937 m，Y方向?yàn)?.012 m，如表1所示。

表1 典型車身段基本參數(shù)

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

車身段側(cè)翻臨界角為39o。參照典型車身段試驗(yàn)?zāi)Ｐ停喕蟮湫蛙嚿矶谓Y(jié)構(gòu)有限元模型如圖7所示。共離散四節(jié)點(diǎn)單元259 976個，節(jié)點(diǎn)258 368個。黃色部分結(jié)構(gòu)為車身骨架結(jié)構(gòu)，綠色部分結(jié)構(gòu)為側(cè)翻翻轉(zhuǎn)支架結(jié)構(gòu)，紅色部分為乘員生存空間。與側(cè)翻試驗(yàn)?zāi)Ｐ拖嗤?jì)算模型骨架及側(cè)翻翻轉(zhuǎn)支架結(jié)構(gòu)材料均選用Q345鋼，材料彈性模量E=2.06×1011Pa，泊松比μ=0.3，密度ρ=7800 kg/m3，屈服強(qiáng)度σs=345MPa；乘員生存空間為剛性材料，密度取值盡量小，ρ=10 kg/m3，以免影響模擬結(jié)果精度［20］。

圖7 典型車身段有限元模型

圖8所示為應(yīng)用所提新型DKT組合薄殼單元改進(jìn)的側(cè)翻一步碰撞算法（以下簡稱“改進(jìn)算法”）模擬的結(jié)構(gòu)最終變形云圖。圖9所示為初始算法模擬的結(jié)構(gòu)最終變形云圖。圖10所示為LS?DYNA仿真的結(jié)構(gòu)最終變形云圖。圖11所示為側(cè)翻試驗(yàn)結(jié)構(gòu)變形結(jié)果圖。

圖8 改進(jìn)算法模擬結(jié)構(gòu)最終變形云圖

圖9 初始算法模擬結(jié)構(gòu)最終變形云圖

圖10 LS?DYNA仿真結(jié)構(gòu)最終變形云圖

圖11 典型車身段側(cè)翻試驗(yàn)結(jié)構(gòu)最終變形

將圖8～圖11進(jìn)行對比可以看出，4種模擬方式所得結(jié)構(gòu)最終變形形態(tài)趨勢基本吻合，應(yīng)用所提新型DKT組合薄殼單元對側(cè)翻一步碰撞算法進(jìn)行改進(jìn)，在實(shí)際工程應(yīng)用中具有一定合理性。

為使上述對比結(jié)論更具說服力，需對4種方式結(jié)構(gòu)最終變形定量對比。車身段側(cè)翻碰撞過程中，底架結(jié)構(gòu)幾乎不產(chǎn)生變形，且與乘員生存空間的侵入量無直接關(guān)聯(lián)。在企業(yè)對該車身段模型側(cè)翻試驗(yàn)中，為獲得結(jié)構(gòu)變形數(shù)據(jù)，在典型車身段封閉環(huán)①和②兩側(cè)立柱內(nèi)表面各選取11個測點(diǎn)進(jìn)行數(shù)據(jù)采集。上述4種方式各立柱變形量如表2所示。

表2 各封閉環(huán)兩側(cè)立柱變形量統(tǒng)計(jì) mm

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

圖12 封閉環(huán)①兩側(cè)立柱變形量對比

圖13 封閉環(huán)②兩側(cè)立柱變形量對比

從圖12和圖13可看出，受試驗(yàn)制備及實(shí)際測量過程誤差的影響，側(cè)翻試驗(yàn)數(shù)據(jù)在穩(wěn)定性方面波動相對較明顯，但對改進(jìn)算法實(shí)際應(yīng)用有效性檢驗(yàn)，具有一定工程參考價(jià)值。通過對比分析發(fā)現(xiàn)，4種方式所獲結(jié)構(gòu)各測點(diǎn)數(shù)據(jù)走勢基本一致，驗(yàn)證了結(jié)構(gòu)最終變形趨勢吻合的結(jié)論。改進(jìn)算法與LS?DYNA平均誤差約為1.9%，精度略有犧牲，但比初始算法提高約1.5%；與側(cè)翻試驗(yàn)平均誤差約為7.1%，小于有限元工程計(jì)算誤差經(jīng)驗(yàn)值范圍（即15%），計(jì)算精度在實(shí)際工程接受范圍內(nèi)。

接著對算法模擬效率進(jìn)行對比。將側(cè)翻一步碰撞改進(jìn)算法模擬效率與初始算法及LS?DYNA進(jìn)行對比，對比結(jié)果如表3所示。

表3 模擬效率對比

從表3可看出，改進(jìn)算法模擬時(shí)間比初始算法稍快4 min，改進(jìn)算法模擬時(shí)間約為LS?DYNA的1/10，初始算法約為LS?DYNA的1/9，改進(jìn)算法模擬效率有一定提升。結(jié)合立柱變形量對比結(jié)論，改進(jìn)算法在模擬效率些許提升的同時(shí)，計(jì)算精度也有一定程度的提高，可在客車結(jié)構(gòu)設(shè)計(jì)初期，更快更準(zhǔn)確獲得側(cè)翻碰撞結(jié)構(gòu)最終變形，對結(jié)構(gòu)側(cè)翻碰撞安全性進(jìn)行快速評價(jià)，縮短產(chǎn)品開發(fā)周期。同時(shí)，可為后續(xù)針對客車側(cè)翻安全性的靈敏度分析、參數(shù)優(yōu)化及拓?fù)鋬?yōu)化算法研究提供必要支撐條件，檢驗(yàn)了本文所提新型單元實(shí)際效果和工程應(yīng)用價(jià)值。

5 結(jié)論

針對客車側(cè)翻碰撞結(jié)構(gòu)變形模擬過程中，因部分結(jié)構(gòu)變形過大導(dǎo)致的單元翹曲問題，先對四節(jié)點(diǎn)薄殼單元進(jìn)行了翹曲修正。接著，在翹曲修正后薄殼單元中性層基礎(chǔ)上，基于塑性全量理論，將四節(jié)點(diǎn)DKT薄板單元和四節(jié)點(diǎn)等參膜單元進(jìn)行組合，獲得新型四節(jié)點(diǎn)DKT組合薄殼單元，替代側(cè)翻一步碰撞初始算法中四節(jié)點(diǎn)Mindlin組合殼單元，對側(cè)翻一步碰撞算法進(jìn)行改進(jìn)。以某12 m公路客車典型車身段作為研究對象，將改進(jìn)算法模擬結(jié)果與初始算法、LS?DYNA及側(cè)翻試驗(yàn)結(jié)果進(jìn)行對比，檢驗(yàn)了本文所提新型單元的實(shí)際中的有效性。

所提新型DKT組合薄殼單元，在厚向及各積分層積分方式和積分點(diǎn)數(shù)量對側(cè)翻一步碰撞算法模擬精度及效率的影響尚未考慮，可能會對模擬結(jié)果有較大影響，在后續(xù)對此做進(jìn)一步深入探討研究，使側(cè)翻一步碰撞算法性能得到進(jìn)一步提升。