气液两相流原理在边坡稳定性分析中的有效性研究

2021-12-23董顺德孙芝玲何玉鹏马文礼

河北工业科技 2021年1期

董顺德,孙芝玲,何玉鹏,马文礼

(1.青海省环境地质勘查局，青海西宁 810001；2. 青海省环境地质重点实验室，青海西宁 810001；3.青海省水利水电科学研究院有限公司，青海西宁 810001)

近年来，极限平衡等分析方法已被广泛应用于雨水渗透导致边坡崩塌的机理研究中[1-4]。降雨使得土壤饱和度增加，致使坡体质量增加，从而导致边坡滑动力增大，出现边坡失稳现象[5]。另一方面，由于坡体中水压上升和基质吸力减小，导致边坡抗剪强度下降，促使滑坡发生。因此对雨水渗流行为和孔隙水压变化的研究是评估边坡稳定性的必要条件。雨水一般通过地表的非饱和区进入边坡中，而非饱和土的渗透系数与其体积含水率或基质吸力之间的关系是非线性的，应用单纯的数值解析法对其研究较为困难。有学者利用有限差分法(FLAC3D)将饱和土流-固耦合计算原理扩展到非饱和土中[6]。还有学者采用有限元法进行了饱和-非饱和渗流分析[7-8]。有限元法主要是对渗透到坡体的雨水流量进行估算，是计算降雨期间孔隙水压力随时间变化的较为有效的方法[9]。

在对降雨边坡稳定性研究中，缺乏对孔隙气压造成影响的分析。针对降雨对非饱和土影响的问题，有学者研究孔隙气压对边坡造成破坏的影响。基础实验证实，在土中发生空气堵塞时，可通过气压的增加导致地表出现局部破坏的现象观察气压变动与地表破坏的关系[10-11]；CHO[12]研究表明，雨水的渗透和水位的上升会引起气压上升和压裂现象，堤防内部的压缩空气会对堤防整体造成一定程度的损伤；SAITO等[13]针对孔隙气压对雨水渗透的影响进行了室内试验研究，分别对两相流和单相流分析法进行了研究，分析此类方法对雨水渗透过程的适用性。室内试验的结果证实了通过气液两相流分析可以重现孔隙气压的上升和排水量增加的现象，但其室内模型试验忽略了对饱和度、基质吸力和渗透系数之间关系的研究。

笔者对气液两相流理论在降雨边坡稳定性分析中的应用进行了研究。首先，利用圆柱形容器进行一维渗流模型实验，采用不考虑孔隙气压影响的单相流分析和考虑孔隙气压影响的气液两相流分析，分别对一维渗流模型实验进行了分析。对比结果证实了两相流分析法的适用性与合理性。针对均匀坡度的风化岩粉质黏土边坡进行了人工降雨实验，通过气液两相流理论对边坡降雨渗透行为进行计算分析。采用气液两相流法进行应力-渗透耦合计算，模拟2012年青海省黄南州同仁县隆务镇西山发生的滑坡，分析降雨诱发边坡失稳的原因。

1 应力-渗透流耦合分析概述

1.1 有限差分法的特征

有限差分法和有限元法是典型的将目标视为连续体的数值分析方法。由于有限差分法能利用结构网格的拓扑优势轻松扩大模板，构造出高精度格式，在基于显式的应力-变形的时间增量计算中，它比有限元法更具优势。本研究使用了有限差分法解析软件FLAC，此软件具有分析多孔介质渗透的功能，可以与应力-变形分开，从而进行独立的渗流分析，也可以进行应力-渗流耦合分析。

1.2 单相流分析的基本方程式

进行单相流分析时，不考虑孔隙气压的影响，仅考虑水的渗透条件，符合达西定律：

(1)

其中：qi为速度矢量；kij为饱和渗透系数；P为水压；gk为重力加速度；ρf为水的密度；κ为比渗透系数(饱和渗透系数与饱和渗透系数的比),用饱和度s表示：

κ=s2(3-2s)。

(2)

流体的平衡方程如式(3)所示：

(3)

其中:qv为水的流入和流出量大小，而分析中，水被视为不可压缩状态物质。

1.3 两相流分析的基本方程

两相流是指气相和液相流体的流动，2种流体流动都遵循达西定律，气相和液相相关的达西定律分别如式(4)和式(5)所示：

