ABSTRACT

In this paper, we firstly derive the stability conditions of high-order staggered-grid schemes for the three-dimensional (3D) elastic wave equation in heterogeneous media based on the energy method. Moreover, the plane wave analysis yields a sufficient and necessary stability condition by the von Neumann criterion in homogeneous case. Numerical computations for 3D wave simulation with point source excitation are given.

Keywords:

3D, Elastic Wave, Inhomogeneous Media, Staggered-Grid Scheme, High-Order, Stability, Energy Estimate

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

$\{\begin{array}{l}\rho \frac{{\partial }^{2}u}{\partial t^2}-\frac{\partial }{\partial x}\left(\left(\lambda +2\mu \right)\frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y}+\lambda \frac{\partial w}{\partial z}\right)-\frac{\partial }{\partial y}\left(\mu \frac{\partial v}{\partial x}+\mu \frac{\partial u}{\partial y}\right)-\frac{\partial }{\partial z}\left(\mu \frac{\partial u}{\partial z}+\mu \frac{\partial w}{\partial x}\right)=f_1,\\ \rho \frac{{\partial }^{2}v}{\partial t^2}-\frac{\partial }{\partial x}\left(\mu \frac{\partial v}{\partial x}+\mu \frac{\partial u}{\partial y}\right)-\frac{\partial }{\partial y}\left(\left(\lambda +2\mu \right)\frac{\partial v}{\partial y}+\lambda \frac{\partial u}{\partial x}+\lambda \frac{\partial w}{\partial z}\right)-\frac{\partial }{\partial z}\left(\mu \frac{\partial v}{\partial z}+\mu \frac{\partial w}{\partial y}\right)=f_2,\\ \rho \frac{{\partial }^{2}w}{\partial t^2}-\frac{\partial }{\partial x}\left(\mu \frac{\partial u}{\partial z}+\mu \frac{\partial w}{\partial x}\right)-\frac{\partial }{\partial y}\left(\mu \frac{\partial v}{\partial z}+\mu \frac{\partial w}{\partial y}\right)-\frac{\partial }{\partial z}\left(\left(\lambda +2\mu \right)\frac{\partial w}{\partial z}+\lambda \frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y}\right)=f_3,\end{array}$ (1)

where $\left(u\mathrm{,}v\mathrm{,}w\right)\left(x\mathrm{,}t\right)$ are the displacement vector at location $x=\left(x\mathrm{,}y\mathrm{,}z\right)$ and time t, $\rho \left(x\right)$ is the density, $\lambda \left(x\right)>0$ and $\mu \left(x\right)\ge 0$ are the Lamé parameters, $f=\left(f_1\mathrm{,}{f}_2\mathrm{,}{f}_3\right)$ is the external force.

Using the stress tensor, we can formulate the above system (1) as a first order in the following ways

