深厚覆盖层上土石坝防渗墙损伤开裂精细化分析及防渗功能评价

邹德高，屈永倩，孔宪京，刘京茂，冉从勇，张 丹

(1.大连理工大学 海岸与近海工程国家重点实验室，辽宁 大连 116024；2.大连理工大学 水利工程学院，辽宁 大连 116024；3.中国电建集团成都勘测设计研究院有限公司，四川 成都 610072)

1 研究背景

深厚覆盖层在我国特别是西南地区的河流中广泛分布，如大渡河、金沙江等，其覆盖层深度可达数百米。我国80%以上的水能资源集中在西部地区，深厚覆盖层上建坝“难以避让”[1-2]。土石坝以其独有的优势，是深厚覆盖层地基上的首选坝型。覆盖层上土石坝工程的地基防渗主要采用竖直防渗墙[3]，如狮子坪、黄金坪、泸定、下坂地、冶勒、旁多等覆盖层深度超过100 m的土石坝工程，均采用混凝土防渗墙的坝基防渗形式，其中旁多、冶勒的覆盖层深度超过400 m，防渗墙深度分别超过100 m和200 m[4-5]。目前已有防渗墙破坏影响大坝工程安全的案例，如Manic3水库的混凝土防渗墙与基岩接触部位发生了挤压破坏[6]；牛头山水库蓄水后，混凝土防渗墙在墙体顶部发生了开裂破坏[6-7]。因此，开展混凝土防渗墙的损伤破坏分析及对防渗性能的影响研究，对科学把握其真实工作状态，进而合理评估深厚覆盖层上土石坝的安全性态至关重要。

覆盖层中混凝土防渗墙作为坝基挡水结构，属于防渗结构而非承载结构。然而目前防渗墙大多采用线弹性模型简化模拟，并采用基于应力的承载能力理念进行评价。线弹性模型计算的混凝土防渗墙应力有时远超材料的抗拉或抗压强度，且超应力区范围可能较大，这些结果虽能给出混凝土防渗墙的相对危险位置，但是不能准确评价墙体的安全状态。实际上，当应力达到强度后，材料就会进入非线性状态(或损伤状态)，应力将发生转移或重新分布，而线弹性模型无法描述这种状态。同时，线弹性模型获得的墙体应力状态很难直接与墙体的防渗性能建立联系。因此，有必要采用更为完备的本构模型模拟墙体混凝土材料，研究覆盖层中混凝土防渗墙的真实状态，进而开展渗流分析，从防渗性能的角度对防渗墙进行评价。

混凝土防渗墙与覆盖层的刚度和尺寸相差百倍甚至千倍，两者相互作用十分明显，导致防渗墙的精细分析尤其是损伤破坏分析十分困难。本课题组于2014年开展了防渗墙的损伤破坏分析[8]，并研究了三维河谷效应、防渗墙墙间夹泥对防渗墙损伤破坏特性的影响规律[5，9]。多家科研机构和研究者们也相继开始采用非线性和弹塑性模型研究防渗墙损伤破坏规律[10-12]。目前的覆盖层上土石坝防渗墙的损伤开裂分析结果相对较少，且局限于覆盖层深度100 m内的封闭式防渗墙，对超深覆盖层上土石坝悬挂式防渗墙的损伤破坏特性还鲜有研究。本构模型方面多采用塑性损伤模型[13]描述混凝土的破损特性，但该模型采用标量损伤因子描述混凝土损伤后的刚度和强度退化特性，导致一个方向发生损伤后，另外两个方向的模量和强度均同程度降低，且拉、压损伤互相干扰，难以精准描述复杂应力状态下的混凝土损伤破坏特性及诱发的材料各向异性问题。

