0 引言
本文提出将PCE与氢脆断裂相场法相结合,建立氢脆的随机断裂相场法理论。输出响应通过PCE建模,并由此构建了完整的理论与数值分析框架,用于分析断裂过程中几种材料因素对破坏的随机化影响,以期为金属受氢脆影响发生断裂破坏提供准确的预测与防控手段。
1 多项式混沌展开理论
1.1 多项式混沌展开及Legendre多项式
系统的随机响应可以通过多项式混沌展开(PCE)表示为均方收敛级数。假设随机输入变量ξ=(ξ1,ξ2,…,ξm)独立且均匀分布在[-1,1],且具有有界二阶矩(即均值和标准差有界)。PCE展开的输出U可以表示为
$U=U\left(\xi \right)=\sum _{i=0}^{\infty }{u}_{i}{\Psi }_{i}\left(\xi \right)$
其中,Ψi(ξ)是多维Legendre正交基函数,ui是展开系数。为简化计算,通常将展开截断至阶数P,并表示为
$\hat{U}=\sum _{i=0}^{P-1}{u}_{i}{\Psi }_{i}\left(\xi \right),\xi =({\xi }_{1},{\xi }_{2},\dots,{\xi }_{m})$
展开的项数P与随机变量个数m及多项式阶数p相关,通常可以通过下式进行计算为
$P=\mathbb{C}_{m+p}^{p}$
Legendre多项式系数{Ln(ξ)${\}}_{n=0}^{\infty }$关于权函数ρ(ξ)=1正交,其中ξ∈[-1,1],并满足正交关系为
${\int }_{-1}^{1}{L}_{i}\left(\xi \right){L}_{j}\left(\xi \right)dx=\frac{2}{2i+1}{\delta }_{ij}$
其中,δij表示 Kronecker函数,其取值为
δij=$\left\{\begin{array}{l}1,i=j\\ 0,i\ne j\end{array}\right.$
Legendre多项式可以通过微分方程求解
$\frac{d}{dx}\left\{(1-{\xi }^{2})\frac{d}{d\xi }{L}_{k}\left\{\xi \right\}\right\}$+k(k+1)Lk{ξ}=0
其解可以通过Rodrigues公式[18]表示为
Lk{ξ}=$\frac{1}{{2}^{k}k!}\frac{d}{d{\xi }^{k}}${(1-ξ2)k}
Legendre多项式的三项递推关系为
(k+1)Lk+1(ξ)=(2k+1)ξLk(ξ)-kLk-1(ξ)
此外,前两阶Legendre多项式为
$\begin{array}{l}{L}_{0}\left(\xi \right)=1\\ {L}_{1}\left(\xi \right)=\xi \end{array}$
根据上述递推关系,可以得到前七阶的一维Legendre多项式如下。
1.2 模型收敛性的判断
通常,预测值$\hat{U}$与精确解U之间存在一定误差。为了量化误差并评估拟合效果,我们引入残差ε,其定义为
$\epsilon =U\left(x\right)-\sum _{i=p}^{\infty }{u}_{i}{\Psi }_{i}\left(\xi \right)$
其中,ξ=(ξ1,ξ2,…,ξm);ui是未知系数;注意到,残差的绝对值越小,代理模型的准确性越高。因此,目标是最小化残差平方和的期望,确定多项式混沌系数向量u为
$\hat{u}$=ArgminE∑(U-$\hat{U}$)2
表1 前七阶一维Legendre多项式Table 1 The first seven one-dimensional Legendre polynomials |
| k | Legendre polynomials:Lk(ξ):ξ∈[-1,1] |
|---|---|
| 0 | 1 |
| 1 | ξ |
| 2 | $\frac{1}{2}$(3ξ2-1) |
| 3 | $\frac{1}{2}$(5ξ3-3ξ) |
| 4 | $\frac{1}{8}$(35ξ4-30ξ2+3) |
| 5 | $\frac{1}{8}$(63ξ5-70ξ3+15ξ) |
| 6 | $\frac{1}{16}$(231ξ6-315ξ4+105ξ2-5) |
| 7 | $\frac{1}{16}$(429ξ7-693ξ5+325ξ3-35ξ) |
通过在标准正态空间中选择N个回归点(ξ1,ξ2,…,ξN)作为输入样本点,多项式混沌展开系数u可以由下式确定为
$\hat{u}$=Argmin$\frac{1}{N}$∑(U-$\hat{U}$)2
最终,多项式混沌展开系数$\hat{u}$=(ΨTΨ)-1ΨTU;代理模型的解可由式公式(2)与多项式混沌展开系数u计算得到。
为评估代理模型的准确性并判断多项式混沌展开的收敛性,使用均方根误差(RMSE)的变异系数CV(RMSE)来衡量[19]。计算公式为
$CV\left(RMSE\right)=\sqrt[ ]{\frac{\sum _{i=1}^{n}({U}_{i}-{\hat{U}}_{i}{)}^{2}}{n}}/\frac{\sum _{i=1}^{n}{U}_{i}}{n}$
其中,n为样本数,Ui表示氢脆相场断裂模型的输出,${\hat{U}}_{i}$是PCE代理模型的输出。当CV(RMSE)小于5%时,认为模型已收敛且拟合精度较高。简写为CV。
2 耦合氢作用的随机相场断裂理论
2.1 相场断裂理论
首先,介绍氢脆的相场模型框架。假设存在一个具有封闭边界δΩ的n维域Ω,其外法向量为 $\overrightarrow{n}$,边界δΩ分为牵引边界δΩt和位移边界δΩu,其中分别施加力边界条件t*(x) 和位移边界条件u*(x)。固体中的裂缝通过裂缝相场d(x)描述,d(x)∈[0,1],其中 d(x)=0表示材料完好,d(x)=1表示形成宏观裂缝。裂缝带的边界为δB,外法线为$\overrightarrow{{n}_{B}}$,如图1所示。
根据Francfort和Marigo提出的裂隙变分框架[20],固体的总能量Π(u)可由应变能Πe(u)、裂隙能量 Πc(d) 和外部力Πext(u) 表示,如下所示为
$ \begin{aligned} \Pi(u) & =\Pi_{e}(\boldsymbol{u})+\Pi_{c}(d)-\Pi_{\text {ext }}(\boldsymbol{u}) \\ & =\int_{\Omega \backslash \Gamma} \psi_{e}(\varepsilon(\boldsymbol{u})) d V+\int_{\Gamma} G_{c} d \Gamma \\ & -\int_{\Omega} \boldsymbol{b} \cdot \boldsymbol{u} d V-\int_{\partial \Omega_{t}} \boldsymbol{t} \cdot \boldsymbol{u} d A \end{aligned} $
其中,ψe(ε(u))表示单位体积的弹性应变能,Gc为能量释放率,b为区域Ω内的体积力,t为边界δΩ上的表面牵引力。将上述总能量正则化为
$\begin{aligned} \Pi(\boldsymbol{u}, d)= & \int_{\Omega} \psi_{e}(\boldsymbol{\varepsilon}(\boldsymbol{u}), d) d V+\int_{\Omega} \boldsymbol{G}_{c} \gamma(d, \nabla d) d V \\ & -\int_{\Omega} \boldsymbol{b} \cdot \boldsymbol{u} d V-\int_{\partial \Omega_{t}} \boldsymbol{t} \cdot \boldsymbol{u} d A \end{aligned}$
其中,γ(d,∇d)为裂缝表面能量密度函数,表示为
γ(d,∇d)=$\frac{1}{{c}_{0}}\left[\frac{1}{l}\alpha \left(d\right)+l{\left|∇ d\right|}^{2}\right]$
其中,c0是与能量密度函数有关的标量参数,l是相场长度尺度。由于损伤的影响,固体的弹性应变能ψe(ε(u))会发生降低,可以表示为
ψe(ε(u),d)=ω(d)ψ0
其中,ω(d)是表示应变能退化程度的退化函数,ψ0是材料完好时的应变能。将ψ0分解为正负部分
ψ0=${\psi }_{0}^{+}$+${\psi }_{0}^{-}$
$\psi_{0}^{ \pm}=\frac{1}{2 \bar{E}_{0}}\left\langle\overline{\boldsymbol{\sigma}}_{1}\right\rangle_{ \pm}^{2}$
其中,<.>±是麦考利括号,有<a>±=<a±|a|>/2;拉伸模量${\overline{E}}_{0}$=λ0+2μ0。因此,分解后的应变能函数可重新表示为
ψe(ε(u),d)=ω(d)${\psi }_{0}^{+}$+${\psi }_{0}^{-}$
得到总能量后,位移场和相场(u,d)通过最小化总能量来确定为
(u,d)=Arg{InfΠ(u,d)},$\stackrel{ ·}{d}$≥0,d∈(0,1)
总能量相对于位移场和相场的变化的控制方程可以写为
$\left\{\begin{array}{l} ∇ ·\sigma +b=0,\\ f-{G}_{c}{\delta }_{d}\gamma =0\end{array}\right.$
$\left\{\begin{array}{l}\sigma ·n=t on\partial {\Omega }_{t}\\ u={u}^{*} on\partial {\Omega }_{u}\\ ∇ d·n=0 on\partial \Omega \end{array}\right.$
其中,能量驱动力f与裂纹密度泛函为
f=-$\frac{\partial \psi }{\partial d}$=-ω'(d)${\psi }_{0}^{+}$
δdγ=$\frac{1}{{c}_{0}}\left[\frac{1}{l}\alpha \text{'}\left(d\right)+2l\Delta d\right]$
为了保证裂纹演化的不可逆性,引入有效裂纹驱动力$\overline{Y}$代替正应变能${\psi }_{0}^{+}$。随着裂纹的发展,$\overline{Y}$也由其最大值替代。有效裂缝驱动力的初始阈值H0为
H0=$\frac{{f}_{t}^{2}}{2{\overline{E}}_{0}}$
2.2 相场特征方程
根据已有的研究[23],相场的特征方程表示为
α(d)=2d-d2,c0=π
ω(d)=$\frac{{(1-d)}^{p}}{{(1-d)}^{p}+{a}_{1}d·(1+{a}_{2}d+{a}_{3}{d}^{2})}$
其中,p≥2,对于脆性断裂的线性软化规律,可以采用以下参数
p=2,a1=2,a2=-0.5,a3=0
2.3 氢浓度场的影响
在区域Ω中,氢的质量扩散过程(由浓度场C描述)受以下方程控制
$\left\{\begin{array}{ll}\stackrel{ ·}{C}+∇ ·J=0& in\Omega \\ J·n={J}^{*}& on\partial {\Omega }_{J}\end{array}\right.$
其中,氢通量J由修正后的菲克定律来确定[24]。
J=-D∇C+$\frac{D{V}_{H}}{RT}$C∇σH
其中,D为扩散系数,VH是偏摩尔体积,σH是静水压力,R代表气体常数,T为温度(300K)。氢浓度变化率为
$\stackrel{ ·}{C}$=DΔC-$\frac{D{V}_{H}}{RT}$∇C·∇σH-$\frac{D{V}_{H}}{RT}$CΔσH
此外,氢浓度C的表面浓度$\stackrel{~}{C}$可通过Langmuir-McLean等温线[25]获得
$\stackrel{~}{C}$=$\frac{C}{C+exp\left(-\frac{\Delta {g}_{b}^{0}}{RT}\right)}$
其中,C的单位是杂质摩尔分数,Δ${g}_{b}^{0}$是吉布斯自由能差。当体积浓度C较低时,以wt.ppm(百万分之一重量)为单位,以5.5×10-5C代替。氢降解效应函数ϕ($\stackrel{~}{C}$)表达式[26]为
ϕ($\stackrel{~}{C}$)=1-β$\stackrel{~}{C}$
氢影响下的失效强度ft和断裂能Gf表达为
ft($\stackrel{~}{C}$)=$\sqrt[ ]{\varphi \left(\stackrel{~}{C}\right)}$ft0,Gf($\stackrel{~}{C}$)=ϕ($\stackrel{~}{C}$)Gf0
其中,ft0、Gf0分别代表未受氢影响时的失效强度与断裂能。
2.4 非侵入式随机断裂相场法实现
如图2所示,氢脆的随机相场断裂模型的计算流程如下:首先,完成相场断裂模型的数值计算;然后,将计算结果作为输入数据集导入后续的多项式混沌展开(PCE)拟合中。这一过程无需对现有的有限元计算流程进行修改,因此具有广泛的适用性。在相场断裂的计算中,我们采用了分离式求解器和交错方法。在图示的三个组成部分中,PCE拟合后的计算结果被视为一个“黑箱”处理单元,其输入为不确定性的材料参数,输出则是代理模型得到的力-位移曲线及相场变量。对于随机空间中的每一个“黑点”,我们通过PCE拟合的模型来预测相应的力学响应。
3 模型验证与收敛性分析
本节介绍了采用多项式混沌展开(PCE)进行拟合的随机相场断裂模拟算例,并与确定性分析结果进行了对比,验证了PCE方法在预测随机断裂方面的优良性能。由于该方法为非侵入式算法,具有较好的适用性,能够方便地集成到现有计算流程中。在计算过程中,首先求解氢脆的相场断裂问题,并将得到的解作为样本输入进行PCE拟合。响应的随机性通过均值和标准差进行表征。我们还开发了一个方法,用于捕捉应力-应变曲线中的关键变化点,模拟应力随应变增长的过程,并提取该阶段的关键数据。接着,我们使用测试集评估PCE拟合模型在氢脆断裂预测中的表现。此外,记录了PCE代理模型的计算时间,以证明该方法在氢脆断裂问题中的应用潜力。
3.1 氢环境下单侧预制缝板的拉裂
表2 氢环境下单侧预制缝板的拉裂:材料参数Table 2 Tensile test of single-side notched plate in hydrogenation environment: Material parameters |
| Parameters | Values | Parameters | Values |
|---|---|---|---|
| ft | 2445.42MPa | Δ${g}_{b}^{0}$ | 30KJ/mol |
| Gc | 30000J/m2 | β | 0.89 |
| E0 | 210GPa | VH | 2×10-6m3/mol |
| v0 | 0.3 | D | 0.0127mm2/s |
我们首先采用确定性材料参数分析缺口板。在氢气作用下,裂纹周围氢浓度增加,导致材料失效强度和能量释放率下降,促进裂纹扩展并最终破坏。图4展示了不同时间下氢在表面的分布。
为评估不确定性材料参数对失效峰值荷载的影响,选择杨氏模量和断裂能作为随机输入变量,并假设其服从均匀分布。图6展示了最大、平均和最小值。结果表明,损伤沿预制裂纹发展,相场变量从0增加到1,标志着材料由完好到完全损坏。随着位移增大,反作用力上升,直到材料破坏,荷载下降。
值得注意的是,随机材料参数导致了不同的峰值应力响应。即使PCE阶数仅为2,构建的代理模型与氢脆断裂相场模拟的结果高度一致,验证了该算法的有效性和准确性。
图6显示了考虑材料参数不确定性后的典型物理力学特征。当随机输入变量为杨氏模量时,力-位移曲线的斜率发生变化,最大值对应最高峰值应力,且上升斜率大于其他情况。通过PCE拟合的代理模型很好地还原了这一特征。材料参数Gc的变化也显示了实际的力学影响:更大的Gc导致裂纹发生时间推迟,这一特征也能在代理模型结果中明显观察到。可以看出,即使PCE阶数较低,仍能得到与蒙特卡洛模拟相似的结果。为了验证这一点,我们在后续工作中计算了PCE与蒙特卡洛模拟结果的均值和标准差。
为说明不确定性材料参数对MCs结果的影响,选取了在相同荷载条件下具有不同材料参数的模拟结果。由于位移以恒定速度施加在边界上,结果选择在相同时间点进行对比。图7展示了在相同位移荷载下的不同裂纹扩展情况。
可以看出,较大的断裂能导致裂纹扩展较晚,这与力-位移曲线的观察一致。在图6中,具有不确定性材料参数的力-位移曲线显示,最小断裂能在位移Δu=2.5μm时发生开裂,而最大断裂能则延迟到Δu=3.3μm。通过PCE拟合的代理模型也能较好地还原这一特征。
在使用相同数据集时,高阶代理模型在提高拟合精度方面表现优异,但同时也会增加计算时间,因此需要在计算效率和模型精度之间找到平衡。图8显示了采用不同阶数时PCE拟合代理模型所需的时间及其对应的CV值。可以发现,随着阶数的增加,拟合时间显著增加,而与较低阶数结果相比,CV值有所减小。
通常,在确保误差可接受的前提下,会选择较低的阶数。为此,我们使用相对误差和CV值评估拟合效果。相对误差包括平均值和标准差的相对误差,标准差作为数据分散度的衡量指标,其公式为:
$\sigma =\sqrt[ ]{{\sigma }^{2}}=\sqrt[ ]{\frac{\sum _{i=1}^{n}({x}_{i}{-\overline{x})}^{2}}{n}}$
相关结果已整理在表3中。可以看出,即便采用较低阶数的PCE,拟合结果依然与蒙特卡洛模拟(MCs)高度一致,验证了其较高的准确性。此外,较低的阶数配置不仅保证了较短的计算时间,还提高了模拟的效率。
表3 氢环境下单侧预制缝板的拉裂:相对误差及CV值Table 3 Tensile test of single-side notched plate in hydrogenation environment: Relative errors and CV values |
| Method | Mean | Relative error(%) | Standard deviation | Relative error(%) | CV (%) |
|---|---|---|---|---|---|
| MCs | 8.96844 | — | 6.72350 | — | — |
| 2-order PCE | 8.96653 | 0.021 | 6.71921 | 0.064 | 0.2309 |
| 3-order PCE | 8.96715 | 0.014 | 6.72148 | 0.03 | 0.0998 |
3.2 预制缺陷板的裂纹扩展分析
材料参数见表4,同时考虑了均匀分布的不确定性材料参数。
表4 预制缺陷板的裂纹扩展分析:材料参数Table 4 Crack growth in prefabricated defective plates: Material parameters |
| Parameters | Values | Parameters | Values |
|---|---|---|---|
| ft | 1778.78MPa | Δ${g}_{b}^{0}$ | 30KJ/mol |
| Gc | 90000J/m2 | β | 0.89 |
| E0 | 200GPa | VH | 2×10-6m3/mol |
| v0 | 0.3 | D | 1×10-8mm2/s |
考虑了材料参数的不确定性,假设其服从均匀分布,随着加载位移增大,材料经历从完好到完全损坏的转变,反作用力逐渐增大,直至最大应力时失效,荷载下降。引入随机材料参数后,得到不同的最大应力值,揭示了材料失效过程中的随机性影响。类似前文,我们分析了不同阶数的PCE拟合代理模型。图11展示了2阶PCE拟合的结果,表明代理模型的预测结果与金属氢脆的相场模拟一致,验证了即使在较低阶数下,PCE代理模型仍能准确捕捉裂纹扩展的两段过程。
图12 预制缺陷板的裂纹扩展分析:随机材料输入参数不同阶段的裂纹扩展路径。(a) Gc取得最大值;(b) Gc取得平均值;(c) Gc取得最小值Fig.12 Crack growth in prefabricated defective plates:Crack propagation paths at different stages under stochastic input material parameters.(a) Gc max; (b) Gc mean; (c) Gc min |
表5 预制缺陷板的裂纹扩展分析:相对误差及CV值Table 5 Crack growth in prefabricated defective plates: Relative errors and CV values |
| Method | Mean | Relative error (%) | Standard deviation | Relative error (%) | CV (%) |
|---|---|---|---|---|---|
| MCs | 7.11479 | / | 5.52506 | / | / |
| 2-order PCE | 7.11747 | 0.0376 | 5.52905 | 0.0722 | 0.3608 |
| 3-order PCE | 7.11602 | 0.0173 | 5.52744 | 0.029 | 0.0363 |
4 结论
提出了一种非侵入式的多项式混沌展开(PCE)拟合氢脆随机相场断裂方法。通过结合PCE拟合后处理与有限元计算工具,构建了完整的理论与数值分析框架,旨在分析断裂过程中多种材料因素对破坏随机化的影响。
研究表明,即便在低阶PCE拟合下,训练后的PCE模型也能有效捕捉力-位移关系,并再现裂纹从起裂到扩展过程中的应力-应变曲线。通过两个经典算例验证了PCE拟合效果,进行了力-位移曲线的不确定性分析,并进一步探讨了PCE的拟合性能。即使在较低阶数的情况下,PCE拟合仍能得到与期望值和标准差相近的结果,证明了该方法在处理不确定性材料参数条件下氢脆问题的有效性。此外,所提出的PCE模型能够显著减少数值模拟的计算时间,并直观展示不确定性因素对结果的影响,从而有效节省计算成本。