spray.py

Go to the documentation of this file.

"""Spray and mixing models (J, We, SMD, x*; TMR and theta are diagnostic only)"""
 
from __future__ import annotations
 
import math
import numpy as np
from typing import Optional
from typing import Tuple, List
from engine.pipeline.config_schemas import SprayConfig
 
 
def momentum_flux_ratio(rho_O: float, u_O: float, rho_F: float, u_F: float) -> float:
    """
    Calculate momentum flux ratio J.
    
    J = (ρ_O × u_O²) / (ρ_F × u_F²)
    
    Parameters:
    -----------
    rho_O : float
        Oxidizer density [kg/m³]
    u_O : float
        Oxidizer velocity [m/s]
    rho_F : float
        Fuel density [kg/m³]
    u_F : float
        Fuel velocity [m/s]
    
    Returns:
    --------
    J : float
    """
    if u_F == 0:
        return np.inf if u_O > 0 else 0.0
    
    J = (rho_O * u_O**2) / (rho_F * u_F**2)
    return float(J)
 
 
def thrust_momentum_ratio(J: float, MR: float) -> float:
    """TMR = J/(1 + MR). DIAGNOSTIC ONLY.
 
    Not a published parameter: the pintle literature's TMR is the radial-to-axial momentum ratio
    (mdot_r v_r)/(mdot_a v_a), which this is not. Nothing in the model reads it; it is reported
    for continuity with old runs.
    """
    TMR = J / (1 + MR) if MR > 0 else J
    return float(TMR)
 
 
def spray_angle_from_J(J: float, k: float, n: float) -> float:
    """
    Calculate spray angle from momentum flux ratio J.
    
    Uses half-angle tangent form: tan(θ/2) = k × J^n
    
    PHYSICS INTUITION:
    ------------------
    Imagine the spray as a cone coming out of the injector:
    
    1. **Momentum Vector Decomposition**:
       - The fuel (RP-1) flows RADIALLY outward (perpendicular to axis)
       - The oxidizer (LOX) flows AXIALLY (along the axis, with some angle)
       - When they collide, the combined spray has both radial and axial momentum
    
    2. **Why Half-Angle?**:
       - The spray is SYMMETRIC about the axis
       - We measure angles FROM THE AXIS, not from the opposite edge
       - A particle at the edge of the spray travels at angle θ/2 from the axis
       - Its momentum vector makes angle θ/2 with the axis (not θ!)
       - When we decompose momentum: p_radial = p × sin(θ/2), p_axial = p × cos(θ/2)
       - Full angle θ = 2×(θ/2) is just a convention - we use θ/2 in calculations
    
    3. **Why Tangent?**:
       - In a conical spray, momentum vectors decompose as:
         * Radial component: p_radial = p_total × sin(θ/2)
         * Axial component: p_axial = p_total × cos(θ/2)
       - The RATIO of radial to axial momentum is:
         p_radial / p_axial = sin(θ/2) / cos(θ/2) = tan(θ/2)
       - This ratio determines how "spread out" the spray is
    
    4. **Connection to J**:
       - J = (ρ_O u_O²) / (ρ_F u_F²) measures relative momentum flux
       - High J → more oxidizer momentum → spray more axial (smaller angle)
       - Low J → more fuel momentum → spray more radial (larger angle)
       - The spray angle emerges from balancing these momentum components
    
    5. **Why Not Direct θ = k×J^n?**:
       - Direct form treats angle as arbitrary scalar
       - But spray angle comes from VECTOR balance of momenta
       - tan(θ/2) captures this vector relationship naturally
       - Also provides natural bounds: tan(θ/2) → ∞ as θ → 90°
    
    Parameters:
    -----------
    J : float
        Momentum flux ratio
    k : float
        Model coefficient
    n : float
        Model exponent
    
    Returns:
    --------
    theta : float [rad]
    """
    if J <= 0:
        return 0.0
    
    # Half-angle tangent form: tan(θ/2) = k × J^n
    # This is the physically meaningful form for spray angles
    tan_half_theta = k * (J ** n)
    theta = 2 * np.arctan(tan_half_theta)
    
    # Clamp to reasonable range [0, π/2] (0-90 degrees)
    # Note: The tan form naturally limits angles, but we clamp for safety
    return float(np.clip(theta, 0, np.pi / 2))
 
 
