Added new flutter solution as implemented by Rodden in Nastran

doc/jcl_template.py

@@ -335,10 +335,13 @@ def __init__(self):
                          'support': [0, 1, 2, 3, 4, 5],
                          # True or False, enables flutter check with k, ke or pk method
                          'flutter': False,
-                         # flutter parameters for k and ke method
+                         # 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
-                         # 'flutter_para': {'method': 'pk', 'Vtas': np.linspace(100.0, 500.0, 100)},
+                         # Flutter parameters for pk method
+                         # There are two implementations of the PK method: 'pk_schwochow', 'pk_rodden'
+                         # Available mode tracking algortihms: 'MAC', 'MACXP', 'MAC*PCC' (recommended)
+                         # 'flutter_para': {'method': 'pk', 'Vtas': np.linspace(100.0, 500.0, 100),
+                         #                  'tracking': 'MAC*PCC'},
                          },
                         ]
         # End

loadskernel/equations/frequency_domain.py

@@ -520,7 +520,17 @@ def calc_eigenvalues(self):
         self.Vtas = np.array(Vtas)
 
 
-class PKMethod(KMethod):
+class PKMethodSchwochow(KMethod):
+    """
+    This PK-Method uses a formulation proposed by Schwochow [1].
+    Summary: The aerodynamic forces are split in a velocity and a deformation dependent part and added to the damping and
+    stiffness term in the futter equation respectively. In this way, the aerodynamic damping and stiffness are treated
+    seperately and in a more physical way. According to Schwochow, this leads to a better approximation of the damping in the
+    flutter solution.
+
+    [1] Schwochow, J., “Die aeroelastische Stabilitätsanalyse - Ein praxisnaher Ansatz Intervalltheoretischen Betrachtung von
+    Modellierungsunsicherheiten am Flugzeug zur”, Dissertation, Universität Kassel, Kassel, 2012.
+    """
 
     def setup_frequence_parameters(self):
         self.n_modes_rbm = 5
@@ -540,14 +550,16 @@ def eval_equations(self):
         self.setup_frequence_parameters()
 
         logging.info('building systems')
-        self.build_AIC_interpolators()  # unsteady
+        self.build_AIC_interpolators()
         logging.info('starting p-k iterations to match k_red with Vtas and omega')
-        # compute initial guess at k_red=0.0 and first flight speed
+        # Compute initial guess at k_red=0.0 and first flight speed
         self.Vtas = self.Vvec[0]
-        eigenvalue, eigenvector = linalg.eig(self.system(k_red=0.0).real)
+        eigenvalue, eigenvector = linalg.eig(self.system(k_red=0.0))
         bandbreite = eigenvalue.__abs__().max() - eigenvalue.__abs__().min()
-        idx_pos = np.where(eigenvalue.__abs__() / bandbreite >= 1e-3)[0]  # no zero eigenvalues
-        idx_sort = np.argsort(np.abs(eigenvalue.imag[idx_pos]))  # sort result by eigenvalue
+        # No zero eigenvalues
+        idx_pos = np.where(eigenvalue.__abs__() / bandbreite >= 1e-3)[0]
+        # Sort initial results by eigenvalue
+        idx_sort = np.argsort(np.abs(eigenvalue.imag[idx_pos]))
         eigenvalues0 = eigenvalue[idx_pos][idx_sort]
         eigenvectors0 = eigenvector[:, idx_pos][:, idx_sort]
         k0 = eigenvalues0.imag * self.macgrid['c_ref'] / 2.0 / self.Vtas
@@ -557,38 +569,47 @@ def eval_equations(self):
         freqs = []
         damping = []
         Vtas = []
-        # loop over modes
+        # Loop over modes
         for i_mode in range(len(eigenvalues0)):
             logging.debug('Mode {}'.format(i_mode + 1))
             eigenvalues_per_mode = []
             eigenvectors_per_mode = []
             k_old = copy.deepcopy(k0[i_mode])
+            eigenvalues_old = copy.deepcopy(eigenvalues0)
             eigenvectors_old = copy.deepcopy(eigenvectors0)
-            # loop over Vtas
+            # Loop over Vtas
             for _, V in enumerate(self.Vvec):
                 self.Vtas = V
                 e = 1.0
                 n_iter = 0
