膨胀土边坡非饱和渗流及渐进性破坏耦合分析

2020-07-21陈亮胜韦秉旭张寒冰

水文地质工程地质 2020年4期

陈亮胜，韦秉旭，廖欢，张寒冰

(1.长沙理工大学交通运输工程学院，湖南长沙 410114；2.湖南科技大学土木工程学院，湖南湘潭 411201)

膨胀土边坡在降雨入渗下易产生膨胀变形和强度衰减[1-2]，与普通黏土受降雨影响的渗流特性及稳定性分析的规律不同，导致膨胀土边坡失稳的现象屡见不鲜，且边坡失稳常常具有浅层性和多重滑动性的特点，对边坡防护及公路运营的威胁不言而喻[3-4]。因此，研究膨胀土经降雨入渗条件下产生的膨胀理论体系，分析降雨诱发膨胀土边坡渐进性破坏的过程，对膨胀土边坡防护及治理能起到至关重要的作用。

目前许多学者采用原位监测和模型试验的手段，来模拟膨胀土边坡的非饱和渗流和渐进性破坏过程。如詹良通等[5]建立人工降雨边坡模型，监测了含水率、孔隙水压力、土体变形和应力状态随降雨入渗的变化规律，发现了吸力消散和膨胀变形致使膨胀土的抗剪强度逐渐软化；梁树等[6]结合原位试验和模型试验，研究了不同降雨条件下膨胀土基坑边坡的入渗规律，并阐述了降雨入渗深度的特征。但原位观测和模型试验具有耗时周期长、受外界条件影响大、步骤繁琐等特点，且在反映不同降雨时刻的边坡破坏特征表现时，难以详细观测边坡内部滑动面渐进性破坏的过程。采用数值模拟法能清晰地反映不同降雨时刻对边坡浅层破坏和渐进性破坏的动态变化过程，以往的研究常采用考虑裂隙和降低强度的数值方法[7-8]，通过非耦合数值方法分析边坡内部非饱和渗流演变规律，模拟分析边坡稳定性，但模拟结果不能合理反映流-固耦合作用下随含水率动态变化的膨胀变形和强度软化。基于温固耦合或二次开发的多场耦合数值分析方法克服了上述缺点，逐渐运用于降雨条件下膨胀土浅层失稳和渐进性破坏的研究。王凯等[9]根据膨胀岩膨胀变形理论和增量理论，通过自编程序引入膨胀力的弹塑性模型，发现膨胀力引起摩尔应力圆偏移，巷道两侧和顶部出现剪切塑性区。李朝阳等[10]通过室内试验考虑变形参数、强度参数和膨胀参数随含水量变化之间的关系，提出了基于Mohr-Coulomb准则的膨胀土弹塑性模型，并通过案例验证了模型的可行性。QI等[11]利用弹塑性力学原理解释了膨胀土吸湿膨胀引起的应力状态变化，研究膨胀土边坡破坏时间和破坏深度与不同初始应力、软化时间和边坡坡度的关系。但上述研究在考虑膨胀软化的基础上，没有实现膨胀土边坡渐进性破坏的过程。为此，丁金华等[12]基于湿度应力理论，利用温度系数等效膨胀系数，实现了温度场模拟膨胀变形场的多场耦合数值模拟，并结合物理模型试验验证了多重滑动性和渐进性的破坏过程，但温度系数只能近似反映膨胀土的膨胀特性。张良以等[13]采用C++编制了膨胀土多场耦合的计算程序，以此探讨膨胀模量和应变软化对膨胀土边坡的渐进性破坏的影响，但膨胀模量的物理意义不明确。

为研究降雨入渗条件下膨胀土边坡的多场耦合数值模拟，采用饱和-非饱和渗流理论、膨胀土弹塑性本构关系和应变软化理论，将膨胀变形引入非饱和流-固耦合模型中，利用应变软化模型、FLAC3D二次开发平台和内置FISH语言，实现非饱和渗流-膨胀变形-应变软化多场耦合作用下膨胀土边坡的数值模拟，并分析降雨入渗条件下膨胀土边坡非饱和渗流、位移响应及渐进性破坏的动态变化过程。

1 非饱和膨胀土理论基础

1.1 饱和-非饱和渗流理论

