張 淼 ，鐘小彥 ，李國(guó)棟

（1.中國(guó)電建集團(tuán)西北勘測(cè)設(shè)計(jì)研究院有限公司，陜西 西安 710065；2.北京清溪環(huán)?？萍加邢薰?，北京 100071；3.西安理工大學(xué)，陜西 西安 710048）

0 引言

在壩工工程、邊坡工程以及地下工程中，滲流對(duì)工程結(jié)構(gòu)的穩(wěn)定性有著很大影響，為了降低滲流對(duì)建筑物的危害并控制滲流，各種各樣的滲控結(jié)構(gòu)如排水孔、井、渠、廊道、防滲帷幕等成為工程設(shè)計(jì)不可分割的一部分。隨著近些年有限元數(shù)值模擬計(jì)算的發(fā)展，數(shù)值分析開始成為滲流計(jì)算中的主流計(jì)算方法，但對(duì)于含有滲控結(jié)構(gòu)的滲流問(wèn)題仍然有較大難度，主要是因?yàn)闈B控結(jié)構(gòu)的比尺梯度變化導(dǎo)致有限元網(wǎng)格劃分困難及復(fù)雜結(jié)構(gòu)中的流體非線性問(wèn)題。為克服有限元建模上的困難，人們提出了用于數(shù)值模擬的排水子結(jié)構(gòu)法[1]、半解析法[2]及復(fù)合單元法[3]等，相關(guān)方法均對(duì)計(jì)算的邊界條件進(jìn)行了一定程度的簡(jiǎn)化，使得有限元的真實(shí)模型并未得到真實(shí)反映，且計(jì)算方法繁復(fù)不宜推廣。既然滲流屬于流體力學(xué)的范疇，也就適用于流體力學(xué)的計(jì)算方法，因此本文主要是應(yīng)用流體力學(xué)平衡方程，將VOF法與多孔介質(zhì)結(jié)合用于模擬滲流場(chǎng)，以FLUENT軟件為抓手展開計(jì)算，過(guò)程不斷與傳統(tǒng)經(jīng)典算例的計(jì)算結(jié)果進(jìn)行比對(duì)分析，并通過(guò)調(diào)整滲控結(jié)構(gòu)尺寸以驗(yàn)證該方法的有效性和魯棒性。

1 計(jì)算方法

本文整體的計(jì)算思路是，通過(guò)有限元求解納維-斯托克斯方程(Navier-Stokes equations)和連續(xù)性方程以計(jì)算出整體的壓力和速度場(chǎng)，以VOF方法來(lái)計(jì)算自由面的位置。

1.1 基于多孔介質(zhì)模型的有限體積法計(jì)算滲流場(chǎng)。

有限體積法（Finite Volume Method)是將計(jì)算區(qū)域劃分成控制單元，借助將控制方程對(duì)控制單元進(jìn)行積分來(lái)導(dǎo)出離散方程，利用其變步長(zhǎng)網(wǎng)格的簡(jiǎn)單離散形式，保證所得差分方程仍然是守恒的。在有限體積法的基礎(chǔ)上，借助多孔介質(zhì)模型，將滲流的應(yīng)用物理概念導(dǎo)入離散元中進(jìn)行求解。由于水是不可壓縮流體，在滲流介質(zhì)為各向同性且當(dāng)滲流流態(tài)為層流時(shí)，忽略加速度、水流阻力、質(zhì)量力和慣性損失，則由一般不可壓縮流體的連續(xù)性方程及那維斯托克斯方程可知多孔介質(zhì)模型控制方程為：

式中：v為流體的運(yùn)動(dòng)粘性系數(shù)；μ是微元斷面的平均速度；p為滲透壓力；為粘性阻力系數(shù)。同時(shí)結(jié)合達(dá)西定律可以得出：

式中：J為滲透坡降；K為滲透系數(shù)；S為滲徑；h為水頭；γ為水的容重；p為滲透壓力。最終可以得到：

在使用該模型計(jì)算滲流場(chǎng)時(shí)，需要注意到粘性阻力系數(shù)與滲流場(chǎng)的地質(zhì)條件關(guān)系最為密切，因此滲透系數(shù)K的選取對(duì)于整體計(jì)算至關(guān)重要，選取時(shí)應(yīng)根據(jù)不同材料的物理特性準(zhǔn)確選取。

1.2 基于VOF方法計(jì)算滲流自由面

對(duì)于自由滲流，浸潤(rùn)線的計(jì)算是關(guān)鍵。運(yùn)用有限元的方法計(jì)算滲流通常是假設(shè)浸潤(rùn)線，借助自由面需要滿足的水頭條件逐漸迭代逼近正確結(jié)果[4]。從流體力學(xué)的角度來(lái)看，浸潤(rùn)線就是滲流的自由面，可以運(yùn)用自由面捕捉體積率的方法獲得。作為在固定歐拉網(wǎng)格下的表面追蹤方法，VOF法可以通過(guò)計(jì)算水相和空氣相的比率來(lái)獲得滲流流體的自由面，通過(guò)對(duì)體積率函數(shù)的求解獲得自由界面，再利用插值進(jìn)行幾何重構(gòu)得到界面，它采用分短線近似的方法來(lái)表示自由水面線,這里為滲流流速。

