水位降落時(shí)土壩壩體滲流特性的交互式有限元模型

臧育櫻，田 淳，趙燕兵

(太原理工大學(xué) 水利科學(xué)與工程學(xué)院，太原 030024)

1 研究背景

土壩滲流和壩坡穩(wěn)定性是影響土壩安全的關(guān)鍵問題。壩前水位降落造成壩體內(nèi)滲流場隨之發(fā)生改變，不同程度地影響壩內(nèi)孔隙水壓力分布，危害壩體穩(wěn)定性，極易造成上游滑坡等重大災(zāi)害[1-3]。

壩坡穩(wěn)定性分析常用方法有極限平衡法和有限元法。目前對庫水位降落時(shí)壩坡穩(wěn)定性分析多采用極限平衡法[4]，例如張大偉等[5]建立二維模型對庫水位降落時(shí)上游壩坡穩(wěn)定性進(jìn)行分析，朱朋等[6]對三峽庫區(qū)滑坡體水位降落過程進(jìn)行計(jì)算，得出滑坡體安全系數(shù)隨水位降落不斷降低。但極限平衡法存在諸多弊端，例如需要假定滑裂面，無法模擬壩坡從局部破壞逐漸擴(kuò)展以至貫通成滑裂面的漸進(jìn)破壞過程等[7]。而有限元法則可較為真實(shí)地反映材料非線性本構(gòu)關(guān)系，實(shí)現(xiàn)各種復(fù)雜邊界條件的求解，其中強(qiáng)度折減法在推求安全穩(wěn)定性系數(shù)方面應(yīng)用廣泛。常見的有限元分析軟件ABAQUS、ANSYS無自帶水位變動(dòng)邊界設(shè)置[8]，多應(yīng)用于穩(wěn)定滲流分析中，對非穩(wěn)定滲流一般將庫水位降落劃分為緩降和驟降，仍按照穩(wěn)定滲流的方式進(jìn)行計(jì)算。對緩降認(rèn)為壩內(nèi)孔隙水壓力隨水位降落同步發(fā)生變化，對驟降則認(rèn)為壩內(nèi)孔隙水壓力隨水位降落不發(fā)生變化，這些都是對實(shí)際問題的一種假設(shè)，與實(shí)際情況不相符，例如劉博[9]按照緩降方式對比加固前、后的心墻壩內(nèi)浸潤線變化，Lane等[10]采用強(qiáng)度折減法確定水位緩降以及水位瞬時(shí)驟降這2種極端情況下的壩坡穩(wěn)定性。目前水位降落過程中壩內(nèi)各部位的變形和孔壓分布規(guī)律的研究較少,故有必要對非穩(wěn)定滲流開展瞬態(tài)有限元分析方法的研究。

2 計(jì)算原理

2.1 流固耦合計(jì)算原理

土壩上下游水位差產(chǎn)生的孔壓以體積力的形式，通過透水介質(zhì)作用于土體，改變壩內(nèi)應(yīng)力場。使壩內(nèi)土體發(fā)生變形，改變土體孔隙比和滲透系數(shù)，造成滲流場變化，二者相互作用，相互影響。耦合模型矩陣形式為

(1)

式中：K為土體剛度矩陣(kN/m3)；Δδ為位移增量(m)；ΔF為由外荷載引起的節(jié)點(diǎn)荷載增量(kN/m2)；ΔFS為變化的滲流體積力引起的節(jié)點(diǎn)荷載增量(kN/m2)；k(σij)為變化的滲透系數(shù)(m/s)；σij為有效應(yīng)力張量場(Pa)；H為水頭(m)；f為滲流場水頭分布函數(shù)。

2.2 強(qiáng)度折減法計(jì)算原理

在彈塑性土體材料計(jì)算中，通過不斷增大的折減系數(shù)Fr，對土體材料的黏聚力c和內(nèi)摩擦角φ進(jìn)行折減(式(2))。通過使土體單元在原荷載作用下因c、φ減小而發(fā)生塑性破壞，形成完整的貫通破壞帶。此時(shí)，對應(yīng)的折減系數(shù)為邊坡穩(wěn)定安全系數(shù)。

