阳 栋，王志亮，2

（1.同济大学土木工程学院，上海200092；2.同济大学岩土与地下工程教育部重点实验室，上海200092）

为了对无限地基进行模拟，必须截取有限范围的地基，然后设置合理的无反射人工边界模拟无限地基辐射阻尼效应.局部人工边界由于其时空解耦性、计算耗时少，在地基辐射阻尼效应模拟中获得广泛应用.Lysmer等［1］基于一维波动理论提出黏性边界，在许多波动问题中得到了广泛应用，也是我国核电站抗震设计规范建议使用的人工边界之一［2］.廖振鹏［3］基于透射理论提出了透射边界，但是计算中可能引起高频失稳，而且编程比较复杂.Deeks等［4］、刘晶波等［2］基于柱面波动方程建立了可以考虑半无限介质弹性恢复能力的二维黏弹性人工边界，能模拟人工边界外无限介质弹性恢复性能，有良好的频率稳定性，但是黏弹性边界不方便考虑初始重力作用等.

无限元概念最早由Ungless［5］提出，后经过Bettess［6］，Beer等［7］和Zienkiewicz［8］的改进和发展，国内学者ZHANG等［9］、葛修润等［10］也对无限元的研究开展了不少工作.燕柳斌［11］利用无限元法模拟半空间弹性地基和重力坝地基；黄胜等［12］用无限元对隧道地震响应进行模拟，但是没有考虑自由波场的作用.相对于黏性边界，无限元不仅能模拟远场地基对地震波能量的吸收作用，还能够准确模拟无穷远处位移为零的边界条件，采用无限元能够大大减少单元数量，节省计算费用［13］.

本文采用有限元与无限元耦合法模拟地基动力响应，为了考虑人工截断边界对波传播影响，应用波场分离技术，在有限元和无限元分界处施加等效节点荷载，以实现地震输入，并对边界和地震输入的有效性进行验证.同时，把该法应用于核电站地震响应计算，对自由场地和考虑土与结构相互作用（SSI）时地震响应进行对比分析.

1 无限元原理与控制方程

ABAQUS中提供了一阶和二阶无限单元，其在连续实体静力分析中能提供刚度，在动态分析中提供静止边界［14］.在进行数值模拟时，由于主要分析近场区域结果，故需用有限元模拟，而远场地基可假定为线性的，用无限元模拟.

假定波沿z轴传播，则一维波动的平衡方程

对于均匀弹性介质平面波

波动方程uz＝f1（z－ct）代表沿z轴正向传播的波，而uz＝f2（z＋ct）代表沿z轴负向传播的波，假设有限元与无限元的分界面为z＝L，因为问题是线性的，所以位移可以叠加，即

当波传播到界面上时，界面上的应力为

式中：f′1，f′2分别为函数f1，f2对z的导数.

如果在有限元边界上设置阻尼器平衡波动力，其阻尼系数为C，并使下式成立：

为了达到不反射的效果，则反射波引起的位移和应力为零，即f2＝0，f′2＝0，由此解得

通过设置参数，就可以避免有限元与无限元边界处波的能量反射.

2 波场分离方法和地震输入

对于内源波动问题，如动力机器、地铁等引起的振动问题，用无限元模拟远场时，可以直接在近场有限元区域内施加振动作用，当波传播到有限元边界时，自动透射进入无限元区域.而对于外源波动问题，如地震作用，存在地震波散射现象，需要在人工边界处进行波场分解.根据不同人工边界上的波场特点，各边界采用不同的分解方案，即：底边界上将总场分解为边界入射场和边界外行场，侧边界上将总场分解为自由场和散射场，边界入射场或自由场可由连续介质力学模型解析计算得到，边界外行场或散射场由人工边界条件模拟并由离散模型依据数值分析方法获得.这种不同边界分区分解的方法，在尽量避免计算自由场的前提下较好地解决了人工边界的平行波模拟问题［15］.

外源作用下，有限元边界任意节点运动方程为

