葉少鋒，周 偉，常曉林，馬 剛

(1. 武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，武漢 430072；2. 武漢大學(xué) 水工巖石力學(xué)教育部重點實驗室，武漢 430072)

0 引 言

在眾多壩型，土石壩因為其具有充分利用當?shù)夭牧?、對壩址地形地質(zhì)條件要求較低、抗震性好等優(yōu)勢而在壩工設(shè)計中得到十分廣泛的運用。其中心墻壩是土石壩常用的壩型之一。眾所周知，心墻水力劈裂是大壩安全的一個重要隱患。在心墻壩的設(shè)計過程中，水力劈裂是一個必須重點研究和防范的關(guān)鍵問題。

許多學(xué)者對心墻水力劈裂問題進行了許多系統(tǒng)深入并且卓有成效的研究。黃文熙[1]在1982年指出水力劈裂是土體中因水壓力抬高導(dǎo)致裂縫產(chǎn)生與擴展的一種物理現(xiàn)象。這一定義很好地反映了水力劈裂現(xiàn)象的本質(zhì)。張丙印等[2]提出一種水壓楔劈模型來解釋水力劈裂裂縫的形成機理。朱俊高等[3]認為高水力梯度是形成水力劈裂的最根本因素。以上三種對水力劈裂形成機理的觀點具有一定的相似之處。在心墻水力劈裂判定方面，殷宗澤等[4]提出將心墻內(nèi)的總應(yīng)力與心墻表面受到的水壓力進行對比，以此分析是否會發(fā)生水力劈裂，其中總應(yīng)力是由土體中有效應(yīng)力與孔隙水壓力疊加得到的。在總結(jié)并比較總應(yīng)力法、綜合法和有效應(yīng)力法三種水力劈裂判別方法的基礎(chǔ)上，陳五一等[5]提出一種新的安全系數(shù)來判定水力劈裂。在數(shù)值模擬方面，李全明等[6]通過綜合比奧固結(jié)理論與彌散裂縫理論，提出了一種有限元計算方法來數(shù)值模擬水力劈裂現(xiàn)象。不過該方法需要預(yù)設(shè)裂縫位置，而且無法模擬裂縫貫穿整個心墻的過程。楊艷等[7]從細觀角度出發(fā)，應(yīng)用離散元方法研究心墻水力劈裂的過程，該研究表明裂縫形成的原因是土體受到的張拉應(yīng)力大于抗拉強度。在水力劈裂室內(nèi)物理試驗方面，袁俊平等[8]自制了一套水力劈裂試驗設(shè)備，在心墻試樣中設(shè)置不同的初始裂縫長度，分析試驗結(jié)果得出：試樣的初始裂縫長度越大，則形成水力劈裂所需的水壓力越小。

綜合以上學(xué)者的研究可以看出，心墻水力劈裂是高水力梯度作用下導(dǎo)致裂縫發(fā)生與擴展的一種流固耦合的物理現(xiàn)象?，F(xiàn)有的研究方法或理論主要關(guān)注于土料應(yīng)力、孔隙水壓力和裂縫等方面的其中一兩點，不夠全面。作為一種能夠有效模擬裂縫萌生與擴展的數(shù)值模型，內(nèi)聚力模型已經(jīng)廣泛運用于巖石等多種材料的開裂行為研究[9-11]，也為分析水力劈裂過程中裂縫的擴展行為提供了新手段。本文采用考慮孔隙水壓力的內(nèi)聚力模型對心墻水力劈裂現(xiàn)象進行了數(shù)值模擬，并且分析了心墻水力劈裂的發(fā)生機理。該方法不僅模擬了流固耦合作用，而且能夠直觀地跟蹤裂縫萌生、擴展直至最終形成貫穿性裂縫的整個過程中孔隙水壓力、應(yīng)力變形等各個方面的狀態(tài)。通過具體分析，該數(shù)值模擬所獲得的成果符合現(xiàn)有的室內(nèi)物理試驗結(jié)果與理論，為今后心墻壩設(shè)計人員分析和解決水力劈裂問題提供了新手段。

1 考慮孔隙水壓力的內(nèi)聚力模型的基本原理

1.1 基本思路

