Repository: ross-rotordynamics/ross
Pull request #278 (2019-09-05 17:48:02)

Add separated DataFrame to seals
This splits bearings from seals in two different dataframes:
* df_bearings for bearings
* df_seals for seals

It aims to better organize the elements and to avoid code interpretation problems such as in static analysis:

This is the free-body diagram of an example while seals and bearings are merged in same dataframe. The code interprets seals like bearings and calculate reaction forces where both seals are placed:
![antigo](https://user-images.githubusercontent.com/45969994/64347715-c6f81000-cfca-11e9-83b1-1acc0cb47c90.PNG)

And now, spliting bearings and seals:
![novo](https://user-images.githubusercontent.com/45969994/64347738-d2e3d200-cfca-11e9-88fb-3e74144f45a4.PNG)

Also, fix a bug when modeling too flexible bearings, displacement were getting larger than acceptable.
It uses axuliar bearings with high stiffness, considering almost zero displacement in bearing nodes.



Patch:
diff --git a/ross/api_report.py b/ross/api_report.py
index 4a1f901b..1f3026c5 100644
--- a/ross/api_report.py
+++ b/ross/api_report.py
@@ -1,7 +1,7 @@
 import numpy as np
 from scipy.signal import argrelextrema
 
-from ross.rotor_assembly import Rotor, rotor_example
+from rotor_assembly import Rotor
 import ross as rs
 
 import bokeh.palettes as bp
@@ -42,7 +42,7 @@ class Report:
         speed_units : str
             String defining the unit for rotor speed.
             Default is "rpm".
-
+            
         Return
         ------
 
@@ -50,20 +50,19 @@ class Report:
         -------
         >>> rotor = rotor_example()
         >>> report = Report(rotor=rotor,
-        ...                 minspeed=400,
-        ...                 maxspeed=1000,
-        ...                 speed_units="rad/s")
+                            minspeed=400,
+                            maxspeed=1000) # doctest: +ELLIPSIS
+        [...]
         """
         self.rotor = rotor
+        self.maxspeed = maxspeed
+        self.minspeed = minspeed
         self.speed_units = speed_units
 
         if speed_units == "rpm":
             minspeed = minspeed * np.pi / 30
             maxspeed = maxspeed * np.pi / 30
 
-        self.maxspeed = maxspeed
-        self.minspeed = minspeed
-                
     @classmethod
     def from_saved_rotors(cls, path):
         rotor = rs.Rotor.load(path)
@@ -81,12 +80,9 @@ class Report:
         Example
         -------
         >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor,
-        ...                 minspeed=400,
-        ...                 maxspeed=1000,
-        ...                 speed_units="rad/s")
-        >>> report.static_forces()
-        array([44.09320349, 44.09320349])
+        >>> report = Report(rotor=rotor, minspeed=400, maxspeed=1000)
+        >>> report.static_forces() # doctest: +ELLIPSIS
+        [...]
         """
         # get reaction forces on bearings
         Fb = self.rotor.run_static().force_data['Bearings Reaction Forces']
@@ -111,12 +107,9 @@ class Report:
         Example
         -------
         >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor,
-        ...                 minspeed=400,
-        ...                 maxspeed=1000,
-        ...                 speed_units="rad/s")
-        >>> report.unbalance_forces(mode=0)
-        [234.5654171598467]
+        >>> report = Report(rotor=rotor, minspeed=400, maxspeed=1000)
+        >>> report.unbalance_forces(mode=0) # doctest: +ELLIPSIS
+        [...]
         """
 
         N = 60 * self.maxspeed / (2 * np.pi)
@@ -168,16 +161,13 @@ class Report:
 
         Example
         -------
-        """
-        """
         >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor,
-        ...                 minspeed=400,
-        ...                 maxspeed=1000,
-        ...                 speed_units="rad/s")
-        >>> clearances = {3:0.001, 5:0.002}
-        >>> report.unbalance_response(clearances=clearances, mode=0)
-        Figure(id='1372', ...)
+        >>> report = Report(rotor=rotor)
+        >>> clearances = dict(3=0.001, 5=0.002)
+        >>> report.unbalance_response(clearances=clearances,
+                                      mode=0,
+                                      ) # doctest: +ELLIPSIS
+        [...]
         """
         maxspeed = self.maxspeed
         minspeed = self.minspeed
@@ -391,12 +381,9 @@ class Report:
         Example
         -------
         >>> rotor = rotor_example()
