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

窦远明 范俊超 王建宁 鞠培东 宋明轩 李景文

摘要 在对地下结构进行抗震分析时，土体边界条件和地震波的施加方法直接关系到运算结果的精准程度。为了使地下结构抗震分析建模更加高效，分析结果更加合理，对边界和地震波的施加算法进行了程序化设计，利用Python语言对ABAQUS进行了二次开发，编写了黏弹性边界和地震波统一自动施加程序，建立了土体—隧洞结构相互作用的三维有限元分析模型。结果表明：该方法可以实现黏弹性边界和地震波的快速自动施加，能够很好地模拟波动在土体中的传播规律;在靠近土体边界附近一定范围内的加速度峰值有3%左右的误差，当模型尺寸取9倍的结构宽度时可以消除这一影响;隧洞结构纵向端部2～3倍结构宽度范围内的计算结果偏大。

关 键 词 地下结构;黏弹性边界;二次开发;有限元分析

中图分类号 TU93     文献标志码 A

Abstract In the seismic analysis of underground structure， the boundary conditions of soil and the seismic wave application methods are directly related to the accuracy of calculation results. In order to make the seismic analysis more effictive and more reasonable， the application algorithm of boundary and seismic waves was programmed and a secondary development of ABAQUS with Python was done. A finite element analysis model with 3D nonlinear soil-tunnel interaction was established by using the unified automatic impose program for viscoelastic boundary and seismic waves. The results show that this method can apply the viscoelastic boundary and the seismic waves automatically and can reflect the seismic law of foundation reasonably. There is a deviation in peak acceleration deviation （about 3%） near the model border and it can be eliminated when the FEM size is 9 times larger than the width of structure. The calculation results in the scope of 2D～3D （D represents the width of structure） width of the tunnel structure end are deviated to be larger.

Key words underground structure; viscoelastic boundary; secondary exploration; finite element analysis

0 引言

在對地下结构进行有限元动力分析的过程中，需要把半无限的土体转换为有边界的空间，而波动在传播的过程中遇到边界会产生反射，这与实际情况并不相符，所以在建立模型时必须对人工边界进行处理，使之符合实际情况[1-2]。当前常用的人工边界条件中，有透射边界[3]、黏性边界[4]、黏弹性边界[5]等，其中透射边界为位移型边界，在多次透射的情况下精度较高，但是容易出现高频震荡[6]，且不容易在有限元软件中实现;黏性边界为应力型人工边界，但是其仅考虑了对散射波能量的吸收，并没有考虑到边界处介质的弹性恢复能力，因此容易导致低频失稳问题[7];而黏弹性边界则克服了上述缺点，能够很好地模拟地基的弹性恢复力和辐射阻尼效应，具有较高的精度[8-9]。除此之外，杜修力等[10-11]结合黏弹性边界条件，将地震作用的总波场进行分解，提出通过荷载的方式在边界处施加地震作用，取得了很好的效果。

虽然人工边界设置和地震波施加的理论取得了很多成果，但是如何将这些方法应用到实际的有限元分析中去依然是目前需要探讨的问题[12]。其中孙奔博等[13]采用combine14单元通过耦合结点的方式在Anasys中实现了黏弹性边界的施加，并对水工隧洞进行了有限元分析;谷音等[2]、刘晶波等[14-15]开发了等效黏弹性边界单元方法，用等效实体单元来替代弹簧—阻尼器原件，在有限元软件中实现了黏弹性边界的施加;章小龙等[6]利用Fortran软件对ABAQUS进行了二次开发，实现了地震动等效荷载的施加。但是对通过Python语言编程实现黏弹性边界和地震动等效荷载的统一施加的方法研究较少。

Python是一种面向对象、解释性、动态数据类型的高级程序设计语言[16-17]，它的语法简洁清晰，易于理解和开发。同时它还具有高效的数据结构，提供了一种简单但很有效的方式进行面向对象的编程[18]，这使其成为大多数平台上应用于各个领域的理想的脚本语言和开发环境，而ABAQUS软件也为用户预留了Python脚本的接口，以便于对软件进行二次开发，这为黏弹性边界和地震动输入问题的研究提供了新的思路。

