HTPB黏弹性微裂纹偏折扩展损伤本构模型

2018-09-29顾志旭郑坚彭威支建庄

航空学报 2018年9期

顾志旭，郑坚，彭威，支建庄

陆军工程大学火炮工程系，石家庄 050003

端羟基聚丁二烯(HTPB)推进剂是一种广泛应用于航空航天和军事武器推进系统的颗粒增强型复合材料，其力学行为直接影响到弹箭武器的使用可靠性[1]。建立准确高效的本构模型，是分析其力学行为的重要前提。

颗粒增强型复合材料的力学行为，不仅取决于组分的性质，更依赖于组分间的细观结构。当加载应力或应变较小时，材料通常表现出基体的黏弹性特征，而当应力或应变较高时，颗粒/基体界面容易出现“脱湿”损伤而形成微空洞，微空洞生长融合形成微裂纹，微裂纹扩展汇合形成宏观裂纹，材料又表现出非线性特征[2]。为了描述上述黏弹性损伤行为，许多学者建立了不同的本构模型。一方面，Xu等[3]基于热力学理论和弹性-黏弹性对应原理建立了热黏弹性损伤本构模型，该模型能够预测宽泛温度范围下推进剂的本构行为。Yun等[4]引入软化和循环加载效应函数，建立了可以反映材料Mullins效应的损伤本构模型。此外，还有Park[5]、Ozupek[6]和Jung[7-8]等的一些工作。上述模型主要是从宏观角度出发，忽略了材料损伤的具体细节，建立的各向同性损伤模型，通常具有较高的计算效率，但模型中往往含有若干物理意义并不十分明确的待定参数，如Xu模型[3]和Park模型[5]中的C(S)线。另一方面，Chen[9-10]，Tohgo[11]，Xu[12]以及Hur[13]等认为颗粒/基体界面的脱湿是材料损伤的主要机制，利用空洞代替脱湿的颗粒，基于细观均匀化理论，建立了不同的细观模型。然而在低应力三轴度时，颗粒/基体界面的脱湿并不完全，颗粒在界面完好的方向仍具有承载能力，上述模型难以反映材料此时表现出的宏观各向异性特征。从另一个角度，彭威等[14]将界面的脱湿视为微裂纹的扩展，采用微裂纹扩展区的概念，建立了一个能够表征各向异性损伤的本构模型，但因需要确定多个损伤内变量，并不方便于应用。同时，由于采用了应力强度因子作为微裂纹扩展的判据，该模型不能反映材料的蠕变损伤特征。

由于HTPB材料的高填充比特性(高达85%)，脱湿形成的微空洞容易融合而形成微裂纹[15]。因此，本文在一个较大的观察尺度上，将材料的主要损伤机制视为微裂纹的扩展和汇合，而将颗粒脱湿、微空洞变形等细观损伤视为裂尖失效区的演化变形。进一步，在Abdel-Tawab宏观本构方程的基础上，考虑微裂纹偏折扩展的情形，建立了一个能够反映材料损伤的应变率和温度依赖性，以及各向异性损伤特征的黏弹性本构模型。

1 宏观黏弹性损伤本构方程

Abdel-Tawab和Weltsman[16]基于不可逆热力学和叠加原理，假定损伤对材料柔量的平衡和瞬态部分具有相同的影响，推导出了如式(1)的本构方程。

(1)

式(1)表明，可以通过有效应力来考虑损伤的影响，但应用的关键在于确定损伤映射张量的具体形式。下文基于微裂纹均匀化理论，推导损伤映射张量的一般的形式。

2 损伤映射张量和损伤内变量

定义损伤内变量为

(2)

式中：ai为微裂纹的半径；ni为裂纹面的单位法向矢量；η为体密度；g(ai)和r(ni)分别为尺寸和取向分布函数；N为微裂纹群数(定义取向和尺寸均相同的微裂纹为一个群)；⊗表示并矢。

