第7章土壤水分运动的数值解课件.ppt

《第7章土壤水分运动的数值解课件.ppt》由会员分享，可在线阅读，更多相关《第7章土壤水分运动的数值解课件.ppt（58页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、第第7 7章章 土壤水分运动的数值解土壤水分运动的数值解 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 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图第第7 7章章 土壤水分运动的数值解土壤水分运动的数值解 主要参考书目主要参考书目

2、主要参考书目主要参考书目Huyakorn,P.S.and G.F.Pinder.Computational methods in subsurface flow.Academic Press,New York,1983.薛禹群，谢春红，水文地质学的数值法，煤炭工业出版社，薛禹群，谢春红，水文地质学的数值法，煤炭工业出版社，1980.土壤水分运动方程的求解方法：土壤水分运动方程的求解方法：解析解（解析解（Analytical solution）数值解（数值解（Numerical solution）(1)有限差分法（有限差分法（Finite Difference method）(2)有限单元法（有

3、限单元法（Finite Element method）7.1 7.1 概述概述7.2 7.2 差分的概念及分类差分的概念及分类 差分的由来（差分的由来（Taylor展开）展开）if =(x)exists,then(1)(2)第第7 7章章 土壤水分运动的数值解土壤水分运动的数值解 when x 0,then 向前差分向前差分 向后差分向后差分 中心差分（中心差分（Eqn.(1)-Eqn.(2)）i.e.截断误差截断误差0(x)截断误差截断误差0(x)截断误差截断误差0(x2)0 x ABCEFxi-1xixi+1如图所示：如图所示：BC表示向前差分；表示向前差分；AB表示向后差分；表示向后差分

4、；AC表示中心差分；表示中心差分；In addition,Eqn.(1)+Eqn.(2)i.e.:其截断误差亦为与其截断误差亦为与 x2同阶的无穷小量。同阶的无穷小量。i=1,z(1)=0i=2i=3.i=n,z(n)=Lz考虚垂直一维问题（如右图所示），考虚垂直一维问题（如右图所示），泛定方程可写为（以泛定方程可写为（以 方程为例）：方程为例）：剖面深度为剖面深度为L，共剖分出共剖分出n个节点（个节点（i=1,2,3,n），），z(1)=0,z(n)=L，已知已知t=tk时刻各节点含水率时刻各节点含水率 的分布，的分布，求求t=t k+1时刻各节点含水率时刻各节点含水率 的分布，的分布，7.

5、37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 7.3.17.3.1显式差分格式显式差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 7.3.17.3.1显式差分格式显式差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 7.3.17.3.1显式差分格式显式差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 若取等步长若取等步长 z，即即 z=zi+1-zi=zi-zi-1，则有：则有：7.3.17.3.1显式差分格式显式差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分

6、 显式差分格式是有条件收敛的显式差分格式是有条件收敛的，一般应满足，一般应满足其中：其中：r多取多取1/2。由于土壤接近饱和时，。由于土壤接近饱和时，Dmax很很大，故大，故 t一般要求取值很小，耗费机时。一般要求取值很小，耗费机时。7.3.17.3.1显式差分格式显式差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 7.3.27.3.2隐式差分格式隐式差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 Letthen (i=2,3,n-1)三对角方程三对角方程7.3.27.3.2隐式差分格式隐式差分格式7.37.3显式差分、隐式差分和中心

7、差分显式差分、隐式差分和中心差分 where隐式差分格式是无条件收敛的。隐式差分格式是无条件收敛的。7.3.27.3.2隐式差分格式隐式差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 7.3.37.3.3中心差分格式中心差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 7.3.37.3.3中心差分格式中心差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 Letthen (i=2,3,n-1)where7.3.37.3.3中心差分格式中心差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差

8、分 中心差分格式也是无条件收敛的。中心差分格式也是无条件收敛的。7.3.37.3.3中心差分格式中心差分格式7.37.3显式差分、隐式差分和中心差分显式差分、隐式差分和中心差分 7.4 7.4 土壤水分运移方程的线性化方法土壤水分运移方程的线性化方法 以基质势水头以基质势水头h为因变量的一维垂直向土壤为因变量的一维垂直向土壤水分运动基本方程为：水分运动基本方程为：7.4 7.4 土壤水分运移方程的线性化方法土壤水分运移方程的线性化方法7.4.17.4.1显式差分格式显式差分格式相应h方程的差分方程为：7.4 7.4 土壤水分运移方程的线性化方法土壤水分运移方程的线性化方法7.4.27.4.2隐

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

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

11、此类似）：(*)求求 ：方法（方法（a）：）：方法方法(b)：算术平均算术平均:几何平均几何平均:调和平均调和平均:求求 ：方法（方法（a）：）：方法方法(b)：算术平均算术平均几何平均几何平均调和平均调和平均 计计算算时时，先先用用显显式式线线性性化化的的方方法法，采采用用隐隐式式差差分分格格式式，求求出出k+1/2时时刻刻各各节节点点的的 m或或（即即 or ）；根根据据这些值可求得相应的这些值可求得相应的 or 等等参参数数，以以此此作作为为计计算算时时始始末末（tk to tk+1）的的平平均均值值，代代入方程求解入方程求解 ，从而完成方程的线性化工作。，从而完成方程的线性化工作。7.

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

