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

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072)

腐蝕化學(xué)環(huán)境會(huì)誘發(fā)金屬結(jié)構(gòu)產(chǎn)生氫脆(HE,hydrogen embrittlement),在與力學(xué)環(huán)境共同作用下,會(huì)更容易造成結(jié)構(gòu)的失效和破壞。相比于單純的斷裂破壞,氫脆斷裂在斷裂前的征兆不易觀察,臨界破壞載荷更小,裂紋擴(kuò)展速度更快,因此帶來(lái)的破壞性更強(qiáng)[1]。氫脆斷裂造成的工程失效案例屢見(jiàn)不鮮,涉及航空、航天、儲(chǔ)能等重大戰(zhàn)略性領(lǐng)域。因此研究氫脆斷裂問(wèn)題有重要的實(shí)際應(yīng)用價(jià)值。

對(duì)于單純的斷裂問(wèn)題,工程上一般采用有限元仿真來(lái)實(shí)現(xiàn)對(duì)斷裂過(guò)程的模擬和預(yù)測(cè),但是傳統(tǒng)的有限元方法(FEM)在處理斷裂問(wèn)題時(shí)需在裂紋擴(kuò)展過(guò)程中不斷追蹤裂紋界面,反復(fù)對(duì)裂尖進(jìn)行網(wǎng)格重新細(xì)分才能精確計(jì)算裂尖奇異場(chǎng)[2],因而在處理斷裂問(wèn)題時(shí)存在計(jì)算復(fù)雜、精度過(guò)低等缺點(diǎn)。

變分?jǐn)嗔严鄨?chǎng)法則可以很好地彌補(bǔ)FEM這些缺點(diǎn):斷裂相場(chǎng)法通過(guò)序參量的自動(dòng)演化來(lái)確定裂紋面,因此可以避免追蹤裂紋面的繁瑣過(guò)程。此外,相場(chǎng)法還可以在不引入額外的失效準(zhǔn)則條件下確定裂紋啟裂和擴(kuò)展的方向[3-4]。

與普通斷裂不同,氫脆斷裂伴隨復(fù)雜的氫擴(kuò)散過(guò)程,需考慮氫濃度場(chǎng)演化對(duì)材料斷裂韌性的影響。值得注意的是,相場(chǎng)斷裂模型的實(shí)質(zhì)就是相場(chǎng)和位移場(chǎng)相互耦合的偏微分方程組,不同場(chǎng)間的影響可以通過(guò)相應(yīng)的控制方程體現(xiàn),因此利用相場(chǎng)法來(lái)處理多場(chǎng)耦合下的斷裂問(wèn)題十分有效。

目前有關(guān)相場(chǎng)氫脆斷裂的主要研究有:利用基本相場(chǎng)斷裂模型,Martínez等[5]提出了相場(chǎng)氫脆斷裂的框架;Wu等[6]基于統(tǒng)一相場(chǎng)理論[7]將相場(chǎng)正則化內(nèi)聚裂縫模型應(yīng)用于氫脆斷裂的模擬;而Mandal等[8]則是對(duì)比研究了不同相場(chǎng)模型下的氫脆斷裂計(jì)算結(jié)果,深入研究了相關(guān)參數(shù)對(duì)于氫脆斷裂計(jì)算結(jié)果的影響。

在基本的相場(chǎng)斷裂模型[5,9-11]中,受拉受壓均會(huì)造成能量的釋放,這在計(jì)算中會(huì)造成一定誤差。因此其應(yīng)用局限于單純拉伸載荷作用下的斷裂模式。在模擬更一般的斷裂時(shí),需要分別考慮拉壓應(yīng)力分量的影響。本文則基于Miehe等[9]提出的應(yīng)變能拉壓譜分解方法,改進(jìn)了相場(chǎng)氫脆斷裂模型,并通過(guò)數(shù)值計(jì)算驗(yàn)證了改進(jìn)模型的可靠性。

