復合固體推進劑顆粒與基體初始界面有無缺陷的細觀模型對比

侯宇菲，許進升，周長省，陳雄，李宏文

(1.南京理工大學 機械工程學院，江蘇 南京 210094；2.晉西工業(yè)集團有限責任公司，山西 太原 030027)

0 引言

端羥基聚丁二烯(HTPB)推進劑作為一種多相(基體相、增強相、界面相)含能材料，廣泛應用于各類火箭武器。HTPB的宏觀力學性能取決于顆粒、基體及顆粒與基體界面的力學性質(zhì)，因此從細觀角度上視為一種非均質(zhì)材料[1]。van Ramshorst等[2]和王陽等[3]基于掃描電鏡(SEM)對HTPB推進劑細觀形態(tài)與宏觀變形及失效之間的關(guān)系進行了研究，發(fā)現(xiàn)在外載荷作用下，顆粒與基體界面發(fā)生脫粘，且裂紋沿著顆粒與基體界面生長，從而引起推進劑力學性能非線性。因此，有必要對復合固體推進劑顆粒與基體界面脫粘進行分析。

Matou?等[4]、常武軍等[5]、Inglis[6]在小應變條件下，采用隨機夾雜顆粒填充模型對推進劑顆粒與基體界面脫粘過程進行了數(shù)值模擬。劉承武等[7]考慮推進劑基體材料大變形的特點，建立了Mori-Tanaka有限元方法，然而隨著界面脫粘程度的加深，基于Mori-Tanaka有限元法計算的應力值偏高，從而無法準確描述推進劑加速損傷的特點。韓龍[8]和Dai等[9]通過修正內(nèi)聚力法則，得到不同應變率下顆粒與基體界面非線性脫粘的損傷演化規(guī)律，然而模型的預測應力值均大于實驗應力值。究其原因，是因為以上研究假設固體推進劑顆粒與基體界面為有效粘結(jié)，并未考慮初始界面缺陷對推進劑力學性能的影響，因此導致預測應力值偏高。Balzer等[10]通過SEM對推進劑顆粒與基體界面在拉伸初始時刻進行觀察，發(fā)現(xiàn)部分顆粒與基體之間存在微小孔洞。這是因為推進劑在固化降溫過程中顆粒與基體的材料屬性不同，導致其形變存在差異，造成顆粒與基體之間形成間隙，使得推進劑內(nèi)部存在隨機分布的初始界面缺陷。劉晉湘等[11]通過顆粒與基體界面含初始缺陷試樣以及顆粒與基體界面無初始缺陷試樣的單軸拉伸實驗發(fā)現(xiàn)，初始界面缺陷決定材料的失效路徑，并且高氯酸銨(AP)顆粒與基體初始界面損傷會降低初始界面的粘結(jié)強度，導致推進劑更容易發(fā)生損傷。黃海峰等[12]和文志杰等[13]在研究巖石斷裂時發(fā)現(xiàn)初始缺陷對細觀非均勻材料宏觀力學性能的影響不容忽略，并利用Weibull函數(shù)對巖石材料的本構(gòu)模型進行修正，能夠更準確地描述這類材料的非線性力學特征。

本文采用通用有限元分析軟件Abaqus的Abaqus/Explicit求解器，結(jié)合VUMAT獲得指數(shù)內(nèi)聚力法則，利用Weibull分布作為顆粒與基體界面粘結(jié)強度的密度函數(shù)，建立考慮初始缺陷的細觀顆粒填充模型，通過對比數(shù)值仿真與實驗結(jié)果驗證所建立模型的準確性。

1 細觀顆粒填充模型

1.1 顆粒填充模型

以HTPB復合固體推進劑為例，采用分子動力學算法[8]生成與實際推進劑具有相同體積分數(shù)的細觀顆粒填充模型，填充比例如表1所示。其中，填充鋁(Al)粉與AP顆粒粒徑分布服從圖1所示的雙峰分布規(guī)律。根據(jù)文獻[14]得知，當代表性體積單元尺寸為最大顆粒粒徑的3～5倍時，認為其有效。因此，本文選用的代表性體積單元尺寸為500 μm×500 μm，如圖2所示。

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

