粘弹性应变能释放率的虚拟网格计算方法①

崔辉如,姜强强,沈 栋,吕 轩,任孝宇

(1.陆军工程大学 国防工程学院,南京 210007;2.陆军工程大学 土木工程博士后科研流动站,南京 210007;3.中国人民解放军91112部队,宁波 315046;4.航天动力先进技术湖北省重点实验室,武汉 430040)

0 引言

当前,粘弹性材料已经被广泛地用于各种工程结构中。伴随着结构尺寸以及载荷环境的日益极端化,粘弹性结构的断裂问题显得尤为重要。以固体火箭发动机为例,其药柱属于典型的粘弹性材料。在日益频繁的战备值班以及实弹演练过程中,药柱结构内部很容易出现微裂纹,进而直接或间接导致药柱结构完整性的失效[1-2]。不同于弹性材料,粘弹性断裂问题尤其是粘弹性断裂参数的数值计算方法还有待深入发掘和研究。

断裂力学参数的计算是裂纹扩展判断的重要依据。应力强度因子以及应变能释放率等断裂力学参数的计算大部分依赖于原始的网格或者构建虚拟的区域进行求解。在原始网格的基础上,最早最直接的方法是基于裂纹尖端前沿的单元应力或裂纹尖端后的节点位移的外推法[3]。王珊[4]利用有限元节点位移确定含裂纹Kirchhoff板辛本征解中的待定系数,计算出了裂纹尖端附近应力场的显式表达式。段静波等[5]将加料有限元法扩展应用于粘弹性介质中裂纹应变能释放率的计算。祁勇峰等[6]在裂纹尖端引入Williams解析解级数,应用解析覆盖与周边数值覆盖联合求解三维穿透直裂纹的应力强度因子。RAJARAM等[7]采用积分点处的摄动梯度与直接解析微分之间的联系,DAIMON等[8]引入一个使得辅助解满足积分平衡关系式的修正项来进行强度因子的求解。JIN等[9]利用弹性-粘弹性对应原理研究了粘弹性梯度材料Ⅰ型及Ⅱ型裂纹的应力强度因子解析解。另外,还有研究采用虚拟网格来进行断裂参数的求解。例如,采用沿裂纹前缘的虚拟圆柱单元,利用极坐标下的高斯点在圆域中进行J积分的计算[10-11]。为避免网格中应力强度因子沿裂纹前沿的振荡,可以利用被积函数的移动最小二乘近似[12]。类似地,还有研究采用盘形域积分法逐点计算J积分和强度因子[17]。GROSSMAN-PONEMON等[14]在虚拟网格中展示了一种基于独立网格的多网格交互作用积分方法。CHOI等[15]利用虚拟网格插值计算虚拟网格裂纹尖端位移,通过计算虚拟网格J积分,进而获取了低质量网格下二维以及三维弹性材料应力强度因子。

裂纹尖端断裂力学参数计算已经有不少的实现办法和实例支撑,但是目前还有以下不足:一是研究对象大多针对线弹性材料,对粘弹性介质断裂问题缺乏深入研究;二是伴随着裂纹扩展过程中非结构且低质量网格的使用,上述方法实现流程相对复杂,并且难以确保仿真精度。本文在前期研究的基础上[15],在裂纹尖端构建了一种新的虚拟网格形式,结合虚拟裂纹闭合方法,发展了一种针对粘弹性材料裂纹尖端应变能释放率的高精度数值计算方法。利用两种典型裂纹案例与文献验证了本文方法的有效性。

1 线粘弹性本构模型

在泊松比为常数的条件下,推进剂三维线粘弹性本构关系表达式为[16]

(1)

其中,

(2)

(3)

式中ν为泊松比;σij、Sij及σkk分别为应力张量、偏应力张量及球应力张量;eij及εkk分别为偏应变张量以及球应变张量;t为加载时间;E(t)为松弛模量。

E(t)Prony级数表达式为

(4)

2 粘弹性裂纹尖端虚拟网格

如图1,以原始网格裂纹尖端P点为中心构建虚拟网格。虚拟网格由2×2个正方形单元构成(D2PKJ、PHNK、EFPD1以及FGHP),虚拟单元的边长为原始网格裂纹尖端平均单元边长的2倍[15]。与此同时,虚拟网格需要在原始裂纹路径上设置同样的裂缝,因此虚拟网格一共有10个虚拟节点。

图1 原始网格(黑色三角形网格)与虚拟网格(蓝色四边形网格)示意图

虚拟节点的位移计算遵循以下原则。假设,原始网格某单元的组成节点为A、B和C,那么虚拟节点N刚好位于该单元区域的判断依据是

SABN+SBCN+SCAN=SABC

(5)

式中SABN、SBCN、,SCAN、SABC分别为三角形ABN、BCN、CAN及ABC的面积。

假设原始网格上单元A、B和C三个节点位移分别是(uA,vA),(uB,vB)以及(uC,vC)。N点、A点、B点以及C点的坐标为分别为(xN,yN)、(xA,yA)、(xB,yB)以及(xC,yC)。

利用面积法进行虚拟节点N位移(uN,vN)的插值计算

