高 宇， 張晨露， 岳 強(qiáng)， 王 橋， 周 偉， 程勇剛， 李通盛

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

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

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

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

1 準(zhǔn)靜態(tài)問(wèn)題的相場(chǎng)模型

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

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

(1)

(2)

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

(3)

式中：?為梯度符號(hào)。

圖1 彈性體中的裂紋面

1.2 正則化方程

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

(4)

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

圖2 s表示的裂紋用相場(chǎng)

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

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

(5)

1.3 控制方程

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

(6)

(7)

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

(8)

(9)

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

divσ=0 onΩ

(10)

(11)

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

(12)

(13)

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

圖3 彈性部分的邊界條件

圖4 相場(chǎng)部分的邊界條件

1.4 歷史場(chǎng)變量

(14)

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

(15)

其中尖括號(hào)的含義為：

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

(16)

式(5)便可以寫為：

(17)

相場(chǎng)方程可改寫為：

(18)

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

(19)

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

(20)

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

2 分步解耦算法

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

(21)

(22)

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

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

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

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

(23)

以二維情況為例，有

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

(24)

(25)

(26)

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

(27)

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

(28)

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

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

其中，對(duì)于平面應(yīng)變問(wèn)題有

(29)

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

(30)

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

(31)

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

(32)

(33)

(34)

4 數(shù)值算例

4.1 含單邊裂紋的受拉方板

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

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

圖6 受拉方板/mm

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

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

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

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

圖10 受拉方板裂紋擴(kuò)展路徑(改進(jìn)的裂紋應(yīng)變-相場(chǎng)建模)

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

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

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

圖12 受剪切方板/mm

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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