高面膜堆石坝防渗结构受力变形数值分析方法——高面膜堆石坝关键技术(四)

姜晓桢，束一鸣

(1．南京水利科学研究院岩土工程研究所，江苏南京 210024;2．水文水资源与水利工程科学国家重点实验室，江苏南京 210098;3．河海大学水利水电学院，江苏南京 210098)

堆石坝是一种较古老和成熟的坝型，多采用混凝土面板和黏土心墙等防渗结构。随着土工合成材料在岩土、水利等工程中广泛应用，土工膜作为一种新型的防渗材料也被应用到堆石坝工程中，发挥其优良的防渗效果。与传统的防渗结构相比，土工膜防渗结构具有适应变形能力强、结构简单、施工方便等优点，是一种值得进一步推广的新型防渗结构［1-2］。目前国际上采用土工膜防渗结构最高的堆石坝是西班牙97m高的Poza de Los Ramos堆石坝，该坝从竣工至今一直运行良好，目前已加高到134 m［3］。我国自20世纪80年代以来，已有数十座新老堆石坝采用土工膜防渗或防渗加固，其中新建最高的是四川仁宗海堆石坝，高56 m［4］。根据国际大坝委员会统计，目前土工膜已经在167座大型堆石坝工程中得到了应用，取得了较好的工程效果［5］。

由于土工膜是一种新型的堆石坝防渗结构，所以采用数值计算方法分析土工膜的受力变形规律时，在模型的建立和分析中都需对传统方法进行改造和完善。国内土工膜防渗结构数值计算起源于20世纪80年代的二维索单元方法［6］，但目前为数不多的与土工膜防渗结构受力变形分析相关的文献中仍存在很多不足之处。束一鸣等［6-7］利用有限元法对深厚覆盖层上土工膜防渗堆石坝的筑坝特性进行了分析，但并未介绍土工膜防渗结构在堆石坝有限元计算中的具体方法和原理。吴兆和等［7-10］利用有限差分软件FLAC3D中的土工格栅单元对堆石坝中土工膜防渗结构的受力变形进行了计算模拟，虽然从原理上基本能反映土工膜防渗结构在堆石坝工作状态下垫层对土工膜的摩擦剪切作用，但无法对锚固处的夹具效应进行模拟。

1 土工膜防渗结构的形式与受力特点

土工膜是一种片状柔性防渗材料，具有良好的适应变形的能力，其防渗结构的基本形式如图1所示，从上到下依次为保护层、土工膜层和垫层。保护层一般为预制的或直接在膜上现浇的混凝土板，其主要作用是保护土工膜免于被水库漂浮物划伤、划破以及直接的紫外线照射引起的老化破坏。垫层材料一般由散粒体材料或其聚合物(无砂混凝土、聚合物拌和混凝土等)组成，其作用一是对土工膜起到承托作用将水压力均匀地传递到坝体和地基;二是利用散粒体材料良好的透水性，起到排水降压的作用。

图1 土工膜防渗结构形式

土工膜防渗结构的施工一般在坝体填筑完成后进行，所以其产生的变形主要是由于蓄水后随着坝体的变形而变形:①膜面由原来的平面凹陷成“锅面”，这种变形会在坝轴线和坝坡两个方向上同时给土工膜带来拉伸变形;②在双向拉伸变形的同时，土工膜还会与两侧的垫层和保护层产生沿着接触面的剪切错动，所以土工膜与两侧材料接触面的力学特性，也影响着土工膜的受力变形特性;③土工膜周边是锚固在河床和岸坡的刚性基岩上的，由于坝体和岸坡基岩之间材料性质差异较大，往往会发生如图2所示的差异变形，由于此时接近锚固处的土工膜在水压力作用下，被紧紧压在坝面上，无法均匀和有效地伸展开来，所以在接近锚固处极小的长度范围内，产生一个较大的拉伸变形，从而引起土工膜的拉伸破坏，这种受力特点也被形象地称为夹具效应。

图2 夹具效应示意图

2 土工膜防渗结构的有限元模型

2．1 土工膜接触面连接单元

为了模拟土工膜与保护层、垫层之间在接触界面的剪切作用，必须在有限元建模过程中在土工膜单元与坝体单元之间建立接触面连接单元。传统的接触面连接单元中接触面两侧单元节点的坐标必须相同，如此则造成了土工膜单元网格的密度与坝体单元网格的密度相同。事实上土工膜在坝体、岸坡基岩以及坝基等处的应力和变形比较复杂，所以往往要求这些部位的土工膜网格尺寸不宜过大，这样势必造成堆石坝坝体单元和土工膜单元在网格尺寸上的不匹配，这也是目前土工膜防渗结构有限元计算的一个难点。

