ABAQUS岩土动力分析的静-动力统一人工边界研究

左得奇，蒋良潍，骆 达，杜美玲，葛学军

（1.西南交通大学土木工程学院，四川成都 610031；2.高速铁路线路工程教育部重点实验室，四川成都 610031；3.中铁二院工程集团有限责任公司，四川成都 610031）

引言

数值分析是岩土静、动力问题不可或缺的研究手段之一，人工边界的设置及动力作用输入方法将直接影响数值计算结果的合理性。因此，各类人工边界和地震动输入方法相继被提出。人工边界可分为全局人工边界和局部人工边界，Lysmer等［1］首先提出了以黏性理论为基础的局部人工边界，其概念清楚、应用方便，但精度较低，容易出现低频失稳问题；Deeks等［2］进一步发展提出了黏弹性人工边界，其引入弹簧单元以模拟地基恢复力，克服了黏性边界低频失稳的问题；刘晶波等［3］推导了三维形式的黏弹性人工边界并将其应用于地下土-结构的相互作用的计算中；廖振鹏等［4］提出了一种对一般无限域模型具有普遍适应性的人工透射边界，入射波由一维单侧波动的叠加构成，其具有实用性广、与有限元方法结合方便的优点，但也存在低频失稳问题；在全局人工边界方面，Ungless［5］首先提出了无限元思想；Kim［6］和Yun［7］利用无限元方法在频域和时域内研究了二维、三维层状土-结构的相互作用，其具有使用方便、前处理工作量小的优势；此外，动力作用输入方式与人工边界条件的设置密切相关，刘晶波等［8］提出了基于黏弹性人工边界的地震动输入方法；杜修力等［9］考虑场地自由场效应推导了下卧刚性基岩条件下的地震动输入方法；马笙杰等［10］则对黏弹性人工边界设置及地震动输入方法进行了详细对比分析；赵武胜等［11］通过分析人工边界相互作用推导了地震波垂直及斜入射条件下地震动输入方法，并以无限元为例进行了正确性验证，上述地震动输入方法多为基于动力人工边界推导得出，适用于动静分算的动力分量计算。

对于线弹性、小变形问题，可分别利用静、动力边界计算静、动力响应分量，再叠加得到完整解。由于土体是一种非线性材料，在力学作用下可能发生非线性响应甚至塑性变形，叠加原理将不再适用，为此研究者提出了静-动力统一人工边界。刘晶波等［12］基于弹性半无限空间中静力问题基本解对黏弹性人工边界进行节点约束刚度修正，提出了黏弹性静-动力统一人工边界，使得静、动力分析能够连续的进行；高峰等［13］通过对比研究发现，该静-动力统一人工边界应用于初始地应力生成误差较大，并提出了一种基于有限元静力边界节点支反力的黏弹性人工边界设置方法。

由上可知，目前人工边界的研究主要以黏弹性人工边界为主，因其在低频问题中计算稳定性好、效率较高，通过边界条件转化［13］、引入边界等效节点力［9-11］等方法，使得其能广泛应用于各类岩土动力问题，而其直接应用于静力计算时精度较低、边界设置较为繁琐，前处理工作量较大。无限元在理论上满足在无穷远处位移为零，波传至无穷远处衰减为零的客观性条件，部分耦合了静、动力计算，且边界设置较为简单、应用方便，大量研究将其直接应用于静力计算（如初始应力状态生成［14-16］）和动力计算（内、外源振动问题［14-16］），王飞等［17］基于ABAQUS无限元的静力计算结果具有可信性的前提，将文献［11］中的地震动输入方法引入无限元人工边界。但实际上，岩土动力计算中，第一步为进行初始应力状态生成即地应力平衡，平衡后位移近似为零，而应力状态的合理性往往被忽略，无限元应用于静力步骤的计算精度会直接影响初始应力状态，由于在体力荷载作用下，静力无限元是通过引入边界节点力以模拟无限域对有限域的约束支承作用，但目前关于该节点力与静力边界是否等效的研究较少；在动力计算中，无限元是一种黏性吸收边界，其无法模拟地基的弹性恢复力及内行波输入，因此，有必要研究基于无限元进行改进的静-动力统一人工边界。

