HTPB推進(jìn)劑粘聚斷裂研究①

韓 波，鞠玉濤，周長(zhǎng)省

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

0 引言

固體推進(jìn)劑藥柱是固體火箭發(fā)動(dòng)機(jī)的重要部件，藥柱結(jié)構(gòu)完整性分析是保證固體火箭發(fā)動(dòng)機(jī)正常內(nèi)彈道性能的必要條件。藥柱表面微小裂紋的存在，對(duì)發(fā)動(dòng)機(jī)能否正常工作有重大影響。因此，研究推進(jìn)劑藥柱表面微小裂紋的起裂、擴(kuò)展規(guī)律具有十分重大現(xiàn)實(shí)意義。

目前，國(guó)內(nèi)外針對(duì)含微小裂紋的固體火箭發(fā)動(dòng)機(jī)藥柱完整性分析，主要采用傳統(tǒng)的斷裂力學(xué)準(zhǔn)則，如應(yīng)力強(qiáng)度因子、能量釋放率等。李九天等［1］使用三維奇異單元，計(jì)算了發(fā)動(dòng)機(jī)藥柱關(guān)鍵部位在內(nèi)壓載荷下的裂尖應(yīng)力強(qiáng)度因子。Earnest［2］使用實(shí)驗(yàn)和數(shù)值仿真結(jié)合的方法，計(jì)算了推進(jìn)劑的能量釋放率，并用于發(fā)動(dòng)機(jī)藥柱的安全性分析。但由于實(shí)際裂紋擴(kuò)展過(guò)程中材料產(chǎn)生了幾何不連續(xù)，因此傳統(tǒng)有限元方法不能準(zhǔn)確模擬出裂紋的擴(kuò)展過(guò)程。粘聚區(qū)模型(Cohesive Zone Model，簡(jiǎn)稱 CZM)最早由 Barenblatt［3］和 Dugdale［4］提出用于研究脆性材料和塑性材料開(kāi)裂過(guò)程。粘聚區(qū)模型為裂紋的起裂、擴(kuò)展分析提供了一種強(qiáng)有力的方法。20世紀(jì)90年代，國(guó)外學(xué)者將粘聚區(qū)模型成功用于有限元數(shù)值仿真中，并取得了成功應(yīng)用［5］。國(guó)內(nèi)劉陸廣等［6］使用粘聚區(qū)模型來(lái)仿真混凝土開(kāi)裂過(guò)程。崔浩等［7］將粘聚區(qū)模型成功用于復(fù)合材料接頭的失效分析中。粘聚區(qū)模型具有明確的物理意義和較高的計(jì)算精度。因此，可用于固體火箭發(fā)動(dòng)機(jī)藥柱裂紋開(kāi)裂和擴(kuò)展分析過(guò)程中。

1 粘聚區(qū)模型

HTPB推進(jìn)劑是一種高固體含量的高聚物復(fù)合材料。其中，包含大量直徑較大的高氯酸銨顆粒和直徑較小的鋁粉。在局部拉伸載荷作用下，HTPB推進(jìn)劑會(huì)在裂尖首先產(chǎn)生顆粒脫濕，進(jìn)而形成微孔洞，隨著孔洞的不斷擴(kuò)大和合并，逐漸形成了宏觀的裂紋。裂尖局部的顆粒脫濕和微孔洞稱為裂尖損傷區(qū)，損傷區(qū)的損傷演化規(guī)律影響著裂紋的起裂和擴(kuò)展過(guò)程。

圖1為粘聚區(qū)單元模型示意圖，(a)為常規(guī)有限元實(shí)體單元，單元之間共節(jié)點(diǎn)，無(wú)法直接模擬材料的宏觀開(kāi)裂過(guò)程;(b)為添加了粘結(jié)單元的有限元網(wǎng)格，實(shí)體單元之間通過(guò)粘結(jié)單元連接，粘結(jié)單元來(lái)模擬材料裂尖的損傷和演化直至宏觀裂紋形成，裂紋的可能擴(kuò)展路徑為實(shí)體單元之間的界面。

圖1 粘聚區(qū)模型示意圖Fig.1 Diagram of cohesive zone model

以平面Ⅰ-Ⅱ型復(fù)合裂紋為例，粘聚區(qū)單元有效位移和有效應(yīng)力定義為

