柏叶口水库面板堆石坝变形数值分析

薛志平 ，孙一清 ，张 震

（1.山西省柏叶口水库建设管理局，山西 吕梁 033000；2.河海大学水利水电学院，江苏 南京 210098）

0 引言

柏叶口水库位于山西省吕梁市交城县会立乡文峪河干流上，控制流域面积 875 km2。坝址区河谷底宽约为 120 m，河谷形状呈 U 形，主要建筑物为面板堆石坝，坝轴线走向为 N 63.23°E。坝顶高程为1 138.30 m，最大坝高为 88.30 m，是山西省最高的大坝，其安全对于吕梁市和汾河中下游的防洪安全具有十分重要的意义。

现场检查发现大坝坝顶及上游面板存在多条裂缝，由于该坝施工速度比较快，前期监测系统不够完备，更新后的监测系统运行时间不长，监测资料匮乏。为清楚地分析大坝变形规律，全面地分析快速施工对大坝安全的影响，有必要通过数值手段全面了解大坝各种不利工况下的变形情况[1-4]。为此结合柏叶口水库大坝实际工程资料，仿真模拟大坝堆石料的填筑、混凝土面板浇筑及蓄水过程，建立考虑具体施工过程分级加载的大坝三维有限元分析模型，运用 COMSOL Multiphysics 数值分析软件获取坝体在竣工期、正常蓄水位、设计及校核洪水位工况下的变形特性。

1 有限元分析模型

1.1 建模过程分析

要建立有限元分析模型，先要在 Auto CAD 中建立准确的混凝土面板堆石坝几何模型。运用 Auto CAD 绘图软件绘制大坝三维几何模型，坐标系的选取如下：X 轴沿坝轴线指向右岸方向为正，零点取在坝轴线左岸处；Y 轴指向水流方向为正，零点取在坝横右 0 + 000.00 m 断面处；Z 轴垂直指向上方为正，与高程一致。为提高计算结果的准确性，计算模型在 X，Y，Z 各轴的边界范围如下：在竖直方向上，底部边界为 Z = 913.5 m 的截面；在顺河向上，上游的边界为 Y = -306.718 m 的截面，下游边界为Y = -298.293 m 的截面；在坝轴线方向上，左边界为 X = 486.6 m 的截面，右边界为 X = -176.6 m 的截面。为避免导入 COMSOL Multiphysics 仿真计算时畸变单元的产生，对大坝几何模型进行优化，避免扭曲面的产生[5-7]。优化过程中不断使用 Live Link功能，连接 Auto CAD 与 COMSOL Multiphysics，实现边修改边导入的数据同步。柏叶口大坝坝体三维有限元网格模型如图 1 所示，剖分后得到 336 632 个单元。

图 1 柏叶口大坝三维有限元网格模型（单位：m）

1.2 荷载分级

真实模拟施工荷载是混凝土面板堆石坝有限元分析中至关重要的一步，可以体现结构体本身随施工过程的变化情况，从而更好地反映材料的非线性，使计算结果更符合实际。通过工程资料可知，坝顶高程为 1 138.30 m，正常蓄水位为 1 133.00 m。在混凝土面板堆石坝的施工中，坝体的填筑过程是分层碾压，具体施工程序如下：1）混凝土趾板浇筑；2）坝体次堆石区高程 1 098.00 m 以下填筑；3）坝体主堆石区高程 1 098.00 m 以下及过渡层、垫层填筑；4）坝体堆石区高程 1 098.00 m 至坝顶及过渡层、垫层填筑；5）面板浇筑；6）铺盖、防浪墙及下游护坡填筑。

初次蓄水时，根据蓄水安全鉴定报告与坝体结构设计和度汛要求，参考已建工程经验，水荷载是分级施加的，结合坝体网格单元的剖分情况，具体加载方式如下：1）水位 1 061.9～1 080.0 m 蓄水；2）水位 1 080.0～1 100.0 m 蓄水；3）水位 1 100.0～1 128.0 m 蓄水；4）水位 1 128.0～1 133.0 m 蓄水到正常蓄水位。

计算工况中还应考虑运行期水位升至设计和校核洪水位的工况，选取竣工期未蓄水、正常蓄水位、设计和校核洪水位 4 种计算工况依次分析面板堆石坝变形分布。

1.3 参数选取及计算工况确定

1.3.1 材料参数的选取

坝体填筑料采用邓肯-张（E-B）模型，相关参数根据现场试验参照类似工程确定，如表 1 所示。基岩地层和混凝土采用线弹性模型，计算参数依据表 2 所示资料。

1.3.2 沉降结果的分析

大坝外部沉降从 2011 年 6 月开始观测（1 078，1 098，1 118 高程观测点），2012 年 6 月坝顶高程开始观测，测点分布如图 2 所示。通过观测资料分析，最大沉降发生在 W12观测点，沉降量为 70.32 mm，坝顶最大沉降发生在 W9观测点，沉降量为 45.63 mm，沉降遵循逐渐下沉规律，没有明显的突变点。观测统计资料[8]如表 3 所示。

