Formulations (Under Development)

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

(1)\dpd{\rho}{t} + \sum_{i=1}^3 \dpd{\rho v_i}{x_i} = 0

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)\bvec{u} \defeq \left(\begin{array}{c} \rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \end{array}\right)

and flux functions

(6)\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)

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)\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)

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)\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

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

(11)\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), \quad \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)

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)\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)

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)\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)

Aided by writing (5) as

\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)

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

(14)\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)

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)\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), \quad \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)

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

(18)\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)

(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)\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)

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)\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)

The left eigenvector matrix is

(24)\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)

Riemann Invariants

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

(25)\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}

(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)\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)

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)\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)