考慮流固耦合的水位降落期壩坡穩(wěn)定分析

田 江 沈長松 王棟良

（1.河海大學 水利水電學院，南京 210098；2.中國水電顧問集團 華東勘測設計研究院，杭州 310014）

我國已建成的各類水庫有8萬多座，但絕大部分興建于20世紀50～70年代，限于當時的技術和經(jīng)濟條件，水庫設計標準偏低、施工質(zhì)量差，經(jīng)過多年運行存在各種安全隱患，急需對其進行除險加固.目前，小（Ⅱ）型病險水庫除險加固工作在全國范圍內(nèi)有條不紊地進行著，水位降落期的土石壩壩坡穩(wěn)定分析是除險加固工程中的重要內(nèi)容之一.因此，選取合理的計算方法來分析水位降落時上游壩坡的穩(wěn)定性具有重要的實際意義.

目前，水位降落期的壩坡穩(wěn)定分析一般采用傳統(tǒng)的極限平衡法［1－2］，該方法主要基于非穩(wěn)定滲流計算原理和極限平衡原理.首先進行水位降落期非穩(wěn)定滲流計算，然后把不同時刻的孔隙水壓力大小應用于極限平衡法，從而計算不同時刻的最小穩(wěn)定系數(shù).該方法除了沒有考慮應力變化對滲流的影響外，還具有極限平衡法自身的缺點，如將土體看作剛體，忽略土體的應力和變形的關系，這給計算結果帶來一定的誤差.

此外，以強度折減法為代表的有限元方法已用于水位降落期壩坡穩(wěn)定分析中［3－4］，該方法基于材料參數(shù)的強度折減原理，并克服了極限平衡法不能考慮土體自身變形對壩坡穩(wěn)定的影響等缺點，但其臨界破壞的判斷準則還有一定的爭議.除了以強度折減法為代表的有限元法外，還有以有限元計算的應力分布結果為基礎，通過張量變換計算出假定滑動面上的穩(wěn)定系數(shù)［5－6］，但這種方法在同時考慮非飽和土流固耦合的影響研究則較少.

水位快速降落過程中，對于滲透系數(shù)較小的均質(zhì)壩，作用于其上游壩坡的水壓力快速消失，而壩體中的超靜孔隙水壓力消散較慢，降落前的飽和區(qū)域變成非飽和區(qū)域，這本質(zhì)上是一個非飽和土的固結問題.本文以某病險水庫為例，以非飽和土流固耦合理論和有限元極限平衡法原理為基礎，對水位降落期上游壩坡進行穩(wěn)定性分析.首先根據(jù)有限元模型進行非飽和土滲流－應力耦合計算，再根據(jù)計算所得各時刻的應力與孔隙水壓力結果進行邊坡穩(wěn)定分析，并與傳統(tǒng)的極限平衡方法進行了對比，得出一些有價值的結論.

1 基本理論

1.1 非飽和土流固耦合理論［7］

1）滲流方程

根據(jù)Darcy定律，非飽和土的二維滲流方程可表示為

式中，kx，ky分別表示x和y方向的滲透系數(shù)；uw為孔隙水壓力；γw為單位水的重量；θw為體積含水率；t為時間.

2）土體本構方程

假設氣壓在整個非飽和區(qū)保持恒定，等于大氣壓力，則非飽和土應力－應變增量形式可以寫成以下形式：

3）有限元平衡方程

應用虛功原理，并分別對公式（1）和公式（2）進行離散，最后耦合方程可以寫成下列形式：

1.2 有限元極限平衡法原理

有限元極限平衡法是以有限元法計算的應力為基礎，通過張量變換方法把高斯點上的應力轉化成假定滑動面上各微段的正應力和剪應力，求解滑動面的穩(wěn)定系數(shù)，運用極限平衡法中滑動面搜索方法可以求出最小穩(wěn)定系數(shù).該方法避免了傳統(tǒng)極限平衡法中的許多假設，計算出的應力更加貼合實際，此時壩坡的穩(wěn)定系數(shù)可按如下公式計算［8］：

式中，∑Sr為沿整個滑面總的最大抗剪力，∑Sm為沿整個滑面總的下滑力.

本文采用基于 Mohr－Coulomb屈服準則的理想彈塑性土體本構模型，以非飽和土流固耦合的應力結果為基礎，求得不同時刻的最小穩(wěn)定系數(shù).在分析中同時考慮了基質(zhì)吸力對壩坡穩(wěn)定的影響，故微小滑弧段的抗剪力和下滑力可表達為：

式中，β為條塊底面長度；c′為有效粘聚力；φ′為有效內(nèi)摩擦角；ua為孔隙氣壓力；uw為孔隙水壓力；θw為體積含水率；θs為飽和體積含水率；θr為殘余體積含水率；σn為條塊底面中心正應力；τm為條塊底面中心剪應力.

2 工程實例及分析

2.1 計算模型及計算參數(shù)

