潘軍，王鴻，趙冰，徐亞星，李濤峰，龍承運(yùn)

（長(zhǎng)沙理工大學(xué) 土木工程學(xué)院，湖南 長(zhǎng)沙410114）

混凝土為重要的建筑材料，其斷裂特性對(duì)工程安全起著關(guān)鍵作用[1]?？煽?、高效的混凝土斷裂模型在橋梁、隧道、大壩等土木工程結(jié)構(gòu)的安全評(píng)估中發(fā)揮著重要作用[2]。對(duì)于混凝土斷裂的研究，尤其是對(duì)其裂紋萌生和擴(kuò)展的研究，引起了國(guó)內(nèi)外學(xué)者越來(lái)越多的關(guān)注。近年發(fā)展起來(lái)的相場(chǎng)斷裂法，通過(guò)跟蹤階參量自動(dòng)演化獲取裂紋路徑，可方便地模擬出裂紋的動(dòng)態(tài)擴(kuò)展過(guò)程[3]，為混凝土斷裂數(shù)值計(jì)算方法提供了新的方法。

傳統(tǒng)的混凝土斷裂模擬方法有單元?jiǎng)h除法、界面單元法和擴(kuò)展有限元法（XFEM）等，是以不連續(xù)位移場(chǎng)來(lái)描述裂紋。其中，單元?jiǎng)h除法和界面單元法只允許裂紋在網(wǎng)格邊界擴(kuò)展，存在明顯的網(wǎng)格依賴性[4]。擴(kuò)展有限元法（XFEM）[5]允許裂紋在網(wǎng)格內(nèi)擴(kuò)展，可以模擬出裂紋擴(kuò)展的任意路徑，但它在處理多裂紋問(wèn)題時(shí)較為繁瑣，大多用來(lái)處理單裂紋問(wèn)題和簡(jiǎn)單的多裂紋問(wèn)題。

相場(chǎng)斷裂法用連續(xù)函數(shù)描述裂紋，其核心是引入階參量描述材料的無(wú)損狀態(tài)和完全破壞狀態(tài)之間的各種狀態(tài)。相場(chǎng)斷裂法與其他方法相比，無(wú)需追蹤裂紋的幾何形狀，在計(jì)算復(fù)雜裂紋擴(kuò)展問(wèn)題的擴(kuò)展路徑時(shí)具有獨(dú)特優(yōu)勢(shì)。相場(chǎng)斷裂方法可分為：①以物理學(xué)Ginzburg-Laudau理論為基礎(chǔ)，可顯示多種相場(chǎng)模型[6]的研究結(jié)果，該類模型能較好地模擬混合型斷裂模型[7]，但Griffith臨界斷裂能本構(gòu)閾值Gc從未進(jìn)入相應(yīng)公式，還需適當(dāng)選取退化（耦合）函數(shù)，其物理意義并不明確。②以脆性斷裂變分理論為基礎(chǔ)，Bourdin 等人[8]結(jié)合相場(chǎng)理論，給出了含階參量斷裂面的彌散表達(dá)式。由材料類型[9]、材料特性[10]、加載方式[11]等因素作用的多種靜態(tài)相場(chǎng)斷裂模型及推廣至動(dòng)態(tài)斷裂[12-13]相場(chǎng)模型的結(jié)果表明：該類相場(chǎng)斷裂模型與傳統(tǒng)斷裂理論的緊密聯(lián)系，力學(xué)概念明確。

目前，相場(chǎng)斷裂法應(yīng)用于混凝土斷裂問(wèn)題的研究并不多見(jiàn)。Schlüter 等人[13]提出了預(yù)測(cè)準(zhǔn)脆性材料斷裂的相場(chǎng)模型，需預(yù)置初始裂紋。Mikeli?等人[14]提出了可解決加壓裂紋的相場(chǎng)模型，但裂紋局部網(wǎng)格較密，沒(méi)有描述網(wǎng)格密度對(duì)裂紋的影響。Hirshikesh 等人[15]提出了準(zhǔn)脆性材料相場(chǎng)模型，未能清楚描述長(zhǎng)度尺度參數(shù)l0與裂紋的關(guān)系。因此，作者擬進(jìn)行混凝土相場(chǎng)斷裂的Abaqus 二次開(kāi)發(fā)，研究單元尺寸和長(zhǎng)度尺度等參數(shù)對(duì)裂紋寬度的影響，分析混凝土破壞機(jī)理，確定合理混凝土斷裂參數(shù)，為混凝土結(jié)構(gòu)的穩(wěn)定性和安全性評(píng)估提供新的理論和計(jì)算方法。

