cese.tex

% (C) 2022 Yung-Yu Chen.  All rights reserved.
% chktex-file 3
% chktex-file 8
% chktex-file 10
% chktex-file 13
% chktex-file 17
% chktex-file 25
% chktex-file 36
% chktex-file 37

\documentclass{turgon}

%\usepackage{lmodern}

\usepackage[printwatermark]{xwatermark}
\newwatermark[allpages,color=black!15,angle=55,scale=5,xpos=0,ypos=0]%
{DRAFT}

%\doublespacing
\linespread{1.2}

\title{
%
Code Development for the Space-Time Conservation Element and Solution Element
Method
%
}

\author{
%
Yung-Yu Chen
%
}

%\date{2008.6.4}

\begin{document}

\maketitle

\tableofcontents

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
\chapter*{Introduction}
\addcontentsline{toc}{chapter}{Introduction}
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Numerical computation, facilitated by the advancement of digital computers,
enables all kinds of simulation.  We use it to solve non-linear hyperbolic
partial differential equations (PDEs), which come from conservation laws
\citep{lax_hyperbolic_1973}.  This note describes a set of libraries that uses
the space-time conservation element and solution element (CESE) method to solve
the equations.  The libraries are written in modern C++ and provide Python
interface for interactive computing and prototyping.

% TODO: Basics of first-order hyperbolic partial differential equations:
% * d'Alembert solution
% * Method of characteristics
% * Non-linear equations
%   * Burger's equations
% * Riemann-invariant
%   * The shallow-water equations
%   * The Euler equations
% * Weak solutions

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
\chapter{The Conservation Element and Solution Element Method}
\label{c:cese}
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The conservation element and solution element (CESE) method solves conservation
laws, which can be written in the following form in one-dimensional space
\begin{align}
  \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0
  .
  \label{e:cese:1d_pde}
\end{align}
$u$ is the dependent solution variable and $f(u)$ is a function.  $(x, t)$ is
the independent variables defining the two axes of two-dimensional Euclidean
space.  Let $\mathbf{h} \defeq (f(u),u)$ and rewrite Eq.~(\ref{e:cese:1d_pde})
to $\nabla\cdot\mathbf{h} = 0$ with the divergence operator.  Assuming
Eq.~(\ref{e:cese:1d_pde}) applies everywhere in the control volume $V$, we
write
\begin{align*}
  \int_V\nabla\cdot\mathbf{h}\dif v = 0
  .
\end{align*}
By using the divergence theorem, the above differential equation is cast into
an integral equation over the control surface $S(V)$ surrounding $V$
\begin{align}
  \oint_{S(V)}\mathbf{h}\cdot\dif\hat{\sigma} = 0
  .  \label{e:cese:1d_integral_form}
\end{align}
$\dif\hat{\sigma}$ is the infinitesimal surface vector.  Equation
(\ref{e:cese:1d_integral_form}) does not require the point-wise divergence free
condition.

The CESE method solves the integral equation, instead of the differential
equation, by enforcing the space-time flux conservation as shown in
Eq.~(\ref{e:cese:1d_integral_form}) \citep{chang_method_1995}.  The solution
variable and its partial derivative in space are independent but solved
together.  It uses a compact stencil that defines two entities: the
conservation element (CE) and solution element (SE).  Space-time invariants are
used to minimize numerical dissipation, but the characteristics-based methods
are not used for obtaining solution.  Ad hoc treatments are avoided as much as
possible.  \cite{chang_new_1991} and \cite{chang_method_1995}.

\begin{figure}[htbp]
\centering
  \includegraphics{cce.eps}
  \caption{A compounded conservation element (CCE), the area enclosed by the
  {\color{red} red} dots of {\color{red} $\square\mathrm{BCEF}$}, contains two
  basic conservation elements (BCEs), the area enclosed by the {\color{blue}
  blue} dots of {\color{blue} $\square\mathrm{ABCD}$} and {\color{blue}
  $\square\mathrm{ADEF}$}.}
  \label{f:cce}
\end{figure}

The CEs discretize the space-time for the integral equation to be solved
(Eq.~\ref{e:cese:1d_integral_form})).  $\mathrm{CE}(j,n)$ denotes a single CE
associated with the grid point $(x_j, t^n)$.  In the CE, the conservation of
$\mathbf{h}$ is approximated as
\begin{align}
  \oint_{S(\mathrm{CE})}\mathbf{h}^*\cdot \dif\hat{\sigma} = 0
  \label{e:conserv_of_approx_h}
\end{align}
where $\mathbf{h}^*$, which will be defined later with SEs, denotes the
approximation of $\mathbf{h}$.  As shown in \figurename~\ref{f:cce}, a CE
defined like that is a compounded conservation element (CCE), consisting of two
adjacent basic conservation elements (BCEs) $\mathrm{CE}_-$ and
$\mathrm{CE}_+$.  Equation (\ref{e:conserv_of_approx_h}) holds in both CCEs and
BCEs, i.e.,
\begin{align*}
  \oint_{S(\mathrm{CE}_\pm)}\mathbf{h}^*\cdot \dif\hat{\sigma} = 0 .
\end{align*}
$S(\mathrm{CE}_{\pm})$ is the bounding surface surrounding $\mathrm{CE}_{\pm}$.

\begin{figure}[hbtp]
  \centering
  \includegraphics{cese_marching.eps}
  \caption{Time-marching the solution by using the cross-shaped solution
  elements.  The {\color{red} red} dotted crosses mark the SEs at $t=t^0$.  The
  {\color{blue} blue} dotted crosses mark the SEs at $t=t^{1/2}$.  The
  {\color{orange} orange} dotted crosses mark the SEs at $t=t^1$.  The bigger
  dots at the horizontal middle points of the crosses are the solution points.
  The arrows show how the solution variable and its spatial derivative, that
  are defined at the solution points, at the previous half time step propagate
  to those at the next half time step.  The solution points at the boundary,
  $x=x_0$ and $x=x_4$, need to be updated by boundary-condition treatments,
  rather than the CESE method scheme.}
  \label{f:cese_marching}
\end{figure}

The SEs determine $\mathbf{h}^*$.  There is more than one way to define SEs,
while an effective and consistent approach is shown in
\figurename~\ref{f:cese_marching}.  Let $\mathrm{SE}(j,n)$ denote the SE
associated with the grid point $(j,n)$, which is the cross-shaped mark enclosed
by the dotted line.  The solution variable approximation is written as
\begin{align*}
  u^*(x,t;j,n) = u_j^n + (u_x)_j^n(x-x_j) + (u_t)_j^n(t-t^n) .
\end{align*}
The grid point $(x_j, t^n)$ is used as the solution point.  $u_j^n$,
$(u_x)_j^n$, and $(u_t)_j^n$ hold constant in $\mathrm{SE}(j,n)$.  It should be
noted that every CE is surrounded by SEs.  Fluxes evaluated through the CE
boundary depends only on the approximation within SEs.  To proceed, write
\begin{align*}
  \frac{\partial u^*(x,t;j,n)}{\partial x} = (u_x)_j^n, \quad
  \frac{\partial u^*(x,t;j,n)}{\partial t} = (u_t)_j^n .
