透射边界结合隐式有限元方法的计算精度探讨

倪同健，丁海平，于彦彦，陈姚杰

（苏州科技大学 江苏省结构工程重点实验室，江苏 苏州215011）

在采用直接法开展土-结构地震反应分析时，需要设定人工边界条件来模拟土体无限域。常用的人工边界很多，如粘性边界[1]、粘-弹性边界[2-3]，叠加边界[4]、旁轴边界[5]、一致边界[6]等，其中多次透射公式（（Multi-Transmitting Formula，MTF））[7]作为一种普适性强，易于实现等优点的局部人工边界条件，特别是将边界节点运动方程和时空解耦的内节点运动方程结合起来，建立了波动数值模拟的解耦技术，大大提高了计算效率和对计算机的内存要求，得到了广泛地关注和应用。随着计算机的高速发展，采用隐式有限元算法与MTF结合进行土-结构相互作用动力分析的计算效率已不成问题。相对地，MTF研究最多的是稳定性问题，特别是针对数值模拟中涉及的高频振荡和低频失稳，不断有新的研究成果产生[8-13]；而MTF的计算精度问题主要是分析与显式有限元方法结合的情形[14-16]，与隐式有限元算法相结合的计算精度的研究并不多，廖振鹏等[7]曾采用频域透射公式做过精度分析，针对SH波散射问题比较了模型计算尺寸L及输入波波长与有限元网格尺寸的关系对精度的影响。本文将针对SV波入射的散射问题，在时域中详细探讨多次透射公式与隐式有限元算法相结合的应用方法，以及时间步长和空间步距对计算精度的影响。

1 系统动力反应分析

如图1所示由人工边界和地表包围的区域结构－地基系统，运动方程为

式中，M为质量矩阵，C为阻尼矩阵，K为刚度矩阵，P为外力矢量。

在动载作用下，本文所讨论结构-地基系统的动力反应是开放系统中的波动问题，其中需要解决好的两个问题：一是边界问题，二是地震波的输入问题，且这两个问题是相关联的。

（1）人工边界。本文的人工边界采用的是多次透射公式（MTF）。设某一入射波以人工波速Ca沿x轴从左侧射向人工边界点0，如图2所示；记点j在p时刻的位移表达式为ujp=u（pΔt，jΔx），j和p为整数，Δt为时间步，Δx为空间步距，若Δt=Δx/Ca，则u0p+1可直接用与人工边界垂直方向上的内部节点的位移ujp+1-j确定[7]

图1 结构－地基系统动力反应

图2 人工边界附近节点编号分析示意图（一维波动模型）

由于MTF模拟的是外行波，即散射波us，其表达式

式中，us为散射波位移，u为全波场位移，ur为参考波场的位移（对底边界，参考波场可直接取输入场；对侧边界，则通常取为自由场）。令多次透射公式（2）中u=us，将式（4）代入式（2）中，则

（2）地震波的输入。假定输入地震波为位移，则输入面为人工边界节点所在位置，人工边界点的全波场位移为

式（6）方程中的u0p+1实际就是结构-地基系统在地震荷载作用下的边界点的位移反应。将式（1）改写成

式中，下标i表示内节点（包括自由表面节点），b表示边界节点。对于散射问题，方程（7）的右边荷载项为0，方程（7）就变成已知边界点ub的位移求解内节点位移问题；对于波源问题，方程（7）的右边Pb为0，方程（7）就变成已知边界点ub的位移和Pb求解内节点位移问题。只要选择一种合适的数值积分格式，两种情形都很容易对方程（7）进行求解。

2 数值模拟

选取分别代表二维均匀水平地形、成层水平地形、凹陷地形的3类场地为算例。假定以SV脉冲波的形式从模型底部入射，其时程曲线和幅值谱如图3所示，输入波的截止频率为20 Hz。

2.1 均匀水平地形

计算区域如图4所示。设无限均匀半空间介质的力学参数：密度为2 000 kg/m3；剪切波速为1 000 m/s；泊松比为0.25。模型尺寸长Lx=100 m，高Ly=50 m，观测点A（0，0）、B（0，50）和C（100，20）的位置见图4。

图3 脉冲波的位移时程和傅里叶谱图

图4 均匀半空间模型

在进行有限元分析时，为了满足数值模拟的精度要求，模型的网格划分一般取每个波长λ具有10个左右单元。为了考察空间步距Δx对计算精度的影响，本文将考虑每个波长采用不同的有限元单元数目，同时采用不同的时间步距，共考虑了8种工况（见表1所列），工况1～6为垂直入射，取一阶MTF；工况7～8为斜入射情形，入射角θ分别取10°和20°，采用一阶和二阶MTF。

表1 均匀半空间计算模型

工况1～6情形下3个观测点位移时程数值模拟结果与精确解的对比见图5（下），图中所有工况的3个观测点的模拟结果与精确解几乎完全重合，很多文献也是通过这样的方法说明数值解的精确度。如果计算二者的误差，即数值解减去精确解，可以发现6种工况的误差还是相差很大的。如图5（上）所示，工况5的误差最小，工况4的误差比工况5的略大，但相差并不明显，二者的Δx相同，Δt分别为0.000 5和0.000 25；工况2比工况1好；工况4比工况3好；工况6较之工况1，同样有所改善但改变不大。

