From 1de1eeaffb65fecf18fc36f2c42a11106ef693ae Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Mon, 13 May 2024 12:09:17 +0200
Subject: [PATCH] Add structural damping via complex proportional stiffness
method.
---
pypk/api_core.py | 3 ++-
pypk/eigensolver.py | 9 +++++----
tests/agard_grad.py | 1 +
3 files changed, 8 insertions(+), 5 deletions(-)
diff --git a/pypk/api_core.py b/pypk/api_core.py
index d5bab2b..4787eb1 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 d0fe1eb..9ff0e6e 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 e44e5c4..77ddfdf 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
--
GitLab