第一作者黄景琦男，博士生，1988年生

通信作者杜修力男，教授,博士生导师，1963年生

近场数值波动分析中地震波输入的一种简化方法

黄景琦，杜修力，赵密，金浏

(北京工业大学城市与工程安全减灾教育部重点实验室,北京100124)

摘要：在显式有限元方法结合黏弹性人工边界的时域波动方法的基础上，建立了地震波垂直输入时的一种简化输入方法。将近场有限元模型沿高度方向进行分层，将地震动的入射运动转化为作用于人工边界底面及分层后每层侧面上的均布力，以实现地震动的输入。与以等效节点力的方式实现的地震动输入相比，施加均布力的方式简化了地震波输入的前处理工作，且又能保证与等效节点力方式相同的精度。自由场数值算例表明：当该方法中的分层高度与波动有限元网格离散要求的最大尺寸相等时，该方法与等效节点力方法具有相同的精度；当局部区域按网格离散要求的最大尺寸进行分层而在其它区域放大分层高度时，局部区域上的近场波动响应仍可保证具有相当高的精度。另外，某隧道结构地震响应算例的计算结果同样说明了该方法的有效性。

关键词：波动输入；地下结构；黏弹性人工边界；时域显式有限元

收稿日期：2013-11-29修改稿收到日期：2014-02-11

中图分类号:TU45文献标志码：A

A simplified seismic wave input method for near-field wave analysis

HUANGJing-qi,DUXiu-li,ZHAOMi,JINLiu(MOE Key Laboratory of Urban Security and Disaster Engineering, Beijing University of Technology, Beijing 100124, China)

Abstract:On the basis of the time-domain wave method combined with the explicit finite element method and the viscous-elastic artificial boundary condition, a simplified seismic wave input method for near-field wave analysis was proposed. With this approach, the near-field model was stratified in the height direction, and then the seismic motion was converted into stresses acted on the viscous-elastic artificial boundary to realize seismic wave input. Compared with the equivalent node force input method, the new method simplified the antecedent processing with an acceptable accuracy. The numerical results obtained from a free ground model indicated that when the stratified size equals the maximum size of the element mesh in dynamical analysis, the accuracy of the simplified method is the same as that of the equivalent node method; when the local parts of the near-field model were stratified with the maximum element mesn size of dynamical analysis while the other parts were stratified with a lager size, the dynamic responses of the local parts have a higher accuracy. Moreover, the dynamic response computation results of a rock tunnel demonstrated the effectiveness of the proposed simplified approach.

Key words:seismic wave input method; underground structure; viscous-elastic artificial boundary condition; explicit finite element method in time domain

许多复杂的工程问题，如高坝、隧道等结构的地震响应分析，均可归结为近场非线性波动问题。这类问题的特点是需要考虑多种复杂的影响因素，计算过程复杂，计算工作量巨大。对于近场波动问题，它涉及到无限域模型有限化处理带来的人工边界及人工边界上的地震动输入问题[1]。目前应用较多的人工边界有黏性边界[2]、透射边界[3]、黏弹性边界[4]等。其中黏弹性人工边界既能吸收远场的散射波又能模拟无限地基的弹性恢复性能，物理意义明确，便于在通用软件中实现，从而得到了广泛的应用。当应用人工边界条件进行结构的地震反应分析时，地震波的输入与所采用的人工边界条件相关联。Joyner等[5]对于一维模型采用将入射运动转化为作用于粘性人工边界上的等效荷载的方法，成功地解决了地震动的波动输入问题；刘晶波等[6]基于此种将入射运动转化为人工边界节点上的等效节点力的思想，建立了一种新的、适用于黏弹性人工边界的波动输入方法，并进行了二维模型的算例验证；王振宇等[7]建立了考虑非线性效应的无限成层地基波动模拟的时域人工边界条件及相应的波动输入方法；杜修力等[8]将一种新的黏弹性人工边界结合显式有限元的时域波动求解方法应用于拱坝-地基系统的地震反应分析中，取得了很好的模拟效果；另外，其他科研工作者[9-11]也在地震动相关的近场波动分析研究中开展了值得重视的工作。上述工作中的地震动输入往往是将人工边界处的入射地震动运动转化为作用于人工边界上的等效荷载，该荷载由抵抗人工边界条件和内域介质两部分构成，从而实现地震动的输入。但在有限元模型中，由于人工边界面上各个人工边界节点的空间位置以及边界节点影响面积的不同，使得不同边界节点的等效节点力时程也不相同。应用通用有限元软件实现波动输入时，在软件中单独对每一个节点施加等效节点力并不容易实现，尤其对于大型结构，单元与节点数目较多。通常的做法是将模型节点和单元信息导出，通过第三方软件进行编程批量的计算等效节点力，然后再导回有限元软件中计算，这使得波动输入的前处理工作较大，也较复杂，不易于波动输入方法的实现及工程应用。

