基于RKDG方法的船體板架爆炸沖擊響應數值模擬

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

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

?

基于RKDG方法的船體板架爆炸沖擊響應數值模擬

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

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

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

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

1基本方程

1.1流體力學方程組

假設爆轟產物及周圍空氣為無粘可壓縮流，則流體力學方程可用下式表示[4]：

Ut+·F(U)=0

(1)

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

(2)

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

1.2Level Set方程

Level Set方程如下：

(3)

(4)

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

(5)

2RKDG方法

2.1空間離散

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

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

(6)

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

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

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

(7)

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

(8)

其中：

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

(9)

由于邊界內側流體的流通量不連續，可將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為勒讓德函數最高階數，ML求解如下：

(15)

ML=diag[m00m11m22]=

Δxdiag[11/34/45]

(16)

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

(17)

2.2時間離散

采用TVD Runge-Kutta方法對流體力學方程組進行時間離散。設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空中爆炸流場數值模擬

3.1自由場空中爆炸流場載荷數值模擬

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

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

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

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

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

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

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

表2　計算工況

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

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

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

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

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

4空中爆炸船體板架結構響應特性

4.1計算結構模型

圖5為船體板架結構有限元計算模型，板架四周剛性固定。

圖5　船體板架結構有限元模型示意圖Fig.5 The FEM of hull grillage

4.2計算工況和考核點

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

設板架中心的坐標為O(0 m,0 m)，選取四個考核點，坐標依次為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，計算得到沖擊波載荷作用大小和范圍，結合LS-DYNA軟件，計算船體板架結構變形，結果如圖7所示。

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

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

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

圖8　陳長海、Langdon G S等給出的試驗結果[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船體板架沖擊響應

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

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

從圖10及表4可以看出，沖擊波傳播過程中，位移、速度響應隨時間逐漸增大，位移及速度響應大小沿板格中心向外逐漸減小，離爆心最近的D點位移、速度響應最大，而A點最小，即板架結構響應與距離爆心距離呈反比。

圖9　船體加筋板架結構應力響應云圖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　位移響應時歷曲線Fig.10 The displacement response time history curve of the inspection point

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

5結論

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

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

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

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

參 考 文 獻

[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年生

摘要：基于歐拉流體力學方程組，對方程組分別進行間斷迦遼金法空間離散和龍格庫塔法時間離散，采用RKDG方法建立數值模型，模擬自由場空中爆炸沖擊波載荷特性，并與經驗公式對比，對RKDG計算程序進行驗證。針對強間斷性的空中爆炸沖擊波與周圍低壓空氣耦合特性問題，給出了流場密度和壓力分布規律。最后編制計算程序，結合大型非線性有限元軟件LS-DYNA，模擬板架結構在空中爆炸載荷作用下的變形和響應特征，計算發現：板架沿著T型材出現塑性變形，板格呈現局部隆起的碟形變形區，響應與距離爆心距離成反比。

關鍵詞：間斷伽遼金;空中爆炸;耦合特性;數值模擬

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

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2015.13.011

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

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

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