Merge pull request #1216 from Raimundovpn/frf

Frf

Author: raphaeltimbo

ross/results.py

@@ -6,20 +6,20 @@
 import copy
 import inspect
 from abc import ABC
-from prettytable import PrettyTable
 from collections.abc import Iterable
 from warnings import warn
 
 import numpy as np
 import pandas as pd
 import toml
+from numba import njit
+from numpy import linalg as la
 from plotly import graph_objects as go
 from plotly.subplots import make_subplots
+from prettytable import PrettyTable
 from scipy.fft import fft
-from numpy import linalg as la
-from numba import njit
 
-from ross.plotly_theme import tableau_colors, coolwarm_r
+from ross.plotly_theme import coolwarm_r, tableau_colors
 from ross.units import Q_, check_units
 from ross.utils import intersection
 
@@ -286,19 +286,16 @@ def calculate_amplitude(self, angle):
             Tuple with (amplitude, phase) value.
             The amplitude units are the same as the ru_e and rv_e used to create the orbit.
         """
-        # find closest angle index
+
         if angle == "major":
             return self.major_axis, self.major_angle
         elif angle == "minor":
             return self.minor_axis, self.minor_angle
 
-        idx = (np.abs(self.angle - angle)).argmin()
-        amplitude = np.sqrt(self.x_circle[idx] ** 2 + self.y_circle[idx] ** 2)
-        phase = self.angle[0] + angle
-        if phase > 2 * np.pi:
-            phase -= 2 * np.pi
-
-        return amplitude, phase
+        amp_complex = self.ru_e * np.cos(angle) + self.rv_e * np.sin(angle)
+        amp_abs = np.abs(self.ru_e * np.cos(angle) + self.rv_e * np.sin(angle))
+        phase = 2 * np.pi - np.mod(np.angle(amp_complex), 2 * np.pi)
+        return amp_abs, phase
 
     def plot_orbit(self, fig=None):
         if fig is None:
@@ -359,7 +356,7 @@ def _init_orbit(ru_e, rv_e):
     nv = nv
     # fmt: off
     T = np.array([[ru * np.cos(nu), -ru * np.sin(nu)],
-                    [rv * np.cos(nv), -rv * np.sin(nv)]])
+                  [rv * np.cos(nv), -rv * np.sin(nv)]])
     # fmt: on
     H = T @ T.T
 
@@ -2625,8 +2622,8 @@ def plot_with_mode_shape(
     ):
         try:
             import random
-            from dash import Dash
-            from dash import dcc, html
+
+            from dash import Dash, dcc, html
             from dash.dependencies import Input, Output
         except ImportError:
             raise ImportError("Please install dash to use this feature.")
@@ -3305,19 +3302,39 @@ def data_magnitude(
                 try:
                     probe_tag = p[2]
                 except IndexError:
-                    probe_tag = f"Probe {i+1} - Node {p[0]}"
+                    probe_tag = f"Probe {i + 1} - Node {p[0]}"
 
             amplitude = []
             for speed_idx in range(len(self.speed_range)):
-                ru_e, rv_e = response[:, speed_idx][
-                    self.rotor.number_dof * node : self.rotor.number_dof * node + 2
-                ]
-                orbit = Orbit(
-                    node=node, node_pos=self.rotor.nodes_pos[node], ru_e=ru_e, rv_e=rv_e
-                )
-                amp, phase = orbit.calculate_amplitude(angle=angle)
-                amplitude.append(amp)
+                if node not in self.rotor.link_nodes:
+                    ru_e, rv_e = response[:, speed_idx][
+                        self.rotor.number_dof * node : self.rotor.number_dof * node + 2
+                    ]
+                    orbit = Orbit(
+                        node=node,
+                        node_pos=self.rotor.nodes_pos[node],
+                        ru_e=ru_e,
+                        rv_e=rv_e,
+                    )
+                    amp, phase = orbit.calculate_amplitude(angle=angle)
+
+                else:
+                    position_in_link = node - len(self.rotor.nodes_pos)
+                    start_index = len(self.rotor.nodes_pos) * self.rotor.number_dof + (
+                        position_in_link * 3
+                    )
+                    ru_e, rv_e = response[:, speed_idx][start_index : start_index + 2]
 
+                    orbit = Orbit(
+                        node=None,
+                        node_pos=None,
+                        ru_e=ru_e,
+                        rv_e=rv_e,
+                    )
+
+                    amp, phase = orbit.calculate_amplitude(angle=angle)
+
+                amplitude.append(amp)
             data[probe_tag] = Q_(amplitude, base_unit).to(amplitude_units).m
 
         df = pd.DataFrame(data)
@@ -3400,21 +3417,42 @@ def data_phase(
                 try:
                     probe_tag = p[2]
                 except IndexError:
-                    probe_tag = f"Probe {i+1} - Node {p[0]}"
+                    probe_tag = f"Probe {i + 1} - Node {p[0]}"
 
             phase_values = []
             for speed_idx in range(len(self.speed_range)):
-                ru_e, rv_e = response[:, speed_idx][
-                    self.rotor.number_dof * node : self.rotor.number_dof * node + 2
-                ]
-                orbit = Orbit(
-                    node=node, node_pos=self.rotor.nodes_pos[node], ru_e=ru_e, rv_e=rv_e
-                )
-                amp, phase = orbit.calculate_amplitude(angle=angle)
+                if node not in self.rotor.link_nodes:
+                    ru_e, rv_e = response[:, speed_idx][
+                        self.rotor.number_dof * node : self.rotor.number_dof * node + 2
+                    ]
+                    orbit = Orbit(
+                        node=node,
+                        node_pos=self.rotor.nodes_pos[node],
+                        ru_e=ru_e,
+                        rv_e=rv_e,
+                    )
+
+                    amp, phase = orbit.calculate_amplitude(angle=angle)
+
+                else:
+                    position_in_link = node - len(self.rotor.nodes_pos)
+                    start_index = len(self.rotor.nodes_pos) * self.rotor.number_dof + (
+                        position_in_link * 3
+                    )
+                    ru_e, rv_e = response[:, speed_idx][start_index : start_index + 2]
+
+                    orbit = Orbit(
+                        node=None,
+                        node_pos=None,
+                        ru_e=ru_e,
+                        rv_e=rv_e,
+                    )
+
+                    amp, phase = orbit.calculate_amplitude(angle=angle)
+
                 phase_values.append(phase)
 
             data[probe_tag] = Q_(phase_values, "rad").to(phase_units).m
-
         df = pd.DataFrame(data)
 
         return df
@@ -5454,7 +5492,7 @@ def data_time_response(
                 try:
                     probe_tag = p[2]
                 except IndexError:
-                    probe_tag = f"Probe {i+1} - Node {p[0]}"
+                    probe_tag = f"Probe {i + 1} - Node {p[0]}"
 
             data[f"angle[{i}]"] = angle
             data[f"probe_tag[{i}]"] = probe_tag

ross/rotor_assembly.py

@@ -5,32 +5,33 @@
 from copy import copy, deepcopy
 from itertools import chain, cycle
 from pathlib import Path
-from methodtools import lru_cache
 
 import numpy as np
 import pandas as pd
 import toml
+from methodtools import lru_cache
 from plotly import express as px
 from plotly import graph_objects as go
 from scipy import io as sio
 from scipy import linalg as la
 from scipy import signal as signal
+from scipy.integrate import cumulative_trapezoid as integrate
+from scipy.linalg import lu_factor, lu_solve
 from scipy.optimize import newton
 from scipy.sparse import linalg as las
-from scipy.integrate import cumulative_trapezoid as integrate
 
 from ross.bearing_seal_element import (
     BallBearingElement,
     BearingElement,
     BearingFluidFlow,
+    CylindricalBearing,
     MagneticBearingElement,
     RollerBearingElement,
     SealElement,
-    CylindricalBearing,
 )
-from ross.faults import Crack, MisalignmentFlex, MisalignmentRigid, Rubbing
-from ross.disk_element import DiskElement
 from ross.coupling_element import CouplingElement
+from ross.disk_element import DiskElement
+from ross.faults import Crack, MisalignmentFlex, MisalignmentRigid, Rubbing
 from ross.materials import Material, steel
 from ross.point_mass import PointMass
 from ross.results import (
@@ -46,17 +47,17 @@
     TimeResponseResults,
     UCSResults,
 )
+from ross.seals.labyrinth_seal import LabyrinthSeal
 from ross.shaft_element import ShaftElement
 from ross.units import Q_, check_units
 from ross.utils import (
-    intersection,
-    newmark,
     assemble_C_K_matrices,
-    remove_dofs,
     convert_6dof_to_4dof,
     convert_6dof_to_torsional,
+    intersection,
+    newmark,
+    remove_dofs,
 )
-from ross.seals.labyrinth_seal import LabyrinthSeal
 
 __all__ = [
     "Rotor",
@@ -1666,7 +1667,7 @@ def _pseudo_modal(self, speed, num_modes):
 
         return matrix_to_modal, vector_to_modal, vector_from_modal
 
-    def transfer_matrix(self, speed=None, frequency=None, modes=None):
+    def transfer_matrix(self, speed=None, frequency=None):
         """Calculate the fer matrix for the frequency response function (FRF).
 
         Paramenters
