考虑颗粒形状的复合固体推进剂细观损伤分析①

侯宇菲，许进升，陈 雄，周长省，朱 亮

(南京理工大学 机械工程学院，南京 210094)

0 引言

复合固体推进剂材料的力学性能在固体火箭发动机药柱的结构完整性中起着重要作用。若想达到固体推进剂在工作过程中的力学指标，必须提前对推进剂的力学性能进行预测，保证固体推进剂的结构完整性。在过去的几十年中，许多研究人员从连续介质力学角度对固体推进剂本构关系进行了大量的探索与研究[1]，这些理论研究均将固体推进剂视为一种连续均匀的介质。然而，复合固体推进剂是一种多颗粒填充的非均匀含能材料，上述方法不能从本质上解释造成固体推进剂力学性能非线性的原因。因此，采用细观方法对复合固体推进剂的力学性能进行研究成为一种有效途径。

采用扫描电镜对丁羟推进剂的端口形貌进行观察，发现颗粒/粘合剂界面脱湿是造成推进剂力学行为非线性的主要原因[2]。由于现阶段试验条件的限制，高填充比颗粒复合材料的细观实验还难以开展，研究者通常采用数值仿真方法进行颗粒/基体界面脱湿的相关力学行为研究分析。彭威[3]使用轴对称圆柱体模型进行数值模拟，发现应力集中现象出现在颗粒/基体界面的极区。Matous[4-5]通过自主开发的软件生成了推进剂的代表性体积单元，研究了不同加载情况下推进剂的细观损伤机理及过程。为进一步探究颗粒/基体界面的力学行为，常武军[6]采用双线性粘聚区模型对颗粒/基体界面脱湿进行数值模拟，得到了不同界面特性及颗粒尺寸对宏观力学的影响规律。封涛[7]在颗粒/基体界面采用Surface-based Cohesive方法对推进剂颗粒/基体脱湿进行了数值模拟。袁崇[8]在研究颗粒/基体界面脱湿时，通过Mori-Tanaka方法得出填充颗粒粒径、填充颗粒体积分数及颗粒界面间最大粘接应力对推进剂的力学性能均有明显影响，并在外载荷作用下，由于颗粒与基体的材料属性不同，推进剂内部的应力分布是不均匀的，当颗粒与基体界面发生脱湿时，推进剂内部的应力将重新分布。韩龙[9]通过推进剂的力学试验与细观数值模拟发现推进剂的力学性能与基体胶片的松弛行为有关，与颗粒的随机分布无关。Zhi[10]着重探讨了颗粒/基体界面脱湿对松弛模量的影响，结果表明推进剂损伤程度越高，其松弛模量越低。刘新国[11]分别从观测及表征方法、理论模型及数值模拟三方面对颗粒/基体脱湿行为进行相关描述，认为解决颗粒/基体脱湿是提高推进剂宏观力学性能的关键。

本文将从细观角度出发，采用分子动力学方法建立圆形颗粒填充模型与多边形颗粒填充模型，在颗粒与基体的界面处使用Python脚本语言嵌入零厚度Cohesive Element，分别使用双线性内聚力模型和指数型内聚力模型对颗粒填充模型进行不同速度载荷下的数值模拟，通过数值仿真计算得到的脱湿点与试验结果对比，验证模型准确性。

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组拉伸速率下的工程应力应变曲线。

Al颗粒作为加速HTPB热分解的添加剂，通常不会对界面的脱湿造成影响。根据文献[10]对复合基体进行应力松弛试验，用Prony级数对松弛曲线进行拟合，其拟合的表达式为

(1)

五阶Prony级数松弛模量主曲线拟合结果为

(2)

从而得到基体胶片的松弛模量应力应变关系：

(3)

图1 HTPB推进剂单轴拉伸应力-应变曲线

2 细观模型建立

2.1 颗粒填充模型

使用扫描电镜对HTPB断面形貌进行观察，扫描倍数为100，如图2所示。

通常HTPB颗粒填充系数较高，颗粒均匀地分布在基体内部，小颗粒散布于大颗粒周围，所有颗粒的形状均类似于球形的多边形，形状及其不规则的颗粒极为少见。从图2可观察到已脱湿颗粒和颗粒脱湿后留下的凹坑，说明在外载荷的作用下，颗粒/基体脱湿率先导致了推进剂的结构破坏。因此，从细观角度出发对推进剂的失效破坏进行分析成为有效途径。

本文采用分子动力学算法生成了随机分布的圆形颗粒填充模型。多边形颗粒是在圆形颗粒基础上使用VB语言生成的。根据圆形颗粒半径及圆心坐标位置确定多边形颗粒的位置及大小：

(4)

式中xp与yp为多边形颗粒的x轴与y轴坐标；xc与yc为圆形颗粒的x轴与y轴坐标；r为圆形颗粒的半径；n为生成随机多边形的边数，本文选用十四边形，即n=14；i为多边形的第i条边。

