Skip to content

[#97] Introduce throttle control; tabulate JICF over momentum-flux ratio #108

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
jjhans wants to merge 1 commit into
base: jicf-cleanup-deadcode-and-mesh
Choose a base branch
from
Open

[#97] Introduce throttle control; tabulate JICF over momentum-flux ratio #108

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 45 additions & 11 deletions examples/04_jicf_heat_release_profile/jicf_heat_release_profile.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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
Expand Down
77 changes: 57 additions & 20 deletions examples/hyshot_ii/case_setup.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand All @@ -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."""
Expand Down
2 changes: 1 addition & 1 deletion examples/hyshot_ii/hyshot_ii_jic.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down
15 changes: 4 additions & 11 deletions src/stanshock/components/combustor.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 (
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -98,16 +98,9 @@ def __init__(
if injector is None:
self.injectors: list[FuelInjector] = []
elif isinstance(injector, list):
self.injectors = [

@TimothyEDawson TimothyEDawson Jul 7, 2026

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not strictly opposed to this change, but the reasoning behind creating the FuelInjector objects internally is to ensure they use the same FluidPhysics and Geometry objects as the Combustor.

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
Expand Down
Loading
Loading