工程科学学报，第37卷，第4期：407-413,2015年4月 Chinese Journal of Engineering,Vol.37,No.4:407-413,April 2015 DOI:10.13374/j.issn2095-9389.2015.04.002:http://journals.ustb.edu.cn 考虑二维降雨入渗的非饱和土边坡瞬态体积含水率 分析 付建新2》，宋卫东12)四，杜建华12) 1)北京科技大学金属矿山高效开采与安全教有部重点实验室，北京1000832)北京科技大学土木与环境工程学院，北京100083 ☒通信作者：E-mail:songwd@usth.cdu.cm 摘要降雨是非饱和土边坡失稳的关键因素之一，非饱和土体的强度与体积含水率有密切的关系.考虑二维降雨入渗条 件，基于非饱和土达西渗流定理及质量守恒定律建立了非饱和土的二维降雨入渗模型方程.基于交替隐式有限差分法，采用 MATLAB编制了计算程序，对方程进行了数值计算，研究了不同降雨条件、边坡内部不同位置的瞬态体积含水率分布.研究 表明，由于考虑了水平方向的入渗，边坡内随着深度和水平距离的增加，体积含水率都发生了变化，在距离边坡表面1m左右 时，体积含水率变化率变小，边坡表层的体积含水率变化远远大于边坡深部.初始体积含水率越小，边坡体积含水率变化越 大，这对边坡的稳定性不利. 关键词非饱和土边坡；边坡稳定性分析：水含量：数学模型：交替方向隐式法：有限差分法 分类号TU42 Transient volume water content analysis of unsaturated soil slopes considering two-dimensional rainfall infiltration FU Jian-xin,SONG Wei-dong,DU Jian-hua 1)Key Laboratory of High-Efficient Mining and Safety of Metal Mines (Ministry of Education),University of Science and Technology Beijing,Beijing 100083.China 2)School of Civil and Environmental Engineering.University of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:songwd@ustb.edu.cn ABSTRACT Rainfall is one of the key factors of unsaturated soil slope instability and the strength of unsaturated soils is closely related to volume water content.Considering two-dimensional rainfall infiltration,on the basis of the Darcy seepage theorem and mass conservation law of unsaturated soils,this paper constructed a two-dimensional infiltration model of unsaturated soils and coded a MATLAB program based on an alternating-direction-implicit (ADI)finite difference method which is used to numerical calculations to study the transient volume water content distribution at different slope locations under different rainfall conditions.The results show that the volume water content changes with the horizontal distance and depth.Moreover,the volume water content at the slope surface changes far greater than that at the deep part.The smaller the initial volume water content is,the greater the change in volume water content of the slope is,which is disadvantageous to slope stability. KEY WORDS unsaturated soil slopes;slope stability analysis;water content:mathematical models:alternating direction implicit method;finite difference method 非饱和土是一种典型的多孔介质四，由于多孔介质中空裂隙的大小、形状复杂，因此地下水在土体中的 收稿日期：2013-10-29 基金项目：国家自然科学基金资助项目(5137403)：教有部科技发展中心资助项目(20120006110022)
工程科学学报，第 37 卷，第 4 期: 407--413，2015 年 4 月 Chinese Journal of Engineering，Vol． 37，No． 4: 407--413，April 2015 DOI: 10． 13374 /j． issn2095--9389． 2015． 04． 002; http: / /journals． ustb． edu． cn 考虑二维降雨入渗的非饱和土边坡瞬态体积含水率 分析 付建新1，2) ，宋卫东1，2) ，杜建华1，2) 1) 北京科技大学金属矿山高效开采与安全教育部重点实验室，北京 100083 2) 北京科技大学土木与环境工程学院，北京 100083 通信作者: E-mail: songwd@ ustb． edu． cn 摘 要 降雨是非饱和土边坡失稳的关键因素之一，非饱和土体的强度与体积含水率有密切的关系． 考虑二维降雨入渗条 件，基于非饱和土达西渗流定理及质量守恒定律建立了非饱和土的二维降雨入渗模型方程． 基于交替隐式有限差分法，采用 MATLAB 编制了计算程序，对方程进行了数值计算，研究了不同降雨条件、边坡内部不同位置的瞬态体积含水率分布． 研究 表明，由于考虑了水平方向的入渗，边坡内随着深度和水平距离的增加，体积含水率都发生了变化，在距离边坡表面 1 m 左右 时，体积含水率变化率变小，边坡表层的体积含水率变化远远大于边坡深部． 初始体积含水率越小，边坡体积含水率变化越 大，这对边坡的稳定性不利． 关键词 非饱和土边坡; 边坡稳定性分析; 水含量; 数学模型; 交替方向隐式法; 有限差分法 分类号 TU42 Transient volume water content analysis of unsaturated soil slopes considering two-dimensional rainfall infiltration FU Jian-xin1，2) ，SONG Wei-dong1，2) ，DU Jian-hua1，2) 1) Key Laboratory of High-Efficient Mining and Safety of Metal Mines ( Ministry of Education) ，University of Science and Technology Beijing，Beijing 100083，China 2) School of Civil and Environmental Engineering，University of Science and Technology Beijing，Beijing 100083，China Corresponding author，E-mail: songwd@ ustb． edu． cn ABSTＲACT Ｒainfall is one of the key factors of unsaturated soil slope instability and the strength of unsaturated soils is closely related to volume water content． Considering two-dimensional rainfall infiltration，on the basis of the Darcy seepage theorem and mass conservation law of unsaturated soils，this paper constructed a two-dimensional infiltration model of unsaturated soils and coded a MATLAB program based on an alternating-direction-implicit ( ADI) finite difference method which is used to numerical calculations to study the transient volume water content distribution at different slope locations under different rainfall conditions． The results show that the volume water content changes with the horizontal distance and depth． Moreover，the volume water content at the slope surface changes far greater than that at the deep part． The smaller the initial volume water content is，the greater the change in volume water content of the slope is，which is disadvantageous to slope stability． KEY WOＲDS unsaturated soil slopes; slope stability analysis; water content; mathematical models; alternating direction implicit method; finite difference method 收稿日期: 2013--10--29 基金项目: 国家自然科学基金资助项目( 5137403) ; 教育部科技发展中心资助项目( 20120006110022) 非饱和土是一种典型的多孔介质［1］，由于多孔介 质中空裂隙的大小、形状复杂，因此地下水在土体中的
·408· 工程科学学报，第37卷，第4期 运动规律复杂，有些地方甚至不连续.可见，研究地下 式中，0为体积含水率，K()和K(山)分别是以土体 水渗流问题不能直接考察水质点的运动轨迹 积含水率和基质势为自变量的非饱和土导水率，业为 非饱和土是一种复杂的三相土，降雨入渗实质上 土水势. 是水分在土体饱气带中的运动6，由于非饱和土水 对于饱和土来说，导水率K为常数，而非饱和土 分渗流问题的复杂性，现有的研究往往只考虑了一维 的导水率K是土的基质势或者体积含水率的函数，并 入渗情况下土体体积含水率的变化规律0，即只考 随着体积含水率和基质势的减小而减小.将式(2)代 虑了水分的垂直入渗，而忽略了水平入渗的影响. 入式(4)中，假设重力势只与深度有关，可得 李宁和许健聪四、詹良通等四采用Fourier积分 q=-K()V(wm±z). (5) 变换分别对不同降雨强度情况下，斜长边坡的降雨入 式中“±”由z来确定，z向上为正时取“+”，z向下为 渗解析解进行了推导，并建立降雨入渗模型，这种方法 正时取“-” 不能直观反映土体含水量随雨水入渗的变化情况，也 1.3非饱和土边坡二维入渗运动方程 没有考虑到土体的渗透系数和扩散率等水分运动参数 非饱和土边坡在二维降雨入渗时包括两个过程， 随含水量的变化而对入渗过程的影响，只是以饱和的 一个是降水从地表垂直向下进入边坡的垂直入渗问 渗透系数进行研究，这与雨水非饱和土中入渗的实际 题，一个则是侧向入渗问题，如图1所示.斜坡在降雨 过程不符. 过程中，边坡表层单元A会同时有垂直通量V和水平 本文综合考虑降雨过程中，水分的垂直和水平入 通量V 渗的相互影响，研究二维入渗情况下土体瞬态体积含 水率的分布规律.不考虑空气压力、土骨架变形、温度 降雨入渗 及水中的溶质等因素，基于水分基本运动方程建立降 雨入渗模型，在此基础上进行网格化，基于交替隐式差 分法，采用MATBLAB编制计算程序，最后以工程实例 垂直通量V 进行验算，研究在不同条件下的土体瞬态体积含水率 在二维降雨入渗条件下的分布情况. 水平通量V 1非饱和土边坡二维入渗模型 1.1非饱和土的土水势 图1边坡二维降雨入渗示意图 非饱和土是由固体颗粒、水及空气构成的三相体 Fig.I Schematie diagram of two-dimensional infiltration into a slope 土体中水的流动是由土水势所决定的.为了建立非饱 取单元A为研究对象，单元体的入渗情况如图2 和土体的水分运动方程，必须要考虑土水势的影响. 所示，单元体水平和竖直分别为dx和dz,另一方向为 一般而言，土水势是指土中某点的水分状态与标 单位长度1，图中所示方向为正. 准状态的势能差团，即物理点体积元内单位数量的土 中水分具有的势能，由下式表示： 中=中。+中。+中m+.+ (1) 式中，山为总土水势，山，为重力势，中。为压力势，中m为 基质势，山，为溶质势，山，为温度势 对于非饱和土，一般不考虑压力势、溶质势及温度 势，因此非饱和土的土水势如下式所示： 少=中m土中 (2) 式中，基质势与土体体积含水率密切相关，随着土体的 饱和程度的增加而降低. 图2单元体 1.2非饱和土达西定律 Fig.2 Element Richards首先将达西定律引入了非饱和土水分运 假设土体单元A的大小不随时间改变，则时间 动规律的研究中，非饱和土中水分运动的方程可表 内，通过单元体的水量为 示为 (3) (ay+ay drd-dr.. (6) q=-K(山)7山 dx dz 或 另一方面，假设单元体内体积含水率为，则单元 q=-K(o)7山. (4) 体内水量为xdz,在单位时间d:内单元体水量变化
工程科学学报，第 37 卷，第 4 期 运动规律复杂，有些地方甚至不连续． 可见，研究地下 水渗流问题不能直接考察水质点的运动轨迹［2--5］． 非饱和土是一种复杂的三相土，降雨入渗实质上 是水分在土体饱气带中的运动［6--8］，由于非饱和土水 分渗流问题的复杂性，现有的研究往往只考虑了一维 入渗情况下土体体积含水率的变化规律［9--10］，即只考 虑了水分的垂直入渗，而忽略了水平入渗的影响． 李宁和许健聪［11］、詹良通等［12］采用 Fourier 积分 变换分别对不同降雨强度情况下，斜长边坡的降雨入 渗解析解进行了推导，并建立降雨入渗模型，这种方法 不能直观反映土体含水量随雨水入渗的变化情况，也 没有考虑到土体的渗透系数和扩散率等水分运动参数 随含水量的变化而对入渗过程的影响，只是以饱和的 渗透系数进行研究，这与雨水非饱和土中入渗的实际 过程不符． 本文综合考虑降雨过程中，水分的垂直和水平入 渗的相互影响，研究二维入渗情况下土体瞬态体积含 水率的分布规律． 不考虑空气压力、土骨架变形、温度 及水中的溶质等因素，基于水分基本运动方程建立降 雨入渗模型，在此基础上进行网格化，基于交替隐式差 分法，采用 MATBLAB 编制计算程序，最后以工程实例 进行验算，研究在不同条件下的土体瞬态体积含水率 在二维降雨入渗条件下的分布情况． 1 非饱和土边坡二维入渗模型 1. 1 非饱和土的土水势 非饱和土是由固体颗粒、水及空气构成的三相体． 土体中水的流动是由土水势所决定的． 为了建立非饱 和土体的水分运动方程，必须要考虑土水势的影响． 一般而言，土水势是指土中某点的水分状态与标 准状态的势能差［13］，即物理点体积元内单位数量的土 中水分具有的势能，由下式表示: ψ = ψg + ψp + ψm + ψs + ψt ． ( 1) 式中，ψ 为总土水势，ψg 为重力势，ψp 为压力势，ψm 为 基质势，ψs 为溶质势，ψt 为温度势． 对于非饱和土，一般不考虑压力势、溶质势及温度 势，因此非饱和土的土水势如下式所示: ψ = ψm ± ψg ． ( 2) 式中，基质势与土体体积含水率密切相关，随着土体的 饱和程度的增加而降低． 1. 2 非饱和土达西定律 Ｒichards 首先将达西定律引入了非饱和土水分运 动规律的研究中［14］，非饱和土中水分运动的方程可表 示为 q = － K( ψm ) Δ ψ ( 3) 或 q = － K( θ) Δ ψ． ( 4) 式中，θ 为体积含水率，K( θ) 和 K( ψm ) 分别是以土体 积含水率和基质势为自变量的非饱和土导水率，ψ 为 土水势． 对于饱和土来说，导水率 K 为常数，而非饱和土 的导水率 K 是土的基质势或者体积含水率的函数，并 随着体积含水率和基质势的减小而减小． 将式( 2) 代 入式( 4) 中，假设重力势只与深度有关，可得 q = － K( θ) Δ ( ψm ± z) ． ( 5) 式中，“± ”由 z 来确定，z 向上为正时取“+ ”，z 向下为 正时取“－ ”． 1. 3 非饱和土边坡二维入渗运动方程 非饱和土边坡在二维降雨入渗时包括两个过程， 一个是降水从地表垂直向下进入边坡的垂直入渗问 题，一个则是侧向入渗问题，如图 1 所示． 斜坡在降雨 过程中，边坡表层单元 A 会同时有垂直通量 Vwy和水平 通量 Vwx ． 图 1 边坡二维降雨入渗示意图 Fig． 1 Schematic diagram of two-dimensional infiltration into a slope 取单元 A 为研究对象，单元体的入渗情况如图 2 所示，单元体水平和竖直分别为 dx 和 dz，另一方向为 单位长度 1，图中所示方向为正． 图 2 单元体 Fig． 2 Element 假设土体单元 A 的大小不随时间改变，则 dt 时间 内， ( 通过单元体的水量为 Vwx x + Vwz ) z dxdzdt． ( 6) 另一方面，假设单元体内体积含水率为 θ，则单元 体内水量为 θdxdz，在单位时间 dt 内单元体水量变化 · 804 ·
付建新等：考虑二维降雨入渗的非饱和土边坡瞬态体积含水率分析 ·409 率为 形成了地表积水或者径流，地面体积含水率可接近饱 a0 dxdzdt. (7) 和体积含水率日，的值并保持不变，此时土体表面处边 at 界条件为： 由质量守恒可得 )()coa (8) R(t）cosa,z=0,0<1<t.: 将式(4)代入式(8)，可得单元体二维渗流方程为 0=0.z=0,1>t 架k@普]+[ko警】冰0 式中，t.为降雨强度超过土的入渗能力的时刻，R()为 降雨强度，0，为饱和体积含水率. (9) (2)下边界条件.假设地下水埋藏较深，土壤深 D()和K()之间的关系为阿 处含水量恒定且分布均匀，即 D(o)=k(0) (10) 0。=0(x,z）,z=L,t>0. C(6) 式中，。为初始体积含水率，L为计算区域的下边界. 式中C(0)为比水容重，定义为 c-总 (3)右边界条件.为了简便的说明问题，假定远 (11 离边坡边界处，没有渗流，即： 则式(9)变为 D(0)0=0,x=L4:D(6)0=0,z=L 翌.(]+o]- ak.(0) 式中，L,为计算区域的右边界坐标，L,为计算区域的 下边界坐标. (12) (4)初始条件. 式(12)即为考虑二维降雨入渗的非饱和土边坡 入渗方程.式中：D(o)为水力扩散系数，cm2·min:0 0l4.=8:(x,z）,z>0,l=0 2.3二维降雨入渗的数值解法 为体积含水率；K(0)为土的导水率，cm·minl;t为时 由上述分析可知，二维降雨入渗过程是一个二阶 间，min. 偏微分方程，采用有限差分法进行求解.为保证精度 2二维降雨入渗模型数值求解 和计算效率，采用交替隐式差分法(ADI). 交替隐式差分法的具体做法是：先在时段，k+ 2.1基本假设 1]内，对x方向取中心隐式差分格式，而z方向则取显 对于二维渗流来说，水的运动方程是一个二阶偏 式差分格式，可列出差分方程并对其求解，得到时段末 微分方程，为求解方便，作如下假设： (t=k+1)的体积含水率分布：然后在时段k+1,k+ (1)假设土体为匀质土体，仅考虑扩散系数D和 2]内对：方向取中心差分分格式，而x方向取显示差 导水系数K随含水量的变化： 分，可列出另一个差分方程，得到时段末(t=k+2)的 (2)为研究方便，在初始状态时，假设土体中各处 体积含水率分布.如此交替进行，便可求出体积含水 的体积含水率处处相等： 率随降雨时间的变化 (3)假设地下水位不变，且埋藏很深，在有限的降 (1)x方向隐式及z方向显式.式(12)中各项导 雨时间内，降雨对地下水位影响不大； 数的差分近似式为： (4)降雨影响深度最大为10m,降雨历时最大为 0_- 24h; dt (13) (5)导水率表达式取K(0)=a ,扩散率取 是{po]- D(8)=b（a. ,系数通过室内实验或原位实验 D)-D( ,(14) △r2 求得： (6)不考虑温度、蒸发等影响因素 是po2]= 2.2边界条件 Dn（t1-)-D-n(,-成-，(15) (1)上边界条件.降雨强度为R()且已知，当降 43 雨开始时，土的入渗能力较大，降雨强度小于土的入渗 aK(0)i-K- (16) 能力，土的实际入渗率就是降雨强度，经历了一定的时 2△z 间后，土的入渗能力减小，降雨强度大于土的入渗率， 式中，Di2K等是D(in)、K(d)等的简写
付建新等: 考虑二维降雨入渗的非饱和土边坡瞬态体积含水率分析 率为 θ t dxdzdt． ( 7) 由质量守恒可得 θ t = Vwx x + Vwz z ． ( 8) 将式( 4) 代入式( 8) ，可得单元体二维渗流方程为 θ t = [ x Kx ( θ) ψm ] x + [ z Kz ( θ) ψm ] z － Kz ( θ) z ． ( 9) D( θ) 和 K( θ) 之间的关系为［15］ D( θ) = K( θ) C( θ) ． ( 10) 式中 C( θ) 为比水容重，定义为 C( θ) = dθ dψm ． ( 11) 则式( 9) 变为 θ t = [ x Dx ( θ) θ ] x + [ z Dz ( θ) θ ] z － Kz ( θ) z ． ( 12) 式( 12) 即为考虑二维降雨入渗的非饱和土边坡 入渗方程． 式中: D( θ) 为水力扩散系数，cm2 ·min － 1 ; θ 为体积含水率; K( θ) 为土的导水率，cm·min － 1 ; t 为时 间，min． 2 二维降雨入渗模型数值求解 2. 1 基本假设 对于二维渗流来说，水的运动方程是一个二阶偏 微分方程，为求解方便，作如下假设: ( 1) 假设土体为匀质土体，仅考虑扩散系数 D 和 导水系数 K 随含水量的变化; ( 2) 为研究方便，在初始状态时，假设土体中各处 的体积含水率处处相等; ( 3) 假设地下水位不变，且埋藏很深，在有限的降 雨时间内，降雨对地下水位影响不大; ( 4) 降雨影响深度最大为 10 m，降雨历时最大为 24 h; ( 5) 导水率表达式取 K( θ) = a1 ( θ θ ) s b1 ，扩散率取 D( θ) = b1 ( θ θ ) s b2 ，系数通过室内实验或原位实验 求得; ( 6) 不考虑温度、蒸发等影响因素． 2. 2 边界条件 ( 1) 上边界条件． 降雨强度为 Ｒ( t) 且已知，当降 雨开始时，土的入渗能力较大，降雨强度小于土的入渗 能力，土的实际入渗率就是降雨强度，经历了一定的时 间后，土的入渗能力减小，降雨强度大于土的入渗率， 形成了地表积水或者径流，地面体积含水率可接近饱 和体积含水率 θs 的值并保持不变，此时土体表面处边 界条件为: － D( θ ( ) θ x sinα + θ z cosα ) + K( θ) cosα = Ｒ( t) cosα，z = 0，0 ＜ t ＜ ta ; θ = θs，z = 0，t ＞ ta { ． 式中，ta 为降雨强度超过土的入渗能力的时刻，Ｒ( t) 为 降雨强度，θs 为饱和体积含水率． ( 2) 下边界条件． 假设地下水埋藏较深，土壤深 处含水量恒定且分布均匀，即 θ0 = θi ( x，z) ，z = L，t ＞ 0． 式中，θ0 为初始体积含水率，L 为计算区域的下边界． ( 3) 右边界条件． 为了简便的说明问题，假定远 离边坡边界处，没有渗流，即: D( θ) θ x = 0，x = Lx ; D( θ) θ z = 0，z = Lz ． 式中，Lx 为计算区域的右边界坐标，Lz 为计算区域的 下边界坐标． ( 4) 初始条件． θ | t = t0 = θi ( x，z) ，z ＞ 0，t = 0． 2. 3 二维降雨入渗的数值解法 由上述分析可知，二维降雨入渗过程是一个二阶 偏微分方程，采用有限差分法进行求解． 为保证精度 和计算效率，采用交替隐式差分法( ADI) ． 交替隐式差分法的具体做法是: 先在时段［k，k + 1］内，对 x 方向取中心隐式差分格式，而 z 方向则取显 式差分格式，可列出差分方程并对其求解，得到时段末 ( t = k + 1) 的体积含水率分布; 然后在时段［k + 1，k + 2］内对 z 方向取中心差分分格式，而 x 方向取显示差 分，可列出另一个差分方程，得到时段末( t = k + 2) 的 体积含水率分布． 如此交替进行，便可求出体积含水 率随降雨时间的变化． ( 1) x 方向隐式及 z 方向显式． 式( 12) 中各项导 数的差分近似式为: θ t = θ k + 1 ij － θ k ij Δt ， ( 13) [ x D( θ) θ ] x = Dk + 1 i + 1 /2，j ( θ k + 1 i + 1，j － θ k + 1 i，j ) － Dk + 1 i － 1 /2，j ( θ k + 1 i，j － θ k + 1 i － 1，j ) Δx 2 ，( 14) [ z D( θ) θ ] z = Dk i，j + 1 /2 ( θ k i，j + 1 － θ k i，j ) － Dk i，j － 1 /2 ( θ k i，j － θ k i，j － 1 ) Δz 2 ， ( 15) K( θ) z = Kk i，j + 1 － Kk i，j － 1 2Δz ． ( 16) 式中，Dk + 1 i + 1 /2，j 、Kk i，j + 1等是 D( θ k + 1 i + 1 /2，j ) 、K( θ k i，j + 1 ) 等的简写． · 904 ·
410 工程科学学报，第37卷，第4期 rDn+n(K-K）. (30) 454 7s △t 25=2△ 等式左边涉及的D均为z方向扩散系数的值，等 式右边涉及的D均为x方向扩散系数的值，K为z方 式(12)可变为： 向的值.同理可以得到如下的三对角矩阵 agc=y (17) 「buC 0 61 ag=rDiing hu (18) a12b12 C12 始 =-+n(Dn+D] he (19) =nDn (20) 01,-2 b1,-2 C1,u-2 h与=-rD-ln,-l+2(D-n+Din）-1]- h1,a-2 0 01,m-l b-LJ nDn+r,（K1-K-i. Lhi.-1 (21) 利用“追赶法”解出每个时段末的体积含水率。对 等式左边涉及的D均为x方向扩散系数的值，等 于i=2这种情况，原理是相同的，只是用的体积含水 式右边涉及的D均为z方向扩散系数的值，K为z方 率是上一个节点求出的体积含水率.以此类推，就可 向的值. 以求出时段末z方向各个节点的体积含水率.这样， 假设上边界条件是第一类边界条件，即体积含水 所有节点的时段末体积含水率都已经求出. 率已知，当=1时对x方向的各个节点写出差分方 当降雨强度大于入渗率时，根据原边界方程补写 程.很明显，我们可以将式(17)写成一个三对角矩阵. 差分方程，如式中的第一个公式，在边界结点j=0处 利用“追赶法”解出每个时段末的体积含水率.对 于=2这种情况原理是相同的，只是用的体积含水率 列差分方程，驰取向前差分，可得差分方程如下： z 是上一个节点求出的体积含水率.这样以此类推，就 可以求出时段末x方向各个节点的体积含水率. :"-:-Ra (31) △z Tou cu 0「1 式中，D。和K分别为z方向的参数D()和 a21 K(:),R1P表示本时段的平均供水强度.式(31) 可写成 0e-2,1 ba-2.1 Cm-2,1 ha-2.1 bo0o+cadi=ho (32) Gn-l,l B- Lhe-1.1 式中，bn=-DAz,cn=-bo,hn=-R1n+K. (2):方向隐式及x方向显式.原方程式(12)中 将上式作为方程组的第一个方程即可组成三对角 各项导数的差分近似式为： 矩阵，用同样的办法可以求得各节点的体积含水率. 0- 基于上述算法，用MATLAB编制了计算程序，可 (22) 迅速的计算出土体内瞬态体积含水率分布. [oo]= 3工程案例分析 De(du-）-2（i-2,(23) 3.1基本参数 △x2 利用本文的程序，对某路基边坡在降雨过程中土 [po0]- 体内部瞬态体积含水率的分布情况（路基边坡坡 脚55)进行如下设定： D ()-D ( 42 ,(24) (1)边坡土体中的水分运动参数为 D,=D,=10.4(0/0.)6,K,=K=0.02(0/8.)28 ak(0) -出 (25) (33) dz 2△z (2)初始体积含水率为0.28，饱和体积含水率 式中，D广nK等是D(n）、K()等的简写. 为0.45： 原方程可变为： (3)计算下边界深度取3m,右边界水平距离取距 ab+c=hg (26) 边坡边缘1.0m,计算步长为5cm,时间步长为1min. ay=rDi (27) 3.2不同工况边坡瞬态体积含水率分析 6,=-0+(Dn+Dn)], (28) (1)设初始体积含水率为0.28，保持降雨强度不 c=aD (29) 变，计算不同时间土体的瞬时体积含水率与深度的关 h=-Di+(Di-i+Di)-1]0j- 系，降雨强度分别选取200、300、400和500mmd,每
工程科学学报，第 37 卷，第 4 期 令 r1 = Δt Δx 2，r2 = Δt Δz 2，r3 = Δt 2Δz ． 式( 12) 可变为: aijθ k + 1 i － 1，j + bijθ k + 1 i，j + cijθ k + 1 i + 1，j = hij， ( 17) aij = r1Dk + 1 i － 1 /2，j ， ( 18) bij = －［1 + r1 ( Dk + 1 i + 1 /2，j + Dk + 1 i － 1 /2，j ) ］， ( 19) cij = r1Dk + 1 i + 1 /2，j ， ( 20) hij = － r2Dk i，j － 1 /2 θ k i，j － 1 +［r2 ( Dk i，j － 1 /2 + Dk i，j + 1 /2 ) － 1］θ k i，j － r2Dk i，j + 1 /2 θ k i，j + 1 + r3 ( Kk i，j + 1 － Kk i，j － 1 ) ． ( 21) 等式左边涉及的 D 均为 x 方向扩散系数的值，等 式右边涉及的 D 均为 z 方向扩散系数的值，K 为 z 方 向的值． 假设上边界条件是第一类边界条件，即体积含水 率已知，当 j = 1 时对 x 方向的各个节点写出差分方 程． 很明显，我们可以将式( 17) 写成一个三对角矩阵． 利用“追赶法”解出每个时段末的体积含水率． 对 于 j = 2 这种情况原理是相同的，只是用的体积含水率 是上一个节点求出的体积含水率． 这样以此类推，就 可以求出时段末 x 方向各个节点的体积含水率． b11 c11 0 a21 b21 c21 … … … an － 2，1 bn － 2，1 cn － 2，1 0 an － 1，1 bn － 1， 1 θ k + 1 11 θ k + 1 21 θ k + 1 n － 2，1 θ k + 1 n － 1， 1 = h11 h21 hn － 2，1 hn － 1， 1 ． ( 2) z 方向隐式及 x 方向显式． 原方程式( 12) 中 各项导数的差分近似式为: θ t = θ k + 1 ij － θ k ij Δt ， ( 22) [ x D( θ) θ ] x = Dk i + 1 /2，j ( θ k i + 1，j － θ k i，j ) － Dk i － 1 /2，j ( θ k i，j － θ k i － 1，j ) Δx 2 ， ( 23) [ z D( θ) θ ] z = Dk + 1 i，j + 1 /2 ( θ k + 1 i，j + 1 － θ k + 1 i，j ) － Dk + 1 i，j － 1 /2 ( θ k + 1 i，j － θ k + 1 i，j － 1 ) Δz 2 ，( 24) K( θ) z = Kk + 1 i，j + 1 － Kk + 1 i，j － 1 2Δz ． ( 25) 式中，Dk i + 1 /2，j 、Kk + 1 i，j + 1等是 D( θ k i + 1 /2，j ) 、K( θ k + 1 i，j + 1 ) 等的简写． 原方程可变为: aijθ k + 1 i，j － 1 + bijθ k + 1 i，j + cijθ k + 1 i，j + 1 = hij， ( 26) aij = r2Dk + 1 i，j － 1， ( 27) bij = －［1 + r2 ( Dk + 1 i，j + 1 /2 + Dk + 1 i，j － 1 /2 ) ］， ( 28) cij = r2Dk + 1 i，j + 1 /2， ( 29) hij = － r1Dk i － 1 /2，jθ k i － 1，j +［r1 ( Dk i － 1 /2，j + Dk i + 1 /2，j ) － 1］θ k i，j － r1Dk i + 1 /2，jθ k i + 1，j + r3 ( Kk + 1 i，j + 1 － Kk － 1 i，j － 1 ) ． ( 30) 等式左边涉及的 D 均为 z 方向扩散系数的值，等 式右边涉及的 D 均为 x 方向扩散系数的值，K 为 z 方 向的值． 同理可以得到如下的三对角矩阵． b11 c11 0 a12 b12 c12 … … … a1，n － 2 b1，n － 2 c1，n － 2 0 a1，n － 1 b1，n － 1 θ k + 1 11 θ k + 1 12 θ k + 1 1，n － 2 θ k + 1 1，n － 1 = h11 h12 h1，n － 2 h1，n － 1 ． 利用“追赶法”解出每个时段末的体积含水率． 对 于 i = 2 这种情况，原理是相同的，只是用的体积含水 率是上一个节点求出的体积含水率． 以此类推，就可 以求出时段末 z 方向各个节点的体积含水率． 这样， 所有节点的时段末体积含水率都已经求出． 当降雨强度大于入渗率时，根据原边界方程补写 差分方程，如式中的第一个公式，在边界结点 j = 0 处 列差分方程，θ z 取向前差分，可得差分方程如下: Dk + 1 i，0 θ k + 1 i，1 － θ k + 1 i，0 Δz － Kk + 1 i，0 = － Ｒk + 1 /2 ． ( 31) 式中，Dk + 1 i，0 和 Kk + 1 i，0 分别 为 z 方 向 的 参 数 D( θ k + 1 i，0 ) 和 K( θ k + 1 i，0 ) ，Ｒk + 1 /2 表示本时段的平均供水强度． 式( 31) 可写成 bi0 θ k + 1 i，0 + ci0 θ k + 1 i，1 = hi0 ． ( 32) 式中，bi0 = － Dk + 1 i，0 /Δz，ci0 = － bi0，hi0 = － Ｒk + 1 /2 + Kk + 1 i，0 ． 将上式作为方程组的第一个方程即可组成三对角 矩阵，用同样的办法可以求得各节点的体积含水率． 基于上述算法，用 MATLAB 编制了计算程序，可 迅速的计算出土体内瞬态体积含水率分布． 3 工程案例分析 3. 1 基本参数 利用本文的程序，对某路基边坡在降雨过程中土 体内部 瞬 态 体 积 含 水 率 的 分 布 情 况 ( 路 基 边 坡 坡 脚 55°) 进行如下设定: ( 1) 边坡土体中的水分运动参数为 Dx = Dz = 10. 4( θ /θs) 3. 6，Kx = Kz = 0. 02( θ /θs) 2. 8 ; ( 33) ( 2) 初始体积含水率为 0. 28，饱和体积含水率 为 0. 45; ( 3) 计算下边界深度取 3 m，右边界水平距离取距 边坡边缘 1. 0 m，计算步长为 5 cm，时间步长为 1 min． 3. 2 不同工况边坡瞬态体积含水率分析 ( 1) 设初始体积含水率为 0. 28，保持降雨强度不 变，计算不同时间土体的瞬时体积含水率与深度的关 系，降雨强度分别选取 200、300、400 和 500 mm·d － 1，每 · 014 ·
付建新等：考虑二维降雨入渗的非饱和土边坡瞬态体积含水率分析 411 个降雨强度下分别取120、240、480、720、960和 1440min进行计算 如图3~图6所示，当降雨强度为300mmd时， 0.46 0.44 降雨时间分别为240、480、960及1440mim时，边坡土 0.42 体内不同位置体积含水率分布图 0.36 是0.34 0.32 0.40 0.30 0.38 0.28 0.36 0.34 050100150200250300350050100150200250300350 深度cm 水平距离em 0.32 未影响区域 0.30 图6降雨时间为1440min时土体体积含水率分布 0.2 00 Fig.6 Distribution of volume water content after 1440 min rainfall 050n00 响区域越来越小，影响深度从1.0m增大到2.5m左 深度em 50300 500 50100150200250300350 水平距离/cm 右，水平影响距离从0.5m增大到1.5m.另外，随着深 度和水平距离增加，曲面上网格逐渐变小，说明降雨对 图3降雨时间为240mim时土体体积含水率分布 Fig.3 Distribution of volume water content after 240 min rainfall 土体体积含水率的影响也越小，由于在边坡表面考虑 了水平方向的入渗，因此表面的体积含水率变化比较 大.对于其他降雨强度下的体积含水率分布规律是相 0.46 同的，这里不再赘述 0.44 图7和图8为降雨时间分别为960min和1440min 0.42 30.40 时边坡土体的体积含水率的等值线分布图.最底部的 0.38 如0.36 一条是体积含水率为0.30的等值线，每条等值线之间 0.34 相隔0.02. 0.32 0.30 0.38 0 0286 0.36 -0.34 0.32 50 1001504h0050300 深度cm 50050100150200250300350 -0.30 水平距离cm 100 图4降雨时间为480min时土体体积含水率分布 150 Fig.4 Distribution of volume water content after 480 min rainfall 200 250 0.46 0 50100150 200 250 300cm 044 0.42 040 图7降雨时间为960min时体积含水率等值线分布图 038 0.36 Fig.7 Contour maps of volume water content after960 min rainfall 0.34 0.32 020 0.28 从各个降雨时间土坡内部的体积含水率等值线可以 50 00 看出：随着降雨时间变长，等值线在不断的下移，并且变得 15 深度em 200 350050100150200230300350 分散，说明降雨影响的深度在不断加深：在斜坡段部分，等 25 3 水平距离/em 值线是向下倾斜，说明了表层的体积含水率大于土体内 部的体积含水率：随着水平距离的变大，等值线变为趋于 图5降雨时间为960min时土体体积含水率分布 平行的直线，说明水平入渗的影响在逐渐变小 Fig.5 Distribution of volume water content after 960 min rainfall (2)边坡初始体积含水率对土体内部瞬态体积含 由图可知，体积含水率是随着深度的增加而逐渐 水率也有很大的影响.选取降雨强度为450mm·d', 减小，逐渐接近土体的初始体积含水率，随着降雨时间 初始体积含水率分别选0.15、0.25、0.30及0.35，降雨 的延长，影响深度及水平距离不断增加，即图3中未影 时间选取480min,距边坡边缘1.0m处作垂直剖面，研
付建新等: 考虑二维降雨入渗的非饱和土边坡瞬态体积含水率分析 个降 雨 强 度 下 分 别 取 120、240、480、720、960 和 1440 min 进行计算． 如图 3 ～ 图 6 所示，当降雨强度为 300 mm·d － 1时， 降雨时间分别为 240、480、960 及 1440 min 时，边坡土 体内不同位置体积含水率分布图． 图 3 降雨时间为 240 min 时土体体积含水率分布 Fig． 3 Distribution of volume water content after 240 min rainfall 图 4 降雨时间为 480 min 时土体体积含水率分布 Fig． 4 Distribution of volume water content after 480 min rainfall 图 5 降雨时间为 960 min 时土体体积含水率分布 Fig． 5 Distribution of volume water content after 960 min rainfall 由图可知，体积含水率是随着深度的增加而逐渐 减小，逐渐接近土体的初始体积含水率，随着降雨时间 的延长，影响深度及水平距离不断增加，即图 3 中未影 图 6 降雨时间为 1440 min 时土体体积含水率分布 Fig． 6 Distribution of volume water content after 1440 min rainfall 响区域越来越小，影响深度从 1. 0 m 增大到 2. 5 m 左 右，水平影响距离从 0. 5 m 增大到 1. 5 m． 另外，随着深 度和水平距离增加，曲面上网格逐渐变小，说明降雨对 土体体积含水率的影响也越小，由于在边坡表面考虑 了水平方向的入渗，因此表面的体积含水率变化比较 大． 对于其他降雨强度下的体积含水率分布规律是相 同的，这里不再赘述． 图7 和图8 为降雨时间分别为960 min 和1440 min 时边坡土体的体积含水率的等值线分布图． 最底部的 一条是体积含水率为 0. 30 的等值线，每条等值线之间 相隔 0. 02． 图 7 降雨时间为 960 min 时体积含水率等值线分布图 Fig． 7 Contour maps of volume water content after 960 min rainfall 从各个降雨时间土坡内部的体积含水率等值线可以 看出: 随着降雨时间变长，等值线在不断的下移，并且变得 分散，说明降雨影响的深度在不断加深; 在斜坡段部分，等 值线是向下倾斜，说明了表层的体积含水率大于土体内 部的体积含水率; 随着水平距离的变大，等值线变为趋于 平行的直线，说明水平入渗的影响在逐渐变小． ( 2) 边坡初始体积含水率对土体内部瞬态体积含 水率也有很大的影响． 选取降雨强度为 450 mm·d － 1， 初始体积含水率分别选 0. 15、0. 25、0. 30 及 0. 35，降雨 时间选取 480 min，距边坡边缘 1. 0 m 处作垂直剖面，研 · 114 ·