韓 波，鞠玉濤，許進(jìn)升，周長省

（南京理工大學(xué) 機(jī)械工程學(xué)院，南京210094）

固體火箭發(fā)動機(jī)裝藥在生產(chǎn)、運輸和使用過程中會承受到各種復(fù)雜載荷的作用，在這些載荷作用下裝藥表面可能會產(chǎn)生微裂紋，裂紋的萌生和發(fā)展會影響發(fā)動機(jī)內(nèi)彈道性能和發(fā)動機(jī)的使用安全性，因此建立推進(jìn)劑裂紋擴(kuò)展過程的數(shù)值仿真方法十分重要.目前國內(nèi)針對推進(jìn)劑斷裂過程的數(shù)值仿真尚不能準(zhǔn)確模擬裂紋擴(kuò)展過程中的應(yīng)力分布、裂紋走向和裂紋終止問題.目前模擬材料裂紋擴(kuò)展過程有虛擬裂紋閉合技術(shù)（VCCT）、擴(kuò)展有限元法（XFEM）和粘聚區(qū)模型（CZM）等多種方法.VCCT需要事先設(shè)定材料的裂紋擴(kuò)展路徑，后兩者可以模擬裂紋走向未知情況下的裂紋擴(kuò)展問題.CZM是一種基于能量平衡和材料損傷的裂紋擴(kuò)展模型，裂紋在擴(kuò)展過程中裂尖材料產(chǎn)生損傷，其力學(xué)性能下降，當(dāng)裂紋消耗能量等于裂紋擴(kuò)展能時，裂尖材料失效擴(kuò)展.裂尖材料的損傷避免了裂尖應(yīng)力的奇異性.國外自20世紀(jì)90年代開始在粘聚區(qū)模型方面開展了大量的研究工作，國內(nèi)相關(guān)學(xué)者近幾年也開始展開相應(yīng)的研究工作[1，2].本文使用粘聚區(qū)模型建立復(fù)合固體推進(jìn)劑裂紋開裂過程的物理和數(shù)學(xué)模型，并結(jié)合ABAQUS用戶自定義單元開發(fā)技術(shù)實現(xiàn)對復(fù)合固體推進(jìn)劑裂紋擴(kuò)展過程的數(shù)值仿真，以期為推進(jìn)劑裝藥結(jié)構(gòu)完整性及安全性分析提供理論支持.

1 物理數(shù)學(xué)模型

1.1 粘聚區(qū)本構(gòu)模型

粘聚區(qū)模型的提出起源于DUGDALE的條狀屈服區(qū)模型和BARENBLAT提出的內(nèi)聚力模型.隨著粘聚區(qū)理論的發(fā)展，之后又提出了各種形式的裂尖粘聚區(qū)應(yīng)力分布模型，并且在瀝青、金屬等材料上取得了成功應(yīng)用.粘聚區(qū)模型（圖1所示）中材料的實際裂尖位于裂尖損傷位移δc處，假設(shè)的損傷裂尖位于損傷應(yīng)力最大處.為了準(zhǔn)確描述裂尖的損傷應(yīng)力變化情況，需要建立裂尖損傷應(yīng)力和損傷位移之間的對應(yīng)關(guān)系——粘聚區(qū)模型本構(gòu).

圖1 粘聚區(qū)模型示意圖

國外有許多學(xué)者提出了不同的粘聚區(qū)本構(gòu)形式，主要有指數(shù)勢函數(shù)形式、多項式形式、線性形式.其中線性粘聚區(qū)本構(gòu)通過合理調(diào)整本構(gòu)參數(shù)可以有效避免人工柔量問題，并且在粘彈性瀝青材料上獲得了成功應(yīng)用[3，4].圖2為線性粘聚區(qū)模型示意圖.在二維情況下，材料裂尖非線性區(qū)的裂紋張開有效位移和有效應(yīng)力定義為

式中，δt和δn為裂尖的切向和法向位移，σt和σn為裂尖的切向和法向應(yīng)力.線性粘聚區(qū)模型中假設(shè)有效位移和有效應(yīng)力關(guān)系呈現(xiàn)為2個線性階段，如圖2所示，其中σmax為材料的應(yīng)力損傷初始值；δcc為材料的損傷初始位移；δc為材料開裂的最終擴(kuò)展位移；Gc為材料的斷裂能，即圖2中三角形面積.當(dāng)裂尖材料的有效應(yīng)力σe達(dá)到損傷應(yīng)力σmax后，材料承載能力下降，當(dāng)裂尖材料總擴(kuò)展位移達(dá)到δc時，其消耗的能量等于材料的斷裂能Gc，之后材料發(fā)生完全斷裂.

圖2 線性粘聚區(qū)本構(gòu)示意圖

線性粘聚區(qū)本構(gòu)的具體表達(dá)形式為

