Repository: ross-rotordynamics/ross
Pull request #697 (2020-11-03 19:46:42)

Building random bearings with freq. dependent coeficients
* These modifications allow the user to generate random coefficients with frequency dependency using `fluidflow`.
* The argument `omega` cannot be considered random. The user inputs a list of frequencies and the code calculate the coefficients' statistics for each frequency value.

codecov[bot]: # [Codecov](https://codecov.io/gh/ross-rotordynamics/ross/pull/697?src=pr&el=h1) Report
> Merging [#697](https://codecov.io/gh/ross-rotordynamics/ross/pull/697?src=pr&el=desc) into [master](https://codecov.io/gh/ross-rotordynamics/ross/commit/4bd744c3d97449fdccb12c03a7da340bb4d33c24?el=desc) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/ross-rotordynamics/ross/pull/697/graphs/tree.svg?width=650&height=150&src=pr&token=SZjY7TrJnS)](https://codecov.io/gh/ross-rotordynamics/ross/pull/697?src=pr&el=tree)

```diff
@@           Coverage Diff           @@
##           master     #697   +/-   ##
=======================================
  Coverage   90.76%   90.76%           
=======================================
  Files          29       29           
  Lines        5218     5221    +3     
=======================================
+ Hits         4736     4739    +3     
  Misses        482      482           
```


| [Impacted Files](https://codecov.io/gh/ross-rotordynamics/ross/pull/697?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [ross/stochastic/st\_bearing\_seal\_element.py](https://codecov.io/gh/ross-rotordynamics/ross/pull/697/diff?src=pr&el=tree#diff-cm9zcy9zdG9jaGFzdGljL3N0X2JlYXJpbmdfc2VhbF9lbGVtZW50LnB5) | `98.86% <100.00%> (+0.04%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/ross-rotordynamics/ross/pull/697?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/ross-rotordynamics/ross/pull/697?src=pr&el=footer). Last update [4bd744c...a10596f](https://codecov.io/gh/ross-rotordynamics/ross/pull/697?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index d059eb0c..8b1fb6d4 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
 -   repo: https://github.com/ambv/black
-    rev: stable
+    rev: 20.8b1
     hooks:
     - id: black
       language_version: python3.7
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index 41e86265..ffea2514 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -1,7 +1,7 @@
 """Bearing Element module.
 
 This module defines the BearingElement classes which will be used to represent the rotor
-bearings and seals. There're 6 different classes to represent bearings options,
+bearings and seals. There are 7 different classes to represent bearings options,
 and 2 element options with 8 or 12 degrees of freedom.
 """
 # fmt: off
@@ -27,6 +27,7 @@ __all__ = [
     "SealElement",
     "BallBearingElement",
     "RollerBearingElement",
+    "BearingFluidFlow",
     "BearingElement6DoF",
     "MagneticBearingElement",
 ]
@@ -912,9 +913,103 @@ class BearingElement(Element):
             color=color,
         )
 
-    @classmethod
-    def from_fluid_flow(
-        cls,
+
+class BearingFluidFlow(BearingElement):
+    """Instantiate a bearing using inputs from its fluid flow.
+
+    This method always creates elements with frequency-dependent coefficients.
+    It calculates a set of coefficients for each frequency value appendend to
+    "omega".
+
+    Parameters
+    ----------
+    n : int
+        The node in which the bearing will be located in the rotor.
+
+    Grid related
+    ^^^^^^^^^^^^
+    Describes the discretization of the problem
+    nz: int
+        Number of points along the Z direction (direction of flow).
+    ntheta: int
+        Number of points along the direction theta. NOTE: ntheta must be odd.
+    nradius: int
+        Number of points along the direction r.
+    length: float
+        Length in the Z direction (m).
+
+    Operation conditions
+    ^^^^^^^^^^^^^^^^^^^^
+    Describes the operation conditions.
+    omega: list
+        List of frequencies (rad/s) used to calculate the coefficients.
+        If the length is greater than 1, an array of coefficients is returned.
+    p_in: float
+        Input Pressure (Pa).
+    p_out: float
+        Output Pressure (Pa).
+    load: float
+        Load applied to the rotor (N).
+
+    Geometric data of the problem
+    ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+    Describes the geometric data of the problem.
+    radius_rotor: float
+        Rotor radius (m).
+    radius_stator: float
+        Stator Radius (m).
+    eccentricity: float
+        Eccentricity (m) is the euclidean distance between rotor and stator centers.
+        The center of the stator is in position (0,0).
+
+    Fluid characteristics
+    ^^^^^^^^^^^^^^^^^^^^^
+    Describes the fluid characteristics.
+    visc: float
+        Viscosity (Pa.s).
+    rho: float
+        Fluid density(Kg/m^3).
+
+    Others
+    ^^^^^^
+    tag : str, optional
+        A tag to name the element
+        Default is None.
+    n_link : int, optional
+        Node to which the bearing will connect. If None the bearing is
+        connected to ground.
+        Default is None.
+    scale_factor : float, optional
+        The scale factor is used to scale the bearing drawing.
+        Default is 1.
+
+    Returns
+    -------
+    bearing: rs.BearingElement
+        A bearing object.
+
+    Examples
+    --------
+    >>> nz = 30
+    >>> ntheta = 20
+    >>> length = 0.03
+    >>> omega = [157.1]
+    >>> p_in = 0.
+    >>> p_out = 0.
+    >>> radius_rotor = 0.0499
+    >>> radius_stator = 0.05
+    >>> load = 525
+    >>> visc = 0.1
+    >>> rho = 860.
+    >>> BearingFluidFlow(0, nz, ntheta, length, omega, p_in,
+    ...                  p_out, radius_rotor, radius_stator,
+    ...                  visc, rho, load=load) # doctest: +ELLIPSIS
+    BearingFluidFlow(n=0, n_link=None,
+     kxx=array([145...
+    """
+
+    def __init__(
+        self,
         n,
         nz,
         ntheta,
@@ -931,99 +1026,8 @@ class BearingElement(Element):
         tag=None,
         n_link=None,
         scale_factor=1.0,
+        color="#355d7a",
     ):
-        """Instantiate a bearing using inputs from its fluid flow.
-
-        This method always creates elements with frequency-dependent coefficients.
-        It calculates a set of coefficients for each frequency value appendend to
-        "omega".
-
-        Parameters
-        ----------
-        n : int
-            The node in which the bearing will be located in the rotor.
-
-        Grid related
-        ^^^^^^^^^^^^
-        Describes the discretization of the problem
-        nz: int
-            Number of points along the Z direction (direction of flow).
-        ntheta: int
-            Number of points along the direction theta. NOTE: ntheta must be odd.
-        nradius: int
-            Number of points along the direction r.
-        length: float
-            Length in the Z direction (m).
-
-        Operation conditions
-        ^^^^^^^^^^^^^^^^^^^^
-        Describes the operation conditions.
-        omega: list
-            List of frequencies (rad/s) used to calculate the coefficients.
-            If the length is greater than 1, an array of coefficients is returned.
-        p_in: float
-            Input Pressure (Pa).
-        p_out: float
-            Output Pressure (Pa).
-        load: float
-            Load applied to the rotor (N).
-
-        Geometric data of the problem
-        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-        Describes the geometric data of the problem.
-        radius_rotor: float
-            Rotor radius (m).
-        radius_stator: float
-            Stator Radius (m).
-        eccentricity: float
-            Eccentricity (m) is the euclidean distance between rotor and stator centers.
-            The center of the stator is in position (0,0).
-
-        Fluid characteristics
-        ^^^^^^^^^^^^^^^^^^^^^
-        Describes the fluid characteristics.
-        visc: float
-            Viscosity (Pa.s).
-        rho: float
-            Fluid density(Kg/m^3).
-
-        Others
-        ^^^^^^
-        tag : str, optional
-            A tag to name the element
-            Default is None.
-        n_link : int, optional
-            Node to which the bearing will connect. If None the bearing is
-            connected to ground.
-            Default is None.
-        scale_factor : float, optional
-            The scale factor is used to scale the bearing drawing.
-            Default is 1.
-
-        Returns
-        -------
-        bearing: rs.BearingElement
-            A bearing object.
-
-        Examples
-        --------
-        >>> nz = 30
-        >>> ntheta = 20
-        >>> length = 0.03
-        >>> omega = [157.1]
-        >>> p_in = 0.
-        >>> p_out = 0.
-        >>> radius_rotor = 0.0499
-        >>> radius_stator = 0.05
-        >>> load = 525
-        >>> visc = 0.1
-        >>> rho = 860.
-        >>> BearingElement.from_fluid_flow(0, nz, ntheta, length, omega, p_in,
-        ...                                p_out, radius_rotor, radius_stator,
-        ...                                visc, rho, load=load) # doctest: +ELLIPSIS
-        BearingElement(n=0, n_link=None,
-         kxx=array([...
-        """
         K = np.zeros((4, len(omega)))
         C = np.zeros((4, len(omega)))
 
@@ -1044,7 +1048,7 @@ class BearingElement(Element):
             )
             K[:, i], C[:, i] = calculate_stiffness_and_damping_coefficients(fluid_flow)
 
-        return cls(
+        super().__init__(
             n,
             kxx=K[0],
             kxy=K[1],
@@ -1055,6 +1059,10 @@ class BearingElement(Element):
             cyx=C[2],
             cyy=C[3],
             frequency=omega,
+            tag=tag,
+            n_link=n_link,
+            scale_factor=scale_factor,
+            color=color,
         )
 
 
diff --git a/ross/defects/abs_defect.py b/ross/defects/abs_defect.py
index 0dba4d29..1d25c164 100644
--- a/ross/defects/abs_defect.py
+++ b/ross/defects/abs_defect.py
@@ -7,9 +7,7 @@ import scipy as sp
 from ross.results import TimeResponseResults
 from ross.units import Q_
 
-__all__ = [
-    "Defect",
-]
+__all__ = ["Defect"]
 
 
 class Defect(ABC):
@@ -31,11 +29,8 @@ class Defect(ABC):
         )
         return results
 
-    def plot_dfft(
-        self, probe, probe_units="rad", range_freq=None, fig=None, **kwargs,
-    ):
-        """
-        """
+    def plot_dfft(self, probe, probe_units="rad", range_freq=None, fig=None, **kwargs):
+        """"""
         if fig is None:
             fig = go.Figure()
 
diff --git a/ross/defects/integrate_solver.py b/ross/defects/integrate_solver.py
index b9bdf230..1020de1a 100644
--- a/ross/defects/integrate_solver.py
+++ b/ross/defects/integrate_solver.py
@@ -12,7 +12,7 @@ import numpy as np
 class Integrator:
 
     """A series of Runge-Kutta time integration algorithms.
-    
+
     Calculates the time response for the rotors input to the routine.
 
     Parameters
diff --git a/ross/defects/misalignment.py b/ross/defects/misalignment.py
index 9a992762..1ec7610b 100644
--- a/ross/defects/misalignment.py
+++ b/ross/defects/misalignment.py
@@ -19,16 +19,13 @@ from ross.units import Q_
 from .abs_defect import Defect
 from .integrate_solver import Integrator
 
-__all__ = [
-    "MisalignmentFlex",
-    "MisalignmentRigid",
-]
+__all__ = ["MisalignmentFlex", "MisalignmentRigid"]
 
 
 class MisalignmentFlex(Defect):
     """A flexible coupling with misalignment of some kind.
-    
-    Calculates the dynamic reaction force of hexangular flexible coupling 
+
+    Calculates the dynamic reaction force of hexangular flexible coupling
     induced by 6DOF's rotor parallel and angular misalignment.
 
     Parameters
@@ -66,15 +63,15 @@ class MisalignmentFlex(Defect):
     print_progress : bool
         Set it True, to print the time iterations and the total time spent.
         False by default.
-    
+
     Returns
     -------
     A force to be applied on the shaft.
 
     References
     ----------
-    .. [1] 'Xia, Y., Pang, J., Yang, L., Zhao, Q., & Yang, X. (2019). Study on vibration response 
-    and orbits of misaligned rigid rotors connected by hexangular flexible coupling. Applied 
+    .. [1] 'Xia, Y., Pang, J., Yang, L., Zhao, Q., & Yang, X. (2019). Study on vibration response
+    and orbits of misaligned rigid rotors connected by hexangular flexible coupling. Applied
     Acoustics, 155, 286-296 ..
 
     Examples
@@ -147,7 +144,7 @@ class MisalignmentFlex(Defect):
         ----------
         rotor : ross.Rotor Object
              6 DoF rotor model.
-                
+
         """
         self.rotor = rotor
         self.radius = rotor.elements[self.n1].odl / 2
@@ -199,7 +196,7 @@ class MisalignmentFlex(Defect):
         self.M = self.rotor.M()
         self.Kst = self.rotor.Kst()
 
-        _, ModMat = scipy.linalg.eigh(self.K, self.M, type=1, turbo=False,)
+        _, ModMat = scipy.linalg.eigh(self.K, self.M, type=1, turbo=False)
         ModMat = ModMat[:, :12]
         self.ModMat = ModMat
 
@@ -278,7 +275,7 @@ class MisalignmentFlex(Defect):
         self.response = self.ModMat.dot(self.displacement)
 
     def _equation_of_movement(self, T, Y, i):
-        """ Calculates the displacement and velocity using state-space representation in the modal domain.
+        """Calculates the displacement and velocity using state-space representation in the modal domain.
 
         Parameters
         ----------
@@ -318,14 +315,14 @@ class MisalignmentFlex(Defect):
 
     def _parallel(self, angular_position):
         """Reaction forces of parallel misalignment.
-        
+
         angular_position : float
                         Angular position of the shaft.
-        
+
         Returns
         -------
         F_mis_p : array
-               Excitation force caused by the parallel misalignment on the entire system. 
+               Excitation force caused by the parallel misalignment on the entire system.
         """
 
         F_mis_p = np.zeros((self.ndof, len(angular_position)))
@@ -434,7 +431,7 @@ class MisalignmentFlex(Defect):
 
         angular_position : float
                 Angular position of the shaft.
-        
+
         Returns
         -------
         F_mis_a : array
@@ -494,7 +491,7 @@ class MisalignmentFlex(Defect):
 
         angular_position : float
                 Angular position of the shaft.
-        
+
         Returns
         -------
         F_misalign : array
@@ -506,8 +503,8 @@ class MisalignmentFlex(Defect):
 
 class MisalignmentRigid(Defect):
     """A rigid coupling with parallel misalignment.
-    
-    Calculates the dynamic reaction force of hexangular rigid coupling 
+
+    Calculates the dynamic reaction force of hexangular rigid coupling
     induced by 6DOF's rotor parallel misalignment.
 
     Parameters
@@ -595,7 +592,7 @@ class MisalignmentRigid(Defect):
         ----------
         rotor : ross.Rotor Object
              6 DoF rotor model.
-                
+
         """
         self.rotor = rotor
         self.ndof = rotor.ndof
@@ -642,7 +639,7 @@ class MisalignmentRigid(Defect):
         self.M = self.rotor.M()
         self.Kst = self.rotor.Kst()
 
-        _, ModMat = scipy.linalg.eigh(self.K, self.M, type=1, turbo=False,)
+        _, ModMat = scipy.linalg.eigh(self.K, self.M, type=1, turbo=False)
         ModMat = ModMat[:, :12]
         self.ModMat = ModMat
 
@@ -742,7 +739,7 @@ class MisalignmentRigid(Defect):
         self.response = self.ModMat.dot(self.displacement)
 
     def _equation_of_movement(self, T, Y, i):
-        """ Calculates the displacement and velocity using state-space representation in the modal domain.
+        """Calculates the displacement and velocity using state-space representation in the modal domain.
 
         Parameters
         ----------
@@ -799,7 +796,7 @@ class MisalignmentRigid(Defect):
 
     def _parallel(self, positions, fir):
         """Reaction forces of parallel misalignment.
-        
+
         Returns
         -------
         Fmis : array
@@ -944,8 +941,8 @@ def base_rotor_example():
 def misalignment_flex_parallel_example():
     """Create an example of a flexible parallel misalignment defect.
 
-    This function returns an instance of a flexible parallel misalignment 
-    defect. The purpose is to make available a simple model so that a 
+    This function returns an instance of a flexible parallel misalignment
+    defect. The purpose is to make available a simple model so that a
     doctest can be written using it.
 
     Returns
@@ -988,8 +985,8 @@ def misalignment_flex_parallel_example():
 def misalignment_flex_angular_example():
     """Create an example of a flexible angular misalignment defect.
 
-    This function returns an instance of a flexible angular misalignment 
-    defect. The purpose is to make available a simple model so that a 
+    This function returns an instance of a flexible angular misalignment
+    defect. The purpose is to make available a simple model so that a
     doctest can be written using it.
 
     Returns
@@ -1032,8 +1029,8 @@ def misalignment_flex_angular_example():
 def misalignment_flex_combined_example():
     """Create an example of a flexible combined misalignment defect.
 
-    This function returns an instance of a flexible combined misalignment 
-    defect. The purpose is to make available a simple model so that a 
+    This function returns an instance of a flexible combined misalignment
+    defect. The purpose is to make available a simple model so that a
     doctest can be written using it.
 
     Returns
@@ -1076,8 +1073,8 @@ def misalignment_flex_combined_example():
 def misalignment_rigid_example():
     """Create an example of a rigid misalignment defect.
 
-    This function returns an instance of a rigid misalignment 
-    defect. The purpose is to make available a simple model so that a 
+    This function returns an instance of a rigid misalignment
+    defect. The purpose is to make available a simple model so that a
     doctest can be written using it.
 
     Returns
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index a3e42073..fb36652b 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -423,28 +423,18 @@ class FluidFlow:
                 for i in range(0, self.nz):
                     for j in range(0, self.ntheta):
                         self.p_mat_analytical[i, j] = (
-                            (
-                                6
-                                * self.viscosity
-                                * self.omega
-                                * (self.radius_rotor / self.radial_clearance) ** 2
-                                * self.eccentricity_ratio
-                                * np.sin(self.dtheta * j)
-                                * (
-                                    2
-                                    + self.eccentricity_ratio * np.cos(self.dtheta * j)
-                                )
-                            )
-                            / (
-                                (2 + self.eccentricity_ratio ** 2)
-                                * (
-                                    1
-                                    + self.eccentricity_ratio * np.cos(self.dtheta * j)
-                                )
-                                ** 2
-                            )
-                            + self.p_in
-                        )
+                            6
+                            * self.viscosity
+                            * self.omega
+                            * (self.radius_rotor / self.radial_clearance) ** 2
+                            * self.eccentricity_ratio
+                            * np.sin(self.dtheta * j)
+                            * (2 + self.eccentricity_ratio * np.cos(self.dtheta * j))
+                        ) / (
+                            (2 + self.eccentricity_ratio ** 2)
+                            * (1 + self.eccentricity_ratio * np.cos(self.dtheta * j))
+                            ** 2
+                        ) + self.p_in
                         if self.p_mat_analytical[i, j] < 0:
                             self.p_mat_analytical[i, j] = 0
         elif self.bearing_type == "medium_size" or self.shape_geometry != "cylindrical":
