基于多介質(zhì)Riemann問題的流體- 固體耦合數(shù)值方法及其在爆炸與沖擊問題中的應(yīng)用

姚成寶， 付梅艷， 韓峰， 閆凱， 雷雨

(西北核技術(shù)研究所， 陜西 西安 710024)

0 引言

地下強(qiáng)爆炸和高速沖擊等物理問題屬于可壓縮多介質(zhì)流動的范疇，通常具有多介質(zhì)、大變形以及高度非線性等特點，物理現(xiàn)象極其復(fù)雜。隨著研究需求的不斷提高，多介質(zhì)流動問題的高置信度數(shù)值模擬逐漸成為必不可少的研究手段，是計算數(shù)學(xué)和計算力學(xué)研究領(lǐng)域中極具挑戰(zhàn)性的課題，具有重要的研究價值。

多介質(zhì)流動數(shù)值方法按其采用的坐標(biāo)系可以分為Lagrange方法、Euler方法、任意Lagrange-Euler方法和Euler-Lagrange耦合方法等。由于Euler方法的坐標(biāo)系固定在空間中，計算網(wǎng)格不隨介質(zhì)的運動而移動，不發(fā)生網(wǎng)格畸變，因此適合處理大變形問題。Euler坐標(biāo)系下的多介質(zhì)流動問題數(shù)值模擬包括兩個關(guān)鍵過程：界面追蹤和介質(zhì)間相互作用。如何描述界面的位置以及混合網(wǎng)格中各介質(zhì)之間的相互作用，是Euler方法求解多介質(zhì)問題的難點和關(guān)鍵。Udaykumar等[1-3]建立了基于Euler坐標(biāo)系下的 Cartesian 網(wǎng)格、具有銳利界面的二維、三維多介質(zhì)彈塑性固體的數(shù)值方法，可用于模擬彈塑性固體在受到高速沖擊和爆炸載荷時的動力學(xué)響應(yīng)。劉凱欣等[4-6]采用時空有限元(CE/SE)方法來求解Euler坐標(biāo)系下的彈塑性固體控制方程組，利用Hybrid Level set方法來捕捉介質(zhì)界面，并采用虛擬流體方法來處理邊界條件和介質(zhì)間的相互作用,對金屬材料受到?jīng)_擊和爆炸等問題進(jìn)行了數(shù)值模擬。丁建許等[7]、王仲琦等[8]采用有限差分方法和流體體積法、水平集方法等界面處理技術(shù)對Euler坐標(biāo)系下彈塑性固體的侵徹和聚能射流等問題進(jìn)行了數(shù)值模擬。

無論采用何種方法來處理介質(zhì)之間的相互作用，如何確定物質(zhì)界面上的物理量參數(shù)是一個關(guān)鍵環(huán)節(jié)。由于物質(zhì)界面實際上是一個接觸間斷，求解物質(zhì)界面上的多介質(zhì) Riemann 問題是處理該類問題的最直接和最有效方法。Riemann問題的精確解能精確捕捉激波和接觸間斷的位置，基于Riemann問題精確解的數(shù)值方法也具有明顯的優(yōu)勢。目前國內(nèi)外已有很多學(xué)者將 Riemann問題的精確解應(yīng)用于多介質(zhì)流動問題的求解，用來提高計算精度。其中，Banks[9]對Jones-Wilkins-Lee(JWL)狀態(tài)方程的 Riemann 問題進(jìn)行了理論分析。Kamm[10]針對滿足凸性的通用狀態(tài)方程，利用二分法迭代求解接觸間斷上的壓力p滿足的代數(shù)方程，得到了Riemann問題的精確解。Despres[11]對非守恒形式亞彈性模型的激波理論進(jìn)行分析，提出了彈性固體的密度ρ與彈性剪切模量μe的乘積βe在可逆的彈性變形中保持為常數(shù)(βe=ρμe)的假設(shè)，獲得了偏應(yīng)力跨過激波時滿足的Rankine-Hugoniot關(guān)系。Gavrilyuk等[12]和Kamm[13]在Despres[11]的基礎(chǔ)上，假定在彈性變形中滿足dβe/dt=0，分析了考慮剪切波時線彈性固體的Riemann問題精確解。Lin等[14]對Riemann不變量進(jìn)行了線性化處理，并對理想彈塑性固體提出了一種基于迭代方法的近似 Riemann算子。Trangenstein等[15]對一維單軸薄壁圓管的理想彈塑性Riemann 問題進(jìn)行了分析。Abuzyarov等[16]和Bazhenov等[17]在忽略內(nèi)能變化時，對等壓形式的理想彈塑性固體的波系結(jié)構(gòu)(激波和稀疏波)進(jìn)行了理論分析。Menshov等[18]基于Gavrilyuk等[12]的假設(shè)，對一維應(yīng)變情形時理想彈塑性固體的Riemann問題進(jìn)行了理論分析。Tang等[19]對具有Murnagham狀態(tài)方程和理想彈塑性本構(gòu)方程的Riemann問題進(jìn)行了分析，并利用Lagrange方法對氣體- 固體、氣體- 液體以及固體- 固體(以下簡稱固固)等多介質(zhì)問題進(jìn)行了一維數(shù)值模擬。Liu等[20]在Tang等[19]的基礎(chǔ)上利用有限體積方法并結(jié)合修正虛擬流體方法對固體在受到強(qiáng)沖擊時的流體- 固體(以下簡稱流固)耦合問題進(jìn)行了理論分析和數(shù)值模擬。Gao等[21-22]對一維理想彈塑性固體的Riemann問題精確解進(jìn)行了理論分析，其中固體在彈性階段和塑性階段的靜水壓力分別采用等熵狀態(tài)方程和剛性氣體狀態(tài)方程，并基于分裂波理論給出了固體處于彈塑性階段的Riemann問題精確解。Berjamin等[23]對一維應(yīng)變下的彈性本構(gòu)方程進(jìn)行了研究，并提出了當(dāng)彈性本構(gòu)模型出現(xiàn)非凸時的Riemann解子，其求解原理和流體狀態(tài)方程出現(xiàn)非凸時基本一致。Hattori[24]研究了考慮相變時的熱彈性本構(gòu)模型的 Riemann問題，并提出了相應(yīng)的Riemann解子。Feng等[25-26]對線性硬化彈塑性固體一維Riemann問題的特征結(jié)構(gòu)進(jìn)行了分析，并結(jié)合剛性氣體狀態(tài)方程，給出了一維流固Riemann問題以及固固Riemann問題的精確解。