非饱和土水力特性主要为土水特征曲线(SWCC)和渗透函数曲线(HCF)，分别用来表示非饱和土的持水性能和导水能力。VG模型[14]是膨胀土常采用的土水特征曲线拟合公式，因其参数容易获取，且拟合效果较好，广泛运用于非饱和土渗流分析。单相渗流模拟的基质吸力忽略了气压力的影响，为负的孔隙水压力，并将体积含水率转化为饱和度，形成饱和度与孔隙水压力的关系表达式：

(1)

式中：s，sr——饱和度和残余饱和度；

uw——孔隙水压力；

a，n，m——拟合参数，满足0

渗透系数函数采用SWCC曲线间接确定：

(2)

有限差分元FLAC3D计算非饱和土内流体扩散作用时，水分运移过程遵循达西定律：

qi=-Kilk(s)(uw-ρfxjgj)il

(3)

式中：qi——单位流量向量；

Kil——渗透系数张量；

k(s)——相对渗透系数；

ρf——流体密度；

xj——笛卡尔坐标分量；

gj——重力加速度分量。

针对非饱和多孔介质结构，其连续性方程为：

(4)

式中：ζ——孔隙水运移过程中发生的流体体积改变量；

M——比奥模量；

n——孔隙率；

t——时间；

α——比奥系数。

当uw≥0时，此时为饱和多孔介质渗流问题，饱和度设置为1，平衡方程更新为：

(5)

FLAC3D内置渗流边界有固定压力边界、固定流量边界和渗漏边界三种，研究降雨入渗多孔介质单相渗流的坡面入渗与产流边界条件时，主要采用前两种边界。其中，压力边界Γw为：

(6)

流量边界Γv为：

(7)

式中：

已知边界孔隙水的压力和流速。

1.2 膨胀土基本理论

降雨入渗作用下膨胀土产生膨胀变形，缪协兴[15]基于湿度应力场理论，提出了膨胀土的膨胀变形与膨胀系数和含水率的变化量有关，其弹性阶段的自由膨胀应变表达式为：

(8)

式中：β——膨胀系数；

Δw——含水率的变化量；

(9)

将式(9)代入式(4)中，得到膨胀土连续性方程：

(10)

非饱和区土体产生负孔隙水压力，引起有效应力增加，采用 Bishop提出的有效应力理论，并认为Bishop参数等于饱和度[16]：

σij=σ′ij+suw

(11)

式中：σij，σ′ij——总应力和有效应力。

根据有效应力原理，将土体结构的静力平衡方程更新为考虑有效应力的格式：

(12)

式中：ρd——干密度；

γw——水的重度。

同时各向同性弹性体的各个方向相同，其应变分量可表示为：

(13)

式中：E——杨氏模量；

ν——泊松比。

根据式(13)，写出总应力分量表达式：

(14)

1.3 应变软化理论

膨胀土在水分运移过程中，同样会引起材料出现强度衰减或应变软化的特征。为描述在不同应变状态下强度的变化规律，以及再现边坡渐进性破坏的全过程，采用应变软化理论模拟膨胀土的应变软化，其剪切屈服准则和张拉屈服准则为：

(15)

Nφ=(1-sinφ)/(1+sinφ)

式中：c，φ——膨胀土黏聚力和内摩擦角；

σt——抗拉强度。

随膨胀土边坡内部含水量和应力应变关系的变化，当应力达到极限状态时，强度参数出现软化劣化效应，逐渐下降至某一残余强度。对计算单元内部进行数值模拟过程中，设置单元剪切应变增量Δκs和张拉应变增量Δkt，表达式为：

(16)

式中：

最大和最小塑性剪切应变增量；

2 膨胀-渗流-软化多场耦合分析法

为了实现降雨条件下膨胀土边坡多场耦合数值模拟，从非饱和渗流、膨胀变形和应变软化三大数值模拟的技术难题着手处理。

(1)采用FISH语言控制FLAC3D渗流模块(Fluid)的方法，可用于实现非饱和渗流的数值模拟。首先，通过设置流体抗压强度允许负孔隙水压力的存在，利用变量zone.pp生成初始孔压场；其次，采用变量gp.sat和zone.property控制公式(1-3)，用以修正各土层的饱和度和渗透系数；再次，利用intface功能获取坡面渗透系数，开发降雨入渗与产流边界；最后，通过FISHCALL命令实现动态时间步下非饱和降雨入渗的数值模拟。其非饱和渗流具体实现过程见文献[17]，此不赘述。