本文數(shù)值模擬是在連續(xù)—離散耦合分析方法(FDEM)的框架下實現(xiàn)的。FDEM方法最早是由Munjiza等提出，其基本思路是在心墻結(jié)構(gòu)內(nèi)劃分有限元網(wǎng)格(本文采用三角形單元)，然后在所有相鄰單元的公共邊上插入無厚度的界面單元，通過界面單元的斷裂來模擬裂縫的萌生與擴展，具體內(nèi)容可以參考文獻[12-17]。其中界面單元是該方法的關(guān)鍵所在，本文采用內(nèi)聚力模型來定義界面單元的力學(xué)行為。由于內(nèi)聚力模型只能模擬裂縫的形成，為了實現(xiàn)對水力劈裂流固耦合作用的模擬，必須在內(nèi)聚力模型中考慮滲流作用，建立考慮孔隙水壓力的內(nèi)聚力模型，簡稱孔壓內(nèi)聚力模型。其模擬滲流的思路是，將流體看作是不可壓縮的，并且流體前緣與裂縫尖端是重合的，計算時不考慮流體的滯后效應(yīng)；有限元網(wǎng)格內(nèi)部孔隙流體流動滿足Darcy定律，界面單元裂縫內(nèi)流體流動包括沿裂縫面的切向流動和垂直裂縫面的法向流動。在水壓力的驅(qū)動下，流體的切向流動驅(qū)使裂縫向前擴展，法向流動使裂縫內(nèi)流體滲透到周邊介質(zhì)中。在計算過程中，水壓力隨著滲流作用進行傳遞，改變模型應(yīng)力狀態(tài)，影響界面單元地張開與閉合狀態(tài)；反之，模型的應(yīng)力狀態(tài)和界面單元地張開與閉合狀態(tài)的變化必然引起孔隙與裂縫內(nèi)流體滲流的變化。兩者形成相互制約和相互作用的關(guān)系。這種關(guān)系就是流固耦合。

在FDEM框架下運用孔隙水壓力內(nèi)聚力模型模擬水力劈裂有兩個優(yōu)勢。一是能夠克服連續(xù)介質(zhì)力學(xué)模型模擬裂縫的困難，有效模擬裂縫的萌生與擴展；二是既能考慮水力劈裂過程中孔隙流體和裂縫流體滲流作用的差異，又能反映裂縫流體與孔隙流體的流量交換。

1.2 孔隙流體滲流方程

假定心墻顆粒骨架和流體均具有不可壓縮性，流體的連續(xù)性方程為：

(1)

式中：εv為體積應(yīng)變；q為滲透流速矢量。

流體滲流的 Darcy 定律為：

q=k▽h

(2)

式中：k為滲透系數(shù)矢量；h為水頭。

1.3 單元應(yīng)力平衡方程

單元的應(yīng)力平衡方程為：

σij,j+fi=0

(3)

式中：σij表示總應(yīng)力；fi表示體力。

本文水力劈裂分析時均采用有效應(yīng)力進行計算，其中應(yīng)力狀態(tài)以拉應(yīng)力為正，孔隙水壓力pw以壓應(yīng)力為正，有效應(yīng)力σ′ij的表達式如下：

σ′ij=σij+δijpw

(4)

由式(3)和式(4)可得出用有效應(yīng)力表示的平衡方程：

σ′ij-pw,j+fi=0

(5)

1.4 內(nèi)聚力模型

內(nèi)聚力模型理論起源于斷裂力學(xué)。在20世紀60年代，Dugdale[18]和Barrenblatt[19]為了解決線彈性斷裂力學(xué)中裂縫尖端處應(yīng)力奇異性問題先后提出了內(nèi)聚力模型(CZM)理論。該理論認為裂縫尖端處應(yīng)力為裂縫開度的函數(shù)。本文通過采用內(nèi)聚力模型定義界面單元受到荷載作用時應(yīng)力與位移的響應(yīng)關(guān)系來描述界面單元的力學(xué)行為。界面單元變形分為彈性變形、起裂、損傷演化和完全失效4個階段；單元起裂之前，發(fā)生彈性變形；起裂后，進入損傷階段，剛度逐漸降低；當剛度降低至0時，單元完全失去承載能力(如圖1所示)。

圖1 內(nèi)聚力模型與界面單元Fig.1 Schematic representation of the cohesive zone model and cohesive interface element

在界面單元出現(xiàn)損傷之前，其線彈性的本構(gòu)關(guān)系如下：

(6)

