基于分段曲线拟合的HTPB固体推进剂损伤数值模拟

开篇：润墨网以专业的文秘视角，为您筛选了一篇基于分段曲线拟合的HTPB固体推进剂损伤数值模拟范文，如需获取更多写作素材，在线客服老师一对一协助。欢迎您的阅读与分享！

摘要： 通过对htpb固体推进剂单轴试验所得的结果进行分段曲线拟合，并对损伤变量进行分段定义，将其引入到HTPB固体推进剂材料的刚度矩阵中，采用Abaqus/Standard中的用户子程序UMAT接口进行编程计算，对材料的损伤特性进行有限元模拟.模拟结果接近于实际结果，可以应用于固体发动机的工程设计中.

关键词： HTPB； 本构方程； 分段曲线拟合； 有限元； 损伤变量； 连续介质损伤力学

中图分类号： V512；TB151.1文献标志码： B

0引言

复合固体推进剂是一种高填充比颗粒增强复合材料，主要由以高分子化合物为母体的黏合剂和掺入其中的大量固体氧化剂颗粒及金属燃料颗粒（如铝粒子等）组成，此外还有少量的附加组分，见图1.黏合剂为黏弹性的橡胶材料，为粒子提供垫子和高氯酸铵，相对于弹性体的黏合剂，主要表现为弹性.铝粉由于具有金属晶体的一些性质，表现出弹塑性.在外载荷作用下，如此复杂的一个复合体不断变形，必然表现出弹性和黏弹性的特征，其中一个力学特性就是拉压与拉剪强度不等.如果考虑大的氧化剂颗粒从黏合剂体系脱开，即所谓的“脱湿”，会形成微裂纹和空洞.事实上，这些微观结构的损伤足以导致推进剂表现出宏观非线.HTPB固体推进剂材料的损伤行为可由BIEGLER和MEHRABADI提出的基于能量的本构模型进行描述[12].

图 1推进剂的主要组分和典型物理模型

对复合固体推进剂损伤性质的研究，必须建立在试验的基础之上，国内外学者在这方面发表过很多文章[36]，但从应用角度看，多是理论强于实践，需要的测试参数很多，在工程上应用起来有很多不便之处.

在单轴拉压条件下的试验结果需通过曲线拟合得到相应的应力应变曲线.若该曲线在一个比较广泛的应变域中，很难用一个统一的函数对其进行精确描述.本文考虑将所描述的应变域进行分段，在每段中以不同的简单函数进行拟合，使固体推进剂的力学行为模拟结果更接近于试验结果.

1本构方程和损伤变量

1.1本构方程描述

根据MEHRABADI和COWIN的表述：（1）对于任意弹性对称体，全应力张量和全应变张量可以分解为6个或少于6个的等同形式的特征张量之和；（2）对于任意弹性对称体，每个特征应力张量直接与其等同形式的特征应变张量成正比；（3）总应变能可分解为6项或少于6项之和，每项都是等同形式的特征应力与特征应变的标量积，表示6个互不作用的能量模态.

所建立材料模型只适用于平面应力情况.在没有损伤影响下，材料在单向应力状态下应力与应变之间的本构关系是线性的.对于正交各向异性的复合材料，应力与应变之间的弹性本构关系为T=C：E（1）式中：C为2阶刚度张量，其矩阵形式可写为

与其他各向异性材料相同，HTPB推进剂表现为拉压和剪切性质的不同，因此式（2）是假定材料为正交各向异性后所得出的沿2个不同方向的初始值.E11和E22可考虑在2个正交方向上取拉伸和压缩试验的平均值，υ12和υ21同样可以通过试验确定.根据文献，要描述复合材料的正交各向异性，应该分解应力张量和应变张量.在平面应力状态下，全应力和全应变分解成对应其特征值的3个特征张量.

由上述第二个性质可知，第i个特征值λi，特征应力与特征应变之间的关系可写为Tci=λiEci（3）由第三个性质，可得2ρ=3iTciEci（4）式中：ρ为总应变能密度.对于这些非相互作用的能量模态，可用于建成立失效准则.

1.2损伤变量的定义

根据连续介质损伤力学可以假定，复合固体推进剂材料的非线性特性完全是由材料的损伤引起的.当损伤变量D=0时，材料是完全弹性的；当D=1时，材料完全损伤.对于线弹性材料，可以证明按材料刚度的改变定义损伤变量与连续介质损伤力学中按横截面损伤的比值定义损伤变量是等价的，因此若损伤之前材料的刚度为C0，损伤后的刚度为C（0

