SV波斜入射下Nam Ngum 5拱坝-库水-地基系统的地震响应分析

刘云辉,陈灯红,潘子悦,赵艺园,刘云龙,陈非凡

(1. 三峡大学 防灾减灾湖北省重点实验室,湖北 宜昌 443002; 2. 三峡大学 土木与建筑学院,湖北 宜昌 443002)

0 引言

由于我国处于两大地震带之间的特殊地理位置,地质构造十分复杂且勘探条件有限,在遭遇强度大、震源浅、分布广的地震时往往引发十分严重的地质灾害。特别是许多大型混凝土重力坝坐落于与高烈度区重合的西部地区,因此,对于大坝抗震安全研究是在水利建设中必不可少的一环。其中利用数值仿真技术来研究大坝抗震问题是非常有效的模拟手段之一。

对于坝体的抗震安全性评价主要有4个重要的环节:坝址地震动输入、高坝地震动反应分析、高坝抗震动力安全评价以及坝体和基岩材料动力特性与破坏机理研究[1]。其中,坝体-地基-库水系统[2-6]地震反应是国内外高坝抗震研究的主要内容。合理的利用地震动输入形式并建立相应的输入方法是大坝抗震分析的重要内容[7]。 JOYNER等[8]考虑了弹性基底中的有限刚度,提出等效节点荷载的理念,很好地解决黏弹性边界上地震波的输入。刘晶波等[9-10]将波动输入方法从一维推广到二维并给出了一种结构-地基相互作用中地震波动输入方法。杜修力等[11-13]在拱坝-地基系统的地震响应分析中,为获得带有人工边界的有限离散动态模型的非耦合方法,将显式有限元的波动求解方法与黏弹性边界相结合,通过设置一组人工边界来模拟出散射波从内部区域到无限外部区域的传播,从而模拟出无限外部区域。

对于实际的工程而言,地震波的入射形式往往不是如假定的垂直入射,而是以一定的角度输入到相应的结构体中。在对于混凝土坝等大体积水工建筑物的地震波斜入射研究中,吴兆营[14]以土石坝为依据探究了坝体在地震波斜入射时的最不利输入。孙奔博等[15]将隐式有限元的求解方式应用于应力型黏弹性人工边界中,将斜入射问题扩展到三维分析模型中。杜修力等[16]以小湾拱坝为研究对象探究了三维情况下SV波斜入射对拱坝的地震反应的影响。李明超等[17]在利用黏弹性人工边界分析斜入射下混凝土重力坝的塑性损伤分析时,发现地基边界材料的复杂性会影响黏弹性边界处等效节点荷载的确定。黄景琦等[18]在文献中指出大型结构计算单元与节点数目过多时,通过相关的编程软件批量处理等效节点荷载导入至有限元计算软件中时,每一个节点施加的等效节点应力不易于实现。SU等[19]提出了一种在ABAQUS中结构整体静动力分析中稳定切换边界条件类型的方案,通过对大坝在静、动荷载作用下的内部和外部波动问题进行了二维有限元分析,验证了该方法能稳定切换边界条件,且提高计算精度。通过以上相关学者的对比研究中发现传统的黏弹性人工边界方法[20]应用于大型结构的仿真分析时,由于大量的等效节点力需要通过三方编程软件处理,使得波动输入前处理工作量大,也比较复杂,不易于波动输入方法的实现及工程应用。且在结构的动力响应分析中,因静、动力分析时边界条件不同,若同时考虑会使得在进行有限元模型计算时,由于静力分析与动力分析人工边界转换过程中会产生一定的误差。尤其是当这种方法应用于地震波斜入射的相关研究问题中,如果利用常规的分析方法进行数值模拟时会使得计算误差进一步增大,因此需对斜入射问题进行更加深入的研究。

本文在以上学者对地震波输入方法研究的基础上提出一种三维黏弹性静-动力统一人工边界对地震波斜入射进行研究。并以Nam Ngum 5重力拱坝为研究对象,建立三维拱坝-库水-地基系统的斜入射模型,推导了任意角度下三维SV波斜入射时的等效节点荷载,构建了三维SV波斜入射下三维黏弹性静-动力统一人工边界,并通过用户子程序接口编写了三维SV波的波动输入子程序,基于大型软件ABAQUS完成地震动输入,探究了三维坝体在SV波斜入射作用下的动力响应。

1 三维SV波斜入射模型构建

1.1 黏弹性人工边界

本文采用刘晶波等[21-22]构建的黏弹性边界,即在边界处通过在有限元模型中边界节点处建立一系列等效弹性弹簧和阻尼器来实现。如图1所示,三维黏弹性人工边界处切向与法向的弹簧刚度与阻尼系数分别为:

图1 黏弹性人工边界示意图Fig. 1 Schematic diagram of viscous-spring artificial boundary

(1)

(2)

式中:cs和cp分别为SV波与P波波速;ρ为介质密度;G为剪切模量;ξ为拉梅常数;r为波源至人工边界点的距离;α、β为无量纲参数,分别取0.8、1.1。

1.2 基于黏弹性边界的等效节点荷载

在杜修力等[23]对地震动输入方法的相关研究中可以得出,地震波输入可以通过在黏弹性人工边界上施加等效节点力的形式来实现对波运动的转化,得到黏弹性边界节点的切线与法线方向的等效节点力可以表示为:

(3)

1.3 入射SV波时等效节点荷载求解

一个入射区域受到以任意角α1入射的SV波,域内任意一点(x,y)的各项位移时程记为ulx(t),左侧人工边界高度为H,下底边界和地表面长度为L,如图2所示。

图2 三维斜入射示意图 图3 三维SV波斜入射反射规律示意图Fig. 2 Schematic diagram of three-dimensional oblique incidence Fig. 3 Schematic diagram of oblique incident reflection law of 3D SV wave

左侧边界处任意一点相对于零时刻入射波的波阵面延迟时间分别为:

(4)

式中:Δt1、Δt2、Δt3分别为入射SV波、反射SV波、反射P波相对于零时刻入射波的波阵面的延迟时间。

当SV波以一定的角度斜入射至立方体截断区域时,根据波动理论中相关理论可知会产生相应的波形转换产生相应的SV波和P波,如图3所示。其中反射角分别为α1、β1,对应的波幅分别为B1、B2。求得左侧边界处的位移场和速度场为:

位移场:

(5)

速度场:

(6)

式(5)与式(6)求解为局部坐标中的位移场与速度场,为求得左侧边界上的应力,则需要将求解出来的位移场与速度场转化为全局坐标系[24]。

(7)

将左侧边界求得的位移场式(5)与速度场式(6)中入射SV波、反射SV波、反射P波任意一点零时刻相对于波阵面的延迟时间分别替换为式(8)即可按照同样的方法求得右侧边界下应力:

(8)

前后边界面应力求解方式与左侧边界求解类似,如图2中结构所示前侧边界与后侧边界处求解出来的应力为相反数。前侧边界入射的SV波、反射的SV波和反射的P波相对于入射波零时刻波阵面的延迟时间与求解应力为:

(9)

(10)

底边界面位移场仅由SV波构成,SV波入射波零时刻波阵面的延迟时间与求解应力分别为:

(11)

(12)

1.4 三维黏弹性静-动力统一人工边界转换方案

在结构动力响应分析中,因静、动力分析时边界条件不同。在进行整体分析中静力分析与动力分析之间会涉及边界条件切换,相关学者[19]研究表明,在边界条件切换过程中会产生一定的误差。为消除边界条件切换产生的误差对结果产生的影响,本文利用一种“单元生死”的技术处理来实现对大坝-库水-地基系统边界条件的转换。

具体实现的步骤总结如下:①地应力平衡:在动力分析时,为使重力荷载作用时地基仍保持初始状态,要进行地应力平衡。因此先对结构进行静力分析,即添加静力边界条件和重力荷载,导出地基的应力场。②求解支反力:将上一步导出的应力场作为初始边界条件,施加静力边界条件、重力荷载等静力荷载进行求解。边界结点的受力可分解为静力约束对地基的作用和地基对静力边界作用的叠加,最后导出支座反力。③动力分析:将得到的支反力作为静力约束的等效替代施加于边界。施加动力边界条件(黏弹性边界)以及静力(包含重力荷载、水压力等)、地震荷载等进行动力分析。

2 算例验证

2.1 SV波斜入射算例验证

本文构造一个三维弹性岩体模型来探究三维SV波斜入射输入方法的计算精度,如图4所示,该岩体计算区域边长为200 m的立方体块,网格尺寸取为5 m的立方体,该模型的上方为半空间自由表面,其余边界全部设置黏弹性边界。分析时时间增量步为0.01 s,计算时长2 s。位移时程曲线如图5所示。岩体介质的力学参数如表1所示。入射的SV波位移时间函数如式(13):

图4 三维弹性岩体有限元模型Fig. 4 Three-dimensional finite element model of elastic rock mass

图5 入射SV波位移时程曲线图Fig. 5 Displacement time history of incident SV wave

表1 介质的力学参数Table 1 Mechanical parameters of the medium

(13)

