Repository: ross-rotordynamics/ross
Pull request #119 (2019-05-18 02:56:20)

Add bearing instantiation from fluid flow
A bearing can now be instantiated from fluid flow by giving fluid flow parameters to a bearing method.


Patch:
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index ca53873d..f25d0ed8 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -9,6 +9,8 @@ import pytest
 import pandas as pd
 import sys
 import xlrd
+from ross.fluid_flow import fluid_flow as flow
+
 
 __all__ = ["BearingElement", "SealElement"]
 bokeh_colors = bp.RdGy[11]
@@ -264,6 +266,7 @@ class BearingElement(Element):
 
             bearing_elements.append(bearing)
         return bearing_elements
+
     def M(self):
         M = np.zeros((4, 4))
 
@@ -350,11 +353,14 @@ class BearingElement(Element):
         to a toml file that can be later loaded by ross.
         Parameters
         ----------
-        :param n : int. The node in which the bearing will be located in the rotor.
-        :param file: str. Path to the file containing the bearing parameters.
+        n : int
+            The node in which the bearing will be located in the rotor.
+        file: str
+            Path to the file containing the bearing parameters.
         Returns
         -------
-        :return: dict that is ready to save to toml and readable by ross.
+        dict
+            A dict that is ready to save to toml and readable by ross.
         """
         b_elem = cls.from_table(n, file)
         data = {
@@ -376,11 +382,13 @@ class BearingElement(Element):
         """Instantiate a bearing using inputs from a table, either excel or csv.
         Parameters
         ----------
-        :param n : int. The node in which the bearing will be located in the rotor.
-        :param file: str. Path to the file containing the bearing parameters.
+        n : int
+            The node in which the bearing will be located in the rotor.
+        file: str
+            Path to the file containing the bearing parameters.
         Returns
         -------
-        :return: A bearing object.
+        A bearing object.
         """
         try:
             df = pd.read_excel(file)
@@ -403,6 +411,70 @@ class BearingElement(Element):
                      "Make sure the table header contains the parameters for the "
                      "BearingElement class.")
 
+    @classmethod
+    def from_fluid_flow(cls, n, nz, ntheta, nradius, length, omega, p_in,
+                        p_out, radius_rotor, radius_stator, visc, rho,
+                        eccentricity=None, load=None):
+        """Instantiate a bearing using inputs from its fluid flow.
+        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: float
+            Rotation of the rotor (rad/s).
+        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).
+        Returns
+        -------
+        A bearing object.
+        """
+        fluid_flow = flow.PressureMatrix(nz, ntheta, nradius, length, omega, p_in,
+                                         p_out, radius_rotor, radius_stator,
+                                         visc, rho, eccentricity=eccentricity, load=load)
+        k = fluid_flow.get_analytical_damping_matrix()
+        c = fluid_flow.get_analytical_stiffness_matrix()
+        return cls(n, kxx=k[0], cxx=c[0], kyy=k[3],
+                   kxy=k[1], kyx=k[2], cyy=c[3],
+                   cxy=c[1], cyx=c[2], w=fluid_flow.omega)
+
 
 class SealElement(BearingElement):
     def __init__(
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index 8cc4c844..3d67805f 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -2,6 +2,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 import sys
 from bokeh.plotting import figure, output_file, show
+import warnings
 
 
 class PressureMatrix:
@@ -33,6 +34,8 @@ class PressureMatrix:
         Input Pressure (Pa).
     p_out: float
         Output Pressure (Pa).
+    load: float
+        Load applied to the rotor (N).
 
     Geometric data of the problem
     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -132,7 +135,7 @@ class PressureMatrix:
     >>> rho = 860.
     >>> my_pressure_matrix = flow.PressureMatrix(nz, ntheta, nradius, length,
     ...                                          omega, p_in, p_out, radius_rotor,
-    ...                                          radius_stator, eccentricity,  visc, rho)
+    ...                                          radius_stator, visc, rho, eccentricity=eccentricity)
     >>> my_pressure_matrix.calculate_pressure_matrix()
     >>> my_pressure_matrix.plot_eccentricity()
     >>> my_pressure_matrix.plot_pressure_z()
@@ -142,8 +145,9 @@ class PressureMatrix:
 
     """
     def __init__(self, nz, ntheta, nradius, length, omega, p_in,
-                 p_out, radius_rotor, radius_stator, eccentricity,
-                 visc, rho):
+                 p_out, radius_rotor, radius_stator, visc, rho, eccentricity=None, load=None):
+        if load is None and eccentricity is None:
+            sys.exit("Either load or eccentricity must be given.")
         self.nz = nz
         self.ntheta = ntheta
         self.nradius = nradius