(2)

式中：c和c′分別為折減前后土體的黏聚力(kPa)；φ和φ′為分別為折減前后土體的內(nèi)摩擦角(°)；Fr為安全系數(shù)。

3 程序二次開發(fā)

3.1 水位變動(dòng)邊界設(shè)置

庫水位降落的非穩(wěn)定滲流場計(jì)算是以上游水位最高時(shí)的穩(wěn)定滲流場為基礎(chǔ)的。對土壩滲流應(yīng)力耦合分析時(shí)，在壩體和壩基的表面應(yīng)同時(shí)定義靜水壓力和設(shè)置孔壓邊界條件[11]。

在穩(wěn)定滲流計(jì)算中，滲流場處于動(dòng)態(tài)平衡過程，靜水壓力、孔隙水壓力改變只與邊界條件的變化相關(guān)，與時(shí)間無關(guān)，僅與垂向坐標(biāo)y呈線性相關(guān)。

即靜水壓力荷載為

F=γw(h(0)-y) 。

(3)

孔隙水壓力荷載為

U1=γw(h(0)-y) 。

(4)

在非穩(wěn)定滲流計(jì)算中，上游水位是時(shí)間t的函數(shù)，設(shè)為h(t)，靜水壓力、孔隙水壓力非線性變化，是時(shí)間t和垂向坐標(biāo)y的函數(shù)，設(shè)為F(t,y)、U1(t,y)。

靜水壓力荷載為

(5)

孔隙水壓力荷載為

U1=γw(h(t)-y) 。

(6)

式中：F為靜水壓力(kN/m2)；U1為孔隙水壓力(kN/m2)；γw為水體重度(kN/m3)；h(0)為上游初始水位(m)；h(t)為上游t時(shí)刻水位(m)；y為垂向坐標(biāo)變量(m)。

3.2 自動(dòng)處理庫水位下降參數(shù)插件程序的實(shí)現(xiàn)

由于ABAQUS中的LOAD模塊僅能施加作用點(diǎn)固定的靜態(tài)荷載，故本文采用Fortran和Python語言聯(lián)合開發(fā)ABAQUS用戶子程序，實(shí)現(xiàn)交互式的非線性動(dòng)態(tài)荷載的施加。該用戶子程序命名為UL插件程序，其可根據(jù)上游初始水位H、水位下降高度Δh和下降歷時(shí)t等水位降落參數(shù)，并基于DLOAD荷載子程序接口和DISP邊界子程序接口與LOAD模塊進(jìn)行實(shí)時(shí)對接，實(shí)現(xiàn)隨時(shí)間和空間位置發(fā)生變化的靜水壓力和孔隙水壓力的賦值功能。UL程序進(jìn)行了人機(jī)交互的可視化界面處理，可與ABAQUS深度兼容，只需在ABAQUS中啟用該程序，根據(jù)提示輸入不同方案下的水位降落參數(shù)，實(shí)現(xiàn)特定庫水位降落過程的滲流場和應(yīng)力場的計(jì)算。

4 模型驗(yàn)證

應(yīng)用上述模型，計(jì)算均質(zhì)土壩水位降落過程中的滲流場和應(yīng)力，分析壩體滲流特性和壩坡穩(wěn)定性，并與理論分析和已有研究結(jié)論進(jìn)行對照。

4.1 滲流場驗(yàn)證

選用模型壩高74 cm，壩頂寬度10 cm，壩底長306 cm，壩面坡比為1∶2，上游水位H1=60 cm，下游水位H2=0 cm的均質(zhì)土壩作為算例，進(jìn)行滲流場的模型驗(yàn)證，模型源于文獻(xiàn)[12]。

4.1.1 穩(wěn)定滲流對比結(jié)果

