For the temporal integration of the governing equations implicit Runge-Kutta schemes are employed. These schemes take the following generic form

"Eqn_IRK-Generic01"

\begin{equation}\label{Eqn_IRK-Generic01} \left[\left(\frac{1}{\Delta \tau} + \frac{1}{\Delta t}\right)\mathbf{I} - a_{ii}\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\right]\Delta \hat{\mathbf{Q}} = -\left[\frac{\hat{\mathbf{Q}}^m - \hat{\mathbf{Q}}^n}{\Delta t} + \sum_{j=1}^{i-1}a_{ij}\mathbf{R}^j + a_{ii}\mathbf{R}^m\right] \end{equation}

where \(\mathbf{R}^j = \mathbf{R}(\hat{\mathbf{Q}}^j)\). As \(\Delta \hat{\mathbf{Q}} \rightarrow 0\), \(\hat{\mathbf{Q}}^{m+1} = \hat{\mathbf{Q}}^{m} = \hat{\mathbf{Q}}^{i}\), where \(\hat{\mathbf{Q}}^{i}\) is the solution of the \(i\)-th stage. For non-deforming meshes Eqn. Eqn_IRK-Generic01 reduces to the following form:

\begin{equation} \left[\left(\frac{1}{J\Delta \tau} + \frac{1}{J\Delta t}\right)\mathbf{I} - a_{ii}\frac{\partial \mathbf{R}}{\partial \mathbf{Q}}\right]\Delta \mathbf{Q} = -\left[\frac{\mathbf{Q}^m - \mathbf{Q}^n}{J\Delta t} + \sum_{j=1}^{i-1}a_{ij} \mathbf{R}^j + a_{ii}\mathbf{R}^m\right] \end{equation}

The flux Jacobian \(\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\) is defined, using Eqn. Eqn_3D-NS-FV, as follows:

\begin{eqnarray} \left(\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\right)_{i,j,k} = &-& \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}} - \frac{\partial}{\partial\xi}\frac{\mathbf{\partial \hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\right)_{i+1/2,j,k} + \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}} - \frac{\partial}{\partial\xi}\frac{\mathbf{\partial \hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\right)_{i-1/2,j,k} \\ &-& \left(\frac{\partial \mathbf{\hat{G}}}{\partial \hat{\mathbf{Q}}} - \frac{\partial}{\partial\eta}\frac{\mathbf{\partial \hat{G}}_v}{\partial \hat{\mathbf{Q}}_{\eta}}\right)_{i,j+1/2,k} + \left(\frac{\partial \mathbf{\hat{G}}}{\partial \hat{\mathbf{Q}}} - \frac{\partial}{\partial\eta}\frac{\mathbf{\partial \hat{G}}_v}{\partial \hat{\mathbf{Q}}_{\eta}}\right)_{i,j-1/2,k} \\ &-& \left(\frac{\partial \mathbf{\hat{H}}}{\partial \hat{\mathbf{Q}}} - \frac{\partial}{\partial\zeta}\frac{\mathbf{\partial \hat{H}}_v}{\partial \hat{\mathbf{Q}}_{\zeta}}\right)_{i,j,k+1/2} + \left(\frac{\partial \mathbf{\hat{H}}}{\partial \hat{\mathbf{Q}}} - \frac{\partial}{\partial\zeta}\frac{\mathbf{\partial \hat{H}}_v}{\partial \hat{\mathbf{Q}}_{\zeta}}\right)_{i,j,k-1/2} \\ &+& \left(\frac{\partial \mathbf{\hat{S}}}{\partial \hat{\mathbf{Q}}}\right)_{i,j,k} \end{eqnarray}

Here \(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\), \(\frac{\partial \mathbf{\hat{G}}}{\partial \hat{\mathbf{Q}}}\) and \(\frac{\partial \mathbf{\hat{H}}}{\partial \hat{\mathbf{Q}}}\) are the inviscid flux Jacobians and \(\frac{\mathbf{\partial \hat{F}}_v}{\partial \hat{\mathbf{Q}}}\), \(\frac{\mathbf{\partial \hat{G}}_v}{\partial \hat{\mathbf{Q}}}\) and \(\frac{\mathbf{\partial \hat{H}}_v}{\partial \hat{\mathbf{Q}}}\) are the viscous flux Jacobians. These are computed analytically in TRACE. The expressions for the various flux Jacobians are given in the following sections.

Inviscid Flux Jacobians

Considering the inviscid flux Jacobian \(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\) at the left face of the \(i\)-th cell, e.g. at \(i+1/2,j,k\), we have

"Eqn_InviscidFluxJacobian"

\begin{eqnarray} \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i+1/2,j,k} &=& \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\hat{\mathbf{Q}}^{m+1} - \frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\hat{\mathbf{Q}}^m\right)_{i+1/2,j,k} \\ &=& \mathbf{\hat{F}}^{m+1}_{i+1/2,j,k} - \mathbf{\hat{F}}^m_{i+1/2,j,k} \end{eqnarray}

Using upwind approximations for the fluxes, we have

"Eqn_RoeFlux"

\begin{equation}\label{Eqn_RoeFlux} \mathbf{\hat{F}}_{i+1/2,j,k} = \frac{1}{2}\left[(\mathbf{F}_L + \mathbf{F}_R) + |\tilde{\mathbf{A}}|(\mathbf{U}_L - \mathbf{U}_R)\right] \end{equation}

