含脱湿的推进剂动态热粘弹性本构模型的构建和应用方法

摘要

本发明公开了一种含脱湿的推进剂动态热粘弹性本构模型的构建和应用方法，模型的构建方法实施步骤包括：建立率相关本构方程，确立脱湿因子、温度修正系数，确立含脱湿因子、温度修正系数的率相关本构模型；模型的应用方法步骤包括构建含脱湿因子、温度修正系数的率相关本构模型并推导模型的增量格式；使用有限元软件的扩展功能接口对增量格式的含脱湿因子、温度修正系数的率相关本构模型进行编程，得到用户子程序，使用有限元软件调用用户子程序。本发明的模型构建方法能更精确的预测粘弹性材料结构在冲击和振动等动态载荷下的力学行为，同时对构建的本构模型进行数值实现，可在有限元软件中对粘弹性材料结构进行准确的动力学分析。

说明书

技术领域

本发明涉及粘弹性材料本构模型研究领域，具体涉及一种含脱湿的推进剂动态热粘弹性本构模型的构建和应用方法，涉及固体推进剂率相关本构模型的构建，用于以固体推进剂为代表的需考虑细观损伤和温度效应的粘弹性材料的本构模型的构建和应用。

背景技术

固体推进剂药柱作为固体火箭发动机的能量来源和主要部分，在其全寿命周期中需要承受温度、内压、振动、冲击等复杂载荷，其中的动态载荷，会使得药柱在一定时间后发生疲劳裂纹或者破坏，影响甚至破坏药柱结构完整性。对于药柱结构的动力学分析，必须采用考虑动态效应的粘弹性本构模型，目前商业有限元软件只包含简单的本构模型，不能精确描述固体推进剂的动态力学响应。

工程经验表明，发动机药柱失效时推进剂往往并未达到其力学性能极限，而与“脱湿”损伤有很强相关性。目前“脱湿”等细观损伤对药柱结构完整性的影响。同时，温度是影响推进剂力学性能的一大因素。复合固体推进剂在动态载荷下呈现出非线性粘弹性特性，其力学响应呈现出较为复杂的率相关性。目前存在的推进剂率相关本构模型，还未有考虑“脱湿”等细观损伤和温度效应的。要实现动态载荷下推进剂药柱的精细结构完整性分析技术，必须有可以精确描述推进剂动态力学性能的本构模型。

现有的商业有限元软件中，对粘弹性材料结构进行动态分析时，采用简单的静态模型，所得结果与实际情况差别较大。但商业有限元软件提供了相应的软件接口，可实现对新型本构模型的编程实现。

发明内容

本发明要解决的技术问题：针对现有技术的上述问题，提供一种针对以固体推进剂为代表的粘弹性材料，在现有率相关本构模型的基础上，考虑脱湿细观损伤和温度对力学性能的影响，添加脱湿因子和温度系数，能更精确的预测粘弹性材料结构在冲击和振动等动态载荷下的力学行为，同时对构建的本构模型进行数值实现，可在有限元软件中对粘弹性材料结构进行准确的动力学分析的含脱湿的推进剂动态热粘弹性本构模型的构建和应用方法。

为了解决上述技术问题，本发明采用的技术方案为：

首先，本发明提供一种含脱湿的推进剂动态热粘弹性本构模型的构建方法，施步骤包括：

1)建立用于描述粘弹性材料动力学行为的率相关本构方程；

2)确立脱湿因子；

3)确立率相关本构方程的温度修正系数；

4)确立含脱湿因子、温度修正系数的率相关本构模型。

优选地，步骤1)中建立的率相关本构方程如式(1)所示；

式(1)中，σ为应力，E₀ε+αε²+βε³为非线性弹性项σ_E，ε为应变，E₀,α,β分别为非线性弹性项σ_E的弹性系数；为低频粘性项σ_VL，E₁,θ₁分别为低频粘性项σ_VL的弹性系数和松弛时间，为高频粘性项σ_VH，E₂,θ₂分别为高频粘性项σ_VH的弹性系数和松弛时间；t和τ均为时间，均为τ时刻的应变率。

优选地，步骤2)中确立脱湿因子如式(2)所示；

式(2)中，D为脱湿因子，ε为应变，ε_th为应变门槛值，为应变率，K_D,a,b均为试验测得的常数。

优选地，所述应变门槛值ε_th取值为推进剂材料的脱湿点伸长率ε_d。

优选地，步骤3)确立率相关本构方程的温度修正系数包括非线性弹性项σ_E的温度修正系数C_E、低频粘性项σ_VL的温度修正系数C_VL、高频粘性项σ_VH的温度修正系数C_VH，函数表达式分别如式(3)～(5)所示；

