基于三场耦合的膨胀土边坡浅层坍滑数值分析*

2022-10-27谭越峰杨和平邹维列陈冠一

工业建筑 2022年7期

肖杰谭越峰童超杨和平常锦,3 邹维列陈冠一, 5

(1.长沙理工大学交通运输工程学院，长沙 410004； 2.中铁科学研究院有限公司，成都 611731；3.长沙学院土木工程学院，长沙 410022； 4.武汉大学土木建筑工程学院，武汉 430072；5.广东鸿高建设集团有限公司，广东东莞 523507)

膨胀土是道路工程、岩土工程和边坡工程常遇的一种区域性特殊土，主要包含蒙脱石和伊利石类强亲水性矿物，具有快速崩解性和胀缩特性[1]。旱季膨胀土坡面失水收缩开裂且裂隙发育，雨季的降水入渗，在边坡体内形成非饱和-饱和渗流以及浅表层的吸湿膨胀，会大大降低其稳定性而造成浅层坍滑，给膨胀土地区工程建设带来极大困扰并产生严重的工程和经济损失[2]。平缓且低矮膨胀土边坡发生的浅层坍滑，按传统的强度理论无法对其破坏原因或机理做出合理解释。为此，大量学者指出在进行膨胀土边坡稳定分析时，除需考虑遇水后土体的强度降低外，亦须考虑膨胀变形的影响[3]。肖杰等探究了干湿循环及低应力条件下膨胀土抗剪强度的衰减及其对边坡浅层坍滑的影响，但未考虑边坡浅、表层土的吸湿膨胀作用对这种衰减、破坏的影响[4]。文献[5-6]介绍了对干湿循环下膨胀土微观结构、体变及力学性能变化规律的研究，指出膨胀土在干湿循环下极易产生裂缝，从而使原土体结构破坏，渗透性显著提高，结构整体性减弱，强度、刚度均发生较大衰减。方瑾瑾等通过真三轴试验得出了相似结论，解释了强度、刚度发生较大衰减原因是干湿循环降低了膨胀土黏聚力及内摩擦角，且发现膨胀土基质吸力与含水率变化路径有关[7]。文献[8-9]指出干湿循环对膨胀土土体持水特性和收缩特性具有显著影响。学者们还基于湿度应力场理论，通过FLAC2D的内置热模块，将数值反演获得的膨胀力类同于重力作用进行计算，一定程度上反映了湿度和变形耦合效应，其实质是采用线性插值将含水率变化等效成温度的均匀变化[10-12]。丁金华等也基于湿度应力场理论，由数值反演有或无荷载作用的室内试验结果，得到饱和度与温度及热膨胀系数间转换关系，能较真实地模拟膨胀土的实时吸湿膨胀效应，并采用多场耦合的方法，较深入地探究了膨胀土吸湿膨胀作用对边坡稳定性的影响，但未考虑干湿循环效应对边坡膨胀土的抗剪强度、渗透系数的影响[13]。张良以等基于膨胀应变与基质吸力增量线性关系，建立了渗流场、应力场、应变场多场耦合数值计算方法，并结合应力软化模型，探究了单次降雨过程膨胀土边坡渐进破坏过程，对比分析了不同膨胀性、不同材料强度参数对边坡破坏形式的影响，但却未考虑多次干湿循环对膨胀土材料特性的影响[14]。降雨入渗显著影响膨胀土边坡的稳定性早已是国内外学者的共识[3,15-16]，但迄今未见同时考虑土体强度衰减和吸湿膨胀变形共同作用对不同降雨强度及历时条件边坡浅层坍滑影响的专门研究。

为此，借助FLAC2D有限差分软件，借鉴多场耦合数值计算思路及方法，综合考虑土体抗剪强度衰减及膨胀土吸湿膨胀两者共同作用，深入研究复杂环境下膨胀土边坡工程性状，进一步分析不同降雨强度和历时条件对边坡稳定性影响，以探究膨胀土边坡浅层坍滑的原因。

1 膨胀土边坡三场耦合实现流程

1.1 热-固耦合分析模型

