基于逐步剖分法求解渗流自由面的方法与流程

本发明涉及一种基于逐步剖分法求解渗流自由面的方法，其属于坝工建设技术领域。

背景技术：

坝体渗流自由面位置的确定在坝工建设中非常重要，坝体渗流自由面也被称为坝体的生命线，只有计算出自由面的位置，才能得到准确的渗流区范围。采用有限元法计算坝体稳定性，通常以坝体自由面为分割线，自由面以下视为饱和区，自由面以上为非饱和区。

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

技术实现要素：

本发明提供一种基于逐步剖分法求解渗流自由面的方法，基于渗流域能量损失率最小求解稳定渗流场自由面的逐步剖分法，将渗流自由面求解转化为求解渗流域能量损失率最小值。

本发明解决上述技术问题的技术方案如下：

一种基于逐步剖分法求解渗流自由面的方法，

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

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

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

采用有限单元法计算渗流问题的关键是运用变分原理确定渗流微分方程的泛函。因此，根据给定的边界条件由(3)式求解渗流水头函数等价于求渗流域ω内泛函的极值问题：

在(2)式中，分别是x、z方向的渗流速度，分别是x、z方向的水力坡降，与渗透力成正比；相当于单位时间渗透力做的功(即功率)，数值上等于渗流的能量损失率。因此，式(2)的积分i(h)与整个渗流域的能量损失率成正比，后文以i(h)表示能量损失率。

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

[k]{h}＝{f}(1)

式(1)中，[k]为总渗透矩阵，[h]是已知常数项列阵；{f}为节点水头列阵。

渗流自由面是渗流域的上部边界，位于其上的部分实际不存在渗流，所以i(h)应不包含这个区域。基于此将整个渗流域分为两部分：渗流实域，即渗流自由面以下的部分bdaf；和渗流虚域，即渗流自由面以上的部分bdc，如图1所示。在求解渗流自由面时，式(2)的泛函只考虑渗流实域，即由上游汇入边界ab、下游流出边界ef、自由渗出边界de，和渗流自由面bd围成的区域。

依据最小功原理，在已确定上、下游及溢出面边界的渗流场中，当渗流达到稳定状态时，域内的能量损失率最小，即满足δi＝0。

以图1所示的矩形坝afcb为例。设上游水位为h1，下游水位为h2，af为不透水边界，d为溢出点。以a为坐标原点，水平向右为x轴的正方向，垂直向上为z轴正方向建立直角坐标系。计算渗流场的能量损失率时，以下游水平面上一点o为坐标原点，水平向右为正方向，表示自由面节点位置不同时能量损失率的大小；垂直向上为正方向，表示不同自由面节点位置，建立如图1所示的直角坐标系。

当渗流场达到稳定状态时，设b3d为真实的渗流自由面，渗流实域b3dfa内的能量损失率为i3；若其他条件不变，将渗流自由面的3点抬高到1点，自由面为b1d时，渗流实域b1dfa的能量损失率为i1，基于以上分析此时应有i1＞i3。若其他条件不变，将渗流自由面的3点压低到2点，自由面为b2d时，计算出渗流实域b2dfa的能量损失率为i2，此时也应有i2＞i3。因此，在上、下游及溢出边界确定的情况下，稳定渗流场中自由面总是位于使得域内能量损失率最小的位置。

逐步剖分法原理：由渗流溢出点向上游入渗点逐层单元逆推，每逆推一层单元，将假定自由面以上单元丢弃，依据实域能量损失率最小确定待求渗流自由面节点的位置，直至逆推到入渗点位置，将每层求得的自由面节点连成光滑的曲线即为真实的自由面曲线。下面以图2a所示的均质矩形坝为例，介绍逐步剖分法求解渗流自由面的具体步骤。

首先，依据横向能量损失率最小计算出渗流溢出点位置，如图2b中节点1所示。