式(3)～(5)中，C_E为非线性弹性项σ_E的温度修正系数，C_VL为低频粘性项σ_VL的温度修正系数，C_VH为高频粘性项σ_VH的温度修正系数，T为当前温度，T_R为参考温度，T₀为初始温度，c_L,k_L,c_H,k_H均为试验测得参数。

优选地，步骤4)确立含脱湿因子、温度修正系数的率相关本构模型如式(6)所示；

式(6)中，σ_d,T为应力，T为当前温度，D为脱湿因子，C_E为非线性弹性项σ_E的温度修正系数，C_VL为低频粘性项σ_VL的温度修正系数，C_VH为高频粘性项σ_VH的温度修正系数，ε为应变，E₀,α,β分别为非线性弹性项的弹性系数；E₁,θ₁分别为应力表达式右侧第四项所示低频粘性项的弹性系数和松弛时间，E₂,θ₂分别为应力表达式右侧第五项所示高频粘性项的弹性系数和松弛时间；t和τ均为时间，均为τ时刻的应变率。

进一步地，本发明还提供一种含脱湿的推进剂动态热粘弹性本构模型的应用方法，实施步骤包括：

S1)采用前述含脱湿的推进剂动态热粘弹性本构模型的构建方法构建含脱湿因子、温度修正系数的率相关本构模型；

S2)推导含脱湿因子、温度修正系数的率相关本构模型的增量格式；

S3)使用有限元软件的扩展功能接口对增量格式的含脱湿因子、温度修正系数的率相关本构模型进行编程，得到用户子程序，使用有限元软件调用用户子程序。

优选地，步骤S2)推导含脱湿因子、温度修正系数的率相关本构模型的增量格式如式(7)所示；

式(7)中，为Δt时间内含脱湿和温度的Kirchhoff应力张量，为t+Δt时刻含脱湿和温度的Kirchhoff应力张量，为t时刻含脱湿和温度的Kirchhoff应力张量，D^t+Δt为t+Δt时刻的脱湿因子，ΔD^t+Δt为t+Δt时刻的脱湿因子的增量。

优选地，步骤S3)中使用有限元软件调用用户子程序的详细步骤包括：

S3.1)有限元软件在调用用户子程序之前，预先为用户子程序传递初始应力σ₀、总应变ε、应变增量Δε以及计算过程中的状态变量；

S3.2)在用户子程序中利用商业有限元软件主程序传递的指定参量获取矩阵，然后返回给主程序，主程序计算应力增量并且得到本次增量结束时的应力

S3.3)用户子程序更新本次增量结束时的状态变量以供下次调用用户子程序使用，所述状态变量包括应变张量和状态变量。

优选地，步骤S3.3)中用户子程序更新本次增量结束时的状态变量的函数式如式(8)所示；

式(8)中，为Δt时间内含脱湿和温度的Kirchhoff应力张量，为t+Δt时刻含脱湿和温度的Kirchhoff应力张量，为t时刻含脱湿和温度的Kirchhoff应力张量。

本发明含脱湿的推进剂动态热粘弹性本构模型的构建方法具有下述优点：本发明针对现有推进剂本构模型中所存在的问题，在现有率相关模型的基础上，考虑推进剂细观损伤和温度效应，添加脱湿因子和温度因子，构建含脱湿损伤和温度效应的率相关粘弹性本构模型(含脱湿因子、温度修正系数的率相关本构模型)，针对以固体推进剂为代表的粘弹性材料，在现有率相关本构模型的基础上，考虑脱湿细观损伤和温度对力学性能的影响，添加脱湿因子和温度系数，能更精确的预测粘弹性材料结构在冲击和振动等动态载荷下的力学行为。

本发明含脱湿的推进剂动态热粘弹性本构模型的应用方法：本发明含脱湿的推进剂动态热粘弹性本构模型的应用方法为采用本发明含脱湿的推进剂动态热粘弹性本构模型的构建方法构建含脱湿因子、温度修正系数的率相关本构模型，因此同样也具有本发明含脱湿的推进剂动态热粘弹性本构模型的构建方法的前述优点，本发明含脱湿的推进剂动态热粘弹性本构同时对构建的本构模型进行数值实现，可在有限元软件中对粘弹性材料结构进行准确的动力学分析，然后对三维本构模型进行推导，得到应力、脱湿因子等的增量格式，并在商业有限元软件中对构建的本构模型进行数值实现，为动态载荷下推进剂药柱的精确结构完整性分析提供理论模型和应用依据。

