# -*- coding: utf-8 -*-
# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
#
# 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
# Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
# Contributors:
# INITIAL AUTHORS - initial API and implementation
# and/or initial documentation
# :author: Sobieski, Agte, and Sandusky
# OTHER AUTHORS - MACROSCOPIC CHANGES
# :author: Damien Guenot
# :author: Francois Gallard
# From NASA/TM-1998-208715
# Bi-Level Integrated System Synthesis (BLISS)
# Sobieski, Agte, and Sandusky
"""
SSBJ Aerodynamics computation
*****************************
"""
from __future__ import division, unicode_literals
import logging
from math import pi
from numpy import array, atleast_2d, cos, sin, sqrt, zeros
LOGGER = logging.getLogger(__name__)
DEG_TO_RAD = pi / 180.0
[docs]class SobieskiAerodynamics(object):
"""Class defining aerodynamical analysis for Sobieski problem and related method to
the aerodynamical problem such as disciplines computation, constraints, reference
optimum."""
DTYPE_COMPLEX = "complex128"
DTYPE_DOUBLE = "float64"
PRESSURE_GRADIENT_LIMIT = 1.04
def __init__(self, sobieski_base):
"""Constructor."""
self.base = sobieski_base
self.constants = self.base.default_constants()
(
self.x_initial,
self.tc_initial,
self.half_span_initial,
self.aero_center_initial,
self.cf_initial,
self.mach_initial,
self.h_initial,
self.throttle_initial,
self.lift_initial,
self.twist_initial,
self.esf_initial,
) = self.base.get_initial_values()
self.dtype = self.base.dtype
self.math = self.base.math
def __set_coeff_drag_f1(self, y_32, x_2):
"""Setting of flags used for determination of polynomial coefficients.
:param y_32: shared variables coming from blackbox_propulsion
- y_32[0]: engine scale factor
:type y_32: numpy array
:param x_2: aero. design variable
- x_2[0]: friction coeff
:type x_2: numpy array
:returns: s_initial, S_new, flag, bound
- s_initial: reference value
- s_new: new values
:rtype: numpy array
"""
s_initial = array([self.esf_initial, self.cf_initial], dtype=self.dtype)
s_new = array([y_32[0], x_2[0]], dtype=self.dtype)
flag = array([1, 1], dtype=self.dtype)
bound = array([0.25, 0.25], dtype=self.dtype)
return s_initial, s_new, flag, bound
def __set_coeff_drag_f2(self, y_12):
"""Setting of flags used for determination of polynomial coefficients.
:param y_12: shared variables coming from blackbox_structure
:type y_12: numpy array
:returns: s_initial, s_new, flag, bound
- s_initial: reference value
- s_new: new values
:rtype: numpy array
"""
s_initial = array([self.twist_initial], dtype=self.dtype)
s_new = array([y_12[1]], dtype=self.dtype)
flag = array([5], dtype=self.dtype)
bound = array([0.25], dtype=self.dtype)
return s_initial, s_new, flag, bound
def __set_coeff_drag_f3(self, x_shared):
"""Setting of flags used for determination of polynomial coefficients.
:param x_shared: global design variables
:type x_shared: numpy array
:returns: s_initial, s_new, flag, bound
- s_initial: reference value
- s_new: new values
:rtype: numpy array
"""
s_initial = array([self.tc_initial], dtype=self.dtype)
s_new = array([x_shared[0]], dtype=self.dtype)
flag = array([1], dtype=self.dtype)
bound = array([0.25], dtype=self.dtype)
return s_initial, s_new, flag, bound
[docs] def compute_k_aero(self, x_shared):
"""Computation of a induced drag coefficient (related to lift)
:param x_shared: global design variables
:type x_shared: numpy array
:returns: k_aero
:rtype: numpy array
"""
# if x_shared[2] >= 1: # Mach
# k = (x_shared[2] ** 2 - 1) * \
# self.math.cos(x_shared[4] * DEG_TO_RAD) / \
# (4.0 * self.math.sqrt(x_shared[4] ** 2 - 1) - 2)
# else:
# k = 1. / (self.math.pi * 0.8 * x_shared[3])
k_aero = (
(x_shared[2] ** 2 - 1)
* self.math.cos(x_shared[4] * DEG_TO_RAD)
/ (4.0 * self.math.sqrt(x_shared[4] ** 2 - 1) - 2)
)
return k_aero
# @staticmethod
# def compute_dk_aero_dAR(self, x_shared):
# # if x_shared[2] >= 1:
# # dk_dAR = 0.0
# # else:
# # dk_dAR = -1.0 / (0.8 * self.math.pi *
# abs(x_shared[3])**2)
# dk_dAR = 0.0
# return dk_dAR
[docs] @staticmethod
def compute_dk_aero_dsweep(x_shared):
"""Computation of a derivative of k_aero wrt sweep.
:param x_shared: global design variables
:type x_shared: numpy array
:returns: dk_aero_dsweep
:rtype: numpy array
"""
# if x_shared[2] >= 1:
# dk_dsweep = -DEG_TO_RAD * (x_shared[2] ** 2 - 1) * \
# self.math.sin(x_shared[4] * DEG_TO_RAD) / \
# (4.0 * self.math.sqrt(x_shared[4] ** 2 - 1) - 2)
# else:
# dk_dsweep = 0.0
u_velo = (x_shared[2] * x_shared[2] - 1.0) * cos(x_shared[4] * DEG_TO_RAD)
up_velo = (
-DEG_TO_RAD
* (x_shared[2] * x_shared[2] - 1.0)
* sin(x_shared[4] * DEG_TO_RAD)
)
v_velo = 4.0 * sqrt(x_shared[4] * x_shared[4] - 1.0) - 2.0
vp_velo = 4.0 * x_shared[4] * (x_shared[4] * x_shared[4] - 1.0) ** -0.5
dk_aero_dsweep = (up_velo * v_velo - u_velo * vp_velo) / v_velo ** 2
return dk_aero_dsweep
[docs] @staticmethod
def compute_dk_aero_dmach(x_shared):
"""Computation of a derivative of k_aero wrt Mach.
:param x_shared: global design variables
:type x_shared: numpy array
:returns: dk_aero_dmach
:rtype: numpy array
"""
# if x_shared[2] >= 1:
# dk_dM = abs(x_shared[3]) * (2.0 * abs(x_shared[2])) * \
# np.cos(x_shared[4] * np.pi / 180.) /
# (4. * abs(x_shared[3]) * \
# np.sqrt(abs(x_shared[4]**2 - 1.) - 2.))
# else:
# dk_dM = 0.0
dk_aero_dmach = (
(2.0 * x_shared[2])
* cos(x_shared[4] * pi / 180.0)
/ (4.0 * sqrt(x_shared[4] ** 2 - 1.0) - 2.0)
)
return dk_aero_dmach
[docs] def compute_dadimcf_dcf(self, x_2):
"""Computation of a derivative of adim friction coeff of polynomial wrt friction
coeff.
:param x_2: local design variables
:type x_2: numpy array
:returns: derivative of adim friction coeff of polynomial
wrt friction coeff
:rtype: numpy array
"""
dadimcf_dcf = self.base.derive_normalize_s(self.cf_initial, x_2[0])
return dadimcf_dcf
[docs] def compute_dadimtwist_dtwist(self, y_12):
"""Computation of a derivative of adim twist of polynomial wrt twist.
:param y_12: coupling variable from weight analysis
:type y_12: numpy array
:returns: derivative of adim twist of polynomial wrt twist
:rtype: numpy array
"""
dadimtwist_dtwist = self.base.derive_normalize_s(self.twist_initial, y_12[1])
return dadimtwist_dtwist
[docs] def compute_dadimesf_desf(self, y_32):
"""Computation of a derivative of adim ESF of polynomial wrt ESF.
:param y_32: coupling variable from propulsion analysis
:type y_32: numpy array
:returns: derivative of adim ESF of polynomial wrt ESF
:rtype: numpy array
"""
return self.base.derive_normalize_s(self.esf_initial, y_32[0])
[docs] def compute_dadimtaper_dtaper(self, x_shared):
"""Computation of a derivative of adim taper-ratio of polynomial wrt taper-
ratio.
:param x_shared: global design variables
:type x_shared: numpy array
:returns: derivative of adim taper-ratio of polynomial wrt taper-ratio
:rtype: numpy array
"""
return self.base.derive_normalize_s(self.tc_initial, x_shared[0])
[docs] def compute_cd_min(self, x_shared, fo1):
"""Computation of a 2D minimum drag coefficient.
:param x_shared: global design variables
:type x_shared: numpy array
:param fo1: coefficient for engine size
:type fo1: numpy array
:returns: CDmin : drag coefficient
:rtype: numpy array
"""
cdmin_incomp = self.constants[4]
return (
cdmin_incomp * fo1
+ 3.05
* x_shared[0] ** (5.0 / 3.0)
* (self.math.cos(x_shared[4] * DEG_TO_RAD)) ** 1.5
)
[docs] @staticmethod
def compute_dcdmin_dsweep(x_shared):
"""Computation of derivative of 2D minimum drag coefficient wrt sweep.
:param x_shared: global design variables
:type x_shared: numpy array
:returns: dCDmin_dsweep
:rtype: numpy array
"""
ang_rad = x_shared[4] * DEG_TO_RAD
return (
-3.05
* 1.5
* x_shared[0] ** (5.0 / 3.0)
* cos(ang_rad) ** 0.5
* DEG_TO_RAD
* sin(ang_rad)
)
[docs] def compute_cd(self, x_shared, lift_coeff, fo1, fo2):
"""Computation of total drag coefficient.
:param x_shared: global design variables
:type x_shared: numpy array
:param lift_coeff: lift coefficient
:type lift_coeff: numpy array
:param fo1: coefficient for engine size
:type fo1: numpy array
:param fo2: coefficient for twist influence on drag
:type fo2: numpy array
:returns: drag_coefficient : drag coefficient
:rtype: numpy array
"""
k_aero = self.compute_k_aero(x_shared)
cdmin = self.compute_cd_min(x_shared, fo1)
return fo2 * (cdmin + k_aero * lift_coeff * lift_coeff)
[docs] def compute_drag(self, drag_coeff, x_shared):
"""Computation of drag from drag coefficient.
:param drag_coeff: lift coefficient
:type drag_coeff: numpy array
:param x_shared: global design variables
:type x_shared: numpy array
:returns: drag as a force
:rtype: numpy array
"""
return self.compute_force(drag_coeff, x_shared[5], x_shared[2], x_shared[1])
[docs] def compute_dcd_dsweep(self, x_shared, lift_coeff, fo2):
"""Computation of derivative of drag coefficient wrt sweep.
:param x_shared: global design variables
:type x_shared: numpy array
:param lift_coeff: lift coefficient
:type lift_coeff: numpy array
:param fo2: coefficient for twist influence on drag
:type fo2: numpy array
:returns: dCd/dsweep
:rtype: numpy array
"""
dcdmin_dsweep = self.compute_dcdmin_dsweep(x_shared)
dk_dsweep = self.compute_dk_aero_dsweep(x_shared)
return fo2 * (dcdmin_dsweep + lift_coeff * lift_coeff * dk_dsweep)
[docs] def compute_dcd_dsref(self, x_shared, y_12, fo2):
"""Computation of derivative of drag coefficient wrt sweep.
:param x_shared: global design variables
:type x_shared: numpy array
:param y_12: lift coefficient
:type y_12: numpy array
:param fo2: coefficient for twist influence on drag
:type fo2: numpy array
:returns: dCd/dsweep
:rtype: numpy array
"""
k_aero = self.compute_k_aero(x_shared)
lift_coeff = self.compute_cl(x_shared, y_12)
dcl_dsref = self.compute_dcl_dsref(x_shared, y_12)
return 2 * k_aero * lift_coeff * dcl_dsref * fo2
[docs] def compute_dcd_dmach(self, x_shared, y_12, fo2):
"""Computation of derivative of drag coefficient wrt Mach.
:param x_shared: global design variables
:type x_shared: numpy array
:param y_12: coupling variable from weight analysis
:type y_12: numpy array
:param fo2: coefficient for twist influence on drag
:type fo2: numpy array
:returns: dcd_dmach
:rtype: numpy array
"""
k_aero = self.compute_k_aero(x_shared)
dk_dmach = self.compute_dk_aero_dmach(x_shared)
lift_coeff = self.compute_cl(x_shared, y_12)
dcl_dmach = self.compute_dcl_dmach(x_shared, y_12)
return (2.0 * k_aero * dcl_dmach + lift_coeff * dk_dmach) * lift_coeff * fo2
[docs] def compute_cl(self, x_shared, y_12):
"""Computation of lift coefficient.
:param x_shared: global design variables
:type x_shared: numpy array
:param y_12: coupling variable from weight analysis
:type y_12: numpy array
:returns: lift coefficient Cl
:rtype: numpy array
"""
return self.compute_adim_coeff(y_12[0], x_shared[5], x_shared[2], x_shared[1])
[docs] def compute_dcl_dh(self, x_shared, y_12):
"""Computation of derivative of lift coefficient wrt altitude.
:param x_shared: global design variables
:type x_shared: numpy array
:param y_12: coupling variable from weight analysis
:type y_12: numpy array
:returns: derivative of lift coefficient wrt altitude
:rtype: numpy array
"""
rhov2 = self.compute_rhov2(x_shared)
drhov2_dh = self.compute_drhov2_dh(x_shared)
return -2 * y_12[0] / x_shared[5] * drhov2_dh / (rhov2 ** 2)
[docs] def compute_dcl_dmach(self, x_shared, y_12):
"""Computation of derivative of lift coefficient wrt reference surface.
:param x_shared: global design variables
:type x_shared: numpy array
:param y_12: coupling variable from weight analysis
:type y_12: numpy array
:returns: derivative of lift coefficient wrt reference surface
:rtype: numpy array
"""
rho, velocity = self.compute_rho_v(x_shared[2], x_shared[1])
dv_dmach = self.compute_dv_dmach(x_shared[1])
dcl_dmach = (
-4.0
* y_12[0]
* dv_dmach
/ (rho * x_shared[5] * velocity * velocity * velocity)
)
return dcl_dmach
[docs] def compute_dcl_dsref(self, x_shared, y_12):
"""Computation of derivative of lift coefficient wrt Mach number.
:param x_shared: global design variables
:type x_shared: numpy array
:param y_12: coupling variable from weight analysis
:type y_12: numpy array
:returns: derivative of lift coefficient wrt Mach number
:rtype: numpy array
"""
return -self.compute_cl(x_shared, y_12) / x_shared[5]
[docs] def compute_rhov2(self, x_shared):
"""Computation of rho * velocity * velocity (2*dynamic pressure)
:param x_shared: global design variables
:type x_shared: numpy array
:returns: rho * velocity * velocity
:rtype: numpy array
"""
rho, velocity = self.compute_rho_v(x_shared[2], x_shared[1])
return rho * velocity * velocity
[docs] def compute_drhov2_dh(self, x_shared):
"""Computation of derivative of rhoV**2 wrt altitude.
:param x_shared: global design variables
:type x_shared: numpy array
:returns: drhov2_dh (derivative of rhoV**2 wrt altitude)
:rtype: numpy array
"""
rho, velocity = self.compute_rho_v(x_shared[2], x_shared[1])
drho_dh, dv_dh = self.compute_drho_dh_dv_dh(x_shared[2], x_shared[1])
return drho_dh * velocity * velocity + 2.0 * rho * dv_dh * velocity
[docs] def compute_drhov2_dmach(self, x_shared):
"""Computation of derivative of rhoV**2 wrt Mach.
:param x_shared: global design variables
:type x_shared: numpy array
:returns: drhov2_dmach (derivative of rhoV**2 wrt Mach)
:rtype: numpy array
"""
rho, velocity = self.compute_rho_v(x_shared[2], x_shared[1])
dv_dmach = self.compute_dv_dmach(x_shared[1])
return 2.0 * rho * dv_dmach * velocity
[docs] def compute_rho_v(self, mach, altitude):
"""From given Mach number and altitude, compute velocity and density.
:param mach: Mach number
:type mach: float
:param altitude: altitude
:type altitude: float
:returns: rho, velocity : density and velocity
:rtype: float
"""
if altitude.real < 36089.0:
velocity = mach * 1116.39 * self.math.sqrt(1 - 6.875e-6 * altitude)
rho = 2.377e-3 * (1 - 6.875e-6 * altitude) ** 4.2561
else:
velocity = mach * 968.1
rho = 2.377e-3 * 0.2971 * self.math.exp((36089.0 - altitude) / 20806.7)
return rho, velocity
[docs] def compute_dv_dmach(self, altitude):
"""Computation of derivative of velocity wrt Mach.
:param altitude: altitude
:type altitude: numpy array
:returns: drhov2_dmach (derivative of rhoV**2 wrt Mach)
:rtype: numpy array
"""
if altitude.real < 36089.0:
dv_dmach = 1116.39 * self.math.sqrt(1 - 6.875e-6 * altitude)
else:
dv_dmach = 968.1
return dv_dmach
[docs] def compute_drho_dh_dv_dh(self, mach, altitude):
"""Compute derivative of velocity and density wrt altitude.
:param mach: Mach number
:type mach: float
:param altitude: altitude
:type altitude: float
:returns: drho_dh, dv_dh
:rtype: float
"""
if altitude.real <= 36089.0:
drho_dh = (
-2.377e-3 * 4.2561 * 6.875e-6 * (1.0 - 6.875e-6 * altitude) ** 3.2561
)
dv_dh = (
-6.875e-6 * 1116.39 * mach * 0.5 * (1.0 - 6.875e-6 * altitude) ** -0.5
)
else:
drho_dh = (
-2.377e-3
* 0.2971
/ 20806.7
* self.math.exp((36089.0 - altitude) / 20806.7)
)
dv_dh = 0.0
return drho_dh, dv_dh
[docs] def compute_adim_coeff(self, force, sref, mach, altitude):
"""Compute adim. force coeff from force (lift or drag)
:param force: Force
:type force: float
:param sref: reference surface
:type sref: float
:param mach: Mach number
:type mach: float
:param altitude: altitude
:type altitude: float
:returns: force coefficient: Cl or Cd according to input force
:rtype: float
"""
rho, velocity = self.compute_rho_v(mach, altitude)
return force / (0.5 * rho * velocity * velocity * sref)
[docs] def compute_force(self, adim_coeff, sref, mach, altitude):
"""Compute force (lift or drag) from adim coeff.
:param adim_coeff: force coefficient
:type adim_coeff: float
:param sref: reference surface
:type sref: float
:param mach: Mach number
:type mach: float
:param altitude: altitude
:type altitude: float
:returns: force: lift or drag
:rtype: float
"""
rho, velocity = self.compute_rho_v(mach, altitude)
# return adim_coeff * 0.5 * rho * V * V * sref
return 0.5 * rho * velocity * velocity * adim_coeff * sref
[docs] def blackbox_aerodynamics(self, x_shared, y_12, y_32, x_2, true_cstr=False):
"""This function calculates drag and lift to drag ratio of A/C.
:param x_shared: shared design variable vector:
- x_shared[0]: thickness/chord ratio
- x_shared[1]: altitude
- x_shared[2]: Mach
- x_shared[3]: aspect ratio
- x_shared[4]: wing sweep
- x_shared[5]: wing surface area
:type x_shared: numpy array
:param y_12: shared variables coming from blackbox_structure:
- y_12[0]: total aircraft weight
- y_12[1]: wing twist
:type y_12: numpy array
:param y_32: shared variables coming from blackbox_propulsion:
- y_32[0]: engine scale factor
:type y_32: numpy array
:param x_2: aero. design variable:
- x_2[0]: friction coeff
:type x_2: numpy array
:param true_cstr: Default value = False)
:returns: y_2, y_21, y_23, y_24, g_2
- y_2: aero. analysis outputs
- y_2[0]: lift
- y_2[1]: drag
- y_2[2]: lift/drag ratio
- y_21: shared variable for blackbox_structure (lift)
- y_23: shared variable for blackbox_propulsion (drag)
- y_24: shared variable for BlackBoxMission (lift/drag ratio)
- g_2: aero constraint (pressure gradient)
:rtype: numpy array, numpy array, numpy array, numpy array, numpy array
"""
y_2 = zeros(3, dtype=self.dtype)
y_23 = zeros(1, dtype=self.dtype)
y_24 = zeros(1, dtype=self.dtype)
y_21 = zeros(1, dtype=self.dtype)
g_2 = zeros(1, dtype=self.dtype)
lift_coeff = self.compute_cl(x_shared, y_12)
# Modification of CDmin for ESF and Cf
s_initial1, s_new1, flag1, bound1 = self.__set_coeff_drag_f1(y_32, x_2)
fo1 = self.base.poly_approx(s_initial1, s_new1, flag1, bound1)
# Modification of drag_coeff for wing twist
s_initial2, s_new2, flag2, bound2 = self.__set_coeff_drag_f2(y_12)
fo2 = self.base.poly_approx(s_initial2, s_new2, flag2, bound2)
drag_coeff = self.compute_cd(x_shared, lift_coeff, fo1, fo2)
drag = self.compute_drag(drag_coeff, x_shared)
y_2[1] = drag
y_2[2] = lift_coeff / drag_coeff
y_2[0] = y_12[0]
y_23[0] = y_2[1]
y_24[0] = y_2[2]
y_21[0] = y_2[0]
# Computation of total drag of A/C
s_initial3, s_new3, flag3, bound3 = self.__set_coeff_drag_f3(x_shared)
# adverse pressure gradient
g_2[0] = self.base.poly_approx(s_initial3, s_new3, flag3, bound3)
# Custom Pgrad function: replace by a linear function
# coeff_dir=(1.04-0.96)/(0.06-0.04)# = 4
# g_2[0] = coeff_dir*Z[0]+0.8
if not true_cstr:
g_2[0] = g_2[0] - self.PRESSURE_GRADIENT_LIMIT
return y_2, y_21, y_23, y_24, g_2
@staticmethod
def __derive_liftoverdrag(cl_cd, lift_jacobian, drag_jacobian, inv_drag):
"""Compute lift over drag jacobian terms.
:param cl_cd: lift over drag
:type cl_cd: numpy array
:param lift_jacobian: d(lift)/d(var)
:type lift_jacobian: numpy array
:param drag_jacobian: d(drag)/d(var)
:type drag_jacobian: numpy array
:returns: jacobian of lift over drag
:rtype: numpy array
"""
return inv_drag * (lift_jacobian - cl_cd * drag_jacobian)
@staticmethod
def __set_coupling_jacobian(jacobian):
"""Set jacobian of coupling variables."""
jacobian["y_21"]["x_2"] = atleast_2d(jacobian["y_2"]["x_2"][0, :])
jacobian["y_21"]["x_shared"] = atleast_2d(jacobian["y_2"]["x_shared"][0, :])
jacobian["y_21"]["y_12"] = atleast_2d(jacobian["y_2"]["y_12"][0, :])
jacobian["y_21"]["y_32"] = atleast_2d(jacobian["y_2"]["y_32"][0, :])
jacobian["y_23"]["x_2"] = atleast_2d(jacobian["y_2"]["x_2"][1, :])
jacobian["y_23"]["x_shared"] = atleast_2d(jacobian["y_2"]["x_shared"][1, :])
jacobian["y_23"]["y_12"] = atleast_2d(jacobian["y_2"]["y_12"][1, :])
jacobian["y_23"]["y_32"] = atleast_2d(jacobian["y_2"]["y_32"][1, :])
jacobian["y_24"]["x_2"] = atleast_2d(jacobian["y_2"]["x_2"][2, :])
jacobian["y_24"]["x_shared"] = atleast_2d(jacobian["y_2"]["x_shared"][2, :])
jacobian["y_24"]["y_12"] = atleast_2d(jacobian["y_2"]["y_12"][2, :])
jacobian["y_24"]["y_32"] = atleast_2d(jacobian["y_2"]["y_32"][2, :])
return jacobian
def __initialize_jacobian(self):
"""Initialization of jacobian matrix.
:returns: jacobian
:rtype: dict(dict(ndarray))
"""
jacobian = {"y_2": {}, "g_2": {}, "y_21": {}, "y_23": {}, "y_24": {}}
jacobian["y_2"]["x_shared"] = zeros((3, 6), dtype=self.dtype)
jacobian["y_2"]["x_2"] = zeros((3, 1), dtype=self.dtype)
jacobian["y_2"]["y_12"] = zeros((3, 2), dtype=self.dtype)
jacobian["y_2"]["y_32"] = zeros((3, 1), dtype=self.dtype)
jacobian["g_2"]["x_2"] = zeros((1, 1), dtype=self.dtype)
jacobian["g_2"]["x_shared"] = zeros((1, 6), dtype=self.dtype)
jacobian["g_2"]["y_12"] = zeros((1, 2), dtype=self.dtype)
jacobian["g_2"]["y_32"] = zeros((1, 1), dtype=self.dtype)
return jacobian
[docs] def derive_blackbox_aerodynamics(self, x_shared, y_12, y_32, x_2):
"""This function calculates drag and lift to drag ratio of A/C.
:param x_shared: shared design variable vector:
- x_shared[0]: thickness/chord ratio
- x_shared[1]: altitude
- x_shared[2]: Mach
- x_shared[3]: aspect ratio
- x_shared[4]: wing sweep
- x_shared[5]: wing surface area
:type x_shared: numpy array
:param y_12: shared variables coming from blackbox_structure:
- y_12[0]: total aircraft weight
- y_12[1]: wing twist
:type y_12: numpy array
:param y_32: shared variables coming from blackbox_propulsion:
- y_32[0]: engine scale factor
:type y_32: numpy array
:param x_2: aero. design variable:
- x_2[0]: friction coeff
:type x_2: numpy array
:returns: jacobian matrix of partial derivatives
:rtype: dict(dict(ndarray))
"""
jacobian = self.__initialize_jacobian()
lift_coeff = self.compute_cl(x_shared, y_12)
# Modification of drag_coeff_min for ESF and Cf
s_initial1, s_new1, flag1, bound1 = self.__set_coeff_drag_f1(y_32, x_2)
fo1, ai_coeff1, aij_coeff1, s_shifted1 = self.base.derive_poly_approx(
s_initial1, s_new1, flag1, bound1
)
drag_coeff_incomp = self.constants[4]
drag_coeff_min = self.compute_cd_min(x_shared, fo1)
k_aero = self.compute_k_aero(x_shared)
# Modification of drag_coeff for wing twist
s_initial2, s_new2, flag2, bound2 = self.__set_coeff_drag_f2(y_12)
fo2, ai_coeff2, aij_coeff2, s_shifted2 = self.base.derive_poly_approx(
s_initial2, s_new2, flag2, bound2
)
drag_coeff = self.compute_cd(x_shared, lift_coeff, fo1, fo2)
cl_cd = lift_coeff / drag_coeff
cl2 = lift_coeff * lift_coeff
drag = self.compute_drag(drag_coeff, x_shared)
rhov2 = self.compute_rhov2(x_shared)
dyn_pressure = 0.5 * rhov2
dyn_force = dyn_pressure * x_shared[5]
inv_drag = 1.0 / drag
# dLift_dCf
# jacobian['y_2']['x_2'][0, 0] = 0.0
# dDrag_dCf
dadimcf_dcf = self.compute_dadimcf_dcf(x_2)
dy2dx2 = dyn_force * fo2 * drag_coeff_incomp * dadimcf_dcf
dy2dx2 *= (
ai_coeff1[1]
+ aij_coeff1[1, 0] * s_shifted1[0]
+ aij_coeff1[1, 1] * s_shifted1[1]
) # dDrag_dCf
jacobian["y_2"]["x_2"][1, 0] = dy2dx2
# d(Lift/Drag)_dCf
dlod = self.__derive_liftoverdrag(
cl_cd, jacobian["y_2"]["x_2"][0, 0], jacobian["y_2"]["x_2"][1, 0], inv_drag
)
jacobian["y_2"]["x_2"][2, 0] = dlod
# dLift/dZ= 0.0
# dDrag/d(t/c)
dy2dxs = dyn_force * fo2 * 3.05 * 5.0 / 3.0 * x_shared[0] ** (2.0 / 3.0)
dy2dxs *= (self.math.cos(x_shared[4] * DEG_TO_RAD)) ** 1.5
jacobian["y_2"]["x_shared"][1, 0] = dy2dxs
# d(Lift/Drag)/d(t/c)
# d(Drag)/dh
drhov2_dh = self.compute_drhov2_dh(x_shared)
dcl_dh = self.compute_dcl_dh(x_shared, y_12)
dy2dxs = 0.5 * x_shared[5] * fo2 * drag_coeff_min * drhov2_dh
add_der = cl2 * drhov2_dh + rhov2 * 2.0 * lift_coeff * dcl_dh
dy2dxs += 0.5 * x_shared[5] * k_aero * fo2 * add_der
jacobian["y_2"]["x_shared"][1, 1] = dy2dxs
# d(Lift/Drag)/dh
# d(Drag)/dM
dcd_dmach = self.compute_dcd_dmach(x_shared, y_12, fo2)
drhov2_dmach = self.compute_drhov2_dmach(x_shared)
jacobian["y_2"]["x_shared"][1, 2] = (
0.5 * x_shared[5] * (drag_coeff * drhov2_dmach + rhov2 * dcd_dmach)
)
# d(Lift/Drag)/dM
# d(Drag)/dAR
# dk_dAR = self.compute_dk_aero_dAR(x_shared) # = 0.0
# jacobian['y_2']['x_shared'][1, 3] = dyn_force * dk_dAR * fo1 * cl2
# d(Lift/Drag)/dAR
# d(Drag)/dsweep
dcd_dsweep = self.compute_dcd_dsweep(x_shared, lift_coeff, fo2)
jacobian["y_2"]["x_shared"][1, 4] = dyn_force * dcd_dsweep
# d(Drag)/dsref
dcd_dsef = self.compute_dcd_dsref(x_shared, y_12, fo2)
dy2dxs = drag / x_shared[5] + dyn_force * dcd_dsef
jacobian["y_2"]["x_shared"][1, 5] = dy2dxs
for i in range(6):
dy2dxs0i = jacobian["y_2"]["x_shared"][0, i]
dy2dxs1i = jacobian["y_2"]["x_shared"][1, i]
dlod = self.__derive_liftoverdrag(cl_cd, dy2dxs0i, dy2dxs1i, inv_drag)
jacobian["y_2"]["x_shared"][2, i] = dlod
# dLift/dWt
jacobian["y_2"]["y_12"][0, 0] = 1.0
# d(Drag)/dWt
jacobian["y_2"]["y_12"][1, 0] = 2.0 * y_12[0] * k_aero * fo2 / dyn_force
# d(Lift/Drag)/dWt
# jacobian['y_2']['y_12'][
# 2, 0] = self.__derive_liftoverdrag(
# cl_cd, jacobian['y_2']['y_12'][0, 0],
# jacobian['y_2']['x_shared'][1, 0], inv_drag)
# dLift/dtwist= 0.0
# d(Drag)/dtwist
dadimtwist_dadim = self.compute_dadimtwist_dtwist(y_12)
jacobian["y_2"]["y_12"][1, 1] = (
dyn_force
* dadimtwist_dadim
* (drag_coeff_min + k_aero * lift_coeff * lift_coeff)
* (ai_coeff2[0] + aij_coeff2[0, 0] * s_shifted2[0])
)
# d(Lift/Drag)/dtwist
# jacobian['y_2']['y_12'][
# 2, 1] = inv_drag * (jacobian['y_2']['y_12'][0, 1] -
# cl_cd * jacobian['y_2']['y_12'][1, 1])
for i in range(2):
dy2dy120 = jacobian["y_2"]["y_12"][0, i]
dy2dy121 = jacobian["y_2"]["y_12"][1, i]
dy2dy122 = self.__derive_liftoverdrag(cl_cd, dy2dy120, dy2dy121, inv_drag)
jacobian["y_2"]["y_12"][2, i] = dy2dy122
# dLift/dESF
# jacobian['y_2']['y_32'][0, 0] = 0.0
# dDrag/dESF
dadimesf_desf = self.compute_dadimesf_desf(y_32)
jacobian["y_2"]["y_32"][1, 0] = (
dyn_force
* fo2
* drag_coeff_incomp
* dadimesf_desf
* (
ai_coeff1[0]
+ aij_coeff1[0, 0] * s_shifted1[0]
+ aij_coeff1[0, 1] * s_shifted1[1]
)
)
# d(Lift/Drag)/dtwist
jacobian["y_2"]["y_32"][2, 0] = self.__derive_liftoverdrag(
cl_cd,
jacobian["y_2"]["y_32"][0, 0],
jacobian["y_2"]["y_32"][1, 0],
inv_drag,
)
# d(dp/dx)/d(t/c)
s_initial3, s_new3, flag3, bound3 = self.__set_coeff_drag_f3(x_shared)
_, ai_coeff3, aij_coeff3, s_shifted3 = self.base.derive_poly_approx(
s_initial3, s_new3, flag3, bound3
)
dadimtaper_dtaper = self.compute_dadimtaper_dtaper(x_shared)
jacobian["g_2"]["x_shared"][0, 0] = dadimtaper_dtaper * (
ai_coeff3[0] + aij_coeff3[0, 0] * s_shifted3[0]
)
# jacobian['g_2']['x_shared'][0, 1:] = 0.0
# jacobian['g_2']['x_2'][0, :] = 0.0
# jacobian['g_2']['y_12'][0, :] = 0.0
# jacobian['g_2']['y_32'][0, :] = 0.0
jacobian = self.__set_coupling_jacobian(jacobian)
return jacobian
