极坐标系下起伏地形的嵌套网格有限差分算法

地震波数值模拟已经成为人们了解地震波在复杂地球介质中传播的重要方法。有限差分因为效率高，方法简单而被广泛的使用，是地震波数值模拟的常见算法之一，基于一阶速度应力方程组的有限差分算法已经被使用于各种模型的计算当中。然而，当对包含较大曲率的模型进行模拟计算时，例如进行全球模型中地震波传播模拟，模型的曲率对地震波传播的影响是不可忽视的，此时使用极坐标系进行模型的处理是非常直观合理的。另外在模拟一些实际模型中将不可避免的引入地形的问题，此时如果忽视地形的影响将会对结果的准确性有较大的影响。多种方法和思路被应用于进行复杂地形的处理，包括使用阶梯网格近似地表形状，或者使用垂向变换网格来处理缓慢起伏地形的情况等。其中精度最高的方法是牵引力镜像法，可以保持空间四阶精度。本文将贴体网格下的牵引力镜像法引入极坐标系，从而希望在考虑地形的情况下可以准确高效的进行全球模型地震波传播的模拟。

在全球模型的数值模拟方面，已经有一些算法被开发出来，例如伪谱法和谱元法，而目前并没有成熟的适用于全球模型的有限差分算法。为了解决极坐标系网格在圆心处存在奇点，以及网格间距过小限制时间步长的问题，本文引入了嵌套网格的方法，通过在圆心处生成笛卡尔网格来避免以上的问题。两套网格分别进行离散，都可以保持优良的正交性。同时，在信息传递的过程中，不需要两套网格点和点空间位置的一一匹配，两套网格之间的信息传递通过高阶拉格朗日插值实现，从而保证了整体算法的精度。另外，极坐标系网格和笛卡尔网格的空间步长基本相同，保证了算法的稳定性和计算效率的充分利用。

本文的极坐标系嵌套网格有限差分算法，首先利用链式求导法则，将经典极坐标系下的弹性力学方程组推导至贴体网格下的速度应力方程组。然后利用自由地表边界条件推导出地表处速度垂向方向导数，并由牵引力镜像反对称条件得出地表处应力垂向方向导数的表达式。同时在自由地表下临界格点采用紧致差分格式以保证空间四阶精度。圆心处的均匀笛卡尔网格对速度应力方程组进行离散，两套网格的交换层通过拉格朗日插值进行信息的传递和更新。另外，本文开发了嵌套网格MPI并行计算模式，并讨论了不同并行模式的加速比和效率。

本文通过与CGFDM和间断伽辽金方法的波形对比验证了本算法在水平地形，起伏地形和复杂介质中的准确性，并模拟和讨论了全火星模型中地形对地震波传播的影响。结果表明，本文的极坐标系嵌套网格有限差分算法是一种可以用于计算包含较大曲率模型的准确高效的有限差分算法。