-                # iteration to match k_red with Vtas and omega of the mode under investigation
+                # Iteration to match k_red with Vtas and omega of the mode under investigation
                 while e >= 1e-3:
-                    eigenvalues_new, eigenvectors_new = self.calc_eigenvalues(self.system(k_old).real, eigenvectors_old)
-                    k_now = np.abs(eigenvalues_new[i_mode].imag) * self.macgrid['c_ref'] / 2.0 / self.Vtas
+                    eigenvalues_new, eigenvectors_new = self.calc_eigenvalues(self.system(k_old),
+                                                                              eigenvalues_old, eigenvectors_old)
+                    # Small switch since this function is reused by class PKMethodRodden
+                    if self.simcase['flutter_para']['method'] in ['pk', 'pk_schwochow']:
+                        # For this implementaion, the reduced frequency may become negative.
+                        k_now = eigenvalues_new[i_mode].imag * self.macgrid['c_ref'] / 2.0 / self.Vtas
+                    elif self.simcase['flutter_para']['method'] in ['pk_rodden']:
+                        # Allow only positive reduced frequencies in the implementation following Rodden.
+                        k_now = np.abs(eigenvalues_new[i_mode].imag) * self.macgrid['c_ref'] / 2.0 / self.Vtas
                     # Use relaxation for improved convergence, which helps in some cases to avoid oscillations of the
                     # iterative solution.
                     k_new = k_old + 0.8 * (k_now - k_old)
                     e = np.abs(k_new - k_old)
                     k_old = k_new
                     n_iter += 1
