One-Dimensional Euler Solver
CFD Julia: A Learning Module Structuring an Introductory Course on Computational Fluid Dynamics
The one-dimensional Euler equations in its conservative form can be written as
\[\cfrac{\partial \mathbf{q}}{\partial t}+\cfrac{\partial \mathbf{F}}{\partial x}=0\]
where
\[\begin{split}\mathbf{q}=\begin{bmatrix}
\rho\\ \rho u \\\rho E
\end{bmatrix}\quad
\mathbf{F}=\begin{bmatrix}
\rho u\\ \rho uu + p \\\rho H u
\end{bmatrix}\end{split}\]
in which
\[\begin{split}\begin{array}{l}
H=E+\cfrac{p}{\rho}\\
E=\cfrac{1}{\gamma-1}\cfrac{p}{\rho}+\cfrac{1}{2}\mathbf{v}^2
\end{array}\end{split}\]
Here, \(\rho\) and \(p\) are the density, and pressure respectively; \(u\) is the horizontal component of velocity;
\(E\) and \(H\) stand for the specific internal energy and specific enthalpy(‘specific’ means per-unit-mass), and \(\gamma\) is the ratio of specific heats. The above equationscan also
be written as
\[\cfrac{\partial \mathbf{q}}{\partial t}+\mathbf{A}\cfrac{\partial \mathbf{q}}{\partial x}=0\]
where \(\mathbf{A}=\cfrac{\partial \mathbf{F}}{\partial \mathbf{q}}\) is the convective flux Jacobian matrix. The matrix
\(\mathbf{A}\) is given as
\[\begin{split}\mathbf{A}=\cfrac{\partial \mathbf{F}}{\partial \mathbf{q}}=\begin{bmatrix}
0& 1 & 0\\
(\gamma-3)\cfrac{u^2}{2}& (3-\gamma)u & (\gamma-1)\\
(\cfrac{\gamma-1}{2}u^{2}-H)u&H+(1-\gamma)u^{2} &\gamma u
\end{bmatrix}\end{split}\]
or
\[\begin{split}\mathbf{A}=\cfrac{\partial \mathbf{F}}{\partial \mathbf{q}}=\begin{bmatrix}
0& 1 & 0\\
\phi ^{2}-u^2& (3-\gamma)u & (\gamma-1)\\
(\phi ^{2}-H)u&H+(1-\gamma)u^{2} &\gamma u
\end{bmatrix}\end{split}\]
where \(\phi^{2}=\cfrac{\gamma-1}{2}q^2\).
in which
\[\begin{split}\begin{array}{l}
q^2=u^{2}\quad\quad\quad\quad\quad\quad(1d)\\
q^2=u^{2}+v^{2}\quad\quad\quad\quad(2d)\\
q^2=u^{2}+v^{2}+w^{2}\ \ \ \quad(3d)\\
\end{array}\end{split}\]
Eigenstructure
\[\mathbf{A}=\mathbf{R}\Lambda \mathbf{L}\]
where
\[\begin{split}\Lambda=\begin{bmatrix}
u & 0 & 0\\
0& u+c & 0\\
0& 0 &u-c
\end{bmatrix}\end{split}\]
\[\begin{split}\mathbf{R}=\begin{bmatrix}
1 & \beta & \beta\\
u& \beta (u+c) & \beta (u-c)\\
\cfrac{\phi^{2}}{\gamma-1}& \beta (H+uc) &\beta (H-uc)
\end{bmatrix}\end{split}\]
\[\begin{split}\mathbf{L}=\begin{bmatrix}
1- \cfrac{\phi^{2}}{c^2}& (\gamma-1)\cfrac{u}{c^{2}} & -\cfrac{(\gamma-1)}{c^{2}}\\
\phi^{2}-uc& -(\gamma-1)u+c & \gamma-1\\
\phi^{2}+uc& -(\gamma-1)u-c &\gamma-1
\end{bmatrix}\end{split}\]
where \(c\) is the speed of sound and is given by \(c^{2}=\gamma p/\rho\), and \(\beta=1/(2c^{2})\).
The semi-discrete form of the Euler equations can be written as
\[\cfrac{\partial \mathbf{q}_{i}}{\partial t}+\cfrac{\mathbf{F}_{i+1/2}-\mathbf{F}_{i-1/2}}{\Delta x}=0\]
where \(\mathbf{q}_{i}\) is cell-centered values stored at nodal points and \(\mathbf{F}_{i-1/2}\) and \(\mathbf{F}_{i+1/2}\) are
the fluxes at left and right cell interface.
We use third-order Runge-Kutta numerical method for time integration. We use WENO-5 reconstruction to compute the left and sight states of the fluxes at the interface.
Roe’s Riemann Solver
\[\mathbf{F}_{i+1/2}=\cfrac{1}{2}(\mathbf{F}^{L}_{i+1/2}+\mathbf{F}^{R}_{i+1/2})
-\cfrac{1}{2}\bar{\mathbf{R}}|\bar{\Lambda}|\bar{\mathbf{L}}(\mathbf{q}^{R}_{i+1/2}-\mathbf{q}^{L}_{i+1/2})\]
First we reconstruct the left and right states of \(\mathbf{q}\) at the interface (i.e., \(\mathbf{q}^{L}_{i+1/2}\) and \(\mathbf{q}^{R}_{i+1/2}\) ) using WENO-5 scheme.
Then, we can compute the left and right state of the fluxes (i.e., \(\mathbf{F}^{L}\) and \(\mathbf{F}^{R}\)).
The Roe averaging formulas to compute approximate values for constructing \(\bar{\mathbf{R}}\) and \(\bar{\mathbf{L}}\)
are given below
\[\begin{split}\begin{array}{c}
\bar{u}=\cfrac{u^{L}\sqrt{\rho^{L}}+u^{R}\sqrt{\rho^{R}}}{\sqrt{\rho^{L}}+\sqrt{\rho^{R}}}\\
\bar{H}=\cfrac{H^{L}\sqrt{\rho^{L}}+H^{R}\sqrt{\rho^{R}}}{\sqrt{\rho^{L}}+\sqrt{\rho^{R}}}
\end{array}\end{split}\]
where the left and right states are computed using WENO-5 reconstruction. The speed of the sound is computed using the below equation
\[\bar{c}=\sqrt{(\gamma-1)\left[\bar{H}-\cfrac{1}{2}\bar{u}^{2}\right]}\]
The eigenvalues of the Jacobian matrix are \(\lambda_{1}=\bar{u}\), \(\lambda_{2}=\bar{u}+\bar{c}\), and \(\lambda_{3}=\bar{u}-\bar{c}\).
\[\begin{split}\begin{array}{l}
\mathbf{q}[0]=\rho\\
\mathbf{q}[1]=\rho u\\
\mathbf{q}[2]=\rho E\\
\end{array}\end{split}\]
\[\begin{split}\begin{array}{l}
\rho E=\cfrac{p}{\gamma-1}+\cfrac{1}{2}\rho (u^2+v^2+w^2)\\
(\gamma-1)\rho E=p+\cfrac{(\gamma-1)}{2}\rho (u^2+v^2+w^2)\\
p=(\gamma-1)\rho E-\cfrac{(\gamma-1)}{2}\rho (u^2+v^2+w^2)\\
p=(\gamma-1)\rho E-\cfrac{(\gamma-1)}{2} \cfrac{((\rho u)^2+(\rho v)^2+(\rho w)^2)}{\rho}\\
p=(\gamma-1)\mathbf{q}[2]-\cfrac{(\gamma-1)}{2} \cfrac{(\mathbf{q}[1]^2)}{\mathbf{q}[0]}\\
p=(\gamma-1)(\mathbf{q}[2]-\cfrac{1}{2} \cfrac{\mathbf{q}[1]^2}{\mathbf{q}[0]})\\
\end{array}\end{split}\]
\[\begin{split}u=\cfrac{\rho u}{\rho}=\mathbf{q}[1]/\mathbf{q}[0]\\\end{split}\]
\[\begin{split}\begin{array}{l}
H=E+\cfrac{p}{\rho}\\
{\rho}H={\rho}E+p\\
u=\cfrac{\rho u}{\rho}=\mathbf{q}[1]/\mathbf{q}[0]\\
{\rho}Hu={\rho}Eu+pu=\mathbf{q}[2]\mathbf{q}[1]/\mathbf{q}[0]+p\mathbf{q}[1]/\mathbf{q}[0]\\
\end{array}\end{split}\]
\[\begin{split}\mathbf{q}=\begin{bmatrix}
\rho\\ \rho u \\\rho E
\end{bmatrix}=\begin{bmatrix}
\mathbf{q}[0]\\ \mathbf{q}[1] \\\mathbf{q}[2]
\end{bmatrix}\end{split}\]
\[\begin{split}\mathbf{F}=\begin{bmatrix}
\rho u\\ \rho uu + p \\\rho H u
\end{bmatrix}= \begin{bmatrix}
\mathbf{q}[1]\\ \cfrac{(\mathbf{q}[1])^{2}}{\mathbf{q}[0]}+p \\
\mathbf{q}[2]\mathbf{q}[1]/\mathbf{q}[0]+p\mathbf{q}[1]/\mathbf{q}[0]
\end{bmatrix}\end{split}\]
\[\begin{split}p=(\gamma-1)(\mathbf{q}[2]-\cfrac{1}{2} \cfrac{\mathbf{q}[1]^2}{\mathbf{q}[0]})\\\end{split}\]