Equivalent Staggered Grid Acoustic Wave Equation

Equivalent Staggered grid  (Di Bartolo et al., 2012) comes from the first order system of equations that come from the continuity equation and Newton's Second Law. First the vectorial linear equation of the motion for the particle velocity field $\mathbf{v}$ is given by:

\begin{equation}
\rho( \mathbf{r}) \frac{\partial \mathbf{v}(\mathbf{r},t)}{\partial t} + \mathbf{\nabla} p(\mathbf{r},t) = \mathbf{f}(\mathbf{r},t)
\label{s1}
\end{equation}

where $\mathbf{r}$ is the position vector, $t$ is the time, $\rho$ os the density of the medium, $\mathbf{f}$ is the density of external body force applied and $\mathbf{\nabla} $ is the gradient operator.
The equivalent scalar equation for pressure $p$ is given by:

\begin{equation}
\frac{1}{\kappa( \mathbf{r})} \frac{\partial p(\mathbf{r},t)}{\partial t} + \mathbf{\nabla} \cdot \mathbf{v}(\mathbf{r},t) = \frac{\partial i_v(\mathbf{r},t)}{\partial t}
\label{s2}
\end{equation}

where $\kappa$ is the adiabatic compression modulus of the medium and $i_v$ , which represents a source, is a volume density injection (an air gun, for example). The $\mathbf{\nabla}$ applies the divergence to the vectorial velocity field $\mathbf{v} $

The Staggered Grid schemes have the great advantage of being able to model regions of water-rock interface, allowing application in the presence of materials with any Poisson's ration (including water), even when density contrasts are considered. In summary are able to model coupling between elastic and acoustic media.

Di Bartolo et. al. 2011 gives an overview of the Staggered Grid methods and its advantages against non-Staggered Methods. Finally its own formulation defined as Equivalent Staggered Grid is based on the following versions of equations (1) and (2), manipulated to remove the vectorial field $\mathbf{v} $.

\begin{equation}
\rho( \mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) - \frac{1}{c^2(\mathbf{r})} \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t} = -s(\mathbf{r},t)
\label{s3}
\end{equation}

with

\begin{equation}
s(\mathbf{r},t) = \rho( \mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t}-\rho( \mathbf{r}) \mathbf{\nabla} \cdot \left(\frac{1}{\rho(\mathbf{r})} \mathbf{f}(\mathbf{r},t)\right)
\label{s4}
\end{equation}

and with the propagation velocity

\begin{equation}
c(\mathbf{r}) = \sqrt{\frac{\kappa(\mathbf{r})}{\rho(\mathbf{r})}}
\label{s5}
\end{equation}

For practical implementation the source term $s(\mathbf{r},t)$ is simplified to $ \rho( \mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t} $ using a null density force $ \mathbf{f}(\mathbf{r},t) $

So the equation (3) can be written as :

\begin{eqnarray}
\rho( \mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) + s(\mathbf{r},t) &=& \frac{1}{c^2(\mathbf{r})} \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t} \nonumber \\
c^2(\mathbf{r}) \rho( \mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) + c^2(\mathbf{r}) \rho( \mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t} &=& \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t} \\
\kappa(\mathbf{r}) \mathbf{\nabla} \cdot \left( \frac{1}{\rho (\mathbf{r})} \mathbf{\nabla} p(\mathbf{r},t) \right) + \kappa(\mathbf{r}) \frac{\partial^2 i_v(\mathbf{r},t)}{\partial^2 t} &=& \frac{\partial^2 p(\mathbf{r},t)}{\partial^2 t}
\end{eqnarray}

From (6) the finite differences marching scheme is defined using the grid properties are $ c(\mathbf{r}) $ and $ \rho( \mathbf{r}) $ and $i_v(\mathbf{r},t)$ can be any combination of punctual sources. The schema is of 2nd order in time using central differences for the second derivative of pressure related to time

For each dimension the divergence term in (6) is written as:

$$ \mathbf{\nabla} \cdot \left( b (\mathbf{r}) \mathbf{\nabla} p(\mathbf{r},t) \right) = \frac{\partial}{\partial x}\left( b (\mathbf{r}) \frac{\partial p}{\partial x}\right)$$

Translated to finite differences and equivalent staggered grid takes the form:

$$ \frac{\partial}{\partial x}\left( b (\mathbf{r}) \frac{\partial p}{\partial x}\right)^n_{i} $$
with $b$ called buoyancy the inverse of density.

$$
=\frac{9}{8} \left[ \frac{b_{i+\frac{1}{2}}}{\Delta x} \left( \frac{9}{8} \frac{ p^n_{i+1}-p^n_i }{\Delta x} - \frac{1}{8} \frac{ p^n_{i+2}-p^n_{i-1} }{3\Delta x}\right) -\frac{b_{i-\frac{1}{2}}}{\Delta x} \left( \frac{9}{8} \frac{ p^n_{i}-p^n_{i-1}}{\Delta x} - \frac{1}{8} \frac{ p^n_{i+1}-p^n_{i-2} }{3\Delta x}\right) \right]
$$
$$
-\frac{1}{8} \left[ \frac{b_{i+\frac{3}{2}}}{3\Delta x} \left( \frac{9}{8} \frac{ p^n_{i+2}-p^n_{i+1} }{\Delta x} - \frac{1}{8} \frac{ p^n_{i+3}-p^n_{i} }{3\Delta x}\right) -\frac{b_{i-\frac{3}{2}}}{3\Delta x} \left( \frac{9}{8} \frac{ p^n_{i-1}-p^n_{i-2}}{\Delta x} - \frac{1}{8} \frac{ p^n_{i}-p^n_{i-3} }{3\Delta x}\right) \right]
$$

better form by Andre

$$
=\frac{9}{64 \Delta x ^2} \left[ b_{i+\frac{1}{2}} \left( 27 \frac{ p^n_{i+1}-p^n_i }{3} - \frac{ p^n_{i+2}-p^n_{i-1} }{3}\right) -b_{i-\frac{1}{2}} \left( 27 \frac{ p^n_{i}-p^n_{i-1}}{3} - \frac{ p^n_{i+1}-p^n_{i-2} }{3}\right) \right]
$$
$$
-\frac{1}{64 \Delta x ^2} \left[ \frac{b_{i+\frac{3}{2}}}{3} \left( 27 \frac{ p^n_{i+2}-p^n_{i+1} }{3} - \frac{ p^n_{i+3}-p^n_{i} }{3}\right) -\frac{b_{i-\frac{3}{2}}}{3} \left( 27 \frac{ p^n_{i-1}-p^n_{i-2}}{3} - \frac{ p^n_{i}-p^n_{i-3} }{3}\right) \right]
$$

$$
=\frac{3}{64 \Delta x ^2} \left[ b_{i+\frac{1}{2}} \left( 27 ( p^n_{i+1}-p^n_i ) - p^n_{i+2}+p^n_{i-1} \right) -b_{i-\frac{1}{2}} \left( 27 (p^n_{i}-p^n_{i-1}) - p^n_{i+1}+p^n_{i-2} \right) \right]
$$
$$
-\frac{1}{576\Delta x ^2} \left[ b_{i+\frac{3}{2}} \left( 27 (p^n_{i+2}-p^n_{i+1} ) - p^n_{i+3}+p^n_{i} \right) -b_{i-\frac{3}{2}} \left( 27 (p^n_{i-1}-p^n_{i-2}) -p^n_{i}+p^n_{i-3} \right) \right]
$$

Then opening $b$ using linear interpolation

$$
=\frac{3}{128 \Delta x ^2} \left[ (b_{i+1}+b_i) \left( 27 ( p^n_{i+1}-p^n_i ) - p^n_{i+2}+p^n_{i-1} \right) -(b_{i-1}+b_i) \left( 27 (p^n_{i}-p^n_{i-1}) - p^n_{i+1}+p^n_{i-2} \right) \right]
$$
$$
-\frac{1}{1152\Delta x ^2} \left[ (b_{i+2}+b_{i+1}) \left( 27 (p^n_{i+2}-p^n_{i+1} ) - p^n_{i+3}+p^n_{i} \right) -(b_{i-2}+b_{i-1}) \left( 27 (p^n_{i-1}-p^n_{i-2}) -p^n_{i}+p^n_{i-3} \right) \right]
$$

 Di Bartolo, L. A new family of finite-difference schemes to solve the heterogeneous acoustic wave equation - Geophysics 2012

Comments