(1.山東電力工程咨詢院有限公司,山東 濟(jì)南 250100; 2.山東黃河勘測設(shè)計研究院,山東 濟(jì)南 250013； 3.山東正元建設(shè)有限責(zé)任公司，山東 濟(jì)南 250100)

研究降雨條件下的非飽和土滲流過程對許多巖土工程或環(huán)境巖土工程的實(shí)踐有重要意義，例如，自然或人工邊坡穩(wěn)定性分析[1-4]、礦床或城市垃圾填埋場設(shè)施的覆蓋層設(shè)計[5-6]、涉及膨脹土和濕陷性土的工程問題[7-8]，等等。解決這些問題的關(guān)鍵和難點(diǎn)在于計算降雨入滲過程在非飽和土體中引起的滲流場[9]。

Richards方程是描述非飽和土體中瞬態(tài)滲流的控制方程。針對該方程的求解，以往的研究多采用數(shù)值模擬手段給出數(shù)值解[10-12]，而解析方法的研究相對較少。實(shí)際上，與數(shù)值解相比，解析解更清晰簡潔，能更直觀地顯示定解問題對邊界條件的響應(yīng)，且便于對各影響因素進(jìn)行參數(shù)敏感性分析，并能對數(shù)值解的結(jié)果進(jìn)行驗(yàn)證[13]。Srivastava 等推導(dǎo)了均質(zhì)、非均質(zhì)土體降雨條件的一維瞬態(tài)入滲解析解[14]，基于此模型，詹良通[13]和李寧[9]等人分別推導(dǎo)了考慮坡度影響和降雨形式的滲流解析解。Chen 等[15-16]提出了非穩(wěn)恒降雨條件下的入滲解析模型。Huang等[17]和Ghotbi等[18]也分別基于Fourier積分變換和同倫分析方法(HAM)提出了降雨入滲過程的解析解。在初始狀態(tài)的處理方面，上述研究或?qū)⒋嬖谇捌诮涤瓴⑦_(dá)到穩(wěn)態(tài)入滲時刻作為初始條件[9,13-14,16]，或假定初始含水率(水頭)分布為線性或特定函數(shù)[15,17-18]，而根據(jù)實(shí)際情況，初始狀態(tài)是可以選擇在任意時刻，相應(yīng)含水率(水頭)分布往往是不規(guī)則的。誠如Zhan 等[19]所述，前期降雨所形成的初始狀態(tài)對水頭分布的影響甚至大于后續(xù)降雨。朱偉等指出在非穩(wěn)定滲流中，能不能正確地考慮初始水分量將會對計算結(jié)果產(chǎn)生非常大的影響[20]。在采用的數(shù)學(xué)算法方面，前人多借鑒Ozisik (1980)關(guān)于偏微分方程的解法[21]，在解析表達(dá)式中存在復(fù)雜的積分變換項(xiàng)[9，15-17]。該類方法在理論上是先進(jìn)的，但數(shù)學(xué)算法稍顯復(fù)雜，不利于直接應(yīng)用。因此，需要進(jìn)一步研究任意初始狀態(tài)的入滲過程解析解，并考量初始狀態(tài)的影響。

本文采用指數(shù)函數(shù)描述非飽和土體的土水特征曲線及其滲透系數(shù)隨水頭的變化，對任意初始狀態(tài)的非飽和土入滲過程解析解進(jìn)行了推導(dǎo)。最后基于所推導(dǎo)模型，分析初始狀態(tài)對含水率分布的影響。

1 理論推導(dǎo)

1.1 基本假設(shè)

(1) 土體為均質(zhì)，非飽和滲流由基質(zhì)吸力差(水頭差)決定，不考慮溶質(zhì)吸力的影響。

(2) 不考慮非飽和土的回滯特性。

(3) 地下水位始終穩(wěn)定，不受土表降雨入滲過程的影響。

(4) 降雨過程雨強(qiáng)穩(wěn)定，即為均勻降雨。

1.2 控制方程

Richards為非線性偏微分方程，目前它已衍生出3種基本格式，即壓力水頭格式(h-form)、混合格式(mixed form)和體積含水率格式(θ-form)，其中混合格式被證實(shí)具有嚴(yán)格遵守質(zhì)量守恒的特性，可表示為

(1)