def spray_angle_from_TMR(TMR: float) -> float:
    """theta = arccos(1/(1 + TMR^0.75)) [rad]. DIAGNOSTIC ONLY.
 
    An invented law on the invented TMR above; no source gives it, and no model reads the angle
    (the injectors publish it as ``theta``). The vaporization march takes its axial drop speed
    from ``spray_axial_velocity``.
    """
    if TMR <= 0:
        return 0.0
    cos_theta = np.clip(1 / (1 + TMR**0.75), -1, 1)
    return float(np.arccos(cos_theta))
 
 
def weber_number(rho: float, u: float, d_char: float, sigma: float) -> float:
    """
    Calculate Weber number.
    
    We = (ρ × u² × d_char) / σ
    
    Parameters:
    -----------
    rho : float
        Density [kg/m³]
    u : float
        Velocity [m/s]
    d_char : float
        Characteristic diameter [m] (port diameter or annulus hydraulic diameter)
    sigma : float
        Surface tension [N/m]
    
    Returns:
    --------
    We : float
    """
    if sigma <= 0:
        return np.inf
    
    We = (rho * u**2 * d_char) / sigma
    return float(We)
 
 
def ohnesorge_number(mu: float, rho: float, sigma: float, d_or: float) -> float:
    """
    Calculate Ohnesorge number.
    
    Oh = μ / √(ρ × σ × d_or)
    
    Parameters:
    -----------
    mu : float
        Dynamic viscosity [Pa·s]
    rho : float
        Density [kg/m³]
    sigma : float
        Surface tension [N/m]
    d_or : float
        Orifice diameter [m]
    
    Returns:
    --------
    Oh : float
    """
    if rho <= 0 or sigma <= 0 or d_or <= 0:
        return 0.0
    
    # FIXED: Ensure sqrt input is positive
    sqrt_arg = rho * sigma * d_or
    Oh = mu / np.sqrt(max(sqrt_arg, 1e-12)) if sqrt_arg > 0 else 0.0
    return float(Oh)
 
 
def smd_lefebvre(
    d_or: float,
    We: float,
    Oh: float,
    C: float,
    m: float,
    p: float
) -> float:
    """
    Calculate Sauter Mean Diameter (SMD) using Lefebvre correlation.
    
    D32 = C × d_or × We^(-m) × (1 + Oh)^p
 
    Raises ValueError when We or d_or is not positive and finite (it used to return d_or).
    
    Parameters:
    -----------
    d_or : float
        Orifice diameter [m]
    We : float
        Weber number
    Oh : float
        Ohnesorge number
    C : float
        Model constant
    m : float
        Weber exponent
    p : float
        Ohnesorge exponent
    
    Returns:
    --------
    D32 : float [m]
    """
    if not (np.isfinite(We) and We > 0 and np.isfinite(d_or) and d_or > 0):
        raise ValueError(f"smd_lefebvre: We={We} and d_or={d_or} must be positive and finite")
 
    # Viscous damping term uses (1 + Oh)^p, NOT Oh^p. Raw Oh^p is unphysical: for a low-viscosity
    # liquid (Oh << 1) it drives D32 -> 0 as Oh -> 0, whereas vanishing viscosity must approach the
    # inviscid limit (no viscous correction => factor -> 1). (1 + Oh)^p -> 1 as Oh -> 0 and grows
    # with viscosity, which is the correct trend (more viscous => larger drops).
    D32 = C * d_or * (We ** (-m)) * ((1.0 + Oh) ** p)
    return float(D32)
 
 
