基于內聚力法則的高能硝酸酯增塑聚醚推進劑開裂過程細觀模型

侯宇菲， 許進升， 古勇軍， 周長省

(1.南京理工大學 機械工程學院, 江蘇 南京 210094; 2.山西北方興安化學工業有限公司, 山西 太原 030008)

0 引言

固體推進劑作為固體火箭發動機的能量來源，其結構完整性直接決定整體火箭系統及作戰平臺的安全。為保證高能的要求，固體推進劑中摻入了大量固體填料，因而形成了一種多相體(基體相、增強相、界面相)含能材料。這種多相材料在載荷作用下，將會引起顆粒與基體界面脫濕、裂紋擴展、變形等缺陷，其宏觀力學性能由細觀多尺度物理所決定[1-2]。因此，從細觀角度出發研究細觀結構對推進劑力學性能的影響，對于推進劑配方的研制具有科學的指導意義[3]。

近年來，國內外學者采用掃描電鏡與數字圖像相關方法對高能材料的細觀結構演化進行觀察，發現顆粒與基體的界面是多相材料失效的起點[4-5]。然而，由于實驗手段的局限性，無法定量分析高能材料的細觀損傷破壞演化過程，因此數值仿真法成為近年來興起的分析多相材料界面破壞的有效方法。陳煜等[6]采用單顆粒軸對稱細觀模型對硝酸酯增塑聚醚(NEPE)推進劑在不同固含量、鍵合劑等作用條件下的應力集中因子及脫濕角進行了研究。Matous等[7]利用分子動力學方法生成推進劑細觀代表性體積單元，研究了在小應變條件下顆粒界面的損傷與失效，證明了界面損傷是造成力學性能非線性的主要原因。Inglis[8]在Matou?的研究基礎上，采用橢圓形顆粒進行數值計算，得知模型損傷演化過程分為3個階段：彈性階段、顆粒與基體界面脫濕和基體大變形。Cui等[9]通過實驗得知推進劑的力學性能與加載速率相關，為更準確地描述推進劑在工作過程中的力學性能，建立了帕克- 保利諾- 勒斯勒爾(PPR)率相關粘聚區準則。Zhang等[10]通過不同溫度及應變速率的單軸拉伸實驗發現，只考慮推進劑顆粒與基體界面脫濕的數值仿真無法準確地預測推進劑的失效行為，以此推測推進劑的失效還包括基體的損傷。周水平等[11]采用原位掃描電鏡得出推進劑的失效包括顆粒與基體界面脫濕與基體拉絲破壞。但以上方法都局限在顆粒與基體界面脫濕對推進劑力學性能的影響，尚未考慮基體失效對推進劑力學性能的影響。

本文針對復合固體推進劑細觀損傷演化規律，建立同時考慮顆粒與基體界面脫濕與基體開裂的細觀顆粒填充模型，并采用VUMAT子程序編譯出多項式- 梯形內聚力法則，分別與雙線性、多項式內聚力法則進行對比研究，得出更符合高能NEPE推進劑開裂過程的細觀顆粒模型。

1 NEPE推進劑細觀顆粒模型

1.1 細觀顆粒填充模型

本文以NEPE推進劑為研究對象，其本構關系通過準靜態單軸拉伸實驗獲取，單軸拉伸實驗在QJ萬能材料試驗機上進行。該NEPE推進劑的主要成分包括鋁(Al)粉、高氯酸銨(AP)及其他輔助制劑，各組分含量如表1所示。通過分子動力學算法[12]生成與復合固體推進劑具有相同體積分數的細觀顆粒填充模型，其中Al粉粒徑為8～20 μm，AP(1)的粒徑為100～120 μm，AP(2)的粒徑為250～300 μm，粒徑分布如圖1所示。

表1 NEPE推進劑基本組分Tab.1 Components of NEPE propellant