各向同性热传导边坡模型如图1所示，左、右两侧边界约束水平向位移，底边界约束水平和竖向位移，边坡模型底边长为21 m，坡底高为6 m，长为6 m，坡顶高为12 m，坡面高度为6 m，坡比1.0∶1.5。因软件无法事先确定边坡的滑裂面，故得不到滑面的法向应力，也就无法直接采用非线性强度拟合函数展开计算。文献[17]指出：不同干湿循环次数及低围压条件下，双直线与广义幂函数法拟合的膨胀土抗剪强度值在高、低围压段之差均不大，可用双直线拟合值替代广义幂函数拟合值做数值计算。肖杰等开展了低应力条件不同干湿循环次数的百色膨胀土饱和三轴固结排水试验，发现经干湿循环作用达到饱和后的土样可压缩性大幅增加，试验中橡皮膜对施低围压土样的环箍效应显著，为使测试值准确可信，重点分析了橡皮膜的约束对试验围压的影响，专门做出相应的强度参数校正，获得强度参数随干湿循环次数的衰减规律[17]。由此，下述数值模拟也将边坡由表及里划分为三层，各层强度参数取文献[17]中双直线法经0次干湿循环确定的低、高围压段拟合值。根据现场调查分析，百色土的一、二层定为强、弱风化层，厚度均取1.0 m，其强度参数都从拟合强度线低应力段取值，余下坡体土均划为第三层，即未风化层，强度参数取高应力段的拟合值。其剪胀角取20°，表1给出各力学参数取值。

表1 百色膨胀土的力学参数取值Table 1 The mechanical parameters of Baise expansive soil

分析中，热传导系数取0.49 W/(m·K)，比热容取2 100 J/(kg ·K)。1)求解初始地应力平衡时，初始状态的温度场取0 ℃，热膨胀系数取零；2)模拟降雨条件下边坡土体的吸湿膨胀时，因边坡的坍滑在浅表层发生，故取地应力σ=6.25 kPa对应的饱和度与热膨胀系数及温度的转换关系[13]：

T=416.7Sr-316.7

(1a)

λ=(20.83Sr-15.83)×10-4

(1b)

式中：Sr为饱和度；T为温度；λ为热膨胀系数。

1.2 非饱和二相流分析模型

二相流模块假定孔隙完全由两种流体填满，且两者间不发生质量转移，故可用于研究非饱和膨胀土的渗流问题。其传导定律可描述为：

(2)

式中：上角w、g为液、气相。

用达西定律可表示为：

(3a)

(3b)

相对渗透系数kr在Van-Genuchten水土特征曲线算式中与有效饱和度Sre的关系为：

(4a)

(4b)

毛细压力pc用Van-Genuchten模型可表示为：

(5a)

(5b)

式中：η为经验系数。

二相流数值计算模型的尺寸大小、网格剖分与热-固耦合模型相同，设地下水位水平，距图形底部为1 m，采用瞬态流降雨入渗法模拟初始状态下的稳定渗流场分布，详见图2。

表2 非饱和-饱和渗流参数取值Table 2 Values for seepage parameters from the unsaturation to saturation

坡面裂隙显著地影响边坡渗流场、应力场及其稳定性，其影响通常由分层等效法考虑：将边坡体由表及里按风化程度划分为三层，各土层饱和渗透系数分别设为2.3×10-6,2.3×10-7,2.3×10-8m/s，为方便分析计算将其换算成FLAC设定形式，即分别为2.346×10-10,2.346×10-11,2.346×10-12m2/(Pa·s)。

位移边界的设置同热力学模型。设底面为不透水边界，水位线以下为压力边界，且饱和度取100%。按国家气象部门规定的降雨量标准，降雨强度分小雨(10 mm/d)、中雨(20 mm/d)、大雨(50 mm/d)及暴雨(100 mm/d)四级，为方便计算分别将它换算为1.157×10-7，2.315×10-7，5.787×10-7，1.157×10-6m/s。因四级降雨强度均小于强风化层土体的饱和渗透系数，故设边坡表面为流量边界。模型底面及左、右两侧均设为不透气边界，表面气压值设为零。

1.3 渗流-位移-应力三场耦合计算流程

