3.3 心墙堆石坝三维仿真计算
3.3.1 三维有限元分析研究
大坝位于大渡河河湾处，河谷深切，坝体空间效应显著，因此非常有必要对心墙堆石坝进行三维应力应变计算分析。下面主要以河海大学岩土所和河海大学水电学院的研究成果为例进行简要介绍。
3.3.1.1 计算参数
河海大学岩土所的坝料计算参数均采用可研阶段参数，河海大学水电学院心墙料、高塑性黏土、过渡料、反滤料、堆石料均采用技施阶段大坝材料复核试验参数，覆盖层采用可研阶段参数，接触面模型参数根据经验拟定，混凝土作为线弹性材料处理，计算参数见表3.8和表3.9。
表3.8 心墙坝坝料E⁃v模型计算参数（河海大学岩土所）
续表
表3.9 心墙坝坝料E⁃v模型计算参数（河海大学水电学院）
3.3.1.2 计算模型及荷载
混凝土与高塑性黏土之间以及防渗墙周边的泥皮均采用接触面单元进行模拟，防渗墙底部设接触面单元模拟沉渣。坝壳堆石料、过渡料、反滤料和心墙砾石土、地基覆盖层等散粒材料均采用邓肯E⁃v模型，混凝土采用线弹性模型，接触面采用Goodman单元。计算中考虑了上游坝壳、过渡料和反滤料等由于水库蓄水而产生的湿化变形。计算采用分级加荷对施工逐级加荷及蓄水过程进行模拟。
3.3.2 蓄水初期应力与变形分析研究
3.3.2.1 河海大学岩土所计算成果
大坝应力变形的最大值和最小值见表3.10，表中，σz/（γh）为1/3～2/3坝高范围内，河谷段心墙各单元的竖向应力σz与自重应力γh之比值的平均值。该值大小能在一定程度上反映由于拱效应影响而导致心墙竖向应力降低的程度，即反映拱效应的强烈程度。蓄水期坝体的变形与应力等值线图见图3.10～图3.22。
表3.10 大坝主要应力变形值一览表
注 水平位移的“－”表示指向上游；沉降的“－”表示向下；应力的“－”表示拉应力。
图3.10 坝体蓄水期最大断面顺河向水平位移等值线（单位：cm）
图3.11 坝体蓄水期最大断面沉降等值线（单位：cm）
图3.12 坝体蓄水期最大断面大主应力等值线（单位：MPa）
图3.13 坝体蓄水期最大断面小主应力等值线（单位：MPa）
图3.14 坝体蓄水期最大断面应力水平等值线（单位：MPa）
图3.15 蓄水期上游防渗墙大主应力等值线（单位：MPa）
图3.16 蓄水期上游防渗墙小主应力等值线（单位：MPa）
图3.17 蓄水期上游防渗墙顺河向水平位移等值线（单位：cm）
图3.18 蓄水期上游防渗墙沉降等值线（单位：cm）
图3.19 蓄水期下游防渗墙大主应力等值线（单位：MPa）
图3.20 蓄水期下游防渗墙小主应力等值线（单位：MPa）
图3.21 蓄水期下游防渗墙顺河向水平位移等值线（单位：cm）
图3.22 蓄水期下游防渗墙沉降轴向水平位移等值线（单位：cm）
1.坝体变形
竣工期，坝体顺河向水平位移基本关于坝轴线呈对称分布，向上游的水平位移最大值为－43.92cm，向下游的水平位移最大值为37.50cm。蓄水后，坝体向上游的水平位移最大值为－31.69cm，向下游的水平位移最大值为52.87cm，占最大坝高（186m）0.28%，但如考虑覆盖层（76m），则只占0.20%。
竣工期沉降最大值为191.54cm，占最大坝高（包括覆盖层）的0.73%；蓄水期沉降最大值为191.74cm。沉降最大值均发生在心墙底部的高塑性土顶部，这主要是高塑性土的变形较大造成的。
2.坝体应力
受拱效应的影响，心墙内大主应力比过渡层的应力明显降低；心墙中小主应力与反滤、过渡及堆石料中的小主应力之间过渡平缓，说明拱效应对小主应力影响不明显。
表3.10显示，河谷段心墙各单元的竖向应力与自重应力（γh）之比值σz/（γh）竣工期为0.695，蓄水后为0.738。该值越小，表明由于拱效应引起的心墙应力减小越严重。
竣工期坝体应力水平总体并不高，堆石体在0.2～0.4范围，心墙大部分在0.4～0.5之间。蓄水期，心墙和下游堆石内应力水平略增大，而上游坝壳的应力水平明显增大，心墙没有剪切破坏区。竣工期和蓄水期大部分高塑土单元的应力水平较高，但大部分范围小于1.0，表明高塑土没有发生剪切破坏。上游防渗墙上游侧的高塑土在蓄水后应力水平达到1.0。这主要是水压力导致防渗墙向下游位移，使得这部分高塑土小主应力减小，发生剪切破坏，这正说明了易发生大剪切变形部位设置高塑土的必要性。因此上游防渗墙上游侧的高塑土是防渗相对薄弱的范围，应予以关注。
3.防渗墙应力变形
（1）上游防渗墙。
1）上游防渗墙竣工期大主应力最大值为37.71MPa，蓄水期大主应力最大值为45.47MPa。上游防渗墙大部分区域在蓄水后大主应力减小，但防渗墙底部由于插入基岩导致弯应力增大，从而蓄水后的最大值比竣工期大。
2）上游防渗墙竣工期和蓄水期防渗墙大部分区域的小主应力为正值，即为压应力，局部边角部位出现了拉应力，但拉应力最大值不超过3MPa。蓄水引起了防渗墙小主应力绝对值增大，对大部分区域，由于小主应力是上下游方向，水压力刚好是这个方向，因此，水压力引起小主应力增大。但对防渗墙下部，同样是由于弯应力，使得小主应力的拉应力更大。
3）上游防渗墙竣工期向上游的水平位移最大值为－5.26cm，位于防渗墙顶部；向下游的水平位移最大值为5.89cm，位于防渗墙中间偏上部的位置。蓄水后防渗墙在其上游面水压力的作用下，向下游方向发生较大位移，蓄水期向下游的水平位移最大值为39.21cm。
4）上游防渗墙沉降最大值位于建基面附近，竣工期为－15.76cm，蓄水期为－15.53cm。
5）上游防渗墙底沉渣单元竣工期最大压缩量为9.02cm，蓄水期最大压缩量为8.64cm。
（2）下游防渗墙。
1）下游防渗墙竣工期大主应力最大值为35.81MPa，蓄水期大主应力最大值为43.60MPa。和上游防渗墙一样，对大部分区域，蓄水是引起防渗墙大主应力减小的，但由于弯应力作用，防渗墙下部压应力增大。
2）下游防渗墙大部分区域的小主应力为正值，即为压应力。竣工期和蓄水期的小主应力和上游防渗墙相近。
3）下游防渗墙竣工期向上游的水平位移最大值为－3.24cm，向下游的水平位移最大值为1.48cm；蓄水后防渗墙在其上游面水压力的作用下，向下游方向发生较大位移，蓄水期向下游的水平位移最大值为34.73cm，位置靠近防渗墙顶部。
4）下游防渗墙沉降最大值位于防渗墙上部中间位置，竣工期为－15.66cm，蓄水期为－15.90cm。
5）下游防渗墙底沉渣单元竣工期最大压缩量为11.65cm，蓄水期最大压缩量为11.94cm。
4.廊道应力
竣工期廊道坝轴向正应力均为压应力，河谷中间断面上的正应力最大，最大值为12.6MPa。廊道底板上边缘的顺河向正应力出现了较大拉应力，且以河谷中心段最大，达到－3.9MPa，靠岸坡段，廊道底板的顺河向拉应力也达到了－2.5MPa。尽管河谷中心段廊道底板拉应力较大，但仅出现在底板的上缘，下缘仍为压应力。竖直向正应力均为压应力，也是河谷中间部位最大，最大值为9MPa。廊道底板剪应力也以河谷中心段最大，达4.8MPa，见图3.23。
图3.23 0＋272断面蓄水期廊道应力（单位：MPa）
蓄水后廊道顺河向应力σy和竖直向应力σz基本规律和竣工期基本一致，大小也接近，但坝轴向应力σx变化明显，主要是上游水压力作用下，廊道河谷中间段向下游位移比两端位移大，导致廊道下游侧σx减小，而上游侧σx增大，最大达17.8MPa。
5.心墙抗水力劈裂研究
计算采用总应力法判定水力劈裂的发生，从上游水压力与心墙竖向应力或中主应力比值的等值线（图3.24和图3.25）可以看出，上游水压力与心墙竖向应力或中主应力的比值均小于1.0，因此，可认为心墙的抗水力劈裂性能是有保证的。
图3.24 蓄水期上游水压力与心墙中应力之比值等值线
图3.25 蓄水期上游水压力与心墙竖向应力之比值等值线
3.3.2.2 河海大学水电学院计算成果
竣工期与蓄水期坝体的变形与应力极值见表3.11，蓄水期坝体的变形与应力等值线图见图3.26～图3.31。
表3.11 三维非线性有限元计算坝体变形与应力极值
图3.26 蓄水期0＋256断面水平向位移（E⁃v模型，单位：cm）
图3.27 蓄水期0＋256断面竖直向位移（E⁃v模型，单位：cm）
图3.28 蓄水期0＋256断面第一主应力（E⁃v模型，单位：MPa）
图3.29 蓄水期0＋256断面第三主应力（E⁃v模型，单位：MPa）
图3.30 蓄水期0＋256断面剪应力水平
图3.31 蓄水期高塑性区黏土的剪应力水平
1.坝体的变形、应力与应力水平
（1）坝体变形。竣工期，坝体向上游最大水平位移为32.5cm，位于上游堆石高程730.00m附近；向下游移动的最大水平位移为37.6cm，位于下游主堆石高程720.00m附近。坝体最大竖向位移为216.2cm，发生在高程740.00m附近的心墙部位。水库蓄水后，坝体向上游最大水平位移减小到21.5cm，位于上游堆石内高程700.00m附近；向下游位移增大到69.9cm，位于高程760.00m附近的心墙内。最大竖向位移为234.5cm，发生在高程740.00m附近的心墙部位。
（2）坝体应力。竣工期坝体第一主应力最大值为3.30MPa，第三主应力最大值为1.32MPa，均发生在坝体与基岩接触部位；蓄水期第一主应力最大值为3.58MPa，第三主应力最大值为1.65MPa，位置与竣工期基本相同。由于水荷载的作用，导致坝体内部心墙、地基防渗墙和廊道共同组成的防渗体的上、下游侧的应力存在不连续现象。
（3）坝体剪应力水平。竣工期上、下游反滤料和心墙结合处的应力水平较高，最大值为0.6；蓄水期上游反滤料区与心墙结合处应力水平较高，最大值为0.93，但没有发生剪切破坏。蓄水期高塑性黏土区的剪应力水平较高，特别是在高塑性黏土与防渗墙与廊道接触区域，单元剪应力水平达到1.0，但是由于较强的塑性特性，该区域黏土与心墙交接部位未发生剪切破坏，说明高塑性黏土区对防渗墙与心墙交接位置的变形与应力具有较强的改善作用。
2.心墙的拱效应与水力劈裂问题
（1）心墙的拱效应。由于心墙的变形模量低于坝壳堆石体的变形模量，从而使坝壳堆石体的应力远大于心墙，发生拱效应。图3.32为蓄水时心墙轴线剖面竖向正应力与该处土柱压力比值［σy/（γsh）］分布云图，γs为心墙砾石土的容重。可见，心墙1/3以上其比值基本处于0.5～1.0之间，在接近坝体顶部拱效应很小，在高程840.00m以下的心墙拱效应现象随着土柱高度的增加而变得越来越明显，拱效应最大处位于心墙底部与高塑性黏土接触位置，处于0.1～0.5之间。
图3.32 蓄水期心墙拱效应σy/（γh）比值分布图
（2）水力劈裂问题。考虑到小主应力方向一般为上下游方向，其作用面与水荷载作用方向垂直，现阶段倾向于以中主应力估算和判断心墙可能发生的水力劈裂。图3.33为蓄水时心墙上游面的水压力与中主应力比值分布云图［图中仅示出σ2/（γwh）＜1.0的区域］，γw为水的容重。当σ2/（γwh）＜1.0的时候，心墙可能发生水力劈裂破坏。由图可见其比值小于1.0主要存在于心墙的两个区域，即两岸坝段的上游水位附近与河床坝段的坝基面附近。
图3.33 蓄水期心墙上游面σ2/（γwh）比值分布图
3.防渗墙的变形与应力
（1）防渗墙的变形。竣工期上游防渗墙顺河向水平位移大部分指向下游，极值为3.9cm，位于右岸防渗墙顶部，而河床段防渗墙顶部顺河向位移指向上游，极值3.5cm；竖直向位移极值为14.6cm，位于河床防渗墙顶部。蓄水期水压作用使得防渗墙的顺河向位移进一步向下游方向发展，极值为27.1cm，竖向位移极值为16.9cm。由于顶部廊道的作用，竣工期下游防渗墙顺河向位移在防渗墙中部指向上游，极值为3.8cm，竖向位移极值为18.2cm。蓄水期水压作用使得防渗墙的顺河向位移进一步向下游方向发展，极值为25.5cm，竖向位移极值为19.4cm，见图3.34～图3.37。
图3.34 蓄水期上游防渗墙水平向位移（E⁃v模型，单位：cm）
图3.35 蓄水期上游防渗墙竖直向位移（E⁃v模型，单位：cm）
图3.36 蓄水期下游防渗墙水平向位移（E⁃v模型，单位：cm）
（2）防渗墙的应力。竣工期上游防渗墙第一主应力都为压应力，极值为27.63MPa。防渗墙的第三主应力极值为1.46MPa，在左、右两岸的端部，防渗墙存在局部的拉应力区，极值达3.38MPa。水库蓄水后，在水压力作用下，防渗墙的大主应力一般大于0，其中第一主应力极值为31.59MPa；第三主应力的压应力极值为1.62MPa，左、右两岸对应的防渗墙，局部拉应力区的数值有所加大，极值为3.78MPa。竣工期下游防渗墙第一主应力都为压应力，极值为35.23MPa。防渗墙的第三主应力为拉应力，极值为1.03MPa。水库蓄水后，第一主应力极值为35.7MPa；而第三主应力的应力状态有所改善，在河床段防渗墙中部出现压应力区，极值为2.42MPa，在防渗墙两端还是有拉应力区，极值为1.57MPa，见图3.38～图3.41。
图3.37 蓄水期下游防渗墙竖直向位移（E⁃v模型，单位：cm）
图3.38 蓄水期上游防渗墙第一主应力（E⁃v模型，单位：MPa）
图3.39 蓄水期上游防渗墙第三主应力（E⁃v模型，单位：MPa）
图3.40 蓄水期下游防渗墙第一主应力（E⁃v模型，单位：MPa）
图3.41 蓄水期下游防渗墙第三主应力（E⁃v模型，单位：MPa）
4.廊道的变形与应力
竣工时，坝内廊道向下游方向移动，最大位移为8.6cm；最大沉陷为19.8cm。水库蓄水后，其水平顺河向、沉陷极值分别增大为31.1cm和21.2cm。廊道在两岸位置发生了扭曲现象，出现了明显的错动变形。竣工期最大错动3.8cm，位于右岸廊道接缝的顶部位置；水库蓄水后，错动值增大到4.1cm。
竣工期，廊道大主应力最大值为12.82MPa，位于河床中央与防渗墙连接的部位；廊道底板和顶拱的小主应力处于受拉状态，但其拉应力基本控制在1.0MPa以内，最大拉应力区分布在廊道的两端，极值达1.97MPa，位于河床廊道与左岸基岩交接的端头位置。水库蓄水后，廊道两岸端头由于岸坡基岩的约束作用，局部出现了明显的应力集中现象。其迎水面侧墙的主拉应力有增大的趋势，最大达3.66MPa，位于左岸廊道端头。背水面侧墙位置，大主应力增大，如左岸廊道下游侧的最大值达18.3MPa。
5.计算成果分析
综合上述三维静力计算成果，可以得出以下主要结论。
（1）大坝的应力和变形分布规律基本合理，成果较为接近，大坝最大沉降均发生在水库满蓄期，坝体沉降最大值为191.7～234.5cm，占大坝与地基覆盖层总高度的0.73%～0.90%；最大水平位移值为52.9～69.9cm。
（2）心墙内没有出现拉应力，应力水平也较低，出现水平拉裂缝和坝肩处的竖向拉裂缝的可能性不大；高塑性土大部分单元的应力水平较高，在心墙、防渗墙、墙顶廊道的接头部位设置高塑性黏土区，改善了大坝的应力状况，使得剪切破坏区大大减小，达到了控制该部位应力与变形的目的。
（3）大坝心墙的拱效应明显，其中坝轴线剖面大部分区域的铅直向应力仅为其土柱重量的70%左右；对蓄水期心墙上游面水力劈裂可能性进行分析得到：在心墙上游面两岸与蓄水位交接面附近，局部单元可能发生水力劈裂。
（4）两道防渗墙的应力和变形状态较为接近，上游防渗墙最大水平位移值为27.1～39.2cm，最大沉降值为15.5～16.9cm，最大主应力值为31.59～45.47MPa；下游防渗墙最大水平位移值为25.5～34.7cm，最大沉降值为15.9～19.4cm，最大主应力值为35.7～43.6MPa。
（5）由于廊道位于防渗墙顶部，导致廊道的拉应力也较大，主拉应力最大值达3.66～3.9MPa。
由前面非线性有限元模拟的结果可得到以下结论。
（1）大坝竣工期和蓄水期的应力和变形分布规律基本合理，大坝最大沉降均发生在水库满蓄期，其中二维固结计算值为246.6cm，三维非线性有限元计算值为234.5cm，三维弹塑性有限元计算值为148.4cm，占大坝与地基覆盖层总高度的比例分别为0.95%、0.9%和0.6%。
（2）坝基防渗墙应力最不利工况都发生在水库满蓄期。其中上游墙在两岸端头伸入基岩，该部位的拉应力最大，如利用邓肯模型计算的主拉应力极值达3.78MPa。下游墙由于墙顶廊道上部心墙土压力的传递作用，河床防渗墙中部的压应力最大，如利用双屈服面模型计算的第一主压应力极值达43.11MPa。且由于廊道位于防渗墙顶部，导致廊道的拉应力也较大，在两岸端头部位设置错缝时主拉应力最大值达3.66MPa。
（3）大坝固结有限元计算的结果表明，考虑心墙中孔隙水压力的产生和消散过程后，蓄水期心墙不会发生水力劈裂。
（4）瀑布沟大坝心墙的拱效应明显，其中坝轴线剖面大部分区域的铅直向应力仅为其土柱重量的70%左右；在心墙上游面两岸与蓄水位交接面附近，局部单元可能发生水力劈裂。
对于坝体应力变形来说，总体上比较均匀，坝壳对心墙存在一定的拱效应，这种拱效应导致心墙上游部位的小主应力减小，但河床部位心墙的小主应力仍然为压应力，表明该部位心墙受力状态良好。岸坡对心墙的约束对心墙的应力影响很大，高塑性黏土良好的塑性是协调心墙与岸坡之间变形的关键。心墙与两岸相接部位小主应力较低，蓄水期上游表面出现一定范围的拉应力，该部位可能出现拉裂缝，但心墙内部偏中间与下游部分小主应力仍然是压应力，这一部位设计厚度大于心墙中部，是合理的设计措施。
根据大坝三维应力应变分析的结果，在技施设计阶段，采取了以下工程措施：
（1）在坝体应变较大部位设置高塑性黏土。尽管根据有限元分析结果，竣工期心墙内的σ1均大于正常蓄水位时的水压力，心墙的水平面上不会产生水力劈裂，但在坝体应力应变比较大的部位仍尽可能设置高塑性黏土。例如在左、右岸坝坡与河床最大剖面的心墙下半部，σ2与水压力相比，当不计心墙土料自身的抗拉强度时，抗水力劈裂能力裕度不大，故在心墙与两岸基岩接触面上铺设3m厚的高塑性黏土。在防渗墙和廊道周围均铺设黏土。
（2）加强心墙上下游的反滤层。在心墙上下游侧及心墙与基础覆盖层间设两层比较宽厚的反滤层，反滤层外再设过渡层。上下游各设两层反滤层，上游两层各厚4m，下游两层厚6m。不仅防止渗流出口的管涌和接触冲刷，一旦出现裂缝还可能促使心墙裂缝自愈，同时在坝壳与心墙间起过渡作用，改善心墙的应力条件，减少心墙拱效应。
（3）在混凝土盖板与砾石土心墙之间设置高塑性黏土区，以协调两者之间的变形。
（4）河床基础覆盖层进行了固结灌浆，以改善坝基接触条件并提高地基变形模量。
（5）砾石土心墙靠两岸部位适当加宽，以增加渗径、改善两岸部位接触条件。
（6）将上游防渗墙上游侧心墙底部5m厚的黏土调整为1m，下游防渗墙下游侧心墙底部5m厚的黏土取消，以减小坝体沉降。
（7）大坝坝顶竣工后的预留沉降超高值不小于2.0m。