式中：tn、ts分別表示法向應(yīng)力與切向應(yīng)力；Kij為剛度矩陣；εn、εs分別表示法向應(yīng)變與切向應(yīng)變，本文計算時未考慮應(yīng)力、應(yīng)變的法向分量與切向分量之間的耦合關(guān)系。

法向應(yīng)變、切向應(yīng)變定義為：

(7)

式中：δn、δt分別表示法向位移和切向位移；T0為界面單元的本構(gòu)厚度。因為本文使用無厚度的界面單元，所以取T0=1，界面單元的應(yīng)變與相應(yīng)方向上的位移大小相等。界面單元的線彈性本構(gòu)方程可寫成：

(8)

式中：kn、ks分別表示法向剛度和切向剛度。

為了模擬裂縫萌生過程，需要定義有效的起裂準則。本文假定材料破壞是由拉應(yīng)力和剪應(yīng)力共同作用下導(dǎo)致的，所以模擬時采用二次應(yīng)力準則判斷界面單元是否開裂。當各個方向的應(yīng)力與其相應(yīng)的臨界應(yīng)力的比值的平方和達到1時，開始形成裂縫，公式如下：

(9)

式中：t0n為界面單元法向應(yīng)力的臨界值，即為心墻土料的抗拉強度；t0s為界面單元切向應(yīng)力的臨界值，即為心墻土料的抗剪強度；〈〉為Macaulay括號：

(10)

界面單元起裂后，即進入裂縫擴展階段，此時采用標量損傷因子來表示界面單元剛度的退化的程度。單元應(yīng)力與位移的表達式為：

(11)

ts=(1-D)ksδs

式中：D表示標量損傷因子，當D=0時，表明界面單元未出現(xiàn)損傷；當D=1時，表明單元完全裂開，失去承載能力。

定義界面單元的等效應(yīng)力和等效位移如下：

(12)

式中：teff為等效應(yīng)力；δm為等效位移。

通過定義界面單元在損傷過程中耗散的能量，即可確定單元的損傷演化過程。耗散能量是指通常所說的斷裂能，它等于應(yīng)力-分離曲線下面所包含面積的大小。本文損傷演化準則采用基于斷裂能的Benzeggagh-Kenane 準則[20]，其表達式如下：

GS=Gs,GT=GS+Gn

(13)

式中：GC為復(fù)合斷裂能；GCn為Ⅰ型斷裂能；GCs為Ⅱ型斷裂能；本文中，η取2。

基于斷裂能的損傷演化準則中界面單元的軟化準則分為線性軟化和指數(shù)軟化。由于指數(shù)軟化曲線更接近心墻料拉伸斷裂的應(yīng)力-應(yīng)變曲線，本文采用指數(shù)軟化準則。指數(shù)軟化準則中，標量損傷因子的表達式如下：

(14)

式中：G0為界面單元開始出現(xiàn)損傷時對應(yīng)的彈性能；δ0m、δfm分別指界面單元開始出現(xiàn)損傷時與完全失效時對應(yīng)的等效位移。

1.5 水力裂縫模型

在考慮孔隙水壓力的內(nèi)聚力模型中，裂縫中流體的流動現(xiàn)象如圖2所示。流體在裂縫中的流動可分為法向流動和切向流動，切向流動驅(qū)使裂縫向前擴展，法向流動反映了流體的濾失效應(yīng)，即部分流體滲透到周邊介質(zhì)中。

圖2 界面單元流體流動Fig.2 Flow within cohesive interface element

流體的切向流動遵守潤滑方程：

qd=-kt▽p

(15)

式中：q為界面單元的體積流率密度矢量；d為界面單元地張開寬度；▽p為沿著界面單元的水壓力梯度；kt為切向滲透系數(shù)，由雷諾方程求得：

(16)

式中：μ為流體的黏滯系數(shù)。

假設(shè)裂縫的上、下表面為可透水層，其滲透系數(shù)即為界面單元的濾失系數(shù)，在水壓力的驅(qū)動下流體產(chǎn)生法向流動，法向流速的計算公式如下：

qt=ct(pi-pt)

qb=cb(pi-pb)

(17)

式中：qt和qb分別為界面單元上、下表面的流速；ct和cb分別為界面單元上、下表面的濾失系數(shù)；pt和pb分別為與界面單元上、下表面相鄰的多孔介質(zhì)的孔隙水壓力；pi為界面單元內(nèi)部的孔隙水壓力。