圖2 細觀顆粒填充模型Fig.2 Microscopic particle packing model

圖2中，AP顆粒與Al顆粒的最大粒徑之比約為20∶1，若直接對該模型進行網(wǎng)格劃分，則將導致網(wǎng)格質(zhì)量降低，進而仿真計算結(jié)果無法收斂。文獻[15]通過實驗證明，由于Al顆粒與基體材料的界面浸潤能力較強，固體推進劑在攪拌固化過程中，顆粒與基體界面無脫粘現(xiàn)象。因此，通過Mori-Tanaka[16]解析法將Al顆粒等效到基體材料中，在數(shù)值模擬過程中只考慮AP顆粒對材料力學性能的影響，并將此類基體稱為等效基體，材料參數(shù)如表2所示。

表2 HTPB推進劑的材料參數(shù)Tab.2 Material parameters of HTPB propellant

1.2 粘結(jié)單元

為模擬顆粒與基體間的界面，本文采用Python腳本語言開發(fā)零厚度粘結(jié)單元，如圖3所示。為方便理解，圖3給出了粘結(jié)單元的幾何厚度，實際計算過程中幾何厚度為0 mm.

圖3 粘結(jié)單元嵌入示意圖Fig.3 Schematic diagram of cohesive element embedding

1.3 指數(shù)內(nèi)聚力法則

內(nèi)聚力法則最早由Dugdale[17]和Barenblatt[18]提出，他們認為在裂紋尖端處存在一個斷裂過程區(qū)，此區(qū)域呈現(xiàn)了裂紋尖端及附近的損傷機制，并引入了內(nèi)聚力和間距的本構(gòu)關(guān)系，以及內(nèi)聚力法則，通過該法則可以模擬斷裂區(qū)域內(nèi)能量的耗散過程。內(nèi)聚力模型根據(jù)其曲線形式可分為雙線性、多項式、指數(shù)型等，每種形式的應力- 位移關(guān)系都有其適用范圍。根據(jù)推進劑在載荷作用下的加速損傷特性，本文采用指數(shù)型內(nèi)聚力法則描述顆粒與基體界面的損傷演化過程。

由于在計算過程中使用了0 mm厚度粘結(jié)單元，指數(shù)內(nèi)聚力法則需通過VUMAT子程序二次開發(fā)得到[19]。指數(shù)內(nèi)聚力法則的斷裂能φ(t)為

(1)

式中：φn為純法向開裂狀態(tài)下的界面斷裂能；tn和tt為界面的法向與切向位移；δn和δt為法向與切向界面開裂時的特征位移；參數(shù)q和r分別由(2)式和(3)式確定：

(2)

(3)

(4)

(5)

式中：Tn與Tt分別為法向應力值與切向應力值。

通過(4)式和(5)式發(fā)現(xiàn)，指數(shù)內(nèi)聚力法則中的Tn和Tt由法向位移和切向位移同時決定，因此更符合實際界面開裂的應力狀態(tài)。

1.4 Weibull分布

復合固體推進劑在制備過程中受工藝和材料特性的限制，內(nèi)部存在大量隨機分布的初始缺陷。其中包括少量氣孔和攪拌過程中由于顆粒碰撞引起的顆粒破碎，更多的缺陷是在固化降溫時，顆粒與基體界面浸潤能力下降導致的顆粒與基體界面間產(chǎn)生空隙，如圖4所示。這些初始缺陷的存在導致HTPB復合固體推進劑在進行單軸拉伸重復性實驗時，應力- 應變曲線的分散性增加，進一步導致材料整體強度降低，因此在預測推進劑的力學性能時需要考慮顆粒與基體界面間的初始缺陷。通常使用Abaqus軟件對顆粒與基體界面材料屬性進行定義時，會默認所有顆粒與基體界面屬性均相同。目前有兩種方法可以表達材料細觀界面的非均質(zhì)性：1)將顆粒與基體界面劃分為不同集合；2)采用自開發(fā)程序，將材料屬性設置為積分點的函數(shù)。方法1由于計算效率低而無法適用于高顆粒填充比材料。方法2可以隨機定義材料屬性，能夠準確模擬界面開裂時的真實狀態(tài)，因此本文采用方法2實現(xiàn)顆粒與基體界面的非均質(zhì)性。文獻[20]認為航空、航天工業(yè)中的復合材料特性符合Weibull分布。因此本文將顆粒與基體界面粘結(jié)強度的密度函數(shù)設置為Weibull函數(shù)，建立考慮初始缺陷的HTPB復合固體推進劑細觀數(shù)值模型。Weibull函數(shù)概率密度f(x)表達式為