$\{\begin{array}{l}\rho \frac{{\partial }^{2}u}{\partial t^2}-\left(\frac{\partial {\tau }^{xx}}{\partial x}+\frac{\partial {\tau }^{xy}}{\partial y}+\frac{\partial {\tau }^{xz}}{\partial z}\right)=0,\\ \rho \frac{{\partial }^{2}v}{\partial t^2}-\left(\frac{\partial {\tau }^{xy}}{\partial x}+\frac{\partial {\tau }^{yy}}{\partial y}+\frac{\partial {\tau }^{yz}}{\partial z}\right)=0,\\ \rho \frac{{\partial }^{2}w}{\partial t^2}-\left(\frac{\partial {\tau }^{xz}}{\partial x}+\frac{\partial {\tau }^{yz}}{\partial y}+\frac{\partial {\tau }^{zz}}{\partial z}\right)=0,\end{array}$ (2)

where

$\{\begin{array}{l}{\tau }^{xx}=\left(\lambda +2\mu \right)\frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y}+\lambda \frac{\partial w}{\partial z},\\ {\tau }^{yy}=\left(\lambda +2\mu \right)\frac{\partial v}{\partial y}+\lambda \frac{\partial u}{\partial x}+\lambda \frac{\partial w}{\partial z},\\ {\tau }^{zz}=\left(\lambda +2\mu \right)\frac{\partial w}{\partial z}+\lambda \frac{\partial u}{\partial x}+\lambda \frac{\partial v}{\partial y},\\ {\tau }^{xy}=\mu \frac{\partial v}{\partial x}+\mu \frac{\partial u}{\partial y},\\ {\tau }^{xz}=\mu \frac{\partial u}{\partial z}+\mu \frac{\partial w}{\partial x},\\ {\tau }^{yz}=\mu \frac{\partial v}{\partial z}+\mu \frac{\partial w}{\partial y}.\end{array}$ (3)

Let $\Delta x$,$\Delta y$ and $\Delta z$ be the spatial steps of $x\mathrm{,}y\mathrm{,}z$ directions respectively. Now discretization of (2) with the second-order accuracy in space gives

$\rho \frac{{\partial }^{2}u}{\partial t^2}\left(i,j,k\right)-\frac{{\tau }_{i+\frac{1}{2},j,k}^{xx}-{\tau }_{i-\frac{1}{2},j,k}^{xx}}{\Delta x}-\frac{{\tau }_{i,j+\frac{1}{2},k}^{xy}-{\tau }_{i,j-\frac{1}{2},k}^{xy}}{\Delta y}-\frac{{\tau }_{i,j,k+\frac{1}{2}}^{xz}-{\tau }_{i,j,k-\frac{1}{2}}^{xz}}{\Delta z}=0,$ (4)

$\rho \frac{{\partial }^{2}v}{\partial t^2}\left(i,j,k\right)-\frac{{\tau }_{i+\frac{1}{2},j,k}^{xy}-{\tau }_{i-\frac{1}{2},j,k}^{xy}}{\Delta x}-\frac{{\tau }_{i,j+\frac{1}{2},k}^{yy}-{\tau }_{i,j-\frac{1}{2},k}^{yy}}{\Delta y}-\frac{{\tau }_{i,j,k+\frac{1}{2}}^{yz}-{\tau }_{i,j,k-\frac{1}{2}}^{yz}}{\Delta z}=0,$ (5)

$\rho \frac{{\partial }^{2}w}{\partial t^2}\left(i,j,k\right)-\frac{{\tau }_{i+\frac{1}{2},j,k}^{xz}-{\tau }_{i-\frac{1}{2},j,k}^{xz}}{\Delta x}-\frac{{\tau }_{i,j+\frac{1}{2},k}^{yz}-{\tau }_{i,j-\frac{1}{2},k}^{yz}}{\Delta y}-\frac{{\tau }_{i,j,k+\frac{1}{2}}^{zz}-{\tau }_{i,j,k-\frac{1}{2}}^{zz}}{\Delta z}=0.$ (6)

For computing this, we need the values of u, v, and w at the grid $\left(i+\frac{1}{2}\mathrm{,}j+\frac{1}{2}\mathrm{,}k+\frac{1}{2}\right)$. One convenient way is to choose averaging the corresponding vales. For example,

$u_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}=\frac{1}{3}\left(\frac{u_{i+1,j,k}+u_{i,j,k}}{2}+\frac{u_{i,j+1,k}+u_{i,j,k}}{2}+\frac{u_{i,j,k+1}+u_{i,j,k}}{2}\right).$

However, such choices have no physical meaning. Another way is to compute $u\mathrm{,}v\mathrm{,}w$ directly on staggered grids. In particular, we replace Equations (4)-(6) with Equations (7)-(9):

$\rho \frac{{\partial }^{2}u}{\partial t^2}\left(i,j,k\right)-\frac{{\tau }_{i+\frac{1}{2},j,k}^{xx}-{\tau }_{i-\frac{1}{2},j,k}^{xx}}{\Delta x}-\frac{{\tau }_{i,j+\frac{1}{2},k}^{xy}-{\tau }_{i,j-\frac{1}{2},k}^{xy}}{\Delta y}-\frac{{\tau }_{i,j,k+\frac{1}{2}}^{xz}-{\tau }_{i,j,k-\frac{1}{2}}^{xz}}{\Delta z}=0,$ (7)

$\begin{array}{l}\rho \frac{{\partial }^{2}v}{\partial t^2}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)-\frac{{\tau }_{i+1,j+\frac{1}{2},k}^{xy}-{\tau }_{i,j+\frac{1}{2},k}^{xy}}{\Delta x}-\frac{{\tau }_{i+\frac{1}{2},j+1,k}^{yy}-{\tau }_{i+\frac{1}{2},j,k}^{yy}}{\Delta y}\\ -\frac{{\tau }_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}^{yz}-{\tau }_{i+\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}}^{yz}}{\Delta z}=0,\end{array}$ (8)

