陳遠強, 鄭 宏, 屈 新

(1. 湖南科技大學 土木工程學院, 湖南 湘潭 411201;2. 北京工業(yè)大學 城市與工程安全減災教育部重點實驗室, 北京 100124;3. 安陽工學院 土木與建筑工程學院, 河南 安陽 455000)

0 引 言

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

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

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

1 控制方程及其定解條件

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

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

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

(1)

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

當動量方程中只保留坡面比降So與摩阻比降Sf時,即可得運動波方程的水力學基礎[29]:

So-Sf=0．

(2)

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

(3)

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

(4)

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

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

(5)

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

1) 初始條件:若將降雨的開始時刻作為時間起點,則此時坡面各處均無徑流產生,即

(6)

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

(7)

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

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

(8)

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

為確定計算域Ω內的滲流場,尚需給定方程(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簡介

NMM由石根華博士于20世紀90年代初首次提出,以拓撲流形和微分流形為基礎,采用兩套獨立的覆蓋系統(tǒng)(即數學覆蓋和物理覆蓋),能夠實現對連續(xù)和非連續(xù)問題的統(tǒng)一求解．數學覆蓋可視為一系列數學片的集合,數學片可以相互重疊且形狀任意,無需與求解區(qū)域的物理邊界重合,但要求為單連通區(qū)域,且其集合要完全覆蓋整個求解區(qū)域．用各種物理邊界(求解域邊界和各種不連續(xù)面)依次切割數學覆蓋中的各個數學片,并拋棄位于求解域之外的部分,就得到相應的物理片,所有物理片的集合就組成了物理覆蓋,物理覆蓋就與求解域邊界精確匹配．而流形單元則是相應物理片之間進一步切割形成的互不重疊的部分．由上可知,NMM的前處理極為靈活,對復雜邊界問題和不連續(xù)性問題具有先天優(yōu)勢．

本文以圖2為例來展示NMM的覆蓋系統(tǒng)．圖中,材料界面將求解域Ω分成兩個子域Ω1和Ω2,兩個矩形數學片MP1和MP2被用來覆蓋整個求解域．數學片MP1經求解域的物理邊界切割,得到兩個物理片:PP1-1和PP1-2;同樣,數學片MP2經求解域的物理邊界切割也得到兩個物理片:PP2-1和PP2-2．四個物理片之間進一步切割,最終生成6個流形單元: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) 數學覆蓋 (c) 數學覆蓋蓋住Ω (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

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

(11)

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

hi(r)=pT(r)di,

(12)

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

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

(13)

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

2.2 控制方程離散

對坡面徑流的連續(xù)方程,考慮鏈式求導法則:

(14)

(15)

本文采用Galerkin加權余量法對相應控制方程進行離散,則構造式(15)的加權余量格式為

(16)

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

對飽和-非飽和滲流控制方程式(8),構造其加權余量格式為

(17)

式中,G(x,y)表示測試函數．

對式(17)中的右端項,可通過分部積分降低階次:

(18)

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

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

(19)

式中,kp為罰因子．

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

(20a)

(20b)

(21a)

(21b)

(21c)

(21d)

(21e)

(21f)

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

2.3 耦合模型建立

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

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

(22a)

(22b)

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

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

(23)

2.4 非線性迭代求解

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

(24)

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

(25)

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

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

(26)

式中,m表示迭代步數;hn+1,m+1表示第n+1個時間步中第m+1個迭代步的壓力水頭列陣,其余下標含義與此類似．

當絕對誤差達到給定的容差ε時,迭代終止,即

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

(27)

式中,‖·‖表示∞-范數．當收斂條件滿足時,令hn+1=hn+1,m+1,然后進入下一時間步的計算．

3 數 值 算 例

3.1 Abdul-Gillham試驗驗證

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

圖3 Abdul-Gillham試驗計算幾何模型Fig. 3 The geometric model for the Abdul-Gillham system

計算時,土-水特征曲線和滲透性函數分別由van Genuchten模型[31]和Mualem模型[32]來描述:

(28)

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

(29)

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

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

(a) 三角形網格形成的數學覆蓋 (b) 四邊形網格形成的數學覆蓋 (a) The mathematical cover generated by triangular grids (b) The mathematical cover generated by quadrilateral grids

圖5 坡腳出流過程對比Fig. 5 Discharges at the base of the slope

3.2 Smith試驗驗證

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

圖6 Smith試驗計算幾何模型Fig. 6 The geometric model for the Smith system

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

圖7 坡腳出流過程對比

3.3 均質邊坡降雨入滲

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

圖9 均質邊坡幾何模型

表1 降雨工況

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

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

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

圖13 降雨結束時刻工況1邊坡壓力水頭分布

4 結 論

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