基于應力-滲流-損傷耦合模型的重力壩三維水力劈裂數(shù)值模擬

沙 莎，張國新

（1.中國電建集團 北京勘測設計研究院有限公司，北京 100024；2.流域水循環(huán)模擬與調控國家重點實驗室 中國水利水電科學研究院，北京 100038）

1 研究背景

隨著我國水資源開發(fā)需求的增長以及筑壩技術的發(fā)展，修建200～300 m的超高壩早已成為現(xiàn)實，混凝土高壩的安全問題成為壩工界一個備受關注的重要課題。由于混凝土材料抗拉強度較低，大體積混凝土壩施工期極容易產生裂縫，尤其是碾壓混凝土壩，還存在薄弱層面，這些裂縫和薄弱面在庫水壓力作用下是否會發(fā)生高壓水劈裂是一個亟需研究的問題。國內外已有幾座重力壩出現(xiàn)了高壓水劈裂的現(xiàn)象，例如：加拿大的雷威爾斯托克壩、美國的德沃歇克壩和盧塞爾壩，施工過程中雖然采用了預冷骨料、水管冷卻、表面保溫等綜合溫度控制措施，壩體上游面仍出現(xiàn)了表面裂縫，水庫蓄水后經過一段時間，表面裂縫突然大范圍擴大，發(fā)展為劈頭裂縫，深入壩內幾十米，有的甚至將整個壩段一分為二，引起嚴重漏水［1］。我國柘溪［2］、桓仁［3］大頭壩，玉石［4］、觀音巖碾壓混凝土重力壩也出現(xiàn)了同樣的問題。因此，對混凝土重力高壩或超高壩進行安全評估時，考慮高壓水劈裂作用的影響是非常重要的。

目前，在混凝土高壓水劈裂的試驗研究方面，已有學者進行了相關研究。Brühwiler和Saouma［5-6］研究了不同級配的混凝土裂縫內的水壓分布變化規(guī)律及裂縫中靜水壓力對混凝土表觀斷裂能、表觀斷裂韌度的影響。Slowik和Saouma［7］研究了裂縫邊緣張開速度不同時裂縫內水壓分布隨時間的變化以及裂縫突然閉合時對水壓分布的影響。王建敏［8］采用大型尺寸楔入式緊湊拉伸試件進行了水力劈裂試驗，研究了裂縫內的靜水壓力對斷裂性能的影響，指出隨著水壓力的增加，試件所能承受的最大荷載減小。賈金生等［9］提出了全級配混凝土試件單軸拉、壓應力作用下高壓水劈裂模擬試驗新方法，并推導了判定重力壩壩踵是否會發(fā)生水力劈裂的分析公式。通過試驗可以揭示水力劈裂的機理，但是由于試驗中采取的試件均是立方體、圓柱體試件，所以試驗不能真實反映混凝土壩的形態(tài)及其受力條件，也無法預測裂縫擴展范圍及危險程度。鑒于試驗的不足，很多學者致力于混凝土壩高壓水劈裂的數(shù)值研究。陳勝宏等［10］采用有限元法研究了高壓水劈裂對小灣高拱壩壩踵開裂的影響，指出考慮高壓水劈裂壩踵裂縫擴展范圍比不考慮時大。Barpi和Valente［11］對混凝土重力壩建基面處裂紋的水力劈裂進行了模擬研究，討論了斷裂過程區(qū)的發(fā)展對裂縫內水壓的影響。董玉文等［12］采用擴展有限元法進行了向家壩壩踵水力劈裂的數(shù)值分析，將裂縫面上的水壓力視為均布荷載。王克峰等［13］采用擴展有限元法研究了考慮流固耦合效應的重力壩的水力劈裂問題。Secchi and Schrefler［14］基于分離裂縫模型和三維自動網(wǎng)格重構技術提出一種模擬飽和孔隙介質三維水力劈裂的方法，并且用于混凝土重力壩在水頭超載作用下壩踵的開裂分析。

水力劈裂是個復雜問題，涉及到四個耦合過程：（1）裂紋面上水壓力引發(fā)的介質變形；（2）裂紋內及周圍介質中流體流動；（3）裂紋的擴展；（4）裂紋內的水滲透到周圍介質中［15］。上述混凝土壩高壓水劈裂的數(shù)值研究中，只有文獻［14］考慮了這4個耦合過程。文獻［10］和［12］將裂縫面上的水壓力視為均布荷載，沒有考慮裂縫擴展過程中水體流動引起的裂縫內的水壓重分布。文獻［11］和［13］則假定裂縫面是不透水的，且沒有考慮周圍介質中流體流動。已有的重力壩水力劈裂的數(shù)值研究大多是二維的，三維方面的十分罕見，特別是重力壩劈頭裂縫方面的。

