Fluid Dynamics

Conservation equations

In fluid domains, the code solves the compressible Navier-Stokes equations in curvilinear coordinates. The basic equations, in a cartesian coordinate space, for the conserved mass density $\rho$ , momentum density $\rho u_i$ , and total energy density $\rho E$ are, in index form with summation convention are given as

$\begin{align} \frac{\partial \rho}{\partial t} + \frac{\partial }{\partial x_j} \rho u_j &= S_\rho \\ \frac{\partial \rho u_i}{\partial t} + \frac{\partial}{\partial x_j}\left(\rho u_i u_j + p\delta_{ij} - \tau_{ij}\right) &= S_{\rho u_i} \\ \frac{\partial \rho E}{\partial t} + \frac{\partial}{\partial x_j}\left(\left\{\rho E + p\right\}u_j + q_j - u_i \tau_{ij}\right) &= S_{\rho E} \end{align}$

where $p$ is the thermodynamic pressure, $\tau_{ij}$ is the viscous stress tensor, and $q_i$ is the heat flux in the $i$ th direction. $S_\rho$ , $S_{\rho u_i}$ , and $S_{\rho E}$ are are mass, momentum, and energy density source terms. These equations can be written in the compact form

$\frac{\partial Q}{\partial t} + \frac{\partial \vec{F}_j}{\partial x_j} = S$

where $Q = [\rho\,\rho \vec{u}\,\rho E]^T$ is the vector of conserved variables, $\vec{F} = \vec{F}^I - \vec{F}^V$ is the flux vector account for both visicd and invisci terms, and $S$ is the source term vector. The so-called RHS $= S - \frac{\partial}{\partial x_j}\vec{F}_j$ is evaluated by the top-level RHS driver function in include/EulerRHS.H.

Viscous stress constitutive relation

The viscous stress tensor is defined as

$\tau_{ij} = \mu \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) + \lambda \frac{\partial u_k}{\partial x_k}\delta_{ij}$

where $\mu$ and $\lambda$ are the first and second coefficients of viscosity, respectively; both may be a function of temperature. Note that Stokes' hypothesis $(\lambda = -\frac{2}{3}\mu)$ is not automatically enforced and that $\lambda$ is related to bulk viscosity $\mu_B$ as $\lambda = \mu_B - (2/3)\mu$ .

$\tau_{ij}$ is computed by ComputeTauBuffer in src/ViscidUtil.C.

Heat flux constitutive relation

The heat flux vector is defined as

$q_i = - \kappa \frac{\partial T}{\partial x_i}$

where $\kappa$ is the thermal conductivity. The heat flux vector, $q_i$ , is computed by ComputeHeatFluxBuffer in src/ViscidUtil.C.

Transport Coefficient Models

The first viscosity coefficient $\mu$ , bulk viscosity coefficient, $\mu_B$ , and the thermal conductivity $k$ depend on the thermodynamic state of the fluid. Currently, only one implementation is available, although others can be easily implemented in ViscidUtil.C.

Power Law

The Power Law transport model in implemented in ComputeTVBufferPower.

The power law model gives the dynamic viscosity, $\mu$ as

$\mu = \beta T^n$

where $\beta$ and $n$ are user specified parameters, typically $n = 0.666$ and $\beta = 4.093 x 10^{-7}$ for air.

The bulk viscosity is defined as

$\mu_B = \alpha \mu$

where $\alpha$ is a user specified parameter, typically $\alpha = 0.6$ for air.

Thus the second coefficient of viscosity can be calculated as

$\lambda = \left(\alpha - 2/3\right) \mu$

The power law model calculates the :w

Equations of state

The equations of state provides closure by relating the intensive state variables, pressure and temperature, to the extensive state variables, specific internal energy and volume. Currently, only one implementation is available, although others can be easily implemented in EulerUtil.C.

Calorically perfect ideal gas

The equation of state currently available is that of an ideal gas, assuming constant specific heats. The equations of state are

$P = \rho R T$

where $R$ is the specific gas constant, defined as $R = R_u / W$ with $R_u$ the universal gas constant, and $W$ the molecular weight.

