Sod's Test Problems¶

Sod, 1978, J.Comp.Phys. Vol 27: 1-31

C. Laney, 1998, Computational Gas Dynamics, Cambridge University Press

You can read more about the Sod Shock Tube on Wikipedia.

The shock tube is a common test for computational solvers and it has an analytical solution for comparison of results.

The shock tube problem set up is as follows. Imagine a tube divided by a massless, rigid membrane separating areas of different pressure and density. If the membrane is instantaneously removed, what is the reaction within the tube?

Test 1:¶

Find the pressure, velocity, speed of sound, density, entropy and Mach number at t = 0.01 s, where¶

$$\underline{w}_L = \left[ \begin{array}{c} \rho_L \\ u_L \\ p_L \\ \end{array}\right] = \left[ \begin{array}{c} 1\ kg/m^3 \\ 0\ m/s \\ 100\ kN/m^2 \\ \end{array}\right]$$$$\underline{w}_R = \left[ \begin{array}{c} \rho_R \\ u_R \\ p_R \\ \end{array}\right] = \left[ \begin{array}{c} 0.125\ kg/m^3 \\ 0\ m/s \\ 10\ kN/m^2 \\ \end{array}\right]$$

Test 2:¶

Find the pressure, velocity, speed of sound, density, entropy and Mach number at t = 0.01 s, where¶

$$\underline{w}_L = \left[ \begin{array}{c} \rho_L \\ u_L \\ p_L \\ \end{array}\right] = \left[ \begin{array}{c} 1\ kg/m^3 \\ 0\ m/s \\ 100\ kN/m^2 \\ \end{array}\right]$$$$\underline{w}_R = \left[ \begin{array}{c} \rho_R \\ u_R \\ p_R \\ \end{array}\right] = \left[ \begin{array}{c} 0.010\ kg/m^3 \\ 0\ m/s \\ 1\ kN/m^2 \\ \end{array}\right]$$

The Conserved Euler Equations¶

The conservative formulation of the Euler equations is required when solutions contain shocks. As you might expect, the Sod's Shock Tube solutions do contain shocks. The conserved variables $\underline{\mathbf{u}}$ are

$$\underline{\mathbf{u}} = \left[ \begin{array}{c} \rho \\ \rho u \\ \rho e_T \\ \end{array} \right]$$

where $e_T = e + \frac{u^2}{2}$ is the specific total energy.

The vector form of the Euler equations is:

$$\frac{\partial }{\partial t} \underline{\mathbf{u}} + \frac{\partial }{\partial t} \underline{\mathbf{f}} = 0$$

where $\underline{\mathbf{f}}$ is the flux vector:

$$\underline{\mathbf{f}} = \left[ \begin{array}{c} \rho u \\ \rho u^2 + p \\ (\rho e_T + p) u \\ \end{array} \right] = \left[ \begin{array}{c} \rho u \\ \rho u^2 + p \\ \rho h_T u \\ \end{array} \right] $$

and where $h_T$ is the total enthalpy:

$$h_T = e_T + \frac{p}{\rho} = h + \frac{u^2}{2}$$

Using the Jacobian matrix, the Euler equations are written in quasi-linear form:

$$\frac{\partial}{\partial t} \underline{\mathbf{u}} + A \frac{\partial}{\partial x}\underline{\mathbf{u}} = 0$$

where

$$A(\underline{\mathbf{u}}) = \frac{\partial \underline{\mathbf{F}}}{\partial \underline{\mathbf{u}}} = \left[ \begin{array}{ccc} \frac{\partial f_1}{\partial u_1} & \frac{\partial f_1}{\partial u_2} & \frac{\partial f_1}{\partial u_3}\\ \frac{\partial f_2}{\partial u_1} & \frac{\partial f_2}{\partial u_2} & \frac{\partial f_2}{\partial u_3}\\ \frac{\partial f_3}{\partial u_1} & \frac{\partial f_3}{\partial u_2} & \frac{\partial f_3}{\partial u_3}\\ \end{array} \right]$$

Euler Jacobian Matrix¶

The flux is given as:

$$\underline{\mathbf{F}} = \left[ \begin{array}{c} \rho u \\ \rho u^2 + p \\ (\rho e_T + p) u \\ \end{array} \right] = \left[ \begin{array}{c} f_1 \\ f_2 \\ f_3 \\ \end{array} \right]$$

To calculate the Jacobian, $\underline{\mathbf{A}}$, we need $\underline{\mathbf{F}}$ in terms of $u_1$, $u_2$ and $u_3$.

Recalling that $\underline{\mathbf{u}} = \left[ \begin{array}{c} \rho \\ \rho u \\ \rho e_T \\ \end{array} \right] = \left[ \begin{array}{c} u_1 \\ u_2 \\ u_3 \\ \end{array} \right]$ we note that $f_1 = \rho u = u_2$

For $f_2$ and $f_3$, however, we must first represent the pressure $p$ in terms of conserved variables. We can relate the pressure to the conserved variables using an equation of state.

Relating pressure to conserved variables¶

For an ideal gas, the equation of state is

$$e = e(\rho, p) = \frac{p}{(\gamma -1) \rho}, \gamma=1.4$$$$\therefore p = (\gamma -1)\rho e $$$$e_T = e+\frac{1}{2} u^2$$$$\therefore e = e_T - \frac{1}{2}u^2$$$$\therefore p = (\gamma -1)\left(\rho e_T - \frac{\rho u^2}{2}\right) = (\gamma -1)\left(u_3 - \frac{1}{2} \frac{u^2_2}{u_1} \right)$$

Substituting for the pressure in terms of $\underline{\mathbf{u}}$ yields

$$\underline{\mathbf{F}} = \left[ \begin{array}{c} f_1 \\ f_2 \\ f_3 \\ \end{array} \right] = \left[ \begin{array}{c} u_2\\ \frac{u^2_2}{u_1} + (\gamma -1)\left(u_3 - \frac{1}{2} \frac{u^2_2}{u_1} \right) \\ \left(u_3 + (\gamma -1)\left(u_3 - \frac{1}{2} \frac{u^2_2}{u_1}\right) \right) \frac{u_2}{u_1}\\ \end{array} \right]$$

Which can be simplied further to

$$\underline{\mathbf{F}} = \left[ \begin{array}{c} f_1 \\ f_2 \\ f_3 \\ \end{array} \right] = \left[ \begin{array}{c} u_2\\ (3-\gamma) \frac{u^2_2}{2u_1} + (\gamma - 1)u_3 \\ \gamma \frac{u_2}{u_1}u_3 - \frac{\gamma -1}{2} \frac{u^3_2}{u^2_1} \\ \end{array} \right]$$

Calculating the Jacobian¶

Recalling that the Jacobian $\underline{\mathbf{A}}$ is of the form:

$$A(\underline{\mathbf{u}}) = \frac{\partial \underline{\mathbf{F}}}{\partial \underline{\mathbf{u}}} = \left[ \begin{array}{ccc} \frac{\partial f_1}{\partial u_1} & \frac{\partial f_1}{\partial u_2} & \frac{\partial f_1}{\partial u_3}\\ \frac{\partial f_2}{\partial u_1} & \frac{\partial f_2}{\partial u_2} & \frac{\partial f_2}{\partial u_3}\\ \frac{\partial f_3}{\partial u_1} & \frac{\partial f_3}{\partial u_2} & \frac{\partial f_3}{\partial u_3}\\ \end{array} \right]$$

Calculate all partial derivatives to determine that the Jacobian matrix is

$$A(\underline{\mathbf{u}}) = \left[ \begin{array}{ccc} 0 & 1 & 0 \\ \frac{1}{2} (\gamma - 3) \left(\frac{u_2}{u_1}\right)^2 & (3- \gamma)\frac{u_2}{u_1} & \gamma - 1 \\ -\gamma \frac{u_2 u_3}{u_1^2} + (\gamma -1)\left(\frac{u_2}{u_1}\right)^3 & \gamma \frac{u_3}{u_1} - \frac{3}{2}(\gamma -1)\left(\frac{u_2}{u_1}\right)^2 & \gamma \frac{u_2}{u_1} \\ \end{array} \right]$$

