渗流力学教程1.doc

渗流力学教程吴梦喜中国科学院力学研究所2010.74目 录1张量的指标符号表示法51.1 指标符号51.2 求和约定51.3 克罗内克(Kronecker )符号：61.4 排列符号61.5 逗号的指标符号71.6 向量、张量及其坐标转换71.6.1 向量71.6.2 张量81.7 Navier-Stokes方程式92 渗流的基本概念与定律102.1多孔介质及其连续性假设102.2 流体的实际平均速度与渗流速度122.3 渗流的基本定律132.3.1 Darcy定律132.3.2 Darcy定律的推广152.4 非Darcy渗流162.5 岩土体的渗流性质162.5.1 岩土中水的存在形式172.5.2 岩土体的水理性质182.5.3 土体的颗粒级配与孔隙率0192.5.4 岩体的透水性0222.6 岩土体的非饱和渗透性232.6.1 非饱和土体的吸力232.6.2 土体的饱和度与吸力和渗透系数的关系243 渗流的偏微分方程与定解条件263.1单相流体渗流的连续性方程263.2两相不溶混渗流的连续性方程283.3 流体与骨架的状态方程283.3.1 不可压缩流体的状态方程293.3.2微可压缩流体的状态方程：293.3.3 气体状态方程293.3.4 岩体的状态方程293.3.5 土体的本构关系303.4单相流体渗流的偏微分方程303.5 定解条件303.5.1 初始条件313.5.2 边界条件314 地下水渗流的理论计算方法324.1 杜布依近似假定与潜水渗流的计算324.2 承压水的计算354.3 地下水向完整井的流动365 A finite-element algorithm for modeling variably saturated flows375.1 Introduction395.2 Basic equations405.3 Finite element formulation415.4 Temporal discretization435.5 Boundary conditions455.6 Discussion on conductivity over a time-step465.7 Algorithm validation and schemes comparison475.7.1 Example 1: one dimensional infiltration into a dry soil475.7.2 Example 2: two-dimensional transient variably saturated water table recharge505.7.3 Example 3: two-dimensional transient unconfined drainage515.8 Conclusion525.9 Acknowledgments535.10 References536 渗透变形与渗透破坏1参考文献1 1张量的指标符号表示法张量分析的主要目的是给人们提供一种数学工具，它可以满足一切物理学定律与坐标系的选择无关的特性。对于同样的一个物理学问题，用张量形式写出的方程与用其他数学形式写出的方程相比，不仅本质上具有普遍性，而且由于符号的对称与简洁，使得方程精练而完美。1.1 指标符号指标符号是指表示一组量的下标（或上标）。如，一组变量，.，可表示为（i=1,2,.,n）。利用指标符号，对三维空间中任意一点的直角笛卡儿坐标系中的坐标不再用x、y、z表示，而用、，或、，表示，计作或，对于三维问题，i=1,2,3，对于二维问题，i=1,2。用指标符号可以使很多繁琐的公式书写大为简化。1.2 求和约定指标与求和约定：(1)除了作特殊的说明外，用作上标或下标的拉丁字母指标，都将取从1到n的值；(2)若一项中有一个指标重复，则意味着要对这个指标遍历范围1,2,n求和。这就是爱因斯坦求和约定。如 （i=1,2,3）(11)显然，求和指标的符号可以任意更换，即(12)根据求和约定，需要记住3条规则：（1）若不是表示求和，不要将同一指标重复两次，即；（2）同一项中不允许同一指标重复两次以上，即没有意义。如果要表示，不能用，只能用来表示；（3）不表示求和的指标称为自由指标，自由指标在公式两侧应采用同一指标符号，如(13)式中，i是自由指标，j是求和指标。在三维笛卡儿坐标系中，式(13)表示一组代数方程式 (14)比较(14)与(13)可以看出，利用指标符号可以将某些公式的书写大为简化。1.3 克罗内克(Kronecker )符号：克罗内克(Kronecker )符号定义为(15)它是数学力学中常用的一个特定符号。采用这一符号可以将空间某一点对坐标原点的距离r的平方表示为(16)其中可以看出，乘以表示i指标被j指标代替而变为，同理，(17)克罗内克(Kronecker )符号还可以定义为单位向量的点积，即式中为i坐标方向的单位向量。1.4 排列符号 为排列符号，三个下标为1、2、3的顺时针排列时称为偶排列，三个下标为1、2、3的逆时针排列称为奇排列(18)由定义可知，,其余排列均为0。利用排列符号可以十分简单地表示行列式的展开式，如行列式(19)1.5 逗号的指标符号在某一个或几个指标符号之间加一逗号“，”，称为逗号指标符号，表示该变量对紧接逗号后面的一个或几个指标符号所表示的自变量取偏导数，如(110)(111)1.6 向量、张量及其坐标转换1.6.1 向量设向量在笛卡儿坐标系中的分量为，用指标符号可以表示为，为第i坐标轴上正方向上的单位向量两个向量的点积为(112)两个向量的叉积为(113)如右图所示，向量在、坐标系中的分量为，在、坐标系中的分量为，显然 (114)写成指标符号的形式 (115)其中 (116)容易证明 (117)式(115)和(117)同样适用于笛卡儿坐标系中的向量坐标转换。势场的梯度为一向量，它可表示为(118)向量场的散度可表示为(119)1.6.2 张量向量的分量仅与一个坐标有关，不同坐标系中其分量满足式式(115)或式(117)的坐标变换关系。将这一概念扩展，可以把张量定义为一组与若干个下标有关的、并满足一定坐标变换关系的量。如（i,j,k=1,2,.,n），n代表张量的维数,下标个数代表张量的阶数。对于2阶张量，应满足(120)对于3阶张量，应满足(121)利用张量的定义，也可以认为向量是一阶张量，但张量与向量不同，既不是向量，也不是标量，而是一个与向量有关的量。弹性力学中一点的应力状态就是典型的二阶张量。二阶张量将其指标对换得到的新张量与原张量相等，即则称这个二阶张量为对称张量，若则称为反对称张量。 若 ，则称为球张量。是一个单位球张量。张量对某一坐标的偏微分可表示为由于则：(122)因而是一个三阶张量。1.7 Navier-Stokes方程式Navier-Stokes方程是有粘性流体运动的控制方程，即(123)若用指标符号表示，式(123)可以简化为 i,j=1,2,3(124)式中，、为流速分量；、为作用力分量；P为水压力；为水的密度，为水的运动粘滞系数。同理，达西定律、水流连续方程、各向异性渗流方程均可用指标符号表示，本文的后续章节将讲到。2 渗流的基本概念与定律渗流是指流体通过多孔介质的流动。渗流现象普遍存在于自然界中。如地下水的渗流；石油、天然气、煤层气的渗流；动物体内的血液微循环和微细支气管的渗流；植物体内的水份和气体输送等。2.1多孔介质及其连续性假设多孔介质是指含有大量空隙的固体。也就是指固体中含有大量相互连通的孔隙、微裂缝或各种类型毛细管体系。图 1为一些天然和人造的多孔介质。图 1 一些天然和人造多孔介质0A海滩沙 B.砂岩 C.石灰岩 D.稞麦面包 E.木材 F.人肺 左.小球 右.砂石从实用的角度考虑，需要对多孔介质的几何性质进行描述。由于介质骨架的复杂性，用数学方程来对骨架的固体颗粒本身的几何特性进行直接描述几乎是不可能的。目前一般采用统计方法进行描述。一是描述宏观的平均特性，如孔隙率，另一种是以骨架的固体颗粒的某些统计性质，如颗粒级配曲线。孔隙率是多孔介质中孔隙的体积在总体积中所占的比重。定义孔隙率为连续函数，绕空间中一点p取一微元体，其孔隙体积为，则孔隙率为。取一系列不断缩小的这样的体积，得到一系列对应的孔隙率，当缩小到某个体积附近时，趋于稳定，如图 2所示，当这个体积进一步缩小，会达到一个临界体积，当再进一步缩小，便开始激烈振荡，这个体积叫做典型单元体；此时，有两种可能：若p点在孔隙上，则孔隙率为1。若在骨架上，则孔隙率为0。多孔介质中p点的孔隙率定义为(21)依据这个定义，孔隙率成为位置的连续函数，通过引入孔隙率和典型单元体的概念，实现了用虚构的连续介质代替了真实的、不连续的多孔介质n。图 2 典型单元体2.2 流体的实际平均速度与渗流速度 连续流体：由于流体力学是研究流体的宏观运动，没有必要对流体进行以分子为单元的微观研究，因而假设流体为连续介质。连续介质场：我们所研究的是流体的宏观运动，即大量流体分子的平均行为，把流体当作连续介质处理，意指任取一个流体微元体积都包含有许多个分子。平均速度：流体通过单位孔隙截面积上的速度，。，渗流速度：流体通过多孔介质截面上的速度，它与平均速度的关系为。流线(Streamline)：某一时刻的一条空间曲线，在它的上面每一点的切线方向都与该处的渗流速度向量重合。流线方程为。2.3 渗流的基本定律2.3.1 Darcy定律1856年，H.Darcy（达西）在解决法国Dijon城的给水问题时，用直立的均质砂柱进行了渗流的实验研究，该实验装置如图 3所示。根据实验结果，Darcy得出结论:流体通过砂柱横截面的体积流量Q与横截面积A和水头差成正比，而与砂柱长度L成反比，即图 3 Darcy实验装置图 （22） 其中：K为水力传导系数或渗透系数，它具有速度量纲，为水力梯度。根据水力学原理，总水头 （23）其中p是对应高度上的水压力，为流体的密度，为重力加速度，z为相对高度，v为流体速度。上式中的速度项相对于压力项和相对高度项一般很小，可以略去。因此达西定律又可以写成（24）其中v为渗流速度。石油工业中，往往是研究沿水平x方向的流动，对于这种情形，达西定律可简化为：（25）对各种不同的单相牛顿流体通过多孔介质流动的研究表明:以上darcy定律中的为绝对渗透率，只与多孔介质本身的结构特性有关，而与单相牛顿流体的特性无关。应当指出:这一结论只适用于牛顿流体，它对非牛顿流体是不适用的。渗透系数与渗透率的换算关系为（26）对于达西定律的一般使用范围的描述是这样的：对于粘性土不但存在起始水力坡降，而且当水力坡降超过起始水力坡降后，渗透速度与水力坡降的规律还偏离达西定律而呈非线性关系。另外，试验也表明，在粗颗粒土中（如砾石、卵石等），只有在小的水力坡降下，此类土的渗透规律才符合达西定律，而在较大的水力坡降下，水在土中的流动即进入紊流状态，渗透速度与水力坡降不符合线性关系。课程中对使用范围的描述是雷诺数，表示作用于流体微团的惯性力与粘性力的比值，雷诺数越小则意味着粘性力影响越显著，雷诺数小，意味着流体流动时各质点间的粘性力占主要地位，流体各质点平行于管路内壁有规则地流动，呈层流流动状态。雷诺数大，意味着惯性力占主要地位，流体呈紊流流动状态。因此，用雷诺数来描述达西定律使用范围是合适的，而不用使用上面文字进行赘述。但是，一般情况下，雷诺数是难以测量并且误差较大。另外，纳吉和卡拉地的实验结果为Re=5，并给出了不同渗流流态的计算公式：层流 过渡区 紊流阻力平方区 因此对于达西定律的使用范围是不很明确的。2.3.2 Darcy定律的推广各向同性介质中单相流体渗流，可以用张量的形式表示,（27）其中，为渗透张量，为流体的运动黏滞系数，i,j为直角坐标轴方向，为流体的密度，为重力加速度，Z为位置高程。渗透张量是元素为实数的对称张量，根据线性代数理论，元素为实数的对称矩阵，其特征值也是实数。也就是说特征方程组，（28）给出的特征向量为渗透张量的主方向，特征值为主方向的渗透系数。若将坐标轴方向与介质的渗透张量主方向一致，则该点渗透张量矩阵具有对角形式（29）用这种形式表示的张量称为对角张量。在这种情况下Darcy定律的推广形式为（210）2.4 非Darcy渗流所有的气体和简单的液体一般是牛顿流体（切应力与速度梯度成正比的流体称为牛顿流体）。聚合物溶液、悬浮液、水泥浆、浆糊等都是非牛顿流体。非牛顿流体一般是微可压缩性流体。唯一已知的明显可压缩非牛顿流体是泡沫流。非牛顿流体的渗流一般不符合Darcy渗流规律。牛顿流体在低速或高速时有时也不满足Darcy渗流规律。石油在低渗特低渗介质中具有启动压力梯度，非线性低速运动方程（211）其中G为启动压力梯度，为渗透率。石油在裂隙溶洞型介质中流动时，遵循高速非线性定律，常用下述运动方程描述 （212）其中n为介于1和2之间的指数。2.5 岩土体的渗流性质岩土渗流性质由作为渗流骨架的岩土性质和其中的流体的性质决定。由于介质的孔隙大小和形状及其分布异常复杂，很难用孔隙形式表征其渗透性，所以常用平均概念和综合性的参数代表其渗流性质。2.5.1 岩土中水的存在形式在岩土体的空隙中存在的水有气态水、液态水、和固态水三种类型。其中液态水又可根据其受力情况分为结合水、毛细水和重力水。此外，还有一种存在于矿物结晶内部及其间的矿物结合水。存在于为饱和岩土体空隙中的水蒸汽称为气态水。气态水可以随空气的流动而移动。它本身也可以由水汽压力大的地方向水汽压力小的地方迁移。当水汽增多达到饱和时，或当气温降低达到露点时，气态水便凝结成为液态水，成为地下水的一种补给来源。结合水通常是指束缚于岩石颗粒表面，不能在重力影响下运动的水。水分子是偶极体，一端带正电荷，一端带负电荷。由于静电引力的作用，带有电荷的岩土颗粒表面，便能吸附水分子形成结合水。依据吸附作用的强弱，结合水又有强结合水和弱结合水之分。紧靠岩石颗粒表面的水叫强结合水（又称吸着水）。这种水不能被植物吸收。结合水的外层水叫做弱结合水（又称薄膜水）。膜薄水一般情况下不能流动，但当施加的外力超过其抗剪强度时，最外层的水分子便发生流动。毛细水是指在表面张力作用下，沿着岩土细小空隙上升的水。毛细水可以从地下水面上升形成毛细水带，也可以脱离地下水面而独立存在，成为悬挂的毛细水。毛细水同时受重力和毛细力作用，能传递静水压力，其上升高度与岩土体空隙大小有关。重力水是指岩土颗粒表面不能吸引，仅受重力影响运动的水。重力水能传递静水压力，并且有溶解盐类的能力。井水、泉水都是重力水。它是水文地质、地下水渗流研究的主要对象。当岩土空隙中的水温低于0C时，液态水便结冰转为固态水。我国东北、青藏高原等地就有部分地下水是以固态的形式存在于岩石空隙之中，形成季节冻结区或多年冻结区。2.5.2 岩土体的水理性质岩土体能容纳一定水量的性能称为容水性。容水性的度量指标为容水度，即 （213）其中为容水度，也称之为含水量，w为岩土中所能容纳的水的体积，V为岩土的总体积。容水度与岩土体的孔隙率n相等。持水性是指在重力作用下释放出水时，由于分子力和表面张力的作用，能在其空隙中保持一定水量的性能。持水性的度量指标为持水度，即 （214）其中为岩土体的持水度，也称之为残余含水量，为在重力作用下保持在岩土体空隙中水的体积，V为岩土的总体积。给水性是指岩土体在重力的作用下能自由排出一定量的水的性能。给水性的度量指标是给水度。 （215）为在重力作用下岩土体排出的水的体积。 容水度等于持水度和给水度之和。 （216）土体中的全部孔隙被水充满时的含水量称之为饱和含水量，用表示。土体中水的体积与孔隙的体积之比称为饱和度,常用s表示。 （217） 当岩土体的孔隙全部被水充满时，饱和度等于1。 岩土体中的含水量等于残余含水量时的饱和度，称为残余饱和度 （218）定义有效饱和度于下 （219）2.5.3 土体的颗粒级配与孔隙率0土分成很多类型，如粘性土、砂性土、砾石、碎石土等，由于生成物质的多样性，沉积条件与生成过程的多变性，往往是不均匀的。由于水或风的搬运，土体多为松散的粒状堆积，属于多孔介质。所以在渗流问题中，经常引用颗粒分布曲线、孔隙率、容重等基本概念来描述土的组成和结构。颗分曲线如图 4所示，横坐标为颗粒粒径，纵坐标为小于对应粒径的土颗粒的重量在总重量中所占的百分比。从颗分曲线，可以看出土的类别，并可查用认为是研究问题的特征粒径或控制粒径，如，等。对应小于总重量10的粒径常被称为有效粒径，称为不均匀系数。这些控制粒径将在管涌和滤层规格中采用。图 4 土颗粒大小分布曲线孔隙率及容重，表示土的密实程度，特别是孔隙率对渗透性及渗透稳定性起着重要的作用。孔隙率n与一般土力学中的孔隙比e（孔隙体积与固项土颗粒体积之比）的关系为 （220） 如图 5所示，孔隙率与颗粒形状排列组合有关。理想化的同样大小球体以最松的立方体排列时的，以最密的菱形体排列时。若小颗粒填在大颗粒之间的孔隙中，则可降低n值。对于粗细混合料，孔隙率与不均匀系数有密切关系。 图 5 土体的孔隙率及其与均匀系数的关系 孔隙率不仅与土颗粒组成有关，而且与土的受压沉积等形成过程有关，如表 1所示，一般天然土的孔隙率介于0.30.6之间，泥炭土的孔隙率则高达0.60.8。虽然孔隙率已大到0.6以上，但由于架空的拱作用，松散堆积体仍能处于稳定状态。表 1 各种土体的孔隙率及渗透系数土体孔隙率n有效孔隙率渗透系数(cm/s)砂质砾0.250.350.200.25310-1510-2砾质砂0.280.350.150.20110-1210-2种砂0.300.380.100.15410-2110-3粉砂0.330.400.080.12210-2110-3砂粘土0.350.450.050.10510-3110-4粉壤土0.400.550.030.08510-4110-6粉粘土0.450.650.020.05110-6孔隙率固然与渗流有关，而孔隙通道的大小更为直接地影响透水性，并对研究管涌起直接的作用。对球体直径为d的立方体排列，其最小孔隙通道直径为，菱形排列时为；对于不均匀的土，若用有效直径表示时，其最小孔隙直径根据皮契肯分析为 （221）因为充填孔隙中的水，除能自由流动的重力水以外，还有受分子力作用的吸着水和薄膜水（或称强结合水和弱结合水），以及毛管水等各种形态的水，而且有水汽和气泡存在，所以真正作为水流动的孔隙必然减小。我们把这种能供给水流动的孔隙称为有效孔隙，它不仅与土粒和其间的孔隙大小有关，也受水头压力大小的影响。有效孔隙影响土的渗透性，其值可参见表 1。内部渗流有效孔隙体积还没有测定的方法，常用在重力作用下自由排出水的体积来近似代替，排水有效孔隙率就是给水度。2.5.4 岩体的透水性0岩石和土的主要区别是破碎程度的不同，甚至没有精确的界限。不过我们在这里要介绍的岩体是属于裂隙介质。岩块本身一般透水性很小，一般岩块的孔隙率约为0.010.1，渗透系数小于110-7cm/s，较之裂隙透水性要小56个数量级，因此可认为不透水，只是互相连通的缝隙发生渗流。当岩体有三组裂隙时，其总的孔隙率表示为 （222）式中，b为裂隙宽度，B为裂隙间距。 形成裂隙的原因主要是构造应力和受大气长期作用所致。垂直裂缝多是构造应力造成的。接近水平方向的层理缝是由沉积造成的。在不太深的表层，由于风化和卸荷作用，而使构造缝和层理缝张开。裂隙介质的岩体具有强烈的各向异性。几种岩体的不同裂隙状况见图 6。岩体中有时还有断层存在，其透水性更强，情况更加复杂。这些裂缝的壁面同样都有吸着水存在，缝开口宽度小于110-6cm时，一般就不会有重力水和自由水流动。图 6 不同岩石的裂隙分布(a)片麻岩 (b)云母片岩 (c)千枚岩 (d)石灰岩 (e) 石英岩 (f) 风化白云岩有时固体介质既有裂隙又是多孔的，即所谓双重介质。例如岩溶地区，由于溶洞的存在就属于多孔性的裂缝岩体，裂缝的粘土则属于裂缝多孔介质。此外对于介质在荷载或孔隙水压力改变的情况下发生体积变形时，则应结合渗流场一并考虑，即所谓的耦合问题。2.6 岩土体的非饱和渗透性2.6.1 非饱和土体的吸力当多孔介质的孔隙中有两种不溶混的流体接触时，两种流体之间的压力存在不连续性。这种压力突变值的大小，取决于该点处界面的曲率。这个压力差称之为毛细压力。 （223）其中，为非湿润相中的压力，为湿润相中的压力。 （224）其中为液面膜的表面张力，为液面膜与孔隙界面的夹角，r为毛细管的半径。当岩土孔隙中的湿项流体为水，非湿项流体为空气时，我们用毛管水头来表示这种压力差，毛管水头即表示为水头高度的孔隙气压力与水压力之差。 （225）其中为毛管水头，为水的密度。 各种土体的毛管水高度表 2所示。表 2 各种土体中毛管水高度土体毛管高度（m）砂0.030.1细砂0.10.5粘质砂0.52.0黄土2.05.0粘质土5.010.0非饱和土中的吸力，就是指非饱和土孔隙之中水与空气两相之间的压力差。主要由毛管压力构成，还包括基质吸力。2.6.2 土体的饱和度与吸力和渗透系数的关系土体的饱和度小于1时，由于孔隙中同时存在孔隙水和气体两种流体，流体界面之间存在表面张力，导致土体中产生吸力。吸力是非饱和孔隙水压力的相反值。由于土体中的孔隙大小不一，当饱和度降低时，重力水首先从大的孔隙中排出。因此，饱和度降低，水和气界面的孔隙半径也减小，因而吸力增加。饱和度或含水量与孔隙水压力之间的关系也称之为土水特征关系，有很多经验公式来描述。如下式所示的Van Genuchtens (1980) (VG)模型， （226）其中，和 n 为形状参数, 。图 7为某黏壤土（clay loam）的土水特征曲线。饱和含水量0.38，残余含水量为0.15。1.82和 n2.43。图 7 某黏壤土的土水特征曲线定义土体非饱和时的渗透系数与饱和渗透系数之比为相对渗透系数，相对渗透系数可以用与VG公式配套的Mualems (1976)公式描述（227）图 8为某黏壤土的渗透系数和含水量的关系。图 8 某黏壤土的渗透系数和含水量的关系 值得注意的是，当土体的含水量低于残余含水量时，由于土体中的不存在重力水，因此，土体中不存在渗流，也不存在吸力，渗透系数也是没有意义的。经验公式在有效饱和度较低时计算出来的吸力和渗透系数，一般情况下都是不符合实际的。3 渗流的偏微分方程与定解条件3.1单相流体渗流的连续性方程流体在多孔介质中流动，遵守质量守恒定律，此质量守恒所满足的方程即为连续性方程。在控制体上，任取一体元，如图 9所示，其表面为，n为边界面的外法线向量。单位时间内通过表面的流体质量：；由于非稳态渗流引起密度随时间变化，设介质的孔隙率为，引起内质量增加量，故整个内地质量增加量为；图 9 微元体若控制体内有源（汇）分布，其强度为q，则单位时间产生的流体质量为，则总质量为；由质量守恒定律得积分形式的连续性方程（31）由奥-高公式（也称散度定理：矢量场通过任意闭合曲面的通量，等于该曲面所包围的体积内矢量场的散度的积分），可以将上式化简为： （32）此即为单相流体的连续性方程。3.2两相不溶混渗流的连续性方程 对于油水两相不溶混的渗流，以上标o和w分别表示油相和水相的量，并用表示密度，s表示相饱和度，则有油相和水相的连续性方程分别为 （33）其中，。但是上式形式的两相流连续性方程在实际使用上不太方便，因为方程中油和水在地层条件下的密度，一般是未知的，而且也不是渗流力学中很关心的变量。另一方面，地层体积系数B在工程上有较成熟的计算方法。可将密度改用地层体积系数B表示， 两相流连续性方程写成： （34）将达西定律代入联立得到： （35）其中：和，分别为油和水的相对渗透率，为解方程还需要补充两个饱和度方程：，。3.3 流体与骨架的状态方程状态方程：指压力、密度、温度或其他热力学参数之间关系的方程3.3.1 不可压缩流体的状态方程不可压缩流体的状态方程： ， (36)3.3.2微可压缩流体的状态方程： 等温条件下微可压缩流体的状态方程可表示为(37)其中为流体的压缩系数。3.3.3 气体状态方程理想气体的状态方程： (38)其中V为气体的体积，T为绝对温度，R为摩尔气体常数，n为气体的摩尔质量。值等于8.31451 J（molK）。它与气体的性质和状态无关，故摩尔气体常数也称为通用气体常数。实际气体的状态方程： (39)z为气体的偏差因子。3.3.4 岩体的状态方程流体赋存于岩体之中，岩体的应力状态因流体压力的改变而变化，因而引起岩体的体积发生改变，其状态方程可以描述为： (310)其中，P为流体的压力，为岩体的体积，为岩体的压缩系数。由于岩体的固体颗粒不可压缩，其体积变化等于孔隙的体积变化，即 (311)3.3.5 土体的本构关系土体的本构关系即土体的应变与应力之间的关系，具有非线性、弹塑性等性质。异常复杂。土体的本构方程应符合的若干原理。根据大量实验数据建立本构方程时，所建立的关系应符合下列原理:1.坐标不变性原理；2.物质客观性原理；3.物质同构性原理；4.相容性原理；5.量纲不变性原理。本构关系的具体内容请参见有关土力学著作。3.4单相流体渗流的偏微分方程微可压缩流体的状态方程可以表示为（39）骨架的状态方程可表示为 （310）设，则忽略高阶分量（311）式（32）中的第一项可改写为 （312）则考虑骨架体积变形的微可压缩流体的单向渗流偏微分方程可表示为： （313）3.5 定解条件 一个非恒定的定解问题，包括初始条件和边界条件。3.5.1 初始条件 在初始时刻，压力P满足 （317）其中为已知函数。3.5.2 边界条件 第一类边界条件： （317）为已知函数 第二类边界条件：（318） 第三类边界条件：（319）其中h为系数。4 地下水渗流的理论计算方法4.1 杜布依近似假定与潜水渗流的计算杜布依(Dupuit,1983)研究地下水缓变流动时，由于渗流自由面的坡度很小，认为可以假定沿深度方向的铅直线是等势线，即测压管水头h=常数，以简化问题。如图 10所示的二维稳定渗流问题，自由面是一条流线，正确的等势线及流速分布如图中(a)所示。按照达西定律，沿自由面上任意点的流速应为（4-1）图 10 杜布依假定的说明因为水平倾斜角Q角很小，杜氏建议用代替，也就是相当于假定等势线是铅直线，即流速是水平向的，而具有静水压力分布，如图(b)所示。基于这种假顶时，平均流速与单宽流量为（4-2）（4-3）杜氏假定能使原来问题的多个独立变量(x,z)减为一个，以代替了，在（43）中的z不再以独立变量出现而使问题变为一元水流问题。杜氏假定实际上是略去了剖面上的垂直流速分量，该值在水平不透水底面上，到自由面上为(等势线与流线正交)。因此自由面变化很陡的地方将产生大的误差。应用杜氏假定与水平不透水层上的缓变无压渗流，如图 11(a)所示，对式（4-3）进行积分可得杜布依公式（4-4）或写成（4-5）式中、为上下游水深；h为从下游起任意距离x处的水深。由上式可解得任意距离x处的水深（4-6） (a) (b)图 11 杜布依假定的渗流水深变化 (a) 水平不透水层；(b) 倾斜不透水层若不透水底面不是水平面而是倾斜面，如图 11(b)所示，按照杜氏假定，任一铅直线上的平均流速为（4-7）这里，J为自由面上的水力坡降。若以H表示倾斜面的水深，并考虑到测压管水头h与水深H的关系，以及底坡，则得（4-8）及（4-9）分析上式就可以得出三种底坡情况，的四种地下水面曲线，分析方法与地表水的稳定缓变流一样，即时水面为回水或壅高曲线，时水面为降落曲线。 最后需要指出，渗流自由面流出下游边界面时常高出下游水面而存在一段自由渗出高度，如图 12中虚线为矩形土坝和井的实际自由面位置（数值计算）与杜氏假定的曲线比较，可知在出渗附近局部急变区发生较大误差。图 12 杜布依假定在急变流区的误差4.2 承压水的计算如图 13所示的厚度变化的承压含水层，其承压水非均匀流的计算式为（4-10）式中、分别为上、下游断面处承压含水层的厚度，、分别为上、下游断面处承压含水层的水头。(a) (b)图 13 含水层厚度变化时的承压水 在地下水水力坡度较大的地区，有时会出现上游是承压水，下游由于水头降至隔水顶板以下而成为无压水的情况，如图 14所示。这种情况可以分段计算，根据流量连续来连接两段的计算公式。图 14 承压无压流4.3 地下水向完整井的流动从井中定量抽水，经过一段时间的非稳定运动后，降落漏斗扩散到边界，周围补给量等于抽水量，则地下水程稳定运动。图 6 潜水完整井不同过水断面流量相等，并且等于井的流量，可得（4-11）即：（4-12）分离变量（4-13）积分得：（4-14）当,时（4-15）则有（4-16）这就是著名的Dupuit稳定井流潜水完整井出水量计算公式。5 A finite-element algorithm for modeling variably saturated flows Abstract: A general numerical algorithm is developed to solve Richards equation, in which a mass-conservative, modified head based scheme (MHB) is proposed to approximate the governing equation, and mass lumping techniques are used to keep the numerical simulation stable. The MHB scheme is compared with the modified Picard iteration scheme (MPI) in a ponding infiltration example. Although the MHB scheme is a little inferior to the MPI scheme in respect of mass balance, it is superior in convergence character and simplicity. Fully implicit, explicit and geometric average conductivity methods are performed and compared, the first one is superior in simulation accuracy and can use large time-step size, but the others are superior in iteration efficiency. The algorithm works well over a wide variety of problems, such as infiltration fronts, steady-state and transient water tables, and transient seepage faces, as demonstrated by its performance against some published experimental data. The algorithm is presented in sufficient detail to facilitate its implementation. Key words: Richards Equation; variably saturated flow; finite element method; transient water table; infiltration. 15.1 IntroductionThe standard approach to model variably saturated flow is using numerical methods to solve Richards Equation (RE). Most common approaches use finite difference or finite element spatial approximation with low order time integration schemes to solve Richards Equation, which are usually expressed in three standard forms: pressure head based, moisture content based and mixed-form where both variables are employed. The pressure based form can be adopted to deal with both saturation and un-saturation cases. However, for highly non-linear problems, such as infiltration into very dry soils, these methods may suffer from mass-balance error, poor iterative efficiency and convergence problems (Celia et al., 1990, Mcbride et al., 2005). The reason of poor mass balance resides in the time derivative term (Celia et al., 1990). While the derivatives of the moisture content (=moisture content, t=time) and (=pore pressure head) are mathematically equivalent in the continuous partial differential equation, their discrete counterparts are not, and the nonequivalence of the discrete form is exacerbated by the highly nonlinear nature of . This leads to serious mass-balance error which grows with the time-step size. Milly (1985) presented a mass-conservative numerical scheme in which a specific soil moisture capacity value (C) averaged over each element during each time steps was used. This approach, coupled with mass lumping (Neuman, 1999), ensures effectively global mass balance in the pressure head based equation. As for the moisture content based form, perfectly mass conservative discrete approximations can be applied, but this form degenerates under fully saturated conditions, since a pressure-saturation relationship no longer exists. It cant also be adopted to simulate layered soils since the moisture content is discontinuous across the contact face of different soils. A numerical technique for approximating the mixed-form equation has been developed to minimize the mass-balance errors and enhance computational efficiency, e.g. Celia et al. (1990) proposed a modified Picard iteration scheme that ensures mass balance by directly evaluating the change of the moisture content in a time step from the water pressure head, which was shown to provide excellent mass balance in modeling unsaturated problems with sharp wetting fronts. The mass balance error problem reported by some other researchers (Kavetski et al., 2001, Hao et al., 2005) in applying this method to free drainage problems arose actually from the error in calculating the unsaturated boundary flow flux, and the error decreased rapidly when a smaller iteration tolerance was employed. A new convergence criterion was introduced in the modified Picard iteration method and was found to be computationally more efficient (Huang et al., 1996). However, the use of mass-conservative methods does not guarantee good iterative efficiency, and thus the main purpose of the present paper is to provide a general, mass-conservative and comput