(2)根据文献[18]提出的测定方法，对南阳膨胀土进行膨胀率试验，得出了膨胀变形εw随初始含水率w0和含水率w的变化关系曲线(图1)，并进行多元非线性回归，膨胀变形经拟合的表达式为：

(17)

图1 膨胀变形随含水率变化Fig.1 Expansion rate varying with water content

基于Visual Studio 2010平台上有限差分FLAC3D本构模型的二次开发接口，采用C++编程语言，在FLAC3D内置模型代码基础上，根据增量理论修正头文件(.h)和源文件(.cpp)，并合成动态数据库文件(.dll)，在本构模型中实现膨胀变形的施加。新定义的本构模型(Plugin)既能体现膨胀土因应力状态而改变的应力应变，也能体现因含水率变化而改变的膨胀变形。

(3)采用FISH语言更新有效应力和非饱和抗剪强度参数，调用修正后的应变软化模型来模拟土体应变软化的过程。当膨胀土边坡内部土体进入塑性屈服状态时，峰值强度τp逐渐转化为残余强度τr，抗剪强度参数c和φ随塑性应变的增加而非线性减小：

(18)

式中：kc，kφ——黏聚力和内摩擦角随塑性应变变化的软化系数，其取值如表1所示；

cp，φp——峰值黏聚力和内摩擦角。

表1 强度参数软化系数取值

以上为三种技术手段的实现思路，其多场耦合数值模拟的计算流程如图2所示。

图2 计算流程Fig.2 Numerical simulation calculation process

3 数值模型及计算方案

3.1 数值模型

河南境内南阳—邓州高速公路某黏土边坡高度为10 m，坡比为1∶1.5，经现场地质勘探发现，该边坡地下水充沛，其出露部分为Q2-3洪积型黏土，由泥灰岩和少量变质岩系风化，再经流水搬运作用洪积而成，是典型的南阳膨胀土结构层。该边坡常年暴露于复杂的气候条件下，边坡深度 3 m范围内分布着不均匀开裂的膨胀土，但裂隙发育的规模及程度较低。

要模拟膨胀土边坡的实际地质情况比较复杂，因此，在保证数值分析可靠性的基础上，根据膨胀土边坡地质情况进行了适当的假定：

(1)根据裂隙随深度自上而下的变化情况，将土层划分为上层土、中层土和下层土三种类型；

(2)视各土层的物理参数为常数，并取各层代表性土样作为试验对象；

(3)各土层为等效渗透系数的连续介质模型，且层间接触条件为完全连续接触；

(4)边坡初始湿度恒定，初始状态下不存在膨胀变形。

图3 数值模型Fig.3 Numerical model

结合实际工程及假定条件，建立如图3所示膨胀土边坡数值模型，上层土、中层土的厚度分别为2，1 m，初始地下水位水平，位于模型底部以上2 m，网格划分为边长为0.25 m的六面体。经土工试验获取各土层的物理参数，其参数取值如表2所示，通过室内压力板试验和双环试验测得各层土的土水特征曲线和饱和渗透系数，采用VG模型进行拟合如图4所示。

表2 各土层物理参数

图4 膨胀土水力特性Fig.4 Hydraulic characteristics of the expansive soil

3.2 边界条件

模型底部为水平及竖向约束的速度边界，模型四周为水平约束的速度边界，并设置模型底部和四周为不透水边界。值得注意的是，在降雨过程中，当降雨强度大于土体的入渗能力时，坡面入渗量由入渗能力决定，此时边界条件取零压力边界；当降雨强度小于土体饱和渗透系数时，坡面入渗量由降雨强度决定，此时边界条件取流量边界。在降雨结束后，一种情况是受坡面出渗的影响，引起坡内水分向外溢出，坡面孔隙水压力逐渐转变为正值，此时边界条件设置为零压力边界；另一种情况是受重力的作用，水分继续向坡内运移，坡面孔隙水压力逐渐转变为负值，此时边界条件设置为零流量边界。

3.3 初始孔压场及降雨方案

在长期自然环境下，坡面不同位置的干湿程度接近，表面基质吸力基本一致。根据坡内地下水位位置和坡面表面基质吸力的大小，将坡面孔隙水压力设置为-180 kPa，地下水位设置为0 kPa，水位至坡面的孔隙水压力非线性分布，生成初始孔压场(图5)。

