1 引言

复合固体推进剂作为目前应用最为广泛的一种固体推进剂, 其初始力学性能及稳定性关系到固体发动机的结构完整性。目前针对复合固体推进剂粘弹性力学的研究, 国内外众多学者已经开发了大量的研究方法及本构模型。Schapery^[1-2]从经典的线粘弹性理论出发, 通过引入折算因子将线性粘弹性本构关系拓展至非线性范围内。Yun^[3]在此基础上考虑了Mullin′s效应及“脱湿”所形成的孔洞对推进剂力学性能的影响, 提出了各向同性含损伤的复合固体推进剂粘弹性本构方程。王哲君^[4]在研究端羟基聚丁二烯(HTPB)推进剂单轴拉伸试验时发现, 其应力-应变关系也表现出明显的超弹特性, 因此为更准确地研究HTPB推进剂在静态及准静态下的力学响应, 需构建同时考虑推进剂粘弹特性与超弹特性的粘超弹本构模型。胡少青^[5]将Yang^[6]等提出的适用于橡胶材料的粘超弹本构模型改进, 用于描述HTPB/IPDI(Isophorone diisocyanate)弹性体的静态及准静态力学响应, 并通过单轴等速拉伸和拉伸松弛试验获得了准确的模型参数。Zhang^[7]在基于Yang^[6]和Pouriayevali^[8-9]的模型基础上引入了率相关温度项, 使得修正后的热-粘超弹本构模型能准确描述不同温度下宽泛应变率范围的NEPE(Nitrate ester plasticized polyether)推进剂的力学性能。王鸿丽^[10]通过一系列的蠕变-回复试验, 推导出改性双基推进剂粘弹-粘塑性本构模型, 并拟合得到了所有模型参数。

以上针对推进剂宏观力学行为的研究都是将推进剂视为宏观均质材料, 基于连续介质力学理论推导得出相应的宏观本构方程。事实上, 复合固体推进剂作为一种非均质的多相复合材料, 在制备过程中由于工艺限制会存在一些初始缺陷, 主要是在固化降温过程中由于颗粒和基体导热性能差异使得二者变形不同而导致的填料-基体界面粘结缺陷, 而考虑填料-基体界面初始缺陷的推进剂本构方程还未见报道。细观实验研究发现初始缺陷的存在会改变推进剂内的细观结构组成, 载荷作用下缺陷的扩展及裂纹的形成很大程度上改变了推进剂的力学性能^[11]。为此, 为定量研究初始缺陷对推进剂力学性能的影响规律, 本研究开展了HTPB推进剂的单轴拉伸及多步松弛试验, 构建了静态及准静态下考虑初始缺陷的HTPB推进剂粘超弹本构模型, 拟合得到了相关的模型参数, 以期为固体火箭发动机推进剂药柱结构完整性分析研究提供理论指导。

2 实验部分 2.1 试件制备及实验方法

采用的HTPB推进剂系由辽宁庆阳化学工业公司提供, 其典型配方组成如表 1所示, 其中Al粉为燃料、高氯酸铵(AP)为氧化剂、HTPB为粘合剂、其余为各种添加剂。为获得含初始界面缺陷的HTPB推进剂, 在制备时分别对质量分数为0%、10%、20%和30%的AP颗粒进行表面处理, 即在其表面滚上一层石蜡, 过筛后在推进剂药浆中混合后进行浇注成型, 使得其对应的颗粒-基体界面处于初始脱粘状态。所有试样经同一批次加工, 参照QJ924-1985标准, 采用刀具冲压法将推进剂制成哑铃型单轴拉伸试件, 尺寸如图 1所示。

在10 mm·min^-1的速率下对不含缺陷的推进剂进行多步松弛实验, 每次将试件拉伸至2%应变, 保持应变不变, 静置600 s; 重复若干次, 直至试件断裂。每次松弛后的应力视为该应变下应力的平衡位置, 依次连接各个平衡应力点, 即可得到HTPB推进剂的平衡响应曲线, 如图 2所示。由于平衡曲线是经过若干次松弛得到的, 因此可将其看成是HTPB推进剂应力应变曲线的超弹性部分^[12]。

