重力坝动力分析黏弹性人工边界及其地震动输入处理方法

谯 雯

(华南农业大学水利与土木工程学院，广州 510642)

0 引言

在进行坝体—地基系统地震动力反应分析时，一个值得关注的问题就是地基辐射阻尼的影响。解决这个问题最有效的方法就是引入人工边界。人工边界的概念最早是由Alterman等[1]于1968年提出的。随后众多人工边界方法被提出，其中由Deeks[2]提出的黏弹性人工边界，因既能够模拟地基辐射阻尼的作用，又能够考虑地基的弹性恢复力，同时还具有良好的稳定性而得到了广泛的研究和应用。刘晶波等[3]在此基础上相继提出了静—动力统一人工边界、一致黏弹性人工边界及黏弹性人工边界单元的概念[4]。苑举卫等[5]、李浩然等[6]在大型有限元商用软件的平台上通过二次开发技术，实现了黏弹性人工边界的添加，取得了良好的模拟效果。黏弹性人工边界单元因便于编程实现，被广泛应用于结构—地基动力有限元问题分析中。

在坝体—地基相互作用问题中，地震动的输入是另一个关键问题[7]。对于大型水工结构而言，在考虑地基辐射阻尼的情况下，地震动输入方式通常按波动输入问题来考虑。刘晶波等[8]通过将地震波转化为人工边界上的等效节点荷载来实现地震动的输入。这种方法既适用于波的竖向输入，又适用于多点输入。

综上原因，本文基于黏弹性人工边界及其相配套的地震动输入方法，引入无厚度Goodman单元模拟黏弹性人工边界，推导了地震动输入方式在有限元中的实现过程。该方法在保证计算精度的前提下，不增加计算规模，保证了计算效率，并且在有限元中易于实现。通过经典数值算例验证了所编程序的正确性，最后将无厚度Goodman黏弹性人工边界单元应用于混凝土重力坝—地基—库水系统的地震动力响应分析中，并与传统的无质量地基模型进行比较，进一步验证了无厚度Goodman黏弹性人工边界单元及其与之相适应的地震动输入在实际工程应用中的强实用性、施加的方便性和计算结果的合理性。

1 黏弹性人工边界及地震动输入

基于球面波动理论推导出的三维黏弹性人工边界，通过在模型边界节点设置连续分布的并列弹簧—阻尼器来模拟无限介质。弹簧刚度系数和阻尼系数根据式(1)和式(2)取得：

(1)

(2)

式中：KT、KN分别为弹簧的切向和法向刚度系数；C为阻尼系数；G、ρ分别为介质的剪切模量、密度；cs、cp分别为剪切波速与压缩波速；αT为弹簧切向修正系数(取值范围0.5～1.0，推荐取0.67)[9]；αN为弹簧法向修正系数(取值范围1.0～2.0，推荐取1.33)[9]；R为人工边界至波源的距离。

黏弹性人工边界是将输入地震动转化为等效节点荷载：

(3)

(4)

2 黏弹性人工边界及其地震动输入在有限元中的处理

2.1 黏弹性人工边界单元在有限元中的处理

为了方便数值模拟，基于矩阵等效原理，本文采用文献[4]推导的利用等效实体边界单元替代弹簧—阻尼器原件，并且这种单元对于不规则边界具有更强的适应性。通过比较证明，单元厚度对计算结果影响甚微，因此本文利用无厚度Goodman单元模拟黏弹性人工边界。三维模型单元类型如图1所示。具体三维黏弹性边界模型的刚度矩阵、阻尼矩阵的推导过程可参见文献[10]。

(a)无厚度六面体单元

(b)无厚度五面体单元图1 Goodman单元的单元类型

2.2 地震动输入在有限元中的处理

假定黏弹性人工边界单元形函数为：

(5)

(6)

将式(6)代入式(3)可得到黏弹性人工边界侧边面上的等效结点荷载列阵为：

(7)

同理，将式(6)代入式(4)可得黏弹性人工边界底边上的等效节点荷载为：

(8)

根据已知输入地震波的位移和速度时程曲线，按式(7)～(8)即可求得人工边界上结点的等效结点荷载，从而实现黏弹性人工边界上地震动输入。

3 Lamb表面源问题的算例验证

为验证自编有限元程序中三维黏弹性人工边界实现的有效性，采用经典Lamb问题对其进行验证。模型长L=1，宽B=1，深度H=0.5，根据对称性原则，在水平方向上取原模型的1/4作为计算模型，如图2所示。详细计算参数见表1。通过在模型的四周和底面设置厚度为0的三维Goodman单元作为等效一致黏弹性人工边界单元，并将人工边界单元的最外层结点固定。在半空间表面处施加垂直集中荷载，荷载激励如图3所示。

图2 Lamb表面源问题计算模型示意图

表1 计算参数

图3 输入荷载时程

取距加载位置0点为0.2和0.4的两点作为观测点，计算结果与对比结果如图4所示。

(a)r=0.20 m处应力时程曲线

(b)r=0.40 m处应力时程曲线图4 计算结果与对比结果图