where \(\mathbf{F}_L = \tilde{\mathbf{A}}_L \mathbf{U}_L\), \(\mathbf{F}_R = \tilde{\mathbf{A}}_R \mathbf{U}_R\) and \(|\tilde{\mathbf{A}}| = |\tilde{\mathbf{A}}(\mathbf{U}_L,\mathbf{U}_R)|\). For the implicit operator first-order upwind approximations are used, and therefore \(\mathbf{U}_L = \mathbf{U}_i\) and \(\mathbf{U}_R = \mathbf{U}_{i+1}\). Eqn. Eqn_RoeFlux then becomes

\begin{equation} \mathbf{F}_{i+1/2,j} = \frac{1}{2}\left[(\mathbf{A}_i\mathbf{U}_i + \mathbf{A}_{i+1}\mathbf{U}_{i+1}) + |\tilde{\mathbf{A}}|_{i+1/2}(\mathbf{U}_i - \mathbf{U}_{i+1})\right] \end{equation}

where \(|\tilde{\mathbf{A}}|_{i+1/2}\) denotes that the Jacobian is computed using the so-called Roe-averaged variables, which are evaluated as a combination of \(\mathbf{U}_i\) and \(\mathbf{U}_{i+1}\). Substituting these expressions into Eqn. Eqn_InviscidFluxJacobian the following expression is obtained for the inviscid flux Jacobian

\begin{eqnarray} \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i+1/2,j,k} &=& \frac{1}{2}\left[(\mathbf{A}_i\hat{\mathbf{Q}}_i + \mathbf{A}_{i+1}\hat{\mathbf{Q}}_{i+1}) + |\tilde{\mathbf{A}}|_{i+1/2} (\hat{\mathbf{Q}}_i - \hat{\mathbf{Q}}_{i+1})\right]^{m+1} \\ &-& \frac{1}{2}\left[(\mathbf{A}_i\hat{\mathbf{Q}}_i + \mathbf{A}_{i+1}\hat{\mathbf{Q}}_{i+1}) + |\tilde{\mathbf{A}}|_{i+1/2} (\hat{\mathbf{Q}}_i - \hat{\mathbf{Q}}_{i+1})\right]^{m} \end{eqnarray}

With the approximation \(\mathbf{A}^{m+1} = \mathbf{A}^m\) the following simplified expression for the inviscid flux Jacobian is obtained

\begin{eqnarray} \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i+1/2,j,k} &=& \frac{1}{2}\left[(\mathbf{A}_i\Delta\hat{\mathbf{Q}}_i + \mathbf{A}_{i+1}\Delta\hat{\mathbf{Q}}_{i+1}) + |\tilde{\mathbf{A}}|_{i+1/2}(\Delta\hat{\mathbf{Q}}_{i} - \Delta\hat{\mathbf{Q}}_{i+1})\right] \\ &=& \mathbf{A}^{+}_{i+1/2,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} + \mathbf{A}^{-}_{i+1/2,j,k}\Delta\hat{\mathbf{Q}}_{i+1,j,k} \end{eqnarray}

where \(\mathbf{A}^+_{i+1/2,j,k} = \frac{1}{2}(\mathbf{A}_{i,j,k} + |\tilde{\mathbf{A}}|_{i+1/2,j,k})\) and \(\mathbf{A}^-_{i+1/2,j,k} = \frac{1}{2}(\mathbf{A}_{i+1,j,k} - |\tilde{\mathbf{A}}|_{i+1/2,j,k})\).

This is the discretization employed in the predictor-corrector solver in TRACE. The individual inviscid fluxes take the following form

"RoeFluxesPC"

