1.5 土壤水分运动基本方程的求解方法

1.5.1 求解方法的总体特征

由于土壤水分运动基本方程是一个偏微分方程，目前仍无法获得其精确解析解，通常采用数值计算和半解析解或近似方法进行求解。目前常用的数值计算方法有有限差分、有限元、边界元、无网格法等方法。

就半解析解或近似分析方法而言，人们通过概化土壤水分运动或水分含量或能量分布特征，建立相应分析方法，目前分析方法主要有以下几种类型。

（1）从土壤水分运动基本方程入手，借助数理方程的基本理论，求解土壤水分运动基本方程。目前常用Boltzmann变换求解水平一维土壤水分运动基本方程，获得土壤水分运动基本方程的半解析解，并发展了土壤水分扩散率测定方法及Philip入渗模型。

（2）从土壤水分运动参数入手，通过假定土壤水分运动参数的变化特征，简化土壤水分运动基本方程。通常假定非饱和土壤导水率是含水量线形函数，或者假定土壤水分扩散系数或导水率是一个特定函数，将土壤水分运动方程的线形化，然后利用数理方程理论获得相应的解析解。

（3）从土壤水分通量入手，通过假定土壤水分通量的变化特征来简化土壤水分运动基本方程，获得土壤水分运动基本方程的近似分析方法。一般假定任一时刻土壤剖面水分通量相等，这样将土壤水分运动基本方程进行简化，获得相应的分析解。如著名的Green-Ampt入渗公式推求过程可以属于这一类型。

（4）从土壤含水量或吸力分布入手，通过假定土壤含水量或吸力剖面分布，进而求解土壤水分运动基本方程。邵明安和王全九通过假定土壤吸力分布，分别获得了推求van Genuchten和Brooks-Corey模型参数的解析表达式，进而建立了相应的简捷方法。

1.5.2 有限差分解算方法

非饱和土壤水分运动方程属于抛物型非线性偏微分方程。目前，在土壤水分运动计算中，多采用直接差分法进行偏微分方程的离散化，其中控制容积法已比较成熟，其推导过程物理概念清晰，且可以保证离散化方程具有守恒特性。以下以直角坐标系下的土壤水分运动的三维和柱坐标系下二维情况为例，说明非饱和土壤水分运动方程离散化的控制容积法。土壤水分一维、二维运动方程的离散化与求解，在此基础上直接简化即可。

1.5.2.1 直角坐标系下的土壤水分运动方程离散化

对计算区域按长方体进行网格剖分，即图1.6所示的实线网格。对于网格节点P，与其相邻的网格节点有E、W、S、N、U、V，取包含P点的控制容积，虚线e、w、s、n、u、v表示控制容积面。

图1.6 直角坐标系下的网格剖分示意图

（a）xy面；（b）xz面；（c）yz面

图1.6中（δx）_e、（δx）_w、（δy）_u、（δy）_v、（δz）_s、（δz）_n为网格间距，对网格采用均匀剖分，即（δx）_e=（δx）_w=Δx，（δy）_u=（δy）_v=Δy，（δz）_s=（δz）_n=Δz。x方向的节点编号为i=1，2，…，nx，y方向的节点编号为j=1，2，…，ny，z方向的节点编号为k=1，2，…，nz。

采用控制容积法对非饱和土壤水分运动基本方程式（1.10）在控制容积中进行积分，时间间隔为［t+Δt］，计算过程如下：

写为差分形式则为

（1）非稳态项。对等式（1.88）左边在时间间隔［t+Δt］内进行积分，取θ在空间变化的型线为阶梯式，即同一控制容积中各处的θ值相同，等于节点P上之值θ_P，故有

（2）扩散项。选取导数随时间作显式阶跃式的变化，得

进一步取θ随x、y、z呈分段线性变化，则式 （1.90）～式 （1.92）中）^、）可表示为

由于网格剖分采用均匀网格，所以上述各式中有（δx）_e=（δx）_w=Δx，（δy）_v=（δy）_u=Δy，（δz）_s=（δz）_n=Δz。