1 基本模型

1.1 相場(chǎng)斷裂基本模型及控制方程

相場(chǎng)法通過(guò)引入一個(gè)輔助標(biāo)量場(chǎng)——相場(chǎng)φ來(lái)表征裂紋拓?fù)鋄9]。在一維情況下,φ的表達(dá)式為

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

(2)

由此可得正則化表面能

(3)

式中，Gc(θ)為臨界能量釋放率。引入應(yīng)力退化函數(shù)g(φ)=(1-φ)2+k,其中k是為了避免計(jì)算奇異而引入的附加參數(shù),在具體計(jì)算中盡量取足夠小,本文中均取k=1×10-6。則原應(yīng)變能ψe(ε)在考慮相場(chǎng)影響下有

(4)

總勢(shì)能泛函則為

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

(5)

在表面力f和體積力b作用以及邊界約束下,可以推導(dǎo)準(zhǔn)靜態(tài)下相場(chǎng)斷裂位移場(chǎng)和相場(chǎng)控制方程

(6)

式中：σ是應(yīng)力張量;n是表面法向量。

1.2 氫離子濃度擴(kuò)散方程

氫脆斷裂發(fā)生時(shí),氫離子會(huì)向材料內(nèi)部應(yīng)力集中部位擴(kuò)散聚集,大大降低構(gòu)件的斷裂韌性。根據(jù)相關(guān)的研究[5],表面氫離子濃度θ對(duì)材料斷裂韌性的影響可用(7)式表征

(7)

Gc(0)為氫離子濃度為0時(shí)材料的臨界能量釋放率。χ是損傷系數(shù),代表氫離子對(duì)斷裂韌性影響。式中表面氫離子濃度θ可由體積氫離子濃度C計(jì)算

(8)

(9)

式中：J是氫離子濃度通量；q是表面氫通量。J為關(guān)于應(yīng)力梯度和氫離子濃度的函數(shù),可以采用如下公式定義

(10)

圖2 濃度擴(kuò)散示意圖

2 改進(jìn)后相場(chǎng)氫脆斷裂模型

2.1 應(yīng)變能拉壓分解

在上述基本相場(chǎng)氫脆斷裂模型中,相場(chǎng)引起的應(yīng)力退化會(huì)引起應(yīng)力張量所有的分量發(fā)生改變,包括拉伸分量和壓縮分量,這與實(shí)際的斷裂機(jī)制有出入。為了將此模型推廣,本文引入Miehe等[9]提出的應(yīng)變能譜分解方法,重新考慮相場(chǎng)對(duì)應(yīng)變能的折減方式。首先對(duì)應(yīng)變張量進(jìn)行譜分解

(11)

式中，〈·〉為麥考利符號(hào)。相應(yīng)地,對(duì)彈性應(yīng)變能進(jìn)行拉壓分解,這樣相場(chǎng)折減只對(duì)正應(yīng)變能有效

(12)

(12)式對(duì)應(yīng)變張量求偏導(dǎo),得到譜分解下的應(yīng)力張量

(13)

λ和μ是拉梅常數(shù),I是二階單位張量。則彈性張量為

(14)

2.2 控制方程的離散化

(15)

N是形函數(shù)矩陣,B是其對(duì)應(yīng)的偏導(dǎo)矩陣。

2.3 求解策略

利用Newton-Raphson算法結(jié)合分步式算法,對(duì)(15)式中的非線性方程組求解

(16)

式中

(17)

時(shí)間導(dǎo)數(shù)項(xiàng)對(duì)應(yīng)的剛度矩陣M的表達(dá)式為

(18)

采用分步算法,可以完成相場(chǎng)和其余兩場(chǎng)的解耦。同時(shí),本文為了提高精確性,在計(jì)算過(guò)程中考慮濃度場(chǎng)和位移場(chǎng)的耦合項(xiàng)剛度矩陣,其推導(dǎo)過(guò)程如下