本文針對具有高度非線性特征、通用形式的Mie-Grüneisen 狀態(tài)方程和流體彈塑性本構(gòu)方程，開發(fā)了一種健壯、高效、誤差可控的多介質(zhì)Riemann問題求解器，并設(shè)計了一種能夠針對強(qiáng)激波和強(qiáng)稀疏波等極端情形的非精確Newton方法來迭代求解Riemann問題所滿足的非線性代數(shù)- 積分組合方程，能有效提高物質(zhì)界面上各物理狀態(tài)量的計算精度。結(jié)合基于Euler坐標(biāo)系的非結(jié)構(gòu)網(wǎng)格、具有銳利界面的守恒型多介質(zhì)流動數(shù)值方法，建立了一套能夠模擬具有高密度比、高壓力比以及復(fù)雜非線性狀態(tài)方程的流固耦合、固固耦合等問題的多介質(zhì)流動計算體系，可用于模擬可壓縮流體和彈塑性固體在極端物理條件下的大變形動力學(xué)行為。最后開展了流固耦合、固固耦合Riemann問題、地下強(qiáng)爆炸、內(nèi)爆壓縮問題以及高速沖擊等大變形問題的數(shù)值模擬，并和理論、實測結(jié)果進(jìn)行比對，結(jié)果表明本文數(shù)值方法能夠有效處理多介質(zhì)流體和彈塑性固體的大變形問題。

1 物理模型

1.1 控制方程

考慮二維計算區(qū)域上的多介質(zhì)流固耦合問題，基于Euler坐標(biāo)系描述可壓縮流體和彈塑性固體在受到爆炸、沖擊等強(qiáng)載荷作用下的大變形動力學(xué)行為，其控制方程組如(1)式所示：

(1)

式中：W表示守恒量；F為通量；H為源項；ρ為密度；u為速度矢量；E為單位體積總能量；p為壓力；S為固體的偏應(yīng)力張量，對于流體S等于0，以下皆同。

1.2 狀態(tài)方程和本構(gòu)模型

除了描述物質(zhì)(流體或固體)的質(zhì)量、動量和能量守恒的控制方程外，物質(zhì)動力學(xué)行為的完整描述還需要提供描述自身熱力學(xué)狀態(tài)的本構(gòu)方程(狀態(tài)方程)。為后續(xù)進(jìn)行統(tǒng)一分析，將流體狀態(tài)方程和固體靜水壓力方程寫成(2)式的統(tǒng)一形式：

p=ωρe+h(ρ)，

(2)

式中：ω為Grüneisen系數(shù)；e為單位質(zhì)量比內(nèi)能；h(ρ)為參考狀態(tài)。本文涉及到的狀態(tài)方程包括：

梯恩梯(TNT)爆轟產(chǎn)物采用JWL狀態(tài)方程，如(3)式所示：

(3)

式中：ρ0為初始密度；e且滿足e=E/ρ-(u2+v2)/2，u,v為速度矢量u的x軸、y軸方向的分量；A1、B1、R1、R2和ω為JWL狀態(tài)方程的參數(shù)，如表1所示。

表1 JWL狀態(tài)方程參數(shù)

空氣采用理想氣體狀態(tài)方程，如(4)式所示：

p=(γ-1)ρe，

(4)

式中：γ為空氣的絕熱指數(shù)，γ=1.4.

剛性氣體狀態(tài)方程(5)式可用于描述金屬、水等物質(zhì)在受到爆炸或沖擊時的動態(tài)響應(yīng)問題，

p=(γ-1)ρe-γp∞，

(5)

式中：p∞表示參考壓力。其中：鋼的靜水壓力狀態(tài)方程參數(shù)為γ= 4.4，p∞= 6×108Pa；水的參數(shù)為γ=7.15，p∞=3.31×108Pa.

Murnagham狀態(tài)方程是一種線性化的Grüneisen狀態(tài)方程，可用于描述金屬等材料在受到?jīng)_擊時的動態(tài)響應(yīng)問題，

(6)

