基于分区有限元与块体界面元混合算法的重力坝深层抗滑稳定分析

2020-06-11李同春

中国农村水利水电 2020年5期

王凯，李同春，程井

(1. 河海大学水利水电学院，南京 210098；2. 河海大学农业工程学院，南京 210098； 3. 河海大学水资源高效利用与工程安全国家工程研究中心，南京 210098)

0 引言

近几十年来，水电事业得到了蓬勃发展，重力坝作为重要坝型得到了广泛应用。随着资源不断开发利用，不可避免地需要在地质条件复杂的岩基上修建重力坝。由于岩体中普遍存在节理、夹层、断层等软弱结构面，特别是缓倾角的结构面在重力和水推力的作用下，有可能形成滑移通道，导致大坝失稳破坏，因此重力坝深层抗滑研究具有实际的工程意义。

岩体中存在的软弱结构面将岩体分割成不连续体，节理、夹层、断层等软弱结构面构成了不连续结构面。坝体失稳往往伴随着岩体中这些不连续结构面地张开、闭合和滑移等现象，属于典型的局部位移非线性的接触问题，因此深层抗滑稳定分析的关键在于如何有效地模拟岩石裂隙的力学行为[1,2]。近年来，主要采用刚体极限平衡法和基于连续介质理论的有限单元法[3-5]进行重力坝深层抗滑稳定分析。刚体极限平衡法由于其概念清晰，计算简便，被大量的工程采用，积累了丰富的工程经验，同时具备一套成熟的设计准则，是一种十分常见的稳定分析方法。但刚体极限平衡法为了求解方便，在力学分析上作了较大的人为假定，与实际情况偏差较大。随着计算机技术的发展，以有限单元法为代表的数值计算方法得到广泛运用。然而有限元法属于连续介质方法，必须满足变形相容性条件，无法真实模拟裂缝地张开滑移。为了解决此类局部位移非连续问题，李同春等[6]提出了分区有限元与块体界面元混合算法，可以较好地模拟不连续结构面的闭合、滑移和张开，同时仅将非线性迭代收缩到接触界面上进行，大大提高了计算效率。本文在此方法基础上，针对某重力坝深层抗滑稳定问题，结合强度折减法，利用接触点对的状态来模拟深层滑动面的力学行为，以失效点对首次贯通作为失稳判据，求解得到深层抗滑稳定安全系数。该方法将重力坝深层抗滑问题考虑为接触问题，为今后的深层抗滑稳定研究提供了一种新的求解思路和方法。

1 深层抗滑稳定分析方法

1.1 刚体极限平衡方法

刚体极限平衡法的求解思路是根据岩基中存在的软弱结构面假设可能的深层滑动面，一般的假设滑面形式包括单滑面、双滑面以及多滑面。然后将滑裂体看作不变形的刚体，不考虑力矩平衡，只通过力的平衡求出滑面上的滑动力和抗滑力，将抗滑力与滑动力比值作为抗滑稳定安全系数。最常见的滑面形式为双滑面，根据设计规范采用等安全系数法，用试算法或迭代法求得整体抗滑稳定安全系数K及抗力Q值。该方法应用十分广泛，但也存在很多缺点，比如计算时需首先假定滑动面，通过试算确定最危险滑面，但对于地质条件复杂的工程，试算工作量巨大；并且将岩体考虑为刚体，不能考虑其受力变形所带来的影响，极限状态与允许的工作状态有较大的出入[3]；不能确定滑面上的应力分布，不能探索破坏的机理及其变化发展过程；等等。但它是经过大量的工程实践验证过的，具有很高的可信度，被规范[7]所采用。

1.2 有限单元法