(4)

(5)

其中：μ为黏性系数；P为孔隙压力；w表示液体；a表示气体。式(6)和式(7)中的比渗透系数κr可以用Van Genuchten模型来表示[14]。

(6)

(7)

其中：a，b，c均为常数；Se为有效饱和度。在本研究中，主要采用基于Van Genuchten模型的水特性曲线：

(8)

流体的平衡方程可表示为

(9)

(10)

式中：Sw为饱和度；ζ为单位体积内流体的体积变化；qv为流入和流出量。流体的本构方程如下：

(11)

(12)

其中:Kw和Kg为体积弹性系数。此外在本研究中，不考虑空气在水中的溶解作用。

2 一维渗流模型实验的再现分析

2.1 实验方法

一维渗流模型的实验装置如图1所示。

图1 一维渗流模型实验装置Fig.1 Experiment device of 1D seepage model

试样是直径为280 mm，高为500 mm的圆柱体，孔隙压力表和张力计分别测量其中心和底部的正、负孔隙压力。采用含水比为10%左右的粉质黏土，土体颗粒密度为2.74 g/cm3，孔隙率为0.4。通过压实达到预定的湿密度制备模型试样。为了得到规定的湿密度，对每层分30次压实，分成8层构建模型。模型完成后，整个模型的初始饱和度将达到预计的60%左右，之后进行变水头实验，得出的渗透系数k=1.8×10-5m/s。

调整降雨强度，并从试样顶部供给。为了便于水从顶部均匀地供给，在试样顶部设置了6 cm的沙层。在砂层与试样的边界处设置了4个直径为2 mm的排水孔，使无法渗入粉质黏土的水从砂层中流出，防止积水对试样造成过大压力。

2.2 数值解析

一维渗流解析数值模型如图2所示。

图2 一维渗流解析数值模型Fig.2 Analytical numerical model of 1D seepage

为了再现一维渗流模型实验，将模型尺寸设定为宽280 mm，高500 mm，渗透系数设定为k=1.8×10-5m/s。中部(No.1)和下部(No.2)的设计与模型实验中的孔隙水压测量点相同，并计算中部和下部的孔隙水压。设定降雨强度为50 mm/h，在底面为不排水条件下进行了一维渗流模型实验和数值解析。在气液两相流分析中，需要确定土壤水分参数，即须确定式(6)-式(8)中的a，b，c和P0值。设定b=1/2，c=1/3[15]，对于a和P0的取值，分别用张力计测量在恒定干密度压实条件下不同饱和度试样的基质吸力，确定饱和度与基质吸力的关系。随后，通过改变式(10)中a和P0的数值，与实验结果进行拟合，拟合结果显示a=0.47，P0=8.5 kPa。

2.3 模型实验和数值模拟结果

图3分别为模型实验、单相流分析和两相流分析中孔隙压力的变化。从实验结果看，模型中心No.1测点的孔隙压力首先开始增大，其次底部测量点No.2测点的孔隙压力增大。这表明，渗透是从顶部向底部进行的。张力计测得的孔隙水负压会很快归零，说明在非饱和状态下，孔隙压力表测量得到的是正孔隙压力。导致这个结果的原因是处于相同深度的孔隙压力表和张力计的受压部分没有设置在相同位置。开始通水后约250 min，No.2测点的孔隙水压力升高，No.1测点的孔隙压力在稍作延迟后也有所增加。这是因为前面的渗透面到达模型底部后，下部饱和，水位从模型底部上升。之后2个测量点在整个模型的饱和状态下，静水压力趋于一致。

图3 模型实验和数值模拟结果Fig.3 Model experiment and numerical simulation results

从单相流分析结果可以看出，孔隙水压先从No.2测点开始增大，出现了孔隙水压力即说明土体饱和。因此，认为在单相流分析中不能考虑孔隙空气的运动，水以一定的速度平稳渗透，从解析结果可以看出试样下部比试样中部更早地达到了饱和状态。在两相流分析中，No.1测点和No.2测点的孔隙压力在大致相同时间内由负变正，但由于有微小的误差，No.2测点测得的结果较快。因此，认为土体饱和是从试样的下部向上部进行的。与单相流分析相比，No.1测点和No.2测点的孔隙压力变为正压所需时间相同，其原因是在两相流分析中，试样在饱和层从底部上升前，中间土层的饱和度较高，使其饱和所需的水量变少。

