溶质运移理论水动力弥散方程的值解法-课件

地下水流数值模拟方法第六章

水动力弥散方程的数值解法中国地质大学环境学院2019春地下水流数值模拟方法第六章水动力弥散方程的数值解法中国地质12一、有限差分法基本思想：

将空间划分成若干网格，把时间分成若干小时段，每一个网格中心点（格点）处的未知变量（H或C）视为该网格平均值。利用差商代替微商。基本步骤：

（1）剖分（2）建立差分方程组（3）求解2一、有限差分法基本思想：23一、有限差分法-导数的有限差分近似图中，去x轴上任意一点i，其坐标为

在改点左右相聚为处分别取（i-1）和（i+1），其坐标分别为和以i为中心，泰勒展开C(x)3一、有限差分法-导数的有限差分近似图中，去x轴上任意一点i34一、有限差分法-导数的有限差分近似整理并略去余项（6-1）-（6-2），再除以略去余项（6-1）-（6-2），再除以略去余项：分别一阶导数的向前差分、向后差分及中心差分二阶中心差分4一、有限差分法-导数的有限差分近似整理并略去余项分别一阶导45一、有限差分法-一维水动力弥散的差分解法一维水动力弥散方程

（6-7）（1）显格式式中仅有一个未知数，解得

式（6-7）中的对流项取中心差分5一、有限差分法-一维水动力弥散的差分解法一维水动力弥散方程56

可以证明，稳定性准则要求

即

即。（1）（2）。由条件（2），格距要求很小；由（2）可知，鉴于较小，导致不能太大，将（2）代入（1）式中，得到

6

可以证明，稳定性准则要求

即67若对流项改为向后差分

解得

稳定性要求不难看出，稳定性限制比对流项取中心差分有所改善。7若对流项改为向后差分解得78（2）隐式格式

整理后得到

隐格式是无条件稳定的。8（2）隐式格式89（2）Grank-Nicolson格式

整理后得到

取隐式和显示的平均，即Grank-Nicolson格式9（2）Grank-Nicolson格式910一、有限差分法-二维水动力弥散的差分解法二维水动力弥散方程

（4-56）（1）显格式式中仅有一个未知数式（4-56）中的对流项取中心差分10一、有限差分法-二维水动力弥散的差分解法二维水动力弥散方1011一、有限差分法-二维水动力弥散的差分解法化简后，有涉及以（i，j）为中心的5个网格点在tn时刻的已知浓度11一、有限差分法-二维水动力弥散的差分解法化简后，有涉及以1112一、有限差分法-二维水动力弥散的差分解法（2）隐格式式（4-56）中右端的对流项取中心差分，右端个C的时阶均取n+1水平12一、有限差分法-二维水动力弥散的差分解法（2）隐格式式（1213一、有限差分法-二维水动力弥散的差分解法（2）隐格式整理后收敛且无条件稳定涉及以（i，j）为中心的5个网格点在tn+1时刻的未知浓度13一、有限差分法-二维水动力弥散的差分解法（2）隐格式整理1314一、有限差分法-二维水动力弥散的差分解法（3）Grank-Nicolson格式将隐式格式的两式相加除以214一、有限差分法-二维水动力弥散的差分解法（3）Grank1415一、有限差分法-二维水动力弥散的差分解法（3）Grank-Nicolson格式整理得15一、有限差分法-二维水动力弥散的差分解法（3）Grank1516一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐式法（ADI法）分两次对三对角矩阵求逆，将一个Δt分成两个Δt/2计算第一个半时间步，对x方向的偏导数采用隐式差分，对y方向的偏导数采用显示差分。16一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐1617一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐式法（ADI法）整理得第二个半时间步，对y方向的偏导数采用隐式差分，对x方向的偏导数采用显示差分。17一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐1718一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐式法（ADI法）整理得收敛且无条件稳定18一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐1819一、有限差分法-求解差分方程的计算机程序举例算例

见教材P81-8519一、有限差分法-求解差分方程的计算机程序举例算例1920二、有限单元法-伽辽金有限单元法伽辽金法

对均匀多孔介质，一维稳定流场中二维水动力弥散方程

取x轴方向与地下水流向相同，记伽辽金法是寻找一个级数形式的试探函数作微分方程的近似解，并使其满足边界条件

（6-26）20二、有限单元法-伽辽金有限单元法伽辽金法（6-26）2021二、有限单元法-伽辽金有限单元法上式一般不会满足方程，因为仅是近似解，得到一个剩余误差函数R(x，y)，在平均意义下迫使误差为0

我们加以权，令剩余的加权积分为0，W是一组权函数。加权剩余法（6-29）21二、有限单元法-伽辽金有限单元法上式一般不会满足方程，因2122二、有限单元法-伽辽金有限单元法在整个研究区D上，基函数NL(x,y)分段定义(1)函数区域连续（2）一阶导数不连续（3）二阶导不容易确定故采用分步积分的方法使二阶导数降阶22二、有限单元法-伽辽金有限单元法在整个研究区D上，基函数2223二、有限单元法-伽辽金有限单元法对一维积分整理后分步积分推广到二维的情况（6-32）代入（6-29）23二、有限单元法-伽辽金有限单元法对一维积分（6-32）代2324二、有限单元法-伽辽金有限单元法有限单元剖分与基函数（1）三角形单元见《地下水流动问题数值模拟》（2）矩形单元将区域记为D，边界记为B。要求各单元均质、等厚，即T、μ为常数。结点（内结点、边界结点）

