考慮顆粒形狀的復合固體推進劑細觀損傷分析①

侯宇菲，許進升，陳 雄，周長省，朱 亮

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

0 引言

復合固體推進劑材料的力學性能在固體火箭發動機藥柱的結構完整性中起著重要作用。若想達到固體推進劑在工作過程中的力學指標，必須提前對推進劑的力學性能進行預測，保證固體推進劑的結構完整性。在過去的幾十年中，許多研究人員從連續介質力學角度對固體推進劑本構關系進行了大量的探索與研究[1]，這些理論研究均將固體推進劑視為一種連續均勻的介質。然而，復合固體推進劑是一種多顆粒填充的非均勻含能材料，上述方法不能從本質上解釋造成固體推進劑力學性能非線性的原因。因此，采用細觀方法對復合固體推進劑的力學性能進行研究成為一種有效途徑。

采用掃描電鏡對丁羥推進劑的端口形貌進行觀察，發現顆粒/粘合劑界面脫濕是造成推進劑力學行為非線性的主要原因[2]。由于現階段試驗條件的限制，高填充比顆粒復合材料的細觀實驗還難以開展，研究者通常采用數值仿真方法進行顆粒/基體界面脫濕的相關力學行為研究分析。彭威[3]使用軸對稱圓柱體模型進行數值模擬，發現應力集中現象出現在顆粒/基體界面的極區。Matous[4-5]通過自主開發的軟件生成了推進劑的代表性體積單元，研究了不同加載情況下推進劑的細觀損傷機理及過程。為進一步探究顆粒/基體界面的力學行為，常武軍[6]采用雙線性粘聚區模型對顆粒/基體界面脫濕進行數值模擬，得到了不同界面特性及顆粒尺寸對宏觀力學的影響規律。封濤[7]在顆粒/基體界面采用Surface-based Cohesive方法對推進劑顆粒/基體脫濕進行了數值模擬。袁崇[8]在研究顆粒/基體界面脫濕時，通過Mori-Tanaka方法得出填充顆粒粒徑、填充顆粒體積分數及顆粒界面間最大粘接應力對推進劑的力學性能均有明顯影響，并在外載荷作用下，由于顆粒與基體的材料屬性不同，推進劑內部的應力分布是不均勻的，當顆粒與基體界面發生脫濕時，推進劑內部的應力將重新分布。韓龍[9]通過推進劑的力學試驗與細觀數值模擬發現推進劑的力學性能與基體膠片的松弛行為有關，與顆粒的隨機分布無關。Zhi[10]著重探討了顆粒/基體界面脫濕對松弛模量的影響，結果表明推進劑損傷程度越高，其松弛模量越低。劉新國[11]分別從觀測及表征方法、理論模型及數值模擬三方面對顆粒/基體脫濕行為進行相關描述，認為解決顆粒/基體脫濕是提高推進劑宏觀力學性能的關鍵。

本文將從細觀角度出發，采用分子動力學方法建立圓形顆粒填充模型與多邊形顆粒填充模型，在顆粒與基體的界面處使用Python腳本語言嵌入零厚度Cohesive Element，分別使用雙線性內聚力模型和指數型內聚力模型對顆粒填充模型進行不同速度載荷下的數值模擬，通過數值仿真計算得到的脫濕點與試驗結果對比，驗證模型準確性。

1 單軸拉伸與松弛試驗

HTPB主要由粘合劑、AP、Al及RDX組成，為降低實驗成本并保證實驗的安全，采用不添加RDX的HTPB-AP-Al定制配方，配方比例如表1所示。

表1 HTPB推進劑基本組分

其中，AP填充顆粒的粒徑比為250∶150∶20=1∶1∶1。對定制配方HTPB推進劑試件采用5組加載速度的單軸拉伸試驗，分別為1、5、20、100、500 mm/min。實驗在電子式萬能試驗機(深圳REGER公司，型號為RGM-X030)上完成。圖1為HTPB在5組拉伸速率下的工程應力應變曲線。

