中信工程设计建设有限公司 湖北武汉 430010
摘要：分析存有自由面的无压渗流的问题时，由于渗流自由面的位置是未知的，需要迭代求解，因而这类问题的分析是非线性的。在饱和—非饱和联合渗流分析中，由于非饱和区的渗透系数函数不再是定值，而是一个与孔隙水压力相关的函数，因而这类问题的分析也是非线性的。目前的通用渗流软件大多数是基于饱和—非饱和联合渗流理论的，在分析无压渗流的问题时，非饱和区的渗透系数函数的取值比较困难。而实际工程中通常重点关注饱和区水的流动，而忽略非饱和区的流量。基于上述考虑，本文基于GEO-SEEP模块分析无压渗流时，将可能存在自由面的土层的渗流系数取值采取与Bathe法近似的方法进行简化处理，将自由面以上的非饱和区域的渗透系数取值为原渗透系数的1/1000，而自由面以下的饱和区保持原来渗透系数不变。并将该简化方法应用到单层均质土层和双层地基的二维无压渗流计算中，与解析解进行了对比。结果表明，这种简化方法的处理具有一定的合理性，并且简化了非饱和区的渗透系数函数的取值方法，为利用通用软件求解存在自由面的无压渗流提供了一定的参考价值。
关键词：自由面；无压渗流；GEO-SEEP模块；非饱和区；渗流系数函数；简化
引 言
在岩土体水力学研究中，为了便于进行数学处理和实际应用，通常采用线性的数学模型（如一般的渗流控制方程是以线性的达西定律推导出来的线性偏微分方程）。但由于受各种因素的综合影响，渗流运动的实际表现呈非线性。文献[1]将岩土体水力学中通常出现的非线性问题按非线性性质划分为以下7类：①渗流区域非线性问题（具有自由面的无压渗流）；②本构关系非线性问题（非达西渗流）；③渗透参数非线性问题（非饱和渗流）；④流态非线性问题（多相体渗流）；⑤介质非线性问题（多重介质渗流）；⑥相互作用非线性问题（渗流与溶质运移、温度、应力等耦合问题）；⑦多重非线性问题（上述两种或两种以上非线性问题的组合）。对于存在自由面的无压渗流的问题，由于渗流自由面的位置是未知的，因而这类问题的分析是非线性的，在数值模拟方法中必须通过迭代求解。无压渗流自由面的确定是渗流计算的主要内容；对于稳定渗流来说，该面上任一点水头 等于该点的位置高程；对于非稳定渗流来说，还应满足第二类边界条件的流量补给关系，即按式⑴计算渗流自由面下降时自由面
流入饱和区的单宽流量
.png)
[2]：
.png)
⑴
式中，为自由面上的水头值；为自由面变动范围内的土体给水度；为自由面外法向与垂线的夹角。
目前求解无压渗流自由面问题的有限元方法有两大类型。第一大类是移动（变动）网格法[3]，这类方法把自由面当作可移动边界处理，在迭代过程中不断改变自由面的位置，使网格发生变化，直到自由面位置稳定为止。采用这类方法虽然取得了不少成功的经验，但也暴露出方法本身的缺陷：①当初始自由面和最终自由面相差较大时，网格过分变形导致网格畸形；②渗流域内有结构物时网格移动会改变结构边界；③自由面附近有不均质介质时，网格移动时会破坏介质区分边界，使计算结果失真；④求解包含有自由面渗流域的应力分布时，求解范围应包括自由面以上区域而移动网格法不能采用同一网格分析渗流场和应力场。
为了避免移动网格法的缺陷，Neuman[4]提出了用固定网格法来进行有自由面渗流分析的伽辽金有限元方法。自此，固定网格法在有自由面渗流区域的研究与应用慢慢的开展起来。固定网格法主要可以归纳为以下四大类型：
第一种类型为流量调整法，如由Desai提出的剩余流量法[5]、张有天在剩余流量法的基础上改进的初流量法[6]、王媛等人提出的改进初流量法[7-10]、速宝玉等提出的节点虚流量法[11]、崔皓东等提出的改进节点虚流量全域迭代法[12]。
第二种类型为渗透矩阵调整法，如Bathe提出的单元渗透矩阵调整法[13]（以下简称为Bathe法）、李春华在Bathe法基础上改进的复合材料单元法[14]、王贤能在Bathe法基础上改进的高斯点法[15]、邓建霞在高斯点法的基础上改进的加密高斯点法[16]。
第三种类型为局部网格调整法，如吴梦喜等提出的虚单元法[17]、梁业国等提出的子单元法[18]、彭华等提出的弃单元法[19]。
第四种类型为截止负压法，如速宝玉等基于变分不等式理论提出的截止负压法[20]、张乾飞等在截止负压法的基础上提出的引用双参数罚函数的改进截止负压法[21]。
另外，王均星等提出了流形单元法[22]，用数值流形方法求解有自由面的渗流问题。
以上的各种类型的求解无压渗流自由面问题的固定网格法均是基于饱和渗流理论的基础上进行计算的。传统的饱和方法将自由面作为活动的边界，必须进行反复的迭代。流量调整法的每一步迭代中都要确定自由面的近似位置，解的稳定性较差；渗透矩阵调整法虽然不需要在每一步迭代中计算自由面位置，但每一步迭代计算中，都要重新集结、分解总体渗透矩阵，计算量比较大；局部网格调整法在迭代中不断变动局部网格，严格地说属于移动网格法；截止负压法以压力场为未知函数，虽然避免了对自由面边界的处理，但自由面附近区域的渗流计算不够精确。
无压渗流自由面的求解既可以利用饱和渗流理论进行，也可以采用饱和—非饱和渗流理论进行。饱和渗流理论只考虑自由面以下饱和区内水的流动；而饱和—非饱和渗流理论同时考虑饱和区和非饱和区内水的流动。饱和—非饱和方法能更好地反映水分运动规律，无需进行自由面调整。因为在饱和—非饱和渗流分析中，自由面不再是运动区域的边界，而是负压区和正压区的分界面，因而在耦合地表水分析、考虑降雨等补给条件方面的独特优势。但由于非饱和渗流的求解也为一非线性求解问题，其方程求解较为困难；另一方面，由于非饱和区的渗透系数函数难以获得，且当计算区域存在渗透系数差异比较大的介质界面时，计算结果不合理，因此饱和—非饱和联合渗流的应用还是受到了很大的限制。
然而，目前的不少通用软件渗流分析模块均是基于饱和—非饱和联合渗流理论建立起来的，比如Geo-studio软件的SEEP模块（以下简称GEO-SEEP）、ABAQUS计算软件等。但是非饱和区的渗透系数函数难以获得；而在工程的实际应用中，多数情况下均只考虑无压渗流自由面以下饱和区的渗流情况，而忽略自由面以上非饱和区的水的流动，即认为自由面以上的区域为干区。基于上述的考虑，本文基于GEO-SEEP模块，在求解存在自由面的无压渗流分析时，将可能存在自由面的土层的渗流系数取值采取与Bathe法近似的方法进行处理，将自由面以上（即负压区）的非饱和区域的单元渗透系数取值为原渗透系数的1/1000，而自由面以下的饱和区（即正压区）的单元保持原来渗透系数不变。经过这样的处理，非饱和区的流量与饱和区的流量相比，可以近似忽略不计，即使得非饱和区近似达到干区的状态；另一方面，这样的处理也可以使得非饱和区的渗透系数参数得以简化，容易获取。
1 二维非饱和土渗流分析原理
在饱和土中，引起水分转移的理的重力和水的压力。在非饱和土中，支配着土壤水在液态下整体转移的是重力和水的表面张力。非饱和土中的渗流与饱和土一样符合达西定律和连续方程，两者控制方程完全相同[23-24]。若将达西定律代入连续方程（忽略渗透过程中总应力的改变和土颗粒骨架的变形），并与总水头 作为未知量，当渗透的主方向与坐标轴一致时，非饱和土渗流的二维控制微分方程可表示为：
.png)
⑵
式中，
.png)
、
.png)
分别为
.png)
、
.png)
方向上的渗透系数；
.png)
为总水头；
.png)
为给定边界流量；
.png)
为体积含水量。
根据土水特征曲线（即体积含水量
.png)
和孔隙水压力
.png)
的关系曲线），体积含水量
.png)
和孔隙水压力
.png)
（以下简称为孔压）有以下关系式：
.png)
⑶
式中，
.png)
为土水特征曲线的斜率。
另外，总水头
.png)
和孔压
.png)
存在以下关系式：
.png)
或
.png)
⑷
式中，
.png)
为位置水头。
将⑷式代入⑶式，可得
.png)
⑸
再将⑸式代入⑵式，即可得到GEO-SEEP 模块中应用的二维饱和—非饱和控制方程式：
.png)
⑹
从⑹式可以看出，非饱和土中的渗透系数不再是常数，而是含水量的一个函数。
因此，利用GEO-SEEP 模块进行饱和—非饱和渗流分析时，首先需要获得非饱和土的渗透系数函数，即土体渗透系数和孔隙水压力的关系。测量非饱和土的渗透系数函数是一个相当复杂的过程，这个过程通常由于过低的渗透系数而时间较长。GEO-SEEP 模块提供了三种渗透系数函数的估计方法，通过测量或估计得体积含水量函数和饱和渗透系数来预测非饱和渗透系数函数。方法一是Fredlund法，由沿着整个体积含水量函数积分来得到非饱和渗透系数的方法组成；方法二是Green和Corey法，需要测定孔隙等级与对应的孔隙压力之间的关系；方法三是Van Genuchten法，提出一个闭合形式的与基质吸力相关的渗透系数函数。这三种方法均是基于实验曲线的拟合而提出的。
目前采用最多的是Fredlund法和Van Genuchten法。采用Fredlund法和Van Genuchten法获取非饱和土的渗透系数函数时，土的渗透系数函数与土体的体积含水量函数（也称储水量函数）有关。体积含水量函数表达了基质吸力变化时土的储水能力的变化情况。体积含水量函数可以通过室内试验获取，也可以通过拟合参数的方法估算。GEO-SEEP模块也提供了几种可参考的估算方法。但这几种估算方法都带有各自特有的拟合参数，需要用户确定，如果不进行相关的试验是难以获取这些参数。
考虑到在工程的实际应用中，多数情况下均只考虑无压渗流自由面以下饱和区的渗流情况，而忽略自由面以上非饱和区的水的流动，即认为自由面以上的区域为干区，这也是饱和渗流理论的观点。基于GEO-SEEP模块进行稳定无压渗流的计算时，本文着力于寻找一种简化非饱和土渗透系数函数的方法，使之既可以跳过获取非饱和土的体积含水量函数和渗透系数函数的过程，而又可以到达满足工程实际应用效果。这种简化方法是基于Bathe[13]法提出的，下面将对简化方法的具体思路作出描述，并将这种简化方法应用到两个实例当中，以验证其合理性。
2 基于Bathe法的简化方法
Bathe[13]将计算区域分为渗流区域（自由面以下区域）和非渗流区域（自由面以上区域），与自由面相交以及自由面以下的单元保持渗透系数 不变，而完全处于非渗流区域的单元渗透系数 取为原来渗透系数的1/1000的小值，表达式如下：
.png)
⑺
式中，
.png)
和
.png)
的定义同上述。
通过不断调整单元传导矩阵，模拟非渗流区域的作用，从而确定真实的渗流饱和区域。Bathe法将自由面的边界非线性问题转化为材料非线性问题，其原理接近于饱和—非饱和联合渗流模型，只是在非渗流区域的渗透系数取值的方法是不相同的。
在饱和渗流分析中，自由面就是总水头等于位置水头的一个法线流量为零的面，也是湿区和干区的分界面。而在饱和—非饱和联合渗流分析中，自由面就是饱和区和非饱和区的分界面，也是正孔压和负孔压的分界面，在自由面上孔压为零。为了使得非饱和区的流量与饱和区的流量相比可以近似不计，即让非饱和区到达近似干区的状态，可基于Bathe法对非饱和区土的渗透系数函数作如下的简化：
.png)
⑻
式中，
.png)
的定义同上述。
简化公式⑻可以在GEO-SEEP模块提供的渗透系数函数（Hyd.Conductivity Fn）中进行设置。渗透系数函数是土的渗透系数和基质吸力
.png)
（Matric Suction）的关系表达式。基质吸力
.png)
，
.png)
为孔隙气压力。当不考虑孔隙气压力的影响时，基质吸力就等于负的孔隙水压力。因此利用GEO-SEEP模块进行稳定无压渗流分析时，只用确定非饱和区土的渗透系数与负孔隙压力的关系即可。通过⑻式的简化方法，可以直接获取非饱和土的渗透系数和与负孔隙压力的相关关系，而不再依赖于体积含水量函数的获取，加快了参数获取的过程。
对自由面穿过的单元，单元高斯区域的渗透系数与其对应的高斯点和自由面的相对位置有关。当高斯点的位置高于自由面时，其在单元内对应的高斯区域的孔压为负孔压，则相应的渗透系数取值为；当高斯点的位置低于自由面时，其在单元被对应的高斯区域的孔压为正孔压，则相应的渗透系数取值为。如图1所示的采用三点积分的三角形单元，位于自由面以下的高斯区域1和2的渗透系数取值为，位于自由面以上的高斯区域3的渗透系数为；对于采用四点积分的四边形单元，位于自由面以下的高斯区域1和2的渗透系数取值为，位于自由面以上的高斯区域3和4的渗透系数为。这种处理与李春华的复合材料单元法[14]有异曲同工之妙。
.png)
图1 高斯区域渗透系数
Fig.1 Conductivity in Gauss regions
以下将以两个二维无压渗流实例验证上述简化方法在GEO-SEEP模块中应用的合理性。
3 单层均质土无压渗流
图2（a）所示位于水平不透水层上的单层均匀各向同性土层，渗透系数为0.001m/s。其底部高程为0m，上游水位
.png)
为12m，下游水位为
.png)
为2m，全长
.png)
为50m，
.png)
为总水头值。该土层的渗流为缓变无压渗流，根据裘布依假定，渗流自由面的近似解为：
.png)
⑼
断面单宽流量
.png)
为：
.png)
⑽
将数据代入⑼和⑽式，即可得到断面单宽流量为0.0014 m2/s；自由面的位置如图2（a）所示。基于GEO-SEEP模块利用上述的简化方法计算该渗流场，其计算结果如图2（b）所示，计算所得的截面单宽流量为0.0013987 m2/s，与上述的近似解结果接近，相对误差仅为0.09%。二者计算所得的自由面的位置也较为接近。
.png)
（a）解析解计算自由面
.png)
（b）GEO-SEEP计算自由面
图2 单层均质土无压渗流对比
Fig.2 Comparison of Unconfined seepage through single-layered homogeneous soil
4 双层地基无压渗流
图3为存在折线分界面的双层地基，上层为强透水层，下层为弱透水层。强透水层的渗透系数为
.png)
，弱透水层渗透系数为
.png)
，
.png)
，则地下水主要在强透水层中流动，
.png)
点为一奇异点，此处的地下水发生转向，水力坡降很大。因为地下水主要在强透水层中流动，可以将该渗流场等效为发生在不规则不透水边界上的无压渗流。
现根据裘布依假定，对该渗流场的解析解作以下的推导：
①水平不透水边界段的自由面位置和截面单宽流量
.png)
的计算公式分别同⑼和⑽式。
②对于斜坡不透水边界段，对任一
.png)
处，
.png)
，
.png)
；则
.png)
。
.png)
为总水头值，
.png)
为过水断面的厚度。
将
.png)
代入单宽流量微分公式：
.png)
，有
.png)
⑾
.png)
图3 双层地基无压渗流
Fig.3 Unconfined seepage through double-layered foundations
求解⑾式的非齐次线性微分方程，可得自由面位置公式：
.png)
⑿
式中，
.png)
为一常数。
将斜坡段左右边界条件分别代入⑿式：
.png)
时，
.png)
.png)
时，
.png)
联立解得
.png)
的表达式和截面单宽流量
.png)
为：
.png)
⒀
.png)
⒁
取特解
.png)
时，上下游水头满足
.png)
；此时自由面位置及单宽流量为：
.png)
⒂
.png)
⒃
单宽流量公式⑽和⒁式都含有未知数
.png)
和
.png)
，其中⒁式为隐式。利用以上两式求解
.png)
和
.png)
时，先在
.png)
范围内假定
.png)
的值，由⑴式可得到对应的
.png)
值。对于⒁式，使误差
.png)
，将假定的
.png)
和对应的
.png)
值代入误差
.png)
公式，即可得到对应的误差
.png)
。
.png)
取得最小值时对应的
.png)
和
.png)
值即为所求的解。计算过程用自编小程序实现。所求得的
.png)
和
.png)
值代入
.png)
表达式，即可得到
.png)
值用⑿式可计算斜坡段自由面的位置。
但⑿式也是一个隐式，使
.png)
。对任一
.png)
值，在
.png)
范围内假定
.png)
的值，并将
.png)
和
.png)
值代入误差
.png)
公式，即可得对应的误差
.png)
。
.png)
取得最小值时对应的
.png)
值即
.png)
处的水头值。计算过程用自编小程序实现。
取上游水位
.png)
为2.0m，下游水位为
.png)
为-3.0m，长度
.png)
为50m。强透水层的渗透系数为
.png)
为0.001m/s，弱透水层渗透系数为
.png)
为1E-6m/s。利用上述的解析解和GEO-SEEP模块求解该渗流场，结果如图4所示。解析解求得的截面单宽流量为4.46E-6 m2/s，GEO-SEEP模块计算所得的截面单宽流量为4.22E-6 m2/s，二者相对误差为5.24%。二者计算所得的自由面位置也较为接近，如图4所示。
结 语
本文基于GEO-SEEP模块，在求解存在自由面的无压渗流分析时，将可能存在自由面的土层的渗流系数取值采取与Bathe法近似的方法进行简化处理，将自由面以上（即负压区）的非饱和区域的单元渗透系数取值为原渗透系数的1/1000，而自由面以下的饱和区（即正压区）的单元保持原来渗透系数不变。并将该简化方法应用到单层均质土层和双层地基的二维无压渗流计算中，将计算所得的截面单宽流量和自由面位置与解析解进行了对比。结果表明，这种简化方法的处理具有一定的合理性，并且简化了非饱和区的
.png)
（a）解析解计算自由面
.png)
（b）GEO-SEEP计算自由面
图4 双层地基无压渗流对比
Fig.4 Comparison of Unconfined seepage through double-layered foundation
渗透系数函数的取值方法。另一方面，简化后的渗流场重点考虑饱和区水的流动情况，而忽略非饱和区的流量，满足实际工程的需求，为利用通用软件求解存在自由面的无压渗流提供了一定的参考价值。
参考文献：
[1] 柴军瑞.岩土体水力学非线性问题[J].岩土力学，2003（S2）：159-162.
[2] 关明芳，陈洪凯.渗流自由面求解方法综述[J].重庆交通大学学报（自然科学版），2005，24（5）：68-73.
[3] 张有天，陈平，王镭.有自由面渗流分析的初流量法[J].水利学报，1988（8）.
[4] 王媛.求解有自由面渗流问题的初流量法的改进[J].水利学报，1998（3）：68-73.
[5] 李新强.无压渗流有限元分析的改进初流量法[J].水利学报，2007，38（8）：961-965.
[6] 潘树来，王金凤，俞缙.利用初流量法分析有自由面渗流问题之改进[J].岩土工程学报，2012（2）：202-209.
[7] 张顺福，丁留谦，刘昌军等.基于子域积分的初流量法改进[J].水电能源科学，2013（3）：15-18.
[8] 速宝玉，朱岳明.不变网格确定渗流自由面的节点虚流量法[J].河海大学学报（自然科学版），1991（5）：113-117.
[9] 崔皓东，朱岳明.有自由面渗流分析的改进节点虚流量全域迭代法[J].武汉理工大学学报（交通科学与工程版），2009，33（2）：238-241.
[10] 李春华.稳定渗流有限元计算时采用固定网格法的初步研究[J].第三届全国渗流力学、第二届水利水电工程渗流学术讨论会，2007.
[11] 王贤能，黄润秋.有自由面渗流分析的高斯点法[J].水文地质工程地质，1997（6）：1-4.
[12 邓建霞.有自由面渗流分析的加密高斯点单元传导矩阵调整法[D].四川大学，2004.
[13] 吴梦喜，张学勤.有自由面渗流分析的虚单元法[J].水利学报，1994（8）：67-71.
[14] 梁业国，熊文林.有自由面渗流分析的子单元法[J].水利学报，1997（8）：34-38.
[15] 彭华，曹定胜.有自由面渗流分析的弃单元法及其应用[J].人民长江，1997（8）：38-40.
[16] 速宝玉，沈振中.用变分不等式理论求解渗流问题的截止负压法[J].水利学报，1996（3）：22-29.
[17] 张乾飞，吴中如.有自由面非稳定渗流分析的改进截止负压法[J].岩土工程学报，2005，27（1）：48-54.
[18] 王均星，吴雅峰，白呈富.有自由面渗流分析的流形单元法[J].水电能源科学，2003，21（4）：23-25.