# Formulations (Under Development)¶

The governing equations of the hydro-acoustic wave include the continuity equation

(1)\begin{align}\begin{aligned}\newcommand{\dpd}[3][]{\mathinner{ \dfrac{\partial{^{#1}}#2}{\partial{#3^{#1}}} }}\\\dpd{\rho}{t} + \sum_{i=1}^3 \dpd{\rho v_i}{x_i} = 0\end{aligned}\end{align}

and the momentum equations

(2)$\dpd{\rho v_j}{t} + \sum_{i=1}^3 \dpd{(\rho v_iv_j + \delta_{ij}p)}{x_i} = \dpd{}{x_j}\left(\lambda\sum_{k=1}^3\dpd{v_k}{x_k}\right) + \sum_{i=1}^3 \dpd{}{x_i} \left[\mu \left( \dpd{v_i}{x_j}+\dpd{v_j}{x_i} \right)\right], \quad j = 1, 2, 3$

where $$\rho$$ is the density, $$v_1, v_2,$$ and $$v_3$$ the Cartesian component of the velocity, $$p$$ the pressure, $$\delta_{ij}, i, j = 1, 2, 3$$ the Kronecker delta, $$\lambda$$ the second viscosity coefficien, $$\mu$$ the dynamic viscosity coefficient, $$t$$ the time, and $$x_1, x_2$$, and $$x_3$$ the Cartesian coordinate axes. Newtonian fluid is assumed.

The above four equations in (1) and (2) have five independent variables $$\rho, p, v_1, v_2$$, and $$v_3$$, and hence are not closed without a constitutive relation. In the bulk package, the constitutive relation (or the equation of state) of choice is

$K = \rho\dpd{p}{\rho}$

where $$K$$ is a constant and the bulk modulus. We chose to use the density $$\rho$$ as the independent variable, and integrate the equation of state to be

(3)$p = p_0 + K \ln\frac{\rho}{\rho_0}$

where $$p_0$$ and $$\rho_0$$ are constants. Substituting (3) into (2) gives

(4)$\dpd{\rho v_j}{t} + \sum_{i=1}^3\dpd{}{x_i} \left[\rho v_iv_j + \delta_{ij}\left(p_0 + K\ln\frac{\rho}{\rho_0}\right) \right] = \sum_{i=1}^3 \dpd{}{x_i} \left[\delta_{ij} \lambda \sum_{k=1}^3 \dpd{v_k}{x_k} + \mu \left( \dpd{v_i}{x_j}+\dpd{v_j}{x_i} \right)\right], \quad j = 1, 2, 3$

## Jacobian Matrices¶

We proceed to analyze the advective part of the governing equations (1) and (4). Define the conservation variables


and flux functions

(6)$\begin{split}\bvec{f}^{(1)} \defeq \left(\begin{array}{c} \rho v_1 \\ \rho v_1^2 + K\ln\frac{\rho}{\rho_0} + p_0 \\ \rho v_1v_2 \\ \rho v_1v_3 \end{array}\right), \quad \bvec{f}^{(2)} \defeq \left(\begin{array}{c} \rho v_2 \\ \rho v_1v_2 \\ \rho v_2^2 + K\ln\frac{\rho}{\rho_0} + p_0 \\ \rho v_2v_3 \end{array}\right), \quad \bvec{f}^{(3)} \defeq \left(\begin{array}{c} \rho v_3 \\ \rho v_1v_3 \\ \rho v_2v_3 \\ \rho v_3^2 + K\ln\frac{\rho}{\rho_0} + p_0 \end{array}\right)\end{split}$

Aided by the definition of conservation variables in (5), the flux functions defined in (6) can be rewritten with $$u_1, u_2, u_3$$, and $$u_4$$

(7)$\begin{split}\bvec{f}^{(1)} = \left(\begin{array}{c} u_2 \\ \frac{u_2^2}{u_1} + K\ln\frac{u_1}{\rho_0} + p_0 \\ \frac{u_2u_3}{u_1} \\ \frac{u_2u_4}{u_1} \end{array}\right), \quad \bvec{f}^{(2)} = \left(\begin{array}{c} u_3 \\ \frac{u_2u_3}{u_1} \\ \frac{u_3^2}{u_1} + K\ln\frac{u_1}{\rho_0} + p_0 \\ \frac{u_3u_4}{u_1} \end{array}\right), \quad \bvec{f}^{(3)} = \left(\begin{array}{c} u_4 \\ \frac{u_2u_4}{u_1} \\ \frac{u_3u_4}{u_1} \\ \frac{u_4^2}{u_1} + K\ln\frac{u_1}{\rho_0} + p_0 \end{array}\right)\end{split}$

By using (5) and (7), the left-hand-side of the governing equations can be cast into the conservative form

(8)$\dpd{\bvec{u}}{t} + \sum_{i=1}^3\dpd{\bvec{f}^{(i)}}{x_i} = 0$

Aided by using the chain rule, (8) can be rewritten in the quasi-linear form

(9)$\dpd{\bvec{u}}{t} + \sum_{i=1}^3\mathrm{A}^{(i)}\dpd{\bvec{u}}{x_i} = 0$

where the Jacobian matrices $$\mathrm{A}^{(1)}, \mathrm{A}^{(2)}$$, and $$\mathrm{A}^{(3)}$$ are defined as

(10)\begin{align}\begin{aligned}\newcommand{\pd}[3][]{ \tfrac{\partial{^{#1}}#2}{\partial{#3^{#1}}} }\\\begin{split}\mathrm{A}^{(i)} \defeq \left(\begin{array}{cccc} \pd{f_1^{(i)}}{u_1} & \pd{f_1^{(i)}}{u_2} & \pd{f_1^{(i)}}{u_3} & \pd{f_1^{(i)}}{u_4} \\ \pd{f_2^{(i)}}{u_1} & \pd{f_2^{(i)}}{u_2} & \pd{f_2^{(i)}}{u_3} & \pd{f_2^{(i)}}{u_4} \\ \pd{f_3^{(i)}}{u_1} & \pd{f_3^{(i)}}{u_2} & \pd{f_3^{(i)}}{u_3} & \pd{f_3^{(i)}}{u_4} \\ \pd{f_4^{(i)}}{u_1} & \pd{f_4^{(i)}}{u_2} & \pd{f_4^{(i)}}{u_3} & \pd{f_4^{(i)}}{u_4} \end{array}\right), \quad i = 1, 2, 3\end{split}\end{aligned}\end{align}

Aided by using (7), the Jacobian matrices defined in (10) can be written out as

(11)$\begin{split}\mathrm{A}^{(1)} &= \left(\begin{array}{cccc} 0 & 1 & 0 & 0 \\ -\frac{u_2^2}{u_1^2} + \frac{K}{u_1} & 2\frac{u_2}{u_1} & 0 & 0 \\ -\frac{u_2u_3}{u_1^2} & \frac{u_3}{u_1} & \frac{u_2}{u_1} & 0 \\ -\frac{u_2u_4}{u_1^2} & \frac{u_4}{u_1} & 0 & \frac{u_2}{u_1} \end{array}\right), \quad \mathrm{A}^{(2)} = \left(\begin{array}{cccc} 0 & 0 & 1 & 0 \\ -\frac{u_2u_3}{u_1^2} & \frac{u_3}{u_1} & \frac{u_2}{u_1} & 0 \\ -\frac{u_3^2}{u_1^2} + \frac{K}{u_1} & 0 & 2\frac{u_3}{u_1} & 0 \\ -\frac{u_3u_4}{u_1^2} & 0 & \frac{u_4}{u_1} & \frac{u_3}{u_1} \end{array}\right), \\ \mathrm{A}^{(3)} &= \left(\begin{array}{cccc} 0 & 0 & 0 & 1 \\ -\frac{u_2u_4}{u_1^2} & \frac{u_4}{u_1} & 0 & \frac{u_2}{u_1} \\ -\frac{u_3u_4}{u_1^2} & 0 & \frac{u_4}{u_1} & \frac{u_3}{u_1} \\ -\frac{u_4^2}{u_1^2} & 0 & 0 & 2\frac{u_4}{u_1} \end{array}\right)\end{split}$

## Hyperbolicity¶

Hyperbolicity is a prerequisite for us to apply the space-time CESE method to a system of first-order PDEs. For the governing equations, (1) and (4), to be hyperbolic, the linear combination of the three Jacobian matrices of their quasi-linear for must be diagonalizable. The spectrum of the linear combination must be all real, too [Warming75], [Chen12].

To facilitate the analysis, we chose to use the non-conservative version of the governing equations. Define the non-conservative variables

(12)$\begin{split}\tilde{\bvec{u}} \defeq \left(\begin{array}{c} \rho \\ v_1 \\ v_2 \\ v_3 \end{array}\right) = \left(\begin{array}{c} u_1 \\ \frac{u_2}{u_1} \\ \frac{u_3}{u_1} \\ \frac{u_4}{u_1} \end{array}\right)\end{split}$

Aided by using (12) and (5), the transformation between the conservative variables and the non-conservative variables can be done with the transformation Jacobian defined as

(13)$\begin{split}\mathrm{P} \defeq \dpd{\tilde{\bvec{u}}}{\bvec{u}} = \left(\begin{array}{cccc} \pd{\tilde{u}_1}{u_1} & \pd{\tilde{u}_1}{u_2} & \pd{\tilde{u}_1}{u_3} & \pd{\tilde{u}_1}{u_4} \\ \pd{\tilde{u}_2}{u_1} & \pd{\tilde{u}_2}{u_2} & \pd{\tilde{u}_2}{u_3} & \pd{\tilde{u}_2}{u_4} \\ \pd{\tilde{u}_3}{u_1} & \pd{\tilde{u}_3}{u_2} & \pd{\tilde{u}_3}{u_3} & \pd{\tilde{u}_3}{u_4} \\ \pd{\tilde{u}_4}{u_1} & \pd{\tilde{u}_4}{u_2} & \pd{\tilde{u}_4}{u_3} & \pd{\tilde{u}_4}{u_4} \end{array}\right) = \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ -\frac{u_2}{u_1^2} & \frac{1}{u_1} & 0 & 0 \\ -\frac{u_3}{u_1^2} & 0 & \frac{1}{u_1} & 0 \\ -\frac{u_4}{u_1^2} & 0 & 0 & \frac{1}{u_1} \end{array}\right) = \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ -\frac{v_1}{\rho} & \frac{1}{\rho} & 0 & 0 \\ -\frac{v_2}{\rho} & 0 & \frac{1}{\rho} & 0 \\ -\frac{v_3}{\rho} & 0 & 0 & \frac{1}{\rho} \end{array}\right)\end{split}$

Aided by writing (5) as

$\begin{split}\bvec{u} = \left(\begin{array}{c} \tilde{u}_1 \\ \tilde{u}_1\tilde{u}_2 \\ \tilde{u}_1\tilde{u}_3 \\ \tilde{u}_1\tilde{u}_4 \end{array}\right)\end{split}$

the inverse matrix of $$\mathrm{P}$$ can be obtained

(14)$\begin{split}\mathrm{P}^{-1} = \dpd{\bvec{u}}{\tilde{\bvec{u}}} = \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ \tilde{u}_2 & \tilde{u}_1 & 0 & 0 \\ \tilde{u}_3 & 0 & \tilde{u}_1 & 0 \\ \tilde{u}_4 & 0 & 0 & \tilde{u}_1 \end{array}\right) = \left(\begin{array}{cccc} 1 & 0 & 0 & 0 \\ v_1 & \rho & 0 & 0 \\ v_2 & 0 & \rho & 0 \\ v_3 & 0 & 0 & \rho \end{array}\right)\end{split}$

and $$\mathrm{P}^{-1}\mathrm{P} = \mathrm{P}\mathrm{P}^{-1} = \mathrm{I}_{4\times4}$$.

The transformation matrix $$\mathrm{P}$$ can be used to cast the conservative equations, (9), into non-conservative ones. Pre-multiplying $$\pd{\tilde{\bvec{u}}}{\bvec{u}}$$ to (9) gives

(15)$\dpd{\tilde{\bvec{u}}}{t} + \sum_{i=1}^3\tilde{\mathrm{A}}^{(i)}\dpd{\tilde{\bvec{u}}}{x_i} = 0$

where

(16)$\tilde{\mathrm{A}}^{(1)} \defeq \mathrm{P}\mathrm{A}^{(1)}\mathrm{P}^{-1}, \quad \tilde{\mathrm{A}}^{(2)} \defeq \mathrm{P}\mathrm{A}^{(2)}\mathrm{P}^{-1}, \quad \tilde{\mathrm{A}}^{(3)} \defeq \mathrm{P}\mathrm{A}^{(3)}\mathrm{P}^{-1}$

To help obtaining the expression of $$\tilde{\mathrm{A}}^{(1)}, \tilde{\mathrm{A}}^{(2)}$$, and $$\tilde{\mathrm{A}}^{(3)}$$, we substitute (5) into (11) and get

(17)$\begin{split}\mathrm{A}^{(1)} &= \left(\begin{array}{cccc} 0 & 1 & 0 & 0 \\ -v_1^2 + \frac{K}{\rho} & 2v_1 & 0 & 0 \\ -v_1v_2 & v_2 & v_1 & 0 \\ -v_1v_3 & v_3 & 0 & v_1 \end{array}\right), \quad \mathrm{A}^{(2)} = \left(\begin{array}{cccc} 0 & 0 & 1 & 0 \\ -v_1v_2 & v_2 & v_1 & 0 \\ -v_2^2 + \frac{K}{\rho} & 0 & 2v_2 & 0 \\ -v_2v_3 & 0 & v_3 & v_2 \end{array}\right), \\ \mathrm{A}^{(3)} &= \left(\begin{array}{cccc} 0 & 0 & 0 & 1 \\ -v_1v_3 & v_3 & 0 & v_1 \\ -v_2v_3 & 0 & v_3 & v_2 \\ -v_3^2 + \frac{K}{\rho} & 0 & 0 & 2v_3 \end{array}\right)\end{split}$

The Jacobian matrices in (16) can be spelled out with the expressions in (13), (14), and (17)

(18)$\begin{split}\tilde{\mathrm{A}}^{(1)} = \left(\begin{array}{cccc} v_1 & \rho & 0 & 0 \\ \frac{K}{\rho^2} & v_1 & 0 & 0 \\ 0 & 0 & v_1 & 0 \\ 0 & 0 & 0 & v_1 \end{array}\right), \quad \tilde{\mathrm{A}}^{(2)} = \left(\begin{array}{cccc} v_2 & 0 & \rho & 0 \\ 0 & v_2 & 0 & 0 \\ \frac{K}{\rho^2} & 0 & v_2 & 0 \\ 0 & 0 & 0 & v_2 \end{array}\right), \quad \tilde{\mathrm{A}}^{(3)} = \left(\begin{array}{cccc} v_3 & 0 & 0 & \rho \\ 0 & v_3 & 0 & 0 \\ 0 & 0 & v_3 & 0 \\ \frac{K}{\rho^2} & 0 & 0 & v_3 \end{array}\right)\end{split}$

(15) is hyperbolic where the linear combination of its Jacobian matrices $$\tilde{\mathrm{A}}^{(1)}$$, $$\tilde{\mathrm{A}}^{(2)}$$, and $$\tilde{\mathrm{A}}^{(3)}$$

(19)$\begin{split}\tilde{\mathrm{R}} \defeq \sum_{i=1}^3 k_i\tilde{\mathrm{A}}^{(i)} = \left(\begin{array}{cccc} \sum_{i=1}^3 k_iv_i & k_1\rho & k_2\rho & k_3\rho \\ \frac{k_1K}{\rho^2} & \sum_{i=1}^3 k_iv_i & 0 & 0 \\ \frac{k_2K}{\rho^2} & 0 & \sum_{i=1}^3 k_iv_i & 0 \\ \frac{k_3K}{\rho^2} & 0 & 0 & \sum_{i=1}^3 k_iv_i \end{array}\right)\end{split}$

where $$k_1, k_2$$, and $$k_3$$ are real and bounded.

The linearly combined Jacobian matrix can be used to rewrite the three-dimensional governing equation, (15), into one-dimensional space

(20)$\dpd{\tilde{\bvec{u}}}{t} + \tilde{\mathrm{R}}\dpd{\tilde{\bvec{u}}}{y} = 0$

where

(21)$y \defeq \sum_{i=1}^3 k_ix_i$

and aided by (19) and the chain rule

$\sum_{i=1}^3 \tilde{\mathrm{A}}^{(i)} \dpd{\tilde{\bvec{u}}}{x_i} = \sum_{i=1}^3 \tilde{\mathrm{A}}^{(i)} \dpd{\tilde{\bvec{u}}}{y} \dpd{y}{x_i} = \sum_{i=1}^3 k_{i} \tilde{\mathrm{A}}^{(i)} \dpd{\tilde{\bvec{u}}}{y} = \tilde{\mathrm{R}}\dpd{\tilde{\bvec{u}}}{y}$

The eigenvalues of $$\tilde{\mathrm{R}}$$ can be found by solving the polynomial equation $$\det(\tilde{\mathrm{R}} - \lambda\mathrm{I}_{4\times4}) = 0$$ for $$\lambda$$, and are

(22)$\lambda_{1, 2, 3, 4} = \phi, \phi, \phi+\sqrt{\frac{K\psi}{\rho}}, \phi-\sqrt{\frac{K\psi}{\rho}}$

where $$\phi \defeq \sum_{i=1}^3 k_iv_i$$, and $$\psi \defeq \sum_{i=1}^3 k_i^2$$. The corresponding eigenvector matrix is

(23)$\begin{split}\mathrm{T} = \left(\begin{array}{cccc} 0 & 0 & \sqrt{\frac{\rho^3\psi}{K}} & \sqrt{\frac{\rho^3\psi}{K}} \\ k_3 & 0 & k_1 & -k_1 \\ 0 & k_3 & k_2 & -k_2 \\ -k_1 & -k_2 & k_3 & -k_3 \end{array}\right)\end{split}$

The left eigenvector matrix is

(24)$\begin{split}\mathrm{T}^{-1} = \left(\begin{array}{cccc} 0 & -\frac{k_1^2-\psi}{k_3\psi} & -\frac{k_1k_2}{k_3\psi} & -\frac{k_1}{\psi} \\ 0 & -\frac{k_1k_2}{k_3\psi} & -\frac{k_2^2-\psi}{k_3\psi} & -\frac{k_2}{\psi} \\ \frac{1}{2\sqrt{\frac{\rho^3\psi}{K}}} & \frac{k_1}{2\psi} & \frac{k_2}{2\psi} & \frac{k_3}{2\psi} \\ \frac{1}{2\sqrt{\frac{\rho^3\psi}{K}}} & -\frac{k_1}{2\psi} & -\frac{k_2}{2\psi} & -\frac{k_3}{2\psi} \end{array}\right)\end{split}$

## Riemann Invariants¶

Aided by (23) and (24), $$\tilde{\mathrm{R}}$$ can be diagonalized as

(25)$\begin{split}\mathrm{\Lambda} \defeq \left(\begin{array}{cccc} \lambda_1 & 0 & 0 & 0 \\ 0 & \lambda_2 & 0 & 0 \\ 0 & 0 & \lambda_3 & 0 \\ 0 & 0 & 0 & \lambda_4 \end{array}\right) = \left(\begin{array}{cccc} \phi & 0 & 0 & 0 \\ 0 & \phi & 0 & 0 \\ 0 & 0 & \phi + \sqrt{\frac{K\psi}{\rho}} & 0 \\ 0 & 0 & 0 & \phi - \sqrt{\frac{K\psi}{\rho}} \end{array}\right) = \mathrm{T}^{-1}\tilde{\mathrm{R}}\mathrm{T}\end{split}$

(25) defines the eigenvalue matrix of $$\tilde{\mathrm{R}}$$. Aach element in the diagonal of the eigenvalue matrix is the eigenvalue listed in (22). Pre-multiplying (20) with $$\mathrm{T}^{-1}$$ gives

$\mathrm{T}^{-1}\dpd{\tilde{\bvec{u}}}{t} + \mathrm{\Lambda}\mathrm{T}^{-1}\dpd{\tilde{\bvec{u}}}{y} = 0$

Define the characteristic variables

(26)$\begin{split}\hat{\bvec{u}} \defeq \left(\begin{array}{c} -\frac{k_1^2-\psi}{k_3\psi}v_1 - \frac{k_1k_2}{k_3\psi}v_2 - \frac{k_1}{\psi}v_3 \\ -\frac{k_1k_2}{k_3\psi}v_1 - \frac{k_2^2-\psi}{k_3\psi}v_2 - \frac{k_2}{\psi}v_3 \\ -\sqrt{\frac{K}{\rho\psi}} + \frac{k_1}{2\psi}v_1 + \frac{k_2}{2\psi}v_2 + \frac{k_3}{2\psi}v_3 \\ -\sqrt{\frac{K}{\rho\psi}} - \frac{k_1}{2\psi}v_1 - \frac{k_2}{2\psi}v_2 - \frac{k_3}{2\psi}v_3 \end{array}\right)\end{split}$

such that

$\mathrm{T}^{-1} = \dpd{\hat{\bvec{u}}}{\tilde{\bvec{u}}}$

Then aided by the chain rule, we can write

(27)$\dpd{\hat{\bvec{u}}}{t} + \mathrm{\Lambda}\dpd{\hat{\bvec{u}}}{y} = 0$

The components of $$\hat{\bvec{u}}$$ defined in (26) are the Riemann invariants.

## Diffusion Term Treatment¶

The momentum equation (4) contains the diffusion term

$\sum_{i=1}^3 \dpd{}{x_i} \left[\delta_{ij} \lambda \sum_{k=1}^3 \dpd{v_k}{x_k} + \mu \left( \dpd{v_i}{x_j}+\dpd{v_j}{x_i} \right)\right], \quad j = 1, 2, 3$

Define

(28)$\begin{split}\bvec{g}^{(1)} \defeq \left(\begin{array}{c} 0 \\ \lambda\sum_{k=1}^3\dpd{v_k}{x_k} + 2\mu\dpd{v_1}{x_1} \\ \mu(\dpd{v_1}{x_2} + \dpd{v_2}{x_1}) \\ \mu(\dpd{v_1}{x_3} + \dpd{v_3}{x_1}) \end{array}\right), \quad \bvec{g}^{(2)} \defeq \left(\begin{array}{c} 0 \\ \mu(\dpd{v_2}{x_1} + \dpd{v_1}{x_2}) \\ \lambda\sum_{k=1}^3\dpd{v_k}{x_k} + 2\mu\dpd{v_2}{x_2} \\ \mu(\dpd{v_2}{x_3} + \dpd{v_3}{x_2}) \end{array}\right), \quad \bvec{g}^{(3)} \defeq \left(\begin{array}{c} 0 \\ \mu(\dpd{v_3}{x_1} + \dpd{v_1}{x_3}) \\ \mu(\dpd{v_3}{x_2} + \dpd{v_2}{x_3}) \\ \lambda\sum_{k=1}^3\dpd{v_k}{x_k} + 2\mu\dpd{v_3}{x_3} \end{array}\right)\end{split}$