附图说明

图1为本发明实施例构建方法的原理示意图。

图2为本发明实施例建立的率相关本构方程模型。

图3为本发明实施例的推进剂“脱湿点”示意图。

图4为本发明实施例模型应用方法的原理示意图。

具体实施方式

如图1所示，本实施例含脱湿的推进剂动态热粘弹性本构模型的构建方法实施步骤包括：

1)建立用于描述粘弹性材料动力学行为的率相关本构方程；

2)确立脱湿因子；

3)确立率相关本构方程的温度修正系数；

4)确立含脱湿因子、温度修正系数的率相关本构模型。

参见图2，步骤1)中建立的率相关本构方程如式(1)所示；

式(1)中，σ为应力，E₀ε+αε²+βε³为非线性弹性项σ_E，ε为应变，E₀,α,β分别为非线性弹性项σ_E的弹性系数；为低频粘性项σ_VL，E₁,θ₁分别为低频粘性项σ_VL的弹性系数和松弛时间，为高频粘性项σ_VH，E₂,θ₂分别为高频粘性项σ_VH的弹性系数和松弛时间；t和τ均为时间，均为τ时刻的应变率。式(1)所示的率相关本构方程为朱-王-唐(ZWT)非线性粘弹性本构方程(简称ZWT模型)。ZWT模型是综合考虑高应变率和低应变率对本构响应的影响，采用两个松弛时间分别为θ₁和θ₂的Maxwell体和一个非线性弹簧并联而成，松弛时间为θ₁的Maxwell体描述低应变率下材料的力学行为，在高应变率条件下，因为实验观察时间比θ₁小很多(通常小3个数量级以上)，来不及松弛，表现为胡克弹簧特性；松弛时间为θ₂的Maxwell体描述高应变率下的粘弹响应，在低应变率条件下，因为实验观察时间比θ₂大很多(通常大3个数量级以上)，早己完全松弛，表现出牛顿粘壶特性，此时的应变率又远小于1S^-1，两者的乘积很小，Maxwell固体在低应变率时表现不出来，相当于不存在。ZWT模型可以用来分析受到从低应变率到高应变率冲击载荷时材料的力学响应。

在外力场作用下，随着推进剂粘合剂基体与固体填料颗粒的相对运动，基体和填料之间的化学键或者物理吸附被拉开，此时推进剂拉伸试件表面颜色开始泛白，这个过程称为“脱湿”。从应力-应变曲线来看，此时斜率发生较大变化，应力随应变增长效果减弱。如图3所示应力-应变曲线，对斜率迅速变化区域两旁曲线取切线，把所成钝角的平分线与曲线交点作为“脱湿点”，其应力为“脱湿点”强度，应变为脱湿点伸长率。本发明使用这种定量的方法，更准确强度地描述“脱湿”这种细观损伤。

根据王礼立等基于热激活机制的损伤演化模型，考虑到脱湿行为与应变关系的非线性，本实施例中，步骤2)中确立脱湿因子(debonding)如式(2)所示；

式(2)中，D为脱湿因子，ε为应变，ε_th为应变门槛值，ε为应变率，K_D,a,b均为试验测得的常数。本实施例中，所述应变门槛值ε_th取值为推进剂材料的脱湿点伸长率ε_d。即ε≥ε_d时，脱湿因子(debonding)如式(2-1)所示；

式(2-1)中，各参量与式(2)完全相同，故在此不再赘述。

多轴应力状态下，将多轴应力状态下的应变退化为等效应变求解脱湿因子项。假设材料为各向同性材料，且认为在三维应力状态下的脱湿损伤也是各向同性的，采用如下式(2-2)形式的等效应变，等效应变率如式(2-3)所示；

式(2-2)和式(2-3)中，ε_eq为等效应变，为等效应变率，ε_ij为Green应变张量的分量，为ε_ij对应的应变率。

本实施例中，步骤3)确立率相关本构方程的温度修正系数包括非线性弹性项σ_E的温度修正系数C_E、低频粘性项σ_VL的温度修正系数C_VL、高频粘性项σ_VH的温度修正系数C_VH，函数表达式分别如式(3)～(5)所示；

式(3)～(5)中，C_E为非线性弹性项σ_E的温度修正系数，C_VL为低频粘性项σ_VL的温度修正系数，C_VH为高频粘性项σ_VH的温度修正系数，T为当前温度，T_R为参考温度，T₀为初始温度，c_L,k_L,c_H,k_H均为试验测得参数。由于对率相关本构模型中粘性项，模量与频率相关，阿伦尼乌斯公式中频率因子与玻璃转换温度相关，可将低频粘性项σ_VL的温度修正系数C_VL设为如式(4)所示的函数表达式；同理，可将高频粘性项σ_VH的温度修正系数C_VH设为如式(5)所示的函数表达式。