本文主要对ABAQUS的Python接口进行了二次开发，并编写了Python程序，在保证模拟精度的前提下实现了黏弹性边界和地震波的快速施加，大大提高了边界施加的效率，并利用该方法对浅埋隧洞的地震动响应进行了分析，为进一步研究地下工程的地震响应提供了新的思路。

1 黏弹性边界和地震波施加原理

1.1 黏弹性边界

在利用有限的空间模拟半无限土体时，局部人工边界条件的施加对于计算的准确性有着很大的影响，许多学者也在这方面进行了大量的研究，其中刘晶波等[12]研究的黏弹性人工边界由于具有比较高的精度[8]和实用性，因而被广泛应用于土—结构的有限元软件计算中。其通过在离散的边界网格结点上施加并联的弹簧和阻尼器来实现从有限域到无限域的外行波吸收，在三维分析中该边界的施加方式如图1所示。其中弹簧的刚度和阻尼器的阻尼由下式计算：

式中：[K1]、[K2]分别为切向（X、Y方向）需要设置弹簧的刚度;[C1]、[C2]分别为切向（X、Y方向）需要设置阻尼器的阻尼;[K3]为法向（Z向）需要设置的弹簧的刚度;[C3]为法向（Z向）需要设置的阻尼器的阻尼值;ρ为介质的密度;R为地下结构中心到边界点的距离;[Ai]为每一个网格结点上的附属面积之和;[cs]、[cp]、G分别为土体介质的剪切波速，压缩波速和剪切模量。

1.2 黏弹性边界下地震波的施加原理

地震波在入射到有限土体中以后，遇到结构或者土体表面会发生散射，所以地震的总波场一般可以分解为自由波场和散射波场[19]。这其中的散射波场由人工边界吸收，而入射波则可以通过作用于人工边界上的等效荷载的方式施加上去。除此之外，还需要根据不同边界面上的波场特点，采用不同的波场分解方案。在底部边界面上可以将总场分解为边界入射场和边界外行场，而在侧面边界处可以将总波场分解为自由场和散射场[12]。其中边界自由场和入射场可以根据连续介质力学模型通过解析计算得到，而边界的外行场或者散射场则可以利用边界条件进行模拟。其中将散射场或者边界外行场作为未知量，并将其用总场与自由場或者入射场的差值表示，就可以得到含外源作用的人工边界面上任一点l在i方向的有限元运动方程[20]

式中：Kli、Cli为边界面上l结点的法向和切向的参数，可由式（1）～（4）计算;fli为入射地震波作用下l结点在i方向的等效结点应力，由式（6）表示

式中，[KliuRli、CliuRli、σRli]分别为克服弹簧、阻尼器和土体介质所需要克服的抗力。

在模型的底面边界垂直入射地震波，设P波和S波的位移时程为[up（t）]和[us（t）]，利用波动的传播规律和波场的应力状态，通过式（6）解析计算边界面上每个点的等效地震荷载[fli（t）]。由此即可以得到任一结点的计算式[20]。

2 ABAQUS-Python的二次开发

2.1 ABAQUS-Python的脚本接口

ABAQUS中的脚本接口能够非常高效地实现ABAQUS/CAE中的所有功能，二者的通信关系如图2所示。当Python程序编写好以后，所有的命令都需要通过Python的特定解释器进行翻译之后才能够进入到ABAQUS/CAE的内核中执行，在执行的过程中会生成记录文件，即扩展名为.rpy的文件。这些命令在进入内核中以后将会被转化为INP文件，最后程序将INP文件输入到求解器中进行求解。

2.2 二次开发主要算法

本次开发使用的主要软件为Python IDLE，首先需要导入ABAQUS相关函数模块，程序会自动读取在ABAQUS中设置好的材料参数，如弹性模量、泊松比、密度等;通过下式可以求得介质的纵向和剪切波速

式中：λ、μ为第一拉密常数和第二拉密常数;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.

[责任编辑 杨 屹]