\begin{eqnarray} \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i+1/2,j,k} &=& \mathbf{A}^+_{i+1/2,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} + \mathbf{A}^-_{i+1/2,j,k}\Delta\hat{\mathbf{Q}}_{i+1,j,k} \\ \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i-1/2,j,k} &=& \mathbf{A}^+_{i-1/2,j,k}\Delta\hat{\mathbf{Q}}_{i-1,j,k} + \mathbf{A}^-_{i-1/2,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} \\ \left(\frac{\partial \mathbf{\hat{G}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i,j+1/2,k} &=& \mathbf{B}^+_{i,j+1/2,k}\Delta\hat{\mathbf{Q}}_{i,j,k} + \mathbf{B}^-_{i,j+1/2,k}\Delta\hat{\mathbf{Q}}_{i,j+1,k} \\ \left(\frac{\partial \mathbf{\hat{G}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i,j-1/2,k} &=& \mathbf{B}^+_{i,j-1/2,k}\Delta\hat{\mathbf{Q}}_{i,j-1,k} + \mathbf{B}^-_{i,j-1/2,k}\Delta\hat{\mathbf{Q}}_{i,j,k} \\ \left(\frac{\partial \mathbf{\hat{H}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i,j,k+1/2} &=& \mathbf{C}^+_{i,j,k+1/2}\Delta\hat{\mathbf{Q}}_{i,j,k} + \mathbf{C}^-_{i,j,k+1/2}\Delta\hat{\mathbf{Q}}_{i,j,k+1} \\ \left(\frac{\partial \mathbf{\hat{H}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i,j,k-1/2} &=& \mathbf{C}^+_{i,j,k-1/2}\Delta\hat{\mathbf{Q}}_{i,j,k-1} + \mathbf{C}^-_{i,j,k-1/2}\Delta\hat{\mathbf{Q}}_{i,j,k} \end{eqnarray}

It is of course possible through the adoption of various assumptions to derive modified expressions for the inviscid flux Jacobians. For example, assuming \(\mathbf{A}_{i,j,k} = \mathbf{A}_{i-1/2,j,k} = \mathbf{A}_{i+1/2,j,k}\) we obtain

\begin{equation} \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i+1/2,j,k} = \mathbf{A}^{+}_{i+1/2,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} + \mathbf{A}^{-}_{i+1/2,j,k}\Delta\hat{\mathbf{Q}}_{i+1,j,k} \end{equation}

where \(\mathbf{A}^+_{i+1/2,j,k} = \frac{1}{2}(\mathbf{A}_{i+1/2,j,k} + |\tilde{\mathbf{A}}|_{i+1/2,j,k})\) and \(\mathbf{A}^-_{i+1/2,j,k} = \frac{1}{2}(\mathbf{A}_{i+1/2,j,k} - |\tilde{\mathbf{A}}|_{i+1/2,j,k})\). Due to the repetition of terms in the expressions for \(A^+\) and \(A^-\), this formulation is computationally less intensive than the previous one. It is, of course, also less accurate and possibly less stable. A third option proposed by Nürnberger for the discretization of the inviscid flux Jacobians is given below

"RoeFluxesSGS"

\begin{eqnarray} \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i+1/2,j,k} &=& \mathbf{A}_{i,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} + \mathbf{A}^-_{i+1/2,j,k}(\Delta\hat{\mathbf{Q}}_{i+1,j,k} - \Delta\hat{\mathbf{Q}}_{i,j,k}) \\ \left(\frac{\partial \mathbf{\hat{F}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i-1/2,j,k} &=& \mathbf{A}_{i,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} - \mathbf{A}^+_{i-1/2,j,k}(\Delta\hat{\mathbf{Q}}_{i,j,k} - \Delta\hat{\mathbf{Q}}_{i-1,j,k}) \\ \left(\frac{\partial \mathbf{\hat{G}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i,j+1/2,k} &=& \mathbf{B}_{i,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} + \mathbf{B}^-_{i,j+1/2,k}(\Delta\hat{\mathbf{Q}}_{i,j+1,k} - \Delta\hat{\mathbf{Q}}_{i,j,k}) \\ \left(\frac{\partial \mathbf{\hat{G}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i,j-1/2,k} &=& \mathbf{B}_{i,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} - \mathbf{B}^+_{i,j-1/2,k}(\Delta\hat{\mathbf{Q}}_{i,j,k} - \Delta\hat{\mathbf{Q}}_{i,j-1,k}) \\ \left(\frac{\partial \mathbf{\hat{H}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i,j,k+1/2} &=& \mathbf{C}_{i,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} + \mathbf{C}^-_{i,j,k+1/2}(\Delta\hat{\mathbf{Q}}_{i,j,k+1} - \Delta\hat{\mathbf{Q}}_{i,j,k}) \\ \left(\frac{\partial \mathbf{\hat{H}}}{\partial \hat{\mathbf{Q}}}\Delta \hat{\mathbf{Q}}\right)_{i,j,k-1/2} &=& \mathbf{C}_{i,j,k}\Delta\hat{\mathbf{Q}}_{i,j,k} - \mathbf{C}^+_{i,j,k-1/2}(\Delta\hat{\mathbf{Q}}_{i,j,k} - \Delta\hat{\mathbf{Q}}_{i,j,k-1}) \end{eqnarray}

where here \(\mathbf{A}^-_{i+1/2,j,k} = \frac{1}{2}(\tilde{\mathbf{A}}_{i+1/2,j,k} - |\tilde{\mathbf{A}}|_{i+1/2,j,k})\) and \(\mathbf{A}^+_{i-1/2,j,k} = \frac{1}{2}(\tilde{\mathbf{A}}_{i-1/2,j,k} + |\tilde{\mathbf{A}}|_{i-1/2,j,k})\), etc. To note here is that the flux on the left and right faces of a cell are computed using different expressions. That is, the expression for the flux at \(i-1/2\) is not simply that for the cell face at \(i+1/2\) offset one index. The central idea behind this approach is to increase the diagonal dominance of the system of equations and thereby the robustness.

Inviscid Flux Jacobians: Formulation I (PC)

Substituting the terms from Eqn. RoeFluxesPC into the expression for the flux Jacobian \(\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\) we obtain the following expression for the inviscid component:

\begin{eqnarray} \left(\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\right)_{i,j,k} &=& \mathbf{D^{I}}\Delta\hat{\mathbf{Q}}_{i,j,k} \\ &-& \mathbf{A}^-_{i+1/2,j,k}\Delta\hat{\mathbf{Q}}_{i+1,j,k} + \mathbf{A}^+_{i-1/2,j,k}\Delta\hat{\mathbf{Q}}_{i-1,j,k} \\ &-& \mathbf{B}^-_{i,j+1/2,k}\Delta\hat{\mathbf{Q}}_{i,j+1,k} + \mathbf{B}^+_{i,j-1/2,k}\Delta\hat{\mathbf{Q}}_{i,j-1,k} \\ &-& \mathbf{C}^-_{i,j,k+1/2}\Delta\hat{\mathbf{Q}}_{i,j,k+1} + \mathbf{C}^+_{i,j,k-1/2}\Delta\hat{\mathbf{Q}}_{i,j,k-1} \end{eqnarray}

where \(\mathbf{D^{I}} = (\mathbf{A}^-_{i-1/2,j,k} - \mathbf{A}^+_{i+1/2,j,k} + \mathbf{B}^-_{i,j-1/2,k} - \mathbf{B}^+_{i,j+1/2,k} + \mathbf{C}^-_{i,j,k-1/2} - \mathbf{C}^+_{i,j,k+1/2})\).

Inviscid Flux Jacobians: Formulation II (SGS)

Substituting the terms from Eqn. RoeFluxesSGS into the expression for the flux Jacobian \(\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\) we obtain the following expression for the inviscid component:

\begin{eqnarray} \left(\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\right)_{i,j,k} &=& \mathbf{D^{II}}\Delta\hat{\mathbf{Q}}_{i,j,k} \\ &-& \mathbf{A}^-_{i+1/2,j,k}\Delta\hat{\mathbf{Q}}_{i+1,j,k} + \mathbf{A}^+_{i-1/2,j,k}\Delta\hat{\mathbf{Q}}_{i-1,j,k} \\ &-& \mathbf{B}^-_{i,j+1/2,k}\Delta\hat{\mathbf{Q}}_{i,j+1,k} + \mathbf{B}^+_{i,j-1/2,k}\Delta\hat{\mathbf{Q}}_{i,j-1,k} \\ &-& \mathbf{C}^-_{i,j,k+1/2}\Delta\hat{\mathbf{Q}}_{i,j,k+1} + \mathbf{C}^+_{i,j,k-1/2}\Delta\hat{\mathbf{Q}}_{i,j,k-1} \end{eqnarray}

where \(\mathbf{D^{II}} = (\mathbf{A}^-_{i+1/2,j,k} - \mathbf{A}^+_{i-1/2,j,k} + \mathbf{B}^-_{i,j+1/2,k} - \mathbf{B}^+_{i,j-1/2,k} + \mathbf{C}^-_{i,j,k+1/2} - \mathbf{C}^+_{i,j,k-1/2})\).

Evaluation of Inviscid Flux Jacobians

The expressions \(\mathbf{A}^-_{i+1/2,j,k}\), \(\mathbf{A}^+_{i-1/2,j,k}\), please see the page about Low-Mach Preconditioning.

Viscous Flux Jacobians

The viscous Jacobians arise from the linearization of the viscous fluxes. These terms are in general not straightforward to obtain and furthermore are complicated by the presence of derivatives of the flow variables, which must be approximated using finite-differences. For these reasons simplified forms of the viscous Jacobians are preferred. In particularly the thin-shear-layer approximation of the Navier-Stokes equations with the assumption of locally constant dynamic viscosity and thermal conductivity coefficients is attractive as the viscous flux Jacobians can then be evaluated analytically.

In their most general form the viscous fluxes are dependent on both the flow variables \(\hat{\mathbf{Q}}\) and their spatial derivatives \(\frac{\partial\hat{\mathbf{Q}}}{\partial \eta}\), etc.:

\begin{eqnarray} \mathbf{\hat{F}}_v &=& \mathbf{\hat{F}}_v\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \xi}, \frac{\partial\hat{\mathbf{Q}}}{\partial \eta}, \frac{\partial\hat{\mathbf{Q}}}{\partial \zeta} \right) \\ \mathbf{\hat{G}}_v &=& \mathbf{\hat{G}}_v\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \xi}, \frac{\partial\hat{\mathbf{Q}}}{\partial \eta}, \frac{\partial\hat{\mathbf{Q}}}{\partial \zeta} \right) \\ \mathbf{\hat{H}}_v &=& \mathbf{\hat{H}}_v\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \xi}, \frac{\partial\hat{\mathbf{Q}}}{\partial \eta}, \frac{\partial\hat{\mathbf{Q}}}{\partial \zeta} \right) \end{eqnarray}

The viscous fluxes can however also be decoupled into components which depend on the vector of the dependent variables and its derivative in either the \(\xi\)-, \(\eta\)- or \(\zeta\)-direction:

\begin{eqnarray} \mathbf{\hat{F}}_v &=& \mathbf{\hat{F}}_{v_{\xi}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \xi}\right) + \mathbf{\hat{F}}_{v_{\eta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \eta}\right) + \mathbf{\hat{F}}_{v_{\zeta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \zeta}\right)\\ \mathbf{\hat{G}}_v &=& \mathbf{\hat{G}}_{v_{\xi}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \xi}\right) + \mathbf{\hat{G}}_{v_{\eta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \eta}\right) + \mathbf{\hat{G}}_{v_{\zeta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \zeta}\right) \\ \mathbf{\hat{H}}_v &=& \mathbf{\hat{H}}_{v_{\xi}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \xi}\right) + \mathbf{\hat{H}}_{v_{\eta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \eta}\right) + \mathbf{\hat{H}}_{v_{\zeta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \zeta}\right) \end{eqnarray}

In its simplest form the thin-layer approximation means that all derivatives not normal to the wall are to be neglected. For example, for a wall located at \(\eta = \) const, the viscous terms would reduce to the following form:

\begin{equation} \mathbf{\hat{G}}_v = \mathbf{\hat{G}}_{v_{\eta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \eta}\right) \end{equation}

where the terms \(\mathbf{\hat{F}}_v\) and \(\mathbf{\hat{H}}_v\) vanish as these appear in the Navier-Stokes equations in the form \(\frac{\partial \mathbf{\hat{F}}_v}{\partial \xi}\) and \(\frac{\partial \mathbf{\hat{H}}_v}{\partial \zeta}\).

In TRACE the use of the thin layer approximation amounts to neglecting the cross derivative terms in the viscous terms. The viscous fluxes therefore take the following form:

\begin{eqnarray} \mathbf{\hat{F}}_v &=& \mathbf{\hat{F}}_{v_{\xi}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \xi}\right)\\ \mathbf{\hat{G}}_v &=& \mathbf{\hat{G}}_{v_{\eta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \eta}\right) \\ \mathbf{\hat{H}}_v &=& \mathbf{\hat{H}}_{v_{\zeta}}\left(\hat{\mathbf{Q}}, \frac{\partial\hat{\mathbf{Q}}}{\partial \zeta}\right) \end{eqnarray}

As explained previously this simplifies the linearization of the flux terms. Considering the term \(\mathbf{\hat{F}}_v\), we now have

"Eqn_LinVisFluxes"

\begin{eqnarray} \mathbf{\hat{F}}_v^{m+1} &=& \mathbf{\hat{F}}_v^{m} + \frac{\partial \mathbf{\hat{F}}_v}{\partial \tau}\Delta\tau + O (\Delta\tau^2)\\ &=& \mathbf{\hat{F}}_v^{m} + \left(\frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}}\frac{\partial \hat{\mathbf{Q}}}{\partial \tau} + \frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\frac{\partial \hat{\mathbf{Q}}_{\xi}}{\partial \tau}\right)\Delta\tau + O (\Delta\tau^2)\\ &=& \mathbf{\hat{F}}_v^{m} + \frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}}\Delta\hat{\mathbf{Q}} + \frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\Delta \hat{\mathbf{Q}}_{\xi} \\ &=& \mathbf{\hat{F}}_v^{m} + \left[\frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}} - \frac{\partial}{\partial \xi}\left(\frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\right)\right]\Delta\hat{\mathbf{Q}} + \frac{\partial}{\partial \xi}\left(\frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\Delta \hat{\mathbf{Q}}\right) \\ &=& \mathbf{\hat{F}}_v^{m} + \frac{\partial}{\partial \xi}\left(\frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\Delta \hat{\mathbf{Q}}\right) \end{eqnarray}

