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

2019-07-31贺元吉赵宏伟江增荣余庆波

固体火箭技术 2019年3期

邓斌，贺元吉，赵宏伟，江增荣，余庆波

(1. 中国人民解放军96901部队，北京 100094；2. 北京理工大学，北京 100081)

0 引言

长期以来，因限于认知、测量手段等因素，通常为简化问题人们往往将火炸药、推进剂等粘弹性材料的泊松比视作常数处理[1-5]。此外，针对火炸药、推进剂等粘弹性材料的泊松比测量，目前国内现有行业标准如QJ 3228—2005[6]和GJB 770B—2005[7]，均将其泊松比视作常数进行处理，并采用了弹性泊松比试验测量手段。

大量研究表明[8-15]，实际的粘弹性材料泊松比为时间或频率的函数，且与温度紧密相关。采用泊松比作为计算参数进行结构分析时，泊松比的微小变化可导致计算结果的显著差异，尤其是对于火炸药、固体推进剂等近似不可压粘弹性材料，以往将其视作常数的处理方式，在某些情况下会导致计算结果的较大误差，尤其是对于某些特定条件如短时冲击载荷以及变温过程，或当粘弹性材料处于状态转变区时，粘弹性材料泊松比的变化往往较大[11-12]。此时用定泊松比模型已无法准确表征这种变化往往导致分析结果的较大误差[12]，难以准确反映粘弹性材料真实力学行为。随着人们对粘弹性结构分析精度要求的不断提高，与时间相关粘弹性泊松比的研究因而受到重视。目前，关于粘弹性泊松比理论和试验研究的成果也较多[9-16]，粘弹性问题计算方法研究也主要集中于定泊松比模型[17-18]，而关于变化泊松比粘弹性结构数值分析方法及其应用研究则鲜见报道。

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

1 粘弹性本构模型

1.1 粘弹性时间相关泊松比形式

(1)

一般地，纵向应变并非阶跃形式，而是随时间的连续变化量。若在任意τ时刻施加一纵向阶跃微应变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)

其中

(6)

(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)

(13)

Δ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)

(22)

(23)

其中

(24)

(25)

(26)

应用中值定理，对式(26)进行积分计算，可得

(27)

其中

(28)

(29)

根据式(8)和式(9)，式(18)和式(20)可写成如下形式：

(30)

(31)

(32)

(33)

(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)

(48)

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级数系数

表2 泊松比Prony级数系数

表3 其他材料参数

3.1 验证算例：粘弹性薄板的应力松弛分析

采用变泊松比线性粘弹性本构方程进行分析。为模拟单轴应力松弛试验过程，采用尺寸为10 mm×50 mm×5 mm的方形薄板模型，其有限元网格模型见图1。

图1 单向拉伸有限元模型

(53)

(54)

(55)

将式(54)两边对t求导并结合式(55)，可得其横向应εy(t)变表达式：

(56)

图2 纵向应力沿随时间变化对比曲线

图3 横向应变沿随时间变化对比曲线

从图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 药柱有限元模型

(a)径向应力 (b)环向应力

(a)径向应变 (b)环向应变

由图5和图6可知，在加载过程的任意时刻，变泊松比所对应的应力应变结果，均位于定泊松比取0.487和0.491时所对应结果之间，且在开始阶段，变泊松比所对应结果接近于定泊松比为0.487时的对应结果，而随着时间的增加，逐渐向定泊松比取0.491时的对应结果靠近，在加载至0.3 s时，其各部位的应力应变结果均已非常接近于定泊松比取0.491时的对应结果(图7和图8)，这与变泊松比值从0.487变化到0.491趋势是一致的，符合力学变化规律。

(a)径向应力 (b)环向应力

(a)径向应变 (b)环向应变

综上，随着分析时间的加长，变泊松比值逐渐趋于接近平衡值，可以预测，对于长时间持续加载分析过程，采用取平衡泊松比值的定泊松比参数与变泊松比参数所得对应结果间几乎无差异，因而此时，在缺乏变泊松比参数时，可采用取平衡泊松比值的定泊松比粘弹性本构模型进行分析。然而，对于泊松比值变化较大的过程，如短时加载、剧烈变温过程，或当粘弹性材料处于状态转变区时，采用泊松比取某一定值所得结果与实际差异会较大，有必要考虑变泊松比粘弹性本构模型进行分析。