分别采用1，5，20，100 mm·min^-1的速率, 对HTPB推进剂进行单轴等速拉伸。取5个试件中一致性较好的3个实验结果的平均值为该速率下的HTPB推进剂单轴拉伸应力-应变曲线, 实验温度为291 K, 相对湿度为48%。将单轴拉伸曲线与平衡响应曲线进行对照, 如图 3所示。

2.2 实验结果分析

由图 3可以看到, 不同拉伸速率下的HTPB推进剂拉伸应力-应变曲线与平衡响应曲线趋势一致, 都是由弹性上升段、应变软化段及应力平台段组成。HTPB推进剂力学行为具有明显的应变率相关性, 即随着加载速率的增加, 材料的初始模量及应力幅值增大, 同一应变下的应力值也增大。因此可以认为HTPB推进剂是一种粘超弹材料。根据周东谟^[13]等人的研究, 可认为在小应变的范围内, 推进剂内部基本无损伤, 因此该阶段应力随应变增加呈现线性增长; 随着应变增加, 推进剂的应力增加放缓, 这主要是因为HTPB粘合剂高分子的键长、键角及小的运动单元的变化引起的。应变进一步增加, 应力逐渐出现平台区, 该阶段内颗粒脱湿导致应变快速增大, 填料-基体界面出现粘结断裂。此时界面已基本丧失承载能力, 随着应变的增加应力几乎保持不变, 这也正是导致应力出现平台的内在原因。

3 粘超弹本构模型

由于各个应变率下的HTPB推进剂拉伸应力-应变曲线均与平衡响应曲线趋势一致, 因此本研究在超弹性本构模型的基础上, 引入应变率参数来描述其率相关特性, 从而建立HTPB推进剂在不同应变率下的粘超弹本构模型, 其表达式可表示为^[14]:

$ \sigma \left( \varepsilon \right) = M\left( {\varepsilon ,\dot \varepsilon } \right)f\left( \varepsilon \right) $

(1)

式中, σ为工程应力; ε为工程应变; ${\dot \varepsilon }$ 为应变率; ${\dot \varepsilon }$ 为应变率相关项; f(ε)为HTPB推进剂超弹部分的本构方程。

3.1 超弹部分

超弹性本构理论认为材料应力应变关系是完全弹性的, 即外力消失后变形可以恢复到原始状态。

HTPB推进剂具体超弹性的力学特性, 其特征是存在某一潜在的能量(应变)函数, 可使用应变能函数W(I₁, I₂, I₃)来表示其本构模型。而后通过第二Piola-Kirchhoff应力张量σ_ij和Green-Lagrange应变张量ε_ij来表示材料的应力-应变关系^[5]:

$ {\sigma _{ij}} = \frac{{\partial W\partial {I_1}}}{{\partial {I_1}\partial {\varepsilon _{ij}}}} + \frac{{\partial W\partial {I_2}}}{{\partial {I_2}\partial {\varepsilon _{ij}}}} + \frac{{\partial W\partial {I_3}}}{{\partial {I_3}\partial {\varepsilon _{ij}}}} $

(2)

式中, I₁、I₂和I₃为左Cauchy变形张量B的应变不变量。材料具有不可压特性时, I₃=1。根据式(2)及Green应变不变量推导得材料应力应变公式:

$ \sigma = 2\left( {\lambda - {\lambda ^{ - 2}}} \right)(\frac{{\partial W}}{{\partial {I_1}}} + {\lambda ^{ - 1}}\frac{{\partial W}}{{\partial {I_2}}}) $

(3)

式中, I₁和I₂为左Cauchy变形张量B的第一、二不变量; λ为主伸长量, 且λ=1+ε, ε是加载方向上的应变。

3.1.1 Mooney-Rivlin模型

Mooney^[15-16]通过物质相变理论及大量试验, 探讨了不可压缩各向同性超弹性材料的有限变形弹性理论, 该模型基于右Cauchy-Green变形张量的第一和第二不变量的应变能函数, 其表达式为^[15]:

$ W = \sum\limits_{i,j = 0}^N {{C_{ij}}{{\left( {{I_1} - 3} \right)}^i}{{\left( {{I_2} - 3} \right)}^j}} $

