开放科学(资源服务)标识码(OSID): 
0 引言
在非饱和入渗领域经典教材或文献中,最广泛的方法是采用Laplace变换求解上述问题[12⇓-14],即通过变换将非线性Richards方程转化为线性微分方程并求解,再通过逆变换得出结果;回顾以往研究成果可以发现,该方法虽已相对成熟,但仍有几个问题一直没有得到充分讨论:①对于该解析解应用的前提,即入渗函数(原函数)发散速度不能超过指数函数这一问题,以往研究皆进行了化简,默认入渗曲线收敛,这种做法在绝大多数工程案例中是正确的,但在理论上却可能出现反例。②在Laplace逆变换这一最为关键步骤上,文献大多进行了省略,而指向查表(例如文献[12]—文献[14]),事实上即便是大的土木类领域,也鲜有对这一逆变换在原理上的解释,它的缺失会导致岩土工程技术人员、尤其是青年学者在看问题时出现“重应用,轻原理”的现象,这对学科的发展是不利的。③在实际计算中,当渗透系数较小时,入渗曲线会出现一个随着深度增加而增大的区间,这种与实际工况不符的现象,在既往文献中也并无指出。
基于此,本文对上述3个可能遗留的问题进行补充讨论:对于问题①,列举了不满足变换条件的一个极端反例,指出在应用Laplace变换求解Richards方程时,先对入渗曲线形态进行估计的必要性;对于问题②,针对Laplace逆变换步骤中历史文献的跳跃性缺失,进行了详细的推导和补充,闭合了土工类文献在Richards方程求解理论中长期缺失的一环;对于问题③,通过实际计算,讨论不同渗透系数下入渗曲线的适用性,并从物理背景和数学理论上分析该解析解所存在的“先天不足”。这些补充和讨论,或可为完善非饱和入渗理论提供一定的借鉴。
1 Laplace正变换的前提及反例
1.1 入渗方程求解的前提
扩散常数为 且不带重力项的Richards方程及其边界条件为
式中:θ为含水率;θ0为边界土体含水率;θi为初始土体含水率;t为入渗时间;z为深度。
根据Laplace变换基本概念,式(1)中以θ(z,t)为原函数,并对其中变量t进行Laplace变换的结果为
式中: (z,p)为象函数;Lt为Laplace变换符号,表示对原函数θ(z,t)中变量t进行变换;p为参数。该式是所有后续推导的基础,通过一系列计算[15],其最终解的形式为
式中erfc(.)为余误差函数。
式(2)在数学上有十分严格的存在条件,即原函数θ(z,t)需收敛或发散速度小于指定函数ept。但如前所述,已有研究皆未讨论这一条件,均默认非饱和入渗曲线的函数表达式是收敛的,事实上这种假设并不全面,稍加考察即可发现,同一时间t,含水率θ随着深度的增加而减小;同一深度z,含水率θ随着时间t的增加而增大。二元函数θ(z,t)同时在t➝∞和z➝∞时,其敛散性实质上是未知的。
1.2 反例
某土体具有幂指数形式入渗曲线θ(z,t)=zNt,其中N为常数,显然其具有上凸的形态,如图1所示。
当t➝∞,z➝∞时有
式(4)表明,该曲线发散速度快于Laplace变换指定函数ept,变换不成立,应采用其他方法求解。但此时若仍以式(3)进行计算,由于该解析解形式的固定性,它最终会收敛为下凹形态,如图2所示。
由图2可以看出,由于解析解和实际曲线的凹凸性相反,此时不论如何提高式(3)精度,都已无法有效模拟含水率真实变化的状态。该现象的数学本质是解析解仅与边界条件和初始条件有关,只能用来模拟某一类入渗特征的曲线,当实际入渗曲线特征与解析解所具备的特征不同时,将会产生较大偏差,这是Laplace变换求解Richards方程的一个缺陷,因此在计算前,对入渗曲线敛散性进行判断,具有实际意义。
2 Laplace逆变换缺失步骤的补充
2.1 象函数微分方程的解
将式(1)经过Laplace变换,可得一个关于象函数 (z,p)的二阶线性常系数非齐次微分方程,该微分方程有固定解的形式为
2.2 逆变换求解原函数
本文逆变换求解涉及到积分变换和变量代换等手段,具体步骤如下:
对余误差函数erfc(z/2 )做拉普拉斯变换,即
对式(6)右端采用分部积分法得
$ \begin{aligned} L\left[\operatorname{erfc}\left(\frac{z}{2 \sqrt{t}}\right)\right] & =-\left.\frac{1}{p} \mathrm{e}^{p t} \operatorname{erfc}\left(\frac{z}{2 \sqrt{t}}\right)\right|_{0} ^{\infty}+ \\ \frac{1}{p} \int_{0}^{\infty} \mathrm{e}^{-p t} \mathrm{efrc}^{\prime}\left(\frac{z}{2 \sqrt{t}}\right) \mathrm{d} t & =\frac{1}{p} \int_{0}^{\infty} \mathrm{e}^{-p t} \operatorname{erfc}^{\prime}\left(\frac{z}{2 \sqrt{t}}\right) \mathrm{d} t \end{aligned}$
其中,
$ \operatorname{erfc}^{\prime}\left(\frac{z}{2 \sqrt{t}}\right)=\frac{2}{\sqrt{\pi}} \frac{\partial}{\partial x}\left(\int_{\frac{z}{2 \sqrt{t}}}^{\infty} \mathrm{e}^{-x^{2}} \mathrm{~d} x\right)=\frac{z}{2 t \sqrt{\pi t}} \mathrm{e}^{-\frac{2^{2}}{4 t}}$
式中x为余误差函数中积分变量。
结合式(6)—式(8)可得
$ L\left[\operatorname{erfc}\left(\frac{z}{2 \sqrt{t}}\right)\right]=\frac{z}{2 p \sqrt{\pi}} \int_{0}^{\infty} \mathrm{e}^{-p t-\frac{z^{2}}{4 t}} \frac{\mathrm{~d} t}{t \sqrt{t}} 。$
对式(9)等号右边积分变量t做t➝ 代换(不影响等号左边)得
$ L\left[\operatorname{erfc}\left(\frac{z}{2 \sqrt{t}}\right)\right]=\frac{z}{2 p \sqrt{\pi}} \int_{0}^{\infty} \mathrm{e}^{-p t^{2}-\frac{z^{2}}{4 t^{2}}} \frac{\sqrt{t} \mathrm{~d} t}{t}$
同理,对式(10)等号右边积分变量t做t➝t2代换得
$ L\left[\operatorname{erfc}\left(\frac{z}{2 \sqrt{t}}\right)\right]=\frac{z}{p \sqrt{\pi}} \int_{0}^{\infty} \mathrm{e}^{-\frac{p}{t^{2}}-\left(\frac{z t}{2}\right)^{2}} \mathrm{~d} t$
此时,问题转化为求解式(11)中等号右边积分,以下对该问题进行求解,对于 、t存在如下所示关系,即
将式(12)整理后得
结合式(12)、式(13),对特殊积分 dt做如下变换,即
$ \begin{array}{c} \int_{0}^{\infty} \mathrm{e}^{-\frac{1}{t^{2}-t^{2}}} \mathrm{~d} t=\int_{0}^{\infty} \mathrm{e}^{-\left(\frac{1}{t} \pm t\right)^{2} \pm 2} \mathrm{~d} t=\mathrm{e}^{ \pm 2} \int_{0}^{\infty} \mathrm{e}^{-\left(\frac{1}{t} \pm t\right)^{2}} \mathrm{~d} t= \\ \mathrm{e}^{2} \int_{0}^{\infty} \mathrm{e}^{-\left(\frac{1}{t}+t\right)^{2}} \mathrm{~d} t=\mathrm{e}^{-2} \int_{0}^{\infty} \mathrm{e}^{\left.-\left(\frac{1}{t}-t\right)\right)^{2}} \mathrm{~d} t \quad \circ \end{array}$
对于 与t,再考虑如下关系,即
结合式(14)和式(15)得
$ \begin{array}{c} \int_{0}^{\infty} \mathrm{e}^{-\frac{1}{t^{2}}-t^{2}} \mathrm{~d} t=\frac{1}{2} \mathrm{e}^{2} \int_{0}^{\infty} \mathrm{e}^{-\left(\frac{1}{t}+t\right)^{2}} \mathrm{~d}\left(t+\frac{1}{t}\right)+ \\ \frac{1}{2} \mathrm{e}^{-2} \int_{0}^{\infty} \mathrm{e}^{\left.-\left(\frac{1}{t}-t\right)\right)^{2}} \mathrm{~d}\left(t-\frac{1}{t}\right) \end{array}$
此时,注意到有以下4个极限存在:
$ \begin{array}{c} \lim _{t \rightarrow 0}\left(t+\frac{1}{t}\right)=+\infty, \quad \lim _{t \rightarrow \infty}\left(t+\frac{1}{t}\right)=+\infty \\ \lim _{t \rightarrow 0}\left(t-\frac{1}{t}\right)=-\infty 、 \quad \lim _{t \rightarrow \infty}\left(t-\frac{1}{t}\right)=+\infty \end{array}$
设x=t± ,并将式(17)代入式(16)得
$ \begin{array}{l} \int_{0}^{\infty} \mathrm{e}^{-\frac{1}{t^{2}}-t^{2}} \mathrm{~d} t=\frac{1}{2} \mathrm{e}^{2} \int_{+\infty}^{+\infty} \mathrm{e}^{-x^{2}} \mathrm{~d} x+ \\ \frac{1}{2} \mathrm{e}^{-2} \int_{-\infty}^{+\infty} \mathrm{e}^{-x^{2}} \mathrm{~d} x=\frac{\sqrt{\pi}}{2} \mathrm{e}^{-2} \end{array}$
同理可得
$ \int_{0}^{\infty} \mathrm{e}^{-\frac{p}{t^{2}}-t^{2}} \mathrm{~d} t=\frac{\sqrt{\pi}}{2} \mathrm{e}^{-2 \sqrt{p}} \quad ;$
再同理可得
$ \int_{0}^{\infty} \mathrm{e}^{-\frac{p}{t^{2}}-\left(\frac{z t}{2}\right){ }^{2}} \mathrm{~d} t=\frac{\sqrt{\pi}}{z} \mathrm{e}^{-z \sqrt{p}}$
将式(20)代入式(11)得
$ L\left[\operatorname{erfc}\left(\frac{z}{2 \sqrt{t}}\right)\right]=\frac{1}{p} \mathrm{e}^{-z \sqrt{p}}$
同理可得
$ L\left[\operatorname{erfc}\left(\frac{z}{2 \sqrt{\bar{D} t}}\right)\right]=\frac{1}{p} \mathrm{e}^{-z \sqrt{\frac{p}{\bar{D}}}} \circ$
至此,式(22)表明,象函数 /p的原函数为erfcz/(2 ),故式(5)中象函数 (z,p)经Laplace逆变换后的原函数θ(z,t)为
$ \theta(z, t)=\left(\theta_{0}-\theta_{\mathrm{i}}\right) \operatorname{erfc}\left(\frac{z}{2 \sqrt{\bar{D} t}}\right)+\theta_{\mathrm{i}}$
该式即为应用Laplace变换所求的Richards方程解,当其他条件不变,考虑重力项且渗透系数为常数 时,式(1)的形式变为
应用本小节同样的方法,可求得第一类边界条件一般形式Richards方程的解为
3 解析解的工程适用性
式(25)在边界条件和初始条件已知情况下,仅是深度z和时间t的函数,应用较为方便,但在实际计算时,会出现一些不符合物理规律的入渗现象。
3.1 计算中易出现的问题
某工况下,土体初始含水率为0.037,边界含水率为0.64,扩散系数 为0.42 cm2/min,渗透系数 为0.02 cm/min,应用式(25)分别计算不同时段t、不同深度z处的含水率,结果如表1所示。限于篇幅,表中只罗列t=70 min时刻、5~80 cm深度处含水率计算数据以供验证,其余数据类似,不再赘述。
表1 Richard方程的Laplace变换解Table 1 Laplace transform solution of Richards Equation |
| θi | θ0 | /(cm2·min-1) | /(cm·min-1) | z/cm | t/min | erfc(v1) | erfc(v2) | 含水率计算值 |
|---|---|---|---|---|---|---|---|---|
| 0.037 | 0.64 | 0.42 | 0.02 | 5 | 70 | 0.638 | 0.404 | 0.742 |
| 0.037 | 0.64 | 0.42 | 0.02 | 10 | 70 | 0.262 | 0.137 | 0.336 |
| 0.037 | 0.64 | 0.42 | 0.02 | 15 | 70 | 0.076 | 0.032 | 0.126 |
| 0.037 | 0.64 | 0.42 | 0.02 | 20 | 70 | 0.015 | 0.005 | 0.055 |
| 0.037 | 0.64 | 0.42 | 0.02 | 25 | 70 | 0.002 | 0.001 | 0.043 |
| 0.037 | 0.64 | 0.42 | 0.02 | 30 | 70 | 0.001 | ≈0 | 0.042 |
| 0.037 | 0.64 | 0.42 | 0.02 | 35 | 70 | ≈0 | ≈0 | 0.041 |
| 0.037 | 0.64 | 0.42 | 0.02 | 40 | 70 | ≈0 | ≈0 | 0.037 |
| 0.037 | 0.64 | 0.42 | 0.02 | 45 | 70 | ≈0 | ≈0 | 0.037 |
| 0.037 | 0.64 | 0.42 | 0.02 | 50 | 70 | ≈0 | ≈0 | 0.037 |
| 0.037 | 0.64 | 0.42 | 0.02 | 80 | 70 | ≈0 | ≈0 | 0.037 |
注:v1=(z- t)/(2 ),v2=(z+ t)/(2 )。 |
为验证解析解在不同分区的适用性,在Ⅰ区选5 cm处不同时刻实测值与计算值进行对比,实际入渗时,当t从10 min增大至90 min时,含水率由0.13增大至0.61,而计算值由0.16增大至0.80,最大误差为32%,从数学角度来看,这一误差值是不可接受的,因此认为解析解在Ⅰ区不适用;同样的,在Ⅱ区选15 cm处不同时刻实测值与计算值对比,此时实测值与计算值最大误差仅为5.8%,解析解在Ⅱ区有良好的适用性。计算结果如图4所示。
3.2 物理现象及归因分析
从图3中可以看出,Ⅰ区中土层上部含水率小于下部含水率,存在“上小下大”这种反常现象,与实际工程认知不符,在何种条件下会出现上述现象,本文对此进行如下验证:
上述现象的成因有其内在物理背景,这是由于在式(25)的解析解中,渗透系数 为常数(扩散系数 也如此),如图6(a)所示,此时t1时刻A点、B点和t2时刻C点的渗透系数相同,导致土体的渗透性不随含水率增加而增大,短时间内,持续增大的含水率和较小的渗透性导致水分难以下渗,从而出现形如Ⅰ区的反常现象。
而在实际下渗中,渗透系数是含水率θ的函数,一般随含水率增加而增大,如图6(b)所示,此时在B点,由时刻t1➝t2时,含水率也从K(θ1)➝K(θ2),随着渗透性增加,水分进一步下渗,湿润锋继续向下运动,这将不会出现“积水”的不合理区域,整个土层符合水分运动规律。
3.3 数学本质与“天然缺陷”
在明确物理背景后,上述问题似乎很容易得到解决,正如前小节所述,只需将式(25)解析解中 、 按照实际水分运动规律替换为K(θ)和D(θ)即可得到“满意”结果。事实上,这种方案是行不通的。因为对于式(1)或式(24),它能进行Laplace变换的本质,就是在变换前,首先假定了渗透系数 和扩散系数 是与积分变量θ无关的常数,正是有了这种假定,才能得出形如式(23)或式(25)的解析解。因此,若在结果中任意代换K(θ)和D(θ),将会产生逻辑错误,这实际上是Laplace变换方案求解Richards方程解析解的“天然缺陷”,无法从根本上消除。
因此,在使用该解析解进行单次计算时, 、 选定后就不能再变动(特别说明,本文3.2节中选定了4种不同的渗透系数,但分别属于4次独立的计算,在单次计算中,参数并未变化),这就使得参数的选择在实际应用中尤为重要,早期研究者曾提出过如下经验公式,即:
式中:θs为饱和含水率;k(θs)、k(θi)分别为饱和渗透系数和初始渗透系数。但为何会有这些经验公式存在于解析解下,一直以来并无系统的背景解答,本文至此给出了完整的阐述。
4 结论
本文系统解答了Richards方程Laplace变换解在理论和应用过程中,前人讨论不充分的3个问题,即入渗曲线敛散性问题、关键求解步骤缺失问题、解析解适用性问题。从数学和物理本质上补全解析解可能存在的瑕疵,本文的主要结论如下:
(1)Laplace正变换在数学中有严格的定义,即原函数发散速度要小于指定指数函数,这在非饱和入渗领域依然有效,当不满足变换条件仍继续计算时,会出现入渗曲线形态与实际不符的现象,它是解析解固有的函数类型造成的,本文反例有力证明了这一点,因此计算前必须首先判断原函数敛散性。
(2)既往研究对于求解步骤,在Laplace正变换和微分方程求解2个阶段论述比较充分,但对关键的逆变换过程,大多省略跳跃,这是由于这一逆变换并不直观,且无通用求解方案,本文给出一种利用积分变换和换元代换的求解方法,完善了这一经典解析解的求解全过程。
(3)该解析解在实际工程中,在渗透系数较小的情况下,会出现土体下层含水率大于上层的悖论,这是由于Laplace变换在求解前,就假定了渗透系数和扩散系数为常数所引起的,这属于该解的天然缺陷,不能通过后期修正渗透系数来调整精度,在使用时需特别注意。