一维非饱和土壤水分运动模型的数值模拟

梁晨璟+李春光+赵文娟+王建华

摘要：分析有限体积法在一维非饱和土壤水分运动模型的数值算法上的适用性及优劣，通过算例对有限体积法进行验证。结果表明，有限体积法同样适用于一维非饱和土壤水分运动，与大多数学者都采用的有限差分法计算结果比较发现两种算法的数值模拟结果相近。但是改变网格疏密后，对有限体积法的计算结果影响不大且与实测值相近，而对有限差分法的计算结果有影响，从而说明在网格稀疏的情况下，有限体积法更好。

关键词：有限体积法；有限差分法；非饱和土壤；水分运动

中图分类号：S152.7； S11+9 文献标识码：A 文章编号：0439-8114（2014）15-3525-04

1-D Numerical Simulation of Water Movement in Unsaturated Soil

LIANG Chen-jing1，LI Chun-guang1，2，ZHAO Wen-juan1，WANG Jian-hua3

（1.Civil Engineering and Water Conservancy， Ningxia University， Yinchuan 750021， China； 2. Beifang University of Nationalities， Yinchuan 750000， China；3. Zhejiang Construction Supervision Co.， Ltd，Hangzhou 310012，China）

Abstract： The applicability and pros or cons of finite volume method in one-dimensional model of water movement in unsaturated soil were analyzed. The finite volume method was verified by a numerical example. The results showed that the finite volume method was appliable to one-dimensional model of water movement in unsaturated soil. Compared with the finite difference method used by most scholars， the numerical simulation results of two algorithms were similar. Changing the mesh density on the finite volume method had little effect on the calculated results similiar to the measured values， while the finite difference method calculation was affected. In the case of sparse grid， the finite volume method was better.

Key words： finite volume method；finite difference method；unsaturated soil； water movement

收稿日期：2014-03-15

基金项目：国家自然科学基金重大研究计划培养项目（91230111）

作者简介：梁晨璟（1988-），女，河北邢台人，在读硕士研究生，研究方向为旱区河流泥沙动力学理论与数值模拟，（电话）15909510629（电子信

箱）liangmeng88-12@163.com；通讯作者，李春光，男，教授，博士，主要从事计算力学和流体力学研究，（电子信箱）

cglizd@hotmail.com。

对于一维非饱和土壤的水分运动模型研究，大多数学者都是采用有限差分法来进行数值模拟。这种方法虽然形式简单，且任意偏微分方程都能写出与之对应的差分方程，但不能从差分方程中体现出各项的物理意义。因此，差分方程只能认为是对微分方程的数学近似，基本上不能反映出物理特征。差分方程在计算结果上也有可能出现不合理的情况[1]。大量的计算实践表明，用通常的差分方法求解土壤水盐运移问题，只有当弥散作用占优时，才能取得较为满意的结果。在求解弥散系数较小、流速相对较大的对流占优问题时，经常遇到数值弥散和数值振荡现象[2]。而有限体积法[3]（Finite volume methed，FVM）是以积分型守衡方程为出发点，通过对流体运动的体积域的离散来构造积分型离散方程。可以看出，这种方法便于用来模拟具有复杂边界区域流体的运动。有限体积法的基本思路易于理解，并能得出直接的物理解释。离散方程的物理意义就是因变量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小的控制体积中的守恒原理一样。有限体积法得出的离散方程对任意一组控制体积都得到满足因变量的积分守恒，因此也满足整个计算区域。正是因为有限体积法有这么多优点，所以在一维非饱和土壤的水分运动模型研究中尝试采用有限体积法。

1 非饱和土壤水分运动数学模型

田间水分运动是非常复杂的，受多种因素影响诸如降雨、畦灌、喷灌及蒸发蒸腾等。因此，有多种数学模型。这里只研究最简单的地表湿润条件下的土壤水分运动问题的数学模型，假定土壤均质，各向同性，根据达西定律和质量守恒方程，可得到如下数学模型[4]：

=

[D（θ）

]-

θ=θa t=0，z=0

θ=θb t=0，z→∞或z=L

式中θ为含水率（cm3/cm3），θa为均匀分布的初始含水率（cm3/cm3），θb为地表因湿润条件而维持不变的含水率（cm3/cm3）， D（θ）为扩散率（cm2/min）， K（θ）为导水率（cm/min）。z为位置坐标，地表取为原点，向下为正方向。

2 数值模型模拟

2.1 建立数值模型

将全部求解区域z=0-L离散化为n个单元，共为n个节点，但上边界点和下边界点并没有计入。距离步长为Δz，时间步长为Δt。

对2到n-1的节点，按有限体积法的完全隐格式可把原方程写出如下形式：

（θik+1-θik）=

-

-

K（

θ）-K（

θ）

上述形式可化简为如下形式：

Aθ+Bθ+Cθ=b

式中，

A=-

B=++

Ci=-

bi=θ-[K（θ）-K（θ）]

i=2，3，4，……， n

对于上边界第一个节点，原方程可以写成如下形式：

