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

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

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

?

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

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

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

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

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

1基本方程

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

假設(shè)爆轟產(chǎn)物及周圍空氣為無粘可壓縮流，則流體力學(xué)方程可用下式表示[4]：

Ut+·F(U)=0

(1)

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

(2)

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

1.2Level Set方程

Level Set方程如下：

(3)

(4)

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

(5)

2RKDG方法

2.1空間離散

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

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

(6)

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

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

在式(6)兩側(cè)同時乘以試探函數(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))進行如下處理[9]：

其中：

(10)

因此式(10)可改寫為

(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)

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

(17)

2.2時間離散

采用TVD Runge-Kutta方法對流體力學(xué)方程組進行時間離散。設(shè)tn時刻流場物理特性為U(0)，則tn+1時刻流場物理特性可由式(18)進行求解[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為勒讓德多項式。

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

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

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

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

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

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

3.2具有強間斷性的空中爆炸流場物理特性

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

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

表2　計算工況

圖3　不同時刻流場密度分布Fig.3 The density distribution of flow field simulation

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

圖4　不同時刻流場壓力分布Fig.4 The pressure distribution of flow field simulation

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

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

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

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

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

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

4.2計算工況和考核點

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

設(shè)板架中心的坐標(biāo)為O(0 m,0 m)，選取四個考核點，坐標(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程序計算時間t=3.0 ms，計算得到?jīng)_擊波載荷作用大小和范圍，結(jié)合LS-DYNA軟件，計算船體板架結(jié)構(gòu)變形，結(jié)果如圖7所示。

圖6　流場初始條件Fig.6 The initial condition ofair explosion

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

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

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

時刻/ms1.52.02.5位移/m0.0870.1370.184

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

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

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

從圖10及表4可以看出，沖擊波傳播過程中，位移、速度響應(yīng)隨時間逐漸增大，位移及速度響應(yīng)大小沿板格中心向外逐漸減小，離爆心最近的D點位移、速度響應(yīng)最大，而A點最小，即板架結(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

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

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

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

5結(jié)論

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

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

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

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

參 考 文 獻

[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] 陳長海,朱錫,侯海量,等.近距空爆載荷作用下固支方板的變形及破壞模式[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é)方程組，對方程組分別進行間斷迦遼金法空間離散和龍格庫塔法時間離散，采用RKDG方法建立數(shù)值模型，模擬自由場空中爆炸沖擊波載荷特性，并與經(jīng)驗公式對比，對RKDG計算程序進行驗證。針對強間斷性的空中爆炸沖擊波與周圍低壓空氣耦合特性問題，給出了流場密度和壓力分布規(guī)律。最后編制計算程序，結(jié)合大型非線性有限元軟件LS-DYNA，模擬板架結(jié)構(gòu)在空中爆炸載荷作用下的變形和響應(yīng)特征，計算發(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

中圖分類號:TH212；TH213.3

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

DOI:10.13465/j.cnki.jvs.2015.13.011

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

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

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