impinging.py

Go to the documentation of this file.

"""Impinging injector model implementation."""
 
from __future__ import annotations
 
import logging
import math
import os
from typing import Tuple, Dict, Any
 
import numpy as np
 
_LOG = logging.getLogger(__name__)
_print_raw = os.environ.get("ENGINE_PRINT_IMPINGING_FEED_CLOSURE", "")
_PRINT_FEED_ORIFICE_CLOSURE = _print_raw != "" and str(_print_raw).lower() not in ("0", "false", "no")
 
_STREAM_ROOT_MAX_EVAL = 200
 
from engine.pipeline.config_schemas import PintleEngineConfig, ImpingingInjectorConfig
from engine.pipeline.feed_loss import delta_p_feed
from engine.pipeline.thermal.regen_cooling import delta_p_regen_channels
from engine.core.discharge import (
    cd_from_re,
    cd_inf_from_inlet_geometry,
    cd_inf_from_orifice_diameter,
    calculate_reynolds_number,
)
from engine.pipeline.assumptions import assume
from engine.core.spray import (
    evaporation_constant_m2_s,
    tau_evap_from_k,
    spray_axial_velocity,
    momentum_flux_ratio,
    thrust_momentum_ratio,
    spray_angle_from_J,
    spray_angle_from_TMR,
    weber_number,
    ohnesorge_number,
    smd_impinging_tn4222,
    tau_evap,
    xstar,
    check_spray_constraints,
)
 
from . import InjectorModel
from engine.core.injectors.layout import effective_discharge
from engine.core.injectors.flow_capacity import (
    effective_flow_areas_from_cd,
    merge_effective_area_warnings,
)
 
 
 
def _stream_flow(Pc, P_tank, rho, A, d_hyd, mu, discharge, T_in, cd_cap, feed_dp,
                 *, cav=None, network=None, detail=None):
    """One stream's closure at a fixed Pc: ṁ = Cd(Re(ṁ), P_inj) A √(2ρ (P_tank − Δp_feed(ṁ) − Pc)).
 
    g(ṁ) = ṁ − Cd A √(2ρ Δp_inj) is negative at ṁ = 0 and non-negative at
    ṁ_max = Cd_cap A √(2ρ (P_tank − Pc)) (Δp_feed ≥ 0, Cd ≤ Cd_cap), so the root is bracketed and
    found by Illinois false position. A tank at or below Pc does not flow.
    Mirrored by ``engine.accel.kernels._stream_flow``.
 
    Returns (ṁ, Cd, P_inj, Δp_feed, Δp_inj, evaluations).
    """
 
    def cd_of(Re, P_in):
        Cd = float(min(cd_from_re(Re, discharge, P_inlet=P_in, T_inlet=T_in, d_hyd_m=d_hyd), cd_cap))
        if cav is not None and P_in > Pc:
            # Cavitating orifice (Nurick 1976): Cd = Cc sqrt(K), K = (P_in - P_v)/(P_in - Pc).
            P_v, Cc = cav
            K = (P_in - P_v) / (P_in - Pc)
            Cd = float(min(Cd, Cc * math.sqrt(max(K, 0.0))))
        return Cd
 
    def state(m):
        dpf = float(feed_dp(m))
        Pi = float(P_tank - dpf)
        dpi = Pi - Pc
        if network is not None:
            # Cd stays the holes' own (flow-weighted); the manifold's loss shows as the gap
            # between the port head and the flow, reported as Cd_eff_manifold.
            m_b, holes, Cd = network.march(m, Pi, Pc, cd_of)
            if detail is not None:
                detail["holes"] = holes
                detail["cd_eff"] = (float(m_b / (A * math.sqrt(2.0 * rho * dpi)))
                                    if dpi > 0 and m_b > 0 else float("nan"))
        else:
            u = m / (rho * A) if A > 0 else 0.0
            Re = calculate_reynolds_number(rho, u, d_hyd, mu)
            Cd = cd_of(Re, Pi)
            m_b = float(Cd * A * np.sqrt(2.0 * rho * dpi)) if dpi > 0 else 0.0
        return m - m_b, Cd, Pi, dpf, max(0.0, dpi)
 
    if not A > 0 or not P_tank > Pc:
        _, Cd, Pi, dpf, dpi = state(0.0)
        return 0.0, Cd, Pi, dpf, dpi, 1
    hi = float(cd_cap * A * np.sqrt(2.0 * rho * (P_tank - Pc)))
    a, b = 0.0, hi
    ga, gb = state(a)[0], state(b)[0]
    n = 2
    tol = 1e-13 * hi
    if gb <= 0.0:
        m = b
    elif ga >= 0.0:  # Cd = 0 at rest: nothing drives the flow
        m = a
    else:
        side = 0
        m = b
        for _ in range(_STREAM_ROOT_MAX_EVAL):
            m = (a * gb - b * ga) / (gb - ga)
            gm = state(m)[0]
            n += 1
            if abs(gm) <= tol or b - a <= tol:
                break
            if gm > 0.0:
                b, gb = m, gm
                if side == -1:
                    ga *= 0.5
                side = -1
            else:
                a, ga = m, gm
                if side == 1:
                    gb *= 0.5
                side = 1
        else:
            _LOG.warning("impinging injector: stream flow root not within tolerance after %d evaluations", n)
    _, Cd, Pi, dpf, dpi = state(m)
    return float(m), Cd, Pi, dpf, dpi, n + 1
 
 
 
def _churchill_f(Re: float, rel_rough: float) -> float:
    """Darcy friction factor, Churchill (1977), laminar through fully rough in one expression."""
    if not Re > 0:
        return 0.0
    A = (2.457 * math.log(1.0 / ((7.0 / Re) ** 0.9 + 0.27 * rel_rough))) ** 16
    B = (37530.0 / Re) ** 16
    return 8.0 * ((8.0 / Re) ** 12 + 1.0 / (A + B) ** 1.5) ** (1.0 / 12.0)
 
 