$\begin{array}{l}\rho \frac{{\partial }^{2}w}{\partial t^2}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)-\frac{{\tau }_{i+1,j,k+\frac{1}{2}}^{xz}-{\tau }_{i,j,k+\frac{1}{2}}^{xz}}{\Delta x}-\frac{{\tau }_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}^{yz}-{\tau }_{i+\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}}^{yz}}{\Delta y}\\ -\frac{{\tau }_{i+\frac{1}{2},j,k+1}^{zz}-{\tau }_{i+\frac{1}{2},j,k}^{zz}}{\Delta z}=0.\end{array}$ (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 $D_x$ on half integer grids with $O\left(\Delta x^{2L}\right)$ order as follows:

${\left(\frac{\partial u}{\partial x}\right)}_{i,j,k}\approx D_x{u}_{i,j,k}={\displaystyle \underset{l=1}{\overset{L}{\sum }}}\frac{{\beta }_l}{\Delta x}\left[u_{i+l-\frac{1}{2},j,k}-u_{i-l+\frac{1}{2},j,k}\right],$ (10)

or

${\left(\frac{\partial u}{\partial x}\right)}_{i+\frac{1}{2},j,k}\approx D_x{u}_{i+\frac{1}{2},j,k}={\displaystyle \underset{l=1}{\overset{L}{\sum }}}\frac{{\beta }_l}{\Delta x}\left[u_{i+l,j,k}-u_{i-l+1,j,k}\right],$ (11)

where ${\beta }_l$ 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 $O\left(\Delta x^{2L}\right)$ order accuracy. For example, when ${\beta }_1=1$ for $L=1$ it has the second-order accuracy $O\left(\Delta x^2\right)$. And when ${\beta }_1=\frac{9}{8}$ and ${\beta }_2=\frac{1}{24}$ for $L=2$ it has the fourth-order accuracy $O\left(\Delta x^4\right)$. The general analytical expression of ${\beta }_l$ is given in Appendix. Similarly, we can define the operators $D_y$ and $D_z$. Here, the subscript of the operator refers to the direction of differentiation.

Now, we can construct the semi-discrete schemes of system (1)

$\begin{array}{l}(\rho \frac{{\partial }^{2}u}{\partial t^2}+D_x\left(\left(\lambda +2\mu \right)D_{x}u+\lambda D_{y}v+\lambda D_{z}w\right)+D_y\left(\mu D_{x}v+\mu D_{y}u\right)\\ \begin{array}{c}\\ \end{array}+D_z\left(\mu D_{z}u+\mu D_{x}w\right))\left(i,j,k\right)=0,\end{array}$ (12)

$\begin{array}{l}(\rho \frac{{\partial }^{2}v}{\partial t^2}+D_x\left(\mu D_{x}v+\mu D_{y}u\right)+D_y\left(\left(\lambda +2\mu \right)D_{y}v+\lambda D_{x}u+\lambda D_{z}w\right)\\ \begin{array}{c}\\ \end{array}+D_z\left(\mu D_{z}v+\mu D_{y}w\right))\left(i+\frac{1}{2},j+\frac{1}{2},k\right)=0,\end{array}$ (13)

$\begin{array}{l}(\rho \frac{{\partial }^{2}w}{\partial t^2}+D_x\left(\mu D_{z}u+\mu D_{x}w\right)+D_y\left(\mu D_{z}v+\mu D_{y}w\right)\\ \begin{array}{c}\\ \end{array}+D_z\left(\left(\lambda +2\mu \right)D_{z}w+\lambda D_{x}u+\lambda D_{y}v\right))\left(i+\frac{1}{2}\mathrm{,}j\mathrm{,}k+\frac{1}{2}\right)=0.\end{array}$ (14)

Applying the central difference approximation for time with the second-order accuracy, we obtain the full-discrete schemes of system (1), we obtain

$\begin{array}{l}{\rho }_{i\mathrm{,}j\mathrm{,}k}{\left(u^{n+1}-2u^n+u^{n-1}\right)}_{i\mathrm{,}j\mathrm{,}k}+\Delta t^2\{D_x\left(\left(\lambda +2\mu \right)D_{x}u+\lambda D_{y}v+\lambda D_{z}w\right)\\ +{D_y\left(\mu D_{x}v+\mu D_{y}u\mathrm{<mo\; stretchy="false">\; )\; </mo>}+D_z\left(\mu D_{z}u+\mu D_{x}w\right)\right)\}}_{i\mathrm{,}j\mathrm{,}k}^n=\mathrm{0,}\end{array}$ (15)