可见，通过 Auto CAD 三维模型导入 COMSOL Multiphysics 进行网格剖分，选取合适的材料及材料参数、本构模型进行模拟，通过数值模拟出的沉降计算值与实际观测的沉降值对比，发现计算值偏小。原因分析如下：1）未将温度及时效对坝体沉降的影响考虑进数值模拟中；2）面板缝与周边缝建模考虑不够全面；3）水渗入面板造成的沉降改变等诸多因素造成。

表 1 非线性材料本构关系的计算参数表

表 2 线弹性材料应力应变关系计算参数表

图 2 柏叶口大坝外部测点分布图

表 3 大坝外部沉降观测资料统计表 mm

2 坝体变形分析

分析过程中，最先模拟基岩的初始地应力场。随着自重和水压力等的施加，堆石料和面板的位移及应力不断变化。到最后加载完成，会较稳定地呈现一定的规律[9-11]。

计算结果中位移的正负与坐标正负一致，在分析该混凝土面板堆石坝的变形时，选取大坝典型断面河床处 1-1 断面（坝 0 + 180.00 断面），左岸 2-2断面（坝 0 + 060.00 断面），右岸 3-3 断面（坝 0 +245.00 断面），分析这 3 个典型断面在竣工期、正常蓄水位、设计及校核洪水位 4 种工况下的变形。

2.1 竣工期变形分析

2.1.1 河床 1-1 断面

坝体河床处 1-1 断面竣工期位移分布如图 3 所示。由于断面上下游相对呈现大致对称结构，计算结果显示，竣工期的位移也具有较好的对称性。由图 3 a 可得，竣工期未蓄水时，坝体上下游最大位移约为 4.50 mm；由图 3 c 可得，坝体的最大沉降约为55.90 mm，由于存在铺盖层（1 A，1 B），沉降最大值出现在坝体的中部偏上 2/3 处。

图 3 坝体河床 1-1 断面在竣工期的位移（单位：mm）

2.1.2 左岸 2-2 断面

竣工期坝体左岸 2-2 断面位移分布如图 4 所示。由图 4 a 可得，竣工期未蓄水时，左岸坝体上下游最大位移约为 2.80 mm；由图 4 c 可得，坝体的最大沉降约为 34.10 mm，由于存在 1 A，1 B 铺盖层，沉降最大值出现在坝体的中部偏上 2/3 处。

图 4 坝体左岸 2-2 断面在竣工期的位移（单位：mm）

2.1.3 右岸 3-3 断面

竣工期坝体右岸 3-3 断面位移分布如图 5 所示。由图 5 a 可得，竣工期未蓄水时，右岸坝体上下游最大位移约为 2.88 mm；由图 5 c 可得，坝体的最大沉降约为 37.50 mm，由于存在 1 A，1 B 铺盖层，沉降最大值出现在坝体的中部偏上 2/3 处。

2.2 正常蓄水位变形分析

正常蓄水位工况下，坝体河床处 1-1、左岸 2-2及右岸 3-3 等断面中，坝体上下游最大位移出现在河床 1-1 断面，且大于左右岸断面，约为 7.30 mm。由于存在 1 A 和 1 B 铺盖层，3 个断面的沉降最大值均出现在坝体的中部偏上 2/3 处，且河床 1-1 断面的坝体沉降最大，约为 60.30 mm。

图 5 坝体右岸 3-3 断面在竣工期的位移（单位：mm）

2.3 设计洪水位变形分析

设计洪水位工况下，坝体河床处 1-1、左岸 2-2及右岸 3-3 等断面中，坝体上下游最大位移出现在河床处 1-1 断面，且大于左右岸断面，约为 7.63 mm。由于存在 1 A 和 1 B 铺盖层，3 个断面的沉降最大值均出现在坝体的中部偏上 2/3 处，且河床 1-1 断面的坝体沉降最大，约为 61.30 mm。

2.4 校核洪水位变形分析

校核洪水位工况下，坝体河床处 1-1、左岸 2-2及右岸 3-3 等断面中，坝体上下游最大位移同样在河床处 1-1 断面，且大于左右岸断面，约为 7.75 mm。由于存在 1 A 和 1 B 铺盖层，3 个断面的沉降最大值均出现在坝体的中部偏上 2/3 处，且河床 1-1 断面的坝体沉降最大，约为 61.60 mm。

3 结语

本研究用非线性有限元法的理论及计算方法、堆石料的本构模型、施工逐级加载方式对柏叶口混凝土面板堆石坝变形进行了研究，解决了以下问题：1）导入 COMSOL Multiphysics 商用软件，使用自由四面体网格，进行了模型的有限元网格剖分；2）竣工期垂直方向坝体的最大位移（沉降）为 55.60 mm，出现在坝体的中上部分，约占最大坝高的 0.33%，向下游的水平位移大于向上游的水平位移。正常蓄水位下，坝体在竖直方向上，坝体的最大沉降为 60.30 mm。随着水位升高坝体沉降不断增大，设计和校核洪水位相对应的沉降量为 61.30和 61.60 mm。由此可知柏叶口混凝土面板堆石坝整体安全，本工程施工速度快对整体变形规律改变不大。下一步将引入考虑蠕变的本构模型，进行深入研究，从而为动态变形监控指标奠定基础。