粘聚區(qū)模型的引入使有限元模型中增加了人工柔量項，造成了整體結(jié)構(gòu)剛度下降.為了便于討論，以圖3所示的結(jié)構(gòu)為例，上下為2個正常的實體單元，中間加入粘結(jié)單元.實體單元初始長度為d，粘結(jié)單元初始厚度為0.在載荷F作用下單元沿上下方向伸長，實體單元伸長Δ，粘結(jié)單元伸長Δc.設(shè)實體單元和粘結(jié)單元剛度分別為Ks和Kc，Es為實體單元的模量.

圖3 粘結(jié)單元柔量分析示意圖

在粘結(jié)單元未損傷前，由力平衡方程可得：

式中，λ＝δcc/δc，該結(jié)構(gòu)的整體剛度為

不加粘結(jié)單元的理論剛度為Ki＝Ks/2，由式（3）、式（4）得：

從式（5）可以發(fā)現(xiàn)，整體剛度與基體材料的模量Es、網(wǎng)格大小d、粘聚區(qū)本構(gòu)參數(shù)δc、σmax和λ相關(guān)，其中Es，δc和σmax為材料的固有屬性，不可改變.增大網(wǎng)格可以增大系統(tǒng)的剛度，但是會造成計算精度的下降.λ表征了粘聚區(qū)本構(gòu)中初始上升段的斜率，通過減小λ可以提高粘結(jié)單元的初始剛度，從而有效避免過大的人工柔量帶來的問題.λ的選擇可以通過加入粘結(jié)單元和未加入粘結(jié)單元情況下的仿真結(jié)果對比確定，λ的取值應(yīng)保證粘結(jié)單元的加入不會對材料未產(chǎn)生損傷前的系統(tǒng)剛度產(chǎn)生較大影響，從而保證計算模型具有合理的精度.

1.2 推進(jìn)劑本構(gòu)模型

HTPB推進(jìn)劑是一種典型的粘彈性材料，在推進(jìn)劑裝藥結(jié)構(gòu)完整性分析中廣泛采用線性粘彈性本構(gòu)模型，該模型可以較好地反映推進(jìn)劑的力學(xué)特性，因此本文在計算中使用線粘彈性本構(gòu)模型.線粘彈性材料的本構(gòu)方程可以寫成：

式中，G（t）和K（t）為剪切模量和體積模量，可表示為

式中，E（t）為楊氏松弛模量，ν為泊松比.E（t）可以寫成Prony級數(shù)的形式：

式中，τi為Prony級數(shù)中的松弛時間.

2 有限元離散及仿真方法研究

ABAQUS提供了豐富的材料和單元庫，但是ABAQUS材料和單元庫中并不包含特定粘聚區(qū)模型和本構(gòu)關(guān)系，需要用戶進(jìn)行開發(fā).ABAQUS提供了用戶自定義單元的接口程序UEL，用戶可以根據(jù)需要自定義各種新單元.下面以二維粘聚區(qū)模型為例，給出ABAQUS二次開發(fā)所需增量形式的粘聚區(qū)單元建立過程.由有限元理論可知，利用Newton－Raphson法求解非線性有限元問題時需要給定單元的切線剛度矩陣KT和節(jié)點平衡矢量列陣R.

Newton法的迭代公式為

式中，

式中，ce為單元選擇矩陣，a為單元位移向量，V為被積單元體積，σ為單元應(yīng)力.

材料的Jacobian矩陣定義為

圖4為2個平面三角形單元和一個粘結(jié)單元變形示意圖.

單元節(jié)點1、4和2、3之間的相對位移為

單元的積分點法向和切向位移表示為

式中，a為圖3中節(jié)點在系統(tǒng)坐標(biāo)下的坐標(biāo)值，R為坐標(biāo)轉(zhuǎn)換矩陣，N為插值形函數(shù).

圖4 粘結(jié)單元

式中，ξ為變換后的坐標(biāo).

本文根據(jù)上述數(shù)學(xué)模型編制了ABAQUS用戶自定義單元開發(fā)程序UEL，建立了復(fù)合推進(jìn)劑裂紋擴(kuò)展數(shù)值仿真方法.

3 仿真與分析

3.1 仿真模型