24二、有限单元法-伽辽金有限单元法有限单元剖分与基函数2425二、有限单元法-伽辽金有限单元法构造单元函数NL基函数取为“双线性插值”

基函数NL（x，y）形状如同一顶高等于1、有4条直线斜边和4条下凹型曲边的尖顶斗笠25二、有限单元法-伽辽金有限单元法构造单元函数NL基函数N2526二、有限单元法-伽辽金有限单元法子区DL以结点L为其共同结点的所有矩形单元（<=4）基函数NL仅在子区上不为0，在非DL上均为0,属于子区DL单元e的单元基函数NL，用NeL表示（>0）26二、有限单元法-伽辽金有限单元法子区DL2627二、有限单元法-伽辽金有限单元法典型矩形单元该单元中心位于坐标原点处，且4条边分别平行x，y轴,结点从左下角开始按逆时针方向编号i,j,m,n。沿x方向长2a，y方向宽2b。27二、有限单元法-伽辽金有限单元法典型矩形单元2728二、有限单元法-伽辽金有限单元法按上述要求所构造的NeL形式为对典型单元，令

则28二、有限单元法-伽辽金有限单元法按上述要求所构造的NeL2829二、有限单元法-伽辽金有限单元法其中即（6-34）29二、有限单元法-伽辽金有限单元法其中即（6-34）2930二、有限单元法-伽辽金有限单元法Nei(x,y)在结点i处为1，其它3处为0根据上式，有平行x或y方向上，Nei(x,y)线性变化故单元基函数为双线性插值函数若结点坐标用表示结点浓度用表示30二、有限单元法-伽辽金有限单元法Nei(x,y)在结点i3031二、有限单元法-伽辽金有限单元法结点浓度CL(t)和单元基函数NeL(x,y)来确定单元内的近似解，即对于区域D，结点L的基函数在子区DL内各单元分块确定。称为矩阵单元e上的基函数（6-36）31二、有限单元法-伽辽金有限单元法结点浓度CL(t)和单元3132二、有限单元法-伽辽金有限单元法伽辽金有限单元法

左端积分范围为D，由于在上NL=0,改写在子区DL的积分，对矩形单元逐个积分、求和。设子区DL上有neL个单元，有（6-39）32二、有限单元法-伽辽金有限单元法伽辽金有限单元法（6-33233二、有限单元法-伽辽金有限单元法由（6-36）得由（6-34）得再由（6-36）得（6-40）（6-42）33二、有限单元法-伽辽金有限单元法由（6-36）得由（6-3334二、有限单元法-伽辽金有限单元法将（6-40）（6-42）代入（6-39）得34二、有限单元法-伽辽金有限单元法将（6-40）（6-43435二、有限单元法-伽辽金有限单元法（6-44）35二、有限单元法-伽辽金有限单元法（6-44）3536二、有限单元法-伽辽金有限单元法用矩阵形式表示一阶微分方程组将一阶导数用差分形式表示按矩阵符号有若﹛C﹜在新时间水平t+Δt上取值，则36二、有限单元法-伽辽金有限单元法用矩阵形式表示一阶微分方3637二、有限单元法-伽辽金有限单元法（1）给定初始条件和边界条件，可求（6-46）各个系数矩阵，于是可解得第一个末时刻各结点的浓度；（2）以此浓度为已知浓度，计算新的系数矩阵；（3）再次求解，得到第二个时间步长末时刻结点浓度（6-46）浓度可在之间任意位置近似表示37二、有限单元法-伽辽金有限单元法（1）给定初始条件和边界3738二、有限单元法-伽辽金有限单元法按单元顺序计算系数矩阵每个矩阵单元涉及4个结点i，j，m，n，涉及4个方程单元e的作用，可用4行与4列的矩阵来表示，称为单元矩阵38二、有限单元法-伽辽金有限单元法按单元顺序计算系数矩阵3839二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体弥散矩阵的关系设研究区划分为8个单元，16个结点，四周为隔水边界编号如下：39二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体3940二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体弥散矩阵的关系计算每个单元弥散矩阵，按双下标已知，放置单元弥散矩阵并叠加，得到总体弥散矩阵。40二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体4041二、有限单元法-伽辽金有限单元法从（6-44）知将坐标原点移到矩形单元中心有41二、有限单元法-伽辽金有限单元法从（6-44）知将坐标原4142二、有限单元法-伽辽金有限单元法通过双下标一致的元素求和，得到总体弥散矩阵诸元素42二、有限单元法-伽辽金有限单元法通过双下标一致的元素求和4243二、有限单元法-伽辽金有限单元法积分求解可用高斯求积方法，详见P97-10543二、有限单元法-伽辽金有限单元法积分求解可用高斯求积方法4344二、有限单元法-伽辽金有限单元法44二、有限单元法-伽辽金有限单元法4445二、有限单元法-伽辽金有限单元法算例：二维水动力弥散的伽辽金法计算机程序

数学模型写成剖分成9个单元共20个结点，空间步长取Δx=Δy=10m，时间步长Δt=5d，渗透流速u=0.01m/d具体代码见P106-11145二、有限单元法-伽辽金有限单元法算例：二维水动力弥散的伽4546三、过程与数值弥散过量现象