在模型实验中，孔隙水压呈阶梯式上升，并趋于静水压力值。在这2种情况下，孔隙水压继续缓慢增加，直至静水压力值。尽管实验中使用的试样是由不均匀性的土壤颗粒组成的，而从数值模拟角度认为模型是一个均质连续体，但是，在两相流分析中孔隙水压的上升趋势和实验结果相同，只是No.1测点的孔隙水压比No.2测点先开始上升。对于单相流分析的结果，试样下部先达到饱和的时间要早于实验结果和两相流分析结果，无法再现水的渗透被孔隙空气阻挡的情况，并且也很难捕捉非饱和状态下的孔隙压力变化情况。

由以上分析可知，在50 mm/h降雨强度下，降雨强度相对于土壤渗透性较大时，两相流分析可以考虑孔隙空气对水渗透的影响。反之，当降雨强度相对于土壤渗透率较小时，认为孔隙空气的影响较小，可以通过单相流分析来评价渗流行为。由于单相流分析无法捕捉到非饱和状态下的孔隙压力变化情况，而通过两相流分析可以得到非饱和状态下孔隙压力和土中空气的影响，所以两相流更适合评估降雨渗透导致的边坡稳定性。

3 边坡渗流行为的室内斜面模型

3.1 斜面模型实验

斜面模型实验如图4所示，斜面模型的形状和仪器的位置如图5所示。实验材料使用的是从同仁县隆务镇西山Ⅲ号滑坡现场采取的表层粉质黏土。边坡模型制作时试样的初始含水量为10%，干密度为1.68 g/cm3，黏聚力为12.7 kPa，内摩擦角为36.9°。坡面模型分8块压实，为了达到预定的湿密度，每段都分为3层并用木槌进行压实。为了确认坡面和坡体重的渗透行为，在No.1—No.3测点设置了张力计，No.2和No.3测点分别安装了孔隙压力计和土壤水分仪，并从模型顶部调整喷水量模拟降雨。为验证不同降雨强度对雨水的渗透行为的影响，进行了3种不同模拟降雨强度条件下的实验，即在45°倾角下降雨强度分别为20，50和100 mm/h。

图4 斜面模型实验Fig.4 Slope model experiment

3.2 边坡模型实验再现数值模型

笔者根据图5所示的模型进行数值模拟。在相同实验条件下，在同一位置测量孔隙水压。在数值模型中，通过对比斜面模型试验中测得的孔隙水压和吸力值与数值模拟得到的孔隙压力和吸力值，对模型进行了验证。应用两相流分析需要的土壤水分参数b=1/2，c=1/3。在模拟中，应用了表示非饱和时的水分特征曲线的式(8)。P0是土壤的特征毛细管压力，它随土壤类型的不同而变化，a和P0的取值是通过下述实验确定的。式(13)是表现Van Genuchten水分特征曲线的模型[14]。式中的ψ表示基质吸力，Se表示有效饱和度，

(13)

为了确定式(9)中的a和P0的取值，在与3.1固结条件相同情形下，将风化泥岩压实到体积为1 000 cm3的模具中，安装张力计和土壤湿度计，从模具上方供水，测量土壤含水量和基质吸力，直至内部饱和。结果如图6所示，4组的结果几乎相同，以图6中Ec5_1&吸力1的值为代表，与SWRC Fit中的土壤水分特征曲线的非线性回归程序拟合，得出如图7所示的表现水分特征曲线模型。从图7得出式(13)中的参数α=0.125，n=1.49。反转横轴和纵轴，将Van Genuchten模型式(8)转化为如图8所示的曲线，得出a=0.33，P0=8.0 kPa。

图6 体积含水率和基质吸力关系Fig.6 Relationship between volumetric water content and matrix suction

图7 式(13)水分特征曲线Fig.7 Water characteristic curve of Eq.(13)

图8 饱和度和基质吸力关系Fig.8 Relationship between saturation and matrix suction

3.3 实验结果和数值模拟结果

