爆炸沖擊下珊瑚砂動態本構模型*

董 凱，任輝啟,，阮文俊，黃 魁，步鵬飛

（1. 南京理工大學能源與動力工程學院，江蘇 南京 210094；2. 中國人民解放軍軍事科學院國防工程研究院，河南 洛陽 471023）

隨著我國海上絲綢之路重大戰略的逐步實施，島礁防護工程的建設需求被密切關注。珊瑚砂廣泛分布于我國南海島礁和瀉湖中，它作為覆土層、填充材料等應用于島礁工程中。在島礁工程面臨突發的動力災變時，珊瑚砂往往先受到來襲彈丸的侵徹、爆炸作用，而成為重要的研究對象[1-2]。

近些年，雖然對于珊瑚砂的本構關系開展了一些研究，但模型主要應用于沉樁、土工建設以及流變分析等低應變率的工程中[3-5]，而對中高應變率加載下珊瑚砂的研究，多見于力學特性分析、顆粒破碎和耗能的研究中，且闡述了顆粒破碎與壓縮特性之間的聯系[6-9]。這些研究成果對推動島礁建設具有重要的作用，然而對于解決復雜的爆炸、侵徹等非線性問題，需要借助于大型商用計算軟件如LS-DYNA、AUTODYN 等進行分析，如徐學勇[10]使用有效應力彈塑性本構模型對飽和鈣質砂在爆炸下的動力響應做了分析。但目前對干燥珊瑚砂動態本構模型的研究十分匱乏，面對大量的實際工程計算需求，急需確定出可嵌入有限元程序中并適用于侵徹、爆炸計算的模型和參數。

砂土的力學參數隨地域分布的不同呈現多樣化，且模型在應用時往往局限于特定的條件，這使沖擊下砂土本構關系具有多種表述形式。目前，廣泛應用的模型大多采用以下兩種框架構建：（1）基于沖擊動力學實驗（如SHPB 實驗、泰勒桿實驗、飛片沖擊等）獲得材料的高應變率壓縮方程，結合描述偏應力的彈塑性關系構成流體彈塑性模型[11-12]；（2）使用靜力學下的基本模型構架，通過增加率相關的參數考慮應變率效應，將模型適用于高應變率加載時的計算[13-14]。兩種框架構建的本構模型在陸源砂土中均得到了廣泛的應用，而在多孔易破碎的珊瑚砂中模型的適用性尚未確定，因此珊瑚砂動力學模型框架適用性的研究應作為模型研究的首要工作。

本文中，根據靜態、動態力學實驗結果，分別基于流體彈塑性模型和Perzyna 黏塑性帽蓋模型確定珊瑚砂的模型參數，借助LS-DYNA 有限元程序實現模型在珊瑚砂中的應用，通過對彈丸侵徹珊瑚砂過程以及平面波加載下珊瑚砂中應力波衰減過程的數值計算驗證模型的有效性，同時對不同壓實密度下珊瑚砂的侵徹規律、爆炸應力波衰減規律進行計算，擬對武器戰斗部在珊瑚砂中的侵徹、爆炸的研究提供重要參考。

1 本構模型和參數

1.1 基于剪切屈服和體積壓縮規律的本構模型

珊瑚砂為脆性松散孔隙巖土材料，在沖擊、爆炸等強動載下具有流體特性，并伴隨著大量的顆粒破碎。高應變率（102s?1以上）下的體積壓縮規律能準確地體現在沖擊受壓過程中材料的宏觀壓密現象，壓實過程在物理上表現為體積不可逆壓縮行為，與靜水壓有關，因此可用不可逆加、卸載路徑下的壓力-體應變函數即物態方程確定。剪切屈服特性采用正應力與偏應力之間的二次函數描述。由于物態方程是在高應變率下測得的，可適用于爆炸沖擊計算，模型中不包含直接體現應變率的參數。LS-DYNA 中的5#材料模型（*MAT_SOIL_AND_FOAM）能很好表達該類模型[15]，其本質屬于流體彈塑性模型，在爆炸、侵徹計算中得到廣泛引用[16-18]，然而模型參數的確定方法未見明確報道，針對珊瑚砂的參數也少有研究。

1.1.1 屈服面函數

假設剪切強度準則可由3 個參數a0、a1、a2構成的二次函數F（J2, p）表達，屈服函數上的偏應力第二不變量J2與平均靜水壓力p 的關系定義為：