基于此，简要介绍了ABAQUS静、动力无限元人工边界的基本理论，结合二维水平地基模型算例，分析无限元直接应用于初始地应力计算单元节点支承刚度的合理性及计算精度，提出引入静力等效节点力以提高无限元静力及内源计算精度的方法；为使得无限元人工边界能应用于需要考虑模型静、动力综合效应的内、外源振动问题，推导内、外源振动问题中需施加的边界等效节点力，编制相关计算程序，提出基于无限元的静-动力统一人工边界，并对该边界的正确性进行模型算例验证。

1 静、动力无限元人工边界作用原理

根据振动源与研究区域的位置关系，可将岩土动力分析问题大致分为内源振动和外源振动，根据波场相对于模型边界的传播方向，可将模型总波场分解为内行波场和外行波场［19］。内源振动问题的振动源位于研究区域内部，如基础振动引起的土中动应力分析等，波源向外部的无限、半无限区域辐射能量，这个过程属于外行波的散射问题；外源振动的振动源位于研究区域之外，如地震引起地基和上部结构的振动分析等，该类问题可分解为从外部无限域的内行波输入和研究区域内部动力响应产生的外行波散射。由上可知，计算模型进行人工边界截断时，内源振动只需考虑外行波的吸收，而外源振动条件下，外行波吸收和内行波输入均需考虑。

ABAQUS动力分析中，将无限单元与有限单元结合使用，所关心的研究区域（近场有限域）以有限元模拟，而采用无限元模拟远场无限域。当非体力荷载（如集中力等）作用时，静力无限元边界基于Zienkiewicz等［18］的静力计算分析理论，设定距模型计算区域无穷远处位移为零，且假设从有限元-无限元交界面至无穷远处，位移呈线性分布。通过引入满足上述假定条件的映射坐标系以及具有二次收敛速度的位移插值函数，使得无限元成为能够应用于静力计算的边界［20］。当体力荷载（如重力等）作用时，ABAQUS程序因无法在无限元中定义体力，以在有限元-无限元边界自动插入节点力的方式体现无限域对有限域的约束支承作用，但插入的节点力只为使模型处于某一力学平衡状态，而程序并不判断插入的节点力是否代表真实状态［20］，即不检验其是否对应于变形体真实位移，故体力荷载作用下直接用无限元进行初始应力分析有可能获得不正确的计算结果。

ABAQUS动力无限元通过在有限元-无限元边界节点设置一系列法向及切向的阻尼器，模拟半无限空间介质的辐射阻尼［20］，其参考了Lysmer等［1］的黏性边界理论，依据弹性介质中S波、P波的波动方程，推导出反射波为零条件下的边界阻尼系数。由此可见，动力计算中无限元人工边界是一种无静刚度的黏性吸收边界，不能模拟地基的弹性恢复力以及内行波的能量输入。

针对ABAQUS静、动力无限元边界所存在的问题，论文开展了静、动力统一人工边界改进方法研究。静力问题从二维水平地基内初始地应力计算精度出发，探讨无限元对体力荷载作用问题的适用性及改进方法，进一步在此基础上推导内、外源振动问题相应的边界等效节点力公式，对应利用半无限地基表面施加动力荷载和底部输入地震动的两种典型情况进行验证。

2 无限元边界初始地应力生成精度探讨

2.1 分析模型

由上节分析可知，将ABAQUS无限元直接应用于体力荷载作用下的初始应力状态生成，可能会带来不正确的结果。为此，建立二维弹性水平地基模型（图1）进行算例验证，模型侧面及底面设置厚度L2=50 m的无限元人工边界，有限元主模型区域50 m×50 m，网格尺寸0.5 m，有限元、无限元区域分别采用CPE4、CINPE4单元，地基土密度ρ=2 000 kg/m3、弹性模量E=200 MPa、泊松比v=0.25［10］。

图1 计算模型示意图Fig.1 Calculation model（unit：m）