根据稀疏估计方法，当忽略摩擦效应时，弹性材料中微裂纹引入的柔度增量可以近似为[17]

ΔSijkl≈

(3)

式中：E0和υ0分别为无损材料的弹性模量和泊松比；ωik为损伤内变量ω的分量；δjl为Kronecker符号。

(4)

式中：S0为无损弹性材料的柔度张量；ΔS(ω)为式(3)定义的柔度增量的张量形式。

根据式(1)中损伤映射张量的定义，有

(5)

将式(3)和式(4)代入式(5)，整理得损伤映射张量为

(6)

记Stest为实验测定的柔度张量，Scurrent为微裂纹扩展至当前状态时的损伤柔度张量，则根据损伤映射张量的定义，有

Stest=S0:Pinitial

(7)

Scurrent=S0:Pcurrent

(8)

式中：Pinitial和Pcurrent分别为初始微裂纹和当前微裂纹状态对应的损伤映射张量，由式(6)求得。显然，S0和Scurrent是未知的。在实际应用时，需要找到由Stest向Scurrent映射的损伤映射张量。根据式(7)和式(8)，不难得到

(9)

3 黏弹性裂纹扩展演化

3.1 能量扩展准则

(10)

根据式(10)，裂纹尖端的张开位移场可以表示为

(11)

式中：ρ为极径；υ为泊松比。

同时，应力场为

(12)

裂纹扩展长度为Δc时所需能量为

(13)

式中：σγy和uγ分别由式(12)和式(11)给出；s为裂纹扩展长度。

考虑Ⅰ型开裂，令φ=0，将式(11)和式(12)代入式(13)，有

(14)

式中：c为裂纹扩展Δc前的尺寸。

裂纹扩展的能量释放率为

(15)

将式(15)推广，考虑复合型开裂时，有

(16)