class _RingManifold:
    """One stream's back channel as a dividing-flow ring, marched hole by hole.
 
    Each of ``n_ports`` ports splits two ways, so ``2 n_ports`` identical branches each feed
    ``n_holes / (2 n_ports)`` holes, the first half a pitch from the port (port between holes).
    Along a branch:
 
    - into the branch from the port (after the feed line's velocity head is dumped there):
      p = P_port - (1 + K_ent) rho u^2 / 2;
    - between holes, Darcy friction (Churchill) on the channel's hydraulic diameter;
    - each hole flows from the LOCAL STATIC pressure, Cd A sqrt(2 rho (p - Pc)) -- with the
      approach flow normal to the hole no velocity head is recovered (Rohde, Richards & Metger,
      NASA TN D-5467, 1969); Cd from the hole's own Reynolds number (and cavitation);
    - past each hole the branch slows and regains C_R rho (u1^2 - u2^2)/2 of static pressure
      (dividing-flow manifold; Acrivos, Babcock & Pigford 1959; Bajura 1971).
 
    The outer root on the stream's total flow closes the branch flow to zero at the dead end.
    Not modelled: the extra lip separation of crossflow over the hole inlet (TN D-5467 Fig. 7;
    Strakey & Talley), which lowers Cd further near the port -- measure it by cold flow.
    """
 
    def __init__(self, *, n_holes, n_ports, A_ch, D_h, r_ring, A_hole, d_hole, rho, mu,
                 K_ent, C_R, roughness):
        self.n_ports = int(max(1, n_ports))
        per = n_holes / (2.0 * self.n_ports)
        self.h = int(round(per))
        if abs(per - self.h) > 1e-9 or self.h < 1:
            from engine.pipeline.assumptions import assume
            assume("injector.plate.channel_inlets", self.n_ports,
                   reason=f"{n_holes} holes do not split evenly over {2 * self.n_ports} branches; "
                          f"{self.h} per branch used")
            self.h = max(1, self.h)
        self.A_ch, self.D_h, self.A_hole, self.d_hole = A_ch, D_h, A_hole, d_hole
        self.rho, self.mu, self.K_ent, self.C_R = rho, mu, K_ent, C_R
        self.s = 2.0 * math.pi * r_ring / float(n_holes)
        self.eps = roughness / D_h if D_h > 0 else 0.0
        self.scale = float(n_holes) / (2.0 * self.n_ports * self.h)
 
    def march(self, m_total, P_port, Pc, cd_of):
        rho, A_ch = self.rho, self.A_ch
        mb = m_total / (2.0 * self.n_ports)
        u = mb / (rho * A_ch)
        p = P_port - (1.0 + self.K_ent) * 0.5 * rho * u * u
        holes, cds = [], []
        for j in range(self.h):
            seg = 0.5 * self.s if j == 0 else self.s
            Re_ch = rho * abs(u) * self.D_h / self.mu
            p -= _churchill_f(Re_ch, self.eps) * (seg / self.D_h) * 0.5 * rho * u * abs(u)
            dp = p - Pc
            q = 0.0
            if dp > 0:
                q = cd_of(0.0, p) * self.A_hole * math.sqrt(2.0 * rho * dp)
                for _ in range(3):   # the hole's Cd at its own Reynolds number
                    Re_h = rho * (q / (rho * self.A_hole)) * self.d_hole / self.mu
                    q = cd_of(Re_h, p) * self.A_hole * math.sqrt(2.0 * rho * dp)
            cds.append(q / (self.A_hole * math.sqrt(2.0 * rho * dp)) if dp > 0 and q > 0 else cd_of(0.0, p))
            mb -= q
            u_new = mb / (rho * A_ch)
            p += self.C_R * 0.5 * rho * (u * u - u_new * u_new)
            u = u_new
            holes.append(q)
        total = 2.0 * self.n_ports * self.scale * sum(holes)
        cd_hole = (sum(c * q for c, q in zip(cds, holes)) / sum(holes)) if sum(holes) > 0 else cds[0]
        return total, holes, cd_hole
 
 
def _ring_manifold_for(config, side, rho, mu, n_holes, d_hole):
    """The stream's ring network when the declared plate has back channels, else None."""
    plate = getattr(config.injector, "plate", None)
    if (plate is None or getattr(plate, "back", None) != "channels"
            or getattr(plate, "manifold_model", "plenum") != "ring_network"):
        return None
    from engine.core.injectors.layout import layout_from_config
    lay = layout_from_config(config, drawings=False, checks=False)
    if lay is None:
        return None
    ch = lay["passages"]["O" if side == "oxidizer" else "F"].get("channel")
    if not ch or not ch.get("flow_area", 0.0) > 0:
        return None
    w, h = float(ch["width"]), float(ch.get("flow_height") or ch["depth"])
    D_h = float(ch.get("hydraulic_diameter") or 4.0 * w * h / (2.0 * (w + h)))   # drawn section, else w x h
    return _RingManifold(
        n_holes=int(n_holes), n_ports=int(getattr(plate, "channel_inlets", 1) or 1),
        A_ch=float(ch["flow_area"]), D_h=D_h, r_ring=float(ch["r_center"]),
        A_hole=math.pi * d_hole ** 2 / 4.0, d_hole=d_hole, rho=rho, mu=mu,
        K_ent=float(getattr(plate, "channel_entry_K", 0.5)),
        C_R=float(getattr(plate, "channel_pressure_regain", 1.0)),
        roughness=float(getattr(plate, "channel_roughness", 3.2e-6)))
 
 
def _manifold_diagnostics(net_O, net_F, detail, mdot_O, mdot_F):
    """Per-hole flow along each channel branch (port first) and the element mixture ratios."""
    if net_O is None and net_F is None:
        return {"manifold_model": "plenum"}
    out = {"manifold_model": "ring_network"}
    per = {}
    for k, net, md in (("O", net_O, mdot_O), ("F", net_F, mdot_F)):
        holes = detail[k].get("holes") if net is not None else None
        if holes:
            mean = md / (2.0 * net.n_ports * net.scale * len(holes)) if md > 0 else float("nan")
            per[k] = list(holes)
            out[f"element_mass_flows_{k}"] = [float(q) for q in holes]
            out[f"element_flow_ratio_min_{k}"] = float(min(holes) / mean) if mean > 0 else float("nan")
            out[f"element_flow_ratio_max_{k}"] = float(max(holes) / mean) if mean > 0 else float("nan")
            out[f"manifold_branch_velocity_{k}"] = float(md / (2.0 * net.n_ports) / (net.rho * net.A_ch))
            out[f"manifold_ports_{k}"] = int(net.n_ports)
            out[f"Cd_eff_manifold_{k}"] = float(detail[k].get("cd_eff", float("nan")))
    if "O" in per and "F" in per and len(per["O"]) == len(per["F"]):
        mrs = [o / f if f > 0 else float("inf") for o, f in zip(per["O"], per["F"])]
        out["element_mixture_ratios"] = mrs
        out["element_mass_flows"] = [o + f for o, f in zip(per["O"], per["F"])]
        out["element_mixture_ratio_min"] = float(min(mrs))
        out["element_mixture_ratio_max"] = float(max(mrs))
    return out
 
 
def impingement_standoff_m(
    n_elements: float,
    spacing_O_m: float,
    spacing_F_m: float,
    angle_O_deg: float,
    angle_F_deg: float,
) -> float:
    """Axial distance from the injector FACE to where a doublet's two jets meet [m].
 
    The two streams of an unlike doublet sit on DIFFERENT pitch circles::
 
        D_pitch = n * spacing / pi          (spacing = arc between adjacent same-stream holes)
        dr      = |D_pitch_O - D_pitch_F| / 2
 
    ``dr`` is the RADIAL gap between the two rings at the face. Each jet is inclined toward
    the other and closes its share of that gap at rate ``tan(theta)`` per unit axial
    distance, so they meet at::
 
        L_imp = dr / (tan(theta_O) + tan(theta_F))
 
    NOTE the offset is ``dr``, not the element spacing. ``spacing`` is the circumferential
    distance between ADJACENT holes of the SAME stream around its own ring -- it says how
    densely one ring is populated and has nothing to do with how far the O and F jets of a
    single doublet are from each other. An earlier version used ``0.5*(s_O + s_F)`` here,
    which is the average ring density: dimensionally a length, physically the wrong one, and
    it under-reported the standoff (2.89 mm vs 4.35 mm on a real design, i.e. L/d 1.59 vs
    2.39) while also shortening the spray length that feeds the residence-time correction.
 
    Returns NaN when the rings coincide or both jets are axial -- in either case there is no
    crossing point.
    """
    import numpy as _np
 
    n = max(1.0, float(n_elements))
    dp_O = n * float(spacing_O_m) / _np.pi
    dp_F = n * float(spacing_F_m) / _np.pi
    dr = 0.5 * abs(dp_O - dp_F)
    tan_sum = float(
        _np.tan(_np.deg2rad(float(angle_O_deg))) + _np.tan(_np.deg2rad(float(angle_F_deg)))
    )
    if not _np.isfinite(tan_sum) or tan_sum <= 1e-9 or dr <= 0.0:
        return float("nan")
    val = dr / tan_sum
    return float(val) if _np.isfinite(val) else float("nan")
 
def momentum_ratio_R_from_bulk_velocities(
    rho_O: float,
    rho_F: float,
    v_O_bulk: float,
    v_F_bulk: float,
) -> float:
    """Jet momentum-flux ratio R = sqrt(rho_O*v_O^2 / (rho_F*v_F^2)), used by Layer 1.
 
    ``v_*_bulk`` are bulk speeds per stream: mdot / (rho * n_elements * A_jet). Since
    mdot = Cd A sqrt(2 rho dp), rho v^2 = 2 Cd^2 dp and R = (Cd_O/Cd_F) sqrt(dp_O/dp_F): the
    injector's Cd-dp split, with no density or hole size in it. It is not Rupe's mixing
    criterion; that is rupe_mixing_ratio.
    """
    if (
        rho_O <= 0
        or rho_F <= 0
        or not np.isfinite(v_O_bulk)
        or not np.isfinite(v_F_bulk)
        or v_F_bulk == 0.0
    ):
        return float("nan")
    num_mom = rho_O * v_O_bulk ** 2
    den_mom = rho_F * v_F_bulk ** 2
    if den_mom <= 0 or num_mom < 0:
        return float("nan")
    return float(np.sqrt(num_mom / den_mom))
 
 
def rupe_mixing_ratio(
    rho_O: float, v_O: float, d_O: float, rho_F: float, v_F: float, d_F: float,
) -> float:
    """Rupe's mixing ratio for an unlike doublet, M = rho_O v_O^2 d_O / (rho_F v_F^2 d_F).
 
    Rupe (JPL PR 20-195), restated by Elverum & Morey, JPL Memo 30-5 (1959) eq. 1: the most
    uniform mixture-ratio distribution is at M = 1 (stream momentum ratio equal to the diameter
    ratio). No angle term; the spray's direction is a separate quantity. M = R^2 d_O / d_F.
    """
    den = rho_F * v_F ** 2 * d_F
    if not (np.isfinite(den) and den > 0 and np.isfinite(v_O) and rho_O > 0 and d_O > 0):
        return float("nan")
    return float(rho_O * v_O ** 2 * d_O / den)
 
 
def injector_face_pressure(Pc: float, gamma: float, contraction_ratio: float) -> float:
    """Static pressure at the injector face for a nozzle-inlet stagnation pressure ``Pc``.
 
    Heat addition in a constant-area chamber (Rayleigh line) costs stagnation pressure between
    the face, where the gas is at rest, and the nozzle inlet at M_c, the subsonic root of
    A_c/A_t: p_face / p0 = (1 + gamma M_c^2) / (1 + (gamma-1)/2 M_c^2)^(gamma/(gamma-1)).
    Sutton & Biblarz ch. 3: ~1.05 at A_c/A_t = 2, under 1 % above ~5. Not applied by the solve
    yet (dp_inj = P_inj - Pc); it needs the CEA gamma the chamber solver holds.
    """
    g, cr = float(gamma), float(contraction_ratio)
    if not (np.isfinite(Pc) and g > 1.0 and cr > 1.0):
        return float(Pc)
    ex = (g + 1.0) / (2.0 * (g - 1.0))
    area = lambda M: ((2.0 / (g + 1.0)) * (1.0 + 0.5 * (g - 1.0) * M * M)) ** ex / M  # noqa: E731
    lo, hi = 1e-12, 1.0
    for _ in range(200):                       # A/A* falls monotonically on the subsonic branch
        mid = 0.5 * (lo + hi)
        lo, hi = (mid, hi) if area(mid) > cr else (lo, mid)
    M2 = (0.5 * (lo + hi)) ** 2
    return float(Pc * (1.0 + g * M2) / (1.0 + 0.5 * (g - 1.0) * M2) ** (g / (g - 1.0)))
 
 
class ImpingingInjector(InjectorModel):
    """Twin-jet impinging injector solver."""
 
    injector_config: ImpingingInjectorConfig
 
    def __init__(self, engine_config: PintleEngineConfig):
        super().__init__(engine_config)
        injector_cfg = engine_config.injector
        if not isinstance(injector_cfg, ImpingingInjectorConfig):
            raise TypeError("ImpingingInjector requires ImpingingInjectorConfig")
        self.injector_config = injector_cfg
 
    def solve(
        self,
        P_tank_O: float,
        P_tank_F: float,
        Pc: float,
    ) -> Tuple[float, float, Dict[str, Any]]:
        config = self.engine_config
        geometry = self.injector_config.geometry
 
        # The land the plate actually gives when the config asks for it (l_over_d_source:
        # plate); otherwise the declared blocks, unchanged.
        discharge_O = effective_discharge(config, "oxidizer")
        discharge_F = effective_discharge(config, "fuel")
        feed_O = config.feed_system["oxidizer"]
        feed_F = config.feed_system["fuel"]
        spray_cfg = config.spray
        fluids = config.fluids
 
        rho_O = fluids["oxidizer"].density
        mu_O = fluids["oxidizer"].viscosity
        sigma_O = fluids["oxidizer"].surface_tension
 
        rho_F = fluids["fuel"].density
        mu_F = fluids["fuel"].viscosity
        sigma_F = fluids["fuel"].surface_tension
 
        # Inlet temperatures for Cd(Re): match pintle — use fluid properties from config when set
        T_tank_O = getattr(fluids["oxidizer"], "temperature", 90.0)
        T_tank_F = getattr(fluids["fuel"], "temperature", 300.0)
 
        A_O = geometry.oxidizer.n_elements * np.pi * (geometry.oxidizer.d_jet / 2.0) ** 2
        A_F = geometry.fuel.n_elements * np.pi * (geometry.fuel.d_jet / 2.0) ** 2
        d_hyd_O = geometry.oxidizer.d_jet
        d_hyd_F = geometry.fuel.d_jet
 
        mdot_O = 0.1
        mdot_F = 0.1
 
        if float(config.solver.closure.Cd_reduction_factor) != 1.0:
            assume("solver.closure.Cd_reduction_factor", 1.0,
                   reason="Cd is orifice geometry and Re; a spray-constraint violation is reported, not traded for Cd")
        for _side, _dc in (("oxidizer", discharge_O), ("fuel", discharge_F)):
            if cd_inf_from_inlet_geometry(_dc) is None and getattr(_dc, "use_geometry_cd", False):
                assume(f"discharge.{_side}.inlet_geometry", "thin plate",
                       reason="no inlet declared: thin-plate Cd, not the drilled L/d passage the layout draws "
                              "(sharp, L/d 4: 0.79, Lichtarowicz 1965)")
        Cd_O_eff = cd_inf_from_orifice_diameter(d_hyd_O, discharge_O)
        Cd_F_eff = cd_inf_from_orifice_diameter(d_hyd_F, discharge_F)
 
        diagnostics = {
            "iterations": 0,
            "constraints_satisfied": False,
            "violations": [],
            "J": None,
            "TMR": None,
            "theta": None,
            "We_O": None,
            "We_F": None,
            "D32_O": None,
            "D32_F": None,
            "x_star": None,
            "injector_type": "impinging",
            "turbulence_intensity_O": None,
            "turbulence_intensity_F": None,
            "turbulence_length_O": None,
            "turbulence_length_F": None,
            "turbulence_intensity_mix": None,
        }
 
        # Typical opposing doublet: each side's ``impingement_angle`` is measured from chamber axis toward
        # the impingement plane; separation between jet centerlines is θ_O + θ_F (schema: included angle per stream).
        _imp_sep_deg = float(geometry.oxidizer.impingement_angle + geometry.fuel.impingement_angle)
        _imp_sep_deg = float(np.clip(_imp_sep_deg, 1.0, 179.0))
        imp_angle_rad = np.deg2rad(_imp_sep_deg)
 
        def _injector_turbulence_fields(uo: float, uf: float) -> Dict[str, float]:
            """Same Reynolds-based turbulence model as pintle (shear-layer mixing intensity)."""
            Re_Oloc = calculate_reynolds_number(rho_O, uo, d_hyd_O, mu_O)
            Re_Floc = calculate_reynolds_number(rho_F, uf, d_hyd_F, mu_F)
            ti_O = 0.16 * (Re_Oloc ** -0.125) if Re_Oloc > 0 else 0.1
            ti_F = 0.16 * (Re_Floc ** -0.125) if Re_Floc > 0 else 0.1
            ti_O = float(np.clip(ti_O, 0.02, 0.3))
            ti_F = float(np.clip(ti_F, 0.02, 0.3))
            tl_O = float(0.07 * d_hyd_O)
            tl_F = float(0.07 * d_hyd_F)
            velocity_total = max(uo + uf, 1e-6)
            ti_mix = (ti_O * uo + ti_F * uf) / velocity_total
            ti_mix = float(np.clip(ti_mix, 0.02, 0.35))
            return {
                "turbulence_intensity_O": ti_O,
                "turbulence_intensity_F": ti_F,
                "turbulence_length_O": tl_O,
                "turbulence_length_F": tl_F,
                "turbulence_intensity_mix": ti_mix,
            }
 
        def _converge_feed_orifice_coupling(
            mdot_o0: float,
            mdot_f0: float,
            cd_eff_o: float,
            cd_eff_f: float,
        ) -> Tuple[
            float,
            float,
            float,
            float,
            float,
            float,
            float,
            float,
            float,
            float,
            int,
        ]:
            """Feed loss, inlet pressure and orifice flow at this Pc, each stream by ``_stream_flow``.
 
            With Pc fixed the streams do not interact, so each is its own bracketed root. The relaxed
            fixed point this replaces stalled at its 150-step cap whenever a probed Pc sat above a
            tank, and stopped ~2e-6 short of the root otherwise.
 
            Returns
            -------
            mdot_o, mdot_f, cd_o, cd_f, Pi_o, Pi_f, dpf_o, dpf_f, dpi_o, dpi_f, evaluations
            """
 
            def _feed_o(m):
                return delta_p_feed(m, rho_O, feed_O, P_tank_O, mu=mu_O)
 
            def _feed_f(m):
                dp = delta_p_feed(m, rho_F, feed_F, P_tank_F, mu=mu_F)
                if config.regen_cooling is not None and config.regen_cooling.enabled:
                    dp += delta_p_regen_channels(m, rho_F, mu_F, config.regen_cooling, P_tank_F)
                return dp
 
            mo, Cdo, Pi_o, dpf_o, dpi_o, n_o = _stream_flow(
                Pc, P_tank_O, rho_O, A_O, d_hyd_O, mu_O, discharge_O, T_tank_O, cd_eff_o, _feed_o,
                cav=cav_O, network=net_O, detail=net_detail["O"],
            )
            mf, Cdf, Pi_f, dpf_f, dpi_f, n_f = _stream_flow(
                Pc, P_tank_F, rho_F, A_F, d_hyd_F, mu_F, discharge_F, T_tank_F, cd_eff_f, _feed_f,
                cav=cav_F, network=net_F, detail=net_detail["F"],
            )
            fp_it = n_o + n_f
 
            mdot_bn_o = float(Cdo * A_O * np.sqrt(2.0 * rho_O * dpi_o)) if dpi_o > 0 else 0.0
            mdot_bn_f = float(Cdf * A_F * np.sqrt(2.0 * rho_F * dpi_f)) if dpi_f > 0 else 0.0
 
            _LOG.debug(
                "impinging feed-orifice closure: fp_iters=%d mdot_solver_O=%s mdot_from_Bernoulli_O=%s "
                "mdot_solver_F=%s mdot_from_Bernoulli_F=%s",
                fp_it,
                repr(mo),
                repr(mdot_bn_o),
                repr(mf),
                repr(mdot_bn_f),
            )
            if _PRINT_FEED_ORIFICE_CLOSURE:
                print(
                    "impinging feed-orifice closure",
                    f"fp_iters={fp_it}",
                    f"mdot_solver_O={mo}",
                    f"mdot_from_Bernoulli_O={mdot_bn_o}",
                    f"mdot_solver_F={mf}",
                    f"mdot_from_Bernoulli_F={mdot_bn_f}",
                )
 
            return mo, mf, Cdo, Cdf, Pi_o, Pi_f, dpf_o, dpf_f, dpi_o, dpi_f, fp_it
 
        # Cavitation limit on every hole (Nurick 1976); a no-op wherever K > (Cd/Cc)^2.
        from engine.core.discharge import contraction_coefficient, inlet_radius_ratio_of
 
        def _cav(side, dc):
            pv = getattr(fluids[side], "vapor_pressure", None)
            if pv is None or not np.isfinite(float(pv)):
                return None
            return float(pv), contraction_coefficient(inlet_radius_ratio_of(dc))
 
        cav_O, cav_F = _cav("oxidizer", discharge_O), _cav("fuel", discharge_F)
        # Back channels as dividing-flow rings (None: the manifold is a plenum at one pressure).
        net_O = _ring_manifold_for(config, "oxidizer", rho_O, mu_O,
                                   geometry.oxidizer.n_elements, geometry.oxidizer.d_jet)
        net_F = _ring_manifold_for(config, "fuel", rho_F, mu_F,
                                   geometry.fuel.n_elements, geometry.fuel.d_jet)
        net_detail = {"O": {}, "F": {}}
 
        feed_orifice_fp_last = 0
 
        # One pass. Cd is orifice geometry and Reynolds number (Lichtarowicz 1965; SP-8089), not a
        # lever for the spray constraints: an x* or We violation is reported in diagnostics. The old
        # loop shrank Cd by Cd_reduction_factor per violation, which also lengthened x*.
        mdot_O, mdot_F, Cd_O, Cd_F, P_inj_O, P_inj_F, delta_p_feed_O, delta_p_feed_F, delta_p_inj_O, delta_p_inj_F, feed_orifice_fp_last = _converge_feed_orifice_coupling(
            mdot_O,
            mdot_F,
            Cd_O_eff,
            Cd_F_eff,
        )
 
        u_O = mdot_O / (rho_O * A_O) if A_O > 0 else 0.0
        u_F = mdot_F / (rho_F * A_F) if A_F > 0 else 0.0
 
        # Impingement relative velocity: law of cosines on the two jet velocity vectors,
        # whose centerlines are separated by the included angle ``imp_angle_rad``. This is the
        # velocity at which the unlike jets collide and shear into a sheet, and it is what
        # drives both breakup (Weber) and the residence/evaporation length (x*).
        u_rel = float(
            np.sqrt(u_O ** 2 + u_F ** 2 - 2 * u_O * u_F * np.cos(imp_angle_rad))
        )
        turb_fields = _injector_turbulence_fields(u_O, u_F)
 
        # Representative chamber gas density for the aerodynamic (breakup) Weber number.
        # rho_g = Pc / (R_gas · T_gas); R_gas and T_gas are configured representative
        # combustion-gas values because the injector solve runs before the CEA chamber state
        # is available. Gas density (not liquid inertia) governs primary atomization.
        rho_gas = float(
            max(Pc / (spray_cfg.smd.chamber_gas_R * spray_cfg.smd.chamber_gas_T), 1e-6)
        )
 
        J = momentum_flux_ratio(rho_O, u_O, rho_F, u_F)
        MR = mdot_O / mdot_F if mdot_F > 0 else np.inf
        TMR = thrust_momentum_ratio(J, MR)
 
        if spray_cfg.spray_angle.model == "J":
            theta = spray_angle_from_J(J, spray_cfg.spray_angle.k, spray_cfg.spray_angle.n)
        else:
            theta = spray_angle_from_TMR(TMR)
 
        Oh_O = ohnesorge_number(mu_O, rho_O, sigma_O, geometry.oxidizer.d_jet)
        Oh_F = ohnesorge_number(mu_F, rho_F, sigma_F, geometry.fuel.d_jet)
 
        # Impinging doublets always use Ingebo (u_rel-based aerodynamic breakup). Legacy Lefebvre
        # on this path used uncalibrated C/m and reported bogus ~1 µm D32 when spray.smd.model was
        # left at "lefebvre" from pintle-era YAML (see dispatch impinging bindings).
        We_O = weber_number(rho_gas, u_rel, geometry.oxidizer.d_jet, sigma_O)
        We_F = weber_number(rho_gas, u_rel, geometry.fuel.d_jet, sigma_F)
        # NACA TN 4222: each stream at its own jet velocity, entering gas at rest (dV = Vj).
        _smd = spray_cfg.smd
        _prop = bool(getattr(_smd, "smd_property_scaling", True))
        _scale = float(getattr(_smd, "smd_scale", 1.0) or 1.0)
        from engine.core.spray import tn4222_transfer_model
        _transfer = tn4222_transfer_model(_smd)
        # A stream that does not flow (tank at or below Pc) has no spray: NaN, not an error.
        D32_O = (_scale * smd_impinging_tn4222(
            geometry.oxidizer.d_jet, u_O, rho_O, mu_O, sigma_O, rho_gas,
            scale_properties=_prop, transfer=_transfer) if u_O > 0 else float("nan"))
        D32_F = (_scale * smd_impinging_tn4222(
            geometry.fuel.d_jet, u_F, rho_F, mu_F, sigma_F, rho_gas,
            scale_properties=_prop, transfer=_transfer) if u_F > 0 else float("nan"))
        # Measured D32 (measurements.d32_*_um) replaces the correlation where it was measured.
        _mO, _mF = getattr(_smd, "d32_measured_O", None), getattr(_smd, "d32_measured_F", None)
        if _mO and u_O > 0:
            D32_O = float(_mO)
        if _mF and u_F > 0:
            D32_F = float(_mF)
 
        # Evaporation constant from propellant properties + chamber state, not a single
        # hardcoded K. Falls back to the legacy tau = K*D32^2 when the derivation cannot run
        # (missing latent heat / boiling point on a custom fluid), so old configs still work.
        _ev = spray_cfg.evaporation
        _use_derived = getattr(_ev, "model", "derived") == "derived"
        k_evap_O = k_evap_F = float("nan")
        if _use_derived:
            k_evap_O = evaporation_constant_m2_s(
                Tc=spray_cfg.smd.chamber_gas_T, Pc=Pc, rho_g=rho_gas, rho_l=rho_O,
                L_vap=float(getattr(fluids["oxidizer"], "latent_heat", 0.0) or 0.0),
                T_boil=float(getattr(fluids["oxidizer"], "boiling_point", 0.0) or 0.0),
                cp_g=float(getattr(_ev, "cp_gas", 2200.0)),
                C_evap=float(getattr(_ev, "C_evap", 1.562)),
            )
            k_evap_F = evaporation_constant_m2_s(
                Tc=spray_cfg.smd.chamber_gas_T, Pc=Pc, rho_g=rho_gas, rho_l=rho_F,
                L_vap=float(getattr(fluids["fuel"], "latent_heat", 0.0) or 0.0),
                T_boil=float(getattr(fluids["fuel"], "boiling_point", 0.0) or 0.0),
                cp_g=float(getattr(_ev, "cp_gas", 2200.0)),
                C_evap=float(getattr(_ev, "C_evap", 1.562)),
            )
        if np.isfinite(k_evap_O) and np.isfinite(k_evap_F):
            tau_evap_O = tau_evap_from_k(D32_O, k_evap_O)
            tau_evap_F = tau_evap_from_k(D32_F, k_evap_F)
        else:
            tau_evap_O = tau_evap(D32_O, _ev.K)
            tau_evap_F = tau_evap(D32_F, _ev.K)
 
        # TRANSPORT velocity, not the jet-to-jet relative velocity. u_rel shears the sheet
        # (it belongs in the Ingebo Weber number above); what carries droplets DOWN the
        # chamber is the momentum-weighted axial resultant of the collision. Using u_rel here
        # over-predicted x* by ~2x on a measured design.
        u_axial = spray_axial_velocity(
            mdot_O=mdot_O, u_O=u_O, theta_O_deg=geometry.oxidizer.impingement_angle,
            mdot_F=mdot_F, u_F=u_F, theta_F_deg=geometry.fuel.impingement_angle,
        )
        u_transport = u_axial if np.isfinite(u_axial) and u_axial > 0 else u_rel
        _xs = [v for v in (xstar(u_transport, tau_evap_O), xstar(u_transport, tau_evap_F)) if np.isfinite(v)]
        x_star = max(_xs) if _xs else float("nan")   # NaN only when neither stream flows
 
        # ---- Impinging-doublet geometry (standoff + ring pitch) -------------------------------
        # Each stream's ``impingement_angle`` is the jet inclination from the chamber axis; the
        # two jets of a doublet lean toward each other and collide on the bisector. The offset
        # they close is RADIAL -- the gap between the two pitch circles --
        #     dr    = |D_pitch_O - D_pitch_F| / 2,   D_pitch = n·spacing/π
        #     L_imp = dr / (tan θ_O + tan θ_F)
        # (purely geometric: each jet closes its half of the offset at rate tan θ). This block
        # used to say ``s_pair`` -- the average RING DENSITY, 0.5*(s_O + s_F). That is
        # dimensionally a length and physically the wrong one; see impingement_standoff_m,
        # which has been computing dr correctly. ``s_pair`` survives below as a diagnostic
        # only. Nothing reads it for the standoff.
        theta_O = float(np.deg2rad(geometry.oxidizer.impingement_angle))
        theta_F = float(np.deg2rad(geometry.fuel.impingement_angle))
        s_O = float(getattr(geometry.oxidizer, "spacing", 0.0) or 0.0)
        s_F = float(getattr(geometry.fuel, "spacing", 0.0) or 0.0)
        s_pair = 0.5 * (s_O + s_F)
        tan_sum = float(np.tan(theta_O) + np.tan(theta_F))
        L_imp = impingement_standoff_m(
            0.5 * (geometry.oxidizer.n_elements + geometry.fuel.n_elements),
            s_O, s_F,
            geometry.oxidizer.impingement_angle,
            geometry.fuel.impingement_angle,
        )
        nO_geom = max(1, int(geometry.oxidizer.n_elements))
        nF_geom = max(1, int(geometry.fuel.n_elements))
        D_pitch_O = float(nO_geom * s_O / np.pi)
        D_pitch_F = float(nF_geom * s_F / np.pi)
        # Circumferential clearance between adjacent same-stream orifices (must stay positive).
        gap_O = float(s_O - geometry.oxidizer.d_jet)
        gap_F = float(s_F - geometry.fuel.d_jet)
        # Axial length the spray needs before it is fully vaporized, measured from the face:
        # standoff to impingement plus the droplet evaporation length x*.
        # Sheet breakup length between impingement and droplet formation. Mass conservation
        # in a radially spreading sheet gives h(r) = d^2/(4r); Kelvin-Helmholtz growth on a
        # thin liquid sheet in gas gives t_b ~ (h/u)*sqrt(rho_l/rho_g), so L_b = u*t_b. At the
        # Weber numbers here (3e4-6e4, far past the ~2e3 regime transition) this comes out
        # small next to L_imp -- the sheet shreds almost immediately -- but it is the step
        # between "jets meet" and "droplets exist", so the spray length is not complete
        # without it.
        _d_avg = 0.5 * (float(geometry.oxidizer.d_jet) + float(geometry.fuel.d_jet))
        _u_sheet = max(float(u_rel), 1e-6)
        _rho_l_avg = 0.5 * (float(rho_O) + float(rho_F))
        # rho_gas carries a 1e-6 floor for solver robustness, and L_b ~ sqrt(rho_l/rho_g),
        # so at that floor the sqrt is ~2.8e4 and L_b explodes -- during Pc bracketing that
        # drove the spray length past the whole chamber, collapsed eta_Lstar, and made the
        # chamber solve report "Supply < Demand at all Pc". Bound it to a physically sane
        # multiple of the jet diameter: at these Weber numbers the sheet shreds in ~1 d_jet,
        # and anything past ~20 is not a doublet sheet any more.
        if np.isfinite(L_imp) and L_imp > 0 and _d_avg > 0 and rho_gas > 1e-4:
            _h_sheet = (_d_avg ** 2) / (4.0 * L_imp)
            L_b = float(_h_sheet * np.sqrt(_rho_l_avg / rho_gas))
            L_b = float(np.clip(L_b, 0.0, 20.0 * _d_avg))
        else:
            L_b = 0.0
        vaporization_length_total = float(
            (L_imp if np.isfinite(L_imp) else 0.0) + L_b + x_star
        )
 
        constraints_ok, violations = check_spray_constraints(We_O, We_F, x_star, spray_cfg)
 
        diagnostics.update(
            {
                "iterations": 1,
                "constraints_satisfied": constraints_ok,
                "violations": violations,
                "J": J,
                "TMR": TMR,
                "theta": theta,
                "u_rel": float(u_rel),
                "rho_gas_breakup": float(rho_gas),
                "Oh_O": float(Oh_O),
                "Oh_F": float(Oh_F),
                "We_O": We_O,
                "We_F": We_F,
                "D32_O": D32_O,
                "D32_F": D32_F,
                "x_star": x_star,
                "u_axial_spray": float(u_transport),
                "k_evap_O": float(k_evap_O),
                "k_evap_F": float(k_evap_F),
                "tau_evap_O": float(tau_evap_O),
                "tau_evap_F": float(tau_evap_F),
                "L_imp": L_imp,
                "L_sheet_breakup": float(L_b),
                "D_pitch_O": D_pitch_O,
                "D_pitch_F": D_pitch_F,
                "s_pair": float(s_pair),
                "element_gap_O": gap_O,
                "element_gap_F": gap_F,
                "vaporization_length_total": vaporization_length_total,
                # The spray-zone tau_res switch travels WITH the spray diagnostics.
                # calculate_combustion_efficiency_advanced() only receives a
                # CombustionEfficiencyConfig, so it used to read this off that object --
                # where the field does not exist -- and a bare except swallowed the
                # mismatch, making the switch permanently False and unreachable.
                "apply_tau_res_correction": bool(
                    getattr(spray_cfg.evaporation, "apply_tau_res_correction", False)
                ),
                "impingement_angle_deg": np.rad2deg(imp_angle_rad),
                "V_rel": float(u_rel),
                "breakup_multiplier": 1.0,
                "penetration_multiplier": 1.0,
                **turb_fields,
                # Discharge coefficients
                "Cd_O": float(Cd_O),
                "Cd_F": float(Cd_F),
            }
        )
 
        # Diagnostics must match last feed–orifice coupling (recomputing from ṁ duplicates float path and drifted Cd).
        delta_p_feed_O_final = float(delta_p_feed_O)
        delta_p_feed_F_final = float(delta_p_feed_F)
        P_inj_O_final = float(P_inj_O)
        P_inj_F_final = float(P_inj_F)
        delta_p_inj_O_final = float(delta_p_inj_O)
        delta_p_inj_F_final = float(delta_p_inj_F)
 
        mdot_from_bernoulli_O = (
            float(Cd_O * A_O * np.sqrt(2.0 * rho_O * delta_p_inj_O_final))
            if delta_p_inj_O_final > 0
            else 0.0
        )
        mdot_from_bernoulli_F = (
            float(Cd_F * A_F * np.sqrt(2.0 * rho_F * delta_p_inj_F_final))
            if delta_p_inj_F_final > 0
            else 0.0
        )
 
        u_O_final = mdot_O / (rho_O * A_O) if A_O > 0 else 0.0
        u_F_final = mdot_F / (rho_F * A_F) if A_F > 0 else 0.0
        # Impingement relative velocity (vector difference of the two jet velocities separated by
        # the included angle); reported as V_rel for the chamber/combustion models.
        u_rel_final = float(
            np.sqrt(
                u_O_final ** 2
                + u_F_final ** 2
                - 2 * u_O_final * u_F_final * np.cos(imp_angle_rad)
            )
        )
        turb_final = _injector_turbulence_fields(u_O_final, u_F_final)
 
        # Momentum metrics on v = mdot / (rho * n_elements * A_jet), A_jet = pi*(d_jet/2)^2:
        # R, the Cd-dp split, and Rupe's M, the mixing criterion.
        n_O = max(1, int(geometry.oxidizer.n_elements))
        n_F = max(1, int(geometry.fuel.n_elements))
        djo = float(geometry.oxidizer.d_jet)
        djf = float(geometry.fuel.d_jet)
        A_jet_O = np.pi * (djo / 2.0) ** 2
        A_jet_F = np.pi * (djf / 2.0) ** 2
        denom_O = rho_O * float(n_O) * A_jet_O
        denom_F = rho_F * float(n_F) * A_jet_F
        v_O_bulk = mdot_O / denom_O if denom_O > 0 else np.nan
        v_F_bulk = mdot_F / denom_F if denom_F > 0 else np.nan
        momentum_ratio_R = momentum_ratio_R_from_bulk_velocities(
            float(rho_O), float(rho_F), float(v_O_bulk), float(v_F_bulk)
        )
        rupe_M = rupe_mixing_ratio(
            float(rho_O), float(v_O_bulk), djo, float(rho_F), float(v_F_bulk), djf
        )
 
        mom_update: Dict[str, Any] = {
            "A_jet_O": float(A_jet_O),
            "A_jet_F": float(A_jet_F),
            "momentum_ratio_n_elements_O": int(n_O),
            "momentum_ratio_n_elements_F": int(n_F),
            "d_jet_O": float(djo),
            "d_jet_F": float(djf),
        }
        if np.isfinite(v_O_bulk):
            mom_update["v_O_bulk"] = float(v_O_bulk)
        if np.isfinite(v_F_bulk):
            mom_update["v_F_bulk"] = float(v_F_bulk)
        if rho_O > 0:
            mom_update["rho_O_momentum"] = float(rho_O)
        if rho_F > 0:
            mom_update["rho_F_momentum"] = float(rho_F)
        if np.isfinite(momentum_ratio_R) and momentum_ratio_R > 0:
            mom_update["momentum_ratio_R"] = momentum_ratio_R
        if np.isfinite(rupe_M) and rupe_M > 0:
            mom_update["rupe_M"] = rupe_M
 
        diagnostics.update(
            {
                # Required by chamber_solver / combustion efficiency (parity with pintle)
                "u_O": float(u_O_final),
                "u_F": float(u_F_final),
                "V_rel": u_rel_final,
                **turb_final,
                "P_injector_O": float(P_inj_O_final),
                "P_injector_F": float(P_inj_F_final),
                "delta_p_injector_O": float(delta_p_inj_O_final),
                "delta_p_injector_F": float(delta_p_inj_F_final),
                "delta_p_feed_O": float(delta_p_feed_O_final),
                "delta_p_feed_F": float(delta_p_feed_F_final),
                "feed_orifice_coupling_iterations": int(feed_orifice_fp_last),
                **_manifold_diagnostics(net_O, net_F, net_detail, mdot_O, mdot_F),
                "mdot_from_bernoulli_O": float(mdot_from_bernoulli_O),
                "mdot_from_bernoulli_F": float(mdot_from_bernoulli_F),
                **mom_update,
            }
        )
 
        # Geometric vs effective flow areas (A_eff = Cd × A_geom); reuse converged Cd_O/Cd_F
        A_eff_O, A_eff_F, eff_warns = effective_flow_areas_from_cd(diagnostics, A_O, A_F)
        diagnostics["A_geom_O"] = float(A_O)
        diagnostics["A_geom_F"] = float(A_F)
        diagnostics["A_eff_O"] = float(A_eff_O)
        diagnostics["A_eff_F"] = float(A_eff_F)
        if eff_warns:
            merge_effective_area_warnings(diagnostics, eff_warns)
 
        return mdot_O, mdot_F, diagnostics