一種改進的相場模型裂紋模擬方法

高 宇， 張晨露， 岳 強， 王 橋， 周 偉， 程勇剛， 李通盛

(1. 大唐宣威水電開發(fā)有限公司， 云南 宣威 655400； 2. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室， 湖北 武漢 430072)

在實際工程應(yīng)用中，對斷裂的數(shù)值計算有較高的要求，目前研究比較多的方法包括內(nèi)聚力模型(Cohesive Zone Model, CZM)[2]、連續(xù)斷裂力學(xué)法(Continuum Damage Mechanics, CDM)[3]、擴展有限元法(Extended Finite Element Method, XFEM)[4]。內(nèi)聚力模型能夠模擬材料中裂紋萌生、演化直至完全失效的連續(xù)過程[5]，但是由于需要事先通過內(nèi)聚力單元來定義斷裂路徑，因此內(nèi)聚力模型很難模擬復(fù)雜受力情況下和含有不同方向的多裂紋的裂紋擴展規(guī)律。連續(xù)斷裂力學(xué)法是另一種常用于有限元法中求解斷裂力學(xué)問題的方法，為了得到較高的計算精度并有效預(yù)測裂紋的擴展規(guī)律，該類型的方法需要在裂紋尖端處加密網(wǎng)格。擴展有限元法是傳統(tǒng)有限元法的改進，擴展有限元法雖然模擬裂紋擴展問題較為有效，但同樣需要將裂紋面處理為不連續(xù)面，在模擬多裂紋交叉時存在一定難度。除了以上研究較多的方法，還有許多其他方法被應(yīng)用于裂紋擴展問題中。例如Tsang等[6]采用離散裂隙網(wǎng)絡(luò)模型(Discrete Fracture Network Model, DFN)能夠綜合考慮裂紋中的斷裂力學(xué)和流體流動，但該方法計算精度較差，并且對裂紋的角度有較大的限制。Zhang等[7]采用了離散元法(Discrete Element Method, DEM)對散體結(jié)構(gòu)中的水力劈裂過程進行了模擬，該方法模擬真實問題時計算量過大，無法實時得到計算結(jié)果。

近年來，基于相場模型(Phase-Field Model, PFM)[8]的裂紋模擬方法得到了廣泛研究。傳統(tǒng)的斷裂力學(xué)中，存在著裂紋面和彈性體之間的不連續(xù)性問題，甚至還會有裂紋尖端的奇異性問題[9]。相場模型通過引入一個相場參數(shù)或函數(shù)來描述被近似為一個裂紋場(Crack Field)的裂紋面，將完全破壞的裂紋面與沒有破壞的彈性體之間光滑化，消除了不連續(xù)性[10]。該方法的優(yōu)勢在于可以非常方便地模擬復(fù)雜裂紋的擴展、貫通、交叉和分叉等行為，能夠預(yù)測裂紋開裂和擴展而不需要額外的判據(jù)[11]。一些不同的相場模型分別由不同的學(xué)者在物理界和力學(xué)界提出，Miehe等[12]提出了一種熱力學(xué)一致性相場模型來模擬斷裂問題，該模型基于連續(xù)力學(xué)和熱力學(xué)。本文應(yīng)用的相場模型主要是基于力學(xué)界中的脆性斷裂變分方程[13]和正則化方程[14]，該方程是Griffith斷裂理論[15]的延伸。對于相場模型的求解，目前大部分的方法均是基于有限元法。

從上文介紹中可以發(fā)現(xiàn)，相場模型是一種求解斷裂問題較為理想的方法，該方法可以非常方便地模擬多裂紋的擴展、貫通、開叉等規(guī)律，并相對容易由二維問題擴展到三維問題，是一種非常有前景的求解斷裂問題的方法。

1 準(zhǔn)靜態(tài)問題的相場模型

1.1 斷裂變分準(zhǔn)則

相場模型是基于Griffith斷裂理論的變分準(zhǔn)則，認(rèn)為彈性體的總勢能包括彈性體能和裂紋表面能[16]。

