Skip to content

Commit

Permalink
Merge pull request #931 from fernandarossi/dev-sync_ucs_map
Browse files Browse the repository at this point in the history 
Dev sync ucs map
  • Loading branch information
raphaeltimbo authored Apr 11, 2024
2 parents 1cc2015 + e69d577 commit d9f7716
Show file tree
Hide file tree
Showing 4 changed files with 685 additions and 113 deletions.
11 changes: 11 additions & 0 deletions docs/references.bib
View file
Original file line number Diff line number Diff line change
Expand Up @@ -720,3 +720,14 @@ @article{genta1988conical
journal = {Journal of Sound and Vibration},
doi = {10.1016/0022-460X(88)90342-2}
}

@article{rouch1980dynamic,
author = {Rouch, K. E. and Kao, J. S.},
title = {Dynamic Reduction in Rotor Dynamics by the Finite Element Method},
journal = {Journal of Mechanical Design},
volume = {102},
number = {2},
pages = {360-368},
year = {1980},
doi = {10.1115/1.3254752}
}
27 changes: 15 additions & 12 deletions ross/disk_element.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,18 +90,21 @@ def __eq__(self, other):
True
"""
false_number = 0
for i in self.__dict__:
try:
if np.allclose(self.__dict__[i], other.__dict__[i]):
pass
else:
false_number += 1

except TypeError:
if self.__dict__[i] == other.__dict__[i]:
pass
else:
false_number += 1
if "DiskElement" in other.__class__.__name__:
for i in self.__dict__:
try:
if np.allclose(self.__dict__[i], other.__dict__[i]):
pass
else:
false_number += 1

except TypeError:
if self.__dict__[i] == other.__dict__[i]:
pass
else:
false_number += 1
else:
false_number += 1

if false_number == 0:
return True
Expand Down
151 changes: 100 additions & 51 deletions ross/rotor_assembly.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -583,7 +583,7 @@ def __eq__(self, other):
return False

@check_units
def run_modal(self, speed, num_modes=12, sparse=True):
def run_modal(self, speed, num_modes=12, sparse=True, synchronous=False):
"""Run modal analysis.
Method to calculate eigenvalues and eigvectors for a given rotor system.
Expand Down Expand Up @@ -615,6 +615,9 @@ def run_modal(self, speed, num_modes=12, sparse=True):
If False, scipy.linalg.eig() is used to calculate all the eigenvalues and
eigenvectors.
Default is True.
synchronous : bool, optional
If True a synchronous analysis is carried out according to :cite:`rouch1980dynamic`.
Default is False.
Returns
-------
Expand All @@ -638,7 +641,9 @@ def run_modal(self, speed, num_modes=12, sparse=True):
>>> mode2 = 1 # Second mode
>>> fig = modal.plot_mode_2d(mode2)
"""
evalues, evectors = self._eigen(speed, num_modes=num_modes, sparse=sparse)
evalues, evectors = self._eigen(
speed, num_modes=num_modes, sparse=sparse, synchronous=synchronous
)
wn_len = num_modes // 2
wn = (np.absolute(evalues))[:wn_len]
wd = (np.imag(evalues))[:wn_len]
Expand Down Expand Up @@ -887,9 +892,15 @@ def convergence(self, n_eigval=0, err_max=1e-02):

return results

def M(self, frequency=None):
def M(self, frequency=None, synchronous=False):
"""Mass matrix for an instance of a rotor.
Parameters
----------
synchronous : bool, optional
If True a synchronous analysis is carried out.
Default is False.
Returns
-------
M0 : np.ndarray
Expand All @@ -913,10 +924,44 @@ def M(self, frequency=None):

for elm in self.elements:
dofs = list(elm.dof_global_index.values())
try:
M0[np.ix_(dofs, dofs)] += elm.M(frequency)
except TypeError:
M0[np.ix_(dofs, dofs)] += elm.M()
if synchronous:
if elm in self.shaft_elements:
try:
M = elm.M(frequency)
except TypeError:
M = elm.M()
G = elm.G()
for i in range(8):
if i == 0 or i == 3 or i == 4 or i == 7:
M[i, 0] = M[i, 0] - G[i, 1]
M[i, 3] = M[i, 3] + G[i, 2]
M[i, 4] = M[i, 4] - G[i, 5]
M[i, 7] = M[i, 7] + G[i, 6]
else:
M[i, 1] = M[i, 1] + G[i, 0]
M[i, 2] = M[i, 2] - G[i, 3]
M[i, 5] = M[i, 5] + G[i, 4]
M[i, 6] = M[i, 6] - G[i, 7]
M0[np.ix_(dofs, dofs)] += M
elif elm in self.disk_elements:
try:
M = elm.M(frequency)
except TypeError:
M = elm.M()
G = elm.G()
M[2, 2] = M[2, 2] - G[2, 3]
M[3, 3] = M[3, 3] + G[3, 2]
M0[np.ix_(dofs, dofs)] += M
else:
try:
M0[np.ix_(dofs, dofs)] += elm.M(frequency)
except TypeError:
M0[np.ix_(dofs, dofs)] += elm.M()
else:
try:
M0[np.ix_(dofs, dofs)] += elm.M(frequency)
except TypeError:
M0[np.ix_(dofs, dofs)] += elm.M()

