隧道地震响应中的人工边界和地震动输入方法研究

2022-07-21丁祖德陈誉升

地震工程与工程振动 2022年3期

丁祖德，陈誉升，资昊

（昆明理工大学建筑工程学院，云南昆明 650500）

引言

动力时程分析法能够反映岩土体－结构动力相互作用过程以及材料的非线性行为，已广泛应用于隧道及地下结构的地震响应分析中［1］。在开展地下结构动力时程分析时，需要考虑无限域地层与结构的相互作用，但受限于计算资源，通常只截取有限域进行研究。为使有限域地层具备模拟无限域地层的辐射阻尼［2］，需要在人工截断的地层周围施加人工边界。此时，地震动输入方法及动力人工边界的选择，将直接影响地下结构动力响应结果精度和可靠性。

目前，国内外学者针对动力人工边界以及地震动输入方法开展了较丰富的研究工作，发展了如粘性边界［3］、粘弹性边界［4］、等位移边界［5］等多种人工边界模型。刘晶波等［6］推导了与波动输入法相匹配的二维粘弹性边界。杨喻声等［7］验证了等位移边界在弹性基岩假设条件下的有效性，简化了常规波动分析方法。李嘉瑞等［8］对比分析了采用不同人工边界对地下结构动力计算结果的影响，建议侧边界采用绑定边界或竖向约束边界，根据下卧基岩场地条件选择底部边界类型。赵武胜等［9］和马笙杰等［10］分别研究了无限元边界、粘弹性边界与地震动输入方法之间的影响机制，结果表明加速度时程输入将导致底部人工边界失效，且导致位移场与理论值产生严重偏差。

已有研究基本上都是从单一人工边界的地震动输入原理展开阐述，且大都基于自由场模型，从位移角度来验证人工边界及地震动输入的有效性［9－11］，而从结构力学响应角度进行验证［8，12］，对地震动输入方法与动力人工边界的匹配性问题的探讨较少。此外，人工边界的施加较为繁琐，在隧道动力响应有限元模型前处理过程中的工作量巨大。基于此，文中依托ABAQUS 有限元软件平台，结合粘弹性边界、无限元边界和等位移边界等动力人工边界理论，基于Python 编程开发相应人工边界下地震动输入程序，以实现模型边界及地震动的自动施加。进一步结合自由场模型和土－隧道相互作用模型，对比振动输入和波动输入方法及常用人工边界的性能，验证文中开发程序的有效性，最后探讨2种地震动输入方法的差异性。

1 人工边界及地震动输入的程序实现

目前，隧道动力响应计算模型中，常用的人工边界有远置边界、无限元边界、粘弹性边界和等位移边界等［8］，常用的地震动输入方法有波动输入法和振动输入法［13］。其中，波动输入法是基于弹性基岩假设，并结合动力人工边界理论，把地震动荷载转换为等效节点力后逐一施加到边界节点上。振动输入法是基于刚性基岩假设，通过在计算模型的底边界上施加位移、速度或加速度时程来实现地震动输入。

波动输入法在施加人工边界时，需要逐一处理边界节点，前处理工作量巨大。鉴于此，文中结合粘弹性边界、无限元边界、自由度绑定边界和竖向约束边界的基本理论，采用Python 编程开发上述人工边界的自动施加程序。为方便使用，进一步将该程序内嵌至ABAQUS 软件中，新增了专门的隧道动力响应辅助计算模块（注：原ABAQUS 软件中无此模块，文中采用Python 语言编程，结合其二次开发功能实现），实现了程序界面操作功能，如图1 所示。接下来，具体阐述不同边界条件下波动输入法的程序实现流程。限于篇幅，人工边界理论推导详见文献［5－6］、［14］。

图1 人工边界及地震动输入自动化程序Fig.1 Automatic program of artificial boundary and ground motion input

1.1 粘弹性边界下的地震动输入

粘弹性边界通过在模型边界处设置并联的阻尼器元件和弹簧元件来模拟无限域地基的辐射阻尼和弹性恢复力，避免了高频失稳及低频漂移的问题。针对粘弹性边界，常采用波动输入法将地震波转化成等效节点荷载输入。其等效节点力计算如式（1）［10］：

式中：Fb为边界节点上的等效荷载；K、C分别为施加到边界上的弹簧—阻尼元件的刚度系数、阻尼系数，其计算公式与取值参考马笙杰等［10］研究；Ab为节点等效面积，分别为边界节点的自由场位移、速度及应力。

