朱悅璐，王義成

(1.南昌工程學(xué)院 水利與生態(tài)工程學(xué)院，江西 南昌 330099；2.中國水利水電科學(xué)研究院，北京 100038)

1 研究背景

如何求土壤水分入滲方程的解析解，一直以來是地下水動(dòng)力學(xué)、土壤水動(dòng)力學(xué)以及非飽和入滲理論中的研究熱點(diǎn)[1-3]，現(xiàn)有的土壤水分入滲方程類型已有很多[4]，以Richards方程為例[5]，這些微分方程的顯示表達(dá)式都具有高度非線性，這種非線性包含兩個(gè)方面：其一是微分方程本身非線性，其二是方程中主要參數(shù)非線性(例如非飽和滲透系數(shù)K(θ)、非飽和擴(kuò)散系數(shù)D(θ)都是含水率θ、深度z、時(shí)間t的非線性函數(shù)[6-7])，這導(dǎo)致入滲方程精確的解析解不易求解，很多情況下甚至不能求解。

目前，最常見的處理手段即為采用各種“變換”對(duì)入滲方程進(jìn)行化簡，例如Laplace變換[8]、Boltzmann變換[9]、最小作用原理變分變換等[10]；這些變換化簡的實(shí)質(zhì)，都是一種線性化方案，即通過損失一定的方程精度，來換取實(shí)際工程中可以接受的近似解、半解析解[11]。但回顧文獻(xiàn)容易發(fā)現(xiàn)[12-14]，既往學(xué)者和工程技術(shù)人員，或因?yàn)閷W(xué)科背景所限，在采用各種線性化手段時(shí)，往往會(huì)忽略其前置的數(shù)學(xué)條件，例如在應(yīng)用最廣泛的Laplace變換線性化Richards方程中，對(duì)于變換的前提，即當(dāng)t→∞時(shí)，原函數(shù)發(fā)散速度不超過指數(shù)函數(shù)這一問題，歷史文獻(xiàn)均未進(jìn)行討論，研究者皆默認(rèn)原函數(shù)θ(z，t)滿足該條件或直接認(rèn)為原函數(shù)收斂。這種做法在大多數(shù)的工程案例中是正確的，但在理論上卻可能存在反例。

非飽和入滲過程中，含水率θ(z，t)是深度和時(shí)間的函數(shù)，一般情況下，當(dāng)z增加時(shí)，θ減小，當(dāng)t增加時(shí)，θ增大，如圖1(a)所示；但在理論上可能存在這樣的曲線，即某種極端的條件下，含水率曲線的函數(shù)式不收斂甚至發(fā)散速度較快，即在數(shù)學(xué)上表現(xiàn)為，二元函數(shù)θ(z，t)的值隨t、z的增大而增大，如圖1(b)所示；這會(huì)造成傳統(tǒng)Laplace變換解在表達(dá)式上看似合理但與實(shí)際工程不符的情況，此時(shí)雖然存在數(shù)值解，但由于函數(shù)θ(z，t)的快速發(fā)散性，Laplace變換的前提已不存在，所化簡的線性方程已不適用，即便能求解出某具體數(shù)值，也已失去了物理意義。

圖1 水在不同土質(zhì)、土壤結(jié)構(gòu)中入滲曲線Fig.1 Water infiltration curve in different soil texture and soil structure

基于此，本文針對(duì)非飽和入滲Richards方程Laplace線性變換方案，提出先通過試驗(yàn)，初步擬合可能出現(xiàn)的含水率曲線形態(tài)，再利用數(shù)學(xué)手段判斷含水率曲線斂散性，最終確定是否能進(jìn)行化簡變換等3個(gè)步驟。從數(shù)學(xué)計(jì)算上完善非飽和土力學(xué)理論，從機(jī)理上明確化簡后的方程，是否能代表入滲曲線最真實(shí)形態(tài)，并給出不適合拉普拉斯變換的極端工況，通過反例說明本研究在實(shí)際工程中的意義和必要性。

2 Richards入滲方程的Laplace變換

(1)

式中：θ0為邊界土體含水率；θi為初始土體含水率；t為入滲時(shí)間。在式(1)中以θ(z，t)為原函數(shù)的Laplace變換為[15]：

(2)

根據(jù)Laplace變換性質(zhì)可得式(3)[16-17]：

(3)

對(duì)式(2)第一個(gè)等號(hào)兩邊求z的二階偏導(dǎo)數(shù)可得式(4)：

(4)

(5)

對(duì)于式(5)，在數(shù)學(xué)上已有成熟的解法求得解析解，具體求解步驟見相關(guān)數(shù)學(xué)教材或工程手冊[18]，本文不再贅述。其解為：

