前坪水库砂砾石黏土心墙坝三维流固耦合分析

历从实，武颖利，皇甫泽华，张兆省，韩 艳

（1.河南省前坪水库建设管理局，河南 郑州450003；2.河南省水利第一工程局，河南 郑州450000；3.南京水利科学研究院，江苏南京210098）

流固耦合作为自然界普遍存在的一种现象,其蕴含着复杂的物理现象和深刻的数学原理［1］。近年来,国内外学者在很多领域对流固耦合数值模拟和算法进行了探索研究,并取得了一定成果［2-6］。在水利工程领域,王晓玲等［7］对考虑流固耦合作用的拱坝抗滑稳定可靠度进行了研究；李冰冰等［8-9］对考虑流固耦合作用的尾矿坝、土石坝渗流稳定进行了研究；刘占涛等［10-11］对基于流固耦合作用的深厚覆盖层和复杂坝基土石坝稳定进行了研究；陈文元等［12］对考虑流固耦合作用的坝体地震动力特性进行了研究；唐克东等［13］对不同开度弧形闸门流固耦合进行了数值模拟,取得了一定的成果。砂砾石黏土心墙坝作为土石坝的一种,相关流固耦合模拟研究文献记载不多。本研究针对前坪水库工程砂砾石黏土心墙坝的特点,结合地质、设计和相关试验资料,采用南水双屈服面本构模型,结合三维流固耦合计算方法,模拟坝体的填筑过程、施工后至蓄水前、蓄水过程以及运行期的坝体变形规律、应力分布规律和工作性态。

1 前坪水库工程坝体结构特点

前坪水库工程是国务院确定的172项节水供水重大水利工程之一。坝型为砂砾石黏土心墙坝,设计坝高90.3 m,坝顶长810 m,坝顶宽10 m,上游坡比1∶2.0～1∶2.5,下游坡比1∶2.0。黏土心墙顶宽为4 m,上下游坡比均为1∶0.3,紧贴心墙外侧为2 m宽粗砂反滤层,粗砂反滤层外侧为2 m宽碎石反滤层,碎石反滤层外侧为砂砾石坝壳料。右岸坝肩坡比为1∶0.75,左岸坝肩坡比为1∶1。

2 三维有限元模型及计算参数

采用基于ABAQUS有限元软件二次开发工具UMAT数据接口编写的南水双屈服面本构模型进行计算。根据设计图纸,建立三维有限元模型。模型深度自建基面向下取80 m,两侧自坝底边线向外延伸350 m作为截断边界。模型底部施加三向约束,四周施加法向约束。坝体三维有限元模型如图1所示。模型共326 648个单元、345 836个节点,单元类型为C3D8P,在反滤料-心墙、防渗墙-坝体间设立薄层单元替代接触面单元。

计算模拟分79个计算步。第1步为地应力平衡步,建立初始地应力场,消除覆盖层初始沉降位移影响；第2～72步为土石坝施工填筑过程；第73～78步为土石坝蓄水至校核水位的过程；第79步模拟校核水位下持续运行10 000 d,研究坝体的工后变形。采用流固耦合方法对填筑施工期、竣工60 d后、蓄水期及蓄水运行20 a四种工况土石坝的受力变形特性进行计算分析。

图1 坝体三维有限元模型

南水双屈服面模型兼顾了邓肯模型和剑桥模型的优点,考虑了土体的剪胀特性,且对各类土体适应性强,砂砾石坝壳料、心墙料、高塑性黏土、砂砾石覆盖层均采用南水双屈服面弹塑性本构模型［14］进行计算,对于基岩、混凝土防渗墙、防渗帷幕采用线弹性模型进行计算。

通过试验和地质勘查资料获取计算参数,见表1。

表1 坝体各区域主要计算参数

3 坝体应力变形计算分析

3.1 施工期坝体应力变形特性

图2为坝体0+280断面施工期应力应变特性图。

由图2可知,施工期0+280断面最大沉降量为118.8 cm,发生在约1/2坝高位置,最大水平位移为16.3 cm,发生在心墙与上游侧反滤料交接面靠近坝顶附近,方向指向下游侧,最小水平位移为-17.3 cm,发生在最大水平位移的对称位置,方向指向上游侧。从大小主应力分布规律来看,施工期坝体大小主应力均为压应力；大小主应力在心墙与反滤层交界位置有明显跌落,心墙内部应力拱比较明显。施工期坝体整体应力水平不高,仅在混凝土防渗墙与高塑性黏土交界位置附近因变形不协调导致应力劣化而出现较高的应力水平,坝体绝大部分位置应力水平较低。

3.2 竣工后60 d坝体应力变形特性

