摘要: 堤坝边坡的稳定性是保障堤坝安全运行的重要因素之一。堤坝边坡的不稳定可能导致堤坝破坏，从而引发严重灾害，造成人员伤亡和财产损失。因此，保证堤坝边坡的稳定性对于保障人民生命财产安全和维护社会稳定具有重要意义。采用有限元COMSOL对堤坝进行建模，建立大坝边坡受各种载荷的二维有限元分析模型，从而进行研究分析，安全系数(FOS)以及重力和静水压力的影响也是预测堤坝边坡稳定性的必要参数。模拟后给出堤坝中的压力水头、失稳前的等效塑性应变和堤坝失稳前的位移大小以及最大位移与FOS的关系示意图。为了控制地表沉降、变形和滑移等几个重要的影响因素提供了重要的参考依据，对堤坝边坡稳定性的研究开展具有重要意义。
Abstract:
The stability of the embankment slope is one of the important factors to ensure the safe operation of the embankment. The instability of the embankment slope may lead to the failure of the embankment, which will lead to serious disasters, resulting in casualties and property losses. Therefore, it is of great significance to ensure the stability of the embankment slope for the safety of people’s life and property and the maintenance of social stability. The finite element COMSOL is used to model the dam slope, and the two-dimensional finite element analysis model of the dam slope subjected to various loads is established for research and analysis. The safety factor (FOS), the influence of gravity and hydrostatic pressure are also necessary parameters to predict the stability of the dam slope. After the analysis, the pressure head, the equivalent plastic strain before the instability, the displacement before the instability of the dam and the relationship between FOS and the maximum displacement are given. In order to control the surface subsidence, deformation and slip of several important factors to provide an important reference basis for the study of embankment slope stability has great significance.
1. 引言
堤坝作为一种重要的水利工程，被广泛应用于水资源的调节和利用。堤坝边坡稳定性问题在我国极具现实意义，随着我国经济的快速发展，我国在地形复杂地区建设的大坝数量越来越多。堤坝边坡不稳定可能导致堤坝破坏，从而引发严重灾害，造成人员伤亡和财产损失。因此，保证堤坝边坡的稳定性对于保障人民生命财产安全和维护社会稳定具有重要意义。堤坝边坡失稳是指堤坝边坡所承受的外力超过了其所能承受的极限，从而导致边坡出现破坏和滑坡等现象。
堤坝失稳前后准确地确定坝体的应力状态和塑性区以及堤坝失稳前的滑弧，将为堤坝建筑的过程以及对之后部分地方的加固等提供了有效的支持。用于边坡稳定性分析的方法主要包括极限平衡方法、极限分析方法和有限元分析方法 [1] [2] [3] 。张鲁渝 [4] 也用有限计算数量级小于10作为边坡失稳的判据，由于采用的有限元程序、收敛容许值和迭代次数的差异，导致计算的安全系数不同。栾茂田 [5] 将抗强度折减基本概念、弹塑性有限元分析原理与计算机结果图形实时显示技术相结合，提出以广义塑性应变及塑性开展区作为边坡失稳的评判依据，可以比较准确地预测边坡潜在破坏面的形状和位置及相应的安全系数。国内学者宋二祥 [6] 以一座土坝为例，利用强度折减系数与土坝顶点水平位移关系曲线标准得出与Bishop法计算的最小安全系数相近的结果。尽管边坡失稳判据不一，但不可否认通过强度折减来分析结构的稳定性，直到临界状态为止。同时有限元强度折减法能分析各种复杂的边坡工程，不需要假设滑动面，考虑了土体的弹塑性本构关系，满足力的平衡方程，计算结果更可靠。
用于模拟仿真和分析堤坝边坡的软件主要有ANSYS、ABAQUS、COMSOL、FLAC、PFC等。大数据统计发现近些年来，采用有限元仿真分析关于堤坝安全稳定性的越来越多，比如说渗流、抗震和边坡等问题备受关注 [7] [8] [9] 。由于COMSOL Multiphysics软件植根于数学物理方法，直面数学物理方程，从最底层的理论出发，揭示物理现象的本质，适于做多物理场仿真分析。它的全耦合求解方法给物理场的耦合问题提供了新的、便捷的解决方案。根据以上特点，本文基于COMSOL Multiphysics有限元软件中的固体力学模块对某一堤坝边坡稳定性进行分析，以了解堤坝在受载荷以及环境作用下的反应 [10] 。
2. 有限元分析模型的建立
2.1. 建立模型
首先，建立几何模型堤坝高度为12 m，长度L1、L2和L3分别为24 m、5 m和24 m。水位为10 m，可能渗流高度为4 m。堤坝总宽度为L1 + L2 + L3。其次，设置材料参数，线弹性材料为各向同性，土体力学参数分别为杨氏模量表达式输入100 Mpa，泊松比为0.4，密度为2000 kg/m3，孔隙率为0.3。最后，设置边界条件，堤坝边界设置除了水库水位方向和渗流高度方向的斜面，其余的均为边界约束。设置两个压力水头，压力水头1为水库水位方向的斜面，压力水头2为渗流高度的斜面。本例中使用平面应变近似在二维模式下对堤坝进行建模。建立模型如图1所示。
Figure 1. Schematic diagram of dam model
图1. 堤坝模型示意图
2.2. 构建网格
Figure 2. Internal network structure diagram of the dam
图2. 大坝内部网络结构图
本例重力和静水压力的影响也包括在内。莫尔–库伦模型的材料属性根据安全系数参数FOS进行参数演化，增加FOS的参数，从而降低了每个参数步骤中的土壤强度。对于1.915以上的参数值，模型不收敛，标明堤坝失稳。
如图2所示，根据网络图分析出网格被自由三角形划分，同时也可以看出，靠近堤坝中间部分的网格很密集，其余区域网格稀疏。
2.3. 研究计算
构建网格后，为了说明孔隙压力和重力载荷对堤坝稳定性的影响，我们建立了三个稳态研究。在第一个研究中，计算达西定律来获得孔隙压力分布图。在第二个研究中，用这种静水压力和重力载荷模拟堤坝受力。在第三个研究步骤中，通过预应力和预应变节点，将第一个研究中产生的孔隙压力和第二个研究中重力载荷产生的预应力考虑在内。在第三个研究中增加了莫尔–库伦准则，研究重力和可变孔隙压力共同作用下的土壤弹塑性破坏。经过计算，得出以下结果，如图3所示，图3表示经过上述三种稳定分析，考虑其堤坝本身的问题，造成的失稳前位移大小，采用三维云图表示，更清晰直观地表示各部位发生位移的大小。
Figure 3. Displacement before dam instability
图3. 堤坝失稳前位移
经上述稳定性分析后，堤坝的压力水头如图4所示：在浸没壁上，它从0 m到10 m不等，而在渗流面上为0 m。正压力水头表示正孔隙压力，表示饱和土，而非饱和土用零压力水头表示，其中零压力水头线是将饱和土和非饱和土分开的潜水面位。
等效塑性应变是用来确定材料经强化后屈服面的位置的物理量，其计算公式与计算等效应力的第四强度理论计算公式是一样的，只是把应力改成应变。在试验中得到的是单向拉伸时材料的弹塑性应力应变曲线，而实际中结构的应力状态往往是空间应力状态，屈服时其也是空间屈服，通过计算等效塑性应变我们可以找到其塑性屈服后应变状态对应的等效应力和此时结构的空间应力状态。
在我们分析堤坝失稳前后的状态时，失稳前的等效塑性应变显示不同的模式，至此塑性区的形状和位移至关重要，图5显示了堤坝失稳前的破坏机理，该区域出现在堤坝的塑性区域一侧，从图中我们也可以看出，该区域后续应该采取相关的措施加以改进。
Figure 4. Pressure Head Change on the Isobaric Line
图4. 等压线上压力水头变化
Figure 5. Equivalent plastic strain before dam instability
图5. 堤坝失稳前的等效塑性应变
堤坝的一侧由于不稳定性，导致土体发生滑移现象。图6中的箭头表示土颗粒的位移方向，由于堤坝下边界为固定约束，因此堤坝右下角土体不会出现滑移。
Figure 6. Sliding arc before dam instability
图6. 堤坝失稳前的滑弧
经过以上的步骤，总结并计算其最大位移与FOS的关系。经计算，结果可以得出图7所示，最大位移是在FOS = 1.9附近显著增加，这标明堤坝开始失稳。
Figure 7. Relationship between maximum displacement and FOS
图7. 最大位移与FOS的关系
3. 结论
采用有限元COMSOL Multiphysics软件模拟堤坝边坡滑移，建立堤坝变形的二维有限元分析模型，从而进行分析。利用有限元强度折减法计算得到安全系数，迭代终止时迎水坡发生的塑性贯通区，同时得出堤坡发生的总位移。计算结果与物理模型试验、数值模拟结果对比数值上会有差异，规律较为合理，证明了考虑静水压力和重力作用下，用有限元强度折减法计算的可靠性和可行性，对深入研究更复杂实际工程、多种外界条件下的多物理场合分析的数值模拟有一定的参考意义。得出的主要结论为：在建造堤坝的时候，堤坝由于压力水头的作用，堤坝的迎水坡会产生塑性区域，因而发生堤坝失稳产生位移。
NOTES
*通讯作者。