\end{align*}
Substitute the approximated solution variable $u^*$ back to the original
differential equation (Eq.~(\ref{e:cese:1d_pde})) and obtain the relation
between $(u_x)_j^n$ and $(u_t)_j^n$ as
\begin{align*}
                  (u_t)_j^n + (f_u)_j^n(u_x)_j^n = 0
  \;\Rightarrow\; (u_t)_j^n = -(f_u)_j^n(u_x)_j^n .
\end{align*}
The approximated solution variable $u^*$ is then rewritten as
\begin{align*}
  u^*(x,t;j,n) = u_j^n + (u_x)_j^n\left[(x-x_j) - (f_u)_j^n(t-t^n)\right] .
\end{align*}
Similarly, the approximated function
\begin{align*}
  f^*(x,t;j,n) = f_j^n + (f_x)_j^n(x-x_j) + (f_t)_j^n(t-t^n)
\end{align*}
is rewritten as
\begin{align*}
  f^*(x,t;j,n) = f_j^n + (f_u)_j^n (u_x)_j^n \left[
    (x-x_j) - (f_u)_j^n(t-t^n)
  \right] .
\end{align*}

To this point, we are ready to write the time-marching formulae for the
$c$-scheme.  It is the most simple time-marching scheme for the CESE method.
The formulae include updating the solution variable and the spatial derivative.
The part for the solution variable $u$ will be obtained by enforcing the
space-time flux conservation over the CCE as shown in Fig.~\ref{f:se_flux}.  It
should be noted that the height of the CEs and SEs is the half time step
$\Delta t/2$, not the full time step $\Delta t$.

\begin{figure}[hbtp]
  \centering
  \includegraphics{se_flux.eps}
  \caption{Space-time flux at the boundary of $\mathrm{CE}(j, n+\frac{1}{2})$
  defined by $\mathrm{SE}(j-\frac{1}{2}, n)$ ({\color{red}red}),
  $\mathrm{SE}(j+\frac{1}{2}, n)$ ({\color{blue} blue}), and $\mathrm{SE}(j,
  n+\frac{1}{2})$ ({\color{orange} orange}).  $x_j$ denotes the grid point of
  the $j$-th SE, $x_j^{\pm}$ the right and left end points, $x_j^s$ the
  solution point.  $\Delta x_j^+ \defeq x_j^+ - x_j$ and $\Delta x_j^- \defeq
  x_j - x_j^-$ are the length of the right and left arms of the $j$-th SE.
  $(\mathbf{h^*})_{j,\pm}^n$ and $(\mathbf{h^*})_{j}^{n,+}$ are the arithmetic
  average of the right, left, and upper arm of $\mathrm{SE}(j,n)$,
  respectively.}
  \label{f:se_flux}
\end{figure}

\begin{figure}[hbtp]
  \centering
  \includegraphics{nonuni_se.eps}
  \caption{The SE definition for a non-uniform one-dimensional grid.  The
  cross-shaped marks are the SEs, and the solid dots are the associated
  solution points.}
  \label{f:nonuni_se}
\end{figure}

The solution point must be at the center of the SE to make first-order
approximation consistent.  In a non-uniform grid, the solution points may not
collocate with grid points (see Fig.~\ref{f:nonuni_se}).  The approximation
formulae should be changed to
\begin{align*}
  u^*(x,t;j,n) &= u_j^n + (u_x)_j^n \left[
    (x-x_j^s) - (f_u)_j^n(t-t^n)
  \right] , \\
  f^*(x,t;j,n) &= f_j^n + (f_u)_j^n (u_x)_j^n \left[
    (x-x_j^s) - (f_u)_j^n(t-t^n)
  \right] .
\end{align*}
The formula for the solution variable $u$ is found to be
\begin{align}
  u_j^{n+\frac{1}{2}}
    = \frac{1}{\Delta x_j}\left\{
      (u^*)_{j-\frac{1}{2},+}^n \Delta x_{j-\frac{1}{2}}^+
    + (u^*)_{j+\frac{1}{2},-}^n \Delta x_{j+\frac{1}{2}}^-
    + \frac{\Delta t}{2} \left[
        (f^*)_{j-\frac{1}{2}}^{n,+}
      - (f^*)_{j+\frac{1}{2}}^{n,+}
      \right]
      \right\}
  \label{e:formula:u}
\end{align}
where
\begin{align*}
  (u^*)_{j\mp\frac{1}{2},\pm}^n
   &= u_{j\mp\frac{1}{2}}^n
    + (u_x)_{j\mp\frac{1}{2}}^n
      \left( x_{j\mp\frac{1}{2}}
           \pm \frac{1}{2} \Delta x_{j\mp\frac{1}{2}}^{\pm}
           - x_{j\mp\frac{1}{2}}^s \right), \\
  (f^*)_{j\mp\frac{1}{2}}^{n,\pm}
   &= f_{j\mp\frac{1}{2}}^n
    + (f_u)_{j\mp\frac{1}{2}}^n(u_x)_{j\mp\frac{1}{2}}^n
      \left[x_{j\mp\frac{1}{2}} - x_{j\mp\frac{1}{2}}^s
          - (f_u)_{j\mp\frac{1}{2}}^n \frac{\Delta t}{4}
      \right] .
\end{align*}
See Fig.~\ref{f:nonuni_se} for the definition of the left-hand side of the
above equations.

The first-order derivative $u_x$ needs another formula.  The one of the
$c$-scheme is
\begin{align}
  (u_x)_j^{n+\frac{1}{2}} = \frac{
    (u')_{j+\frac{1}{2}}^n - (u')_{j-\frac{1}{2}}^n
  }{\Delta x_j}
  \label{e:formula:ux:c}
