# Sod’s Shock Tube (Under Development)¶

A classic example to verify whether a CFD algorithm the Sod shock tube problem [Sod78]. We will introduce this problem in what follows.

In short, a shock tube problem is a Riemann problem with the Euler equations. This is a good benchmark to compare different CFD algorithm results.

The Euler equations consist of conservation of mass (Eq. (1)), of momentum (Eq. (2)), and of energy (Eq. (3)).

(1)\begin{align}\begin{aligned}\newcommand{\dpd}[]{\mathinner{ \dfrac{\partial{^{#1}}#2}{\partial{#3^{#1}}} }}\\\dpd{\rho}{t} + \dpd{{\rho}{v}}{x} = 0\end{aligned}\end{align}
(2)$\dpd{\rho{v}}{t} + \dpd{(p+\rho{v^2})}{x} = 0$
(3)$\dpd{}{t}\left(\frac{p}{\gamma-1} + \frac{\rho{v^2}}{2}\right) + \dpd{}{x}\left(\frac{\gamma}{\gamma-1}pv+\frac{1}{2}\rho{v^3}\right) = 0$

By defining

(4)\begin{align}\begin{aligned}\newcommand{\bvec}{\mathbf{#1}} \newcommand{\defeq}{\buildrel{\text{def}}\over{=}}\\\begin{split}\bvec{u} = \left(\begin{array}{c} u_1 \\ u_2 \\ u_3 \end{array}\right) \defeq \left(\begin{array}{c} \rho \\ \rho v \\ \rho\left(\frac{1}{\gamma-1}\frac{p}{\rho} + \frac{v^2}{2}\right) \end{array}\right)\end{split}\end{aligned}\end{align}
(5)$\begin{split}\bvec{f} = \left(\begin{array}{c} f_1 \\ f_2 \\ f_3 \end{array}\right) \defeq \left(\begin{array}{c} u_2 \\ (\gamma-1)u_3 - \frac{\gamma-3}{2}\frac{u_2^2}{u_1} \\ \gamma\frac{u_2u_3}{u_1} - \frac{\gamma-1}{2}\frac{u_2^3}{u_1^2} \end{array}\right)\end{split}$

we can rewrite Eqs. (1), (2), and (3) in a general form for nonlinear hyperbolic PDEs:

(6)$\dpd{\bvec{u}}{t} + \dpd{\bvec{f}(\bvec{u})}{x} = 0$

The initial condition of the Riemann problem is defined as:

(7)$\begin{split}\bvec{u} = \left(\begin{array}{c} \rho_L \\ u_L \\ p_L \end{array}\right) \text{ for } x <= 0 \text{ and } \bvec{u} = \left(\begin{array}{c} \rho_R \\ u_R \\ p_R \end{array}\right) \text{ for } x > 0\end{split}$

By using Eq. (7), Sod’s initial conditions can be set as:

(8)$\begin{split}\bvec{u} = \left(\begin{array}{c} 1 \\ 0 \\ 1 \end{array}\right) \defeq \bvec{u}_L \text{ for } x <= 0 \text{ and } \bvec{u} = \left(\begin{array}{c} 0.125 \\ 0 \\ 0.1 \end{array}\right) \defeq \bvec{u}_R \text{ for } x > 0 \text{at } t=0\end{split}$

We divide the solution of the problem in “5 zones”. From the left ($$x<0$$) to the right ($$x>0$$) of the diaphragm.

• Region I
• There is no boundary of the tube. The status is always $$\bvec{u}_L$$.
• Region II
• Rarefaction wave. The status is continuous from the region 1 to the region 3.
• Region III
• In the shock “pocket”, there is “no more shock” and the hyperbolic PDE (6) told us $$u_{\mathrm{III}}=u_{\mathrm{IV}}$$ are Riemann invariants. Together with Rankine-Hugoniot conditions, we know $$p_{\mathrm{III}}=p_{\mathrm{IV}}$$ and the density is not continuous.
• Region IV
• Because of the expansion of the shock, there is shock discontinuity. The discontinuity status could be determined by Rankine-Hugoniot conditions [Wesselling01].
• Region V
• There is no boundary of the tube, so the status is always $$\bvec{u}_R$$

To derive the analytic solution, we will begin from the region $$\mathrm{IV}$$ to get $$\bvec{u}_{\mathrm{IV}}$$, then $$\bvec{u}_{\mathrm{III}}$$ and finally $$\bvec{u}_{\mathrm{II}}$$.

## Derivation of $$\mathbf{u}_{\mathrm{IV}}$$¶

Rankine-Hugoniot conditions gives that the jump conditions must hold across a shock:

(9)$u_{shock}(\rho_{2} - \rho_{1}) = m_2 - m_1$
(10)$u_{shock}(m_2 - m_1) = \frac{{m_2}^2}{\rho_2} + p_2 - \frac{{m_1}^2}{\rho_1} - p_1$
(11)$u_{shock}(\rho_{2} E_2 - \rho_{1} E_1) = m_{2} H_{2} - m_{1} H_{1}$

If this is a stationary shock, $$u_{shock} = 0$$. (9) tells us $$m_2 = m_1$$.

Because $$u_{shock} = 0$$ and (9), from (10) we get:

$\begin{split}\frac{{m_2}^2}{\rho_2} + p_2 - \frac{{m_1}^2}{\rho_1} - p_1 = 0 \\ \text{divided by } m_1 \Rightarrow \frac{{m_2}^2}{\rho_{2}m_1} + \frac{p_2}{m_1} - \frac{{m_1}^2}{\rho_1{m_1}} - \frac{p_1}{m_1} = 0 \\ \text{please remember } m_1 = m_2 \Rightarrow \frac{m_{2}}{\rho_2} + \frac{p_2}{m_2} - \frac{m_{1}}{\rho_1} - \frac{p_1}{m_1}=0 \\ \text{use } m_1 = \rho_{1}{u_1}, m_2 = \rho_{2}{u_2} \Rightarrow u_2 + \frac{{\gamma}{p_2}}{\gamma{\rho_{2}{u_2}}} - u_1 - \frac{{\gamma}{p_1}}{\gamma{\rho_{1}{u_1}}} \\ \text{because } c_1 = \sqrt{\frac{{\gamma}{p_1}}{\rho_1}}, c_2 = \sqrt{\frac{{\gamma}{p_2}}{\rho_2}} \Rightarrow u_2 + \frac{{c_2}^2}{u_2} - u_1 - \frac{{c_1}^2}{u_1} = 0\end{split}$

Thus, we get

(12)$u_1 - u_2 = \frac{{c_2}^2}{u_2} - \frac{{c_1}^2}{u_1}$

Since $$u_{shock} = 0$$, $$m_2 = m_1$$ and (11), we get $$H_1 = H_2$$. Use $$H=h+\frac{{c}^2}{2}$$, namely $$H_1=h_1+\frac{{c_1}^2}{2}$$ and $$H_2=h_2+\frac{{c_2}^2}{2}$$, and we could rewrite $$H_1 = H_2$$ as

$\begin{split}& H_1 = h_1+\frac{{u_1}^2}{2} = h_2+\frac{{u_2}^2}{2} = H_2 \\ & \text{Use } h = c_{p}T = \frac{c^2}{\gamma - 1} \\ & \text{that is } \quad h_1 = c_{p}T_1 = \frac{{c_1}^2}{\gamma - 1} \quad h_2 = c_{p}T_2 = \frac{{c_2}^2}{\gamma - 1} \\ \Rightarrow \quad & h_1 + \frac{{u_1}^2}{2} = \frac{{c_1}^2}{\gamma - 1} + \frac{{u_1}^2}{2} \\ \quad & h_2 + \frac{{u_2}^2}{2} = \frac{{c_2}^2}{\gamma - 1} + \frac{{u_2}^2}{2} \\ \Rightarrow \quad & \frac{{c_1}^2}{\gamma - 1} + \frac{{u_1}^2}{2} = \frac{{c_2}^2}{\gamma - 1} + \frac{{u_2}^2}{2}\end{split}$

Assume $$u_1 > \text{sonic speed} c_{*} > u_2$$. Because of continuity, there must be a point with the speed $$u_{*}$$ equal to the sound speed $$c_{*}$$ which satisfies:

$\frac{{c_1}^2}{\gamma - 1} + \frac{{u_{*}}^2}{2} = \frac{{c_1}^2}{\gamma - 1} + \frac{{c_{*}}^2}{2} = \frac{(\gamma-1)+2}{2(\gamma-1)}{c^2_{*}} = \frac{(\gamma+1)}{2(\gamma-1)}{c^2_{*}}$

And

(13)$\frac{{c_1}^2}{\gamma - 1} + \frac{{u_1}^2}{2} = \frac{{c_2}^2}{\gamma - 1} + \frac{{u_2}^2}{2} = \frac{(\gamma+1)}{2(\gamma-1)}{c^2_{*}}$

Now let’s try to get $$c_{*}$$ represented by $$u_{1}$$ and $$u_{2}$$. Because of (13)

$\begin{split}& \frac{{c_1}^2}{\gamma - 1} + \frac{{u_1}^2}{2} = \frac{(\gamma+1)c_{*}}{2(\gamma-1)} \\ & \frac{{c_2}^2}{\gamma - 1} + \frac{{u_2}^2}{2} = \frac{(\gamma+1)c_{*}}{2(\gamma-1)} \\ \text{multipled by } \frac{(2\gamma-1)}{\gamma{u_1}} & \text{ and multipled by } \frac{(2\gamma-1)}{\gamma{u_2}} \text{ seperately} \\ \Rightarrow & \frac{2{c_1}^2}{\gamma{u_1}} + \frac{{u_1}(\gamma-1)}{\gamma} = \frac{(\gamma+1)c_{*}}{\gamma{u_1}} \\ & \frac{2{c_2}^2}{\gamma{u_2}} + \frac{{u_2}(\gamma-1)}{\gamma} = \frac{(\gamma+1)c_{*}}{\gamma{u_2}} \\ \Rightarrow & \frac{2{c_1}^2}{\gamma{u_1}} = \frac{(\gamma+1)c_{*}}{\gamma{u_1}} - \frac{{u_1}(\gamma-1)}{\gamma} \\ & \frac{2{c_2}^2}{\gamma{u_2}} = \frac{(\gamma+1)c_{*}}{\gamma{u_2}} - \frac{{u_2}(\gamma-1)}{\gamma} \\ \Rightarrow & \frac{{c_1}^2}{\gamma{u_1}} = \frac{(\gamma+1)c_{*}}{2\gamma{u_1}} - \frac{{u_1}(\gamma-1)}{2\gamma} = [\frac{(\gamma+1)c_{*}}{\gamma-1}+{u_1}^2] (\frac{(\gamma-1)}{2\gamma{u_1})}) \\ & \frac{{c_2}^2}{\gamma{u_2}} = \frac{(\gamma+1)c_{*}}{2\gamma{u_2}} - \frac{{u_2}(\gamma-1)}{2\gamma} = [\frac{(\gamma+1)c_{*}}{\gamma-1}+{u_2}^2] (\frac{(\gamma-1)}{2\gamma{u_2})}) \\ \Rightarrow & \frac{{c_1}^2}{\gamma{u_1}} - \frac{{c_2}^2}{\gamma{u_2}} = \frac{(\gamma+1)c_{*}}{2\gamma{u_1}} - \frac{{u_1}(\gamma-1)}{2\gamma} - \frac{(\gamma+1)c_{*}}{2\gamma{u_2}} + \frac{{u_2}(\gamma-1)}{2\gamma}\end{split}$

$\begin{split}u_2 - u_1 & = \frac{(\gamma+1)c_{*}}{2\gamma{u_1}} - \frac{{u_1}(\gamma-1)}{2\gamma} - \frac{(\gamma+1)c_{*}}{2\gamma{u_2}} + \frac{{u_2}(\gamma-1)}{2\gamma} \\ \Rightarrow & c_{*}(\frac{(\gamma+1)}{2\gamma{u_1}} - \frac{(\gamma+1)}{2\gamma{u_2}}) = u_2\frac{\gamma+1}{2\gamma} + u_1\frac{\gamma+1}{2\gamma} \\ \Rightarrow & c_{*}(\frac{1}{u_1}-\frac{1}{u_2}) = u_2 - u_1 \\ \Rightarrow & c_{*} = {u_1}{u_2}\end{split}$
(14)$c_{*} = {u_1}{u_2}$
Now defining the Mach number $$M_1 = \frac{u_1}{c_1}$$ and $$M^{*} = \frac{u}{c_*}$$, we get from (13):
$\begin{split}& \frac{{c_1}^2}{(\gamma - 1)} + \frac{{u_1}^2}{2} = \frac{(\gamma+1)}{2(\gamma-1)}{c^2_{*}} \\ \Rightarrow & \frac{{c_1}^2}{(u_1)^2}\frac{(u_1)^2}{\gamma - 1} + \frac{{u_1}^2}{2} = \frac{(\gamma+1)}{2(\gamma-1)}{c^2_{*}} \\ \text{dividied by } c^2_{*} \Rightarrow & \frac{{c_1}^2}{(u_1)^2}\frac{(u_1)^2}{c^2_{*}(\gamma - 1)} + \frac{{u_1}^2}{(c^2_{*})2} = \frac{(\gamma+1)}{2(\gamma-1)} \\ & \frac{M^*_1}{M_1(\gamma-1)} + \frac{M^{*2}_{1}}{2} = \frac{{u_1}^2}{(c^2_{*})2}\end{split}$