式中 δt和δn為裂尖粘結(jié)單元面上的切向和法向位移;σt和σn為裂尖的切向和法向力。

σe和δe之間的變化關(guān)系，即粘聚區(qū)的本構(gòu)模型，它預(yù)示著材料裂尖的損傷和演化過(guò)程。圖2為冪指數(shù)粘聚區(qū)本構(gòu)的示意圖。

圖2中，橫坐標(biāo)代表材料裂尖損傷區(qū)變形量;縱坐標(biāo)代表裂尖損傷應(yīng)力。本構(gòu)關(guān)系中，上升段是為了數(shù)值模擬的需要，代表著裂尖材料未損傷前的響應(yīng);下降段代表著材料的損傷軟化形式。

式(2)為冪指數(shù)形式的粘聚區(qū)具體形式，可通過(guò)調(diào)節(jié)參數(shù)α來(lái)模擬不同形式的軟化形式。

圖2 冪指數(shù)粘聚區(qū)模型本構(gòu)示意圖Fig.2 Power law CZM constitution diagram

2 粘聚區(qū)本構(gòu)的實(shí)驗(yàn)測(cè)定

粘聚區(qū)本構(gòu)模型中包含3個(gè)重要參數(shù):粘聚區(qū)斷裂能、粘聚區(qū)斷裂強(qiáng)度和粘聚區(qū)函數(shù)形式。國(guó)外有相關(guān)研究學(xué)者認(rèn)為，粘聚區(qū)斷裂能等于材料的臨界擴(kuò)展J積分［8］。因此，本文使用推進(jìn)劑的臨界擴(kuò)展J積分作為粘聚區(qū)斷裂能。Begley和Landes［9］提出了多試樣法，通過(guò)一組初始裂紋長(zhǎng)度相近的試件來(lái)測(cè)量材料的J積分。

多試樣法需多個(gè)試樣進(jìn)行測(cè)量，實(shí)驗(yàn)量較大，實(shí)驗(yàn)過(guò)程較為繁瑣。因此，有學(xué)者提出了單試樣法的測(cè)量方法。Sumpter［10］提出了η因子的概念求解材料的J積分，并得到了廣泛應(yīng)用。

式中 U為載荷位移曲線積分;B為試樣厚度;W為試樣寬度;a為初始裂紋長(zhǎng)度;A為裂紋體初始斷裂韌帶面積。

本文將使用單試樣法來(lái)測(cè)定粘聚區(qū)斷裂能，但單試樣法公式(4)中的η因子需使用多試樣法來(lái)進(jìn)行標(biāo)定［11］。

對(duì)比式(3)、式(4)，得η因子的公式為

η因子的標(biāo)定在20℃下展開(kāi)，使用QJ 211B萬(wàn)能材料試驗(yàn)機(jī)在20 mm/min的拉伸速率進(jìn)行拉伸實(shí)驗(yàn)，試樣幾何尺寸B=5 mm、W=30 mm、裂紋初始長(zhǎng)度a不同。實(shí)驗(yàn)過(guò)程中記錄載荷位移曲線，同時(shí)使用CCD同步測(cè)量裂紋開(kāi)裂過(guò)程。使用式(5)對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行處理，可獲得η因子隨裂紋初始長(zhǎng)度的關(guān)系，如圖3示。從圖3可發(fā)現(xiàn)，η隨裂紋初始尺寸稍有變化，平均值變化范圍包含在測(cè)量標(biāo)準(zhǔn)差之內(nèi)。因此，可認(rèn)為η在所取的裂紋初始尺寸范圍內(nèi)為恒定值。

圖3 η因子與裂紋初始尺寸關(guān)系Fig.3 η factor vs initial crack size

J積分的測(cè)量使用單試樣法在20℃下展開(kāi)，試樣的初始裂紋長(zhǎng)度a/W=0.5，共進(jìn)行5次重復(fù)性實(shí)驗(yàn)。實(shí)驗(yàn)過(guò)程中采集試樣的載荷-位移曲線。對(duì)載荷-位移曲線積分，獲得裂紋試樣起裂時(shí)的裂尖輸入能量U，結(jié)合測(cè)量的η因子代入式(4)，得到HTPB推進(jìn)劑在20 mm/min的拉伸速度下的粘聚區(qū)斷裂能為(1.007±0.050)N/mm。粘聚區(qū)斷裂強(qiáng)度的測(cè)量采用標(biāo)準(zhǔn)單軸拉伸實(shí)驗(yàn)獲取。在20℃環(huán)境溫度下，進(jìn)行一組20 mm/min的等速拉伸試驗(yàn)，得到HTPB推進(jìn)劑的極限斷裂強(qiáng)度為(0.548±0.033)MPa。粘聚區(qū)損傷函數(shù)的最終軟化參數(shù)在第4章詳述。

