水下爆炸船體結(jié)構(gòu)響應(yīng)間斷伽遼金法數(shù)值模擬

于福臨，郭 君，姚熊亮，任少飛

(哈爾濱工程大學(xué) 船舶工程學(xué)院, 150001 哈爾濱)

?

水下爆炸船體結(jié)構(gòu)響應(yīng)間斷伽遼金法數(shù)值模擬

于福臨，郭 君，姚熊亮，任少飛

(哈爾濱工程大學(xué) 船舶工程學(xué)院, 150001 哈爾濱)

為求解水下爆炸強(qiáng)間斷流場(chǎng)，采用Level Set方法定位多相流界面位置，應(yīng)用虛擬流體方法處理鄰近界面兩側(cè)物理量，并用RKDG方法進(jìn)行空間和時(shí)間的離散，求解流場(chǎng)的Euler方程，并進(jìn)行一維、二維評(píng)價(jià)，計(jì)算結(jié)果能夠較好地反映水下爆炸沖擊波產(chǎn)生、傳播、反射和爆炸產(chǎn)物的膨脹等現(xiàn)象.最后，結(jié)合大型非線(xiàn)性有限元軟件ABAQUS，模擬了船體板在水下爆炸載荷作用下的變形和響應(yīng)特征.模擬結(jié)果表明,間斷迦遼金法能夠?qū)崿F(xiàn)對(duì)水下爆炸船體結(jié)構(gòu)響應(yīng)精確模擬，板架結(jié)構(gòu)響應(yīng)與爆心距離成反比.

水下爆炸;間斷伽遼金;載荷計(jì)算;結(jié)構(gòu)響應(yīng);有限元;數(shù)值模擬

水面艦船在作戰(zhàn)中，會(huì)受到水雷、魚(yú)雷等水下武器的攻擊，對(duì)艦船生命力構(gòu)成威脅.國(guó)內(nèi)外學(xué)者對(duì)水下爆炸載荷及結(jié)構(gòu)響應(yīng)方面開(kāi)展了廣泛的研究.試驗(yàn)研究是最有效的手段，國(guó)內(nèi)外開(kāi)展了一些實(shí)船[1]、艙段以及模型試驗(yàn)[2]，但水下爆炸試驗(yàn)耗資巨大.解析法只能對(duì)簡(jiǎn)單邊界的水下爆炸載荷進(jìn)行求解，文獻(xiàn)[3-4]給出了水下爆炸載荷的經(jīng)驗(yàn)公式，但是經(jīng)驗(yàn)公式無(wú)法對(duì)過(guò)程進(jìn)行描述，而且難以對(duì)結(jié)構(gòu)響應(yīng)進(jìn)行求解.隨著計(jì)算機(jī)性能的提高，數(shù)值技術(shù)成為研究水下爆炸的一種有效手段，發(fā)展出了有限元、有限體積、邊界元和無(wú)網(wǎng)格等數(shù)值方法.間斷伽遼金方法是介于有限元法和有限體積法之間的一種算法[5]，具備兩者的優(yōu)點(diǎn).一方面，本質(zhì)上是一種高精度的微分方程空間離散方法，在離散過(guò)程中使用弱形式積分，具有有限維解，在單元內(nèi)使用多項(xiàng)式逼近物理量的真實(shí)值，通過(guò)簡(jiǎn)單的增加多項(xiàng)式次數(shù)，即可直接提高單元階數(shù)，構(gòu)造高階格式，與有限元方法相似;另一方面，間斷伽遼金方法引入數(shù)值通量的概念，在單元與單元間計(jì)入流場(chǎng)間斷，使其可用于沖擊波等強(qiáng)間斷現(xiàn)象的捕捉，同時(shí)還具備多種非線(xiàn)性的限制器，在保證其精度的前提下，抑制非物理振蕩，這與有限體積法相似.間斷伽遼金法本身就具有高分辨率捕捉的特性，能有效處理水下爆炸中出現(xiàn)的壓力、密度、速度等物理量的間斷.

本文基于間斷伽遼金方法，結(jié)合虛擬流體方法(ghost fluid method)和Level Set方法，編制計(jì)算程序，對(duì)一維和二維算例進(jìn)行了計(jì)算，并模擬了近壁面條件下的水下爆炸流場(chǎng)特性，最后結(jié)合ABAQUS，對(duì)水下爆炸載荷作用下船體結(jié)構(gòu)的變形和響應(yīng)特征進(jìn)行了模擬.

