摘要 : 引力位是地球自由振荡数值模拟中不可或缺的一部分，也是重力异常研究中需要计算的对象。由于引力位满足无界的泊松-拉普拉斯方程，其在无穷远边界处为零，这对数值模拟造成了困扰。针对此，采用无限元-谱元混合法，直接对无穷远边界进行拟合，不再对边界条件做近似；同时考虑到三维地球模型计算效率问题，采用2.5维控制方程；最后，通过数值试验，验证了方法的准确性。 Application of infinite-spectral hybrid method in 2.5-dimensional gravitational potential calculation Abstract : Gravitational potential is an inseparable part of the earth's free oscillation simulation. It is also the object of calculation in the study of gravity anomalies. Since the gravitational potential satisfies the unbounded Poisson-Laplace equation, which is zero at infinity, this is frustrating for numerical simulations. The purpose of numerical simulation is usually achieved by limiting the solution domain and approximating its boundary conditions. We attempt to use the infinite-spectral hybrid method to fit the infinity boundary directly and no longer approximate the boundary conditions. Considering the computational efficiency of the 3D Earth model, this study uses a 2.5-dimensional governing equation. Finally, numerical experiments verify the factual accuracy of this method.

引力位从地球内部到地球外部无处不在。当大地震发生时，地球内部物质产生位移，引起引力位扰动。因此，在地球自由振荡以及低频大尺度(如全球尺度)地震波场数值模拟中，引力位计算是一个不可避免的问题。但无论引力位还是引力位扰动，在地球内部均满足泊松方程，在地球外满足拉普拉斯方程，其边界条件均为在无穷远处衰减为零 ^{[ 1 ]} ，这迫使我们需要利用数值方法求解一个无界微分方程。长期以来，为了避免求解该无界微分方程而采用Cowling近似 ^{[ 2 ]} ，舍弃掉引力位扰动项。这会造成模拟波形的不准确，从而导致精细频谱分析困难，特别是对于地球自由振荡 ^{[ 3 ]} 、固体潮汐 ^{[ 4 ]} 等问题。

针对此类问题，常规的处理方法是将求解域设定在有限的范围内。例如将地球内部网格向外延伸到足够远的地球外(如数倍地球半径处)，然后在这一人为指定的边界处设置边界条件，只要延伸得足够远，则对地球内部及近地表的引力位影响较小。该边界条件的选取有多种方法：例如Cai和Wang ^{[ 5 ]} 提出的 r ^-4 近似假设，用于近似表示位置 r 处的函数值；Chaljub和Valette ^{[ 6 ]} 采用球谐函数展开构建DtN算子，用于该边界条件的处理；Xu和Huo ^{[ 7 ]} 使用Legendre多项式近似表示边界处的值。这些处理总体上都是通过对边界条件取值的近似来满足所需的计算精度。事实上，对于无界泊松问题，数学上很早就有研究 ^{[ 8 ]} ，只是一些方法没有在地学中得到使用。

本文在延拓网格的基础上，运用无限元映射 ^{[ 9 ]} 的思想，直接通过构造插值多项式，自然地满足无穷远趋于零的边界条件，从而避免了对边界条件取值的近似处理，因而结果更准确。

本文采用谱元法 ^{[ 10 ]} 作为数值求解方法。谱元法是一种综合型方法，既保留了传统有限元方法处理复杂网格和边界的便利性，又确保了类似谱配置方法的高精度。考虑到全球尺度模拟的计算量和模型的复杂度(如S40RTS ^{[ 11 ]} 地幔模型)，本文采用2.5维 ^{[ 12 - 13 ]} 的控制方程。在2.5维假设下，引力位的计算简化到二维D形区域。由于对称轴的原因，轴单元的数值积分需要进行改进。本文采用Gerritsma和Phillips ^{[ 14 ]} 提出的归一化方法，在轴单元引入GRL(Gauss-Radau-Legendre)数值积分。

为分析此模拟方法的有效性，本文开展了数值试验，以三维球对称实心球为例，数值模拟其引力位，并与解析解做对比。

1 方法 1.1 控制方程

