戴前偉， 孔重陽， 雷 軼， 張 彬， 韓行進

(1.中南大學 地球科學與信息物理學院， 湖南 長沙 410012； 2.中南大學 有色金屬成礦預測與地質(zhì)環(huán)境監(jiān)測教育部重點實驗室， 湖南 長沙 410012； 3.五凌電力有限公司， 湖南 長沙 410004)

1 研究背景

水利工程中水壩的修建能調(diào)節(jié)控制下泄水量，從而降低洪澇等災害帶來的影響，并且提高水資源的利用率，滿足人們生產(chǎn)和生活的需要。無壓滲流問題在巖土工程、水利工程、地下水等領域廣泛存在，滲流或者孔隙壓力過大會導致水壩失穩(wěn)[1]。土石壩作為最常見的一種水壩類型，其壩體具有一定的透水性，在水的滲流作用下容易發(fā)生滲漏而遭到破壞[2-4]，因此需要對土石壩的滲流安全問題進行更深入的研究。

目前，壩體的工程設計、施工、檢測等常常要求探測出自由面的位置，從而對壩體的滲漏情況進行初步判斷并提前做好相應的防滲措施，然而自由面的位置和形狀是事先未知的，并且在求解過程中會不斷發(fā)生變化[5]，如何能高效準確地確定自由面以及滲流溢出點的位置正是滲流分析中需要解決的關鍵問題[6]。滲流分析中最常用的方法是數(shù)值模擬，不同的數(shù)值模擬方法已經(jīng)被國內(nèi)外的水文地質(zhì)學者們應用于相應的研究之中[7]。Neumann[8]率先利用固定網(wǎng)格的有限元法對滲流自由面問題進行了求解，其主要特點是在求解滲流自由面的迭代過程中不對原始單元網(wǎng)格進行修改和調(diào)整。付延玲等[9]將整個求解域的網(wǎng)格單元劃分為3種類型，在處理過渡單元的滲透矩陣時，引入了雙參數(shù)罰函數(shù)，利用流量等效法對滲透系數(shù)進行調(diào)整，提高了計算過程的效率以及結(jié)果的精度。唐紅影[10]對棄單元法進行改進，使計算域逐漸靠近實際滲流區(qū)域，并對復合單元進行二次細分，降低了網(wǎng)格單元對計算精度的影響。Darbandi等[11]提出的有限體積法確保了單元的質(zhì)量守恒，但要求網(wǎng)格與自由面相匹配。Chugh等[12]的邊界元法只對域邊界進行離散，通過移動邊界節(jié)點可以簡單地改變域的幾何結(jié)構(gòu)，雖然可以降低問題的維數(shù)，簡化求解，但對于非均勻區(qū)域和非線性問題，邊界元法很難求解。Zheng等[13]將數(shù)值流形方法與無網(wǎng)格Galerkin 方法相結(jié)合，對非均質(zhì)壩體的自由面進行了數(shù)值模擬，通過增加積分點支撐域的場節(jié)點，計算精度得到了很大提高。Zhang 等[14]用移動Kriging 無網(wǎng)格法和蒙特卡羅積分法分析了無約束滲流問題，但在迭代過程中會移動場節(jié)點，這給有限元分析中自適應網(wǎng)格法的使用帶來了一定的困難。Dai等[15]提出了一種新的綜合方法，將自適應網(wǎng)格策略與Galerkin有限元法相結(jié)合，對自由曲面的位置進行搜索，取得了很好的效果。

本文將光滑有限元法應用于無壓滲流問題的分析中，建立二維無壓滲流問題的光滑有限元模型，在飽和滲流區(qū)域計算求解自由面的位置。通過對矩形均質(zhì)壩和直角梯形均質(zhì)壩經(jīng)典算例的計算，并與其他數(shù)值方法的結(jié)果進行對比來驗證本文方法的有效性。同時模擬了土石壩中不同的滲透異常體對滲流速度、流體壓力、水頭分布等滲流參數(shù)的具體影響，為之后的滲流反演提供相應的依據(jù)。

