一維土壤水分運動模型的有限體積法求解及其影響因素分析

李怡航

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

1 概述

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

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

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

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

2.1 一維入滲運動的控制方程

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

(1)

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

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

(2)

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

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

圖1 入滲模型FVM示意圖

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

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

2.3 邊界條件

2.3.1上邊界條件的處理

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

2.3.2下邊界條件的處理

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

2.4 模型求解

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

(3)

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

(4)

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

(5)

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

(6)

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

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

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

(7)

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

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

(8)

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

圖2 入滲模型程序流程圖

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

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

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

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

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

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

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

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

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

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

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

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

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

表1 誤差分析表

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

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

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

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

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

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

Ⅲ號曲線的含水率變化非常明顯，該段曲線幾乎水平，且該段曲線即為入滲模型的濕潤鋒。

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

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

(1)降雨持時相同而降雨強度不同

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

圖5 不同降雨強度下含水率空間分布曲線

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

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

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

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

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

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

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

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

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

5 總結(jié)

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

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