碎片云SPH方法數值模擬中的材料失效模型*

邸德寧，陳小偉

(1.中國工程物理研究院總體工程研究所，四川 綿陽 621999; 2.北京理工大學前沿交叉科學研究院，北京 100081; 3.北京理工大學爆炸科學與技術國家重點實驗室，北京 100081)

光滑粒子流體動力學法(smoothed particle hydrodynamics, SPH)是一種Lagrange無網格粒子算法，超高速撞擊中的強沖擊波使材料表現出如流體一樣的性質，伴隨著材料大變形甚至破碎，傳統網格方法難以求解，而基于流體動力學控制方程的SPH方法非常適用[1]。在SPH數值模擬前的處理中，必須設定材料模型，包括狀態方程、強度模型和失效模型等。材料狀態方程不可或缺，描述材料密度、壓力和內能之間的關系；材料強度一般不可忽略，引入強度模型可計算材料偏應力張量和材料當前應變、應變率、溫度等對屈服強度的影響。薄板超高速撞擊中材料破碎形成碎片云的過程，被認為是稀疏波導致的拉應力超過材料斷裂應力后材料不斷層裂破壞的過程，而狀態方程和強度模型均不能表征材料斷裂方面的特征，因此需要考慮在計算中引入失效模型。

利用SPH方法對碎片云進行數值模擬，計算結果與實驗吻合很好，但同時發現材料模型的選擇對計算結果有影響。Groenenboom[2]發現SESAME狀態方程在壓力和徑向速度方面稍小于多項式狀態方程；張偉等[3]發現Tillotson和Puff狀態方程相比于Shock狀態方程，碎片云徑向速度稍大，碎片云顆粒更細小且分散，而強度模型的選擇對結果影響不大；Higashide等[4]發現不包含材料相變的狀態方程計算得到的碎片最大，考慮材料液化和氣化的熱力學完全狀態方程對相態分布計算較準確，其中考慮了亞穩態相態情況的狀態方程計算結果最準確；李寶寶等[5]通過比較GRAY三相物態方程與Tillotson物態方程，也發現材料相變對碎片云外形尺寸影響較大。然而，SPH方法失效模型的選擇尚無具體結論，有研究者認為SPH方法中不需要材料失效模型，另一些研究者則采用最大拉應力或Grady失效模型。另外，現有對材料模型方案的研究大多針對碎片云外形或運動狀態，未能很好地討論對碎片質量和數量分布的影響，以及影響程度與撞擊工況的相關性，難以實際參考應用。

本文中，利用AUTODYN軟件中的SPH模塊，針對無失效模型、Grady失效模型和最大拉應力失效模型3種方案，考察失效模型對碎片云外形尺寸及速度、粒子點失效表現和碎片數量分布的影響，討論失效模型的影響程度隨撞擊工況的變化，并比較Grady失效模型下不同失效閾值導致的碎片云差異。本文旨在補充SPH方法中材料模型的選擇，對碎片云數值模擬特別是碎片分布進行深入分析。

1 碎片云幾何表現的差異

1.1 計算模型設置

本文中所用算例共2部分，第1部分參照實驗情況[6]，取直徑為9.53 mm的球形彈丸，材料為2017-T4鋁合金，采用Shock狀態方程和Johnson-Cook強度模型；薄板材料為6061-T6鋁合金，采用Shock狀態方程和Steinberg-Guinan強度模型；彈丸正撞擊薄板，表1中列出了共3種工況7個算例的薄板厚度、撞擊速度及所采用的失效模型方案。

最大拉應力模型即材料拉應力超過設定的失效閾值后材料失效，本文中通過主應力失效實現，設定6061-T6鋁合金的失效閾值為2.6 GPa[3]，2017-T4鋁合金的失效閾值為2.5 GPa[6]，不考慮最大剪切力失效情況。

表1 各算例的設置Table 1 Settings of various cases

Grady失效模型同樣基于主應力判斷失效，即主應力超過失效閾值ps后材料失效，不同的是失效閾值ps非固定值，而由下式計算得到：

