固體推進劑低溫細觀損傷仿真研究

周紅梅，袁 軍，賴建偉，袁 嵩

(1.海軍航空工程學院七系，煙臺 264001；2.海軍西安地區導彈發動機專業軍代室，西安 710025)

0 引言

復合固體推進劑是一種顆粒填充復合材料，主要由金屬、氧化劑顆粒和高聚物粘合劑基體組成。粘合劑基體為粘彈性連續相，具有強烈的時間和載荷歷史效應；固體顆粒作為填料分布于基體中，對基體起增強作用，并顯著提高固體推進劑承載能力。固體推進劑的這種細觀結構特點滿足了各方面性能要求，但同時也導致力學損傷的復雜性。從細觀力學的角度來看，作為一種非均勻材料，固體顆粒的形狀、尺寸和分布存在明顯差異，隨著外界載荷的作用，易在微細觀程度上產生基體裂紋和顆粒/基體界面脫濕等損傷。特別是在低溫條件下，基體、顆粒及界面的力學性能發生變化，從而引起損傷發生的時機、方式和位置發生改變。因此，有必要針對低溫條件下推進劑的各種損傷模式進行分析，以期找出環境溫度影響機理及推進劑宏觀力學性能發生變化的原因。

顆粒/基體界面脫濕和基體斷裂是復合固體推進劑材料的兩種主要損傷模式[1]。目前研究大多針對界面脫濕一種模式[2-3]，對于界面脫濕和基體斷裂同時存在的分析較少涉及。本文采用有限元方法對固體推進劑進行細觀建模，利用Cohesive單元分析脫濕損傷，編制基體VUMAT本構子程序分析基體斷裂，考慮與溫度有關的參數設置，可較為真實地描述固體推進劑的低溫細觀損傷過程。

1 物理模型及計算方法

1.1 細觀建模與自適應網格劃分

為考察不同顆粒數量和直徑大小對推進劑破壞過程的影響，建立單顆粒、4顆粒和8顆粒三種推進劑二維細觀模型。假設高氯酸銨(AP)顆粒為圓形顆粒，不考慮顆粒形狀、取向的影響。其中，單顆粒模型含一個AP大顆粒，4顆粒模型含2個AP大顆粒和2個AP小顆粒，8顆粒模型含2個AP大顆粒和6個AP小顆粒。模型中AP大顆粒直徑為330 μm，AP小顆粒直徑為130 μm。4顆粒模型中顆粒呈對稱分布，8顆粒模型中顆粒呈隨機分布。

采用以Delaunay法為基礎的“力平衡”網格優化算法[4-6]，利用距離函數來描述劃分的區域，采用“力平衡”法更新節點，從而實現重置拓撲結構和網格優化，具有較強的自適應性，可滿足復合固體推進劑細觀結構的網格劃分，三種模型細觀網格如圖1所示。由網格劃分局部放大圖(圖1(d))可看出，推進劑顆粒邊界處的網格根據曲率自動調整疏密，顆粒內部的網格尺寸較邊界處的要大，大顆粒邊界處的網格粗，小顆粒邊界處的網格細，顆粒間距小的區域算法能自動進行加密處理。由此可見，劃分的網格充分體現了網格的自適應特點，較好的考慮到細觀結構的幾何特征及在力學載荷作用下的力學特征，為后續仿真分析的模型提供精度保證。

1.2 Cohesive單元模型

ABAQUS有限元分析軟件的Cohesive單元已經廣泛用于模擬推進劑材料的顆粒/基體界面的力學響應，本文采用二維4節點Cohesive單元來表征顆粒/基體界面脫粘情況，如圖2所示。

固體推進劑顆粒/基體界面力學行為可用純法向和切向的雙線性粘結模型表示，如圖3所示。界面主要經歷彈性段、損傷段和完全失效段三個階段[7]：

彈性段：界面受力較小，界面應力與張開位移成線性增加關系。

損傷段：界面應力達到σmax，界面張開位移達到臨界值δ0，界面開始出現損傷，界面應力與張開位移成線性減小關系。

完全失效段：界面張開位移大于δf。界面完全脫粘不再承力。

1.3 基體開裂模擬

文獻[8-9]采用編制VUMAT子程序的方法開展了裂紋擴展過程研究，取得了較為滿意的效果，驗證了方法的可行性。為模擬界面脫粘過程中推進劑基體的裂紋擴展情況，使用ABAQUS中的子程序接口，開發包含失效準則且具有單元刪除功能的基體本構VUMAT子程序，通過采用最大應力失效準則直接處理材料積分點，當距離裂紋尖端附近的某個積分點的應力值達到斷裂臨界應力值時，此單元就開裂，程序就會刪除此單元。該處理方法在網格劃分足夠精細的情況下，可在未知裂紋擴展路徑時，模擬任意未知方向裂紋開裂和擴展，更加接近推進劑破壞的真實情況。