对模型施加重力荷载，利用ABAQUS的Geostatic分析步进行地应力平衡，计算模型尺寸与计算精度呈正相关，采用模型尺寸调整系数ξ探讨边界位置对计算结果的影响，ξ定义为主模型调整后尺寸（L+2L1）与原尺寸（L0）之比即ξ=（L+2L1）/L0，分别按照ξ=1，2，3，4，5，10，15确定有限元计算区域，按0.5 m网格尺寸对模型进行剖分。考察主模型中心沿路径BE的竖向应力，并与理论解答σyy=γ（50-y）作对比（γ为地基土重度，y为计算点纵坐标），记录计算耗时以研究模型尺寸对计算结果及效率的影响，具体方案如表1所示，其中K0为ξ=1并施加静力边界条件（两侧和底部分别采用滚轴、固定边界）的计算结果。

表1 无限元边界条件下模型尺寸对计算结果影响工况表Table 1 Effect size model calculation condition table membered infinite boundary conditions m

2.2 计算结果

经地应力平衡后，K0～K7模型位移量值均小于10-7m，仅从位移大小角度来看，地应力平衡效果均较好。然而考察主模型区域竖向应力分布（图2，图3）、计算精度及耗时（图4）随ξ的变化可知：

图2 不同模型尺寸下主模型区域竖向应力云图Fig.2 Vertical stress contour of the main model with different model sizes

（1）施加传统静力边界条件的K0模型计算结果与理论解完全一致，无限元边界条件下计算精度受模型边界远近影响，若地应力平衡结果要满足工程要求，边界需要达到主模型区域15～20倍以上的远置距离。当ξ=1时，竖向应力明显偏离自重应力理论分布，呈现出沿中轴线BE向模型两侧边界逐渐减小的形式，BE路径竖向应力与理论值的最大相对误差达-33.6%；随着ξ增大，边界逐渐远置，竖向应力分布趋于横向均匀，逐渐接近理论（静力边界条件下）结果；当ξ=15时，模型尺寸已达750 m×750 m，竖向应力与理论值最大相对误差仍有-6.1%，仅近似满足工程精度要求。

（2）竖向应力计算误差与深度基本呈正相关，即ABAQUS静力无限元自动插入节点力的等效刚度与静力边界刚度的差异随深度增加而变大。若需将误差控制在工程允许误差（5%以内），ξ=1时，仅2.5 m深度范围内的计算结果可信（图3）；而按传统经验取3～5倍主模型区域尺寸设置边界进行建模（图4中粉红色区域ξ=3～5），计算结果可信深度范围为8.5～14 m，仍难以满足主模型全区域的精度需要。

图3 各模型沿路径BE竖向应力分布情况Fig.3 Vertical stress distribution along path BE

用图4中所示的三次多项式对相对误差曲线进行拟合，得到尺寸调整系数ξ与误差百分比ω的拟合函数（决定系数0.982）。根据拟合函数进行外延，若误差控制在5%，ξ应不小于20，计算量急剧增加，因此有必要提出一种适用于主模型区域高精度、高效率分析的静力无限元边界改进方法并使之能应用于动力计算的初始状态生成。

图4 竖向应力误差和计算耗时与尺寸调整系数关系Fig.4 The relationship between vertical stress error，calculation time and size adjustment coefficient

3 基于无限元的静-动力统一人工边界改进方法

3.1 静力无限元边界改进

3.1.1 基于静力边界等效节点力的无限元刚度修正

体力荷载作用下，ABAQUS无限元静力计算产生误差的根本原因是程序引入的节点力等效刚度与静力状态的实际边界刚度不匹配，因此，若考虑在施加无限元人工边界的同时，对有限元-无限元边界节点施加传统静力边界（滚轴、固定）计算得到的边界支反力，由于边界支承刚度相等，主模型区域受力状态应与静力边界条件下完全一致。下面对其进行理论验证：

静力计算中，总体平衡方程为

式中：K为刚度矩阵；U为位移向量；F为节点力向量。

将式（1）写成如下形式：

