基于相場法的氫脆斷裂模擬

陳鵬程, 馬玉娥, 彭帆, 周玲瓏

(西北工業大學 航空學院, 陜西 西安 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的表達式為

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,較為合理。

圖10 濃度為0.5×10-6時濃度擴散示意圖(Ⅱ型)

圖13 不同濃度下的位移載荷曲線(Ⅱ型)

圖14 濃度為0.5×10-6時不同加載步計算結果(Ⅱ型)

4 結 論

1) 本文引入應變能拉壓分解方法,改進了相場斷裂氫脆模型。改進模型相較于改進前的模型,不僅能模擬Ⅰ型氫脆斷裂問題,并且能有效模擬Ⅱ型氫脆斷裂過程。將濃度為0.5×10-6的Ⅰ型斷裂計算的結果與相關文獻中的計算結果對比,改進模型精度可靠,且改進模型能精確模擬脆彈性材料的瞬斷過程。

2) 數值計算結果表明:氫離子濃度越高,斷裂發生越快,臨界破壞載荷越小;同時可以看出:在斷裂過程中,氫離子會發生顯著的聚集現象即向著應力集中區域聚集,這均符合氫脆過程的特征。

3) 通過加載步對計算結果的影響分析可見,分步計算時加載步對計算結果有著一定的影響:加載步越大,計算結果的穩定性和精度較差,加載步愈小,計算結果趨向于穩定值,愈精確。