土石壩心墻裂縫發展的擴展有限元模擬

趙曉龍,卞漢兵,鄭 威,孫兆輝,章賽澤,邱秀梅*

1.河海大學 巖土力學與堤壩工程教育部重點實驗室,江蘇 南京 210098

2.河海大學 巖土工程科學研究所,江蘇 南京 210098

3.山東農業大學 水利土木工程學院,山東 泰安 271018

4.LEM3,CNRS 7239，洛林大學,梅茲 法國 57045

5.徐州市水利建筑設計研究院,江蘇 徐州 221112

土石壩水力劈裂是亟待解決的問題之一[1,2]。目前關于水力劈裂的數值模擬，學者們做了很多工作，如有限元、邊界元和無單元法等，其中有限元應用最為廣泛和成熟。李全明等[3]引入裂縫彌散理論來描述裂縫發展過程，從而實現水力劈裂的有限元模擬。王俊杰[4]等引入四節點等參單元建立裂縫及影響區的有限元模型，實現對裂縫的模擬。上述研究成果對于水力劈裂數值研究有很好的啟發意義，但由于有限元法本身的缺陷，使其模擬裂縫開展遇到很多問題。

有限元在處理裂紋問題時，需要根據裂縫的位置分布網格，在裂尖區域需要分布密集的網格來反映應力的變化，當裂縫擴展時網格需要重新劃分，因此計算效率很低，對于復雜的裂縫擴展問題甚至無法解答[5]。

在此背景下，Belytschko教授于1999年提出了擴展有限元法XFEM[6,7]，用以解決有限元難以解決的不連續問題，如裂紋擴展等。與常規有限元不同，擴展有限元分析裂縫時，網格無需重新劃分，而是在包含裂縫的單元內增加自由度來描述不連續體的位移間斷，位移間斷通過富集函數來描述，可以顯示裂紋的擴展路徑，反映土體在流固耦合作用下變形特性，因而具有優越性[5,8]。

近年來，國內一些學者[9-12]開始將擴展有限元應用到混凝土壩體材料的開裂分析中，但土石壩粘土心墻水力劈裂的擴展有限元模擬鮮有報道。

本文基于擴展有限元，以損傷力學為理論基礎，參照混凝土等材料的建模思路，考慮到基質吸力對黏土力學性能和裂縫發育的影響，建立了壓實黏土的孔隙擴展模型和裂縫張開模型，對裂縫的發展進行了模擬。該思路對于土石壩的水力劈裂數值研究有一定的借鑒性。

1 心墻黏土水力劈裂模型

1.1 黏土孔隙擴展模型

1.1.1 模型本構關系 根據損傷力學理論，將土體內部的孔隙看作微裂紋（細觀裂紋），土體由于微裂紋發展而產生類似損傷的演化，將土體模擬為連續介質，引入損傷變量，以熱力學為基礎，將孔隙的擴展視為能量轉化的過程，由自由能和耗散勢導出孔隙擴展的本構關系和演化方程，描述材料從孔隙擴展到出現宏觀裂紋的過程。

壓實黏土微裂紋演化的動力學特征在張拉和壓縮應力作用下是不同的，為說明其演化對加載路徑的依賴性，使用標量ωc、ωt表示微裂紋在受壓和受拉情況下的演變。構造總體裂紋影響系數ω來表征兩種裂紋對黏土力學性能的影響：

式中αt和(1-αt)用來描述張拉和壓縮產生的微裂紋對材料力學性能總體影響的貢獻率。αt的定義見下式：

于是，壓實黏土的本構關系可以直接通過自由能函數，即方程導出：

通過自由能函數，可推出兩個損傷變量，即在受壓和受拉狀態的關聯共軛熱力學力：

1.1.2 孔隙擴展準則的建立 基于不可逆轉的熱力學過程，孔隙擴展需要一個適當的準則來決定。準則可以是損傷共軛力(見式)的函數。但這種準則難以通過試驗測定，因而現實多采納以物理學為依據的方法。由于微裂紋的形成和張拉應變之間存在固然聯系，引入等效張拉應變的定義，有其中，iε代表i方向的主應變。括號表示只有x為正值才會被考慮，若x為負值或零，則〈x〉=0。將孔隙擴展驅動力定義為等效張拉應變加載歷史中達到的最大值：Yω=max(εeq，Yhis)。Yω取等效張拉應變和加載歷史中達到的最高值Yhis之間的最大值。這表明微裂紋的形成和擴展是一個不可逆的過程。

