曾 俊，李 星，凌曉鳴，姚孟迪，盧文由頁

(1.武漢大學(xué)水資源與水電工程科學(xué)國家重點實驗室，武漢430072；2.武漢大學(xué)水工巖石力學(xué)教育部重點實驗室， 武漢430072；3.中國三峽建設(shè)管理有限公司，成都 610000)

0 引 言

在水利水電工程施工過程中，滲流問題是影響工程安全至關(guān)重要的因素之一。有限單元法是模擬滲流問題的一種有效工具，它經(jīng)過多年的發(fā)展，其理論愈加成熟。固定網(wǎng)格有限單元法中變分不等式法以其嚴(yán)謹(jǐn)?shù)睦碚摶A(chǔ)，計算工作量小等優(yōu)點得到了廣泛的運用。Brezis等[1]提出了一種非線性相對較弱的擴展壓力法，在一定程度上克服了最初變分不等式非線性較強的缺點；鄭宏等[2]吸收初流量法的思想，將達西定律延拓至全域，建立起了有自由面穩(wěn)定滲流問題的橢圓形Signorini變分不等式法，在理論上解決了出滲點的奇異性以及網(wǎng)格依賴性；陳益峰等[3]將鄭宏所提出的方法運用于工程實際中，并引入了逐步放寬定義的自適應(yīng)罰Heaviside函數(shù)替代原有罰函數(shù)，提高計算數(shù)值的穩(wěn)定性的同時，確保了計算方法的魯棒性，在諸多實際工程之中有良好的運用效果。

在實際工程措施中，為了對滲流進行有效控制，常常布置有排水孔幕等滲控措施。因此，在有限元計算過程中，能否準(zhǔn)確地反映排水孔的滲流控制效果對工程設(shè)計以及運行都具有重要的意義。最初Fipps等[4]提出以點代孔法，但是此方法難以解決無壓滲流時自由面與排水孔的界面問題。隨著理論的完善，諸多學(xué)者提出了一系列更為成熟的排水孔模擬方法：桿單元法[5]、空氣單元法[7]、子結(jié)構(gòu)法[8]等。在上述方法中，桿單元既有較為明確的物理意義，同時不必在空間結(jié)構(gòu)上反映排水孔，計算簡便，經(jīng)濟實用，在工程計算中得到了廣泛的運用。例如王建等[9]在桿單元法的基礎(chǔ)上，進一步在地下水運動的物理意義方面進行改進，建立了匯線單元法，對于單井與復(fù)雜井群的滲流問題進行了深度剖析。

本文在穩(wěn)定滲流框架下，基于含自適應(yīng)罰Heaviside函數(shù)的橢圓形Signorini變分不等式，提出排水孔模擬的桿單元法。并通過一個簡單模型算例，初步驗證本文方法的準(zhǔn)確性與可靠性。同時，針對實際工程的滲流問題也開展了深度分析，證明該方法在模擬含復(fù)雜排水結(jié)構(gòu)的工程滲流問題中也是可靠的，進一步驗證其準(zhǔn)確性。

1 穩(wěn)定滲流分析理論

1856年，法國水力學(xué)學(xué)者、工程師達西(H.Darcy)對水在均質(zhì)砂槽中的流動作了一系列實驗，建立了如下關(guān)系式：

(1)

假定巖土體為非均質(zhì)各向異性可壓縮材料，利用達西定律，則三維穩(wěn)定滲流的控制方程可表示為：

(2)

式中：h=h(x,y,z)為待求水頭函數(shù)；kij為滲透張量。

在實際穩(wěn)定滲流問題中，針對式(2)求解，需要以下4種邊界條件：

(1)水頭邊界條件：

(3)

(2)流量邊界條件：

(4)

(3)出滲面Signorini型互補邊界條件：

(5)

式中：Γs為潛在出滲邊界。

(4)自由面邊界條件：

qn|Ωw=qn|Ωd=0 (onΓf)

(6)

式中：Γf={(x,y,z)|φ=z}為自由面，即濕區(qū)與干區(qū)的分界面。

2 排水孔幕模擬方法

2.1 以管代孔

在滲流分析中，假定水流是一種連續(xù)流體，根據(jù)其連續(xù)性條件，當(dāng)水流經(jīng)過滲透性不同的2種材料的界面時，會產(chǎn)生水流的折射現(xiàn)象，見圖1。

