用于牛顿噪声分析的高效有限元公式

《Computer Physics Communications》：An efficient finite element formulation for Newtonian noise analysis

　　

摘要
爱因斯坦望远镜是第三代地下引力波观测站，旨在实现前所未有的灵敏度，可达 3 Hz。由人为或自然振动源引起的在土壤中传播的波会产生密度波动，从而引起引力吸引，导致激光干涉仪的镜面产生运动，这被称为牛顿噪声。后者是通过积分测试质量周围土壤区域中由地震波场引起的密度波动来计算的。
本文提出了一种有限元公式，该公式利用高斯积分，根据定义在有限元网格上的地震波场，评估总牛顿噪声以及体积和表面贡献。支持线性和二次四面体（4 节点和 10 节点）及六面体（8 节点和 20 节点）有限元。该方法计算总体、体积和表面贡献。相应的体积和表面积分用有限元耦合矩阵表示，这些矩阵仅取决于几何形状和材料属性。这使得可以在不重新计算积分的情况下，针对不同地震波场高效评估牛顿噪声。
该公式经过验证，适用于平面 P 波和 S 波在弹性均匀全空间中传播的情况，假设镜面悬挂在球形腔体中，且波长远大于腔体半径，因此可以忽略波的散射。结果与解析解获得了极好的吻合度。对于瑞利波在均匀弹性半空间自由表面上传播，且测试质量悬挂在自由表面上方有限距离处的牛顿噪声，也报告了类似的良好吻合度。
该方法已实现在 ANNA 牛顿噪声分析工具箱中，该工具箱在 MATLAB 编程和数值计算平台上运行，并兼容开源 GNU Octave 科学编程语言；同时也提供 Python 版本。所提出的有限元框架为计算非均匀介质中的引力 - 地震耦合提供了一种物理上一致且计算高效的方法。

引言
爱因斯坦望远镜是第三代引力波探测器，计划在接下来的十年内建造 [1]。它由一个地下三角形激光干涉仪组成，旨在观测引力波——即 4D 时空连续体中的涟漪——例如由双黑洞或双中子星合并引起。激光干涉仪的镜面安装在悬挂塔中，这些塔放置在干涉仪角点的大型洞室内。
为了观测由更远距离或更大质量合并引起的引力波，LIGO 和 Virgo 等第二代引力波探测器的灵敏度和工作频率范围（10 Hz 以上）是不够的 [1]。为了将灵敏度提高一个数量级并将工作频率降低到 3 Hz，爱因斯坦望远镜最好在 seismic 安静的地区地下建造 [2]，以便最大程度地减少道路和铁路交通、风力涡轮机、采石场、工业和建筑活动等人为振动源的干扰。
地震牛顿噪声在低频段占主导地位：在土壤中传播地震波会产生密度波动，从而导致引力吸引和激光干涉仪镜面的相应运动。由于这些引力无法屏蔽，唯一的缓解措施是通过波场重构从地震测量中估计牛顿噪声，并将其从干涉仪数据中减去 [3]。
对于涉及平面谐波 P 波或 S 波在均匀线弹性全空间中传播的简化问题，已有牛顿噪声的解析解 [4]，假设波长远大于腔体尺寸，因此可以忽略波的散射。Harms [4] 还推导了一个解析解，用于计算悬挂在均匀线弹性半空间自由表面上方有限距离处的镜面在谐波瑞利波激发下的牛顿噪声。然而，在实际中，场地的地形和地质更为复杂，涉及山丘、山谷、沉积过程引起的分层，以及土壤和岩层的褶皱。地震波场由多种人为振动源引起，并受到土壤地形和分层的影响，以及干涉仪的洞室、隧道和竖井的影响。因此，需要数值方法来计算地震波场和相应的牛顿噪声。
有几种方法可用于计算由人为振动源引起的三维 (3D) 地震波场，并考虑动态土 - 结构相互作用。直接刚度公式或薄层法 [5]，如 MATLAB 的 ElastoDynamics 工具箱 [6] 中所实现的那样，用于模拟水平分层介质中由于谐波或瞬态激发引起的表面波或强制振动。道路和铁路交通引起的振动可以使用 MATLAB 中的 TRAFFIC 工具箱计算，该工具箱将半解析道路或轨道模型与分层土壤的二维半 (2.5D) 边界元 (BE) 模型耦合 [7], [8], [9]；或者也可以使用 MOTIV [10], [11], [12] 和 MEFISSTO [13] 程序包。配备无限元、吸收边界条件或完美匹配层 [14], [15], [16], [17] 的 3D 有限元公式、耦合有限元 - 边界元公式 [13], [18], [19] 以及 SPECFEM3D 和 SPEED [20] 等谱元模型，可用于解决更复杂的动态土 - 结构相互作用问题，例如由桥上运行的列车或风力涡轮机引起的振动，以及干涉仪洞室、隧道和竖井引起的波散射。
近期的数值研究通过首先模拟地震波场，然后直接数值评估由此产生的牛顿噪声体积积分来评估牛顿噪声。例如，在 [21] 中，利用直接刚度法求解水平分层介质的弹性动力学方程，随后通过高斯积分计算牛顿噪声。类似地，在 [22] 中，采用谱元模拟来建模地震波场。相应的牛顿噪声随后通过基于偶极子的体积积分进行评估。在这些方法中，牛顿噪声积分被视为依赖于特定地震实现的 后处理步骤，需要对每个地震波场进行完整的重新评估。假设地震波场是在 3D 有限元网格上计算的，且网格足够细化以准确解析波传播，本文提出了一种有效的有限元公式，通过高斯积分数值积分其总体、体积和表面贡献来评估牛顿噪声。该公式产生取决于几何形状的矩阵，这些矩阵仅组装一次，随后应用于任何地震波场。由此产生的框架能够在非均匀三维介质中进行计算高效的牛顿噪声分析，同时保持解析测试质量附近引力波动所需的数值精度。所提出的有限元公式通过上述解析公式进行了验证，考虑了均匀线弹性全空间中悬挂在球形腔体内的镜面在平面 P 波或 S 波激发下的牛顿噪声 [4]，以及均匀线弹性半空间自由表面上方有限距离处悬挂的镜面在瑞利波激发下的牛顿噪声 [4]。
本文的组织结构如下。第 2 节给出了引力势总导数的表达式，以及分解为体积和表面贡献的总牛顿噪声。第 3 节开发了一种有限元公式，用于评估定义在有限元网格上的地震波场的牛顿噪声。推导了有限元矩阵的表达式，当与地震位移向量相乘时，产生总牛顿噪声以及体积和表面贡献。第 4 节通过具有解析解的问题展示了所提出的有限元公式的验证。第 5 节总结了本文的结论。