(19)

由剛度矩陣的表達(dá)式可以看出,在計(jì)算過(guò)程中需要確定積分點(diǎn)的靜水壓力值以及其梯度值。本文采用同樣形式的形函數(shù)Ni來(lái)計(jì)算靜水壓力

(20)

式中,σH,i是節(jié)點(diǎn)i處的靜水壓力,其表達(dá)式為

(21)

d是單元節(jié)點(diǎn)個(gè)數(shù),則有

(22)

利用Voigt標(biāo)記法將四階彈性張量Ds寫(xiě)成矩陣形式,則平面應(yīng)變問(wèn)題下的耦合項(xiàng)剛度矩陣公式為

(23)

式中，L的表達(dá)式為

L=

(24)

為了驗(yàn)證改進(jìn)方法的有效性,基于Matlab程序開(kāi)發(fā)平臺(tái),編寫(xiě)了多場(chǎng)耦合斷裂的有限元求解程序,并設(shè)計(jì)算例進(jìn)行驗(yàn)證。

3 數(shù)值算例

為了驗(yàn)證改進(jìn)模型的可行性以及對(duì)Ⅰ型氫脆斷裂的模擬情況，選取含單條裂紋的正方形板作為研究對(duì)象,分別施加拉伸位移載荷和剪切位移載荷,如圖3～4所示。試樣的邊長(zhǎng)為1 mm,預(yù)制裂紋長(zhǎng)度均為0.5 mm,位于試樣的正中間,底部為固定位移邊界。

圖3 單裂縫正方形板 圖4 單裂縫正方形板受拉伸位移載荷 受剪切位移載荷

3.1 準(zhǔn)靜態(tài)Ⅰ型裂紋問(wèn)題

在單純的拉伸位移載荷下,整個(gè)區(qū)域內(nèi)設(shè)置固定氫離子濃度條件C,分別取4種不同的濃度邊界條件C=0,0.1×10-6,0.5×10-6,1×10-6來(lái)研究氫離子濃度對(duì)斷裂過(guò)程的影響。圖5對(duì)比了不同濃度下的位移-載荷曲線,可知斷裂過(guò)程對(duì)氫離子濃度十分敏感:隨著濃度的增加,臨界破壞載荷不斷減小。以濃度為0.5×10-6時(shí)為例,如圖6所示,可以看出在斷裂過(guò)程中,氫離子不斷向應(yīng)力集中處擴(kuò)散;在相同位移加載下,濃度越大則相場(chǎng)演化越快,如圖7和圖8所示。同時(shí)為了驗(yàn)證本文模型的準(zhǔn)確性,選取了濃度為0.5×10-6時(shí)的計(jì)算結(jié)果與Martínez等[5]的研究結(jié)果進(jìn)行對(duì)比,如圖9a)所示。結(jié)果表明,改進(jìn)模型所計(jì)算臨界破壞載荷與原模型相差僅為3.954%,因此,改進(jìn)后模型能滿足計(jì)算精度;同時(shí)可以看出,改進(jìn)后模型能捕捉脆彈性材料的瞬斷過(guò)程。

圖5 不同氫濃度下的位移-載荷曲線

此外,為了研究位移加載步對(duì)計(jì)算結(jié)果的影響,在C=0.5×10-6的濃度邊界條件下,分別選取了Δu為8×10-5,4×10-5,2×10-5,1×10-5,5×10-6mm 5種不同位移加載步值,計(jì)算了相應(yīng)加載步下的破壞過(guò)程曲線,如圖9b)所示。在本文中采用分步式算法求氫脆相場(chǎng)斷裂模型時(shí),位移場(chǎng)和濃度場(chǎng)完全耦合,它們和相場(chǎng)的求解在前后相鄰的加載步間進(jìn)行。因此,大的加載步將影響求解的精度。如圖9b)所示,當(dāng)加載步逐漸減小,計(jì)算結(jié)果趨近于穩(wěn)定,但是計(jì)算時(shí)間成本不斷增加。因此在選取加載步時(shí)需要權(quán)衡計(jì)算精度和計(jì)算效率,而本算例所選取的加載步,其對(duì)應(yīng)的計(jì)算結(jié)果如圖9b)紅色粗實(shí)線所示,綜合考慮其計(jì)算精度和計(jì)算效率,是個(gè)較理想的選擇。