2 推進劑拉伸破壞過程仿真

2.1 模型及參數設置

推進劑細觀模型底部采用固支約束，頂部施加位移載荷。在界面處引入接觸算法，避免界面完全脫粘而導致顆粒與基體之間相互滲透。

在進行有限元分析時需使用與溫度相關的粘合劑基體、顆粒和顆粒/基體界面的相關參數。假設顆粒為線彈性材料，粘合劑基體為粘彈性材料，基體本構模型采用開發的VUMAT子程序；顆粒/基體界面采用Cohesive單元模擬具有3個獨立參數[10]，即界面強度σmax、彈性模量E和失效張開位移δf。低溫條件下參數在參照常溫參數的基礎上，根據相關參數隨溫度變化的特征規律進行確定，計算中需要使用的參數如表1和表2所示。

表1 固體推進劑顆粒/基體界面參數

表2 固體推進劑AP顆粒參數

HTPB推進劑基體為粘彈性材料，其松弛模量表達式為[11]

Em(t)= 1.04e-t/2×10-6+1.50e-t/2×10-5+1.07e-t/2×10-4+

1.10e-t/2×10-3+1.40e-t/2×10-2+1.13e-t/2×10-1+

1.07e-t/2+0.62e-t/2×101+0.58e-t/2×102+

0.23e-t/2×103+0.27e-t/2×104+0.09e-t/2×105+

0.16e-t/2×106+1.19e-t/2×107

基體的時溫等效方程為

在ABAQUS進行顯式分析時需要設置模型中各材料的密度值，如表3所示。

表3 基體/AP顆粒材料密度

2.2 固體推進劑破壞過程仿真結果分析

圖4為包含一個AP顆粒和粘合劑的固體推進劑細觀模型，顆粒體積分數為35%，在模型上方垂直方向上施加拉伸位移載荷。由圖4可知，在顆粒的上極區拉應力最大的地方首先產生脫粘，隨著位移載荷的增加，脫粘面逐漸增大，下極區也逐漸開始脫粘，最終顆粒上下極區完全脫粘。當脫粘發展到一定程度時，顆粒赤道附近基體出現應力集中，在大致沿顆粒赤道方向上，裂紋由脫粘面向粘合劑基體擴展，并最終沿水平方向貫穿基體內部。在整個拉伸過程中顆粒始終保持完好，沒有產生裂紋。通過計算發現，增加顆粒直徑或提高顆粒體積分數，其破壞現象相同；將界面強度增加至高于基體的強度時，界面脫粘仍然先發生，然后發生基體斷裂，與文獻報道的PBX炸藥破壞過程一致[12]；將基體的強度增大時，裂紋由界面向基體擴展的時間延長，顆粒上下脫粘面形成的孔穴變大。從圖4中模擬結果還可看出，拉伸作用下發生了界面脫粘和基體開裂，而顆粒沒有破壞，這是因為AP顆粒的模量比基體模量高出幾個數量級，AP顆粒的內部應力不足以引起顆粒發生破壞。經過仿真計算發現，低溫和常溫條件下單顆粒模型破壞過程仿真結果類似。

為考察顆粒相互作用和不同半徑的影響，建立了由4個AP顆粒和粘合劑基體組成的模型，模型中顆粒半徑相差較大且成對稱分布，破壞過程計算結果如圖5(a)～(c)所示。由圖5(a)～(c)可知，模型中上方兩AP顆粒首先發生脫粘，其中大顆粒仍然首先脫粘，之后下方顆粒開始脫粘，脫粘完全后下方大、小顆粒赤道附近基體先后開始發生斷裂直至破壞。結果表明，有不同大小的顆粒存在時，大顆粒首先發生脫粘，與顯微鏡觀察結果也是一致的[1，7]。由于顆粒之間的影響，導致多個顆粒之間的基體應力較大且分布復雜，該處基體首先發生斷裂，方向具有一定的隨機性，之后顆粒赤道附近基體才開始斷裂。這說明由于顆粒大小不同，應力分布發生變化，顆粒之間的相互作用對顆粒的損傷破壞產生了重要的影響。

