復合材料薄膜充氣床墊的多物質ALE數值模擬

楊顏志，金先龍，2，張偉偉

(1.上海交通大學 機械與動力工程學院，上海 200240;2.上海交通大學 機械系統與振動國家重點實驗室，上海 200240)

氣囊結構由于其良好的減震作用，已廣泛應用于交通安全、產品防震、航天器登陸與回收等存在較大沖擊過載的領域［1-3］。傳統的氣囊設計多采用物理試驗對其進行性能分析，隨著計算機技術的發展，仿真技術在氣囊設計中的作用越來越大，已成為氣囊設計的重要手段［4-5］。

目前，國內外在氣囊仿真研究中常用的方法有兩種。一種是基于熱力學原理的控制體積方法(CV)，該方法比較簡便計算需求較少，但無法準確計算氣囊變形過程，計算精度只能滿足一般的工程需要［6-7］。另一種是基于有限元技術的ALE方法，它可以準確計算不同時刻氣囊的變形，同時精確描述充氣氣體、外部空氣與氣囊織物表面相互作用。該方法一般用于氣囊特性的詳細仿真，它具有計算準確的優勢，但比較復雜計算需求較大［8］。

本文所研究某型充氣床墊是典型的氣囊結構，該床墊由新型柔性薄膜材料制成，力學性能復雜。床墊充氣物理試驗過程中當充氣內壓達到25 kPa時發生局部爆破，本文仿真的主要目的包括兩部分:第一是模擬充氣內壓25 kPa時床墊爆破過程，在此基礎上提出結構改進方案，并作進一步仿真。第二是模擬床墊載人工作狀況，通過參數分析方法判斷床墊的最佳工作內壓。因爆破工況和載人工作工況的復雜性，為準確模擬床墊的變形與受力過程，本文采用多物質ALE方法對其進行詳細數值模擬。計算在上海超級計算中心“魔方”超級計算機上采用并行計算方法完成，計算平臺為Ls-dyna軟件和Optimus優化平臺。

1 ALE有限元方法

ALE有限元方法是針對傳統Lagrangian描述和Eulerian描述的不足，而提出的新的描述方法，該方法綜合了Lagrangian描述與Eulerian描述的優點，網格以給定的方式運動，而材料的運動并不與網格一致，因此即可追蹤自由表面，又可保證大變形問題的計算精度與計算穩定性［9］。ALE方法的計算網格獨立于變形體和空間運動，其根據需要自由選擇運動狀態，ALE方法具體描述如下。

如圖1所示，連續介質在t0時刻的初始構形記為ΩX，t時刻的現時構形記為ΩY，同時引入獨立于初始構形和現時構形的參照構形ΩZ。物質、空間和參考點分別由Lagrangian坐標系中X、Eulerian坐標系中Y和參考坐標系中Z描述。

在時間t時刻，空間點既可以通過材料運動形式Y=φ(X，t)來反映，又可以表現為網格運動形式Y=φ(Z，t)。而在現時構形中物質點的運動速度v和參考點的運動速度u可表示為:

圖1 ALE方法的運動描述Fig.1 Motion description in ALE

當u=v時，ALE描述轉化為Lagrangian描述;當u=0時，ALE描述轉化為Eulerian描述;只有當u≠v≠0時，對應一般的ALE描述。

ALE方法的控制方程(質量、動量和能量守恒方程)可表示為:

式中:ρ為密度，v為物質速度，w為對流速度，b為單位體積力，σ為柯西應力張量，E為能量。

ALE方法分為單物質ALE方法和多物質ALE方法。多物質ALE方法容許在一個網格中包含多種物質材料，通過跟蹤每種材料的邊界，在相應的單元中進行物質交換和輸送。本文利用多物質ALE方法處理充氣床墊內外空氣與床墊薄膜間的相互作用問題。流體與結構體的相互作用通過耦合來定義，耦合方式有多種，一般選擇罰函數耦合方式［10］。

2 問題描述

