基于有限元法和響應面法的閘基滲流反演分析

劉麗君，范雪楓，吳震宇

(1.四川大學 水力學與山區河流開發保護國家重點實驗室 水利水電學院，四川 成都 610065；2.昆明勘測設計研究院，云南 昆明 650216)

我國西南山區大多數水閘建在覆蓋層地基上。這些水閘的地基在材料組成、滲流特性、防滲方式和滲流監測方面具有一定的相似性。山區河流的河床覆蓋層普遍較為深厚，結構和級配變化大，局部存在架空現象，透水性強，滲透穩定性差。滲壓監測數據可以真實反映閘基滲流性態，但缺點是測點數量有限，只能觀測局部滲流情況。采用有限元法等數值模擬技術進行滲流分析，可以比較真實地模擬地基材料組成、各種材料的滲流特性、防滲設施和滲流邊界條件等，得到整個計算模型區域內的滲流場，因此，與基于監測數據方法相比，數值模擬更加全面地分析地基滲流狀態[1- 10]。滲透系數等材料滲流參數的取值對地基滲流數值模擬結果影響較大，如果參數選取不合理，則滲流計算結果與實際滲流情況會存在較大差異。根據滲流監測數據率定材料的滲流參數(即滲流參數反演分析)，然后采用率定后的參數進行滲流數值模擬，可使數值模擬結果更接近實際情況。

某水閘位于四川岷江干流，主要功能為發電，地基主要由漂卵石層和砂壤土層相間組成，閘基滲流安全是該工程重要問題之一。在閘基防滲墻后鉆孔埋設了8個滲壓計，用于監測閘基滲壓變化情況。本文首先根據該水閘2010—2016年期間的滲壓監測數據，分析閘基滲壓變化規律及其原因。然后根據滲壓監測數據，采用基于響應面法的參數反演分析方法，對覆蓋層材料的滲透系數進行率定。最后采用率定后的滲透系數，對枯期和汛期的閘基滲流場進行有限元數值模擬，驗證了滲壓變化的原因，并對閘基整體滲流情況進行了分析，獲得了監測數據未能反映的閘基滲透坡降情況，為閘基滲流安全評估提供了參考。

1 閘基滲壓監測數據分析

1.1 工程概況

水閘于20世紀70年代建成。該水閘從左岸至右岸由取水口攔污柵閘、攔河閘、擋水壩等建筑物組成。攔河閘橫河向長度為156m，順河向長度為25m，最大閘高21.4m，設有6孔泄洪沖沙閘。水庫總庫容93萬m3，正常蓄水位945.00m，設計洪水位945.60m，校核洪水位948.87m。閘基覆蓋層最大深度45m，采用上游混凝土鋪蓋和懸掛式防滲墻防滲，防滲墻位于閘室上游20m，厚度0.7m，最大深度27m。閘室上游鋪蓋長75～95m，閘室下游護坦長56m。閘址處兩岸岸坡較緩，臨河坡高300～500m，河谷較開闊，河床寬80～100m。閘址處河床覆蓋層主要由漂卵石層和砂壤土層相間組成，厚度約40～62m。1#閘孔順河向剖面如圖1所示。

圖1 1#閘孔順河向剖面圖

為監測閘基滲壓，以及上游鋪蓋和懸掛式防滲墻的防滲效果，在防滲墻后沿壩軸線方向鉆孔埋設了兩排滲壓計，其中上游側布置了5支滲壓計(編號為UP1～UP5)，下游側布置了3支滲壓計(編號為UP6～UP8)，總共8支滲壓計。閘基滲流監測布置如圖2所示。

圖2 閘基滲壓監測平面布置圖

1.2 閘基滲壓監測數據分析

閘基滲壓分析采用2010年1月—2016年12月期間的監測數據。為了便于和庫水位進行比較，滲壓測值已轉化為滲壓水頭高程。各測點滲壓水頭高程隨時間變化規律一致，受篇幅限制，只給出了UP2測點的滲壓測值過程線，如圖3所示。

圖3 閘基滲壓實測過程線

