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

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

?

基于RKDG方法的船體板架爆炸沖擊響應(yīng)數(shù)值模擬

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

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

RKDG方法(簡(jiǎn)稱(chēng)間斷伽遼金方法)由Reed等[1]提出，在很長(zhǎng)一段時(shí)間內(nèi)僅用于線(xiàn)性方程的求解。Cockburn等[2]結(jié)合顯式Runge-Kutta方法，利用間斷伽遼金法求解了雙曲型守恒律方程，間斷伽遼金法作為一種新型的計(jì)算流體方法開(kāi)始迅速發(fā)展。間斷伽遼金方法是介于有限元法和有限體積法之間的一種CFD算法[3]，具備兩者的優(yōu)點(diǎn)。一方面，間斷伽遼金法本質(zhì)上是一種高精度的微分方程空間離散方法，在離散過(guò)程中使用弱形式積分，具有有限維解，在單元內(nèi)使用多項(xiàng)式逼近物理量的真實(shí)值，通過(guò)簡(jiǎn)單的增加多項(xiàng)式次數(shù)，即可直接提高單元階數(shù)，構(gòu)造高階格式，也與有限元方法相似。此外，相比其它可用于流場(chǎng)間斷捕捉的空間離散方法，如WENO[4]格式，ENO[5]格式等，間斷伽遼金方法不需要大規(guī)模的模板，緊湊型極好，求解過(guò)程中僅需要當(dāng)前單元及其直接相鄰的單元即可。另一方面，間斷伽遼金方法引入數(shù)值通量的概念，在單元與單元間計(jì)入流場(chǎng)間斷，使其可用于沖擊波等強(qiáng)間斷現(xiàn)象的捕捉，同時(shí)還具備多種非線(xiàn)性的限制器，在保證其精度的前提下，抑制非物理振蕩，這與有限體積法相似。間斷伽遼金法本身就具有高分辨率捕捉的特性，能有效處理水下爆炸中出現(xiàn)的壓力、密度、速度等物理量的間斷。

本文編制RKDG方法Fortran計(jì)算程序，對(duì)空中爆炸流場(chǎng)物理特性進(jìn)行數(shù)值模擬，求解沖擊波載荷特性，并結(jié)合大型非線(xiàn)性有限元軟件LS-DYNA，模擬船體加筋板架結(jié)構(gòu)的爆炸沖擊響應(yīng)特征。

1基本方程

1.1流體力學(xué)方程組

假設(shè)爆轟產(chǎn)物及周?chē)諝鉃闊o(wú)粘可壓縮流，則流體力學(xué)方程可用下式表示[4]：

Ut+·F(U)=0

(1)

二維無(wú)粘可壓縮流的歐拉流體力學(xué)方程組如下式所示：

(2)

其中：ρ為流體密度，p為壓力，u,v,w為速度，E為單位體積總能量。

1.2Level Set方程

Level Set方程如下：

(3)

(4)

為方便求解將符號(hào)函數(shù)S(Φ0)光滑化為[6]：

(5)

2RKDG方法

2.1空間離散

由于爆轟流體運(yùn)動(dòng)速度相當(dāng)快，所以可被視為無(wú)粘可壓縮流。一維無(wú)粘可壓縮流的DG方法空間離散如圖1所示[7]，歐拉方程如下式所示：

U(x,t)t+F(U(x,t))x=0

(6)

將計(jì)算區(qū)域Ω劃分為N個(gè)單元，每個(gè)單元長(zhǎng)度為dxi=dxi+1/2-dxi-1/2。

圖1　一維DG方法空間離散Fig.1 Discontinuous Galerkin method in 1D spatial discretization

在式(6)兩側(cè)同時(shí)乘以試探函數(shù)w(x)，積分得[8]

(7)

式(7)分部積分得[8]：

(8)

其中：

F(U(xi-1/2,t)t)w(xi-1/2)

(9)

由于邊界內(nèi)側(cè)流體的流通量不連續(xù)，可將F(U(xi±1/2,t))進(jìn)行如下處理[9]：

其中：

(10)

因此式(10)可改寫(xiě)為

(11)

U(t)+Ux(t)P1ξ+Uxx(t)(ξ2-1/3)

(13)

將式(13)代入式(12)中得

(14)

其中:k=0,1,…,p，p為勒讓德函數(shù)最高階數(shù)，ML求解如下：

(15)

ML=diag[m00m11m22]=

