陳遠(yuǎn)強(qiáng), 鄭 宏, 屈 新

(1. 湖南科技大學(xué) 土木工程學(xué)院, 湖南 湘潭 411201;2. 北京工業(yè)大學(xué) 城市與工程安全減災(zāi)教育部重點(diǎn)實(shí)驗(yàn)室, 北京 100124;3. 安陽(yáng)工學(xué)院 土木與建筑工程學(xué)院, 河南 安陽(yáng) 455000)

0 引 言

降雨入滲引起的水力環(huán)境變化,是邊坡失穩(wěn)的重要誘因之一．精確獲取降雨過(guò)程中邊坡水分運(yùn)移規(guī)律是準(zhǔn)確評(píng)價(jià)邊坡穩(wěn)定性的前提．大量專家學(xué)者對(duì)邊坡降雨入滲過(guò)程中的穩(wěn)定性進(jìn)行了深入研究:姚海林等[1]對(duì)降雨情況下影響非飽和膨脹土邊坡穩(wěn)定性的參數(shù)進(jìn)行了研究;董建軍等[2]研究了降雨入滲過(guò)程中非飽和滲流場(chǎng)與應(yīng)力場(chǎng)耦合作用下的邊坡穩(wěn)定性;Xiong等[3]對(duì)降雨和庫(kù)水共同作用下三峽庫(kù)區(qū)某邊坡的穩(wěn)定性進(jìn)行了研究．然而, 上述關(guān)于降雨入滲的研究, 均未考慮雨水沿地表的運(yùn)動(dòng)過(guò)程, 與真實(shí)滲流場(chǎng)往往存在一定的差異[4-5]．

為反映地表徑流對(duì)邊坡滲流場(chǎng)的影響,就必然需要建立邊坡降雨條件下徑流與滲流的耦合求解模型．目前,降雨坡面徑流與坡體滲流的各種耦合求解模型中,在坡面土體飽和之前,對(duì)坡面入滲邊界的處理方式基本一致,即將其視作流量邊界,流量大小與降雨強(qiáng)度相等;當(dāng)坡面土體飽和之后,徑流產(chǎn)生,降雨入滲邊界即由流量邊界轉(zhuǎn)換為壓力水頭邊界,壓力水頭在數(shù)值上與徑流水深相等．對(duì)于徑流產(chǎn)生之后,邊坡徑流場(chǎng)與滲流場(chǎng)之間的耦合求解大體可分為兩種:第一種是以坡面流量交換、水深相等為依據(jù),交替迭代依次求解坡面徑流和坡體滲流,直至滿足迭代收斂條件．此方法交替求解坡面徑流與坡體滲流,迭代計(jì)算量大,且難以確保徑流與滲流間流量交換完全相等,但由于其理論成熟,易于編程實(shí)現(xiàn),仍得到廣泛應(yīng)用,如Morita和Yen[6]基于地表二維擴(kuò)散波方程和地下三維非飽和滲流方程,建立了相應(yīng)的耦合求解模型;Panday等[7]發(fā)展了一種分布式地表水、地下水耦合模型,考慮了地形、植被和建筑物的影響;張培文等[8]將降雨條件下坡面徑流和降雨入滲的模擬互為條件,采用交替迭代的耦合分析方法較好地模擬了降雨坡面徑流;Erduran[9]構(gòu)建了一個(gè)綜合數(shù)值模型,模擬了植被覆蓋下地表徑流與飽和地下水流之間的相互作用;朱磊等[10]基于地表徑流物理機(jī)制和非飽和土壤水分運(yùn)動(dòng)理論,分別采用雙層節(jié)點(diǎn)耦合法和整合離散方程的整體法,建立了降雨情況下地表徑流與地下滲流的耦合模型;李根等[5]將地表徑流模型與土壤水流模型進(jìn)行耦合分析,計(jì)算了土坡降雨入滲,研究表明地表徑流對(duì)土體入滲有較大提高;李尚輝等[11]建立了降雨邊坡徑流與大孔隙流的耦合模型,并基于有限元軟件COMSOL實(shí)現(xiàn)了其數(shù)值求解．第二種求解方法是將降雨入滲面作為徑流場(chǎng)與滲流場(chǎng)的內(nèi)部域,實(shí)現(xiàn)坡面徑流控制方程與坡體非飽和滲流控制方程的統(tǒng)一聯(lián)立求解,該方法避免了求解過(guò)程中的流量交換,不需要進(jìn)行交替迭代求解,在保證計(jì)算精度的同時(shí)極大地提高了計(jì)算效率．Kollet等[12]提出了一種新的耦合模型,將地表水系統(tǒng)視作地下水系統(tǒng)的上部邊界,不需要考慮流量交換;童富果等[13]建立了降雨入滲與坡面徑流直接耦合計(jì)算模型與求解方法,避免了兩者間的迭代計(jì)算與流量交換,提高了計(jì)算效率和計(jì)算精度;Weill等[14]將擴(kuò)散波方程作適當(dāng)變換,使其在數(shù)學(xué)上與飽和-非飽和滲流Richards方程具有相同形式,并在多孔介質(zhì)表面引入徑流層概念,實(shí)現(xiàn)滲流區(qū)域和虛擬徑流層的統(tǒng)一求解,提高了計(jì)算效率;Tian等[15-16]在童富果等的研究基礎(chǔ)上進(jìn)一步深入,并將耦合模型拓展至二維坡面徑流與三維非飽和滲流的耦合;王樂(lè)等[17]建立了邊坡滲流與坡面徑流聯(lián)合求解的三維模型,并對(duì)產(chǎn)流后的入滲邊界流量進(jìn)行了修正;Liu等[18]將水氣兩相流融入地表徑流與地下水滲流的耦合模型,實(shí)現(xiàn)了邊坡降雨-入滲-產(chǎn)流的數(shù)值模擬．