從圖3可以看出，滲壓測值呈明顯周期性變化規律，即汛期(5—10月)水位高，枯期水位低。各滲壓計的最大測值基本出現在7月，在938.34～945.15m之間(受2013年7月9號泥石流影響，下游水位壅高超945m高程，各滲壓計測值也相應偏高)；最小測值基本出現在12月～次年4月，在935.07～939.34m之間。滲壓測值與庫水位的相關性弱，如2015—2016年滲壓測值的年變幅為2.74～4.25m，而相應庫水位的年變幅僅為1.09～1.64m，因此，庫水位變化不是引起滲壓測值周期性變化的主要原因。

在枯期，攔河閘關閉，通過左岸取水口引水發電，下游河道基本無水。在汛期，開閘泄洪，下游河道水位較高。地基采用上游混凝土鋪蓋和懸掛式防滲墻防滲，防滲墻位于閘室上游20m。從地基防滲布置來看，其主要作用是控制上游庫水向下游滲透，而不能有效控制下游河水向上游滲透。所以，滲壓測值的周期性變化估計是汛、枯期下游河道水位大幅變化而行成較大水位差變化所引起的。

采用靠近防滲墻的一排滲壓計(UP1～UP5)的枯期測值(上、下游水位差較大)計算閘基滲壓折減系數

α=(HSP-HRB)/(HR-HRB)

(1)

式中，HSP—滲壓水頭高程；HR—庫水位；HRB—下游河床高程。

計算滲壓折減系數基本在0.1～0.3之間，表明閘基防滲系統的防滲效果較好。

由于沒有下游河道水位測值，故滲壓測值周期性變化的原因還需進一步驗證。此外，滲壓計的數量有限，只能觀測局部滲流情況，且只能獲得滲壓測值，無法確定閘基滲透坡降。因此，采用有限元法進行閘基滲流分析。為了使有限元計算結果更接近實際情況，先根據滲壓監測數據率定閘基覆蓋層材料的滲透系數。

2 閘基覆蓋層滲透系數反演分析

2.1 閘基滲流有限元模型

采用ANSYS軟件[11]建立1#閘孔和2#閘孔的閘基平面滲流有限元模型。模型包含閘底板、上游鋪蓋、下游護坦(含排水孔)、防滲墻、閘基覆蓋層及下部花崗巖，如圖4—5所示。護坦中的排水孔通過增大相應位置處的單元的滲透系數的方式模擬。有限元模型的模擬范圍為：自上游鋪蓋末端向上游延伸75m，自下游護坦末端向下游延伸80m，自覆蓋層底部向下延伸50.4m。為了對比計算滲壓和實測滲壓，根據滲壓計的埋設位置，在有限元模型中相應位置布置結點(1#閘孔閘基中的滲壓計為UP2和UP6，2#閘孔閘基中的滲壓計為UP3和UP7)。

由于滲流分析和熱分析的基本方程相似，因此采用ANSYS軟件的熱分析模塊計算閘基滲流。采用PLANE55熱分析單元離散有限元模型，1#閘孔閘基的有限元模型劃分為6922個單元和7095個結點，2#閘孔閘基的有限元模型劃分為6354個單元和6527個節點。模型上下游邊界及底部邊界為不透水邊界。模型頂部邊界水頭取枯期和汛期的上、下游水位，上游水位采用實測值，枯期下游按無水考慮，汛期的下游水位根據水力學計算結果確定。

2.2 滲透系數反演分析方法

混凝土閘底板、混凝土鋪蓋、混凝土護坦、混凝土防滲墻和覆蓋層下部花崗巖的滲透系數均有試驗值。因為難以取到閘基覆蓋層原狀材料，以及材料顆粒尺寸大難以進行室內試驗，所以，閘基覆蓋層材料滲透系數的設計值是通過現場勘察并結合工程經驗確定的，具有不確定性。而材料的滲透系數對閘基滲流計算結果影響較大，因此，根據滲壓監測數據對覆蓋層材料的滲透系數進行反演分析。

滲透系數的反演分析即尋找一組滲透系數使有限元計算的閘基滲壓與滲壓計測到的滲壓最吻合。需要采用滲壓的計算值和實測值建立目標函數，求解一個以滲透系數為未知變量的優化問題。采用有限元法計算滲壓，既耗時又無法建立目標函數，給優化問題的求解造成困難。本文采用滲壓計算值的響應面方程(即以滲壓計算值為因變量和需要率定的滲透系數為自變量的函數)代替有限元計算，可以建立目標函數，提高參數反演分析的效率。

