基于相场法的氢脆断裂模拟

2022-07-21陈鹏程马玉娥彭帆周玲珑

西北工业大学学报 2022年3期

陈鹏程, 马玉娥, 彭帆, 周玲珑

(西北工业大学航空学院, 陕西西安 710072)

腐蚀化学环境会诱发金属结构产生氢脆(HE,hydrogen embrittlement),在与力学环境共同作用下,会更容易造成结构的失效和破坏。相比于单纯的断裂破坏,氢脆断裂在断裂前的征兆不易观察,临界破坏载荷更小,裂纹扩展速度更快,因此带来的破坏性更强[1]。氢脆断裂造成的工程失效案例屡见不鲜,涉及航空、航天、储能等重大战略性领域。因此研究氢脆断裂问题有重要的实际应用价值。

对于单纯的断裂问题,工程上一般采用有限元仿真来实现对断裂过程的模拟和预测,但是传统的有限元方法(FEM)在处理断裂问题时需在裂纹扩展过程中不断追踪裂纹界面,反复对裂尖进行网格重新细分才能精确计算裂尖奇异场[2],因而在处理断裂问题时存在计算复杂、精度过低等缺点。

变分断裂相场法则可以很好地弥补FEM这些缺点:断裂相场法通过序参量的自动演化来确定裂纹面,因此可以避免追踪裂纹面的繁琐过程。此外,相场法还可以在不引入额外的失效准则条件下确定裂纹启裂和扩展的方向[3-4]。

与普通断裂不同,氢脆断裂伴随复杂的氢扩散过程,需考虑氢浓度场演化对材料断裂韧性的影响。值得注意的是,相场断裂模型的实质就是相场和位移场相互耦合的偏微分方程组,不同场间的影响可以通过相应的控制方程体现,因此利用相场法来处理多场耦合下的断裂问题十分有效。

目前有关相场氢脆断裂的主要研究有:利用基本相场断裂模型,Martínez等[5]提出了相场氢脆断裂的框架;Wu等[6]基于统一相场理论[7]将相场正则化内聚裂缝模型应用于氢脆断裂的模拟;而Mandal等[8]则是对比研究了不同相场模型下的氢脆断裂计算结果,深入研究了相关参数对于氢脆断裂计算结果的影响。

在基本的相场断裂模型[5,9-11]中,受拉受压均会造成能量的释放,这在计算中会造成一定误差。因此其应用局限于单纯拉伸载荷作用下的断裂模式。在模拟更一般的断裂时,需要分别考虑拉压应力分量的影响。本文则基于Miehe等[9]提出的应变能拉压谱分解方法,改进了相场氢脆断裂模型,并通过数值计算验证了改进模型的可靠性。

1 基本模型

1.1 相场断裂基本模型及控制方程

相场法通过引入一个辅助标量场——相场φ来表征裂纹拓扑[9]。在一维情况下,φ的表达式为

γ(φ)=12l0φ·φ+1l0φ2]

(2)

由此可得正则化表面能

(3)

式中，Gc(θ)为临界能量释放率。引入应力退化函数g(φ)=(1-φ)2+k,其中k是为了避免计算奇异而引入的附加参数,在具体计算中尽量取足够小,本文中均取k=1×10-6。则原应变能ψe(ε)在考虑相场影响下有

(4)

总势能泛函则为

Ψ(φ,u)=Ψs+Ψb

(5)

在表面力f和体积力b作用以及边界约束下,可以推导准静态下相场断裂位移场和相场控制方程

(6)

式中：σ是应力张量;n是表面法向量。

1.2 氢离子浓度扩散方程

氢脆断裂发生时,氢离子会向材料内部应力集中部位扩散聚集,大大降低构件的断裂韧性。根据相关的研究[5],表面氢离子浓度θ对材料断裂韧性的影响可用(7)式表征

(7)

Gc(0)为氢离子浓度为0时材料的临界能量释放率。χ是损伤系数,代表氢离子对断裂韧性影响。式中表面氢离子浓度θ可由体积氢离子浓度C计算

(8)

(9)

式中：J是氢离子浓度通量；q是表面氢通量。J为关于应力梯度和氢离子浓度的函数,可以采用如下公式定义

(10)

图2 浓度扩散示意图

2 改进后相场氢脆断裂模型

2.1 应变能拉压分解

在上述基本相场氢脆断裂模型中,相场引起的应力退化会引起应力张量所有的分量发生改变,包括拉伸分量和压缩分量,这与实际的断裂机制有出入。为了将此模型推广,本文引入Miehe等[9]提出的应变能谱分解方法,重新考虑相场对应变能的折减方式。首先对应变张量进行谱分解

(11)

式中，〈·〉为麦考利符号。相应地,对弹性应变能进行拉压分解,这样相场折减只对正应变能有效

(12)

(12)式对应变张量求偏导,得到谱分解下的应力张量

(13)

λ和μ是拉梅常数,I是二阶单位张量。则弹性张量为

