碾压式沥青混凝土心墙坝静力非线性有限元分析

许 涛，侯爰冰，杨旭亮

(1.甘肃省水利水电勘测设计研究院有限责任公司，甘肃 兰州 730030；2.青海民族大学土木与交通工程学院，青海 西宁 810007)

因非连续变形分析与有限元耦合算法既可模拟块体间的不连续变形又可提升块体内应力场的计算精度，对分缝重力坝进行非线性分析[1]。徐远杰等[2]对混凝土面板堆石坝进行三维非线性有限元计算并分析坝体和面板的应力及变形等情况。肖亚子等[3]对沥青混凝土心墙坝采用邓肯E-μ模型进行三维静力非线性分析，将数值模拟的大坝沉降变形结果与原型观测数据比较，判断二者的一致性。赵廷红等[4]运用ABAQUS软件中的混凝土损伤塑性模型对碾压混凝土重力坝进行三维非线性静力计算，并与线性结果比较，验证非线性分析在重力坝应力位移分析中的合理性及适用性。黄焜等[5]采用邓肯-张E-B模型模拟堆石体，利用ANSYS软件建立面板堆石坝三维有限元模型，分析大坝蓄水期的应力和变形，对比坝体位移的计算值与监测值，验证有限元计算的可靠性。孙明权等[6]利用 ANSYS 软件中的APDL开发应用邓肯-张E-B模型，对沥青混凝土心墙堆石坝进行静力有限元计算，从竣工期和蓄水期阶段分别进行应力及变形分析。前人采用不同有限元软件选择邓肯-张本构模型对重力坝和堆石坝的应力和变形进行非线性分析，主要研究坝体的应力与变形，并验证计算结果的合理性，但对于心墙坝的非线性有限元分析相对较少，因此本文结合沥青混凝土心墙坝的坝体分区、施工及加载条件，通过静力非线性有限元分析坝体和心墙工作性态，研究坝体及心墙的应力及变形。

1 静力有限元分析原理

1.1 单元选择

堆石体和心墙采用八结点等参单元，沥青混凝土与过渡层之间采用接触面单元模拟，并采用三棱柱单元及四面体单元以适应边界的不规则变化[7]。

1.2 静力分析方法

按位移求解时，非线性有限元法的基本平衡方程采用增量初应变法迭代求解[8]，则基本平衡方程式为

[K]{Δu}={ΔR}+{ΔR0}

(1)

式中，[K]—整体劲度矩阵；{Δu}—结点位移增量列阵；{ΔR}—结点荷载增量列阵；{ΔR0}—初应变的等效结点荷载列阵。

用中点增量法求解非线性方程组，按筑坝顺序将自重荷载分级，据此进行有限元网格划分。若第i级荷载增量为{ΔR}i，先根据上一级荷载计算终了的{δ}i-1确定弹性常数Ei-1和vi-1，并以此形成劲度矩阵[K]i-1，施加本级荷载增量的一半{ΔR}i/2于结构，用式(2)求位移增量：

[K]i-1{Δδ}i-1/2={ΔR}i/2

(2)

再计算应变增量，累加到上一级终了的应变上，可得本级中点应变{ε}i-1/2，由此确定Ei-1/2和vi-1/2。由式(3)求解本级全荷载对应的位移增量{Δδ}i，据此求解本级的应力与应变增量[9]。

[K]i-1/2{Δδ}i={ΔR}i

(3)

1.3 材料的静力本构关系

1.3.1堆石料及棱体排水料

筑坝堆石料为非线性材料，是沥青混凝土心墙堆石坝的主体，其变形与荷载大小、加荷路径相关，应力应变关系呈非线性。三轴试验研究结果表明，邓肯-张模型对土体应力应变非线性特性能较好的反映，选择邓肯-张E-B模型作为坝体堆石料、过渡料、棱体排水等的本构模型[10]。模型以切线弹性模量Et和切线体积模量Bt作为计算参数，其中切线弹性模量Et为