考虑到等效节点荷载是基于边界节点逐个施加的，且每个节点处均需要计算对应的等效节点力时程，文中采用Python语言编制粘弹性边界及地震动输入程序，程序流程如图2所示。

图2 粘弹性边界实现及地震动输入流程Fig.2 Implementation of viscoelastic boundary and seismic input process

1.2 无限元边界下的地震动输入

无限元边界是基于无限单元实现的静动力人工边界，在静力分析时通过有限长的无限单元来模拟半无限地基，能反映地基无穷远处位移为零的条件，在动力分析时通过在无限单元中内置阻尼器来实现粘性边界的吸收效果［14］。

无限元边界与粘弹性边界的地震动施加原理相似，同样采用Python 语言编程实现无限元边界及地震动等效节点力时程的输入，程序流程如图3所示。等效节点荷载计算式如下：

图3 无限元边界的地震动输入流程Fig.3 Flow chart of ground motion input with infinite element boundary

式中：Fb为边界节点上的等效荷载；Ab为节点的等效面积；C为阻尼系数分别为边界节点的自由场速度及应力。

1.3 等位移边界下的地震动输入

等位移边界分为自由度绑定边界和竖向约束边界。其中，自由度绑定边界通过在模型两侧设置运动绑定约束，约束两侧边界上节点的自由度防止其旋转，迫使土层侧边同一高度位置的节点同时运动，从而避免侧边界对波的反射［5］。竖向约束边界通过约束两侧截断边界的竖向自由度，以模拟地基模型在剪切波作用下的剪切动力响应。

等位移边界下的地震动输入可以采用波动输入法，即在两侧边设置等位移边界，底边设置无限元或粘弹性边界来实现，也可以直接采用振动输入法。结合上述粘弹性边界、无限元边界的Python 程序，同样开发了等位移边界及地震动输入（图1）。

1.4 远置边界下的地震动输入

在结构动力响应分析时间内，由于远置边界距离结构足够远，模型边界反射波不会对结构响应产生影响，保证了动力响应分析的可靠性。模型边界宽度根据下式确定［15］：

式中：L为模型边界距离结构所在位置的距离；T为结构动力响应分析时间；C为地震波在土体介质中的传播速度。

根据式（3）可知，当地震波持时较长、土体介质波速较大时，需要巨大的模型尺寸才可避免边界对结构位置的影响，并不具备实用性。远置边界在动力分析时间内是准确的，常被作为波动问题的解析解，以验证其他边界的有效性。

2 波动输入法人工边界有效性验证

本节从自由场模型、土层－隧道相互作用模型2 个方面，分别建立二维地震响应计算模型，对无限元边界、粘弹性边界、自由度绑定边界、竖向约束边界及相应的地震动输入方式进行系统性验证。针对波动输入方式下的模型验证时，各人工边界的选择见表1。

表1 波动输入方式下的边界选择Table 1 Boundary selection under wave input methods

为方便对比各边界条件下的计算结果，从模型底部垂直输入SV波，其表达式如下［16］：

2.1 自由场模型验证

考虑不同人工边界的二维土体自由场模型如图4所示。模型尺寸取200 m×50 m，单元网格尺寸为1 m×1 m。土体单元采用线弹性的4节点平面应变单元（CPE4R）模拟，计算参数见表2。计算中，以远置边界自由场模型计算结果及理论值来验证人工边界模型的有效性。为了满足远置边界条件，根据土体参数，由式（3）计算的远置边界自由场模型宽度至少为589 m，故取该边界下模型宽度为1 000 m。按上述参数建立有限元模型后，采用文中开发的动力计算辅助程序施加各人工边界及输入地震动，开展自由场模型动力响应分析。

表2 土体物理力学参数Table 2 Soil physical mechanics parameters

图4 不同人工边界的自由场模型Fig.4 Free-field model considering different artificial boundaries

为分析动力响应结果，在模型底部及顶部中点设置A、B两个监测点。图5 为各边界下监测点的水平位移时程结果。图5中的理论解显示，不考虑土体介质阻尼的情况下，垂直入射的SV波经过0.51 s后到达自由表面和反射波叠加，使自由表面监测点A的位移振幅相比入射波放大2倍，而波形保持不变。反射波经0.51 s后传至模型底部，引起模型底部振动，之后被底部人工边界吸收而不再引发有限域的振动。各人工边界条件下，模型监测点处的水平位移时程曲线均与理论曲线及远置边界模型结果曲线吻合较好。其中，粘弹性边界在波动结束后有一个小幅震荡，与文献［10］结果类似，在实际计算中对结果产生的误差可忽略不计。可见，上述人工边界均能较好地吸收地表反射波，具有良好的计算精度，也说明所开发的人工边界及地震动输入程序是有效性的。

