diff --git a/pypk/api_core.py b/pypk/api_core.py
index d5bab2b2c88f8426f6b62935965d1bf79e609c3a..4787eb10579d712050fa867719dad899a1688cf0 100644
--- a/pypk/api_core.py
+++ b/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
diff --git a/pypk/eigensolver.py b/pypk/eigensolver.py
index d0fe1eb8f684f27461cf142f9fa3971ef1e0781b..9ff0e6e787ce202567799d10ab85244b85e5654a 100644
--- a/pypk/eigensolver.py
+++ b/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]
diff --git a/tests/agard_grad.py b/tests/agard_grad.py
index e44e5c495abab311e0ba21c86532839b052275f3..77ddfdfdac2f83f626d19f3602b84849dd8f4416 100644
--- a/tests/agard_grad.py
+++ b/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