圖1 NEPE固體推進劑顆粒粒徑分布Fig.1 Particle size distribution of NEPE propellant

圖2 NEPE固體推進劑顆粒填充幾何模型Fig.2 Particle packing model of NEPE propellant

根據文獻[13]得知，當代表性體積單元的尺寸為最大粒徑尺寸的3～5倍時，則認為代表性體積單元有效。因此，本文根據圖1填充顆粒粒徑分布規律得到推進劑細觀顆粒模型的代表性體積單元，如圖2所示，其尺寸為1 000 μm×1 000 μm. 由于Al粉顆粒的粒徑小且數量多，若直接對該模型進行網格劃分，將會降低網格質量并影響計算結果的收斂性。文獻[14]通過實驗證明，Al顆粒作為加速復合推進劑熱分解的添加劑，不會引起顆粒與基體界面脫濕。因此采用Mori-Tanaka解析法[15]將Al顆粒對推進劑力學性能的影響等效到基體中，在計算時只考慮AP顆粒的影響。材料參數如表2所示。

表2 細觀材料力學屬性Tab.2 Mechanical properties of meso-components

1.2 粘結單元

粘結單元的基本思想是構造一種可以賦予材料損傷特性的特殊單元，其力學行為遵循內聚力法則的力學響應，當粘結單元與有限元模型結合后，實現模型的失效破壞。本文利用Python腳本語言在顆粒與基體界面及基體內部嵌入零厚度粘結單元，實現NEPE推進劑顆粒與基體界面脫濕至基體斷裂的過程。具體如下：

1) 將顆粒填充模型導入Abaqus軟件中，分別建立基體單元集合與顆粒單元集合，采用四節點的平面應變單元對以上集合進行網格劃分，其網格類型為CPE4，并以Abaqus軟件inp格式輸出。

2) 遍歷整個模型集合，劃分顆粒單元與基體單元的公共節點及基體單元內部公共節點，并建立新的單元集合。如圖3(a)所示，深色部分為AP顆粒單元，淺色部分為基體單元。

圖3 粘結單元插入示意圖Fig.3 Schematic diagram of embedded cohesive elements

3) 將所有公共節點拆分并重新編號，保留新舊節點坐標對應關系，如圖3(b)所示，節點4被拆分為節點4和節點13，節點6被拆分為節點7和節點10.

4) 基于新舊節點坐標對應關系，對顆粒單元和基體單元重新編號并生成粘結單元，粘結單元編號時必須按逆時針方向編號，其網格類型為COH2D4. 粘結單元分布如圖3(b)所示，在實際計算過程中粘結單元的幾何厚度為0 mm.

5) 生成新的inp文件，為方便后續添加材料屬性，建立4個集合放置4種不同材料屬性的單元：AP顆粒單元集合，基體單元集合，界面單元集合，基體粘結單元集合。

1.3 邊界條件

為模擬NEPE推進劑在單軸等速拉伸下的細觀失效過程，需對所建立的細觀顆粒填充模型選用合理的邊界條件，為滿足細觀模型周期性及材料均勻性假設，通常對代表性體積單元施加周期性邊界條件。然而，周期性邊界條件在實現過程中依賴于周期性網格，不適用于 NEPE推進劑的復雜細觀結構。文獻[16]表明，只要代表性體積單元的尺寸適當，均勻位移邊界條件可替代周期性邊界條件。因此本文選用均勻位移邊界條件，具體施加方式如圖4所示。

圖4 NEPE推進劑細觀模型邊界條件施加Fig.4 Boundary conditions of mesoscopic model of NEPE propellant

2 多項式- 梯形內聚力法則的建立

2.1 雙線性內聚力法則

20世紀60年代，Barenblatt[17]和Dugdale[18]針對斷裂力學中裂紋尖端的應力奇異性問題，首次提出了內聚力法則，內聚力模型如圖5所示。內聚力區域的力學狀態通過最大斷裂強度σmax、界面臨界張開位移δ0和斷裂能φ確定[19]。其中，雙線性內聚力法則最為常見，如(1)式所示：