The specific heat capacity at constant volume and pressure are defined as

$\begin{align} C_v &= \left(\frac{\partial E}{\partial T}\right)_v \\ C_p &= \left(\frac{\partial H}{\partial T}\right)_p \end{align}$

Then, by substitution into the equation of state we get the following relation

$R = C_p - C_v$

By defining the specific heat ratio, $\gamma = \frac{C_p}{C_v}$ , the following expressions give the realtionship between specific energy, pressure, and temperature as implemented by ComputeDVBufferIdeal.

$\begin{align} P &= (\gamma -1) \rho e \\ T &= \frac{\gamma-1}{R} e \end{align}$

Non-dimensionalization

PlasCom2 can run in either a dimensional or non-dimensional mode. The code uses the following variables to define the non-dimensional scaling:

$\rho^*_\infty$ , $P^*_\infty$ , $T^*_\infty$ , and $L^*$ , a length scale. Where $*$ denotes a dimensional value and $\infty$ denotes the reference state. There are two optional non-dimensional spaces available to the user, as shown in the table below.

PlasCom2 Non-Dimensional spaces
Standard (nonDimensional=1)	Legacy PlasComCM (nonDimensional=2)
$u^_\infty = \sqrt \frac{P^_\infty}{\rho^*_\infty}$	$u^_\infty = \sqrt \frac{\gamma P^_\infty}{\rho^*_\infty}$
$e^_\infty = u^_\infty^2 = \frac{P^_\infty}{\rho^_\infty}$	$e^_\infty = u^_\infty^2 = \frac{\gamma P^_\infty}{\rho^_\infty}$
$\rho = \rho^* /\rho^*_\infty$	$\rho = \rho^* /\rho^*_\infty$
$P = P^* /P^*_\infty$	$P = P^* /(\rho^_\infty u^_\infty^2)$
$T = T^* /T^*_\infty$	$T = T^* /((\gamma-1)T^*_\infty)$
$u_i = u^_i /u^_\infty$	$u_i = u^_i /u^_\infty$
$e = e^* /e^*_\infty$	$e = e^* /e^*_\infty$
$t = t^* /(L^* / u^*_\infty)$	$t = t^* /(L^* / u^*_\infty)$

Substitution into the dimensional form of the Navier-Stokes equations yields the non-dimensional equivalent

$\begin{align} \frac{\partial \rho}{\partial t} + \frac{\partial }{\partial x_j} \rho u_j &= S_\rho \\ \frac{\partial \rho u_i}{\partial t} + \frac{\partial}{\partial x_j}\left(\rho u_i u_j + p\delta_{ij} - \tau_{ij}\right) &= S_{\rho u_i} \\ \frac{\partial \rho E}{\partial t} + \frac{\partial}{\partial x_j}\left(\left\{\rho E + p\right\}u_j + q_j - u_i \tau_{ij}\right) &= S_{\rho E} \end{align}$

with the following non-dimensionalization for the source terms

$\begin{align} S_\rho &= \frac{S^*_\rho L^*}{\rho^*_\infty U^*_\infty} \\ S_{\rho u_i} &= \frac{S^*_{\rho u_i } L^*}{\rho^*_\infty U^*_\infty^2 } \\ S_{\rho E} &= \frac{S^*_{\rho E} L^*}{\rho^*_\infty U^*_\infty^3} \end{align}$

by choosing the following non-dimensionalizations for the transport coefficients

$\begin{align} \mu &= \mu^* /\mu^*_\infty \\ \lambda &= \lambda^* /\lambda^*_\infty \\ \kappa &= \kappa^* /\kappa^*_\infty \\ \end{align}$

the non-dimensional viscous stress tensor and heat flux vector can be written as

$\begin{align} \tau_{ij} &= \frac{\mu}{\RE} \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) + \frac{\lambda}{\RE} \frac{\partial u_k}{\partial x_k}\delta_{ij} \\ q_i &= - \frac{\mu}{\RE \Pr} \frac{\partial T}{\partial x_i} \end{align}$

