2.6 心墙堆石坝张拉裂缝及水力劈裂计算方法

2.6.1 张拉裂缝工程实用估算方法——变形倾度有限元法

通过大量土坝变形观测资料分析表明，由于可压缩土层厚度和土料性质的不一致，坝体不同区域会产生不均匀沉降，坝体横断面因不均匀沉降，会在近于铅直方向上相互错动，当错动达到一定的限度时，坝体就沿错动方向发生破坏；同时由于坝体纵向不均匀水平位移的发展，使破坏面两侧坝体相互分离，形成一定宽度的裂缝。因此，对于一般坝体，因沉降差而造成的破坏面，一般都表现为不同宽度的裂缝，为此我国的一些学者提出了根据坝体沉降观测资料来预测坝体裂缝的变形倾度法（图2.6-1）。

对基于现场监测变形的变形倾度法进行扩展，通过在有限元计算程序中嵌入变形倾度计算模块，发展了基于有限元变形计算的变形倾度有限元法。该法简洁实用，在坝体变形有限元计算分析结果的基础上，可以得到坝体在施工期或运行期任意时刻的变形倾度，据此可以进行坝体发生裂缝可能性的判断，是一种分析和判别土石坝是否会发生表面张拉裂缝的工程实用方法。采用土工离心机模型试验验证了变形倾度有限元法的适用性，试验结果表明，目前工程中采用临界倾度值γ_c取1%是基本合适的。

2.6.2 心墙压实黏土张拉断裂特性及本构模型

压实黏土，作为一种天然材料或者经过简单人工混合的半人工材料，主要包括黏土颗粒、石子、孔隙水以及孔隙气等。压实黏土的变形特性和破坏发展过程与材料中的微裂隙、孔隙的变形等都有密切的联系。可将黏土材料的裂缝扩展过程划分为微裂、亚临界扩展和失稳扩展3个阶段。微裂阶段对应拉应力水平较低的情况，此时压实黏土中微孔隙刚刚开始发生扩展和连通，土体宏观变形模量降低不明显，宏观上呈现近似线弹性。亚临界扩展阶段为压实黏土极限抗拉强度两侧附近的阶段，该阶段发生了数目较多的孔隙连通、裂纹合并和交叉，在宏观上表现为承载能力的缓慢增加或者缓慢降低。当孔隙继续连通时，会逐步形成宏观连通的裂缝，此时土体的抗拉能力会发生快速地降低直至完全丧失，这一阶段称失稳扩展阶段。

压实黏土在受拉状态下裂缝的发展可以认为是微孔隙不断扩张和沿着土粒之间黏聚薄弱部位逐步扩展的过程。在扩展过程中由于石子自身的强度较高，其“镶嵌”效果往往使裂缝的发展路径绕过石子，沿着石子和黏土颗粒结合的薄弱面向前扩展。

图2.6-2给出了拉伸试验中断裂发生后试样的照片。从图2.6-2上可见，对于压实黏土，从宏观裂缝向土样内部延伸时存在着一个损伤软化的过程区，裂缝扩展过程就是损伤过程区在土体中的延伸和扩展的过程。

在糯扎渡心墙堆石坝工程中，针对心墙土料进行了抗拉抗裂单轴拉伸和三轴拉伸试验研究，分析了压实土料的张拉断裂变形特性、破坏形式及强度准则，建议了压实黏土在受拉条件下的本构模型。

1.压实黏土单轴拉伸条件下的本构关系

对单轴拉伸试验，在土样达到极限抗拉强度发生断裂的过程中，断裂区土体和非断裂区土体表现出完全不同的应力变形特性。只有断裂区土体的应力应变关系才真正是土体抗拉应力应变的全过程曲线。根据由试验得到的全过程曲线的特性，建议了压实黏土在单轴受拉条件下应力应变全过程曲线的整理方法。对拉伸曲线的上升段建议了如下的数学表达式

式中 E——初始拉伸变形模量；

μ₀——初始切线变形模量折减系数，表示土体在承受拉应力前的初始损伤程度；

A₁和B₁——材料常数；

ε_f——峰值拉应变。

对拉伸曲线的软化段建议采用如下的负指数函数进行描述

式中 ε_u——极限拉应变；

α——材料参数，描述软化段曲线的倾斜程度。

2.压实黏土三轴拉伸条件下的本构关系

对处于受压工作状态的土石结构，当研究局部土体由受压状态变化到受拉状态过程中的应力变形特性时，要求其本构模型能够综合考虑土体的拉压变形特性。根据压实黏土的常规三轴试验和三轴拉伸试验结果，将邓肯-张EB模型发展为能够统一考虑黏性土拉压特性的邓肯-张扩展本构模型。该模型主要包括如下的部分：

