地震动波动输入方法在高土石坝动力分析中的应用研究

魏匡民 陈生水 李国英 米占宽 沈 婷

(1.南京水利科学研究院岩土工程研究所，南京 210024；2.水利部土石坝破坏机理与防控技术重点实验室，南京 210029)

地震动输入是大坝地震安全评价的首要前提，也是大坝工程抗震分析中的一个薄弱环节[1].目前土石坝有限元动力分析中，大多采用一致性均匀输入，即将土石坝动力反应作为一个能量封闭的振动问题，计算中不计结构与地基的相互作用，将地震惯性力作用在坝体建基面上，则动力计算结果中加速度、速度、位移均为相对于建基面运动的相对值.一致性地震动输入方法不能反映河谷不同部位地震动差异和地震动传播的“行波效应”，该方法仅适用于尺寸小，刚度小，输入频率低的建筑物[2].对于200～300m 级特高土石坝结构，其尺寸、跨度与质量均非常大，在地震分析中应考虑坝体与地基的相互作用，这种相互作用既包括半无限地基对坝体动力特性的影响又包括坝体对动输入的影响.若不考虑坝体与地基的相互作用，则坝体与地基不存在能量交换，尤其不能反映地震能量向无限域的逸散现象，又称为“辐射阻尼效应”.辐射阻尼效应和自由场地震动输入机制密切相关，地基-坝体动力系统的地震响应包含自由场入射地震波以及由坝体和河谷产生的散射外行波，外行波在向无限山体和地基传播过程中由于几何扩散和地基阻尼逐步逸散，而实际数值模拟中计算范围仅能取1～2倍坝高，若在数值模型截断边界处采用固定约束或自由端，则一部分原本逸散到无限地基中的散射波会经截断边界反射回坝体-地基系统内，显著影响坝体结构的地震响应.

为了在数值计算中合理模拟地基的辐射阻尼现象，众多学者开展了截断边界处人工边界的研究，比较有代表性的有：廖振鹏[3]提出的人工透射边界，Lysmer[4]提出的粘性人工边界以及以Deeks[5]、刘晶波[6-9]、杜修力[9-12]等学者为代表的粘弹性人工边界.相较于粘性边界和透射边界，粘弹性边界在模拟地基的弹性恢复能力以及低频稳定性方面更有优势[7]，该方法物理概念清晰，数值稳定性高，具有很强的实用性.近些年，粘弹性人工边界在大坝地震分析中，尤其是混凝土高坝地震分析中，得到了较为广泛的应用，其中杜修力等采用粘弹性人工边界对小湾拱坝进行了动力反应分析[11]，张伯艳等[13]采用粘弹性人工边界对乌东德拱坝进行了非线性动力分析.在土石坝地震分析领域，孔宪京等[14-15]将非一致性波动输入方法引入了土石坝动力分析程序，研究了地震波频谱特性、坝体高度和地基模量对高土石坝地震反应的影响.总体来说，非一致波动输入方法在土石坝中的研究应用仍很少，目前土石坝工程规模不断增大，一些拟建的大坝如大石峡砂砾石面板坝、如美心墙堆石坝其坝高已达到250～300m 级，土石坝抗震计算中采用更合理的地震动输入方法有助于更准确评估大坝实际的抗震性能.本次研究在三维土石坝静、动力分析程序TOSS3D 基础上，引入了基于粘弹性人工边界的非一致性地震动波动输入方法，通过与解析结果对比验证了程序开发的正确性.以新疆大石峡面板砂砾石坝为例，比较研究了一致性和非一致性地震动输入下大坝在加速度反应、面板动应力、震后残余变形等方面的差异，评价了不同地震动输入方法对大坝地震反应的影响.

1 粘弹性人工边界单元

人工粘弹性边界是在地基边界处人为设置一种连续的应力型边界条件，其基本思想是将外行散射波作为波源问题，从理想介质中平面、柱面、球面波的标准方程推求结构边界面上法向和切向应力与位移、速度的关系式，以获得动态阻抗矩阵的局部阻尼和刚度集中系数，相当于设置阻尼器和弹簧的粘弹性边界.Deeks[5]提出了轴对称平面应变问题中的时域粘弹性边界.刘晶波等[16]提出了三维静动力统一人工边界的简洁表达，法向弹簧弹性系数KN，切向弹簧弹性系数KT，法向阻尼器阻尼系数CN，切向阻尼器阻尼系数CT的表达如式：

式 中 ，G为 地 基 剪 切 模 量 ；E为 地 基 弹 性 模 量 ；ρ为 地基 密 度 ；Cs为S波 波 速 ；Cp为P波 波 速 ；rb为 散 射 源到人工边界的距离；αN，αT分别为边界法向修正系数和切向修正系数，与波源问题中所采用的波型有关，从地基底部入射的地震波，若其竖向和水平向分量取以P波和S波传播的平面波时，αN，αT均取0.5.

