ShockTube module

The ShockTube module provides functions for analyzing shock tube flows and gas-gas interfaces using the Rankine-Hugoniot relations and exact Riemann solvers.

Overview

Shock tubes are fundamental devices in experimental fluid dynamics, used to generate controlled shock waves and study high-speed gas dynamics. This module provides tools to:

Calculate shock jump conditions across moving shock waves
Determine required driver section pressures for desired shock strengths
Solve complete shock tube problems including reflected shocks
Analyze gas-gas interface dynamics using exact Riemann solvers

Basic shock tube functions

Shock jump conditions

Calculate the post-shock state for a gas traversed by a shock wave:

using PyThermo
using PyThermo.ShockTube

# Define the initial state
driven = Species("N2", P=100u"kPa", T=300u"K")

# Calculate shock jump for Mach 2.0 shock
shocked, velocity = shockjump!(driven, 2.0)

pressure(shocked)    # post-shock pressure
temperature(shocked) # post-shock temperature
velocity             # flow velocity

The shockjump! variant modifies the input state in-place:

driven = Species("N2")
shocked, velocity = shockjump!(driven, 2.0)
# driven is now modified to the shocked state
@assert shocked === driven

Driver pressure calculation

Calculate the required driver section pressure to achieve a desired shock strength:

driver = Species("He")
driven = Mixture(["He" => 0.95, "acetone" => 0.05], T=18u"°C", P=85u"kPa")
Ms = 2.2

# Calculate required driver pressure
driver_calc = driverpressure!(driver, driven, Ms)
pressure(driver_calc) # required driver pressure (~3.73 MPa)

The driverpressure! variant modifies the driver state in-place:

driver = Species("He")
driver_calc = driverpressure!(driver, driven, Ms)
@assert driver_calc === driver

Complete shock tube calculation

Perform a complete shock tube analysis including reflected shock:

driver = Species("He")
driven = Mixture(["He" => 0.95, "acetone" => 0.05], T=18u"°C", P=85u"kPa")
Ms = 2.2

result = shockcalc!(driver, driven, Ms)

# Access all states
result.driven    # initial driven state
result.shocked   # shocked state
result.reflected # reflected shock state
result.driver    # driver state

# Access flow properties
result.Ms # incident shock Mach
result.Mr # reflected shock Mach
result.u2 # post-shock velocity

# Density ratio (useful for test time estimation)
density(result.shocked) / density(result.driven) # ~2.64

Riemann interface solver

The exact Riemann solver analyzes the interaction between two adjacent gases with differing properties.

Basic usage: Direct specification

Solve for the interface state between two gases at rest:

using PyThermo.ShockTube: riemann_interface

# Define left and right states
left = Species("N2", P=500u"kPa", T=600u"K")
right = Species("SF6", P=100u"kPa", T=300u"K")

# Solve the Riemann problem
sol = riemann_interface(left, right)

# Interface properties
sol.p_star    # interface pressure
sol.u_star    # interface velocity
sol.S_contact # contact discontinuity speed

# Post-wave states
sol.rho_star_L # left star density
sol.rho_star_R # right star density
sol.T_star_L   # left star temperature
sol.T_star_R   # right star temperature

# Wave speeds
sol.S_L # left wave speed
sol.S_R # right wave speed

# Gamma values
sol.gamma_L # left gas isentropic exponent
sol.gamma_R # right gas isentropic exponent

With initial velocities

Specify non-zero initial velocities for each gas:

left = Species("N2", P=200u"kPa", T=400u"K")
right = Species("Ar", P=100u"kPa", T=300u"K")

# Left gas moving right at 100 m/s, right gas moving left at 50 m/s
sol = riemann_interface(left, right, 100.0u"m/s", -50.0u"m/s")

Shock tube convenience method

Automatically apply shock jump and solve interface problem:

# Define initial states
driven = Species("N2", P=100u"kPa", T=300u"K")
test_gas = Species("SF6", P=100u"kPa", T=300u"K")

# Shock Mach number
Ms = 2.0

# Automatically calculates shock jump in driven gas,
# then solves interface with test gas
sol = riemann_interface(driven, test_gas, Ms)

sol.p_star # interface pressure
sol.u_star # interface velocity

This is equivalent to:

# Manual approach
shocked, u_shocked = shockjump!(driven, Ms)
sol = riemann_interface(shocked, test_gas, u_shocked, 0.0)

Understanding the results

RiemannSolution structure

The RiemannSolution struct contains all information about the solved interface:

Interface Properties (same on both sides):

p_star: Interface pressure
u_star: Interface velocity
S_contact: Contact discontinuity speed (equals u_star)

Left Star Region (between left wave and contact):

rho_star_L: Density
T_star_L: Temperature
S_L: Left wave speed (shock or rarefaction tail)

Right Star Region (between contact and right wave):

rho_star_R: Density
T_star_R: Temperature
S_R: Right wave speed (shock or rarefaction tail)

Original State Information:

gamma_L: Specific heat ratio of left gas
gamma_R: Specific heat ratio of right gas

Wave structure

The Riemann solution consists of three waves:

Left Wave (speed S_L): Propagates into left gas
- If p_star > p_left: Shock wave (compression)
- If p_star < p_left: Rarefaction wave (expansion)
Contact Discontinuity (speed S_contact = u_star): Material interface
- Pressure and velocity continuous across contact
- Density and temperature discontinuous (different gases)
Right Wave (speed S_R): Propagates into right gas
- If p_star > p_right: Shock wave (compression)
- If p_star < p_right: Rarefaction wave (expansion)

Physical interpretation

left = Species("N2", P=400u"kPa", T=500u"K")
right = Species("Ar", P=100u"kPa", T=300u"K")

sol = riemann_interface(left, right)

# High pressure on left causes expansion wave into left gas
# and compression (shock) into right gas
sol.p_star < pressure(left)  # left wave is rarefaction
sol.p_star > pressure(right) # right wave is shock

# Interface moves toward lower pressure region
sol.u_star > 0.0u"m/s" # interface moving right if true, left if false

# Wave speeds bracket the contact
@assert sol.S_L < sol.S_contact < sol.S_R

Display and printing

The RiemannSolution struct has a custom display format:

sol = riemann_interface(left, right)
println(sol)

Output:

RiemannSolution:
  Interface: p* = 229 kPa, u* = 175 m/s
  Left  state: ρ = 1.81 kg/m³, T = 427 K, S_L = -454 m/s
  Right state: ρ = 2.59 kg/m³, T = 425 K, S_R = 460 m/s

Examples

Example 1: Helium-driven shock tube

Complete analysis of a helium-driven shock tube with acetone seeding:

using PyThermo
using PyThermo.ShockTube

# Define gases
driver = Species("He")
driven = Mixture(["He" => 0.95, "acetone" => 0.05], T=18u"°C", P=85u"kPa")

# Target shock Mach number
Ms = 2.2

# Complete calculation
result = shockcalc!(driver, driven, Ms)

# Shock tube results
result.Ms                                        # incident shock Mach
result.u2                                        # post-shock velocity
pressure(result.shocked)                         # post-shock pressure
temperature(result.shocked)                      # post-shock temperature
density(result.shocked) / density(result.driven) # density ratio

# Reflected shock
result.Mr                     # reflected shock Mach number
pressure(result.reflected)    # reflected shock pressure
temperature(result.reflected) # reflected shock temperature

# Driver requirements
pressure(result.driver) # required driver pressure

Example 2: Gas-gas interface analysis

Analyze the interface between shocked nitrogen and sulfur hexafluoride:

using PyThermo
using PyThermo.ShockTube: riemann_interface

# Initial states
driven = Species("N2", P=100u"kPa", T=300u"K")
test_gas = Species("SF6", P=100u"kPa", T=300u"K")

# Shock strength
Ms = 1.8

# Solve interface problem
sol = riemann_interface(driven, test_gas, Ms)

# Interface analysis
sol.p_star # interface pressure
sol.u_star # interface velocity

# Left (shocked N2) properties
sol.rho_star_L # density
sol.T_star_L   # temperature
sol.S_L        # wave speed

# Right (SF6) properties
sol.rho_star_R # density
sol.T_star_R   # temperature
sol.S_R        # wave speed

sol.S_contact # contact discontinuity speed

Example 3: Different gas combinations

Compare interface behavior for different gas pairs:

using PyThermo.ShockTube: riemann_interface

# Common left state
left = Species("He", P=300u"kPa", T=500u"K")