土力學中，通常認為無黏性砂的黏聚力為零。而大量實驗結果表明，珊瑚砂存在一定的黏聚力[19-20]，主要源于不規則棱角的砂顆粒間的相互咬合，剪切時需要破壞此咬合結構才能達到破壞。Fasanella 等[17]給出a0=a1=0 來表征黏聚力為零，這與實際是不相符的。Wright[21]給出了確定屈服面參數的方法，但公式中存在一些錯誤，這里修正如下。

在剪切屈服面上，珊瑚砂符合Mohr-Coulomb準則，強度包絡線如圖1 所示。在破壞面上有關系式：

圖1 確定屈服參數的摩爾圓Fig.1 Mohr’s circle geometry used to determine yield surface parameters

式中：σx為測得的軸向應力，在試樣均勻化加載的條件下與正應力相等，在屈服時 σf=σx；σr為側限壓力；ξ為側壓力系數，定義為ξ =σr/σx，并通過SHPB 實驗確定ξ =0.495。為便于公式推導，這里近似取p=2σf/3，可依據式(4)得到：

并可得到：

則屈服面參數可表示為：

文獻[16]中的土壤模型參數被廣泛應用于數值計算中，根據其中參數a0=3 400 kPa2，a2=0.30，利用上述方法計算得到擬合參數a1=63.96 kPa，與文獻[16]中a1=70.33 kPa 進行對比，誤差約為9.06%。可見，該方法可對屈服面參數進行確定。

根據直剪實驗結果，珊瑚砂的黏聚力C=7.24 kPa，內摩擦因數tan φ=0.462，由式(8)可得到：a0=84.77 kPa2，a1=16.23 kPa，a2=0.777。由于巖土材料對溫度的敏感性較低，模型沒有考慮溫度項的影響。當靜水拉應力超過了拉伸截斷閾值時，則設置拉應力為截斷閾值且偏應力張力為零。由于松散的珊瑚砂無法承受拉伸狀態，可取拉伸截斷閾值為零。

1.1.2 沖擊時平均壓力-體應變方程

在高荷載作用下，多孔珊瑚砂被壓縮密實，引起體積模量的增大，也稱為材料的壓硬性。由于珊瑚砂的顆粒破碎特性，加載初期卸載后基本不發生回彈，但在壓實后材料凝聚為塊狀，卸載時可認為彈性卸載。因此，使用分段物態方程分別描述加、卸載行為：

式中：p 為平均壓力，?V為體積應變。

高應變率壓縮曲線對于模型的確定至關重要，然而通過SHPB 實驗獲得的曲線應變僅達到約0.12，對應的應力也往往不超過20 MPa，無法描述高孔隙率的珊瑚砂動態壓實特性。通過比較靜態與動態的應力應變曲線特征，可發現：（1）動態屈服應變與靜態近似相等；（2）珊瑚砂沒有明顯的應變強化拐點εh，表明在壓實過程中顆粒破碎比較平穩，變形機制未發生根本變化；（3）珊瑚砂在屈服后強化因數σd/σs隨應變的增加表現出線性關系[23]。

基于上述結論，假設線性關系持續至壓實狀態，根據文獻[23] 中相對密實度Dr為0.90（密度為1.260 g/cm3）的靜態、高應變率下的一維壓縮應力應變曲線，以準靜態應力應變為基準，考慮強化因數的影響（σd/σs=1.663），結合測得的圍壓將準靜態一維應力應變曲線換算成高應變率壓力-體應變曲線，如圖2所示。文祝等[22]通過預壓方法確定了珊瑚砂的壓縮方程，為p-?V關系的擬合提供了有效的方法，與本文確定的壓縮曲線相比，整體趨勢保持一致，但由于級配的不同導致曲線存在一定差別。兩種方法中最大應力均達到100 MPa 以上，使用應變率強化因數擬合壓實段物態方程的方法，有如下優點：（1）能夠準確描述動態下初始階段的屈服特性；（2）避免了預壓階段未計及的應變率效應問題。

根據實驗結果[23]確定了3 種不同密度的珊瑚砂的p-?V曲線，如圖3 所示。隨著密度的增大，在相同體應變下受到的壓力更大，為了準確描述屈服和硬化特性，宜采用分段函數表達，擬合的物態方程為：

式中：a、m、b、n 均為密實度Dr的線性函數，模型的適用性因參數的線性特征而增強。它們可分別表示為：

此外，模型中加載體積模量K 可由SHPB 實驗確定，取彈性條件下的初始值，砂土材料泊松比μ通常為0.3，可確定剪切模量G=3K(1-μ)/[2(1+μ)]，卸載體積模量Ku可根據壓實段的模量確定[21]。綜合以上，得到基于5#模型的3 種密度珊瑚砂參數，見表1～3。