$\begin{array}{l}l{\rho }_{i+\frac{1}{2},j+\frac{1}{2},k}{\left(v^{n+1}-2v^n+v^{n-1}\right)}_{i+\frac{1}{2},j+\frac{1}{2},k}+\Delta t^2\{D_x\left(\mu D_{x}v+\mu D_y^{0}u\right)\\ +D_y\left(\left(\lambda +2\mu \right)D_{y}v+\lambda D_{x}u+\lambda D_{z}w\right)+{D_z\left(\mu D_{z}v+\mu D_{y}w\right)\}}_{i+\frac{1}{2},j+\frac{1}{2},k}^n=0,\end{array}$ (16)

$\begin{array}{l}{\rho }_{i+\frac{1}{2},j,k+\frac{1}{2}}{\left(w^{n+1}-2w^n+w^{n-1}\right)}_{i+\frac{1}{2},j,k+\frac{1}{2}}+\Delta t^2\{D_x\left(\mu D_{z}u+\mu D_{x}w\right)\\ +{D_y\left(\mu D_{z}v+\mu D_{y}w\right)+D_z\left(\left(\lambda +2\mu \right)D_{z}w+\lambda D_{x}u+\lambda D_{y}v\right)\}}_{i+\frac{1}{2},j,k+\frac{1}{2}}^n=0,\end{array}$ (17)

where n denotes the time index and $\Delta t$ 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:

$E=E_c+E_p,$ (18)

with

$E_c=\frac{1}{2}{\displaystyle \underset{R^3}{\int }}\rho \left[{\left(\frac{\partial u}{\partial t}\right)}^2+{\left(\frac{\partial v}{\partial t}\right)}^2+{\left(\frac{\partial w}{\partial t}\right)}^2\right]\text{d}x\text{d}y\text{d}z,$

$\begin{array}{l}{E}_p=\frac{1}{2}{\displaystyle \underset{R^3}{\int }}\lambda {\left(\frac{\partial u}{\partial x}+\frac{\partial u}{\partial y}+\frac{\partial u}{\partial z}\right)}^2+2\mu \left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial z}\right)}^2\right]\\ +\mu \left[{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2+{\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)}^2\right]\text{d}x\text{d}y\text{d}z.\end{array}$

In a source-free infinite medium, the energy is conservative, i.e., $\frac{\text{d}E}{\text{d}t}=0$. The discrete energy at time $\left(n+\frac{1}{2}\right)\Delta t$ is $E^{n+\frac{1}{2}}=E_c^{n+\frac{1}{2}}+E_p^{n+\frac{1}{2}}$. In order to compute $E_c^{n+\frac{1}{2}}$ and $E_p^{n+\frac{1}{2}}$ and analyze the stability, we define the following functional spaces

$\begin{array}{l}{L}_{000}^3=\{\phi \in L^3\left(R^3\right)|\phi ={\displaystyle \underset{i,j,k=-\infty }{\overset{\infty }{\sum }}}{\phi }_{i,j,k}\\ \times I_{\left[\left(i-\frac{1}{2}\right)\Delta x,\left(i+\frac{1}{2}\right)\Delta x\right]\times \left[\left(j-\frac{1}{2}\right)\Delta y,\left(j+\frac{1}{2}\right)\Delta y\right]\times \left[\left(k-\frac{1}{2}\right)\Delta z,\left(k+\frac{1}{2}\right)\Delta z\right]}\left(x,y,z\right)\},\end{array}$

$\begin{array}{l}{L}_{\ast \ast \ast }^3=\{\phi \in L^3\left(R^3\right)|\phi ={\displaystyle \underset{i,j,k=-\infty }{\overset{\infty }{\sum }}}{\phi }_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}}\\ \begin{array}{c}\\ \end{array}\times I_{\left[i\Delta x,\left(i+1\right)\Delta x\right]\times \left[j\Delta y,\left(j+1\right)\Delta y\right]\times \left[k\Delta z,\left(k+1\right)\Delta z\right]}\left(x,y,z\right)\},\end{array}$