Δxdiag[11/34/45]

(16)

對(duì)于多維流體力學(xué)方程組采用同樣方法，可以得到下式[8]。

(17)

2.2時(shí)間離散

采用TVD Runge-Kutta方法對(duì)流體力學(xué)方程組進(jìn)行時(shí)間離散。設(shè)tn時(shí)刻流場(chǎng)物理特性為U(0)，則tn+1時(shí)刻流場(chǎng)物理特性可由式(18)進(jìn)行求解[11]：

U(1)=U(0)+ΔtR(U(0))

(18)

其中:R(Uk)可由下式求解[12]:

[Pk(1)F(Uk(1,t))i-

Pk(-1)F(Uk(-1,t))i]}

(19)

Pk為勒讓德多項(xiàng)式。

3空中爆炸流場(chǎng)數(shù)值模擬

3.1自由場(chǎng)空中爆炸流場(chǎng)載荷數(shù)值模擬

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

圖2　沖擊波超壓計(jì)算值與經(jīng)驗(yàn)值對(duì)比圖Fig.2 The comparison of the peak overpressure of blast wave

距離/m0.500.650.750.901經(jīng)驗(yàn)公式解/MPa36.8726.9222.3817.5315.18間斷迦遼金法解/MPa32.1226.4421.7918.8515.34誤差/%12.881.772.627.511.02

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

3.2具有強(qiáng)間斷性的空中爆炸流場(chǎng)物理特性

藥包為球形TNT裝藥，空氣與爆轟產(chǎn)物采用理想氣體狀態(tài)方程，r=1.4。上下邊界設(shè)置為固壁邊界，左右邊界條件設(shè)置為流體流出。由RKDG方法計(jì)算得到不同時(shí)刻整個(gè)區(qū)域流場(chǎng)的密度、壓力分布。

空中爆炸流場(chǎng)密度計(jì)算結(jié)果如圖3所示。

表2　計(jì)算工況

圖3　不同時(shí)刻流場(chǎng)密度分布Fig.3 The density distribution of flow field simulation

空中爆炸流場(chǎng)壓力計(jì)算結(jié)果如圖4所示。

圖4　不同時(shí)刻流場(chǎng)壓力分布Fig.4 The pressure distribution of flow field simulation

由圖3可以看出，流場(chǎng)密度隨著時(shí)間的推移迅速衰減，密度由圓心向四周逐漸減小，當(dāng)沖擊波到達(dá)剛性壁面時(shí)，壁面密度增大，且邊界中心密度大于兩側(cè)密度。

由圖4可以看出，與密度變化規(guī)律相似，沖擊波超壓隨著時(shí)間推移迅速衰減，沖擊波傳播至上下固壁邊界時(shí)，壓力迅速上升，在邊界中心形成局部高壓，甚至在t=4.0 ms時(shí)高于爆心壓力。

4空中爆炸船體板架結(jié)構(gòu)響應(yīng)特性

4.1計(jì)算結(jié)構(gòu)模型

圖5為船體板架結(jié)構(gòu)有限元計(jì)算模型，板架四周剛性固定。

圖5　船體板架結(jié)構(gòu)有限元模型示意圖Fig.5 The FEM of hull grillage

4.2計(jì)算工況和考核點(diǎn)

爆炸流場(chǎng)初始條件如圖6所示，上邊界為船體加筋板架結(jié)構(gòu)。應(yīng)用RKDG方法計(jì)算整個(gè)流場(chǎng)載荷，并與LS-DYNA軟件結(jié)合，計(jì)算船體板架響應(yīng)特性。

設(shè)板架中心的坐標(biāo)為O(0 m,0 m)，選取四個(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.3船體板架毀傷變形

RKDG程序計(jì)算時(shí)間t=3.0 ms，計(jì)算得到?jīng)_擊波載荷作用大小和范圍，結(jié)合LS-DYNA軟件，計(jì)算船體板架結(jié)構(gòu)變形，結(jié)果如圖7所示。

圖6　流場(chǎng)初始條件Fig.6 The initial condition ofair explosion

圖7　船體板架毀傷變形Fig.7 The damage deformation of hull grillage

沖擊波傳播首先達(dá)到板架中心，如強(qiáng)間斷性空中爆炸沖擊波載荷模擬中展示的規(guī)律，板架中心處的沖擊波壓力最大，沿中心向四周壓力不斷衰減，從圖8可以看出，板架中心變形最大，沿板架中心向四周變形逐漸減小，這與沖擊波超壓分布特性一致，板架結(jié)構(gòu)呈現(xiàn)局部隆起的碟形變形區(qū)，且中心處位移最大。這與陳長(zhǎng)海、Langdon等[14-15]試驗(yàn)變形結(jié)果是一致的。