由计算结果图4显示，采用固定边界模拟时，脉冲波往复反射，观测点的位移产生了大幅度的波动，模拟结果与解析解相差甚大。而在黏弹性人工边界模型中，当脉冲波到达黏弹性人工边界后，能量几乎完全被边界的黏滞阻尼吸收，模拟结果与解析解非常一致。证明了本文基于无厚度Goodman单元的黏弹性人工边界能够很好地模拟无限地基辐射阻尼，同时也验证了所编程序的可靠性和有效性。

4 工程应用

4.1 有限元计算模型

由于实际工程中，重力坝在沿坝轴线方向设置了横缝，强震作用下横缝传递的剪力远小于惯性力。加之相邻坝段会由于坝高、自振周期的不同而形成单个坝段的振动。因此，对于重力坝的动力分析，一般不考虑沿着坝轴线方向的振动，常见做法是选取单个坝段进行动力分析。因此，本文选取福建省水口水电站19th坝段作为分析对象，对重力坝—地基—库水系统进行地震动力反应分析，分别对地基被施加黏弹性人工边界(模型Ⅰ)和固定边界(模型Ⅱ)两种情况建立了三维有限元模型，其中，施加固定边界的模型的地基范围较大。分析时坝体和地基均取为线弹性模型。假定库水是不可压缩的无阻尼理想流体，忽略振动中库水表面的影响，不计库底竖向位移对水体的影响。采用Westergard附加质量法实现库水动水压力对大坝地震响应的影响，采用直接滤频法求解特征方程。坝体地基网格划分如图5所示，库水位网格图如图6所示，大部分坝体、地基均采用三维八结点六面体等参单元离散，少部分坝身用六结点五面体单元过渡。在地基左右和底侧截断边界处各延伸一层厚度为0的Goodman单元作为等效边界单元来模拟黏弹性人工边界。对于施加了黏弹性人工边界的模型Ⅰ，共划分为5 852个结点，2 890个单元，而施加固定边界的模型II则划分为了6 972个结点，3 364个单元。

4.2 计算参数

重力坝的坝体以及地基假设为均匀各向同性体。在进行动力分析之前，需对重力坝进行静力计算，以确定坝体和坝基的初始应力状态。本文结合坝体和地基材料参数以及文献[11]的建议，参数取值见表2。

(a)模型Ⅰ：坝—地基系统有限元网格

(b)模型Ⅱ：坝—地基系统有限元网格图5 地基系列有限元网络

图6 库水位61.00 m的水体有限元网格图

4.3 地震波的选择与输入

这里选用位于Ⅲ类场地的EI Centro波实际地震记录，根据本工程相关资料，福州地区场地类别为Ⅲ类，处于抗震设防烈度7度区，现将EI Centro波峰值加速度调整到0.1g，时间步距0.02 s，调整后的EI Centro波加速度时程曲线如图7所示。

图7 调整后的EI Centro波加速度时程

4.4 动力计算成果分析

由于篇幅限制，本文只给出了坝顶65号结点和坝底3 109号结点的动位移时程，坝顶4号单元和坝底543号单元第一主应力和第三主应力时程，并同时给出了传统无质量地基固定边界的结果作为比较。

从图8～13对比可以看出：采用本文方法的计算模型进行计算，重力坝坝顶和坝基处的几个主要部位的顺河向动位移响应和动主应力响应与采用传统无质量地基固定边界条件下的计算模型对比，其计算的结果有明显的降低，其降低幅度大致在19%～55%之间。主要原因是由于无质量地基未考虑地基辐射阻尼的作用，导致反射波的存在，从而使得动力响应值偏大。从图中还可以看出：坝顶结点动位移响应较大，沿建基面位置附近的动力响应相对较小，其主要原因是地基及坝体对地震波的放大作用造成的。

表2 计算参数

图8 坝顶65号结点顺河向位移时程曲线

图9 坝底3 109号结点顺河向位移时程曲线

图10 坝顶4号单元第一主应力时程曲线

图11 坝底543号单元第一主应力时程曲线

图12 坝顶4号单元第三主应力时程曲线

图13 坝底543号单元第三主应力时程曲线

从以上结果比较可以看出，本文方法地基边界处对地震波动具有很好的吸收作用，波动在边界处的反射作用得到了明显的减弱，能够考虑地基辐射阻尼的影响，使得计算模型的计算结果更加符合实际，再次验证了本文动力有限元计算程序的可靠性和黏弹性人工边界施加的正确性。

5 结语

黏弹性人工边界不仅能够吸收外行散射地震波能量，而且能够模拟人工边界半无限地基介质的弹性恢复性能，被广泛地应用于实际工程。本文探讨了无厚度Goodman单元模拟黏弹性人工边界的可行性，推导了地震动输入转化为等效节点荷载在有限元中的实现过程。得出以下结论：

1)经典Lamb问题算例表明：使用一致黏弹性人工边界时，在接近加载点附近的结点和靠近边界的结点上，其位移计算结果与解析解很接近，计算的精度高。同时也可以看出，在对黏弹性人工边界的模拟上，采用本文的无厚度Goodman单元代替弹簧—阻尼器的方法，将使黏弹性人工边界的施加方式更为简单，保证计算效率。

2)重力坝—地基—库水系统进行地震动力反应分析实例表明：合理选取和模拟人工边界条件将直接影响动力分析结果的准确性，黏弹性人工边界条件原理简单、易于运用，具有良好的模拟效果，稳定性较好，而且可以在边界不规则的情况下使用，从而可为各类工程的抗震分析提供参考。