where the term in square brackets is zero if the transport coefficients are approximated as locally constant. This derivation and the definitions of \(\frac{\partial \mathbf{\hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\), \(\frac{\partial \mathbf{\hat{G}}_v}{\partial \hat{\mathbf{Q}}_{\eta}}\) and \(\frac{\partial \mathbf{\hat{H}}_v} {\partial \hat{\mathbf{Q}}_{\zeta}}\) used in TRACE are taken from the work of Tysinger and Caughey [59] .

The viscous flux Jacobian is then

\begin{equation} \left.\frac{\partial}{\partial\xi}\left(\frac{\mathbf{\partial \hat{F}}_v}{\partial \hat{\mathbf{Q}}_{\xi}}\Delta \hat{\mathbf{Q}}\right)\right|_{i+1/2,j,k} = \left(\mathbf{A}_v \Delta \hat{\mathbf{Q}}\right)_{i+1,j,k} - \left(\mathbf{A}_v \Delta \hat{\mathbf{Q}}\right)_{i,j,k} \end{equation}

Viscous Flux Jacobians: Formulation

Substituting the terms from Eqn. Eqn_LinVisFluxes into the expression for the flux Jacobian \(\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\) and approximating the spatial derivatives on the cell faces with second-order central difference approximations we obtain the following expression for the viscous component:

\begin{eqnarray} \left(\frac{\partial \mathbf{R}}{\partial \hat{\mathbf{Q}}}\right)_{i,j,k} &=& \frac{1}{Re_0}\left[-\left((2\mathbf{A}_v + 2\mathbf{B}_v + 2\mathbf{C}_v)\Delta\hat{\mathbf{Q}}\right)_{i,j,k} \right. \\ &+& \left(\mathbf{A}_v\Delta\hat{\mathbf{Q}}\right)_{i+1,j,k} + \left(\mathbf{A}_v\Delta\hat{\mathbf{Q}}\right)_{i-1,j,k} \\ &+& \left(\mathbf{B}_v\Delta\hat{\mathbf{Q}}\right)_{i,j+1,k} + \left(\mathbf{B}_v\Delta\hat{\mathbf{Q}}\right)_{i,j-1,k} \\ &+& \left.\left(\mathbf{C}_v\Delta\hat{\mathbf{Q}}\right)_{i,j,k+1} + \left(\mathbf{C}_v\Delta\hat{\mathbf{Q}}\right)_{i,j,k-1}\right] \end{eqnarray}

Evaluation of Viscous Flux Jacobians

The viscous flux jacobians \(\mathbf{A}_v\), \(\mathbf{B}_v\) and \(\mathbf{C}_v\) are computed as follows following the work of Tysinger and Caughey [59] :

\begin{equation} \mathbf{A}_v, \mathbf{B}_v, \mathbf{C}_v = J\left( \begin{array}{ccccc} 0 & 0 & 0 & 0 & 0 \\ e_{21} & e_{22} & e_{23} & e_{24} & 0 \\ e_{31} & e_{32} & e_{33} & e_{34} & 0 \\ e_{41} & e_{42} & e_{43} & e_{44} & 0 \\ e_{51} & e_{52} & e_{53} & e_{54} & e_{21} \end{array} \right) \end{equation}

where

\begin{eqnarray} e_{21} &=& -\alpha_{xx}\left(\frac{u}{\rho}\right) - \alpha_{xy}\left(\frac{v}{\rho}\right) - \alpha_{xz}\left(\frac{w}{\rho}\right) \\ e_{22} &=& \alpha_{xx}\left(\frac{1}{\rho}\right) \\ e_{23} &=& \alpha_{xy}\left(\frac{1}{\rho}\right) \\ e_{24} &=& \alpha_{xz}\left(\frac{1}{\rho}\right) \end{eqnarray}

\begin{eqnarray} e_{31} &=& -\alpha_{xy}\left(\frac{u}{\rho}\right) - \alpha_{yy}\left(\frac{v}{\rho}\right) - \alpha_{yz}\left(\frac{w}{\rho}\right) \\ e_{32} &=& \alpha_{xy}\left(\frac{1}{\rho}\right) \\ e_{33} &=& \alpha_{yy}\left(\frac{1}{\rho}\right) \\ e_{34} &=& \alpha_{yz}\left(\frac{1}{\rho}\right) \end{eqnarray}

\begin{eqnarray} e_{41} &=& -\alpha_{xz}\left(\frac{u}{\rho}\right) - \alpha_{yz}\left(\frac{v}{\rho}\right) - \alpha_{zz}\left(\frac{w}{\rho}\right) \\ e_{42} &=& \alpha_{xz}\left(\frac{1}{\rho}\right) \\ e_{43} &=& \alpha_{yz}\left(\frac{1}{\rho}\right) \\ e_{44} &=& \alpha_{zz}\left(\frac{1}{\rho}\right) \end{eqnarray}

\begin{eqnarray} e_{51} &=& -\alpha_{xx}\left(\frac{u^2}{\rho}\right) - \alpha_{yy}\left(\frac{v^2}{\rho}\right) - \alpha_{zz}\left(\frac{w^2}{\rho}\right) \\ && - 2\alpha_{xy}\left(\frac{uv}{\rho}\right) - 2\alpha_{xz}\left(\frac{uw}{\rho}\right) - 2\alpha_{yz}\left(\frac{vw}{\rho}\right) \\ && + \alpha_e \left(-\frac{e}{\rho^2} + \frac{u^2 + v^2 + w^2}{\rho}\right) \\ e_{52} &=& -\alpha_{e}\left(\frac{u}{\rho}\right) - e_{21}\\ e_{53} &=& -\alpha_{e}\left(\frac{v}{\rho}\right) - e_{31} \\ e_{54} &=& -\alpha_{e}\left(\frac{w}{\rho}\right) - e_{41} \\ e_{55} &=& \alpha_{e}\left(\frac{1}{\rho}\right) \end{eqnarray}

with

\begin{eqnarray} \alpha_{xx} &=& (2\mu + \lambda)\kappa^2_x + \mu\kappa^2_y + \mu\kappa^2_z \\ \alpha_{yy} &=& (2\mu + \lambda)\kappa^2_y + \mu\kappa^2_z + \mu\kappa^2_x \\ \alpha_{zz} &=& (2\mu + \lambda)\kappa^2_z + \mu\kappa^2_x + \mu\kappa^2_y \\ \alpha_{xy} &=& (2\mu + \lambda)\kappa_x\kappa_y \\ \alpha_{xz} &=& (2\mu + \lambda)\kappa_x\kappa_z \\ \alpha_{yz} &=& (2\mu + \lambda)\kappa_y\kappa_z \\ \alpha_e &=& \gamma\mu {Pr}^{-1}\left(\kappa^2_x + \kappa^2_y + \kappa^2_z\right) \end{eqnarray}

and \(\kappa = \xi\), \(\eta\) or \(\zeta\) for \(\mathbf{A}_v\), \(\mathbf{B}_v\) or \(\mathbf{C}_v\), respectively.

Source Term

As a final step the source terms

\begin{equation} \mathbf{\hat{S}} = \frac{1}{J}\left[ \begin{array}{c} 0 \\ 0 \\ \rho \omega_1 (y \omega_1 + 2w) \\ \rho \omega_1 (z \omega_1 - 2v) \\ 0 \end{array} \right] = \frac{1}{J}\left[ \begin{array}{c} 0 \\ 0 \\ \omega_1 (Q_1y \omega_1 + 2Q_4) \\ \omega_1 (Q_1z \omega_1 - 2Q_3) \\ 0 \end{array} \right] \end{equation}

must be linearized and included in the left hand side. The linearization requires the evaluation of source Jacobian:

\begin{equation} \frac{\partial \mathbf{S}}{\partial \mathbf{Q}} = \left(\begin{array}{ccccc} \frac{\partial S_1}{\partial Q_1} & \frac{\partial S_1}{\partial Q_2} & \frac{\partial S_1}{\partial Q_3} & \frac{\partial S_1}{\partial Q_4} & \frac{\partial S_1}{\partial Q_5} \\ \frac{\partial S_2}{\partial Q_1} & \frac{\partial S_1}{\partial Q_2} & \frac{\partial S_1}{\partial Q_3} & \frac{\partial S_1}{\partial Q_4} & \frac{\partial S_1}{\partial Q_5} \\ \frac{\partial S_3}{\partial Q_1} & \frac{\partial S_1}{\partial Q_2} & \frac{\partial S_1}{\partial Q_3} & \frac{\partial S_1}{\partial Q_4} & \frac{\partial S_1}{\partial Q_5} \\ \frac{\partial S_4}{\partial Q_1} & \frac{\partial S_1}{\partial Q_2} & \frac{\partial S_1}{\partial Q_3} & \frac{\partial S_1}{\partial Q_4} & \frac{\partial S_1}{\partial Q_5} \\ \frac{\partial S_5}{\partial Q_1} & \frac{\partial S_1}{\partial Q_2} & \frac{\partial S_1}{\partial Q_3} & \frac{\partial S_1}{\partial Q_4} & \frac{\partial S_1}{\partial Q_5} \end{array}\right) \end{equation}

where \(\mathbf{Q} = [\rho\:\:\rho u\:\:\rho v\:\:\rho w\:\:\rho E]^T\). Evaluating the individual terms we obtain the follow

\begin{equation} \frac{\partial \mathbf{S}}{\partial \mathbf{Q}} = \mathbf{A}^{s} = \frac{1}{J}\left( \begin{array}{ccccc} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ y {\omega_1}^2 & 0 & 0 & 2\omega_1 & 0 \\ z{\omega_1}^2 & 0 & -2\omega_1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{array} \right) \end{equation}

LHS Matrix

The implicit matrix takes the following general form:

"Eqn_LHS-MatrixGeneral"

\begin{equation}\label{Eqn_LHS-MatrixGeneral} \left[ \begin{array}{cccccccccccccc} \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} & 0 & \cdots & 0 & \mathbf{G} & 0 & 0 & 0 & 0 \\ \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} & 0 & \cdots & 0 & \mathbf{G} & 0 & 0 & 0 \\ 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} & 0 & \cdots & 0 & \ddots & 0 & 0 \\ \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} & 0 & \cdots & 0 & \mathbf{G} & 0 \\ 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} & 0 & \cdots & 0 & \mathbf{G} \\ \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} & 0 & \cdots & 0 \\ 0 & \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} & 0 & \cdots \\ \cdots & 0 & \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} & 0 \\ 0 & \cdots & 0 & \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 & \mathbf{E} \\ \mathbf{F} & 0 & \cdots & 0 & \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots & 0 \\ 0 & \mathbf{F} & 0 & \cdots & 0 & \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 & \cdots \\ 0 & 0 & \mathbf{F} & 0 & \cdots & 0 & \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} & 0 \\ 0 & 0 & 0 & \mathbf{F} & 0 & \cdots & 0 & \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} & \mathbf{C} \\ 0 & 0 & 0 & 0 & \mathbf{F} & 0 & \cdots & 0 & \mathbf{D} & 0 & \cdots & 0 & \mathbf{A} & \mathbf{B} \\ \end{array} \right] \end{equation}