(1)

式中：T為斷裂強度；σ為張開位移；δf為特征位移。

圖5 內聚力模型示意圖Fig.5 Schematic diagram of real fracture process of cohesive zone model

最大斷裂強度σmax與界面臨界張開位移δ0的比值為模量K. 先前研究者認為內聚力法則中的參數最關鍵，而曲線形式是次要的。但是，文獻[20]在研究鉻與鋯鈦酸鉛(Cr/PZT)界面開裂時發現，雙線性內聚力法則能夠與實驗結果匹配，而指數型的仿真結果較差。

2.2 多項式- 梯形內聚力法則

針對NEPE推進劑復雜的力學行為，所選用的內聚力法則必須同時滿足高斷裂能和張力位移非線性損傷特點。對于雙線性內聚力法則而言，若想增加斷裂能，則必須提高斷裂強度或臨界開裂位移，而斷裂強度或臨界開裂位移的增加會引起應力值增大或界面無法開裂的現象。為解決張力位移非線性問題，Tvergaard等[21]在研究彈塑性材料破壞時提出梯形內聚力法則，這種內聚力法則有效地提高了斷裂能，更好地描述了材料的塑性損傷階段。然而梯形內聚力法則的上升段、下降段均為線性函數，并不適用于表征NEPE推進劑材料非線性張力位移關系(通常只有5%的彈性階段)。

Needleman等[22]為解決張力位移非線性問題，建立了多項式內聚力法則，與雙線性內聚力法則、梯形內聚力法則不同，其張力位移關系為連續性方程，需要通過斷裂能控制方程求偏導得到。二維平面狀態下的多項式內聚力模型斷裂能方程為

(2)

式中：un和ut為法向與切向的臨界界面位移分離值；α為法向剛度與切向剛度的比例系數。

(2)式求導，可得

(3)

(4)

式中：Tn和Tt分別為法向和切向的斷裂強度。

為同時滿足推進劑材料的延展性損傷特點和張力位移非線性變化關系，本文建立了多項式- 梯形內聚力法則，將多項式與梯形內聚力法則用分段函數表示，應力的上升段和下降段為多項式函數，應力屈服段為常數，如圖6所示，并通過VUMAT子程序完成了該法則的數值開發，可更準確地描述NEPE推進劑的力學性能。

圖6 梯形- 多項式內聚力法則Fig.6 Trapezoidal-polynomial cohesive law

多項式- 梯形內聚力法則的位移與張力控制方程為

(5)

(6)

斷裂能φ控制方程為

(7)

(8)

由(8)式計算各向應力σ：

(9)

式中：T0為粘結單元計算厚度。

subroutine vumat()

Input:

Smax !斷裂強度

Disnn !延性損傷位移

Distn !特征位移

Do 100i=1, nblock

Compute(stress, strain, displacement, energy)

Improvement(stress, strain, displacement, energy)

If (stateNew(error)==0) then

stateNew(delete)=0

else

stateNew(delete)=1

end if

100 continue

Return

End

圖7 多項式- 梯形內聚力法則算法

Fig.7 Trapezoidal-polynomial cohesive law algorithm

為實現粘結單元的損傷演化至失效過程，通過設置狀態變量SDV(stateNew)判斷粘結單元的失效，當狀態變量由1變為0時，粘結單元進行不可逆刪除。

2.3 參數反演識別

對于內聚力法則參數獲取方法，現階段很難直接由實驗獲得。本文采用反演分析法獲得內聚法則的材料參數。以實驗曲線上數據點的位移值為基準，在仿真曲線上找到最靠近該基準點的最小區間，隨后在最小區間的兩個端點間線性插值得到仿真值，最后利用目標函數得到實驗值與仿真值的最小誤差R(具體步驟見文獻[23])：

