core.equations

machwave.core.equations

This module groups pure functions that return the right-hand side of differential equations (DEs) used throughout Machwave.

They are intentionally stateless and side-effect-free so that any numerical integrator (RK4, RK45, implicit Euler, JAX-based solvers, etc.) can consume them without modification.

Variable names may not be according to PEP8, since the common notation in academic literature is preferred for readability.

compute_chamber_pressure_mass_balance_lre(P0, R, T0, V0, At, k, m_dot_ox, m_dot_fuel)

Calculates the chamber pressure by solving the mass balance equation for liquid rocket engines.

Parameters:

Name	Type	Description	Default
P0	float	Chamber pressure [Pa].	required
R	float	Gas constant per molecular weight [J/(kg·K)].	required
T0	float	Flame temperature [K].	required
V0	float	Chamber free volume [m^3].	required
At	float	Nozzle throat area [m^2].	required
k	float	Isentropic exponent of the mix.	required
m_dot_ox	float	Instantaneous oxidizer mass flow rate [kg/s].	required
m_dot_fuel	float	Instantaneous fuel mass flow rate [kg/s].	required

Returns:

Type	Description
tuple[float]	Derivative of chamber pressure with respect to time.

Source code in machwave/core/equations/lre_mass_balance.py

def compute_chamber_pressure_mass_balance_lre(
    P0: float,
    R: float,
    T0: float,
    V0: float,
    At: float,
    k: float,
    m_dot_ox: float,
    m_dot_fuel: float,
) -> tuple[float]:
    """
    Calculates the chamber pressure by solving the mass balance equation for liquid
    rocket engines.

    Args:
        P0: Chamber pressure [Pa].
        R: Gas constant per molecular weight [J/(kg·K)].
        T0: Flame temperature [K].
        V0: Chamber free volume [m^3].
        At: Nozzle throat area [m^2].
        k: Isentropic exponent of the mix.
        m_dot_ox: Instantaneous oxidizer mass flow rate [kg/s].
        m_dot_fuel: Instantaneous fuel mass flow rate [kg/s].

    Returns:
        Derivative of chamber pressure with respect to time.
    """
    m_dot_out = (
        At
        * P0
        * k
        * (np.sqrt((2 / (k + 1)) ** ((k + 1) / (k - 1))))
        / (np.sqrt(k * R * T0))
    )
    m_dot_in = m_dot_ox + m_dot_fuel
    dP0_dt = (R * T0 / V0) * (m_dot_in - m_dot_out)
    return (dP0_dt,)

compute_chamber_pressure_mass_balance_srm(P0, Pe, Ab, V0, At, pp, k, R, T0, r, Cd=1.0)

Calculates the chamber pressure by solving Hans Seidel's differential equation.

This differential equation was presented in Seidel's paper named "Transient Chamber Pressure and Thrust in Solid Rocket Motors", published in March, 1965.

Parameters:

Name	Type	Description	Default
P0	float	Chamber pressure [Pa].	required
Pe	float	External pressure [Pa].	required
Ab	float	Burn area [m^2].	required
V0	float	Chamber free volume [m^3].	required
At	float	Nozzle throat area [m^2].	required
pp	float	Propellant density [kg/m^3].	required
k	float	Isentropic exponent of the mix.	required
R	float	Gas constant per molecular weight [J/(kg·K)].	required
T0	float	Flame temperature [K].	required
r	float	Propellant burn rate [m/s].	required
Cd	float	Discharge coefficient, default is 1.0.	1.0

Returns:

Type	Description
tuple[float]	Derivative of chamber pressure with respect to time.

Source code in machwave/core/equations/srm_mass_balance.py

def compute_chamber_pressure_mass_balance_srm(
    P0: float,
    Pe: float,
    Ab: float,
    V0: float,
    At: float,
    pp: float,
    k: float,
    R: float,
    T0: float,
    r: float,
    Cd: float = 1.0,
) -> tuple[float]:
    """
    Calculates the chamber pressure by solving Hans Seidel's differential equation.

    This differential equation was presented in Seidel's paper named "Transient Chamber
    Pressure and Thrust in Solid Rocket Motors", published in March, 1965.

    Args:
        P0: Chamber pressure [Pa].
        Pe: External pressure [Pa].
        Ab: Burn area [m^2].
        V0: Chamber free volume [m^3].
        At: Nozzle throat area [m^2].
        pp: Propellant density [kg/m^3].
        k: Isentropic exponent of the mix.
        R: Gas constant per molecular weight [J/(kg·K)].
        T0: Flame temperature [K].
        r: Propellant burn rate [m/s].
        Cd: Discharge coefficient, default is 1.0.

    Returns:
        Derivative of chamber pressure with respect to time.

    """
    critical_pressure_ratio = get_critical_pressure_ratio(k=k)
    Pr = Pe / P0

    if Pr <= critical_pressure_ratio:  # choked
        H = (k**0.5) * (2 / (k + 1)) ** ((k + 1) / (2 * (k - 1)))
    else:  # sub-critical
        H = ((k / (k - 1)) ** 0.5) * Pr ** (1 / k) * (1 - Pr ** ((k - 1) / k)) ** 0.5

    m_dot_gen = pp * r * Ab
    m_dot_exit = Cd * P0 * At * H / (R * T0) ** 0.5

    dP0_dt = (R * T0 / V0) * (m_dot_gen - m_dot_exit)
    return (dP0_dt,)

compute_point_mass_trajectory(y, v, T, D, M, g)

Returns the derivatives of elevation and velocity.

Parameters:

Name	Type	Description	Default
y	float	Instant elevation [m].	required
v	float	Instant velocity [m/s].	required
T	float	Instant thrust [N].	required
D	float	Instant drag constant (Cd * A * rho / 2) [kg/m].	required
M	float	Instant total mass [kg].	required
g	float	Instant acceleration of gravity [m/s^2].	required

Returns:

Type	Description
tuple[float, float]	Derivatives of elevation and velocity.

Source code in machwave/core/equations/point_mass_trajectory.py

def compute_point_mass_trajectory(
    y: float, v: float, T: float, D: float, M: float, g: float
) -> tuple[float, float]:
    """
    Returns the derivatives of elevation and velocity.

    Args:
        y: Instant elevation [m].
        v: Instant velocity [m/s].
        T: Instant thrust [N].
        D: Instant drag constant (Cd * A * rho / 2) [kg/m].
        M: Instant total mass [kg].
        g: Instant acceleration of gravity [m/s^2].

    Returns:
        Derivatives of elevation and velocity.
    """
    if v < 0:
        x = -1
    else:
        x = 1

    dv_dt = (T - x * D * (v**2)) / M - g
    dy_dt = v

    return (dy_dt, dv_dt)