现有FLAC软件各版本均未提供热力学与二相流模块的实时耦合计算功能，为传递两模块间信息，自行编制了相应FISH程序，最终实现考虑地下水位、渗透系数变化及强度衰减三因素下的膨胀土边坡三场耦合。具体流程如下：

1)分别赋予其相应的参数、边界及初始条件，以实现两模块间信息相互调用。

2)采用热力学模块求边坡初始地应力平衡并输出其应力场、位移场。基于瞬态流降雨入渗法，编FISH程序，将上述应力场、位移场作为初始条件，由二相流模块模拟天然状态下的稳定渗流场，用INITIAL命令清零位移场，输出其应力及饱和度场。据此开展初始状态下非饱和渗流计算，其信息的输出用INFO命令。

3)由式(1)将二相流模拟所得饱和度场换算成温度场及热膨胀系数，输出其应力场，编FISH程序将三者赋予热力学模块，以实时模拟各结点吸湿膨胀效应并求解热、力学平衡。

4)编FISH程序将热力学模块所得应力场赋予非饱和二相流模型作下一降雨时段渗流分析初始条件，重复步骤2)、3)至分析结束，最终实现三因素下的膨胀土坡的三场耦合。

2 膨胀土吸湿膨胀对边坡稳定性的影响

以下分析除直接读取数值原图外，还提取三处控制断面处的孔隙水压、安全系数及潜在滑面深度、坡体变形绘制的相应图表，以方便直观比对吸湿膨胀效应对边坡稳定性影响。

2.1 边坡渗流场影响分析

图4～6分别为是否考虑膨胀土吸湿膨胀，降中雨时边坡的孔隙水压云，以及边坡三断面处孔隙水压随降雨历时的变化状况。结合两图分析：1)考虑与未考虑吸湿膨胀作用两坡体内的孔隙水压增长趋势大致相同，均随降雨历时增长而增大，但考虑吸湿膨胀作用的增幅及影响范围均远大于未考虑的同一边坡。2)考虑吸湿膨胀的边坡，降中雨5 d后，坡脚(x=5 m)处强、弱两土层基本呈饱和状态，饱和区范围达2.2 m，而坡中(x=10.5 m)、坡顶(x=16 m)两断面处分别约为2.0，0.8 m。3)未考虑吸湿膨胀的边坡，除地下水引起的饱和区外，三控制断面处一直都处于非饱和状态，孔隙水压的变化幅度不大。

2.2 边坡位移矢量分析

图7、8为是否考虑膨胀土吸湿膨胀，降中雨时边坡的水平、竖向位移云。图9～12为考虑膨胀土吸湿膨胀与否，随降雨历时变化的边坡三断面处水平及竖向位移的对比。分析图7～12可得：1)考虑吸湿膨胀作用下，边坡三控制断面处水平、竖向位移均比不考虑的同一处大两个数量级。显然，仅考虑雨水入渗后土体强度衰减而忽略其吸湿膨胀效应，分析所得的坡体位移与实际浅层存在的吸湿膨胀后的相差甚远。2)考虑膨胀土吸湿膨胀，坡顶处最大竖向位移均比坡脚(x=5 m)及坡体下部1/3处大，约52.18 mm，但其影响深度却小于后者。原因是坡顶处的水平和竖向均临空，故其竖向位移在一定深度范围发展较快。未考虑吸湿膨胀作用时，坡体下部1/3处及坡脚(x=5 m)处的变化幅度基本一致，且均比坡顶处大。3)无论考虑吸湿膨胀与否，边坡的水平位移变化趋势基本一致，变幅最大均发生在坡体下部约1/3处，坡顶次之，坡脚处最小，究其原因是该处虽水平向临空，但因受到侧向土体的挤压，其水平位移变化幅度仍较坡顶处大。4)有或无吸湿膨胀条件的水平、竖向位移均只能影响边坡的浅表层，且由表及里该影响渐减，这与膨胀土边坡坍滑的浅层性吻合。

2.3 边坡剪应力分析

