core.compressible_flow

machwave.core.compressible_flow

Compressible flow theory and analysis.

apply_thrust_coefficient_correction(ideal_thrust_coefficient, nozzle_correction_factor)

Apply nozzle efficiency correction to thrust coefficient.

Parameters:

Name	Type	Description	Default
ideal_thrust_coefficient	float	Ideal thrust coefficient.	required
nozzle_correction_factor	float	Nozzle efficiency (0-1).	required

Returns:

Type	Description
float	Corrected thrust coefficient.

Source code in machwave/core/compressible_flow/nozzle.py

def apply_thrust_coefficient_correction(
    ideal_thrust_coefficient: float,
    nozzle_correction_factor: float,
) -> float:
    """Apply nozzle efficiency correction to thrust coefficient.

    Args:
        ideal_thrust_coefficient: Ideal thrust coefficient.
        nozzle_correction_factor: Nozzle efficiency (0-1).

    Returns:
        Corrected thrust coefficient.
    """
    return ideal_thrust_coefficient * nozzle_correction_factor

get_boundary_layer_percentage_loss(chamber_pressure_psi, throat_diameter_inch, expansion_ratio, time, c_1, c_2)

Boundary layer correction factor accounts for the decrement in performance due to the viscous and heat transfer effects in the nozzle walls. It is time dependent.

Valid for liquid, solid, and hybrid propellants.

The time dependence is exponential due to the transient heat up, important in motors with short burn durations (less than 4 seconds). Dependence on expansion ratio represents the effect of a the amount of nozzle surface area.

Time constant C2 comes from analysis of a the transient heating of a BATES motor.

Time constant C1 was obtained from a direct measurement of the heat loss in a BATES motor, among other things.

Ordinary nozzle: C1 = 0.003650 C2 = 0.000937

Solid steel nozzle with relatively thick walls: C1 = 0.005060 C2 = 0.000000

Parameters:

Name	Type	Description	Default
chamber_pressure_psi	float	The chamber pressure [psi].	required
throat_diameter_inch	float	The throat diameter [in].	required
expansion_ratio	float	The expansion ratio of the nozzle.	required
time	float	The time in seconds [s].	required
c_1	float	Coefficient for the boundary layer correction factor.	required
c_2	float	Coefficient for the boundary layer correction factor.	required

Returns: The boundary layer correction factor.

Source code in machwave/core/compressible_flow/losses.py

@decorators.check_bounds(lower=0.0, upper=100.0)
@decorators.warn_if_outside_range(**TYPICAL_RANGES["boundary_layer_loss"])
def get_boundary_layer_percentage_loss(
    chamber_pressure_psi: float,
    throat_diameter_inch: float,
    expansion_ratio: float,
    time: float,
    c_1: float,
    c_2: float,
) -> float:
    """
    Boundary layer correction factor accounts for the decrement in performance due to
    the viscous and heat transfer effects in the nozzle walls. It is time dependent.

    Valid for liquid, solid, and hybrid propellants.

    The time dependence is exponential due to the transient heat up, important in
    motors with short burn durations (less than 4 seconds). Dependence on expansion
    ratio represents the effect of a the amount of nozzle surface area.

    Time constant C2 comes from analysis of a the transient heating of a BATES motor.

    Time constant C1 was obtained from a direct measurement of the heat loss in a BATES
    motor, among other things.

    Ordinary nozzle:
    C1 = 0.003650
    C2 = 0.000937

    Solid steel nozzle with relatively thick walls:
    C1 = 0.005060
    C2 = 0.000000

    Args:
        chamber_pressure_psi: The chamber pressure [psi].
        throat_diameter_inch: The throat diameter [in].
        expansion_ratio: The expansion ratio of the nozzle.
        time: The time in seconds [s].
        c_1: Coefficient for the boundary layer correction factor.
        c_2: Coefficient for the boundary layer correction factor.
    Returns:
        The boundary layer correction factor.
    """
    term_1 = c_1 * (chamber_pressure_psi**0.8) / (throat_diameter_inch**0.2)
    term_2 = 1 + 2 * np.exp(
        (-c_2 * chamber_pressure_psi**0.8 * time) / (throat_diameter_inch**0.2)
    )
    term_3 = 1 + 0.016 * (expansion_ratio - 9)

    return term_1 * term_2 * term_3

