泡沫填充夾芯墻多胞結構的耐撞性多目標優化設計

第一作者文桂林男，教授，博士生導師，1970年生

通信作者尹漢鋒男，講師, 碩士生導師，1982年生

泡沫填充夾芯墻多胞結構的耐撞性多目標優化設計

文桂林1,2，孔祥正1,2，尹漢鋒1,2，肖久如2

(1.湖南大學汽車車身先進設計制造國家重點實驗室，長沙410082；2.湖南大學特種裝備先進設計技術與仿真教育部重點實驗室，長沙410082)

摘要：泡沫填充管由于其優秀的能量吸收性能引起越來越多的關注。根據對絲瓜結構的研究，提出一種與絲瓜橫截面類似的泡沫填充夾芯墻多胞結構?；诜蔷€性有限元軟件LS-DYNA，建立此種泡沫填充夾芯墻多胞結構的軸壓分析有限元模型?；诙囗検酱砟Ｐ停捎枚嗄繕肆Ｗ尤簝灮惴▽Σ煌孛嫘螤畹亩喟Y構進行了優化設計。根據優化結果對不同結構的選擇給出指導性意見，使結構在滿足峰值力限制的情況下擁有最好的吸能特性。提出的新型結構在車輛工程及航空航天等領域具有很好的應用前景。

關鍵詞：泡沫填充；絲瓜結構；耐撞性；多目標優化

基金項目：國家杰出青年科學基金(11225212);國家自然科學基金(青年項目) (11302075);國家科技支撐計劃(2012BAH09B02);高等學校博士學科點專項科研基金(20120161130001,20120161120009);汽車車身先進設計制造國家重點實驗室基金(自主課題71275003; 開放基金31275006)

收稿日期：2013-09-30修改稿收到日期：2014-03-07

中圖分類號:TG302文獻標志碼：A

Multi-objective crashworthiness optimization design of foam-filled sandwich wall multi-cell structures

WENGui-lin1,2,KONGXiang-zheng1,2,YINHan-feng1,2,XIAOJiu-ru2(1．State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body, Hunan University, Changsha 410082, China; 2．The MOE Key Laboratory of Advanced Design and Simulation Techniques for Special Equipment, Hunan University, Changsha 410082, China)

Abstract:Foam-filled tubes gain more and more attention due to their excellent energy absorption capacity. According to the study of luffa structures, a kind of foam-filled sandwich wall multi-cell structures similar to the luffa cross-section was proposed. By using the nonlinear finite element code LS-DYNA, the axial crushing finite element model of the foam-filled sandwich wall multi-cell structure was established. Based on the polynomial surrogate model, the structures with different sizes were optimized by adopting multiobjective particle swarm optimization algorithm. According to the result of optimization, suggestions for the choice of reasonable structures were given to make the selected structure have the best energy absorption characteristics under the peak force constraint. The new type of structures proposed here had a good application prospect in vehicle, aerospace and other fields.

Key words:foam-filled sandwich; luffa structure; crashworthiness; multiobjectie optimization

泡沫填充薄壁管結構由于其優秀的吸能特性和非凡的輕量化性能，在汽車、航空、軍事裝備以及其他產業擁有廣泛的應用前景[1]。泡沫填充薄壁管結構相比于相應的非泡沫填充薄壁管結構在不增加太多總質量的情況下，能吸收更多的能量[2]。因此，近年來國內外許多研究人員采用理論分析、試驗和數值仿真等方法對其吸能特性進行了研究。通過研究發現，在加載固定的情況下泡沫填充薄壁管的吸能特性主要取決于薄壁管的幾何參數和泡沫的密度[3]。

根據上述情況，有必要對泡沫填充薄壁管結構進行耐撞性優化來尋找最佳的薄壁管的幾何尺寸和泡沫的密度。Zarei等[4-5]應用多準則優化設計流程來得到泡沫填充鋁管的最大吸收能量和最小質量。Hou等[6-7]設計了一種多目標優化框架對泡沫填充薄壁管進行優化。在他們的優化過程中，使用代理模型來減少計算耗時。Bi等[8]通過約束平均壓縮力來優化單胞和三胞的泡沫填充薄壁金屬管以得到最大比吸能。Zhang等[9]使用Kriging代理模型以及遺傳算法和非支配排序遺傳算法對泡沫填充方管的耐撞性能進行優化。