Al顆粒作為加速HTPB熱分解的添加劑，通常不會對界面的脫濕造成影響。根據文獻[10]對復合基體進行應力松弛試驗，用Prony級數對松弛曲線進行擬合，其擬合的表達式為

(1)

五階Prony級數松弛模量主曲線擬合結果為

(2)

從而得到基體膠片的松弛模量應力應變關系：

(3)

圖1 HTPB推進劑單軸拉伸應力-應變曲線

2 細觀模型建立

2.1 顆粒填充模型

使用掃描電鏡對HTPB斷面形貌進行觀察，掃描倍數為100，如圖2所示。

通常HTPB顆粒填充系數較高，顆粒均勻地分布在基體內部，小顆粒散布于大顆粒周圍，所有顆粒的形狀均類似于球形的多邊形，形狀及其不規則的顆粒極為少見。從圖2可觀察到已脫濕顆粒和顆粒脫濕后留下的凹坑，說明在外載荷的作用下，顆粒/基體脫濕率先導致了推進劑的結構破壞。因此，從細觀角度出發對推進劑的失效破壞進行分析成為有效途徑。

本文采用分子動力學算法生成了隨機分布的圓形顆粒填充模型。多邊形顆粒是在圓形顆粒基礎上使用VB語言生成的。根據圓形顆粒半徑及圓心坐標位置確定多邊形顆粒的位置及大小：

(4)

式中xp與yp為多邊形顆粒的x軸與y軸坐標；xc與yc為圓形顆粒的x軸與y軸坐標；r為圓形顆粒的半徑；n為生成隨機多邊形的邊數，本文選用十四邊形，即n=14；i為多邊形的第i條邊。

按照表1給出的HTPB各顆粒組分含量，生成了如圖3所示的隨機顆粒填充模型。

圖2 HTPB的斷口形貌

(a)圓形顆粒

(b)多邊形顆粒

由于Al和AP粒徑相差較大，會降低細觀模型的網格質量，并影響數值仿真結果。為提高網格質量，本文將Al與基體粘結劑視為一體，其初始模量為1.232 8 MPa。根據文獻[5]獲得AP顆粒的相關參數，彈性模量為32 447 MPa，泊松比為0.143 3 MPa。

2.2 Cohesive element二次開發

Cohesive element有幾何厚度和計算厚度。先前使用VB語言在基體/顆粒之間嵌入了有限厚度Cohesive element進行數值仿真，然而由于粘結層的厚度很小，大約為0.000 1 mm，這種方法通常會造成網格質量下降，進而影響仿真結果收斂性。為了解決這個問題，現使用Python腳本語言開發零厚度Cohesive element。下面詳述零厚度Cohesive element的生成：

(1)將建立好的隨機顆粒填充模型導入ABAQUS中，分別建立顆粒單元集合與基體單元集合，采用四節點的應變單元對所有單元進行網格劃分，并以ABAQUS.inp文件格式輸出。

(2)找出顆粒單元與基體單元的公共點，命名為界面集合。

(3)對所有公共節點進行拆分，拆分后節點編號進行重新排列。需注意的是，粘結單元的編號必須使用逆時針編號，否則ABAQUS不能識別。其網格類型為COH2D4。

(4)將新生成的文件以ABAQUS.inp的文件格式重新導入ABAQUS中，進行數值模擬。

2.2.1 雙線性內聚力模型

雙線性內聚力模型廣泛應用于材料與結構的開裂破壞研究中，最早在研究脆性材料的斷裂時被提出[12]，后來將內聚力模型與有限元方法相結合，實現數值模擬的材料斷裂[13]。典型的雙線型內聚力模型見圖4。

雙線性內聚力模型張力位移關系的控制方程為

(5)

(6)

(a)法向

(b)切向

(7)

在雙線性張力位移關系中，除了最大應力值和臨界能值作為參數必須給出外，還需給界面臨界張開位移或應力上升階段斜率K。

2.2.2 指數型內聚力模型

