使用新型物態方程的超高速碰撞物質點法模擬*

李依瀟，王生捷

（1. 中國航天科工集團第二研究院研究生院，北京 100854；2. 北京機械設備研究所，北京 100854）

超高速碰撞是指“材料強度遠小于自身所受慣性力”[1]的碰撞。超高速碰撞現象廣泛存在于航天、兵器、天體物理等領域。超高速碰撞問題具有復雜性，涉及固體力學、流體力學、熱力學等多門學科，理論研究存在較大困難，實驗研究受到成本、技術等方面的限制，因此數值模擬成為研究超高速碰撞問題的重要手段。

物質在超高速碰撞過程中大變形、高應變率、破碎等表現，給數值模擬帶來巨大挑戰。數值模擬的核心是數值計算方法，即對連續介質的離散近似。拉格朗日法將計算網格固連在物體上，在求解超高速碰撞問題時會發生網格畸變，無法有效模擬材料破碎等現象；歐拉法將計算網格固定在空間中，材料相對網格運動，不存在網格畸變問題，但控制方程存在對流項，求解較困難；針對以上兩種方法的優缺點，任意拉格朗日-歐拉法(arbitrary Lagrange-Euler, ALE)將兩者結合使用，可解決一大批只用拉格朗日法或歐拉法解決不了的問題。為克服網格法存在的不足，無網格法應運而生。光滑粒子流體動力學(smoothed particle hydrodynamics, SPH)方法將計算域離散為一系列有相互作用的粒子，這些粒子承載質量、速度、能量等物理量，通過核函數發生相互作用。SPH 方法在超高速碰撞問題數值模擬中得到廣泛應用，但仍存在固有的缺陷，如計算效率低、拉伸不穩定、邊界條件難以處理等。物質點法(material point method, MPM)是一種采用拉格朗日質點和歐拉網格雙重描述的數值計算方法，綜合了網格法和無網格法的優點，非常適合分析超高速碰撞問題。物質點法的基本思想是將連續體離散成一組帶有質量的質點，并在物質運動區域建立背景網格，質點攜帶質量、速度、應力、應變、能量等物質信息，網格僅用于動量方程的求解與空間導數的計算，在每一個時間步中，質點與網格完全固連，在時間步結束時丟棄已變形的背景網格[2-3]。

在超高速碰撞問題中，材料的強度與其所受的壓力相比可以忽略不計，材料處于流動狀態，材料中的壓力需要使用物態方程來計算。物態方程用于描述物質壓力、內能、體積的函數關系，在物質點法模擬的每一時間步內，通過將所求得的質點體積、內能代入物態方程，解得當前時刻質點的壓力。Mie-Grüneisen 物態方程形式簡單，可以很好地描述絕大多數金屬固體在沖擊載荷作用下的熱力學行為，常用于沖擊動力學問題的數值模擬[2,4]。Tillotson 物態方程考慮材料的凝聚態和氣態，沒有考慮熔化，主要優點是適合于高壓區域[5]，在沖擊動力學問題數值模擬中的應用也很廣泛，碰撞速度高于4 km/s 的問題就應該使用Tillotson 物態方程計算[6]。由于超高速碰撞問題涉及材料液化、汽化等相變現象，在數值模擬中使用Mie-Grüneisen 物態方程或Tillotson 物態方程，會導致計算結果與實驗結果存在偏差。本文中結合Grover 定標律方程與分子動力學方法，計算得到金屬鋁的新型物態方程，并代入MATLAB 自編柱坐標物質點法計算程序，將計算所得碎片云與實驗結果，使用Mie-Grüneisen 物態方程、Tillotson 物態方程的計算結果進行對比，證明新型物態方程的有效性。

1 物質點法基本原理

1.1 控制方程

在超高速碰撞問題中，碰撞點附近區域的壓力遠高于材料本身強度，可以認為該區域的材料是可壓縮流體。材料在超高速碰撞過程中經歷高應變率的變形，變形過程很短暫，可忽略熱傳導效應[7]。軸對稱沖擊動力學分析常使用柱坐標形式的更新拉格朗日格式控制方程[8-9]。

質量守恒：

動量方程：

能量方程：

控制方程是一組偏微分方程，求解此類方程的方法分兩類：一類是直接求解；另一類是先建立和原微分方程及其定解條件等價的弱形式，再以此為基礎建立近似解法[2]。物質點法的求解采用后一種方法。

1.2 離散形式

物質點法將連續體離散為一系列質點，質點的運動代表了物體的運動。物質點法采用質點積分，即將虛功方程中的各項積分轉化為被積函數在各質點處的值與該質點所代表的體積乘積的和。在求解動量方程時，質點和背景網格完全固連，隨背景網格一起運動，質點和網格結點之間通過形函數NIp建立信息的映射。質點攜帶的位移、應力等信息的導數，可以由結點信息插值得到。質點與結點之間的信息映射公式[2]如下。

