第07期（总第446期）吉林水利
2019年07
月
［文章编号］1009-2846(2019)07-0022-05
［收稿日期］
2018-01-29
［作者简介］杨旭（1987-），男，辽宁铁岭人，本科，助理工程师，主要从事水利水电工程等工作。
我国是农业用水大国，而沟灌作为我国目前应用于耕作物的主要地表灌溉方式之一，灌水条件复杂[1]
。灌水过程数学模型都是以非线性偏微分
方程表述，地表水流运动及土壤水分运动等条件复杂，很难求出精确解析解，采用模型数值方法求解其数值解是目前最有效的办法。而Matlab建模和仿真功能强大且直观易用，数值计算能力突出，在工程科学研究与计算方面得到广泛运用。对土壤水分运动和地表水流运动进行有效的耦合是研究复杂条件下灌水模型的主要手段。Gandolfi等[2]运用MacCormark数值方法对地表径流渗流耦合模型进行求解并验证了模型的精度；董勤各等[3]建立了一维畦灌的地表水流－土壤水耦合动力学模型并采用混合数值方法进行求解；章少辉等[4]在描述沟灌地表水流运动的基础上运用有限体积法对沟灌耦合模型进行数值模拟。虽然前人的研究取得了长足的进展，但仍需对此进行进一步研究。本文以辽宁某灌区玉米覆盖垄为研究对象，对于影响因素多、条件复杂的沟灌，研究沟灌过程中地表水流运动模型和土壤水分运动方程，将地表水流运动和土壤水分运动模型进行耦合，借助Matlab软件强大的数值计算功能进行数值模拟，研究成果可为类似工程提供参考。
1工程概况
在作物生长过程中，灌溉方式会影响其生长及产量。当前我国农业沟灌模式仍在北方地区大面积应用，研究其地表-土壤水动力学耦合关系具
有重要意义。辽宁某区玉米覆盖垄采用沟灌灌水方式进行灌溉，灌区土质为粉土、黏土，现场调研可知，单沟为梯形槽状明渠，沟深15cm左右，沟顶宽约
50cm，沟底宽约20cm，沟底纵坡0.4%。
现场实测土壤水力参数如表1所示，地表水流运动的基本参数如表2所示。对现场试验大田
进行灌水试验，得到的实测数据与模拟结果进行对比。
地表-土壤水动力学耦合模型在复杂沟灌条件下的应用
杨旭（葫芦岛平山供水有限责任公司，辽宁葫芦岛125000）
［摘要］反映地表水流运动和土壤水分运动的数学模型大多都是非线性偏微分方程，通常采用数值计算方法来求解较为复杂的地表水流运动模型和土壤水分运动方程。本文以辽宁某灌区玉米覆盖垄为例，利用零惯量模型和土壤水分运动方程（Richards方程）来描述沟灌灌水时水分的运动特征，构建沟灌地表-土壤水动力学的耦
合模型，通过迭代法实现有效耦合，借助Matlab软件进行数值模拟与计算，将模拟结果与实测结果相对比，模拟结果符合实际情况，能够满足实际工作要求。［关键词］沟灌；零惯量模型；Richards方程；地表-土壤水动力学的耦合模型
［中图分类号］S275.3［文献标识码］
B
土壤质地滞留含水率θr饱和含水率θs经验参数α1经验参数nn饱和导水率
K
S
粉土0.004220.090.0040.00220.49
土壤水力参数表1
22--
DOI:10.15920/j.cnki.22-1179/tv.2019.07.006图1沟灌模型示意图及其几何特征
地表-土壤水动力学耦合模型在复杂沟灌条件下的应用吉林水利杨旭2019年07月2沟灌地表水流运动模型
运动波模型、零惯量模型、水量平衡模型和完全水动力学模型是当前用于描述地表水流运动过程的4种主要模型。其中零惯量模型适用于一般
情况下畦灌和沟灌的地表水水流运动过程[5]
，求解
精度足够且求解过程简便，故本文采用零惯量模型对地表水流运动进行研究。沟灌地表水流运动属明渠非稳定流，满足圣维南方程[6]
。沟灌模型几何特征示意图如图1所示。在灌水过程工程实践中，沟灌地表水深h0与水流流速v都较小，此时圣维南运动方程中的加速度项和惯性项可忽略不计，即简化为零惯性模型[7]，其基本方程为：连续方程：鄣A鄣t+鄣Q鄣y+鄣Z鄣t=0（1）运动方程：鄣h0鄣y=S0-Sf（2）其中，A为地表水流的过水断面面积，m2；Q为地表水流的过水断面流量，m3/s；Z为累积入渗量，m2；y为水流推进距离，m；t为时间，s；h0为地表水深，m，h0=δ1Aδ2，δ1、δ2为拟合参数，对于沟灌时可取δ1=0.782，δ2=0.536[7]；S0为田面坡度；Sf为地表水流阻力坡降，Sf=Q2n2A2R43=Q2n2ρ1Aρ2，ρ1和ρ2为经验系数，沟灌可取ρ1=1、ρ2=3.3333，n为田面曼宁糙率。初始条件及边界条件分别为Q（y，t）=0，t=0h0（y，t）=0，t=鄣0
（3）
Q（y，t）=Q0y=0，0≤t≤t1
Q（y，t）=0y=0，t1≤t≤t2
A（y，t）=0y=yr，t2≤t≤t4
≤≤≤≤≤≤
≤≤≤≤≤≤
≤
（4）
A（y，t）=0y=ya，0≤t≤t3
Q（y，t）=0y=L，t3≤t≤t4
≤（5）
其中，Q0为沟首进水流量，m3/s；t1为沟首的断水时间，s；t2为沟首的垂直消退时间，s；t3为水流前
锋推进到沟尾端时间，s；t4为水流尾锋退水到沟尾
端时间，s；ya为水流前锋位置，m；yr为水流退水尾锋位置，m；L为沟长，m。
3沟灌土壤水分运动Richards方程
Richards方程常用于描述土壤水分运动[8]。为
简化计算，将土壤设为各向同性均匀多孔介质，忽略空气阻力、温度和蒸发量对入渗的影响。根据非饱和达西定律和质量守恒定律，可以建立二维Richards方程以描述二维非饱和土壤水分运动方
程[9]
。以负压水头h为独立变量的Richards方程基
本形式为：
C（h）鄣h鄣t=鄣鄣xK（h）鄣h鄣x≤≤+鄣鄣zK（h）鄣h鄣z≤≤-鄣K（h）
鄣z（6）
其中，h为负压水头，cm；C（h）为比水容量，
cm-1；K（h）为非饱和土壤的导水率，cm·min-1；x为
空间水平坐标，cm；z为空间垂直坐标，规定向下为正，cm
。
沟长（m）入沟流量（m3/s
）田面糙率地面坡度入渗系数入渗指数
1200.004220.090.0040.00220.49
地表水流运动基本参数表2
xyAFO
HDB
ah0
CB
z
23--图3沟灌迭代耦合过程流程图图2求解区域初始条件为：
h（x，z，t）＝hin（x，z）0≤x≤X，0≤z≤Z，t=0（7）其中，X，Z为求解区域水平和垂直方向的最大
距离，cm；hin（x，z）为初始负压水头分布，cm。
求解区域Ω如图2所示，则Richards方程边界条件为：h=h0-H，0＜t＜T，DE边
h=（h0-H）cosα，0＜t＜T，EG边
鄣h鄣x
=0，t＞0，AB和CD边
-K（h）鄣h鄣x-K（h）鄣h鄣z---1=0，t＞0，GF边
-K（h）鄣h鄣z+K（h）=0，t＞0，FA边
h=hin，t＞0，BC
----------------
----------------
-
边
（8）
其中，t为时间，min；T为停水时间，min；h0为
地表水深，cm；H为沟深，cm；α为EF边与x轴负方向的夹角。
4沟灌地表-土壤水动力学耦合模型及数
值模拟
在沟灌灌水过程中，地表水流沿田块沟槽方向向前流动，同时水分通过地表表面发生入渗。地表水流运动可利用零惯量模型来描述其运动过程，土壤水分运动可通过Richards方程来描述其
运动过程。将两模型进行迭代耦合求解，借助Matlab软件进行数值计算，实现对沟灌地表-土壤
水动力学耦合模型的数值模拟。两种模型耦合流程图如图3所示。
通过上述耦合过程对沟灌地表-土壤水运动
过程进行耦合之后，即可得到地表水流运动和土壤水分运动的过程。对零惯量模型，取空间步长Δy=1，水流运动基本参数如表2所示。对Richards方程取Ω=［0，100］×［-135，0］，α=45°，沟深
H=
OD=15，AB=135，AF=75，DE=10，BC=100，Δt=1，Δz=-5，Δx1=Δx2=Δx3=5。
耦合过程中累积入渗量Z采用Kostiakov入渗
公式[10]
即式（7）计算。
Z=krtαr（9）
其中，kr为入渗系数，t为入渗时间，αr为入渗指数。对实验大田进行沟灌试验，沟灌时长达2h，断
水阶段不计，并测出水流推进到不同沟长位置的
时间实测值，利用自制取芯杆于T=3000s时迅速
沿y=0和y=60处土壤剖面中心线插入（杆长
1.5m，截面直径3cm，插入深度1.2m，上端突出水
面），取出土芯，利用烘干法测出不同埋深处含水
地表-土壤水动力学耦合模型在复杂沟灌条件下的应用吉林水利杨旭2019年07月AFO
HGED
BC
z
Ω
xh0
利用零惯量模型求解到某结点处的初始地表水深h0
将初始地表水深h0作为Richards方程的上边界条件代入，对Richards方程差分求解得到该结点处初始负压水头h的二维分布情况
根据
θ（h）~h关系式求得土壤剖面各结点处的含水率θ（h）
对含水率积分可求得土壤水分累积入渗量Z将累积入渗量Z代入零惯量模型再一次求解得到该结点处更新后的地表水深h0（1）
将更新后的地表水深h0（1）作为Richards方程的上边界条件，对Richards方程差分求解得到该结点处新的负压水头h的分布情况
根据θ（h）~h关系式求得土壤剖面各结点处含水率
对含水率进行积分得到累积入渗量Z（1）
否
是该结点处地表水深h0和入渗量Z
是否满足收敛要求
结束本结点的耦合计算进入下一结点的计算
24--