按照表1给出的HTPB各颗粒组分含量，生成了如图3所示的随机颗粒填充模型。

图2 HTPB的断口形貌

(a)圆形颗粒

(b)多边形颗粒

由于Al和AP粒径相差较大，会降低细观模型的网格质量，并影响数值仿真结果。为提高网格质量，本文将Al与基体粘结剂视为一体，其初始模量为1.232 8 MPa。根据文献[5]获得AP颗粒的相关参数，弹性模量为32 447 MPa，泊松比为0.143 3 MPa。

2.2 Cohesive element二次开发

Cohesive element有几何厚度和计算厚度。先前使用VB语言在基体/颗粒之间嵌入了有限厚度Cohesive element进行数值仿真，然而由于粘结层的厚度很小，大约为0.000 1 mm，这种方法通常会造成网格质量下降，进而影响仿真结果收敛性。为了解决这个问题，现使用Python脚本语言开发零厚度Cohesive element。下面详述零厚度Cohesive element的生成：

(1)将建立好的随机颗粒填充模型导入ABAQUS中，分别建立颗粒单元集合与基体单元集合，采用四节点的应变单元对所有单元进行网格划分，并以ABAQUS.inp文件格式输出。

(2)找出颗粒单元与基体单元的公共点，命名为界面集合。

(3)对所有公共节点进行拆分，拆分后节点编号进行重新排列。需注意的是，粘结单元的编号必须使用逆时针编号，否则ABAQUS不能识别。其网格类型为COH2D4。

(4)将新生成的文件以ABAQUS.inp的文件格式重新导入ABAQUS中，进行数值模拟。

2.2.1 双线性内聚力模型

双线性内聚力模型广泛应用于材料与结构的开裂破坏研究中，最早在研究脆性材料的断裂时被提出[12]，后来将内聚力模型与有限元方法相结合，实现数值模拟的材料断裂[13]。典型的双线型内聚力模型见图4。

双线性内聚力模型张力位移关系的控制方程为

(5)

(6)

(a)法向

(b)切向

(7)

在双线性张力位移关系中，除了最大应力值和临界能值作为参数必须给出外，还需给界面临界张开位移或应力上升阶段斜率K。

2.2.2 指数型内聚力模型

ABAQUS提供了三种内聚力子程序开发方法：自定义界面本构子程序、自定义单元子程序及自定义材料子程序。自定义材料子程序方法包括UMAT和VUMAT两种方法。其中，UMAT不能和零厚度的粘结单元共同使用。因此，本文选用VUMAT来实现内聚力模型的张力位移关系。VUMAT通过更新状态变量矩阵与材料常数矩阵PROPS(*)实现材料的本构关系。状态变量矩阵为StateOld(*)和StateNew(*)。应力张量矩阵为StressOld(*)和StressNew(*)。StateOld(*)和StressOld(*)代表上一增量步的状态变量值与应力张量；StateNew(*)和StressNew(*)代表增量步结束时的状态变量值与应力张量。例如，在第n个增量步中有：

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

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

(8)

式中n代表第n个增量步；i代表第i个状态变量或应力张量。

如图5所示，指数内聚力模型是连续性方程，其应力变化是连续的。应力在减小的过程中呈非线性地渐进为零。指数型内聚力模型选用断裂能作为判断单元损伤的条件：

(9)

其中，Δn、Δt为界面在法向和切向的位移值，应力最大点对应的位移值为特征位移值，分别为δn和δt，φn为纯法向开裂时所需的断裂能; 参数q、r分别为

(10)

(11)

(a)法向

(b)切向

2.3 内聚力模型的参数获取

内聚力模型参数的准确性直接影响数值模拟结果的正确性。在目前的工程应用和理论分析中，内聚力参数是很难直接由实验获得的，大多数取值来源于文献和经验，但是这种方法得到的模型参数往往不精确。本文通过Hooke-Jeeves的反演方法[14]得到模型参数，其基本原理是通过实验结果拟合出模型参数，如图6所示。在拟合过程中通过不断改变参数取值来调整仿真结果，使得仿真结果与实验结果的误差小于规定值。

图6 界面反演流程图

双线性内聚力模型需确定3个参量，分别为界面刚度K、内聚强度σmax和破坏位移δf。首先在模型中设定3个参数的初始值：界面刚度K=200 MPa/mm，内聚强度σmax=0.2 MPa，破坏位移δf=0.01 mm。依据仿真应力应变曲线minR得到应力值，并与试验曲线获得应力值进行比较。通过目标函数minR得知，试验值与仿真值的误差大小为

(12)

通过上式对函数值进行迭代，不断更新界面刚度、内聚强度和破坏位移，使得minR小于规定的容差Rlim。当反演结果与实际实验结果的minR符合规定值时，认为反演参数具备可行性。指数内聚力模型的参数反演和双线性内聚力模型参数获得途径一致。指数内聚力模型只需确定内聚强度σmax和破坏位移δf，设定初始内聚力强度σmax=0.2 MPa和破坏位移δf=0.01 mm。反演后得到的参数如表2所示。