图5 初始孔隙水压力Fig.5 Initial pore water pressure

为考虑最不利降雨情况对膨胀土边坡非饱和渗流特性的影响，依据河南省南阳地区气象观测极值确定降雨工况，选取降雨强度为1.65×10-6m/s，并拟定降雨历时3 d，停雨历时6 d。

4 数值结果及分析

4.1 非饱和渗流分析

暂态饱和区指在降雨入渗过程中，坡体饱和度随降雨历时的变化发生改变，逐渐在坡内形成饱和度大于0.9的区域[19]。暂态饱和区厚度和面积的动态变化体现着土质边坡的非饱和渗流时空分布规律。为此，首先对膨胀-渗流-应力多场耦合作用下暂态饱和区的时空变化规律展开分析。

图6所示为不同降雨时刻饱和度的演化规律云图，降雨1 d时，坡面出现暂态饱和区，随降雨持续到第3天，饱和度为0.75以下的非饱和区已完全消散，暂态饱和区随着时间的推移呈悬挂型向坡体内部扩展，暂态饱和区的范围与降雨入渗深度呈现明显的正相关关系。这是由于降雨强度大于土体入渗能力，湿润锋呈饱和状态向下扩展。降雨结束后，一部分水分由于重力作用湿润锋继续向下扩展，坡面的水分向坡体内部渗透，补充着湿润锋下缘附近的水源；另一部分水分由于坡面出渗向坡外溢出，但水分向坡外溢出的速度落后于向坡内渗透的速度，且下层土的入渗能力限制了湿润锋向下扩展的速度，此时坡面基质吸力的消散速度低于悬挂型暂态饱和区的增长速度，导致停雨0～6 d内的暂态饱和区大于降雨期间暂态饱和区的范围，暂态饱和区呈显著的时间滞后现象。但湿润锋被均化，停雨后暂态饱和区内的整体饱和度有所缩减。

图6 饱和度变化规律Fig.6 Evolution of saturation

膨胀土非饱和渗流的过程伴随着膨胀变形的产生，为体现膨胀土的膨胀特性，分别模拟不考虑膨胀变形和考虑膨胀变形的渗流耦合分析，将前者视为普通黏土，后者视为膨胀土。图7为截面x=10 m和x=20 m 处基质吸力随高程的变化规律。从图7中可看出：降雨3 d时，湿润锋呈饱和状态向坡内扩展，膨胀土的湿润锋深度大于普通黏土；停雨6 d时，湿润锋被均化，受湿润锋影响的区域扩大，在考虑膨胀变形的基础上，基质吸力变化幅度的范围缩小。造成上述现象的原因是：降雨入渗期间，膨胀土吸水产生膨胀变形，使土体受压而变得密实，阻碍了湿润锋向坡内扩展，对膨胀土的降雨入渗过程起抑制作用，造成湿润锋的扩展深度减小。当降雨结束后，湿润锋上缘基质吸力有所恢复，湿润锋下缘继续扩展，暂态饱和区具有时间滞后性，受湿润锋影响任意位置的基质吸力均大于初始状态的基质吸力，膨胀变形对水分运移的过程继续发挥着抑制作用，湿润锋上缘被均化的速度和下缘扩展的速度受膨胀变形的调节。

图7 基质吸力随高程的变化Fig.7 Matrix suction changes with the slope depth

对比图7(a)和(b)中坡脚附近和坡顶附近基质吸力变化可以发现，降雨3 d时，坡面吸力由降雨前的180 kPa，逐渐消散至0，且膨胀土和普通黏土坡脚和坡顶处基质吸力的变化规律相同。停雨6 d时，由于膨胀土的膨胀变形限制了水分运移的速度，膨胀土基质吸力在坡顶和坡脚处基质吸力变化不明显，而普通黏土在坡顶附近的基质吸力明显高于坡脚附近。此外，坡脚附近湿润锋深度约3.5 m，坡顶附近湿润锋深度约3.25 m，坡脚附近湿润锋的影响程度高于坡顶附近。这是由于水分难以透过渗透率低的下层土，汇集于中下土层交界处，加上水分沿坡面方向由坡顶向坡脚运移，造成坡脚处基质吸力的消散速度高于坡顶处。

4.2 位移响应及渐进性破坏分析