式中：uli，，分别为l节点在i方向的位移、速度和加速度；Cl为节点l的阻尼系数；Kl为节点l的刚度系数；和分别为模拟边界自由场和散射场在l节点沿i方向施加的等效荷载；ml为节点l质量.

在模拟散射场时，有限元边界节点上施加的应力由无限元提供，无需在边界节点上附加额外的荷载；模拟自由场时，在有限元边界节点上需要施加包括克服无限单元阻尼所需的应力和自由场在有限元边界处的应力.利用式（7）进行波动输入的前提是首先确定自由场，而自由场一般采用解析法或数值法确定.杜修力等［16］给出了避免计算自由场前提下的黏弹性边界等效荷载计算方法，在此基础上去掉弹簧力作用，即可得到无限元边界等效荷载的计算公式.下面以平面问题S波从有限元底部边界垂直入射为例，其等效荷载计算式为在底边界上

在侧边界处

式中：Δt1，Δt2分别为l节点处入射S波和地表反射S波的时间延迟，即Δt1＝h／cS，Δt2＝（2 H－h）／cS，其中H为地表面到有限元区域底部边界的距离，h为l节点到有限元底边界的距离，cS为S波波速；Al为有限元边界节点l的影响面积；ClN，ClT为法向和切向阻尼系数，其值为ClN＝ρcS，ClT＝ρcP，cP为纵波波速.等效地震荷载下标代表节点号和分量方向，上标代表节点所在有限元边界界面的外法线方向，与坐标轴方向一致为正.

在ABAQUS计算时，先用FORTRAN90编制程序计算等效荷载，然后选取边界节点，通过调用子程序在边界节点上施加等效荷载.

3 边界及地震输入的可行性验证

为了验证无限元和等效荷载输入模拟自由场地震的有效性，本文通过一个简单的模型计算，验证半无限地基中波动传播规律.对于二维弹性半空间模型，从底部垂直入射一S波，入射波方程如式（13）所示，求二维弹性半平面内各点的位移反应

已知岩石地基的剪切模量G＝4.233 6GPa，泊松比ν＝0.25，密度ρ＝2 100kg·m－3，波速cS＝1 400m·s－1，cP＝2 425m·s－1.选取有限元的计算区域宽为762m，高381m，用平面应变四边形减缩积分单元剖分网格，单元尺寸为Δx＝Δy＝19.05 m，有限单元数为800个，节点数为861个.对于无限元区域，选取顶部中点（C点）作为辐射中心，延长C点与有限元边界节点的连线，并使无限单元的长度等于有限元边界节点到辐射中心的长度，建立无限元网格.总的计算时间为2.0s，取固定增量步长为0.005s，计算网格模型如图1a所示.

图2 A，B和C点水平位移时程曲线Fig.2 Time histories of horizontal displacement at nodes A，Band C

由一维波动理论可知，波分别在t＝0.136s和t＝0.272s时到达模型中部（B点）和顶部（C点）.由于顶部是自由表面，故顶部的最大振幅相比入射波放大2倍，波在顶部经自由表面反射后向下传播，在t＝0.408s时反射波到达中部，在t＝0.545s时反射波到达底部（A点），在t＝0.945s时波完全经过底部有限元边界并透射进入无限元区域，此时有限元区域中所有节点的位移均变为零.从图2可看出，模拟结果与理论值吻合较好，这也说明上述采用无限元结合等效荷载输入来模拟波在半无限地基中传播的方法是正确的.

4 工程应用实例

4.1 计算模型和地震输入

设有一核电站建在上述场地上（见图1b），把核电站简化为一矩形结构，尺寸114.3m×38.1m，单元尺寸取为19.05m.混凝土结构密度为2 643 kg·m－3，弹性模量为31.027GPa，泊松比为0.15.在ABAQUS中采用隐式积分方法求解.

本例以Koyna地震的S波［14］作为输入波，对自由场及考虑SSI效应下的地震响应进行分析.由于实际记录的是地表地震波，其中已包含地基的放大效应，在应用到底部基岩时，可取其幅值的一半.通常记录的地震波数据中混杂了各种干扰，直接对加速度时程进行二次积分得到的位移时程会产生漂移失真，因此要对加速度时程进行校准和滤波，达到消除噪声的目的.此处用地震分析软件Seismosignal对原始的加速度时程数据进行处理，图3是经校准和滤波后的加速度时程曲线、位移时程曲线及加速度反应谱.