（1）压缩变形段，仍采用原邓肯-张模型的计算公式

（2）直接拉伸状态下，曲线上升段采用和受压段邓肯-张模型相同的计算公式

其中，R_ft为三轴拉伸条件下的破坏比；对于三轴拉伸状态，σ₁为周围压力形成的径向应力，σ₃为轴向应力；破坏强度（σ₁-σ₃）_f根据三轴拉伸破坏联合强度准则确定；E_it为起始变形模量，具体形式如下

其中，K和n与式（2.6-3）中的K和n相等。

（3）三轴拉伸的软化段采用负指数函数形式，计算公式如下

式中 σ_rt——三轴拉伸试验确定的残余强度；

ε_f——拉伸强度对应的峰值拉应变；

α——材料参数，描述软化段曲线的倾斜程度。

（4）压缩卸载-拉伸耦合段，从压缩卸载到反向拉伸仍采用双曲线模拟应力应变关系曲线。当反向拉伸达到破坏强度时，也采用负指数函数对软化段进行模拟，压缩卸载的起始点根据加载的应力水平确定。此时，卸载初始模量为

上述邓肯-张扩展模型能够统一考虑黏性土拉压特性，其对典型试验结果的拟合情况见图2.6-3和图2.6-4。

2.6.3 土石坝张拉裂缝的有限元数值仿真算法

2.6.3.1 有限元法中裂缝的处理方法

目前，在岩石和混凝土等研究领域模拟裂缝的模型很多，有限元法中常用的处理方法有：①单元边界的单独裂缝，将裂缝处理为单元边界，一旦出现新的裂缝就增加新的节点，重新划分单元，使裂缝总是处于单元和单元之间的边界；②弥散裂缝，在弥散裂缝理论中，使用一条含平行、密集裂缝的断裂带来描述结构体中的宏观裂缝。裂缝带内的材料表现为脆性或应变软化的材料属性，而在裂缝带以外的材料保持该材料在承载状态下的一般拉伸特性；③单元内部的单独裂缝，该方法将裂缝的不连续变形特性引入单元的形函数和本构关系中。

压实黏土中的裂缝扩展过程就是断裂区在土体结构中延伸和发展的过程，故本书中引入弥散裂缝理论来模拟压实黏土的开裂行为。将裂缝弥散于实体单元，其扩展过程通过调整单元的刚度矩阵来实现，因而在模拟裂缝的扩展过程时，无需改变单元的拓扑关系和进行单元网格的重新划分。

2.6.3.2 张拉裂缝有限元-无单元耦合模拟计算方法

1.基于径向基函数的点插值无单元法

近年来，无单元法作为一种新兴的数值计算方法已经成功地应用于很多领域。无单元法只需要节点信息，无需单元信息，克服了有限元计算中网格畸变和重新生成带来的困难，可以方便地在开裂或大变形区域增加节点以提高计算精度，故其在分析裂缝扩展和局部大变形等问题方面具有优势。其中，基于径向基的点插值无单元法（RPIM），其形函数具有Delta函数性质，便于施加本质边界条件，而且形函数及其导数形式简单，计算效率高，便于与有限元法耦合，从而发挥无单元法和有限元法各自的优势。

2.点插值无单元法与有限元耦合法

为了发挥无单元法和有限元法各自的优势，清华大学提出了径向基点插值无单元法与有限元直接耦合的计算方法，并用数值算例验证了该耦合方法的有效性和适用性。

由于点插值无单元法的形函数φ(x)具有Kronecker δ属性性质、单位分解属性和重构属性，使得点插值无单元法可以与有限元法直接进行耦合。耦合界面处无需进行任何处理，耦合界面上的节点既是有限元的节点，也可作为无单元法的影响节点，见图2.6-5。在对点插值无单元法与有限元直接进行耦合时，无需保证有限元单元与无单元域的背景积分单元重合，两者可分别在各自的域内进行积分计算。为了提高耦合面处的计算精度，无单元域内节点的影响域往往延伸至有限元区域内，即有限元节点可以作为无单元法的影响节点。

点插值无单元法与以往无单元法的一个本质区别就是其形函数具备插值特性。对于点插值无单元法来说，本质边界条件可以如有限元法一样直接施加。因此，点插值无单元法与有限元进行耦合的意义更主要是在于减少计算量，提高计算效率。在计算效率方面，虽然点插值无单元法较以往的无单元方法有了较大的提高，但仍不及传统的有限元法。在较大规模的计算中，例如高土石坝三维计算中，在精度要求高的区域或者需要模拟裂缝扩展的区域采用无单元法，其他区域仍采用有限元法，可以发挥两者各自的优势，从而达到较好的计算效果。