（3）源项。对等式 （1.88）右边最后一项）_{按照源项进行处理，对该项随x、}y、t取阶梯式型线，则

将离散化后的各项代入式（1.88），整理后得隐格式的离散化方程为

简化之，可以写为如下形式

式中：；S=-（k_s-k_n）ΔxΔy；ΔV=ΔxΔyΔz；a_P=+a_E+a_W+a_V+a_U+a_S+a_N；θ_P、θ_E、θ_W、θ_S、θ_N、θ_U、θ_V分别为各节点的未知含水率 （即t+Δt时刻的含水率）；为P点在前一个时刻的含水率 （即t时刻的含水率）；D_e、D_w、D_s、D_n、D_u、D_v分别为控制容积面e、w、s、n、u、v处的扩散率；K_s、K_n分别为控制容积面s、n处的导水率；Δx、Δy、Δz分别为网格间距。

由离散化方程式（1.95）可看出，节点P在t+Δt时刻的含水率是其上一时刻的含水率及其相邻各点在本时刻的含水率θ_P、θ_E、θ_W、θ_S、θ_N、θ_U、θ_V以不同权重叠加的结果。离散化方程式 （1.95）含有7个未知变量，因此对整个计算区域离散化，可得到每行含有7个未知元素的离散化方程组。

1.5.2.2 柱坐标系下的土壤水分运动方程离散化

对计算区域按矩形网格进行剖分，如图1.7所示的实线网格。对于网格节点P，与其相邻的网格节点有E、W、S、N。取包含P点的控制容积，如图1.7中的阴影部分，虚线e、w、s、n表示控制容积面。图中（δr）_e、（δr）_w、（δz）_s、（δz）_n为网格间距，如果网格为均匀剖分，则（δr）_e=（δr）_w=Δr，（δz）_s=（δz）_n=Δz。为简单起见，以下按均匀网格处理。

图1.7 柱坐标系下的网格剖分示意图