式中：K表示體積模量，K=2.225×1011Pa；p0表示初始壓力。其中，鋼的靜水壓力狀態(tài)方程參數(shù)為ρ0=7 800 kg/m3，p0=1.0×105Pa，γ=3.7.

與流體不同的是，固體材料除了能夠承受壓縮或膨脹帶來的體積變化外，還能夠承受一定程度的剪切變形(形狀改變)。流體彈塑性模型將固體的力學(xué)響應(yīng)過程分解為兩個部分：體積變形和剪切變形，其中體積變形仍然由狀態(tài)方程(2)式來確定，而剪切變形在受到強(qiáng)沖擊或者爆炸載荷作用時可能經(jīng)歷3個過程：彈性變形階段、塑性變形階段以及流體動力學(xué)階段，如(7)式所示。

(7)

2 流固耦合Riemann問題求解

2.1 流固耦合Riemann問題

由于多維Riemann問題的復(fù)雜性，相關(guān)的理論研究十分困難，常見的做法是將多維Riemann問題簡化為界面法向上的局部一維Riemann問題[27]，利用界面兩側(cè)介質(zhì)的壓力和法向速度的相容性條件來求解一維Riemann問題。

對于流體與彈塑性固體耦合的一維多介質(zhì)Riemann問題，當(dāng)不考慮黏性時，Riemann問題的解由初值條件以及相應(yīng)介質(zhì)的狀態(tài)方程或本構(gòu)方程決定。不失一般性，假設(shè)流體位于物質(zhì)界面的左側(cè)，彈塑性固體位于物質(zhì)界面的右側(cè)，此時一維流固耦合的多介質(zhì)Riemann問題定義為

Uτ+Ψξ(U)=0，

(8)

圖1 一維流固耦合Riemann問題的波系結(jié)構(gòu)Fig.1 Typical wave structures of one-dimensional fluid-solid Riemann problem

10 kV配電網(wǎng)不停電作業(yè)綜合絕緣抱桿也獲得了行業(yè)內(nèi)的高度認(rèn)可，于2017年7月獲實用新型專利(專利號ZL 201720051536.5)，12月獲國網(wǎng)四川省電力公司工人技術(shù)創(chuàng)新一等獎。

2.2 流體中的波系結(jié)構(gòu)

對于流體，當(dāng)p>pl時，非線性波的類型為激波，且激波前、后的物理量狀態(tài)滿足Rankine-Hugoniot關(guān)系：

當(dāng)p≤pl時，非線性波的類型為稀疏波，且稀疏波前、后的Riemann不變量滿足：

式中：c表示聲速，且

綜上所述，fl(p,Ul)=-(u-ul)滿足：

(9)

對于JWL、多項式等形式復(fù)雜的狀態(tài)方程，fl(p,Ul)具有高度非線性，解析求解(9)式非常困難。本文采用數(shù)值方法來近似求解fl(p,Ul)，并通過在數(shù)值求解過程中進(jìn)行誤差控制來保證近似解收斂到(9)式的精確解。

2.2.1 激波

當(dāng)p>pl時，該非線性波的類型為激波。對(2)式進(jìn)行變換，

el(p,ρ)=(p-hl(ρ))/ωlρ.

可以建立激波曲線上p和ρ之間滿足的關(guān)系(Hugoniot函數(shù))：

(10)

利用牛頓迭代法對(10)式進(jìn)行迭代求解：

(11)

2.2.2 稀疏波

當(dāng)p≤pl時，該非線性波的類型為稀疏波，p和ρ滿足等熵關(guān)系dρ/dp=1/c2. 結(jié)合(9)式的稀疏波關(guān)系，可得

(12)

利用自適應(yīng)步長的Runge-Kutta-Fehlberg方法[28-29]對(12)式進(jìn)行聯(lián)立求解，可完成稀疏波分支fl(p,Ul)的計算。

2.3 固體中的波系結(jié)構(gòu)

由于固體在彈性和塑性階段具有不同的狀態(tài)方程和本構(gòu)模型，當(dāng)固體經(jīng)歷從彈性到塑性的相變過程時，稀疏波或激波的前、后狀態(tài)跨過了固體的彈性、塑性屈服極限，稀疏波的等熵線或激波的Rankine-Hugoniot曲線的斜率將發(fā)生間斷，產(chǎn)生明顯不同于流體的波系結(jié)構(gòu)。下面依次對各變形階段的波系結(jié)構(gòu)進(jìn)行分析。

2.3.1 彈性變形階段的波系結(jié)構(gòu)

2.3.1.1 彈性激波

(13)

且靜水壓力p和偏應(yīng)力s分別滿足

(14)

與(10)式類似，建立彈性固體在跨過激波時靜水壓力p和密度ρ之間的關(guān)系為

結(jié)合(14)式，建立彈性激波曲線上σ和ρ之間滿足的關(guān)系(Hugoniot函數(shù))：

2.3.1.2 彈性稀疏波

2.3.2 塑性變形階段的波系結(jié)構(gòu)

2.3.2.1 塑性激波

2.3.2.2 塑性稀疏波

2.3.3 流體變形階段的波系結(jié)構(gòu)

2.3.3.1 激波

2.3.3.2 稀疏波

2.3.4 彈塑性變形階段的波系結(jié)構(gòu)