(10)

本文分別對雙線性內聚力法則、多項式內聚力法則和多項式- 梯形內聚力法則進行參數反演識別，并假設推進劑為各向同性材料，內聚力法則的法向參數與切向參數相同，反演后的參數如表3～表5所示。

表3 雙線性內聚力模型參數Tab.3 Parameters of bilinear cohesive zone model

表4 多項式內聚力模型參數Tab.4 Parameters of polynomial cohesive zone model

表5 多項式- 梯形內聚力模型參數Tab.5 Parameters of polynomial-trapezoidal zone model

3 推進劑開裂過程細觀仿真及失效分析

3.1 網格收斂性

為消除有限元模型的物理離散誤差，有必要對其進行網格收斂性檢查。通過減小單元尺寸即提高插值函數的階次，使有限元模型的解收斂于精確解。同時為確保粗網格的計算結果是定性正確的，粗網格的尺寸Δy不宜過大，并且不同網格尺寸Δy的比值c保持一致，即

(11)

圖8給出了4種網格尺寸Δy的有限元模型。下面對顆粒與基體界面為雙線性內聚力法則的細觀顆粒填充模型施加微小載荷，采用有限元計算獲得初始彈性模量，與實驗彈性模量作對比，獲得最精準的網格尺寸。

圖8 NEPE推進劑網格示意圖Fig.8 Schematic diagram of NEPE propellant mesh

如圖9所示，隨著網格尺寸的減小，有限元模型的初始彈性模量逐漸趨于穩定并接近于實驗值，表明實體單元階次對計算結果的影響隨網格尺寸減小而逐漸降低。當網格尺寸Δy=0.02 mm時，初始模量E0的計算結果收斂，因此采用該網格尺寸作為本文相關模型的網格尺寸標準。

3.2 推進劑脫濕細觀仿真及分析

為證明NEPE推進劑細觀應力場的分布與內聚力法則的曲線形式有關，下面采用不同曲線形式的內聚力模型對NEPE推進劑細觀顆粒填充模型的損傷演化過程進行分析。其中：模型1為顆粒與基體界面采用雙線性內聚力法則；模型2為顆粒與基體界面采用多項式內聚力法則；模型3顆粒與基體界面采用多項式- 梯形內聚力法則。

表6分別給出了模型1、模型2、模型3應變為8%、24%、34%、50%時的應力云圖及損傷演化情況。當應變為8%時，顆粒與基體之間未出現明顯的位移分離，界面粘結性能良好。但由于顆粒與基體剛度嚴重不匹配，導致模型內部應力分布極不均勻，顆粒應力明顯高于基體應力，表明模型的整體應力強度因AP顆粒而提高。當應變增加到24%時，大顆粒與基體之間發生位移分離，這是因為顆粒剛度遠大于基體剛度，隨著應變的增加，基體的變形遠大于顆粒的變形，導致顆粒與基體界面成為整個模型中最薄弱的環節，即顆粒與基體之間出現位移分離的現象。當應變為34%時，基體變形增大，顆粒與基體的位移分離量增加，導致顆粒從孔洞中裸露出來進而形成脫濕。當應變達到50%時，模型中顆粒與基體界面脫濕數量持續增加，直至所有顆粒與基體的界面都消失。

由表6可知，3個模型的損傷演化過程相同，但應力值分布不同，表明模型應力分布和內聚力法則的曲線形式有關，但單從應力云圖無法判斷哪一種內聚力法則更適合描述NEPE推進劑的顆粒與基體界面失效過程。

圖9 NEPE推進劑網格- 初始彈性模量曲線Fig.9 NEPE propellant mesh-initial elastic modulus

3.3 推進劑脫濕損傷分析

