含能材料变泊松比粘弹性本构及其数值实现方法

含能材料变泊松比粘弹性本构及其数值实现方法DENG Bin;HE Yuanji;ZHAO Hongwei;JIANG Zengrong;YU Qingbo
【摘要】为研究含能材料变泊松比粘弹力学行为,基于热流变简单材料假设,建立了一种采用时间相关泊松比的三维热粘弹性本构方程;利用积分算法对该本构进行数值离散,给出了相应的应力更新方法;基于Abaqus软件UMAT接口技术编制了相应材料子程序,完成了算例分析.结果表明,数值解与对应解析解吻合良好,并通过了合理性算例验证,分析方法及算法程序有效,可为后续三维装药精细结构完整性分析提供支持;粘弹材料泊松比具有时间相关性,其较小变化对计算结果影响显著,尤其是在短时加载或状态转变过程.
【期刊名称】《固体火箭技术》
【年(卷),期】2019(042)003
【总页数】7页(P290-296)
【关键词】含能材料;时间相关泊松比;有限元分析;材料子程序
【作者】DENG Bin;HE Yuanji;ZHAO Hongwei;JIANG Zengrong;YU Qingbo 【作者单位】;;;;
【正文语种】中文
【中图分类】V512
0 引言
长期以来，因限于认知、测量手段等因素，通常为简化问题人们往往将火炸药、推
进剂等粘弹性材料的泊松比视作常数处理[1-5]。

此外，针对火炸药、推进剂等粘
弹性材料的泊松比测量，目前国内现有行业标准如QJ 3228—2005[6]和GJB
770B—2005[7]，均将其泊松比视作常数进行处理，并采用了弹性泊松比试验测
量手段。

大量研究表明[8-15]，实际的粘弹性材料泊松比为时间或频率的函数，且与温度紧密相关。

采用泊松比作为计算参数进行结构分析时，泊松比的微小变化可导致计算结果的显著差异，尤其是对于火炸药、固体推进剂等近似不可压粘弹性材料，以往将其视作常数的处理方式，在某些情况下会导致计算结果的较大误差，尤其是对于某些特定条件如短时冲击载荷以及变温过程，或当粘弹性材料处于状态转变区时，粘弹性材料泊松比的变化往往较大[11-12]。

此时用定泊松比模型已无法准确表征
这种变化往往导致分析结果的较大误差[12]，难以准确反映粘弹性材料真实力学行为。

随着人们对粘弹性结构分析精度要求的不断提高，与时间相关粘弹性泊松比的研究因而受到重视。

目前，关于粘弹性泊松比理论和试验研究的成果也较多[9-16]，粘弹性问题计算方法研究也主要集中于定泊松比模型[17-18]，而关于变化泊松比
粘弹性结构数值分析方法及其应用研究则鲜见报道。

为此，考虑到粘弹性材料泊松比的时间相关性以及泊松比参数变化对计算结果的显著影响，本文采用一种时间相关粘弹性材料泊松比形式建立三维线性粘弹性本构模型，研究相应工程应用的数值分析方法，旨在为火炸药、推进剂等含能材料装药精细结构完整性分析提供支持。

1 粘弹性本构模型
1.1 粘弹性时间相关泊松比形式
对于各向同性线粘弹性材料，泊松比定义为在给定阶跃纵向应变条件下，H(t)为Heaviside函数，与时间相关的横向应变εy(t)和纵向应变的负比值[8]，即
一般地，纵向应变并非阶跃形式，而是随时间的连续变化量。

若在任意τ时刻施
加一纵向阶跃微应变dεx(τ)，则持续作用至t时刻的横向应变响应为dεy(t)=-ν(t-τ)dεx(τ)，当τ在整个[0,t]上连续变化时，根据线性系统的Boltzmann叠加原理，可得t时刻的总横向应变响应：
(2)
可知，横向应变εy(t)要滞后于纵向应变εx(t)响应历程，且粘弹性泊松比是横向变形的记忆函数，而非简单的代数关系。

为区分传统的定泊松比，后续也称时间相
关泊松比ν(t)为变泊松比。

1.2 采用时间相关泊松比的粘弹性本构方程
对于各向同性线弹性材料，采用杨氏模量E和泊松比ν表示的本构形式如下：
(3)
其中
(4)
根据弹性-粘弹性对应原理，得到式(3)、式(4)在相域内的对应线粘弹性关系式，再对其作Laplace逆变换，并基于热流变简单材料假设，可得时域内的线热粘弹性
本构方程：
(5)
其中
(7)
其中，泊松比ν(t)和松弛模量E(t)Prony级数形式分别为
(8)
(9)
其中，为温差；ξ和ξ′表达式为
(10)
其中，aT为温度平移因子，满足如下WLF方程：
(11)
式中 C1、C2为材料常数；T为当前温度；T0为参考温度。

