Governing Equations

The compressible Navier-Stokes equations, written in strong conservation law (divergence) form are

\[\begin{split}\begin{gather*} \frac{\partial \rho}{\partial t} + \frac{\partial \rho u_j}{\partial x_j} = 0\\[0.1in] \frac{\partial \rho u_i}{\partial t} + \frac{\partial \rho u_i u_j}{\partial x_j} = -\frac{\partial p}{\partial x_i} + \frac{\partial \sigma_{ij}}{\partial x_j}\\[0.1in] \frac{\partial \rho E}{\partial t} + \frac{\partial (\rho E + p )u_j}{\partial x_j} = \frac{\partial u_i\sigma_{ij}}{\partial x_j} - \frac{\partial q_j}{\partial x_j} \end{gather*}\end{split}\]

where \(\sigma_{ij}\) and \(q_j\) are the viscous stress and thermal conduction, defined as

\[\begin{split}\begin{eqnarray} \sigma_{ij} & = & \mu\left[\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_i}{\partial x_j}\right) - \frac{2}{3}\left(\frac{\partial u_k}{\partial x_k}\right) \delta_{ij}\right],\\[0.1in] q_j & = & -k\frac{\partial T}{\partial x_j}, \end{eqnarray}\end{split}\]

and where \(\mu\) is the dynamic viscosity, \(k\) is the thermal conductivity, and \(\delta_{ij}\) is the Kronecker delta.

The solution variables are the density, \(\rho\), the momentum per unit volume \(\rho u_i\), and the total energy per unit volume \(\rho E = \rho(e + 1/2u_ku_k)\), where \(e=c_vT\) is the internal energy per unit mass. Here \(c_v\) is the heat capacity at constant volume and \(T\) is the temperature. The pressure \(p\) appears as an auxiliary variable and it is related to the density and temperature through the ideal gas law

\[p = \rho R T,\]

Where \(R\) is the gas constant.

In the finite volume framework we start with integral forms of the conservation laws, not the differential forms given above. The integral forms can be derived in a variety of ways, including simply integrating the differential forms over a control volume. Doing this and making use of Gauss' theorem leads to

\[\begin{split}\begin{gather*} \frac{d}{dt}\int_\mathcal{V}\rho\: d\mathcal{V} \:+\; \int_\mathcal{S} \rho u_j n_j\: d\mathcal{S} \;=\; 0 \\[0.1in] \frac{d}{dt}\int_\mathcal{V}\rho u_i\: d\mathcal{V} \:+\; \int_\mathcal{S} \rho u_iu_j n_j\: d\mathcal{S} \;= \; \int_\mathcal{S} (-p\delta_{ij} \:+\; \sigma_{ij})n_j\: d\mathcal{S}\\[0.1in] \frac{d}{dt}\int_\mathcal{V}\rho E\: d\mathcal{V} \:+\; \int_\mathcal{S} (\rho E + p)u_j n_j\: d\mathcal{S} \;=\; \int_\mathcal{S} u_i\sigma_{ij}n_j\: d\mathcal{S} \:-\: \int_\mathcal{S} q_jn_j\: d\mathcal{S} \end{gather*}\end{split}\]

where \(\mathcal{V}\) is a control volume and \(\mathcal{S}\) is its bounding with a local outward-pointing unit normal vector \(n_j\).

The solution variables are stored at the computational cell centroids and we assume that they take on linear variations between points. When products of solution variables are considered for the flux terms, we interpolate the appropriate factors of the product independently and then retain only the net linear variations in the plane of the face. Using these approximations, the linear variations integrate to zero when volume or surface integrals are considered, thereby producing simple discrete forms. For example, the discrete form of the momentum equation is

\[\frac{d\rho u_i}{dt} = \frac{1}{\mathcal{V}}\sum_{l=1}^{6~faces} \left[-(\rho u_i)_l(u_j)_l - (p)_l\delta_{ij} + (\sigma_{ij})_l\right] \left(n_j\Delta S\right)_l\]

where \(l\) is an index that covers the six faces of a hexahedral computational cell. Here the linear variations have been integrated away, leaving just the cell-center value of \(\rho u_i\) on the left and the face-center values of the flux terms on the right. The latter are determined via simple interpolations from the cell center values, i.e. \((\rho u_i)_l\) means the quantity \(\rho u_i\) interpolated to face \(l\).