式中：KI、KIB、KBI、KBI为刚度矩阵的分块子矩阵；UI、UB分别表示模型内部与边界节点的位移；FI、FB分别表示模型内部与边界节点的节点力。

由于边界及内部节点力一部分由外荷载引起，一部分由约束反力引起，式（2）可以进一步表示为

式中：上标f表示外荷载引起的节点力分量；而上标r表示约束引起的节点力分量，在施加静力边界时即为静力边界约束引起的节点力，而施加无限元边界时则为无限元边界引起的节点力。

施加静力边界UB=0时，由于边界节点在约束或指定位移方向不会对边界以内节点产生节点力［13］，故F rI=0，于是

可解得

将F rB代回式（3），仍可解得UB=0，且原UI的解也不变，上述推导从物理角度可解释为对模型施加约束条件和对模型施加与约束条件UB=0相对应的等效节点力F rB是等效的，因此可在需要去除静力约束的情况下，施加相应的约束反力以达到不改变原静力平衡状态的目的，下文将F rB称为静力等效节点力。

而当施加无限元人工边界时，程序自动插入的节点力FINFB并不一定满足FINFB=F rB，故将FINFB代回式（4）不一定得到真实位移边界条件UB=0，计算所得的应力状态亦不一定与静力边界一致。若在施加无限元人工边界的同时，对主模型区域与无限元交界面处一并施加静力等效节点力F rB，则主模型区域的应力场将保持与静力边界条件下一致。

3.1.2 模型验证

由上节推导可知，可引入静力等效节点力以提高无限元的初始地应力场计算精度，具体实现步骤为：（1）建立有限元模型，对模型施加静力边界条件，计算完成后，提取模型左、右及底部支座反力即静力等效节点力。（2）去除静力边界，施加无限元人工边界、导入步骤（1）中计算所得的应力场并施加有限元-无限元静力等效节点力即可完成初始地应力生成，静力等效节点力在后续分析过程中保持不变。

为方便对比，取第2.1节中尺寸调整系数ξ=1的模型进行验证。将施加静力等效节点力前后的计算结果、误差连同理论解绘入图5，由图可知，仅设置无限元边界的模型竖向应力误差可达-33.6%，而采用文中引入静力等效节点力方法与理论解完全一致，说明引入静力等效节点力的方法是合理的。

图5 引入静力等效节点力方法与仅设置无限元边界计算所得竖向应力Fig.5 Vertical stress difference between the static equivalent nodal force and the infinite boundary

3.2 静-动力统一人工边界等效节点力分析

3.2.1 内源振动问题

对于需要考虑静、动力综合效应（即静动力合算）的内源振动问题，动力计算只需考虑外行波吸收，无需考虑内行波输入，而无限元本身可较好地模拟半无限空间的辐射阻尼以吸收外行波，因此仅需保证静力计算精度，故内源振动问题的等效节点力等于静力等效节点力：式中：下标α为主模型区域边界位置；β为等效节点力方向。

3.2.2 外源振动问题

对于静动力合算的外源振动问题，则需满足静力计算精度并同时考虑外行波吸收和内行波输入。其中，外行波吸收可由无限元本身解决，对于内行波输入，需克服人工边界的能量吸收才能由域外透射至域内，文献［11］考虑无限域与有限域之间的相互作用力和克服人工边界的阻尼力，导出了外源振动条件下用于地震动输入的等效节点力公式（下文称之为动力等效节点力）；静力计算部分可采用上文第3.2.1节静力等效节点力方法以保证精度，故外源振动条件下，静-动力统一人工边界等效节点力应为模拟内行波输入的动力等效节点力和保证静力计算精度的静力等效节点力共同组成。

将动力等效节点力［11］与静力等效节点力（式（6））叠加即得到外源振动条件下静-动力统一人工边界的等效节点力，下面以（以图1为例）给出平面剪切S波和压缩P波沿y方向从模型底部垂直向上传播时，有限元-无限元边界等效节点力计算公式。

平面P波垂直向上传播时，考虑底部入射波f1和地表反射波f2，则无限元静-动力统一人工边界等效节点力，在主模型左侧AD边界为：

