重力壩壩踵裂縫水力劈裂耦合分析

鄭安興，潘國勇

(1.浙江水利水電學(xué)院 水利與環(huán)境工程學(xué)院，浙江 杭州 310018； 2.杭州市富陽區(qū)水利水電工程質(zhì)量安全服務(wù)保障中心，浙江 杭州 311400)

壩踵裂縫在高水頭壓力作用下有可能進一步擴展，裂縫擴展會嚴重降低重力壩的承載力，威脅大壩安全[1]。目前從理論上解釋重力壩壩踵裂縫水力劈裂耦合作用是非常困難的，然而數(shù)值模擬方法為壩踵裂縫水力劈裂耦合分析提供了有效途徑[2-3]。其中擴展有限元法（XFEM）是目前解決斷裂等不連續(xù)問題較為有效的數(shù)值模擬方法[4]。該方法處理斷裂問題時結(jié)構(gòu)內(nèi)部幾何或物理邊界和計算網(wǎng)格獨立，因而能方便分析水力劈裂問題。Wang 等[5]利用擴展有限元和有限體積法相結(jié)合的形式研究了混凝土重力壩的水壓致裂機理；Ren 等[6]基于非耦合模型用擴展有限元模擬了混凝土中的水力劈裂問題；Paul 等[7]建立基于黏聚區(qū)模型的三維耦合HM-XFEM，并將其應(yīng)用于非平面水力裂縫傳播和多重水力裂縫干擾中；Shi 等[8]建立了模擬三維非平面流體驅(qū)動裂縫擴展的全耦合流體力學(xué)XFEM 模型；Zheng 等[9]基于擴展有限元法研究了水力裂縫與天然裂縫的相互作用；董玉文等[10]采用擴展有限元法研究了混凝土重力壩水力劈裂問題；Tan 等[11]采用基于擴展有限元的CZM 方法，研究了具有過渡帶的橫觀各向同性層狀頁巖地層中水力裂縫垂向擴展行為。

綜上可見，采用擴展有限元法對混凝土水力劈裂問題進行了廣泛研究，而重力壩水力劈裂機理研究相對較少。本文考慮裂縫內(nèi)的水壓力作用，建立重力壩壩踵裂縫水力劈裂耦合分析的擴展有限元模擬方法，進行重力壩壩踵裂縫水力劈裂耦合分析。

1 擴展有限元法的基本理論

依據(jù)單位分解方法，在傳統(tǒng)有限元近似位移函數(shù)式中加入反映裂縫尖端奇異性和裂縫面不連續(xù)性的附加函數(shù)，得到擴展有限元法的近似位移形式[12]如下：

式中：Ni為節(jié)點形函數(shù)； ?、 ?Γ與 ?Λ分別為區(qū)域內(nèi)所有離散節(jié)點的集合、裂縫貫穿單元的節(jié)點集及裂縫尖端單元的節(jié)點集；ui、ai與分別為節(jié)點自由度向量、裂縫貫穿單元的附加自由度向量及裂尖單元的附加自由度向量；H(x)為階躍函數(shù)，反映裂縫兩側(cè)位移的不連續(xù)性，可表示為：

Fl(x)為裂縫尖端附加函數(shù)，反映裂縫尖端奇異性，對各向同性彈性體而言，其表達式[13]為：

式中：r與θ為裂尖局部極坐標，如圖1 所示。

圖1 裂尖局部極坐標Fig.1 Local polar coordinates of crack tip

對于裂縫貫穿單元和裂縫尖端單元，裂縫面上同一點的相對位移可由式（1）求得：

式中：w為裂縫面上同一點的相對位移。

擴展有限元法的近似位移形式構(gòu)造后，通過虛功原理推導(dǎo)其控制方程：

式中： σ為柯西應(yīng)力張量； δu為虛位移； δε為相應(yīng)的虛應(yīng)變；b、t和p分別為計算域 ?上的體力、邊界 Γt上的外力和裂縫面上 Γc的水壓力。

由擴展有限元近似位移函數(shù)式（1）和平衡方程（5），聯(lián)合式（4），可得到擴展有限元的控制方程：

式中：K為整體勁度矩陣，由單元勁度矩陣集合而成：

fi為整體荷載向量，可表示為：

式中：l=1，2，3，4。

d為節(jié)點位移未知量向量，可表示為：

2 裂縫斷裂準則及裂縫應(yīng)力強度因子

裂縫斷裂準則采用最大周向拉應(yīng)力強度因子準則，裂縫開裂角 θc由下式確定：

式中：KI為I 型應(yīng)力強度因子；KII為II 型應(yīng)力強度因子。

裂縫斷裂判據(jù)為：

式中：KIeq為周向拉應(yīng)力強度因子臨界值。