The entries in the matrix are themselves 5x5 matrices. The precise form of these depends on the discretization scheme used and the approximations used for the Jacobians. Using the two formulations outlined previously for the inviscid Jacobians the following two expressions are obtained for the terms in Eqn. Eqn_LHS-MatrixGeneral.

Formulation I: (PC)

\begin{eqnarray} \mathbf{A}_{i,j,k} &=& -a_{ii}\left(\mathbf{A}^{v}_{i-1,j,k} + \mathbf{A}^+_{i-1/2,j,k}\right) \\ \mathbf{B}_{i,j,k} &=& \left(\frac{1}{J\Delta \tau} + \frac{1}{J\Delta t}\right)\mathbf{I} \\ &-& a_{ii}\left(\mathbf{D^{I}}_{i,j,k} - 2(\mathbf{A}^{v}_{i,j,k} + \mathbf{B}^{v}_{i,j,k} + \mathbf{C}^{v}_{i,j,k}) + \mathbf{A}^{s}_{i,j,k}\right) \\ &=& \left(\frac{1}{J\Delta \tau} + \frac{1}{J\Delta t}\right)\mathbf{I} \\ &-& \mathbf{A}_{i+1,j,k} - \mathbf{C}_{i-1,j,k} - \mathbf{D}_{i,j+1,k} \\ &-& \mathbf{E}_{i,j-1,k} - \mathbf{F}_{i,j,k+1} - \mathbf{G}_{i,j,k-1} - a_{ii}\mathbf{A}^{s}_{i,j,k} \\ \mathbf{C}_{i,j,k} &=& -a_{ii}\left(\mathbf{A}^{v}_{i+1,j,k} - \mathbf{A}^-_{i+1/2,j,k}\right) \\ \mathbf{D}_{i,j,k} &=& -a_{ii}\left(\mathbf{B}^{v}_{i,j-1,k} + \mathbf{B}^+_{i,j-1/2,k}\right) \\ \mathbf{E}_{i,j,k} &=& -a_{ii}\left(\mathbf{B}^{v}_{i,j+1,k} - \mathbf{B}^-_{i,j+1/2,k}\right) \\ \mathbf{F}_{i,j,k} &=& -a_{ii}\left(\mathbf{C}^{v}_{i,j,k-1} + \mathbf{C}^+_{i,j,k-1/2}\right) \\ \mathbf{G}_{i,j,k} &=& -a_{ii}\left(\mathbf{C}^{v}_{i,j,k+1} - \mathbf{C}^-_{i,j,k+1/2}\right) \end{eqnarray}

