库水压缩性及分缝布置对拱坝非线性地震响应的影响

陈灯红 赵艺园 杨紫辉 洪海丰 曹文昱

(1.三峡大学 土木与建筑学院, 湖北 宜昌 443002;2.三峡大学 水利与环境学院, 湖北 宜昌 443002)

中国的西南、西北地区是全国水能资源最为丰富的地区,也是已建、在建、拟建高混凝土坝最多的区域.地震科学数据共享中心[1]的数据显示,过去10年这一区域发生5级以上地震次数占中国境内5级以上地震总数的60.75%,占中国大陆5级以上地震总数的91.12%.频发的地震会对区域内的水库大坝造成不同程度的损伤,因此对混凝土大坝进行抗震性能分析,研究其在地震作用下的动态响应,对于安全、合理、充分地利用丰富的水能资源具有实际的工程意义[2-3].

目前研究混凝土大坝在地震过程中动态响应的方法主要有:通过实际观测和记录大坝地震反应、室内振动台模型试验、数值计算及理论分析.由于模型试验存在一定的局限性,因此近年来在国内外学者的研究和推广下,数值模拟计算在应用于混凝土大坝动力响应分析方面取得了长足进步.杜修力等[4]建立了拱坝-可压缩库水-地基模型进行拱坝地震响应分析.邱奕翔等[5]以拉西瓦拱坝为例,建立三维有限元模型,采用附加质量与流固耦合法,研究库水模拟对拱坝动力特性的影响.杨柳等[6]采用ADINA 研究考虑库水可压缩性对拱坝动力响应的影响.龙渝川等[7]采用限制接触面切向滑移的弹簧单元,模拟横缝在地震动荷载作用下的非线性力学行为,并利用有限元软件的隐士迭代算法对高拱坝进行动力响应的求解,并将计算结果与清华大学学者扩展的ADAP-88 程序计算结果进行比较,验证了该方法的合理性.张宇等[8]以沙牌碾压混凝土拱坝为例,通过物理模拟与已观测的实际震害对比,研究其考虑结构缝的地震破坏机理和失效模式,为类似的工程设计提供了一定的参考依据.

综合上述学者的研究结果,可压缩库水能真实反映库水的特性及库水对坝体动力特性的影响;数值模拟计算与试验验证均表明设置横缝是必要的.但对于这两类问题的研究大都是独立的,这样建模时考虑的层面较为单一,对此国内外学者已经发展了多个拱坝-库水-地基非线性动力分析模型,为研究拱坝地震损伤特性方面提供了一定的参考依据.Chopra等[9-10]通过考虑库水的压缩性与坝基岩石相互作用,建立三维非线性大坝分析模型对大坝的地震响应进行研究.张楚汉等[11]分析总结了高混凝土坝抗震安全评价方面最新的研究成果与高拱坝抗震研究的发展方向,并提出了一种综合的大坝抗震安全评价体系研究思路.陈健云等[12]通过考虑横缝的高拱坝在不同强度地震荷载作用下响应的分析,提出了拱坝抗震性能的3个水平阶段,为解决超静定拱坝在强震下的整体抗震性能评价指标方面提供了参考.范书立等[13]采用考虑键槽咬合作用、横缝设置、Westergaard 附加质量的非线性动力分析模型,对白鹤滩拱坝进行地震易损性分析.Kadkhodayan 等[14]通过考虑坝体收缩缝、库水-地基和库水-坝体的相互作用,建立三维有限元模型,并基于IDA 方法对拱坝进行抗震性能研究.

在上述研究的基础上,构建了一种较为综合全面的高拱坝-库水-地基动力相互作用三维模型,该模型通过采用声学单元模拟库水,以考虑库水与坝体及地基的耦合;采用动接触-联结混合模型中能够实现接触面间法向光滑过渡的指数型动接触力模型模拟横缝接触非线性;在混凝土受拉(压)本构模型中引入动力放大系数实现混凝土损伤塑性模型的率型化;并且,对锦屏Ⅰ级高拱坝-库水-地基系统进行了地震响应分析,研究了库水压缩性及分缝布置对该拱坝非线性地震响应的影响.

1 计算模型

