竖向成层介质中标量波传播问题的高精度人工边界条件

李会芳，赵 密，杜修力

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

在地下结构抗震分析中，有限元法因其灵活高效的特点，被广泛用于有限域模拟。为了更有效地反映截断无限土体的辐射阻尼效应，需要引入人工边界条件[1]。人工边界的适用性和精度直接影响计算结果的精度，并通过影响计算区域的大小而影响计算效率。城市抗震防灾和重大工程项目建设中经常遇到竖向成层介质以及不规则地形场地。在我国西南部山岭地区中这种地形地质条件较为常见，如滇中引水工程中香炉山隧道穿越的含多条断裂带的复杂地形地质条件[2]。厦门市轨道交通2 号线跨海盾构隧道穿越了基岩、基岩强风化、基岩全风化等多个地层[3]。不均匀介质和不规则地形场地条件引起复杂的波动反射和透射，会造成与均匀半空间场地明显不同的动力反应。复杂场地中结构动力反应研究表明，不均匀地质和不规则地形的存在会对场地和结构的动力反应产生显著影响[4-8]，在地下结构抗震分析领域引起越来越多的关注。建立竖向成层介质中的高精度人工边界条件具有研究和实际意义。国内外研究者已提出大量人工边界条件[9]，其中包括：粘性边界[10]、粘弹性边界[11]、透射边界[12]、无限元法[13-14]、边界元法[15-16]、Engquist-Majda 边界[17]、Higdon边界[18]、Bayliss-Turkel 边界[19]、Dirichlet-to-Neumann(DtN) 边界[20-23]、一致边界(薄层法)[24-26]、比例边界有限元法[27]、吸收层法如完美匹配层法[28-29]等。以上人工边界条件主要用于均匀介质的半空间或全空间场地以及底部固定的单层或水平成层场地。

近年来适用于水平成层半空间场地的人工边界条件得到进一步发展。2009 年，蒋通和田治见宏[30]提出一致边界结合粘性边界来处理水平成层半空间介质的波传播问题，粘性边界中阻尼器的吸能效果对整体精度影响较大。2011 年，Lee 和Tassoulas[31]提出一致边界结合连分式吸收层的方法；2012 年，João 等[32]提出采用一致边界结合完美匹配离散层来解决水平成层半空间介质的波传播问题，分析表明连分式吸收层可等效转化为完美匹配离散层，其中匹配层厚度为基于频率和外行波入射角的纯虚数。2015 年，Hamdan 等[33]提出采用一致边界结合旁轴边界来模拟水平成层半空间介质中波传播问题，旁轴边界采用二阶泰勒展开代替半空间的解析刚度精度较低。2020 年，李会芳等[34]提出H 形高精度人工边界处理水平成层半空间标量波传播问题，此边界应用简单方便且能得到高精度、高效率的求解。

基于水平成层半空间场地中人工边界条件的研究工作，本文针对含竖向成层介质以及不规则地形场地中标量波传播问题，建立了频域下折线形高精度人工边界条件。竖向成层介质的底边界采用广义一致边界来拟合，通过适当的空间变换并精确模拟两侧的半空间刚度，将原来适用于表面自由、底部固定的一致边界扩展应用于两侧开放的竖向成层介质底边界；成层介质侧边界采用基于连分式的高精度边界。由于提出边界的高精度特性，有限域可以取得尽量小，对于无内域结构仅关心地表反应的情况，可以直接将人工边界加在地表，极大地减少自由度，提高计算效率。通过在每个竖层边界上引入斜角变换，广义一致边界可以为任意折线形，适用于竖向成层、不规则地表的场地反应。

1 问题描述

地震作用下重大基础设施结构的动力响应需考虑与周围介质的相互作用，形成无限域介质中的波传播问题。图1(a)所示为二维竖向成层半空间场地中的标量波传播问题，场地中包含任意不规则地形以及广义结构。为高效求解出平面波动问题，引入人工边界分别将竖向地层和两侧无限域截断，如图1(b)所示计算模型1 由人工边界条件与内域有限元部分组成，其中有限域中可以包含任意不均匀、非线性材料以及不规则结构，采用有限元法模拟。由于提出人工边界条件的高精度，有限域尺寸可以尽量小。对于不含地下结构仅关心地面反应的情况如图1(c)所示，可以直接将人工边界加在地表，仅用折线形高精度人工边界条件来模拟竖向成层半空间场地，如图1(d)中计算模型2，最大限度减少计算自由度高效求解地表反应。

图1 竖向成层且地表起伏的半空间场地中标量波传播问题Fig. 1 Scalar wave propagation in half space including vertical stratified media and irregular topography

