基于能量损失率最小原理求解稳定渗流场的自由面曲线

王凤梅，侯兴民，郑珊珊

(烟台大学土木工程学院，山东烟台264005)

0 引 言

坝体渗流自由面又被称为坝体的生命线，只有计算出自由面的位置，才能得到准确的渗流域范围。采用有限元法计算坝体稳定性，通常以坝体自由面为分割线。自由面以下视为饱和区，自由面以上为非饱和区[1]。用有限元法计算渗流场时，自由面的求解是关键。

在稳定渗流求解中，自由面应同时满足水头边界条件和流量边界条件，而自由面位置事先是未知的。国内外学者在解决该问题时，多根据水头边界条件或流量边界条件逐步迭代逼近计算渗流自由面[2-5]，如自由面适应网格法[6]是同时逼近水头和流量边界条件的迭代计算；虚单元法[7- 8]、丢单元法[9]则是通过逼近水头边界条件的迭代计算；初流量法[10]、改进初流量法[11-13]、改进截止负压法[14]是逼近流量边界条件的迭代计算。采用以上方法计算渗流自由面时，为减少误差，使计算结果更加接近实际的自由面，需提高迭代次数，从而增加了计算量。

本文提出了一种基于渗流场能量损失率最小原理确定稳定渗流场自由面的方法。首先对渗流场做有限单元划分，计算渗流逸出点；然后由逸出边界向上游边界逐层推进有限单元，并基于能量损失率最小原理计算每层单元的自由面点，直至上游汇入边界；最后将每一层的自由面点以及逸出点连线得到完整的自由面曲线。采用该方法计算了有电模拟试验解的矩形坝、有甘油模型试验解的矩形坝、有解析解的梯形坝，并与试验结果或解析解作了对比，结果表明该方法具有很高的计算精度。

1 渗流域能量损失率最小原理

1.1 渗流域能量损失率

根据达西定律和地下水运动的连续性条件，不考虑土体和水体的压缩性，均质各向异性土体的二维稳定渗流满足以下控制微分方程：

∂∂xkx∂h∂x+∂∂zkz∂h∂z=0

(1)

式中，h为水头函数；kx、kz分别为x方向、z方向的渗透系数。

对稳定渗流场，边界条件有：①上、下游的水头边界条件；②自由渗出段的水头边界条件和流量边界条件；③底部不透水边界的流量边界条件；④渗流自由面的水头边界条件和流量边界条件。其中，②、④边界是未知的。

根据给定的边界条件，由式(1)求解渗流水头函数等价于求渗流域Ω内泛函的极值问题[1]：

I(h)=∬Ω12kx∂h∂x2+12kz∂h∂z2dxdz

(2)

式中，kx∂h∂x、kz∂h∂z分别为x、z方向的渗流速度；∂h∂x、∂h∂z分别为x、z方向的水力坡降，与渗透力成正比；γwkx∂h∂x2+kz∂h∂z2相当于单位时间渗透力做的功，数值上等于渗流域的能量损失率。因此，式(2)的积分I(h)与整个渗流域的能量损失率成正比。

对式(2)求极小值，将所得线性微分方程组汇总成如下的矩阵形式

[K]{h}={f}

(3)

式中，[K]为总渗透矩阵；{h}为节点水头列阵；{f}为已知常数项列阵。式(3)为采用有限单元法计算渗流问题的基础。

渗流计算域包括完全实域单元、自由面穿过单元、完全虚域单元3部分。本文利用侯兴民等[18]对自由面穿过单元处理的方法，计算能量损失率时，考虑去除渗流计算域中的完全虚域单元与自由面穿过单元中的虚域部分所求的能量损失率最小值，与仅去除计算域中的完全虚域单元(后文将仅去除完全虚域单元后的计算域称为渗流全域[18])所算得的能量损失率最小值，两最小值对应的自由面点位置相同，故无需对自由面穿过单元作特殊处理。图1为矩形坝渗流计算中的两组坐标系：①以不透水层和矩形坝左边界交点A为坐标原点，水平向右为x轴的正方向，垂直向上为z轴正方向的直角坐标系；②以下游水平面上任一点O为坐标原点，水平向右为正方向，垂直向上为正方向的直角坐标系。求解渗流自由面时，采用点A为原点的坐标系，设坝上游水位为h1，下游水位为h2，式(2)的泛函只考虑渗流全域，即由底部不透水边界AF、上游汇入边界AB、下游流出边界EF、自由渗出边界DE和渗流自由面BD所围成的区域。计算渗流场的能量损失率时，采用点O为坐标原点的直角坐标系，其中，水平方向表示渗流逸出点(自由面点)位置不同时渗流全域能量损失率I(h)的大小，垂直方向表示渗流逸出点(自由面点)的坐标h。