對(duì)于邊坡徑流與滲流耦合模型的求解,通常采用的數(shù)值方法包括有限差分法[6-7,12]、有限單元法[8,11,13-18]和有限體積法[5,9]等．然而,上述數(shù)值方法在遇到復(fù)雜求解域時(shí)存在網(wǎng)格生成困難的問(wèn)題,而石根華博士[19]提出的數(shù)值流形法(NMM)以其獨(dú)特的覆蓋系統(tǒng),對(duì)復(fù)雜邊界具有極強(qiáng)的適應(yīng)性,目前已被廣泛應(yīng)用于滲流計(jì)算領(lǐng)域[20-22]．本文基于前人的研究基礎(chǔ),對(duì)邊坡徑流與滲流的耦合模型進(jìn)一步優(yōu)化,采用壓力水頭形式的Richards方程描述地下水非飽和滲流,用運(yùn)動(dòng)波方程描述坡面徑流,建立了耦合模型的變分提法,并基于NMM推導(dǎo)出離散求解格式,通過(guò)程序編制實(shí)現(xiàn)了邊坡降雨-入滲-產(chǎn)流的全過(guò)程數(shù)值模擬,數(shù)值算例驗(yàn)證了算法及程序的有效性與正確性．

1 控制方程及其定解條件

邊坡降雨-入滲-產(chǎn)流涉及三個(gè)方面的研究?jī)?nèi)容:降雨入滲、壤中流和坡面徑流,是一個(gè)“三水”轉(zhuǎn)換問(wèn)題,其示意圖如圖1所示．降雨入滲到坡體之中,形成壤中流,其運(yùn)移規(guī)律由飽和-非飽和滲流Richards方程[23-25]描述;當(dāng)坡表土體達(dá)到飽和狀態(tài)或者滲透能力小于降雨強(qiáng)度時(shí),坡面徑流產(chǎn)生,其運(yùn)移規(guī)律由坡面徑流Saint-Venant方程[6,9,26]描述．鑒于Saint-Venant方程求解的復(fù)雜性, 目前多采用其簡(jiǎn)化形式: 運(yùn)動(dòng)波或擴(kuò)散波方程[27-29]．因此, “三水”問(wèn)題的求解即可轉(zhuǎn)換為Richards方程與Saint-Venant方程或者其簡(jiǎn)化形式的聯(lián)立求解．