2应力应变曲线的分段拟合

通过对复合固体推进剂材料的试样进行拉伸、剪切或者压缩试验，可以得到固体推进剂材料的力学性质.若能对应力应变试验值进行较简单的拟合，得到一个统一表达的函数，则可以使上述编程过程变得相对简单.但是，在一个较大的应变域中，可能很难找到一个能够描述试验形态的统一函数.因此，可以考虑用分段拟合的方法得到相应的应力应变关系曲线.

设第r（r=1，2，…）组试验数据，取线性无关组φ（i）ij（ε）=εi-1r

i=1，2，…，m，j=1，2…，n（21）由m次最小二乘拟合，可得多项式σr（ε）=mk=0Crkφrk（ε）=mk=0Crkεk-1r（22）将式（22）代入式（15）中，可以得到在不同区域的范围内，对损伤变量的不同定义为Dr=I-1CroTr（E）E（23）式（23）中需要确定Cro，在增量计算中Cro为分段区域端点上2条拟合曲线所共同具有的切线斜率，因此在分段拟合中，保持2条拟合曲线在端点处满足条件Tr-1（E-0）=Tr（E+0）=k0

Tr-1（E）EE=E-0=Tr（E）EE=E+0（24）k0需满足在计算过程中的稳定性要求.将计算所得的Dr代入（21）式，通过修改Abaqus/Standard中UMAT的DDSDDE进行计算.

某配方常温下HTPB固体推进剂应力应变曲线见图2.

图 2HTPB固体推进剂典型应力应变曲线

局部放大后，当ε

表 1HTPB推进剂材料性质参数εmaxab数值102.13.2

但是，当ε>10%以后，用双曲正切函数并不能很好地拟合试验曲线，因此在此区间，可采用最小二乘法进行拟合，见表2.

表 2HTPB固体推进剂的材料性质第二段试验数据应变102555应力/MPa0.410.690.92

由表2 拟合出如下函数.σ（ε）=0.187 9ε0.6-0.024 92ε（26）根据式（26）可得出在ε>10%后的损伤变量定义为D2=1-175.40（0.112 7ε-0.6+0.024 92）（27）3有限元模拟和结果讨论

推进剂材料使用圆柱型试件，其尺寸为Φ25 mm×50 mm圆柱型裂纹试件，两端采用专用夹具与拉伸机连接.程序开头片段见图3.

图 3程序开头片段

计算结果均通过von Mises应变εe和应力σe进行分析.试件的整体应力云图见图4；推进剂材料在破坏前的整体应力云图见图5；载荷时间分别为推进剂在1，2.5，7.5和15 s时的剖面应力云图见图6.图 4含夹具试件的整体应力云图

图 5推进剂的整体应力云图

图 6推进剂在载荷时间分别为1，2.5，7.5和

15 s时的剖面应力云图

根据应力计算分析的结果，推进剂发生破坏的位置在推进剂材料的中心位置，随着载荷时间的增加，推进剂内部的应力逐渐增大.若要发生破坏，则从推进剂最中心的位置产生，从这个位置扩展，向外层表面延伸，直到试件完全破坏.

由图6可知，推进剂表现出明显的非线性黏弹性特征.随着位移的增加，应力呈曲线上升；因为拉伸速率为定值，所以应变为线性发展趋势，这与试验情况一致.材料单轴加载试验表明，在ε

4结束语

依据试验结果，通过对单轴试验曲线的适当分段拟合，用有限元法进行模拟，可以对HTPB推进剂材料的损伤行为做出较好解释，其模拟结果与试验也有较好的一致性.参考文献：

[1]COWIN S C， MEHRABADI M M. The structure of the linear anisotropic elastic symmetries[J]. Mech Phys Solids， 1992， 40（7）： 14591471.

[2]ARRAMON Y P， MEHRABADI M M， MARTIN D W， et al. A multidimentional anistropic strength criterion based on kelvin modes[J]. Int J Solids & Struct， 2000， 37（21）： 29152935.

[3]沈观林， 胡更开. 复合材料力学[M]. 北京： 清华大学出版社， 2006： 2944.

[4]李灏. 损伤力学基础[M]. 济南： 山东科学技术出版社， 1992： 4349.

[5]LEMAITRE J. A continuous damage mechanics model for ductile fracture[J]. J Eng Mat & Technol， 1983， 107（1）： 8389.

[6]LOCKYEAR S A. Mechanical analysis of transversely loaded wood/plastic sections[D]. Pullman： Washington State University， 1999.