第五章--降雨和灌水入渗条件下土壤水分运动2Word格式.doc

第五章--降雨和灌水入渗条件下土壤水分运动2Word格式.doc

《第五章--降雨和灌水入渗条件下土壤水分运动2Word格式.doc》由会员分享，可在线阅读，更多相关《第五章--降雨和灌水入渗条件下土壤水分运动2Word格式.doc（23页珍藏版）》请在冰点文库上搜索。

R（t）——降雨或灌水入渗强度。

（3）当降雨或灌水强度大于土壤入渗强度，地表形成积水，成为压力入渗。

即（一类边界）

（2-5-4）

H（t）——地表积水深度。

当地表积水而没有产生径流时，地表水深为H（t）；

若产生地表径流，积水深度H（t）可根据来水强度R（t）、土壤入渗强度i（t）及地表径流量Q（t）求得。

2．下边界条件

（1）地下水埋深较小，以地下水位作边界。

当地下水位变化很小或基本保持不变时，则地下水面处土壤含水率为饱和含水率（地下水面离地面距离为d），故

（2-5-5）

当地下水面随时间而变化时，即地下水埋深d为时间t函数d（t），则地下水面处负压为零，即

（2-5-6）

（2）地下水埋深较大的情况，在计算范围内，下边界土壤剖面含水率保持初始含水率，即

（2-5-7）

在上述条件下，如初始含水率上下一致，，得则

（2-5-8）

k（θi）––––离地表距离d处断面通量。

（3）不透水边界。

下边界为流量等于零的边界，即

（2-5-9）

上述表明，研究入渗时边界条件是较为复杂的，所以，计算方法也较为复杂。

第二节土壤水入渗线性化方程的近似解

在垂直入渗情况下，一维土壤水分运动的基本方程可写作：

（2-5-10）

如降雨或灌水前的初始含水率（在土壤剖面上含水率均匀分布）为θi，则初始条件为

（2-5-11）

在地表有一薄水层时，表层含水率等于饱和含水率θS；

在地下水埋深较大时，计算时段内入渗水不会到达下边界。

为此，下边界土壤含水率不变，等于初始含水率，则边界条件可以写作以下形式：

（2-5-12）

由于式（2-5-10）为非线性方程（因为扩散度D（θ）及水力传导度k（θ）均为待求含水率θ的函数），求解比较困难，为了简化计算，近似地以平均扩散度代替D（θ），并以代替，则式中（2-5-10）可简化为

（2-5-13）

式中（2-5-13）为常系数线性方程，可以用拉普拉斯变换求解。

对式（2－5－13）采用拉普拉斯变换后可得象函数方程：

（2-5-14）

式（2-5-14）的通解为

（2-2-15）

式（2-5-12）经拉氏变换后，得：

（2-2-16）

（2-2-17）

根据边界条件式（2－5－16）、式（2一5－17）确定常数：

代入式（2-5-15），得象函数的表达式为

（2-5-18）

进行逆变换后，得含水率的表达式为

（2-5-19）

补余误差函数可自表1－2－2查得。

式（2－5－19）中可用下式计算：

（2-5-20）

若已知D与θ的关系式，代入式（2－5－20）积分，即可求得。

采用式（2-5-19）求得的土壤剖面上含水率分布如示意图2-5-2所示。

由于地表的入渗强度，为了推求入渗强度，首先根据的象函的表达式求：

（2-5-51）

地表处，z=0，则

（2-5-21’）

在入渗初期，，相当于时，，，则式（2-5-21）可近似写成：

（2-5-21’’）

经逆变换得：

（2-5-22）

入渗初期：

［当D（θ）取平均值D时］

（2-5-23）

入渗时间较久，即当之，相当干时

o

t（min）

i（cm/min）

Ks

图2-5-3入渗率随时间变化图

，代入（2-5-21’）式，则，所以

则入渗时间久时，入渗强度（i→ks）为

（2-2-24）

自式（2－5－23）、式（2－5－24）得入渗速度在时间上的变化过程如图（2-5-3）所示。

第三节Green－Ampt模型的入渗解

Green－Ampt模型［50］是1911年提出的一种简化的入渗模型，它是建立在毛管理论基础上的一种入渗模型。

假定土壤是由一束直径不相同的毛管组成，水在土壤入渗过程中，湿润锋面几乎是水平锋面，且在锋面上各点的吸力水头均为Sm。