非饱和渗流具有显著的时空分布不均匀性，引起的膨胀变形和应变软化对膨胀土边坡具有不同的位移响应程度。为探究膨胀-渗流-应力多场耦合作用下膨胀土边坡的位移响应特征，综合计算了4种位移监测方案：(1)考虑膨胀变形和应变软化；(2)考虑膨胀变形、不考虑应变软化；(3)考虑应变软化、不考虑膨胀变形；(4)不考虑膨胀变形和应变软化。在x=10 m和x=20 m的坡顶位置分别设置2个A、B监测位移的特征点，图8为特征点A和B位移的时程变化。

图8 不同位置位移时程变化Fig.8 Displacement response of different slope locations

针对膨胀土边坡位移响应的影响，从图8中不难发现，4种方案敏感性分析的先后顺序为：考虑膨胀和软化、考虑膨胀、考虑软化，以及不考虑膨胀和软化。模拟结果显示了膨胀变形和应变软化对膨胀土位移响应的效果，体现了多场耦合作用下膨胀土边坡的非饱和渗流、膨胀变形和应变软化之间的相互联系。另外，特征点A各方案最大位移分别为75.6，60.4，41.1和24.3 mm，特征点B各方案最大位移分别为69.8，55.9，37.3和21.5 mm，各方案最大位移发生的时间分别为9，8，6，5 d附近，均处于停雨阶段，且坡脚附近整体位移响应程度大于坡顶附近。这是由于湿润锋在降雨结束后会继续向边坡内部和坡脚扩展，造成坡脚剪应力集中和基质吸力消散速度变快。

膨胀变形和应变软化对边坡位移响应产生一定影响，但影响效果不同：

(1)降雨初期坡体受雨水入渗引起应变软化效应不明显，应变软化对边坡变形无显著影响，此时方案1和方案2的位移响应相近。但随着降雨持续和湿润锋进一步向下扩展，受应变软化影响的区域增加，应变软化对位移响应程度逐渐加深。当不考虑应变软化时，难以如实反映降雨入渗条件下边坡变形和强度衰减的变化。

(2)边坡土体受降雨入渗影响而湿化，产生的膨胀变形对边坡变形的影响贯彻水分运移的整个阶段，膨胀变形对位移响应的影响显著大于应变软化。当不考虑膨胀变形时，膨胀土将退化为普通黏土，边坡位移响应程度急剧减小，导致边坡稳定性被高估。

根据不同时刻塑性应变的变化规律(图9)，可以看出：第一次滑动发生在降雨期间，塑性破坏区垂直于坡面向坡内扩展，产生局部塑性破坏；第二次滑动发生在停雨阶段，塑性破坏区由坡脚向坡顶发展，逐渐形成贯通的滑动面，最终导致整体性失稳。其渐进性破坏过程呈现出多重滑动性和后退牵引性的破坏特征，且整体性破坏发生在雨后，具有显著的时间滞后性。

图9 塑性应变演化规律Fig.9 Evolution of maximum shear strain

产生上述现象的原因是：在降雨过程中，坡面形成悬挂型暂态饱和区，暂态饱和区的形成引起饱和区内基质吸力消散、抗剪强度降低、产生膨胀变形和土体重度增加，导致饱和区内竖向应力和水平应力分布状态增大，大主应力的方向发生偏转，逐渐平行于坡面方向。此外，浅层土受上覆荷载作用小，沿坡面法向方向的约束能力差，在变形协调作用影响下，会产生较大的剪应力，形成第一条平行于坡面向坡内扩展的剪切滑动带。当降雨结束后，随着湿润锋被均化，湿润锋上缘含水率下降，下缘抵达中下土层分界面上，逐渐形成饱和-非饱和的交界带。饱和带停滞时间长，饱和带附近的膨胀变形迅速增加，且土层上下表现出显著的强度差异，下层土的抗剪强度远高于中层土，产生大主应力集中的现象。在第一次滑动基础上，坡面土体已遭受塑性破坏，不能有效支撑边坡内部土体。此外，坡脚附近的含水率高于坡顶附近，坡脚承受更大的膨胀变形和应变软化双重作用，加之边坡的形态结构，率先在坡脚附近形成塑性破坏区，并逐渐沿中下土层分界面由坡脚向坡顶扩展，形成第二条向上贯通的剪切滑动带，最终导致边坡产生后退牵引式的整体性破坏。