(4)

切线体积变形模量Bt为：

(5)

对于卸荷或再加荷情况，采用回弹模量Eur进行计算：

(6)

在三维计算中，以广义剪应力q代替(σ1-σ3)，以平均主应力p代替σ3，即

(7)

(8)

破坏偏应力(σ1-σ3)f则根据三维问题的摩尔-库伦准则，表示为

(9)

1.3.2接触面

因过渡层的堆石料与沥青混凝土心墙材料刚度差异较大，为了模拟不同材料之间的相互作用，在有限元分析时设置Goodman接触面单元处理位移不协调问题[11]。

2 坝体静力非线性有限元模型

2.1 三维数值模型

甘肃酒泉某综合利用水库工程为碾压式沥青混凝土心墙砂砾石坝，总库容4910万m3，坝顶高程为1965.40m，坝顶宽10.0m，最大坝高82.40m，坝长556m。上游坝坡为1∶2.5，下游坝坡为1∶2.0和1∶1.8，心墙顶厚度0.5m，至底部渐变为2m。

结合大坝分区、施工和蓄水过程，建立三维有限元模型进行非线性静力分析，将堆石体、心墙、上下游过渡层、排水棱体、覆盖层、接缝等区域剖分为25386个单元。碾压式沥青混凝土砂砾石坝三维有限元模型如图1所示。

2.2 静力计算参数

混凝土护坡在达到破坏强度之前线性关系一般较好，故按线弹性材料处理，弹性材料参数γ=25.0kN/m3，E=30000MPa，ν=0.167；心墙、过渡层、砂砾石层等材料的本构模型采用邓肯-张E-B模型。依据勘察资料、坝料试验资料并参考相似工程，有限元计算中各材料参数取值见表1—2。

表1 接触面参数

表2 坝料的邓肯-张(E-B)模型参数

2.3 坝体填筑加载过程

因坝体施工分层填筑，自重分为10级逐级加载模拟，具体顺序为上游围堰填筑至1905m及1916m；坝体填筑至1912.87、1915.40、1926.00、1930.40、1935.40、1941.40、1950.40、1965.40m。根据度汛及蓄水设计，水荷载分为7级逐级施加模拟，分别为水库蓄水至1899.1、1912.87、1926(度汛水位)、1938、1951.00、1961.60(正常、设计水位)、1962.25m(校核水位)。

3 坝体非线性静力计算结果与分析

3.1 坝体变形

3.1.1典型断面竣工期坝体变形

由竣工期0+430、0+640、0+910横截面的坝体水平位移及竖向位移的等值线可得，横剖面上的水平位移分布规律呈现上游土体位移指向上游，下游土体位移指向下游，与竣工期均质土坝上下游位移分布规律相一致；3个断面最大沉降分别为 25.71、38.95、28.46cm，均发生在坝体的2/3坝高处，具体高程为1940.96、1941.31、1941.49m。图2—3为桩号0+640剖面位移等值线图，向上游最大位移为13.79cm，向下游最大位移为15.68cm，发生在高程1929.26m处。

图2 竣工期0+640截面顺河向位移等值线(单位：cm)

图3 竣工期0+640截面竖直向位移等值线(单位：cm)

图4—5为竣工期坝轴线断面上的顺河向位移和竖向位移等值线，竖向位移分布规律由两岸向河谷方向逐渐增大，最大值为39.32cm，在坝体2/3坝高处。

图4 竣工期坝轴线截面顺河向位移等值线(单位：cm)

图5 竣工期坝轴线截面竖向位移等值线(单位：cm)

3.1.2典型断面满蓄时坝体变形

结合蓄水期0+430、0+640、0+910横截面和坝轴线截面坝体的水平位移及竖向位移分布等值线得水库蓄水后，在水荷载的作用下，上游侧土体向上游的位移减小、下游侧土体向下游的位移增大，坝体沉降变化不大，蓄水期坝体的位移统计见表3。

表3 蓄水期坝体的位移统计