-                    if n_iter > 80:
-                        logging.warning('poor convergence for mode {} at Vtas={:.2f} with k_red={:.5f} and e={:.5f}'.format(
+                    # If no convergence is achieved, stop and issue a warning. Typically, the iteration converges in less than
+                    # ten loops, so 50 should be more than enough and prevents excessive calculation times.
+                    if n_iter > 50:
+                        logging.warning('No convergence for mode {} at Vtas={:.2f} with k_red={:.5f} and e={:.5f}'.format(
                             i_mode + 1, V, k_new, e))
-                    if n_iter > 100:
-                        logging.warning('p-k iteration NOT converged after 100 loops.')
                         break
+                eigenvalues_old = eigenvalues_new
                 eigenvectors_old = eigenvectors_new
                 eigenvalues_per_mode.append(eigenvalues_new[i_mode])
                 eigenvectors_per_mode.append(eigenvectors_new[:, i_mode])
-            # store
+            # Store results
             eigenvalues_per_mode = np.array(eigenvalues_per_mode)
             eigenvalues.append(eigenvalues_per_mode)
             eigenvectors.append(np.array(eigenvectors_per_mode).T)
@@ -605,22 +626,29 @@ def eval_equations(self):
                     }
         return response
 
-    def calc_eigenvalues(self, A, eigenvector_old):
+    def calc_eigenvalues(self, A, eigenvalues_old, eigenvectors_old):
         eigenvalue, eigenvector = linalg.eig(A)
-        # idx_pos = np.where(eigenvalue.imag >= 0.0)[0]  # nur oszillierende Eigenbewegungen
-        # idx_sort = np.argsort(np.abs(eigenvalue.imag[idx_pos]))  # sort result by eigenvalue
-        MAC = fem_helper.calc_MAC(eigenvector_old, eigenvector, plot=False)
+        # To match the modes with the previous step, use a correlation cirterion as specified in the JCL.
+        # The MAC matrix is used as default.
+        if 'tracking' not in self.simcase['flutter_para']:
+            MAC = fem_helper.calc_MAC(eigenvectors_old, eigenvector)
+        elif self.simcase['flutter_para']['tracking'] == 'MAC':
+            MAC = fem_helper.calc_MAC(eigenvectors_old, eigenvector)
+        elif self.simcase['flutter_para']['tracking'] == 'MACXP':
+            MAC = fem_helper.calc_MACXP(eigenvectors_old, eigenvector)
+        elif self.simcase['flutter_para']['tracking'] == 'MAC*PCC':
+            # This is a combination of modal assurance criterion and pole correlation and improves the handling of complex
+            # conjugate poles.
+            MAC = fem_helper.calc_MAC(eigenvectors_old, eigenvector) * fem_helper.calc_PCC(eigenvalues_old, eigenvalue)
         idx_pos = self.get_best_match(MAC)
-        eigenvalues = eigenvalue[idx_pos]  # [idx_sort]
-        eigenvectors = eigenvector[:, idx_pos]  # [:, idx_sort]
+        eigenvalues = eigenvalue[idx_pos]
+        eigenvectors = eigenvector[:, idx_pos]
         return eigenvalues, eigenvectors
 
     def get_best_match(self, MAC):
         """
-        Before: idx_pos = [MAC[x, :].argmax() for x in range(MAC.shape[0])]
-        With purely real eigenvalues it happens that the selection (only) by highest MAC value does not work.
-        The result is that from two different eigenvalues one is take twice. The solution is to keep record
-        of the matches that are still available so that, if the bets match is already taken, the second best match is selected.
+        It is important that no pole is dropped or selected twice. The solution is to keep record of the matches that are
+        still available so that, if the best match is already taken, the second best match is selected.
         """
         possible_matches = [True] * MAC.shape[1]
         possible_idx = np.arange(MAC.shape[1])
@@ -641,26 +669,17 @@ def calc_Qhh_2(self, Qjj_unsteady):
         return self.PHIlh.T.dot(self.aerogrid['Nmat'].T.dot(self.aerogrid['Amat'].dot(Qjj_unsteady).dot(self.Djh_2)))
 
     def build_AIC_interpolators(self):
-        # do some pre-multiplications first, then the interpolation
         Qhh_1 = []
         Qhh_2 = []
-        if self.jcl.aero['method'] in ['freq_dom']:
-            for Qjj_unsteady in self.aero['Qjj_unsteady']:
-                Qhh_1.append(self.calc_Qhh_1(Qjj_unsteady))
-                Qhh_2.append(self.calc_Qhh_2(Qjj_unsteady))
-        elif self.jcl.aero['method'] in ['mona_unsteady']:
-            ABCD = self.aero['ABCD']
-            for k_red in self.aero['k_red']:
-                D = np.zeros((self.aerogrid['n'], self.aerogrid['n']), dtype='complex')
-                j = 1j  # imaginary number
-                for i_pole, beta in zip(np.arange(0, self.aero['n_poles']), self.aero['betas']):
-                    D += ABCD[3 + i_pole, :, :] * j * k_red / (j * k_red + beta)
-                Qjj_unsteady = ABCD[0, :, :] + ABCD[1, :, :] * j * k_red + ABCD[2, :, :] * (j * k_red) ** 2 + D
-                Qhh_1.append(self.calc_Qhh_1(Qjj_unsteady))
-                Qhh_2.append(self.calc_Qhh_2(Qjj_unsteady))
-
-        self.Qhh_1_interp = MatrixInterpolation(self.aero['k_red'], Qhh_1)
-        self.Qhh_2_interp = MatrixInterpolation(self.aero['k_red'], Qhh_2)
+        # Mirror the AIC matrices with respect to the real axis to allow negative reduced frequencies.
+        Qjj_mirrored = np.concatenate((np.flip(self.aero['Qjj_unsteady'].conj(), axis=0), self.aero['Qjj_unsteady']), axis=0)
+        k_mirrored = np.concatenate((np.flip(-self.aero['k_red']), self.aero['k_red']))
+        # Do some pre-multiplications first, then the interpolation
+        for Qjj in Qjj_mirrored:
+            Qhh_1.append(self.calc_Qhh_1(Qjj))
+            Qhh_2.append(self.calc_Qhh_2(Qjj))
+        self.Qhh_1_interp = MatrixInterpolation(k_mirrored, Qhh_1)
+        self.Qhh_2_interp = MatrixInterpolation(k_mirrored, Qhh_2)
 
     def system(self, k_red):
         rho = self.atmo['rho']
@@ -674,5 +693,40 @@ def system(self, k_red):
         lower_part = np.concatenate((-Mhh_inv.dot(self.Khh - rho / 2.0 * self.Vtas ** 2.0 * Qhh_1),
                                      -Mhh_inv.dot(self.Dhh - rho / 2.0 * self.Vtas * Qhh_2)), axis=1)
         A = np.concatenate((upper_part, lower_part))
+        return A
 
+
+class PKMethodRodden(PKMethodSchwochow):
+    """
+    This PK-Method uses a formulation as implemented in Nastran by Rodden and Johnson [2], Section 2.6, Equation (2-131).
+    Summary: The matrix of the aerodynamic forces Qhh includes both a velocity and a deformation dependent part. The real and
+    the imaginary parts are then added to the damping and stiffness term in the futter equation respectively.
+
+    [2] Rodden, W., and Johnson, E., MSC.Nastran Version 68 Aeroelastic Analysis User’s Guide. MSC.Software Corporation, 2010.
+    """
+
+    def build_AIC_interpolators(self):
+        # Same formulation as in K-Method, but with custom, linear matrix interpolation
+        Qhh = []
+        for Qjj_unsteady, k_red in zip(self.aero['Qjj_unsteady'], self.aero['k_red']):
+            Qhh.append(self.PHIlh.T.dot(self.aerogrid['Nmat'].T.dot(self.aerogrid['Amat'].dot(Qjj_unsteady).dot(
+                self.Djh_1 + complex(0, 1) * k_red / (self.macgrid['c_ref'] / 2.0) * self.Djh_2))))
+        self.Qhh_interp = MatrixInterpolation(self.aero['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. In addition, limit k_red to the smallest
+        # frequency the AIC matrices were calculated for.
+        k_red = np.max([k_red, np.min(self.aero['k_red'])])
+
+        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), dtype='complex128'),
+                                     np.eye(self.n_modes, dtype='complex128')), 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/state_space.py