某小（Ⅱ）型水庫，大壩為均質(zhì)土壩，壩頂高程為84.20m，壩頂寬度為3.0m，上下游壩坡均為1∶2.5，上游正常蓄水位80.6m，死水位為72.6m，設計洪水位81.76m，校核洪水位82.96m，下游79m高程處設有6m寬平臺.壩體料為粉質(zhì)粘土，壩基為強風化千枚巖.大壩無排水措施，下游出溢點較高，滲漏現(xiàn)象嚴重，故需對其進行壩坡穩(wěn)定分析，本文只對其進行水位降落期上游壩坡的穩(wěn)定分析.

選取大壩最大斷面作為計算典型斷面，大壩分區(qū)和網(wǎng)格剖分如圖1所示.

圖1 典型剖面及網(wǎng)格剖分圖

耦合計算時需要給定位移邊界條件和水力邊界條件，位移邊界條件為地基左右兩側水平向約束，底部固向約束，壩坡面為自由面；水力邊界條件為上下游均為已知水頭.為了解上游水位降落速度對上游壩坡穩(wěn)定的影響，分析了水庫水位分別以不同的速度從校核洪水位82.96m降落至死水位72.60m后水位保持不變不同工況下的壩坡穩(wěn)定性.工況1、工況2和工況3依次為1.0m／d、2.0m／d和3.0m／d的速度從校核洪水位降落至死水位.

大壩各個分區(qū)的材料參數(shù)根據(jù)地勘資料并參考類似工程選取，具體見表1.

表1 壩體物理力學參數(shù)

壩體填筑料為粉質(zhì)粘土，本次計算采用的粉質(zhì)粘土的土－水特征曲線是根據(jù)其飽和體積含水率和殘余體積含水率，由GeoStudio軟件中內(nèi)置的粉質(zhì)粘土的土－水特征曲線實驗數(shù)據(jù)擬合得到，如圖2所示；基質(zhì)吸力與滲透系數(shù)的關系曲線則根據(jù)Van Genuchten模型估算得出，如圖3所示.

圖2 土－水特征曲線

圖3 基質(zhì)吸力與滲透系數(shù)的關系曲線

2.2 計算結果及其分析

1）上游水位降落速度對壩坡穩(wěn)定系數(shù)的影響

圖4為水位降落至死水位時刻3種工況下壩體內(nèi)浸潤線的位置，壩體內(nèi)孔隙水壓力的消散速度明顯低于上游水位的下降速度，這就形成了“逆流”現(xiàn)象，不利于上游壩坡的穩(wěn)定.圖5反映了不同工況下上游壩坡最小穩(wěn)定系數(shù)隨時間的變化規(guī)律，都呈現(xiàn)先下降后慢慢上升的趨勢.上游水位降落越快，壩坡最小穩(wěn)定系數(shù)降落也越快，最小穩(wěn)定系數(shù)隨著降落速度的增大而減小，主要是因為水位降落速度越快，降落至死水位時間越短，壩體內(nèi)的孔隙水消散時間相對縮短，“逆流”現(xiàn)象越嚴重.

2）基質(zhì)吸力對壩坡穩(wěn)定的影響

基質(zhì)吸力對壩坡穩(wěn)定性的影響程度與壩體中非飽和區(qū)域的大小有關，從表2和圖4分析可知，由于初始狀態(tài)壩體內(nèi)的飽和區(qū)域很大，當水位降落過快，壩體滲透系數(shù)較小時，壩體內(nèi)的超靜孔隙消散很慢，水位降落至死水位時壩體內(nèi)形成的相對非飽和區(qū)域很小，此時，降落期基質(zhì)吸力對壩坡穩(wěn)定的影響相對較小.因此，是否考慮基質(zhì)吸力對壩坡穩(wěn)定的影響可根據(jù)實際工程在水位降落工程中壩體內(nèi)非飽和區(qū)域的大小而定.

表2 最小穩(wěn)定系數(shù)對比

3）耦合與非耦合情況下的計算結果對比

如何精確確定壩體中孔隙水壓力的分布情況是水位降落期壩坡穩(wěn)定分析的一個重要問題.以工況1為例，將忽略外部荷載變化對滲流場影響的非飽和非穩(wěn)定滲流計算所得不同時刻的孔隙水壓力應用于穩(wěn)定分析模塊進行穩(wěn)定分析得出的計算結果作為非耦合工況，比較耦合與非耦合兩種情況下孔隙水壓力分布規(guī)律的差異.

圖6為上述兩種情況下相同滑動面上的孔隙水壓力分布情況（橫坐標對應計算模型的x軸），用非飽和土流固耦合計算出的孔隙水壓力比忽略應力場變化對滲流影響的非飽和非穩(wěn)定滲流計算所得的結果較小，且更加符合實際，這與文獻［9］中得出的結論是相似的.這是因為水位降落過程是一個卸載過程，上游壩體土體孔隙體積在一定程度上有所膨脹，從而在一定程度上減小了孔隙水壓力.

圖6 滑動面孔隙水壓力分布

4）本次有限元極限平衡法與傳統(tǒng)極限平衡法的結果對比

為了說明該方法計算結果的準確性，以工況1并考慮基質(zhì)吸力影響為例，分別采用Morgentern－Price法、Bishop法、Spencer法、Janbu法和Ordinary法對其進行了壩坡穩(wěn)定性分析.極限平衡法計算所采用的不同時刻孔隙水壓力大小同樣由非飽和土耦合計算得到，各種方法計算出的穩(wěn)定系數(shù)見表3.