本文通過對絲瓜結構的研究，設計一種新的泡沫填充夾芯墻多胞結構。使用非線性有限元軟件LS-DYNA建立不同截面形狀泡沫填充夾芯墻多胞結構的軸壓分析有限元模型。隨后，基于多項式代理模型，采用多目標粒子群優化算法對不同截面形狀的多胞結構進行了優化設計，并根據優化結果對不同結構的選擇給出指導性意見，使結構在滿足峰值力限制的情況下擁有最大比吸能。泡沫填充夾芯墻多胞結構的優化設計對類似新型結構耐撞性的設計與優化具有一定的指導意義。

1防撞性準則及結構

1.1防撞性準則

一個緩沖性能好的吸能結構首先必須滿足吸能的要求，即充分吸收物體的沖擊動能。同時吸能結構的質量也是設計人員關心的重要指標，吸能結構質量越小，相關的成本也就越低。所以吸能結構單位質量的吸能量-比吸能SEA(Special Energy Absorption)是衡量吸能結構耐撞性的重要指標，其大?。?/p>

(1)

式中:M為吸能結構的質量，EA(Energy Absorbed)為吸能結構所能吸收的能量，其大小可用下式計算：

(2)

式中:d為碰撞吸能過程中的壓縮位移，F為碰撞沖擊力。典型薄壁結構的軸向壓縮力-位移曲線如圖1所示。對于一個給定的變形，平均壓縮力MCF(Mean Crushing Force)由下式計算可得：

(3)

另外一個衡量能量吸收性能的重要指標是沖擊負載效率CLE(Crash Load Efficiency),可由下式計算得到：

(4)

式中:PCF為薄壁結構的峰值壓縮力，如圖1所示。

圖1　典型薄壁結構的軸向壓縮力—位移曲線 Fig.1 Axial crushing force-displacement curves of a typical thin-walled structure

1.2泡沫填充夾芯墻多胞結構

Xie等[10]通過對絲瓜海綿材料的剛度、強度以及吸能特性的試驗研究，認為其可以成為有多種實際用途的可持續性機械材料。目前絲瓜海綿材料受其原材料的獲取及加工的限制，還無法大量的用于實際工程設計中，但是其形狀結構可以供我們參考用于其他材料。

圖2為絲瓜結構實物圖，本文根據其結構設計出一種新的泡沫填充金屬薄壁管結構，其截面形狀如圖3所示。泡沫的截面形狀設計成與絲瓜實際截面形狀類似，外徑長度取80 mm，寬度用d表示，并且根據胞數的多少分別設計出二胞、三胞、四胞、五胞和六胞5種截面形狀。將此種截面形狀的泡沫填充在鋁薄板中間，泡沫密度用ρ表示，鋁板厚度用t表示。

圖2　絲瓜結構的實物圖(左圖：內部表面， 中圖：外部表面，右圖：橫截面) Fig.2 Physical map of luffa fibres(left:inner surface; middle:outer surface;right:cross section)

圖3　泡沫填充夾芯墻多胞結構的截面示意圖(n表示結構的胞數) Fig.3 Cross-section diagram of foam-filled sandwich wall multi-cell structures(n:the number of cell)

2泡沫填充夾芯墻多胞結構的有限元模型

2.1有限元模型

泡沫填充夾芯墻多胞結構軸向壓縮分析的有限元模型如圖4所示。結構長度取200 mm，其中一端與剛性墻接觸，另一端用質量為500 kg的剛體墻以15 m/s的初速度進行軸向加載。

本文中，采用顯式非線性有限元軟件LS-DYNA來對泡沫填充夾芯墻多胞結構的軸向壓縮進行仿真計算。鋁板結構的單元為Belytschko-Tsay-4-node薄殼單元，在厚度方向上有3個積分點，在面內有1個積分點。對于泡沫結構，采用有1個積分點的8節點實體單元?；趧偠鹊纳陈┛刂朴脕肀苊鈧瘟隳茏冃文Ｊ?，降階積分用來避免體積鎖死。