粘弹性边界在有限元方法中通过粘弹性边界单元实现，文献[6]提出了二维的一致粘弹性人工边界单元，并推导了二维单元的刚度矩阵、阻尼矩阵，基于Nastran软件进行了结构的动力分析.本研究在此基础上开发了无厚度的三维的粘弹性人工边界单元，空间坐标系中z'方向为单元法向.

局部坐标系x'y'z'中，单元结点位移{ω}、速度与应力{σ}之间的关系可以写为

由于该单元是无厚度的 下盘结点1234坐标分别与上盘结点5，6，7，8 重合，单元形函数可选择为：

则局部坐标中单元刚度矩阵与阻尼矩阵为

其中，应变矩阵B和雅克比行列式|J|表达式为

式(4)和(5)均是在局部坐标系建立的，实际计算中应转换至整体坐标系.

2 人工粘弹性边界地震动输入

在数值计算中，需要首先确定半空间自由波场，计算在所有人工边界上需要施加的荷载.近些年，研究人员推导了地震P波，S波垂直和斜入射时的人工边界地震动输入形式[17]，由于地壳介质的密度由地表往下随地层深度而增大，按物理学中波在不同介质中传播的折射和反射定律，由地壳深部往地表传播的地震波，其入射方向将逐渐接近于垂直水平地表的竖向[2].所以，工程计算中一般采用垂直入射的地震动输入方式.即将地表面作为半无限体三维上表面，距离表面较远的平行于表面的一个假想截面，地震波从该面输入，设地震波不随坐标X和Y变化，即将地震波作为一维平面波从底部输入.如图1所示，在垂直向输入P波，沿X轴和Y轴方向分别输入S波.

图1 地震动输入的物理模型

本次研究采用文献[1]给出的各人工边界地震荷载表达.

3 非一致地震动输入程序正确性验证

本次研究基于南京水利科学研究院开发的TOSS3D 软件平台，实现了人工粘弹性边界和地震动的波动输入，为了验证本次开发程序正确性，采用经典的自由场地基P波，S波传播案例验证[10]，计算模型截取边长2000m 的正方体，模型及有限元网格如图2(a)所示，粘弹性边界单元剖分如图2(b)所示.地基密度2630kg/m3，弹性模量E=32.5GPa，泊松比=0.22，动力计算时取时间步长为 0.006s.

图2 自由场有限元模型

地震波从三维模型底部垂直输入，输入地震波位移、速度和加速度时程表达式如下

位移时程：

速度时程：

加速度时程：

半无限自由场底部输入位移过程、速度过程、加速度过程分别为u0(t)，u0(t)，¨u0(t)，则任 意 一 点 任意时刻的位移解析解为：

显然地表位移、加速度解析解为考虑行波延迟后放大2倍的入射位移速度和加速度时程.图3为P波入射时模型表面中心点(图1中A 点)的位移、加速度时程，可以看出，模型顶部中心点动位移和加速度幅值是输入值的2 倍，解析解与数值解能够较好吻合.图4 为S 波入射时模型表面中心点(图1 中A点)的位移、加速度时程，可以看出，模型顶部中心点动位移和加速度幅值也均为输入值的2倍，解析解与数值解能够较好吻合.可见，本次研究开发的软件具有良好的模拟精度.

图3 P波输入时地表点解析解与数值解对比

图4 S波输入时地表结点解析解与数值解比较

4 大石峡砂砾石面板坝地震分析研究

4.1 工程概况

大石峡砂砾石面板坝坝高247m，正常蓄水位1700m，为拟建的世界第一高面板坝，场区地震基本烈度7度，拦河坝设防烈度为8度，设计地震100年超越概率2%地震动峰值加速度为387gal.典型剖面材料分区如图5所示，有限元模型如图6所示.单元总数73167个，结点总数71554个.填筑和蓄水过程分51个荷载级模拟.

图5 大石峡砂砾石面板坝材料分区

图6 大石峡面板砂砾石坝三维有限元网格

4.2 计算参数与输入地震动

数值模拟中，筑坝料静力计算模型采用土石坝工程中广泛使用的“南水”模型[18]，动力模型采用沈珠江动力模型[19]，计算参数见表1～2，面板和趾板均为C30混凝土.

本次计算采用100年超越概率2%地震动输入，地震安评部门提供了3条地震动时程，由于并未直接规定这3条地震动的输入方向，所以计算过程中有6种荷载组合，计算结果表明各荷载组合计算坝体动反应规律基本相同，但量值上有差异，本文选最不利荷载组合为例进行分析，该组合中坝轴向、顺河向、竖直向分别输入的加速度时程如图7所示，其中竖直向加速度峰值为水平向的2/3.一致性地震动输入计算中地震动从基岩面输入，非一致性波动地震动计算中，地震动以面力形式输入.

表1 大石峡筑坝材料南水模型参数

表2 沈珠江动力模型参数

图7 输入地震动时程

4.3 计算结果分析

1)加速度反应