圖10所示為實驗應力- 應變曲線。從圖10中可以發現，NEPE推進劑的力學特性曲線可分為4個階段：線性段、軟化段、下降段和斷裂破壞段。在線性段，應力隨應變的增加而線性增加，此時推進劑性能保持完整，尚未發生損傷。當到達軟化段時，應力- 應變曲線出現拐點(即A點處)，曲線斜率緩慢降低，材料的模量開始下降，結合應力云圖發現此時顆粒與基體界面發生了位移分離。在下降段，應力隨應變的增加而快速降低，曲線的斜率進一步降低。在斷裂破壞段，應力隨應變增加而緩慢增加，當跨過應力平臺區后，應力快速下降，形成斷裂。

表6 NEPE推進劑Mises應力云圖

圖10 NEPE推進劑應力- 應變曲線Fig.10 Stress-strain curves of NEPE propellant

對比預測結果與實驗結果得知，在應力- 應變曲線的線性段，3種模型的仿真曲線與實驗曲線完全擬合。當到達軟化段時，模型2與模型3預測結果較好，模型1的應力值偏高。這是因為雙線性內聚力法則與多項式內聚力法則的損傷變量形式不同，在雙線性內聚力法則中，只有當材料點的應力值到達斷裂強度時模型才開始損傷，而多項式內聚力法則的損傷是從載荷施加的起點開始迭代，因此多項式內聚力法則表征的顆粒與基體界面剛度軟化較快。當3種模型處于下降段和斷裂破壞段時，其預測結果均偏離實驗結果，模型1和模型2的應力隨著應變的增加而增加，模型Ⅲ經過一個微小的應力平臺區后，其應力值也隨即增大，最終3種模型的應力值均為發生下降。由此可知，顆粒與基體界面脫濕只能造成推進劑模量的下降，并不會引起推進劑的失效破壞。因此，只考慮顆粒與基體界面損傷的細觀顆粒模型并不能完整地預測推進劑的力學性能。

3.4 推進劑開裂過程細觀仿真及分析

通過3.2節得知，僅考慮顆粒與基體界面脫濕的細觀顆粒模型并不能完整地預測NEPE推進劑的宏觀力學性能，因此下面采用考慮基體失效的細觀顆粒填充模型進行數值仿真，其中： 模型4為顆粒與基體界面及基體內部采用雙線性內聚力法則；模型5為顆粒與基體界面及基體內部采用多項式內聚力法則；模型6為顆粒與基體界面及基體內部采用多項式- 梯形內聚力法則；模型7為顆粒與基體界面采用多項式內聚力法則，基體內部采用多項式- 梯形內聚力法則。

圖11 應變分別為8%、24%、34%、50% 時NEPE 推進劑Mises應力云圖Fig.11 Mises stress nephograms of NEPE propellant at the strains of 8%, 24%, 34%, and 50%

由于模型數量較多，并且各個模型的云圖變化過程相似，下面只給出模型7在應變為 8%、24%、34%、50%時的應力云圖，如圖11所示。當應變為8%、24%、34%時，其仿真結果與3.2節的仿真結果相同，表現為未損段、顆粒與基體產生位移分離段和顆粒與基體界面脫濕段。然而，當應變達到50%時，未考慮基體失效與考慮基體失效的細觀顆粒填充模型其應力云圖產生明顯區別，如圖10與圖11所示。考慮基體失效的細觀顆粒填充模型其顆粒與基體界面脫濕數量維持水平，部分孔洞發生匯聚，導致基體內部形成微裂紋。而未考慮基體失效的細觀顆粒填充模型，其顆粒與基體界面脫濕數量持續增加，模型內部未出現微裂紋。通過文獻[10]得知，考慮基體失效的細觀顆粒填充模型其預測結果和實際情況相符，并且裂紋方向與載荷加載方向相垂直。這是因為在未考慮基體損傷的模型中，雖然顆粒與基體界面脫濕導致基體內部形成高應力區，但由于基體未考慮損傷，高應力區域無法形成裂紋，模型整體即可持續拉伸，直至所有顆粒與基體界面都脫濕。而考慮基體損傷的模型中，當顆粒與基體界面脫濕后，基體成為載荷承擔主體，在承載的初始階段應力隨著應變增加而小幅度增加，但由于顆粒與基體界面脫濕后產生了大量的孔洞，在與載荷方向相垂直的孔洞處形成了高應力區，隨著載荷的持續加載，高應力區迅速形成裂紋。當裂紋產生后未脫濕區的應力值轉而降低，因此顆粒與基體界面脫濕數量維持恒定。

