Skip to content

core.compressible_flow

Compressible flow relations for convergent-divergent nozzle analysis.

  • Isentropic flow — Critical pressure ratio, exit Mach from expansion ratio, exit pressure, choked-flow detection.
  • Nozzle performance — Ideal and corrected thrust coefficient (Cf), optimal expansion ratio for a given altitude, effective exit conditions under overexpanded flow separation, thrust from Cf.

machwave.core.compressible_flow

Compressible flow theory and analysis.

FlowBranch

Bases: StrEnum

Branch of the area-Mach relation that a solution is taken from.

Source code in machwave/core/compressible_flow/isentropic.py 
 9
10
11
12
13
class FlowBranch(enum.StrEnum):
    """Branch of the area-Mach relation that a solution is taken from."""

    SUBSONIC = "subsonic"
    SUPERSONIC = "supersonic"

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 
22
23
24
25
26
27
28
29
30
31
32
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, branch=FlowBranch.SUPERSONIC)

Get exit Mach number from expansion ratio.

Every expansion ratio above unity satisfies the area-Mach relation at one subsonic and one supersonic Mach number. The supersonic root is returned by default, which is the started nozzle flowing full. A nozzle that is not started carries a normal shock in its divergent section and leaves subsonically, so its physical root is the subsonic one and has to be asked for explicitly.

Parameters:

Name Type Description Default
k float

Isentropic exponent.

 required
expansion_ratio float

Expansion ratio (A_exit / A_throat).

 required
branch FlowBranch

Branch of the area-Mach relation to solve on.

 SUPERSONIC

Returns:

Type Description
float

Exit Mach number.

Raises:

Type Description
ValueError

If the expansion ratio lies outside the range the branch covers.
Source code in machwave/core/compressible_flow/isentropic.py 
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def get_exit_mach_from_expansion_ratio(
    k: float,
    expansion_ratio: float,
    branch: FlowBranch = FlowBranch.SUPERSONIC,
) -> float:
    """
    Get exit Mach number from expansion ratio.

    Every expansion ratio above unity satisfies the area-Mach relation at one
    subsonic and one supersonic Mach number. The supersonic root is returned by
    default, which is the started nozzle flowing full. A nozzle that is not started
    carries a normal shock in its divergent section and leaves subsonically, so its
    physical root is the subsonic one and has to be asked for explicitly.

    Args:
        k: Isentropic exponent.
        expansion_ratio: Expansion ratio (A_exit / A_throat).
        branch: Branch of the area-Mach relation to solve on.

    Returns:
        Exit Mach number.

    Raises:
        ValueError: If the expansion ratio lies outside the range the branch covers.
    """
    mach_limit = MACH_LIMITS[branch]
    maximum_expansion_ratio = get_maximum_expansion_ratio(k, branch)
    if not 1.0 <= expansion_ratio <= maximum_expansion_ratio:
        raise ValueError(
            f"Expansion ratio {expansion_ratio} has no {branch} solution for "
            f"isentropic exponent {k}. The {branch} branch covers expansion ratios "
            f"from 1 at the sonic throat to {maximum_expansion_ratio:.6g} at the "
            f"Mach {mach_limit:g} limit of the solver."
        )

    if expansion_ratio <= get_expansion_ratio_from_exit_mach(SONIC_MACH, k):
        return SONIC_MACH  # Both roots meet at the throat.

    return cast(
        float,
        scipy.optimize.brentq(
            lambda m: get_expansion_ratio_from_exit_mach(m, k) - expansion_ratio,
            a=min(SONIC_MACH, mach_limit),
            b=max(SONIC_MACH, mach_limit),
        ),
    )

get_exit_pressure(k_exhaust, expansion_ratio, chamber_pressure, branch=FlowBranch.SUPERSONIC)

Get exit pressure from isentropic relations.

Parameters:

Name Type Description Default
k_exhaust float

Isentropic exponent at exit.

 required
expansion_ratio float

Expansion ratio.

 required
chamber_pressure float

Chamber pressure [Pa].

 required
branch FlowBranch

Branch of the area-Mach relation to solve the exit Mach number on.

 SUPERSONIC

Returns:

Type Description
float

