近些年来, 定量研究溶质通过非饱和土壤的运移, 已成为土壤物理中的热点领域之一.由于非饱和土壤的非均质性和优势流等因素的影响, 溶质在其中的运移往往是对流占优势的.有关对流占优问题(大网格Peclet数) 的数值算法, 在地下水的研究中已有许多例子, 但在非饱和土壤中的研究还不多见.通量校正运移(flux-corrected transport) 方法是由Boris等^[1]和Book等^[2]提出的, 其后, Zalesak^[3]、Hills等^[4]在这方面又做了很多研究.FCT方法是一种混合有限差分格式, 它利用了低阶近似和高阶校正格式, 这种方法可以减少数值弥散和数值振荡误差.本文是在Hills研究的基础上, 结合用FUCG算法解Richards方程, 以数值例子验证了FCT方法的有效性.

1.   方程的分解

考虑溶质在非饱和土壤中的二维运移:

式中, c是溶质的浓度, θ是水体积分数, q_x, q_z为x方向和z方向上的通量, D_xx, D_zz分别是x方向和z方向上的水动力弥散系数.

边界条件是:

式中, α, β是参数, α=1, β=0说明是第一类边界条件; α=0, β≠0则是第二类边界条件; n_j是方向指向Γ外侧并与γ垂直的单位向量; j表示x或z; θ是水体积分数.

把(1) 式中的时间导数项用差分近似得出:

引进中间变量(θc)^M, 把(3) 式分成两部分:

式(4), (5) 组成一纯对流问题, 可用FCT方法来解.把式(4), (5) 代入式(2), (3) 得

(6), (7) 式构成了一纯弥散问题, 可以用常规的交替方向隐式差分方法来求解.

2.   二维FCT算法

对于二维溶质运移问题, (4) 式可以写成:

式中, f=q_xc, g=q_zc.

FCT算法包括两个主要阶段: 低阶对流或运移阶段和高阶反扩散或校正阶段.

2.1   低阶对流阶段的解

低阶对流阶段的解可通过预测—校正算法来实现.下列各式中的和分别是(i, j) 和(i+1, j) 的算术平均值, 以及(i, j) 和(i, j+1) 的算术平均值.

预测步:

校正步:

(10), (12) 是两个显式差分方程, 可方便地求出未知值c_{i, j}^*和c_{i, j}^L.

2.2   高阶校正阶段的解

在第一阶段的解中, 可能会引进过量的数值弥散, 因此, FCT算法的第二阶段是通过反扩散通量来校正第一阶段的解.

首先, 由下列步骤求出高阶通量和:

高阶通量表示为:

那么, 浓度的临时值(即预测值) c_{i, j}^p可用下面的方程显式求出:

用c_{i, j}^p来修正通量f_{i, j}和g_{i, j},

这样, 再用下列方程得出校正的高阶通量:

其次确定反扩散通量和:

假定c_{i, j} > c_{i+1, j}, 反扩散通量的作用总是趋于增大c_{i, j}, 而减小c_{i+1, j}.反扩散通量的作用与此类似.这种效应与粘滞性的作用刚好相反, 粘滞性使数值解趋向比较均一.

如果把反扩散通量直接加到低阶通量中去并用它来解方程, 会使方程的解产生振荡.为了消除这种振荡, 必须对反扩散通量加以校正(限制).具体做法如下:

同样地,

在完成了上述情况下反扩散通量和的校正后, 还须对其他情况下的反扩散通量进行校正.校正后的反扩散通量为:

式中, 为校正因子, 它们都是小于1的正数.

最后, 用校正后的反扩散通量求得低阶中间浓度c_{i, j}^Μ:

至此, FCT算法就算完成了.

3.   整个问题的解

当用(21) 式求得中间值c_{i, j}^Μ后, 把c_{i, j}^Μ代入(6) 式, 并结合其边界条件(7) 式, 用隐式差分解法即可求得c_{i, j}ⁿ⁺¹, 这样就完成了一个时段Δt的计算, 连续进行下去, 就可达到解得所指定的某一时刻t时的浓度的目的.

即在x方向:

在z方向上:

(22), (23) 式都是三对角方程, 可用追赶法求解.为了保持质量守恒, 我们用Kirkland等^[5]提出的新通量迭代共轭梯度(flux-updating iterative conjugate gradient, FUCG) 算法来求水体积分数, 方法如下:

以压力形式表示的Richards方程可写为:

用FUCG算法对(24) 进行差分近似:

式中, C=∂θ/∂h为比水容量, h^*是新的时间步长对h的估计. (25) 式可用交替方向隐式差分法求解.

h_{i, j}^*求出后, 通量用下列方程计算:

则水体积分数

4.   数值例子

考虑1个100 cm×100 cm的垂向(X-Z) 土壤剖面, 该土壤的非饱和水力传导率K (h) 适合van Genechten-Mualem模型.模型中的参数K_s, θ_s, θ_r, α和n采用Tseng等^[6]给出的值, 即K_s=0.3 cm/min, θ_r=0.05, θ_s=0.4, α=0.019, n=1.59.另外, 土壤的初始水体积分数给定为: θ (x, z, 0) =0.3, 弥散系数为D=0.01 v (cm²/min), v为水分的运动速率.

假定在研究的土壤剖面顶部连续注入质量浓度为1 mg·L^-1的溶质1 h, 那么, 水分和溶质运移问题的初边值边界条件分别为:

取空间步长Δx=Δz=10 cm, 时间步长Δt=60 min, 用前面提到的方法解上述定解问题, 则在剖面的任一垂直方向上, 溶质随深度的变化如图 1所示.

图  1  质量浓度随深度的变化曲线

Fig.  1.  Variation curves of mass concentration with depth

下载: 全尺寸图片 幻灯片

从图 1可以看出, 在t=2 d, 5 d, 10 d, 20 d和30 d的质量浓度锋面都近于陡立, 说明用FCT方法解非饱和土壤中溶质运移问题是可行的.

5.   结论

用FCT算法解对流占优(大网格Peclet数) 问题, 实质上是把原对流弥散问题分裂成对流和弥散二部分, 对流部分用FCT法模拟, 弥散部分用常规的交替方向隐式差分法求解.在一个时段内, 用FCT法求出的浓度作为模拟弥散部分时的初值.而解得弥散部分的值, 即为整个问题在该时段内的解.尽管推导FCT算法的过程比较烦琐, 但由于各式都是一些显式方程, 并不需要通过形成方程组来求解, 求解非常简单, 且容易编程, 提高了计算效率.数值例子表明FCT算法可有效地消除数值弥散和数值振荡.另外, 用FUCG方法解Richards方程, 得到水体积分数, 保持了质量守恒, 为下一步用FCT算法解对流弥散问题做了有效的保证.