聶相田, 莊濮瑞, 焦延濤, 王博

(1.華北水利水電大學(xué) 水利學(xué)院,河南 鄭州 450046; 2.水資源高效利用與保障工程河南省協(xié)同創(chuàng)新中心,河南 鄭州 450046; 3.河南省水環(huán)境模擬與治理重點實驗室,河南 鄭州 450046)

滲流問題是影響土壩安全的重要因素之一。只有深入分析土壩的滲流分布,才能確保土壩運行安全[1]。而土壩的滲流場與應(yīng)力場相互作用,是一種耦合關(guān)系[2]。在分析土壩邊坡的穩(wěn)定性時,若只考慮單獨的滲流場或應(yīng)力場進(jìn)行計算,得到的結(jié)果必然與實際情況存在差異。近年來,已有許多學(xué)者在土壩邊坡或者其他巖質(zhì)邊坡的穩(wěn)定性分析中,考慮滲流場與應(yīng)力場的耦合作用,并取得了一些研究成果。如苗麗等[3]基于多孔介質(zhì)的滲流特性和土的非線性本構(gòu)關(guān)系,研究了滲流場與應(yīng)力場的耦合作用,并對某非均質(zhì)土壩的邊坡穩(wěn)定情況進(jìn)行了計算分析。李宗坤等[4]探討了滲流場與應(yīng)力場相互作用的機理,提出了考慮滲流與應(yīng)力直接耦合作用的邊坡穩(wěn)定分析方法。賈善坡等[5]將損傷力學(xué)理論引入到滲流-應(yīng)力耦合計算分析中,對比利時核廢料庫開挖過程中圍巖的損傷演化以及滲流場和應(yīng)力場耦合過程進(jìn)行了計算分析,得到了圍巖損傷特性、孔隙壓力以及滲透性的變化規(guī)律。柳厚祥等[6]考慮應(yīng)力場與滲流場的耦合特性,對某尾礦壩的非穩(wěn)定滲流進(jìn)行了計算分析。宋傳旺等[7]采用有限元軟件MIDAS建立了某尾礦壩的計算模型,探究了應(yīng)力場和滲流場耦合作用對該尾礦壩穩(wěn)定性的影響。沈振中等[8]應(yīng)用迭代解耦技術(shù)建立了基于無單元法的應(yīng)力場與滲流場耦合分析計算模型,并研制開發(fā)了計算程序。

綜上所述,目前關(guān)于巖土體邊坡穩(wěn)定性分析的研究成果中,考慮滲流-應(yīng)力耦合影響的已有不少,而且多數(shù)是結(jié)合有限元強度折減法進(jìn)行計算的。采用有限元強度折減法進(jìn)行邊坡穩(wěn)定性計算時,邊坡失穩(wěn)判據(jù)大致有3種：判據(jù)一,數(shù)據(jù)計算是否收斂;判據(jù)二,塑性區(qū)是否貫通;判據(jù)三,特征點位移突變[9]。根據(jù)上述判據(jù),目前大多采用Drucker-Prager模型或Mohr-Coulomb模型進(jìn)行巖土體的位移和應(yīng)力計算,如文獻(xiàn)[4][9][10]等。盡管采用Drucker-Prager模型或Mohr-Coulomb模型等理想彈塑性模型能很好地描述巖土體的非線性特性,但此類模型存在一些缺陷,如屈服面存在尖角，導(dǎo)致計算煩瑣,且收斂速度緩慢,甚至不收斂[4]。因此,當(dāng)采用此類模型進(jìn)行邊坡穩(wěn)定性分析,并結(jié)合判據(jù)一和判據(jù)二進(jìn)行邊坡失穩(wěn)破壞判斷時,容易誤判,甚至結(jié)果失真。而鄧肯-張(Duncan-Chang)模型是一種基于三軸試驗數(shù)據(jù)得出的增量彈性模型,不僅適用于描述土體的非線性變形特征,而且可以在一定程度上反映土體的彈塑性變形,同時,該模型所需參數(shù)不多,各參數(shù)的物理意義明確,獲得參數(shù)的途徑亦較簡單,計算也很容易收斂。因此,鄧肯-張模型在巖土體的應(yīng)力變形計算分析中得到了廣泛的應(yīng)用?；卩嚳?張模型的邊坡穩(wěn)定有限元強度折減法也有學(xué)者進(jìn)行過研究[11-12],但考慮滲流-應(yīng)力相互耦合、相互影響,并基于鄧肯-張模型的邊坡穩(wěn)定有限元強度折減法的研究則相對較少。因此,本文基于ANSYS軟件,利用其APDL命令流語言并結(jié)合Fortran語言,編制模塊化的滲流-應(yīng)力耦合計算程序,并基于鄧肯-張模型的邊坡穩(wěn)定有限元強度折減法對某均質(zhì)土壩的邊坡穩(wěn)定性進(jìn)行計算分析,以期為土石壩的邊坡穩(wěn)定分析提供一些參考。