2 基于內(nèi)聚力模型的心墻水力劈裂數(shù)值模擬

2.1 內(nèi)聚力模型力學(xué)參數(shù)選取

模型參數(shù)的選取是數(shù)值模擬的一個關(guān)鍵因素。筆者查閱了大量文獻，尚未找到將內(nèi)聚力模型運用于土石壩心墻材料研究的文獻資料，因而缺乏心墻材料的內(nèi)聚力模型力學(xué)參數(shù)的參考數(shù)據(jù)。為此，筆者在進行心墻水力劈裂的數(shù)值模擬前，先采用內(nèi)聚力模型對文獻[21]中的心墻土料試樣的巴西圓盤劈裂試驗進行數(shù)值模擬，初步確定心墻材料內(nèi)聚力模型參數(shù)。由于堆石壩橫河向的寬度一般遠大于順河向的寬度，結(jié)構(gòu)變形特點更接近于平面應(yīng)變，所以本文以平面應(yīng)變的角度來數(shù)值模擬巴西圓盤劈裂和水力劈裂試驗。

巴西圓盤劈裂物理試驗采用直徑為39.1 mm，高度為80 mm的圓盤試樣。試驗分別選取5組不同含水率的試樣和5組不同密度的試樣。本文對其中密度最大且比較接近工程實際情況的一組試樣進行模擬，其密度為1.65 g/cm3，含水率為16%。數(shù)值模型為直徑39.1 mm的圓形平面模型，劃分為4 464 個界面單元和3 024 個三角形實體單元，具體模型如圖3所示。數(shù)值模擬時通過上下兩個剛性板以一定的速率向圓盤中心移動來進行加載，本文剛性板的加載速率為9×10-3m/s。隨著剛性板的移動，圓盤逐漸形成上下貫通的裂縫，如圖4所示。計算結(jié)束后，通過提取剛性板受到的反力與圓盤裂縫方向的徑向位移，繪制出拉應(yīng)力與徑向位移的關(guān)系曲線并與文獻[21]的物理試驗結(jié)果進行比較。其中由反力計算拉應(yīng)力的公式如下：

(18)

式中：σt為拉應(yīng)力；F為剛性板受到的反力；w和l分別表示試樣直徑和厚度。其中，對于二維數(shù)值模擬，試樣的厚度取為1。

圖3 巴西圓盤劈裂試樣網(wǎng)格劃分Fig.3 Mesh configuration of Beazilian samples

圖4 巴西圓盤劈裂數(shù)值試樣裂縫圖Fig.4 Cracks of numerical Beazilian samples

運用內(nèi)聚力模型所需的力學(xué)參數(shù)如表1所示，包括密度ρ、彈性模量E、泊松比υ、抗拉強度tn、抗剪強度ts、法向剛度kn、切向剛度ks、Ⅰ型斷裂能GCn、Ⅱ型斷裂能GCs等。本文數(shù)值模型中，沒有考慮到土料的含水量這一因素，但是含水量不同，試樣的力學(xué)性能也會不同，含水量對試樣的影響會通過力學(xué)參數(shù)體現(xiàn)出來。數(shù)值模擬時，不斷調(diào)整各個參數(shù)進行試算，使試樣最大拉應(yīng)力與其對應(yīng)的徑向位移與物理試樣結(jié)果相吻合。最終各個參數(shù)的取值如表1所示，數(shù)值試驗與物理試驗的拉應(yīng)力和徑向位移曲線如圖5所示。從圖5可以看出，數(shù)值試驗曲線在拉應(yīng)力達到最大值后，急劇下降至0，而物理試驗曲線則沒有這一過程。這主要是由于物理試驗最后壓裂試樣的過程很短暫，儀器來不及記錄數(shù)據(jù)，試樣就已經(jīng)完全失去承載能力，而數(shù)值試驗則沒有這個問題，因而數(shù)值試驗記錄了最后壓裂的過程曲線。除了上面這一差異，該數(shù)值試驗的應(yīng)力變形曲線基本與物理實驗結(jié)果相符合，因此初步確定的內(nèi)聚力模型參數(shù)能夠模擬心墻試樣的劈裂力學(xué)行為。

