降雨和灌水入渗条件下土壤水分运动2Word格式.docx
《降雨和灌水入渗条件下土壤水分运动2Word格式.docx》由会员分享，可在线阅读，更多相关《降雨和灌水入渗条件下土壤水分运动2Word格式.docx（31页珍藏版）》请在冰点文库上搜索。
（2-5-5）
$（d,t）=乞，z=d,t=0
h（d,t）=0,z=d,t>
当地下水面随时间而变化时，即地下水埋深d为时间t函数d（t）,则地下水面处负压为零，即
（2-5-6）
h（d（t）,t）=0,z=d（t）,t0
（2）地下水埋深较大的情况，在计算范围内，下边界土壤剖面含水率保持初始含水率，
即
在上述条件下，如初始含水率上下一致，耳（z）-耳,得「i（z）=0则
胡i
q=—D（旳【k（m）二k（R）z=d,tO（2-5-8）
cz
k（0i）--离地表距离d处断面通量。
（3）不透水边界。
下边界为流量等于零的边界，即
hh
q--k（h）
（1）=0,1,z=d,t，0（2-5-9）
dzdz
上述表明，研究入渗时边界条件是较为复杂的，所以，计算方法也较为复杂。
第二节土壤水入渗线性化方程的近似解
在垂直入渗情况下，一维土壤水分运动的基本方程可写作：
化行d（日严（2-5-10）
:
z:
z
如降雨或灌水前的初始含水率（在土壤剖面上含水率均匀分布）为0i，则初始条件为
二（乙0）-（2-5-11）
在地表有一薄水层时，表层含水率等于饱和含水率0s；
在地下水埋深较大时，计算时
段内入渗水不会到达下边界。
为此，下边界土壤含水率不变，等于初始含水率，则边界条件
可以写作以下形式：
E（0,t）“sZ=O,t>
0丿（2-5-12）
日（°
°
t）=©
ZT旳，tA0
由于式（2-5-10）为非线性方程（因为扩散度D（0）及水力传导度k（0）均为待求
平均扩散度D代替D（0）,
含水率0的函数），求解比较困难，为了简化计算，近似地以并以⑴孚严代替字，则式中（25）可简化为
2•
（2-5-13）
（2-5-14）
otj—C0少
D「-N-:
z:
式中（2-5-13）为常系数线性方程，可以用拉普拉斯变换求解对式（2-5-13）采用拉普拉斯变换后可得象函数方程：
2——
dNZP-
—.—U
dz2DdzD
式（2-5-14）的通解为
（2-2-15）
（2-2-16）
（2-2-17）
根据边界条件式（2—5—16）、式（2一5—17）确定常数:
代入式（2-5-15），得象函数的表达式为
N~
--iE
（2-5-18）
二Z,P'
s-
PP
进行逆变换后，得含水率的表达式为
eJ皿〕+e畔erfc
7十Nt
2
l2^Dt丿
l2』Dt/
1—2—2查得。
式（2—5-19）
（2-5-19）
补余误差函数可自表
中D可用下式计算:
5/3-
D」D…f®
入-71。
3
（2-5-20）
若已知D与0的关系式，代入式（2—5—20）积分,
即可求得D。
采用式（2-5-19）求得的土壤剖面上含水率分布如示意图2-5-2所示。
由于地表的入渗强度i二-Dkb〕，为了推
圉2-^-2人邃条件下，剖血上含水率分布示总图
求入渗强度，首先根据二的象函二的表达式求—
（2-5-51）
-Z
2D
'
2
NP
D（4D
在入渗初期，t:
：
0.2-Dy,
N2
相当于P204D
（2-5-21'
P4：
4D
，则式（2-5-21）可近似写成:
宀基7PU1
；
P■.DD、P
经逆变换得:
1
胡_（兀-3）严
z、、二D
入渗初期：
［当D（0）取平均值D时］
）'
（2-5-22）
i=-Dk忑…7
Sz
（2-5-23）
入渗时间较久，即当
P」也之,
204D
t时
cQcQ
）式，贝U=0,所以=0
.1-..
zz
则入渗时间久时，入渗强度（iTks）为
i—Dk忑二k忑（2-2-24）
自式（2—5-23）、式（2—5-24）得入渗速度在时间上的变化过程如图（2-5-3）所示。
第三节Green—Ampt模型的入渗解
Green—Ampt模型5^是1911年提出的一种简化的入渗模型，它是建立在毛管理论基础上的一种入渗模型。
假定土壤是由一束直径不相同的毛管组成，水在土壤入渗过程中，湿
润锋面几乎是水平锋面，且在锋面上各点的吸力水头均为
Sm。
锋面后面的土壤含水率为均
一的，如图（2-5-4）所示。
所以k（0）也为常数，这种模型又称活塞模型。
根据达西定律：
H+S+7
q=—krJ=m——（2-5-25）
H――地面以上水层厚度；
Sm――锋面处土壤负压；
z一锋面推进距离。
图2-5-4Green---AmpiA.
根据水量平衡原理，应等于土体内增加的水量q=dA-,
dt
式（2-5-26）积分:
所以
式（2-5-27）为z〜t关系式，原则上可以求得任何时刻当然也就不难求得该时刻的累计入渗量：
（2-5-28）
W\~-^0z
Ht0时，式（2-5-27）可写作:
t时入渗总量
Z=,2门厂〒0—厂Smt
（2-5-29）
I对t求导，得:
式（2-5-27）'
表明，入渗初期，入渗深度z与•..t成正比。
将
（2-5-30）
当t大时，式（2-5-26）中
Z>
>
H+Sm,因此H—Sm一Z:
1，则由式（2-5-25）可
知:
（2-5-31）
第四节水平入渗条件下的Philip解法
水平入渗条件下的Philip解【51］是一种半解析法，即前半部用解析法，利用博茨曼
（Boltzmann）变换，将偏微分方程转换为常微分方程；
后半部采用迭代计算，求解常微分方程。
由于求解过程中未作过分简化，求得结果较为严密。
、水平入渗的常微分方程推求
水平入渗的基本方程
c9
ft
limy-补
t=0,x0,
t0,x=0,
t0,x—'
（2-5-32）
将式（2-5-32）中基本方程改写为以坐标x（0,t）为变量的方程，根据第二章中方程
（2-5-34）
（2-5-35）
（2-5-36）
（2-5-37）
（2-2-17），在水平入渗时应为
（2-5-33）
采用Boltzmann变换，引入变量从入（0）,且令
v-X（V,t）t2
则XV,t-'
Vt2
m丄V门
.:
t2
将式（2—5-36）、式（2—5—37）代入式（2-5-33）得:
经整理后得微分方程
（2-5-38）
由边界条件已知：
为了求得入〜0的关系式，将式（2—5—38）常微分方程自0i至0积分得：
八dV
-二dv--2D*L（2-5-39）
<
d'
因•二Xt3，即可由入〜B求得B（x,t），从而求得剖面上任何时间，任何距离的含水率分布。
若能实测入（B）,则代入上式可求得D（0）关系曲线。
、迭代计算法
在D（0）已知时，式（2-5—39）为入〜0关系式，它的关系曲线如图2—5—5所示,
00^0i等份成n份，步长为=（00-0i）/n,若等分点按顺序其含
0n-1,0n。
相应各等分点相应入值为入1,入2,…入n,两个相邻等分
1/2,D1+1/2，…D（n-2）
亠_入_r宀（2-5-40）
Dr+1/2取平均值，定义
式（2-5-43）表示入〜0关系曲线所包含的面积。
式（2-5-44）表示0=0r+1/2时，0r+1/2〜0n之间入〜0曲线所包含面积的近似值。
由图可知面积为
r矩形面积后的一小
（2-5-45'
S「+1—-为入〜0的曲线与入r,入r+1/2之间所包含面积减去块面积。
由式（2—5一43）及式（2-5-45）得：
当△0取得足够小时，则Sr+1可以忽略，所以
—J"
则式（2一5一42）可以改写为
J1:
J1丄
J21
=J1-'
1L.：
l
7丄
由图形可知:
A9
2J1
（2-5-46）
（2-5-49）
J1
r-
2"
根据式（2—5—46）和式（2—5—49），若已知J1/2及D（0）则可交替使用此二式，求得入1,入2…入r…，入n-1。
因此，为了求得入〜0关系值必须首先求得J1/2及D〜0关系曲线。
由上述可知J1/2代表了入〜0曲线所包含面积的近似值，由式（2—5—49）可知：
若将面积看作一系列的长为（0一0n）,高为d入的水平方向的矩形，则
od
J1=v一片d'
（2-5-50）
入---和D（0）的函数。
对于初始值J1/2，我们可将D（0）考虑为在整个含水率变化范围内的平均值，由式（2
—5—50）可知：
若用加权平均值来代替上式中D值，以D'
表示：
（2-5-51）
其中等号右端系数2是为使等式衡等，D'
=D'
（若将D'
代替D（0）即可算得该系数）。
将式（2-5-51）写成有限差的形式：
因D'
为常数，直接利用入渗问题线性化解，以，二xt^代人得:
（2-5-58）
（J）=2n.cD2
2/
有了Ji/2就可以由式（2一5—56）和式（2-5-49）重复交替计算求得入〜0关系曲线，但由于J/2是估算值，存在偏差，求得的入〜0曲线必然也有偏差，为了使迭代计算尽快收敛，采用两种不同方法计算Jn-1/2,不断地修正Ji/2，以使两种方法求得的Jn-1/2值足够
的接近，以便得到一个比较精确的入〜0曲线。
第一种方法采用
jFj^J1-■nJ^（2-5-59）[即式（2-5-49）]
勺n一21
第二种方法采用
t
2y2
式中A（y）定义为
（2-5-64）
二2erfy
A（y）的函数表见表（2-5-1）。
对y较大时，A（y）可以渐近级数表示:
（2-5-68）
当厶0足够小时，且假定J1/2是正确的，则JFn-1/2与JSn-1/2应该是相等的，但由于计算
中采用了近似值两者误差为.訂；
，若在所要求精度范围内即可，否则进行第二、第三、……,
N--式-（2—5-70）两边相等（或误差在所要求的精确度范围内）时送代次数。
为了加速收敛，Jl/2还有另一改进公式：
（2-5-71）
P=2,3,,N
由于在其他著作中采用了与上述不同符号,代公式。
这里介绍一下用其他符号表示的Philip的迭
Philip
定义Ir+1/2为
所以,
（2-5-72）
r：
其他以I为参数的各式中均比J为参数的多除一个A0。
—r丄小
22
Dr计算平均扩散度。
（3）根据公式J”2=2n•△B（D'
n）1/2计计算第一次送代初值（Ji/2）1,并由公式（2—5—46）计算入1,由（J1/2）i和入i即可计算（Ji+i/2），再代入公式（2-5-49）求入2,依次类推，最后求得（丿尸“人。
n
（4）由以上计算同时可求得入n-1和Dn-1/2，代入公式（2-5-60）计算（JS1）1。
n一
（5）计算误差△i=（JF1）1-（JS1）1，若△“满足要求，以上计算结束，否则要进行第
匕nP
二、三、…次的送代，直至满足要求为止。
（6）进行多次送代时，重复以上
（2）〜（5）步骤，直至迭代误差△很小，且入1,入2,…、不再随Ji/2改变时为止。
若计算结果不理想，可能是选用△B不合适，则可重新选用
△B值，重复以上步骤进行迭代计算。
11
（7）计算出入〜B关系曲线后，因-=xt2，故X='
t2，对于给定的时间t，由入可计算x值，即得0〜x关系曲线。
也即求得了B在剖面上的分布。
Philip水平入渗解的理论可用于水平土柱法测定土壤水扩散度。
土壤水分在较长的（理
论上说是半无限边界）均质土柱中发生水平运动时，一维土壤水平运动方程如式（2-5-32）
所示，经Boltzmann变换可得[如前述方程（2-5—39）]:
屮ed日fde
——
1心丿
入（0）―--,t的函数。
1当水平土柱首端供水，根据观测水分运移资料，即可绘出入〜0曲线，（■C）二xt_2）,
Q
吕-（Rdr为由初始含水率0i至0x之间入〜0曲线所包含的面积，如图2-5-5所示。
d■'
（）d为在入〜0曲线上对应含水率0X点的斜率。
因此，通过水平土柱试验即可求得土壤
d，
水扩散度D值。
第五节垂直入渗条件下的Philip解法
一维垂直入渗基本方程写成z（0,t）为函数时为
D®
）
&
」
对一类边界定解条件为
■---iz■?
0,t=0,
炉=日0Z=0,t>
0（2-5-75）
0=3ZT°
O,t：
Philip级数解的形式［52］
13
=送・（日）2（2-5-76）
i二
ZRt「二t2•2二t1•3一厂t2■■-
（2-5-77）
相应边界条件为
1航=2卞一…一%=0
9时的位置z。
以下我们以待
（2-5-78）
若求得式（2-5-77）中系数ni（9）值，则可求得在一定定系数法求解ni（9）。
式（2-5-74）可改写为
将式（2-5-80）代入式（2-5-79）归并后得：
丄m：
丄
1二t22二t3二存八i二存=0（2-5-81）
im
式中前四项系数n1,n2,n3,n4分别为
D
d（旳
n3=0得:
（2-5-85）
n4=0得:
閉（日）丿d日
「d日徑凹丫,2屮3（日）_屮2（日）［23（0）八护1（日）丿］小2（日）d=（0）」
（2-5-86）
根据以上各式，若已知n1（0）则可逐步求得外n2（9）,n3（B）,…，即可求得方程（2—5—74）的解。
由于收敛较快，求得前四项就足够精确了。
上述n1（0）,n2（0）,n3（0）,n4（0）都较复杂，n1（0）即为水平入渗时的入，为了求解这些值并不容易，为此Philip采用了递推公式求解。
以下以求得n2为例
说明求解方法。
2-5-6关系福线
DFAC中h可以下式表示:
为了求解系数，我们可以写出与式
（2-5-39）类似的近似系数一般方程，则
I■-df
fdv-:
—一：
（2-5-87）
直d日
a,3—均为9的已知函数。
式（2—5-87）左边可以写作式（2—5
—88），如图2—5—6所示。
1；
fd=「fdv‘2fd：
（2—5
r'
~n雹r
—88）
式（2一5—88）中右端第二个积分为
由图2—5—6可知，该积分近似等于
h，其中h为四边形BCDE的平均高度。
△
9取得足够小时，DE可以看作为一直线，在梯形
（2-5-89）
1FX1f
hr二f「[fri—f「fr13fr
44
将式（2—5—89）、式（2—5—91）代入式（2—5—88）,经整理后得
（2-5-94）
将式（2一5—94）代入式（2—5—93）经整理后得:
式（2-5-95）为Philip的第一个递推公式。
第二个递推公式，将r改为r—l，则式（2-5-94）为
A0
fd—1f
由式（2—5—94），若r=n—I，则
（2-5-101）
当r=n时，fn并不趋于无穷大，则由式（2—5-101）和式（2-5—99）可知，一=0。
以上两递推公式如果已知初始和边界条件即可交替使用，求得f〜0曲线。
如果边界条件为0=00时，f=0,且f0=0，则可由式（2-5-95）求得11/2，由式（2—5—99）又可求得I，当r=1时，代入式（2-5-94）求得f2。
重复以上过程即可求得f〜0曲线。
当f=入，n,^,3等不同项的系数时，都可求得对应的与0关系曲线。
在某一t时，可以求得不同
0的z值，则可求得对应于不同t时的0〜z关系曲线，即求得任何时间剖面上含水率分布0（乙t）。
入渗量计算如下：
渗入土壤的水量应等于流出水量与剖面上含水量变化之和，在半无限情况下，底部含水
率将保持初始含水率，即为0n,导水率也为常量k，因此由下部流出水量
如将土柱剖分为n等份，步长为△z,取横断面为1个单位时，经t时土柱顶部含水率
为B0,对半无限土壤剖面含水量增量近似表示为
精确值可用积分表示：
上式表示0=0n以上B〜z曲线所包含的面积，这面积也可表示为
Wzd）（2-5-103）
式（2-5-102）与式（2-5—103）相加即为t时间内渗入剖面总水量，即
I二.zdrknt（2-5-104）
r
1234
将Philip解的表达式•2t八3t"
4住•…代入式（2—5—104）,并
加以整理得，
I=At2kntBtCt2Dt（2-5-105）
也」6
A二q「drB「冷2rdr
A、B、C、D等均为常数。
式（2—5—105）中右边第一项相当于水平入渗，所以水平入渗量为
（2-5-106）
丨k=At2
入渗率的计算如下：
水进入土壤的速率即为z=0处的水流通量，即为
I以式（2—5—105）代入：
按式（2—5—107）求得的入渗率当t较小时误差较小，但时间长后，计算值与实验资
料偏离较大，Philip对其进行了修正，但修正式并不适用于t小时情况。
Miller和Klute用
Hagan等人的结果，用下式计算与实验值接近。
1丄
iAt2if（2-5-108）
A,if——均为土壤特性常数。
A值大小决定于土壤初始含水率，一般称为吸水率。
if为稳定入渗速度，相当于土壤饱
和时的水力传导度（即渗透系数）ks。
入渗初期，入渗速度很大，ks远较1/2At-1/2为小，可
忽略不计。
随着时间增大，入渗速度迅速减小，在入渗时间很久，上78,式（2—5—108）
中右端第一项趋近于0,i=if,即为稳定入渗速度。
入渗速度在时间上变化如图2-5-3所示。
A值可通过初期入渗总量Ih确定：
lh=At2（2-5-109）
故A=1ht2（2-5-109）
入渗率理论公式中的常数需通过试验确定。
在生产中常采用经验公式计算入渗率i和入
渗量I。
第六节土壤水入渗的经验公式
前节介绍的理论公式，既
升级会员 