邸德寧，陳小偉

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

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

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

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

1 碎片云幾何表現(xiàn)的差異

1.1 計算模型設(shè)置

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

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

表1 各算例的設(shè)置Table 1 Settings of various cases

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

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

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

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

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

圖1 算例01～03數(shù)值模擬結(jié)果Fig.1 Simulation results of cases 01-03

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

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

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

2 材料失效細觀表現(xiàn)的差異

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

2.1 有、無失效模型計算結(jié)果的差異

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

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

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

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

2.2 最大拉應(yīng)力及Grady失效模型結(jié)果比較

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

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

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

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

綜上，從粒子點壓力表現(xiàn)來看，在本文所采用失效參數(shù)下，最大拉應(yīng)力模型下材料更難失效。

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

3 碎片數(shù)量分布差異

3.1 碎片數(shù)量分布的整體比較

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

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

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

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

3.2 大質(zhì)量碎片分布差異比較

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

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

工況Grady模型最大拉應(yīng)力模型實驗值12.822.512.9526.235.926.3631.301.281.45

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

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

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

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

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

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

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

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

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

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

Grady模型中臨界應(yīng)變常數(shù)εc直接影響失效閾值計算值，其對所有材料被建議取值0.15[7]，對于鋁合金材料已被驗證該取值較為適用。考慮到銅材料實驗失效應(yīng)力為1.0～2.5 GPa，而εc=0.15時Grady模型計算值僅為1.0 GPa[7]，因此本文以O(shè)FHC為例，嘗試ε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狀態(tài)方程和Steinberg-Guinan強度模型。

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

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

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

圖9 碎片云碎片識別結(jié)果側(cè)視圖Fig.9 Side view of debris recognition result of debris clouds

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

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

6 結(jié) 論

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

在本文所用失效參數(shù)下，最大拉應(yīng)力模型下材料更難失效，導(dǎo)致碎片云擴散程度小幅減弱，小碎片聚集為大碎片，大碎片質(zhì)量變大，碎片云侵徹極限提高。增大失效閾值同樣導(dǎo)致材料更難失效，亦有上述表現(xiàn)。