Dev sync ucs map

docs/references.bib

@@ -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}
+}

ross/disk_element.py

@@ -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

ross/rotor_assembly.py

@@ -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.
@@ -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
         -------
@@ -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]
@@ -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
@@ -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
 
@@ -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
@@ -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
         -------
@@ -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
@@ -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.
 
@@ -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
         -------
@@ -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
@@ -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
@@ -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]