From ccbcdc6ed480235f5000e65899ff25a94a116419 Mon Sep 17 00:00:00 2001 From: James Hansen <109622320+jjhans@users.noreply.github.com> Date: Fri, 3 Jul 2026 10:59:42 -0700 Subject: [PATCH] [#97] Introduce throttle control; tabulate JICF over momentum-flux ratio Replace the pre-specified mass flow rate schedule with a throttle-driven, dynamically-updatable injector, and tabulate the JICF model over a range of momentum-flux ratios J (at a user-specified resolution) rather than a schedule of mass flow rates. The original fixed-schedule behavior is recovered by setting the throttle as a function of time. - JICModel: throttle-agnostic library keyed on a J grid instead of a time schedule. Takes J + per-J (rho_inj, u_inj, T_inj); sorts J ascending and reorders the injected-fluid state to match. All tables and interpolators are keyed on J. J is stored in the HDF5 groups and guarded on load (_cached_J_matches) so a stale mdot-indexed cache is not silently reused. - FuelInjector: now holds the runtime throttle schedule and derives J(t) from the live crossflow inflow and jet state, keeping lookups correct if the inflow is later allowed to vary. Convected fluid tips carry [x, mdot, J]. - Combustor: injector= now takes pre-built FuelInjector(s). - get_injectors_fpv: builds the J grid by sweeping the equivalence-ratio range at the requested resolution, plus a separate runtime throttle schedule, and returns a FuelInjector. plot.py, hyshot_ii_jic.py, and example 04 updated to the new interfaces. --- .../jicf_heat_release_profile.py | 56 ++++++-- examples/hyshot_ii/case_setup.py | 77 ++++++++--- examples/hyshot_ii/hyshot_ii_jic.py | 2 +- src/stanshock/components/combustor.py | 15 +- src/stanshock/models/jicf/generate.py | 129 +++++++++--------- src/stanshock/models/jicf/source.py | 91 ++++++++++-- src/stanshock/processing/plot.py | 4 +- 7 files changed, 253 insertions(+), 121 deletions(-) diff --git a/examples/04_jicf_heat_release_profile/jicf_heat_release_profile.py b/examples/04_jicf_heat_release_profile/jicf_heat_release_profile.py index 321443e..4d0336e 100644 --- a/examples/04_jicf_heat_release_profile/jicf_heat_release_profile.py +++ b/examples/04_jicf_heat_release_profile/jicf_heat_release_profile.py @@ -10,7 +10,7 @@ import numpy as np from stanshock.components.combustor import Combustor -from stanshock.models.jicf import JICModel +from stanshock.models.jicf import FuelInjector, JICModel from stanshock.numerics.boundary_conditions import BCInput, SpecifiedFace from stanshock.physics.flamelet import FPVTable from stanshock.processing.csv_writer import CSVWriter @@ -152,21 +152,28 @@ def build_fpv_table(paths: dict[str, Path]) -> FPVTable: def build_injector( geometry: Box, physics: FPVTable, paths: dict[str, Path], t_end: float -) -> JICModel: - t_inj = np.array([0.0, t_end]) - rho_inj = np.array([RHO_F, RHO_F]) - +) -> FuelInjector: cache_dir = paths["cache_dir"] - return JICModel( - fuel="H2", + + # Momentum-flux ratio at the (single) operating condition, and a small + # grid bracketing it for the throttle-agnostic table. Under choked + # injection the jet velocity/temperature are ~constant, so only the + # injected density varies across the grid. + J0 = RHO_F * U_F**2 / (RHO_IN * U_IN**2) + J_grid = np.linspace(0.5 * J0, 1.5 * J0, 5) + rho_grid = J_grid * RHO_IN * U_IN**2 / U_F**2 + u_grid = np.full_like(J_grid, U_F) + T_grid = np.full_like(J_grid, T_F) + + jicf = JICModel( x_inj=X_INJ, x_noz=L_DOMAIN, n_inj=N_INJ, d_inj=D_F, - t_inj=t_inj, - rho_inj=rho_inj, - u_inj=U_F, - T_inj=T_F, + J=J_grid, + rho_inj=rho_grid, + u_inj=u_grid, + T_inj=T_grid, rho=RHO_IN, u=U_IN, T=T_IN, @@ -176,6 +183,33 @@ def build_injector( physics=physics, ) + # Constant-throttle runtime schedule at the operating condition. + A_f = np.pi * (D_F / 2.0) ** 2 + mdot_f = N_INJ * RHO_F * U_F * A_f + A_in = float(geometry.area(0.0, geometry.xf[0])) + mdot_ox = RHO_IN * U_IN * A_in + gas = physics.gas + gas.TDX = T_IN, RHO_IN, X_OX + Y_ox = gas.Y.copy() + gas.TDX = T_F, RHO_F, X_F + Y_f = gas.Y.copy() + E_f = float(gas.int_energy_mass) + 0.5 * U_F**2 + w_f = mdot_f / (mdot_f + mdot_ox) + gas.Y = (1.0 - w_f) * Y_ox + w_f * Y_f + phi0 = float(gas.equivalence_ratio(X_F, X_OX)) + + t_inj = np.array([0.0, t_end]) + return FuelInjector( + jicf, + t_inj=t_inj, + phi_inj=np.array([phi0, phi0]), + mdot_inj=np.array([mdot_f, mdot_f]), + u_inj=np.array([U_F, U_F]), + E_inj=np.array([E_f, E_f]), + geometry=geometry, + physics=physics, + ) + def collect_profiles(combustor: Combustor) -> dict[str, np.ndarray]: idx = combustor.geometry.idx_cells diff --git a/examples/hyshot_ii/case_setup.py b/examples/hyshot_ii/case_setup.py index 7402e08..fe0731c 100644 --- a/examples/hyshot_ii/case_setup.py +++ b/examples/hyshot_ii/case_setup.py @@ -9,7 +9,7 @@ from stanshock.components.combustor import Combustor from stanshock.models.inlet_diffuser import InletDiffuser -from stanshock.models.jicf import JICModel +from stanshock.models.jicf import FuelInjector, JICModel from stanshock.models.wall_models import ( CompressibleHeatFlux, CompressibleReactingSkinFriction, @@ -323,7 +323,9 @@ def get_injectors_fpv( U_in: float, mdot: Literal["constant", "schedule"] = "constant", fpv_dir: Path = Path("./data"), -) -> JICModel: + phi_range: tuple[float, float] = (0.1, 0.7), + n_throttle: int = 8, +) -> FuelInjector: assert isinstance(geometry, AsymmetricBox) # Freeze after constant cross section region of the combustor L_const = 300.0e-3 # m @@ -390,29 +392,42 @@ def get_injectors_fpv( ] ) - t_f = np.zeros(t_phi_gl_schedule.shape[0]) - rho_f = np.zeros(t_phi_gl_schedule.shape[0]) - U_f = np.zeros(t_phi_gl_schedule.shape[0]) - T_f = np.zeros(t_phi_gl_schedule.shape[0]) - for i in range(t_phi_gl_schedule.shape[0]): - t_f[i] = t_phi_gl_schedule[i, 0] - rho_f[i], U_f[i], T_f[i] = fuel_props_from_phi( - physics, t_phi_gl_schedule[i, 1], mdot_ox, T0_f, P_in, A_f_tot - ) - # # NOTE: Assuming perfect gas & isentropic choked flow, only rho_f changes with phi/mdot - # U_f = U_f[-1] - # T_f = T_f[-1] + def _props(phi_gl: float) -> tuple[float, float, float, float, float, float]: + """Injected-fluid state for a global equivalence ratio ``phi_gl``. - return JICModel( + Returns ``(rho_f, U_f, T_f, mdot_f, E_f, J)``, where ``J`` is the + momentum-flux ratio relative to the (fixed) crossflow inflow. + """ + rho_f, U_f, T_f = fuel_props_from_phi( + physics, phi_gl, mdot_ox, T0_f, P_in, A_f_tot + ) + mdot_f = N_f * rho_f * U_f * A_f + J = rho_f * U_f**2 / (rho_in * U_in**2) + gas = physics.gas + gas.TDX = T_f, rho_f, physics.fuel_def + E_f = float(gas.int_energy_mass) + 0.5 * U_f**2 + return rho_f, U_f, T_f, mdot_f, E_f, J + + # Tabulation grid: discretize the throttle range as a range of momentum-flux + # ratios J (the throttle-agnostic table dimension) at the requested + # resolution. The runtime throttle schedule below maps onto this table. + phi_grid = np.linspace(phi_range[0], phi_range[1], n_throttle) + rho_g = np.zeros(n_throttle) + U_g = np.zeros(n_throttle) + T_g = np.zeros(n_throttle) + J_g = np.zeros(n_throttle) + for i, phi_gl in enumerate(phi_grid): + rho_g[i], U_g[i], T_g[i], _, _, J_g[i] = _props(phi_gl) + + jicf = JICModel( x_inj=x_inj, x_noz=L_const, n_inj=N_f, d_inj=2 * r_f, - t_inj=t_f, - phi_inj=t_phi_gl_schedule[:, 1], - rho_inj=rho_f, - u_inj=U_f, - T_inj=T_f, + J=J_g, + rho_inj=rho_g, + u_inj=U_g, + T_inj=T_g, rho=rho_in, u=U_in, T=T_in, @@ -422,6 +437,28 @@ def get_injectors_fpv( datadir=fpv_dir, ) + # Runtime throttle schedule: recover the original fixed-schedule behaviour by + # evaluating the injected-fluid state at each scheduled equivalence ratio. + t_f = t_phi_gl_schedule[:, 0] + phi_s = t_phi_gl_schedule[:, 1] + n_s = t_phi_gl_schedule.shape[0] + mdot_s = np.zeros(n_s) + U_s = np.zeros(n_s) + E_s = np.zeros(n_s) + for i in range(n_s): + _, U_s[i], _, mdot_s[i], E_s[i], _ = _props(phi_s[i]) + + return FuelInjector( + jicf, + t_inj=t_f, + phi_inj=phi_s, + mdot_inj=mdot_s, + u_inj=U_s, + E_inj=E_s, + geometry=geometry, + physics=physics, + ) + class Hyshot2Interface: """Wrapper around HyShot-II scramjet to facilitate calling from 6-DOF trajectory simulations.""" diff --git a/examples/hyshot_ii/hyshot_ii_jic.py b/examples/hyshot_ii/hyshot_ii_jic.py index 92094da..81c6793 100644 --- a/examples/hyshot_ii/hyshot_ii_jic.py +++ b/examples/hyshot_ii/hyshot_ii_jic.py @@ -55,7 +55,7 @@ # ["Y_H2", "Y_OH", "Y_H2O"], # ] ss.xt_diagrams = [XTDiagram(ss, variable, skip_steps=10) for variable in plot_variables] -ss.advance_simulation(ss.injectors[0].jicf.t_inj[-1]) +ss.advance_simulation(ss.injectors[0].t_inj[-1]) for diagram in ss.xt_diagrams: diagram.plot(figdir=fig_dir) diff --git a/src/stanshock/components/combustor.py b/src/stanshock/components/combustor.py index b11738f..251adc0 100644 --- a/src/stanshock/components/combustor.py +++ b/src/stanshock/components/combustor.py @@ -4,7 +4,7 @@ from stanshock.models.area_change import AreaChange from stanshock.models.boundary_layer import BoundaryLayer -from stanshock.models.jicf import FuelInjector, JICFChemistrySource, JICModel +from stanshock.models.jicf import FuelInjector, JICFChemistrySource from stanshock.models.pseudoshock import Pseudoshock from stanshock.models.wall_models import HeatFlux, SkinFriction from stanshock.numerics.boundary_conditions import ( @@ -69,7 +69,7 @@ def __init__( source_terms: RightHandSide | list[RightHandSide] | None = None, # Catch-all source term(s) - injector: list[JICModel] | JICModel | None = None, # injector model + injector: list[FuelInjector] | FuelInjector | None = None, # injector model(s) flux_function: RiemannSolver = hllc_flux_vectorized, inviscid_face_extrapolator: type[FaceExtrapolator] = FifthOrderWeno, viscous_face_extrapolator: type[FaceExtrapolator] = FirstOrder, @@ -98,16 +98,9 @@ def __init__( if injector is None: self.injectors: list[FuelInjector] = [] elif isinstance(injector, list): - self.injectors = [ - FuelInjector(inj, geometry=inj.geometry, physics=inj.physics) - for inj in injector - ] + self.injectors = injector else: - self.injectors = [ - FuelInjector( - injector, geometry=injector.geometry, physics=injector.physics - ) - ] + self.injectors = [injector] self.optimization_iteration = optimization_iteration self.physics: FluidPhysics = physics self.initialization: Initialization = initialization diff --git a/src/stanshock/models/jicf/generate.py b/src/stanshock/models/jicf/generate.py index 3bc9aee..ae1ac3c 100644 --- a/src/stanshock/models/jicf/generate.py +++ b/src/stanshock/models/jicf/generate.py @@ -29,8 +29,7 @@ def __init__( x_noz: float, n_inj: int, d_inj: float, - t_inj: Array, - phi_inj: Array, + J: Array, rho_inj: Array, u_inj: Array, T_inj: Array, @@ -47,6 +46,13 @@ def __init__( """ This method initializes the Jet-in-Crossflow model with the following parameters: + + The model is tabulated over a range of momentum-flux ratios ``J`` + (the throttle-agnostic table dimension) rather than a specific + schedule of mass flow rates. The runtime throttle schedule lives in the + :class:`~stanshock.models.jicf.source.FuelInjector`, which maps the + live throttle/inflow state onto this table via ``J``. + x_inj: float The x-coordinate of the injection point x_noz: float @@ -57,16 +63,16 @@ def __init__( The diameter of the injected jet theta_inj: float The angle of the jet relative to the x axis (rads) - t_inj: np.ndarray - Time array for the injection profile - phi_inj: np.ndarray - Scheduled equivalence ratio of injected jet + J: np.ndarray + The grid of momentum-flux ratios to tabulate over. Sets the table + resolution and range; each entry corresponds to one injected-fluid + state ``(rho_inj[i], u_inj[i], T_inj[i])``. rho_inj: np.ndarray - The density of the injected jet, as a function of time + The density of the injected jet at each ``J`` grid point u_inj: np.ndarray - The velocity of the injected jet + The velocity of the injected jet at each ``J`` grid point T_inj: np.ndarray - The temperature of the injected jet + The temperature of the injected jet at each ``J`` grid point rho: float The density of the crossflow u: float @@ -106,10 +112,6 @@ def __init__( self.d_inj = d_inj self.theta_inj = theta_inj if theta_inj is not None else 0.0 - self.rho_inj = rho_inj - self.u_inj = u_inj - self.T_inj = T_inj - self.rho = rho self.alpha = alpha if alpha is not None else 1e6 @@ -136,13 +138,19 @@ def __init__( self.M = self.u / self.c self.Y_ox = gas.Y - # Properties of the injected fluid - self.t_inj = t_inj - self.phi_inj = phi_inj - sol = ct.SolutionArray(gas, shape=self.t_inj.shape) + # Momentum-flux-ratio grid (the table dimension). Sort ascending and + # carry the per-grid injected-fluid state along so the tables and the + # RegularGridInterpolators built below share a monotone axis. + order = np.argsort(np.asarray(J, dtype=float)) + self.J = np.asarray(J, dtype=float)[order] + self.rho_inj = np.asarray(rho_inj, dtype=float)[order] + self.u_inj = np.asarray(u_inj, dtype=float)[order] + self.T_inj = np.asarray(T_inj, dtype=float)[order] + + # Properties of the injected fluid at each grid point + sol = ct.SolutionArray(gas, shape=self.J.shape) sol.TDX = self.T_inj, self.rho_inj, self.fuel_def self.p_inj = sol.P - self.E_inj = sol.int_energy_mass + 0.5 * self.u_inj**2 self.W_inj = sol.mean_molecular_weight self.gamma_inj = sol.cp / sol.cv self.c_inj = sol.sound_speed @@ -150,25 +158,6 @@ def __init__( self.M_inj = 1.0 self.mdot_inj = self.n_inj * self.rho_inj * self.u_inj * self.A_inj self.mdot_inj[np.isnan(self.mdot_inj)] = 0.0 - self.mdot_inj_unique, self.mdot_inj_unique_idx = np.unique( - self.mdot_inj, return_index=True - ) - self.mdot_inj_unique_idx = self.mdot_inj_unique_idx[ - np.argsort(self.mdot_inj_unique) - ] - - self.mdot_inj_unique = self.mdot_inj[self.mdot_inj_unique_idx] - self.rho_inj_unique = self.rho_inj[self.mdot_inj_unique_idx] - self.u_inj_unique = self.u_inj[self.mdot_inj_unique_idx] - self.p_inj_unique = self.p_inj[self.mdot_inj_unique_idx] - - # Mass flow rate and equivalence ratio schedules - self.mdot_f_interp = interpolate.interp1d( - self.t_inj, self.mdot_inj, bounds_error=False, fill_value=0.0 - ) - self.phi_f_interp = interpolate.interp1d( - self.t_inj, self.phi_inj, bounds_error=False, fill_value=0.0 - ) # Set up the analytic JICF model self.analytic = AnalyticJICF( @@ -177,16 +166,18 @@ def __init__( h=self.h, n_inj=n_inj, d_inj=d_inj, - rho_inj=self.rho_inj_unique, - u_inj=self.u_inj_unique, + rho_inj=self.rho_inj, + u_inj=self.u_inj, rho=rho, u=u, physics=self.physics, theta_inj=theta_inj, ) - # Precompute a 3D array of the mixture fraction and generate an interpolator - if self._h5_has("Z_3D/Z"): + # Precompute a 3D array of the mixture fraction and generate an interpolator. + # Only reuse a cached table if it was generated on the same J grid; + # otherwise it belongs to a different throttle range and is recomputed. + if self._h5_has("Z_3D/Z") and self._cached_J_matches("Z_3D"): with h5py.File(self.model_file, "r") as f: group = f["Z_3D"] self.x_3D_data = group["x"][:] @@ -199,7 +190,7 @@ def __init__( # Precompute the axial mean and variance profiles of Z, along with the # mesh they were generated on. - if self._h5_has("Z_profiles/Z_avg"): + if self._h5_has("Z_profiles/Z_avg") and self._cached_J_matches("Z_profiles"): with h5py.File(self.model_file, "r") as f: group = f["Z_profiles"] self.x_profile = group["x"][:] @@ -208,19 +199,19 @@ def __init__( else: self.calc_Z_avg_var_profiles(write=True) - # Map mdot -> Z mean/variance on the *generation* mesh (self.x_profile). + # Map J -> Z mean/variance on the *generation* mesh (self.x_profile). # Building the interpolators on the stored mesh rather than the current # simulation mesh decouples the profiles from the simulation grid, so the # model does not need to be regenerated when the mesh changes. Out-of-range # query points (e.g. ghost cells) are linearly extrapolated. self.Z_avg_profile_interp = interpolate.RegularGridInterpolator( - (self.mdot_inj_unique, self.x_profile), + (self.J, self.x_profile), self.Z_avg_profile, bounds_error=False, fill_value=None, ) self.Z_var_profile_interp = interpolate.RegularGridInterpolator( - (self.mdot_inj_unique, self.x_profile), + (self.J, self.x_profile), self.Z_var_profile, bounds_error=False, fill_value=None, @@ -247,6 +238,20 @@ def _h5_has(self, key: str) -> bool: with h5py.File(self.model_file, "r") as f: return key in f + def _cached_J_matches(self, group: str) -> bool: + """Return True if the cached group was generated on the current J grid. + + Guards against silently reusing a table tabulated over a different + throttle range (momentum-flux-ratio grid). + """ + if not self.model_file.exists(): + return False + with h5py.File(self.model_file, "r") as f: + if group not in f or "J" not in f[group]: + return False + J_cached = f[group]["J"][:] + return J_cached.shape == self.J.shape and np.allclose(J_cached, self.J) + def _h5_write(self, group: str, data: dict[str, Array]) -> None: """Write (overwriting if present) a group of named arrays to the model file.""" with h5py.File(self.model_file, "a") as f: @@ -277,7 +282,7 @@ def calc_Z_3D_interp(self, write=False): self.xc[0], self.xc[-1], dx, 1.1, self.x_inj ) Nx = len(self.x_3D_data) - self.Z_3D_data = np.zeros([len(self.mdot_inj_unique), Nx, Ny, Nz]) + self.Z_3D_data = np.zeros([len(self.J), Nx, Ny, Nz]) for i in tqdm(range(Nx)): for j in range(Ny): for k in range(Nz): @@ -286,12 +291,13 @@ def calc_Z_3D_interp(self, write=False): self.y_3D_data[j], self.z_3D_data[k], ) - self.Z_3D_data[np.isnan(self.rho_inj_unique)] = 0.0 + self.Z_3D_data[np.isnan(self.rho_inj)] = 0.0 if write: self._h5_write( "Z_3D", { + "J": self.J, "x": self.x_3D_data, "y": self.y_3D_data, "z": self.z_3D_data, @@ -303,7 +309,7 @@ def calc_Z_3D_interp(self, write=False): def _build_Z_3D_interp(self) -> None: self.Z_3D_interp = [] - for i_m in range(len(self.mdot_inj_unique)): + for i_m in range(len(self.J)): interp = interpolate.RegularGridInterpolator( (self.x_3D_data, self.y_3D_data, self.z_3D_data), self.Z_3D_data[i_m], @@ -312,17 +318,17 @@ def _build_Z_3D_interp(self) -> None: self.Z_3D_interp.append(interp) def eval_Z_3D_interp(self, x, y, z): - Z_arr = np.zeros_like(self.mdot_inj_unique) - for i_m in range(len(self.mdot_inj_unique)): + Z_arr = np.zeros_like(self.J) + for i_m in range(len(self.J)): Z_arr[i_m] = self.Z_3D_interp[i_m]((x, y, z)) return Z_arr def Z_avg_var(self, x): - Z_avg = np.zeros_like(self.mdot_inj_unique) - Z_var = np.zeros_like(self.mdot_inj_unique) + Z_avg = np.zeros_like(self.J) + Z_var = np.zeros_like(self.J) x_local = x - self.x_inj - for i_m in range(len(self.mdot_inj_unique)): - if np.isnan(self.rho_inj_unique[i_m]): + for i_m in range(len(self.J)): + if np.isnan(self.rho_inj[i_m]): Z_avg[i_m] = 0.0 Z_var[i_m] = 0.0 continue @@ -351,10 +357,10 @@ def func(z, y, i_m=i_m): return Z_avg, Z_var def Z_avg_var_adjusted(self, x): - Z_avg = np.zeros_like(self.mdot_inj_unique) - Z_var = np.zeros_like(self.mdot_inj_unique) - for i_m in range(len(self.mdot_inj_unique)): - if np.isnan(self.rho_inj_unique[i_m]): + Z_avg = np.zeros_like(self.J) + Z_var = np.zeros_like(self.J) + for i_m in range(len(self.J)): + if np.isnan(self.rho_inj[i_m]): Z_avg[i_m] = 0.0 Z_var[i_m] = 0.0 continue @@ -389,9 +395,9 @@ def calc_Z_avg_var_profiles(self, write=False): # Record the mesh the profiles are generated on so the interpolators can # be rebuilt independently of the simulation mesh. self.x_profile = np.asarray(self.xc, dtype=float) - n_mdot = len(self.mdot_inj_unique) - self.Z_avg_profile = np.zeros([n_mdot, len(self.x_profile)]) - self.Z_var_profile = np.zeros([n_mdot, len(self.x_profile)]) + n_J = len(self.J) + self.Z_avg_profile = np.zeros([n_J, len(self.x_profile)]) + self.Z_var_profile = np.zeros([n_J, len(self.x_profile)]) for i in tqdm(range(len(self.x_profile))): if self.x_profile[i] > self.x_noz: # Freeze the profiles in the nozzle @@ -406,6 +412,7 @@ def calc_Z_avg_var_profiles(self, write=False): self._h5_write( "Z_profiles", { + "J": self.J, "x": self.x_profile, "Z_avg": self.Z_avg_profile, "Z_var": self.Z_var_profile, diff --git a/src/stanshock/models/jicf/source.py b/src/stanshock/models/jicf/source.py index aa34d61..93ebd39 100644 --- a/src/stanshock/models/jicf/source.py +++ b/src/stanshock/models/jicf/source.py @@ -1,6 +1,7 @@ from __future__ import annotations import numpy as np +from scipy import interpolate from stanshock.models.jicf.generate import JICModel from stanshock.physics.flamelet import FPVTable @@ -11,14 +12,40 @@ class FuelInjector(RightHandSide): - """Source terms for a fuel injector based on a JICF model.""" + """Source terms for a fuel injector based on a JICF model. + + The injected mass flow rate is driven by a runtime *throttle schedule* + (``t_inj`` -> ``mdot_inj``/``u_inj``/``E_inj``), which is decoupled from the + throttle-agnostic :class:`JICModel` tables. The original fixed-schedule + behaviour is recovered by passing that schedule here. Table lookups are + performed on the momentum-flux ratio ``J``, computed from the injected jet + state and the crossflow inflow so that they remain correct if the inflow is + allowed to vary. + """ def __init__( - self, jicf: JICModel, **precompute_steps: Unpack[PrecomputeSteps] + self, + jicf: JICModel, + t_inj: Array, + phi_inj: Array, + mdot_inj: Array, + u_inj: Array, + E_inj: Array, + **precompute_steps: Unpack[PrecomputeSteps], ) -> None: """ jicf: JICModel - The underlying jet-in-crossflow model + The underlying jet-in-crossflow model (tabulated over ``J``) + t_inj: Array + Time stamps of the throttle schedule + phi_inj: Array + Scheduled global equivalence ratio (used for reporting/plotting) + mdot_inj: Array + Scheduled total injected mass flow rate at each time stamp + u_inj: Array + Scheduled injected jet velocity at each time stamp + E_inj: Array + Scheduled injected total energy per unit mass at each time stamp geometry: Box The geometry object describing the mesh and cross-section physics: FPVTable @@ -29,14 +56,46 @@ def __init__( assert isinstance(self.physics, FPVTable) self.jicf = jicf + # Runtime throttle schedule + self.t_inj = np.asarray(t_inj, dtype=float) + self.phi_inj = np.asarray(phi_inj, dtype=float) + self.mdot_inj = np.asarray(mdot_inj, dtype=float) + self.u_inj = np.asarray(u_inj, dtype=float) + self.E_inj = np.asarray(E_inj, dtype=float) + + # Interpolants over the throttle schedule (also used for reporting/plotting) + self.mdot_f_interp = interpolate.interp1d( + self.t_inj, self.mdot_inj, bounds_error=False, fill_value=0.0 + ) + self.phi_f_interp = interpolate.interp1d( + self.t_inj, self.phi_inj, bounds_error=False, fill_value=0.0 + ) + # Extract some information about the geometry self.xc = self.geometry.xc[self.idx_input] self.cell_volumes: Array | None = None if self.geometry.dlnA_dt is None: self.cell_volumes = self.geometry.volume()[:, None] - # Position of the first injected fluid particle - self.fluid_tips = np.array([[self.jicf.x_inj, self.jicf.mdot_inj[0]]]) + # Position of the first injected fluid particle: [x, mdot, J]. Both the + # mass flow rate (for reporting) and the momentum-flux ratio (the table + # lookup axis) are convected with each parcel. + self.fluid_tips = np.array( + [[self.jicf.x_inj, self.mdot_inj[0], self._J(self.t_inj[0])]] + ) + + def _J(self, t: float) -> float: + """Momentum-flux ratio ``J = rho_inj u_inj**2 / (rho u**2)`` at time ``t``. + + The crossflow inflow (``rho``, ``u``) is currently fixed on the model; + when it becomes a live flow-field state this is where it would be read. + """ + mdot = float(np.interp(t, self.t_inj, self.mdot_inj)) + u_inj = float(np.interp(t, self.t_inj, self.u_inj)) + if mdot <= 0.0 or u_inj <= 0.0: + return 0.0 + rho_inj = mdot / (self.jicf.n_inj * u_inj * self.jicf.A_inj) + return rho_inj * u_inj**2 / (self.jicf.rho * self.jicf.u**2) def update_fluid_tip_positions(self, dt, t, u): """ @@ -54,9 +113,9 @@ def update_fluid_tip_positions(self, dt, t, u): # Update the fluid tip positions self.fluid_tips[:, 0] += dt * np.interp(self.fluid_tips[:, 0], self.xc, u) - # Emit a new fluid tip - mdot = np.interp(t, self.jicf.t_inj, self.jicf.mdot_inj) - next_tip = np.array([[self.jicf.x_inj, mdot]]) + # Emit a new fluid tip, carrying its mass flow rate and momentum-flux ratio + mdot = np.interp(t, self.t_inj, self.mdot_inj) + next_tip = np.array([[self.jicf.x_inj, mdot, self._J(t)]]) self.fluid_tips = np.concatenate([self.fluid_tips, next_tip], axis=0) # Drop fluid tips that have passed the end of the domain @@ -75,9 +134,9 @@ def source( state_array_local = np.reshape(state_array_local, self.shape_input) rhs = np.zeros_like(state_array_local) - mdot = np.interp(time, self.jicf.t_inj, self.jicf.mdot_inj) - u_inj = np.interp(time, self.jicf.t_inj, self.jicf.u_inj) - E_inj = np.interp(time, self.jicf.t_inj, self.jicf.E_inj) + mdot = np.interp(time, self.t_inj, self.mdot_inj) + u_inj = np.interp(time, self.t_inj, self.u_inj) + E_inj = np.interp(time, self.t_inj, self.E_inj) L_src = 3e-2 xf = self.geometry.xf idx = ( @@ -138,13 +197,15 @@ def source_implementation( assert state.normalized_progress_variable is not None factor = self.physics.get_source_progress_variable_compressibility_factor(state) - # Get the mixture fraction variance profile - mdot_inj = np.interp( + # Get the mixture fraction variance profile. The variance table is + # indexed by momentum-flux ratio J (column 2 of the fluid tips), which + # is convected with each injected parcel. + J_inj = np.interp( self.injector.xc, np.flip(self.injector.fluid_tips, axis=0)[:, 0], - np.flip(self.injector.fluid_tips, axis=0)[:, 1], + np.flip(self.injector.fluid_tips, axis=0)[:, 2], ) - Zvar = self.injector.jicf.Z_var_profile_interp((mdot_inj, self.injector.xc)) + Zvar = self.injector.jicf.Z_var_profile_interp((J_inj, self.injector.xc)) Zvar = np.maximum(Zvar, 10 ** self.injector.jicf.logsigma2_vec.min()) return ( diff --git a/src/stanshock/processing/plot.py b/src/stanshock/processing/plot.py index ebf3c84..d6d27b5 100644 --- a/src/stanshock/processing/plot.py +++ b/src/stanshock/processing/plot.py @@ -232,7 +232,7 @@ def update(self, domain: Combustor) -> None: if domain.injectors is not None: mdot_f = 0.0 for inj in domain.injectors: - mdot_f += float(inj.jicf.mdot_f_interp(domain.t)) + mdot_f += float(inj.mdot_f_interp(domain.t)) self.mdot.append(mdot_f) def plot(self, figdir: Path | str = ".") -> None: @@ -489,7 +489,7 @@ def plot_state( inj.fluid_tips[:, 1] * 1e3 * inj.jicf.n_inj, s=1, ) - phi_tot += inj.jicf.phi_f_interp(domain.t) + phi_tot += inj.phi_f_interp(domain.t) ax.set_title(rf"$\phi={phi_tot:.2f}$") ax.set_ymargin(0.1) ax.set_ylabel(r"$\dot{m}_f$ [g/s]")