Introduction

For the numerical simulations of near-field wave motions and the response of geological structures, the control equations of different media should be determined to obtain the reliable wave propagation characteristics [1–5]. Moreover, we need to truncate models of media in a finite-computational domain by introducing artificial boundaries. Inappropriately set artificial boundary conditions might incur spurious reflections, which not only affect the computational precision at inner grid nodes and boundary nodes and the resolution of wave-field simulation, but also interfere with the response of geological structure [6–8]. Numerous studies have been conducted on artificial boundary conditions since the late 1960s [6, 9–19]. Specifically, the Multi-Transmitting Formula (MTF), which is a local artificial boundary condition proposed by Liao, et al. [20, 21], has been favored because of its simple physical concept, wide adaptability, as well as easy implementation of decoupled high-precision numerical simulations of wave motions and for the wave scattering problem, its unique advantages are more obvious [20, 22–26].

Similar to other local artificial boundary conditions, the MTF is subjected to numerical instability problems in implementation of MTF into the numerical simulation by time-step integration, such as the high-frequency oscillation instability and drift instability. In regards to the high-frequency instability, its mechanism has been clarified, and measures that attempt to eliminate this instability have been proposed [27–30]. For drift instability, the mechanism of high-order drift instability has been explored by using numerical experiments, and a measure has been suggested to suppress it [31]. By reducing the order of MTF to order 1, the drift instability can be suppressed to a certain extent. However, the precision of MTF-1 is not enough to meet the needs of engineering application. Especially, a small computational domain which brings a relatively large incident angle between artificial boundary and input will amplify the drift instability. Therefore, it is worthwhile to further elaborate on the mechanism of drifting instability and to propose corresponding elimination measures. The numerical stability of local artificial boundary conditions is mathematically equivalent to that of the initial-boundary value problems of the hyperbolic partial differential equations. Moreover, for the latter the necessary and sufficient condition for one-dimensional numerical stability, known as the Gustafsson, Kreiss and Sundström (GKS) criterion [32], has been identified. To increase the applicability of this criterion to multi-dimensional case and simplify this complex mathematical theorem, the existing study interprets physical implication of this criterion, i.e., the internal traveling waves which satisfy both artificial boundary conditions and internal motion equations are not allowed [32–34].

In this study, based on above interpretation and combined with decoupled numerical simulations [21, 35], the mechanism of drift instability was clarified, and a simple measure for eliminating drift instability in the numerical simulation of wave motions was proposed, which is to introduce a modified operator $$ \gamma B_0^0$$ into the MTF. Moreover, the physical implication of $$ \gamma B_0^0$$ was explained by the spherical wave propagation principle and the medium damping effect; then, the effect of the $$ \gamma B_0^0$$ on numerical simulation accuracy was analyzed. Finally, the validity of the proposed stability measure was verified by numerical experiments of wave source problem and scattering problem.

Mechanism of and eliminating measure for MTF drift instability

The x-axis is set to be the outer normal of the artificial boundary, and the origin, o, coincides with a boundary node, as shown in Fig 1. The displacement of a one-way wave motion propagating forward along the x-axis, denoted as u(t,x), is a function of time t and coordinate x.

Fig 1

Schematic diagram of the artificial boundary node o and MTF computational points.

Based on the simulation of wave motion, the MTF can be expressed as [24]:

The backward moving operator $$ B_m^n$$ is defined as [24]:

The operator satisfies the following operation rules:

Using the backward moving operator, the MTF can be expressed as:

In the discrete model, the harmonic form of the interior node solution that satisfies the numerical integral stability condition of the nodal motion equation in the computational domain of x < 0 can be expressed as [33]:

The numerical stability of local artificial boundary conditions is mathematically equivalent to that of the initial-boundary value problems of the hyperbolic partial differential equations. The necessary and sufficient condition for one-dimensional numerical stability of the latter is known as the GKS criterion. According to the GKS criterion, in order to achieve the stability of MTF against any wave motions and prevent internal traveling waves from entering the computational domain, the MTF boundary conditions (Eq (4) should be not satisfied for the internal traveling wave solution (Eq (5)) [33, 34]. Substituting Eq (5) into Eq (4) and combining it with moving operator $$ B_m^n$$, the stability requirement of MTF for numerical integration can be expressed as:

Because Δt and Δx are positive real numbers, ω and k are real numbers with the same signs, the condition (Eq (6)) holds for all non-zero wavenumbers and non-zero frequencies, with the exception of zero frequency and zero wavenumber. In the case of zero frequency and zero wavenumber, the components of the numerical solution enter the computational domain through the boundary, which caused the violation of GKS criterion and leads to the drift instability. Based on this explanation, we proposed to introduce a small positive parameter into the boundary conditions in the numerical implementation of MTF, and Eq (4) should be rewritten as:

In the next section, the spherical wave propagation principle and the medium damping effect were utilized to analyze the physical implication of the modified operator, $$ \gamma B_0^0$$.

Physical implication of $$ \mathit{\gamma }{\mathit{B}}_{\mathbf{0}}^{\mathbf{0}}$$

The spherical wave generated by the point source S is assumed to propagate outward at a wave velocity c (as shown in Fig 2). The distance between the observation point A and point source S is r₀, and the distance between the observation point o and point source S is r₀+Δr. Assuming the displacement of point A is known to be u_A(t), the displacement of point o can be obtained according to the spherical wave propagation principle:

It is known through comparison with the MTF that according to the MTF with $$ \gamma B_0^0$$, the displacements of the inner nodes should be multiplied by the corresponding factor when extrapolating the displacement of the artificial boundary nodes with the displacements of the inner nodes.

Fig 2

Schematic diagram of the source and observation points.

For the sake of simplicity but without a loss of generality, let the MTF order N = 2 in Eq (9), and thus:

Eq (12) is another expression of the second-order MTF with $$ \gamma B_0^0$$, which is equivalent to “Eq (22)” on page 178 of Reference [25]. It can be seen that both the incident wave and the error wave propagate in the form of spherical waves with α = 1/1+γ. Therefore, it can be concluded that the MTF with $$ \gamma B_0^0$$ considers the geometric diffusion characteristics of the medium in the numerical simulation of wave motion.

In the following, the second-order MTF is still used as an example to provide another physical explanation of the MTF with $$ \gamma B_0^0$$. For the condition of N = 2, we can obtain [6]:

For the condition μ = ln(1+γ), we can obtain:

Through comparison of Eqs (16) and (10), it is easy to recognize that they are equivalent. Accordingly, it is explained that the MTF with $$ \gamma B_0^0$$ also introduces damping characteristics of the medium in the numerical simulation of wave motion.

Based on the two explanations of physical of the MTF with $$ \gamma B_0^0$$ described above, the MTF with $$ \gamma B_0^0$$ either consider the geometric diffusion characteristics of the medium or introduce the damping characteristics of the medium, thereby the two explanations of physical implications indicate that the MTF with $$ \gamma B_0^0$$ introduces the absorption mechanism of the medium to wave motion. Moreover, inspired by the introduction of the damping characteristics of the medium, it is worthwhile to explore how much damping that is proportional to the velocity of the particle motion should be introduced to eliminate the drift instability of the MTF, which is an issue that will be studied in the future.

Accuracy analysis of modified operator $$ \mathit{\gamma }{\mathit{B}}_{\mathbf{0}}^{\mathbf{0}}$$

The accuracy of the artificial boundary is usually expressed by the reflection coefficient. In the steady-state case, the reflection coefficient of the artificial boundary is generally defined as [24, 25]:

Next, the reflection coefficients of the MTF with $$ \gamma B_0^0$$ are discussed in two limit cases, Case A and Case B, in which the total wave field is composed only of incident waves in Case A and total wave field is composed of both incident waves and fully developed reflected waves in Case B.

Case A

For steady-state wave motions, Eq (9) can be expressed as:

where U_j(j = 0,1,2,3,…) denotes the amplitude of motion at the computational point j (Fig 1), a = exp(iωΔt), $$ i=\sqrt{‐1}$$, and ω is the circular frequency.

In the calculation of the reflection coefficient, the term related to time appears in expression of the total displacement and the incident displacement at the same time, which can be eliminated finally during the calculation progress. Considering that the time factor has no influence on the result of reflection coefficient in this case, the time factor (exp(iωpΔt)) in the discrete form of the incident plane harmonic wave is ignored (here and hereinafter), and the discrete form can be expressed as:

where $$ U_j^I$$ and $$ U_0^I$$ are amplitudes of the incident harmonic wave at the discrete nodes, x = jc_aΔt, and the artificial boundary node, x = 0, respectively. Moreover:

where c_x is the apparent wave velocity propagating along the x-axis, and a_x represents the phase change of displacement when the traveling distance of incident plane wave is c_aΔt on the x-axis. Assuming the wave velocity of the incident plane wave to be c and the incident angle to be θ, thus c_x = c/cosθ. Assuming c_a = c, we can obtain:

In this case, $$ U_j=U_j^I$$. According to Eq (20), we can have:

Substituting Eq (23) into Eq (19) results in:

where $$ b_I={\displaystyle \sum _{j=1}^N{\left(-1\right)}^{j+1}{C}_j^N{a}^{-j}{a}_x^j/{\left(1+\gamma \right)}^j}$$. Using the binomial Equation $$ {\left(1+x\right)}^N=1+{\displaystyle \sum _{j=1}^N{C}_j^N{x}^j}$$, b_I can be expressed as:

Substitution of Eq (24) into Eq (18) results in:

Substitution of a and a_x into Eq (26) leads to:

where period T = 2π/ω.

Case B

For Case B, the incidence of a plane scalar wave is discussed here. c_a is assumed to be equal to the scalar wave velocity c, and the incident angle is θ. Therefore, c_x = c/cosθ. Because the total displacement at the discrete node, x = −jc_aΔt, on the x-axis, is composed of the incident wave and the fully developed reflected wave, we can derive:

where $$ U_j^I$$ is the incident wave displacement, given by Eq (20), and $$ U_j^R$$ is the reflected wave displacement, which can be determined according to:

Substituting Eqs (28), (20) and (29) into Eq (19) results in:

where b_I is determined by Eq (25), and b_R is determined by:

Substituting $$ U_0=U_0^I+U_0^R$$ into Eq (30) results in:

According to the definition of the reflection coefficient, namely Eq (17), substituting b_I and b_R into Eq (31) results in:

Substituting a and a_x into Eq (32) results in:

To more intuitively explain the effect of $$ \gamma B_0^0$$ on the MTF reflection coefficient, we plot the relationship curves between the reflection coefficient and the incident angle corresponding to different γ, according to Eqs (27) and (33) for N = 2 and Δt/T = 1/10 and Δt/T = 1/20. As shown in Fig 3, when Δt/T was 1/10 or 1/20, the reflection coefficients with the modified operator γ of 0.02 became similar to the result that without γ, which demonstrated that the modified operator γ could have little influence on the computational precision if its value is small. However, with the increase of γ to 0.05, the error of the reflection coefficient after adding the γ became larger. Moreover, it was seen that when γ was 0.05, the overall error of reflection coefficients was larger with the decrease of Δt/T from 1/10 to 1/20. As for the case B shown in Fig 4, the effect of different value of modified operator on the reflection coefficient is consistent with that shown in case A (shown in Fig 3), which indicated that if the added modified operator is small enough, its effect on MTF can be almost ignored.

Fig 3

Reflection coefficient in case A: (a) Δt/T = 1/10; (b) Δt/T = 1/20.

Fig 4

Reflection coefficient in case B: (a) Δt/T = 1/10; (b) Δt/T = 1/20.

Numerical experiments

The decoupling method for the numerical simulation of wave motion [21, 35], which is a combination of the MTF and Lumped mass explicit finite element method, is used to verify the proposed measure of the MTF stably implementation by numerical experiments. Typical examples of the three-dimensional wave source problem and scattering problem are used to examine how this measure effectively eliminates drift instability and ensures the computational precision of numerical simulation.

Wave source problem

Wave motions generated by a concentrated source of force in a homogeneous, isotropic linear elastic infinite medium are considered. With a Cartesian coordinate system ox₁x₂x₃, the origin of the coordinate system is the same as the action point of the concentrated force. Assumed to act in the direction along the x₃ axis, the concentrated force, p(t), has a time function of an approximate δ-pulse [25] (Fig 5). The numerical simulation of wave motion is performed in the homogeneous box domain, and all of boundaries are the artificial boundary, the sizes of which are represented by B₁, B₂ and B₃ (Fig 6). The computational domain is discretized by the cubic finite elements with a spatial step size of Δx, whose value is determined by the computational precision of the finite element numerical simulation of the wave motion [25], and a system of ordinary differential equations relative to inner nodes can be formed on this basis. The numerical integration of ordinary differential equations is carried out by the central difference method, and the time step Δt of numerical integration is determined by the stability criterion of the central difference method. The discrete equation of nodes on the artificial boundary is a second-order MTF, in which the displacements of the computational points are obtained by the three-point interpolation method with second-order precision based on displacements of finite element nodes [25], the expression is as follows:

where

in which $$ {\overline{u}}_j^p$$ and $$ {\overline{u}}_j^{p-1}\left(j=0,1,2\right)$$ represent the displacement responses of the finite element nodes at the times of pΔt and (p−1)Δt, respectively.

Fig 5

Time function of the concentrated force.

Fig 6

Computational domain in the wave source problem (Quadrant I).

The calculation performed in the wave source problem with dimensionless parameters and numerical simulation results of wave motion corresponded to the case of B₁ = B₂ = B₃ = 0.4, medium mass density ρ = 1, longitudinal wave velocity $$ c_1=\sqrt{3}$$, shear wave velocity c₂ = 1, Δx = 0.02, Δt = 0.01, and artificial wave velocity c_a = 1. The displacements of inner observation point A (0.2, 0.2, 0.2) and boundary observation point B (0.4, 0.4, 0.4) under different values of γ were shown in Figs 7 and 8 and compared with the analytical solution [36] respectively.

Fig 7

Displacement time-histories of observation point A.

Fig 8

Displacement time-histories of observation point B.

As shown in the Fig 7, all three components of displacement responses of the observation point A which is inside the computational domain showed significant drift instability when the MTF is without the modified operator. After the addition of the modified operator, the drift instability was effectively suppressed, and as the value of γ increased from 0.005 to 0.02, the suppression effect seemed to become better. However, when the value of γ continues to increase to 0.05, the error by the addition of modified operator appeared. In the Fig 7C, it can be seen that when the value of γ was 0.05, the trough of the component of displacement response in the numerical simulation was obviously deviated from the exact solution, which means that the suppression effect of drift instability became poor at this time. As shown in the Fig 8, it was found that the influence of the value of γ on the displacement response of the boundary point B was similar to the influence on the inner point A. Besides, the addition γ had a worse suppression effect on the component of Z-direction of displacement response, which was due to the fact that the incident angle of the vertical component is larger than that of the horizontal component. In summary, in the wave source problem, a small value of modified operator γ can produce a good suppression effect on the drift instability. Whereas, there was an optimal value of γ, and the optimal solution of this numerical case seemed to be around 0.02 and less than 0.05.

Scattering problem

The validity of the proposed measure for the elimination of MTF drift instability is examined by numerical simulation of wave motion on the scattering problem. Scattering problem corresponds to scattering field generated by a concave domain (as shown in Fig 9) under vertical incidence plane S-wave, which is located on the free surface of a homogeneous, isotropic semi-infinite elastic space with a medium density ρ = 2000 kg/m³, shear wave velocity c_s = 1000 m/s, and Poisson's ratio ν = 0.25. The size of the concave domain is 250 m by 250 m by 250 m, and incidence plane S-wave with vibration direction of particles along the x-axis of the Cartesian coordinate system oxyz is an approximate δ-pulse (Fig 5) with 1 s pulse width. The numerical simulation of wave motion is carried out in a computational domain of the homogeneous box with top free surface (Fig 9), and both side and bottom boundaries of the box domain are the artificial boundary, and the solution method is the decoupling numerical simulation method used for solving the wave source problem. The computational domain is discretized by cubic finite elements with a spatial step size of Δx, and a system of ordinary differential equations relative to inner nodes can be formed on this basis. The numerical integration of ordinary differential equations is carried out by the central difference method, and the time step Δt of numerical integration is determined by the stability criterion of the central difference method. The discrete equation of nodes on the artificial boundaries is a second-order MTF. In numerical simulation of wave motion, the space step Δx is 12.5 m, the time step Δt is 0.005 s, and the artificial wave velocity c_a is c_s.

Fig 9

Computational domain in the scattering problem.

Figs 10 and 11 showed the displacement responses at observation points A and B in Fig 9 respectively, and the comparison with the numerical analytical solution, i.e. artificial boundaries are far enough from the domain of interest so that the numerical solution of this domain is not affected by artificial boundaries during whole analytical time. In the scattering problem, all three components of displacement responses of the observation point A and point B showed obvious drift instability when the MTF is without the modified operator. As the value of γ increased from 0.001 to 0.01, the displacement responses of the observation point A and point B become closer to the solution of extended computational domain. However, when the values of γ exceeded 0.01, it showed a worse suppression effect on the drift instability for observation point A and point B, especially in the peak of the displacement responses. The variation of the suppression effect with the change of the value of γ in scattering problem showed the same trend as that of the wave source problem. But the optimal value of γ in the scattering problem seemed to be around 0.01, which was different with that in the source wave problem. From the numerical results of the wave source problem and the scattering problem, it can be seen that the addition of γ can effectively suppress the drift instability of MTF, but the value of γ should not be too small or too large.

Fig 10

Displacement time-histories of observation point A.

Fig 11

Displacement time-histories of observation point B.

It can be seen from the results of two numerical experiments that the proposed measure is effective in eliminating the drift instability of MTF. The curves in the smaller time windows of Figs 7, 8, 10 and 11 also show that introducing $$ \gamma B_0^0$$ into MTF has no influence on the computational precision of numerical simulation of wave motion before the drift instability of MTF appears. Theoretically, γ can be any small positive number, but it should not be too small when considering the existence of error noise in numerical simulation, and neither should it be too large considering the computational precision of the numerical simulation. It can be seen that the addition of γ can effectively suppress the drift instability of MTF, and the range of which is recommend from 0.01 to 0.05 referring to the calculation cases of our work. Though there is no further study on the optimal value of γ, the optimal value should be related to the scale of the model computational domain. The general principle is to minimize the value of γ under the premise of stable implementation of the MTF in numerical simulations of the wave motions.

Conclusions

The Mechanism of drift instability in the numerical implementation of MTF is theoretically analyzed, and it reveals that drift instability is caused by the reason that the GKS criterion is violated in the MTF with a zero frequency and zero wavenumber. To eliminate this instability, a simple measure in the numerical simulation of wave motions, introducing a modified operator $$ \gamma B_0^0$$ into the MTF, was proposed. Due to introduction of the modified operator $$ \gamma B_0^0$$, the absorption mechanism of the medium to the wave in the numerical simulation of wave motion is drawn into the MTF. In addition, through the accuracy analysis based on the calculation of the reflection coefficient and the verification of numerical experiments, the proposed measure is effective in eliminating drift instability of MTF and has slight influence on the accuracy of numerical simulation under the appropriate small value of γ. It is recommended that the value of γ should not be too large or too small, the optimal value varies according to the different calculation models. Based on the numerical simulations we have carried, the trial range of value of is suggested to be in 0.01 to 0.05. Besides, our strategy for eliminating the drift instability of MTF is also suitable for all versions of MTF.