2 本构方程的数值模型
2.1 本构方程的增量形式
针对式(5)所示的本构方程，采用增量有限元法对其进行数值离散，给出其增量形式。

将分析时间[0,t]划分为[0,t1],[t1,t2], …,[tm,tm+1],…,[tM-1,tM]共M个子时间增量步，则可将其在任意增量步[tm,tm+1]内的增量形式写成：
(12)
其中
Δσij(tm+1)=σij(tm+1)-σij(tm)
ΔSij(tm+1)=Sij(tm+1)-Sij(tm)
(14)
Δσkk(tm+1)=σkk(tm+1)-σkk(tm)
(15)
针对式(6)所示本构方程偏张量部分，结合Stieltjes卷积定理，可得其在[tm,tm+1]内的增量形式：
(16)
其中
(17)
(18)
(19)
(20)
(21)
综合式(8)、式(17)，可得的形式如下：
(22)
假设各增量步内应力随ε呈线性变化，可得式(22)中表达式：
(23)
其中
(24)
(25)
将式(9)代入式(19)可得表达式：
(26)
应用中值定理，对式(26)进行积分计算，可得
(27)
其中
(28)
(29)
根据式(8)和式(9)，式(18)和式(20)可写成如下形式：
(30)
(31)
其中，和求解递推公式：
(32)
(33)
根据式(9)和式(21)，易得的计算式：
(34)
对于式(7)所示本构部分，类似式(6)处理方法，可得其在[tm,tm+1]内增量形式：
(35)
其中
(36)
(37)
(38)
(39)
(40)
针对式(36)～式(40)，可采用类似式(17)～式(21)的方法进行数值离散，得到式
(36) ～式(40)的表达式：
(41)
(42)
(43)
(44)
(45)
其中
(46)
(47)
针对式(25)中和式(33)中考虑计算精度和稳定性，采用如下通用计算式βn(Δξ)，即
(48)
若将式(48)中的τn先后替换成则分别得到的计算式。

2.2 本构方程的切线刚度
利用一致切线刚度可保证有限元分析所采用的Newton法具有二阶收敛速率。

根据切线刚度的定义，对于小变形或者小体变的大变形问题，其一致切线刚度的形式为
(49)
联立式(12)和式(49)，可得在第tm+1时刻的切线刚度张量各个分量：
(50)
(51)
(52)
3 算例分析
针对上述变泊松比粘弹性本构增量模型，基于Abaqus软件平台二次开发技术，采用Fortran语言编写相应UMAT材料子程序，并依托软件平台前后处理和求解功能，实现该本构的工程应用。

以某HTPB推进剂为例，取式(11)中的C1=4.97，C2=156.1，T0=293.15 K；式(9)所示松弛模量取4项时的各参数见表1；根据某型推进剂泊松比试验数据，拟合得到式(8)所示取5项时的相关系数见表2；其他部件材料性能参数见表3。

表1 松弛模量Prony级数系数Table 1 Prony series coefficients for relaxation modulusn01234τn/s—27.3222.732 20.273 220.027 322En/MPa15.8681.154 84.058 61.532 06.469 7
表2 泊松比Prony级数系数Table 2 Prony series coefficients of viscoelastic Poisson's ration012τνn/s—30.130 7301.307νn0.486 90.001 2740.011
239n345τνn/s3013.0730 130.7301 307νn-0.012 240.001 050.011 486
表3 其他材料参数Table 3 Other material parameters材料参数构件壳体绝热层推进剂E/MPa2.0×10512.19E(t)αT/(1/K)1.1×10-59.2×10-51.0×10-
4ν0.30.495ν or ν(t)ρ/(kg/m3)780012001860
3.1 验证算例：粘弹性薄板的应力松弛分析
采用变泊松比线性粘弹性本构方程进行分析。

为模拟单轴应力松弛试验过程，采用尺寸为10 mm×50 mm×5 mm的方形薄板模型，其有限元网格模型见图1。

图1 单向拉伸有限元模型Fig.1 Finite element model of thin-slab structure 对于上述粘弹性薄板模型，分析时约束方板左侧和下侧表面的法向位移，保持温度恒定，并在模型的右侧表面施加纵向阶跃位移则对应的应变响应为此时的横向应力σy(t)为0。

