基于差分法及試驗聯(lián)合確定非飽和土滲透系數(shù)的方法

江 明 王世梅 張?zhí)m慧 劉 凡

(三峽大學(xué) 土木與建筑學(xué)院， 湖北 宜昌 443002)

非飽和土滲透系數(shù)是非飽和土研究的重要參數(shù)之一，對于非飽和土滲透系數(shù)的求解，國內(nèi)外學(xué)者從試驗和理論兩個角度展開了研究．試驗方法[1-2]耗時過長，且其測得的非飽和土滲透系數(shù)的解的范圍相對有限．因此大多數(shù)學(xué)者更傾向于理論研究，F(xiàn)rendlund等[3]在土水特征曲線的基礎(chǔ)上建立了非飽和土滲透函數(shù)經(jīng)驗?zāi)Ｐ?，缺點是計算過程復(fù)雜；孫大松等[4]建立了用分維值和進氣值表示的土水特征曲線和滲透系數(shù)分形模型；Campbell模型[5]、Gardner模型[6]是以土的粒徑分布、孔隙率等特性指標確定的土壤轉(zhuǎn)換函數(shù)來間接估計非飽和土滲透系數(shù)的經(jīng)驗?zāi)Ｐ停欢鲜瞿Ｐ投疾痪哂衅毡樾?，目前Van Genuchten[7]提出的VG模型在預(yù)測非飽和土滲透系數(shù)方面?zhèn)涫軐W(xué)者青睞．該模型簡單且具有普遍適用性，其不足之處在于VG模型是經(jīng)驗公式，與工程實際略有偏差，且其表達的是滲透系數(shù)與基質(zhì)吸力之間的關(guān)系，忽視了干濕循環(huán)作用下的滯后效應(yīng)[8]，相關(guān)研究表明[9]，以體積含水量為自變量求得的滲透系數(shù)不受滯后效應(yīng)的影響．因此本文基于有限差分法和一維垂直滲流試驗，建立了一種以體積含水量為自變量求解非飽和土滲透系數(shù)的新方法，為今后非飽和土體滲透系數(shù)的預(yù)測提供一定的依據(jù)．

1 基于差分法的公式推導(dǎo)

1.1 一維垂直滲流控制方程

非飽和土水分運動基本方程的建立需要滿足廣義達西定律和質(zhì)量守恒定律兩個條件．廣義達西定律是指單位時間內(nèi)通過單位面積的水量即土體中液相的流速與其水力梯度成正比．

(1)

式中，V為水的流速，Kθ為滲透系數(shù)，ψ為總土水勢，?ψ/?z為z方向的水力梯度．

理想狀態(tài)下飽和土體中各位置體積含水量和基質(zhì)吸力均相同，故飽和土體中的滲透系數(shù)K始終為常量．而非飽和土體中的滲透系數(shù)Kθ則隨基質(zhì)吸力或體積含水量的變化而變化．

非飽和土瞬態(tài)流形式中的水分流動滿足質(zhì)量守恒定律，可描述為單位體積土流入水的質(zhì)量等于儲存在該單位體積土中水的質(zhì)量和流出水的質(zhì)量之和．微元體內(nèi)水的儲存量可以用單位時間內(nèi)的土體體積含水量的變化率表示．故瞬態(tài)水流可表示為：

(2)

式中，θ為體積含水量，t為時間．將(1)式代入(2)式可得一維非飽和土水分運動方程：

(3)

對非飽和土來說，溶質(zhì)勢和溫度勢極其微小，總土水勢ψ忽略溶質(zhì)勢和溫度勢對土壤水分運動的影響．如果以原點為重力勢能的零點，由于坐標正向向下，則重力勢可表示為-z，總土水勢為：

ψ=h-z

(4)

式中，h為基質(zhì)勢；-z為重力勢．

令

(5)

則式(3)變換為如下形式：

(6)

(6)式表達的即是一維垂直滲流控制方程，該方程僅僅只是適用于土體骨架不變形，且沒有考慮生物或者化學(xué)作用對水流動影響的情況．

1.2 一維垂直滲流控制方程離散

由泰勒級數(shù)可得：

f(x+Δx)≈

(7)

f(x-Δx)≈

(8)

由(7)式得向前差商公式：

(9)

由(8)式得向后差商公式：

(10)

聯(lián)立(7)、(8)式得中心差商公式：

(11)

利用有限差分的思想對(6)式進行離散．將研究區(qū)域劃分為矩形網(wǎng)格(如圖1所示)，設(shè)i代表位置，j代表時間，時間步長為Δt，空間步長為Δz，對(6)式空間i和時間j兩個維度進行離散．

