**1. Introduction**

Elastic wave numerical simulation is an important means by which to study the law of seismic wave propagation within the ground, and it plays a vital role in the whole process of seismic exploration, including data acquisition, processing, and interpretation. Many numerical simulation methods, including the finite element [1,2], finite difference [3–5], spectral element [6], and pseudo spectral [7] methods have been successfully applied to elastic wave forward modeling. The staggered grid finite difference method (SGFD) is widely used in the forward modeling of elastic wave equations owing to its high calculation efficiency, high accuracy, and convenient implementation process. However, there are two problems with SGFD. One is that the fixed grid step may discretize the interface to the wrong depth, resulting in an inaccurate travel time for the reflected waves. The second is that a stepped interface will appear when the undulating interface is discretized with a regular grid, which will generate false diffracted waves. To eliminate grid diffraction, many scholars have carried out extensive research, including using variable grid algorithms [8], vertical grid mapping based on coordinate transformation [9], body fitted grids [10], and other methods. These methods effectively suppress stepped diffraction to a certain extent, but they also have problems. For example, the variable mesh algorithm essentially uses a smaller mesh step size to discretize the undulating interface, so it cannot fundamentally

**Citation:** Liu, S.; Zhou, Z.; Zeng, W. Simulation of Elastic Wave Propagation Based on Meshless Generalized Finite Difference Method with Uniform Random Nodes and Damping Boundary Condition. *Appl. Sci.* **2023**, *13*, 1312. https://doi.org/ 10.3390/app13031312

Academic Editor: Emanuel Guariglia

Received: 8 December 2022 Revised: 10 January 2023 Accepted: 13 January 2023 Published: 18 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

eliminate the stepped diffraction; the difference method based on coordinate transformation requires a corresponding transformation of the wave equation; and the implementation process is relatively complex.

In recent years, the meshless method, which avoids mesh generation, discretizes the solution area into a series of independent nodes and constructs approximate functions on these discrete nodes to obtain linear equations [11]. Since the nodes change flexibly with the velocity model, the meshless method can effectively avoid grid diffraction [12]. The generalized finite difference method is a meshless method with a simple principle, flexible node arrangement, and high calculation accuracy, and it has been widely used to solve a variety of mathematical and engineering problems [13–15]. Jensen (1972) first used the Taylor series expansion of several adjacent nodes around a center point, based on a distance function, to solve differential equations [16]. Later, Benito et al. (2001) developed the explicit formula of the generalized finite difference method [17]. Based on the Taylor series expansion and weighted least square fitting, the partial derivative of the unknown quantity in the differential equation is expressed as a linear combination of the function values of several adjacent nodes. Thus, the basic theory of the generalized finite difference method was formed. Since then, many scholars have improved the theory of this method and its implementation technology. For example, Benito and Ureña systematically analyzed various factors that affect the accuracy of the GFDM calculations, such as node generation [18,19], star of nodes shape [20], and weight function [17], and they pointed out that as the discretization (cloud of nodes) becomes more regular, the results obtained become more stable [21].

Ureña et al. (2011) first applied the generalized finite difference method to the forward modeling of elastic wave equations, and proposed an adaptive method to minimize the effect of the irregularity of node distributions on dispersion [22]; on this basis, Ureña et al. (2012) discussed the dispersion and stability of elastic wave forward modeling using the generalized finite difference method under regular and irregular node conditions [23]. Benito et al. (2013) further studied the generalized finite difference method for solving the problem of determining seismic wave propagation in homogeneous isotropic media [24]. Benito et al. (2015) discussed the influence of node settings on simulation accuracy using circular hole models with Dirichlet uniform boundary conditions and square hole model scatterers with free boundary conditions [25]. Benito et al. (2017) applied the generalized finite difference method (GFDM) to solve the problem of elastic wave propagation, and analyzed the influence of the type of node clouds (regular and irregular) on discretization [17]. Salete et al. (2017) put forward a generalized finite difference scheme to solve the two-dimensional seismic wave propagation problem with a perfectly matched layer absorption boundary, and discussed the stability of PML. They also pointed out that the stability condition at the boundary of PML is stricter than that in the internal computational domain [26].