表2 双线性内聚力模型及指数型内聚力模型界面参数

3 结果与讨论

3.1 仿真结果

使用Abaqus软件对圆形颗粒与多边形颗粒模型进行5组加载速度下的仿真计算，由于图形较多，只给出采用指数型内聚力模型的圆形颗粒填充模型与多边形颗粒填充模型在拉伸速度为1 mm/min时的Mises应力云图，如图7和图8所示。

(a)应变3% (b)应变10% (c)应变15%

(a)应变3% (b)应变10% (c)应变15%

图7(a)和图8(a)为应变为3%时的Mises应力云图变化。此时，颗粒与基体均处于弹性变化范围内，其形状没有发生大的改变，只是云图的颜色发生了改变，说明颗粒与基体的应力分布不均匀，颗粒与基体之间正在进行力的传递。

图7(b)和图8(b)为应变为10%时Mises应力云图变化。在大粒径圆形颗粒与多边形颗粒的极区位置，均发生了界面位移分离。这是因为颗粒与基体的材料属性不同，其模量也不同。颗粒模量大于基体模量。在作用力相同的情况下，颗粒变形小于基体变形，使得界面成为最薄弱的部位，故发生位移分离现象。

当应变达到15%时，如图7(c)和图8(c)所示，在颗粒聚集较密集的地方也发生了界面位移分离。已经脱湿的颗粒从孔洞中裸露出来，基体明显伸长。此时，颗粒与基体的界面不再承担外载荷作用，即到达了颗粒的最大脱湿点。

3.2 结果分析

通过图7和图8发现，在外载荷的作用下，颗粒与基体之间的界面会产生分离现象，进而造成颗粒/基体界面脱湿和孔洞的生成。为进一步探究颗粒与基体界面分离对推进剂力学性能的影响并验证模型准确性，将对实验获得的初始模量及脱湿点模量与仿真模型得出的初始模量及脱湿点模量进行比较，以验证各模型的预测能力和确定适用范围。

如图9所示，随着加载速率的提高，HTPB固体推进剂的初始模量增大，表现出显著的粘弹性效应。使用双线性内聚力圆形颗粒填充模型和多边形颗粒填充模型的初始模量预测结果明显大于使用指数型内聚力圆形颗粒填充模型与多边形颗粒填充模型。这是因为指数型内聚力模型与双线性内聚力模型在计算法则不同，双线性内聚力模型更适用于描述脆性材料的本构，而指数型内聚力模型适合描述非线性材料的本构。因此，指数型内聚力模型的预测值与试验结果误差更小。

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

当应力应变曲线的斜率明显降低时，即发生了颗粒的脱湿，这一点称为颗粒/基体界面脱湿点。如图10所示，和初始模量的预测结果类似，指数型内聚力模型比双线性内聚力模型的预测结果更为准确。因为两种模型的损伤变量形式不同，指数型内聚力模型的损伤模量要高于双线性内聚力模型的损伤模量。所以，指数型内聚力模型的损伤速度较快，使得界面能更快软化，得到脱湿点模量值更加准确。多边形颗粒填充模型较圆形颗粒填充模型得到脱湿点模量更准确，这是因为AP颗粒的形状并不是真的圆形，为计算方便，将AP颗粒简化成圆形进行数值模拟。因此，多边形颗粒填充模型比圆形颗粒填充模型更接近推进剂内部的真实情况。

图10 不同拉伸速率下的脱湿点模量

从图11可看出，拉伸速率越大，模量下降率越大，并且不同加载速率的脱湿点模量近等相似，说明HTPB固体推进剂这种材料在高速率外载荷的作用下，更容易出现损伤。为保证推进剂的力学性能，下一步应探究在高速率的工作载荷下，延缓推进剂损伤的措施。

图11 模量下降率

数值模拟得出的模量值均略高于试验值，这是因为本文所建立的模型均是基于二维平面应变，而推进剂真实形态是三维状态的，且本文在建立模型时，认为模型的内部是不存在微小缺陷的。在实际情况中，推进剂的内部存在微孔洞等缺陷，这会造成应力值的下降。以上原因造成了模型的计算值大于试验值。

4 结论

(1)分子动力学方法生成的多边形颗粒填充模型比圆形颗粒填充模型能更好地描述推进剂细观结构，更准确地判断损伤位置。

(2)由于目前细观试验手段的局限性，引入Hook-Jeeves优化反演方法，获得了更准确的内聚力模型参数。通过不同损伤形式内聚力模型的比较，指数内聚力模型更适用于延展性材料的损伤，能更好地预测推进剂的力学性能，使得仿真结果更接近于试验值。

(3)通过5组加载速率的试验结果与仿真结果比较，发现推进剂的损伤和加载速率相关，载荷加载速率越大，推进剂受损程度越大，模量下降率越高。

(4)下一步将开发三维推进剂颗粒填充模型，以提高模型的准确性。