類似于混凝土材料中的Mazars模型，采用指數形式描述受拉和受壓狀態下微裂紋的演化過程：

式中，參數Bc和Bt用來控制在壓縮和張拉荷載條件下壓實黏土微裂紋的演化速度，很明顯Bt＞Bc，即在壓實黏土中，受拉條件下微裂紋更容易形成和擴展。Yc0和Yt0定義為微裂紋的初始閾值。此外，ωc和ωt的演化法則應當驗證中的能量耗散條件。

1.1.3 基質吸力對壓實黏土孔隙演化的影響 在非飽和土中，還要考慮到負空壓引起的基質吸力對材料力學性能及微裂紋發展的影響。在非飽和條件下，線彈性孔隙介質的本構關系可以表達為：

式中π是等效孔隙壓力，b是比奧系數，由多孔介質理論可以給出：

式中，Kb是材料的不排水體積模量，Km是矩陣體積模量，E是楊氏模量，v是泊松比。楊氏模量E通過函數E(ω)給出，說明比奧系數也受到孔隙擴展的影響。對于壓實黏土，b可以直接取1。非飽和土中，基質吸力對壓實粘土的力學性質及微裂紋的擴展有顯著影響，因此需要將其影響考慮進來。如在土樣失水干燥過程中，可以觀察到干燥裂縫。這主要是由于土樣固有的非均質性和結構上的影響。目前工作中只認為干燥裂縫是由于結構的影響。為此，在損傷演化中，要考慮到由于受限制的收縮變形而引起的局部張拉應力。提出了等效張拉應變的一個新形式?εeq，該形式將力學荷載和基質荷載引起的應變結合起來。由于基質荷載而引起的應變計算公式如下：

式第二項為基質吸力產生的等效拉應變。很明顯，基質吸力引起的等效拉應變與等效孔隙壓力π、比奧系數b、土體的體積模量K0有直接關系，取決于飽和度，因而?εeq受土樣吸水飽和或失水干燥過程的控制。在等效應變的新形式中，式中的損傷驅動力是?εeq而不是εeq的函數。

1.1.4 模型參數的確定 在提出的壓實黏土的力學模型中，有六個參數需要確定：兩個未損傷材料的彈性參數(K0,G0)，及四個損傷特征參數(Bc，Bt，Yt0，Yc0)。參數確定方法如下：初始彈性常數，即初始體積模量K0和剪切模量G0，與初始楊氏模量E0和泊松比v0相關的，可從初始損傷前單軸張拉(或單軸壓縮)試驗應力應變曲線的線性部分獲得。與張拉損傷有關的參數Bt和Yt0和與壓縮損傷有關的參數Bc和Yc0，可分別通過拉伸和壓縮試驗確定。

1.2 壓實黏土的裂縫本構模型

在壓實黏土未形成宏觀裂紋時，假定為均質，此時采用1.1中建立的黏土孔隙擴展模型。在圖1中，應力達到峰值前，采用上述模型，這是典型的非線性過程。當孔隙發展到臨界值，即圖中標注的ω0時，應力也即將達到峰值，宏觀裂縫出現。此時連續介質力學不再適用，需要將不連續體和宏觀裂紋單獨考慮。

在出現宏觀裂紋的應力軟化階段需要引入軟化曲線。以試驗測量結果為依據，給出了一個經驗公式如下：

式中，ft是材料的單軸抗拉強度，[un]m是對應應力為0的最大裂縫寬度，c1是模型參數，用來控制應力—裂紋寬度曲線的形狀。實際中，由于裂紋展開而耗散的斷裂能Gf可以直接通過峰后曲線來獲得，如圖2所示，其大小代表了曲線下的面積。上述三個參數表述了裂紋出現后的力學行為特征，可通過單軸張拉試驗確定。

圖1 巖土材料單軸拉伸應力應變全曲線及孔隙演化過程Fig.1 Complete stress-strain curve of geotechnical material under uniaxial tension and its pore evolution process

圖2 巖土材料宏觀裂縫的應力—裂縫寬度關系Fig.2 Relationship between stress and crack width of macroscopic crack of geotechnical material

