三角网格有限元法波动模拟的数值频散及稳定性研究

三角网格有限元法能够准确模拟复杂构造和复杂介质条件下的地震波场,数值频散和稳定性条件是地震波数值模拟中参数选择的主要依据.基于均匀的线性三角网格单元,根据结构刚度矩阵的组装原理以及平面波理论,推导了集中质量矩阵下两种网格结构的声波频散函数以及稳定性条件,并对数值频散特性以及稳定性进行了详细研究:三角网格单元中波动的数值频散除了受到空间采样间隔、单元网格纵横比和波传播方向等常规因素的影响外,还受到网格布局的影响,过锐或过钝的三角单元会对波动数值频散产生不良的影响,不同类型的单元网格、单元纵横比对应着不同的稳定性条件,正三角单元中的波动具有较好的数值频散特性,其数值各向异性(频散随波传播方向的变化)效应最弱,稳定性条件也较为宽松.最后通过数值模拟直观地验证了以上分析结果,为有限元正演三角网格的剖分和参数的设置提供一定的理论依据.

The study for numerical dispersion and stability of wave motion with triangle-based finite element algorithm

The forward modeling of wave equation plays a very important role in seismic data acquisition, processing and interpretation. In order to accurately simulate the propagation of seismic wave in underground medium, not only requiring geophysical model is consistent with the actual formation, but also needing high precision numerical simulation method. Finite element method (FEM) with triangular grid can divide arbitrary complex boundary effectively, so it can accurately simulate the propagation of seismic wave field in complex medium. The spatial discretization of continuous medium by FEM introduces dispersion errors to the solution of wave equation, which causes phase velocity (numerical phase velocity) varies with the frequency of seismic wave that do not actually exist in the continuum, and this artefact is called numerical dispersion. In general, numerical dispersion not only has no practical significance, but also affects the understanding of the real fluctuation phenomenon. Therefore, it is very necessary to clarify the influence factors of the numerical dispersion, which can help to improve the accuracy of the numerical simulation.This paper focuses on studying the numerical dispersion and stability of wave motion with triangle-based finite element algorithm. Considering two kinds of commonly used grid structures (denoted by I mesh structure and II mesh structure, respectively), according to the assemble principle of structural stiffness matrix and plane wave theory, we obtain dispersion functions under the assumption that computational area is uniform and without borders. The stability conditions are obtained in the condition of second order middle difference, to provide reference for the selection of time step. In order to fully understand the characteristics of numerical dispersion in triangular meshes, we analysed quantitatively the effect of spatial sampling interval, the propagation direction of seismic wave, the ratio of vertical to horizontal (mesh shape) on numerical dispersion. To obtain optimal mesh generation, the influences of mesh shape are studied mainly. In a practical simulation, we should consider the numerical calculation accuracy. Therefore, the maximum dispersion errors of the two kinds of grid structures are analyzed comparatively. Finally, we verify our conclusions by wave-field simulations. The authors hope that the study in this paper can provide some useful suggestions for the division of triangular mesh and parameter setting.From the theoretical analyses and numerical experiments, the following results can be gained: (1) The numerical phase velocity varies significantly with the propagation direction, and this variation is periodic, the period of the numerical dispersion is π except regular triangle mesh, the period of which is π/3, because the mesh structure is invariant under π/3 rotation. (2) When the ratio of vertical to horizontal is equal to 1, the maximum dispersion error and minimum dispersion error of I mesh structure appear in the direction of θ=0° and θ=45°, respectively; and the maximum dispersion error and minimum dispersion error of II mesh structure appear in the direction of θ=90°and θ=0°, respectively. (3) In terms of I mesh structure, the dispersion error in horizontal direction has nothing to do with the ratio of vertical to horizontal, which is also the maximum dispersion error, in other words, by reducing the ratio of vertical to horizontal to suppress the numerical dispersion is not an effective method. (4) In terms of II mesh structure, the maximum dispersion error can be suppressed effectively by reducing the ratio of vertical to horizontal as long as it does not appear obtuse triangle mesh; it is worth pointing out that regular triangle mesh almost removes the directional dependence of the phase velocity, which means the numerical anisotropy is the weakest. (5) In the case of no visible numerical dispersion, the number of sampling point in the wavelength corresponding to peak frequency of source wavelet of II mesh structure is less than that of I mesh structure, and this can help improve the calculation efficiency and reduce the occupation of memory. (6) The stability factor of II mesh structure is also larger than that of I mesh structure.The dispersion characteristics analysed in the above can provide some theoretical guidances for mesh generation. The arrangements of elements should be chosen with care; unreasonable meshes will reduce the precision of numerical simulation, and even lead to erroneous results. The best arrangement of triangles appears to be hexagonal, which has the further benefit that it almost removes the directional dependence of the phase velocity, which can be useful in maintaining wave front definition.