地下水溶质运移理论及模型读书报告

《地下水溶质运移理论及模型》读书报告

0、前言

本书作者陈崇希教授，浙江温州人士，生于1933年10月，1956年毕业于北京地质学院水文地质及工程地质专业。教授、博士生导师。现任中国地质大学环境地质研究所所长。《水文地质、工程地质》、《勘察科学技术》与《地球科学》杂志编委。湖北省地质学会名誉理事。长期从事地下水渗流理论，地下水数值模拟技术，地下水资源评价与管理，地质环境保护及地质灾害防治等方面的教学与科研工作；主持过“五五”、“八五”、“九五”国家科技攻关，国家自然科学基金等各类的科研项目31项；以第一作者身份获奖的有：省部级科技进步三等奖2项，二等奖3项，部优秀教材2等奖1项，国家科技进步三等奖1项。独著或以第一作者合作的著作有《地下水动力学》、《地下水流动问题数值方法》、《地下水溶质运移理论及模型》、《地下水混合井流的理论及应用》和《地下水不稳定井流计算方法》等七部，发表的论文有《用数值一解析法预测毛里塔尼亚伊迪尼水源地地下水开采动态》、《滨海多含水层系统地下水开采――水环境系统若干问题》等40余篇。1997年获地质矿产部“八五”科技工作有突出贡献先进个人。1992年获政府特殊津贴。培养硕士生18名，博士生13名。是中国地质大学211工程建设《地质环境保护及地质灾害防治》学科群的首席科学家。

该书是陈教授在“多孔介质水动力弥散理论及水质模型”的讲稿基础上修正补充而成，分7章进行叙述。《地下水溶质运移理论及模型》书中主要分三个方向：水动力弥散（微分）方程、水动力弥散（微分）方程的解及水动力弥散系数的计算。通过对本书的阅读，现就关于水动力弥散方程的解的做一下简单介绍。

水动力弥散方程的解法主要分为水动力弥散方程的解析解法和数值解法。本书第四章，重点介绍以基本解为基础，给出一、二、三维的解析解；第六章介绍了数值解，主要是利用有限差分和有限元法来求一、二维的解，并给出一些修

正过的数值解法。

1、 对水动力弥散方程的解析解法的理解 尽管解析解法在求解复杂的水动力弥散方程定解中存在一定缺陷，但仍然不可忽略它所起的作用。室内或野外试验都要根据解析解的实用条件来进行设计，并用解析解去拟合观测资料以求得水动力弥散系数。解析解中将瞬时注入点源问题的解称为基本解。由基本解出发，利用叠加原理导出线源、面源、多点源及连续注入问题的解。因此，点源问题的解是一切解的根本，需十分重视。

（1）空间瞬时点源的解