半空间场地中，左右两侧无限域以及竖向无限延伸的地层内为均匀线弹性材料。频域下标量波传播问题的控制方程为：

地表为自由边界条件，各竖向地层和无限域之间的交界面以及人工边界两侧均满足位移和应力连续条件。另外，满足无限远处辐射条件以及初始静止条件。

2 高精度人工边界条件

以图1(b)中计算模型1 为例给出本文高精度折线形人工边界条件的实现过程。不含广义结构时，人工边界条件可加在地表处，此时计算模型中无内域和竖层侧边界，计算模型1 简化为图1(d)中计算模型2。另外，计算模型1 中竖层底边界可以根据结构形式、计算效率等需要设置为任意角度折线形边界。

2.1 两侧无限域连分式边界

则侧边无限域y方向的动力刚度为：

以上为左侧无限域的人工边界条件，同理可以得到右侧无限域的人工边界条件。

2.2 竖层底部广义一致边界

2.3 有限域与人工边界条件耦合

将有限域侧边力位移关系的连分式展开式(11)沿竖向x轴离散，通过线性插值得到侧边人工边界条件有限元方程：

将有限域与底部人工边界条件式(17)以及两侧人工边界条件式(18)组装，得到频域下耦合系统有限元方程：

式中，下标I、B、C分别对应内域自由度、边界自由度以及连分式引入的辅助自由度。

3 数值算例

本节给出不同地形、地质条件下四种场地中波传播问题的数值算例来验证提出高精度折线形人工边界条件的有效性和精度。参考解为大区域有限元模型的时域解，模型尺寸足够大能够保证在观测时间内结果不受到截断边界处反射波的影响。同时采用了简单常用的粘弹性人工边界条件[11]进行了对比分析，说明提出的折线形人工边界条件适用于竖向成层不规则地形场地，且在计算精度上有显著改进。本文人工边界条件中连分式阶数统一取为5 阶，阶数为零时可退化为旁轴近似。采用Ricker 波作为荷载，其时程和频谱曲线如图2 所示。

图2 Ricker 波的时程和频谱曲线Fig. 2 Time history and Fourier spectrum of Ricker wavelet impulse

3.1 水平地表场地

平直地表半空间场地中含有一个竖向无限延伸地层，如图3 所示。中间地层的剪切波速c1=100 m/s，两侧无限域介质剪切波速c2= 200 m/s。坐标原点取地表中心点A。计算中选取了两种形式的人工边界，如图3 中点划线和虚线所示。加载点为点A，观察点为A(0, 0)、B(0, -100)、C(100, 0)、D(100, -100)。有限元网格尺寸为2 m，计算时间步长为0.01 s。

图3 含竖直地层半空间场地中标量波传播问题Fig. 3 Scalar wave propagation in half space field with one vertical layer

观察点的位移时程结果如图4 所示。将采用新边界的两种计算模型和采用粘弹性边界的计算模型1 与采用大区域有限元模型得到的参考解进行比较。其中计算模型1 由人工边界条件与内域有限元部分组成，如图3 中点划线所示区域；计算模型2 为直接将人工边界加在地表，如图3 中虚线所示。下文中图例意义相同。从图中可以看出，无论是加载点还是边界处节点以及材料交界面上节点，采用新人工边界条件的两种模型计算结果均与参考解吻合较好，而粘弹性边界模型在边界及材料交界面节点上计算误差较大。验证了新提出人工边界条件的有效性和精度。

图4 水平地表场地中观察点位移时程结果Fig. 4 Time histories of displacement solutions at observation points in a flat surface site

3.2 折线地表场地

折线形地表半空间场地中含有两个竖向地层，如图5 所示。左侧地层的剪切波速c1=100 m/s，右侧地层的剪切波速c2= 300 m/s，两侧无限域介质剪切波速c3= 200 m/s。图5 中各点坐标为A(0, 0)、B(50, 50)、C(150, 0)、D(0, -50)、E(50, -50)、F(50,-50)，其中加载点为B，观察点为B、C、E、F。有限元网格尺寸为1 m～3 m，计算时间步长为0.002 s。

图5 含竖向成层介质且地表起伏的半空间场地中标量波传播问题Fig. 5 Scalar wave propagation in half space field with vertical stratified media and irregular topography

观察点的位移时程结果如图6 所示。观察点包括地表节点、人工边界上的节点以及材料交界面处的节点。从图6 中可以看出，两种新边界计算模型的位移解与参考解吻合较好，而粘弹性边界模型在边界及材料交界面节点上计算误差较大。表明提出的人工边界条件可以很好的模拟边界处和材料变化带来的反射和散射波，以及随时间空间变化的辐射效应，验证了新人工边界条件的有效性和精度。

