推进剂“脱湿”损伤研究的内聚力单元方法

崔辉如，吕 轩，许玉荣

(1. 陆军工程大学 国防工程学院, 江苏 南京 210007； 2. 湖北航天技术研究院总体设计所， 湖北 武汉 430040)

近年来，细观力学的发展为研究推进剂非线性力学性能提供了新的研究手段。在此基础上，可以结合数值损伤模型，开展推进剂细观力学性能研究，揭示推进剂“脱湿”损伤机理，对推进剂非线性力学性能的预示具有较大的推动意义[1-7]。

通过颗粒填充手段模拟生成的推进剂细观结构为推进剂松弛模量的预测[4, 8]、推进剂细观损伤的数值模拟研究[5]提供了便利条件。推进剂细观结构中，颗粒和基体之间黏接界面的破坏是“脱湿”产生的主要原因。Zhang[9]以及Han[10]等建立了固体推进剂的细观力学模型，对大变形下的推进剂内部界面“脱湿”现象进行模拟，采用内聚单元对黏接界面进行建模。为了模拟不同固体含量的推进剂中细观损伤的发展历程，职世君等[6]通过分子动力学的方法得到了推进剂的颗粒装填模型，其基于内聚接触手段结合双线性内聚力本构关系，建立了AP颗粒和基体之间的损伤模型；基于有限元方法，对比分析了不同固体含量对细观损伤形态的影响并指出高体分比状态下颗粒排列非常致密，界面损伤的相互影响较大是加快固体推进剂产生损伤至形成宏观裂纹的过程的主要因素；此外，研究还表明，内聚强度越小，复合固体推进剂的损伤起始时间越早，损伤程度越大，相应的拉伸应力越小。韩龙等[2]开展了复合固体推进剂细观界面率相关参数的反演分析并揭示了黏接界面性能随加载速率的变化规律。封涛等[1]采用双线性内聚力模型来描述颗粒与基体黏接单元的力学响应，深入研究了细观结构对端羟基聚丁二烯(hydroxyl terminated polybutadiene)推进剂力学性能的影响。周红梅等[11]同样采用了内聚力单元的方式，开展了推进剂细观结构在常温以及低温下的损伤机理研究，研究表明低温下的推进剂细观结构破坏形式更为复杂。此外，还有研究[12]指出，在细观结构下，推进剂一般首先在极区产生脱黏，颗粒的大小对推进剂的力学性能有较大影响。内聚力单元方法是推进剂细观“脱湿”损伤研究的重要途径之一[13]。在推进剂细观几何的处理、内聚力单元模型的构建以及边界条件的设计方面，上述文献中的处理方式还不是特别全面和完善。通过内聚力单元方式进行推进剂“脱湿”损伤机理的研究还有待进一步深入。

为了进一步研究推进剂细观“脱湿”损伤机理，本文在合理的代表体积单元尺寸的基础上，采用内聚力单元方法对颗粒和基体之间的黏接界面进行建模。结合周期几何和周期边界条件，采用单轴拉伸和剪切试验的手段，开展推进剂细观“脱湿”损伤的机理研究，并对不同推进剂体分比、应变速率以及内聚强度进行“脱湿”损伤力学行为影响分析。

1 推进剂细观有限元模型

1.1 细观几何模型的构建

推进剂是典型的高填充型复合材料，颗粒的体积分数往往会达到60%，甚至更高。目前分子动力学方法是行之有效的方法，Kochevets[14]、申柳雷[15]的博士学位论文均给出了利用分子动力学进行推进剂颗粒装填模型建立的具体细节。

为了描述一个真实的细观结构的力学响应，细观模型必须要满足周期几何以及周期边界条件的要求。在进行颗粒装填时，如果发生颗粒的“压线”事件，那么对应几何边(或者三维模型中的面)上就要有相应的颗粒，才能满足周期几何的限制条件。此外，在进行颗粒装填之前，颗粒的数量以及相应的几何尺寸是按照材料的配方计算出来的固定的值，并不能因为装填的过程而改变。为了解决上述的两个问题，以模型的中心构建直角坐标系，第i个颗粒的坐标(颗粒边界顶点的坐标)为ηi(xi,yi,zi)，其中无量纲的坐标范围为[-1,1]。上述问题可以转化为

(1)

通过这样的手段，可以使得所有“压线”的颗粒都在模型的对应几何边界上有对称的颗粒，尽管颗粒的数量比初始计算的数量增加，但这并不影响总的装填分数，如图1所示。

图1 颗粒发生“压线”事件以及周期几何的处理Fig.1 "Line pressing" event of particles and periodic geometry processing

1.2 推进剂细观“脱湿”行为仿真