恒定庫水位下，對滲透系數(shù)為0.001 7 cm/s的模型壩體分別采用數(shù)值模擬和水力學(xué)三段法計(jì)算壩內(nèi)浸潤線，計(jì)算結(jié)果對比如圖1所示，模型計(jì)算浸潤線逸出點(diǎn)高度為26.94 cm，水力學(xué)法計(jì)算值為25.67 cm，模型計(jì)算結(jié)果較水力學(xué)法高4.43%，兩種方法計(jì)算所得浸潤線位置基本一致，證明采用本文模型對壩體進(jìn)行穩(wěn)定滲流場計(jì)算可行。

圖1 壩體穩(wěn)定滲流浸潤線對比

4.1.2 非穩(wěn)定滲流對比結(jié)果

分別對滲透系數(shù)為0.001 7、0.000 2 cm/s的兩種不同土質(zhì)的壩體水位降落過程進(jìn)行計(jì)算。隨著庫水位降落，壩體內(nèi)浸潤線隨之降落并出現(xiàn)分水嶺，分水嶺凸點(diǎn)高程逐漸降低，與已有成果得到的水位降落下浸潤線變化規(guī)律一致。陶同康[12]通過沙槽試驗(yàn)，劉潔等[13]、沙金煊[14]通過攝動(dòng)法分別對各時(shí)段自由面分水嶺凸點(diǎn)高程進(jìn)行計(jì)算，同本文數(shù)值模擬計(jì)算結(jié)果對比見圖2，各方法計(jì)算偏差見表1，說明模型對非穩(wěn)定滲流場計(jì)算合理。

表1 浸潤線凸點(diǎn)高程計(jì)算偏差表

圖2 非穩(wěn)定滲流浸潤線凸點(diǎn)高程對比

4.2 應(yīng)力場驗(yàn)證

模型選用壩高21.3 m，壩頂寬度7.3 m，基礎(chǔ)底厚7.3 m，上游水位H1=17.1 m，下游水位H2=0 m，壩體、壩基采用同種材料的均質(zhì)土壩[15]，示意圖見圖3。

圖3 土壩示意圖

計(jì)算上游水位自17.1 m降落至壩底過程中壩坡瞬態(tài)安全系數(shù)，0 h時(shí)刻壩體為17.1 m形成的穩(wěn)定滲流，以位移-折減系數(shù)曲線出現(xiàn)明顯拐點(diǎn)時(shí)為失穩(wěn)條件。采用強(qiáng)度折減法分別計(jì)算水位下降0、85 h和90 h時(shí)的下游壩坡安全系數(shù)。分別同賈蒼琴等[15]計(jì)算結(jié)果進(jìn)行比較，見表2。本文計(jì)算結(jié)果同極限平衡法平均偏差為0.99%，同賈蒼琴等有限元法平均偏差為4.85%，模型計(jì)算水位瞬態(tài)降落過程安全系數(shù)的方法較為合理。圖4、圖5分別為水位未發(fā)生降落時(shí)、水位降落85 h時(shí)刻滑動(dòng)面對比。

表2 3種方法的安全系數(shù)比較

圖4 水位未發(fā)生降落0 h時(shí)各方法計(jì)算的臨界破壞機(jī)制對比

圖5 水位下降85 h時(shí)各方法計(jì)算的臨界破壞機(jī)制對比

5 工程算例

5.1 工程概況

十里河水庫位于山西左云縣西南，該壩為均質(zhì)土壩，壩體全長為1 230 m，壩頂高程1 325.9 m，頂寬4 m，壩高16 m，上游壩坡1∶2.75。下游馬道以上壩坡1∶2.25，馬道以下壩坡1∶2.5。壩基上層為3 m厚的強(qiáng)透水的卵石混合土層，壩軸線上游馬道附近碎石層中加設(shè)坡比為1的黏土截水槽。根據(jù)碾壓土石壩規(guī)范，確定非常運(yùn)行條件水位變動(dòng)范圍，自校核洪水位1 325.54 m降至死水位1 318.2 m。具體工程布置剖面見圖6。

圖6 十里河水庫大壩橫剖面圖

5.2 模型設(shè)置

