摘要
为分析四组元HTPB推进剂不同拉伸速率组合对并行流变框架模型的精度影响,基于并行流变框架建立了推进剂的非线性黏弹性本构模型;通过对不同拉伸速率实验进行组合构建了本构模型,得到相应的模型参数。将有限元模型和数值计算结果与实验结果进行对比,将不同速率组别构建的本构模型误差进行了对比分析。结果表明当采用高速率与低速率组别相结合的方式即可精确建立本构模型,无须进行中间速率的大量实验。且高速率推进剂拉伸实验的速率可以截至3000 mm/min,无须再增加速率。此分析为简化推进剂材料实验提供了合理建议,提高了实验效率。同时,为快速预测推进剂材料的力学性能开辟了一条有效途径。
Abstract
In order to analyze the effect of different combinations of stretching rates of four-component HTPB propellant on the accuracy of the parallel rheological framework model, a nonlinear viscoelastic constitutive model of the propellants was developed based on the parallel rheological framework method. A constitutive model was constructed by combining experiments with different tensile rates, the corresponding constitutive model parameters were obtained and compared with experimental results using finite element models and numerical calculations. The constitutive model errors of the calibration for different rate groups were compared and analyzed. The results show that the model can be calibrated more accurately when a combination of high rate and low rate groups is used, without the need for extensive experiments at intermediate rates. And the high rate propellant tensile test rate can be up to 3000 mm/min, no need to increase the rate. This analysis provides justified suggestions for simplifying propellant material tests and improving the efficiency of the experiments. Concurrently, it paves an effective way to rapidly predict the mechanical properties of propellant materials.
固体火箭发动机由于其存储时间较长、机动性强和操作简便等特点,被广泛用于航天领域特别是战术导弹方面。其中固体火箭推进剂作为发动机的主要能量来源,在保障导弹完成技战术指标方面起着至关重要的作用。而推进剂燃烧前,在极短的时间内需达到稳定的压力状态,整个过程中压力梯度较大,高燃气体易对推进剂材料形成较强的冲击效应[1]。张明等[2]在对点火压强和增压速率的研究中发现,点火瞬态应变响应结果明显高于准静态计算结果,也从一个侧面表现出增压速率对于推进剂的力学响应有着至关重要的影响。而在推进剂点火建压的过程中,在围压环境下推进剂周向出现较大的拉应变是推进剂失效的主要原因。所以如何设计推进剂的力学性能实验和利用实验结果建立准确的本构模型是当前研究的热点。
在设计实验时,目前主要考虑的参数包括温度、压力和反应速率。王玉峰等[3]对端羟基聚丁二烯(hydroxyl-terminated polybutadiene,HTPB)推进剂进行了不同应变率条件下的实验,发现随着拉伸速率的提高,推进剂的抗拉强度和伸长率呈增加的趋势。申志彬等[4]针对三组元HTPB推进剂在宽温和变压力条件下进行了实验研究,得到了推进剂的相关力学性能结果。马浩等[5]对HTPB推进剂也进行了不同温度下的快慢拉组合的实验研究,发现拉伸速率和温度对于推进剂的“脱湿点”有较大的影响。通过以上的研究可见,不同的实验设计目的都在于对推进剂在不同条件下的力学性能进行反映,明晰推进剂材料内部的变化机理,从而为更精确地建立非线性本构模型打下基础。
关于准确地建立本构模型方面,推进剂是一种典型的黏弹性材料,不同于其他的金属材料,黏弹性材料的力学性能随温度、时间、速率的变化很大。经典的线弹性本构模型虽然可以针对黏弹性材料的蠕变和松弛现象给出一定的描述,但是与实验数据对比发现准确度并不是很高。为了更加精确地描述推进剂的黏弹性行为,大量的专家学者对非线性本构模型进行了研究[6-7]。Henriques等[8]用分数阶导数本构模型来描述两种聚合物泡沫的黏弹性行为。Zhang等[9]通过分数阶黏弹性对材料的大变形进行描述。现阶段使用最多的非线性本构模型是朱-王-唐(zhu-wang-tang,ZWT)本构模型。张亮[10]对HTPB推进剂的围压性能进行了研究,在ZWT本构模型中加入损伤函数和围压因子,发现模型可以很好地预示推进剂的脱湿现象和围压现象,并且对推进剂曲线出现双峰现象进行了解释。周风华等[11]在对有机玻璃的研究中基于ZWT本构模型加入损伤项,对高应变率下的有机玻璃非线性力学行为进行了描述。ZWT本构模型中拟合的方式较为复杂,参数较多,并且适用应变率范围较窄的情况,故而建立宽应变率范围内精确的本构模型是黏弹性材料研究的重点。而并行流变框架模型结合了橡胶领域中常用的超弹性本构模型与流变学中的黏弹性本构模型,其具体形式将在1.2节中做详细的介绍。其思想最早由Bergström等[12]通过应力松弛实验发现,实验中发现不同拉伸和压缩状态下的橡胶材料在卸载后都有着向一种平衡状态靠近的行为,故而将橡胶的这种非线性力学行为表达为超弹性和黏弹性项之和。其中超弹性模型为八链模型,黏弹性模型是以朗之万统计函数模型为基础得到微观分子与宏观流变特性的模型,研究发现此模型的适应性广,而且随着应变率的变化,模型对于黏弹性材料的力学特征有很好的捕捉效果。赵华等[13]用Prony级数表示的线性黏弹性材料模型和用五项Mooney-Rivlin应变能函数表示的超弹性材料模型相结合建立聚氨酯弹性体的非线性黏弹本构模型,数值分析结果表明该模型能准确地表征聚氨酯弹性体材料在压缩变形时的响应特性。非线性黏弹性部分描述橡胶材料在振动、冲击载荷下的动态响应,模型预测结果与实验结果吻合度较高。周梦雨等[14]通过对橡胶材料运用三阶Ogden超弹性和幂率硬化流动法则相结合的并行流变框架模型,与传统的显示动力学对比分析了轮胎在滚动过程中的温升变化,结果表明并行流变框架模型更能体现材料的黏弹性特征响应。
本文借助前期研究的结果,基于并行流变框架模型对四组元HTPB推进剂材料展开研究。重点设计不同拉伸速率组别实验,分析不同拉伸速率组合对本构模型精度的影响。提出合理的实验建议,为后续简化推进剂实验数据和快速预测推进剂力学性能提供方法,为分析固体火箭推进剂结构完整性奠定基础。
1 本构模型与实验
HTPB推进剂材料的典型黏弹性特性导致传统线弹性本构模型对其力学性能的描述能力较弱。为了构建能够表征其非线性力学特性的模型,本文借用橡胶工业中常使用的并行流变框架模型将超弹性和黏弹性模型并联,用以表征其力学特性。
1.1 实验设计
1.1.1 实验开展
为了减少实验次数并且获得相对准确的本构模型,达到对HTPB推进剂材料非线性力学行为的准确模拟,实验共采用 6 种不同速率(分别为10 mm/min、100 mm/min、500 mm/min、2 000 mm/min、3 000 mm/min和4 000 mm/min)对HTPB推进剂试件进行单轴拉伸实验,试件采用GJB 770B—2005《火药试验方法》中制备的B型试样,试样尺寸:长度为(120±0.5)mm,工程标距为(70±0.5)mm,厚度为(10±0.5)mm,宽度为(25±0.5)mm,具体如图1所示。
图1试件示意图
Fig.1Schematic diagram of test piece
实验主要通过电子万能试验机进行,用与试件相一致的夹具对试件进行固定。而后采取引伸计测量试件的标距部分(70 mm)的应变变化。实验的加载速率主要通过电子万能试验机进行控制,实验每0.5 s记录一次数据,主要包括加载的时间、加载的速率、夹具下端的负载拉力和位移。
通过对所测数据进行工程应力应变向真实应力应变的转化可以得到6种不同拉伸速率下的HTPB推进剂材料力学性能曲线,具体如图2所示。
图26种拉伸速率下推进剂真实应力应变曲线
Fig.2True stress-strain curves of propellant at six different tensile rates
1.1.2 实验组别设计
实验测试速率共有10 mm/min、100 mm/min、500 mm/min、2 000 mm/min、3 000 mm/min、4 000 mm/min六种不同的数据。由于2 000 mm/min的数据需要用来对模型的准确性进行验证,所以将其从校核参数的数据中提出,以剩下的5组速率校核模型。5组速率共有120种不同的组合方式,但是根据速率分布的规律,采用速率从高到低依次递减选择了8个组别,具体见表1。
1)第1~4组的设置是为了探究减少低速率组别对本构模型精度的影响。
2)第5组的设置是为了探究由最高和最低速率构建出的本构模型能否达到所需求的精度要求,为了验证猜测——低速率对于反应推进剂缓慢变化过程(蠕变或松弛)有较强的作用,不能轻易将其舍去。故设置了最高速率与最低速率组合(第5组)和次高速率与最低速率组合(第6组)来验证此种猜想,舍弃了中间速率100 mm/min、500 mm/min、2 000 mm/min。
3)第6组的设置是为了与第5组形成对照实验,验证3 000 mm/min与4 000 mm/min速率所构建本构模型的差异。此中低速率不取100 mm/min而取10 mm/min,一方面是因为相比100 mm/min单轴拉伸实验,10 mm/min单轴拉伸实验更加容易实现,方便以后实验的开展;另一方面本文认为推进剂在长期贮藏时蠕变与松弛是很典型的变形方式,而蠕变与松弛的速率更接近10 mm/min,所以低速率组别选取10 mm/min更符合研究目的。
4)第7组、第8组与第5组、第6组形成对照组别,验证在高、低速率中加入中间速率是否有利于提高模型的精度。
表1拉伸速率组别设置
Tab.1 Tensile rates group settings
1.2 基于并行流变框架非线性本构模型
并行流变框架模型是广义Maxwell模型的变形,广义Maxwell模型如图3所示。不同于广义Maxwell模型的是并行流变框架中的弹性元件与黏性元件是非线性的。并行流变框架模型如图4所示,弹性元件(图4中A、B中的非线性弹簧)作为主承载元件,黏性元件(图4中B中的非线性黏壶)捕捉材料的黏性流动但不承载。将两者结合实现对材料的非线性力学行为准确的描述。
1.2.1 超弹性本构模型
超弹性模型的本构关系一般使用应变能密度函数对变形不变量求导表述。本文中弹性元件(图4的A、B网络中的非线性弹性元件)选用Yeoh模型,其应变能密度函数为
图3广义Maxwell模型
Fig.3Generalized Maxwell model
图4并行流变框架模型
Fig.4Parallel rheological framework model
(1)
式中:C10、C20、C30是材料的相关参数;κ为材料的体模量,需要通过实验测定。
(2)
其中,F是变形梯度,I1可以表示为
(3)
C是右Cauchy-Green变形张量,λi是材料变形主方向的变形量。
将式(2)和式(3)代入式(1)中,并使式(1)对I1求导可以得出应力的表达式为
(4)
其中,b*是左Cauchy-Green变形张量的应变偏量矩阵,I是二阶单位张量。下面考虑一种特殊形式:单轴拉伸。假设推进剂为近似不可压缩材料,单轴拉伸的变形梯度矩阵可以写为
(5)
1.2.2 黏弹性本构模型
本文中黏性元件(图4的B网络中的非线性黏性元件)使用Boyce流动法则,其黏性流动遵循下述表达式:
(6)
其中,γB为材料的黏性应变,
(7)
τB是B网络的有效应力测度,表示为
(8)
ξ、C、和m都是材料相关的参数,也是后文中待校准的参数。C的物理意义为材料变形快慢,主要对本构模型中变形梯度进行修正;是黏性流动的驱动应力,即超过这个应力,材料出现黏性流动;m主要反映材料流动速度的快慢;因为ξ是一个很小的参数,只是为了防止表达式出现数值奇异,所以一般ξ=0.001。
黏性变形率可以表示为
(9)
联立式(6)~(9)可以获得黏性变形,从而可以获得弹性变形。则并行流变框架的本构模型可以表达为
(10)
(11)
(12)
其中,s是B网络的应力与A网络的应力的刚度系数之比。非线性弹簧使用超弹性本构模型,非线性黏壶采用黏弹性本构模型。如上文所述非线性黏壶不承载,故而产生应力的只有A、B中的非线性弹簧,其表达式如式(11)~(12)所示,而并行流变框架的总应力表达式如式(10)所示,简单来说即对并行流变框架中的A、B两网络中的承载部分进行叠加。A、B两网络都采用非线性超弹性本构模型,B网络中更添加了非线性的流动模型,使得总体模型更能精确地捕捉材料的非线性力学性能变化。
2 数值计算与分析
为了提高计算的效率并且达到所要求的精度,基于上述并行流变框架理论和不同速率下的实验对推进剂的本构模型进行拟合,通过分析得到相关的结论。
2.1 参数的确定
描述一个2节点并行流变框架非线性本构模型需要7个参数,具体如表2所示。
表2并行流变框架参数
Tab.2 Parameters of parallel rheological framework
上述参数的初值可以根据文献[15]获得,具体如表3所示。
表3优化参数初值
Tab.3 Initial parameter values for optimization
校核的具体流程如图5所示。
拟合参数的评判标准为实验值与数值计算出的值方差之和达到最小,如式(13)所示。
(13)
其中:N为每种速率实验中的取点个数;n为所取的速率组别个数,即表1中速率的个数,在拟合中考虑了不同速率对本构模型参数的影响。
图5本构模型构建流程
Fig.5Construction process of the constitutive model
至此,不同种拉伸速率组合构建出的本构模型可以用以上的7个参数表示出来。
2.2 不同速率的计算结果
针对1.1.2节中的不同的8组速率实验,运用粒子种群算法对本构模型进行优化和参数的获取。粒子种群算法不容易陷入局部解,善于找寻整体最优解。经过计算,每组速率校核出的参数如表4所示。
表4不同速率组别构建的本构模型参数
Tab.4 Construct constitutive model parameters for different rate groups
通过编写程序进行一维计算,而三维通过有限元商业软件建模进行计算,并设置与一维对应的材料属性,具体过程可见前期研究成果文献[15],计算结果如图6所示。
从图6可见,基于并行流变框架的非线性本构模型对于HTPB推进剂材料非线性力学性能有很好的预示能力,具体的误差将在2.3节中进行详细的讨论。但是由图6可见,非线性本构模型对低速率拉伸推进剂力学性能的预示能力并没有其对高速率拉伸推进剂力学性能预示能力强。分析原因可能是当推进剂材料在低速拉伸时,其内部的承载较均匀,应力波速率远大于拉伸速率,所以推进剂的曲线变化不大,基本呈现线性加载的现象(如图6中10 mm/min拉伸速率下,推进剂在10%的应变内力学性能基本为线性)。而非线性本构模型是基于宏观力学提出的本构模型,为了求得最优解,只能以大部分速率曲线的形态确定本构模型,故会出现低速率实验曲线出现偏差的情况。
图6不同速率组别本构模型计算结果
Fig.6Calculated results of different rates combinations of constitutive model
2.3 误差分析
取不同组速率模型一维理论计算数据与实验数据方差之和为评判模型精度的标准,即式(13)。表5为不同速率组别本构模型方差之和(误差)。
表5不同速率组别本构模型方差之和
Tab.5 Sum of variance of constitutive model for different rate groups
2.3.1 低速率组别对比影响分析
对比前四组实验结果,发现逐渐去掉低速率组别对实验误差的影响较大,结果如图7所示。原因是低速率组别反映了材料长期缓慢变化下的力学特性,大量的研究表明应变率变化对推进剂材料力学性能的影响较大[16]。低速率下的推进剂材料破坏主要以推进剂内部界面“脱湿”为主,“脱湿”现象导致推进剂颗粒从基体中析出,反映至宏观力学中即为材料的刚度下降,所以低速拉伸主要影响力学性能曲线的线性段部分,使得推进剂屈服强度偏小,初始模量偏小。而本构方程超弹性部分C10参数对于材料的初始模量较为敏感,所以低速率对于本构方程的影响不可忽略。
图7前四组速率组合本构模型误差结果
Fig.7Error results for the first four sets of rate combination constitutive model
2.3.2 高速率组别对比影响分析
对比第5组与第6组实验的误差值,3 000 mm/min的拉伸速率与4 000 mm/min的拉伸速率对模型的精度影响只有1%左右,故后续的实验可以将拉伸速率截断到3 000 mm/min,此时的实验数据可以用来建立精度较好的本构模型。由图2可见,3 000 mm/min与4 000 mm/min速率力学性能曲线已经较为接近,线性段(应变小于0.1)内两者基本重合。因为在高速拉伸的情况下,推进剂材料以基体的撕裂和颗粒的脆断为主要损伤形式,当速度达到一定时推进剂材料的失效形式已相对确定,其力学性能也会较为相似。
2.3.3 中间速率组别对比影响分析
后四组速率组合本构模型误差结果如图8所示,对比第5组、第6组、第7组与第8组实验的误差值,可以得出添加500 mm/min拉伸速率对模型的精度会有较大的影响。从图6(e)~(h)四图对比中可见,第7、8组大多数计算曲线比第5、6组计算曲线和实验数据曲线要低,意味着当添加500 mm/min拉伸速率时,所得本构模型的一维计算值与三维计算值会出现偏低的情况,即500 mm/min拉伸速率的添加使得本构模型出现偏低速的力学性能,导致所得本构模型的精度下降,也从侧面反映出在校核本构模型时速率组别的选择和组合方式十分重要。
图8后四组速率组合本构模型误差结果
Fig.8Error results for the last four sets of rate combination constitutive model
2.3.4 速率组别波动对于误差的影响
为了验证高低速率组合可以较好地构建模型,并且分析所选速率在最优速率周围波动时对模型精度的影响,在以上速率组别的基础上增加两组速率分别为:3 000 mm/min和100 mm/min、4 000 mm/min和100 mm/min。增设两组对照组意在与上文中最优组合3 000 mm/min和10 mm/min、4 000 mm/min和10 mm/min进行对比。将3 000 mm/min和100 mm/min组别设为第9组,同理将4 000 mm/min和100 mm/min设为第10组,图9为四组误差对比图。由图9可见,当把3 000 mm/min和10 mm/min组别(组别6)作为标准,其他三组误差百分比如表6所示。
从表6中可以看出当速率组别出现波动时,对于模型的精度有较大的影响。当速率组别在高速率周围波动时,模型的精度相较于速率组别在低速率周围波动要更高;当高低速率组别同时偏离最佳速率组别时,模型的精度会出现较大的下降;从第9组校核结果来看,用100 mm/min的拉伸速率数据代替10 mm/min拉伸速率数据可以提高模型的精度,从一种程度上说明实验速率的选择不是绝对的,但是高低速率组合的方式已能满足模型的校核需求。
图9速率组合波动对本构模型误差影响结果
Fig.9Results of the effect of rate combination fluctuations on the error of the constitutive model
表6不同速率组合相对最佳组合本构模型的计算误差
Tab.6 Computational errors of different rate combinations relative to the best combination of constitutive model
3 结论
本文通过对固体火箭推进剂不同应变率组别的对比分析,设置8组实验进行对照,并根据并行流变框架建立推进剂材料的本构模型,以理论模型与实验数据之间的方差之和为模型精度的评判标准,梳理了各种速率组合的精度情况,得出以下结论:
1)低速率组别如10 mm/min对于推进剂本构模型的影响不可忽略,建模时必须考虑推进剂在低速拉伸下出现的“脱湿”损伤;
2)高速率组别经过实验现象和建模误差对比分析,发现实验速率在3 000 mm/min时即可实现对本构模型较精确的校准,无须进行更高速率的实验;
3)高低速率组合可以实现对本构模型的精确校核,为推进剂材料实验给出了简化建议,提高了实验的效率,为快速预测材料的力学性能给出了一种方法。