3 數(shù)值算法

本文使用 ABAQUS的用戶自定義單元子程序UEL開(kāi)發(fā)出粘聚區(qū)單元，用來(lái)模擬HTPB推進(jìn)劑的起裂和開(kāi)裂過(guò)程。用ABAQUS Standard求解時(shí)，需提供粘結(jié)單元的切線剛度矩陣和力矢量，下面將詳細(xì)推導(dǎo)其具體形式。

根據(jù)有限元虛功原理，不考慮體力情況下，包含粘結(jié)單元的平衡方程和邊界的等效弱積分形式為

式中 δεij、δi、δui分別為實(shí)體單元應(yīng)變、粘結(jié)單元相對(duì)位移、邊界位移;σij、Tic、Tie分別代表實(shí)體單元應(yīng)力、粘結(jié)單元粘聚力、力邊界;dV、dSc、dSe分別代表實(shí)體單元控制體、粘結(jié)單元表面、外邊界。

式(6)中的第二項(xiàng)代表粘結(jié)單元的貢獻(xiàn)部分。圖4為粘結(jié)單元變形示意圖。

圖4 粘結(jié)單元Fig.4 Cohesive element

粘結(jié)單元相對(duì)位移為

其中:

式中 a為粘結(jié)單元節(jié)點(diǎn)在系統(tǒng)坐標(biāo)下的坐標(biāo)值;R為坐標(biāo)轉(zhuǎn)換矩陣。

插值形函數(shù)為

粘結(jié)單元的切線剛度矩陣可表示為

式中 C為粘結(jié)單元的Jacobian矩陣。

粘結(jié)單元所提供的力矢量表示為

ABAQUS開(kāi)發(fā)過(guò)程中需提供式(9)、式(10)，系統(tǒng)建立起總體剛度矩陣，按Newton迭代法進(jìn)行迭代求解。式(11)為系統(tǒng)剛度矩陣，式(12)為牛頓迭代法的計(jì)算格式。

4 仿真和模型驗(yàn)證

為了確定所使用的冪指數(shù)形式的粘聚區(qū)模型中的本構(gòu)參數(shù)α，使用所建立的粘聚區(qū)有限元仿真方法，對(duì)包含I型裂紋的HTPB推進(jìn)劑試樣進(jìn)行數(shù)值仿真。仿真模型為包含I型裂紋的板條形推進(jìn)劑試樣。試樣初始裂紋長(zhǎng)度為15 mm，上下表面施加20 mm/min等速位移載荷。HTPB推進(jìn)劑是一種典型的粘彈性材料，其力學(xué)特性和加載歷史存在相關(guān)性。在仿真過(guò)程中，使用線粘彈性本構(gòu)來(lái)描述HTPB推進(jìn)劑的力學(xué)行為。圖5中，灰色區(qū)域?yàn)槔鞂?shí)驗(yàn)獲得的載荷-時(shí)間變化范圍。圖5中，3條曲線分別為不同粘聚區(qū)本構(gòu)參數(shù)α下通過(guò)仿真獲得的載荷-時(shí)間曲線。從圖5可發(fā)現(xiàn)，隨著α的減小，仿真的最大載荷逐漸增大，斷裂時(shí)間后延，但曲線之間的上升段初期基本重合，且變化趨勢(shì)基本一致。通過(guò)對(duì)比發(fā)現(xiàn)，當(dāng)α=0.55時(shí)，仿真曲線基本位于實(shí)驗(yàn)結(jié)果的中間部分。因此，可確定出HTPB推進(jìn)劑粘聚區(qū)本構(gòu)參數(shù)α≈0.55。

圖5 單邊裂紋載荷-時(shí)間曲線仿真值和實(shí)驗(yàn)值Fig.5 Numerical and experimental load-time curve of signal edge notched sample