圖8　陳長(zhǎng)海、Langdon G S等給出的試驗(yàn)結(jié)果[14-15]Fig.8 The test results provided by CHEN Changhai and Langdon G S

時(shí)刻/ms1.52.02.5位移/m0.0870.1370.184

4.4船體板架沖擊響應(yīng)

圖9為t=0.05 ms時(shí)船體板架結(jié)構(gòu)應(yīng)力響應(yīng)云圖，可以看出，沖擊波作用于板架后，T型材區(qū)域應(yīng)力響應(yīng)明顯大于板格應(yīng)力響應(yīng)，且在交叉梁處產(chǎn)生了應(yīng)力集中，應(yīng)力沿中心向四周逐漸減小。此時(shí)最大應(yīng)力值為5.1×108Pa，在空中爆炸載荷作用下板架沿著T型材出現(xiàn)塑性變形。這與London等[14]觀(guān)察到的試驗(yàn)結(jié)果是一致的。

圖10、圖11為不同考核點(diǎn)的位移和速度響應(yīng)曲線(xiàn)。

從圖10及表4可以看出，沖擊波傳播過(guò)程中，位移、速度響應(yīng)隨時(shí)間逐漸增大，位移及速度響應(yīng)大小沿板格中心向外逐漸減小，離爆心最近的D點(diǎn)位移、速度響應(yīng)最大，而A點(diǎn)最小，即板架結(jié)構(gòu)響應(yīng)與距離爆心距離呈反比。

圖9　船體加筋板架結(jié)構(gòu)應(yīng)力響應(yīng)云圖Fig.9 The stress nephogram of hull grillage under far field air explosion

考核點(diǎn)ABCD位移值/cm0.0280.0310.1490.163速度/(m·s-1)16.520.449.8107.6

圖10　位移響應(yīng)時(shí)歷曲線(xiàn)Fig.10 The displacement response time history curve of the inspection point

圖11　速度響應(yīng)時(shí)歷曲線(xiàn)Fig.11 The velocity response time history curve of the inspection point

5結(jié)論

本文將RKDG方法應(yīng)用于空中爆炸流場(chǎng)物理特性求解中，模擬了自由場(chǎng)空中爆炸載荷以及具有強(qiáng)間斷性的空中爆炸沖擊波流場(chǎng)，給出了流場(chǎng)載荷特性。結(jié)合有限元方法，對(duì)爆炸載荷作用下典型船體板架結(jié)構(gòu)的響應(yīng)就行了分析，結(jié)論如下：

(1) RKDG方法計(jì)算得到的空中爆炸沖擊波載荷計(jì)算結(jié)果與經(jīng)驗(yàn)公式吻合較好，不同爆距下的超壓峰值最小誤差為1.02%，最大誤差不超過(guò)13%，可以有效模擬空中爆炸沖擊波載荷特性。

(2) 模擬具有強(qiáng)間斷性的空中爆炸流場(chǎng)載荷特性時(shí)發(fā)現(xiàn)，流場(chǎng)密度隨著時(shí)間的推移迅速衰減，以爆心為圓心的圓形區(qū)域內(nèi)密度始終最大；而沖擊波超壓也隨著時(shí)間推移迅速衰減，沖擊波傳播至上下固壁邊界時(shí)，壓力迅速上升，在邊界中心形成局部高壓，甚至在t=4.0 ms時(shí)邊界中心沖擊波超壓高于爆心壓力。

(3) 分析船體板架結(jié)構(gòu)在空中爆炸載荷作用下的響應(yīng)特征時(shí)，發(fā)現(xiàn)板架結(jié)構(gòu)呈現(xiàn)局部隆起的碟形變形，且中心處位移最大，在空中爆炸載荷作用下板架沿著T型材出現(xiàn)塑性變形，與試驗(yàn)結(jié)果一致，并發(fā)現(xiàn)板架結(jié)構(gòu)響應(yīng)與距離爆心距離呈反比。

參 考 文 獻(xiàn)

[1]Reed W H, Hill T R. Triangular mesh methods for the neutron transport equation [J]. Los Alamos Report LA-UR-73-479, 1973.