1.1 坡面徑流控制方程及其定解條件

坡面徑流由Saint-Venant方程描述,具體包含連續(xù)方程和動(dòng)量方程．對(duì)連續(xù)方程,考慮坡度影響,可以表述為

(1)

式中,h為坡面徑流水深;q為單寬流量;I=I(t)為降雨強(qiáng)度;f為地表入滲率;θ為坡面傾角;t為時(shí)間;l為坡頂沿坡面的長(zhǎng)度．

當(dāng)動(dòng)量方程中只保留坡面比降So與摩阻比降Sf時(shí),即可得運(yùn)動(dòng)波方程的水力學(xué)基礎(chǔ)[29]:

So-Sf=0．

(2)

一般地,若考慮水流為紊流狀態(tài),則摩阻比降Sf可表示為

(3)

式中,κ為無(wú)量綱的摩擦因子;m為一實(shí)數(shù)．因此,由式(3)可推導(dǎo)出坡面水流的速度公式為

(4)

通常,對(duì)式(4)中參數(shù)κ和m的取值有兩種方式:Manning公式或Chézy公式．若考慮為Manning公式,則有κ=1/nm,m=2/3,其中nm為Manning粗糙度系數(shù);若考慮為Chézy公式,則有κ=Cc,m=1/2,其中Cc為Chézy系數(shù)．當(dāng)然,參數(shù)κ和m也可以根據(jù)試驗(yàn)得到[30]．

由坡面單寬流量q=uh,并考慮為Manning公式,則一維坡面徑流的運(yùn)動(dòng)波模型可用如下方程組描述:

(5)

方程組(5)的定解條件包括初始條件和邊界條件:

1) 初始條件:若將降雨的開(kāi)始時(shí)刻作為時(shí)間起點(diǎn),則此時(shí)坡面各處均無(wú)徑流產(chǎn)生,即

(6)

2) 邊界條件:坡頂處水深和流量均為0,即

(7)

1.2 飽和-非飽和滲流控制方程及其定解條件

水在坡體中的運(yùn)移規(guī)律由飽和-非飽和滲流Richards方程描述,本文采用其壓力水頭格式(h-form),其二維表達(dá)式如下:

(8)

式中,h為壓力水頭;C(h)=?θ/?h為容水度,其中θ為體積含水量;?為梯度算子;k(h)為土體的滲透系數(shù),且k(h)=kskr(h),其中,ks為飽和滲透系數(shù),kr為相對(duì)滲透系數(shù)．需要注意的是,體積含水量θ和相對(duì)滲透系數(shù)kr均可表示為壓力水頭h的函數(shù),目前已建立了多種通用模型來(lái)描述相應(yīng)的函數(shù)關(guān)系,如Genuchten-Mualem模型[31-32]、Gardner模型[33]和Brooks-Corey模型[34]等．

為確定計(jì)算域Ω內(nèi)的滲流場(chǎng),尚需給定方程(8)的初始條件及邊界條件:

1) 初始條件:

h(x,y,0)=h0(x,y,t0), inΩ．

(9)

2) 邊界條件:

(a) 壓力水頭邊界條件Γh

(10a)

(b) 流量邊界條件Γq

(10b)

(c) 材料界面條件Γm

(10c)

2 控制方程離散及耦合模型建立

2.1 NMM簡(jiǎn)介

