特高混凝土面板堆石坝的应力变形比较研究

陈 平，朱家奇，王庆祥

(1.中国电建集团中南勘测设计研究院有限公司，湖南 长沙 410014；2.河海大学水利水电学院，江苏 南京 210024)

近年来，我国水利水电建设发展迅速，在土石坝填筑设计理论和筑坝技术等方面均取得了长足的进步，面板堆石坝正成为当今坝工建设中的主要坝型之一。由于水电工程的复杂性和建设运行经验的欠缺，DL/T 5016—2011《混凝土面板堆石坝设计规范》规定对高混凝土面板堆石坝应进行变形和应力有限元计算复核；同时，考虑到堆石料本构模型尚未成熟，对于特高混凝土面板堆石坝要求采用不同本构模型进行对比验证，以提高计算成果的精度和适用性。因此，选择合适的本构模型对特高混凝土面板堆石坝的有限元分析十分重要。

学者们一直致力于堆石料本构关系的研究并提出了相应模型[1]。在众多堆石体计算的本构模型当中，以邓肯-张模型[2]为代表的弹性非线性模型，由于其参数简单易懂、物理意义明确、适用性强、并且能够能充分反映出土体的非线性特性而得到广泛应用，但该模型不能反映堆石料高紧密度条件下的剪胀性[3]。沈珠江[4]采用邓肯-张模型应力应变双曲线将硬化规律通过显式表示，提高了双屈服面模型的实用性，应用较广；但该模型高估了堆石料的剪胀性，使得堆石体的计算变形一般偏小。以往传统的塑性理论常常需要提前假设出土体的塑性势面和关联流动法则，并不考虑其应力主轴的旋转，而广义塑性模型[5- 6]不需要采用塑性势面的显式表达式及屈服面，可以通过定义微分形式的加载方向和塑性流动方向，并简单方便地定义出不同加载情况下的塑性模量，以完成各种加载条件下土体应力变形的预测任务。陈生水等[7]提出了基于广义塑性理论的本构模型，模型通过改变摩擦角反映堆石料颗粒破碎导致峰值应力随平均应力减小。朱晟等[8]采用更具普适性的剪胀方程，提出了适应于复杂应力路径的静动力耦合统一广义塑性模型。

拉哇水电站混凝土面板堆石坝目前为世界上在建的最高的混凝土面板堆石坝，坝顶总长398.0 m，坝顶宽15.0 m，坝顶高程2 709.0 m，最大坝高239.0 m。河床深槽部位趾板设置基座混凝土结构，拟最大厚度为5 m，基座顶部高程2 475.0 m，上游面板底部坝坡坡比1∶1.45，下游坝坡在2 590.0 m高程以上布置“之”字形上坝道路，综合坡比为1∶1.628。本文以拉哇水电站的特高混凝土面板堆石坝为例，通过采用邓肯E-B模型和统一广义塑性模型两种计算方法，计算其坝体变形及应力分布情况，以此进一步分析堆石体和面板的变形协调特性，为坝体分区结构设计提供相应依据。

1 堆石体计算本构模型

1.1 邓肯E-B模型[2]

堆石体在自重荷载及蓄水期水压力的作用下往往会表现出明显的非线性特征，因此可以采用充分反映土体非线性的邓肯E-B模型，切线弹性模量Et和体积变形模量Bt的计算公式为

(1)

(2)

(3)

(4)

式中，pa为大气压力；K、n、Rf、Kb、m、φ0、Δφ为模型参数。

1.2 河海统一广义塑性模型[8]

统一广义塑性模型的弹塑性矩阵为

(5)

式中，De为弹性矩阵；Dp为塑性矩阵；ngL/U为加卸载过程中对应的塑性流动方向；HL/U为加卸载过程中对应的塑性模量；n为加载时的方向。

根据等向压缩特性，堆石体的弹性模量为

(6)

式中，υ为泊松比，一般取0.3。

堆石体的塑性流动方向，在p-q平面内也就是它对应的剪胀方程为

(7)

