基于多介質Riemann問題的流體- 固體耦合數值方法及其在爆炸與沖擊問題中的應用

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

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

0 引言

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

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

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

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

1 物理模型

1.1 控制方程

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

(1)

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

1.2 狀態方程和本構模型

除了描述物質(流體或固體)的質量、動量和能量守恒的控制方程外，物質動力學行為的完整描述還需要提供描述自身熱力學狀態的本構方程(狀態方程)。為后續進行統一分析，將流體狀態方程和固體靜水壓力方程寫成(2)式的統一形式：

p=ωρe+h(ρ)，

(2)

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

梯恩梯(TNT)爆轟產物采用JWL狀態方程，如(3)式所示：

(3)

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

表1 JWL狀態方程參數

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

p=(γ-1)ρe，

(4)

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

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

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

(5)

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

Murnagham狀態方程是一種線性化的Grüneisen狀態方程，可用于描述金屬等材料在受到沖擊時的動態響應問題，

(6)

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

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

(7)

2 流固耦合Riemann問題求解

2.1 流固耦合Riemann問題

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

對于流體與彈塑性固體耦合的一維多介質Riemann問題，當不考慮黏性時，Riemann問題的解由初值條件以及相應介質的狀態方程或本構方程決定。不失一般性，假設流體位于物質界面的左側，彈塑性固體位于物質界面的右側，此時一維流固耦合的多介質Riemann問題定義為

Uτ+Ψξ(U)=0，

(8)

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

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

2.2 流體中的波系結構

對于流體，當p>pl時，非線性波的類型為激波，且激波前、后的物理量狀態滿足Rankine-Hugoniot關系：

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

式中：c表示聲速，且

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

(9)

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

2.2.1 激波

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

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

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

(10)

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

(11)

2.2.2 稀疏波

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

(12)

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

2.3 固體中的波系結構

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

2.3.1 彈性變形階段的波系結構

2.3.1.1 彈性激波

(13)

且靜水壓力p和偏應力s分別滿足

(14)

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

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

2.3.1.2 彈性稀疏波

2.3.2 塑性變形階段的波系結構

2.3.2.1 塑性激波

2.3.2.2 塑性稀疏波

2.3.3 流體變形階段的波系結構

2.3.3.1 激波

2.3.3.2 稀疏波

2.3.4 彈塑性變形階段的波系結構

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

2.3.4.1 彈性屈服極限狀態的確定

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

式中：S表示非線性波后的偏應力張量。

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

2.3.4.2 彈塑性激波

相應的Hugoniot函數Φr(σ,ρ)為

式中：

2.3.4.3 彈塑性稀疏波

2.3.5 塑性- 流體階段的波系結構

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

2.3.5.1 塑性屈服極限狀態的確定

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

2.3.5.2 塑性- 流體激波關系

(15)

相應的Hugoniot函數為

(16)

式中：

2.3.5.3 塑性- 流體稀疏波

2.3.6 彈性- 塑性- 流體變形階段的波系結構

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

2.3.6.1 彈性- 塑性- 流體激波

相應的Hugoniot函數為

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

(17)

(18)

2.3.7 彈塑性固體Riemann問題的統一形式

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

(19)

(20)

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

對于稀疏波情形，滿足：

聯立求解稀疏波方程

(21)

和等熵方程

(22)

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

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

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

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

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

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

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

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

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

4)更新得到m+1迭代步時，物質界面上的應力：

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

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

3 多介質流動數值方法

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

(23)

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

(24)

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

3.1 界面幾何重構

圖2 物質界面幾何結構示意圖Fig.2 Geometry of material interface

3.2 單元邊界數值通量

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

(25)

3.3 物質界面數值通量

物質界面數值通量是界面兩側的介質由于相互作用力帶來的守恒量變化，滿足：

(26)

3.4 守恒量更新

當計算得到Ki,n的單元邊界數值通量和物質界面數值通量后，可對該單元中每種物質的守恒量進行如下更新：

(27)

4 數值算例

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

4.1 流固Riemann問題

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

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

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

4.2 理想彈塑性固體Riemann問題

計算一個理想彈塑性模型的固固Riemann問題，仍取自文獻[32]。固體采用和4.1節相同的本構模型和材料參數，計算區域和網格尺寸同4.1節，仍采用cm-g-ms單位制。初值條件如下：

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

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

4.3 流體彈塑性固體Riemann問題

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

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

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

4.4 地下強爆炸應力波問題

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

圖6 地下強爆炸應力波參數Fig.6 Stress wave parameters from underground intense explosion

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

如何建立合理的巖石材料本構模型，能較為準確地模擬大跨度壓力范圍內爆炸產生的近、遠區力學效應，長期以來一直是一個極具挑戰性的課題。本算例重點關注的部分在于爆炸產物與巖石之間的流固耦合作用過程，后續將在現有基礎上進一步開展考慮斷裂、損傷的巖石本構模型研究，提高地下爆炸遠區應力波參數計算的準確度。

4.5 空中強爆炸沖擊波問題

計算一個當量為1 kt TNT、爆炸高度為50 m的空中強爆炸沖擊波傳播問題。強爆炸產物采用等溫、等壓球模型，取爆炸總當量的85%作為強爆炸的力學初始能量[34-35]。空氣采用理想氣體狀態方程，初始壓力為101.3 kPa，初始密度為1.29 kg/m3. 計算區域在徑向距離r=0 m和距地面高度z=0 m設為固壁反射邊界條件，其余邊界設為無反射邊界條件。

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

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

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

圖8 地面不同距離處的沖擊波參數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的鋁材料靶體。兩種金屬材料的靜水壓力均采用剛性氣體狀態方程，偏應力部分采用理想彈塑性材料模型。其中，彈體(鋼)的材料參數為：絕熱指數γ=4.075,初始密度ρ0=7 840 kg/m3,彈性剪切模量μe=75.8 GPa，屈服強度Ye=1.6 GPa；靶體(鋁)的材料參數為：絕熱指數γ=2.75，初始密度ρ0=2 790 kg/m3，彈性剪切模量μe=27.4 GPa，屈服強度Ye=0.34 GPa.

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

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

表2 侵徹深度和侵徹半徑的數值結果與實驗結果[37]的對比

5 結論

本文通過爆炸與沖擊動力學問題的數值模擬方法與應用，得出了以下主要結論：

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

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

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

參考文獻(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] 陳千一. 層裂與成型裝藥金屬射流問題的數值分析[D]. 北京：北京大學, 2011.

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

[7] 丁建許. 爆炸與沖擊問題的高精度界面處理及數值研究[D]. 北京：北京理工大學, 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] 王仲琦. 面向對象的爆炸力學 Euler 型多物質數值方法及其應用研究[D]. 北京：北京理工大學, 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]. 北京：國防工業出版社, 2002.

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

[34] 姚成寶, 王宏亮, 浦錫鋒,等. 空中強爆炸沖擊波地面反射規律數值模擬研究[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]. 北京：國防工業出版社, 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] 江松，王瑞利. 多介質可壓縮流體力學的若干測試基準問題[R]. 北京：北京應用物理與計算數學研究所, 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)