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

(5)\[ \begin{align}\begin{aligned}\newcommand{\bvec}[1]{\mathbf{#1}} \newcommand{\defeq}{\buildrel{\text{def}}\over{=}}\\\begin{split}\bvec{u} \defeq \left(\begin{array}{c} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \end{array}\right)\end{split}\end{aligned}\end{align} \]

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}\]