@@ -1675,9 +1676,6 @@ def transfer_matrix(self, speed=None, frequency=None, modes=None):
             Excitation frequency. Default is rotor speed.
         speed : float, optional
             Rotating speed. Default is rotor speed (frequency).
-        modes : list, optional
-            List with modes used to calculate the matrix.
-            (all modes will be used if a list is not given).
 
         Returns
         -------
@@ -1693,33 +1691,17 @@ def transfer_matrix(self, speed=None, frequency=None, modes=None):
         if frequency is None:
             frequency = speed
 
-        lti = self._lti(speed=speed)
-        B = lti.B
-        C = lti.C
-        D = lti.D
+        I = np.eye(self.M().shape[0])
 
-        # calculate eigenvalues and eigenvectors using la.eig to get
-        # left and right eigenvectors.
-        evals, psi = self._eigen(speed=speed, frequency=frequency)
-
-        psi_inv = la.inv(psi)
+        lu, piv = lu_factor(
+            -(frequency**2) * self.M(frequency=frequency)
+            + 1j * frequency * (self.C(frequency=frequency) + frequency * self.G())
+            + self.K(frequency=frequency)
+        )
+        H = lu_solve((lu, piv), I)
 
-        if modes is not None:
-            n = self.ndof  # n dof -> number of modes
-            m = len(modes)  # -> number of desired modes
-            # idx to get each evalue/evector and its conjugate
-            idx = np.zeros((2 * m), int)
-            idx[0:m] = modes  # modes
-            idx[m:] = range(2 * n)[-m:]  # conjugates (see how evalues are ordered)
-            evals = evals[np.ix_(idx)]
-            psi = psi[np.ix_(range(2 * n), idx)]
-            psi_inv = psi_inv[np.ix_(idx, range(2 * n))]
-
-        with np.errstate(divide="ignore", invalid="ignore"):
-            diag = np.diag([1 / (1j * frequency - lam) for lam in evals])
-            diag[np.isnan(diag)] = 0
-
-        H = C @ psi @ diag @ psi_inv @ B + D
+        if np.isnan(H).any():
+            H = np.zeros((H.shape))
 
         return H
 