圖6 濃度為0.5×10-6時(shí)濃度擴(kuò)散示意圖(Ⅰ型)

圖7 濃度為0.5×10-6時(shí)相場(chǎng)演化示意圖(Ⅰ型)

圖8 濃度為0時(shí)相場(chǎng)演化示意圖(Ⅰ型)

圖9 位移-載荷曲線

3.2 準(zhǔn)靜態(tài)Ⅱ型裂紋問(wèn)題

引入拉壓分解后,可以更精確地模擬Ⅱ型斷裂過(guò)程。本例中,在剪切位移載荷下,分別取2種不同的濃度條件C=0,0.5×10-6來(lái)考慮氫離子濃度對(duì)斷裂過(guò)程的影響。以0.5×10-6為例,濃度擴(kuò)散和斷裂(相場(chǎng))演化的過(guò)程分別如圖10～11所示。由圖可見(jiàn),引入拉壓分解后的模型可以有效地模擬Ⅱ型氫脆斷裂過(guò)程,相較于原模型,采取應(yīng)變能拉壓分解后的模型不會(huì)出現(xiàn)裂紋的偽分叉[12]。同時(shí),圖12給出了濃度為0×10-6同位移載荷下的相場(chǎng)演化模擬圖。結(jié)合圖13不同濃度下斷裂過(guò)程的位移-載荷曲線,可見(jiàn)在Ⅱ型斷裂模式下,整個(gè)斷裂過(guò)程對(duì)氫離子濃度同樣很敏感。

本節(jié)同樣進(jìn)行了位移加載增量步對(duì)Ⅱ型相場(chǎng)氫脆斷裂計(jì)算結(jié)果的影響分析。如圖14所示,加載步愈小,最終計(jì)算結(jié)果愈精確?？梢钥闯?在本算例中,綜合計(jì)算精度和計(jì)算效率考慮后所選取的加載步Δu=1×10-5mm,較為合理。

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

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

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

4 結(jié) 論

1) 本文引入應(yīng)變能拉壓分解方法,改進(jìn)了相場(chǎng)斷裂氫脆模型。改進(jìn)模型相較于改進(jìn)前的模型,不僅能模擬Ⅰ型氫脆斷裂問(wèn)題,并且能有效模擬Ⅱ型氫脆斷裂過(guò)程。將濃度為0.5×10-6的Ⅰ型斷裂計(jì)算的結(jié)果與相關(guān)文獻(xiàn)中的計(jì)算結(jié)果對(duì)比,改進(jìn)模型精度可靠,且改進(jìn)模型能精確模擬脆彈性材料的瞬斷過(guò)程。

2) 數(shù)值計(jì)算結(jié)果表明:氫離子濃度越高,斷裂發(fā)生越快,臨界破壞載荷越小;同時(shí)可以看出:在斷裂過(guò)程中,氫離子會(huì)發(fā)生顯著的聚集現(xiàn)象即向著應(yīng)力集中區(qū)域聚集,這均符合氫脆過(guò)程的特征。

3) 通過(guò)加載步對(duì)計(jì)算結(jié)果的影響分析可見(jiàn),分步計(jì)算時(shí)加載步對(duì)計(jì)算結(jié)果有著一定的影響:加載步越大,計(jì)算結(jié)果的穩(wěn)定性和精度較差,加載步愈小,計(jì)算結(jié)果趨向于穩(wěn)定值,愈精確。