基于<bold>2.5</bold>维有限元-完全匹配层技术的隧道振动建模研究

赵立财; ZHAO Licai

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

摘要

为解决隧道受周围区域复杂影响的三维隧道动力模拟问题，基于有限单元法，采用2.5D技术，应用完全匹配层（PML）处理由有限单元网格截断的边界，提出一种针对超长隧道结构动态性能模拟的高效、精确解法。结合2.5D技术、有限单元法（FEM）及完全匹配层（PML）的方法并不常见，推导2.5D技术和有限单元法（FEM）结合法的方程，着重说明其与2.5D技术和完全匹配层（PML）结合法的兼容性。通过与常用数值模型（MEF-MEC）和半解析解（PIP）法的结果验证，证明了模型的准确性，并进行了参数研究，以评估隧道和土体性质对土体振动的影响。结果表明：提出的2.5D MEF-PML方法能够解决有限元法不能处理无限域的缺点，与常规PML法相比,精度更高，在隧道动力分析时，能准确考虑隧道周围区域因素的影响。

关键词

2.5D技术; 有限元—完全匹配层法; 波传播; 隧道—土体相互作用

近年来，城市铁路建设引起了一系列振动问题，已经影响到了附近居民的生活质量。由于铁路基础设施与具有住宅或服务功能的建筑之间距离较近，因此，该问题在城市隧道建设中尤为突出。

影响系统响应的因素较多，特别是列车动态激励的产生机制、隧道与周围土体的动态相互作用、隧道的几何形状和位置、岩土土体的几何异质性和固有材料等，因此，隧道振动的产生和传播问题十分复杂。Gupta等^[

1]在开发MEF-MEC类型的混合三维数值模型基础上，运用了动态隧道模拟周期模型，该模型利用有限元法（FEM）模拟复杂的几何体，在隧道中，使用轮廓元法（MEC）模拟与岩土质量相对应的不确定介质。尽管这类模型的精度很高，但其固有的数值计算量很大，不利于实际应用。此外，由于2D模型不能模拟波在三维空间的传播，而在计算模型大、计算步长短的情况下，3D模型时间成本又非常高。假定隧道结构、围岩的几何形状及材料特性沿纵向均匀不变，利用2.5D公式将隧道纵向振动频率等变量进行傅立叶变换，将其分解成一系列沿该方向传播的简谐波，进而得到波数和频域下的应变方程，将3D问题变为平面应变问题。在计算每个简谐波时，只需对隧道二维横截面结构与围岩进行有限元离散化，得出三维解，大大缩减了模型的自由度，使计算效率较3D模型有所提高。2.5D数值模拟技术已被应用于土体和隧道等铁路运输基础设施的研究中^{[参考文献 2-3}2-3]，该方法既可以应用于基于有限元法和轮廓元法模型^{[参考文献 4

百度学术}4]，也可以应用于基于这两种方法的混合模型^{[参考文献 5-6}5-6]，甚至是无网格法与FEM的结合，如Amado-Mendes等^{[参考文献 7

百度学术}7]进行的初步研究。在静态问题中，地基的作用体现在刚度上，可以在距离观察区域足够远的地方截断域，使变形足够小，从而忽略不计。然而，在动态分析中，地基所采用的模型不仅要满足典型土体的动态刚度要求，还要满足索末菲尔德^{[参考文献 8-9}8-9]的辐射条件，该条件假定在各向同性和(半)不确定性的均质环境中辐射能量不可能返回到观察领域。这一要求需要对边界条件进行特殊处理，因为当波到达有限元网格的极限时，不能出现杂散反射。为实现这一目标，针对隧道之间的相互作用问题，已经提出并实施了一些策略。

在使用MEF时，由于观察领域的几何限制，处理人工边界的不同方法和手段可以分为局部程序和整体程序。总体程序在概念上比较一致，但其复杂程度和通用性较低。这类方法包括在有限元背景下精确考虑波传播方程和索末菲尔德辐射条件的方法，如边界元法、基本解法^[

10-11]。这些程序由基于波传播理论的近似解组成，通过处理人工边界来避免波反射。尽管这些方法的效率取决于离散域的大小（相对于产生的波长），但易用性加上与MEF方法的良好兼容性，使其可以广泛运用于实际工程中。局部程序包括吸收边界法、无限单元法和吸收层方法，吸收边界法是其中最简单的方法。在局部程序中，还有另一种处理波传播问题的方法：完全匹配层（PML）方法，也称为PML公式，通过有限大小的层（称为PML层）替换无限域来处理人工边界。该介质能实现两个主要目的：确保当波到达其边界时没有反射；吸收影响它的波。

介绍2.5D技术在有限元法中的运用，强调通过大型隧道进行系统模拟时需要考虑的具体内容；提出2.5D技术与PML的结合法，并强调了该方法与有限元法的异同点。在提出两种方法之后，通过与不确定介质的分析解决方案进行比较来验证这两种模型，将插入不确定介质中隧道的解与PIP（管套管热交换器）半解析程序提供的解进行了比较。最后，基于观察区域横截面上瞬时波的传播参考变化函数,分别进行了隧深度、衬砌厚度及土体材料阻尼特性对振动传播的影响研究。