本文在显式有限元法结合黏弹性人工边界的时域波动方法的基础上，建立了地震波垂直波动输入的简化输入方法。在该方法中，首先将有限元模型沿高度方向进行分层，然后将地震动的入射运动转化为作用于人工边界底面及分层后每层侧面上的均布力，由通用有限元程序来计算每个节点上的等效集中力。与以等效节点力形式实现的波动输入相比，本文的方法只要在人工边界面上分层施加均布力即可，不需在有限元软件外另外单独进行编程批量求解等效节点力，可在通用有限元软件中直接实现地震动输入，从而简化了地震动的波动输入在通用有限元软件中的实现。

1以等效节点力方式实现的地震动波动输入

在显式有限元法结合黏弹性边界的时域整体分析方法中，人工边界面l节点i方向的集中质量有限元总波场运动方程为[1]

(1)

将人工边界处的总波场分解为自由场(上标f表示)和散射场(上标S表示)，人工边界节点l方向i的总位移和总作用应力可以分别写为

(2)

(3)

散射场采用黏弹性人工边界条件模拟[12]，人工边界节点l方向i的应力-运动关系可以写为

(4)

式中：Kli和Cli为节点l方向i的人工边界参数，可参照文献[4]、[12]选取。

将式(2)的散射场位移及其时间导数代入式(4)，进一步将式(4)代入式(3)，最终将式(3)代入式(1)，整理得到考虑无限域辐射阻尼和地震波输入条件下人工边界节点的集中质量有限元运动方程

(5)

与荷载为零的非人工边界节点有限元方程相比，式(5)等号左边增加的两项是为了模拟无限域辐射阻尼而施加的黏弹性人工边界条件，物理上相当于在节点l方向i施加一个另一端固定的并联弹簧-阻尼器单元，如图1所示。

图1　三维黏弹性人工边界条件 Fig.1 Three-dimensionalviscous-spring artificial boundary conditions

式(5)等号右边表示地震荷载，即通过自由波场得到的人工边界节点l方向i处的自由场所对应的等效荷载，其中前两项表示产生自由场反应所需抵抗人工边界物理元件的节点力，第三项表示产生自由场反应所需抵抗近场介质的节点力。

对于以等效节点力形式实现的地震动波动输入而言，输入的关键问题是由已知入射平面波时程确定人工边界节点处的自由场位移时程和速度时程，从而计算边界节点处的等效节点力时程，然后将等效节点力加到模型的人工边界节点上，以此实现地震动的输入。但对于大型结构而言，有限元模型的节点数目较多，需要通过在第三方软件中另外编写程序批量求解等效节点力，这使得波动输入的前处理工作较大也较复杂，不便于波动输入方法的实现及工程应用。

2地震动波动输入的简化方法