為了驗證本文所建立的仿真計算方法的可行性，利用文獻(xiàn)[6]中的模型，對 HTPB推進(jìn)劑Ⅰ－Ⅱ型裂紋進(jìn)行了數(shù)值仿真計算，圖5為仿真模型示意圖.模型寬度W＝50 mm，長度H＝100 mm，中心裂紋長度l＝20mm.藥柱下表面固定，上表面施加60mm/min等速拉伸載荷.HTPB推進(jìn)劑松弛模量通過松弛實驗獲取，Prony級數(shù)參數(shù)見表1，初始模量E0＝15 MPa，泊松比取0.499.粘聚區(qū)本構(gòu)使用式（3）所示的形式，本文主要研究 HTPB復(fù)合推進(jìn)劑的裂紋擴(kuò)展有限元計算方法，粘聚區(qū)本構(gòu)參數(shù)的具體實驗獲取不在本文研究范圍之內(nèi).根據(jù)單軸拉伸實驗，裂尖損傷應(yīng)力σmax大致取0.5 MPa，推進(jìn)劑斷裂能Gc根據(jù)文獻(xiàn)[5]大致取500J/m2.

圖5 仿真模型示意圖

表1 松弛模量數(shù)據(jù)

使用粘聚區(qū)模型模擬裂紋擴(kuò)展方向未知情況下的材料開裂過程需要在正常實體單元之間加入粘結(jié)單元，通過ABAQUS CAE無法實現(xiàn).本文采用MATLAB編程語言生成包含粘結(jié)單元和二維實體單元的有限元網(wǎng)格.如圖6所示，生成的網(wǎng)格中包含10 475個粘結(jié)單元和7 052個三角形實體單元.為了準(zhǔn)確模擬裂紋擴(kuò)展過程中的應(yīng)力、應(yīng)變變化情況，在預(yù)測的裂紋擴(kuò)展路徑四周進(jìn)行了網(wǎng)格細(xì)化.

圖6 有限元網(wǎng)格

3.2 仿真結(jié)果分析

為了確定合理的λ以避免人工柔量的影響，本文對比了無粘結(jié)單元情況下和加入了粘結(jié)單元情況下的仿真計算結(jié)果，如圖7所示，圖中，F(xiàn)s為圖6中有限元模型計算的上下表面拉力.

圖7 不同λ下的載荷－時間曲線

從圖7可以發(fā)現(xiàn)λ的值明顯影響了計算所得的載荷－時間曲線.當(dāng)不插入粘結(jié)單元時，載荷隨時間持續(xù)上升，推進(jìn)劑不會產(chǎn)生斷裂.粘結(jié)單元的引入可以計算出推進(jìn)劑從裂尖產(chǎn)生局部損傷直至斷裂的整個過程.當(dāng)λ＝0.01時，從圖上可以看出粘聚區(qū)模型計算結(jié)果和無粘聚區(qū)模型下的計算結(jié)果相差很大，其裂紋未擴(kuò)展前的曲線斜率明顯低于后者，這是由于過大的人工柔量導(dǎo)致了系統(tǒng)剛度的下降.當(dāng)λ＝0.001時，在3s之前粘聚區(qū)模型的使用與否對計算結(jié)果影響不大，說明該值可以較好地反映粘聚區(qū)本構(gòu)模型和裂紋擴(kuò)展的基本規(guī)律.在3s之后由于推進(jìn)劑裂尖損傷的產(chǎn)生，材料承載能力開始下降，此時常規(guī)有限元方法已經(jīng)不能反映出裂紋的擴(kuò)展過程.在4.56s時載荷達(dá)到最大值，之后推進(jìn)劑失穩(wěn)快速開裂.

圖8和圖9為在λ＝0.001情況下，3.07s和4.56s時的有限元網(wǎng)格變形情況.在3.07s時推進(jìn)劑裂尖開始產(chǎn)生了應(yīng)力損傷情況，此時裂紋將要產(chǎn)生擴(kuò)展，損傷裂尖位于裂尖的初始位置.從圖7上可以看出4.56s時仿真模型所受到的載荷達(dá)到最大值.對比圖9，發(fā)現(xiàn)推進(jìn)劑產(chǎn)生了明顯的擴(kuò)展裂紋，損傷裂尖位置基本上達(dá)到了推進(jìn)劑的邊緣，此時推進(jìn)劑的真實裂尖并不位于損傷裂尖，推進(jìn)劑仍能承受外載荷的作用，但是與未受損傷前相比其承載能力已經(jīng)明顯下降.

圖8 3.07s有限元網(wǎng)格

圖9 4.56s有限元網(wǎng)格

圖10為仿真和文獻(xiàn)[6]中實驗獲得的拉伸破壞后的HTPB推進(jìn)劑裂紋擴(kuò)展路徑.通過實驗獲得的初始裂紋起裂角平均值為51.9°，之后裂紋擴(kuò)展方向轉(zhuǎn)變?yōu)槠街绷鸭y，而仿真所獲得的初始裂紋起裂角約為62°，之后裂紋轉(zhuǎn)變?yōu)槠街绷鸭y.由于粘聚區(qū)有限元仿真中裂紋沿單元界面開裂，裂尖的網(wǎng)格細(xì)密程度決定了仿真獲得的初始裂紋擴(kuò)展角精度.本文仿真所獲得的初始裂紋擴(kuò)展角與實驗結(jié)果的差距部分主要是由裂尖網(wǎng)格劃分精度導(dǎo)致.