1 混凝土的相場(chǎng)斷裂模型

1.1 斷裂總勢(shì)能

混凝土的斷裂總勢(shì)能由體積儲(chǔ)存能、裂紋表面能和外力做功三部分組成[8]：

式中：П為斷裂總勢(shì)能；Y為裂紋表面能；U為體積儲(chǔ)存能；W為外力功；Gc為Griffith 臨界斷裂能本構(gòu)閾值；ψ為體積儲(chǔ)存能密度；V為被積體積；Γ為完全斷裂的裂紋面；Sσ為應(yīng)力邊界；ui為位移；Ti為面力；fi為體力。

Bourdin 等人[8]在斷裂變分理論的基礎(chǔ)上，定義一個(gè)在[0,1]區(qū)間變化的標(biāo)量φ為裂紋相場(chǎng)（階參量）。當(dāng)φ=1 時(shí)，表示有裂紋；當(dāng)φ=0 時(shí)，表示材料完好。相場(chǎng)形式的裂紋表面能為：

對(duì)于二維問(wèn)題，單位體積的裂紋表面密度函數(shù)γ(φ)的具體形式為：

式中：l0為控制裂紋“擴(kuò)散”量的長(zhǎng)度尺度參數(shù)。

體積儲(chǔ)存能密度ψ由應(yīng)變?chǔ)舏j和階參量φ共同確定：

式中：H(εij)為彈性應(yīng)變能密度；Dijkl為彈性張量；εij、εkl為應(yīng)變張量；g(φ)為應(yīng)力退化函數(shù)。

g(φ)的表達(dá)式為：

式中：k≈10-6，使剛度矩陣在斷裂模擬時(shí)保持良好的非奇異性。

由式（4）～（5），體積儲(chǔ)存能U的表達(dá)式為：

1.2 相場(chǎng)斷裂變分理論

不考慮面力的情況下，將式（2）、（3）、（6）代入式（1），斷裂總勢(shì)能的表達(dá)式為：

斷裂總勢(shì)能的一階變分為：

式中：nj為邊界?V外法線的方向余弦；為彈性應(yīng)變能密度。

斷裂變分理論認(rèn)為，任意位置裂紋的任意動(dòng)態(tài)擴(kuò)展?fàn)顟B(tài)都會(huì)使總勢(shì)能最小。即在平衡狀態(tài)時(shí)，Dijklεkl為有效應(yīng)力；δП=0對(duì)任意的δui和δφ都成立，則由式(8)得：

斷裂問(wèn)題的控制方程為：

σij，j+fi= 0。 （9）

斷裂問(wèn)題的相場(chǎng)演化方程為：

式中：?為梯度。

力的邊界條件為：

σijnj= 0。 （11）

相場(chǎng)的自然邊界條件為：

相場(chǎng)斷裂變分理論繼承和發(fā)展了傳統(tǒng)的Griffith 理論，并解決了Griffith 理論無(wú)法解決的裂紋萌生、擴(kuò)展路徑和失穩(wěn)分叉等問(wèn)題，為研究斷裂問(wèn)題提供了新的視角。

1.3 相場(chǎng)斷裂模型的有限元形式

在平面應(yīng)力狀態(tài)下，單元內(nèi)部一點(diǎn)的位移u與階參量φ可以離散為：

u=Nuae，φ=Nφφe。 （13）

其中，ae={u1,v1,u2,v2,…,un,vn}T為單元節(jié)點(diǎn)位移向量；φe={φ1,φ2,…,φn}T為單元節(jié)點(diǎn)階參量向量；n為單元節(jié)點(diǎn)數(shù)。

Nu=，為位移的形函數(shù)矩陣；Nφ=[ ]N1N2…Nn為階參量的形函數(shù)矩陣。

單元內(nèi)任意一點(diǎn)的應(yīng)變?chǔ)排c階參量的梯度?φ可由ae和φe表達(dá)式為：

ε=Buae，?φ=Bφφe。 （14）

其中，

單元節(jié)點(diǎn)位移向量ae與單元節(jié)點(diǎn)階參量向量φe用結(jié)構(gòu)節(jié)點(diǎn)位移a與節(jié)點(diǎn)階參量φ表達(dá)式為：

ae=Ga，φe=Gφ。 （15）

式中：G為單元集成的組裝矩陣。

將式（13）～（15）代入式（7）中，結(jié)構(gòu)總勢(shì)能П的有限元形式為：

式中：ne為單元數(shù)。

由變分原理，總勢(shì)能П的一階變分為零，即：

δП=0對(duì)任意的δaT和δφT都成立，則：

