基于ABAQUS的水平成层场地地震反应无限元分析

李琦

(泰山学院建筑与机械工程系，山东泰安 271021)

1 引言

对于地下结构抗震设计的拟静力计算方法，无论是反应位移法还是反应加速度法，都要先进行一维土层的地震反应分析.传统的方法都是采用动力有限元法进行分析，它能够出来土体材料的非线性、成层性等问题.但是，由于地下结构处于半无限的空间中，如果在进行动力分析时，对选取的计算区域设置简单的截断边界是不合适的.通常采用耦合的方法在设置合理的边界条件，可以在截断处采用边界元或者设置粘滞边界［1］.但是在有限元计算中，这种边界的设置和使用过于复杂，工程设计人员难以掌握，会带来较高的计算成本.

为实现保证计算精度的前提下尽可能简化计算过程和提高计算效率，本文结合文献［2］中的一维土体的无限元计算方法，基于ABAQUS有限元分析平台中的无限元模型来模拟岩体的半无限边界条件.这种方法的基本原理和过程与有限元分析十分接近，可以实现简便、快速、准确的一维岩体的动力反应计算［2］.

2 ABAQUS中的无限单元

ABAQUS中可以提供轴对称平面无限元和三维无限元，其收敛阶数为一阶或二阶，其中平面无限元相应的单元编号如图2－1所示.

ABAQUS╲CAE并不支持直接定义无限单元，需要在input文件中加以定义.以四结点的无限元为例，在input文件中定义无限元，首先要确定虚拟点r0，然后从r0开始向需要定义的方向延伸.虚拟的点的选择可以根据解的先验信息进行，如对于半空间边界上左右的点荷载问题，可以选加载点作为虚拟点.然后，定义无限元与之耦合的有限元的结点r1.之后从r1开始沿着r0r1的方向延伸出去，进而定义r2.如图2－2所示，用*NCOPY，POLE命令可以方便地在某一无限方向上定义r2.在定义r2时要注意r1与r2的配置，r2到r0的距离一定要是r1到r0的两倍，并务必使无限元的两条无限边界的延长线在远处不能相交.

图2 －1 平面无限单元的编号

图2 －2 在某一无限方向上定义r2

3 民治地下洞室厂址区一维岩体地震动分析

3.1 模型参数及地震动输入

对于一维成层场地地震分析，本文结合民治地下洞室厂址的现场围岩测试数据，按照有关规程规范要求，结合建筑物布置和已建工程类比，确定场地共分三层，由此选取模型参数:岩体物理模型采用各向同性线弹性模型，参数见表3－1.

表3 －1 材料参数表［3］

有限元模型:对厂区岩体进行二维有限元离散，有限元计算范围为垂直于厂房轴线向取323.0米，高程方向自+1230.0米至山体表面，此处将山体表面处理为水平表面，其高程为+1650.0米，即切取自岩体表面向下420米深，宽度为323米的计算区域.有限元模型共分256个四结点四边形线性减缩积分平面应变单元(CPE4R)，一共289个结点.在有限元计算区域两侧设置34个4结点线性平面应变无限单元(CINPE4)与有限单元进行耦合，如图3－1.

地震波选用印度的Koyna重力坝地震时记录的水平向加速度时程谱，在ABAQUS的inp文件里用右图曲线坐标的点来给ABAQUS动力反应计算提供反应谱曲线，计算时间10秒，在底部基岩处输入地震波进行动力响应计算.

图3 －2 一维岩体自振一阶振型

图3 －1 一维场地地震分析计算简图

3.2 阻尼的确定

模型采用瑞利阻尼，令α=0，并设定模态阻尼中的阻尼比为ξ=0.3，则假设β只与场地围岩的第一阶临界阻尼和特征频率有关，对计算模型提取自振振型及固有频率，得到结果如图3－2，通过计算可得所取厂区岩体的一阶特征值为80.664，与之对应的一阶自振频率为ω1=8.982(RAD/S)，代入公式β= 2ξ/ω1，可以求得模型的瑞利阻尼系数为α=0和β=0.00668.