1 2.5维MEF-PML法数值模型的概念和计算公式

1.1 概况

假设一个纵向截面不变的结构，如图1所示。通过运用2.5D公式在纵向坐标（x方向）进行傅立叶变换，将所有变量变换为波数k₁，然后对截面进行离散化即可得出问题的三维解，最后，在2.5D PML单元层辅助的同时，通过使用2.5D有限元离散观察领域。其中，域截断发生在PML单元之后。因此，该数值装置可以吸收以任意入射角到达2.5D MEF和2.5D PML方法观察域之间边界的波，而且不会出现杂散反射。

图1 2.5D有限元模型

Fig. 1 2.5D finite element method

通过改变复数域PML单元的节点y和z坐标，使x坐标变换到频率—波数域，进而满足边界吸收条件，实现波沿PML单元传播人为衰减的目的，如图2所示。

图2 PML层内波衰减示意图

Fig. 2 Schematic diagram of wave attenuation in PML layer

该方法的关键在于所采用参考函数的变化，因为PML层不仅要发挥吸收材料的作用，还要发挥保证波在观察领域和PML边界处不发生反射的作用。只有满足这两个要求，才能在观察领域内获得准确解。

1.2 2.5D MEF-PML结合法公式

1.2.1 MEF 2.5D结合法的通式

为了减少三维模拟极高的计算成本并减少计算量，对于几何和力学特性无限延伸的三维结构，可以通过2.5D计算公式来处理。由于结构的主要属性由其横截面表示，故被视为二维结构。假定结构响应为线性，则可以在频率/波数域中进行分析，并在无限发展结构方向上对时间和空间进行傅立叶变换。转换后的量（作用量和响应量）成为空间方向 $x$ 和时间方向 $t$ 的傅里叶图像函数，傅里叶图像分别由波数 $k_{1}$ 和角频率 $ω$ 表示，在未处理域中只保留两个方向来表征三维环境，并通过有限元对这些方向定义的平面进行离散化，如图1所示。

基于有限元法的一般步骤，即强形式和弱形式，可以为三维域中任何无限体积建立平衡方程

\int_{V} δ ε σ d V + \int_{V} δ u ρ \ddot{u} d V = \int_{V} δ u p d S

（1）

式中： $δ ε$ 为虚拟扩展场； $σ$ 为应力场； $δ u$ 为虚拟位移场； $u$ 为位移场； $ρ$ 为介质密度； $p$ 为施加在表面 $S$ 上的外部压力场。

应用傅立叶变换后，域的横截面保留在非变换域中并被离散化为有限元，使得根据节点变量重写方程式（1）具有可能性。虚功定理假设为

\int δ f (x) p (x) d x = \int δ f (- k_{1}) p (k_{1}) d k_{1}

（2）

由式（1）、式（2）可知，变换域内力和惯量所做的虚功分别为

\begin{array}{l} \int_{V} δ ε σ d V = \int_{k_{1}} δ u_{n}^{T} (- k_{1}, ω) \int_{x} \int_{y} B^{T} (- k_{1}) \cdot \\ D B (k_{1}) d y d z u_{n} (k_{1}, ω) d k_{1} \end{array}

（3）

\begin{array}{l} \int_{V} δ u ρ \ddot{u} d V = - ω^{2} \int_{k_{1}} δ u_{n}^{T} (- k_{1}, ω) \cdot \\ \int_{Z} \int_{y} N^{T} ρ N d y d z u_{n} (k_{1}, ω) d k_{1} \end{array}

（4）

式中： $B$ 为形式函数导数的矩阵； $N$ 为形式函数矩阵（仅在横截面平面上定义）； $D$ 为弹性矩阵（应力与变形相关）； $n$ 为节点位移矢量（集合变换域内单元截面节点3个自由度对应的位移）；其中激发源被认为在角频率 $ω$ 下随时间呈谐波变化。

利用几何体在ZY平面上仅为离散的事实，计算外部荷载所做的虚功。因此，考虑平行于施加动作元素的坐标s，外部荷载所做的虚功由式（5）得出。

\int_{S} δ u p d S = \int_{k_{1}} δ u_{n}^{T} (- k_{1}, ω) p_{n} (k_{1}, ω) d k_{1}

（5）

式中：向量 $p_{n}$ 汇集了沿有限元边界施加的压力 $p$ 所产生的等效节点力。

由于式（3）、式（4）、式（5）适用于任意虚位移 $δ u_{n} (- k_{1})$ ，因此，可以取 $k_{1}$ 的所有积分，使每个有限元的平衡关系为

(\int_{x} \int_{y} B^{T} (- k_{1}) D B (k_{1}) d y d z - ω^{2} \int_{z} \int_{y} N^{T} ρ N d y d z) u_{n} (k_{1}, ω) = p_{n} (k_{1}, ω)

（6）

通过采用有限元法的常用命名法，厚度矩阵和质量矩阵可分别定义为式（7）和式（8）。

K = \int_{z} \int_{y} B^{T} (- k_{1}) D B (k_{1}) d y d z

（7）

M = \int_{z} \int_{y} N^{T} ρ N d y d z

（8）