3.2 坝体应力

竣工期坝体0+430、0+640、0+910三个截面的第一、第三主应力等值线与坝坡基本平行，从外向内呈逐渐加大的趋势；蓄水后受水荷载作用坝体应力极值稍增大，但增幅均较小。统计了四个断面竣工期和蓄水后坝体的应力如图6所示。

图6 坝体的应力极值

3.3 坝体剪应力水平

由竣工期和蓄水期坝体剪应力水平等值线分布得，3个剖面中0+640剖面局部剪应力水平最大，分别为0.94和0.96；竣工期与蓄水期0+640截面剪应力水平等值线如图7—8所示，坝体单元大部分区域的剪应力水平合理，均小于1.0，未发现剪切破坏单元。

图7 竣工期0+640截面剪应力水平等值线图

图8 蓄水期0+640截面剪应力水平等值线图

3.4 心墙变形

对心墙坝轴线剖面的变形进行分析，竣工期，心墙上游向最大位移为0.53cm，位于1943.40m高程处；下游向最大位移为1.00cm，位于心墙坝顶处；心墙最大竖向位移为 39.32cm，位于1937.90m高程处。蓄水后，由于受水推力作用心墙顺河向位移增大，其中下游向最大位移为19.40cm，位于1937.90m高程对应2/3坝高处。图9为0+640截面心墙挠曲线，最大挠度与心墙高度的比值为0.245%，小于弯曲破坏时的最大挠跨比7.70%；最大竖向位移为39.40cm，位于1943.40m高程处，因此蓄水对心墙竖向位移影响不大。

图9 0+640位置心墙变形挠曲线

图10 0+640位置心墙垂直应力与水压对比曲线

3.5 心墙应力

竣工期，心墙第一主应力与第三主应力沿高程大致平行分布，且从坝基向坝顶逐渐减小，极值分别为1.58MPa和1.04MPa。蓄水期的应力等值线分布规律与竣工期的应力分布相似，在水荷载作用下心墙应力有所增大，但增幅均较小，蓄水期心墙截面第一主应力和第三主应力等值线见图11—12，极值分别为1.95MPa和1.10MPa。

图11 蓄水期心墙的第一主应力等值线(单位：MPa)

图12 蓄水期心墙的第三主应力等值线(单位：MPa)

心墙的最大主应力均大于单轴压缩试验的抗压强度1.45MPa，但小于三轴试验的最大破坏应力3.65 MPa，考虑到最大应力出现在心墙底部，且心墙两侧受坝壳挤压处于三轴受压状态，应力虽超过了单轴抗压强度，但比三轴试验抗压强度小很多，因此心墙不会发生压缩破坏。图10为0+640位置满蓄期心墙的垂直应力与水压对比曲线，心墙各部位垂直应力均大于蓄水位下的水压，心墙不会发生水力劈裂。

3.6 心墙剪应力水平

心墙剪应力水平等值线分布呈心墙岸坡段较大，竣工期和蓄水期最大值分别为0.44和0.64。水库蓄水后，心墙剪应力水平有所升高，心墙单元的剪应力水平都小于1.0，未发现剪切破坏单元。心墙的位移和应力极值见表4。

表4 心墙的计算结果

4 结论与建议

(1)坝体：竣工期堆石料向心墙上、下游侧变形；蓄水期上游位移减小，下游位移增加。最大压应力在0+640剖面底部中央，最大剪应力小于1.0，不会发生剪切破坏。

(2)心墙：竣工期顺河向位移很小，轴向位移由河谷向中央挤压；蓄水期顺河向位移增大，轴向位移变化不大。应力均小于三轴试验抗压强度3.65 MPa，不产生压缩破坏，不发生剪切破坏及水力劈裂。

通过对大坝应力、位移分布规律研究，可为大坝施工质量控制起指导作用，为同类型工程设计和结构分析提供参考。建议对应力、位移进行实时监测，实时掌握大坝健康状况，确保建设、运行安全。