渗流场与应力场耦合下的边坡稳定性分析

郎 颖

（抚顺市水利勘测设计研究院，辽宁抚顺113008）

渗流场与应力场耦合下的边坡稳定性分析

郎 颖

（抚顺市水利勘测设计研究院，辽宁抚顺113008）

摘要：为了解水库库水位变化对岸坡稳定性的影响情况，建立了有限元计算模型，利用P ha s e 2有限元分析软件，对渗流场和应力场耦合作用下的库岸边坡稳定性模拟计算，其结果为水库安全的运行和管理提供了依据，

关键词：库水位；渗流场；应力场；稳定

D 0I：10.3969/j.i s s n.1672-2469.2016.06.041

水对岩土边坡的影响很大，一般边坡的破坏都受地下水和降雨等水的影响［1]。但是起初人们在对岩土体边坡稳定性进行分析时往往忽略了水的作用，忽略了渗流场对边坡稳定性的影响。直到1970年，随着岩体渗流耦合理论的提出，人们才开始对渗流场和应力场的耦合影响下边坡的稳定性进行研究［2-3]。从开始研究边坡岩土体的应力与其渗透系数、孔隙率、裂隙结构等的关系，到研究围压、孔压与渗透系数的关系，基于以上成果，有学者提出了钉床模型、洞穴模型、洞穴—凸起模型、球形颗粒排列模型和毛细管模型等。这些物理模型可以模拟不同孔隙大小及分布的岩土体中渗流场和应力场的耦合特征。随着计算机技术的发展，有学者建立了数学模型，如等效连续介质模型、网格线素法、裂隙网格模型、电阻网裂隙岩体渗流网格模型、三维裂隙网格模型、双重介质模型以及有限单元法等［4]。经过几十年的发展，渗流场和应力场的理论和研究方法已比较成熟，数值分析方法的作用和影响也越来越大［5]。但也存在着对非连续介质模型计算困难、等效连续介质模型的有限元分析精度和效率不协调、未考虑实际岩体的实际和参数的随机性等问题。本文在对渗流场和应力场耦合的理论分析的基础上，以某水库岸坡为例，研究了库水升降影响下的渗流场和应力场的等效连续介质数值模型，对库岸边坡稳定性进行了数值分析，为水库安全运行和管理提供了依据。

1 渗流场与应力场耦合方法

1.1 洞穴—凸起结合模型

洞穴—凸起结合模型是指用一定形式排列钉状物来表示裂隙面上的凸起，岩体受力时引起的裂隙面的变形情况可以通过钉状物的压缩来表示的物理模型。其裂隙面由两侧凸起和之间的洞穴构成，应力的增加会使裂隙面产生变形，适用于单裂隙岩体的计算。

1.2 球形颗粒排列模型

球形颗粒排列模型是假定岩土体为一个立方体，将均匀分布的球形颗粒看做岩土体中的孔隙，将球形颗粒堆砌在立方体中，球形颗粒的分布表示岩土体中孔隙的分布。模型的渗透系数与孔隙的几何特征符合公式1所示的关系［6]。式中：N为与颗粒几何位置常数；d为与孔隙尺寸有关的参数；μ为粘滞系数，p为围压，g为重力加速度。

将上式中的参数N、d、μ、p、g等采用康采尼理论和铁木辛柯方法进行转换，消除孔隙的几何特性后，得到了介质未固结情况下的球形颗粒排列模型中的渗透系数与应力的关系，如公式2所示［7]。

收稿日期：2015-12-30

作者简介：郎 颖（1975年—），女，工程师。

式中：正负号表示压缩和拉伸加载；v为泊松比；E为固体颗粒有效模量。

1.3 毛细管模型

毛细管模型是用许多毛细管组成不同直径的毛细管束，来代表岩土体中的不同种类的孔隙以及该种类孔隙的渗透特性。模型中的渗透系数与毛细管的多少和直径的大小有关，如公式3所示。式中：Ni为单位面积上第i根毛细管束中毛细管的数目；m为不同直径的毛细管种类数量，φi为第i根毛细管的直径。

根据弹性力学的理论，毛细管的直径会随着岩土体内围压的不同而发生不同的变化，记录这种变化可以得到毛细管模型中岩土体渗透系数与其围压之间的关系，如公式4所示。式中：p为岩土体所受的围压；v为泊松比；E=f （p）为岩体弹性模量。

2 岩体稳定性分析

2.1 工程概况

选取某水库右坝肩一处山体单薄、地下水埋藏深、库水位升降会对边坡稳定造成影响的一段边坡。该边坡长约460m，顶高程约1470m，迎水坡在左侧，底高程约1340m，背水坡位于右侧，底高程约1240m。背水坡较陡，迎水坡较缓。如图1所示。

图1 稳定计算边坡剖面图

边坡岩性自上而下依次为灰岩、大理岩、花岗闪长岩，在花岗闪长岩与大理岩中间夹有一层蚀变形成的石英片岩。岩体的风化基本与地面线平行，从外由内依次为强风化带、弱风化带和微风化带，其中强风化带厚度约有11m，弱风化带厚度约有33m。从左侧迎水坡到右侧背水坡中的初始水位线是先稍微升高，后开始降低，有明显的分水岭，最低水位高程约为1292m。水库蓄水时的正常蓄水位为1392m，每天水位变幅在30m左右。本次主要针对水位自初始水位上升到正常蓄水位后再降低30m，这一过程进行边坡稳定性分析。

