固体发动机推进剂/绝热层界面I型脱粘力学行为试验与仿真研究

马晓琳，申志彬，崔辉如

(国防科技大学 空天科学学院 应用力学系，长沙 410073)

0 引言

作为航天运载系统的动力装置，固体发动机具有结构简单、工作可靠、性能稳定等特点，在导弹与航天技术等领域得到了广泛应用。贴壁浇注的固体发动机多采用壳体/绝热层/衬层/推进剂的多界面结构形式[1]，包括壳体/绝热层、绝热层/衬层、衬层/推进剂三类粘接界面。据统计，国外失败的固体发动机中有近1/3是因界面脱粘造成的[2]。因此，粘接界面是固体发动机结构中的薄弱环节。其中，衬层/推进剂界面是最为薄弱的部分，一旦出现脱粘现象，极易发生蹿火甚至爆炸等事故。粘接界面的力学行为直接关乎固体发动机的结构完整性，进而影响到其贮存寿命。因此，对界面脱粘的力学行为进行研究与评价具有重要的工程意义。

为了对界面粘接性能进行合理的预测及评判，有必要寻找合理的模型对粘接界面的力学行为进行描述，研究者通常选取传统断裂力学模型。然而，在实际应用中，断裂力学模型不适用于完好粘接界面的开裂情况以及开裂较大的非奇异裂纹场。此外，粘接结构实际上包含两个界面层和一个粘接层，不符合断裂力学中对裂纹的定义。因此，基于损伤力学的内聚力模型由此出现[3-4]，内聚力模型用内聚力模拟裂纹尖端的非线性断裂区，克服了线弹性断裂力学中产生的应力奇异缺陷，能较好地模拟界面的开裂过程。

本文通过商业有限元软件ABAQUS用户子程序UEL对PPR内聚力单元进行开发，设计了固体发动机推进剂/绝热层Ⅰ型界面脱粘试验，并通过基于试验的反演分析方法，获得PPR内聚力模型参数，对不同加载速率下粘接界面的力学行为进行了相关研究。

1 内聚力模型

内聚力模型从损伤力学角度出发，假设牵引力与界面开裂位移满足连续性关系，随着界面的开裂，牵引力逐渐增加至最大值，开始产生损伤，随后牵引力逐渐减小至零，损伤演变，直到位移达到临界值，此时界面已完全失效。牵引力-分离位移关系简洁明了，便于有限元实现和开发，已成为界面数值仿真的主流选择。本章将对商业有限元软件ABAQUS内置的双线性内聚力模型以及PPR内聚力模型进行介绍。

1.1 双线性内聚力模型

典型双线性内聚力模型如图1所示，符合式(1)所示的牵引力-分离位移关系：

(1)

商业有限元软件ABAQUS内置多种损伤演化类型，基于有效位移的损伤演化定义如下：

(2)

图1 典型双线性内聚力模型的牵引力-分离位移关系

商业有限元软件ABAQUS内置有4类损伤起始准则，包括最大名义应力准则、最大名义应变准则、二次名义应力准则以及二次名义应变准则。以二次名义应变准则为例，当包含名义应变的二次函数值达到1时，界面开始发生损伤。

(4)

1.2 PPR内聚力模型

PPR内聚力模型由Park、Paulino、Roesler三人[5]共同提出，其牵引力-分离位移关系可由基于势函数的内聚力模型获得，如图2所示。

内聚强度的势函数表示如下：

(5)

式中Γn、Γt、α、β、m和n分别为PPR内聚力模型的特征参数；φn和φt分别表示法向和切向的内聚能。

(6)

(7)

图2 PPR内聚力模型的牵引力-分离位移关系

基于势函数的特性，PPR势函数分别对法向和切向位移求导，即可得到法向和切向的内聚强度。

(8)

(9)

Γn和Γt表示与断裂能相关的能量常数，当界面沿法向和切向的断裂能不同时，能量常数表示为

(10)

若界面法向和切向的断裂能相同，能量常数可简化为

(11)

指数m和n表示PPR内聚力模型的初始斜率。

(12)

其中，α和β为形状参数，用于表示不同材料所产生的软化响应。当α，β≈2时，PPR内聚力模型的软化曲线为近似线性；当1<α，β<2时，PPR内聚力模型的软化曲线为凹形；当α>2，β>2时，PPR内聚力模型的软化曲线为凸形。λn和λt为初始刚度的指示因子，定义为损伤起始位移和最终破坏位移的比值。

(13)

因此，界面的最终破坏位移可表示如下：

(14)

(15)

1.3 比较和讨论

与PPR内聚力模型相比，双线性内聚力模型具有以下两点局限性：

(1)双线性模型在软化阶段提供正刚度。换言之，内聚力随分离位移的增大而增大，这类牵引力-分离位移关系通常与实际情况不符(除非材料随分离位移的增大表现出强化行为)。

