黏弹性边界的二次开发及其在地下结构抗震分析中的应用

材料参数，如弹性模量、泊松比、密度等;通过下式可以求得介质的纵向和剪切波速

式中：λ、μ为第一拉密常数和第二拉密常数;G、ρ为介质的剪切模量和密度。

将式（7）、（8）以代码的形式编写到Python程序中去。获取边界面上所有的网格结点，设成集合，并将总数记为h。通过式（1）～（4）分别计算该点需要施加的弹簧和阻尼元件的具体参数，利用reader函数将地震波的位移和速度函数分别以表格的形式读入到程序中来，每次从集合中取出一个结点，并利用for...in语句遍历表格中所有数据，以元组的形式表示出来。每次从边界的集合中取出一个结点，计算该点的相应参数，然后利用attend函数，通过while语句创建列表函数u（t-Δti），并通过函数的运算计算出该点对应的荷载函数，接着判断该点处于人工截面的具体位置，将随时间变化的荷载利用Python命令施加到该结点上去，同时将该点需要施加的边界条件，即弹簧和阻尼器施加上去，依次循环，直到每个节点都施加上相应的荷载和边界后，结束运行。具体流程如图3所示。

2.3 程序核心代码

该程序的主要功能是在有限的土体介质中施加地震波。只要给出土体介质的弹性模量、密度、泊松比等基本参数以后，便可以实现黏弹性边界和地震波的自动施加，程序部分核心代码如表1所列。

将程序写好后保存为.py格式的文件，当建立好模型并划分好网格以后，通过菜单栏的文件→运行脚本→选择相应Python文件运行即可。

3 程序验证

为了验证施加程序的有效性，在这里建立简单模型进行验证。在ABAQUS中建立长宽高均为50 m的立方土体，模型介质参数如下：弹性模量E=92.1 MPa，泊松比为0.25，密度为1 880 kg/m3。在底部输入加速度峰值为0.2g的天津波，地震波时长为19.19 s，计算时间取为22 s，取介质自由表面中心点的位移和加速度作为参考，时间增量取为0.01 s。为保证精度，网格大小为1 m × 1 m，土体外侧边界共有12 500个结点，程序仅运行80 s即将所有边界条件和地震波元件施加完毕。模型计算结果如图4、5所示。其中图4为土体入射面处和土体表面中心点处的加速度曲线对比图，图5为二者的位移对比图。

由波动的传播原理可知，对于单一的弹性介质而言，当波动传播到介质自由表面时其位移和加速度为入射时的2倍[3]。从图中可以看出，土体自由表面处的位移放大了两倍，同时加速度峰值也变为0.4g，同样被放大了两倍，由此可见利用此种方式施加的地震波具有良好的精度和效率，可以应用于地震波的施加和地下结构抗震的分析。

4 地下浅埋隧洞的有限元分析

为了考察本次二次开发在地下结构抗震中的效果，建立土体—隧洞结构有限元模型，对地震动在土体中的传播规律和土体—结构的相互作用进行研究。在这里选用某工程实地土层，具体参数如表2所列，根据经验公式计算，场地卓越周期为1.24 s。为了研究边界和地震波的施加效果，土体深度取为60 m，横向宽度取为90 m，纵向长度取为120 m，隧洞直径9 m，埋深6 m（浅埋），具体模型如图6所示。在模型底部输入水平横向地震动，加速度峰值分别为0.1g、0.2g、0.3g的天津波，其中0.1g天津波的加速度时程曲线和傅里叶谱如图7所示。同时建立同等参数下的自由场模型，并在底部输入水平横向加速度为0.1g的天津波作为对比。

4.1 土体竖直方向上的加速度变化和频谱特性分析

取0.1g天津波作用下自由场模型自由表面中心点的加速度，并绘出相应的傅里叶谱如图8所示。从图中可以看出，相比于输入时的地震波和傅氏谱（见图7），地震波的加速度峰值由0.1g变为0.18g，波形变化不大，这说明地震时土体本身对地震波的加速度有放大作用;从傅氏谱中可以看出，地表处接收到的地震波的低频成分增多，高频成分减少，傅里叶幅值的峰值频率由1.22 Hz变为0.81 Hz。这一数值正好与土体的基本频率吻合（f=1/T=1/1.24≈0.8 Hz），这是由于土体本身的基本频率比较小，放大了地震波中的低频成分，同时过滤了其中的高频成分。计算结果符合地震动在土体中的传播规律[21]。

