Merge pull request #481 from fast-aircraft-design/takeoff_segments_an…

…d_polar

Added only takeoff segments and polar.

src/fastoad/model_base/flight_point.py

@@ -90,6 +90,7 @@ class FlightPoint:
     CL: float = None  # pylint: disable=invalid-name
     #: Drag coefficient.
     CD: float = None  # pylint: disable=invalid-name
+    lift: float = None  #: Aircraft lift in Newtons
     drag: float = None  #: Aircraft drag in Newtons.
     thrust: float = None  #: Thrust in Newtons.
     thrust_rate: float = None  #: Thrust rate (between 0. and 1.)
@@ -101,6 +102,8 @@ class FlightPoint:
     sfc: float = None  #: Specific Fuel Consumption in kg/N/s.
     slope_angle: float = None  #: Slope angle in radians.
     acceleration: float = None  #: Acceleration value in m/s**2.
+    alpha: float = None  #: angle of attack in radians
+    slope_angle_derivative: float = None  #: slope angle derivative in rad/s
     name: str = None  #: Name of current phase.
 
     _units = dict(
@@ -113,12 +116,15 @@ class FlightPoint:
         mach="-",
         CL="-",
         CD="-",
+        lift="N",
         drag="N",
         thrust="N",
         thrust_rate="-",
         sfc="kg/N/s",
         slope_angle="rad",
         acceleration="m/s**2",
+        alpha="rad",
+        slope_angle_derivative="rad/s",
     )
 
     def __post_init__(self):

src/fastoad/models/performances/mission/polar.py

@@ -1,4 +1,5 @@
 """Aerodynamic polar data."""
+
 #  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
 #  Copyright (C) 2021 ONERA & ISAE-SUPAERO
 #  FAST is free software: you can redistribute it and/or modify
@@ -11,28 +12,39 @@
 #  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 <https://www.gnu.org/licenses/>.
-
 from numpy import ndarray
 from scipy.interpolate import interp1d
 from scipy.optimize import fmin
 
 
 class Polar:
-    def __init__(self, cl: ndarray, cd: ndarray):
-
+    def __init__(self, cl: ndarray, cd: ndarray, alpha: ndarray = None):
         """
         Class for managing aerodynamic polar data.
 
         Links drag coefficient (CD) to lift coefficient (CL).
         It is defined by two vectors with CL and CD values.
+        If a vector of angle of attack (alpha) is given, it links alpha and CL
 
         Once defined, for any CL value, CD can be obtained using :meth:`cd`.
+        For any alpha given, CL is obtained using :meth:'cl'.
 
         :param cl: a N-elements array with CL values
         :param cd: a N-elements array with CD values that match CL
+        :param alpha: a N-elements array with angle of attack corresponding to CL values
+
         """
+
         self._definition_CL = cl
-        self._cd = interp1d(cl, cd, kind="quadratic", fill_value="extrapolate")
+        self._definition_CD = cd
+
+        # Interpolate cd
+        self._cd_vs_cl = interp1d(cl, cd, kind="quadratic", fill_value="extrapolate")
+
+        # CL as a function of AoA
+        self._definition_alpha = alpha
+        if alpha is not None:
+            self._cl_vs_alpha = interp1d(alpha, cl, kind="linear", fill_value="extrapolate")
 
         def _negated_lift_drag_ratio(lift_coeff):
             """Returns -CL/CD."""
@@ -45,6 +57,16 @@ def definition_cl(self):
         """The vector that has been used for defining lift coefficient."""
         return self._definition_CL
 
+    @property
+    def definition_cd(self):
+        """The vector that has been used for defining drag coefficient."""
+        return self._definition_CD
+
+    @property
+    def definition_alpha(self):
+        """The vector that has been used for defining AoA."""
+        return self._definition_alpha
+
     @property
     def optimal_cl(self):
         """The CL value that provides larger lift/drag ratio."""