目前模擬水力劈裂的數(shù)值方法有：相場法［16］，近場動力學［17］，元胞自動機法［18］，離散元法［19］，流形元法［20］，無單元法［21］，有限元法［14］，邊界元法［22］等。有限元法是最成熟，應用最廣的數(shù)值方法，其本質上是一種連續(xù)介質力學方法，必須加以改進才能被用來分析裂縫擴展這樣的不連續(xù)問題，改進的方法分為兩類：變網(wǎng)格法和固定網(wǎng)格法。變網(wǎng)格法隨著裂縫的產生，需要不斷地重新劃分網(wǎng)格，裂紋面必須與單元邊保持一致，在裂尖區(qū)域需要細化單元或改用特殊形式的縫端奇異元，產生極大的前處理工作量。固定網(wǎng)格法則是保持網(wǎng)格不變，通過修改開裂單元插值關系、本構關系反映裂縫的存在，比如擴展有限元（XFEM）、連續(xù)損傷模型。相比較而言，固定網(wǎng)格法的應用更為方便。本文基于有限元采用一種各向同性損傷模型進行水力劈裂的模擬，將混凝土視為飽和的孔隙介質，根據(jù)孔隙介質有效應力原理，提出一種應力-滲流-損傷耦合模型，考慮水力劈裂過程中的耦合效應。采用該耦合模型，對國內某重力壩上游面垂直裂縫水力劈裂進行數(shù)值模擬。

2 應力-滲流-損傷耦合模型

假定大壩混凝土是飽和的孔隙介質。應力-滲流-損傷耦合模型包括了以下非線性行為：（1）單元損傷后的應力-應變關系；（2）損傷引起的孔隙水壓影響系數(shù)的改變；（3）應力、損傷引起的單元滲透系數(shù)的改變。計算中假定受壓時，應力應變關系是線彈性的，因為重力壩壓應力一般不會超過其抗壓強度。以下介紹這種耦合模型的特性。

2.1 損傷模型 混凝土的斷裂特性表現(xiàn)為在真實的縫端前面存在一個斷裂過程區(qū)，在斷裂過程區(qū)中，通過骨料的咬合和界面間的摩擦力，相互之間仍有應力作用，斷裂過程區(qū)的存在使得混凝土呈現(xiàn)出應變軟化特性。由Hillerborg等［23］提出的虛擬裂縫模型可以很好地模擬斷裂過程區(qū)的力學特性。目前，關于混凝土拉伸曲線的數(shù)學表達式，大多數(shù)學者主張上升段采用直線，主要區(qū)別在下降段，包括：單直線下降、分段下降、曲線下降。但無論采用何種形式，均應保持拉伸曲線的斷裂能相同［24］。本文下降段采用由江見鯨［24］提出的指數(shù)下降型表達式。因此，單軸受拉條件下，應力-應變關系表達為：

式中：E0為材料無損時的彈性模量；ft0為抗拉強度；εt為單元的單軸拉伸應變；ε0為開裂時的應變，ε0=ft0/E0；α為控制下降段的軟化系數(shù)。

斷裂能GF的含義如圖1所示，為應力-裂縫寬度曲線與坐標軸圍成的面積，可表達為式（2）；單位開裂寬度的斷裂能gf的含義如圖2中陰影所示，為應力-應變曲線下降段與ε=ε0、橫坐標軸ε圍成的面積，可表達為式（3）。

圖1 應力-裂縫寬度曲線和斷裂能Gf

圖2 單位開裂寬度的斷裂能gf

式中：w為過程區(qū)內所有微裂縫張開的位移量之和，h為混凝土微裂縫分布區(qū)寬度。

根據(jù)式（1）

斷裂能GF與單位開裂寬度的斷裂能gf存在如下關系：

式中：lt為混凝土單元特征尺寸。對于平面單元，取積分點區(qū)域面積的平方根；對于實體單元，取積分點區(qū)域體積的立方根。

由式（4）、式（5），本構曲線下降段的軟化系數(shù)應滿足

根據(jù)Lemaitre［25］提出的等效應變原理，假定Cauchy應力σ作用在受損材料上的應變與有效應力作用在無損材料上的應變等價，即

得到

式中：E0為材料無損時的彈性模量；為材料受損后的彈性模量；d為損傷變量，d=0對應于無損傷狀態(tài)，d=1對應于完全損傷狀態(tài)，0＜d＜1對應于不同程度的損傷狀態(tài)。