然后，以节点1(即溢出点)为起点，向上游方向逆推一层单元，设节点1左上方的相邻的节点2为此层单元的待求自由面节点。假定线段1-2为局部渗流自由面，赋予其自由面的边界条件，并将其以上单元丢弃，形成新的渗流计算模型，如图2b所示。

通过依次调整假定自由面节点2的高程，基于能量损失率最小求解自由面点原理计算出实域能量损失率最小所对应的节点位置，即为该层单元真实的自由面点。如图3所示，当自由面节点分别位于2a、2b、2c位置时，实域能量损失率分别为i2a(h)、i2b(h)、i2c(h)。因为图3中节点位于2a处的实域能量损失率最小，所以2a即为所求自由面位置。

依据上述步骤，依次逆推求解出各层单元渗流自由面节点的位置。将求解出的每层自由面节点、溢出点以及入渗点连成一条光滑曲线，即为真实的渗流自由面。

逐步剖分法的网格划分要求，选取六节点三角形单元计算二维稳定渗流问题，对渗流域进行横向与竖向的规则剖分，所有网格均划分为四边形；然后将每个四边形单元沿对角线划分为两个三角形单元。

为方便计算过程中自由面节点的调整，避免自由面节点出现跨单元调整造成单元畸形，竖向单元划分的尺寸取

式(4)中，h1为上游水位，h3为溢出点的高程，l为上游汇入点与溢出点之间的水平距离，l0为单元横向尺寸。

如图2所示，横向网格剖分时，距离溢出点1m的初始段内，应至少划分两层单元网格。因为以溢出点为起始点进行逐步剖分推进计算，在初始段，计算实域能量损失时渗流虚域的面积存在较大误差，对渗流自由面的计算结果影响较大，当此段内横向网格仅划分一个单元时，渗流虚域的面积误差对渗流自由面的结果集中在此段内，计算误差过大，也会增加后续计算中自由面节点的误差积累；当该段横向网格划分加密，在推进求解时，渗流虚域在数学模型中所占比例的逐渐减少会消减一部分误差，该段内自由面节点的计算误差会减小，后续计算中的积累误差也会减小。因此加密初始段横向单元划分，可提高整体自由面的计算精度。

根据上述原理，采用fortran语言编制了求解平面稳定渗流自由面的有限元程序，计算步骤如下：

①选取六节点三角形单元对计算模型进行有限元划分，设在x轴方向共划分了m层单元。依据横向能量损失率最小确定溢出点位置；

②由节点1(即溢出点)位置往上游逆推一层网格，如图2b所示，假定自由面节点的位置为高于溢出点位置的节点2，连接1-2并丢弃其上的单元：引入单元函数zb(i)以区分1-2自由面以上或以下的单元，zb(i)＝0表示以上的单元，将此类单元丢弃；zb(i)＝1表示以下的单元，保留此类单元；

③将已知的上、下游和溢出面各节点水头值以及假定的自由面1-2的水头值代入到式(3)中，求解计算域内所有节点的水头值。引入单元分类函数pb(i)，若单元完全位于渗流虚区，即单元内三个角节点的水头值与高程的差值均小于零，则pb(i)＝0；若单元完全位于渗流实区，即单元内三个角节点的水头值与高程的差值均大于零，则pb(i)＝1；

④依据式(2)，计算渗流实域的能量损失i(h)，即pb(i)＝1单元与跨自由面单元的能量损失总和；

⑤将假定的自由面节点高程调整为zq＝zq-1+nα，其中q为逆推的次数，n为第q次逆推计算的次数(n＝1,2,3…)，α为每次计算抬高的高度，重复步骤②、③、④，计算出相应的渗流实域能量损失in(h)。若in(h)＜in-1(h)，则将假定的自由面节点抬高α，重复计算步骤⑤，直至in(h)＞in-1(h)；

⑥由步骤⑤确定的自由面节点往上游方向逆推一层网格，如图2c所示，重复步骤③、④、⑤，直至确定出全部自由面节点的位置。