# Different right gases
gases = ["Ar", "N2", "SF6", "CO2"]

# Gas pair comparisons (left = He at 300 kPa, 500 K)
for gas in gases
    right = Species(gas, P=100u"kPa", T=300u"K")
    sol = riemann_interface(left, right)
    
    # He / gas results:
    sol.p_star                      # interface pressure
    sol.u_star                      # interface velocity
    sol.rho_star_L                  # left star density
    sol.rho_star_R                  # right star density
    sol.rho_star_L / sol.rho_star_R # density jump
end

API reference

PyThermo.ShockTube.shockjump! — Function

shockjump!(gas::Chemical, Mach::Real) -> (Chemical, Unitful.Velocity)
shockjump(gas::Chemical, Mach::Real) -> (Chemical, Unitful.Velocity)

Calculate shock jump conditions across a normal shock wave using the Rankine-Hugoniot relations.

The mutating version shockjump! modifies the input gas object in place, while shockjump creates a copy and leaves the original unchanged.

Parameters

gas: A Species or Mixture object representing the pre-shock gas state
Mach: Shock Mach number (ratio of shock velocity to sound speed in pre-shock gas)

Returns

A tuple containing:

Modified gas object with post-shock temperature and pressure
Post-shock gas velocity relative to the laboratory frame (with units m/s)

Physics

The function calculates:

Pressure ratio (PR) across the shock
Temperature ratio (TR) across the shock
Post-shock gas velocity using the Rankine-Hugoniot jump conditions

Examples

julia> driven = Mixture(["He" => 0.95, "acetone" => 0.05], T = 18u"°C", P = 85u"kPa")
Mixture(95% He, 5% acetone, 291.1 K, 8.500e+04 Pa)

julia> shocked, velocity = shockjump(driven, 2.2);

julia> round(Int, ustrip(u"K", temperature(shocked)))
624

julia> round(Int, ustrip(u"m/s", velocity))
1024

See Also

shockcalc!: Comprehensive shock tube calculation
driverpressure!: Calculate required driver pressure

source

PyThermo.ShockTube.driverpressure! — Function

driverpressure!(driver::Chemical, driven::Chemical, Ms::Real) -> Chemical
driverpressure(driver::Chemical, driven::Chemical, Ms::Real) -> Chemical

Calculate the required driver gas pressure to achieve a specified shock Mach number in a shock tube.

The mutating version driverpressure! modifies the input driver object in place, while driverpressure creates a copy and leaves the original unchanged.

Parameters

driver: A Species or Mixture object representing the driver gas (high-pressure section)
driven: A Species or Mixture object representing the driven gas (low-pressure test section)
Ms: Desired shock Mach number in the driven section

Returns

The modified driver object with updated pressure that will produce the specified shock Mach number.

Physics

Uses the shock tube theory relating driver and driven gas properties through the contact surface and expansion wave. The calculation accounts for:

Isentropic exponents (γ) of both gases
Sound speeds in both sections
The relationship between shock strength and pressure ratio across the contact surface

Examples

julia> driver = Species("He");

julia> driven = Mixture(["He" => 0.95, "acetone" => 0.05], T = 18u"°C", P = 85u"kPa");

julia> driver_calc = driverpressure(driver, driven, 2.2);

julia> round(typeof(1.0u"MPa"), pressure(driver_calc), digits=2)
3.73 MPa

See Also

shockjump!: Calculate shock jump conditions
shockcalc!: Comprehensive shock tube calculation

source

PyThermo.ShockTube.shockcalc! — Function

shockcalc!(driver::Chemical, driven::Chemical, Ms::Real) -> ShockCalc
shockcalc(driver::Chemical, driven::Chemical, Ms::Real) -> ShockCalc

Perform a comprehensive shock tube calculation for a given shock Mach number.

This function calculates all relevant states in a shock tube experiment:

Required driver pressure to achieve the specified Mach number
Post-shock conditions in the driven section (region 2)
Reflected shock conditions at the end wall (region 5)
Wave velocities and gas velocities

The mutating version shockcalc! modifies the input driver object in place, while shockcalc creates a copy and leaves the original unchanged.

Parameters

driver: A Species or Mixture object representing the driver gas (high-pressure section)
driven: A Species or Mixture object representing the driven gas (low-pressure test section)
Ms: Desired incident shock Mach number

