李怡航

(四川大學(xué)水利水電學(xué)院，四川 成都 610065)

1 概述

降雨入滲是自然界中水循環(huán)的重要環(huán)節(jié)。降雨入滲是指降雨后水分從地表垂直向下進(jìn)入土壤，促使土壤從不飽和狀態(tài)向飽和狀態(tài)轉(zhuǎn)變的過程。水體滲流不僅可使粘土質(zhì)巖石發(fā)生泥化、軟化現(xiàn)象，還可促使膨脹土邊坡或巖質(zhì)邊坡的體積膨脹，增大自重壓力[1]，產(chǎn)生滑坡?？梢妼?duì)于滲流模型的建立是十分必要的。

在模擬入滲運(yùn)動(dòng)數(shù)值模擬方面，文獻(xiàn)[2]中采用有限差分法對(duì)其進(jìn)行求解，但此法對(duì)數(shù)值解的守恒性難以保證，在遇到邊界條件和復(fù)雜區(qū)域時(shí)，會(huì)出現(xiàn)數(shù)值振蕩和彌散現(xiàn)象。建立降雨入滲模型是分析邊坡的穩(wěn)定性較為常用的研究方法[3]，如陳守義用極限平衡分析法模擬的入滲模型[4]，于玉貞、劉金龍等利用有限元滲流分析的模型[5- 6]，而有限體積法結(jié)合了有限差分和有限元，其積分方程表示的是控制體積的通量平衡，因此數(shù)值解的守恒性得到了保障，相比于其他方法能有效地克服數(shù)值彌散和振蕩現(xiàn)象[7]。

本文對(duì)一維Richards方程進(jìn)行有限體積離散[8]，基于MATLAB語言編寫一維入滲的程序，結(jié)合了木拉提等的中壤土入滲試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，從而來驗(yàn)證程序的可靠性。本文還分析了“降雨持時(shí)相同而降雨強(qiáng)度不同”和“降雨總量相同而降雨強(qiáng)度不同”兩方面因素對(duì)一維入滲模型的空間分布規(guī)律進(jìn)行分析。

2 數(shù)學(xué)模型及其求解

2.1 一維入滲運(yùn)動(dòng)的控制方程

本文研究的土壤水分垂直運(yùn)動(dòng)是針對(duì)單一土壤的，因此可假設(shè)其為均質(zhì)、各向同性且骨架不變形的多孔介質(zhì)。不考慮土壤中根系吸水的影響，Richards通過實(shí)驗(yàn)驗(yàn)證了非飽和土壤水運(yùn)動(dòng)符合Darcy定律[9- 10]，將其引入非飽和土壤水運(yùn)動(dòng)，結(jié)合Buckingham的能量法提出了不飽和土壤水運(yùn)動(dòng)的控制方程，土壤水分入滲運(yùn)動(dòng)用Richards方程描述為：

(1)

式中，θ—土壤體積含水率；K(θ)—土壤導(dǎo)水率；D(θ)—非飽和土壤水?dāng)U散率。

D(θ)根據(jù)其定義是導(dǎo)水率K(θ)和比水容量C(θ)的比值，即：

(2)

2.2 計(jì)算網(wǎng)格的處理

圖1為一維垂直方向上入滲模型網(wǎng)格示意圖。

圖1 入滲模型FVM示意圖

在距地面垂直距離L這段空間中分為n個(gè)節(jié)點(diǎn)，空間步長(zhǎng)選取Δz，并分別編號(hào)j=1～n，其中j=1和j=n的兩處節(jié)點(diǎn)分別為計(jì)算模型的靠近上邊界和下邊界的節(jié)點(diǎn)，且距離邊界距離Δz/2。時(shí)間步長(zhǎng)選取Δt，時(shí)間層編序號(hào)依次為k=0，1，2，…。