(1)

(2)

式中：Ω為彈性體的求解域；Γc?Ω為彈性體內(nèi)的裂隙面(圖1)；u為位移向量；ε為應(yīng)變張量；Gc為材料的斷裂韌度；ψ0為彈性體的能量密度；λ，μ為Lame常數(shù)。應(yīng)變張量和位移向量的關(guān)系為：

(3)

式中：?為梯度符號。

圖1 彈性體中的裂紋面

1.2 正則化方程

式(1)中包含對裂紋面的積分，在時間數(shù)值計算中存在兩大困難，一是位移場在裂紋處存在不連續(xù)性，二是從能量方程中無法直接得到最佳的裂紋擴展路徑。為了解決這一問題，一些學(xué)者提出了相應(yīng)的正則化方程[14]：

(4)

式中：s為用來表征裂紋場的場變量，其值在0～1之間變化，如果是0，表示完全破壞，如果是1，表示沒有破壞；參數(shù)l0>0用來控制s的過渡區(qū)域的寬度(圖2)；η為大于0且遠(yuǎn)遠(yuǎn)小于1的無量綱參數(shù)，用來模擬斷裂時，即s=0，產(chǎn)生的殘余剛度。

圖2 s表示的裂紋用相場

于是應(yīng)力張量可以定義為：

=(s2+η)(λtr(ε)I+2με)

(5)

1.3 控制方程

為了得到相場模型的最終控制方程，可以通過最小勢能原理建立的泛函推導(dǎo)出相場模型的弱積分形式：

(6)

(7)

式(6)(7)分別對應(yīng)彈性力學(xué)部分和相場部分，分別對式(6)左端域積分和式(7)的第一項域積分進行分部積分，可得：

(8)

(9)

式中：n是邊界?Ω上的外法相向量。由于式(8)對任意δu都成立，于是式(8)(9)對應(yīng)的控制方程分別為：

divσ=0 onΩ

(10)

(11)

相應(yīng)的邊界條件分別為：

(12)

(13)

式中：?Ωt為面力已知的邊界(圖3)；?Ωs′為相場模型的Neumann邊界(圖4)。

圖3 彈性部分的邊界條件

圖4 相場部分的邊界條件

1.4 歷史場變量

(14)

εi和ni分別為主應(yīng)變和主應(yīng)變方向，則

(15)

其中尖括號的含義為：

〈x〉±=(x± |x|)/2

(16)

式(5)便可以寫為：

(17)

相場方程可改寫為：

(18)

式中：?為歷史常變量，其定義為：

(19)

最終的控制方程和邊界條件可以歸納為：

(20)

彈性部分的邊界條件和普通無裂紋的彈性力學(xué)問題的邊界條件一樣(圖3)，?Ωu為位移已知邊界，相場部分的邊界條件為在整個求解域的邊界?Ω上為Neumann邊界條件，即在?Ω=?Ωs′上有sini=0，且在裂紋面Γc上s=0(圖4)。

2 分步解耦算法

在采用有限元法數(shù)值計算過程中，對于每一加載時間步，可以直接對式(20)進行聯(lián)立求解，為避免直接求解非線性方程組，本文采用解耦的方法分別求解彈性方程和相場方程。對于每一時間步，將式(20)解耦為彈性方程和相場方程：

(21)

(22)

式中：σn，un，sn，?n分別為第n時間步的應(yīng)力張量值、位移場、相場和歷史場變量。分步解耦算法的具體流程圖見圖5。

圖5 分步解耦算法的具體流程圖

3 應(yīng)變-相場法預(yù)設(shè)裂紋的改進方法

本文擬采用有限元法求解相場模型，位移場通過有限元法形函數(shù)近似可得

(23)

以二維情況為例，有

u=[u1(x)u2(x)]T

(24)

(25)

(26)

能量密度函數(shù)ψ0可以寫為：

(27)

當(dāng)模型有初始裂紋時，一般的方法是采用如下的應(yīng)變-相場法來預(yù)設(shè)初始裂紋[1]：

