地下水数值模拟技术

地下水数值模拟技术朱国荣编著南京大学水科学系2006年9月1第一篇有限单元法有限单元法（FiniteElementMethod—FEM）是求解数学物理问题的一种行之有效的数值方法，基本思想诞生于1943年（Courant）。二十世纪五十年代，有限单元法从结构分析的基础上发展起来，六十年代在工程计算领域得到了相当的推广。二十世纪六十年代中期，有限单元法被用来求解地下水运动问题，而昀早将该方法引入我国地下水研究的时间是六十年代末，其中，肖树铁、孙纳正、谢春红、陈崇希、薛禹群等人在这方面起到了重要的作用。二十世纪七十年代，以南京大学、武汉水利电力学院、山东大学和中国地质科学院为首的教学研究单位，在推广和应用有限单元法求解地下水问题上发挥了重大的作用，形成了目前广泛使用这种方法开展地下水流模拟的局面，也确立了有限单元法在地下水资源评价中不可动摇的地位。应该说，用有限单元法建立地下水数值模型是现阶段研究地下水运动规律，地下水资源评价、地下水管理、地下水溶质运移、包气带中地下水的运动，地面变形等方面的基础工作。显然，掌握这种基本的地下水模拟技术的重大意义也就不言而喻了。目前，与有限单元法具备同等价值的常用地下水数值模拟技术还有有限差分法（DefineElementMethod—DEM）、边界单元法（BoundaryElementMethod—BEM），这些方法的数学原理不同，离散方法各异，但昀终都将对微分方程的求解转化为求解一个水头方程组。本篇专门讨论有限单元法的原理和应用，通过对一组理想问题的讨论，由浅入深地介绍给予有限元法的地下水模拟技术。2第一章有限单元法的基本概念1.1有限单元法的基本分类从本质上说，有限单元法昀终把对定解问题的求解转化为对现行方程组得求解，而根据建立方程组得过程不同，有限单元法可以分为若干种，目前昀普遍使用的是雷利—里茨（Rayleigh-Ritz）法、迦辽金（Galerkin）法和单元均衡法。以下用求解地下水运动问题来简单介绍这三种方法的基本思想。1.1.1雷利—里茨有限单元法雷利—里茨有限单元法简称里茨法，它以变分原理和剖分插值为基础，把对描述地下水运动的定界问题的求解转化为对等价于该定解问题的泛函求极值问题，然后在各离散元确立水头插值函数以便建立水头方程组。1.1.2迦辽金有限单元法迦辽金法从昀小二乘法原理出发，根据剩余理论引入对应于地下水运动的定解问题的权函数，通过剖分后，确定各离散元的水头插值函数进而建立水头方程组。1.1.3单元均衡法单元均衡法简称均衡法，它以质量守恒为基础，首先把研究区域D离散成有限多个足够小的子域Di（图1-1中的阴影处），然后建立每个子域的流量平衡关系。根据这个思想，某一个研究时段内流进和流出子域的水量差应该等于子域内的含水层在该时段内由于水头下降或上升所释放或贮存的水量，据此可利用达西定律建立一个满足于定解问题所描述的水头方程组。1.2空间离散空间离散的实质就是将一个被研究的连续空间按照某种规则离散为一系列彼此相依的子空间，这些子空间称为单元，所以又将这些单元称为离散单元，这个过程又成为单元剖分。单元的形态根据研究区域的几何形状可以有不同的形态。对于二维问题，昀常用的单图1-1单元均衡法中的均衡子域3元是三角形单元，有时候也用矩形或四边形单元（图1-2中的b，c）。前者可以看作是后者的特殊情况。当然，在需要的时候还可以是其它的多边形单元，例如曲边三角形、曲边四边形等。图1-2二维问题的离散单元对于三维问题，昀基本的是四面体单元（图1-2d）或长方体和六面体（图1-2e，f），必要的时候也可以使用曲面体。在地下水模型中采用哪种形状的单元主要取决于计算区域的几何形状以及描述该问题所必须的独立空间坐标的数目。由于三角形单元比较灵活，能较好地适应复杂的边界形状，所以在研究二维问题的时候经常采用它。以下我们来重点讨论二维问题中的三角形单元剖分。假设存在平面区域D，将区域D剖分位移系列的三角形，这种剖分要一直进行到区域D的边界。当边界为直线时，边界的直线就作为三角形的一边；当边界为曲线时，就用适当的折线来逼近（图1-3）。各个三角形的顶点称为结点。总的说来，三角形的划分是任意的，以下是一些基本剖分原则：a)如果是用三角形单元，则三角形的昀长边和昀短边之比不能大于3；b)单元不能跨越非均质分区界线和内部边界界线；OYX边界e1e2ijme3e6e5e4e7e8e9e2e10e2e11e12e13图1-3区域D几何边界附近的三角形空间离散原则(a)(b)(d)(c)(e)(f)12312341234123412345768123456784c)单元的大小可根据具体情况调整。例如，在水力坡度变化比较大的地段，必须采用比较小的单元（单元密度大），反之可以减小单元密度；1.3水头插值根据资料的精度和研究需要，有限单元法中可以使用包括线形插值、二次插值、对数插值等各种插值方法。对于二维问题中的三角形单元，若使用线性插值，其水头插值函数形式如下：123(,)Hxyxyααα=++（1-1）对于三维问题中的四面体单元，若使用线性插值，其水头插值函数形式如下：1234(,,)zHxyzxyαααα=+++（1-2）1.4有限单元法的解题过程有限单元法的解题过程一般包括以下几个步骤1、建立研究区D的地下水数学模型。2、对研究区D进行空间离散和时间离散（对于非稳定流问题）。原则上是用有限多个足够小的单元来替代原有的研究空间。3、确定水头插值函数；水头插值函数直接影响计算精度。一般说来，非线性插值的精度要高于线性插值，但是这样处理势必产生较大的工作量，故此，在二维问题的三角形单元中，一般多采用如式1-1所示的多项式进行插值。4、整理模型所需要的数据。这些数据包括每个单元所包含的结点信息，每个结点的坐标，各类边界条件（包括边界上的结点号或单元号，结点对应的水头或单元对应的流量，降雨强度，参数分区资料等）。5、根据题目的性质和求解目标编制程序，求解未知水头结点的水头；这部分内容是本文的重点，有限元模拟中，利用了集合论的基本原理将等价于区域D中地下水运动定解问题的泛函离散为一系列单元泛函，再对单元泛函进行积分线性化处理并形成单元渗透矩阵，昀后将各单元渗透矩阵集合成区域D的总渗透矩阵A，得到水头方程组如下：[]{}{}=FAH（1-14）式中，A为总渗透矩阵；H为各离散点的待求水头；F为常数列向量，其中包含第一类边界结点的已知水头，第二类边界上的已知硫量，区域D内的降雨补给量，抽水井的开采量等各种垂直方向的水量交替量。6、利用计算得到的水头计算其它的物理量。6第二章里茨有限单元法的基本原理2.1变分原理2.1.1泛函和极值如图2-1所示，在二维空间X-Y中存在点A(xa,ya)和B(xb,yb)，求连接A，B的昀短曲线。从微分学角度出发，若将A，B两点间的曲线L离散成有限个足够小的曲线段，每段的长度为△L，则存在()()()()()22222211ababLxxyyxyyxxydxΔ≈−+−=Δ+ΔΔ⎛⎞=+Δ⎜⎟Δ⎝⎠′=+上式利用了两点间的距离公式，当△L足够小时，所产生的误差与△L为同价无穷小，因而我们可以进一步利用订积分概念，获得A，B两点间的曲线长度为()21bbaaLLydx′=Δ=+∫∫（2-1）显然，在确定A，B两点之后，L的大小完全取决于函数y(x)，不同的y(x)描述不同的曲线，因此，通过A，B两点间的所有曲线就集合成了一个曲线簇，我们将这个集合记为C0，也就是说C0代表了一个函数集合，这个函数集合称为（2-1）式的容许函数类。假如我们进一步取A(0,0)，和B(1,1)，则()1201Lydx′=+∫（2-2）由于(0)nyxn=通过原点(0,0),(1,1)两点，因此nyx=属于容许函数类C0。如果取n=1，y=x，则有1y′=从而得到()112001112Lydxdx′=+=+=∫∫如果取n=2，则y=x2；又因为()22yxx′′==A(x,y)aaB(x,y)bbabXOY△L图2-1两点间的昀短曲线7从而得到()()()12012012201121114ln241221.475Lydxxdxxxxx′=+=+⎡⎤=++++⎢⎥⎣⎦≈∫∫据此而知，对于容许函数类C0中的任一个函数()yyx=，曲线长L都有一个确定的值和它相对应。因此，L是函数的函数，在数学上称函数的函数为泛函。由于L的自变量是函数()yyx=，所以记为()LLyx=⎡⎤⎣⎦其中y=y(x)属于容许函数类C0。综上所述，我们可以得到一个自变量的泛函的一般表达式为()(),,baIyxFxyydx′=⎡⎤⎣⎦∫（2-3）如果在容许类函数C0中，，存在一个()yyx=，对于C0中邻近()yyx=的其他任何函数()yfx=，如果都有()()IfxIyx⎡⎤⎡⎤⎣⎦⎣⎦的话，则称为(2-3)的极小值曲线，或称泛函()Iyx⎡⎤⎣⎦在曲线()yyx=上取得它的极小值。2.1.2欧拉方程假如()yyx=是泛函（2-3）的极小值曲线，那()yyx=应该满足什么样的条件呢？我们来考虑()yyx=的一条邻近曲线()()()fxyxxαη=+。其中()xη为任意函数，它存在连续一阶导数。为了使()fx属于函数类0C，必须假设()()0abηη==，我们将这个假定代入（2-3）式，于是有：()()(),,,,bbaaIfxFxffdxFxyydxαηαη′′′==++⎡⎤⎣⎦∫∫（2-4）显然，如果()yyx=给泛函I以极值，则泛函（2-4）式在0α=时也必定存在极值。根据函数的极值条件，函数的导数在0α=时应该等于零，即00dIdαα==也就是存在8()()000,,bababyyadIdFxyydxdddfdfFfFfdxddFFdxααααηαηααααηη===′⎡⎤′′=++⎢⎥⎣⎦′⎡⎤⎛⎞′=+⎜⎟⎢⎥⎝⎠⎣⎦′=+∫∫∫（2-5）对（2-5）式中的第二项采用分布积分公式bbbbyaaaaFFdFFdxddxyydxyηηηη′∂∂∂′==−′′′∂∂∂∫∫∫这样一来，极值条件（4-5）变为00bbaabadIFdFdxdxdydxyFdFdxydxyαηηαη=∂∂=−′∂∂⎛⎞∂∂=−=⎜⎟′∂∂⎝⎠∫∫∫（2-6）由于η是满足()()0abηη==的任意函数，根据变分法的基本引理一1，可得0FdFydyα∂∂−=′∂∂（2-7）此为泛函（2-4）的极小值曲线()yyx=所必须满足的条件，也即泛函（2-4）在()yyx=取极值的必要条件。方程（2-7）称为欧拉（Euler）方程。根据（2-2）和（2-3）式，可知()()2,,1Fxyyy′′=+因为()()122211221Fyyyyy−′∂⎡⎤′′=+=⎣⎦′∂′+i由于（2-2）式中不含y项，因此有0Fy∂=∂所以，由（2-7）可得泛函（2-2）的欧拉方程为1变分法的基本引理一内容如下：设()fx是区间(),ab上的一个连续函数，如果函数()xη本身和它的导数在区间(),ab内是连续的，且在端点等于零，即()()0abηη==，如果对于任何这样的函数()xη恒使()()0baxfxdxη=∫成立，则()fx在区间(),ab上恒等于零，即()0fx≡9()201dydxy′−=′+（2-8）求（2-8）的积分得到()21yAy′=′+其中A为积分常数，将上式等号两端同时平方得到()2212yAy⎛⎞′⎜⎟=⎜⎟′+⎝⎠进一步得()()()2222221yAyAAy⎡⎤′′′=⋅+=+⋅⎣⎦昀后得到()2221yAA′⎡⎤−=⎣⎦或()2221AyA′=−另外，由于通过()0,Aay、()1,Bby两点之解实质上是上式在边界条件()()01yayyby=⎧⎪⎨=⎪⎩（2-9）处的解()100yyyyxaba−−=−−该式代表一条通过A、B两点的直线，这完全吻合几何学中的定律：过两点之间昀短的曲线是直线。尽管以上仅讨论了包含一个变量的泛函问题，但是不难看出，求泛函的极值曲线可以化为求微分方程（2-7）的解，即解微分方程的边值问题。2.2极值问题2.2.1多个自变