2.2.1滲壓計算值的響應面方程

本文通過構造如下響應面方程[12- 13]，代替耗時的有限元計算，快速確定不同滲透系數對應的滲壓計算值，極大提高參數反演分析的效率。

(2)

對需要率定的滲透系數進行組合。然后對每一組滲透系數，采用有限元法進行閘基滲流計算，提取滲壓計埋設處的滲壓計算值。最后通過回歸分析確定響應面方程(2)的系數。

本文采用正交試驗設計[14- 15]進行滲透系數的組合。以三種材料的滲透系數組合為例進行說明。每種材料的滲透系數取3個值，即μ、μ(1-kδ)、μ(1+kδ)，其中μ為滲透系數均值，δ為滲透系數變異系數，k為系數。然后根據正交試驗表，對3種材料的9個滲透系數進行組合，得到9個試驗點，如圖6中立方體網格中的黑色圓點所示。由圖6可見，立方體網格共有6個平面，每個平面上都有三個試驗點，每條直線上都有1個試驗點，并且這些試驗點是均勻分散的，能很好地代表大量試驗點的情況。

圖6 正交設計示意圖

2.2.2滲透系數反演分析的目標函數與求解

利用滲壓計算值的響應面方程(1)和滲壓的實測值建立如下滲透系數反演分析的目標函數：

(3)

(4)

滲透系數的反演分析就是在各種材料的滲透系數的上、下限范圍內搜尋一組滲透系數使式(3)表達的目標函數值最小。因此，滲透系數的反演分析可以轉化為帶約束條件的非線性規劃問題。本文利用LINGO軟件[16]求解式(3)—(4)表達的非線性規劃問題。

2.3 覆蓋層滲透系數反演分析結果

根據UP2、UP3、UP6和UP7四個滲壓計對2015年汛期和枯期的滲壓監測數據，采用2.2節所述方法對漂卵石、卵礫石和砂壤土三種材料的滲透系數進行反演分析。需要率定滲透系數的材料有三種，故滲壓計算值的響應面式(2)共有10個系數。采用回歸分析確定這些系數需要10個以上的試驗點。因此，進行兩次正交試驗設計，共得到18個試驗點。第一次正交試驗設計每種材料的滲透系數取μ、μ(1-δ)、μ(1+δ)3個值，第二次取μ、μ(1-3δ)、μ(1+3δ)3個值，其中μ采用滲透系數的設計值，見表1，δ取0.3[17]。采用SPSS軟件[18]的正交試驗設計模塊對漂卵石、卵礫石和砂壤土三種材料的滲透系數進行組合。然后對每一組滲透系數，采用圖4—5所示的有限元模型進行汛期和枯期兩種工況下的閘基滲流計算，提取UP2、UP3、UP6和UP7四個滲壓計埋設處的滲壓計算值。最后利用SPSS軟件的回歸分析模塊確定汛期和枯期工況下UP2、UP3、UP6和UP7四個測點處滲壓計算值的響應面方程(總共8個)系數。8個響應面方程的回歸分析的復相關系數均大于0.9，同時采用響應面方程計算的滲壓和有限元計算結果的差別很小，最大相對誤差小于0.34%，大部分相對誤差小于0.1%。

表1 材料滲透系數的設計值和率定值

利用滲壓計算值響應面方程和滲壓實測值，按照式(2)—(3)建立滲透系數反演分析目標函數，滲透系數的上限設為μ(1+3δ)，下限設為μ(1-3δ)。采用LINGO軟件搜索漂卵石、卵礫石和砂壤土滲透系數的最優組合，使目標函數值最小，結果見表1。由表1可見，漂卵石和卵礫石滲透系數的率定值與設計值有較大差別。采用覆蓋層材料滲透系數的率定值進行閘基滲流計算。圖7顯示了UP2、UP3、UP6及UP7四個測點處的滲壓計算值和實測值的對比。可以看出，滲壓的計算值和實測值總體比較接近。計算值與實測值的差異可能是滲壓計的測量誤差和有限元模型誤差(沒有采用三維模型，此外，用于建模的閘基地質剖面是根據有限鉆孔數據繪制的，與實際地層情況可能存在差別)造成的。