式中：ρ為材料密度，由材料連續性方程或核函數更新；Y為屈服強度，由材料強度模型根據材料當前應變、應變率和溫度計算更新，變化幅度較大，對ps計算值影響嚴重；c0為材料體積聲速，由材料體積模量和當前密度決定；εc為臨界失效應變[7]，材料常數，鋁合金材料一般取εc=0.15?？梢姡诓煌瑓^域、不同時刻，ps的取值均可能不同。

1.2 不同失效模型方案下碎片云幾何表現的比較

為分析失效準則對碎片云計算結果的影響，設計了算例01～03作對比，算例01不引入失效模型，算例02采用Grady模型，算例03采用最大拉應力模型。另外，文中碎片云圖像及分布統計均基于撞擊后20 μs時的計算結果。

如圖1對照所示，有、無失效模型對碎片云外形影響嚴重，特別是對彈丸碎片云部分；而Grady失效模型和最大拉應力失效模型的選擇并無明顯差別。測量碎片云整體尺寸，發現：算例01徑向尺寸為算例02的81%，軸向尺寸為算例02的96%；算例03和算例02碎片云整體尺寸相差很小。

圖1 算例01～03數值模擬結果Fig.1 Simulation results of cases 01-03

針對彈丸碎片云部分，圖2為算例01～07彈丸碎片主體碎片識別結果放大圖，其中碎片識別為只顯示包含有未失效粒子的碎片。對比各圖可見，算例01中粒子聚集緊密形成一個大碎片，包含了彈丸中94.8%的粒子點；算例02和算例03結果相近，碎片孤立并高度分散，但算例02中碎片更多，且最前端平整，更接近實驗圖像。算例01中碎片云尺寸明顯偏小。分析算例02～07碎片云特征速度[6]發現，最大拉應力模型下碎片擴散程度不及Grady模型，即前端和徑向速度偏小但后端速度偏大，不過兩者相差不大，都和實驗很接近。

綜上，不引入失效準則的計算結果明顯與實驗不符，有失效模型的計算結果均與實驗符合，但最大拉應力模型下碎片云的擴散程度稍弱于Grady模型下。

圖2 算例01～07彈丸碎片云碎片識別結果與實驗對照[6]Fig.2 Main parts of projectile fragments of cases 01-07 compared with experiments[6]

2 材料失效細觀表現的差異

失效模型方案不同導致碎片云計算結果有差異，為分析其內在機理，在計算模型中設置測量點，提取粒子點壓力-時間歷程作對比分析。

2.1 有、無失效模型計算結果的差異

圖3 測量點A、B壓力時間歷程Fig.3 Pressure-time curves of gauge points A and B

首先對比有、無失效模型的情況，在算例01和算例02中設置測量點A和B，圖3繪制了測量點壓力-時間歷程。在沖擊波作用段內同一測量點曲線完全重合，表明材料被沖擊過程與失效模型無關。稀疏波到達一段時間后開始出現小的分歧，進入負壓后曲線出現明顯分歧。觀察右上的放大圖可見，算例01中測量點負壓均超過6 GPa，并在達到峰值后圓滑回落；算例02中A點負壓在到達2.86 GPa后突然跳到0 GPa并保持，B點負壓在到達峰值2.12 GPa后圓滑回落，并在后期繼續變化。上述差異正是由失效模型導致的，不引入失效模型時粒子點在強拉力下不會失效，從而完整地傳遞稀疏波，但這和工程材料拉伸斷裂的實際表現明顯不符；引入失效模型后，拉伸主應力超過失效閾值即會判為失效，粒子不再能承受拉應力或剪切力，如算例02中A點表現所示，負壓跳變為0 GPa。

