HTPB推進劑粘聚斷裂研究①

韓 波，鞠玉濤，周長省

(南京理工大學機械工程學院，南京 210094)

0 引言

固體推進劑藥柱是固體火箭發動機的重要部件，藥柱結構完整性分析是保證固體火箭發動機正常內彈道性能的必要條件。藥柱表面微小裂紋的存在，對發動機能否正常工作有重大影響。因此，研究推進劑藥柱表面微小裂紋的起裂、擴展規律具有十分重大現實意義。

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

1 粘聚區模型

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

圖1為粘聚區單元模型示意圖，(a)為常規有限元實體單元，單元之間共節點，無法直接模擬材料的宏觀開裂過程;(b)為添加了粘結單元的有限元網格，實體單元之間通過粘結單元連接，粘結單元來模擬材料裂尖的損傷和演化直至宏觀裂紋形成，裂紋的可能擴展路徑為實體單元之間的界面。

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

以平面Ⅰ-Ⅱ型復合裂紋為例，粘聚區單元有效位移和有效應力定義為

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

σe和δe之間的變化關系，即粘聚區的本構模型，它預示著材料裂尖的損傷和演化過程。圖2為冪指數粘聚區本構的示意圖。

圖2中，橫坐標代表材料裂尖損傷區變形量;縱坐標代表裂尖損傷應力。本構關系中，上升段是為了數值模擬的需要，代表著裂尖材料未損傷前的響應;下降段代表著材料的損傷軟化形式。

式(2)為冪指數形式的粘聚區具體形式，可通過調節參數α來模擬不同形式的軟化形式。

圖2 冪指數粘聚區模型本構示意圖Fig.2 Power law CZM constitution diagram

2 粘聚區本構的實驗測定

粘聚區本構模型中包含3個重要參數:粘聚區斷裂能、粘聚區斷裂強度和粘聚區函數形式。國外有相關研究學者認為，粘聚區斷裂能等于材料的臨界擴展J積分［8］。因此，本文使用推進劑的臨界擴展J積分作為粘聚區斷裂能。Begley和Landes［9］提出了多試樣法，通過一組初始裂紋長度相近的試件來測量材料的J積分。

多試樣法需多個試樣進行測量，實驗量較大，實驗過程較為繁瑣。因此，有學者提出了單試樣法的測量方法。Sumpter［10］提出了η因子的概念求解材料的J積分，并得到了廣泛應用。

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

本文將使用單試樣法來測定粘聚區斷裂能，但單試樣法公式(4)中的η因子需使用多試樣法來進行標定［11］。

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

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

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

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

3 數值算法

本文使用 ABAQUS的用戶自定義單元子程序UEL開發出粘聚區單元，用來模擬HTPB推進劑的起裂和開裂過程。用ABAQUS Standard求解時，需提供粘結單元的切線剛度矩陣和力矢量，下面將詳細推導其具體形式。

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

式中 δεij、δi、δui分別為實體單元應變、粘結單元相對位移、邊界位移;σij、Tic、Tie分別代表實體單元應力、粘結單元粘聚力、力邊界;dV、dSc、dSe分別代表實體單元控制體、粘結單元表面、外邊界。

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

圖4 粘結單元Fig.4 Cohesive element

粘結單元相對位移為

其中:

式中 a為粘結單元節點在系統坐標下的坐標值;R為坐標轉換矩陣。

插值形函數為

粘結單元的切線剛度矩陣可表示為

式中 C為粘結單元的Jacobian矩陣。

粘結單元所提供的力矢量表示為

ABAQUS開發過程中需提供式(9)、式(10)，系統建立起總體剛度矩陣，按Newton迭代法進行迭代求解。式(11)為系統剛度矩陣，式(12)為牛頓迭代法的計算格式。

4 仿真和模型驗證

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

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

為了驗證所建立的HTPB推進劑粘聚區模型的準確性和所建立的有限元仿真方法的精確性，下面將數值仿真預測結果和實驗結果進行對比。如圖6所示，仿真模型左右兩側開有10 mm初始裂紋，左右裂紋上下相距20 mm。由于裂紋擴展路徑未知，需在所有網格單元之間嵌入粘結單元，本文使用MATLAB編程，生成了ABAQUS識別的網格文件。圖6為27 s時的仿真和實驗過程中的斷裂形貌。可發現，2條初始裂紋已出現擴展，試件的整體變形、裂紋擴展形貌基本一致。圖6(a)給出了仿真中的Mises應力分布情況。可看到，裂紋尖端應力分布上下不對稱，此時裂紋呈現出平面Ⅰ-Ⅱ型復合裂紋。

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

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

圖7 仿真和實驗裂紋開裂路徑Fig.7 Numerical and experimental crack path

圖8為仿真和實驗獲得的載荷時間曲線。圖8中，灰色區域為實驗測得的載荷時間變化情況，曲線為仿真預測結果。從圖8中發現，仿真曲線的最高大值與實驗測量結果基本一致。曲線前半部分在實驗測量范圍之內，說明所建立的模型對裂紋的起裂預測較為準確。曲線后半部分下降段曲線有所提前，這可能是因為實際情況下推進劑裂尖局部區域外也產生了脫濕等損傷，造成了材料整體軟化、模量下降。因此，實際載荷時間曲線比預測值后移。

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

綜合圖7、圖8可發現，所建立的粘聚區本構模型和仿真方法，為HTPB推進劑裂紋起裂和裂紋擴展過程提供了一個較先進的物理模型和數值仿真方法。傳統有限元計算方法在計算裂紋起裂問題時，由于裂尖應力奇異性造成了裂尖應力計算的不可信。使用粘聚區模型來模擬裂尖應力時，粘聚區本構模型的損傷特性避免了裂尖應力的奇異性，使裂尖應力更加合乎物理實際。裂紋起裂擴展后，造成了物理模型在裂尖的幾何不連續性，這是傳統有限元單元所無法實現的。本文建立的仿真計算方法，通過合理地構建粘結單元來模擬裂紋開裂這一物理現象，可有效地解決由于裂紋擴展所造成的幾何不連續性。

實驗中發現，HTPB推進劑斷裂性能和斷裂強度隨加載速率變化，HTPB推進劑粘聚區本構也存在一定的粘彈性效用。因此，建立相關的HTPB推進劑粘聚區斷裂模型將是下一步的研究重點。

5 結論

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

(2)基于ABAQUS建立了模擬裂紋擴展的粘結單元，并建立了模擬復合型裂紋開裂過程的嵌入粘結單元法，使用該方法成功預測了HTPB推進劑裂紋擴展路徑。

［1］李九天，雷勇軍，唐國金，等.固體火箭發動機藥柱表面裂紋分析［J］.固體火箭技術，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］劉陸廣，歐卓成，段卓平，等.混凝土動態斷裂數值模擬［J］.兵工學報，2010，31(6):741-745.

［7］崔浩，李玉龍，劉元鏞，等.基于粘聚區模型的含填充區復合材料接頭失效數值模擬［J］.復合材料學報，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.