矩阵B是由微分算子L（在变换域内）与矩阵N的乘积决定，矩阵N以元素的形式将函数分组。由于方向 $x$ 对波数 $k_{1}$ 进行变换，因此，对同一方向的导数进行解析计算时，矩阵L为

L = {[\begin{matrix} i k_{1} & 0 & 0 & \frac{\partial}{\partial y} & 0 & \frac{\partial}{\partial z} \\ 0 & \frac{\partial}{\partial y} & 0 & i k_{1} & \frac{\partial}{\partial z} & 0 \\ 0 & 0 & \frac{\partial}{\partial z} & 0 & \frac{\partial}{\partial y} & i k_{1} \end{matrix}]}^{T}

（9）

由于在模型中采用了材料阻尼效应的形式，因此，在考虑复杂刚度参数基础上使用了滞后阻尼模型。通过将刚度矩阵K分解为子矩阵，使相同矩阵中的所有项都与波数无关，提高了计算效率。为此，矩阵B被认为是两个矩阵的和，从而将数值导数与解析导数分离，如式（10）所示。

B = B_{1} + i k_{1} B_{2}

（10）

将式（10）代入式（7）中并对其进行分组，就能将式（6）重写为

(K_{1} + i k_{1} K_{2} + k_{1}^{2} K_{3} - ω^{2} M) u_{n} (k_{1}, ω) = p_{n} (k_{1}, ω)

（11）

1.2.2 2.5D与PML结合法公式

为了克服MEF本身固有的局限性，将采用2.5D和PML相结合的吸收边界法。其原理是在观察区域周围引入一层PML单元，该区域也将由2.5D有限元离散，这意味着它能吸收入射能量，而且不会在两个区域之间的边界发生反射。因此，PML层应该具有吸收性，而不是反射性，如图2所示。

为了满足上述条件，需考虑PML域中的复数坐标，维持问题的控制方程，从而避免在观察域中的边界上出现虚假反射。鉴于上述情况，需将PML所描述的域更改为复域。由于方向x被变换到波数域（参见图1），因此，该关系由式（12）和式（13）仅适用于y和z方向的参考变化函数给出。

\tilde{y} = \int_{0}^{y} λ_{y} (y) d y

（12）

\tilde{z} = \int_{0}^{y} λ z (z) d z

（13）

式中： $λ_{y}$ 和 $λ_{z}$ 分别为 $y$ 和 $z$ 方向的变换函数。

在制定有限元类型的公式时，还需在物理和复数这两个空间的节点坐标导数之间建立关系。这些关系可以通过式（14）、式（15）建立。

\frac{\partial}{\partial \tilde{y}} = \frac{1}{λ_{y} (y)} \frac{\partial}{\partial y}

（14）

\frac{\partial}{\partial \tilde{z}} = \frac{1}{λ_{z} (z)} \frac{\partial}{\partial z}

（15）

由于PML域内的解满足观察域内微分方程的要求，因此，如果式（1）符合复数空间的坐标变化关系，则可以使用定义2.5D和MEF结合法的程序来定义2.5D和PML结合法。因此，在运用式（12）~式（15）的基础上，结合Garlekin程序，可以为PML区域建立刚度矩阵K^*和质量矩阵M^*

K^{*} = \int_{z} \int_{y} B * (- k_{1}) D B^{*} (k_{1}) λ_{y} λ_{z} d y d z

（16）

M^{*} = \int_{z} \int_{y} N^{T} ρ N λ_{y} λ_{z} d y d z

（17）

矩阵B^*为微分算子L^*与矩阵N的乘积，矩阵N为元素2.5D PML形式函数的组合。为了遵循物理空间坐标和复数空间坐标之间的关系，新的微分算子由式（18）给出。

L^{*} = {[\begin{matrix} i k_{1} & 0 & 0 & \frac{1}{λ_{y} (y)} \frac{\partial}{\partial y} & 0 & \frac{1}{λ_{z} (Z)} \frac{\partial}{\partial z} \\ 0 & \frac{1}{λ_{y} (y)} \frac{\partial}{\partial y} & 0 & i k_{_{1}} & \frac{1}{λ_{z} (Z)} \frac{\partial}{\partial z} & 0 \\ 0 & 0 & \frac{1}{λ_{z} (Z)} \frac{\partial}{\partial z} & 0 & \frac{1}{λ_{y} (y)} \frac{\partial}{\partial y} & i k_{_{1}} \end{matrix}]}^{T}

（18）

PML区域可以在两个方向、一个方向或没有方向上吸收入射能量。后一种情况对应于 $λ_{y} = λ_{z} = 1$ 。将该关系代入式（16）和式（17）表示的刚度和质量矩阵中，可以分别得出式（7）和式（8）。

最后，在代入有限元和PML元的动态刚度矩阵及诺伊曼和狄利克雷边界条件之后，方程组的整体形式即可被确定，如式（19）所示。

\begin{array}{l} \{[K_{F E M}^{g} (k_{1})] + [K_{F E M}^{g} (k_{1}, ω)] - \\ ω^{2} ([M_{F E M}^{g}] + [M_{F E M}^{g} (k_{1}, ω)])\} u_{n} (k_{1}, ω) = p_{n} (k_{1}, ω) \end{array}

（19）

方程组的解位于变换域中，需进行两次傅立叶逆变换才能将解转换到时空域中。