圖1 有限差分網(wǎng)格劃分

因為θ=θ(t，z)，對于公式(6)，令

(12)

則(6)式化簡為：

(13)

對(13)式進行中心差商：

其中：

(14)

(15)

(16)

(17)

將式(14)～(17)代入式(13)得滲流控制方程的差分格式(18)：

(18)

將式(18)寫成非飽和土滲透系數(shù)的迭代格式(19)：

(19)

將式(19)寫成簡化模式(20)，p，q，r由(19)式因子算出．

(20)

2 試驗及方法驗證

開展了非飽和土一維垂直滲流試驗及土水特征曲線試驗，基于試驗數(shù)據(jù)分別采用本方法及VG模型求解的方法獲得了土樣的滲透系數(shù)與體積含水量關(guān)系曲線，將兩者進行對比，驗證了該方法的可靠性．

2.1 一維垂直滲流試驗

2.1.1 試驗裝置

模擬非飽和土一維垂直入滲的試驗裝置由供水系統(tǒng)、入滲土柱及采集系統(tǒng)組成，如圖2所示．

圖2 一維垂直滲流試驗原理圖

供水系統(tǒng)主要是由潛水泵將水箱內(nèi)的蒸餾水抽入到高位水箱中，再通過溢流法將水緩慢輸入土樣中．入滲土柱試樣裝置內(nèi)徑300 mm，高600 mm，分別在50 mm、250 mm、350 mm、450 mm、550 mm處安置了1個體積含水量傳感器(試驗前已率定)．為了使水能夠均勻入滲，在土樣頂部與底部都裝有反濾層，反濾層與土樣之間用濾紙隔開，試樣底部的天平能夠?qū)崟r稱量流出水的質(zhì)量．

2.1.2 試驗土樣

試驗所用粘土過2 mm篩后，通過控制初始含水量與干密度獲得土樣的質(zhì)量，將稱好的土等分為3層，從下至上逐層均勻壓實，裝好的土樣如圖3所示．通過常規(guī)物性試驗及常水頭試驗測得的指標見表1．

圖3 裝好的人滲土柱

表1 土常規(guī)物性指標

2.1.3 試驗方案

在模擬恒定水頭下非飽和土一維垂直入滲過程時，考慮到不同的水頭高度只會影響土樣中體積含水量的變化速率，為了避免高水頭條件下體積含水量變化過快而捕捉不到詳細的變化趨勢，將水頭高度設(shè)置為2.5 cm，控制數(shù)據(jù)采集間隔Δt為1 min，待試驗數(shù)據(jù)在6 h內(nèi)不再變化時，說明入滲土柱已經(jīng)達到飽和．

整個試驗過程持續(xù)時間為1 200 min，試驗測得5個位置處體積含水量隨時間的變化曲線如圖4所示．

圖4 不同深度體積含水量隨時間變化關(guān)系曲線

從圖4中可以分析得到以下結(jié)果：①深度越淺即離水土接觸面越近的位置，含水量傳感器的讀數(shù)越先開始變化．②隨著深度的增加，傳感器的讀數(shù)維持初始狀態(tài)的時間越長．這是由于水分在土柱中的運移是從上部到下部的一個連續(xù)的由非飽和接近飽和的過程．土柱越往下，由非飽和狀態(tài)到飽和狀態(tài)所需要的時間越長．③隨著時間不斷延長，整個土柱的含水量都逐漸增加，達到穩(wěn)定后整個試驗土柱的含水量基本上一致，并且接近飽和．④5個位置的試樣在試驗達到穩(wěn)定狀態(tài)時的含水量都接近飽和，但是相比飽和含水量還有些許差距，這是由于試樣并不是固、液兩相，土體中存在沒有隨水流排出的氣體．⑤隨著深度的增加，不同位置飽和狀態(tài)的體積含水量是逐漸減小的，但是減小的程度不是很明顯．主要原因是土樣制備是從下到上逐層壓密的，從上到下土體的孔隙逐漸減小，飽和體積含水量不斷減小，這是符合實際的．

2.2 土水特征曲線試驗

為了與一維垂直滲流試驗所用土樣保持一致，控制土水特征曲線試驗所用土樣干密度為1.5 g/cm3，使用LAB523型壓力膜儀(Tempe儀)分多級壓力對3組平行試樣進行試驗，多級壓力下測得的3組含水量求平均值，試驗結(jié)果如圖5所示．