where

\begin{eqnarray} \mathbf{D^{I}}_{i,j,k} &=& \mathbf{A}^-_{i-1/2,j,k} - \mathbf{A}^+_{i+1/2,j,k} + \mathbf{B}^-_{i,j-1/2,k} \\ &-& \mathbf{B}^+_{i,j+1/2,k} + \mathbf{C}^-_{i,j,k-1/2} - \mathbf{C}^+_{i,j,k+1/2} \end{eqnarray}

\begin{eqnarray} \mathbf{A}^{+}_{i+1/2,j,k} &=& \frac{1}{2}(\mathbf{A}_{i,j,k} + |\tilde{\mathbf{A}}|_{i+1/2,j,k}) \\ \mathbf{A}^{+}_{i-1/2,j,k} &=& \frac{1}{2}(\mathbf{A}_{i-1,j,k} + |\tilde{\mathbf{A}}|_{i-1/2,j,k}) \\ \mathbf{A}^{-}_{i+1/2,j,k} &=& \frac{1}{2}(\mathbf{A}_{i+1,j,k} - |\tilde{\mathbf{A}}|_{i+1/2,j,k}) \\ \mathbf{A}^{-}_{i-1/2,j,k} &=& \frac{1}{2}(\mathbf{A}_{i,j,k} - |\tilde{\mathbf{A}}|_{i-1/2,j,k}) \end{eqnarray}