(6)

這樣即求出Richards入滲方程經(jīng)過Laplace變換后的近似入滲曲線θ(z，t)，由于erfc(·)函數(shù)值不超過1[19]，結(jié)合θ0、θi定義和運(yùn)算規(guī)則，易得式(6)值恒小于1，含水率計(jì)算結(jié)果收斂。

2.2 經(jīng)典變換存在的問題以上討論了經(jīng)典變換后入滲方程有形如式(6)的近似解，并知道由式(6)中高斯誤差函數(shù)的值域，可約束含水率計(jì)算結(jié)果收斂。但事實(shí)上，在實(shí)際工程中，入滲曲線θ(z，t)是否收斂在事先是未知的，而上述步驟的實(shí)質(zhì)是，在式(2)右端積分符號(hào)內(nèi)，研究者事先默認(rèn)θ(z，t)是連續(xù)且收斂的，以確保后續(xù)運(yùn)算，而由于二階線性常系數(shù)非齊次微分方程解的固定性，總能得出形如式(6)解的形式，這在物理上相當(dāng)于用式(6)這種唯一的曲線形態(tài)，來擬合所有可能出現(xiàn)的入滲情況，顯然，這種做法在數(shù)學(xué)上并不嚴(yán)謹(jǐn)。

因此本文認(rèn)為，在進(jìn)行Laplace變換之前，需進(jìn)行先驗(yàn)試驗(yàn)，估計(jì)某種土樣在實(shí)際入滲過程中曲線的初步形態(tài)，并判斷其斂散性，只有確定式(2)中的先驗(yàn)θ(z，t)函數(shù)發(fā)散速度慢于Laplace變換條件所規(guī)定的指數(shù)函數(shù)時(shí)，才能應(yīng)用Laplace變換對(duì)方程進(jìn)行線性化，其結(jié)果用式(6)擬合才是具有物理意義的。以下小節(jié)，本文將針對(duì)這一問題進(jìn)行討論。

3 幾種常見的入滲曲線斂散性判斷

比較不同函數(shù)發(fā)散程度的快慢，有多種數(shù)學(xué)方案，對(duì)于本研究的實(shí)際問題，含水率函數(shù)值恒大于零，因此可采用“作比取極限法”進(jìn)行判斷[20]：任取兩函數(shù)f(t)、g(t)，當(dāng)t→∞時(shí)，若lim(f(t)/g(t))→∞，則f(t)發(fā)散速度快于g(t)；若lim(f(t)/g(t))→0，則f(t)發(fā)散速度慢于g(t)。

3.1 多項(xiàng)式形態(tài)入滲曲線斂散性當(dāng)θ(z，t)曲線具有多項(xiàng)式形態(tài)時(shí)，以三次多項(xiàng)式為例，其標(biāo)準(zhǔn)方程為：

θ(z，t)=a1z+a2z2+a3z3+b1t+b2t2+b3t3+c1zt2+c2z2t+c3zt+d

(7)

Laplace變換性質(zhì)表明，要式(2)成立，需有θ(z，t)≤Heβt，其中H、β為常數(shù)，這等價(jià)于證明t→∞時(shí)，θ(z，t)發(fā)散速度慢于Heβt。顯然對(duì)于式(7)當(dāng)t→∞時(shí)，上述命題成立，此時(shí)Laplace解滿足變換條件。在這種情況下，應(yīng)用Laplace變換將入滲方程線性化，并求形如式(6)的近似解是可行的。

3.2 指數(shù)形態(tài)入滲曲線的斂散性當(dāng)θ(z，t)具有指數(shù)函數(shù)形態(tài)時(shí)，一種常見的入滲方程表達(dá)式為：

θ(z，t)=zNt

(8)

式中N為常數(shù)。在實(shí)際入滲過程中，若z由工程條件限制取常數(shù)時(shí)，易知當(dāng)t→∞時(shí)，zNt與Heβt具有同階斂散性，因此，當(dāng)入滲深度固定時(shí)，指數(shù)形態(tài)入滲曲線發(fā)散速度不會(huì)超過Heβt，該種情況亦可用Laplace變換化簡方程。

3.3 冪指數(shù)形態(tài)入滲曲線的斂散性若在某些條件下，深度z持續(xù)變化，式(8)由指數(shù)函數(shù)轉(zhuǎn)化為冪指數(shù)函數(shù)(此時(shí)z，t均為自變量)，當(dāng)z→∞，t→∞時(shí)，有：

(9)

這表明入滲曲線θ(z，t)發(fā)散速度快于Heβt，此時(shí)，式(2)Laplace變換最右邊等號(hào)的積分不收斂，因此變換不能成立，在這種情況下，不能用Laplace變換將Richards方程線性化，因?yàn)槠浣庖咽ノ锢硪饬x。

