基于Lagrange 及SPH 算法的花崗巖侵徹數(shù)值模擬

靳紹虎，劉科偉，黃 進，楊家彩，靳少博

（中南大學(xué)資源與安全工程學(xué)院，湖南 長沙 410083）

自海灣戰(zhàn)爭以來，隨著精確制導(dǎo)武器和小型鉆地核武器的不斷發(fā)展，地下防護空間的開發(fā)利用也在世界各國的國防建設(shè)中得到普遍發(fā)展，并受到高度重視[1]。巖石作為一種天然的防護材料，具有強度高、硬度大、彈性模量大等優(yōu)點，被廣泛應(yīng)用于防護工程和重要軍事設(shè)施的構(gòu)筑。在過去的幾十年里，國內(nèi)外專家學(xué)者對于巖石靶體的侵徹效應(yīng)進行了深入研究，依據(jù)不同研究手段，分為實驗[2-5]、理論[6-9]和數(shù)值模擬研究[10-13]。

近年來，隨著計算機技術(shù)的發(fā)展，數(shù)值模擬憑借其經(jīng)濟性、高效性和準(zhǔn)確性，逐漸被學(xué)者應(yīng)用于侵徹問題的研究。Warren 等[10]建立了彈丸和石灰?guī)r靶體的有限元模型，對斜侵徹過程中彈道軌跡進行預(yù)測，與室內(nèi)實驗結(jié)果基本吻合。Tham[11]采用Lagrange、Euler-Lagrange 耦合和SPH-Lagrange 耦合3 種方法模擬彈體侵徹過程，預(yù)測混凝土靶體的最大侵徹深度以及徑向應(yīng)力-時間的關(guān)系。強洪夫等[12]進一步優(yōu)化了SPH（Smooth particle hydrodynamics）算法，用于花崗巖目標(biāo)的碰撞侵徹研究，給出了成坑深度與侵徹速度的關(guān)系式。在侵徹過程中，材料往往處于高應(yīng)力、大變形狀態(tài)，采用Lagrange 有限元方法的數(shù)值計算會面臨網(wǎng)格畸變、計算精度下降等問題[13]，即使引入合理的單元刪除準(zhǔn)則，也會帶來系統(tǒng)能量不守恒、結(jié)構(gòu)剛度降低等問題[14]；而SPH 方法屬于無網(wǎng)格方法，避免了強沖擊載荷下網(wǎng)格扭曲和單元失效，但是計算每個粒子的物理量比較費時，其計算效率遠(yuǎn)低于有限元法[15]，而且SPH 方法在固體力學(xué)的應(yīng)用中存在邊界拉伸不穩(wěn)定等問題[16]。即便有學(xué)者提出采用SPH-Lagrange 耦合的方法來提高計算效率和避免邊界拉伸的不穩(wěn)定性，但是這種方式導(dǎo)致系統(tǒng)總能量嚴(yán)重?fù)p耗[17]。

在上述研究的基礎(chǔ)之上，依據(jù)Li 等[18]的侵徹實驗，利用有限元分析軟件LS-DYNA，選取合適的網(wǎng)格尺寸和模型參數(shù)，建立彈體和花崗巖靶體的計算模型，采用Lagrange 算法、SPH-Lagrange 耦合算法以及SPH 算法分別進行數(shù)值模擬，基于計算效率、侵徹深度、速度衰減、靶體損傷、Mises 應(yīng)力分布等方面，對比分析3 種算法在彈體侵徹、貫穿巖石靶體數(shù)值模擬中的優(yōu)勢及不足。

1 模型與算法

1.1 計算模型