從損傷力學的角度看，巖石、混凝土宏觀應力-應變的非線性，是由于其受力后不斷損傷引起的微裂紋萌生和擴展造成的。這種脆性在拉伸作用下更加明顯，因此采用彈性的損傷力學本構方程來描述其力學特性是合適的。經過研究和實驗證明，采用彈塑性損傷本構模型模擬的結果和彈性損傷本構模型的結果相差不大［26］。因此，本文采用彈性損傷模型進行計算。

由式（1）、式（8），單軸拉伸作用下，混凝土損傷演化方程可表示為：

2.2 孔隙介質有效應力原理 有效應力概念應用于被單相流體浸潤的巖石、混凝土材料，可以看成太沙基有效應力原理應用于土體的推廣。基于有效應力原理，Biot首次引進一個標量參數(shù)，即Biot系數(shù)，反映孔隙水壓對有效應力的影響（Biot，1941，1955，1977）［27-29］。在文獻［28］中，孔隙率 bi被定義為有效孔隙面積與橫截面Ai的比值。有效孔隙面積則被定義為在垂直于橫截面Ai的方向，單位長度所有微小孔隙面積的總和?？紫堵蔮i同樣可以代表孔隙體積（Vp）與代表體體積（Vb）的比值。在混凝土壩彈性分析中，通常假定含有孔隙的大壩混凝土是均質、各向同性的，因而經常采用bi的一個各向同性值bx=by=bz=b0進行混凝土壩的彈性分析。對于土體的應力計算，通常取b0≈1。然而，對于彈性范圍內混凝土的應力計算，取b0≈1則是不合理的，因為混凝土即使在即將破壞的狀態(tài)下，b0仍小于1。

圖3 飽和孔隙介質的應力分解

文中應力符號以拉應力為正。飽和孔隙介質的應力分解如圖3所示，作用在單元表面的總應力可以分解成兩部分：與內部水壓力p相互平衡的外力bp和有效應力，單元的平均變形只與有效應力有關，可表示為：

2.3 孔隙水壓影響系數(shù)與損傷關系 單元未損傷時假定bx=by=bz=b0，單元損傷后，垂直于裂縫方向的孔隙水壓影響系數(shù)與平行于裂縫方向孔隙水壓影響系數(shù)理應是不一樣的。圖4定性的展示了開裂對孔隙水壓影響系數(shù)的影響，1、2、3分別代表第一、二、三主應力的方向，建立以主應力方向為坐標軸的局部坐標系。開裂對垂直于裂縫面方向的b1影響較大，b2和b3在一個范圍內，即b0＜b2=b3＜b1［30］。Bary and Bournazel［31］通過實驗得出了一個各向異性的Biot張量，其與各個方向的損傷值和孔隙水壓有關。由于缺乏開裂對平行于裂縫方向的孔隙水壓影響系數(shù)影響的試驗資料，本文假定混凝土損傷之后，孔隙水壓影響系數(shù)仍各向同性，即b1=b2=b3。彈性狀態(tài)下，b1=b2=b3=b0，b0為初始Biot系數(shù)；完全損傷后，b1=b2=b3=1，孔隙水壓影響系數(shù)隨損傷的演化方程可表示為：

2.4 滲透系數(shù)與應力及損傷關系 混凝土在微觀上是由骨架和孔隙組成的，這種構造特征使得其在受荷載或擾動作用后，其微觀幾何形態(tài)、孔隙的結構發(fā)生改變，從而導致孔隙率和滲透性發(fā)生改變?？紫堵实淖兓饕蓛刹糠纸M成：一是由于結構形變引起孔隙體積變化；二是微裂紋等缺陷的萌生、擴展及貫通，使得材料內部的孔隙結構和大小發(fā)生改變，即材料發(fā)生損傷而引起的孔隙結構變化［32］。

當單元未發(fā)生拉損傷，處于彈性狀態(tài)，即損傷變量d=0時，滲透系數(shù)與有效應力成指數(shù)關系，表示如下：

圖4 開裂對孔隙水壓影響系數(shù)及滲透性的影響

式中：k0為初始滲透系數(shù)；a0為耦合系數(shù)；σ1、σ2、σ3為有效主應力，以受拉為正；k11、k22、k33為方向與有效主應力一致的主滲透系數(shù)。

當單元第一主應力達到抗拉強度，發(fā)生損傷，即損傷變量d＞0時，單元中產生裂縫，圖4定性地展示了開裂對單元滲透性的影響?；炷林械乃饕刂芽p流動，方向2、方向3的滲透系數(shù)大大提高，此時近似按照裂隙滲流考慮，采用在裂隙滲流研究中運用得最為廣泛的單裂隙平行板水力模型，則滲透系數(shù)表示如下：