主模型右侧C F边界：

主模型底部DF边界：

同理可得剪切S波垂直向上传播时等效节点力为：

主模型左侧A D边界：

主模型右侧CF边界：

主模型底部BE边界：

式中：下标l、r、b分别表示模型左侧、右侧和底面；x、y分别表示x，y方向；A i为对应节点的控制面积；λ为第一Lamé常数。

3.3 实现流程及适用性验证

3.3.1 实现流程

静、动力等效节点力可分别基于式（6）和式（11）～式（22）进行编程计算。对于内源振动问题，等效节点力即为静力等效节点力，对于外源振动问题，等效节点力则由静力等效节点力和动力等效节点力共同组成。完整的静、动力统一人工边界实现流程如图6所示，文中研究以Python语言编制等效节点力输入程序。

图6 静-动力统一人工边界计算实现流程Fig.6 The implementation method of static-dynamic unified artificial boundary

3.3.2 内源振动问题验证

仍考虑第2.1中的二维水平地基模型，静力荷载仅为重力，在模型地表中心B点施加如图7所示的动荷载，峰值200 kN，作用时间0.5 s，计算时长4 s，对模型施加内源振动问题对应的等效节点力（式6），求解时间步长为0.001 s。为验证计算结果准确性，另建立一远置边界模型作为近似理论解答，根据文献［19］的建议，监测点到人工边界的最近距离应大于求解时间与波速乘积的一半，设置模型尺寸为1 280 m×640 m，采用传统静力边界条件，材料参数、单元尺寸及类型均与主模型区域一致。首先进行地应力平衡，随后进行动力计算。

图7 输入动荷载时程Fig.7 Dynamic load time history

考虑到地应力平衡后模型位移近似为零，位移响应主要取决于动力荷载，应同时考察模型应力与位移响应，监测模型B点的竖向位移及中心G点（深度24.75 m）的竖向应力时程，分别如图8（a）和（b）所示。由图可知，采用静-动力统一人工边界方法计算结果与远置边界基本一致，精度良好；除此之外，采用文中方法耗时305 s，仅为远置边界模型耗时（86 631 s）的3.52%，因此该方法应用于静、动力合算问题可极大提高计算效率。

图8 冲击荷载作用下B点位移及G点竖向应力响应时程Fig.8 Displacement response of point B and the vertical stress response of point G

3.3.3 外源振动问题验证

（1）计算模型及工况

基于第2.1中的二维水平地基模型，静力荷载仅为重力，底部x方向根据地震位移时程曲线（图9，最大幅值0.01 m），利用式（17）确定的等效节点力输入地震波，计算时长4 s，计算对比3种边界形式下模型动力响应差异，其中，工况Ⅰ为传统无限元人工边界；工况Ⅱ为主模型侧面施加文献［11］所推得的动力等效节点力，但不考虑静力等效节点力；工况Ⅲ为主模型侧面施加本文所提出的静-动力统一人工边界。

图9 输入波位移时程曲线Fig.9 The displacement time history of seismic wave

首先进行地应力平衡，随后进行动力计算。监测模型顶部A、B、C点和底部E、F、G点位移时程以及内部G点的竖向应力时程。

（2）计算结果分析

由弹性波动理论可知，当地震波垂直向上传播至地表时，反射波与入射波幅值大小及方向均应一致，地表位移幅值为输入波的2倍；剪切波在传播过程中，不会引起单元正应力而仅引起切应力，因此单元正应力应基本保持不变，即等于前期的静力计算结果。可由上述2项规律，判断各边界条件下位移及应力计算结果的合理性。

顶部A、B、C点及底部D、E、F点的位移响应情况分别如图10（a）～（f）所示，由图可知：

（1）工况Ⅰ由于仅设置了无限元人工边界而未考虑静、动力等效节点力，地表位移（图10（a））响应显著偏小，同时模型底部（图10（d））未出现地表反射波，均与理论响应曲线差别较大，地表A、B、C三点间同步性较差，模型底部D、E、F点间也出现类似情况，说明该边界条件会造成波动场的不均匀；

