基于渗流与应力耦合的封闭式坝体防渗墙的应力研究
高江林
天津大学建筑工程学院
江西省水利科学研究院
为探讨土石坝加固工程中封闭式坝体防渗墙的应力状态及其适用性，利用ABAQUS有限元软件建立渗流与应力耦合数值模型，实现了同时考虑墙体与坝体相互接触、应力与渗流耦合共同作用的模拟。通过对典型土石坝工程的计算，在比较流固耦合作用对计算结果影响的基础上，对不同工程条件下的封闭式坝体防渗墙的应力进行了分析。分析结果表明，在蓄水工况下，是否考虑耦合作用对防渗墙应力的计算结果影响明显；不同墙体弹性模量、坝体土变形模量及坝基透水率等条件下的墙体应力分布有明显区别。本文分析结论对于土石坝除险加固工程中混凝土防渗墙的应用及设计计算方法具有一定的参考意义。
1 引言
在近几年国内大规模开展的病险水库除险加固工程建设中，混凝土防渗墙得到了大量应用。对于坝基覆盖层不厚的土石坝，一般采用直接嵌入基岩的封闭式坝体防渗墙进行防渗处理。在江西省近年来的土石坝加固工程中，绝大多数采用了这种形式的防渗墙。
目前，关于新建土石围堰堰体防渗墙和深厚覆盖层中坝基防渗墙的应力变形研究较多［1-4］。针对土石坝除险加固工程中增建混凝土防渗墙的应力研究相对较少。由于封闭式坝体防渗墙主要承受水平荷载，且建设防渗墙时一般在大坝运行几十年之后，土体固结对墙体应力的影响基本可以忽略不计［5］。因而，土石坝加固工程中增建的封闭式坝体防渗墙的承载形式及土体固结的影响，与一般新建土石围堰中的堰体防渗墙及深厚覆盖层中的坝基防渗墙均有明显区别。因此，针对目前土石坝加固工程中的封闭式坝体防渗墙，有必要探讨墙体在不同工程条件下的应力状况及混凝土防渗墙的适用性。
研究方法方面，早期一般采用简化方法来计算墙体应力，如采用简化地基梁进行计算［6］。简化方法应用相对简单，但由于假设条件的偏差，往往导致结果与实际情况相差较大。随着有限元分析方法在岩土工程领域的发展，国内学者采用数值计算方法对防渗墙应力变形开展了一系列研究，取得很多有益成果。但这些研究中，较少考虑渗流与应力耦合作用。而土石坝属于典型的流固耦合问题，为反映渗流对防渗墙应力的影响，往往通过水压力和湿容重（或干容重）来实现［7］，或将渗流计算出的渗透力作为荷载加到墙体上再进行应力计算［8］。间接方法需先假定或计算渗流场，分析过程繁琐，且无法从本质上反映坝体渗流与墙体应力的直接作用关系及相互影响。因此，考虑渗流与应力耦合作用是防渗墙应力变形研究方法进一步发展的必然要求。
为分析土石坝加固工程中增建的坝体防渗墙的受力性状，笔者利用ABAQUS有限元软件建立二维有限元耦合数值模型，同时考虑坝体土与防渗墙接触、渗流与应力耦合的共同作用。在对比耦合作用对墙体应力计算结果影响的基础上，重点分析不同条件下的墙体应力分布及各条件对墙体应力的影响规律，为土石坝加固工程中的防渗墙应用和设计计算提供一些有益参考。
本文发表于2013年。
2 防渗墙与坝体相互作用的数值模拟
土石坝中防渗墙应力的有限元分析，关键是处理好防渗墙与坝体土相互作用的模拟。主要包括三个方面的计算：初始应力状态、渗流与应力耦合、墙体与坝体接触（包括渗流的计算方法）。
2.1 初始应力状态
土石坝除险加固工程中增建混凝土防渗墙一般是在大坝运行几十年以后进行，此时坝体土大多已固结完成，土体固结对墙体应力的影响很小。另一方面，防渗墙施工时一般空库状态，即大坝上游无水压力作用。因此，进行防渗墙与坝体相互作用分析的初始状态应该是加固完后水库蓄水前的天然状态，坝体和墙体存在初始应力但初始位移为零。
为获得初始应力状态，须进行初始应力平衡。在ABAQUS中实现的方法为：首先将重力荷载施加于土体，并施加符合工程实际情况的边界条件，计算出在重力荷载下的应力场；根据得到的初始应力场，计算出各单元节点应力，并将其和重力荷载一起施加到原始有限元模型中计算，即可得出初始状态的应力，同时实现初始位移为零。需要注意的是，由于混凝土防渗墙施工多为水下浇筑混凝土，且施工时一般为空库状态，初始应力状态计算时不应考虑墙体与坝体的接触，否则会导致墙体初始应力偏大。
2.2 渗流与应力的耦合分析
土石坝中的渗流与应力耦合分析，需统一考虑饱和与非饱和渗流计算。本文假定非饱和与饱和土体均服从Darcy，两者的区别在于孔隙水压力和渗透系数的不同［9］。因此，非饱和渗流与饱和渗流可以用统一的方程形式，此时可将非饱和与饱和区域作为统一区域进行渗流计算。渗流与应力耦合计算需同时满足应力平衡方程和渗流连续方程。
某体积域内的平衡方程可用某一时刻的虚功原理表述：
式中：σ为柯西应力；δε为虚应变；δv为虚位移；fs为单位面积上的面力；f为单位体积上流体重力之外的体积力；s为土体饱和度；n为土体孔隙率；ρw为流体的密度；g为重力加速度。
计算域用有限元网格离散后，单元网格中流体运动满足连续性，即单位时间内流入土体的流体体积等于流体体积的增加速率，其连续方程可表示为：
式中：n为渗流面S的法线方向；vw为流体的平均流速；ρow为流体参照密度。
在上述应力平衡方程和渗流连续方程的基础上，采用Galerkin有限元格式，将节点位移和孔隙水压力作为节点自由度进行空间离散，形成应力平衡方程和渗流连续方程的矩阵形式。在同时满足位移边界和渗流边界的条件下，直接对耦合控制方程进行求解。
2.3 墙体-坝体的接触单元
2.3.1 接触模拟
防渗墙与坝体土的接触采用ABAQUS提供的接触单元进行模拟，接触行为包含接触面的法向作用和切向作用。接触面的法向模拟采用硬接触［10］，即两物体只有在压紧状态时才能传递法向压力p，若两物体之间有间隙时则不能传递法向应力。
接触面的切向作用采用库仑摩擦模型，用摩擦系数μ来模拟在两个表面之间的摩擦行为。当接触面处于闭合状态（即有法向接触压力）时，接触面可以传递切向应力（摩擦力）。设定极限剪应力：τcrit=μp。若摩擦力小于某一极限值τcrit时，认为接触面处于黏结状态；若摩擦力大于τcrit后，接触面开始出现相对滑移变形，称为滑移状态。
2.3.2 渗流模拟
为实现整个计算域的渗流与应力耦合计算，必须采用能够考虑渗流与应力耦合作用的接触单元。本文采用ABAQUS提供的能够考虑孔压-位移耦合的接触单元，除满足上述传力机理外，还满足以下渗流条件：
图1 接触面单元渗流模式示意图
（1）孔压在接触对应的实体上完全连续，用公式表示为pA-pB=0，其中pA、pB为接触面两侧节点上的孔隙水压力。
（2）接触面的切线方向无流体流动。
（3）接触面压力为有效压力，不包含孔隙水压力的作用。
接触面单元渗流模式如图1所示。
3 计算模型
本文以一般地基上的均质土坝作为分析对象。坝高30.0m，坝基覆盖层厚5.0m，防渗墙自坝顶嵌入至基岩0.5m，墙深35.5m，墙厚0.6m，强风化和弱风化基岩厚度均取为10.0m，正常蓄水位25.0m，下游水位为坝脚位置。具体计算条件如图2所示。
有限元计算网格如图3所示。
土体采用Mohr-Coulomb弹塑性模型，采用变形模量进行计算。根据本文计算结果，考虑到混凝土压应力水平较低及抗拉状态时的应力应变曲线接近线性关系［5，11］，且整个基岩应力水平较低，本文计算中防渗墙和基岩材料采用弹性模型。对于封闭式防渗墙而言，坝基岩层一般为渗透系数较小或进行过灌浆处理的基岩，可作为相对不透水层看待。部分有限元计算参数见表1所示。
图2 计算模型图（单位：m）
图3 有限元网格图
表1 主要计算参数
4 耦合作用对防渗墙应力结果的影响
为分析耦合作用对墙体应力计算结果的影响，分别对是否考虑耦合时墙体的最大、最小主应力进行比较。不考虑耦合时，仅考虑重力荷载和迎水面水压力荷载作用，无坝体渗流作用。由于ABAQUS中应力以拉为正，因此，最大主应力反映的是墙体最大拉应力，而最小主应力则反映墙体的最大压应力。以刚性混凝土防渗墙为例进行分析。
4.1 墙体应力
是否考虑耦合时墙体的最大、最小主应力分布如图4所示。
由图4可知，考虑耦合作用时防渗墙最大拉应力和压应力分布为4.51MPa和3.56MPa，均出现基岩面位置附近。最大拉应力大于压应力，相对混凝土的极限承载能力而言，墙体破坏由拉应力控制。拉应力和压应力呈对称分布，说明墙体承受典型的弯矩作用，墙侧土摩阻力引起的墙体压应力贡献相对较小。
图4 防渗墙的大、小主应力
图5 刚性墙墙身弯矩分布
不考虑耦合作用时墙体的最大主应力和最小主应力均为负值，说明此时防渗墙整体处于受压状态，此时墙侧摩阻力引起的压力作用更加突出。最大压应力为1.18MPa，远小于墙体材料的允许值，说明不考虑渗流与应力耦合会导致墙体应力计算结果偏小，尤其是控制墙体破坏的拉应力。
4.2 墙身弯矩分布
为探讨正常蓄水工况下防渗墙的受力机理，进一步分析墙身弯矩分布。利用墙体各单元节点上的应力结果，分别计算出防渗墙迎水面和背水面的竖向应力，进一步可推算出墙体弯矩。墙体弯矩沿深度分布如图5所示。弯矩方向以顺时针为正。
根据图5可以看出，考虑耦合作用时墙体弯矩在中部（21.0m深左右）和底部出现两处极值，分别为-84.6kN·m和245.8kN·m，之间弯矩发生反向，说明弯矩主要有中下部水平荷载引起。最大弯矩值出现在墙体底部。弯矩的大小与分布规律与应力相对应。可以推断，对于高坝或更高水位工况下墙体会出现更加不利的应力状态。
不考虑渗流与应力耦合时，最大弯矩值为26.7kN·m，出现在墙体底部，分布规律类似，但弯矩值远小于考虑耦合作用的计算结果。
4.3 坝体应力场和渗流场的比较
为说明耦合作用对墙体应力影响，进一步分析耦合作用下的坝体渗流场和应力场。
4.3.1 坝体孔压分布
渗流与应力耦合作用下坝体孔压分布如图6所示。
图6 坝体孔压分布
图6所示的灰色区域为非饱和土，分界线为坝体浸润线位置。由于防渗墙的截渗作用，墙体前后水头有明显突变，水头差约12.5m。
4.3.2 坝体应力场
考虑耦合时和不考虑耦合时的坝体竖向有效应力分布如图7和图8所示。为方便比较，仅列出坝体部分的土体有效应力分布。
图7 坝体竖向有效应力分布（考虑耦合）
图8 坝体竖向有效应力分布（不考虑耦合）
由图7可知，考虑耦合作用时的最大竖向有效应力为0.44MPa，相对集中出现在防渗墙底部附近，且防渗墙前后土体的有效应力相差明显，且整体上墙后土体的有效应力大于墙前土体，说明防渗墙形成向下游挤压的趋势。根据墙体前后土体有效应力差和墙体弯矩分布，可以进一步推断，引起墙体向下游侧挤压的主要作用力来自于上游土体中的孔隙水压力。
由图8可知，不考虑耦合作用时的最大坝体竖向有效应力为0.54MPa，由于无孔隙水压力作用，最大有效应力值大于考虑耦合时的计算结果。此时防渗墙前后坝体的有效应力基本相当。由此看出，不考虑渗流与应力耦合时，无法反映墙体前后水头差对墙体应力的影响。
5 不同条件下的墙体应力分析
结合实际工程情况，为探讨不同条件下的混凝土防渗墙应力情况及其在土石坝中的适用性，本文重点分析不同混凝土弹模、坝体填筑材料、坝基透水率等条件下的墙体应力及变化规律。根据前面分析可知，封闭式坝体防渗墙主要承受水平荷载作用，墙体破坏由拉应力控制，下文仅对不同条件下的最大主应力（拉应力）进行分析。
5.1 不同混凝土弹模条件下的墙体应力
在近年来江西省已建成的土石坝坝体防渗墙中，墙体混凝土的抗压强度从几兆帕的塑性混凝土到二十几兆帕刚性混凝土都有，相应弹性模量差别也很大。在工程质量控制上，一般仅把混凝土抗压强度和渗透系数作为墙体材料的质量控制指标，很少对墙体材料的弹性模量进行控制。
为分析混凝土弹性模量对墙体应力的影响，分别对刚性和塑性混凝土防渗墙的应力进行计算。刚性混凝土弹性模量Ec分别取为2.5×104MPa和2.0×104MPa，塑性混凝土弹性模量取为1500MPa和500MPa。不同弹模对应的最大主应力如图9所示。
根据图9可知，墙体应力随着墙体材料弹性模量的增大而增大，且刚性混凝土防渗墙的拉应力的远大于塑性墙，并超出材料本身的极限值。而塑性混凝土防渗墙拉应力均在允许范围内，当Ec=500MPa时，墙体基本无拉应力存在。由此看出，对于承受水平荷载为主的封闭式坝体防渗墙，墙体弹模对应力影响显著。墙体弹性模量越小（与坝体土弹模越接近），越有利于墙体受力，即对墙体受力而言，混凝土抗压强度不是越大越好。因此，在实际工程中应增加对墙体材料弹性模量的上限控制。
图9 不同混凝土弹模时的墙体最大主应力
5.2 不同土体模量条件下的墙体应力
由上述分析可知，墙体和坝体土弹性模量相差越大，对于墙体受力越不利。在实际工程中，由于坝体填筑材料种类的不同，坝体土本身的变形模量同样存在一定的差异。为探讨混凝土防渗墙在不同坝体填筑材料中的应力情况及其适用性，进一步分析不同坝体土变形模量对刚性和塑性墙体应力的影响。坝体土变形模量Es分别为50MPa、100MPa、300MPa、600MPa。
不同土体变形模量时刚性混凝土防渗墙的最大主应力分布如图10所示。
图10 不同土体模量时刚性墙的最大主应力
由图10可以看出，坝体土的变形模量对刚性混凝土墙体应力影响明显，土体模量越小，墙体拉应力越大。当坝体填筑材料为软黏土时（Es<100MPa），墙体应力急剧增加，最大拉应力达到6.0MPa，大大超出混凝土材料的极限值。当坝体填筑材料为模量较大的硬质黏土或砂性土时（Es>300MPa），对墙体受力更为有利。因此，进行刚性混凝土防渗墙设计时，尤其是对于坝体填筑材料为变形模量较小的软黏土时，应充分考虑墙体材料弹模过大对墙体应力的不利影响。
进一步分析不同坝体土变形模量时塑性墙体的应力情况，如图11所示。
由图11可知，增建塑性混凝土防渗墙时，墙体最大拉应力同样随坝体土变形模量的减小而增大，但塑性混凝土墙的最大拉应力远远小于刚性混凝土墙。一般情况下，塑性墙体的最大拉应力均在材料允许值范围内。从墙体受力角度考虑，对于不同的坝体填筑材料，塑性墙较刚性墙具有更好的适用性。
图11 不同坝体土模量时塑性墙的最大主应力
5.3 不同坝基透水率时的墙体应力
由于地质条件和所采取工程措施的不同，实际土石坝工程中的坝基透水率存在一定差异。风化程度较低或进行了灌浆处理的坝基透水率往往较小，而对于一般强风化的基岩透水率则较大。一般情况下，坝基渗透系数在1.0×10-4～1.0×10-6cm/s之间。为分析不同坝基透水率时混凝土防渗墙的应力情况，坝基渗透系数Kr分别取为1.0×10-4cm/s、1.0×10-5cm/s、1.0×10-6cm/s，此时模型中的强风化和弱风化基岩采用相同的渗透系数。不同坝基渗透系数对应的刚性混凝土防渗墙的最大主应力分布如图12所示。
图12 不同坝基透水率时刚性墙的最大主应力
图13 不同坝基透水率时塑性墙的最大主应力
根据图13可以看出，坝基渗透系数越小，墙体最拉应力越大。这主要是由于坝基透水率越小，墙体前后形成的水头差越大（三个渗透系数对应的水头差分别为7.5m、9.5m、12.5m），墙体所承受的水平荷载也越大，相应引起的拉应力也更大。当坝基渗透系数Kr=1.0×10-6cm/s时，最大拉应力为4.51MPa，已超出混凝土极限值。当Kr=1.0×10-4cm/s时，最大拉应力为2.84MPa，不同地基透水率的最大拉应力值相差近2倍，说明坝基透水性能对刚性墙体应力的影响不能忽略。
不同坝基渗透系数对应的塑性混凝土防渗墙的最大主应力分布如图13所示。
由图13可知，塑性混凝土防渗墙的最大拉应力同样随坝基渗透系数减小而增大，当Kr=1.0×10-6cm/s时，最大拉应力为0.15MPa；当Kr=1.0×10-4cm/s时，最大主应力为-0.17MPa（整体为受压状态）。不同坝基透水率时的塑性混凝土防渗墙最大拉应力相差很小。说明不同坝基透水率对塑性混凝土防渗墙应力的影响较小。
可以看出，墙体应力与坝基透水性本身是一对矛盾，从基础防渗角度，坝基透水率越小越好，而从墙体受力而言，透水率越大，对墙体受力越有利。而建设防渗墙的根本目的就是防渗，针对坝基透水率小的情况，应首先考虑建设塑性混凝土防渗墙。
6 结语
（1）本文基于渗流与应力耦合分析，实现了防渗墙与坝体之间相互接触、渗流与应力耦合的共同作用的模拟。分析结果表明，所建立的耦合数值模型能较好地反映防渗墙与坝体的相互作用，可为土石坝中混凝土防渗墙的应力变形计算提供方法参考。
（2）是否考虑渗流与应力的耦合作用，对于防渗墙的应力结果影响明显。正常蓄水工况下，封闭式坝体防渗墙主要承受墙体前后水头差引起的水平荷载作用。仅考虑上游坝坡水压力作用（不考虑渗流）进行墙体应力计算是偏于不安全的。
（3）防渗墙墙体材料的弹模对于墙体应力影响明显，墙体材料弹模越大，墙体应力越大。水平荷载作用下，塑性混凝土应力远小于刚性混凝土墙，此时塑性混凝土防渗墙具有明显的承载优势。为保证封闭式坝体防渗墙的承载能力和工程安全，建议在防渗墙设计时明确墙体材料弹性模量的上限值。
（4）坝体填筑材料的变形模量对刚性混凝土墙体应力影响明显，坝体材料模量越小，墙体应力越大。当坝体材料为模量很小的软黏土时，对于刚性混凝土防渗墙受力极为不利。不同坝体填筑材料变形模量对塑性混凝土墙应力的影响相对较小，对墙体受力而言，塑性墙具有更广泛的适应性。因此，在土体模量相对较小的土石坝中增建混凝土防渗墙时，一定要充分考虑墙体材料弹模过大对墙体应力的不利影响。
（5）坝基透水率对刚性混凝土防渗墙的应力影响明显，坝基透水率越小，水平荷载越大，刚性混凝土防渗墙的墙体应力越大。相对而言，坝基透水率对塑性混凝土防渗墙的应力影响较小，不同渗透系数时的墙体应力均在材料允许范围内，且相差不大。因此，针对坝基透水率小的情况，应优先考虑建设塑性混凝土防渗墙。
参考文献
［1］沈新慧.防渗墙及周围土体的应力探讨［J］.水利学报.1995（11）：39-45.
［2］卢廷浩，汪荣大.瀑布沟土石坝防渗墙应力变形分析［J］.河海大学学报.1998，26（2）：41-44.
［3］王刚，张建民，濮家骝.坝基混凝土防渗墙应力位移影响因素分析［J］.土木工程学报.2006，39（4）：73-77.
［4］熊欢，王清友，高希章，等.沙湾水电站一期围堰塑性混凝土防渗墙应力变形分析［J］.水力发电学报，2010，29（2）：197-203.
［5］王清友，孙万功，熊欢.塑性混凝土防渗墙［M］.北京：中国水利水电出版社，2008.
［6］郑秀培.土石坝地基混凝土防渗墙设计与计算［M］.北京：水利电力出版社，1979.
［7］沈珠江，左元明.三峡工程二期混凝土防渗墙围堰初步设计方案的应力应变计算［J］.长江科学院院报.1987（3）：31-37.
［8］段伟，骆原，陈珂，等.黄壁庄水库副坝渗流及防渗墙应力分析［J］.人民黄河.2011，33（5）：124-127.
［9］费康，张建伟.ABAQUS在岩土工程中的应用［M］.北京：中国水利水电出版社，2010.
［10］ABAQUS Documentation［CP/DK］.ABAQUS Theory Manual（Version 6.10），2010.
［11］高丹盈，王四巍，宋帅奇，等.塑性混凝土单向受压应力—应变关系的试验研究［J］.水利学报.2009，40（1）：82-87.
本文发表于2013年。