如一维流动水动力弥散中单位浓度梯度函数注入的情况。

对比解析解与数值解C-x曲线，在浓度封面附近数值计算的浓度超过最大浓度值1和小于最小浓度0，失去物理意义。数值弥散

若纯对流方程采用有限差分逼近，其结果呈现似有弥散的存在。46三、过程与数值弥散过量现象数值弥散4647三、过程与数值弥散47三、过程与数值弥散4748三、过程与数值弥散过量现象-解析

由于时间步长与空间尺度匹配不佳，因而含水层在数值上不能“吸收”住污染物所引起。解决办法：

将时间增量选择为几何级数的子项；

做试验校正时间步长和差分网格；48三、过程与数值弥散过量现象-解析4849三、过程与数值弥散数值弥散-解析

由截断误差引起。

对流项一阶空间导数采用向后差分逼近，浓度Ci-1泰勒展开变形得

向后差分的截断误差主部分为相当于弥散系数的弥散项49三、过程与数值弥散数值弥散-解析相当于弥散系数的弥散项4950三、过程与数值弥散数值弥散-解析

由故纯对流方程有限差分近似表示后，结果如同对流-弥散方程的结果，记用数值弥散系数与物理弥散系数的比值来评价数值弥散对物理弥散的干涉程度此条件下的数值弥散系数50三、过程与数值弥散数值弥散-解析此条件下的数值弥散系数5051三、过程与数值弥散数值弥散-解析

定义为Peclet数，时，弥散占优；

时，对流占优。速度u与物理弥散系数DL比越大，越大。对于小数，数值解与精确解非常接近，对于大数，数值解过渡带变宽，浓度锋面变平缓，在锋面的前后还出现解的振动—过量现象。控制Δx的大小从而减少数值弥散51三、过程与数值弥散数值弥散-解析控制Δx的大小从而减少数5152三、过程与数值弥散对纯对流方程

52三、过程与数值弥散对纯对流方程5253三、过程与数值弥散对纯对流方程

结果也造成数值弥散此时，总的数值弥散是两者之和53三、过程与数值弥散对纯对流方程结果也造成数值弥散此时，总5354四、对流占主导地位时的数值方法上游加权法

对流项的有限差分化过程中乘上一个适当的权因子，则波动和过量得到改善或消除一维水动力弥散中对流项取中心隐式差分若相邻格点i-1或（和）i+1的浓度越大，则格点i的浓度也越大，反映到差分方程上，系数Gi-1和Gi+1不得与Gi同号，上式要求

不能满足时，对流占优势54四、对流占主导地位时的数值方法上游加权法对流项的有限差分5455四、对流占主导地位时的数值方法上游加权法

假定DL和u为常数，取等格距差分网格。一维对流-弥散方程在积分，有对流-弥散通量用有限差分近似表示

55四、对流占主导地位时的数值方法上游加权法假定DL和u为常5556四、对流占主导地位时的数值方法上游加权法

ω为权因子；当ω=1/2时，对流项相当于中心差分；当ω=1时，对流项相当于取向后差分，

当ω=0时，对流项相当于向前差分。速度不同，

ω取值不同56四、对流占主导地位时的数值方法上游加权法5657四、对流占主导地位时的数值方法上游加权法

上游浓度权因子大于1/2。物理意义：加强上游格点浓度值在对流项中的作用可通过浓度的线性（或二次）插值具体确定在Δt期间C在格边上随时间的变化，从而确定ω值。57四、对流占主导地位时的数值方法上游加权法上游浓度权因子大5758四、对流占主导地位时的数值方法

仍用Gi-1、Gi和Gi+1分别表示上式左端第一、二、三项括号内系数。则采用上游加权法有限差分格式，可消除过量和波动，但截断误差增大，数值弥散反而增大。中心差分格式，截断误差为二阶上游加权有限差分截断误差为一阶58四、对流占主导地位时的数值方法仍用Gi-1、Gi5859四、对流占主导地位时的数值方法

上游权因子的物理实质：对流项浓度C在某断面，例如断面的中心差分意味着对流项的浓度在Δt时段内均以i与i-1两格点的平均浓度，即通过断面；实际为：Δt时段的初始时刻是该平均浓度值，随后在对流的作用下，该断面的浓度显然会向上游格点的浓度Ci-1(u>0时)的变化。59四、对流占主导地位时的数值方法上游权因子的物理实5960四、对流占主导地位时的数值方法特征值方法

令源汇项为0，展开对流-弥散方程

60四、对流占主导地位时的数值方法特征值方法6061四、对流占主导地位时的数值方法特征值方法

假定流体密度值ρ不变（不可压缩），含水层不可压缩，则，（6-75）简化为基本思路：物质导数表示一个沿轨迹（特征曲线）

运动的研究对象浓度变化率。由弥散作用引起。若存在运会想，还将受到沿轨迹源与汇、化学、生物化学反应、吸附与解析等。对流项用沿特征线的示踪质点的运动描述

（6-76）61四、对流占主导地位时的数值方法特征值方法（6-76）6162四、对流占主导地位时的数值方法特征值方法

特征点均匀分布研究区，赋一个初始浓度。

含水层中未受污染部分特征点的初始浓度为零，特征点浓度眼轨迹的瞬时浓度变化，借助矩形网格差分格式计算。62四、对流占主导地位时的数值方法特征值方法6263四、对流占主导地位时的数值方法63四、对流占主导地位时的数值方法6364四、对流占主导地位时的数值方法模型精确度取决于网格距大小。可以在同一网格上计算速度场。如果采用较少的网格距，可以改善流场的描述。