章节片段
理论背景
地震牛顿噪声是由土壤中的密度波动引起激光干涉仪镜面的引力吸引而产生的不期望的运动（加速度），这些密度波动是由地震波场引起的。
总牛顿噪声
为了评估方程 (17) 中的体积积分，将控制体积 Ω? 离散化为 ne 个有限元 (FE)，体积为 Ω?e 且密度 ρe 恒定，从而得出以下求和：
δa_t(x?, ω) ? Σ (从 e=1 到 ne) G ∫ (在 Ω?e 上) ρe / ||x - x?||3 * (I? - 3 e_r ? e_r) * ?_e(x, ω) dv
其中 ?_e(x, ω) 是单元位移。
局部坐标系中的单元位移 ?_e(ξ, ω) 使用形函数 N_e^k(ξ) 进行近似：
?_e(ξ, ω) ? Σ (从 k=1 到 n) N_e^k(ξ) * ū_e^k = N_e(ξ) * ū_e
其中 ū_e = {ū_e^1T, ū_e^2T, ?, ū_e^nT}T 是一个列向量

验证
所得的数值公式使用现有的牛顿噪声解析解 [4] 进行了验证。考虑了三种配置。首先，检查平面谐波 P 波和 S 波在具有半径 r? 的腔体的均匀线弹性全空间中传播的情况。接下来，考虑瑞利波在均匀线弹性半空间表面上传播，并计算距离自由表面高度为 h 处的测试质量上的牛顿噪声。

结论
本文提出了一种有效的有限元公式，用于从 3D 地震波场计算腔体内镜面上的牛顿噪声。使用 3D 有限元网格上的高斯积分数值评估定义总牛顿噪声的体积和表面积分，以及体积和表面贡献。由于有限元矩阵仅取决于问题域的几何形状和密度，它们只计算一次。

CRediT 作者贡献声明
P. Reumers：写作 – 审查与编辑，写作 – 原始草稿，软件，方法论，调查，形式分析，概念化。X. Kuci：写作 – 审查与编辑，写作 – 原始草稿，可视化，验证，软件，方法论，调查，形式分析。S. Fran?ois：写作 – 审查与编辑，监督，软件，资源，项目管理，方法论，资金获取，形式分析，概念化。G. Degrande：写作 – 审查与编辑，

利益冲突声明
作者声明以下财务利益/个人关系，这可能被视为潜在的利益冲突：Xhorxha Kuci 和 Pieter Reumers 和 Stijn Francois 和 Geert Degrande 报告获得了佛兰德斯研究基金会的财务支持。如果有其他作者，他们声明没有已知的竞争财务利益或个人关系，这些关系可能影响本文报告的工作。

P. Reumers|X. Kuci|S. Fran?ois|G. Degrande
荷语鲁汶大学，土木工程系，结构力学系，Kasteelpark Arenberg 40, Leuven B-3001, 比利时