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', 'MAC*PCC' (recommended), 'MAC*HDM'
+                         # 'flutter_para': {'method': 'pk', 'Vtas': np.linspace(100.0, 500.0, 100),
+                         #                  'tracking': 'MAC*PCC'},
                          },
                         ]
         # End

doc/tutorials/DC3_model/JCLs/jcl_dc3_flutter.py

@@ -239,5 +239,5 @@ def __init__(self):
                          # True or False, enables flutter check with k, ke or pk method
                          'flutter': True,
                          # flutter parameters for pk method
-                         'flutter_para': {'method': 'pk', 'Vtas': np.linspace(20.0, 300.0, 20)}}]
+                         'flutter_para': {'method': 'pk', 'Vtas': np.linspace(20.0, 300.0, 20), '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,35 @@ def eval_equations(self):
                     }
         return response
 
-    def calc_eigenvalues(self, A, eigenvector_old):
+    def calc_eigenvalues(self, A, eigenvalues_old, eigenvectors_old):
+        # Find all eigenvalues and eigenvectors
         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)
-        idx_pos = self.get_best_match(MAC)
-        eigenvalues = eigenvalue[idx_pos]  # [idx_sort]
-        eigenvectors = eigenvector[:, idx_pos]  # [:, idx_sort]
+        # To match the modes with the previous step, use a correlation cirterion as specified in the JCL.
+        if 'tracking' not in self.simcase['flutter_para']:
+            # Set a default.
+            tracking_method = 'MAC'
+        else:
+            tracking_method = self.simcase['flutter_para']['tracking']
+        # Calculate the correlation bewteen the old and current modes.
+        if tracking_method == 'MAC':
+            # Most simple, use only the modal assurance criterion (MAC).
+            correlation = fem_helper.calc_MAC(eigenvectors_old, eigenvector)
+        elif tracking_method == 'MAC*PCC':
+            # Combining MAC and pole correlation cirterion (PCC) for improved handling of complex conjugate pairs.
+            correlation = fem_helper.calc_MAC(eigenvectors_old, eigenvector) * fem_helper.calc_PCC(eigenvalues_old, eigenvalue)
+        elif tracking_method == 'MAC*HDM':
+            # Combining MAC and hyperboloic distance metric (HDM) for improved handling of complex conjugate pairs.
+            correlation = fem_helper.calc_MAC(eigenvectors_old, eigenvector) * fem_helper.calc_HDM(eigenvalues_old, eigenvalue)
+        # Based on the correlation matrix, find the best match and apply to the modes.
+        idx_pos = self.get_best_match(correlation)
+        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 +675,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 +699,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,67 @@ 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')
-
-        return MAC, plt
+        plot_correlation_matrix(MAC, 'MAC')
     return MAC
 
 
+def calc_PCC(lam1, lam2, plot=False):
+    """
+    This is the most simple pole correlation criterion I could think of. It calculates a value between 0.0 and 1.0 where
+    1.0 indicates identical poles. What is important is that the criterion works for complex eigenvalues.
+    """
+    PCC = np.zeros((len(lam1), len(lam2)))
+    for jj in range(len(lam2)):
+        # Calculate the delta with respect to all other poles.
+        delta = np.abs(lam1 - lam2[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:
+        plot_correlation_matrix(PCC, 'PCC')
+    return PCC
+
+
+def calc_HDM(lam1, lam2, plot=False):
+    """
+    Hyperbolic distance metric. Proposed by [1], implementation and application shown in [2], section 6.2.5.
+
+    [1] Luspay, T., Péni, T., Gőzse, I., Szabó, Z., and Vanek, B., “Model reduction for LPV systems based on approximate modal
+    decomposition”, International Journal for Numerical Methods in Engineering, vol. 113, no. 6, pp. 891–909, 2018,
+    https://doi.org/10.1002/nme.5692.
+    [2] Jelicic, G., “System Identification of Parameter-Varying Aeroelastic Systems using Real-Time Operational Modal
+    Analysis”, Deutsches Zentrum für Luft- und Raumfahrt  e. V., 2022, https://doi.org/10.57676/P9QV-CK92.
+    """
+    HDM = np.zeros((len(lam1), len(lam2)))
+    # Map to unit circle.
+    fs = 2.56 * np.max(np.abs(np.concatenate((lam1, lam1))) / 2.0 / np.pi)
+    z1 = np.exp(lam1 / fs)
+    z2 = np.exp(lam2 / fs)
+    # Mirror unstable poles about unit circle.
+    pos1 = lam1.real > 0.0
+    pos2 = lam2.real > 0.0
+    z1[pos1] = 1.0 / z1[pos1].conj()
+    z2[pos2] = 1.0 / z2[pos2].conj()
+    # Implementation of equations (4) and (A4) in [1] per row
+    for jj in range(len(lam2)):
+        HDM[:, jj] = 1.0 - np.abs((z1 - z2[jj]) / (1.0 - z1 * z2[jj].conj()))
+    # Optionally, visualize the results
+    if plot:
+        plot_correlation_matrix(HDM, 'HDM')
+    return HDM
+
+
+def plot_correlation_matrix(matrix, name='Correlation Matrix'):
+    plt.figure()
+    plt.pcolor(matrix, cmap='hot_r')
+    plt.colorbar()
+    plt.grid('on')
+    plt.title(name)
+    plt.show()
+
+
 def force_matrix_symmetry(matrix):
     return (matrix + matrix.T) / 2.0