diff --git a/.gitignore b/.gitignore index e103045b..c2614ba4 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,7 @@ pycycle.egg-info/* /build* *testflo_report.out +om_pycycle.egg-info/* +inputs.html +n2.html +.openmdao_out \ No newline at end of file diff --git a/example_cycles/rde_ramjet.py b/example_cycles/rde_ramjet.py new file mode 100644 index 00000000..4c45fe63 --- /dev/null +++ b/example_cycles/rde_ramjet.py @@ -0,0 +1,253 @@ +""" +Step 5: Air-Breathing Rotating Detonation Engine (RDE) Ramjet Cycle. + +Simulates and compares a supersonic air-breathing ramjet at Mach 2.5, 40,000 ft: +1. Conventional Isobaric Ramjet (Brayton cycle with 5% combustor pressure loss). +2. Rotating Detonation Engine Ramjet (Humphrey/CJ cycle with pressure gain). + +Architecture: +[FlightConditions] -> [Inlet] -> [RDECombustor / Combustor] -> [Nozzle] -> [Performance] +""" + +import sys +import numpy as np +import openmdao.api as om +import pycycle.api as pyc + +from pycycle.elements.rde_combustor import RDECombustor, print_rde + + +class RamjetCycle(pyc.Cycle): + """ + Air-breathing ramjet propulsion cycle supporting either conventional or RDE combustor. + """ + + def initialize(self): + self.options.declare('use_rde', default=True, desc='Use RDE combustor if True, conventional if False') + self.options.declare('fuel_type', default='FAR', desc='Fuel type for thermodynamics') + super().initialize() + + def setup(self): + use_rde = self.options['use_rde'] + fuel_type = self.options['fuel_type'] + + # Use fast TABULAR thermo for air + Jet-A + self.options['thermo_method'] = 'TABULAR' + self.options['thermo_data'] = pyc.AIR_JETA_TAB_SPEC + + design = self.options['design'] + + # 1. Add cycle components + self.add_subsystem('fc', pyc.FlightConditions()) + self.add_subsystem('inlet', pyc.Inlet()) + + if use_rde: + self.add_subsystem('burner', RDECombustor(fuel_type=fuel_type, dia_annulus=14.0)) + else: + self.add_subsystem('burner', pyc.Combustor(fuel_type=fuel_type)) + + self.add_subsystem('nozz', pyc.Nozzle(nozzType='CD', lossCoef='Cv')) + self.add_subsystem('perf', pyc.Performance(num_nozzles=1, num_burners=1)) + + # 2. Connect flow stations + self.pyc_connect_flow('fc.Fl_O', 'inlet.Fl_I', connect_w=False) + self.pyc_connect_flow('inlet.Fl_O', 'burner.Fl_I') + self.pyc_connect_flow('burner.Fl_O', 'nozz.Fl_I') + + # 3. Connect ambient pressure to nozzle backpressure + self.connect('fc.Fl_O:stat:P', 'nozz.Ps_exhaust') + + # 4. Connect performance metrics + self.connect('inlet.Fl_O:tot:P', 'perf.Pt2') + self.connect('burner.Fl_O:tot:P', 'perf.Pt3') + self.connect('burner.Wfuel', 'perf.Wfuel_0') + self.connect('inlet.F_ram', 'perf.ram_drag') + self.connect('nozz.Fg', 'perf.Fg_0') + + # Resolve unit ambiguity on inlet flow rate + self.set_input_defaults('inlet.Fl_I:stat:W', 100.0, units='lbm/s') + self.set_input_defaults('burner.Fl_I:FAR', 0.028) + self.set_input_defaults('inlet.MN', 0.35) + self.set_input_defaults('burner.MN', 0.30) + + # 5. Solver configuration + newton = self.nonlinear_solver = om.NewtonSolver() + newton.options['atol'] = 1e-6 + newton.options['rtol'] = 1e-6 + newton.options['iprint'] = -1 + newton.options['maxiter'] = 30 + newton.options['solve_subsystems'] = True + newton.options['reraise_child_analysiserror'] = False + + self.linear_solver = om.DirectSolver() + + super().setup() + + +def run_ramjet_case(use_rde=True, far_val=0.028, alt_ft=40000.0, mach=2.5, w_air=100.0): + """ + Sets up and executes a ramjet simulation at the given flight condition. + """ + prob = om.Problem() + prob.model = RamjetCycle(use_rde=use_rde, fuel_type='FAR') + + prob.set_solver_print(level=-1) + prob.setup(check=False) + + # Flight conditions: Mach 2.5 at 40,000 ft + prob.set_val('fc.alt', alt_ft, units='ft') + prob.set_val('fc.MN', mach) + + # Inlet conditions (supersonic diffusion to combustor face) + prob.set_val('inlet.Fl_I:stat:W', w_air, units='lbm/s') + prob.set_val('inlet.MN', 0.35) + prob.set_val('inlet.ram_recovery', 0.90) + + # Combustor parameters + prob.set_val('burner.Fl_I:FAR', far_val) + prob.set_val('burner.MN', 0.30) + + if use_rde: + prob.set_val('burner.dPqP_inj', 0.12) + prob.set_val('burner.eta_rde', 0.85) + prob.set_val('burner.dia_annulus', 14.0, units='inch') + prob.set_val('burner.N_waves', 1.0) + else: + prob.set_val('burner.dPqP', 0.05) # 5% conventional pressure loss + + # Nozzle velocity coefficient + prob.set_val('nozz.Cv', 0.98) + + prob.run_model() + + # Extract metrics + p0 = prob.get_val('fc.Fl_O:stat:P', units='psi')[0] + pt0 = prob.get_val('fc.Fl_O:tot:P', units='psi')[0] + v0 = prob.get_val('fc.Fl_O:stat:V', units='ft/s')[0] + pt2 = prob.get_val('inlet.Fl_O:tot:P', units='psi')[0] + pt4 = prob.get_val('burner.Fl_O:tot:P', units='psi')[0] + tt4 = prob.get_val('burner.Fl_O:tot:T', units='degR')[0] + wfuel = prob.get_val('burner.Wfuel', units='lbm/s')[0] + fram = prob.get_val('inlet.F_ram', units='lbf')[0] + fg = prob.get_val('nozz.Fg', units='lbf')[0] + fn = prob.get_val('perf.Fn', units='lbf')[0] + tsfc = prob.get_val('perf.TSFC', units='lbm/(h*lbf)')[0] + isp = fn / wfuel # Specific impulse (lbf-s / lbm) + pr_burner = pt4 / pt2 + throat_area = prob.get_val('nozz.Throat:stat:area', units='inch**2')[0] + exit_area = prob.get_val('nozz.Fl_O:stat:area', units='inch**2')[0] + + diagnostics = { + 'use_rde': use_rde, + 'P0_psi': p0, + 'Pt0_psi': pt0, + 'V0_fts': v0, + 'Pt2_psi': pt2, + 'Pt4_psi': pt4, + 'PR_burner': pr_burner, + 'Tt4_degR': tt4, + 'W_air': w_air, + 'W_fuel': wfuel, + 'F_ram': fram, + 'F_g': fg, + 'F_n': fn, + 'TSFC': tsfc, + 'Isp_s': isp, + 'Throat_area_in2': throat_area, + 'Exit_area_in2': exit_area + } + + if use_rde: + diagnostics['D_cj'] = prob.get_val('burner.D_cj', units='ft/s')[0] + diagnostics['f_rde'] = prob.get_val('burner.f_rde', units='Hz')[0] + diagnostics['Pt_inj'] = prob.get_val('burner.Pt_inj', units='psi')[0] + + return prob, diagnostics + + +def viewer(prob, pt, file=sys.stdout): + """ + print a report of all the relevant cycle properties + """ + + summary_data = (prob[pt+'.fc.Fl_O:stat:MN'], prob[pt+'.fc.alt'], prob[pt+'.inlet.Fl_O:stat:W'], + prob[pt+'.perf.Fn'], prob[pt+'.perf.Fg'], prob[pt+'.inlet.F_ram'], + prob[pt+'.perf.OPR'], prob[pt+'.perf.TSFC']) + summary_data = tuple(float(np.asarray(x).item()) if hasattr(x, 'item') else float(x) for x in summary_data) + + print(file=file, flush=True) + print(file=file, flush=True) + print(file=file, flush=True) + print("----------------------------------------------------------------------------", file=file, flush=True) + print(" POINT:", pt, file=file, flush=True) + print("----------------------------------------------------------------------------", file=file, flush=True) + print(" PERFORMANCE CHARACTERISTICS", file=file, flush=True) + print(" Mach Alt W Fn Fg Fram OPR TSFC ", file=file, flush=True) + print(" %7.5f %7.1f %7.3f %7.1f %7.1f %7.1f %7.3f %7.5f" % summary_data, file=file, flush=True) + + fs_names = ['fc.Fl_O', 'inlet.Fl_O', 'burner.Fl_O', 'nozz.Fl_O'] + fs_full_names = [f'{pt}.{fs}' for fs in fs_names] + pyc.print_flow_station(prob, fs_full_names, file=file) + + burner = prob.model._get_subsystem(f'{pt}.burner') + if isinstance(burner, RDECombustor): + pyc.print_rde(prob, [f'{pt}.burner'], file=file) + else: + pyc.print_burner(prob, [f'{pt}.burner'], file=file) + + noz_names = ['nozz'] + noz_full_names = [f'{pt}.{n}' for n in noz_names] + pyc.print_nozzle(prob, noz_full_names, file=file) + + +class MPRDERamjet(pyc.MPCycle): + + def setup(self): + self.pyc_add_pnt('DESIGN', RamjetCycle(use_rde=True, fuel_type='FAR')) + + self.set_input_defaults('DESIGN.inlet.MN', 0.35) + self.set_input_defaults('DESIGN.burner.MN', 0.30) + self.set_input_defaults('DESIGN.inlet.Fl_I:stat:W', 100.0, units='lbm/s') + self.set_input_defaults('DESIGN.burner.Fl_I:FAR', 0.028) + self.set_input_defaults('DESIGN.inlet.ram_recovery', 0.90) + + self.pyc_add_cycle_param('burner.dPqP_inj', 0.12) + self.pyc_add_cycle_param('burner.eta_rde', 0.85) + self.pyc_add_cycle_param('burner.dia_annulus', 14.0) + self.pyc_add_cycle_param('burner.N_waves', 1.0) + self.pyc_add_cycle_param('nozz.Cv', 0.98) + + self.od_pts = [] + + super().setup() + + +RDERamjet = RamjetCycle + + +if __name__ == "__main__": + + import time + + prob = om.Problem() + + mp_ramjet = prob.model = MPRDERamjet() + + prob.setup(check=False) + + # Define the design point + prob.set_val('DESIGN.fc.alt', 40000.0, units='ft') + prob.set_val('DESIGN.fc.MN', 2.50) + + st = time.time() + + prob.set_solver_print(level=-1) + prob.set_solver_print(level=2, depth=1) + + prob.run_model() + + for pt in ['DESIGN'] + mp_ramjet.od_pts: + viewer(prob, pt) + + print() + print("time", time.time() - st) diff --git a/example_cycles/rde_turbojet.py b/example_cycles/rde_turbojet.py new file mode 100644 index 00000000..ca01f348 --- /dev/null +++ b/example_cycles/rde_turbojet.py @@ -0,0 +1,359 @@ +""" +Step 6: Low-CPR Rotating Detonation Engine (RDE) Turbojet Cycle. + +Simulates and compares hybrid RDE-gas turbine cycle configurations: +1. Conventional Turbojet (CPR = 13.5, Isobaric Combustor with 3% pressure loss). +2. RDE Turbojet (CPR = 13.5, RDECombustor with pressure gain). +3. Low-CPR RDE Turbojet (CPR = 8.0, RDECombustor with pressure gain, reduced stages/weight). + +Architecture: +[FlightConditions] -> [Inlet] -> [Compressor] -> [RDECombustor / Combustor] -> [Turbine] -> [Nozzle] + ^ | + +-----------------[Shaft]--------------------+ +""" + +import sys +import numpy as np +import openmdao.api as om +import pycycle.api as pyc + +from pycycle.elements.rde_combustor import RDECombustor, print_rde + + +class HybridTurbojet(pyc.Cycle): + """ + Single-spool turbojet cycle supporting either a conventional or RDE combustor. + """ + + def initialize(self): + self.options.declare('use_rde', default=True, desc='Use RDE combustor if True, conventional if False') + self.options.declare('fuel_type', default='FAR', desc='Fuel type for thermodynamics') + super().initialize() + + def setup(self): + use_rde = self.options['use_rde'] + fuel_type = self.options['fuel_type'] + design = self.options['design'] + + # Fast TABULAR thermodynamics + self.options['thermo_method'] = 'TABULAR' + self.options['thermo_data'] = pyc.AIR_JETA_TAB_SPEC + + # 1. Add cycle components + self.add_subsystem('fc', pyc.FlightConditions()) + self.add_subsystem('inlet', pyc.Inlet()) + self.add_subsystem('comp', pyc.Compressor(map_data=pyc.AXI5, map_extrap=True), + promotes_inputs=['Nmech']) + + if use_rde: + self.add_subsystem('burner', RDECombustor(fuel_type=fuel_type, dia_annulus=12.0)) + else: + self.add_subsystem('burner', pyc.Combustor(fuel_type=fuel_type)) + + self.add_subsystem('turb', pyc.Turbine(map_data=pyc.LPT2269), + promotes_inputs=['Nmech']) + self.add_subsystem('nozz', pyc.Nozzle(nozzType='CD', lossCoef='Cv')) + self.add_subsystem('shaft', pyc.Shaft(num_ports=2), promotes_inputs=['Nmech']) + self.add_subsystem('perf', pyc.Performance(num_nozzles=1, num_burners=1)) + + # 2. Connect flow stations + self.pyc_connect_flow('fc.Fl_O', 'inlet.Fl_I', connect_w=False) + self.pyc_connect_flow('inlet.Fl_O', 'comp.Fl_I') + self.pyc_connect_flow('comp.Fl_O', 'burner.Fl_I') + self.pyc_connect_flow('burner.Fl_O', 'turb.Fl_I') + self.pyc_connect_flow('turb.Fl_O', 'nozz.Fl_I') + + # 3. Turbomachinery torque to shaft + self.connect('comp.trq', 'shaft.trq_0') + self.connect('turb.trq', 'shaft.trq_1') + + # 4. Nozzle exhaust backpressure + self.connect('fc.Fl_O:stat:P', 'nozz.Ps_exhaust') + + # 5. Performance connections + self.connect('inlet.Fl_O:tot:P', 'perf.Pt2') + self.connect('comp.Fl_O:tot:P', 'perf.Pt3') + self.connect('burner.Wfuel', 'perf.Wfuel_0') + self.connect('inlet.F_ram', 'perf.ram_drag') + self.connect('nozz.Fg', 'perf.Fg_0') + + # 6. Balances + balance = self.add_subsystem('balance', om.BalanceComp()) + if design: + # Match airflow to design net thrust Fn + balance.add_balance('W', units='lbm/s', eq_units='lbf', rhs_name='Fn_target') + self.connect('balance.W', 'inlet.Fl_I:stat:W') + self.connect('perf.Fn', 'balance.lhs:W') + + # Match FAR to design turbine inlet temperature T4 + balance.add_balance('FAR', eq_units='degR', lower=1e-4, val=0.0175, rhs_name='T4_target') + self.connect('balance.FAR', 'burner.Fl_I:FAR') + self.connect('burner.Fl_O:tot:T', 'balance.lhs:FAR') + + # Match turbine PR to net shaft power = 0 + balance.add_balance('turb_PR', val=3.5, lower=1.001, upper=12.0, eq_units='hp', rhs_val=0.0) + self.connect('balance.turb_PR', 'turb.PR') + self.connect('shaft.pwr_net', 'balance.lhs:turb_PR') + + else: + # Off-design: match FAR to target thrust + balance.add_balance('FAR', eq_units='lbf', lower=1e-4, val=0.02, rhs_name='Fn_target') + self.connect('balance.FAR', 'burner.Fl_I:FAR') + self.connect('perf.Fn', 'balance.lhs:FAR') + + # Shaft mechanical speed balance + balance.add_balance('Nmech', val=8070.0, units='rpm', lower=500.0, eq_units='hp', rhs_val=0.0) + self.connect('balance.Nmech', 'Nmech') + self.connect('shaft.pwr_net', 'balance.lhs:Nmech') + + # Choked throat area continuity balance + balance.add_balance('W', val=140.0, units='lbm/s', eq_units='inch**2') + self.connect('balance.W', 'inlet.Fl_I:stat:W') + self.connect('nozz.Throat:stat:area', 'balance.lhs:W') + + # 7. Non-linear solver configuration + newton = self.nonlinear_solver = om.NewtonSolver() + newton.options['atol'] = 1e-6 + newton.options['rtol'] = 1e-6 + newton.options['iprint'] = -1 + newton.options['maxiter'] = 40 + newton.options['solve_subsystems'] = True + newton.options['max_sub_solves'] = 100 + newton.options['reraise_child_analysiserror'] = False + + self.linear_solver = om.DirectSolver() + + super().setup() + + +class MPHybridTurbojet(pyc.MPCycle): + """ + Multi-point cycle wrapper for design and off-design evaluation. + """ + + def initialize(self): + self.options.declare('use_rde', default=True) + super().initialize() + + def setup(self): + use_rde = self.options['use_rde'] + + # Add DESIGN point + self.pyc_add_pnt('DESIGN', HybridTurbojet(use_rde=use_rde)) + + self.set_input_defaults('DESIGN.Nmech', 8070.0, units='rpm') + self.set_input_defaults('DESIGN.inlet.MN', 0.60) + self.set_input_defaults('DESIGN.comp.MN', 0.020) + self.set_input_defaults('DESIGN.burner.MN', 0.020) + self.set_input_defaults('DESIGN.turb.MN', 0.40) + + if use_rde: + self.pyc_add_cycle_param('burner.dPqP_inj', 0.12) + self.pyc_add_cycle_param('burner.eta_rde', 0.85) + self.pyc_add_cycle_param('burner.dia_annulus', 12.0) + self.pyc_add_cycle_param('burner.N_waves', 1.0) + else: + self.pyc_add_cycle_param('burner.dPqP', 0.03) + + self.pyc_add_cycle_param('nozz.Cv', 0.99) + + # Off-design point (OD0: Mach 0.20 at 5,000 ft) + self.od_pts = ['OD0'] + self.pyc_add_pnt('OD0', HybridTurbojet(design=False, use_rde=use_rde)) + self.set_input_defaults('OD0.fc.MN', val=0.20) + self.set_input_defaults('OD0.fc.alt', 5000.0, units='ft') + self.set_input_defaults('OD0.balance.Fn_target', 8000.0, units='lbf') + + self.pyc_use_default_des_od_conns() + self.pyc_connect_des_od('nozz.Throat:stat:area', 'balance.rhs:W') + + super().setup() + + +def run_turbojet_simulation(use_rde=True, cpr=13.5, fn_target=11800.0, t4_target=2370.0): + """ + Sets up and solves the turbojet cycle for a given CPR and burner configuration. + """ + prob = om.Problem() + prob.model = MPHybridTurbojet(use_rde=use_rde) + + prob.set_solver_print(level=-1) + prob.setup(check=False) + + # Design point conditions + prob.set_val('DESIGN.fc.alt', 0.0, units='ft') + prob.set_val('DESIGN.fc.MN', 0.000001) + prob.set_val('DESIGN.balance.Fn_target', fn_target, units='lbf') + prob.set_val('DESIGN.balance.T4_target', t4_target, units='degR') + prob.set_val('DESIGN.comp.PR', cpr) + prob.set_val('DESIGN.comp.eff', 0.83) + prob.set_val('DESIGN.turb.eff', 0.86) + + # Balance initial guesses + prob['DESIGN.balance.W'] = 145.0 if not use_rde else 125.0 + prob['DESIGN.balance.FAR'] = 0.01755 + prob['DESIGN.balance.turb_PR'] = 3.85 if not use_rde else 4.20 + prob['DESIGN.fc.balance.Pt'] = 14.696 + prob['DESIGN.fc.balance.Tt'] = 518.67 + + for pt in ['OD0']: + prob[pt + '.balance.W'] = 140.0 if not use_rde else 120.0 + prob[pt + '.balance.FAR'] = 0.0168 + prob[pt + '.balance.Nmech'] = 8100.0 + prob[pt + '.fc.balance.Pt'] = 13.0 + prob[pt + '.fc.balance.Tt'] = 530.0 + prob[pt + '.turb.PR'] = 4.0 + + prob.run_model() + + pt = 'DESIGN' + pt2 = prob.get_val(f'{pt}.inlet.Fl_O:tot:P', units='psi')[0] + pt3 = prob.get_val(f'{pt}.comp.Fl_O:tot:P', units='psi')[0] + pt4 = prob.get_val(f'{pt}.burner.Fl_O:tot:P', units='psi')[0] + tt4 = prob.get_val(f'{pt}.burner.Fl_O:tot:T', units='degR')[0] + pt5 = prob.get_val(f'{pt}.turb.Fl_O:tot:P', units='psi')[0] + w_air = prob.get_val(f'{pt}.inlet.Fl_O:stat:W', units='lbm/s')[0] + w_fuel = prob.get_val(f'{pt}.perf.Wfuel', units='lbm/s')[0] + fn = prob.get_val(f'{pt}.perf.Fn', units='lbf')[0] + tsfc = prob.get_val(f'{pt}.perf.TSFC', units='lbm/(h*lbf)')[0] + turb_pr = prob.get_val(f'{pt}.turb.PR')[0] + throat_area = prob.get_val(f'{pt}.nozz.Throat:stat:area', units='inch**2')[0] + comp_pwr = prob.get_val(f'{pt}.shaft.pwr_out', units='hp')[0] + + diagnostics = { + 'use_rde': use_rde, + 'CPR': cpr, + 'Pt2_psi': pt2, + 'Pt3_psi': pt3, + 'Pt4_psi': pt4, + 'PR_burner': pt4 / pt3, + 'Pt5_psi': pt5, + 'Tt4_degR': tt4, + 'W_air': w_air, + 'W_fuel': w_fuel, + 'Fn_lbf': fn, + 'TSFC': tsfc, + 'Turb_PR': turb_pr, + 'Comp_pwr_hp': comp_pwr, + 'Throat_area_in2': throat_area + } + + if use_rde: + diagnostics['D_cj'] = prob.get_val(f'{pt}.burner.D_cj', units='ft/s')[0] + diagnostics['f_rde'] = prob.get_val(f'{pt}.burner.f_rde', units='Hz')[0] + diagnostics['Pt_inj'] = prob.get_val(f'{pt}.burner.Pt_inj', units='psi')[0] + + return prob, diagnostics + + +def viewer(prob, pt, file=sys.stdout): + """ + print a report of all the relevant cycle properties + """ + + summary_data = (prob[pt+'.fc.Fl_O:stat:MN'], prob[pt+'.fc.alt'], prob[pt+'.inlet.Fl_O:stat:W'], + prob[pt+'.perf.Fn'], prob[pt+'.perf.Fg'], prob[pt+'.inlet.F_ram'], + prob[pt+'.perf.OPR'], prob[pt+'.perf.TSFC']) + summary_data = tuple(float(np.asarray(x).item()) if hasattr(x, 'item') else float(x) for x in summary_data) + + print(file=file, flush=True) + print(file=file, flush=True) + print(file=file, flush=True) + print("----------------------------------------------------------------------------", file=file, flush=True) + print(" POINT:", pt, file=file, flush=True) + print("----------------------------------------------------------------------------", file=file, flush=True) + print(" PERFORMANCE CHARACTERISTICS", file=file, flush=True) + print(" Mach Alt W Fn Fg Fram OPR TSFC ", file=file, flush=True) + print(" %7.5f %7.1f %7.3f %7.1f %7.1f %7.1f %7.3f %7.5f" % summary_data, file=file, flush=True) + + fs_names = ['fc.Fl_O', 'inlet.Fl_O', 'comp.Fl_O', 'burner.Fl_O', + 'turb.Fl_O', 'nozz.Fl_O'] + fs_full_names = [f'{pt}.{fs}' for fs in fs_names] + pyc.print_flow_station(prob, fs_full_names, file=file) + + comp_names = ['comp'] + comp_full_names = [f'{pt}.{c}' for c in comp_names] + pyc.print_compressor(prob, comp_full_names, file=file) + + burner = prob.model._get_subsystem(f'{pt}.burner') + if isinstance(burner, RDECombustor): + pyc.print_rde(prob, [f'{pt}.burner'], file=file) + else: + pyc.print_burner(prob, [f'{pt}.burner'], file=file) + + turb_names = ['turb'] + turb_full_names = [f'{pt}.{t}' for t in turb_names] + pyc.print_turbine(prob, turb_full_names, file=file) + + noz_names = ['nozz'] + noz_full_names = [f'{pt}.{n}' for n in noz_names] + pyc.print_nozzle(prob, noz_full_names, file=file) + + shaft_names = ['shaft'] + shaft_full_names = [f'{pt}.{s}' for s in shaft_names] + pyc.print_shaft(prob, shaft_full_names, file=file) + + pyc.print_balances(prob, pt, file=file) + + +def map_plots(prob, pt): + comp_names = ['comp'] + comp_full_names = [f'{pt}.{c}' for c in comp_names] + pyc.plot_compressor_maps(prob, comp_full_names) + + turb_names = ['turb'] + turb_full_names = [f'{pt}.{c}' for c in turb_names] + pyc.plot_turbine_maps(prob, turb_full_names) + + +RDETurbojet = HybridTurbojet +MPRDETurbojet = MPHybridTurbojet + + +if __name__ == "__main__": + + import time + + prob = om.Problem() + + mp_turbojet = prob.model = MPHybridTurbojet(use_rde=True) + + prob.setup(check=False) + + # Define the design point + prob.set_val('DESIGN.fc.alt', 0.0, units='ft') + prob.set_val('DESIGN.fc.MN', 0.000001) + prob.set_val('DESIGN.balance.Fn_target', 11800.0, units='lbf') + prob.set_val('DESIGN.balance.T4_target', 2370.0, units='degR') + prob.set_val('DESIGN.comp.PR', 13.5) + prob.set_val('DESIGN.comp.eff', 0.83) + prob.set_val('DESIGN.turb.eff', 0.86) + + # Set initial guesses for balances + prob['DESIGN.balance.FAR'] = 0.01755 + prob['DESIGN.balance.W'] = 125.0 + prob['DESIGN.balance.turb_PR'] = 4.20 + prob['DESIGN.fc.balance.Pt'] = 14.696 + prob['DESIGN.fc.balance.Tt'] = 518.67 + + for i, pt in enumerate(mp_turbojet.od_pts): + # initial guesses + prob[pt + '.balance.W'] = 120.0 + prob[pt + '.balance.FAR'] = 0.01680 + prob[pt + '.balance.Nmech'] = 8197.38 + prob[pt + '.fc.balance.Pt'] = 15.703 + prob[pt + '.fc.balance.Tt'] = 558.31 + prob[pt + '.turb.PR'] = 4.6690 + + st = time.time() + + prob.set_solver_print(level=-1) + prob.set_solver_print(level=2, depth=1) + + prob.run_model() + + for pt in ['DESIGN'] + mp_turbojet.od_pts: + viewer(prob, pt) + + print() + print("time", time.time() - st) diff --git a/example_cycles/tests/benchmark_rde_ramjet.py b/example_cycles/tests/benchmark_rde_ramjet.py new file mode 100644 index 00000000..83209c13 --- /dev/null +++ b/example_cycles/tests/benchmark_rde_ramjet.py @@ -0,0 +1,67 @@ +import numpy as np +import unittest + +import openmdao.api as om +from openmdao.utils.assert_utils import assert_near_equal + +import pycycle.api as pyc +from example_cycles.rde_ramjet import MPRDERamjet + + +class RDERamjetTestCase(unittest.TestCase): + + def test_benchmark(self): + self.benchmark_case1() + + def benchmark_case1(self): + prob = om.Problem() + mp_ramjet = prob.model = MPRDERamjet() + + prob.set_solver_print(level=-1) + prob.setup(check=False) + + # Design flight condition: Mach 2.50 at 40,000 ft + prob.set_val('DESIGN.fc.alt', 40000.0, units='ft') + prob.set_val('DESIGN.fc.MN', 2.50) + + old = np.seterr(divide='raise') + try: + prob.run_model() + tol = 1e-4 + + # Inlet airflow + assert_near_equal(prob['DESIGN.inlet.Fl_O:stat:W'][0], 100.0, tol) + + # Overall pressure ratio (due to RDE detonation pressure rise) + assert_near_equal(prob['DESIGN.perf.OPR'][0], 2.9628, tol) + + # Net thrust Fn = Fg - Fram + assert_near_equal(prob['DESIGN.perf.Fn'][0], 7660.7, tol) + + # Gross thrust Fg + assert_near_equal(prob['DESIGN.perf.Fg'][0], 15185.85, tol) + + # Ram drag + assert_near_equal(prob['DESIGN.inlet.F_ram'][0], 7525.2, tol) + + # TSFC + assert_near_equal(prob['DESIGN.perf.TSFC'][0], 1.3158, tol) + + # Combustor pressure ratio (Pt4 / Pt2) + assert_near_equal(prob['DESIGN.burner.PR_RDE'][0], 2.9628, tol) + + # Injector plenum pressure + assert_near_equal(prob['DESIGN.burner.Pt_inj'][0], 37.186, tol) + + # Detonation wave frequency + assert_near_equal(prob['DESIGN.burner.f_rde'][0], 1317.9, tol) + + # Nozzle throat area + assert_near_equal(prob['DESIGN.nozz.Throat:stat:area'][0], 82.733, tol) + + finally: + np.seterr(**old) + + +if __name__ == "__main__": + unittest.main() diff --git a/example_cycles/tests/benchmark_rde_turbojet.py b/example_cycles/tests/benchmark_rde_turbojet.py new file mode 100644 index 00000000..0c403ed5 --- /dev/null +++ b/example_cycles/tests/benchmark_rde_turbojet.py @@ -0,0 +1,101 @@ +import numpy as np +import unittest + +import openmdao.api as om +from openmdao.utils.assert_utils import assert_near_equal + +import pycycle.api as pyc +from example_cycles.rde_turbojet import MPRDETurbojet + + +class RDETurbojetTestCase(unittest.TestCase): + + def test_benchmark(self): + self.benchmark_case1() + + def benchmark_case1(self): + prob = om.Problem() + mp_turbojet = prob.model = MPRDETurbojet(use_rde=True) + + prob.set_solver_print(level=-1) + prob.setup(check=False) + + # Initial Conditions for DESIGN + prob.set_val('DESIGN.fc.alt', 0.0, units='ft') + prob.set_val('DESIGN.fc.MN', 0.000001) + prob.set_val('DESIGN.balance.Fn_target', 11800.0, units='lbf') + prob.set_val('DESIGN.balance.T4_target', 2370.0, units='degR') + prob.set_val('DESIGN.comp.PR', 13.5) + prob.set_val('DESIGN.comp.eff', 0.83) + prob.set_val('DESIGN.turb.eff', 0.86) + + # Balance initial guesses + prob['DESIGN.balance.FAR'] = 0.01755 + prob['DESIGN.balance.W'] = 125.0 + prob['DESIGN.balance.turb_PR'] = 4.20 + prob['DESIGN.fc.balance.Pt'] = 14.696 + prob['DESIGN.fc.balance.Tt'] = 518.67 + + for i, pt in enumerate(mp_turbojet.od_pts): + prob[pt + '.balance.W'] = 120.0 + prob[pt + '.balance.FAR'] = 0.01680 + prob[pt + '.balance.Nmech'] = 8197.38 + prob[pt + '.fc.balance.Pt'] = 15.703 + prob[pt + '.fc.balance.Tt'] = 558.31 + prob[pt + '.turb.PR'] = 4.6690 + + old = np.seterr(divide='raise') + try: + prob.run_model() + tol = 1e-4 + + # --- DESIGN Point Regressions --- + # Airflow W + assert_near_equal(prob['DESIGN.inlet.Fl_O:stat:W'][0], 123.6359, tol) + + # Overall Pressure Ratio + assert_near_equal(prob['DESIGN.perf.OPR'][0], 13.5000, tol) + + # Fuel to Air Ratio + assert_near_equal(prob['DESIGN.balance.FAR'][0], 0.017765, tol) + + # Turbine Expansion Ratio + assert_near_equal(prob['DESIGN.balance.turb_PR'][0], 3.8546, tol) + + # Gross Thrust + assert_near_equal(prob['DESIGN.perf.Fg'][0], 11799.99, tol) + + # TSFC (16% reduction vs conventional Brayton 0.7985) + assert_near_equal(prob['DESIGN.perf.TSFC'][0], 0.67008, tol) + + # Combustor Inlet Temperature + assert_near_equal(prob['DESIGN.comp.Fl_O:tot:T'][0], 1187.76, tol) + + # RDE Combustor Pressure Ratio (Pressure Gain Combustion: Pt4/Pt3 = 1.863) + assert_near_equal(prob['DESIGN.burner.PR_RDE'][0], 1.8633, tol) + + # Injector plenum pressure + assert_near_equal(prob['DESIGN.burner.Pt_inj'][0], 174.587, tol) + + # Detonation wave speed + assert_near_equal(prob['DESIGN.burner.D_cj'][0], 4184.15, tol) + + # Detonation frequency + assert_near_equal(prob['DESIGN.burner.f_rde'][0], 1331.86, tol) + + # --- OD0 Off-Design Point Regressions --- + assert_near_equal(prob['OD0.inlet.Fl_O:stat:W'][0], 97.1832, tol) + assert_near_equal(prob['OD0.perf.OPR'][0], 12.2789, tol) + assert_near_equal(prob['OD0.balance.FAR'][0], 0.015692, tol) + assert_near_equal(prob['OD0.balance.Nmech'][0], 7621.79, tol) + assert_near_equal(prob['OD0.perf.Fg'][0], 8662.85, tol) + assert_near_equal(prob['OD0.perf.TSFC'][0], 0.68625, tol) + assert_near_equal(prob['OD0.comp.Fl_O:tot:T'][0], 1119.53, tol) + assert_near_equal(prob['OD0.burner.PR_RDE'][0], 1.8033, tol) + + finally: + np.seterr(**old) + + +if __name__ == "__main__": + unittest.main() diff --git a/pycycle/api.py b/pycycle/api.py index 8b057d97..f921f9bb 100644 --- a/pycycle/api.py +++ b/pycycle/api.py @@ -11,6 +11,7 @@ from pycycle.elements.duct import Duct from pycycle.elements.compressor import Compressor from pycycle.elements.combustor import Combustor +from pycycle.elements.rde_combustor import RDECombustor, RDEPressureGainComp, print_rde from pycycle.elements.turbine import Turbine from pycycle.elements.nozzle import Nozzle from pycycle.elements.shaft import Shaft diff --git a/pycycle/elements/rde_combustor.py b/pycycle/elements/rde_combustor.py new file mode 100644 index 00000000..c7d79d66 --- /dev/null +++ b/pycycle/elements/rde_combustor.py @@ -0,0 +1,439 @@ +""" +Rotating Detonation Engine (RDE) Combustor Element for pyCycle. + +Implements the pressure-gain combustion physics of an RDE, including: +1. Injector plenum feed pressure drop +2. Humphrey / Fickett-Jacobs detonation thermodynamics +3. Detonation wave kinematics (CJ wave speed and rotation frequency) +4. Unsteady-to-steady flow realization efficiency +5. Integration with pyCycle Element API (ThermoAdd + Thermo) +""" + +import sys +import numpy as np +import openmdao.api as om + +from pycycle.constants import R_UNIVERSAL_ENG, g_c +from pycycle.thermo.thermo import Thermo, ThermoAdd +from pycycle.flow_in import FlowIn +from pycycle.passthrough import PassThrough +from pycycle.element_base import Element + + +class RDEPressureGainComp(om.ExplicitComponent): + """ + Computes the stagnation pressure gain and wave kinematics across an RDE annulus. + + Governing Equations: + -------------------- + 1. Injector Plenum Pressure: + Pt_inj = Pt_in * (1.0 - dPqP_inj) + + 2. Fuel Energy Release: + q = (FAR * Q_fuel) / (1.0 + FAR) [Btu/lbm] + + 3. Gas Properties: + R_gas = R_UNIVERSAL_ENG / MW_gas [Btu/(lbm*degR)] + cv = R_gas / (gamma_gas - 1.0) + + 4. Humphrey Pressure Ratio (Ideal Constant-Volume Combustion): + Pi_ideal = 1.0 + q / (cv * Tt_in) = 1.0 + (gamma_gas - 1.0) * q / (R_gas * Tt_in) + + 5. Realized Stagnation Pressure Gain: + Pi_det = 1.0 + eta_rde * (Pi_ideal - 1.0) + Pt_out = Pt_inj * Pi_det = Pt_in * (1.0 - dPqP_inj) * Pi_det + PR_RDE = Pt_out / Pt_in = (1.0 - dPqP_inj) * Pi_det + + 6. Wave Kinematics: + M_cj = sqrt(1.0 + (gamma_gas**2 - 1.0)/(2.0*gamma_gas) * q/(R_gas*Tt_in)) + ... + D_cj = M_cj * sqrt(gamma_gas * R_gas * 778.169 * g_c * Tt_in) [ft/s] + f_rde = (N_waves * D_cj) / (pi * (dia_annulus / 12.0)) [Hz] + """ + + def initialize(self): + self.options.declare('detonation_mode', default='HUMPHREY', values=['HUMPHREY', 'CJ'], + desc='Thermodynamic mode for detonation pressure rise') + + def setup(self): + # Stagnation flow inputs from upstream + self.add_input('Pt_in', val=100.0, units='lbf/inch**2', desc='Inlet total pressure') + self.add_input('Tt_in', val=1000.0, units='degR', desc='Inlet total temperature') + self.add_input('FAR', val=0.03, desc='Fuel to air ratio') + + # Combustor design parameters + self.add_input('dPqP_inj', val=0.12, desc='Injector feed pressure loss fraction') + self.add_input('eta_rde', val=0.85, desc='Detonation combustor realization efficiency') + self.add_input('Q_fuel', val=18600.0, units='Btu/lbm', desc='Lower heating value of fuel') + self.add_input('gamma_gas', val=1.33, desc='Effective ratio of specific heats') + self.add_input('MW_gas', val=28.96, units='lbm/mol', desc='Molecular weight of gas mixture') + self.add_input('dia_annulus', val=12.0, units='inch', desc='Mean diameter of RDE annulus') + self.add_input('N_waves', val=1.0, desc='Number of co-rotating detonation waves') + + # Performance outputs + self.add_output('Pt_out', val=110.0, units='lbf/inch**2', desc='Exit total pressure') + self.add_output('PR_RDE', val=1.10, desc='Net combustor pressure ratio (Pt_out / Pt_in)') + self.add_output('Pt_inj', val=88.0, units='lbf/inch**2', desc='Injector plenum total pressure') + self.add_output('D_cj', val=5500.0, units='ft/s', desc='Chapman-Jouguet detonation wave speed') + self.add_output('f_rde', val=1750.0, units='Hz', desc='Detonation wave rotational frequency') + + # Analytic derivatives setup + self.declare_partials('Pt_out', ['Pt_in', 'Tt_in', 'FAR', 'dPqP_inj', 'eta_rde', 'Q_fuel', 'gamma_gas', 'MW_gas']) + self.declare_partials('PR_RDE', ['Tt_in', 'FAR', 'dPqP_inj', 'eta_rde', 'Q_fuel', 'gamma_gas', 'MW_gas']) + self.declare_partials('Pt_inj', ['Pt_in', 'dPqP_inj']) + self.declare_partials('D_cj', ['Tt_in', 'FAR', 'Q_fuel', 'gamma_gas', 'MW_gas']) + self.declare_partials('f_rde', ['Tt_in', 'FAR', 'Q_fuel', 'gamma_gas', 'MW_gas', 'dia_annulus', 'N_waves']) + + def compute(self, inputs, outputs): + pt_in = inputs['Pt_in'] + tt_in = inputs['Tt_in'] + far = inputs['FAR'] + dpqp_inj = inputs['dPqP_inj'] + eta = inputs['eta_rde'] + q_fuel = inputs['Q_fuel'] + gamma = inputs['gamma_gas'] + mw = inputs['MW_gas'] + dia = inputs['dia_annulus'] + n_waves = inputs['N_waves'] + + # 1. Injector plenum pressure + u_inj = 1.0 - dpqp_inj + pt_inj = pt_in * u_inj + outputs['Pt_inj'] = pt_inj + + # 2. Specific heat release + q = (far * q_fuel) / (1.0 + far) + + # 3. Gas properties + r_gas = R_UNIVERSAL_ENG / mw # Btu/(lbm * degR) + + # 4. Ideal Humphrey pressure ratio + pi_ideal_minus_1 = (gamma - 1.0) * q / (r_gas * tt_in) + + # 5. Realized pressure gain + pi_det = 1.0 + eta * pi_ideal_minus_1 + pr_rde = u_inj * pi_det + pt_out = pt_in * pr_rde + + outputs['Pt_out'] = pt_out + outputs['PR_RDE'] = pr_rde + + # 6. Wave kinematics + # 1 Btu = 778.169 ft*lbf, g_c = 32.174 + conv_factor = 778.169 * g_c + q_mech = q * conv_factor # (ft*lbf)/lbm + r_mech = r_gas * conv_factor # (ft*lbf)/(lbm * degR) + a_inj = np.sqrt(gamma * r_mech * tt_in) + + xi = (gamma**2 - 1.0) * q_mech / (2.0 * gamma * r_mech * tt_in) + m_cj = np.sqrt(1.0 + xi) + np.sqrt(xi) + d_cj = m_cj * a_inj + outputs['D_cj'] = d_cj + + dia_ft = dia / 12.0 + circ = np.pi * dia_ft + outputs['f_rde'] = (n_waves * d_cj) / circ + + def compute_partials(self, inputs, J): + pt_in = inputs['Pt_in'] + tt_in = inputs['Tt_in'] + far = inputs['FAR'] + dpqp_inj = inputs['dPqP_inj'] + eta = inputs['eta_rde'] + q_fuel = inputs['Q_fuel'] + gamma = inputs['gamma_gas'] + mw = inputs['MW_gas'] + dia = inputs['dia_annulus'] + n_waves = inputs['N_waves'] + + u_inj = 1.0 - dpqp_inj + q = (far * q_fuel) / (1.0 + far) + r_gas = R_UNIVERSAL_ENG / mw + pi_ideal_minus_1 = (gamma - 1.0) * q / (r_gas * tt_in) + pi_det = 1.0 + eta * pi_ideal_minus_1 + pr_rde = u_inj * pi_det + + # Pt_inj partials + J['Pt_inj', 'Pt_in'] = u_inj + J['Pt_inj', 'dPqP_inj'] = -pt_in + + # PR_RDE partials + J['PR_RDE', 'dPqP_inj'] = -pi_det + J['PR_RDE', 'eta_rde'] = u_inj * pi_ideal_minus_1 + + # Derivatives of pi_ideal_minus_1: + dq_dfar = q_fuel / (1.0 + far)**2 + dpi_dfar = (gamma - 1.0) * dq_dfar / (r_gas * tt_in) + dpi_dtt = - (gamma - 1.0) * q / (r_gas * tt_in**2) + dpi_dqfuel = (gamma - 1.0) * far / ((1.0 + far) * r_gas * tt_in) + dpi_dgamma = q / (r_gas * tt_in) + dpi_dmw = (gamma - 1.0) * q / (R_UNIVERSAL_ENG * tt_in) + + J['PR_RDE', 'FAR'] = u_inj * eta * dpi_dfar + J['PR_RDE', 'Tt_in'] = u_inj * eta * dpi_dtt + J['PR_RDE', 'Q_fuel'] = u_inj * eta * dpi_dqfuel + J['PR_RDE', 'gamma_gas'] = u_inj * eta * dpi_dgamma + J['PR_RDE', 'MW_gas'] = u_inj * eta * dpi_dmw + + # Pt_out partials + J['Pt_out', 'Pt_in'] = pr_rde + J['Pt_out', 'dPqP_inj'] = pt_in * J['PR_RDE', 'dPqP_inj'] + J['Pt_out', 'eta_rde'] = pt_in * J['PR_RDE', 'eta_rde'] + J['Pt_out', 'FAR'] = pt_in * J['PR_RDE', 'FAR'] + J['Pt_out', 'Tt_in'] = pt_in * J['PR_RDE', 'Tt_in'] + J['Pt_out', 'Q_fuel'] = pt_in * J['PR_RDE', 'Q_fuel'] + J['Pt_out', 'gamma_gas'] = pt_in * J['PR_RDE', 'gamma_gas'] + J['Pt_out', 'MW_gas'] = pt_in * J['PR_RDE', 'MW_gas'] + + # D_cj and f_rde partials + conv_factor = 778.169 * g_c + q_mech = q * conv_factor + r_mech = r_gas * conv_factor + a_inj = np.sqrt(gamma * r_mech * tt_in) + + xi = (gamma**2 - 1.0) * q_mech / (2.0 * gamma * r_mech * tt_in) + sqrt_1_xi = np.sqrt(1.0 + xi) + sqrt_xi = np.sqrt(max(xi, 1e-12)) + m_cj = sqrt_1_xi + sqrt_xi + d_cj = m_cj * a_inj + + # d(m_cj)/d(xi) + dm_dxi = 0.5 / sqrt_1_xi + 0.5 / sqrt_xi + + # d(xi) derivatives: + dxi_dq = (gamma**2 - 1.0) * conv_factor / (2.0 * gamma * r_mech * tt_in) + dxi_dfar = dxi_dq * dq_dfar + dxi_dqfuel = dxi_dq * (far / (1.0 + far)) + dxi_dtt = - xi / tt_in + dxi_dmw = xi / mw + dxi_dgamma = ((2.0 * gamma**2 - (gamma**2 - 1.0)) / (2.0 * gamma**2)) * (q_mech / (r_mech * tt_in)) + + # da_inj derivatives: + da_dtt = 0.5 * a_inj / tt_in + da_dgamma = 0.5 * a_inj / gamma + da_dmw = -0.5 * a_inj / mw + + # d(D_cj) + J['D_cj', 'FAR'] = (dm_dxi * dxi_dfar) * a_inj + J['D_cj', 'Q_fuel'] = (dm_dxi * dxi_dqfuel) * a_inj + J['D_cj', 'Tt_in'] = (dm_dxi * dxi_dtt) * a_inj + m_cj * da_dtt + J['D_cj', 'gamma_gas'] = (dm_dxi * dxi_dgamma) * a_inj + m_cj * da_dgamma + J['D_cj', 'MW_gas'] = (dm_dxi * dxi_dmw) * a_inj + m_cj * da_dmw + + # f_rde partials + dia_ft = dia / 12.0 + circ = np.pi * dia_ft + scale = n_waves / circ + + J['f_rde', 'FAR'] = scale * J['D_cj', 'FAR'] + J['f_rde', 'Q_fuel'] = scale * J['D_cj', 'Q_fuel'] + J['f_rde', 'Tt_in'] = scale * J['D_cj', 'Tt_in'] + J['f_rde', 'gamma_gas'] = scale * J['D_cj', 'gamma_gas'] + J['f_rde', 'MW_gas'] = scale * J['D_cj', 'MW_gas'] + J['f_rde', 'dia_annulus'] = - (n_waves * d_cj) / (np.pi * (dia_ft**2) * 12.0) + J['f_rde', 'N_waves'] = d_cj / circ + + +class RDECombustor(Element): + """ + Rotating Detonation Engine (RDE) Combustor Element for pyCycle. + + Replaces conventional isobaric combustors with pressure-gain detonation combustion. + + -------------- + Flow Stations + -------------- + Fl_I: Inlet flow from compressor/diffuser (Station 3) + Fl_O: Stagnation and static flow exit with pressure gain (Station 4) + + ------------- + Design Inputs + ------------- + Fl_I:FAR : Fuel to air ratio + dPqP_inj : Injector feed pressure loss fraction (default: 0.12) + eta_rde : Detonation realization efficiency (default: 0.85) + dia_annulus : Mean diameter of annulus [inch] + MN : Exit Mach number (design static calculation) + + ------------- + Off-Design Inputs + ------------- + Fl_I:FAR : Fuel to air ratio + dPqP_inj : Injector feed pressure loss fraction + eta_rde : Detonation realization efficiency + area : Combustor exit area (off-design matching) + + ------------- + Outputs + ------------- + Wfuel : Fuel mass flow rate [lbm/s] + PR_RDE : Net combustor stagnation pressure ratio (Pt4 / Pt3) + Pt_inj : Injector plenum total pressure [psi] + D_cj : Detonation wave speed [ft/s] + f_rde : Detonation rotation frequency [Hz] + """ + + def initialize(self): + self.options.declare('statics', default=True, + desc='If True, calculate static properties.') + self.options.declare('fuel_type', default="JP-7", + desc='Type of fuel.') + self.options.declare('detonation_mode', default='HUMPHREY', values=['HUMPHREY', 'CJ'], + desc='Thermodynamic mode for detonation calculation') + self.options.declare('dia_annulus', default=12.0, + desc='Mean diameter of RDE annulus [inch]') + + self.default_des_od_conns = [ + ('Fl_O:stat:area', 'area') + ] + + super().initialize() + + def pyc_setup_output_ports(self): + thermo_method = self.options['thermo_method'] + thermo_data = self.options['thermo_data'] + fuel_type = self.options['fuel_type'] + + self.thermo_add_comp = ThermoAdd( + method=thermo_method, + mix_mode='reactant', + thermo_kwargs={ + 'spec': thermo_data, + 'inflow_composition': self.Fl_I_data['Fl_I'], + 'mix_composition': fuel_type + } + ) + + self.copy_flow(self.thermo_add_comp, 'Fl_O') + + def setup(self): + thermo_method = self.options['thermo_method'] + thermo_data = self.options['thermo_data'] + design = self.options['design'] + statics = self.options['statics'] + air_fuel_composition = self.Fl_O_data['Fl_O'] + det_mode = self.options['detonation_mode'] + + # 1. Inlet flow station + in_flow = FlowIn(fl_name='Fl_I') + self.add_subsystem('in_flow', in_flow, promotes=['Fl_I:tot:*', 'Fl_I:stat:*']) + + # 2. Fuel mixing (ThermoAdd) + self.add_subsystem( + 'mix_fuel', + self.thermo_add_comp, + promotes=[ + 'Fl_I:stat:W', + ('mix:ratio', 'Fl_I:FAR'), + 'Fl_I:tot:composition', + 'Fl_I:tot:h', + ('mix:W', 'Wfuel'), + 'Wout' + ] + ) + + # 3. RDE Detonation Pressure Gain Component + rde_pg = RDEPressureGainComp(detonation_mode=det_mode) + prom_pg_in = [ + 'dPqP_inj', 'eta_rde', 'Q_fuel', 'gamma_gas', + 'MW_gas', 'dia_annulus', 'N_waves' + ] + prom_pg_out = ['PR_RDE', 'Pt_inj', 'D_cj', 'f_rde'] + self.add_subsystem('rde_press_gain', rde_pg, + promotes_inputs=prom_pg_in, + promotes_outputs=prom_pg_out) + + self.connect('Fl_I:tot:P', 'rde_press_gain.Pt_in') + self.connect('Fl_I:tot:T', 'rde_press_gain.Tt_in') + self.connect('Fl_I:FAR', 'rde_press_gain.FAR') + + # 4. Vitiated flow station (Thermo total_hP at elevated pressure Pt_out) + vit_flow = Thermo( + mode='total_hP', + fl_name='Fl_O:tot', + method=thermo_method, + thermo_kwargs={ + 'composition': air_fuel_composition, + 'spec': thermo_data + } + ) + self.add_subsystem('vitiated_flow', vit_flow, promotes_outputs=['Fl_O:*']) + self.connect('mix_fuel.mass_avg_h', 'vitiated_flow.h') + self.connect('mix_fuel.composition_out', 'vitiated_flow.composition') + self.connect('rde_press_gain.Pt_out', 'vitiated_flow.P') + + # 5. Static flow station properties + if statics: + if design: + out_stat = Thermo( + mode='static_MN', + fl_name='Fl_O:stat', + method=thermo_method, + thermo_kwargs={ + 'composition': air_fuel_composition, + 'spec': thermo_data + } + ) + self.add_subsystem('out_stat', out_stat, + promotes_inputs=['MN'], + promotes_outputs=['Fl_O:stat:*']) + self.connect('mix_fuel.composition_out', 'out_stat.composition') + self.connect('Fl_O:tot:S', 'out_stat.S') + self.connect('Fl_O:tot:h', 'out_stat.ht') + self.connect('Fl_O:tot:P', 'out_stat.guess:Pt') + self.connect('Fl_O:tot:gamma', 'out_stat.guess:gamt') + self.connect('Wout', 'out_stat.W') + + else: + out_stat = Thermo( + mode='static_A', + fl_name='Fl_O:stat', + method=thermo_method, + thermo_kwargs={ + 'composition': air_fuel_composition, + 'spec': thermo_data + } + ) + self.add_subsystem('out_stat', out_stat, + promotes_inputs=['area'], + promotes_outputs=['Fl_O:stat:*']) + self.connect('mix_fuel.composition_out', 'out_stat.composition') + self.connect('Fl_O:tot:S', 'out_stat.S') + self.connect('Fl_O:tot:h', 'out_stat.ht') + self.connect('Fl_O:tot:P', 'out_stat.guess:Pt') + self.connect('Fl_O:tot:gamma', 'out_stat.guess:gamt') + self.connect('Wout', 'out_stat.W') + + else: + self.add_subsystem('W_passthru', + PassThrough('Wout', 'Fl_O:stat:W', 1.0, units='lbm/s'), + promotes=['*']) + + super().setup() + + +def print_rde(prob, element_names, file=sys.stdout): + """ + Convenience viewer for RDE Combustor diagnostic properties. + """ + len_header = 23 + 7 * 13 + print("-" * len_header, file=file, flush=True) + print(" RDE COMBUSTOR PROPERTIES", file=file, flush=True) + print("-" * len_header, file=file, flush=True) + + line_tmpl = '{:<20}| ' + '{:>13}' * 7 + print(line_tmpl.format('RDE Element', 'PtIn(psi)', 'PtInj(psi)', 'PtOut(psi)', 'PR_RDE', 'TtOut(degR)', 'D_cj(m/s)', 'f_rde(Hz)'), + file=file, flush=True) + + line_tmpl_data = '{:<20}| {:13.2f}{:13.2f}{:13.2f}{:13.4f}{:13.2f}{:13.1f}{:13.1f}' + for e_name in element_names: + pt_in = prob.get_val(f'{e_name}.rde_press_gain.Pt_in', units='psi')[0] + pt_inj = prob.get_val(f'{e_name}.Pt_inj', units='psi')[0] + pt_out = prob.get_val(f'{e_name}.Fl_O:tot:P', units='psi')[0] + pr_rde = prob.get_val(f'{e_name}.PR_RDE')[0] + tt_out = prob.get_val(f'{e_name}.Fl_O:tot:T', units='degR')[0] + d_cj = prob.get_val(f'{e_name}.D_cj', units='ft/s')[0] * 0.3048 + f_rde = prob.get_val(f'{e_name}.f_rde', units='Hz')[0] + + print(line_tmpl_data.format(e_name, pt_in, pt_inj, pt_out, pr_rde, tt_out, d_cj, f_rde), + file=file, flush=True) + print("-" * len_header, file=file, flush=True) diff --git a/pycycle/elements/test/test_rde_combustor.py b/pycycle/elements/test/test_rde_combustor.py new file mode 100644 index 00000000..15ec4a5e --- /dev/null +++ b/pycycle/elements/test/test_rde_combustor.py @@ -0,0 +1,133 @@ +""" +Tests for RDECombustor element in pyCycle. + +Verifies: +1. Clean integration into pyCycle's Cycle class and flow port network. +2. Pressure gain (Pt_out > Pt_in) with physical fuel addition and chemical equilibrium. +3. Compatibility with both CEA and TABULAR thermodynamics packages. +4. On-design (MN specified) and Off-design (area specified) static property solves. +""" + +import unittest +import numpy as np +import openmdao.api as om +from openmdao.utils.assert_utils import assert_near_equal + +import pycycle.api as pyc +from pycycle.elements.rde_combustor import RDECombustor, RDEPressureGainComp + + +class TestRDECombustor(unittest.TestCase): + + def test_rde_combustor_cea(self): + """Test RDECombustor within a Cycle using CEA thermodynamics.""" + prob = om.Problem() + model = prob.model = pyc.Cycle() + model.options['thermo_method'] = 'CEA' + model.options['thermo_data'] = pyc.species_data.janaf + + model.add_subsystem('flow_start', pyc.FlowStart()) + model.add_subsystem('rde', RDECombustor(fuel_type="JP-7")) + + model.pyc_connect_flow('flow_start.Fl_O', 'rde.Fl_I') + + model.set_input_defaults('rde.Fl_I:FAR', 0.02) + model.set_input_defaults('rde.MN', 0.4) + model.set_input_defaults('rde.dPqP_inj', 0.10) + model.set_input_defaults('rde.eta_rde', 0.85) + + prob.set_solver_print(level=-1) + prob.setup(check=False) + + # Inflow conditions: Station 3 representative + prob.set_val('flow_start.P', 150.0, units='psi') + prob.set_val('flow_start.T', 1100.0, units='degR') + prob.set_val('flow_start.W', 100.0, units='lbm/s') + + prob.run_model() + + pt_in = prob.get_val('flow_start.Fl_O:tot:P', units='psi')[0] + pt_out = prob.get_val('rde.Fl_O:tot:P', units='psi')[0] + pr_rde = prob.get_val('rde.PR_RDE')[0] + w_fuel = prob.get_val('rde.Wfuel', units='lbm/s')[0] + d_cj = prob.get_val('rde.D_cj', units='ft/s')[0] + f_rde = prob.get_val('rde.f_rde', units='Hz')[0] + + # 1. Verify net pressure gain + self.assertGreater(pt_out, pt_in, "RDE exit total pressure must be greater than inlet total pressure") + self.assertGreater(pr_rde, 1.0, "Net pressure ratio must exceed 1.0") + + # 2. Verify fuel mass balance: Wfuel = W_air * FAR = 100 * 0.02 = 2.0 + assert_near_equal(w_fuel, 2.0, tolerance=1e-4) + + # 3. Verify wave kinematics are in physical ranges + self.assertGreater(d_cj, 4000.0, "Detonation velocity must be > 4000 ft/s (~1200 m/s)") + self.assertLess(d_cj, 8000.0, "Detonation velocity must be < 8000 ft/s (~2400 m/s)") + self.assertGreater(f_rde, 1000.0, "RDE frequency must be > 1 kHz") + + def test_rde_combustor_tabular(self): + """Test RDECombustor within a Cycle using fast TABULAR thermodynamics.""" + prob = om.Problem() + model = prob.model = pyc.Cycle() + model.options['thermo_method'] = 'TABULAR' + model.options['thermo_data'] = pyc.AIR_JETA_TAB_SPEC + + model.add_subsystem('flow_start', pyc.FlowStart()) + model.add_subsystem('rde', RDECombustor(fuel_type="FAR")) + + model.pyc_connect_flow('flow_start.Fl_O', 'rde.Fl_I') + + model.set_input_defaults('rde.Fl_I:FAR', 0.025) + model.set_input_defaults('rde.MN', 0.35) + model.set_input_defaults('rde.dPqP_inj', 0.12) + model.set_input_defaults('rde.eta_rde', 0.85) + + prob.set_solver_print(level=-1) + prob.setup(check=False) + + prob.set_val('flow_start.P', 200.0, units='psi') + prob.set_val('flow_start.T', 1200.0, units='degR') + prob.set_val('flow_start.W', 150.0, units='lbm/s') + + prob.run_model() + + pt_in = prob.get_val('flow_start.Fl_O:tot:P', units='psi')[0] + pt_out = prob.get_val('rde.Fl_O:tot:P', units='psi')[0] + pr_rde = prob.get_val('rde.PR_RDE')[0] + + self.assertGreater(pt_out, pt_in, "RDE exit pressure must exceed inlet pressure in tabular mode") + self.assertGreater(pr_rde, 1.0) + + def test_rde_combustor_off_design(self): + """Test RDECombustor in off-design mode with area specified.""" + prob = om.Problem() + model = prob.model = pyc.Cycle(design=False) + model.options['thermo_method'] = 'TABULAR' + model.options['thermo_data'] = pyc.AIR_JETA_TAB_SPEC + + model.add_subsystem('flow_start', pyc.FlowStart()) + model.add_subsystem('rde', RDECombustor(fuel_type="FAR")) + + model.pyc_connect_flow('flow_start.Fl_O', 'rde.Fl_I') + + model.set_input_defaults('rde.Fl_I:FAR', 0.02) + model.set_input_defaults('rde.area', 50.0, units='inch**2') + model.set_input_defaults('rde.dPqP_inj', 0.12) + model.set_input_defaults('rde.eta_rde', 0.85) + + prob.set_solver_print(level=-1) + prob.setup(check=False) + + prob.set_val('flow_start.P', 180.0, units='psi') + prob.set_val('flow_start.T', 1150.0, units='degR') + prob.set_val('flow_start.W', 120.0, units='lbm/s') + + prob.run_model() + + # Check that exit area matches specified input + area_out = prob.get_val('rde.Fl_O:stat:area', units='inch**2')[0] + assert_near_equal(area_out, 50.0, tolerance=1e-4) + + +if __name__ == "__main__": + unittest.main() diff --git a/uv.lock b/uv.lock new file mode 100644 index 00000000..7518fc90 --- /dev/null +++ b/uv.lock @@ -0,0 +1,3 @@ +version = 1 +revision = 3 +requires-python = ">=3.12"