时间域边界元方法中对数空间奇点的解析积分处理
赵峰,
李晓坤,
罗俊成,
雷卫东,
李洪军
《Mathematical and Computational Applications》:Analytical Integration for Logarithmic Spatial Singularities in the Time Domain Boundary Element Method
Feng Zhao,
Xiaokun Li,
Juncheng Luo,
Weidong Lei and
Hongjun Li
【字体:
大
中
小
】
时间:2026年05月01日
来源:Mathematical and Computational Applications 2.1
编辑推荐:
摘要 对数空间奇异积分的处理是影响时间域边界元方法(TD-BEM)求解弹性动力学问题结果可靠性的关键挑战。为了解决这个问题,本文基于Heaviside函数的基本性质,推导出了一组分析上严谨的对数空间奇异积分公式,使得在TD-BEM中可以直接进行时空分析求解
摘要 对数空间奇异积分的处理是影响时间域边界元方法(TD-BEM)求解弹性动力学问题结果可靠性的关键挑战。为了解决这个问题,本文基于Heaviside函数的基本性质,推导出了一组分析上严谨的对数空间奇异积分公式,使得在TD-BEM中可以直接进行时空分析求解这类奇异积分。这些公式填补了TD-BEM中弹性动力学问题中对数空间奇异积分缺乏直接分析解的研究空白,并丰富了TD-BEM处理奇异积分的理论体系。选择了三个典型的弹性动力学工程问题进行数值验证,包括在均匀突然载荷作用下的固定-固定梁、在边界爆破载荷作用下的无限域带单个空腔的模型,以及位于山谷地形下的双隧道受到地铁振动载荷影响的模型。将所提方法的计算结果与参考解进行比较,结果表明计算结果与参考解吻合良好,有效验证了这些分析积分公式的正确性和工程适用性。
1. 引言
随着科学计算和工程仿真需求的增长,时间域边界元方法(TD-BEM)作为一种高效的数值分析方法,在求解波动传播和瞬态动力学问题方面展现了独特的优势。与有限元方法(FEM)和有限差分方法(FDM)相比,边界元方法(BEM)具有降维特性,显著减少了计算自由度。此外,TD-BEM的边界积分方程是基于完整的时空基本解构建的,这使得该方法天然适合处理涉及无限或半无限边界的问题。然而,由基本解的奇异性引起的奇异积分的处理仍然是TD-BEM数值实现中的关键挑战。这些奇异性不仅会降低数值解的可靠性,还可能导致不稳定的甚至错误的解[1,2],从而限制了TD-BEM的进一步应用和普及。
在TD-BEM框架内,奇异性主要来源于两个来源。首先是波前奇异性,表现为由于应力波的到达导致波前附近的位移和应力场发生突变。其次是强奇异或对数空间奇异性,当单位脉冲作用在无限小面积上时,在源点产生无限应力。这两种奇异性的存在被称为双重奇异性。在过去的几十年中,许多学者提出了各种方法来解决TD-BEM中的奇异积分问题。例如,Guiggiani等人采用了局部坐标展开来处理奇异积分[3,4],而Niu和Zhou应用正则化技术来解决势问题中的强奇异积分[5]。Xie等人提出了sinh变换方法,用于求解薄壁结构BEM分析中的奇异和近奇异积分问题[6,7]。Giannopoulos和Anifantis[8]为基于BEM的界面现象建模和奇异行为处理提供了有价值的见解。Mansur和Brebbia率先使用Hadamard主值积分来处理波前奇异性[9,10,11,12,13,14],这有效提高了计算精度和效率,从而显著推动了TD-BEM的发展。对于空间奇异性,曾经广泛使用刚体位移方法,结合动态和静态奇异性特性进行处理。然而,这种方法是一种间接的数值方法,可能导致计算误差的累积。为了实现弹性动力学中空间奇异性的直接求解,Lei等人[15,16]、Xie等人[17]和Zhong等人[18]开发了独立于刚体位移方法的分析积分方法。然而,这些分析方法仅适用于强奇异性,不适用于对数空间奇异积分的积分。在弹塑性静态分析领域,常用常应变场方法[19,20]和初始应力展开技术[21]来处理空间奇异性。正如Telles等人[19,20]和Banerjee等人[21]所指出的,这两种方法需要全域离散化,这损害了BEM在解决无限域问题时的固有优势。Deng等人[22]将互补理论引入弹塑性边界积分方程,将域积分转换为边界积分,并将域奇异性转化为边界奇异性,从而实现了强奇异域积分的求解。这种方法保留了BEM的固有优势,但仅适用于弹塑性静态问题。对于弹塑性动态问题,Li等人[23,24]应用分析积分直接处理TD-BEM中的奇异性,获得了稳定且准确的计算结果,其中只需离散化边界和局部塑性区域。不幸的是,这些分析方法仍然不适用于对数空间奇异性。在早期,Johnston和Elliott[25]提出了Telles方法的通用形式[26],通过任意奇数阶多项式坐标变换优化了BEM中对数奇异积分的数值解,显著提高了计算精度和稳定性。Smith[27]为等参元素(线性、二次和三次元素)开发了直接的高斯求积公式,使得无需分离奇异项即可高效计算对数奇异性。这两种方法都对该领域做出了重要贡献。然而,这两种方法本质上仍然依赖于高斯求积进行数值积分,并没有建立完全分析的解决方案路径。
作为TD-BEM数值实现的核心基石,完善奇异积分处理的理论体系至关重要,因为它是确保计算结果可靠性和准确性的关键前提。本文首先研究了TD-BEM中対数空间奇异性的产生机制及其对数值解的影响。在此基础上,利用Heaviside函数的基本性质,提出了一种直接分析积分方法来求解対数空间奇异积分。推导出了一组TD-BEM中对数空间奇异积分的解析公式。与数值策略[26,27]不同,所提出的方法完全避免了数值积分误差,且不依赖于网格密度和时间步长的细化,从而在动态分析中提供了更高的稳定性和可靠性。设计了三个涵盖有限域、无限域和半无限域的典型数值示例来验证所提公式的正确性和适用性。通过一个无限域示例也验证了结合所提公式的TD-BEM的鲁棒性。
2. 位移边界积分方程
本文提出的公式适用于没有体力作用的平面应变弹性动力学问题。一旦将弹性模量E和泊松比ν替换为(1 + 2ν)E/(1 + ν)2和ν/(1 + ν),就可以得到平面应力问题的公式。笛卡尔坐标系O ? x1x2中的位移边界积分方程如方程(1)所示,该方程也在[11]中给出。
(1)
其中P和τ分别代表单位脉冲的源点和激励时刻;Q和t分别代表场点和观测时刻;cik代表源点P的位置系数,如方程(2)所示;
δik是克罗内克德尔塔;φ表示弹性静力学中的表面牵引力[20];Γε是以点P为中心、半径为ε的圆弧,如图1所示;θ1和θ2分别是x1轴正方向与点P两侧切线外法线之间的角度,逆时针方向为正值,顺时针方向为负值。
(2)
vd和vs分别表示P波(主波)和S波(剪切波)的速度;r表示源点P与场点Q之间的距离;方程(3)和(4)中的其他参数如下:
μ和λ是拉梅常数,φ = λ/2μ;n是内部域边界的外向单位法线;下标i、k、l = 1、2分别对应x1和x2方向;w用“d”和“s”表示P波和S波。
图1. P位于边界的情况。
3. 数值实现
弹性动力学的边界元积分方程需要在时间和空间上进行离散化。网格长度和时间步长是根据Mansur和Brebbia[10]的研究结果确定的。
3.1. 时间离散化
时间域[0, t]被划分为M个间隔,时间步长为Δt = t/M。时间节点定义为tm = mΔt,其中m = 0, 1, 2, …, M。在每个间隔内,假设表面牵引力是恒定的,而位移是线性变化的。这种假设适用于具有时间跳变的载荷[28,29],当时间步长较小时也适用于没有时间跳变的载荷[15]。场点Q在间隔[tm?1, tm]内的表面牵引力和位移可以使用时间节点变量通过插值方法表示,如方程(5)所示。
(5)
其中
(6)
3.2. 空间离散化
假设边界被离散为ne个线性边界元素,对于边界元素e内的场点Q,表面牵引力和位移在元素内变化。它们可以用元素节点变量表示。
(7)
其中N1(ξ)和N2(ξ)是形状函数,如方程(8)所示。
(8)
3.3. 边界元积分方程的离散化
方程(1)中积分的离散表达式详细见方程(9)。
(9)
其中
和分别表示场点Q在时刻tm对源点P在时刻t的位移的影响系数,如方程(10)所示。
(10)
4. 元素影响系数的求解
影响系数已在参考文献[16]中计算;因此,本研究仅求解。
中的奇异性来源于两个主要来源。首先是波前奇异性,其数学形式为[vw(t ? τ) ? r]?n/2 → ∞。这种奇异性对应于波前附近响应的突变,由应力波的到达引起。其次是対数空间奇异性,表示为lnr → ∞,这是由于单位脉冲作用在无限小面积上产生的。波前奇异性和对数空间奇异性都发生在节点与源点P重合的元素中,而在其他元素中只出现波前奇异性。
根据影响系数积分中包含的奇异性类型,采用不同的积分方案。对于仅包含波前奇异性的积分,采用Hadamard主值积分进行时间积分,高斯求积进行空间积分。对于包含对数空间奇异性的积分,直接分析积分方法用于时间和空间积分。对于非奇异积分,采用Riemann积分进行时间积分,高斯求积进行空间积分。
4.1. 无空间奇异性的影响系数
当源点P与场点Q不重合时,我们有r ≠ 0,如图2所示。在这种情况下,不存在对数空间奇异性;即lnr保持有限。只出现波前奇异性,表示为[vw(t ? τ) ? r]?n/2 → ∞。
图2. 源点P与元素e的场点Q不重合。
4.1.1. 核函数的表达式
可以通过将方程(3)代入方程(10)得到的形式。
(11)
其中
和分别是与时间和空间相关的核函数。
4.1.2. 核函数的时间积分
引入系数Aw1~Aw6来简化表达式,这些系数依赖于M ? m。根据Heaviside函数的性质,奇异时间trw(trw = t ? r/vw)与时间间隔[tm?1, tm]之间的关系可以分为以下四种类型。核函数可以相应地计算得出。对于 和 ,有 和 ;对于 和 ,有 和 ;对于 和 ,有 和 。4.1.3. 影响系数的空间积分将核函数代入方程(11)后,使用高斯求积法进行影响系数的空间积分。计算公式如方程(12)所示。其中 ξg 和 ωg 分别代表高斯求积点和权重;ng 是高斯积分点的数量;r(ξg) 表示场点 Q 与对应于 ξg 的源点 P 之间的距离,如方程(13)所示。4.2. 具有对数空间奇异性的影响系数当点 P 与元素 e 的节点重合时,距离 r 从 0 变化到 Le。图 3 展示了 P 与元素的首个节点重合的情况。当 r → 0 时,会出现对数空间奇异性,即 lnr → ∞。图 3. 源点 P 与元素 e 的首个节点重合。4.2.1. 奇异积分部分的分离在图 3 所示的条件下,局部坐标 ξ 的表达式可以推导为方程(14)。将方程(14)代入形状函数方程(8),得到:将方程(15)和(16)代入方程(11),影响系数变为:在 中没有空间奇异性,而在 中存在对数空间奇异性。第 4.1 节中介绍的方法可以用来处理 ,而 需要分解为非奇异部分 和对数奇异部分。由于奇异积分具有不同的解特性,因此可以将其分解为奇异部分和非奇异部分。非奇异部分可以通过高斯求积等常规数值积分方法直接计算,而奇异部分不能通过普通数值求积方法评估,因此需要专门的处理方法。在这项工作中,提出了一种用于奇异部分的解析解方法。 的表达式为 (19)。当源点 P 与元素 e 的第二个节点重合时,也可以得到类似的结果。一旦解决了奇异部分 ,就可以得到所有的影响系数。4.2.2. 奇异积分部分的计算时空积分域和核函数图 4 中用点填充的区域是核函数的时空积分域,其数学形式由方程(20)给出。图 4. 三种类型的时空积分域。在图 4 中,定义了 、 和 。积分域 Dτr 可以分为以下三种类型:矩形域 ①:;梯形域 ②:;混合域 ③:。这三种类型的积分域都适用于先对 r 进行积分,然后再对 τ 进行积分的顺序。通过交换方程(19)中的积分顺序,我们得到 (21)。每种积分域的解决方法如下。对于矩形域和梯形域,首先确定 τ 的积分范围为 tm?1 到 tm。在每个 τ 时间间隔内,r 的下限固定在 0,上限对于矩形域为 Le,对于梯形域为 vw(t ? τ)。完成对 r 的积分后,依次对 τ 进行积分。对于混合积分域,首先将整个区域划分为独立的矩形域和梯形域子域,tLw 作为边界。每个子域分别按照矩形域和梯形域的相应积分规则进行计算。就奇异性特征而言,只有矩形域存在对数空间奇异性。相比之下,梯形域和混合域同时包含对数空间奇异性和波前奇异性,导致更强的奇异行为。这种高阶奇异性对数值结果的稳定性和准确性有更显著的影响。因此,适当处理梯形域和混合域对于确保 TD-BEM 解的可靠性至关重要。2. 核函数的空间积分核函数的空间积分根据积分域的类型,使用解析积分方法进行评估。对于 r2 ≥ Le 的情况,Dτr 成为矩形域 ①,其中 r ∈ [0, Le], 和 分别表示 bw 的上下积分限,其表达式如下。由于只需考虑 bw 的上积分限,进一步推导得到:对于 r1 ≤ Le 的情况,Dτr 成为梯形域 ②,其中 r ∈ [0, vw(t ? τ)];对于 r2 < Le 且 r1 > Le 的情况,Dτr 成为混合域 ③。对 τ ∈ [tm?1, tLw) 的积分可以使用矩形域的方法进行评估,对 τ ∈ [tLw, tm] 的积分可以使用梯形域的方法进行评估。3. 核函数的时间积分在获得空间积分结果后,可以对 τ ∈ [tm?1, tm] 的核函数进行时间积分。引入系数 Bw1~Bw6 来简化表达式,这些系数取决于 M ? m,具体如下。对于矩形域 ①,其中 r ∈ [0, Le];对于梯形域 ②,其中 r1 ≤ Le;对于混合域 ③,其中 r2 < Le 且 r1 > Le。对于 τ ∈ [t1, tLw),可以通过在情况 ① 的公式中将上限 t2 替换为 tLw 来简单地进行时间积分;对于 τ ∈ [tLw, t2],可以通过在情况 ② 的公式中将下限 t1 替换为 tLw 来进行时间积分。可以通过将核函数和系数 、 和 代入方程(21)来获得 。5. 组装和求解5.1. 组装通过依次将源点 P 与每个边界节点对齐,可以确定所有元素的元素影响系数。在时间和空间上组装元素影响系数后,形成与边界积分方程对应的代数方程组。代数方程组的矩阵形式如方程(22)所示。其中 HMM 是位移影响系数矩阵,GMM 是表面牵引影响系数矩阵;uM 是所有节点的节点位移列向量,pM 是节点表面牵引列向量;列向量 aM 仅与前 M ? 1 个时间步的响应相关。5.2. 代数方程组的求解uM 和 pM 中可能存在未知量。因此,需要将方程(22)重新排列成标准形式,将未知量和已知量分别放在方程组的左侧和右侧。求解过程如下。代数方程组的标准形式将 uM 中的已知量和 pM 中的未知量交换位置,同时交换 HMM 和 GMM 的相应列,并改变它们的符号。然后,uM 和 pM 分别转换为列向量 xM(包含所有边界未知量)和列向量 (包含所有已知量);同时 HMM 和 GMM 分别转换为 和 。设 ,则方程(22)的标准形式为方程(24)。2. 边界节点的未知位移和表面牵引可以通过方程(25)求解未知向量 ,然后结合边界条件,可以得到所有边界节点的位移和表面牵引。3. 内部点的位移和应力在求解所有边界节点的位移和表面牵引后,可以通过内部点边界积分方程(见 [30])得到内部点的位移和应力。由于内部点无法与边界节点重合,内部点边界积分方程不包含空间奇异性。6. 验证选择了三个典型的弹性动力学问题进行验证,即固定-固定梁在均匀施加的突然载荷下、具有单个空腔的无限域在其边界受到爆破载荷的作用,以及位于山谷地形下的双隧道受到地铁振动载荷的作用。这些问题涉及有限域、无限域和半无限域的典型场景,以及简单和复杂的几何边界配置。由于所提出的对数奇异性处理方法仅应用于源点与元素节点重合的少数奇异元素,而大多数非奇异元素与传统 TD-BEM 相同,因此额外的计算成本可以忽略不计。因此,不需要单独进行效率研究。通过将结果与从特征线方法 [31] 和 FEM 获得的参考解进行比较,验证了所提出方法的鲁棒性和正确性。对于有解析解的情况,使用特征线方法;当没有解析解时,使用 FEM 结果作为参考。6.1. 鲁棒性分析为了验证所提出方法的数值鲁棒性,建立了一个具有单个空腔的无限域模型,并在空腔边界施加爆破载荷,如图 5 所示。空腔的半径 r0 = 1 m,爆破载荷 p(t) = 1.435p0 × (e?1279t ? e?12792t),其中 p0 = 50 MPa,如图 6 所示。材料的泊松比、杨氏模量和密度分别为 ν = 0.25、E = 4.0 × 1010 N/m2 和 ρ = 2.8 × 103 kg/m3。P 波速度相应地确定为 vd = 4140 m/s。空腔边界分别由 20、28、40 个线性边界元素离散化。空腔边界的空间网格离散化如图 7 所示,边界元素按顺时针方向依次离散化。时间步长设置为 5 × 10?5 s、3.5 × 10?5 s、2.5 × 10?5 s。离散化系数 β 的范围从 0.33 到 1.32 [10]。图 5. 受到均匀爆破载荷作用的无限域模型示意图。图 6. 爆破载荷的时间历史曲线。图 7. 空间网格离散化示意图。提取了 r = 2 m 的圆形轮廓处的峰值径向位移。通过以大约 0.7 的缩放因子依次减小时间步长和空间网格尺寸,分析了动态响应的变化特性。通过与特征线方法的解析解进行比较,计算了相对误差。表 1 显示了半径 r = 2 m 的圆上的峰值径向位移和相对误差。从表 1 可以看出,r = 2 m 处的峰值径向位移逐渐趋近于特征线方法的解析解,且相对误差随着空间网格的细化和时间步长的减小而持续降低。值得注意的是,随着网格尺寸的减小,数值精度显著提高,而由时间步长调整引起的收敛速率相对较慢。当总网格数从 20 增加到 40,时间步长从 5 × 10?5 s 减少到 2.5 × 10?5 s 时,最大相对误差从 2.15% 降至仅 0.009%。在参数变化过程中,没有出现数值振荡或发散,数值结果保持稳定和单调收敛。上述结果表明,所提出的 TD-BEM 算法对空间和时间离散化变化具有良好的鲁棒性。当空间网格数量为 20 且时间步长设置为 5 × 10?5 s,即 β = 0.66 时,计算相对误差已经控制在约 2% 以内,实现了令人满意的高计算精度。因此,图 8 展示了在这种离散方案下 r = 2 m 的圆形轮廓上的径向位移和应力的时间历史曲线,进一步验证了所提出算法的全时域计算精度。“特征线”标记的曲线表示通过特征线方法 [31] 获得的解析解。圆环的应力与位移曲线(半径r = 2米,每个三角形标记代表两个离散时间步长):(a) 径向应力曲线;(b) 切向应力曲线;(c) 位移曲线。6.2. 准确性分析 6.2.1. 受均匀突然载荷作用的固定-固定梁 固定-固定梁的示意图如图9所示,其中a = 1米,突然载荷p(t) = 200 MPa(t ≥ 0),如图10所示。材料的泊松比、杨氏模量和密度分别为ν = 0.3、E = 2.1 × 10^11 N/m^2和ρ = 7.9 × 10^3 kg/m^3。相应的P波速度可确定为vd = 5982 m/s。为了降低计算复杂度,采用了半结构模型,网格和相应的边界约束条件如图11所示。共使用了30个线性边界元素(每个元素长度为0.2米)来离散化边界域,时间步长设置为2.5 × 10^-5秒。离散化系数β = 0.75 [10]。图12显示了中间跨度点A的位移结果,其中u0代表点A的静态解 [32]。图9. 受突然载荷作用的固定-固定梁示意图。图10. 突然载荷的时间历史曲线。图11. 半结构的网格和约束条件。图12. 中间跨度点A的位移曲线(每个三角形标记代表两个离散时间步长)。6.2.2. 受地铁振动载荷作用的山谷地形下的双隧道 山谷地形下的双隧道示意图如图13所示,两个隧道的中心分别为O1(?8米, 0)和O2(8米, 0),半径r0 = 3米,中心角θ = π/3。地面上关键点的坐标如下(单位:米):B1(100, 16)、B2(52, 16)、B3(10, 16)、B4(9, 18)、B5(7, 18)、B6(1, 12)、B7(?4, 12)、B8(?7, 18)、B9(?9, 18)、B10(?10, 16)、B11(?52, 16)、B12(?100, 16)。地铁运行引起的振动载荷p(t)的时间历史曲线如图14所示。材料的泊松比、杨氏模量和质量密度分别为ν = 0.25、E = 2.2 × 10^7 N/m^2和ρ = 1.8 × 10^3 kg/m^3。相应的P波速度可确定为vd = 121 m/s。网格如图15所示。B1B2和B11B12方向的元素大小为3米,而B2B3和B10B11方向的元素大小为2米。山谷边界B3B10方向的网格大小从1米到1.5米不等。每个隧道的内边界被离散化为24个元素,每个元素的长度约为0.8米。半无限远场的辐射条件通过两个自适应半无限边界元素(ASIBE)[33]来模拟。时间步长设置为1.25 × 10^-2秒。离散化系数β = 0.76 [10]。图16显示了地面上点O3(0, 16)的位移。图13. 受地铁振动载荷作用的半无限域中的双圆柱形空腔示意图。图14. 地铁运行引起的振动载荷。图15. 山谷地形下的双隧道网格图。图16. 点B8的位移曲线(每个三角形标记代表四个离散时间步长):(a) 水平位移曲线;(b) 垂直位移曲线。6.2.3. 讨论 本文有针对性地选择了三种典型的动态工程问题作为验证示例,涵盖了有限域、无限域和半无限域。这些示例包括:受均匀突然载荷作用的固定-固定梁、受边界爆破载荷作用的无限域中的空腔,以及受地铁振动载荷作用的山谷地形下的双隧道。这三个示例分别对应于结构动力学、爆破动力学和城市轨道交通环境振动工程中的典型场景,可以全面验证在不同空间域类型下建立的对数空间奇异性解析积分公式的适用性和正确性。数值计算结果表明,所提出方法在三个示例中获得的结果与参考解非常吻合:在有限域的固定-固定梁示例中,中间跨度节点的位移时间历史曲线与参考解高度一致,验证了这些公式在处理有限域结构动力响应计算中的对数空间奇异性方面的准确性。在无限域中的空腔爆破示例中,距离空腔中心2米(r = 2米)的圆环上的径向应力、切向应力和位移的计算结果与参考解一致,表明所提出的公式保留了TD-BEM的固有特性,这些特性本质上满足无限域的远场辐射条件。同时,它们保证了无限域弹性动力学问题中高精度的奇异积分计算。在受地铁振动载荷作用的山谷地形下的双隧道示例中,指定地面测量点的位移时间历史结果与参考解吻合良好,证明了所提出的公式仍然可以有效解决具有复杂多腔边界半无限域工程问题中的对数空间奇异性问题,并且计算结果稳定。上述三个不同空间域类型的示例一致验证结果充分证明了本文建立的对数空间奇异性解析积分公式的数学正确性和工程适用性。这些公式可以有效处理TD-BEM中的对数空间奇异积分,为涉及不同空间域类型的TD-BEM弹性动力学问题的数值计算提供了一种可靠的奇异积分处理方法。7. 结论 为了应对TD-BEM中处理对数空间奇异积分的挑战,本文基于Heaviside函数的特性推导并建立了一套解析积分公式,实现了此类奇异积分的直接时空解析解。通过涵盖有限域、无限域和半无限域的三个典型弹性动力学示例验证了这些公式的有效性。主要结论如下:这些公式的推导过程在数学上是严谨的,为TD-BEM中的对数空间奇异积分提供了直接的解析解;这些公式从根本上避免了数值积分方法在处理奇异问题时影响准确性的因素,进一步改进了TD-BEM中奇异积分处理的理论体系,并为工程应用提供了可靠的方法;所提出的解析积分公式的稳健性、正确性和工程适用性得到了充分验证,计算结果与参考解非常吻合。这些公式可以有效处理不同空间域类型下的TD-BEM中对数空间奇异积分,适用于结构动力学、爆破动力学和城市轨道交通环境振动等各种场景下的弹性动力学问题计算。所提出的对数空间奇异性解析积分公式为TD-BEM在弹性动力学问题中的数值应用提供了重要的理论支持。在本研究中,该公式是为使用线性元素的2D线性弹性动力学问题开发的。通过推导相应的奇异积分解析表达式,该方法可以扩展到更高阶的等参元素和3D弹性动力学问题。这将是我们未来工作的主要方向。