1.2.3 参考系的转换

在观察区域之外，通过解决方案的连续性推导出这些函数，与此同时，提出其衰减演化时的平稳解决方案，以避免在PML域内出现杂散反射。通过采用式（20）、式（21）所示函数来满足要求。

λ_{y} (y) = \frac{2 π}{a b s (k)} \frac{y}{H_{y}} - i \frac{k_{0}}{k} {(\frac{y}{H_{y}})}^{2}

（20）

λ_{z} (z) = \frac{2 π}{a b s (k)} \frac{z}{H_{z}} - i \frac{k_{0}}{k} {(\frac{z}{H_{z}})}^{2}

（21）

式中： $k_{0}$ 是一个合适的常数值（本研究的结论是 $k_{0} = 20$ 时将产生较好的结果）； $H_{y}$ 和 $H_{z}$ 分别为y和z方向上吸收层的厚度；k为问题横截面上传播波长对应的波数，由式（22）给出。

k = \sqrt[]{{(\frac{ω}{C_{s}})}^{2} - k_{1}^{2}}

（22）

式中： $C_{s}$ 为切割波在场内的传播速度。

为了调整网格，以获得足够传播波长的PML域厚度，引入参考变化函数的实际部分。因此，在所提出的方法中，PML网格总是以1 m的厚度生成，并被分为5或6个子层，该数字对应于与结果敏锐度相适应的最小离散值，其适应性由参考变化函数决定。实际上，PML域的网格是通过波数自适应的（拉伸或收缩），这种效果通过变换函数实现。通过将所提出模型得到的结果与理论方案的结果进行比较，验证了所提出的参考变化函数的正确性。

应当注意的是，如果式（22）的根为负值，则参考变化函数没有虚部，即，沿着问题横截面传播的波是瞬变的，因此，参考变化函数式（20）和式（21）是纯实数值，波没有衰减，为了更好地模拟域的无限特性，需要对网格进行调整。

2 理论验证

2.1 不确定介质中的垂直载荷

为了验证模型的可靠性和准确性，对一个简单例子进行模拟，并将结果与Kausel等^[

12]提出的显式分析解决方案进行比较。图3显示了问题的主要特征及采用的2.5D有限元网格的几何形状。在数值分析中，利用问题的对称性和反对称性条件，确定2.5D有限元网格的边界条件，如图3（b）所示。

（a）不确定介质性质

（b） 2.5D有限元网格模型

图3 验证问题的描述

Fig. 3 Verify the description of the problem

图3中观察区域为3 m×3 m的截面，由2.5D有限元单元格表示，每个单元的宽度 $l_{e}$ 为0.3 m。沿着该区域的外部轮廓，即人工边界处施加一层厚1 m的PML单元，其尺寸可根据在YOZ平面上传播的波长进行调整。不确定介质由垂直电荷激发，该垂直电荷相对于x坐标的纵轴为谐波，波数为 $k_{1}$ 。在空间频域中，荷载由式（23）给出。

p_{j} (x, y, z, ω) = δ_{i j} δ (y) δ (z) c o s (k_{1} x)

（23）

由于电荷是谐波，由此产生的位移也是谐波

u_{j} (x, y, z, ω) = u_{j} (k_{1}, y, z, ω) c o s (k_{1} x)

（24）

图4和图5分别显示了 $k_{1}$ 波数（定义为 $k_{1} = k_{1} C_{s} / ω$ ）在频率为10、75 Hz时沿y轴的垂直位移。

（a）波数维度k₁=0.5

（b）波数维度 k₁=1.0

（c）波数维度k₁=1.5

图4 垂直载荷频率为10 Hz时沿不确定介质y轴垂直位移 $u_{z}$ $(x, y, z, ω)$ 的实部和虚部

Fig. 4 When the vertical load frequency is 10 Hz, the real and imaginary parts of the vertical displacement along the y-axis of the uncertain medium are presented

（a）波数维度k₁=0.5

（b）波数维度k₁=1.0

（c）波数维度k₁=1.5

图5 垂直载荷频率为75 Hz时沿不确定介质y轴垂直位移 $u_{z} (x, y, z, ω)$ 的实部和虚部

Fig. 5 When the vertical load frequency is 75 Hz, the real and imaginary parts of the vertical displacement along the y-axis of the uncertain medium are presented

由图4和图5可以看出，无论波数或频率如何，数值解与理论解的一致性都很好。需要注意的是，当k₁大于1时，沿横截面仅产生大于1的瞬时波，如式（22）所示。因为所考虑的材料阻尼非常小，所以，当k₁大于1时，响应的虚分量几乎为零，说明了响应的易逝性，数值模型也很好地再现了这种效果，数值结果与理论结果完全吻合。从图4（a）可以看出，该模型能够准确地再现传播波，即使波长大大超过了观察区域。在图4（a）所示的条件下，传播的切割波波长约为14.5 m，相当于观察区域尺寸的5倍左右。该模型具有显著的再现能力是因为PML所采用的变化参考函数配置式（20）、式（21）在截面上引入了实数分量作为传播波数的函数，使模型能够自动调整PML单元的尺寸，从而模拟大波长。事实上，如果k₁=0，2.5 D问题就转变成了二维平面变形分析问题。对于这种情况，假设无限空间具有图3（a）所示的特性及图3（b）所示的模型。考虑采用Ricker子波脉冲相对应的瞬态垂直载荷如图6所示，动作的应用点对应于参考点的原点，如图3（a）所示，Ricker子波脉冲引起的垂直移动如图7所示。