$\begin{array}{l}{L}_{0\ast \ast }^3=\{\phi \in L^3\left(R^3\right)|\phi ={\displaystyle \underset{i,j,k=-\infty }{\overset{\infty }{\sum }}}{\phi }_{i,j+\frac{1}{2},k+\frac{1}{2}}\\ \times I_{\left[\left(i-\frac{1}{2}\right)\Delta x,\left(i+\frac{1}{2}\right)\Delta x\right]\times \left[j\Delta y,\left(j+1\right)\Delta y\right]\times \left[k\Delta z,\left(k+1\right)\Delta z\right]}\left(x,y,z\right)\},\end{array}$

where

$I_{\left[a,b\right]\times \left[c,d\right]\times \left[e,f\right]}\left(x,y,z\right)=(\begin{array}{ll}1,\hfill & \left(x,y,z\right)\in \left[a,b\right]\times \left[c,d\right]\times \left[e,f\right],\hfill \\ 0,\hfill & \left(x,y,z\right)\notin \left[a,b\right]\times \left[c,d\right]\times \left[e,f\right],\hfill \end{array}$

where 0 represents the integer grid $\left(i\mathrm{,}j\mathrm{,}k\right)$ and $\ast$ the half integer grid $i+\frac{1}{2}$ or $j+\frac{1}{2}$ or $k+\frac{1}{2}$. The other functional spaces $L_{00\ast }^3$,$L_{\ast \ast 0}^3$,$L_{\ast 00}^3$ $L_{\mathrm{<mo>*</mo>00}}^3$ , $L_{\ast 0\ast }^3$,$L_{0\ast \ast }^3$ can be defined similarly. For saving space, we omit their definitions. Let the scalar inner product be defined in $L_{000}^3$ by

${\left(f,g,h\right)}_{000}={\displaystyle \underset{i,j,k=-\infty }{\overset{\infty }{\sum }}}{f}_{i,j,k}{g}_{i,j,k}{h}_{i,j,k}\Delta x\Delta y\Delta z.$

Other inner products such as ${\left(f\mathrm{,}g\mathrm{,}h\right)}_{\ast \ast \ast }$,${\left(f\mathrm{,}g\mathrm{,}h\right)}_{0\ast \ast }$ and so on have similar meaning. In the following, we compute $E_c^{n+\frac{1}{2}}$ and $E_p^{n+\frac{1}{2}}$ respectively.

$\begin{array}{c}{E}_c^{n+\frac{1}{2}}=\frac{1}{2}\{{\left(\rho \frac{u^{n+1}-u^n}{\Delta t}\mathrm{,}\frac{u^{n+1}-u^n}{\Delta t}\right)}_{000}+{\left(\rho \frac{v^{n+1}-v^n}{\Delta t}\mathrm{,}\frac{v^{n+1}-v^n}{\Delta t}\right)}_{\ast \ast 0}\\ +{\left(\rho \frac{w^{n+1}-w^n}{\Delta t}\mathrm{,}\frac{w^{n+1}-w^n}{\Delta t}\right)}_{0\ast 0}-\frac{\Delta t^2}{4}[(\lambda D_x\frac{u^{n+1}-u^n}{\Delta t}\\ +\lambda D_y\frac{v^{n+1}-v^n}{\Delta t}+\lambda D_z\frac{w^{n+1}-w^n}{\Delta t}\mathrm{,}{D}_x\frac{u^{n+1}-u^n}{\Delta t}+D_y\frac{v^{n+1}-v^n}{\Delta t}\\ +{D_z\frac{w^{n+1}-w^n}{\Delta t})}_{\ast 00}+(\mu D_x\frac{v^{n+1}-v^n}{\Delta t}+\mu D_y\frac{u^{n+1}-u^n}{\Delta t}\mathrm{,}\\ D_x\frac{v^{n+1}-v^n}{\Delta t}+{D_y\frac{u^{n+1}-u^n}{\Delta t})}_{0\ast 0}+(\mu D_z\frac{u^{n+1}-u^n}{\Delta t}\end{array}$

$\begin{array}{c}+\mu D_x\frac{w^{n+1}-w^n}{\Delta t}\mathrm{,}{D}_z\frac{u^{n+1}-u^n}{\Delta t}+{D_x\frac{w^{n+1}-w^n}{\Delta t})}_{00\ast }\\ +(\mu D_y\frac{w^{n+1}-w^n}{\Delta t}+\mu D_z\frac{v^{n+1}-v^n}{\Delta t}\mathrm{,}\mu D_y\frac{w^{n+1}-w^n}{\Delta t}+{\mu D_z\frac{v^{n+1}-v^n}{\Delta t})}_{\ast \ast \ast }\\ +{\left(2\mu D_x\frac{u^{n+1}-u^n}{\Delta t}\mathrm{,}{D}_x\frac{u^{n+1}-u^n}{\Delta t}\right)}_{\ast 00}+(2\mu D_y\frac{v^{n+1}-v^n}{\Delta t}\mathrm{,}\\ {D_y\frac{v^{n+1}-v^n}{\Delta t})}_{\ast 00}+{\left(2\mu D_z\frac{w^{n+1}-w^n}{\Delta t}\mathrm{,}{D}_z\frac{w^{n+1}-w^n}{\Delta t}\right)}_{\ast 00}]\}\mathrm{,}\end{array}$

$\begin{array}{c}{E}_p^{n+\frac{1}{2}}=\frac{1}{2}\{(\lambda D_x\frac{u^{n+1}+u^n}{2}+\lambda D_y\frac{v^{n+1}+v^n}{2}+\lambda D_z\frac{w^{n+1}+w^n}{2},\\ {D_x\frac{u^{n+1}+u^n}{2}+D_y\frac{v^{n+1}+v^n}{2}+D_z\frac{w^{n+1}+w^n}{2})}_{\ast 00}\end{array}$

$\begin{array}{c}+{\left(2\mu D_x\frac{u^{n+1}+u^n}{2},D_x\frac{u^{n+1}+u^n}{2}\right)}_{\ast 00}+(2\mu D_y\frac{v^{n+1}+v^n}{2},\\ {D_y\frac{v^{n+1}+v^n}{2})}_{0*0}+{\left(2\mu D_z\frac{w^{n+1}+w^n}{2},D_z\frac{w^{n+1}+w^n}{2}\right)}_{\ast 00}\end{array}$

$\begin{array}{c}+{\left(\mu D_y\frac{u^{n+1}+u^n}{2}+\mu D_x\frac{v^{n+1}+v^n}{2},D_y\frac{u^{n+1}+u^n}{2}+D_x\frac{v^{n+1}+v^n}{2}\right)}_{0\ast 0}\\ +{\left(\mu D_z\frac{u^{n+1}+u^n}{2}+\mu D_x\frac{w^{n+1}+w^n}{2},D_z\frac{u^{n+1}+u^n}{2}+D_x\frac{w^{n+1}+w^n}{2}\right)}_{00\ast }\\ +{\left(\mu D_z\frac{v^{n+1}+v^n}{2}+\mu D_y\frac{w^{n+1}+w^n}{2},D_z\frac{v^{n+1}+v^n}{2}+D_y\frac{w^{n+1}+w^n}{2}\right)}_{\ast \ast \ast }\}.\end{array}$

We have conservation of the discrete energy, that is: $\frac{E^{n+\frac{1}{2}}-E^{n-\frac{1}{2}}}{\Delta t}=0$. The stability of the scheme will be proven if the potential energy $E_p^{n+\frac{1}{2}}$ and the kinetic energy $E_c^{n+\frac{1}{2}}$ are positive. Since $E_p^{n+\frac{1}{2}}$ is obviously positive, we need to find out under what conditions $E_c^{n+\frac{1}{2}}$ is positive.

The problem can be reformulated as: $\forall u\in L_{000}^3$,$\forall v\in L_{\ast \ast 0}^3$ and $\forall w\in L_{\ast 0\ast }^3$ with

$\begin{array}{c}I={\left(\lambda D_{x}u+\lambda D_{y}v+\lambda D_{z}w\mathrm{,}{D}_{x}u+D_{y}v+D_{z}w\right)}_{\ast 00}\\ +{\left(\mu D_{x}v+\mu D_{y}u\mathrm{,}{D}_{x}v+D_{y}u\right)}_{0\ast 0}+{\left(\mu D_{z}u+\mu D_{x}w\mathrm{,}{D}_{z}u+D_{x}w\right)}_{00\ast }\\ +{\left(\mu D_{y}w+\mu D_{z}v\mathrm{,}{D}_{y}w+D_{z}v\right)}_{\ast \ast \ast }+{\left(2\mu D_{x}u\mathrm{,}{D}_{x}u\right)}_{\ast 00}\\ +{\left(2\mu D_{y}v\mathrm{,}{D}_{y}v\right)}_{\ast 00}+{\left(2\mu D_{z}w\mathrm{,}{D}_{z}w\right)}_{\ast 00}\end{array}$

we look for the corresponding conditions. Since

$\frac{\Delta t^2}{4}I\le {\left(\rho u\mathrm{,}u\right)}_{000}+{\left(\rho v\mathrm{,}v\right)}_{\ast \ast 0}+{\left(\rho w\mathrm{,}w\right)}_{\ast 0\ast }\mathrm{,}$ (19)

we can bound I as follows:

$\begin{array}{c}I\le 2[{\left(\lambda D_{x}u,D_{x}u\right)}_{\ast 00}+{\left(\lambda D_{y}v,D_{y}v\right)}_{\ast 00}+{\left(\lambda D_{z}w,D_{z}w\right)}_{\ast 00}\\ +{\left(\mu D_{x}v,D_{x}v\right)}_{0\ast 0}+{\left(\mu D_{y}u,D_{y}u\right)}_{0\ast 0}+{\left(\mu D_{x}w,D_{x}w\right)}_{00\ast }\\ +{\left(\mu D_{z}u,D_z^{0}u\right)}_{00\ast }+{\left(\mu D_{z}v,D_{z}v\right)}_{\ast \ast \ast }+{\left(\mu D_{y}w,D_{y}w\right)}_{\ast \ast \ast }\\ +{\left(\mu D_{x}u,D_{x}u\right)}_{\ast 00}+{\left(\mu D_{y}v,D_y^{*}v\right)}_{\ast 00}+{\left(\mu D_{z}w,D_{z}w\right)}_{\ast 00}],\end{array}$

or

$\begin{array}{c}I\le 2[\underset{I_1}{\underset{︸}{{\left(\left(\lambda +\mu \right)D_{x}u,D_{x}u\right)}_{\ast 00}+{\left(\mu D_{y}u,D_{y}u\right)}_{0\ast 0}+{\left(\mu D_{z}u,D_{z}u\right)}_{00\ast }}}\\ +\underset{I_2}{\underset{︸}{{\left(\left(\lambda +\mu \right)D_{y}v,D_{y}v\right)}_{\ast 00}+{\left(\mu D_{x}v,D_{x}v\right)}_{0\ast 0}+{\left(\mu D_{z}v,D_{z}v\right)}_{\ast \ast \ast }}}\end{array}$

$+\underset{I_3}{\underset{︸}{{\left(\left(\lambda +\mu \right)D_{z}w,D_{z}w\right)}_{\ast 00}+{\left(\mu D_{x}w,D_{x}w\right)}_{00\ast }+{\left(\mu D_{y}w,D_{y}w\right)}_{\ast \ast \ast }}}].$

Obviously, $I\le 2\left(I_1+I_2+I_3\right)$. In the following we estimate $I_1$,$I_2$ and $I_3$ respectively. Since

$\begin{array}{c}{I}_1={\displaystyle \underset{i,j,k}{\sum }}{\left(\lambda +\mu \right)}_{i+\frac{1}{2},j,k}{\left[{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\frac{{\beta }_l}{\Delta x}\left(u_{i+l,j,k}-u_{i-l+1,j,k}\right)\right]}^2\Delta x\Delta y\Delta z\\ +{\displaystyle \underset{i,j,k}{\sum }}{\mu }_{i,j+\frac{1}{2},k}{\left[{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\frac{{\beta }_l}{\Delta y}\left(u_{i,j+l,k}-u_{i,j-l+1,k}\right)\right]}^2\Delta x\Delta y\Delta z\\ +{\displaystyle \underset{i,j,k}{\sum }}{\mu }_{i,j,k+\frac{1}{2}}{\left[{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\frac{{\beta }_l}{\Delta z}\left(u_{i,j,k+l}-u_{i,j,k-l+1}\right)\right]}^2\Delta x\Delta y\Delta z,\end{array}$

we have the following estimates for $I_1$ :

$\begin{array}{c}{I}_1\le 2\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta x}\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta x}\right|{\left(\lambda +\mu \right)}_{i+\frac{1}{2},j,k}\left(u_{i+l,j,k}^2+u_{i-l+1,j,k}^2\right)\Delta x\Delta y\Delta z\\ +2\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta y}\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta y}\right|{\mu }_{i,j,k+\frac{1}{2}}\left(u_{i,j,k+l}^2+u_{i,j,k-l+1}^2\right)\Delta x\Delta y\Delta z\\ +2\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta z}\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta z}\right|{\mu }_{i,j+\frac{1}{2},k}\left(u_{i,j+l,k}^2+u_{i,j-l+1,k}^2\right)\Delta x\Delta y\Delta z,\end{array}$