当土体中有地下结构存在时，地震波在土体中的传播规律会发生改变，在这里选取0.1g天津波作用下的自由场模型和带隧洞结构模型，考察土体中心附近沿竖直方向的加速度放大系数变化，将土体自由表面處记为0 m，则基岩处为―60 m，每隔3 m取一个记录点，曲线如图9所示。

从图中可以看出，带隧洞土体和自由场土体的加速度放大系数基本一致，自底部至土体表面大致呈放大趋势，但是由于不同土层的性质有所不同，所以放大程度并不一致。其中在浅层0～15 m埋深范围内，由于隧洞结构的存在使得围岩中的加速度放大系数有所增加，同时也使得地表处的加速峰值变大。通过以上对比分析说明利用该方法施加的边界条件和地震波合理，地下浅埋隧洞结构的存在对围岩的加速度变化有一定影响。

4.2 加速度在水平方向上的变化分析

为了分析边界条件施加对模型的影响程度，考察0.1g天津波输入情况下自由场和带结构时土体自由表面中心处沿水平横向（共90 m）和水平纵向（共120 m）的加速度峰值变化趋势，每隔6 m取一个记录点，曲线如图10和图11所示。

对于土体表面水平横向加速度峰值变化而言，从图10中可以看出，在距离土体两侧边界的20～30 m范围内加速度峰值波动比较大，可见在这一范围内模型会受到边界效应的影响。而在中心附近的30 m范围内加速度峰值比较平稳，此时边界效应的影响减小。其中带隧洞土体表面加速度最小值为0.208g，最大值为0.215g，仅比最小值增大3%;自由场表面加速度最小值为0.200g，最大值为0.206g，同样仅比最小值增大3%，且数值波动出现的位置在距离土体边界2D～3D（D为结构直径）范围内;这一规律在水平纵向曲线（见图11）中同样体现出来，由此可以看出在靠近边界附近加速度峰值有轻微波动，但是误差较小，在可接受范围之内，在中间段则比较平稳，效果较好。经过多次试算结果表明，当将模型尺寸取为9倍的结构宽度时，可以消除两侧的误差带来的影响。

4.3 隧洞结构在地震动下的响应分析

隧洞结构在地震动输入时会产生变形甚至破坏，为了研究该种边界条件下地下结构的动力响应，取0.1g～0.3g天津波作用下隧洞结构纵向长度上顶部和底部的相对位移值和跨中的横截面环向应力值进行分析。

0.2g天津波作用下隧洞纵向长度上顶部和底部的相对位移沿隧洞纵向的变化曲线如图12所示，从图中可以看出，隧洞两端2D～3D（D为隧洞直径）范围内的相对位移值较大，在中部趋于平稳，这是由于隧洞结构在人为被截断以后两端的约束减弱，使得位移变大，这与实际的情况并不相符，在离开端部一段距离后这种影响才逐渐减弱。所以在建模计算时需要将隧洞纵向端部长度尺寸额外延长2～3倍的结构宽度，中间部分的稳定值才可以认为符合平截面假定。由于篇幅限制，本文仅给出0.1g～0.3g天津波作用下隧洞纵向跨中处的拱顶和拱底位移以及二者相对位移，如表3所列。从表中可以看出，隧洞的拱顶和拱底位移值稍大，但是相对位移比较小，这说明圆形浅埋隧洞在地震动作用下的地震响应虽然比较大，但是整体性比较好，抗震性能可以保证[22]，计算结果符合已有研究规律[23]。

为了研究隧洞横截面上的应力分布规律，分别取0.1g～0.3g天津波作用下隧洞跨中处横断面，将拱顶处记为0°，沿着顺时针方向每30°取一个观测点，即顺时针方向分别为0°（拱顶）、30°、60°、90°（拱底）、……330°、360°（与0°重合），横断面上各观测点沿顺时针方向的应力变化曲线如图13所示。

从图中可以看出，3种情况下隧洞的应力最大值均出现在120°～150°之间和210°～240°之间，而拱顶（0°和360°）和拱底（180°）处的应力值非常小，说明浅埋隧洞在V类围岩下衬砌的顶部和底部的应力值较小，斜下角处的应力值较大，这与文献[24]中的结论一致，分析结果合理，符合浅埋隧洞的地震响应规律。

5 结语

本文对ABAQUS的Python脚本接口进行了研究，成功开发出基于黏弹性边界的地震波统一施加程序，并对其精度进行了验证，最后对分层土体浅埋隧洞的动力响应进行了有限元分析，得到如下结论：

1）该方法可以实现黏弹性边界和地震波的统一施加，大大提高建模效率，同时可以有效消除波动在人工界面上的反射，计算结果符合波动在土体中的传播规律，可以应用于地下结构抗震分析中。