针对上述问题，本文建立了有限元-比例边界元耦合的超深覆盖层上土石坝防渗墙受力性态的跨尺度精细分析模型，采用混凝土塑性损伤模型和内聚力模型分别描述受压破坏和受拉开裂特性，提出了深厚覆盖层上土石坝防渗墙破坏演化分析方法，避免了传统采用标量损伤因子难以描述混凝土损伤后诱发的材料各向异性问题。开展了深厚覆盖层上土石坝防渗墙的损伤开裂分析，阐明了超深覆盖层上悬挂式混凝土防渗墙的应力和变形机制，揭示了混凝土防渗墙的损伤开裂模式，定位了防渗墙的薄弱区域。在此基础上，开展了防渗墙破损后的渗流分析，评价了防渗墙损伤开裂对防渗性能的影响，提出了相应的措施并量化了其防渗效果。

2 理论方法

2.1 塑性损伤模型Lee和Fenves[13]在Barcelona塑性损伤模型[14]的基础上，分别采用拉、压两个损伤变量描述混凝土的损伤断裂特性，并将屈服函数推至有效应力空间，提出了Lee-Fenves塑性损伤模型。该模型在模拟混凝土的力学特性方面具有较大优势且便于数值实现，目前已在混凝土坝工程中得到广泛应用[15-16]。文献[17-18]将该模型引入土石坝领域并进行了改进，发展了可考虑不同纤维含量的钢纤维混凝土塑性损伤模型，并应用于面板坝混凝土面板的动力损伤演化分析。本构模型的介绍和验证参见文献[13，17-18]。但是该模型采用标量损伤因子描述混凝土损伤后的刚度和强度退化特性，导致一个方向发生损伤后，另外两个方向的模量和强度均同程度降低，且发生拉损伤后会影响受压性能，难以精准描述复杂应力状态下的混凝土损伤破坏特性。

2.2 内聚力模型该模型认为在裂纹尖端存在一个微小的内聚力区(Fracture process zone)，以内聚力抵抗界面间的相对分离。在开裂过程中，界面上的应力为开裂位移的函数，从而避免了线弹性断裂力学中的裂纹尖端应力奇异性问题[19]。采用内聚力模型模拟开裂时，在裂缝可能发生和扩展的部位布置界面单元(cohesive interface element，CIE)，界面单元与周围的实体单元相连，见图1。随荷载增大，界面单元的应力状态达到破坏准则，刚度和承载能力降低，当刚度降低到0时，界面单元失效，新的裂缝面出现。鉴于双线型模型适用于描述脆性材料的断裂，且拥有较高的计算效率[20]，因此本文采用双线型模型描述内聚力单元的牵引力分离关系。该模型已应用于混凝土坝[21-22]，但在土石坝领域应用较少。文献[23-24]将内聚力模型应用于混凝土面板坝的面板动力开裂模拟，并量化了抗震措施效果，本构模型的介绍和验证参见文献[23-24]。但是该模型仅能描述混凝土的拉伸、剪切和拉(压)剪破坏，无法描述受压破坏。

图1 内聚力模型和界面单元

2.3 有限元-比例边界元跨尺度分析方法比例边界有限元法(SBFEM)[25-26]由Song和Wolf于1997年提出，是一种半解析方法，并继承了边界元的优点。该方法支持任意边数多边形和任意面数的多面体单元，能灵活地处理更复杂几何形状单元的求解问题，在跨尺度网格的疏密网格过渡处理中具有无法比拟的优势。Chen等[27-29]基于该方法原有边界高斯积分点，首先通过SBFEM理论构造用于非线性分析的单元插值函数和位移应变矩阵，采用材料常刚度矩阵方法计算相关系数矩阵，规避传统耗时繁琐的计算求解过程，改善分析效率。其次，通过增加单元内部高斯点，采用内部积分方案求解位移应变协调矩阵、单元刚度矩阵和应力积分等，发展了可用于弹塑性分析的多边形和多面体单元，解决了SBFEM受限于弹性分析的问题。并在同一程序框架内集成了FEM-SBFEM耦合算法，发展了FEM-SBFEM的跨尺度分析方法，并应用于土石坝的三维静动力分析[23-24，27-30]。

3 深厚覆盖层上土石坝防渗墙的精细化分析

