Add structural damping via complex proportional stiffness method.

pypk/api_core.py

@@ -35,7 +35,8 @@ def init_pypk(cfg):
     else:
         raise RuntimeError('"fluid" entry must be "unmatched" or "matched_isa"!\n')
     # Eigen solver
-    eig = EigenSolver(cfg['l_ref'], cfg['n_modes'])
+    g = cfg['g_struct'] if 'g_struct' in cfg else 0.0
+    eig = EigenSolver(cfg['l_ref'], cfg['n_modes'], g)
     # Aggregator
     agg = KSAggregator(cfg['rho_ks'])
     # p-k method

pypk/eigensolver.py

@@ -20,9 +20,10 @@ import scipy.linalg as spla
 class EigenSolver:
     """Eigen solver
     """
-    def __init__(self, l, n):
+    def __init__(self, l, n, g=0):
         self._l = l # reference length
         self._n = n # number of modes
+        self._g = (1 + 1j * g) # structural damping (complex proportional stiffness)
 
     def compute(self, rho, u, m, k, q):
         """Compute eigenvalues and eigenvectors
@@ -32,7 +33,7 @@ class EigenSolver:
         k : stiffness matrix
         q : aerodynamic matrix
         """
-        p, v = spla.eig(k - 0.5 * rho * u**2 * q, -u**2 / self._l**2 * m) # det((p*u/l)^2*M + K - q_inf*Q) = 0
+        p, v = spla.eig(self._g * k - 0.5 * rho * u**2 * q, -u**2 / self._l**2 * m) # det((p*u/l)^2*M + K - q_inf*Q) = 0
         p = np.sqrt(p)
         # Sort
         p = np.where(np.imag(p) < 0, -p, p) # change sign of p if imag(p) < 0
@@ -58,10 +59,10 @@ class EigenSolver:
         a = np.zeros((self._n + 1, self._n + 1), dtype=complex)
         b = np.zeros((self._n + 1, 1), dtype=complex)
         a[:-1,0] = -2 * p * u**2 / self._l**2 * m @ v
-        a[:-1,1:] = p**2 * u**2 / self._l**2 * m + k - 0.5 * rho * u**2 * q
+        a[:-1,1:] = p**2 * u**2 / self._l**2 * m + self._g * k - 0.5 * rho * u**2 * q
         a[-1,0] = 0
         a[-1,1:] = v.T
-        b[:-1,0] = (p**2 * u**2 / self._l**2 * dm + dk - 0.5 * rho * u**2 * dq) @ v
+        b[:-1,0] = (p**2 * u**2 / self._l**2 * dm + self._g * dk - 0.5 * rho * u**2 * dq) @ v
         b[-1,0] = 0
         x = spla.lu_solve((spla.lu_factor(a)), b)
         return x[0,0]

tests/agard_grad.py

@@ -32,6 +32,7 @@ def get_cfg():
     'mu': 143.920, # mass ratio
     'f_ref': 40.3511, # frequency of first torsional mode
     'n_modes': 4, # number of modes
+    'g_struct': 0.02, # structural damping
     'rho_ks': 1e6, # KS aggregation parameter for KS
     'method': 'nipk', # method type
     'fluid': 'unmatched', # fluid type