-        >>> report = Report(rotor=rotor,
-        ...                 minspeed=400,
-        ...                 maxspeed=1000,
-        ...                 speed_units="rad/s")
-        >>> report.mode_shape(mode=0)
-        ([], array([3.]))
+        >>> report = Report(rotor=rotor)
+        >>> report.mode_shpae(mode=0) # doctest: +ELLIPSIS
+        [...]
         """
         nodes_pos = self.rotor.nodes_pos
 
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index 7383f3f7..fc02a322 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -168,7 +168,7 @@ class PressureMatrix:
         radius_stator,
         viscosity,
         density,
-        beta=None,
+        beta=2.356,
         eccentricity=None,
         load=None,
     ):
@@ -202,6 +202,7 @@ class PressureMatrix:
         else:
             self.bearing_type = "medium_size"
         self.eccentricity = eccentricity
+        self.beta = beta
         self.eccentricity_ratio = None
         self.load = load
         if self.eccentricity is None:
@@ -211,12 +212,8 @@ class PressureMatrix:
         self.eccentricity_ratio = self.eccentricity / self.difference_between_radius
         if self.load is None:
             self.load = self.get_rotor_load()
-        if beta is None:
-            self.beta = self.calculate_attitude_angle()
-        else:
-            self.beta = beta
-        self.xi = self.eccentricity * np.cos(2*np.pi - self.beta)
-        self.yi = self.eccentricity * np.sin(2*np.pi - self.beta)
+        self.xi = self.eccentricity * np.cos(2*np.pi - beta)
+        self.yi = self.eccentricity * np.sin(2*np.pi - beta)
         self.re = np.zeros([self.nz, self.ntheta])
         self.ri = np.zeros([self.nz, self.ntheta])
         self.z = np.zeros([1, self.nz])
@@ -237,20 +234,6 @@ class PressureMatrix:
         self.analytical_pressure_matrix_available = False
         self.numerical_pressure_matrix_available = False
 
-    def calculate_attitude_angle(self):
-        """Calculates the attitude angle based on the eccentricity.
-        Returns
-        -------
-        float
-            Attitude angle
-        Examples
-        --------
-        >>> my_fluid_flow = pressure_matrix_example()
-        >>> my_fluid_flow.calculate_attitude_angle() # doctest: +ELLIPSIS
-        1.5...
-        """
-        return np.arctan(np.pi*(1 - self.eccentricity_ratio**2)/(4*self.eccentricity_ratio))
-
     def calculate_pressure_matrix_analytical(self, method=0, force_type=None):
         """This function calculates the pressure matrix analytically.
         Parameters
diff --git a/ross/results.py b/ross/results.py
index 56ebdff6..838a2a07 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -439,6 +439,143 @@ class FrequencyResponseResults:
         mag_plot.xaxis.axis_label_text_font_size = "14pt"
         mag_plot.yaxis.axis_label_text_font_size = "14pt"
 
