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

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

0 引言

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

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

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

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組拉伸速率下的工程應(yīng)力應(yīng)變曲線。

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

(1)

五階Prony級數(shù)松弛模量主曲線擬合結(jié)果為

(2)

從而得到基體膠片的松弛模量應(yīng)力應(yīng)變關(guān)系：

(3)

圖1 HTPB推進劑單軸拉伸應(yīng)力-應(yīng)變曲線

2 細觀模型建立

2.1 顆粒填充模型

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

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

本文采用分子動力學算法生成了隨機分布的圓形顆粒填充模型。多邊形顆粒是在圓形顆?；A(chǔ)上使用VB語言生成的。根據(jù)圓形顆粒半徑及圓心坐標位置確定多邊形顆粒的位置及大?。?/p>

(4)

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

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

圖2 HTPB的斷口形貌

(a)圓形顆粒

(b)多邊形顆粒

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

2.2 Cohesive element二次開發(fā)

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

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

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

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

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

2.2.1 雙線性內(nèi)聚力模型

雙線性內(nèi)聚力模型廣泛應(yīng)用于材料與結(jié)構(gòu)的開裂破壞研究中，最早在研究脆性材料的斷裂時被提出[12]，后來將內(nèi)聚力模型與有限元方法相結(jié)合，實現(xiàn)數(shù)值模擬的材料斷裂[13]。典型的雙線型內(nèi)聚力模型見圖4。

雙線性內(nèi)聚力模型張力位移關(guān)系的控制方程為

(5)

(6)

(a)法向

(b)切向

(7)

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

2.2.2 指數(shù)型內(nèi)聚力模型

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

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

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

(8)

式中n代表第n個增量步；i代表第i個狀態(tài)變量或應(yīng)力張量。

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

(9)

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

(10)

(11)

(a)法向

(b)切向

2.3 內(nèi)聚力模型的參數(shù)獲取

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

圖6 界面反演流程圖

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

(12)

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

表2 雙線性內(nèi)聚力模型及指數(shù)型內(nèi)聚力模型界面參數(shù)

3 結(jié)果與討論

3.1 仿真結(jié)果

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

(a)應(yīng)變3% (b)應(yīng)變10% (c)應(yīng)變15%

(a)應(yīng)變3% (b)應(yīng)變10% (c)應(yīng)變15%

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

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

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

3.2 結(jié)果分析

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

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

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

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

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

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

圖11 模量下降率

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

4 結(jié)論

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

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

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

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