某型號充氣床墊，由新型柔性薄膜材料制成。具有重量輕、可折疊、方便攜帶等特點，適用于遠距離長期野外活動。該充氣床墊分為床墊和枕頭兩大部分，外部由床墊面膜封閉，內部由多種設計的拉筋結構與面膜相連，可分擔一定的面膜拉力。床墊外部面膜為植絨TPU薄膜材料，厚0.2 mm。內部拉筋材料為TPU薄膜，厚0.3 mm。

本次計算包括兩個部分:第一是針對25 kPa充氣內壓時的爆破工況，采用空氣-床墊耦合的多物質ALE方法進行模擬，并結合試驗結果，分析充氣床墊失效原因。同時提出改進方案，并作進一步仿真計算。該部分計算采用16個CPU，計算需時約10小時。第二是針對載人工作工況，模擬假人-床墊-地面相互作用下床墊的變形與受力狀況。此時假人、床墊和地面通過接觸的方式相互作用，空氣和床墊通過多物質ALE耦合方法相互作用，參見圖4。該工況在Optimus優化分析平臺聯合Ls-dyna仿真，通過參數分析的方法，重點分析不同內壓下床墊的變形和應力情況，并判斷床墊最佳工作內壓。該部分計算采用64個CPU，計算需時約136小時。

3 計算模型

3.1 有限元模型

根據充氣床墊設計資料，采用三維建模方法建立了充氣床墊三維有限元模型。由于結構均為柔性薄膜材料構成，因此床墊面膜與內部拉筋均采用膜單元模擬。如圖2為充氣床墊三維有限元模型，該模型單元數20 956，節點數19 509。圖3為三維有限元假人模型，假人模型來源于三維人體造型軟件Poser。通過模型轉換，將Poser中的幾何人體模型轉換為Hypermesh中的有限元模型，從而引入Ls-dyna中進行計算。由于本文研究的重點為充氣床墊結構載人工況時的應力變形情況，假人模型只是作為一個載荷施加于床墊結構，對假人本身不作分析，因此只需保證假人正確的平臥姿態即可。模型在Poser中首先進行姿態調整，使人體模型為通常的平臥睡姿，從而保證了假人模型的合理性。有限元假人模型采用三節點殼單元模擬，材料為剛體材料。單元數83 991，節點數42 381;假人高1.75 m，重80 kg。

圖2 充氣床墊三維有限元模型Fig.2 3-D finite element model of air-charge mattress

圖3 三維有限元假人模型Fig.3 3-D finite element model of dummy

采用ALE方法模擬空氣床墊的充氣和工作過程，空氣與床墊通過耦合方式實現交互作用;同時假人、床墊和地板三者間建立了基于罰函數的接觸模型［11］。圖5為系統整體有限元模型，與床墊耦合的空氣分為內外層空氣，其中內層空氣為高壓空氣，外層空氣為外界大氣。內層空氣壓力可根據不同工況調整，外層空氣在周圍邊界施加無反射邊界模擬無限遠空氣特征?？諝鈫卧獢?27 200，節點數239 112。

圖4 系統整體有限元模型Fig.4 Global finite element model of system

3.2 材料模型

床墊薄膜為TPU植絨復合材料，拉筋采用TPU材料，由于TPU材料具有復雜的力學特性，通過簡單拉伸試驗測定其力學參數。拉伸試樣尺寸為200mm×50mm，經試驗測定TPU植絨復合材料具有明顯的各向異性，且在經緯向分別呈近似線彈性特性，本文選用Ls-dyna各向異性彈性材料模型(*MAT_ORTHOTROPIC_ELASTIC)模擬床墊面膜材料。密度 ρ=1 200 kg/m3，徑向彈性模量Ex=2.77e2 MPa，緯向彈性模量Ey=1.45e2 MPa，泊松比 υ =0.48。

圖5 TPU材料應力應變關系Fig.5 Stress vs strain curve of TPU