or

$\begin{array}{c}{I}_1\le 4\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta x}\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta x}\right|\left[\frac{{\left(\lambda +\mu \right)}_{i+l-\frac{1}{2},j,k}}{2{\rho }_{i,j,k}}+\frac{{\left(\lambda +\mu \right)}_{i-l+\frac{1}{2},j,k}}{2{\rho }_{i,j,k}}\right]{\rho }_{i,j,k}{u}_{i,j,k}^2\Delta x\Delta y\Delta z\\ +4\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta z}\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta z}\right|\left[\frac{{\mu }_{i,j,k+l-\frac{1}{2}}}{2{\rho }_{i,j,k}}+\frac{{\mu }_{i,j,k-l+\frac{1}{2}}}{2{\rho }_{i,j,k}}\right]{\rho }_{i,j,k}{u}_{i,j,k}^2\Delta x\Delta y\Delta z\\ +4\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta y}\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|\frac{{\beta }_l}{\Delta y}\right|\left[\frac{{\mu }_{i,j+l-\frac{1}{2},k}}{2{\rho }_{i,j,k}}+\frac{{\mu }_{i,j-l+\frac{1}{2},k}}{2{\rho }_{i,j,k}}\right]{\rho }_{i,j,k}{u}_{i,j,k}^2\Delta x\Delta y\Delta z,\end{array}$

