考慮滲流—應力耦合的黏土心墻壩邊坡穩定分析

黃 香，崔志瑞，郭瑩瑩，李 穎

(1.三峽金沙江云川水電開發有限公司，云南 昆明 650224；2.華北水利水電大學， 河南 鄭州 450046)

1 概述

對于黏土心墻壩而言，滲流與穩定問題關系到大壩能否安全運行[1]。而邊坡穩定分析是黏土心墻壩設計的基礎，通過邊坡穩定分析可反映出大壩設計方案的合理性，為大壩結構設計提供科學依據[2]。目前關于黏土心墻壩邊坡穩定的研究成果較多，如文獻[3-4]等均采用傳統的邊坡穩定計算方法，對一些黏土心墻壩的邊坡穩定情況進行了計算分析。許多研究成果[5-6]表明，在天然狀態下，滲流與土體孔隙空間相互關聯，滲流場與應力場是一種耦合關系。對土石壩(包括黏土心墻壩)進行邊坡穩定分析時，只考慮單獨的滲流場或應力場進行計算，得到的結果會與實際情況存在差異，采用傳統的邊坡穩定計算方法，會存在一些不足。另一些文獻[7-8]考慮了滲流—應力雙場耦合的作用，并對黏土心墻壩的應力變形分布規律和工作性態進行了計算分析，但對雙場耦合作用對邊坡穩定的影響并沒有進行更深入的研究。綜合分析文獻可知，目前考慮滲流—應力雙場耦合的土石壩邊坡穩定分析，較多研究偏重于尾礦壩或均質土壩[9-10]，而考慮雙場耦合的黏土心墻壩邊坡穩定的研究則比較少。

另外，從研究方法層面，目前大多數研究采用有限元強度折減法進行土石壩邊坡穩定計算分析。通過此法進行邊坡穩定計算時，邊坡失穩判據大致有三種：判據一，數據計算是否收斂；判據二，塑性區是否貫通；判據三，特征點位移突變[11]。根據判據進行土石壩邊坡穩定計算分析時，大多數研究采用Drucker-Prager模型等彈塑性模型進行壩體的位移和應力計算[9-11]。盡管上述理想彈塑性模型能很好描述巖土體的非線性特性，但此類模型如屈服面存在尖角導致計算繁瑣，且收斂速度緩慢，甚至不收斂的缺陷[9]。因此，當采用此類模型進行邊坡穩定計算分析，并結合判據一和判據二進行邊坡失穩破壞判斷時，容易造成誤判。而鄧肯—張模型做為一種增量型彈性模型，可克服上述彈塑性模型計算時不易收斂的缺陷，而且可以在一定程度上反映土體的彈塑性變形，目前鄧肯—張模型在巖土體的應力變形計算分析中得到了廣泛的應用，但考慮滲流—應力雙場相互耦合，相互影響，并基于鄧肯—張模型的邊坡穩定有限元強度折減法的研究則相對尚較少。

綜上所述，本文利用ANSYS軟件的APDL命令流語言并結合Fortran語言，編制模塊化的滲流—應力耦合計算程序，并基于鄧肯—張模型的邊坡穩定有限元強度折減法，對某黏土心墻壩的邊坡穩定進行計算分析，以期為黏土心墻壩的邊坡穩定分析提供借鑒和參考。

2 計算原理

2.1 滲流—應力耦合計算原理

目前，滲流—應力雙場耦合計算方法主要為直接耦合法和間接耦合法，考慮直接耦合計算時需建立復雜的數學模型且求解困難，因此，本文采用間接耦合方法進行雙場的耦合計算。

1)滲流場對應力場的影響

滲流場對應力場的影響有動水壓力(滲透體積力)和靜水壓力(滲透壓力)兩種形式。靜水壓力(滲透壓力)為作用與接觸面上的力，計算時可按面力作用形式計入節點荷載，即：

(1)

式中NT(x)為形函數；Γ為滲流邊界；p為滲透壓力，其可表示為：

p=γw(H-y)

(2)

式中γw為水的容重；H為水頭；y為位置y向坐標。

動水壓力(滲透體積力)可視為與自重類似，計算時按體力計入節點荷載，即：

fs=?ΩNT(x)fpdΩ

(3)