diff --git a/ross/fluid_flow/fluid_flow_coefficients.py b/ross/fluid_flow/fluid_flow_coefficients.py
index 7d11423b..82621679 100644
--- a/ross/fluid_flow/fluid_flow_coefficients.py
+++ b/ross/fluid_flow/fluid_flow_coefficients.py
@@ -423,7 +423,7 @@ def find_equilibrium_position(fluid_flow_object, print_equilibrium_position=Fals
         )
         bearing.geometry_description()
         bearing.calculate_pressure_matrix_numerical()
-        (_, _, fx, fy,) = calculate_oil_film_force(bearing, force_type="numerical")
+        (_, _, fx, fy) = calculate_oil_film_force(bearing, force_type="numerical")
         return np.array([fx, (fy - bearing.load)])
 
     if fluid_flow_object.load is None:
@@ -437,15 +437,9 @@ def find_equilibrium_position(fluid_flow_object, print_equilibrium_position=Fals
     move_rotor_center_abs(fluid_flow_object, x0[0], x0[1])
     fluid_flow_object.geometry_description()
     fluid_flow_object.calculate_pressure_matrix_numerical()
-    (_, _, fx, fy,) = calculate_oil_film_force(
-        fluid_flow_object, force_type="numerical"
-    )
+    (_, _, fx, fy) = calculate_oil_film_force(fluid_flow_object, force_type="numerical")
     result = least_squares(
-        residuals,
-        x0,
-        args=[fluid_flow_object],
-        jac="3-point",
-        bounds=([0, -1], [1, 0]),
+        residuals, x0, args=[fluid_flow_object], jac="3-point", bounds=([0, -1], [1, 0])
     )
     move_rotor_center_abs(
         fluid_flow_object,
diff --git a/ross/results.py b/ross/results.py
index 157cdee1..7928eec9 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -2140,7 +2140,7 @@ class ForcedResponseResults:
 
         fig.update_xaxes(title_text=f"Rotor Length ({rotor_length_units})")
         fig.update_yaxes(
-            title_text=f"Bending Moment ({moment_units})", title_font=dict(size=12),
+            title_text=f"Bending Moment ({moment_units})", title_font=dict(size=12)
         )
         fig.update_layout(**kwargs)
 
@@ -2254,10 +2254,10 @@ class ForcedResponseResults:
                 domain=dict(x=[0.47, 1]),
             ),
             title=dict(
-                text=f"Deflected Shape<br>Speed = {speed_str} {frequency_units}",
+                text=f"Deflected Shape<br>Speed = {speed_str} {frequency_units}"
             ),
             legend=dict(
-                orientation="h", xanchor="center", yanchor="bottom", x=0.5, y=-0.3,
+                orientation="h", xanchor="center", yanchor="bottom", x=0.5, y=-0.3
             ),
             width=800,
             height=600,
diff --git a/ross/rotor_assembly.py b/ross/rotor_assembly.py
index 88023764..f96b8f0f 100644
--- a/ross/rotor_assembly.py
+++ b/ross/rotor_assembly.py
@@ -2377,7 +2377,7 @@ class Rotor(object):
             Coupling type. The avaible types are: flex, by default; and rigid.
 
         **kwargs: dictionary
-        
+
             In the case of coupling = "flex", **kwargs receives:
                 dt : float
                     Time step.
diff --git a/ross/stochastic/st_bearing_seal_element.py b/ross/stochastic/st_bearing_seal_element.py
index 4f67c6b1..b72a1c7f 100644
--- a/ross/stochastic/st_bearing_seal_element.py
+++ b/ross/stochastic/st_bearing_seal_element.py
@@ -464,36 +464,37 @@ class ST_BearingElement:
         >>> import ross.stochastic as srs
         >>> nz = 8
         >>> ntheta = 64
-        >>> length = 0.01
-        >>> omega = 157.1
+        >>> length = 1.5 * 0.0254
+        >>> omega = [157.1, 300]
         >>> p_in = 0.
         >>> p_out = 0.
-        >>> radius_rotor = 0.08
-        >>> radius_stator = 0.1
-        >>> viscosity = np.random.uniform(0.01, 0.02, 5)
-        >>> density = 860.
-        >>> eccentricity = 0.001
-        >>> attitude_angle = np.pi
-        >>> elms = ST_BearingElement.from_fluid_flow(
-        ...     0, nz, ntheta, length,
-        ...     omega, p_in, p_out, radius_rotor,
-        ...     radius_stator, viscosity, density, attitude_angle,
-        ...     eccentricity, is_random=["viscosity"]
+        >>> radius_rotor = 3 * 0.0254
+        >>> radius_stator = 3.003 * 0.0254
+        >>> viscosity = np.random.uniform(2.4e-03, 2.8e-03, 5)
+        >>> density = 860
+        >>> load = 1244.1
+        >>> elms = srs.ST_BearingElement.from_fluid_flow(
+        ...     0, nz=nz, ntheta=ntheta, length=length,
+        ...     omega=omega, p_in=p_in, p_out=p_out, radius_rotor=radius_rotor,
+        ...     radius_stator=radius_stator, viscosity=viscosity, density=density,
+        ...     load=load, is_random=["viscosity"]
         ... )
         >>> len(list(iter(elms)))
         5
         """
         attribute_dict = locals()
+        attribute_dict.pop("cls")
+        attribute_dict.pop("omega")
         size = len(attribute_dict[is_random[0]])
         args_dict = {
-            "kxx": [],
-            "kxy": [],
-            "kyx": [],
-            "kyy": [],
-            "cxx": [],
-            "cxy": [],
-            "cyx": [],
-            "cyy": [],
+            "kxx": np.zeros((len(omega), size)),
+            "kxy": np.zeros((len(omega), size)),
+            "kyx": np.zeros((len(omega), size)),
+            "kyy": np.zeros((len(omega), size)),
+            "cxx": np.zeros((len(omega), size)),
+            "cxy": np.zeros((len(omega), size)),
+            "cyx": np.zeros((len(omega), size)),
+            "cyy": np.zeros((len(omega), size)),
         }
 
         for k, v in attribute_dict.items():
@@ -502,33 +503,34 @@ class ST_BearingElement:
             else:
                 attribute_dict[k] = np.asarray(v)
 
-        for i in range(size):
-            fluid_flow = flow.FluidFlow(
-                attribute_dict["nz"][i],
-                attribute_dict["ntheta"][i],
-                attribute_dict["length"][i],
-                attribute_dict["omega"][i],
-                attribute_dict["p_in"][i],
-                attribute_dict["p_out"][i],
-                attribute_dict["radius_rotor"][i],
-                attribute_dict["radius_stator"][i],
-                attribute_dict["viscosity"][i],
-                attribute_dict["density"][i],
-                attribute_dict["attitude_angle"][i],
-                attribute_dict["eccentricity"][i],
-                attribute_dict["load"][i],
-                attribute_dict["omegap"][i],
-                attribute_dict["immediately_calculate_pressure_matrix_numerically"][i],
-            )
-            k, c = calculate_stiffness_and_damping_coefficients(fluid_flow)
-            args_dict["kxx"].append(k[0])
-            args_dict["kxy"].append(k[1])
-            args_dict["kyx"].append(k[2])
-            args_dict["kyy"].append(k[3])
-            args_dict["cxx"].append(c[0])
-            args_dict["cxy"].append(c[1])
-            args_dict["cyx"].append(c[2])
-            args_dict["cyy"].append(c[3])
+        for j, frequency in enumerate(omega):
+            for i in range(size):
+                fluid_flow = flow.FluidFlow(
+                    attribute_dict["nz"][i],
+                    attribute_dict["ntheta"][i],
+                    attribute_dict["length"][i],
+                    frequency,
+                    attribute_dict["p_in"][i],
+                    attribute_dict["p_out"][i],
+                    attribute_dict["radius_rotor"][i],
+                    attribute_dict["radius_stator"][i],
+                    attribute_dict["viscosity"][i],
+                    attribute_dict["density"][i],
+                    attribute_dict["attitude_angle"][i],
+                    attribute_dict["eccentricity"][i],
+                    attribute_dict["load"][i],
+                    attribute_dict["omegap"][i],
+                )
+                k, c = calculate_stiffness_and_damping_coefficients(fluid_flow)
+
+                args_dict["kxx"][j, i] = k[0]
+                args_dict["kxy"][j, i] = k[1]
+                args_dict["kyx"][j, i] = k[2]
+                args_dict["kyy"][j, i] = k[3]
+                args_dict["cxx"][j, i] = c[0]
+                args_dict["cxy"][j, i] = c[1]
+                args_dict["cyx"][j, i] = c[2]
+                args_dict["cyy"][j, i] = c[3]
 
         return cls(
             n,
@@ -540,7 +542,7 @@ class ST_BearingElement:
             cyy=args_dict["cyy"],
             cxy=args_dict["cxy"],
             cyx=args_dict["cyx"],
-            frequency=[fluid_flow.omega],
+            frequency=omega,
             tag=tag,
             n_link=n_link,
             scale_factor=scale_factor,
diff --git a/ross/stochastic/st_disk_element.py b/ross/stochastic/st_disk_element.py
index fa057d80..822e1c78 100644
--- a/ross/stochastic/st_disk_element.py
+++ b/ross/stochastic/st_disk_element.py
@@ -55,9 +55,7 @@ class ST_DiskElement:
     5
     """
 
-    def __init__(
-        self, n, m, Id, Ip, tag=None, color="Firebrick", is_random=None,
-    ):
+    def __init__(self, n, m, Id, Ip, tag=None, color="Firebrick", is_random=None):
         attribute_dict = dict(n=n, m=m, Id=Id, Ip=Ip, tag=tag, color=color)
 
         self.is_random = is_random
@@ -218,7 +216,7 @@ class ST_DiskElement:
         >>> # fig.show()
         """
         label = dict(
-            m="Mass", Id="Diametral moment of inertia", Ip="Polar moment of inertia",
+            m="Mass", Id="Diametral moment of inertia", Ip="Polar moment of inertia"
         )
         if var_list is None:
             var_list = self.is_random
@@ -234,9 +232,7 @@ class ST_DiskElement:
         )
 
     @classmethod
-    def from_geometry(
-        cls, n, material, width, i_d, o_d, tag=None, is_random=None,
-    ):
+    def from_geometry(cls, n, material, width, i_d, o_d, tag=None, is_random=None):
         """Random disk element.
 
         Creates an object containing a list with random instances of
@@ -297,7 +293,7 @@ class ST_DiskElement:
             o_d = np.array(o_d)
 
         attribute_dict = dict(
-            n=n, material=material, width=width, i_d=i_d, o_d=o_d, tag=tag,
+            n=n, material=material, width=width, i_d=i_d, o_d=o_d, tag=tag
         )
         size = len(attribute_dict[is_random[0]])
 
@@ -340,6 +336,6 @@ def st_disk_example():
     2
     """
     elm = ST_DiskElement(
-        n=1, m=[30, 40], Id=[0.2, 0.3], Ip=[0.5, 0.7], is_random=["m", "Id", "Ip"],
+        n=1, m=[30, 40], Id=[0.2, 0.3], Ip=[0.5, 0.7], is_random=["m", "Id", "Ip"]
     )
     return elm
diff --git a/ross/stochastic/st_materials.py b/ross/stochastic/st_materials.py
index d0f95f2a..8adde6df 100644
--- a/ross/stochastic/st_materials.py
+++ b/ross/stochastic/st_materials.py
@@ -91,7 +91,7 @@ class ST_Material:
             Poisson = np.asarray(Poisson)
 
         attribute_dict = dict(
-            name=name, rho=rho, E=E, G_s=G_s, Poisson=Poisson, color=color,
+            name=name, rho=rho, E=E, G_s=G_s, Poisson=Poisson, color=color
         )
         self.is_random = is_random
         self.attribute_dict = attribute_dict
diff --git a/ross/stochastic/st_point_mass.py b/ross/stochastic/st_point_mass.py
index 86205b0a..1613280a 100644
--- a/ross/stochastic/st_point_mass.py
+++ b/ross/stochastic/st_point_mass.py
@@ -57,9 +57,9 @@ class ST_PointMass:
     """
 
     def __init__(
-        self, n, m=None, mx=None, my=None, tag=None, color="DarkSalmon", is_random=None,
+        self, n, m=None, mx=None, my=None, tag=None, color="DarkSalmon", is_random=None
     ):
-        attribute_dict = dict(n=n, m=m, mx=mx, my=my, tag=tag, color=color,)
+        attribute_dict = dict(n=n, m=m, mx=mx, my=my, tag=tag, color=color)
 
         self.is_random = is_random
         self.attribute_dict = attribute_dict
@@ -217,7 +217,7 @@ class ST_PointMass:
         >>> # fig.show()
         """
         label = dict(
-            mx="Mass on the X direction", my="Mass on the Y direction", m="Mass",
+            mx="Mass on the X direction", my="Mass on the Y direction", m="Mass"
         )
         if var_list is None:
             var_list = self.is_random
diff --git a/ross/stochastic/st_results.py b/ross/stochastic/st_results.py
index 0af9c744..d0752c32 100644
--- a/ross/stochastic/st_results.py
+++ b/ross/stochastic/st_results.py
@@ -375,7 +375,7 @@ class ST_FrequencyResponseResults:
         self.phase = phase
 
     def plot_magnitude(
-        self, percentile=[], conf_interval=[], units="mic-pk-pk", **kwargs,
+        self, percentile=[], conf_interval=[], units="mic-pk-pk", **kwargs
     ):
         """Plot amplitude vs frequency.
 
@@ -611,7 +611,7 @@ class ST_FrequencyResponseResults:
         return fig
 
     def plot_polar_bode(
-        self, percentile=[], conf_interval=[], units="mic-pk-pk", **kwargs,
+        self, percentile=[], conf_interval=[], units="mic-pk-pk", **kwargs
     ):
         """Plot polar forced response using Plotly.
 
@@ -718,7 +718,7 @@ class ST_FrequencyResponseResults:
                     linecolor="black",
                     linewidth=2.5,
                 ),
-            ),
+            )
         )
 
         return fig
@@ -1167,7 +1167,7 @@ class ST_TimeResponseResults:
                 xaxis=dict(title=dict(text="<b>Rotor Length</b>"), showspikes=False),
                 yaxis=dict(title=dict(text="<b>Amplitude - X</b>"), showspikes=False),
                 zaxis=dict(title=dict(text="<b>Amplitude - Y</b>"), showspikes=False),
-            ),
+            )
         )
         return fig
 
@@ -1644,7 +1644,7 @@ class ST_ForcedResponseResults:
             polar=dict(
                 radialaxis=fig2.layout.polar.radialaxis,
                 angularaxis=fig2.layout.polar.angularaxis,
-            ),
+            )
         )
 
         return fig
diff --git a/ross/stochastic/st_shaft_element.py b/ross/stochastic/st_shaft_element.py
index d735ddf5..9963da2e 100644
--- a/ross/stochastic/st_shaft_element.py
+++ b/ross/stochastic/st_shaft_element.py
@@ -338,6 +338,6 @@ def st_shaft_example():
     from ross.materials import steel
 
     elm = ST_ShaftElement(
-        L=[1.0, 1.1], idl=0.0, odl=[0.1, 0.2], material=steel, is_random=["L", "odl"],
+        L=[1.0, 1.1], idl=0.0, odl=[0.1, 0.2], material=steel, is_random=["L", "odl"]
     )
     return elm


Test patch:
diff --git a/ross/tests/test_bearing_seal_element.py b/ross/tests/test_bearing_seal_element.py
index 2951e5fb..ced115bf 100644
--- a/ross/tests/test_bearing_seal_element.py
+++ b/ross/tests/test_bearing_seal_element.py
@@ -8,7 +8,7 @@ import pytest
 from numpy.testing import assert_allclose
 
 from ross.bearing_seal_element import (BallBearingElement, BearingElement,
-                                       BearingElement6DoF,
+                                       BearingElement6DoF, BearingFluidFlow,
                                        MagneticBearingElement,
                                        RollerBearingElement)
 
@@ -388,3 +388,42 @@ def test_save_load(bearing0, bearing_constant, bearing_6dof):
     bearing_6dof.save(file)
     bearing_6dof_loaded = BearingElement6DoF.load(file)
     assert bearing_6dof == bearing_6dof_loaded
+
+
+def test_bearing_fluid_flow():
+    nz = 30
+    ntheta = 20
+    length = 0.03
+    omega = [157.1]
+    p_in = 0.0
+    p_out = 0.0
+    radius_rotor = 0.0499
+    radius_stator = 0.05
+    load = 525
+    visc = 0.1
+    rho = 860.0
+    bearing = BearingFluidFlow(
+        0,
+        nz,
+        ntheta,
+        length,
+        omega,
+        p_in,
+        p_out,
+        radius_rotor,
+        radius_stator,
+        visc,
+        rho,
+        load=load,
+    )
+
+    # fmt: off
+    K = np.array([[14547442.70620538, 15571505.36655864],
+                  [-25596382.88167763, 12526684.40342712]])
+
+    C = np.array([[ 263025.76330117, -128749.90335233],
+                  [ -41535.76386708,  309417.62615761]])
+    # fmt: on
+
+    assert_allclose(bearing.K(0), K, rtol=1e-3)
+    assert_allclose(bearing.C(0), C, rtol=1e-3)
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
index 40402f70..755f4d4e 100644
--- a/ross/tests/test_fluid_flow.py
+++ b/ross/tests/test_fluid_flow.py
@@ -320,7 +320,7 @@ def test_plots():
 
 def test_find_equilibrium_position():
     bearing = fluid_flow_example2()
-    (n, t, force_x, force_y,) = calculate_oil_film_force(bearing)
+    (n, t, force_x, force_y) = calculate_oil_film_force(bearing)
     assert math.isclose(force_x, 0, abs_tol=1e-4)
     assert math.isclose(force_y, bearing.load, abs_tol=1e-2)
 
diff --git a/ross/tests/test_materials.py b/ross/tests/test_materials.py
index 6e9a23d7..74c444e8 100644
--- a/ross/tests/test_materials.py
+++ b/ross/tests/test_materials.py
@@ -7,7 +7,7 @@ from ross.materials import *
 @pytest.fixture
 def AISI4140():
     return Material(
-        name="AISI4140", rho=7850, E=203200000000.0, G_s=80000000000.0, color="#525252",
+        name="AISI4140", rho=7850, E=203200000000.0, G_s=80000000000.0, color="#525252"
     )
 
 
diff --git a/ross/tests/test_stochastic_elements.py b/ross/tests/test_stochastic_elements.py
index c8e8df3d..f800b4ad 100644
--- a/ross/tests/test_stochastic_elements.py
+++ b/ross/tests/test_stochastic_elements.py
@@ -33,7 +33,7 @@ def rand_shaft():
 @pytest.fixture
 def rand_disk_from_inertia():
     elm = ST_DiskElement(
-        n=1, m=[30, 40], Id=[0.2, 0.3], Ip=[0.5, 0.7], is_random=["m", "Id", "Ip"],
+        n=1, m=[30, 40], Id=[0.2, 0.3], Ip=[0.5, 0.7], is_random=["m", "Id", "Ip"]
     )
     return elm
 
diff --git a/ross/tests/test_stochastic_rotor_assembly.py b/ross/tests/test_stochastic_rotor_assembly.py
index ef8e5ab3..536c0b12 100644
--- a/ross/tests/test_stochastic_rotor_assembly.py
+++ b/ross/tests/test_stochastic_rotor_assembly.py
@@ -48,7 +48,7 @@ def rotor1():
 ###############################################################################
 def test_st_shaft_elements_odd_length():
     tim0 = ST_ShaftElement(
-        L=[1, 1.1], idl=0, odl=[0.1, 0.2], material=steel, is_random=["L", "odl"],
+        L=[1, 1.1], idl=0, odl=[0.1, 0.2], material=steel, is_random=["L", "odl"]
     )
     tim1 = ST_ShaftElement(
         L=[1, 1.1, 1.2],
@@ -78,17 +78,17 @@ def test_st_disk_elements_odd_length():
 
 
 def test_st_bearing_elements_odd_length():
-    tim0 = ShaftElement(L=0.25, idl=0, odl=0.05, material=steel,)
-    tim1 = ShaftElement(L=0.25, idl=0, odl=0.05, material=steel,)
+    tim0 = ShaftElement(L=0.25, idl=0, odl=0.05, material=steel)
+    tim1 = ShaftElement(L=0.25, idl=0, odl=0.05, material=steel)
     shaft_elm = [tim0, tim1]
 
     disk0 = DiskElement(n=1, m=20, Id=1, Ip=1)
 
     brg0 = ST_BearingElement(
-        n=0, kxx=[1e6, 2e6], cxx=[1e3, 2e3], is_random=["kxx", "cxx"],
+        n=0, kxx=[1e6, 2e6], cxx=[1e3, 2e3], is_random=["kxx", "cxx"]
     )
     brg1 = ST_BearingElement(
-        n=2, kxx=[1e6, 2e6, 3e6], cxx=[1e3, 2e3, 3e3], is_random=["kxx", "cxx"],
+        n=2, kxx=[1e6, 2e6, 3e6], cxx=[1e3, 2e3, 3e3], is_random=["kxx", "cxx"]
     )
 
     with pytest.raises(ValueError) as ex:
@@ -118,7 +118,7 @@ def test_st_point_mass_elements_odd_length():
 
 def test_elements_odd_length():
     tim0 = ST_ShaftElement(
-        L=[1, 1.1], idl=0, odl=[0.1, 0.2], material=steel, is_random=["L", "odl"],
+        L=[1, 1.1], idl=0, odl=[0.1, 0.2], material=steel, is_random=["L", "odl"]
     )
     shaft_elm = [tim0, tim0]
 
@@ -127,10 +127,10 @@ def test_elements_odd_length():
     disks = [disk0, disk1]
 
     brg0 = ST_BearingElement(
-        n=0, kxx=[1e6, 2e6], cxx=[1e3, 2e3], is_random=["kxx", "cxx"],
+        n=0, kxx=[1e6, 2e6], cxx=[1e3, 2e3], is_random=["kxx", "cxx"]
     )
     brg1 = ST_BearingElement(
-        n=2, kxx=[1e6, 2e6], cxx=[1e3, 2e3], is_random=["kxx", "cxx"],
+        n=2, kxx=[1e6, 2e6], cxx=[1e3, 2e3], is_random=["kxx", "cxx"]
     )
     bearings = [brg0, brg1]