(28)

式中：L為裂紋曲線或裂紋面；d(x,L)為點x到L的最短距離；H0為初始?xì)v史場變量的值。

但參照已有的數(shù)值算例表明，采用該方法的計算結(jié)果與幾何不連續(xù)法在某些情況下差別較大，主要原因為幾何不連續(xù)法可以認(rèn)為是將裂紋面上下單元弱化，裂紋面上下單元完全分離。本文在此基礎(chǔ)上，提出了改進的應(yīng)變-相場法來預(yù)設(shè)裂紋，期望得到與幾何不連續(xù)法預(yù)設(shè)裂紋相近似的效果，其基本思想是假設(shè)在裂紋面附近的單元一直弱化，無論其應(yīng)變狀態(tài)如何。

其中，對于平面應(yīng)變問題有

(29)

在改進的應(yīng)變 - 相場法中，首先通過式(28)來計算裂紋面附近積分點的初始?xì)v史場變量，且在所有計算過程中，采用下式來計算弱化后的D矩陣

(30)

于是式(21)的弱積分形式可以寫為

(31)

最終可得的彈性方程的有限元離散格式為：

(32)

(33)

(34)

4 數(shù)值算例

4.1 含單邊裂紋的受拉方板

考慮如圖6所示包含初始裂紋的受拉方板，邊長為1 mm，劃分為如圖7所示的54040個三角形單元，在預(yù)計裂紋擴展路徑處局部加密，該處網(wǎng)格大小約為0.001 mm。設(shè)該問題為平面應(yīng)變問題，彈性模量為210 kN/mm2，泊松比為0.3，斷裂韌度為Gc=2.7×10-3kN/mm2。采用位移加載法，前500步位移增量為Δu=1×10-5mm，之后位移增量調(diào)整為Δu=1×10-6mm直到完全破壞。在計算過程中，取l0=0.004 mm，η=0。

首先采用幾何不連續(xù)法來預(yù)設(shè)初始裂紋，即在裂紋處將上下單元在幾何上分離，計算得到的裂紋擴展過程如圖8所示。

圖6 受拉方板/mm

圖7 受拉方板網(wǎng)格劃分

圖8 受拉方板裂紋擴展路徑(裂紋幾何建模)

若采用應(yīng)變 - 相場法來預(yù)設(shè)初始裂紋，得到的裂紋擴展路徑如圖9所示。可以發(fā)現(xiàn)，采用該方法來預(yù)設(shè)裂紋時，模型會早一些開裂，主要原因在于對于該模型，采用應(yīng)變 - 相場法預(yù)設(shè)裂紋時，把裂紋附近在l0/2范圍內(nèi)的單元均弱化了，而采用幾何不連續(xù)來預(yù)設(shè)裂紋時，并未弱化任何單元，或者說只將包含裂紋面的單元弱化了，因此對于該模型而言，采用應(yīng)變 - 相場法來預(yù)設(shè)裂紋時，更多的單元被弱化，模型會更早開裂。進一步采用改進的應(yīng)變 - 相場法來預(yù)設(shè)裂紋，其得到的裂紋擴展路徑如圖10所示。通過比照受拉方板的力 - 位移曲線(圖11)，可以發(fā)現(xiàn)采用改進的方法會進一步提前裂紋的開裂，但該計算結(jié)果是合理的，與采用幾何不連續(xù)來建立的模型結(jié)果差別較大，根本原因如前所述，數(shù)值模型有一些差別。

圖9 受拉方板裂紋擴展路徑(裂紋應(yīng)變-相場建模)

圖10 受拉方板裂紋擴展路徑(改進的裂紋應(yīng)變-相場建模)

圖11 受拉方板的力-位移曲線

4.2 含單邊裂紋的受剪切方板

采用與算例4.1同樣的材料參數(shù)研究方板受剪切作用下的裂紋擴展過程，其計算模型如圖12所示，所有的邊界上豎向位移均約束為零。采用位移加載法，每一加載步位移增量為Δu=1×10-4mm， 采用非均勻網(wǎng)格，在預(yù)計裂紋擴展的位置加密網(wǎng)格，最小網(wǎng)格大小在0.001 mm左右，網(wǎng)格模型如圖13所示。