2 無壓滲流問題的數(shù)學描述

為了方便描述無壓滲流問題，以圖1為例示意土石壩中的無壓滲流。圖1中的問題域被自由面AE分為兩部分，自由面以上的部分ΩU為不飽和域，以下的部分ΩS為飽和域。假設問題域中滲流僅僅流經(jīng)飽和區(qū)域。無壓滲流分析中最重要的問題即為自由面的求解，只要尋找到自由面的位置，飽和域與不飽和域便可區(qū)分開來。

圖1 土石壩壩體無壓滲流示意圖

在壩體飽和域中，任意一點的水頭值的數(shù)學表達式為：

(1)

式中：φ為滲流測壓水頭，m；p為流體壓力，Pa；ρ為流體密度，kg/m3；g為重力加速度，m/s2；y為豎直方向坐標，m。

飽和域中流體連續(xù)性方程為：

(2)

根據(jù)達西定律，飽和域中的滲流速度v為：

(3)

式中：vx、vy分別為水平向、豎直向的滲流速度，m/s；kx、ky分別為水平向、豎直向的滲透系數(shù)，m/s。

將公式(3)代入公式 (2)中，可得到二維無壓穩(wěn)定滲流的微分方程[16]：

(4)

在圖1中，相應邊界條件如下：

(1) 上下游庫水AB、尾水CD水頭邊界條件：

φAB=H1

(5)

φCD=H2

(6)

(2) 底部邊界BC為不透水邊界，其外法向方向流量qn為0，即滿足流量邊界條件：

qn=v·n=0

(7)

式中：n為單位外法向向量。

(3) 自由面AE需同時滿足兩個邊界條件：

qn=0,φ=y

(8)

(4) 滲流面DE邊界條件：

qn≤0,φ=y

(9)

即在滲流面上有流量流出。

3 無壓滲流場的光滑有限元模型

本部分簡要介紹應用在無壓滲流問題中的光滑有限元法。利用固定網(wǎng)格將整個問題域進行離散化，每個節(jié)點的水頭值可以通過插值得到：

uh(x)=N·H

(10)

式中：N為形函數(shù)矩陣；H為網(wǎng)格節(jié)點的水頭值。

求解無壓滲流問題控制方程的過程中，需要得到水頭梯度，而該值的獲得需要計算出形函數(shù)梯度。在經(jīng)典有限元法中，形函數(shù)梯度矩陣由形函數(shù)求導得到，而在光滑有限元法中，利用梯度光滑技術得到形函數(shù)梯度矩陣：

(11)

本文中采用的光滑函數(shù)為Heaviside-type型函數(shù)：

(12)

利用高斯散度定理(或者是分部積分法)，將面積分轉(zhuǎn)化為線積分,并將公式(12)代入公式(11)得：

(13)

由公式(13)可知，光滑梯度技術將原先二維的面積分轉(zhuǎn)換成沿著邊界進行的線積分。將公式(13)代入節(jié)點水頭公式(10)中得到：

(14)

(15)

其中：

(16)

將控制方程和邊界條件改寫為積分弱形式，并將水頭插值代入得到無壓滲流問題的Garlerkin弱形式：

(17)

(18)

弱形式中，R與自然邊界條件相關，由前文邊界條件可知邊界沿外法向方向流量為0，即R=0[17]。

利用固定網(wǎng)格法求解無壓滲流問題，可以降低對網(wǎng)格的依賴性，并且不需要在迭代過程中進行網(wǎng)格的修改。由于網(wǎng)格是固定的，每次迭代過程中得到的自由面會將整個求解域分為上下兩部分，產(chǎn)生3種類型的光滑單元，如圖2所示，自由面以上的單元稱為外部光滑單元，以下的稱為內(nèi)部光滑單元，自由面經(jīng)過的單元稱為相交光滑單元。由于采用的是四邊形單元，自由面經(jīng)過相交光滑單元時會將該單元切分成不同的形狀，可能出現(xiàn)的切分形狀如圖3所示。

