侯宇菲，許進(jìn)升，陳 雄，周長(zhǎng)省，朱 亮

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

0 引言

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

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

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

1 單軸拉伸與松弛試驗(yàn)

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

表1 HTPB推進(jìn)劑基本組分

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

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

(1)

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

(2)

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

(3)

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

2 細(xì)觀模型建立

2.1 顆粒填充模型

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

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

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

(4)

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

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

圖2 HTPB的斷口形貌

(a)圓形顆粒

(b)多邊形顆粒

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

2.2 Cohesive element二次開發(fā)

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

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

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

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

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

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

雙線性內(nèi)聚力模型廣泛應(yīng)用于材料與結(jié)構(gòu)的開裂破壞研究中，最早在研究脆性材料的斷裂時(shí)被提出[12]，后來將內(nèi)聚力模型與有限元方法相結(jié)合，實(shí)現(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來實(shí)現(xiàn)內(nèi)聚力模型的張力位移關(guān)系。VUMAT通過更新狀態(tài)變量矩陣與材料常數(shù)矩陣PROPS(*)實(shí)現(xiàn)材料的本構(gòu)關(guān)系。狀態(tài)變量矩陣為StateOld(*)和StateNew(*)。應(yīng)力張量矩陣為StressOld(*)和StressNew(*)。StateOld(*)和StressOld(*)代表上一增量步的狀態(tài)變量值與應(yīng)力張量；StateNew(*)和StressNew(*)代表增量步結(jié)束時(shí)的狀態(tài)變量值與應(yīng)力張量。例如，在第n個(gè)增量步中有：

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

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

(8)

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

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

(9)

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

(10)

(11)

(a)法向

(b)切向

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

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

圖6 界面反演流程圖

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

(12)

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

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

3 結(jié)果與討論

3.1 仿真結(jié)果

使用Abaqus軟件對(duì)圓形顆粒與多邊形顆粒模型進(jìn)行5組加載速度下的仿真計(jì)算，由于圖形較多，只給出采用指數(shù)型內(nèi)聚力模型的圓形顆粒填充模型與多邊形顆粒填充模型在拉伸速度為1 mm/min時(shí)的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)?%時(shí)的Mises應(yīng)力云圖變化。此時(shí)，顆粒與基體均處于彈性變化范圍內(nèi)，其形狀沒有發(fā)生大的改變，只是云圖的顏色發(fā)生了改變，說明顆粒與基體的應(yīng)力分布不均勻，顆粒與基體之間正在進(jìn)行力的傳遞。

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

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

3.2 結(jié)果分析

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

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

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

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

圖10 不同拉伸速率下的脫濕點(diǎn)模量

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

圖11 模量下降率

數(shù)值模擬得出的模量值均略高于試驗(yàn)值，這是因?yàn)楸疚乃⒌哪Ｐ途腔诙S平面應(yīng)變，而推進(jìn)劑真實(shí)形態(tài)是三維狀態(tài)的，且本文在建立模型時(shí)，認(rèn)為模型的內(nèi)部是不存在微小缺陷的。在實(shí)際情況中，推進(jìn)劑的內(nèi)部存在微孔洞等缺陷，這會(huì)造成應(yīng)力值的下降。以上原因造成了模型的計(jì)算值大于試驗(yàn)值。

4 結(jié)論

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

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

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

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