"""Driver: regenerates every number in report SR-2026-009. Usage: python driver.py [--quick] (default) the full 45.2M-candidate design-space search. On one 192-core node about five minutes; on two cores roughly 50 hours. --quick a coarser grid, a few minutes on two cores, which reproduces the structural conclusions and the recommended design but not the exact feasible count. Requires sensitivity.py, campaign.py and optimize4d.py alongside; all four files are reproduced in full as listings in the report appendix. Deterministic: no random draws are used anywhere. Writes: design_space.csv, pareto_5d.csv, design_space_meta.json, de_haller_boundary.csv, reconciliation.csv, verification.csv, validation.csv. The figures are drawn from these files by plotting code that is not part of this archive. """ import json, math, os, sys from multiprocessing import Pool sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import numpy as np, pandas as pd import campaign as cp import optimize4d as op from sensitivity import BASE QUICK = "--quick" in sys.argv U0, PHI0, PSI0 = BASE["U_rms"], BASE["phi"], BASE["psi"] GC, CPC = 1.4, 1005.0 FT = 0.3048 DHMIN = op.DE_HALLER_MIN if QUICK: NS = list(range(8, 41, 2)); UR = np.arange(0.45, 1.85, 0.10) PR = np.arange(0.25, 2.05, 0.10); PSI = np.arange(0.12, 0.52, 0.04) NU = (1, 3, 5, 8) else: NS = list(range(6, 41)); UR = np.arange(0.45, 1.85, 0.025) PR = np.arange(0.25, 2.05, 0.025); PSI = np.arange(0.12, 0.52, 0.01) NU = (1, 2, 3, 4, 5, 6, 8, 10) # ---- 1. the diffusion boundary, in closed form and numerically ------------------- def phi_min_closed(psi, dh=DHMIN): """Lowest flow coefficient the diffusion limit permits at this stage loading. cos(arctan x) = 1/sqrt(1+x^2), so with tan(b1) = 1/phi and tan(b2) = (1-psi)/phi the de Haller number satisfies DH^2 = (phi^2 + (1-psi)^2) / (phi^2 + 1), which inverts to phi_min = sqrt((DH^2 - (1-psi)^2) / (1 - DH^2)). The argument is negative, and no flow coefficient is forbidden, when DH >= 1 - psi, i.e. psi <= 1 - DH. """ num = dh ** 2 - (1 - psi) ** 2 return math.sqrt(num / (1 - dh ** 2)) if num > 0 else float("nan") def boundary(): """The (phi, psi) feasibility curve, closed form checked against bisection.""" rows = [] for psi in np.arange(0.10, 0.52, 0.005): closed = phi_min_closed(psi) # independent numerical root of de_haller(phi, psi) = DHMIN lo, hi, num = 1e-3, 5.0, float("nan") if (op.de_haller(lo, psi) - DHMIN) * (op.de_haller(hi, psi) - DHMIN) < 0: for _ in range(200): mid = 0.5 * (lo + hi) if op.de_haller(mid, psi) < DHMIN: lo = mid else: hi = mid num = 0.5 * (lo + hi) rows.append(dict(psi=psi, phi_min=num, phi_min_closed=closed, free=not np.isfinite(closed))) BND = pd.DataFrame(rows) BND.to_csv("de_haller_boundary.csv", index=False) return BND # ---- 2. the design-space search --------------------------------------------------- def evaluate(arg): ns, ur, pr, ps, n = arg phi, U = PHI0 * pr, U0 * ur if op.de_haller(phi, ps) < DHMIN: return None try: o = op.objectives(ns, n, phi, U, psi=ps) except Exception: return None if o is None: return None return {k: o[k] for k in ("n_stages", "n_units", "phi", "psi", "U_rms", "de_haller", "p_bar", "mass_total", "fuel_kgh", "h_min_mm", "worst_margin", "active_at_limit")} def search(): ncore = int(os.environ.get("SLURM_CPUS_ON_NODE", os.cpu_count() or 1)) args = [(ns, ur, pr, ps, n) for ns in NS for ur in UR for pr in PR for ps in PSI for n in NU] with Pool(ncore) as pool: out = [r for r in pool.imap_unordered(evaluate, args, chunksize=2000) if r] D = pd.DataFrame(out) D.to_csv("design_space.csv", index=False) P, best = [], math.inf for _, r in D.sort_values("mass_total").iterrows(): if r.fuel_kgh < best - 1e-12: P.append(r) best = r.fuel_kgh PF = pd.DataFrame(P) PF.to_csv("pareto_5d.csv", index=False) b_m = D.loc[D.mass_total.idxmin()] b_f = D.loc[D.fuel_kgh.idxmin()] meta = dict(candidates=len(args), feasible=len(D), quick=QUICK, pareto_points=len(PF), psi_range=[float(D.psi.min()), float(D.psi.max())], phi_range=[float(D.phi.min()), float(D.phi.max())], stage_range=[int(D.n_stages.min()), int(D.n_stages.max())], mass_range=[float(D.mass_total.min()), float(D.mass_total.max())], fuel_range=[float(D.fuel_kgh.min()), float(D.fuel_kgh.max())], lightest={k: (float(b_m[k]) if k != "active_at_limit" else b_m[k]) for k in D.columns}, least_fuel={k: (float(b_f[k]) if k != "active_at_limit" else b_f[k]) for k in D.columns}) json.dump(meta, open("design_space_meta.json", "w"), indent=1) return D, PF, meta # ---- 2b. the frozen-stage-loading comparison --------------------------------------- def frozen_search(): """Repeat the search with psi held at the parent package's value. This is the comparison that settles what freeing the stage loading actually buys. An earlier draft of the report claimed that freezing psi is what made a previous phase conclude the requirement infeasible; this search refutes that, and shipping it lets a reader check the refutation rather than take it on trust. Same grid in the other four variables, same objective function. """ ncore = int(os.environ.get("SLURM_CPUS_ON_NODE", os.cpu_count() or 1)) args = [(ns, ur, pr, PSI0, n) for ns in NS for ur in UR for pr in PR for n in NU] with Pool(ncore) as pool: out = [r for r in pool.imap_unordered(evaluate, args, chunksize=2000) if r] F = pd.DataFrame(out) F.to_csv("frozen_psi_search.csv", index=False) return F, len(args) # ---- 3. verification --------------------------------------------------------------- def verification(best): """Residuals of the balances the solver imposes, plus three independent checks: the de Haller closed form against explicit velocity triangles, the boundary closed form against numerical bisection, and the optimiser's system mass against a hand re-solve at the sizing corner.""" alt_hi = 60000 * FT d = op.point(int(best.n_stages), int(best.n_units), float(best.phi), float(best.U_rms), alt_hi, 0.70, 10.0, psi=float(best.psi)) amb = op.mission(alt_hi, 0.70) gh, cph = BASE["gam_h"], BASE["cp_h"] M = 0.70 Tt2 = amb["T_inf"] * (1 + 0.2 * M ** 2) pt2 = amb["p_inf"] * (1 + 0.2 * M ** 2) ** 3.5 * BASE["pi_inlet"] w_c = CPC * Tt2 * (d["PR"] ** ((GC - 1) / GC) - 1) / BASE["eta_c"] Tt3 = Tt2 + w_c / CPC dh = cph * (BASE["TIT"] - Tt3) f = dh / (BASE["eta_b"] * BASE["LHV"] - dh) m5 = (1 - d["b"]) * (1 + d["f"]) dTt = w_c / (BASE["eta_mech"] * m5 * cph) rows = [("Fuel-air ratio from the burner balance", abs(f - d["f"]) / d["f"]), ("Delivery temperature from the compressor work", abs(Tt3 - d["T_del"]) / d["T_del"]), ("Shaft power balance", abs(BASE["eta_mech"] * m5 * cph * dTt - w_c) / w_c), ("Inlet continuity", abs(d["m_cap_unit"] / (d["rho1"] * d["Vx"]) - d["A_in"]) / d["A_in"])] def march(nsub): Tt, pt = Tt2, pt2 dT = (Tt3 - Tt2) / nsub n_ = GC / (GC - 1) * BASE["eta_c"] for _ in range(nsub): pt = pt * (1 + dT / Tt) ** n_ Tt += dT return pt p1 = march(1) for nsub in (16, 256): rows.append((f"Polytropic march refinement, {nsub} substeps", abs(march(nsub) - p1) / p1)) def dh_triangle(phi, psi): Vx = 1.0 U = Vx / phi return math.hypot(Vx, U - psi * U) / math.hypot(Vx, U) rows.append(("de Haller, closed form against velocity triangles", abs(op.de_haller(float(best.phi), float(best.psi)) - dh_triangle(float(best.phi), float(best.psi))) / dh_triangle(float(best.phi), float(best.psi)))) BND = pd.read_csv("de_haller_boundary.csv") ok = BND[np.isfinite(BND.phi_min) & np.isfinite(BND.phi_min_closed)] rows.append(("Diffusion boundary, closed form against numerical bisection", float((ok.phi_min_closed - ok.phi_min).abs().max() / ok.phi_min.max()))) o = op.objectives(int(best.n_stages), int(best.n_units), float(best.phi), float(best.U_rms), psi=float(best.psi)) d_s = max([op.point(int(best.n_stages), int(best.n_units), float(best.phi), float(best.U_rms), a, M_, 10.0, want_mass=True, psi=float(best.psi)) for a in op.ALT_M for M_ in op.MACH], key=lambda c: c["m_unit"]) # NOT an independent check: the expression below is the same one # objectives() evaluates, over the same corner set, with the same arguments, # so its residual is zero by construction. It is retained as bookkeeping and # labelled as such in the report, not credited as a measurement. rows.append(("System mass, optimiser against same-expression re-solve (identity)", abs(o["mass_total"] - d_s["m_unit"] * int(best.n_units)) / o["mass_total"])) V = pd.DataFrame(rows, columns=["test", "result"]) V.to_csv("verification.csv", index=False) return V # ---- 4. validation against the frozen parent package -------------------------------- def validation(): """The package is a 20 km, Mach 0.70 machine and is reproduced at ITS condition. The fuel-air ratio is reported but not scored: the package ran a hotter turbine inlet than the limit imposed here, and that temperature is solved for.""" d = cp.evaluate({"n_stages": 12, "n_units": 1, "U_rms": 296.147, "m_prod_total": 6.14582727 * 0.4519778}) PKG = dict(b=0.4519778, f=0.0230033, A_in=0.3807215, A_out=0.0355616) V = pd.DataFrame([(k, d[k], PKG[k], abs(d[k] - PKG[k]) / abs(PKG[k])) for k in PKG], columns=["quantity", "this_work", "package", "rel_diff"]) Tt2p = 216.65 * 1.098 w = CPC * Tt2p * (d["PR"] ** ((GC - 1) / GC) - 1) / BASE["eta_c"] Tt3p = Tt2p + w / CPC lo, hi = Tt3p + 1.0, 2500.0 for _ in range(200): mid = 0.5 * (lo + hi) dh = BASE["cp_h"] * (mid - Tt3p) fm = dh / (BASE["eta_b"] * BASE["LHV"] - dh) if fm < PKG["f"]: lo = mid else: hi = mid V["note"] = ["", "package TIT %.1f K against %.1f K here; not scored" % (0.5 * (lo + hi), BASE["TIT"]), "", ""] V.loc[V.quantity == "f", "rel_diff"] = np.nan V.to_csv("validation.csv", index=False) return V, 0.5 * (lo + hi) # ---- 5. reconciliation -------------------------------------------------------------- def reconciliation(D, BND, best): def n_needed(p_del, p_inf, M, U, psi): pt2 = p_inf * (1 + 0.2 * M ** 2) ** 3.5 * BASE["pi_inlet"] Tt2 = 216.65 * (1 + 0.2 * M ** 2) return ((p_del / pt2) ** ((GC - 1) / GC) - 1) * Tt2 * CPC / (BASE["eta_c"] * psi * U ** 2) _, p60, _ = cp.isa(60000 * FT) psi_b = float(best.psi) rows = [ dict(quantity="Stage loading below which the limit never binds", estimate=1 - DHMIN, computed=float(BND[BND.free].psi.max()), band=0.02), # At a stage loading at or below 1 - DH the limit is inactive and the # closed form has no root: there is no "lowest legal phi" to predict, and # the row is skipped rather than scored against a NaN. *([dict(quantity="Lowest legal flow coefficient at the recommended loading", estimate=phi_min_closed(psi_b), computed=float(D[abs(D.psi - psi_b) < 1e-9].phi.min()), band=0.05)] if np.isfinite(phi_min_closed(psi_b)) else []), dict(quantity="Lowest legal flow coefficient at the package loading", estimate=phi_min_closed(PSI0), computed=PHI0, band=0.01), dict(quantity="Stages for 2 bar at 60 kft, recommended loading and blade speed", estimate=n_needed(2.0e5, p60, 0.70, float(best.U_rms), psi_b), computed=float(best.n_stages), band=0.05), dict(quantity="Delivery pressure of the least-fuel design (bar)", estimate=2.0, computed=float(D.loc[D.fuel_kgh.idxmin()].p_bar), band=0.02), ] R = pd.DataFrame(rows) R["ratio"] = R.computed / R.estimate R["verdict"] = np.where(np.abs(R.ratio - 1) <= R.band, "within band", "outside band") R.to_csv("reconciliation.csv", index=False) return R if __name__ == "__main__": BND = boundary() free = BND[BND.free] print("diffusion limit inactive at any flow coefficient for psi <= %.4f" % free.psi.max()) print("closed form phi_min(psi) = sqrt((DH^2-(1-psi)^2)/(1-DH^2)); " "at the package loading %.4f -> %.6f" % (PSI0, phi_min_closed(PSI0))) D, PF, meta = search() print("\nFEASIBLE DESIGNS: %d of %d candidates (%.3f%%)" % (meta["feasible"], meta["candidates"], 100 * meta["feasible"] / meta["candidates"])) b = D.loc[D.mass_total.idxmin()] print("lightest: %d stages, %d units, phi %.3f, psi %.3f, U %.1f -> " "%.1f kg, %.1f kg/h, span %.2f mm, DH %.4f" % (b.n_stages, b.n_units, b.phi, b.psi, b.U_rms, b.mass_total, b.fuel_kgh, b.h_min_mm, b.de_haller)) bf = D.loc[D.fuel_kgh.idxmin()] print("least fuel: %d stages, %d units, phi %.3f, psi %.3f -> %.1f kg, %.1f kg/h" % (bf.n_stages, bf.n_units, bf.phi, bf.psi, bf.mass_total, bf.fuel_kgh)) print("\nPareto front (%d points):" % len(PF)) print(PF[["n_stages", "n_units", "phi", "psi", "U_rms", "de_haller", "p_bar", "mass_total", "fuel_kgh", "h_min_mm", "active_at_limit"]] .round(3).to_string(index=False)) print("\nactive constraint census:\n%s" % D.active_at_limit.value_counts().to_string()) FZ, nfz = frozen_search() print("\nFROZEN STAGE LOADING (psi = %.7f), %d nodes: %d feasible" % (PSI0, nfz, len(FZ))) if len(FZ): bz = FZ.loc[FZ.mass_total.idxmin()] print(" lightest: %d stages, %d units, phi %.3f, U %.1f -> %.1f kg, %.1f kg/h, span %.4f mm" % (bz.n_stages, bz.n_units, bz.phi, bz.U_rms, bz.mass_total, bz.fuel_kgh, bz.h_min_mm)) print(" => freeing psi buys %.0f%% of mass, NOT feasibility." % (100 * (1 - b.mass_total / bz.mass_total))) print(" feasible designs within 0.01 of the package loading: %d" % int((D.psi - PSI0).abs().le(0.01).sum())) o_rec = op.objectives(int(b.n_stages), int(b.n_units), float(b.phi), float(b.U_rms), psi=float(b.psi)) print("\nfuel at full duty is a CORNER MEAN: %.3f kg/h over corners spanning " "%.2f to %.2f kg/h (%.1f%%)" % (o_rec["fuel_kgh"], o_rec["fuel_full_duty_min"], o_rec["fuel_full_duty_max"], 100 * (o_rec["fuel_full_duty_max"] / o_rec["fuel_full_duty_min"] - 1))) ties = D[(D.fuel_kgh - D.fuel_kgh.min()).abs() < 1e-9] print("designs tied for least fuel: %d, masses %.0f to %.0f kg" % (len(ties), ties.mass_total.min(), ties.mass_total.max())) print("\nspan and mass belong to different corners:") for a_ in op.ALT_M: dd = op.point(int(b.n_stages), int(b.n_units), float(b.phi), float(b.U_rms), a_, 0.70, 10.0, want_mass=True, psi=float(b.psi)) print(" %5.0f ft: span %.3f mm, r_tip %.1f mm, rpm %.0f, m_unit %.2f kg" % (a_ / FT, dd["h_last_mm"], 1000 * dd["r_tip"], dd["rpm"], dd["m_unit"])) print("\nbleed fraction at fixed delivery pressure varies with flow coefficient:") for pr_ in (0.31, 0.45, 0.59, 0.80): dd = op.point(int(b.n_stages), int(b.n_units), pr_, float(b.U_rms), 60000 * FT, 0.70, 10.0, psi=float(b.psi)) if dd: print(" phi %.2f -> b %.6f at %.4f bar" % (pr_, dd["b"], dd["p_del_bar"])) print("\ncensored range endpoints (feasible designs sitting on a search bound):") print(" psi at %.3f: %d | phi at %.4f: %d | stages at %d: %d" % (D.psi.min(), int((D.psi == D.psi.min()).sum()), D.phi.min(), int((D.phi == D.phi.min()).sum()), D.n_stages.max(), int((D.n_stages == D.n_stages.max()).sum()))) print("\nverification residuals:\n%s" % verification(b).to_string(index=False)) VD, TIT_pkg = validation() print("\nvalidation against the frozen package:\n%s" % VD.to_string(index=False)) print("\npredictions against the computed space:\n%s" % reconciliation(D, BND, b).round(4).to_string(index=False))