锋面后面的土壤含水率为均一的，如图（2-5-4）所示。

所以k（θ）也为常数，这种模型又称活塞模型。

根据达西定律：

（2-5-25）

H——地面以上水层厚度；

Sm——锋面处土壤负压；

z一锋面推进距离。

式（2-5-25）为单位时间，单位面积流入土体的水量。

根据水量平衡原理，应等于土体内增加的水量，即

（2-5-26）

式（2-5-26）积分：

所以

（2-5-27）

式（2－5－27）为z～t关系式，原则上可以求得任何时刻t时入渗锋面所达到的位置，当然也就不难求得该时刻的累计入渗量：

（2-5-28）

H→0时，式（2-5-27）可写作：

（2-2-27’）

或由式（2-5-26），当t很小时，该式的H+Sm+z项中z略去，所以。

此时

积分得（2-2-27’’）

t时入渗总量

（2-5-29）

式（2-5-27’’）表明，入渗初期，入渗深度z与成正比。

将I对t求导，得：

（2-5-30）

当t大时，式（2-5-26）中Z>

>

H+Sm，因此，则由式（2－5－25）可知：

（2-5-31）

即入渗强度近似等于土壤饱和渗透系数。

第四节水平入渗条件下的Philip解法

水平入渗条件下的Philip解[51]是一种半解析法，即前半部用解析法，利用博茨曼（Boltzmann）变换，将偏微分方程转换为常微分方程；

后半部采用迭代计算，求解常微分方程。

由于求解过程中未作过分简化，求得结果较为严密。

一、水平入渗的常微分方程推求

水平入渗的基本方程

（2-5-32）

将式（2-5-32）中基本方程改写为以坐标x（θ，t）为变量的方程，根据第二章中方程（2-2-17），在水平入渗时应为

（2-5-33）

采用Boltzmann变换，引入变量从λ（θ），且令

（2-5-34）

则（2-5-35）

故（2-5-36）

（2-5-37）

将式（2－5-36）、式（2－5－37）代入式（2－5－33）得：

经整理后得微分方程

（2-5-38）

由边界条件已知：

为了求得λ～θ的关系式，将式（2－5－38）常微分方程自θi至θ积分得：

（2-5-39）

式（2-5-39）为λ～θ关系式，若已知D（θ）关系式代入上式即可求得λ～θ关系式，因，即可由λ～θ求得θ（x，t），从而求得剖面上任何时间，任何距离的含水率分布。

若能实测λ（θ），则代入上式可求得D（θ）关系曲线。

二、迭代计算法

在D（θ）已知时，式（2－5－39）为λ～θ关系式，它的关系曲线如图2－5－5所示，Philip的送代法是将θ0→θi等份成n份，步长为Δθ=（θ0-θi）/n，若等分点按顺序其含水率为θ1，θ2，…，θn-1，θn。

相应各等分点相应λ值为λ1，λ2，…λn，两个相邻等分点中点的含水率为θ1/2，θ1+1/2，…，θ（n-2）+1/2，其相应扩散率值为D1/2，D1+1/2，…D（n-2）+1/2，任一点含水率θr，可用下式计算

（2-5-40）

Dr+1/2取平均值，定义

（2-5-41）

写出θ1/2，θ1+1/2，…θ（n-2）+1/2各点的式（3-5-39）：

（2-5-42）

令（2-5-43）

当r=0时，（2-5-44）

式（2－5－43）表示λ～θ关系曲线所包含的面积。

式（2-5-44）表示θ=θr+1/2时，θr+1/2～θn之间λ～θ曲线所包含面积的近似值。

由图可知面积为

（2-5-45）

δr+1—––为λ～θ的曲线与λr，λr+1/2之间所包含面积减去矩形面积后的一小块面积。

由式（2－5一43）及式（2－5－45）得：

当Δθ取得足够小时，则δr+1可以忽略，所以

（2-5-45’）

则式（2一5一42）可以改写为

（2-5-46）

由图形可知：

Δr––––矩形面积与λ～θ曲线与λr－1/2、λr之间包含的面积之差。

上式可写为

（2-5-47）

当Δθ足够小时，同样，（Δr-δr）可忽略。

（2-5-48）

根据式（2－5－47），我们可以得到：