本文選取控制體l-r為對(duì)象進(jìn)行研究，由于點(diǎn)l、r并非模型中的節(jié)點(diǎn)，其K、D值不能直接獲得，但可通過相鄰節(jié)點(diǎn)L、M、R取近似值。

2.3 邊界條件

2.3.1上邊界條件的處理

(1)當(dāng)上邊界的土壤處于飽和狀態(tài)時(shí)，為第一類邊界條件。取控制體cv中包含節(jié)點(diǎn)j=1，且l即與上邊界重合，因此：

2.3.2下邊界條件的處理

入滲模型的下邊界條件：θ(L，t)=θ0，即視作在x=L處的體積含水率不隨時(shí)間的變化而發(fā)生變化，即該值為一恒定值。在節(jié)點(diǎn)j=n處截取控制體cv，則r與下邊界重合，即：

2.4 模型求解

式(1)的有限體積法表述為：

(3)

關(guān)于導(dǎo)水率，本文采用能適用于水流對(duì)流擴(kuò)散的Up-wind格式，即：

(4)

關(guān)于擴(kuò)散率，本文采用算術(shù)平均法，即：

(5)

因此，對(duì)于非邊界節(jié)點(diǎn)(j=2，3，…，n-1)，式(5)最終求解可得下式：

(6)

α的取值不同會(huì)產(chǎn)生不同的計(jì)算格式，其中最常用的有顯式格式(α=0)、C-N格式(α=0.5)以及全隱格式(α=1)。其中全隱格式的各項(xiàng)系數(shù)均為正數(shù)，且對(duì)時(shí)間步長(zhǎng)的選取是無條件穩(wěn)定的，本文即選取全隱格式進(jìn)行計(jì)算，即選取α=1，每一時(shí)間層的含水率分布均是通過多次迭代得到。

令r=Δt/Δz2，采用全隱格式后，式(6)可簡(jiǎn)化為：

aLθL+aMθM+aRθR=b

(7)

將上下邊界條件進(jìn)行有限體積離散，亦可得到式(7)，但其系數(shù)稍有不同。

由此，將各節(jié)點(diǎn)上含水量的計(jì)算表達(dá)式以矩陣形式表達(dá)：

(8)

其中對(duì)于每一節(jié)點(diǎn)j=2，3，…，L-1來說，aL=Aj、aM=Bj、aR=Cj。式(8)左邊的第一個(gè)矩陣為為三對(duì)角矩陣，左邊第二個(gè)矩陣與含水率有關(guān)，是未知矩陣，其中θ1、θL均是通過上下邊界條件確定。求解上式可利用Matlab工具求解，計(jì)算流程如圖2所示。

圖2 入滲模型程序流程圖

3 數(shù)學(xué)模型的驗(yàn)證

程序的有效性需要借助試驗(yàn)對(duì)模型求解的空間含水率分布進(jìn)行驗(yàn)證。下面將通過木拉提·胡塞因等[11]對(duì)中壤土(其干容重為γd=1.4×103kg/m3)進(jìn)行測(cè)試的室內(nèi)試驗(yàn)數(shù)據(jù)進(jìn)行程序驗(yàn)證。其試驗(yàn)的條件如下：

(1)初始條件：θ0=0.09(初始含水率)。

(2)邊界條件：θs=0.42(上邊界條件)θ|z=60=0.09(下邊界條件)。

(3)參數(shù)選?。?/p>

計(jì)算時(shí)，空間步長(zhǎng)選取Δz=1cm，時(shí)間步長(zhǎng)Δt=1min。程序運(yùn)算t=600，1230，1890min時(shí)各節(jié)點(diǎn)的體積含水量，最后程序運(yùn)算結(jié)果與試驗(yàn)數(shù)據(jù)如圖3所示，圖3中計(jì)算結(jié)果同試驗(yàn)結(jié)果較為吻合。

圖3 入滲模型含水率空間分布圖

