殼單元和實體單元模擬爆炸荷載作用鋼板動力響應的差異性比較

朱俊杰，盧永剛，馮星奎，陶俊林

(1. 西南科技大學 土木工程與建筑學院，四川 綿陽 621010；2. 中國工程物理研究院 總體工程研究所，四川 綿陽 621000；3. 火箭軍駐中國工程物理研究院 型號辦事處，四川 綿陽 621900)

殼單元和實體單元模擬爆炸荷載作用鋼板動力響應的差異性比較

朱俊杰1，盧永剛2，馮星奎3，陶俊林1

(1. 西南科技大學 土木工程與建筑學院，四川 綿陽 621010；2. 中國工程物理研究院 總體工程研究所，四川 綿陽 621000；3. 火箭軍駐中國工程物理研究院 型號辦事處，四川 綿陽 621900)

以艦船艙室鋼板為研究對象，利用動力分析有限元程序 Ls-dyan，分別采用殼單元和實體單元 2 種建模方式對爆炸荷載作用下鋼板的動力響應進行數值模擬研究。對比鋼板的應力分布云圖、中線上各點最大位移曲線圖，發現 2 種方法計算結果與試驗值大致吻合，采用殼單元建模在保證計算結果準確性的前提下可以大大提高計算效率。如何正確地選擇單元進行分析是數值模擬很重要的一個環節，通過對殼單元和實體單元計算結果的初步對比，為艦船艙室數值模擬研究如何選取單元提供參考。

船用鋼板；爆炸荷載；數值模擬；殼單元

0 引 言

艦艇是由縱橫構架和殼板結構組成，其中殼板主要承受水平荷載。艦船艙室結構承受的主要載荷之一是戰斗部爆炸產生的沖擊波載荷[1]，鋼板在受到具有傳播速度快、峰值荷載大、作用時間短等特征的爆炸沖擊波載荷作用時，其受力特性和動態響應顯著異于靜載情況下，且十分復雜，故而研究鋼板在爆炸荷載作用下的動力響應很有必要，是研究艙室結構的基礎。目前國內對艦船艙室方面的研究大多集中在艙室結構的仿真分析上，朱建方等[2]對艦艇艙室在爆炸荷載下的動態響應過程進行了數值模擬，得到了艦艇艙室在爆炸沖擊波作用下的響應特性和破壞機理。鄢順偉等[3]仿真計算了不同壁厚艙室結構在反艦導彈戰斗部內爆作用下的毀傷效應。但是，對船用鋼板在爆炸荷載作用下的響應研究較少。本文采用數值模擬的方法對爆炸載荷下四邊固支的鋼板進行研究。在利用數值分析研究時，單元類型的選取對計算效率有很大的影響。鋼板屬于板結構，實體單元與殼單元均可以對其進行模擬，但是 2 種建模方法得到結果的差異值得探討。通過分析兩者的差別，可以更好地為鋼板使用何種單元在保證其結果準確性的基礎上又能提高計算效率提供參考。為此，本文以文獻中試驗為基礎[4–5]，在驗證數值計算方法可靠的基礎上，分別采用殼單元與實體單元對四邊固支的方形板在爆炸荷載作用下的動力響應進行數值模擬研究，考察使用不同的建模單元類型對計算結果的影響[6]，并對結果的差異性進行了比較與分析。

1 Shell163 單元介紹

Shell163 是一個具有 4 節點單元，有彎曲和膜特征，可施加平面和法向載荷。單元在每個節點上有 12個自由度，分別為在節點 x，y 和 z 方向的平動、加速度、速度和繞 x，y 和 z 軸的轉動。Shell163 單元一共有 12 種算法，這里選第 10 種算法，對于翹曲情況可以得到正確結果。和實體單元一樣，積分點的個數直接影響著 CPU 的時間。典型地，對于彈性材料沿厚度方向需要 2 個積分點，而對于塑性材料則需要 3 個或更多的積分點[7]。

2 數值計算模型

2.1 算法選取

分析主要涉及到流體流動、固體大變形以及流體與固體相互作用等問題，Ls-dyna 程序中提供了 Lagrange，Euler 和 ALE（任意拉格郎日歐拉算法）3 種不同的算法，其中 ALE 算法兼具 Lagrange 算法和 Euler算法的優點，可以描述流體與固體的相互耦合作用；同時在求解過程中網格還可以根據定義的參數適當調整位置，克服網格嚴重畸變的問題[8]。因此，本文中選取 Ls-dyna 程序提供的 ALE 算法建立流體與固體間的聯系，實現流固耦合動態分析，模擬方形鋼板在爆炸荷載作用下的響應和破壞。

2.2 材料本構與狀態方程

結構材料在爆炸沖擊載荷作用下會產生巨大的變形，與受靜載荷作用的情況相比，動態屈服應力比靜態屈服應力有明顯提高，材料力學性能表現相對也更為復雜。固支方板采用低碳鋼，鋼板選用與應變率相關的各向同性塑性隨動硬化模型，材料參數見表 1。