本研究采用Lagrange、SPH-Lagrange 及SPH 算法分別模擬彈體侵徹花崗巖的全過程，采用文獻(xiàn)[5]中的彈體侵深及靶體損傷，驗證數(shù)值模擬的有效性，實驗數(shù)據(jù)見表1。實驗中，采用二級輕氣炮，使子彈以1 200～1 600 m/s 的速度正侵徹靶體，子彈和靶體試樣見圖1[5,18]。其中，子彈直徑10.8 mm，長54.0 mm，彈頭形狀系數(shù)CRH（Caliber-radius-head）值為3.0，彈體材料采用30CrMnSiNi2A 高強度合金，密度7 850 kg/m3，初始質(zhì)量32.45 g，楊氏模量211 GPa；花崗巖靶體試樣長300 mm，寬300 mm，高800 mm，四周用外徑480 mm、厚10 mm 的鋼桶包裹，并用10 mm 厚鋼板封底，在鋼桶和試樣之間用混凝土填充，確保對試樣的約束，整體可視為半無限靶體。靶體材料選取山東省五蓮縣花崗巖，密度2 670 kg/m3，單軸抗壓強度約150 MPa，楊氏模量54.6 GPa，泊松比0.3，縱波速度4 200 m/s。

表1 侵徹實驗數(shù)據(jù)Table 1 Test results of projectile penetration

圖1 實驗試件Fig. 1 Experiment sample

數(shù)值模擬中，考慮彈體和目標(biāo)靶體的軸對稱性，在對稱面設(shè)置對稱邊界，建立1/4 模型，以提高計算效率，并在靶體模型外圍設(shè)置無反射邊界，模擬半無限靶體。為研究彈體貫穿靶體的情況，將靶體模型高度由800 mm 設(shè)置為100 mm（實驗中子彈速度為1 196 m/s 時的侵深為118.80 mm）。為更好區(qū)分兩類有限元模型，如圖2 所示，用于模擬實驗的數(shù)值模型稱為800 mm 模型，用于研究彈體貫穿靶體的數(shù)值模型稱為100 mm 模型。

圖2 計算模型Fig. 2 Model of projectile and granite target

初始時刻建立Lagrange 有限元模型，模型網(wǎng)絡(luò)劃分采用Solid 164 八節(jié)點六面體單元，利用侵蝕面面接觸算法模擬彈靶接觸，并添加單元刪除準(zhǔn)則。然后利用LS-PrePost 將單元轉(zhuǎn)化為SPH 粒子，用侵蝕點面接觸代替侵蝕面面接觸，不再使用單元刪除準(zhǔn)則。如圖2 所示，800 mm 模型中，由于單元數(shù)眾多，將其全部轉(zhuǎn)化為SPH 粒子后計算困難，所以只將彈靶作用近區(qū)內(nèi)的單元轉(zhuǎn)化為SPH 粒子，建立SPH-Lagrange耦合模型，SPH 粒子與Lagrange 網(wǎng)格之間設(shè)置點面綁定接觸；100 mm 模型中，單元轉(zhuǎn)化為SPH 粒子，分為兩種形式，一是建立SPH-Lagrange 耦合模型，二是將靶體單元全部轉(zhuǎn)化為SPH 粒子，建立SPH 模型。

1.2 模型選擇

本研究中，子彈采用LS-DYNA 數(shù)據(jù)庫中*MAT_PLASTIC_ KINEMATIC 模型，該模型屬于考慮應(yīng)變率效應(yīng)的雙線性模型，模型屈服后會沿線性硬化，可描述各向同性硬化和隨動硬化塑性，能夠很好地模擬金屬等彈塑性材料的力學(xué)響應(yīng)[19]。詳細(xì)參數(shù)[20]見表2，其中 ρp為質(zhì)量密度，VP 表示考慮應(yīng)變率效應(yīng)的計算方法。當(dāng)VP = 0 時，通過縮放屈服應(yīng)力進行計算；當(dāng)VP = 1時，調(diào)用黏塑性公式執(zhí)行計算。

表2 子彈材料參數(shù)Table 2 Properties of test projectile

模型本構(gòu)方程為

根據(jù) β值的不同，可以用來描述不同的硬化模型，如圖3(a)所示： β=0 時，隨動硬化，屈服面大小不變，沿塑性應(yīng)變方向移動； β=1 時，各向同性硬化，屈服面位置不變，大小隨應(yīng)變而變化； 0 ＜β ＜1 時，混合硬化。

圖3 模型描述Fig. 3 Model description