考慮失效模型對碎片云宏觀表現的影響，當粒子點失效后不再承受拉應力，效果即為斷裂點，會阻礙稀疏波的傳播，從而導致斷裂點兩側粒子運動狀態差異較大，碎片高度分散。大量失效粒子會構成斷裂面，一個孤立碎片的表面即為失效粒子組成的斷裂面，碎片內部粒子表現如同算例02中B點，應力波在碎片內部振蕩導致其壓力值不為零。反之，不引入失效模型，則不會出現斷裂點，通過SPH算法核函數近似，近鄰粒子之間無障礙地相互影響，導致粒子點速度矢量比較接近從而集中分布，出現圖2中算例01所示聚集現象。但核函數近似的影響域有限大，若2個粒子速度差異較大而距離較遠，則難以相互影響，也會產生相對孤立的碎片，形式上表現為碎片云。

綜上，不引入失效模型也能得到相對孤立的碎片，但只是形式上表現為碎片云，材料表現與實際不符；引入失效模型后粒子點失效表現才符合物理實際，得到的碎片分布符合實驗結果。

2.2 最大拉應力及Grady失效模型結果比較

進一步考慮，最大拉應力模型中材料失效閾值恒定，但Grady模型中隨著材料當前狀態的改變，閾值ps的計算值也會變化。圖4示例了2種失效模型下某粒子點處失效閾值變化及主應力失效規則的觸發，可見Grady模型下失效閾值在沖擊波到達后開始劇烈變化，事實上主要隨屈服強度變化，密度和體積聲速隨時間變化不大，而粒子點在最大主應力超過失效閾值后立即被判為失效，應力值歸零。

圖4 材料失效閾值變化及對應失效規則觸發的示例Fig.4 An instance of material failure regulation and its failure threshold-time curve

在算例02和算例03中取不同位置處C、D、E等3個測量點，如圖5所示為3個點處壓力時間歷程，比較2種失效模型下粒子點失效表現的差異。C點粒子在Grady和最大拉應力模型下同時失效，但最大拉應力模型下負壓更大；D點粒子在最大拉應力模型下負壓也更大，而且比Grady模型下更晚失效；E點粒子在Grady模型下失效，但在最大拉應力模型下未失效。

由圖3可見算例02中B點未失效，但其負壓峰值遠低于無失效模型的算例01下，可見其他粒子的失效導致的壓力變化會影響B點壓力，失效準則的引入降低了粒子點負壓；無失效模型情況可視為粒子最難失效情況(不失效)，結合圖5中最大拉應力模型下粒子點負壓大于Grady模型下，由此推斷：材料更難失效的失效模型下粒子負壓更大。從粒子失效時間可得到相近結論，更難失效的模型下粒子更晚失效甚至不會失效。

綜上，從粒子點壓力表現來看，在本文所采用失效參數下，最大拉應力模型下材料更難失效。

圖5 測量點C、D、E壓力時間歷程Fig.5 Pressure-time curves of gauge points C, D and E

3 碎片數量分布差異

3.1 碎片數量分布的整體比較

粒子失效后成為斷裂點，孤立碎片的邊界必然全是已失效粒子，AUTODYN軟件基于此進行碎片識別，獲得SPH粒子和碎片的對應關系及碎片特征。本文中利用AUTODYN碎片識別功能，統計數值模擬結果中正向運動的彈丸碎片數量、質量。由于不引入失效模型的方案明顯與實驗不符，此處只對算例02～07進行比較，分析最大拉應力模型和Grady模型下結果的差異。

統計模擬結果中所有碎片，發現3個工況中Grady模型下碎片總數分別比最大拉應力模型下多1.78%、4.48%、0.42%，碎片越多意味著失效粒子越多，這說明Grady模型下材料更易失效。因此，從碎片總數來看，在本文所采用失效參數下，最大拉應力模型下材料更難失效。

統計質量大于0.01 mg的彈丸碎片，得圖6所示對數坐標下碎片無量綱質量m/M與累計數量Nc的關系，其中Nc為質量不小于m的碎片累計數量，M為0.01 mg以上碎片總質量，3種工況下依次取683、936、340 mg。需要注意，對數坐標下坐標分布非線性，小質量碎片段差異被壓縮。不同模型下曲線變化速度有差異，導致不同質量處碎片累計數量相對大小不一，lg(m/M)=-2.9標記線起不同模型下曲線基本一致，直到lgNc=1.1所標記的大碎片段位置。