NMM由石根華博士于20世紀(jì)90年代初首次提出,以拓?fù)淞餍魏臀⒎至餍螢榛A(chǔ),采用兩套獨(dú)立的覆蓋系統(tǒng)(即數(shù)學(xué)覆蓋和物理覆蓋),能夠?qū)崿F(xiàn)對(duì)連續(xù)和非連續(xù)問(wèn)題的統(tǒng)一求解．?dāng)?shù)學(xué)覆蓋可視為一系列數(shù)學(xué)片的集合,數(shù)學(xué)片可以相互重疊且形狀任意,無(wú)需與求解區(qū)域的物理邊界重合,但要求為單連通區(qū)域,且其集合要完全覆蓋整個(gè)求解區(qū)域．用各種物理邊界(求解域邊界和各種不連續(xù)面)依次切割數(shù)學(xué)覆蓋中的各個(gè)數(shù)學(xué)片,并拋棄位于求解域之外的部分,就得到相應(yīng)的物理片,所有物理片的集合就組成了物理覆蓋,物理覆蓋就與求解域邊界精確匹配．而流形單元?jiǎng)t是相應(yīng)物理片之間進(jìn)一步切割形成的互不重疊的部分．由上可知,NMM的前處理極為靈活,對(duì)復(fù)雜邊界問(wèn)題和不連續(xù)性問(wèn)題具有先天優(yōu)勢(shì)．

本文以圖2為例來(lái)展示NMM的覆蓋系統(tǒng)．圖中,材料界面將求解域Ω分成兩個(gè)子域Ω1和Ω2,兩個(gè)矩形數(shù)學(xué)片MP1和MP2被用來(lái)覆蓋整個(gè)求解域．?dāng)?shù)學(xué)片MP1經(jīng)求解域的物理邊界切割,得到兩個(gè)物理片:PP1-1和PP1-2;同樣,數(shù)學(xué)片MP2經(jīng)求解域的物理邊界切割也得到兩個(gè)物理片:PP2-1和PP2-2．四個(gè)物理片之間進(jìn)一步切割,最終生成6個(gè)流形單元:E1～E6(E1源自PP1-1,E2源自PP1-1和PP2-1,E3源自PP2-1,E4源自PP1-2,E5源自PP1-2和PP2-2,E6源自PP2-2)．

(a) 求解區(qū)域Ω (b) 數(shù)學(xué)覆蓋 (c) 數(shù)學(xué)覆蓋蓋住Ω (a) Solution domain Ω (b) The mathematical cover (c) The mathematical cover on Ω

(d) MP1與Ω相互切割 (e) MP2與Ω相互切割 (d) MP1 and Ω cutting each other (e) MP2 and Ω cutting each other

(f) 流形單元生成(f) Generation of manifold elements

基于上述理論,定義在流形單元上的壓力水頭近似函數(shù)可表示為

(11)

式中,r=(x,y)代表位置向量;Np表示覆蓋該流形單元的物理片個(gè)數(shù),若數(shù)學(xué)覆蓋由三角形網(wǎng)格形成,則每個(gè)流形單元是3個(gè)物理片的交集,即Np=3,同樣,若數(shù)學(xué)覆蓋由四邊形網(wǎng)格構(gòu)成,則每個(gè)流形單元是4個(gè)物理片的交集,有Np=4;wi(r)表示第i個(gè)物理片上的單位分解函數(shù),源自生成該物理片的數(shù)學(xué)片上的權(quán)函數(shù);hi(r)表示定義于第i個(gè)物理片上的局部近似函數(shù),可表達(dá)為

hi(r)=pT(r)di,

(12)

式中,di表示定義在第i個(gè)物理片上的廣義自由度向量;p(r)為完全多項(xiàng)式基函數(shù),有零階、一階和二階等多種形式,本文采用零階基函數(shù),則其數(shù)學(xué)表達(dá)式為p(r)=1T．

考慮采用零階基函數(shù),則壓力水頭近似函數(shù)可重新表示為