3.1 有限元模型深厚覆盖层上的沥青混凝土心墙坝分析模型如图2所示。坝高148 m。坝体上游坡度1∶2.5，下游坡度1∶2.0，坝顶宽10 m。覆盖层厚度超过500 m，分为4层，从上至下，依次为①层、②层、③层、④层，底部基岩为第⑤层。防渗墙深140 m，延伸至第③层；宽1.4 m，顶部有一扩大端，并与心墙底部相连。

图2 深厚覆盖层上沥青混凝土心墙坝的三维分析模型

采用作者发展的有限元-比例边界元-增量迭代法-预设薄层单元的跨尺度非线性分析方法[30]，精细描述防渗墙的受力状态。防渗墙底部采用有限元-比例边界元耦合跨尺度分析方法，实现防渗墙局部网格的多级细分，将该区域网格加密至2.0 m；通过比例边界元的多面体单元连接交界处网格，实现不同网格尺寸的跨越；在心墙-防渗墙接头区和防渗墙底部设置薄层单元模拟局部土体呈带状的大应变剪切变形特性，以精细描述防渗墙的受力性态[30]。最终建立三维的有限元-比例边界元耦合的防渗墙精细分析模型，并采用增量迭代法进行求解。

计算采用大连理工大学自主研发的大型岩土工程高性能分析软件系统 GEODYNA，该软件已集成了有限元-比例边界元-无网格-离散元耦合的多数值分析方法，有中点增量法、增量迭代法、等效线性方法等隐式求解方法和显式求解方法，已广泛应用于土石坝的静、动力分析[23-24，27-33]。

土体采用非线性弹性邓肯张E-μ模型，坝体和覆盖层材料的计算参数见表1。在防渗墙与覆盖层土体、心墙与过渡料、心墙与防渗墙之间设置Goodman接触单元，采用双曲线模型描述其接触特性，参数见表2。

表1 坝料参数表

表2 防渗墙与覆盖层间双曲线接触面模型材料参数

3.2 悬挂式防渗墙的受力变形模式为便于讨论防渗墙的受力变形特性，首先采用线弹性模型描述防渗墙。防渗墙的受力示意见图3，竣工期防渗墙的变形模式见图4，主应力及其方向见图5和图6。竣工期，由于上、下游坝体基本对称，防渗墙顺河向位移很小。受坝体传来的上覆土压力和覆盖层摩擦力的作用，悬挂式防渗墙发生向下的竖向位移，由于河谷存在一定坡度，防渗墙向下位移的同时，会出现指向河谷的坝轴向压缩变形。(1)防渗墙两岸嵌入基岩，起到了支撑作用，因此防渗墙两岸底部压应力最大，压应力方向近似垂直于两岸岸坡。(2)两岸底部受基岩的支撑作用，防渗墙发生类似于“外伸梁”变形模式，底部支撑对应的顶部拉应力最大，其方向为坝轴向。支撑部位向河谷延伸，过“反弯点”后，防渗墙顶部出现高压应力，底部出现高拉应力，方向均为坝轴向。(3)河谷中部防渗墙的主应力方向基本为竖向和坝轴向，均为压应力。(4)重力作用下，覆盖层主要发生竖向压缩变形，而防渗墙的模量远大于覆盖层土体，其本身的压缩变形很小，防渗墙将会刺入两端土体，因此防渗墙中上部区域受向下的摩擦力，而防渗墙中下部受向上的摩擦力，见图7。因此，对于河谷中部防渗墙，在摩擦力方向转换的位置，防渗墙压应力最大。

图3 防渗墙受力示意图

图4 竣工期防渗墙变形模式(变形放大30倍)

图5 防渗墙主应力(受压为正)

图6 防渗墙主应力方向

图7 防渗墙受到的覆盖层摩擦力方向

