渗流计算中的三维非结构化网格下的自由面适应网格方法

第42卷第4期四川大学学报（工程科学版）Vol．42No．4 2010年7月JOURNAL OF SICHUAN UNIVERSITY（ENGINEERING SCIENCE EDITION）
July2010文章编号：1009-3087（2010）04-0035-07
渗流计算中的三维非结构化网格下的
自由面适应网格方法
李连侠1，廖华胜1*，郑小玉2，刘达3，邹俊4
（1．四川大学水力学及山区河流开发保护国家重点实验室，四川成都610065；2．中国水电顾问集团成都勘测设计研究院，四川成都610072；
3．广东省水利水电勘测设计研究院，广东广州510610；4．珠江水利科学研究院，广东广州510611）摘要：在2维正交网格系统下自由面适应网格方法的基础上，提出了3维非结构化网格下的自由面适应网格方法。

该方法在迭代求解渗流自由面时，无需改动网格，而是用映射的方法去让自由面适应既定网格系统，也不需要对出逸边界进行特殊处理，其方法简单，易于程序实现；通过将数值计算结果与甘油试验及室内砂槽模拟结果相对比，表明了该方法在3维非结构化网格系统中捕捉浸润面的精度较高，拓宽了自由面适应网格方法在3维带自由面渗流问题中的应用范围。

关键词：渗流；自由面；自由面适应网格；3维非结构化网格
中图分类号：TV131文献标识码：A
Approach of Water Free Surface Adapting Grid System to Simulating Phreatic Surface Seepage Flow with3D Unstructured Grids
LI Lian-xia1，LIAO Hua-sheng1*，ZHENG Xiao-yu2，LIU Da3，ZOU Jun4
（1．State Key Lab．of Hydraulics and Mountain River Eng．，Sichuan Univ．，Chengdu610065，China；
2．HydroChina Chengdu Eng．Co．，Chengdu610072，China；
3．Guangdong Provincial Investigation Design and Research Inst．of Water Conservancy and Electric Power，Guangzhou510610，China；
4．Pearl River Hydraulics Research Inst．，Guangzhou510611，China）
Abstract：The method of water free surface adapting grid system with2D structured grids was extended to3D un-structured grids to simulate phreatic surface in seepage flow．In this approach，the grid system does not need to be modified or regenerated during nonlinear iterations，whereas free surface is mapped to grid nodes to determinate the final phreatic surface without handling seepage face specifically．The good agreement of computations and experi-mental results indicated that this method can capture the phreatic surface in3D unstructured grids with well preci-sion．This study extended the application range of adapting grid method in unconfined seepage flow．
Key words：seepage flow；phreatic surface；free surface adaping grid system；3D unstructured grids
在地下水输运问题中，潜水面问题经常遇到，这
收稿日期：2009－11－29
基金项目：国家科技支撑计划资助项目（2007BAC18B01）；教育部博士点新教师基金资助项目（20090181120013）；
四川大学青年基金资助项目（校青2008037）
作者简介：李连侠（1978－），男，讲师．研究方向：地下水数值模拟及工程水力学．E-mail：lianxiali@yahoo．com．cn
*通讯联系人可以归结为自由面的确定问题，至今仍是众多学者研究的重要课题［1－5］。

由于在计算中自由面边界是未知的，使得求解问题变得复杂和困难，通常用迭代逼近的方法来求其近似解。

目前为止，求解该类问题，主要有3种方法［6］：变网格法，固定网格法和无网格法。

其中，固定网格法以其诸多优势，最受青睐。

文献［1］在2维结构化网格系统中提出了一种新的固定网格方法———自由面适应网格方法，它在
自由面迭代求解过程中不改变网格而是让浸润线不断地去适应网格，采用网格映射的方法去确定浸润线，
直至最后达到允许精度，得到要求的浸润线分布；该方法具有计算简单，编程方便，精度较高，易于推广到3维网格系统下等优点。

虽然该方法从结构化网格推广到非结构化网格
的原理很简单，但由于网格系统变得复杂了，其几何关系也不像结构化网格下那样易于确定。

因此，作者对该方法进行了进一步的研究，并通过程序实现了3维自由面适应网格方法，
将其计算结果与甘油试验结果和砂槽模型试验结果进行对比分析，以期发挥网格映射方法的优势，拓宽其在3维非结构化网格系统下模拟复杂渗流自由面的能力。