花崗巖靶體選用LS-DYNA 數(shù)據(jù)庫中的*MAT_RHT 模型。該模型由德國Ernst-Mach 研究所（EMI）的 Riedel、Thoma 和Hiermaier 于1998 年首次提出[19]，經(jīng)過近20 年的發(fā)展，被廣泛應(yīng)用于爆炸沖擊等領(lǐng)域的數(shù)值模擬中，能夠較好地描述混凝土、巖石等脆性材料受侵徹爆炸作用時的動態(tài)力學(xué)行為。RHT 模型充分考慮了彈性極限面、失效面以及損傷軟化段殘余強度，用于研究混凝土在動載作用下的初始屈服強度、失效強度和殘余強度的變化[21]，還可以較準(zhǔn)確地描述混凝土、巖石等脆性材料應(yīng)力-應(yīng)變曲線彈性、線性強化以及損傷軟化3 個階段，如圖3(b)所示，其中， σˉ 、 εˉ表示等效應(yīng)力、等效應(yīng)變，σfail、σelastic和σresidual分別為失效應(yīng)力、彈性極限應(yīng)力和殘余應(yīng)力。同時，RHT 模型采用了Herrmann 提出的p-α狀態(tài)方程，用于描述在高速沖擊等問題中當(dāng)材料處于流動狀態(tài)時靜水壓力、材料密度和內(nèi)能之間的關(guān)系[22]。

對于材料壓密后的狀態(tài)方程，用多項式函數(shù)表示

在測得花崗巖基本物理力學(xué)參數(shù)的基礎(chǔ)上，通過理論計算和參數(shù)研究[23]確定用于模擬花崗巖的RHT 模型參數(shù)，具體參數(shù)見表3。

表3 RHT 模型參數(shù)Table 3 Parameters of RHT material model

1.3 SPH 算法

SPH 算法是一種流體力學(xué)模擬的無網(wǎng)格光滑粒子算法，以插值理論為基礎(chǔ)，通過任意分布的光滑粒子來計算各種邊界條件下的物理力學(xué)方程。在算法中，并不是單獨的計算粒子，而是通過核函數(shù)來逼近區(qū)域內(nèi)任意粒子周圍的場變量?？紤]材料的彈塑性效應(yīng)時，核函數(shù)W和插值公式可分別表示為[13]

2 數(shù)值驗證

2.1 網(wǎng)格收斂性分析

選擇合適的網(wǎng)格尺寸是精確模擬的關(guān)鍵，為確定本研究中花崗巖模型的網(wǎng)格大小，首先建立正方形網(wǎng)格，并采用無量綱比例參數(shù) λ，其定義為彈體半徑與網(wǎng)格長度之比。如表4 所示，當(dāng) λ為 1、2、3、4、5、6 時，構(gòu)造相應(yīng)的數(shù)值模型，并進行彈體侵徹深度計算。在數(shù)值模擬中，子彈選用*MAT_PLASTIC_KINEMATIC 模型，彈體材料的屈服強度高，侵徹結(jié)束后，彈體只發(fā)生輕微形變，彈徑未發(fā)生變化，可以將其視為具有剛性性能的高屈服強度彈體，因此，在固定彈徑條件下進行靶體網(wǎng)格收斂性分析時，子彈采用該雙線性模型是合理可行的。圖4 展示了彈體速度為1 196 m/s 時，不同比例參數(shù) λ對應(yīng)的侵徹深度-時間曲線。隨著λ 的增加，即隨著網(wǎng)格尺寸變小，侵徹深度逐漸趨于收斂， λ為 4、5、6 時，侵徹深度基本一致，且與實驗結(jié)果較接近，滿足精度要求。因此，為減少計算資源占用，提高計算效率，后續(xù)仿真模型采用 λ = 4 的比例參數(shù)，即網(wǎng)格尺寸為1.35 mm。

圖4 不同網(wǎng)格尺寸下侵深-時間關(guān)系曲線Fig. 4 Time history of penetration depth with different mesh sizes

表4 網(wǎng)格尺寸與侵徹深度Table 4 Mesh size and penetration depth

2.2 數(shù)值驗證