式中NT(x)為形函數；Ω為滲流計算區域；fp為動水壓力，其可表示為：

(4)

考慮滲流場影響的應力場控制方程，則可表示為：

(5)

2)應力場對滲流場的影響

在荷載作用下，應力場的變化，會導致巖土體的體積及孔隙發生變化，從而導致其滲透系數發生改變。應力與滲透系數的關系可用文獻[6]推薦的關系式表達，即：

K=K0exp[-β(σ-p)]

(6)

式中σ為土體應力；K0為初始滲透系數；β為經驗系數，取為3e-7；p為靜水壓力。

考慮應力場影響的滲流場控制方程，可表示為：

(7)

式中H為水頭分布；K(σij)為滲透系數，可由式(6)表示；Ω為滲流計算區域；Γ1為已知流量邊界；Γ2為自由面邊界；Γ3為已知水頭邊界；n2和n3為對應邊界的法線方向。

2.2 基于鄧肯—張模型的強度折減法原理

利用有限元強度折減法時，黏聚力c、內摩擦角φ按下式進行計算：

cf=c/Fs,φf=arctan(tanφ/Fs)

(8)

式中Fs為折減系數(即安全系數)；cf、φf分別為折減后的材料內聚力和內摩擦角。

對于鄧肯—張模型，其切線模量Et可由下式表示：

(9)

式中Ei為初始彈性模量；Rf為破壞比；兩者的具體表達式可參見文獻[5]。而上式中(σ1-σ3)f表示破壞時的主應力差，可由Mohr-Coulomb定律推導得到，具體為：

(10)

將式(10)及Ei的表達式代入式(9)，可得切線模量的表達式為：

(11)

從式(11)可以看出，只要對式中的黏聚力c和內摩擦角φ進行折減，即可實現與Drucker-Prager模型或Mohr-Coulomb模型相似的材料強度折減效果，另外本文在進行強度折減時，還考慮將初始切線模量k進行折減，以充分滿足土體的變形性質，即：

kf=k/Fs

(12)

將式(8)和式(12)代入式(11)，即可得到考慮有限元強度折減法的鄧肯—張模型切線模量計算表達式，如式(13)所示：

(13)

3 工程實例及結果分析

3.1 工程實例

某黏土心墻壩，工程總庫容為110.87萬m3，工程規模為小(1)型，工程等別為Ⅳ等，永久性水工建筑物級別為4級，臨時性水工建筑物級別為5級。跨河布置總長為185.00 m，壩頂高程為1 804.30 m，壩頂設高1.2m混凝土防浪墻，墻頂高程為1 805.50 m，壩頂寬度為5.00 m。大壩位于河谷處，最大壩高斷面及材料分區如圖1所示。限于篇幅，本文中只給出正常蓄水工況的計算結果，大壩正常蓄水位為1 800.08 m。根據實驗資料[12]，本次計算分析時，壩體材料采用鄧肯—張E-v模型，具體材料參數見表1，壩基及帷幕各材料假定為線彈性材料(見表2)。

圖1 黏土心墻壩計算斷面及材料分區示意(單位：m)

表1 壩體鄧肯—張E-v模型材料參數

表2 壩基彈性模型材料參數

3.2 計算結果分析

進行滲流—應力耦合計算時，計算斷面的有限元網格如圖2所示。有限元模型的頂邊界截至壩頂高程為1 840.30 m，底邊界截取至壩基帷幕最深處以下，至高程為1 680.00 m處。本次有限元模擬計算，滲流計算邊界條件取為：壩基底部為不透水邊界，左側為上游水頭邊界，右側無水，因此不加水頭，但下游側斜邊視為可能溢出邊界；應力計算部分邊界條件為：壩基底部為x、y向位移約束，壩基上下游側x向位移約束。模型的坐標系選用笛卡爾直角坐標系，x軸指向下游為正，y軸向上為正。

圖2 計算斷面有限元網格示意

1)滲流計算結果分析

從圖3可以看出:考慮和不考慮耦合作用時，壩體的水頭等值線分布規律大致相同。但從中也可看出考慮耦合作用時，水頭為90.737 m和92.884 m的等值線更靠近下游，這是由于考慮耦合作用后，受應力場的影響使材料的滲透系數有所增大，心墻削減水頭的能力減弱。