1
自由面适应网格方法
1．1
自由面迭代方法
具有自由面三维渗流问题的数学模型可用下式
表示［7］
：
·（K H ）+Q H =S S H t ，t t 0，（x ，y ，z ）∈Ω
H （x ，y ，z ，t 0）=H 0（x ，y ，z ），t =t 0，（x ，y ，z ）∈ΩH （x ，y ，z ，t ）=H 1（x ，y ，z ，t ），t t 0，（x ，y ，z ）∈Γ1K n 2 H
n 2=q （x ，
y ，z ），t t 0，（x ，y ，z ）∈Γ2H （x ，y ，z ，t ）=z （t ），t t 0，（x ，y ，z ）∈Γ3μ H t = ·（K H ）－ H z （K z
+ε）+ε，
t t 0，（x ，y ，z ）∈Γ
3
式中，
K 为渗透系数张量；H =H （x ，y ，z ，t ）为水头；Ω为渗流计算区域；Q H 为源（汇）项；S S 为贮水
率；t 0为计算初始时刻；H 0（x ，
y ，z ）为t 0时刻Ω域的水头分布；Γ1为第1类（给定水头）边界，
H 1（x ，y ，z ，t ）为Γ1上的水头分布；Γ2为第2类（给定流量）边界；n 2为Γ2的法向方向；q （x ，y ，z ，t ）为Γ2上的
流量分布；Γ3为第3类边界，
即渗流自由面边界；ε为Γ3上的补给强度；μ为给水度。

2维正交网格下自由面适应网格方法在文献［1］中做了详细的介绍，该方法利用自由面上H =z 的特点将计算浸润线abcdef 映射在固定网格系统中的A 1A 2BC 1C 2DE 1E 2F 各节点上（图1），无需对出逸面做特殊处理，
很容易推广到3维正交网格系统中。

在非正交如非结构化网格系统下，自由面的映
射方法则相对复杂些，
但其迭代过程与2维正交网格下基本相同，整个迭代过程如下：
图1
自由面适应网格方法实现原理示意图
Fig．1
Free surface adapting to grid
1）假定一初始自由面，简单起见，可以取一水平面。

2）求出自由面与计算区域内哪些单元相交，即
定出自由面单元，假定其满足H =Z ，Z 为单元形心的高程坐标。

3）求解渗流场，求出渗流单元的水头分布。

4）如果所有自由面单元均满足|Z －H |小于给定误差，则迭代完成，到第7）步，否则继续第5）步。

5）将这些单元的形心投影到水平面上，将其三角化，生成三角形网格。

6）将生成的三角形网格顶点的Z 坐标变成计算出的水头值，则变成立体网格，该网格就是所谓的自由面；网格节点所在的单元就是调整后的自由面
单元，
假定其满足H =Z ，回到第3步继续迭代。

7）满足收敛控制标准，迭代结束，得出结果。

1．2
涉及的几何问题
虽然从正交网格到非正交网格的原理很简单，但由于网格系统变得复杂了，其几何关系也不像正交网格下那样易于确定。

从1．1节的方法介绍中可以看出，在确定自由面单元的时候，需要用到面与体的相交；在将投影后的自由面单元形心三角化时，需要用到平面三角形网格生成技术。

下面对解决这些几何问题的方法进行论述［8］。

1．2．1
平面与凸多面体的相交
要解决凸多面体与一个平面的相交，就必须首
先解决多边形与平面或三角形相交的问题。

1）三角形和平面的相交
假设平面P 用点法式来表示，
即P ：（P ，n ），其中，
P 为其上的一点，n 为其法向向量；三角形T 用3个顶点来表示，即T ：（Q 0，Q 1，Q 2）。

如果平面与三角形相交，那么三角形的一个顶点将与其它的两个顶
6
3四川大学学报（工程科学版）第42卷
点分别位于平面的两边。

如果计算每一个顶点与平面之间的有符号距离，并比较它们的符号，那么就可以立即判断出三角形与平面是否相交。

因为只需要知道三角形网格是否与平面相交，而且不需要知道其交点，所以实现起来较为简单。

2）三角形与三角形的相交
三角形与三角形的相交判断，与三角形和平面相交判断相似，不同的是，需要先求出交点，看交点是否在另一三角形之内。

这需要用到点是否在三维三角形内的判断，这里给出点是否包含于一般凸多边形的判断方法。

首先要求该点位于多边形所在的平面内，否则点必在多边形外。

设多边形的边为e i=V i+1－V i，并且其外边法向为n i=Perp（e i）。

当n i·（P－V t）＜0对所有的i都成立时，点P位于多边形内。

如果n
i
·
（P－V
t
）＞0对某些i成立，则该点位于多边形外。

