Stability of High-Order Staggered-Grid Schemes for 3D Elastic Wave Equation in Heterogeneous Media ()
1. Introduction
Numerical simulation of wave propagation has important applications in many scientific fields such as geophysics and seismic inversion. There are several types of numerical methods to solve the wave equations, for example, the finite difference method, the finite element method [1] [2] [3] [4] , the spectral element method [5] , the discontinuous Galerkin method [6] [7] and the finite volume method [8] [9] . Each of the above numerical methods has its own advantages and disadvantages. In this paper, we consider the finite difference method.
The finite difference method is a very popular method because of high computational efficiency. In fact, it has been applied to wave simulation for several decades [10] [11] [12] [13] . Since perfect numerical simulation depends on both stability and the order of accuracy, the high-order schemes and the corresponding stability are an important research topic of this field. In particularly, we may list a few here. In [14] , Cohen and Joly construct and analyses a family of fourth-order schemes for the acoustic wave equation in nonhomogeneous media. In [15] , Sei analysis the stability of high-order difference schemes for the 2D elastic wave equation in heterogeneous media. The stable difference approximation for the 3D elastic wave equation in the second-order formulation in heterogeneous media has been investigated in [16] . In [17] , a new family of locally one-dimensional schemes with fourth-order accuracy both in space and time for the 3D elastic wave equation is constructed and the stability is derived. The constructed new schemes in [17] only involve a three-point stencil in each spatial direction to achieve fourth-order accuracy. In this paper, based on the energy method, we study the stability analysis for the high-order staggered-grid schemes of the 3D elastic wave equation in heterogeneous media. To our knowledge, there is no work in this respect and our result is new.
The reminder of the paper is organized as follows. In Section 2, we present the governing equation and the high-order difference schemes in heterogeneous on staggered-grid grids. In Section 3, the stability analysis for the high-order difference schemes in heterogeneous is presented. In Section 4, the plane wave analysis in homogeneous case is investigated. In Section 5, we present numerical comparisons for 3D elastic wave simulation. Finally the conclusion and discussions are given in Section 6.
2. High-Order Spatial Discretization
We consider the following three-dimensional (3D) elastic wave equations in isotropic heterogeneous media
(1)
where
are the displacement vector at location
and time t,
is the density,
and
are the Lamé parameters,
is the external force.
Using the stress tensor, we can formulate the above system (1) as a first order in the following ways
(2)
where
(3)
Let
,
and
be the spatial steps of
directions respectively. Now discretization of (2) with the second-order accuracy in space gives
(4)
(5)
(6)
For computing this, we need the values of u, v, and w at the grid
. One convenient way is to choose averaging the corresponding vales. For example,
However, such choices have no physical meaning. Another way is to compute
directly on staggered grids. In particular, we replace Equations (4)-(6) with Equations (7)-(9):
(7)
(8)
(9)
Obviously, the schemes (7)-(9) are the second-order accuracy in space. In order to construct high-order accuracy scheme in space, we first define the following functional spaces. Now we introduce the differentiation operator
on half integer grids with
order as follows:
(10)
or
(11)
where
is difference coefficients on the staggered grids, which can be calculated by a fast algorithm [18] [19] [20] by Matlab tool. Obviously, the approximation (10) or (11) has
order accuracy. For example, when
for
it has the second-order accuracy
. And when
and
for
it has the fourth-order accuracy
. The general analytical expression of
is given in Appendix. Similarly, we can define the operators
and
. Here, the subscript of the operator refers to the direction of differentiation.
Now, we can construct the semi-discrete schemes of system (1)
(12)
(13)
(14)
Applying the central difference approximation for time with the second-order accuracy, we obtain the full-discrete schemes of system (1), we obtain
(15)
(16)
(17)
where n denotes the time index and
the time step.
3. Stability Analysis
We now turn to the study of the numerical stability of the schemes (15)-(17). We are going to proceed by the energy method in analogy with continuous energy given by:
(18)
with
In a source-free infinite medium, the energy is conservative, i.e.,
. The discrete energy at time
is
. In order to compute
and
and analyze the stability, we define the following functional spaces
where
where 0 represents the integer grid
and
the half integer grid
or
or
. The other functional spaces
,
,
,
,
can be defined similarly. For saving space, we omit their definitions. Let the scalar inner product be defined in
by
Other inner products such as
,
and so on have similar meaning. In the following, we compute
and
respectively.
We have conservation of the discrete energy, that is:
. The stability of the scheme will be proven if the potential energy
and the kinetic energy
are positive. Since
is obviously positive, we need to find out under what conditions
is positive.
The problem can be reformulated as:
,
and
with
we look for the corresponding conditions. Since
(19)
we can bound I as follows:
or
Obviously,
. In the following we estimate
,
and
respectively. Since
we have the following estimates for
:
or
or
Thus we obtain
Set
(20)
Then for
, we have
(21)
Similarly, we have
(22)
where
(23)
And
(24)
where
(25)
Substituting (21), (22) and (24) into (19), we have
(26)
where
. Thus we obtain the sufficient stability condition for the numerical scheme (15)-(17). If the grid is uniform, i.e.
, then (26) gives
Therefore we summarize the conclusion above into the following theorem.
Theorem 1. A sufficient stability condition for the numerical schemes (15)-(17) is
(27)
If
, then it reduces to
(28)
where
, and
,
and
are given by (20), (23) and (25) respectively.
4. Plane Wave Analysis
We turn to Fourier analysis [21] and we will derive the dispersion relation and by the von Neumann criterion we will get a necessary and sufficient stability condition. In homogeneous case for (12)-(14), the full-discrete schemes can be written as
(29)
(30)
(31)
We assume that
is a solution of Equation (29)-(31), where
,
is the angular frequency,
amplitude, and
is the wave vector. Here
is the propagation angle and
the propagation azimuth. The two angles determine the movement direction of the plane wave in the 3D space.
Substituting the plane wave solution into Equations (29)-(31), we obtain the following relations:
(32)
(33)
(34)
By introducing the matrix B with elements
defined by
we can write the relations (32)-(34) as the following matrix form
(35)
The eigenvalues of B then express
as a function of
, which is the dispersion relation. There are three eigenvalues for matrix B. One eigenvalue is corresponding to the longitudinal or compressional wave, the double eigenvalues are corresponding to the transverse or shear wave. Thus we have the following two different relations
(36)
where
here
and
are the velocities of compressional and shear waves. Note that
is always larger than
. With the dispersion relations (36), we can apply the von Neumann stability criterion. A necessary stability is that the eigenvalues of B must be lower than 1. Thus we have
(37)
It is easy to verify that
(38)
Therefore we obtain the following theorem.
Theorem 2. In the homogeneous case, a sufficient and necessary stability condition for the numerical schemes (29)-(31) is given by
(39)
Proof Combining Equations (37) and (38), we obtain (39). Moreover, the matrix B is symmetric in homogeneous case. So the condition (39) is a sufficient and necessary condition. The proof is complete.
We now define the normalized phase error
as follow:
(40)
which is a function of
, where
indicates
or
which is related to different kinds of compressional wave and shear wave.
The stability condition
is defined by Courant-Friedrichs-Lewy (CFL) condition which bounds the interval for stability. We plot some dispersion curves based on Equation (40). Without loss of generality, we present dispersion curves for some special propagation angle and azimuth. Figure 1 is the normalized phase error for fixed
and
with different values of CFL condition and it shows that the phase error drops as increasing the order of accuracy. Figure 2 shows the normalized phase error for
and different values of
for different order or L. The figures for other propagation angle
and azimuth
are similar we omit them for saving space.
5. Numerical Computations
Wave simulation ignited by a point source is usually adopted in geophysical applications. For convenience, we simulate 3D elastic wave propagation in a homogeneous cubic model. The computational domain is
. The source is located in the center of the model and its time function is given by
(41)
which is loaded on the u component. The compressional velocity is 4000 m/s and the shear velocity 2500 m/s. The time step is
and the space step is
. Figure 3 shows the 3D snapshot of u component at propagation time 0.2 s. For brevity, we present some 2D slices of the 3D snapshots of u, v, and w components. The xz sections of 3D snapshots of u, v, and w components at propagation time 0.2 s are shown in Figures 4-6 respectively. We omit other sections for space. In our computations the scheme with fourth-order accuracy in space is applied. We remark that the comparisons between the numerical solution and the exact solution can be found in [17] . From Figures 4-6, we can clearly see the two types of waves, i.e. the compressional wave and the shear wave, which is consistent with the physical phenomenon.
Figure 3. The 3D snapshot of u component at propagation time 0.2 s.
Figure 4. The xz-section of 3D snapshot of u component at propagation time 0.2 s.
Figure 5. The xz-section of 3D snapshot of v component at propagation time 0.2 s.
Figure 6. The xz-section of 3D snapshot of w component at propagation time 0.2 s.
6. Conclusion
The staggered-grid difference method is a very important technique to solve wave equations numerically because of its high efficiency and the character of energy preservation. It has been well applied to seismic wave propagation for more than two decades. Based on the energy estimate method, we implement the stability analysis for the high-order staggered-grid schemes of the inhomogeneous 3D elastic wave equation. The stability result is controlled by the space varying parameters and the difference coefficients. The plane wave analysis in homogeneous media is completed and by the von Neumann criterion a necessary and sufficient stability condition is obtained. The analysis is helpful to design the computational parameters such as the time step and the space steps. Numerical computations are given to verify the effectiveness of the schemes. The key point of this paper is the theoretical analysis. In the future, we will consider more numerical computations for inhomogeneous media.
Acknowledgements
This work is supported by National Natural Science Foundation of China under grant numbers 11471328 and 51739007. It is also partially supported by National Center for Mathematics and Interdisciplinary Sciences, Chinese Academy of Sciences.
Appendix: Expression of the Coefficient
The calculation of difference coefficients on both regular and staggered-grid grids has been investigated by several authors [15] [18] [19] [20] . In this appendix we present the analytical expression of the difference coefficient
in (10) or (11). In order to calculate the coefficient
, we consider an explicit staggered-grid difference expression for a function
:
where h is the step size, L is a positive number and
are the difference coefficients. Now consider
and
, where k is the wave number,
and
is a constant then we have
Then Taylor’s series expansion gives
Now equating the coefficient of
both sides we get
We can rewrite this equation in the following form
Now solving the above system, we get the following solutions