圖7 汛期和枯期滲壓計算值與實測值的對比

3 閘基滲流有限元分析

3.1 上下游水位變化對閘基滲壓的影響分析

根據1.2節的分析，汛、枯期下游河道水位大幅變化是引起滲壓計測值周期性變化的主要原因。本節采用有限元數值模擬，研究上下游水位變化對閘基滲壓的影響，驗證滲壓計測值周期性變化的原因。首先將上游水位固定在944.5m，將下游水位從936m以0.5m的間隔升高至938.5m，計算UP2、UP3、UP6及UP7四個測點處的滲壓，分析下游水位變化對滲壓的影響。然后將下游水位固定在936m，將上游水位從942.5m以0.5m的間隔升高至945m，計算各測點處的滲壓，分析上游水位變化對滲壓的影響。測點處滲壓隨上下游水位的變化情況如圖8所示，可以看出滲壓隨上游或下游水位的升高而增大，并呈線性關系。滲壓隨下游水位變化的斜率明顯大于隨上游水位變化的斜率，上游水位升高2.5m引起的滲壓增大值不超過1m，而下游水位升高2.5m引起的滲壓增大值為2～2.5m。因此，說明滲壓計測值的變化主要受下游水位變化影響。

圖8 上下游水位對滲壓的影響

3.2 閘基滲流特性分析

采用圖4—5所示的有限元模型(其中覆蓋層材料滲透系數采用率定值)進行汛期(上游水位944.60m，下游水位938.50m)和枯期(上游水位944.59m，下游水位936.00m)兩種工況下的閘基滲流計算。由于枯期和汛期兩種工況下閘基滲壓分布規律相似，只是量值不同。故僅給出枯期時1#閘孔的閘基滲流等勢線和滲透坡降，分別如圖9—10所示。

圖9 枯期1#閘孔閘基滲流等勢線

圖10 枯期1#閘孔閘基滲透坡降

由圖9—10可以看出，閘基滲壓從上游至下游逐漸降低，且在防滲墻兩側滲壓下降較快，說明防滲墻發揮了降低閘基滲壓的作用。閘基三個部位的滲透坡降相對較大，即防滲墻底端附近、混凝土鋪蓋上游端附近和護坦上游排水孔附近，其中以防滲墻底端附近的滲透坡降最大。由于防滲墻靠上游布置，不能控制下游河水向上游滲透，因此，汛期下游水位大幅升高后，防滲墻后的滲壓比枯期時明顯增大。由于枯期上下游水位差較大，因此，相應的閘基滲透坡降也較大。上述三個部位在汛期的最大滲透坡降分別為0.302、0.129和0.148，枯期分別為0.425、0.181和0.209。除上述幾個部位外，閘基其余部分的滲透坡降都很小，閘基滲流總體比較穩定。由有限元計算結果可見，防滲墻底端、混凝土鋪蓋上游端和護坦上游排水孔附近的閘基土體發生滲透破壞的風險較大。特別是在枯期運行時，需要加強巡視檢查，觀察是否存在閘前水庫旋渦、護坦排水孔中析出細顆粒等異常現象。到目前為止，該水閘閘基無異常滲流現象。

4 結論

與基于監測數據的方法相比，數值模擬可以分析地基整體滲流狀態。根據滲壓監測數據率定材料滲透系數，再采用率定后的滲透系數進行滲流計算，可使計算結果更接近實際情況。

本文提出基于響應面法的滲透系數反演分析方法。通過構造滲壓計算值的響應面方程，代替耗時的有限元計算，而且可以建立顯式表達的目標函數，因此，能夠提高滲透系數反演分析的效率。

結合監測數據和數值模擬對水閘的地基滲流進行了分析。監測數據表明該水閘地基滲壓呈汛期高、枯期低的周期性變化規律，原因是汛、枯期下游河道水位大幅變化。根據滲壓監測數據，對閘基覆蓋層材料滲透系數進行了反演分析。采用率定后的滲透系數，對枯期和汛期的閘基滲流場進行了有限元數值模擬，驗證了滲壓周期性變化的原因，并揭示了閘基中滲透坡降較大的部位。本文提出的方法可用于類似水閘地基的滲流分析與安全評估。