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

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

?

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

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

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

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

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

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

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

1 數(shù)值模擬方法

1.1 控制方程

對于無黏可壓縮流場，歐拉方程組可表示[6]為

(1)

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

(2)

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

1.2 狀態(tài)方程

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

(3)

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

(4)

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

1.3 空間與時間離散

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

(5)

(6)

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

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

(7)

式中U(0)為tn時刻流場物理特性.

1.4 界面追蹤

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

(8)

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

(9)

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

1.5 虛擬流動方法

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

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

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

(10)

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

(11)

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

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

(12)

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

(13)

2 算例驗(yàn)證

2.1 一維算例

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

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

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

2.2 二維激波管算例

選取文獻(xiàn)[15]中的算例，計(jì)算區(qū)域區(qū)取為[0,325]×[-44.5,44.5]，一Mach數(shù)為1.22的激波通過一個氣泡，氣泡中心位于(175,0)，半徑25，激波位于x=225，向左傳播，網(wǎng)格為320×80.無量綱化初始條件:氣泡內(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 流場不同時刻密度

圖2 流場不同時刻壓力

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

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

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

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

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

圖3 近壁面水下爆炸流場物理特性

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

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

建立船體板架結(jié)構(gòu)有限元模型，模型尺寸為6 m×4 m，厚度20 mm，綜合考慮計(jì)算精度和計(jì)算耗時，選取網(wǎng)格數(shù)量120×80，編制水下爆炸載荷計(jì)算程序并結(jié)合非線性有限元軟件ABAQUS，對船體板架結(jié)構(gòu)響應(yīng)進(jìn)行數(shù)值模擬.上邊界為船體板結(jié)構(gòu)，板四周剛性固定.應(yīng)用RKDG方法計(jì)算整個流場載荷，并與ABAQUS軟件結(jié)合，計(jì)算船體板架響應(yīng)特性.設(shè)板架中心的坐標(biāo)為O(0 m,0 m)，選取4個考核點(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)時，板架承受壓力從0.1 MPa增大GPa至級，達(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)生圓盤狀拱起變形，與板格形成塑性鉸.

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

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

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

5 結(jié) 論

1)能夠?qū)崿F(xiàn)對流場強(qiáng)間斷的高分辨率捕捉，對激波具有較高的分辨率，通過在間斷附近構(gòu)造Riemann問題的方式，能有效的抑制界面附近非物理震蕩，并能消除間斷附近的數(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)生圓盤狀拱起變形，與板格形成塑性鉸，板架結(jié)構(gòu)響應(yīng)與距離爆心距離成反比.

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

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

[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] 崔杰. 近場水下爆炸氣泡載荷及對結(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.

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

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

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

姚熊亮，yaoxiongliang@163.com.

U663.2

A

0367-6234(2016)04-0139-05