\begin{eqnarray} \mathbf{B}^{+}_{i,j+1/2,k} &=& \frac{1}{2}(\mathbf{B}_{i,j,k} + |\tilde{\mathbf{B}}|_{i,j+1/2,k}) \\ \mathbf{B}^{+}_{i,j-1/2,k} &=& \frac{1}{2}(\mathbf{B}_{i,j-1,k} + |\tilde{\mathbf{B}}|_{i,j-1/2,k}) \\ \mathbf{B}^{-}_{i,j+1/2,k} &=& \frac{1}{2}(\mathbf{B}_{i,j+1,k} - |\tilde{\mathbf{B}}|_{i,j+1/2,k}) \\ \mathbf{B}^{-}_{i,j-1/2,k} &=& \frac{1}{2}(\mathbf{B}_{i,j,k} - |\tilde{\mathbf{B}}|_{i,j-1/2,k}) \end{eqnarray}

\begin{eqnarray} \mathbf{C}^{+}_{i,j,k+1/2} &=& \frac{1}{2}(\mathbf{C}_{i,j,k} + |\tilde{\mathbf{C}}|_{i,j,k+1/2}) \\ \mathbf{C}^{+}_{i,j,k-1/2} &=& \frac{1}{2}(\mathbf{C}_{i,j,k-1} + |\tilde{\mathbf{C}}|_{i,j,k-1/2}) \\ \mathbf{C}^{-}_{i,j,k+1/2} &=& \frac{1}{2}(\mathbf{C}_{i,j,k+1} - |\tilde{\mathbf{C}}|_{i,j,k+1/2}) \\ \mathbf{C}^{-}_{i,j,k-1/2} &=& \frac{1}{2}(\mathbf{C}_{i,j,k} - |\tilde{\mathbf{C}}|_{i,j,k-1/2}) \end{eqnarray}