图5 各人工边界模型监测点A、B处的水平位移时程曲线（波动输入）Fig.5 Time-history curve of horizontal displacement at monitoring points A and B of each artificial boundary models（wave input）

2.2 隧道地震响应模型验证

为验证各人工边界对隧道地震响应计算的有效性，在自由场模型基础上，进一步建立考虑土层－隧道动力相互作用的横向地震响应模型见图6。其中，将远置边界模型结果作为解析解。模型尺寸和土层条件同自由场模型，隧道埋深取为10 m，隧道外半径取为3 m、衬砌厚度为0.35 m。衬砌采用线弹性梁单元（B21）模拟，计算参数见表3。

表3 衬砌物理力学参数Table 3 Lining physical mechanical parameters

图6 不同人工边界的自由场模型Fig.6 Free-field model considering different artificial boundaries

同时，引入圆形隧道在横向剪切地震作用下的4种拟静力解析解作为衬砌内力数值解的对比（表4）［17－20］。剪切波作用下的圆形隧道解析解，通常是基于土－隧道系统周围的远场地层处于纯剪状态的假设［12］，将剪应力等效为沿45o方向产生的拉、压主应力，此时衬砌发生椭圆化的剪切变形，将土体的变形施加到衬砌结构上，即可确定衬砌结构在横向地震动作用下的内力。其中，土体变形由土体所受剪应力的大小以及土－隧道结构的相对刚度来决定。对于半径为R 的圆形隧道，隧道周边的远场土层受到的剪应力通过式（5）计算得到，而土－隧道结构的相对刚度则采用式（6）中的压缩系数C和柔度系数F量化表示［20］：

表4 横向地震作用下圆形隧道结构内力解析解Table 4 Analytical solution of internal force of circular tunnel structure under transverse earthquake

式中：R为隧道外半径；Es、νs分别为土体弹模、泊松比；γc为隧道周围2R范围内土体平均剪切应变峰值；τ为周围土层在剪切地震作用下产生的剪应力；C、F分别为压缩系数和柔度系数；El、νl、Al、Il分别为衬砌的弹模、泊松比、截面积和惯性矩。

需要指出的是，隧道结构与周围土体介质的接触面是比较粗糙的，不易出现完全滑移的极端条件［12］，因此，本文验证计算仅考虑土与隧道结构无相对滑移的情况。

图7 为各边界条件下衬砌拱肩（θ= 45o）处的内力时程曲线。从图7 可看出，各人工边界模型的计算结果在0～1.5 s 之间与远置边界模型的结果基本吻合，其中，粘弹性边界的计算内力结果在1.5～2 s 之间有较小波动；自由度绑定边界和竖向约束边界的计算内力结果在2.5～3 s之间有较小波动。总体而言，各人工边界计算结果的整体误差较小，在实际分析中可忽略不计。

图7 不同人工边界典型截面内力时程曲线（波动输入）Fig.7 Time-history curve of internal force of typical section with different artificial boundaries（wave input）

结合土－隧道相互作用模型的参数，通过式（4）可计算得到土－隧道系统的压缩系数C和柔度系数F分别为0.002 5、1.476。其次，利用自由场模型进行地震响应分析，得到隧道所处地层范围内的平均剪切应变峰值γc为0.002 18。据此，可计算出地震作用下隧道所处地层平均应变峰值时刻衬砌截面内力的解析解，并与动力时程分析中的峰值内力结果对比，如图8所示。

由图8中可知，Wang解、Bobet解和Park 解的全截面动轴力结果比较接近，而Penzien 解严重低估了截面动轴力。3种典型人工边界下的衬砌截面动轴力分布完全一致。与3种典型人工边界下的动轴力相比，解析解在45°及135°截面及其附近约45°范围内偏小，最大相差8.5%，其它截面位置的动轴力分布与人工边界下的基本一致。

图8 波动输入的峰值时刻衬砌截面动内力分布（各人工边界对比）Fig.8 Distribution of dynamic internal force of lining section at peak time of wave input（comparison of artificial boundaries）