坐标为的特征点所在的单元，可用下标i,j确定

精确度受到每一网格内所选择的最初的特征点数目影响。64四、对流占主导地位时的数值方法模型精确度6465四、对流占主导地位时的数值方法缺点：（1）实际程序冗长，需记录特征点的踪迹；

（2）抽水点或出流边界的特征点需要消除，入渗边界须有特征点生成，不透水边界特征点须反射回来；

（3）要通过消去特征点来避免滞留点附近特征点堆积；（4）特征点耗尽的单元，必须产生新的特征点。65四、对流占主导地位时的数值方法缺点：6566四、对流占主导地位时的数值方法动坐标系方法与网格变形方法能消除过量现象，还保持较陡的浓度锋面。

对式

设坐标变换

得，

取空间步长，则

66四、对流占主导地位时的数值方法动坐标系方法与网格变形方法6667四、对流占主导地位时的数值方法动坐标系方法与网格变形方法

在固定坐标系中取格距为，有

动坐标系中格点运动与固定坐标系格点运动重合。

设固定坐标系中格点xj在tn+1时刻与动坐标系中Xi重合，有

动坐标系下浓度结点的表达式67四、对流占主导地位时的数值方法动坐标系方法与网格变形方法6768四、对流占主导地位时的数值方法动坐标系方法与网格变形方法

式变成

对半无限长砂柱，当DL=0时68四、对流占主导地位时的数值方法动坐标系方法与网格变形方法6869四、对流占主导地位时的数值方法网格变形方法

设tn时刻有限元结点在空间中的分布是对应浓度分布为

69四、对流占主导地位时的数值方法网格变形方法6970四、对流占主导地位时的数值方法网格变形方法

主要步骤：结点位置随时间变化，基函数及有限元方程的系数都依赖于时间，每推进一个Δt都需要重新计算。二维问题还须进行畸形判断

70四、对流占主导地位时的数值方法网格变形方法7071四、对流占主导地位时的数值方法结合动坐标的网格变形方法

一维对流-弥散在某个动点在t时刻位置x可写成

其中，x0是动点的初始位置，有取时间的一阶偏导数对流-弥散方程写成

若控制动点速度dx/dt并让它接近流速u，则可转换成弥散为主的的方程式71四、对流占主导地位时的数值方法结合动坐标的网格变形方法若7172四、对流占主导地位时的数值方法若使用伽辽金法，基函数因依赖于结点位置而成为时间的函数。若用表示随时间变化的基函数，则对任意t，有

故72四、对流占主导地位时的数值方法若使用伽辽金法，基函数因依7273四、对流占主导地位时的数值方法随机步行法

采用示踪剂描述污染物运移。只有被污染的特征点在流场中运动。每个特征点都给定一个不变的污染物质量，此质量的和等于排进含水层的总污染物质量。取许多单个特征点轨迹（随机步行）的平均值，可得到一个有弥散的特征点的分布。要得到一个浓度分布，就要先叠加一个网格并算每个网格单元所含有的特征点数。

73四、对流占主导地位时的数值方法随机步行法7374四、对流占主导地位时的数值方法随机步行法-以一维运移为例在x=0处，一个瞬时注入的、质量为Δ

M的理想失踪及，在t时刻浓度分布对一个固定时间t，看做正态分布74四、对流占主导地位时的数值方法随机步行法-以一维运移为例7475四、对流占主导地位时的数值方法随机步行法-以一维运移为例在时间t=0和位置x=0处，放置具有污染物质量为ΔM/N的N个特征点，每一个特征点移动距离x，而在时刻t到达它们各自位置，即Z是一个平均值为0，标准差为1的正态分布随机变量

所得的频率分布，通过一个标准化因子，有对有限数量的颗粒，可近似去顶C(x,t)，将空间变量划分为长度Δx的若干区间，区间中点浓度75四、对流占主导地位时的数值方法随机步行法-以一维运移为例7576四、对流占主导地位时的数值方法随机步行法-以一维运移为例若孔隙平均流速是时间和位置的函数，那么时间为0开始的颗粒P在时间间隔Δt内的随机轨迹xp(t)为

76四、对流占主导地位时的数值方法随机步行法-以一维运移为例7677四、对流占主导地位时的数值方法随机步行法-二维流场

首先假定流动方向平行于x轴，得到Z和Z’是正态分布的随机变量的两个值。对任意方向的流动，需考虑弥散张量特性。77四、对流占主导地位时的数值方法随机步行法-二维流场7778四、对流占主导地位时的数值方法随机步行法-二维流场

首先假定流动方向平行于x轴，得到Z和Z’是正态分布的随机变量的两个值。对任意方向的流动，需考虑弥散张量特性。对流运动可确定在水流流动方向及其正交方向的弥散运动78四、对流占主导地位时的数值方法随机步行法-二维流场对流运7879四、对流占主导地位时的数值方法随机步行法

79四、对流占主导地位时的数值方法随机步行法7980四、对流占主导地位时的数值方法随机步行法

数值结果的质量主要取决于特征点数目N，要考虑网格的选择。大网格间距会产生极平均的结果，小网格间距会产生十分粗糙的分布。