如果n i·（P－V t） 0对所有i成立，且至少有1个i使得等号成立，则该点位于多边形的边界上；如果等式仅对1个i成立，那么P位于1个边上，但不是顶点；如果等式对两个i成立，那么该点是1个顶点；等式不能出现3次或者3次以上。

1．2．2离散点的三角化
1）凸包算法
由于在自由面迭代过程得到的投影点是离散的，无序的，必须首先找到边界点才能将这些离散点三角化，这又用到了2维凸包算法。

考虑一个有限点集Q，如果该集合的所有点都共线，那么其凸包就是一条线段；当至少有三个点不共线时，凸包是一个成为凸多边形所包围的区域；凸包的顶点必须是原点集的一个子集，将Q内的凸多边形的顶点标识出来就建立了凸包。

计算点集凸包的算法很多，比较常用的有Graham算法等［9］。

2）德洛奈三角剖分
德洛奈三角剖分（Delaunay Trianglation），也叫德洛奈三角化，是将平面上1个点集联接成三角形的方法［10］。

这些三角形具有如下特点：互不重叠；覆盖所有的点；每个三角形的外接圆不包含除构成该三角形的3个点的其它点集中的点。

对于给定的有限点集，德洛奈三角剖分开始于构建一个能包含这个点集的足够大的超级三角形，使得它能够包含最终的三角剖分的外接圆的并集；之后采用增量构建方法将每一个输入点P都被插入这个三角剖分内。

注意，上述过程并不一定保证所生成的三角形
满足第三个条件，但不影响自由面适应网格的迭代。

实际上，在网格自适应迭代过程中，生成的网格只需要满足前两条即可，但如果同时满足第三点则会使得自由面更光滑。

另外，对自由面单元形心投影后的离散点，不一定非要用三角剖分将其联接起来，可以是任意多边形剖分，只要能够满足不漏点不重叠，则该方法仍然是可用的，这也说明该方法的适应性较好。

1．2．3点在凸多面体内的检测
将三角化后的投影点Z坐标换成计算出的水头值后，要确定自由面单元，需要判断网格点在哪一个单元内。

要检测一个点P和一个具有不共面顶点V i的多面体之间的关系，与检测一个点与多边形的关系有相似之处，具体可以用如下方法来判断：设多面体的面包含于n·（X－V i）=0内，其中，V i是这个面的一个顶点，n i是该面的外法向向量。

如果n·（P－V
i
）＜0对所有i均成立，则点P位于多面体内。

如果
n·（P－V
i
）＞0对某些i成立，则该点在体外。

n·（P
－V
i
） 0对所有i都成立且至少有1个i使等式成立，那么该点位于多面体的1个面上。

如果等式出现3次或3次以上，那么点P位于1个顶点上，等号出现的次数就是共用该顶点的面数。

2计算实例验证
为了对3维自由面适应网格方法的计算精度进行验证，分别进行了均质模型和非均质模型两组计算，并将模拟结果和试验结果进行了对比分析。

2．1均质模型
选取一上游水位为6．0m、下游水位为1．0m的均质矩形坝进行验证计算［11］，其剖面图见图2（a）。

