复合固体推进剂脱湿过程细观建模与损伤定量表征①

张 超，侯俊玲，李 群

(西安交通大学 航天航空学院 机械结构强度与振动国家重点实验室，西安 710049)

0 引言

复合固体推进剂作为一种含能材料，在航天航空领域得到了非常广泛的应用，是固体火箭发动机不可维修的最薄弱环节之一。固体推进剂的结构完整性直接影响着导弹武器的贮存和使用情况，必须在使用前就对其力学性能进行准确预测。早期关于固体推进剂力学行为的研究，主要是从连续介质力学角度基于唯象理论建立其本构关系[1-2]，这些理论均将固体推进剂作为连续均匀介质来处理。然而，在细观尺度上，固体推进剂是一种高填充比颗粒增强复合粘弹性材料，主要由以高分子聚合物为基体的粘合剂HTPB和大量的AP颗粒及金属燃料(Al)颗粒组成，其力学行为主要由固体填料的体积分数、粒径大小、粘合剂基体的粘弹性质以及粘合剂基体与固体填料之间的界面性质所决定。因此，从细观角度出发研究固体推进剂的力学行为具有重要意义。

要从细观角度出发研究固体推进剂的力学行为，首先要建立能够有效表征固体推进剂细观结构的颗粒填充模型。在已发表的文献中，关于颗粒复合材料填充模型的建立，主要分为两类：第一类是通过实验观测的内部微结构重构，建立材料的真实细观结构[3]，其能够准确反映材料的真实细观结构，但对试验条件和技术手段等均有较高的要求，在相关图像和数据的处理过程中工作量巨大，成本十分高昂。关于颗粒填充复合材料细观结构建模，学者们大多采用第二类方法，即通过数值模拟仿真方法，对材料典型细观结构进行建模[4-8]。例如，基于蒙特卡罗[9]算法的随机模拟算法和基于分子动力学[10]算法的碰撞模拟方法建立二维颗粒填充模型。近年来，为建立高体积分数三维颗粒填充模型，学者们进行了大量的研究，并取得了不错的成果。例如，Zhang等[4]提出的下落堆积法、Yu等[6]提出的空间压缩法和传统的随机序列吸附法[7]。但是，采用这些方法，在颗粒填充过程中，当颗粒体积分数高于50%时，往往会遇到颗粒投放重叠检测迭代次数大、建模时间过长甚至失效的问题，无法达到材料配方所要求的体积分数。因此，亟需要开发一种建立高体积分数三维颗粒填充模型的算法，以便开展高体积分数颗粒填充复合材料的数值仿真研究。

此外，大量研究结果表明[11-22]，颗粒/基体界面脱湿是固体推进剂在外界载荷作用下的主要失效形式。针对推进剂颗粒/基体界面脱湿问题，国内外学者们开展了大量工作，探究了细观结构损伤演化规律及其对宏观力学行为的影响。袁嵩等[23]建立 了简化的单颗粒和多颗粒细观结构模型，研究了颗粒与粘合剂基体在粘结完好和存在脱粘两种情况下固体推进剂内部的应力分布。Matous K等[24]建立了推进剂代表体积单元的细观模型，数值模拟了颗粒/基体界面损伤的产生及发展。李高春等[25]综合考虑固体推进剂细观损伤的特点，研究了推进剂的界面脱湿过程及其对宏观力学性能的影响。张建伟[26]、职世君[27]等研究了固体填料的体积分数、颗粒位置的随机分布、粒径大小对推进剂细观损伤及宏观力学性能的影响。目前，研究大多集中于二维模型或三维空间低颗粒体积分数模型，开展高体积分数的固体推进剂三维细观数值模拟，能更好为地反映受载条件下固体推进剂的细观损伤特性，进而为固体推进剂的使用与配方设计提供指导。

本文基于随机序列吸附法、热膨胀原理及分步填充的思想，实现了高体积分数颗粒填充，并通过几何合并的方式，加入粘合剂基体和颗粒/基体界面，建立固体推进剂三维细观结构模型。引入内聚力界面单元，基于双线性损伤内聚力模型，考虑粘合剂基体的粘弹特性，研究不同应变率加载下固体推进剂材料内部的界面脱湿行为，并基于内聚力模型定义损伤变量，对其细观损伤进行了定量表征。

1 三维细观结构模型

本文以HTPB推进剂[11]为研究对象，其主要成分包括HTPB基体、AP颗粒、Al颗粒以及其他助剂，其基本组分的比例见表1。HTPB推进剂是一个典型的颗粒填充复合材料，其颗粒的粒径分布情况如图1所示[10]。