(13)

式中,N=[N1,…,NNp]為形函數(shù)矩陣,其中,Ni=wi;d為廣義自由度向量,d=[d1,…,dNp]T．

2.2 控制方程離散

對(duì)坡面徑流的連續(xù)方程,考慮鏈?zhǔn)角髮?dǎo)法則:

(14)

(15)

本文采用Galerkin加權(quán)余量法對(duì)相應(yīng)控制方程進(jìn)行離散,則構(gòu)造式(15)的加權(quán)余量格式為

(16)

式中,g(x)表示測(cè)試函數(shù);Γs代表降雨徑流面．

對(duì)飽和-非飽和滲流控制方程式(8),構(gòu)造其加權(quán)余量格式為

(17)

式中,G(x,y)表示測(cè)試函數(shù)．

對(duì)式(17)中的右端項(xiàng),可通過(guò)分部積分降低階次:

(18)

式中,?Ω為求解區(qū)域Ω的外邊界,?Ω=Γh+Γq;ni為邊界外法線向量的方向余弦．

將式(18)代入式(17),結(jié)合流量邊界條件式(10b),并采用罰函數(shù)法施加本質(zhì)邊界條件和界面連續(xù)性條件,則式(17)可進(jìn)一步表示為

(19)

式中,kp為罰因子．

用式(13)中的近似場(chǎng)函數(shù)分別逼近式(16)與式(19)中的真實(shí)場(chǎng)函數(shù),并采用Galerkin加權(quán)余量法,則可得到坡面徑流和飽和-非飽和滲流的總體方程分別為

(20a)

(20b)

(21a)

(21b)

(21c)

(21d)

(21e)

(21f)

式中,B=LN,L=[?/?x,?/?y]T;k為滲透系數(shù)矩陣．

2.3 耦合模型建立

當(dāng)對(duì)坡面徑流與坡體滲流進(jìn)行單獨(dú)分析時(shí),對(duì)坡面徑流而言,需要給定坡面處的下滲分布情況;對(duì)于坡體的飽和-非飽和滲流而言,需要給定坡面處的入滲分布情況．因此,耦合模型需要解決二者間流量交換的問(wèn)題．從現(xiàn)實(shí)角度來(lái)看,坡面處的下滲分布和入滲分布是一致的,如果將坡面邊界作為內(nèi)部域,就無(wú)須考慮二者之間的流量交換．

當(dāng)坡體表面未出現(xiàn)徑流時(shí),降雨全部滲進(jìn)邊坡,此時(shí)僅需考慮飽和-非飽和滲流總體控制方程式(20b);當(dāng)坡面出現(xiàn)徑流之后,坡面處的水深應(yīng)滿足式(20a),整個(gè)坡體的壓力水頭應(yīng)滿足式(20b),此時(shí),坡面處的壓力水頭與徑流水深在數(shù)值上應(yīng)近似相等,即式(20a)與式(20b)求解變量相同．又考慮坡面徑流與坡體滲流共有坡面邊界,為敘述方便,將坡面邊界記為Γgs,假定坡面土體的入滲率為f,則式(20a)、(20b)中單元上的等效節(jié)點(diǎn)流量列陣可分別重寫(xiě)為

(22a)

(22b)

式中,F0為飽和-非飽和滲流問(wèn)題中除降雨坡面邊界外的其他邊界形成的等效節(jié)點(diǎn)流量列陣．

將式(20a)與(20b)相加,并考慮式(22a)與(22b),即可得到降雨條件下坡面徑流與坡體滲流全耦合模型的總體控制方程:

(23)

2.4 非線性迭代求解

對(duì)耦合模型的總體控制方程式(23),采用差分法對(duì)時(shí)間域進(jìn)行離散化處理,有

(24)

將式(24)代入式(23),可得耦合模型的流形元迭代求解格式:

(25)