where $\ensuremath{\mathit{Re}}$ is defined as the code Reynolds number, $\RE = \frac{\rho^*_\infty U^*_\infty L^*}{\mu^*_\infty}$ and $\ensuremath{\mathit{Pr}}$ is defined as the Prandtl number, $\PR = \frac{C^*_p_\infty\mu^*_\infty}{k^*_\infty} = \frac{C_p\mu}{k}$ which define the dimensional reference values $\mu^*_\infty$ and $\kappa^*_\infty$ respectively.

Non-dimensional equation of state

There are no special modifications to the calorically perfect gas equation of state, with the exception of the specific gas constant. The reference gas constant is calculated and non-dimensionalized as follows

$\begin{align} R^*_\infty &= \frac{P^*_\infty}{\rho^*_\infty T^*_\infty} \\ R &= R^* /R^*_\infty \\ \end{align}$

For the standard non-dimensionalization, $R$ is exactly 1.0. For the legacy non-dimensionalization, $R = \frac{\gamma-1}{\gamma}$ .

Curvilinear coordinate systems

It is possible to express the Navier-Stokes equations in any other coordinate system $\xi_i$ provided the mapping $\Xi$ from $x_i$ to $\xi_i$ is given. The Cartesian coordinates $(\vec{x},t)$ can be mapped to another coordinate system $(\vec{\xi},\tau)$ via the time-dependent mappings

$\vec{x} = {X}(\vec{\xi},\tau) \qquad \mbox{ with inverse } \qquad \vec{\xi} = \Xi(\vec{x},t)$

where $X^{-1} = \Xi$ and we only consider non-singular mappings such that $X^{-1}$ exists and is well defined. Moreover we take $t = \tau$ . The Jacobian of the transformation is defined as $J = \mathrm{det}(\partial \Xi_i/\partial x_j)$ and is strictly positive.

Under these conditions and with simple application of the chain rule it can be shown[2] that the vector (compact) form tranforms to

$\frac{\partial}{\partial \tau}\left(\frac{Q}{J}\right) + \frac{\partial \hat{\vec{F}}^I_i}{\partial \xi_i} - \frac{\partial \hat{\vec{F}}^V_i}{\partial \xi_i}= \frac{S}{J}$

after using the identities

$\begin{split} \frac{\partial}{\partial \xi_j}\left(\frac{1}{J}\frac{\partial \xi_j}{\partial x_i}\right) &= 0 \qquad \mbox{ for $i = 1, \dots, N$} \\ \frac{\partial}{\partial \tau}\left(\frac{1}{J}\right) + \frac{\partial}{\partial \xi_j}\left(\frac{1}{J}\frac{\partial \xi_j}{\partial t}\right) &= 0, \end{split}$

where $N$ is the number of dimensions. If we define the weighted metric $\hat{\xi_i} = J^{-1}(\partial \xi/\partial x_i)$ and contravariant velocity $\hat{U} = u_j \hat{\xi}_j + \hat{\xi}_t$ , with similar expressions for the remaining components, then the inviscid fluxes $\hat{F}^I_i$ are

$\hat{\vec{F}}^I_1 = \begin{bmatrix} \rho \hat{U} \\ \rho u \hat{U} + p\hat{\xi}_x \\ \rho v \hat{U} + p \hat{\xi}_y \\ (\rho E + p)\hat{U} - \hat{\xi}_t p \end{bmatrix}\qquad \mbox{ and } \qquad \hat{\vec{F}}^I_2 = \begin{bmatrix} \rho \hat{V} \\ \rho u \hat{V} + p\hat{\eta}_x \\ \rho v \hat{V} + p \hat{\eta}_y \\ (\rho E + p)\hat{V} - \hat{\eta}_t p \end{bmatrix}$

in two dimensions and