13、数值，然后对参数先求出计算时段始末的参数值，然后对参数取平均。取平均。平均值的计算可采用算术平均、几何平均、调和平平均值的计算可采用算术平均、几何平均、调和平均等方法。均等方法。7.5.17.5.1线性化方法线性化方法7.5.1.27.5.1.2预报校正法预报校正法（2）线性外推法（已知线性外推法（已知 时刻，外推时刻，外推 时刻）时刻）假定含水率在相邻时段内是线性变化的，根据前假定含水率在相邻时段内是线性变化的，根据前一个时段末的含水率，用直线外推近似求出计算时段一个时段末的含水率，用直线外推近似求出计算时段末的含水率：末的含水率：7.5.17.5.1线性化方法线性化方法7.5.1.27.5

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

15、参数值（预报值），再求解隐式差分方程，得到时段末的（预报值），再求解隐式差分方程，得到时段末的含水率的值（校正值）。含水率的值（校正值）。计计算算时时，先先假假定定本本时时段段末末（tk+1）的的 m或或（一一般般做做法是先取时段初的函数值法是先取时段初的函数值 or ），or 等参数的迭代初值等参数的迭代初值 求解线性代数方程组得求解线性代数方程组得 or 直至前后两次迭代值的误差满足精度要求为止。直至前后两次迭代值的误差满足精度要求为止。7.5.1.3迭代法迭代法（a）若取绝对误差，判断条件为：若取绝对误差，判断条件为：or7.5.17.5.1线性化方法线性化方法7.5.1.3迭代法迭代法

16、（b）若取相对误差，判断条件为：若取相对误差，判断条件为：or7.5.17.5.1线性化方法线性化方法7.5.1.3迭代法 1与与 2 的取值可根据实际问题的精度要求而定。的取值可根据实际问题的精度要求而定。此此法法计计算算稍稍复复杂杂，但但精精度度较较高高，数数值值计计算算时时经经常常采采用用。为为避避免免迭迭代代次次数数太太多多，一一般般可可设设置置一一个个判判断断语语句句：当当迭迭代代次次数数m5时时，将将时时间间步步长长 减半，再重新计算。减半，再重新计算。7.5.17.5.1线性化方法线性化方法7.5.1.3迭代法7.5 线性化方法与土壤水分参数的取值线性化方法与土壤水分参数的取值7

17、.5.2.1直接法直接法先求出两结点的参数，然后先求出两结点的参数，然后(1)取两结点参数的算术平均值：取两结点参数的算术平均值：7.5.2参数取值参数取值(2)取两结点参数的调和平均值：取两结点参数的调和平均值：(3)取两结点参数的几何平均值：取两结点参数的几何平均值：7.5 线性化方法与土壤水分参数的取值线性化方法与土壤水分参数的取值7.5.2 2间接法间接法先求出两结点处含水率的算术、调和、几何平均值，先求出两结点处含水率的算术、调和、几何平均值，然后再求出相应的参数值。然后再求出相应的参数值。(1)算术平均值：算术平均值：(2)调和平均值：调和平均值：(3)几何平均值：几何平均值：7.

18、5 线性化方法与土壤水分参数的取值线性化方法与土壤水分参数的取值7.5.2 3其他方法其他方法（1）求出的参数和两结点处的参数和，三点参）求出的参数和两结点处的参数和，三点参数取平均，故可称为数取平均，故可称为“三点式三点式”：（2）（3）由邻近结点的参数代替，如）由邻近结点的参数代替，如7.6 7.6 边界条件的处理及追赶法求解三边界条件的处理及追赶法求解三对角方程对角方程前前已已叙叙及及，采采用用显显式式差差分分格格式式可可直直接接求求解解，采采用用隐隐式式或或中中心心差差分分格格式式时时总总可可以以得得到到如如下下形形式式的的三三对对角方程组（以角方程组（以 方程为例）：方程为例）：共共

