考虑峰后应变软化的蓄水期坝肩岩体变形数值模拟研究

汪维彬

（辽宁润中供水有限责任公司,辽宁 沈阳 110166）

1 工程概况

随着岫岩县经济社会的迅速发展，城区面积不断扩大，城镇供水安全问题也日益凸显出来。为了保证城区供水安全，当地政府计划实施第二水源建设工程。根据相关规划，拟在岫岩县前营子镇西大营村与狄家堡村交界处新建一座蓄水调节中型水库。前营子调节水库的设计库容为3350 万m3。水库大坝为混凝土面板堆石坝，坝轴线长253.00 m，坝顶高程248.00 m，最大坝高34.20 m。水库按照百年一遇洪水标准设计，千年一遇洪水标准校核，设计水位244.23 m，校核水位246.33 m。水库的配套设施为日处理3.5 万m3的净水厂一座，输水管线5 km。工程建成后可以为岫岩县城区、前营子农高区、雅河园区提供生活生产用水，项目计划投资5 亿元。

水库的坝址区属于典型的高山峡谷地貌，左岸山峰高度约450 m，为顺层岸坡，岸坡岩体的岩性主要为二叠纪玄武岩。鉴于坝址区地质条件十分复杂，存在较多的断层、层内和层间错动带、原生和构造裂隙等不良地质构造，而存在这些不良地质构造的左坝肩岩体在蓄水过程中的大坝推力以及渗透压的综合作用下，必将产生机制复杂的应力变形[1]。显然，对上述变形进行合理预测和变形机制探讨对工程蓄水过程的安全性具有重要意义[2]。在岩体应变数值模拟过程中，一般采用理想弹塑性本构模型[3]。由于该模型难以反映岩石峰后应变软化特征，因此模拟计算结果往往不够准确[4]。基于此本文研究基于试验室试验成果数据，构建既能反映岩体强度和变形参数，又能反映峰后劣化变软特征的本构模型，利用应变软化本构模型和摩尔—库伦屈服准则展开对蓄水期坝肩岩体变形数值模拟研究。

2 FLAC3D计算模型

2.1 模型的构建

利用FLAC3D软件构建前营子调节水库左坝肩岸坡三维数值网格模型。模型以坝轴线指向左岸的方向为X轴正方向，以垂直于X轴指向下游的方向为Y轴正方向，以竖直向上的方向为Z轴正方向。模型在X轴和Y轴方向的宽度分别为103 m和45 m，Z轴方向的高度为370 m，其中模型底面高程为80 m，最高处高程为450 m。在建模过程中，将岩体按照风化程度分为Ⅳ类、Ⅲ类和Ⅱ类岩体，分别代表强风化、中风化和微风化玄武岩。针对研究区内结构面的影响，主要考虑F14、F15、F16、F17、F18等大型断层以及 C2、C3、C3-1、C4、C5、C6、C7、C8、C9、C10 等层间错动带。鉴于岩体内各个岩层在形状并不规则，因此选择四面体单元进行模型的网格单元剖分，最终获得60964个计算单元13599个节点[5]。有限元模型示意图见图1。

图1 有限元模型示意图

2.2 模型的边界条件

受到水流侵蚀和地质构造双重作用，坝址区存在十分显著的应力集中现象，因此模型计算要考虑岩体的自重应力[6]。因此，模型的X轴方向的两侧需要施加应力边界条件；Y轴方向的两个侧面与模型的底面施加法向约束条件；模型的上部为自由边界条件；根据现场的实测地质数据资料确定初始应力场。

2.3 本构模型和屈服准则

FLAC3D软件中现有的岩体应变软化模型是基于物理力学参数恒定不变的理想弹塑性模型开发的[7]。事实上，在岩体发生屈服破坏过程中，其内摩擦角、粘聚力以及抗拉强度等物理力学参数是不断变化的[8]。因此，研究中基于岩体材料的应力和应变之间的关系，将上述物理力学参数定义为剪应变或塑性拉应变的分段线性函数，并以此为基础建立应变软化本构模型。在该模型下，岩体的剪切应变和塑性拉应变可以通过增加每一计算步长的硬化参数进行测量。研究中对坝肩岩体进行不同渗透压力条件下的单轴、三轴压缩实验，获得岩石试样的应力—应变全过程曲线，结果见图2。当岩体处于弹性阶段时，其应变大致相当于弹性应变，岩石内部尚未发生塑性变形。在软化阶段岩石的屈服应力逐渐下降，其主要物理力学参数逐渐变弱，岩体逐渐发生塑性变形和屈服破坏，直至下降至残余强度。在残余阶段，岩体的塑性变形持续发展，应变变形不断增大，但是岩体的屈服强度不再发生变化。

图2 应变软化本构模型应力应变关系简图

2.4 计算方案

为了研究前营子调节水库蓄水期左坝肩岩体的应变特征，同时对比考虑和不考虑岩体劣化作用数值模拟结果，本文设计了两种不同的计算方案：方案一不考虑峰后应变软化对岩体劣化作用的影响，主要考虑的荷载为水库蓄水过程中的上、下游水位、扬压力、淤沙荷载。方案二考虑峰后应变软化对岩体劣化作用的影响，主要荷载为水库蓄水过程中的上、下游水位、扬压力、淤沙荷载，同时考虑峰后应变软化对岩体的强度和变形参数的影响。在该方案下，库水为每上升5 m，均需要对相应参数进行重新赋值，直至库水位达到285 m。