a 考慮耦合作用

b 不考慮耦合作用

2)應力變形計算結果分析

從圖4可以看出，考慮耦合后壩體的y向位移與不考慮耦合時相比變大，這是由于耦合作用滲透水壓力的豎向分力使壩體的土體有下沉趨勢，導致考慮耦合效應后，壩體的y向位移增大。

從圖5可以看出，考慮耦合和不考慮耦合情況下，壩體的y向應力均符合一般規律。但不考慮耦合情況下，心墻y向應力小于兩側壩殼料的應力，具有拱效應，這是因為心墻的變形模量較兩側壩殼料小，心墻位移大，兩側壩殼料位移小，心墻部分自重傳遞到壩殼，反映了心墻內部應力小于兩側壩殼料的應力。考慮耦合情況下，心墻部位也產生了明顯的拱效應，但考慮耦合情況下，心墻上下游兩側應力分布不成對稱分布，且心墻上下游兩側同一高程處，上游側y向應力大于下游側，這是由于耦合情況下，考慮了滲透水壓力的影響，上游水頭較高，滲透水壓力較大，導致上游土體壓力大于下游土體的壓力。

a 考慮耦合作用

b 不考慮耦合作用

a 考慮耦合作用

b 不考慮耦合作用

3)邊坡穩定計算結果分析

考慮邊坡失穩判據1和判據2的特點及局限性以及失穩判據與鄧肯—張模型的相容性，本文采用判據3對邊坡的失穩破壞進行判定。采用判據3時，需確定特征點后將水平或豎向位移作為位移突變的判據。因為黏土心墻變形比兩側壩殼料大，更易導致邊坡的失穩破壞，故本文選取心墻頂的9 672號節點作為特征點(圖2中A點)。在開始計算時，折減系數以0.1為間隔進行變化，以找出安全系數的大致范圍，再以0.025或更小的間隔進行折減，直到找出更精確的安全系數。根據上述方法得到的特征點豎向位移與安全系數的關系曲線如圖6所示。

圖6 特征點豎向位移與安全系數的關系曲線示意

從圖6可以看出，當考慮耦合時，折減系數由1.2增大為1.3時，特征點的位移有明顯的突變。對折減系數進行加密計算，當安全系數在1.250～1.275之間時，有明顯的的位移突變，因此，考慮耦合計算時，本文所研究黏土心墻壩的邊坡穩定安全系數為1.250。同樣可以看出，當不考慮耦合計算時，當折減系數由1.3增大為1.4時，特征點的位移有明顯的突變。同樣，對折減系數進行加密，可以看出安全系數在1.300～1.400之間時，有明顯的的位移突變，因此，不考慮耦合計算時，可確定本文所研究黏土心墻壩的邊坡穩定安全系數為1.375。從圖6還可以看出，考慮耦合后計算所得壩體邊坡的穩定安全系數比不考慮滲流—應力耦合時小，究其原因是，一方面考慮耦合作用時，滲透水壓力豎向分力使壩體的土體有下沉趨勢，滲流力加大了壩坡的滑移力，從而增大了壩體邊坡土體的下滑趨勢，另一方面考慮耦合作用時，受滲流場影響，土體的變形模量變小，土體軟化，更易產生變形，對壩坡抗滑穩定不利。

4 結語

本文利用ANSYS軟件，建立了某黏土心墻壩的有限元模型，計算結果表明：

1)考慮耦合作用后，受應力的影響，壩體土料的滲透系數增大后會使心墻削減水頭的能力有所減弱。

2)考慮耦合作用后，滲透水壓力豎向分力使壩體的土體有明顯下沉趨勢，豎向位移較不考慮耦合時明顯增大，而且滲透水壓力的作用使得壩體上游側土體的應力明顯比下游大。

3)考慮耦合作用進行邊坡穩定計算分析時，滲流場作用不僅會影響壩坡的滑移力，而且會影響土體的變形模量，對壩坡的抗滑穩定不利，從而使計算得到邊坡穩定安全系數明顯比不考慮耦合作用的小。為偏安全，在進行黏土心墻壩有限元計算分析及邊坡穩定計算時，建議考慮滲流場—應力場的雙場耦合效應。