[2]Cockburn B, Shu C W. TVB Runge-Kutta local projection discontinuous galerkin finite element method for conservation laws II: General framework[J]. Mathematics of Computation, 1989, 52(186): 411-435.

[3]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.

[4]Zhu J, Liu T, Qiu J, et al. RKDG methods with WENO limiters for unsteady cavitating flow[J]. Computers & Fluids, 2012, 57: 52-65.

[5]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.

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

[7]Cheng Y, Gamba I M, Majorana A, et al. Discontinuous galerkin solver for Boltzmann-poisson transients[J]. Journal of Computational Electronics, 2008, 7(3): 119-123.

[8]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.

[9]Wang Z J. High-order computational fluid dynamics tools for aircraft design[J]. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 2014, 372(2022): 20130318..

[10] Du J, Shu C W, Zhang M. A simple weighted essentially non-oscillatory limiter for the correction procedure via reconstruction (CPR) framework[J]. Applied Numerical Mathematics, 2014,279:1-26.

[11] 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.

[12] Pesch L, van der Vegt J J W. A discontinuous Galerkin finite element discretization of the Euler equations for compressible and incompressible fluids [J].Journal of Computational Physics, 2008, 227: 5426-5446.

[13] Henrych J, Major R. The dynamics of explosion and its use[M]. Amsterdam: Elsevier, 1979.

[14] 陳長(zhǎng)海,朱錫,侯海量,等.近距空爆載荷作用下固支方板的變形及破壞模式[J]. 爆炸與沖擊, 2008, 32(4): 368-375.

CHEN Chang-hai, ZHU Xi, HOU Hai-liang, et al. Deformation and failure modes of clamped squared plates under close-range air blast loads[J]. Explosion and Shock Waves, 2008, 32(4): 368-375.

第一作者 于福臨 男，博士生，1988年生

摘要：基于歐拉流體力學(xué)方程組，對(duì)方程組分別進(jìn)行間斷迦遼金法空間離散和龍格庫(kù)塔法時(shí)間離散，采用RKDG方法建立數(shù)值模型，模擬自由場(chǎng)空中爆炸沖擊波載荷特性，并與經(jīng)驗(yàn)公式對(duì)比，對(duì)RKDG計(jì)算程序進(jìn)行驗(yàn)證。針對(duì)強(qiáng)間斷性的空中爆炸沖擊波與周?chē)蛪嚎諝怦詈咸匦詥?wèn)題，給出了流場(chǎng)密度和壓力分布規(guī)律。最后編制計(jì)算程序，結(jié)合大型非線(xiàn)性有限元軟件LS-DYNA，模擬板架結(jié)構(gòu)在空中爆炸載荷作用下的變形和響應(yīng)特征，計(jì)算發(fā)現(xiàn)：板架沿著T型材出現(xiàn)塑性變形，板格呈現(xiàn)局部隆起的碟形變形區(qū)，響應(yīng)與距離爆心距離成反比。

關(guān)鍵詞：間斷伽遼金;空中爆炸;耦合特性;數(shù)值模擬

Numerical simulation of hull grillage responses under blast loading based on RKDG method

YUFu-lin,GUOJun,YAOXiong-liang,RENShao-fei(College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China)

Abstract:In order to simulate the characteristics of far field air explosion shock load, the discretized Euler equations were solved with RKDG method. The numerical model was established and the far field air explosion shock wave load characteristics were simulated. The results were compared with those of the empirical formula. The density and pressure distribution of high pressure shock wave with strong discontinuity were simulated with a Fortran program. Using LS-DYNA, the impulse responses of hull plates under shock load were simulated. The results indicated that the plastic deformation shape profiles appear along the stiffener direction, and the plate grillages reveal local dishing deformation, and the responses are inversely proportional to the distance from the explosion center.

Key words:discontinuous Galerkin method; air explosion; coupling characteristics; numerical simulation

中圖分類(lèi)號(hào):TH212；TH213.3

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2015.13.011

通信作者郭君 男，博士，副教授，1981年生

收稿日期：2014-09-23修改稿收到日期：2015-01-20

基金項(xiàng)目：國(guó)家自然科學(xué)基金資助項(xiàng)目(51109042)；中國(guó)博士后科學(xué)基金資助項(xiàng)目(2012M520707)；黑龍江省自然科學(xué)基金資助項(xiàng)目(E201124)；黑龍江省博士后資助項(xiàng)目(LBH-Z12064)