圖4 固體推進劑初始缺陷示意圖Fig.4 Initial defect of solid propellant

(6)

式中：x表示隨機函數(shù)；m為形狀參數(shù)；λ為比例參數(shù)。形狀參數(shù)m反映了f(x)的離散程度，m越大，分散性越小，材料的可靠度越高，如圖5所示。

圖5 Weibull分布密度函數(shù)Fig.5 Weibull distribution density function

1.5 內(nèi)聚力參數(shù)及邊界條件

為模擬HTPB推進劑在單軸等速拉伸下顆粒與基體界面脫粘的過程，對所建立的細觀模型選用合理邊界條件進行數(shù)值模擬，如圖6所示。圖6(a)、圖6(b)、圖6(d)的下邊界施加固定約束載荷、上邊界作為載荷承載端，其中：圖6(a)左邊界和右邊界為無約束邊界；圖6(b)左邊界和右邊界為豎直邊界，即在變形過程中邊界保持豎直；圖6(d)左邊界和右邊界為周期性邊界；圖6(c)上、下、左、右邊界均為周期性邊界。文獻[21]認為若代表性體積單元尺寸選取適當，則圖6(a)、圖6(b)、圖6(c)的邊界條件對數(shù)值計算結(jié)果影響較小。圖6(d)中，由于軸約束不匹配，導致計算結(jié)果與實際情況不符。本文所考慮的顆粒分布狀態(tài)呈非周期性，故選用圖6(b)中的邊界條件進行數(shù)值計算。

圖6 邊界條件Fig.6 Boundary conditions

由于目前實驗手段的局限性，很難獲得微米級顆粒與基體界面間的粘結(jié)能，普遍采用的方法是通過實驗反演法獲得內(nèi)聚力模型的材料參數(shù)。主要包括如下3個步驟：1)采用材料實驗機對HTPB復合固體推進劑試樣進行單軸拉伸實驗，實驗環(huán)境溫度為20 ℃，濕度保持在40%±5%范圍內(nèi)，實驗機加載速率為1 mm/min，實驗重復3次，從而獲得HTPB復合固體推進劑的單軸拉伸應力- 應變曲線，如圖7所示；2)建立考慮初始缺陷的固體推進劑細觀模型，假設試樣3無顆粒與基體初始界面缺陷，并將其應力- 應變曲線定為目標函數(shù)；3)根據(jù)文獻[22-23]中的Hook-Jeeves優(yōu)化算法使仿真應力- 應變曲線盡可能接近目標函數(shù)，當仿真曲線與目標函數(shù)誤差值小于設定閾值時，確定內(nèi)聚力模型材料參數(shù)，具體反演過程見文獻[22]。

圖7 實驗應力- 應變曲線Fig.7 Experimental stress-strain curves of propellant

指數(shù)內(nèi)聚力模型只需確定粘結(jié)強度σmax和特征位移δc即可，反演后的參數(shù)如表3所示。

表3 指數(shù)內(nèi)聚力參數(shù)Tab.3 Exponential cohesive parameters

2 結(jié)果分析

2.1 無缺陷模型仿真結(jié)果

為消除有限元模型網(wǎng)格尺寸所帶來的計算誤差，對模型的網(wǎng)格收斂性進行計算，取網(wǎng)格尺寸y分別為0.040 mm、0.020 mm、0.010 mm、0.005 mm，如圖8所示。4種網(wǎng)格尺寸的計算模型，其顆粒與基體界面采用指數(shù)內(nèi)聚力法則，各點粘結(jié)強度均相同，不考慮初始缺陷，稱這種模型為無缺陷模型。由于計算結(jié)果較多，下面僅給出網(wǎng)格尺寸為0.01 mm的模型應力計算云圖，如圖9所示。

