VOF數值模擬在含復雜滲控結構滲流問題中的應用

張 淼 ，鐘小彥 ，李國棟

（1.中國電建集團西北勘測設計研究院有限公司，陜西 西安 710065；2.北京清溪環保科技有限公司，北京 100071；3.西安理工大學，陜西 西安 710048）

0 引言

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

1 計算方法

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

1.1 基于多孔介質模型的有限體積法計算滲流場。

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

式中：v為流體的運動粘性系數；μ是微元斷面的平均速度；p為滲透壓力；為粘性阻力系數。同時結合達西定律可以得出：

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

在使用該模型計算滲流場時，需要注意到粘性阻力系數與滲流場的地質條件關系最為密切，因此滲透系數K的選取對于整體計算至關重要，選取時應根據不同材料的物理特性準確選取。

1.2 基于VOF方法計算滲流自由面

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

圖1 幾何重構自由水面示意圖

2 實例驗證計算方法的準確性和穩定性

2.1 不同滲透系數的兩個矩形壩體

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

圖2 不同滲透系數的兩個矩形壩體自由面計算結果

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

2.2 尾水的矩形三維壩體實例

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

圖3 實例幾何形體自由面

圖4 截面計算結果對比

2.3 非均質矩形三維壩體計算實例

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

圖5 非均質壩體滲流自由面計算結果

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

同時圖7-(a)與圖8的計算結果與文獻[5]幾乎完全相同，也驗證了VOF法的準確性。

3 在復雜滲控結構的應用及敏感性分析

3.1 工程概況

本文主要是為了驗證VOF法計算滲流的穩定性及魯棒性，故為方便進行結果對比，采用文獻[5]中的典型算例進行，即采用的壩體體型參數、邊界條件及計算條件也均與文獻條件相同，整體工況如下圖7所示。一混凝土重力壩典型壩段，壩高和壩長分別為170 m和30 m，頂部和底部寬度分別為12 m和123 m，大壩上、下游水位分別為168.0 m和28.5 m。在計算空間上，向壩踵和壩趾上、下游各取200 m，從建基面向下深度取170 m。在壩上游面考慮0.8 m的混凝土防滲層，如圖所示共布置7條截面尺寸為2 m×2 m的水平排水廊道。壩基布置一排間距5 m，孔徑12.73 cm的垂直排水孔，排水孔深入壩基40 m，與廊道相通進行系統性排水。壩基設有深60 m，寬2.5 m的防滲帷幕。滲透系數分別為：壩體混凝土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 含排水孔幕和排水廊道的重力壩段示意圖

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

3.2 計算結果分析

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

圖8 有或無排水孔幕條件下壩基揚壓力對比

4 驗證復雜滲控結構中VOF法的敏感性分析

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

4.1 排水孔孔距對計算結果的影響

采用控制變量法，先固定排水孔直徑不變的情況下，改變孔間距s分別為2 m、3 m、4 m和5 m，計算得到大壩滲流自由面位置如圖10-a所示。可以看出，當排水孔間距小于4 m時，自由面位置隨著間距等額減小明顯降低，排水孔間距對滲流場的影響趨勢增強；然而當排水孔間距大于4 m時，間距等額增大對滲流場的影響趨勢減弱。然而，排水孔間距對單寬滲流量的影響卻很小，究其原因主要是因為伴隨著排水孔間距等額減小，排水孔數量不斷增加，排水量按常理也應該增加，但考慮到整個滲控結構周邊的滲透壓逐漸變小，導致排水孔排水的強度減小了，兩者之間有一個抵消作用，因此排水孔間距對滲流量的影響較小。

圖9 不同排水孔間距&孔徑工況對應的滲流自由面分布

4.2 排水孔孔徑對計算結果的影響分析

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

表1 不同排水孔間距&孔徑工況對應的單寬滲流量

5 結論

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