非线性粘弹性本构模型在固体火箭发动机装药结构完整性分析中应用

李 毅,田维平,董新刚,姚 东

(中国航天科技集团公司四院四十一所，西安 710025)

0 引言

复合固体推进剂是一种颗粒增强的高聚型含能材料，广泛用于运载火箭、导弹等装备的动力系统，其性能直接影响装备的运载能力、射程、贮存性等关键性能。复合固体推进剂力学特性的粘弹属性与其自身的配方特点有关，也与服役环境如温度、加载应变率、应力应变状态等因素高度相关。在变形较小时，可认为其本构模型为线粘弹性；随着变形量增大，推进剂内部微裂纹、微孔洞等损伤逐渐演化，出现明显的脱湿现象，此时线性理论不再适用，应选取合适的非线性粘弹性本构模型。

早期基于复合推进剂的线粘弹性本构模型开展了工程领域的结构分析，结合产品研制的经验修正取得了良好的应用效果[1]。Scahpery[2]及其合作者一直致力于含损伤非线性粘弹性本构模型的研究。Ulrich Mandel等[3]对纤维增强复合材料进行研究，提出一种三维的非线性粘弹性本构模型，并指出此模型可适用于任意载荷情况下。Muliana[4]从细观角度入手，提出一种非线性粘弹性本构模型。张晓等[5]从应变能出发推导出HTPB推进剂非线性粘弹性本构模型。常新龙等[6]研究了HTPB高应变率下的非线性粘弹性本构模型。王哲君等[7]对HTPB压应力状态下的本构模型进行了研究。姚东等[8]从损伤效应的拟合出发，提出了考虑应力状态的非线性本构模型。唐志平、朱兆祥等[9]提出了被后来的学者称为朱-王-唐非线性本构模型的成果，均为复合固体推进剂非线性粘弹性本构模型的研究作出了贡献。

本文基于朱-王-唐非线性本构模型，针对其积分项繁多、本构参数获取难度较大的特点，提出一种形式简单、便于工程应用的本构模型。

1 推进剂本构模型的提出

朱-王-唐非线性本构模型如式(1)所示：

σ=E0ε+αε2+βε3+

(1)

式中E0、E1、E2、α、β、θ1、θ2均为材料常数，前3项表示非线性弹性部分，后2项积分项表示粘性部分。

式(1)表明，非线性粘弹性本构模型可写成弹性部分和粘性部分之和的形式。参考该结论，本文将推进剂本构模型表述为如下形式：

σ(t,ε)=fE(t,ε)+fv(t,ε)

(2)

式中fE(t,ε)为弹性项；fv(t,ε)为粘性项。

Staverman和Schwarzl[10]在大量实验观测基础上指出，采用式(3)：

σ(t)=E(t)f(ε)

(3)

代替σ(t)=E(t)ε可较好地描述实验结果。赵荣国[11]由此假设过去不同时刻的应变对现时应力的影响可忽略不计，并对其进行了验证。联想到应力松弛实验，本文认为推进剂的受力分为应变变化引起的应力及松弛现象导致的应力耗散。对于应变从0～ε的变化，都可认为是以上2种应力叠加产生的效果，这样应力的表达形式就可简单地写成：

σ(t,ε)=fE(ε)-σ松弛

(4)

式中fE(ε)、σ松弛分别为弹性、粘性应力项，其函数形式及其中的材料参数还有待通过实验数据进行确定。

2 基于松弛及拉伸实验的本构模型建模

推进剂在变形较小时，服从线粘弹性规律，当变形达到一定值后，损伤的逐渐演化导致其非线性显著增加，本文假设应力耗散项不受损伤变量的影响。因此，式(4)中的σ松弛就可用式(5)表示：

σ松弛=[E(0)-E(t)]ε

(5)

式中E(0)为初始时刻的松弛模量值；E(t)为t时刻松弛模量值；ε为t时刻应变。

粘性应力耗散项是状态量，仅与当前状态有关，如应变、时间等，不受以往应力-应变状态的影响。

弹性项的应力-应变关系包括大应变状态的数据，因此，还需进行一定的拉伸实验。

2.1 松弛实验及其分析

推进剂松弛实验采用哑铃模型，如图1所示。

试验件拉伸到应变为5%，随后测其应力值，并将数据转化为松弛模量。为得到完备的结果，需进行多组温度下的松弛实验。

根据时温等效原理E(t,T)=E(t/αT,T0)，对不同温度下的松弛模量进行处理，使所有数据基本重合为一条曲线，即松弛模量主曲线。随后，将松弛模量主曲线拟合为6阶Prony级数形式：