图3 输入波的加速度时程、位移时程曲线及加速度反应谱Fig.3 Acceleration and displacement time histories and acceleration response spectrum of input wave

4.2 自由场分析

图4 无限元自由场地震加速度反应谱Fig.4 Acceleration response spectrum of free－field modeled with infinite element

此处以图1a的自由场为例，分别进行施加重力与不施加重力的地震模拟，前者在地震输入前单独建立一个时间步长为1s的重力加载步，总的分析时间为12.0s.图4a，b分别为无重力与考虑重力作用下，自由场中A，B，C三点的加速度反应谱，图5a为自由场顶部C点的位移时程曲线.从图4和5a可看出，重力对地震加速度和位移响应都无影响，这是因为重力不影响水平方向力的平衡.在周期为0.25s和0.70s时，场地顶部的加速度反应谱放大超过2倍，而周期为1.00s时，场地对加速度反应谱的放大作用最不明显，由于I类核电场地的特征周期约为0.25s，故0.25s处的加速度反应谱值最大.图4a中A点曲线为有限元底部地震波入射点处的加速度反应谱，由于其中包括反射波的影响，故其与图3c中入射波的加速度反应谱有所区别.C点位移是入射波和反射波共同作用的结果，对比图5a和图3b，可以看出自由场位移响应值为入射位移波的2倍.通过与一维地震分析程序SHAKE91［17］的计算结果作对比（如图4a、图5a所示），也可以看出本文采用的边界和地震动输入是有效的.

4.3 SSI效应分析

对图1b中模型输入上述Koyna波进行计算，并选取节点A1～G1对结果进行分析.图6为考虑SSI效应后的地震加速度反应谱，对比图6a和4a可以看出，考虑SSI后，结构底部处C点的加速度反应谱最大值从1.92g下降到0.94g（其中g为重力加速度），且在周期为0.25s附近，加速度反应谱下降，而在周期为0.30s到4.00s范围内，加速度反应谱增大，即SSI效应使短周期段的加速度反应谱值降低，而长周期段的加速度反应谱值增加.对于结构底部标高的地基节点C1，E1，F1，G1，从图6b可看出，在结构与有限元边界的中间点F1处，加速度反应谱最大值为2.36g，大于自由场中的加速度反应谱最大值1.92g，即在S波的激励下，SSI效应使结构附近地基的加速度反应谱降低，而远离结构处的加速度反应谱峰值相比自由场反而略有升高.图5b为自由场顶部C点和考虑SSI效应后相同位置处C1点的位移时程曲线图，可以看出S波作用下，SSI效应对地基的水平位移无影响.图7为地表C点在自由场和考虑SSI效应时最大主应力和最小主应力的时程曲线，由于结构底部和地基之间将产生较大的摩擦力，从图中容易发现SSI效应使得地基中C点的最大主应力和最小主应力的幅值均增大2倍之多.

5 结论

（1）采用近场地基由有限元模拟，远场土体由无限元模拟，并使用本文所提出的等效荷载计算公式计算等效节点力，在有限元和无限元交界处输入地震作用，此法能够有效仿真地震波在半无限地基中的传播过程.

图7 自由场C点及考虑SSI效应时C1点的主应力时程曲线Fig.7 Time histories of principal stress at point Cin free－field and point C1with a consideration of SSI effect

（2）对于岩石地基，当地震波从底部基岩传播到地表时，加速度反应谱峰值和位移振幅都放大2倍，通过自由场地施加重力和无重力的对比分析，可以看出S波作用下，重力对地震作用时地基的加速度反应谱和水平位移响应无明显影响.

（3）在S波作用下，SSI效应使短周期段的加速度反应谱值降低，而长周期段的加速度反应谱值增加；结构附近地基的加速度反应谱峰值相比自由场降低，而远离结构处加速度反应谱峰值相比自由场略有提高；SSI作用对地基的水平位移无影响，但使结构底部地基的最大主应力和最小主应力幅值都有增大趋势.

