基于总势能极小原理的渗流溢出点位置计算方法

李远东，侯兴民，郑珊珊，孙伟建

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



基于总势能极小原理的渗流溢出点位置计算方法

李远东，侯兴民，郑珊珊，孙伟建

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

在考虑渗流影响的工程中，自由面的计算是渗流计算的主要内容,精确确定溢出点位置是准确计算自由面的关键。基于总势能极小原理的渗流溢出点位置计算方法，是在溢出点位置调整的过程中，通过寻找全域总势能的极小值确定溢出点位置的一种方法。通过算例分析，得出如下结论：基于总势能极小原理的渗流溢出点位置计算方法弥补了先前溢出点计算方法精度低、需迭代计算、不易收敛的缺点，精确、高效地计算出溢出点位置。

溢出点；渗流自由面；总势能极小值

渗流计算在土石坝稳定性[1]、基坑稳定[2]计算中占有重要地位，其中自由面的确定是渗流计算的主要内容之一。以土石坝稳定性计算为例，通常以坝体的自由面为分割线，自由面以下为饱和区，取浮重度计算；自由面以上为非饱和区，取天然重度计算。

确定溢出点位置是计算自由面的关键。在渗流场中，上游最高水位点为自由面的起始点，溢出点为自由面的终止点。上游水位已知，若能准确确定溢出点位置，无论在应用无网格法[3-5]、虚单元法[6-8]、初流量法[9-11]、单元渗透矩阵调整法[12-14]、或其他方法[15-19]求解渗流自由面时，均可更准确地计算自由面位置。

目前，已经提出的确定溢出点的方法主要有沿坡面滑动法和二次曲线相交法。沿坡面滑动法是将溢出点和渗流自由面上的节点一同进行迭代计算，沿坡面滑动调整其位置[17]。然而此方法溢出点需要与自由面一同迭代，计算量大且会存在不收敛的情况；渗流自由面形状一般可用二次曲线来描述，二次曲线相交法通过确定自由面上三点做二次曲线，与坡面的交点即为溢出点[17]。但此方法中对三点的精度要求太高，误差往往比较大。

针对上述方法的不足，学者们已经做了一些改进与研究。吴梦喜等[6]提出了在可能的溢出面范围内，先假设一排较低的点为溢出点。溢出点以下结点为已知节点，以上结点为未知节点，计算区域各节点的势，反复迭代直到收敛为止。比较溢出点上部一点的势与其位势大小，若势小于位势，则假定的溢出点位置合理，否则应把其上一点作为溢出点，重新计算直至满足要求为止的方法。此方法为局部变网格法，需将自由面与溢出点共同迭代，计算量大。朱军等[13]指出，可将溢出点边界作为第二类边界条件，在求出流量后，通过调整自由项，将其转化为第一类边界条件迭代计算。此方法受起始条件的影响较大，且需反复迭代，计算量大。黄蔚等[15]提出在处理溢出段时，根据水头与高程的关系将溢出面上的点定义为固定点与活动点，通过迭代与线性插值最终找到溢出点位置。此方法也是将自由面与溢出点一同迭代，需重复构造总渗透矩阵与自由项。谢国海等[18]在等效渗透系数法中，根据水头与高程的差值，将溢出点与自由面一起迭代，直至满足水头与高程的差的绝对值小于收敛值。此方法将溢出点与自由面一起迭代，收敛判别条件复杂，有时可能不收敛。

本文从全域总势能的概念入手，提出了一种确定溢出点位置的新方法，不进行迭代，也不重新构造总渗透矩阵，依据全域总势能的极小值即可精确确定溢出点位置，并通过与有解析解的算例对比验证了该方法的精度。

1　渗流的基本理论

1.1支配方程和边界条件

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

(1)

其中kx，ky为x，y方向的渗透系数，h为水头函数。

第一类边界条件又称水头边界条件，为边界上给定的位势函数或水头分布函数，可写为：

(2)

在稳定渗流场中，上、下游的已知水头值，就属于此类边界条件。

第二类边界条件又称流量边界条件，为在边界上给出位势函数或水头的法向导数，其表达式为：

(3)

1.2应用有限单元法求解渗流场

在上述边界条件下，依据式(1)求解水头函数是十分复杂的，对边界的要求十分苛刻。故常采用对泛函求极值的方法求解微分方程，构造如下泛函：

(4)

将全域单元划分后，在每个单元中对节点水头求导，得到单元渗透矩阵，组装为总渗透矩阵，得到求解水头的线性方程组：

[K]{h}={f}

(5)

[K]为总渗透矩阵，{f}为自由项，当没有流量流入与流出时，{f}中的元素为0。求解式(5)即可得到渗流场内各节点的水头值。

2　基于全域总势能确定溢出点

2.1全域总势能

在二维渗流场中，以上游坡面底部为坐标原点，水平方向向左为x轴正方向，铅垂向上方向为y轴正方向，建立直角坐标系，对全域划分有限单元，定义每延米上的单元总势能为：

(6)

其中：Ee为每延米上的单元总势能，J/m。ρ为水的密度；g为重力加速度；Ni为单元形函数；n为单元节点数。