4 反例及對(duì)比

表1 線性化Laplace變換解

將表1計(jì)算結(jié)果繪制成θ～z函數(shù)曲線，如圖2所示。該曲線是一簇下凹的函數(shù)圖形，表示隨深度增加，含水率減少，同一深度時(shí)，隨入滲時(shí)間增加，含水率增大。以圖中藍(lán)色豎線為分隔，在計(jì)算時(shí)間內(nèi)，線性化方程計(jì)算結(jié)果顯示，入滲深度約為15 cm，在z>15 cm區(qū)域，含水率為初始含水率0.037，這表明入滲尚未到達(dá)該區(qū)域。在很多工程案例中，該種形態(tài)的入滲曲線經(jīng)常出現(xiàn)[21-23]，這一計(jì)算結(jié)果“似乎”并無問題。

圖2 線性化方程解的入滲曲線Fig.2 Infiltration curve of linearized equation solution

采用同樣數(shù)據(jù)，以算例中實(shí)際入滲曲線θ(z，t)=z0.8t代入計(jì)算，其結(jié)果如圖3所示，由于冪指數(shù)函數(shù)增長過快，圖3采用單對(duì)數(shù)坐標(biāo)系以便考察，但這并不影響函數(shù)的增減性；可以看出，在計(jì)算時(shí)間內(nèi)，θ(z，t)函數(shù)值隨著z，t增加而增加，所有深度均為滲透區(qū)，該結(jié)果同時(shí)也表明，這種極端情況在數(shù)學(xué)上客觀存在。

圖3 實(shí)際入滲曲線Fig. 3 Actual infiltration curve

為更直觀展示上述結(jié)果，將線性化方程計(jì)算結(jié)果和原方程計(jì)算結(jié)果繪制在三維坐標(biāo)系θ-z-t中，如圖4、圖5所示；此時(shí)可以清晰的看出，兩種情況下，含水率入滲曲面隨z、t的變化，呈現(xiàn)相反的凹凸性，這表明，此時(shí)無論采用何種方式提高線性化方程式(6)本身計(jì)算精度，都已不能有效模擬含水率變化情況，而需要采用其他方法；造成這一現(xiàn)象的數(shù)學(xué)本質(zhì)是傳統(tǒng)計(jì)算方案將非飽和入滲問題化簡為一個(gè)泛定方程，進(jìn)而將這類方程又轉(zhuǎn)化為僅與初始條件和邊界條件有關(guān)的定解問題，而這種不包含滲流場本身信息的微分方程，在某些應(yīng)用場合存在理論瑕疵。

圖5 三維坐標(biāo)下實(shí)際入滲曲線Fig. 5 Actual infiltration curve in 3-dimensional coordinates

上述反例即可證明，Laplace變換化簡Richards方程并非適用于所有情況，事實(shí)上，在實(shí)際工程中，由于土壤性質(zhì)、孔隙結(jié)構(gòu)以及工程條件的不同，真實(shí)非飽和入滲曲線形態(tài)差異很大，所對(duì)應(yīng)入滲曲線方程，在數(shù)學(xué)上能滿足嚴(yán)格收斂條件的亦是少數(shù)，因此在進(jìn)行化簡變換前，本文所提出的原函數(shù)斂散性判斷，不論在數(shù)學(xué)上，還是物理上都十分必要。

5 結(jié)論

本研究對(duì)非飽和入滲Richards方程的Laplace變換過程中，原函數(shù)的斂散性進(jìn)行了討論，并得出以下結(jié)論：

(1)Laplace變換存在性定理，在非飽和土力學(xué)中仍然適用，凡在變換過程中，直接默認(rèn)原函數(shù)收斂或其發(fā)散速度小于Heβt的方案，在數(shù)學(xué)上都不嚴(yán)謹(jǐn)。

(2)在實(shí)際工程中，入滲曲線θ(z，t)具有多項(xiàng)式形式、指數(shù)形式時(shí)，原函數(shù)發(fā)散速度都小于Heβt，可用化簡后的線性化方程求解入滲曲線。入滲曲線若具有冪指數(shù)形式，則原函數(shù)發(fā)散速度快于Heβt，此時(shí)Laplace變換積分發(fā)散，該變換不成立，因此需用其他方法進(jìn)行求解。

(3)凡是先驗(yàn)試驗(yàn)中，含水率曲線θ(z，t)的發(fā)散速度大于Heβt的，Laplace線性化求解方案均不成立，本文中的反例有效驗(yàn)證了這一結(jié)論。