Repository: ross-rotordynamics/ross (python)

Issue:
Include Sommerfeld number in the fluid flow class
According to Friswell et al. in Dynamics of Rotating Machines (2010, page 178), if the magnitude of the load is known, then the bearing eccentricity may be obtained by solving a quartic equation including the Sommerfeld number. The program must be able to calculate the Sommerfeld number and give the eccentricity of the bearing.

Pull request description:
Add Sommerfeld number function
Add a function to calculate the Sommerfeld number for a given load. Also add two more functions: one to calculate eccentricity for a given load using the Sommerfeld number, and one to calculate the load using the eccentricity given when the class is instantiated.
Close #101 

Patch:
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index d35952ca..f6414f35 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -111,6 +111,8 @@ class PressureMatrix:
     plot_counter: int
         for each plot, it increments this counter, in order to plot
         figures with unique IDs.
+    radial_clearance: float
+        Difference between both stator and rotor radius, regardless of eccentricity.
 
     Examples
     --------
@@ -177,6 +179,7 @@ class PressureMatrix:
         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
@@ -258,6 +261,51 @@ class PressureMatrix:
 
         return radius_external, xre, yre
 
+    def get_rotor_load(self):
+        """Returns the load applied to the rotor.
+        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 sommerfeld_number(self, f):
+        """Return the sommerfeld number.
+        Parameters
+        ----------
+        f: float
+            Load applied to the rotor.
+        Returns
+        -------
+        float
+            The sommerfeld number.
+        """
+        modified_s = (self.radius_stator*2*self.omega*self.visc*(self.length**3)) / \
+                     (8*f*(self.radial_clearance**2))
+        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.
+        Returns
+        -------
+        float
+            The eccentricity ratio.
+        """
+        s = self.sommerfeld_number(f)
+        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)):
+            if 0 <= roots[i] <= 1:
+                return np.sqrt(roots[i])
+        sys.exit("Eccentricity ratio could not be calculated.")
+
     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.


Test patch:
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
new file mode 100644
index 00000000..fddbfdaf
--- /dev/null
+++ b/ross/tests/test_fluid_flow.py
@@ -0,0 +1,28 @@
+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)
+    """
+    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
+    eccentricity = 0
+    visc = 0.1
+    rho = 860.
+    bearing = flow.PressureMatrix(nz, ntheta, nradius, length, omega, p_in,
+                                  p_out, radius_rotor, radius_stator,
+                                  eccentricity, visc, rho)
+    s = bearing.sommerfeld_number(525)
+    assert math.isclose(s, 3.571, rel_tol=0.001)
+