在誤差結(jié)果分析中，抽取試驗(yàn)數(shù)據(jù)在對(duì)應(yīng)時(shí)刻對(duì)應(yīng)深度的含水率計(jì)算值，與試驗(yàn)數(shù)據(jù)進(jìn)行誤差分析，得到表1。

由表1分析可知，采集的計(jì)算點(diǎn)其相對(duì)誤差維持在15%以內(nèi)，且絕對(duì)誤差也控制在0.2以內(nèi)，因此可以認(rèn)為該程序運(yùn)算結(jié)果是可靠的。

4 空間分布規(guī)律與因素影響分析

在驗(yàn)證章節(jié)2中的有限體積法模型的可靠性后，利用該模型對(duì)空間含水率分布進(jìn)行進(jìn)一步分析。在分析中采用文獻(xiàn)[9]的試驗(yàn)，土壤類型為壤土，其參數(shù)有：塑性指數(shù)Ip=11.8，而干容重γd=1.365g/cm3。

初始條件：初始含水量沿垂直方向恒為0.028。

邊界條件：上邊界最高含水量達(dá)0.445。

表1 誤差分析表

D(θ)和K(θ)的關(guān)系式是經(jīng)驗(yàn)公式，根據(jù)文獻(xiàn)得到：

將編寫的上述算例計(jì)算程序在Matlab R2017a軟件上運(yùn)行。將結(jié)果進(jìn)行處理后對(duì)降雨入滲后的含水率空間分布曲線特征進(jìn)行分析，取時(shí)間步長(zhǎng)Δt=1min，空間步長(zhǎng)Δz=1cm，降雨持時(shí)tm=1600min，降雨強(qiáng)度I=0.0500cm/min為例，其計(jì)算后得到的含水率空間分布曲線如圖4所示。

圖4 含水率空間分布曲線及分段

該曲線與Green-Ampt[12- 15]入滲模型近似，濕潤(rùn)鋒可近似看作水平線。該曲線可分成四段：

Ⅰ號(hào)曲線是從z=0cm出開始，該段曲線的含水率接近或達(dá)到飽和含水率，該段曲線近乎垂直。

Ⅱ號(hào)曲線的含水率隨深度的變化已較為明顯，但是該段曲線隨深度增加，含水率值的減小幅度較小，且各點(diǎn)的切線角度已明顯不是垂直或者水平。

Ⅲ號(hào)曲線的含水率變化非常明顯，該段曲線幾乎水平，且該段曲線即為入滲模型的濕潤(rùn)鋒。

Ⅳ號(hào)曲線處于濕潤(rùn)鋒前，其含水率隨深度未發(fā)生變化，維持在初始含水率狀態(tài)。

以“降雨持時(shí)相同而降雨強(qiáng)度不同”和“降雨總量相同而降雨強(qiáng)度不同”兩方面因素對(duì)入滲模型的空間分布規(guī)律進(jìn)行分析。

(1)降雨持時(shí)相同而降雨強(qiáng)度不同

從直觀感覺上分析，降雨強(qiáng)度越大，入滲量越大，濕潤(rùn)鋒深度zf越深。該觀點(diǎn)是否正確，對(duì)程序輸入時(shí)間步長(zhǎng)Δt=1min，空間步長(zhǎng)Δz=1cm，迭代次數(shù)為tm=1600min，降雨強(qiáng)度I分別為0.1000，0.0500，0.0250，0.0057cm/min，運(yùn)算得到各降雨強(qiáng)度下最終含水率空間分布圖如圖5所示。

圖5 不同降雨強(qiáng)度下含水率空間分布曲線

圖5中四條線的分布符合“降雨強(qiáng)度越大，入滲量越大，濕潤(rùn)鋒深度zf越深”的觀點(diǎn)。但是，降雨強(qiáng)度0.1000cm/min、0.0500cm/min、0.0250cm/min之間的數(shù)量差距較大，但最終的入滲量與濕潤(rùn)鋒深度差距卻不是很大，這與降雨強(qiáng)度為0.0057cm/min的結(jié)果曲線卻有極大差異。而且，0.0057cm/min強(qiáng)度的結(jié)果曲線含水率并未達(dá)到飽和狀態(tài)。分析原因如下：