get_critical_pressure_ratio(k)

Get critical pressure ratio for choked flow.

Parameters:

Name	Type	Description	Default
k	float	Isentropic exponent.	required

Returns:

Type	Description
float	Critical pressure ratio.

Source code in machwave/core/compressible_flow/isentropic.py

def get_critical_pressure_ratio(k: float) -> float:
    """Get critical pressure ratio for choked flow.

    Args:
        k: Isentropic exponent.

    Returns:
        Critical pressure ratio.
    """
    return (2 / (k + 1)) ** (k / (k - 1))

get_exit_mach_from_expansion_ratio(k, expansion_ratio)

Get exit Mach number from expansion ratio.

Parameters:

Name	Type	Description	Default
k	float	Isentropic exponent.	required
expansion_ratio	float	Expansion ratio (A_exit / A_throat).	required

Returns:

Type	Description
float	Exit Mach number.

Raises:

Type	Description
ValueError	If solver fails to converge.

Source code in machwave/core/compressible_flow/isentropic.py

def get_exit_mach_from_expansion_ratio(k: float, expansion_ratio: float) -> float:
    """Get exit Mach number from expansion ratio.

    Args:
        k: Isentropic exponent.
        expansion_ratio: Expansion ratio (A_exit / A_throat).

    Returns:
        Exit Mach number.

    Raises:
        ValueError: If solver fails to converge.
    """
    try:
        exit_mach = cast(
            float,
            scipy.optimize.brentq(
                lambda m: get_expansion_ratio_from_exit_mach(m, k) - expansion_ratio,
                a=1.001,  # Just above sonic
                b=20.0,  # High supersonic
            ),
        )
        return exit_mach
    except ValueError as e:
        raise ValueError(
            f"Failed to converge for expansion_ratio={expansion_ratio}, k={k}"
        ) from e

get_exit_pressure(k_ex, expansion_ratio, chamber_pressure)

Get exit pressure from isentropic relations.

Parameters:

Name	Type	Description	Default
k_ex	float	Isentropic exponent at exit.	required
expansion_ratio	float	Expansion ratio.	required
chamber_pressure	float	Chamber pressure [Pa].	required

Returns:

Type	Description
float	Exit pressure [Pa].

Source code in machwave/core/compressible_flow/isentropic.py

def get_exit_pressure(
    k_ex: float, expansion_ratio: float, chamber_pressure: float
) -> float:
    """Get exit pressure from isentropic relations.

    Args:
        k_ex: Isentropic exponent at exit.
        expansion_ratio: Expansion ratio.
        chamber_pressure: Chamber pressure [Pa].

    Returns:
        Exit pressure [Pa].
    """
    exit_mach = get_exit_mach_from_expansion_ratio(k_ex, expansion_ratio)
    return chamber_pressure * (1 + 0.5 * (k_ex - 1) * exit_mach**2) ** (
        -k_ex / (k_ex - 1)
    )

get_expansion_ratio_from_exit_mach(mach, k)

Get expansion ratio from exit Mach number.

Parameters:

Name	Type	Description	Default
mach	float	Mach number.	required
k	float	Isentropic exponent.	required

Returns:

Type	Description
float	Expansion ratio (A / A_throat).

Source code in machwave/core/compressible_flow/isentropic.py

def get_expansion_ratio_from_exit_mach(mach: float, k: float) -> float:
    """Get expansion ratio from exit Mach number.

    Args:
        mach: Mach number.
        k: Isentropic exponent.

    Returns:
        Expansion ratio (A / A_throat).
    """
    term1 = (2 / (k + 1)) * (1 + 0.5 * (k - 1) * mach**2)
    term2 = (k + 1) / (2 * (k - 1))
    return (1 / mach) * (term1**term2)

get_ideal_thrust_coefficient(chamber_pressure, exit_pressure, external_pressure, expansion_ratio, k_ex)

Get ideal thrust coefficient for DeLaval nozzle.

Parameters:

Name	Type	Description	Default
chamber_pressure	float	Chamber pressure [Pa].	required
exit_pressure	float	Exit pressure [Pa].	required
external_pressure	float	External pressure [Pa].	required
expansion_ratio	float	Expansion ratio.	required
k_ex	float	Isentropic exponent at exit.	required

Returns:

Type	Description
float	Ideal thrust coefficient.

References

https://www.nakka-rocketry.net/th_thrst.html

Source code in machwave/core/compressible_flow/nozzle.py

def get_ideal_thrust_coefficient(
    chamber_pressure: float,
    exit_pressure: float,
    external_pressure: float,
    expansion_ratio: float,
    k_ex: float,
) -> float:
    """Get ideal thrust coefficient for DeLaval nozzle.

    Args:
        chamber_pressure: Chamber pressure [Pa].
        exit_pressure: Exit pressure [Pa].
        external_pressure: External pressure [Pa].
        expansion_ratio: Expansion ratio.
        k_ex: Isentropic exponent at exit.

    Returns:
        Ideal thrust coefficient.

    References:
        https://www.nakka-rocketry.net/th_thrst.html
    """
    pressure_ratio = exit_pressure / chamber_pressure
    return (
        np.sqrt(
            (2 * (k_ex**2) / (k_ex - 1))
            * ((2 / (k_ex + 1)) ** ((k_ex + 1) / (k_ex - 1)))
            * (1 - (pressure_ratio ** ((k_ex - 1) / k_ex)))
        )
        + expansion_ratio * (exit_pressure - external_pressure) / chamber_pressure
    )

get_kinetics_percentage_loss(i_sp_th_frozen, i_sp_th_shifting, chamber_pressure_psi)

The kinetics correction factor accounts for the decrement in performance due to incomplete heat transfer of latent heat to sensible heat caused by the finite time required for the gas phase chemical reactions to occur.

Valid for liquid, solid, and hybrid propellants. The expansion ratio of the i_sp_th_frozen and i_sp_th_shifting should be the same.

Pressure correction is applied for chamber pressures above 1.379 MPa (200 psi), in order to dampen the effect of the kinetics correction factor.

Parameters:

Name	Type	Description	Default
i_sp_th_frozen	float	The specific impulse of the frozen flow [s].	required
i_sp_th_shifting	float	The specific impulse of the shifting flow [s].	required
chamber_pressure_psi	float	The chamber pressure [psi].	required

Returns: The kinetics correction factor.

Source code in machwave/core/compressible_flow/losses.py

@decorators.check_bounds(lower=0.0, upper=100.0)
@decorators.warn_if_outside_range(**TYPICAL_RANGES["kinetics_loss"])
def get_kinetics_percentage_loss(
    i_sp_th_frozen: float, i_sp_th_shifting: float, chamber_pressure_psi: float
) -> float:
    """
    The kinetics correction factor accounts for the decrement in performance due to
    incomplete heat transfer of latent heat to sensible heat caused by the finite time
    required for the gas phase chemical reactions to occur.

    Valid for liquid, solid, and hybrid propellants.
    The expansion ratio of the i_sp_th_frozen and i_sp_th_shifting should be the same.

    Pressure correction is applied for chamber pressures above 1.379 MPa (200 psi), in
    order to dampen the effect of the kinetics correction factor.

    Args:
        i_sp_th_frozen: The specific impulse of the frozen flow [s].
        i_sp_th_shifting: The specific impulse of the shifting flow [s].
        chamber_pressure_psi: The chamber pressure [psi].
    Returns:
        The kinetics correction factor.
    """
    i_sp_th_ratio = i_sp_th_frozen / i_sp_th_shifting

    if chamber_pressure_psi < KINETICS_LOSS_PRESSURE_THRESHOLD_PSI:
        pressure_correction = 1.0
    else:
        pressure_correction = (
            KINETICS_LOSS_PRESSURE_THRESHOLD_PSI / chamber_pressure_psi
        )

    return 33.3 * (1 - i_sp_th_ratio) * pressure_correction

get_nozzle_divergent_percentage_loss(divergent_angle)

Calculates the divergent nozzle correction factor given the half angle. NOTE: only applicable for a conical convergent-divergent nozzle.

Parameters:

Name	Type	Description	Default
divergent_angle	float	The half angle of the divergent nozzle.	required

Returns:

Type	Description
float	The divergent correction factor.

Source code in machwave/core/compressible_flow/losses.py