根據(jù)大壩剖面圖，模型范圍取壩體和壩基兩部分。壩基深度取20 m，水平方向向上、下游各取約30 m，總長150 m，采用CPE4P四節(jié)點(diǎn)單元，網(wǎng)格總數(shù)15 877個(gè)，節(jié)點(diǎn)16 255個(gè)，如圖7。壩基約束左右兩側(cè)水平位移，底部同時(shí)約束水平、豎向位移，上下游壩坡、河床和壩基兩側(cè)均為透水邊界，壩基底部為不透水邊界。模型整體施加重力荷載，上游壩坡施加非線性的靜水壓力荷載及孔隙水壓力邊界。模型分區(qū)中①為壩體素填土區(qū)；②為強(qiáng)透水的卵石混合土層；③為截水槽黏土層；④為壩基素填土層，具體材料參數(shù)見表3。

圖7 十里河水庫大壩模型設(shè)置及分區(qū)示意圖

表3 材料參數(shù)

5.3 計(jì)算結(jié)果分析

5.3.1 浸潤線分析

上游水位自1 325.54 m分別以速率0.05、0.10、0.50、1.00、1.50 m/d下降146.8、73.4、14.68、7.34、4.89 d至1 318.2 m，取初始時(shí)刻、水位降落完成時(shí)和水位每下降1.5 m(時(shí)間間隔分別為30、15、3、1.5、1 d)典型時(shí)刻浸潤線，繪制成浸潤線隨時(shí)間變化圖，見圖8。水位下降過程中滲流場3個(gè)主要特征：①壩體中浸潤線的下降同庫水位下降相比，明顯滯后。壩前水位下降速度越大，滯后現(xiàn)象越明顯。②壩體上游臨水側(cè)浸潤線呈凸形，出現(xiàn)分水嶺且降速越大，浸潤線上曲率越大。③浸潤線最高點(diǎn)位置隨著下降高度的增加向下游移動(dòng)。這些特征均對上游壩坡不利，增大上游壩坡下滑的趨勢。

圖8 不同降水速率下壩體浸潤線隨時(shí)間變化

5.3.2 速度矢量分析

滲流場對壩體影響主要在上游部位，圖9為上游水位以0.5 m/d速率下降至1 323.54、1 320.54、1 318.2 m時(shí)上游壩坡速率矢量圖。在水位降落初期，壩體內(nèi)孔隙水壓力的消散速度與水位的降落速度相適應(yīng)如圖9(a)、圖9(b)所示，表現(xiàn)為浸潤面降落相比于庫水位降落稍有滯后，仍為從上游流向下游。隨庫水位持續(xù)降低，浸潤面滯后現(xiàn)象顯著并產(chǎn)生分水嶺，當(dāng)水位降至1 320.54 m，如圖9(c)所示。臨水側(cè)流場轉(zhuǎn)變?yōu)榭拷嫌尾糠值乃鞒蛏嫌瘟鲃?dòng)，而其余部分仍朝下游流動(dòng)，當(dāng)水位持續(xù)下降現(xiàn)象更為明顯，如圖9(d)所示。

圖9 0.5 m/d降水速率下典型水位壩體流速矢量圖

5.3.3 瞬時(shí)抗滑穩(wěn)定性分析

上游水位分別以0.05、0.10、0.50、1.00、1.5 m/d速度下降時(shí)，上、下游壩坡的安全系數(shù)隨上游水位變化關(guān)系如圖10所示，降水完成時(shí)不同降速的上、下游壩坡安全系數(shù)見表4。

圖10 不同降速下壩體上下游安全系數(shù)

表4 降水完成時(shí)上、下游邊坡安全系數(shù)

圖10、表4表明：

(1)穩(wěn)定滲流狀態(tài)下，十里河水庫大壩強(qiáng)透水層結(jié)合黏土截水槽的工程措施，能有效降低壩體內(nèi)浸潤線、逸出點(diǎn)高程，減少壩內(nèi)滲徑，加快壩內(nèi)孔隙水壓力的消散，增大下游壩坡安全系數(shù)，從而提高壩體的整體穩(wěn)定性，使得下游壩坡安全系數(shù)受庫水位降落影響不大。