将上游入渗点、溢出点以及确定的各自由面节点连成一条光滑曲线，即为所求的渗流自由面曲线。其程序流程图如图4所示。

本发明的有益效果是：提出了一种基于渗流域能量损失率最小求解稳定渗流场自由面的逐步剖分法，将渗流自由面求解问题转化为求解渗流域能量损失率最小值的问题。对有确定上、下游边界和自由渗出边界的渗流域进行有限元剖分，求出溢出点位置，然后由溢出边界逐层单元向汇入边界推进，每推进一层单元，将该层的渗流虚域切掉，并基于能量损失率最小求解该层单元的自由面点，直到形成完整的渗流自由面。本方法具有很高的计算精度，避免了现有自由面曲线求解方法中需要迭代计算的缺陷。

附图说明

图1为不同自由面位置的实域能量损失变化示意图；

图2a为原计算模型单元划分；

图2b为求解局部自由面1-2；

图2c为求解局部自由面6-7；

图3为不同自由面点高度的实域能量损失率变化趋势示意图；

图4为有限元计算程序框图；

图5为实施例一的网格划分；

图6为实施例二的网格划分；

图7为实施例一和二的自由面位置；

图8为实施例三的网格划分；

图9为实施例三的自由面位置；

图10为梯形坝的网格划分。

具体实施方式

以下结合附图对本发明的原理和特征进行描述，所举实例只用于解释本发明，并非用于限定本发明的范围。

一种基于逐步剖分法求解渗流自由面的方法，采用fortran语言编制求解平面稳定渗流自由面的有限元程序，其计算步骤如下：

①选取六节点三角形单元对计算模型进行有限元划分，设在x轴方向共划分了m层单元，依据横向能量损失率最小确定溢出点位置；

所述有限元划分采用对渗流域进行横向与竖向的规则剖分，所有网格均划分为四边形，然后将每个四边形单元沿对角线划分为两个三角形单元；

②由溢出点即节点1位置往上游逆推一层网格，假定自由面节点的位置为高于节点1位置的节点2，连接节点1和节点2得到自由面1-2，并丢弃其上的单元：引入单元函数zb(i)以区分节点1-2面以上或以下的单元，zb(i)＝0表示以上的单元，将此类单元丢弃，zb(i)＝1表示以下的单元，保留此类单元；

③将已知的上、下游和溢出面各节点水头值以及假定的自由面1-2的水头值代入到由线性微分方程组汇总的矩阵形式

[k]{h}＝{f}(1)

中，求解计算域内所有节点的水头值；在式(1)中，[k]为总渗透矩阵、[h]是已知常数项列阵、{f}为节点水头列阵；

引入单元分类函数pb(i)，若单元完全位于渗流虚区，即单元内三个角节点的水头值与高程的差值均小于零，则pb(i)＝0；若单元完全位于渗流实区，即单元内三个角节点的水头值与高程的差值均大于零，则pb(i)＝1；

④依据式

计算渗流实域的能量损失i(h)，即pb(i)＝1单元与跨自由面单元的能量损失总和；

在式(2)中，分别是x、z方向的渗流速度，分别是x、z方向的水力坡降，与渗透力成正比；相当于单位时间渗透力做的功，数值上等于渗流的能量损失率；

⑤将假定的自由面节点高程调整为zq＝zq-1+nα，其中q为逆推的次数，n为第q次逆推计算的次数(n＝1,2,3…)，α为每次计算抬高的高度，重复步骤②、③、④，计算出相应的渗流实域能量损失in(h)；

若in(h)＜in-1(h)，则将假定的自由面节点抬高α，重复计算步骤⑤，直至in(h)＞in-1(h)；

⑥由步骤⑤确定的自由面节点往上游方向逆推一层网格，重复步骤③、④、⑤，直至确定出全部自由面节点的位置；

⑦将上游入渗点、溢出点以及确定的各自由面节点连成一条光滑曲线，即为所求的渗流自由面曲线。

