A Coupling Analysis of Rainfall Infiltration and Slope Surface Runoff Based on the Numerical Manifold Method
-
摘要: 降雨时坡地的入渗-产流分析，是降雨型滑坡、泥石流等地质灾害机理研究中的重要课题之一. 为实现边坡降雨-入渗-产流的全过程数值模拟，进一步提高计算效率，考虑将降雨入渗面视作坡面径流与坡体渗流的内部域，基于一维运动波方程和二维压力水头格式的Richards方程建立耦合模型，并推导出其总体控制方程，采用数值流形法(numerical manifold method, NMM)实现其数值求解，通过编制相应的计算程序分析了边坡降雨产流过程. 数值分析结果表明：所建模型的计算结果与试验数据及前人模拟结果吻合良好，验证了该文模型及计算方法的有效性与可靠性；降雨强度越大，产流时间越早，坡面积水深度越大，对坡体内的水分分布影响范围越广. 研究表明，所建模型能真实反映边坡降雨-入渗-产流全过程，可为降雨诱发的各类地质灾害分析提供计算依据.Abstract: The infiltration-runoff processes of slopes during rainfall were of significance in the mechanism study of rainfall-induced landslides, debris flows, and other geological disasters. To realize the numerical simulation of the whole rainfall infiltration-runoff process and further improve the computation efficiency, a coupling model and its total governing equations were derived from the 1D kinematic wave equation and the 2D h-based Richards' equation, with the rainfall infiltration surface deemed as the internal domain of the runoff and the seepage. Then the numerical manifold method (NMM) was used to solve the total governing equations, and the computation program was compiled to simulate the rainfall infiltration-runoff processes of slopes. The numerical analysis results showed that, the coupling model solution is in good agreement with the experimental data and previous results, verifying the validity and reliability of the proposed model. The higher the rainfall intensity is, the earlier the runoff producing time will be; the deeper the ponding depth is, and the wider the influence range on the water distribution within the slope will be. The proposed model can truly reflect the whole rainfall infiltration-runoff process of the slope, and provides a calculation basis for the analysis of various geological disasters induced by rainfall.
-
Key words:
- rainfall infiltration /
- slope runoff /
- saturated-unsaturated seepage /
- numerical manifold method
edited-byedited-by1) (我刊编委郑宏来稿) -
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 边坡降雨-入渗-产流示意图Figure 1. The diagrammatic sketch of rainfall-infiltration-runoff of the slope1.1 坡面径流控制方程及其定解条件
坡面径流由Saint-Venant方程描述，具体包含连续方程和动量方程. 对连续方程，考虑坡度影响，可以表述为
∂h∂t+∂q∂l=Icosθ−f, (1) 式中，h为坡面径流水深；q为单宽流量；I=I(t)为降雨强度；f为地表入渗率；θ为坡面倾角；t为时间；l为坡顶沿坡面的长度.
当动量方程中只保留坡面比降So与摩阻比降Sf时，即可得运动波方程的水力学基础[29]：
So−Sf=0. (2) 一般地，若考虑水流为紊流状态，则摩阻比降Sf可表示为
Sf=(uκhm)2, (3) 式中，κ为无量纲的摩擦因子；m为一实数. 因此，由式(3)可推导出坡面水流的速度公式为
u=κhmS1/2f=κhmS1/2o. (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公式，则一维坡面径流的运动波模型可用如下方程组描述：
{∂h∂t+∂q∂l=Icosθ−f,q=1nmh5/3S1/2o. (5) 方程组(5)的定解条件包括初始条件和边界条件：
1) 初始条件：若将降雨的开始时刻作为时间起点，则此时坡面各处均无径流产生，即
{h(l,t)|t=0=0,q(l,t)|t=0=0,0⩽l⩽L. (6) 2) 边界条件：坡顶处水深和流量均为0，即
{h(l,t)|l=0=0,q(l,t)|l=0=0,t>0. (7) 1.2 饱和-非饱和渗流控制方程及其定解条件
水在坡体中的运移规律由饱和-非饱和渗流Richards方程描述，本文采用其压力水头格式(h-form)，其二维表达式如下：
C(h)∂h∂t=∇[k(h)∇(h+y)], (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
h(x,y,t)=ˉh(x,y,t), on Γh. (10a) (b) 流量边界条件Γq
-k \frac{\partial(h+y)}{\partial \boldsymbol{n}}=\bar{q}(x, y, t), \quad \text { on } \varGamma_{\mathrm{q}} \text {. } (10b) (c) 材料界面条件Γm
\left\{\begin{array}{l} h^{+}(x, y, t)=h^{-}(x, y, t), \\ -k^{+} \frac{\partial(h+y)^{+}}{\partial \boldsymbol{n}}=-k^{-} \frac{\partial(h+y)^{-}}{\partial \boldsymbol{n}}, \quad \text { on } \varGamma_{\mathrm{m}} . \end{array}\right. (10c) 式(9)、(10)中，t0为初始时刻；h0为初始时刻的压力水头值；\bar{h}和\bar{q}分别为给定的压力水头值和流量值；n为边界单位外法线向量；上标“+”和“-”表示材料界面Γm两侧不同材料上的取值.
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).
基于上述理论，定义在流形单元上的压力水头近似函数可表示为
h^{\mathrm{h}}(\boldsymbol{r})=\sum\limits_{i=1}^{N_{\mathrm{p}}} w_{i}(\boldsymbol{r}) h_{i}(\boldsymbol{r}), (11) 式中，\boldsymbol{r}=(x, y)代表位置向量；Np表示覆盖该流形单元的物理片个数，若数学覆盖由三角形网格形成，则每个流形单元是3个物理片的交集，即Np=3，同样，若数学覆盖由四边形网格构成，则每个流形单元是4个物理片的交集，有Np=4；wi(r)表示第i个物理片上的单位分解函数，源自生成该物理片的数学片上的权函数；hi(r)表示定义于第i个物理片上的局部近似函数，可表达为
h_{i}(\boldsymbol{r})=\boldsymbol{p}^{\mathrm{T}}(\boldsymbol{r}) \boldsymbol{d}_{i}, (12) 式中，di表示定义在第i个物理片上的广义自由度向量；p(r)为完全多项式基函数，有零阶、一阶和二阶等多种形式，本文采用零阶基函数，则其数学表达式为p(r)=1T.
考虑采用零阶基函数，则压力水头近似函数可重新表示为
h^{\mathrm{h}}(\boldsymbol{r})=\sum\limits_{i=1}^{N_{\mathrm{p}}} w_{i}(\boldsymbol{r}) \boldsymbol{p}^{\mathrm{T}}(\boldsymbol{r}) \boldsymbol{d}_{i}=\boldsymbol{N} \boldsymbol{d}, (13) 式中, N=[N1, …, NNp]为形函数矩阵，其中，Ni=wi；d为广义自由度向量，d =[d1, …, dNp]T.
2.2 控制方程离散
对坡面径流的连续方程，考虑链式求导法则：
\frac{\partial q}{\partial l}=\frac{\partial q}{\partial h} \frac{\partial h}{\partial x} \frac{\partial x}{\partial l}=\frac{5 S_{o}^{1 / 2} h^{2 / 3} \cos \theta}{3 n_{\mathrm{m}}} \frac{\partial h}{\partial x}. (14) 记k_{\mathrm{o}}(h)=\frac{5 S_{\mathrm{o}}^{1 / 2} h^{2 / 3} \cos \theta}{3 n_{\mathrm{m}}}，则连续方程可重新表示为
\frac{\partial h}{\partial t}+k_{0}(h) \frac{\partial h}{\partial x}=I \cos \theta-f. (15) 本文采用Galerkin加权余量法对相应控制方程进行离散，则构造式(15)的加权余量格式为
\int_{\varGamma_{\mathrm{s}}}\left(\frac{\partial h}{\partial t}+k_{\mathrm{o}}(h) \frac{\partial h}{\partial x}\right) g \mathrm{~d} \varGamma=\int_{\varGamma_{\mathrm{s}}}(I \cos \theta-f) g \mathrm{~d} \varGamma, (16) 式中，g(x)表示测试函数；Γs代表降雨径流面.
对饱和-非饱和渗流控制方程式(8)，构造其加权余量格式为
\iint_{\varOmega} C(h) \frac{\partial h}{\partial t} G \mathrm{~d} \varOmega=\iint_{\varOmega} \nabla[k(h) \nabla(h+y)] G \mathrm{~d} \varOmega, (17) 式中，G(x, y)表示测试函数.
对式(17)中的右端项，可通过分部积分降低阶次：
\begin{gathered} \iint_{\varOmega} \nabla[k(h) \nabla(h+y)] G \mathrm{~d} \varOmega=\iint_{\varOmega} \frac{\partial}{\partial x_{i}}\left[k(h) \frac{\partial(h+y)}{\partial x_{i}}\right] G \mathrm{~d} \varOmega= \\ \oint_{\partial \varOmega} G k(h) \frac{\partial(h+y)}{\partial x_{i}} n_{i} \mathrm{~d} \varGamma-\iint_{\varOmega} \frac{\partial G}{\partial x_{i}} k(h) \frac{\partial(h+y)}{\partial x_{i}} \mathrm{~d} \varOmega, \end{gathered} (18) 式中，\partial \varOmega为求解区域Ω的外边界，\partial \varOmega=\varGamma_{\mathrm{h}}+\varGamma_{\mathrm{q}}; n_{i}为边界外法线向量的方向余弦.
将式(18)代入式(17)，结合流量边界条件式(10b)，并采用罚函数法施加本质边界条件和界面连续性条件，则式(17)可进一步表示为
\begin{gathered} \iint_{\varOmega} C(h) G \frac{\partial h}{\partial t} \mathrm{~d} \varOmega+\iint_{\varOmega} \nabla G k(h) \nabla h \mathrm{~d} \varOmega+k_{\mathrm{p}} \int_{\varGamma_{\mathrm{h}}} G(h-\bar{h}) \mathrm{d} \varGamma+k_{\mathrm{p}} \int_{\varGamma_{\mathrm{m}}} G\left(h^{+}-h^{-}\right) \mathrm{d} \varGamma= \\ \quad \iint_{\varOmega} \nabla G k(h) \nabla y \mathrm{~d} \varOmega+\int_{\varGamma_{\mathrm{q}}} G \bar{q} \mathrm{~d} \varGamma, \end{gathered} (19) 式中，kp为罚因子.
用式(13)中的近似场函数分别逼近式(16)与式(19)中的真实场函数，并采用Galerkin加权余量法，则可得到坡面径流和饱和-非饱和渗流的总体方程分别为
\boldsymbol{A} \frac{\partial \boldsymbol{h}}{\partial t}+\boldsymbol{D} \boldsymbol{h}=\boldsymbol{Q}, (20a) \boldsymbol{E} \frac{\partial \boldsymbol{h}}{\partial t}+\boldsymbol{K} \boldsymbol{h}=\boldsymbol{F}, (20b) 式中，h为广义自由度列阵；\frac{\partial \boldsymbol{h}}{\partial t}为广义自由度对时间的导数列阵；系数矩阵 A，D，E，K和等效节点流量列阵 Q，F均由相应的单元矩阵组装而成，其在相应单元上的具体表达式可表示如下：
\boldsymbol{A}^{e} =\int_{\varGamma_{\mathrm{s}}^{e}} \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N} \mathrm{d} \boldsymbol{\varGamma}, (21a) \boldsymbol{D}^{e} =\int_{\varGamma_{\mathrm{s}}^{e}} k_{\mathrm{o}} \boldsymbol{N}^{\mathrm{T}} \frac{\partial \boldsymbol{N}}{\partial x} \mathrm{~d} \boldsymbol{\varGamma}, (21b) \boldsymbol{Q}^{e} =\int_{\varGamma_{\mathrm{s}}^{e}} \boldsymbol{N}^{\mathrm{T}}(I \cos \theta-f) \mathrm{d} \boldsymbol{\varGamma}, (21c) \boldsymbol{E}^{e} =\int_{\varOmega^{e}} C \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N} \mathrm{d} \varOmega, (21d) \boldsymbol{K}^{e} =\int_{\varOmega^{e}} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{k} \boldsymbol{B} \mathrm{d} \varOmega+\int_{\varGamma_{\mathrm{h}}} k_{\mathrm{p}} \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N} \mathrm{d} \boldsymbol{\varGamma}+\int_{\varGamma_{\mathrm{m}}^{e}} k_{\mathrm{p}}\left(\boldsymbol{N}^{+}-\boldsymbol{N}^{-}\right)^{\mathrm{T}}\left(\boldsymbol{N}^{+}-\boldsymbol{N}^{-}\right) \mathrm{d} \boldsymbol{\varGamma}, (21e) \boldsymbol{F}^{e} =\int_{\varOmega^{e}} \boldsymbol{B}^{\mathrm{T}} \boldsymbol{k} \boldsymbol{B} \mathrm{d} \boldsymbol{\varOmega} \cdot \boldsymbol{y}+\int_{\varGamma_{\mathrm{q}}^{e}} \boldsymbol{N}^{\mathrm{T}} \bar{q} \mathrm{~d} \boldsymbol{\varGamma}+\int_{\varGamma_{\mathrm{h}}} k_{\mathrm{p}} \boldsymbol{N}^{\mathrm{T}} \bar{h} \mathrm{~d} \boldsymbol{\varGamma}, (21f) 式中，\boldsymbol{B}=\boldsymbol{L} \boldsymbol{N}, \boldsymbol{L}=[\partial / \partial x, \partial / \partial y]^{\mathrm{T}}; \boldsymbol{k}为渗透系数矩阵.
2.3 耦合模型建立
当对坡面径流与坡体渗流进行单独分析时，对坡面径流而言，需要给定坡面处的下渗分布情况；对于坡体的饱和-非饱和渗流而言，需要给定坡面处的入渗分布情况. 因此，耦合模型需要解决二者间流量交换的问题. 从现实角度来看，坡面处的下渗分布和入渗分布是一致的，如果将坡面边界作为内部域，就无须考虑二者之间的流量交换.
当坡体表面未出现径流时，降雨全部渗进边坡，此时仅需考虑饱和-非饱和渗流总体控制方程式(20b)；当坡面出现径流之后，坡面处的水深应满足式(20a)，整个坡体的压力水头应满足式(20b)，此时，坡面处的压力水头与径流水深在数值上应近似相等，即式(20a)与式(20b)求解变量相同. 又考虑坡面径流与坡体渗流共有坡面边界，为叙述方便，将坡面边界记为Γgs，假定坡面土体的入渗率为f，则式(20a)、(20b)中单元上的等效节点流量列阵可分别重写为
\boldsymbol{Q}^{e} =\int_{\varGamma_{\mathrm{gs}}^{e}} \boldsymbol{N}^{\mathrm{T}} I \cos \theta \mathrm{d} \boldsymbol{\varGamma}-\int_{\varGamma_{\mathrm{gs}}^{e}} \boldsymbol{N}^{\mathrm{T}} f \mathrm{~d} \boldsymbol{\varGamma}, (22a) \boldsymbol{F}^{e} =\boldsymbol{F}_{0}^{e}+\int_{\varGamma_{\mathrm{gs}}} \boldsymbol{N}^{\mathrm{T}} f \mathrm{~d} \boldsymbol{\varGamma}, (22b) 式中，F0为饱和-非饱和渗流问题中除降雨坡面边界外的其他边界形成的等效节点流量列阵.
将式(20a)与(20b)相加，并考虑式(22a)与(22b)，即可得到降雨条件下坡面径流与坡体渗流全耦合模型的总体控制方程：
\boldsymbol{M} \frac{\partial \boldsymbol{h}}{\partial t}+\boldsymbol{O} \boldsymbol{h}=\boldsymbol{P}, (23) 式中，\boldsymbol{M}=\boldsymbol{A}+\boldsymbol{E}; \boldsymbol{O}=\boldsymbol{D}+\boldsymbol{K}; \boldsymbol{P}=\boldsymbol{F}_{0}+\int_{\varGamma_{\mathrm{g}}^{e}} \boldsymbol{N}^{\mathrm{T}} I \cos \theta \mathrm{d} \boldsymbol{\varGamma}. 需要注意的是，对于耦合模型中的一维坡面径流问题，并不需要单独划分网格，而是直接将坡体渗流问题网格中的坡面降雨入渗边界作为径流计算网格. 当采用三角形或四边形网格形成的数学覆盖去覆盖整个求解区域时，则坡面上的每个径流单元也相应被三个或四个物理片所覆盖，类似于渗流问题的边界条件.
2.4 非线性迭代求解
对耦合模型的总体控制方程式(23)，采用差分法对时间域进行离散化处理，有
\left\{\begin{array}{l} \frac{\partial \boldsymbol{h}}{\partial t}=\frac{\boldsymbol{h}_{n+1}-\boldsymbol{h}_{n}}{\Delta t}, \\ \boldsymbol{h}=\theta \boldsymbol{h}_{n+1}+(1-\theta) \boldsymbol{h}_{n}. \end{array}\right. (24) 将式(24)代入式(23)，可得耦合模型的流形元迭代求解格式：
\left(\frac{\boldsymbol{M}}{\Delta t}+\theta \boldsymbol{O}\right) \boldsymbol{h}_{n+1}=\left[\frac{\boldsymbol{M}}{\Delta t}-(1-\theta) \boldsymbol{O}\right] \boldsymbol{h}_{n}+\boldsymbol{P}, (25) 式中，下标n表示时间步数；Δt为时间步长；θ为权重参数，0≤θ≤1，θ的取值不同，对应着不同的时间差分格式，也将直接影响解的精度和稳定性. 研究表明，当θ=1，即对时间域的离散采用向后差分公式时，式(25)在求解理论上是无条件稳定的，因此本文亦采用向后差分公式.
由于总体控制方程式(25)具有强烈的非线性，需要采用非线性迭代方法进行求解. 本文采用Picard迭代法进行迭代求解，取θ=1，则迭代求解格式(25)可进一步表示为
\left(\frac{\boldsymbol{M}_{n+1, m}}{\Delta t}+\boldsymbol{O}_{n+1, m}\right) \boldsymbol{h}_{n+1, m}=\frac{\boldsymbol{M}_{n+1, m}}{\Delta t} \boldsymbol{h}_{n}+\boldsymbol{P}_{n+1, m}, (26) 式中，m表示迭代步数；hn+1, m+1表示第n+1个时间步中第m+1个迭代步的压力水头列阵，其余下标含义与此类似.
当绝对误差达到给定的容差ε时，迭代终止，即
\left\|\boldsymbol{h}_{n+1, m+1}-\boldsymbol{h}_{n+1, m}\right\| \leqslant \boldsymbol{\varepsilon}, (27) 式中，‖ · ‖表示∞-范数. 当收敛条件满足时，令\boldsymbol{h}_{n+1}=\boldsymbol{h}_{n+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-5 m/s，持续时间为20 min. 坡脚处安装监测系统，用于记录坡脚处的出流情况，监测总时长为25 min.
计算时，土-水特征曲线和渗透性函数分别由van Genuchten模型[31]和Mualem模型[32]来描述：
\theta=\theta_{\mathrm{r}}+\frac{\theta_{\mathrm{s}}-\theta_{\mathrm{r}}}{\left[1-(-a h)^{b}\right]^{c}}, (28) k_{\mathrm{r}}=\varTheta^{1 / 2}\left[1-\left(1-\varTheta^{1 / c}\right)^{c}\right]^{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可以看出，数值模拟结果的产流时间均要早于试验结果，分析其原因为数值模型中均没有考虑土壤中空气的压缩性.
图 4 Abdul-Gillham试验几何模型及数学覆盖Figure 4. Geometric models and mathematical covers of the Abdul-Gillham system3.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×105 cm-1 · s-1.
图 6 Smith试验计算几何模型注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同.Figure 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等的模拟结果也存在类似问题，分析其原因可能是初始时刻土壤含水率的实测数据较少，并不能真实反映实际含水情况，导致建立的土-水特征曲线和渗透系数函数模型并不能准确表征土体的物理特性. 但从饱和度随时间的演化过程来看，其符合降雨入渗规律. 因此，综上所述，可认为本文计算结果是客观合理的.
图 8 x=5.6 m剖面土体饱和度演化过程Figure 8. Evolutions of the saturation degree along the profile for x=5.6 m3.3 均质边坡降雨入渗
本小节对均质边坡的降雨入渗展开模拟. 假定边坡水平长度为100 m，竖直厚度为40 m，坡角为15°，其几何模型如图 9所示. 土体的饱和渗透系数为4.332×10－4 m/min，土-水特征曲线及渗透性函数如图 10所示. 假定初始时刻土体的体积含水量为0.045，坡面Manning粗糙度系数为0.035 m－1/3 · min. 计算时，分别采用表 1中5种不同工况的降雨强度进行模拟分析，以探究降雨强度对坡面产流的影响，降雨历时10 h. 本小节仅采用四边形网格形成的数学覆盖去覆盖整个计算区域，共生成1 458个物理片和1 374个流形单元.
图 10 土-水特征曲线及渗透性函数Figure 10. The soil-water characteristic curve and the permeability function表 1 降雨工况Table 1. The rainfall intensitiescase number 1 2 3 4 5 rainfall intensity I/(m/min) 2.502×10－3 2.085×10－3 1.668×10－3 1.251×10－3 8.340×10－4 图 11和12分别为5种不同降雨工况计算得到的平衡条件坡面水位线图和坡脚出流情况. 由图可知，降雨强度越大，坡面产流时间越早，达到平衡条件时的坡面积水越深，坡脚径流量也越大. 工况1降雨强度最大，产流时间最早，约为25 min，达到平衡条件时坡脚积水深度为8.5 cm；而工况5降雨强度最小，降雨约79 min之后坡面才开始产流，达到平衡条件时的坡脚积水深度也仅有4.3 cm.
图 13和14分别给出了工况1情况下降雨结束时刻坡体和坡面的压力水头分布. 由图可以看出，由于初始时刻边坡的体积含水量较低，降雨只对边坡表面产生影响，最大影响深度约为1 m左右，而边坡底层土体的体积含水量基本不发生变化.
图 13 降雨结束时刻工况1边坡压力水头分布Figure 13. The pressure head distribution of the slope in case 1 at the end of rainfall
图 14 降雨结束时刻工况1边坡表层压力水头分布(单位: m)Figure 14. The pressure head distribution at the surface of the slope in case 1 at the end of rainfall (unit: m)4. 结论
本文采用一维运动波方程描述坡面径流，用二维压力水头形式的Richards方程描述土体非饱和渗流，考虑将降雨入渗面作为径流与渗流的内部域，推导出边坡降雨径流与渗流全耦合模型的总体控制方程. 基于NMM对其进行数值离散，建立了对应的迭代求解格式，并通过编程实现了边坡降雨-入渗-产流的全过程数值模拟. 研究表明，本文模型及计算方法准确可靠，能够较为真实地反映边坡的降雨入渗与产流过程，可为降雨型滑坡的稳定性评价及排水治理设计提供技术参考.
-
图 1 边坡降雨-入渗-产流示意图
Figure 1. The diagrammatic sketch of rainfall-infiltration-runoff of the slope
图 4 Abdul-Gillham试验几何模型及数学覆盖
Figure 4. Geometric models and mathematical covers of the Abdul-Gillham system
图 6 Smith试验计算几何模型
注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同.
Figure 6. The geometric model for the Smith system
图 8 x=5.6 m剖面土体饱和度演化过程
Figure 8. Evolutions of the saturation degree along the profile for x=5.6 m
图 10 土-水特征曲线及渗透性函数
Figure 10. The soil-water characteristic curve and the permeability function
图 13 降雨结束时刻工况1边坡压力水头分布
Figure 13. The pressure head distribution of the slope in case 1 at the end of rainfall
图 14 降雨结束时刻工况1边坡表层压力水头分布(单位: m)
Figure 14. The pressure head distribution at the surface of the slope in case 1 at the end of rainfall (unit: m)
表 1 降雨工况
Table 1. The rainfall intensities
case number 1 2 3 4 5 rainfall intensity I/(m/min) 2.502×10－3 2.085×10－3 1.668×10－3 1.251×10－3 8.340×10－4 
下载: 导出CSV
-
[1] 姚海林, 郑少河, 李文斌, 等. 降雨入渗对非饱和膨胀土边坡稳定性影响的参数研究[J]. 岩石力学与工程学报, 2002, 21(7): 1034-1039.YAO Hailin, ZHENG Shaohe, LI Wenbin, et al. Parametric study on the effect of rain infiltration on stability of unsaturated expansive soil slope[J]. Chinese Journal of Rock Mechanics and Engineering, 2002, 21(7): 1034-1039. (in Chinese) [2] 董建军, 王思萌, 杨晓萧, 等. 基于非饱和-饱和渗流的降雨入渗边坡稳定性分析[J]. 水利与建筑工程学报, 2018, 16(6): 99-104.DONG Jianjun, WANG Simeng, YANG Xiaoxiao, et al. Analysis of the rainfall infiltration on slope stability based on unsaturated-saturated seepage[J]. Journal of Water Resources and Architectural Engineering, 2018, 16(6): 99-104. (in Chinese) [3] XIONG X, SHI Z, XIONG Y, et al. Unsaturated slope stability around the Three Gorges Reservoir under various combinations of rainfall and water level fluctuation[J]. Engineering Geology, 2019, 261: 105231. doi: 10.1016/j.enggeo.2019.105231 [4] WALLACH R, GRIGORIN G, RIVLIN J. The errors in surface runoff prediction by neglecting the relationship between infiltration rate and overland flow depth[J]. Journal of Hydrology, 1997, 200: 243-259. doi: 10.1016/S0022-1694(97)00017-6 [5] 李根, 韩同春, 吴俊杨, 等. 基于有限体积法的地表径流与土壤水流耦合分析[J]. 浙江大学学报(工学版), 2022, 56(5): 947-955.LI Gen, HAN Tongchun, WU Junyang, et al. Coupled analysis on surface runoff and soil water movement by finite volume method[J]. Journal of Zhejiang University (Engineering Science), 2022, 56(5): 947-955. (in Chinese) [6] MORITA M, YEN B C. Modeling of conjunctive two-dimensional surface-three-dimensional subsurface flows[J]. Journal of Hydraulic Engineering, 2002, 128(2): 184-200. doi: 10.1061/(ASCE)0733-9429(2002)128:2(184) [7] PANDAY S, HUYAKORN P S. A fully coupled physically-based spatially-distributed model for evaluating surface/subsurface flow[J]. Advances in Water Resources, 2004, 27(4): 361-382. doi: 10.1016/j.advwatres.2004.02.016 [8] 张培文, 刘德富, 郑宏, 等. 降雨条件下坡面径流和入渗耦合的数值模拟[J]. 岩土力学, 2004, 25(1): 109-113.ZHANG Peiwen, LIU Defu, ZHENG Hong, et al. Coupling numerical simulation of slope runoff and infiltration under rainfall conditions[J]. Rock and Soil Mechanics, 2004, 25(1): 109-113. (in Chinese) [9] 埃尔杜拉K S. 一个综合植被地表径流与饱和地下水流之间相互作用的数值模型[J]. 应用数学和力学, 2012, 33(7): 828-844. doi: 10.3879/j.issn.1000-0887.2012.07.004ERDURAN K S. An integrated numerical model for vegetated surface and saturated subsurface flow interaction[J]. Applied Mathematics and Mechanics, 2012, 33(7): 828-844. (in Chinese) doi: 10.3879/j.issn.1000-0887.2012.07.004 [10] 朱磊, 田军仓, 孙骁磊. 基于全耦合的地表径流与土壤水分运动数值模拟[J]. 水科学进展, 2015, 26(3): 322-330.ZHU Lei, TIAN Juncang, SUN Xiaolei. A fully couple numerical simulation for surface runoff and soil water movement[J]. Advances in Water Science, 2015, 26(3): 322-330. (in Chinese) [11] 李尚辉, 刘代文, 阙云, 等. 坡面径流与大孔隙流耦合作用下边坡水分场参数敏感性分析[J]. 济南大学学报(自然科学版), 2023, 37(3): 264-272.LI Shanghui, LIU Daiwen, QUE Yun, et al. Parameter sensitivity analysis on slope water fields under slope runoffs and macropore flows coupling[J]. Journal of University of Jinan (Science and Technology), 2023, 37(3): 264-272. (in Chinese) [12] KOLLET S J, MAXWELL R M. Integrated surface-groundwater flow modeling: a free-surface overland flow boundary condition in a parallel groundwater flow model[J]. Advances in Water Resources, 2006, 29(7): 945-958. doi: 10.1016/j.advwatres.2005.08.006 [13] 童富果, 田斌, 刘德富. 改进的斜坡降雨入渗与坡面径流耦合算法研究[J]. 岩土力学, 2008, 29(4): 1035-1040.TONG Fuguo, TIAN Bin, LIU Defu. A coupling analysis of slope runoff and infiltration under rainfall[J]. Rock and Soil Mechanics, 2008, 29(4): 1035-1040. (in Chinese) [14] WEILL S, MOUCHE E, PATIN J. A generalized Richards equation for surface/subsurface flow modelling[J]. Journal of Hydrology, 2009, 366(1/4): 9-20. [15] TIAN Dongfang, LIU Defu. A new integrated surface and subsurface flows model and its verification[J]. Applied Mathematical Modelling, 2011, 35: 3574-3586. doi: 10.1016/j.apm.2011.01.035 [16] TIAN Dongfang, ZHENG Hong, LIU Defu. A 2D integrated FEM model for surface water-groundwater flow of slopes under rainfall condition[J]. Landslides, 2017, 14(2): 577-593. doi: 10.1007/s10346-016-0716-4 [17] 王乐, 田东方. 边坡渗流与坡面径流联合求解三维有限元模型[J]. 人民珠江, 2019, 40(5): 24-29.WANG Le, TIAN Dongfang. 3D FEM integrated model for seepage and runoff process[J]. Pearl River, 2019, 40(5): 24-29. (in Chinese) [18] LIU Gang, TONG Fuguo, TIAN Bin. A finite element model for simulating surface run-off and unsaturated seepage flow in the shallow subsurface[J]. Hydrological Processes, 2019, 33(26): 3378-3390. doi: 10.1002/hyp.13564 [19] SHI Genhua. Manifold method of material analysis[C]//Transactions of the 9th Army Conference on Applied Mathematics and Computing. US Army Research Office, 1991. [20] 陈远强, 杨永涛, 郑宏, 等. 饱和-非饱和渗流的数值流形法研究与应用[J]. 岩土工程学报, 2019, 41(2): 338-347.CHEN Yuanqiang, YANG Yongtao, ZHENG Hong, et al. Saturated-unsaturated seepage by numerical manifold method[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(2): 338-347. (in Chinese) [21] LIN Shan, ZHENG Hong, ZHANG Zhihong, et al. Evaluation of permeability of soil & rock aggregate using meshless numerical manifold method[J]. Computers and Geotechnics, 2022, 151: 104953. doi: 10.1016/j.compgeo.2022.104953 [22] CHEN Yuanqiang, ZHENG Hong, YIN Boyuan, et al. The MLS-based numerical manifold method for Darcy flow in heterogeneous porous media[J]. Engineering Analysis With Boundary Elements, 2023, 148: 220-242. doi: 10.1016/j.enganabound.2022.12.030 [23] RICHARDS L A. Capillary conduction of liquids through porous mediums[J]. Journal of Applied Physics, 1931, 1: 318-333. http://ci.nii.ac.jp/naid/10006144118 [24] 朱帅润, 李绍红, 钟彩尹, 等. 时间分数阶的非饱和渗流数值分析及其应用[J]. 应用数学和力学, 2022, 43(9): 966-975. doi: 10.21656/1000-0887.420334ZHU Shuairun, LI Shaohong, ZHONG Caiyin, et al. Numerical analysis of time fractional-order unsaturated flow and its application[J]. Applied Mathematics and Mechanics, 2022, 43(9): 966-975. (in Chinese) doi: 10.21656/1000-0887.420334 [25] 刘能胜, 曹恒明. 通用函数边界下潜水非稳定流模型的解及应用[J]. 应用数学和力学, 2020, 41(9): 1048-1056. doi: 10.21656/1000-0887.400371LIU Nengsheng, CAO Hengming. Solution and application of the transient phreatic flow motion model under general function boundary[J]. Applied Mathematics and Mechanics, 2020, 41(9): 1048-1056. (in Chinese) doi: 10.21656/1000-0887.400371 [26] 郑素佩, 李霄, 赵青宇, 等. 求解二维浅水波方程的旋转混合格式[J]. 应用数学和力学, 2022, 43(2): 176-186. doi: 10.21656/1000-0887.400124ZHENG Supei, LI Xiao, ZHAO Qingyu, et al. A rotated mixed scheme for solving 2D shallow water equation[J]. Applied Mathematics and Mechanics, 2022, 43(2): 176-186. (in Chinese) doi: 10.21656/1000-0887.400124 [27] LIGHTHILL M J, WHITHAM G B. On kinematic waves: flood movement in long rivers[J]. Proceedings of the Royal Society of London (Series A): Mathematical and Physical Sciences, 1955, 229(1178): 281-316. doi: 10.1098/rspa.1955.0088 [28] PONCE V M, SIMONS D B, LI R M. Applicability of kinematic and diffusion models[J]. Journal of the Hydraulics Division, 1978, 104(3): 353-360. doi: 10.1061/JYCEAJ.0004958 [29] 文康, 金管生, 李蝶娟. 地表径流过程的数学模拟[M]. 北京: 水利电力出版社, 1991.WEN Kang, JIN Guansheng, LI Diejuan. Mathematical Simulation of Surface Runoff Processes[M]. Beijing: China Water & Power Press, 1991. (in Chinese) [30] SMITH R E, WOOLHISER D A. Mathematical simulation of infiltrating watersheds[R]. Colorado: Colorado State University, 1971. [31] VAN GENUCHTEN M T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5): 892-898. doi: 10.2136/sssaj1980.03615995004400050002x [32] MUALEM Y. A new model for predicting the hydraulic conductivity of unsaturated porous media[J]. Water Resources Research, 1976, 12(3): 513-522. doi: 10.1029/WR012i003p00513 [33] GARDNER W R. Some steady-state solutions of the unsaturated moisture flow equation with application to evaporation from a water table[J]. Soil science, 1958, 85(4): 228-232. doi: 10.1097/00010694-195804000-00006 [34] BROOKS R A, COREY A T. Hydraulic properties of porous media[R]. Hydrology Papers, 1964. [35] ABDUL A S, GILLHAM R W. Laboratory studies of the effects of the capillary fringe on streamflow generation[J]. Water Resources Research, 1984, 20(6): 1-18. [36] VANDERKWAAK J E. Numerical simulation of flow and chemical transport in integrated surface-subsurface hydrologic systems[D]. Waterloo: University of Waterloo, 1999. [37] SINGH V, BHALLAMUDI S M. Conjunctive surface-subsurface modeling of overland flow[J]. Advances in Water Resources, 1998, 21(7): 567-579. doi: 10.1016/S0309-1708(97)00020-1 期刊类型引用(0)
其他类型引用(1)
-
下载:
计量
- 文章访问数: 490
- HTML全文浏览量: 166
- PDF下载量: 52
- 被引次数: 1
渝公网安备50010802005915号