式中，β，α，Mg为粗粒土参数。

加载时的塑性流动方向为

(8)

为了模拟堆石体的“卸载体缩”现象，将卸载时的塑性流动方向定义为

(9)

采用非关联流动法则，加载方向为

(10)

加载时的塑性模量为

(11)

由于颗粒破碎的影响，堆石料的峰值应力比Mf并不是一个常数，而是与平均压力或者是围压相关的，将堆石料的峰值强度统一描述为

(12)

图1 河床典型剖面的剖分网格示意

表1 坝体填筑分级加载过程

式中，Mf0与nf为材料参数，可通过常规三轴试验峰值强度拟合得到。

1.3 流变本构模型[9]

目前混凝土面板堆石坝常采用碾压的施工方式，在分级填筑碾压的过程中，堆石体的流变往往会呈现出“继效特性”。如果考虑到流变过程中的遗传效应，可以采用朱晟等[8]提出的堆石料增量流变模型来进行分析。

假设混凝土面板堆石坝在ζi时刻加载第i级的应力增量，则在该级应力增量下大坝的累积流变量可以表示为(t>ζl)

(13)

(14)

式中，εv(t)、εs(t)分别为i级应力增量作用下的累积体积流变量、累计剪切流变量；Δεvfi、Δεsfi分别为第i级应力增量作用下的最终体积流变量和最终剪切流变量；c、α均为试验常数。

切线流变体积模量Kt和剪切模量Gt分别为

(15)

(16)

式中，kv、nv、ks、ns、Rsf为模型参数，可根据室内三轴流变试验确定。

2 特高混凝土面板堆石坝的三维有限元计算

2.1 坝体网格剖分及施工加载模拟

针对拉哇特高混凝土面板堆石坝，按照前期地质勘察所得到的地质水文条件，可以将其地基作为刚性来考虑。因此依据河床断面的典型特征及坝体分区，采用自动剖分程序进行有限元前处理的过程，同时结合面板分缝特点，沿坝轴线方向共计设置55个剖面，其中河床典型剖面0+188如图1a所示，大坝三维整体剖分网格如图1b所示。坐标系按照常规方向进行设置，以上游指向下游为x轴正方向，以铅直向上的高程增加方向为y轴正方向，以沿坝轴线左岸指向右岸的方向为z轴正方向。拉哇混凝土面板堆石坝总计剖分得到18 518个结点，17 102个单元，其中包括14 154个堆石体单元，964个接触面、面板及竖缝单元，56个周边缝单元。

由于采用考虑流变遗传效应的增量流变模型，故按照大坝实际填筑的分级施工顺序，计算中将荷载施加的过程与设计的填筑分期过程相对应，按逐级加载的方式共分为45个加载级来进行模拟。坝体填筑加载过程及填筑程序分别见表1和图2。

图2 坝体填筑分期(单位：m)

2.2 本构模型参数

结合室内三轴试验结果及拟合过程，确定坝料广义塑性模型计算参数如表2所示。利用上述坝料的室内三轴试验资料，整理邓肯E-B模型参数如表3所示。同时根据室内三轴流变试验结果，整理坝料的增量流变模型参数于表4。

表2 坝料河海统一广义塑性模型参数

表3 坝料邓肯E-B模型参数

表4 坝料增量流变模型计算参数

表5 不同本构模型的堆石体有限元计算极值

2.3 不同本构模型大坝计算结果对比分析

2.3.1 堆石体的应力和变形

选取河床典型0+188纵向剖面进行不同模型计算结果的对比分析，为便于比较，将相同试验资料的应力变形有限元计算极值列于表5。图3为满蓄期该典型剖面的主应力等值线，由图3可知，两种模型计算得到的坝体应力结果数值相近，第一主应力的等值线与坝坡平行，最大值分别为3.82 MPa和3.78 MPa，产生在坝底的轴线附近。统一广义塑性模型计算的第三主应力最大值较大为1.70 MPa，是邓肯E-B模型的1.4倍。