基于以上讨论，选取PPR内聚力模型描述固体发动机推进剂/绝热层粘接界面的力学行为。

(a)法向内聚强度随切向位移增加而变化 (b)切向内聚强度随法向位移增加而变化

(a)法向内聚强度随切向位移增加而变化 (b)切向内聚强度随法向位移增加而变化

2 PPR内聚力单元的二次开发

商业有限元软件ABAQUS可提供用户材料(UMAT)以及用户子程序单元(UEL)两种内聚力模型的开发方式。UMAT将单元的真实分离位移转化为名义应变，基于牵引力-分离位移关系建立内聚力单元的本构模型，商业有限元软件ABAQUS内置的双线性内聚力模型采用的就是这种基于位移的方法。而UEL采用基于力的方法，将用户单元所受的牵引力转化为等效节点力。研究表明，UEL能够准确预测界面的开裂和大裂纹的扩展现象[6]，UMAT能够连续预测裂纹的产生，但不能准确预测裂纹的扩展。综上所述，本文选择UEL对PPR内聚力单元进行二次开发。

2.1 基本方程

与标准平面单元相比，粘接界面内聚力单元的有限元构成方式有所不同：一方面，内聚力单元的初始厚度为零；另一方面，在力的平衡关系中，只考虑单元上表面和下表面所受的牵引力。本文以初始厚度为零的线性二维内聚力单元为例进行讨论，如图5所示。基于虚功原理，域(Ω)内虚应变能和界面的内聚断裂能(Γc)的和与外力(Γext)在边界(Γe)上做的虚功等效，控制方程的弱形式为

(16)

式中ε、u、Δ分别为柯西应变、界面开裂位移以及单元的位移；Tc为沿断裂面分布的内聚力，其值取决于开裂位移的大小。

图5 局部坐标系下的二维内聚力单元

2.2 用户子程序实现

商业有限元软件ABAQUS中自带用户子程序接口，便于特殊有限元的实现。在一般分析步骤中，用户单元对模型的主要贡献是提供节点残差量，残差的定义为

(17)

(18)

(19)

其中，R由坐标转换矩阵组成。

(20)

(21)

其中，θ为总体坐标系和局部坐标系之间的夹角，如图6所示。

图6 局部坐标系下的二维内聚力单元及其位移

sinθ和cosθ的具体定义如下：

(22)

其中

(23)

(24)

(25)

(26)

图7 内聚力单元沿法向和切向的分离位移

式(26)可用矩阵形式表示为

(27)

其中，L为对应转换矩阵。

(28)

局部坐标系下沿内聚力单元表面的分离位移可由节点局部位移插值得到。

(29)

其中，N为形函数。

(30)

式中N1、N2为局部坐标ξ的线性函数。

(31)

(32)

其中，转换矩阵为B=NLR。

(33)

式中T0为粘接界面的真实厚度。

此外，用户子程序还需要提供单元刚度矩阵，其定义为

(34)

3 粘接界面试验及仿真分析

3.1 Ⅰ型界面脱粘试验

根据ASTM关于金属材料I型界面开裂试验规范[7]中所介绍的双悬臂梁(DCB)试验，考虑到固体推进剂和绝热层材料的刚度，参照Zhou等[8]对DCB试验的改进设计，本文设计双悬臂夹层梁(DSCB)试验件。如图8所示，铝制夹具(110 mm×5 mm)与推进剂片材(80 mm×10 mm)以及绝热层片材(80 mm×2.5 mm)通过高强度环氧树脂胶粘接，推进剂和绝热层之间为衬层丁羟胶(70 mm×0.2 mm)。为保证试验件由衬层处脱粘，在夹具左端留有10 mm空隙不浇注丁羟胶，形成预制脱粘区。整个模型的面外厚度为10 mm，左侧通过销钉与试验机连接，如图9所示。

图8 推进剂/绝热层粘接界面DSCB试验件的几何构成

(a)脱粘前 (b)脱粘后

采用微机控制的电子万能试验机进行推进剂/绝热层Ⅰ型界面脱粘试验，试验机记录下拉伸速率分别为10、50、100 mm/min时加载点的力-位移曲线，如图10所示。

图10 不同速率下DSCB试验获取的加载力-位移曲线

由试验获取的加载力-位移曲线可观察到，加载力随位移的变化可分为两个阶段：

(1)强化阶段。加载力随位移的增加逐渐增大，曲线线性上升，随后斜率稍有衰减，曲线到达峰值，粘接界面开始产生损伤。

(2)损伤演化阶段。加载力随位移的增加逐渐减小，裂纹逐渐扩展，直至完全破坏。此外，随着加载速率的增大，粘接界面的内聚能和强度均明显增大，表明粘接界面的力学行为具有显著的速率相关性。但在强化阶段，不同加载速率下的加载力-位移曲线刚度相差不大，表明粘接界面力学行为的速率相关性主要体现在损伤演化阶段。