+        minspeed = 320
+        maxspeed = 450
+
+        idx_max = argrelextrema(mag[inp, out, :], np.greater)
+        wn = frequency_range[idx_max[0]]
+
+        AF_table = []
+        SM_table = []
+        SM_ref_table = []
+
+        for i, peak in enumerate(mag[inp, out, :][idx_max[0]]):
+            peak_n = 0.707 * peak
+            peak_aux = np.linspace(peak_n, peak_n, len(frequency_range))
+
+            idx = np.argwhere(
+                np.diff(np.sign(peak_aux - mag[inp, out, :]))
+            ).flatten()
+            idx = np.sort(np.append(idx, idx_max[0][i]))
+
+            # if speed range is not long enough to catch the magnitudes
+            try:
+                idx_aux = [
+                    list(idx).index(idx_max[0][i]) - 1,
+                    list(idx).index(idx_max[0][i]) + 1,
+                ]
+                idx = idx[idx_aux]
+            except IndexError:
+                idx = [list(idx).index(idx_max[0][i]) - 1,
+                       len(frequency_range) - 1
+                ]
+
+            # Amplification Factor (AF)
+            AF = wn[i] / (
+                frequency_range[idx[1]] - frequency_range[idx[0]]
+            )
+
+            # Separation Margin (SM)
+            if AF > 2.5 and wn[i] < minspeed:
+                SM = min([16, 17 * (1 - 1 / (AF - 1.5))]) / 100
+                SMspeed = wn[i] * (1 + SM)
+                SM_ref = (minspeed - wn[i]) / wn[i]
+                source = ColumnDataSource(
+                    dict(
+                        top=[max(mag[inp, out, :][idx_max[0]])],
+                        bottom=[0],
+                        left=[wn[i]],
+                        right=[SMspeed],
+                        tag1=[wn[i]],
+                        tag2=[SMspeed],
+                    )
+                )
+
+            elif AF > 2.5 and wn[i] > maxspeed:
+                SM = min([26, 10 + 17 * (1 - 1 / (AF - 1.5))]) /100
+                SMspeed = wn[i] * (1 - SM)
+                SM_ref = (wn[i] - maxspeed) / maxspeed
+                source = ColumnDataSource(
+                    dict(
+                        top=[max(mag[inp, out, :][idx_max[0]])],
+                        bottom=[0],
+                        left=[SMspeed],
+                        right=[wn[i]],
+                        tag1=[wn[i]],
+                        tag2=[SMspeed],
+                    )
+                )
+
+            else:
+                SM = None
+                SM_ref = None
+                SMspeed = None
+
+            AF_table.append(AF)
+            SM_table.append(SM)
+            SM_ref_table.append(SM_ref)
+
+            mag_plot.quad(
+                top="top",
+                bottom="bottom",
+                left="left",
+                right="right",
+                source=source,
+                line_color=bokeh_colors[8],
+                line_width=0.8,
+                fill_alpha=0.2,
+                fill_color=bokeh_colors[8],
+                legend="Separation Margin",
+                name="SM2",
+            )
+            hover = HoverTool(names=["SM2"])
+            hover.tooltips = [
+                ("Critical Speed :", "@tag1"),
+                ("Speed at 0.707 x peak amplitude :", "@tag2"),
+            ]
+
+        table_source = ColumnDataSource(
+                dict(
+                    Wd=frequency_range[idx_max[0]],
+                    SM_table=SM_table,
+                    AF_table=AF_table,
+                    SM_ref_table=SM_ref_table,
+                )
+        )
+        form1 = NumberFormatter(format='0.00')
+        form2 = NumberFormatter(format='0.00%')
+        columns = [
+                TableColumn(
+                    field="Wd", title="Critical Speed", formatter=form1
+                ),
+                TableColumn(
+                    field="SM_table", title="Required Separation Margin", formatter=form2
+                ),
+                TableColumn(
+                    field="SM_ref_table", title="Available Separation Margin", formatter=form2
+                ),
+                TableColumn(
+                    field="AF_table", title="Amplification Factor", formatter=form1
+                )
+            ]
+        data_table = DataTable(
+                source=table_source, columns=columns, width=700, height=450
+        )
+        table = widgetbox(data_table)
+
+        mag_plot.quad(
+            top=max(mag[inp, out, :][idx_max[0]]),
+            bottom=0,
+            left=minspeed,
+            right=maxspeed,
+            line_color="green",
+            line_width=0.8,
+            fill_alpha=0.2,
+            fill_color="green",
+            legend="Operation Speed Range",
+        )
+        mag_plot.add_tools(hover)
+
         source = ColumnDataSource(dict(x=frequency_range, y=mag[inp, out, :]))
         mag_plot.line(
             x="x",
@@ -449,7 +586,7 @@ class FrequencyResponseResults:
             line_width=3,
         )
 
-        return mag_plot
+        return mag_plot, table
 
     def plot_phase_matplotlib(self, inp, out, ax=None, **kwargs):
         """Plot frequency response.
@@ -621,11 +758,11 @@ class FrequencyResponseResults:
         --------
         """
         # bokeh plot axes
-        bk_ax0 = self.plot_magnitude_bokeh(inp, out, ax=ax0)
+        bk_ax0, table = self.plot_magnitude_bokeh(inp, out, ax=ax0)
         bk_ax1 = self.plot_phase_bokeh(inp, out, ax=ax1)
 
         # show the bokeh plot results
-        grid_plots = gridplot([[bk_ax0], [bk_ax1]])
+        grid_plots = gridplot([[bk_ax0, table], [bk_ax1]])
         show(grid_plots)
 
         return bk_ax0, bk_ax1
@@ -972,39 +1109,6 @@ class ForcedResponseResults:
 
 
 class ModeShapeResults:
-    """Evaluates the mode shapes for the rotor.
-
-    This analysis presents the vibration mode for each critical speed.
-
-    Parameters
-    ----------
-    modes : list
-        list of eigenvectors
-    ndof : int
-        Number of degrees of freedom of the system
-    nodes : list
-        list of node numbers
-    nodes_pos : list
-        list of nodes positions
-    elements_length : list
-        list with length of each shaft element
-    w : float, list
-        rotor speed
-    wd : list
-        list with damped natural frequency
-    log_dec : list
-        list with logarithmic decrements
-    kappa_modes : list
-        list with values of kappa
-
-    Returns
-    -------
-        A graphic object
-
-    Example
-    -------
-    """
-
     def __init__(
         self,
         modes,
@@ -1028,37 +1132,6 @@ class ModeShapeResults:
         self.kappa_modes = kappa_modes
 
     def calc_mode_shape(self, mode=None, evec=None):
-        """
-        Method that calculate the arrays describing the mode shapes.
-
-        Parameters
-        ----------
-        mode : int
-            The n'th vibration mode
-            Default is None
-        evec : array
-            Array containing the system eigenvectors
-
-        Returns
-        -------
-        xn : array
-            absolut nodal displacement - X direction
-        yn : array
-            absolut nodal displacement - Y direction
-        zn : array
-            absolut nodal displacement - Z direction
-        x_circles : array
-            orbit description - X direction
-        y_circles : array
-            orbit description - Y direction
-        z_circles_pos : array
-            axial location of each orbit
-        nn : int
-            number of points to plot lines
-
-        Example
-        -------
-        """
         evec0 = self.modes[:, mode]
         nodes = self.nodes
         nodes_pos = self.nodes_pos
@@ -1128,31 +1201,6 @@ class ModeShapeResults:
         return xn, yn, zn, x_circles, y_circles, z_circles_pos, nn
 
     def plot(self, mode=None, evec=None, fig=None, ax=None):
-        """
-        Method that calculate the arrays describing the mode shapes.
-
-        Parameters
-        ----------
-        mode : int
-            The n'th vibration mode
-            Default is None
-        evec : array
-            Array containing the system eigenvectors
-        fig : matplotlib figure
-            The figure object with the plot.
-        ax : matplotlib axes
-            The axes object with the plot.
-
-        Returns
-        -------
-        fig : matplotlib figure
-            Returns the figure object with the plot.
-        ax : matplotlib axes
-            Returns the axes object with the plot.
-
-        Example
-        -------
-        """
         if ax is None:
             fig = plt.figure()
             ax = fig.gca(projection="3d")
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index d512dd90..066a72f1 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -227,7 +227,12 @@ class Rotor(object):
 
         df_shaft = pd.DataFrame([el.summary() for el in self.shaft_elements])
         df_disks = pd.DataFrame([el.summary() for el in self.disk_elements])
-        df_bearings = pd.DataFrame([el.summary() for el in self.bearing_seal_elements])
+        df_bearings = pd.DataFrame(
+                [el.summary() for el in self.bearing_seal_elements if (
+                        el.__class__.__name__ == 'BearingElement')])
+        df_seals = pd.DataFrame(
+                [el.summary() for el in self.bearing_seal_elements if (
+                        el.__class__.__name__ == 'SealElement')])
 
         nodes_pos_l = np.zeros(len(df_shaft.n_l))
         nodes_pos_r = np.zeros(len(df_shaft.n_l))
@@ -245,15 +250,15 @@ class Rotor(object):
 
         df_shaft["nodes_pos_l"] = nodes_pos_l
         df_shaft["nodes_pos_r"] = nodes_pos_r
-        # bearings
 
-        df = pd.concat([df_shaft, df_disks, df_bearings], sort=True)
+        df = pd.concat([df_shaft, df_disks, df_bearings, df_seals], sort=True)
         df = df.sort_values(by="n_l")
         df = df.reset_index(drop=True)
 
         self.df_disks = df_disks
         self.df_bearings = df_bearings
         self.df_shaft = df_shaft
+        self.df_seals = df_seals
 
         # check consistence for disks and bearings location
         if df.n_l.max() > df_shaft.n_r.max():
@@ -492,11 +497,11 @@ class Rotor(object):
                [ 0., -6.,  1.,  0.],
                [ 6.,  0.,  0.,  1.]])
         """
+        K0 = np.zeros((self.ndof, self.ndof))
+
         if frequency is None:
             frequency = self.w
 
-        K0 = np.zeros((self.ndof, self.ndof))
-
         for elm in self.elements:
             dofs = elm.dof_global_index()
             n0 = dofs[0]
@@ -529,11 +534,11 @@ class Rotor(object):
                [0., 0., 0., 0.],
                [0., 0., 0., 0.]])
         """
+        C0 = np.zeros((self.ndof, self.ndof))
+
         if frequency is None:
             frequency = self.w
 
-        C0 = np.zeros((self.ndof, self.ndof))
-
         for elm in self.elements:
             dofs = elm.dof_global_index()
             n0 = dofs[0]
@@ -1871,8 +1876,25 @@ class Rotor(object):
             for node_y in range(int(len(self.M()) / 4)):
                 grav[4 * node_y + 1] = -9.8065
 