图8和图9分别给出了地震动一致性和波动输入下坝体的在坝轴向、顺河向、垂直向的加速度放大倍数.一致性输入条件下坝体在轴向、顺河向、竖直向的加速度放大倍数最大值分别为2.83，2.92，3.03，波动输入条件下，坝体在轴向、顺河向、竖直向的加速度放大倍数最大值分别为2.11，2.70，2.87，可以看出，一致性输入方法计算的坝体反应要大于波动输入方法.两种方法计算的坝体加速度极值位置均位于坝体，但分布位置略有不同，从放大倍数分布还可以看出，采用波动输入坝顶的“鞭梢效应”较一致性输入弱，此现象与粘弹性边界吸收了外行波有关.

图8 一致性地震动输入加速度放大倍数

图9 地震动波动输入加速度放大倍数

2)动位移反应

图10和图11分别给出了地震动一致性和波动输入下坝体的在坝轴向、顺河向、垂直向动位移极值分布.一致性输入条件下轴向、顺河向、竖直向动位移最大值为22.6cm，40.7cm，14.3cm，波动输入条件下轴向、顺河向、竖直向动位移最大值为42.5cm，40.9cm，31.2cm.由于波动方法计算所得的动位移为坝体的绝对位移，一致性输入方法计算的动位移是坝体相对于基岩面的相对位移，即基岩面的动位移为零.所以，从量值上波动输入条件下动位移更大，且分布规律与一致性地震动输入有较大差异.

图10 一致性地震动输入动位移极值

图11 地震动波动输入动位移极值

3)震后残余变形

一致性地震动输入和波动输入，大坝残余变形分布规律基本一致，即坝轴向残余变形表现为两岸向河床方向的挤压变形，顺河向残余变形基本表现为指向下游，大坝震陷极值出现在坝顶.篇幅所限，图12给出了波动输入下坝体的在坝轴向、顺河向、垂直向的残余变形分布.从计算结果表明，一致性地震动输入，坝轴向残余变形指向右岸和左岸最大值分别为10.1 cm 和13.7cm，顺河向位移指向下游最大值为62.7 cm，震陷最大值为92.1cm；波动输入方法坝轴向残余变形指向右岸和左岸最大值分别为9.5cm 和9.2 cm，顺河向位移指向下游最大值为57.1cm，震陷最大值为76.5cm，可见，波动输入方法计算的坝体残余变形较一致性输入小.

图12 地震动波动输入坝体残余变形

4)防渗体应力变形

地震过程中混凝土面板内会因为地震位移产生一定的动拉应力和动压应力.本次计算结果表明，一致性和波动输入方法下，面板动应力分布规律基本一致，即顺坡向压应力、拉应力最大值出现在面板中部偏向坝顶位置，坝轴向压、拉应力最大值出现在面板靠近两岸偏向坝顶位置.一致性地震动输入情况下，面板最大轴向动压、动拉应力分别为5.65MPa、5.26 MPa；顺坡向动压和动拉应力最大值分别为8.64 MPa，7.98MPa.地震动波动输入情况下，面板最大轴向动压和动拉应力分别为5.24MPa、3.94MPa；顺坡向动压、动拉应力最大值分别为7.44MPa，6.22 MPa.一致性地震动输入情况下地震期面板周边缝三向变位分别为错动39.8mm、沉陷43.4mm、张开21.6mm.地震动波动输入情况下震后面板周边缝三向变位分别为错动38.6mm、沉陷40.8mm、张开19.9mm.

表3为一致性输入和波动输入条件下大坝动力响应和残余变形汇总，可以看出，地震动波动输入面板轴向压、拉应力分别较一致性输入减小了约7.3%，25.1%，顺河向压、拉应力分别减小了13.8%，22.1%.坝体震陷，波动输入较一致性输入减小了约17%，地震期周边缝三向变位错动、沉陷、张拉分别减小了1.2mm，2.6mm，1.7mm.可见，一致性输入偏大地预测了坝体的动力反应.

表3 波动输入与一致性输入计算结果汇总

5 结 论

1)波动输入方法在模拟地基辐射阻尼和地震波传播的“行波效应”方面较传统的一致性输入方法更为合理.

2)将粘弹性边界和均质场地地震动输入方法嵌入了TOSS3D 有限元分析软件，通过与解析解比较验证了开发软件的正确性.

3)以拟建的247m 大石峡砂砾石面板坝为例，对比研究了地震动一致性输入方法和波动输入方法坝体加速度反应、动位移、残余变形以及防渗体应力变形方面的差异，结果表明，由于考虑了半无限地基的辐射阻尼，波动输入方法计算所得的坝体动力反应和残余变形整体较一致性输入小，且分布规律也有所差异.从工程角度来说，传统的一致性地震动输入方法是偏于安全的，但地震动波动输入方法更有助于客观评价大坝的抗震能力.

4)目前拟建的土石坝坝高已达到250～300m级，由于其质量大、跨度长坝体动力分析中需建立坝体与地基的相互作用的系统，以更合理模拟坝体与地基的能量交换行为.