質點坐標：

質點位移：

質點位移導數：

結點質量：

結點動量：

結點內力：

式中：下標p表示質點攜帶的變量，下標I表示背景網格結點攜帶的變量；i=r，z；j=r，z；np表示與結點I相鄰的質點總數。

1.3 求解格式

在物質點法中，物體的所有物質信息由質點攜帶，背景網格結點不存儲任何物質信息。在求解動量方程時，需要將質點在當前時刻的質量、動量、應力等信息通過形函數映射至背景網格結點，在背景網格上進行求解。之后，將網格結點的速度、位置增量映射回質點，使質點速度、位置得到更新。

物質點法的顯式求解方法按照應力更新方式的不同分為先更新應力(update stress first，USF)和后更新應力(update stress last, USL)和改進型USL (modified USL, MUSL)3 種，其中USF 在處理超高速碰撞問題時具有效率高、精度高的特點[10]。

物質點法的USF 求解格式如下：（1）將質點的質量和動量通過形函數映射到背景網格結點，求得結點質量和動量；（2）對結點動量施加邊界條件，進行修正；（3）由結點的質量和動量求得結點速度，由此計算質點速度梯度、應變增量和旋量增量，更新質點的體積、應力偏量、內能；（4）通過物態方程求解質點的壓力；（5）計算背景網格結點的合力，并根據邊界條件進行修正，更新背景網格結點動量；（6）將背景網格結點位置、速度的變化量映射回質點，更新質點的位置和速度；（7）丟棄已變形的網格，在下一時間步使用未變形的新網格。

2 材料模型

材料模型描述了材料在外力作用下的響應。材料模型包括本構模型、失效模型和物態方程。

本構模型用于描述材料的偏應力與偏應變的關系，偏應力更新算法如下式所示[2]：

材料在超高速碰撞過程中的屈服應力可用Johnson-Cook 模型表示[2]：

超高速碰撞問題涉及材料的沖擊破壞，本文中使用聯合失效模型進行描述。當質點的拉應力大于給定值或溫度高于熔點時，質點失效。失效質點的應力偏量為零，且不能承受拉應力。

3 新型物態方程的建立

3.1 分子動力學基本原理

分子動力學模擬是一種用來計算多體體系平衡和傳遞性質的一種確定性方法。分子動力學方法將原子視為經典粒子，體系的多體相互作用由包含經驗參數的解析函數直接給出[11]，粒子的運動遵循牛頓運動方程：

式中：mi為粒子i的質量，ri為粒子i的位置矢量，fi為粒子i所受的力，Ep為勢函數。

在分子動力學模擬中，勢函數表征原子間的相互作用，嵌入原子(embedded atom method，EAM)勢在描述金屬的性質時與實際情況能夠較好地吻合。EAM 勢的基本思想是將組成體系的原子看成一個個嵌入由其他所有原子構成的有效介質中的客體原子的集合，從而將系統的總能量表達為嵌入能和相互作用的對勢之和[12]。常用的EAM 勢有Johnson 分析型EAM 勢、Cai-Ye EAM、Zhou EAM 等，本文中采用Cai-Ye EAM 勢[12]。EAM 勢的總能表達式為：

式中：Fi為嵌入項，Φ 為對勢項，N為體系粒子數，ρi為粒子i的電荷密度，ri為粒子i的位置矢量，rj為粒子j的位置矢量。

在分子動力學方法中，粒子的微觀量與系統的宏觀量通過統計物理聯系起來。體系的總能量E由下式給出：

式中：Ek為體系的熱能，Ec為體系的冷能，ec為質量冷能，vi為粒子i的速率。

體系的冷壓pc可表示為：

式中：V為體系的體積，rij為粒子i與粒子j間的距離，fij為粒子i與粒子j間的作用力。

實際使用中，首先利用分子動力學方法計算得到冷能、冷壓與體積比的關系，如圖1～2 所示。在超高速碰撞問題的物質點法模擬中，每一時間步內，根據質點當前時刻的體積，利用插值法求得質點的冷能和冷壓。

圖1 冷能與體積比的關系Fig. 1 Relation between cold energy and volume ratio

圖2 冷壓與體積比的關系Fig. 2 Relation between cold pressure and volume ratio

3.2 新型物態方程

新型物態方程基于Grover 定標律方程和分子動力學方法構建，利用分子動力學方法計算冷能、冷壓，避免了使用Hugoniot 曲線數據所導致的物態方程形式繁瑣。新型物態方程在固-液相區采用Grover 定標律方程的形式，在氣相區采用分子動力學中常用的維里方程的形式[13]。針對超高速碰撞問題的物質點法模擬，單個物質點處于單相區或混合相區對整體的影響可忽略不計，因此新型物態方程不考慮混合相區。新型物態方程的形式如下。

