Fluid Equations
Primary field variables for the fluid equations in ChiDG are Density,
Momentum-1, Momentum-2, Momentum-3, and Energy
\[\begin{split}Q_{cart} =
\begin{pmatrix}
\rho \\ \rho v_x \\ \rho v_y \\ \rho v_z \\ \rho E
\end{pmatrix}
\quad\quad
Q_{cyl} =
\begin{pmatrix}
\rho \\ \rho v_r \\ r \rho v_\theta \\ \rho v_z \\ \rho E
\end{pmatrix}\end{split}\]
Note
In Cylindrical coordinates, angular momentum, \(r \rho v_\theta\),
is being conserved. Not tangential momentum, \(\rho v_\theta\).
The transformation of the momentum vector between Cartesian and
Cylindrical coordinates is
\[\begin{split}\begin{bmatrix}
\rho v_x \\ \rho v_y \\ \rho v_z
\end{bmatrix}
=
\begin{bmatrix}
cos(\theta) & -sin(\theta) & 0 \\
sin(\theta) & cos(\theta) & 0 \\
0 & 0 & 1
\end{bmatrix}
\begin{bmatrix}
\rho v_r \\ \rho v_\theta \\ \rho v_z
\end{bmatrix}\end{split}\]
The velocity vector in Cartesian and Cylindrical coordinates is
\[\vec{v} = [v_x, v_y, v_z] \quad\quad \vec{v} = [v_r, v_\theta, v_z]\]
Some distinctions are made here between absolute, relative, and frame velocities.
The absolute velocity(\(v\)) is that which is perceived in the inertial frame of
reference. The relative(or advection)(\(w\)) velocity is the velocity, which is perceived
in the non-inertial frame of reference. The frame velocity(\(u\)) is the velocity of
the moving coordinate system. These are defined in the Cartesian and Cylindrical
coordinate systems as
\[\begin{split}\vec{v} &= [v_x, v_y, v_z] \quad\quad\:\:\:\: \vec{v} = [v_r, v_\theta, v_y] \\
\vec{w} &= [w_x, w_y, w_z] \quad\quad \vec{w} = [w_r, w_\theta, w_y] \\
\vec{u} &= [u_x, u_y, u_z] \quad\quad\:\:\: \vec{u} = [u_r, u_\theta, u_y]\end{split}\]
For calculations using an inertial frame, the frame velocity is zero.
\[\vec{u} = [0, 0, 0]\]
For a rotating Cylindrical coordinate system, the frame velocity is
\[\vec{u}_{cyl} = [0, \omega r, 0]\]
where \(\omega\) is the rotational rate in \(rad/s\). The relation
between absolute, relative, and frame velocities is
\[\vec{v} = \vec{w} + \vec{u}\]
At times, components of the velocity vectors are represented in the same manner
for both coordinate systems
\[\vec{v} = [v_1,v_2,v_3]\]
Within the framework, effort has been invested in representing quantities
such that they are consistent with vector calculus. As such, representations
of a gradient in Cartesian and Cylindrical coordinates as
\[\begin{split}\nabla f &= \bigg[\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}\bigg]
\\
\\
\nabla f &= \bigg[\frac{\partial f}{\partial r}, \frac{1}{r}\frac{\partial f}{\partial \theta}, \frac{\partial f}{\partial z}\bigg]\end{split}\]
are represented in terms of the components of the gradient vector as
\[\nabla f = [\nabla_1 f, \nabla_2 f, \nabla_3 f]\]
Fluid Advection
The fluxes representing advective transport for the fluid equations
are given compactly as
\[\begin{split}\vec{F}_a(Q) =
\begin{pmatrix}
\rho \vec{w} \\
\mathcal{S} (\rho \vec{v} \otimes \vec{w} + \overline{\overline{I}}p) \\
\rho H \vec{w} + \vec{u} p
\end{pmatrix}
\quad\quad
\mathcal{S}_{cart} =
\begin{pmatrix}
1 \\ 1 \\ 1
\end{pmatrix}
\quad\quad
\mathcal{S}_{cyl} =
\begin{pmatrix}
1 \\ r \\ 1
\end{pmatrix}\end{split}\]
where \(\mathcal{S}\) is a scaling vector that is sued to provide the correct
dimensional scaling for the angular momentum equation. Expanded forms for the
fluxes in Cartesian and Cylindrical coordinates are given as follows.
Cartesian
\[\begin{split}\vec{F}_a(\vec{x},Q) =
\begin{pmatrix}
\begin{bmatrix}
\rho w_x \\
\rho v_x w_x + p \\
\rho v_y w_x \\
\rho v_z w_x \\
\rho H w_x + u_x p
\end{bmatrix}
,
\begin{bmatrix}
\rho w_y \\
\rho v_x w_y \\
\rho v_y w_y + p \\
\rho v_z w_y \\
\rho H w_y + u_y p
\end{bmatrix}
,
\begin{bmatrix}
\rho w_z \\
\rho v_x w_z \\
\rho v_y w_z \\
\rho v_z w_z + p \\
\rho H w_z + u_z p
\end{bmatrix}
\end{pmatrix}
\quad\quad
S_a(\vec{x},Q) =
\begin{pmatrix}
0 \\ 0 \\ 0 \\ 0 \\ 0
\end{pmatrix}\end{split}\]
Cylindrical
\[\begin{split}\vec{F}_a(\vec{r},Q) =
\begin{pmatrix}
\begin{bmatrix}
\rho w_r \\
\rho v_r w_r + p \\
r \rho v_\theta w_r \\
\rho v_z w_r \\
\rho H w_r + u_r p
\end{bmatrix}
,
\begin{bmatrix}
\rho w_\theta \\
\rho v_r w_\theta \\
r \rho v_\theta w_\theta + r p \\
\rho v_z w_\theta \\
\rho H w_\theta + u_\theta p
\end{bmatrix}
,
\begin{bmatrix}
\rho w_z \\
\rho v_r w_z \\
r \rho v_\theta w_z \\
\rho v_z w_z + p \\
\rho H w_z + u_z p
\end{bmatrix}
\end{pmatrix}
\quad\quad
S_a(\vec{r},Q) =
\begin{pmatrix}
0 \\ \frac{\rho v_\theta w_\theta - p}{r} \\ 0 \\ 0 \\ 0
\end{pmatrix}\end{split}\]
\(H = \frac{\rho E + p}{\rho}\) is the total enthalpy and \(p\) is the static pressure.
Fluid Diffusion
The fluxes representing diffusive transport for the fluid equations are given compactly as
\[\begin{split}\vec{F}_d(Q,\nabla Q) =
-
\begin{pmatrix}
0 \\ \mathcal{S} \overline{\overline{\tau}} \\ k \nabla T + \overline{\overline{\tau}} \cdot \vec{v}
\end{pmatrix}
\quad\quad
\mathcal{S}_{cart} =
\begin{pmatrix}
1 \\ 1 \\ 1
\end{pmatrix}
\quad\quad
\mathcal{S}_{cyl} =
\begin{pmatrix}
1 \\ r \\ 1
\end{pmatrix}\end{split}\]
where \(\mathcal{S}\) is a scaling vector that is used to provide the correct dimensional
scaling for the angular momentum equation. Expanded forms for the fluxes in Cartesian
and Cylindrical coordinates are given as follows.
Cartesian
\[\begin{split}\vec{F}_d(\vec{x},Q,\nabla Q) =
-
\begin{pmatrix}
\begin{bmatrix}
0 \\
\tau_{11} \\
\tau_{21} \\
\tau_{31} \\
k \nabla_1 T + \overline{\overline{\tau}} \cdot \vec{v}
\end{bmatrix}
,
\begin{bmatrix}
0 \\
\tau_{12} \\
\tau_{22} \\
\tau_{32} \\
k \nabla_2 T + \overline{\overline{\tau}} \cdot \vec{v}
\end{bmatrix}
,
\begin{bmatrix}
0 \\
\tau_{13} \\
\tau_{23} \\
\tau_{33} \\
k \nabla_3 T + \overline{\overline{\tau}} \cdot \vec{v}
\end{bmatrix}
\end{pmatrix}
\quad
S_d(\vec{x},Q) =
\begin{pmatrix}
0 \\ 0 \\ 0 \\ 0 \\ 0
\end{pmatrix}\end{split}\]
Cylindrical
\[\begin{split}\vec{F}_d(\vec{r},Q,\nabla Q) =
-
\begin{pmatrix}
\begin{bmatrix}
0 \\
\tau_{11} \\
r\tau_{21} \\
\tau_{31} \\
k \nabla_1 T + \overline{\overline{\tau}} \cdot \vec{v}
\end{bmatrix}
,
\begin{bmatrix}
0 \\
\tau_{12} \\
r\tau_{22} \\
\tau_{32} \\
k \nabla_2 T + \overline{\overline{\tau}} \cdot \vec{v}
\end{bmatrix}
,
\begin{bmatrix}
0 \\
\tau_{13} \\
r\tau_{23} \\
\tau_{33} \\
k \nabla_3 T + \overline{\overline{\tau}} \cdot \vec{v}
\end{bmatrix}
\end{pmatrix}
\quad
S_d(\vec{r},Q) =
\begin{pmatrix}
0 \\ -\frac{\tau_{22}}{r} \\ 0 \\ 0 \\ 0
\end{pmatrix}\end{split}\]
Models
Equations of State
Ideal Gas
| Model name: |
Ideal Gas |
| Model fields: |
Pressure Temperature |
The ideal gas equation of state computes Pressure and Temperature as
\[\begin{split}p &= (\gamma - 1)\bigg(\rho E - \frac{1}{2}\frac{\vec{\rho v} \cdot \vec{\rho v}}{\rho} \bigg)
\\
\\
T &= \frac{p}{\rho R}\end{split}\]
where \(R\) is the specific gas constant
\[R = 287.15 \quad\quad \bigg[ \frac{\text{J}}{\text{kg} \cdot \text{K}} \bigg]\]
Viscosity
Sutherland’s Law
| Model name: |
Sutherlands Law |
| Model fields: |
Dynamic Viscosity |
Sutherland’s Law computes Dynamic Viscosity as a function of temperature
using
\[\mu = \mu_0 \bigg(\frac{T}{T_0}\bigg)^{3/2} \frac{T_0 + S}{T + S}\]
where the model constants are
\[\begin{split}\mu_0 &= 1.7894e^{-5} \quad\quad &\bigg[\frac{\text{kg}}{\text{m} \cdot \text{s}}\bigg] \\
T_0 &= 273.11 \quad\quad &\big[K\big] \\
S &= 110.56 \quad\quad &\big[K\big]\end{split}\]
Constant Viscosity
Velocity Gradients
Gradients of velocity are computed using the chain rule. From the ChiDG framework, we have gradients of the
primary field variables. Here, we have gradients of the components of
momentum: \((\nabla \rho v_1,\nabla \rho v_2, \nabla \rho v_3)\). Gradients
of velocity, \((\nabla v_1, \nabla v_2, \nabla v_3)\) are computed by recognizing that in general
\[\nabla (\phi f) = \nabla(\phi) f + \phi \nabla (f)\]
So the gradient of the velocity components can be computed as
\[\begin{split}\nabla (v_1) &= \frac{\nabla(\rho v_1)}{\rho} - \frac{v_1 \nabla(\rho)}{\rho}
\\
\nabla (v_2) &= \frac{\nabla(\rho v_2)}{\rho} - \frac{v_2 \nabla(\rho)}{\rho}
\\
\nabla (v_3) &= \frac{\nabla(\rho v_3)}{\rho} - \frac{v_3 \nabla(\rho)}{\rho}\end{split}\]
Note
In Cylindrical coordinates, we have \(\nabla(r \rho v_\theta)\) instead
of \(\nabla(\rho v_\theta)\). The gradient of tangential momentum is
computed from the angular momentum gradient as
\[\nabla(\rho v_\theta) =
\begin{bmatrix}
\frac{\nabla_1(r \rho v_\theta)}{r} - \frac{v_\theta}{r}, &
\frac{\nabla_2(r \rho v_\theta)}{r}, &
\frac{\nabla_3(r \rho v_\theta)}{r}
\end{bmatrix}\]
Shear Stress
| Model name: |
Shear Stress |
| Model fields: |
Shear-11, Shear-22, Shear-33, Shear-12, Shear-13, Shear-23 |
The shear stress tensor is defined as
\[\overline{\overline{\tau}} = \mu(\nabla \vec{v} + \nabla \vec{v}^T) + \lambda \overline{\overline{I}} \nabla \cdot \vec{v}\]
The tensor compnents are
\[\begin{split}\overline{\overline{\tau}} =
\begin{bmatrix}
\tau_{11} & \tau_{12} & \tau_{13} \\
\tau_{21} & \tau_{22} & \tau_{23} \\
\tau_{31} & \tau_{32} & \tau_{33} \\
\end{bmatrix}\end{split}\]
Note
The stress tensor is symmetric. So, only the upper triangular components of the tensor
are computed.
The components of the stress tensor are computed as
Cartesian
\[\begin{split} \tau_{11} &= 2 \mu \bigg(\nabla_1 v_x \bigg) + \lambda(\nabla \cdot \vec{v}) \\
\tau_{22} &= 2 \mu \bigg(\nabla_2 v_y \bigg) + \lambda(\nabla \cdot \vec{v}) \\
\tau_{33} &= 2 \mu \bigg(\nabla_3 v_z \bigg) + \lambda(\nabla \cdot \vec{v}) \\
\\
\tau_{12} &= \mu \bigg( \nabla_2 v_x + \nabla_1 v_y \bigg) \\
\tau_{13} &= \mu \bigg( \nabla_3 v_x + \nabla_1 v_z \bigg) \\
\tau_{23} &= \mu \bigg( \nabla_2 v_z + \nabla_3 v_y \bigg)\end{split}\]
\[\nabla \cdot \vec{v} = \bigg( \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} + \frac{\partial v_z}{\partial z} \bigg) = \bigg[ \nabla_1 v_x + \nabla_2 v_y + \nabla_3 v_z \bigg]\]
Cylindrical
\[\begin{split} \tau_{11} &= 2 \mu \bigg(\nabla_1 v_r \quad\quad \bigg) + \lambda(\nabla \cdot \vec{v}) \\
\tau_{22} &= 2 \mu \bigg(\nabla_2 v_\theta + \frac{v_r}{r} \bigg) + \lambda(\nabla \cdot \vec{v}) \\
\tau_{33} &= 2 \mu \bigg(\nabla_3 v_z \quad\quad \bigg) + \lambda(\nabla \cdot \vec{v}) \\
\\
\tau_{12} &= \mu \bigg( \nabla_2 v_r + \nabla_1 v_\theta - \frac{v_\theta}{r} \bigg) \\
\tau_{13} &= \mu \bigg( \nabla_3 v_r + \nabla_1 v_z \quad\quad \bigg) \\
\tau_{23} &= \mu \bigg( \nabla_2 v_z + \nabla_3 v_\theta \quad\quad \bigg)\end{split}\]
\[\nabla \cdot \vec{v} = \bigg( \frac{1}{r}\frac{\partial r v_r}{\partial r} + \frac{1}{r}\frac{\partial v_\theta}{\partial \theta} + \frac{\partial v_z}{\partial z}\bigg) = \bigg( \frac{\partial v_r}{\partial r} + \frac{1}{r}\frac{\partial v_\theta}{\partial \theta} + \frac{\partial v_z}{\partial z} + \frac{v_r}{r} \bigg) = \bigg[\nabla_1 v_r + \nabla_2 v_\theta + \nabla_3 v_z + \frac{v_r}{r}\bigg]\]
Temperature Gradient
| Model name: |
Temperature Gradient |
| Model fields: |
Temperature Gradient - 1 Temperature Gradient - 2 Temperature Gradient - 3 |
Gradients of temperature are computed using the chain rule. From the ChiDG framework,
we have gradients of the primary field variables. The gradient of the scalar
temperature field \(\nabla T(Q)\) can be computed by expanding
\[\nabla T(Q) = \frac{\partial T}{\partial Q} \nabla Q\]
as
\[\nabla T(\rho, \rho v_1, \rho v_2, \rho v_3, \rho E) =
\frac{\partial T}{\partial \rho} \nabla \rho +
\frac{\partial T}{\partial \rho v_1} \nabla \rho v_1 +
\frac{\partial T}{\partial \rho v_2} \nabla \rho v_2 +
\frac{\partial T}{\partial \rho v_3} \nabla \rho v_3 +
\frac{\partial T}{\partial \rho E} \nabla \rho E\]
Note
In Cylindrical coordinates, we have \(\nabla(r \rho v_\theta)\) instead
of \(\nabla(\rho v_\theta)\). The gradient of tangential momentum is
computed from the angular momentum gradient as
\[\nabla(\rho v_\theta) =
\begin{bmatrix}
\frac{\nabla_1(r \rho v_\theta)}{r} - \frac{v_\theta}{r}, &
\frac{\nabla_2(r \rho v_\theta)}{r}, &
\frac{\nabla_3(r \rho v_\theta)}{r}
\end{bmatrix}\]
Vorticity
| Model name: |
Vorticity |
| Model fields: |
Vorticity-1 Vorticity-2 Vorticity-3 |
Vorticity, used in the Spalart-Allmaras turbulence model, is defined as the Curl of
velocity as
\[\vec{\omega} = \nabla \times \vec{v}\]
In Cartesian coordinates, this is computed as
\[\begin{split}\vec{\omega} = \nabla \times \vec{v} &= \bigg(\frac{\partial v_z}{\partial x} - \frac{\partial v_y}{\partial z}\bigg) \hat{x} + \bigg(\frac{\partial v_x}{\partial z} - \frac{\partial v_z}{\partial x}\bigg) \hat{y} + \bigg( \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y}\bigg) \hat{z} \\
&= \bigg[ \bigg(\nabla_2 v_3 - \nabla_3 v_2\bigg), \bigg(\nabla_3 v_1 - \nabla_1 v_3\bigg), \bigg(\nabla_1 v_2 - \nabla_2 v_1\bigg)\bigg]\end{split}\]
In Cylindrical coordinates, the Curl of a vector is given as
\[\begin{split}\vec{\omega} = \nabla \times \vec{v} &= \bigg(\frac{1}{r}\frac{\partial v_z}{\partial \theta} - \frac{\partial v_\theta}{\partial z}\bigg) \hat{r} + \bigg(\frac{\partial v_r}{\partial z} - \frac{\partial v_z}{\partial r}\bigg) \hat{\theta} + \frac{1}{r}\bigg( \frac{\partial r v_\theta}{\partial r} - \frac{\partial v_r}{\partial \theta}\bigg) \hat{z} \\
&= \bigg(\frac{1}{r}\frac{\partial v_z}{\partial \theta} - \frac{\partial v_\theta}{\partial z}\bigg) \hat{r} + \bigg(\frac{\partial v_r}{\partial z} - \frac{\partial v_z}{\partial r}\bigg) \hat{\theta} + \bigg( \frac{\partial v_\theta}{\partial r} - \frac{1}{r}\frac{\partial v_r}{\partial \theta} + \frac{v_\theta}{r}\bigg) \hat{z} \\
&= \bigg[ \bigg(\nabla_2 v_3 - \nabla_3 v_2\bigg), \bigg(\nabla_3 v_1 - \nabla_1 v_3\bigg), \bigg(\nabla_1 v_2 - \nabla_2 v_1 + \frac{v_2}{r}\bigg)\bigg]\end{split}\]
In the non-inertial frame for Cylindrical coordinates, the relative
vorticity is accounted for as
\[\omega_3 = \omega_3 - 2\omega\]
Note, that here \(\omega_3\) is the third component of the vorticity
vector, while \(\omega\) is the rate of rotation for the non-inertial frame.