當(dāng)固體經(jīng)歷從彈性到塑性的相變過程時，激波的Rankine-Hugoniot曲線或稀疏波等熵線的斜率將發(fā)生間斷，并產(chǎn)生以彈性屈服極限時的狀態(tài)作為臨界點的分裂激波(如圖1(b)所示)。

2.3.4.1 彈性屈服極限狀態(tài)的確定

在求解彈塑性固體的Riemann問題前，首先需要確定固體材料處于壓縮和拉伸屈服極限時的狀態(tài)。根據(jù)Von-Mises屈服條件和(14)式，可得

式中：S表示非線性波后的偏應(yīng)力張量。

彈性極限時的其他物理量可通過相應(yīng)的計算得到：

2.3.4.2 彈塑性激波

相應(yīng)的Hugoniot函數(shù)Φr(σ,ρ)為

式中：

2.3.4.3 彈塑性稀疏波

2.3.5 塑性- 流體階段的波系結(jié)構(gòu)

當(dāng)固體經(jīng)歷從塑性到流體的相變過程時，激波Rankine-Hugoniot曲線或稀疏波等熵線的斜率將發(fā)生間斷，并產(chǎn)生以塑性屈服極限狀態(tài)作為臨界點的分裂激波(類似于圖1(b)所示)。

2.3.5.1 塑性屈服極限狀態(tài)的確定

利用與彈性屈服極限類似的方法，可以計算塑性屈服極限時的各物理量為

2.3.5.2 塑性- 流體激波關(guān)系

(15)

相應(yīng)的Hugoniot函數(shù)為

(16)

式中：

2.3.5.3 塑性- 流體稀疏波

2.3.6 彈性- 塑性- 流體變形階段的波系結(jié)構(gòu)

當(dāng)固體經(jīng)歷從彈性、塑性到流體的相變過程時，激波Rankine-Hugoniot曲線或稀疏波曲線的斜率將發(fā)生間斷，并產(chǎn)生分別以彈性、屈服極限時的狀態(tài)作為臨界點的3個分裂波。

2.3.6.1 彈性- 塑性- 流體激波

相應(yīng)的Hugoniot函數(shù)為

2.3.6.2 彈性- 塑性- 流體稀疏波

(17)

(18)

2.3.7 彈塑性固體Riemann問題的統(tǒng)一形式

定義fr(p,Ur)=u-ur，將2.3.1節(jié)～2.3.6節(jié)中固體各變形階段的激波關(guān)系和稀疏波關(guān)系寫成統(tǒng)一的形式。其中，激波情形滿足：

(19)

(20)

與流體類似，利用牛頓迭代法對(20)式進(jìn)行迭代求解：

對于稀疏波情形，滿足：

聯(lián)立求解稀疏波方程

(21)

和等熵方程

(22)

可完整得到沿稀疏波分支的接觸間斷狀態(tài)。與流體類似，利用自適應(yīng)步長的Runge-Kutta-Fehlberg方法對(21)式和(22)式進(jìn)行聯(lián)立求解，可完成稀疏波分支fr(σ,Ur)的求解。

2.4 流固耦合Riemann問題的求解步驟

綜上所述，本文提出的流固耦合Riemann問題的整體求解步驟如下：

1)提供物質(zhì)界面應(yīng)力的初始估計σ0和整體誤差ε0，其中

2)假設(shè)第m步迭代的值σm已知，確定左、右非線性波的類型：

①如果σm>max {σl,σr}，則兩側(cè)非線性波均為激波；

②如果min {σl,σr}<σm≤max {σl,σr},則一側(cè)非線性波為激波，另一側(cè)為稀疏波；

③如果σm≤min {σl,σr}，則兩側(cè)非線性波均為稀疏波。

3)根據(jù)非線性波的類型和當(dāng)前迭代步的局部求值誤差εm，控制激波分支的非線性代數(shù)方程(10)式、(11)式和(19)式、(20)式的迭代殘差，以及稀疏波分支的常微分方程組(12)式和(21)式、(22)式的局部求值誤差，計算得到當(dāng)前步的fl(σm,Ul)和fr(σm,Ur)。

4)更新得到m+1迭代步時，物質(zhì)界面上的應(yīng)力：

5)當(dāng)應(yīng)力的變化達(dá)到指定的整體誤差ε0時，迭代終止，將得到的收斂解σm近似作為物質(zhì)界面的壓力σ*，否則返回步驟2.

6)計算物質(zhì)界面上的法向速度：

3 多介質(zhì)流動數(shù)值方法

在前期工作[30]的基礎(chǔ)上，進(jìn)一步建立了流固耦合問題的多介質(zhì)流動數(shù)值方法。其中，利用水平集(Level Set)方法來追蹤物質(zhì)界面，即定義距離函數(shù)φ(x,t)的零等值面來表示物質(zhì)界面Γ(t)。利用特征線方法[31]求解φ(x,t)的演化方程

(23)

為保持Level Set函數(shù)的符號距離性質(zhì)，采用顯式正系數(shù)格式[31]求解重新初始化方程

(24)

式中：φ0為重新初始化前φ(x,t)的值；ε為一常數(shù)。

3.1 界面幾何重構(gòu)

圖2 物質(zhì)界面幾何結(jié)構(gòu)示意圖Fig.2 Geometry of material interface

