基于COMSOL Multiphysics的高面板堆石坝应力变形分析

邓七生

(鄱阳县水利局 江西上饶 333100)

混凝土面板堆石坝因其对环境适应性强、施工方便、造价相对较低、抗震性能好等优点，在世界范围内得到广泛应用。大坝变形是面板堆石坝的一个严重问题，因此对混凝土面板堆石坝的应力变形进行分析具有重要意义[1]。目前，对混凝土面板堆石坝的应力和变形分析进行了一些研究。黄向志采用三维非线性静力有限元法分析坝体和防渗体的应力和变形。李福忠等采用邓肯E-B非线性本构模型计算和分析了面板堆石坝的有限元应力和变形。张双梅等人使用ANSYS软件和有限元法对大坝进行了模拟和分析[2]。刘红等人使用ABAQUS 建立了有限元模型，以分析和研究大坝位移、大坝应力和板应力。李斌通过数值模拟分析了复杂覆盖层上面板堆石坝的应力变形。部分研究对象仍然是100 m 以下的面板堆石坝，随着施工技术的不断进步和完善，许多面板堆石大坝的高度已经超过100 m。因此，分析高面板堆石坝的应力和变形具有重要意义。该文采用Duncan-Zhang E-B 非线性弹性双曲线模型，利用有限元软件COMSOL Multiphysics 建立模型，利用MATLAB函数进行二次开发，对某水库混凝土面板堆石坝的应力和变形进行了分析[3]。该水库流域面积1 915.00 km2，总库容6.08 亿m3，水库死水位208.00 m，正常水位254.00 m。大坝为混凝土面板堆石坝，坝顶高程275.70 m，最大坝高102.20 m。

1 数值模拟模型

1.1 堆石体本构模型

在有限元数值模拟中，混凝土面板堆石坝坝体采用邓肯-张E-B 非线性弹性双曲线模型。非线性弹性双曲线模型的计算公式如下。

式（1）和（2）中，Pa为实际标准大气压力；K为弹性模量；n为弹性模量；σ1为第一主应力；σ3为第三主应力（围压）。

根据莫尔-库仑强度准则：

式（3）中，C和φ均为抗剪强度的指标。体积变形模量的计算公式为

式（4）中，Kb为体积模量系数；m为体积模量指数。摩擦φ随围压σ3变化的公式如下：

式（5）中，φ0为σ3=Pa条件下的内摩擦角。

1.2 接触模型

采用无厚度的Goodman单元模拟面板与衬垫之间的接触状态。Goodman单元以两侧对应节点的相对位移为变量，可以很好地模拟接触面的位错、滑移或张开，并可以考虑接触变形的非线性特征[4]。

2 初始应力场计算

根据以下方法确定新填充层的初始应力状态：

式（6）和（7）中，σ1和σ3分别为新填充层接合处的主应力和次主应力；γ为新填土层的重力；h为元素质心土壤表面以下的深度。式（8）中，K0为土壤的静态侧压力系数；φ为材料的内摩擦角。

3 基本方程的解

采用分步加载法模拟施工过程中的实际应力应变状态。基本方程采用荷载分级中点增量法求解。将载荷分为若干载荷增量，并对每个增量进行两次有限元计算。首先取荷载增量的一半，以求平均应力和相应的切线模量。然后，将其作为该类荷载增量的平均值，并求解在全荷载增量作用下的位移、应力和应变增量，该增量叠加在总位移、应力和应变上[5]。

4 计算软件及其二次开发

为了有效模拟和评估混凝土面板堆石坝施工和运行期间的结构安全状态，考虑到前处理和后处理的需要以及堆石坝非线性本构模型的二次开发，选择COMSOL Multiphysics 有限元计算软件进行研究。在邓肯拉伸E-B 模型中，堆石材料的切向弹性模量和体积变形模量通过加载和卸载问题计算，并且应力状态随每次加载而变化。基于COMSOL和MATLAB联合开发函数，使用MATLAB 函数计算每个加载步骤的切线弹性模量和体积变形模量。通过在COMSOL计算中自动调用MATLAB 函数，可以获得相应的切线弹性模量和体积变形模量[6]。对于初始应力场的计算，也基于COMSOL 与MATLAB 联合开发功能，并用MATLAB 进行相应的二次开发。

5 有限元模型

考虑盘石头水库混凝土面板堆石坝的实际断面结构和材料分区，建立了坝体的二维有限元模型。首先在CAD中建立超元模型（如图1所示），然后将CAD模型导入COMSOL Multiphysics 有限元软件，并采用极精细网格细分模式自动细分有限元网格（如图2 所示）。大坝底部固定。细分单元总数为20 792个。大坝填筑分为11个阶段。水负荷分为5个阶段。

图1 CAD模型图

图2 有限元网格

6 计算工况和参数

不同材料的混凝土和堆石体计算参数如表1所示。

表1 计算参数

7 结果分析

堆石体的水平和垂直位移云图以及主要和次要主应力云图如图3～图6所示。

图3 正常蓄水条件下的垂直位移和水平位移云图

图4 正常蓄水条件下的主要主应力和次要主应力云图

图5 校核洪水位条件下的垂直位移和水平位移云图

图6 校核洪水位条件下的主要主应力和次要主应力云图

大坝各种条件下堆石的最大位移和最大应力如表2所示。

表2 最大位移和最大应力

完工期的板挠度较小，为10.3 cm，出现在覆盖层顶部的高程处。在正常蓄水条件下，板的最大挠度出现在板的中间，值为34.2 cm，周边接缝的沉降为1.04 cm，周边缝的开口位移为1.18 cm。在设计洪水位条件下，板的最大挠度出现在板的中间，值为31.4 cm，周边接缝的沉降为1.22 cm，周边缝的开口位移为1.30 cm。在校核洪水位条件下，板的最大挠度出现在板的中间，值为53.9 cm，周边接缝的沉降为1.28 cm，周边缝的开口位移为1.36 cm。当水压作用在面板上时，面板法线方向的应力基本一致，约等于面板前水压的大小，其值小于混凝土的抗压强度。面板主要承受沿坝坡的拉应力。在正常蓄水条件下，面板上的最大拉应力为2.58 MPa，但周边接缝附近的位置除外，该位置出现在194.5 m高程的面板表面。从258.37 m 标高开始，标高以下的拉应力超过C25混凝土抗拉强度的设计值。在设计洪水位条件下，面板上的最大拉应力为2.86 MPa，但周边接缝附近的位置除外，该位置出现在194.5 m高程的面板表面。从261.92 m 标高开始，标高以下的拉应力超过C25 混凝土抗拉强度的设计值。在校核洪水位条件下，面板上的最大拉应力为3.14 MPa，出现在194.5 m高程的面板表面。从标高262.27 m 开始，标高以下的拉应力超过C25混凝土抗拉强度的设计值。

8 结语

该文基于COMSOL 与MATLAB 的联合开发函数，利用MATLAB 函数进行二次开发，模拟分析了不同工况下堆石坝的应力和变形特性。计算结果与实际相符，表明模型和计算方法科学合理，为相关工程提供了一种新的计算方法。有限元计算和分析表明，由于下游堆石区的软岩材料，大坝的沉降明显从大坝位移分布向下游倾斜。大坝沉降不大，但上下游大坝沉降差较大。周边接缝在正常范围内变形。大坝中没有明显的剪切破坏区。根据有限元计算结果，在水压作用下，面板的小主应力会出现超过混凝土抗拉强度的拉应力。因此，应加强该小组的日常维护和监测。