（2-5-49）

根据式（2－5－46）和式（2－5－49），若已知J1/2及D（θ）则可交替使用此二式，求得λ1，λ2…λr…，λn-1。

因此，为了求得λ～θ关系值必须首先求得J1/2及D～θ关系曲线。

由上述可知J1/2代表了λ～θ曲线所包含面积的近似值，由式（2－5－49）可知：

若将面积看作一系列的长为（θ一θn），高为dλ的水平方向的矩形，则

（2-5-50）

λ––––θ和D（θ）的函数。

对于初始值J1/2，我们可将D（θ）考虑为在整个含水率变化范围内的平均值，由式（2－5－50）可知：

若用加权平均值来代替上式中D值，以D’表示：

（2-5-51）

——的加权因子。

其中等号右端系数2是为使等式衡等，D’=D’（若将D’代替D（θ）即可算得该系数）。

将式（2-5-51）写成有限差的形式：

（2-5-52）

因D’为常数，直接利用入渗问题线性化解，以代人得：

（2-5-53）

将式（2-5-53）代入式（2-5-50）得：

（2-5-54）

式（2-5-54）为的初始近似值第一次迭代值，因此可写作。

将α=λ/2（D’）1/2代人上式，且令积分上限为β0=β→∞以n·

Δθ代替（θ0－θn）则式（2-5-54）可写作：

（2-5-55）

为了求解式（2－5－55）积分式，利用erfcα=1－erfα，则

已知

得（2-5-56）

当时,，所以式（2－5－56）为

（2-5-57）

将式（2-5-57）代入式（2-5-55）得：

（2-5-58）

有了J1/2就可以由式（2一5－56）和式（2－5－49）重复交替计算求得λ～θ关系曲线，但由于J1/2是估算值，存在偏差，求得的λ～θ曲线必然也有偏差，为了使迭代计算尽快收敛，采用两种不同方法计算Jn－1/2，不断地修正J1/2，以使两种方法求得的Jn－1/2值足够的接近，以便得到一个比较精确的λ～θ曲线。

第一种方法采用

（2－5－59）［即式（2－5－49）］

第二种方法采用

（2－5－59）［即式（2－5－49）

JS式中，由于θ变化于θn，θn-1之间，若Δθ取得足够小，则D可以看作为一常量Dn－1/2：

（2-5-61）

式（2-5-61）表明，在θn和θn－1范围内将D视作常量，当Δθ取得足够小时，不会引起明显误差。

上述问题可为求解半无限长水平土柱的入渗问题，即

（2-5-62）

根据土壤水入渗的线性化方程的近似解法，式（2－5－62）定解问题的解为

（2-5-63）

式中A（y）定义为

（2-5-64）

其中（2-5-65）

A（y）的函数表见表（2-5-1）。

对y较大时，A（y）可以渐近级数表示：

（2-5-66）

λn-1由式（2－5－59）求得，对某一个λ～θ曲线λn-1及Dn-1/2为常量，所以，A（y）也为常数，可采用式（2－5－64）计算，根据y求A（y）。

首先，已知

（2-5-67）

式中的β为虚拟变量，可以改写为

（2-5-68）

由式（2-5-68）求得erfc（y），代入式（2-5-64）即可求得A（y）、A（y）～y已列于表（2-5-1）中，所以根据某Dn－1/2，λn－1求得y，查表2－5－1求得。

将式（2－5－63）代人式（2－5－60）得：

（2-5-69）

当Δθ足够小时，且假定J1/2是正确的，则JFn-1/2与JSn-1/2应该是相等的，但由于计算中采用了近似值两者误差为；

，若在所要求精度范围内即可，否则进行第二、第三、……，更多次的选代。

由与（J1/2）1；

可以得到一个更接近于实际的初始值（J1/2）2。

重复上述过程进行第二个循环的计算。

（J1/2）2一般可写作：

迭代P次后：

（2-5-71）

N––––式（2－5－70）两边相等（或误差在所要求的精确度范围内）时送代次数。

为了加速收敛，J1/2还有另一改进公式：

（2-5-71）

由于在其他著作中采用了与上述不同符号，这里介绍一下用其他符号表示的Philip的迭代公式。

Philip定义Ir+1/2为

（2-5-72）

所以，其他以I为参数的各式中均比J为参数的多除一个Δθ。