坝体0+280断面竣工结束60 d的应力变形特性为：坝体最大沉降量为123.1 cm,较施工期增大了4.3 cm,最大值发生的位置与施工期一致。最大水平位移为16.3 cm,最小水平位移为-17.4 cm,水平位移分布规律与施工期一致,最小水平位移略有减小。竣工结束60 d坝体大小主应力分布规律与施工期基本一致,心墙内部大小主应力和应力水平随着孔压的消散略有调整。

4.3 蓄水期坝体应力变形特性计算分析

图3为坝体0+280断面蓄水期的应力、变形特性图。

由图3可知,蓄水期0+280断面最大沉降量为126.7 cm,较竣工期沉降量增大了3.6 cm,沉降分布规律与施工期基本一致。最大水平位移为34.6 cm,发生在上游坝体坝壳料上部区域,方向指向下游；最小水平位移发生在下游坝坡中上部位置,最小值为-6.0 cm,方向指向上游。蓄水期大小主应力分布有所调整,心墙内部大小主应力变化比较明显。蓄水期坝体应力水平整体不高。蓄水期坝体渗流场分布合理,符合一般心墙坝的孔压分布规律。

图2 坝体0+280断面施工期不考虑流固耦合作用下的应力变形特性

图3 坝体0+280断面蓄水期应力变形特性

3.4 满蓄状态运行20 a坝体应力变形特性

图4 为坝体0+280断面蓄水运行20 a后的应力应变特性图。

由图4可知,坝体满蓄状态运行20 a后0+280断面坝体沉降量为128.9 cm,较蓄水期增大了2.2 cm,沉降分布规律基本不变。坝体应力分布和应力水平分布与蓄水期基本保持一致。

3.5 坝体耦合（蓄水）前后最大应力变形特性对比

选取坝体典型断面0+280、0+550、0+650,对不同工况计算数据进行统计分析,见表2。

坝体大部分沉降在施工期完成,施工期黏土心墙残余超孔压未完全消散,在竣工60 d后,坝体各典型断面的最大沉降量均略有增大；蓄水后,受水荷载影响以及心墙内部孔压分布的调整,坝体最大沉降量继续增大；在大坝完成蓄水,保持满库运行20 a后,坝体最大沉降量较蓄水期有所增大,并趋于稳定。竣工期坝体最大断面位置水平位移基本保持上下游侧对称分布,其他断面受地形和坝体不对称的影响水平位移极值有一定差异；蓄水期受水荷载影响,蓄水对坝体水平变形比较敏感；满库运行20 a后,坝体水平位移极值受坝体内部孔压调整影响有微小变化,极值为35.2 cm,整体规律与蓄水期保持一致。黏土心墙在竣工期和蓄水期的大小主应力最大值受孔压调整的影响逐渐减小,黏土心墙应力水平最大值也表现为逐渐减小。

图4 坝体0+280蓄水运行20 a后的应力变形特性

表2 坝体耦合（蓄水）前后最大应力、变形比较

3.6 坝体施工期实际监测最大沉降量与理论计算最大沉降量对比

坝体施工期埋设多种安全监测仪器,施工期统计至2019年2月的沉降量监测数据见表3。

由表3可知坝体施工期黏土心墙最大沉降率为1.35%,由表2可知理论计算的坝体施工期黏土心墙最大沉降率为1.25%,二者吻合较好。

4 结 论

通过对前坪水库砂砾石黏土心墙坝采用三维流固耦合计算分析,得出如下结论：

（1)通过对比监测资料和计算结果可知,坝体施工期实际监测最大沉降率与理论计算沉降率比较吻合,施工期完成了大部分沉降。

（2)黏土心墙大小主应力分布较为合理,不同时期心墙的应力水平最大值均出现在高塑性黏土与混凝土防渗墙交界区域,心墙区域绝大部分应力水平在0.6以下,表明黏土心墙有较高的强度安全储备,心墙在不同时期不会发生剪切破坏。

（3)从心墙水平位移分布规律来看,竣工期坝体最大最小水平位移分别为17.7 cm和-17.3 cm,水平位移基本沿坝轴线左右对称；竣工60 d水平位移较竣工期略有变化,量值很小；蓄水期心墙最大、最小水平位移分别为34.6 cm和-7.4 cm,发生在心墙顶部,心墙表现为向下游侧变形；满蓄运行20 a后,心墙水平位移的分布规律基本与蓄水期一致。

（4)采用流固耦合分析计算方法,通过考虑施工时间效应、蓄水时间效应以及运行期长期变形的模拟,分析了坝体在各个时期的变形特征和应力分布,体现了心墙固结特性对坝体变形和应力分布的影响。根据计算结果,工程建设管理部门在对坝顶道路、防浪墙以及围护栏施工时,应选择合适的施工时期和分缝方式,以消除坝体工后变形对建筑物的影响。

表3 坝体施工期实际监测沉降量统计