選取表1 中的3 組實驗（No.1、No.2 和No.4）進行數(shù)值模擬，子彈選用*MAT_PLASTIC_KINEMATIC 模型，花崗巖靶體選用*MAT_RHT 模型，網(wǎng)格比例參數(shù)取 λ = 4，模擬獲得的侵深-時間曲線如圖5 所示。當(dāng)子彈速度為1 196、1 400、1 630 m/s 時，對應(yīng)的侵深分別為119.07、149.17、171.80 mm，與表1 中實驗數(shù)據(jù)較為接近。然而隨著侵徹速度不斷增加，模擬中侵深與實驗之間的差距逐漸變大，主要原因是隨著侵徹速度增加，實驗中子彈的磨損不斷增大，而模擬中子彈形變不大，從而導(dǎo)致模擬中侵徹深度略高于實驗。圖6 顯示了侵徹速度為1 400 m/s 時，模擬中靶體損傷情況及實驗中靶體宏觀破壞形態(tài)，其中數(shù)值模擬通過云圖的形式顯示靶體的破壞損傷情況，0 表示巖石未損傷，1 表示巖石已完全破壞。

圖5 不同侵徹速度下侵深-時間關(guān)系曲線Fig. 5 Time history of penetration depth with different velocities

圖6 侵徹速度為1 426 m/s 時靶體的損傷Fig. 6 Target damage at an impact velocity of 1 426 m/s

通過實驗與數(shù)值計算的結(jié)果對比，數(shù)值計算獲得的侵徹深度及靶體損傷與實驗結(jié)果吻合度高，模擬重現(xiàn)了花崗巖靶體在彈體高速侵徹下的損傷與破壞，因此建立的計算模型適用于彈體高速侵徹花崗巖靶體的相關(guān)研究，尤其是在預(yù)測彈體侵深和靶體損傷方面。

3 數(shù)值計算及結(jié)果分析

基于已驗證的本構(gòu)模型和網(wǎng)格比例參數(shù) λ，在子彈初速度為1 196 m/s 的條件下，分別對800 mm 模型和100 mm 模型進行數(shù)值計算。其中，800 mm 模型分為兩類：一是有限元模型，命名FEM-1（Finite element model-1），單元數(shù)3 647 448；二是耦合模型，命名SPH-FEM-1（Smooth particle hydrodynamics-finite element model-1），單元數(shù)3 341 923，SPH 粒子數(shù)305 525。100 mm 模型分為3 類：一是有限元模型，下稱FEM-2，單元數(shù)690 408；二是耦合模型，下稱SPH-FEM-2，單元數(shù)521 008，SPH 粒子數(shù)169 400；三是SPH 模型，下稱SPH-2，SPH 粒子數(shù)690 408。

3.1 800 mm 模型

圖7 給出了采用Lagrange 算法和SPH-Lagrange 耦合算法數(shù)值模擬得到的侵徹路徑和靶體損傷的對比。如圖7 所示，Lagrange 算法導(dǎo)致彈體侵徹路徑周邊的單元發(fā)生畸變，耦合算法則無單元畸變，能夠很好地模擬子彈擊打靶體時粒子濺射飛散的過程，這是因為彈體高速侵徹時Lagrange 網(wǎng)格產(chǎn)生大變形而SPH 粒子可以任意流動。同時，與Lagrange 算法相比，耦合算法中侵徹會形成較大前坑區(qū)，而且侵徹路徑周邊嚴(yán)重?fù)p傷區(qū)域（0.8～1.0）直徑更大，與實驗結(jié)果吻合度更高，這是由于Lagrange 算法中需要添加刪除準(zhǔn)則，彈體打擊產(chǎn)生的應(yīng)力波在向四周傳遞過程中由于單元刪除導(dǎo)致更多能量損耗，因此前坑區(qū)和嚴(yán)重?fù)p傷區(qū)域比較窄。但是，耦合算法中，由于應(yīng)力波在不同算法模型中傳播速度不同以及點面綁定接觸算法在接觸力傳遞方面有一定缺陷，導(dǎo)致在SPH 粒子和Lagrange 網(wǎng)格接觸面的兩側(cè)，損傷變化非線性，存在明顯的損傷差異。

圖7 靶體損傷云圖Fig. 7 Target damage of 800 mm model