@@ -160,11 +164,27 @@ class PressureMatrix:
         self.p_out = p_out
         self.radius_rotor = radius_rotor
         self.radius_stator = radius_stator
-        self.eccentricity = eccentricity
-        self.xi = np.sqrt(2)*eccentricity/2
-        self.yi = -self.xi
         self.visc = visc
         self.rho = rho
+        self.radial_clearance = self.radius_stator - self.radius_rotor
+        self.difference_between_radius = radius_stator - radius_rotor
+        self.bearing_type = ''
+        if self.length/self.radius_stator <= 1/8:
+            self.bearing_type = 'short_bearing'
+        elif self.length/self.radius_stator > 4:
+            self.bearing_type = 'long_bearing'
+        else:
+            self.bearing_type = 'medium_size'
+        self.eccentricity = eccentricity
+        self.eccentricity_ratio = None
+        self.load = load
+        if self.eccentricity is None:
+            self.eccentricity = self.calculate_eccentricity_ratio()*self.difference_between_radius
+        self.eccentricity_ratio = self.eccentricity / self.difference_between_radius
+        if self.load is None:
+            self.load = self.get_rotor_load()
+        self.xi = np.sqrt(2)*self.eccentricity/2
+        self.yi = -self.xi
         self.re = np.zeros([self.nz, self.ntheta])
         self.ri = np.zeros([self.nz, self.ntheta])
         self.z = np.zeros([1, self.nz])
@@ -173,13 +193,9 @@ class PressureMatrix:
         self.yre = np.zeros([self.nz, self.ntheta])
         self.yri = np.zeros([self.nz, self.ntheta])
         self.p_mat = np.zeros([self.nz, self.ntheta])
-        self.bearing_type = ''
         self.plot_counter = 0
         self.calculate_coefficients()
         self.pressure_matrix_available = False
-        self.difference_between_radius = radius_stator - radius_rotor
-        self.eccentricity_ratio = self.eccentricity/self.difference_between_radius
-        self.radial_clearance = self.radius_stator - self.radius_rotor
 
     def calculate_pressure_matrix(self):
         """This function calculates the pressure matrix
@@ -199,12 +215,6 @@ class PressureMatrix:
         of the discrete pressure (central differences in the second
         derivatives). It is executed when the class is instantiated.
         """
-        if self.length/self.radius_stator <= 1/8:
-            self.bearing_type = 'short_bearing'
-        elif self.length/self.radius_stator > 4:
-            self.bearing_type = 'long_bearing'
-        else:
-            self.bearing_type = 'medium_size'
         for i in range(self.nz):
             zno = i * self.dz
             self.z[0][i] = zno
@@ -265,65 +275,64 @@ class PressureMatrix:
         yre: float
             The position y of the returned external radius.
         """
-        betha = -np.pi/4
-        alpha = betha - gama
-        radius_external = np.sqrt(self.radius_stator**2 - (self.eccentricity * np.sin(alpha))**2) + self.eccentricity * np.cos(alpha)
+        beta = -np.pi/4
+        alpha = beta - gama
+        radius_external = (np.sqrt(self.radius_stator**2 - (self.eccentricity * np.sin(alpha))**2) +
+                           self.eccentricity * np.cos(alpha))
         xre = radius_external * np.cos(gama)
         yre = radius_external * np.sin(gama)
 
         return radius_external, xre, yre
 
     def get_rotor_load(self):
-        """Returns the load applied to the rotor.
+        """Returns the load applied to the rotor, based on the eccentricity ratio.
+        Suitable only for short bearings.
         Returns
         -------
         float
             Load applied to the rotor.
         """
