地下水流有限元模型

地下水数值模拟的有限元法的原理与有限差分法有所不同，有限差分法以偏微分方程的差分近似为基础，而有限元法以积分方程的离散近似为基础。在以往的研究中，推导有限元方程有两种方法，分别为迦辽金(Galerkin)法与里茨(Ritz)法。

Galerkin法采用加权余量的积分公式(以平面模型为例)为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:R(H')为近似解H'形成的余量函数;w为权函数;x为模型的某个子域。

Ritz法是利用与地下水偏微分方程等价的函数极限值函数进行求解，即

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:J(H)是等价范函，积分内公式F定义为使J(H)取极小值的函数H(x，y，t)恰好满足地下水流的偏微分方程。于是，地下水方程的近似解通过以下极值条件方程:

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

得到，式中Hp是模型子域上的节点水头。已有证明，当有限元网格单元采用相同的插值函数及形函数时，上述两种方法实际上是完全等价的，所形成的代数方程组相同。为此，本书中不再对有限元法的积分方程进行推导。

最常见的平面有限单元网格是三角形网格，而四边形网格也有较多的应用。三位有限元模型一般采用分层网格建立6节点或8节点等参单元。有限元的积分方程就是在网格单元的基础上近似求解的。

首先讨论平面三角形单元(图2-3)。单元e的3个顶点编号按照逆时针顺序分别为i，j，k，水头在单元内任意坐标点(x，y)的分布特征采用线性插值函数进行近似描述:

图2-3 三角形单元及其顶点编号

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:β1、β2、β3为几何系数。把节点i，j，k处的坐标和水头代入式(2-28)得

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

这是线性方程组，求解得

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:Ae为三角形单元面积。将式(2-33)～式(2-35)代入式(2-29)，则插值函数可改写为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

令:

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

是节点i、j、k在单元e内的形参数。对于固定坐标点(x，y)，形函数具有以下的性质:

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:上标e表示求和只能在同一三角形单元内进行。当所有节点的水头已知时，就可以把形参函数代入式(2-40)计算三角形单元内任意坐标位置的水头。

四边形单元也可以采用类似的方法建立插值函数和形函数，但是在形式上复杂一些。首先，任意的四边形都可以影射为一个正方形，图2-4a中(x，y)坐标系下的四边形1－2－3－4，可通过变换，形成2-4b中(η，ξ)坐标系下的正方形1－2－3－4。这个正方形单元被称为等参单元，其形函数定义为

图2-4 任意四边形单元的坐标变换

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:(ηi，ξi)为节点i的映射坐标。等参单元内的任意一点(η，ξ)，通过上述形参数与实际四边形单元内的对应坐标点建立如下关系，即

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:xi，yi为节点i的实际坐标。因此，等参单元内形函数与实际单元函数具有一一对应的关系:

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

而等参单元内的任意一点(η，ξ)的水头，也可以通过形函数变换为坐标下(x，y)处的水头

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

在建立有限元方程时，需要用到一下类型的偏导数

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

对于三角形网格中的一个节点p(图2-5)，建立有限元方程的积分区域包括了与p点相邻的所有节点p－1，p－2等，即以p点为共同顶点的所有三角形单元所覆盖的区域。积分采用分片法，即在相邻的三角形单元内分别积分，然后求总和，得到对应于p点的有限元方程。

图2-5 三角形网格中的节点及其水均衡控制区(阴影部分)

通过建立有限元网格节点的水均衡控制区，采用水量均衡原理，也可推导出地下水流的有限元方程，与积分法得到的结果相同。这种推导方法更加清楚地反映了地下水流的物理机制。如图2-5所示，节点p的水均衡控制区是一系列三角形重心与侧边中点的连线所围成的区域。每个三角形单元中，属于节点p水均衡控制区的面积是其单元面积的1/3。因此，节点p水均衡控制区的面积为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:e－n为第n个相邻三角形单元的整体编号;Ae－n为其面积。

地下水通过控制区周边流入控制区内的侧向流量QH，可以用流速在控制区边界线上的线积分得到

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:V为流速向量;n为边界线S上微分线段dS的内法线矢量;Vx和Vy分别为x方向和y方向的流速分量;m为含水层的厚度。设渗透系数的主轴与坐标轴一致，则流速向量可根据达西定律获取:

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

代入式(2-51)有

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

确定上述积分的方法是在节点p相邻的三角形单元内分别计算，然后求和，即

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

其中，每个三角形单元e－n流量贡献为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:S－n为在三角形单元e－n内的积分路径。以图2-5中的单元e－1为例，节点p控制区在该单元内的积分路径为A－B－C，其中A为点p和点p－1连线的中点，B为点p和点p－2连线的中点，C为单元e－1的重心。

根据三角形单元内水头的插值函数式(2-40)，有

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

根据式(2-41)得

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

同理有

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

可见水力梯度在三角形单元内是一常量，可以提取到式(255)的积分号外部，即

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

这样在线积分中只剩下积分路径在x轴和y轴上的投影，有

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

根据点A和点C与单元e－1的3个顶点之间的关系，有

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

将式(2-61)和式(2-61)代入式(2-60)，得

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

将式(2-63)代入式(2-59)，得

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

三角形单元e－1中的节点p、节点p－1以及节点p－2分别与图4-9中普遍单元的顶点i，j，k相互对应，有了式(2-64)之后，引入如下的系数来反映单元e内节点与节点之间的关系:

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:GpL为节点L与节点p的控制区之间的关系。利用式(2-65)可以把式(2-64)改写为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

代入式(2-64)得到流向节点p控制区的总侧向流量:

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:p－j和p－k分别表示沿着单元e－i从节点p逆时针移动遇到的第一个节点与第二个节点。

以图25中的节点p为例，展开式(2-67)，得到

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

在整个网格中，节点之间在各个单元中的关联系数可以合并为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:GpL为节点p与节点L之间在整体网络上的关联系数;EpL为共同拥有节点p和L的单元的数目，说明节点p和L之间没有直接关系。关联系数的重要特征是对称性:EpL=GpL。运用式(2-69)，流向节点p控制区的侧向流量可表示为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:N为网格中的总节点数。

对于承压含水层，在时间步长Δtk内，节点p控制区内水体积的增加量为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:k为时间步长。标准的有限元法是把差值函数引入上述积分，并分别在相邻单元内积分再求和，导出ΔVW中包含全部相邻节点的水头变化信息。目前，地下水流有限元模型中，普遍采用储量集中处理法，即用节点p的水头变化表示控制区的平均水头变化，有

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:Ap按照式(2-50)计算，这种简单的做法反而可以改善有限元的计算结果。根据水均衡原理，控制区水体积的增量与外部流量补给的水量相等，即

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:Qa是其他外部不计流量。利用式(2-70)和式(2-72)，可以展开式(2-73)为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

这就是承压含水层模型中对应节点p的隐式有限元方程。

潜水含水层的有限元方程可写为

典型煤矿地下水运动及污染数值模拟:Feflow及Modflow应用

式中:Sy为给水度。关联系数GpL中的导水系数TX，TY与当前时刻的带球水头有关，因此属于非线性方程。

按照引力加速度和惯性加速度的分离方法，航空重力测量分为航空重力测量和航空重力梯度测量（吴美平，张开东，2007）。

1航空重力测量

航空重力测量的基本原理：采用两个不同的加速度测量系统，其中一个系统的输出含有引力加速度，也就是测量的是比力；而另一个系统输出的是不含引力的加速度，于是在同一个坐标系内对两组加速度输出进行求差，即可消除共有的载体运动加速度。剩下的差值中包含引力加速度和传感器的系统误差。

根据牛顿第二定律，在惯性坐标系i中，质点的动力学方程为：

航空重力勘探理论方法及应用

式中： 为质点的惯性加速度；Gi为引力加速度；fi为比力，为重力传感器的观测值。

对于静态重力测量，此时惯性加速度 ，通过调整可使得重力仪保持水平，从而使重力仪敏感轴对准重力矢量方向，如果不考虑测量误差则重力仪的观测值就是重力加速度的大小。静态重力测量比较简单，精度可以达到很高；动态下的航空重力测量，情况就要复杂得多。一方面由于航空重力测量惯性加速度 不再为零，此时重力传感器的观测量是惯性加速度和引力加速度的共同影响；另一方面重力传感器敏感轴的指向难以保持稳定。爱因斯坦广义相对论中的等效原理指出：在一个封闭系统内的观测者不能区分作用于它的力是引力还是它所在的系统正在作加速运动，也就是说惯性加速度所造成的“重量感”和牛顿万有引力的效应是完全一样的。在运动载体上，重力仪就是这个封闭系统内的观测者，它感受引力和惯性力的作用，但不能对二者进行区别。因此，原理上航空重力测量需要解决两个基本问题（孙中苗，2004）：①运动状态下如何保持重力传感器的稳定指向或精确地获得重力传感器的指向值？②如何分离出引力加速度和惯性加速度？

航空重力测量指的是在飞机等运动平台上对地球重力加速度的测量，按照测量对象的不同，可分为航空重力矢量测量和标量测量。航空重力矢量测量又称比力测量，可同时获得重力加速度矢量g沿三个垂直方向上的分量gx、gy、gz。航空重力标量测量指的是对重力加速度沿铅垂线方向分量的测量，航空标量重力测量只需要测量重力扰动矢量垂直分量的大小（重力异常）。目前航空标量重力测量已经进入实用阶段，通常航空重力测量指的就是标量重力测量。航空矢量重力测量需要进一步获取重力扰动矢量的两个水平分量。

对当地测量坐标系而言（图3-1-1），原点0在仪器传感器中心，X沿子午圈方向指向北（N）,Y沿卯酉圈方向指向东（E）,Z沿椭球内法线方向指向地心（Down）。假设n为重力等位面的外法矢方向（即垂直向上的方向），则航空重力标量测量的值g为：

航空重力勘探理论方法及应用

航空重力标量测量的重力g值实质上是重力位沿重力梯度方向的空间变化率，或重力位沿向下的铅垂线方向的梯度值gz，也就是重力加速度沿铅垂线方向的分量，它与距离的平方成反比。

2航空重力梯度测量

重力梯度测量的基本原理：将两个共基线的加速度计的输出求差，以共模方式消除载体运动加速度 的影响；如果共用基线是旋转稳定的，则由该差值观测量可得到重力梯度分量。航空重力梯度测量对重力异常的高频分量非常敏感，因此在资源勘探等领域具有广阔的应用前景。

航空重力梯度测量指的是对重力加速度沿铅垂线的垂向分量gz对相应坐标轴方向的一次导数或梯度值Wzx、Wzy、Wzz的测量（图3-1-1）（张贵宾等，2007）。因此，航空重力梯度测量的是航空重力梯度张量的部分分量，即全部或部分航空重力梯度分量gzx、gzy、gzz。由于重力梯度测量对加速度计的精度要求非常高，其技术难度很大，所以对航空重力梯度测量不再作进一步展开。

图3-1-1 重力场矢量及重力梯度张量分量示意图