在降雨初期入滲能力是大于降雨強(qiáng)度的，因此在初期入滲強(qiáng)度等于降雨強(qiáng)度。隨著持續(xù)的降雨，土壤含水率不斷上升直至飽和，表層土壤的入滲能力卻隨著土壤含水率增加而減小。待其入滲能力等于降雨強(qiáng)度，坡面便開始產(chǎn)流(定義產(chǎn)流時(shí)刻的入滲能力為產(chǎn)流入滲能力，符號(hào)fp)，不過此刻表層土壤的含水率未必在飽和狀態(tài)，入滲能力仍有可能繼續(xù)下降。飽和后，土壤的入滲強(qiáng)度維持不變且小于降雨強(qiáng)度，因此表面能形成積水。

(2)降雨總量相同而降雨強(qiáng)度不同

取時(shí)間步長(zhǎng)Δt=1min，空間步長(zhǎng)Δz=1cm，以降雨總量為480mm設(shè)定了四個(gè)降雨事件，其降雨持時(shí)和降雨強(qiáng)度見表3。

表3 降雨總量為480mm的不同降雨事件表

按照表3中事件所列述的降雨強(qiáng)度與降雨持時(shí)帶入程序中，運(yùn)算并繪制了降雨總量為480mm的不同降雨事件下含水率的空間分布曲線，如圖6所示。

圖6 不同降雨事件的含水率空間分布曲線

如圖6所示，在降雨總量不變的情況下，隨著降雨強(qiáng)度的減小，入滲量反而增加，濕潤(rùn)鋒深度加深。當(dāng)降雨強(qiáng)度為0.0057cm/min，降雨持時(shí)為8640min時(shí)，降雨入滲已達(dá)到計(jì)算深度，但是其表層無法形成積水，上層邊界土也未達(dá)到飽和狀態(tài)，降雨量全部入滲到土體中。

降雨強(qiáng)度越大，表層土壤含水率增加速率上升，產(chǎn)流時(shí)間縮短。因此，產(chǎn)流的時(shí)間以及產(chǎn)流量明顯增加，造成這一現(xiàn)象的發(fā)生。這說明，降雨強(qiáng)度越低，形成的徑流高度越小甚至無法形成徑流。

5 總結(jié)

(1)本文推導(dǎo)了一維狀態(tài)下降雨入滲模型控制方程的有限體積格式，結(jié)合上下邊界的離散利用MATLAB程序?qū)δＰ瓦M(jìn)行求解。并通過已知文獻(xiàn)數(shù)據(jù)驗(yàn)證了入滲模型計(jì)算程序的可靠性，程序計(jì)算數(shù)據(jù)結(jié)果與文獻(xiàn)試驗(yàn)記錄結(jié)果接近，兩者比較得到的誤差均控制在15%以內(nèi)，因此程序的計(jì)算結(jié)果是可靠的，可用于初步降雨入滲的預(yù)測(cè)與分析。

(2)從降雨強(qiáng)度和空間的關(guān)系對(duì)影響入滲模型的因素進(jìn)行分析。當(dāng)降雨持時(shí)相同而降雨強(qiáng)度不同，服從“降雨強(qiáng)度越大，入滲量越大，濕潤(rùn)鋒深度zf越深”的規(guī)律，但若降雨強(qiáng)度小到一定程度則土壤各節(jié)點(diǎn)均不能達(dá)到飽和狀態(tài)；當(dāng)降雨總量相同而降雨強(qiáng)度不同，降雨強(qiáng)度越大，其產(chǎn)流時(shí)間越短，入滲量越少，濕潤(rùn)鋒深度越淺。