（a）时域表示

（b）频域表示（频率含量）

图6 Ricker子波脉冲对应的瞬态垂直载荷

Fig.6 Transient vertical load corresponding to the Ricker subwale pulse

（a） A点

(b) B点

图7 Ricker子波脉冲引起的垂直移动

Fig. 7 Vertical movement caused by Ricker wavelet pulse

图7显示了图3（b）所示位于观察区域边界的点A、B在时空域内的响应，在数值技术结果的基础上，给出了理论计算结果（使用Kausel等^[

12]提出的Green 2.5D函数计算）。同样，数值解和理论解之间具有极高的一致性，即使在观察区域边界上的点，所获得的结果之间也没有区别，可以体现出所提出模型的高效性。围绕相同的验证示例，图8给出了更全面的计算结果，描述了一些时间段的垂直位移场。

（a） t=0.214 s

（b） t=0.219 s

（c） t=0.224 s

（d） t=0.229 s

图8 Ricker子波脉冲在不同时间点引起的垂直位移场（单位：cm）

Fig. 8 Vertical displacement field caused by Ricker wavelet pulse at different time points (Unit: cm)

在图8中仅第二象限用数值方法计算，其余象限的结果为Kausel等^[

12]提出的理论解。可以看出，Ricker子波脉冲在不同时间点引起的垂直位移场大小分布形式基本一致，对于各时刻，数值解和理论解基本一致，说明提出的理论可信，可用于进一步研究。

2.2 不确定介质中的隧道

选择Gupta等^[

1]研究中给出的一组结果作为对比组,所考虑的案例用于验证该研究中的数值模型，为此，将数值解与Forrest等^{[参考文献 13

百度学术}13]开发的PIP（管套管热交换器）模型提供的半解析解结果进行比较。比较了3种分析不确定介质中隧道动态响应的方法：本文提出的2.5D MEF-PML模型、Gupta等^{[参考文献 1

百度学术}1]提出的周期性MEF-MEC模型；半解析解的PIP模型。

图9所示为不确定介质性质的隧道,隧道内半径为2.75 m，混凝土壁厚为0.25 m，在其底部施加垂直谐波点载荷。观测点由问题截面上的极坐标和正交方向上的矩形坐标定义，x=0 m，载荷施加在坐标为（0，2.75，0）的点上。有限元网格宽10 m，高40 m，即相对于隧道中心在z方向延伸20 m（A点和C点位于有限元域和PML域之间的界面）。将该域分解为2.5D有限元，每个单元的宽度小于0.6 m，以获得频率高达80 Hz的精确解。为了避免波反射到达人工边界，沿人工边界铺设了一层厚1 m、由5层2.5D有限元组成的PML单元。图10为由在隧道底部施加不同频率谐波点载荷引起的沿YOZ平面（即由载荷和观测点组成的平面）位移场标准的实际部分。

（a）问题描述

（b） 2.5D有限元网格

图9 不确定介质中的隧道

Fig. 9 Tunnel in uncertain medium

（a） 20 Hz

（b） 40 Hz

（c） 60 Hz

（d） 80 Hz

图10 由在隧道底部施加不同频率谐波点荷载引起的沿YOZ平面的实际位移场

Fig. 10 Actual displacement field along the YOZ plane caused by applying harmonic point loads with different frequencies at the bottom of the tunnel

尽管图10中显示为二维形式，但其模拟结果对应于三维分析，解决了波数维度k₁连续值的问题。不同模型结果的比较如图11所示，图11显示了图9（a）中所示的每个观测点垂直位移传递函数的变化规律。

（a） A点

（b） B点

（c） C点

图11 不同位置垂直位移的传递函数

Fig. 11 Transfer function of vertical displacement at different positions

由图11中可以看出，提出的模型具有很高的准确性，它提供的数值结果与使用PIP模型所获得的结果几乎完全一致。PIP模型是一个基于厚壁管道中波传播理论的半解析解模型，可以在当前条件下（隧道在不确定介质中）获得精确解。此外，提出的数值模型具有很强的通用性，可以模拟复杂的几何形状，这在半解析模型（如PIP模型）中不可能实现。

从上述结果可以看出，2.5D MEF-PML模型在精度方面可以取代Gupta等^[

1]提出的MEF-MEC模型，如图11（a）所示的函数波动在2.5D MEF-PML模型中得到了很好的呈现，但在MEF-MEC模型中却不那么明显。

同样，提出的模型即使在观察域（即MEF 2.5D所描述的域）和PML域之间边界上的点也具有极高的准确性。

3 隧道振动的参数

3.1 参考方案说明

