基于数值流形法的降雨入渗与坡面径流耦合算法研究

陈远强, 郑 宏, 屈 新

(1. 湖南科技大学 土木工程学院, 湖南 湘潭 411201;2. 北京工业大学 城市与工程安全减灾教育部重点实验室, 北京 100124;3. 安阳工学院 土木与建筑工程学院, 河南 安阳 455000)

0 引 言

降雨入渗引起的水力环境变化,是边坡失稳的重要诱因之一．精确获取降雨过程中边坡水分运移规律是准确评价边坡稳定性的前提．大量专家学者对边坡降雨入渗过程中的稳定性进行了深入研究:姚海林等[1]对降雨情况下影响非饱和膨胀土边坡稳定性的参数进行了研究;董建军等[2]研究了降雨入渗过程中非饱和渗流场与应力场耦合作用下的边坡稳定性;Xiong等[3]对降雨和库水共同作用下三峡库区某边坡的稳定性进行了研究．然而, 上述关于降雨入渗的研究, 均未考虑雨水沿地表的运动过程, 与真实渗流场往往存在一定的差异[4-5]．

为反映地表径流对边坡渗流场的影响,就必然需要建立边坡降雨条件下径流与渗流的耦合求解模型．目前,降雨坡面径流与坡体渗流的各种耦合求解模型中,在坡面土体饱和之前,对坡面入渗边界的处理方式基本一致,即将其视作流量边界,流量大小与降雨强度相等;当坡面土体饱和之后,径流产生,降雨入渗边界即由流量边界转换为压力水头边界,压力水头在数值上与径流水深相等．对于径流产生之后,边坡径流场与渗流场之间的耦合求解大体可分为两种:第一种是以坡面流量交换、水深相等为依据,交替迭代依次求解坡面径流和坡体渗流,直至满足迭代收敛条件．此方法交替求解坡面径流与坡体渗流,迭代计算量大,且难以确保径流与渗流间流量交换完全相等,但由于其理论成熟,易于编程实现,仍得到广泛应用,如Morita和Yen[6]基于地表二维扩散波方程和地下三维非饱和渗流方程,建立了相应的耦合求解模型;Panday等[7]发展了一种分布式地表水、地下水耦合模型,考虑了地形、植被和建筑物的影响;张培文等[8]将降雨条件下坡面径流和降雨入渗的模拟互为条件,采用交替迭代的耦合分析方法较好地模拟了降雨坡面径流;Erduran[9]构建了一个综合数值模型,模拟了植被覆盖下地表径流与饱和地下水流之间的相互作用;朱磊等[10]基于地表径流物理机制和非饱和土壤水分运动理论,分别采用双层节点耦合法和整合离散方程的整体法,建立了降雨情况下地表径流与地下渗流的耦合模型;李根等[5]将地表径流模型与土壤水流模型进行耦合分析,计算了土坡降雨入渗,研究表明地表径流对土体入渗有较大提高;李尚辉等[11]建立了降雨边坡径流与大孔隙流的耦合模型,并基于有限元软件COMSOL实现了其数值求解．第二种求解方法是将降雨入渗面作为径流场与渗流场的内部域,实现坡面径流控制方程与坡体非饱和渗流控制方程的统一联立求解,该方法避免了求解过程中的流量交换,不需要进行交替迭代求解,在保证计算精度的同时极大地提高了计算效率．Kollet等[12]提出了一种新的耦合模型,将地表水系统视作地下水系统的上部边界,不需要考虑流量交换;童富果等[13]建立了降雨入渗与坡面径流直接耦合计算模型与求解方法,避免了两者间的迭代计算与流量交换,提高了计算效率和计算精度;Weill等[14]将扩散波方程作适当变换,使其在数学上与饱和-非饱和渗流Richards方程具有相同形式,并在多孔介质表面引入径流层概念,实现渗流区域和虚拟径流层的统一求解,提高了计算效率;Tian等[15-16]在童富果等的研究基础上进一步深入,并将耦合模型拓展至二维坡面径流与三维非饱和渗流的耦合;王乐等[17]建立了边坡渗流与坡面径流联合求解的三维模型,并对产流后的入渗边界流量进行了修正;Liu等[18]将水气两相流融入地表径流与地下水渗流的耦合模型,实现了边坡降雨-入渗-产流的数值模拟．

对于边坡径流与渗流耦合模型的求解,通常采用的数值方法包括有限差分法[6-7,12]、有限单元法[8,11,13-18]和有限体积法[5,9]等．然而,上述数值方法在遇到复杂求解域时存在网格生成困难的问题,而石根华博士[19]提出的数值流形法(NMM)以其独特的覆盖系统,对复杂边界具有极强的适应性,目前已被广泛应用于渗流计算领域[20-22]．本文基于前人的研究基础,对边坡径流与渗流的耦合模型进一步优化,采用压力水头形式的Richards方程描述地下水非饱和渗流,用运动波方程描述坡面径流,建立了耦合模型的变分提法,并基于NMM推导出离散求解格式,通过程序编制实现了边坡降雨-入渗-产流的全过程数值模拟,数值算例验证了算法及程序的有效性与正确性．

1 控制方程及其定解条件

边坡降雨-入渗-产流涉及三个方面的研究内容:降雨入渗、壤中流和坡面径流,是一个“三水”转换问题,其示意图如图1所示．降雨入渗到坡体之中,形成壤中流,其运移规律由饱和-非饱和渗流Richards方程[23-25]描述;当坡表土体达到饱和状态或者渗透能力小于降雨强度时,坡面径流产生,其运移规律由坡面径流Saint-Venant方程[6,9,26]描述．鉴于Saint-Venant方程求解的复杂性, 目前多采用其简化形式: 运动波或扩散波方程[27-29]．因此, “三水”问题的求解即可转换为Richards方程与Saint-Venant方程或者其简化形式的联立求解．