在稳定流动和一个源或数个源都有固定强度的情况下，可以节省很多特征点。80四、对流占主导地位时的数值方法随机步行法8081818182四、对流占主导地位时的数值方法随机步行法

nδ是在时间τ内，有在时间t=0时瞬时注入所产生的分布。所得的特征点分布加权叠加到最终的特征点分布上。用迟滞流速u/R代替平均孔隙流速u，并用迟滞系数去除全部注入质量，可引入线性吸附作用。分配给特征点一个按时间变化的污染物质量82四、对流占主导地位时的数值方法随机步行法8283四、对流占主导地位时的数值方法随机步行法-特点

不会出现数值弥散，但在污染带边缘的得到满意的计算结果需要大量的特征点。

完整的，易于变成并能用于每一个水流模型的方法，易推广到三维模型。

灵敏度较小的参数的变化引起的浓度分布变化会被随机性浓度的变化覆盖。

83四、对流占主导地位时的数值方法随机步行法-特点8384四、对流占主导地位时的数值方法引入人工扩散量方法

对流-弥散方程简写成简写成脚标注释的形式引入人工扩散系数并赋权得

（6-107）84四、对流占主导地位时的数值方法引入人工扩散量方法（6-18485四、对流占主导地位时的数值方法引入人工扩散量方法

当θu=0时，导出对称的正定矩阵，若不校正，格式是不稳定的。在时间步长中点泰勒展开

（6-109）85四、对流占主导地位时的数值方法引入人工扩散量方法（6-18586四、对流占主导地位时的数值方法引入人工扩散量方法

（6-107）对时间求偏导并代入（6-109）得在时，是具有二阶截断误差的Grank-Nicolson格式

86四、对流占主导地位时的数值方法引入人工扩散量方法8687四、对流占主导地位时的数值方法引入人工扩散量方法

为评价对流项引起的误差，将上式各项归并后得当人工扩散项与弥散权重

采用如下形式，对流项取非中心权而引起的误差可消除：87四、对流占主导地位时的数值方法引入人工扩散量方法8788四、对流占主导地位时的数值方法引入人工扩散量方法

只有当时，方可导出对称正定系数矩阵，上两式可写成

（6-116）（6-117）88四、对流占主导地位时的数值方法引入人工扩散量方法（6-18889四、对流占主导地位时的数值方法引入人工扩散量方法

（6-116）所定义的人工扩散是一种有效校正。

（6-117）弥散项取隐式对数值弥散没有影响，却提高求解精度。89四、对流占主导地位时的数值方法引入人工扩散量方法89地下水流数值模拟方法第六章

水动力弥散方程的数值解法中国地质大学环境学院2019春地下水流数值模拟方法第六章水动力弥散方程的数值解法中国地质9091一、有限差分法基本思想：

将空间划分成若干网格，把时间分成若干小时段，每一个网格中心点（格点）处的未知变量（H或C）视为该网格平均值。利用差商代替微商。基本步骤：

（1）剖分（2）建立差分方程组（3）求解2一、有限差分法基本思想：9192一、有限差分法-导数的有限差分近似图中，去x轴上任意一点i，其坐标为

在改点左右相聚为处分别取（i-1）和（i+1），其坐标分别为和以i为中心，泰勒展开C(x)3一、有限差分法-导数的有限差分近似图中，去x轴上任意一点i9293一、有限差分法-导数的有限差分近似整理并略去余项（6-1）-（6-2），再除以略去余项（6-1）-（6-2），再除以略去余项：分别一阶导数的向前差分、向后差分及中心差分二阶中心差分4一、有限差分法-导数的有限差分近似整理并略去余项分别一阶导9394一、有限差分法-一维水动力弥散的差分解法一维水动力弥散方程

（6-7）（1）显格式式中仅有一个未知数，解得

式（6-7）中的对流项取中心差分5一、有限差分法-一维水动力弥散的差分解法一维水动力弥散方程9495

可以证明，稳定性准则要求

即

即。（1）（2）。由条件（2），格距要求很小；由（2）可知，鉴于较小，导致不能太大，将（2）代入（1）式中，得到

6

可以证明，稳定性准则要求

即9596若对流项改为向后差分

解得

稳定性要求不难看出，稳定性限制比对流项取中心差分有所改善。7若对流项改为向后差分解得9697（2）隐式格式

整理后得到

隐格式是无条件稳定的。8（2）隐式格式9798（2）Grank-Nicolson格式

整理后得到

取隐式和显示的平均，即Grank-Nicolson格式9（2）Grank-Nicolson格式9899一、有限差分法-二维水动力弥散的差分解法二维水动力弥散方程