3.3 防渗墙损伤破坏分析由上节的弹性分析结果可知，防渗墙的压应力达到了150 MPa，拉应力达到了44 MPa，远超混凝土强度，这是不符合实际的。因此本节采用混凝土塑性损伤模型和内聚力模型描述防渗墙的损伤破坏特性。其中，塑性损伤模型描述受压破坏，内聚力模型描述受拉破坏。由弹性分析结果可知，防渗墙的拉应力方向主要为坝轴向。同时防渗墙施工时存在竖向分布的槽段接缝，见图8，该接缝界面的强度和模量均要低于防渗墙混凝土。因此，综合考虑拉应力方向和防渗墙施工工艺，在防渗墙的槽段接缝处预设竖向的内聚力界面单元并采用内聚力模型描述其受拉破坏特性，防渗墙墙体采用实体单元和塑性损伤模型描述其受压破坏特性，见图9。塑性损伤模型参数见表3，槽段界面处的内聚力模型参数相比墙体考虑一定的弱化，由于缺少试验结果，本文参考了新旧混凝土结合面的强度折减系数[34]。由于槽段界面无法做凿毛、清洗等处理，且槽段界面可能存在墙间夹泥，这将进一步降低结合面强度，因此折减系数取较小值，本文在计算中取为0.25。需要说明的是，防渗墙槽段界面为弧形，为简化建模与计算，本文偏保守地将槽段界面简化为平面。

表3 混凝土塑性损伤模型参数

图8 防渗墙槽段接缝

图9 防渗墙损伤开裂模拟的模型说明示意图

竣工期的防渗墙损伤开裂分析结果见图10，防渗墙应力见图11。考虑损伤开裂之后，防渗墙的超压应力和超拉应力均被释放，而在防渗墙两岸顶部和底部拉应力较大的区域发生开裂，其中，防渗墙左岸中上部轴向150 m范围内出现竖向裂缝，最大裂缝宽度为2.6 cm，单条裂缝最大长度为70 m，裂缝总长度为841 m；防渗墙右岸中上部轴向190 m范围内出现竖向裂缝，最大裂缝宽度为2.4 cm，单条裂缝最大长度为75 m，裂缝总长度为930 m；防渗墙左岸底部轴向320 m范围内出现竖向裂缝，最大裂缝宽度为0.8 cm，单条裂缝最大长度为30 m，裂缝总长度为501 m；防渗墙右岸底部轴向350 m范围内出现竖向裂缝，最大裂缝宽度为1.6 cm，单条裂缝最大长度为30 m，裂缝总长度为375 m。防渗墙两岸发生压损伤，损伤面积约700 m2，约占防渗墙面积的0.3%，其中两岸底部最大损伤因子超过0.95。

图10 防渗墙损伤开裂分析结果

图11 考虑损伤开裂后的防渗墙应力(受压为正)

弹性分析和损伤开裂分析的防渗墙变形对比见图12，考虑混凝土损伤开裂后，两岸底部的混凝土发生严重的压损伤，刚度降低，防渗墙竖向位移增大约0.7 m，高压应力释放，但是在稍微远离防渗墙与基岩接触区域，防渗墙依然存在高压应力，超过了24 MPa。两岸“支撑”对应的防渗墙顶部不再出现反弯变形，而河谷中部防渗墙的竖向位移差别不大。防渗墙上部依然存在较大的压应力，达到了45 MPa，这是因为塑性损伤模型采用多轴强度进行损伤判别，该区域防渗墙在坝轴向、竖向和顺河向三个方向均处于受压状态，因此最大压应力超过单轴抗压强度时，混凝土并不一定会发生压损伤。

图12 防渗墙变形模式对比(变形放大15倍)

3.4 渗流分析及防渗性能评价防渗墙完好状态和考虑损伤开裂后的防渗墙压力水头分布见图13和图14，渗漏量见表4。防渗墙完好时，总渗漏量为74 842 m3/d，其中防渗墙部分的渗漏量为326 m3/d，仅占总渗漏量的0.4%。而防渗墙损伤开裂后，防渗墙下游面压力水头在损伤区域明显增大，其中两岸顶部损伤区的压力水头30 m增加到70 m，两岸底部损伤区的压力水头由190 m增加至210 m。大坝总的渗漏量增加至114 389 m3/d，增大52.8%，而防渗墙区域的渗漏量由326 m3/d增加至50 430 m3/d，增大153.6倍，占总渗漏量的44.1%。

表4 渗漏量