為了驗(yàn)證所建立的HTPB推進(jìn)劑粘聚區(qū)模型的準(zhǔn)確性和所建立的有限元仿真方法的精確性，下面將數(shù)值仿真預(yù)測(cè)結(jié)果和實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。如圖6所示，仿真模型左右兩側(cè)開(kāi)有10 mm初始裂紋，左右裂紋上下相距20 mm。由于裂紋擴(kuò)展路徑未知，需在所有網(wǎng)格單元之間嵌入粘結(jié)單元，本文使用MATLAB編程，生成了ABAQUS識(shí)別的網(wǎng)格文件。圖6為27 s時(shí)的仿真和實(shí)驗(yàn)過(guò)程中的斷裂形貌?？砂l(fā)現(xiàn)，2條初始裂紋已出現(xiàn)擴(kuò)展，試件的整體變形、裂紋擴(kuò)展形貌基本一致。圖6(a)給出了仿真中的Mises應(yīng)力分布情況?？煽吹剑鸭y尖端應(yīng)力分布上下不對(duì)稱，此時(shí)裂紋呈現(xiàn)出平面Ⅰ-Ⅱ型復(fù)合裂紋。

圖6 t=27 s斷裂形貌Fig.6 Fracture morphology at 27 s

圖7(a)為仿真結(jié)束后的網(wǎng)格，裂紋沿初始網(wǎng)格界面進(jìn)行擴(kuò)展。由于網(wǎng)格布置存在一定的不對(duì)稱性和計(jì)算中的舍入誤差原因，仿真結(jié)果中右上部裂紋貫穿，左下角裂紋未貫穿。圖7(b)為實(shí)驗(yàn)結(jié)束后裂紋擴(kuò)展路徑。由于裂紋體的不對(duì)稱性，所做的3次重復(fù)性實(shí)驗(yàn)均呈現(xiàn)出如圖7所示的單條裂紋貫穿情況，貫穿裂紋在擴(kuò)展過(guò)程中逐漸轉(zhuǎn)向。通過(guò)對(duì)比可發(fā)現(xiàn)，仿真預(yù)測(cè)的裂紋擴(kuò)展路徑和實(shí)驗(yàn)結(jié)果吻合良好，證明所建立的嵌入粘結(jié)單元的仿真方法，可較好地模擬出復(fù)合裂紋的擴(kuò)展路徑。本文仿真過(guò)程中，裂紋表面在開(kāi)裂后處于自由表面狀態(tài)。在真實(shí)狀況下，發(fā)動(dòng)機(jī)藥柱在裂紋擴(kuò)展過(guò)程中，伴隨著燃?xì)獾母Z入過(guò)程。因此，開(kāi)裂后的裂紋面將處于壓力載荷作用之下。如果需考慮到裂紋面內(nèi)存在壓力載荷的情況，需在式(10)所示的單元力矢量列陣中，加入外界壓力載荷的作用項(xiàng)。

圖7 仿真和實(shí)驗(yàn)裂紋開(kāi)裂路徑Fig.7 Numerical and experimental crack path

圖8為仿真和實(shí)驗(yàn)獲得的載荷時(shí)間曲線。圖8中，灰色區(qū)域?yàn)閷?shí)驗(yàn)測(cè)得的載荷時(shí)間變化情況，曲線為仿真預(yù)測(cè)結(jié)果。從圖8中發(fā)現(xiàn)，仿真曲線的最高大值與實(shí)驗(yàn)測(cè)量結(jié)果基本一致。曲線前半部分在實(shí)驗(yàn)測(cè)量范圍之內(nèi)，說(shuō)明所建立的模型對(duì)裂紋的起裂預(yù)測(cè)較為準(zhǔn)確。曲線后半部分下降段曲線有所提前，這可能是因?yàn)閷?shí)際情況下推進(jìn)劑裂尖局部區(qū)域外也產(chǎn)生了脫濕等損傷，造成了材料整體軟化、模量下降。因此，實(shí)際載荷時(shí)間曲線比預(yù)測(cè)值后移。

圖8 驗(yàn)證試樣仿真和實(shí)驗(yàn)載荷-時(shí)間曲線Fig.8 Numerical and experimental load-time curve of verifiable sample