圖2 平均壓力-體應變的擬合曲線Fig.2 Average pressure-volumetric strain fitting curves

圖3 在不同相對密實度下平均壓力與體應變的關系Fig.3 Average pressure-volumetric strain curves under different compactness levels

表1 當 D r=0.30時珊瑚砂的5#材料模型參數Table 1 Parameters of 5# constitutive model for coral sand when Dr=0.30

表2 當 D r=0.60時珊瑚砂的5#材料模型參數Table 2 Parameters of 5# constitutive model for coral sand when Dr=0.60

表3 當 D r=0.90時珊瑚砂的5#材料模型參數Table 3 Parameters of 5# constitutive model for coral sand when Dr=0.90

1.2 基于Perzyna 形式的黏塑性帽蓋模型

定義帽蓋模型的屈服函數f 為應力張量第一不變量I1和偏應力第二不變量J2的關系，屈服面可分為3 部分（見圖4），以下簡單介紹。

圖4 黏塑性帽蓋模型的屈服面Fig.4 Yield surface for viscoplastic cap model

（1）當 I1≤?T時，定義為拉伸失效面部分，可表示為：

式中：?T 為材料拉伸截止閾值。

（2）當 ? T＜I1＜L(k)時，定義為強度失效面區，表示為：

破裂面為非硬化的，因此采用不考慮硬化的改進Drucker-Prager 形式，為：

式中：α、β、γ、θ 為材料特性相關的參數。

（3）當 I1≥L(k)時，定義為帽蓋面部分，表示為：

觀察患者呼吸機相關性肺炎的發生率,若通氣治療48h后符合:①X線胸片示新的或進行性肺浸潤;②發熱;③外周血白細胞計數>20.0×109/L或C反應蛋白>8mg/L;④氣道分泌物細菌培養陽性。基礎條件為X線胸片所示改變,若另外3條中2條符合,即可診斷患者患有呼吸機相關性肺炎[4]。

帽蓋面定義為半橢圓形，為：

式中：X(k)位于帽蓋面與橫軸I1軸的交點，為：

L(k)為帽蓋在開始位置時的靜水壓值，并符合：

由此，帽蓋面可以表示為：

式中：W 為最大塑性體應變參數，D 為體積變化率參數，X0與帽蓋面在橫軸的初始位置有關。

Perzyna 黏塑性帽蓋模型需確定14 個參數：W、D、R、X0為帽蓋面函數的參數，α、β、γ、θ 為強度失效面函數中的參數，T 為拉伸截止面中的參數，η、f0、N 為描述黏塑性流動法則的參數，還有，描述珊瑚砂的彈性響應階段的加載體積模量K 和剪切模量G。

丁育青[14]給出了該模型參數的詳細確定方法，并應用于非飽和黏土中的爆炸計算，結合珊瑚砂的三軸圍壓實驗[20]、準靜態壓縮實驗和SHPB 實驗結果[23]擬合得到了密實度0.30 珊瑚砂的Perzyna 模型參數，見表4。由于Perzyna 黏塑性模型參數較多，需依靠大量的力學實驗確定，這里基于力學實驗結果僅擬合得到了一組參數用于分析計算，模型需使用LS-DYNA 中的二次開發模塊自定義本構，具體算法可參考文獻[13-14, 24]。

表4 當 D r=0.30時珊瑚砂Perzyna 黏塑性帽蓋模型參數Table 4 Perzyna viscoplastic cap model parameters of coral sand when Dr=0.30

2 數值模擬結果和分析

2.1 算例和結果

苗偉偉等[25]針對珊瑚砂在彈丸侵徹下的規律進行了一系列的實驗研究，獲得了初速為351～972 m/s的彈丸侵徹深度，同時使用理論方法進行了侵徹深度的計算并獲得了理想的結果。目前，對于砂土的侵徹模型大都基于空腔膨脹理論進行研究[26]，為了確定所建立模型的適用性，使用LS-DYNA 對實驗的侵徹過程進行數值計算。實驗中彈丸形狀和尺寸如圖5 所示，彈丸材料為35CrMnSiA 合金鋼，質量為80 g，在侵徹計算時采用剛體模型（*MAT_RIGID），密度ρ=7.85 g/cm3，彈性模量E=210 GPa，泊松比μ=0.29。珊瑚砂采用5#模型（Dr=0.90，ρ=1.260 g/cm3）。模型中的密度參數與侵徹實驗一致，彈丸與珊瑚砂采用侵蝕接觸，采用1/2 對稱三維模型，計算靶體為厚10 cm、寬30 cm、長160 cm 的立方體，網格劃分如圖6 所示，其中彈丸網格數量為4 194，靶體網格數量為600 萬。