圖8 網(wǎng)格尺寸Fig.8 Mesh size

從圖9中可以看出，顆粒與基體材料的應力變化分為3個階段。當應變ε為0.03時，模型整體未出現(xiàn)損傷，此時推進劑處于彈性階段。當應變ε為0.05時，顆粒作為彈性體，其彈性模量遠大于基體材料彈性模量，當應變相同時，顆粒承受的應力大于基體承受的應力，導致顆粒與基體界面成為推進劑結(jié)構(gòu)中的最薄弱區(qū)域，在載荷的持續(xù)作用下，顆粒與基體界面發(fā)生位移分離。對于單一粒徑配方，模型中部的顆粒與基體界面最先發(fā)生損傷。當應變ε為0.10時，顆粒與基體界面發(fā)生脫粘，部分顆粒從微孔洞中裸露出來，此時顆粒兩極不再受力，其赤道處出現(xiàn)高應力區(qū)。最終，高應力區(qū)形成微裂紋，進而導致推進劑失效破壞。

下面通過彈性模量- 應變曲線分析推進劑細觀失效時的力學響應，如圖10所示。由圖10可見：拉伸初始階段，初始彈性模量恒定，模型處于彈性階段；當應變繼續(xù)增大時，彈性模量開始下降，結(jié)合應力云圖可以發(fā)現(xiàn)顆粒與基體界面發(fā)生脫粘，由此得知顆粒與基體界面脫粘會導致推進劑彈性模量下降。

圖10 無損傷模型彈性模量- 應變曲線Fig.10 Modulus-strain curves of damage-free model

對比圖10中不同網(wǎng)格尺寸的無缺陷模型預測結(jié)果與試樣3的實驗結(jié)果發(fā)現(xiàn)，無缺陷模型對網(wǎng)格尺寸的敏感性較高，模型預測結(jié)果隨著網(wǎng)格尺寸的減小而接近實驗結(jié)果。當網(wǎng)格尺寸為0.020 mm時，模型對網(wǎng)格尺寸的敏感性大幅度降低，計算結(jié)果收斂。為考察網(wǎng)格尺寸對細觀模型應力場的影響，下面給出模型整體應變?yōu)?.06時，模型網(wǎng)格尺寸為0.005 mm與0.040 mm的應力云圖，如圖11所示。由圖11可見：粗網(wǎng)格模型的應力過渡區(qū)有大量鋸齒狀區(qū)域出現(xiàn)，且連續(xù)性較差，過度不均勻；細網(wǎng)格模型的應力過度區(qū)光滑，其應力場分布連續(xù)性較好。以上分析表明，為提高無缺陷模型計算結(jié)果的準確性，需采用合適的網(wǎng)格尺寸。

圖11 不同網(wǎng)格尺寸的von Mises應力云圖Fig.11 von Mises stress nephograms with different mesh sizes

2.2 初始缺陷損傷模型結(jié)果分析

采用2.1節(jié)的4種網(wǎng)格尺寸，對初始缺陷模型進行數(shù)值仿真，其顆粒與基體界面粘結(jié)強度的密度函數(shù)呈Weibull分布。從圖5得知，粘結(jié)強度分布密度與形狀參數(shù)m有關(guān)，為確定符合復合固體推進劑材料粘結(jié)強度的分布規(guī)律，將m分別取值為2.0、2.5、3.0、3.5、4.0、4.5、5.0代入計算模型中進行數(shù)值模擬。由于模型數(shù)量較多，下面只給出m=3.0時網(wǎng)格尺寸為0.010 mm的初始缺陷模型應力云圖，如圖12所示。

圖12 初始缺陷損傷模型的von Mises應力云圖Fig.12 von Mises stress nephograms of initial damage model

