沈 辉，罗先启，郑安兴，毕金锋，翁永红

（1.上海交通大学 船舶海洋与建筑工程学院，上海 200240；2.长江水利委员会 长江勘测规划设计研究院，武汉 430010）

1 引 言

拱坝坝肩岩体稳定性分析是一个受多种因素影响的地质力学问题。通常需要在地质勘探、现场调查研究的基础上，通过地质条件综合分析建立地质体的概念模型，通过地质力学模型试验或数值模拟方法模拟地质体力学环境演化的全过程，以达到预测分析的目的。

随着理论分析的成熟和计算机技术的发展，在大型水利水电工程及岩土工程等领域，数值模拟方法以其高效低成本的特点，成为地质力学问题分析的关键技术手段。从目前的研究情况来看，地勘资料的详细程度及研究问题的不同，决定了数值模型的精细度及规模，也使得分析人员根据实际情况选择合理的建模及分析软件[1]。由于地质模型建立、空间单元划分等前处理工作烦琐，往往不得不采用简化地质模型的做法，对工程建筑物近域地基中复杂的地质条件以及坝址河谷两岸实际地形的影响考虑较少。文献[2]采用ANSYS建立拱坝及概化密集节理模型，探讨密集节理对坝肩稳定的影响问题；文献[3]以大型通用有限元程序ADINA为计算平台,对比分析坝肩边坡卸荷后，不同地震加速度及不同锚固长度和锚固吨位在卸荷条件下关键点的位移和应力以及塑性区的变化规律；文献[4]采用考虑开挖过程的有限元强度折减法分析锦屏一级水电站左岸坝肩开挖边坡在施工期的稳定性。也有许多学者[1,5]将三维地学建模与数值模拟进行耦合，让数值模拟直接继承和利用三维地质详细的地勘数据。但需要掌握三维地学建模系统的格式和算法，编制计算程序复杂，不利于推广。

本文以乌东德超高拱坝坝肩稳定性分析数值模型为例，对可研阶段地质勘探数据、二维地质剖面和工程设计数据进行预处理，提取包括地形地貌、区域地层、断层构造、岩溶系统等自然地质对象界限，以及开挖面、拱坝等人工对象边界。基于Hypermesh软件对各类对象进行插值、拟合和几何建模，利用其强大的网格参数交互调整功能以及独立单元与节点的处理功能，建立拱坝联合坝址区岩体的三维数值精细模型，为坝址区地应力反演分析，坝体载荷、工程开挖等因素影响下坝肩稳定性分析提供了良好的分析平台。

2 工程地质概况

乌东德水电站大坝为混凝土双曲拱坝，坝顶高程为988 m，坝底高程为723 m，最大坝高为265 m。坝址处河谷形态呈狭窄的“V”型，两岸地形对称，高陡险峻，坡形完整。拱座岩体为落雪组第3段第1～4亚段（P t32-11～P t32-14）变质灰岩、大理岩，以厚层为主，夹中厚层、巨厚层及少量薄层；地层产状走向为260°～280°，倾向SW，倾角为70°～85°，岩层走向与河流流向交角为 60°～80°，为斜横向谷。拱座岩体中没有斜切边坡、贯穿性强的断层和裂隙性断层发育。

根据可研阶段资料，大坝拱肩槽呈弧形开挖，左岸开挖边坡走向从62°转至100°，开挖走向与岩层走向接近平行，拱肩槽正面坡为横向坡；右岸开挖边坡走向从68°转至 76°，拱肩槽正面坡为斜交顺向坡。左岸开挖边坡坡顶高程约为1163 m，最大坡高约为440 m；右岸开挖边坡坡顶高程约为1120 m，最大坡高约为397 m。

左坝肩岩体主要为 Pt2l3-1厚层灰岩夹大理岩，局部有少量（高程为810 m以下）Pt2l2-3互层大理岩化白云岩。右坝肩中下部（高程为875～905 m以下）为Pt2l3-1厚层灰岩夹大理岩；上部（高程875～905 m以上）为Pt2l3-2-1～Pt2l2-2-6巨厚层白云岩夹中厚层灰岩、厚层石英岩。两岸坝肩建基岩体以Ⅱ1级为主，岩体质量优良，左坝肩主要为Ⅱ1(95.2%)，少量为Ⅲ1(4.8%)，右坝肩主要为Ⅱ1(82.7%)，少量为Ⅲ1(17.3%)。

3 拱坝联合坝址区岩体三维模型

3.1 拱坝坝址区自然地质条件三维模型

拱坝坝址区自然地质条件的三维建模主要包括计算范围的选择、地形地貌、关键地质结构（地层、断层、溶洞系统、软弱夹层等）。