1 數(shù)值模擬方法

1.1 控制方程

對(duì)于無(wú)黏可壓縮流場(chǎng)，歐拉方程組可表示[6]為

(1)

對(duì)于二維計(jì)算區(qū)域，式(1)可表示為

(2)

式中：ρ為流體密度；p為壓力； u為沿x軸速度；v為沿y軸速度；E為單位體積總能量.

1.2 狀態(tài)方程

氣體采用理想氣體狀態(tài)方程描述[7]

(3)

對(duì)于水而言，采用Tait狀態(tài)方程描述[8]為

(4)

式中：ρ為密度；B取3.31×108Pa；γ=7.15；A=105Pa；ρ0=1 000 kg/m3.

1.3 空間與時(shí)間離散

對(duì)于水下爆炸無(wú)黏可壓縮性流場(chǎng)，在空間上采用間斷伽遼金法進(jìn)行離散，設(shè)流體計(jì)算區(qū)域記為Ω，在式(1)兩側(cè)同時(shí)乘以試探函數(shù)w(x)，引入流通量得[9]

(5)

(6)

式中,Pm(x)、Pn(y)分別為L(zhǎng)egendre多項(xiàng)式，將式(6)帶入式(5)即得到離散方程.

采用3階TVDRunge-Kutta方法對(duì)流體力學(xué)方程組進(jìn)行時(shí)間離散[11]，具體形式為

(7)

式中U(0)為tn時(shí)刻流場(chǎng)物理特性.

1.4 界面追蹤

應(yīng)用Level Set方法定位多相流界面位置，Level Set方程如下[12]

(8)

式中：V=(u,v,w)為運(yùn)動(dòng)速度；Φ為距離函數(shù).

(9)

式中：Φ0為初始時(shí)刻Φ值;S(Φ0)為光滑函數(shù).

1.5 虛擬流動(dòng)方法

原始的虛擬流動(dòng)方法簡(jiǎn)單的采用虛擬流體物理場(chǎng)的特性對(duì)應(yīng)真實(shí)流體的物理場(chǎng)特性，方法較為簡(jiǎn)單，在處理強(qiáng)間斷問(wèn)題時(shí)，無(wú)法給出精確的界面狀態(tài)[13]，誤差較大，甚至無(wú)法計(jì)算.

本文使用修正的真實(shí)虛擬流動(dòng)方法重構(gòu)界面狀態(tài)，使用兩相流黎曼問(wèn)題的求解方式，沿掃掠方向拾取狀態(tài)參數(shù).為了方便描述，沿x向描述這種算法.界面位于xL和xR之間，線(xiàn)段l1和l2代表界面.

界面法向n使用以下近似[14]:

(10)

設(shè)xL和xR處的狀態(tài)參數(shù)為

(11)

速度分解成法向和切向分量為

使用投影態(tài)給出的參量求解黎曼問(wèn)題為

(12)

最后，通過(guò)重新結(jié)合黎曼解得到的中間狀態(tài)和切向速度得到了界面狀態(tài)參數(shù)為

(13)

2 算例驗(yàn)證

2.1 一維算例

采用RKDG方法編制計(jì)算程序，模擬自由場(chǎng)空中爆炸流場(chǎng)，并與Henrych經(jīng)驗(yàn)公式進(jìn)行對(duì)比.表1和圖1,2分別為不同爆距條件下沖擊波超壓峰值計(jì)算值與經(jīng)驗(yàn)值對(duì)比結(jié)果.

表1 沖擊波超壓峰值計(jì)算結(jié)果

從表1可以看出，RKDG方法計(jì)算得到的沖擊波超壓峰值與經(jīng)驗(yàn)公式吻合很好，最小誤差僅為1.02%，最大誤差不超過(guò)13%，說(shuō)明本文提出的RKDG計(jì)算程序能夠有效模擬空中爆炸沖擊波載荷特性.

2.2 二維激波管算例