(14)

2.2 控制方程的离散化

(15)

N是形函数矩阵,B是其对应的偏导矩阵。

2.3 求解策略

利用Newton-Raphson算法结合分步式算法,对(15)式中的非线性方程组求解

(16)

式中

(17)

时间导数项对应的刚度矩阵M的表达式为

(18)

采用分步算法,可以完成相场和其余两场的解耦。同时,本文为了提高精确性,在计算过程中考虑浓度场和位移场的耦合项刚度矩阵,其推导过程如下

(19)

由刚度矩阵的表达式可以看出,在计算过程中需要确定积分点的静水压力值以及其梯度值。本文采用同样形式的形函数Ni来计算静水压力

(20)

式中,σH,i是节点i处的静水压力,其表达式为

(21)

d是单元节点个数,则有

(22)

利用Voigt标记法将四阶弹性张量Ds写成矩阵形式,则平面应变问题下的耦合项刚度矩阵公式为

(23)

式中，L的表达式为

(24)

为了验证改进方法的有效性,基于Matlab程序开发平台,编写了多场耦合断裂的有限元求解程序,并设计算例进行验证。

3 数值算例

为了验证改进模型的可行性以及对Ⅰ型氢脆断裂的模拟情况，选取含单条裂纹的正方形板作为研究对象,分别施加拉伸位移载荷和剪切位移载荷,如图3～4所示。试样的边长为1 mm,预制裂纹长度均为0.5 mm,位于试样的正中间,底部为固定位移边界。

图3 单裂缝正方形板图4 单裂缝正方形板受拉伸位移载荷受剪切位移载荷

3.1 准静态Ⅰ型裂纹问题

在单纯的拉伸位移载荷下,整个区域内设置固定氢离子浓度条件C,分别取4种不同的浓度边界条件C=0,0.1×10-6,0.5×10-6,1×10-6来研究氢离子浓度对断裂过程的影响。图5对比了不同浓度下的位移-载荷曲线,可知断裂过程对氢离子浓度十分敏感:随着浓度的增加,临界破坏载荷不断减小。以浓度为0.5×10-6时为例,如图6所示,可以看出在断裂过程中,氢离子不断向应力集中处扩散;在相同位移加载下,浓度越大则相场演化越快,如图7和图8所示。同时为了验证本文模型的准确性,选取了浓度为0.5×10-6时的计算结果与Martínez等[5]的研究结果进行对比,如图9a)所示。结果表明,改进模型所计算临界破坏载荷与原模型相差仅为3.954%,因此,改进后模型能满足计算精度;同时可以看出,改进后模型能捕捉脆弹性材料的瞬断过程。

图5 不同氢浓度下的位移-载荷曲线

此外,为了研究位移加载步对计算结果的影响,在C=0.5×10-6的浓度边界条件下,分别选取了Δu为8×10-5,4×10-5,2×10-5,1×10-5,5×10-6mm 5种不同位移加载步值,计算了相应加载步下的破坏过程曲线,如图9b)所示。在本文中采用分步式算法求氢脆相场断裂模型时,位移场和浓度场完全耦合,它们和相场的求解在前后相邻的加载步间进行。因此,大的加载步将影响求解的精度。如图9b)所示,当加载步逐渐减小,计算结果趋近于稳定,但是计算时间成本不断增加。因此在选取加载步时需要权衡计算精度和计算效率,而本算例所选取的加载步,其对应的计算结果如图9b)红色粗实线所示,综合考虑其计算精度和计算效率,是个较理想的选择。

图6 浓度为0.5×10-6时浓度扩散示意图(Ⅰ型)

图7 浓度为0.5×10-6时相场演化示意图(Ⅰ型)

图8 浓度为0时相场演化示意图(Ⅰ型)

图9 位移-载荷曲线

3.2 准静态Ⅱ型裂纹问题

引入拉压分解后,可以更精确地模拟Ⅱ型断裂过程。本例中,在剪切位移载荷下,分别取2种不同的浓度条件C=0,0.5×10-6来考虑氢离子浓度对断裂过程的影响。以0.5×10-6为例,浓度扩散和断裂(相场)演化的过程分别如图10～11所示。由图可见,引入拉压分解后的模型可以有效地模拟Ⅱ型氢脆断裂过程,相较于原模型,采取应变能拉压分解后的模型不会出现裂纹的伪分叉[12]。同时,图12给出了浓度为0×10-6同位移载荷下的相场演化模拟图。结合图13不同浓度下断裂过程的位移-载荷曲线,可见在Ⅱ型断裂模式下,整个断裂过程对氢离子浓度同样很敏感。

本节同样进行了位移加载增量步对Ⅱ型相场氢脆断裂计算结果的影响分析。如图14所示,加载步愈小,最终计算结果愈精确。可以看出,在本算例中,综合计算精度和计算效率考虑后所选取的加载步Δu=1×10-5mm,较为合理。