+            aux_brg = []
+            for n in self.df_bearings["n"]:
+                aux_brg.append(BearingElement(n=n, kxx=1e14, cxx=0))
+
+            aux_rotor = Rotor(self.shaft_elements, self.disk_elements, aux_brg)
+            aux_K = np.zeros((aux_rotor.ndof, aux_rotor.ndof))
+
+            for elm in aux_rotor.elements:
+                if elm.__class__.__name__ != "SealElement":
+                    dofs = elm.dof_global_index()
+                    n0 = dofs[0]
+                    n1 = dofs[-1] + 1  # +1 to include this dof in the slice
+                    try:
+                        aux_K[n0:n1, n0:n1] += elm.K(0)
+                    except TypeError:
+                        aux_K[n0:n1, n0:n1] += elm.K()
+
             # calculates x, for [K]*(x) = [M]*(g)
-            disp = (la.solve(self.K(0), self.M() @ grav)).flatten()
+            disp = (la.solve(aux_K, self.M() @ grav)).flatten()
 
             # calculates displacement values in gravity's direction
             # dof = degree of freedom
@@ -1890,10 +1912,10 @@ class Rotor(object):
             BrgForceToReturn = []
             for i, node in enumerate(self.df_bearings["n"]):
                 BrgForce[node] = (
-                    -disp_y[node] * self.df_bearings.loc[i, "kyy"].coefficient[0]
+                    -disp_y[node] * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0]
                 )
                 BrgForceToReturn.append(
-                    np.around(-disp_y[node] * self.df_bearings.loc[i, "kyy"].coefficient[0], decimals=1)
+                    np.around(-disp_y[node] * aux_rotor.df_bearings.loc[i, "kyy"].coefficient[0], decimals=1)
                 )
 
             # Disk Forces
@@ -1973,7 +1995,7 @@ class Rotor(object):
                 self.Bm,
                 self.df_shaft,
                 self.df_disks,
-                self.df_bearings,
+                aux_rotor.df_bearings,
                 self.nodes,
                 self.nodes_pos,
                 Vx_axis,


Test patch:
diff --git a/ross/tests/test_rotor_assembly.py b/ross/tests/test_rotor_assembly.py
index 059cc3ae..fee2a570 100644
--- a/ross/tests/test_rotor_assembly.py
+++ b/ross/tests/test_rotor_assembly.py
@@ -711,13 +711,13 @@ def test_static_analysis_rotor3(rotor3):
         rotor3.disp_y,
         (
             [
-                -0.00061784,
-                -0.00108199,
-                -0.00143205,
-                -0.00157464,
-                -0.00147798,
-                -0.00115111,
-                -0.00069521,
+                -4.94274533e-12,
+                -4.51249085e-04,
+                -7.88420867e-04,
+                -9.18114192e-04,
+                -8.08560219e-04,
+                -4.68788888e-04,
+                -5.56171636e-12
             ]
         ),
         decimal=6,
@@ -779,17 +779,17 @@ def test_static_analysis_rotor5(rotor5):
         rotor5.disp_y,
         (
             [
-                0.00026382,
-                -0.00015021,
-                -0.00056947,
-                -0.00098566,
-                -0.00130592,
-                -0.00143686,
-                -0.00134669,
-                -0.00104446,
-                -0.00063136,
-                -0.00021282,
-                0.0002005,
+                8.12651626e-04,
+                4.08939282e-04,
+                -5.69465378e-12,
+                -4.05876595e-04,
+                -7.15824882e-04,
+                -8.36443708e-04,
+                -7.35964234e-04,
+                -4.23416398e-04,
+                -6.31362481e-12,
+                4.28859620e-04,
+                8.52492302e-04
             ]
         ),
         decimal=6,
@@ -868,17 +868,17 @@ def test_static_analysis_rotor6(rotor6):
         rotor6.disp_y,
         (
             [
-                -2.58459386e-06,
-                -8.83514255e-05,
-                -1.79345202e-04,
-                -2.82280131e-04,
-                -3.83451012e-04,
-                -4.71329601e-04,
-                -5.55753491e-04,
-                -7.08192594e-04,
-                -1.02148266e-03,
-                -1.55824814e-03,
-                -2.22220183e-03,
+                -1.03951876e-04,
+                -4.93624668e-05,
+                -1.79345202e-12,
+                3.74213098e-05,
+                7.66066703e-05,
+                1.29084322e-04,
+                1.85016673e-04,
+                1.72933811e-04,
+                -1.02148266e-11,
+                -3.96409257e-04,
+                -9.20006704e-04
             ]
         ),
         decimal=6,