Expressing the Jacobian in terms of the velocity, $u$ and the speed of sound (for an ideal gas) $$a = \sqrt{\gamma \frac{p}{\rho}}$$ yields

$$A(\underline{\mathbf{u}}) = \left[ \begin{array}{ccc} 0 & 1 & 0 \\ \frac{1}{2} (\gamma - 3) u^2 & (3 - \gamma)u & \gamma - 1 \\ \frac{1}{2}(\gamma -2)u^3 - \frac{a^2 u}{\gamma -1} & \frac{3-2 \gamma}{2}u^2+\frac{a^2}{\gamma -1} & \gamma u\\ \end{array} \right]$$

Laying the groundwork¶

There's a lot to keep track of in the Sod tests, so it is a good idea to build several functions which can take care several common operations.

Initial Conditions¶

The initial conditions for both Sod tests are below.

Test1¶

$$\underline{w}_L = \left[ \begin{array}{c} \rho_L \\ u_L \\ p_L \\ \end{array}\right] = \left[ \begin{array}{c} 1\ kg/m^3 \\ 0\ m/s \\ 100\ kN/m^2 \\ \end{array}\right]$$$$\underline{w}_R = \left[ \begin{array}{c} \rho_R \\ u_R \\ p_R \\ \end{array}\right] = \left[ \begin{array}{c} 0.125\ kg/m^3 \\ 0\ m/s \\ 10\ kN/m^2 \\ \end{array}\right]$$

Test 2¶

$$\underline{w}_L = \left[ \begin{array}{c} \rho_L \\ u_L \\ p_L \\ \end{array}\right] = \left[ \begin{array}{c} 1\ kg/m^3 \\ 0\ m/s \\ 100\ kN/m^2 \\ \end{array}\right]$$$$\underline{w}_R = \left[ \begin{array}{c} \rho_R \\ u_R \\ p_R \\ \end{array}\right] = \left[ \begin{array}{c} 0.010\ kg/m^3 \\ 0\ m/s \\ 1\ kN/m^2 \\ \end{array}\right]$$

Instead of writing all of this out for every test, we can create functions to "reset" all of our initial conditions.

We're also using a standard discretization throughout, so we can roll that into our reset functions.

Discretization¶

Test 1¶

• N = 50 points in [-10m, 10m]
• $\Delta x$ = 20m / 50 = 0.4 m
• initial CFL = 0.4

Test 2¶

• N = 50 points in [-10m, 15m]
• $\Delta x$ = 25m / 50 = 0.5 m
• initial CFL = 0.3

Initial maximum wave speed = 374.17 m/s

Conserved Variables Vector¶

Recall that

$$\underline{\mathbf{u}} = \left[ \begin{array}{c} \rho \\ \rho u \\ \rho e_T \\ \end{array} \right]$$

where $e_T = e + \frac{u^2}{2}$ is the specific total energy and

$$e = e(\rho, p) = \frac{p}{(\gamma -1) \rho}, \gamma=1.4$$

The Flux Vector¶

The Jacobian¶

Discretizing the Euler Equations¶

Lax Friedrichs scheme (1st order)¶

$$\underline{u}^{n+1}_i = \frac{1}{2} \left(\underline{u}^n_{i+1} + \underline{u}^n_{i-1} \right) - \frac{\Delta t}{2 \Delta x}\left[ \underline{f}\left(\underline{u}^n_{i+1}\right) - \underline{f}\left(\underline{u}^n_{i-1} \right) \right]$$

Making sense of the results¶

We've solved for the vector $\underline{\mathbf{u}}$ using the conditions specified, but the goal is solve for pressure, velocity, sound speed, density, entropy and Mach number.

$\underline{\mathbf{u}}$ needs to be decomposed into those elements in order to compare results with the analytical solution.

Richtmyer method¶

2-step method

$$\underline{\mathbf{u}}^{n+\frac{1}{2}}_{i+\frac{1}{2}} = \frac{1}{2} \left( \underline{\mathbf{u}}^n_{i+1} + \underline{\mathbf{u}}^n_i \right) - \frac{\Delta t}{2 \Delta x} \left( \underline{\mathbf{f}}^n_{i+1} - \underline{\mathbf{f}}^n_i\right)$$$$\underline{\mathbf{u}}^{n+1}_i = \underline{\mathbf{u}}^n_i - \frac{\Delta t}{\Delta x} \left(\underline{\mathbf{f}}^{n+\frac{1}{2}}_{i+\frac{1}{2}} - \underline{\mathbf{f}}^{n+\frac{1}{2}}_{i-\frac{1}{2}} \right)$$