式中,下標(biāo)n表示時(shí)間步數(shù);Δt為時(shí)間步長(zhǎng);θ為權(quán)重參數(shù),0≤θ≤1,θ的取值不同,對(duì)應(yīng)著不同的時(shí)間差分格式,也將直接影響解的精度和穩(wěn)定性．研究表明,當(dāng)θ=1,即對(duì)時(shí)間域的離散采用向后差分公式時(shí),式(25)在求解理論上是無(wú)條件穩(wěn)定的,因此本文亦采用向后差分公式．

由于總體控制方程式(25)具有強(qiáng)烈的非線性,需要采用非線性迭代方法進(jìn)行求解．本文采用Picard迭代法進(jìn)行迭代求解,取θ=1,則迭代求解格式(25)可進(jìn)一步表示為

(26)

式中,m表示迭代步數(shù);hn+1,m+1表示第n+1個(gè)時(shí)間步中第m+1個(gè)迭代步的壓力水頭列陣,其余下標(biāo)含義與此類似．

當(dāng)絕對(duì)誤差達(dá)到給定的容差ε時(shí),迭代終止,即

‖hn+1,m+1-hn+1,m‖≤ε,

(27)

式中,‖·‖表示∞-范數(shù)．當(dāng)收斂條件滿足時(shí),令hn+1=hn+1,m+1,然后進(jìn)入下一時(shí)間步的計(jì)算．

3 數(shù) 值 算 例

3.1 Abdul-Gillham試驗(yàn)驗(yàn)證

Abdul-Gillham[35]試驗(yàn)經(jīng)常被用來(lái)檢驗(yàn)坡面徑流與飽和-非飽和滲流耦合求解模型的性能．試驗(yàn)在一個(gè)長(zhǎng)為140 cm,高為120 cm,寬為8 cm的玻璃箱內(nèi)進(jìn)行．箱內(nèi)鋪設(shè)一定厚度的中細(xì)砂,使其形成坡角為12°的斜坡,坡腳高度為76 cm,模型幾何尺寸如圖3所示．初始時(shí)刻,水位面與坡腳高度一致,位于76 cm處,坡體底部和兩側(cè)均為不透水邊界．坡面上方設(shè)置降雨模擬器,模擬的降雨強(qiáng)度為1.2×10-5m/s,持續(xù)時(shí)間為20 min．坡腳處安裝監(jiān)測(cè)系統(tǒng),用于記錄坡腳處的出流情況,監(jiān)測(cè)總時(shí)長(zhǎng)為25 min．

圖3 Abdul-Gillham試驗(yàn)計(jì)算幾何模型Fig. 3 The geometric model for the Abdul-Gillham system

計(jì)算時(shí),土-水特征曲線和滲透性函數(shù)分別由van Genuchten模型[31]和Mualem模型[32]來(lái)描述:

(28)

kr=Θ1/2[1-(1-Θ1/c)c]2,

(29)

式中,θs為飽和體積含水率;θr為殘余體積含水率;Θ為有效飽和度,具體表示為Θ=(θ-θr)/(θs-θr);a,b,c均為與土體性質(zhì)有關(guān)的參數(shù),其中c=1-1/b．本算例中,參數(shù)a=2.3 m-1,b=5.5．土體的飽和滲透系數(shù)為3.5×10-5m/s,飽和體積含水率取為0.34,殘余含水率取為1×10-6．坡面的Manning粗糙度系數(shù)為0.185 m-1/3·s．

