基于双重逾渗模型的裂隙多孔介质连通性研究

基于双重逾渗模型的裂隙多孔介质连通性研究 DOI :10. 3880/ j . issn. 1006Ο7647. 2012. 02. 006 基于水2气二相流模型的土坡稳定性分析 1 2 3 4张晓悦,王 栋,张晓乐,沈跃军 (1. 浙江水利水电专科学校水利工程系 ,浙江 杭州 310018 ; 2. 信息产业电子第十一设计研究院科技工程股份有限公司杭州分院 ,浙江 杭州 310020 ; ) 3. 浙江华东工程安全技术有限公司 ,浙江 杭州 310014 ; 4. 浙江省钱塘江管理局勘测设计院 ,浙江 杭州 310016摘要 :为了研究孔隙气压力和负孔隙水压力在土坡稳定分析中所起的作用 ,基于多相流理论建立水 2气二相流模型 ,利用此模型对稳定渗流情况和降雨情况下的土体边坡内水相和气相的渗流状态进 行模拟 ,根据模拟得到的孔隙气压力和孔隙水压力值求出土坡安全系数 。计算结果表明 ,土坡非饱 和区负孔隙水压力的存在较大地提高了土坡稳定性 ,在稳定渗流情况下 ,孔隙气压力对土坡稳定的影响可忽略 ,在降雨情况下 ,非饱和区产生的孔隙气压力使得安全系数减小 ,滑动面与地下水位的 距离越大 ,非饱和区的孔压对土坡稳定的影响越大 。 关键词 :非饱和土 ;二相流 ;相态 ;降雨 ;土坡稳定 + () 文献标识码 :A中图分类号 : TU441 . 35 文章编号 :1006Ο7647 201202Ο0023Ο05 1 2 Sta bility analysis of soil slope based on water2a ir t wo2phase flo w model/ / ZHANG Xiao2yue, WANG Dong, ZHANG Xiao2 3 4 ( le, SHEN Yue2jun1. Department of Hydraulic Engineering , Zhejiang Water Conservancy and Hydropower College , Hangzhou 310018 , China ; 2. Hangzhou B ranch of the IT Electronics Eleventh Design & Research Institute , Scientif ic and Technological Engineering Co. , Ltd. , Hangzhou 310020 , China ; 3. Zhejiang Huadong S afety Technology Engineering Limited Company , Hangzhou 310014 , China ; 4. S urveying and Design Institute of Qiantang River Administration of Zhejiang Province , Hangzhou )310016 , China Abstract : Based on the multi2phase flow theory , a water2air two2phase flow model is established to study the effects of the pore air pressure and negative pore water pressure in the soil slope stability analysis. Using the model , the water flow and air flow in soil slope under steady seepage and rainfall conditions are simulated , and safety coefficients of the soil slope are calculated. The results show that the soil slope stability can evidently increase under the effect of the negative pore water pressure of the unsaturated zone . The pore air pressure in the unsaturated zone can be neglected under steady seepage condition , which should be taken into consideration under rainfall condition as it can reduce the safety coefficient . Influence of the pore pressure of unsaturated zone on soil slope stability will increase with the distance between the slip surface and the underground water level . Key words : unsaturated soil ; two2phase flow ; phase state ; rainfall ; slope stability 3 型的一种空间和时间离散方法 ,White 等分析了多 土坡稳定分析是估计天然及人工边坡安全系数 4 孔介质饱和度和相对渗透系数关系 ,Laroche 等对 ,一般进行边坡稳定分析时 ,由于测 的一个常用工具 量负孔隙水压力和孔隙气压力并将其纳入稳定分析 多相流模型中毛细压力和饱和度的关系做了研究 。 中比较困难 ,对于地下水位以上非饱和区中由负孔 本文建立水2气二相流模型来模拟土体边坡在 稳定渗流情况和降雨情况下的水相和气相渗流场 , 隙水压力和孔隙气压力所提供的部分抗剪强度通常 予以忽略不计 。近年来 ,随着人们对多相流的研究 全面考虑土体中的水相和气相流动以及水相和气相 越来越多 ,在土坡稳定分析中全面考虑孔隙水压力 成分的相互转换 ,利用二相流计算结果求得边坡安 1 和孔隙气压力成为可能 。例如 ,Dijke 等给出了一 全系数 ,分析孔隙气压力和负孔隙水压力在土坡稳 2 定中的作用 。 个较完整的多相流模型 ,Matthew 等提出二相流模 ( ) ()基金项目 :科技部中奥两国科技合作项目 CN 01/ 2007;国家自然科学基金 50539020 () 作者简介 :张晓悦 1984 —,女 ,山西忻州人 ,讲师 ,博士 ,主要从事水工结构研究 。E2mail :zhangxy @zjwchc . com 阶项 , 在第 m + 1 次迭代步 、k + 1 次时间步 : 1 水2气二相流模型 ( ) ( )F x?+ F x k +1 , m +1k +1 , m 1. 1 定义以及基本假设 5 F )( ( - x ) 4 x k +1 , m +1 k +1 , m 5 x k +1 , m 非饱和土是四相混合体 ,包括固相 、气相 、液相 5 ( ) () 令 F x= 0 , 将式 4改写为 k + 1 , m + 1以及称为收缩膜的水气分界面。在土体中可能有 ) ( )( )( 5 B xuF x= - k +1 , m k +1 , m 多种液相 ,本文只考虑水相 。非饱和土中的渗流为 式中 : B 为系数矩阵 , B = 5 F/ 5 x ; u 为主要变量增 包括气相和液相流动的二相流 ,气相和液相都包括 量 , u= x-) ( x; F x为第 k + 1 时 j k + 1 , m + 1 k + 1 , m k + 1 , m 水和空气两种组分 。液相和气相通过蒸发与溶解相 5 间步 、m 迭代步的误差值 。系数矩阵 B 用数值差分 互转换,相态及其相互转换关系见图 1 。 法来求 , 系数 b值为ij 5 Fi , k +1 , m ( Δb= ) ? [ Fx, , x, x+x, x, - ij i 1 j - 1 j j j +1 5 xj , k +1 , m ) Δ( Δ/ 2x Fx, , x, x- x, x, i 1 j - 1 j j j +1 j )( (), n j = 1 , 2 , 6 图 1 相态及相态转换 Δ:x为一微小增量 ; n 为主要变量数乘以单元式中 j 模型假设土体骨架不变形 ,即孔隙率为常数 ,孔 数的值 。隙流体不可压缩 ,忽略化学生物学反应 ,恒温 ,不考 Newton2Raphson 法具体计算步骤如下 : 第 k + 1虑本构关系中渗透系数与饱和度之间的滞后作用 , 气相中所有组分都遵守理想气体定律 ,气相中水组 时 间 步 , 初 始 值 x为 上 一 时 间 步 的 解 , 当 k + 1 , 0 分分压力等于此温度下的饱和蒸汽压力 。 ( ( ) ) ‖F xεB x ‖ >时 , 求 解 u = k + 1 , m 2k + 1 , m [8 ] 1. 2 模型控制方程及其求解) η(η( , 令 x= x + u < 1 ) - F x , k + 1 , m k + 1 , m + 1 k + 1 , m j 模型考虑气相和液相的流动以及各组分之间的质( ) 为迭代的步长系数 减幅系数, 进行下一迭代步计 ) ε( 算 , 直到 ‖F x‖<, 进行下一时间步计算 。 k + 1 , m 2 6κ量传递 ,模型的控制方程即组分的质量守恒方程 : 1. 3 主要变量及相态转换 ρ5 Sx在给定的控制体中 , 相态并非固定不变 , 对应不 αα ακm ,,? α= g ,w 同相态需采用不同的主要变量 , 在计算中单元的相 - < 5 t 态和主要变量需 在 每 一 迭 代 步 后 重 新 调 整 。根 据 [ 6 ] ρkxαα ακ rm ,,Gibbsian 相态 定 律, 一 个 多 组 分 多 相 流 系 统 中 自 (ρ ) div Kgrad p- g - Q= 0 α ακ ? μ αα= g ,w ( ) 由度的数量为 F = C - P + 2 + P - 1= C + 1 , 其中 ()C 为组 分 数 , P 为 相 数 , 本 模 型 中 , 由 于 温 度 为 常 1 量 , 只需另外选择两个独立的主要变量 。只有气相 土体固有渗透系数 : 时 , x和 p为主要变量 ; 只有液相时 , x和 p为 g ,wg w ,gg μk sww( ) 2 K = 主要变量 ; 水相气相同时存在时 , S和 p为主要变 w g ρg w 量 。相态转换的执行程序见图 2 。 ακ式中 :下标 和分别代表相态和组分 , 取 g 时为气 相或空气组分 , 取 w 时为水相或水组分 ; < 为土体孔 ρα α 隙率 ; S为相的饱和度 ;= n/ V, n为相 α α αα α m ,3 ρα 的摩尔质量 , V为 n对应的体积 , m;为相的 α α α 密 ακ度 ; x为相中组分的摩尔分数 ; t 为时间 ; k 为土α,κsw α为pα 的饱和渗透系数 ; k 为相的相对渗透系数 ; α rα μα相的压力 ;为相的绝对黏滞性 , 温度一 定的情况α κ 下为常数 ; g 为重力加速度 ; Q为组分 κ 的源汇项 , kg/ s 。 控制方程的变量分为主要变量和次要变量两类 , 次要变量为主要变量的函数 。采用有限差分法对控 [7 ] 制方程进行空间和时间的离散, 用 Newton2Raphson [6 ]() 法求解模型控制方程 。将式 1简化为函数形式 : ( ) ()F x= 0 3 图 2 相态转换执行程序 ( ) 式中 x 为主要变量 。泰勒级数展开式 3, 忽略高 1. 4 次要变量条件单元赋初值实现 ; Neumann 边界条件 , 即流量边 11411 摩尔分数和密度界条件 , 包括水相流量和气相流量边界 , 给边界条件 ( ( ) ) ακ单元施加一个对应于组分质量守恒方程 式 1的 相组分的摩尔分数 x满足 :ακ, 源汇项 。在稳定渗流情况下源汇项为 0 , 在降雨情 ( )x= 17 ακ ,? [5 ]κ= w ,g 况下 , 水相源汇项计算公式为其中气相中水组分的摩尔分数 : ρ()Q = 15 A q w w t()8 xp/ p= g ,w w ,sat g式中 : A 为与降雨方向垂直的土体表层有效面积 ; q t 式中 p为此温度下的饱和蒸汽压力 。 w ,sat为 t 时刻的降雨强度 。 根据 Henry 定律和 Dalton 分压定律 , 得到 水 相 [9 ] 2 土坡稳定分析公式 中空气组分的摩尔分数: ( ))(x9 = Kp- p 采用圆弧滑动面 ,根据简化 Bishop 法 ,假设 X- H g w ,satw ,gR X= 0 ,其中 X与 X分别为土条右侧和左侧的条 L R L ( 式中 : K为 Henry 系 数 , K= 01894 2 + 1147exp H H [11 ] 间竖向剪力, 已知水 2气二相流模型计算出的单 - 10 - 1 ( ) ) - 01043 94 T×10 Pa ; T 为绝对温度 。 忽略元孔隙压力值 , 要将该计算结果应用于边坡稳定分 ρ压力影响 , 在常温下 , 水相密度 取常数 w 析中 , 就要根据水 2气二相流模型网格节点的孔隙压 3 1 000 kg/ m。由理想气体定律可得力插值得到滑动面上各点的孔隙压力值 , 可采用等 ()ρ= p/ R T 10 gm ,g [12 ] 参单元的有限元逆变换来实现 , 本文采用 Taylor ρρ( )= xM11 κκgm ,g g ,? κ= w ,g 级数展开的线性项进行逆变换 , 插值求得各土条底 κ式中 : R 为气体常数 ; M为组分的分子质量 。 κ 边中点的孔隙压力值 , 然后根据力矩平衡来计算边 11412 毛细压力和相对渗透系数 收缩膜两侧承受 坡安全系数 。 的水压力和空气压力之差称为 根据竖向力平衡 , 作用于第 i 土条底面上的总 毛细压力 , 即 p= p- p?0 Pa , 毛细压力和相对渗 c w g 法向力 N 计算公式为 i 透系数都为饱和度的函数 , 采用 Genuchten 提出的近 βαβα csinsin 1i ii iφ[10 ] N W-+ ptan+= ii w , ib 似公式: mF F α, i ξ λ- 1/ 1/ ( ) ()( ) 12 pS= PS - 1 c wβ α 0 e sin i i (φ φ))(ptan- tan16 g , ib/ 2 1ξξ/21F ( )() () = S [ 1 - 1 - S ] 13k S erw we ξξαφ1/ 3 /2 sintan 1i( )( ) ( )14 = 1 - S [ 1 - S ]kS e erg wαm = co s 其中 + α, i i F ξλλξ式中 :,为经验参数 , 取 = 1 - 1/ ; P为土的进 0 β式中 : φ 为 i 土条底面斜向长度 ; c 为有效黏聚力 ; i ( ) ( ) 气值 ; S 为有效饱和度 , S = S- S / S - S , e e w rwsw rwφ为有效内摩擦角 , 而参数 代表抗剪 强 度 随 毛 细 b 其中 S 为剩余水饱和度 , S 为最大水饱和度 。 rwswφ, 当土条底面位于饱 和 区 内 时 ,压力的增加情况 b 1. 5 定解条件αφ 采用 值 ;为 i 土条底面中点的切线与水平面的 夹i a . 边界条件的实现 :在渗流场边界外围加一圈 角 ; F 为安全系数 ; W为 i 土条重 , 浸润线以下按 i 薄层单元 , 称之为边界条件单元 , 通过赋初值给边界 饱和密度计算 ; p, p分别为 i 土条底面中点的孔 w , i g , i 条件单元来设定渗流场的边界条件 , 假设其体积为 隙水压力和孔隙气压力 , 为相对压力 。 50 3 一个极大值 , 如 10m, 以保证在计算过程中边界条 由力矩平衡可得安全系数为 件不发生改变 。 φ1 tan bβc+β F = i N - p - i w , i i ? φtan b. 稳定渗流计算的初始条件 :假设边坡内部土αi Wsini i? i 体单元水相饱和 , 主要变 量 为x和p, 令x= w ,g g w ,g φtan b- 10 φ()βtan 17 p1 - g , i i10 , p等于大气压 ; 上下游水位以上的边界 g φtan 条件单元气相饱和 , 主要变量为 x和 p, 令 x=g ,gg g ,g () ( ) 联合式 16和 17, 迭代求解安全系数 , 若不考 01999 9 , p等于标准大气压 ; 上下游水位以下的边界 g 虑土体的负孔隙水压力 p和负孔隙气压力 p, 则 w , i g , i 条件单元水相饱和 , 主要变量为 x和 p, 令 x= w ,g g w ,g () ( ) 式 16和式 17变为边坡稳定性分析的饱和土计算 - 10 10 , p等于该边界单元形心点处的水压力值 。降 公式 。g 雨情况的初始条件为稳定渗流计算的结果 。 3 土坡稳定分析算例 c . 两种边界条件 :Dirichlet 边界条件 , 即已知气 取一均质土坡进行计算 ,左侧水位 11 m ,右侧水 压力边界条件和已知水压力边界条件 , 通过给边界 位 16 m ,坡度 27?,底边不透水 ,边坡截面见图 3 。取两 个滑动面进行研究 ,滑动面 1 为深层滑动面 ,大多数 () 土条底面位于饱和区 ,滑动圆心为 16100 m ,25100 m, 滑动半径为 16116 m ; 滑动面 2 为浅层滑动面 ,大多 ( 数土 条 底 面 位 于 非 饱 和 区 , 滑 动 圆 心 为 20100 m , ) 27100 m,滑动半径为 13174 m 。 ()图 3 边坡剖面 单位 :m 边坡土体材料参数 : 饱和渗透系数 k = 110 ×sw - 6ξ10 m/ s , 经验参数 = 01457 , 进气值 P= 20 kPa , 孔 0 隙率 < = 01421 ; 剩余水饱和度 S = 0115 , 最大水饱 rw 和度 S = 110 , 有效黏聚力 c = 1011 kPa , 有效内摩 sw φ φ擦角 = 4216,? 参数 = 3510,? 温度 T = 15 ?, 土体 密度b 3 3 为 1 930 kg/ m, 土体饱和密度为2 000 kg/ m。 图 5 降雨情况水 2气二相流模型计算结果 利用水 2气二相流模型进行模拟 , 得到土坡在稳 润线以上的非饱和区 , 孔隙气压力基本为 0 Pa , 存在 定渗流情况及降雨强度为 712 mm/ h 、降雨时间为 5 h 负孔隙水压力 ; 降雨后 , 土体表层形成暂态饱和区 , 的降雨结束瞬时的孔压及饱和度如图 4 和图 5 所示 非饱和区的水相饱和度增加 , 使得负孔隙水压力减 ( 小 , 由于雨水入渗 , 孔隙中空气尚未全部消散 , 空气 图中孔隙气压力和孔隙水压力分别用相应的压力 所产生的阻力阻碍水相渗入 , 造成孔隙气压力明显 ( ) ρ( ) ρp- p/ g 和 p- p/ g 表示 , 其中水头 增大 , 在浸润线附 近 , 孔 隙 气 压 力 形 成 一 条 基 本 为 p g 0 w w 0 w 0 ) 为大气压力 。图 6 和图 7 同。 0 Pa 的过渡带 。 [13 ] 仅考虑水相流动用有限元法计算, 得到土坡 在稳定渗流和降雨后的孔隙水压力等值线如图 6 和 图 7 所示 , 所以通过对比水 2气二相流和单相流计算 的结果 , 可验证水 2气二相流模型的正确性 。 ( )图 6 单相稳定渗流情况孔隙水压力 m单位 : ( )图 7 单相降雨情况孔隙水压力 m单位 : 根据水 2气 二 相 流 模 型 计 算 出 的 孔 压 值 , 由 式 图 4 稳定渗流情况水 2气二相流模型计算结果 () ( ) 16和式 17求出滑动面 1 和 2 在以下 4 种情况下 由图 4 和图 5 可知 , 边坡浸润线以下的饱和区 ,: ?仅考虑 p > 0 Pa 部分 , 即按饱和土公的安全系数 w 孔隙气压力等于孔隙水压力 。稳定渗流情况下 , 浸式计算 ; ?同时考虑正孔隙水压力 p > 0 Pa 和孔隙 w 气压力 p; ?同时考虑正负孔隙水压力 p> 0 Pa 和 g w 参考文献 : p< 0 Pa ; ?同时考虑正负孔隙水压力 p及孔隙气w w 1 van DIJ KE M I J ,van der ZEE S E A T M ,van DUIJN C J . 压力 p。由 水 2气 二 相 流 模 型 求 得 的 安 全 系 数 见 g Multi2phase flow modeling of air sparging J . Advances in 表 1 。情况 2 和 4 与情况 1 相比较 , 安全系数的增加 () Water Resources ,1995 ,18 6:3192333. 百分比见表 2 。 2 MATTHEW W F , CHRISTOPHER E K , CASS T M. Mixed 表 1 水 2气二相流模型安全系数计算结果 finite element methods and higher order temporal approximations for variably saturated groundwater flow J . 工 况 滑动面 情况 1 情况 2 情况 3 情况 4 Advances in Water Resources ,2003 ,26 :3732394. 1 11997 7 1997 5 1172 7 1172 6 122稳定渗流 2 21981 3 21981 2 41028 7 41028 5 3 WHITE M D , OOSTROM M , IMHARD R J . Modeling fluid 11991 4 11985 0 21118 0 21111 6 1 flow and transport in variably saturated porous media with the 降 雨 2 21973 9 21926 3 31722 0 31674 4 STOMP simulator J . Advances in Water Resources ,1995 ,18 () 6:3532364. % 表 2 不同计算情况下安全系数增加百分比 4 LAROCHE C ,VIZIKA O. Two2phase flow properties prediction 工 况 滑动面 情况 2 与情况 1 相比 情况 2 与情况 4 相比 from small2scale data using pore2network modeling J . - 01010 1755 1 8Transport in Porous Media ,2005 ,61 :77291. 稳定渗流 - 01003 351126 2 5 FREDLUND D G , AHARDIO H. Soil mechanics for - 01321 61036 1 降 雨 unsaturated soilsM . New York :John & Sons ,1993. - 11601 231555 2 6 CLASS H , HELMIG R ,BASTIAN P. Numerical simulation of non2isothermal multiphase multicomponent processes in porous 由表 1 可知 , 考虑非饱和区的孔隙气压力和负 ( ) mediaJ . Advances in Water Resources , 2002 , 25 5: 5332 ( ) 孔隙水压力 情况 4之后 , 边坡安全系数比按饱和 550. ( ) 土公式计算 情况 1的结果有了明显增加 ; 稳定渗 7 李强 ,封建湖 ,蔡体敏 ,等 . 二相流计算的一种差分算法 流情况下 , 孔隙气压力对安全系数的影响可以忽略 , () J . 应用数学和力学 ,2004 ,25 5:4882496.主要是负孔隙水压力对安全系数的增加起作用 , 降 8 韩林 ,张子明 ,倪志强 . 应用改进算法的对称逐步超松 雨后 , 孔隙气压力对安全系数的影响变大 , 孔隙气压 弛预处理共轭梯度法进行大体积混凝土仿真计算 J . 河力起减小安全系数的作用 , 但与负孔隙水压力的作 () 海大学学报 :自然科学版 ,2010 ,38 3:2782283. 用相 比 , 孔 隙 气 压 力 对 稳 定 的 影 响 依 然 较 小 。由 9 KOLDITZ O , de JONGE J . Non 2isothermal two2phase flow in low2permeable porous media J . Computational Mechanics , 表 2 可知 , 滑动面与地下水位的距离越大 , 负孔隙水 2004 ,33 :3452364. 压力和孔隙气压力对边坡稳定的影响就越大 。 10 van GENUCHTEN M T. A closed2form equation for predicting 4 结 论 the hydraulic conductivity of unsaturated soilsJ . Soil Science Society of America Journal ,1980 ,44 :8922898. a . 在稳定渗流情况下 ,非饱和区孔隙气压力基 11 李凯 ,陈国荣 . 基于滑移线场理论的边坡稳定性有限元 本为 0 Pa ,负孔隙水压力不为 0 Pa ,安全系数的计算 ( ) 分析 J . 河海大学学报 : 自然科学版 ,2010 ,38 2: 1912 公式中孔隙气压力的作用可以忽略 ;降雨情况下 ,非 195. 饱和区的负孔隙水压力减小 ,由于空气阻力的作用 , 12 徐燕萍 ,项阳 ,刘泉声 ,等 . 等参元逆变换插值法的研究 孔隙气压力增加 ,此时负孔隙水压力和孔隙气压力 () 及其应用 J . 岩土力学 ,2001 ,22 2:2262228 都对边坡的稳定有影响 。 13 黄耀英 ,吴中如 ,顾璇 . 基于区间参数摄动法的黏弹性b. 在考虑了非饱和区孔隙气压力与负孔隙水 区间有限元研究 J . 河海大学学报 :自然科学版 ,2008 ,压力提供的抗剪强度之后 ,边坡安全系数比按饱和 () 36 5:7022706. 土公式计算的结果有了较大提高 ,其中孔隙气压力 ( ) 收稿日期 :2011 Ο12 Ο06 编辑 :熊水斌使安全系数减小 ,负孔隙水压力使安全系数增加 ,滑 动面与地下水位距离越大 ,则负孔隙水压力与孔隙 气压力对边坡稳定的影响越明显 。 c . 在边坡稳定分析时 ,往往忽略孔隙气压力和 负孔隙水压力的作用 ,由本文分析可知负孔隙水压 力在提高土坡稳定方面的作用不可忽视 ,而在降雨情 况下孔隙气压力对土坡稳定的影响也应引起注意 。