平衡方程可改寫(xiě)為：

式中：Ku和Kφ分別為結(jié)構(gòu)整體的位移和相場(chǎng)剛度，Pu和Pφ分別為結(jié)構(gòu)位移和相場(chǎng)的整體荷載向量。

對(duì)式（18）～（19）做微分運(yùn)算，得：

增量平衡方程可寫(xiě)為：

其中，

由于相場(chǎng)斷裂有限元模型的方程組是非線性，因此，必須采用增量迭代來(lái)計(jì)算解。有限元軟件Abquas 能解決許多復(fù)雜的非線性問(wèn)題，被廣泛應(yīng)用工程計(jì)算中。選擇在Abaqus 軟件中建立模型，以便利用其內(nèi)置的非線性求解器，采用牛頓-拉夫森算法和自動(dòng)時(shí)間步進(jìn)方案。先使用Abaqus 的用戶自定義單元子程序UEL，定義新的二維四節(jié)點(diǎn)等參單元，該單元每個(gè)節(jié)點(diǎn)具有1個(gè)位移自由度和1 個(gè)相場(chǎng)自由度，共3 個(gè)自由度，UEL 還會(huì)計(jì)算單元的切線剛度矩陣和節(jié)點(diǎn)力向量。再編寫(xiě)1個(gè)用戶自定義材料子程序（UMAT），該子程序用于公共塊傳遞信息，并在積分點(diǎn)之間進(jìn)行插值。本研究基于Abquas 的UEL 和UMAT 實(shí)現(xiàn)混凝土相場(chǎng)斷裂的Abquas二次開(kāi)發(fā)。

2混凝土相場(chǎng)斷裂的模擬結(jié)果

2.1相場(chǎng)斷裂模型的基本驗(yàn)證

通過(guò)單向拉伸試驗(yàn)，驗(yàn)證混凝土相場(chǎng)斷裂模型的可靠性，得到應(yīng)力應(yīng)變曲線為[11]：

式中：B、C的取值分別為0.499 8、1.001 0。

試驗(yàn)材料的彈性模量E=25.5×103MPa，抗拉強(qiáng)度f(wàn)t=1.5 MPa，泊松比μ=0.2，受拉峰值應(yīng)變?chǔ)舉=1×10-4，本構(gòu)閾值為Gc=15.500 50 N/m，長(zhǎng)度尺度參數(shù)l0=0.018 48 m?；炷羻蜗蚶煸囼?yàn)與相場(chǎng)斷裂模型的應(yīng)力應(yīng)變曲線如圖1 所示。從圖1 中可以看出，點(diǎn)B 為峰值應(yīng)變?chǔ)舉對(duì)應(yīng)的峰值應(yīng)力點(diǎn)。當(dāng)曲線處于AB 段時(shí)，ε<εe，單向拉伸試驗(yàn)與相場(chǎng)斷裂模型的應(yīng)力應(yīng)變曲線吻合程度較好。當(dāng)曲線處于BC 段時(shí)，曲線呈下凹趨勢(shì)，此時(shí)材料處于軟化階段，裂紋開(kāi)始萌生擴(kuò)展，裂紋寬度隨著應(yīng)力的減小而不斷增加，單向拉伸試驗(yàn)與相場(chǎng)斷裂模型的應(yīng)力應(yīng)變曲線略有差別。當(dāng)曲線處于CD 段時(shí)，材料處于裂紋貫通階段，隨著應(yīng)變的增加，應(yīng)力緩慢減小，單向拉伸試驗(yàn)與相場(chǎng)斷裂模型的應(yīng)力應(yīng)變曲線吻合。表明：相場(chǎng)斷裂模型能較好地模擬出混凝土的拉伸變形過(guò)程。

圖1混凝土應(yīng)力應(yīng)變曲線Fig.1Stress strain curves of concrete

2.2 單元尺寸對(duì)裂紋寬度的影響

為研究單元尺寸對(duì)裂紋寬度模擬結(jié)果的影響，采用Abaqus的二維相場(chǎng)斷裂模型，利用5組不同單元尺寸的網(wǎng)格對(duì)混凝土梁三點(diǎn)彎曲試驗(yàn)進(jìn)行模擬。三點(diǎn)彎曲試樣幾何模型如圖2所示，材料參數(shù)為：彈性模量E=25.5×103MPa，抗拉強(qiáng)度f(wàn)t=1.5 MPa，泊松比μ=0.2，受 拉 峰 值應(yīng) 變?chǔ)舉=1×10-4，本 構(gòu) 閾 值Gc=15.500 5 N/m，長(zhǎng)度尺度參數(shù)l0=0.018 48 m。5組不同單元尺寸網(wǎng)格單元的最小長(zhǎng)度、寬度見(jiàn)表1。

