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

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

1 概述

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

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

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

2 計算原理

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

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

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

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

(1)

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

p=γw(H-y)

(2)

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

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

fs=?ΩNT(x)fpdΩ

(3)

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

(4)

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

(5)

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

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

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

(6)

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

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

(7)

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

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

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

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

(8)

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

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

(9)

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

(10)

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

(11)

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

kf=k/Fs

(12)

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

(13)

3 工程實例及結果分析

3.1 工程實例

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

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

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

表2 壩基彈性模型材料參數(shù)

3.2 計算結果分析

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

圖2 計算斷面有限元網(wǎng)格示意

1)滲流計算結果分析

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

a 考慮耦合作用

b 不考慮耦合作用

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

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

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

a 考慮耦合作用

b 不考慮耦合作用

a 考慮耦合作用

b 不考慮耦合作用

3)邊坡穩(wěn)定計算結果分析

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

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

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

4 結語

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

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

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

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