19、n-2个个方方程程，再再加加上上上上、下下边边界界条条件件，即即可可组组成成n n 阶阶的的n元元一一次次方方程程组组，线线性性化化后后求求解解该该方方程程组组可可获获得得原定解问题的解。原定解问题的解。7.6.17.6.1边界条件的处理边界条件的处理 上、下边界不论为第一、第二或第三类边界，上、下边界不论为第一、第二或第三类边界，总可以化为：总可以化为：(2)(1)7.6 7.6 边界条件的处理及追赶法求解三边界条件的处理及追赶法求解三对角方程对角方程7.6.17.6.1边界条件的处理边界条件的处理 验验 证：证：当上、下边界为第一类边界时当上、下边界为第一类边界时=+)(,0,1)2.(E

20、qn)(,0,1)1.(Eqn12n11111knnktfHEFtfHGF7.6 7.6 边界条件的处理及追赶法求解三边界条件的处理及追赶法求解三对角方程对角方程7.6.17.6.1边界条件的处理边界条件的处理 当上、下边界为第二、三类边界时当上、下边界为第二、三类边界时=+-=+-+-+-+-+)(),(122/111112111112312111223knnnknknnkkkktqKzzDtqKzzDqqqqq经线性化后，原定解问题可化为经线性化后，原定解问题可化为n n阶线性代数方程组，即：阶线性代数方程组，即：以矩阵的形式表示即为：以矩阵的形式表示即为：=-+-+-nnknknkkkn

21、nnnnHHHHHFEGFEGFEGFEGF132111113121111133322211LLLLLLLLLLqqqqq7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程1原理原理(1)(1a)7.6 7.6 边界条件的处理及追赶法求解三边界条件的处理及追赶法求解三对角方程对角方程7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程1原理原理(2)(2a)7.6 7.6 边界条件的处理及追赶法求解三边界条件的处理及追赶法求解三对角方程对角方程.(ia)(i)7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程(n-1)a)(n-1)7.6.27.6.2 追赶法

22、求解三对角方程追赶法求解三对角方程7.6 7.6 边界条件的处理及追赶法求解三边界条件的处理及追赶法求解三对角方程对角方程(na)(n)7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程7.6 7.6 边界条件的处理及追赶法求解三边界条件的处理及追赶法求解三对角方程对角方程 解出解出 后，利用式（后，利用式（(n-1)a）,，(ia),，(2a),(1a)往前递推，依次即可解出往前递推，依次即可解出 ，于是得到，于是得到k+1时刻各节点的时刻各节点的 值（值（i=1,2,，n）。）。若若 均为已知（第一类边界条件），则可直均为已知（第一类边界条件），则可直接从式（接从式（2）开始求

23、解，依次递推至式（）开始求解，依次递推至式（n-1）即可。即可。7.6.27.6.2 追赶法求解三对角方程追赶法求解三对角方程7.6 7.6 边界条件的处理及追赶法求解三边界条件的处理及追赶法求解三对角方程对角方程7.7 7.7 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图(1)输入输入D、K含水量含水量SITA关系（建议采用关系（建议采用函数函数子程序调用）；子程序调用）；节节点点剖剖分分总总数数n及及对对应应的的空空间间坐坐标标z(1:n)（若若为为等等步步长长，则则输入空间步长输入空间步长DZ）；）；初始含水量剖面分布初始含水量剖面分布SITA0（1:n）；）；上上、

24、下下边边界界条条件件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=17.7 7.7 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图(

25、3)判断计算时间是否已超过需模拟时间判断计算时间是否已超过需模拟时间if(T(KT)=DTMAX)then DT=DTMAX T(KT)=T(KT-1)+DT endif7.7 7.7 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图(5)开始计算系数矩阵，求解线性代数方程组开始计算系数矩阵，求解线性代数方程组 赋迭代计次数赋迭代计次数 MP1=07.7 7.7 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图(6)if(MP1=5)then (else goto(12)(7)计算本时刻相应的边界值计算本时刻相应的边界值QSUBC(T(KT)与与QSLBC(T(

26、KT)*（以下计算三对角方程系数（以下计算三对角方程系数（i=1,n）F(1),G(1),H(1)E(n),F(n),H(n)(8)计算三对角方程组各系数计算三对角方程组各系数 do i=2,n-1由由DT,SITAM(1:n)及及SITA0（1:n）求求 D，K及及E(i),F(i),G(i),H(i)enddo7.7 7.7 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图(9)求解三对角方程得求解三对角方程得SITA(1:n)(10)if then (else goto(11)MP1=MP1+1 do i=1,n SITAM(i)=SITA(i)enddo goto(6)

27、7.7 7.7 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图(11)else（输出结果，进入下一时段计算）输出结果，进入下一时段计算）print SITA(1:n)KT=KT+1T(KT)=2.25*T(KT-1)-1.25*T(KT-2)do i=1,n SITA0(i)=SITA(i)enddogoto(3)endif(10)7.7 7.7 垂直一维非饱和土壤水流计算流程图垂直一维非饱和土壤水流计算流程图(12)else（迭迭代代次次数数超超过过5次次时时，将将时时间间步步长长减减半半，再重新计算）再重新计算）MP1=0 DT=0.5*DT T(KT)=T(KT-1)+

28、DT do i=1,n SITAM(i)=SITA0(i)赋迭代初值赋迭代初值 enddo goto(6)(13)else（计算时间超过模拟时间）计算时间超过模拟时间）stop end endif(3)第第第第7 7章章章章 作作作作 业业业业1、用差分法求解、用差分法求解取取x=2，对对0 x10进行剖分（手算，有步骤）。进行剖分（手算，有步骤）。2、算至稳定为止（先手算，然后编程计算）。算至稳定为止（先手算，然后编程计算）。3、Philip解作业解作业4、Warrick,A.W.,et al.Simultaneous solute and water transport for an unsaturated soil.W.R.R.,7(5):1216-1225,1971.