3.5 推進劑斷裂損傷分析

通過上述分析得知，顆粒與基體界面脫濕引起基體內部產生高應力區是造成推進劑破壞的關鍵，因此基體損傷不容忽略。與3.3節的結果相似，以上4個模型的損傷演化過程相同，但應力值分布不同，因此下面通過細觀模型預測應力- 應變曲線與實驗應力- 應變曲線來驗證模型的準確性。

如圖12所示，考慮基體失效的細觀顆粒填充模型可以完整地模擬出NEPE推進劑力學曲線變化的4個階段。在線性段和軟化段，采用多項式內聚力法則的細觀顆粒填充模型，其預測結果與實驗結果更為接近。在下降段，模型6、模型7較模型4、模型5的應力值下降更為緩慢，更接近實驗結果。在斷裂破壞段，模型4和模型5的應力隨應變的增加而增加，當到達應力峰值后快速斷裂，期間并未出現應力平臺區；而模型6和模型7的曲線走勢基本和實驗結果相一致，當跨越一個緩慢上升的應力平臺區后快速斷裂。

圖12 NEPE推進劑應力- 應變曲線Fig.12 Stress-strain curves of NEPE propellant

通過上述分析得知，模型6和模型7的仿真結果與實驗結果最為接近。這是因為：1)多項式內聚力法則中的張力位移為非線性變化，更適合描述NEPE推進劑顆粒與基體界面脫濕的實際損傷過程；2)NEPE推進劑基體材料變形較大，采用多項式- 梯形內聚力法則可以有效地模擬出推進劑在拉伸過程中的應力平臺區，并提高模型的斷裂能，而雙線性內聚力法則更適合描述脆性材料的損傷，無法體現推進劑材料延展性損傷的特點。

在圖12中，模型6和模型7的仿真結果略高于實驗結果，分析誤差產生的原因，有以下3種：1)本文推進劑的細觀模型是二維的，網格的單元類型均基于平面應變，與真實推進劑的三維應力與應變狀態存在一定區別；2)目前，為簡化模型，將填充顆粒均設置為圓形顆粒，但在真實結構中顆粒的形狀是類似于圓形的；3)真實推進劑在固化攪拌過程中，顆粒與基體界面存在一定缺陷，而在本文模型中，認為顆粒與基體之間是完全粘結在一起的，不存在任何缺陷。這些原因均使得仿真結果偏高。

4 結論

本文基于Abaqus有限元軟件建立了NEPE推進劑細觀顆粒填充模型，分析了顆粒與基體界面脫濕與基體開裂對推進劑力學性能的影響。得出如下主要結論：

1) 本文建立的多項式- 梯形內聚力法則較多項式內聚力法則具有更高的斷裂能，較雙線性內聚力法則實現了延展性損傷特點，更符合NEPE推進劑材料的力學性能。

2) 通過顆粒與基體界面脫濕和基體開裂細觀模型的計算結果得知，顆粒與基體界面脫濕引起基體內部形成高應力區，高應力區在載荷作用下形成裂紋，進而導致推進劑開裂。因此，推進劑開裂由顆粒與基體界面脫濕和基體開裂共同造成。

3) 本文建立的同時考慮顆粒與基體界面脫濕和基體開裂細觀模型，其預測應力- 應變曲線與實驗應力- 應變曲線擬合較好，表明模型有效。