@decorators.check_bounds(lower=0.0, upper=100.0)
@decorators.warn_if_outside_range(**TYPICAL_RANGES["divergent_loss"])
def get_nozzle_divergent_percentage_loss(divergent_angle: float) -> float:
    """
    Calculates the divergent nozzle correction factor given the half angle.
    NOTE: only applicable for a conical convergent-divergent nozzle.

    Args:
        divergent_angle: The half angle of the divergent nozzle.

    Returns:
        The divergent correction factor.
    """
    return 50 * (1 - np.cos(np.deg2rad(divergent_angle)))

get_optimal_expansion_ratio(k, chamber_pressure, atmospheric_pressure)

Get optimal expansion ratio for DeLaval nozzle.

Parameters:

Name	Type	Description	Default
k	float	Isentropic exponent.	required
chamber_pressure	float	Chamber pressure [Pa].	required
atmospheric_pressure	float	External pressure [Pa].	required

Returns:

Type	Description
float	Optimal expansion ratio.

Source code in machwave/core/compressible_flow/nozzle.py

def get_optimal_expansion_ratio(
    k: float, chamber_pressure: float, atmospheric_pressure: float
) -> float:
    """Get optimal expansion ratio for DeLaval nozzle.

    Args:
        k: Isentropic exponent.
        chamber_pressure: Chamber pressure [Pa].
        atmospheric_pressure: External pressure [Pa].

    Returns:
        Optimal expansion ratio.
    """
    return (
        (((k + 1) / 2) ** (1 / (k - 1)))
        * ((atmospheric_pressure / chamber_pressure) ** (1 / k))
        * np.sqrt(
            ((k + 1) / (k - 1))
            * (1 - (atmospheric_pressure / chamber_pressure) ** ((k - 1) / k))
        )
    ) ** -1

get_overall_nozzle_efficiency(eta_div, eta_kin, eta_bl, eta_2p, other_losses=OTHER_LOSSES_DEFAULT)

Calculates the overall nozzle efficiency by combining the correction factors.

Parameters:

Name	Type	Description	Default
eta_div	float	The divergent nozzle correction factor.	required
eta_kin	float	The kinetics correction factor.	required
eta_bl	float	The boundary layer correction factor.	required
eta_2p	float	The two-phase flow correction factor.	required

Returns:

Type	Description
float	The overall nozzle efficiency.

Source code in machwave/core/compressible_flow/losses.py

@decorators.check_bounds(lower=0.0, upper=1.0)
def get_overall_nozzle_efficiency(
    eta_div: float,
    eta_kin: float,
    eta_bl: float,
    eta_2p: float,
    other_losses: float = OTHER_LOSSES_DEFAULT,
) -> float:
    """
    Calculates the overall nozzle efficiency by combining the correction factors.

    Args:
        eta_div: The divergent nozzle correction factor.
        eta_kin: The kinetics correction factor.
        eta_bl: The boundary layer correction factor.
        eta_2p: The two-phase flow correction factor.

    Returns:
        The overall nozzle efficiency.
    """
    return 1 - (eta_div + eta_kin + eta_bl + eta_2p + other_losses) / 100

get_thrust_from_thrust_coefficient(thrust_coefficient, chamber_pressure, nozzle_throat_area)

Get thrust from thrust coefficient.

Parameters:

Name	Type	Description	Default
thrust_coefficient	float	Thrust coefficient.	required
chamber_pressure	float	Chamber stagnation pressure [Pa].	required
nozzle_throat_area	float	Nozzle throat area [m^2].	required

Returns:

Type	Description
float	Thrust [N].

Source code in machwave/core/compressible_flow/nozzle.py

def get_thrust_from_thrust_coefficient(
    thrust_coefficient: float, chamber_pressure: float, nozzle_throat_area: float
) -> float:
    """Get thrust from thrust coefficient.

    Args:
        thrust_coefficient: Thrust coefficient.
        chamber_pressure: Chamber stagnation pressure [Pa].
        nozzle_throat_area: Nozzle throat area [m^2].

    Returns:
        Thrust [N].
    """
    return thrust_coefficient * chamber_pressure * nozzle_throat_area

get_two_phase_flow_percentage_loss(chamber_pressure_psi, mole_fraction_of_condensed_phase, expansion_ratio, throat_diameter_inch, characteristic_length_inch)

Two-phase flow correction factor accounts for the decrement in performance due to the presence of a condensed phase in the combustion products.

Valid for solid, and hybrid propellants.

