基于粘聚區(qū)模型的推進劑開裂數(shù)值仿真

韓 波，鞠玉濤，許進升，周長省

（南京理工大學 機械工程學院，南京210094）

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

1 物理數(shù)學模型

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

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

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

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

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

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

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

粘聚區(qū)模型的引入使有限元模型中增加了人工柔量項，造成了整體結(jié)構(gòu)剛度下降.為了便于討論，以圖3所示的結(jié)構(gòu)為例，上下為2個正常的實體單元，中間加入粘結(jié)單元.實體單元初始長度為d，粘結(jié)單元初始厚度為0.在載荷F作用下單元沿上下方向伸長，實體單元伸長Δ，粘結(jié)單元伸長Δc.設實體單元和粘結(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和λ相關，其中Es，δc和σmax為材料的固有屬性，不可改變.增大網(wǎng)格可以增大系統(tǒng)的剛度，但是會造成計算精度的下降.λ表征了粘聚區(qū)本構(gòu)中初始上升段的斜率，通過減小λ可以提高粘結(jié)單元的初始剛度，從而有效避免過大的人工柔量帶來的問題.λ的選擇可以通過加入粘結(jié)單元和未加入粘結(jié)單元情況下的仿真結(jié)果對比確定，λ的取值應保證粘結(jié)單元的加入不會對材料未產(chǎn)生損傷前的系統(tǒng)剛度產(chǎn)生較大影響，從而保證計算模型具有合理的精度.

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

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

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

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

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

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

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

Newton法的迭代公式為

式中，

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

材料的Jacobian矩陣定義為

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

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

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

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

圖4 粘結(jié)單元

式中，ξ為變換后的坐標.

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

3 仿真與分析

3.1 仿真模型

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

圖5 仿真模型示意圖

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

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

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

3.2 仿真結(jié)果分析

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

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

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

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

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

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

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

圖10 仿真和實驗開裂路徑

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

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

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

圖11 裂尖擴展路徑上不同位置的應力－時間曲線

4 結(jié)論

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

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

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

[1]劉陸廣，歐卓成，段卓平，等.混凝土動態(tài)斷裂數(shù)值模擬[J].兵工學報，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ū)復合材料接 頭 失 效 數(shù) 值 模 擬 [J].復 合 材 料 學 報，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]張亞，強洪夫，楊月誠.國產(chǎn) HTPB復合固體推進劑Ⅰ－Ⅱ型裂紋斷裂性能實驗研究[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）