圖1 幾何重構(gòu)自由水面示意圖

2 實(shí)例驗(yàn)證計(jì)算方法的準(zhǔn)確性和穩(wěn)定性

2.1 不同滲透系數(shù)的兩個(gè)矩形壩體

有一個(gè)混凝土壩體，由兩個(gè)滲透系數(shù)不同的矩形區(qū)域組成，下游區(qū)域的滲透系數(shù)是上游的10倍（即k2/k1=10），基于此，數(shù)值震蕩有可能在兩個(gè)滲透系數(shù)不同的材料交界面上出現(xiàn)，一般這種情況下有可能出現(xiàn)計(jì)算精度和計(jì)算有效性的明顯下降。數(shù)值計(jì)算方法計(jì)算的結(jié)果如圖2中的紅色線條所示。

圖2 不同滲透系數(shù)的兩個(gè)矩形壩體自由面計(jì)算結(jié)果

可以看出，由VOF法計(jì)算的自由面的位置和前面文獻(xiàn)中用EP方法計(jì)算的結(jié)果相吻合，證明VOF法在滲透系數(shù)變化劇烈的工況中仍然具有很高的精度。

2.2 尾水的矩形三維壩體實(shí)例

這個(gè)實(shí)例是一個(gè)經(jīng)典算例，是一個(gè)有尾水的矩形三維均質(zhì)壩體的滲流情況。均質(zhì)壩體的尺寸長(zhǎng)寬高分別為9 m×5 m×10 m，上游初始最大水頭為10 m，下游水頭為2 m，如圖3所示的研究實(shí)例幾何形體。計(jì)算過(guò)程中，三個(gè)方向的網(wǎng)格劃分均為0.3 m，滲透系數(shù)K=0.0005 m/day,松弛因子w=1.2，收斂因子e=0.5。剖面A-A上的計(jì)算結(jié)果在圖4中顯示，由圖表中可以看出VOF方法和其他方法得到的結(jié)果很好的吻合。圖5中給出了三維自由面。

圖3 實(shí)例幾何形體自由面

圖4 截面計(jì)算結(jié)果對(duì)比

2.3 非均質(zhì)矩形三維壩體計(jì)算實(shí)例

在這個(gè)例子中，和圖5中的算例有著相同的幾何體，邊界條件和計(jì)算參數(shù)的假設(shè)也都相同，但在該算例中考慮非均勻的介質(zhì)。水力滲透系數(shù)值如下：K1=0∶0005 m/day，K2=0∶005 m/day，K3=0 ∶00005 m/day。如下圖 6所示，同樣用VOF求解出壩體滲流的自由面。可以看出，非均勻介質(zhì)很大程度的影響了自由面，而自由面的計(jì)算結(jié)果也真實(shí)地反映了實(shí)際情況。

圖5 非均質(zhì)壩體滲流自由面計(jì)算結(jié)果

計(jì)算出的穿排水孔中心線平面水頭等值線如下圖8所示。從圖中可以看出，當(dāng)壩段未布置排水孔和帷幕時(shí)，大壩下游斜邊界的出滲情況明顯好于不設(shè)置排水孔和帷幕的工況，且從圖9可以看出，有或無(wú)排水孔幕條件下壩基揚(yáng)壓力的影響也非常大，排水孔幕對(duì)滲流流場(chǎng)的分布影響很大,由于排水孔幕的存在導(dǎo)致建基面揚(yáng)壓力和滲流自由面急劇降低。因此也可以看出，排水孔、帷幕等滲控結(jié)構(gòu)對(duì)大壩的滲流情況起到控制性影響，在滲流計(jì)算中，對(duì)排水孔的模擬非常重要。

同時(shí)圖7-(a)與圖8的計(jì)算結(jié)果與文獻(xiàn)[5]幾乎完全相同，也驗(yàn)證了VOF法的準(zhǔn)確性。

3 在復(fù)雜滲控結(jié)構(gòu)的應(yīng)用及敏感性分析

3.1 工程概況

本文主要是為了驗(yàn)證VOF法計(jì)算滲流的穩(wěn)定性及魯棒性，故為方便進(jìn)行結(jié)果對(duì)比，采用文獻(xiàn)[5]中的典型算例進(jìn)行，即采用的壩體體型參數(shù)、邊界條件及計(jì)算條件也均與文獻(xiàn)條件相同，整體工況如下圖7所示。一混凝土重力壩典型壩段，壩高和壩長(zhǎng)分別為170 m和30 m，頂部和底部寬度分別為12 m和123 m，大壩上、下游水位分別為168.0 m和28.5 m。在計(jì)算空間上，向壩踵和壩趾上、下游各取200 m，從建基面向下深度取170 m。在壩上游面考慮0.8 m的混凝土防滲層，如圖所示共布置7條截面尺寸為2 m×2 m的水平排水廊道。壩基布置一排間距5 m，孔徑12.73 cm的垂直排水孔，排水孔深入壩基40 m，與廊道相通進(jìn)行系統(tǒng)性排水。壩基設(shè)有深60 m，寬2.5 m的防滲帷幕。滲透系數(shù)分別為：壩體混凝土k=4.09×10-6m/d；混凝土防滲層 k=1.62×10-6m/d；壩基巖體 k=2.30×10-2m/d；防滲帷幕k=7.18×10-3m/d。

