The transformation of the governing equations (NSE) to curvlinear coordinates
\begin{align} \tau = t, \qquad \xi = \xi(x,y,z,t), \qquad \eta = \eta(x,y,z,t), \qquad \zeta = \zeta(x,y,z,t) \end{align}
yields the following transformed Navier-Stokes equations in conservative form
\begin{align} \newcommand{\ddelt } [1] {\frac{\partial #1}{\partial t}} \newcommand{\ddelxi } [1] {\frac{\partial #1}{\partial \xi}} \newcommand{\ddeleta } [1] {\frac{\partial #1}{\partial \eta}} \newcommand{\ddelzeta } [1] {\frac{\partial #1}{\partial \zeta}} \ddelt{\mv{\hat{Q}}} + \ddelxi{\mv{\hat{E}}} + \ddeleta{\mv{\hat{F}}} + \ddelzeta{\mv{\hat{G}}} = 0 \end{align}
where
\begin{align} \mv{\hat{Q}} = J\mv{Q} \end{align}
and
\begin{align} \mv{\hat{E}} &= J\left(\xi_t\mv{Q} + \xi_x\mv{E} + \xi_y\mv{F} + \xi_z\mv{G}\right)\\ \mv{\hat{F}} &= J\left(\eta_t\mv{Q} + \eta_x\mv{E} + \eta_y\mv{F} + \eta_z\mv{G}\right)\\ \mv{\hat{G}} &= J\left(\zeta_t\mv{Q} + \zeta_x\mv{E} + \zeta_y\mv{F} + \zeta_z\mv{G}\right) \end{align}
The quantities \xi_i, \eta_i\, and \zeta_i denote metric terms and J denotes the Jacobian.
To compute the transformed fluxes \mv{\hat{E}}, \mv{\hat{F}}, and \mv{\hat{G}} with second-order accuracy central differences around the cell center are computed via
\begin{align} \left.\frac{\partial \mv{\hat{E}}}{\partial \xi}\right|_{i,j,k} &= \mv{\hat{E}}_{i+\frac{1}{2},j,k} - \mv{\hat{E}}_{i-\frac{1}{2},j,k}\\ \left.\frac{\partial \mv{\hat{F}}}{\partial \eta}\right|_{i,j,k} &= \mv{\hat{F}}_{i,j+\frac{1}{2},k} - \mv{\hat{F}}_{i,j-\frac{1}{2},k}\\ \left.\frac{\partial \mv{\hat{G}}}{\partial \zeta}\right|_{i,j,k} &= \mv{\hat{G}}_{i,j,k+\frac{1}{2}} - \mv{\hat{G}}_{i,j,k-\frac{1}{2}}\,. \end{align}
The convective flux is computed with the Advection Upstream Splitting Method (AUSM) by [LiouSteffen1993]. In this method the inviscid flux is split into a convective part, computed by an upwind formulation, and a pressure part, computed by central differences. The inviscid flux at a cell interface, i.e., the surface that connects two cell volumes, by
\begin{align} \newcommand{\vecfive } [5] {\begin{pmatrix} #1\\ #2\\ #3\\ #4\\ #5 \end{pmatrix}} \mv{\hat{E}}_{inv} = J\Big[ \vecfive{\rho U}{\rho Uu}{\rho Uv}{\rho Uw}{\rho U(E + p/\rho)} + \vecfive{0}{p\xi_x}{p\xi_y}{p\xi_z}{p\xi_t}\Big] \end{align}
where the contravariant velocity is given by
\begin{align} U = u\xi_x + v\xi_y + w\xi_z + \xi_{t}\,. \end{align}
The introduction of the speed of sound a yields
\begin{align} \newcommand{\vecfive } [5] {\begin{pmatrix} #1\\ #2\\ #3\\ #4\\ #5 \end{pmatrix}} \mv{\hat{E}}_{L/R} = \underbrace{ \frac{U}{\left|\nabla \xi \right|a} }_{M_{L/R}} \underbrace{ J\left|\nabla \xi\right| \vecfive{\rho a}{\rho au}{\rho av}{\rho aw}{\rho a\left(E + p/\rho\right)} }_{\mv{\hat{E}}^c_{L/R}} + \underbrace{ J\vecfive{0}{p\xi_x}{p\xi_y}{p\xi_z}{p\xi_t} }_{\mv{\hat{E}}^p_{L/R}} \end{align}
where L and R denote the left and right state of the surface. That is, the left or right state the interface Mach number
\begin{align} M_{L/R} = \frac{1}{2}\left[M_L + M_R\right] \end{align}
where
\begin{align} M_{L,R} = \frac{U_{L,R}}{\left|\nabla \xi \right|a_{L,R}} \end{align}
and
\begin{align} a_{L,R} = \sqrt{\gamma\frac{p_{L,R}}{\rho_{L,R}}} \end{align}
Whether the left or right state is chosen for convection depends on the Mach number, i.e.,
\begin{align} \mv{\hat{E}}^c_{L/R} = \left\{ \begin{array}{ll} (.)_L, &M_{L/R} \geq 0\\ (.)_R, &M_{L/R} < 0 \end{array} \right.\,. \end{align}
As an example, the total discretized inviscid flux in the \xi-direction is computed by
\begin{align} \mv{\hat{E}}_{L/R} = \frac{1}{2}\left[M_{L/R}\left(\mv{\hat{E}}^c_L + \mv{\hat{E}}^c_R\right) + |M_{L/R}| \left(\mv{\hat{E}}^c_L - \mv{\hat{E}}^c_R\right)\right] + \mv{\hat{E}}^p_{L/R} \end{align}
where
\begin{align} \mv{\hat{E}}^p_{L/R} = \frac{1}{2}\left[p_L + p_R\right]J\vecfive{0}{\xi_x}{\xi_y}{\xi_z}{\xi_t}\,. \end{align}
is the pressure part of the flux. To achieve a second-order accuracy the left and right states are computed by a monotonic upstream-centered scheme for conservation laws (MUSCL) [vanLeer1979], such that
\begin{align} \newcommand{\im } [0] {{i - \frac{1}{2}}} \newcommand{\ip } [0] {{i + \frac{1}{2}}} \left.\mv{P}_L\right|_{i+\frac{1}{2}} &= \mv{P}_i + \frac{l_\ip}{\left(l_\ip + l_\im\right)^2} \left[l_\im \left(\mv{P}_{i+1} -\mv{P}_{i}\right) + l_\ip\left(\mv{P}_{i} -\mv{P}_{i-1} \right)\right]\\ \left.\mv{P}_R\right|_{i+\frac{1}{2}} &= \mv{P}_{i-1} - \frac{l_\ip}{\left(l_\ip + l_{i+\frac{3}{2}}\right)^2} \left[l_\ip \Delta_{i+\frac{3}{2}} \left(\mv{P}_{i+2} - \mv{P}_{i+1}\right) + l_{i+\frac{3}{2}}\left(\mv{P}_{i+1} - \mv{P}_{i}\right)\right] \end{align}
where
\begin{align} l_{i\pm\frac{1}{2}} = \sqrt{\left(\pm\mv{x}_{i\pm 1} \mp \mv{x}_{i}\right)\cdot\left(\pm\mv{x}_{i\pm 1} \mp \mv{x}_{i}\right)} \end{align}
is the distance function between the cell centers of neighbour cells in the chosen direction.
The viscous flux is computed by a modified cell-vertex (MCV) scheme [Meinke1993], where the flux through each interface, e.g., \left(i+\frac{1}{2},j,k\right) and \left(i-\frac{1}{2},j,k\right) is determined by an average of the flux on the four cell interface corner points. For instance, for the \mv{E}-flux
\begin{align} \left.\frac{\partial\mv{\hat{E}}_{vis}}{\partial \xi}\right|_{i,j,k} = &\frac{1}{4}\left( \mv{\hat{E}}_{vis,i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}} + \mv{\hat{E}}_{vis,i+\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}} + \mv{\hat{E}}_{vis,i+\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} + \mv{\hat{E}}_{vis,i+\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}} \right)\\ - &\frac{1}{4}\left( \mv{\hat{E}}_{vis,i-\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}} + \mv{\hat{E}}_{vis,i-\frac{1}{2},j-\frac{1}{2},k+\frac{1}{2}} + \mv{\hat{E}}_{vis,i-\frac{1}{2},j+\frac{1}{2},k-\frac{1}{2}} + \mv{\hat{E}}_{vis,i-\frac{1}{2},j-\frac{1}{2},k-\frac{1}{2}} \right) \end{align}
The corner fluxes are computed by central differences using the variables on the surrounding eight cell centers, e.g., for \frac{\partial u}{\partial \xi} at the corner point \left(i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}\right)
\begin{align} \left.u_{\xi}\right|_{i+\frac{1}{2},j+\frac{1}{2},k+\frac{1}{2}} =&\frac{1}{4}\left(u_{i+1,j,k} + u_{i+1,j+1,k} + u_{i+1,j,k+1} + u_{i+1,j+1,k+1} \right) -\\ &\frac{1}{4}\left(u_{i,j,k} + u_{i,j+1,k} + u_{i,j,k+1} + u_{i,j+1,k+1}\right)\,. \end{align}
Combining the discretized spatial derivatives into the discrete right-hand-side operator RHS\left(t;\cdot \right) the semi-discrete equation
\begin{align} \frac{\partial \mv{\hat{Q}}}{\partial t} = RHS\left(t;\mv{\hat{Q}}\right) \end{align}
is then integrated in time by a Runge-Kutta scheme, such that
\begin{align} \begin{array}{ccl} \mv{\hat{Q}}^{(0)} &=& \mv{\hat{Q}}^n\\ &\vdots& \\ \mv{\hat{Q}}^{(k)} &=& \mv{\hat{Q}}^{(0)} + \alpha_k\Delta t RHS\left(t^n + \alpha_{k-1}\Delta t; \mv{\hat{Q}}^{(k-1)}\right)\quad k = 1,\dots ,s\,,\\ &\vdots& \\ \mv{\hat{Q}}^{n+1} &=& \mv{\hat{Q}}^{(k)} \end{array} \end{align}
In most computations the coefficients are set to \alpha_k = (1/4, 1/6, 3/8, 1/2, 1) such that a second-order accurarcy in time is achieved.
1.9.5