图9表示了实验和数值模拟过程中孔隙水压的变化情况。实验开始时孔隙水压大小不同的原因是由于模型制作时压实过程产生的振动，导致模型中的孔隙水向下移动。通过调整模型中的初始饱和度，使分析实验在降雨开始后立即获得与实验中相同的孔隙水压值。在降雨开始后，因为No.2测点比No.1测点更靠近坡面，No.2测点的孔隙水压立即开始增加，随后No.1测点的孔隙水压才有上升趋势。当浸润面到达计测点时，孔隙水压迅速增加，在饱和度维持一段时间后，模型下方上升的水面到另一测点，孔隙水压都是先增加一次，然后在短时间内又开始增加，呈现阶梯变化。数值分析表明，孔隙水压并不像实验中那样呈阶梯式增加，而是连续增加。这可能是因为和第2节的一次元渗透流模型实验相同，无法充分再现现阶段土壤孔隙中的水渗透行为。但当浸润面已经达到计测点时，孔隙水压突然上升的行为与实验中的是一样的。其次，孔隙水压开始增加后，即浸润面到达计测点的顺序和时间与实验结果差别也不大。此外，降雨强度为100 mm/h的孔隙水压结果如图9 c)所示，与其他2种情况不同，计测点No.2测点的孔隙水压也同样比No.3测点的孔隙水压增加得快。从这些结果看，虽然无法准确再现实际边坡渗透的复杂行为，但总体上可以再现渗流行为。

图9 孔隙水压变化Fig.9 Changes of pore pressures

4 基于暴雨引起的滑坡数值模拟

4.1 边坡崩塌现场的描述

笔者通过气液两相流分析对2012年青海省黄南州同仁县隆务镇西山发生的滑坡进行了模拟。滑坡全貌如图10所示。

图10 黄南州同仁县隆务镇西山滑坡全貌Fig.10 General view of Xishan Landslide in Longwu Town, Tongren County, Huangnan prefecture

以西山Ⅲ号滑坡为研究对象，滑坡岩性以强风化泥岩堆积的粉质黏土为主，土质不均匀，泥岩风化较严重，局部含少量砾石。强风化泥岩节理裂隙发育，风化程度较高，且形成的粉质黏土不但自身强度低，而且对含水量变化极为敏感，遇水后易软化，抗剪强度大大降低，易形成软弱滑动面(带)，为易滑地层。滑坡体平均坡度约35°，由浅层滑坡形成，滑体轴向长为280～310 m，横向宽为50～180 m，平均厚度为21.5 m。建议采用抗滑桩和桩间挡土墙的组合形式进行边坡防治。

4.2 基于两相流分析的渗透行为评价

图11表示滑坡的数值模型。由于滑坡是强风化泥岩土层的表层塌陷，因此设定作为基层土的弱风化岩是弹性体，而表层的强风化泥岩土是按照Mohr-Coulomb断裂标准的弹塑性体。在边界条件方面，模型的右侧和底部设计为不排水，而左侧和坡面设计为允许排水。底面为垂直和水平方向约束，两侧为水平方向约束。弱风化岩土层和强风化岩土层的初始饱和度分别设定为1.0和0.3，降雨强度设定为100 mm/h。分析开始时，对模型进行自重模拟，将计算后的稳定状态定义为初始应力状态。

图11 滑坡再现数值模型Fig.11 Numerical model of landslide reproduction

模拟结果显示：降雨开始2 h后，坡面附近饱和度增加，可以确认到坡脚孔隙水压开始转为正值。降雨开始3 h后，边坡上部孔隙水压转为正值，强风化泥岩土层孔隙气压增大。降雨开始4 h后，上部的孔隙压力继续增大，孔隙水压增加的位置孔隙气压也在增加。降雨开始6 h后，孔隙水压力的增加领域的面积逐渐扩大，在剪切应变增大的地方，逐渐失稳，上部坡面浅部最大剪应变增大，可以确定为表层滑坡。结果显示，塌方发生在坡面，该处塌方最大深度约3 m，现场滑坡发生的位置与数值模拟确定的划坡位置高度一致。

数值模拟结果表明，滑坡形态和滑坡位置与实际滑坡高度相似。边坡滑动的机理是由于雨水的渗透，孔隙水压力和孔隙压力增大，而有效应力减小。所以表层滑坡不一定是由地下水位上升导致滑坡面有效应力下降引起的，在非饱和条件下也可能发生剪切应变，即在降雨量达到足以抬高地下水之前就有发生滑坡的可能性。