（4-56）（1）显格式式中仅有一个未知数式（4-56）中的对流项取中心差分10一、有限差分法-二维水动力弥散的差分解法二维水动力弥散方99100一、有限差分法-二维水动力弥散的差分解法化简后，有涉及以（i，j）为中心的5个网格点在tn时刻的已知浓度11一、有限差分法-二维水动力弥散的差分解法化简后，有涉及以100101一、有限差分法-二维水动力弥散的差分解法（2）隐格式式（4-56）中右端的对流项取中心差分，右端个C的时阶均取n+1水平12一、有限差分法-二维水动力弥散的差分解法（2）隐格式式（101102一、有限差分法-二维水动力弥散的差分解法（2）隐格式整理后收敛且无条件稳定涉及以（i，j）为中心的5个网格点在tn+1时刻的未知浓度13一、有限差分法-二维水动力弥散的差分解法（2）隐格式整理102103一、有限差分法-二维水动力弥散的差分解法（3）Grank-Nicolson格式将隐式格式的两式相加除以214一、有限差分法-二维水动力弥散的差分解法（3）Grank103104一、有限差分法-二维水动力弥散的差分解法（3）Grank-Nicolson格式整理得15一、有限差分法-二维水动力弥散的差分解法（3）Grank104105一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐式法（ADI法）分两次对三对角矩阵求逆，将一个Δt分成两个Δt/2计算第一个半时间步，对x方向的偏导数采用隐式差分，对y方向的偏导数采用显示差分。16一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐105106一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐式法（ADI法）整理得第二个半时间步，对y方向的偏导数采用隐式差分，对x方向的偏导数采用显示差分。17一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐106107一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐式法（ADI法）整理得收敛且无条件稳定18一、有限差分法-二维水动力弥散的差分解法（4）交替方向隐107108一、有限差分法-求解差分方程的计算机程序举例算例

见教材P81-8519一、有限差分法-求解差分方程的计算机程序举例算例108109二、有限单元法-伽辽金有限单元法伽辽金法

对均匀多孔介质，一维稳定流场中二维水动力弥散方程

取x轴方向与地下水流向相同，记伽辽金法是寻找一个级数形式的试探函数作微分方程的近似解，并使其满足边界条件

（6-26）20二、有限单元法-伽辽金有限单元法伽辽金法（6-26）109110二、有限单元法-伽辽金有限单元法上式一般不会满足方程，因为仅是近似解，得到一个剩余误差函数R(x，y)，在平均意义下迫使误差为0

我们加以权，令剩余的加权积分为0，W是一组权函数。加权剩余法（6-29）21二、有限单元法-伽辽金有限单元法上式一般不会满足方程，因110111二、有限单元法-伽辽金有限单元法在整个研究区D上，基函数NL(x,y)分段定义(1)函数区域连续（2）一阶导数不连续（3）二阶导不容易确定故采用分步积分的方法使二阶导数降阶22二、有限单元法-伽辽金有限单元法在整个研究区D上，基函数111112二、有限单元法-伽辽金有限单元法对一维积分整理后分步积分推广到二维的情况（6-32）代入（6-29）23二、有限单元法-伽辽金有限单元法对一维积分（6-32）代112113二、有限单元法-伽辽金有限单元法有限单元剖分与基函数（1）三角形单元见《地下水流动问题数值模拟》（2）矩形单元将区域记为D，边界记为B。要求各单元均质、等厚，即T、μ为常数。结点（内结点、边界结点）

24二、有限单元法-伽辽金有限单元法有限单元剖分与基函数113114二、有限单元法-伽辽金有限单元法构造单元函数NL基函数取为“双线性插值”

基函数NL（x，y）形状如同一顶高等于1、有4条直线斜边和4条下凹型曲边的尖顶斗笠25二、有限单元法-伽辽金有限单元法构造单元函数NL基函数N114115二、有限单元法-伽辽金有限单元法子区DL以结点L为其共同结点的所有矩形单元（<=4）基函数NL仅在子区上不为0，在非DL上均为0,属于子区DL单元e的单元基函数NL，用NeL表示（>0）26二、有限单元法-伽辽金有限单元法子区DL115116二、有限单元法-伽辽金有限单元法典型矩形单元该单元中心位于坐标原点处，且4条边分别平行x，y轴,结点从左下角开始按逆时针方向编号i,j,m,n。沿x方向长2a，y方向宽2b。27二、有限单元法-伽辽金有限单元法典型矩形单元116117二、有限单元法-伽辽金有限单元法按上述要求所构造的NeL形式为对典型单元，令

则28二、有限单元法-伽辽金有限单元法按上述要求所构造的NeL117118二、有限单元法-伽辽金有限单元法其中即（6-34）29二、有限单元法-伽辽金有限单元法其中即（6-34）118119二、有限单元法-伽辽金有限单元法Nei(x,y)在结点i处为1，其它3处为0根据上式，有平行x或y方向上，Nei(x,y)线性变化故单元基函数为双线性插值函数若结点坐标用表示结点浓度用表示30二、有限单元法-伽辽金有限单元法Nei(x,y)在结点i119120二、有限单元法-伽辽金有限单元法结点浓度CL(t)和单元基函数NeL(x,y)来确定单元内的近似解，即对于区域D，结点L的基函数在子区DL内各单元分块确定。称为矩阵单元e上的基函数（6-36）31二、有限单元法-伽辽金有限单元法结点浓度CL(t)和单元120121二、有限单元法-伽辽金有限单元法伽辽金有限单元法