1 計算原理

1.1 滲流計算原理

ANSYS軟件本身并沒有集成滲流計算模塊,但集成了溫度場分析模塊,鑒于溫度場的計算原理與滲流場的相似,故可以借用其溫度場分析模塊來進(jìn)行滲流場的計算分析。關(guān)于滲流場與溫度場計算原理的相似性,文獻(xiàn)[1]和[13]均從基本理論、微分方程、初始邊界條件3個方面進(jìn)行了詳細(xì)的論證,本文不再贅述。

1.2 應(yīng)力變形計算原理

目前,描述巖土體應(yīng)力-應(yīng)變關(guān)系的數(shù)學(xué)模型很多,主要包括彈性模型和彈塑性模型兩大類,但真正在工程實踐中得到廣泛應(yīng)用的卻為數(shù)不多。其中,鄧肯-張模型因其具有多種優(yōu)點而獲得了最為廣泛的應(yīng)用,它包括初始的E-v模型和修正后的E-B模型。目前ANSYS軟件的各版本中尚未嵌入鄧肯-張模型,但由于其計算公式簡單,可以通過APDL語言進(jìn)行編程實現(xiàn)。本文即通過APDL語言創(chuàng)建鄧肯-張E-v模型的宏文件,然后通過調(diào)用該宏文件來實現(xiàn)土體的應(yīng)力計算。關(guān)于鄧肯-張模型的詳細(xì)介紹,可參考文獻(xiàn)[1]和[13]。

1.3 滲流-應(yīng)力耦合計算原理

在滲流過程中,滲流給土體介質(zhì)施加滲透作用力,應(yīng)力或應(yīng)變引起其滲透性質(zhì)變化,繼而影響滲流場的分布,從而反過來影響巖土體的應(yīng)力場[14]。目前,進(jìn)行滲流-應(yīng)力耦合計算大致有兩種途徑：直接耦合和間接耦合。直接耦合計算需求出包含有滲流場與應(yīng)力場相關(guān)未知量的復(fù)雜數(shù)學(xué)模型的解析解,這是很難實現(xiàn)的。因此,本文采用間接耦合方法,通過迭代計算來求得滲流-應(yīng)力耦合作用下的滲流場和應(yīng)力場。

1.3.1 滲流場對應(yīng)力場的影響

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

(1)

p=γw(H-y)。

(2)

式中：fw為節(jié)點所受面力;NT(x)為形函數(shù);Γ為滲流邊界;p為節(jié)點靜水壓力列陣,γw為水的容重;H為節(jié)點水頭列陣;y為節(jié)點y向坐標(biāo)列陣。

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

fs=?ΩNT(x)fpdΩ，

(3)

(4)

式中：fs為節(jié)點所受體力;Ω為滲流計算區(qū)域;fp為動水壓力。

1.3.2 應(yīng)力場對滲流場的影響

應(yīng)力場的變化會引起巖土體的體積及孔隙發(fā)生變化,從而導(dǎo)致其滲透系數(shù)發(fā)生改變。應(yīng)力與滲透系數(shù)的關(guān)系用文獻(xiàn)[14]推薦的關(guān)系式來表達(dá),即：

K=K0exp[-β(σ-pe)]。

(5)

式中：K為土體單元的滲透系數(shù)列陣;σ為土體單元的應(yīng)力列陣;pe為土體單元的靜水壓力列陣,本文取每個單元4個節(jié)點靜水壓力的平均值作為土體單元的靜水壓力值;K0為土體的初始滲透系數(shù);β為經(jīng)驗系數(shù),取為3e-7。

1.3.3 耦合分析求解步驟

步驟1選取一個初始滲透系數(shù),利用有限元方法求解滲流場,得區(qū)域內(nèi)的水頭分布H(x,y,z)。

步驟2將滲流場分析得到的水頭H(x,y,z)分別代入式(2)與式(4)中,求解各節(jié)點的靜水壓力p和動水壓力fp。