式中：u為裂縫張開的寬度，近似按u=（ε1-εt）lt計算，εt、lt的定義見2.1節(jié)；υ為水的運動黏滯系數(shù)，g為重力加速度。由k11、k22、k33形成局部坐標系下的滲透矩陣［k′］，經坐標轉換將其轉換到整體坐標系中，則整體坐標系下的滲透矩陣［k］可表示為：

2.5 滲流場基本微分方程 假定水不可壓縮，根據(jù)Darcy定律，可得無源非穩(wěn)定情況下滲流連續(xù)方程：

式中：S為飽和度，H為水頭勢，kij為滲透系數(shù)張量；

對于非恒定滲流，其定解條件包括如下邊界條件和初始條件：

水頭邊界：

流量邊界：

隔水邊界：

初始條件：

式中：f為已知水頭邊界值，q為已知流量邊界，H0為已知函數(shù)。

如果按照恒定滲流計算，即式（15）的右側為0，沒有初始條件式（19）。

2.6 程序實現(xiàn) SAPTIS軟件是中國水利水電科學研究院結構材料所獨立開發(fā)的大型結構多場仿真與非線性分析軟件，該軟件具備混凝土壩溫度、滲流、變形、應力等多場耦合仿真分析功能［33-35］。本文在SAPTIS原有功能基礎上加入了應力-滲流-損傷耦合非線性功能，進行水力劈裂的數(shù)值模擬計算。

3 實驗驗證

為了驗證本文應力-滲流-損傷耦合模型的正確性，對賈金生等［9］做的圓柱體混凝土試件水力劈裂試驗進行了數(shù)值模擬，試件尺寸及受力條件如圖5所示，本文計算模型如圖6所示。選取文獻［9］中3個試件進行數(shù)值模擬，試件力學參數(shù)、受力狀態(tài)與文獻［9］一致，見表1。水壓加載過程與文獻［9］一致，如圖7所示，即采用梯級加載方式，水壓小于1.0 MPa時，每小時增加0.5 MPa；水壓1 MPa至2 MPa時，每小時增加0.2 MPa，水壓大于2 MPa時，每小時增加0.1 MPa。滲流場計算按照非恒定滲流考慮，混凝土滲透系數(shù)k0=5×10-9m/s，耦合系數(shù)a0=0.01。

圖5 高壓水劈裂試驗原理（單位：mm）

圖6 計算模型

表1 圓柱體混凝土試件力學參數(shù)及受力狀態(tài)

試件Ⅳ發(fā)生水力劈裂破壞時的損傷分布如圖8所示，由圖8可以看出，試件損傷區(qū)域幾乎和預制縫平行。不同水壓作用下，試件Ⅳ預制縫所在截面損傷分布如圖9所示，孔隙水壓分布如圖10所示。由圖9可知，水壓1.4 MPa，縫端開始損傷，隨著水壓的增加，損傷區(qū)域越來越大，水壓3.4 MPa，損傷區(qū)域面積占總面積55%，與文獻［9］描述的劈裂面積十分接近，水壓超過3.4 MPa，截面全部損傷，本文認為水壓3.4 MPa時，試件發(fā)生了水力劈裂破壞。由圖10可以看出，隨著水壓增加，最大孔隙水壓分布區(qū)域逐漸增大，其邊界與圖9中損傷區(qū)域的邊界十分吻合，這是由混凝土單元損傷之后，滲透系數(shù)顯著增加，水不斷滲入引起的。不同試件水力劈裂水壓的試驗值、計算值見表2，計算值與試驗值吻合很好，可以看出，本文耦合模型可以很好地模擬水力劈裂現(xiàn)象。

4 國內某重力壩水力劈裂數(shù)值模擬

圖7 水壓加載過程

圖8 P=3.4MPa，試件損傷分布

圖9 不同水壓作用下，預制縫所在截面損傷分布等值線圖

該水電站為西南某水電基地的中游河段，水庫正常蓄水位1134 m，庫容約20.72億m3。施工過程中，16號壩段上游面出現(xiàn)了表面裂縫，蓄水過程中這些表面裂縫擴展為劈頭裂縫。本文以16號壩段為研究對象，在其上游面設置了一個豎直的初始縫面，模擬了其在蓄水過程中的水力劈裂現(xiàn)象。該壩段壩頂高程1139 m，建基面高程990 m，壩體及壩體-地基模型如圖11所示，模型中考慮了壩體排水、廊道、地基排水、帷幕。由于缺乏施工期形成的表面裂縫的形狀、范圍，本文假定初始裂縫面形狀為三角形，底高程1005 m，頂高程1025 m，最大深度為5 m。壩基的計算范圍為向上游、下游和基巖深度方向各取1.5倍壩高，壩基底邊界為豎向約束，四個側邊界均為水平約束。