图13为考虑膨胀土吸湿膨胀与否，降中雨时边坡剪应力云。不难发现：1)两者的剪应力集中区均先在边坡中、下部出现，考虑吸湿膨胀的边坡其集中区随非饱和区-饱和区的扩展向上、深部发展，而不考虑的边坡其剪应力区几乎不变化。究其原因是计及吸湿土膨胀时，坡体内的含水率变化使其浅表层土体产生非均匀变形，从而改变了土中的应力状态。

2.4 边坡安全系数及潜在滑裂面深度分析

受热-固耦合模块限制，不能考虑降雨的影响。为使求解准确，在FISH程序中将热力学模块所得应力场赋予二相流模块后再进行相关计算，求得边坡安全系数。

图14、15分别为考虑膨胀土是否吸湿膨胀，降中雨时边坡的安全系数及潜在滑面深度云。图16为不同降中雨历时边坡的安全系数大小及潜在滑面深度。可知：1)考虑吸湿膨胀的边坡，降中雨4 d后安全系数由1.46降至0.89，减幅39.0%，已发生浅层破坏；而未考虑吸湿膨胀的边坡，降中雨6 d后安全系数仍为1.23，仅减小15.6%，仍处于稳定状态。

2)考虑吸湿膨胀的边坡，降中雨4 d后的潜在滑面深度由1.8 m降至1.36 m，降幅24.4%；而未考虑吸湿膨胀的同一边坡，降中雨6 d后的潜在滑面仍为1.6 m深。因此，膨胀土的吸湿膨胀对降雨入渗边坡的稳定性及潜在滑面深度的影响均不容忽视。

3 降雨条件对边坡三场及浅层坍滑的影响

限于篇幅，仅对比降小、中雨时边坡三控制断面处渗流、位移、应力三场差异，并分析各级降雨强度、历时下边坡的安全性及潜在滑面深度，以获得降雨条件对膨胀土边坡浅层坍滑的影响。

3.1 降雨对边坡渗流场的影响

图17、18为降小、中雨时三断面处孔隙水压的分布。可知：1)降小雨6 d后，边坡坡脚处开始出现正孔隙水压，局部出现0.25 m深度饱和区，与文献[4，19]模拟的相同降雨条件的边坡一直处于非饱和状态的结果不一致，这是因为后两者均未计及吸湿膨胀效应的影响。尽管降雨渗流路径较长，若考虑膨胀土的吸湿膨胀、土变松散，坡体内也可能产生正孔隙水压。2)降雨强度越大，历时越长，孔隙水压的变化幅度及雨水入渗的深度会越大，且坡脚处孔隙水压比坡中及坡顶处的大。

3.2 降雨对边坡位移场的影响

图19～22分别为降小、中雨三断面处经不同降雨历时的水平、竖向位移对比。可知：降雨条件下，水平位移主要集中于边坡下部约1/3处范围。降雨强度越大，历时越长，边坡浅表层土体的水平、竖向位移显著增加，其影响也由表及里逐渐减小，且二者的增长速率在不同降雨条件下也不一致。

3.3 降雨对边坡应力场的影响

图23为降小雨不同历时的膨胀土边坡剪应力云。结合图13、22分析可知：边坡土中剪应力随降雨历时增长显著增大，坡体下部约1/3处出现应力集中，在中下部及坡脚区域先形成集中区，然后逐渐向上、深部发展。且考虑吸湿膨胀与否两条件下的剪应力均随降雨强度增大而增大。

3.4 降雨对边坡安全性及潜在滑面深度的影响

图24为不同降雨条件下边坡的安全系数及潜在滑面深度的对比。不难发现：1)降小雨6 d后，安全系数仍为1.08，降中雨5 d后降至0.89，降大雨30 h后为0.88，降暴雨24 h后仅剩0.66，其降幅分别是39.0%、39.7%、54.8%，边坡均已破坏。2)同降雨历时，降雨强度越大时边坡安全系数及其潜在滑面深越小，安全系数衰减的速率也越快。3)同降雨强度，边坡的安全系数随降雨历时增长而减小。

边坡降雨入渗后，开裂土吸湿、饱和、膨胀，沿临空向非均匀变形并使土中应力重分布，当所受剪力大于土体已衰减强度时局部出现塑性区，恶劣降雨条件下其扩张并贯通，最终导致边坡坍滑。