Shock tube calculator

This module calculates analytical solutions for a Sod shock tube problem.

Standard test case

A shock tube consists of a long tube filled with the same gas in two different physical states. The tube is divided into two parts, separated by a diaphragm. The initial state is defined by the values for density, pressure and velocity.

When the membrane bursts, discontinuity between the two initial states breaks into leftward and rightward moving waves, separated by a contact surface.

• Region 5 is the low pressure gas which is not disturbed by the shock wave.
• Regions 3 and 4 (divided by the contact surface) contain the gas immediately behind the shock traveling at a constant speed.
• The contact surface across which the density and the temperature are discontinuous lies within this zone.
• The zone between the head and the tail of the expansion fan is noted as Region 2. In this zone, the flow properties gradually change since the expansion process is isentropic.
• Region 1 denotes the undisturbed high pressure gas.

About the code

Sod solver returns after time of evolution the following variables: 1. Positions of head and foot of rarefation wave, contact discontinuity and shock 2. Constant pressure, density and velocity for all regions except rarefaction region 3. Pressure, density and velocity sampled across the domain of interest

The usage should be straightforward:

---

```python from shocktubecalc import sod

leftstate and rightstate set p, rho and u

geometry sets left boundary on 0., right boundary on 1

and initial position of the shock xi on 0.5

t is the time evolution for which positions and states in tube should be calculated

gamma denotes specific heat

npts is number of points to be calculated for rarefaction wave

Note that gamma and npts are default parameters (1.4 and 500) in solve function.

positions, regions, values = sod.solve(leftstate=(1, 1, 0), rightstate=(0.1, 0.125, 0.), geometry=(0., 1., 0.5), t=0.2, gamma=1.4, npts=500) ```

---

Return value of a solve function

The positions is a dictionary of region names and their positions in the pipe after t=0.2 seconds.

Printing positions for each and every region: python print('Positions:') for desc, vals in positions.items(): print('{0:10} : {1}'.format(desc, vals))

The regions is a dictionary of constant pressure, density and velocity for all regions except for the rarefaction wave. python print('States:') for desc, vals in regions.items(): print('{0:10} : {1}'.format(desc, vals))

The values is a dictionary of pressure , density , velocity and x arrays. If one would want to plot said values as a function of x : ```python import matplotlib.pyplot as pl plt.figure(1) plt.plot(values['x'], values['p'], linewidth=1.5, color='b') plt.ylabel('pressure') plt.xlabel('x') plt.axis([0, 1, 0, 1.1])

plt.figure(2) plt.plot(values['x'], values['rho'], linewidth=1.5, color='r') plt.ylabel('density') plt.xlabel('x') plt.axis([0, 1, 0, 1.1])

plt.figure(3) plt.plot(values['x'], values['u'], linewidth=1.5, color='g') plt.ylabel('velocity') plt.xlabel('x')

plt.show() ```

Calculating other values

One can easily calculate energy, temperature, speed of sound and a Mach number:

• Energy e = p/(gamma-1)
• Temperature T = p/rho
• Speed of sound c = (gamma p /rho)^0.5
• Mach number M = u/c

python e = values['p']/(gamma -1) T = values['p']/values['rho'] c = np.sqrt(gamma *values['p']/values['rho']) M = values['p']/c