3.2 商业有限元软件ABAQUS仿真

为了简化建模，忽略铝制夹具和推进剂以及绝热层之间的环氧树脂胶，不考虑二者间的界面。如图11所示，运用商业有限元软件ABAQUS对DSCB试验件建模仿真，铝制夹具、推进剂和绝热层均采用四节点四边形线性缩减积分平面应力单元(CPS4R)，衬层采用四节点二维内聚力单元(COH2D4)，其初始厚度设置为0.2 mm。约束试验件左上端节点沿X、Y方向的位移，以及左下端节点沿X方向的位移，对左下端节点以恒定速度施加位移。

铝制夹具的杨氏模量为70 000 MPa，泊松比为0.33。推进剂松弛模量用Prony级数表示为

(35)

其中，E0为初始模量；n为Prony级数的阶数；t为时间；Ei为τi时间下Prony级数的各阶系数，表 1为5阶Prony级数拟合所得系数。此外，推进剂泊松比假设为0.498。

图11 DCSB试验的有限元模型及其边界条件

n012345τi/s—0.0030.0701.88450.911376Ei/MPa55.733.18.724.472.692.51

本文采用基于应变能势函数的Mooney-Rivlin模型描述绝热层的超弹性力学行为。

W=C10(I1-3)+C01(I2-3)

(36)

其中，I1和I2分别为左Cauthy-Green应变张量的第一、第二不变量；C10=-1.162 1和C01=1.672 3为Mooney-Rivlin模型的系数，可通过最小二乘法对试验曲线进行回归得到。

对于衬层材料，采用UEL实现PPR内聚力模型φn、φt、Tn、Tt、α、β、m和n8个参数的设置。仿真结果如图12所示，内聚力单元达到破坏条件时即可删除。

图12 DCSB试验的有限元仿真结果(应力云图)

3.3 内聚力模型参数的获取

根据试验获取的加载力-位移曲线为宏观结果，仅能对界面内聚能和内聚强度的值进行粗略预估，然而，对于粘接界面问题，界面的起裂以及损伤演化特性才是主要关注点。因此，有必要基于试验结果对仿真模型的内聚力参数进行相应修正，以获取更为准确的内聚力模型参数，即为基于试验的反演分析方法。构建仿真曲线与试验曲线之间的误差，作为反演优化的目标函数 。

(37)

考虑到目标函数不存在梯度，且参数空间非连续，故本文采用Hooke-Jeeves算法对PPR内聚力模型的8个参数进行直接搜索优化，反演优化获取的加载力-位移曲线同试验曲线对比如图13所示。从图13可看出，仿真曲线和试验曲线重合度较好，说明PPR内聚力模型能够较好的描述界面的脱粘过程。

表2展示了不同加载速率下反演优化得到的PPR内聚力模型参数，根据已获得的PPR内聚力模型参数，可得到不同加载速率下的内聚力模型的法向牵引力-分离位移关系，如图14所示。结合表2分析可知，随着加载速率的增大，PPR内聚力模型对应的法向内聚能、内聚强度均增大，法向初始刚度、损伤起始位移均减小，具有明显的速率相关性。此外，10 mm/min时的最终开裂位移介于50 mm/min和100 mm/min之间，这是由曲线形状不同导致的，50 mm/min对应的α=2.035，近似于2，PPR内聚力模型的软化曲线为近似线性，而10 mm/min和100 mm/min对应的α均大于2，PPR内聚力模型的软化曲线为凹形。

(a)10 mm/min (b)50 mm/min (c)100 mm/min

图14 不同加载速率下反演优化获取的PPR内聚力模型

PPR内聚力模型参数加载速率/(mm/min)1050100ϕn/(N/mm)0.761 11.3883.428ϕt/(N/mm)19.015.8952.545Tn/MPa4.22412.5719.53Tt/MPa1.79214.3133.83α9.9982.0356.050β9.4628.6331.047λn0.170 10.093 760.002 676λt0.151 00.003 5140.483 0

4 结论

(1)推进剂/绝热层Ⅰ型界面脱粘试验表明，加载力随位移的变化可分为强化阶段和损伤演化阶段，且粘接界面的力学行为具有显著的速率相关性，主要体现在损伤演化阶段。

(2)通过商业有限元软件ABAQUS用户子程序对PPR内聚力单元进行开发，仿真结果与试验结果重合度较好，表明PPR内聚力模型能较好地描述界面脱粘过程。

(3)通过基于试验的反演分析方法获得PPR内聚力模型参数，随着加载速率的增大，法向内聚能、内聚强度均增大，法向初始刚度、损伤起始位移均减小，表明粘接界面的力学行为具有显著的速率相关性。