传统的接触面连接单元接触面两侧单元节点的坐标必须相同，即以面对面的方式进行连接，如图3所示。如果将传统的面对面连接方式改为点对面的连接方式，那么就可以实现土工膜网格与坝体网格的独立划分，如图4所示。三角形土工膜单元的每个节点分别隶属于不同的(或者相同的)的坝体网格单元内(面上)，所以1个土工膜接触面连接单元中包含了1个土工膜单元的一个节点以及该节点所在的坝体三维实体单元的所有节点。

图3 传统接触面连接单元

图4 土工膜接触面连接单元

在节点力Fe的作用下，土工膜接触面连接单元内产生平行于接触面的两个相互垂直的相对切向位移ws1、ws2和1个垂直接触面的法向相对位移wr，即

相应地，在土工膜与坝体接触面上产生两个切向应力 τs1、τs2和，1 个法向应力 σr，即

在线弹性假定下，应力σ与相对位移w成正比，即

式中ks、kr分别为切向、法向的单位面积劲度系数，N/m3。法向劲度系数kr可取为与坝体弹性模量相近的一个数，而切向的劲度系数则可根据不同的接触面本构关系，由试验确定。

根据有限元法的基本原理［11］，坝体单元内任意一点的位移都可以由坝体单元节点的位移插值而得，所以当土工膜接触面连接单元中的土工膜节点位于坝体单位内(或其表面上)时，坝体单元在该点上的位移可表示为

式中:NFE为坝体单元的插值形函数;ue为坝体单元的节点位移。

如此土工膜接触面连接单元的相对位移可由土工膜接触面连接单元内各个节点的位移表示，并通过坐标转化，将土工膜单元局部坐标系下相对位移转化为三维整体坐标下的相对位移:

式中:I为三阶单元矩阵(矩阵对角线元素为1，其余为零);M为土工膜接触面单元局部坐标系(r，s1，s2)与整体坐标系(x，y，z)之间的转化矩阵;ae为土工膜接触面连接单元各个节点的位移列阵。

综合式(1)～(6)，根据虚位移原理可推得

式中:ke为土工膜连接单元的单元劲度矩阵;D为土工膜连接单元中土工膜节点所代表的面积，本文中其值为围绕该土工膜节点的所有土工膜单元按面积大小进行权值平均后的面积。

根据上述的分析，可以发现计算土工膜接触面连接单元的单元劲度矩阵时与传统的接触面单元相类似，唯一不同的是，计算之前需确认连接单元中土工膜节点在坝体单元内的位置，并以此确定插值形函数的具体数值。当坝体三维实体单元为四面体单元时，可根据四面体单元的体积坐标确定4个插值形函数的具体数值:

式中:V为四面体单元体积;x、y、z分别为连接单元中土工膜单元节点的坐标值;xj、yj、zj、xm、ym、zm、xp、yp、zp分别为坝体四面体单元节点坐标。

当坝体三维实体单元为六面体等参单元或五面体等参单元时，由于进行了等参变换，所以形函数数值并不是三维整体坐标的显式表达式，无法直接求得，但还是可以通过迭代的数值方法得出形函数数值［12］。

2．2 土工膜夹具效应的有限元模型

根据土工膜夹具效应的受力特点，如果从有限元数值计算模型的角度来看，在建立网格模型时需做出一定的改变，如图5所示。传统的方法是将土工膜节点与基岩节点以共节点的形式来模拟锚固约束，所以网格节点是连续的，如此则无法体现坝体、岸坡以及坝基等处的差异变形。如果要对差异变形进行模拟，那就必须在土工膜单元节点和基岩单元节点之间再建立新的单元才有可能计算出差异变形的大小。本文将原来锚固处的节点分离开来，在土工膜节点与基岩节点之间建立新的三维空间弹簧单元，并对锚固处土工膜的受力变形情况进行分析。

图5 夹具效应有限元网格模型示意图

三维空间弹簧单元包含1个基岩节点和1个土工膜节点，总共有6个自由度。在节点力Fe的作用下，两个节点间会产生1个三维空间的相对位移:

相应地，在三维空间弹簧单元内产生1个三维空间内的弹簧内力:

弹簧内力σ'与相对位移w'可根据弹簧的劲度有如下的关系:

式中kg为弹簧的劲度系数，N/m，可由土工膜的弹性模量与膜厚度的乘积来计算。

三维空间弹簧单元的相对位移可由三维空间弹簧单元两个节点位移值求得

综合式(9)～(13)，根据虚位移原理可推得

式中:ke为三维空间弹簧单元的单元劲度矩阵;Ls为三维空间弹簧单元中土工膜锚固处节点所代表的锚固段长度，本文中取为围绕该锚固节点的所有土工膜单元锚固边按边长权值平均后的长度。

3 土工膜防渗结构的有限元计算