步驟3將各節(jié)點的靜水壓力p和動水壓力fp分別代入式(1)和式(3),根據(jù)有限元法求解應(yīng)力場,獲取應(yīng)力分布σ(x,y,z)。

步驟4根據(jù)各單元4個節(jié)點的靜水壓力值求得各單元的靜水壓力值pe,將σ(x,y,z)和pe代入式(5),求解新的滲透系數(shù),并替代初始的滲透系數(shù),重新進(jìn)行滲流場的求解。

步驟5重復(fù)步驟1—4，進(jìn)行迭代計算,直至滿足如下精度要求：

|Hn+1(x,y,z)-Hn(x,y,z)|≤εH。

(6)

式中：Hn+1(x,y,z)和Hn(x,y,z)分別為第n+1次和第n次迭代求得的滲流場水頭分布;εH為水頭的求解精度。

1.4 基于鄧肯-張模型的強度折減法

有限元強度折減法得到的安全系數(shù)是基于強度儲備概念的安全系數(shù),黏聚力c、內(nèi)摩擦角φ按下式進(jìn)行折減[10]：

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

(7)

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

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

(8)

(9)

式中：Ei為初始彈性模量;Rf為破壞比;(σ1-σ3)f為破壞時的主應(yīng)力差,可由Mohr-Coulomb定律推導(dǎo)得到。

將式(9)及Ei的表達(dá)式代入式(8),可得切線模量的最終表達(dá)式為：

(10)

式中：k為初始切線模量;n為模量指數(shù);Pa為大氣壓強,取101.4 kPa。

從式(10)可以看出,只要對式中的黏聚力c和內(nèi)摩擦角φ進(jìn)行折減,即可實現(xiàn)與Drucker-Prager模型或Mohr-Coulomb模型相似的材料強度折減效果,從而保證有限元強度折減法的順利實現(xiàn),也就是說,基于鄧肯-張模型的有限元強度折減法在原理上是可行的。

2 工程實例分析

2.1 工程實例

王莊水庫大壩為均質(zhì)連續(xù)的各向同性土壩,壩高22.50 m,壩頂高程613.95 m,壩底高程591.00 m。壩頂總長128.00 m,壩頂寬4.00 m,大壩上游坡比為1∶2.5,下游坡比為1∶2,下游設(shè)貼坡排水。校核洪水位612.50 m,設(shè)計洪水位611.40 m,正常蓄水位609.50 m,3種工況下下游均無水。限于篇幅,本文中只給出運行期正常蓄水位工況的計算結(jié)果,所謂運行期即指水庫蓄水后壩體和壩基形成穩(wěn)定的滲流場。因此,本次計算假定土體為飽和土體,土料的飽和滲透系數(shù)取為2.23×10-8m/s。對于穩(wěn)定滲流計算而言,浸潤面以上節(jié)點對于流場計算貢獻(xiàn)不大,可稱之為虛點,本文參考相關(guān)文獻(xiàn)[1]的經(jīng)驗,采用ANSYS提供的死活單元技術(shù),將處于浸潤線上部的單元“殺死”,使其不參與計算,而只“激活”浸潤線下部的單元網(wǎng),當(dāng)進(jìn)行應(yīng)力場分析時重新“激活”浸潤線上部單元,使其參與計算。參考相關(guān)文獻(xiàn)[3,15]及試驗資料,本文采用鄧肯-張E-v模型進(jìn)行應(yīng)力計算,所用的飽和土體材料參數(shù)見表1。

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

2.2 滲流-應(yīng)力耦合計算結(jié)果分析

進(jìn)行滲流-應(yīng)力耦合計算時,為節(jié)省計算時間,本文只建立壩體的有限元模型。滲流計算邊界條件取為：壩體底部視為不透水邊界；左側(cè)為上游水頭邊界；右側(cè)無水,因此不加水頭,但右側(cè)斜邊視為可能溢出面。應(yīng)力計算部分邊界條件為：壩體底部為x、y向位移約束。模型的坐標(biāo)系選用笛卡爾直角坐標(biāo)系,x軸指向下游為正,y軸向上為正。

在正常蓄水位工況下,滲流-應(yīng)力耦合計算時壩體y向位移云圖與非耦合計算時壩體y向位移云圖的對比如圖1所示。從圖1中可以看出,考慮滲流-應(yīng)力耦合作用后,壩體的y向位移明顯增大。這是因為考慮滲流-應(yīng)力耦合作用時,一方面,滲透水壓力的豎向分力使壩體的土體有下沉趨勢;另一方面,地下水動水壓力作用于土體骨架,改變了原有土體顆粒結(jié)構(gòu)的排列,從而導(dǎo)致壩體的y向位移明顯增大。

