高土石围堰复合土工膜防渗三维有限元分析

徐 晗,陈 云,饶锡保,潘家军

(长江科学院 水利部岩土力学与工程重点实验室，武汉 430010)

高土石围堰复合土工膜防渗三维有限元分析

徐 晗,陈 云,饶锡保,潘家军

(长江科学院 水利部岩土力学与工程重点实验室，武汉 430010)

如何在三维有限元计算中反映复合土工膜的力学特性是亟待解决的技术难题。总结了现有的土工膜数值模拟方法，在三维有限元计算中采用薄膜单元来模拟土工膜，建立了某高65 m的高土石围堰土工膜防渗三维有限元模型，土工膜与上、下游坝料接触部位均设置接触。为模拟施工填筑与蓄水全过程，分别计算了土工膜边界柔性约束与刚性约束2种工况。结果表明:薄膜单元较好地反映了土工膜柔性抗拉材料的力学特性，能适应三维有限元数值模拟中常用的分步填筑施工模拟法，且后处理中能直接展示土工膜的位移、应力云图，具有方便高效的特点；土工膜边界柔性约束位移分布均与刚性约束相同，只是最大主应力及拉应力分布区域略小于刚性约束，计算结果符合土石坝工程的一般规律。

复合土工膜；薄膜单元；三维有限元；高土石围堰；变形特性

1 研究背景

复合土工膜(简称土工膜)具有防渗效果好、铺设方便等优点，可铺设在坝体上游作斜墙或在土石坝中间作心墙，已广泛应用于水利工程建设中[1-2]。防渗体系是土石坝工程的生命线，如果土工膜拉裂或拉断，将给整个工程的防渗带来安全隐患。现有规范[3]规定，高水头(>50 m)挡水建筑物采用土工膜防渗应经过计算论证，由于土工膜厚度薄，不具备抗压和抗弯特性，与坝料的力学特性差异较大，需要在土工膜与坝体之间设置非连续接触面，因而在计算中如何考虑土工膜的实际作用存在较大的技术难度。

虽然学者们发展了解析法[4]、有限差分法[5]等多种方法来计算土工膜的应力变形，但更多的是理论探讨，较少运用于实际土石坝工程。土石坝应力变形计算多采用有限元法，在二维计算中可以将土工膜简化为柔索单元[6]或不设置土工膜单元[7]，但土工膜三维有限元数值模拟研究甚少，所采用的力学模型仍需要进一步研究论证，因而在三维有限元计算中考虑土工膜的作用是亟待解决的技术难题。

2 三维有限元土工膜模拟方法

在三维有限元计算中，主要采用如下几种方法来考虑土工膜的作用。

(1) 土工膜内嵌入实体单元：该方法将土工膜内嵌入实体单元，并采用切向弹簧模拟土工膜与坝体之间的切向摩擦作用;在法向方向假定土工膜随所嵌入的坝体单元一起运动。由于实际工程中蓄水时水压力作用下土工膜与上、下游坝体之间有可能脱开，而采用该方法无法模拟此脱开性状，同时水压力也只能施加在所嵌入的单元面上，因而无法反映土工膜的真实受力性状。

(2) 采用实体单元建立土工膜：该方法需要对土工膜建立实体单元，并划分网格，土工膜上、下游侧均与上、下游坝体之间设置接触面，来模拟坝体与土工膜之间的接触。虽然该方法可以在土工膜上加水压力，并模拟土工膜与坝体之间的张拉脱开、剪切滑移等接触力学特性，但是实际土工膜厚度只有1～2 mm，在三维模型中建模难度极大，且实体单元无法反映土工膜单元不具备抗压和抗弯特性的特点。

(3) 土工膜依附在单元表面：该方法假定土工膜依附在单元表面上，不单独建立土工膜单元，而根据加载前后土工膜所在位置单元表面的节点位移变化，来换算土工膜应变，此计算方法非常粗略，仅能反映土工膜大致的受力性状。

上述的几种三维有限元计算中模拟土工膜的方法均存在明显的缺点，本文采用ABAQUS提供的薄膜单元来模拟土工膜[8]。薄膜单元是一种面单元，只能传递面内的力，不承受弯矩，即它没有弯曲刚度，在ABAQUS中常用的薄膜单元为M3D4R (三维4节点减缩积分薄膜单元)，单元的正、负法向面均可以设置非连续接触面。虽然有学者对薄膜单元在土石坝中的应用进行了初步的探讨[9]，但并没有模拟土工膜随堰体施工填筑而逐步增高的过程，因而需要对薄膜单元是否适用于土石坝数值模拟中常见的填筑过程进行进一步的论证。

上述薄膜单元属于柔性抗拉材料，可以直接设置土工膜的厚度，但不具备抗压和抗弯特性；同时在土工膜与上、下游坝料之间均可设置接触，接触面的力学模型又可分为法向模型和切向模型，可以模拟土工膜与坝料之间的法向张拉脱开、切向剪切滑移等接触力学特性；在此方法中，水压力也直接加在土工膜上，因而在三维有限元计算中较好地反映了土工膜自身的基本力学特性。