有限单元法可以考虑岩体的非线性弹塑性本构关系，以及变形对应力的影响，结合强度折减法，可以模拟岩体渐进破坏到整体剪切破坏过程。通过有限元计算得到坝基的塑性区图，基于最大最小值理论及有限元的最小势能原理，可以近似认为塑性区图上塑性应变值最大的点的连线(对平面问题)即为临界滑动面[8]，即可得到重力坝的失稳模式。为了获得抗滑稳定安全系数，一般需要结合失稳判据确定。常用的失稳判据包括特征点位移突变和塑性区贯通[5]，但这两种失稳判据均存在较大的人为性，比如位移突变拐点常常不明确，得到的抗滑稳定安全系数为一个区间，以及塑性值达到多大才算贯通一直没有统一标准，判据还存在较大争议[9]。同时，有限元法基于连续介质模型，需满足位移协调条件，无法模拟岩石裂隙地张开滑移，不能真实反映滑动面上的力学特性，而且弹塑性有限元法还存在计算不容易收敛、接近临界失稳状态无法求解等问题，因此有必要寻找一种既能反映岩体的连续-非连续位移特征，又兼具较高计算效率的分析方法。

1.3 分区有限元与块体界面元混合算法

分区有限元与块体界面元混合算法将求解系统分解为若干块体与非连续界面，每个块体作为系统有限元的一个分区，将不连续结构面视为接触界面，在接触界面上生成接触点对，通过接触点对的闭合、滑移、张开模拟界面的破坏情况。以块体内的位移作为分区有限元求解的基本变量，以块体形心的刚体位移作为刚体运动变量，通过界面上力与变形之间的关系和刚体平衡方程，形成以界面上的接触内力和块体形心刚体位移为混合变量与界面结点总位移进行迭代求解的混合方程[10,11]。具体的方程推导见文献[6]。

分区有限元与块体界面元混合算法计算结果只能得到接触点对的状态、结点的应力应变和块体的刚体位移。因此将混合算法与强度折减法相结合，不断折减滑动面上的抗剪断强度参数，以接触点对首次全部失效作为失稳判据，将此时的折减系数作为抗滑稳定安全系数。

针对某典型重力坝工程的深层滑动问题，采用分区有限元与块体界面元混合算法进行分析。基于该混合算法的思路，首先需要将求解系统分解为若干块体和非连续界面，需预先确定可能的深层滑动面位置，故求解该重力坝深层滑动问题的主要步骤如下：①首先根据岩基的地质构造建立有限元模型，采用弹塑性有限元强度折减法确定最可能的滑动面。②修改初始有限元模型，在滑动面位置建立接触点对。③结合强度折减法，采用分区有限元与块体界面元混合算法对修改后的模型进行计算，以接触点对全部失效作为失稳判据，求出重力坝抗滑稳定安全系数。

2 典型工程应用

2.1 工程概况

图1 挡水坝段剖面图(单位：m)Fig.1 Retaining dam section profile

2.2 材料参数及计算工况

混凝土坝体分三期施工，一期、二期和三期分别采用C30、C25、C20混凝土进行浇筑。岩基主要为大理岩夹云母石英片岩。在坝踵附近发现有一条倾向下游的夹层，编号f511，宽度约为15cm。坝体、地基及软弱夹层的物理力学参数如表1所示。

选取正常蓄水位工况进行分析计算。重力坝整体结构所受的荷载主要包括坝体自重、泥沙压力、上下游水压力以及地基受到的渗透压力。初始地应力仅考虑为地基自重所引起的地应力场。上下游正常蓄水位高程分别为2 702.0和2 544.32 m，水重度取9.81 kN/m3。淤沙高程2 595.35 m，浮容重7.8 kN/m3，内摩擦角12°。对于渗透压力，考虑坝体为不透水材料，仅考虑地基受渗流力作用，先根据各分区渗透系数等计算参数，进行稳定渗流场的计算，得到渗流体积力，然后将计算得到的渗透压力作为结点荷载施加到地基中[12]。

2.3 数值计算

2.3.1 弹塑性有限元强度折减法

建立重力坝初始有限元模型，如图 2所示。坝基范围分别向上下游、底部延伸2倍坝高。整体采用4结点单元进行离散，共计2488个节点，2368个单元。坝体材料采用线弹性本构模型，软弱夹层采用MC屈服准则，地基岩石采用DP屈服准则。坝基两侧采用法向约束，底部采用固定约束。