3.2 單元邊界數(shù)值通量

單元邊界數(shù)值通量是指單元Ki,n內(nèi)每種物質(zhì)與相鄰單元的同種物質(zhì)在單元邊界上通過流入流出或相互作用力而帶來的守恒量變化，滿足：

(25)

3.3 物質(zhì)界面數(shù)值通量

物質(zhì)界面數(shù)值通量是界面兩側(cè)的介質(zhì)由于相互作用力帶來的守恒量變化，滿足：

(26)

3.4 守恒量更新

當(dāng)計算得到Ki,n的單元邊界數(shù)值通量和物質(zhì)界面數(shù)值通量后，可對該單元中每種物質(zhì)的守恒量進(jìn)行如下更新：

(27)

4 數(shù)值算例

利用建立的數(shù)值方法對流固Riemann問題、固固Riemann問題、地下強(qiáng)爆炸問題、空中強(qiáng)爆炸問題和高速侵徹等問題進(jìn)行數(shù)值模擬，并將計算結(jié)果與理論和實測數(shù)據(jù)進(jìn)行對比，驗證數(shù)值方法的正確性。

4.1 流固Riemann問題

考慮一個氣體- 彈塑性固體Riemann問題[32]，氣體采用理想氣體狀態(tài)方程來描述，絕熱指數(shù)取γ=2. 采用cm-g-ms單位制，固體的靜水壓力采用Murnagham狀態(tài)方程，偏應(yīng)力采用理想彈塑性模型，狀態(tài)方程參數(shù)為K=2.225×1011Pa，ρ0=7.8 g/cm3，γ=3.7，彈性剪切模量μe=8.53×1010Pa，屈服極限Ye=6.5×109Pa. 計算區(qū)域為[0 cm,1 cm]×[0 cm,0.01 cm]，網(wǎng)格尺寸0.002 5 cm. 初值條件如下：

計算時間t=7.17×10-4ms，計算結(jié)果與文獻(xiàn)[32]的解析結(jié)果對比如圖3所示。從圖3中可知，二者計算結(jié)果比較一致，固體在高壓氣體作用下發(fā)生了彈塑性變形，且前驅(qū)彈性波和后驅(qū)塑性波均被正確的捕捉。

圖3 氣體- 彈塑性固體Riemann問題Fig.3 Gas-elastoplastic Riemann problem

4.2 理想彈塑性固體Riemann問題

計算一個理想彈塑性模型的固固Riemann問題，仍取自文獻(xiàn)[32]。固體采用和4.1節(jié)相同的本構(gòu)模型和材料參數(shù)，計算區(qū)域和網(wǎng)格尺寸同4.1節(jié)，仍采用cm-g-ms單位制。初值條件如下：

計算時間t=6.751×10-5ms. 計算結(jié)果與文獻(xiàn)[32]的解析結(jié)果對比如圖4所示。由圖4對比可知二者計算結(jié)果比較一致，在相界面兩側(cè)的固體中均產(chǎn)生了彈性、塑性激波，且在相界面和激波附近沒有產(chǎn)生非物理震蕩。

圖4 彈塑性- 彈塑性固體Riemann問題 Fig.4 Elastoplastic-elastoplastic Riemann problem

4.3 流體彈塑性固體Riemann問題

計算一個流體彈塑性模型的固固Riemann問題。固體的靜水壓力狀態(tài)方程和參數(shù)與4.1節(jié)相同，偏應(yīng)力采用流體彈塑性模型，且彈性剪切模量μe=8.53×1010Pa，塑性剪切模量μp=4.27×1010Pa，彈性屈服極限Ye=6.5×108Pa，塑性屈服極限Yp=9.75×108Pa. 計算區(qū)域和網(wǎng)格尺寸同4.1節(jié)，初值條件如下：

計算時間t=6.751×10-5ms. 計算結(jié)果如圖5所示，在相界面兩側(cè)的固體中均依次產(chǎn)生了彈性激波、塑性激波和流體激波，且在相界面和激波附近沒有產(chǎn)生非物理震蕩。

圖5 流體彈塑性固體Riemann問題Fig.5 Hydro-elastoplastic Riemann problem

4.4 地下強(qiáng)爆炸應(yīng)力波問題

計算了1 kt TNT當(dāng)量的地下填實強(qiáng)爆炸產(chǎn)生的應(yīng)力波傳播過程。計算區(qū)域為1 000 m，強(qiáng)爆炸產(chǎn)物采用真實氣體狀態(tài)方程[34]，地介質(zhì)采用流體彈塑性材料模型，且μe=58.9×109Pa，μp=18.344×109Pa，Ye=1.05× 108Pa,Yp=1.55×108Pa. 計算得到地下強(qiáng)爆炸不同距離r處的應(yīng)力波峰值速度ur和峰值應(yīng)力σr，并與由實驗結(jié)果擬合的經(jīng)驗公式[33]對比。由圖6可知，數(shù)值計算得到的峰值徑向速度與經(jīng)驗公式吻合較好，峰值徑向應(yīng)力在200 m范圍內(nèi)與經(jīng)驗公式基本符合，200 m以外比經(jīng)驗公式偏高，最大相差約1倍。

圖6 地下強(qiáng)爆炸應(yīng)力波參數(shù)Fig.6 Stress wave parameters from underground intense explosion