3.基于点插值无单元法的弥散裂缝模型

在基于有限元法的弥散裂缝模型中，将裂缝弥散于整个实体单元。裂缝的扩展过程通过调整开裂单元的刚度矩阵来实现，弥散裂缝模型示意见图2.6-6（a）。弥散裂缝模型假定开裂应变均匀分布于一定的区域材料内，通过材料应力应变关系调整来反映开裂引起的材料力学性能的变化。

清华大学发展了基于无单元法的弥散裂缝模型，其将裂缝弥散于积分点的影响域内，裂缝的扩展过程通过调整开裂积分点所在影响域的刚度矩阵来实现，无单元法弥散裂缝模型的示意图见图2.6-6(b)。当积分点发生开裂后，开裂应变根据影响权重分布于影响域内，开裂过程释放的能量也根据影响权重作用在影响节点的裂缝面法向上。距离插值点较近的影响节点，其张开位移和释放的节点力都较大，随着与插值点距离的增大，张开位移和释放的节点力逐渐减小。具体的开裂计算原则如下：①若为首次开裂，则不仅本级的应力增量被释放，以前由逐级荷载累计的应力也被释放；②若已处于开裂状态，则在裂缝法向的任何拉应力增量将在每次迭代中被释放；③在某级荷载下，若在垂直于开裂方向产生了压应力，则裂缝将闭合。

基于无单元法的弥散裂缝模型模拟张拉断裂过程的示意见图2.6-7。开裂前，结构体处理成各向同性材料。当最小主应力达到抗拉强度后，结构体在垂直于最小主应力的方向产生张拉裂缝。开裂后，结构体在裂缝面的法向逐渐丧失抗拉刚度，但在其他方向仍然可以承受荷载的作用，所以可将发生裂缝后的结构体处理成各向异性材料，以模拟其在裂缝法向抗拉刚度的丧失和垂直方向的承载能力。

清华大学根据所发展的基于无单元法的弥散裂缝模型，利用上一节试验研究得出的压实黏土张拉断裂本构关系，推导了压实黏土三维脆性断裂模型和三维钝断裂带模型的无单元计算模式，并编制了相应的计算程序。该算法程序对于土石坝表面张拉裂缝问题具有较好的适用性，可用于土石坝坝体发生张拉裂缝和裂缝发生规模的计算分析。

2.6.4 土石坝水力劈裂发生过程的数值仿真算法

2.6.4.1 土石坝水力劈裂发生机制研究

1.心墙的拱效应

在通常条件下，由于土的容重远大于水的容重，竖向土压力总是大于同样深度处的水压力，所以不会发生水力劈裂。在土石坝里存在堆石体对心墙拱效应，减少了心墙中的应力，从而使得水力劈裂的发生成为可能。心墙的拱效应被认为是心墙发生水力劈裂最重要的因素之一。

在糯扎渡心墙堆石坝工程可行性研究阶段，针对直心墙和斜心墙两种坝型拟定的坝体分区方案，进行了二维非线性有限元计算分析，研究了堆石体对心墙的拱效应。为了确定一个具体的比较标准，对两种坝型分别假想了相应的“均质坝”，其心墙单元的计算应力σ₁作为判定拱效应大小的基准。拱效应G由下式计算

表2.6-1给出了总体的计算结果，可见由于堆石体对心墙拱效应的存在，使得心墙的垂直向应力显著降低，对直心墙坝尤其如此，心墙底部、中部和上部拱效应系数G的大小在坝料原参数的条件下分别达到42.7%、61%和56.4%。从这个角度讲，堆石体对心墙拱效应的存在确实是导致水力劈裂发生的重要条件。然而根据计算结果，尽管堆石体对心墙的拱效应十分明显，心墙的竖向应力仍近似为大主应力的作用方向，且压应力数值仍有一定的量级，所以心墙对发生水平向水力劈裂仍有较高的安全度。

2.渗透弱面水压楔劈效应