参考方案所采用的几何形状和力学性能如图12所示。模型假定为质量均匀的浅埋隧道，通过在隧道底部施加单位振幅（1 N）且频率可变的谐波垂直载荷产生动态激励。图12展示了土体上不同位置的观测点及荷载作用点，观测点和作用点应记录在同一横截面上（静载条件）。在数值模拟时，利用所分析问题的对称性条件，建立了2.5D MEF-PML模型。观察区域为20 m宽，也就是说，图12中所示F点的垂直方向对应于该区域的边界。有限元法可以分析高达80 Hz的频率范围。与观察区域相结合，使用了厚1 m的2.5D PML层，并将其分为6个子层。图13为参考方案所采用的有限元网格。

图12 几何形状和力学性能

Fig. 12 Geometry and aesthetic performance used in the reference scheme

图13 参考方案中的有限元网格

Fig. 13 References in the finite element mesh

3.2 隧道深度的影响

从土体动态响应的角度分析，隧道深度是一个重要方面。除参考方案中H=9 m外，还假定另外两个深度，分别为H=12 m和H=15 m。图14显示了位于土体表面的参考点和3种不同情况下的垂直位移。

（a） A点

（b） B点

（c） C点

（d） D点

（e） E点（f ） F点

图14 不同隧道深度下观测点的垂直位移

Fig. 14 Vertical displacement of observation points at different tunnel depths

由图14可以看出，土体表面的动态响应趋势相当复杂。造成这种复杂性的因素有很多，包括压缩波、截止波和瑞利波在自由表面的相互作用，这种相互作用的结果是曲线上出现波动和反向峰（下降），反向峰之间的频率间距随隧道深度的不同而变化，原因是它们的发生取决于源与接收器之间的距离。

研究显示，隧道越深，振动水平越低，在距离隧道路线最近的点（图14（a）、（b））和远离隧道路线的点（图14（e）、（f））体现得非常明显，当频率大于20 Hz时，后者的影响更加明显。这种依赖于隧道深度的衰减效应主要归结于两个原因：1）阻尼。因为从源到接收器的距离随着隧道深度的增加而增加；2）激励产生射线波的频率范围减小。第1个原因的影响主要体现在离隧道最近的点上，接收源距离受隧道深度的影响最大；第2个影响因素主要体现在离隧道最远的点上。对此，Gupta等^[

1]研究表明，瑞利波的频率含量被限制在

C_{R} / d

频率范围内，其中

d

为源深，

C_{R}

为瑞利波的速度，因此，隧道深度增加导致频率范围缩小的原因主要是瑞利波。当H=9 m时，瑞利波的截止频率约为20~30 Hz，而当隧道深度较深时，截止频率约为13~20 Hz。

另一方面，图14中显示了隧道深度与产生最小响应值频率的关系。隧道的深度越深，该影响发生的频率就越低，低频与高频范围内的频率差异就越小，并且，随着观测点距离的增加，这种差异也会明显减小，如图14（f ）所示。

对于施加荷载作用下点的动态响应，隧道深度的影响可以忽略不计。如图15所示，不同隧道深度荷载作用点所获得的位移幅值只有轻微的差异，尤其是在更深的隧道中。

图15 不同隧道深度荷载作用点的垂直位移幅值

Fig. 15 Vertical displacement amplitude of the different tunnel depth load points

综上可知，从铁路交通振动的角度来看，该结论非常重要。事实上，列车运行过程中产生的动态荷载取决于轨道底部的动态性能，由于不同深度隧道的动态响应十分相似，特别是当隧道上方的土体厚度超过其直径的1.5倍时，列车与系统其他部分(轨固隧道)之间的动态相互作用问题与隧道的深度无关。因此，铁路交通在某一深度隧道中产生的动力相互作用载荷可以作为激励源（输入数据），用于研究不同深度隧道中产生的振动。

3.3 隧道衬砌厚度的影响

隧道衬砌厚度是影响隧道在平面内的动力性能和纵向抗弯刚度^[

14]从而影响隧道与土体相互作用机制的参数之一。为了研究隧道衬砌厚度对系统整体动态性能的影响，设计了两种新方案，即衬砌厚度d=0.5 m和d=0.7 m，隧道外半径为常数，其值如图12所示。与隧道直径相比，衬砌厚度的假设值似乎有些夸张，但从理论角度来看，重要的是要考虑衬砌厚度的高对比度，以突出该参数的影响。