表1 內(nèi)聚力模型力學(xué)參數(shù)Tab.1 Mechanical parameters of cohesive zone model

圖5 物理試驗與數(shù)值試驗的拉伸應(yīng)力-位移曲線對比[21]Fig.5 Experimental tensile stress-displacement response comparised to numerical predictions[21]

2.2 心墻水力劈裂數(shù)值模擬

2.2.1 建立計算模型

本文選取文獻[2]中的室內(nèi)水力劈裂試驗數(shù)據(jù)資料進行水力劈裂二維數(shù)值模擬。該實驗物理模型如圖6所示，心墻試樣為邊長20 cm的立方體，試樣中部有一層長10 cm、厚0.5 cm的砂層，表示心墻的初始裂縫或者缺陷。進行試驗時，先在模型的左側(cè)和上側(cè)通過氣囊施加初始應(yīng)力σ3和σ1使試樣固結(jié)，然后在左側(cè)施加水壓力，通過觀察窗口和測量滲流量的變化判斷是否形成水力劈裂。

圖6 物理模型示意圖[2]Fig.6 Sketch of physical model[2]

數(shù)值模型取20 cm×20 cm的正方形，模型中部有一道長為10 cm的初始裂縫或缺陷。模型由3 628 個實體單元和5 362 個界面單元組成，如圖7所示。由于水力劈裂數(shù)值試驗是采用考慮孔隙水壓力作用的內(nèi)聚力模型，所以計算時所用參數(shù)除了表1中力學(xué)參數(shù)外，還有流體計算參數(shù)，包括滲透系數(shù)k、黏滯系數(shù)μ、重度γw、孔隙比e和濾失系數(shù)ct與cb等，參數(shù)的選取參照文獻[22]，具體數(shù)值見表2。

表2 流體計算參數(shù)Tab.2 Computational parameters for fluid

圖7 數(shù)值模型示意圖Fig.7 Sketch of numerical model

2.2.2 數(shù)值模擬結(jié)果分析

本文模擬心墻水力劈裂的過程如下：首先對模型上表面和上游表面分別施加σ1和σ3的面荷載進行固結(jié)；然后在模型上游施加水壓力；一開始取一較小的水壓力值，如果計算結(jié)果沒有形成裂縫，則設(shè)置比較大的水壓力值重新進行計算；如此重復(fù)上述步驟，直到計算結(jié)果出現(xiàn)貫穿上下游的裂縫為止(如圖8所示)，此時施加的水壓力就是劈裂水壓力值。

圖8 水力劈裂后的貫穿性裂縫Fig.8 Penetrating cracks after hydraulic fracturing

本文根據(jù)室內(nèi)物理試驗σ1、σ3的3組初始應(yīng)力組合方案(見表3)分別進行了水力劈裂數(shù)值模擬，數(shù)值模擬與物理試驗的劈裂水壓力值對比如表3與圖9所示。根據(jù)圖9所呈現(xiàn)的規(guī)律可以得出，隨著初始豎向應(yīng)力的增大，數(shù)值模擬與物理試驗的劈裂水壓力值均逐漸增大，劈裂水壓力值與豎向應(yīng)力大致呈線性相關(guān)。另外，數(shù)值模擬的劈裂水壓力值要大于室內(nèi)試驗的劈裂水壓力。筆者認為這可能源于3個原因。第一，室內(nèi)心墻材料的物理特性受含水率、密度等因素的影響很大，兩個試驗并非采用同一心墻土料進行實驗，水力劈裂試樣的力學(xué)參數(shù)與巴西圓盤劈裂試樣可能有所差異，本文直接采用巴西圓盤劈裂試樣的參數(shù)來模擬水力劈裂，可能導(dǎo)致結(jié)果有一定誤差。第二，室內(nèi)物理試驗在達到水力劈裂破壞前，經(jīng)過了多次水的滲透作用，心墻土在水的反復(fù)作用下發(fā)生軟化作用，強度有所下降，而數(shù)值模擬沒能考慮到這一因素的影響。第三，數(shù)值模擬在考慮實際情況的基礎(chǔ)上必須做一定合理假定和簡化，這也會使數(shù)值試驗與物理試驗結(jié)果之間產(chǎn)生一些無法避免的誤差。

表3 數(shù)值模擬結(jié)果[2]Tab.3 Numerical results[2]