表1 堰体材料的邓肯-张E-μ模型参数Table 1 Parameters of Duncan-Zhang E-μ model of cofferdam material

3 有限元模型及模拟的填筑过程

某围堰堰顶宽15 m，最大堰高65 m，堰体由石渣和砂砾石组成，采用混凝土防渗墙上接土工膜防渗，建立的三维有限元模型如图1所示，防渗墙+土工膜防渗纵断面如图2所示，其中土工膜最大高度为45 m。

图2 防渗墙+土工膜纵断面Fig.2 Longitudinal section of impervious wall and geotechnical membrane

图3 土工膜接触面 设置示意Fig.3 Schematic diagram of the contact interface between geotechnical membrane and dam

三维有限元计算中采用薄膜单元来模拟土工膜，对于土工膜与坝料接触部位也设置接触单元，反映土工膜与坝料之间的接触特性。土工膜与坝体之间的相互接触如图3所示。根据室内试验成果，土工膜与堰体填料间的界面摩擦系数为0.29，防渗墙与覆盖层之间的摩擦系数参照三峡二期围堰实测值取0.19。复合土工膜采用线弹性模型，其厚度取2.5 mm，弹性模量为100 MPa。

模拟施工的全过程，荷载分级按堰体填筑次序进行，模拟的施工程序为：围堰地基初应力分析→分级填筑截流堤→水下分级填筑上游堆石→水下分级填筑下游堆石→施工防渗墙→水上分级填筑碾压堆石→铺设过渡料、坝料及土工膜→围堰分级蓄水，下游降水到围堰底。

堰体材料采用的邓肯-张E-μ模型参数见表1，其中：c为内聚力；φ为摩擦角；Rf为破坏比；K为初始变形模量的基数；n反映变形模量随围压增长的关系；G为初始泊松比的基数；F反映初始泊松比随围压增长而降低的速率；D反映初始泊松比随轴向应变增长的关系。

4 土工膜应力变形分析

4.1 土工膜变形分析

图4为土工膜沉降等值线，可知土工膜变形在堰体最大横断面处位移值最大。表2为蓄水期堰体与土工膜变形极值，可知土工膜与坝体竖向位移量值基本一致，沉降最大值发生在堰体中部，坝轴向位移量值也一致，只是顺河向位移极值有所差异，这是由于堰体水平位移极值发生在下游堰体内，如图5所示。通过分析图5可知，在土工膜与上、下游堰体交接处，竖向位移与水平位移等值线均有明显的突变，表明该模型较好地反映了土工膜与上、下游堰体之间的非连续接触变形特性。

图4 土工膜变形等值线Fig.4 Deformation contours of geomembrane表2 蓄水期堰体与土工膜变形极值Table 2 Extreme values of the deformation of cofferdam and geomembrane impoundment stage

cm

图5 最大横断面位移等值线Fig.5 Displacement contours at the maximum cross section

4.2 土工膜应力分析

土工膜为柔性构件，只能承受拉应力，不能承受压应力，因而计算结果显示了土工膜所受到的拉应力数值及区域。

表3为刚性约束条件下土工膜拉应力、应变极值，图6为土工膜最大主应力、最大主应变等值线。由表3可知，土工膜的拉应力可以直观地显现出来，刚性约束下的最大主应力0.67 MPa，最大主应变0.71%。土工膜刚性约束后的拉应力主要集中在土工膜的上部及周边约束地区，这是由于该地区变形较大或者受约束所致；而由于土工膜与防渗墙交接地带处于坝体中下部，该部位埋深较大，下游堰体的抗力较强，土工膜变形较小，因而没有拉应力产生。

表3 土工膜应力、应变极值Table 3 Extreme values of the stress and strain of geomembrane

图6 土工膜最大主应力、最大主应变等值线Fig.6 Contours of maximum principal stress and maximum principal strain of geomembrane

综上所述，土工膜不具备抗压和抗弯特性，因而主要呈现拉应力，没有压应力产生，计算结果较好地体现了土工膜柔性抗拉材料的力学特性。

5 边界约束对土工膜应力变形的影响

工程实践表明，在土工膜与周边结构刚性约束时，在接头处容易产生应力集中，致使土工膜拉裂[10]。为了验证土工膜边界约束施加方法的正确性，在薄膜单元与结构连接处，分别施加了刚性约束与柔性约束2种边界条件，刚性约束为土工膜边界各位移分量全限制，柔性约束则释放部分位移分量，土工膜的应力变形结果见表4。

表4 边界约束对土工膜应力变形的影响结果Table 4 Influence of boundary constraint on the stress and deformation of geomembrane