裂縫應(yīng)力強度因子計算采用相互作用積分方法?？紤]兩種應(yīng)力狀態(tài)，狀態(tài)為真實狀態(tài)，取擴展有限元數(shù)值解；狀態(tài)為輔助狀態(tài)，取裂縫尖端的漸近場[11]。狀態(tài)1、2 的相互作用積分I(1,2)表達式如下：

式中：Q為權(quán)函數(shù)； δ1j為科羅內(nèi)克爾符號；A為裂縫尖端周圍的積分區(qū)域；pj為裂縫面上的水壓力。

相互作用應(yīng)變能密度W(1,2)為：

狀態(tài)1、2 的相互作用積分和應(yīng)力強度因子之間存在如下關(guān)系：

式中：對于平面應(yīng)力，E?=E；對于平面應(yīng)變，E*=E/(1－v2)，v為泊松比。

3 水力劈裂耦合數(shù)值模型

假定水流在裂縫內(nèi)的運動滿足層流條件，裂縫內(nèi)水流的流量q應(yīng)滿足立方定律：

式中： μ為水流的黏性系數(shù)；w為裂縫的水力開度；p為水壓力。

對于低滲透介質(zhì)體，可以忽略裂縫面上的濾失量，裂縫中水流還應(yīng)滿足質(zhì)量守恒方程：

根據(jù)式（16）和（17），通過控制體積法推導(dǎo)得到裂縫內(nèi)水壓的差分形式：

圖2 為裂縫內(nèi)水流運動模型，控制體積內(nèi)的四邊形單元由裂縫上的高斯點概化而來，結(jié)合式（18）～（20）按順序聯(lián)立求解每個流體單元，可得到整條裂縫內(nèi)的水壓力。

圖2 裂縫內(nèi)水流分析模型Fig.2 Analysis model of flow in crack

4 試件水力劈裂驗證

試件的幾何尺寸如圖3 所示，裂縫的初始長度為1 m，斷裂韌度為3.8×107N/m3/2，彈性模量E=40 MPa，泊松比ν=0.2，試件的右側(cè)邊界為固定約束，初始裂縫邊界水壓力為3 MPa，水力劈裂計算時不考慮試件自重。

采用四邊形等參單元進行擴展有限元分析，該計算區(qū)域共劃分625 個單元，676 個節(jié)點，有限元網(wǎng)格劃分如圖3 所示。水力劈裂作用下裂縫擴展步數(shù)共5 步，裂縫擴展步長為0.4 m。不同裂縫擴展長度時裂縫內(nèi)水壓力分布如圖4 所示。可見，在裂縫擴展前，裂縫內(nèi)水壓力基本與邊界水壓力相同，當裂縫開始擴展時，裂縫內(nèi)水壓會產(chǎn)生劇烈下降，隨著裂縫不斷擴展，裂縫內(nèi)水壓力又會發(fā)展成邊界全水頭。上述計算結(jié)果與Brühwiler 等[15]得到的試驗結(jié)論一致，從而驗證了計算程序的準確性。

圖3 計算網(wǎng)格(單位：m)Fig.3 Computational grid (unit: m)

圖4 不同裂縫長度下裂縫內(nèi)水壓力Fig.4 Water pressures in crack under different crack lengths

5 重力壩水力劈裂耦合分析

某混凝土重力壩的壩高、壩頂寬和壩底寬分別為192.0、10.0 和159.8 m。考慮最不利工況，假定上游水頭192.0 m，上游坡面坡度為3%，下游坡面坡度為70%。壩踵處有一初始裂縫，裂縫長度為10 m。壩基底部和壩基左右邊界為法向約束，采用四邊形等參單元進行擴展有限元分析，計算模型共劃分5 920 個單元，6 109 個節(jié)點及有限元網(wǎng)格劃分如圖5 所示。壩體和壩基的材料計算參數(shù)分別為：重度24 和28 kN/m3，彈性模量22 和24 GPa，泊松比0.167 和0.200，斷裂韌度21 500 和25 800 kN/m3/2。

圖5 重力壩幾何尺寸及裂縫擴展路徑（單位：m）Fig.5 Geometric dimensions and crack propagation path of gravity dam (unit: m)

重力壩所受荷載主要有自重與水壓力。水壓力考慮重度超載方式，保持水位在192 m，只提高水重度，水壓增加前后比值為超載系數(shù)[14]，本算例的重度超載系數(shù)取為2.5。

考慮裂縫張開寬度和裂縫面水壓力之間的耦合關(guān)系，應(yīng)用最大周向拉應(yīng)力斷裂準則，裂縫每步擴展步長取為3 m，裂縫擴展分9 步。采用擴展有限元方法計算得到的考慮水力劈裂耦合作用的重力壩裂縫擴展如圖5 所示。從圖5 中可以看出，壩踵初始逐漸向壩基底部擴展，且裂縫擴展朝下游方向。