@@ -1788,12 +1770,12 @@ def run_freq_response(
         --------
         >>> import ross as rs
         >>> rotor = rs.rotor_example()
-        >>> speed = np.linspace(0, 1000, 101)
+        >>> speed =np.linspace(0, 1000, 101)
         >>> response = rotor.run_freq_response(speed_range=speed)
 
         Return the response amplitude
         >>> abs(response.freq_resp) # doctest: +ELLIPSIS
-        array([[[1.00000000e-06, 1.00261725e-06, 1.01076952e-06, ...
+        array([[[0.00000000e+00, 1.00261725e-06, 1.01076952e-06, ...
 
         Return the response phase
         >>> np.angle(response.freq_resp) # doctest: +ELLIPSIS
@@ -1810,7 +1792,7 @@ def run_freq_response(
         Selecting the disirable modes, if you want a reduced model:
         >>> response = rotor.run_freq_response(speed_range=speed, modes=[0, 1, 2, 3, 4])
         >>> abs(response.freq_resp) # doctest: +ELLIPSIS
-        array([[[2.00154633e-07, 2.02422522e-07, 2.09522044e-07, ...
+        array([[[0.00000000e+00, 1.00261725e-06, 1.01076952e-06, ...
 
         Plotting frequency response function:
         >>> fig = response.plot(inp=13, out=13)
@@ -1877,14 +1859,12 @@ def _run_freq_response(
         accl_resp = np.empty((self.ndof, self.ndof, len(speed_range)), dtype=complex)
 
         if free_free:
-            transfer_matrix = lambda s, m: self.transfer_matrix(
-                speed=0, modes=m, frequency=s
-            )
+            transfer_matrix = lambda s: self.transfer_matrix(speed=0, frequency=s)
         else:
-            transfer_matrix = lambda s, m: self.transfer_matrix(speed=s, modes=m)
+            transfer_matrix = lambda s: self.transfer_matrix(speed=s)
 
         for i, speed in enumerate(speed_range):
-            H = transfer_matrix(speed, modes)
+            H = transfer_matrix(speed)
             freq_resp[..., i] = H
             velc_resp[..., i] = 1j * speed * H
             accl_resp[..., i] = -(speed**2) * H
@@ -2006,7 +1986,6 @@ def run_forced_response(
             forced_resp[:, i] = freq_resp.freq_resp[..., i] @ force[..., i]
             velc_resp[:, i] = freq_resp.velc_resp[..., i] @ force[..., i]
             accl_resp[:, i] = freq_resp.accl_resp[..., i] @ force[..., i]
-
         forced_resp = ForcedResponseResults(
             rotor=self,
             forced_resp=forced_resp,
@@ -2210,7 +2189,7 @@ def run_unbalance_response(
 
         Return the response phase
         >>> np.angle(response.forced_resp) # doctest: +ELLIPSIS
-        array([[ 0.00000000e+00, ...
+        array([[ 0.        ,  0.        ,  0.        , ...
 
         Using clustered points option.
         Set `cluster_points=True` and choose how many modes the method must search and

ross/tests/test_results.py

@@ -195,7 +195,7 @@ def test_orbit_calculate_amplitude():
     # create orbit with major axis at 45deg
     orb = Orbit(node=0, node_pos=0, ru_e=(2 + 1j), rv_e=(2 - 1j))
 
-    assert_allclose(orb.calculate_amplitude(Q_(0, "deg"))[0], 1.7949401413591568)
+    assert_allclose(orb.calculate_amplitude(Q_(0, "deg"))[0], 2.23606797749979)
     assert_allclose(orb.calculate_amplitude(Q_(45, "deg"))[0], 2.8284271247461903)
     assert_allclose(
         orb.calculate_amplitude(Q_(135, "deg"))[0], 1.4142135623730947, rtol=1e-3