由于泡沫填充夾芯墻多胞結構由泡沫和鋁板兩種不同的材料組成，在利用LS-DYNA軟件仿真分析過程中必須使用不同的接觸算法。泡沫與鋁板之間采用可失效的固連接觸，在這種接觸定義中，如果拉伸或剪切應力的值超過了其極限強度，那么粘接將會失效接觸自動變為面面接觸，本文中拉伸與剪切的極限強度均取為150 MPa[11]。同時為了模擬結構在軸向壓縮過程中由于折疊引起的自身接觸，泡沫與鋁板分別定義自動的單面接觸。質量塊與泡沫、鋁板之間定義自動的點面接觸，靜摩擦系數和動摩擦系數均為0.2。

為了減小在數值計算中網格大小對結果精確度的影響，本文做了一系列收斂實驗來決定網格的大小。以兩胞結構為例，ρ，t和d分別取0.35 g/cm3，1.0 mm和10 mm，取鋁板的網格大小分別為1 mm×1 mm、2 mm×2 mm和3 mm×3 mm，泡沫的網格大小分別為3 mm×3 mm×3 mm、4 mm×4 mm×4 mm和5 mm×5 mm×5 mm，網格大小兩兩組合建立兩胞結構的軸向壓縮有限元模型。為了更直觀，我們只取網格最小組合、網格中間大小組合和網格最大組合得到的平均壓縮力與位移曲線，如圖5所示。從圖中可以看出鋁板網格在2 mm×2 mm、泡沫網格在4 mm×4 mm×4 mm時已經可以很好地收斂，考慮到計算效率，最后確定鋁板和泡沫的網格尺寸分別為2 mm×2 mm和4 mm×4 mm×4 mm，如圖4所示。仿真中泡沫的網格尺寸要根據寬度的變化做適當調整。

圖4　泡沫填充夾芯墻多胞結構軸壓 分析有限元模型(以2胞為例) Fig.4 Axial crushing finite element model of foam-filled sandwich wall multi-cell structures(here take n=2)

2.2材料屬性

本文中薄壁墻的材料為鋁合金AA6060 T4，其密度為ρ=2.7×103kg/m3，彈性模量E=68.2 GPa，屈服強度σy=80 MPa，極限強度σu=80 MPa，泊松比μ=0.3，且硬化指數n=0.23[12]。薄壁墻的本構模型是基于各項同性塑性法則的彈塑性材料模型。在這個本構模型中，材料的塑性硬化是通過一些分段線來定義的。為了定義這些分段線，塑性應變與對應的真實壓力值如表2所示。由于鋁合金材料對應變速率不靈敏，在這個模型中忽略應變速率的影響。

表1　材料AA 6060 T4的應變硬化數據

表2　泡沫鋁材料常數

泡沫模型是由Deshpande和Fleck提出的[13]，泡沫材料的屈服準則通過下式來定義：

Φ=σe-σy≤0

(5)

這里σy為屈服應力，等效應力σe通過下式可得：

(6)

式中：σv和σm分別表示范式有效應力和平均應力。參數α控制著屈服表面的形狀，是關于塑性泊松比vp的函數，可通過下式定義：

(7)

通過上式可以推斷出，當vp= 0時，α= 2.12。

這個材料模型的應變硬化準則可通過下式表達：

圖5　兩胞結構不同網格大小平均壓縮力與位移曲線 Fig.5 Mean crushing forces versus displacement with different mesh size of two-cell structure

(8)

式中：εe為等效應變，σp，α2，γ，β和εD為材料參數，其與泡沫的密度有以下關系：

(9)

式中:ρf和ρf0分別代表泡沫的密度和基材料的密度。C0，C1和q為常數，見表3。從式(9)可以看出，泡沫材料的彈性模量Ep也是關于泡沫密度的函數。

2.3有限元模型驗證

圖6所示為泡沫填充夾芯墻兩胞結構在軸向動態加載過程中的動能、內能、總能量和沙漏能，從圖中可以得到內能的增加等于動能的減少而總能量保持不變，沙漏能不到總能量的1%。沙漏能是衡量網格質量好壞的一個重要指標，一般情況下要控制沙漏能不到系統總能量的5%。因此，此處采用的網格尺寸能夠準確的模擬泡沫填充夾芯墻多胞結構的軸壓試驗。

圖6　系統的動能、內能、總能量及沙漏能 Fig.6 Kinetic,internal,total and hourglass energy of the system

3多目標優化設計

3.1優化方法