SV波分别在0°、15°及30°斜入射下地表中心点A的水平位移时程曲线如图6所示,与解析解的位移时程曲线进行对比发现:SV波斜入射三维算例计算结果的位移时程曲线与解析解的位移时程曲线吻合良好,且几乎不存在计算结果在位移归零时出现漂移的情况,可以说明本文中建立的模型及计算程序具有较好的模拟精度。

图6 斜入射时A点水平位移时程曲线Fig. 6 Horizontal displacements at point A under oblique incidence

2.2 静动力统一人工边界算例验证

为验证本文中三维黏弹性静-动力统一人工边界的模拟精度,选取50 m×50 m×50 m范围的有限元模型,网格尺寸取2 m的立方体,如图7所示,其力学参数为:弹性模量E=2.4×1010Pa,泊松比ν=0.25,质量密度ρ=2450 kg/m3,在顶面x向施加一个切向为q=1×108Pa的均布荷载,有限元模型侧面和底面施加黏弹性边界,入射SV波时程曲线与2.1节中输入一致,其位移时间函数如式(13),时间步长取0.01 s。

图7 岩体模型及参考点Fig. 7 Rock mass model and reference point

传统人工边界方案中对于静、动力分析过程中均采用黏弹性边界来实现结构的动力响应与本文中黏弹性静-动力统一人工边界求得顶面中点A处与底面中点B处位移时程曲线图,如图8所示,相比于传统方案,本文中采取的三维黏弹性静-动力统一人工边界转换方案消除边界切换过程中产生的误差影响,减小了与理论值之间的误差,其精度优势能够得到突显。说明本文中提出的计算方案与传统方案相比具有较高的模拟精度。

图8 模型顶部及底部中点X向位移时程曲线Fig. 8 X displacement time histories of the top and bottom midpoints

3 Nam Ngum 5拱坝地震响应分析

3.1 工程概况

本节选取Nam Ngum 5重力拱坝为研究对象。Nam Ngum 5重力拱坝的几何尺寸以及相应的设计资料采用文献[25]依据,并参照GB 51247—2018《水工建筑物抗震设计标准》[26]确定,最大坝高99.00 m,坝顶宽度6.00 m,拱冠梁底宽42.00 m,心线弧长234.838 m,弧高比2.42,拱顶中心角度92.794°。以此为依据建立Nam Ngum 5重力拱坝-库水-地基系统的三维有限元模型,如图9所示,共204538个节点,190505个单元。其中,坝体有189921个节点,6440个单元;库水14616个节点,12880个单元;地基189921个节点,171184个单元。

图9 Nam Ngum 5重力拱坝系统有限元网格图Fig. 9 Finite element mesh of Nam Ngum 5 gravity arch dam system

3.2 模型材料参数

坝体混凝土考虑为一种混凝土材料,并假定为线弹性,动弹模取为静弹模的1.5倍[26],阻尼比取5%;地基假定为均质线弹性介质不考虑材料阻尼。材料参数如表2所示,具体坝体、库水、地基相应的有限元模型以及具体尺寸如图10所示。

表2 材料力学参数Table 2 Mechanical parameter of materials参数坝体库水地基静态弹性模量/GPa33.15—20密度/(kg/m3)240010002700泊松比0.16—0.30图10 三维系统示意图Fig. 10 Three-dimensional system

3.3 计算荷载与边界条件

对于拱坝-库水-地基系统的作用荷载主要分为静荷载与动荷载两部分,其中静荷载为自重、静水压力等,动荷载为动水压力与地震荷载。参照规范可知该地区地震基本烈度为VI度,地震动峰值加速度为0.15g,本次计算中Δt取0.01 s,地震波是依据规范谱生成的随机地震波如图11所示。

图11 SV波加速度、位移、速度时程曲线Fig. 11 Acceleration, displacement and velocity of SV wave

采用三角级数法生成人工地震波,对于给定的功率谱密度函数Sx(ω),按式(14)可以方便的生成以Sx(ω)为功率谱密度函数、均值为0的高斯平稳过程a(t)[27]。同时为了更好地描述地面运动的非稳定性,采用包络函数f(t)乘以平稳过程a(t),具体如式(15)所示的人工地震波型。

(14)

x(t)=f(t)·a(t)

(15)

在系统模型中,截断边界采用黏弹性边界以模拟无限地基的辐射阻尼效应,采用声学单元及无反射边界条件模拟坝-库动力相互作用[28]。在进行结构整体静动力分析中,通过一种“单元生死”技术,实现静态约束边界到动态黏弹性边界的稳定切换,从而实现静动力统一计算。

4 结果分析