@@ -4,10 +4,10 @@
 
 from scipy import linalg
 
-from loadskernel.equations.frequency_domain import PKMethod
+from loadskernel.equations.frequency_domain import PKMethodSchwochow
 
 
-class StateSpaceAnalysis(PKMethod):
+class StateSpaceAnalysis(PKMethodSchwochow):
 
     def eval_equations(self):
         self.setup_frequence_parameters()
@@ -19,19 +19,20 @@ def eval_equations(self):
         bandbreite = eigenvalue.__abs__().max() - eigenvalue.__abs__().min()
         idx_pos = np.where(eigenvalue.__abs__() / bandbreite >= 1e-3)[0]  # no zero eigenvalues
         idx_sort = np.argsort(np.abs(eigenvalue.imag[idx_pos]))  # sort result by eigenvalue
-        # eigenvalues0 = eigenvalue[idx_pos][idx_sort]
+        eigenvalues0 = eigenvalue[idx_pos][idx_sort]
         eigenvectors0 = eigenvector[:, idx_pos][:, idx_sort]
 
         eigenvalues = []
         eigenvectors = []
         freqs = []
         damping = []
         Vtas = []
+        eigenvalues_old = copy.deepcopy(eigenvalues0)
         eigenvectors_old = copy.deepcopy(eigenvectors0)
         # loop over Vtas
         for _, V in enumerate(self.Vvec):
             Vtas_i = V
-            eigenvalues_i, eigenvectors_i = self.calc_eigenvalues(self.system(Vtas_i), eigenvectors_old)
+            eigenvalues_i, eigenvectors_i = self.calc_eigenvalues(self.system(Vtas_i), eigenvalues_old, eigenvectors_old)
 
             # store
             eigenvalues.append(eigenvalues_i)
@@ -40,6 +41,7 @@ def eval_equations(self):
             damping.append(eigenvalues_i.real / np.abs(eigenvalues_i))
             Vtas.append([Vtas_i] * len(eigenvalues_i))
 