從曲線整體來看，2種模型下曲線變化趨勢一致，各工況中Grady模型下曲線均接近直線，而最大拉應力模型下曲線波動大，偏離直線。不同工況下材料強度及沖擊波強度等不同，Grady模型能夠很好地與撞擊工況自適應，對應的調整失效閾值而在對數坐標下呈現出更好地線性分布，表現更佳。

3.2 大質量碎片分布差異比較

圖6中2種模型下曲線lgNc=1.1的大碎片段差異明顯，3種工況中均呈現Grady模型下碎片質量先更大后更小的規律。分析碎片位置發現Grady模型下將最大的若干個大碎片的部分材料分給了稍小的大碎片，使大碎片段質量分布更平緩。最大拉應力模型下粒子更聚集，導致最大的個別碎片偏大，其他大碎片偏小。另外，各工況中均為Grady模型下碎片總數更多，但在分界線處碎片數量基本一致，說明最大拉應力模型下是小碎片數量偏少。這同樣是最大拉應力模型下材料難失效導致，小碎片附著在大碎片上未剝離，小碎片數量大幅減小，同時增加大碎片質量。

圖6 0.01 mg以上碎片累計數量分布Fig.6 Distribution of cumulative number of debris above 0.01 mg

工況Grady模型最大拉應力模型實驗值12.822.512.9526.235.926.3631.301.281.45

最大碎片的侵徹性能可代表碎片云的侵徹極限，比較發現：各工況中Grady模型下最大碎片的質量分別比最大拉應力模型下小9.84%、18.5%、19.2%，而軸向速度基本一致。綜上可推斷：Grady模型下得到的碎片云侵徹極限更小，對后續靶板的損傷更均勻；最大拉應力模型下后續靶板損傷更嚴重。

測量計算結果中最大碎片在各坐標方向尺寸并計算其球等效直徑，如表2所示為最大碎片球等效直徑計算值與實驗值的對比。雖然最大拉應力模型下最大碎片的質量較大，但由于其粒子聚集度高，碎片等效直徑卻稍小。不難看到，所有工況中Grady模型下最大碎片等效直徑計算值更接近實驗值，Grady模型計算結果更符合實驗。

4 失效模型影響隨撞擊工況變化

根據各工況板厚及撞擊速度可知，工況1下材料初步破碎，工況3下材料充分破碎，工況2介于前兩者之間。觀察圖6中各工況下兩曲線間的差異，材料破碎越充分，曲線間差異越小，失效模型方案對計算結果影響越弱，從碎片總數也能得到此結論。

從材料失效考慮，當材料破碎不充分時，稀疏波強度與失效閾值相當，因此失效閾值上細微的差別就能導致很大的碎片分布差異；材料破碎充分時，稀疏波強度遠大于失效閾值，失效模型影響便較弱。另外，Grady模型中失效閾值主要隨材料屈服強度改變，撞擊較弱時應變、應變率等較小，導致失效閾值計算值偏小。Grady模型下材料本就更易失效，更加偏小的失效閾值就會和最大拉應力模型下恒定的閾值差距更大，加劇兩種模型下碎片分布差異。

撞擊速度對稀疏波強度影響顯著，工況2撞擊速度遠低于工況1，而工況3的稍大于工況1。因此，工況2相比于工況1和3，不同失效模型下碎片數量曲線之間差異明顯，且最大碎片顯著減小。板厚主要改變強稀疏波作用時長和作用區域大小，嚴重影響小碎片數量，但彈丸材料總量不變，導致工況3下碎片總數顯著增大，曲線斜率明顯變化，但對曲線間差異影響不明顯。

從細觀上分析失效模型方案對粒子點失效的影響，在彈丸和薄板模型中軸線遍布測量點，圖7為只在一種失效模型下失效而在另一種模型下不失效的測量點的軸向坐標Z，0點到彈丸尾端和薄板背面等距。圖7中不同模型下失效表現相反的最遠粒子的坐標，表征了失效模型對材料失效表現產生明顯影響的起始位置。工況2中表現相反的最遠粒子坐標最大，工況3中則坐標最小，正是上述的稀疏波強度和失效閾值的比較。工況2下初始稀疏波強度即和失效閾值相當，邊緣位置處失效模型的影響已展現出；工況3下需要稀疏波不斷衰弱直至中心附近才能表現出失效模型的影響；工況1下最遠粒子坐標與工況3接近而與工況2差距較大，主要是由撞擊速度的差異導致的。