推进剂细观“脱湿”损伤仿真的核心是颗粒和基体界面分离行为的模拟。应用比较成功的主要有内聚接触模型和内聚力单元模型两种。内聚接触是ABAQUS自带的功能，其本质是对接触行为的模拟。当然，这种接触是一种广义的“接触”，既可以描述界面之间的相互挤压行为，又可以描述界面之间的相互分离行为。内聚力单元会将单元上分布的内聚力等效到单元的节点上，代入总体刚度阵的计算。然而，内聚接触是分布在界面上的一连串的离散的弹簧，被黏接的物体是通过一系列离散的弹簧相连接，至于没有弹簧连接的地方，是不存在任何相互作用的；内聚力单元则是界面上分布的连续的弹簧，被黏接物之间没有任何的间隙，黏接界面之间任何微小的分离都将被捕捉，并被代入单元节点力的计算中。正因为两种不同的建模理念，内聚接触方法会比内聚力单元方法更容易收敛，但是精度会大大降低，甚至会导致错误的结果，毕竟分段式的连接并不符合黏接的实际情况。因此，本文将采用更符合实际的内聚力单元进行界面损伤的建模。图2(b)给出了图2(a)中颗粒的有限元模型以及界面上分布的零厚度内聚力单元，可以发现，界面上分布的内聚力单元和颗粒单元的边界是完全吻合的。

图2 颗粒网格以及零厚度内聚力单元分布Fig.2 Distribution of mesh of particles and zero thickness cohesive elements

PPR(Park-Paulino-Roesler)内聚力模型是考虑物理断裂参数以及连续断裂边界的势相关内聚力模型[16]，模型的势函数表达式为

(2)

其中：Δn和Δt分别为法向和切向分离位移；Γn和Γt为能量常数；m和n为无量纲参数；〈·〉为MacAuley算子；αk和βk为模型的形状参数，分别控制法向和切向的内聚力和分离位移曲线下降段的凹凸特性；δn和δt为模型的法向和切向的临界位移；φt和φn分别代表切向和法向内聚能，其表达式为

(3)

本文中，黏接界面初始内聚强度σmax和τmax为0.2 MPa、初始临界位移δn和δt为50 μm、初始无量纲参数m和n为0.01、形状参数αk和βk为2.10。此外，颗粒采用线弹性材料属性，基体采用线黏弹性材料属性。

1.3 单元平均应力应变的后处理方法

商业有限元软件并不能直接给出模型的等效应力应变关系。为了便于进一步分析模型的力学响应特性，需要对计算结果进行处理。对于三角形网格模型，其平均应力应变关系有如下表达式

(4)

(5)

考虑到网格数量达到一定的规模时，网格对计算结果的影响差异不大，进行有限元的分析过程中，节点的位移是计算出来的直接的结果。但是考虑到部分节点既属于基体又属于颗粒夹杂，节点处应力的处理方法并不统一。因此，本文直接采用积分点出的应力和应变结果代入式(4)～(5)进行计算，即

(6)

(7)

2 推进剂“脱湿”损伤机理

2.1 单轴拉伸下的推进剂细观“脱湿”损伤分析

从单轴拉伸试验出发，着重分析考虑和不考虑损伤状态下，推进剂拉伸阶段的应力应变曲线。“脱湿”现象和应变速率具有较大的关联性，因为基体和黏接界面的力学性能对应变速率较为敏感。本文暂时不考虑基体和黏接界面的率相关力学性能。在细观推进剂结构拉伸过程中，结构整体的应变速率为0.2%/s。

图3(a)首先给出了结构在定速拉伸阶段，纵向应力应变曲线随时间的关系。由图可以发现，颗粒和基体黏接界面的损伤效应会明显地降低推进剂的刚度和承载能力。图3(b)为拉伸过程中纵向“脱湿”因子的变化曲线。“脱湿”因子具体的定义为

(8)

其中，D代表“脱湿”因子，σwithout和σwith分别表示不考虑损伤和考虑损伤状态下的应力。

在加载开始阶段，结构的“脱湿”因子在2 s的时间内就增加到了0.11左右。这说明，在加载的开始阶段，界面的损伤效应就已经存在。值得注意的是，在“脱湿”因子随加载时间稳定增长之前，曲线存在一个“凹陷”。本文认为“凹陷”的存在是局部黏接界面处界面刚度符号的改变导致的：在加载的初始阶段，黏接界面处的分离位移都比较小，界面刚度基本处于上升阶段。随着时间的推移，局部黏接界面处的分离位移加速增长，界面刚度开始出现下降，这一局部刚度的突然改变导致了“脱湿”因子曲线上“凹陷”的存在。局部黏接界面处的损伤效应向外传递时，大部分界面处的刚度还处于增长阶段，局部的负刚度效应会被整体的刚度抵消，因此，在很短的时间内，“脱湿”因子曲线继续随时间不断增长。

