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}[3][]{\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}[1]{\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}\]

please recall (12), thus

\[\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}\]

The relation

(14)\[c_{*} = {u_1}{u_2}\]

is called Prandel-Meyer relation. It means the flow at one side of a shock must be supersonic, and the other side must be subsonic.

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}\]