圖5 彈丸形狀和尺寸[25]Fig.5 Projectile geometry[25]

圖6 侵徹計算模型網格劃分（靶體為部分顯示）Fig.6 Finite element mesh of calculated model (target is partially displayed)

計算得到的彈丸在不同入射速度下的最終侵徹深度如圖7 所示，可見基于5#模型的計算結果與實驗結果良好吻合，表明模型適用于剛性彈丸侵徹珊瑚砂的數值模擬。最終侵徹深度隨入射速度的變化呈非線性關系，在不同初速侵徹時，彈丸速度隨深度的下降規律如圖8 所示。在侵徹初期由于過載非常大，隨著入射速度的增大，彈丸速度下降非常明顯，呈指數形式降低；隨著速度的降低，過載明顯減小。該結論與Omidvar 等[27]通過實驗測得的土中彈丸侵徹過載的變化規律一致，表明砂土在剛性彈丸高速侵徹時呈流體響應。

圖7 最終侵徹深度與入射速度的關系Fig.7 Final penetration depth versus initial velocity

圖8 不同入射速度時速度與深度的關系Fig.8 Velocity versus penetration depth at different initial velocities

在文獻[25]的實驗工況C27 中，彈丸以入射速度972 m/s 侵徹珊瑚砂，在侵徹初期彈丸速度遠高于珊瑚砂波速（通常低于300 m/s）。用5#模型計算得到的珊瑚砂高壓區呈錐形分布，具有典型的流體特征，如圖9 所示。隨著速度的下降，高壓損傷區的邊界輪廓逐漸由錐形向圓形過渡演化，侵徹過程中彈丸頭部始終處于高壓力作用狀態，因此在頭部出現了明顯的磨蝕損傷，這與文獻[25]的實驗相符。在侵徹計算過程中，彈體發生偏轉現象，文獻[25]中雖未對該現象進行闡述，但彈道偏轉情況在侵徹過程中難以避免，尤其在顆粒材料的侵徹中更加普遍[28]。通過以上結論與計算結果表明，基于5#模型對于彈丸侵徹珊瑚砂的數值計算是有效的。

圖9 珊瑚砂在不同時刻的壓力場和彈丸產生的磨蝕區Fig.9 Pressure fields of coral sand at diffident times and scratch area of projectile

為了探索不同壓實密度對珊瑚砂抗侵徹規律的影響，使用了1.1 節中的模型計算了彈丸對不同密實度珊瑚砂的侵徹，選擇彈丸速度510 m/s 的情形進行對比，獲得的速度隨深度的變化關系如圖10 所示。對比密實度0.30 的珊瑚砂，密實度0.90 的密度提高了約6.96%，侵徹深度降低了約5.03%。結果表明，提高壓實度雖然可以降低侵徹深度，但3 種壓實密度下侵徹深度相差很少。另外，本文的珊瑚砂為不良級配，其最大、最小干密度間的差很小，因此對侵徹深度的影響也很小。

圖10 不同壓實密度時速度與深度的關系Fig.10 Velocity versus penetration depth under different compactness levels

2.2 珊瑚砂中爆炸應力波衰減規律

武器戰斗部在砂土中爆炸時，產生的爆炸波沿砂土向周邊傳播，對地下結構造成破壞，這是防護工程中最關注的問題。砂土中應力波的衰減規律是決定傳播至結構表面處壓力的重要因素，而目前對珊瑚砂中的地沖擊效應研究較少。

趙章泳等[29]通過實驗確定了平面波加載下的珊瑚砂中應力波衰減規律，使用第1 節中確定的模型參數對該爆炸過程進行了數值計算，得到的應力波峰值衰減曲線如圖11 所示。可見，5#模型計算得到的衰減曲線與實驗較好吻合，而Perzyna 模型則略顯剛硬，表現為應力波衰減緩慢。