［1］ Lysmer J，Kulemeyer R L.Finite dynamic model for infinite media［J］.Journal of Engineering Mechanics Division，ASCE，1969，95（4）：759.

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

LIU Jingbo，GU Yin，DU Yixin.Consistent viscous－spring artificial boundaries and viscous－spring boundary elements［J］.Journal of Geotechnical Engineering，2006，28（9）：1070.

［3］ 廖振鹏.工程波动理论导论［M］.北京：科学出版社，2002.

LIAO Zhenpeng.Introduction to wave motion theories in engineering［M］.Beijing：Science Press，2002.

［4］ Deeks A J，Randolph M F.Axisymmetric time－domain transmitting boundaries［J］.Journal of Engineering Mechanics Division，ASCE，1994，120（1）：25.

［5］ Ungless R F.An infinite finite element［D］.Vancouver：University of British Columbia，1973.

［6］ Bettess P.More on infinite element［J］.International Journal for Numerical Methods in Engineering，1980，15（11）：1613.

［7］ Beer G，Meek J L.Infinite domain element［J］.International Journal for Numerical Methods in Engineering，1981，17（3）：43.

［8］ Zienkiewicz O C.A novel boundary infinite element［J］.International Journal for Numerical Methods in Engineering，1983，19（3）：393.

［9］ ZHANG Chuhan，ZHAO Chongbin.Coupling method of finite and infinite elements for strip foundation wave problems［J］.Earthquake Engineering and Structural Dynamics，1987，15（7）：839.

［10］ 葛修润，谷先荣，丰定祥.三维无限元和节理无界元［J］.岩土工程学报，1986，8（5）：9.

GE Xiurun，GU Xianrong，FENG Dingxiang.Threedimensional infinite domain elements and joint infinite domain elements［J］.Journal of Geotechnical Engineering，1986，8（5）：9.

［11］ 燕柳斌.用三维映射无限元模拟重力坝地基［J］.水利学报，1991，2（10）：7.

YAN Liubin.Modeling of gravity dam foundation using three dimensional mapped infinite element［J］.Journal of Hydraulic Engineering，1991，2（10）：7.

［12］ 黄胜，陈卫忠，杨建平，等.地下工程地震动力响应及抗震研究［J］.岩石力学与工程学报，2009，28（3）：483.

HUANG Sheng，CHEN Weizhong，YANG Jianping，et al.Research on earthquake－induced dynamic response and aseismic measures for underground engineering［J］.Chinese Journal of Rock Mechanics and Engineering，2009，28（3）：483.

［13］ 姜袁，陈灯红，姚艳华.有限元－无限元法在重力坝应力分析中的应用［J］.武汉大学学报：工学版，2009，42（3）：322.

JIANG Yuan，CHEN Denghong，YAO Yanhua.Application of finite element and infinite element method to stress analysis of gravity dam［J］.Engineering Journal of Wuhan University，2009，42（3）：322.

［14］ SIMULIA Inc.ABAQUS theory manual［EB／OL］.［2011－04－20］.USA：https：／／www.sharcnet.ca／Software／Abaqus610／Documentation／ddoc／v6.10／books／stm／default.htm.

［15］ 赵建锋，杜修力，韩强，等.外源波动问题数值模拟的一种实现方式［J］.工程力学，2007，24（4）：52.

ZHAO Jianfeng，DU Xiu1i，HAN Qiang，et al.An approach to numerical simulation for external source wave motion［J］.Engineering Mechanics，2007，24（4）：52.

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

DU Xiuli，ZHAO Mi.Analysis method for seismic response of arch dams in time domain based on viscous－spring artificial boundary condition［J］.Journal of Hydraulic Engineering，2006，37（9）：1063.

［17］ Idriss I M，Sun J I.SHAKE 91：a computer program for conducting equivalent linear seismic response analyses of horizontally layered soils deposits［EB／OL］.［2011－04－20］.USA：http：／／nees.ucdavis.edu／publications／Shak91＿manual.pdf.