連續介質力學和宏觀裂隙的力學行為通過臨界損傷值緊密地聯系在一起，形成一個完整的分析，得到了連續的應力應變全曲線。如圖3所示，在上升曲線部分，為連續介質力學的分析結果，下降曲線則是運用了裂縫張開模型分析的結果。

圖3 巖土材料單軸拉伸試驗曲線及其數值模擬Fig.3 Uniaxial tensile test curve of geotechnical material and its numerical simulation

圖4 帶有裂紋的不連續體Fig.4 Discontinuity with crack

2 裂縫擴展數值方法

2.1XEFM方法

擴展有限元的基本觀點是“擴展”被裂縫分割的單元節點自由度。“擴展”的自由度用來描述不連續體的“位移間斷”。

2.1.1 XEFM位移模式 擴展有限元通過對裂縫分割單元結點增加額外自由度來表征不連續體，位移模式的一般形式：

式中Ω為帶有裂紋的面積域(見圖4)，Γ為裂紋，Ni為單元插值函數，ai為單元節點位移，Hj(x)為Heaviside函數，在裂紋面的上下兩側分別取值+1和-1，作為富集函數用以表示裂紋處位移的間斷，

bi為“擴展”自由度的節點位移。式中第一項即為有限元的公式，而第二項為擴展自由度的貢獻。

2.1.2 XEFM支配方程 利用虛功原理，建立擴展有限元法的支配方程：

式中，da是常規有限元節點位移增量，db是擴展自由度的節點位移增量，為總體勁度矩陣，fext為外界施加在節點上的節點荷載，為由內力引起的常規自由度和擴展自由度的節點荷載列陣。實際上，上述支配方程可以簡寫為KU=F，這一表達和有限元法完全一致。

總體勁度矩陣由Kaa、Kab、Kba和Kbb四個單元勁度矩陣構成，其中與常規有限元單元勁度矩陣一致；考慮了常規和擴展自由度之間的耦合；描述了不連續體對勁度矩陣的貢獻，其中T為裂縫的勁度矩陣，為下標n和t代表裂縫的法線和切線方向。

2.1.3 裂紋積分方法 被裂縫分割的單元，由于跳躍方程在整個方程區域內是不連續的，故在計算單元剛度矩陣和單元節點力時不像有限元法中那樣進行Gauss積分，而是需要進行分區處理。被裂縫分割以后，分區則是不規則的四邊形或三角形。為了解決這一問題，需將裂縫兩邊的區域全部離散為三角形。值得注意的是，這些三角形并不是新生成的單元，而是積分過程中使用的積分域，不形成額外的節點和單元。為了考慮裂縫對單元剛度矩陣和節點荷載的貢獻，裂縫表面也進行了線積分。

2.2 裂縫擴展準則和擴展方向

關于黏土開裂，目前常用的準則主要有拉應力開裂準則和剪切開裂準則。拉應力開裂準則，即當土體中拉主應力達到土體抗拉強度時，土體開裂，裂縫的方向垂直于該拉主應力方向。該準則物理意義明確，但土石壩心墻內部通常都處于受壓狀態，并不存在拉應力。很多學者都認為心墻發生開裂是由于其滿足拉應力開裂準則[13-15]，但即便存在“拱效應”，也很難使心墻的應力為零或負值[16]。剪切開裂準則即莫爾—庫倫強度理論。在這里，還可以采用拉應變開裂準則，即當土中拉應變達到土體抗拉強度對應的拉應變時，土體開裂。相對于拉應力準則，對于心墻土體，拉應變準則更加適合，因為即便土體為壓應力的條件下，拉應變也可能存在。

在本文計算中，極限拉應變、拉應力、抗剪強度、臨界損傷同時作為裂縫開裂和擴展的準則。在每一個加載步后，那些沒有被宏觀裂紋分割的單元內的每個節點的損傷值都會被校核并和開裂準則比較。若開裂準則滿足，那么有兩種被區分的情況。第一種情況，若該單元在裂紋尖端范圍內，這屬于裂紋擴展情況。另一種情況，若該單元周圍是普通的單元，那么在普通單元內會產生新的裂紋尖端。為了便于計算，允許裂紋擴展以確保尖端總在單元的邊界上。同時，將富集函數添加到普通單元節點上以描述單元的位移間斷。隨著新裂紋的形成，在當前加載步下，整個結構的平衡會重新進行計算和校核，直到達到新的平衡。