Returns

A ShockCalc struct containing:

All four gas states (driver, driven, shocked, reflected)
Incident and reflected shock Mach numbers
Post-shock gas velocity

The returned object displays as a formatted table with pressure, temperature, density, and sound speed for each region.

Examples

julia> driver = Species("He")
Species(He, 298.1 K, 1.013e+05 Pa)

julia> driven = Mixture(["He" => 0.95, "acetone" => 0.05], T = 18u"°C", P = 85u"kPa")
Mixture(95% He, 5% acetone, 291.1 K, 8.500e+04 Pa)

julia> result = shockcalc(driver, driven, 2.2);

julia> round(density(result.shocked) / density(result.driven), digits=2)
2.65

Physics

The calculation follows standard shock tube theory:

Determines driver pressure needed for specified shock Mach number
Calculates incident shock jump conditions using Rankine-Hugoniot relations
Computes reflected shock conditions at the end wall
Returns complete state information for analysis

See Also

shockjump!: Calculate shock jump conditions only
driverpressure!: Calculate required driver pressure only
ShockCalc: Result container struct

source

PyThermo.ShockTube.ShockCalc — Type

ShockCalc

Container struct for comprehensive shock tube calculation results.

This struct stores the complete state information for all regions in a shock tube experiment, including the driver section, driven section, post-shock (region 2), and reflected shock regions.

Fields

driver::Chemical: Driver gas state (high-pressure section, region 4)
driven::Chemical: Initial driven gas state (low-pressure test section, region 1)
shocked::Chemical: Post-incident-shock gas state (region 2)
reflected::Chemical: Post-reflected-shock gas state (region 5)
Ms::Float64: Incident shock Mach number
Mr::Float64: Reflected shock Mach number
u2::Unitful.Velocity: Post-shock gas velocity (region 2)

Display

When displayed, ShockCalc presents results in a formatted table showing pressure, temperature, density, and sound speed for each region, along with incident and reflected wave velocities.

Examples

julia> driver = Species("He");

julia> driven = Mixture(["He" => 0.95, "acetone" => 0.05], T = 18u"°C", P = 85u"kPa");

julia> result = shockcalc(driver, driven, 2.2);

julia> result.Ms
2.2

julia> round(typeof(1.0u"kPa"), pressure(result.shocked), digits=1)
481.8 kPa

See Also

shockcalc!: Function that creates ShockCalc objects

source

PyThermo.ShockTube.riemann_interface — Function

riemann_interface(left::Chemical, right::Chemical, u_left=0.0u"m/s", u_right=0.0u"m/s") -> RiemannSolution
riemann_interface(left::Chemical, right::Chemical, Ms::Real) -> RiemannSolution

Solve the exact Riemann problem for two gases in contact at an interface.

This function uses the exact Riemann solver to determine the post-shock state at a gas interface, including pressure, velocity, densities, temperatures, and wave speeds for both gases. The solution accounts for different isentropic exponents (γ) on each side of the interface.

Parameters

Method 1: Direct specification of states and velocities

left::Chemical: Left gas state (Species or Mixture)
right::Chemical: Right gas state (Species or Mixture)
u_left::Unitful.Velocity: Velocity of left gas (default 0.0u"m/s")
u_right::Unitful.Velocity: Velocity of right gas (default 0.0u"m/s")

Method 2: Shock tube convenience method

left::Chemical: Unshocked driven gas (will be passed through shock)
right::Chemical: Right gas state (Species or Mixture, assumed at rest)
Ms::Real: Incident shock Mach number

When using Method 2, the function automatically:

Applies shock jump conditions to the left gas at the given Mach number
Computes the post-shock velocity
Solves the Riemann problem with the shocked left state
Assumes right gas is at rest (u_right = 0)

Returns

A RiemannSolution struct containing:

Interface pressure and velocity
Star region densities and temperatures for both gases
Wave speeds (left wave, right wave, contact discontinuity)
Isentropic exponents for reference

Physics

The exact Riemann solver:

Iteratively solves for the star region pressure using Newton-Raphson
Computes the star region velocity from momentum conservation
Calculates star region densities from Rankine-Hugoniot relations (shock) or isentropic relations (rarefaction)
Determines temperatures from the ideal gas law
Computes all wave speeds

Examples