1.1 库水与地基、坝体的流固耦合

考虑到库水黏性很小的特性,将其假定为均质、无黏、无旋、小扰动的理想流体.设流体中质点(x,y,z)的运动位移为(u,v,w),由牛顿定律可得:

若不考虑库水的可压缩性,可令式(3)右侧等于0[5].应用Galerkin法,并乘以声压的变分δp,在流体区域V内积分[13],经过计算得

式中:u为S面上位移向量;L=∇().

将流体方程离散化,分成若干个有限单元,单元内任意一点声压和质点的位移及其对时间的各阶导数均可由该单元节点上相应值插值表示,并将声压变分约去,得到完全耦合的结构流体运动方程为:

式中:MS、CS、KS分别为结构质量矩阵、结构阻尼阵和刚度矩阵;Mf、Cf、Kf分别为流体质量矩阵、声阻尼矩阵和流体刚度矩阵;R为流体和结构的耦合矩阵;U、P为节点位移向量和声压向量;FS为结构载荷向量.

基于以上的运动方程,并结合本文研究内容设置如图1所示库水模型的边界条件.

图1 库水模型边界条件

式中:n为截断库水边界内法向;θ为与内法向之间的入射角[15].

1.2 动接触边界及键槽咬合作用的模拟

采用单一力学模型不足以描述横缝接触中的非线性.例如,仅使用联结单元模型,如何选取既保证精度又保证收敛性的单元刚度系数,依赖于研究人员的经验;仅使用动接触力模型,则会忽略键槽的咬合作用.因此,以动接触-联结混合模型模拟考虑键槽咬合作用的横缝接触非线性.非线性接触横缝的相互作用分为法向和切向,法向采用动接触边界中的指数型模型;切向采用库仑摩擦模型描述混凝土接触面的摩擦,键槽的咬合作用则采用附加切向刚度弹簧的点对联结单元模型,以此构建了一种带附加切向刚度的平缝模型.该模型与实际的键槽横缝相比,不仅能改善由于引入联结单元模型导致的模型计算收敛性较差的问题,也不要求接触两侧的网格完全匹配,可在一定程度上提高求解效率,也使建模过程更为便利.

考虑了如图2所示的接触边界,该接触边界由主面和从面构成,从面上的节点S′在主面上有且仅有一个确定的节点M′(锚点)与之对应.锚点与对应结点的连线方向即为接触的法线方向;从面结点与锚点的法向和切向相对距离,分别表征动接触的横缝开度与切向位移.通过从面结点与主面锚点之间的相对位移可以确定变形过程中横缝的接触状态[7].

图2 接触边界示意图

从面结点S′到主面区间(M-1,M)的法向和切向距离可由式(11)确定.当h＞0时,主从面接触;当h＜0时,主从面分离;当h=0时,主从面接触且相互之间没有力的作用.

式中:c为初始间隙(取c=5×10-7m);h为接触面之间的相对位移(以嵌入为正);p为接触点对上的接触压力;p0为特征接触力(取p0=5×109Pa).

考虑横缝两端坝体接触是键槽的咬合作用,相邻坝段间的顺河向滑动受到约束.因此,采用联结单元模型中的点对弹簧单元,即附加切向弹簧的方式来模拟键槽的咬合作用(如图3所示).沿缝面布置切向刚度为KS的线性弹簧单元,通过调整KS可以达到控制切向位移的目的.

1.3 率相关本构方程构建及损伤因子的转化

1.3.1 率相关本构方程的构建

高烈度地震区修建的混凝土高坝,其混凝土材料在高应变率下体现的动力特性,会直接影响到坝体的地震响应.为了更合理地对混凝土坝进行抗震性能分析,计算模型应考虑混凝土材料的率相关性.

Lubliner等[16-17]提出的混凝土塑性损伤模型(CDP)可以较好地模拟动静荷载作用下混凝土的损伤演化过程,在混凝土坝地震响应分析中得到较多应用.在此,将选取的混凝土受拉(压)本构模型嵌入到该CDP模型中,并引入动力放大系数(即在静态的等效单轴拉压本构关系中考虑应变率),完成三维动态混凝土率相关本构模型的构建.