图1 逸出点和自由面点的全域能量损失率变化示意

1.2 能量损失率最小原理

在已确定上、下游边界的渗流场中，已知上下游以及逸出面的边界条件，当渗流达到稳定状态时，渗流全域的能量损失率最小，此时应满足δI(h)=0[1]。

如图1a所示，当渗流场达到稳定状态时，设点D为真实的渗流逸出点，渗流全域内的能量损失率为I0；若其他条件不变，将渗流逸出点D抬高到D1点，渗流全域的能量损失率为I1，基于以上分析此时应有I1>I0。若其他条件不变，将渗流逸出点D压低到D2点，计算出渗流全域的能量损失率为I2，此时也应有I2>I0。因此，在上、下游边界确定的情况下，稳定渗流场中渗流逸出点总使得全域内能量损失率，也就是式(2)中的泛函I(h)最小。

同样地，在已确定上下游边界和逸出面边界的条件下，当渗流达到稳定状态时，图1b中自由面曲线上的任一点都满足使得全域能量损失率最小。设B1D为真实的渗流自由面，渗流全域B1DFA内的能量损失率为I1；若其他条件不变，将渗流自由面的1点抬高到3点，即自由面为B3D时，渗流全域B3DFA的能量损失率为I3，基于以上分析，此时应有I3>I1。若其他条件不变，将渗流自由面的1点压低到2点，自由面为B2D时，计算出渗流全域B2DFA的能量损失率为I2，此时也应有I2>I1。因此，在上、下游及逸出边界确定的情况下，稳定渗流场中使得渗流全域内能量损失率最小的位置即为所求的自由面点，这样，求解自由面曲线，实质上成为求解I(h)最小值的问题。

2 基于渗流场能量损失率最小原理求解渗流自由面

图3 有限元计算程序流程

2.1 基于能量损失率最小求解渗流自由面的原理

由能量损失率最小原理求解渗流逸出点位置，再由所求出的渗流逸出点向上游入渗点逐层单元逆推，每逆推一层单元，确定该层渗流自由面节点的位置，直至逆推到入渗点位置，将每层求得的自由面节点及逸出点连接成曲线即为所求的自由面曲线。以图2所示的均质矩形坝为例，介绍基于能量损失率最小求解渗流自由面的原理。

首先，由能量损失率最小原理计算出渗流逸出点1，然后设点1左上方相邻的节点2为此层单元的待求自由面节点，假定线段1-2为局部渗流自由面，赋予其自由面的边界条件。不计其上虚域单元，计算全域的渗流能量损失率。依次调整假定自由面节点2的高程，按照1.2所述方法计算出渗流全域能量损失率最小所对应的节点位置，即为该层单元的自由面点。当自由面节点分别位于2、2a、2b位置时，域内能量损失率分别为I2、I2a、I2b。节点位于2处的全域能量损失率最小，所以2点即为该层单元处的自由面位置。依据上述原理，依次求解出各层单元渗流自由面节点的位置。将求解出的每层自由面节点、逸出点以及入渗点连成一条曲线，即为自由面。

图2 剖面2处自由面位置的计算

2.2 计算步骤

根据上述原理，采用Fortran语言编制求解平面稳定渗流场自由面的有限元程序，程序流程如图3所示，具体计算步骤为:

表1 有电模拟试验解的矩形坝自由面计算结果比较

(1)选取六节点三角形单元划分计算模型，假设在x轴方向共划分了m层单元。

(2)假定逸出点初始水头值为z0=h2+α。其中，h2为下游水位；α为最小误差限(α的取值详见本文算例1)。

(3)将已知的上下游和假定的逸出面各点水头值代入式(3)中，求解计算域内所有节点的水头值。引入单元分类函数pb(i)，若单元位于渗流虚区，即单元内三个角节点的水头值与高程的差值均小于零，则pb(i)=0；反之，pb(i)=1。

(4)依据式(2)，计算渗流全域的能量损失率，即pb(i)=1单元的能量损失率总和。

(5)将假定的逸出点高程调整为z0=h2+nα，n为逸出点试算次数(n=1， 2， 3…)。重复步骤(3)～(5)，计算出相应的渗流全域能量损失In(h)。若In(h)In-1(h)，并依据图1a计算渗流逸出点。

