[WIP] Feature GAFs from CFD

doc/jcl_template.py

@@ -14,7 +14,7 @@ class jcl:
 
     def __init__(self):
         # Give your aircraft a name and set some general parameters
-        self.general = {'aircraft': 'DLR F-19-S',
+        self.general = {'aircraft': 'MyAircraft',
                         # Reference span width (from tip to tip)
                         'b_ref': 15.375,
                         # Reference chord length
@@ -30,7 +30,7 @@ def __init__(self):
         to be implemented as a python module.
         """
         # Electronic flight control system
-        self.efcs = {'version': 'mephisto', # Name of the corresponding python module
+        self.efcs = {'version': 'MyEFCS', # Name of the corresponding python module
                      # Path where to find the EFCS module
                      'path': '/path/to/EFCS',
                      }
@@ -70,10 +70,11 @@ def __init__(self):
         self.aero = {'method': 'mona_steady',
                      # 'mona_steady'      - steady trim and quasi-steady time domain simulations
                      # 'mona_unsteady'    - unsteady time domain simulation based on the RFA, e.g. for gust
-                     # 'freq_dom'         - frequency domain simulations, e.g. gust, continuous turbulence, flutter, etc
+                     # 'mona_freq_dom'    - frequency domain simulations, e.g. gust, continuous turbulence, flutter, etc
                      # 'nonlin_steady'    - steady trim and quasi-steady time domain simulations with some non-linearities
                      # 'cfd_steady'       - steady trim
                      # 'cfd_unsteady'     - unsteady time domain simulation, e.g. for gust
+                     # 'cfd_freq_dom'     - frequency domain simulations, e.g. flutter
                      #
                      # True or False, aerodynamic feedback of elastic structure on aerodynamics can be deactivated.
                      # You will still see deformations, but there is no coupling.
@@ -101,12 +102,14 @@ def __init__(self):
                      # number of poles for rational function approximation (RFA)
                      'n_poles': 4,
                      # Additional parameters for CFD
-                     'para_path': '/scratch/tau/',
-                     'para_file': 'para',
+                     'para_path': '/path/to/CFD-workspace',
+                     'para_file': 'my_baseline_para_file',
                      # Currently implemented interfaces: 'tau' or 'su2'
                      'cfd_solver': 'tau',
                      'tau_solver': 'el',
                      'tau_cores': 16,
+                     # Name of job that was used to precompute CFD-based GAFs
+                     'job_name_gafs': 'jcl_dc3_gafs',
                      # --- Start of experimental section, only for special cases ---
                      # Correction coefficient at CG, negativ = destabilizing
                      'Cn_beta_corr': [-0.012],
@@ -126,7 +129,7 @@ def __init__(self):
         # General CFD surface mesh information
         self.meshdefo = {'surface': {'fileformat': 'netcdf', # implemented file formats: 'cgns', 'netcdf', 'su2'
                                      # file name of the CFD mesh
-                                     'filename_grid': 'tau.grid',
+                                     'filename_grid': 'CFD-mesh.grid',
                                      # list of markers [1, 2, ...] or ['upper', 'lower', ...] of surfaces to be included in
                                      # deformation
                                      'markers': [1, 3],
@@ -211,8 +214,8 @@ def __init__(self):
                        'n_blades': [2, 2],
                        # Mach number for VLM4Prop
                        'Ma': [0.25],
-                       # Input-file ('.yaml') for PyPropMAt and VLM4Prop
-                       'propeller_input_file': 'HAP_O6_PROP_pitch.yaml',
+                       # Input-file ('.yaml') for PyPropMat and VLM4Prop
+                       'propeller_input_file': 'Prop.yaml',
                        }
         # CFD-specific
         # In case a pressure inlet boundary is modeled in the CFD mesh, the boundary condition will be
@@ -338,10 +341,17 @@ def __init__(self):
                          # Flutter parameters for k and ke method
                          'flutter_para': {'method': 'k', 'k_red': np.linspace(2.0, 0.001, 1000)},
                          # Flutter parameters for pk method
-                         # There are two implementations of the PK method: 'pk_schwochow', 'pk_rodden'
-                         # Available mode tracking algortihms: 'MAC', 'MAC*PCC' (recommended), 'MAC*HDM'
-                         # 'flutter_para': {'method': 'pk', 'Vtas': np.linspace(100.0, 500.0, 100),
-                         #                  'tracking': 'MAC*PCC'},
+                         # 'flutter_para': {'method': 'pk_schwochow',  # Available implementations: 'pk_schwochow', 'pk_rodden'
+                         #                  'Vtas': np.linspace(100.0, 500.0, 100),
+                         #                  'tracking': 'MAC*PCC',  # Available: 'MAC', 'MAC*PCC' (recommended), 'MAC*HDM'
+                         #                  },
+                         # True or False, enables computation of CFD-based generalized aerodynamic forces (GAFs)
+                         'gaf': True,
+                         'gaf_para': {'df': 1.0,  # The frequency resolution governs the time length.
+                                      # The max. frequency governs the dt of the time domain simulation.
+                                      # Rule of thumb: 20x the highest frequency of interest.
+                                      'fmax': 500.0
+                                      },
                          },
                         ]
         # End

loadskernel/cfd_interfaces/meshdefo.py

@@ -10,26 +10,30 @@ class Meshdefo():
     """
 
     def Ux2(self, Ux2):
-        logging.info('Apply control surface deflections to cfdgrid.')
-        Ujx2 = np.zeros(self.aerogrid['n'] * 6)
-        if 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'y':
-            hingeline = 'y'
-        elif 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'z':
-            hingeline = 'z'
-        else:  # default
-            hingeline = 'y'
-        for i_x2 in range(len(self.x2grid['key'])):
-            logging.debug('Apply deflection of {} for {:0.4f} [deg].'.format(
-                self.x2grid['key'][i_x2], Ux2[i_x2] / np.pi * 180.0))
-            if hingeline == 'y':
-                Ujx2 += np.dot(self.Djx2[i_x2], [0, 0, 0, 0, Ux2[i_x2], 0])
-            elif hingeline == 'z':
-                Ujx2 += np.dot(self.Djx2[i_x2], [0, 0, 0, 0, 0, Ux2[i_x2]])
-        self.transfer_deformations(self.aerogrid, Ujx2, '_k', rbf_type='wendland2', surface_spline=False, support_radius=1.5)
+        if np.any(Ux2):
+            logging.debug('Apply control surface deflections to cfdgrid.')
+            Ujx2 = np.zeros(self.aerogrid['n'] * 6)
+            if 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'y':
+                hingeline = 'y'
+            elif 'hingeline' in self.jcl.aero and self.jcl.aero['hingeline'] == 'z':
+                hingeline = 'z'
+            else:  # default
+                hingeline = 'y'
+            for i_x2 in range(len(self.x2grid['key'])):
+                logging.debug('Apply deflection of {} for {:0.4f} [deg].'.format(
+                    self.x2grid['key'][i_x2], Ux2[i_x2] / np.pi * 180.0))
+                if hingeline == 'y':
+                    Ujx2 += np.dot(self.Djx2[i_x2], [0, 0, 0, 0, Ux2[i_x2], 0])
+                elif hingeline == 'z':
+                    Ujx2 += np.dot(self.Djx2[i_x2], [0, 0, 0, 0, 0, Ux2[i_x2]])
+            self.transfer_deformations(self.aerogrid, Ujx2, '_k', rbf_type='wendland2',
+                                       surface_spline=False, support_radius=1.5)
+        else:
+            logging.debug('Apply NO control surface deflections to cfdgrid.')
 
     def Uf(self, Uf):
-        if 'flex' in self.jcl.aero and self.jcl.aero['flex']:
-            logging.info('Apply flexible deformations to cfdgrid.')
+        if 'flex' in self.jcl.aero and self.jcl.aero['flex'] and np.any(Uf):
+            logging.debug('Apply flexible deformations to cfdgrid.')
             # set-up spline grid
             if self.jcl.spline['splinegrid']:
                 # make sure that there are no double points in the spline grid as this would cause a singularity of the
@@ -45,4 +49,4 @@ def Uf(self, Uf):
 
             self.transfer_deformations(splinegrid, Ug_f_body, '', rbf_type='tps', surface_spline=False)
         else:
-            logging.info('Apply NO flexible deformations to cfdgrid.')
+            logging.debug('Apply NO flexible deformations to cfdgrid.')

loadskernel/cfd_interfaces/su2_interface.py

@@ -113,7 +113,7 @@ def set_deformations(self):
         Communicate the change of coordinates of the fluid interface to the fluid solver.
         Prepare the fluid solver for mesh deformation.
         """
-        logging.info('Sending surface deformations to SU2.')
+        logging.debug('Sending surface deformations to SU2.')
         for x in range(self.local_mesh['n']):
             self.FluidSolver.SetMarkerCustomDisplacement(self.local_mesh['MarkerID'][x],
                                                          self.local_mesh['VertexIndex'][x],
@@ -199,7 +199,7 @@ def init_solver(self):
             self.get_local_mesh()
 
     def run_solver(self, i_timestep=0):
-        logging.info('Waiting until all processes are ready to perform a coordinated start...')
+        logging.debug('Waiting until all processes are ready to perform a coordinated start...')
         self.comm.barrier()
         logging.info('Launch SU2 for time step {}.'.format(i_timestep))
         # start timer
@@ -510,16 +510,23 @@ def update_timedom_para(self):
                 os.symlink(filename_steady, filename_unsteady1)
             except FileExistsError:
                 pass
-
-            """
-            In the time domain, simply rely on the the density residual to establish convergence.
-            This is because SU2 only needs a few inner iterations per time step, which are too few for a meaningful
-            cauchy convergence.
-            """
+            # Set-up inner iterations
+            if 'gaf' in self.simcase and self.simcase['gaf']:
+                # In case of GAF computation, no convergence criterion can be used because the reference / zero solution
+                # (without a pulse) has to have exactly the same number of steps as the pulse solution.
+                # Possibly, the number of inner iterations depends on the aircraft, mesh size, etc.
+                config['INNER_ITER'] = 4
+                if 'CONV_RESIDUAL_MINVAL' in config:
+                    config.pop('CONV_RESIDUAL_MINVAL')
+            else:
+                # In the time domain, simply rely on the the density residual to establish convergence.
+                # This is because SU2 only needs a few inner iterations per time step, which are too few for a meaningful
+                # cauchy convergence.
+                config['INNER_ITER'] = 30
+                config['CONV_RESIDUAL_MINVAL'] = -6
+            # Remove other convergence criteria
             if 'CONV_FIELD' in config:
                 config.pop('CONV_FIELD')
-            config['INNER_ITER'] = 30
-            config['CONV_RESIDUAL_MINVAL'] = -6
 
             # There is no need for restart solutions, they only take up storage space. Write plotting files only.
             config['OUTPUT_FILES'] = ['TECPLOT', 'SURFACE_TECPLOT']
@@ -561,7 +568,6 @@ def update_gust_para(self, Vtas, Vgust):
             # Note: In SU2 this is the full gust length, not the gust gradient H (half gust length).
             config['GUST_WAVELENGTH'] = 2.0 * self.simcase['gust_gradient']
             config['GUST_PERIODS'] = 1.0
-            config['GUST_AMPL'] = Vgust
             config['GUST_BEGIN_TIME'] = 0.0
             config['GUST_BEGIN_LOC'] = -2.0 * self.simcase['gust_gradient'] - self.simcase['gust_para']['T1'] * Vtas
 

loadskernel/equations/cfd_frequency_domain.py

@@ -0,0 +1,129 @@
+import logging
+import numpy as np
+
+from scipy.interpolate import interp1d
+from scipy.fftpack import fft
+
+from loadskernel.interpolate import MatrixInterpolation
+from loadskernel.equations.mona_frequency_domain import KMethod as MonaKMethod
+from loadskernel.equations.mona_frequency_domain import KEMethod as MonaKEMethod
+from loadskernel.equations.mona_frequency_domain import PKMethodRodden as MonaPKMethodRodden
+from loadskernel.equations.mona_frequency_domain import GustExcitation as MonaGustExcitation
+
+
+class GustExcitation(MonaGustExcitation):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        # Initialize additional attributes to avoid defining them outside __init__
+        # Interpolators
+        self.Qhh_interp = None
+        self.Qhk_interp = None
+        self.Qk_gust_interp = None
+
+    def build_AIC_interpolators(self):
+        # Similar as in the fultter solutions, but not scaled with the dynamic pressure q_dyn.
+        Qhh = []
+        for i_k, _ in enumerate(self.GAFs['k_red']):
+            Qhh.append(self.PHIkh.T.dot(self.GAFs['Qhk'][:, :, i_k]))
+        self.Qhh_interp = MatrixInterpolation(self.GAFs['k_red'], Qhh)
+        # Reorder matrices (freq, aero panels, modes) because the interpolator works along the first axis.
+        Qhk = np.moveaxis(self.GAFs['Qhk'], -1, 0)
+        Qk_gust = np.moveaxis(self.GAFs['Qk_gust'], -1, 0)
+        self.Qhk_interp = MatrixInterpolation(self.GAFs['k_red'], Qhk)
+        self.Qk_gust_interp = MatrixInterpolation(self.GAFs['k_red'], Qk_gust)
+
+    def transfer_function(self, f):
+        omega = 2.0 * np.pi * f
+        Qhh = self.Qhh_interp(self.f2k(f))
+        TF = np.linalg.inv(-self.Mhh * omega ** 2 + complex(0, 1) * omega * self.Dhh + self.Khh - Qhh)
+        return TF
+
+    def calc_gust_excitation(self, freqs, t):
+        gust_f = fft(self.one_m_cosine_gust(t))
+        Ph_fourier = np.zeros((self.n_modes, len(freqs)), dtype='complex128')
+        Pk_fourier = np.zeros((self.aerogrid['n'] * 6, len(freqs)), dtype='complex128')
+        for i, f in enumerate(freqs):
+            Qk_gust = self.Qk_gust_interp(self.f2k(f))
+            Pk_fourier[:, i] = Qk_gust.dot(gust_f[i])
+            Ph_fourier[:, i] = self.PHIkh.T.dot(Pk_fourier[:, i])
+        return Ph_fourier, Pk_fourier
+
+    def one_m_cosine_gust(self, t):
+        # This is "only" the gust signal; unlike with panel methods, there is no relationship with the aicraft geometry here.
+        # The effect of the aircraft penetrating into the gust was already captured during the GAF computations.
+        tw = self.simcase['gust_gradient'] * 2.0 / self.Vtas
+        gust = self.Vtas * self.WG_TAS * 0.5 * (1 - np.cos(2.0 * np.pi * t / tw))
+        gust[np.where(t > tw)] = 0.0
+        return gust
+
+    def calc_aero_response(self, freqs, Uh, dUh_dt):
+        # Because the motion is included in the GAF computations, matrix Qhk is multiplied "only" by the generalized
+        # deformations Uh.
+        Ph_fourier = np.zeros((self.n_modes, len(freqs)), dtype='complex128')
+        Pk_fourier = np.zeros((self.aerogrid['n'] * 6, len(freqs)), dtype='complex128')
+        for i, f in enumerate(freqs):
+            Qhk = self.Qhk_interp(self.f2k(f))
+            Pk_fourier[:, i] = Qhk.dot(Uh[:, i])
+            Ph_fourier[:, i] = self.PHIkh.T.dot(Pk_fourier[:, i])
+        return Ph_fourier, Pk_fourier
+
+
+class KMethod(MonaKMethod):
+
+    def __init__(self, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        # initialize frequently used attributes to satisfy linters/static analyzers
+        self.n_freqs = None
+        self.n_modes = None
+        self.k_reds = None
+
+    def build_AIC_interpolators(self):
+        Qhh = []
+        for i_k, _ in enumerate(self.GAFs['k_red']):
+            Qhh.append(self.PHIkh.T.dot(self.GAFs['Qhk'][:, :, i_k]) / self.GAFs['q_dyn'])
+        self.Qhh_interp = interp1d(self.GAFs['k_red'], Qhh, kind='cubic', axis=0, fill_value="extrapolate")
+
+    def setup_frequence_parameters(self):
+        self.n_modes = self.model['mass'][self.trimcase['mass']]['n_modes'][()] + 5
+        self.k_reds = self.simcase['flutter_para']['k_red']
+        self.n_freqs = len(self.k_reds)
+
+        if self.k_reds.max() > np.max(self.GAFs['k_red']):
+            logging.warning('Required reduced frequency = %0.3f but GAFs given only up to %0.3f',
+                            self.k_reds.max(), np.max(self.GAFs['k_red']))
+
+
+class KEMethod(MonaKEMethod, KMethod):
+    """
+    The CFD-based KE-Method uses the combined formulations of the CFD-based K-Method (see above)
+    and the Mona-based KE-Method (imported as MonaKMethod). This is achieved by inheriting twice.
+    """
+
+
+class PKMethodRodden(MonaPKMethodRodden):
+
+    def build_AIC_interpolators(self):
+        # Same formulation as in K-Method, but with custom, linear matrix interpolation
+        Qhh = []
+        for i_k, _ in enumerate(self.GAFs['k_red']):
+            Qhh.append(self.PHIkh.T.dot(self.GAFs['Qhk'][:, :, i_k]) / self.GAFs['q_dyn'])
+        self.Qhh_interp = MatrixInterpolation(self.GAFs['k_red'], Qhh)
+
+    def system(self, k_red):
+        rho = self.atmo['rho']
+        # Make sure that k_red is not zero due to the division by k_red. If k_red=0.0, set to a small value.
+        # This line is the only difference to the mona-based PKMethodRodden, because GAFs from CFD are also
+        # calculated for k_red=0.0.
+        k_red = np.max([k_red, 0.001])
+
+        Qhh = self.Qhh_interp(k_red)
+        Mhh_inv = np.linalg.inv(self.Mhh)
+
+        upper_part = np.concatenate((np.zeros((self.n_modes, self.n_modes)),
+                                     np.eye(self.n_modes)), axis=1)
+        lower_part = np.concatenate((-Mhh_inv.dot(self.Khh - rho / 2 * self.Vtas ** 2.0 * Qhh.real),
+                                     -Mhh_inv.dot(self.Dhh - rho / 4 * self.Vtas * self.macgrid['c_ref'] / k_red * Qhh.imag)),
+                                    axis=1)
+        A = np.concatenate((upper_part, lower_part))
+        return A

loadskernel/equations/cfd.py → loadskernel/equations/cfd_time_domain.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from loadskernel.equations.steady import Steady
+from loadskernel.equations.mona_time_domain import Steady
 from loadskernel.solution_tools import gravitation_on_earth
 
 
@@ -127,7 +127,7 @@ def equations(self, X, t, modus):
         # aerodynamics
         self.cfd_interface.update_general_para()
         self.cfd_interface.update_timedom_para()
-        if self.simcase['gust']:
+        if 'gust' in self.simcase and self.simcase['gust']:
             self.cfd_interface.update_gust_para(Vtas, self.WG_TAS * Vtas)
         self.cfd_interface.init_solver()
         self.cfd_interface.set_euler_transformation(delta_XYZ, PhiThetaPsi)