Acoustic Wave Equation

Derivation

The acoustic wave equation in a heterogeneous medium with variable density \(\rho(\mathbf{x})\) and bulk modulus \(K(\mathbf{x})\) is derived from the linearized equations of momentum conservation and mass conservation (continuity equation).

Momentum Conservation (Equation of Motion):

(1)\[\rho \frac{\partial \mathbf{v}}{\partial t} = -\nabla p + \mathbf{f},\]

where:

  • \(\mathbf{v}\) is the particle velocity vector.

  • \(p\) is the acoustic pressure.

  • \(\mathbf{f}\) is the body force source term (force per unit volume).

  • \(\rho\) is the mass density.

Mass Conservation (Continuity Equation):

(2)\[\frac{1}{K} \frac{\partial p}{\partial t} = -\nabla \cdot \mathbf{v} + s,\]

where:

  • \(K = \rho c^2\) is the adiabatic bulk modulus, with \(c\) being the acoustic wave speed.

  • \(s\) is a source term (volume injection rate).

Implementation

The system is discretized using a finite-difference time-domain (FDTD) method on a staggered grid.

Staggered Grid Layout

For a 2D domain in the \((y, x)\) plane (using \(y\) for the first dimension and \(x\) for the second to match the array indexing convention):

  • Pressure \(p\) and bulk modulus \(K\) are defined at grid nodes \((y, x)\).

  • Vertical velocity \(v_y\) and density \(\rho_y\) are defined at half-integer steps in \(y\): \((y+1/2, x)\).

  • Horizontal velocity \(v_x\) and density \(\rho_x\) are defined at half-integer steps in \(x\): \((y, x+1/2)\).

Time Integration

Time stepping is performed using an explicit second-order leapfrog scheme. Pressure and particle velocities are staggered in time by half a time step \(h_t/2\).

  1. Update Particle Velocities:

    (3)\[\mathbf{v}^{t+1/2} = \mathbf{v}^{t-1/2} - h_t \frac{1}{\rho} \nabla p^t + h_t \frac{1}{\rho} \mathbf{f}^t.\]

  2. Update Pressure:

    (4)\[p^{t+1} = p^t - h_t K (\nabla \cdot \mathbf{v}^{t+1/2} - s^{t+1/2}).\]

Absorbing Boundaries (CPML)

To mitigate artificial reflections from the domain boundaries, a Convolutional Perfectly Matched Layer (CPML) is employed. This method introduces memory variables to implement the coordinate stretching.

For the velocity update, we replace the gradient \(\nabla p\) with \(\nabla p + \mathbf{\psi}\). The memory variable \(\mathbf{\psi}\) is updated as:

\[\psi_i^t = a_i \psi_i^{t-1} + b_i \partial_i p^t\]

where \(i \in \{y, x\}\).

For the pressure update, we replace the divergence \(\nabla \cdot \mathbf{v}\) with \(\nabla \cdot \mathbf{v} + \sum_i \phi_i\). The memory variable \(\mathbf{\phi}\) is updated as:

\[\phi_i^{t+1/2} = a_i \phi_i^{t-1/2} + b_i \partial_i v_i^{t+1/2}\]

The coefficients \(a\) and \(b\) are determined by the PML damping profile.

Matrix Form

We express the update steps in matrix form, which facilitates the derivation of the adjoint operators. We use \(B\) for buoyancy (\(1/\rho\)) and \(h_t\) for the time step \(h_t\). The notation \(y'\) and \(x'\) denotes locations shifted by half a grid cell (e.g., \(y+1/2\)).

1. Velocity Update

This step updates the particle velocities and the gradient memory variables.

\[\begin{split}\begin{pmatrix} v_y^{t+\frac{1}{2}} \\ v_x^{t+\frac{1}{2}} \\ \psi_y^t \\ \psi_x^t \end{pmatrix} = \begin{pmatrix} -h_t B^{y'x} (1+b^{y'}) \partial_{y'} & 1 & 0 & -h_t B^{y'x} a^{y'} & 0 & h_t B^{y'x} & 0 \\ -h_t B^{yx'} (1+b^{x'}) \partial_{x'} & 0 & 1 & 0 & -h_t B^{yx'} a^{x'} & 0 & h_t B^{yx'} \\ b^{y'} \partial_{y'} & 0 & 0 & a^{y'} & 0 & 0 & 0 \\ b^{x'} \partial_{x'} & 0 & 0 & 0 & a^{x'} & 0 & 0 \end{pmatrix} \begin{pmatrix} p^t \\ v_y^{t-\frac{1}{2}} \\ v_x^{t-\frac{1}{2}} \\ \psi_y^{t-1} \\ \psi_x^{t-1} \\ f_y^t \\ f_x^t \end{pmatrix} \label{acoustic_matrix_v}\end{split}\]

2. Pressure Update

This step updates the pressure and the divergence memory variables.

\[\begin{split}\begin{pmatrix} p^{t+1} \\ \phi_y^{t+\frac{1}{2}} \\ \phi_x^{t+\frac{1}{2}} \end{pmatrix} = \begin{pmatrix} -h_t K^{yx} (1+b^y) \partial_y & -h_t K^{yx} (1+b^x) \partial_x & 1 & -h_t K^{yx} a^y & -h_t K^{yx} a^x & h_t K^{yx} \\ b^y \partial_y & 0 & 0 & a^y & 0 & 0 \\ 0 & b^x \partial_x & 0 & 0 & a^x & 0 \end{pmatrix} \begin{pmatrix} v_y^{t+\frac{1}{2}} \\ v_x^{t+\frac{1}{2}} \\ p^t \\ \phi_y^{t-\frac{1}{2}} \\ \phi_x^{t-\frac{1}{2}} \\ s^{t+\frac{1}{2}} \end{pmatrix} \label{acoustic_matrix_p}\end{split}\]