#: NACA TN 4222 (Ingebo 1958) test conditions: n-heptane jet pairs in an airstream at 82 F and
#: 29.3 in Hg abs. The correlation carries no property terms, so these are its reference state.
TN4222_RHO_AIR = 29.3 * 3386.389 / (287.05 * ((82.0 - 32.0) * 5.0 / 9.0 + 273.15))   # 1.149 kg/m^3
TN4222_HEPTANE_SIGMA = 0.0197      # N/m, n-heptane at 28 C
TN4222_HEPTANE_MU = 3.8e-4         # Pa s
TN4222_HEPTANE_RHO = 684.0         # kg/m^3
#: Nukiyama-Tanasawa distribution TN 4222 fits (number ~ D^2 exp(-bD)): D30 = 3.915/b,
#: D32 = 5/b, so D32/D30 = 5/3.915.
TN4222_D32_OVER_D30 = 5.0 / 3.915
 
 
#: How TN 4222's heptane/air D32 is carried to another liquid and gas (SMDConfig.smd_property_transfer).
TN4222_PROPERTY_TRANSFERS = ("tn4087", "dombrowski_johns", "none")
 
 
def tn4222_property_factor(transfer: str, *, rho_liq: float, mu_liq: float, sigma: float,
                           rho_gas: float) -> float:
    """D32(liquid, gas) / D32(n-heptane, air at TN 4222's conditions), same jets and speeds.
 
    - "tn4087": Ingebo's crossflow exponents (NACA TN 4087: D30 ~ (We Re)^(-1/4)), so
      (sigma mu_l / rho_l)^(1/4) rho_g^(-1/4).
    - "dombrowski_johns": the inviscid sheet-breakup length scale of Dombrowski & Johns (Chem. Eng.
      Sci. 18, 1963), D ~ (sigma^2 / (rho_g rho_l U^4))^(1/6) at fixed sheet geometry and speed, so
      (sigma^2 / (rho_g rho_l))^(1/6) relative to heptane/air.
    - "none": 1 (heptane in air as measured).
    The spread between them is the uncertainty of carrying a single-liquid correlation to
    LOX/ethanol at chamber density; a measured D32 (smd_scale) replaces it.
    """
    if transfer not in TN4222_PROPERTY_TRANSFERS:
        raise ValueError(f"TN 4222 property transfer {transfer!r} is not one of {TN4222_PROPERTY_TRANSFERS}")
    if transfer == "none":
        return 1.0
    for name, v in (("rho_liq", rho_liq), ("mu_liq", mu_liq), ("sigma", sigma), ("rho_gas", rho_gas)):
        if not (np.isfinite(v) and v > 0):
            raise ValueError(f"tn4222_property_factor: {name}={v} must be positive and finite")
    if transfer == "tn4087":
        prop = (sigma * mu_liq / rho_liq) / (TN4222_HEPTANE_SIGMA * TN4222_HEPTANE_MU / TN4222_HEPTANE_RHO)
        return float((rho_gas / TN4222_RHO_AIR) ** -0.25 * prop ** 0.25)
    ref = TN4222_HEPTANE_SIGMA ** 2 / (TN4222_RHO_AIR * TN4222_HEPTANE_RHO)
    return float((sigma ** 2 / (rho_gas * rho_liq) / ref) ** (1.0 / 6.0))
 
 
def tn4222_transfer_model(smd_cfg) -> str:
    """The property transfer a spray.smd config asks for: 'none' when smd_property_scaling is off,
    else smd_property_transfer (default 'tn4087')."""
    if not bool(getattr(smd_cfg, "smd_property_scaling", True)):
        return "none"
    return str(getattr(smd_cfg, "smd_property_transfer", "tn4087") or "tn4087")
 
 