圖10 仿真和實驗開裂路徑

文獻(xiàn)[6]中觀察到“裂紋初始擴(kuò)展后都有向橫向裂紋的轉(zhuǎn)變趨勢”，仿真結(jié)果也出現(xiàn)裂紋初始擴(kuò)展之后轉(zhuǎn)變?yōu)闄M向裂紋的現(xiàn)象.圖10中的裂紋擴(kuò)展路徑基本重合，表明粘聚區(qū)模型可以較好地預(yù)測HTPB復(fù)合裂紋的裂紋擴(kuò)展路徑.

圖11為裂紋擴(kuò)展路徑上距離初始裂尖不同位置處的Von Mises應(yīng)力－時間曲線，所示的幾條曲線均呈現(xiàn)出先增大后減小的形狀，曲線的峰值對應(yīng)損傷裂尖到達(dá)該點的時間.從圖上可以看出裂尖位置隨時間的變化情況，圖中4個位置達(dá)到應(yīng)力最大值的時間間隔呈現(xiàn)出逐漸縮小的趨勢，這說明裂紋呈現(xiàn)出加速擴(kuò)展的趨勢.

本文建立了一種復(fù)合固體推進(jìn)劑開裂過程數(shù)值仿真方法，為了能更加準(zhǔn)確地描述出復(fù)合推進(jìn)劑的裂紋擴(kuò)展過程，需要通過大量的實驗建立起準(zhǔn)確的推進(jìn)劑粘聚區(qū)模型，這也是筆者下一階段的工作目標(biāo).

圖11 裂尖擴(kuò)展路徑上不同位置的應(yīng)力－時間曲線

4 結(jié)論

本文采用粘聚區(qū)模型理論針對復(fù)合固體推進(jìn)劑裂紋擴(kuò)展過程進(jìn)行了研究，建立了針對推進(jìn)劑斷裂過程的物理和數(shù)學(xué)模型，并編制了ABAQUS二次開發(fā)程序，實現(xiàn)了對固體推進(jìn)劑Ⅰ－Ⅱ型復(fù)合裂紋擴(kuò)展過程的數(shù)值仿真計算，得到了如下結(jié)論：

①使用粘聚區(qū)模型可以很好地模擬出復(fù)合固體推進(jìn)劑裂尖的損傷應(yīng)力場、裂紋擴(kuò)展路徑和裂紋體開裂過程.

②粘聚區(qū)模型可以為固體推進(jìn)劑裝藥完整性和安全性分析提供一種可靠的分析計算方法.

[1]劉陸廣，歐卓成，段卓平，等.混凝土動態(tài)斷裂數(shù)值模擬[J].兵工學(xué)報，2010，31（6）：741－745.LIU Lu－guang，OU Zhuo－cheng，DUAN Zhuo－ping，et al.Simulation of concrete dynamic fracture[J].Acta Armamentarii，2010，31（6）：741－745.（in Chinese）

[2]崔浩，李玉龍，劉元鏞，等.基于粘聚區(qū)模型的含填充區(qū)復(fù)合材料接 頭 失 效 數(shù) 值 模 擬 [J].復(fù) 合 材 料 學(xué) 報，2010，27（2）：161－168.CUI Hao，LI Yu－long，LIU Yuan－yong，et al.Numerical simulation of composites joints failure based on cohesive zone model[J].Acta Materiae Composite Sinica，2010，27（2）：161－168.（in Chinese）

[3]SONG S H，PAULINO G H，BUTTLAR W G.A bilinear cohesive zone model tailored for fracture of asphalt concrete considering viscoelastic bulk material[J].Engineering Fracture Mechanics，2006，73（18）：2 829－2 848.

[4]SONG S H，PAULINO G H，BUTTLAR W G.Simulation of crack propagation in asphalt concrete using an intrinsic cohesive zone model[J].Journal of Engineering Mechanics，2006，132（11）：1 215－1 223.

[5]MAROM G，HAREL H，ROSNER J.Fracture energies of composite propellants[J].Jounal of Applied Polymer Science，1977，21（6）：1 629－1 634.

[6]張亞，強(qiáng)洪夫，楊月誠.國產(chǎn) HTPB復(fù)合固體推進(jìn)劑Ⅰ－Ⅱ型裂紋斷裂性能實驗研究[J].含能材料，2007，15（4）：359－361.ZHANG Ya，QIANG Hong－fu，YANG Yue－cheng.Fracture behavior of HTPB composite propellant inⅠ－Ⅱ mixed mode crack[J].Chinese Journal of Energetic Materials，2007，15（4）：359－361.（in Chinese）