3 计算结果与分析

3.1 应力计算结果与分析

利用上节构建的模型，对两种计算方案下的左坝肩主应力进行模拟计算，获得最大主应力和最小主应力分布曲线见图3、图4。由图3可知，在258 m～263 m以及268 m～285 m高程范围内，两种方案的最大主应力计算结果基本一致，在246 m～258 m以及263 m～268 m高程范围内，两种计算方案的差别比较明显，方案二，也就是考虑峰后应变软化特征方案的最大主应力值下降比较显著，257.5 m高程部位的最大主应力值下降幅度最大，为2.48 MPa。

由图4可知，在 246 m～251 m 以及 271 m～285 m 高程范围内，两种方案的最小主应力计算结果基本一致，在251 m～271 m高程范围内，两种计算方案的结果有比较明显的差别，但是方案二相比方案一的下降幅度并不大，仅在252 m、268 m和270 m高程部位下降幅度较大，较方案一分别下降了0.71 MPa、0.54 MPa、0.57 MPa。

图3 最大主应力分布图

图4 最小主应力分布图

3.2 位移计算结果与分析

利用模型计算结果绘制出坝肩岩体位移曲线，见图5～图8。由图可知，坝肩岩体的顺河向位移比较明显，相较于横河向位移和竖向位移明显偏大，同时顺河向位移在不同高程部位的变化曲线与总位移具有基本一致的变化特征。具体而言，高程257 m部位的顺河向位移最大，为20.60 mm，而该部位的横河向位移为8.12 mm，竖向位移为10.71 mm。从坝肩岩体的总位移来看，高程283 m部位由于受到大坝推力的作用产生的位移最小，仅为5.57 mm。在276 m～285 m高程范围内，坝肩岩体的总位移量随着高程的增加而降低；在246 m～276 m高程范围内，坝肩岩体的总位移量呈现出随着高程的增加而先增大后减小的特征，变化范围为11.82 mm～24.73 mm。

图5 横河向位移分布图

图6 顺河向位移分布图

图7 竖向位移分布图

图8 总位移位移分布图

3.3 岩体屈服破坏情况

为了进一步了解前营子调节水库蓄水期左坝肩岩体的屈服破坏情况，对A-A′典型剖面在两种计算方案下的岩体塑性区进行模拟计算，A-A′剖面的位置示意图见图9，计算结果见图10、图11。由图可知，在没有考虑峰后应变软化特征的情况下，坝肩岩体的塑性区主要分布在250 m～290 m高程的坝肩岩体的表面以及边坡内部存在不良结构面的部位，同时，这些塑性区并不连贯，没有向边坡内部扩展的趋势，不会对坝肩岩体的整体稳定造成显著的不良影响。在考虑坝肩岩体峰后应变软化特征的情况下，其塑性区的分布特征并没有显著的变化，但是边坡表面岩体的塑性区范围和分布略有增加，岩体内部不良结构面部位的塑性区范围有明显扩大，但是没有继续向内部扩展的趋势。

图9 A-A＇剖面的位置示意图

图10 方案一塑性区分布图

图11 方案二塑性区分布图

总之，在水库初期蓄水过程中，由于大坝的推力作用，部分坝肩岩体会发生屈服破坏。鉴于坝肩岩体存在峰后应变软化，因此强度和变形参数会有一定程度的降低，这一方面表现为岩体承载力的降低和变形的增大，另一方面会导致周边未破坏岩体分担荷载量的增加，进而造成塑性区范围的扩大，模拟计算结果正好印证了这一点。在工程建成后的蓄水过程中，对坝肩位移变形进行了实际监测，各监测点的数据与本文方案二的计算结果更为接近。因此，采用峰后应变软化模型对前营子调节水库蓄水过程中坝肩岩体的变形情况进行预测分析时，获得的预测结果更为准确。

4 结论

本次研究根据前营子调节水库的地质勘查资料，构建起基于峰后应变软化的三维数值网络模型，进而对坝肩岩体的应力、位移以及塑性区的变化情况进行分析计算，获得主要结论如下：

（1）与未考虑峰后应变软化情况相比，在相同荷载作用下，岩体的最大主应力和最小主应力具有相似的分布特征，但是最大主应力和最小主应力值略有减小。

（2）与未考虑峰后应变软化情况相比，坝肩岩体的位移量具有相似的分布特征，但是顺河向位移、横河向位移、竖向位移以及总位移值均有明显增大。

（3）在考虑坝肩岩体峰后应变软化特征的情况下，其塑性区的分布特征并没有显著的变化，但是边坡表面岩体的塑性区范围和分布略有增加，岩体内部不良结构面部位的塑性区范围有明显扩大。

（4）采用峰后应变软化模型对前营子调节水库蓄水过程中坝肩岩体的变形情况进行预测分析时，获得的预测结果更为准确。