具体实施例：

实施例一有电模拟试验解的矩形坝

尺寸为10m×10m的均质土坝，上、下游水位分别为10m和2m，无垂向补给，其底部为不透水边界，下游渗出面边界取为h＝z，由横向能量损失率最小求得溢出点高度为4.75m。建立两种有限元计算模型：①将矩形坝划分为800个两直角边为0.5m的六节点三角形单元，共1681个节点，如图5所示；②将矩形坝划分为220个六节点三角形单元，在x方向0m-9m段三角形单元两直角边为1m，9m到10m段三角形单元两直角边为0.5m，共有273个节点，如图6所示。然后利用逐步剖分法分别求解两种模型，本文方法与电模拟试验解计算结果的对比见图7，与节点虚流量法、初流量法、电模拟试验解的精度对比见表1。

表1有电模拟试验解的矩形坝自由面计算结果比较(单位：m)

由表1可看出，本文方法计算结果与电模拟试验解对比，最大绝对误差为0.27m，最大相对误差为4.94％，比节点虚流量法、初流量法的计算精度高。本文方法通过依次逆推求解渗流自由面，避免了以往有限元法计算渗流自由面时，需不断迭代且易不收敛的缺点。

实施例二有甘油模型试验解的矩形坝

均质土坝尺寸6m×4m，，上、下游水位分别为6m和1m，无垂向补给，其底部为不透水边界，下游渗出面边界取为h＝z，由横向能量损失率最小求得溢出点高度为3.2m。将矩形坝划分为100个六节点三角形单元，在x方向0m-3m段三角形单元两直角边为1m，3m到4m段三角形单元两直角边为0.5m，共有231个节点，如图8所示。本文方法与甘油模型试验解计算结果的对比见图9，精度对比见表2。

表2有甘油模型试验解的矩形坝自由面计算结果比较(单位：m)

由表2可知，本文方法计算结果与甘油模型试验解对比，最大绝对误差为0.08m，最大相对误差为1.83％，本文提出的方法具有很高的计算精度。

实施例三有解析解的梯形坝

图10所示为不透水地基上均质梯形土石坝，上游坡面边坡系数m1＝3，下游坡面边坡系数m2＝2。上游水位h1＝15m，下游水位h2＝4m。坝体高度hc为17m，坝顶fe长度为12m，梯形坝底bc总长为97m，由横向能量损失率最小求得溢出点高度为6.3m。梯形坝体竖向三角形有限单元的划分尺寸为1m，由b到c三段(bg段、gh段、hc段)的横向网格尺寸分别为3m、2m、2m。利用逐步剖分法求解梯形坝上、下游水位之间的渗流自由面位置，与解析解计算结果的对比见表3。

表3有解析解的梯形坝自由面计算结果对比(单位：m)

由表3可知，本文方法计算结果与解析解对比，最大绝对误差为0.37m，最大相对误差为2.85％，计算精度很高，且避免了以往有限元法计算渗流自由面时，需不断迭代且不易收敛的缺点。

综上所述，本文提出了一种依据实域能量损失率最小求解渗流自由面的逐步剖分法，利用该方法求解了有电模拟试验解的矩形坝、有甘油模型试验解的矩形坝、有解析解的梯形坝的渗流自由面，并与试验解和解析解作了对比分析，得到如下结论：

⑴从泛函的物理意义出发，将稳定渗流场自由面的求解问题转化为：在已知上、下游边界条件和自由渗出边界条件下，满足渗流实域能量损失率最小条件的渗流自由面位置；

⑵逐步剖分法对有确定上、下游边界以及自由渗出边界的渗流场，由渗流溢出边界逐段向渗流汇入边界推进，每推进一层单元，依据实域能量损失率最小，求解该层的自由面点，并切掉该层的渗流虚域，直到得到完整的渗流自由面；

⑶该方法具有很高的计算精度，计算结果稳定。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。