衬砌截面动弯矩的4 种解析解略有差异，其中，在±45°轴线的截面上差别最明显，弯矩值最大相差约16.02%，Wang解最大，Penzien解次之，Bobet及Park解最小且完全一致。各人工边界下的衬砌截面动弯矩分布完全一致，不同边界下的动弯矩在0°～180°范围内与Penzien解基本一致，而在180°～360°范围内则与Bobet解及Park解吻合。对比各边界下的数值解与解析解可知，各人工边界下的内力分布及内力值均与解析解接近，表明本文开发的人工边界及地震动施加程序可用于隧道地震响应分析。对比上述边界还可看出，相较于无限元边界和粘弹性边界，等位移边界有相当的计算精度，而且边界施加简单，更实用便捷。

3 振动输入法人工边界有效性验证

3.1 自由场模型验证

同样采用式（4）所示的剪切波进行验证。考虑到振动输入会使地震波在有限域中来回反射，为了反映这种特点，将动力分析时长调整为6 s。并以远置边界计算结果为解析解，按式（3）计算得到远置边界模型宽度需大于1 176.7 m，因此，取自由场模型宽度为2 000 m。结合上述分析结果可知，等位移边界（如自由度绑定边界）与其他人工边界相比，施加更简便。为节约篇幅，本节采用自由度绑定边界开展振动输入法的有效性验证，其自由场模型和隧道地震响应模型的尺寸及参数设置分别沿用2.1节和2.2节。

自由场模型计算结果如图9 所示，与波动输入法不同，振动输入法的模型底部由输入的位移时程控制，不因自由表面反射波的到来而改变，但下行反射波会经过底部边界再次向自由表面反射，引起模型自由表面位移出现反复的波动。

图9 不同人工边界模型监测点A、B处的水平位移时程曲线（振动输入）Fig.9 Time-history curves of horizontal displacement（vibration input）at monitoring points A and B of different artificial boundary models

3.2 隧道地震响应模型验证

图10为2种边界条件下的衬砌拱肩（θ=45o）部位内力时程曲线。与波动输入法明显不同，由于地震波在模型中的来回反射，土层－隧道体系受到了循环的地震动荷载作用，导致隧道结构内力也出现了循环响应特征。

图10 不同人工边界下的衬砌典型截面内力时程曲线（振动输入）Fig.10 Time-history curve of internal force of typical lining section under different artificial boundaries（vibration input）

图11为解析解与振动输入法数值解的对比结果。从图中可以看出，对于衬砌截面动轴力，除Penzien解以外的3种解析解与数值解都比较接近。对于衬砌截面动弯矩，数值解与解析解均相差较小。由此说明，自由度绑定边界对振动输入法是有效性的。

图11 峰值时刻衬砌截面内力分布曲线（振动输入）Fig.11 The internal force distribution curve of the lining section at the peak moment（vibration input）

4 波动输入与振动输入差异性

在地下结构抗震分析中，由于输入的地震波持时较长，导致结构的地震响应可能由多个波场共同决定。因此，为探讨这种多波场条件下振动输入与波动输入的区别，将土体弹性模量调整为1 GPa，其余模型参数不变，则土体介质剪切波速cs=438.5 m/s，由此可计算出入射波从底边界入射后，仅需0.228 s便可经过自由表面反射而传回底部，此时原先的入射波场振动还未结束，随即产生了多波场叠加效应。本节同样以自由度绑定边界为例，如图12（a）所示，波动输入模型底部A点在0.228～0.5 s间的位移便是由入射波与反射波叠加形成的。而振动输入模型底部受位移时程控制，不因自由表面反射波的到来而改变。如图12（b）所示，振动输入的模型底部反射波还未再次到达自由表面时（0.342 s之前），自由表面仅发生一次反射，其位移较入射位移放大2倍；当自由表面再次发生反射（0.342 s之后），此时多波场叠加，位移场放大倍数超过2倍，土体地震响应加剧。由于底部人工边界对反射波的具有吸收功能，采用波动输入的模型顶部B点仅发生一次反射作用，位移场放大2倍。

图12 2种地震动输入下监测点A、B的位移时程曲线Fig.12 Displacement time-history curves of monitoring points A and B under two types of ground motion inputs

如图13、图14所示，由于地层中多波场叠加效应使得位移场增强，振动输入时的隧道衬砌的峰值内力大于波动输入模型的计算结果。可见，振动输入与波动输入对结构的地震响应有明显的影响，振动输入会使地层结构产生更剧烈的循环加卸载行为，会放大土层及结构的动力响应。因此，在进行地下结构地震响应计算时应根据具体的场地条件选择合理的地震动输入方法，当下卧基岩满足刚性条件假定可考虑使用振动输入法，而绝大多数地下结构的场地条件宜采用波动输入法。