分析其原因，主要是由于當(dāng)前的流體彈塑性本構(gòu)方程和模型參數(shù)無法完全反映真實的地下強(qiáng)爆炸應(yīng)力波傳播過程而導(dǎo)致的。在地下發(fā)生強(qiáng)爆炸時，高溫、高壓的強(qiáng)爆炸產(chǎn)物會劇烈壓縮周圍巖石，在不同距離處會依次產(chǎn)生汽化區(qū)、液化區(qū)、破碎區(qū)、塑性區(qū)和彈性區(qū)，且不同區(qū)域的覆蓋范圍通常很大、巖石的力學(xué)性質(zhì)差異較大，使得地下強(qiáng)爆炸效應(yīng)的完全數(shù)值模擬十分困難。其中，在汽化和液化區(qū)，由于沖擊波強(qiáng)度非常高，巖石吸收大量能量，并很快轉(zhuǎn)變?yōu)榱黧w，此時材料強(qiáng)度對應(yīng)力波的影響基本可以忽略，采用流體彈塑性模型能夠近似模擬該階段的巖石動力學(xué)特性。在破碎區(qū)和塑性區(qū)，由于巖石是一種典型的各向異性材料，通常包含著各種尺度的缺陷，當(dāng)受到一定強(qiáng)度的載荷作用時，巖石容易產(chǎn)生裂紋并導(dǎo)致破裂。此外，由于巖石樣品的不均質(zhì)性，也使實驗數(shù)據(jù)也具有很大的離散性，使得其本構(gòu)關(guān)系的建立十分困難。

如何建立合理的巖石材料本構(gòu)模型，能較為準(zhǔn)確地模擬大跨度壓力范圍內(nèi)爆炸產(chǎn)生的近、遠(yuǎn)區(qū)力學(xué)效應(yīng)，長期以來一直是一個極具挑戰(zhàn)性的課題。本算例重點關(guān)注的部分在于爆炸產(chǎn)物與巖石之間的流固耦合作用過程，后續(xù)將在現(xiàn)有基礎(chǔ)上進(jìn)一步開展考慮斷裂、損傷的巖石本構(gòu)模型研究，提高地下爆炸遠(yuǎn)區(qū)應(yīng)力波參數(shù)計算的準(zhǔn)確度。

4.5 空中強(qiáng)爆炸沖擊波問題

計算一個當(dāng)量為1 kt TNT、爆炸高度為50 m的空中強(qiáng)爆炸沖擊波傳播問題。強(qiáng)爆炸產(chǎn)物采用等溫、等壓球模型，取爆炸總當(dāng)量的85%作為強(qiáng)爆炸的力學(xué)初始能量[34-35]?？諝獠捎美硐霘怏w狀態(tài)方程，初始壓力為101.3 kPa，初始密度為1.29 kg/m3. 計算區(qū)域在徑向距離r=0 m和距地面高度z=0 m設(shè)為固壁反射邊界條件，其余邊界設(shè)為無反射邊界條件。

圖7 空中強(qiáng)爆炸典型時刻的壓力云圖Fig.7 Pressure contours of air blast at typical time

計算得到典型時刻的沖擊波壓力等值線如圖7所示，從中可看出，當(dāng)空中強(qiáng)爆炸產(chǎn)生的沖擊波傳播到地面時，最早在爆心投影點發(fā)生正反射，然后以逐漸增大的入射角發(fā)生斜反射，產(chǎn)生雙激波結(jié)構(gòu)的規(guī)則反射(見圖7(a))。隨著沖擊波繼續(xù)向外傳播，入射角不斷增大。當(dāng)入射角增大到臨界角附近時將發(fā)生馬赫反射，此時反射波陣面和入射波陣面的交點離開地面，形成接近垂直于地面的馬赫桿(見圖7(b))。隨著沖擊波進(jìn)一步向外傳播，馬赫桿快速增長(見圖7(c)～圖7(e))，并最終形成半球反射，此時沖擊波以接近球面波的形式繼續(xù)向外傳播。

圖8給出了數(shù)值計算得到的地面沖擊波載荷分布(峰值超壓p和沖量I)與實測結(jié)果[36]的對比情況。由圖8對比可知，計算得到的地面沖擊波參數(shù)與實測結(jié)果符合較好，最大誤差不超過20%，驗證了本文數(shù)值方法的正確性。

圖8 地面不同距離處的沖擊波參數(shù)Fig.8 Parameters of blast wave on the ground at different distances

4.6 高速侵徹問題

計算一個高速沖擊下的彈塑性固體侵徹問題。模型如圖9(a)所示，模型中上方、直徑為0.5 cm的圓柱體鋼材料彈體以700 m/s速度撞擊下方靜止、直徑為4 cm、厚度為1 cm的鋁材料靶體。兩種金屬材料的靜水壓力均采用剛性氣體狀態(tài)方程，偏應(yīng)力部分采用理想彈塑性材料模型。其中，彈體(鋼)的材料參數(shù)為：絕熱指數(shù)γ=4.075,初始密度ρ0=7 840 kg/m3,彈性剪切模量μe=75.8 GPa，屈服強(qiáng)度Ye=1.6 GPa；靶體(鋁)的材料參數(shù)為：絕熱指數(shù)γ=2.75，初始密度ρ0=2 790 kg/m3，彈性剪切模量μe=27.4 GPa，屈服強(qiáng)度Ye=0.34 GPa.