圖2 三點(diǎn)彎曲試驗(yàn)試件幾何參數(shù)（單位：m）Fig.2 Geometry parameters of specimen in three point bending test(unit:m)

表1 網(wǎng)格單元尺寸Table 1 Grid element sizes m

在試件跨中頂部采用增量位移加載，初始位移為u0=10-6m，加載至u=7.70×10-5m，加載步Δu=10-6m，模擬結(jié)果如圖3 所示?；炷亮褐行詫拥碾A參數(shù)（裂紋）分布及不同的單元長(zhǎng)度下的裂紋寬度如圖4所示。

從圖3 和圖4(a)中可以看出，裂紋寬度隨著單元尺寸的減小而變窄。當(dāng)單元尺寸≤Mesh3 時(shí)，裂紋寬度趨于穩(wěn)定。從圖4(b)中可以看出，當(dāng)裂紋完全貫穿部分，寬度的模擬結(jié)果隨單元長(zhǎng)度的減小而線性減小。當(dāng)單元的長(zhǎng)度足夠?。ㄚ呌诹悖r(shí)，裂紋寬度收斂至0.003 92 m?？梢哉J(rèn)為收斂值是由理論模型決定的，與單元尺寸無(wú)關(guān)的物理量。

圖3 不同單元尺寸下相場(chǎng)斷裂模型三點(diǎn)彎曲試樣的裂紋寬度Fig.3 Crack widths of three point bending specimens in phase field fracture model under different element sizes

圖4 相場(chǎng)斷裂模型的三點(diǎn)彎曲試驗(yàn)值Fig.4 Values of three point bending tests obtained by phase field fracture model

2.3 長(zhǎng)度參數(shù)對(duì)裂紋寬度的影響

為了研究長(zhǎng)度尺度參數(shù)l0對(duì)裂紋寬度的影響，基于Abaqus 的二維相場(chǎng)斷裂模型，利用Mesh4 網(wǎng)格，研究不同l0時(shí)，混凝土梁三點(diǎn)彎曲試驗(yàn)的裂紋擴(kuò)展，其余的材料參數(shù)不變，模擬結(jié)果如圖5 所示?；炷亮褐行詫拥碾A參數(shù)（裂紋）分布及不同l0下的裂紋寬度如圖6所示。

從圖5 和圖6(a)中可以看出，裂紋寬度隨l0減小而變小，裂紋的彌散程度隨l0減小而顯著變小。從圖6(b)還可以看出，裂紋寬度及彌散程度的模擬結(jié)果隨l0減小而線性減小。當(dāng)l0足夠?。ㄚ呌诹悖r(shí)，裂紋的寬度收斂至零，表明：裂紋的寬度與l0呈正比例關(guān)系，可以認(rèn)為裂紋寬度及裂紋彌散程度受l0控制。因此，在試驗(yàn)測(cè)得裂紋寬度的前提下，可以利用本數(shù)值模型，采用數(shù)值打靶法確定材料l0取值。

圖5 不同長(zhǎng)度尺度的相場(chǎng)斷裂模型三點(diǎn)彎曲試驗(yàn)的裂紋寬度Fig.5 Crack widths of three point bending specimens in phase field fracture model under different length scale parameters

圖6 不同長(zhǎng)度的相場(chǎng)斷裂模型三點(diǎn)彎曲試驗(yàn)值Fig.6 Values of three point bending tests obtained by phase field fracture model under different length scale parameters

3 結(jié)論

使用UEL 子程序?qū)崿F(xiàn)相場(chǎng)斷裂模型的Abaqus二次開(kāi)發(fā)，編寫(xiě)了附加的UMAT 子程序，實(shí)現(xiàn)了相場(chǎng)斷裂模擬結(jié)果的可視化。利用數(shù)值算例分析了單元尺寸和長(zhǎng)度尺度參數(shù)對(duì)三點(diǎn)彎曲混凝土梁裂紋擴(kuò)展寬度的影響。得到結(jié)論為：

1）混凝土相場(chǎng)斷裂法可以較好地模擬出單向拉伸試驗(yàn)的全過(guò)程。

2）當(dāng)單元尺寸足夠小時(shí)，裂紋寬度是一個(gè)由理論模型決定、與單元尺寸無(wú)關(guān)的物理量。

3）裂紋寬度和裂紋彌散程度受長(zhǎng)度尺度參數(shù)控制。