圖1 正常蓄水位工況下壩體y向位移

2.3 邊坡穩(wěn)定計算結(jié)果分析

根據(jù)前文所述的基于鄧肯-張模型的有限元強度折減法原理,進(jìn)行壩體邊坡的穩(wěn)定計算?？紤]到邊坡失穩(wěn)判據(jù)一和判據(jù)二的特點及局限性,并考慮到失穩(wěn)判據(jù)與鄧肯-張模型的相容性,本文采用判據(jù)三(即特征點位移突變判據(jù))來進(jìn)行邊坡的失穩(wěn)破壞判定。采用判據(jù)三時,需確定特征點,并將其水平或豎向位移作為位移突變的判據(jù)。為此,本文選取壩頂?shù)?2號節(jié)點作為特征點,用以監(jiān)控壩體位移的改變。72號節(jié)點具體位置如圖2中A點所示。

采用強度折減法計算時,首先可使折減系數(shù)以0.1為間隔進(jìn)行變化,大概確定安全系數(shù)范圍;然后再以0.01或更小的間隔折減土體的黏聚力和內(nèi)摩擦角,直到找出更精確的安全系數(shù)。

根據(jù)上述方法得到的特征點豎向位移與安全系數(shù)的關(guān)系曲線如圖3所示。從圖3可以看出：

1)考慮滲流-應(yīng)力耦合計算時,若以0.1的間隔改變折減系數(shù),當(dāng)折減系數(shù)介于1.2和1.3之間時，特征點的豎向位移有明顯的突變。為了找出更準(zhǔn)確的安全系數(shù),在1.2與1.3之間以0.025的間隔加密折減系數(shù),安全系數(shù)在1.225與1.250之間時有明顯的位移突變。因此,當(dāng)考慮滲流-應(yīng)力耦合作用時,可確定本文所研究的均質(zhì)土壩的邊坡穩(wěn)定安全系數(shù)為1.225。

圖3 特征點豎向位移與安全系數(shù)的關(guān)系曲線

2)不考慮滲流-應(yīng)力耦合計算時,當(dāng)折減系數(shù)介于1.3和1.4之間時特征點的豎向位移有明顯的突變。同樣,對折減系數(shù)以0.025的間隔進(jìn)行加密,安全系數(shù)在1.325與1.350之間時特征點有明顯的豎向位移突變。因此,當(dāng)不考慮滲流-應(yīng)力耦合作用時,可確定本文所研究的均質(zhì)土壩的邊坡穩(wěn)定安全系數(shù)為1.325。

3)考慮滲流-應(yīng)力耦合作用后,壩體邊坡的穩(wěn)定安全系數(shù)較不考慮滲流-應(yīng)力耦合作用時有所減小。原因應(yīng)為考慮耦合作用時,一方面,滲流場對土體產(chǎn)生的滲流力在計算時直接加入到壩體原有的應(yīng)力場里面進(jìn)行計算,滲流力加大了壩坡的滑移力,從而降低了大壩的抗滑穩(wěn)定性;另一方面,受滲流場的影響,土體出現(xiàn)軟化,更容易發(fā)生變形,對壩坡抗滑穩(wěn)定不利。因此,這也說明了考慮滲流-應(yīng)力耦合作用時,滲流場對壩體邊坡的穩(wěn)定是不利的。

3 結(jié)論

基于滲流場與應(yīng)力場的相互作用原理,本文利用ANSYS軟件,提出了考慮滲流-應(yīng)力耦合作用下基于鄧肯-張模型的有限元強度折減法的計算方法。計算結(jié)果表明：

1)滲透水壓力的豎向分力使壩體的土體有明顯的下沉趨勢,因此對土石壩進(jìn)行應(yīng)力-變形計算分析時,應(yīng)充分考慮滲透水壓力的影響。

2)進(jìn)行邊坡穩(wěn)定計算分析時,考慮滲流-應(yīng)力耦合作用計算得到的邊坡穩(wěn)定安全系數(shù)明顯比不考慮滲流-應(yīng)力耦合作用的小,滲流場對壩體邊坡的穩(wěn)定是不利的。因此,在進(jìn)行土石壩邊坡穩(wěn)定計算時,應(yīng)盡量考慮滲流-應(yīng)力的耦合作用,使計算結(jié)果更為準(zhǔn)確。