圖2 3種類型的光滑單元示意圖

圖3 邊界相交單元被自由面可能切分的形狀

(19)

式中：Ω1為內(nèi)部光滑單元區(qū)域；Ω2為相交光滑單元區(qū)域。

在每個單元上運用梯度光滑技術得到單元的水頭值梯度，公式(19)可改寫為單元積分和的形式：

(20)

在計算過程中，認為每個網(wǎng)格中的滲透系數(shù)為常數(shù)，而且光滑形函數(shù)梯度矩陣也為常數(shù)矩陣，則總體剛度矩陣為：

(21)

4 迭代求解算法

利用光滑有限元法求解無壓滲流問題，根本目的是尋找自由面的位置，根據(jù)自由面滿足的邊界條件，不斷在迭代過程中更新自由面位置和形狀，直至達到一定的收斂條件。具體的迭代過程如圖4所示。

圖4 自由面迭代求解過程示意圖

(1) 選擇1個線段作為初始自由面，如圖4中AD。自由面的溢出點需要在CE上，而初始點位于BF上。

(2) 在CE上選擇一點D作為自由面溢出點的初始解。

(3) 連接點A與D，形成的線段AD作為自由面的初始解。

(4) 將初始自由面AD切分成一系列控制點H1，H2，…，Hn。選擇如圖4所示的BC線段做為基線，并將基線BC劃分為一系列基點G1，G2，…，Gn。

(5) 分別連接H1點與G1點，H2點與G2點，…，Hn點與Gn點，得到一系列射線H1G1、H2G2、…、HnGn。在迭代過程中，控制點H1，H2，…，Hn組成了自由面，并且這些控制點會沿著相應的射線移動，形成新的自由面。

上述準備工作完成后，利用假設的初始自由面AD解上述方程，可以得到求解域各個節(jié)點的水頭值。在這一步驟中，控制點的位置必須沿著射線進行修改以滿足自由面上的邊界條件?？刂泣c的垂直坐標表示為：

yinew=yiold+λσi

(22)

σi=hi-yiold

(23)

式中：yiold和yinew分別為第i個控制點在前一次和下一次迭代的高程值；hi為在控制點處計算得到的測壓水頭值；λ為給定的步長，本文取λ=0.8。若步長過小，則迭代次數(shù)增多，計算消耗較大；若步長過大，則難以到達收斂狀態(tài)。在能達到收斂狀態(tài)的情況下，選擇適合的步長來減少迭代次數(shù)，降低計算消耗。每次進行迭代前進行如下判斷：

(24)

式中：n為控制點個數(shù)；ε為非常小的參數(shù)，本文中取值為ε=10-3，該值越小，說明計算結(jié)果中控制點水頭值與高程值越接近，但迭代次數(shù)會增加、計算消耗增大，也有可能出現(xiàn)迭代難以收斂的問題；該值越高，則計算結(jié)果的誤差越大，難以滿足計算精度的要求。如果判斷式(24)得以滿足，則迭代過程終止，否則根據(jù)公式(22)、(23)對控制點位置進行更新。

5 算例分析

算例分析即利用本文提出的算法來求解無壓滲流問題經(jīng)典模型和異常體模型的自由面位置，通過與其他方法[11,15,18-20]在標準模型應用得到的結(jié)果進行對比來驗證本文算法的有效性，在異常體模型中研究異常體的存在對滲流速度、水頭分布、流體壓力等滲流場參數(shù)的影響。

5.1 矩形均質(zhì)壩模型(示例1)

示例1采用的經(jīng)典模型為高11 m，寬5 m的矩形均質(zhì)壩模型。上游庫水的水頭值H1=10 m，下游水頭值H2=2 m。根據(jù)文獻[17]、[21]中選取的滲透系數(shù)值，本部分滲透系數(shù)設置為1 m/s，底部邊界設置為不透水邊界。所采用的網(wǎng)格如圖5(a)所示，初始自由面和計算得到的自由面位置如圖5(b)所示。假設該模型的滲流區(qū)域為均質(zhì)各向同性介質(zhì)，初始自由面可以由Dupuit公式得到：