由表4可知柔性约束条件下，蓄水期土工膜水平位移最大值为11.6 cm(向下游)，竖向位移最大值为43.2 cm，坝轴向位移最大值为5.9 cm，与刚性约束边界位移值相比，蓄水期土工膜各向位移均与刚性约束下相同；只是最大主应力为0.59 MPa，最大主应变为0.57%，略小于刚性约束条件。

图7为蓄水期不同边界约束条件土工膜拉应力区域云图，可知在土工膜与两侧山体搭接部位容易产生应力集中，柔性约束条件下土工膜的拉应力区域及数值明显小于刚性约束条件，因而在土石坝工程中，土工膜边界应尽量采用柔性约束。

6 结 论

(1) 采用ABAQUS提供的薄膜单元来模拟土工膜，可在土工膜与上、下游坝料接触部位设置接触，薄膜单元能较好地反映土工膜柔性抗拉材料的力学特性，在土石坝三维计算中能适应常用的分步填筑施工模拟法，且后处理中能直接展示土工膜的位移、拉应力与应变等值线，具有方便高效的特点。

(2) 当土工膜边界为刚性约束时，在蓄水期水压力的作用下，在土工膜与两侧山体搭接部位易产生应力集中；土工膜边界柔性约束下各向位移分布均与刚性约束相同，只是最大主应力及拉应力分布区域小于刚性约束，计算结果符合土石坝的一般规律。

[1] 顾淦臣.复合土工膜或土工膜堤坝的实例述评(续)[J].水利水电技术, 2003, 34(1): 56-59.

[2] 包承纲. 三峡二期围堰创新及其对我国高围堰建设的意义[J]. 长江科学院院报, 2014, 31(9): 33-42.

[3] SL/T 225—98, 水利水电工程土工合成材料应用技术规范[S]. 北京: 中国水利水电出版社, 1998.

[4] 高 登, 朱 斌, 陈云敏. 圆形构筑物周边土工膜沉降时的褶皱产生机理及分析模型[J]. 岩土工程学报, 2008, 30(8): 1126-1132.

[5] 吴兆和, 沈长松, 瞿忠烈, 等. 复合土工膜堆石坝膜受力变形特性分析[J]. 水电能源科学, 2010, 28(3): 91-94.

[6] 花加凤, 束一鸣, 张贵科, 等. 土石坝坝面防渗膜中的夹具效应[J]. 水利水电科技进展，2007, 27(2): 66-68.

[7] 姜海波. 复合土工膜心墙与斜墙高土石坝应力应变研究[J]. 长江科学院院报, 2014, 31(1): 53-57.

[8] 庄 茁. ABAQUS有限元软件6.4版入门指南[M]. 北京: 清华大学出版社, 2004.

[9] 尚 层.土石坝复合土工膜防渗斜墙应力变形分析[D].乌鲁木齐: 新疆农业大学, 2012.

[10]李 波, 程永辉, 程展林. 围堰防渗墙与复合土工膜联接型式离心模型试验研究[J]. 岩土工程学报, 2012, 34(11): 2081-2086.

(编辑：黄 玲)

Three-dimensional Finite Element Analysis on Composite Geomembraneof High Earth Rockfill Cofferdam

XU Han，CHEN Yun，RAO Xi-bao，PAN Jia-jun

(Key Laboratory of Geotechnical Mechanics and Engineering of Ministry of Water Resources,Yangtze River Scientific Research Institute, Wuhan 430010, China)

To reflect the mechanical properties of composite geomembrane in 3-D finite element calculation is a technical problem to be solved. On the basis of summarizing current numerical methods of geomembrane simulation, we propose to simulate the geomembrane using membrane element in 3-D finite element calculation. By building a finite element model of high earth-rock cofferdam of 65 m height with geomembranes contacted with dam in the upstream and downstream, we simulated the flexible and rigid constraint conditions of geomembrane boundary respectively in filling construction and impounding stages. Results suggest that the membrane element could well reflect the flexible tensile properties of geomembrane. The proposed method could also adapt to the stepped filling in 3-D finite element simulation. Moreover, it could conveniently and effectively display the displacements and stresses of geomembranes in post-processing. The displacement distribution of geomembrane under flexible boundary restraint was the same with that under rigid boundary constraint; while the principal tensile stress and its distributed region under flexible boundary constraint were slightly smaller than those under rigid constraint. The calculation results are consistent with the general rules of earth-rockfill dam.

composite geomembrane; membrane element; 3D finite element; high earth rockfill cofferdam；deformation characteristics

2016-06-23;

2016-08-23

国家自然科学基金项目(51309029)；长江科学院创新团队项目(CKSF2015051/YT)；长江科学院院所基金项目(CKSF2015036/YT)

徐 晗(1978-)，男，湖北公安人，教授级高级工程师，博士，主要从事岩土工程数值分析研究，(电话)027-82829743(电子信箱)mechanics007@aliyun.com。

10.11988/ckyyb.20161056

TV641.1

A

1001-5485(2017)02-0110-04

2017,34(2):110-113