(a) 应力应变随时间分布(a) Stress and strain versus time

(b) “脱湿”因子随时间分布(b) "Dewetting" factor versus time图3 定速拉伸阶段纵向应力应变和“脱湿”因子变化曲线Fig.3 Longitudinal stress/strain and "dewetting" factor versus time in stretch period with constant speed

图4分别给出了四种典型时刻推进剂细观结构的变形。在t=20 s时刻，整个结构并没有发生相对于结构尺寸而言较为明显的裂纹。结合“脱湿”因子曲线，可以判断这一时刻，结构内部细微的裂纹早就已经让结构的承载能力产生了较大幅度的削减。在t=40 s时刻，结构内部已经开始出现明显的裂纹。如图4(b)所示，在图中的实线框，裂纹呈现“成对”出现的现象，排列顺序和拉伸方向一致。也就是说，某个界面处的裂纹并不是“孤立”存在的，在裂纹的附近，沿着结构受力的方向，总能找到另一个明显的裂纹。图4(c)中，t=60 s时刻，裂纹“成对”出现的现象更加明显。这样的“裂纹对”可以在两个颗粒与基体之间的黏接界面发生，也可以在三个颗粒与基体之间的黏接界面发生。此外，这样的“裂纹对”的产生还需要界面与界面之间有足够近的距离，图中粗虚线框的部分，由于在界面附近没有距离足够近的另一个界面，“裂纹对”现象并没有出现。图4(d)中，t=70 s时刻，“裂纹对”已经布满整个结构表面了。部分“裂纹对”的参与界面甚至从两个发展到了三个。此时，黏接界面处的裂纹已经有扩展到基体中的趋势。

(a) t=20 s (b) t=40 s

(c) t=60 s (d) t=70 s图4 定速纵向阶段推进剂变形Fig.4 Diagram of deformed propellant in stretch period with constant speed

2.2 纯剪试验下的推进剂细观“脱湿”损伤分析

本小节将进一步利用纯剪试验对推进剂细观结构的“脱湿”机理进行研究。图5分别给出了两种剪切方向(模式A和模式B)下的变形示意，切应变的应变速率均为0.02%/s。可以明显地看到，“裂纹对”在沿着剪切的方向(双箭头线)上有序分布。在垂直于剪切方向的位置上，几乎看不到明显的裂纹。此外，部分界面由于在沿着剪切方向上的黏接界面数目较少，距离较远，同样没有“裂纹对”的出现。

(a) 模式A(a) Model A (b) 模式B(b) Model B图5 t=40 s时刻模式A和B的变形Fig.5 Diagram of deformed model A and B at t=40 s

图6(a)所示为剪切阶段切应力和应变随时间的变化。以沿着二四象限的剪切模式为例，在考虑和不考虑界面损伤状态下，结构在xy以及yx方向上的应力应变曲线是一致的。这就从细观层面上验证了切应力互等原理，进一步说明了细观结构和仿真手段的准确性。另外，考虑界面损伤效应的应力应变曲线明显低于无损伤状态下的应力应变曲线。这说明界面结构的损伤会削弱推进剂的抗剪能力。图6(b)给出了两种加载方向下的“脱湿”因子变化曲线，同一模式(模式A或模式B)下不同方向上的“脱湿”因子变化曲线几乎是重合的。尽管不同模式下(模式A和模式B)的“脱湿”因子变化曲线不一致，但是误差保持在1%左右。此外，类似于单轴拉伸，“凹陷”现象同样存在，这里不再重复解释。

(a) 应力应变随时间分布(a) Stress and strain versus time

(b) “脱湿”因子随时间分布(b) "Dewetting" factor versus time图6 纯剪作用下推进剂的应力应变和“脱湿”因子变化曲线Fig.6 Stress/strain and "dewetting" factor versus time under pure shear effect

3 考虑“脱湿”损伤的推进剂力学行为分析

3.1 推进剂体分比对“脱湿”损伤的影响规律

推进剂中颗粒的体积与推进剂总体积的比值称为体分比。按照同一级配(粒径分布)生成四种不同体分比的推进剂细观有限元模型，对应的体分比Vf分别为60%、65.3%、70%以及75%。采用同样的单轴拉伸试验，纵向应变速率固定为0.2%/s。图7分别给出了这些有限元模型在考虑和不考虑界面“脱湿”损伤效应下的应力应变曲线。可以发现，随着体分比的增加，结构中固体填料增多，整体刚度变大，应力应变曲线整体水平获得了提升。但是，随着体分比的增大，考虑界面“脱湿”损伤效应下的应力应变水平比不考虑损伤下的应力应变曲线之间差异愈发明显。