尽管堆石体对心墙的拱效应十分明显，但根据有限元计算经验，单凭堆石体对心墙的拱效应，不可能使得心墙在垂直方向变为小主应力方向并产生拉应力。由于在工程实践中多发生水平向的水力劈裂裂缝，联系到土石坝的水平向填筑过程可能在心墙内产生水平向渗透弱面的事实，提出了渗透弱面水压楔劈效应作用模型。认为在心墙中可能存在的渗水弱面以及在水库快速蓄水过程中所产生的弱面水压楔劈效应，应是心墙发生水力劈裂的另一个重要条件。所谓的渗透弱面是指在心墙上游面由于土料和施工不均匀、偶然掺入的堆石料、未充分压实的局部土层或者由偶然因素产生的初始细小裂缝等所造成的渗透系数相对较大的区域。由于坝体和心墙是水平向逐层填筑碾压施工的，所以出现水平向渗透弱面的可能性相对较大。

图2.6-8所示为水平向渗透弱面水压楔劈效应的产生过程以及其对发生水力劈裂诱导作用的机理。在水库蓄水时，水压力首先会沿这些水平渗透弱面快速渗入心墙，使得在心墙内产生竖直方向的水压力梯度，从而产生竖直向的渗透压力。该渗透压力作用在水平渗透弱面的上下两个边壁上，使得水平渗透弱面有张开的趋势，也即会减少心墙上游面竖直方向的应力。本书称这种由渗透弱面导致的可使心墙上游面竖直应力减少的渗透压力为渗透弱面的水压楔劈效应。

堆石体对心墙的拱效应和上述水压楔劈效应综合作用的结果可使心墙在渗透弱面处产生竖直向应力局部降低的现象。当该综合效应较大时，可使得竖向应力变成小主应力，甚至拉应力，从而导致劈裂裂缝的发生。

根据工程经验，水力劈裂一般均同水库的快速蓄水过程相联系。使用渗透弱面水压楔劈效应作用模型可以很好地解释这种现象。当库水位上升速度较慢时，一部分水压力已可渗入渗透软弱面周围的土体，从而使得水压力梯度降低，水压楔劈效应减少。反之在水库快速蓄水时，水压力来不及渗入渗透软弱面周围的土体，会产生较大的水压力梯度，从而导致较大的水压楔劈效应，增加发生水力劈裂的可能性。另外，在接近坝顶附近，由于上部的坝体压重较小，初始竖向应力较小，发生水力劈裂所需要的水压楔劈效应较小，所以是通常发生水力劈裂的危险部位。

2.6.4.2 水力劈裂发生过程的计算程序系统

清华大学将弥散裂缝理论和所建立的压实黏土脆性断裂模型引入水力劈裂问题的研究中，扩展了弥散裂缝的概念并与比奥固结理论相结合，推导和建立了用于描述水力劈裂发生和扩展过程的有限元数值仿真模型和有限元-无单元耦合数值仿真算法。

在水力劈裂的数值模拟中，一个重要的问题是模拟当裂缝发生后水压力沿劈裂裂缝的快速渗入过程，也即土体在渗透特性方面的各向异性。为此本书中将弥散裂缝的概念进行了推广，将其应用于开裂土体各向异性渗透特性的描述。通过增大单元在裂缝方向的渗透系数，模拟单元开裂后水压力沿裂缝方向的渗入过程。为此，假定平行于裂缝面方向的渗透系数分量k_t与裂缝面的法向有效应力σ_y′之间存在以下的指数方程

式中 k₀——压实黏土受压状态下的渗透系数；

α——耦合参数。

在计算中，一旦发现单元发生张拉裂缝，则需要计算裂缝张开方向，修改开裂单元的刚度矩阵并根据裂缝法向有效应力值，计算裂缝及其前缘单元的渗透矩阵，同时还需修正初始孔压场和外荷载，修改计算时间，进行迭代计算，直到无单元开裂为止。图2.6-9给出了计算程序系统的流程。

2.6.4.3 水力劈裂破坏过程的工程实例仿真模拟

利用所建立的水力劈裂发生与扩展过程的数学模型和仿真算法，对已发生水力劈裂破坏的挪威Hyttejuvet坝水力劈裂的发生与扩展过程进行了仿真模拟。

现场观测和计算结果均表明，Hyttejuvet坝施工期心墙内存在较高的超静孔隙水压力。由于该坝为变宽度狭窄心墙坝，且心墙土料和坝壳堆石料的模量相差较大，堆石料对心墙的拱效应较大。两者都使得该坝在初次蓄水时较容易发生水力劈裂破坏。

为了计算模拟该坝发生水力劈裂的过程，在心墙上游侧变宽度位置设置了3处初始渗透弱面，研究了在拱效应和渗透弱面水压楔劈效应共同作用下Hyttejuvet坝水力劈裂发生的过程。图2.6-10所示为水力劈裂发生引起的裂缝张开图。计算得到的水力劈裂发生的时间和位置同监测结果基本一致。