左端积分范围为D，由于在上NL=0,改写在子区DL的积分，对矩形单元逐个积分、求和。设子区DL上有neL个单元，有（6-39）32二、有限单元法-伽辽金有限单元法伽辽金有限单元法（6-3121122二、有限单元法-伽辽金有限单元法由（6-36）得由（6-34）得再由（6-36）得（6-40）（6-42）33二、有限单元法-伽辽金有限单元法由（6-36）得由（6-122123二、有限单元法-伽辽金有限单元法将（6-40）（6-42）代入（6-39）得34二、有限单元法-伽辽金有限单元法将（6-40）（6-4123124二、有限单元法-伽辽金有限单元法（6-44）35二、有限单元法-伽辽金有限单元法（6-44）124125二、有限单元法-伽辽金有限单元法用矩阵形式表示一阶微分方程组将一阶导数用差分形式表示按矩阵符号有若﹛C﹜在新时间水平t+Δt上取值，则36二、有限单元法-伽辽金有限单元法用矩阵形式表示一阶微分方125126二、有限单元法-伽辽金有限单元法（1）给定初始条件和边界条件，可求（6-46）各个系数矩阵，于是可解得第一个末时刻各结点的浓度；（2）以此浓度为已知浓度，计算新的系数矩阵；（3）再次求解，得到第二个时间步长末时刻结点浓度（6-46）浓度可在之间任意位置近似表示37二、有限单元法-伽辽金有限单元法（1）给定初始条件和边界126127二、有限单元法-伽辽金有限单元法按单元顺序计算系数矩阵每个矩阵单元涉及4个结点i，j，m，n，涉及4个方程单元e的作用，可用4行与4列的矩阵来表示，称为单元矩阵38二、有限单元法-伽辽金有限单元法按单元顺序计算系数矩阵127128二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体弥散矩阵的关系设研究区划分为8个单元，16个结点，四周为隔水边界编号如下：39二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体128129二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体弥散矩阵的关系计算每个单元弥散矩阵，按双下标已知，放置单元弥散矩阵并叠加，得到总体弥散矩阵。40二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体129130二、有限单元法-伽辽金有限单元法从（6-44）知将坐标原点移到矩形单元中心有41二、有限单元法-伽辽金有限单元法从（6-44）知将坐标原130131二、有限单元法-伽辽金有限单元法通过双下标一致的元素求和，得到总体弥散矩阵诸元素42二、有限单元法-伽辽金有限单元法通过双下标一致的元素求和131132二、有限单元法-伽辽金有限单元法积分求解可用高斯求积方法，详见P97-10543二、有限单元法-伽辽金有限单元法积分求解可用高斯求积方法132133二、有限单元法-伽辽金有限单元法44二、有限单元法-伽辽金有限单元法133134二、有限单元法-伽辽金有限单元法算例：二维水动力弥散的伽辽金法计算机程序

数学模型写成剖分成9个单元共20个结点，空间步长取Δx=Δy=10m，时间步长Δt=5d，渗透流速u=0.01m/d具体代码见P106-11145二、有限单元法-伽辽金有限单元法算例：二维水动力弥散的伽134135三、过程与数值弥散过量现象

如一维流动水动力弥散中单位浓度梯度函数注入的情况。

对比解析解与数值解C-x曲线，在浓度封面附近数值计算的浓度超过最大浓度值1和小于最小浓度0，失去物理意义。数值弥散

若纯对流方程采用有限差分逼近，其结果呈现似有弥散的存在。46三、过程与数值弥散过量现象数值弥散135136三、过程与数值弥散47三、过程与数值弥散136137三、过程与数值弥散过量现象-解析

由于时间步长与空间尺度匹配不佳，因而含水层在数值上不能“吸收”住污染物所引起。解决办法：

将时间增量选择为几何级数的子项；

做试验校正时间步长和差分网格；48三、过程与数值弥散过量现象-解析137138三、过程与数值弥散数值弥散-解析

由截断误差引起。

对流项一阶空间导数采用向后差分逼近，浓度Ci-1泰勒展开变形得

向后差分的截断误差主部分为相当于弥散系数的弥散项49三、过程与数值弥散数值弥散-解析相当于弥散系数的弥散项138139三、过程与数值弥散数值弥散-解析

由故纯对流方程有限差分近似表示后，结果如同对流-弥散方程的结果，记用数值弥散系数与物理弥散系数的比值来评价数值弥散对物理弥散的干涉程度此条件下的数值弥散系数50三、过程与数值弥散数值弥散-解析此条件下的数值弥散系数139140三、过程与数值弥散数值弥散-解析

定义为Peclet数，时，弥散占优；

时，对流占优。速度u与物理弥散系数DL比越大，越大。对于小数，数值解与精确解非常接近，对于大数，数值解过渡带变宽，浓度锋面变平缓，在锋面的前后还出现解的振动—过量现象。控制Δx的大小从而减少数值弥散51三、过程与数值弥散数值弥散-解析控制Δx的大小从而减少数140141三、过程与数值弥散对纯对流方程

52三、过程与数值弥散对纯对流方程141142三、过程与数值弥散对纯对流方程

结果也造成数值弥散此时，总的数值弥散是两者之和53三、过程与数值弥散对纯对流方程结果也造成数值弥散此时，总142143四、对流占主导地位时的数值方法上游加权法

对流项的有限差分化过程中乘上一个适当的权因子，则波动和过量得到改善或消除一维水动力弥散中对流项取中心隐式差分若相邻格点i-1或（和）i+1的浓度越大，则格点i的浓度也越大，反映到差分方程上，系数Gi-1和Gi+1不得与Gi同号，上式要求

不能满足时，对流占优势54四、对流占主导地位时的数值方法上游加权法对流项的有限差分143144四、对流占主导地位时的数值方法上游加权法

假定DL和u为常数，取等格距差分网格。一维对流-弥散方程在积分，有对流-弥散通量用有限差分近似表示