兩種模型計算得到的球形裝藥觸地爆時爆炸波衰減的一般規律，如圖12 所示，其中監測點1 位于比例距離0.2 m/kg1/3處，監測點2 位于比例距離0.3 m/kg1/3處。可以看出，兩種模型在相同時刻達到監測點，表明計算得到的應力波速度一致。5#模型計算得到的應力波在衰減過程中表現出上升沿漸緩、長歷時、低幅值的現象，于瀟等[30]、Yu 等[31]通過SHPB 實驗測量珊瑚砂中應力波衰減規律，也發現了該現象，認為顆粒破碎導致高頻分量被過濾，Perzyna 黏塑性模型無法描述因顆粒破碎產生的波前觀測弛豫現象，因此在易破碎的珊瑚砂中的計算存在較大誤差。而5#模型的壓縮曲線基于SHPB 實驗測得，有關破碎引起的高壓縮特性已經耦合于方程內部，能夠表達珊瑚砂的宏觀力學本質；同時，參考王禮立等[32]對塑性本構關系和流動型本構關系的研究可以推斷，基于流動率確定的本構模型更符合干燥珊瑚砂在爆炸、侵徹下的變形規律。

圖11 峰值壓力的衰減Fig.11 Calculated and experiment results of peak pressure attenuation

圖12 兩種模型壓力波Fig.12 Pressure waves calculated by two models

珊瑚砂堆積密度對應力波衰減的影響同樣具有舉足輕重的作用，然而由于模型參數的局限，對壓實密度這個參量影響爆炸波傳播規律的研究較少。使用1.1 節中的5#模型參數，對3 種密度的珊瑚砂中的爆炸波衰減規律進行數值計算，采用觸地爆炸模型（見圖13）作為研究對象，炸藥為200 g 的塊狀TNT，尺寸為10 cm×5 cm×2.5 cm，模型采用1/4 對稱結構，在對稱面施加對稱邊界，其余采用非反射邊界條件，空氣與珊瑚砂采用Euler 網格，均采用ALE 多物質單元，材料可在網格內流動，TNT 與空氣的材料模型及相關算法可參考文獻[33-34]。

圖13 計算模型Fig.13 Numerical model

塊狀炸藥爆炸時產生近似球形波，介質表面爆炸時處于炸藥正下方的壓力最大，通過計算得到壓力波峰值壓力隨比例距離的衰減關系如圖14 所示。可見，松散的珊瑚砂具有更好的消波能力，相比密實度0.90，密實度0.30 的珊瑚砂應力峰值最高可降低26.1%。因此，工程應用中，在承載力允許的條件下，使用松散珊瑚砂更利于爆炸波的衰減。

圖14 不同相對密度時爆炸峰值壓力與比例距離的關系Fig.14 Peak pressure versus scaled distance under different compactness levels

綜上所述，珊瑚砂的壓實密度對于彈丸侵徹的應力波衰減存在一定的影響，為提高承載力和抗侵徹能力，往往需要提高壓實密度；而面對爆炸波的沖擊作用，又需要降低密實度以增強消波能力。因此，不論是評價武器對島礁的毀傷效應還是島礁防護工程的設計，珊瑚砂的動態本構模型研究都十分重要。

本文中研究的珊瑚砂為原狀砂，在實驗時保留了原生顆粒級配的完整性。工程中有時將砂土作夯實處理，密度超過原生砂的最大干密度但同時會引起顆粒破碎，此時，需對本構模型另做改進。

3 結 論

研究的兩種珊瑚砂模型中：5#模型更適合于爆炸沖擊下的數值計算；Perzyna 黏塑性帽蓋模型中的參數需依靠較多的實驗獲得，在工程應用中略顯不足，但其模型理論性較強，開發出考慮顆粒破碎的黏塑性帽蓋模型將更符合珊瑚砂的變形規律。通過基本力學實驗結果確定了珊瑚砂的本構參數，并通過數值計算獲得了以下結論。

（1）對于應變率效應較敏感的珊瑚砂，使用動態增強系數確定壓實段應力應變曲線的方法可以彌補SHPB 實驗中無法將珊瑚砂加載至密實狀態的不足，為擬合高應變率下的物態方程提供方法。

（2）在彈丸侵徹珊瑚砂的計算中，5#模型能表現出高應變率下的流動特征，對于珊瑚砂的應力場響應描述符合客觀規律。

（3）流體彈塑性模型更能反映干燥珊瑚砂在爆炸作用下應力波衰減時出現的歷時增長現象，Perzyna 黏塑性帽蓋模型因對顆粒破碎描述的缺失，計算應力波衰減時得到的誤差較大。

（4）不良級配的珊瑚砂由于最大、最小干密度相差較小，相對密實度對彈丸侵徹深度影響較小，但對于爆炸波的壓力峰衰減影響較大。