In this study, the elastic wave equation was solved using the meshless generalized finite difference method. This method first discretizes the solution area into a series of independent nodes, and then constructs an approximate function for these discrete nodes based on the Taylor series expansion and weighted least squares fitting. Finally, it solves the linear system to obtain an approximate solution for the boundary value of the elastic wave equation. An algorithm combining the Poisson disk node generation algorithm and the centroid Voronoi node adjustment algorithm is suggested to improve the stability of the node "star". The Cerjan damping boundary condition is introduced to avoid the instability caused by the absorbing boundary conditions. In some worked examples, the GFDM is compared with SGFD, and analytical solutions based on homogeneous, horizontal layered, undulating interface, and fault models are also tested.

### **2. Methodology**

### *2.1. Elastic Wave Equation*

The two-dimensional elastic wave equation can be expressed as:

$$\begin{cases} \rho \frac{\partial^2 u\_x}{\partial t^2} = (\lambda + 2\mu) \frac{\partial^2 u\_x}{\partial x^2} + (\lambda + \mu) \frac{\partial^2 u\_z}{\partial x \partial z} + \mu \frac{\partial^2 u\_x}{\partial z^2} + f\_x\\ \rho \frac{\partial^2 u\_z}{\partial t^2} = \mu \frac{\partial^2 u\_z}{\partial x^2} + (\lambda + \mu) \frac{\partial^2 u\_x}{\partial x \partial z} + (\lambda + 2\mu) \frac{\partial^2 u\_z}{\partial z^2} + f\_z \end{cases} \tag{1}$$

where *t* is the time, *ρ* is the density of the medium, *λ* and *μ* are Lame constants, *λ* = *ρ v*2 *<sup>p</sup>* − 2*v*<sup>2</sup> *s* , *μ* = *ρv*<sup>2</sup> *<sup>s</sup>*, *vp*, and *vs* are the velocities of compressional and shear waves, *ux* and *uz* represent displacement in the *x* and *z* directions, and *fx* and *fz* represent source terms in the *x* and *z* directions.

### *2.2. Central Difference for Time Partial Derivative Approximation*

The displacement at time *t* + Δ*t* and *t* − Δ*t* can be expanded at time *t* via the Taylor series:

$$u(t + \Delta t) = u(t) + \Delta t \frac{\partial u(t)}{\partial t} + \frac{\left(\Delta t\right)^2}{2!} \frac{\partial^2 u(t)}{\partial t^2} + \frac{\left(\Delta t\right)^3}{3!} \frac{\partial^3 u(t)}{\partial t^3} + \dots \tag{2}$$

$$u(t - \Delta t) = u(t) - \Delta t \frac{\partial u(t)}{\partial t} + \frac{\left(\Delta t\right)^2}{2!} \frac{\partial^2 u(t)}{\partial t^2} - \frac{\left(\Delta t\right)^3}{3!} \frac{\partial^3 u(t)}{\partial t^3} + \dots \tag{3}$$

Truncating the 2M-order obtains the central difference expression:

$$u(t + \Delta t) = 2u(t) - u(t - \Delta t) + 2\sum\_{m=1}^{M} \frac{\Delta t^{2m}}{2m!} \frac{\partial^{2m} u(t)}{\partial t^{2m}} + O\left(\Delta t^M\right) \tag{4}$$

When *M* = 1, then:

$$u(t + \Delta t) = 2u(t) - u(t - \Delta t) + \Delta t^2 \frac{\partial^2 u(t)}{\partial t^2} + O\left(\Delta t^2\right) \tag{5}$$