$\hat{\vec{F}}^I_1 = \begin{bmatrix} \rho \hat{U} \\ \rho u \hat{U} + p\hat{\xi}_x \\ \rho v \hat{U} + p \hat{\xi}_y \\ \rho w \hat{U} + p\hat{\xi}_z \\ (\rho E + p)\hat{U} - \hat{\xi}_t p \end{bmatrix}, \quad \hat{\vec{F}}^I_2 = \begin{bmatrix} \rho \hat{V} \\ \rho u \hat{V} + p\hat{\eta}_x \\ \rho v \hat{V} + p \hat{\eta}_y \\ \rho w \hat{V} + p \hat{\eta}_z \\ (\rho E + p)\hat{V} - \hat{\eta}_t p \end{bmatrix}, \quad \mbox{and} \quad \hat{\vec{F}}^I_3 = \begin{bmatrix} \rho \hat{W} \\ \rho u \hat{W} + p\hat{\zeta}_x \\ \rho w \hat{W} + p \hat{\zeta}_y \\ \rho w \hat{W} + p \hat{\zeta}_z \\ (\rho E + p)\hat{W} - \hat{\zeta}_t p \end{bmatrix}$

The inviscid fluxes are computed dimension-at-a-time by Euler::Flux1D in kernels/Euler.f90.

Forms of the viscous terms

The viscous fluxes may be expressed in a number of forms, depending on the particular goal. The main difference between the particular forms is how second derivatives are taken; namely, either $\partial^2/\partial x^2$ is taken directly or as two repeated derivatives, $\partial/\partial x\left(\partial/\partial x\right)$ . The use of $\partial^2/\partial x^2$ is advantageous in that it allows for the most physical dissipation in the code, as determined by the modified wavenumber. (See, e.g., Lele[1] for a discussion of the modified wavenumber and its meaning.) However, it is also expensive, especially in three or more dimensions. Using repeated derivatives in advantageous for two reasons: (i) it keeps the equations in conservative form which is useful for shock capturing and (ii) it is less expensive (by a factor around 2.5 in three dimensions) than the fully expanded form. However, it provides less physical dissipation, especially at the highest wavenumbers, exactly where it is most needed.

Strong form of the viscous terms in xi-coordinates

Following Anderson, Tanehill, and Pletcher (1984), the strong form of the viscous terms is as follows

$\begin{align} \frac{\partial}{\partial t}\left(\frac{\rho u_1}{J}\right) &= \cdots \frac{\partial}{\partial \xi}\left(\hat{\xi}_i\tau_{i1}\right) + \frac{\partial}{\partial \eta}\left(\hat{\eta}_i\tau_{i1}\right) + \frac{\partial}{\partial \zeta}\left(\hat{\zeta}_i\tau_{i1}\right) \\ \frac{\partial}{\partial t}\left(\frac{\rho u_2}{J}\right) &= \cdots \frac{\partial}{\partial \xi}\left(\hat{\xi}_i\tau_{i2}\right) + \frac{\partial}{\partial \eta}\left(\hat{\eta}_i\tau_{i2}\right) + \frac{\partial}{\partial \zeta}\left(\hat{\zeta}_i\tau_{i2}\right) \\ \frac{\partial}{\partial t}\left(\frac{\rho u_3}{J}\right) &= \cdots \frac{\partial}{\partial \xi}\left(\hat{\xi}_i\tau_{i3}\right) + \frac{\partial}{\partial \eta}\left(\hat{\eta}_i\tau_{i3}\right) + \frac{\partial}{\partial \zeta}\left(\hat{\zeta}_i\tau_{i3}\right) \\ \frac{\partial}{\partial t}\left(\frac{\rho E}{J}\right) &= \cdots \frac{\partial}{\partial \xi}\left(\hat{\xi}_i [ u_j \tau_{ij} - q_i ]\right) + \frac{\partial}{\partial \eta}\left(\hat{\eta}_i [ u_j \tau_{ij} - q_i ]\right) + \frac{\partial}{\partial \zeta}\left(\hat{\zeta}_i [ u_j \tau_{ij} - q_i ]\right) \end{align}$

This form is much faster than the other forms, but does not have any dissipation at the highest wavenumbers. The viscous flux components are computed dimension-at-a-time by Viscid::StrongFlux1D in kernels/Viscid.f90.

Selected by the user with METRICTYPE = NONORTHOGONAL_STRONGFORM in .