图13 防渗墙完好状态下的总水头分布(单位：m)

图14 考虑损伤开裂后的防渗墙下游面水头分布(单位：m)

由表4可知，相比于防渗墙底部损伤开裂区，防渗墙左右岸顶部的裂缝宽度较大，渗漏量较大，分别达到了17 518 m3/d和21 271 m3/d，占防渗墙区域渗漏量的76.9%。由于顶部的破坏区域更小、破坏更为集中且埋深浅，防渗处理相对容易，因此应重点对防渗墙两岸顶部区域采取措施。本工程采用悬挂式防渗墙，防渗墙底高程以下不设防渗结构，覆盖层的渗漏量已达56 328 m3/d，而防渗墙底部的损伤开裂程度相对较小，其渗漏量约为覆盖层渗漏量的20%，且损伤范围大、埋深大、处理困难，因此防渗墙两岸底部区域采取措施的性价比并不高。

综合上述分析结果，可以在防渗墙两岸的上游侧设置辅助防渗墙(图15)，厚度和深度与主防渗墙一致，间距为2 m，左岸和右岸辅助防渗墙的坝轴向长度分别为160 m和200 m，均比防渗墙顶部开裂区长度大10 m。主防渗墙和辅助防渗墙间采用灌浆处理，渗透系数取为3.0×10-5cm/s[36]。

考虑处理措施的防渗墙下游面压力水头分布见图16，上部的压力水头与防渗墙未损伤工况下的结果基本一致，而防渗墙底部未处理区的压力水头分布基本不变。考虑处理措施的渗漏分析结果见表4，与未处理工况相比，两岸顶部损伤区的渗漏量减小97.8%，大坝总渗漏量降低31.7%；与设计目标(防渗墙完好状态)相比，大坝总渗漏量仅增加了4.4%，防渗处理措施起到了很好的效果。

图16 损伤开裂后考虑处理措施的防渗墙下游面水头分布(单位：m)

4 结论

本文提出了比例边界元-有限元耦合跨尺度离散、塑性损伤模型和内聚力模型分离描述压损伤和受拉开裂、破损后防渗功能目标评价的精细化分析方法，开展了深厚覆盖层上土石坝的应力变形分析和渗流分析，阐明了超深覆盖层上悬挂式防渗墙的应力和变形机制、揭示了混凝土防渗墙的损伤开裂模式，定位了防渗墙薄弱区域，评价了防渗墙损伤开裂对防渗性能的影响，提出了相应的处理措施并量化了其性能效果。主要结论如下：

(1)两岸基岩的支撑作用使防渗墙两岸底部出现近似垂直于岸坡的高压应力；类“外伸梁”的面内弯曲变形使得防渗墙两岸底部支撑对应的顶部区域出现坝轴向的高拉应力，随着向河谷延伸，经过“反弯点”后，防渗墙顶部出现坝轴向的高压应力，底部出现坝轴向的高拉应力；河谷中部防渗墙出现竖向和坝轴向的主压应力。

(2)靠近防渗墙两岸的顶部和底部区域的坝轴向高拉应力导致槽段间出现竖向裂缝；防渗墙两岸底部因近似垂直于岸坡的高压应力发生压损伤；防渗墙两岸损伤后刚度降低，沉降增加，而河谷中部防渗墙的竖向位移基本不变。发展的损伤开裂分析方法可揭示防渗墙的破坏模式、精准定位防渗墙的薄弱位置、量化破坏程度，并与防渗性能建立直接关系，为防渗墙的性能评价提供基础。

(3)在防渗墙两岸的上游侧局部设置辅助防渗墙，并在主防渗墙和辅助防渗墙间采用灌浆处理，可大幅降低防渗墙损伤开裂后的渗漏量，起到了很好的防渗效果。

(4)本文提出的防渗墙精细化分析和评价方法，可从防渗性能角度评价防渗墙破坏的影响，并量化措施的防渗效果，实现了防渗墙从传统承载能力评价到功能性态评价的跨越，可为深厚覆盖层上土石坝防渗墙的安全评价和设计优化提供理论依据和技术支持。