对于引力位 Ψ 以及引力位扰动 ϕ ，均满足泊松-拉普拉斯方程 ^{[ 1 ]} ：

$\nabla^2 \mathit{\Psi}= \begin{cases}4 {\rm{ \mathsf{ π} }} g \rho, & \text { in V }, \\ 0, & \text { outside V }.\end{cases}$ $\nabla^2 \phi(d)= \begin{cases}-4 {\rm{ \mathsf{ π} }} g \nabla \cdot(\rho d), & \text { in V, } \\ 0, & \text { outside V }.\end{cases}$ $\int_{\mathit{\Omega}} \nabla u \cdot \nabla \psi \mathrm{d} V=\int_{\mathit{\Omega}} s \cdot \psi \mathrm{d} V.$ $x(\boldsymbol{\xi})=\sum\limits_{\alpha=1}^{n_\alpha} \boldsymbol{x}_\alpha N_\alpha(\boldsymbol{\xi}) .$ $f(x(\xi, \eta)) \approx \sum\limits_{\alpha=0}^N \sum\limits_{\beta=0}^N f\left(\xi_\alpha, \eta_\beta\right) \ell_\alpha(\xi) \ell_\beta(\eta).$ $\left\{\begin{array}{l} \frac{\partial f}{\partial \xi}(x(\xi, \eta))=\sum\limits_{\alpha, \beta=0}^N f^{\alpha \beta} \ell^{\prime}{ }_\alpha(\xi) \ell_\beta(\eta), \\ \frac{\partial f}{\partial \eta}(x(\xi, \eta))=\sum\limits_{\alpha, \beta=0}^N f^{\alpha \beta} \ell_\alpha(\xi) \ell^{\prime}{ }_\beta(\eta) . \end{array}\right.$ $\int_{\mathit{\Omega}_e} f(x(\xi, \eta)) \mathrm{d} \xi \mathrm{d} \eta \approx \sum\limits_{i, j=0}^N \omega_i \omega_j f^{i j} .$ $\boldsymbol{K}_e=\int_{\mathit{\Omega}_e} \nabla \psi \cdot \nabla u r \mathrm{~d} r \mathrm{~d} z=\int_{\mathit{\Omega}_e} \nabla \psi \cdot \nabla u\left|\mathcal{J}_e\right| r \mathrm{~d} \xi \mathrm{d} {\eta} .$ $\boldsymbol{F}_e=\int_{\mathit{\Omega}_e} s \cdot \psi r \mathrm{~d} r \mathrm{~d} z=\int_{\mathit{\Omega}_e} s \cdot \psi\left|\mathcal{J}_e\right| r \mathrm{~d} \xi \mathrm{d} \eta .$ $\mathcal{J}_e=\left\|\begin{array}{ll} \frac{\mathrm{d} r}{\mathrm{~d} \xi} & \frac{\mathrm{d} z}{\mathrm{~d} \xi} \\ \frac{\mathrm{d} r}{\mathrm{~d} \eta} & \frac{\mathrm{d} z}{\mathrm{~d} \eta} \end{array}\right\|.$ $\boldsymbol{K}_e=\int_{\mathit{\Omega}_e} \nabla u \cdot \nabla \psi \frac{r}{\omega} \omega \mathrm{d} r \mathrm{~d} z .$ $\lim\limits_{\xi \rightarrow-1} \frac{r}{\omega}=\lim\limits_{\xi \rightarrow-1} \frac{r(\xi, \eta)}{1+\xi}=\left.\frac{\partial r}{\partial \xi}\right|_{\xi=-1} .$ $\tilde{r}= \begin{cases}\frac{r}{\omega}, & \text { if } r>0, \\ \frac{\partial r}{\partial \xi}, & \text { if } r=0 .\end{cases}$ $\mathit{\Psi}= \begin{cases}2 {\rm{ \mathsf{ π} }} g \rho\left(R^2-\frac{1}{3} r^2\right), & r \leqslant R, \\ \frac{4}{3} \frac{{\rm{ \mathsf{ π} }} \rho g R^3}{r}, & r>R .\end{cases}$ 3 结论与讨论

通过数值试验，验证了本文提出的无限元-谱元混合法在引力位计算模拟中的真实准确性。只需要有限的几层网格延拓，并不需要真的把求解域延伸到无穷远，即可实现对无穷远边界条件的把控。同时引入的2.5维近似假设可以很好地解决3维真实地球模型计算的复杂度。对于传统的PREM模型 ^{[ 12 ]} 可以直接应用此方法。对于更复杂的地球模型(如S40RTS ^{[ 11 ]} )，可以考虑保留其有效的球谐展开阶次，利用傅里叶变换转化为本文的2.5维问题 ^{[ 13 ]} 。最后，本文通过无限元-谱元混合法的引入，弥补了传统引力位边界条件近似处理的缺陷，为地球自由振荡和低频大尺度(如全球尺度)地震波场数值模拟提供有效支持。

Cowling T G. The non-radial oscillations of polytropic stars[J]. Monthly Notices of the Royal Astronomical Society, 1941, 101(8): 367-375. Doi:10.1093/mnras/101.8.367 Yan R, Woith H, Wang R J, et al. Earth's free oscillations excited by the 2011 Tohoku Mw 9.0 earthquake detected with a groundwater level array in mainland China[J]. Geophysical Journal International, 2016, 206(3): 1457-1466. Doi:10.1093/gji/ggw213 Dahlen F A, Tromp J. Theoretical global seismology[M]. Princeton: Princeton University Press, 2021. Doi:10.1515/9780691216157 Cai Y E, Wang C Y. Fast finite-element calculation of gravity anomaly in complex geological regions[J]. Geophysical Journal International, 2005, 162(3): 696-708. Doi:10.1111/j.1365-246X.2005.02711.x Chaljub E, Valette B. Spectral element modelling of three-dimensional wave propagation in a self-gravitating Earth with an arbitrarily stratified outer core[J]. Geophysical Journal International, 2004, 158(1): 131-141. Doi:10.1111/j.1365-246X.2004.02267.x Xu C Y, Huo Z J. Numerical simulation of gravity anomaly based on the unstructured element grid and finite element method[J]. Mathematical Problems in Engineering, 2020, 2020: 3604084. Doi:10.1155/2020/3604084 Hejlesen M M, Rasmussen J T, Chatelain P, et al. A high order solver for the unbounded Poisson equation[J]. Journal of Computational Physics, 2013, 252: 458-467. Doi:10.1016/j.jcp.2013.05.050 O.C Z, R.L T, Zhu J Z. The finite element method: Its basis and fundamentals[M]. Sixth. Burlington (Mass.): Elsevier, 2005. Chaljub E, Komatitsch D, Vilotte J P, et al. Spectral-element analysis in seismology[J]. Advances in Geophysics. Elsevier, 2007, 48: 365-419. Doi:10.1016/S0065-2687(06)48007-9 Ritsema J, Deuss A, van Heijst H J, et al. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements[J]. Geophysical Journal International, 2011, 184(3): 1223-1236. Doi:10.1111/j.1365-246X.2010.04884.x Nissen-Meyer T, Fournier A, Dahlen F A. A two-dimensional spectral-element method for computing spherical-earth seismograms-Ⅰ. Moment-tensor source[J]. Geophysical Journal International, 2007, 168(3): 1067-1092. Doi:10.1111/j.1365-246X.2006.03121.x Nissen-Meyer T, Fournier A, Dahlen F A. A 2-D spectral-element method for computing spherical-earth seismograms-Ⅱ. Waves in solid-fluid media[J]. Geophysical Journal International, 2008, 174(3): 873-888. Doi:10.1111/j.1365-246X.2008.03813.x Gerritsma M I, Phillips T N. Spectral element methods for axisymmetric stokes problems[J]. Journal of Computational Physics, 2000, 164(1): 81-103. Doi:10.1006/jcph.2000.6574 Shen J, Tang T. Spectral and high-order methods with applications[M]. Beijing: Science Press, 2006. 李瑞浩. 重力学引论[M]. 北京: 地震出版社, 1988.