Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove rigid body modes from modal analysis #1079

Merged
merged 2 commits into from
Jul 11, 2024
Merged

Remove rigid body modes from modal analysis #1079

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
09853c3
Remove rigid body modes from modal analysis
jguarato Jul 5, 2024
d82ba29
Fix test results using the eigen solver
jguarato Jul 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 19 additions & 7 deletions ross/rotor_assembly.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1327,6 +1327,13 @@ def _eigen(
----------
speed : float, pint.Quantity
Rotor speed. Default unit is rad/s.
num_modes : int, optional
The number of eigenvalues and eigenvectors to be calculated using ARPACK.
If sparse=True, it determines the number of eigenvalues and eigenvectors
to be calculated. It must be smaller than Rotor.ndof - 1. It is not
possible to compute all eigenvectors of a matrix with ARPACK.
If sparse=False, num_modes does not have any effect over the method.
Default is 12.
frequency: float, pint.Quantity
Excitation frequency. Default units is rad/s.
sorted_ : bool, optional
Expand Down Expand Up @@ -1372,15 +1379,20 @@ def _eigen(
try:
evalues, evectors = las.eigs(
A,
k=num_modes,
k=2 * num_modes,
sigma=1,
ncv=2 * num_modes,
ncv=4 * 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))

# Disregard rigid body modes:
idx = np.where(np.abs(evalues) > 0.1)[0]
evalues = evalues[idx]
evectors = evectors[:, idx]
except las.ArpackError:
evalues, evectors = la.eig(A)
else:
Expand Down Expand Up @@ -2500,14 +2512,14 @@ def run_ucs(
The bearing frequency range used to calculate the intersection points.
In some cases bearing coefficients will have to be extrapolated.
The default is None. In this case the bearing frequency attribute is used.
num : int
Number of steps in the range.
Default is 20.
num_modes : int, optional
Number of modes to be calculated. This uses scipy.sparse.eigs method.
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.
num : int
Number of steps in the range.
Default is 20.
synchronous : bool, optional
If True a synchronous analysis is carried out according to :cite:`rouch1980dynamic`.
Default is False.
Expand Down Expand Up @@ -4268,9 +4280,9 @@ def rotor_example_6dof():
>>> rotor_speed = 100.0 # rad/s
>>> modal6 = rotor6.run_modal(rotor_speed)
>>> print(f"Undamped natural frequencies: {np.round(modal6.wn, 2)}") # doctest: +ELLIPSIS
Undamped natural frequencies: [ 0. 47.62 ...
Undamped natural frequencies: [ 47.62 91.84 96.36 274.44 ...
>>> print(f"Damped natural frequencies: {np.round(modal6.wd, 2)}") # doctest: +ELLIPSIS
Damped natural frequencies: [ 0. 47.62 ...
Damped natural frequencies: [ 47.62 91.84 96.36 274.44 ...
>>> # Plotting Campbell Diagram
>>> camp6 = rotor6.run_campbell(np.linspace(0, 400, 101), frequencies=6)
>>> fig = camp6.plot()
Expand Down
18 changes: 9 additions & 9 deletions ross/tests/test_6dof_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,12 +94,12 @@ def test_run_modal(rotor_6dof):

# fmt:off
wn = np.array([
7.84532831e-05, 7.84532831e-05, 9.17864452e+01, 9.62950609e+01,
2.74045342e+02, 2.96967970e+02, 7.17111118e+02
91.78677 , 96.29604 , 274.059892, 296.994125, 717.35166 ,
770.203624, 774.349678
])
wd = np.array([
0. , 0. , 91.78644519, 96.29506086,
274.04534155, 296.9679697, 717.11111643
91.78677 , 96.29604 , 274.059892, 296.994125, 717.35166 ,
770.203624, 774.349678
])
# fmt:on

Expand All @@ -109,11 +109,11 @@ def test_run_modal(rotor_6dof):

def test_modal_results_equality(rotor_6dof, rotor_4dof):
speed = 100.0
modal_6dof = rotor_6dof.run_modal(speed, num_modes=14)
modal_6dof = rotor_6dof.run_modal(speed, num_modes=10)
modal_4dof = rotor_4dof.run_modal(speed, num_modes=10)

assert_allclose(modal_6dof.wn[2:], modal_4dof.wn, rtol=5e-2, atol=1)
assert_allclose(modal_6dof.wd[2:], modal_4dof.wd, rtol=5e-2, atol=1)
assert_allclose(modal_6dof.wn, modal_4dof.wn, rtol=5e-2, atol=1)
assert_allclose(modal_6dof.wd, modal_4dof.wd, rtol=5e-2, atol=1)


def test_campbell(rotor_6dof):
Expand All @@ -132,15 +132,15 @@ def test_campbell(rotor_6dof):
])
# fmt:on

assert_allclose(campbell.wd[:, 2], wd, rtol=1e-3)
assert_allclose(campbell.wd[:, 0], wd, rtol=1e-3)


def test_campbell_equality(rotor_6dof, rotor_4dof):
speed_range = np.linspace(315, 1150, 31)
campbell1 = rotor_6dof.run_campbell(speed_range)
campbell2 = rotor_4dof.run_campbell(speed_range)

assert_allclose(campbell1.wd[:, 2], campbell2.wd[:, 0], rtol=1e-3)
assert_allclose(campbell1.wd[:, 0], campbell2.wd[:, 0], rtol=1e-3)


def test_run_freq(rotor_6dof):
Expand Down
20 changes: 4 additions & 16 deletions ross/tests/test_rotor_assembly.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -314,12 +314,6 @@ def test_evals_sorted_rotor2(rotor2):
1.1024641658e-11 + 598.0247411456j,
4.3188161105e-09 + 3956.2249777612j,
2.5852376472e-11 + 3956.2249797838j,
4.3188161105e-09 - 3956.2249777612j,
2.5852376472e-11 - 3956.2249797838j,
7.4569772223e-11 - 598.0247411492j,
1.1024641658e-11 - 598.0247411456j,
1.4667459679e-12 - 215.3707255735j,
3.9623200168e-12 - 215.3707255733j,
]
)

Expand All @@ -331,24 +325,18 @@ def test_evals_sorted_rotor2(rotor2):
8.482603e-08 + 3470.897616j,
4.878990e-07 + 3850.212629j,
4.176291e01 + 3990.22903j,
4.176291e01 - 3990.22903j,
4.878990e-07 - 3850.212629j,
8.482603e-08 - 3470.897616j,
5.045245e-01 - 215.369011j,
-5.045245e-01 - 215.369011j,
-4.838034e-14 - 34.822138j,
]
)
modal2_0 = rotor2.run_modal(speed=0)
rotor2_evals, rotor2_evects = rotor2._eigen(speed=0)
assert_allclose(rotor2_evals, evals_sorted, rtol=1e-3)
assert_allclose(modal2_0.evalues, evals_sorted, rtol=1e-3)
assert_allclose(rotor2_evals[:6], evals_sorted, rtol=1e-3)
assert_allclose(modal2_0.evalues[:6], evals_sorted, rtol=1e-3)
modal2_10000 = rotor2.run_modal(speed=10000)
assert_allclose(modal2_10000.evalues, evals_sorted_w_10000, rtol=1e-1)
assert_allclose(modal2_10000.evalues[:6], evals_sorted_w_10000, rtol=1e-1)

# test run_modal with Q_
modal2_10000 = rotor2.run_modal(speed=Q_(95492.96585513721, "RPM"))
assert_allclose(modal2_10000.evalues, evals_sorted_w_10000, rtol=1e-1)
assert_allclose(modal2_10000.evalues[:6], evals_sorted_w_10000, rtol=1e-1)


@pytest.fixture
Expand Down
Loading