即其计算区域变为：〔XˑYˑZ〕=〔（0，4）ˑ（0，0．6）ˑ（0，6．5）〕；边界条件为：
H=6．0m，X=0，Z∈［0，6］；
H=1．0m，X=4，Z∈［0，1］；
其余物理边界为第二类计算边界
{。

用非结构化网格剖分计算区域，共包含35083个四面体单元（见图2（b）），初始浸润面假定为一与上游水面齐平的水平面。

该计算实例其实是1个3维网格的2维流动问题，旨在验证方法的可行性和精度。

计算的水头等值线图示于图2（c）中。

可见，流动结果是合理的。

73
第4期李连侠，等：渗流计算中的三维非结构化网格下的自由面适应网格方法
图2矩形断面坝体
Fig．2
Computational model
在Y =0．3m 处，作一竖直剖面，其与3维浸润面相交即该剖面沿X 方向的浸润线，该浸润线与甘油试验值和2维（网格尺寸为Δx =Δz =0．1m ）计算对比结果见图3，除了3维计算结果出逸处附近略微偏高外，
3者基本吻合。

图3矩形坝浸润线计算和试验值对比图Fig．3
Phreatic surfaces of computing and experiment
图4是自由面迭代过程。

图中的0号面为初
设自由面，
1，2，……分别代表第1，2，……次迭代后得到的中间结果（自由面），
20号面为最终迭代结果。

可以看出，任意给定一个初始自由面（图4中0
号自由面），迭代一步后得到1号自由面，与最终结
果相差很大，继续迭代调整，逐步向最终结果（图4中第20号自由面）收敛，这是个自由面适应网格自动调整的过程。

在迭代过程中，自由面是由三角网格构成的。

图4矩形3维算例自由面迭代过程Fig．4
Iteration of phreatic surfaces
图5给出了迭代过程中的一个自由面三角网
格及其水头等值云图。

可见，
自由面迭代过程中，对德洛奈三角剖分质量要求并不高，只需能覆盖整个
计算区域且无重叠即可。

图5
自由面三角网格在平面上的投影及其水头等值云图
Fig．5
Triangulation and head contour of phreatic surface during iteration
2．2
非均质模型
针对图6所示非均质渗流模型，进行了砂槽试验（装置见图7和图8）和数值计算。

砂槽的一侧为
灰塑料板，装有量测各部位的孔隙水压力的测压管
（纵向间距20cm ，垂向间距8cm ），另一侧为有机玻璃以便观察。

矩形槽长3．2m （试验中渗流区长度
8
3四川大学学报（工程科学版）第42卷
1．79m ），宽0．5m ，高2m 。

矩形槽由槽首、槽身和
槽尾三段组成。

各段之间用可移动的过滤网隔开，使槽身具有伸缩性，以便适应设计模型大小的变化。

通过试验率定，试验中采用的介质的渗透系数选定为：砂0．398cm /s ，碎石7．72cm /s 。

其中浸润面位置采用土壤水分仪率定获取，测点布置情况见图6（圆圈黑点标记），渗流量通过称重法得到。

图6不同介质分块3维试验平面示意图及自由面测点布置Fig．6Plan view of model and monitor points of phreatic
surface
图7渗流模型示意图
Fig．7
Layout of seepage experimental
model
图8
砂槽实景照片
Fig．8
Layout of sand simulation model
试验和模拟主要针对上游水位保持0．842m 不变，下游水位分别为（1）0．515m 和（2）0．262m 两个工况进行。

图9给出了工况1计算中所采用的非结构化网格剖分情况。

两个工况下水头等值线及渗流场计算结果见图10。

由图10可见，水头场的分布沿横向与纵向并不均匀，这正是3维流场的特点；在模型右侧的两块介质中，等势线从稀疏而近于垂直逐渐变密变弯曲，在
右侧砂块（第2块介质）的前部等值线最密，等势线向后发展越加弯曲；模型左侧的两块介质中，等势
图9非结构化网格剖分Fig．9
Unstructured grids used
图10工况1及工况2水头场及渗流场计算结果Fig．10
Computing hydraulic head and seepage field of case 1and 2
线的大致趋势与右侧相同，但此时等势线在左侧砂块（第1块介质）的尾部最密集。

这是因为，一般情况下，渗透系数越小，其浸润线越陡，水头变化越大，与试验测量结果吻合。

不同介质分块模型中，明显存在Y 方向上的渗流流动，使得流场扭曲。

各种工况下自由面计算结果见图11。

由图11可见，自由面适应网格方法正确地模拟出了分区非均质3维渗流的特点，自由面的分布沿横向和纵向均有变化。

在模型的前半部份，自由面从最初的近乎相平变化到呈明显右高左低的趋势；在模型的后半部份，自由面从右高左低逐步变化到近于平齐。

自由面在模型的中部，可见明显扭曲。

9
3第4期李连侠，等：渗流计算中的三维非结构化网格下的自由面适应网格方法
图12 14为3种工况下沿3个纵剖面（Y =0．490m ，Y =0．375m 和Y =0．125m ）的浸润线计算结果和试验结果比较情况。

计算得到的浸润面结
果随比试验实测结果略低，
但总体来说与试验结果均能较好吻合。

图11工况1及工况2三维自由面等值线计算结果Fig．11Computing 3D phreatic surfaces of case 1and
2
图12工况1及工况2在Y =0．490m 剖面处浸润线计算结果与试验结果对比Fig．12Phreatic surfaces of computing and experiment of case 1and 2（Y =0．490m
）
图13工况1及工况2职在Y =0．375m 剖面处浸润线计算结果与试验结果对比Fig．13Phreatic surfaces of computing and experiment of case 1and 2（Y =0．375m
）
图14工况1及工况2在Y =0．125m 剖面处浸润线计算结果与试验结果对比Fig．14Phreatic surfaces of computing and experiment of case 1and 2（Y =0．125m ）
4四川大学学报（工程科学版）第42卷
表1给出了计算与试验所得的渗流量与出逸面位置结果，二者基本一致，误差均在8%以内，工况2出逸位置试验结果与计算结果相差略大，与边界对试验观测造成的误差有关。

由上可见，作者提出的自由面适应网格方法能够在3维非结构化网格下较为精确地模拟出复杂渗流自由面。

表1渗流量及出逸面高程计算和试验结果对比Tab．1Seepage fluxes and exit positions from computing and experimental results
工况
渗流量/（L·s－1）出逸面高程/m
计算试验误差/%计算试验误差/%
11．561．531．960．520．515 0．5200 1
21．821．782．250．380．35 0．405．2 7．9 3结论
提出了3维非结构化网格下的自由面适应网格方法，在迭代求解变动的自由面时，无需改动网格，而是用映射的方法去让自由面适应既定网格系统，不需要对出逸边界进行特殊处理，其方法简单，易于程序实现；通过将数值计算结果与甘油试验及砂槽模拟结果相对比，表明了该方法在3维非结构化网格系统中捕捉浸润面的精度较高，拓宽了网格映射方法在带自由面渗流问题中的应用范围。

本文研究中在自由面适应网格时采用的是1阶线形映射方法，如果采用2阶或更高阶映射或加密网格的方法则可以进一步提高计算精度。

参考文献：
［1］Li Lianxia，Liao Huasheng，Liu Da，et al．A method of water free surface adapting grid system to simulate phreatic surface in seepage flow［J］．Journal of Sichuan Univeristy：Engineering Science Edition，2006，38（5）：76－81．［李连侠，廖华胜，刘达，等．渗流计算中求解浸润线的自由面适应网格方法［J］．四川大学学报：工程科学版，2006，38（5）：76－81．］
［2］Ning Bo，Long Haiyan，Chai Junrui，et al．Brief description of numerical method research of free surface seepage［J］．Advances in Science and Technology of Water Resources，2009，29（1）：85－89．［宁博，龙海燕，柴军瑞，等．具有自由面渗流问题的数值方法研究简述［J］．水利水电科技进展，2009，29（1）：85－89．］
［3］Chen Jiang．Method for seepage analysis with free surface based on ansys［J］．Geotechnical Engineering Technique，
2008，22（5）：236－240．［陈江．基于Ansys的渗流自由面计算方法［J］．岩土工程技术，2008，22（5）：236－240．］［4］Shu Zhongying，Deng Jianxia，Li Longguo．Application of element conduct matrix regulative method by encrypting gauss point in seepage analysis with free surface［J］．Journal of Sichuan University：Engineering Science Edition，2007，39（1）：48－52．［舒仲英，邓建霞，李龙国．加密高斯点单元传导矩阵调整法在有自由面渗流分析中的应用［J］．四川大学学报：工程科学版，2007，39（1）：48－52．］［5］Chen Yeying，Tang Hongmei，Chen Hongkai．Seepage free surface solution of the Three Gorges Reservoir bank and its application［J］．Port and Waterway Engineering，2006，11（395）：16－19．［陈野鹰，唐红梅，陈洪凯．三峡水库岸坡渗流自由面求解方法及应用［J］．水运工程，2006，11（395）：16－19．］
［6］Guan Mingfang，Chen Hongkai．Review of solution methods of seepage free surface of slope［J］．Journal of Chongqing Jiaotong University，2005，24（5）：68－73．［关明芳，陈洪凯．渗流自由面求解方法综述［J］．重庆交通学院学报，2005，24（5）：68－73．］
［7］Bear J．Dynamics of fluids in porous media［M］．New York：Dover Publications，INC，1972．
［8］Schneider P J，Eberly D H（美）．计算机图形学几何工具算法详解［M］．周长发，译．北京：电子工业出版社，2005．
［9］Liu Runtao．A new linear algorithm for computing the con-vex hull of simple polygon［J］．Journal of Engineering Graphics，2002（2）：120－126．［刘润涛．一种简单多边形凸包的新线性算法［J］．工程图学学报，2002（2）：120－126．］
［10］Liu Jianxin，Lu Xinming，Yue Hao．Fast algorithm for delaunay triangulation of simple polygon based on maximum triangle weights［J］．Computer Technology and Develop-ment，2006，16（7）：126－128．［刘建新，卢新明，岳昊．简单多边形快速Delaunay三角剖分算法［J］．计算机技术与发展，2006，16（7）：126－128．］
［11］柴军瑞．大坝工程渗流力学［M］．西藏：西藏人民出版社，2001．（编辑张琼）
14
第4期李连侠，等：渗流计算中的三维非结构化网格下的自由面适应网格方法。