图6 折线地表场地中观察点位移时程结果Fig. 6 Time histories of displacement solutions at observation points in a zigzag terrain site

3.3 阶梯地形场地

阶梯地形半空间场地中含有三个竖向地层，如图7 所示。从左到右介质刚度依次增大，左侧无限域介质剪切波速cl= 100 m/s，三个竖向地层的剪切波速分别为c2= 150 m/s、c3= 200 m/s、c4=250 m/s，右侧无限域介质剪切波速c5= 300 m/s。图7 中定位点坐标为A(0, 0)、B(50, 50/3)、C(100,100/3)、D(150, 50)、E(0, -50)、F(50, -50)、G(100,-50)、H(150,-50)，其中加载点为B，观察点为A、C、E、G。有限元网格尺寸为1 m～2.2 m，计算时间步长为0.002 s。

图7 含竖向成层介质的阶梯地形半空间场地中标量波传播问题Fig. 7 Scalar wave propagation in step-shaped half space field with vertical stratified media

观察点的位移时程结果如图8 所示。观察点包括地表节点、人工边界上的节点以及材料交界面处的节点。从图8 中可以看出，两种新边界计算模型的位移解与参考解吻合较好，而粘弹性边界模型在边界及材料交界面节点上计算误差较大。表明新人工边界条件可以很好的吸收边界以及材料变化产生的反射和散射波，模拟无限场地的辐射效应，验证了新人工边界条件的有效性和精度。

图8 阶梯地形场地中观察点位移时程结果Fig. 8 Time histories of displacement solutions at observation points in a stepped terrain site

3.4 含地下结构场地

如图9 所示，地下结构建立在含竖向成层介质且两侧不等高地形的半空间场地中，其地表高低起伏可模拟河谷与山坡同时存在的复杂地形，本算例分析其中的标量波传播问题。从左到右5 种介质刚度依次增大，左侧无限域介质剪切波速cl= 100 m/s，三个竖向层的剪切波速分别为c2=150 m/s，c3= 200 m/s，c4= 250 m/s，右侧无限域介质剪切波速c5= 300 m/s。因含有地下结构，计算模型2 不再适用或者说误差较大，本节仅采用计算模型1 与参考解比较，人工边界位置如图9中点划线所示，分别施加新边界条件和粘弹性边界条件[11]。图9 中定位点坐标为A(0, 0)、B(50, -20)、C(150, 30)、D(200, 10)、E(0, -50)、F(50, -50)、G(150, -50)、H(200, -50)，地下结构圆心坐标为(100, -30)，半径10 m。其中加载点为B，观察点为B、C、E、H以及地下结构上节点C1(90, -30)、C2(100, -20)。有限元网格最大尺寸小于4 m，计算时间步长为0.002 s。

图9 含地下结构的复杂地形地质半空间场地中标量波传播问题Fig. 9 Scalar wave propagation in complex topographical geological half space field including an underground structure

观察点的位移时程结果如图10 所示。观察点包括地表节点、人工边界和材料界面处的节点以及地下结构上的节点。从图10 中可以看出，新人工边界条件可以很好地模拟边界处、材料界面以及地下结构带来的反射波和散射波，以及随时间空间变化的辐射效应。新边界计算模型的场地反应和结构反应均与参考解吻合较好，而粘弹性人工边界误差较大。验证了新人工边界条件的有效性和精度。

图10 含地下结构场地中观察点位移时程结果Fig. 10 Time histories of displacement solutions at observation points in a site containing underground structure

4 结论

本文针对竖向成层介质中的标量波传播问题，基于连分式展开和扩展的一致边界，建立了含竖向成层介质和地表不规则场地的折线形高精度人工边界条件。半空间场地划分为竖向成层介质以及两侧无限域。两侧无限域的人工边界条件基于动力刚度的连分式展开，引入辅助变量来模拟无限域辐射效应；竖向成层介质为广义的一致边界，两端开放与侧边连分式边界耦合，并引入斜角变换使得提出的人工边界条件可以为任意折线形。侧边界与底边界耦合后的折线形人工边界可用于竖向成层介质中标量波传播分析。数值算例结果表明，提出的高精度人工边界条件适用于多种竖向成层介质中标量波传播问题，具有较高的计算精度，极大提高计算效率。

与水平成层半空间中标量波传播问题的H 形人工边界[34]相比，本文折线形人工边界条件主要是对原有方法的扩展应用。水平成层半空间模型转换到竖向成层半空间模型的过程中，多层和半空间区域的无限延伸方向以及物理边界均发生变化，人工边界的截断位置也随之改变。因此针对水平成层半空间中的人工边界条件不再直接适用，但具有重要参考价值。为进一步扩展人工边界条件的适用范围，适用于倾斜地层的人工边界条件正在研究中。