如圖8 所示，在SPH 粒子與單元接觸面左右13.5 mm 的范圍內(nèi)，獲取靶體自由面和對稱面交界面上單元或SPH 粒子的PPV（Peak particle velocity），編號1～20。從圖8 中看出，相較于Lagrange 算法，耦合算法中粒子區(qū)域的PPV 較低，而且在粒子與單元接觸面的PPV 驟增，這與彈體和SPH 粒子及粒子和單元之間的點面接觸有關(guān)。但是，在有限單元區(qū)域，兩種算法中的PPV 逐漸趨于一致。由此，基于PPV 準(zhǔn)則[24]，應(yīng)用耦合算法研究花崗巖靶體的侵徹?fù)p傷時，要注意中心區(qū)域的應(yīng)力及損傷可能會低于預(yù)期，粒子和單元的過渡區(qū)域損傷會高于預(yù)期，而一定范圍外花崗巖靶體的損傷和穩(wěn)定性可以基于PPV 準(zhǔn)則進行較準(zhǔn)確的評估。

圖8 質(zhì)點峰值振動速度Fig. 8 PPV of target points

圖9 展示了兩種算法得到的靶體的Mises 應(yīng)力分布情況。圖10 為兩種算法得到的侵徹速度-時間關(guān)系曲線。由圖9 可以看出，與Lagrange 算法中靶體Mises 應(yīng)力環(huán)形分布相比，耦合算法中，中間正方形區(qū)域為SPH 粒子，在粒子與單元接觸表面的4 個角出現(xiàn)應(yīng)力滯后，接觸面比尖角處應(yīng)力大，主要是由粒子與單元之間的接觸算法導(dǎo)致的。由圖10 可以看出，Lagrange 算法中侵深為119.07 mm，耦合算法中侵深為114.79 mm，均與實驗結(jié)果118.80 mm接近，考慮到數(shù)值模擬中彈體無磨損，因此Lagrange算法更符合實際情況。但是，從侵徹速度-時間關(guān)系曲線中發(fā)現(xiàn)，與耦合算法相比，Lagrange 算法得到的速度衰減不平滑，這是由算法中添加的刪除準(zhǔn)則造成的。同時，兩種算法中侵徹速度衰減至零后均出現(xiàn)速度反向增加的情況，主要原因是彈體采用*MAT_PLASTIC_ KINEMATIC 模型描述，侵徹過程中產(chǎn)生形變，儲存部分變形能，在侵徹結(jié)束后，能量釋放，導(dǎo)致彈體出現(xiàn)反彈。

圖9 0.04 ms 時Mises 應(yīng)力分布Fig. 9 Von Mises stress at 0.04 ms

圖10 侵深-時間和速度-時間關(guān)系曲線Fig. 10 Time history of penetration depth and velocity

計算效率方面，設(shè)置同樣的計算總時間和時間步長，應(yīng)用Lagrange 算法數(shù)值模擬需要3.0 h，耦合算法需要17.1 h。綜合考慮靶體損傷、彈體速度和侵深、Mises 應(yīng)力分布以及計算時間等方面，800 mm 模型的計算中，如果只考慮侵徹路徑近區(qū)的靶體損傷和粒子飛濺效果，耦合算法更合適；如果考慮計算效率、侵徹深度、應(yīng)力波傳播等，推薦Lagrange 算法。

3.2 100 mm 模型

圖11 為Lagrange 算法、SPH-Lagrange 耦合算法和SPH 算法數(shù)值模擬得到的侵徹路徑和靶體損傷的對比。由圖11 可知，Lagrange 算法導(dǎo)致彈體侵徹路徑周邊的單元發(fā)生畸變，其他兩種算法則無單元畸變，能夠很好地模擬子彈擊打靶體時粒子濺射飛散過程，這得益于SPH 粒子的任意流動性。在耦合算法和SPH 算法中，侵徹會形成較大的前坑區(qū)和后坑區(qū)，而且侵徹路徑周邊嚴(yán)重?fù)p傷（0.8～1.0）范圍更廣，與彈體侵徹貫穿靶體的工況相吻合；同時Lagrange 算法只有較小的前坑區(qū)，不存在后坑區(qū)，這主要與添加的刪除準(zhǔn)則有關(guān)。在耦合算法中，由于應(yīng)力波在不同算法模型中傳播速度不同以及點面綁定接觸算法在接觸力傳遞方面有一定缺陷，導(dǎo)致在SPH 粒子和Lagrange 網(wǎng)格接觸面的兩側(cè)，損傷變化非線形，存在明顯的損傷差異，而其他兩種算法則不存在這種差異。