-        return -((np.pi*self.radius_stator*2*self.omega*self.visc*(self.length**3)*self.eccentricity_ratio)
-                 / (8*(self.radial_clearance**2)*((1 - self.eccentricity_ratio**2)**2))) \
-            * np.sqrt((16/np.pi - 1)*self.eccentricity_ratio + 1)
-
-    def modified_sommerfeld_number(self, f):
-        """Return the modified sommerfeld number.
-        Parameters
-        ----------
-        f: float
-            Load applied to the rotor.
+        if not self.bearing_type == 'short_bearing':
+            warnings.warn("Function get_rotor_load suitable only for short bearings. "
+                          "The ratio between the bearing length and its radius should be less or "
+                          "equal to 0.125. Currently we have "+str(self.length/self.radius_stator)+".")
+        return ((np.pi*self.radius_stator*2*self.omega*self.visc*(self.length**3)*self.eccentricity_ratio)
+                / (8*(self.radial_clearance**2)*((1 - self.eccentricity_ratio**2)**2))) \
+            * (np.sqrt((16/(np.pi**2) - 1)*self.eccentricity_ratio**2 + 1))
+
+    def modified_sommerfeld_number(self):
+        """Returns the modified sommerfeld number.
         Returns
         -------
         float
             The modified sommerfeld number.
         """
         return (self.radius_stator*2*self.omega*self.visc*(self.length**3)) / \
-               (8*f*(self.radial_clearance**2))
+               (8*self.load*(self.radial_clearance**2))
 
-    def sommerfeld_number(self, f):
-        """Return the sommerfeld number.
-        Parameters
-        ----------
-        f: float
-            Load applied to the rotor.
+    def sommerfeld_number(self):
+        """Returns the sommerfeld number.
         Returns
         -------
         float
             The sommerfeld number.
         """
-        modified_s = self.modified_sommerfeld_number(f)
+        modified_s = self.modified_sommerfeld_number()
         return (modified_s/np.pi)*(self.radius_stator*2/self.length)**2
 