圖9 數(shù)值模擬結(jié)果與室內(nèi)試驗結(jié)果對比[2]Fig.9 Numerical and laboratory results[2]

為了從水壓力角度探究心墻水力劈裂的發(fā)生機理，從第二組應(yīng)力組合的數(shù)值試驗分別提取水平裂縫內(nèi)和垂直裂縫路徑上的水壓力值。圖10為裂縫擴展至距上游面0.15 m處時刻的水平裂縫內(nèi)水壓力分布圖。從圖10 可以看出，水壓力分布曲線可分為4段。第1段曲線橫坐標區(qū)間為距離上游面0～0.10 m，此處為模型的初始裂縫或者缺陷，水壓力等于上游面的水壓力值0.33 MPa。第2段曲線橫坐標區(qū)間為距離上游面0.10～0.14 m，此處界面單元損傷因子已經(jīng)到達1，界面單元完全裂開形成張開裂縫，隨著水流進入裂縫，水壓力由初始裂縫向新生裂縫處傳播，由于水的黏滯作用，水壓力在向前傳播的過程中逐漸減弱。第3段曲線橫坐標區(qū)間為0.14～0.15 m，水壓力急劇下降，呈現(xiàn)很大的水力梯度。此處為裂縫的最前端，在高水力梯度的作用下，界面單元出現(xiàn)損傷，損傷因子介于0～1，裂縫由閉合狀態(tài)逐漸張開。第4段橫坐標區(qū)間為距離上游面0.15～0.20 m，水壓力值接近于0。此處界面單元完好未出現(xiàn)損傷，沒有裂縫出現(xiàn)，水壓力以Darcy定律在介質(zhì)中傳播，水壓力值接近于下游面的初始水壓力值。圖11為垂直穿過距離上游面0.12 m處水平裂縫的路徑上的水壓力值分布圖，提取水壓力時刻與圖10一樣，其中橫坐標等于0表示裂縫的位置，橫坐標大于0表示位于裂縫的下方，橫坐標小于0表示位于裂縫的上方。從圖11中可以看出，裂縫中的水壓力要遠大于裂縫上下方介質(zhì)中的孔隙水壓力值。這是由于裂縫形成后，初始裂縫內(nèi)的水壓力迅速向裂縫傳播，而裂縫上下方介質(zhì)的孔隙水壓力則是上游面初始水壓力滿足Darcy定律向下游傳播的結(jié)果；在心墻材料較低滲透性與裂縫內(nèi)水壓力未能及時向周邊介質(zhì)消散的條件下，裂縫內(nèi)水壓力擴散速率要大于孔隙水壓傳播速度，因而裂縫內(nèi)水壓力大于上下方介質(zhì)中的孔隙水壓力。在裂縫內(nèi)與上下方介質(zhì)中的水壓力差的作用下，上下方的介質(zhì)分別向上與向下移動，裂縫逐漸張開，直到裂縫內(nèi)水壓力向周邊介質(zhì)消散達到水壓平衡為止。綜上所述，水力劈裂過程中，高水力梯度和裂縫與周邊介質(zhì)的水壓差是劈裂裂縫前端的兩個重要特征。結(jié)合文獻[2]描述的“水壓楔劈效應(yīng)”可以看出，高水力梯度和裂縫與周邊介質(zhì)的水壓差是裂縫萌生與擴展的兩個重要原因。

圖10 裂縫內(nèi)的水壓力Fig.10 Pressure in the fractures

圖11 垂直水平裂縫路徑上的水壓力Fig.11 Pore pressure in a vertical path