ABAQUS提供了三種內聚力子程序開發方法：自定義界面本構子程序、自定義單元子程序及自定義材料子程序。自定義材料子程序方法包括UMAT和VUMAT兩種方法。其中，UMAT不能和零厚度的粘結單元共同使用。因此，本文選用VUMAT來實現內聚力模型的張力位移關系。VUMAT通過更新狀態變量矩陣與材料常數矩陣PROPS(*)實現材料的本構關系。狀態變量矩陣為StateOld(*)和StateNew(*)。應力張量矩陣為StressOld(*)和StressNew(*)。StateOld(*)和StressOld(*)代表上一增量步的狀態變量值與應力張量；StateNew(*)和StressNew(*)代表增量步結束時的狀態變量值與應力張量。例如，在第n個增量步中有：

StateOld(n,i)=StressNew(n-1,i)

StressOld(n,i)=StressNew(n-1,i)

(8)

式中n代表第n個增量步；i代表第i個狀態變量或應力張量。

如圖5所示，指數內聚力模型是連續性方程，其應力變化是連續的。應力在減小的過程中呈非線性地漸進為零。指數型內聚力模型選用斷裂能作為判斷單元損傷的條件：

(9)

其中，Δn、Δt為界面在法向和切向的位移值，應力最大點對應的位移值為特征位移值，分別為δn和δt，φn為純法向開裂時所需的斷裂能; 參數q、r分別為

(10)

(11)

(a)法向

(b)切向

2.3 內聚力模型的參數獲取

內聚力模型參數的準確性直接影響數值模擬結果的正確性。在目前的工程應用和理論分析中，內聚力參數是很難直接由實驗獲得的，大多數取值來源于文獻和經驗，但是這種方法得到的模型參數往往不精確。本文通過Hooke-Jeeves的反演方法[14]得到模型參數，其基本原理是通過實驗結果擬合出模型參數，如圖6所示。在擬合過程中通過不斷改變參數取值來調整仿真結果，使得仿真結果與實驗結果的誤差小于規定值。

圖6 界面反演流程圖

雙線性內聚力模型需確定3個參量，分別為界面剛度K、內聚強度σmax和破壞位移δf。首先在模型中設定3個參數的初始值：界面剛度K=200 MPa/mm，內聚強度σmax=0.2 MPa，破壞位移δf=0.01 mm。依據仿真應力應變曲線minR得到應力值，并與試驗曲線獲得應力值進行比較。通過目標函數minR得知，試驗值與仿真值的誤差大小為

(12)

通過上式對函數值進行迭代，不斷更新界面剛度、內聚強度和破壞位移，使得minR小于規定的容差Rlim。當反演結果與實際實驗結果的minR符合規定值時，認為反演參數具備可行性。指數內聚力模型的參數反演和雙線性內聚力模型參數獲得途徑一致。指數內聚力模型只需確定內聚強度σmax和破壞位移δf，設定初始內聚力強度σmax=0.2 MPa和破壞位移δf=0.01 mm。反演后得到的參數如表2所示。

表2 雙線性內聚力模型及指數型內聚力模型界面參數

3 結果與討論

3.1 仿真結果

使用Abaqus軟件對圓形顆粒與多邊形顆粒模型進行5組加載速度下的仿真計算，由于圖形較多，只給出采用指數型內聚力模型的圓形顆粒填充模型與多邊形顆粒填充模型在拉伸速度為1 mm/min時的Mises應力云圖，如圖7和圖8所示。

(a)應變3% (b)應變10% (c)應變15%

(a)應變3% (b)應變10% (c)應變15%

圖7(a)和圖8(a)為應變為3%時的Mises應力云圖變化。此時，顆粒與基體均處于彈性變化范圍內，其形狀沒有發生大的改變，只是云圖的顏色發生了改變，說明顆粒與基體的應力分布不均勻，顆粒與基體之間正在進行力的傳遞。

圖7(b)和圖8(b)為應變為10%時Mises應力云圖變化。在大粒徑圓形顆粒與多邊形顆粒的極區位置，均發生了界面位移分離。這是因為顆粒與基體的材料屬性不同，其模量也不同。顆粒模量大于基體模量。在作用力相同的情況下，顆粒變形小于基體變形，使得界面成為最薄弱的部位，故發生位移分離現象。