表 1 鋼板材料參數Tab. 1 Steel plate material parameters

炸藥在空氣中爆炸時，會在極短的時間內把化學能全部轉變為熱能，導致爆轟產物向周圍迅速膨脹而做功，形成空氣沖擊波。文中采用 *HIGH_EXPLOSIVE_BURN 模型作為炸藥的材料模型，采用 Jones-Wilkins-Lee（JWL）狀態方程來描述，即

式中：A ，B，R1，R2和 ω 為 JWL 狀態方程參數；E為炸藥的內能；V 為當前相對體積。

對空氣采用 *MAT_NULL 材料模型以及線性多項式 *EOS_LINEAR_POLYNOMI 狀態方程加以描述，即

式中：C0～C6為常量；e 為初始體積內能；μ = 1/V0–1；V0為相對體積；C1～C4均為 0。

2.3 有限元模型

1） 空氣網格的驗證

在數值模擬中，網格的劃分至關重要，網格的大小不僅關系著數值模擬計算的精度，也關系著數值模擬計算所用的時間。一般而言，網格越小則計算的結果越接近于真實情況，但隨著網格的減少計算時間也急劇增大，所以選擇合理的網格尺寸很有必要。網格驗證幾何模型尺寸為空氣 120 mm × 120 mm × 60 mm，炸藥為 27 mm × 27 mm × 27 mm。利用對稱性取 1/4 模型進行構建，模型如圖 2 所示。

表 2 炸藥與狀態方程參數Tab. 2 Explosive and equation of state parameters

表 3 多項式狀態方程參數Tab. 3 Polynomial equation of state parameters

空氣與炸藥均采用 8 節點實體單元，TNT 當量為 20 g，炸藥中心位置起爆。在對稱面上施加對稱約束，四周（除對稱面）施加固定邊界條件。為了模擬出無限空氣域，空氣模型邊界施加透射邊界。工況如表 4 所示。

表 4 空氣網格驗證工況表Tab. 4 Air grid authentication working table

上述工況在距離炸藥中心 6 cm 和 8 cm 位置處的超壓如圖 3～圖 4 所示。

由圖 5～圖6 可以看出，隨著網格的減少，壓力的峰值越來越大，當網格取 0.3 cm 和 0.2 cm 時壓力峰值的差值約為 5%，考慮到計算精度和時間的因素空氣網格取 0.3 cm。

2） 炸藥網格的驗證

根據上面網格確定炸藥網格的驗證工況如表 5 所示。

表 5 炸藥網格驗證工況表Tab. 5 Explosive grid authentication working table

由圖 5～圖 6 可知，隨著炸藥網格的減少，壓力峰值有所增加，當網格取為 0.2 cm 和 0.1 cm 時的壓力峰值差值約為 2%，炸藥網格選取 0.2 cm 較為合適。

3） 2 種單元計算模型

采用 Ansys/Ls-dyna 建立有限元模型，對于固支鋼板在爆炸載荷的作用過程進行了數值模擬研究。固支鋼板尺寸為 63 mm × 63 mm，厚度為 1.6 mm，建立 1/4模型。方板尺寸如圖 7 所示。

數值模型由炸藥、空氣、鋼板 3 部分組成，采用殼單元和實體單元分別建立板模型，2 種模擬方式除了方板的單元類型不同，其他參數均一致。空氣單元尺寸為 3 mm，炸藥單元尺寸為 2 mm，板單元為 0.4 mm，有限元模型如圖 8～圖9 所示。

2.4 失效準則

采用 *MAT_ADD_EROSION 對鋼板材料失效準則進行定義，塑性失效應變 εmax= 0.30，考慮到在沖擊荷載作用下薄壁材料容易出現斷裂、絕熱剪切等破壞情況，在計算模型中增加最大剪應變失效準則，剪切失效應變 τmax= 0.20[9]。利用單元刪除技術來模擬鋼板的破壞情況，即當單元變形引起的塑性應變或剪切應變超過臨界值時，就認為該單元破壞。

3 殼單元與實體單元模擬結果比較

3.1 方形鋼板的動態響應

1）t = 0.15 ms 時刻，方形鋼板 Von Mises 應力云圖對比。

從圖 10 可知，無論是采用殼單元還是實體單元，方形板中的應力分布大致相同，最大值均出現在炸藥最上方中心位置處。同一時刻，實體單元中最大應力值為 503.6 MPa，殼單元中最大應力值為 504.7 MPa。

從圖 11 可知，方形鋼板的變形大致呈現為中心處的外鼓，呈“金字塔”狀，開始先發生中心處的翹曲，隨著時間的增加，翹曲越來越嚴重，而后慢慢下降，2 種方法模擬的鋼板變形發展趨勢與試驗結果變形趨勢相同。

2）鋼板中線上各點處最大位移值對比