(6)由逸出点位置往上游方向逆推一层网格，如图2所示，假定自由面节点位置为高于逸出点位置的节点2，连接1-2，不计其上的有限单元：可引入单元函数zb(i)以区分已求得的自由面(包括待求的前一层自由面)以上及以下的单元，zb(i)=0表示面上单元并丢弃；zb(i)=1表示面下的单元并保留。

(7)将已知的上、下游和所求逸出面各节点水头值以及图2中假定的待求自由面1-2的水头值代入式(3)中，重复步骤(4)，求解计算域的全域能量损失率。

(8)将假定自由面节点高程调整为zq=zq-1+nα。其中，q为逆推单元层数；n为第q次逆推的试算次数。重复步骤(7)，计算出相应渗流全域的能量损失率In(h)。若In(h)In-1(h)；依据图1b确定第q层的自由面点。

(9)由第(8)步确定的自由面节点往上游方向递推一层网格，重复步骤(7)、(8)，直至确定出全部自由面点。将上游入渗点、逸出点以及确定的各自由面节点连成一条曲线，即为所求的渗流自由面曲线。

3 算例

3.1 计算模型1、2(与电模拟试验解对比)

尺寸为10 m×10 m的均质土坝，上、下游水位分别为10 m和2 m，无垂向补给，其底部为不透水边界，建立如图4所示的2种有限元计算模型：①将矩形坝划分为200个两直角边为1 m的六节点三角形单元，共有441个节点；②将矩形坝划分为800个两直角边为0.5 m的六节点三角形单元，共1 681个节点。利用能量损失率最小原理分别求解两种模型的自由面，与电模拟试验解的精度进行对比，结果见表1。

图4 有限元计算模型

求解逸出点及自由面点时，α取0.01 m。由表1可看出，本文方法计算精度比传统数值计算方法精度高。尽管也需要做多次调整自由面节点的位置以获取能量损失率最小值，但是编制的有限元计算程序的计算速度快。不存在其他有限元法计算方法中可能出现的不收敛情况。

3.2 计算模型3(与甘油模型试验解对比)

均质土坝尺寸6 m×4 m，上、下游水位分别为6 m和1 m，无垂向补给，其底部为不透水边界。建立如图5所示的计算模型3，将矩形坝划分为48个两直角边长分别为1 m的六节点三角形单元，共有117个节点。本文方法与甘油模型试验解计算结果的对比见图6，精度对比见表2。

图5 计算模型3网格划分

图6 自由面位置比较

横坐标/m自由面高度/m甘油模型试验解[16]计算模型3绝对误差/m相对误差/%0.006.006.000.000.001.005.635.50-0.13-2.302.005.105.04-0.06-1.173.004.384.18-0.20-4.564.003.253.10-0.15-4.61

由表2可知，本文方法计算结果与甘油模型试验解对比的最大绝对误差为-0.20 m，最大相对误差-4.61%，尽管单元划分较少，该方法仍保持了很高的计算精度。

3.3 计算模型4(与解析解对比)

图7所示为不透水地基上均质梯形土石坝，上游坡面边坡系数m1=3，下游坡面边坡系数m2=2。上游水位h1=15 m，下游水位h2=4 m。坝体高hHC=17 m，坝顶长lHG=12 m，坝底长lBE=97 m。坝体竖向三角形有限单元的划分尺寸为1 m，其中BC、CD、DE段的横向网格尺寸分别为3、2、2 m。基于能量损失率最小原理求解梯形坝上、下游水位之间的渗流自由面位置，与解析计算结果的对比见表3。

图7 梯形坝的网格划分

表3 有解析解的梯形坝自由面计算结果对比

由表3可知，本文方法计算结果与解析解对比，最大绝对误差为-0.42 m，最大相对误差为-3.89%，计算精度很高。

4 结 论

从泛函的物理意义出发，将稳定渗流场自由面的求解问题转化为求解已知上、下游边界条件和自由渗出边界条件下，满足渗流全域能量损失率最小条件的渗流自由面位置。依据全域能量损失率最小原理求出逸出点，基于该原理由逸出边界逐层单元推至上游汇入边界，将所求各层自由面点与逸出点连线得渗流自由面曲线。利用该方法求解了有电模拟试验解的矩形坝、有甘油模型试验解的矩形坝、有解析解的梯形坝的渗流自由面，并与试验解或解析解作了对比分析，该方法具有很高的计算精度，计算结果稳定。

本文只验证了二维稳定渗流场自由面的求解，应进一步论证所提方法在更复杂渗流问题中的适用性，所提出的逸出点求解及渗流自由面求解方法的理论研究以及工程验证值得做进一步研究的课题。