此问题已退化为平面应力问题，此时有σz=0，σxz=0，σyz=0，再根据式(5)可得
(53)
(54)
根据初始应变条件并结合式(53)，可得其纵向应力σx(t)表达式:
(55)
将式(54)两边对t求导并结合式(55)，可得其横向应εy(t)变表达式：
(56)
针对上述问题的分析，取计算时间t=100 s，取并分别采用本文有限元解和相应解析解进行求解。

得到纵向应力、横向应变随时间的变化曲线，分别见图2和图3。

图2 纵向应力沿随时间变化对比曲线Fig.2 Comparison of longitudinal stress vs time
图3 横向应变沿随时间变化对比曲线Fig.3 Comparison of transverse strain vs time
从图2和图3中可看出，采用变泊松比粘弹性本构模型进行分析时，本文有限元解与对应解析解吻合良好，且具有很高的计算精度，这进一步验证了本文算法和程序的有效性。

3.2 合理性算例：内压载荷下的装药结构分析
相比于定泊松比本构模型，同等条件下的变泊松比粘弹性本构模型相关问题的求解要复杂得多，一般难以给出其解析解，且因限于条件暂无法开展相关试验验证。

下面就算法程序做一些合理性验证分析，其基本思路：将变泊松比粘弹性本构模型所得结果与相应的定泊松比本构模型所得结果进行比较，并根据力学规律来判断结果的合理性来验证算法程序合理性。

根据式(8)及表2中的系数可知，该推进剂泊松比为时间递增函数，其初始值约为0.487，在0.3 s内其值位于0.487～0.491之间，若本文算法程序有效，则对于变泊松比模型所得结果，应位于定泊松比为0.487和0.491时所对应结果之间，且开始时刻结果接近于0.487对应值，而在0.3 s时结果接近于0.491对应值。

以受均匀内压载荷下的细长固体火箭发动机为例，模型的几何尺寸为内径200 mm，外径397 mm，取其一段所建计算模型见图4。

在同等条件下进行分析，当本构(5)退化至定泊松比线粘弹性本构时，所得结果与Marc线粘弹性本构解吻合良好，限于篇幅，此处不详述。

在其内表面施加均布压力载荷P(t)=6(1-e-20t)MPa，计算时间t=0.3 s，分析时取定泊松比ν=0.487和0.491分别进行计算，采用有限元网格模型、载荷及位移边界条件，先后考虑本文变泊松比和定泊松比(取为0.487和0.491)粘弹性本构模型，分别对上述问题进行
分析，其中该问题在采用定泊松比时所得环向、径向的应力应变关系式见文献[19]。

得到圆管内表面的应力、应变随时间变化曲线分别见图5和图6所示；加载至0.3 s时的应力、应变沿径向变化规律分别如图7和图8所示。

图4 药柱有限元模型Fig.4 Finite element model of grain
(a)径向应力 (b)环向应力图5 应力随时间变化曲线Fig.5 Curves of stress vs time
(a)径向应变 (b)环向应变图6 应变随时间变化曲线Fig.6 Curves of strain vs time
由图5和图6可知，在加载过程的任意时刻，变泊松比所对应的应力应变结果，
均位于定泊松比取0.487和0.491时所对应结果之间，且在开始阶段，变泊松比
所对应结果接近于定泊松比为0.487时的对应结果，而随着时间的增加，逐渐向
定泊松比取0.491时的对应结果靠近，在加载至0.3 s时，其各部位的应力应变结果均已非常接近于定泊松比取0.491时的对应结果(图7和图8)，这与变泊松比值从0.487变化到0.491趋势是一致的，符合力学变化规律。

(a)径向应力 (b)环向应力图7 应力沿径向变化曲线Fig.7 Curves of stress vs radial path
(a)径向应变 (b)环向应变图8 应变沿径向变化曲线Fig.8 Curves of strain along radial path
综上，随着分析时间的加长，变泊松比值逐渐趋于接近平衡值，可以预测，对于长时间持续加载分析过程，采用取平衡泊松比值的定泊松比参数与变泊松比参数所得对应结果间几乎无差异，因而此时，在缺乏变泊松比参数时，可采用取平衡泊松比值的定泊松比粘弹性本构模型进行分析。

然而，对于泊松比值变化较大的过程，如短时加载、剧烈变温过程，或当粘弹性材料处于状态转变区时，采用泊松比取某一定值所得结果与实际差异会较大，有必要考虑变泊松比粘弹性本构模型进行分析。

4 结论
(1)采用变泊松比粘弹性本构的本文数值解与对照解析解吻合良好，合理性算例结果符合力学规律，验证了本文算法程序的正确性。