def smd_impinging_tn4222(
    d_jet: float,
    v_jet: float,
    rho_liq: float,
    mu_liq: float,
    sigma: float,
    rho_gas: float,
    *,
    dv: Optional[float] = None,
    scale_properties: bool = True,
    transfer: Optional[str] = None,
) -> float:
    """Sauter mean diameter of an impinging-jet pair, Ingebo, NACA TN 4222 (1958).
 
        Dj / D30 = 2.64 (Dj Vj)^(1/2) + 0.97 Dj dV          (Dj, D30 in inches; Vj, dV in ft/s)
        D32 = (5 / 3.915) D30                                 (Nukiyama-Tanasawa, as fitted)
 
    ``v_jet`` is the stream's own jet velocity and ``dv`` the jet-to-gas velocity difference; a
    face-injected jet enters gas at rest, dV = Vj (TN 4222 eq. 11), the default.
 
    TN 4222 measured n-heptane only. ``transfer`` (see ``tn4222_property_factor``) carries it to
    other liquids and gas densities; None means "tn4087" when ``scale_properties`` else "none".
    That transfer is an assumption: measure D32 and set ``spray.smd.smd_scale``.
 
    Raises ValueError on a non-positive or non-finite input (it used to return d_jet, a drop the
    size of the orifice, without saying so).
 
    Replaces an earlier function that cited this report but ran TN 4087's crossflow form with the
    We*Re product on the jets' closing speed, and read its D30 as D32: about half the drop size
    TN 4222 measured.
    """
    for name, v in (("d_jet", d_jet), ("v_jet", v_jet), ("rho_liq", rho_liq), ("mu_liq", mu_liq),
                    ("sigma", sigma), ("rho_gas", rho_gas)):
        if not (np.isfinite(v) and v > 0):
            raise ValueError(f"smd_impinging_tn4222: {name}={v} must be positive and finite")
    if transfer is None:
        transfer = "tn4087" if scale_properties else "none"
    dj_in = d_jet / 0.0254
    vj = v_jet / 0.3048
    dV = (v_jet if dv is None else abs(float(dv))) / 0.3048
    d30 = dj_in / (2.64 * math.sqrt(dj_in * vj) + 0.97 * dj_in * dV) * 0.0254
    d32 = TN4222_D32_OVER_D30 * d30
    return float(d32 * tn4222_property_factor(transfer, rho_liq=rho_liq, mu_liq=mu_liq,
                                              sigma=sigma, rho_gas=rho_gas))
 
 
def evaporation_constant_m2_s(
    *,
    Tc: float,
    Pc: float,
    rho_g: float,
    rho_l: float,
    L_vap: float,
    T_boil: float,
    cp_g: float = 2200.0,
    D0: float = 2.0e-5,
    C_evap: float = 1.562,
) -> float:
    """d²-law evaporation constant k_evap [m²/s], derived from propellant properties.
 
    Replaces a single hardcoded ``K`` that was identical in all 14 configs and for every
    propellant. That constant's VALUE was defensible for ethanol at ~3100 K (it implies
    k_evap = 3.33e-6 m²/s, within 10% of a T^2 extrapolation of measured ethanol data), but its
    SHAPE was wrong: it could not distinguish ethanol from methane from RP-1, and it did not
    move with chamber temperature or pressure at all -- so a 200 psi design and a 600 psi design
    evaporated at identical rates.
 
        k_evap = C * (8 * rho_g * D_v / rho_l) * ln(1 + B_M)
        B_M    = cp_g * (Tc - T_boil) / L_vap          (Spalding mass transfer number)
        D_v    = D0 * (Tc/300)^1.75 * (101325/Pc)      (vapor diffusivity)
 
    Every input except C and D0 already exists in the propellant presets, so this needs no new
    tables. ``C_evap`` is the one calibration constant, anchored so ethanol at 3094 K / 450 psi
    reproduces the 3.33e-6 m²/s the old fixed K implied -- i.e. this is a re-shaping of a value
    that was already right for the propellant it was tuned on, not a change of magnitude.
 
    The same form is already written (and disabled) in combustion_physics.py under
    ``SPALDING_DIAGNOSTIC_ENABLED``, where it also hardcodes ``fuel="RP-1"``.
    """
    if not all(np.isfinite(v) for v in (Tc, Pc, rho_g, rho_l, L_vap, T_boil)):
        return float("nan")
    if Tc <= 0 or Pc <= 0 or rho_g <= 0 or rho_l <= 0 or L_vap <= 0:
        return float("nan")
    D_v = float(D0) * (float(Tc) / 300.0) ** 1.75 * (101325.0 / float(Pc))
    B_M = float(cp_g) * max(0.0, float(Tc) - float(T_boil)) / float(L_vap)
    k = float(C_evap) * (8.0 * float(rho_g) * D_v / float(rho_l)) * float(np.log1p(B_M))
    return float(k) if (np.isfinite(k) and k > 0) else float("nan")
 
 