(6)

其中,

(7)

以此类推,不难计算出任意时间步上虚拟网格的位移场。

考虑到粘弹性材料具有典型的时间相关特性,为此在进行虚拟网格应力场计算的时候,不仅需要代入上一步计算中的应力、应变增量等数据还要将当前步产生的应力、应变增量等数据传递到下一计算步中。有关粘弹性本构模型的数值化计算方法可以参考相关文献,这里不再详细展开[17-18]。

虚拟网格应力计算以及原始网格的应力计算都需要进行相关状态变量的数据传递,上述操作会给主程序的设计和编程带来很大的不便。考虑虚拟网格中节点和单元总数较少,可以在每一步计算中利用相同的位移插值计算方案,计算出虚拟网格在当前步以及前序时间步中的位移场,进而直接计算出当前载荷步上虚拟网格的应力以及应变场。上述操作可以有效降低时间步之间状态变量传递的工作量,简化编程任务。为了方便表示,将该虚拟网格应力重构的方法记作(Stress Recovery Technique,SRT)。图2给出了本文中粘弹性应变能释放率数值计算的流程示意。

图2 应变能释放率数值计算流程示意图

本文中的虚单元均采用常规一阶四节点四边形单元。虚单元的单元积分方案有两种,分别是2×2高斯积分点的完全积分以及1×1高斯积分点的缩减积分。如无特殊交代,虚单元采用完全积分方案。此外,本文中的原始网格均采用三节点三角形单元,积分点数为1。

3 虚拟裂纹闭合方法

Irwin在假设“势能改变与裂纹闭合一个扩展增量所需的功等效”基础上,利用裂纹闭合积分计算裂纹尖端的能量释放率[3]:

(8)

式中(见图3)a和Δa分别为裂纹长度和裂纹扩展增量;σyy和τxy分别为沿着闭合裂纹面上的法向和切向应力;T为裂纹体厚度;Δv和Δu分别为当闭合裂纹张开时闭合裂纹面上的法向和切向位移。

图3 闭合的虚拟裂纹示意图

上式在实际的计算过程中需要开展两次分析计算,这给裂纹扩展问题的研究带来很大的不便。为此,RYBICKI等[19]首次提出了应变能释放率计算的虚拟裂纹闭合方法。该方法假设虚拟裂纹尖端后面的张开位移和实际裂纹尖端后面的张开位移近似相等。此外,为了避免积分号中出现不精确的应力值,采用节点力代替应力的积分,即

(9)

式中(见图4)Fx1和Fy1分别为节点1处内力沿着水平和垂直方向的分力;Δv2,3和Δu2,3分别为节点2和3之间垂直和水平方向的位移差。

图4 虚拟裂纹闭合方法计算应变能释放率示意图

在上节虚拟网格的基础上,利用编程很容易实现节点1处支反力以及节点2和3处位移的计算,进而计算出任意时刻裂纹尖端的应变能释放率。

4 验证算例

图5 三参数线粘弹性模型

对于上述粘弹性结构受阶跃载荷时,张淳源[21]给出了结构裂纹尖端应变能释放率的解析表达式:

(10)

式中C(t)为材料的蠕变柔量;平面应变时,χ=1.0-v2,平面应力时,χ=1.0;KⅠ和KⅡ分别为结构所受阶跃载荷变成定常载荷时裂纹尖端的Ⅰ型和Ⅱ型裂纹强度因子。

对于经典的三参数线粘弹性模型,其蠕变柔量的形式为

(11)

4.1 单边裂纹受拉测试

图6(a)为一单边含裂纹缺陷的粘弹性平板,板的长度和宽度分别2h=20和b=10。平板的厚度为1,板中间位置设置长度为a=1的水平裂纹。平板两端受均布拉应力σ(t)=σ0H(t)作用,H(t)代表单位阶跃函数,σ0=1,作用时间为200 s。对于定常均布拉应力σ0作用,裂纹尖端的Ⅰ型强度因子满足[22]:

(a) Geometric model of tensile specimen with single edge crack (b) Semi-structured low-quality grid model

(12)

其中,

(13)

为方便,定义相对误差用以表示计算结果的准确性。各个时刻相对误差表达式如下

(14)

图6(b)给出了平板的非结构网格划分情况,网格划分均采用一阶三角形单元。在远离裂纹尖端的平板边缘,单元尺寸为2.0,在裂纹尖端附近对网格进行了加密,单元平均尺寸为0.05。

结合裂纹尖端的虚拟网格技术(虚拟单元尺寸为0.1)与虚拟裂纹闭合方法,分别开展泊松比为0.3时,平面应力以及平面应变条件下的Ⅰ型应变能释放率的计算。从图7可以明显地看到,本文提出的虚拟网格技术比文献[20]中的加料单元方法相比精度大幅提高,相对误差均在5%以内。加料单元方法的本质是将断裂参数作为单元的附加自由度,即所谓的“加料”。该方法需要对裂纹尖端相邻单元的位移模式进行精细化设计。含裂纹结构有限元分析结束后,主程序将会同时获得断裂参数和位移场的结果。但是,加料后的单元并没有改善裂纹尖端应力奇异的现状,断裂参数的计算还是在裂纹尖端奇异应力的基础上获得的。这是本文提出的虚拟网格技术比加料单元方法精度更高的主要原因。以上算例说明,本文提出的方法在应对低泊松比条件下的Ⅰ型应变能释放率计算是准确可行的。