(2)此粘弹材料泊松比为时间递增函数，在任意给定时间内存在相应上下限值，且该变泊松比模型所得结果，均位于定泊松比模型分别取该上下限值时所对应结果之间，且其变化趋势符合规律，验证了本文算法程序的合理性。

(3)利用本文算法及材料子程序方法，有效实现了时间相关泊松比粘弹性本构方程的数值分析与工程应用，可为后续火含能材料装药精细结构分析提供支持。

参考文献：
【相关文献】
[1] Park S W,Schapery R A.A viscoelastic constitutive model for particulate composites with growing damage[J].International Journal of Solids Structures,1997,34(8):931-947. [2] Jung G D,Youn S K,Kim B K.A three-dimensional nonlinear viscoelastic constitutive model of solid propellant[J].International Journal of Solids and
Structures,2000,37(34):4715-4732.
[3] Ho S Y.High strain-rate constitutive models for solid rocket propellants[J].Journal of Propulsion and Power,2002,18(5):1106-1111.
[4] Hinterhoelzl R M,Schapery R A.FEM implementation of a three-dimensional viscoelastic constitutive model for particulate composites with damage
growth[J].Mechanics of Time-Dependent Materials,2004,8(1):65-94.
[5] Xu F,Aravas N,Sofronis P.Constitutive modeling of solid propellant materials with evolving microstructural damage[J].Journal of the Mechanics and Physics of
Solids,2008,56(5):2050-2073.
[6] 中国航天标准化研究所.复合固体推进剂泊松比试验方法:QJ 3228—2005 [S].2006.
China Institute of Aerospace Standardization.Poisson's ratio test method for composite solid propellants:QJ 3228—2005[S].2006.
[7] 国防科学技术委员会.火药试验方法.方法414.4 泊松比引伸计法:GJB 770B—2005 [S].2006. National Defense Science and Technology Commission.Test methods for
explosives.Method 414.4:GJB 770B—2005 [S].2006.
[8] Christensen R M.Theory of viscoelasticity:An introduction[M].New York:Academia Press,1982.
[9] 赵伯华.粘弹性泊松比与动态复数泊松比的研究 [J].推进技术,1995,16(3):1-7.
ZHAO Bohua.An investigation on viscoelastic Poisson's ratio and dynamic complex Poisson's ratio[J].Journal of Propulsion Technology,1995,16(3):1-7.
[10] Farruggia F,Perre P.Microscopic tensile test in the transverse plane of earlywood and latewood parts of spruce[J].Wood Science and Technology,2000,34(2):65-82.
[11] Qi H Y,Zhou L Z,Yang X G.Measurement of Young's modulus and Poisson's ratio of thermal barrier coatings[J].Chinese Journal of Aeronautics,2005,18(2):180-184.
[12] Lakes R S.The time-dependent Poisson's ratio of viscoelastic materials can increase or decrease[J].Cellular Polymers,1992,11(1):466-469.
[13] Tschoegl N W,Knauss W G,Emri I.Poisson's ratio in viscoelasticity:A critical
review[J].Mechanics of time-dependent materials,2002,6(1):3-51.
[14] Michaeli M,Shtark A,Grosbein H,et al.Analytical,experimental and computational viscoelastic material characterizations absent Poisson's ratios[R].AIAA 2011-1809. [15] Hilton H H.Clarifications of certain ambiguities and failings of Poisson's ratios in linear viscoelasticity[J].Journal of Elasticity,2011,104(1):5155-5161.
[16] 张海联,周建平.固体推进剂药柱的近似不可压缩粘弹性增量有限元法 [J].固体火箭技
术,2001,24(2):36-40.
ZHANG Hailian,ZHOU Jianping.Near-incompressible viscoelastic incremental finite element analysis of solid propellant grains [J].Journal of Solid Rocket
Technology,2001,24(2):36-40.
[17] 孟红磊,鞠玉涛.含损伤非线性粘弹性本构模型及其数值仿真应用[J].固体火箭技
术,2012,35(6):764-768.
MENG Honglei,JU Yutao.Nonlinear viscoelastic equation with cumulative damage and its application on numerical simulation [J].Journal of Solid Rocket
Technology,2012,35(6):764-768.
[18] Zienkiewicz O C,Taylor R L.The finite element method for solid and structural mechanics (Sixth edition)[M].Oxford:Elsevier Butterworth-Heinemann publications,2005.
[19] 王元有.推进剂药柱在内压力载荷下的应力、应变的粘弹性分析 [J].兵工学报,1983,4(3):20-33. WANG Yuanyou.Viscoelastic analysis of stresses and strains in solid propellant grains under pressure loadings[J].Acta Armamentarii,1983,4(3):20-33.。