基于两种本构模型的堆石坝应力变形比较研究

朱俊高，单一峰*，郑惠峰，刘 忠

(1.河海大学 岩土力学与堤坝工程教育部重点实验室，江苏 南京 210098；2.河海大学 岩土工程科学研究所，江苏 南京 210098；3.中国电建集团华东勘测设计研究院有限公司，浙江 杭州 330110；4.黄河水利委员会黄河水利科学研究院，河南 郑州 450003)

近年来，随着堆石坝筑坝高度的增加，有限元分析在坝体设计中起到了不可或缺的作用。有限元计算对选取的本构模型及其模型参数具有较大的依赖性，不同的本构模型在计算堆石坝的应力变形时，计算结果相差较大[1-3]。因此，在有限元计算中，本构模型的合理选取至关重要。在我国堆石坝的有限元分析中，弹性非线性模型(如邓肯模型、K-G模型)和弹塑性模型(如南水模型、椭圆-抛物双屈服面模型)都得到了不同程度的应用。邓肯E-B模型[4]，参数确定简单，模型参数的物理意义明确，为弹性非线性模型的代表模型，被广泛应用于堆石坝的有限元分析中[5-7]，但该模型基于广义胡克定律，不能反映剪胀性。殷宗泽提出的椭圆-抛物双屈服面弹塑性模型[8]，综合了单屈服面模型的优点，能较好反映土的剪缩剪胀性及复杂应力路径对土体变形的影响，具有较好的实用性。

目前，已有学者比较了这几类模型在有限元计算中的差异[9-10]，如赵晓龙[11]等比较了邓肯-张模型和椭圆-抛物双屈服面模型土石坝应力变形的计算结果。朱俊高[12]等人结合糯扎渡心墙坝实例，分析比较了邓肯-张模型和邓肯E-B模型在有限元计算中的差异。结果表明使用不同模型计算所得到的坝体变形和应力分布均有所区别。尽管上述研究对比众多，但少有学者对邓肯E-B模型和椭圆-抛物双屈服面模型计算的堆石坝应力变形进行全面的分析比较。基于此，本文以两岔河心墙堆石坝为例，分别采用邓肯E-B模型与椭圆-抛物双屈服面模型进行三维有限元计算，并对计算结果进行了对比分析，探究两种模型在应力和变形的差异，为心墙堆石坝有限元计算本构模型的选取提供了参考依据。

1 工程概况和有限元模型

两岔河水库位于四川省凉山州会东县。挡水建筑物为粘土心墙堆石坝，坝顶长度247.00 m，宽10 m，最大坝高74.50 m。心墙顶宽4.0 m，上、下游坡比1∶0.25，心墙最大底宽39.25 m，心墙与两岸坝肩混凝土垫层接触部位采用厚度2.0 m的高塑性粘土过渡。坝基覆盖层防渗采用一道混凝土防渗墙，墙厚1.2 m。防渗墙顶通过廊道与防渗心墙连接，底部嵌入基岩。坝体分区见图1。

坝体各分区的邓肯E-B模型和椭圆-抛物双屈服面模型的计算参数分别见表1和表2。其中邓肯E-B模型参数通过常规三轴试验确定，双屈服面模型参数通过河海大学岩土所研制程序[13]由E-B模型参数转换得到，具体如下：由已知E-B模型参数计算出对应土石料应力应变曲线族，再用最优化方法拟合，从而获得椭圆-抛物双屈服面模型参数。

大坝共划分59 480个结点、56 848个单元。单元网格划分时，对土石料、混凝土等单元，大部分用8结点6面体单元，少数用6结点5面体、4结点4面体等单元过渡，边界条件设置为底面xy双向约束，坝体两岸yz双向约束，同时采用接触面单元模拟各种材料边界力与变形的协调关系。为了更好地模拟现场实际工程情况，计算采用分级加荷对大坝逐级施工及蓄水过程进行有限元计算，分19级进行加载，其中前15级为施工加荷，后4级为蓄水加荷，不考虑渗流的影响，仅作为附加水荷载。大坝三维有限元网格见图2，河谷段最大断面网格见图3。

表1 邓肯E-B模型计算参数Tab.1 Parameters of Duncan E-B model

表2 椭圆-抛物双屈服面模型计算参数Tab.2 Parameters of ellipse-parabolic double yield surface model

图2 大坝三维有限元网格Fig.2 Three-dimension finite element mesh

图3 河谷段最大断面网格Fig.3 The maximum cross section of valley segment

2 有限元计算结果及分析比较

本文采用河海大学岩土工程研究所研制的TDAD三维有限元软件，进行两种不同本构模型的应力变形计算，计算整理得到的最大沉降和水平位移见表3。为研究应力变形分布规律，整理出了(正常高水位)满蓄期坝体及防渗墙的沉降和顺河向水平位移等值线，见图4—图7，以及满蓄期大主应力和小主应力等值线，见图8—图11。

2.1 坝体位移分析

