《第7章土壤水分运动的数值解课件.ppt》由会员分享，可在线阅读，更多相关《第7章土壤水分运动的数值解课件.ppt（58页珍藏版）》请在三一办公上搜索。
1、第7章 土壤水分运动的数值解,7.1 概述7.2 差分的概念及分类7.3 显式差分、隐式差分和中心差分7.4 土壤水分运移方程的线性化方法7.5 线性化方法与土壤水分参数的取值7.6 边界条件的处理及追赶法求解三对角方程7.7 垂直一维非饱和土壤水流计算流程图,第7章 土壤水分运动的数值解,土壤水分运动方程的求解方法：解析解（Analytical solution）数值解（Numerical solution）(1)有限差分法（Finite Difference method）(2)有限单元法（Finite Element method）,7.1 概述,7.2 差分的概念及分类,差分的由来（T
2、aylor展开）if=(x)exists,then,第7章 土壤水分运动的数值解,when x 0,then,向前差分,向后差分,中心差分（Eqn.(1)-Eqn.(2)）,i.e.,如图所示：BC表示向前差分；AB表示向后差分；AC表示中心差分；,In addition,Eqn.(1)+Eqn.(2),i.e.:,其截断误差亦为与x2同阶的无穷小量。,考虚垂直一维问题（如右图所示），泛定方程可写为（以 方程为例）：,剖面深度为L，共剖分出n个节点（i=1,2,3,n），z(1)=0,z(n)=L，已知t=tk时刻各节点含水率 的分布，求t=t k+1时刻各节点含水率 的分布，,7.3显式差分
3、、隐式差分和中心差分,7.3.1显式差分格式,7.3显式差分、隐式差分和中心差分,7.3.1显式差分格式,7.3显式差分、隐式差分和中心差分,7.3.1显式差分格式,7.3显式差分、隐式差分和中心差分,若取等步长z，即z=zi+1-zi=zi-zi-1，则有：,7.3.1显式差分格式,7.3显式差分、隐式差分和中心差分,显式差分格式是有条件收敛的，一般应满足,其中：r多取1/2。由于土壤接近饱和时，Dmax很大，故t一般要求取值很小，耗费机时。,7.3.1显式差分格式,7.3显式差分、隐式差分和中心差分,7.3.2隐式差分格式,7.3显式差分、隐式差分和中心差分,Let,then,(i=2,3
4、,n-1),7.3.2隐式差分格式,7.3显式差分、隐式差分和中心差分,where,隐式差分格式是无条件收敛的。,7.3.2隐式差分格式,7.3显式差分、隐式差分和中心差分,7.3.3中心差分格式,7.3显式差分、隐式差分和中心差分,7.3.3中心差分格式,7.3显式差分、隐式差分和中心差分,Let,then,(i=2,3,n-1),where,7.3.3中心差分格式,7.3显式差分、隐式差分和中心差分,中心差分格式也是无条件收敛的。,7.3.3中心差分格式,7.3显式差分、隐式差分和中心差分,7.4 土壤水分运移方程的线性化方法,以基质势水头h为因变量的一维垂直向土壤水分运动基本方程为：,7
5、.4 土壤水分运移方程的线性化方法,7.4.1显式差分格式相应h方程的差分方程为：,7.4 土壤水分运移方程的线性化方法,7.4.2隐式差分格式相应h方程的差分方程为,7.4 土壤水分运移方程的线性化方法,7.4.3中心差分格式相应h方程的差分方程为,7.5 线性化方法与土壤水分参数的取值,土壤水分运移方程中的各参数均依赖于变量 或m，从而使得方程呈现出较强的非线性性，求解前必须将其线性化，得到n元一次线性代数方程组，以便于求解。（1）显式线性化 计算过程中，方程中的各参数如：C(m)，K(m)or D()，K()等均以时段初的值（即前一时段的值）代入。适用条件：此法只能用于m或 变化缓慢的情
6、况；当变化剧烈时，此法会导致较大的偏差。,7.5.1线性化方法,例如：已知k-1与k时刻各节点的 值，需求k+1时刻各节点的 值，应用显式外推线性化参数的方法如下（以导水率K的求法为例，其他各参数的求法与此类似）：,(*),求：方法（a）：,方法(b)：算术平均:几何平均:调和平均:,求：方法（a）：,方法(b)：算术平均几何平均调和平均,计算时，先用显式线性化的方法，采用隐式差分格式，求出k+1/2时刻各节点的m或（即 or）；根据这些值可求得相应的 or 等参数，以此作为计算时始末（tk to tk+1）的平均值，代入方程求解，从而完成方程的线性化工作。,7.5 线性化方法与土壤水分参数的
7、取值,7.5.1线性化方法,7.5.1.2预报校正法,（1）显式外推法（先求显式差分方程，得 作为预报值）计算参数时，据前一时段始末的函数值m或，用线性外推近似求出本时段末的函数值，然后,方法(a)：求出计算时段始末的平均值，再通过各参数与m或 的关系，求出各参数；或 方法(b)：先求出计算时段始末的参数值，然后对参数取平均。,平均值的计算可采用算术平均、几何平均、调和平均等方法。,7.5.1线性化方法,7.5.1.2预报校正法,（2）线性外推法（已知 时刻，外推 时刻）假定含水率在相邻时段内是线性变化的，根据前一个时段末的含水率，用直线外推近似求出计算时段末的含水率：,7.5.1线性化方法,
8、7.5.1.2预报校正法,对于第一时段，无法应用线性外推法，可用其他线性化方法。,7.5.1线性化方法,7.5.1.2预报校正法,（3）显式预测法（两次求解隐式：第一次预报值；第二次得校正值）,首先利用显式线性化方法，由时段初的含水率及相应的参数值，按隐式差分格式求解方程，求得本时段末（或时段中间）的含水率及相应的参数值（预报值），再求解隐式差分方程，得到时段末的含水率的值（校正值）。,计算时，先假定本时段末（tk+1）的m或（一般做法是先取时段初的函数值 or），or 等参数的迭代初值 求解线性代数方程组得 or 直至前后两次迭代值的误差满足精度要求为止。,7.5.1.3迭代法,（a）若取绝
9、对误差，判断条件为：,or,7.5.1线性化方法,7.5.1.3迭代法,（b）若取相对误差，判断条件为：,or,7.5.1线性化方法,7.5.1.3迭代法,1与2 的取值可根据实际问题的精度要求而定。此法计算稍复杂，但精度较高，数值计算时经常采用。为避免迭代次数太多，一般可设置一个判断语句：当迭代次数m5时，将时间步长 减半，再重新计算。,7.5.1线性化方法,7.5.1.3迭代法,7.5 线性化方法与土壤水分参数的取值,7.5.2.1直接法先求出两结点的参数，然后(1)取两结点参数的算术平均值：,7.5.2参数取值,(2)取两结点参数的调和平均值：,(3)取两结点参数的几何平均值：,7.5
10、线性化方法与土壤水分参数的取值,7.5.2 2间接法先求出两结点处含水率的算术、调和、几何平均值，然后再求出相应的参数值。,(1)算术平均值：,(2)调和平均值：,(3)几何平均值：,7.5 线性化方法与土壤水分参数的取值,7.5.2 3其他方法（1）求出的参数和两结点处的参数和，三点参数取平均，故可称为“三点式”：,（2）,（3）由邻近结点的参数代替，如,7.6 边界条件的处理及追赶法求解三对角方程,前已叙及，采用显式差分格式可直接求解，采用隐式或中心差分格式时总可以得到如下形式的三对角方程组（以 方程为例）：共n-2个方程，再加上上、下边界条件，即可组成n n 阶的n元一次方程组，线性化后
11、求解该方程组可获得原定解问题的解。,7.6.1边界条件的处理,上、下边界不论为第一、第二或第三类边界，总可以化为：,(2),(1),7.6 边界条件的处理及追赶法求解三对角方程,7.6.1边界条件的处理,验 证：当上、下边界为第一类边界时,=,=,=,=,=,=,+,+,),(,0,1,),2,.(,Eqn,),(,0,1,),1,.(,Eqn,1,2,n,1,1,1,1,1,k,n,n,k,t,f,H,E,F,t,f,H,G,F,7.6 边界条件的处理及追赶法求解三对角方程,7.6.1边界条件的处理,当上、下边界为第二、三类边界时,=,+,-,-,-,=,+,-,-,-,+,-,-,+,-,
12、+,-,+,+,+,+,),(,),(,1,2,2,/,1,1,1,1,1,2,1,1,1,1,1,2,3,1,2,1,1,1,2,2,3,k,n,n,n,k,n,k,n,n,k,k,k,k,t,q,K,z,z,D,t,q,K,z,z,D,q,q,q,q,q,经线性化后，原定解问题可化为n n阶线性代数方程组，即：,以矩阵的形式表示即为：,=,-,+,+,-,+,+,+,-,-,-,n,n,k,n,k,n,k,k,k,n,n,n,n,n,H,H,H,H,H,F,E,G,F,E,G,F,E,G,F,E,G,F,1,3,2,1,1,1,1,1,3,1,2,1,1,1,1,1,3,3,3,2,2,2
13、,1,1,L,L,L,L,L,L,L,L,L,L,q,q,q,q,q,7.6.2 追赶法求解三对角方程,1原理,(1),(1a),7.6 边界条件的处理及追赶法求解三对角方程,7.6.2 追赶法求解三对角方程,1原理,(2),(2a),7.6 边界条件的处理及追赶法求解三对角方程,.,.,(ia),(i),7.6.2 追赶法求解三对角方程,(n-1)a),(n-1),7.6.2 追赶法求解三对角方程,7.6 边界条件的处理及追赶法求解三对角方程,(na),(n),7.6.2 追赶法求解三对角方程,7.6 边界条件的处理及追赶法求解三对角方程,解出 后，利用式（(n-1)a）,，(ia),，(2
14、a),(1a)往前递推，依次即可解出，于是得到k+1时刻各节点的 值（i=1,2,，n）。若 均为已知（第一类边界条件），则可直接从式（2）开始求解，依次递推至式（n-1）即可。,7.6.2 追赶法求解三对角方程,7.6 边界条件的处理及追赶法求解三对角方程,7.7 垂直一维非饱和土壤水流计算流程图,(1)输入D、K含水量SITA关系（建议采用函数子程序调用）；节点剖分总数n及对应的空间坐标z(1:n)（若为等步长，则输入空间步长DZ）；初始含水量剖面分布SITA0（1:n）；上、下边界条件UBC、LBC（=1，2，or3）及相应的边界值QSUBC、QSLBCT（建议采用子程序调用）；初始时间
15、步长DT0、最大时间步长DTMAX，计算结束时间TEND，迭代控制标准DET。,(2)赋初值T(0)=0 初始时间T(1)=DT0 第一次计算时间*（以下赋含水量计算中间值）do i=1,nSITAM(i)=SITA0(i)enddo时段计数KT=1,7.7 垂直一维非饱和土壤水流计算流程图,(3)判断计算时间是否已超过需模拟时间if(T(KT)=TEND)then(Note:else goto(13),(4)DT=T(KT)-T(KT-1)*(以下判断时间步长是否超过允许最大时间步长)if(DT=DTMAX)then DT=DTMAX T(KT)=T(KT-1)+DT endif,7.7 垂
16、直一维非饱和土壤水流计算流程图,(5)开始计算系数矩阵，求解线性代数方程组 赋迭代计次数 MP1=0,7.7 垂直一维非饱和土壤水流计算流程图,(6)if(MP1=5)then(else goto(12)(7)计算本时刻相应的边界值QSUBC(T(KT)与QSLBC(T(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)enddo,7.7 垂直一维非饱和土壤水流计算流程图,(9)求解三对角方程
17、得SITA(1:n),(10)if then(else goto(11)MP1=MP1+1 do i=1,n SITAM(i)=SITA(i)enddo goto(6),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 垂直一维非饱和土壤水流计算流程图,(12)else（迭代次数超过5次时，将时间步长减半，再重新计算）MP1=0 DT=0.5*DT T(KT)=T(KT-1)+DT do i=1,n SITAM(i)=SITA0(i)赋迭代初值 enddo goto(6)(13)else（计算时间超过模拟时间）stop end endif(3),第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.,
