1 引言

复合固体推进剂不仅是固体火箭发动机的能量物质, 也是其重要的结构材料, 在成型、贮存与使用中, 往往存在各种机械与温度载荷, 引起各种不同形式的微损伤, 主要有粒子断裂、基体损伤、粒子与基体的界面脱湿^[1]。文献[2]与[3]分别从试验与理论上, 确定在复合固体推进剂中主要损伤形式是粒子与基体的界面脱湿。界面脱湿是一种细观的力学行为, 是引起复合固体推进剂本构关系从线性行为到非线性行为的主要决定因素, 并有可能发展成宏观裂纹, 从而导致复合固体推进剂失效。在实验中, 单独地分析界面脱湿对复合固体推进剂本构关系的影响, 存在诸多影响因素, 难以展开。在理论上, 唯象的宏观本构关系, 一般不考虑复合固体推进剂中的细观结构, 如粒子的尺寸, 粒子的体积分数, 及粒子与基体的界面粘结情况等等。因此, 需要发展具有细观特征的宏细观结合模型, 分析界面脱湿对宏观本构关系的影响, 并正确预测复合固体推进剂的力学性能。

复合固体推进剂是一种粒子增强体复合材料, 主要由聚合物基体与刚性的氧化剂粒子组成^[3]。国外, 针对复合固体推进剂已开展了一些研究, 如Schapery^[4]从不可逆热力学角度, 引入两个内变量描述复合固体推进剂的损伤, 提出了损伤本构关系。Vratsanos与Farris^[5]从能量守恒角度, 假定粒子的脱湿是从大粒子到小粒子依次脱湿, 认为粒子的脱湿是导致复合固体推进剂非线性的主要原因, 并确定了宏观本构关系。Ravichandran与Liu^[6]基于Eshelby等效夹杂理论与Mori-Tanaka方法, 考虑了复合固体推进剂中的界面损伤, 提出了一个率无关的唯象本构关系。Tan^{[3, 7]}针对PBX9501高能炸药, 使用数字图像相关(DIC)技术, 获取了基体与粒子的界面粘性定律, 并确定了模型参数, 并在RVE(Representative Volume Element)上, 使用平均化与Mori-Tanaka法, 确定了高三轴应力下的线弹性宏观本构关系, 利用有限元技术数值模拟了宏观本构关系, 显示了增强体粒径大小对PBX9501高能炸药的界面脱湿具有重大影响, 并得出了临界脱湿粒径。之后, Tan^[8]在线弹性本构关系的基础上, 使用Laplace变换原理, 将线弹性宏观本构关系发展到线粘弹性的宏观本构关系。国内, 彭威^[1]针对复合固体推进剂, 考虑粒子增强作用与界面脱湿损伤, 建立了含损伤变量的粘弹性宏观本构关系, 并与拉伸曲线对比, 结果吻合较好。李丹与胡更开^[9]针对高粒子体积分数聚合物材料, 基于Laplace变换和双夹杂相互作用的弹性模型, 建立了细观力学模型, 计算了玻璃微珠/ED-6复合材料的有效松弛模量及常应变率下的本构关系。

可以看出, 国内外针对复合固体推进剂, 建立宏细观结合的损伤本构关系已经展开研究, 但是对于细观结构对宏观本构关系的影响难于从定量上确定。因此, 本文针对复合固体推进剂, 对于其内部细观结构的界面脱湿, 假定含粒子与空泡的三相物理模型描述了经历界面脱湿的复合固体推进剂, 利用细观理论与能量守恒原理从定量上确定不同载荷作用下, 增强粒子的脱湿程度对宏观力学行为的影响, 并确定了宏观有效割线模量。最后, 针对粒子粒径服从对数正态分布的复合固体推进剂, 展开了数值模拟, 获取计算结果, 与实验数据对比, 结果吻合较好^[10]; 并针对高粒子体积分数的复合固体推进剂, 分析了粒子体积分数对复合固体推进剂宏观本构关系的影响。

2 能量守恒定律的表达形式

为了确定复合固体推进剂含损伤的本构关系, Vratsanos与Farris根据热力学第一定律—能量守恒定律, 假定复合固体推进剂中的损伤完全由界面脱湿决定, 从而确定能量平衡方程^[5]:

$ \delta {U_{{\rm{strain}}}} + \delta {U_{{\rm{surf}}}} = \delta W + \delta Q $

(1)

式中, δU_strain是复合固体推进剂中粒子与基体存储的应变能, δU_surf是由脱湿而引起的界面耗散能, δQ是外界导入系统的净热量, δW是外界对系统做的净功, 单位均是J; δ是微小量符号。假定系统是绝热的, 并使用虚功原理, 能量平衡方程用应力与应变如下表达^[5]:

$ 2{G_c}\frac{{\delta A}}{{{V_0}}} = {\sigma _{ij}}\delta {\varepsilon _{ij}}-\delta {\sigma _{ij}}{\varepsilon _{ij}} $

(2)

式中, G_c是基体与粒子界面的单位面积耗散能, 单位是J·m^-2; A是基体与粒子脱湿界面的面积, m²; σ_ij, ε_ij分别是复合固体推进剂的应力与应变, 下标满足张量运算; 由于粒子与基体界面脱湿是成对的发生, 因此左端存在一个系数2, V₀是复合固体推进剂的初始体积, m³。

在单轴拉伸下, 给定应力边界条件, 得到控制方程:

$ \left\{ \begin{array}{l} 2{G_c}\frac{{\delta A}}{{{V_0}}} = {\sigma _{ij}}\delta {\varepsilon _{ij}}-\delta {\sigma _{ij}}{\varepsilon _{ij}}\\ {\sigma _{11}} = {\sigma _0};{\sigma _{22}} = {\sigma _{33}} = 0; \end{array} \right. $

(3)

利用复合固体推进剂的本构关系:

$ {\varepsilon _{ij}} =-\frac{\nu }{E}{\sigma _{kk}}{\delta _{ij}} + \frac{{1 + \nu }}{E}{\sigma _{ij}} $

(4)

式中, E为宏观的弹性模量, ν为宏观的泊松比, δ_ij是克罗内克符号; 最后, 可以确定临界脱湿应变方程:

$ -2\frac{{{G_c}}}{{{V_0}}}\left( {\frac{{{\rm{d}}A}}{{{\rm{d}}c}}} \right) = {\varepsilon _{11}}^2\left( {\frac{{{\rm{d}}E}}{{{\rm{d}}c}}} \right) $

(5)

式中,d是微分符号; 因此, 只需要确定粒子与基体的界面脱湿面积A随粒子体积分数c的变化关系, 及复合固体推进剂中弹性模量E随粒子体积分数c的变化关系, 就得到临界应变, 并确定应力。利用几何关系可以确定界面脱湿面积A的变化; 利用细观力学理论Mori-Tanaka方法确定弹性模量的变化。

假定复合固体推进剂中的粒子是球形, 那么一个半径为R的粒子发生完全脱湿, 面积增加4πR², 体积分数减少$\frac{4}{3}{\rm{ \mathsf{ π} }}{R^3}/{V_0}$, 因此,

$ \frac{{{\rm{d}}A}}{{{\rm{d}}c}} =-6{V_0}/R $

(6)

3 Mori-Tanaka方法的有效场理论

粒子增强体复合材料的细观力学有效场理论, 常见方法有广义自洽法(GSM), 微分法和Mori-Tanaka方法(M-T)。Christensen^[11]针对高体积分数的球形粒子增强体复合材料, 使用这三种方法分别计算了复合材料的弹性性能, 比较结果得出, M-T方法的精度与GSM相当, 微分法次之, 但GSM方法用于多相复合材料具有局限性。另外, Weng^[12]证明对于多相复合材料, M-T方法的精度较高。

假设复合固体推进剂是增强粒子与聚合物基体组成的两相复合材料, 但是由于粒子的脱湿, 使用空泡代替粒子, 从而具有与粒子不同弹性属性的第三相, 因此复合固体推进剂属于三相复合材料。本文使用M-T方法计算包含空泡的复合材料有效属性。引用文献^[13]的推导结果, 设复合材料的平均刚度矩阵是[C], 基体的刚度矩阵是[C⁰], 增强粒子的刚度矩阵是[Cⁱ], 以及由粒子脱湿而形成的空泡刚度矩阵是[C^ν], 那么有如下关系:

$ \begin{array}{l} \left[{\bar C} \right] = \left[{{C^0}} \right] + \\ \;\;\;\;\;\;\;\;\;\;[{C^0}]*\{ {c^i}\left[{{{\mathit\Gamma }^i}} \right]({c^i}\left[{I-{S^i}-{{\mathit\Gamma}^i}} \right]+\left[{{S^i}} \right] + \left[A \right] + \\ \;\;\;\;\;\;\;\;\;\;{c^\nu }\left[{I-{S^\nu }-{\mathit\Gamma ^v}} \right]{\left[{{S^\nu } + B} \right]^{ - 1}}{\left[{{S^\nu } + B} \right]^{ - 1}}\left[{{S^i} + A} \right]{)^{ - 1}}\} + \\ \;\;\;\;\;\;\;\;\;\;[{C^0}]*\{ {c^\nu }\left[{{{\mathit\Gamma}^\nu }} \right]({c^\nu }\left[{I-{S^\nu }-{\mathit\Gamma ^\nu }} \right] + \left[{{S^\nu }} \right] + \left[B \right] + \\ \;\;\;\;\;\;\;\;\;\;{c^i}\left[{I-{S^i}-{\mathit\Gamma ^i}} \right]{\left[{{S^i} + A} \right]^{ - 1}}\left[{{S^\nu } + B} \right]{)^{ -1}}\} \end{array} $

(7)

式中, 方括号表示方阵, [A]=[Cⁱ-C⁰]^-1[C⁰], [B]=[C^ν-C⁰]^-1[C⁰], cⁱ表示粒子的体积分数, c^ν表示由粒子脱湿而产生的空泡体积分数。[S]是Eshelby矩阵, 依赖于基体的泊松比及增强粒子的形状, 对于球形粒子, 有^[12]:

$ {S_{ijkl}} = \frac{1}{{15\left( {1-{\nu _0}} \right)}}\{ \left( {5{\nu _0}-1} \right){\delta _{ij}}{\delta _{kl}} + \left( {4-5{\nu _0}} \right)({\delta _{ik}}{\delta _{jl}} + {\delta _{il}}{\delta _{jk}})\} $

(8)

在单轴拉伸下, 粒子不同部位脱湿程度不同, 一般是从极区开始脱湿, 因此, 引入因子F_b表示不同的脱湿程度, 那么空泡的刚度矩阵, 表示为^[10]:

$ {\left[{{C^\nu }} \right] = m\left[{\begin{array}{*{20}{c}} {{F_b}{E_{11}}(1-{\nu _{23}}^2)}&0&0\\ 0&{{E_{22}}}&{{E_{22}}{\nu _{32}}}\\ 0&{{E_{33}}{\nu _{23}}}&{{E_{33}}} \end{array}} \right]} $

(9)

式中, m=(1-ν₂₃ν₃₂)^-1, E_ii(1≤i≤3)是弹性粒子的弹性模量, ν_ij(1≤i, j≤3)是弹性粒子的泊松比。为了考虑相与相之间的相互影响, 公式(7)中的[Γ^r]表示相互影响矩阵, 当r=i表示粒子; 当r=ν表示空泡, 相互影响矩阵表示为^[10]:

$ \left[{{\mathit\Gamma ^r}} \right] = \left[I \right] + \frac{{5{c^r}}}{{4{\beta ^2}}}Y\left[{{W^r}} \right] $

(10)

式中, [I]表示单位方阵, c^r表示相r的体积分数; Y=Y_m(1-c^r), Y_m为相互影响系数; [W^r]=ξ₁δ_ijδ_kl+ξ₂(δ_ikδ_jl+δ_ilδ_jk), 参数ξ₁, ξ₂, β的取值依据文献[13]。

4 数值模拟与计算分析

为了验证本文算法的正确性, 针对复合固体推进剂HTPB展开数值模拟。颗粒粒径频度服从对数正态分布, 采用文献[10]的推进剂试件数据, 确定了平均粒径r =32 μm, 对数标准差lnσ=0.167, 及初始粒子体积分数fⁱ=30%与初始的空泡体积分数f ^ν=0;本文预估了界面粘结能G_c=2.241 J·m^-2, 部分脱湿因子F_b=10^-5, 及相互影响系数Y_m=1.18。依据文献[10]确定了复合固体推进剂组分的物理属性, 玻璃珠增强粒子的剪切模量G_i=30 GPa与泊松比ν_i=0.16;基体的泊松比ν₀=0.495, 考虑非线性弹性基体, 弹性模量与应变相关, 单位MPa, 表示如下^[10]:

$ {E_0} = 1.554865-0.497499\varepsilon + 0.321452{\varepsilon ^2} $

(11)

采用公式(5)与(7)模拟计算, 得到计算结果如图 1, 为了比较, 同时给出了文献[10]的实验数据。可以看出, 计算结果与试验数据较为吻合。易知, 复合固体推进剂的力学行为主要分为两个阶段, 第一阶段是线弹性阶段, 粒子没有发生脱湿, 粒子对复合固体推进剂起到增强作用; 第二阶段是非线性阶段, 由于粒子的脱湿, 导致复合固体推进剂中空泡的产生, 粒子的增强作用减弱, 引起复合固体推进剂弹性模量的下降, 验证了复合固体推进剂中的非线性行为主要是界面脱湿, 与文献^[2]的结论相一致。

针对目前国内使用的复合固体推进剂HTPB, 主要由高氯酸铵(AP)、铝粒子与粘合剂基体组成, 固体含量体积分数高达90%, 其中AP粒子体积分数70%。铝粒子的粒径较AP粒子的粒径小得多, 在拉伸试件破坏之前, 铝粒子的界面基本粘结完好, 对复合固体推进剂主要起到增强作用, 因此作为简单考虑, 作者针对含AP粒子与粘合剂基体的复合推进剂展开数值模拟, 研究了不同AP粒子体积分数, 界面脱湿对复合固体推进剂力学行为的影响。

从图 2中可以看出, AP粒子体积分数高的HTPB比AP粒子体积分数低的HTPB具有更高的强度, 但是AP粒子发生脱湿后, 导致宏观本构关系的非线性行为出现, 同时应力下降的趋势更明显, 粒子的增强作用转化为软化作用。这是因为粒子体积分数含量越高, 对HTPB起到越高的增强作用, 但是高的粒子体积分数, 导致HTPB更容易发生脱湿, 粒子的增强作用减弱较快, 表现出更明显的应力下降趋势。

图 3给出了HTPB宏观有效弹性模量变化, 可以看出有效弹性模量随应变的增大而逐渐降低, 但在非线性阶段, 弹性模量下降速率加快。这是因为, 一方面, 本文使用非线性弹性基体, 基体弹性模量随应变增大而降低, 因此在粒子粘结时, 宏观有效模量随载荷的增大而降低; 另一方面, 复合固体推进剂HTPB内粒子脱湿数量随应变增大而增加, 也就是损伤值增大, 粒子的增强作用减弱, 从而宏观有效弹性模量随损伤值增大而降低。同时易知, HTPB中包含越多的粒子, 就具有更高的初始弹性模量, 在图 3中粒子含量的宏观有效模量是粒子含量的4倍, 但是粒子含量的宏观有效模量退化速率也更快。图 4给出了HTPB体积膨胀应变的变化趋势, 可以看出体积膨胀应变随应变的增大而增大, 在非线性阶段, 由于粒子脱湿而产生空泡, 导致体积膨胀应变的增长, 也就是HTPB内部损伤值的上升。由于粒子体积分数高的HTPB更容易脱湿, 导致在HTPB中空泡更快增长, 从而体积膨胀应变更大。

5 结论

针对HTPB推进剂, 建立了数学模型, 模拟了界面脱湿对其力学性能的影响, 可以推广至一般的粒子增强体复合固体推进剂。并从以上的数值模拟与分析中, 可以得出:

第一, 对于弱粘结界面的复合固体推进剂, 容易发生界面脱湿及空泡出现, 粒子的增强作用减弱, 导致宏观模量的下降, 但损伤值与体积膨胀都增大。

第二, 增加复合固体推进剂的粒子含量, 可以明显地提高宏观模量与强度, 但更容易发生界面脱湿。