式中，h為孔隙水壓力水頭；θ為體積含水率；k為滲透系數(shù)；t為時間；z為縱坐標(biāo)，正方向向下(如圖1所示)。

圖1 入滲過程示意Fig.1 Schematic illustration of infiltration process

Yanful和Mousavi[22]指出，在非飽和土滲流中相較于大的的壓力水頭差，式(1)中的重力項(xiàng)可以忽略。同時，忽略重力項(xiàng)可以極大簡化推導(dǎo)過程，避免解析解中復(fù)雜的積分或級數(shù)項(xiàng)[23]。因此，Richards公式可以簡化為

(2)

參照文獻(xiàn)中常用的指數(shù)函數(shù)形式[13-17,24]，土水特征曲線和滲透系數(shù)方程可描述為

θ(h)=θr+(θs-θr)eαh

(3a)

k(h)=kseαh

(3b)

式中，θs和θr分別為飽和和殘余體積含水率；ks為飽和滲透系數(shù)；α為減飽和系數(shù)。將式(3)代入式(2)可得

(4)

為了簡便計算，定義土的擴(kuò)散系數(shù)D和無量綱的體積含水率Θ分別為

(5a)

Θ=(θ-θr)/(θs-θr)

(5b)

將式(3a)、(5a)和(5b)代入式(4)，并進(jìn)一步化簡，可得：

(6)

1.3 入滲過程解析解

降雨強(qiáng)度與土體飽和滲透系數(shù)的相互關(guān)系會對入滲過程產(chǎn)生影響。對于降雨強(qiáng)度小于飽和滲透系數(shù)的情況，即q≤ks，降雨將全部入滲，土表不會形成積水，也即整個過程均處于定流量入滲階段。而對于降雨強(qiáng)度大于飽和滲透系數(shù)的情況，即q>ks，其起始為定流量入滲階段，降雨一段時間后，表層土體達(dá)到飽和，降雨不能完全入滲，土表存在積水，定流量入滲階段轉(zhuǎn)化為土表積水階段。

1.3.1定流量入滲階段

在定流量入滲階段，上邊界為恒定的流量邊界，也即降雨強(qiáng)度。此時的初始和邊界條件可表示為

Θ(z,t)|z=L=1

(7a)

(7b)

Θ(z,t)|t=0=f(z)

(7c)

式中，L為地下水位深度；q為降雨強(qiáng)度；f(z)為任意初始狀態(tài)含水率分布函數(shù)。式(7a)表示在地下水位處，下邊界土體始終飽和，Θ始終為1；式(7b)表示上邊界為恒定流量邊界條件；式(7c)表示初始狀態(tài)。

首先求穩(wěn)態(tài)入滲時的含水率分布，設(shè)為w(z),則w(z)滿足：

(8a)

w(z)|z=L=1

(8b)

-D(θs-θr)w′(z)|z=0=q

(8c)

聯(lián)立式(8)可以求得w(z)的表達(dá)式為

(9)

由于原問題的邊界條件為非齊次(如式(7))，因此為方便求解，引入一變量使其齊次化，令：

(10)

將式(9)和式(10)代入式(6)和式(7)，可得：

(11a)

(11b)

(11c)

(11d)

這里，式(11a)為控制方程，式(11b)和式(11c)為邊界條件，式(11d)為初始條件。根據(jù)附錄中的推導(dǎo)過程，式(11)的解可以寫為

(12)

由于式(12)中存在積分項(xiàng)，較難定量求得。為簡便應(yīng)用，用一個二次函數(shù)來擬合初始含水率分布，以消除積分變量，簡化計算。由于z=L時，Θ始終為1，因此令：

(13)

式中，a和b為擬合參數(shù)。將式(9)和式(13)代入式(12)，并進(jìn)一步化簡，可得定流量入滲階段土體含水率瞬態(tài)分布的解析解表達(dá)式：

(14)

1.3.2土表積水階段

對于q>ks的情況，假設(shè)土表積水開始時間為tp。不考慮積水造成的土表水頭堆積的影響，則當(dāng)t>tp時，上邊界條件為恒定含水率邊界條件。

首先對于積水時刻tp，可令式(14)等于1求得。設(shè)此時的含水率分布為p(z)，則p(z)的表達(dá)式為

(15)

當(dāng)t>tp時，定義一個新的時間變量T=t-tp，將T代入式(6)可得：