(4)

Tschoegl^[17-18]认为, 含高阶项的Mooney-Rivlin模型能更好描述填充与非填充橡胶材料, 因此本研究采用五项Mooney-Rivlin模型, 其应变能函数为^[17]:

$ \begin{array}{l} W = {C_{10}}\left( {{I_1} - 3} \right) + {C_{01}}\left( {{I_2} - 3} \right) + {C_{20}}{\left( {{I_1} - 3} \right)^2} + \\ \quad {C_{02}}{\left( {{I_2} - 3} \right)^2} + {C_{11}}\left( {{I_1} - 3} \right)\left( {{I_2} - 3} \right) \end{array} $

(5)

式中, C₁₀、C₀₁、C₂₀、C₀₂、C₁₁为材料超弹参数。

将式(5)代入式(2), 得到五项Mooney-Rivlin模型名义应力-应变表达式为^[14]:

$ \begin{array}{l} \sigma = 2\left[ {\left( {1 + \varepsilon } \right) - {{\left( {1 + \varepsilon } \right)}^{ - 2}}} \right]\left\langle {{C_{10}} + 2{C_{20}}[{{\left( {1 + \varepsilon } \right)}^2} + } \right.\\ \quad 2{\left( {1 + \varepsilon } \right)^{ - 1}} - 3] + {C_{11}}\left[ {2\left( {1 + \varepsilon } \right) + {{\left( {1 + \varepsilon } \right)}^{ - 2}} - 3} \right] + \\ \quad {\left( {1 + \varepsilon } \right)^{ - 1}}\{ {C_{01}} + 2{C_{02}}\left[ {2\left( {1 + \varepsilon } \right) + {{\left( {1 + \varepsilon } \right)}^{ - 2}} - 3} \right] + \\ \quad \left. {{C_{11}}\left[ {{{\left( {1 + \varepsilon } \right)}^2} + 2{{\left( {1 + \varepsilon } \right)}^{ - 1}} - 3} \right]\} } \right\rangle \end{array} $

(6)

该模型形式简单、便于实现, 广泛应用于工程研究之中, 可较好地拟合中等应变范围材料的力学实验结果。

3.1.2 Ogden模型

Ogden^[19-20]为了克服采用变形张量不变量引起的关系式复杂, 提出了以主伸长率表征的应变能函数, 其形式为^[19]:

$ W = \sum\limits_{n = 1}^\infty {\frac{{2{\mu _n}}}{{\alpha _n^2}}\left( {\lambda _1^{{\alpha _n}} + \lambda _2^{{\alpha _n}} + \lambda _3^{{\alpha _n}} - 3} \right)} $

(7)

式中, μ_n和α_n是材料参数, 可通过非线性最小二乘法拟合得到。级数的项数可根据实验结果进行选择, 所以Ogden模型具有很大的灵活性。将式(7)代入式(2)可得名义应力-应变关系为^[20]:

$ \sigma = \sum\limits_{n = 1}^\infty {\frac{{2{\mu _n}}}{{{\alpha _n}}}\left[ {{{\left( {1 + \varepsilon } \right)}^{{\alpha _n} - 1}} - {{\left( {1 + \varepsilon } \right)}^{ - {\alpha _n}/2 - 1}}} \right]} $

(8)

3.1.3 参数拟合与分析

分别采用五项Mooney-Rivlin模型、Ogden(n=2)模型和Ogden(n=3)模型对HTPB推进剂超弹性的拉伸应力-应变曲线进行拟合, 拟合结果见图 4和表 2所示。

由表 2可以看出, 3种超弹性本构模型对实验结果的拟合都比较好, 考虑到Ogden模型本身的稳定性好, 并且选用伸长率表示的本构模型比选用应变不变量表示的本构模型更为直观, 加之2阶的Ogden模型参数少, 形式简便, 因此本研究选用Ogden(n=2)模型作为HTPB推进剂超弹性部分f(ε)的本构方程。

3.2 粘弹部分

根据式(1)中提出的粘超弹本构模型, 应变率参数的求取采用各应变率下, 任一应变时应力与该应变下应力的超弹性部分相除得到, 即:

$ M\left( {\varepsilon ,\dot \varepsilon } \right) = \frac{{{\sigma _{{\rm{test}}}}\left( {\varepsilon ,\dot \varepsilon } \right)}}{{{f_{{\rm{test}}}}\left( \varepsilon \right)}} $

(9)

式中, ${{\sigma _{{\rm{test}}}}\left({\varepsilon, \dot \varepsilon } \right)}$ 为HTPB推进剂单轴拉伸实验中各应变率下的应力值; f_test(ε)为多步松弛所得的推进剂的平衡应力。

计算结果列于图 5, 明显看出各个应变率下, 应变率参数 $M\left({\varepsilon, \dot \varepsilon } \right)$ 与应变表现出一种线性关系, 因此对不同应变率下的曲线进行线性拟合, 拟合结果如表 3所示。

将表 3中各参数随应变率 ${\dot \varepsilon }$ 的变化规律画于图 6中, 不难发现, 参数b和k与应变率成对数关系, 因此选用对数函数对 $b -\dot \varepsilon $ 和 $k -\dot \varepsilon $ 进行拟合, 即令:

$ \begin{array}{l} b = {p_1} + {q_1}{\rm{ln}}\left( {\dot \varepsilon } \right)\\ k = {p_2} + {q_2}{\rm{ln}}\left( {\dot \varepsilon } \right) \end{array} $

(10)

拟合结果如图 6和表 4所示。可以看出, 采用对数函数对参数b和k拟合的结果较好, 因此参数b和k可表示为应变率对数的函数。

3.3 粘超弹本构模型

将式(8)~式(10)代入方程(1), 即可得到HTPB推进剂粘超弹本构方程为:

$ \begin{array}{l} \sigma = \left[ {({p_2} + {q_2}{\rm{ln}}\left( {\dot \varepsilon } \right)) \times \varepsilon + ({p_1} + {q_1}{\rm{ln}}\left( {\dot \varepsilon } \right))} \right] \times \\ \quad \sum\limits_{n = 1}^2 {\frac{{2{\mu _n}}}{{{\alpha _n}}}\left[ {{{\left( {1 + \varepsilon } \right)}^{{\alpha _n} - 1}} - {{\left( {1 + \varepsilon } \right)}^{\left( { - {\alpha _n}/2} \right) - 1}}} \right]} \end{array} $

(11)

将表 2和表 4中的模型参数代入式(11)即可得到只关于应变和应变率的本构方程, 该方程可用于描述HTPB推进剂在有限范围内任一应变率下的应力应变响应。将单轴拉伸实验的应变率代入式(11), 将计算得到的曲线与实验结果进行对比, 结果见图 7。可以看出, 该模型对实验结果的拟合较好, 可用于描述HTPB推进剂在准静态单轴拉伸下的力学性能。

为验证模型的适用性, 补做了拉伸速率为10 mm·min^-1的验证试验, 并通过式(11)计算得到该拉伸速率下的预测应力-应变曲线, 结果见图 7。从图 7中可以看出, 预测曲线和实验结果吻合较好, 同时在获取应变率参数 $M\left({\varepsilon, \dot \varepsilon } \right)$ 时并没有用到该应变率下的实验数据, 这说明构建的粘超弹本构模型能用于预测HTPB推进剂应变率介于3.33×10^-4 s^-1和3.33×10^-2 s^-1之间时的应力应变响应。

4 含缺陷粘超弹本构模型

参照文献[21]中对含孔洞缺陷的铝合金本构模型的构建方法, 在前节模型基础上, 引入缺陷损伤因子, 构建了含缺陷HTPB推进剂粘超弹本构模型:

$ \sigma \left( \varepsilon \right) = \left[ {1 - f\left( {\eta ,\varepsilon } \right)} \right] \cdot {\sigma _0}\left( {\varepsilon ,\dot \varepsilon } \right) $

(12)

式中, ${\sigma _0}\left({\varepsilon, \dot \varepsilon } \right)$ 为不含缺陷的HTPB推进剂粘超弹本构方程; f(η, ε)为初始缺陷损伤因子; η为初始缺陷含量。