or

$\begin{array}{c}{I}_1\le \frac{4}{\Delta x^2}\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|{\beta }_l\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|{\beta }_l\right|\left[\frac{{\left(\lambda +\mu \right)}_{i+l-\frac{1}{2},j,k}}{2{\rho }_{i,j,k}}+\frac{{\left(\lambda +\mu \right)}_{i-l+\frac{1}{2},j,k}}{2{\rho }_{i,j,k}}\right]{\left(\rho u,u\right)}_{000}\\ +\frac{4}{\Delta z^2}\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|{\beta }_l\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|{\beta }_l\right|\left[\frac{{\mu }_{i,j,k+l-\frac{1}{2}}}{2{\rho }_{i,j,k}}+\frac{{\mu }_{i,j,k-l+\frac{1}{2}}}{2{\rho }_{i,j,k}}\right]{\left(\rho u,u\right)}_{000}\\ +\frac{4}{\Delta y^2}\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|{\beta }_l\right|\right){\displaystyle \underset{i,j,k}{\sum }}{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|{\beta }_l\right|\left[\frac{{\mu }_{i,j+l-\frac{1}{2},k}}{2{\rho }_{i,j,k}}+\frac{{\mu }_{i,j-l+\frac{1}{2},k}}{2{\rho }_{i,j,k}}\right]{\left(\rho u,u\right)}_{000}.\end{array}$

Thus we obtain

$\begin{array}{l}{I}_1\le \left(\frac{4}{\Delta x^2}+\frac{4}{\Delta y^2}+\frac{4}{\Delta z^2}\right)\left({\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|{\beta }_l\right|\right)\underset{i,j,k}{\max }{\displaystyle \underset{l=1}{\overset{L}{\sum }}}\left|{\beta }_l\right|\times [\frac{{\left(\lambda +\mu \right)}_{i+l-\frac{1}{2},j,k}}{2{\rho }_{i,j,k}}\\ +\frac{{\left(\lambda +\mu \right)}_{i-l+\frac{1}{2},j,k}}{2{\rho }_{i,j,k}}+\frac{{\mu }_{i,j,k+l-\frac{1}{2}}}{2{\rho }_{i,j,k}}+\frac{{\mu }_{i,j,k-l+\frac{1}{2}}}{2{\rho }_{i,j,k}}+\frac{{\mu }_{i,j+l-\frac{1}{2},k}}{2{\rho }_{i,j,k}}+\frac{{\mu }_{i,j-l+\frac{1}{2},k}}{2{\rho }_{i,j,k}}]{\left(\rho u,u\right)}_{000}.\end{array}$