每延米上全域总势能为：

E=∑Ee

(7)

在式(6)中，ρ、g均为常数，故单元总势能值的变化受积分式的影响，而积分式在数值上等于单元中所有点水头值的和。故单元总势能受单元中水头总值的影响，全域总势能受全域水头总值的影响。

2.2基于总势能极小原理的求解溢出点方法

图1给定真实溢出点情况下的边界分布情况

由于真实溢出点的位置是未知的，计算时需要对溢出点的位置进行设置，即将h=y的边界条件赋予可能的溢出段。当设置的溢出点低于真实的溢出点时，部分真实溢出段边界将被定义为不透水边界(如图2中FF0所示)，这阻止了部分溢出段边界处的流量外泄，导致此时的全域水头总值比真实情况的全域水头总值高，使得此时全域总势能比真实情况的全域总势能高；当设置的溢出点高于真实的溢出点时，原本位于BF边界之上，应满足h

图2　给定溢出点低于真实的溢出点情况下的边界分布情况

图3给定溢出点高于真实的溢出点情况下的边界分布情况

因此，全域总势能在赋予溢出点位置升高的过程中，会呈现先减小后增大的变化过程，当全域总势能取极小值时，对应的溢出点位置就是真实的溢出点位置(如图4中F点所示)。

图4全域总势能与溢出点位置关系示意图

2.3计算单元的选取

根据上述全域总势能计算方法及溢出点求解方式，进行有限元数值计算时需要注意以下几点：

(1)在计算全域总势能时，需在求出单元的水头函数后，将水头函数对单元面积积分，得到单元的总势能，再对所有单元的总势能求和。为提高计算效率，推荐采用水头函数及单元形函数较为简单的单元类型；

(2) 为溢出段及上、下游水位赋予第一类边界条件时，应尽可能的准确。对于直线边界，边界上水头函数为线性时最佳。

综上所述，本文选取三角形单元作为计算区域划分的基本单元。水头函数的表达式为

h(x,y)=α1+α2x+α3y

(8)

α1，α2，α3为与单元相关的常数。

此类单元可以很好地模拟第一类边界条件，水头函数对面积的积分式也较为简单。

2.4计算步骤

根据全域总势能确定溢出点位置的程序框图如图5所示。

图5溢出点计算程序框图

下面以斜坡溢出面边界为例，具体解释精确计算溢出点位置的程序步骤：

(1) 全域划分有限单元，释放上、下游水位边界条件，依据式(5)，求解各节点水头值，再依据式(6)、式(7)求解在没有赋予溢出点情况下的全域总势能E0。

(2) 将溢出点升高一个节点，依据式(5)、式(6)、式(7)求解全域总势能Ep，p为执行此步骤的序号。

(3) 判断Ep与Ep-1的大小，若Ep>Ep-1，执行下一步；若Ep-1>Ep，回到步骤(2)。

(4) 设此时的p值为m，当p=m-2，m-1，m时，溢出点的坐标分别为(x1，y1)，(x2，y2)，(x3，y3)，令E10=Em-2。

(5) 设最小误差为α，令y2=y1+qα，x2=x1+(x3-x1)qα/(y3-y1)，根据(x2,y2)的新坐标，微调总渗透矩阵[K]，并将(x2，y2)作为新的溢出点，依据式(5)、式(6)、式(7)求解全域总势能E1q，q为执行此步骤的序号。

(6) 判断E1q与E1q-1的大小，若E1q>E1q-1，执行下一步；若E1q-1>E1q，回到步骤(5)。

(7) 取y=y2-α，x=x2-(x3-x1)α/(y3-y1)，(x，y)即为所求的溢出点坐标。

3　算　例

3.1算例1

以10m×10m均质土坝为例进行分析，土坝上下游水位分别为10m和2m，底部为不透水边界，该问题的自由面解析表达式为y2=100-8x。其溢出点位置的解析解为4.472m。采取两种不同的单元剖分进行计算，计算模型1共200个单元，以直角边长为1m的三角形进行划分，如图6(a)；计算模型2共50个单元，以直角边长为2m的三角形进行划分，如图6(b)。应用上文所述方法求解溢出点位置。全域总势能的变化情况如表1所示。

图6　算例1网格剖分图

由表1知，在计算模型1中，求出的溢出点位置为4.55m，与解析解相差0.08m，相对误差为1.74%；在计算模型2中，求出的溢出点位置为4.35m，与解析解相差-0.12m，相对误差为-2.73%，满足规范的误差要求，在大网格划分下，仍然能得到接近实际位置的溢出点。

由表2知，同样采取计算模型1时，本文方法算得的溢出点，精度比节点虚流量法、初流量法、改进的初流量法、改进的截至负压法、改进的丢单元法等方法都高。

表2　溢出点精度对比

3.2算例2

为了进一步验证本文方法，选取带有坡面的土石坝坝体进行溢出点位置的计算，坝体上、下游水位分别为15m和4m，上、下游坡比分别为1∶3和1∶2，最高处高17m，底部横向宽度为97m，底部为不透水边界，解析解求得的溢出点的位置为6.23m。

