A stable implementation measure of multi-transmitting formula in the numerical simulation of wave motion

The Multi-Transmitting Formula (MTF) proposed by Liao et al. is a local artificial boundary condition widely used in numerical simulations of wave propagation in an infinite medium, while the drift instability is usually caused in its numerical implementation. In view of a physical interpretation of the Gustafsson, Kreiss and Sundström criterion on numerical solutions of initial-boundary value problems in the hyperbolic partial differential equations, the mechanism of the drift instability of MTF was discussed, and a simple measure for eliminating the drift instability was proposed by introducing a modified operator into the MTF. Based on the theory of spherical wave propagation and damping effect of medium, the physical implication on modified operator was interpreted. And the effect of the modified operator on the reflection coefficient of MTF was discussed. Finally, the validity of the proposed stable implementation measure was verified by numerical tests of wave source problem and scattering problem.


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][2][3][4][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][7][8]. Numerous studies have been conducted on artificial boundary conditions since the late 1960s [6,[9][10][11][12][13][14][15][16][17][18][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][23][24][25][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 a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 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][28][29][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][33][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 gB 0 0 into the MTF. Moreover, the physical implication of gB 0 0 was explained by the spherical wave propagation principle and the medium damping effect; then, the effect of the gB 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.
Based on the simulation of wave motion, the MTF can be expressed as [24]: where N is the MTF order, u p j ¼ uðpDt; À jc a DtÞ, p and j are arbitrary integers, Δt is the time step, c a is the artificial wave velocity, and C N j ¼ N!=½ðN À jÞ!j!�. The backward moving operator B n m is defined as [24]: The operator satisfies the following operation rules: Using the backward moving operator, the MTF can be expressed as: where B 0 0 is the unit operator.
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]: where ω and k are the circular frequency of the steady-state one-way wave and the apparent wavenumber along the x-axis, respectively, and they are real numbers with the same signs. 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 n m , 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: where γ is a small positive value close to 0. It is easy to understand that this measure of numerical implementation in the MTF guarantees that the GKS criterion is met in the case with a zero frequency and zero wavenumbers. At the same time, the small value of γ incurs only a slight effect on the computational precision of the numerical simulations over the entire computational domain.
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, gB 0 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 0 , and the distance between the observation point o and point source S is r 0 +Δ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: where α is the medium geometry diffusion factor, and α = r 0 /(r 0 +Δr) = 1/(1+Δr/r 0 ). Assuming � g ¼ Dr=r, then we had a ¼ 1=ð1 þ � gÞ. According to operation rules of the moving operator B n m , Eq (7) can be expanded as: It is known through comparison with the MTF that according to the MTF with gB 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.
For the sake of simplicity but without a loss of generality, let the MTF order N = 2 in Eq (9), and thus: Eq (10) can be rewritten as: After simplification, Eq (11) can also be expressed as: where:ũ Eq (12) is another expression of the second-order MTF with gB 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

PLOS ONE
wave propagate in the form of spherical waves with α = 1/1+γ. Therefore, it can be concluded that the MTF with gB 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 gB 0 0 . For the condition of N = 2, we can obtain [6]: where exp(−μp) is the introduced damping factor. According to the assumption in Eq (14), Equation (13) can be rewritten as: 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 gB 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 gB 0 0 described above, the MTF with gB 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 gB 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 γB 0 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]: where U I 0 and U R 0 are the amplitudes of the incident plane harmonic wave and the reflected plane harmonic wave, respectively, at the boundary node. If the total displacement of the boundary node, U 0 , is expressed as Next, the reflection coefficients of the MTF with gB 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), , 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 I j and U I 0 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 I j . According to Eq (20), we can have: Substituting Eq (23) into Eq (19) results in: where ðÀ 1Þ jþ1 C N j a À j a j x =ð1 þ gÞ j . Using the binomial Equation 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 I j is the incident wave displacement, given by Eq (20), and U R j 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: 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 gB 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.

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 1 x 2 x 3 , 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 3 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 1 , B 2 and B 3 (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: 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

PLOS ONE
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 semiinfinite elastic space with a medium density ρ = 2000 kg/m 3 , 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

PLOS ONE
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.
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 gB 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 gB 0 0 into the MTF, was proposed. Due to introduction of the modified operator gB 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.