分別采用三角形和四邊形網(wǎng)格形成的數(shù)學(xué)覆蓋去覆蓋整個(gè)求解區(qū)域,由于降雨時(shí)坡面處水力梯度變化劇烈,為提高計(jì)算精確度,本文對(duì)坡面處網(wǎng)格進(jìn)行了加密,其幾何模型和數(shù)學(xué)覆蓋如圖4所示．三角形和四邊形網(wǎng)格形成的數(shù)學(xué)覆蓋生成的物理片個(gè)數(shù)分別為1 196和1 132, 流形單元數(shù)目分別為2 256和1 204．?dāng)?shù)值模擬得到的坡腳出流情況如圖5所示,同時(shí)圖中也給出了Abdul、Gillham[35]的試驗(yàn)結(jié)果與VanderKwaak[36]、Kollet等[12]、Tian等[16]的模擬結(jié)果．由圖5可知,本文采用三角網(wǎng)格和四邊形網(wǎng)格形成數(shù)學(xué)覆蓋的坡腳出流過(guò)程線幾乎完全重合,并與VanderKwaak、Kollet等和Tian等的模擬結(jié)果基本一致,均與試驗(yàn)結(jié)果吻合良好,從而驗(yàn)證了本文所建立的耦合求解模型是正確可靠的．此外,由圖5可以看出,數(shù)值模擬結(jié)果的產(chǎn)流時(shí)間均要早于試驗(yàn)結(jié)果,分析其原因?yàn)閿?shù)值模型中均沒(méi)有考慮土壤中空氣的壓縮性．

(a) 三角形網(wǎng)格形成的數(shù)學(xué)覆蓋 (b) 四邊形網(wǎng)格形成的數(shù)學(xué)覆蓋 (a) The mathematical cover generated by triangular grids (b) The mathematical cover generated by quadrilateral grids

圖5 坡腳出流過(guò)程對(duì)比Fig. 5 Discharges at the base of the slope

3.2 Smith試驗(yàn)驗(yàn)證

Smith等[30]曾對(duì)降雨產(chǎn)流進(jìn)行過(guò)試驗(yàn)研究．試驗(yàn)土體長(zhǎng)為12.2 m,厚度為1.22 m,寬為5.1 cm,土體按照密度不同大致分為四層,各層厚度分別為1.27 cm,6.35 cm,22.86 cm和91.52 cm,土槽底部坡度為0.01,其幾何示意圖如圖6所示．其中,頂層土體和第三層土體密度相同,飽和滲透系數(shù)為0.254 cm/min,第一層土體和底層土體飽和滲透系數(shù)分別為0.394 cm/min和0.186 cm/min．土-水特征曲線和滲透系數(shù)函數(shù)由Smith等經(jīng)試驗(yàn)測(cè)定,并采用Brooks-Corey模型[34]描述,Singh等[37]對(duì)Simth等的試驗(yàn)數(shù)據(jù)進(jìn)行了更為詳細(xì)地分析,得到了更為完善的模型,其具體參數(shù)見(jiàn)文獻(xiàn)[37]．模擬降雨強(qiáng)度為0.42 cm/min,歷時(shí)15 min,坡面徑流速度計(jì)算公式采用Smith等得到的經(jīng)驗(yàn)公式u=CSoh2,其中參數(shù)C=7.87×105cm-1·s-1．

圖6 Smith試驗(yàn)計(jì)算幾何模型Fig. 6 The geometric model for the Smith system

模擬時(shí),分別采用三角形和四邊形網(wǎng)格形成的數(shù)學(xué)覆蓋去覆蓋整個(gè)求解區(qū)域,兩種覆蓋方法生成的物理片數(shù)目均為3 233,生成的流形單元數(shù)目分別為5 600和2 800．計(jì)算得到的坡腳出流情況對(duì)比如圖7所示,由圖可以看出,本文兩種數(shù)學(xué)覆蓋的計(jì)算結(jié)果及Smith等[30]、Morita等[6]模擬結(jié)果與實(shí)測(cè)出流情況基本一致,再次驗(yàn)證了本文耦合模型的正確性．圖8給出了x=5.6 m剖面上的土體飽和度隨時(shí)間變化過(guò)程,模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)存在一定差異,Smith等和Morita等的模擬結(jié)果也存在類似問(wèn)題,分析其原因可能是初始時(shí)刻土壤含水率的實(shí)測(cè)數(shù)據(jù)較少,并不能真實(shí)反映實(shí)際含水情況,導(dǎo)致建立的土-水特征曲線和滲透系數(shù)函數(shù)模型并不能準(zhǔn)確表征土體的物理特性．但從飽和度隨時(shí)間的演化過(guò)程來(lái)看,其符合降雨入滲規(guī)律．因此,綜上所述,可認(rèn)為本文計(jì)算結(jié)果是客觀合理的．