表3 計算結果對比

有限元極限平衡法計算出的最小穩(wěn)定系數(shù)為1.485，雖與Bishop法計算出的最小穩(wěn)定系數(shù)相差較大，為10.1%，但與理論相對嚴謹?shù)臉O限平衡計算公式Morgentern－Price法、Spencer法計算的穩(wěn)定系數(shù)都比較接近，相對誤差在2.2%之內(nèi)，說明了基于非飽和土流固耦合理論的有限元極限平衡法在計算水位降落下的壩坡穩(wěn)定是合理可行的.

5）材料參數(shù)對壩坡穩(wěn)定性的影響分析

本文采用基于 Mohr－Coulomb屈服準則的理想彈塑性土體本構模型，以工況1為例分析研究理想彈塑性本構模型中主要材料參數(shù)（彈性模量、泊松比、粘聚力和內(nèi)摩擦角）對壩坡穩(wěn)定性的影響.

從圖7分析可知，基于非飽和土固結理論的有限元極限平衡法中，在所取的參數(shù)變化范圍內(nèi)，穩(wěn)定系數(shù)與材料參數(shù)呈線性變化.彈性模量對穩(wěn)定系數(shù)幾乎無影響，穩(wěn)定系數(shù)隨著泊松比的增加而增加，但影響程度較小，而粘聚力和內(nèi)摩擦角對穩(wěn)定系數(shù)的影響較大，最大影響值在7%～15%之間.

圖7 材料參數(shù)與穩(wěn)定系數(shù)的關系曲線

3 結 論

基于非飽和土流固耦合理論和有限元極限平衡法原理，對某水庫在水位降落過程中上游壩坡的穩(wěn)定性進行了分析研究，得出以下結論：

1）采用基于非飽和土流固耦合理論的有限元法分析簡單均質(zhì)壩水位降落過程中壩坡的穩(wěn)定性是合理可行的，其計算得出孔隙水壓力分布和穩(wěn)定穩(wěn)定系數(shù)更符合實際情況，具有較好的推廣前景.因本文僅對材料簡單的均質(zhì)壩進行了分析，對于壩體材料分區(qū)多、地基條件復雜的土石壩，使用該方法進行穩(wěn)定分析時，其結果的合理性有待進一步研究.

2）基于非飽和土固結理論的有限元極限平衡法中，材料的粘聚力和內(nèi)摩擦角對壩坡穩(wěn)定系數(shù)的影響最大，泊松比影響相對較小，彈性模量幾乎無影響，實際工程應根據(jù)實驗獲得準確而合理的計算參數(shù).

3）基質(zhì)吸力對壩坡穩(wěn)定性的影響程度與壩體中非飽和區(qū)域的大小有關，當初始狀態(tài)壩體內(nèi)飽和區(qū)域較大，降落過程中壩體內(nèi)形成的非飽和區(qū)域又相對較小，此時基質(zhì)吸力對壩坡穩(wěn)定的影響相對較小，是否考慮基質(zhì)吸力對壩坡穩(wěn)定的影響可根據(jù)實際工程在水位降落過程中壩體內(nèi)非飽和區(qū)域的大小而定.

4）對于滲透系數(shù)較小的均質(zhì)壩，上游水位快速降落時壩體中的孔隙水往往來不及消散而形成逆流現(xiàn)象，這將大大降低上游壩坡的抗滑穩(wěn)定性甚至會出現(xiàn)滑坡，水位降落越快越不利于穩(wěn)定，所以在大壩運行管理過程中要合理調(diào)度，控制上游水位下降速度.

［1］ 陳祖煜.土質(zhì)邊坡穩(wěn)定分析——原理、方法和程序［M］.北京：中國水利水電出版社，2003.

［2］ 謝定義，姚仰平，黨發(fā)寧.高等土力學［M］.北京：高等教育出版社，2007.

［3］ 賈蒼琴，黃茂松，王貴和.非飽和非穩(wěn)定滲流作用下土坡穩(wěn)定分析的強度折減有限元方［J］.巖石力學與工程學報，2007，26（6）：1290－1296.

［4］ 周貴云，李同春.基于非飽和土固結理論的有限元強度折減法［J］.巖土力學，2008，29（4）：1133－1137.

［5］ 曾亞武，田偉明.邊坡穩(wěn)定性分析的有限元法與極限平衡法的結合［J］.巖石力學與工程學報，2005，24（增2）：5355－5359.

［6］ 王學武，許尚杰，黨發(fā)寧，等.水位驟降時的非飽和壩坡穩(wěn)定分析［J］.巖土力學，2010：31（9）：2760－2764.

［7］ GEO－SLOPE International Ltd.Stress－Deformation Modeling with SIGMA／W 2007［M］.Canada：GEOSLOPE International Ltd，March，2008：183－192，338－339.

［8］ GEO－SLOPE International Ltd.SIGMA／W Example file：Rapid Drawdown［EB／OL］.Http：／／www.geoslope.com，2012－5－16.