对于控制方程式 （1.15）右边的最后一项 [即按源项处理，则通过在控制容积中对控制方程进行积分，其中时间间隔由t到t+Δt，可得全隐格式的离散化方程为

式中：；；a_P=a_E+a_W+a_S+a_N+；ΔV=+r_W）ΔrΔz，（ΔV为控制容积的体积）；θ_P、θ_E、θ_W、θ_S、θ_N分别为各节点的未知含水率（即t+Δt时刻的含水率）；为P点在前一个时刻的含水率 （即t时刻的含水率）；D_e、D_w、D_s、D_n分别为控制容积面e、w、s、n处的扩散率；K_s、K_n分别为控制容积面s、n处的导水率；r_E、r_W、r_P分别为节点E、W、P处的径向坐标。

由离散化方程式（1.96）可以清晰地看出，节点P在t+Δt时刻的含水率是其上一时刻的含水率及其相邻各点在本时刻的含水率θ_E、θ_W、θ_S、θ_N以不同的权重叠加的结果。离散化方程式 （1.96）含有5个未知变量，因此对整个计算区域离散化，可得每行含有5个未知元素的离散化方程组。

1.5.2.3 非线性系数的处理

以上述柱坐标系下的离散化方程为例，控制容积面e、w、s、n处的扩散率D_e、D_w、D_s、D_n和控制容积面s、n处的导水率K_s、K_n，要利用该控制容积面两侧节点在t+Δt时刻的值进行计算，一般可取两点的算术平均值。

扩散率D（θ）和导水率K（θ）都是含水率的函数，然而节点E、W、S、N在t+Δt时刻的含水率是待求的未知量。对于该非线性问题可以采用显式线性化法、预报校正法、迭代法等来处理。处理非线性问题最常用的方法是对每一个时间步的计算采用迭代法来进行处理。先以各点在t时刻的含水率θ^t计算其扩散率D（θ）和导水率K（θ），从而得到所需的控制容积面e、w、s、n处的扩散率D_e、D_w、D_s、D_n和控制容积面s、n处的导水率K_s、K_n，通过对离散化方程的解算，可以得到各界点含水率在t+Δt时刻的第1次迭代结果θ^t+Δt；再利用该结果计算D_e、D_w、D_s、D_n、K_s、K_n等值，进行第2次迭代……直至前后两次迭代计算出的含水率之差小于某一规定值。

1.5.2.4 离散化方程组的求解

求解方程组的方法可以分为直接法和迭代法两类。对于该非线性问题，离散方程中的系数也是未知量θ的函数，这样整个问题的求解必然是迭代性质的，经过多次迭代直到获得收敛的解。在这一计算过程中，每一次求解代数方程时，其系数都是临时的，若采用直接方法求解，则所得到的是关于这一组临时系数的解，但既然代数方程本身的系数是有待改进的，就没有必要把相应的真解求出来；若采用迭代法，则可控制在适当时候中止迭代，以在改进代数方程系数后再求解，故对非线性问题来说，直接解法也是不经济的。迭代法包括点迭代法、块迭代法和交替方向隐式迭代法（ADI方法），其中ADI方法计算简单有效，在土壤入渗计算中常被采用。

1.直角坐标系下的土壤水分运动

ADI方法的具体实施方式很多，最简单的就是采用Jacobi方式的按行与列的交替迭代，对于离散化方程式（1.95）可表示为

在k方向上：

在i方向上：

在j方向上：

在一个时间步长的计算过程中，先按式（1.97）在k方向上扫描，结合边界点的边界条件得三对角系数矩阵的方程组，按追赶法解算，得到的含水率作为对j方向扫描的初始含水率；在k方向扫描完毕后，按式（1.98）对i方向进行扫描计算，得到的含水率作为对j方向扫描的初始含水率；最后按式（1.99）对j方向进行扫描计算，完成一次完整扫描。对i、j、k三个方向进行交替扫描，直至前后两次扫描结果之差满足精度要求为止。然后逐步进行下一步长的计算。

在计算过程中，对于列或行的扫描方向（即，在按列扫描的情况下，是从左向右逐列扫描，还是从右向左逐列扫描；在按行扫描的情况下，是从上向下逐行扫描，还是从下向上逐行扫描）是一个影响收敛速度的重要问题。原则上是从对整个计算区域影响较大的边界处开始。

时间步长与空间步长是影响计算精度的重要因素。若步长取得过大，则会引起较大误差，过小又会占用较大内存、花费较长的计算时间。在对非线性偏微分方程进行数值求解时，目前尚无成熟的理论来选择适宜的时间步长与空间步长。对一个具体的计算问题，一般采用试取的方法来确定其时间步长与空间步长。

2.柱坐标系下的土壤水分运动

设通过对计算区域的网格剖分，在z方向共有m条网格线，在r方向共有n条网格线，所以共有网格节点m×n个。于是，通过离散化方程式（1.96）并考虑边界条件，可得含有（m×n-m-n+4）个未知变量的方程组，该方程组的系数矩阵为大型稀疏矩阵，对其直接进行解算（即直接解法），占用计算机内存多，计算较复杂。

对于一维抛物型偏微分方程，其全隐格式的离散化方程组的系数矩阵均为三对角矩阵，采用追赶法求解非常方便。对于高维抛物型偏微分方程采用某些迭代格式，便可以对每一个时间步通过多次求解三对角方程组来得到问题的解（迭代解法）。其中交替方向迭代法（ADI方法）计算简单有效，在土壤水分运动计算中常被采用。

ADI方法的具体实施方式很多，最简单的就是采用Jacobi方式的按列与行交替迭代，对于离散化方程式（1.96）可表示为

式中：带有上标t的量为在时刻t的值；带有上标t+Δt/2的量为在时刻t与时刻t+Δt的中间值，在式（1.100）中为未知变量，在式（1.101）中为已知变量；带有上标t+Δt的量为在时刻t+Δt的值。

在计算过程中，先对列按式（1.100）列方程组，得三对角系数矩阵的方程组，按追赶法解算，图1.8为按列扫描的图像。逐列扫描完毕后，再对行按式（1.101）进行逐行扫描计算，于是完成一次完整扫描。对行、列交替进行扫描，直至前后两次扫描结果之差满足精度要求为止。

图1.8 按列扫描的图像