Shocks and detonations module#

In this section, you will find the documentation of the routines implemented to solves processes that involve strong changes in dynamic pressure, such as steady state shock and detonation waves in either normal or oblique stream configurations within the limits of regular shock reflections.

ShockSolver#

class ShockSolver#

Bases: handle

The ShockSolver() class is used to solve shock waves problems under different flow configurations. The class provides methods to compute the post-shock state of the gas mixture, including incident and reflected shocks, oblique shocks, and shock polar diagrams.

The ShockSolver() object can be initialized as follows:

solver = ShockSolver('problemType', problemType, ...)

Here problemType represents the acronym of the problem to be solved (see below). Additional optional parameters can be provided to customize the solver’s behavior.

Problem types:

  • SHOCK_I: Incident shock

  • SHOCK_R: Reflected shock

  • SHOCK_OBLIQUE: Oblique shock

  • SHOCK_OBLIQUE_R: Oblique reflected shock

  • SHOCK_POLAR: Shock polar diagrams

  • SHOCK_POLAR_R: Shock polar diagrams for incident and reflected states

  • SHOCK_POLAR_LIMITRR: Shock polar within the limit of regular reflection

  • SHOCK_PRANDTL_MEYER: Prandtl-Meyer expansion wave

See also: Mixture(), EquilibriumSolver(), DetonationSolver(), solve(), solveArray(), report()

FLAG_CACHE = 'true'#

Flag to clear cache after calculations

FLAG_REPORT = 'false'#

Flag to print predefined plots

FLAG_RESULTS = 'true'#

Flag to print results

FLAG_TIME = 'true'#

Flag to print elapsed time

equilibriumSolver = None#

EquilibriumSolver object

itLimitRR = '10'#

Max number of iterations - limit of regular reflections

itMax = '50'#

Max number of iterations - shocks and detonations

itOblique = '20'#

Max number of iterations - oblique shocks

machThermo = '2'#

Pre-shock Mach number above which T2_guess will be computed considering h2 = h1 + u1^2 / 2

numPointsPolar = '100'#

Number of points to compute shock/detonation polar curves

numPointsPrandtlMeyer = '100'#

Number of points to compute Prandtl-Meyer curves

plotConfig = None#

PlotConfig object

problemType = None#

Problem type

time = None#

Elapsed time

tol0 = '1e-5'#

Tolerance of shocks/detonations kernel

tolLimitRR = '1e-4'#

Tolerance to calculate the limit of regular reflections

tolOblique = '1e-3'#

Tolerance oblique shocks algorithm

shockIncident(obj, mix1, u1, varargin)#

Compute pre-shock and post-shock states of a planar incident shock wave

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • u1 (float) – Pre-shock velocity [m/s]

Optional Args:

mix2 (Mixture): Properties of the mixture in the post-shock state (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2 (Mixture): Properties of the mixture in the post-shock state

Examples

  • [mix1, mix2] = shockIncident(ShockSolver(), mix1, u1)

  • [mix1, mix2] = shockIncident(ShockSolver(), mix1, u1, mix2)

shockIncidentIdeal(obj, gamma, M1)#

Compute jump conditions assuming a thermochemically frozen gas (calorically perfect gas)

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • gamma (float) – Adiabatic index [-]

  • M1 (float) – Pre-shock Mach number [-]

Returns:

Tuple containing

  • R (float): Density ratio [-]

  • P (float): Pressure ratio [-]

  • T (float): Temperature ratio [-]

  • Gammas (float): Rankine-Hugoniot slope parameter [-]

  • M1 (float): Pre-shock Mach number [-]

Example

[R, P, T, Gammas, M1] = shockIncidentIdeal(ShockSolver(), 1.4, 2.0)

shockObliqueBeta(obj, mix1, u1, beta, varargin)#

Compute pre-shock and post-shock states of an oblique shock wave given the wave angle (one solution)

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • u1 (float) – Pre-shock velocity [m/s]

  • beta (float) – Wave angle [deg] of the incident oblique shock

Optional Args:

mix2 (Mixture): Properties of the mixture in the post-shock state (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2 (Mixture): Properties of the mixture at the post-shock state

Examples

  • [mix1, mix2] = shockObliqueBeta(ShockSolver(), mix1, u1, beta)

  • [mix1, mix2] = shockObliqueBeta(ShockSolver(), mix1, u1, beta, mix2)

shockObliqueReflectedTheta(obj, mix1, u2, theta, mix2, varargin)#

Compute pre-shock and post-shock states of an oblique reflected shock wave given the deflection angle.

Two solutions:

  • Weak shock

  • Strong shock

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state of the incident shock

  • u2 (float) – Post-shock velocity [m/s] of the incident shock

  • theta (float) – Deflection angle [deg]

  • mix2 (Mixture) – Properties of the mixture in the post-shock state of the incident shock

Optional Args:

  • mix5_1 (Mixture): Properties of the mixture in the post-shock state of the reflected shock - weak shock (previous calculation)

  • mix5_2 (Mixture): Properties of the mixture in the post-shock state of the reflected shock - strong shock (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state of the incident shock

  • mix2 (Mixture): Properties of the mixture in the post-shock state of the incident shock

  • mix5_1 (Mixture): Properties of the mixture in the post-shock state of the reflected shock - weak shock

  • mix5_2 (Mixture): Properties of the mixture in the post-shock state of the reflected shock - strong shock

Examples

  • [mix1, mix2, mix5_1, mix5_2] = shockObliqueReflectedTheta(ShockSolver(), mix1, u2, theta, mix2)

  • [mix1, mix2, mix5_1, mix5_2] = shockObliqueReflectedTheta(ShockSolver(), mix1, u2, theta, mix2, mix5_1, mix5_2)

shockObliqueTheta(obj, mix1, u1, theta, varargin)#

Compute pre-shock and post-shock states of an oblique shock wave given the deflection angle.

Two solutions:

  • Weak shock

  • Strong shock

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • u1 (float) – Pre-shock velocity [m/s]

  • theta (float) – Deflection angle [deg]

Optional Args:

  • mix2_1 (Mixture): Properties of the mixture in the post-shock state - weak shock (previous calculation)

  • mix2_2 (Mixture): Properties of the mixture in the post-shock state - strong shock (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2_1 (Mixture): Properties of the mixture in the post-shock state - weak shock

  • mix2_2 (Mixture): Properties of the mixture in the post-shock state - strong shock

Examples

  • [mix1, mix2_1, mix2_2] = shockObliqueTheta(ShockSolver(), mix1, u1, theta)

  • [mix1, mix2_1, mix2_2] = shockObliqueTheta(ShockSolver(), mix1, u1, theta, mix2_1, mix2_2)

shockPolar(obj, mix1, u1)#

Compute shock polar diagrams

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • u1 (float) – Pre-shock velocity [m/s]

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2 (Mixture): Properties of the mixture at the post-shock state with the shock polar results

Example

  • [mix1, mix2] = shockPolar(ShockSolver(), mix1, u1)

shockPolarLimitRR(obj, mix1, u1)#

Obtain polar curves for the given pre-shock conditions using Broyden’s method

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • u1 (float) – Pre-shock velocity [m/s]

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2 (Mixture): Properties of the mixture in the post-shock state - polar diagrams from mix1 (incident)

  • mix2_1 (Mixture): Properties of the mixture in the post-shock state - weak shock

  • mix3 (Mixture): Properties of the mixture in the post-shock state - polar diagrams from mix2_1 (reflected)

Example

[mix1, mix2, mix2_1, mix3] = shockPolarLimitRR(ShockSolver(), mix1, u1)

shockPrandtlMeyer(obj, mix1, u1, theta2, varargin)#

Computes the downstream state associated with an isentropic expansion through a given deflection angle theta2 by marching along the SV (defined entropy and specific volume) expansion path. The objective function thetaGuess - theta2 = 0 is solved using a Newton-Raphson method

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-expansion state

  • u1 (float) – Pre-expansion velocity [m/s]

  • theta2 (float) – Deflection angle [deg]

Optional Args:

  • mix2 (Mixture): Properties of the mixture in the post-expansion state (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the initial state

  • mix2 (Mixture): Properties of the final state after expansion

Notes

  • mix2.polar stores the path history (states along the expansion)

Example

[mix1, mix2] = shockPrandtlMeyer(ShockSolver(), mix1, 15.0)

shockReflected(obj, mix1, mix2, varargin)#

Compute pre-shock and post-shock states of a planar reflected shock wave

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:

  • obj (ShockSolver) – ShockSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • mix2 (Mixture) – Properties of the mixture at the post-shock state of the incident shock

Optional Args:

mix5 (Mixture): Properties of the mixture in the post-shock state of the reflected shock (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state of the incident shock

  • mix2 (Mixture): Properties of the mixture at the post-shock state of the incident shock

  • mix5 (Mixture): Properties of the mixture in the post-shock state of the reflected shock

Examples

  • [mix1, mix2, mix5] = shockReflected(ShockSolver(), mix1, mix2)

  • [mix1, mix2, mix5] = shockReflected(ShockSolver(), mix1, mix2, mix5)

DetonationSolver#

class DetonationSolver(varargin)#

Bases: handle

The DetonationSolver() class is used to solve detonation waves problems under different flow configurations. The class provides methods to compute the post-detonation state of the gas mixture, including incident and reflected detonations in over-driven and under-driven conditions, oblique detonations, and detonation polar diagrams.

The DetonationSolver() object can be initialized as follows:

solver = DetonationSolver('problemType', problemType, ...)

Here problemType represents the acronym of the problem to be solved (see below). Additional optional parameters can be provided to customize the solver’s behavior.

Problem types:

  • DET: Chapman-Jouguet detonation

  • DET_R: Chapman-Jouguet reflected detonation

  • DET_OVERDRIVEN: Over-driven detonation

  • DET_OVERDRIVEN_R: Reflected over-driven detonation

  • DET_UNDERDRIVEN: Under-driven detonation

  • DET_UNDERDRIVEN_R: Reflected under-driven detonation

  • DET_OBLIQUE: Oblique detonation

  • DET_OBLIQUE_R: Oblique reflected detonation

  • DET_POLAR: Detonation polar diagrams

  • DET_POLAR_R: Detonation polar diagrams for incident and reflected states

See also: Mixture(), EquilibriumSolver(), ShockSolver(), solve(), solveArray(), report()

DetonationSolver(varargin)#

Constructor

FLAG_CACHE = 'true'#

Flag to clear cache after calculations

FLAG_REPORT = 'false'#

Flag to print predefined plots

FLAG_RESULTS = 'true'#

Flag to print results

FLAG_TIME = 'true'#

Flag to print elapsed time

equilibriumSolver = None#

EquilibriumSolver object

itGuess = '5'#

Max number of iterations - guess detonation

itLimitRR = '10'#

Max number of iterations - limit of regular reflections

itMax = '50'#

Max number of iterations - shocks and detonations

itOblique = '20'#

Max number of iterations - oblique shocks

machThermo = '2'#

Pre-shock Mach number above which T2_guess will be computed considering h2 = h1 + u1^2 / 2

numPointsPolar = '100'#

Number of points to compute shock/detonation polar curves

plot(mixArray1, mixArray2, varargin)#

Plot results

Parameters:

  • obj (DetonationSolver) – DetonationSolver object

  • mixArray1 (Mixture) – Array of Mixture objects (pre-shock state)

  • mixArray2 (Mixture) – Array of Mixture objects (post-shock state)

Optional Args:

  • mixArray_i (Mixture): Array of Mixture objects

Examples

  • plot(DetonationSolver(), mixArray1, mixArray2);

  • plot(DetonationSolver(), mixArray1, mixArray2, mixArray3);

plotConfig = None#

PlotConfig object

printTime()#

Print execution time

Parameters:

obj (EquilibriumSolver) – Object of the class EquilibriumSolver

problemType = None#

Problem type

report(mixArray1, mixArray2, varargin)#

Postprocess all the results with predefined plots

Parameters:

  • obj (DetonationSolver) – DetonationSolver object

  • mixArray1 (Mixture) – Array of Mixture objects (pre-shock state)

  • mixArray2 (Mixture) – Array of Mixture objects (post-shock state)

Optional args:

  • mixArray_i (Mixture): Array of Mixture objects

Examples

  • report(DetonationSolver(), mixArray1, mixArray2);

  • report(DetonationSolver(), mixArray1, mixArray2, mixArray3);

static setMixtureShockPolar(mix, mixPolar)#

Assign values from the polar curves

shockSolver = None#

ShockSolver object

solve(mix1, varargin)#

Solve detonations problems

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mix1 (Mixture) – initial Mixture object

Optional Args:

  • mixGuess_i (Mixture): Mixture object from previous calculation

Returns:

varargout (Mixture) – Updated Mixture objects depending on the problem type

Examples

  • [mix1, mix2] = solve(DetonationSolver(), mix1);

  • [mix1, mix2] = solve(DetonationSolver(), mix1, mix2Guess);

  • [mix1, mix2, mix3] = solve(DetonationSolver(), mix1, mix2Guess, mix3Guess);

solveArray(mixArray1, varargin)#

Solve a set of detonation problems

Parameters:

  • obj (EquilibriumSolver) – EquilibriumSolver object

  • mixArray1 (Mixture) – Array of initial Mixture objects

Returns:

varargout (Mixture) – Updated arrays of Mixture objects depending on the shock problem type

Examples

  • [mixArray1, mixArray2] = solveArray(DetonationSolver(), mixArray1);

  • [mixArray1, mixArray2, mixArray3] = solveArray(DetonationSolver(), mixArray1);

time = None#

Elapsed time [s]

tol0 = '1e-5'#

Tolerance of shocks/detonations kernel

tolLimitRR = '1e-4'#

Tolerance to calculate the limit of regular reflections

tolOblique = '1e-3'#

Tolerance oblique shocks algorithm

detonationCJ(obj, mix1, varargin)#

Compute pre-shock and post-shock states of a Chapman-Jouguet detonation

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:

  • obj (DetonationSolver) – DetonationSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

Optional Args:

mix2 (Mixture): Properties of the mixture in the post-shock state (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2 (Mixture): Properties of the mixture in the post-shock state

Examples

  • [mix1, mix2] = detonationCJ(DetonationSolver(), mix1)

  • [mix1, mix2] = detonationCJ(DetonationSolver(), mix1, mix2)

detonationGuess(obj, mix1)#

Obtain guess of the jump conditions for a Chapman-Jouguet detonation as in NASA’s CEA code (see Sec. 8.3 of [1])

[1] Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:

  • obj (DetonationSolver) – DetonationSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

Returns:

Tuple containing

  • P (float): Pressure ratio [-]

  • T (float): Temperature ratio [-]

  • STOP (float): Relative error [-]

Example

[P, T, STOP] = detonationGuess(obj, mix1)

detonationObliqueBeta(obj, mix1, driveFactor, beta, varargin)#

Compute pre-shock and post-shock states of an oblique detonation wave given the wave angle.

Two solutions:

  • Underdriven detonation

  • Overdriven detonation

Parameters:

  • obj (DetonationSolver) – Detonation solver object

  • mix1 (struct) – Properties of the mixture in the pre-shock state

  • driveFactor (float) – Driven factor [-]

  • beta (float) – Wave angle [deg] of the incident oblique detonation

Optional Args:

mix2 (Mixture): Properties of the mixture in the post-shock state (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-detonation state

  • mix2_1 (Mixture): Properties of the mixture in the post-detonation state - underdriven detonation

  • mix2_2 (Mixture): Properties of the mixture in the post-detonation state - overdriven detonation

Examples

  • [mix1, mix2_1, mix2_2] = detonationObliqueBeta(DetonationSolver(), mix1, 2, 60)

  • [mix1, mix2_1, mix2_2] = detonationObliqueBeta(DetonationSolver(), mix1, 2, 60, mix2)

detonationObliqueTheta(obj, mix1, driveFactor, theta, varargin)#

Compute pre-shock and post-shock states of an oblique detonation wave given the deflection angle (only overdriven solution).

Two solutions:

  • Weak detonation

  • Strong detonation

Parameters:

  • obj (DetonationSolver) – Detonation solver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • driveFactor (float) – Drive factor [-]

  • theta (float) – Deflection angle [deg]

Optional Args:

  • mix2_1 (Mixture): Properties of the mixture in the post-shock state - weak detonation (previous calculation)

  • mix2_2 (Mixture): Properties of the mixture in the post-shock state - strong detonation (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-detonation state

  • mix2_1 (Mixture): Properties of the mixture in the post-detonation state - weak detonation

  • mix2_2 (Mixture): Properties of the mixture in the post-detonation state - strong detonation

Examples

  • [mix1, mix2_1, mix2_2] = detonationObliqueTheta(DetonationSolver(), mix1, 2, 30)

  • [mix1, mix2_1, mix2_2] = detonationObliqueTheta(DetonationSolver(), mix1, 2, 30, mix2)

detonationOverdriven(obj, mix1, driveFactor, varargin)#

Compute pre-shock and post-shock states of an overdriven planar detonation

Parameters:

  • obj (DetonationSolver) – DetonationSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • driveFactor (float) – Overdriven factor [-]

Optional Args:

mix2 (Mixture): Properties of the mixture in the post-shock state (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2 (Mixture): Properties of the mixture in the post-shock state

Examples

  • [mix1, mix2] = detonationOverdriven(DetonationSolver(), mix1, 1.5)

  • [mix1, mix2] = detonationOverdriven(DetonationSolver(), mix1, 1.5, mix2)

detonationPolar(obj, mix1, driveFactor, varargin)#

Compute detonation polar diagrams

Parameters:

  • obj (DetonationSolver) – DetonationSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • driveFactor (float) – Driven factor

Optional Args:

mix2 (Mixture): Properties of the mixture in the post-shock state (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2 (Mixture): Properties of the mixture at the post-shock state with the shock polar results

Examples

  • [mix1, mix2] = detonationPolar(DetonationSolver(), mix1, 3000)

  • [mix1, mix2] = detonationPolar(DetonationSolver(), mix1, 3000, mix2)

detonationReflected(obj, mix1, mix2, varargin)#

Compute pre-shock and post-shock states of a planar reflected shock wave

This method is based on the method outlined in Gordon, S., & McBride, B. J. (1994). NASA reference publication, 1311.

Parameters:

  • obj (DetonationSolver) – DetonationSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • mix2 (Mixture) – Properties of the mixture at the post-shock state of the incident shock

Optional Args:

mix5 (Mixture): Properties of the mixture in the post-shock state of the reflected shock (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state of the incident shock

  • mix2 (Mixture): Properties of the mixture at the post-shock state of the incident shock

  • mix5 (Mixture): Properties of the mixture in the post-shock state of the reflected shock

Examples

  • [mix1, mix2, mix5] = detonationReflected(DetonationSolver(), mix1, mix2)

  • [mix1, mix2, mix5] = detonationReflected(DetonationSolver(), mix1, mix2, mix5)

detonationUnderdriven(obj, mix1, driveFactor, varargin)#

Compute pre-shock and post-shock states of an overdriven planar detonation

Parameters:

  • obj (DetonationSolver) – DetonationSolver object

  • mix1 (Mixture) – Properties of the mixture in the pre-shock state

  • driveFactor (float) – Underdriven factor [-]

Optional Args:

mix2 (Mixture): Properties of the mixture in the post-shock state (previous calculation)

Returns:

Tuple containing

  • mix1 (Mixture): Properties of the mixture in the pre-shock state

  • mix2 (Mixture): Properties of the mixture in the post-shock state

Examples

  • [mix1, mix2] = detonationUnderdriven(DetonationSolver(), mix1, 1.5)

  • [mix1, mix2] = detonationUnderdriven(DetonationSolver(), mix1, 1.5, mix2)

JumpConditionsSolver#

class JumpConditionsSolver(varargin)#

Bases: handle

The JumpConditionsSolver() class computes the Rankine–Hugoniot jump conditions and associated amplification parameters for gas mixtures across a shock wave. It supports calorically perfect, thermally perfect, and calorically imperfect models. In addition to post-shock thermodynamic properties, the class evaluates the dimensionless slopes of the Hugoniot curve—Gammas2 (\(\Gamma\)), Gammas1 (\(\Gamma_\rho\)), and Gammas3 (\(\Gamma_p\))—which quantify the local sensitivity of the post-shock density to changes in downstream pressure (\(\Gamma\)), and the influence of upstream density and pressure perturbations on the post-shock state (\(\Gamma_\rho\), \(\Gamma_p\)), as formalized in Ref. [1, 2].

The JumpConditionsSolver() object can be initialized as follows:

solver = JumpConditionsSolver('equilibriumSolver', eqSolver, ...)

Optional name-value arguments can be specified to customize tolerances, interpolation settings, plotting options, and solver behavior. Default objects for EquilibriumSolver() and ShockSolver() are created automatically if not provided.

Flags:

  • FLAG_FAST: Use previous state as initial guess for equilibrium solver

  • FLAG_PLOT: Plot results

  • FLAG_PAPER: Apply the normalization and notation conventions described in Cuadra2025a

  • FLAG_RESULTS: Print computed values to console

  • FLAG_TIME: Print elapsed computation time

  • FLAG_REPORT: Generate predefined plots for diagnostic analysis

See also: ShockSolver(), EquilibriumSolver(), solve(), report(), plot(), plotGammas()

References

[1] Cuadra, A., Di Renzo, M., Hoste, J. J. O., Williams, C. T., Vera, M., & Huete, C. (2025).

Review of shock-turbulence interaction with a focus on hypersonic flow. Physics of Fluids, 37(4). DOI: 10.1063/5.0255816.

[2] Cuadra, A., Williams, C. T., Di Renzo, M. & Huete, C. (2024). Compressibility

and vibrational-excitation effects in hypersonic shock-turbulence interaction. Tech. Rep. Summer Program Proceedings, Center for Turbulence Research, Stanford University.

FLAG_PAPER = 'false'#

Flag to apply the normalization and notation conventions described in Cuadra2025a

FLAG_PLOT = 'false'#

Flag to plot results

FLAG_REPORT = 'false'#

Flag to print predefined plots

FLAG_RESULTS = 'true'#

Flag to print results

FLAG_TIME = 'true'#

Flag to print elapsed time

JumpConditionsSolver(varargin)#

Constructor

Optional Args:

  • param1 (type): Description of param1

  • param2 (type): Description of param2

Returns:

obj (JumpConditionsSolver) – Instance of the solver

Example

obj = JumpConditionsSolver(‘param1’, value1, ‘param2’, value2, …);

dp2dM1 = None#

Derivative of p2 with respect to M1

drho2dM1 = None#

Derivative of rho2 with respect to M1

equilibriumSolver = None#

EquilibriumSolver object

static getDensityRatioPerfect(gamma, mach)#

Get the density ratio for a calorically perfect gas

Parameters:

  • gamma (float) – Adiabatic index

  • mach (float) – mach number

Returns:

Rratio (float) – Density ratio

getGammas1(jumpConditions)#

Get the dimensionless slope of the Hugoniot curve Gammas1

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • jumpConditions (struct) – Structure containing the jump conditions

Returns:

Gammas1 (float) – Dimensionless slope of the Hugoniot curve (partial derivative at constant rho2, p1)

getGammas2(jumpConditions)#

Get the dimensionless slope of the Hugoniot curve Gammas2

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • jumpConditions (struct) – Structure containing the jump conditions

Returns:

Gammas2 (float) – Dimensionless slope of the Hugoniot curve (partial derivative at constant rho1, p1)

getGammas3(jumpConditions)#

Get the dimensionless slope of the Hugoniot curve Gammas3

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • jumpConditions (struct) – Structure containing the jump conditions

Returns:

Gammas3 (float) – Dimensionless slope of the Hugoniot curve (partial derivative at constant rho1, rho2)

getGammasPerfect(gamma, mach)#

Get the dimensionless slope of the Hugoniot curve for a calorically perfect gas

Parameters:

  • gamma (float) – Adiabatic index

  • mach (float) – mach number

Returns:

Tuple containing –

  • Gammas1 (float): Dimensionless slope of the Hugoniot curve (partial derivative at constant rho2, p1)

  • Gammas2 (float): Dimensionless slope of the Hugoniot curve (partial derivative at constant rho1, p1)

  • Gammas3 (float): Dimensionless slope of the Hugoniot curve (partial derivative at constant rho1, rho2)

getJumpData(mixArray1)#

Get jump conditions

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • mixArray1 (Mixture) – Mixture object for the pre-shock state

Returns:

Tuple containing –

  • jumpConditions (struct): Structure containing the jump conditions

  • mixArray1 (Mixture): Mixture object for the pre-shock state

  • mixArray2 (Mixture): Mixture object for the post-shock state

static getPressureRatioPerfect(gamma, mach)#

Get the pressure ratio for a calorically perfect gas

Parameters:

  • gamma (float) – Adiabatic index

  • mach (float) – mach number

Returns:

Pratio (float) – Pressure ratio

mach = None#

Mach number array of the array of Mixture objects

mix = None#

Mixture object for initialization

plot(jumpConditions)#

Plot results

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • jumpConditions (struct) – Structure containing the jump conditions

Example

  • plot(JumpConditionsSolver(), results);

plotGammas(jumpConditions)#

Plot dimensionless slopes of the Hugoniot curve

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • jumpConditions (struct) – Structure containing the jump conditions

Returns:

ax (axes) – Axes handle for the plot

Example

plotGammas(JumpConditionsSolver(), jumpConditions);

pressure = None#

Pressure array of the array of Mixture objects

printTime()#

Print execution time

Parameters:

obj (EquilibriumSolver) – Object of the class EquilibriumSolver

report(jumpConditions)#

Generate report for the jump conditions

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • jumpConditions (struct) – Structure containing the jump conditions

Example

report(JumpConditionsSolver(), jumpConditions);

static reshapeVector(x, y, varargin)#

Reshape vectors to be column vectors of the form

[x(1), y(1), x(2), y(2), …]

Parameters:

  • x (float) – First vector

  • y (float) – Second vector

Optional Args:

  • varargin (float): Additional vectors

shockSolver = None#

ShockSolver object

solve(mixArray1, varargin)#

Solve the jump conditions for the given mixture

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • mixArray (Mixture) – Array of Mixture objects

Returns:

obj (JumpConditionsSolver) – JumpConditionsSolver object with updated properties

Example

jumpConditions = solve(JumpConditionsSolver(), mixArray);

solveShock(T1, p1)#

Solve shock problem for the given temperature and pressure

Parameters:

  • obj (JumpConditionsSolver) – JumpConditionsSolver object

  • T1 (float) – Temperature of the mixture

  • p1 (float) – Pressure of the mixture

Returns:

Tuple containing –

  • rho1 (float): Density of the mixture in the pre-shock state

  • rho2 (float): Density of the mixture in the post-shock state

  • p2 (float): Pressure of the mixture in the post-shock state

temperature = None#

Temperature array of the array of Mixture objects

time = None#

Elapsed time

tolGammas1 = '1e-4'#

Tolerance for the calculation of the dimensionless RH slope Gammas1

tolGammas2 = '1e-4'#

Tolerance for the calculation of the dimensionless RH slope Gammas2

tolGammas3 = '1e-4'#

Tolerance for the calculation of the dimensionless RH slope Gammas3