混凝土单轴受压状态下的本构方程采用王乾峰[17]修正后的单轴受压本构模型来定义单轴受压曲线,并对其试验数据进行拟合,得出受压状态下的动力放大系数(CDIF),二者结合确定混凝土单轴受压率相关本构方程为:

1.3.2 损伤因子的转化

本文构建的混凝土单轴受压(拉)率相关本构方程是基于弹性损伤模型建立的,若将损伤因子dc、dt直接嵌入到混凝土塑性损伤(CDP)模型中,会出现由于卸载路径交叉而导致计算不收敛的情况,因此需要将损伤因子dc、dt转化成CDP模型中统一的格式.

Sidoroff根据能量等效原理所提出的混凝土损伤因子计算方法,可直接嵌入到CDP 模型中进行计算,且计算结果能较好的符合实际情况[21].

无损材料的弹性余能为:

将式(15)、(18)和式(25)进行对比分析,可得出如下关系式:

由式(26)进一步可推出两种损伤因子之间的换算公式为:

将式(16)、式(19)得出的受拉(压)损伤因子经由式(27),转换为塑性损伤因子后,可用于CDP模型中计算.

2 工程算例

选取锦屏一级混凝土双曲拱坝为工程算例,其坝顶高程1 885 m,建基面高程1 580 m,最大坝高305 m.电站正常蓄水位1880 m,死水位1800 m,拱冠梁顶厚16 m,拱冠梁底厚63 m,最大中心角93.12°,顶拱中心线弧长552.23 m,厚高比0.207,弧高比1.811.

2.1 模型参数

将坝体按照高程从下至上分为3个区域(如图4所示).

图4 坝体混凝土分区示意图

坝体混凝土材料的基本属性为:密度2500 kg/m3;泊松比0.167;A 区混凝土弹性模量29.8 GPa、B 区混凝土弹性模量27.9 GPa、C 区混凝土弹性模量25.4 GPa.各分区混凝土材料的动态应力控制标准按照《水工建筑物抗震设计标准》(GB 51247—2018)取值[22],动态抗拉强度的标准值可取为其动态抗压强度标准值的10%.模型参数采用CDP模型参数(见表1),据此来模拟坝体混凝土的动态力学性能.为突出研究重点,地基按传统的无质量地基模型考虑[12,23],其静弹性模量为30 GPa、泊松比为0.2.

表1 CDP模型参数取值

2.2 有限元模型

为研究库水压缩性及坝体横缝布置对大坝抗震性能的影响,建立以声学单元模拟库水的三维坝体-库水-地基模型(如图5所示).结合本文研究目的设立了3种计算工况:工况1:设有13条横缝的有限元模型、库水可压缩;工况2:设有13条横缝的有限元模型、库水不可压缩;工况3:设有0条横缝的有限元模型、库水可压缩.其中不设横缝模型共有节点11 324个,单元8 935个;设置13条横缝模型共有节点12 158个,单元8935个.模型中地基单元3 900个;坝体单元2 205个;库水单元2 830个.

图5 三维坝体-库水-地基有限元模型

模型的计算静荷载为坝体自重与静水压力;动荷载为动水压力与地震动荷载,基线调零后人工拟合的地震动加速度时程曲线如图6所示,其最大可信地震的水平、竖直峰值加速度代表值分别为0.323g、0.215g.

图6 基线调零后人工拟合地震动加速度时程曲线

3 库水压缩性对拱坝非线性地震响应的影响

3.1 库水压缩性对比模型的选取

流固耦合模型中库水压缩性,可通过控制声学介质的体积模量来实现.本节取设有13条横缝的三维坝体-库水-地基模型作为研究对象,通过改变库水单元的体积模量,达到对库水压缩性的模拟.其中,库水可压缩工况(工况1)体积模量为2.07 GPa,库水不可压缩工况(工况2)体积模量为2.07×1020Pa.

3.2 库水压缩性对于拱坝结构模态参数的影响