由式(10)和式(16)可见，在蠕变加载下，裂尖的能量释放率为加载时间的递增函数。若在加载初期稳定的裂纹(G

3.2 裂纹偏折

假定初始缺陷半径为a的钱币形裂纹，则裂尖前缘任一点Q处的应力强度因子为[19]

(17)

将式(17)代入式(16)，得Q点处的能量释放率为

GQ(φ)=A1+A2sin2φ+A3cos2φ+A4sin(2φ)

(18)

式中：

假定在裂尖法向平面内，裂纹沿周向应力最大的方向进行扩展，则裂纹前缘不同点处的偏折角θk满足:

(19)

式(19)表明，不同点处的偏折角通常并不相同。同样，A和A′点的偏折角为最大，而B和B′点的偏折角为最小0(发生自相似扩展)。故偏折扩展后的币形裂纹如图2所示，裂纹前缘为一椭圆周线，长轴为CC′，短轴为DD′。

3.3 等效裂纹

对于图2所示的三维偏折裂纹，分析其对材料柔量的影响时，存在着较大的困难。简化起见，采用Francois和Dascalu[20]的方法，将其等效为长轴为CC′，短轴为DD′的椭圆形裂纹，并且有

(20)

(21)

式中：an和an+1分别为扩展前后的尺寸，如图3所示；l为偏折扩展的长度，由3.4节扩展速率方程的确定；θkmax为币形裂纹在A点的偏折角。

为了沿用损伤内变量式(2)的定义，根据Sevostianov和Kachanov[21]的研究，用币形裂纹来近似上述椭圆形裂纹，两者的几何关系为

(22)

式中：L为等效币形裂纹的半径；k=f/e,e和f分别为椭圆的长半轴和短半轴；Ξ(k)为第一类完全椭圆积分。

3.4 扩展速率方程

Schapery[22]基于图4所示的裂纹模型，得到的裂纹扩展速率方程为

(23)

式中：Π为裂纹扩展所需的断裂能；σm为失效应力；KI为Ⅰ型应力强度因子；C1和n为平面蠕变柔量Cv(t)≈C1tn时对应的系数。图4所示中，P为真实裂尖，P′为表观裂尖，α为裂尖失效区长度，d为裂纹张开位移。

假定上一t1时刻，真实裂尖到达H点，表观裂尖在I点，如图5所示。当前时刻t，裂纹向前扩展α长度，真实裂尖和表观裂尖分别到达J和H点。

在Δt=t-t1时间段内，H点的张开位移由0增加到dm[24]，

(24)

假定KI与t呈幂函数关系，即KI～tμ，并注意到Cv(t)≈C1tn，则有

(25)

式中：r=2μ；κn=Γ(r)Γ(1+n)/Γ(1+n+r)。

裂纹扩展长度为α时所需断裂能为[24]

(26)

将式(26)代入式(23)，可得

(27)

(28)

式(28)即为裂纹扩展的速率方程。

4 模型验证

利用ABAQUS二次开发接口UMAT，将上述模型嵌入其中，进行分析计算。本构方程式(1)的离散方法参见文献[25]，而速率演化方程式(28)采用Euler方法求解，以节省计算量。损伤内变量的计算更新流程如图6所示。

假定材料初始状态为各向同性，初始微裂纹空间均匀分布，尺寸为其统计平均值a0。模型参数a0、σm和N对本构曲线的影响，同微裂纹自相似扩展的情形[26]，此处不再赘述。

T/℃Gc/(mJ·mm-2)σm/MPaε=5 min-1ε=0.5 min-1ε=0.05 min-1ε=0.005 min-1-400.013.652.321.150.88 200.0010.780.510.420.28 600.000 10.450.370.300.20

表2 临界能量释放率Gc的拟合值Table 2 Fitted values of critical energy release rate Gc

不同温度和应变率下，临界能量释放率Gc的拟合值见表2。可见在同一温度下，Gc受应变率的影响有限，最大变幅为17%。故在表1中将其近似取为常数，并不影响预测结果的准确性，如图7所示。然而，其受温度的影响较大，随着温度的升高而减小。

为便于说明，定义损伤应力-应变曲线偏离无损曲线的起始点为临界损伤点(图7中虚线与实线的交叉点)，该点对应的应力为临界损伤应力，记为σd。不同温度和应变率下，σd的近似估计值见表3。结合表1可见，失效应力σm取值约为临界损伤应力σd的1～1.3倍，并且随着应变率的降低和温度的升高而减小。这是由于，一方面本文采用了稀疏估计的方法，微裂纹周围的平均应力就为远场加载应力。另一方面，受裂尖应力放大效应和脱湿卸载效应的共同影响，裂尖失效区内的平均应力大于微裂纹周围的平均应力。

表3 临界损伤应力σd的近似值Table 3 Approximate values of critical damage stress σd

图8给出了20 ℃，5/min应变率加载下，部分微裂纹的扩展情况。由图8可见，在单轴x2向加载下，初始θ=90°的微裂纹不发生扩展，而初始θ=7.95°的微裂纹扩展尺寸最大，介于两者之间的微裂纹，扩展尺寸随着初始θ角的增大而减小，等效微裂纹向垂直于加载轴的方向偏转。

(29)

(30)

最后，检验3.4节中的基本假设KI～tμ。图10给出了初始θ=7.95°,58.59°的2个微裂纹，分别对应图1中A和B点的KI-t曲线。其中离散点为UMAT计算的KI值，实线为KI～tμ关系的拟合曲线。由图10可见，对于不同取向的微裂纹，假设KI～tμ近似成立，但μ的取值并不相同。为便于应用，进一步假定不同微裂纹的μ值相同。当μ取0.5时，式(28)与文献[24]中的式(76)相同，此时的KI～t0.5拟合曲线如图9中虚线所示。可见，在一定的误差范围内，可以用KI～t0.5来近似KI～tμ，而不影响模型的预测精度，见图7。