@@ -59,5 +81,17 @@ def cd(self, cl=None):
         :return: CD values for each provide CL values
         """
         if cl is None:
-            return self._cd(self._definition_CL)
-        return self._cd(cl)
+            return self._cd_vs_cl(self._definition_CL)
+        return self._cd_vs_cl(cl)
+
+    def cl(self, alpha):
+        """
+        The lift coefficient corresponding to alpha (rad)
+
+        :param alpha: the angle of attack at which CL is evaluated
+        :return: CL value for each alpha.
+        """
+        if self._definition_alpha is None:
+            raise ValueError("Polar was instantiated without AoA vector.")
+
+        return self._cl_vs_alpha(alpha)

src/fastoad/models/performances/mission/polar_modifier.py

@@ -0,0 +1,112 @@
+""" Aerodynamics polar modifier."""
+#  This file is part of FAST-OAD : A framework for rapid Overall Aircraft Design
+#  Copyright (C) 2021 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 <https://www.gnu.org/licenses/>.
+
+from abc import ABC, abstractmethod
+from dataclasses import dataclass
+
+from fastoad.model_base.flight_point import FlightPoint
+from .polar import Polar
+
+
+@dataclass
+class AbstractPolarModifier(ABC):
+
+    """
+    Base class to implement a change to the polar during the mission computation
+    """
+
+    @abstractmethod
+    def modify_polar(self, polar: Polar, flight_point: FlightPoint) -> Polar:
+
+        """
+        :param polar: an instance of Polar
+        :param flight_point: an intance of FlightPoint containg only floats
+        :return: the modified polar for the flight point
+        """
+
+
+@dataclass
+class UnchangedPolar(AbstractPolarModifier):
+    """
+    Default polar modifier returning the polar without changes
+    """
+
+    def modify_polar(self, polar: Polar, flight_point: FlightPoint) -> Polar:
+        """
+        :param polar: a polar instance
+        :param flight_point: a FlightPoint instance
+        :return: the polar instance
+        """
+        return polar
+
+
+@dataclass
+class GroundEffectRaymer(AbstractPolarModifier):
+
+    """
+    Evaluates the drag in ground effect, using Raymer's model:
+        'Aircraft Design A conceptual approach', D. Raymer p304
+    """
+
+    #: Wingspan
+    span: float
+
+    #: Main landing gear height
+    landing_gear_height: float
+
+    #: Induced drag coefficient, multiplies CL**2 to obtain the induced drag
+    induced_drag_coef: float
+
+    #: Winglet effect tuning coefficient
+    k_winglet: float
+
+    #: Total drag tuning coefficient
+    k_cd: float
+
+    #: Altitude of ground w.r.t. sea level
+    ground_altitude: float = 0.0
+
+    def modify_polar(self, polar: Polar, flight_point: FlightPoint) -> Polar:
+        """
+        Compute the ground effect based on altitude from ground and return an updated polar
+
+        :param polar: a Polar instance used as basis to apply ground effect
+        :param flight_point: a flight point containing the flight conditions
+        for calculation of ground effect
+        :return: a copy of polar with ground effect
+        """
+
+        h_b = (
+            self.span * 0.1
+            + self.landing_gear_height
+            + flight_point.altitude
+            - self.ground_altitude
+        ) / self.span
+        k_ground = 33.0 * h_b ** 1.5 / (1 + 33.0 * h_b ** 1.5)
+        cd_ground = (
+            self.induced_drag_coef
+            * polar.definition_cl ** 2
+            * self.k_winglet
+            * self.k_cd
+            * (k_ground - 1)
+        )
+
+        # Update polar interpolation
+        modified_polar = Polar(
+            polar.definition_cl,
+            polar.definition_cd + cd_ground,
+            polar.definition_alpha,
+        )
+
+        return modified_polar

src/fastoad/models/performances/mission/segments/base.py

@@ -34,6 +34,7 @@
 from .exceptions import FastUnknownMissionSegmentError
 from ..base import IFlightPart
 from ..exceptions import FastFlightSegmentIncompleteFlightPoint
+from ..polar_modifier import AbstractPolarModifier, UnchangedPolar
 
 _LOGGER = logging.getLogger(__name__)  # Logger for this module
 
@@ -202,6 +203,7 @@ def compute_from(self, start: FlightPoint) -> pd.DataFrame:
         # Let's ensure we do not modify the original definitions of start and target
         # during the process
         start_copy = deepcopy(start)
+
         if start_copy.altitude is not None:
             try:
                 self.complete_flight_point(start_copy)
@@ -317,6 +319,10 @@ class AbstractTimeStepFlightSegment(
     #: The Polar instance that will provide drag data.
     polar: Polar = MANDATORY_FIELD
 
+    #: A polar modifier that can apply dynamic changes to the original polar
+    # (the default value returns a polar without change)
+    polar_modifier: AbstractPolarModifier = UnchangedPolar()
+
     #: The reference area, in m**2.
     reference_area: float = MANDATORY_FIELD
 
@@ -330,7 +336,7 @@ class AbstractTimeStepFlightSegment(
 
     #: Minimum and maximum authorized mach values. If computed Mach gets beyond these limits,
     #: computation will be interrupted and a warning message will be issued in logger.
-    mach_bounds: tuple = (0.0, 5.0)
+    mach_bounds: tuple = (-1.0e-6, 5.0)
 
     #: If True, computation will be interrupted if a parameter stops getting closer to target
     #: between two iterations (which can mean the provided thrust rate is not adapted).
@@ -390,8 +396,9 @@ def complete_flight_point(self, flight_point: FlightPoint):
             )
 
             if self.polar:
+                modified_polar = self.polar_modifier.modify_polar(self.polar, flight_point)
                 flight_point.CL = flight_point.mass * g / reference_force
-                flight_point.CD = self.polar.cd(flight_point.CL)
+                flight_point.CD = modified_polar.cd(flight_point.CL)
             else:
                 flight_point.CL = flight_point.CD = 0.0
             flight_point.drag = flight_point.CD * reference_force
@@ -623,3 +630,80 @@ def get_distance_to_target(
     ) -> float:
         current = flight_points[-1]
         return target.time - current.time
+
+
+@dataclass
+class TakeOffSegment(AbstractManualThrustSegment, ABC):
+    """
+    Class for computing takeoff segment
+    """
+
+    # Default time step for this dynamic segment
+    time_step: float = 0.1
+
+    def compute_from_start_to_target(self, start: FlightPoint, target: FlightPoint) -> pd.DataFrame:
+        self.polar_modifier.ground_altitude = start.altitude
+        return super().compute_from_start_to_target(start, target)
+
+
+@dataclass
+class GroundSegment(TakeOffSegment, ABC):
+    """
+    Class for computing accelerated segments on the ground with wheel friction.
+    """
+
+    # Friction coefficient considered for acceleration at take-off.
+    # The default value is representative of dry concrete/asphalte
+    wheels_friction: float = 0.03
+
+    # The angle of attack of the aircraft at the beginning of the segment
+    alpha: float = 0.0
+
+    def get_gamma_and_acceleration(self, flight_point: FlightPoint):
+        """
+        For ground segment, gamma is assumed always 0 and wheel friction
+        (with or without brake) is added to drag
+        """
+        mass = flight_point.mass
+        drag_aero = flight_point.drag
+        lift = flight_point.lift
+        thrust = flight_point.thrust
+
+        drag = drag_aero + (mass * g - lift) * self.wheels_friction
+
+        # edit flight_point fields
+        flight_point.drag = drag
+
+        acceleration = (thrust - drag) / mass
+
+        return 0.0, acceleration
+
+    def complete_flight_point(self, flight_point: FlightPoint):
+        """
+        Computes data for provided flight point using AoA and apply polar modification if any
+
+        :param flight_point: the flight point that will be completed in-place
+        """
+        self._complete_speed_values(flight_point)
+
+        # Ground segment may force engine setting like reverse or idle
+        flight_point.engine_setting = self.engine_setting
+
+        atm = self._get_atmosphere_point(flight_point.altitude)
+        reference_force = 0.5 * atm.density * flight_point.true_airspeed ** 2 * self.reference_area
+
+        if self.polar:
+            alpha = flight_point.alpha
+            modified_polar = self.polar_modifier.modify_polar(self.polar, flight_point)
+            flight_point.CL = modified_polar.cl(alpha)
+            flight_point.CD = modified_polar.cd(flight_point.CL)
+        else:
+            flight_point.CL = flight_point.CD = 0.0
+
+        flight_point.drag = flight_point.CD * reference_force
+        flight_point.lift = flight_point.CL * reference_force
+
+        self.compute_propulsion(flight_point)
+        flight_point.slope_angle, flight_point.acceleration = self.get_gamma_and_acceleration(
+            flight_point
+        )