由表3可知，邓肯E-B模型和双屈服面模型计算的最大沉降量很接近，双屈服面模型计算的上下游水平位移小于邓肯E-B模型计算结果，坝轴向整体水平位移(左右岸位移之和)较邓肯E-B模型也偏小。总体上，两种模型的沉降差异很小，水平位移差异稍大。说明了两种模型在反映土石料性质能力上存在差异。

表3 满蓄期坝体沉降位移结果Tab.3 Settlement and displacement of dam after water storage

注：表中正值代表位移指向坝体下游或沿坝轴指向右岸，负值代表位移指向坝体上游或沿坝轴指向左岸。

图4给出了满蓄期坝体沉降等值线，可以看出，两种模型计算的等值线图分布规律基本相似，坝体沉降分布呈对称形状。总体上坝体沉降不大，坝体最大沉降在同一区域，位于心墙自底向上约1/3坝高处中心位置。满蓄期两种模型计算的最大沉降量分别是113.6、116.2 cm，分别占最大坝高119.5 m(含覆盖层)的0.95%和0.97%。

经收集各大坝体的实际沉降数据资料发现[14]，一般情况下，满蓄期的堆石坝最大沉降值约为坝高的0.52%～1.36%，最大位移比为0.029～0.749，邓肯E-B模型和双屈服面模型计算的位移比分别是0.294与0.240，均在合理范围内。同时需要注意的是，在使用邓肯E-B模型进行有限元计算的过程中，若泊松比μ大于0.5，劲度矩阵发生异常，无法进行下一步计算，因此该模型不能反映土体的剪胀性，并且，该模型也无法反映土体压缩及剪切变形的交叉影响。综上，使用可以反映剪胀性的双屈服面模型进行有限元计算所得到的应力变形结果在理论上更加合理。

从图4还可以看出，廊道与高塑性黏土区内的等值线分布比较密集，说明此处沉降梯度较大，这是由于廊道及防渗墙的支撑顶托作用导致该处产生较大的剪切变形。从沉降等值线图还可以看出，覆盖层的沉降等值线分布较为稀疏，相比于坝体沉降并不大。

图5给出了满蓄期坝体顺河向水平位移等值线，其中，正值表示向下游位移，负值表示向上游位移。可以看出两个模型计算的水平位移分布规律相近，由于土体自重作用，坝体在沉降的同时会产生向上下游的水平位移，但蓄水后由于水压力向下游最大水平位移都明显大于向上游的位移，蓄水后坝体向上游的最大水平位移位于靠近原地面线的上游坝壳内。满蓄期坝体向下游的最大位移，用邓肯E-B模型计算为30.2 cm，用双屈服面模型计算为27.9 cm；满蓄期坝体向上游的最大位移，用邓肯E-B模型计算为9.3 cm，用双屈服面模型计算为7.1 cm。相比邓肯E-B模型，双屈服面模型计算的上下游最大位移均更靠近中轴线。

图6、图7分别给出了满蓄期防渗墙的沉降等值线和满蓄期顺河向水平位移等值线。图中可见，防渗墙在河谷中央顶部沉降最大，往左右两边及底部逐渐减小，邓肯E-B模型和双屈服面模型计算的最大沉降分别为1.8和2.4 cm，双屈服面模型计算的沉降较大，等值线变化梯度更大。蓄水后，防渗墙受覆盖层土体侧向挤压均表现为向下游的位移，呈凹向上游变形。防渗墙最大水平位移位于近河谷中央的位置，往左右两岸、顶部及底部位移逐渐减小。比较两种模型计算结果，邓肯E-B模型计算结果稍大，最大顺河向水平位移为5.8 cm。

图4 满蓄期沉降等值线图(单位：cm)Fig.4 Contour lines of settlement after water storage

图5 满蓄期顺河向水平位移等值线(单位：cm)Fig.5 Contour lines of horizontal displacement along the river after water storage

图6 满蓄期防渗墙沉降等值线(单位：cm)Fig.6 Contour lines of settlement of cutoff wall after water storage

图7 满蓄期防渗墙顺河向水平位移等值线(单位：cm)Fig.7 Contour lines of horizontal displacement along the river of cutoff wall after water storage

2.2 坝体应力分析

图8给出了满蓄期坝体典型横断面大主应力等值线。无论邓肯E-B模型还是双屈服面模型，满蓄期坝壳内大主应力的等值线分布基本平行于坝坡，由外向内大主应力逐渐增加。受拱效应的影响，心墙内大主应力等值线比过渡层的大应力等值线更为稀疏[15]，尤其在心墙的中上部更为显著。在心墙下部由于心墙及高塑性粘土受到下部廊道及防渗墙的顶托作用，心墙向下的沉降受到影响，使得这种拱效应作用下降，因此在心墙下部，心墙内大主应力较过渡层降低不明显。由于双屈服面模型考虑了塑性变形，能较为全面地反映土体的剪胀剪缩性及剪切变形的交叉影响，该模型下计算得到的坝体应力分布比较合理。