return M0

Expand Down Expand Up @@ -1053,7 +1098,7 @@ def G(self):

return G0

def A(self, speed=0, frequency=None):
def A(self, speed=0, frequency=None, synchronous=False):
"""State space matrix for an instance of a rotor.
Parameters
Expand All @@ -1063,6 +1108,9 @@ def A(self, speed=0, frequency=None):
Default is 0.
frequency : float, optional
Excitation frequency. Default is rotor speed.
synchronous : bool, optional
If True a synchronous analysis is carried out.
Default is False.
Returns
-------
Expand All @@ -1089,7 +1137,7 @@ def A(self, speed=0, frequency=None):
# fmt: off
A = np.vstack(
[np.hstack([Z, I]),
np.hstack([la.solve(-self.M(frequency), self.K(frequency)), la.solve(-self.M(frequency), (self.C(frequency) + self.G() * speed))])])
np.hstack([la.solve(-self.M(frequency, synchronous=synchronous), self.K(frequency)), la.solve(-self.M(frequency,synchronous=synchronous), (self.C(frequency) + self.G() * speed))])])
# fmt: on

return A
Expand Down Expand Up @@ -1240,7 +1288,14 @@ def _index(eigenvalues):

@check_units
def _eigen(
self, speed, num_modes=12, frequency=None, sorted_=True, A=None, sparse=True
self,
speed,
num_modes=12,
frequency=None,
sorted_=True,
A=None,
sparse=True,
synchronous=False,
):
"""Calculate eigenvalues and eigenvectors.
Expand All @@ -1264,6 +1319,9 @@ def _eigen(
sparse : bool, optional
If sparse, eigenvalues will be calculated with arpack.
Default is True.
synchronous : bool, optional
If True a synchronous analysis is carried out.
Default is False.
Returns
-------
Expand All @@ -1280,20 +1338,34 @@ def _eigen(
91.796...
"""
if A is None:
A = self.A(speed=speed, frequency=frequency)
A = self.A(speed=speed, frequency=frequency, synchronous=synchronous)

if sparse is True:
try:
evalues, evectors = las.eigs(
A, k=num_modes, sigma=0, ncv=2 * num_modes, which="LM", v0=self._v0
)
# store v0 as a linear combination of the previously
# calculated eigenvectors to use in the next call to eigs
self._v0 = np.real(sum(evectors.T))
except las.ArpackError:
evalues, evectors = la.eig(A)
else:
if synchronous:
evalues, evectors = la.eig(A)
idx = np.where(np.imag(evalues) != 0)[0]
evalues = evalues[idx]
evectors = evectors[:, idx]
idx = np.where(np.abs(np.real(evalues) / np.imag(evalues)) < 1000)[0]
evalues = evalues[idx]
evectors = evectors[:, idx]
else:
if sparse is True:
try:
evalues, evectors = las.eigs(
A,
k=num_modes,
sigma=0,
ncv=2 * num_modes,
which="LM",
v0=self._v0,
)
# store v0 as a linear combination of the previously
# calculated eigenvectors to use in the next call to eigs
self._v0 = np.real(sum(evectors.T))
except las.ArpackError:
evalues, evectors = la.eig(A)
else:
evalues, evectors = la.eig(A)

if sorted_ is False:
return evalues, evectors
Expand Down Expand Up @@ -2176,9 +2248,8 @@ def run_ucs(
Default is 16. In this case 4 modes are plotted, since for each pair
of eigenvalues calculated we have one wn, and we show only the
forward mode in the plots.
synchronous : bool
If True a synchronous analysis is carried out and the frequency of
the first forward model will be equal to the speed.
synchronous : bool, optional
If True a synchronous analysis is carried out.
Default is False.
Returns
Expand Down Expand Up @@ -2214,32 +2285,10 @@ def run_ucs(
for i, k in enumerate(stiffness_log):
bearings = [b.__class__(b.n, kxx=k, cxx=0) for b in bearings_elements]
rotor = self.__class__(self.shaft_elements, self.disk_elements, bearings)
speed = 0
if synchronous:

def wn_diff(x):
"""Function to evaluate difference between speed and
natural frequency for the first mode."""
modal = rotor.run_modal(speed=x, num_modes=num_modes)
# get first forward mode
if modal.whirl_direction()[0] == "Forward":
wn0 = modal.wn[0]
else:
wn0 = modal.wn[1]

return wn0 - x

speed = newton(wn_diff, 0)
modal = rotor.run_modal(speed=speed, num_modes=num_modes)

# if sync, select only forward modes
if synchronous:
rotor_wn[:, i] = modal.wn[modal.whirl_direction() == "Forward"]
# if not sync, with speed=0 whirl direction can be confusing, with
# two close modes being forward or backward, so we select one mode in
# each 2 modes.
else:
rotor_wn[:, i] = modal.wn[::2]
modal = rotor.run_modal(
speed=0, num_modes=num_modes, synchronous=synchronous
)
rotor_wn[:, i] = modal.wn[::2]

bearing0 = bearings_elements[0]

Expand Down
Loading

0 comments on commit d9f7716

Please sign in to comment.