+            eigenvalues_old = eigenvalues_i
             eigenvectors_old = eigenvectors_i
 
         response = {'eigenvalues': np.array(eigenvalues),
@@ -85,7 +87,7 @@ def system(self, Vtas):
         return A
 
 
-class JacobiAnalysis(PKMethod):
+class JacobiAnalysis(PKMethodSchwochow):
 
     def __init__(self, response):
         self.response = response

loadskernel/fem_interfaces/fem_helper.py

@@ -19,17 +19,80 @@ def calc_MAC(X, Y, plot=False):
     # This loop is necessary in case the size of the the eigenvectors differs.
     for jj in range(Y.shape[1]):
         MAC[:, jj] = np.real(np.abs(q3[:, jj]) ** 2.0 / q1 / q2[jj])
-    # Optionally, vizualize the results
+    # Optionally, visualize the results
     if plot:
         plt.figure()
         plt.pcolor(MAC, cmap='hot_r')
         plt.colorbar()
         plt.grid('on')
+        plt.title('MAC')
 
         return MAC, plt
     return MAC
 
 
+def calc_MACXP(lam_1, my_1, lam_2, my_2, plot=False):
+    """
+    This is function vectorizes the calculation of the pole-weighted modal assurance criterion (MACXP), see equation 28 in [1].
+    Notation: eigenvalues = lam_1,2 and eigenvectors = my_1,2
+
+    [1] Vacher, P., Jacquier, B., and Bucharles, A., “Extensions of the MAC criterion to complex mode”, in Proceedings of the
+    international conference on noise and vibration engineering., Leuven, Belgium, 2010.
+    """
+    # For the MACXP criterion, the number of modes has to be the same
+    assert len(lam_1) == len(lam_2), 'Number of eigenvalues not equal: {} vs. {}'.format(len(lam_1), len(lam_2))
+    assert my_1.shape[0] == my_1.shape[1], 'Size of eigenvector 1 not square: {}'.format(my_1.shape)
+    assert my_2.shape[0] == my_2.shape[1], 'Size of eigenvector 2 not square: {}'.format(my_2.shape)
+
+    # Helper functions
+    def nominator(lam_1, my_1, lam_2, my_2):
+        return np.abs(my_1.T.dot(my_2)) / np.tile(np.abs(lam_1 + lam_2), (len(lam_1), 1))
+
+    def denom(lam, my):
+        return np.diag(my.conj().T.dot(my) / (2.0 * np.abs(np.real(lam))) + np.abs(my.T.dot(my)) / (2.0 * np.abs(lam)))
+
+    # Pre-compute the terms in the nominator and denominator
+    nomin1 = nominator(lam_1.conj(), my_1.conj(), lam_2, my_2)
+    nomin2 = nominator(lam_1, my_1, lam_2, my_2)
+    denom1 = denom(lam_1, my_1)
+    denom2 = denom(lam_2, my_2)
+    # Assemble everything
+    MACXP = (nomin1 + nomin2) ** 2.0 / (denom1 * denom2)
+    # The results should be completely real, but an imaginary part that is numerically zero remains and is discarded here.
+    MACXP = np.abs(MACXP)
+    # Optionally, visualize the results
+    if plot:
+        plt.figure()
+        plt.pcolor(MACXP, cmap='hot_r')
+        plt.colorbar()
+        plt.grid('on')
+        plt.title('MACXP')
+        return MACXP, plt
+    return MACXP
+
+
+def calc_PCC(lam_1, lam_2, plot=False):
+    """
+    This is the most simple pole correlation criterion I could think of. It calculated a value between 0.0 and 1.0 where
+    1.0 indicates identical poles. The cirterion works also for complex eigenvalues.
+    """
+    PCC = np.zeros((len(lam_1), len(lam_2)))
+    for jj in range(len(lam_2)):
+        # Calculate the delta with respect to all other poles.
+        delta = np.abs(lam_1 - lam_2[jj])
+        # Scale the values to the range of 0.0 to 1.0 and store in matrix.
+        PCC[:, jj] = 1.0 - delta / delta.max()
+    # Optionally, visualize the results
+    if plot:
+        plt.figure()
+        plt.pcolor(PCC, cmap='hot_r')
+        plt.colorbar()
+        plt.grid('on')
+        plt.title('PCC')
+        return PCC, plt
+    return PCC
+
+
 def force_matrix_symmetry(matrix):
     return (matrix + matrix.T) / 2.0
 

loadskernel/plotting_standard.py

@@ -440,9 +440,9 @@ def plot_eigenvalues(self):
             h = self.model['atmo'][trimcase['altitude']]['h'][()]
 
             # this kind of plot is only feasible for methods which iterate over Vtas, e.g. not the K- or KE-methods
-            if 'flutter' in simcase and simcase['flutter_para']['method'] not in ['pk', 'statespace']:
-                logging.info(
-                    'skip plotting of eigenvalues and -vectors for {}'.format(trimcase['desc']))
+            if 'flutter' in simcase and simcase['flutter_para']['method'] not in ['pk_rodden', 'pk_schwochow',
+                                                                                  'pk', 'statespace']:
+                logging.info('skip plotting of eigenvalues and -vectors for {}'.format(trimcase['desc']))
                 continue
 
             # Plot boundaries