(6)

式中各参数的值如表1所示。在对不同温度下松弛模量进行处理时，得到不同温度对应的时温等效因子的值，时温等效因子lgαT仅和温度有关。根据WLF方程(7)，拟合得到的C1和C2的值分别为8.719 0、144.695 8。拟合后的松弛模量主曲线与实验数据的对比如图2所示。

(7)

2.2 拉伸实验及其分析

拉伸试件同图1所示模型，实验在Instron 4210型试验机上完成，测试温度 25 ℃，应变率 0.023 8 s-1，测得数据结果如图3所示。

式(4)中的应力弹性项不包括粘性应力。因此，根据拉伸实验得到的松弛模量及式(5)，对拉伸数据进行修正，得到弹性应力项的表达式：

(8)

fE(ε)=σ实验+σ松弛

修正后的应力用式(9)形式的函数进行拟合效果如图4所示。其中，E0和Es均为材料参数，不同温度下的E0和Es值如表2所示。

(9)

将E0和Es表示为温度T的函数关系式，拟合的函数关系式如式(10)、式(11)所示，拟合曲线如图5所示。

(10)

(11)

3 本构模型及其参数的数值验证

3.1 本构模型一维MATLAB验证

本文采用MATLAB编写代码对本构模型进行一维数值验证，MATLAB中仅模拟边界条件的变化，无需定义单元类型。在应变率0.023 8 s-1、25 ℃的情况下，计算出的应力-应变走势与实验数据的对比如图6所示。可看出，应力随应变变化的趋势基本和实验数据吻合，脱湿点对应的应变约为20%，脱湿后的误差明显大于线性段，相对偏差达8%。本文假设应力耗散与损伤无关，损伤的增加有可能会导致粘性应力耗散加大，具体规律需进一步研究。

粘弹性材料的力学性能具有温度相关性和载荷率相关性，需进行多组不同边界条件的数值仿真，现将不同条件下的数值仿真结果同实验的对比汇总如图7所示。从图7可看出，初始阶段应力应变曲线差别不大，符合线粘弹性规律；在相同拉伸速率条件下，随着温度的升高，应力值有一个减小的趋势，温度越低，应力水平越高；当温度不变时，拉伸速率越大，应力值越大。数值拟合结果与实验、理论都具有相同的规律，相对偏差不超过10%。

用MATLAB对松弛实验进行数值仿真，结果如图8所示。结果表明，从拉伸实验数据推导出的本构模型同样适用于非拉伸情况，且相对偏差不超过10%。

3.2 本构模型模型三维ABAQUS验证

利用ABAQUS用户材料子程序UMAT对本构模型进行二次开发，将本文提出的本构模型编写成UMAT子程序，在ABAQUS中进行本构模型的三维算例验证。

将本构模型从一维扩展到三维[12-13]如式(12)。

Aijkl(E(0)-E(t))εkl

(12)

其中，Aijkl为一个四阶张量，此处可表示为矩阵形式：

(13)

(14)

粘性项的应力计算可按照线弹性部分进行三维计算，需注意一点：粘性项是状态量，不可按照增量形式进行计算。

应力更新过程分为2步：(1)计算增量步结束时的弹性应力。t时刻传入的应力σ(t)为真实应力，欲得到增量步结束时的弹性应力，需先得到增量步开始时的弹性应力值，加上上一步减去的粘性应力项，才可得到增量步结束时的弹性应力。(2)计算增量步结束时的真实应力。用步骤(1)中得到的弹性应力减去增量步结束时的粘性应力，便可得到增量结束时的真实应力。增量步结束时真实应力为

(15)

则Δt时间段内应力增量为

Δσij(t+Δt)=σij(t+Δt)-σij(t)

=[Dij-E(0)+E(t+Δt)]AijklΔεkl(t+Δt)+

[E(t+Δt)-E(t)]Aijklεkl(t)

(16)

令：

(17)

依据上述有限元增量形式，将其编写为UMAT在ABAQ中进行数值仿真计算。

数值仿真模型尺寸与图1所示试件尺寸相同，在ABAQUS中建立模型如图9所示。模型包括试验件和夹具两部分。试验件网格尺寸1 mm，网格数量23 230，单元类型为C3D20H。夹具网格尺寸2 mm，与图9所示垂直纸面方向尺寸为1 mm，网格数量5410，材质为钢。计算过程中固定一个夹具，另一夹具匀速向另一方向移动。

推进剂泊松比取为0.496，单元类型选择杂交单元(Hybrid Formulation)以消除沙漏效应，应力更新过程有UMAT实现。