def tau_evap_from_k(D32: float, k_evap_m2_s: float) -> float:
    """Standard d²-law evaporation time: τ = D32² / k_evap. Note k is in m²/s, NOT its reciprocal."""
    if not (np.isfinite(D32) and np.isfinite(k_evap_m2_s)) or k_evap_m2_s <= 0 or D32 <= 0:
        return float("nan")
    return float((D32 ** 2) / k_evap_m2_s)
 
 
def spray_axial_velocity(
    *,
    mdot_O: float, u_O: float, theta_O_deg: float,
    mdot_F: float, u_F: float, theta_F_deg: float,
) -> float:
    """Axial velocity of the combined spray after the doublet collides [m/s].
 
    Momentum balance on the pair: the two jets arrive canted at theta_O / theta_F from the
    chamber axis, and what survives downstream is the momentum-weighted axial component
 
        u_axial = (mdot_O*u_O*cos(theta_O) + mdot_F*u_F*cos(theta_F)) / (mdot_O + mdot_F)
 
    This is the velocity that TRANSPORTS droplets down the chamber, and therefore the one that
    belongs in x* = u * tau_evap. The code previously used ``u_rel``, the jet-to-jet RELATIVE
    velocity -- correct for atomisation (it is what shears the sheet, and it belongs in the
    Ingebo Weber number) but wrong for transport: droplets do not travel downstream at the
    closing speed of two jets. Measured on a converged design, u_rel was 28.4 m/s against a true
    axial 13.8 m/s, so x* -- and hence the required chamber length -- was over-predicted 2.06x.
    """
    mt = float(mdot_O) + float(mdot_F)
    if not np.isfinite(mt) or mt <= 0:
        return float("nan")
    ax = (float(mdot_O) * float(u_O) * np.cos(np.deg2rad(float(theta_O_deg)))
          + float(mdot_F) * float(u_F) * np.cos(np.deg2rad(float(theta_F_deg))))
    u = ax / mt
    return float(u) if (np.isfinite(u) and u > 0) else float("nan")
 
 
def tau_evap(D32: float, K: float) -> float:
    """
    Calculate evaporation time.
 
    τ_evap = K × D32²
    
    Parameters:
    -----------
    D32 : float
        Sauter Mean Diameter [m]
    K : float
        Evaporation constant [s/m²]
    
    Returns:
    --------
    tau : float [s]
    """
    tau = K * (D32 ** 2)
    return float(tau)
 
 
def xstar(U_rel: float, tau_evap: float) -> float:
    """
    Calculate evaporation length x*.
    
    x* = U_rel × τ_evap
    
    Parameters:
    -----------
    U_rel : float
        Relative velocity [m/s]
    tau_evap : float
        Evaporation time [s]
    
    Returns:
    --------
    x_star : float [m]
    """
    x_star = U_rel * tau_evap
    return float(x_star)
 
 