坝体典型横断面蓄水工况下的小主应力等值线见图9。图中可见，小主应力在坝体内的分布也近似呈平行于坝坡的形式，由外向内逐渐增加。由于拱效应，心墙内小主应力与过渡层的应力相比数值也有所下降，但与大主应力相比拱效应作用不明显。邓肯E-B模型计算的坝壳小主应力相比双屈服面模型的结果降低得更为显著。上游坝壳的小主应力由于蓄水影响明显减小，心墙内靠近上游侧小主应力也有所减小，但仍为正值(此处正值代表压应力，负值代表拉应力)，说明心墙内并未出现拉应力。

对河谷段设有廊道和防渗墙的大坝，由于蓄水时在廊道和防渗墙上下游面存在水头差，廊道和防渗墙将产生向下游顺河向水平位移，而廊道两端和防渗墙周边受到基岩的约束，有可能在防渗墙和廊道部分区域产生拉应力，故需要特别关注。图10给出了满蓄期防渗墙中线处大主应力等值线图，可以看出，大主应力等值线基本呈对称分布，由防渗墙顶往左右两边及底部逐渐减小。大主应力均为压力，防渗墙墙顶中部区域应力值最大，邓肯E-B模型和双屈服面模型计算值分别为15.5和20.5 MPa。

图11为满蓄期防渗墙中轴面小主应力等值线，由于受到周边基岩的约束作用，小主应力在左右岸坡边角部出现一定的拉应力。邓肯E-B模型和双屈服面模型计算的拉应力值分别为1.56和1.89 MPa。双屈服面模型计算的拉应力的值及变化梯度均大于邓肯E-B模型，且左岸拉应力差异较右岸更为明显。

图8 满蓄期大主应力等值线图(单位：MPa)Fig.8 Contour lines of major principal stress after water storage

图9 满蓄期小主应力等值线图(单位：MPa)Fig.9 Contour lines of minor principal stress after water storage

图11 满蓄期防渗墙中线处小主应力等值线图(单位：MPa)Fig.11 Contour lines of minor principal stress at the midline of cutoff wall after water storage

由大小主应力等值线图可以看出，两种模型计算的应力差异较大，故整理了截断面中轴线的应力计算结果，对应力作进一步分析比较。图12给出了心墙与防渗墙在最大断面中轴线位置上的竖向应力σz和顺河向应力σy的分布，可以看出，双屈服面计算的应力均比E-B模型结果大。这是由于双屈服面模型计算的沉降值大，而水平位移又较小，坝体材料挤压效果更明显。对于心墙部分，其双屈服面模型计算的内部应力越大，对其防渗效果的发挥更好，如果心墙内应力太低，则可能产生心墙裂缝，甚至水力劈裂。对防渗墙，双屈服面模型计算的应力比E-B模型的结果要大5 MPa左右，尽管都没有超过C30混凝土强度，但对于高坝或更深厚的覆盖层，还是应该引起注意，E-B模型计算的防渗墙应力可能偏小。

3 结论

1)邓肯E-B模型与椭圆-抛物双屈服面模型，两种模型计算的坝体沉降十分接近，坝体水平位移有所差异。与双屈服面模型相比，邓肯E-B模型计算的上下游水平位移稍大，其坝轴向水平位移也比双屈服面模型计算结果大。邓肯E-B模型与双屈服面模型计算的最大沉降分别占最大坝高(含覆盖层)的0.95%和0.97%，位移比分别是0.294和0.240，与多个堆石坝工程监测数据比较，两种模型计算的结果均符合实际工程沉降规律。

图12 满蓄期河谷中心坝轴线位置应力分布对比(单位：MPa)Fig.12 Stress comparison diagram of dam body truncated surface during full storage period

2)对于两种模型计算得到的大小主应力，邓肯E-B模型受拱效应的影响较为明显，相同高程下，心墙内的大小主应力比过渡层的应力下降得更为明显。心墙内的小主应力均为正，心墙处于受压状态。

3)双屈服面模型计算的沉降较大，水平位移较小，使得坝体侧向挤压效果比邓肯E-B模型更明显，导致防渗墙大小主应力及应力变化梯度均大于邓肯E-B模型。受周边基岩的约束作用，防渗墙中轴面小主应力在左右岸坡区域均出现拉应力。两种模型计算的防渗墙应力都在其混凝土的容许应力范围内，其中左岸拉应力差异较为明显，双屈服面模型计算的拉应力的值及变化梯度较邓肯E-B模型大。

邓肯E-B模型和椭圆-抛物双屈服面模型计算的心墙及防渗墙处的应力变形有明显差异，可供设计综合参考。基于双屈服面模型能反映土体压缩及剪切变形的交叉影响，且采用其计算结果设计更为保守。综合考虑，在满蓄期工况条件下，采用双屈服面模型进行有限元计算更加合理准确。本文研究成果可为类似工程设计时的模型选用提供一定的参考依据。