对不同温度、不同加载速率条件进行数值仿真，取试件中心节点为研究对象，查看其应力-应变曲线并与实验数据进行对比，结果如图10所示。三维数值仿真结果符合推进剂拉伸实验结果，不同温度、拉伸速率下的变化规律也符合理论预测。

从图10可看出，三维结果较一维结果误差有所增大，相对偏差最大为14.8%，这是因为三维计算在ABAQUS下进行，所用应力应变为真实应力应变，较名义应力应变误差增大。

误差原因分析：

(1)单向拉伸试验得到的应力-应变均为名义值，转换为真实值，进而拟合本构参数时存在一定偏差；

(2)泊松比按经验取为0.496，与真实值存在一定偏差；

(3)粘性应力耗散项的计算由松弛模量确定，文中假设这一项与损伤无关，而损伤对推进剂本构模型影响极大。因此，在损伤较大的情况下，松弛模量的形式不变必然导致误差的存在。

4 模拟发动机热-结构分析

根据本文提出的非线性粘弹性本构模型，对贴壁浇注固体火箭发动机硫化降温过程进行分析。计算模型如图11所示，主要结构参数采取无量纲表示。

在ABAQUS/CAE中建立模型。药柱、绝热层、壳体之间采取布尔加运算进行合并。为控制网格对分析结果的影响，对关键部位采取了不同的种子控制策略(尺寸及过渡要求)，并对药柱、绝热结构等对应材料泊松比较大的部位采用杂交单元技术以消除沙漏效应，如图12所示。在药柱两侧施加对称边界条件，在后裙端面施加固支条件以消除刚体位移。

推进剂、绝热层、壳体各部件的材料属性如表3所示。推进剂模量部分本文提出的本构模型在UMAT进行计算。

热边界：药柱硫化降温过程是放在一定的温度环境中，该环境温度为10 ℃；药柱、绝热层、壳体初始温度均为60 ℃，推进剂和绝热层与空气的对流换热系数为1 K/(m2·W)，壳体与空气的对流换热系数为10 K/(m2·W)，时间历程设置为30 d。

表3 推进剂、绝热层、壳体的材料属性

药浆浇注后的硫化过程涉及粘合剂基体高分子网络的交联、粘合剂基体与固体颗粒组分的粘接等一系列化学-热力学耦合过程，反应机理复杂、影响因素众多。本文所述硫化降温过程，以药浆硫化反应达到“零应力温度”为起点，并假设：

(1)“零应力温度”直至目标环境温度的过程中，药柱热物理属性及力学性能均可由推进剂试件的测试结果进行有效表征；

(2)在“零应力温度”直至目标环境温度范围内，推进剂热物理属性不随温度变化。

受限于实验测试条件，在进行硫化降温实验中，仅可测试药柱中孔处的位移变化，因此选取此作为本构模型的验证标准。

取模型上3个典型位置如图13所示，其输出温度随时间变化的数据如图14所示。外壁面与目标环境的对流换热系数最大，故温度下降快；药柱内孔处与目标环境的对流换热系数较小，因此温度下降较为缓慢；肉厚中间处热损失最少，故温度下降最为缓慢。

硫化降温过程结构的变化主要是因为温度场的变化引起的，初始阶段温度变化剧烈，所以中孔位移变化也会较快，当温度场平衡后，位移也不会发生变化。从图14可看出，约200 h后，温度场趋于平衡。因此，中孔位移必然也有类似的规律。将中孔位移随时间变化的数据绘制在图15中，数值仿真结果与分析一致。

数值仿真计算得到的中孔位移最大为0.014 48R(R为壳体半径)。采用相同结构尺寸的发动机进行了10 ℃环境下的温度平衡试验，8 d后测得药柱内孔处温度与环境温度一致，表明药柱温度达到了平衡；此时进行药柱内径测量，并与设计值对比，发现药柱内孔径向位移约为0.013 8R。数值仿真与实验数据的相对偏差仅为4.9%，验证了本构模型的准确性。

5 结论

(1)本文构建的考虑损伤、温度的非线性粘弹性本构模型，由弹性项和粘性项组成；弹性项仅为应变的函数，结合松弛实验和拉伸实验对参数进行了拟合；粘性项为状态量，且不受损伤的影响，可由松弛模量和应变计算得到。

(2)通过MATLAB和QBAQUS对本构模型进行了数值仿真验证，仿真结果与实验数据吻合较好，一维仿真结果与实验数据相对偏差不超过10%，三维仿真结果相对偏差略有增大，最大为14.8%。