Parameters:

Name	Type	Description	Default
chamber_pressure_psi	float	The chamber pressure [psi].	required
mole_fraction_of_condensed_phase	float	The mole fraction of the condensed phase [moles/100g].	required
expansion_ratio	float	The expansion ratio of the nozzle.	required
throat_diameter_inch	float	The throat diameter [in].	required
characteristic_length_inch	float	The characteristic length [in].	required

Returns: The two-phase flow correction factor.

Source code in machwave/core/compressible_flow/losses.py

@decorators.check_bounds(lower=0.0, upper=100.0)
@decorators.warn_if_outside_range(**TYPICAL_RANGES["two_phase_flow_loss"])
def get_two_phase_flow_percentage_loss(
    chamber_pressure_psi: float,
    mole_fraction_of_condensed_phase: float,
    expansion_ratio: float,
    throat_diameter_inch: float,
    characteristic_length_inch: float,
) -> float:
    """
    Two-phase flow correction factor accounts for the decrement in performance due to
    the presence of a condensed phase in the combustion products.

    Valid for solid, and hybrid propellants.

    Args:
        chamber_pressure_psi: The chamber pressure [psi].
        mole_fraction_of_condensed_phase: The mole fraction of the condensed phase
            [moles/100g].
        expansion_ratio: The expansion ratio of the nozzle.
        throat_diameter_inch: The throat diameter [in].
        characteristic_length_inch: The characteristic length [in].
    Returns:
        The two-phase flow correction factor.
    """
    particle_size_um: float = _get_two_phase_phase_loss_particle_size(
        chamber_pressure_psi,
        mole_fraction_of_condensed_phase,
        throat_diameter_inch,
        characteristic_length_inch,
    )
    xi: float = mole_fraction_of_condensed_phase

    if xi >= 0.09:
        c_4 = 0.5
        if throat_diameter_inch < 1.0:
            c_3, c_5, c_6 = 9.0, 1.0, 1.0
        elif throat_diameter_inch < 2.0:
            c_3, c_5, c_6 = 9.0, 1.0, 0.8
        else:  # throat_diameter_inch >= 2
            if particle_size_um < 4.0:
                c_3, c_5, c_6 = 13.4, 0.8, 0.8
            elif particle_size_um <= 8.0:
                c_3, c_5, c_6 = 10.2, 0.8, 0.4
            else:
                c_3, c_5, c_6 = 7.58, 0.8, 0.33
    else:  # xi < 0.09
        c_4 = 1.0
        if throat_diameter_inch < 1.0:
            c_3, c_5, c_6 = 30.0, 1.0, 1.0
        elif throat_diameter_inch < 2.0:
            c_3, c_5, c_6 = 30.0, 1.0, 0.8
        else:  # throat_diameter_inch >= 2
            if particle_size_um < 4.0:
                c_3, c_5, c_6 = 44.5, 0.8, 0.8
            elif particle_size_um <= 8.0:
                c_3, c_5, c_6 = 34.0, 0.8, 0.4
            else:
                c_3, c_5, c_6 = 25.2, 0.8, 0.33

    numerator = (xi**c_4) * (particle_size_um**c_5)
    denominator = (
        (chamber_pressure_psi**0.15)
        * (expansion_ratio**0.08)
        * (throat_diameter_inch**c_6)
    )

    return c_3 * numerator / denominator

is_flow_choked(chamber_pressure, external_pressure, critical_pressure_ratio)

Check if flow is choked.

Parameters:

Name	Type	Description	Default
chamber_pressure	float	Chamber pressure [Pa].	required
external_pressure	float	External pressure [Pa].	required
critical_pressure_ratio	float	Critical pressure ratio.	required

Returns:

Type	Description
bool	True if flow is choked, False otherwise.

Source code in machwave/core/compressible_flow/isentropic.py

def is_flow_choked(
    chamber_pressure: float,
    external_pressure: float,
    critical_pressure_ratio: float,
) -> bool:
    """Check if flow is choked.

    Args:
        chamber_pressure: Chamber pressure [Pa].
        external_pressure: External pressure [Pa].
        critical_pressure_ratio: Critical pressure ratio.

    Returns:
        True if flow is choked, False otherwise.
    """
    return chamber_pressure >= external_pressure / critical_pressure_ratio