如图7所示，分别以高为1m的三角形与高为2m的三角形对此坝体进行单元剖分，应用本文方法编制有限元程序，对溢出点位置进行计算，全域总势能如表3所示。

图7　算例2网格剖分图

由表3知，1m网格划分下，求出溢出点位置为6.40m，与解析解的差为0.17m，相对误差为2.72%；在2m网格划分下，求出溢出点位置为6.00m，与解析解的差为-0.23m，相对误差为-3.69%，满足工程精度要求。

4　结　论

本文提出了一种以全域总势能极小值精确计算溢出点的方法，应用有限元程序求解了二维计算模型与土石坝坝体的溢出点位置并与其他方法作了对比，主要结论如下：

(1) 应用本文基于全域总势能计算溢出点的方法，相比以前的一些方法，可更精确地确定溢出点位置。

(2) 本文方法避免了同类研究由于将溢出点与自由面一同迭代，结果易不收敛的问题；在求解溢出点时，不必重新构造总渗透矩阵，减少计算量，提高计算效率。

(3) 在较大单元网格划分中，本文方法仍然具有很高的精度。

[1]刘娟奇，王志强，梁收运.库水位下降对新集水库均质土坝渗流及稳定性影响分析[J]. 水利与建筑工程学报,2014,12(6):38-43.

[2]位俊俊，张利伟，孔德志.基坑渗流稳定分析[J].水利与建筑工程学报,2012,10(3):79-87.

[3]DarbaniM,OuahsineA,VillonP,etal.Meshlessmethodforshallowwaterequationswithfreesurfaceflow[J].AppliedMathematicsandComputation, 2011,217(11):5113-5124.

[4]JieYX,LiuWJ,LiGX.ApplicationofNEMinseepageanalysiswithafreesurface[J].MathematicsandComputersinSimulation, 2013,89(2):23-37.

[5]TangJing,LiuYongbiao.Penaltyfunctionelementfreemethodtosolvecomplexseepagefieldofearthfilldam[J].IERIProcedia, 2012(1):117-123.

[6]吴梦喜,张雪勤.有自由面渗流分析的虚单元法[J].水利学报,1994(8):67-71.

[7]凌道盛.有自由面渗流分析的虚节点法[J].浙江大学学报(工学版),2002,36(3):243-246.

[8]崔皓东,朱岳明.有自由面渗流分析的改进节点虚流量全域迭代法[J].武汉理工大学学报(交通科学与工程版),2009,33(2):238-241.

[9]张有天,陈平,王镭.有自由面渗流分析的初流量法[J].水利学报,1988(8):18-26.

[10]王媛.求解有自由面渗流问题的初流量法的改进[J].水利学报,1998(3):68-73.[11]潘树来,王全凤,俞 缙.利用初流量法分析有自由面渗流问题之改进[J].岩土工程学报,2012,34(2):202-209.[12]王贤能.有自由面渗流分析的高斯点法[J].水文地质工程地质,1997(6):1-4.

[13]朱军,刘光廷.改进的单元渗透矩阵调整法求解无压渗流场[J].水利学报,2001(8):49-52.

[14]黄蔚,刘迎曦,周承芳.三维无压渗流场的有限元算法研究[J].水利学报,2001(6):33-36.

[15]梁业国,熊文林,周创兵.有自由面渗流分析的子单元法[J].水利学报,1997(8):34-38.

[16]王均星,吴雅峰,白呈富.有自由面渗流分析的流行单元法[J].水电能源科学,2003,21(4):23-26.

[17]毛昶熙,段祥宝,李祖贻，等.渗流数值计算与程序应用[M].南京:河海大学出版社,1999.

[18]谢国海,章广成.二维无压稳定渗流分析的改进方法[J].地球与环境,2005,33(S1):52-56.

[19]张宜虎.采用等效渗透系数法搜索渗流自由面[J].工程地质学报,2004,12(2):177-192.

Calculation Method of Seepage Overflow Point Based on the Principle of Minimum Total Potential Energy

LI Yuandong, HOU Xingmin, ZHENG Shanshan, SUN Weijian

(SchoolofCivilEngineering,YantaiUniversity,Yantai,Shandong264005,China)

In engineering that needs to consider the effect of seepage, determining the overflow point precisely is the key to the calculation of free surface, which is the main content of seepage calculation. Based on the principle of minimum total potential energy to calculate the overflow point is a method that the location of overflow point is determined through calculating the minimum of global total potential energy by adjusting its position. Example analysis showed that this method can calculate the location of overflow point precisely and efficiently, which overcomes the shortcomings of previous methods such as low accuracy, requiring iterative calculation and difficult to converge.

overflow point; seepage free surface; minimum total potential energy

10.3969/j.issn.1672-1144.2016.04.017

2016-04-17

2016-05-14

国家自然科学基金项目(51479174)

李远东(1990—)，男，山东青岛人，硕士研究生，研究方向为建筑与土木工程。 E-mail：wddd2013@163.com

侯兴民(1970—)，男，山东寿光人，博士，教授，主要从事岩土工程方面教学与科研工作。 E-mail：houxm@ytu.edu.cn

TV139.14

A

1672—1144(2016)04—0084—05