Sod shock tube

From Wikipedia, the free encyclopedia
Jump to: navigation, search
Density Plot after time evolution of t=0.2[-]

The Sod shock tube problem, named after Gary A. Sod, is a common test for the accuracy of computational fluid codes, like Riemann solvers, and was heavily investigated by Sod in 1978.

The test consists of a one-dimensional Riemann problem with the following parameters, for left and right states of an ideal gas.

\left( \begin{array}{c}\rho_L\\P_L\\v_L\end{array}\right)
\left( \begin{array}{c}1.0\\1.0\\0.0\end{array} \right)
\left( \begin{array}{c}\rho_R\\P_R\\v_R\end{array}\right)
\left( \begin{array}{c}0.125\\0.1\\0.0\end{array}\right)


  • \rho is the density
  • P is the pressure
  • v is the velocity

The time evolution of this problem can be described by solving the Euler equations. Which leads to three characteristics, describing the propagation speed of the various regions of the system. Namely the rarefaction wave, the contact discontinuity and the shock discontinuity. If this is solved numerically, one can test against the analytical solution, and get information how well a code captures and resolves shocks and contact discontinuities and reproduce the correct density profile of the rarefaction wave.

Analytic derivation[edit]

The different states of the solution are separated by the time evolution of the three characteristics of the system, which is due to the finite speed of information propagation. Two of them are equal to the speed of sound of the both states

cs_1 = \sqrt{\gamma \frac{P_L}{\rho_L}}
cs_5 = \sqrt{\gamma \frac{P_R}{\rho_R}}

The first one is the position of the beginning of the rarefaction wave while the other is the velocity of the propagation of the shock.


\Gamma = \frac{\gamma - 1}{\gamma + 1}, \beta = \frac{\gamma - 1}{2 \gamma}

The states after the shock are connected by the Rankine Hugoniot shock jump conditions.

\rho_4 = \rho_5 \frac{P_4 + \Gamma P_5}{P_5 + \Gamma P_4}

But to calculate the density in Region 4 we need to know the pressure in that region. This is related by the contact discontinuity with the pressure in region 3 by

P_4 = P_3

Unfortunately the pressure in region 3 can only be calculated iteratively, the right solution is found when u_2 equals u_4

u_4 = \left(P_3' - P_5\right)\sqrt{\frac{1-\Gamma}{\rho_R(P_3'+\Gamma P_5)}}
u_2 =\left(P_1^\beta-P_3'^\beta\right) \sqrt{\frac{(1-\Gamma^2)P_1^{1/\gamma}}{\Gamma^2 \rho_L}}
u_2 - u_4 = 0

This function can be evaluated to an arbitrary precision thus giving the pressure in the region 3

P_3 = \operatorname{calculate}(P_3,s,s,,)

finally we can calculate

u_3 = u_5 + \frac{(P_3 - P_5)}{\sqrt{\frac{\rho_5}{2}((\gamma+1)P_3 +(\gamma-1)P_5)}}
u_4 = u_3

and \rho_3 follows from the adiabatic gas law

\rho_3 = \rho_1 \left(\frac{P_3}{P_1}\right)^{1/\gamma}


See also[edit]