為了初步探討水力劈裂過程中界面單元的應(yīng)力狀態(tài)，從第二組數(shù)值試驗選取開裂路徑上任意一個界面單元，提取其切向應(yīng)力、法向應(yīng)力以及損傷因子的歷時曲線，如圖12所示。圖中，法向應(yīng)力正值表示拉伸應(yīng)力，切向應(yīng)力正值表示逆時針方向，0～1 s的時間段為施加初始應(yīng)力的階段，1～3 s為施加水壓力階段。從圖中可以看出，法向應(yīng)力在施加初始應(yīng)力階段由0逐漸增加到0.3 MPa的壓應(yīng)力；在施加水壓力階段的前期，還未產(chǎn)生裂縫，孔隙水壓力抵消了部分初始應(yīng)力，壓應(yīng)力逐漸變小到0.2 MPa左右；在2.2 s時刻界面單元開始出現(xiàn)損傷，上游水壓力傳遞到裂縫中，法向應(yīng)力由壓應(yīng)力減小至0并轉(zhuǎn)變成拉應(yīng)力；當損傷因子達到1時，界面單元完成裂開失去承載能力，拉應(yīng)力全部釋放，法向應(yīng)力降低至0。界面單元在施加水壓力的前期受滲流作用存在部分切向應(yīng)力，在界面單元裂開過程中，切向應(yīng)力逐漸降低至0。從圖中的曲線可以發(fā)現(xiàn)，界面單元在開裂過程中，法向應(yīng)力的變化幅度以及起裂時的應(yīng)力值要遠高于切向應(yīng)力。從以上的分析中可以看出，在二維水力劈裂數(shù)值模擬中，水力劈裂的破壞機理是以拉伸破壞為主。值得注意的是，在三維情況下，心墻結(jié)構(gòu)的應(yīng)力狀態(tài)要復(fù)雜得多，其水力劈裂的破壞機理還有待進一步研究。

圖12 界面單元切向應(yīng)力、法向應(yīng)力和損傷因子歷時曲線Fig.12 Time-history curves of normal stress, shear stress and the damage variable of cohesive elements

為了進一步驗證數(shù)值模擬與物理試驗具有相同的一般性規(guī)律，本文針對第二組初始應(yīng)力組合分別建立4組具有不同初始裂縫長度的模型進行水力劈裂數(shù)值模擬，其中將初始裂縫長度相對模型邊長的百分比定義為裂縫相對長度。除了裂縫相對長度不同外，4組數(shù)值試驗所施加的初始應(yīng)力和約束條件等均相同。最終計算結(jié)果如圖13所示，從圖13中可以看出，劈裂水壓力值隨著初始裂縫相對長度的增加逐漸降低，此規(guī)律與文獻[8]的物理試驗結(jié)果相同。該規(guī)律表明在受到相同外界荷載和邊界條件下，如果有較長的初始裂縫或者更大的缺陷，那么心墻更容易形成水力劈裂。因為初始裂縫越長，裂縫尖端處孔隙水壓力集中程度越劇烈，導(dǎo)致形成更大的水力梯度，促進了裂縫的萌生與擴展。在心墻壩建設(shè)過程中，控制心墻上游面附近初始裂縫或缺陷的大小，將有效降低水力劈裂發(fā)生的幾率。

圖13 初始裂縫相對長度與劈裂水壓力的關(guān)系Fig.13 Relation between initial crack length and fracturing pressure

綜上所述，基于考慮孔隙水壓力內(nèi)聚力模型的心墻水力劈裂數(shù)值模擬，雖然計算出的劈裂水壓力值偏大，但是，能夠反映心墻水力劈裂的一般規(guī)律，并且合理地揭示了水力劈裂的發(fā)生機理。

3 結(jié) 語

(1)本文運用內(nèi)聚力模型模擬了心墻試樣巴西圓盤劈裂試驗，計算得出的拉應(yīng)力—徑向位移曲線與物理實驗基本相符，初步確定了心墻材料的內(nèi)聚力模型力學(xué)參數(shù)。

(2)將考慮孔隙水壓力的內(nèi)聚力模型運用于心墻水力劈裂數(shù)值模擬，雖然由于一些不可避免的原因，與室內(nèi)物理試驗相比，結(jié)果中的劈裂水壓力值大于物理試驗；但是，劈裂水壓力值與豎向初始應(yīng)力的關(guān)系與物理試驗結(jié)果相同，即隨著豎向初始應(yīng)力的增大，劈裂水壓力值線性增大。

(3)通過探究數(shù)值試驗中水壓力在裂縫及其周邊介質(zhì)的分布，可以得出：高水力梯度和裂縫與周邊介質(zhì)的水壓差是水力劈裂裂縫萌生與擴展的兩個重要原因。

(4)提取開裂路徑上一個界面單元的法向應(yīng)力、切向應(yīng)力以及損傷因子，并繪制出歷時曲線。分析曲線的規(guī)律可以發(fā)現(xiàn)，在二維水力劈裂數(shù)值模擬中，水力劈裂的破壞機理是以拉伸破壞為主。