圖1 滲流折射示意圖Fig.1 Seepage refraction

由水流的連續(xù)性可得如下滲流的折射定律：

(7)

式中：K1、K2為介質(zhì)Ⅰ、Ⅱ的滲透系數(shù)；V1、V2為水流在2介質(zhì)的滲透速度，它們與界面法向夾角分別為θ1和θ2。

圖2 孔底式排水孔 h～K、Q～K曲線Fig.2 Drop dowm type drain hole h～K、Q～K graph

圖3 孔口式排水孔h～K、Q～K曲線Fig.3 Over flow type drain hole h～K、Q～K graph

2.2 桿單元法及其有限元格式

使用桿單元法模擬排水孔時，可直接利用排水孔周邊巖體網(wǎng)格節(jié)點建立排水孔桿單元。設(shè)[Kw]為桿單元滲透矩陣，其表達形式可由下式確定：

(8)

式中：ΔF為桿單元的橫截面積；ΔL0為桿單元長度；K為桿單元滲透系數(shù)，其取值按照圖2、圖3曲線確定。

圖4 桿單元示意圖Fig.4 Rod drainage elements

假定一桿單元(見圖4)，i,j2節(jié)點的總水頭分別為Hi,Hj。除i,j2節(jié)點之外，其余各個節(jié)點的總水頭矩陣為{Ha}。根據(jù)滲流控制的基本方程的有限元格式，在沒有排水孔時，可以建立如下方程：

(9)

引入排水孔之后，先只需將排水孔桿單元的滲透矩陣合計入有限元網(wǎng)格的總滲透矩陣中，在滲流計算中，直接求解整合后的滲控方程即可。其有限元格式轉(zhuǎn)變?yōu)椋?/p>

(10)

式中：λ為判定節(jié)點位置的參數(shù)；Zi為i節(jié)點的位置高程，當(dāng)Hi>Zi時，節(jié)點位于自由面以上，λ=1；當(dāng)Hi

3 算例驗證

以含若干個排水孔的穩(wěn)定滲流問題為例，驗證本文方法的可靠性與準(zhǔn)確性。計算模型見圖5。巖體滲透系數(shù)K0=1×10-7m/s，排水孔滲透系數(shù)K=350K0。模型長100 m，寬60 m，高50 m。為便于全面了解和掌握模型內(nèi)的滲流狀況，選取2個典型剖面進行滲流計算成果的整理分析。

圖5 三維計算模型示意圖Fig.5 3D computing model

各典型剖面的具體位置說明如下：A-A為沿排水孔幕橫剖面；B-B為模型中心縱剖面。

人人都知道“創(chuàng)業(yè)難，守業(yè)更難”，李志勇也不例外。隨著野生菌的銷量越來越可觀，他也發(fā)現(xiàn)了很多問題，其中最顯著的問題是，野生菌生長周期短，保鮮時間更短，可是用冰柜將新鮮野生菌冰凍起來，并不能解決問題，一經(jīng)過冰凍的野生菌，解凍后會發(fā)生性狀改變，一下鍋就都軟爛掉了。

模型的主要邊界條件見圖6，與B-B剖面垂直的2個界面取為定水頭邊界，上游定水頭為49.00 m，下游定水頭為25.00 m，其余的外部邊界面為均取為隔水邊界，同時排水孔于廊道溢出點取為Signorini型出滲邊界。同時，由于本文僅研究排水孔幕的效果，為避免廊道取出滲邊界而對排水孔出滲流量產(chǎn)生影響，廊道面除了與排水孔交界處取Signorini型出滲邊界外，其余均取為隔水邊界。

圖6 三維計算模型邊界條件示意圖Fig.6 the boundary conditions of the 3D computing model

為了探究在不同的排水孔間距下本文方法的適用性，分別取排水孔間距為15、10、6 m 3種工況，驗證上述方法的準(zhǔn)確性與可靠性。

(1)當(dāng)排水孔間距為15 m時，模型共劃分單元131 336個，節(jié)點138 776個，共有排水孔3根。本文所計算的各個排水孔流量，以及同工況下鄧琦[10]等使用基于棄單元網(wǎng)格法的空氣單元法計算各排水孔所對應(yīng)的流量比對結(jié)果見圖7。圖8為A-A橫剖面等水頭線分布，圖9為B-B縱剖面自由面及等水頭線分布，圖10為B-B縱剖面滲透坡降矢量圖。