圖7 坡腳出流過(guò)程對(duì)比

3.3 均質(zhì)邊坡降雨入滲

本小節(jié)對(duì)均質(zhì)邊坡的降雨入滲展開(kāi)模擬．假定邊坡水平長(zhǎng)度為100 m,豎直厚度為40 m,坡角為15°,其幾何模型如圖9所示．土體的飽和滲透系數(shù)為4.332×10-4m/min,土-水特征曲線及滲透性函數(shù)如圖10所示．假定初始時(shí)刻土體的體積含水量為0.045,坡面Manning粗糙度系數(shù)為0.035 m-1/3·min．計(jì)算時(shí),分別采用表1中5種不同工況的降雨強(qiáng)度進(jìn)行模擬分析,以探究降雨強(qiáng)度對(duì)坡面產(chǎn)流的影響,降雨歷時(shí)10 h．本小節(jié)僅采用四邊形網(wǎng)格形成的數(shù)學(xué)覆蓋去覆蓋整個(gè)計(jì)算區(qū)域,共生成1 458個(gè)物理片和1 374個(gè)流形單元．

圖9 均質(zhì)邊坡幾何模型

表1 降雨工況

圖11和12分別為5種不同降雨工況計(jì)算得到的平衡條件坡面水位線圖和坡腳出流情況．由圖可知,降雨強(qiáng)度越大,坡面產(chǎn)流時(shí)間越早,達(dá)到平衡條件時(shí)的坡面積水越深,坡腳徑流量也越大．工況1降雨強(qiáng)度最大,產(chǎn)流時(shí)間最早,約為25 min,達(dá)到平衡條件時(shí)坡腳積水深度為8.5 cm;而工況5降雨強(qiáng)度最小,降雨約79 min之后坡面才開(kāi)始產(chǎn)流,達(dá)到平衡條件時(shí)的坡腳積水深度也僅有4.3 cm．

圖11 平衡條件下的坡面水位線

圖13和14分別給出了工況1情況下降雨結(jié)束時(shí)刻坡體和坡面的壓力水頭分布．由圖可以看出,由于初始時(shí)刻邊坡的體積含水量較低,降雨只對(duì)邊坡表面產(chǎn)生影響,最大影響深度約為1 m左右,而邊坡底層土體的體積含水量基本不發(fā)生變化．

圖13 降雨結(jié)束時(shí)刻工況1邊坡壓力水頭分布

4 結(jié) 論

本文采用一維運(yùn)動(dòng)波方程描述坡面徑流,用二維壓力水頭形式的Richards方程描述土體非飽和滲流,考慮將降雨入滲面作為徑流與滲流的內(nèi)部域,推導(dǎo)出邊坡降雨徑流與滲流全耦合模型的總體控制方程．基于NMM對(duì)其進(jìn)行數(shù)值離散,建立了對(duì)應(yīng)的迭代求解格式,并通過(guò)編程實(shí)現(xiàn)了邊坡降雨-入滲-產(chǎn)流的全過(guò)程數(shù)值模擬．研究表明,本文模型及計(jì)算方法準(zhǔn)確可靠,能夠較為真實(shí)地反映邊坡的降雨入滲與產(chǎn)流過(guò)程,可為降雨型滑坡的穩(wěn)定性評(píng)價(jià)及排水治理設(shè)計(jì)提供技術(shù)參考．