三维计算模型以坝体为中心，沿顺河向、横河向及坝底深部基岩以下一定的范围内生成。下游模拟范围及左右两岸模拟范围均大于2.5倍坝高。整体模拟范围为1225 m×1225 m×1440 m。

坝址区会理群（Pt2H）地层主要出露因民组（Pt2y）和落雪组（Pt2l）两组地层，岩性有8类。根据地质勘探数据及二维地质剖面信息，将设计单位提供的 CAD格式的工程地质图、地形图、各高程平切图等进行预处理，剔除非几何边界的冗余信息，并直接输出为 DXF矢量数据格式，包含地表及岩性分界线的DXF数据文件直接输入Hypermesh软件，对地表、地层及断层进行插值和拟合处理，坝址区岩层几何模型如图1所示。

坝址区地表实测断层有35条，规模较大。对坝肩稳定性有影响的断层有4条，从上游至下游分别为红沟断层F3、梁子断层F7、白沟断层F14、雷家湾沟断层F15。F15走向NNE，与河谷斜交。其余断层与河谷大角度相交，坝线选择段地表断层不发育，可见9条小断层，模型考虑距离坝体较近的次一级断层 f40和f42。右岸存在已停止发育的充填古岩溶系统(K25)，岩溶系统发育于Pt2l3-3厚层灰岩与Pt2l3-2厚层白云岩分界面附近 Pt2l3-3厚层灰岩中，基本上顺层发育，高程向下具向坡外发育和规模逐渐减小的趋势，平切面上呈椭圆状，长轴与岩层走向方面近于一致，短轴与岩层倾向近于一致。断层和岩溶与坝体相对位置如图3所示。

图1 岩层三维地质模型Fig.1 3D geological model for rock mass

图2 断层和岩溶与坝体的位置Fig.2 Positions of dam,faults and karst

图3 整体有限元网格模型Fig.3 Integral finite element model

借助Hypermesh软件在几何体的布尔（Boolean）操作运算方面的优越性，模型三维实体模型采用图2结构面切割图1岩层几何体完成。异形几何体的空间复杂关系使得有限元网格的划分工作十分繁琐。利用HyperMesh的网格划分模块交互调整特性，方便调整曲面或边界的网格参数，特别是对单个单元及节点的调整功能，提高不规则几何体网格划分效率。

本文岩体全部采用8节点六面体和6节点五面体实体单元，F3、F14、F15及 K25采用实体单元模拟。小断层f40和f42采用无厚度节理单元模拟，计算网格如图 3所示。有限元网格划分后单元总数为294992，节点总数为 300272（坝体单元总数为10040）。根据岩体力学参数建议值，地层岩性计算力学参数见表1。

表1 地层岩性计算力学参数Table 1 Mechanical parameters of materials for calculation

3.2 坝址区初始地应力反演分析三维模型

假定现今初始地应力场主要是远古地应力场在后期的地表剥蚀与河流侵蚀下切等外动力地质作用下逐渐形成。根据坝址区实测地应力资料，假定应力边界条件，反演最优地表剥蚀和河流侵蚀分层厚度[6－7]。计算需模拟从远古模拟夷平面地表剥蚀至现今坝址区河谷地貌，反演构造应力边界条件。将经过神经网络反演得到的边界条件和重力加速度代入进行正算，从而获取坝址区初始地应力场。

本文坝址区夷平面的高程选取1800 m高程作为地表剥蚀和河流侵蚀下切的起点，以坝址区现今河谷高程最低的点（高程值为723 m）为地表侵蚀下切的终点，河谷的下切过程分为6次完成。现今地表与远古夷平面间平均插值5个侵蚀平面，利用Hypermesh软件网格投影功能，地表面网格分别在5个侵蚀面及远古夷平面进行二维投影。侵蚀面间实体网格形态由侵蚀面二维预置网格控制，由此方法可快速构建地应力反演分析数值模型，如图4所示。

图4 初始地应力反演分析网格模型Fig.4 Finite element model for initial geostress field inversion

3.3 坝肩槽开挖分析三维模型

模型高边坡包括坝肩坝顶高程以上边坡及拱肩槽边坡。左岸高边坡开口线最高高程为1155 m，平均开挖坡度为75°；右岸高边坡开口线最高高程为1135 m，平均开挖坡度为65°。拱肩槽边坡平面形态似弧形，略凸向上游侧。根据上、下游边坡不同形态，可进一步分为拱肩槽上游边坡、正面边坡及拱肩槽下游边坡。左右岸开挖边坡大致呈对称分布，两岸一般每隔20 m设置一马道。开挖后地形地貌基本顺原坡体走向，退坡深度不大。