-    def calculate_eccentricity_ratio(self, f):
-        """Calculate the eccentricity ratio using the sommerfeld number.
-        Parameters
-        ----------
-        f: float
-            Load applied to the rotor.
+    def calculate_eccentricity_ratio(self):
+        """Calculate the eccentricity ratio for a given load using the sommerfeld number.
+        Suitable only for short bearings.
         Returns
         -------
         float
             The eccentricity ratio.
         """
-        s = self.modified_sommerfeld_number(f)
+        if not self.bearing_type == 'short_bearing':
+            warnings.warn("Function calculate_eccentricity_ratio suitable only for short bearings. "
+                          "The ratio between the bearing length and its radius should be less or "
+                          "equal to 0.125. Currently we have "+str(self.length/self.radius_stator)+".")
+        s = self.modified_sommerfeld_number()
         coefficients = [1, -4, (6 - (s**2)*(16 - np.pi**2)), -(4 + (np.pi**2)*(s**2)), 1]
         roots = np.roots(coefficients)
         for i in range(0, len(roots)):
@@ -331,6 +340,57 @@ class PressureMatrix:
                 return np.sqrt(roots[i])
         sys.exit("Eccentricity ratio could not be calculated.")
 
+    def get_analytical_stiffness_matrix(self):
+        """Returns the stiffness matrix calculated analytically.
+        Suitable only for short bearings.
+        Returns
+        -------
+        list of floats
+            A list of length four including stiffness floats in this order: kxx, kxy, kyx, kyy
+        """
+        if not self.bearing_type == 'short_bearing':
+            warnings.warn("Function get_analytical_stiffness_matrix suitable only for short bearings. "
+                          "The ratio between the bearing length and its radius should be less or "
+                          "equal to 0.125. Currently we have "+str(self.length/self.radius_stator)+".")
+        f = self.get_rotor_load()
+        h0 = 1.0/(((np.pi**2)*(1 - self.eccentricity_ratio**2) + 16*self.eccentricity_ratio**2)**1.5)
+        a = f/self.radial_clearance
+        kxx = a*h0*4*((np.pi**2)*(2 - self.eccentricity_ratio**2) + 16*self.eccentricity_ratio**2)
+        kxy = (a*h0*np.pi*((np.pi**2)*(1 - self.eccentricity_ratio**2)**2 -
+               16*self.eccentricity_ratio**4) /
+               (self.eccentricity_ratio*np.sqrt(1 - self.eccentricity_ratio**2)))
+        kyx = (-a*h0*np.pi*((np.pi**2)*(1 - self.eccentricity_ratio**2) *
+               (1 + 2*self.eccentricity_ratio**2) +
+               (32*self.eccentricity_ratio**2)*(1 + self.eccentricity_ratio**2)) /
+               (self.eccentricity_ratio*np.sqrt(1 - self.eccentricity_ratio**2)))
+        kyy = (a*h0*4*((np.pi**2)*(1 + 2*self.eccentricity_ratio**2) +
+                                  ((32*self.eccentricity_ratio**2) *
+                                   (1 + self.eccentricity_ratio**2))/(1 - self.eccentricity_ratio**2)))
+        return [kxx, kxy, kyx, kyy]
+
+    def get_analytical_damping_matrix(self):
+        """Returns the damping matrix calculated analytically.
+        Suitable only for short bearings.
+        Returns
+        -------
+        list of floats
+            A list of length four including stiffness floats in this order: cxx, cxy, cyx, cyy
+        """
+        if not self.bearing_type == 'short_bearing':
+            warnings.warn("Function get_analytical_damping_matrix suitable only for short bearings. "
+                          "The ratio between the bearing length and its radius should be less or "
+                          "equal to 0.125. Currently we have "+str(self.length/self.radius_stator)+".")
+        f = self.get_rotor_load()
+        h0 = 1.0/(((np.pi**2)*(1 - self.eccentricity_ratio**2) + 16*self.eccentricity_ratio**2)**1.5)
+        a = f/(self.radial_clearance * self.omega)
+        cxx = (a*h0*2*np.pi*np.sqrt(1 - self.eccentricity_ratio**2) *
+               ((np.pi**2)*(1 + 2*self.eccentricity_ratio**2) - 16*self.eccentricity_ratio**2)/self.eccentricity_ratio)
+        cxy = (-a*h0*8*((np.pi**2)*(1 + 2*self.eccentricity_ratio**2) - 16*self.eccentricity_ratio**2))
+        cyx = cxy
+        cyy = (a*h0*(2*np.pi*((np.pi**2)*(1 - self.eccentricity_ratio**2)**2 + 48*self.eccentricity_ratio**2)) /
+               (self.eccentricity_ratio*np.sqrt(1 - self.eccentricity_ratio**2)))
+        return [cxx, cxy, cyx, cyy]
+
     def plot_eccentricity(self, z=0):
         """This function assembles pressure graphic along the z-axis.
         The first few plots are of a different color to indicate where theta begins.
diff --git a/ross/results.py b/ross/results.py
index e8a21608..d46f3171 100644
--- a/ross/results.py
+++ b/ross/results.py
@@ -484,8 +484,6 @@ class ForcedResponseResults(Results):
         -------
         ax : matplotlib.axes
             Matplotlib axes with phase plot.