圖7 孔間距15 m時各孔流量比對結(jié)果Fig.7 Flow ratio comparison of each hole (hole pitch:15 m)

圖8 孔間距15 m時A-A剖面等水頭線分布(單位：m)Fig.8 The water head lines distribution of A-A section (hole pitch:15 m)

圖9 孔間距15 m時B-B剖面自由面及等水頭線分布(單位：m)Fig.9 The free surface and water head lines distribution of B-B section (hole pitch:15 m)

圖10 孔間距15 m時B-B剖面滲透坡降矢量圖Fig.10 The seepage gradient vector of B-B section (hole pitch:15 m)

(2)當(dāng)排水孔間距為10 m時，模型共劃分單元144 808個，節(jié)點153 164個，共有排水孔5根。本文所計算的各個排水孔流量，以及同工況下使用基于棄單元網(wǎng)格法的空氣單元法計算各排水孔所對應(yīng)的流量比對結(jié)果見圖11。圖12為A-A橫剖面等水頭線分布，圖13為B-B縱剖面自由面及等水頭線分布，圖14為B-B縱剖面滲透坡降矢量圖。

圖11 孔間距10 m時各孔流量比對結(jié)果Fig.11 Flow ratio comparison of each hole (hole pitch:10 m)

圖12 孔間距10 m時A-A剖面等水頭線分布(單位：m)Fig.12 The water head lines distribution of A-A section (hole pitch:10 m)

圖13 孔間距10 m時B-B剖面自由面及等水頭線分布(單位：m)Fig.13 The free surface and water head lines distribution of B-B section (hole pitch:10 m)

圖14 孔間距10 m時B-B剖面滲透坡降矢量圖Fig.14 The seepage gradient vector of B-B section (hole pitch:10 m)

(3)當(dāng)排水孔間距為6 m時，模型共劃分單元155 024個，節(jié)點163 556個，共有排水孔共有9根。本文所計算的各個排水孔流量，以及同工況下使用基于棄單元網(wǎng)格法的空氣單元法計算各排水孔所對應(yīng)的流量比對結(jié)果見圖15。圖16為A-A橫剖面等水頭線分布，圖17為B-B縱剖面自由面及等水頭線分布，圖18為B-B縱剖面滲透坡降矢量圖。

圖15 孔間距6 m時各孔流量比對結(jié)果Fig.15 Flow ratio comparison of each hole (hole pitch:6 m)

圖16 孔間距6 m時A-A剖面等水頭線分布(單位：m)Fig.16 The water head lines distribution of A-A section (hole pitch:6 m)

圖17 孔間距6 m時B-B剖面自由面及等水頭線分布(單位：m)Fig.17 The free surface and water head lines distribution of B-B section (hole pitch: 6 m)

圖18 孔間距6 m時B-B剖面滲透坡降矢量圖Fig.18 The seepage gradient vector of B-B section (hole pitch:6 m)

從上述結(jié)果分析可知，無論排水孔間距取為何值，滲流在經(jīng)過排水孔幕位置處自由面均有明顯下降，排水孔周邊巖體形成了明顯的降落漏斗，說明排水孔起到了良好的降壓排水效果。各個工況下滲透坡降最大值出現(xiàn)在排水孔與廊道的交界面處，當(dāng)排水孔間距分別為15、10、6 m時，其對應(yīng)的最大滲透坡降值分別為13、7、6。說明巖體滲流確實通過排水孔幕進入到了排水廊道，排水孔發(fā)揮了其效用，說明本文所提出的模擬方法是有效的。

在同一條件下，使用基于Signorini變分不等式的桿單元法的計算流量結(jié)果與使用空氣單元法進行計算的結(jié)果十分接近。當(dāng)排水孔間距分別為15、10、6 m時，其最大相對偏差分別為0.11%、0.91%、4.70%。相對偏差均較小，僅當(dāng)排水孔間距為6 m時相差略大。說明本文所提出的模擬方法，在各排水孔工況下均能準(zhǔn)確模擬，適用性強、較為可靠。針對在模擬過程中所產(chǎn)生偏差，經(jīng)過深度分析，產(chǎn)生這些偏差的原因有以下幾點。