（θ-θ）=

-

-

[K（θ）-K（θ）]

可化简为Bθ+Cθ=b这样的形式

对于下边界第n个节点，原方程可以写成如下形式：

（θ-θ）=

-

-

[K（θ）-K（θ）]

同样也可以可化简为Aθ+Bθ=b

从而第一个节点到第n个节点都可把原方程转化成a+aθ+aθ=b这样的形式：

B1 C1 0 0 0

A2 B2 C2 0 0

0 ? ? ?

0 0 An-1 Bn-1 Cn-1

0 0 0 An Bn ·θ

θ

[…]

θ

+

θ

+=b1

b2

[…]

bn-1

bn-2

即三对角矩阵[A][θ]k+1=[b]，此时只需利用追赶法便可进行求解。

其中，对于扩散项中D（θ）这项采用算数平均法，即：

D（θ）=

D（θ）=

对于对流项中，K（θ）这项采用迎风格式，即：

K（θ）=K（θ）

K（θ）=K（θ）

2.2模型的实现

上述模型的运算是基于Matlab数学软件。程序采用有限体积法进行计算，经过方程的离散后，方程便化写成一个三对角矩阵。对于三对角矩阵采用的是TDMA法（追赶法）进行求解。TAMA算法是一种直接解法，是由高斯消元法应用于这种特殊方程而得到的一种解法。对于数学模型中的土壤扩散率D（θ）和土壤导水率K（θ）均是由试验测得。在程序中一般是定步长，根据需要改变步长。

2.3算例验证

对于一个程序的有效性是需要算例验证的。该程序的验证算例为地表湿润条件下的土壤水分运动问题[5]。土壤的垂向深度是100 cm，地表含水率θa=0.38 cm3/cm3，深度在100 cm时含水率θb=0.28 cm3/cm3，初始含水率θ0=0.28 cm3/cm3，饱和含水率θs=0.4 cm3/cm3，则可以表示成如下形式：

=

[D（θ）

]-

θa=0.38 t=0，z≥0

θb=0.28 t=0，z=100 cm

K（θ）=10.25（θ/θs）3.6

D（θ）=0.0379（θ/θs）3.487

式中，D（θ）为扩散率（cm2/min），K（θ）为导水率（cm/min）。设置时间步长为1 min，空间步长1 cm，用Matlab编程序进行数值模拟，其计算结果如图1。

从图1可以看出，土壤含水率随着时间的推移，逐步向下入渗。当t=5 min时入渗到20 cm处，在t=25 min时入渗到35 cm左右，最后到t=55 min时，入渗到了65 cm左右处。从整个模拟的趋势上看是符合客观规律的。从而也能说明，有限体积法也可以对一维非饱和土壤水分运动进行数值模拟。

3 有限体积法和有限差分法比较

在实验室模拟土壤水分运动，只考虑地表湿润条件，不考虑其他情况。选取的土壤为中壤土，初始含水率为0.09 cm3/cm3，土样容重为1.4 g/cm3。为保证土柱初始含水率均匀和密度均一，先对土样进行了风干和过筛，然后土样按控制的容重装入直径为10 cm，高为80 cm的有机玻璃管，土柱长60 cm。在土柱的顶部加水，形成一个薄水层（1 cm），由自动加水器供水，使水位保持不变，并计时。然后用γ射线测得1 230 min时刻的土柱剖面含水率。

数学模型如下：

=

[D（θ）

]-

θa=0.42 t=0，z≥0 （1）

θb=0.09 t=0，z=60 cm （2）

D（θ）=0.0002e19.902θ cm2/min

K（θ）=555.22θ13.181 cm/min

（1）式为上边界条件，即在表层含水率0.42 cm3/cm3。（2）式为下边界条件，垂直深度60 cm时含水率为0.09 cm3/cm3。初始含水率为0.09 cm3/cm3。D（θ）扩散率（cm2/min）是通过水平土柱吸渗法测定得出。K（θ）导水率（cm/min）由K（θ）=C·D（θ）导出，其中C为比水容量，通过土壤水分特征曲线得出。

当时间步长为30 min，空间步长为2 cm，网格数为30时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图2。图2中，可明显看出有限体积法模拟效果优于有限差分法，此时的网格数为30。

当时间步长不变仍为30 min，空间步长为1cm，网格数为60时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图3。图3中，有限体积法和有限差分法两者的模拟的效果基本相同，并没有太大差别。此时的网格数为60，网格加密了一倍。与图2比较，有限体积法并没太大的改变，而有限差分法模拟效果明显变好。

当时间步长不变仍为30 min，空间步长为2/3cm，网格数为90时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图4。图4中，有限体积法和有限差分法模拟结果基本重合在一起。此时的网格数为90，较之图2网格加密了两倍，有限差分法模拟效果明显优于图2。

4 结论

