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

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

固體火箭發(fā)動(dòng)機(jī)裝藥在生產(chǎn)、運(yùn)輸和使用過程中會(huì)承受到各種復(fù)雜載荷的作用，在這些載荷作用下裝藥表面可能會(huì)產(chǎn)生微裂紋，裂紋的萌生和發(fā)展會(huì)影響發(fā)動(dòng)機(jī)內(nèi)彈道性能和發(fā)動(dòng)機(jī)的使用安全性，因此建立推進(jìn)劑裂紋擴(kuò)展過程的數(shù)值仿真方法十分重要.目前國內(nèi)針對(duì)推進(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ò)展能時(shí)，裂尖材料失效擴(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ù)實(shí)現(xiàn)對(duì)復(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所示）中材料的實(shí)際裂尖位于裂尖損傷位移δc處，假設(shè)的損傷裂尖位于損傷應(yīng)力最大處.為了準(zhǔn)確描述裂尖的損傷應(yīng)力變化情況，需要建立裂尖損傷應(yīng)力和損傷位移之間的對(duì)應(yīng)關(guān)系——粘聚區(qū)模型本構(gòu).

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

國外有許多學(xué)者提出了不同的粘聚區(qū)本構(gòu)形式，主要有指數(shù)勢函數(shù)形式、多項(xiàng)式形式、線性形式.其中線性粘聚區(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個(gè)線性階段，如圖2所示，其中σmax為材料的應(yīng)力損傷初始值；δcc為材料的損傷初始位移；δc為材料開裂的最終擴(kuò)展位移；Gc為材料的斷裂能，即圖2中三角形面積.當(dāng)裂尖材料的有效應(yīng)力σe達(dá)到損傷應(yīng)力σmax后，材料承載能力下降，當(dāng)裂尖材料總擴(kuò)展位移達(dá)到δc時(shí)，其消耗的能量等于材料的斷裂能Gc，之后材料發(fā)生完全斷裂.

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

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

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

圖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)的剛度，但是會(huì)造成計(jì)算精度的下降.λ表征了粘聚區(qū)本構(gòu)中初始上升段的斜率，通過減小λ可以提高粘結(jié)單元的初始剛度，從而有效避免過大的人工柔量帶來的問題.λ的選擇可以通過加入粘結(jié)單元和未加入粘結(jié)單元情況下的仿真結(jié)果對(duì)比確定，λ的取值應(yīng)保證粘結(jié)單元的加入不會(huì)對(duì)材料未產(chǎn)生損傷前的系統(tǒng)剛度產(chǎn)生較大影響，從而保證計(jì)算模型具有合理的精度.

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

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

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

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

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

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

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

Newton法的迭代公式為

式中，

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

材料的Jacobian矩陣定義為

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

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

單元的積分點(diǎn)法向和切向位移表示為

式中，a為圖3中節(jié)點(diǎn)在系統(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 仿真模型

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

圖5 仿真模型示意圖

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

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

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

3.2 仿真結(jié)果分析

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

圖7 不同λ下的載荷－時(shí)間曲線

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

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

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

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

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

圖10 仿真和實(shí)驗(yàn)開裂路徑

文獻(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)力－時(shí)間曲線，所示的幾條曲線均呈現(xiàn)出先增大后減小的形狀，曲線的峰值對(duì)應(yīng)損傷裂尖到達(dá)該點(diǎn)的時(shí)間.從圖上可以看出裂尖位置隨時(shí)間的變化情況，圖中4個(gè)位置達(dá)到應(yīng)力最大值的時(shí)間間隔呈現(xiàn)出逐漸縮小的趨勢，這說明裂紋呈現(xiàn)出加速擴(kuò)展的趨勢.

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

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

4 結(jié)論

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

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

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

[1]劉陸廣，歐卓成，段卓平，等.混凝土動(dòng)態(tài)斷裂數(shù)值模擬[J].兵工學(xué)報(bào)，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é) 報(bào)，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)劑Ⅰ－Ⅱ型裂紋斷裂性能實(shí)驗(yà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）