图10 模型顶部及底部测点位移响应时程Fig.10 Displacement response of the top and bottom points

（2）分别考虑了动力等效节点力和静、动力等效节点力的工况Ⅱ、工况Ⅲ，顶部及底部位移响应均与理论值基本相同，且地表A、B、C点及底部D、E、F点响应时程曲线同步性较好，说明波的传播是均匀的，仅从位移角度来看，2种边界条件下的计算精度均较好。

3种工况下，G点竖向正应力响应分别如图11（a）～（c）所示，由图可知：

图11 模型中部G点单元竖向应力响应时程Fig.11 Vertical normal stress response of element in position G

（1）工况Ⅰ和工况Ⅱ单元竖向正应力响应量值（-326.2 kPa）均与理论值（-495.0 kPa）差异较大，其应力波动基准值与理论值相比均较小，精度较低；

（2）工况Ⅲ单元竖向正应力响应量值均稳定在理论值附近，波动最大幅值仅为理论值的0.011%，可见，采用本文静-动力统一人工边界获得的应力响应是可靠的。

由以上分析可知，虽然工况Ⅱ边界条件下的位移响应与理论值基本一致，但应力响应出现较大失真。岩土动力计算中，第一步往往是地应力平衡，动力加载前模型位移很小（＜10-5m），动位移响应主要取决于动力荷载，然而由于初始地应力的误差影响后续动应力响应结果精度。

3.3.4 边界条件适用性探讨

从监测点位移及应力响应与理论解对比可知，工况Ⅰ位移、应力响应量值及规律均明显偏离理论解，因此传统直接设置无限元的边界条件不适用于外源振动问题；工况Ⅱ位移响应与理论解基本一致但初始地应力出现较大系统误差，因此得到的动应力响应并不可靠，工况Ⅱ边界条件仅适用于静、动力分算的外源振动问题；采用文中的无限元静-动力统一人工边界（工况Ⅲ）的位移及应力响应均与理论解相符，且可根据求解问题的类型考虑按需引入静力、动力等效节点力，使得其对静力问题和静、动力分算（合算）的内、外源振动问题均具有较强适用性。综上，3种工况对应的边界条件对应的适用范围如表2所示。

表2 常用人工边界设置方法及其对岩土静、动力问题适用性Table 2 The common artificial boundary′s applicability to rock and soil static and dynamic problems

4 结论

通过分析ABAQUS静、动力无限元理论，利用二维水平地基模型，分析讨论了无限元直接应用于体力荷载初始应力计算时的精度较差且对模型尺寸要求大的局限性，针对性的提出了引入静力等效节点力的改进方法，建立了无限元静-动力统一人工边界并进行了算例验证，得出以下结论：

（1）因ABAQUS静力无限元引入节点力等效刚度与静力边界的差异，无限元边界直接应用于体力荷载下初始地应力的计算结果误差较大。算例表明，50 m×50 m水平地基的竖向自重应力计算结果严重偏小，最大相对误差可达-33.6%；边界远置方法虽能一定程度改善精度，但所需的模型尺寸过大、计算效率极低。故有必要研究适用于无限元人工边界的高精度、高效率的改进方法。

（2）以无限元应用于静力工况及内、外源振动的普遍适用性为目标，基于模型约束条件和对应节点力的等效性，提出引入静力等效节点力提高初始应力计算精度的方法，并通过模型计算验证了其具有与传统静力边界一致的效果；进一步基于弹性波动理论考虑有限域-无限域相互作用和边界阻尼效应，推导了内、外源振动条件下的等效节点力，提出了无限元静-动力统一人工边界并编制了相应的等效节点力输入程序。

（3）针对半无限地基表面的动力加载和地基深部振动向地表传播的两种典型内、外源振动情况，以模型位移与应力响应为监测量验证了无限元静-动力统一人工边界的有效性，分析了传统无限元与仅考虑动力等效节点力的边界条件对静、动力合算问题的局限性，论文提出的无限元静-动力统一人工边界对内、外源振动的静、动力分（合）算问题均可具有较好的适用性。