裂縫擴展方向的確定是宏觀裂紋擴展模型中最難解決的問題之一。在本文，采用了以平均張拉應力為基礎的準則。宏觀裂紋從它的尖端以垂直最大張拉應力的方向進行擴展。考慮到局部相互作用的影響，后期通過非局部平均應力張量來確定。

式中，w是Gauss型權函數，定義如下：

式中，r代表當前點和裂紋尖端之間的距離，lc是材料的特征長度，該特征長度決定了裂紋尖端周圍相互作用域的范圍。lc的數值取決于特征元的尺寸，在二維問題中，可以表示為：其中，Ae是二維網格中的平均單元面積。

3 數值模擬算例

工程背景為蒙陰縣岸堤水庫大壩，壩高29.80 m，頂寬7 m，為混凝土路面。主河槽壩段為粘土心墻砂殼壩。該壩修建于上世紀50年代末，當時建設標準偏低，施工質量較差。心墻局部土料含砂粒較高，滲透系數偏大，且土料整體壓實不均勻，因而存在水力劈裂破壞的風險。

模擬試樣為長20 cm高5 cm的長方形，其邊界條件如圖5所示。為了模擬土石壩心墻的邊界條件，限制水平面前后方向的側向變形，實際上就是一個平面應變問題，在試件的上部施加恒定的初始壓力。在試件的左右兩邊水平位移全部限定，試件下部的豎向位移限定。試件的左邊施加隨時間變化的水壓力，而試件右邊水壓力一直設定為0，試件的上部和下部皆為不透水條件。在試件中央設定一條預設的裂縫，其寬度為0.5 cm，長度為10 cm。庫水壓力的大小和變化可通過調整水壓力的大小來模擬。

整個試件被劃分1600個單元，1701個節點（圖6），預設裂縫為藍色區域部分，用于模擬土石壩中的軟弱滲透面，豎向應力為0.5 MPa，左邊的水壓力在0.05 s內升高到0.2 MPa，然后保持這一水壓力直到試驗結束。土樣力學參數通過三軸試驗，并參考大壩試驗報告相關成果綜合確定，其中粘土的滲透系數為3.99×10-9m/s，砂土部分的滲透系數為1×10-3m/s。因為土樣一直處于恒定壓力狀態下，采用開裂準則為拉應變準則。在裂縫擴展過程中，裂縫穿越單元的滲透系數也隨之增加，采用和砂土一致的滲透系數。這是一個流固耦合加裂縫擴展的復雜過程。

圖5 水力劈裂模擬試件邊界條件Fig.5 Boundary condition of hydraulic fracture simulation specimen

圖6 試件網格劃分Fig.6 Specimen mesh generation

圖7 試件不同階段水壓力分布及裂縫擴展示意圖Fig.7 Distribution of hydraulic pressure and crack extension of specimen in different stages

圖7 給出了對應于不同階段的試件的水壓力分布和變形網格豎向位移的分布。可以看出，在0.5 MPa的豎向壓力下，在0.2 MPa的水壓力沖擊下，土樣被劈裂。在0.05 s時，由于土樣的滲透率較低，試件內部的水壓力還基本保留為初始水壓力。預設滲透弱面的滲透系數遠遠大于試件的滲透系數，其水壓力在加壓的同時達到了0.2 MPa。在這一水壓力沖擊下，預設的滲透弱面被水壓力沖開，呈張開趨勢，如圖7(b.1)示意。使得裂縫前端的單元內部產生較大的拉應變，雖然在該時刻，試件仍然處于受壓狀態，但是預設裂縫的前部單元的拉應變已經達到了裂縫開裂準則規定的拉應變，裂縫開始擴展。

在裂縫擴展過程中，由于裂縫的張開，水壓沿著裂縫向試件深部，即試件右邊擴展，在水壓力作用下，裂縫被撐開，使得裂縫前部單元一直處于拉伸狀態，因此裂縫一直快速向前擴展，如圖7(a.2)和(b.2)示意。在裂縫即將貫通時，裂縫被撐開達到最大程度。在裂縫貫通的一瞬間，由于試件右側處于排水狀態，裂縫中的水壓力，特別是試件右邊，水壓力立即降到0。對應的裂縫同時閉合，隨著時間的增加，裂縫中的水壓力基本保持不變，相應的裂縫開度變化也不大。但是由于裂縫的擴展，土樣與水的接觸面加大，土樣中的水壓力很快達到穩定狀態。