在地震动的波动输入方法中，人工边界面上需要施加产生自由场反应所需抵抗人工边界物理元件的应力及产生自由场反应所需抵抗近场介质的应力。因此在截断的人工边界面上，地震动的作用形式如图2(a)所示，人工边界面上的不同点具有不同的应力时程。以等效节点力形式实现的波动输入实际上是忽略人工边界节点影响区域内人工边界应力的差异，将区域内的人工边界应力均匀的等效为集中力施加在人工边界节点上(图2(b))。同样，也可将实际作用在人工边界面上一定区域内的非均布应力等效为均布力施加在人工边界面上(图2(c))，模型计算时，由有限元程序将施加的分布力转化为节点上的集中力，这样就可以省去对此区域内的每个节点上的集中力进行人工求解，从而在建模的前处理工作上进行了简化。

图2　人工边界面上应力分布及其等效化方式 Fig.2 The stress distribution on the viscous-spring artificial boundary surface and the equivalent method

图3　振动方向简化的波动输入方法 Fig.3The simplified approximate seismic wave input method

在有限元软件的具体实现中，如图3中建立的有限元区域，首先将有限元模型沿高度方向进行分层，然后将求解的均布力施加在相应的边界面上便可实现地震动的波动输入。

2.1垂直入射地震波的均布力求解

对于垂直入射的S波，取图4的盒子状人工边界面，振动方向与X方向平行。对于底边界面的输入，由于模型底边界点的运动时程沿整个底面都相同，所以底边界的输入为无任何均匀等效的精确输入，施加均布力为

(6)

对于模型侧边界，由于侧边界上进行高度分层处理，对于每层侧边界都需要施加均布力。对于某层的左右侧面，施加均布力为

(7)

对于某层的前后侧面，施加均布力为

(8)

图4　简化波动输入方法中的均布力施加方法 Fig.4 Applied method for stress insimplified approximate input method

3分层高度与复合分层方法讨论

图5　自由场模型 Fig.5 Free ground model

在上述简化输入方法中，由于侧边界上沿高度进行分层输入，假定每层侧面上地震动运动时程相同，即忽略了在每层侧边上波动延迟效应，因此分层的高度尺寸会影响输入的精度。在波动有限元分析方法中，为了能反映地震波传播过程中的波形特性，有限元网格离散的尺寸要满足不大于地震波波长的1/8～1/10的要求，过大的网格会产生波形失真的问题。对于应力型人工边界面而言，地震波输入是将入射运动转化为人工边界面上的应力进行输入，如图2(a)，为了能够反应入射波的波形特性，在简化输入方法中，分层的高度尺寸也应满足不大于地震波波长的1/8～1/10的要求。等效节点力法是分层高度等于有限单元尺度时的本文方法特例，当分层高度包括几个单元尺度时，本文方法的精度需要讨论。

下面将以自由场模型作进一步讨论。建立的自由

场模型如图5所示，模型高度200 m，长、宽各100 m，材料的弹模为6×109MPa，泊松比0.3，密度2 308 kg/m3，剪切波速约为1 000 m/s。输入地震波为经典Kobe波，图6为Kobe波的加速度时程曲线。

3.1分层高度讨论

分别以100 m、50 m及20 m均匀的沿高度方向对有限元模型进行分层，其中20 m尺度满足小于地震波波长的1/10要求。

图7为Mises应力峰值随深度的变化，从中可以看出100 m尺度的计算结果与等效节点力方法差距较大。随着分层高度的减小，计算结果逐渐向等效节点力方法的结果靠近，当分层高度为20 m时，简化输入方法与等效节点方法的计算结果完全重合。图8为水平位移峰值随深度的变化，同样随着分层高度的减小，简化输入方法与等效节点力方法的计算结果逐渐靠近，当分层高度为20 m时，简化输入方式与等效节点力方法的计算结果完全相同。

图6 输入的Kobe波加速度时程曲线Fig.6AccelerationtimehistoryofimputedKobewave图7 Misses应力峰值随深度变化Fig.7DistributionofMissesstress图8 水平位移峰值随深度变化Fig.7Distributionofhorizontaldisplacementamplitude

图7和图8的分析结果说明：采用本文建立的简化输入方法时，计算的精度取决于分层的高度。分层高度越小，计算的结果将越为精确。当分层高度与波动有限元网格离散要求的最大尺寸相等时，本文方法与等效节点力方法具有相同的精度。