圖11 靶體損傷云圖Fig. 11 Target damage of 100 mm model

圖12 給出了3 種算法中靶體的Mises 應(yīng)力分布情況。圖13 給出了3 種算法中彈體侵徹速度-時間關(guān)系曲線。由圖12 可知，在耦合算法中，中間正方形區(qū)域為SPH 粒子，在粒子與單元接觸表面的4 個角會出現(xiàn)應(yīng)力滯后，接觸面比尖角處應(yīng)力大，而在其他兩種算法中靶體的Mises 應(yīng)力與實際情況一致，呈環(huán)形分布，這種應(yīng)力分布的差異主要是由粒子與單元之間的接觸算法造成的。由圖13 可以看出，在耦合算法和SPH 算法中，兩種彈體侵徹速度-時間曲線平滑且剩余速度一致。但是在Lagrange 算法中，由于添加的單元刪除準(zhǔn)則，導(dǎo)致侵徹速度衰減存在波動且剩余速度較大。

圖12 0.04 ms 時Mises 應(yīng)力分布Fig. 12 Von Mises stress at 0.04 ms

圖13 采用不同算法得到的侵徹速度-時間關(guān)系曲線Fig. 13 Time history of velocity with different algorithm

計算效率方面，設(shè)置同樣的計算總時間和時間步長，采用Lagrange 算法進行數(shù)值模擬需要0.65 h，耦合算法需要3.60 h，SPH 算法需要8.30 h。綜合考慮靶體損傷、彈體速度和侵深、Mises 應(yīng)力分布以及計算時間等方面，100 mm 模型的數(shù)值模擬中，如果只考慮彈體侵徹貫穿靶體形成的前/后坑區(qū)、粒子飛濺效果以及子彈的剩余速度，耦合算法更合適；如果只考慮時間效率，則推薦Lagrange 算法；如果考慮數(shù)值模擬的準(zhǔn)確性、應(yīng)力波傳播等，則SPH 算法更優(yōu)。

4 結(jié) 論

基于Lagrange 及 SPH 算法，對彈體高速侵徹巖石靶體進行了數(shù)值模擬，對比分析了Lagrange算法、SPH-Lagrange 耦合算法以及SPH 算法的優(yōu)劣性。

（1）在彈體侵徹巖石靶體的數(shù)值模擬中，Lagrange 算法的計算效率高，計算結(jié)果較準(zhǔn)確，通過設(shè)置無反射邊界可以模擬彈體侵徹半無限靶體，但是存在無法展示靶體濺射、前坑區(qū)小、無后坑區(qū)等問題，利用RHT 模型顯示的損傷云圖，則可以獲取較準(zhǔn)確的前坑區(qū)和后坑區(qū)范圍，進而間接解決Lagrange 算法存在前坑區(qū)小、無后坑區(qū)的問題。

（2）SPH 算法可以清晰顯示侵徹時靶體粒子的飛散，計算結(jié)果較準(zhǔn)確，但是計算效率低，粒子過多時還可能造成計算停止。

（3）SPH-Lagrange 耦合算法兼具SPH 和Lagrange 算法的優(yōu)勢，但是由于粒子與單元之間接觸的問題，會產(chǎn)生應(yīng)力滯后和應(yīng)力波不穩(wěn)定傳播的現(xiàn)象。在應(yīng)用耦合算法研究侵徹爆炸問題時，關(guān)于侵徹爆炸中心區(qū)域的巖體破壞判定，需要考慮應(yīng)力滯后、應(yīng)力波不穩(wěn)定傳播等的影響。關(guān)于一定范圍外地下結(jié)構(gòu)的穩(wěn)定性評估，可基于PPV 準(zhǔn)則進行較準(zhǔn)確的研究。

（4）在單元數(shù)較少的小模型中，采用SPH 算法可以更準(zhǔn)確地研究侵徹問題。在單元數(shù)較多的大模型中，Lagrange 算法是比較合理的選擇，耦合算法則需要進一步研究。