专利名称：一种基于等势面的地下水渗流量计算方法
技术领域：
本发明属于水利工程地下水控制技术领域，涉及地下水渗流工程数值模拟技术， 尤其是涉及一种利用有限元求解地下水渗流量的计算方法。
背景技术：
渗流是一种地下水运动，当建筑物及其地基内存在水位差时，地下水将从水位高 的地方向水位低的地方产生流动，形成渗流。渗流是一门与水力学和岩土力学有着密切关 系的学科，随着近代科学技术的不断发展，渗流在基本理论、试验手段、计算方法和应用等 方面都得到了极大发展，逐渐成为一门专门的学科，已能解决各种复杂的工程问题。各种复 杂情况下的土石坝渗流，都可以在计算机上模拟出来。实践证明，利用有限元法求解均质或 非均质、各向异性或各向同性以及复杂边界条件的土石坝渗流问题，虽然得到的解是近似 的，然而却是满意的解答，对于土石坝渗流问题，已基本上可取代模拟试验。渗流计算中最重要的内容之一就是分析和预测渗流量，渗流量计算的准确性和精 度对于渗流问题的分析往往是至关重要的，特别是对大坝防渗系统的布置、防渗效果评价、 坝后排水设施的分布与设置以及排水井的排水效果等具有重要的参考价值。例如，在土石 坝病险的诊断中，需要根据渗流量、渗透水透明度和水质观测资料及巡查结果，结合渗透比 降，判别有无管涌、流土、接触冲刷、接触流失等渗透变形现象。对渗流量计算目前主要采用中断面法和等效结点流量法。中断面法的主要原理是，在求得渗流场水头函数H的有限单元法数值解后，对于 任意单元(如8节点等参单元)，选择其一对面之间4条棱的中点构成的截面(中断面)作为 过流断面；在二维问题中是以中线作为过流断面的。有限单元法求得的水头函数数值解精 度较高，一般能满足工程应用的要求。中断面法计算原理简单且很容易通过程序实现。但 是，由于水头函数数值解为数值离散解，并且，实际选用的过流断面为各个单元的中断面， 因此当计算区域材料分区和地质条件复杂时，单元形状很不规则，其中断面也是极不规则 的扭曲面，所计算出的渗流量的准确性大大降低，有时并不能满足工程应用的要求。如图1 所示，图中弧形实线为自由面4，虚线表示中面法的三个过水断面剖面1、2、3，因为剖分网 格的问题，只有剖面3是正常的有效过水断面，剖面1和2都没有覆盖整个渗流区域，因而 不可能计算准确。等效结点流量法将任一过流断面上的渗流量表示成相关单元的传导系数与相应 结点水头的乘积的代数和。它避免了对水头离散解的进一步求导运算，所求得的渗流量计 算精度与水头解的计算精度同阶。但是缺点是渗流量缺乏明确的物理意义，过水断面的选 择不太容易，所得到的断面渗流量为前后单元所谓等效结点流量的代数和，计算时需要将 前后单元的贡献分开计算，否则计算结果为0。而将前后单元的贡献分开计算需要给定条件 或人工干涉，工作量大且容易出错。如图2中所示，图中实线表示自由面4，虚线所包围的结 点构成该过水断面，在计算该断面流量时，如果简单的将所有相关单元的传导系数与相应 结点水头的相乘，则代数和为0，因此必须将左右两侧的单元分开单独计算才可以。
有学者在渗流有限元数值计算的基础上，研究计算任意断面渗流量的方法，以解 决复杂渗流场中渗流量的计算问题(相关文献祁书文，《基于有限元法的复杂三维渗流场 渗流量计算方法研究》，河海大学硕士论文，2007. 5)。如图3所示，该方法在三维渗流场中取相互靠近的两个平面5、6并进行剖分，求出 两个面5、6相邻单元形心处的水头，然后由两个面之间的距离为L，计算出通过该流管的渗 透流量。该方法需要假定两平面形心距离L，需要将计算断面重新剖分，需要用到迭代法求 解截面单元的形心在整体区域内所经过单元的单元坐标系中的局部坐标，计算时间较长， 从理论上来说假定的平面形心距离L也影响流量的计算精度。
发明内容
本发明的目的是提供一种计算渗流区域渗流量的便利计算方法，从而提高计算精 度，简化用户操作过程，避免了复杂渗流场中某一流量断面既有正向流动又有反向流动时 所带来的流量计算误差。本发明的技术方案是
一种基于等势面的地下水渗流量计算方法， 计算步骤如下
(1. 1)求解渗流场。用常规的有限元法求解渗流场，求出各单元结点水头值；
(1.2)确定等势面的水头值；指定所要确定的等势面的水头值h0；等水头结点构成等 势面；
(1.3)对于渗流场的任意一个单元或第ie个单元，做以下工作 (1. 3. 1)判断各结点水头值与h0的关系；
各结点水头值与ho之间的差值有以下四种情况
①该单元上的所有结点的水头值都大于ho；
②该单元上的所有结点的水头值都小于ho；
③该单元上的结点的水头值有大于ho的，也有小于ho的；
④该单元上的结点的水头值有等于ho的，也有不等于ho的； (1.3.2)内插水头，求出等水头结点；
对于(1. 3. 1)中前两种情况①②，因不可能出现大小为h0的等势面，所以不再计算； 而对于(1.3. 1)中后两种情况③④，使用内插方法求出值为h0的水头等势面片；单元 内等水头结点构成的面片称为等势面片，
(1. 3. 3)求出该单元所有的等水头结点；剔除等水头结点中重复的点，还剩N个点，有 六种情况，N的值在广6之间取值，
其中，N=I时只是一个点、N=2时只是一条线，形不成等水头面片，通过流量为零，对流 量计算结果没有影响，因此不再处理；
当N大于2时，这N个点所构成的空间N边形就是该单元上的水头等势面片；
(1. 3. 4)水头等势面片形成后，当N=4时等势面片是空间四边形，此时依如下步骤处
理
(1.3. 4. 1)此时利用常规空间二维有限元法中的等参单元法，建立局部坐标系 (1. 3. 4. 2)用常规的内插法求出此二维单元积分点即单元内部的高斯积分点在三维单元即第ie个单元中的位置坐标；
(1. 3. 4. 3)用常规三维单元的计算方法求出第(1. 3. 4. 2)步中高斯积分点在三维单元 中的水力梯度，水力梯度有X、1、ζ三个方向的分量；
(1. 3. 4. 4)用常规二维单元的计算方法求出平面局部坐标系积分点所代表的面积，该 面积在整体坐标三个坐标平面中也有X、1、ζ三个方向的分量；
(1. 3. 4. 5)相应的面积分量与水力梯度的分量相乘，相乘后相加求和，其和值即为通过 该面片的渗流量；
(1. 3. 5)水头等势面片形成后，当N=3时等势面片是空间三角形，将三角形视为退化的 空间四边形，此时可按第(1. 3. 4)步的方法计算通过该面片的渗流量； (1. 3. 6)水头等势面片形成后，当N=5时等势面片是空间五边形； 连接空间五边形上不相邻的两个点，将此五边形分解为一个四边形和一个三角形； 对于四边形，按第(1. 3. 4)步的方法计算通过该小面片的渗流量； 对于三角形，按第(1. 3. 5)步的方法计算通过该小面片的渗流量； 将两者求和，和值即为通过该五边形面片的渗流量； (1. 3. 7)水头等势面片形成后，当N=6时等势面片是空间六边形； 连接空间六边形上间隔两个点的两点，将此六边形分解为两个四边形； 对于两个四边形，均按第(1. 3. 4)步的方法计算通过该小面片的渗流量； 将两者求和，其和值即为通过该六边形面片的渗流量； (1.4)对整个渗流域的每一个单元均做第(1. 3)步； (1. 5)将各面片的流量求和，即可求出经过指定水头等势面的流量； (1. 6)重新指定等势面的水头值，重复第(1. 2)、( 1. 3)、( 1. 4)、( 1. 5)步，求出经过新指 定水头等势面的流量。所述第(1.3. 2)中的内插方法的步骤为
(1. 3. 2. 1)取单元的任意一条边Li为对象，判断该边上各结点水头值与ho的关系，会 出现以下五种情况
(1. 3. 2. 1. a)该边上的两个结点的水头值都大于h0，不再计算； (1.3.2.1. b)该边上的两个结点的水头值都小于h0，不再计算； (1.3.2.1. c)该边上的两个结点的水头值一个大于h0，一个小于h0，则根据坐标和水 头关系内插出水头值为h0的点的位置；
设定两个结点的坐标分别为(xl，yl，zl)和(口^？，^，水头分别为！^和！^，则所求 水头值为h0的点的坐标(X，y, ζ)可由下式求出 x=xl+(x2-xl)/(hl-h2) X (hl-hO) y=yl+(y2-yl)/(hl-h2) X (hl_h0) z=zl+(z2-zl)/(hl-h2) X (hl-hO)
(1.3.2.1. d)该边上的两个结点的水头值一个等于h0，一个不等于h0，则(x，y，z)的 值就取等于h0的结点的坐标；
(1.3.2.1. e)该边上的两个结点的水头值都等于h0，则两个结点都是等水头结点。(1.3. 2.2)重复第(1. 3. 2. 1)步，对单元的各条边进行循环，判断每一条边是否存 在等水头结点，若有则求出等水头结点，将各边循环完后，求出该单元所有的等水头结点。
本方法为利用有限元求解地下水渗流量问题时的新的方法，优点
1.直接在等势面上积分，过水剖面与流向垂直，没有分量的干扰，计算精度高。2.由于通过水头等势面的流量只有一个方向，因此避免了中断面法中不但需要 确定水头梯度的方向还要确定面积向量的方向。3.避免了复杂渗流场中某一流量断面既有正向流动又有反向流动时所带来的流
量计算误差。4.避免了复杂渗流场中不同有限元剖分时流量断面不好确定的困难。5.不用对过水断面重新剖分。6.可以同时计算多个等势面的流量。
图1是中面法过水断面示意图
图2是等效结点流量法过水断面示意图； 图3是任意过水断面示意图4是等势面与单元相交时交点的数量(第ie个单元所有的等水头结点中剔除重复的
点)N=I时的示意图5是等势面与单元相交时交点的数量N=2时的示意图； 图6是等势面与单元相交时交点的数量N=3时的示意图； 图7是等势面与单元相交时交点的数量N=4时的示意图； 图8是等势面与单元相交时交点的数量N=5时的示意图； 图9是等势面与单元相交时交点的数量N=6时的示意图； 图10是等势面片局部坐标系； 图11是将三角形视为退化的空间四边形的示意图； 图12是五边形分解方法的示意图； 图13是六边形分解方法的示意图； 图14是基于等势面(线)的渗流量计算方法流程图； 图15是本发明的一个实施例的实验说明的示意图。
具体实施例方式
本发明所采用的方法是基于等势面的渗流量计算方法，其原理是
由于满足达西定律的渗流运动满足拉普拉斯方程
‘ 3H
所以，在渗流场中通过某断面s的渗流量的计算式为q =-Un^ds ,
方程式中S为过水断面；η为过流剖面的正法线方向的单位向量。
由相关理论可知在各向同性场中渗流流速和坡降的方向一致，形成流线，流线与 水头等势面正交，因此水头等势面的法线方向就是流速方向。
根据渗流场这一特性，如果计算过水断面采用水头等势面，渗流量计算正好满足 渗流量计算公式中的过水断面的法向与水流流速方向一致的要求，不会产生其他方向的分 量，因而计算精度最高。这也正是本发明基于等势面的渗流量计算的基本思想。本发明的计算过程如下计算流程如图14所示，
(1. 1)求解渗流场。用常规的有限元法求解渗流场，求出各单元结点水头值。(1.2)确定等势面的水头值。指定所要确定的等势面的水头值h0。等水头结点构 成等势面
(1. 3)对于第ie个单元(指任意一个单元)，做以下工作 (1. 3. 1)判断各结点水头值与h0的关系。各结点水头值与h0之间的差值有以下四种情况
①该单元上的所有结点的水头值都大于h0；
②该单元上的所有结点的水头值都小于ho；
③该单元上的结点的水头值有大于ho的，也有小于ho的；
④该单元上的结点的水头值有等于ho的，也有不等于(大于或小于)h0的。(1. 3. 2)内插水头，求出等水头结点。对于前两种情况(①该单元上的所有结点的水头值都大于h0 ；②该单元上的所有 结点的水头值都小于h0)，因不可能出现大小为h0的等势面，所以不再计算。而对于后两种情况(③该单元上的结点的水头值有大于h0的，也有小于h0的；④ 该单元上的结点的水头值有等于ho的，也有不等于ho的)，需要内插求出值为ho的水头等 势面片。具体内插方法
(1.3. 2. 1)取单元的任意一条边Li为对象，判断该边上各结点水头值与h0的关系，也 会出现以下几种情况
(1.3.2. La)该边上的两个结点的水头值都大于h0，不再计算； (1.3.2. l.b)该边上的两个结点的水头值都小于h0，不再计算； (1. 3. 2. 1. c)该边上的两个结点的水头值一个大于h0，一个小于h0，则根据坐标和水 头关系内插出水头值为h0的点的位置；
设定两个结点的坐标分别为(xl，yl，zl)和(口^？，^，水头分别为！^和！^，则所求 水头值为h0的点的坐标(X，y, ζ)可由下式求出 x=xl+(x2-xl)/(hl-h2) X (hl-h0) y=yl+(y2-yl)/(hl-h2) X (hl_h0) z=zl+(z2-zl)/(hl-h2) X (hl-h0)
(1.3.2. l.d)该边上的两个结点的水头值一个等于h0，一个不等于h0，则(x, y, ζ)的 值就取等于h0的结点的坐标。(1. 3. 2. 1. e)该边上的两个结点的水头值都等于h0，则两个结点都是等水头结
点ο(1.3. 2.2)重复第(1. 3. 2. 1)步，对单元的各条边进行循环，判断每一条边是否存 在等水头结点，若有则求出等水头结点，将各边循环完后，求出该ie单元所有的等水头结
点ο
上述内插方法仅仅是其中一种，也可以采用其它的内插方法，比如反距离加权 法、多元回归法、Kriging法、样条函数法等。(1. 3. 3)剔除上述结点中重复的点，还剩N个点，N的值在广6之间取值，即有如图 4 图9所示的六种情况，N的值在Γ6之间取值，有六种情况:N=1、N=2、N=3、4、N=5、N=6 ；
单元内等水头结点构成的面片称为等水头面片(或称等势面片)。其中，N=I时只是一个点、N=2时只是一条线，形不成等水头面片，通过流量为零， 对流量计算结果没有影响，因此不再处理。当N大于2时，这N个点所构成的空间N边形就是该单元上的水头等势面片。(1. 3. 4)水头等势面片形成后，当N=4时等势面片是空间四边形，此时依如下步骤 处理
(1.3.4. 1)此时利用常规空间二维有限元法中的等参单元法，等参单元法为公知技 术，建立局部坐标系，如图10所示
(1. 3. 4. 2)用常规的内插法求出此二维单元积分点(单元内部的高斯积分点)在三维 单元(即第ie个单元)中的位置坐标。(1. 3. 4. 3)用常规三维单元的计算方法求出第(1. 3. 4. 2)步中高斯积分点在三维 单元中的水力梯度，水力梯度有χ、y、ζ三个方向的分量。(1. 3. 4. 4)用常规二维单元的计算方法求出平面局部坐标系积分点所代表的面 积，该面积在整体坐标三个坐标平面中也有x、y、ζ三个方向的分量。(1. 3. 4. 5)相应的面积分量与水力梯度的分量相乘，相乘后相加求和，其和值即 为通过该面片的渗流量。(1. 3. 5)水头等势面片形成后，当N=3时等势面片是空间三角形，如图11所示，将 三角形视为退化的空间四边形，此时可按第(1. 3. 4)步的方法计算通过该面片的渗流量。(1. 3. 6)水头等势面片形成后，当N=5时等势面片是空间五边形。如图12所示，连接1点和4点，将此五边形分解为四边形1234和三角形145 ； 对于四边形1234，按第(1. 3. 4)步的方法计算通过该小面片的渗流量；
对于三角形145，按第(1. 3. 5)步的方法计算通过该小面片的渗流量； 将两者求和，和值即为通过该五边形面片的渗流量。当然，也可以连接2点和5点，将此五边形分解为四边形2345和三角形125，然后 对四边形2345和三角形125分别计算，求出其各自的渗流量。或者连接3点和5点，将此五边形分解为四边形2345和三角形125 ；或者连接1点 和3点，将此五边形分解为四边形1345和三角形123。将两者求和，和值即为通过该五边形 面片的渗流量。(1. 3. 7)水头等势面片形成后，当N=6时等势面片是空间六边形。如图13所示，连接1点和4点，将此六边形分解为四边形1234和四边形1456 ； 对于四边形1234，按第(1. 3. 4)步的方法计算通过该小面片的渗流量；
对于四边形1456，按第(1. 3. 4)步的方法计算通过该小面片的渗流量； 将两者求和，其和值即为通过该六边形面片的渗流量。也可以连接3点和6点，将此六边形分解为四边形1236和四边形3456 ；或者连接 2点和5点，将此六边形分解为四边形1256和四边形2345。将两者求和，其和值即为通过该六边形面片的渗流量。( 1. 4)对整个渗流域的每一个单元做第(1. 3)步。(1. 5)将各面片的流量求和，即可求出经过指定水头等势面的流量。(1.6)重新指定等势面的水头值，重复(1.2) (1.5)步。
下面是本发明的一个实施例的实验说明。对于一个矩形坝，当剖分网格如图1时，采用常规的中面法的断面可能就像虚线3 所示那样没有覆盖住流断面，计算出的渗流量结果不正确。而采用本发明(如图15所示)所给的出渗段之前的各个断面7均为与水流垂直的 过水断面，计算精度自然较高，且这些断面7是由水头值确定的，能够自动计算，且不需要 人工干涉。
权利要求
1. 一种基于等势面的地下水渗流量计算方法，其特征在于 计算步骤如下(1.1)求解渗流场；用常规的有限元法求解渗流场，求出各单元结点水头值；(1.2)确定等势面的水头值；指定所要确定的等势面的水头值h0；等水头结点构成等 势面；(1.3)对于渗流场的任意一个单元或第ie个单元，做以下工作 (1. 3. 1)判断各结点水头值与h0的关系；各结点水头值与ho之间的差值有以下四种情况①该单元上的所有结点的水头值都大于ho；②该单元上的所有结点的水头值都小于ho；③该单元上的结点的水头值有大于ho的，也有小于ho的；④该单元上的结点的水头值有等于ho的，也有不等于ho的； (1.3.2)内插水头，求出等水头结点；对于(1. 3. 1)中前两种情况①②，因不可能出现大小为h0的等势面，所以不再计算； 而对于(1.3. 1)中后两种情况③④，使用内插方法求出值为h0的水头等势面片；单元 内等水头结点构成的面片称为等势面片，(1. 3. 3)求出该单元所有的等水头结点；剔除等水头结点中重复的点，还剩N个点，有 六种情况，N的值在广6之间取值，其中，N=I时只是一个点、N=2时只是一条线，形不成等水头面片，通过流量为零，对流 量计算结果没有影响，因此不再处理；当N大于2时，这N个点所构成的空间N边形就是该单元上的水头等势面片；(1. 3. 4)水头等势面片形成后，当N=4时等势面片是空间四边形，此时依如下步骤处理(1.3. 4. 1)此时利用常规空间二维有限元法中的等参单元法，建立局部坐标系 (1. 3. 4. 2)用常规的内插法求出此二维单元积分点即单元内部的高斯积分点在三维单 元即第ie个单元中的位置坐标；(1. 3. 4. 3)用常规三维单元的计算方法求出第(1. 3. 4. 2)步中高斯积分点在三维单元 中的水力梯度，水力梯度有X、1、ζ三个方向的分量；(1. 3. 4. 4)用常规二维单元的计算方法求出平面局部坐标系积分点所代表的面积，该 面积在整体坐标三个坐标平面中也有χ、1、ζ三个方向的分量；(1. 3. 4. 5)相应的面积分量与水力梯度的分量相乘，相乘后相加求和，其和值即为通过 该面片的渗流量；(1. 3. 5)水头等势面片形成后，当N=3时等势面片是空间三角形，将三角形视为退化的 空间四边形，此时可按第(1. 3. 4)步的方法计算通过该面片的渗流量； (1. 3. 6)水头等势面片形成后，当N=5时等势面片是空间五边形； 连接空间五边形上不相邻的两个点，将此五边形分解为一个四边形和一个三角形； 对于四边形，按第(1. 3. 4)步的方法计算通过该小面片的渗流量； 对于三角形，按第(1. 3. 5)步的方法计算通过该小面片的渗流量；将两者求和，和值即为通过该五边形面片的渗流量； (1. 3. 7)水头等势面片形成后，当N=6时等势面片是空间六边形； 连接空间六边形上间隔两个点的两点，将此六边形分解为两个四边形； 对于两个四边形，均按第(1. 3. 4)步的方法计算通过该小面片的渗流量； 将两者求和，其和值即为通过该六边形面片的渗流量； (1.4)对整个渗流域的每一个单元均做第(1. 3)步； (1. 5)将各面片的流量求和，即可求出经过指定水头等势面的流量； (1. 6)重新指定等势面的水头值，重复第(1. 2)、( 1. 3)、( 1. 4)、( 1. 5)步，求出经过新指 定水头等势面的流量。
2.根据权利要求1所述的基于等势面的地下水渗流量计算方法，其特征在于 所述第(1.3. 2)中的内插方法的步骤为(1. 3. 2. 1)取单元的任意一条边Li为对象，判断该边上各结点水头值与h0的关系，会 出现以下五种情况(1. 3. 2. 1. a)该边上的两个结点的水头值都大于h0，不再计算； (1.3.2.1. b)该边上的两个结点的水头值都小于h0，不再计算； (1.3.2.1. c)该边上的两个结点的水头值一个大于h0，一个小于h0，则根据坐标和水 头关系内插出水头值为h0的点的位置；设定两个结点的坐标分别为(xl，yl，zl)和(口^？，^，水头分别为！^和！^，则所求 水头值为h0的点的坐标(X，y, ζ)可由下式求出 x=xl+(x2-xl)/(hl-h2) X (hl-hO) y=yl+(y2-yl)/(hl-h2) X (hl_h0) z=zl+(z2-zl)/(hl-h2) X (hl-hO)(1.3.2.1. d)该边上的两个结点的水头值一个等于h0，一个不等于h0，则(x，y，z)的 值就取等于h0的结点的坐标；(1. 3. 2. 1. e)该边上的两个结点的水头值都等于h0，则两个结点都是等水头结点； (1.3. 2.2)重复第(1. 3. 2. 1)步，对单元的各条边进行循环，判断每一条边是否存在等 水头结点，若有则求出等水头结点，将各边循环完后，求出该单元所有的等水头结点。
全文摘要
本发明公开了一种基于等势面的地下水渗流量计算方法，计算步骤如下(1.1)求解渗流场。用常规的有限元法求解渗流场，求出各单元结点水头值；(1.2)确定等势面的水头值；指定所要确定的等势面的水头值h0；等水头结点构成等势面；(1.3)对于渗流场的任意一个单元或第ie个单元，做以下工作(1.3.1)判断各结点水头值与h0的关系；本方法为利用有限元求解地下水渗流量问题时的新的方法，优点直接在等势面上积分，过水剖面与流向垂直，没有分量的干扰，计算精度高。由于通过水头等势面的流量只有一个方向，因此避免了中断面法中不但需要确定水头梯度的方向还要确定面积向量的方向。
文档编号G06F19/00GK102063577SQ20111000644
公开日2011年5月18日 申请日期2011年1月13日 优先权日2011年1月13日
发明者宋志宇, 宋海亭, 景来红, 李斌 申请人:黄河勘测规划设计有限公司