柱坐标系下地下滴灌土壤水分运动规律研究

柱坐标系下地下滴灌土壤水分运动规律研究

目前，中国的地下灌溉设计主要参照《微灌技术技术规范》。该系统根据土壤灌溉和简单孔口流量的原理进行规划设计。但地下滴灌的灌水器埋于地下,实际田间的滴头流量无法直接观测,其土壤水分运动规律也有别于地面滴灌。研究地下滴灌土壤水分运动规律对提高其灌溉理论水平,指导地下滴灌工程的规划设计有着重要的意义。数值模拟法由于具有无需建立复杂的专门设备即可在计算机上进行各种试验处理的优点,所以自20世纪60年代以来,它已成为研究土壤水分运动的主要方法。仵峰等(1996)将地下滴灌的土壤水分运动简化为垂直平面的二维流动,建立的h方程采用了交替隐式差分解法(AlternativeDirectionImplicit,简称ADI法),李光永等(1996)对建立的轴对称条件下的θ方程进行了ADI法求解,许迪等(2001)对柱坐标系下的h方程进行了Galerkin求解,模型都得到了较好的验证。事实上,地下滴灌土壤水运动属于典型的三维饱和-非饱和流动,θ方程并不适宜。另外,ADI法原理相对简单,它能弥补隐式或中心差分法以及Galerkin法计算工作量大的缺点,具有占用内存少、计算速度快、稳定性较好等优点。基于此,本文对建立的柱坐标系下地下滴灌土壤水分运动的h方程进行了ADI法求解,模型进行了室内试验的验证。1下u3000下灌地水分运动h方程若假定土壤均质、各向同性,地下滴灌的三维饱和-非饱和流动可作为轴对称的二维问题进行处理。建立柱坐标系下地下滴灌土壤水分运动的h方程如下(z轴向下为正):C(h)∂h∂t=∂∂r[Κ(h)∂h∂r]+Κ(h)r∂h∂r+∂∂z[Κ(h)∂h∂z]-∂Κ(h)∂z-Sr(z,t)(1)式中:C(h)和K(h)分别为土壤容水度和非饱和土壤水力传导度,均为土壤负压h的函数;r,z分别表示柱半径和垂直坐标;Sr为根系吸水率。2模型的确定和分解2.1计算区域在上述假定条件下,可以只研究如图1中阴影部分的半个领域。其中,线段de为滴头直径;Ze为滴头埋深;Z、R为研究区域纵横向长度。2.2饱和土壤侵蚀面为0.假定在灌水初始时刻,土壤负压水头均匀分布,计算区域内各点土壤负压都等于某一固定值h0,即h(r,z,0)=h0(r,z)∈除去弧段的整个区域(2)对于de弧边界,可以设灌水初期地下滴头与土壤接触面上为饱和状态,即h=0,当(z-Ze)2+r2=r20,r>0时(3)式中,r0为滴头半径。实际上,这个饱和弧区域大于滴头弧边界,但这对土壤水后期分布没有太大影响。2.3排水边界条件如图1,设滴头半径为r0,对称性边界oa(除去滴头部分)及bc均属于给定流量边界,可将其作为虚拟的不透水边界处理,即-Κ(h)∂(h-z)∂z=0,所以进一步可得:①对于oa边界(除去de段):∂h∂z=1,当r=0,0<z<Ze-r0或Ze+r0<z<Z时(4)②对于bc边界:∂h∂z=1,当r=R,0<z<Z时(5)③计算区域的下边界ab,考虑地下水埋深较大的情况,采用自由排水边界条件,则∂h∂z=0,当z=Z,0≤r≤R时(6)④计算区域的上边界oc,垂直方向的水流通量应与蒸发强度或降雨强度相等,即-Κ(h)∂(h-z)∂z=ε,当z=0,0≤r≤R时(7)同时还需满足ha≤h≤hs=0,ha和hs分别为表层土壤的最小和最大负压水头(一般不考虑表层积水情况);另外,式(7)中,当ε表示蒸腾强度时取负值;当ε表示降雨时取正值;当降雨和蒸腾都不考虑时ε取零。⑤对于当t>0时的de弧段,由于其周围的饱和区域会随灌溉时间的增加而逐渐增大,所以de弧段实际上是一个动边界条件。3时间步长及整体系统以Δr,Δz为距离步长将区域划分为矩形网格,时间步长为Δt,结点记为(i,j,n),其中,i=1,2,…,I+1;j=1,2,…,J+1;n=1,2,…,N。3.1方程的离散度j.2e-snri,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,kni,t32Cni,j*hn+1/2i,j-hni,jΔt=Κn+1/2i+1/2,j(*hn+1/2i+1,j-*hn+1/2i,j)-Κn+1/2i-1/2,j(*hn+1/2i,j-*hn+1/2i-1,j)(Δr)2+Κn+1/2i,jiΔr(*hn+1/2i+1,j-*hn+1/2i-1,j)2Δr+Κni,j+1/2(hni,j+1j-hni,j)-Κni,j-1/2(hni,j-hni,j-1)(Δz)2-Κni,j+1-Κni,j-12Δz-Snri,j(8)取r1=Δt2(Δr)2,r2=Δt2(Δz)2,r3=Δt2Δz,经整理得:(-2r1Kn+1/2i-1/2,j+r1Kn+1/2i,j/i)*hn+1/2i-1,j+(Cni,j+2r1Kn+1/2i+1/2,j+2r1Kn+1/2i-1/2,j)*hn+1/2i,j+(-2r1Kn+1/2i+1/2,j-r1Kn+1/2i,j/i)*hn+1/2i+1,j=2r2Kni,j+1/2hni,j+1+(Cni,j-2r2Kni,j+1/2-2r2Kni,j-1/2)hni,j+2r2Kni,j-1/2hni,j-1-r3(Kni,j+1-Kni,j-1)-Snri,jΔt若令ai,j=-2r1Kn+1/2i-1/2,j+r1Kn+1/2i,j/ibi,j=Cni,j+2r1Kn+1/2i+1/2,j+2r1Kn+1/2i-1/2,jci,j=-2r1Kn+1/2i+1/2,j-r1Kn+1/2i,j/ifi,j=2r2Kni,j+1/2hni,j+1+(Cni,j-2r2Kni,j+1/22r2Kni,j-1/2)hni,j+2r2Kni,j-1/2hni,j-1-r3(Kni,j+1Kni,j-1)-Snri,jΔt上述差分方程可简化为:ai,j*hn+1/2i-1,j+bi,j*hn+1/2i,j+ci,j*hn+1/2i+1,j=fi,j(9)若初始tn时刻水头已知,则通过求解式(9)即可求得带*号的负压水头——在tn+1/2时刻的过渡解。/2i,10102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,102i,Cn+1/2i,jhn+1i,j-*hn+1/2i,jΔt=Κn+1/2i+1/2,j(*hn+1/2i+1,j-hn+1/2i,j)-Κn+1/2i-1/2,j(hn+1/2i,j-hn+1/2i-1,j)(Δr)2+Κn+1/2i,jiΔr(*hn+1/2i+1,j-*hn+1/2i-1,j)2Δr+Κn+1i,j+1/2(hn+1i,j+1-hn+1i,j)-Κn+1i,j-1/2(hn+1i,j-hn+1i,j-1)(Δz)2-Κn+1/2i,j+1-Κn+1/2i,j-12Δz-Sn+1/2ri,j(10)经整理得:-2r2Kn+1i,j-1/2hn+1i,j-1+(2Cn+1/2i,j+2r2Kn+1i,j+1/2+2r2Kn+1i,j-1/2)hn+1i,j+2r2Kn+1i,j+1/2hn+1i,j+1=2r1Kn+1/2i+1/2,j+r1Kn+1/2i,j/i)*hn+1/2i+1,j+(2Cn+1/2i,j-2r1Kn+1/2i+1/2,j-2r1Kn+1/2i-1/2,j)*hn+1/2i,j+(2r1Kn+1/2i-1/2,j-r1Kn+1/2i,j/i)*hn+1/2i-1,jr3(Kn+1/2i,j+1-Kn+1/2i,j-1)-Sn+1/2ri,jΔt若令a′i,j=-2r2Kn+1i,j-1/2b′i,j=-2Cn+1/2i,j+2r2Kn+1i,j+1/2+2r2Kn+1i,j-1/2c′i,j=2r2Kn+1i,j+1/2f′i,j=(2r1Kn+1/2i+1/2,j+r1Kn+1/2i,j/i)*hn+1/2i+1,j+(2Cn+1/2i,j2r1Kn+1/2i+1/2,j-2r1Kn+1/2i-1/2,j)*hn+1/2i,j+(2r1Kn+1/2i-1/2,j-r1Kn+1/2i,j/i)*hn+1/2i-1,jr3(Kn+1/2i,j+1-Kn+1/2i,j-1)-Sn+1/2ri,jΔt上述差分方程可简化为:a′i,jhn+1i-1,j+b′i,jhn+1i,j+c′i,jhn+1i+1,j=f′i,j(11)将已求得的tn+1/2时刻的过渡解当作已知代入(11),即可求得tn+1时刻各结点的负压水头。3.2分散条件下的分辨率与等式2相同，初始条件的差分方程为h0i,j=h0(12)弧边界可近似各个点z,则jN=jM+5,其中jM,jN分别为初始饱和区上、下缘结点编号。则de弧边界可近似有10个结点(如图2),分别为:h02,jΜ=h03,jΜ=h03,jΜ+1=0h04,jΜ+1=h04,jΜ+2=h04,jΜ+3=h04,jΜ+4=0h03,jΜ+4=h03,jΝ=h02,jΝ=0(13)与等式4相同，即从方程de中取出的有限方程为hn2,j-hn1,j=Δz(14)与等式5相同，即边界市场的离散方程为hΙ+1,jn-hΙ,jn=Δz(15)对应于等式6，即从方程ab的角度来看的离散方程为hi,J+1n-hi,Jn=0(16)与等式7相同，即限制oc的有限方程为-Κi,1nhi,2n-hi,1nΔz=ε⇒hi,2n-hi,1n-εΚi,1nΔz(17)饱和区内控制边界条件为计算方便,以滴头边界作为初始饱和边界计算得第一时段末的负压水头。此时,若没有负压水头接近零(根据精度要求,可取h≈0.0001),在下一时段仍按初始时段边界条件进行计算;若发现已有负压水头接近零,则饱和区已经扩展到新的位置,下一时段应以新的饱和边界进行运算。而饱和区内部结点负压水头都已为零,可不再计算。如此每个时段进行数值计算,每个计算时段结束后均需按饱和区新达到的位置来确定后续时段的边界条件。3.3生长分区方程对于差分方程(9),假定nΔt为初始时刻,则由初始条件知相应的hi,jn应为已知,也即fi,j为已知。若在第j行沿水平r方向依次列出各结点差分方程,则可组成如下三对角方程组:(a2,jb2,jc2,ja3,jb3,jc3,j⋱⋱⋱aΙ-1,jbΙ-1,jcΙ-1,jaΙ,jbΙ,jcΙ,j)(Ι-1)×(Ι+1)×(*h1,jn+1/2*h2,jn+1/2⋮*hΙ,jn+1/2*hΙ+1,jn+1/2)(Ι+1)×1=(f2,jf3,j⋮fΙ-1,jfΙ,j)(Ι-1)×1(18)又由oa和bc上的边界条件离散式(14)和(15)补充两个方程:-*h1,jn+1/2+*h2,jn+1/2=Δz和-*hΙ,jn+1/2+*hΙ+1,jn+1/2=Δz于是,共有I+1个方程,联立以上3式可求解I+1个未知水头。采用多次迭代法将差分方程线性化,结点间土壤水力特性参数的取值采用几何平均值法,即对r方向上的半结点有:Ki±1/2,j=(Ki,jKi±1,j)1/2,Ci±1/2,j=(Ci,jCi±1,j)1/2,对于z方向上半结点处的处理可采用同样的方法。依次求解的各对角方程组,便可以求得在过渡时刻(n+1/2)Δt计算区域内所有结点的负压水头*hi,jn+1/2。同理,由差分方程(11),将求得的过渡解*hi,jn+1/2作为初始条件,则f′i,j相当于已知。在第i列沿垂直z方向依次列出各结点差分方程,可组成如下三对角方程组:(a′i,2b′i,2c′i,2a′i,3b′i,3c′i,3⋱⋱⋱a′i,J-1b′i,J-1c′i,J-1a′i,Jb′i,Jc′i,J)(Ι-1)×(Ι+1)×(hi,1n+1hi,2n+1⋮hi,Jn+1hi,J+1n+1)(Ι+1)×1=(f′i,2f′i,3⋮f′i,J-1f′i,J)(Ι-1)×1(19)再由oc和ab上的边界条件离散式(16)和(17)补充两个方程:-hi,Jn+1+hi,J+1n+1=0和hi,2n+1-hi,1n+1=-εKi,1n+1Δz同样,共有I+1个方程,可联立求解I+1个未知水头。由此,用迭代法可求出计算区域中所有结点在(n+1)Δt时刻的负压水头。反复交替逐时段地使用并求解上述两个三对角方程组,就可以求出计算区域中全部结点在各个时段末的负压水头,直到所规定的时刻为止。4数值模型的检验作者于2003年4月在武汉大学农田水利与水环境重点试验室开展了室内地埋点源条件下的地下滴灌灌水试验。4.1土壤湿性边界试验在自制的40cm×40cm×50cm(长×宽×高)有机玻璃箱中进行,装置如图3。试验采用粘壤土,滴头埋在地表面