圖 12給出了在 20 g 炸藥量產生的爆炸荷載作用下鋼板中線上各點的最大位移曲線。由圖中可以看出，2種方法得到的各點處最大位移曲線基本一致，與試驗值較為吻合，各點處最大位移差值均在 10% 以內。最大的位移發生在鋼板中心位置處，2 種模擬方法得到的結果與試驗值差值均在 5% 以內，其他位置處的差值較小。

3.2 兩種方法計算時間對比

2 種模型的炸藥和空氣域單元數目相同，而使用殼單元建模使得鋼板單元數目減少 50%。由圖 13 可以看出使用殼單元建模計算時間降低約 50%，顯著提高計算效率。尤其對于單元數目比較多的模型，殼單元的優勢更為突出。

無論是從鋼板的動力響應歷程還是鋼板的變形特征層面上的對比驗證，其結構均表明：本次數值模擬的建模技術和參數選取合理、可靠，其數值模擬的結果可以艦船艙室內爆選取單元提供依據。

4 結 語

對方形鋼板在爆炸荷載作用下的響應，分別采用殼單元和實體單元進行模擬。對比 2 種方法中板的應力分布和中線處各點的最大位移值，得到的結果大致相同，差值在 10% 以內。結果表明使用實體單元建模和使用殼單元建模對于等厚度方板的計算結果影響不大。而采用殼單元建模相比采用實體單元建模可以顯著提高計算效率。

數值分析方法在爆炸荷載作用下艦船艙室的動力響應研究中有著廣泛的應用，但如何正確地選擇單元進行分析是數值模擬中很重要的一個環節。本文通過對殼單元和實體單元模擬方形鋼板計算結果的初步對比，為艦船艙室模型中如何選擇單元提供參考。

[1]陳攀, 劉志忠. 艙室內爆沖擊波載荷特性及影響因素[J]. 艦船科學技術, 2016, 38(2): 43–48. CHEN Pan, LIU Zhi-zhong. Research on loading of explosive and influencing factors inside closed cabin [J]. Ship Science and Technology, 2016, 38(2): 43–48.

[2]朱建方, 王偉力, 曾亮. 艦艇艙室內爆毀傷的建模與仿真分析[J]. 系統仿真學報, 2009, 21(22): 7066–7068.

[3]杜志鵬, 李曉彬, 夏利娟, 等. 反艦導彈攻擊艦船舷側防護結構過程數值仿真[J]. 哈爾濱工程大學學報, 2006. 27(4): 484–487.

[4]LONGDON G S, CHUNG Ki Y S, NURICK G N. Experimental and numercial studiesion the response of quadrangular stiffened plates. Part I: subjected to uniform blast load [J]. In-ternational Journal of Impact Engineering. 2005, (31): 55–83.

[5]LONGDON G S, CHUNG Ki Y S, NURICK G N. Experimental and numercial studiesion the response of quadrangular stiffened plates. Part II: localist blast loading [J]. International Journal of Impact Engineering. 2005, (31): 85–11.

[6]徐新光. 梁單元和實體單元模擬框架結構的差異性比較[J].四川建材, 2014, 40(5): 46–48.

[7]ANSYS工程結構數值分析[M]. 北京. 人民交通出版社, 2012: 11–13.

[8]ANSYS/LSDYNA在爆炸與沖擊領域內的工程應用[M]. 北京:中國建筑工業出版社, 2011.

[9]李帆. 船艦室內爆炸破壞的數值模擬和試驗研究[D]. 北京: 國防科技大學, 2012.

Comparision of differentiation with the steel plates' dynamic response subjected to blast loads using shell element and solid element modeling methods

ZHU Jun-jie1, LU Yong-gang2, FENG Xing-kui3, TAO Jun-lin1
(1. School of Civil Engineering and Architecture, Southwest University of Science and Technology, Mianyang 621010, China; 2. Institute of Structural Mechanics,China Academy of Engineering Physics, Mianyang 621900, China; 3. The Office of Rocket Force, China Academy of Eugineering Physicies, Mianyang 621900, China)

Take the steel plates of warship cabins as research objects, numerical simulation study on the response of the quadrangular steel plates subjected to blast loads is presented by dynamic analysis finite element procedure Ls-dyan3d, using modeling methods of shell element and solid element respectively. Comparing pressure distribution cloud picture and each point of midcourt line miximum displacement curve graph, it can be found that transformation trend of steel plates is similar for different thickness while each point of midcourt line miximum displacement is discriminating. However results from two methods are alike for equal thickness. Meanwhile, it can improve calculative effciency. How to choose which kind of element types correctly to analysis is an important link. Through the preliminary contrast study of calculative results of shell element and solid element, it can provide reference for choosing element type in warship cabins numerical simulation.

marine steel plate；blast loads；numerical simulation；shell element

TJ03

A

1672–7619(2017)03–0018–05

10.3404/j.issn.1672–7619.2017.03.004

2016–09–18；

2016–10–27

國家自然科學基金資助項目(11502258)

朱俊杰(1991–)，女，碩士研究生，研究方向為結構沖擊與爆炸。

陶俊林(1972–)，男，博士，教授，E-mail:junlintao@126.com。