(1)理論基礎(chǔ)不同。本文使用的Signorini型變分不等式法是屬于固定網(wǎng)格法中的變分不等式法，而棄單元法屬于固定網(wǎng)格法中的直覺化方法。直覺化方法通過自由面迭代來區(qū)分干濕區(qū)，其理論基礎(chǔ)不同。

(2)排水孔模擬方法不同。桿單元法不能在空間構(gòu)造上反映排水孔的實際布置情況，而空氣單元法是對排水孔具體結(jié)構(gòu)的直接模擬，在一定程度上空氣單元法能夠更好地還原排水孔的實際結(jié)構(gòu)，其模擬結(jié)果更為精確。

(3)模型網(wǎng)格劃分的差異。模型的網(wǎng)格劃分往往與數(shù)值模擬計算結(jié)果有著密切聯(lián)系，而在實際建模中，很難完全保證各模型網(wǎng)格劃分與網(wǎng)格形態(tài)完全一致，故網(wǎng)格形態(tài)不一致會導(dǎo)致結(jié)果存在有一定的差異。

4 工程應(yīng)用

小灣水電站位于云南省大理州南澗縣和臨滄市鳳慶縣的瀾滄江中游，是瀾滄江中下游河段8個梯級電站的第2級。工程完成蓄水后水位高程為1 215.40 m，達到正常蓄水位高程為1 240.00 m。為了驗證本文方法解決實際工程中含排水孔滲流問題的實用性與可靠性，利用已有的工程資料建立小灣22號壩段的精細(xì)有限元模型，并使用本文所提出方法對該模型中132根排水孔進行模擬。模型見圖19。

圖19 小灣22號壩段有限元模型Fig.19 Finite element model of Xiaowan 22 dam section

模型共劃分單元83 784個，節(jié)點100 839個。模型在計算過程中各個材料滲透參數(shù)取值參考文獻[11]。計算過程中模型各邊界條件如下：各個廊道取為Signorini型出滲邊界，模型左右側(cè)、上下游側(cè)以及底部均取為隔水邊界，拱壩后至二道壩前定水頭邊界取為1 004.00 m，二道壩后定水頭邊界取為991.09 m。同樣，為了驗證本文所提方法在各水頭下的適用性，拱壩上游定水頭根據(jù)工況不同分別取為：1 215.40 m、1 240.00 m。

圖20分別為拱壩上游水位取為1 215.40、1 240.00 m之時，模型的自由面以及等水頭線的分布圖。圖21給出了在各個工況下，流經(jīng)壩體壩基廊道的滲流量的計算模擬值與實測值的對比。

圖20 模型自由面與等水頭線分布(單位：m)Fig.20 The free surface and water head lines distribution of model

圖21 壩基底部廊道流量對比Fig.21 Flow ratio comparison of foundation gallery

對研究結(jié)果進行分析后可知，在數(shù)值計算過程中，小灣拱壩壩內(nèi)排水孔將大部分壩內(nèi)滲流控制在排水孔幕前，且在壩基排水孔附近形成了降落漏斗?？傮w而言，在計算過程中排水孔起到了良好的排水降壓的作用，能夠有效的控制壩區(qū)滲流場。另外，將通過壩基底部廊道的總滲流量與實際工程監(jiān)測值進行對比，其相對偏差較小，相對偏差分別為2.7%、3.3%，準(zhǔn)確性較高。通過上述分析可印證本文所提出方法在工程實例中可得到良好運用，進一步證明方法的準(zhǔn)確性與可靠性。

5 結(jié) 論

本文通過將含自適應(yīng)罰Heaviside函數(shù)的Signorini變分不等式法與桿單元法結(jié)合，提出了一種物理意義比較明確，計算收斂性好，數(shù)值穩(wěn)定性強，網(wǎng)格依賴性小，建模以及計算工作量小的含復(fù)雜排水結(jié)構(gòu)的滲流分析方法。通過對不同排水孔間距下一個簡單滲流問題以及多種水位條件下小灣工程22號壩段滲流問題的模擬，驗證了本文所提出方法的有效性、準(zhǔn)確性與可靠性，為解決含有復(fù)雜排水結(jié)構(gòu)的巖土體的三維滲流問題提供了一種可行方案。

□