1.1 坡面径流控制方程及其定解条件

坡面径流由Saint-Venant方程描述,具体包含连续方程和动量方程．对连续方程,考虑坡度影响,可以表述为

(1)

式中,h为坡面径流水深;q为单宽流量;I=I(t)为降雨强度;f为地表入渗率;θ为坡面倾角;t为时间;l为坡顶沿坡面的长度．

当动量方程中只保留坡面比降So与摩阻比降Sf时,即可得运动波方程的水力学基础[29]:

So-Sf=0．

(2)

一般地,若考虑水流为紊流状态,则摩阻比降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 饱和-非饱和渗流控制方程及其定解条件

水在坡体中的运移规律由饱和-非饱和渗流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年代初首次提出,以拓扑流形和微分流形为基础,采用两套独立的覆盖系统(即数学覆盖和物理覆盖),能够实现对连续和非连续问题的统一求解．数学覆盖可视为一系列数学片的集合,数学片可以相互重叠且形状任意,无需与求解区域的物理边界重合,但要求为单连通区域,且其集合要完全覆盖整个求解区域．用各种物理边界(求解域边界和各种不连续面)依次切割数学覆盖中的各个数学片,并抛弃位于求解域之外的部分,就得到相应的物理片,所有物理片的集合就组成了物理覆盖,物理覆盖就与求解域边界精确匹配．而流形单元则是相应物理片之间进一步切割形成的互不重叠的部分．由上可知,NMM的前处理极为灵活,对复杂边界问题和不连续性问题具有先天优势．

本文以图2为例来展示NMM的覆盖系统．图中,材料界面将求解域Ω分成两个子域Ω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) 求解区域Ω (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 控制方程离散

对坡面径流的连续方程,考虑链式求导法则:

(14)

(15)

本文采用Galerkin加权余量法对相应控制方程进行离散,则构造式(15)的加权余量格式为

(16)

式中,g(x)表示测试函数;Γs代表降雨径流面．

对饱和-非饱和渗流控制方程式(8),构造其加权余量格式为

(17)

式中,G(x,y)表示测试函数．

对式(17)中的右端项,可通过分部积分降低阶次:

(18)

式中,∂Ω为求解区域Ω的外边界,∂Ω=Γh+Γq;ni为边界外法线向量的方向余弦．

将式(18)代入式(17),结合流量边界条件式(10b),并采用罚函数法施加本质边界条件和界面连续性条件,则式(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)中单元上的等效节点流量列阵可分别重写为

(22a)

(22b)

式中,F0为饱和-非饱和渗流问题中除降雨坡面边界外的其他边界形成的等效节点流量列阵．

将式(20a)与(20b)相加,并考虑式(22a)与(22b),即可得到降雨条件下坡面径流与坡体渗流全耦合模型的总体控制方程:

(23)

2.4 非线性迭代求解

对耦合模型的总体控制方程式(23),采用差分法对时间域进行离散化处理,有

(24)

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

(25)

式中,下标n表示时间步数;Δt为时间步长;θ为权重参数,0≤θ≤1,θ的取值不同,对应着不同的时间差分格式,也将直接影响解的精度和稳定性．研究表明,当θ=1,即对时间域的离散采用向后差分公式时,式(25)在求解理论上是无条件稳定的,因此本文亦采用向后差分公式．

由于总体控制方程式(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,持续时间为20 min．坡脚处安装监测系统,用于记录坡脚处的出流情况,监测总时长为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．

分别采用三角形和四边形网格形成的数学覆盖去覆盖整个求解区域,由于降雨时坡面处水力梯度变化剧烈,为提高计算精确度,本文对坡面处网格进行了加密,其几何模型和数学覆盖如图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

模拟时,分别采用三角形和四边形网格形成的数学覆盖去覆盖整个求解区域,两种覆盖方法生成的物理片数目均为3 233,生成的流形单元数目分别为5 600和2 800．计算得到的坡脚出流情况对比如图7所示,由图可以看出,本文两种数学覆盖的计算结果及Smith等[30]、Morita等[6]模拟结果与实测出流情况基本一致,再次验证了本文耦合模型的正确性．图8给出了x=5.6 m剖面上的土体饱和度随时间变化过程,模拟结果与实测数据存在一定差异,Smith等和Morita等的模拟结果也存在类似问题,分析其原因可能是初始时刻土壤含水率的实测数据较少,并不能真实反映实际含水情况,导致建立的土-水特征曲线和渗透系数函数模型并不能准确表征土体的物理特性．但从饱和度随时间的演化过程来看,其符合降雨入渗规律．因此,综上所述,可认为本文计算结果是客观合理的．

图7 坡脚出流过程对比

3.3 均质边坡降雨入渗

本小节对均质边坡的降雨入渗展开模拟．假定边坡水平长度为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．本小节仅采用四边形网格形成的数学覆盖去覆盖整个计算区域,共生成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左右,而边坡底层土体的体积含水量基本不发生变化．

图13 降雨结束时刻工况1边坡压力水头分布

4 结 论

本文采用一维运动波方程描述坡面径流,用二维压力水头形式的Richards方程描述土体非饱和渗流,考虑将降雨入渗面作为径流与渗流的内部域,推导出边坡降雨径流与渗流全耦合模型的总体控制方程．基于NMM对其进行数值离散,建立了对应的迭代求解格式,并通过编程实现了边坡降雨-入渗-产流的全过程数值模拟．研究表明,本文模型及计算方法准确可靠,能够较为真实地反映边坡的降雨入渗与产流过程,可为降雨型滑坡的稳定性评价及排水治理设计提供技术参考．