3.2复合分层方法讨论

上面讨论了地震动波动输入简化方法的分层高度的标准，即当分层高度与波动有限元网格离散要求的最大尺寸相等时，简化方法与等效节点力的方法具有相同的精度。按此原则进行分层时，对于埋深较浅，模型竖向尺度不大时，分层较少。当埋深较深时，模型的竖向尺度相对与地震波的波长较大时，则划分的层数将会增多。对于我们关心的模型区域，比如工程隧道所在的深度范围附近，可以选择按网格离散的最大尺度进行严格分层，对于远离工程结构的其它高度区域可选择将划分高度放大，以减小工作量。为了验证上述复合分层方法的可行性，将有限元模型在50 m～150 m深度范围内进行严格分层，即每20 m划分一层，其余部位按50 m进行放大分层。

图9　复合分层时Mises应力峰值随深度变化 Fig.9 Distribution of Misses stress

图9与图10分别为应力计算结果与位移计算结果随深度的变化情况。从图中可以看出，虽然其它部位进行了按50 m的放大分层处理，但严格分层区域内的自由场地震响应仍然与等效节点力方式的计算结果较为吻合，说明进行较大尺寸的工程计算时，局部区域可按离散要求的最大尺寸值进行严格分层而在其它区域放大分层高度，以此进一步的减小地震输入的工作量。

图10　复合分层时水平位移峰值随深度变化 Fig.10 Distribution of horizontal displacement amplitude

4工程算例

4.1有限元模型

以某岩体隧道洞身段为研究对象，如图11所示，横向(x方向)取200 m，竖向(y向)取200 m，纵向(z向)取200 m，隧道埋深100 m。衬砌厚度取50 cm。图11中为衬砌结构中布置的A、B、C、D四个监测点。本构模型中，岩体、衬砌均采用线弹性模型，材料参数见表1。在有限元模型的4个侧面和底面施加黏弹性人工边界。

图11　计算模型及监测点分布 Fig.11 Computational model and location of monitoring points

材料密度/(kg·m-3)弹性模量/GPa泊松比剪切波速/(m·s-1)围岩230860.31000衬砌2500300.2—

4.2地震波输入

图12　有限元模型分层示意图 Fig.12 Hierarchical information of model

在模型底部垂直输入S波，输入地震波为经典Kobe波(图6)。地震动输入分别按本文建立的简化输入方法和等效节点力方法进行输入，图12为有限元模型进行的分层情况。在隧道周围按分层高度小于地震波波长1/10的尺寸进行分层，分层高度为20 m，其余部分进行了放大分层,分层高度为70 m。

4.3结果分析

图12为监测点A、B、C和D点加速度时程曲线。从图12中可以看出简化输入方法的隧道监测点加速度反应与等效节点力方法的加速度反应基本相同。图13和14为监测点A、B、C和D点Mises应力时程曲线，如同加速度反应曲线的对比，两种输入方法下的隧道Mises应力反应也基本相同。另外，下文也给出了两种输入方法下的隧道监测点第一主应力峰值与第三主应力峰值，见表2。表2同样能对比出在两种输入方式下，岩体隧道的地震动反应较为吻合。

图13　各监测点加速度曲线 Fig.13 Acceleration time history of A,B,C,D points under seismic wave input

图14　各监测点Mises应力曲线 Fig.14 Misesstress of A,B,C,D points under seismic wave input

通过上述岩体隧道地震响应模拟的对比分析，本文建立的简化输入方法与等效节点力方法相比，在岩体隧道地震反应数值计算中具有较好的吻合度，说明了简化输入方法的有效性与准确性。

表2　各监测点主应力幅值(MPa)

附注：法一为等效节点力方法；法二为本文简化方法

5结论

本文在显式有限元方法结合黏弹性人工边界的时域波动方法的基础上，建立了近场数值波动分析中地震动输入的一种简化方法。与等效节点力的方法相比，本文方法简化了地震动输入的前处理工作，且又能保证与等效节点力方式相同精度。另外文中探讨了简化输入方法时分层高度的标准与复合分层的方法，并对比了等效节点力方法与本文简化方法下的隧道地震响应结果，得到以下结论：