def check_spray_constraints(
    We_O: float,
    We_F: float,
    x_star: float,
    config: SprayConfig
) -> Tuple[bool, List[str]]:
    """
    Check if spray constraints are satisfied.
    
    Constraints:
    - We_i ≥ We_min
    - x* < x_limit
    
    Parameters:
    -----------
    We_O : float
        Oxidizer Weber number
    We_F : float
        Fuel Weber number
    x_star : float
        Evaporation length [m]
    config : SprayConfig
        Spray configuration
    
    Returns:
    --------
    constraints_satisfied : bool
    violations : list[str]
        List of violated constraint names
    """
    violations = []
    
    # Check Weber number constraints
    We_min = config.weber.get("We_min", 15.0)
    if We_O < We_min:
        violations.append(f"We_O < We_min ({We_O:.2f} < {We_min:.2f})")
    if We_F < We_min:
        violations.append(f"We_F < We_min ({We_F:.2f} < {We_min:.2f})")
    
    # Check evaporation length constraint
    if config.evaporation.use_constraint:
        if x_star >= config.evaporation.x_star_limit:
            violations.append(
                f"x* >= x_limit ({x_star:.4f} >= {config.evaporation.x_star_limit:.4f} m)"
            )
    
 
    constraints_satisfied = len(violations) == 0
    return constraints_satisfied, violations
 
 
def smd_pintle(
    L_open: float,
    V_rel: float,
    rho_f: float,
    mu_f: float,
    sigma_f: float,
    C: float,
    B: float,
    n: float,
    p: float,
) -> float:
    """
    Calculate Sauter Mean Diameter (SMD) for pintle injector using relative velocity physics.
    
    Formula: SMD = C * L_open * We_rel^(-n) * (1 + B * Oh_f)^p
    
    Parameters:
    -----------
    L_open : float
        Characteristic length (Pintle Opening / Gap Height) [m]
    V_rel : float
        Relative velocity magnitude between streams [m/s]
    rho_f : float
        Fuel (sheet) density [kg/m³]
    mu_f : float
        Fuel (sheet) dynamic viscosity [Pa·s]
    sigma_f : float
        Fuel surface tension [N/m]
    C, B, n, p : float
        Correlation constants
        
    Returns:
    --------
    D32 : float [m]
    """
    if L_open <= 0:
        raise ValueError(f"smd_pintle: L_open must be positive, got {L_open}")
    if rho_f <= 0:
        raise ValueError(f"smd_pintle: rho_f must be positive, got {rho_f}")
    if sigma_f <= 0:
        raise ValueError(f"smd_pintle: sigma_f must be positive, got {sigma_f}")
    if mu_f < 0:
        raise ValueError(f"smd_pintle: mu_f must be non-negative, got {mu_f}")
    
    if C <= 0:
        raise ValueError(f"smd_pintle: C must be positive, got {C}")
    if n <= 0:
        raise ValueError(f"smd_pintle: n must be positive, got {n}")
    if B < 0:
        raise ValueError(f"smd_pintle: B must be non-negative, got {B}")
    if p < 0:
        raise ValueError(f"smd_pintle: p must be non-negative, got {p}")
 
    if V_rel < 0:
        raise ValueError(f"smd_pintle: V_rel must be non-negative, got {V_rel}")
    
    if V_rel == 0:
        import warnings
        warnings.warn("smd_pintle: V_rel is 0.0, We=0. SMD calculation relying on Ohnesorge term only.")
 
    # Calculate Weber number based on relative velocity
    # We = rho * V^2 * L / sigma
    We_rel = (rho_f * (V_rel ** 2) * L_open) / sigma_f
    
    # Calculate Ohnesorge number for the liquid sheet
    # Oh = mu / sqrt(rho * sigma * L)
    denom = np.sqrt(rho_f * sigma_f * L_open)
    Oh_f = mu_f / denom if denom > 0 else 0.0
    
    # Apply correlation
    # SMD = C * L_open * We^(-n) * (1 + B * Oh)^p
    factor_we = (We_rel ** (-n)) if We_rel > 0 else 1.0
    factor_oh = (1.0 + B * Oh_f) ** p
    
    D32 = C * L_open * factor_we * factor_oh
        
    return float(D32)