2）距离土体边界附近一定距离内的加速度峰值会出现波动，但是误差范围很小，当将模型尺寸取为9倍的结构宽度时，可以消除两侧的误差带来的影响。

3）隧洞结构纵向端部2～3倍宽度度范围内的计算结果偏大，需要通过额外增加相应的隧洞纵向长度取值来解决这一问题。

4）该方法为地下结构抗震分析提供了新的思路，也为进一步研究地下结构地震响应规律奠定了基础。

参考文献：

[1]    解祎. 基于黏弹性人工边界的地铁隧道地震响应分析[D]. 天津：天津大学，2013.

[2]    谷音，刘晶波，杜义欣. 三维一致粘弹性人工边界及等效粘弹性边界单元[J]. 工程力学，2007，24（12）：31-37.

[3]    廖振鵬. 工程波动理论导论 [M]. 2版. 北京：科学出版社，2002.

[4]    LYSMER J，KULEMEYER R L. Finite dynamic model for infinite media[J]. Journal of Engineering Mechanics，ASCE，1969，95（4）：859-877.

[5]    DEEKS A J，RANDOLPH M F. Axisymmetric time-domain transmitting boundaries[J]. Journal of Engineering Mechanics，ASCE，1994，120（1）：25-42.

[6]    章小龙，李小军，陈国兴，等. 黏弹性人工边界等效荷载计算的改进方法[J]. 力学学报，2016，48（5）：1126-1135.

[7]    刘晶波，吕彦东. 结构-地基动力相互作用问题分析的一种直接方法[J]. 土木工程学报，1998，31（3）：55-64.

[8]    刘晶波，李彬. 三维黏弹性静-动力统一人工边界[J]. 中国科学E辑：工程科学 材料科学 2005，35（9）：966-980.

[9]    LIU J B，DU Y X，DU X L，et al. 3D viscous-spring artificial boundary in time domain[J]. Earthquake Engineering and Engineering Vibration，2006，5（1）：93-102.

[10]  LIU J B，LIU S，DU X L. A method for the analysis of dynamic response of structure containing non-smooth contactable interfaces[J]. Acta Mechanica Sinica，1999，15（1）：63-72.

[11]  DU X L，WANG J T. Seismic response analysis of arch dam-water-rock foundation systems[J]. Earthquake Engineering and Engineering Vibration，2004，3（2）：283-291.

[12]  邱流潮，金峰. 地震分析中人工边界处理与地震动输入方法研究[J]. 岩土力学，2006，27（9）：1501-1504.

[13]  孙奔博，胡良明，付晓龙. 围岩类型及衬砌厚度对水工隧洞地震动力响应的影响分析[J]. 水利水电技术，2017，48（3）：19-24，45.

[14]  刘晶波，谷音，杜义欣. 一致粘弹性人工边界及粘弹性边界单元[J]. 岩土工程学报，2006，28（9）：1070-1075.

[15]  LIU J B，WANG V，LI B. Efficient procedure for seismic analysis of soil-structure interaction system[J]. Tsinghua Science and Technology，2006，11（6）：625-631.

[16]  张强，马永，李四超. 基于Python的ABAQUS二次开发方法与应用[J]. 舰船电子工程，2011，31（2）：131-134.

[17]  王辉，李廷春，王清标，等. 基于Python的ABAQUS二次开发及其在浅埋偏压隧道分析中的应用[J]. 现代隧道技术，2015，52（2）：160-165.

[18]  李勇. Python web开发学习实例[M]. 北京：清华大学出版社，2011.

[19]  黄胜，陈卫忠，伍国军，等. 地下工程抗震分析中地震动输入方法研究[J]. 岩石力学与工程学报，2010，29（6）：1254-1262.

[20]  杜修力，赵密. 基于黏弹性边界的拱坝地震反应分析方法[J]. 水利学报，2006，37（9）：1063-1069.

[21]  王国波，徐海清，于艳丽. “群洞效应”对紧邻交叠盾构隧道及场地土地震响应影响的初步分析[J]. 岩土工程学报，2013，35（5）：968-973.

[22]  GB 50011-2010，建筑抗震设计规范[S]. 北京：中国建筑工业出版社，2011.

[23]  王国波，陈梁，徐海清，等. 紧邻多孔交叠隧道抗震性能研究[J]. 岩土力学，2012，33（8）：2483-2490.

[24]  李凤稳，孙常新. 水工隧洞地震作用计算模型的适用性研究[J]. 水利水电技术，2017，48（4）：42-46.

[责任编辑 杨 屹]

推荐访问： 抗震 边界 弹性 地下 结构

---