[1]Geli L, Bard P-Y,Jullien B.The effect of topography on earthquake ground motion:A review and new results[J]. Bulletin of the Seismological Society of America, 1988,78(1):42-63.
[2]Campillo M, Bouchon M. Synthetic SH seismograms in a laterally varying medium by the discrete wavenumber method[J]. Geophysical Journal International,1985, 83(1):307-317.
[3]Chen X. Seismogram Synthesis for Multi-Layered Media with Irregular Interfaces by Global Generalized Reflection/Transmission Matrices Method.I. Theory of Two-Dimensional SH Case[J]. Bulletin of the Seismological Society of America, 1988, 80(6A):1696-1724.
[4]Chen X. Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method.I. Applications for 2D SH case[J]. Bulletin of the Seismological Society of America, 1995, 85(4):1094-1106.
[5]Chen X. Seismogram Synthesis for Multi-Layered Media with Irregular interfaces by Global Generalized Reflection/Transmission Matrices Method. . Theory of 2D P-SV Case[J]. Bulletin of the Seismological Society of America, 1996, 86(2):389-405.
[6]臧楠.复杂地形地震波模拟的对接网格和嵌套网格有限差分算法研究[D].安徽.中国科学技术大学, 2022.
[7]AltermanZ, KaralF C. Propagation of elastic waves in layered media by finite difference methods[J]. Bulletin of the Seismological Society of America, 1968, 58(1):367-398.
[8]KellyKR,WardR W,Treitel S,et al. Synthetic seismograms; a finite- difference approach[J]. Geophysics, 1976, 41(1):2-27.
[9]Zahradnik J. Simple elastic finite-difference scheme[J]. Bulletin of the Seismological Society of America, 1995, 85(6):1879-1887.
[10]Oprsal I,Zahradnik J. Elastic finite-difference method for irregular grids[J].Geophysics, 1999, 64(1):240-250.
[11]Moczo, P., Kristek, J., Vavrycuk,et al. 3D heterogeneous staggered-grid Finite-Difference Modeling of Seismic Motion with Volume Harmonic and Arithmetic Averaging of Elastic Moduli and Densities[J]. Bulletin of the Seismological Society of America, 2002, 92(8):3042-3066.
[12]Tessmer E,Kosloff D. 3-D elastic modeling with surface topography by a Chebychev spectral method[J]. Geophysics, 1994, 59(3):464-473.
[13]Lysmer J,DrakeL A. A finite element method for seismology[J] .Methods of computational physics, 1972, (11):181-216.
[14]SerónFJ,SanzFJ,Kindelán M,et al. Finite-element method for elastic wave propagation[J]. Communications in Applied Numerical Methods, 1990, 6(5): 359-368.
[15]Komatitsch D, TrompJ. Introduction to the spectral element method for three-dimensional seismic wave propagation[J].Geophysical Journal International, 1999, 139(3):806-822.
[16]Komatitsch D,Tromp J. Spectral-element simulations of global seismic wave propagation—I.Validation[J].Geophysical Journal International, 2002, 149(2):390-412.
[17]Komatitsch D,TrompJ.Spectral-element simulations of global seismic wave propagation—II. Three-dimensional models, oceans, rotation and self-gravitation[J].Geophysical Journal International, 2002, 150(1):303-318.
[18]Komatitsch D,Vilotte J-P. The spectral element method:An efficient tool to simulate the seismic response of 2D and 3D geological structures[J].Bulletin of the Seismological Society of America, 1998, 88(2):368-392.
[19]Kosloff D,Baysal E. Forward modeling by a Fourier method[J]. Geophysics, 1982, 47(10):1402-1412.
[20]Chaljub, Emmanuel, Yann Capdeville and Jean-Pierre Vilotte. “Solving elastodynamics in a fluid-solid heterogeneous sphere: a parallel spectral element approximation on non-conforming grids.” Journal of Computational Physics, 2003, 187: 457-491.
[21]Dimitri Komatitsch, Jeroen Tromp, Introduction to the spectral element method for three-dimensional seismic wave propagation, Geophysical Journal International,1999, 139: 806-822.
[22]Yanbin Wang, Hiroshi Takenaka, Takashi Furumura, Modelling seismic wave propagation in a two dimensional cylindrical whole earth model using the pseudospectral method, Geophysical Journal International, 2001, 145: 689–708.
[23]Steger JL,Dougherty FC,Benek JA. A chimera grid scheme[J]. 1983. DOI:http://dx. doi.org/.
[24]Xavier Juvigny, Elodie Canonne and ChristopheBenoit. "Multigrid Algorithms for the Chimera Method," AIAA 2004-758.
[25]Shih TI-P.Overset grids:fundamentals and practical issues.AIAA 2002-3259,2002.
[26]阎超.计算流体力学方法及应用.北京：北京航空航天大学出版社，2006.
[27]Benek JA,Steger JL,Dougherty F.A flexible grid em- bedding technique with application to the Euler equations. AIAA Paper 83-1944,1983.
[28]Benek JA,Steger JL.Chimera:A grid- embedding technique.AEDC-TR-85-64,1985.
[29]LaBoozetta WF,Gatzke TD.MACGS-towards the complete grid generation system. AIAA Paper 94-1923,1994.
[30]Chiu IT, Meakin RL.On automating domain connectivity or overset grids.AIAA Paper 95-0854,1995.
[31]刘周，周伟江，杨云军.使用非结构重叠网格方法模拟多体分离[C]//第十五届全国计算流体力学会议.2012.
[32]ZhangW,ChenX.Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[J].Geophysical Journal International, 2006, 167(1):337-353.
[33]Appelo D,Petersson N A. A Stable Finite Difference Method for the Elastic Wave Equation on Complex Geometries with Free Surfaces[J].Communications in Computational Physics, 2009, 5(1):84-107.
[34]Robertsson J O A. A numerical free-surface condition for elastic/visco elastic finite-difference modeling in the presence of topography[J]. Geophysics, 1996, 61(6):1921-1934.
[35]Lee S-J,Chen H-W,Liu Q,et al. Three-Dimensional Simulations of Seismic-Wave Propagation in the Taipei Basin with Realistic Topography Based upon the Spectral-Element Method[J]. Bulletin of the Seismological Society of America, 2008, 98(1):253-264.
[36]Zhang W,Shen Y,Chen X. Numerical simulation of strong ground motion for the M s8.0 Wenchuan earthquake of 12 May 2008[J].Science in China Series D:Earth Sciences, 2008, 51(12):1673-1682.
[37]Kaser M,Igel H,Sambridge M,et al. A Comparative Study of Explicit Differential Operators on Arbitrary Grids[J]. Journal of Computational Acoustics, 2001, 09(03):1111-1125.
[38]Nilsson S,Petersson N,Sjogreen B,et al. Stable Difference Approximations for the Elastic WaveEquationin SecondOrder Formulation[J].SIAM Jounal on Numerical Analysis, 2007, 45(5):1902-1936.
[39]Zahradník, J., Hron, F. Robust finite-difference scheme for elastic waves on coarse grids. Stud Geophys Geod ,1992, 36, 1–19.
[40]Takao Ohminato, Bernard A. Chouet. A free-surface boundary condition for including 3D topography in the finite-difference method. Bulletin of the Seismological Society of America, 1997, 87 (2): 494–515.
[41]Ekkehart Tessmer, Dan Kosloff, Alfred Behle. Elastic wave propagation simulation in the presence of surface topography, Geophysical Journal International, 1992, 108: 621–632.
[42]Eric T. Chung, Chi Yeung Lam, Jianliang Qian. A staggered discontinuous Galerkin method for the simulation of seismic waves with surface topography.GEOPHYSICS, 2015, 80:4, T119-T135.
[43]徐剑侠,张伟,陈晓非.极坐标系有限差分中起伏地表边界条件处理[J].地球物理学报, 2021.DOI:10.6038/cig2021M0316.
[44]Zhang W,Chen X. Traction image method for irregular free surface boundaries in finite difference seismic wave simulation[J]. Geophysical Journal of the Royal Astronomical Society, 2010, 167(1):337-353.
[45]Cerjan C,Kosloff D,Kosloff R,et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equations[J].Geophysics, 1985, 50(4):705-708.
[46]Wang ZJ.A fuly conservative structured/unstructured chimera grid scheme AIAA-95-0671,1995.
[47]Chicheportich J,Gloerfelt X.Study of interpolation meth- ods for high-accuracy computations on overlapping grids. Computers & Fluids,2012,68:112-133.
[48]常兴华,马戎,王年华,等.非结构重叠网格并行化隐式装配技术研究[C]//第四届全国非定常空气动力学学术会议.2018.
[49]刘鑫,陆林生.重叠区域找重策略和插值方法的研究[J].计算机应用研究, 2006, 23(7):3.DOI:10.3969/j.issn.1001-3695.2006.07.007.
[50]吴杰,赵慧勇,贺元元.重叠动网格算法及其实现[C]//全国激波与激波管学术交流会.2012.
[51]田书玲.基于非结构网格方法的重叠网格算法研究[D].南京航空航天大学
[2024-04-22].DOI:10.7666/d.d052041.
[52]Marstin CW,McConnaughey HV.Computational problems on composite grids.AIAA ,1984,84-1611.
[53]Farrell, P., Piggott, M.D., Pain, C.C., Gorman, G. Conservative interpolation between unstructured meshes via supermesh construction. Computer Methods in Applied Mechanics and Engineering, 2009, 2632-2642.
[54]Moon YJ, Liou MS. Conservative treatment of boundary interfaces for overlaid grids and multi - level grid adaptations.AIAA, 1989, Paper 89-1980.
[55]Zhao Xiang, Guan Huanwen, Yang Zhi, et al. An implicit and globally conservative unstructured chimera grid method. AIAA, 2011, 2011 - 777.
[56]Chesshire G, Henshaw W. Composite overlapping meshes for the solution of partial diferential equations. J Comput Phys, 1990, 90: 1- 64.
[57]Zhang Xing,NiSaizhen,He Guowei.A pressure-correction method and ts applications on an unstructured chimera grid. Computers& Fluids, 2008,37:993-1010.
[58]Zhang Xing. Computation of viscous incompressible flow using pressure correction method on unstructured chimera grid. International Journal of Computational Fuid Dynamics, 2006, 20(9): 637-650.
[59]Anderson D L, Miller W F, Latham G V, et al. Seismology on Mars. Journal of Geophysical Research, 1977, 82(28): 4524‒4546.
[60]Clinton J F, Giardini D, Böse M, et al. The marsquake service — building a Martian seismicity catalogue for InSight. Space Sci Rev, 2018, 214: 1‒33.
[61]Knapmeyer-Endrun B, Ceylan S, van Driel M. Crustal S-wave velocity from apparent incidence angles: a case study in preparation for InSight. Space Sci Rev, 2018, 214: 1‒40.
[62]Zheng Y, Nimmo F, Lay T. Seismological implications of a lithospheric low seismic velocity zone in Mars. Phys Earth Planet Inter, 2015, 240: 132–141.
[63]Khan A, van Driel M, Böse M, et al. Single-station and single-event marsquake location and inversion for structure using synthetic Martian waveforms. Phys Earth Planet Inter, 2016, 258: 28–42.
[64]Lognonné P, Karakostas F, Rolland L, et al. Modeling of atmospheric-coupled Rayleigh waves on planets with atmosphere: from Earth observation to Mars and Venus perspectives. J Acoust Soc Am, 2016, 140(2): 1447–1468.
[65]Bissig F, Khan A, van Driel M, et al. On the detectability and use of normal modes for determining interior structure of Mars. Space Sci Rev, 2018, 214: 1‒28.
[66]Bozdag E, Ruan Y, Metthez N, et al. Simulations of seismic wave propagation on Mars. Space Sci Rev, 2017, 211(1/2/3/4): 571–594.
[67]Sohl F, Spohn T. The interior structure of Mars: implications from SNC meteorites. Journal of Geophysical Research, 1997, 102: 1613‒1635.
[68]Stahler S C, Khan A, Banerdt W B, et al. Seismic detection of the martian core[J]. Science, 2021, 373(6553): 443-448.