基于相场法的氢脆断裂模拟
陈鹏程, 马玉娥, 彭帆, 周玲珑     
西北工业大学 航空学院, 陕西 西安 710072
摘要: 通过引入应变能拉压分解, 对相场氢脆断裂模型进行了改进, 并给出了改进模型的数值求解公式, 推导了浓度场和位移场耦合项刚度矩阵。编写Matlab相场氢脆断裂计算程序, 分别模拟了Ⅰ型和Ⅱ型2种不同模式的氢脆断裂过程。数值计算表明: 在断裂过程中氢离子会向裂尖处扩散集中, 且氢浓度越高, 临界破坏载荷越小; 对比已有模型计算结果, 改进模型能精确计算Ⅰ型断裂问题的临界破坏载荷, 准确捕捉脆弹性材料的瞬断过程, 且能更有效地模拟Ⅱ型的氢脆断裂过程。
关键词: 相场法    多场耦合    氢脆    拉压分解    脆弹性断裂    

腐蚀化学环境会诱发金属结构产生氢脆(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]。在一维情况下, ϕ的表达式为

(1)

式中,a表示尖锐裂纹发生位置。当ϕ=0时, 表示构件在该处完整无裂纹; 当ϕ=1时, 则表示完全开裂。l0是表征裂纹弥散程度的参数。由图 1可见, l0越大, 相场弥散范围越宽。

图 1 标量场ϕ的分布以及l0对标量场分布的影响

由此推广到二维和三维, 引入裂纹表面密度函数

(2)

由此可得正则化表面能

(3)

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

(4)

总势能泛函则为

(5)

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

(6)

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

1.2 氢离子浓度扩散方程

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

(7)

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

(8)

式中:Δgb0为吉布斯自由能;R是气体常数;T是温度。浓度扩散过程如图 2所示, 根据守恒定律, 扩散控制方程可以表示为

图 2 浓度扩散示意图
(9)

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

(10)

式中:D是浓度扩散系数;VH是固溶体中氢的偏摩尔体积;σH是静水压力。

2 改进后相场氢脆断裂模型 2.1 应变能拉压分解

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

(11)

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

(12)

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

(13)

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

(14)

式中, H(x)为阶跃函数: 当x < 0时H(x)=0, 当x≥0时, H(x)=1; 为球形四阶单位张量; 投影算子。由(14)式可知, 弹性张量会随着加载步的增加而不断更新。

2.2 控制方程的离散化

考虑氢离子浓度对断裂韧性的影响(见(7)式),引入2.1节中拉压分解下的应变能、应力以及弹性张量的新形式。将1.1节和1.2节中的控制方程弱化并且离散后得到位移场残差riu, 相场残差riϕ以及浓度场残差riC分别如下

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

(23) 式中κ为材料的体积模量。(24)式中Dijd表示节点d处的弹性矩阵第ij列值。Ni, xd表示节点d处的Nix的偏导值; 而耦合项

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

3 数值算例

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

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

在数值计算过程中, 参数取值如下: ①材料弹性模量取E=210 GPa, 泊松比ν=0.3。②相场参数: 弥散参数l0=0.05 mm, k=1×10-6。临界能量释放率取Gc(0)=2.7×10-3 kN/mm。③损伤系数χ=0.89, 扩散系数D=0.012 7 mm2/s, VH=2 000 mm3/mol。④数值模拟的室温取T=300 K, 位移增量步取Δu=1×10-5mm。

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 不同氢浓度下的位移-载荷曲线
图 6 浓度为0.5×10-6时浓度扩散示意图(Ⅰ型)
图 7 浓度为0.5×10-6时相场演化示意图(Ⅰ型)
图 8 浓度为0时相场演化示意图(Ⅰ型)
图 9 位移-载荷曲线

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

3.2 准静态Ⅱ型裂纹问题

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

图 10 浓度为0.5×10-6时浓度扩散示意图(Ⅱ型)
图 11 浓度为0.5×10-6相场演化示意图(Ⅱ型)
图 12 浓度为0时相场演化示意图(Ⅱ型)
图 13 不同浓度下的位移载荷曲线(Ⅱ型)

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

图 14 浓度为0.5×10-6时不同加载步计算结果(Ⅱ型)
4 结论

1) 本文引入应变能拉压分解方法, 改进了相场断裂氢脆模型。改进模型相较于改进前的模型, 不仅能模拟Ⅰ型氢脆断裂问题, 并且能有效模拟Ⅱ型氢脆断裂过程。将浓度为0.5×10-6的Ⅰ型断裂计算的结果与相关文献中的计算结果对比, 改进模型精度可靠, 且改进模型能精确模拟脆弹性材料的瞬断过程。

2) 数值计算结果表明: 氢离子浓度越高, 断裂发生越快, 临界破坏载荷越小; 同时可以看出: 在断裂过程中, 氢离子会发生显著的聚集现象即向着应力集中区域聚集, 这均符合氢脆过程的特征。