根据式(12)提出的本构方程, 初始缺陷损伤因子f(η, ε)的求取采用不同缺陷含量下应力与不含缺陷时的应力之比得到, 即:

$ f\left( {\eta ,\varepsilon } \right) = 1 - \frac{{{\sigma _{{\rm{test}}}}\left( {\eta ,\varepsilon } \right)}}{{{\sigma _{{\rm{test}}}}\left( {{\eta _0},\varepsilon } \right)}} $

(13)

式中, σ_test(η, ε)为不同缺陷含量的HTPB推进剂单轴拉伸试验的应力值; σ_test(η₀, ε)为不含缺陷的HTPB推进剂试验应力值。

对不同缺陷含量的HTPB推进剂进行速率为5 mm·min^-1的单轴拉伸试验, 依式(13)对实验结果进行处理得到初始缺陷损伤因子f(η, ε), 如图 8所示。观察曲线发现其走势与反正切函数一致, 故选用反正切函数对f(η, ε)进行拟合, 即令:

$ f = A \times {\rm{arctan}}\left( {B \times \varepsilon + C} \right) + D $

(14)

拟合结果如图 8、表 5所示。

经观察, 发现参数A、B、C、D均随缺陷含量线性变化, 因此可选用线性函数对其进行拟合, 同时考虑到当缺陷含量为零时, 损伤因子f(η, ε) =0, 具体拟合结果见图 9。

将参数函数表达式代入式(12)中, 得到常温下含初始缺陷的HTPB推进剂粘超弹本构方程表达式:

$ \begin{array}{l} \sigma = \left\{ {1 - \left[ {A \times {\rm{arctan}}\left( {B \times \varepsilon + C} \right) + D} \right]} \right\} \times \\ \quad \left[ {({p_2} + {q_2}{\rm{ln}}\left( {\dot \varepsilon } \right)) \times \varepsilon + ({p_1} + {q_1}{\rm{ln}}\left( {\dot \varepsilon } \right))} \right] \times \\ \quad \sum\limits_{n = 1}^2 {\frac{{2{\mu _n}}}{{{\alpha _n}}}\left[ {{{\left( {1 + \varepsilon } \right)}^{{\alpha _n} - 1}} - {{\left( {1 + \varepsilon } \right)}^{\left( { - {\alpha _n}/2} \right) - 1}}} \right]} \end{array} $

(15)

最后分别用推进剂缺陷含量为10%及30%、拉伸速率为1, 10, 100 mm·min^-1的实验数据对本构方程进行验证, 结果如图 10所示。可见, 实验结果与预测结果吻合得较好, 二者最大偏差仅为4.4%, 说明了构建的本构模型具有一定的可靠性, 对含缺陷的HTPB推进剂力学性能的研究具有一定的参考意义。

5 结论

(1) 利用万能试验机, 对HTPB推进剂进行了不同速率的准静态单轴拉伸试验和多步松弛试验。实验结果表明, HTPB推进剂的在松弛阶段, 应力初始急剧下降, 而后松弛趋势逐渐变缓, 最后趋于稳定, 拉伸曲线率相关性明显。因此可以认为HTPB推进剂的力学特性具有典型的粘超弹特性。

(2) 在典型的超弹性本构模型基础上, 引入应变率参数函数, 构建了描述及预测HTPB推进剂单轴拉伸力学性能的粘超弹本构模型, 并拟合得到相应模型参数。通过10 mm·min^-1的验证实验, 验证了模型的准确性。

(3) 在粘超弹本构模型的基础上, 引入初始缺陷损伤因子, 并拟合得到了损伤因子的变化规律, 从而得到可描述含初始缺陷的HTPB推进剂力学性能的本构方程。采用不同缺陷含量不同拉伸速率下的实验结果对本构方程进行验证, 预测结果与试验结果最大偏差仅为4.4%, 验证了模型的可靠性。

需要说明的是, 本研究所建的含初始缺陷的HTPB推进剂粘超弹本构模型, 只考虑了推进剂内填料-基体界面的初始脱粘, 至于真实推进剂中可能含有的微孔洞、微裂纹以及颗粒夹杂等, 这些缺陷对推进剂力学性能的影响形式是否也如此仍需进一步研究。