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

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

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)

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)