當泡沫填充夾芯墻多胞結構作為緩沖吸能結構時，我們希望其單位質量的吸能量越大越好，所以在耐撞性優化中，比吸能SEA將被作為耐撞性指標之一。同時，為了保證被緩沖設備或人員的安全，我們也希望泡沫填充夾芯墻多胞結構的緩沖峰值力PCF越小越好，這樣才能使得被緩沖設備或人員所承受的過載盡可能的小。所以，為了滿足這兩方面的設計要求，泡沫填充夾芯墻多胞結構的耐撞性優化問題就是典型的多目標優化問題，該多目標優化問題可以描述為：

(10)

式中:ρ，d和t分別代表填充泡沫的密度、厚度和鋁板的厚度。

本文將采用圖7所示的多目標優化設計流程求解式(10)的多目標優化問題。其中，本研究中采用的試驗設計方法為全因子法，全因子法可以全面地反映設計空間的信息。在試驗設計過程中，我們采用4水平的全因子法進行采樣，這樣對于ρ、d和t三個設計變量，總共在設計空間中均勻采集了4×4×4=64個樣本點，同時我們采用LS-DYNA計算了這些樣本點對應的目標函數值。根據這些樣本點的信息，我們采用多項式函數[14]來構建目標函數的代理模型。

(11)

相對誤差越小，代理模型越準確。此外，我們還采用均方根誤差RMSE(Root Mean Square Error)、最大絕對誤差MAX和R2值來衡量近似模型的精度，分別采用下式進行計算：

(12)

(13)

(14)

式中，k為驗證代理模型的試驗點個數，SSE和SST分別采用一下公式進行計算：

(15)

(16)

在代理模型精確足夠的基礎上，泡沫填充夾芯墻多胞結構的耐撞性多目標優化設計是可以實現的。與傳統的多目標優化算法NSGA[15]和PEAS[16-17]相比，多目標粒子群算法由于其收斂速度快、產生Pareto前沿分布均勻、算法簡單等優點，目前得到廣泛的關注。多目標粒子群算法的細節可參考文獻[16]。本文中，多目標粒子群算法被用來獲取兩個矛盾變量SEA和PCF的Pareto前沿。實際上，Patero前沿就是一系列優化方案的需求點，工程師可以根據實際需求在所得Pareto前沿中得到他們想要的優化點。

3.2泡沫填充夾芯墻多胞結構的耐撞性優化

根據多項式函數代理模型的構建方法，本文分別構建了泡沫填充夾芯墻多胞結構的SEA和PCF的一階到四階多項式函數模型。同時，根據式(11)~(16)計算了這些多項式函數代理模型的相對誤差RE、均方根誤差RSME、最大絕對誤差MAX和R2值，從得到的結果可以看出四階多項式函數的精度最高，各胞結構四

階多項式函數代理模型的精度如表4所示。求解過程中，為了求解精確將變量取值進行了歸一化處理。由于篇幅有限，本文只給出四胞結構的SEA和PCF關于變量ρ、t和d的四階多項式，其函數表達式如下：

SEA(ρ,t,d)=16.998 7+4.902 6ρ-7.836 2t+

1.591 2d+0.997 2ρ2+1.434 5ρt-3.101 5ρd+

13.908 9t2+3.287 3td-0.491 9d2-

3.391 1ρ2t+5.196 4ρ2d-6.309 9ρt2+

0.576 9ρtd+1.708 5ρd2-16.392 4t2d+

1.489 5td2+0.548 5ρ4+0.796 3ρ3t-2.327 8ρ3d-

0.708 7ρ2t2+4.334 6ρ2td-4.040 6ρ2d2+7.330 2ρt3-

7.454 7ρt2d+3.614 4ρtd2+0.441 1ρd3-6.309 8t4+

9.564 6t3d+1.926 0t2d2+0.312 3td3-1.806 7d4

PCF(ρ,t,d)=81.077 5+4.743 7ρ+19.707 3t-

5.132 1d+42.050 8ρ2-1.788 3ρt+

75.477 9ρd+9.981 1t2-0.354 6td+6.476 1d2-

29.704 7ρ2t-89.635 1ρ2d-10.594 6ρt2-

26.152 3ρtd-58.192 2ρd2-18.532 6t2d+

49.079 6td2-19.320 0ρ4+12.253 7ρ3t+

58.109 8ρ3d+20.076 5ρ2t2b-10.095 6ρ2td+