-        mag_plot : bokeh axes
-            bokeh axes with magnitude plot
         Examples
         --------
         """
@@ -497,13 +495,10 @@ class ForcedResponseResults(Results):
 
         if units == "m":
             ax.set_ylabel("Amplitude $(m)$")
-            y_axis_label = "Amplitude (m)"
         elif units == "mic-pk-pk":
             mag = 2 * mag * 1e6
             ax.set_ylabel("Amplitude $(\mu pk-pk)$")
-            y_axis_label = "Amplitude $(\mu pk-pk)$"
 
-        # matplotlib plotting
         ax.plot(frequency_range, mag[dof], **kwargs)
 
         ax.set_xlim(0, max(frequency_range))
@@ -513,26 +508,7 @@ class ForcedResponseResults(Results):
         ax.set_xlabel("Frequency (rad/s)")
         ax.legend()
 
-        # bokeh plot - create a new plot
-        mag_plot = figure(
-                tools="pan, box_zoom, wheel_zoom, reset, save",
-                width=900,
-                height=400,
-                title="Forced Response - Magnitude",
-                x_axis_label="Frequency",
-                x_range=[0, max(frequency_range)],
-                y_axis_label=y_axis_label
-        )
-        source = ColumnDataSource(dict(x=frequency_range, y=mag[dof]))
-        mag_plot.line(x='x',
-                      y='y',
-                      source=source,
-                      line_color=bokeh_colors[0],
-                      line_alpha=1.0,
-                      line_width=3,
-                      )
-
-        return ax, mag_plot
+        return ax
 
     def plot_phase(self, dof, ax=None, **kwargs):
         """Plot frequency response.
@@ -550,10 +526,10 @@ class ForcedResponseResults(Results):
             the plot (e.g. linestyle='--')
         Returns
         -------
-        ax : matplotlib.axes
+        ax0 : matplotlib.axes
+            Matplotlib axes with amplitude plot.
+        ax1 : matplotlib.axes
             Matplotlib axes with phase plot.
-        phase_plot : bokeh axes
-            Bokeh axes with phase plot
         Examples
         --------
         """
@@ -571,28 +547,10 @@ class ForcedResponseResults(Results):
 
         ax.set_ylabel("Phase")
         ax.set_xlabel("Frequency (rad/s)")
-        ax.legend()
 
-        # bokeh plot - create a new plot
-        phase_plot = figure(
-                tools="pan, box_zoom, wheel_zoom, reset, save",
-                width=900,
-                height=400,
-                title="Forced Response - Magnitude",
-                x_axis_label="Frequency",
-                x_range=[0, max(frequency_range)],
-                y_axis_label="Phase"
-        )
-        source = ColumnDataSource(dict(x=frequency_range, y=phase[dof]))
-        phase_plot.line(x='x',
-                        y='y',
-                        source=source,
-                        line_color=bokeh_colors[0],
-                        line_alpha=1.0,
-                        line_width=3,
-                        )
+        ax.legend()
 
-        return ax, phase_plot
+        return ax
 
     def plot(self, dof, ax0=None, ax1=None, **kwargs):
         """Plot frequency response.
@@ -617,34 +575,22 @@ class ForcedResponseResults(Results):
             Matplotlib axes with amplitude plot.
         ax1 : matplotlib.axes
             Matplotlib axes with phase plot.