图1 固体推进剂颗粒粒径分布Fig.1 Particle size distribution of solid propellants

本文基于随机序列吸附法、热膨胀原理及分步填充的思想，实现了高体积分数颗粒填充，详细颗粒填充流程如图2所示。在三维细观结构颗粒填充过程中，首先将颗粒粒径等比例缩小，基于随机序列吸附法，运用Matlab编程语言，向胞元内随机投放小颗粒，记录小颗粒的粒径、热膨胀系数和球心坐标。然后，使用Python脚本语言进行Abaqus二次开发建立初始模型，通过热膨胀分析提高颗粒粒径。建模过程中，为提高效率，每一轮填充的颗粒数不宜过多，需要经过多步填充和热膨胀过程，才能够建立高体积分数三维细观结构模型。需要注意的是，若颗粒初始粒径过小，会造成热膨胀分析过程中有限元网格畸变和收敛的困难。因此，小颗粒粒径的选取应以颗粒初始体积分数不小于实际体积分数的33%为宜。

图2 三维细观结构颗粒填充流程图Fig.2 Particle filling flowchart of three dimensional mesostructure

三维细观结构建模过程详细说明如下：

(1)第一轮颗粒填充：考虑到大颗粒对固体推进剂的力学行为影响较大，首先基于随机序列吸附法，将大颗粒以真实粒径由大到小依次填充进胞元内，直至单个颗粒的重叠检测次数等于预设值，颗粒填充效果如图3(a)所示，颗粒数为6，颗粒体积分数为31.49%。因为大颗粒是以真实粒径填充的，因此不需要进行热膨胀分析。

(2)第二轮颗粒填充：在第一轮填充颗粒的基础上，进行第二轮颗粒填充时，为保证热膨胀分析的收敛性并减小计算量，第二轮填充进12个颗粒，此时胞元内共有18个颗粒，体积分数达到40.72%，第二轮颗粒填充初始状态及网格划分如图3(b)所示。特别需要注意的是，在热膨胀分析过程中，装配体不包含粘合剂基体和颗粒/基体界面部件，在颗粒填充完成后，通过几何合并的方式将粘合剂基体和颗粒/基体界面加入装配体即可，同时在模型中引入通用接触以防止颗粒重叠。

(3)第二轮填充颗粒热膨胀分析：在热膨胀分析过程中，建立6个薄壁面并完全固定，以保证颗粒在热膨胀分析过程中始终处于胞元所在的空间内。颗粒热膨胀分析后的应变和位移云图分别如图3(c)、(d)所示，热膨胀分析后颗粒体积分数为45.05%，较第二轮填充的初始状态提高了4.33%。由图3(c) 、(d)可见，在颗粒热膨胀分析过程中，随着温度逐渐升高，第二轮填充的颗粒粒径不断增大，并可通过颗粒与薄壁面、颗粒与颗粒之间的相互作用，使颗粒产生刚性位移，将每一个颗粒重新分配到合适的位置。但当其中某些颗粒不存在理论膨胀空间时，即使在其他颗粒仍有膨胀空间的情况下，热膨胀分析过程也会出现不收敛的情况，导致热膨胀分析终止。

(4)第二轮部分颗粒热膨胀分析：基于Python脚本语言进行Abaqus二次开发，提取所有颗粒中心节点的坐标，并根据热膨胀系数和热膨胀分析时间计算颗粒当前粒径，重新建立模型。同时，将第二轮没有理论膨胀空间的颗粒热膨胀系数置为零，仅对有膨胀空间的颗粒进行热膨胀分析，再次提高颗粒体积分数，图3(e)、(f)分别为仅对某一个颗粒进行热膨胀分析的应变和位移云图。同理，对其他有膨胀空间的颗粒进行热膨胀分析，直至所有颗粒都不存在理论膨胀空间时，第二轮热膨胀分析结束，此时颗粒体积分数为47.95%，再次提高2.9%。

(a)Particle filling effect of the first round (b)Initial state and grid division of particle filling

(c)Strain distribution of particle thermal expansion analysis (d)Displacement distribution of particle thermal expansion

(5)重复步骤(2)～(4)工作，将更多颗粒填充进胞元内，最终可得到胞元尺寸为0.613 3 mm3，颗粒数为150，颗粒体积分数为63.8%的三维细观结构模型，如图4所示。

图4 颗粒填充细观模型(体积分数为63.8%)Fig.4 Particle filled mesoscopic model (volume fraction is 63.8%)

2 界面脱湿数值模拟