9.179 1ρ2d2+30.005 2ρt3-56.672 8ρt2d+

74.065 7ρtd2+16.280 6ρd3-0.112 6t4-5.687 6t3d+

20.251 1t2d2-50.975 4td3-3.948 8d4

圖7　多目標耐撞性優化設計流程圖 Fig.7 Flowchart of the crashworthiness multiobjective optimal design

胞數RMSEMAXR2RE(%)n=20.18090.61180.9927[-3.2990,3.7147]n=30.23650.66990.9893[-2.6700,2.8225]SEAn=40.13820.38500.9955[-2.0405,1.5425]n=50.14630.43380.9946[-2.0369,1.9325]n=60.13560.41070.9955[-1.5559,0.2320]n=21.22482.79120.9922[-2.1912,2.9138]n=31.83075.97660.9912[-4.7577,3.4719]PCFn=41.55274.10050.9917[-3.2010,3.3930]n=51.33504.89680.9937[-2.8558,3.3804]n=61.15063.23000.9961[-3.1402,1.8643]

各胞結構根據精度最高的四階多項式代理模型，采用多目標粒子群算法對式(10)的多目標優化問題進行計算。在采用粒子群算法優化過程中，慣性權重c0=0.4，加速因子c1=c2=1.0，種群大小p=100，解集合A數量為50。當迭代次數取100時Pareto前沿已經收斂穩定，Pareto前沿即為式(10)描述的多目標優化問題的解，各胞結構的Pareto前沿如圖8所示。從圖中可以看出，相同比吸能所對應的峰值力隨著泡沫填充夾芯墻多胞結構的胞數增加而增加；而在所考慮變量范圍內，最大比吸能先隨著胞數的增加而增加，在四胞結構時比吸能達到最大，隨后隨著胞數的增加，比吸能會變小，但五胞結構和六胞結構的比吸能變化范圍差別不是很大，說明四胞之后再增加結構的胞數對結構的比吸能影響不是很大。在所考慮的五種結構中，若不考慮峰值力的限制，四胞結構可取得最大比吸能，如圖8中點a所示，優化設計的參數如表5工況a所示；若實際問題中要求緩沖結構的PCF值必須小于105 kN，要想得到最大比吸能，最優設計點選擇圖8中的b點，此時為三胞結構，優化設計的參數如表5工況b所示；若實際問題中要求緩沖結構的PCF值必須小于90 kN，要想得到最大比吸能，最優設計點選擇圖8中的c點，此時為兩胞結構，優化設計的參數如表4工況c所示。

圖8　泡沫填充夾芯墻多胞結構的的Pareto前沿 Fig.8 Pareto front of foam-filled sandwich wall multi-cell structures

分別建立這三種優化設計點對應的軸向壓縮有限元模型，采用LS-DYNA計算得到SEA和PCF值如表4所示。從表中可以看出，有限元仿真結果與代理模型的誤差小于3%，說明代理模型的精度很高可以用于實際工程設計。圖9為各優化工況對應的峰值壓縮力與位移曲線，圖中可以看出優化工況b和工況c的最大峰值力非常接近105 kN和90 kN，但是仍然小于其對應值。另外在出現最大峰值力之后，壓縮力趨于平緩而且靠近最大峰值力。圖10為各優化工況結構60%壓縮時的變形圖，從圖中可以看出三種工況對應結構在壓縮過程中折疊形狀非常有規律，這說明泡沫填充夾芯墻多胞結構具有良好的碰撞穩定性，如果作為能量吸收器將擁有非常吸引人的吸能特性。

圖9　各優化設計工況對應峰值壓縮力與位移曲線 Fig.9 Forces-displacement curves of the different optimal design conditions

圖10　各優化設計工況結構60%壓縮時變形圖 Fig.10 60% deformation figures of the different optimal design conditions

工況ρ/(g·cm-3)t/mmd/mmSEA/(kJ·kg-1)PCF/kN近似模型仿真誤差/%近似模型仿真誤差/%a0.40.87.957623.610323.8645-1.0652115.7285116.9200-1.0191b0.40.87.397623.259622.82521.9032104.6009103.69000.8776c0.38470.8039621.396021.12451.285288.964588.26300.7948

4結論

