Repository: ross-rotordynamics/ross
Pull request #659 (2020-08-21 23:43:55)

Find equilibrium position
In this Pull Request, the equilibrium position for long bearings is calculated by means of an iterative method. It should work well for bearings with non-cylindrical geometry.


Patch:
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index e1d52ae7..7fce4f16 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -918,7 +918,6 @@ class BearingElement(Element):
         n,
         nz,
         ntheta,
-        nradius,
         length,
         omega,
         p_in,
@@ -989,26 +988,24 @@ class BearingElement(Element):
         --------
         >>> 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 = (radius_stator - radius_rotor)*0.2663
+        >>> load = 525
         >>> visc = 0.1
         >>> rho = 860.
-        >>> BearingElement.from_fluid_flow(0, nz, ntheta, nradius, length, omega, p_in,
+        >>> BearingElement.from_fluid_flow(0, nz, ntheta, length, omega, p_in,
         ...                                p_out, radius_rotor, radius_stator,
-        ...                                visc, rho, eccentricity=eccentricity) # doctest: +ELLIPSIS
+        ...                                visc, rho, load=load) # doctest: +ELLIPSIS
         BearingElement(n=0, n_link=None,
          kxx=[...
         """
         fluid_flow = flow.FluidFlow(
             nz,
             ntheta,
-            nradius,
             length,
             omega,
             p_in,
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index 27b8760b..1bfcd572 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -4,6 +4,7 @@ import sys
 import numpy as np
 import scipy as sp
 
+from ross.fluid_flow.fluid_flow_coefficients import find_equilibrium_position
 from ross.fluid_flow.fluid_flow_geometry import (calculate_attitude_angle,
                                                  calculate_eccentricity_ratio,
                                                  calculate_rotor_load,
@@ -198,9 +199,12 @@ class FluidFlow:
         load=None,
         omegap=None,
         immediately_calculate_pressure_matrix_numerically=True,
+        bearing_type=None,
     ):
         if load is None and eccentricity is None:
             sys.exit("Either load or eccentricity must be given.")
+        calculate_equilibrium_position = load is not None
+
         self.nz = nz
         self.ntheta = ntheta
         self.n_interv_z = nz - 1
@@ -219,65 +223,83 @@ class FluidFlow:
         self.density = density
         self.characteristic_speed = self.omega * self.radius_rotor
         self.radial_clearance = self.radius_stator - self.radius_rotor
-        self.bearing_type = ""
-        if self.length / (2 * self.radius_stator) <= 1 / 4:
-            self.bearing_type = "short_bearing"
-        elif self.length / (2 * self.radius_stator) > 4:
-            self.bearing_type = "long_bearing"
-        else:
-            self.bearing_type = "medium_size"
+        self.bearing_type = bearing_type
+        if bearing_type is None:
+            if self.length / (2 * self.radius_stator) <= 1 / 4:
+                self.bearing_type = "short_bearing"
+            elif self.length / (2 * 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:
-            modified_s = modified_sommerfeld_number(
-                self.radius_stator,
-                self.omega,
-                self.viscosity,
-                self.length,
-                self.load,
-                self.radial_clearance,
-            )
-            self.eccentricity = (
-                calculate_eccentricity_ratio(modified_s) * self.radial_clearance
-            )
+
         self.omegap = omegap
         if self.omegap is None:
             self.omegap = self.omega
         else:
             self.omegap = omegap
-        self.eccentricity_ratio = self.eccentricity / self.radial_clearance
-        if self.load is None:
-            self.load = calculate_rotor_load(
-                self.radius_stator,
-                self.omega,
-                self.viscosity,
-                self.length,
-                self.radial_clearance,
-                self.eccentricity_ratio,
-            )
-        if attitude_angle is None:
-            self.attitude_angle = calculate_attitude_angle(self.eccentricity_ratio)
-        else:
-            self.attitude_angle = attitude_angle
-        self.xi = self.eccentricity * np.cos(3 * np.pi / 2 + self.attitude_angle)
-        self.yi = self.eccentricity * np.sin(3 * np.pi / 2 + self.attitude_angle)
+        self.z_list = np.zeros(self.nz)
         self.re = np.zeros([self.nz, self.ntheta])
         self.ri = np.zeros([self.nz, self.ntheta])
         self.xre = np.zeros([self.nz, self.ntheta])
         self.xri = np.zeros([self.nz, self.ntheta])
         self.yre = np.zeros([self.nz, self.ntheta])
         self.yri = np.zeros([self.nz, self.ntheta])
-        self.z_list = np.zeros(self.nz)
         self.gama = np.zeros([self.nz, self.ntheta])
         self.t = 0
         self.xp = 0
         self.yp = 0
+        if self.bearing_type == "short_bearing":
+            if self.eccentricity is None:
+                modified_s = modified_sommerfeld_number(
+                    self.radius_stator,
+                    self.omega,
+                    self.viscosity,
+                    self.length,
+                    self.load,
+                    self.radial_clearance,
+                )
+            if calculate_equilibrium_position is True:
+                self.eccentricity = (
+                    calculate_eccentricity_ratio(modified_s) * self.radial_clearance
+                )
+            self.eccentricity_ratio = self.eccentricity / self.radial_clearance
+            if attitude_angle is None:
+                self.attitude_angle = calculate_attitude_angle(self.eccentricity_ratio)
+            else:
+                self.attitude_angle = attitude_angle
+            self.xi = self.eccentricity * np.cos(3 * np.pi / 2 + self.attitude_angle)
+            self.yi = self.eccentricity * np.sin(3 * np.pi / 2 + self.attitude_angle)
+            self.geometry_description()
+            if self.load is None:
+                self.load = calculate_rotor_load(
+                    self.radius_stator,
+                    self.omega,
+                    self.viscosity,
+                    self.length,
+                    self.radial_clearance,
+                    self.eccentricity_ratio,
+                )
+        else:
+            if calculate_equilibrium_position is True:
+                find_equilibrium_position(self)
+            else:
+                if attitude_angle is None:
+                    sys.exit("Attitude angle must be given.")
+                else:
+                    self.attitude_angle = attitude_angle
+            self.geometry_description()
+            self.eccentricity_ratio = self.eccentricity / self.radial_clearance
+            self.xi = self.eccentricity * np.cos(3 * np.pi / 2 + self.attitude_angle)
+            self.yi = self.eccentricity * np.sin(3 * np.pi / 2 + self.attitude_angle)
+
         self.p_mat_analytical = np.zeros([self.nz, self.ntheta])
         self.p_mat_numerical = np.zeros([self.nz, self.ntheta])
-        self.geometry_description()
         self.analytical_pressure_matrix_available = False
         self.numerical_pressure_matrix_available = False
+
         if immediately_calculate_pressure_matrix_numerically:
             self.calculate_pressure_matrix_numerical()
 
@@ -675,5 +697,55 @@ def fluid_flow_example2():
         visc,
         rho,
         load=load,
+        immediately_calculate_pressure_matrix_numerically=True,
+    )
+
+
+def fluid_flow_example3():
+    """This function returns a different instance of a simple fluid flow.
+    The purpose is to make available a simple model
+    so that doctest can be written using it.
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+    An instance of a fluid flow object.
+
+    Examples
+    --------
+    >>> my_fluid_flow = fluid_flow_example3()
+    >>> my_fluid_flow.eccentricity
+    0.0001
+    """
+    from ross.fluid_flow import fluid_flow as flow
+
+    nz = 8
+    ntheta = 32 * 4
+    omega = 100.0 * 2 * np.pi / 60
+    p_in = 0.0
+    p_out = 0.0
+    radius_rotor = 1
+    h = 0.000194564
+    radius_stator = radius_rotor + h
+    length = 8 * 2 * radius_stator
+    visc = 0.015
+    rho = 860.0
+    eccentricity = 1e-4
+    attitude_angle = np.pi / 4
+    return flow.FluidFlow(
+        nz,
+        ntheta,
+        length,
+        omega,
+        p_in,
+        p_out,
+        radius_rotor,
+        radius_stator,
+        visc,
+        rho,
+        eccentricity=eccentricity,
+        attitude_angle=attitude_angle,
         immediately_calculate_pressure_matrix_numerically=False,
     )
diff --git a/ross/fluid_flow/fluid_flow_coefficients.py b/ross/fluid_flow/fluid_flow_coefficients.py
index 9813904c..767f391d 100644
--- a/ross/fluid_flow/fluid_flow_coefficients.py
+++ b/ross/fluid_flow/fluid_flow_coefficients.py
@@ -1,10 +1,15 @@
-import warnings
+import sys
 from math import isnan
 
 import numpy as np
 from scipy import integrate
+from scipy.optimize import least_squares
 
-from ross.fluid_flow.fluid_flow_geometry import move_rotor_center
+# fmt: off
+from ross.fluid_flow.fluid_flow_geometry import (move_rotor_center,
+                                                 move_rotor_center_abs)
+
+# fmt: on
 
 
 def calculate_oil_film_force(fluid_flow_object, force_type=None):
@@ -367,192 +372,78 @@ def calculate_short_damping_matrix(fluid_flow_object):
     return [cxx, cxy, cyx, cyy]
 
 
-def find_equilibrium_position(
-    fluid_flow_object,
-    print_along=True,
-    tolerance=1e-05,
-    increment_factor=1e-03,
-    max_iterations=10,
-    increment_reduction_limit=1e-04,
-    return_iteration_map=False,
-):
-    """This function returns an eccentricity value with calculated forces matching the load applied,
-    meaning an equilibrium position of the rotor.
-    It first moves the rotor center on x-axis, aiming for the minimum error in the force on x (zero), then
-    moves on y-axis, aiming for the minimum error in the force on y (meaning load minus force on y equals zero).
+def find_equilibrium_position(fluid_flow_object, print_equilibrium_position=False):
+    """This function finds the equilibrium position of the rotor such that the fluid flow
+    forces match the applied load.
     Parameters
     ----------
     fluid_flow_object: A FluidFlow object.
-    print_along: bool, optional
-        If True, prints the iteration process.
-    tolerance: float, optional
-    increment_factor: float, optional
-        This number will multiply the first eccentricity found to reach an increment number.
-    max_iterations: int, optional
-    increment_reduction_limit: float, optional
-        The error should always be approximating zero. If it passes zeros (for instance, from a positive error
-        to a negative one), the iteration goes back one step and the increment is reduced. This reduction must
-        have a limit to avoid long iterations.
-    return_iteration_map: bool, optional
-        If True, along with the eccentricity found, the function will return a map of position and errors in
-        each step of the iteration.
+    print_equilibrium_position: bool, optional
+        If True, prints the equilibrium position.
     Returns
     -------
-    None, or
-    Matrix of floats
-        A matrix [4, n], being n the number of iterations. In each line, it contains the x and y of the rotor
-        center, followed by the error in force x and force y.
-    Examples
+    None
     --------
     >>> from ross.fluid_flow.fluid_flow import fluid_flow_example2
     >>> my_fluid_flow = fluid_flow_example2()
-    >>> find_equilibrium_position(my_fluid_flow, print_along=False,
-    ...                           tolerance=0.1, increment_factor=0.01,
-    ...                           max_iterations=5, increment_reduction_limit=1e-03)
+    >>> find_equilibrium_position(my_fluid_flow)
+    >>> (my_fluid_flow.xi, my_fluid_flow.yi) # doctest: +ELLIPSIS
+    (2.25...
     """
-    fluid_flow_object.calculate_coefficients()
+
+    def residuals(x, *args):
+        """Calculates x component of the forces of the oil film and the
+        difference between the y component and the load.
+        Parameters
+        ----------
+        x: array
+            Rotor center coordinates
+        *args : dict
+            Dictionary instantiating the ross.FluidFlow class.
+        Returns
+        -------
+        array
+            Array with the x component of the forces of the oil film and the difference
+            between the y component and the load.
+        """
+        bearing = args[0]
+        move_rotor_center_abs(
+            bearing,
+            x[0] * fluid_flow_object.radial_clearance,
+            x[1] * fluid_flow_object.radial_clearance,
+        )
+        bearing.geometry_description()
+        bearing.calculate_pressure_matrix_numerical()
+        (_, _, fx, fy,) = calculate_oil_film_force(bearing, force_type="numerical")
+        return np.array([fx, (bearing.load - fy)])
+
+    if fluid_flow_object.load is None:
+        sys.exit("Load must be given to calculate the equilibrium position.")
+    x0 = np.array([0, -1e-3 * fluid_flow_object.radial_clearance])
+    move_rotor_center_abs(fluid_flow_object, x0[0], x0[1])
+    fluid_flow_object.geometry_description()
     fluid_flow_object.calculate_pressure_matrix_numerical()
-    r_force, t_force, force_x, force_y = calculate_oil_film_force(
+    (_, _, fx, fy,) = calculate_oil_film_force(
         fluid_flow_object, force_type="numerical"
     )
-    increment = increment_factor * fluid_flow_object.eccentricity
-    error_x = abs(force_x)
-    error_y = abs(force_y - fluid_flow_object.load)
-    error = max(error_x, error_y)
-    k = 1
-    map_vector = []
-    while error > tolerance and k <= max_iterations:
-        increment_x = increment
-        increment_y = increment
-        iter_x = 0
-        iter_y = 0
-        previous_x = fluid_flow_object.xi
-        previous_y = fluid_flow_object.yi
-        infinite_loop_x_check = False
-        infinite_loop_y_check = False
-        if print_along:
-            print("\nIteration " + str(k) + "\n")
-        while error_x > tolerance:
-            iter_x += 1
-            move_rotor_center(fluid_flow_object, increment_x, 0)
-            fluid_flow_object.calculate_coefficients()
-            fluid_flow_object.calculate_pressure_matrix_numerical()
-            (
-                new_r_force,
-                new_t_force,
-                new_force_x,
-                new_force_y,
-            ) = calculate_oil_film_force(fluid_flow_object, force_type="numerical")
-            new_error_x = abs(new_force_x)
-            move_rotor_center(fluid_flow_object, -increment_x, 0)
-            if print_along:
-                print("Iteration in x axis " + str(iter_x))
-                print("Force x: " + str(new_force_x))
-                print("Previous force x: " + str(force_x))
-                print("Increment x: ", str(increment_x))
-                print("Error x: " + str(new_error_x))
-                print("Previous error x: " + str(error_x) + "\n")
-            if new_force_x * force_x < 0:
-                infinite_loop_x_check = False
-                increment_x = increment_x / 10
-                if print_along:
-                    print("Went beyond error 0. Reducing increment. \n")
-                if abs(increment_x) < abs(increment * increment_reduction_limit):
-                    if print_along:
-                        print("Increment too low. Breaking x iteration. \n")
-                    break
-            elif new_error_x > error_x:
-                if print_along:
-                    print("Error increased. Changing sign of increment. \n")
-                increment_x = -increment_x
-                if infinite_loop_x_check:
-                    break
-                else:
-                    infinite_loop_x_check = True
-            else:
-                infinite_loop_x_check = False
-                move_rotor_center(fluid_flow_object, increment_x, 0)
-                error_x = new_error_x
-                force_x = new_force_x
-                force_y = new_force_y
-                error_y = abs(new_force_y - fluid_flow_object.load)
-                error = max(error_x, error_y)
-
-        while error_y > tolerance:
-            iter_y += 1
-            move_rotor_center(fluid_flow_object, 0, increment_y)
-            fluid_flow_object.calculate_coefficients()
-            fluid_flow_object.calculate_pressure_matrix_numerical()
-            (
-                new_r_force,
-                new_t_force,
-                new_force_x,
-                new_force_y,
-            ) = calculate_oil_film_force(fluid_flow_object, force_type="numerical")
-            new_error_y = abs(new_force_y - fluid_flow_object.load)
-            move_rotor_center(fluid_flow_object, 0, -increment_y)
-            if print_along:
-                print("Iteration in y axis " + str(iter_y))
-                print("Force y: " + str(new_force_y))
-                print("Previous force y: " + str(force_y))
-                print("Increment y: ", str(increment_y))
-                print(
-                    "Force y minus load: " + str(new_force_y - fluid_flow_object.load)
-                )
-                print(
-                    "Previous force y minus load: "
-                    + str(force_y - fluid_flow_object.load)
-                )
-                print("Error y: " + str(new_error_y))
-                print("Previous error y: " + str(error_y) + "\n")
-            if (new_force_y - fluid_flow_object.load) * (
-                force_y - fluid_flow_object.load
-            ) < 0:
-                infinite_loop_y_check = False
-                increment_y = increment_y / 10
-                if print_along:
-                    print("Went beyond error 0. Reducing increment. \n")
-                if abs(increment_y) < abs(increment * increment_reduction_limit):
-                    if print_along:
-                        print("Increment too low. Breaking y iteration. \n")
-                    break
-            elif new_error_y > error_y:
-                if print_along:
-                    print("Error increased. Changing sign of increment. \n")
-                increment_y = -increment_y
-                if infinite_loop_y_check:
-                    break
-                else:
-                    infinite_loop_y_check = True
-            else:
-                infinite_loop_y_check = False
-                move_rotor_center(fluid_flow_object, 0, increment_y)
-                error_y = new_error_y
-                force_y = new_force_y
-                force_x = new_force_x
-                error_x = abs(new_force_x)
-                error = max(error_x, error_y)
-        if print_along:
-            print("Iteration " + str(k))
-            print("Error x: " + str(error_x))
-            print("Error y: " + str(error_y))
-            print(
-                "Current x, y: ("
-                + str(fluid_flow_object.xi)
-                + ", "
-                + str(fluid_flow_object.yi)
-                + ")"
-            )
-        k += 1
-        map_vector.append(
-            [fluid_flow_object.xi, fluid_flow_object.yi, error_x, error_y]
+    result = least_squares(
+        residuals,
+        x0,
+        args=[fluid_flow_object],
+        jac="3-point",
+        bounds=([-1, -1], [1, 1]),
+    )
+    move_rotor_center_abs(
+        fluid_flow_object,
+        result.x[0] * fluid_flow_object.radial_clearance,
+        result.x[1] * fluid_flow_object.radial_clearance,
+    )
+    fluid_flow_object.geometry_description()
+    if print_equilibrium_position is True:
+        print(
+            "The equilibrium position (x0, y0) is: (",
+            result.x[0] * fluid_flow_object.radial_clearance,
+            ",",
+            result.x[1] * fluid_flow_object.radial_clearance,
+            ")",
         )
-        if previous_x == fluid_flow_object.xi and previous_y == fluid_flow_object.yi:
-            if print_along:
-                print("Rotor center did not move during iteration. Breaking.")
-            break
-
-    if print_along:
-        print(map_vector)
-    if return_iteration_map:
-        return map_vector
diff --git a/ross/fluid_flow/fluid_flow_geometry.py b/ross/fluid_flow/fluid_flow_geometry.py
index 98a7db86..501f627b 100644
--- a/ross/fluid_flow/fluid_flow_geometry.py
+++ b/ross/fluid_flow/fluid_flow_geometry.py
@@ -1,9 +1,10 @@
 import numpy as np
-from scipy.optimize import root
+from scipy.optimize import least_squares, root
 
 
 def calculate_attitude_angle(eccentricity_ratio):
     """Calculates the attitude angle based on the eccentricity ratio.
+    Suitable only for short bearings.
     Parameters
     ----------
     eccentricity_ratio: float
@@ -334,7 +335,7 @@ def move_rotor_center(fluid_flow_object, dx, dy):
     >>> my_fluid_flow = fluid_flow_example2()
     >>> move_rotor_center(my_fluid_flow, 0, 0)
     >>> my_fluid_flow.eccentricity # doctest: +ELLIPSIS
-    2.6627...
+    2.54...
     """
     fluid_flow_object.xi = fluid_flow_object.xi + dx
     fluid_flow_object.yi = fluid_flow_object.yi + dy
@@ -347,3 +348,36 @@ def move_rotor_center(fluid_flow_object, dx, dy):
     fluid_flow_object.attitude_angle = np.arccos(
         abs(fluid_flow_object.yi / fluid_flow_object.eccentricity)
     )
+
+
+def move_rotor_center_abs(fluid_flow_object, x, y):
+    """Moves the rotor center to the coordinates (x, y) and calculates new eccentricity,
+    attitude angle, and rotor center.
+    Parameters
+    ----------
+    fluid_flow_object: A FluidFlow object.
+    x: float
+        Coordinate along x axis.
+    y: float
+        Coordinate along y axis.
+    Returns
+    -------
+    None
+    --------
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
+    >>> move_rotor_center_abs(my_fluid_flow, 0, -1e-3*my_fluid_flow.radial_clearance)
+    >>> my_fluid_flow.eccentricity # doctest: +ELLIPSIS
+    1.99...
+    """
+    fluid_flow_object.xi = x
+    fluid_flow_object.yi = y
+    fluid_flow_object.eccentricity = np.sqrt(
+        fluid_flow_object.xi ** 2 + fluid_flow_object.yi ** 2
+    )
+    fluid_flow_object.eccentricity_ratio = (
+        fluid_flow_object.eccentricity / fluid_flow_object.radial_clearance
+    )
+    fluid_flow_object.attitude_angle = np.arccos(
+        abs(fluid_flow_object.yi / fluid_flow_object.eccentricity)
+    )


Test patch:
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
index 68808ab8..28a341de 100644
--- a/ross/tests/test_fluid_flow.py
+++ b/ross/tests/test_fluid_flow.py
@@ -7,6 +7,7 @@ from numpy.testing import assert_allclose
 from plotly import graph_objects as go
 
 from ross.fluid_flow import fluid_flow as flow
+from ross.fluid_flow.fluid_flow import fluid_flow_example2
 from ross.fluid_flow.fluid_flow_coefficients import (
     calculate_oil_film_force, calculate_short_damping_matrix,
     calculate_short_stiffness_matrix,
@@ -75,6 +76,7 @@ def fluid_flow_short_friswell(set_load=True):
             rho,
             load=load,
             immediately_calculate_pressure_matrix_numerically=False,
+            bearing_type="short_bearing",
         )
     else:
         return flow.FluidFlow(
@@ -90,6 +92,7 @@ def fluid_flow_short_friswell(set_load=True):
             rho,
             eccentricity=eccentricity,
             immediately_calculate_pressure_matrix_numerically=False,
+            bearing_type="short_bearing",
         )
 
 
@@ -188,7 +191,6 @@ def fluid_flow_short_numerical():
     radius_stator = 0.1
     visc = 0.015
     rho = 860.0
-    attitude_angle = None
     eccentricity = 0.001
     return flow.FluidFlow(
         nz,
@@ -201,7 +203,6 @@ def fluid_flow_short_numerical():
         radius_stator,
         visc,
         rho,
-        attitude_angle=attitude_angle,
         eccentricity=eccentricity,
         immediately_calculate_pressure_matrix_numerically=False,
     )
@@ -220,6 +221,7 @@ def fluid_flow_long_numerical():
     visc = 0.015
     rho = 860.0
     eccentricity = 0.0001
+    attitude_angle = (5 * np.pi) / 100
     return flow.FluidFlow(
         nz,
         ntheta,
@@ -232,7 +234,8 @@ def fluid_flow_long_numerical():
         visc,
         rho,
         eccentricity=eccentricity,
-        immediately_calculate_pressure_matrix_numerically=False,
+        attitude_angle=attitude_angle,
+        immediately_calculate_pressure_matrix_numerically=True,
     )
 
 
@@ -314,20 +317,10 @@ def test_plots():
 
 
 def test_find_equilibrium_position():
-    bearing = flow.fluid_flow_example2()
-    find_equilibrium_position(
-        bearing,
-        print_along=False,
-        tolerance=0.1,
-        increment_factor=0.01,
-        max_iterations=2,
-        increment_reduction_limit=1e-03,
-    )
-    assert_allclose(
-        bearing.eccentricity,
-        (bearing.radius_stator - bearing.radius_rotor) * 0.2663,
-        atol=0.001,
-    )
+    bearing = fluid_flow_example2()
+    (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)
 
 
 def test_move_rotor_center():