1、地下水流数值模拟方法地下水流数值模拟方法第六章 水动力弥散方程的数值解法2一、有限差分法一、有限差分法3一、有限差分法一、有限差分法-导数的有限差分近似导数的有限差分近似xxi ixx1-x1-i ix1x1ii4一、有限差分法一、有限差分法-导数的有限差分近似导数的有限差分近似x2分别一阶导数的向前差分、向后差分及中心差分分别一阶导数的向前差分、向后差分及中心差分2x二阶中心差分二阶中心差分5一、有限差分法一、有限差分法-一维水动力弥散的差分解法一维水动力弥散的差分解法（6-76-7）1niC11122221()2()()2nnnnLLLiiiiDtu tDtDtu tCCCCxxxxx（1
2、）显格式式中仅有一个未知数，解得 1111122()2nnnnnnniiiiiiiLCCCCCCCDutxx式（6-7）中的对流项取中心差分 6可以证明，稳定性准则要求可以证明，稳定性准则要求221()LDtx2()2LxtD 22()Lu tDtxx()22LLLLDa uDxau 即 即 。（1）（2）。2Lxa xt22()222LLLLLDa uxDatDuu由条件（由条件（2），格距要求很小；由（，格距要求很小；由（2）可知，鉴于）可知，鉴于较小，导致较小，导致不能太大，将（不能太大，将（2）代入（）代入（1）式中，得到）式中，得到 7111122()nnnnnnniiiiiiiLC
3、CCCCCCDutxx11122212()()()nnnnLLLiiiiDtu tu tDtDtCCCCxxxxx2(1)22(1)LLLxxxDtDDuxuuxx 若对流项改为向后差分若对流项改为向后差分 解得解得 稳定性要求稳定性要求不难看出，稳定性限制比对流项取中心差分有所改善。不难看出，稳定性限制比对流项取中心差分有所改善。8（2 2）隐式格式）隐式格式111111111122()2nnnnnnniiiiiiiLCCCCCCCDutxx 1111122221()2()()2nnnnLLLiiiiDtu tDtDtu tCCCCxxxxx整理后得到整理后得到 隐格式是无条件稳定的。9（2
4、 2）Grank-NicolsonGrank-Nicolson格式格式 整理后得到整理后得到 取隐式和显示的平均，即取隐式和显示的平均，即Grank-Nicolson格式xCCxCCCDtCninniLn2u-221C,1in12,1i,in1,i1ni,，（），xCCxCCCDninniiL2u-21,1i1n121,1i1,1-n11,121,21,122-t-122t-niLniLniLCxtuxDCxtDCxtuxDniLniLniLCxtuxDCxtDCxtuxD,1-2,2,122-t-122t10一、有限差分法一、有限差分法-二维水动力弥散的差分解法维水动力弥散的差分解法（4-5
5、64-56）（1）显格式式中仅有一个未知数式（4-56）中的对流项取中心差分 2,1i,in,1,i1nj,i,2CxCCCDtCnjnjjiLnj，xCCyCCCDnjjinjnjjiT2u-2,1,in,12,1,i,in,1，11一、有限差分法一、有限差分法-二维水动力弥散的差分解法维水动力弥散的差分解法化简后，有njiTLnjiLnjiCyDxDCxDC,22,121,t2t21x2tu-tnjiTnjiTnjiLCyDCyDCxD,1-,2,1,2,1-2ttx2tu-t涉及以（i，j）为中心的5个网格点在tn时刻的已知浓度12一、有限差分法一、有限差分法-二维水动力弥散的差分解法维
6、水动力弥散的差分解法（2）隐格式式（4-56）中右端的对流项取中心差分，右端个C的时阶均取n+1水平 21,1i,1,in,11,i1nj,i,2CxCCCDtCnjnjjiLnj，xCCyCCCDnjjinjnjjiT2u-21,i1n,121,1,i1,in,1，13一、有限差分法一、有限差分法-二维水动力弥散的差分解法维水动力弥散的差分解法（2）隐格式整理后收敛且无条件稳定涉及以（i，j）为中心的5个网格点在tn+1时刻的未知浓度14一、有限差分法一、有限差分法-二维水动力弥散的差分解法维水动力弥散的差分解法（3）Grank-Nicolson格式将隐式格式的两式相加除以215一、有限差分
7、法一、有限差分法-二维水动力弥散的差分解法维水动力弥散的差分解法（3）Grank-Nicolson格式整理得16一、有限差分法一、有限差分法-二维水动力弥散的差分解法维水动力弥散的差分解法（4）交替方向隐式法（ADI法）分两次对三对角矩阵求逆，将一个t分成两个t/2计算 第一个半时间步，对x方向的偏导数采用隐式差分，对y方向的偏导数采用显示差分。2,1,i,in,1221,1i21,i21n,1-22yCCCDxCCCDnjnjjiTnjnjjiL，2tu-j,i21nj21j,1i21-njniniCCxCC，17一、有限差分法一、有限差分法-二维水动力弥散的差分解法维水动力弥散的差分解法（
8、4）交替方向隐式法（ADI法）整理得 第二个半时间步，对y方向的偏导数采用隐式差分，对x方向的偏导数采用显示差分。21,1,i1,i1n,1221,1i21,i21n,122yCCCDxCCCDnjnjjiTnjnjjiL，18一、有限差分法一、有限差分法-二维水动力弥散的差分解法维水动力弥散的差分解法（4）交替方向隐式法（ADI法）整理得收敛且无条件稳定2tu-j,i21nj21j,1i21-njniniCCxCC，19一、有限差分法一、有限差分法-求解差分方程的计算机程序举例求解差分方程的计算机程序举例算例算例 见教材见教材P81-85P81-8520二、有限单元法二、有限单元法-伽辽金有
9、限单元法（6-6-2626）21二、有限单元法二、有限单元法-伽辽金有限单元法C（6-6-2929）22二、有限单元法二、有限单元法-伽辽金有限单元法23二、有限单元法二、有限单元法-伽辽金有限单元法（6-6-3232）代入（代入（6-296-29）24二、有限单元法二、有限单元法-伽辽金有限单元法25二、有限单元法二、有限单元法-伽辽金有限单元法基函数NL（x，y）形状如同一顶高等于1、有4条直线斜边和4条下凹型曲边的尖顶斗笠26二、有限单元法二、有限单元法-伽辽金有限单元法27二、有限单元法二、有限单元法-伽辽金有限单元法28二、有限单元法二、有限单元法-伽辽金有限单元法29二、有限单元法
10、二、有限单元法-伽辽金有限单元法（6-6-3434）30二、有限单元法二、有限单元法-伽辽金有限单元法nnmmjjiiyCyCyCyC,x,x,x,x、nnmmjjiiyyyy,x,x,x,x、31二、有限单元法二、有限单元法-伽辽金有限单元法txC,y,eyNyNyNyNeeee,x,x,x,xnmji、（6-6-3636）32二、有限单元法二、有限单元法-伽辽金有限单元法LDxy,4neL33二、有限单元法二、有限单元法-伽辽金有限单元法34二、有限单元法二、有限单元法-伽辽金有限单元法35二、有限单元法二、有限单元法-伽辽金有限单元法36二、有限单元法二、有限单元法-伽辽金有限单元法37
11、二、有限单元法二、有限单元法-伽辽金有限单元法1t 1tC lC B 2tt C Cttt与38二、有限单元法二、有限单元法-伽辽金有限单元法39二、有限单元法二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体弥散矩阵的关系算例：单元弥散矩阵与整体弥散矩阵的关系设研究区划分为设研究区划分为8 8个单元，个单元，1616个结点，四周为隔水边个结点，四周为隔水边界编号如下：界编号如下：40二、有限单元法二、有限单元法-伽辽金有限单元法算例：单元弥散矩阵与整体弥散矩阵的关系算例：单元弥散矩阵与整体弥散矩阵的关系计算每个单元弥散矩阵，按双下标已知，放置单元弥散矩计算每个单元弥散矩阵，按双下标已知
12、，放置单元弥散矩阵并叠加，得到总体弥散矩阵。阵并叠加，得到总体弥散矩阵。41二、有限单元法二、有限单元法-伽辽金有限单元法42二、有限单元法二、有限单元法-伽辽金有限单元法43二、有限单元法二、有限单元法-伽辽金有限单元法积分求解可用高积分求解可用高斯求积方法，详斯求积方法，详见见P97-105P97-10544二、有限单元法二、有限单元法-伽辽金有限单元法45二、有限单元法二、有限单元法-伽辽金有限单元法算例：二维水动力弥散的伽辽金法计算机程序算例：二维水动力弥散的伽辽金法计算机程序 数学模型写成数学模型写成46三、过程与数值弥散三、过程与数值弥散47三、过程与数值弥散三、过程与数值弥散48
13、三、过程与数值弥散三、过程与数值弥散49三、过程与数值弥散三、过程与数值弥散相当于弥散系数的弥散项相当于弥散系数的弥散项50三、过程与数值弥散三、过程与数值弥散此条件下的数值弥散系数此条件下的数值弥散系数LeeLLDxuPPDxuDD即22n51三、过程与数值弥散三、过程与数值弥散eP2eP2ePePePeP控制控制x x的大小从而减少数值弥散的大小从而减少数值弥散52三、过程与数值弥散三、过程与数值弥散53三、过程与数值弥散三、过程与数值弥散结果也造成数值弥散结果也造成数值弥散54四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法对流项的有限差分化过程中乘上一个适当的对流项的有限差
14、分化过程中乘上一个适当的权因子，则波动和过量得到改善或消除权因子，则波动和过量得到改善或消除55四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法2/1i2/1ix,x56四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法57四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法58四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法59四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法2/1ix2C1iiC2/1ix60四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法61四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法0t0
15、eH或0div tyt,x62四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法63四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法64四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法 ttppyx、65四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法66四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法67四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法动坐标系下浓度结点的表达式动坐标系下浓度结点的表达式68四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法69四、对流占主导地位时的数值方法四、对流占主导地位时的数
16、值方法nnXCnX70四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法71四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法若控制动点速度若控制动点速度dx/dtdx/dt 并让它接近流速并让它接近流速u u，则可转换成弥散为主的的方程式，则可转换成弥散为主的的方程式72四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法t,xNi txNtCtxCNii,i173四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法74四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法t4exp4,20LLDutxtDCtxC75四、对流占主导地位时的数值方法四
17、、对流占主导地位时的数值方法N76四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法77四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法78四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法对流运动可确定在水流流动方对流运动可确定在水流流动方向及其正交方向的弥散运动向及其正交方向的弥散运动79四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法80四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法8182四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法83四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法84四、对流占主导
18、地位时的数值方法四、对流占主导地位时的数值方法0t,ijijiixCDxxCuC0t,ijjiiiCDCuC1ji1jij1i1nj,DDundnidnidCCCtCi0D1D1u1jijij1ij,ndidniunCCCtCi85四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法tCtDtCtujiijdii2/12/1ujiijjiijiiiCDCDCtCu02/1-2ijdttCtDji86四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法87四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法ijDd88四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法0u89四、对流占主导地位时的数值方法四、对流占主导地位时的数值方法