圖12 受剪切方板/mm

圖13 受剪切方板不均勻網(wǎng)格

首先采用幾何不連續(xù)法來預(yù)設(shè)裂紋，其計算得到的裂紋擴展路徑如圖14所示，若采用應(yīng)變-相場法來預(yù)設(shè)裂紋，計算得到的裂紋擴展路徑如圖15所示，可以發(fā)現(xiàn)與圖14差別較大，材料很久才開裂，且開裂路徑也不一致，主要原因是二者計算模型在該算例下差別較大，采用應(yīng)變 - 相場法預(yù)設(shè)裂紋時，初始裂紋附近并未完全弱化，而采用幾何不連續(xù)法預(yù)設(shè)裂紋時，裂紋面上下節(jié)點是分離的，可以有錯位。

圖14 受剪切方板裂紋擴展路徑(裂紋幾何建模,非均勻網(wǎng)格)

如果采用本文改進的應(yīng)變 - 相場法來預(yù)設(shè)裂紋，其得到的裂紋擴展路徑如圖16所示，可以發(fā)現(xiàn)，該計算結(jié)果與圖14的結(jié)果非常接近，證明了本文方法的有效性。

圖15 受剪切方板裂紋擴展路徑(裂紋應(yīng)變 - 相場建模，非均勻網(wǎng)格)

圖16 受剪切方板裂紋擴展路徑(裂紋改進的應(yīng)變 - 相場建模，非均勻網(wǎng)格)

為了進一步驗證本文方法的有效性，采用158404個均勻四邊形網(wǎng)格進行再次計算，如圖17所示，消除由于網(wǎng)格不均勻帶來的影響，取l0=0.005 mm，其計算得到的裂紋擴展路徑分別如圖18～20所示，可以發(fā)現(xiàn)，本文提出方法的結(jié)果與采用幾何不連續(xù)建模得到計算結(jié)果接近。

圖17 受剪切方板均勻網(wǎng)格

圖18 受剪切方板裂紋擴展路徑(裂紋幾何建模，均勻網(wǎng)格)

圖19 受剪切方板裂紋擴展路徑(裂紋應(yīng)變-相場建模，均勻網(wǎng)格)

圖20 受剪切方板裂紋擴展路徑(裂紋改進的應(yīng)變-相場建模，均勻網(wǎng)格)

計算得到的力 - 位移曲線如圖21示，可以發(fā)現(xiàn)，無論是均勻網(wǎng)格還是非均勻網(wǎng)格，本文提出的改進方法得到的荷載位移曲線與采用幾何不連續(xù)建模得到的計算結(jié)果吻合較好，而采用初始的應(yīng)變 - 相場建模得到的荷載要大很多。

圖21 受剪切方板荷載 - 位移曲線

5 結(jié) 論

本文通過改進的應(yīng)變 - 相場模型，對四種情況的準(zhǔn)靜態(tài)裂紋擴展問題進行了數(shù)值模擬，得出的主要結(jié)論如下：

(1)采用改進的應(yīng)變 - 相場法計算單邊裂紋受拉方板的開裂時間較幾何不連續(xù)法計算結(jié)果提前，這是由于改進的應(yīng)變 - 相場法弱化了更多單元造成的，但其相較于傳統(tǒng)的應(yīng)變-相場法模擬精度有大幅度提高；

(2)無論是均勻網(wǎng)格還是非均勻網(wǎng)格，三種模型計算單邊裂紋受剪切方板的結(jié)果表明改進的應(yīng)變-相場法計算出的裂紋擴展路徑與幾何不連續(xù)法的計算結(jié)果吻合較好，驗證了本文提出的改進的應(yīng)變-相場模型的合理性、可靠性，為模擬巖石或混凝土工程中的裂縫復(fù)雜擴展模式提供了一種新的思路。