(1)可能滑动面搜索。首先进行坝基渗流场计算，考虑防渗帷幕的作用，上下游河谷分别施加相应水头，其他作为不透水边界，得到渗透压力分布如图 3所示。将渗透压力作为节点荷载施加到地基中，并施加其他荷载，不断折减地基材料的抗剪断强度参数，塑性区不断扩展并最终贯通。根据之前的分析，将塑性区图上塑性应变值最大的点连接起来即为最危险滑动面，求得的可能滑动面标记如图 4所示。可以看到软弱夹层自然构成深层滑移的主滑面，滑移剪出面角度约20°左右，从坝趾处附近剪出。滑面为典型的双滑面形式。

图3 岩基渗透压力云图(单位：kPa)Fig.3 Pore-pressure nephogram of rock bed

图4 岩基塑性区通道示意图(折减系数Fs=2.4)Fig.4 Plastic zone of rock bed (reduction factor Fs=2.4)

(2)基于弹塑性有限元强度折减法的抗滑稳定安全系数。选取坝顶和坝趾为特征点，图 5为特征点水平位移随折减系数的过程线。从曲线图可以看出，两特征点位移拐点皆不明显，无法直接准确获得抗滑稳定安全系数，塑性值的选取也没有统一的标准，该方法的失稳判据在确定抗滑稳定安全系数数值时具有不可避免的人为因素。对于本工程，结合位移拐点和塑性区贯通判据，综合判断抗滑稳定安全系数约在2.3～2.5之间。

图5 特征点位移过程线Fig.5 Curve of x displacement-reduction factor for feature points

2.3.2 分区有限元与块体界面元混合算法

基于分区有限元与块体界面元混合算法，在搜索得到的可能滑动面上建立接触点对，修改后的有限元模型见图6。模型一共包含2 931个节点，2 804个单元，滑动面上一共生成53对接触点对。坝体和部分地基作为滑动体，滑动面以下的地基作为基础。计算时仅考虑滑动面上的强度特性，滑动块体和基础块体按弹性体考虑。采用分区有限元与块体界面元混合算法对系统进行计算，得到接触点对的状态。接触点对所处状态定义有3种：0为张开，1为闭合，2为滑移，计算时认为接触点对处于张开或滑移状态判定为失效。折减系数从1.00开始逐渐递增，初始每次递增0.10，在接触点对几乎全部失效，即将贯通时，变为每次递增0.01，统计不同折减系数下的接触点对的状态，结果如表2和图7所示。可以看出随着折减系数不断增大，接触点对的失效比例逐渐上升。从图8失效点对的位置可以发现，软弱夹层抗剪断强度参数低，其上的接触点对首先失效。折减系数不断增大时，剪出面上的接触点对开始失效，并最终全部失效贯通，此时判定大坝失稳。

图6 重力坝修改后的有限元模型(含滑动面)Fig.6 Modified FEM meshes of gravity dam (sliding surface included)

序号折减系数接触点对失效组数/接触点对总组数11.0024/5321.3025/5331.4030/5341.5037/5351.6043/5361.8046/5371.9049/5382.0150/5392.1051/53102.1353/53

图7 接触点对失效比例随折减系数变化过程线Fig.7 Curve of failure ratio of nodal pairs-reduction factor

图8 不同折减系数下接触点对状态图(Δ为失效)Fig.8 States of nodal pairs of different reduction factor (Δ represents failure)

根据统计结果，当折减系数为2.13时，接触点对首次全部失效，认为重力坝失稳，求得抗滑稳定安全系数为2.13。结果与弹塑性有限元强度折减法求得结果相比，数值上十分接近，证明了此方法的正确性。但前者比后者结果要小，原因在于分区有限元与块体界面元混合算法很好地考虑了滑动面上的力学行为，模拟了裂缝的滑移与张开。而有限单元法基于连续介质模型，必须满足位移协调性，滑动面无法张开或滑移，无法及时进行应力重分配，导致法向压应力相对增加和切向合力减小[1]，滑动面的安全裕度较大。