圖5 示例1中網(wǎng)格尺寸及初始自由面與最終自由面結(jié)果對比

(25)

式中：H1、H2分別為上游庫水、尾水的水頭值，m；L為底部邊界的長度，m。

示例1中本文方法對自由面計算結(jié)果與其他方法計算結(jié)果的對比如圖6(a)所示，本文方法計算得到的水頭等值圖、滲流速度圖及流體壓力等值圖如圖6(b)～6(d)所示。由圖6(a)可看出，本文方法計算得到自由面位置與其他方法得到的自由面位置具有很好的一致性，因此本文方法可以很好地預測自由面以及溢出點的位置。圖6(d)中的流體壓力可通過公式(1)計算得到， 其零壓力等值線與自由面位置相吻合。

圖6 示例1中自由面本文方法與其他方法計算結(jié)果對比及滲流場各參數(shù)計算結(jié)果

5.2 具有單異常體矩形壩模型(示例2和示例3)

在土石壩中，由于介質(zhì)滲透系數(shù)的不同，自由面的位置會受到一定的影響[21]。實際中常常需要根據(jù)自由面的位置和形狀推測滲透異常區(qū)域的位置，因此研究異常體對自由面所產(chǎn)生的影響是十分有必要的。示例2和示例3的模型尺寸、邊界條件和網(wǎng)格與示例1相同，但在模型中加入了不同滲透系數(shù)的異常體，將具有不同異常體位置(近自由面和遠自由面)的兩個模型分別設為示例2和示例3。分別對示例2和示例3計算表1中列出的4種不同的滲透系數(shù)情況下的自由面位置及流場分布，示例2和示例3背景區(qū)域的滲透系數(shù)設置為1 m/s。

表1 示例2和示例3模型中異常體不同滲透率的溢出點高程

與示例1相同。示例2近自由面滲透異常體位置及異常體不同滲透情況下的自由面計算結(jié)果見圖7，模型流場分布見圖8；示例3遠自由面滲透異常體位置及異常體不同滲透情況下的自由面計算結(jié)果見圖9，模型流場分布見圖10。

圖7 示例2近自由面滲透異常體位置及異常體不同滲透情況下的自由面計算結(jié)果

圖8 示例2異常體不同滲透情況下的模型滲流場

圖9 示例3遠自由面滲透異常體位置及異常體不同滲透情況下的自由面計算結(jié)果

圖10 示例3異常體不同滲透情況下的模型滲流場

由圖7(b)可以看出，當異常體位置靠近自由面時，自由面受異常體影響較大。當異常體滲透系數(shù)高于背景區(qū)域滲透系數(shù)時，自由面向異常體方向彎曲；而異常體滲透系數(shù)低于背景區(qū)域時，自由面向異常體相反方向彎曲。且異常體滲透系數(shù)與背景區(qū)域滲透系數(shù)差值越大，則自由面受影響程度越大。由圖9(b)可以看出，當異常體位置距離自由面較遠時，自由面受異常體影響較小。由表1可知，溢出點位置受異常體的影響，異常體離自由面越近，對其位置的影響越大，異常體位置較遠時，溢出點位置的變化不大。

圖8和10顯示了異常體不同滲透系數(shù)影響下的滲流速度變化。當異常體的滲透系數(shù)由低于背景區(qū)域逐漸增大至高于背景區(qū)域時，滲流速度逐漸增大，流向由被低滲透異常體排斥變?yōu)橄蚋邼B透異常體匯聚。

通過示例2和示例3計算結(jié)果的比較，發(fā)現(xiàn)異常體位置對自由面和滲流速度均有影響。異常體離自由面越近，對自由面影響越大。然而，由于自由面附近的滲流速度較小，滲流速度受異常位置的影響相對較小。異常體離自由面越遠，對自由面影響越小，異常體位置對滲流速度的影響越大。