经过算例的验证，有限体积法同样适用于一维非饱和土壤水分运动。与有限差分法比对发现，两种方法的数值模拟结果相近。但是在网格稀疏的时候，有限体积法模拟的效果更好。在数值计算中，网格的疏密会影响计算精度和计算时间。网格越密计算精度越高，但是运算时间也会随之增长。通过改变网格的疏密后有限体积法的计算值与实测值并无太大误差，表明可选择有限体积法在网格相对较疏时进行数值运算，效果较好。

参考文献：

[1] 李人宪.有限体积法基础[M].北京：国防工业出版社，2008.

[2] 吕岁菊，李春光.土壤水盐运移规律数值模拟研究综述[J].农业科学研究，2005，26（1）：80-84.

[3] VERSTEEG H K， MALALASEKERA W. An Introduction to Computational Fluid Dynamics： The Finite V Method[M]. Edinburgh， England： Pearson Education Limited，1995.

[4] 雷志栋，杨师秀，谢森传.土壤水动力学[M].北京：清华大学出版社，1988.

[5] 赵慧丽，张 弥.降雨入渗对基坑非饱和土瞬态含水率影响的数值计算[J].石家庄铁道学院学报，2000，12（13）：43-46.

当时间步长为30 min，空间步长为2 cm，网格数为30时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图2。图2中，可明显看出有限体积法模拟效果优于有限差分法，此时的网格数为30。

当时间步长不变仍为30 min，空间步长为1cm，网格数为60时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图3。图3中，有限体积法和有限差分法两者的模拟的效果基本相同，并没有太大差别。此时的网格数为60，网格加密了一倍。与图2比较，有限体积法并没太大的改变，而有限差分法模拟效果明显变好。

当时间步长不变仍为30 min，空间步长为2/3cm，网格数为90时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图4。图4中，有限体积法和有限差分法模拟结果基本重合在一起。此时的网格数为90，较之图2网格加密了两倍，有限差分法模拟效果明显优于图2。

4 结论

经过算例的验证，有限体积法同样适用于一维非饱和土壤水分运动。与有限差分法比对发现，两种方法的数值模拟结果相近。但是在网格稀疏的时候，有限体积法模拟的效果更好。在数值计算中，网格的疏密会影响计算精度和计算时间。网格越密计算精度越高，但是运算时间也会随之增长。通过改变网格的疏密后有限体积法的计算值与实测值并无太大误差，表明可选择有限体积法在网格相对较疏时进行数值运算，效果较好。

参考文献：

[1] 李人宪.有限体积法基础[M].北京：国防工业出版社，2008.

[2] 吕岁菊，李春光.土壤水盐运移规律数值模拟研究综述[J].农业科学研究，2005，26（1）：80-84.

[3] VERSTEEG H K， MALALASEKERA W. An Introduction to Computational Fluid Dynamics： The Finite V Method[M]. Edinburgh， England： Pearson Education Limited，1995.

[4] 雷志栋，杨师秀，谢森传.土壤水动力学[M].北京：清华大学出版社，1988.

[5] 赵慧丽，张 弥.降雨入渗对基坑非饱和土瞬态含水率影响的数值计算[J].石家庄铁道学院学报，2000，12（13）：43-46.

当时间步长为30 min，空间步长为2 cm，网格数为30时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图2。图2中，可明显看出有限体积法模拟效果优于有限差分法，此时的网格数为30。

当时间步长不变仍为30 min，空间步长为1cm，网格数为60时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图3。图3中，有限体积法和有限差分法两者的模拟的效果基本相同，并没有太大差别。此时的网格数为60，网格加密了一倍。与图2比较，有限体积法并没太大的改变，而有限差分法模拟效果明显变好。

当时间步长不变仍为30 min，空间步长为2/3cm，网格数为90时，计算1 230 min后的入渗情况。有限体积法和有限差分法的计算结果如图4。图4中，有限体积法和有限差分法模拟结果基本重合在一起。此时的网格数为90，较之图2网格加密了两倍，有限差分法模拟效果明显优于图2。

4 结论

经过算例的验证，有限体积法同样适用于一维非饱和土壤水分运动。与有限差分法比对发现，两种方法的数值模拟结果相近。但是在网格稀疏的时候，有限体积法模拟的效果更好。在数值计算中，网格的疏密会影响计算精度和计算时间。网格越密计算精度越高，但是运算时间也会随之增长。通过改变网格的疏密后有限体积法的计算值与实测值并无太大误差，表明可选择有限体积法在网格相对较疏时进行数值运算，效果较好。

参考文献：

[1] 李人宪.有限体积法基础[M].北京：国防工业出版社，2008.

[2] 吕岁菊，李春光.土壤水盐运移规律数值模拟研究综述[J].农业科学研究，2005，26（1）：80-84.

[3] VERSTEEG H K， MALALASEKERA W. An Introduction to Computational Fluid Dynamics： The Finite V Method[M]. Edinburgh， England： Pearson Education Limited，1995.

[4] 雷志栋，杨师秀，谢森传.土壤水动力学[M].北京：清华大学出版社，1988.

[5] 赵慧丽，张 弥.降雨入渗对基坑非饱和土瞬态含水率影响的数值计算[J].石家庄铁道学院学报，2000，12（13）：43-46.