3.3 计算结果分析

应用ABAQUS计算出一维岩体在地震作用下的加速度分布，计算时间为10秒，每个分析步的时间为0.01秒，每十个分析步输出一次计算结果，图3－3列出岩体分别在0.2、0.4、0.6、0.8秒时的加速度分布.由图3－3可以看出，岩体在受到地震作用时水平加速度出现明显的成层分布，而入射到两侧的地震波都被无限元边界吸收，没有产生反射效应，这与实际情况是相符的.

图3 －3 岩体在地震过程中的水平加速度分布图(单位:m/s2)

图3 －4 结点151、411、563、617的水平加速度时程(time:s;A1:m/s2)

在进行计算时设置了对应历史变量的输出，下面选取地面、基岩以及洞室所在位置的几个关键结点，将其10秒内的加速度时程显示如图3－4.其中，结点151为模型地面处结点，最大加速度(绝对值，以下同)为0.325764m/s2，结点411为模型主厂房所在区域关键结点，最大加速度为1.14749m/s2，结点563为模型主变室所在区域关键结点，最大加速度为1.47996m/s2，结点617为模型基岩处结点，最大加速度为4.26931m/s2.

图3 －6 剪切应力沿高程的分布(t=5.5s)

图3 －5 相对位移沿高程的分布(t=5.5s)

图3 －7 水平加速度沿高程的分布(t=5.5s)

图3 －8 水平向速度沿高程的分布(t=5.5s)

场地岩体上下底位置的最大相当位移发生在t=5.5s时，此时由一维场地地震动分析得到的岩体各层上的相对位移、剪切应力、加速度和速度的数值沿高程的分布情况见图3－5～图3－8.并且，取各个结点加速度时程中的最大值，列出其沿高程方向的分布见图3－9.

将计算所得结果与文献［2］中:将土层沿深度详细划分，用SHAKE程序对自由土层进行计算，得到土层上下底位置发生最大相对位移时(t=4.59 s)一维自由土层的加速度、速度、相对位移和剪应力沿深度方向的分布数据;逐一进行比较，在相对应的深度范围内其分布情况基本吻合，各个分布数据的数值相差为相对位移偏小25%～30%，剪切应力偏小20%～33%，加速度偏小60%～70%，速度偏小75%左右.

鉴于文献［2］研究的是自由土体，而本文研究的是自由岩体，其弹性模量、泊松比以及自身密度都有较大差异，因此可以得出结论:本文中应用ABAQUS程序进行的一维岩体地震动分析所得到的数据是可靠的，可以作为进行进一步动力分析和抗震分析的参考数据，并为实际工程设计提供依据.

图3 －9 场地岩体一维地震动反应分析水平加速度峰值

4 小结

本文基于动力无限元分析理论，采用ABAQUS大型有限元软件进行对于弹性动本构岩体的一维地震动反应分析.应用ABAQUS中的无限元模型模拟半无限边界，使用软件中的瑞利阻尼理论模拟岩石对地震波的耗散，并采用多点输入地震加速度时程，实现了在强震情况下的一维岩体动力分析.进而同基于SHAKE程序进行的一维土体地震动反应分析结果作了对比，说明本文在ABAQUS软件中实现的一维岩体动力反应分析方法是可行并且具有合理性的.

［1］陈健云，胡志强，林皋.超大型地下洞室群的三维地震响应分析［J］.岩土工程学报，2001，23(4).

［2］谷拴成，朱彬，马德梅.地下结构抗震计算中拟静力法的地震荷载施加方法研究［J］.西安科技大学学报，2005，25(2).

［3］邵国建，苏静波.民治水电站地下厂房系统围岩稳定、变形、开挖次序及支护参数分析研究［R］.南京:河海大学工程力学系，2007.