5.3 具有多異常體矩形壩模型(示例4)

示例4旨在探討高滲透和低滲透多種不同滲透異常體對水頭分布、滲流速度和流體壓力的影響。各異常體的具體位置如圖11(a)所示。背景區(qū)滲透系數(shù)仍設為1 m/s，低滲透系數(shù)和高滲透系數(shù)分別設為0.1和10 m/s，高滲和低滲異常體的滲透系數(shù)均與背景區(qū)域相差10倍。模型的邊界條件和網(wǎng)格與示例1一致。

具有多種不同滲透異常體的自由面計算結(jié)果如圖11(b)所示，水頭等值圖、滲流速度圖及流體壓力等值圖如圖12所示。由圖12(a)、12(b)可看出，高滲透區(qū)域滲流速度增大，水頭等值線變稀疏，低滲透區(qū)域排斥周圍滲流，通過低滲透區(qū)域的滲流速度減小，低滲透區(qū)域水頭等值線變密集。兩個低滲透區(qū)域?qū)λ^和滲流流速的影響程度也不相同，深部的低滲透區(qū)域?qū)χ車^和滲流流速的影響更加明顯。在圖12(c)中，流體壓力隨異常滲透率的不同而發(fā)生變化，流體壓力等值線在高滲透區(qū)域變得稀疏，而在低滲透區(qū)域變得密集。

圖11 示例4各異常體位置及自由面計算結(jié)果

圖12 示例4滲流場各參數(shù)計算結(jié)果

通過研究不同滲透系數(shù)的異常體對滲流場參數(shù)分布產(chǎn)生的具體影響，有利于后續(xù)開展通過滲流場參數(shù)分布來反演異常體位置的研究。

5.4 梯形均質(zhì)壩模型(示例5)

為了驗證本文方法在其他模型上的有效性，示例5采用了直角梯形斷面形狀的土石壩模型。壩底長為7 m，壩寬為2 m，上游庫水高程為5 m，下游尾水高程為1 m。該模型中的網(wǎng)格劃分如圖13(a)所示，求解域的滲透系數(shù)設置為1 m/s。以本文方法計算得出的最終自由面如圖13(b)所示。

圖13 示例5模型網(wǎng)格劃分及最終自由面計算結(jié)果

示例5自由面中本文方法計算結(jié)果與其他方法結(jié)果對比及模型滲流場見圖14。由圖14可看出，利用光滑有限元法計算得到的結(jié)果與其他方法具有很好的一致性；隨著深度的增加，滲流速度大小逐漸增大，且溢出點位置以上和以下的滲流均有流向溢出點的趨勢。

圖14 示例5自由面中本文方法計算結(jié)果和其他方法結(jié)果對比及模型滲流場

6 結(jié) 論

本文采用光滑有限元法對土石壩無壓滲流場進行數(shù)值模擬，通過在經(jīng)典算例模型中和其他方法結(jié)果的對比驗證了文中方法的有效性。研究了不同滲透系數(shù)異常體對水頭分布、滲流速度和流體壓力等滲流場參數(shù)的影響，得到如下結(jié)論：

(1) 通過光滑梯度技術將單元面積分優(yōu)化為沿著邊界進行的線積分，建立了無壓滲流問題的光滑有限元模型，不需要對形函數(shù)進行求導運算，減少了計算量，消除了對網(wǎng)格形狀的依賴性，在處理因自由面穿過而產(chǎn)生的畸形單元時更加靈活。

(2) 對經(jīng)典模型的自由面計算結(jié)果與其他方法得到的自由面結(jié)果基本一致，驗證了本文方法的有效性。異常體模型的數(shù)值模擬結(jié)果表明，自由面和溢出點的位置會受異常體位置和滲透系數(shù)的影響。

(3) 通過研究異常體對水頭分布、流體壓力和滲流速度等滲流場參數(shù)的具體影響，為利用滲流場參數(shù)分布來反演滲透異常體位置的研究奠定基礎。