圖5 土水特征曲線

由于VG模型是目前采用最為普遍的土水特征曲線模型，故用Matlab中非線性擬合函數(shù)lsqcurvefit對土水特征曲線數(shù)據(jù)進行擬合，得到VG模型的相關(guān)參數(shù)，擬合曲線相關(guān)系數(shù)R2為0.930 37，擬合曲線如圖5所示，VG模型參數(shù)見表2．

表2 VG模型相關(guān)參數(shù)

2.3 非飽和土滲透系數(shù)的求解

1)基于差分法及試驗聯(lián)合求解：

圖6 迭代計算流程圖

2)VG模型求解

利用土水特征曲線擬合得到的VG模型參數(shù)即可求解出非飽和土滲透系數(shù)Kvg，VG模型表達式如下：

(21)

(22)

Kvg=KsΦ0.5[1-(1-Φ1/m)m]2

(23)

式中，θ為計算時段土壤體積含水量，h為基質(zhì)勢，α，m，n為經(jīng)驗擬合系數(shù)，θr為殘余體積含水量，θs為飽和體積含水量，Ks為飽和滲透系數(shù)，Kvg為相對應(yīng)θ的滲透系數(shù)．

將差分格式求出的計算值Kt與VG模型所求值Kvg進行對比分析，如圖7所示，兩種數(shù)據(jù)變化的趨勢基本相同．

圖7 計算值與VG模型求解值對比圖

2.4 誤差分析

對兩種方法求出的非飽和土滲透系數(shù)進行趨勢對比后還需要分析兩者之間的誤差．

由于非飽和土滲透系數(shù)值的變化范圍過大，跨越6個數(shù)量級，一般的均方差無法較好的判定二者之間的差異情況．為了定量分析兩種方式所求值的誤差，定義差分格式計算值與VG模型計算值比值平均值[10]為GMER．

(24)

式中，n為分析樣本數(shù)，Kt為差分格式計算值，Kvg為VG模型模擬值．按照式(24)對差分格式計算值進行誤差分析，評判標準為GMER越接近于1.0，誤差越小，擬合的精度越高．

引用王成華等[11]對幾種經(jīng)驗公式模型的誤差分析數(shù)據(jù)與差分格式的誤差分析作對比，見表3．

表3 8種經(jīng)驗公式模型GMER值

由表3可知，Vereecken、Rawls(2)、Wosten(1)、Wosten(2)模型的GMER值均小于1.0，Campbell、GHM、Rawls(1)模型的GMER值均大于1.0，Wosten(1)模型計算得到的GMER值0.94與1.0最為接近，故其擬合誤差最?。疚腉MER值3.84，相比于上述3種GMER值大于1.0的模型，其擬合誤差最?。畬Ρ劝l(fā)現(xiàn)，通過差分格式求解非飽和土滲透系數(shù)的方法具有一定的可靠性．

由于差分格式求出的是數(shù)據(jù)點，不能完整的表示滲透系數(shù)值，因此對差分格式所求值進行擬合得到了相關(guān)函數(shù)曲線．?dāng)M合函數(shù)表達式為：

lg(K)=a+b*θ

(25)

式中，a、b均為參數(shù)，相關(guān)系數(shù)R2為0.999 2．

圖8 滲透系數(shù)與體積含水量的關(guān)系曲線

由圖8表明a，b參數(shù)均與土的顆粒組成與內(nèi)部結(jié)構(gòu)有關(guān)，同一體積含水量狀態(tài)下，a參數(shù)越大，滲透系數(shù)越大，b參數(shù)與滲透系數(shù)的增長速率有關(guān)，b參數(shù)越大，滲透系數(shù)增長越快．

3 結(jié) 論

本文提出了基于一維垂直滲流控制方程差分求解及非飽和滲透性試驗聯(lián)合確定非飽和土滲透系數(shù)的新方法．主要得到了以下結(jié)論：

1)推導(dǎo)了一維垂直滲流控制方程的差分格式，該差分格式只需由試驗測得相關(guān)參數(shù)就可以通過Matlab編程求解出非飽和土滲透系數(shù)．

2)試驗粘土滲透系數(shù)的對數(shù)lg(K)與體積含水量θ的分布呈直線變化規(guī)律．

3)用差分格式與VG模型兩種方法求解的非飽和土滲透系數(shù)較為一致，用GMER算出的誤差較小，認為本文提出的求解非飽和土滲透系數(shù)的新方法具有一定的可靠性，對非飽和土體滲透系數(shù)的預(yù)測具有一定的借鑒意義．