The time partial derivatives *∂*2*ux*(*t*)/*∂t* <sup>2</sup> and *∂*2*uz*(*t*)/*∂t* <sup>2</sup> in Equation (1) can be converted using Equation (5) into the following form:

$$\begin{cases} \begin{aligned} \boldsymbol{u}\_{x}(t+\Delta t) &= 2\boldsymbol{u}\_{x}(t) - \boldsymbol{u}\_{x}(t-\Delta t) \\ &+ \Delta t^{2} \left[ \frac{\lambda+2\mu}{\rho} \frac{\partial^{2}\boldsymbol{u}\_{x}}{\partial\mathbf{x}^{2}} + \frac{\mu}{\rho} \frac{\partial^{2}\boldsymbol{u}\_{x}}{\partial\boldsymbol{z}^{2}} + \frac{\lambda+\mu}{\rho} \frac{\partial^{2}\boldsymbol{u}\_{x}}{\partial\boldsymbol{z}\partial\boldsymbol{z}} \right] + f\_{x}(t+\Delta t) \\ \boldsymbol{u}\_{z}(t+\Delta t) &= 2\boldsymbol{u}\_{z}(t) - \boldsymbol{u}\_{z}(t-\Delta t) \\ &+ \Delta t^{2} \left[ \frac{\lambda+2\mu}{\rho} \frac{\partial^{2}\boldsymbol{u}\_{z}}{\partial\boldsymbol{x}^{2}} + \frac{\mu}{\rho} \frac{\partial^{2}\boldsymbol{u}\_{z}}{\partial\boldsymbol{z}^{2}} + \frac{\lambda+\mu}{\rho} \frac{\partial^{2}\boldsymbol{u}\_{x}}{\partial\boldsymbol{x}\partial\boldsymbol{z}} \right] + f\_{z}(t+\Delta t) \end{aligned} \tag{6}$$

### *2.3. GFDM for Spatial Partial Derivative Approximation*

In the GFDM, through utilizing the Taylor series expansion and weighted least squares fitting, the derivatives of unknown variables can be approximated by a linear combination of function values of some neighboring nodes, which are located inside a star. First, a regular or irregular cloud of points is generated in the computation domain and along the boundary. For a given node *i*, which is denoted as a central node, the *m* nearest nodes surrounding the node *i* will be found. The concept of the star refers to a group of established nodes in relation to the central node, as shown in the black circle in Figure 1. Each node is assigned an associated star.

**Figure 1.** An irregular cloud of points and the selection of star via distance criterion.

If *u*<sup>0</sup> is the function value at the central node *x*<sup>0</sup> of the star and *ui*(*i* = 1, 2, ··· , *m*) is the function value at one of the rest of the nodes, then the Taylor series expansion around the central node can be expressed in the following form:

$$u\_i = u\_0 + h\_i \frac{\partial u\_0}{\partial \mathbf{x}} + l\_i \frac{\partial u\_0}{\partial z} + \frac{h\_i^2}{2} \frac{\partial^2 u\_0}{\partial \mathbf{x}^2} + \frac{l\_i^2}{2} \frac{\partial^2 u\_0}{\partial z^2} + h\_i l\_i \frac{\partial^2 u\_0}{\partial \mathbf{x} \partial z} + \cdots \tag{7}$$

where the coefficients *hi* = *xi* − *x*0, *li* = *zi* − *z*0, and (*x*0, *z*0) are the coordinates of the central node, and (*xi*, *zi*) are the coordinates of the node *i* in the star. By truncating the Taylor series with the second derivative, the residual function *B*(*u*) can be defined by the following equation:

$$B(\boldsymbol{u})\_{2^{\text{nd}}} = \sum\_{i=1}^{m} \left[ \left( u\_0 - u\_i + h\_i \frac{\partial u\_0}{\partial \mathbf{x}} + l\_i \frac{\partial u\_0}{\partial \boldsymbol{z}} + \frac{h\_i^2}{2} \frac{\partial^2 u\_0}{\partial \mathbf{x}^2} + \frac{l\_i^2}{2} \frac{\partial^2 u\_0}{\partial \boldsymbol{z}^2} + h\_i l\_i \frac{\partial^2 u\_0}{\partial \mathbf{x} \partial \mathbf{z}} \right) \mathbf{u}\_i \right]^2 \tag{8}$$

By truncating the Taylor series with the fourth derivative, the residual function *B*(*u*) can be defined by the following:

*<sup>B</sup>*(*u*)4*th* <sup>=</sup> *<sup>m</sup>* ∑ *i*=1 *u*<sup>0</sup> <sup>−</sup> *ui* <sup>+</sup> *hi ∂u*<sup>0</sup> *<sup>∂</sup><sup>x</sup>* + *li ∂u*<sup>0</sup> *<sup>∂</sup><sup>z</sup>* <sup>+</sup> *<sup>h</sup>*<sup>2</sup> *i* 2 *∂*2*u*<sup>0</sup> *<sup>∂</sup>x*<sup>2</sup> <sup>+</sup> *<sup>l</sup>* 2 *i* 2 *∂*2*u*<sup>0</sup> *<sup>∂</sup>z*<sup>2</sup> + *hili ∂*2*u*<sup>0</sup> *∂x∂z* +*h*<sup>3</sup> *i* 6 *∂*3*u*<sup>0</sup> *<sup>∂</sup>x*<sup>3</sup> <sup>+</sup> *<sup>l</sup>* 3 *i* 6 *∂*3*u*<sup>0</sup> *<sup>∂</sup>z*<sup>3</sup> <sup>+</sup> *<sup>h</sup>*<sup>2</sup> *i li* 2 *∂*3*u*<sup>0</sup> *<sup>∂</sup>x*2*∂<sup>z</sup>* <sup>+</sup> *hil* 2 *i* 2 *∂*3*u*<sup>0</sup> *<sup>∂</sup>x∂z*<sup>2</sup> <sup>+</sup> *<sup>h</sup>*<sup>4</sup> *i* 24 *∂*4*u*<sup>0</sup> *<sup>∂</sup>x*<sup>4</sup> <sup>+</sup> *<sup>l</sup>* 4 *i* 24 *∂*4*u*<sup>0</sup> *∂z*<sup>4</sup> +*h*<sup>3</sup> *i li* 6 *∂*4*u*<sup>0</sup> *<sup>∂</sup>x*3*∂<sup>z</sup>* <sup>+</sup> *hil* 3 *i* 6 *∂*4*u*<sup>0</sup> *<sup>∂</sup>x∂z*<sup>3</sup> <sup>+</sup> *<sup>h</sup>*<sup>2</sup> *i l* 2 *i* 4 *∂*4*u*<sup>0</sup> *∂x*2*∂z*<sup>2</sup> *wi* 2 (9)

In Equations (8) and (9), *wi* = *w*(*hi*, *li*) represents the distance weight function of the star. In all the examples considered in this paper, the weighting function was chosen as per [27]:

$$w\_i = \begin{cases} 1 - 6\left(\frac{d\_i}{d\_{\text{max}}}\right)^2 + 8\left(\frac{d\_i}{d\_{\text{max}}}\right)^3 - 3\left(\frac{d\_i}{d\_{\text{max}}}\right)^4, & d\_i \le d\_{\text{max}}\\ 0, & d\_i > d\_{\text{max}} \end{cases} \tag{10}$$

where *di* denotes the distance between nodes (*xi*, *zi*) and (*x*0, *z*0), and *d*max is the distance between the central node and the farthest node in the star. Let the following terms be defined:

$$\mathbf{D}\_{\mathbf{u}} = \left[ \frac{\partial u\_0}{\partial \mathbf{x}}, \frac{\partial u\_0}{\partial z}, \frac{\partial^2 u\_0}{\partial \mathbf{x}^2}, \frac{\partial^2 u\_0}{\partial z^2}, \frac{\partial^2 u\_0}{\partial \mathbf{x} \partial z'}, \cdots \right] \tag{11}$$

$$\mathbf{p}\_{i} = \left\{ h\_{i}, l\_{i}, \frac{h\_{i}^{2}}{2}, \frac{l\_{i}^{2}}{2}, h\_{i}l\_{i}, \dots \right\}, i = 1, 2, \dots, m \tag{12}$$

$$\mathbf{P} = \begin{bmatrix} \mathbf{P}\_1 \\ \mathbf{P}\_2 \\ \vdots \\ \mathbf{P}\_m \end{bmatrix} = \begin{bmatrix} h\_1 & l\_1 & \cdots & h\_1 l\_1 & \cdots \\ h\_2 & l\_2 & \cdots & h\_2 l\_2 & \cdots \\ \vdots & \vdots & \ddots & \vdots & \cdots \\ h\_m & l\_m & \ddots & h\_m l\_m & \cdots \end{bmatrix} \tag{13}$$

The residual function defined in Equations (8) and (9) can be expressed in matrix form:

$$B(\boldsymbol{\mu}) = \left(\mathbf{PD}\_{\mathbf{u}} + \mathbf{u}\_0 - \mathbf{u}\right)^T \mathbf{W} (\mathbf{PD}\_{\mathbf{u}} + \mathbf{u}\_0 - \mathbf{u})\tag{14}$$

where **u** = [*u*1, *u*2, ··· , *um*] *<sup>T</sup>*, **u**<sup>0</sup> = [*u*0, *u*0, ··· , *u*0] *<sup>T</sup>*, and **W** = *diag*- *w*2 <sup>1</sup>, *<sup>w</sup>*<sup>2</sup> <sup>2</sup>, ··· , *<sup>w</sup>*<sup>2</sup> *m* .

To minimize the residual function *B*(*u*) with respect to the unknown derivatives at the central point (*x*0, *z*0):

$$\frac{\partial B(\mu)}{\partial \mathbf{D}\_{\mathbf{u}}} = 0 \tag{15}$$

yields the following linear equation system,

$$\mathbf{AD}\_{\mathbf{u}} = \mathbf{b} \tag{16}$$

and

$$\mathbf{A} = \mathbf{P}^T \mathbf{W} \mathbf{P} \tag{17}$$

$$\mathbf{b} = \mathbf{P}^{\mathrm{T}} \mathbf{W}(\mathbf{u} \cdot \mathbf{u}\_{\mathbf{0}}) \tag{18}$$

According to Equations (16)–(18), the partial derivative vector **Du** can be expressed as:

$$\mathbf{D}\_{\mathbf{u}} = \mathbf{A}^{-1} \mathbf{b} = \mathbf{A}^{-1} \mathbf{P}^{T} \mathbf{W} (\mathbf{u} - \mathbf{u}\_{0}) = -\mathbf{a}u\_{0} + \sum\_{i=1}^{m} \mathbf{a}\_{i}u\_{i} \tag{19}$$

where **<sup>a</sup>** <sup>=</sup> *<sup>m</sup>* ∑ *i*=1 **a***i*, **a** = [*a*1, *a*2, *a*3, *a*4, *a*5] *<sup>T</sup>*, and *<sup>m</sup>* ∑ *i*=1 **a***<sup>i</sup>* = *m* ∑ *i*=1 *a*1*i*, *m* ∑ *i*=1 *a*2*i*, *m* ∑ *i*=1 *a*3*i*, *m* ∑ *i*=1 *a*4*i*, *m* ∑ *i*=1 *a*5*i T* . Expanding the partial derivative vector **Du**, we can obtain:

$$\begin{cases} \frac{\partial u\_0}{\partial x} = -a\_1 u\_0 + \sum\_{i=1}^m a\_{1i} u\_i\\ \frac{\partial u\_0}{\partial z} = -a\_2 u\_0 + \sum\_{i=1}^m a\_{2i} u\_i\\ \frac{\partial^2 u\_0}{\partial x^2} = -a\_3 u\_0 + \sum\_{i=1}^m a\_{3i} u\_i\\ \frac{\partial^2 u\_0}{\partial z^2} = -a\_4 u\_0 + \sum\_{i=1}^m a\_{4i} u\_i\\ \frac{\partial^2 u\_0}{\partial x \partial z} = -a\_5 u\_0 + \sum\_{i=1}^m a\_{5i} u\_i \end{cases} \tag{20}$$

Substituting Equation (20) into the Equation (6), the discrete formula for solving the elastic wave equation can be obtained:

$$\begin{cases} \begin{aligned} u\_{x,0}^{n+1} &= 2u\_{x,0}^{n} - u\_{x,0}^{n-1} + \Delta t^{2} \left[ \frac{\lambda + 2\mu}{\rho} \left( -a\_{3}u\_{x,0}^{n} + \sum\_{i=1}^{m} a\_{3i}u\_{x,i}^{n} \right) \right. \\ & \left. + \frac{\lambda + \mu}{\rho} \left( -a\_{5}u\_{z,0}^{n} + \sum\_{i=1}^{m} a\_{5i}u\_{z,i}^{n} \right) + \frac{\mu}{\rho} \left( -a\_{4}u\_{x,0}^{n} + \sum\_{i=1}^{m} a\_{4i}u\_{x,i}^{n} \right) \right] + f\_{x}^{n+1} \\\ u\_{z,0}^{n+1} &= 2u\_{z,0}^{n} - u\_{z,0}^{n-1} + \Delta t^{2} \left[ \frac{\lambda + 2\mu}{\rho} \left( -a\_{4}u\_{z,0}^{n} + \sum\_{i=1}^{m} a\_{4i}u\_{z,i}^{n} \right) \right. \\ & \left. + \frac{\lambda + \mu}{\rho} \left( -a\_{5}u\_{x,0}^{n} + \sum\_{i=1}^{m} a\_{5i}u\_{x,i}^{n} \right) + \frac{\mu}{\rho} \left( -a\_{3}u\_{z,0}^{n} + \sum\_{i=1}^{m} a\_{3i}u\_{z,i}^{n} \right) \right] + f\_{z}^{n+1} \end{aligned} \tag{21}$$

For the stability analysis, the aim is to construct a harmonic decomposition of the approximated solution at the grid points at a given time level (n). The amplification factor must be less than or equal to one in order to determine the stability limit. This has been studied by [23] and the condition for the stability of the star has been established as:

$$\Delta t \le \sqrt{\frac{4}{\left(v\_p^2 + v\_s^2\right)\left(|a\_3| + |a\_4|\right) + \left(v\_p^2 - v\_s^2\right)\sqrt{\left(a\_3 - a\_4\right)^2 + 4a\_5^2}}}\tag{22}$$

where *vp* and *vs* are the velocity of compressional and shear waves, respectively.