Richtmyer, Test 1¶

Using the values for test 1, let's see how the Richtmyer scheme handles the shock tube.

Damping¶

Some schemes are prone to overshoots, which produce the sharp spikes in the plots above. One way to smooth out these overshoots is to provide a damping term.

For the Richtmyer scheme, the damping term is added when the $\underline{\mathbf{u}}^{n+1}_i$ term is calculated.

The damping term is simply added on to the existing equation, with the level of damping controlled by a damping coefficient $\epsilon$.

The full damping term is defined as:

$$\epsilon (\underline{\mathbf{u}}^{n}_{i+1} - 2\underline{\mathbf{u}}^{n}_{i} + \underline{\mathbf{u}}^{n}_{i-1})$$

You can see that this is implemented in the Richtmyer function above with the line

damp * (U[:,2:] - 2*U[:,1:-1] + U[:,:-2])

where damp is another parameter passed to the function. damp is given a default value of 0, so to add damping, simply pass a non-zero value of damp to the same function.

Try experimenting with different values of $\epsilon$ and see what happens. Below, a value of $\epsilon = 0.1$ is used.

Look at the difference! Just a small damping value and all of the spikes smooth out nicely.

Richtmyer, Test 2¶

As with the first set of conditions, there are some spikes formed by overshoots. Applying damping with $\epsilon = 0.125$, yields a similar smoothing effect.

MacCormack Method¶

This method is also called the 'predictor-corrector' method, where we predict the next value of $\underline{\mathbf{u}}$, $\underline{\mathbf{u}}^*_i$, and then update it.

$$\underline{\mathbf{u}}^*_i = \underline{\mathbf{u}}^n_i - \frac{\Delta t}{\Delta x} \left( \underline{\mathbf{f}}^n_{i+1} - \underline{\mathbf{f}}^n_i \right)$$$$\underline{\mathbf{u}}^{n+1}_i = \frac{1}{2}\left(\underline{\mathbf{u}}^n_i + \underline{\mathbf{u}}^*_i \right) - \frac{\Delta t}{2 \Delta x} \left(\underline{\mathbf{f}}^*_i - \underline{\mathbf{f}}^*_{i-1} \right) $$

Damping¶

Damping can also be used to control overshoots in the MacCormack method. Unlike the Richtmyer damping, however, the damping used with the MacCormack method is applied to the 'predictor' step of the scheme. Again, the default value of damp is 0, so the undamped behavior can also be observed.

MacCormack, Test 1¶

Catastrophic negative overshoots¶

What happened? If the overshoots are very large, the value of $\underline{\mathbf{u}}$ can become negative which will introduce some large errors into the calculation. It will also create some complex numbers when we try to calculate anything involving square roots.

Remember that the Sod tests are very restrictive, in the sense that the initial conditions we are using are very 'coarse.' There are a few ways to mitigate negative overshoots available to us.

1. Check and 'correct' negative velocity, density and pressure values at each time iteration.

2. Introduce damping and see if it can overcome the negative overshoots.

3. Cheat. Or rather, try out a third set of test values with a much finer mesh and see if the problem persists.

4. Cheat a different way. Instead of using a fixed timestep, introduce an adaptive timestep. Remember that we define out timestep nt in terms of the initial sound speed. More on this later.

For now, let's try introducing some damping, with $\epsilon = .125$, and see what happens.

That worked much better. How does Test 2 fare with damping added?

MacCormack, Test 2¶

MacCormack, 'Test 3'¶

What about the first 'cheat'? If we change the problem and give it a much finer mesh to work with, will the overshoots still cause trouble? Below is a third set of initial conditions, with the same physical parameters as test 2, but with nx = 400.

Even with a much finer mesh, the overshoots have apparently caused some trouble. But maybe we can get away with a much smaller damping coefficient?

---

The cell below executes the style of this notebook.