(a) Vf=60%

(b) Vf=65.3%

(c) Vf=70%

(d) Vf=75%图7 不同体分比下的应力应变随时间分布Fig.7 Stress/strain versus time under different volume factors

图8 不同体分比下的推进剂“脱湿”因子随时间变化曲线Fig.8 "Dewetting" factor versus time under different volume factors

图8分别给出了这些不同体分比下的推进剂细观结构“脱湿”因子随时间的变化关系。可以看到，“凹陷”现象在初始阶段已经变成异常的上下波动，但是这样的波动之后，“脱湿”因子曲线会继续平稳发展。不同体分比下的“脱湿”因子曲线进入稳定发展阶段的时刻并没有特别明显的差异。此外，尽管体分比越高，“脱湿”因子曲线提升现象越明显，但是曲线的斜率之间并没有较大的差异。

3.2 应变速率对“脱湿”损伤的影响规律

推进剂在不同的应变速率下往往会表现出不同的力学性能。图9分别给出了体分比为65.3%、内聚强度为0.2 MPa时，不同应变速率v下的推进剂细观结构应力应变曲线。可以看到，随着应变速率的提升，应力应变曲线不断提升，相应的“脱湿”损伤效应愈明显。

图10给出了不同应变速率下的“脱湿”损伤因子变化曲线。这里，高应变速率下的“脱湿”因子曲线明显大于低应变速率下的“脱湿”因子曲

(a) v=0.1%/s

(b) v=0.2%/s

(c) v=0.3%/s

(d) v=0.4%/s图9 不同应变速率下的应力应变曲线Fig.9 Stress/strain versus time under different tensile rates

线。此外，高应变速率下的“脱湿”因子曲线斜率明显大于低应变率下的“脱湿”因子曲线。这是因为，应变速率提升的状态下，结构中基体和颗粒间的“脱湿”来不及“愈合”，那么相应的损伤效应就会加剧。

图10 应变速率对推进剂“脱湿”因子的影响Fig.10 Effects of tensile rates on "dewetting" factor

3.3 内聚强度对“脱湿”损伤的影响规律

(a) σmax=τmax=0.1 MPa

(b) σmax=τmax=0.2 MPa

(c) σmax=τmax=0.3 MPa

(d) σmax=τmax=0.4 MPa图11 不同内聚强度下应力应变随时间分布Fig.11 Stress/strain versus time under different cohesive strength

为了探究不同的内聚强度对应力应变曲线的影响，图11给出了不同内聚强度下推进剂细观结构在0.2%/s的应变速率下的应力应变曲线。在考虑界面“脱湿”损伤状态时，内聚强度越小，细观结构的“脱湿”损伤效应越明显。如果界面的内聚强度接近无穷大，结构的计算结果应该和不考虑“脱湿”损伤的状态一致。这一经验理论和仿真预示是吻合的。

图12给出了拉伸阶段，不同内聚强度下的“脱湿”因子变化曲线。可以看到，当内聚强度达到0.4 MPa时，“脱湿”因子曲线在“凹陷”点之后一直呈现下降的趋势。在PPR内聚关系中，内聚强度的提升会使得内聚刚度增大，因此在内聚强度很高时，随着分离位移的提高，结构整体的刚度加大，界面的“脱湿”损伤效应就会得到抑制。因此，在内聚强度为0.3 MPa和0.4 MPa时，结构整体的“脱湿”因子曲线会出现持续下降的趋势。也就是说，在考虑界面“脱湿”损伤状态时，应变速率越大、内聚强度越小，细观结构的“脱湿”损伤效应越明显。

图12 内聚强度对推进剂“脱湿”因子的影响Fig.12 Effects of cohesive strength on “dewetting” factor

4 结论

从推进剂细观分析模型出发，基于分子动力学的推进剂颗粒装填方法、周期几何以及周期边界的处理手段，在单轴拉伸和纯剪试验下，进行了基于内聚力模型的推进剂“脱湿”损伤机理研究，开展了推进剂体分比、应变速率、内聚强度参数对推进剂细观结构仿真结果的影响分析。主要结论如下：

1)在外载荷作用下，沿着细观结构的受力方向，“裂纹对”现象的出现是推进剂发生“脱湿”损伤的内因。

2)随着体分比的增大，推进剂细观结构“裂纹对”现象增多，“脱湿”因子曲线提升现象明显。体分比由60%提升至75%的过程中，“脱湿”因子增加近一倍。

3)应变速率的减小和内聚强度的提升会使得“脱湿”因子减小，相应地削弱推进剂细观结构的“脱湿”损伤效应。在内聚强度为0.3 MPa和0.4 MPa时，结构整体的“脱湿”因子曲线会出现持续下降的趋势。