3) 通过加载步对计算结果的影响分析可见, 分步计算时加载步对计算结果有着一定的影响: 加载步越大, 计算结果的稳定性和精度较差, 加载步愈小, 计算结果趋向于稳定值, 愈精确。

参考文献
[1] 刘德林, 陶春虎, 刘昌奎, 等. 钢氢脆失效的新现象与新认识[J]. 失效分析与预防, 2015, 10(6): 376-383.
LIU Delin, TAO Chunhu, LIU Changkui, et al. New phenomenons and knowledge of steel hydrogen embrittlement[J]. Failure Analysis and Prevention, 2015, 10(6): 376-383. (in Chinese) DOI:10.3969/j.issn.1673-6214.2015.06.009
[2] BELYTSCHKO T, BLACK T. Elastic crack growth in elements with minimal remeshing[J]. International Journal for Numerical Method in Engineering, 1999, 45: 601-620. DOI:10.1002/(SICI)1097-0207(19990620)45:5<601::AID-NME598>3.0.CO;2-S
[3] YANG L, YANG Y, ZHENG H. A phase field numerical manifold method for crack propagation in quasi-brittle materials[J]. Engineering Fracture Mechanics, 2021, 241: 107427. DOI:10.1016/j.engfracmech.2020.107427
[4] FRANCFORT G A, MARIGO J J. Revisiting brittle fracture as an energy minimization problem[J]. Journal of the Mechanics and Physics of Solids, 1998, 46(8): 1319-1342. DOI:10.1016/S0022-5096(98)00034-9
[5] MARTÍNEZ-PAÑEDA E, GOLAHMAR A, NIORDSON C F. A phase field formulation for hydrogen assisted cracking[J]. Computer Methods in Applied Mechanics and Engineering, 2018, 342: 742-761. DOI:10.1016/j.cma.2018.07.021
[6] WU J Y, MANDAL T K, NGUYEN V P. A phase-field regularized cohesive zone model for hydrogen assisted cracking[J]. Computer Methods in Applied Mechanics and Engineering, 2020, 358: 112614. DOI:10.1016/j.cma.2019.112614
[7] 吴建营. 固体结构损伤破坏统一相场理论、算法和应用[J]. 力学学报, 2021, 53(2): 301-329.
WU Jianying. On the theoretical and numerical aspects of the unified phase-field theory for damage and failure in solids and structures[J]. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(2): 301-329. (in Chinese)
[8] MANDAL T K, NGUYEN V P, WU J Y. Comparative study of phase-field damage models for hydrogen assisted cracking[J]. Theoretical and Applied Fracture Mechanics, 2021, 111: 102840. DOI:10.1016/j.tafmec.2020.102840
[9] MIEHE C, WELSCHINGER F, HOFACKER M. Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations[J]. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1273-1311. DOI:10.1002/nme.2861
[10] PENG F, HUANG W, MA Y E, et al. Phase field modeling of brittle fracture based on the cell-based smooth FEM by considering spectral decomposition[J]. International Journal of Computational Methods, 2020, 18(2): 2050016.
[11] MSEKH M A, SARGADO J M, JAMSHIDIAN M, et al. Abaqus implementation of phase-field model for brittle fracture[J]. Computational Materials Science, 2015, 96: 472-484. DOI:10.1016/j.commatsci.2014.05.071
[12] LIU G, LI Q, MSEKH M A, et al. Abaqus implementation of monolithic and staggered schemes for quasi-static and dynamic fracture phase-field model[J]. Computational Materials Science, 2016, 121: 35-47.
Simulating hydrogen embrittlement fracture based on phase field method
CHEN Pengcheng, MA Yu'e, PENG Fan, ZHOU Linglong     
School of Aeronautics, Northwestern Polytechnical University, Xi'an 710072, China
Abstract: The phase field hydrogen embrittlement fracture model is improved by introducing tension-compression split of strain energy. The numerical formulas of the model are provided, besides, the coupling term of concentration field and displacement field is deduced. The matlab software is used to compile the numerical program of phase field hydrogen embrittlement fracture. The modes I and II cracks of hydrogen embrittlement are simulated respectively. The simulation results show that hydrogen ions concentrate at the crack tip where stress concentration happens, and that the hydrogen concentration reduces the critical failure load of the square plate. Compared with the numerical results of the existing models, the improved model can accurately calculate the critical failure load in the mode I crack and capture the embrittlement fracture phenomenon when the phase field and the concentration field are accumulated near the crack tip. Moreover, the improved model can effectively simulate the mode II crack with hydrogen embrittlement.
Keywords: phase field method    multi-field coupling    hydrogen embrittlement    tension-compression split    brittle fracture    
西北工业大学主办。
0

文章信息

陈鹏程, 马玉娥, 彭帆, 周玲珑
CHEN Pengcheng, MA Yu'e, PENG Fan, ZHOU Linglong
基于相场法的氢脆断裂模拟
Simulating hydrogen embrittlement fracture based on phase field method
西北工业大学学报, 2022, 40(3): 504-511.
Journal of Northwestern Polytechnical University, 2022, 40(3): 504-511.

文章历史

收稿日期: 2021-08-12

相关文章

工作空间