第7章土壤水分运动的数值解

1、7.1 7.1 概述概述7.2 7.2 差分的概念及分类差分的概念及分类7.3 7.3 显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分7.4 7.4 土壤水分运移方程的线性化方法土壤水分运移方程的线性化方法7.5 7.5 线性化方法与土壤水分参数的取值线性化方法与土壤水分参数的取值7.6 7.6 边界条件的处理及追赶法求解三对角方程边界条件的处理及追赶法求解三对角方程7.7 7.7 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图Huyakorn, P. S. and G. F. Pinder. Computational methods in subsurface

2、flow. Academic Press, New York, 1983.薛禹群，谢春红，水文地质学的数值法，煤炭工业出版社，薛禹群，谢春红，水文地质学的数值法，煤炭工业出版社，1980.土壤水分运动方程的求解方法：土壤水分运动方程的求解方法： 解析解（解析解（Analytical solution） 数值解（数值解（Numerical solution）(1) 有限差分法（有限差分法（Finite Difference method）(2) 有限单元法（有限单元法（Finite Element method）7.1 7.1 概述概述 差分的由来（差分的由来（Taylor展开）展开） if =

3、 (x) exists, thenkkkxkxxxxxx!)()(kkkkxkxxxxxx!) 1()()(1)(2)(0)()(xxxxxx)(0)()(xxxxxxwhen x 0, thenxxxxxxxxxxx)()(lim)()(lim00 向前差分向前差分xxxxx)()( 向后差分向后差分xxxxx)()( 中心差分（中心差分（Eqn.(1) - Eqn. (2)）xxxxxxxx)(02)()(2xxxxxx2)()(i.e.截断误差截断误差0 ( x)截断误差截断误差0 ( x)截断误差截断误差0 ( x2)0 x ABCEFxi-1xixi+1如图所示：如图所示： BC表示

4、向前差分；表示向前差分； AB表示向后差分；表示向后差分； AC表示中心差分；表示中心差分；In addition, Eqn.(1) + Eqn.(2) )(0)(2)()(22222xxxxxxxxx)(0)()(2)(2222xxxxxxxxi.e.:其截断误差亦为与其截断误差亦为与 x2同阶的无穷小量。同阶的无穷小量。i=1, z(1)=0i=2i=3.i=n, z(n)=Lz考虚垂直一维问题（如右图所示），考虚垂直一维问题（如右图所示），泛定方程可写为（以泛定方程可写为（以 方程为例）：方程为例）：剖面深度为剖面深度为L，共剖分出，共剖分出n个节点（个节点（i = 1, 2, 3, ,

5、 n），），z (1) = 0, z (n) = L，已，已知知t = tk时刻各节点含水率时刻各节点含水率 的分布，求的分布，求t = t k+1时刻各节点含水率时刻各节点含水率 的分布，的分布， kkttt1zKzDzt 7.3.17.3.1显式差分格式显式差分格式kiizKzDztttkikii121212121iikikikizzzDzDzDz iikikikikizzDzD112121 112121iikikikikizzDzD7.3.17.3.1显式差分格式显式差分格式 2112111211121 zzzzDzzDzDziiikikikiiikikikiki21212121iiki

6、kikizzKKzK21212121112111211iiiiiikikiiiikikiikikizzKKzzDzzDt)(21112121iiiizzzz) 1, 3, 2( 2212111211121111niKKzzDzzDzztiiiikikiiiikikiiiikiki7.3.17.3.1显式差分格式显式差分格式 若取等步长若取等步长 z，即，即 z =zi+1 - zi = zi - zi-1，则有：，则有：) 1, 3, 2(212112112121niKKzDDztiikikiikikiikiki7.3.17.3.1显式差分格式显式差分格式 显式差分格式是有条件收敛的显式差分格

7、式是有条件收敛的，一般应满足，一般应满足max2Dzrt其中：其中：r多取多取1/2。由于土壤接近饱和时，。由于土壤接近饱和时，Dmax很很大，故大，故 t一般要求取值很小，耗费机时。一般要求取值很小，耗费机时。7.3.17.3.1显式差分格式显式差分格式 1kiiKzDzt) 1, 3, 2( 22/12/1111121111121111niKKzzDzzDzztiiiikikiiiikikiiiikiki7.3.27.3.2隐式差分格式隐式差分格式 )(1111iiiizzzztr)(1112iiiizzzztrLetikiikiikiiHGFE11111then ( i = 2, 3,

8、, n - 1)三对角方程三对角方程7.3.27.3.2隐式差分格式隐式差分格式 where;211212iiiiDrGDrE11211212iiiiiDrDrGEF2121112iiiikiiKKzztH隐式差分格式是无条件收敛的。隐式差分格式是无条件收敛的。7.3.27.3.2隐式差分格式隐式差分格式 kikikiiKzDzKzDzKzDzt12/1217.3.37.3.3中心差分格式中心差分格式 ) 1, 3, 2()(2)(22212111111211111121111niKKzzDzzDzztiiiikikikikiiiikikikikiiiikiki7.3.37.3.3中心差分格式

9、中心差分格式 )(1111iiiizzzztr)(1112iiiizzzztrLetikiikiikiiHGFE11111then ( i = 2, 3, , n - 1)where;211212iiiiDrGDrE7.3.37.3.3中心差分格式中心差分格式 11211212iiiiiDrDrGEF中心差分格式也是无条件收敛的。中心差分格式也是无条件收敛的。212111112) 1(iiiikiikiiikiiiKKzztGGEEH7.3.37.3.3中心差分格式中心差分格式 7.4 7.4 土壤水分运移方程的线性化方法土壤水分运移方程的线性化方法 以基质势水头以基质势水头h为因变量的一维垂

10、直向土壤为因变量的一维垂直向土壤水分运动基本方程为：水分运动基本方程为： K hhhC hK htzzz7.4 7.4 土壤水分运移方程的线性化方法土壤水分运移方程的线性化方法7.4.17.4.1显式差分格式显式差分格式相应h方程的差分方程为：11/211/21211kkkkkkkkkiiiiiiiiiChhKhhKhhtz1/21/2kkiiKKz7.4 7.4 土壤水分运移方程的线性化方法土壤水分运移方程的线性化方法7.4.27.4.2隐式差分格式隐式差分格式相应h方程的差分方程为1111111111/211/2111/21/22kkkkkkkkkkiiiiiiiikiiihhKhhKhh

11、KKCtzz7.4 7.4 土壤水分运移方程的线性化方法土壤水分运移方程的线性化方法7.4.37.4.3中心差分格式中心差分格式相应h方程的差分方程为1111111221221122111111122212222kkkkkkkkiiiiiiiikkkkkkiikiiiiihhhhhhCKtzKKhhhhKzz 土壤水分运移方程中的各参数均依赖于变量土壤水分运移方程中的各参数均依赖于变量 或或 m，从而使得方程呈现出较强的非线性性，求解前必须将其线从而使得方程呈现出较强的非线性性，求解前必须将其线性化，得到性化，得到n元一次线性代数方程组，以便于求解。元一次线性代数方程组，以便于求解。（1）显式

12、线性化）显式线性化 计算过程中，方程中的各参数如：计算过程中，方程中的各参数如：C ( m)，K( m) or D( )，K( )等均以时段初的值（即前一时段的值）代入。等均以时段初的值（即前一时段的值）代入。 适用条件：此法只能用于适用条件：此法只能用于 m或或 变化缓慢的情况；当变化缓慢的情况；当变化剧烈时，此法会导致较大的偏差。变化剧烈时，此法会导致较大的偏差。7.5.17.5.1线性化方法线性化方法例如：已知例如：已知k - 1与与k时刻各节点的时刻各节点的 值，需求值，需求k + 1时刻各时刻各节点的节点的 值，应用显式外推线性化参数的方法如下（以值，应用显式外推线性化参数的方法如下

13、（以导水率导水率K的求法为例，其他各参数的求法与此类似）：的求法为例，其他各参数的求法与此类似）：1111)(kkkkkikikikitttt(*) 求求 ： 方法（方法（a）：）：21kiK)(21121kikiki2121kikiKK方法方法(b)：算术平均算术平均:几何平均几何平均:调和平均调和平均:121)(21kikikiKKK)()(121kikikiKKK)()()()(21121kikikikikiKKKKK 求求 ： 方法（方法（a）：）：kiK21)(21121kikikikikiKK2121方法方法(b)：算术平均算术平均几何平均几何平均调和平均调和平均kikikiKKK

14、121)(21)()(121kikikiKKK)()()()(21121kikikikikiKKKKK 计算时，先用显式线性化的方法，采用隐式差分格式，计算时，先用显式线性化的方法，采用隐式差分格式，求出求出k + 1/2时刻各节点的时刻各节点的 m或或 （即（即 or ）；根据）；根据这些值可求得相应的这些值可求得相应的 or 等参数，以此作为计算时始末（等参数，以此作为计算时始末（tk to tk+1）的平均值，代）的平均值，代入方程求解入方程求解 ，从而完成方程的线性化工作。，从而完成方程的线性化工作。21kmi21ki2121kmikiKK21kiK1ki7.5.17.5.1线性化方法

15、线性化方法7.5.1.27.5.1.2预报校正法预报校正法（1）显式外推法（先求显式差分方程，得）显式外推法（先求显式差分方程，得 作为预报作为预报值）值） 计算参数时，据前一时段始末的函数值计算参数时，据前一时段始末的函数值 m或或 ，用，用线性外推近似求出本时段末的函数值，然后线性外推近似求出本时段末的函数值，然后 方法方法(a)：求出计算时段始末的平均值，再通过各参数与：求出计算时段始末的平均值，再通过各参数与 m或或 的关系，求出各参数；或的关系，求出各参数；或 方法方法(b)：先求出计算时段始末的参数值，然后对参数取：先求出计算时段始末的参数值，然后对参数取平均。平均。 平均值的计算

16、可采用算术平均、几何平均、调和平平均值的计算可采用算术平均、几何平均、调和平均等方法。均等方法。7.5.17.5.1线性化方法线性化方法7.5.1.27.5.1.2预报校正法预报校正法1ki（2）线性外推法（已知线性外推法（已知 时刻，外推时刻，外推 时刻）时刻） 假定含水率在相邻时段内是线性变化的，根据前假定含水率在相邻时段内是线性变化的，根据前一个时段末的含水率，用直线外推近似求出计算时段一个时段末的含水率，用直线外推近似求出计算时段末的含水率：末的含水率：7.5.17.5.1线性化方法线性化方法7.5.1.27.5.1.2预报校正法预报校正法,1k k 1k 112kkkiii11231

17、22kkkiii 对于第一时段，无法应用线性外推法，可用其他对于第一时段，无法应用线性外推法，可用其他线性化方法。线性化方法。 7.5.17.5.1线性化方法线性化方法7.5.1.27.5.1.2预报校正法预报校正法（3 3）显式预测法（两次求解隐式：第一次预报值；）显式预测法（两次求解隐式：第一次预报值；第二次得校正值）第二次得校正值） 首先利用显式线性化方法，由时段初的含水率首先利用显式线性化方法，由时段初的含水率及相应的参数值，按隐式差分格式求解方程，求及相应的参数值，按隐式差分格式求解方程，求得本时段末（或时段中间）的含水率及相应的参得本时段末（或时段中间）的含水率及相应的参数值（预报

18、值），再求解隐式差分方程，得到时数值（预报值），再求解隐式差分方程，得到时段末的含水率的值（校正值）。段末的含水率的值（校正值）。计算时，先假定本时段末（计算时，先假定本时段末（tk+1）的）的 m或或 （一般做（一般做法是先取时段初的函数值法是先取时段初的函数值 or ） ， or 等参数的迭代初值等参数的迭代初值 求解线性代数方程组得求解线性代数方程组得 or 直至前后两次迭代值的误差满足精度要求为止。直至前后两次迭代值的误差满足精度要求为止。kmiki)1(1kiK)1(1kiD)1(1kiC)1(1kmi)1 (1ki)2(1)2(1)2(1or ,kikikiCDK)2(1)2(1o

19、r kikmi)(1)(1or pkipkmi7.5.1.3迭代法迭代法（a）若取绝对误差，判断条件为：）若取绝对误差，判断条件为：), 2, 1(max1) 1( 1)( 1nipkipkii), 2, 1(max1)1(1)(1nipkmipkmiior7.5.17.5.1线性化方法线性化方法7.5.1.3迭代法迭代法（b）若取相对误差，判断条件为：）若取相对误差，判断条件为：), 2, 1(max2)(1)1(1)(1nipkipkipkii), 2, 1(max2)(1)1(1)(1nipkmipkmipkmiior7.5.17.5.1线性化方法线性化方法7.5.1.3迭代法 1与与

20、2 的取值可根据实际问题的精度要求而定。的取值可根据实际问题的精度要求而定。 此法计算稍复杂，但精度较高，数值计算时此法计算稍复杂，但精度较高，数值计算时经常采用。为避免迭代次数太多，一般可设置一经常采用。为避免迭代次数太多，一般可设置一个判断语句：当迭代次数个判断语句：当迭代次数m5时，将时间步长时，将时间步长 减半，再重新计算。减半，再重新计算。kkkttt117.5.17.5.1线性化方法线性化方法7.5.1.3迭代法7.5 线性化方法与土壤水分参数的取值线性化方法与土壤水分参数的取值7.5.2.1直接法直接法先求出两结点的参数，然后先求出两结点的参数，然后(1)取两结点参数的算术平均值

21、：取两结点参数的算术平均值：7.5.2参数取值参数取值11212iiiDDD(2)取两结点参数的调和平均值：取两结点参数的调和平均值：11122/iiiiiDDDDD(3) 取两结点参数的几何平均值：取两结点参数的几何平均值：112iiiDDD7.5 线性化方法与土壤水分参数的取值线性化方法与土壤水分参数的取值7.5.2 2间接法间接法先求出两结点处含水率的算术、调和、几何平均值，先求出两结点处含水率的算术、调和、几何平均值，然后再求出相应的参数值。然后再求出相应的参数值。(1)算术平均值：算术平均值：11212iii1122iiDD(2)调和平均值：调和平均值：11122/i iiii112

22、2iiDD(3)几何平均值：几何平均值：112i ii1122iiDD7.5 线性化方法与土壤水分参数的取值线性化方法与土壤水分参数的取值7.5.2 3其他方法其他方法（1）求出的参数和两结点处的参数和，三点参）求出的参数和两结点处的参数和，三点参数取平均，故可称为数取平均，故可称为“三点式三点式”：11122124iiiiDDDD（2） 11121iiiiiDDd （3）由邻近结点的参数代替，如）由邻近结点的参数代替，如11122,iiiiDD DD前已叙及，采用显式差分格式可直接求解，采前已叙及，采用显式差分格式可直接求解，采用隐式或中心差分格式时总可以得到如下形式的三对用隐式或中心差分格

23、式时总可以得到如下形式的三对角方程组（以角方程组（以 方程为例）：方程为例）：共共n - 2个方程，再加上上、下边界条件，即可组成个方程，再加上上、下边界条件，即可组成n n 阶的阶的n元一次方程组，线性化后求解该方程组可获得元一次方程组，线性化后求解该方程组可获得原定解问题的解。原定解问题的解。) 1, 3 , 2( 11111niHGFEikiikiikii7.6.17.6.1边界条件的处理边界条件的处理 上、下边界不论为第一、第二或第三类边界，上、下边界不论为第一、第二或第三类边界，总可以化为：总可以化为：1121111HGFkknknnknnHFE111(2)(1)7.6.17.6.1

24、边界条件的处理边界条件的处理 验验 证：证： 当上、下边界为第一类边界时当上、下边界为第一类边界时)()()(),()(), 0(121111121kknkkztftftftLtft)(, 0, 1)2.(Eqn)(, 0, 1) 1.(Eqn12n11111knnktfHEFtfHGF7.6.17.6.1边界条件的处理边界条件的处理1121111HGFkknknnknnHFE111 当上、下边界为第二、三类边界时当上、下边界为第二、三类边界时)4()()3(), 0(210tqKzDttqKzDzLzz)(),(122/111112111112312111223knnnknknnkkkktq

25、KzzDtqKzzD1121111HGFkknknnknnHFE111211212123111111112231)(,),(,nknnnnnnnkkKtqHEFzzDEKtqHFGzzDFnknnknnnknnknnknnkkkkkkkkHFEHGFEHGFEHGFEHGF111111111121314313312321321221121121111经线性化后，原定解问题可化为经线性化后，原定解问题可化为n n阶线性代数方程组，即：阶线性代数方程组，即：以矩阵的形式表示即为：以矩阵的形式表示即为：nnknknkkknnnnnHHHHHFEGFEGFEGFEGF13211111312111113

26、33222117.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程1原理原理1121111HGFkk(1)(1a) ,( 111111111111HBFAAGABkk7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程1原理原理2132122112HGFEkkk1122132121212ABEHGAEGFkk13222212kkAGAB(2)(2a)1122212122,ABEHBAEGFA.1111111iiiikiikiiiiiABEHGAEGF1111,iiiiiiiiiiABEHBAEGFA.111kiiiiikiAGAB(ia)(i)ikiikiikiiHGFE1

27、11117.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程111111121nknnknnknnHGFE221111112121nnnnknnknnnnnABEHGAEGF1111111knnnnnknAGAB2211121211,nnnnnnnnnnABEHBAEGFA(n-1)a)(n-1)7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程nknnknnHFE11111111nnnnknnnnnABEHAEGF11111nnnnnnnnknAEGFABEH(na)(n)7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程 解出解出 后，利用式（后，利用式（

28、(n-1)a）,，(ia),，(2a), (1a)往前递推，依次即可解出往前递推，依次即可解出 ，于，于是得到是得到k+1时刻各节点的时刻各节点的 值（值（i = 1, 2, ，n）。）。 若若 均为已知（第一类边界条件），则可直均为已知（第一类边界条件），则可直接从式（接从式（2）开始求解，依次递推至式（）开始求解，依次递推至式（n - 1）即可。）即可。1kn11121211 , , , , kkknkn111 , knk7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程(1) 输入输入D、K含水量含水量SITA关系（建议采用关系（建议采用函数函数子程序调用）；子程序调用）；节

29、点剖分总数节点剖分总数n及对应的空间坐标及对应的空间坐标z(1:n)（若为等步长，则（若为等步长，则输入空间步长输入空间步长DZ）；）；初始含水量剖面分布初始含水量剖面分布SITA0（1:n）；）；上、下边界条件上、下边界条件UBC、LBC（=1，2，or3）及相应的边界）及相应的边界值值QSUBC、QSLBCT（建议采用子程序调用）；（建议采用子程序调用）；初始时间步长初始时间步长DT0、最大时间步长、最大时间步长DTMAX，计算结束时，计算结束时间间TEND，迭代控制标准，迭代控制标准DET。(2) 赋初值赋初值T(0)=0 初始时间初始时间T(1)=DT0 第一次计算时间第一次计算时间*（以下赋含水量计算中间值）（以下赋含水量计算中间值） do i=1, nSITAM(i)=SITA0(i) enddo时段计数时段计数KT=1(3) 判断计算时间是否已超过需模拟时间判断计算时间是否已超过需模拟时间if (T(KT)=DTMAX) then DT=DTMAX T(KT)=T(KT-1)+DT endif(5) 开始计算系数矩阵，求解线性代数方程组开始计算系数矩阵，求解线性代数方程组 赋迭代计次数赋迭代计次数 MP1=0(6) if (MP1=5)then (else goto (12)(7) 计算本时