選取文獻(xiàn)[15]中的算例，計(jì)算區(qū)域區(qū)取為[0,325]×[-44.5,44.5]，一Mach數(shù)為1.22的激波通過(guò)一個(gè)氣泡，氣泡中心位于(175,0)，半徑25，激波位于x=225，向左傳播，網(wǎng)格為320×80.無(wú)量綱化初始條件:氣泡內(nèi)ρ=0.138，u=0，v=0，p=1;激波前ρ=1，u=0，v=0，p=1;激波后ρ=1.376 4，u=-0.394，v=0，p=1.569 8.

圖1 流場(chǎng)不同時(shí)刻密度

圖2 流場(chǎng)不同時(shí)刻壓力

由圖1可以看出，平面沖擊波與氣泡相互作用后，氣泡向左運(yùn)動(dòng)，當(dāng)沖擊波掃過(guò)氣泡后，氣泡的迎波面(右側(cè))向內(nèi)凹陷，而背波面(左側(cè))形態(tài)基本保持不變.隨著時(shí)間推移，迎波面繼續(xù)向內(nèi)凹陷，迎波面形成渦環(huán)狀結(jié)構(gòu)，這是由于沖擊波誘導(dǎo)氣泡發(fā)生射流而形成的漩渦.由圖2可以看出，沖擊波傳播至氣泡處時(shí)，氣泡內(nèi)部壓力升高，并在氣泡上下靠近壁面處形成局部高壓，在t=100時(shí)，沖擊波在傳播方向的最前方形成高壓區(qū)，隨著時(shí)間推移，沖擊波誘導(dǎo)氣泡發(fā)生射流并在射流后方形成高壓區(qū).

由一維和二維算例可以看出本文的數(shù)值方法計(jì)算沖擊波壓力值較精確，能夠?qū)崿F(xiàn)對(duì)流場(chǎng)間斷的高分辨率捕捉，對(duì)激波具有較高的分辨率，通過(guò)在間斷附近構(gòu)造Riemann問(wèn)題的方式，能有效的抑制界面附近非物理震蕩，并能消除間斷附近的數(shù)值耗散.

3 近壁面水下爆炸沖擊波載荷數(shù)值模擬

初始時(shí)刻，爆炸流場(chǎng)內(nèi)密度為1 630 kg/m3，壓力為7.8×109Pa，外部流場(chǎng)水的密度為1 000 kg/m3，壓力為1×105Pa，爆炸初始流場(chǎng)半徑為1 m，上邊界為固壁邊界條件.

初始時(shí)刻，沖擊波和材料界面都迅速向外膨脹，并向周?chē)鲌?chǎng)中釋放一道球面沖擊波，稱(chēng)為主波，同時(shí)形成一道朝向氣泡中心匯聚的內(nèi)聚沖擊波，如圖3(a).當(dāng)沖擊波傳播至壁面后，被剛性壁面完全反射，并仍以沖擊波的形式朝向氣泡傳播.當(dāng)沖擊波傳播至氣泡表面后，由于氣泡外部介質(zhì)的阻抗高于氣泡內(nèi)部，因此一部分沖擊波將以稀疏波的形式反射回流場(chǎng)中，而另一部分將以沖擊波的形式透射入氣泡內(nèi)部，如圖3(b)、3(c).由于主波相對(duì)壁面來(lái)說(shuō)是以大角度入射的斜沖擊波，因此開(kāi)始出現(xiàn)馬赫反射現(xiàn)象.主波脫離壁面表面，形成馬赫桿結(jié)構(gòu)，與壁面反射波一起交匯于“三點(diǎn)波”處，并在三波結(jié)構(gòu)后出現(xiàn)一道滑移線(xiàn).主波、壁面反射波以及馬赫桿均為沖擊波，波前波后的密度、壓力和Ma均間斷，而滑移線(xiàn)內(nèi)側(cè)的壓力連續(xù)，密度和Ma間斷，如圖3(d)、(e).當(dāng)稀疏波傳播至壁面，即從低阻抗介質(zhì)入射至高阻抗介質(zhì)，又將以稀疏波的形式被反射，傳播方向相反的兩道稀疏波相互疊加，最終使得該區(qū)域內(nèi)的壓力低于臨界壓力以下，并形成空化，如圖3(f)所示.

圖3 近壁面水下爆炸流場(chǎng)物理特性

4 船體板架結(jié)構(gòu)響應(yīng)數(shù)值模擬