圖9 高速侵徹的密度云圖Fig.9 Density contours of high speed penetration problem

圖9給出了不同時刻彈體和靶體內(nèi)的密度云圖。由圖9計算結(jié)果表明，當(dāng)彈體侵徹靶體時，在彈體和靶體中均產(chǎn)生強(qiáng)烈的沖擊波，并分別沿著各自介質(zhì)向外傳播。由于彈體比靶體具有更高的密度和強(qiáng)度，在彈體侵徹靶體過程中，彈體的形狀變化較小，而靶體發(fā)生了劇烈的變形。

表2 侵徹深度和侵徹半徑的數(shù)值結(jié)果與實驗結(jié)果[37]的對比

5 結(jié)論

本文通過爆炸與沖擊動力學(xué)問題的數(shù)值模擬方法與應(yīng)用，得出了以下主要結(jié)論：

1)提出一種能夠處理高度非線性狀態(tài)方程和彈塑性相變本構(gòu)模型的多介質(zhì)Riemann問題通用求解方法，能有效提高物質(zhì)界面上各物理量的計算精度。

2)結(jié)合基于Euler坐標(biāo)系的互不相溶、具有清晰銳利界面的可壓縮多介質(zhì)流動數(shù)值方法，建立了一套能夠模擬具有高密度比、高壓力比以及復(fù)雜狀態(tài)方程的流固耦合、固固耦合問題的多介質(zhì)計算體系。

3)依次對一維流固Riemann問題、固固Riemann問題、地下強(qiáng)爆炸問題、空中強(qiáng)爆炸問題和高速侵徹等問題進(jìn)行數(shù)值模擬，計算結(jié)果與理論分析和實測數(shù)據(jù)比較一致，表明數(shù)值方法能夠有效應(yīng)用于地下強(qiáng)爆炸和高速侵徹等具有多介質(zhì)、大變形的實際工程應(yīng)用問題。

參考文獻(xiàn)(References)

[1] KAPAHI A, SAMBASIVAN S, UDAYKUMAR H. A three-dimensional sharp interface Cartesian grid method for solving high speed multi-material impact, penetration and fragmentation problems[J]. Journal of Computational Physics, 2013, 241: 308-332.

[2] SAMBASIVAN S, KAPAHI A, UDAYKUMAR H. Simulation of high speed impact, penetration and fragmentation problems on locally refined Cartesian grids[J]. Journal of Computational Physics, 2013, 235: 334-370.

[3] SAMBASIVAN S, UDAYKUMAR H. A sharp interface method for high-speed multi-material flows: strong shocks and arbitrary materialpairs[J]. International Journal of Computational Fluid Dynamics, 2011, 25(3): 139-162.

[4] CHEN Q Y, LIU K X. A high-resolution Eulerian method for numerical simulation of shaped charge jet including solid-fluid coexistence and interaction[J]. Computers & Fluids, 2012, 56: 92-101.

[5] CHEN Q Y, WANG J T, LIU K X. Improved CE/SE scheme with particle level set method for numerical simulation of spall fracture due to high-velocity impact[J]. Journal of Computational Physics, 2010, 229(19): 7503-7519.

[6] 陳千一. 層裂與成型裝藥金屬射流問題的數(shù)值分析[D]. 北京：北京大學(xué), 2011.

CHEN Q Y. Numerical analysis of spall fracture and shaped charge jet[D]. Beijing: Peking University, 2011. (in Chinese)

[7] 丁建許. 爆炸與沖擊問題的高精度界面處理及數(shù)值研究[D]. 北京：北京理工大學(xué), 2017.

DING J X. The high order inerfacial treatment and numerical investigation of explosion and impact problem[D]. Beijing: Beijing Institute of Technology, 2017. (in Chinese)

[8] 王仲琦. 面向?qū)ο蟮谋W(xué) Euler 型多物質(zhì)數(shù)值方法及其應(yīng)用研究[D]. 北京：北京理工大學(xué), 2000.

WANG Z Q. Object-oriented Euler multi-material numerical met-hods on and its applications[D]. Beijing: Beijing Institute of Technology, 2000. (in Chinese)

[9] BANKS J R. On exact conservation for the Euler equations with complex equations of state[J]. Communications in Computational Physics, 2010, 8(5): 995.

[10] KAMM J. Solution of the 1D Riemann problem with a general EOS in ExactPack[C]∥Proceedings of the 4th ASME Conference on Verification and Validation of Simulations. Las Vegas,NV, US：LANL, 2015.

[11] DESPRES B. A geometrical approach to nonconservative shocks and elastoplastic shocks[J]. Archive for Rational Mechanics and Analysis, 2007, 186(2): 275-308.

[12] GAVRILYUK S, FAVRIE N, SAUREL R. Modelling wave dynamics of compressible elastic materials[J]. Journal of Computational Physics, 2008, 227(5): 2941-2969.

[13] KAMM J. FLAG simulations of the elasticity test problem of Gavrilyuk et al[R]. Los Alamos, NV,US: Los Alamos National Laboratory,2014.