當V＜VJ且Ek＜Em（固相）時：

當V＜VJ且Em≤Ek≤EG（液相）時：

當V＜VJ且Ek＞EG（熱液相）時：

當V≥VJ（氣相）時：

3.3 新型物態方程計算流程

在使用新型物態方程的超高速碰撞物質點法模擬中，分子動力學方法用于求解質點的冷能和冷壓，質點的熱能用于判斷所處的相區，具體計算流程如下：

（1）根據物質點法解得的當前時刻質點體積V，得到當前時刻體積V與初始體積V0之比，利用插值法求得質點在當前時刻的冷能Ec和冷壓pc，并將體積比代入公式求得Tm、TG、Em、EG；

（2）根據物質點法解得的當前時刻質點內能E，與上一步已求得的冷能Ec，做差求得質點在當前時刻的熱能Ek；

（3）根據V、VJ、Ek、Em、EG判斷質點所處的相區，利用式(16)～(19)中的能量關系式求得質點溫度T，將冷壓pc、溫度T代入所屬相區的壓力關系式，求得質點壓力p。

4 鋁彈鋁靶超高速碰撞數值模擬結果

直徑為9.5 mm 的球形鋁彈丸以6 640 m/s 的超高速撞擊厚度為2.2 mm 的鋁靶。碰撞6.6 μs 后的實驗結果[15]如圖3(a)所示, 碎片云長37.7 mm，寬34.2 mm。

彈體和靶體的離散尺寸為0.05 mm，質點總數為58 010，背景網格采用0.1 mm 進行計算，分別使用新型物態方程、Mie-Grüneisen 物態方程、Tillotson 物態方程進行計算，碰撞6.6 μs 后的計算結果如圖3(b)～(d)所示。

圖3 數值模擬結果與實驗結果[15]的對比Fig. 3 Comparison between experimental[15] and numerical results

使用新型物態方程計算所得碎片云(圖3(b))長36.8 mm，寬33.3 mm，與實驗結果[15]符合較好。計算所得碎片云前部呈現出與實驗結果相符的階梯狀，碎片云后部反濺部分的形態也與實驗結果符合較好。

使用Mie-Grüneisen 物態方程計算所得碎片(圖3(c))長36.8 mm，寬32.4 mm。計算所得碎片云前部呈現出階梯狀，但與實驗結果[15]存在一定差異，碎片云后部反濺部分的形態與實驗結果[15]符合較好。

使用Tillotson 物態方程計算所得碎片(圖3(d))長38.7 mm，寬34.1 mm。計算所得碎片云前部沒有出現階梯狀，與實驗結果[15]不符；碎片云后部反濺部分的形態與實驗結果[15]符合較好。

可見，在碎片云尺寸方面，使用3 種物態方程計算結果均與實驗結果[15]符合較好。碎片云形態方面，使用新型物態方程計算結果與實驗結果最接近；使用Mie-Grüneisen 物態方程計算所得的碎片云前部雖然呈現出階梯狀，但與實驗結果[15]仍有差距；使用Tillotson 物態方程計算所得的碎片云前部沒有出現階梯狀，形態與實驗結果[15]差距最大。

5 結 論

（1）使用柱坐標物質點法能夠對軸對稱沖擊動力學問題進行準確度較高的數值模擬，計算所得碎片云尺寸與實驗結果符合較好。（2）盡管Mie-Grüneisen 物態方程只適用于描述金屬固體在沖擊載荷下的熱力學行為，但對于碰撞速度達到6 640 m/s 的超高速碰撞問題，使用Mie-Grüneisen 物態方程計算所得的碎片云尺寸與實驗結果相比誤差較小，且形態優于使用Tillotson 物態方程的計算結果，因此不必換用Tillotson 物態方程。（3）新型物態方程基于Grover 定標律方程和分子動力學方法建立，可以有效處理材料的相變問題，適用于超高速碰撞問題的分析。使用新型物態方程計算所得的碎片云形態、尺寸與實驗結果吻合很好，成功再現鋁彈鋁靶超高速碰撞碎片云前部的階梯狀形態。（4）使用新型物態方程和Mie-Grüneisen 物態方程計算所得的碎片云尺寸小于實驗結果和使用Tillotson 物態方程的計算結果，在使用新型物態方程和Mie-Grüneisen 物態方程對航天器防護結構超高速碰撞進行數值模擬時，有可能因碎片云分布范圍小于實際情況而導致不正確的穿透，但有助于提高彈道極限方程的保守性。