(3)基于本文所建立的本构模型对发动机硫化降温过程进行了数值仿真，得到中孔径向位移为0.144 8R，与实验测得的数据(0.013 8R)相对偏差仅为4.9%，说明本构模型在发动机热-结构分析具有较好的精度。

[1] 刘中兵,利凤祥,李越森,等.高过载条件下固体推进剂药柱结构完整性分析计算[J].固体火箭技术,2003,26(2):12-16.

LIU Zhongbing,LI Fengxiang,LI Yuesen,et al.A calculation analysis for structural integrity of solid propellant grains under high overload[J].Journal of Solid Rocket Technology,2003,26(2):12-16.

[2] Ha K,Schapery R A.A three dimensional viscoelastic constitutive model for particulate composites with growing damage and its experimental validation[J].International Journal of Solid Structures,1998,35(26-27):3497-3517.

[3] Ulrich Mandel,Robin Taubert,Roland Hinterh?lzl.Three-dimensional nonlinear constitutive model for composites[J].Composite Structures,2016,142:78-86.

[4] Muliana A,Rajagopal K R,scharnuterD T,et al.A nonlinear viscoelastic constitutive model for polymeric solids based on multiple natural configuration theory[J].International Journal of Solid Structures,2016,100-101:95-100.

[5] 张晓,郑坚,彭威,等.HTPB复合固体推进剂粘弹性应变能及非线性本构模型[J].固体火箭技术,2015,38(6):827-832.

ZHANG Xiao,ZHENG Jian,PENG Wei,et al.Viscoelastic strain energy and nonlinear constitutive model for HTPB composite solid propellant[J].Journal of Solid Rocket Technology,2015,38(6):827-832.

[6] 常新龙,赖建伟,张晓军,等.HTPB推进剂高应变率粘弹性本构模型研究[J].推进技术,2014,35(1):123-127.

CHANG Xinlong,LAI Jianwei,ZHANG Xiaojun,et al.High strain-rate viscoelastic constitutive model for HTPB propellant[J].Journal of Propulsion Technology,2014,35(1):123-127

[7] 王哲君,强洪夫,王广,等.中应变率下HTPB推进剂压缩力学性能和本构模型研究[J].推进技术,2016,37(4):776-782.

WANG Zhejun,QIANG Hongfu,WNAG Guang,et al.Mechanical properties and constitutive model for HTPB propellant under intermediate strain rate compression[J].Journal of Propulsion Technology,2016,37(4):776-782.

[8] 姚东,张光喜,高波.考虑应力状态的HTPB/AP推进剂含损伤热-粘弹性本构方程[J].固体火箭技术,2014,34(4):496-499.

YAO Dong,ZHANG Guangxi,GAO Bo.Constitutive equations involving damage for HTPB/AP propellant considering stress state[J].Journal of Solid Rocket Technology,2014,34(4):496-499.

[9] 唐志平,田兰桥,朱兆祥.高应变率下环氧树脂的力学性能[C]//全国第二届爆炸力学会议文集.扬州,1981.

TANG Zhiping,TIAN Lanqiao,ZHU Zhaoxiang.Mechanical properties of epoxy resin under high strain rate[C]//Anthology of the Second National Conference on Mechanics of Explosion.Yangzhou,1981.

[10] Staverman A J,Schwarzl F.Nonlinear deformation behavior of high polymers[M].Stuart H A,Ed.Die Phyaile der Hochpolymeren,1956.

[11] 赵荣国,张淳源.一个考虑损伤的非线性粘弹性本构关系[J].湘潭大学自然学科学报,2001,23(3):28-33.

ZHAO Rongguo,ZHANG Chunyuan.A nonlinear viscoelastic constitutive relation with damage[J].National Science Journal of Xiangtan University,2001,23(3):28-33.

[12] 冯震宙,王新军,王富生,等.朱-王-唐非线性粘弹性本构模型在有限元分析中的实现及其应用[J].材料科学与工程学报,2007,25(2):269-272.

FENG Zhenzhou,WANG Xinjun,WANG Fusheng,et al.Implementation and its application in finite element analysis of constitutive model for ZWT nonlinear viscoelastic material[J].Joural of Materials Sciense and Engineering,2007,25(2):269-272.

[13] 彭云.基于 ABAQUS 的非线性粘弹性本构模型二次开发[J].南昌航空大学学报,2011,25(2):38-41.

PENG Yun.Developing of nonlinear viscoelastic constitutive model based on ABAQUS[J].Joural of Nanchang Hangkong University (Natural Sciences),2011,25(2):38-41.