對比低溫損傷破壞過程圖5(f)～(h)發現，由于受顆粒之間相互影響，相同加載條件下，低溫使損傷表現更加復雜。低溫條件下，大顆粒仍然首先發生脫粘，其次小顆粒發生脫粘。但基體斷裂與常溫表現不同，在脫粘進行到一定程度時，右上方大顆粒赤道附近基體發生斷裂，左下方大顆粒由于受其他顆粒影響，在上極區處附近基體出現斷裂。隨著載荷的繼續，下方小顆粒赤道附近也開始產生裂紋。說明低溫條件下，基體斷裂更加容易發生在大顆粒附近。

圖6(a)～(c)為8個隨機分布的顆粒拉伸作用下損傷破壞的計算結果，顆粒體積分數為25%。與4顆粒模型結果相同，脫粘首先發生在大顆粒上，然后再發展到其他小顆粒上。顆粒/基體界面脫粘、基體斷裂及裂紋的擴展都不同程度地受到顆粒之間相互作用的影響，在顆粒分布較密的中間小顆粒由于受附近幾個大、小顆粒的影響，上極區的孔穴充分發展，并首先引起基體的裂紋擴展，之后附近顆粒基體也逐漸發生斷裂。與4顆粒模型不同的是，在8顆粒模型中基體斷裂主要發生在顆粒分布較密處，其他部分基體未發生斷裂，見圖6(c)。由圖6(d)～(f)可知，低溫條件下，基體破壞也出現在顆粒分布較密處，但由于低溫下基體變硬，導致基體斷裂裂尖較為銳利，整體斷裂面更加光滑平整。

綜上所述，低溫會對推進劑界面脫粘和基體斷裂的時機、位置和損傷程度造成影響。

3 結論

(1)通過建立三種不同顆粒尺寸、數量和分布的推進劑細觀模型，開發模擬基體開裂的VUMAT本構子程序，并采用Cohesive單元描述推進劑材料的顆粒/基體界面，可以較好地模擬常溫和低溫條件下其損傷破壞過程。

(2)低溫條件下，單顆粒推進劑脫濕首先發生在大顆粒表面，基體破壞沿顆粒赤道附近發生，裂紋有脫粘面向粘合劑基體擴展，并最終貫穿基體內部。多顆粒推進劑界面脫粘首先發生在大顆粒，其次小顆粒開始發生脫粘，而基體斷裂更容易發生在大顆粒附近，斷裂方向沒有明顯規律。外界低溫環境對推進劑界面脫粘和基體斷裂的時機、位置和損傷程度造成顯著影響，導致推進劑破壞形式更加復雜。

(3)由于本文建立的模型尺寸較小，且主要集中在現象分析上，所揭示出的推進劑低溫損傷規律未能量化，后續可對模型進行合理優化，擴大模型尺度，開展低溫條件下推進劑細觀損傷規律與宏觀力學性能相關性的研究。

[1] 劉著卿,李高春,邢耀國,等.復合固體推進劑細觀損傷掃描電鏡實驗及數值模擬[J].推進技術,2011,32(3):412-416.

[2] 袁嵩,湯衛紅,李高春.復合推進劑細觀失效機理分析[J].固體火箭技術,2006,29(1):48-51.

[3] 劉承武,陽建紅,鄧凱,等.基于Mori-Tanaka有限元法的粘彈復合推進劑非線性界面脫粘[J].推進技術,2011,32(2):225-229.

[4] Persson P O,Strang G.A simple mesh generator in matlab[J].SIAM Review,2004,46(2): 329-345.

[5] Persson P O.Mesh generation for implicit geometries[D].PhD thesis,Massachusetts Institute of Technology,2005.

[6] Watson D.Computing then-dimensional delaunay tessellation with application to voronoi polytopes[J].The Computer Journal,1981,24(2): 167-172.

[7] Matous K,Inglis H M,Gu X,el al.Multiscale modeling of solid propellants: from particle packing to failure[J].Composites Science and Technology,2007,67(8): 1694-1708.

[8] 郭素娟,康國政.金屬延性裂紋擴展過程的有限元模擬[C]//中國力學學會2009學術大會論文集,2009:200-205.

[9] 彭云.基于ABAQUS的三維裂紋動態擴展有限元模擬[J].南昌航空大學學報,2010,24(2):41-45.

[10] 趙玖玲,強洪夫.基于粘附功的復合推進劑AP/基體界面損傷宏細觀仿真 [J].固體火箭技術,2011,34(5):614-618.

[11] Bardenhagen S G,Brackbill J U,Sulsky D L.Shear deformation in granular material[J].Eleventh International Detonation Symposium,1998:547-555.

[12] 職世君,張建偉,張澤遠.復合固體推進劑細觀損傷形貌數值模擬[J].固體火箭技術,2015,38(2):239-244.