-        bk_ax0 : bokeh axes
-            Bokeh axes with amplitude plot.
-        bk_ax1 : bokeh axes
-            Bokeh axes with phase plot.
         Examples
         --------
         """
         if ax0 is None and ax1 is None:
             fig, (ax0, ax1) = plt.subplots(2)
 
-        ax0 = self.plot_magnitude(dof, ax=ax0, **kwargs)[0]
+        ax0 = self.plot_magnitude(dof, ax=ax0, **kwargs)
         # remove label from phase plot
         kwargs.pop("label", None)
         kwargs.pop("units", None)
-        ax1 = self.plot_phase(dof, ax=ax1, **kwargs)[0]
+        ax1 = self.plot_phase(dof, ax=ax1, **kwargs)
 
         ax0.set_xlabel("")
         ax0.legend()
 
-        # bokeh plot axes
-        bk_ax0 = self.plot_magnitude(dof, ax=ax0, **kwargs)[1]
-        bk_ax1 = self.plot_phase(dof, ax=ax1, **kwargs)[1]
-
-        # show the bokeh plot results
-        grid_plots = gridplot([[bk_ax0], [bk_ax1]])
-        show(grid_plots)
-
-        return ax0, ax1, bk_ax0, bk_ax1
+        return ax0, ax1
 
 
 class ModeShapeResults(Results):


Test patch:
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
index 674c57dc..e2ea9ed8 100644
--- a/ross/tests/test_fluid_flow.py
+++ b/ross/tests/test_fluid_flow.py
@@ -2,12 +2,7 @@ from ross.fluid_flow import fluid_flow as flow
 import math
 
 
-def test_sommerfeld_number():
-    """
-    This function instantiate a bearing using the fluid flow class and test if it matches the
-    expected results for the sommerfeld number.
-    Taken from example 5.5.1, page 181 (Dynamics of rotating machine, FRISSWELL)
-    """
+def fluid_flow_eccentricity():
     nz = 30
     ntheta = 20
     nradius = 11
@@ -17,12 +12,75 @@ def test_sommerfeld_number():
     p_out = 0.
     radius_rotor = 0.0499
     radius_stator = 0.05
-    eccentricity = 0
+    eccentricity = (radius_stator - radius_rotor)*0.2663
     visc = 0.1
     rho = 860.
-    bearing = flow.PressureMatrix(nz, ntheta, nradius, length, omega, p_in,
-                                  p_out, radius_rotor, radius_stator,
-                                  eccentricity, visc, rho)
-    e = bearing.calculate_eccentricity_ratio(525)
-    assert math.isclose(e, 0.2663, rel_tol=0.001)
+    return flow.PressureMatrix(nz, ntheta, nradius, length, omega, p_in,
+                               p_out, radius_rotor, radius_stator,
+                               visc, rho, eccentricity=eccentricity)
+
+
+def fluid_flow_load():
+    nz = 30
+    ntheta = 20
+    nradius = 11
+    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.
+    return flow.PressureMatrix(nz, ntheta, nradius, length, omega, p_in,
+                               p_out, radius_rotor, radius_stator,
+                               visc, rho, load=load)
+
+
+def test_sommerfeld_number():
+    """
+    This function instantiate a bearing using the fluid flow class and test if it matches the
+    expected results for the sommerfeld number and eccentricity ratio.
+    Taken from example 5.5.1, page 181 (Dynamics of rotating machine, FRISSWELL)
+    """
+    bearing = fluid_flow_load()
+    assert math.isclose(bearing.eccentricity_ratio, 0.2663, rel_tol=0.001)
+
+
+def test_get_rotor_load():
+    """
+    This function instantiate a bearing using the fluid flow class and test if it matches the
+    expected results for the load over the rotor, given the eccentricity ratio.
+    Taken from example 5.5.1, page 181 (Dynamics of rotating machine, FRISSWELL)
+    """
+    bearing = fluid_flow_eccentricity()
+    assert math.isclose(bearing.load, 525, rel_tol=0.1)
+
 
+def test_stiffness_matrix():
+    """
+    This function instantiate a bearing using the fluid flow class and test if it matches the
+    expected results for the stiffness matrix, given the eccentricity ratio.
+    Taken from example 5.5.1, page 181 (Dynamics of rotating machine, FRISSWELL)
+    """
+    bearing = fluid_flow_eccentricity()
+    kxx, kxy, kyx, kyy = bearing.get_analytical_stiffness_matrix()
+    assert math.isclose(kxx/10**6, 12.81, rel_tol=0.01)
+    assert math.isclose(kxy/10**6, 16.39, rel_tol=0.01)
+    assert math.isclose(kyx/10**6, -25.06, rel_tol=0.01)
+    assert math.isclose(kyy/10**6, 8.815, rel_tol=0.01)
+
+
+def test_damping_matrix():
+    """
+    This function instantiate a bearing using the fluid flow class and test if it matches the
+    expected results for the damping matrix, given the eccentricity ratio.
+    Taken from example 5.5.1, page 181 (Dynamics of rotating machine, FRISSWELL)
+    """
+    bearing = fluid_flow_load()
+    cxx, cxy, cyx, cyy = bearing.get_analytical_damping_matrix()
+    assert math.isclose(cxx/10**3, 232.9, rel_tol=0.01)
+    assert math.isclose(cxy/10**3, -81.92, rel_tol=0.01)
+    assert math.isclose(cyx/10**3, -81.92, rel_tol=0.01)
+    assert math.isclose(cyy/10**3, 294.9, rel_tol=0.01)