在此需要说明的是,使用的动接触-联结混合模型中,从面结点仅在与主面接触时产生切向约束力,当横缝张开时,主面与从面之间就没有任何切向约束即切向摩擦力为零,只有附加切向力发生作用[5].从而在计算自振频率时切向非线性弹簧会退化成线性弹簧,则可以通过计算自振频率并使用瑞利阻尼考虑结构阻尼对拱坝动力响应的影响.阻尼系数α=2ξω1ω2/(ω1+ω2),β=2ξ/(ω1+ω2),其中ξ为阻尼比,ω1和ω2为选取的模态频率上下限.

据表2所示结果,对三维坝体-库水-地基有限元模型工况1和工况2计算后的自振频率及阵型进行分析与整理,并对比工况1 和工况2 前8 阶振型得出:拱坝前4阶振型的对称形式相同,但第5、7、8阶振型两种工况的对称形式则有所不同;除第3阶工况2较工况1自振频率的减少量达到34.33%外,其他自振频率减少量变动幅度相对较小.

表2 两种工况下拱坝的前8阶自振频率

3.3 结果分析

3.3.1 应力分析

表3给出了两种工况下坝体主要位置的拉压应力峰值.据表中的数据得出,可压缩工况计算得到的横河向与竖直向的拉应力与压应力值均小于不可压缩工况,其中拉应力峰值相差较大的位置位于竖直向坝顶拱冠处,可达43.33%.这表明不可压缩工况过大地估计了拱坝抗震主要部位的应力.顺河向结果恰好相反,这是因为可压缩库水基频与拱坝坝体正对称阵型所对应的基频接近,库水高阶频率对坝体的反应远不如不可压缩工况大,从而导致考虑库水可压缩性时所得的拱坝坝体顺河向应力值较不可压缩工况大.

表3 不同工况下坝体主要位置应力峰值

3.3.2 损伤分析

如图7所示,两种工况下坝体受拉损伤因子分布区域大致相同,坝体与地基接触面的受损情况要远大于坝体其他部位.将横缝按上游坝面左岸至右岸依次编为1～13号,其中,工况1拉损极值(0.933)出现在左岸坝肩1号横缝缝端;工况2拉损极值(0.905)出现在拱冠梁附近的9号横缝缝端.如图8所示,两种工况拱冠梁的拉损极值均出现在坝踵区域,工况2的拉损(0.797)情况比工况1(0.526)要严重一些.出现这一现象的原因是工况2库水体积模量的增加,使得库水单元的刚度变大,在地震动荷载作用下能够产生更大且垂直于坝面的动水压力.

图7 坝体受拉损伤分布

图8 拱冠梁受拉损伤分布

3.3.3 位移分析

图9给出了坝体上游面顺河向位移等值线图.从图中可以得出两种工况下:坝体位移由坝肩向拱冠梁逐渐增大,由于库水刚度不同,导致坝面动水压力的差异,使得相同区域内工况1 的位移部分大于工况2.在最大可信地震作用下拱冠梁的残余相对位移,工况1为8.00 cm,工况2为7.00 cm.说明了地震动荷载作用下工况2拱冠梁运动幅度大,但其地震动荷载作用后的永久位移要小于工况1.

图9 坝体上游面顺河向位移云图(cm)

3.3.4 横缝开度

如图10所示为两种工况下横缝开度极值包络图,可以得出在地震动作用下的横缝张开情况:基本呈左右对称分布,横缝开度较大的区域分布在左、右岸坝肩以及拱冠梁附近.其中工况1横缝最大开度为15.41 mm(12 号横缝),工况2 横缝最大开度为11.96 mm(1号横缝),工况1横缝最大开度较工况2增加了28.85%.

图10 横缝开度极值包络图

库水体积模量会影响横缝开度的最值,库水体积模量的增加会使得横缝开度最值的降低,横缝开度最值会出现在左、右岸坝肩区域内,工程实际中对于位于该区域应该予以重视.

3.3.5 动水压力分析

如图11给出了两个不同工况拱冠梁的动水压力分布图和包络曲线图.从图中对比得出:等高程处库水可压缩工况拱冠梁的动水压力要小于不可压缩工况的动水压力.坝踵处动水压力时程曲线如图12所示,在相同位置,可压缩工况动水压力最大值为0.619 MPa,不可压缩工况动水压力最大值为0.913 MPa,前者比后者小约32.20%.通过纵向与横向对比得出不考虑拱坝库水的压缩性会夸大动水压力作用.