(1)本文方法的计算精度取决于分层的高度。分层高度越小，计算的结果将越为精确；当分层高度与波动有限元网格离散要求的最大尺寸相等时，本文方法与等效节点力方法具有相同的精度；

(2)在进行较大尺寸的工程计算时，局部区域按可按网格离散要求的最大尺寸值进行严格分层而在其它区域放大分层高度，计算精度依然可以获得保证;

(3)简化输入方法与等效节点力方法相比，在岩体隧道地震反应数值计算中具有较好的吻合度。

参考文献

[1]杜修力. 工程波动理论与方法[M]. 北京:科学出版社，2009: 15-21.

[2]Lysmer J, Kuhlemeyer R L. Finite dynamic model forinfinite media[J].Journal of the Engineering Mechanics Division ASCE 95(EM4),1969:869-877.

[3]廖振鹏. 工程波动理论导论[M]. 北京: 科学出版社, 2002.

[4]Deeks A J, Randolph M F. Axisymmetric time-domain transmitting boundaries[J]. Journal of Engineering Mechanics, ASCE, 1994, 120(1): 25-42.

[5]Joyner W B,Chen A T F. Calculation of nonlinear ground response in earthquake[J]. Bull. Seism. Soc. Amer., 1975, 65(5): 1315-1336.

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

LIU Jing-bo, LU Yan-dong. A direct method for analysis of dynamic soil-structure interaction[J]. China Civil Engineering Journal, 1998, 31(3): 55-64.

[7]王振宇, 刘晶波. 成层地基非线性波动问题人工边界与波动输入研究[J]. 岩石力学与工程学报, 2004, 23(7): 1169-1173.

WANG Zhen-yu, LIU Jing-bo. Study of wave motion input and artificial boundary for problem of nonlinear wave motion in layered soil[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(7): 1169-1173.

[8]杜修力, 赵密. 基于黏弹性边界的拱坝地震反应分析方法[J]. 水利学报, 2006, 37(9): 1063-1069.

DU Xiu-li, ZHAO Mi. Analysis method for seismic response of arch dams in time domain based on viscous-spring artificial boundary conditions[J]. Shuli Xuebao, 2006, 37(9): 1063-1069.

[9]张顺福, 王海波, 杨会臣, 等效黏性边界单元及黏性边界波动输入方法[J]. 水利学报, 2008, 39(10): 1248-1255.

ZHANG Shun-fu, WANG Hai-bo, YANG Hui-chen. Equivalent viscous boundary elements and the method of wave input for viscous boundary[J].Shuli Xuebao, 2008, 39(10): 1248-1255.

[10]张波, 李术才, 杨学英, 等. 三维黏弹性人工边界地震波动输入方法研究[J]. 岩土力学, 2009, 30(3): 774-778.

ZHANG Bo, LI Shu-cai, YANG Xue-ying, et al. Research on seismic wave input with three-dimensional viscoelastic artificial boundary[J]. Rock and Soil Mechanics, 2009, 30(3): 774-778.

[11]赵武胜, 陈卫忠, 黄胜, 等. 地下结构抗震分析中地震动输入方法研究[J]. 土木工程学报, 2010, 43(S): 554-559.

ZHAO Wu-sheng, CHEN Wei-zhong, HUANG Sheng, et al. Research on method of seismic motion input in a seismic analysis for underground structure[J]. China Civil Engineering Journal, 2010, 43(S): 554-559.

[12]刘晶波, 王振宇, 杜修力, 等. 波动问题中的三维时域粘弹性人工边界[J]. 工程力学, 2005, 22(6): 46-51.

LIU Jing-bo, WANG Zhen-yu, DU Xiu-li, et al.Three-dimensional viscos-elastic artificial boundaries in time domain for wave motion problems [J]. Engineering Mechanics, 2005, 22(6): 46-51.