图4为满蓄期典型纵向剖面堆石体的竖向位移等值线。由图4可知，两种模型计算的最大沉降分别为209.9 cm和202.8 cm，约为坝高的0.88%和0.85%，位置均出现在2 600 m高程附近，且位移分布规律保持一致。由表5可得，无论是竣工期还是满蓄期，邓肯E-B模型计算的坝体水平向下游的变形明显较大，满蓄期水平向下游的位移分别为32.2 cm和21.7 cm，约为统一广义塑性模型计算结果的1.5倍。

图4 满蓄期典型纵向剖面竖向位移等值线(单位：cm)

表6为不同工程满蓄期的堆石体沉降极值情况。由表6可知，拉哇混凝土面板堆石坝采用两种不同堆石体本构模型的计算结果与类似工程的实际监测数据数值接近，坝高占比符合一般规律。从堆石体的沉降变形来看，两种模型均可以较好地计算得到大坝变形的分布规律。

表6 不同工程满蓄期的堆石体沉降极值比较

2.3.2 混凝土面板的应力和变形

不同本构模型面板有限元计算极值如表7所示。

图5为满蓄期面板轴向应力的等值线分布。由图5可知，两种模型计算得到的轴向应力结果分布规律类似，面板在水荷载作用下基本处于受压状态，两岸有局部区域存在较小的拉应力。从数值上来看，邓肯E-B模型计算的应力最大值为5.66 MPa，而统一广义塑性模型计算最大值为2.82 MPa，两者均处于2 550 m高程附近。总体而言，统一广义塑性模型计算的面板轴向应力均比邓肯E-B模型小。

图6为满蓄期面板挠度的等值线分布。满蓄期时，邓肯E-B模型和广义塑性模型计算出的面板整体分布规律相近，与类似工程的变形经验相一致[10]。邓肯E-B模型计算的挠度最大值产生在河床典型剖面的2 620 m高程附近，为85.7 cm；而统一广义塑性模型计算的面板中下部坝高部位数值较大，最大值为50.0 cm，可以看出两者计算得到的面板挠度极值差异较大。

图7为满蓄期面板轴向位移的等值线分布。在蓄水后，两种模型计算的面板轴向位移分布规律较为一致，面板左侧向右岸移动，右侧向左岸移动，整体呈现由外部向中心挤压的趋势。两者的最大变形均出现在左右侧中部位置，邓肯E-B模型的轴向位移最大值分别为21.1、23.5 cm，而统一广义塑性模型的轴向位移最大值分别为10.7、11.0 cm，量值上邓肯E-B模型计算值约为统一广义塑性模型的两倍。

3 结 论

利用基于分形理论的堆石料室内缩尺试验成果，分别采用邓肯E-B模型及统一广义塑性模型，对拉哇特高混凝土面板堆石坝进行三维有限元分析计算，研究不同堆石体本构模型的坝体内部应力变形分布规律，得到如下结论：

(1)采用邓肯E-B模型和统一广义塑性模型计算得到的满蓄期堆石体最大沉降分别为209.9 cm和202.8 cm，占坝高的百分比分别为0.88%和0.85%，两者与水布垭、三板溪及洪家渡等类似工程的实际监测数据比较一致，两种模型均可较好地描述面板堆石坝堆石体的变形规律。

表7 不同本构模型的面板有限元计算极值

图5 满蓄期面板轴向应力等值线分布(单位：MPa)

图6 满蓄期面板挠度等值线分布(单位：cm)

图7 满蓄期面板轴向位移等值线分布(单位：cm)

(2)两种本构模型计算的面板变形与应力差异较大，如满蓄时面板挠度最大值分别为85.7 cm和50.0 cm，统一广义塑性模型计算面板的拉应力区域更小，分析原因主要是邓肯模型不能考虑堆石的剪胀性，计算的水平变形更大引起。

(3)两种本构模型计算的坝体应力变形规律相似，计算应力与变形分布规律良好，说明了两种模型应用于特高混凝土面板堆石坝有限元计算分析的可靠性和合理性，计算结果可为后续特高坝工程的设计施工提供依据。