3．1 模型网格

以某已建堆石坝为例，对其采用的面膜防渗结构进行有限元建模和计算，并对本文所提出的数值计算方法的合理性和适用性进行分析。

该堆石坝坝高178m，坝顶宽12 m，坝长1104 m，上、下游坝坡度均为1∶1.4，两岸岸坡基本对称，坡度为25°。大坝整体有限元网格和最大断面坝体分区网格如图6所示，面膜网格如图7所示。

图6 坝体有限元网格

图7 面膜网格

三维坐标系建立如下:顺水流指向下游为x轴正向;沿坝轴线指向左岸为y轴正向;竖直向上为z轴正向。坝基基岩厚度为150 m，水平方向上分别从坝踵和坝趾向上下游延伸150 m。整个模型实体单元全部划分为三维等参六面体单元，共计3024个单元， 3865个节点;土工膜三维平面膜单元全部划分为三角形单元，共计2104个单元，1143个节点。

计算共分为施工和蓄水两种工况进行。整个计算过程分为3个步骤:①建立坝基及两岸基岩模型，并进行初始地应力计算，然后位移清零，并保留单元应力信息;②分9层建立坝体模型，并模拟坝体施工填筑过程;③在土工膜节点上分9步逐级施加170m的水头，模拟水库蓄水过程，所以总共进行了18个荷载增量步的计算。

3．2 材料本构模型和计算参数

坝体堆石料为非线性材料，其应力-应变关系呈非线性，坝体各分区堆石料采用邓肯E-B模型，其参数根据相关工程资料选取［13-14］，如表1所示。土工膜材料的应力-应变关系也具有非线性特性，采用土工膜非线性弹性模型，其切线弹性模量Et和切线泊松比 νt按下式［15］计算:

表1 坝体各分区堆石料邓肯E-B模型参数

式中:σ1、σ2分别为土工膜单元的大、小主应力;C1、C2、…、C6分别为土工膜非线性弹性模型的参数，各参数取值如表2所示。

表2 土工膜非线性弹性模型参数

根据土工膜接触面直剪试验的结果发现其剪切位移与剪切应力关系基本符合克拉夫-邓肯的双曲线模型［16］，所以本文土工膜接触面连接单元的切向劲度系数ks的取值按下式计算:

式中:K1、n、Rf为通过试验确定的参数;ρw为水的密度;Pa为大气压力;δ为接触面材料摩擦角;c为接触面黏聚力;τs为接触面上的切向应力。各参数取值如表3所示。

表3 土工膜接触面双曲线模型参数

3．3 计算结果分析

土工膜防渗结构在有限元建模过程中运用了3种不同的单元来模拟土工膜防渗结构的受力变形特性:①在土工膜四周锚固处的三维空间弹簧单元;②土工膜与坝面的接触面连接单元;③由平面三角形单元组成的面膜单元。这3种单元各自反映了土工膜在不同部位的受力特点，三维弹簧单元的计算结果能反映出锚固处差异变形对土工膜张拉作用的大小，土工膜接触面连接单元的计算结果能反映出土工膜和坝面之间沿着接触面相对错动的大小以及错动所带来的摩擦剪切应力的大小，而坝面上的平面三角形土工膜单元则能综合反映出面膜内部的应变和应力分布。

图8为面膜的位移分布情况。从图8可知，水库蓄水后土工膜防渗结构在水压力的作用下，随着坝面一起发生变形，由原来的平面凹陷成了“锅面”，凹陷最大处位于坝面中下部1/3坝高处。坝轴线方向的位移基本呈对称分布，两侧靠近岸坡处的土工膜随着凹陷的“锅面”向中心变形。

图8 面膜位移分布情况 (单位:cm)

为了防止这种夹具效应的发生，可将土工膜沿着差异变形的反方向铺设一段距离，即为将来可能发生的差异变形留下了足够大的预留量。相关试验结果表明，这两种铺膜方式下，无论是坝基处还是岸坡处的夹具效应的影响都被大大削弱。预留量的大小可根据数值计算或结构模型试验的结果，同时结合工程实际进行灵活设定。

4 结语

土工膜防渗结构在水压力作用下，由原来的平面凹陷成“锅面”，同时坝面上土工膜向坝面凹陷中心发生位移，在这种变形规律下，坝面上土工膜基本处于受拉状态，越靠近坝基和岸坡锚固处部位的土工膜拉应力越大，而坝面中部土工膜几乎不受拉应力作用。在这种应力分布规律下，提高锚固处的土工膜网格的密度，降低坝面中部的网格密度，就显得很有必要，而土工膜接触面连接单元恰好能够解决坝面网格和土工膜网格在尺寸大小上不一致的问题。