根据坝肩边坡设计及坝体体型方案，在 CAD软件中按照开挖线高程及坝体横剖面高程，构建人工对象边界模型。导入到Hypermesh软件，快速生成曲面，与坝址区三维地质模型进行体与面之间的布尔运算，将开挖岩体部分切割分离。

图5 拱肩开挖边坡网格模型Fig.5 Finite element model for abutment slope excavation

3.4 坝肩岩体滑动楔块三维滑面模型

地勘裂隙统计资料表明，拱座内无定位块体，控制拱座滑动稳定的主要是短小裂隙形成的非连通的滑动边界。裂隙可分为SN组和EW组，其中SN组与河流平行，两岸均发育，占总数的54.1%，EW组占总数的 15.4%，与河流近垂直，两岸均发育。裂隙以中倾角为主，次为陡倾角，缓倾角裂隙不发育，中倾角裂隙以SN组裂隙为主，优势产状为倾向273°、倾角 51°[8]。根据裂隙产状和规模及岩层分布情况，确定9个可能滑动楔块，典型楔块滑面如图6所示。

图6 典型楔块滑面Fig.6 Sliding surface of typical wedge block

本文研究的问题中存在多个可能滑动楔块的情况下，若三维建模时考虑所有滑面情况，则使得拱座内岩体支离破碎，单元网格质量差，操作繁琐，不利于推广。

本文采用的是采用滑面多重网格的方法，即将有限元计算的应力结果插值到重建的滑面网格上，从而适应复杂滑面的组合[9]。在Hypermesh软件中确定楔块滑面后，利用卸荷带界线切割楔块滑面并进行网格划分。拱座内滑动楔块由新鲜岩体、弱卸荷岩体及强卸荷岩体组成，如图7所示。

图7 楔块岩体卸荷带界线Fig.7 Classification of unloading zones of wedge blocks

4 计算结果

将三维数值模型通过接口程序进行单元数据格式转换[10]，模型导入到FLAC3D程序中进行三维弹塑性数值计算。根据地应力反演结果，坝址区地应力分布可分为3个区域，沿坝轴线剖面初始地应力s1分布如图8所示。

图8 坝轴线剖面地应力分布Fig.8 Stress distribution at dam axis section

两侧河谷近地表地应力带：受地表剥蚀和金沙江长期冲蚀下切作用，最大水平主应力量值一般为10 MPa左右，以小于10 MPa的低地应力为主；总体上，左岸较右岸稍大，岸坡地形对地应力有一定影响。河床深部地应力区：两岸深部高程在 750～850 m范围内，应力量值属中等地应力水平；河床深部最大主应力达20 MPa左右，局部达到高应力水平。构造应力区：在断层区域，受断层挤压剪切影响，断层两侧地应力突变，以左岸雷家沟断层F15最为明显。因此，断层穿越地段受断层控制，为应力集中的构造应力带。

计算模拟过程按照初始地应力反演、高边坡及坝肩槽开挖、坝体施工、坝体载荷的步骤进行。考虑荷载不同组合，计算模拟分为4个方案。方案1荷载组合为正常蓄水+自重+泥沙压力+浪压力+冰压力。方案2荷载组合为方案1荷载组合叠加温升。方案3荷载组合为方案1荷载组合叠加温降。方案4荷载组合为方案2荷载组合叠加地震荷载。坝肩槽开挖完成后，坝肩附近岩体顺河向位移如图9所示，右岸拱肩槽 800 m高程附近位移值最大达到12.9 mm。左右岸可能滑动楔块安全系数计算如表2所示。

正常蓄水后（方案 1），各楔块具有较高的安全裕度，整体上左岸稳定性强于右岸，安全系数计算结果如表2所示。右岸滑面通过岩溶区域的楔块相对稳定性较差，楔块 RM1安全系数最小值为5.399。基本组合2条件下，温降时拱圈缩短并偏向右岸下游方向变形，左岸楔块相对方案 1，楔块稳定性有所增加。右岸部分楔块稳定性降低，最大降低幅度7.4%。地震荷载作用下所有楔块安全系数均大幅下降，左岸楔块LM2-1安全系数最小为5.511。

图9 拱肩开挖边坡位移分布(沿河向)Fig.9 Displacement distribution of excavation slope

表2 楔块安全系数Table 2 Safety factors of blocks

5 结 论