[14] LIN X, BALLMANN J. A Riemann solver and a second-order Godunov method for elasticplastic wave propagation in solids[J]. International Journal of Impact Engineering, 1993, 13(3): 463-478.

[15] TRANGENSTEIN J, COLELLA P. A higher-order Godunov method for modeling finite deformation in elastic-plastic solids[J]. Communications on Pure and Applied Mathematics, 1991, 44(1): 41-100.

[16] ABUZYAROV M, BAZHENOV V, KOTOV V. A Godunov-type method in dvnamics of elastoplastic media[J]. Computational Mathematics and Mathematical Physics, 2000, 40(6): 900-913.

[17] BAZHENOV V, KOTOV V. Modification of Godunov’s numerical scheme for solving problems of pulsed loading of soft soils[J]. Journal of Applied Mechanics and Technical Physics, 2002, 43(4): 603-611.

[18] MENSHOV I, MISCHENKO A, SEREJKIN A. Numerical mo-deling of elastoplastic flows by the Godunov method on moving Eulerian grids[J]. Mathematical Models and Computer Simulations, 2014, 6(2): 127-141.

[19] TANG H S, HUANG D. A second-order accurate capturing scheme for 1D inviscid flows of gas and water with vacuum zones[J]. Journal of Computational Physics, 1996, 128(2): 301-318.

[20] LIU T G, XIE W F, KHOO B C. The modified ghost fluid method for coupling of fluid and structure constituted with hydro-elasto-plastic equation of state[J]. SIAM Journal on Scientific Computing, 2008, 30(3): 1105-1130.

[21] GAO S, LIU T G. 1D exact elastic-perfectly plastic solid Riemann solver and its multimaterial application[J]. Advances in Applied Mathematics and Mechanics, 2017, 9(3): 621-650.

[22] GAO S, LIU T G, YAO C B. A complete list of exact solutions for one-dimensional elasticperfectly plastic solid Riemann problem without vacuum[J]. Communications in Nonlinear Science and Numerical Simulation, 2018, 63: 205-227.

[23] BERJAMIN H. Analytical solution to 1D nonlinear elastodynamics with general constitutive laws[J]. Wave Motion, 2017, 74: 35-55.

[24] HATTORI H. The Riemann problem for thermoelastic materials with phase change[J]. Journal of Differential Equations, 2004, 205(1): 229-252.

[25] FENG Z, KABOUDIAN A, RONG J. The simulation of compressible multi-fluid multi-solid interactions using the modified ghost method[J]. Computers & Fluids, 2017, 154(5): 12-26.

[26] FENG Z, RONG J, KABOUDIAN A. The modified ghost method for compressible multimedium interaction with elastic-plastic solid[J]. Communications in Computational Physics, 2017, 22(5): 1258-1285.

[27] TORO E F. Riemann solvers and numerical methods for fluid dynamics[M]. Berlin, Germany: Springer, 2009: 102-200.

[28] FEHLBERG E. Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems: R-315[R]. Washington, DC,US: NASA, 1969.

[29] MATHEWS J H, FINK K D. Numerical methods using MATLAB[M]. Upper Saddle River, NJ,US: Pearson Prentice Hall, 2004.

[30] GUO Y H, LI R, YAO C B. A numerical method on Eulerian grids for two-phase compressible flow[J]. Advances in Applied Mathematics and Mechanics, 2016, 8(2):187-212.

[31] DI Y N, LI R, TANG T. Level set calculations for incompressible two-phase flows on a dynamically adaptive grid[J]. Journal of Scientific Computing, 2007, 31(1):75-98.

[32] LIU T G, XIE W F, KHOO B C. The modified ghost fluid method for coupling of fluid and structure constituted with hydro elastoplastic equation of state [J]. SIAM Journal on Scientific Computing, 2008, 30(3): 1105-1130.

[33] 郝寶田. 地下核爆炸及其作用[M]. 北京：國防工業(yè)出版社, 2002.

HAO B T. Underground explosions and application [M]. Beijing: National Defense Industry Press, 2002. (in Chinese)

[34] 姚成寶, 王宏亮, 浦錫鋒,等. 空中強(qiáng)爆炸沖擊波地面反射規(guī)律數(shù)值模擬研究[J]. 爆炸與沖擊, 2019, 39(11): 112201-1-112201-8.

YAO C B, WANG H L, PU X F, et al. Numerical simulation of intense blast wave reflected on rigid ground[J]. Explosion and Shock Waves, 2019, 39(11): 112201-1-112201-8. (in Chinese)

[35] 喬登江. 核爆炸物理概論[M]. 北京：國防工業(yè)出版社, 2003:51-55.

QIAO D J. An introduction to nuclear explosion physics[M]. Beijing: National Defense Industry Press, 2003:51-55. (in Chinese)

[36] GLASSTONE S, DOLAN P J. The effects of nuclear weapons [M]. US: United States Department of Defense and the Energy Research and Development Administration, 1977: 453-501.

[37] 江松，王瑞利. 多介質(zhì)可壓縮流體力學(xué)的若干測試基準(zhǔn)問題[R]. 北京：北京應(yīng)用物理與計算數(shù)學(xué)研究所, 2013.

JIANG S, WANG R L. Several benchmark problems of compressible multi-material fluid dynamics[R]. Beijing:Institute of Applied Physics and Computational Mathematics, 2013. (in Chinese)