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方法在复杂渗控结构数值计算中的鲁棒性，也为工程实际中的渗流数值模拟提供了参考。