試驗測定TPU材料為各向同性材料，但具有明顯的非線性力學特征，為準確模擬其非線性特征，本文選用Ls-dyna分段彈塑性材料模型(*MAT_PIECEWISE_LINEAR_PLASTICITY)模擬拉筋材料。密度ρ=1 200 kg/m3，泊松比υ=0.48，應力應變關系直接按實驗所得曲線模擬，圖5描述了TPU材料的應力應變關系。

空氣采用空材料模型(Null)。該模型在求解中不計算材料的應力和應變，而只是求解材料的運動情況。內外層空氣壓力采用Ideal Gas狀態方程計算，具體形式為:

式中:P為壓力，γ為比熱比，Ev，0為初始內能，vr為比體積。氣囊內空氣密度ρ=1.18 kg/m3，體積模量Mu=0.14 MPa，γ =1.4，Ev，0=63 328 J/m3，vr=1.0。外部空氣密度 ρ=1.18 kg/m3，體積模量Mu=0.14 MPa，γ =1.4，Ev，0=0.0 J/m3，vr=1.0。

4 計算結果與試驗驗證

4.1 25 kPa內壓爆破工況分析

(1)結構薄弱位置分析

圖6～圖7分別為爆破工況下充氣床墊內部空氣變形狀況，以及不同時刻充氣床墊的應力分布云圖?？芍?.1 s時充氣床墊變形與應力分布較均勻，最大應力出現在床墊側面四角位置，1.0 s時充氣床墊結構出現形狀突變和應力集中，并且有進一步增大的趨勢。最大變形和應力均出現在縱拉筋與橫拉筋相交的床墊面膜位置。此時床墊面膜出現嚴重隆起狀畸變并出現應力集中，初步判斷結構已發生失穩，因此可認為此處為爆破工況薄弱位置。與爆破試驗對比分析，該結論與試驗結果一致，從而證明了仿真分析的可靠性。

圖6 充氣床墊內部空氣變形示意圖Fig.6 Deformation of interior air in air-charge mattress

圖7 充氣床墊應力分布Fig.7 Stress distribution of air-charge mattress

(2)結構優化改進

參考圖2充氣床墊內部結構模型，分析充氣床墊出現局部畸變的原因，由于縱拉筋和橫拉筋相交處有多種拉筋結構存在，且臨近枕頭氣囊與床墊氣囊的連接位置，結構突變明顯，因此極容易出現應力集中并造成局部畸變。在上述仿真的基礎上，對床墊結構進行局部修改，如圖8所示，為強化縱橫拉筋相交位置，將充氣床墊縱橫拉筋間的間距縮小一半，其它結構保持不變。

圖9為局部結構改進后爆破工況下充氣床墊不同時刻的應力分布云圖，從圖中可知，不同時刻充氣床墊變形與應力都分布均勻，最大應力出現在床墊側面四角位置?？v拉筋與橫拉筋相交的床墊面膜位置不再出現應力集中和局部隆起狀畸變。說明對床墊縱橫拉筋的局部結構改進，很好的改善了床墊該位置的受力特性，提高了床墊整體的力學性能。

圖8 充氣床墊結構改進Fig.8 Structure improvement of air-charge mattress

圖9 結構改進后充氣床墊應力分布Fig.9 Stress distribution of improved air-charge mattress

4.2 載人工況分析

本次計算的主要目的是通過分析得出床墊載人工況下的最佳工作內壓，即在該內壓條件下，床墊的變形和應力最小。設定的約束條件為床墊充氣內壓，范圍為2～8 kPa。設定的目標條件為最小的變形和應力狀況。在Optimus優化平臺下聯合Ls-dyna仿真，通過參數分析方法進行分析，計算定義了50組設計實驗，每組實驗對應不同的充氣內壓，通過對比分析得出載人工況下的最佳充氣內壓。圖10給出了參數分析的計算流程。

圖10 參數分析的計算流程Fig.10 The calculation flow of parameter analysis

圖11～圖12分別給出了工作工況下假人、床墊和地板三者相互作用示意，以及床墊的變形情況?？芍矇|在與假人臀部及背部接觸的位置發生較大變形，并且出現局部應力集中現象。

