改性双基推进剂非线性本构模型及其数值实现.pdf

第 卷第 期 固 体 火 箭 技 术 改性双基推进剂非线性本构模型及其数值实现张 旭，刘向阳，王士欣，闫晴霄，张广龙，张 源（北京理工大学 宇航学院，北京；中国船舶重工集团公司第七一三研究所，郑州）摘要：为分析改性双基（）推进剂在复杂载荷下的力学响应，建立了一种含损伤变量的 推进剂非线性粘弹性本构模型。鉴于 推进剂的拉压不对称性，通过单轴拉伸、单轴压缩以及对应的松弛试验获得本构模型中的关键参数。当使用含应变和应变率的双指数形式表征与应变率呈“负相关”的损伤变量时，推进剂的非线性力学行为能被准确地描述。同时，引入体积应变作为拉伸和压缩条件的判定变量，发展出基于 （）子程序的非线性本构模型数值实现方式。将该考虑拉压不对称性的本构模型及数值实现方式应用于燃气发生器药柱在点火冲击载荷作用下的力学响应分析，并与单一拉伸参数分析结果进行对比。结果表明，考虑拉压不对称性时，结构的等效应变和损伤变量与使用单一拉伸参数时的分析结果存在显著差异。因此，判断结构的拉压状态并选用正确参数具有重要意义。关键词：改性双基推进剂；本构模型；拉压不对称性；体积应变；中图分类号：文献标识码：文章编号：（）：，（，；，）：（），：；引言复合改性双基推进剂（）是一种经典的固体推进剂，在导弹等武器系统中有广泛的使用。在硝化纤维素、硝化甘油中加入了固体氧化剂等组分后，其性能有了明显提升，但同时其内部组分间的多相界面使其力学性能变得更为复杂。刘远祥等进行了不同拉伸速率下 推进剂的单轴拉伸试验，结果表明，推进剂的力学行为明显受到应变率的影响，收稿日期：；修回日期：。作者简介：张旭，男，博士生，研究方向为固体推进剂本构模型。通讯作者：刘向阳，男，教授，研究方向为固体火箭发动机寿命预估。随着应变率的增加，推进剂的应力应变曲线逐渐出现“强化现象”。孙朝翔在进行不同应变率下的 压缩性能试验后发现，的压缩初始模量和屈服应力随着应变率的升高而升高。王鸿丽等在分别进行 拉伸和压缩试验后提出，具有明显的应力和应变拉压不对称性，其认为拉压不对称性是由材料初始缺陷的扩展等因素造成的。综上可知，推进剂的力学行为存在率相关性和拉压不对称性。在此基础上，发展出针对 推进剂力学行为的本构模型和数值实现方式具有重要的意义。针对 推进剂的本构模型及其数值实现方式，国内外研究者展开了很多研究。孟红磊等建立了 推进剂的非线性粘弹性本构模型，利用 将本构模型二次开发后，其分析了压力、星孔参数及包覆层材料对装药结构完整性的影响。沙宝林等建立了 推进剂统一本构，其中损伤水平由恢复应变能密度的方法预测，在 二次开发后，对燃烧室力学行为进行了预测。以上本构模型具有参考价值。在加入固体填料后，推进剂与其他固体推进剂的力学性能较为类似。因此，在建立适用于 推进剂的本构模型时，也可以参考其他类型固体推进剂等粘弹性材料的相关研究。等通过引入应变软化函数构建了高伸长率下的唯象型非线性本构模型，该软化函数类似于宏观损伤演化函数。许进升基于端羟基聚丁二烯（）推进剂单轴试验和松弛试验，发展出基于 级数的损伤粘弹性本构模型并应用于发动机结构完整性分析。等使用了一种简单但形象的损伤变量形式来构建 推进剂非线性本构模型，其损伤定义为衰减应力与线性粘弹性应力的比值。以上非线性本构模型均通过引入损伤变量（或软化函数）构建，需要在此基础上建立适用于 材料的率相关性和拉压不对称性的损伤形式。损伤变量作为构建非线性粘弹性本构模型的关键参数，其演化规律对推进剂力学行为的描述意义重大。但目前还缺乏对 推进剂损伤变量与应变率和拉压条件相关性的研究。本文建立了改性双基推进剂含损伤的非线性粘弹性本构模型并推导了其增量形式。在通过单轴和松弛试验拟合得到本构模型的关键参数后，基于拉压不对称性发展出了本构模型的数值实现方式。最后，利用该本构模型和数值实现方式对点火冲击载荷下的燃气发生器药柱力学响应做出分析。本构模型根据文献，固体推进剂的非线性本构模型可定义为（）（）（）（）（）式中（）为应力；（）为应变；（）为松弛模型；为积分变量。损伤变量 定义为衰减应力所占线粘弹性应力的比值，其取值范围为，当 时表示材料无损伤，当 时表示材料破坏。当使用由广义 模型得到的 级数形式松弛模量作为积分型本构模型中的松弛模量时，得到的非线性本构模型为（）（）()（）（）式中（）为 级数形式松弛模量；为松弛平衡模量；为第 组 模型的弹性模量；为松弛时间。如果应变率为恒定值，其可表示为（）（）代入式（）得到积分后恒应变率下的本构模型：（）（）()（）为后续分析便捷，定义线性粘弹性应力为 ()（）以上得到一维条件下全量型非线性粘弹性本构模型。在材料的变形过程中，很难保持恒定的应变率。因此，需要推导本构模型的增量形式，并对应变率进行更新。在 时刻，材料的损伤变量为（）；在 时刻，材料的损伤变量为（）。则在时间增量 内，材料的应力增量略去高阶项后表示为（）（）（）（）（）（）（）（）（）（）当时间增量 足够小时，在 内，材料的应变率可表示为（）（）（）（）线性粘弹性应力增量为（）（）（）（）（）()（）因此，式（）可表示为（）（）（）（）()（）（）（）()（）年 月固体火箭技术第 卷 以上得到的增量形式本构模型中，应变率可随着时间不断更新，故该增量型本构模型能预测复杂应变率下材料的大变形粘弹性力学行为。试验与参数拟合 单轴试验及松弛试验针对 推进剂存在的拉压不对称性问题，为准确分析拉伸和压缩条件下材料的力学响应，本研究中分别进行拉伸和压缩条件下的单轴和松弛试验来获得所建立非线性粘弹性本构模型的参数。单轴拉伸和拉伸松弛试验使用哑铃板材试样，单轴压缩和压缩松弛使用圆柱棒材试样，哑铃试样和圆柱试样尺寸如图 所示，单位为。单轴拉伸和压缩试验在万能试验机中进行，试验温度为 。拉伸速率分别选取、，对应的应变率分别为、；压缩速率分别选取、，对应的应变率分别为、。每种条件进行 组试验，选取 组试验真实应力和应变的平均值绘制应力应变曲线，拉伸和压缩条件下应力应变曲线如图 所示。R12（）2030（）图 试件尺寸图 图 中，当应变较小时，材料处于线性粘弹性范围；当应变达到某一数值时，应力应变曲线出现一个模量突减的现象，该模量突减的点被称为“损伤点”，该点对应的应变值被称为“损伤阈值”。且拉伸和压缩应力应变曲线中，随着应变率的增加，应力有明显的增加，即 推进剂存在明显率效应。拉伸和压缩松弛试验中，试验温度为 ，初始恒定应变为，松弛时间为 ，进行 组试验。一般认为，采用 阶 级数能准确描述 的松弛性能，将松弛模量和松弛时间按照式（）拟合为 阶 级数形式，拟合结果如表 所示。6.6710 s-4-1-3-13.3310 s-2-11.6710 s-2-13.3310 s（）1.3910 s-3-1-3-15.5510 s-2-12.7810 s-2-15.5510 s（）图 单轴试验应力应变曲线 表 松弛模量参数 损伤变量拟合利用式（）计算得到 加载速率下拉伸及压缩条件下损伤变量演化曲线如所图 所示。由图 可知，在未产生损伤时，线性粘弹性应力应变曲线与真实应力应变曲线重合度较好，此时损伤值为。在产生损伤后，真实应力明显小于线性粘弹性应力。损伤变量随着应变的增加而增加，直至达到一个数值并保持稳定。在这一模量稳定的区间内，损伤变量也保持稳定。其后，当应变达到最大伸长率后，损伤变量迅速增加至。年 月张旭，等：改性双基推进剂非线性本构模型及其数值实现第 期损伤变量随应变近似满足指数分布。拉伸条件下 推进剂的损伤变量在初始阶段随应变的增加速率小于压缩条件，且拉伸损伤变量的最终稳定值也小于压缩损伤变量的稳定值。由应力应变数据计算得到的不同应变率下的损伤变量如图 散点所示。由图 可知，随着应变率的增加，损伤阈值有所增加，但增加的幅值逐渐减小，呈指数型演化。同时，在不同应变率下，最大伸长率处拉伸损伤变量数值（图 最终值）变化的范围宽于压缩情况。整体而言，推进剂的损伤变量与应变率呈“负相关”。且拉伸条件下相关性强于压缩条件。Stree/MPaStree/MPa（）（）图 加载速率下损伤变量演化曲线 Exp.1.3910 s-fit-3-1-3-15.5510 s-fit-2-12.7810 s-fit-2-15.5510 s-fitExp.6.6710 s-fit-3-1-3-13.3310 s-fit-2-11.6710 s-fit-2-13.3310 s-fit （）（）图 双指数型损伤变量与试验结果对比 由式（）相关理论可知，当应变率增加时，线性粘弹性应力呈比例增加，损伤变量计算的分母更大。此外，当应变率较高时，推进剂的模量更高，更不容易发生损伤。以上两点原因共同导致损伤变量与应变率呈“负相关”特性。在应变较大时，随着应变的增加，损伤变量数值有所减小，这是由于推进剂材料在非线性段存在一定的“强化现象”，即模量随应变增加而增加的现象。当模量增加时，损伤变量数值会减小。文献中，将与应变率呈对数关系、与围压强度呈指数关系的损伤变量表征为对数与指数结合的形式。本文中，将与变形状态与应变率均呈指数关系的损伤变量表征为双指数形式。当 时，损伤变量定义为 （）（）（）式中 为损伤阈值；、和 为损伤变量中的参数。当 时，损伤变量 。其中，包含应变率的指数项用于表征损伤变量与应变率的“负相关”性，即 应为负值。损伤阈值表达式为 （）（）式中、和 为损伤阈值参数。由双指数损伤变量模型拟合得到的损伤变量及损伤阈值参数结果如表 所示。年 月固体火箭技术第 卷双指数形式的损伤变量计算结果如图 实线所示。由图 对比可知，双指数形式的损伤变量与实验结果具有很好的一致性（拉伸损伤变量拟合相关系数为 ，压缩损伤变量拟合相关系数为 ）。表 损伤变量及损伤阈值参数 数值实现利 用中 的（）用户材料子程序建立本构模型的数值实现方法，子程序中进行应力状态判定和应力更新的流程如图 所示。MisesthMisesth图 应力更新流程图 应力更新过程中，关键的部分有：（）变形状态判别因为 材料存在拉压不对称性，因此，必须对单元所处的拉伸和压缩状态进行区分。对单元受力状态的判断采用体积应变来判断，体积应变指物体单位体积的改变量。三维应变单元的体积应变为 （）式中 为体积应变，其数值与第一应变不变量 相同；为变形前体积；为变形后体积；、和 分别为沿、和 方向的应变。（）应变率计算在本文的研究中，选用 等效应变作为损伤状态的判断量，选用 等效应变率用于损伤阈值、损伤变量和松弛修正系数的计算。等效应变基于 子程序获得。使用等效应变率可以避免因为选取不同方向的应变而导致节点出现变形状态的判断差别。等效应变率的计算表达式为（）（）（）式中 为 应变率；、和 分别是、和 面内的剪应变。本构模型应用 有限元模型在建立非线性本构模型及其数值实现方式后，本节对某型燃气发生器单根药柱结构在点火冲击载荷下的力学响应进行分析。燃气发生器单根药柱结构如图 所示，包括药柱、包覆层、挡药板和弹性垫圈等部件。BaffleCoatingGasketGrain图 燃气发生器单根药柱结构 该型燃气发生器单根药柱长度约为 。点火冲击载荷基于燃气发生器常温测试试验得到，在 内，药柱内外表面压强线性增加至 。分析中使用到的材料参数如表 所示。表 燃气发生器药柱组件材料参数 （）（）药柱结构划分 网格，共 个单元；其他结构划分 网格，共 个单元。挡药板外表面施加固定约束边界条件，点火冲击载荷作用在所有部件的外表面，不同部件接触面施加绑定约束。年 月张旭，等：改性双基推进剂非线性本构模型及其数值实现第 期 结果分析点火冲击载荷作用后，药柱结构体积应变云图如图 所示。由图 可知，在点火冲击载荷作用下，整个药柱的体积应变均为负值，即药柱处于压缩状态。因此，在对结构进行力学响应分析时，自动选择压缩本构参数。同时，由体积应变云图可得，药柱两端面处变形最大。点火冲击载荷作用后，药柱端面等效应力应变云图如图 所示。可见，由于与挡药板的接触压力，挡药板端的应力分布较为复杂。挡药板端最大等效应力为，最大等效应变为。图 药柱体积应变 3.7043.4013.0992.7972.4942.1921.8901.5871.2850.9830.6800.3780.0760.0100.0090.0080.0070.0060.0050.0040.0030.0020.0010.000-0.001（）（）8.3217.6346.9476.2605.5734.8864.1983.5112.8242.1371.4500.7630.0760.0310.0290.0260.0240.0210.0190.0160.0140.0110.0080.0060.0030.001（）（）图 药柱等效应力应变 相较于挡药板端，垫圈端药柱的等效应力较为规整，药柱的应力呈圆环状分布。中间部位应力应变数值较大，表面和内孔处数值较小，这是由于内外表面共同向中间部位挤压造成。垫圈端最大等效应力为，最大等效应变为。药柱的最大力学响应在垫圈端，故对垫圈端的损伤做出分析，含包覆层药柱垫圈端损伤变量分布如图 所示。由图 分析可知，药柱结构仅有垫圈端药柱中心区域发生了损伤，最大损伤变量数值为，即该部分区域进入非线性粘弹性范围。除靠近端面部位外，其余大部分部位仍处于线性粘性粘弹性范围，即未发生损伤。为验证考虑拉压不对称性的必要性，采用拉伸本构参数对点火冲击载荷作用下药柱的力学响应进行分析并与前序结果对比。当不考虑拉压不对称性，仅使用拉伸本构参数用于结构力学响应分析时，计算得到药柱的等效应力应变分布情况与考虑拉压不对性时相同，但数值有所差别。0.1570.1440.1310.1180.1050.0920.0780.0650.0520.0390.0260.0130.000D图 垫圈端损伤变量 药柱的最大等效应力和等效应变分别为 和。分析得到的等效应力数值没有明显差别，但最大等效应变显著大于考虑拉压不对称性时分析结 年 月固体火箭技术第 卷果（最大等效应力和等效应变分别为 和，见图）。这是由于在应变较小时，拉伸条件下需要更大的应变才能达到压缩条件下较小应变的应力水平。拉伸参数分析得到的垫圈端损伤云图如图 所示。由图 可知，使用拉伸参数计算得到的损伤变量最大值为，远小于考虑拉压不对称性时的损伤变量最大值（，见图）。0.0230.0210.0190.0170.0150.0130.0120.0100.0080.0060.0040.0020.000D图 基于拉伸参数的垫圈端损伤变量 以上现象由两个因素导致：（）在较高应变率下，拉伸损伤阈值大于相同应变率下的压缩损伤阈值；（）在发生损伤后，拉伸损伤变量的增加幅度慢于压缩损伤变量（见图）。综上所述，在对复杂载荷作用下的结构进行力学响应分析时，考虑 推进剂的拉压不对称性，判断结构的变形状态并选择准确的参数具有必要性。结论本文构建了改性双基（）推进剂非线性本构模型并发展了其数值实现方式，得到的主要结论如下：（）通过单轴试验和松弛试验构建了含损伤变量的非线性本构模型。试验结果表明，在应变小于损伤阈值时，损伤变量接近于；当应变大于损伤阈值时，损伤变量随应变呈指数型增加。（）推进剂的损伤变量与应变率呈“负相关”，拉伸条件下相关性强于压缩条件。使用包含应变和应变率的双指数可以准确描述与变形状态与应变率相关的损伤变量。（）引入将体积应变作为拉压状态的判定变量，基于 子程序完成了考虑拉压不对称的 推进剂本构模型数值实现。（）将数值实现方法应用于燃气发生器药柱在点火冲击下的力学响应分析。与使用单一拉伸参数时进行力学响应分析的结果相比，考虑拉压不对称性时结构的等效应变更小，但损伤变量更大。考虑拉压不对称性并选用准确参数具有重要意义。参考文献：，：，（）：，（）：刘远祥，胡少青，张亚俊，等不同拉伸速率下改性双基推进剂的力学特性推进技术，（）：，（）：孙朝翔宽泛应变率和温度下改性双基推进剂本构模型及应用研究南京：南京理工大学，：，王鸿丽，许进升，周长省，等改性双基推进剂拉压不对称力学性能研究固体火箭技术，（）：，（）：孟红磊改性双基推进剂装药结构完整性数值仿真方法研究南京：南京理工大学，：，孟红磊，鞠玉涛含损伤非线性粘弹性本构模型及数值仿真应用固体火箭技术，（）：，（）：沙宝林，侯晓压力环境下固体推进剂含损伤的统一本构研究强度与环境，（）：，年 月张旭，等：改性双基推进剂非线性本构模型及其数值实现第 期 ，（）：，（）：许进升复合推进剂热黏弹性本构模型实验及数值仿真研究南京：南京理工大学，：，：周东谟定应变下 推进剂老化机理及发动机寿命预估研究北京：北京理工大学，：，：，：赵伯华，沈月萍火箭装药体积形变力学特性的研究推进技术，（）：，（）：（编辑：吕耀辉）年 月固体火箭技术第 卷