當應變達到15%時，如圖7(c)和圖8(c)所示，在顆粒聚集較密集的地方也發生了界面位移分離。已經脫濕的顆粒從孔洞中裸露出來，基體明顯伸長。此時，顆粒與基體的界面不再承擔外載荷作用，即到達了顆粒的最大脫濕點。

3.2 結果分析

通過圖7和圖8發現，在外載荷的作用下，顆粒與基體之間的界面會產生分離現象，進而造成顆粒/基體界面脫濕和孔洞的生成。為進一步探究顆粒與基體界面分離對推進劑力學性能的影響并驗證模型準確性，將對實驗獲得的初始模量及脫濕點模量與仿真模型得出的初始模量及脫濕點模量進行比較，以驗證各模型的預測能力和確定適用范圍。

如圖9所示，隨著加載速率的提高，HTPB固體推進劑的初始模量增大，表現出顯著的粘彈性效應。使用雙線性內聚力圓形顆粒填充模型和多邊形顆粒填充模型的初始模量預測結果明顯大于使用指數型內聚力圓形顆粒填充模型與多邊形顆粒填充模型。這是因為指數型內聚力模型與雙線性內聚力模型在計算法則不同，雙線性內聚力模型更適用于描述脆性材料的本構，而指數型內聚力模型適合描述非線性材料的本構。因此，指數型內聚力模型的預測值與試驗結果誤差更小。

圖9 不同拉伸速率下的初始模量

當應力應變曲線的斜率明顯降低時，即發生了顆粒的脫濕，這一點稱為顆粒/基體界面脫濕點。如圖10所示，和初始模量的預測結果類似，指數型內聚力模型比雙線性內聚力模型的預測結果更為準確。因為兩種模型的損傷變量形式不同，指數型內聚力模型的損傷模量要高于雙線性內聚力模型的損傷模量。所以，指數型內聚力模型的損傷速度較快，使得界面能更快軟化，得到脫濕點模量值更加準確。多邊形顆粒填充模型較圓形顆粒填充模型得到脫濕點模量更準確，這是因為AP顆粒的形狀并不是真的圓形，為計算方便，將AP顆粒簡化成圓形進行數值模擬。因此，多邊形顆粒填充模型比圓形顆粒填充模型更接近推進劑內部的真實情況。

圖10 不同拉伸速率下的脫濕點模量

從圖11可看出，拉伸速率越大，模量下降率越大，并且不同加載速率的脫濕點模量近等相似，說明HTPB固體推進劑這種材料在高速率外載荷的作用下，更容易出現損傷。為保證推進劑的力學性能，下一步應探究在高速率的工作載荷下，延緩推進劑損傷的措施。

圖11 模量下降率

數值模擬得出的模量值均略高于試驗值，這是因為本文所建立的模型均是基于二維平面應變，而推進劑真實形態是三維狀態的，且本文在建立模型時，認為模型的內部是不存在微小缺陷的。在實際情況中，推進劑的內部存在微孔洞等缺陷，這會造成應力值的下降。以上原因造成了模型的計算值大于試驗值。

4 結論

(1)分子動力學方法生成的多邊形顆粒填充模型比圓形顆粒填充模型能更好地描述推進劑細觀結構，更準確地判斷損傷位置。

(2)由于目前細觀試驗手段的局限性，引入Hook-Jeeves優化反演方法，獲得了更準確的內聚力模型參數。通過不同損傷形式內聚力模型的比較，指數內聚力模型更適用于延展性材料的損傷，能更好地預測推進劑的力學性能，使得仿真結果更接近于試驗值。

(3)通過5組加載速率的試驗結果與仿真結果比較，發現推進劑的損傷和加載速率相關，載荷加載速率越大，推進劑受損程度越大，模量下降率越高。

(4)下一步將開發三維推進劑顆粒填充模型，以提高模型的準確性。