Direct specification of states:

julia> left = Species("N2", P=500u"kPa", T=600u"K");

julia> right = Species("SF6", P=100u"kPa", T=300u"K");

julia> sol = riemann_interface(left, right);

julia> round(ustrip(u"kPa", sol.p_star), digits=1)
319.6

julia> round(ustrip(u"m/s", sol.u_star), digits=1)
155.8

With velocities:

julia> left = Species("N2", P=200u"kPa", T=400u"K");

julia> right = Species("Ar", P=100u"kPa", T=300u"K");

julia> sol = riemann_interface(left, right, 100.0u"m/s", -50.0u"m/s");

julia> sol.u_star > 0.0u"m/s"
true

Shock tube convenience method:

julia> driven = Species("N2");

julia> test_gas = Species("SF6");

julia> sol = riemann_interface(driven, test_gas, 1.5);

julia> round(ustrip(u"m/s", sol.S_contact), digits=1)
159.2

See Also

RiemannSolution: Result container struct
shockjump!: Calculate shock jump conditions

source

riemann_interface(left::Chemical, right::Chemical, Ms::Real) -> RiemannSolution

Convenience method for shock tube problems: automatically applies shock jump conditions to the left gas and solves the Riemann problem.

This method assumes:

left is the unshocked driven gas
right is at rest (u_right = 0)
A shock of Mach number Ms passes through the left gas

The function computes the shocked state and velocity of the left gas, then solves the exact Riemann problem at the interface.

Parameters

left::Chemical: Unshocked driven gas (will be passed through shock at Ms)
right::Chemical: Test gas at rest
Ms::Real: Incident shock Mach number

Returns

A RiemannSolution with the interface solution

Examples

julia> driven = Species("N2");

julia> test_gas = Species("SF6");

julia> sol = riemann_interface(driven, test_gas, 1.5);

julia> round(ustrip(u"m/s", sol.S_contact), digits=1)
159.2

source

PyThermo.ShockTube.RiemannSolution — Type

RiemannSolution

Container struct for exact Riemann solver results at a gas interface.

This struct stores the complete solution of the Riemann problem for two gases in contact, including the star region properties (interface conditions), star densities for both gases, and all relevant wave speeds. All properties can be accessed using dot notation.

Fields

p_star::Unitful.Pressure: Pressure at the interface (star region)
u_star::Unitful.Velocity: Velocity at the interface (star region)
rho_star_L::Unitful.Density: Density of left gas in star region
rho_star_R::Unitful.Density: Density of right gas in star region
T_star_L::Unitful.Temperature: Temperature of left gas in star region
T_star_R::Unitful.Temperature: Temperature of right gas in star region
S_L::Unitful.Velocity: Left wave speed (shock or rarefaction head)
S_R::Unitful.Velocity: Right wave speed (shock or rarefaction head)
S_contact::Unitful.Velocity: Contact discontinuity speed
gamma_L::Float64: Isentropic exponent of left gas
gamma_R::Float64: Isentropic exponent of right gas

Examples

julia> left = Species("N2", P=500u"kPa", T=600u"K");

julia> right = Species("SF6", P=100u"kPa", T=300u"K");

julia> sol = riemann_interface(left, right);

julia> round(u"kPa", sol.p_star, digits=2)
319.59 kPa

julia> round(u"m/s", sol.u_star, digits=1)
155.8 m s^-1

See Also

shockcalc!: Function that creates ShockCalc objects

source

Theory background

Rankine-Hugoniot relations

The shock jump conditions relate pre-shock and post-shock states through conservation of mass, momentum, and energy:

Mass: ρ₁u₁ = ρ₂u₂
Momentum: p₁ + ρ₁u₁² = p₂ + ρ₂u₂²
Energy: h₁ + u₁²/2 = h₂ + u₂²/2

Where subscript 1 denotes pre-shock and subscript 2 denotes post-shock states.

Riemann problem

The Riemann problem considers the evolution of two initially separated gas regions. The solution consists of three waves:

A wave propagating into the left state
A contact discontinuity (material interface)
A wave propagating into the right state

The interface pressure p* and velocity u* are found by solving:

fL(p, u) = fR(p, u)

where f represents the relations across the left and right waves. The solver iteratively finds p* and u* that satisfy continuity of pressure and velocity at the interface.