4.1 結(jié)構(gòu)模型及工況設(shè)置

建立船體板架結(jié)構(gòu)有限元模型，模型尺寸為6 m×4 m，厚度20 mm，綜合考慮計(jì)算精度和計(jì)算耗時(shí)，選取網(wǎng)格數(shù)量120×80，編制水下爆炸載荷計(jì)算程序并結(jié)合非線(xiàn)性有限元軟件ABAQUS，對(duì)船體板架結(jié)構(gòu)響應(yīng)進(jìn)行數(shù)值模擬.上邊界為船體板結(jié)構(gòu)，板四周剛性固定.應(yīng)用RKDG方法計(jì)算整個(gè)流場(chǎng)載荷，并與ABAQUS軟件結(jié)合，計(jì)算船體板架響應(yīng)特性.設(shè)板架中心的坐標(biāo)為O(0 m,0 m)，選取4個(gè)考核點(diǎn)，坐標(biāo)依次為A(2.0 m,0 m)、B(1.5 m,0 m)、C(1.0 m,0 m)、D(0.5 m,0 m).

4.2 船體板架結(jié)構(gòu)毀傷變形分析

當(dāng)沖擊波傳播作用到板架結(jié)構(gòu)時(shí)，板架承受壓力從0.1 MPa增大GPa至級(jí)，達(dá)到材料屈服極限，結(jié)構(gòu)出現(xiàn)塑性變形，計(jì)算得到的船體板架結(jié)構(gòu)變形如圖4所示.從圖4可以看出，沖擊波作用范圍內(nèi)船體板產(chǎn)生塑性變形，且沖擊波作用的圓形區(qū)域內(nèi)部結(jié)構(gòu)應(yīng)變較小，結(jié)構(gòu)產(chǎn)生圓盤(pán)狀拱起變形，與板格形成塑性鉸.

圖4 船體板架結(jié)構(gòu)塑性變形

4.3 船體板架結(jié)構(gòu)響應(yīng)分析

圖5給出了1.2 ms時(shí)船體板架結(jié)構(gòu)響應(yīng)云圖.圖6給出了A、B、C、D這4個(gè)考核點(diǎn)的板架結(jié)構(gòu)響應(yīng)計(jì)算結(jié)果.從圖6可以看出，沖擊波傳播過(guò)程中，位移及加速度響應(yīng)大小沿板格中心向外逐漸減小，離爆心最近的D點(diǎn)位移、加速度響應(yīng)最大，而A點(diǎn)最小，即板架結(jié)構(gòu)響應(yīng)與距離爆心距離呈反比.

5 結(jié) 論

1)能夠?qū)崿F(xiàn)對(duì)流場(chǎng)強(qiáng)間斷的高分辨率捕捉，對(duì)激波具有較高的分辨率，通過(guò)在間斷附近構(gòu)造Riemann問(wèn)題的方式，能有效的抑制界面附近非物理震蕩，并能消除間斷附近的數(shù)值耗散.

2)RKDG-GFM-LSM方法能很好的模擬近壁面水下爆炸沖擊波產(chǎn)生、傳播、反射和爆炸產(chǎn)物的膨脹等現(xiàn)象，并觀察到?jīng)_擊波傳播中的馬赫桿以及氣泡空化現(xiàn)象.船體板架結(jié)構(gòu)水下爆炸模擬中發(fā)現(xiàn)，結(jié)構(gòu)產(chǎn)生圓盤(pán)狀拱起變形，與板格形成塑性鉸，板架結(jié)構(gòu)響應(yīng)與距離爆心距離成反比.

圖5 船體板架結(jié)構(gòu)響應(yīng)云圖

圖6 考核點(diǎn)位移、加速度響應(yīng)時(shí)歷曲線(xiàn)

[1] SHIN Y S, SCHNEIDER N A. Ship shock trial simulation of USS Winston S. Churchill (DDG 81): modeling and simulation strategy and surrounding fluid volume effects[C]//Proceedings of the 74thShock and Vibration Symposium. San Diego, CA, Oct. : [s.n.], 2003: 27-31.

[2] 崔杰, 張阿漫, 郭君, 等. 艙段結(jié)構(gòu)在氣泡射流作用下的毀傷效果[J]. 爆炸與沖擊, 2012, 32(4): 355-361.