Formulation II: (SGS)

\begin{eqnarray} \mathbf{A}_{i,j,k} &=& -a_{ii}\left(\mathbf{A}^{v}_{i-1,j,k} + \mathbf{A}^+_{i-1/2,j,k}\right) \\ \mathbf{B}_{i,j,k} &=& \left(\frac{1}{J\Delta \tau} + \frac{1}{J\Delta t}\right)\mathbf{I} \\ &-& a_{ii}\left(\mathbf{D^{II}}_{i,j,k} - 2(\mathbf{A}^{v}_{i,j,k} + \mathbf{B}^{v}_{i,j,k} + \mathbf{C}^{v}_{i,j,k}) + \mathbf{A}^{s}_{i,j,k}\right) \\ \mathbf{C}_{i,j,k} &=& -a_{ii}\left(\mathbf{A}^{v}_{i+1,j,k} - \mathbf{A}^-_{i+1/2,j,k}\right) \\ \mathbf{D}_{i,j,k} &=& -a_{ii}\left(\mathbf{B}^{v}_{i,j-1,k} + \mathbf{B}^+_{i,j-1/2,k}\right) \\ \mathbf{E}_{i,j,k} &=& -a_{ii}\left(\mathbf{B}^{v}_{i,j+1,k} - \mathbf{B}^-_{i,j+1/2,k}\right) \\ \mathbf{F}_{i,j,k} &=& -a_{ii}\left(\mathbf{C}^{v}_{i,j,k-1} + \mathbf{C}^+_{i,j,k-1/2}\right) \\ \mathbf{G}_{i,j,k} &=& -a_{ii}\left(\mathbf{C}^{v}_{i,j,k+1} - \mathbf{C}^-_{i,j,k+1/2}\right) \end{eqnarray}

where

\begin{eqnarray} \mathbf{D^{II}}_{i,j,k} &=& \mathbf{A}^-_{i+1/2,j,k} - \mathbf{A}^+_{i-1/2,j,k} + \mathbf{B}^-_{i,j+1/2,k} \\ &-& \mathbf{B}^+_{i,j-1/2,k} + \mathbf{C}^-_{i,j,k+1/2} - \mathbf{C}^+_{i,j,k-1/2} \\ \mathbf{A}^{\pm}_{i+1/2,j,k} &=& \frac{1}{2}(\tilde{\mathbf{A}}_{i+1/2,j,k} \pm |\tilde{\mathbf{A}}|_{i+1/2,j,k}) \\ \mathbf{A}^{\pm}_{i-1/2,j,k} &=& \frac{1}{2}(\tilde{\mathbf{A}}_{i-1/2,j,k} \pm |\tilde{\mathbf{A}}|_{i-1/2,j,k}) \\ \mathbf{B}^{\pm}_{i,j+1/2,k} &=& \frac{1}{2}(\tilde{\mathbf{B}}_{i,j+1/2,k} \pm |\tilde{\mathbf{B}}|_{i,j+1/2,k}) \\ \mathbf{B}^{\pm}_{i,j-1/2,k} &=& \frac{1}{2}(\tilde{\mathbf{B}}_{i,j-1/2,k} \pm |\tilde{\mathbf{B}}|_{i,j-1/2,k}) \\ \mathbf{C}^{\pm}_{i,j,k+1/2} &=& \frac{1}{2}(\tilde{\mathbf{C}}_{i,j,k+1/2} \pm |\tilde{\mathbf{C}}|_{i,j,k+1/2}) \\ \mathbf{C}^{\pm}_{i,j,k-1/2} &=& \frac{1}{2}(\tilde{\mathbf{C}}_{i,j,k-1/2} \pm |\tilde{\mathbf{C}}|_{i,j,k-1/2}) \end{eqnarray}