利用合成的三维SV波实现对Nam Ngum 5拱坝-库水-地基系统的位移、应力、动水压力的分析,由于SV波入射时存在临界角[29]且公式为α=arcsin(cs/cp)=32.3°,故本次选取0°、5°、10°、15°、20°、25°来进行工程实例分析。

4.1 位移分析

坝顶点横河向与顺河向位移曲线,如图12所示。由图可知,随着入射角度的增加横河向位移在逐渐减小,25°较 0°斜入射减小了8.32%,顺河向位移在逐渐增加,其中在地震波25°斜入射下顺河向位移达到0.15 m,这主要是因为输入的地震荷载为SV波(横河向位移)且本文中的地震动输入仅仅考虑平动分量,使得地震波在地表发生了反射使得地表结构产生了顺河向位移,导致水平位移逐渐减小顺河向位移逐渐增大。

图12 不同角度斜入射下坝顶点横河向、顺河向位移Fig. 12 Cross-river and along-river displacements with different angles of oblique incidence at the dam crest

4.2 应力分析

不同入射角度时坝体最大与最小主应力等值线图如图13和图14所示,由图可知,上游面主拉应力与下游面主拉应力均成对称分布,且随着入射角度的增大坝体主拉、压主应力逐渐增大,压应力极值集中出现在坝趾处,拉应力极值集中在坝踵处,规律性均与实际相符。分别提取坝踵、坝趾2个关键点的应力峰值曲线如图15所示。随着入射角度的逐渐增大,坝踵点、坝趾点2个关键点的应力均呈现增大的趋势,且在地震波25°斜入射下相应的主拉应力与主压应力达到了1.47 MPa与3.00 MPa。而入射角度为0°(即垂直入射)与其他入射角度相比,最大或最小主应力是最小的,这说明三维情况下SV波随着入射角度的增加对坝体动力响应的影响越大。

图13 SV波斜入射时坝体最大主拉应力等值线图Fig. 13 Contour map of the maximum principal tensile stress in the dam at oblique incidence of SV waves

图14 SV波斜入射时坝体最大主压应力等值线图Fig. 14 Contour map of the maximum principal compressive stress of dam at oblique incidence of SV wave

图15 坝体主应力峰值Fig. 15 Principal peak stress of dam body

4.3 动水压力分析

坝踵点处动水压力时程曲线如图16所示,由图可知:在相同位置随着入射角度的增加动水压力也在呈现出逐步增大的线性关系,当入射角度为25°时达到最大,此时最大动水压力峰值为0.33 MPa,最小动水压力峰值为-0.29 MPa。入射角度为0°时最小,此时最大动水压力峰值为0.05 MPa,最小动水压力峰值为-0.044 MPa。这说明三维情况下SV波随着入射角度的增加对大坝的动水压力也会逐渐增大。

图16 坝踵点动水压力Fig. 16 Hydrodynamic pressure at the dam heel

5 结论

本文将黏弹性人工边界与相应地震动输入方法相结合,推导了任意角度下三维SV波倾斜入射的等效节点荷载,构建了地震波斜入射的三维黏弹性静-动力统一人工边界,并且通过对ABAQUS的二次开发,编写了三维SV波斜入射的等效节点荷载子程序;构建了三维SV波斜入射拱坝-库水-地基系统分析模型,采用所建立的倾斜入射的模型与分析方法,将其应用于Nam Ngum 5重力拱坝的地震响应分析中,实现对模型中坝体在地震波倾斜入射波时的位移、应力、动水压力的分析,得到以下几点结论:

1)三维SV波斜入射与垂直入射下相比,重力拱坝坝体关键点地震响应规律存在明显差异。且SV波斜入射时坝体关键点的顺和向位移、主应力及动水压力与入射角度呈正相关,其中25°入射情况下的坝体关键点位移、主应力、动水压力最大。

2)由于输入的地震荷载为SV波(横河向位移),地震波在地表发生了反射使得地表结构产生了顺河向位移,导致水平位移逐渐减小顺河向位移逐渐增大。其中25°入射与0°入射相比,随着入射角度的增加横河向位移呈现逐渐减小的趋势,顺河向位移在逐渐增加,其中在地震波25°斜入射下顺河向位移达到最大,为0.15 m。

3)随着入射角度的增加坝体关键点的动水压力及主应力均在不断增加,其中25°与0°斜入射相比,在地震波25°斜入射下主拉应力与主压应力达到最大分别为1.47 MPa与3.00 MPa;动水压力也是25°入射时达到最大,最大动水压力峰值为0.33 MPa,最小动水压力峰值为-0.29 MPa。这表明三维SV波斜入射情况下,随着入射角度的增加对大坝的主应力与动水压力均会产生一定的影响。