2.2 模型建立及边界条件

利用P ha s e 2有限元软件，嵌入渗流场与应力场耦合方程组，对上面的边坡进行二维数值模拟稳定性分析。分析过程中假设各岩层均为各项同性，采用总应力法来分析平面应变问题，岩块和结构面符合摩尔-库伦准则，岩层的抗剪强度参数采用强度折减法进行折减，最终建立等效连续介质模型。如图2所示。模型宽460m，高230m，采用三角形网格进行单元划分，共划分计算单元约10000个。

图2 边坡的有限元计算模型

模拟三种工况下边坡的稳定性：工况1初始水位1293m；工况2正常蓄水位1392m；工况3放水后水位1362m。岩体的抗剪参数为试验室提供和参考《工程地质手册》取得，见表1。层状结晶灰岩层面抗剪强度为c=0.08MPa，φ=31°；层状大理岩层面抗剪强度为c=0.02MPa，φ=31°。结晶灰岩层面倾角为24°，模拟时层厚设定为2.0m；大理岩层面倾角为21°，模拟时层厚设定为2.5m。根据岩体的风化程度将其透水性划分为强透水性、弱透水性和微透水性，其渗透系数见表2。

模型的边界条件有位移边界条件和水力边界条件。各工况的位移边界条件均为边坡两侧边界受x方向约束，底部边界受x、y方向约束。水力边界条件受水位的影响有所区别：工况1左侧水头为1293m，工况2左侧水头为1392m，工况3左侧水头为1362m，各工况的其余部分均为未知地下水边界条件。

表1 岩体抗剪强度参数

表2 岩体渗透系数

3 计算结果分析

3.1 应力场分析

对上面的有限元模型进行上述三种工况下应力场变化模拟，结果如图3所示。

图3 各工况主应力分布云图

由图3中a、b、c图可知：整体上看，最大主应力均与库水位的变化一致，库水位高时最大主应力大，库水位降低时最大主应力减小，各工况的最大主应力分布情况基本一致。从工况1到工况2，随着库水位的升高，岩体内的孔隙水压力增大，最大主应力随之增大；从工况2到工况3，随着库水位的降低，孔隙水压力也减小，最大主应力也随之降低。从工况1到工况3的过程中，最大主应力呈现先升高后降低的趋势。由图3中d、e、f图可知：整体上看，最小主应力也与库水位的变化一致，库水位高时最小主应力大，库水位降低时最小主应力减小，迎水坡的最小主应力大于背水坡。整个边坡在蚀变带内的最小主应力均较大，迎水坡一侧的大于背水坡侧，在背水坡的中后部最小主应力与上下岩体中的区别更明显，明显比上下岩层中的要大，蚀变带上部岩层中的最小主应力在靠近坡体表面时较小，蚀变带处背水坡坡脚处的最小主应力在蚀变带内最小。

3.2 破坏位置分析

各工况下库岸边坡总位移矢量如图4所示。

图4 各工况总位移分布矢量图

由图4可知：各工况下，库岸边坡的破坏位置均在背水坡，背水坡坡脚蚀变带处位移量一直较大，随库水的升高位移值不断增大，但增幅较小，说明库水位的升降对总位移的影响较小，渗流作用及其变化不会引起边坡破坏位置的改变，但是会导致坡体变形增大，加剧坡体的破坏，对边坡的稳定不利。

4 结论

通过建立有限元计算模型，利用P ha s e 2有限元分析软件模拟计算了渗流场和应力场耦合作用下的库岸边坡稳定性，结果表明：各工况的最大主应力分布情况基本一致，最大主应力均与库水位的变化一致，最小主应力也与库水位的变化一致，迎水坡的最小主应力大于背水坡。各工况下，库岸边坡的破坏位置均在背水坡，库水位的升降对总位移的影响较小，渗流作用及其变化不会引起边坡破坏位置的改变，但是会导致坡体变形增大，对边坡的稳定不利。

参考文献

［1]何玉红.降雨情况下排土场边坡稳定性数值模拟研究［J].长江科学院院报，2015（02）：44-47.

［2]郑跃军.乌鲁瓦提水利枢纽大坝渗流监测与分析［J].水利规划与设计，2015（03）：60-63.

［3]雷子祥.复杂地质条件下的洞挖施工技术［J].水利技术监督，2015（03）：87-91.

［4]李宗坤，陈丽刚，韩立炜.基于A B A Q U S渗流与应力耦合的边坡稳定性分析［J].人民黄河，2011（02）：103-104.

［5]原先凡，邓华锋，蔡健，等.库水位升降条件下邻近坝址区堆积体边坡稳定性分析［J].水电能源科学，2013（09）：135-137.

［6]白阳.挡渣坝坝体渗流分析与施工设计［J].水利规划与设计，2014（11）：80-82.

［7]来金生.巩哈泉一库渗漏分析及治理方案［J].水利技术监督，2013（06）：51-53.

中图分类号：T V 871

文献标识码：B

文章编号：1672-2469（2016）06-0122-04