55四、对流占主导地位时的数值方法上游加权法假定DL和u为常144145四、对流占主导地位时的数值方法上游加权法

ω为权因子；当ω=1/2时，对流项相当于中心差分；当ω=1时，对流项相当于取向后差分，

当ω=0时，对流项相当于向前差分。速度不同，

ω取值不同56四、对流占主导地位时的数值方法上游加权法145146四、对流占主导地位时的数值方法上游加权法

上游浓度权因子大于1/2。物理意义：加强上游格点浓度值在对流项中的作用可通过浓度的线性（或二次）插值具体确定在Δt期间C在格边上随时间的变化，从而确定ω值。57四、对流占主导地位时的数值方法上游加权法上游浓度权因子大146147四、对流占主导地位时的数值方法

仍用Gi-1、Gi和Gi+1分别表示上式左端第一、二、三项括号内系数。则采用上游加权法有限差分格式，可消除过量和波动，但截断误差增大，数值弥散反而增大。中心差分格式，截断误差为二阶上游加权有限差分截断误差为一阶58四、对流占主导地位时的数值方法仍用Gi-1、Gi147148四、对流占主导地位时的数值方法

上游权因子的物理实质：对流项浓度C在某断面，例如断面的中心差分意味着对流项的浓度在Δt时段内均以i与i-1两格点的平均浓度，即通过断面；实际为：Δt时段的初始时刻是该平均浓度值，随后在对流的作用下，该断面的浓度显然会向上游格点的浓度Ci-1(u>0时)的变化。59四、对流占主导地位时的数值方法上游权因子的物理实148149四、对流占主导地位时的数值方法特征值方法

令源汇项为0，展开对流-弥散方程

60四、对流占主导地位时的数值方法特征值方法149150四、对流占主导地位时的数值方法特征值方法

假定流体密度值ρ不变（不可压缩），含水层不可压缩，则，（6-75）简化为基本思路：物质导数表示一个沿轨迹（特征曲线）

运动的研究对象浓度变化率。由弥散作用引起。若存在运会想，还将受到沿轨迹源与汇、化学、生物化学反应、吸附与解析等。对流项用沿特征线的示踪质点的运动描述

（6-76）61四、对流占主导地位时的数值方法特征值方法（6-76）150151四、对流占主导地位时的数值方法特征值方法

特征点均匀分布研究区，赋一个初始浓度。

含水层中未受污染部分特征点的初始浓度为零，特征点浓度眼轨迹的瞬时浓度变化，借助矩形网格差分格式计算。62四、对流占主导地位时的数值方法特征值方法151152四、对流占主导地位时的数值方法63四、对流占主导地位时的数值方法152153四、对流占主导地位时的数值方法模型精确度取决于网格距大小。可以在同一网格上计算速度场。如果采用较少的网格距，可以改善流场的描述。

坐标为的特征点所在的单元，可用下标i,j确定

精确度受到每一网格内所选择的最初的特征点数目影响。64四、对流占主导地位时的数值方法模型精确度153154四、对流占主导地位时的数值方法缺点：（1）实际程序冗长，需记录特征点的踪迹；

（2）抽水点或出流边界的特征点需要消除，入渗边界须有特征点生成，不透水边界特征点须反射回来；

（3）要通过消去特征点来避免滞留点附近特征点堆积；（4）特征点耗尽的单元，必须产生新的特征点。65四、对流占主导地位时的数值方法缺点：154155四、对流占主导地位时的数值方法动坐标系方法与网格变形方法能消除过量现象，还保持较陡的浓度锋面。

对式

设坐标变换

得，

取空间步长，则

66四、对流占主导地位时的数值方法动坐标系方法与网格变形方法155156四、对流占主导地位时的数值方法动坐标系方法与网格变形方法

在固定坐标系中取格距为，有

动坐标系中格点运动与固定坐标系格点运动重合。

设固定坐标系中格点xj在tn+1时刻与动坐标系中Xi重合，有

动坐标系下浓度结点的表达式67四、对流占主导地位时的数值方法动坐标系方法与网格变形方法156157四、对流占主导地位时的数值方法动坐标系方法与网格变形方法

式变成

对半无限长砂柱，当DL=0时68四、对流占主导地位时的数值方法动坐标系方法与网格变形方法157158四、对流占主导地位时的数值方法网格变形方法

设tn时刻有限元结点在空间中的分布是对应浓度分布为

69四、对流占主导地位时的数值方法网格变形方法158159四、对流占主导地位时的数值方法网格变形方法

主要步骤：结点位置随时间变化，基函数及有限元方程的系数都依赖于时间，每推进一个Δt都需要重新计算。二维问题还须进行畸形判断

70四、对流占主导地位时的数值方法网格变形方法159160四、对流占主导地位时的数值方法结合动坐标的网格变形方法

一维对流-弥散在某个动点在t时刻位置x可写成

其中，x0是动点的初始位置，有取时间的一阶偏导数对流-弥散方程写成

若控制动点速度dx/dt并让它接近流速u，则可转换成弥散为主的的方程式71四、对流占主导地位时的数值方法结合动坐标的网格变形方法若160161四、对流占主导地位时的数值方法若使用伽辽金法，基函数因依赖于结点位置而成为时间的函数。若用表示随时间变化的基函数，则对任意t，有

故72四、对流占主导地位时的数值方法若使用伽辽金法，基函数因依161162四、对流占主