為了分析裂縫擴展對重力壩應(yīng)力分布的影響，由擴展有限元法計算得到的考慮水力劈裂耦合效應(yīng)的重力壩應(yīng)力分布見圖6。由圖6 可知，裂縫尖端區(qū)域和壩踵附近都出現(xiàn)應(yīng)力集中現(xiàn)象，且裂縫尖端區(qū)域的主應(yīng)力值明顯高于壩踵附近的主應(yīng)力值。

圖6 應(yīng)力云圖Fig.6 Stress distribution diagram

圖7 為裂縫張開寬度與裂縫長度之間的關(guān)系，可見，裂縫張開寬度隨著裂縫長度的增加而增大。不同裂縫長度下的裂縫內(nèi)水壓力分布見圖8?？梢?，在裂縫擴展前，裂縫內(nèi)水壓力基本與邊界水壓力相同；當裂縫開始擴展時，裂縫內(nèi)水壓下降；而后裂縫張開寬度不斷增大，裂縫內(nèi)水壓力又變成邊界全水頭。

圖7 不同裂縫長度下的裂縫張開寬度Fig.7 Crack opening widths under different crack lengths

圖8 不同裂縫長度下的裂縫內(nèi)水壓力Fig.8 Hydraulic pressures along crack face under different crack lengths

不考慮水力劈裂（裂縫內(nèi)無水壓）、不考慮水力劈裂耦合（水壓為邊界全水頭）和考慮水力劈裂耦合作用下的裂縫擴展路徑如圖9 所示?？梢?，有、無水力劈裂作用下裂縫都向壩基底部擴展，無水力劈裂作用下的裂縫開裂角大于水力劈裂作用下的，無水力劈裂耦合作用下的裂縫開裂角小于水力劈裂耦合作用下的。

圖9 有無水力劈裂作用下裂縫擴展路徑Fig.9 Crack propagation path with or without hydraulic fracturing

不考慮水力劈裂（裂縫內(nèi)無水壓）、不考慮水力劈裂耦合（水壓為邊界全水頭）和考慮水力劈裂耦合作用下的裂尖應(yīng)力強度因子隨裂縫擴展長度的變化規(guī)律見圖10?？梢?，不考慮裂縫水力劈裂時的Ⅰ型應(yīng)力強度因子曲線比較光滑，隨著裂縫擴展長度不斷增大，Ⅰ型應(yīng)力強度因子不斷減小。不考慮水力劈裂耦合和考慮水力劈裂耦合作用下的裂尖Ⅰ型應(yīng)力強度因子變化規(guī)律比較一致，具有波動性，隨著裂縫擴展長度不斷增加，無水力劈裂作用下的裂尖Ⅰ型應(yīng)力強度因子要小于水力劈裂作用下的裂尖Ⅰ型應(yīng)力強度因子，這說明裂縫水力劈裂導(dǎo)致裂尖Ⅰ型應(yīng)力強度因子增大，降低了重力壩裂縫的穩(wěn)定性。不考慮裂縫水力劈裂時的Ⅱ型應(yīng)力強度因子比較平緩，隨著裂縫擴展長度不斷增大，Ⅱ型應(yīng)力強度因子變化不大。不考慮水力劈裂耦合和考慮水力劈裂耦合作用下的裂尖Ⅱ型應(yīng)力強度因子變化規(guī)律比較一致，具有波動性。

圖10 有無水力劈裂作用下裂縫長度與Ⅰ、Ⅱ型應(yīng)力強度因子之間的關(guān)系曲線Fig.10 Relation curve between crack length and types Ⅰ, Ⅱ stress intensity factors with or without hydraulic fracturing

6 結(jié) 語

本文考慮裂縫內(nèi)的水壓力作用，建立擴展有限元法框架下重力壩壩踵裂縫水力劈裂耦合數(shù)值模型，采用擴展有限元法對重力壩壩踵裂縫水力劈裂耦合進行了數(shù)值模擬分析，得到以下結(jié)論：

（1）自重、上游庫水壓力和水力劈裂作用下壩踵初始裂縫逐漸向壩基底部擴展，且裂縫擴展朝下游方向，裂縫尖端區(qū)域和壩踵附近都出現(xiàn)應(yīng)力集中現(xiàn)象，且裂縫尖端區(qū)域的主應(yīng)力值明顯高于壩踵附近的。

（2）壩踵初始裂縫擴展前，裂縫內(nèi)水壓力基本與邊界水壓力相同；當裂縫開始擴展時，裂縫內(nèi)水壓力下降；而后裂縫張開寬度不斷增大，裂縫內(nèi)水壓力又會變成邊界全水頭。

（3）無水力劈裂作用下的裂縫開裂角大于水力劈裂作用下的，無水力劈裂耦合作用下的裂縫開裂角小于水力劈裂耦合作用下的。

（4）裂縫水力劈裂會導(dǎo)致裂尖Ⅰ型應(yīng)力強度因子增大，降低重力壩裂縫的穩(wěn)定性，影響大壩安全。