图11 拱冠梁剖面动水压力分布图与包络曲线图

图12 库水压缩性对比工况的坝踵处动水压力时程曲线

4 分缝布置对于拱坝非线性地震响应的影响

4.1 分缝布置对比模型的选取

本节选取设有13条横缝(工况1)、设有0条横缝(工况3)的三维坝体-库水-地基模型作为研究对象(库水可压缩,体积模量为2.07 GPa),研究设置横缝对坝体动力响应的影响.

4.2 结果分析

4.2.1 应力分析

经运算得到了不同工况下坝体主要位置各方向拉压应力峰值,见表4.对表中数据对比得出,两个工况的结果存在一定的差异,导致这一现象产生的原因是横缝在地震动作用下会出现张开闭合现象,破坏了结构的整体性并降低了拱坝的刚度,引起应力重分布进而影响拱坝的抗震安全.

表4 不同工况下坝体主要位置应力峰值(单位:MPa)

4.2.2 位移分析

工况1、3拱冠梁顺河向相对位移时程曲线如图13所示,正(负)位移值表示拱冠梁沿y轴向下(上)游运动.工况1曲线峰值为0.262 m、-0.114 m;工况3曲线峰值为0.317 m、-0.391 m.表明在拱坝有限元分析中,若不考虑横缝的影响,拱冠梁在地震动作用下将发生剧烈运动,且工况3的位移曲线在12.54 s后出现了明显的漂移,而位移曲线的漂移常常作为拱冠梁失效的指标之一.即在拱坝的有限元分析中,忽略横缝的作用和影响,会使计算结果失真.

图13 工况1、3拱冠梁顺河向相对位移时程曲线

4.2.3 损伤分析

坝体上游面、坝基的受拉损伤云图如图14、15所示.可以得出对于拱坝而言,考虑横缝与否,不会影响损伤区域集中在左右坝肩、拱冠梁顶以及坝踵区域.工况1的损伤状况明显好于工况3,说明了横缝的存在使得拱向应力向梁向应力转化,改善了坝肩的受力状况.如图14中工况3在左岸坝肩处出现了贯穿坝体的损伤区域,工况1 在此区域内虽有损伤超过80%,但损伤区域扩展长度小于坝厚的1/3.表明横缝的存在改善了由于“角缘效应”导致的左岸坝肩损伤集中现象.

图14 坝体上游面受拉损伤云图

图15 坝基受拉损伤云图

5 结 论

本文研究了库水可压缩性及分缝布置对高拱坝非线性地震响应的影响,得出以下结论:

1)在横缝模型中考虑库水压缩性,可压缩工况的横河向与竖直向的拉应力与压应力值均小于不可压缩工况,其中拉应力峰值(可压缩工况0.68 MPa,不可压缩工况1.20 MPa)相差较大的位置位于竖直向坝顶拱冠处,达到43.33%.表明库水不可压缩工况过大地估计了拱坝主要部位的应力.

2)库水压缩性对横缝开度的峰值也有一定影响,不仅影响横缝开度峰值大小也影响其产生的位置,其中可压缩工况横缝开度最值较不可压缩工况相差28.85%,且压缩工况横缝开度峰值出现在左岸坝肩,而不可压工况峰值出现在右岸坝肩.

3)在横缝模型中考虑库水压缩性,通过对比两个工况的动水压力值发现,同一位置(坝踵处)可压缩工况动水压力最大值(0.619 MPa)较不可压缩工况的最大值(0.913 MPa)小32.20%;拱冠梁等高程处库水可压缩工况的动水压力值均小于不可压缩工况,即不考虑拱坝库水的压缩性会夸大动水压力作用.

4)在考虑库水可压缩性时,不设横缝模型的各向拉应力虽未超过各区域动态应力控制标准,但其在左岸坝肩处的损伤区域已经贯穿坝体,且拱冠梁的相对位移也出现了明显的漂移;设置13条横缝模型的损伤区域分布及拱冠梁的相对位移,均未出现拱坝运行状态不良的迹象,说明了在拱坝的地震响应分析中,考虑横缝的影响,能够对拱坝的地震响应进行更好的描述.