对于夹具效应，本文改变了传统基岩处土工膜单元与基岩单元共节点的连接方式，在基岩节点和土工膜节点之间增加了一组三维弹簧单元，通过弹簧单元的三向变形来反映差异变形，从而对锚固处土工膜的受力变形进行了有效的模拟，计算结果反映的规律符合一般工程实际。

［1］顾淦臣．复合土工膜或土工膜堤坝实例述评［J］．水利水电技术，2002，33(12):26-32．(GU Ganchen．Review ofdams and dikes with composite geomembrane impervious structure［J］．Water Resources and Hydropower Engineering，2002，33(12):26-32．(in Chinese))

［2］顾淦臣．复合土工膜或土工膜堤坝实例述评(续)［J］．水利水电技术，2003，34(1):55-61．(GU Ganchen．Review of dams and dikes with composite geomembrane impervious structure(Ⅱ)［J］．Water Resources and Hydropower Engineering，2003，34(1):55-61．(in Chinese))

［3］《土工合成材料工程应用手册》编写委员会．土工合成材料工程应用手册［M］．2版．北京:中国建筑工业出版社，2000．

［4］忘海洋．仁宗海水库电站大坝上游坡面复合土工膜防渗工程施工技术［J］．吉林水利，2010，2(12):66-70．(WANG Haiyang．The construction technology of the compound genmembrane against seepage applied in the upstream slope protection of Renzonghai reservoir power station dam［J］．Jilin Water Resources，2010，2(12):66-70．(in Chinese))

［5］Geomembrane sealing systems for dams-Design principles and return of experience(Bulletin 135)［R］．Paris:The International Commission on Large Dams，2010．

［6］束一鸣，顾淦臣，土工薄膜中央防渗土石坝有限元分析［J］．河海大学学报，1988，16(增刊 1):79-91．(SHU Yiming，GU Ganchen．The FEM analysis of the rockfill dam using genmembrane as central barrier［J］．Journal of Hohaiuniversity，1988，16 (Sup1):79-91．(in Chinese))

［7］岑威钧，沈长松，童建文，等． 深厚覆盖层上复合土工膜防渗堆石坝筑坝特性研究［J］． 岩土力学，2009，30 ( 1) : 175-180． ( CEN Weijun，SHEN Changsong，TONG Jianwen， et al． Study of construction behavior feasibility of composite geomembrane rockfill dam on thick alluvium deposit［J］． Rock and Soil Mechanics， 2009， 30( 1) : 175-180． ( in Chinese) )

［8］吴兆和，沈长松，瞿忠烈，等． 复合土工膜堆石坝膜受力变形特性分析［J］． 水电能源科学， 2010， 28( 3) : 91-94．( WU Zhaohe，SHEN Changsong，QU Zhonglie，et al．Analysis of mechanical behaviors of composite geomembrane used in rockfill dam ［J］． Water Resources and Power， 2010， 28( 1) : 91-94． ( in Chinese) )

［9］花加凤，束一鸣，张贵科，等．土石坝坝面防渗膜中的夹具效应［J］．水利水电科技进展，2007，27(2):66-68．(HUA Jiafeng，SHU Yiming，ZHANG Guike，et al．Clamping effect in geomembranes used for seepage control of upstream surface of earth-rock dams［J］．Advances in Science and Technology of Water Resources，2007，27(2):66-68．(in Chinese))

［10］束一鸣，李永红． 较高土石坝膜防渗结构设计方法探讨［J］． 河海大学学报: 自然科学版，2006，34( 1) : 60-64．( SHU Yiming，LI Yonghong． Design of impervious structure with geomembrane for high earth-rock dams［J］．Journal of Hohai University: Natural Sciences，2006，34( 1) : 60-64． ( in Chinese) )

［11］陈国荣．有限单元法原理及应用［M］．北京:科学出版社，2008．

［12］徐燕萍，项阳，刘泉声，等．等参元逆变换插值法的研究及其应用［J］．岩土力学，2001，22(2):226-228．(XU Yanping，XIANG Yang，LIU Quanshen，et al．Research on numerical inverse isoparametric mapping interpolation and its application ［J］．Rock and Soil Mechanics，2001，22(2):226-228．(in Chinese))

［13］DUNCAN J M，CHANG C Y．Nonlinear analysis of stress and strain in soils［J］．Journal of Soil Mechanics and Foundation Division，1970，96(5):1629-1653．

［14］单宏伟．高面板堆石坝模型参数反演及应力变形分析［D］．北京:清华大学，2008．

［15］姜晓桢．土石坝坝面防渗土工膜运行性态及其分析方法［D］．南京:河海大学，2013．

［16］ClOUGH G W，DUNCAN J M．Finite element analysis of retaining wall behaviour［J］．Journal of Soil mechanics and Foundation Engineering，1971，12(10)，1657-1672．