其基本条件是：①均质各向同性介质；②静止流场0=u ，弥散系数为常数，流体密度为常数（ρ=常数）；③0=t 时，在原点处瞬时注入溶质的质量为m 。 以瞬时点源的位置为原点，可以得出浓度C 是相对于原点对称的。可简化出纯弥散方程：

)(222222z

C y C x C

D t C ??+??+??=?? 式中，D 代表多孔介质的分子扩散系数。该式可看出，是球对称的，有利于纯弥散方式的应用讨论。

取半径为R 和R+dR 的两个球面所构成的单元体为均衡段，根据质量均衡有：

t

C V J n W J n W V dR R

D R D ??=??-??+ 式中，W 为球面积；n 为有效孔隙率；J D 为弥散通量，且R C D

J D ??-=,V v 为均衡段空隙体积。

忽略高阶微量，化简后得：

)(122R C R R R

D t C ??????=?? 于是该点源的定解问题可以写成：

??? ?

??????=??R C R R R D t C 22 (R ≧0,t>0) 0),(0==t t R C (R>0)

0),(=∞→R t R C (t>0)

0),(0==R t R C (t>0) m dR R n C =???∞

024π (t>0) （该式将点源处浓度

限制在有限区域）

通过Boltzmann 变换，将原来的偏微分方程定解问题转变为常微分方程定解问题，可求得空间瞬时点源的解为：

Dt R e Dt n m

t R C 423

2)(8),(-=π

从上式可得出：①等浓度面为圆心位于原点处的球面；②任何时候的浓度最大值都在原点处，且随着时间的增加，原点处的浓度减小。

⑵空间瞬时无限线源解

空间瞬时无限线源的作用可看着点源的连续分布，因考虑到点源基本解的微分方程是线性的，故采用叠加的方法，即积分法，可得空间无限线源的基本解为：

Dt r e nDt m t r C 412

4),(-=π 从上式可看出，浓度C 与z 无关，即在z 方向不产生弥散问题。也就是说我们可以将空间上的无限线源弥散问题转化成xy 平面上的二维弥散问题。于是，该解也可为平面瞬时点源问题的基本解。

⑶空间瞬时无限面源的解

根据点、线、面的构成原理，同理，可将空间无限面源看成是无数连续排列的无限线源组成，通过对无限线源的积分，可以得出空间无限面源的基本解为：

DT

z f e Dt n m t x C 42

2),(-=

π

从上式可看出：y 与z 无关，也就是说上述定解问题实质上是一维弥散问题。 以上解都是没有边界限制的，若加上边界，便成了有限空间问题。若边界简单，则可利用类似于水流问题中的反映法，将其变成无界问题，然后再采用叠加方法求出所需求的解。

（4）一维稳定流下水动力弥散问题的解

本章中水动力弥散问题都是在一维稳定流情况下讨论的，分为一、二、三维水动力弥散问题的解。

I 、一维水动力弥散问题与一维瞬时点源问题相近，初始条件与边界条件都相同。只是在示踪剂瞬时注入时，设其原有溶液浓度00=C ,并有速度u 稳定流动，求浓度),(t x C 的分布，从而造成一维水动力弥散问题比之多了一个对流项。本书中，采取坐标转换（按ut x X -=），利用一维瞬时点源问题的解，消去对流项，令,X x ut T t =-=。将新变量X 、T 反变换后得到：t D ut x L L e t D n m t x C 4)(22),(--=πω 当一维水动力弥散问题里初始浓度成阶梯状分布，即形成一维稳定流动一维水动力弥散问题，其数学模型可写成：

x C u x

C D t C L ??-??=??22

0(,)C t C -∞= 0>t

1(,)C t C ∞= 0>t

我们可以通过利用点源的基本解进行积分，再令ut x X -=，用换元法对它进行简化，得出：t D ut x L L e t D n m t x C 4)(22),(--=πω

而在半无限，一端为定浓度边界的限定情况下，一维水动力弥散问题的数学模型为：

x C u x

C D t C L ??-??=??22 0)0,(=x C 0≥x

0),0(C t C = 0>t

0),(=∞t C 0>t

该模型通过Laplace 变换，并利用边界条件，换元法可得出该定解问题的

解：01122L kx

D C erfc e erfc C =+，当 x 足够大或t 足够长时，该式为)2(210t

D ut x erfc C C L -=。 II 、二维水动力弥散问题中，注入平面瞬时点源时，同样可利用平面瞬时点源的基本解，通过换元等一系列转化、积分求得所求之解。只是必须清楚该问题在假定条件上有新的变化：①0≠u ，为一定值，流体非静止②水动力弥散系数为各向异性。通过一定关系的转化，得出该问题的解：???

?????+--=t D y t D ut x T L M T L e t D D M

m t y x C 44)(224),,(π。

当注入平面连续点源时，可将连续点源的作用看为无数瞬时点源之和，通过叠加原理，积分求得解：

220(,,)2()(,)

4L xu D L u t C x y t K W D ββ??=-????

式中，β=0()K β为第二类零阶修正贝赛尔函数，2(,)4L

u t W D β为第一类越流系统井函数。 当时间足够较长时，上式可简化为：

)44(2),(22222021

T L L D xu T L D D y u D x u K e D D Mn m y x C L +=π，此式也是计算水动力弥散系数常用的公式之一。

对于拟稳定条件下示踪剂的径向弥散，通常以井为中心，通过达西定律求出其平均速度，在极坐标下建立二维弥散方程，并利用复变函数理论求出其精确解，贝尔给出该定解问题的近似解为：

)342(21),(3

20r At r erfc C t r C L ??-=α（常用于确定实测的纵向弥散度L α）

III 、三维水动力弥散问题中，对于空间瞬时点源，其弥散系数D 是各向异性，且属于二度各向异性，若要利用前面基本解的结果（各向同性），就需进行相应的坐标变换，得到该问题的解为：

()22244313222(,,,)8L T x ut y z D t D t L T m

C x y z t e n

D D t π??-+??--????=

可看出:三维水动力弥散问题中，浓度C 于时间t 的二分之三次方成反比。对比

一、二维情况，不难看出，随着弥散维数的增加，浓度C 的衰减度也加快。根据浓度空间分布的时间函数，等浓度面是一个旋转的椭球面，其长轴沿x 方向。 注入空间连续点源时，假定注入的是理想示踪剂。将连续点源视为无数的瞬时点源之和，直接利用空间瞬时点源的解，利用积分得出解。当时间足够长时，该问题的解为：()12(,,)4u R x m C x y z e n DR π-=

2、 对水动力弥散方程的数值解法的理解 数值解法可以应用于复杂的情况，在实际应用中起着很好的效果。现就书上出现的数值解法做一些简单介绍。

（1）有限差分法

有限差分法的基本思想是：将研究空间划分成许多小的网格，把时间分成许多小段t ?，每个网格中心点处的未知变量视为该网格上的平均值，然后利用差商近似代替微商，形成研究区域上离散分布的有限个代数方程，求解方程组便可得该t ?时刻上各格点上的取值。然后按照一个个的t ?逐个往前求解。

其基本原理：就是将某点处的浓度函数的导数用该点处和其n 个相邻点处的浓度值及其间距近似表示，常通过泰勒展开式建议浓度导数的近似。

其求解地下水动力弥散问题的基本步骤：

①剖分渗流区，确定离散点；

②建立水的动力弥散问题的差分方程；

③求解差分方程。

I、一维水动力弥散的差分解法

根据一维水动力弥散方程，根据差分原理，采用向前、向后和中心差分3种不同的差分格式进行差分，可分别得到显式、隐式和Crank-Nicolson等3种不同的差分格式。

显式是有稳定条件要求的，而隐式和Crank-Nicolson式都是无条件稳定的，都可以利用“追赶法”解三对角方程组求

t时候的浓度值

1

n

II、二维水动力弥散的差分解法

以一维流动二维水动力弥散方程为依据，用差分近似代替方程中的偏导数，同理得到二维水动力弥散方程的显式、隐式和Crank-Nicolson 3种不同的差分格式。III、一、二维水动力弥散的差分解法的比较

相同点：都采用类似的差分原理进行差分，得到的差分格式的基本类型一致。区别：一维条件下3种格式采用“追赶法”求解，的三对角线方程组。而二维条件下所给出的三种格式组成的方程组是五对角线方程组。为了避免解五对角线方程组的困难，特提出交替方向隐式法，简称ADI方法。它的优点是：不是一次对整体矩阵求逆，而是分两次对三对角线矩阵求逆，这样就把二维问题简化为多次解一维问题。

⑵有限单元法

与有限差分法相同，有限单元法也是根据区域剖分和插值方法将水动力弥散的定解问题化为代数方程组进行求解的。本书主要介绍一维稳定流场中二维水动力弥散问题伽辽金有限单元法，伽辽金有限单元法在伽辽金法的基础上发展而来，是伽辽金法与有限单元法的结合。

I、迦辽金法是寻找一个级数形式的试探函数作为微分方程的近似解，并使其满足给定的边界条件，并令剩余的加权积分为零，经过积分，最终得出伽辽金方程：

Γ=??+??+????+???????Γd N J dxdy N t C N x C u y N y C D x N x C D L D D L L L T L L

)(~

~~~ n L ,...,3,2,1=

迦辽金法属于加权剩余法，且由于其他加权剩余法，应用更普遍。 II 、有限单元法

用有限单元法求解地下水水动力弥散问题，首先要将区域单元均质剖分，常用的单元形状有三角形单元和矩形单元。本书中介绍的矩形有限元。在单元剖分后，随即构造单元基函数N ’L ，且矩形单元满足：在L 结点处，N L =1；在其他结点处，N L =0.基函数N L 在各个矩形单元中都按双线性规律变化。

III 、伽辽金有限单元法

将迦辽金方程与有限元剖分思想结合起来，就建立了矩形单元的迦辽金有限元方程。最终得出：

①当L 为内结点，则N L =0，故F L =0;

②当L 为第二类结点时，,,1122

L D iL D Lj F iLJ LjJ =+。 式中，,D iL J 为边界iL 段上的弥散通量；,D Lj J 为边界Lj 段上是弥散通量。

（3）修正方法

在求解弥散方程中常出现过量和数值弥散，为克服它们引起的误差，提高数值解法的精确性和稳定性，本书特提出几种修正的数值方法（上游加权法、特征值方法、动坐标系方法方法与网格变形方法、随机步行法、引入人工扩散量的方法），这里将不再累赘。

3、结论

本书中提到的水动力弥散方程的解法都是现目前在地下水溶质运移里常用的两种方法，解析解法和数值解法都有着它们各自的优点，应合理利用它们的优势，力争求得与实际更为接近的值，为之后的研究提供精确有力的依据。当然，我们的专业研究人员更应该致力于将这些弥散方程的解法发展拓宽，与时俱进地与计算机技术更好地融合在一起，以便更方便更快捷地得出结果，更有效率地促

进水文地质学的发展。