使用新型物態(tài)方程的超高速碰撞物質(zhì)點法模擬*

李依瀟，王生捷

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

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

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

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

1 物質(zhì)點法基本原理

1.1 控制方程

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

質(zhì)量守恒：

動量方程：

能量方程：

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

1.2 離散形式

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

質(zhì)點坐標：

質(zhì)點位移：

質(zhì)點位移導數(shù)：

結(jié)點質(zhì)量：

結(jié)點動量：

結(jié)點內(nèi)力：

式中：下標p表示質(zhì)點攜帶的變量，下標I表示背景網(wǎng)格結(jié)點攜帶的變量；i=r，z；j=r，z；np表示與結(jié)點I相鄰的質(zhì)點總數(shù)。

1.3 求解格式

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

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

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

2 材料模型

材料模型描述了材料在外力作用下的響應(yīng)。材料模型包括本構(gòu)模型、失效模型和物態(tài)方程。

本構(gòu)模型用于描述材料的偏應(yīng)力與偏應(yīng)變的關(guān)系，偏應(yīng)力更新算法如下式所示[2]：

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

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

3 新型物態(tài)方程的建立

3.1 分子動力學基本原理

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

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

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

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

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

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

體系的冷壓pc可表示為：

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

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

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

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

3.2 新型物態(tài)方程

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

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

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

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

當V≥VJ（氣相）時：

3.3 新型物態(tài)方程計算流程

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

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

（2）根據(jù)物質(zhì)點法解得的當前時刻質(zhì)點內(nèi)能E，與上一步已求得的冷能Ec，做差求得質(zhì)點在當前時刻的熱能Ek；

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

4 鋁彈鋁靶超高速碰撞數(shù)值模擬結(jié)果

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

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

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

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

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

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

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

5 結(jié) 論

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