在考虑的3种情况下，根据激励频率，观测点在地块表面的垂直位移情况如图16所示。结果可分为两大类：1）8 m以下距离的动态响应，如图16（a）、（b）、（c）所示；2）8 m以上距离的动态响应，如图16（d）、（e）、（f ）所示。从图16（b）、(c）可以看出，随着隧道刚度的增加，第一个反向峰向更高的频率移动，这可能是因为在平面上自由隧道弯曲模态固有频率的增加由其衬砌厚度的增加引起。8 m以上距离的动态响应总体趋势是，随着衬砌厚度的增加，垂直位移减少，尤其是在频率超过35 Hz的情况下，这种影响随着源与接收器距离的增加而变得更加明显，这在很大程度上是由于隧道纵向弯曲刚度的增加使所施加的荷载进一步退化，从而使响应最小化。显然，这种影响在较短的波长，即较高的频率下更为明显。

（a） A点

（b） B点

（c） C点

（d） D点

（e） E点（f ） F点

图16 不同隧道衬砌厚度下观测点的垂直位移

Fig. 16 Vertical displacement of observation points at different tunnel lining thicknesses

另一方面，如图16（d）、（e）、（f ）所示，随着衬砌厚度的变化，反向峰频率增加，这说明了隧道刚度对隧道与土体相互作用机制的影响。衬砌厚度的增加导致隧道在平面和纵向的固有频率增加，从而导致这种影响发生的频率增加。

综上所述，当源与接收器之间的距离超过隧道直径的1.5倍时，增加隧道衬砌的厚度是一种减少交通诱发振动的措施，但对于较短的距离则效果不大。关于衬砌厚度对隧道自身动态响应的影响，图17比较了3种情况下荷载作用点的垂直位移。

图17 不同隧道衬砌厚度下荷载作用点的垂直位移

Fig. 17 Vertical displacement of load points under different tunnel lining thicknesses

由图17可以看出，在3种情况下，荷载作用点垂直位移值的总体趋势相似，然而，衬砌厚度的增加将导致位移幅度的显著减小，因此，隧道衬砌厚度可以在列车—轨道—隧道的相互作用机制中发挥重要作用。衬砌厚度的增加会导致轨道支撑刚度的增加，这会对轨道交通引起的振动形成机制产生影响，特别是该机制是由车辆和其余系统之间的动态相互作用产生的情况，而这种影响的大小取决于轨道的弹性特性。

3.4 土体物质阻尼的影响

隧道内部引起的振动是一个复杂问题，不仅涉及隧道的力学性能，而且涉及岩体的刚度、分层和阻尼等特性。土体阻尼参数会随着震源与接收器之间距离的变化而变化，很难实现现场量化，因此，应该受到特别关注^[

15]。为了说明土体材料阻尼的影响，提出一项参数化研究，除参考方案阻尼系数

ξ

= 0.04外，还分析了阻尼系数为0.02和0.06的两种情况。振动幅度的衰减由两种机制引起：几何阻尼和材料阻尼，就几何阻尼而言，并不代表系统的能量损失，而是波的传播过程中导致的系统扩散，因此，只取决于波前的几何形状、波的类型及从观测点到发电源的距离。笔者只讨论材料阻尼的影响，是因为几何阻尼只是系统自身结构的一个固有特征，而材料阻尼与有效的能量耗散有关。激励频率越高（波长越短），源与接收器之间的距离越远，其影响就越明显，这种效应在图18中体现得十分明显，图18所示为以质量表面各点垂直位移振幅作为激励频率的函数。

（a） A点

（b） B点

（c） C点

（d） D点

（e） E点（f ） F点

图18 不同材料阻尼值下参考的垂直位移

Fig. 18 Vertical displacement of the reference at different material damping values

由图18可看出：1）阻尼的增加与振动水平的降低有关；2）无论源—接收器距离如何，激励频率越高，阻尼效应越明显；3）当源—接收器距离较大时，材料阻尼引起的振幅明显减小（见图18（f））。这些表现是合理的，因为阻尼效应与材料所经历的振动周期数量有关^[

16-17]，对于较高的频率，即较短的波长，当源—接收器的距离较远时，即材料所承受的循环次数较多时，振动水平将大大降低。关于土体阻尼对隧道阈值点动态响应的影响，结果如图19所示。

图19 不同材料阻尼值下荷载作用点的垂直位移值

Fig. 19 The vertical displacement value of the load action point under different material damping values

由图19可知，在评价隧道本身的位移时，材料阻尼的影响几乎可以忽略不计。通过结论可以发现，材料阻尼对动态列车与线路相互作用机制的影响很小，因此，与传播过程中的情况相反，它对振动产生过程的影响并不明显。

4 结论

1）介绍并验证了2.5D MEF-PML模型，该计算方法适用于分析隧道等交通基础设施无限发展结构的动态响应。

2）由于PML依据的方程与有限元分析中常用的方程非常相似，因此，将2.5D PML结合法代入基于2.5D MEF的数值计算中相对简单，这种合并的基本方法是通过定义和建立参考变化函数来实现的，能够处理域横截面上瞬时波和波的传播参考变化函数。

3）通过比较发现，2.5D MEF-PML模型解与PIP模型结果几乎完全重合，并且，较常用的MEF-MEC模型，其更加简洁，且结果更加准确。事实上，即使使用与观察区域大小相同的有限元网格，2.5D MEF-PML模型也可以获得非常精确的结果，并且该方法克服了有限元程序不能处理无限域的缺点，2.5D MEF-PML方法通过使用一个有限大小的层（称为PML层）替换无限域来处理人工边界，相比无限单元法，具有计算代价更小的优势。不过，由于提出的2.5D MEF-PML中空间域的物理空间坐标是基于笛卡尔坐标，因此，该模型只适用于处理直线型观察域，不能处理曲线型观察域。

4）目前的数值模型是由之前的2.5D MEF-MEC模型发展而来，能够模拟移动载荷及列车与线路相互作用的问题，从这个角度看，2.5DMEF-PML是解决隧道内轨道交通诱发振动问题的综合数值工具。

参考文献

GUPTA S, HUSSEIN M F M, DEGRANDE G, et al. A comparison of two numerical models for the prediction of vibrations from underground railway traffic [J]. Soil Dynamics and Earthquake Engineering, 2007, 27(7): 608-624. [百度学术]

徐斌, 雷晓燕, 徐满清, 等. 饱和土体中空沟对移动荷载被动隔振的2.5D边界元法分析[J]. 岩土力学, 2012, 33(4): 1079-1086, 1102. [百度学术]

XU B, LEI X Y, XU M Q, et al. Analysis of effectiveness of passive isolation for vibration due to moving loads on saturated soil by using open trench with 2.5D boundary element method [J]. Rock and Soil Mechanics, 2012, 33(4): 1079-1086, 1102. (in Chinese) [百度学术]

高广运, 姚哨峰, 孙雨明, 等. 2.5维有限元分析高铁荷载诱发非饱和土地面振动[J]. 同济大学学报(自然科学版), 2019, 47(7): 957-966. [百度学术]

GAO G Y, YAO S F, SUN Y M, et al. Unsaturated ground vibration induced by high-speed train loads based on 2.5D finite element method [J]. Journal of Tongji University (Natural Science), 2019, 47(7): 957-966. (in Chinese) [百度学术]

LOMBAERT G, DEGRANDE G. Ground-borne vibration due to static and dynamic axle loads of InterCity and high-speed trains [J]. Journal of Sound and Vibration, 2009, 319(3/4/5): 1036-1066. [百度学术]

罗振东. 混合有限元法基础及其应用[M]. 北京: 科学出版社, 2006. [百度学术]

LUO Z D. Foundation and application of hybrid finite element method [M]. Beijing: Science Press, 2006. (in Chinese) [百度学术]

王志军, 李先枝. 非线性Sobolev方程的非协调H¹-Galerkin混合有限元方法[J]. 数学的实践与认识, 2018, 48(12): 193-199. [百度学术]

WANG Z J, LI X Z. A nonconforming H¹-Galerkin mixed finite element method for nonlinear Sobolev equations [J]. Mathematics in Practice and Theory, 2018, 48(12): 193-199. (in Chinese) [百度学术]

AMADO-MENDES P, GODINHO L M C, COSTA P A. 2.5D modeling of soil-structure interactionusing a coupled MFS-FEM formulation [C]// 11th World Congress on Computational Mechanics (WCCM XI), 2014: Barcelona. [百度学术]

朱志辉, 刘禹兵, 王力东, 等. 基于2.5维有限元法和虚拟激励法的地铁交通场地随机振动分析[J]. 中国铁道科学, 2020, 41(4): 29-39. [百度学术]

ZHU Z H, LIU Y B, WANG L D, et al. Analysis of ground random vibration induced by subway transit based on 2.5-dimensional finite element and pseudo excitation methods [J]. China Railway Science, 2020, 41(4): 29-39. (in Chinese) [百度学术]

JOHNSON S G. Notes on perfectly matched layers (PMLs) [J]. ArXiv Eprints, 2021: ArXiv: 2108. 0538. [百度学术]

邓琴.弹塑性边界元法的若干基础研究及应用[D].北京：中国科学院大学, 2011. [百度学术]

DENG Q. Some basic research and application of elastoplastic boundary element method [D]. Beijing: University of Chinese Academy of Sciences, 2011. [百度学术]

申光宪. 边界元法[M]. 北京: 机械工业出版社, 1998. [百度学术]

SHEN G X. Boundary meta-method[M]. Beijing: China Machine Press, 1998. (in Chinese) [百度学术]

KAUSEL E, DE OLIVEIRA BARBOSA J M. PMLs: A direct approach [J]. International Journal for Numerical Methods in Engineering, 2012, 90(3): 343-352. [百度学术]

FORREST J A, HUNT H E M. A three-dimensional tunnel model for calculation of train-induced ground vibration [J]. Journal of Sound and Vibration, 2006, 294(4/5): 678-705. [百度学术]

YU H T, CAI C, BOBET A, et al. Analytical solution for longitudinal bending stiffness of shield tunnels [J]. Tunnelling and Underground Space Technology, 2019, 83: 27-34. [百度学术]

FACCIORUSSO J, MADIAI C. On cohesive soil damping estimation by free vibration method in resonant column test [J]. Geotechnical Testing Journal, 2020, 43(6): 20180241. [百度学术]

SALEH ASHEGHABADI M, RAHGOZAR M A. Finite element seismic analysis of soil-tunnel interactions in clay soils [J]. Iranian Journal of Science and Technology, Transactions of Civil Engineering, 2019, 43(4): 835-849. [百度学术]

REN X W, WU J F, TANG Y Q, et al. Propagation and attenuation characteristics of the vibration in soft soil foundations induced by high-speed trains [J]. Soil Dynamics and Earthquake Engineering, 2019, 117: 374-383. [百度学术]

基于2.5维有限元-完全匹配层技术的隧道振动建模研究 PDF

摘要

关键词

1 2.5维MEF-PML法数值模型的概念和计算公式

1.1 概况

1.2 2.5D MEF-PML结合法公式

2 理论验证

2.1 不确定介质中的垂直载荷

2.2 不确定介质中的隧道

3 隧道振动的参数

3.1 参考方案说明

3.2 隧道深度的影响

3.3 隧道衬砌厚度的影响

3.4 土体物质阻尼的影响

4 结论

参考文献