采用通用软件平台，建立用于高拱坝坝肩稳定分析所需的三维数值模型的方法，并以乌东德工程为例，构建适用于坝肩稳定性分析的三维精细数值模型。利用地勘资料和工程设计数据，分析了影响坝肩稳定的地质特征和工程措施，采用面向对象技术、多种 CAD软件协同作业，建立自然地质条件和工程对象的几何模型。以坝肩部位为中心，基于CAE领域功能强大的Hypermesh软件，较好地完成了复杂地质结构与工程对象的网格划分。通过接口程序进行三维数值模型单元格式转换，在 FLAC3D完成三维弹塑性计算。

本文将地应力反演分析、坝体载荷、工程开挖等多因素影响下坝肩静动力稳定性分析的问题统一到独立模型中。计算结果表明，精细的三维模型最大程度地模拟了工程岩体复杂结构，为高拱坝坝肩岩体稳定性分析提供可靠的基础。

[1]李明超,钟登华,秦朝霞,等.基于三维地质模型的工程岩体结构精细数值建模[J].岩石力学与工程学报,2007,26(9): 1893－1898.LI Ming-chao,ZHONG Deng-hua,QIN Zhao-xia,et al.Refined modeling for numerical simulation of engineering rock mass structures based on 3D geological model[J].Chinese Journal of Rock Mechanics and Engineering,2007,26(9): 1893－1898.

[2]杨宝全,张林,陈建叶,等.复杂地质条件下拱坝坝肩稳定分析及密集节理影响研究[J].岩石力学与工程学报,2010,29(增刊2): 3972－3978.YANG Bao-quan,ZHANG Lin,CHEN Jian-ye,et al.Analysis of arch dam’s abutment stability under complicated geological conditions and study of influence of dense joints[J].Chinese Journal of Rock Mechanics and Engineering,2010,29(Supp.2): 3972－3978.

[3]刘杰，李建林，宛良朋，等.基于卸荷和抗震分析理论的大岗山坝肩边坡锚固优化研究[J].岩土力学,2012,33(增刊 2): 275－282.LIU Jie,LI Jian-lin,WAN Liang-peng,et al.Research on Dagangshan dam abutment slope anchorage optimization based on theory of unloading and seismic analysis[J].Rock and Soil Mechanics,2012,33(Supp.2): 275－282.

[4]漆祖芳,姜清辉,唐志丹,等.锦屏一级水电站左岸坝肩边坡施工期稳定分析[J].岩土力学,2012,33(2): 531－538.QI Zu-fang,JIANG Qing-hui,TANG Zhi-dan,et al.Stability analysis of abutment slope at left bank of Jinping-I hydropower project during construction[J].Rock and Soil Mechanics,2012,33(2): 531－538.

[5]侯恩科,吴立新,李建民,等.三维地学模拟与数值模拟的耦合方法研究[J].煤炭学报,2002,27(4): 388－392.HOU En-ke,WU Li-xin,LI Jian-min,et al.Study on the coupling of 3D geoscience modeling with numerical simulation[J].Journal of China Coal Society,2002,27(4): 388－392.

[6]张传庆,周辉,冯夏庭,等.局域地应力场获取的插值平衡方法[J].岩土力学,2008,29(8): 2016－2024.ZHANG Chuan-qing,ZHOU Hui,FENG Xia-ting,et al.Method for obtaining local stress field based on interpolation and numerical simulation method[J].Rock and Soil Mechanics,2008,29(8): 2016－2024.

[7]张振华.深切河谷岸坡开挖过程动态预警方法研究[博士学位论文 D].武汉: 中国科学院武汉岩土力学研究所,2008.

[8]沈辉,罗先启,李野,等.乌东德拱坝坝肩三维抗滑稳定分析[J].岩石力学与工程学报,2012,31(5): 1026－1033.SHEN Hui,LUO Xian-qi,LI Ye,et al.Three-dimensional stability analysis of dam abutment of Wudongde arch dam[J].Chinese Journal of Rock Mechanics and Engineering,2012,31(5): 1026－1033.

[9]杨强,朱玲,翟明杰.基于三维非线性有限元的坝肩稳定刚体极限平衡法机理研究[J].岩石力学与工程学报,2005,24(19): 3403－3409.YANG Qiang,ZHU Ling,ZHAI Ming-jie.On mechanism of rigid-body limit equilibrium method for abutment stability analysis of arch dam[J].Chinese Journal of Rock Mechanics and Engineering,2005,24(19): 3403－3409.

[10]刘心庭,唐辉明.FLAC3D复杂网格模型的构建及其工程应用[J].金属矿山,2010,(11): 108－111.LIU Xin-ting,TANG Hui-ming.FLAC3Dcomplex grid model construction and its engineering application[J].Metal Mine,2010,(11): 108－111.