diff --git a/src/fastoad/models/performances/mission.py b/src/fastoad/models/performances/mission.py new file mode 100644 index 000000000..5972d68e2 --- /dev/null +++ b/src/fastoad/models/performances/mission.py @@ -0,0 +1,252 @@ +""" +Simple module for performances +""" +# This file is part of FAST : A framework for rapid Overall Aircraft Design +# Copyright (C) 2020 ONERA & ISAE-SUPAERO +# FAST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import numpy as np +from scipy.constants import g +from scipy.integrate import cumtrapz +from scipy.interpolate import interp1d +import openmdao.api as om + +from fastoad.constants import FlightPhase +from fastoad.utils.physics import Atmosphere + +CLIMB_MASS_RATIO = 0.97 # = mass at end of climb / mass at start of climb +DESCENT_MASS_RATIO = 0.98 # = mass at end of descent / mass at start of descent +RESERVE_MASS_RATIO = 0.06 # = (weight of fuel reserve)/ZFW +CLIMB_DESCENT_DISTANCE = 500 # in km, distance of climb + descent + + +class _CruiseTimeSpeedDistance(om.ExplicitComponent): + """ + Estimation of time, speed and distance vectors for a given cruise distance + """ + + def initialize(self): + self.options.declare("flight_point_count", 50, types=(int, tuple)) + + def setup(self): + shape = self.options["flight_point_count"] + self.add_input("data:mission:sizing:cruise:time:initial", np.nan, units="s") + self.add_input("data:TLAR:cruise_mach", np.nan, shape=shape) + self.add_input("data:mission:sizing:cruise:altitude", np.nan, shape=shape, units="m") + self.add_input("data:mission:sizing:cruise:distance:initial", np.nan, units="m") + self.add_input("data:mission:sizing:cruise:distance:final", np.nan, units="m") + + self.add_output("data:mission:sizing:cruise:time", shape=shape, units="s") + self.add_output("data:mission:sizing:cruise:time:final", units="s") + self.add_output("data:mission:sizing:cruise:speed", shape=shape, units="m/s") + self.add_output("data:mission:sizing:cruise:distance", shape=shape, units="m") + + self.declare_partials("data:mission:sizing:cruise:time", "*", method="fd") + self.declare_partials("data:mission:sizing:cruise:speed", "*", method="fd") + self.declare_partials("data:mission:sizing:cruise:time:final", "*", method="fd") + self.declare_partials("data:mission:sizing:cruise:distance", "*", method="fd") + + def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None): + flight_point_count = self.options["flight_point_count"] + t_0 = inputs["data:mission:sizing:cruise:time:initial"] + mach = inputs["data:TLAR:cruise_mach"] + initial_distance = inputs["data:mission:sizing:cruise:distance:initial"] + final_distance = inputs["data:mission:sizing:cruise:distance:final"] + + total_distance = final_distance - initial_distance + + dx_distance = total_distance / (flight_point_count - 1) + + atmosphere = Atmosphere( + inputs["data:mission:sizing:cruise:altitude"], altitude_in_feet=False + ) + + speed = atmosphere.speed_of_sound * mach + + distance = np.full((flight_point_count - 1), dx_distance) + distance = np.concatenate((initial_distance, distance)) + distance = np.cumsum(distance) + + time = distance / speed + t_0 + + outputs["data:mission:sizing:cruise:time"] = time + outputs["data:mission:sizing:cruise:time:final"] = time[-1] + outputs["data:mission:sizing:cruise:speed"] = speed + outputs["data:mission:sizing:cruise:distance"] = distance + + +class _CruiseAltitude(om.ExplicitComponent): + """ + Estimation of altitude vector + """ + + def initialize(self): + + altitude_ref = 0.3048 * np.linspace(0.0, 60e3, num=1000) + + rho_ref = Atmosphere(altitude_ref, altitude_in_feet=False).density + + self.options.declare("flight_point_count", 50, types=(int, tuple)) + self.options.declare( + "altitude_interpolation", interp1d(rho_ref, altitude_ref), types=interp1d + ) + + def setup(self): + + shape = self.options["flight_point_count"] + # TODO: is it necessary to keep initial altitude ? + self.add_input("data:mission:sizing:cruise:altitude:initial", np.nan, units="m") + self.add_input("data:mission:sizing:cruise:speed", np.nan, shape=shape, units="m/s") + self.add_input("data:mission:sizing:cruise:weight", np.nan, shape=shape, units="kg") + self.add_input("data:aerodynamics:aircraft:cruise:optimal_CL", np.nan) + self.add_input("data:geometry:aircraft:wing:area", np.nan, units="m**2") + + # TODO: solver needs an initial guess, how to not hard code it ? + self.add_output( + "data:mission:sizing:cruise:altitude", + 0.3048 * np.full(shape, 35000), + shape=shape, + units="m", + ) + + self.declare_partials("data:mission:sizing:cruise:altitude", "*", method="fd") + + def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None): + altitude_interpolation = self.options["altitude_interpolation"] + initial_altitude = inputs["data:mission:sizing:cruise:altitude:initial"] + speed = inputs["data:mission:sizing:cruise:speed"] + weight = inputs["data:mission:sizing:cruise:weight"] + optimum_cl = inputs["data:aerodynamics:aircraft:cruise:optimal_CL"] + wing_area = inputs["data:geometry:aircraft:wing:area"] + + rho = 2 * weight * g / optimum_cl / wing_area / speed ** 2 + + altitude = altitude_interpolation(rho) + altitude[0] = initial_altitude + + outputs["data:mission:sizing:cruise:altitude"] = altitude + + +class _CruiseThrust(om.ExplicitComponent): + """ + Estimation of thrust vector + """ + + def initialize(self): + self.options.declare("flight_point_count", 50, types=(int, tuple)) + + def setup(self): + + shape = self.options["flight_point_count"] + self.add_input("data:mission:sizing:cruise:altitude", np.nan, shape=shape, units="m") + self.add_input("data:TLAR:cruise_mach", np.nan, shape=shape) + self.add_input("data:mission:sizing:cruise:weight", np.nan, shape=shape, units="kg") + self.add_input("data:aerodynamics:aircraft:cruise:optimal_CL", np.nan) + self.add_input("data:aerodynamics:aircraft:cruise:optimal_CD", np.nan) + + self.add_output("data:propulsion:phase", FlightPhase.CRUISE, shape=shape) + self.add_output("data:propulsion:use_thrust_rate", False, shape=shape) + self.add_output( + "data:propulsion:required_thrust_rate", 0.0, shape=shape, lower=0.0, upper=1.0 + ) + self.add_output("data:propulsion:required_thrust", shape=shape, units="N") + self.add_output("data:propulsion:altitude", shape=shape, units="m") + self.add_output("data:propulsion:mach", shape=shape) + + self.declare_partials("data:propulsion:phase", "*", method="fd") + self.declare_partials("data:propulsion:use_thrust_rate", "*", method="fd") + self.declare_partials("data:propulsion:required_thrust_rate", "*", method="fd") + self.declare_partials("data:propulsion:required_thrust", "*", method="fd") + self.declare_partials("data:propulsion:altitude", "*", method="fd") + self.declare_partials("data:propulsion:mach", "*", method="fd") + + def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None): + weight = inputs["data:mission:sizing:cruise:weight"] + optimum_cl = inputs["data:aerodynamics:aircraft:cruise:optimal_CL"] + optimum_cd = inputs["data:aerodynamics:aircraft:cruise:optimal_CD"] + + outputs["data:propulsion:altitude"] = inputs["data:mission:sizing:cruise:altitude"] + outputs["data:propulsion:mach"] = inputs["data:TLAR:cruise_mach"] + + thrust = weight * g / (optimum_cl / optimum_cd) + + outputs["data:propulsion:required_thrust"] = thrust + + +class _CruiseWeight(om.ExplicitComponent): + """ + Estimation of weight vector + """ + + def initialize(self): + self.options.declare("flight_point_count", 50, types=(int, tuple)) + + def setup(self): + + shape = self.options["flight_point_count"] + self.add_input("data:mission:sizing:cruise:consumption:initial", np.nan, units="kg") + self.add_input("data:mission:sizing:cruise:time", np.nan, shape=shape, units="s") + self.add_input("data:mission:sizing:cruise:weight:initial", np.nan, units="kg") + self.add_input("data:propulsion:SFC", np.nan, shape=shape, units="kg/N/s") + self.add_input("data:propulsion:required_thrust", np.nan, shape=shape) + + self.add_output( + "data:mission:sizing:cruise:weight", + np.full(shape, 77037 * CLIMB_MASS_RATIO), + shape=shape, + units="kg", + ) + self.add_output("data:mission:sizing:cruise:consumption", shape=shape, units="kg") + + self.declare_partials("data:mission:sizing:cruise:weight", "*", method="fd") + self.declare_partials("data:mission:sizing:cruise:consumption", "*", method="fd") + + def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None): + initial_consumption = inputs["data:mission:sizing:cruise:consumption:initial"] + time = inputs["data:mission:sizing:cruise:time"] + initial_weight = inputs["data:mission:sizing:cruise:weight:initial"] + sfc = inputs["data:propulsion:SFC"] + thrust = inputs["data:propulsion:required_thrust"] + + consumption = cumtrapz(sfc * thrust, time) + consumption = np.concatenate((initial_consumption, consumption)) + weight = initial_weight - consumption + + outputs["data:mission:sizing:cruise:weight"] = weight + outputs["data:mission:sizing:cruise:consumption"] = consumption + + +class Cruise(om.Group): + """ + Complete Cruise + """ + + def initialize(self): + self.options.declare("flight_point_count", 50, types=(int, tuple)) + + def setup(self): + + shape = self.options["flight_point_count"] + + self.add_subsystem( + "time_speed_distance", + _CruiseTimeSpeedDistance(flight_point_count=shape), + promotes=["*"], + ) + self.add_subsystem("altitude", _CruiseAltitude(flight_point_count=shape), promotes=["*"]) + self.add_subsystem("thrust", _CruiseThrust(flight_point_count=shape), promotes=["*"]) + self.add_subsystem("weight", _CruiseWeight(flight_point_count=shape), promotes=["*"]) + + self.nonlinear_solver = om.NonlinearBlockGS() + self.nonlinear_solver.options["iprint"] = 2 + self.nonlinear_solver.options["maxiter"] = 100 + self.linear_solver = om.LinearBlockGS() diff --git a/src/fastoad/models/performances/tests/test_mission.py b/src/fastoad/models/performances/tests/test_mission.py new file mode 100644 index 000000000..9d42e942d --- /dev/null +++ b/src/fastoad/models/performances/tests/test_mission.py @@ -0,0 +1,178 @@ +# This file is part of FAST : A framework for rapid Overall Aircraft Design +# Copyright (C) 2020 ONERA & ISAE-SUPAERO +# FAST is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import openmdao.api as om +from numpy.testing import assert_allclose +import numpy as np + +from tests.testing_utilities import run_system +from ..breguet import BreguetFromMTOW, BreguetFromOWE +from ..mission import ( + _CruiseTimeSpeedDistance, + _CruiseAltitude, + _CruiseThrust, + Cruise, + CLIMB_MASS_RATIO, +) +from ...propulsion.fuel_engine.rubber_engine import OMRubberEngine + + +def test_cruise_time_speed_distance(): + flight_points_count = 50 + ivc = om.IndepVarComp() + ivc.add_output( + "data:mission:sizing:cruise:altitude", np.full(flight_points_count, 35000), units="ft" + ) + ivc.add_output("data:TLAR:cruise_mach", np.full(flight_points_count, 0.78)) + ivc.add_output("data:mission:sizing:cruise:time:initial", 0.0) + ivc.add_output("data:mission:sizing:cruise:distance:initial", 0.0) + ivc.add_output("data:mission:sizing:cruise:distance:final", 500, units="NM") + + problem = run_system(_CruiseTimeSpeedDistance(flight_point_count=flight_points_count), ivc) + + assert_allclose( + problem["data:mission:sizing:cruise:speed"], + np.full(flight_points_count, 231.298488), + rtol=1e-3, + ) + assert_allclose(problem["data:mission:sizing:cruise:time:final"], 4003, rtol=1) + + +def test_cruise_altitude(): + flight_points_count = 5 + ivc = om.IndepVarComp() + ivc.add_output("data:mission:sizing:cruise:altitude:initial", 35000, units="ft") + ivc.add_output( + "data:mission:sizing:cruise:speed", np.full(flight_points_count, 231.298488), units="m/s" + ) + ivc.add_output( + "data:mission:sizing:cruise:weight", + np.linspace( + 77037 * CLIMB_MASS_RATIO, 77037 * CLIMB_MASS_RATIO * 0.5, num=flight_points_count + ), + units="kg", + ) + ivc.add_output("data:aerodynamics:aircraft:cruise:optimal_CL", 0.52) + ivc.add_output("data:geometry:aircraft:wing:area", 130.233, units="m**2") + + problem = run_system(_CruiseAltitude(flight_point_count=flight_points_count), ivc) + + assert_allclose( + problem["data:mission:sizing:cruise:altitude"], + np.array([10668.0, 11176.0, 12153.0, 13310.0, 14725.0]), + atol=1, + ) + + +def test_cruise_thrust(): + flight_points_count = 5 + ivc = om.IndepVarComp() + ivc.add_output("data:TLAR:cruise_mach", np.full(flight_points_count, 0.78)) + ivc.add_output( + "data:mission:sizing:cruise:altitude", np.full(flight_points_count, 35000), units="ft" + ) + ivc.add_output( + "data:mission:sizing:cruise:weight", + np.linspace( + 77037 * CLIMB_MASS_RATIO, 77037 * CLIMB_MASS_RATIO * 0.5, num=flight_points_count + ), + units="kg", + ) + ivc.add_output("data:aerodynamics:aircraft:cruise:optimal_CL", 0.52) + ivc.add_output("data:aerodynamics:aircraft:cruise:optimal_CD", 0.03414) + ivc.add_output("data:geometry:aircraft:wing:area", 130.233, units="m**2") + + problem = run_system(_CruiseThrust(flight_point_count=flight_points_count), ivc) + + assert_allclose( + problem["data:propulsion:required_thrust"], + np.array([48112.0, 42098.0, 36084.0, 30070.0, 24056.0]), + atol=1, + ) + + +def test_cruise_no_engine(): + flight_points_count = 5 + ivc = om.IndepVarComp() + ivc.add_output("data:mission:sizing:cruise:consumption:initial", 0.0) + ivc.add_output("data:TLAR:cruise_mach", np.full(flight_points_count, 0.78)) + ivc.add_output("data:mission:sizing:cruise:time:initial", 0.0) + ivc.add_output("data:mission:sizing:cruise:distance:initial", 0.0) + ivc.add_output("data:mission:sizing:cruise:distance:final", 500, units="NM") + ivc.add_output("data:mission:sizing:cruise:altitude:initial", 10161.0 / 0.3048, units="ft") + ivc.add_output("data:mission:sizing:cruise:weight:initial", 77037, units="kg") + ivc.add_output("data:aerodynamics:aircraft:cruise:optimal_CL", 0.52) + ivc.add_output("data:aerodynamics:aircraft:cruise:optimal_CD", 0.03414) + ivc.add_output("data:geometry:aircraft:wing:area", 130.233, units="m**2") + ivc.add_output("data:propulsion:SFC", np.full(flight_points_count, 1.706583172872032e-05)) + + group = om.Group() + group.add_subsystem("cruise", Cruise(flight_point_count=flight_points_count), promotes=["*"]) + + problem = run_system(group, ivc) + + assert_allclose( + problem["data:mission:sizing:cruise:altitude"], + np.array([10161.0, 10131.0, 10202.0, 10274.0, 10345.0]), + atol=1, + ) + assert_allclose( + problem["data:mission:sizing:cruise:weight"], + np.array([77037.0, 76201.0, 75372.0, 74551.0, 73737.0]), + atol=1, + ) + + +def test_cruise(): + flight_points_count = 5 + ivc = om.IndepVarComp() + ivc.add_output("data:mission:sizing:cruise:consumption:initial", 0.0) + ivc.add_output("data:TLAR:cruise_mach", np.full(flight_points_count, 0.78)) + ivc.add_output("data:mission:sizing:cruise:time:initial", 0.0) + ivc.add_output("data:mission:sizing:cruise:distance:initial", 0.0) + ivc.add_output("data:mission:sizing:cruise:distance:final", 500, units="NM") + ivc.add_output("data:mission:sizing:cruise:altitude:initial", 10161.0 / 0.3048, units="ft") + ivc.add_output("data:mission:sizing:cruise:weight:initial", 77037, units="kg") + ivc.add_output("data:aerodynamics:aircraft:cruise:optimal_CL", 0.52) + ivc.add_output("data:aerodynamics:aircraft:cruise:optimal_CD", 0.03414) + ivc.add_output("data:geometry:aircraft:wing:area", 130.233, units="m**2") + + ivc.add_output("data:propulsion:rubber_engine:bypass_ratio", 5) + ivc.add_output("data:propulsion:rubber_engine:maximum_mach", 0.95) + ivc.add_output("data:propulsion:rubber_engine:design_altitude", 35000, units="ft") + ivc.add_output("data:propulsion:MTO_thrust", 100000, units="N") + ivc.add_output("data:propulsion:rubber_engine:overall_pressure_ratio", 30) + ivc.add_output("data:propulsion:rubber_engine:turbine_inlet_temperature", 1500, units="K") + + group = om.Group() + group.add_subsystem( + "engine", OMRubberEngine(flight_point_count=flight_points_count), promotes=["*"] + ) + group.add_subsystem("cruise", Cruise(flight_point_count=flight_points_count), promotes=["*"]) + group.nonlinear_solver = om.NonlinearBlockGS() + group.nonlinear_solver.options["iprint"] = 0 + group.nonlinear_solver.options["maxiter"] = 100 + group.linear_solver = om.LinearBlockGS() + + problem = run_system(group, ivc) + + assert_allclose( + problem["data:mission:sizing:cruise:altitude"], + np.array([10161.0, 10152.0, 10242.0, 10330.0, 10415.0]), + atol=1, + ) + assert_allclose( + problem["data:mission:sizing:cruise:weight"], + np.array([77037.0, 75954.0, 74906.0, 73900.0, 72930.0]), + atol=1, + )