綜合圖7、圖8可發(fā)現(xiàn)，所建立的粘聚區(qū)本構(gòu)模型和仿真方法，為HTPB推進(jìn)劑裂紋起裂和裂紋擴(kuò)展過(guò)程提供了一個(gè)較先進(jìn)的物理模型和數(shù)值仿真方法。傳統(tǒng)有限元計(jì)算方法在計(jì)算裂紋起裂問(wèn)題時(shí)，由于裂尖應(yīng)力奇異性造成了裂尖應(yīng)力計(jì)算的不可信。使用粘聚區(qū)模型來(lái)模擬裂尖應(yīng)力時(shí)，粘聚區(qū)本構(gòu)模型的損傷特性避免了裂尖應(yīng)力的奇異性，使裂尖應(yīng)力更加合乎物理實(shí)際。裂紋起裂擴(kuò)展后，造成了物理模型在裂尖的幾何不連續(xù)性，這是傳統(tǒng)有限元單元所無(wú)法實(shí)現(xiàn)的。本文建立的仿真計(jì)算方法，通過(guò)合理地構(gòu)建粘結(jié)單元來(lái)模擬裂紋開(kāi)裂這一物理現(xiàn)象，可有效地解決由于裂紋擴(kuò)展所造成的幾何不連續(xù)性。

實(shí)驗(yàn)中發(fā)現(xiàn)，HTPB推進(jìn)劑斷裂性能和斷裂強(qiáng)度隨加載速率變化，HTPB推進(jìn)劑粘聚區(qū)本構(gòu)也存在一定的粘彈性效用。因此，建立相關(guān)的HTPB推進(jìn)劑粘聚區(qū)斷裂模型將是下一步的研究重點(diǎn)。

5 結(jié)論

(1)建立的粘聚區(qū)斷裂模型能很好地反映出HTPB推進(jìn)劑裂紋起裂和裂紋擴(kuò)展的失效過(guò)程。

(2)基于ABAQUS建立了模擬裂紋擴(kuò)展的粘結(jié)單元，并建立了模擬復(fù)合型裂紋開(kāi)裂過(guò)程的嵌入粘結(jié)單元法，使用該方法成功預(yù)測(cè)了HTPB推進(jìn)劑裂紋擴(kuò)展路徑。

［1］李九天，雷勇軍，唐國(guó)金，等.固體火箭發(fā)動(dòng)機(jī)藥柱表面裂紋分析［J］.固體火箭技術(shù)，2008，31(5):471-474.

［2］Earnest T E.RSRM TP-H1148 main grain propellant crack initiation evaluation［R］.AIAA 2005-3061.

［3］Barenblatt G I.The formation of equilibrium cracks during brittle fracture:general ideas and hypotheses axially-symmetric cracks［J］.Journal of Applied Mathematics and Mechanics，1959，23(3):622-636.

［4］Dugdale D.Yielding of steel sheets containing slits［J］.Journal of the Mechanics and Physics of Solids，1960，8(2):100-104.

［5］Xu X，Needleman A.Numerical simulation of fast crack growth in brittle solids［J］.Journal of the Mechanics and Physics of Solids，1994，42(9):1397-1434.

［6］劉陸廣，歐卓成，段卓平，等.混凝土動(dòng)態(tài)斷裂數(shù)值模擬［J］.兵工學(xué)報(bào)，2010，31(6):741-745.

［7］崔浩，李玉龍，劉元鏞，等.基于粘聚區(qū)模型的含填充區(qū)復(fù)合材料接頭失效數(shù)值模擬［J］.復(fù)合材料學(xué)報(bào)，2010，27(2):161-168.

［8］Tabakovic A，Karac A，Ivankovic A，et al.Modelling the quasi-static behaviour of bituminous material using a cohesive zone model［J］.Engineering Fracture Mechanics，2010，77(13):2403-2418.

［9］Begley J A，Landes J D.The J integral as a fracture criterion［R］.ASTM STP 514，1972.

［10］Sumpter J D G.Elastic plastic fracture analysis and design using the finite element method［D］.London:University of London，1973.

［11］Ramorino G，Agnelli S，De R Santis，et al.Investigation of fracture resistance of natural rubber/clay nanocomposites by J-testing［J］.Engineering Fracture Mechanics，2010，77(10):1527-1536.