[3] 崔杰. 近場(chǎng)水下爆炸氣泡載荷及對(duì)結(jié)構(gòu)毀傷試驗(yàn)研究[D]. 哈爾濱: 哈爾濱工程大學(xué), 2013.

[4] SUSANTO A, IVAN L, De STERCK H, et al. High-order central ENO finite-volume scheme for ideal MHD[J]. Journal of Computational Physics, 2013, 250: 141-164.

[5] JEONG W, SEONG J. Comparison of effects on technical variances of computational fluid dynamics (CFD) software based on finite element and finite volume methods[J]. International Journal of Mechanical Sciences, 2014, 78: 19-26.

[6] QIU Jianxian, LIU Tiegang, KHOO B C. Simulations of compressible two-medium flow by runge-kutta discontinuous galerkinmethods with the ghost fluid method[J]. Communications in Computational Physics, 2008, 3(3): 479-504.[7] 張學(xué)瑩，趙寧，王春武. 可壓縮多介質(zhì)流體數(shù)值模擬中的Level-Set間斷跟蹤方法[J]. 計(jì)算物理,2006, 23(5): 518-524.[8] JOHNSEN E. Numerical simulations of non-spherical bubble collapse: with applications to shockwave lithotripsy[D]. Pasadena: California Institute of Technology, 2008.

[9] 任少飛. 基于間斷迦遼金法的艦船水下及空中爆炸模擬[D]. 哈爾濱：哈爾濱工程大學(xué), 2012.

[10]ZHU Jun, LIU Tiegang, QIU Jianxian, et al. RKDG methods with WENO limiters for unsteady cavitating flow[J]. Computers & Fluids, 2012, 57: 52-65.

[11]BILLET G, RYAN J, BORREL M. A Runge Kutta Discontinuous Galerkin approach to solve reactive flows on conforming hybrid grids: the parabolic and source operators[J]. Computers & Fluids, 2014, 95: 98-115.

[12]OSHER S, FEDKIW R P. Level set methods: an overview and some recent results[J]. Journal of Computational Physics, 2001, 169(2): 463-502.

[13]陳榮三. 大密度比和大壓力比可壓縮流的數(shù)值計(jì)算[J]. 應(yīng)用數(shù)學(xué)和力學(xué), 2008, 29(5): 609-617.

[14]BO W, GROVE J W. A volume of fluid method based ghost fluid method for compressible multi-fluid flows[J]. Computers & Fluids, 2014, 90: 113-122.

[15]HAAS J F, STURTEVANT B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J]. Journal of Fluid Mechanics, 1987, 181: 41-76.

(編輯 張 紅)

Numerical simulation of the response for Hull plates subjected to underwater explosion based on RKDG method

YU Fulin, GUO Jun, YAO Xiongliang, REN Shaofei

(College of Shipbuilding Engineering, Harbin Engineering University, 150001 Harbin，China)

In order to solve the underwater explosion flow field with large discontinuities, Level Set method was applied to track the interface position of the multi-medium flow, Ghost Fluid method was used to calculate the physical parameter of both sides of the interface, time and space were discretized by Runge-Kutta Discontinuous Galerkin Method, Euler equations of the flow field were solved. One-dimensional and two-dimensional assessments were conducted by RKDG approach. The results reflect the phenomena of underwater explosion shock wave generation, propagation, reflection and explosion products expansion. Finally, the shock responses and damage characteristics of hull plates under shock load were simulated with the nonlinear FEM softeware ABAQUS. The RKDG method can be applied to simulate the hull plates response with high accuracy. The response of hull plates is inversely proportional to the blast center distance.

underwater explosion; RKDG; shock load calculation; hull plates response; FEM; numerical simulation

10.11918/j.issn.0367-6234.2016.04.023

2014-11-26.

國(guó)家自然科學(xué)基金(51109042)；中國(guó)博士后科學(xué)基金(2012M520707)；黑龍江省自然科學(xué)基金 (E201124).

于福臨(1988—)，男，博士研究生；

姚熊亮(1963—)，男，教授，博士生導(dǎo)師.

姚熊亮，yaoxiongliang@163.com.

U663.2

A

0367-6234(2016)04-0139-05