I wrote a program in Fortran which calculates Sod shock tube numerically. Now I want to compare it with exact (analytical) solution, but I don't know how to get it.

Can I find it somewhere or do I have to write a new program which calculates it? If it's the latter, where can I find the equations I need?

## Best Answer

For the Sod shock tube, there are 5 states, two of which come from initial conditions (states I & V in the linked Wiki entry). Two of the remaining 3 can be calculated with relative ease given the linked Wikipedia entry; the most difficult state (III) requires a root solver to find the value of the pressure.

What's not covered in the Wiki entry is state II, which is the rarefaction wave region. In this region, we have that: \begin{align} u_{2}(x)&=\frac{2}{\gamma+1}\cdot\left(c_1+\frac{\left(x-x_\text{mid}\right)}{t}\right)\\ \rho_{2}(x)&=\rho_1\cdot\left(1-\frac{\gamma-1}{2}\cdot\frac{u_{2}(x)}{c_1}\right)^{2/(\gamma-1)}\\ P_{2}(x)&=P_1\cdot\left(1-\frac{\gamma-1}{2}\cdot\frac{u_{2}(x)}{c_1}\right)^{2\gamma/(\gamma-1)} \end{align} where $x_\text{mid}$ is the point that divides the two states (often taken to be 0.5 with boundaries at 0 & 1), $t$ the time since the simulation started and $c_1$ the speed of sound in state I (the left state)--note that this assumes that the left state is the over-pressure region--all other terms take their normal meaning.

Building a function or two to find the states should be pretty straight-forward, once the above is known.