本实施例中，步骤4)确立含脱湿因子、温度修正系数的率相关本构模型如式(6)所示；

式(6)中，σ_d,T为应力，T为当前温度，D为脱湿因子，C_E为非线性弹性项σ_E的温度修正系数，C_VL为低频粘性项σ_VL的温度修正系数，C_VH为高频粘性项σ_VH的温度修正系数，ε为应变，E₀,α,β分别为非线性弹性项的弹性系数；E₁,θ₁分别为应力表达式右侧第四项所示低频粘性项的弹性系数和松弛时间，E₂,θ₂分别为应力表达式右侧第五项所示高频粘性项的弹性系数和松弛时间；t和τ均为时间，均为τ时刻的应变率。

如图4所示，本实施例推进剂动态热粘弹性本构模型的应用方法的实施步骤包括：

S1)采用本实施例前述含脱湿的推进剂动态热粘弹性本构模型的构建方法构建含脱湿因子、温度修正系数的率相关本构模型；

S2)推导含脱湿因子、温度修正系数的率相关本构模型的增量格式；

S3)使用有限元软件的扩展功能接口对增量格式的含脱湿因子、温度修正系数的率相关本构模型进行编程，得到用户子程序，使用有限元软件调用用户子程序。

本实施例中，步骤S2)推导含脱湿因子、温度修正系数的率相关本构模型的增量格式如式(7)所示；

式(7)中，为Δt时间内含脱湿和温度的Kirchhoff应力张量，为t+Δt时刻含脱湿和温度的Kirchhoff应力张量，为t时刻含脱湿和温度的Kirchhoff应力张量，D^t+Δt为t+Δt时刻的脱湿因子，ΔD^t+Δt为t+Δt时刻的脱湿因子的增量。

本实施例中，步骤S2)推导含脱湿因子、温度修正系数的率相关本构模型的增量格式的详细步骤包括：

2.1)将步骤1)中建立的如式(1)所示的率相关本构方程由一维形式改写为式(7-1)所示的三维形式；

式(7-1)中，S_ij为第二类Piola-Kirchhoff应力张量的分量，E₀,α,β分别为非线性弹性项σ_E的弹性系数，[A]为转换矩阵，E_kl为原构型应变张量分量，E_ij为Green-Langrage应变张量的分量，E₁,θ₁分别为低频粘性项σ_VL的弹性系数和松弛时间，E₂,θ₂分别为高频粘性项σ_VH的弹性系数和松弛时间，t和τ均为时间；其中，转换矩阵[A]的函数表达式如式(7-2)所示；

式(7-2)中，[A]为转换矩阵，υ为泊松比，转换矩阵[A]的左下角与右上角为对称结构，故式(7-2)所示函数表达式省略留空。

2.2)根据式(7-1)，得到：非线性弹性项σ_E的增量格式如式(7-3)所示，低频粘性项σ_VL的应力张量VL_ij如式(7-4)所示，低频粘性项σ_VL的应力张量VL_ij的增量格式如式(7-5)所示，高应变率粘性项应力张量如式(7-6)所示，高应变率粘性项应力张量的增量格式如式(7-7)所示；

式(7-3)中，为t时刻的非线性弹性项应力张量的增量，ΔE_kl为原构型应变分量的增量，E_ij为Green-Langrage应变张量的分量，E₀,α,β分别为非线性弹性项σ_E的弹性系数，[A]为转换矩阵。

式(7-4)和式(7-5)中，为t时刻的低频粘性项σ_VL的应力张量，E₁,θ₁分别为低频粘性项σ_VL的弹性系数和松弛时间，[A]为转换矩阵，t和τ均为时间，E_kl为原构型应变张量分量，为低频粘性项应力张量在t时刻的增量，为在t+Δt时刻的低频粘性项应力张量的分量，为在t时刻的低频粘性项应力张量的分量，为在t时刻的原构型应变张量分量的增量，Δt为时间增量；

式(7-6)和式(7-7)中，为高应变率粘性项应力张量，E₂,θ₂分别为高频粘性项σ_VH的弹性系数和松弛时间，[A]为转换矩阵，t和τ均为时间，E_kl为原构型应变张量分量，为高频粘性项应力张量在t时刻的增量，为在t时刻的原构型应变张量分量的增量，Δt为时间增量；

2.3)综上式(7-1)～式(7-7)，可知Kirchhoff应力张量增量如式(7-8)所示，应力累加式的函数表达式如式(7-9)所示；