Exit pressure [Pa].
Source code in machwave/core/compressible_flow/isentropic.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def get_exit_pressure(
    k_exhaust: float,
    expansion_ratio: float,
    chamber_pressure: float,
    branch: FlowBranch = FlowBranch.SUPERSONIC,
) -> float:
    """
    Get exit pressure from isentropic relations.

    Args:
        k_exhaust: Isentropic exponent at exit.
        expansion_ratio: Expansion ratio.
        chamber_pressure: Chamber pressure [Pa].
        branch: Branch of the area-Mach relation to solve the exit Mach number on.

    Returns:
        Exit pressure [Pa].
    """
    exit_mach = get_exit_mach_from_expansion_ratio(k_exhaust, expansion_ratio, branch)
    return chamber_pressure * (1 + 0.5 * (k_exhaust - 1) * exit_mach**2) ** (
        -k_exhaust / (k_exhaust - 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 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
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_terms(chamber_pressure, exit_pressure, external_pressure, expansion_ratio, k_exhaust)

Get the momentum and pressure terms of the ideal thrust coefficient.

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_exhaust float

Isentropic exponent at exit.

 required

Returns:

Type Description
float

Momentum term then pressure term. The pressure term is positive when under
float

expanded and negative when over expanded.
References

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

Source code in machwave/core/compressible_flow/nozzle.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
def get_ideal_thrust_coefficient_terms(
    chamber_pressure: float,
    exit_pressure: float,
    external_pressure: float,
    expansion_ratio: float,
    k_exhaust: float,
) -> tuple[float, float]:
    """
    Get the momentum and pressure terms of the ideal thrust coefficient.

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

    Returns:
        Momentum term then pressure term. The pressure term is positive when under
        expanded and negative when over expanded.

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

get_maximum_expansion_ratio(k, branch)

Get the largest expansion ratio the exit Mach solver resolves on a branch.

Parameters:

Name Type Description Default
k float

Isentropic exponent.

 required
branch FlowBranch

Branch of the area-Mach relation.

 required

Returns:

Type Description
float

Expansion ratio at the Mach limit of the branch.
Source code in machwave/core/compressible_flow/isentropic.py 
51
52
53
54
55
56
57
58
59
60
61
62
def get_maximum_expansion_ratio(k: float, branch: FlowBranch) -> float:
    """
    Get the largest expansion ratio the exit Mach solver resolves on a branch.

    Args:
        k: Isentropic exponent.
        branch: Branch of the area-Mach relation.

    Returns:
        Expansion ratio at the Mach limit of the branch.
    """
    return get_expansion_ratio_from_exit_mach(MACH_LIMITS[branch], k)

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 
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
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_separated_exit_conditions(k_exhaust, expansion_ratio, chamber_pressure, external_pressure, separation_pressure_ratio)

Get the effective exit conditions accounting for flow separation.

The model implemented is based on the work of Summerfield et al. (1954) and assumes that flow separation occurs when the exit pressure is below a certain fraction of the ambient pressure.

Parameters:

Name Type Description Default
k_exhaust float

Isentropic exponent at exit.

 required
expansion_ratio float

Geometric expansion ratio.

 required
chamber_pressure float

Chamber pressure [Pa].

 required
external_pressure float

Ambient pressure [Pa].

 required
separation_pressure_ratio float

Separation-to-ambient pressure ratio.

 required

Returns:

Type Description
float

Effective expansion ratio and effective exit pressure [Pa]. The effective
float

expansion ratio never falls below the throat value of unity.

Raises:

Type Description
ValueError

If the geometric expansion ratio lies outside the range the supersonic branch of the area-Mach relation covers.
References

Summerfield, M., Foster, C. R., & Swan, W. C. (1954). Flow separation in overexpanded supersonic exhaust nozzles.

Source code in machwave/core/compressible_flow/nozzle.py 
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def get_separated_exit_conditions(
    k_exhaust: float,
    expansion_ratio: float,
    chamber_pressure: float,
    external_pressure: float,
    separation_pressure_ratio: float,
) -> tuple[float, float]:
    """
    Get the effective exit conditions accounting for flow separation.

    The model implemented is based on the work of Summerfield et al. (1954) and assumes
    that flow separation occurs when the exit pressure is below a certain fraction of
    the ambient pressure.

    Args:
        k_exhaust: Isentropic exponent at exit.
        expansion_ratio: Geometric expansion ratio.
        chamber_pressure: Chamber pressure [Pa].
        external_pressure: Ambient pressure [Pa].
        separation_pressure_ratio: Separation-to-ambient pressure ratio.

    Returns:
        Effective expansion ratio and effective exit pressure [Pa]. The effective
        expansion ratio never falls below the throat value of unity.

    Raises:
        ValueError: If the geometric expansion ratio lies outside the range the
            supersonic branch of the area-Mach relation covers.

    References:
        Summerfield, M., Foster, C. R., & Swan, W. C. (1954). Flow separation in
        overexpanded supersonic exhaust nozzles.
    """
    exit_pressure = isentropic.get_exit_pressure(
        k_exhaust, expansion_ratio, chamber_pressure
    )
    separation_pressure = separation_pressure_ratio * external_pressure
    if exit_pressure >= separation_pressure:
        return expansion_ratio, exit_pressure

    sonic_pressure = chamber_pressure * isentropic.get_critical_pressure_ratio(
        k_exhaust
    )
    if separation_pressure >= sonic_pressure:
        return 1.0, sonic_pressure

    # The exit pressure falls monotonically from the sonic throat to the geometric
    # exit, so bracketing from the throat always contains the separation point.
    effective_expansion_ratio = cast(
        float,
        scipy.optimize.brentq(
            lambda ratio: (
                isentropic.get_exit_pressure(k_exhaust, ratio, chamber_pressure)
                - separation_pressure
            ),
            a=1.0,
            b=expansion_ratio,
        ),
    )
    return effective_expansion_ratio, separation_pressure

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 
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
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

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 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
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