根據絲瓜結構，設計了一種與絲瓜橫截面類似的泡沫填充夾芯墻多胞結構，有二胞、三胞、四胞、五胞和六胞五種截面形狀；基于非線性有限元軟件LS-DYNA，建立了不同胞數結構軸壓分析的有限元模型，并對模型網格大小進行了收斂性驗證。

基于多項式代理模型，采用粒子群優化算法對不同截面形狀的多胞結構進行了優化設計，并根據優化結果對不同結構的選擇給出指導性意見。在考慮變量范圍內，若無峰值力限制取四胞結構可得最大比吸能；若峰值力要小于105 kN或90 kN，則分別取三胞結構或兩胞結構可得最大比吸能。以此建立的有限元模型與代理模型結果的誤差在3%以內，說明基于多項式代理模型的泡沫填充夾芯墻多胞結構的耐撞性優化是可行且有效的，可指導類似薄壁管結構的設計與優化。

參考文獻

[1]Ahmad Z, Thambiratnam D P. Application of foam-filled conical tubes in enhancing the crashworthiness performance of vehicle protective structures[J]. Int J Crashworthines, 2009, 14(4): 349-363.

[2]Seitzberger M, Rammerstorfer F G, Degiseher HP, et al. Crushing of axially compressed steel tubes filled with aluminium foam[J]. Acta Mech,1997,125(1-4): 93-105.

[3]Yin H F, Wen G L, Hou S J, et al. Multiobjective crashworthiness optimization of functionally lateral graded foam-filled tubes[J]. Materials and Design,2013,44: 414-428.

[4]Zarei H R, Kr?ger M. Optimization of the foam-filled aluminum tubes for crush box application[J]. Thin Wall Structure,2008,46(2): 214-221.

[5]Zarei H R, Kr?ger M. Crashworthiness optimization of empty and filled aluminum crash boxes[J]. Int J Crashworthiness,2007, 12(3): 255-264.

[6]Hou S J, Li Q, Long S Y,et al. Crashworthiness design for foam filled thin-wall structures[J]. Mater Des, 2009,30(6): 2024-2032.

[7]Hou S J, Han X, Sun G Y,et al.Multiobjective optimization for tapered circular tubes[J]. Thin Wall Struct,2011, 49(7): 855-863.

[8]Bi J, Fang H B, Wang Q,et al. Modeling and optimization of foam-filled thin-walled columns for crashworthiness designs[J]. Finite Elem Anal Des,2010,46(9): 698-709.

[9]Zhang Y, Sun G Y, Li G Y,et al. Optimization of foam-filled bitubal structures for crashworthiness criteria[J]. Mater Des,2012,38: 99-109.

[10]Shen J H, Xie Y M, Huang X D, et al. Mechanical properties of luffa sponge[J]. Journal of the Mechanical Behavior of Biomedical Materials,2012,15: 141-152.

[11]Santosa S, Wierzbicki T. Crash behavior of box columns filled with aluminium honeycomb or foam[J]. Comput Struct,1998, 68(4): 343-367.

[12]Santosa S P, Wierzbicki T, Hanssen A G,et al. Experiment and numerical studies of foam-filled sections[J]. Int J Impact Eng,2000,24(5): 509-534.

[13]Deshpande V S, Fleck N A. Isotropic constitutive models for metallic foams[J]. J Mech Phys Solids,2000,48(6-7): 1253-1283.

[14]Hou S J, Li Q, Long S Y, et al. Multiobjective optimization of multi-cell sections for the crashworthiness design[J]. International Journal of Impact Engineering,2007,35(2008): 1355-1367.

[15]Liao X T, Li Q, Yang X J, et al. Multiobjective optimization for crash safety design of vehicles using stepwise regression model[J]. Structural and Multidisciplinary Optimization,2008,35(6):561-569.

[16]Raquel C, Naval P. An effective use of crowding distance in multiobjective particle swarm optimization[C]//Proceedings of the 2005 conference on Genetic and evolutionary computation, Washington(DC), USA, 2005:257-264.

[17]Liu D S, Tan C K, Goh C K. A multiobjective memetic algorithm based on particle swarm optimization[J]. IEEE Transactions on Systems Man and Cybernetics Part B-Cybernetics,2007,37(1): 42-50.

[18]白金澤．LS-DYNA3D理論基礎與實例分析[M].北京：科學出版社，2005.