固体推进剂低温细观损伤仿真研究

周红梅，袁 军，赖建伟，袁 嵩

(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.