由圖12可見，初始缺陷模型與無缺陷模型的應力云圖變化過程類似，其應力云圖變化同樣分為3個階段，但二者的第2階段存在明顯區(qū)別。在無缺陷模型中，位于模型中部的顆粒最先發(fā)生脫粘，并且由于顆粒與基體界面的粘結(jié)強度相同，顆粒兩級處應力最大，顆粒與基體界面脫粘從極點沿兩側(cè)緩慢進行。而在初始缺陷模型中，由于界面粘結(jié)強度分布的非均勻性，顆粒與基體界面脫粘位置發(fā)生改變，不再是位于模型中部的顆粒率先發(fā)生界面脫粘，而是顆粒與基體界面損傷含量較高的位置率先發(fā)生脫粘；并且位于極點周圍的區(qū)域均為初始脫粘區(qū)，如圖12(b)所示。初始缺陷模型更好地解釋了推進劑在進行重復拉伸實驗時試件斷裂位置的不可預估性。

為了更直觀地分析初始缺陷對復合固體推進劑力學性能的影響，本文給出了初始缺陷模型的彈性模量- 應變曲線，如圖13所示。

圖13 初始缺陷損傷模型的彈性模量- 應變曲線Fig.13 Modulus-strain curves of initial damage model

從圖13(a)～圖13(g)中可以看出，初始彈性模量隨形狀參數(shù)m值降低而減小。這是因為m值決定了顆粒與基體界面粘結(jié)強度的分布規(guī)律，m值越小，粘結(jié)強度分布越分散，細觀模型內(nèi)部的初始缺陷越多，顆粒與基體界面發(fā)生脫粘損傷的概率越大，從而界面?zhèn)鬟f載荷能力下降，導致模型的整體強度下降。在相同應變條件下對應的應力越低，最終導致初始彈性模量下降。相反，m值越大，粘結(jié)強度分布越集中，數(shù)值計算結(jié)果得到的初始彈性模量越高。為了進一步證明初始缺陷損傷模型是否符合真實推進劑的力學狀態(tài)，將所有模型的計算結(jié)果求均值并與實驗結(jié)果比較，如圖13(h)所示。從圖13(h)中可以看出，仿真結(jié)果與實驗結(jié)果基本吻合，表明本文所建模型有效，同時也表明推進劑在制備過程中，推進劑內(nèi)部存在隨機分布的初始缺陷，因此需要采用考慮初始缺陷的細觀模型對其進行力學性能分析。

對比圖10與圖13(h)發(fā)現(xiàn)，初始缺陷損傷模型對網(wǎng)格尺寸的敏感性大幅度降低，計算所得初始彈性模量最大值與最小值的相對誤差小于4%，提高了網(wǎng)格收斂性。眾所周知，網(wǎng)格數(shù)量越多，模型的計算時間越長，導致計算成本增加。因此，可以采用初始缺陷損傷模型代替無缺陷損傷模型進行數(shù)值計算，從而降低計算成本。

3 結(jié)論

本文基于Abaqus有限元軟件建立了HTPB復合固體推進劑顆粒填充模型，分析了顆粒與基體初始界面缺陷對推進劑力學性能的影響。得出如下主要結(jié)論：

1) 針對HTPB復合固體推進劑顆粒與基體初始界面缺陷問題，本文以Weibull分布作為顆粒與基體界面粘結(jié)強度的密度函數(shù)，分別建立了無缺陷模型與初始缺陷模型。數(shù)值仿真結(jié)果表明：初始彈性模量隨著顆粒與基體初始界面缺陷含量的增多而降低。

2) HTPB復合固體推進劑細觀損傷分為3個階段：彈性階段、顆粒與基體之間產(chǎn)生位移分離、顆粒與基體界面脫粘并形成孔洞。結(jié)合仿真結(jié)果與應力應變曲線得知：顆粒與基體界面脫粘引起推進劑彈性模量降低，改善顆粒與基體之間的浸潤能力是提高復合固體推進劑力學性能的關(guān)鍵。

3) 網(wǎng)格尺寸敏感性分析表明：無缺陷模型與初始缺陷模型的數(shù)值計算結(jié)果隨著網(wǎng)格尺寸減小而收斂，并且初始缺陷模型的網(wǎng)格尺寸敏感性較低，在后續(xù)工作中可采用初始缺陷模型來降低計算成本。