此算例雖然簡單，但卻驗證了水力劈裂發生的兩個基本條件，即“缺陷”和“快速蓄水的初期”[17,18]。缺陷即指心墻中應該存在薄弱面，算例中用滲透系數較大的單元來表示；而快速蓄水指需要一定的水力梯度，心墻內部的水壓力本來很小，幾乎為零，0.2 MPa的水頭很快地進入到薄弱面中，才導致了裂縫的擴展，從而發生了水力劈裂。從圖7中還可以看到，隨著水力劈裂的發展，裂縫周圍土樣內部的孔隙水應力不斷升高，即紅色區域不斷增大，這在一定程度上抵消裂縫和心墻的水頭差，從而減緩裂縫擴展程度。當裂縫貫穿后，裂縫內的高水壓力突然消失，裂縫閉合，土樣內的水壓力重新分布并達到穩定。正是由于水力劈裂后裂縫的快速閉合，才導致實際中的水力劈裂現象難以觀察。算例不僅可以反映裂縫連續的擴展狀態，還能反映裂縫周圍孔隙水應力的連續變化，這是擴展有限元的優越性。

4 結論

本文建立了孔隙擴展模型和裂縫張開模型，并用擴展有限元法模擬水力劈裂過程，其模型參數可通過試驗測得，且模擬過程嚴格服從熱力學的不可逆定律。相比文獻[4]，考慮了心墻中裂縫周圍水壓力隨裂縫擴展過程的變化，也是一種進步。本方法對水力劈裂的數值研究提供了一種思路，但研究成果還有待室內試驗和工程實例的驗證。

[1]張丙印,于玉貞,張建民.高土石壩的若干關鍵技術問題[C].北京：中國土木工程學會第九屆土力學及巖土工程學術會議,2003：163-186

[2]郭沖,邱秀梅,趙曉龍,等.土石壩心墻填土恒壓吸濕變形試驗研究[J].山東農業大學學報：自然科學版,2013,44(1)：76-80

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

[4]王俊杰,朱俊高.堆石壩心墻抗水力劈裂性能研究[J].巖石力學與工程學報,2007,26(增1)：2880-2886

[5]李錄賢,王鐵軍.擴展有限元法XFEM及其應用[J].力學進展,2005,35(1)：5-20

[6]Belytschko T,Black T.Elastic crack growth in finite elements with minimal remeshing[J].International Journal for Numerical Methods in Engineering,1999,45：601-620

[7]Mo?s N,Dolbow J,Belytschko T.A finite element method for crack growth without remeshing[J].International Journal for Numerical Methods in Engineering,1999,46：131-150

[8]郭歷倫,陳忠富,羅景潤,等.擴展有限元方法及應用綜述[J].力學季刊,2011,32(4)：612-625

[9]杜效鵠,段云嶺,王光綸.重力壩斷裂數值分析研究[J].水利學報,2005,36(9)：1035-1042

[10]方修君,金 峰,王進廷.基于擴展有限元法的Koyna重力壩地震開裂過程模擬[J].清華大學學報：自然科學版,2008,48(12)：2065-2069

[11]方修君,金 峰.裂隙水流與混凝土開裂相互作用的耦合模型[J].水利學報,2007,38(12)：1466-1474

[12]董玉文,任青文.重力壩水力劈裂分析的擴展有限元法[J].水利學報,2011,42(11)：1361-1367

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

[14]孫亞平.水力劈裂機理研究[D].北京：清華大學,1985

[15]馮曉瑩,徐澤平.心墻水力劈裂機理的離心模型試驗研究[J].水利學報,2009,40(10)：1259-1263

[16]王俊杰,朱俊高,張 輝.關于土石壩心墻水力劈裂研究的一些思考[J].巖石力學與工程學報,2005,24(增2)：5664-5668

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

[18]趙曉龍,邱秀梅,韓慧敏,等.土石壩帶裂縫黏土心墻破壞機理試驗研究[J].中國農村水利水電,2016(2)：134-138