(a)State of plane stress (b)State of plane strain

对于推进剂这一类典型的粘弹性材料,其泊松比一般大于0.495,甚至接近0.5。在利用虚拟网格技术时,分别采用两种数值积分手段,进行积分点上应力应变的数值计算。一种是采用全积分形式,即采用4个积分点进行计算。考虑到虚拟单元方法的本质是对裂纹尖端应力场的重构,为此将全积分形式的虚拟网格方法记作(Stress Recovery Technique-Fully Integration,SRT-F)。另一种方案采用缩减积分的形式,即采用一个高斯积分点,将该方法记作(Stress Recovery Technique-Reduced Integration, SRT-R)。

图8给出了平面应变状态时两种积分方案下,应变能释放率的计算结果与文献解的对比情况。由图8可以看到,虚拟网格应力的全积分计算方案效果较差,相对误差接近100%。与同样采用缩减积分方案的文献解相比,本文的虚拟网格技术获得的应变能释放率精度更高,在加载后期,相对误差低于1%。由此可见,虚拟网格缩减积分技术可以应对高泊松比条件下的Ⅰ型应变能释放率的数值计算。

图8 平面应变状态阶跃加载下应变能释放率相对误差随时间变化(v=0.495)

4.2 单边裂纹受剪测试

图9(a)为一单边含裂纹缺陷的粘弹性平板,板的长度和宽度分别2h=20和b=10。平板的厚度为1,板中间位置设置长度为a=3的水平裂纹。平板左侧边缘中间位置分别受一对相互平行集中力Q(t)=Q0H(t)作用,Q0=1,作用时间为200 s。对于定常集中力Q0作用,裂纹尖端的Ⅱ型强度因子满足[23]:

(a) Geometric model of shear specimen with single edge crack (b) Semi-structured low-quality grid model

(15)

其中,fⅡ=2.5935。

图9(b)给出了平板的非结构网格划分情况,网格划分同样采用一阶三角形单元。在远离裂纹尖端的平板边缘,单元尺寸为2.0,在裂纹尖端附近对网格进行了加密,单元平均尺寸同样为0.05。

图10给出了平面应变状态时两种积分方案下,应变能释放率的计算结果与文献解的对比情况。虚拟网格应力的全积分计算方案效果较差,相对误差在90%左右。本文的虚拟网格技术获得的应变能释放率相对误差低于7%。另外一方面,文献解的精度较差,加载后期一度超过40%。可以认为,本文提出的虚拟网格技术能够应对高泊松比条件下的Ⅱ型应变能释放率的高精度数值计算。

图10 平面应变状态阶跃加载下应变能释放率

进一步地,为了检验复杂载荷历程下,虚拟网格技术的计算能力,按照文献[20],设计了如图11两种载荷历程形式,其中相互平行集中力Q(t)=Q0f(t)。

(a)Load history I (b)Load history II

(a)Load history I (b)Load history II

因此,文献解与本文提出的虚拟网格全积分求解方案相对误差较大。第50 s时刻,本文提出的虚拟网格缩减积分求解方案计算出的应变能释放率为0.499 2,可以认为这是三种方案中,精度最高的求解手段。类似地,在载荷历程二中,应变能释放率的最大值不会超过0.538 8,文献解与虚拟网格全积分求解方案存在较大偏差。

5 结论

本文提出了一种粘弹性应变能释放率的虚拟网格求解方案,通过两种典型含裂纹缺陷,结合解析表达式,验证了求解方案的准确性和有效性,主要结论如下:

(1)虚拟网格技术结合虚拟裂纹闭合方法可以有效地利用低质量网格实现粘弹性裂纹尖端应变能释放率的高精度数值计算。与文献提出的加料单元方法相比,实现思路简洁、计算精度更高并且容易编程实现。此外,由于不需要单独设计网格,本文提出的方案更加适用于粘弹性裂纹扩展过程中的应变能释放率计算。

(2)对于泊松比接近于0.5的粘弹性材料,仿真结果表明,本文提出的虚拟网格缩减积分方案同样可以高效地应对。

(3)本文提出的二维粘弹性应变能释放率求解方案可以移植和扩展到任意材料二维裂纹前缘的应变能释放率计算问题中。对于三维应变能释放率的求解,可以参考本文方法,在裂纹面前缘构建三维虚拟网格并结合虚拟裂纹闭合方法,实现三维裂纹面上断裂参数的数值计算。此外,本文方法将应用到粘弹性材料的裂纹扩展仿真中,为任意低质量网格条件下的裂纹扩展提供更加精确和科学的判断。

(4)与文献[5]求解方案相比,本文提出的粘弹性应变能释放率虚拟网格计算方法实现方便、精度高并且易于在后期的裂纹扩展仿真研究中实现。