(2)庫水位降落完成時(shí)，下游壩坡安全系數(shù)稍有增大，上游壩坡安全系數(shù)持續(xù)減小且水位下降速度越大安全系數(shù)越低。庫水位在1 325.54 m未降落時(shí)，上游壩坡安全系數(shù)為3.685，下游壩坡為2.459，降水完成后，上游壩坡安全系數(shù)降低率為30%左右，下游壩坡變化率僅為0.3%～2%左右，當(dāng)降速為1.5 m/d時(shí)，上游壩坡安全系數(shù)降低率為33.81%，下游壩坡增加僅為0.25%。

(3)當(dāng)上游水位降落4.5 m至H=1 321.04 m時(shí)，如圖10虛線所示，水位下降速度越大，降落至同一位置時(shí)，上游壩坡安全系數(shù)越小。當(dāng)降速為1.5 m/d時(shí)，此時(shí)上游壩坡安全系數(shù)為2.69，降速為0.05 m/d時(shí)，安全系數(shù)為3.009。

(4)水位降落初期，上游壩坡安全系數(shù)隨水位下降速度的增加而大幅減小，隨著降落時(shí)間的持續(xù)，安全系數(shù)變化幅度受水位下降速度的影響減小。

(5)在初始水位形成的穩(wěn)定滲流情況下，壩體下游坡安全系數(shù)較小，較為危險(xiǎn)。水位下降速度較小時(shí)，上游壩坡安全系數(shù)遠(yuǎn)大于下游壩坡，壩體的整體安全系數(shù)由下游壩坡決定。當(dāng)水位下降速度增大，隨降水歷時(shí)的增加，壩體的整體安全系數(shù)逐漸轉(zhuǎn)為由上游壩坡決定。如當(dāng)以1.5 m/d降至1 318.6 m附近，上游壩坡安全系數(shù)小于下游壩坡，壩體整體穩(wěn)定性由上游壩坡決定。

6 結(jié) 論

本文旨在結(jié)合算例探討二次開發(fā)的庫水位降落土壩模型的可靠性，并對水位降落瞬態(tài)過程中壩體滲流及穩(wěn)定性進(jìn)行分析，結(jié)果如下。

(1)借助Fortran語言編寫UL子程序，Python 語言對UL子程序進(jìn)行可視化處理，實(shí)現(xiàn)了交互式的非線性動(dòng)態(tài)靜水壓力和孔隙水壓力荷載的施加。

(2)利用二次開發(fā)模型，對壩體滲流和穩(wěn)定性進(jìn)行計(jì)算，對比理論分析和已有研究成果中浸潤線位置、壩坡瞬態(tài)安全系數(shù)，結(jié)果均較為合理，驗(yàn)證了模型的可靠性。

(3)針對十里河水庫大壩進(jìn)行計(jì)算分析得出：①庫水位降落時(shí)，浸潤線降落較庫水位降落存在滯后現(xiàn)象且出現(xiàn)分水嶺，降落速度越大，浸潤線降落的滯后現(xiàn)象越明顯；②隨水位降落歷時(shí)的增加，臨水側(cè)流場自上游流向下游轉(zhuǎn)變?yōu)榭拷嫌尾糠值乃鞒蛏嫌瘟鲃?dòng)；③在降水初期，上游坡安全系數(shù)降幅隨水位下降速度的增大而增大，隨降水歷時(shí)的增加，上游安全系數(shù)降幅受水位下降速度的影響減小。且降落至任意水位時(shí)，水位下降速度越大上游壩坡安全系數(shù)越低。水位自1 325.54 m降至1 318.2 m時(shí)，0.05 m/d降速時(shí)上游坡安全系數(shù)減小25.06%，1.5 m/d降速時(shí)減小33.81%。

(4)二次開發(fā)的UL子程序拓展了ABAQUS有限元法在非穩(wěn)定滲流中的應(yīng)用，可將其推廣應(yīng)用于各類降水工程中非穩(wěn)定滲流場、應(yīng)力場的計(jì)算。