計算中考慮的荷載為自重、上下游水荷載，先施加自重荷載，再施加水壓力，上游水位按每荷載步增加10 m的方式從1020 m抬升至1120 m，最后一個荷載步水位從1120 m抬升至正常蓄水位1134 m，下游水位1010 m保持不變。壩體混凝土、地基力學特性見表3。滲流場計算按照恒定滲流考慮，混凝土滲透系數(shù)k0=5×10-9m/s，耦合系數(shù)a0=0.01。

不同水位壩體損傷分布如圖12所示。由圖12可以看出，（1）蓄水至1070 m時，壩體出現(xiàn)損傷；（2）隨著水位增加，損傷區(qū)域從裂縫面尖端逐步向壩體上部、內部、下部擴展；（3）水位從1070 m上升至1110 m過程中，損傷區(qū)域緩慢擴展，從1110 m上升至1120 m，損傷區(qū)域明顯增大，向內部擴展至壩體排水管處；（4）水位由1120 m上升至正常蓄水位1134 m時，損傷區(qū)域顯著增大，上部達到高程1060 m處，下部達到高程992 m處，最大裂縫深度33 m，最終裂縫面呈圓弧狀。

圖10 不同水壓作用下，預制縫所在截面孔隙水壓分布等值線圖（單位：m）

表2 水力劈裂水壓計算值、試驗值對比

圖11 壩體、壩體-地基模型

表3 國內某重力壩材料參數(shù)

圖12 不同水位壩體損傷分布等值線

設計院根據(jù)安全監(jiān)測實測結果得到的正常蓄水位下的裂縫范圍如圖13所示。由圖13可以看出，正常蓄水位下，裂縫范圍上部達到高程1070 m處，底部達到高程997 m處，裂縫向壩體內部擴展的最大深度達到40 m左右。比較圖12（f）及圖13可以看出，本文數(shù)值模擬得到的正常水位下裂縫的擴展范圍與設計院根據(jù)安全監(jiān)測實測結果得到的范圍基本吻合，設計院的裂縫范圍比本文向內部擴展的更深一點，向上擴展的更高一點。經過初步分析差異主要與初始裂縫面形狀、范圍，壩體力學特性參數(shù)，蓄水過程等有關，計算中當初始裂縫面的最大深度增加時，正常水位下裂縫擴展范圍增大，壩體抗拉強度降低時，正常水位下裂縫擴展范圍也增大，所以掌握真實的初始裂縫面形狀、范圍，壩體力學特性參數(shù)，蓄水過程等資料是采用本文耦合模型準確預測裂縫是否擴展、擴展過程、最終擴展范圍的關鍵。

圖13 裂縫預估范圍

5 結論

本文提出了一種應力-滲流-損傷耦合模型，并用于混凝土重力壩三維水力劈裂的模擬。

該耦合模型具有以下特點：（1）考慮了混凝土斷裂過程區(qū)的應變軟化特性；（2）基于斷裂能守恒原理將損傷模型與斷裂力學相結合，使得斷裂能的消散不會受到網(wǎng)格大小的影響；（3）考慮了水力劈裂過程中的4個耦合過程。

采用該耦合模型，對內置裂縫的圓柱體試件的水力劈裂進行了數(shù)值模擬，水力劈裂水壓的計算值與試驗值吻合很好，驗證了所提耦合模型的合理性。

在此基礎上，進行了國內某混凝土重力壩三維水力劈裂的模擬：（1）蓄水至1070 m時，壩體出現(xiàn)損傷；（2）隨著水位增加，損傷區(qū)域從裂縫面尖端逐步向壩體上部、內部、下部擴展，蓄水至正常蓄水位1134 m，損傷區(qū)域上部達到高程1060 m處，下部達到高程992 m處，最大裂縫深度33 m，最終裂縫面呈圓弧狀；（3）數(shù)值模擬得到的損傷區(qū)域與設計院根據(jù)安全監(jiān)測實測結果得到的范圍基本吻合。研究結果表明：采用該耦合模型可以十分便捷有效地進行混凝土重力壩三維水力劈裂的模擬，評估表面裂縫的危險程度。然而水力劈裂是一個復雜的問題，損傷對孔隙水壓影響系數(shù)、滲透系數(shù)的影響還需要進一步研究。