圖6 含排水孔幕和排水廊道的重力壩段示意圖

本算例計(jì)算一排排水孔的工況，因此應(yīng)用對(duì)稱邊界到壩身及壩基的前、后兩側(cè)；壩基左、右與底邊界均設(shè)為不透水條件；壩基排水孔視為出滲邊界，水頭值設(shè)為與其相通的排水廊道底高程；壩體排水廊道和垂直排水孔均視為出滲邊界。

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

圖7 順河向排水孔剖面的水頭等值線和自由面

圖8 有或無(wú)排水孔幕條件下壩基揚(yáng)壓力對(duì)比

4 驗(yàn)證復(fù)雜滲控結(jié)構(gòu)中VOF法的敏感性分析

為了進(jìn)一步驗(yàn)證計(jì)算方法的敏感性，本研究進(jìn)一步對(duì)排水孔的尺寸及其布置形式的變化進(jìn)行敏感性分析，計(jì)算結(jié)果見圖9。

4.1 排水孔孔距對(duì)計(jì)算結(jié)果的影響

采用控制變量法，先固定排水孔直徑不變的情況下，改變孔間距s分別為2 m、3 m、4 m和5 m，計(jì)算得到大壩滲流自由面位置如圖10-a所示?？梢钥闯?，當(dāng)排水孔間距小于4 m時(shí)，自由面位置隨著間距等額減小明顯降低，排水孔間距對(duì)滲流場(chǎng)的影響趨勢(shì)增強(qiáng)；然而當(dāng)排水孔間距大于4 m時(shí)，間距等額增大對(duì)滲流場(chǎng)的影響趨勢(shì)減弱。然而，排水孔間距對(duì)單寬滲流量的影響卻很小，究其原因主要是因?yàn)榘殡S著排水孔間距等額減小，排水孔數(shù)量不斷增加，排水量按常理也應(yīng)該增加，但考慮到整個(gè)滲控結(jié)構(gòu)周邊的滲透壓逐漸變小，導(dǎo)致排水孔排水的強(qiáng)度減小了，兩者之間有一個(gè)抵消作用，因此排水孔間距對(duì)滲流量的影響較小。

圖9 不同排水孔間距&孔徑工況對(duì)應(yīng)的滲流自由面分布

4.2 排水孔孔徑對(duì)計(jì)算結(jié)果的影響分析

如果固定各排水孔之間距離為2 m不變,孔徑d分別取10 cm、15 cm、20 cm、25 cm時(shí)，可計(jì)算出滲流自由面位置如上圖10-b所示。伴隨著孔徑增大，單寬排水量增大，滲流自由面隨之逐漸降低,但單寬排水量和自由面降低的幅度均不大?？梢缘贸觯?dāng)上下游水位差固定時(shí)，排水孔間距比孔徑對(duì)壩基揚(yáng)壓力的分布影響更大。如果增大排水孔間距到5 m，由從圖10-c可以看出，設(shè)置排水孔的排水效果已經(jīng)不太明顯了，因此在實(shí)際工程應(yīng)用中，排水孔的間距不宜過(guò)大，結(jié)果見表1。

表1 不同排水孔間距&孔徑工況對(duì)應(yīng)的單寬滲流量

5 結(jié)論

滲流計(jì)算本身就是流體力學(xué)計(jì)算的分支，從理論上講，利用流體數(shù)值計(jì)算方法可以計(jì)算線性和非線性滲流，具有更廣泛適用性。本文結(jié)合有限體積法，多孔介質(zhì)模型、VOF模型、以及有限體積自由面捕捉技術(shù)，進(jìn)一步研究了基于多孔介質(zhì)模型和自由面捕捉技術(shù)的滲流計(jì)算方法的有效性及敏感性。通過(guò)對(duì)多個(gè)算例的結(jié)果對(duì)比表明，VOF自由面捕捉技術(shù)能夠準(zhǔn)確地對(duì)滲流自由面及出滲點(diǎn)進(jìn)行捕捉，同時(shí)通過(guò)復(fù)雜滲控結(jié)構(gòu)中排水孔孔距及孔徑的對(duì)比計(jì)算，驗(yàn)證了VOF方法在復(fù)雜滲控結(jié)構(gòu)數(shù)值計(jì)算中的魯棒性，也為工程實(shí)際中的滲流數(shù)值模擬提供了參考。