圖11 工作工況下假人、床墊和地板相互作用Fig.11 Interactive action between dummy，mattress and ground in working condition

圖12 工作工況下床墊變形和應力分布Fig.12 Stress and Displacement distribution of air-charger mattress in working condition

圖13 床墊最大應力隨充氣內壓的變化曲線Fig.13 The curve of maximal stress versus pressure

圖13～圖14分別給出了載人工作工況下床墊最大應力隨充氣內壓的變化曲線，以及床墊最大豎向位移隨充氣內壓的變化曲線?？芍?，當充氣內壓為4.3 kPa左右時，床墊結構的最大應力處于最小狀態。當充氣內壓小于4.3 kPa時，床墊最大豎向位移隨內壓的增加顯著減小;當充氣內壓大于4.3 kPa時，床墊最大豎向變形隨內壓的增加減小量不明顯。綜合分析，床墊充氣內壓為4.3 kPa左右時，應力水平最低，變形也處于較低水平。較小的內壓將導致床墊變形過大，并產生過大的應力。較大的內壓對抑制床墊變形效果有限，而床墊的應力卻較大。因此，可認為床墊最佳工作氣壓為4.3 kPa左右。

圖14 床墊最大豎向變形隨充氣內壓的變化曲線Fig.14 The curve of maximal vertical displacement versus pressure

5 結論

(1)采用多物質ALE方法對某型薄膜結構充氣床墊進行了25 kPa內壓下爆破工況的數值仿真。數值仿真通過與物理試驗對比驗證，判斷出縱拉筋與橫拉筋相交的床墊面膜位置為薄弱位置，并對該局部結構進行改進，經過進一步的數值計算，證明結構改進很好的改善了床墊整體力學性能。

(2)采用多物質ALE方法對某型薄膜結構充氣床墊進行了載人工況的數值模擬。在Optimus優化平臺下聯合Ls-dyna仿真，通過參數分析方法進行不同設計實驗下床墊應力和變形的對比分析，結果表明4.3 kPa左右充氣內壓為最佳工作內壓。最終的數值結果可作為該產品設計定型的重要依據，并為同類問題提供了參考。

［1］Nicholas W B，Matthew G D.Advanced airbag landing systems for planetary landers［C］.18th AIAA Aerodynamic Decelerator Systems Technology Conference and Seminar，AIAA，2005，1615.

［2］溫金鵬，李 斌，譚德偉，等.考慮織布彈性的軟著陸氣囊緩沖特性研究［J］.振動與沖擊，2010，29(2):79-83.

［3］陳 帥，李 斌，溫金鵬，等.軟著陸氣囊緩沖特性與參數設置的理論研究［J］.振動與沖擊，2009，28(4):25-28.

［4］ Marklund P O，Nilsson L.Simulation of airbag inflation processes using a coupled fluid structure approach［J］.Computational Mechanics，2002，29:289-297.

［5］何 文，鐘志華，楊濟匡.汽車安全氣囊工作性能仿真試驗驗證技術研究［J］.機械工程學報，2002，38(4):126-129.

［6］王亞偉，楊春信，柯 鵬.貨臺空投系統氣囊緩沖過程仿真［J］.系統仿真學報，2007，19(14):3176-7179.

［7］王曉冬，張金換，黃世霖.氣囊模擬計算技術綜述［J］.設計.計算.研究，2000，9:4-7.

［8］代小芳，宗 智，王喜軍，等.折疊氣囊展開過程CV和ALE數值模擬方法對比研究［J］.汽車工程，2008，30(8):676-680.

［9］杜新光，金先龍，陳向東.基于流固耦合的新拌混凝土流動及變形仿真方法［J］.上海交通大學學報，2008，42(12):1993-1996.

［10］張曉敏，金先龍，陳向東.風荷載作用下大型柔性儲液結構動態響應仿真方法研究［J］.振動與沖擊，2009，28(5):115-118.

［11］郭毅之，金先龍，丁峻宏.沉管隧道地震響應分析中的三維接觸模型與算法研究［J］.應用力學學報，2006，23(1):48-52.