值得關注的是，圖6中大碎片段碎片質量出現明顯分歧，但開始分歧位置均在lgNc=1.1標記線附近，說明撞擊工況對大碎片段差異大小影響弱。大碎片多靠近彈丸中心，材料不斷層裂破碎同時降低稀疏波強度，導致大碎片位置時稀疏波強度均較低，因此撞擊工況對大碎片段差異影響不顯著。

圖7 模型中軸線上不同失效模型下失效表現不同的粒子點Fig.7 Particles with different failure performance under Grady and max-tension models

5 失效閾值對碎片云的影響

Grady模型中臨界應變常數εc直接影響失效閾值計算值，其對所有材料被建議取值0.15[7]，對于鋁合金材料已被驗證該取值較為適用。考慮到銅材料實驗失效應力為1.0～2.5 GPa，而εc=0.15時Grady模型計算值僅為1.0 GPa[7]，因此本文以OFHC為例，嘗試εc=0.15,0.30,0.45和0.60，考察失效閾值對碎片云的影響，同時確認銅材料合適的εc取值。

第2部分算例參照實驗[8]取直徑7.72 mm、厚2.36 mm的OFHC圓盤正撞擊1 mm厚的6061-T6鋁薄板，彈速為6.39 km/s，圓盤軸線與速度方向呈4.6°夾角，采用Shock狀態方程和Steinberg-Guinan強度模型。

圖8為圓盤中心粒子失效時刻壓力的變化曲線，隨著εc增大，粒子負壓增大且更晚失效。臨界應變常數εc越大，失效閾值越大，材料越難失效。比較碎片云特征速度發現：隨著εc增大，碎片云最前端速度減小1.3%，最外側徑向速度減小3.5%，失效閾值的增大導致碎片云擴散程度下降。但由于邊緣位置稀疏波強度很高，失效閾值對碎片云速度影響非常小。事實上，各εc取值下模擬的碎片云擴散程度均不及實驗，因此εc=0.15時速度結果最接近實驗，但不同取值時碎片云特征速度變化不超過3%，合適的εc取值需要根據具體算例中其他材料參數和人工黏性等確定。

圖8 圓盤中心粒子失效時刻壓力的變化Fig.8 Pressure-time curves of the center particles of the disc projectile

εc碎片總數最大碎片/mg0.15177 83214.430.30161 61622.230.45152 95631.210.60145 48034.14

圖9 碎片云碎片識別結果側視圖Fig.9 Side view of debris recognition result of debris clouds

圖10 0.01 mg以上碎片累計數量分布Fig.10 Distribution of cumulative number of debris above 0.01 mg

綜上，增大Grady模型下臨界應變常數εc直接導致材料更難失效，宏觀表現為碎片云擴散程度小幅減弱，碎片總數減少但0.01 mg以上碎片增加，最大碎片質量明顯增大，和前文中采用最大拉應力模型導致的材料更難失效表現基本一致。

6 結 論

比較發現，不引入失效模型時，材料在強負壓下不斷裂的表現不符合物理實際，彈丸碎片中粒子高度聚集，與實驗結果嚴重不符。因此，碎片云數值模擬中必須引入材料失效模型。有失效模型時，從碎片識別圖像、碎片分布曲線及最大碎片尺寸來看，Grady模型計算結果更符合實驗。不過失效模型對碎片云的影響程度與撞擊工況相關，材料破碎越充分，不同失效模型下碎片分布曲線之間差異越小。

在本文所用失效參數下，最大拉應力模型下材料更難失效，導致碎片云擴散程度小幅減弱，小碎片聚集為大碎片，大碎片質量變大，碎片云侵徹極限提高。增大失效閾值同樣導致材料更難失效，亦有上述表現。