(16)

在土表積水階段，初始和邊界條件可以表示為

Θ(z,T)|z=L=1

(17a)

Θ(z,T)|z=0=1

(17b)

Θ(z,T)|T=0=p(z)

(17c)

根據(jù)Carslaw和Jaeger(1959)[25]的解法，式(16)在邊界條件式(17)下的解為

(18)

由于式(18)含p(z)函數(shù)未解積分項(xiàng)，若直接將式(15)代入求解，算法復(fù)雜，且不能得到其解析表達(dá)式，限制其直接應(yīng)用。因此，我們用二次函數(shù)來代替式(15)，以實(shí)現(xiàn)簡化。因?yàn)楫?dāng)T=0時，在z=0和z=L處，Θ=1，所以二次函數(shù)的表達(dá)式可寫為

(19)

其中，c為擬合參數(shù)，利用式(19)擬合式(15)來可以求得其值。將式(19)代入式(18)可求得存在土表積水階段的含水率分布解析解：

(20)

2 算例分析和討論

本文針對降雨強(qiáng)度和飽和滲透系數(shù)的關(guān)系，分別進(jìn)行如下兩個算例討論。如表1所示，在算例1中降雨強(qiáng)度q小于飽和滲透系數(shù)ks；而在算例2中，q大于ks。對算例1和算例2，均進(jìn)行兩種初始條件(線性和非線性含水率分布)的分析。

表1 邊界條件和土性參數(shù)Tab.1 The boundary conditions and soil properties

圖2為降雨強(qiáng)度小于飽和滲透系數(shù)時，不同時刻土體內(nèi)含水率的分布，其值均由式(14)計算得出。圖2(a)中，a和b分別賦值0和1來擬合線性初始含水率分布，而圖2(b)中，a和b分別為2和-1.5。比較圖2(a)和(b)，不難發(fā)現(xiàn)初始狀態(tài)對含水率分布產(chǎn)生非常大的影響。另外，圖2(b)的初始狀態(tài)可認(rèn)為存在較強(qiáng)的前期降雨，由于后期降雨強(qiáng)度小于飽和滲透系數(shù)，土體內(nèi)部滲流較快，使得表層土體含水率在t=5 h小于初始含水率。

圖3展示了降雨強(qiáng)度大于飽和滲透系數(shù)時，不同時刻土體內(nèi)含水率的分布。其初始狀態(tài)設(shè)定與圖2情況相同。對圖3(a)和(b)，其積水時刻tp分別為9.69 h和9.39 h，土表邊界在ttp時上邊界始終維持在飽和含水率，含水率分布由式(20)計算求得，且以tp時刻的含水率分布為初始條件。對應(yīng)圖3中的兩種

圖2 算例1中不同時刻土體內(nèi)含水率分布Fig.2 Water content redistribution during infiltration for Case 1

圖3 算例2中不同時刻土體含水率分布Fig.3 Water content redistribution during infiltration for Case 2

情況，式(20)中的擬合參數(shù)c分別為1.53和1.68。隨著時間推移，含水率分布的變化幅度逐漸變小，并最后趨于穩(wěn)態(tài)。由圖3可見，初始條件不僅影響含水率的大小，而且影響含水率分布形態(tài)。

3 結(jié) 論

(1) 基于指數(shù)函數(shù)來描述土水特征曲線和滲透系數(shù)隨水頭的變化，得到降雨強(qiáng)度大于或小于飽和滲透系數(shù)的入滲過程解析解。該模型可以考慮任意初始狀態(tài)，并且避免了復(fù)雜的數(shù)學(xué)算法表達(dá)，無需借助編程來實(shí)現(xiàn)數(shù)學(xué)計算，易于應(yīng)用。計算結(jié)果表明初始狀態(tài)對入滲過程的計算結(jié)果有顯著影響。

(2) 提出的模型能夠揭示特定情況下的入滲過程機(jī)理，形式簡單，能克服數(shù)值模擬的復(fù)雜數(shù)學(xué)算法及計算收斂問題，便于考察入滲過程的機(jī)理。

(3) 但同時需要考慮到，解析解的得出是基于一定的假設(shè)和簡化。實(shí)際應(yīng)用中，應(yīng)將解析解手段和數(shù)值模擬方法結(jié)合應(yīng)用。