\end{align}
where
\begin{align*}
  (u')_{j\pm\frac{1}{2}}^n &\defeq
      u_{j\pm\frac{1}{2}}^n
    + (u_x)_{j\pm\frac{1}{2}}^n \left[
        x_{j\pm\frac{1}{2}} - x_{j\pm\frac{1}{2}}^s
      - (f_u)_{j\pm\frac{1}{2}}^n \frac{\Delta t}{2} \right]
\end{align*}
are the Taylor expansion to $t^{n+\frac{1}{2}}$ with respect to
$\mathrm{SE}(j\pm\frac{1}{2}, n)$.  Equations (\ref{e:formula:u}) and
(\ref{e:formula:ux:c}) together form the $c$-scheme.

A weighting function should be introduced to treat discontinuity in space.
Equations (\ref{e:formula:u}) and (\ref{e:formula:ux:c}) are called the
$c$-scheme because Eq.~(\ref{e:formula:ux:c}) uses central-differencing to
approximate $(u_x)_j^{n+\frac{1}{2}}$.  It is second-order accurate but doesn't
give correct result when the solution variable $u$ is discontinuous between
$x_{j-\frac{1}{2}}$ and $x_{j+\frac{1}{2}}$.  To define the weighting function,
let
\begin{align*}
  (u_{x\pm})_j^{n+\frac{1}{2}} \defeq
    \frac{
      (u')_{j\pm\frac{1}{2}}^n - u_j^{n+\frac{1}{2}}
    }{x_{j\pm\frac{1}{2}} - x_j}
\end{align*}
be the spatial differences in the two intervals $[x_j, x_{j+\frac{1}{2}}]$ and
$[x_{j-\frac{1}{2}}, x_j]$, respectively.  The weighted average of the spatial
difference can then be calculated with
\begin{align*}
  (u_x^w)_j^{n+\frac{1}{2}} =
  W\left(
    (u_{x-})_j^{n+\frac{1}{2}}, (u_{x+})_j^{n+\frac{1}{2}}
  \right)
\end{align*}
where $W$ is the weighting function.  When there is discontinuity in $u$
between $x_{j-\frac{1}{2}}$ and $x_{j+\frac{1}{2}}$, the weighting function
should detect it and approximate the spatial differencing properly.  A choice
of the weighting function is
\begin{align}
  W_{\alpha} =
    \frac{|u_{x+}|^{\alpha}u_{x-} + |u_{x-}|^{\alpha}u_{x+}}
         {|u_{x-}|^{\alpha} + |u_{x+}|^{\alpha}},
  \label{e:formula:wa}
\end{align}
where $\alpha \in \mathbb{R}$ is a constant parameter and $|u_{x-}|^{\alpha} +
|u_{x+}|^{\alpha} > 0$.  $\alpha$ is often picked as a positive integer, e.g.,
2, for saving computation cycles.  For non-linear equations or discontinuous
initial conditions, a weighting function, e.g., Eq.~(\ref{e:formula:wa}), is
necessary to keep the solution from diverging.

\section{CFL Insensitive Scheme}
\label{s:cese:ctau}

\section{Local Time-Stepping}
\label{s:cese:lts}

\section{Boundary-Condition Treatment}
\label{s:cese:bc}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
\chapter{Solving Systems of Equations}
\label{c:syseq}
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{The Euler Equations}
\label{s:syseq:euler}

\subsection{Governing Equations of Fluid Flows}

The governing equations of fluid flow consist of the conservation of mass,
momentum, and energy
%
\begin{gather}
  \frac{\partial\rho}{\partial t} + \frac{\partial\rho v_j}{\partial x_j}
  = 0
  \label{e:euler:mass}
  \\
  \frac{\partial\rho v_i}{\partial t}
  + \frac{\partial\rho v_iv_j}{\partial x_j}
  = \frac{\partial p}{\partial x_j} + \rho b_i
  \label{e:euler:momentum}
  \\
  \frac{\partial}{\partial t}
  \left[\rho\left( e + \frac{v_k^2}{2} \right)\right]
  + \frac{\partial}{\partial x_j}
  \left[\rho\left( e + \frac{v_k^2}{2} \right)v_j\right]
  = \rho \dot{q} - \frac{\partial pv_j}{\partial x_j} + \rho b_jv_j
  \label{e:euler:energy}
\end{gather}
%
Einstein's index summation convention is used.
The unknowns are density $\rho$, velocity $\mathbf{v}$, pressure $p$, and
internal energy $e$.
Body force $\mathbf{b}$ and heat generation $\dot{q}$ are given.

There are 5 equations and 6 unknowns ($\rho$, $\mathbf{v}$, $p$, and $e$).
The equation of state is used to close the system of equation:
%
\begin{align}
  p = \rho RT \label{e:euler:eos}
\end{align}
where $R$ is the ideal gas constant and $T$ the temperature.
%
The temperature $T$ is related to the internal energy $e$ by
%
\begin{align}
  e = c_vT = \frac{RT}{\gamma-1} = \frac{1}{\gamma-1}\frac{p}{\rho}
  \label{e:euler:internal_energy}
\end{align}
%
where $c_v$ is the specific heat at constant volume and $\gamma$ the specific
heat ratio.

The 5 governing equations (\ref{e:euler:mass}) (\ref{e:euler:momentum})
(\ref{e:euler:energy}) are closed by using the two additional equations
(\ref{e:euler:eos}) and (\ref{e:euler:internal_energy}) and the additional
variable $T$.

The conservation variables:
%
\begin{align}
  \mathbf{u} \defeq \left(
  \begin{array}{c}
    u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5
  \end{array}\right) = \left(
  \begin{array}{c}
    \rho \\ \rho v_1 \\ \rho v_2 \\ \rho v_3 \\
    \rho\left(e+\frac{v_k^2}{2}\right)
  \end{array}\right)
  \label{e:euler:unknown}
\end{align}
%
Pressure is an important quantity and the representation by the conservation
variables is (subscripts are expanded explicitly):
\begin{align*}
  p = (\gamma-1)\left(u_5 - \frac{u_2^2+u_3^2+u_4^2}{2u_1}\right)
\end{align*}
%
Rewrite the 5 governing equations by using the conservation variables:
%
\begin{align}
  \frac{\partial u_1}{\partial t}
  + \frac{\partial u_2}{\partial x_1}
  + \frac{\partial u_3}{\partial x_2}
  + \frac{\partial u_4}{\partial x_3} = 0
  \label{e:euler:gov1}
\end{align}
%
\begin{align}
  \begin{aligned} &\frac{\partial u_2}{\partial t}
  + \frac{\partial}{\partial x_1}\left(\frac{u_2^2}{u_1}\right)
  + \frac{\partial}{\partial x_2}\left(\frac{u_2u_3}{u_1}\right)
  + \frac{\partial}{\partial x_3}\left(\frac{u_2u_4}{u_1}\right) = \\
  &\quad -\frac{\partial}{\partial x_1}\left[
    (\gamma-1)\left(u_5 - \frac{u_2^2+u_3^2+u_4^2}{2u_1}\right)
    \right] + b_1u_1
  \end{aligned}
  \label{e:euler:gov2}
\end{align}
%
\begin{align}
  \begin{aligned} &\frac{\partial u_3}{\partial t}
  + \frac{\partial}{\partial x_1}\left(\frac{u_2u_3}{u_1}\right)
  + \frac{\partial}{\partial x_2}\left(\frac{u_3^2}{u_1}\right)
  + \frac{\partial}{\partial x_3}\left(\frac{u_3u_4}{u_1}\right) = \\
  &\quad -\frac{\partial}{\partial x_2}\left[
    (\gamma-1)\left(u_5 - \frac{u_2^2+u_3^2+u_4^2}{2u_1}\right)
    \right] + b_2u_1
  \end{aligned}
  \label{e:euler:gov3}
\end{align}
%
\begin{align}
  \begin{aligned} &\frac{\partial u_4}{\partial t}
  + \frac{\partial}{\partial x_1}\left(\frac{u_2u_4}{u_1}\right)
  + \frac{\partial}{\partial x_2}\left(\frac{u_3u_4}{u_1}\right)
  + \frac{\partial}{\partial x_3}\left(\frac{u_4^2}{u_1}\right) = \\
  &\quad -\frac{\partial}{\partial x_3}\left[
    (\gamma-1)\left(u_5 - \frac{u_2^2+u_3^2+u_4^2}{2u_1}\right)
    \right] + b_3u_1
  \end{aligned}
  \label{e:euler:gov4}
\end{align}
%
\begin{align}
  \begin{aligned} &\frac{\partial u_5}{\partial t}
  + \frac{\partial}{\partial x_1}\left(\frac{u_2u_5}{u_1}\right)
  + \frac{\partial}{\partial x_2}\left(\frac{u_3u_5}{u_1}\right)
  + \frac{\partial}{\partial x_3}\left(\frac{u_4u_5}{u_1}\right) = \\
  &\quad - \frac{\partial}{\partial x_1}\left[
    (\gamma-1)\left(u_5 - \frac{u_2^2+u_3^2+u_4^2}{2u_1}\right)
    \frac{u_2}{u_1}
    \right] \\
  &\quad - \frac{\partial}{\partial x_2}\left[
    (\gamma-1)\left(u_5 - \frac{u_2^2+u_3^2+u_4^2}{2u_1}\right)
    \frac{u_3}{u_1}
    \right] \\
  &\quad - \frac{\partial}{\partial x_3}\left[
    (\gamma-1)\left(u_5 - \frac{u_2^2+u_3^2+u_4^2}{2u_1}\right)
    \frac{u_4}{u_1}
    \right]
  + \rho\dot{q} + b_1u_2 + b_2u_3 + b_3u_4
  \end{aligned}
  \label{e:euler:gov5}
\end{align}

\subsection{Vector Flux Functions}

Reorganize the above 5 equations (\ref{e:euler:gov1}) (\ref{e:euler:gov2})
(\ref{e:euler:gov3}) (\ref{e:euler:gov4}) (\ref{e:euler:gov5}) into a vector
form
%
\begin{align}
  \frac{\partial\mathbf{u}}{\partial t}
  + \sum_{\mu=1}^3 \frac{\partial\mathbf{f}^{(\mu)}}{\partial x_{\mu}}
  = \mathbf{s}
  \label{e:euler:vec}
\end{align}
%
the symbol $\mathbf{s}$ at the right-hand side is the lumped source term.
There are 5 equations in Eq.~(\ref{e:euler:vec}).
The flux function $\mathbf{f}^{(1)}$ is defined as
%
\begin{align}
  \begin{aligned}
    f^{(1)}_1 &= u_2 \\
    f^{(1)}_2 &= (\gamma-1)u_5 - \frac{\gamma-3}{2}\frac{u_2^2}{u_1}
    - \frac{\gamma-1}{2}\frac{u_3^2}{u_1}
    - \frac{\gamma-1}{2}\frac{u_4^2}{u_1} \\
    f^{(1)}_3 &= \frac{u_2u_3}{u_1} \\
    f^{(1)}_4 &= \frac{u_2u_4}{u_1} \\
    f^{(1)}_5 &= \gamma\frac{u_2u_5}{u_1}
    - \frac{\gamma-1}{2}\frac{u_2^2+u_3^2+u_4^2}{u_1}\frac{u_2}{u_1}
  \end{aligned}
  \label{e:euler:flux1}
\end{align}
%
$\mathbf{f}^{(2)}$ as
%
\begin{align}
  \begin{aligned}
    f^{(2)}_1 &= u_3 \\
    f^{(2)}_2 &= \frac{u_2 u_3}{u_1} \\
    f^{(2)}_3 &= (\gamma-1)u_5 - \frac{\gamma-1}{2}\frac{u_2^2}{u_1}
    - \frac{\gamma-3}{2}\frac{u_3^2}{u_1}
    - \frac{\gamma-1}{2}\frac{u_4^2}{u_1} \\
    f^{(2)}_4 &= \frac{u_3 u_4}{u_1} \\
    f^{(2)}_5 &= \gamma\frac{u_3u_5}{u_1}
    - \frac{\gamma-1}{2}\frac{u_2^2+u_3^2+u_4^2}{u_1}\frac{u_3}{u_1}
  \end{aligned}
  \label{e:euler:flux2}
\end{align}
%
$\mathbf{f}^{(3)}$ as
%
\begin{align}
  \begin{aligned}
    f^{(3)}_1 &= u_4 \\
    f^{(3)}_2 &= \frac{u_2 u_4}{u_1} \\
    f^{(3)}_3 &= \frac{u_3 u_4}{u_1} \\
    f^{(3)}_4 &= (\gamma-1)u_5 - \frac{\gamma-1}{2}\frac{u_2^2}{u_1}
    - \frac{\gamma-1}{2}\frac{u_3^2}{u_1}
    - \frac{\gamma-3}{2}\frac{u_4^2}{u_1} \\
    f^{(3)}_5 &= \gamma\frac{u_4u_5}{u_1}
    - \frac{\gamma-1}{2}\frac{u_2^2+u_3^2+u_4^2}{u_1}\frac{u_4}{u_1}
  \end{aligned}
  \label{e:euler:flux3}
\end{align}
%
The lumped source term $\mathbf{s}$ is
\begin{align}
  \begin{aligned}
    s_1 &= 0 \\
    s_2 &= b_1 u_1 \\
    s_3 &= b_2 u_1 \\
    s_4 &= b_3 u_3 \\
    s_5 &= \dot{q}u_1 + b_1 u_2 + b_2 u_3 + b_3 u_4
  \end{aligned}
  \label{e:euler:sterm}
\end{align}

\subsection{Quasi-Linear System Equations}

Expand Eq.~(\ref{e:euler:vec}) to an index form:
%
\begin{align}
  \frac{\partial u_m}{\partial t}
  + \sum_{\mu=1}^3 \frac{\partial f^{(\mu)}_m}{\partial x_{\mu}}
  = s_m, \quad m = 1, \ldots, 5
  \label{e:euler:idx_full}
\end{align}
%
The Euler equations are inviscid.
The source term on the right-hand side of the above equation
(\ref{e:euler:idx_full}) is dropped for the Euler equations:
%
\begin{align}
  \frac{\partial u_m}{\partial t}
  + \sum_{\mu=1}^3 \frac{\partial f^{(\mu)}_m}{\partial x_{\mu}}
  = 0, \quad m = 1, \ldots, 5
  \label{e:euler:idx}
\end{align}
%
Aided by the notation
%
\begin{align*}
  \begin{aligned}
    u_{mt} &\defeq \frac{\partial u_m}{\partial t} \\
    u_{mx_{\mu}} &\defeq \frac{\partial u_m}{\partial x_{\mu}} \\
    f^{(\mu)}_{m,l} &\defeq \frac{\partial f^{(\mu)}_m}{\partial u_l}
  \end{aligned},
  \quad m, l = 1, 2, 3, 4, 5
\end{align*}
%
the Euler equations can be written in the quasi-linear form:
%
\begin{align}
  \frac{\partial\mathbf{u}}{\partial t} + \sum_{\mu=1}^3
  \mathrm{A}^{(\mu)} \frac{\partial\mathbf{u}}{\partial x_{\mu}} = 0
  \label{e:euler:qlinear}
\end{align}
%
where
\begin{align*}
  \mathrm{A}^{(\mu)} = \left[f^{(\mu)}_{m,l}\right],
  \quad \mu = 1, 2, 3 \;\mathrm{and}\; m, l = 1, 2, 3, 4, 5
\end{align*}

List the elements of $\mathrm{A^{(1)}}$, $\mathrm{A^{(2)}}$, and
$\mathrm{A^{(3)}}$ as
%
\begin{gather}
  \begin{gathered}
    \mathrm{A}^{(1)} = \left(
    \begin{array}{ccccc}
      0             & 1             & 0             & 0             & 0 \\
      f^{(1)}_{2,1} & f^{(1)}_{2,2} & f^{(1)}_{2,3} & f^{(1)}_{2,4} &
      \gamma - 1 \\
      f^{(1)}_{3,1} & f^{(1)}_{3,2} & f^{(1)}_{3,3} & 0             & 0 \\
      f^{(1)}_{4,1} & f^{(1)}_{4,2} & 0             & f^{(1)}_{4,4} & 0 \\
      f^{(1)}_{5,1} & f^{(1)}_{5,2} & f^{(1)}_{5,3} & f^{(1)}_{5,4} &
      f^{(1)}_{5,5}
    \end{array}
    \right)
    \\
    \begin{aligned}
      f^{(1)}_{2,1} &= \frac{\gamma-3}{2}\frac{u_2^2}{u_1^2}
      + \frac{\gamma-1}{2}\frac{u_3^2}{u_1^2}
      + \frac{\gamma-1}{2}\frac{u_4^2}{u_1^2}, \\
      f^{(1)}_{2,2} &= -(\gamma-3)\frac{u_2}{u_1}, \quad
      f^{(1)}_{2,3} = -(\gamma-1)\frac{u_3}{u_1}, \quad
      f^{(1)}_{2,4} = -(\gamma-1)\frac{u_4}{u_1}, \\
      f^{(1)}_{3,1} &= -\frac{u_2 u_3}{u_1^2}, \quad
      f^{(1)}_{3,2} = \frac{u_3}{u_1}, \quad
      f^{(1)}_{3,3} = f^{(1)}_{4,4} = \frac{u_2}{u_1}, \\
      f^{(1)}_{4,1} &= -\frac{u_2 u_4}{u_1^2}, \quad
      f^{(1)}_{4,2} = \frac{u_4}{u_1}, \quad
      f^{(1)}_{4,4} = f^{(1)}_{3,3} = \frac{u_2}{u_1}, \\
      f^{(1)}_{5,1} &= -\gamma\frac{u_2 u_5}{u_1^2}
      + (\gamma-1)\frac{u_2^2+u_3^2+u_4^2}{u_1^2}\frac{u_2}{u_1}, \\
      f^{(1)}_{5,2} &= \gamma\frac{u_5}{u_1}
      - \frac{\gamma-1}{2}\frac{3u_2^2 + u_3^2 + u_4^2}{u_1^2}, \\
      f^{(1)}_{5,3} &= -(\gamma-1)\frac{u_2 u_3}{u_1^2}, \quad
      f^{(1)}_{5,4} = -(\gamma-1)\frac{u_2 u_4}{u_1^2}, \quad
      f^{(1)}_{5,5} = \gamma\frac{u_2}{u_1}
    \end{aligned}
  \end{gathered}
  \label{e:euler:jaco1}
\end{gather}
%
\begin{gather}
  \begin{gathered}
    \mathrm{A}^{(2)} = \left(
    \begin{array}{ccccc}
      0             & 0             & 1             & 0             & 0 \\
      f^{(2)}_{2,1} & f^{(2)}_{2,2} & f^{(2)}_{2,3} & 0             & 0 \\
      f^{(2)}_{3,1} & f^{(2)}_{3,2} & f^{(2)}_{3,3} & f^{(2)}_{3,4} &
      \gamma - 1 \\
      f^{(2)}_{4,1} & 0             & f^{(2)}_{4,3} & f^{(2)}_{4,4} & 0 \\
      f^{(2)}_{5,1} & f^{(2)}_{5,2} & f^{(2)}_{5,3} & f^{(2)}_{5,4} &
      f^{(1)}_{5,5}
    \end{array}
    \right)
    \\
    \begin{aligned}
      f^{(2)}_{2,1} &= -\frac{u_2 u_3}{u_1^2}, \quad
      f^{(2)}_{2,2} = f^{(2)}_{4,4} = \frac{u_3}{u_1}, \quad
      f^{(2)}_{2,3} = \frac{u_2}{u_1}, \\
      f^{(2)}_{3,1} &= \frac{\gamma-1}{2}\frac{u_2^2}{u_1^2}
      + \frac{\gamma-3}{2}\frac{u_3^2}{u_1^2}
      + \frac{\gamma-1}{2}\frac{u_4^2}{u_1^2}, \\
      f^{(2)}_{3,2} &= -(\gamma-1)\frac{u_2}{u_1}, \quad
      f^{(2)}_{3,3} = -(\gamma-3)\frac{u_3}{u_1}, \quad
      f^{(2)}_{3,4} = -(\gamma-1)\frac{u_4}{u_1}, \\
      f^{(2)}_{4,1} &= -\frac{u_3 u_4}{u_1^2}, \quad
      f^{(2)}_{4,3} = \frac{u_4}{u_1}, \quad
      f^{(2)}_{4,4} = f^{(2)}_{2,2} = \frac{u_3}{u_1}, \\
      f^{(2)}_{5,1} &= -\gamma\frac{u_3 u_5}{u_1^2}
      + (\gamma-1)\frac{u_2^2+u_3^2+u_4^2}{u_1^2}\frac{u_3}{u_1}, \\
      f^{(2)}_{5,2} &= -(\gamma-1)\frac{u_2 u_3}{u_1^2}, \\
      f^{(2)}_{5,3} &= \gamma\frac{u_5}{u_1}
      - \frac{\gamma-1}{2}\frac{u_2^2 + 3u_3^2 + u_4^2}{u_1^2}, \\
      f^{(2)}_{5,4} &= -(\gamma-1)\frac{u_3 u_4}{u_1^2}, \quad
      f^{(2)}_{5,5} = \gamma\frac{u_3}{u_1}
    \end{aligned}
  \end{gathered}
  \label{e:euler:jaco2}
\end{gather}
%
\begin{gather}
  \begin{gathered}
    \mathrm{A}^{(3)} = \left(
    \begin{array}{ccccc}
      0             & 0             & 0             & 1             & 0 \\
      f^{(3)}_{2,1} & f^{(3)}_{2,2} & 0             & f^{(3)}_{2,4} & 0 \\
      f^{(3)}_{3,1} & 0             & f^{(3)}_{3,3} & f^{(3)}_{3,4} & 0 \\
      f^{(3)}_{4,1} & f^{(3)}_{4,2} & f^{(3)}_{4,3} & f^{(3)}_{4,4} &
      \gamma - 1 \\
      f^{(3)}_{5,1} & f^{(3)}_{5,2} & f^{(3)}_{5,3} & f^{(3)}_{5,4} &
      f^{(3)}_{5,5}
    \end{array}
    \right)
    \\
    \begin{aligned}
      f^{(3)}_{2,1} &= -\frac{u_2 u_4}{u_1^2}, \quad
      f^{(3)}_{2,2} = f^{(3)}_{3,3} = \frac{u_4}{u_1}, \quad
      f^{(3)}_{2,4} = \frac{u_2}{u_1}, \\
      f^{(3)}_{3,1} &= -\frac{u_3 u_4}{u_1^2}, \quad
      f^{(3)}_{3,3} = f^{(3)}_{2,2} = \frac{u_4}{u_1}, \quad
      f^{(3)}_{3,4} = \frac{u_3}{u_1}, \\
      f^{(3)}_{4,1} &= \frac{\gamma-1}{2}\frac{u_2^2}{u_1^2}
      + \frac{\gamma-1}{2}\frac{u_3^2}{u_1^2}
      + \frac{\gamma-3}{2}\frac{u_4^2}{u_1^2}, \\
      f^{(3)}_{4,2} &= -(\gamma-1)\frac{u_2}{u_1}, \quad
      f^{(3)}_{4,3} = -(\gamma-1)\frac{u_3}{u_1}, \quad
      f^{(3)}_{4,4} = -(\gamma-3)\frac{u_4}{u_1}, \\
      f^{(3)}_{5,1} &= -\gamma\frac{u_4 u_5}{u_1^2}
      + (\gamma-1)\frac{u_2^2+u_3^2+u_4^2}{u_1^2}\frac{u_4}{u_1}, \quad
      f^{(3)}_{5,4} = \gamma\frac{u_5}{u_1}
      - \frac{\gamma-1}{2}\frac{u_2^2 + u_3^2 + 3u_4^2}{u_1^2}, \\
      f^{(3)}_{5,2} &= -(\gamma-1)\frac{u_2 u_4}{u_1^2}, \quad
      f^{(3)}_{5,3} = -(\gamma-1)\frac{u_3 u_4}{u_1^2}, \quad
      f^{(3)}_{5,5} = \gamma\frac{u_4}{u_1}
    \end{aligned}
  \end{gathered}
  \label{e:euler:jaco3}
\end{gather}

Equation~(\ref{e:euler:idx}) and the chain rule provide the following
properties:
%
\begin{gather}
  u_{mt} =
  -\sum_{\mu=1}^3
  \left(
  \sum_{l=1}^5 f^{(\mu)}_{m,l} u_{lx_{\mu}}
  \right),
  \quad m = 1, 2, 3, 4, 5
  \label{e:euler:ut}
  \\
  f^{(\mu)}_{mt} =
  \sum_{l=1}^5 f^{(\mu)}_{m,l} u_{lt},
  \quad \mu = 1, 2, 3, m = 1, 2, 3, 4, 5
  \label{e:euler:ft}
\end{gather}

\subsection{One-Dimension Euler Equations}

Reduce the governing equations (\ref{e:euler:mass}), (\ref{e:euler:momentum}),
and(\ref{e:euler:energy}) to one dimensional:
%
\begin{gather}
  \frac{\partial\rho}{\partial t} + \frac{\partial\rho v}{\partial x}
  = 0
  \label{e:euler1d:mass}
  \\
  \frac{\partial\rho v}{\partial t}
  + \frac{\partial\rho v^2}{\partial x}
  = \frac{\partial p}{\partial x} + \rho b
  \label{e:euler1d:momentum}
  \\
  \frac{\partial}{\partial t}
  \left[\rho\left( e + \frac{v^2}{2} \right)\right]
  + \frac{\partial}{\partial x}
  \left[\rho\left( e + \frac{v^2}{2} \right)v\right]
  = \rho \dot{q} - \frac{\partial pv}{\partial x} + \rho b v
  \label{e:euler1d:energy}
\end{gather}
%
Equations~(\ref{e:euler:eos}) and (\ref{e:euler:internal_energy}) are also used
to close the equations.

The conservation variables:
%
\begin{align}
  \mathbf{u} \defeq \left(
  \begin{array}{c}
    u_1 \\ u_2 \\ u_3
  \end{array}\right) = \left(
  \begin{array}{c}
    \rho \\ \rho v \\ \rho\left(e+\frac{v^2}{2}\right)
  \end{array}\right)
  \label{e:euler1d:unknown}
\end{align}
%
Important physical quantities
%
\begin{gather*}
  \rho = u_1 \\
  v = \frac{u_2}{u_1} \\
  p = (\gamma - 1)(u_3 - \frac{u_2^2}{2u_1}) \\
  T = \frac{\gamma - 1}{R}
  \left(
  \frac{u_3}{u_1} - \frac{1}{2} \frac{u_2^2}{u_1^2}
  \right)
\end{gather*}
%
\begin{gather}
  \frac{\partial\mathbf{u}}{\partial t} + \frac{\partial\mathbf{f}}{\partial x}
  = 0
  \label{e:euler1d:vec}
\end{gather}
%
The flux function is
%
\begin{align}
  \begin{aligned}
    f_1 &= u_2 \\
    f_2 &= (\gamma-1)u_3 - \frac{\gamma-3}{2}\frac{u_2^2}{u_1} \\
    f_3 &= \gamma\frac{u_2 u_3}{u_1}
    - \frac{\gamma-1}{2}\frac{u_2^3}{u_1^2}
  \end{aligned}
  \label{e:euler1d:flux}
\end{align}

The index form is expanded from Eq.~(\ref{e:euler1d:vec})
%
\begin{align}
  \frac{\partial u_m}{\partial t} + \frac{\partial f_m}{\partial x}
  = 0, \quad m = 1, 2, 3
  \label{e:euler1:idx}
\end{align}
%
The properties derived from the index form
%
\begin{gather}
  u_{mt} = -\sum_{l=1}^3 f_{m,l} u_{lx},  \quad m = 1, 2, 3
  \label{e:euler1:ut}
  \\
  f_{mt} = \sum_{l=1}^3 f_{m,l} u_{lt},  \quad m = 1, 2, 3
  \label{e:euler1:ft}
\end{gather}
%
The quasi-linear form
%
\begin{align}
  \frac{\partial\mathbf{u}}{\partial t}
  + \mathrm{A} \frac{\partial\mathbf{u}}{\partial x} = 0
  \label{e:euler1d:qlinear}
\end{align}
%
has the elements in the matrix
%
\begin{gather}
  \begin{gathered}
    \mathrm{A} = \left(
    \begin{array}{ccc}
      0       & 1       & 0          \\
      f_{2,1} & f_{2,2} & \gamma - 1 \\
      f_{3,1} & f_{3,2} & f_{3,3}
    \end{array}
    \right)
    \\
    \begin{aligned}
      f_{2,1} &= \frac{\gamma-3}{2}\frac{u_2^2}{u_1^2}, \quad
      f_{2,2} = -(\gamma-3)\frac{u_2}{u_1}, \\
      f_{3,1} &= -\gamma\frac{u_2 u_3}{u_1^2}
      + (\gamma-1)\frac{u_2^3}{u_1^3}, \\
      f_{3,2} &= \gamma\frac{u_3}{u_1}
      - \frac{3}{2}(\gamma-1)\frac{u_2^2}{u_1^2}, \\
      f_{3,3} &= \gamma\frac{u_2}{u_1}
    \end{aligned}
  \end{gathered}
  \label{e:euler1d:jaco}
\end{gather}

\subsection{Shock Tube Problem}

Consider a (one-dimensional) tube filled with two gases that are separated by a
diaphragm (\figurename~\ref{f:tube_at_rest}).
A high-pressure gas is at the left-hand side, and a low-pressure gas is at
the right-hand side.
$p$ denotes the pressure, $\rho$ the mass density, $\gamma$ the ratio of
specific heat, and $v$ the velocity.
The gases is at rest initially ($t = t_0$).
%
\begin{align}
  p_1 > p_5 , \quad
  \rho_1 > \rho_5, \quad
  v_1 = v_5
\end{align}
%
The gas at the high-pressure side is called the driver gas, while the gas at
the low-pressure side is called the driven gas.
When the diaphragm is removed, the driver gas pushes toward the driven gas
and the gases around the diaphragm starts to move to right.
See \figurename~\ref{f:tube_move_right}.

\begin{figure}[h]
  \centering
  \includegraphics{tube_at_rest.eps}
  \caption{Gases are at rest in the tube.}
  \label{f:tube_at_rest}
  \includegraphics{tube_move_right.eps}
  \caption{Gases move to right after the diaphragm rupture.}
  \label{f:tube_move_right}
\end{figure}

\begin{figure}[h]
  \centering
  \includegraphics{tube_zones.eps}
  \caption{The zones of flow after the diaphragm disrupts in the tube.}
  \label{f:tube_zones}
\end{figure}

The rupture of the diaphragm generates a right-moving normal shock wave and a
left-moving expansion wave.
The flow in the tube looks like \figurename~\ref{f:tube_zones}.
The flow properties in zones 1 and 5 remain unchanged.
The expansion wave is in zone 2.
The contact surface is between zones 3 and 4.
The entropy across the contact surface is discontinuous, but the condition
$v_3 = v_4 = v_f$ and $p_3 = p_4$ hold.
The right-moving normal shock wave is between zones 4 and 5.
Analytical solution of the problem can be obtained by solving the problem of the
moving shock wave and the expansion wave~\citep{anderson_modern_2003}.

\subsubsection{Moving Normal Shock}

Across zones 4 and 5 there is a normal shock moving at the speed $v_s$.
The gas velocity $v_5$ in zone 5 is 0.
The conservation of mass, momentum, and energy across the normal shock moving
at the speed $v_s$ are written as
%
\begin{align}
  \rho_4 (v_s - v_4) &= \rho_5 v_s
  \label{e:nshock:mass}
  \\
  p_4 + \rho_4 (v_s - v_4)^2 &= p_5 + \rho_5 v_s^2
  \label{e:nshock:momentum}
  \\
  h_4 + \frac{(v_s - v_4)^2}{2} &= h_5 + \frac{v_s^2}{2}
  \label{e:nshock:energe}
\end{align}
%
Obtain the relationship between $v_s$ and $v_s-v_4$ by using
Eq.~(\ref{e:nshock:mass}) (conservation of mass)
%
\begin{gather*}
  v_s = \frac{\rho_4}{\rho_5}(v_s - v_4) \\
  v_s - v_4 = \frac{\rho_5}{\rho_4}v_s
\end{gather*}
%
Substitute the equation above to Eq.(\ref{e:nshock:momentum}) to have
%
\begin{gather*}
  p_4 + \rho_4 \left(\frac{\rho_5}{\rho_4}v_s\right)^2 = p_5 + \rho_5 v_s^2 \\
  \Rightarrow \quad
  p_4 - p_5 = \rho_5 v_s^2 - \frac{\rho_5^2}{\rho_4} v_s^2
  = \rho_5 v_s^2 \left(1 - \frac{\rho_5}{\rho_4}\right) \\
  \Rightarrow \quad
  v_s^2 = \frac{p_4 - p_5}{\rho_5\left(1 - \frac{\rho_5}{\rho_4}\right)}
  = \left(\frac{\rho_4}{\rho_5}\right)
  \left(\frac{p_4 - p_5}{\rho_4 - \rho_5}\right) \\
  \Rightarrow \quad
  (v_s - v_4)^2 = \left(\frac{\rho_5}{\rho_4}\right)^2 v_s^2
  = \left(\frac{\rho_5}{\rho_4}\right)
  \left(\frac{p_4 - p_5}{\rho_4 - \rho_5}\right) \\
\end{gather*}
%
Recall the relationship $h = e + \frac{p}{\rho}$ and use the above expressions
of $v_s^2$ and $(v_s - v_4)^2$.
Substitute them to Eq.~(\ref{e:nshock:mass}) to have
%
\begin{gather*}
  e_4 + \frac{p_4}{\rho_4}
  + \frac{1}{2} \left(\frac{\rho_5}{\rho_4}\right)
  \left(\frac{p_4 - p_5}{\rho_4 - \rho_5}\right)
  = e_5 + \frac{p_5}{\rho_5}
  + \left(\frac{\rho_4}{\rho_5}\right)
  \left(\frac{p_4 - p_5}{\rho_4 - \rho_5}\right)
  \\
  \Rightarrow
  e_4 - e_5 = \frac{p_5}{\rho_5} - \frac{p_4}{\rho_4}
  + \frac{1}{2} \left(\frac{\rho_4}{\rho_5} - \frac{\rho_5}{\rho_4}\right)
  \left(\frac{p_4 - p_5}{\rho_4 - \rho_5}\right)
\end{gather*}
%
Rewrite to obtain
%
\begin{align*}
  e_4 - e_5 &= \frac{2 (\rho_4 p_5 - \rho_5 p_4) (\rho_4 - \rho_5)}
  {2 \rho_4 \rho_5 (\rho_4 - \rho_5)}
  + \frac{(\rho_4^2 - \rho_5^2) (p_4 - p_5)}
  {2 \rho_4 \rho_5 (\rho_4 - \rho_5)}
  \\
  &= \left(\frac{1}{2 \rho_4 \rho_5}\right)
  \left(
  \cancel{2} \rho_4 p_5 - \cancel{2} \rho_5 p_4
  + \rho_4 p_4 - \cancel{\rho_4 p_5} + \cancel{\rho_5 p_4} - \rho_5 p_5
  \right)
  \\
  &= \left(\frac{p_4 + p_5}{2}\right)
  \left(\frac{\rho_4 - \rho_5}{\rho_4 \rho_5}\right)
\end{align*}
%
The Hugoniot equation is obtained for the moving shock
%
\begin{align}
  e_4 - e_5 = \left(\frac{p_4 + p_5}{2}\right)
  \left(\frac{1}{\rho_5} - \frac{1}{\rho_4}\right)
  \label{e:nshock:hugoniot}
\end{align}

By assuming calorically ideal gas, $e = c_v T$, $\gamma = \frac{R}{c_v} + 1$,
and $p = \rho RT$, the Hugoniot equation leads to
%
\begin{gather}
  \frac{T_4}{T_5} = \frac{p_4}{p_5}
  \left(
  \frac{\dfrac{\gamma + 1}{\gamma - 1} + \dfrac{p_4}{p_5}}
  {1 + \dfrac{\gamma + 1}{\gamma - 1} \dfrac{p_4}{p_5}}
  \right)
  \label{e:nshock:T}
  \\
  \frac{\rho_4}{\rho_5} =
  \left(
  \frac{1 + \dfrac{\gamma + 1}{\gamma - 1} \dfrac{p_4}{p_5}}
  {\dfrac{\gamma + 1}{\gamma - 1} + \dfrac{p_4}{p_5}}
  \right)
  \label{e:nshock:rho}
\end{gather}
%
Define the Mach number of the moving normal shock
%
\begin{gather*}
  M_s = \frac{v_s}{a_5}
\end{gather*}
%
The relation of the pressure across the shock wave is
%
\begin{gather*}
  \frac{p_4}{p_5} = 1 + \frac{2\gamma}{\gamma + 1}(M_s^2 - 1) \\
  \Rightarrow \quad
  M_s = \sqrt{\frac{\gamma + 1}{2\gamma}\left(\frac{p_4}{p_5} - 1\right) + 1}
\end{gather*}

Aided by the definition of $M_s$, the velocity of the shock can be written as
%
\begin{gather}
  v_s = a_5
  \sqrt{\frac{\gamma + 1}{2\gamma}\left(\frac{p_4}{p_5} - 1\right) + 1}
  \label{e:nshock:v_s}
\end{gather}
%
Plug Eq.~(\ref{e:nshock:rho}) and the above equation (\ref{e:nshock:v_s}) to
Eq.~(\ref{e:nshock:mass}) (conservation of mass) to have
%
\begin{align*}
  v_4 &= \left(1 - \frac{\rho_5}{\rho_4}\right) v_s
  \\
  &= a_5 \left(1 - \frac{\rho_5}{\rho_4}\right)
  \sqrt{\frac{\gamma + 1}{2\gamma}\left(\frac{p_4}{p_5} - 1\right) + 1}
  \\
  &=
  \frac{a_5}{\gamma} \left(\frac{p_4}{p_5} - 1\right)
  \left[
    \left(\frac{\gamma+1}{2\gamma}\right) \left(\frac{p_4}{p_5} - 1\right) + 1
    \right]^{-\frac{1}{2}}
\end{align*}
%
Further simplify to get
%
\begin{align}
  v_4 = \frac{a_5}{\gamma} \left(\frac{p_4}{p_5} - 1\right)
  \left(
  \dfrac{\dfrac{2\gamma}{\gamma+1}}
  {\dfrac{p_4}{p_5} + \dfrac{\gamma-1}{\gamma+1}}
  \right)^{\frac{1}{2}}
  \label{e:nshock:v_4}
\end{align}
%
The value of the speed of sound (assuming calorically ideal gas) can be
obtained by
%
\begin{align*}
  a = \sqrt{\gamma R T} = \sqrt{\frac{\gamma p}{\rho}}
\end{align*}

\subsubsection{Expansion Wave}

By using the method of characteristics, the solution of the expansion wave in
zone 2 is obtained.
The speed of sound is
%
\begin{align}
  \frac{a}{a_1} = 1 - \frac{\gamma - 1}{2}\frac{v}{a_1}
  \label{e:ewave:speedofsound}
\end{align}
%
Aided by the relation $a = \sqrt{\gamma R T}$
%
\begin{align}
  \frac{T}{T_1} =
  \left[1 - \frac{\gamma - 1}{2}\left(\frac{v}{a_1}\right)\right]^2
  \label{e:ewave:temperature}
\end{align}
%
\begin{align}
  \frac{p}{p_1} =
  \left[1 - \frac{\gamma - 1}{2}\left(\frac{v}{a_1}\right)\right]^{
    \frac{2\gamma}{\gamma-1}}
  \label{e:ewave:pressure}
\end{align}
%
\begin{align}
  \frac{\rho}{\rho_1} =
  \left[1 - \frac{\gamma - 1}{2}\left(\frac{v}{a_1}\right)\right]^{
    \frac{2}{\gamma-1}}
  \label{e:ewave:density}
\end{align}

Plugging Eq.~(\ref{e:ewave:speedofsound}) into the characteristic (passing
origin)
%
\begin{align*}
  \frac{\dif x}{\dif t} = v - a
  \quad \Rightarrow \quad x = (v - a) t
\end{align*}
%
and obtain
%
\begin{align}
  v = \frac{2}{\gamma + 1}\left(a_1 + \frac{x}{t}\right)
\end{align}
%
The head of the left-running expansion wave moves at $-a_1$, the speed of sound
in zone 1.

\subsubsection{Connect the Shock Wave and Expansion Wave Solutions}

The relationship between the shock wave and expansion wave in the shock tube
can be obtained by the condition across the contact surface, $v_3 = v_4 = v_f$
and $p_3 = p_4$.
Recall Eq.~(\ref{e:nshock:v_4}) for the flow speed in zone 4
%
\begin{align}
  v_4 = \frac{a_5}{\gamma} \left(\frac{p_4}{p_5} - 1\right)
  \left(
  \dfrac{\dfrac{2\gamma}{\gamma+1}}
  {\dfrac{p_4}{p_5} + \dfrac{\gamma-1}{\gamma+1}}
  \right)^{\frac{1}{2}}
  \tag{\ref{e:nshock:v_4}}
\end{align}
%
Recall Eq.~(\ref{e:ewave:pressure})
%
\begin{align}
  \frac{p_3}{p_1} =
  \left[1 - \frac{\gamma - 1}{2}\left(\frac{v_3}{a_1}\right)\right]^{
    \frac{2\gamma}{\gamma-1}}
  \tag{\ref{e:ewave:pressure}}
\end{align}
%
The flow speed in zone 3 and 4 can be expressed by using the pressure ratio as
%
\begin{align*}
  v_3 = \frac{2 a_1}{\gamma - 1}
  \left[
  1 - \left(\frac{p_3}{p_1}\right)^{\frac{\gamma-1}{2\gamma}}
  \right]
\end{align*}
%
Because $p_3 = p_4$
%
\begin{align}
  v_3 = \frac{2 a_1}{\gamma - 1}
  \left[
  1 - \left(\frac{p_4}{p_1}\right)^{\frac{\gamma-1}{2\gamma}}
  \right]
  \label{e:tube:v_3}
\end{align}
%
Because $v_3 = v_4$, by combining Eq.~(\ref{e:nshock:v_4}) and the above
Eq.~(\ref{e:tube:v_3}), we have
%
\begin{align*}
  \frac{a_5}{\gamma} \left(\frac{p_4}{p_5} - 1\right)
  \left(
  \dfrac{\dfrac{2\gamma}{\gamma+1}}
  {\dfrac{p_4}{p_5} + \dfrac{\gamma-1}{\gamma+1}}
  \right)^{\frac{1}{2}}
  = \frac{2 a_1}{\gamma - 1}
  \left[
  1 - \left(\frac{p_4}{p_1}\right)^{\frac{\gamma-1}{2\gamma}}
  \right]
\end{align*}
%
It can be rearranged to relate $\frac{p_1}{p_5}$ to $\frac{p_4}{p_5}$ as
%
\begin{align}
  \frac{p_1}{p_5} = \frac{p_4}{p_5}
  \left\{
  1 - \frac
  {(\gamma-1)\dfrac{a_5}{a_1}\left(\dfrac{p_4}{p_5} - 1\right)}
  {\sqrt{2\gamma
  \left[
    2\gamma + (\gamma+1)\left(\dfrac{p_4}{p_5}-1\right)
    \right]}}
  \right\}^{-\frac{2\gamma}{\gamma-1}}
\end{align}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
\chapter{Two- and Three-Dimensional Unstructured Mesh}
\label{c:ustm}
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
\clearpage
\addcontentsline{toc}{chapter}{Bibliography}
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\bibliographystyle{myunsrtnat} % no sort (order in appearance)
\bibliographystyle{myplainnat} % sort by author
\bibliography{turgon_main}

\end{document}

% vim: set ff=unix fenc=utf8 et sw=2 ts=2 tw=79: