diff --git a/ceasiompy/ThermoData/README.md b/ceasiompy/ThermoData/README.md new file mode 100644 index 000000000..77e55d5ff --- /dev/null +++ b/ceasiompy/ThermoData/README.md @@ -0,0 +1 @@ +# Thermo_Data \ No newline at end of file diff --git a/ceasiompy/ThermoData/__specs__.py b/ceasiompy/ThermoData/__specs__.py new file mode 100644 index 000000000..86c68fb89 --- /dev/null +++ b/ceasiompy/ThermoData/__specs__.py @@ -0,0 +1,69 @@ +from pathlib import Path +from ceasiompy.utils.moduleinterfaces import CPACSInOut +from ceasiompy.utils.commonxpath import ( + REF_XPATH, + CLCALC_XPATH, + SU2_FIXED_CL_XPATH, + SU2_TARGET_CL_XPATH, + ENGINE_TYPE_XPATH, + RANGE_XPATH, +) + +# ===== Module Status ===== +# True if the module is active +# False if the module is disabled (not working or not ready) +module_status = True + +# ===== Results directory path ===== + +RESULTS_DIR = Path("Results", "Thermodata") + +# ===== CPACS inputs and outputs ===== + +cpacs_inout = CPACSInOut() + +# ===== Input ===== + + +cpacs_inout.add_input( + var_name="net_force", + var_type=float, + default_value=2000, + unit="1", # AGGIUNGERE UNITA DI MISURA + descr="Engine net force", + xpath=RANGE_XPATH + "/NetForce", + gui=True, + gui_name="NetForce", + gui_group="Cruise", +) + +cpacs_inout.add_input( + var_name="engine type", + var_type=list, + default_value=[0, 1], + unit=None, + descr="0: TBJ, 1: TBF ", + xpath=ENGINE_TYPE_XPATH, + gui=True, + gui_name="0 for Turbojet 1 for Turbofan", + gui_group="User inputs", +) + + +# ===== Output ===== + +cpacs_inout.add_output( + var_name="target_cl", + default_value=None, + unit="1", + descr="Value of CL to achieve to have a level flight with the given conditions", + xpath=SU2_TARGET_CL_XPATH, +) + +cpacs_inout.add_output( + var_name="fixed_cl", + default_value=None, + unit="-", + descr="FIXED_CL_MODE parameter for SU2", + xpath=SU2_FIXED_CL_XPATH, +) diff --git a/ceasiompy/ThermoData/func/run_func.py b/ceasiompy/ThermoData/func/run_func.py new file mode 100644 index 000000000..eb0a6b58a --- /dev/null +++ b/ceasiompy/ThermoData/func/run_func.py @@ -0,0 +1,162 @@ +""" +CEASIOMpy: Conceptual Aircraft Design Software + +Developed by CFS ENGINEERING, 1015 Lausanne, Switzerland + +Function to run the PyCycle code in order to obtain turbojet and turbofan boundary conditions +at given altitude, mach number and net force = thrust of the engine + +Python version: >=3.8 + +| Author: Giacomo Benedetti and Francesco Marcucci +| Creation: 2023-12-12 + +""" + +from pathlib import Path + +from ceasiompy.ThermoData.func.turbofan_func import ( + turbofan_analysis, + write_hbtf_file, +) + +from ceasiompy.ThermoData.func.turbojet_func import ( + turbojet_analysis, + write_turbojet_file, +) + +from ceasiompy.utils.ceasiomlogger import get_logger + +from ceasiompy.utils.commonxpath import ( + ENGINE_TYPE_XPATH, + ENGINE_BC, + RANGE_XPATH, + SU2_AEROMAP_UID_XPATH, +) +from ceasiompy.utils.commonnames import ( + ENGINE_BOUNDARY_CONDITIONS, +) +from cpacspy.cpacsfunctions import ( + get_value_or_default, + add_float_vector, +) +from cpacspy.cpacsfunctions import create_branch, get_value_or_default +from cpacspy.cpacspy import CPACS + +log = get_logger() + +MODULE_DIR = Path(__file__).parent +MODULE_NAME = MODULE_DIR.name + + +def thermo_data_run(cpacs_path, cpacs_out_path, wkdir): + """Running the PyCycle code by choosing between turbojet or turbofan engine + + Args + cpacs_path (Path): Path to CPACS file + cpacs_out_path (Path):Path to CPACS output file + wkdir (str): Path to the working directory + """ + + if not wkdir.exists(): + raise OSError(f"The working directory : {wkdir} does not exit!") + + cpacs = CPACS(cpacs_path) + tixi = cpacs.tixi + + Fn = get_value_or_default(tixi, RANGE_XPATH + "/NetForce", 2000) + + aeromap_list = cpacs.get_aeromap_uid_list() + + if aeromap_list: + aeromap_default = aeromap_list[0] + log.info(f"The aeromap is {aeromap_default}") + aeromap_uid = get_value_or_default( + cpacs.tixi, SU2_AEROMAP_UID_XPATH, aeromap_default + ) + activate_aeromap = cpacs.get_aeromap_by_uid(aeromap_uid) + alt_list = activate_aeromap.get("altitude").tolist() + mach_list = activate_aeromap.get("machNumber").tolist() + T_tot_out_array = [] + P_tot_out_array = [] + + for case_nb in range(len(alt_list)): + alt = alt_list[case_nb] + MN = mach_list[case_nb] + case_dir_name = f"Case{str(case_nb).zfill(2)}_alt{alt}_mach{round(MN, 2)}" + case_dir_path = Path(wkdir, case_dir_name) + + if not case_dir_path.exists(): + case_dir_path.mkdir() + + EngineBC = Path(case_dir_path, ENGINE_BOUNDARY_CONDITIONS) + + f = open(EngineBC, "w") + + engine_type = get_value_or_default(tixi, ENGINE_TYPE_XPATH, 0) + create_branch(cpacs.tixi, ENGINE_BC) + + if engine_type == 0: + ( + T_tot_out, + V_stat_out, + MN_out, + P_tot_out, + massflow_stat_out, + T_stat_out, + P_stat_out, + ) = turbojet_analysis(alt, MN, Fn) + + T_tot_out_array.append(T_tot_out) + P_tot_out_array.append(P_tot_out) + + f = write_turbojet_file( + file=f, + T_tot_out=T_tot_out, + V_stat_out=V_stat_out, + MN_out=MN_out, + P_tot_out=P_tot_out, + massflow_stat_out=massflow_stat_out, + T_stat_out=T_stat_out, + P_stat_out=P_stat_out, + ) + + else: + + ( + T_tot_out_byp, + V_stat_out_byp, + MN_out_byp, + P_tot_out_byp, + massflow_stat_out_byp, + T_stat_out_byp, + T_tot_out_core, + V_stat_out_core, + MN_out_core, + P_tot_out_core, + massflow_stat_out_core, + T_stat_out_core, + ) = turbofan_analysis(alt, MN, Fn) + + T_tot_out_array.append(T_tot_out_core) + P_tot_out_array.append(P_tot_out_core) + + f = write_hbtf_file( + file=f, + T_tot_out_byp=T_tot_out_byp, + V_stat_out_byp=V_stat_out_byp, + MN_out_byp=MN_out_byp, + P_tot_out_byp=P_tot_out_byp, + massflow_stat_out_byp=massflow_stat_out_byp, + T_stat_out_byp=T_stat_out_byp, + T_tot_out_core=T_tot_out_core, + V_stat_out_core=V_stat_out_core, + MN_out_core=MN_out_core, + P_tot_out_core=P_tot_out_core, + massflow_stat_out_core=massflow_stat_out_core, + T_stat_out_core=T_stat_out_core, + ) + add_float_vector(tixi, ENGINE_BC + "/temperatureOutlet", T_tot_out_array) + add_float_vector(tixi, ENGINE_BC + "/pressureOutlet", P_tot_out_array) + + cpacs.save_cpacs(cpacs_out_path, overwrite=True) diff --git a/ceasiompy/ThermoData/func/turbofan_func.py b/ceasiompy/ThermoData/func/turbofan_func.py new file mode 100644 index 000000000..52521c47f --- /dev/null +++ b/ceasiompy/ThermoData/func/turbofan_func.py @@ -0,0 +1,436 @@ +""" +CEASIOMpy: Conceptual Aircraft Design Software + +Developed by CFS ENGINEERING, 1015 Lausanne, Switzerland + +Function to run the PyCycle code for the turbofan engine + +Python version: >=3.8 + +| Author: Giacomo Benedetti and Francesco Marcucci +| Creation: 2023-12-12 + +""" +import sys + +import openmdao.api as om + +import pycycle.api as pyc + +from scipy.constants import convert_temperature + +from ceasiompy.utils.ceasiomlogger import get_logger + +log = get_logger() + +# ================================================================================================= +# FUNCTIONS +# ================================================================================================= + + +def turbofan_analysis(alt, MN, Fn): + class HBTF(pyc.Cycle): + def setup(self): + + # Setup the problem by including all the relevant components here - comp, burner, turbine etc + + # Create any relevant short hands here: + + design = self.options["design"] + + USE_TABULAR = False + if USE_TABULAR: + self.options["thermo_method"] = "TABULAR" + self.options["thermo_data"] = pyc.AIR_JETA_TAB_SPEC + FUEL_TYPE = "FAR" + else: + self.options["thermo_method"] = "CEA" + self.options["thermo_data"] = pyc.species_data.janaf + FUEL_TYPE = "Jet-A(g)" + + # Add subsystems to build the engine deck: + self.add_subsystem("fc", pyc.FlightConditions()) + self.add_subsystem("inlet", pyc.Inlet()) + + self.add_subsystem( + "fan", + pyc.Compressor(map_data=pyc.FanMap, bleed_names=[], map_extrap=True), + promotes_inputs=[("Nmech", "LP_Nmech")], + ) + self.add_subsystem("splitter", pyc.Splitter()) + self.add_subsystem("duct4", pyc.Duct()) + self.add_subsystem( + "lpc", + pyc.Compressor(map_data=pyc.LPCMap, map_extrap=True), + promotes_inputs=[("Nmech", "LP_Nmech")], + ) + self.add_subsystem("duct6", pyc.Duct()) + self.add_subsystem( + "hpc", + pyc.Compressor( + map_data=pyc.HPCMap, + bleed_names=["cool1", "cool2", "cust"], + map_extrap=True, + ), + promotes_inputs=[("Nmech", "HP_Nmech")], + ) + self.add_subsystem("bld3", pyc.BleedOut(bleed_names=["cool3", "cool4"])) + self.add_subsystem("burner", pyc.Combustor(fuel_type=FUEL_TYPE)) + self.add_subsystem( + "hpt", + pyc.Turbine( + map_data=pyc.HPTMap, bleed_names=["cool3", "cool4"], map_extrap=True + ), + promotes_inputs=[("Nmech", "HP_Nmech")], + ) + self.add_subsystem("duct11", pyc.Duct()) + self.add_subsystem( + "lpt", + pyc.Turbine( + map_data=pyc.LPTMap, bleed_names=["cool1", "cool2"], map_extrap=True + ), + promotes_inputs=[("Nmech", "LP_Nmech")], + ) + self.add_subsystem("duct13", pyc.Duct()) + self.add_subsystem("core_nozz", pyc.Nozzle(nozzType="CV", lossCoef="Cv")) + + self.add_subsystem("byp_bld", pyc.BleedOut(bleed_names=["bypBld"])) + self.add_subsystem("duct15", pyc.Duct()) + self.add_subsystem("byp_nozz", pyc.Nozzle(nozzType="CV", lossCoef="Cv")) + + # Create shaft instances. Note that LP shaft has 3 ports! => no gearbox + self.add_subsystem( + "lp_shaft", + pyc.Shaft(num_ports=3), + promotes_inputs=[("Nmech", "LP_Nmech")], + ) + self.add_subsystem( + "hp_shaft", + pyc.Shaft(num_ports=2), + promotes_inputs=[("Nmech", "HP_Nmech")], + ) + self.add_subsystem("perf", pyc.Performance(num_nozzles=2, num_burners=1)) + + # Now use the explicit connect method to make connections -- connect(, ) + + # Connect the inputs to perf group + self.connect("inlet.Fl_O:tot:P", "perf.Pt2") + self.connect("hpc.Fl_O:tot:P", "perf.Pt3") + self.connect("burner.Wfuel", "perf.Wfuel_0") + self.connect("inlet.F_ram", "perf.ram_drag") + self.connect("core_nozz.Fg", "perf.Fg_0") + self.connect("byp_nozz.Fg", "perf.Fg_1") + + # LP-shaft connections + self.connect("fan.trq", "lp_shaft.trq_0") + self.connect("lpc.trq", "lp_shaft.trq_1") + self.connect("lpt.trq", "lp_shaft.trq_2") + # HP-shaft connections + self.connect("hpc.trq", "hp_shaft.trq_0") + self.connect("hpt.trq", "hp_shaft.trq_1") + # Ideally expanding flow by conneting flight condition static pressure to nozzle exhaust pressure + self.connect("fc.Fl_O:stat:P", "core_nozz.Ps_exhaust") + self.connect("fc.Fl_O:stat:P", "byp_nozz.Ps_exhaust") + + balance = self.add_subsystem("balance", om.BalanceComp()) + if design: + balance.add_balance("W", units="lbm/s", eq_units="lbf") + # Here balance.W is implicit state variable that is the OUTPUT of balance object + self.connect( + "balance.W", "fc.W" + ) # Connect the output of balance to the relevant input + self.connect( + "perf.Fn", "balance.lhs:W" + ) # This statement makes perf.Fn the LHS of the balance eqn. + self.promotes("balance", inputs=[("rhs:W", "Fn_DES")]) + + balance.add_balance("FAR", eq_units="degR", lower=1e-4, val=0.017) + self.connect("balance.FAR", "burner.Fl_I:FAR") + self.connect("burner.Fl_O:tot:T", "balance.lhs:FAR") + self.promotes("balance", inputs=[("rhs:FAR", "T4_MAX")]) + + # Note that for the following two balances the mult val is set to -1 so that the NET torque is zero + balance.add_balance( + "lpt_PR", + val=1.5, + lower=1.001, + upper=8, + eq_units="hp", + use_mult=True, + mult_val=-1, + ) + self.connect("balance.lpt_PR", "lpt.PR") + self.connect("lp_shaft.pwr_in_real", "balance.lhs:lpt_PR") + self.connect("lp_shaft.pwr_out_real", "balance.rhs:lpt_PR") + + balance.add_balance( + "hpt_PR", + val=1.5, + lower=1.001, + upper=8, + eq_units="hp", + use_mult=True, + mult_val=-1, + ) + self.connect("balance.hpt_PR", "hpt.PR") + self.connect("hp_shaft.pwr_in_real", "balance.lhs:hpt_PR") + self.connect("hp_shaft.pwr_out_real", "balance.rhs:hpt_PR") + + # Set up all the flow connections: + self.pyc_connect_flow("fc.Fl_O", "inlet.Fl_I") + self.pyc_connect_flow("inlet.Fl_O", "fan.Fl_I") + self.pyc_connect_flow("fan.Fl_O", "splitter.Fl_I") + self.pyc_connect_flow("splitter.Fl_O1", "duct4.Fl_I") + self.pyc_connect_flow("duct4.Fl_O", "lpc.Fl_I") + self.pyc_connect_flow("lpc.Fl_O", "duct6.Fl_I") + self.pyc_connect_flow("duct6.Fl_O", "hpc.Fl_I") + self.pyc_connect_flow("hpc.Fl_O", "bld3.Fl_I") + self.pyc_connect_flow("bld3.Fl_O", "burner.Fl_I") + self.pyc_connect_flow("burner.Fl_O", "hpt.Fl_I") + self.pyc_connect_flow("hpt.Fl_O", "duct11.Fl_I") + self.pyc_connect_flow("duct11.Fl_O", "lpt.Fl_I") + self.pyc_connect_flow("lpt.Fl_O", "duct13.Fl_I") + self.pyc_connect_flow("duct13.Fl_O", "core_nozz.Fl_I") + self.pyc_connect_flow("splitter.Fl_O2", "byp_bld.Fl_I") + self.pyc_connect_flow("byp_bld.Fl_O", "duct15.Fl_I") + self.pyc_connect_flow("duct15.Fl_O", "byp_nozz.Fl_I") + + # Bleed flows: + self.pyc_connect_flow("hpc.cool1", "lpt.cool1", connect_stat=False) + self.pyc_connect_flow("hpc.cool2", "lpt.cool2", connect_stat=False) + self.pyc_connect_flow("bld3.cool3", "hpt.cool3", connect_stat=False) + self.pyc_connect_flow("bld3.cool4", "hpt.cool4", connect_stat=False) + + # Specify solver settings: + newton = self.nonlinear_solver = om.NewtonSolver() + newton.options["atol"] = 1e-8 + + # set this very small, so it never activates and we rely on atol + newton.options["rtol"] = 1e-99 + newton.options["iprint"] = 2 + newton.options["maxiter"] = 50 + newton.options["solve_subsystems"] = True + newton.options["max_sub_solves"] = 1000 + newton.options["reraise_child_analysiserror"] = False + # ls = newton.linesearch = BoundsEnforceLS() + ls = newton.linesearch = om.ArmijoGoldsteinLS() + ls.options["maxiter"] = 3 + ls.options["rho"] = 0.75 + # ls.options['print_bound_enforce'] = True + + self.linear_solver = om.DirectSolver() + + super().setup() + + class MPhbtf(pyc.MPCycle): + def setup(self): + + self.pyc_add_pnt( + "DESIGN", HBTF(thermo_method="CEA") + ) # Create an instance of the High Bypass ratio Turbofan + + super().setup() + + def viewer(prob, pt, file=sys.stdout): + """ + print a report of all the relevant cycle properties + """ + + if pt == "DESIGN": + MN = prob["DESIGN.fc.Fl_O:stat:MN"] + LPT_PR = prob["DESIGN.balance.lpt_PR"] + HPT_PR = prob["DESIGN.balance.hpt_PR"] + FAR = prob["DESIGN.balance.FAR"] + else: + MN = prob[pt + ".fc.Fl_O:stat:MN"] + LPT_PR = prob[pt + ".lpt.PR"] + HPT_PR = prob[pt + ".hpt.PR"] + FAR = prob[pt + ".balance.FAR"] + + summary_data = ( + 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"], + prob[pt + ".splitter.BPR"], + ) + + 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 BPR ", + file=file, + flush=True, + ) + print( + " %7.5f %7.1f %7.3f %7.1f %7.1f %7.1f %7.3f %7.5f %7.3f" % summary_data, + file=file, + flush=True, + ) + + fs_names = [ + "fc.Fl_O", + "core_nozz.Fl_O", + "byp_nozz.Fl_O", + ] + fs_full_names = [f"{pt}.{fs}" for fs in fs_names] + pyc.print_flow_station(prob, fs_full_names, file=file) + + import time + + prob = om.Problem() + + prob.model = mp_hbtf = MPhbtf() + + prob.setup() + + prob.set_val("DESIGN.fan.PR", 1.685) + prob.set_val("DESIGN.fan.eff", 0.8948) + + prob.set_val("DESIGN.lpc.PR", 1.935) + prob.set_val("DESIGN.lpc.eff", 0.9243) + + prob.set_val("DESIGN.hpc.PR", 9.369) + prob.set_val("DESIGN.hpc.eff", 0.8707) + + prob.set_val("DESIGN.hpt.eff", 0.8888) + prob.set_val("DESIGN.lpt.eff", 0.8996) + + prob.set_val("DESIGN.fc.alt", alt, units="ft") + prob.set_val("DESIGN.fc.MN", MN) + + prob.set_val("DESIGN.T4_MAX", 2857, units="degR") + prob.set_val("DESIGN.Fn_DES", Fn, units="lbf") + + # Set initial guesses for balances + prob["DESIGN.balance.FAR"] = 0.1 + prob["DESIGN.balance.W"] = 10.0 + prob["DESIGN.balance.lpt_PR"] = 2 + prob["DESIGN.balance.hpt_PR"] = 2.0 + prob["DESIGN.fc.balance.Pt"] = 2 + prob["DESIGN.fc.balance.Tt"] = 500.0 + + prob.set_solver_print(level=-1) + prob.set_solver_print(level=2, depth=1) + + flight_env = [(0.9, 6000)] + + viewer_file = open("hbtf_des_view.out", "w") + + for PC in [1, 0.9]: + viewer(prob, "DESIGN", file=viewer_file) + prob.run_model() + + print() + + # Obtaining variables names + + # variable_names = open("variables_hptf.txt", "w") + # print( + # prob.model.list_outputs(val=True, units=True, implicit=False), + # file=variable_names, + # ) + + # BYPASS VARIABLES + T_tot_out_byp = convert_temperature( + (prob.get_val("DESIGN.byp_nozz.throat_total.flow.Fl_O:tot:T")), + "Rankine", + "Celsius", + ) + V_stat_out_byp = prob.get_val("DESIGN.byp_nozz.mux.Fl_O:stat:V") * 0.3048 + MN_out_byp = prob.get_val("DESIGN.byp_nozz.mux.Fl_O:stat:MN") + P_tot_out_byp = ( + prob.get_val("DESIGN.byp_nozz.throat_total.flow.Fl_O:tot:P") * 6894.7573 + ) # Pa + massflow_stat_out_byp = ( + prob.get_val("DESIGN.byp_nozz.mux.Fl_O:stat:W") * 0.45359237 + ) # kg/s + T_stat_out_byp = convert_temperature( + (prob.get_val("DESIGN.byp_nozz.mux.Fl_O:stat:T")), "Rankine", "Celsius" + ) # celsius + + # CORE VARIABLES + T_tot_out_core = convert_temperature( + (prob.get_val("DESIGN.core_nozz.throat_total.flow.Fl_O:tot:T")), + "Rankine", + "Kelvin", + ) + V_stat_out_core = prob.get_val("DESIGN.core_nozz.mux.Fl_O:stat:V") * 0.3048 + MN_out_core = prob.get_val("DESIGN.core_nozz.mux.Fl_O:stat:MN") + P_tot_out_core = ( + prob.get_val("DESIGN.core_nozz.throat_total.flow.Fl_O:tot:P") * 6894.7573 + ) # Pa + massflow_stat_out_core = ( + prob.get_val("DESIGN.core_nozz.mux.Fl_O:stat:W") * 0.45359237 + ) # kg/s + T_stat_out_core = convert_temperature( + (prob.get_val("DESIGN.core_nozz.mux.Fl_O:stat:T")), "Rankine", "Kelvin" + ) # celsius + + return ( + T_tot_out_byp, + V_stat_out_byp, + MN_out_byp, + P_tot_out_byp, + massflow_stat_out_byp, + T_stat_out_byp, + T_tot_out_core, + V_stat_out_core, + MN_out_core, + P_tot_out_core, + massflow_stat_out_core, + T_stat_out_core, + ) + + +def write_hbtf_file( + file, + T_tot_out_byp, + V_stat_out_byp, + MN_out_byp, + P_tot_out_byp, + massflow_stat_out_byp, + T_stat_out_byp, + T_tot_out_core, + V_stat_out_core, + MN_out_core, + P_tot_out_core, + massflow_stat_out_core, + T_stat_out_core, +): + + file.write(f"T_tot_out_core = {T_tot_out_core} [K]\n") + file.write(f"V_stat_out_core = {V_stat_out_core} [m/s]\n") + file.write(f"MN_out_core = {MN_out_core} [adim]\n") + file.write(f"P_tot_out_core = {P_tot_out_core} [Pa]\n") + file.write(f"massflow_out_core = {massflow_stat_out_core} [kg/s]\n") + file.write(f"T_stat_out_core = {T_stat_out_core} [K]\n") + file.write(f"T_tot_out_byp = {T_tot_out_byp} [K]\n") + file.write(f"V_stat_out_byp = {V_stat_out_byp} [m/s]\n") + file.write(f"MN_out_byp = {MN_out_byp} [adim]\n") + file.write(f"P_tot_out_byp = {P_tot_out_byp} [Pa]\n") + file.write(f"massflow_stat_out_byp = {massflow_stat_out_byp} [kg/s]\n") + file.write(f"T_stat_out_byp = {T_stat_out_byp} [K]\n") + + log.info("hbtf.dat file generated!") + + return file diff --git a/ceasiompy/ThermoData/func/turbojet_func.py b/ceasiompy/ThermoData/func/turbojet_func.py new file mode 100644 index 000000000..371195bad --- /dev/null +++ b/ceasiompy/ThermoData/func/turbojet_func.py @@ -0,0 +1,228 @@ +""" +CEASIOMpy: Conceptual Aircraft Design Software + +Developed by CFS ENGINEERING, 1015 Lausanne, Switzerland + +Function to run the PyCycle code for the turbojet engine + +Python version: >=3.8 + +| Author: Giacomo Benedetti and Francesco Marcucci +| Creation: 2023-12-12 + +""" + +import openmdao.api as om + +import pycycle.api as pyc + +from scipy.constants import convert_temperature + +from ceasiompy.utils.ceasiomlogger import get_logger + +log = get_logger() + +# ================================================================================================= +# FUNCTIONS +# ================================================================================================= + + +def turbojet_analysis(alt, MN, Fn): + class Turbojet(pyc.Cycle): + def setup(self): + + USE_TABULAR = True + + if USE_TABULAR: + self.options["thermo_method"] = "TABULAR" + self.options["thermo_data"] = pyc.AIR_JETA_TAB_SPEC + FUEL_TYPE = "FAR" + + design = self.options["design"] + + # Add engine elements + 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"], + ) + 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)) + + # 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") + + # Make other non-flow connections + # Connect turbomachinery elements to shaft + self.connect("comp.trq", "shaft.trq_0") + self.connect("turb.trq", "shaft.trq_1") + + # Connnect nozzle exhaust to freestream static conditions + self.connect("fc.Fl_O:stat:P", "nozz.Ps_exhaust") + + # Connect outputs to perfomance element + 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") + + # Add balances for design and off-design + balance = self.add_subsystem("balance", om.BalanceComp()) + if design: + + 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") + + balance.add_balance( + "FAR", eq_units="degR", lower=1e-4, val=0.017, rhs_name="T4_target" + ) + self.connect("balance.FAR", "burner.Fl_I:FAR") + self.connect("burner.Fl_O:tot:T", "balance.lhs:FAR") + + balance.add_balance( + "turb_PR", val=1.5, lower=1.001, upper=8, eq_units="hp", rhs_val=0.0 + ) + self.connect("balance.turb_PR", "turb.PR") + self.connect("shaft.pwr_net", "balance.lhs:turb_PR") + + newton = self.nonlinear_solver = om.NewtonSolver() + newton.options["atol"] = 1e-6 + newton.options["rtol"] = 1e-6 + newton.options["iprint"] = 2 + newton.options["maxiter"] = 15 + 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 MPTurbojet(pyc.MPCycle): + def setup(self): + self.pyc_add_pnt("DESIGN", Turbojet()) + + 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) # .2 + self.set_input_defaults("DESIGN.burner.MN", 0.020) # .2 + self.set_input_defaults("DESIGN.turb.MN", 0.4) + + self.pyc_add_cycle_param("burner.dPqP", 0.03) + self.pyc_add_cycle_param("nozz.Cv", 0.99) + + super().setup() + + import time + + # defining the optimization problem + + prob = om.Problem() + + mp_turbojet = prob.model = MPTurbojet() + prob.setup(check=False) + + # define input values + prob.set_val("DESIGN.fc.alt", alt, units="ft") + prob.set_val("DESIGN.fc.MN", MN) + prob.set_val("DESIGN.balance.Fn_target", Fn, 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) + + if Fn > 2700: + prob["DESIGN.balance.FAR"] = 0.0175506829934 + prob["DESIGN.balance.W"] = 75.453135137 + prob["DESIGN.balance.turb_PR"] = 4.46138725662 + prob["DESIGN.fc.balance.Pt"] = 14.6955113159 + prob["DESIGN.fc.balance.Tt"] = 518.665288153 + + elif 2700 <= Fn < 850: + prob["DESIGN.balance.FAR"] = 0.0175506829934 + prob["DESIGN.balance.W"] = 50.453135137 + prob["DESIGN.balance.turb_PR"] = 4.46138725662 + prob["DESIGN.fc.balance.Pt"] = 14.6955113159 + prob["DESIGN.fc.balance.Tt"] = 518.665288153 + + else: + prob["DESIGN.balance.FAR"] = 0.0175506829934 + prob["DESIGN.balance.W"] = 10.453135137 + prob["DESIGN.balance.turb_PR"] = 4.46138725662 + prob["DESIGN.fc.balance.Pt"] = 14.6955113159 + prob["DESIGN.fc.balance.Tt"] = 518.665288153 + + prob.set_solver_print(level=-1) + # prob.set_solver_print(level=2, depth=1) + + prob.run_model() + + print() + + # command to visualize all the output of the system + # a = prob.model.list_outputs(val=True, units=True) + + T_tot_out = convert_temperature( + (prob.get_val("DESIGN.nozz.throat_total.flow.Fl_O:tot:T")), "Rankine", "Kelvin" + ) + V_stat_out = (prob.get_val("DESIGN.nozz.mux.Fl_O:stat:V")) * 0.3048 + MN_out = prob.get_val("DESIGN.nozz.mux.Fl_O:stat:MN") + T_stat_out = convert_temperature( + prob.get_val("DESIGN.nozz.mux.Fl_O:stat:T"), "Rankine", "Kelvin" + ) # celsius + massflow_stat_out = prob.get_val("DESIGN.nozz.mux.Fl_O:stat:W") * 0.45359237 # kg/s + P_tot_out = ( + prob.get_val("DESIGN.nozz.throat_total.flow.Fl_O:tot:P") * 6894.7573 + ) # Pa + P_stat_out = prob.get_val("DESIGN.nozz.mux.Fl_O:stat:P") * 6894.7573 # Pa + + return ( + T_tot_out, + V_stat_out, + MN_out, + P_tot_out, + massflow_stat_out, + T_stat_out, + P_stat_out, + ) + + +def write_turbojet_file( + file, + T_tot_out, + V_stat_out, + MN_out, + P_tot_out, + massflow_stat_out, + T_stat_out, + P_stat_out, +): + + file.write(f"T_tot_out = {T_tot_out} [K]\n") + file.write(f"V_stat_out = {V_stat_out} [m/s]\n") + file.write(f"MN_out = {MN_out} [adim]\n") + file.write(f"P_tot_out = {P_tot_out} [Pa]\n") + file.write(f"massflow_out = {massflow_stat_out} [kg/s]\n") + file.write(f"T_stat_out = {T_stat_out} [K]\n") + file.write(f"P_stat_out = {P_stat_out} [Pa]\n") + + log.info("turbojet.dat file generated!") + + return file diff --git a/ceasiompy/ThermoData/testing_pycycle/des_hbtf_func_test.py b/ceasiompy/ThermoData/testing_pycycle/des_hbtf_func_test.py new file mode 100644 index 000000000..dbc960816 --- /dev/null +++ b/ceasiompy/ThermoData/testing_pycycle/des_hbtf_func_test.py @@ -0,0 +1,390 @@ +import sys + +import numpy as np + +import openmdao.api as om + +import pycycle.api as pyc + +import re + +from scipy.constants import convert_temperature + + +def run_turbofan_analysis_test(alt, MN): + class HBTF(pyc.Cycle): + def setup(self): + + # Setup the problem by including all the relevant components here - comp, burner, turbine etc + + # Create any relevant short hands here: + + design = self.options["design"] + + USE_TABULAR = False + if USE_TABULAR: + self.options["thermo_method"] = "TABULAR" + self.options["thermo_data"] = pyc.AIR_JETA_TAB_SPEC + FUEL_TYPE = "FAR" + else: + self.options["thermo_method"] = "CEA" + self.options["thermo_data"] = pyc.species_data.janaf + FUEL_TYPE = "Jet-A(g)" + + # Add subsystems to build the engine deck: + self.add_subsystem("fc", pyc.FlightConditions()) + self.add_subsystem("inlet", pyc.Inlet()) + + self.add_subsystem( + "fan", + pyc.Compressor(map_data=pyc.FanMap, bleed_names=[], map_extrap=True), + promotes_inputs=[("Nmech", "LP_Nmech")], + ) + self.add_subsystem("splitter", pyc.Splitter()) + self.add_subsystem("duct4", pyc.Duct()) + self.add_subsystem( + "lpc", + pyc.Compressor(map_data=pyc.LPCMap, map_extrap=True), + promotes_inputs=[("Nmech", "LP_Nmech")], + ) + self.add_subsystem("duct6", pyc.Duct()) + self.add_subsystem( + "hpc", + pyc.Compressor( + map_data=pyc.HPCMap, + bleed_names=["cool1", "cool2", "cust"], + map_extrap=True, + ), + promotes_inputs=[("Nmech", "HP_Nmech")], + ) + self.add_subsystem("bld3", pyc.BleedOut(bleed_names=["cool3", "cool4"])) + self.add_subsystem("burner", pyc.Combustor(fuel_type=FUEL_TYPE)) + self.add_subsystem( + "hpt", + pyc.Turbine( + map_data=pyc.HPTMap, bleed_names=["cool3", "cool4"], map_extrap=True + ), + promotes_inputs=[("Nmech", "HP_Nmech")], + ) + self.add_subsystem("duct11", pyc.Duct()) + self.add_subsystem( + "lpt", + pyc.Turbine( + map_data=pyc.LPTMap, bleed_names=["cool1", "cool2"], map_extrap=True + ), + promotes_inputs=[("Nmech", "LP_Nmech")], + ) + self.add_subsystem("duct13", pyc.Duct()) + self.add_subsystem("core_nozz", pyc.Nozzle(nozzType="CV", lossCoef="Cv")) + + self.add_subsystem("byp_bld", pyc.BleedOut(bleed_names=["bypBld"])) + self.add_subsystem("duct15", pyc.Duct()) + self.add_subsystem("byp_nozz", pyc.Nozzle(nozzType="CV", lossCoef="Cv")) + + # Create shaft instances. Note that LP shaft has 3 ports! => no gearbox + self.add_subsystem( + "lp_shaft", + pyc.Shaft(num_ports=3), + promotes_inputs=[("Nmech", "LP_Nmech")], + ) + self.add_subsystem( + "hp_shaft", + pyc.Shaft(num_ports=2), + promotes_inputs=[("Nmech", "HP_Nmech")], + ) + self.add_subsystem("perf", pyc.Performance(num_nozzles=2, num_burners=1)) + + # Now use the explicit connect method to make connections -- connect(, ) + + # Connect the inputs to perf group + self.connect("inlet.Fl_O:tot:P", "perf.Pt2") + self.connect("hpc.Fl_O:tot:P", "perf.Pt3") + self.connect("burner.Wfuel", "perf.Wfuel_0") + self.connect("inlet.F_ram", "perf.ram_drag") + self.connect("core_nozz.Fg", "perf.Fg_0") + self.connect("byp_nozz.Fg", "perf.Fg_1") + + # LP-shaft connections + self.connect("fan.trq", "lp_shaft.trq_0") + self.connect("lpc.trq", "lp_shaft.trq_1") + self.connect("lpt.trq", "lp_shaft.trq_2") + # HP-shaft connections + self.connect("hpc.trq", "hp_shaft.trq_0") + self.connect("hpt.trq", "hp_shaft.trq_1") + # Ideally expanding flow by conneting flight condition static pressure to nozzle exhaust pressure + self.connect("fc.Fl_O:stat:P", "core_nozz.Ps_exhaust") + self.connect("fc.Fl_O:stat:P", "byp_nozz.Ps_exhaust") + + balance = self.add_subsystem("balance", om.BalanceComp()) + if design: + balance.add_balance("W", units="lbm/s", eq_units="lbf") + # Here balance.W is implicit state variable that is the OUTPUT of balance object + self.connect( + "balance.W", "fc.W" + ) # Connect the output of balance to the relevant input + self.connect( + "perf.Fn", "balance.lhs:W" + ) # This statement makes perf.Fn the LHS of the balance eqn. + self.promotes("balance", inputs=[("rhs:W", "Fn_DES")]) + + balance.add_balance("FAR", eq_units="degR", lower=1e-4, val=0.017) + self.connect("balance.FAR", "burner.Fl_I:FAR") + self.connect("burner.Fl_O:tot:T", "balance.lhs:FAR") + self.promotes("balance", inputs=[("rhs:FAR", "T4_MAX")]) + + # Note that for the following two balances the mult val is set to -1 so that the NET torque is zero + balance.add_balance( + "lpt_PR", + val=1.5, + lower=1.001, + upper=8, + eq_units="hp", + use_mult=True, + mult_val=-1, + ) + self.connect("balance.lpt_PR", "lpt.PR") + self.connect("lp_shaft.pwr_in_real", "balance.lhs:lpt_PR") + self.connect("lp_shaft.pwr_out_real", "balance.rhs:lpt_PR") + + balance.add_balance( + "hpt_PR", + val=1.5, + lower=1.001, + upper=8, + eq_units="hp", + use_mult=True, + mult_val=-1, + ) + self.connect("balance.hpt_PR", "hpt.PR") + self.connect("hp_shaft.pwr_in_real", "balance.lhs:hpt_PR") + self.connect("hp_shaft.pwr_out_real", "balance.rhs:hpt_PR") + + # Set up all the flow connections: + self.pyc_connect_flow("fc.Fl_O", "inlet.Fl_I") + self.pyc_connect_flow("inlet.Fl_O", "fan.Fl_I") + self.pyc_connect_flow("fan.Fl_O", "splitter.Fl_I") + self.pyc_connect_flow("splitter.Fl_O1", "duct4.Fl_I") + self.pyc_connect_flow("duct4.Fl_O", "lpc.Fl_I") + self.pyc_connect_flow("lpc.Fl_O", "duct6.Fl_I") + self.pyc_connect_flow("duct6.Fl_O", "hpc.Fl_I") + self.pyc_connect_flow("hpc.Fl_O", "bld3.Fl_I") + self.pyc_connect_flow("bld3.Fl_O", "burner.Fl_I") + self.pyc_connect_flow("burner.Fl_O", "hpt.Fl_I") + self.pyc_connect_flow("hpt.Fl_O", "duct11.Fl_I") + self.pyc_connect_flow("duct11.Fl_O", "lpt.Fl_I") + self.pyc_connect_flow("lpt.Fl_O", "duct13.Fl_I") + self.pyc_connect_flow("duct13.Fl_O", "core_nozz.Fl_I") + self.pyc_connect_flow("splitter.Fl_O2", "byp_bld.Fl_I") + self.pyc_connect_flow("byp_bld.Fl_O", "duct15.Fl_I") + self.pyc_connect_flow("duct15.Fl_O", "byp_nozz.Fl_I") + + # Bleed flows: + self.pyc_connect_flow("hpc.cool1", "lpt.cool1", connect_stat=False) + self.pyc_connect_flow("hpc.cool2", "lpt.cool2", connect_stat=False) + self.pyc_connect_flow("bld3.cool3", "hpt.cool3", connect_stat=False) + self.pyc_connect_flow("bld3.cool4", "hpt.cool4", connect_stat=False) + + # Specify solver settings: + newton = self.nonlinear_solver = om.NewtonSolver() + newton.options["atol"] = 1e-8 + + # set this very small, so it never activates and we rely on atol + newton.options["rtol"] = 1e-99 + newton.options["iprint"] = 2 + newton.options["maxiter"] = 50 + newton.options["solve_subsystems"] = True + newton.options["max_sub_solves"] = 1000 + newton.options["reraise_child_analysiserror"] = False + # ls = newton.linesearch = BoundsEnforceLS() + ls = newton.linesearch = om.ArmijoGoldsteinLS() + ls.options["maxiter"] = 3 + ls.options["rho"] = 0.75 + # ls.options['print_bound_enforce'] = True + + self.linear_solver = om.DirectSolver() + + super().setup() + + class MPhbtf(pyc.MPCycle): + def setup(self): + + self.pyc_add_pnt( + "DESIGN", HBTF(thermo_method="CEA") + ) # Create an instance of the High Bypass ratio Turbofan + + self.set_input_defaults("DESIGN.inlet.MN", 0.751) + self.set_input_defaults("DESIGN.fan.MN", 0.4578) + self.set_input_defaults("DESIGN.splitter.BPR", 5.105) + self.set_input_defaults("DESIGN.splitter.MN1", 0.3104) + self.set_input_defaults("DESIGN.splitter.MN2", 0.4518) + self.set_input_defaults("DESIGN.duct4.MN", 0.3121) + self.set_input_defaults("DESIGN.lpc.MN", 0.3059) + self.set_input_defaults("DESIGN.duct6.MN", 0.3563) + self.set_input_defaults("DESIGN.hpc.MN", 0.2442) + self.set_input_defaults("DESIGN.bld3.MN", 0.3000) + self.set_input_defaults("DESIGN.burner.MN", 0.1025) + self.set_input_defaults("DESIGN.hpt.MN", 0.3650) + self.set_input_defaults("DESIGN.duct11.MN", 0.3063) + self.set_input_defaults("DESIGN.lpt.MN", 0.4127) + self.set_input_defaults("DESIGN.duct13.MN", 0.4463) + self.set_input_defaults("DESIGN.byp_bld.MN", 0.4489) + self.set_input_defaults("DESIGN.duct15.MN", 0.4589) + self.set_input_defaults("DESIGN.LP_Nmech", 4666.1, units="rpm") + self.set_input_defaults("DESIGN.HP_Nmech", 14705.7, units="rpm") + + # --- Set up bleed values ----- + + self.pyc_add_cycle_param("inlet.ram_recovery", 0.9990) + self.pyc_add_cycle_param("duct4.dPqP", 0.0048) + self.pyc_add_cycle_param("duct6.dPqP", 0.0101) + self.pyc_add_cycle_param("burner.dPqP", 0.0540) + self.pyc_add_cycle_param("duct11.dPqP", 0.0051) + self.pyc_add_cycle_param("duct13.dPqP", 0.0107) + self.pyc_add_cycle_param("duct15.dPqP", 0.0149) + self.pyc_add_cycle_param("core_nozz.Cv", 0.9933) + self.pyc_add_cycle_param("byp_bld.bypBld:frac_W", 0.005) + self.pyc_add_cycle_param("byp_nozz.Cv", 0.9939) + self.pyc_add_cycle_param("hpc.cool1:frac_W", 0.050708) + self.pyc_add_cycle_param("hpc.cool1:frac_P", 0.5) + self.pyc_add_cycle_param("hpc.cool1:frac_work", 0.5) + self.pyc_add_cycle_param("hpc.cool2:frac_W", 0.020274) + self.pyc_add_cycle_param("hpc.cool2:frac_P", 0.55) + self.pyc_add_cycle_param("hpc.cool2:frac_work", 0.5) + self.pyc_add_cycle_param("bld3.cool3:frac_W", 0.067214) + self.pyc_add_cycle_param("bld3.cool4:frac_W", 0.101256) + self.pyc_add_cycle_param("hpc.cust:frac_P", 0.5) + self.pyc_add_cycle_param("hpc.cust:frac_work", 0.5) + self.pyc_add_cycle_param("hpc.cust:frac_W", 0.0445) + self.pyc_add_cycle_param("hpt.cool3:frac_P", 1.0) + self.pyc_add_cycle_param("hpt.cool4:frac_P", 0.0) + self.pyc_add_cycle_param("lpt.cool1:frac_P", 1.0) + self.pyc_add_cycle_param("lpt.cool2:frac_P", 0.0) + self.pyc_add_cycle_param("hp_shaft.HPX", 250.0, units="hp") + + super().setup() + + import time + + prob = om.Problem() + + prob.model = mp_hbtf = MPhbtf() + + prob.setup() + + prob.set_val("DESIGN.fan.PR", 1.685) + prob.set_val("DESIGN.fan.eff", 0.8948) + + prob.set_val("DESIGN.lpc.PR", 1.935) + prob.set_val("DESIGN.lpc.eff", 0.9243) + + prob.set_val("DESIGN.hpc.PR", 9.369) + prob.set_val("DESIGN.hpc.eff", 0.8707) + + prob.set_val("DESIGN.hpt.eff", 0.8888) + prob.set_val("DESIGN.lpt.eff", 0.8996) + + prob.set_val("DESIGN.fc.alt", alt, units="ft") + prob.set_val("DESIGN.fc.MN", MN) + + prob.set_val("DESIGN.T4_MAX", 2857, units="degR") + prob.set_val("DESIGN.Fn_DES", 5900.0, units="lbf") + + # Set initial guesses for balances + prob["DESIGN.balance.FAR"] = 0.025 + prob["DESIGN.balance.W"] = 100.0 + prob["DESIGN.balance.lpt_PR"] = 4.0 + prob["DESIGN.balance.hpt_PR"] = 3.0 + prob["DESIGN.fc.balance.Pt"] = 5.2 + prob["DESIGN.fc.balance.Tt"] = 440.0 + + prob.set_solver_print(level=-1) + # prob.set_solver_print(level=2, depth=1) + + prob.run_model() + + print() + + # Obtaining variables names + + # variable_names = open("variables_hptf.txt", "w") + # print( + # prob.model.list_outputs(val=True, units=True, implicit=False), + # file=variable_names, + # ) + + # BYPASS VARIABLES + T_tot_out_byp = convert_temperature( + (prob.get_val("DESIGN.byp_nozz.throat_total.flow.Fl_O:tot:T")), + "Rankine", + "Celsius", + ) + V_stat_out_byp = prob.get_val("DESIGN.byp_nozz.mux.Fl_O:stat:V") * 0.3048 + MN_out_byp = prob.get_val("DESIGN.byp_nozz.mux.Fl_O:stat:MN") + P_tot_out_byp = ( + prob.get_val("DESIGN.byp_nozz.throat_total.flow.Fl_O:tot:P") * 6894.7573 + ) # Pa + massflow_stat_out_byp = ( + prob.get_val("DESIGN.byp_nozz.mux.Fl_O:stat:W") * 0.45359237 + ) # kg/s + T_stat_out_byp = convert_temperature( + (prob.get_val("DESIGN.byp_nozz.mux.Fl_O:stat:T")), "Rankine", "Celsius" + ) # celsius + + # CORE VARIABLES + T_tot_out_core = convert_temperature( + (prob.get_val("DESIGN.core_nozz.throat_total.flow.Fl_O:tot:T")), + "Rankine", + "Kelvin", + ) + V_stat_out_core = prob.get_val("DESIGN.core_nozz.mux.Fl_O:stat:V") * 0.3048 + MN_out_core = prob.get_val("DESIGN.core_nozz.mux.Fl_O:stat:MN") + P_tot_out_core = ( + prob.get_val("DESIGN.core_nozz.throat_total.flow.Fl_O:tot:P") * 6894.7573 + ) # Pa + massflow_stat_out_core = ( + prob.get_val("DESIGN.core_nozz.mux.Fl_O:stat:W") * 0.45359237 + ) # kg/s + T_stat_out_core = convert_temperature( + (prob.get_val("DESIGN.core_nozz.mux.Fl_O:stat:T")), "Rankine", "Kelvin" + ) # celsius + + res = np.array( + [ + T_tot_out_byp, + V_stat_out_byp, + MN_out_byp, + P_tot_out_byp, + massflow_stat_out_byp, + T_stat_out_byp, + T_tot_out_core, + V_stat_out_core, + MN_out_core, + P_tot_out_core, + massflow_stat_out_core, + T_stat_out_core, + ] + ) + print(f"T_tot_out_core = {T_tot_out_core} [K]") + print(f"V_stat_out_core = {V_stat_out_core} [m/s]") + print(f"MN_out_core = {MN_out_core} [adim]") + print(f"P_tot_out_core = {P_tot_out_core} [Pa]") + print(f"massflow_out_core = {massflow_stat_out_core} [kg/s]") + print(f"T_stat_out_core = {T_stat_out_core} [K]") + print(f"T_tot_out_byp = {T_tot_out_byp} [K]") + print(f"V_stat_out_byp = {V_stat_out_byp} [m/s]") + print(f"MN_out_byp = {MN_out_byp} [adim]") + print(f"P_tot_out_byp = {P_tot_out_byp} [Pa]") + print(f"massflow_stat_out_byp = {massflow_stat_out_byp} [kg/s]") + print(f"T_stat_out_core = {T_stat_out_core} [K]") + + return res + + +if __name__ == "__main__": + + alt = float(input("insert altitude [ft]: ")) + + MN = float(input("insert mach number[adim]: ")) + + run_turbofan_analysis_test(alt, MN) + print("temperature [celsius], stat velocity [m/s], MN for bypass and core") diff --git a/ceasiompy/ThermoData/testing_pycycle/des_turbojet_func.py b/ceasiompy/ThermoData/testing_pycycle/des_turbojet_func.py new file mode 100644 index 000000000..d36447fef --- /dev/null +++ b/ceasiompy/ThermoData/testing_pycycle/des_turbojet_func.py @@ -0,0 +1,173 @@ +import sys + +import openmdao.api as om + +import pycycle.api as pyc + +import numpy as np + +import re + +from scipy.constants import convert_temperature + + +def run_turbojet_analysis(alt, MN, Fn): + class Turbojet(pyc.Cycle): + def setup(self): + + USE_TABULAR = True + + if USE_TABULAR: + self.options["thermo_method"] = "TABULAR" + self.options["thermo_data"] = pyc.AIR_JETA_TAB_SPEC + FUEL_TYPE = "FAR" + + design = self.options["design"] + + # Add engine elements + 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"], + ) + 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)) + + # 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") + + # Make other non-flow connections + # Connect turbomachinery elements to shaft + self.connect("comp.trq", "shaft.trq_0") + self.connect("turb.trq", "shaft.trq_1") + + # Connnect nozzle exhaust to freestream static conditions + self.connect("fc.Fl_O:stat:P", "nozz.Ps_exhaust") + + # Connect outputs to perfomance element + 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") + + # Add balances for design and off-design + balance = self.add_subsystem("balance", om.BalanceComp()) + if design: + + 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") + + balance.add_balance( + "FAR", eq_units="degR", lower=1e-4, val=0.017, rhs_name="T4_target" + ) + self.connect("balance.FAR", "burner.Fl_I:FAR") + self.connect("burner.Fl_O:tot:T", "balance.lhs:FAR") + + balance.add_balance( + "turb_PR", val=1.5, lower=1.001, upper=8, eq_units="hp", rhs_val=0.0 + ) + self.connect("balance.turb_PR", "turb.PR") + self.connect("shaft.pwr_net", "balance.lhs:turb_PR") + + newton = self.nonlinear_solver = om.NewtonSolver() + newton.options["atol"] = 1e-6 + newton.options["rtol"] = 1e-6 + newton.options["iprint"] = 2 + newton.options["maxiter"] = 15 + 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 MPTurbojet(pyc.MPCycle): + def setup(self): + self.pyc_add_pnt("DESIGN", Turbojet()) + + 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) # .2 + self.set_input_defaults("DESIGN.burner.MN", 0.020) # .2 + self.set_input_defaults("DESIGN.turb.MN", 0.4) + + self.pyc_add_cycle_param("burner.dPqP", 0.03) + self.pyc_add_cycle_param("nozz.Cv", 0.99) + + super().setup() + + import time + + # defining the optimization problem + + prob = om.Problem() + mp_turbojet = prob.model = MPTurbojet() + prob.setup(check=False) + + # define input values + prob.set_val("DESIGN.fc.alt", alt, units="ft") + prob.set_val("DESIGN.fc.MN", MN) + prob.set_val("DESIGN.balance.Fn_target", Fn, 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) + + # define initial guesses for the balances + prob["DESIGN.balance.FAR"] = 0.0175506829934 + prob["DESIGN.balance.W"] = 168.453135137 + prob["DESIGN.balance.turb_PR"] = 4.46138725662 + prob["DESIGN.fc.balance.Pt"] = 14.6955113159 + prob["DESIGN.fc.balance.Tt"] = 518.665288153 + + prob.set_solver_print(level=-1) + # prob.set_solver_print(level=2, depth=1) + + prob.run_model() + + print() + + # command to visualize all the output of the system + # a = prob.model.list_outputs(val=True, units=True) + + T_tot_out = convert_temperature( + (prob.get_val("DESIGN.nozz.throat_total.flow.Fl_O:tot:T")), "Rankine", "Kelvin" + ) + V_stat_out = (prob.get_val("DESIGN.nozz.mux.Fl_O:stat:V")) * 0.3048 + MN_out = prob.get_val("DESIGN.nozz.mux.Fl_O:stat:MN") + T_stat_out = convert_temperature( + prob.get_val("DESIGN.nozz.mux.Fl_O:stat:T"), "Rankine", "Kelvin" + ) # celsius + massflow_stat_out = prob.get_val("DESIGN.nozz.mux.Fl_O:stat:W") * 0.45359237 # kg/s + P_tot_out = ( + prob.get_val("DESIGN.nozz.throat_total.flow.Fl_O:tot:P") * 6894.7573 + ) # Pa + + res = np.array( + [T_tot_out, V_stat_out, MN_out, P_tot_out, massflow_stat_out, T_stat_out] + ) + print(f"T_tot_out = {T_tot_out} [K]") + print(f"V_stat_out = {V_stat_out} [m/s]") + print(f"MN_out = {MN_out} [adim]") + print(f"P_tot_out = {P_tot_out} [Pa]") + print(f"massflow_out = {massflow_stat_out} [kg/s]") + print(f"T_stat_out = {T_stat_out} [K]") + return res diff --git a/ceasiompy/ThermoData/thermodata.py b/ceasiompy/ThermoData/thermodata.py new file mode 100644 index 000000000..b7b621fde --- /dev/null +++ b/ceasiompy/ThermoData/thermodata.py @@ -0,0 +1,89 @@ +""" +CEASIOMpy: Conceptual Aircraft Design Software + +Developed by CFS ENGINEERING, 1015 Lausanne, Switzerland + +Calculate outlet conditions fot turbojet and turbofan engines by using the open source code PyCycle +and saving those conditions in a text file + +| Author: Giacomo Benedetti and Francesco Marcucci +| Creation: 2023-12-12 +""" + +# ================================================================================================= +# IMPORTS +# ================================================================================================= + +import shutil +from pathlib import Path + +from ceasiompy.ThermoData.func.run_func import ( + thermo_data_run, +) + +from ceasiompy.utils.ceasiomlogger import get_logger + +from ceasiompy.utils.moduleinterfaces import ( + check_cpacs_input_requirements, + get_toolinput_file_path, + get_tooloutput_file_path, +) +from ceasiompy.utils.commonxpath import ( + REF_XPATH, + CLCALC_XPATH, + ENGINE_TYPE_XPATH, +) +from ceasiompy.utils.commonnames import ( + ENGINE_BOUNDARY_CONDITIONS, +) +from cpacspy.cpacsfunctions import ( + get_value_or_default, +) +from cpacspy.cpacsfunctions import create_branch, get_value, get_value_or_default +from cpacspy.cpacspy import CPACS +from markdownpy.markdownpy import MarkdownDoc +from ceasiompy.utils.ceasiompyutils import get_results_directory + + +log = get_logger() + +MODULE_DIR = Path(__file__).parent +MODULE_NAME = MODULE_DIR.name + + +# ================================================================================================= +# CLASSES +# ================================================================================================= + + +# ================================================================================================= +# FUNCTIONS +# ================================================================================================= + + +# ================================================================================================= +# MAIN +# ================================================================================================= + + +def main(cpacs_path, cpacs_out_path): + + log.info("----- Start of " + MODULE_NAME + " -----") + + results_dir = get_results_directory("ThermoData") + + thermo_data_run(cpacs_path, cpacs_out_path, results_dir) + + folder_name = "reports" + + shutil.rmtree(folder_name) + + log.info("----- End of " + MODULE_NAME + " -----") + + +if __name__ == "__main__": + + cpacs_path = get_toolinput_file_path(MODULE_NAME) + cpacs_out_path = get_tooloutput_file_path(MODULE_NAME) + + main(cpacs_path, cpacs_out_path)