2.1 推进剂组分参数

HTPB推进剂的主要组分如表1所示，在数值模拟过程中，主要参数为固体填充颗粒的弹性模量和泊松比，以及粘合剂基体的松弛模量。由于AP颗粒的弹性模量远大于粘合剂基体的弹性模量，在相同载荷的作用下，相较于基体的变形，AP颗粒的变形可忽略不计，故假设AP颗粒为弹性体，其弹性模量EAP=32 450 MPa，泊松比vAP=0.143 3[3]。同时，为了进一步简化计算模型，将Al颗粒对固体推进剂力学性能的影响等效到基体中考虑，文中Al颗粒和粘合剂基体混合物的初始弹性模量设为7.39 MPa[3]。

HTPB推进剂粘合剂基体具备粘弹性材料的基本性质，是导致推进剂具有粘弹性的根本原因，粘合剂基体松弛模量依据文献[28]中推进剂应力松弛实验结果选取，对松弛曲线用Prony级数的形式进行拟合，其拟合的表达式为

(1)

式中t为时间；E∞为平衡模量，是t→∞时模量的稳态模量；Ei和τi分别为第i个Maxwell单元的模量和松弛时间。

所得Prony级数表达式的各项系数见表2，进而得到粘合剂基体的松弛模量应力应变关系为

(2)

表2 HTPB推进剂松弛模量Prony级数表达式系数Table 2 Prony series expression coefficient of relaxation modulus of HTPB propellant

2.2 内聚力模型及参数

图5 双线性内聚力模型Fig.5 Bilinear cohesion model

本研究中，损伤初始阶段采用最大应力准则，其控制方程为

(3)

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

(4)

其中，D为损伤因子，定义为

(5)

双线性内聚力模型主要包括3个参数：初始线性模量、最大应力值和最终开裂位移值。本研究中，参考文献[25]中双线性内聚力模型数据，初始线性模量取为基体模量的10倍；界面最大应力取为0.73 MPa；最终开裂位移值取为37 μm。

2.3 有限元网格及边界条件

在有限元网格划分过程中，AP颗粒采用结构化网格划分技术，单元类型为C3D8R，共20 424个单元；界面采用结构化网格划分技术，单元类型为COH3D8，共10 320个单元；由于颗粒的存在，导致基体结构非常复杂，采用自由网格划分技术，单元类型为C3D10，共83 676个单元。AP颗粒和界面网格划分如图6(a)所示，粘合剂基体网格划分如图6(b)所示。

(a)Mesh of AP particles and interface

(b)Mesh of adhesive matrix图6 细观模型网格划分Fig.6 Mesh of mesoscopic model

为模拟固体推进剂在均布位移载荷条件下的力学行为，采用如下边界条件：面EFGH保持固定，约束面ADHE、面BCGF、面ABFE和面DCGH的法向位移，在面ABCD沿x轴正方向施加均布位移载荷，如图7所示。同时，在模型中引入通用接触以防止界面脱湿后导致的界面之间相互渗透。

图7 细观模型边界条件Fig.7 Boundary conditions of mesoscopic model

3 结果与分析

3.1 脱湿过程数值模拟结果

基于双线性损伤内聚力模型，考虑粘合剂基体的粘弹特性，对所建三维细观结构模型进行不同应变率下的单轴拉伸数值模拟，HTPB推进剂在不同应变率下的单轴拉伸应力-应变曲线如图8所示。由图8可见，固体推进剂在不同应变率载荷下表现出不同的力学行为，应变率越高，应力-应变曲线的斜率越高，材料表现出更高的刚度。图9和图10分别给出应变率为3.33×10-3s-1，应变分别为3%、8%和13%时，对应的固体推进剂内部Mises应力和刚度衰减率SDEG损伤云图，其中刚度衰减率SDEG是Abaqus后处理结果中表征材料刚度衰减率的变量，当SDEG=0时表示材料没有损伤，SDEG=1时表示材料已经完全损坏。

图8 HTPB推进剂单轴拉伸应力-应变曲线Fig.8 Uniaxial tensile stress-strain curves of HTPB propellant

图9(a)和图10(a)分别为应变为3%时的Mises应力和SDEG损伤云图。在应变为3%时，应力-应变呈典型线性阶段(如图8所示)。此时，在均布位移载荷的作用下，由于颗粒与基体、颗粒与颗粒之间的相互作用，在大颗粒的极区位置产生很高的局部应力场，但颗粒与基体仍然处于弹性变化范围，材料内部尚未出现界面脱湿损伤。

(a)Strain is 3%

(b)Strain is 8%

(c)Strain is 13%图9 固体推进剂Mises应力云图Fig.9 Mises stress distribution in solid propellants

随着应变增加到8%，如图9(b)和图10(b)所示，在大颗粒的极区位置应力集中明显，出现了一定的界面损伤。这是因为颗粒与基体的材料属性不同，颗粒模量远大于基体模量，在作用力相同情况下，颗粒变形小于基体变形，使得界面成为最薄弱的部位，故界面最先出现损伤，且因界面损伤导致颗粒承载能力下降，材料等效模量降低。从图8所示的应力-应变曲线也可看出，在8%应变处，应力应变呈非线性关系，表现出一定程度的软化，这表明材料内部出现了一定的界面损伤。

由图9(c)和图10(c)可知，随着加载继续进行到应变为13%时，在比较密集的颗粒附近，也出现了大量的界面损伤，且在颗粒界面损伤严重的区域，Mises应力进一步降低。这说明，由于界面损伤导致大颗粒承载能力进一步下降，材料等效模量迅速减小。对比图8所示的应力-应变曲线，在应变为13%时，应力-应变呈强非线性关系，且软化加剧，表明材料内部出现了大量的界面损伤。

(a)Strain is 3%

(b)Strain is 8%

(c)Strain is 13%%图10 固体推进剂SDEG损伤云图Fig.10 SDEG distribution in solid propellants

3.2 基于界面刚度衰减率的脱湿损伤定量表征

为了更好表征固体推进剂材料的界面脱湿行为，本文提出一种基于界面刚度衰减率的损伤定量表征方法。定义固体推进剂材料内部脱湿面积总和(DA)为

(6)

式中ED(i)为第i个界面单元的损伤体积；d为界面单元的厚度；N为所有界面单元的个数。

ED(i)定义为

(7)

式中SDEG(j)为第i个界面单元第j个高斯积分点的刚度衰减率；V(j)为第i个界面单元第j个高斯积分点的几何体积；n为第i个界面单元的高斯积分点个数。

选取颗粒体积分数为63.8%的固体推进剂细观模型，图11展示了不同应变率加载情况下，脱湿面积DA随载荷的变化情况。

(a)Dewetting area-strain

(b)Enlarged view partially图11 损伤-应变曲线Fig.11 Damage-strain curves

从图11(a)和12(a)中可看出，在低应变载荷情况下，损伤面积DA=0，表明材料处于弹性范围内，内部没有出现界面损伤。而当应变超过6%附近时，随时载荷的增大，脱湿面积DA的值显著增大，这表明材料内部出现界面损伤。图11(a)所示的脱湿面积DA的拐点，可认为是材料的初始脱湿点，其应变载荷值大约在6%附近，如图11(b)和12(b)所示。随着载荷增大，DA值逐渐急剧增大，表明其损伤状态严重，如图11(c)和12(c)所示。

此外，从图11(b)给出大应变载荷的损伤面积局部放大图中可以看出，在相同应变的情况下，应变率越高，材料内部脱湿面积越大。这是由于粘合剂基体具有与时间相关的粘弹特性，在相同应变的情况下，应变率越高，加载时间越短，材料内部应力越大(如图8所示)，更易造成界面损伤。

(a)Strain is 3%

(b)Strain is 6%

(c)Strain is 13%图12 刚度衰减率SDEG损伤云图 (strain rate=3.33×10-3 s-1)Fig.12 SDEG damage distribution (strain rate=3.33×10-3 s-1)

4 结论

(1)基于随机序列吸附法、热膨胀原理和分步填充的思想，实现了高体积分数颗粒填充，并通过几何合并的方式，加入粘合剂基体和颗粒/基体界面，建立了高体积分数固体推进剂三维细观结构模型。所建模型能有效表征固体推进剂的三维细观结构。

(2)基于双线性损伤内聚力模型，在颗粒/基体界面嵌入内聚力单元，考虑粘合剂基体与时间相关的粘弹特性，对不同应变率下固体推进剂内部脱湿过程进行了数值模拟。结果表明，由于颗粒与基体材料属性相差较大，脱湿首先出现在大颗粒及颗粒比较密集的区域，界面损伤导致材料承载能力下降；且应变率越高，材料内部越容易出现损伤。

(3)使用Python脚本语言进行Abaqus二次开发，提取所有内聚力界面单元高斯积分点的几何体积和刚度衰减率SDEG数据，定义固体推进剂材料内部脱湿面积，可准确捕捉固体推进剂颗粒/基体界面的脱湿点，并对其细观损伤进行定量表征。