即（2-5-73）

迭代计算步骤：

（1）确定θ0至θn的分段数n，步长Δθ=（θ0-θn）。

（2）以计算平均扩散度。

（3）根据公式J1/2=2n·

Δθ（D’/π）1/2计计算第一次送代初值（J1/2）1，并由公式（2－5－46）计算λ1，由（J1/2）1和λ1即可计算（J1+1/2），再代入公式（2-5-49）求λ2，依次类推，最后求得。

（4）由以上计算同时可求得λn-1和Dn-1/2，代入公式（2-5-60）计算。

（5）计算误差Δ1=－，若Δ1满足要求，以上计算结束，否则要进行第二、三、…次的送代，直至满足要求为止。

（6）进行多次送代时，重复以上

（2）～（5）步骤，直至迭代误差Δ很小，且λ1，λ2，…、不再随J1/2改变时为止。

若计算结果不理想，可能是选用Δθ不合适，则可重新选用Δθ值，重复以上步骤进行迭代计算。

（7）计算出λ～θ关系曲线后，因，故，对于给定的时间t，由λ可计算x值，即得θ～x关系曲线。

也即求得了θ在剖面上的分布。

Philip水平入渗解的理论可用于水平土柱法测定土壤水扩散度。

土壤水分在较长的（理论上说是半无限边界）均质土柱中发生水平运动时，一维土壤水平运动方程如式（2－5－32）所示，经Boltzmann变换可得[如前述方程（2－5－39）]：

（2-5-39）

λ（θ）—––x，t的函数。

当水平土柱首端供水，根据观测水分运移资料，即可绘出λ～θ曲线，（），为由初始含水率θi至θx之间λ～θ曲线所包含的面积，如图2-5-5所示。

为在λ～θ曲线上对应含水率θx点的斜率。

因此，通过水平土柱试验即可求得土壤水扩散度D值。

第五节垂直入渗条件下的Philip解法

一维垂直入渗基本方程写成z（θ，t）为函数时为

（2-5-74）

对一类边界定解条件为

（2-5-75）

Philip级数解的形式[52]

（2-5-76）

相应边界条件为

（2-5-77）

（2-5-78）

若求得式（2-5-77）中系数ηi（θ）值，则可求得在一定θ时的位置z。

以下我们以待定系数法求解ηi（θ）。

式（2－5－74）可改写为

（2-5-79）

将z的解式（2-5-76）分别代入以上各项：

（2-5-80）

将式（2-5-80）代入式（2-5-79）归并后得：

（2-5-81）

式中前四项系数η1，η2，η3，η4分别为

（2-5-82）

由于式（2－5－8l）左端等于零，所以η1=η2=η3=η4…=0

η1=0，即

全式除以得

即

式（2－5－83）表明，η1（θ）即为λ（θ）为水平入渗时解。

同理，以η2=0求得：

（2-5-84）

η3=0得：

（2-5-85）

η4=0得：

（2-5-86）

根据以上各式，若已知η1（θ）则可逐步求得外η2（θ），η3（θ），…，即可求得方程（2－5－74）的解。

由于收敛较快，求得前四项就足够精确了。

上述η1（θ），η2（θ），η3（θ），η4（θ）都较复杂，η1（θ）即为水平入渗时的λ，为了求解这些值并不容易，为此Philip采用了递推公式求解。

以下以求得η2为例说明求解方法。

为了求解系数，我们可以写出与式（2-5-39）类似的近似系数一般方程，则

（2-5-87）

α，β—均为θ的已知函数。

式（2－5－87）左边可以写作式（2－5－88），如图2－5－6所示。

（2－5－88）

式（2一5－88）中右端第二个积分为

由图2－5－6可知，该积分近似等于，其中h为四边形BCDE的平均高度。

Δθ取得足够小时，DE可以看作为一直线，在梯形DFAC中h可以下式表示：

（2-5-89）

（2-5-90）

式（2-5-87）中，若θ=θr+1/2，则方程为

（2-5-91）

αr+1/2为θr+1<

θ<

θr范围内平均值，且以来近似（df/dθ）r+1/2（其负号表示θ值增加的方向与f值增加方向相反），代入式（2－5－91）得：

（2-5-92）

将式（2－5－89）、式（2－5－91）代入式（2－5－88），经整理后得

令（2-5