式(7-8)中，为应力张量分量的增量，E₀,α,β分别为非线性弹性项σ_E的弹性系数，[A]为转换矩阵，ΔE_kl为原构型应变分量的增量，E_kl为原构型应变分量，E_ij为Green-Langrage应变张量的分量，E₁,θ₁分别为低频粘性项σ_VL的弹性系数和松弛时间，E₂,θ₂分别为高频粘性项σ_VH的弹性系数和松弛时间，t和τ均为时间，为在t时刻的原构型应变张量分量的增量，Δt为时间增量，为t时刻的低频粘性项σ_VL的应力张量，为高应变率粘性项应力张量；

式(7-9)中，为t+Δt时刻的Kirchhoff应力张量，为t时刻的Kirchhoff应力张量，为Δt时间内的Kirchhoff应力张量；

2.4)对于如式(2-1)所示的脱湿因子(debonding)，t时刻等效应变如式(7-10)所示；由式(7-11)所示应力的增量形式可得式(7-12)，可知等效应变率如式(7-13)所示；在较短时间内，有式(7-14)成立，故可推导得出式(7-15)；

式(7-10)中，为t时刻等效应变，为t时刻的Green应变张量的分量；

式(7-11)～式(7-13)中，为t+Δt时刻的Green应变张量的分量，为t时刻的Green应变张量的分量，为t+Δt时刻的Green应变张量分量的增量，为t+Δt时刻的等效应变，为t+Δt时刻的Green应变张量分量，为t时刻的等效应变率，为t+Δt时刻的等效应变率，为t时刻的Green应变率张量分量，为t+Δt时刻的Green应变率张量分量；

2.5)当推进剂材料发生脱湿时，即ε≥ε_d时，脱湿因子的增量格式可以为为式(7-16)，并可表达为如式(7-17)所示；认为脱湿损伤是不可恢复的，则脱湿因子D的累加形式为式(7-18)所示，t时刻含脱湿应力张量如式(7-19)所示；

式(7-16)和式(7-17)中，ΔD^t为t时刻脱湿因子增量，为t+Δt时刻的等效应变率，ε_d为脱湿点伸长率，为t时刻的等效应变率，K_D,a,b均为试验测得的常数；

D^t+Δt＝D^t+ΔD^t(7-18)

式(7-18)和式(7-19)中，D^t+Δt为t+Δt时刻的脱湿因子，D^t为t时刻的脱湿因子，ΔD^t为t时刻的脱湿因子的增量，为t时刻含脱湿的Kirchhoff应力张量，为t时刻的Kirchhoff应力张量；

2.6)由式(7-19)及式(7-9)的应力累加式，则含脱湿应力张量的增量为式(7-20)；

式(7-20)中，为Δt时间内含脱湿的Kirchhoff应力张量；为t+Δt时刻含脱湿的Kirchhoff应力张量，为t时刻含脱湿的Kirchhoff应力张量，D^t+Δt为t+Δt时刻的脱湿因子，为t+Δt时刻的Kirchhoff应力张量，D^t为t时刻的脱湿因子，为t时刻的Kirchhoff应力张量，为Δt时间内的Kirchhoff应力张量，ΔD^t+Δt为t+Δt时刻的脱湿因子的增量；

2.7)在式(7-20)的基础上，针对温度效应，针对非线性弹性项σ_E添加温度修正系数C_E、为低频粘性项σ_VL添加温度修正系数C_VL、为高频粘性项σ_VH添加温度修正系数C_VH，最终可得到含脱湿因子、温度修正系数的率相关本构模型的增量格式如式(7)所示。

本实施例中，步骤S3)中使用有限元软件调用用户子程序的详细步骤包括：

S3.1)有限元软件在调用用户子程序之前，预先为用户子程序传递初始应力σ₀、总应变ε、应变增量Δε以及计算过程中的状态变量；

S3.2)在用户子程序中利用商业有限元软件主程序传递的指定参量获取矩阵，然后返回给主程序，主程序计算应力增量并且得到本次增量结束时的应力

S3.3)用户子程序更新本次增量结束时的状态变量以供下次调用用户子程序使用，所述状态变量包括应变张量和状态变量。

本实施例中，步骤S3.3)中用户子程序更新本次增量结束时的状态变量的函数式如式(8)所示；

式(8)中，为Δt时间内含脱湿和温度的Kirchhoff应力张量，为t+Δt时刻含脱湿和温度的Kirchhoff应力张量，为t时刻含脱湿和温度的Kirchhoff应力张量。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。