(5)進一步模擬設(shè)置不同初始裂縫長度的水力劈裂試驗，得到如下規(guī)律：初始裂縫長度越大，心墻越容易形成水力劈裂，控制心墻上游面附近初始裂縫或缺陷的大小，將有效降低水力劈裂發(fā)生的幾率。

□

[1] 黃文熙. 對土石壩科研工作的幾點看法[J]. 水利水電技術(shù), 1982,(4):23-28.

[2] 張丙印, 李 娜, 李全明. 土石壩水力劈裂發(fā)生機理及模型試驗研究[J]. 巖土工程學(xué)報, 2005,27(11):1 277-1 281.

[3] 朱俊高, 王俊杰, 張 輝. 土石壩心墻水力劈裂機制研究[J]. 巖土力學(xué), 2007,28(3):487-492.

[4] 殷宗澤，朱俊高，袁俊平，等. 心墻堆石壩的水力劈裂分析[J]. 水利學(xué)報，2006，37(11)：1 348-1 353.

[5] 陳五一, 趙顏輝. 土石壩心墻水力劈裂計算方法研究[J]. 巖石力學(xué)與工程學(xué)報, 2008, 27(7):1 380-1 386.

[6] 李全明，張丙印，于玉貞，等. 土石壩水力劈裂發(fā)生過程的有限元數(shù)值模擬[J]. 巖土工程學(xué)報，2007，29(2)：212-217.

[7] 楊 艷, 周 偉, 常曉林, 等. 高心墻堆石壩心墻水力劈裂的顆粒流模擬[J]. 巖土力學(xué), 2012,33(8):2 513-2 520.

[8] 袁俊平, 王啟貴. 水力劈裂楔劈效應(yīng)試驗研究[J]. 巖土力學(xué), 2015,36(12):3 360-3 364.

[9] 王永剛, 胡劍東, 王禮立. 金屬材料層裂破壞的內(nèi)聚力模型[J]. 固體力學(xué)學(xué)報, 2012,33(5):465-470.

[10] 胡 超, 周 偉, 常曉林, 等. 基于內(nèi)聚力模型的高心墻堆石壩壩頂裂縫模擬及其成因分析[J]. 中南大學(xué)學(xué)報:自然科學(xué)版, 2014,45(7):2 303-2 310.

[11] 馬 剛, 周創(chuàng)兵, 常曉林, 等. 巖石破壞全過程的連續(xù)-離散耦合分析方法[J]. 巖石力學(xué)與工程學(xué)報, 2011,30(12):2 444-2 455.

[12] Munjiza A. The combined finite-discrete element method[M]. Cichester：John Wiley and Son，Inc.，2004：233-243.

[13] Munjiza A, Bangash T, John N W M. The combined finite-discrete element method for structural failure and collapse[J]. Engineering Fracture Mechanics, 2004,71(46)：469-483.

[14] Munjiza A, Owen D R J, Bicanic N. A combined finite-discrete element method in transient dynamics of fracturing solids[J]. Engineering Computations, 1995,12(2)：145-174.

[15] Owen D R J, Feng Y T. The modelling of multi-fracturing solids and particulate media[J]. International Journal for Numerical Methods in Engineering, 2004, 60(1)：317-339.

[16] Zhou Wei, Yuan Wei, Ma Gang, et al. Combined finite-discrete element method modeling of rockslides[J], Engineering Computations, 2016,33(5):1-27.

[17] MA Gang, Zhou Wei, Chang Xiaolin, et al. Combined FEM/DEM modeling of triaxial compression tests for rockfills with polyhedral particles[J]. International Journal of Geomechanics, 2014,14(4):04014014.

[18] Dugdale D S. Yielding of steel sheets containing slits[J]. Journal of the Mechanics and Physics of Solids, 1960,8(2):100-104.

[19] Barenblatt G I. The mathematical theory of equilibrium cracks in brittle fracture[J]. Advances in applied mechanics, 1962,7(1):55-129.

[20] Benzeggagh M L, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus[J]. Composites science and technology, 1996,56(4):439-449.

[21] 王惠敏. 黏性土抗拉特性試驗研究[D]. 南京: 南京大學(xué), 2011.

[22] 龔迪光, 曲占慶, 李建雄, 等. 基于ABAQUS平臺的水力裂縫擴展有限元模擬研究[J]. 巖土力學(xué), 2016,37(5):1 510-1 520.