比较误差发现，对于最简单的半空间模型，λ/Δx越大，数值模拟精度越高；Δt越小，精度更高，只是Δt的改变对精度的改善与λ/Δx的改变对精度的影响相比较小；当取λ/Δx=40时，一阶MTF就有相当高的精度。采用工况5的Δt和λ/Δx取值考察斜入射的情况，即工况7～8。根据斜入射的结果，即图6和图7，二阶MTF的精度明显高于一阶MTF，因此需要采用二阶MTF用以提高斜入射的数值模拟精度。

2.2 均匀水平成层地形

设无限均匀水平2层介质的力学参数为：上部土层密度为2 000 kg/m3，剪切波速为500 m/s，泊松比为0.25；下部土层密度为2 000 kg/m3，剪切波速为1 000 m/s，泊松比为0.25。计算区域如图8所示，模型尺寸Lx=100 m，高Ly=50 m；观测点A（0，0）、B（0，50）和C（100，25）的位置见图8。

图5 θ=0°时观测点位移时程和误差比较（均匀水平地形）

图6 θ=10°时位移时程和误差比较（均匀水平地形）

图7 θ=20°时位移时程和误差比较（均匀水平地形）

根据2.1节的均匀水平地形分析，本节考虑的工况9～12（垂直入射）与工况2～5对应，工况13～14（斜入射）与工况7～8对应。因为2层土的剪切波速相差2倍，当采用相同的空间步距划分有限元网格时，2种介质中的每个波长的有限元单元数目也是2倍关系，如工况9中采用Δx=2.5 m，则对于上部土层，每个波长的有限元单元数目等于10，而对于下部土层，每个波长的有限元单元数目等于20。见表2所列。

工况9～12情形下3个观测点位移时程数值模拟结果与精确解的对比见图9（下），图中所有工况的3个观测点的模拟结果的精度看起来都相同，但通过误差对比，如图9（上）所示，可以发现4种工况的误差还是相差很大。

工况11和工况12的误差明显比工况9和10的误差小，其原因是工况11和12中每个波长的有限元单元数目，即λ/Δx，比工况9和10中的多；工况12和工况11之间的误差改变并不明显，同样地，对工况10和工况9而言，其误差虽有所改善也相差不大，可能是由于所选取的Δt对于所探讨的模型而言，其精度已经足够，因此Δt改变后精度虽然有所提高但改善不是很明显。故采用工况12的Δt和λ/Δx取值考察斜入射的情况，即工况13和14。根据斜入射的结果见图10和图11，同样发现采用二阶MTF的精度明显高于一阶MTF。

图8 均匀水平成层模型

表2 均匀水平成层计算模型

图9 θ=0°时观测点位移时程和误差比较（均匀水平成层地形）

图10 θ=10°时位移时程和误差比较（均匀水平成层地形）

图11 θ=20°时位移时程和误差比较（均匀水平成层地形）

2.3 凹陷地形

设在无限均匀半空间表面有一凹坑，土体介质的力学参数为：密度2 000 kg/m3，剪切波速1 000 m/s，泊松比0.25。计算区域如图12所示，模型尺寸Xb=Yb=150 m，Xa=Ya=25 m，高Ly=50 m，模型关于AD轴对称，观察点A（150，125）、B（175，125）和C（175，150）的位置见图12。

根据上述均匀水平地形和均匀水平成层地形的分析，凹陷地形模型考虑15～17三种工况（见表3所列），主要是考察隐式有限元与一阶和二阶MTF结合的精度。

工况15的3个观测点位移时程数值模拟结果与精确解（即远置解）的对比见图13（下），图中二阶MTF的模拟结果与远置解吻合很好；同时通过误差对比图13（上）发现，二阶MTF的精度明显高于一阶MTF。对于斜入射情形，即工况16-17，根据图14和图15，采用二阶MTF的精度也明显高于一阶MTF。

图12 凹陷地形模型

表3 凹陷地形计算模型

图13 θ=0°时观测点位移时程和误差比较（凹陷地形）

图15 θ=20°时位移时程和误差比较（凹陷地形）

3 结论

本文分别选取二维均匀水平地形、成层水平地形和凹陷地形等3种场地模型，采用隐式有限元与MTF相结合的方法，进行了地震波作用下的场地响应模拟，对比了单位网格内包含不同单元数目的数值模拟结果与精确解的差异，同时也考虑了时间步距的影响，得出下列结果。

（1）对于有限元分析中网格大小常常采用每个波长λ有10个单元的情况，一般可以得到较好的计算精度，如果提高λ/Δx的数目，数值模拟的精度将会更高；

（2）时间步距Δt的减小，对于精度同样有所改善，只是可能文中所选取的时间步距已经可以满足所需要的精度要求，因此时间步距的改变对于精度的影响与空间步长的改变对精度的影响相比没有那么明显；

（3）对于复杂地形，若取λ/Δx=20，并采用二阶MTF，可以很大地提高地震作用下场地响应的计算精度。