Repository: ross-rotordynamics/ross
Pull request #349 (2019-10-17 16:28:42)

Improve oil film force calculations


Patch:
diff --git a/ross/bearing_seal_element.py b/ross/bearing_seal_element.py
index d8b14f82..7bf54963 100644
--- a/ross/bearing_seal_element.py
+++ b/ross/bearing_seal_element.py
@@ -895,7 +895,7 @@ class BearingElement(Element):
         BearingElement(n=0,
          kxx=[...
         """
-        fluid_flow = flow.PressureMatrix(
+        fluid_flow = flow.FluidFlow(
             nz,
             ntheta,
             nradius,
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index b216a240..34d38897 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -1,17 +1,12 @@
 import sys
 import numpy as np
-from ross.fluid_flow.fluid_flow_geometry import (
-    calculate_attitude_angle,
-    internal_radius_function,
-    external_radius_function,
-    modified_sommerfeld_number,
-    calculate_eccentricity_ratio,
-    calculate_rotor_load,
-)
+from ross.fluid_flow.fluid_flow_geometry import calculate_attitude_angle, internal_radius_function,\
+    external_radius_function, modified_sommerfeld_number, calculate_eccentricity_ratio, calculate_rotor_load
 
 
-class PressureMatrix:
-    r"""This class calculates the pressure matrix of the object with the given parameters.
+class FluidFlow:
+    r"""This class intends to calculate the pressure matrix and the stiffness and damping matrices
+    of a fluid flow with the given parameters.
 
     It is supposed to be an attribute of a bearing element,
     but can work on its on to provide graphics for the user.
@@ -52,8 +47,8 @@ class PressureMatrix:
     eccentricity: float, optional
         Eccentricity (m) is the euclidean distance between rotor and stator centers.
         The center of the stator is in position (0,0).
-    beta: float, optional
-        Attitude angle. Angle between the origin and the eccentricity (rad).
+    attitude_angle: float, optional
+        Attitude angle. Angle between the load line and the eccentricity (rad).
 
     Fluid characteristics
     ^^^^^^^^^^^^^^^^^^^^^
@@ -65,7 +60,7 @@ class PressureMatrix:
 
     Returns
     -------
-    An object containing the pressure matrix and its data.
+    An object containing the fluid flow and its data.
 
     Attributes
     ----------
@@ -118,8 +113,8 @@ class PressureMatrix:
     bearing_type: str
         type of structure. 'short_bearing': short; 'long_bearing': long;
         'medium_size': in between short and long.
-        if length/radius_stator <= 1/4 it is short.
-        if length/radius_stator > 8 it is long.
+        if length/diameter <= 1/4 it is short.
+        if length/diameter > 8 it is long.
     radial_clearance: float
         Difference between both stator and rotor radius, regardless of eccentricity.
     analytical_pressure_matrix_available: bool
@@ -145,11 +140,11 @@ class PressureMatrix:
     >>> viscosity = 0.015
     >>> density = 860.
     >>> eccentricity = 0.001
-    >>> beta = np.pi
-    >>> my_fluid_flow = flow.PressureMatrix(nz, ntheta, nradius, length,
+    >>> attitude_angle = np.pi
+    >>> my_fluid_flow = flow.FluidFlow(nz, ntheta, nradius, length,
     ...                                          omega, p_in, p_out, radius_rotor,
     ...                                          radius_stator, viscosity, density,
-    ...                                          beta=beta, eccentricity=eccentricity)
+    ...                                          attitude_angle=attitude_angle, eccentricity=eccentricity)
     >>> my_fluid_flow.calculate_pressure_matrix_analytical() # doctest: +ELLIPSIS
     array([[-0.00000...
     >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
@@ -158,25 +153,25 @@ class PressureMatrix:
     >>> # show(my_fluid_flow.plot_eccentricity())
     >>> # show(my_fluid_flow.plot_pressure_theta(z=int(nz/2)))
     >>> matplot_pressure_theta(my_fluid_flow, z=int(nz/2)) # doctest: +ELLIPSIS
-    <matplotlib.axes._subplots...
+    <matplotlib.axes...
     """
 
     def __init__(
-        self,
-        nz,
-        ntheta,
-        nradius,
-        length,
-        omega,
-        p_in,
-        p_out,
-        radius_rotor,
-        radius_stator,
-        viscosity,
-        density,
-        beta=None,
-        eccentricity=None,
-        load=None,
+            self,
+            nz,
+            ntheta,
+            nradius,
+            length,
+            omega,
+            p_in,
+            p_out,
+            radius_rotor,
+            radius_stator,
+            viscosity,
+            density,
+            attitude_angle=None,
+            eccentricity=None,
+            load=None,
     ):
         if load is None and eccentricity is None:
             sys.exit("Either load or eccentricity must be given.")
@@ -201,9 +196,9 @@ class PressureMatrix:
         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 / 4:
+        if self.length / (2 * self.radius_stator) <= 1 / 4:
             self.bearing_type = "short_bearing"
-        elif self.length / self.radius_stator >= 8:
+        elif self.length / (2 * self.radius_stator) >= 8:
             self.bearing_type = "long_bearing"
         else:
             self.bearing_type = "medium_size"
@@ -212,33 +207,21 @@ class PressureMatrix:
         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.radius_stator, self.omega, self.viscosity, self.length,
+                       self.load, self.radial_clearance)
             self.eccentricity = (
-                calculate_eccentricity_ratio(modified_s)
-                * self.difference_between_radius
+                    calculate_eccentricity_ratio(modified_s) * self.difference_between_radius
             )
         self.eccentricity_ratio = self.eccentricity / self.difference_between_radius
         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 beta is None:
-            self.beta = calculate_attitude_angle(self.eccentricity_ratio)
+            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.beta = beta
-        self.xi = self.eccentricity * np.cos(2 * np.pi - self.beta)
-        self.yi = self.eccentricity * np.sin(2 * np.pi - self.beta)
+            self.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.re = np.zeros([self.nz, self.ntheta])
         self.ri = np.zeros([self.nz, self.ntheta])
         self.z_list = np.zeros(self.nz)
@@ -279,11 +262,11 @@ class PressureMatrix:
             Pressure matrix of size (nz x ntheta)
         Examples
         --------
-        >>> my_fluid_flow = pressure_matrix_example()
+        >>> my_fluid_flow = fluid_flow_example()
         >>> my_fluid_flow.calculate_pressure_matrix_analytical() # doctest: +ELLIPSIS
         array([[...
         """
-        if self.bearing_type == "short_bearing" or force_type == "short":
+        if self.bearing_type == "short_bearing" or force_type == 'short':
             if method == 0:
                 for i in range(0, self.nz):
                     for j in range(0, self.ntheta):
@@ -310,41 +293,24 @@ class PressureMatrix:
                         # fmt: on
                         if self.p_mat_analytical[i][j] < 0:
                             self.p_mat_analytical[i][j] = 0
-        elif self.bearing_type == "long_bearing" or force_type == "long":
+        elif self.bearing_type == "long_bearing" or force_type == 'long':
             if method == 0:
                 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.ri[i][j] / self.difference_between_radius) ** 2
-                                * self.eccentricity_ratio
-                                * np.sin(self.gama[i][j])
-                                * (
-                                    2
-                                    + self.eccentricity_ratio * np.cos(self.gama[i][j])
-                                )
-                            )
-                            / (
-                                (2 + self.eccentricity_ratio ** 2)
-                                * (
-                                    1
-                                    + self.eccentricity_ratio * np.cos(self.gama[i][j])
-                                )
-                                ** 2
-                            )
-                            + self.p_in
-                        )
+                        self.p_mat_analytical[i][j] = (6 * self.viscosity * self.omega *
+                                                       (self.ri[i][j] / self.difference_between_radius) ** 2 *
+                                                       self.eccentricity_ratio * np.sin(self.gama[i][j]) *
+                                                       (2 + self.eccentricity_ratio * np.cos(self.gama[i][j]))) / (
+                                                              (2 + self.eccentricity_ratio ** 2) *
+                                                              (1 + self.eccentricity_ratio * np.cos(
+                                                                  self.gama[i][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":
-            raise ValueError(
-                "The pressure matrix for a bearing that is neither short or long can only be calculated "
-                "numerically. Try calling calculate_pressure_matrix_numerical or setting force_type "
-                "to either 'short' or 'long' in calculate_pressure_matrix_analytical."
-            )
+            raise ValueError("The pressure matrix for a bearing that is neither short or long can only be calculated "
+                             "numerically. Try calling calculate_pressure_matrix_numerical or setting force_type "
+                             "to either 'short' or 'long' in calculate_pressure_matrix_analytical.")
         self.analytical_pressure_matrix_available = True
         return self.p_mat_analytical
 
@@ -354,7 +320,7 @@ class PressureMatrix:
         derivatives). It is executed when the class is instantiated.
         Examples
         --------
-        >>> my_fluid_flow = pressure_matrix_example()
+        >>> my_fluid_flow = fluid_flow_example()
         >>> my_fluid_flow.calculate_coefficients()
         >>> my_fluid_flow.c0w # doctest: +ELLIPSIS
         array([[...
@@ -365,11 +331,11 @@ class PressureMatrix:
             eccentricity_error = False
             for j in range(0, self.ntheta):
                 # fmt: off
-                self.gama[i][j] = j * self.dtheta + (np.pi - self.beta)
+                self.gama[i][j] = j * self.dtheta + np.pi / 2 + self.attitude_angle
                 [radius_external, self.xre[i][j], self.yre[i][j]] = \
                     external_radius_function(self.gama[i][j], self.radius_stator)
                 [radius_internal, self.xri[i][j], self.yri[i][j]] = \
-                    internal_radius_function(self.gama[i][j], self.beta, self.radius_rotor,
+                    internal_radius_function(self.gama[i][j], self.attitude_angle, self.radius_rotor,
                                              self.eccentricity)
                 self.re[i][j] = radius_external
                 self.ri[i][j] = radius_internal
@@ -402,7 +368,7 @@ class PressureMatrix:
                 # fmt: on
                 if not eccentricity_error:
                     if abs(self.xri[i][j]) > abs(self.xre[i][j]) or abs(
-                        self.yri[i][j]
+                            self.yri[i][j]
                     ) > abs(self.yre[i][j]):
                         eccentricity_error = True
             if eccentricity_error:
@@ -415,7 +381,7 @@ class PressureMatrix:
         """This function assembles the matrix M and the independent vector f.
         Examples
         --------
-        >>> my_fluid_flow = pressure_matrix_example()
+        >>> my_fluid_flow = fluid_flow_example()
         >>> my_fluid_flow.mounting_matrix()
         >>> my_fluid_flow.M # doctest: +ELLIPSIS
         array([[...
@@ -480,7 +446,7 @@ class PressureMatrix:
         """This function resolves the linear system [M]{P}={f}.
         Examples
         --------
-        >>> my_fluid_flow = pressure_matrix_example()
+        >>> my_fluid_flow = fluid_flow_example()
         >>> my_fluid_flow.mounting_matrix()
         >>> my_fluid_flow.resolves_matrix()
         >>> my_fluid_flow.P # doctest: +ELLIPSIS
@@ -496,7 +462,7 @@ class PressureMatrix:
             Pressure matrix of size (nz x ntheta)
         Examples
         --------
-        >>> my_fluid_flow = pressure_matrix_example()
+        >>> my_fluid_flow = fluid_flow_example()
         >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
         array([[...
         """
@@ -513,8 +479,8 @@ class PressureMatrix:
         return self.p_mat_numerical
 
 
-def pressure_matrix_example():
-    """This function returns an instance of a simple pressure matrix.
+def fluid_flow_example():
+    """This function returns an instance of a simple fluid flow.
     The purpose is to make available a simple model
     so that doctest can be written using it.
 
@@ -523,27 +489,20 @@ def pressure_matrix_example():
 
     Returns
     -------
-    An instance of a pressure matrix object.
+    An instance of a fluid flow object.
 
     Examples
     --------
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> my_fluid_flow = fluid_flow_example()
     >>> my_fluid_flow.eccentricity
     0.001
     """
-    my_pressure_matrix = PressureMatrix(
-        nz=8,
-        ntheta=64,
-        nradius=11,
-        length=0.01,
-        omega=100.0 * 2 * np.pi / 60,
-        p_in=0.0,
-        p_out=0.0,
-        radius_rotor=0.08,
-        radius_stator=0.1,
-        viscosity=0.015,
-        density=860.0,
-        eccentricity=0.001,
-        beta=np.pi,
-    )
+    my_pressure_matrix = FluidFlow(nz=8, ntheta=64, nradius=11, length=0.01, omega=100. * 2 * np.pi / 60,
+                                   p_in=0., p_out=0., radius_rotor=0.08, radius_stator=0.1,
+                                   viscosity=0.015, density=860., eccentricity=0.001, attitude_angle=np.pi)
     return my_pressure_matrix
+
+
+
+
+
diff --git a/ross/fluid_flow/fluid_flow_coefficients.py b/ross/fluid_flow/fluid_flow_coefficients.py
index a7512012..a138c53e 100644
--- a/ross/fluid_flow/fluid_flow_coefficients.py
+++ b/ross/fluid_flow/fluid_flow_coefficients.py
@@ -1,5 +1,6 @@
 import numpy as np
 from scipy import integrate
+from math import isnan
 
 
 def calculate_analytical_stiffness_matrix(load, eccentricity_ratio, radial_clearance):
@@ -20,8 +21,8 @@ def calculate_analytical_stiffness_matrix(load, eccentricity_ratio, radial_clear
         A list of length four including stiffness floats in this order: kxx, kxy, kyx, kyy
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> load = my_fluid_flow.load
     >>> eccentricity_ratio = my_fluid_flow.eccentricity_ratio
     >>> radial_clearance = my_fluid_flow.radial_clearance
@@ -66,8 +67,8 @@ def calculate_analytical_damping_matrix(load, eccentricity_ratio, radial_clearan
         A list of length four including damping floats in this order: cxx, cxy, cyx, cyy
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> load = my_fluid_flow.load
     >>> eccentricity_ratio = my_fluid_flow.eccentricity_ratio
     >>> radial_clearance = my_fluid_flow.radial_clearance
@@ -92,81 +93,89 @@ def calculate_analytical_damping_matrix(load, eccentricity_ratio, radial_clearan
     return [cxx, cxy, cyx, cyy]
 
 
-def calculate_oil_film_force(pressure_matrix_object, force_type=None):
+def calculate_oil_film_force(fluid_flow_object, force_type=None):
     """This function calculates the forces of the oil film in the N and T directions, ie in the
     opposite direction to the eccentricity and in the tangential direction.
     Parameters
     ----------
-    pressure_matrix_object: A PressureMatrix object.
+    fluid_flow_object: A FluidFlow object.
     force_type: str
         If set, calculates the oil film force matrix analytically considering the chosen type: 'short' or 'long'.
         If set to 'numerical', calculates the oil film force numerically.
     Returns
     -------
-    normal_force: float
+    radial_force: float
         Force of the oil film in the opposite direction to the eccentricity direction.
     tangential_force: float
         Force of the oil film in the tangential direction
+    f_x: float
+        Components of forces in the x direction
+    f_y: float
+        Components of forces in the y direction
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> calculate_oil_film_force(my_fluid_flow) # doctest: +ELLIPSIS
     (...
     """
-    if force_type != 'numerical' and (force_type == 'short' or pressure_matrix_object.bearing_type == 'short_bearing'):
-        normal_force = 0.5 * \
-           pressure_matrix_object.viscosity * \
-           (pressure_matrix_object.radius_rotor / pressure_matrix_object.difference_between_radius) ** 2 * \
-           (pressure_matrix_object.length ** 3 / pressure_matrix_object.radius_rotor) * \
-           ((2 * pressure_matrix_object.eccentricity_ratio ** 2 * pressure_matrix_object.omega) /
-            (1 - pressure_matrix_object.eccentricity_ratio ** 2) ** 2)
-
-        tangential_force = 0.5 * pressure_matrix_object.viscosity * \
-            (pressure_matrix_object.radius_rotor / pressure_matrix_object.difference_between_radius) ** 2 * \
-            (pressure_matrix_object.length ** 3 / pressure_matrix_object.radius_rotor) * (
-             (np.pi * pressure_matrix_object.eccentricity_ratio * pressure_matrix_object.omega) /
-             (2 * (1 - pressure_matrix_object.eccentricity_ratio ** 2) ** (3. / 2)))
-    elif force_type != 'numerical' and (force_type == 'long' or pressure_matrix_object.bearing_type == 'long_bearing'):
-        normal_force = 6 * pressure_matrix_object.viscosity * \
-            (pressure_matrix_object.radius_rotor / pressure_matrix_object.difference_between_radius) ** 2 * \
-            pressure_matrix_object.radius_rotor * pressure_matrix_object.length * \
-            ((2 * pressure_matrix_object.eccentricity_ratio ** 2 * pressure_matrix_object.omega) /
-             ((2 + pressure_matrix_object.eccentricity_ratio ** 2) *
-              (1 - pressure_matrix_object.eccentricity_ratio ** 2)))
-        tangential_force = 6 * pressure_matrix_object.viscosity * (pressure_matrix_object.radius_rotor / pressure_matrix_object.difference_between_radius) ** 2 * \
-            pressure_matrix_object.radius_rotor * pressure_matrix_object.length * \
-            ((np.pi * pressure_matrix_object.eccentricity_ratio * pressure_matrix_object.omega) /
-             ((2 + pressure_matrix_object.eccentricity_ratio ** 2) *
-             (1 - pressure_matrix_object.eccentricity_ratio ** 2) ** 0.5))
+    if force_type != 'numerical' and (force_type == 'short' or fluid_flow_object.bearing_type == 'short_bearing'):
+        radial_force = 0.5 * \
+           fluid_flow_object.viscosity * \
+           (fluid_flow_object.radius_rotor / fluid_flow_object.difference_between_radius) ** 2 * \
+           (fluid_flow_object.length ** 3 / fluid_flow_object.radius_rotor) * \
+           ((2 * fluid_flow_object.eccentricity_ratio ** 2 * fluid_flow_object.omega) /
+            (1 - fluid_flow_object.eccentricity_ratio ** 2) ** 2)
+
+        tangential_force = 0.5 * fluid_flow_object.viscosity * \
+            (fluid_flow_object.radius_rotor / fluid_flow_object.difference_between_radius) ** 2 * \
+            (fluid_flow_object.length ** 3 / fluid_flow_object.radius_rotor) * (
+             (np.pi * fluid_flow_object.eccentricity_ratio * fluid_flow_object.omega) /
+             (2 * (1 - fluid_flow_object.eccentricity_ratio ** 2) ** (3. / 2)))
+    elif force_type != 'numerical' and (force_type == 'long' or fluid_flow_object.bearing_type == 'long_bearing'):
+        radial_force = 6 * fluid_flow_object.viscosity * \
+            (fluid_flow_object.radius_rotor / fluid_flow_object.difference_between_radius) ** 2 * \
+            fluid_flow_object.radius_rotor * fluid_flow_object.length * \
+            ((2 * fluid_flow_object.eccentricity_ratio ** 2 * fluid_flow_object.omega) /
+             ((2 + fluid_flow_object.eccentricity_ratio ** 2) *
+              (1 - fluid_flow_object.eccentricity_ratio ** 2)))
+        tangential_force = 6 * fluid_flow_object.viscosity * (fluid_flow_object.radius_rotor / fluid_flow_object.difference_between_radius) ** 2 * \
+            fluid_flow_object.radius_rotor * fluid_flow_object.length * \
+            ((np.pi * fluid_flow_object.eccentricity_ratio * fluid_flow_object.omega) /
+             ((2 + fluid_flow_object.eccentricity_ratio ** 2) *
+             (1 - fluid_flow_object.eccentricity_ratio ** 2) ** 0.5))
     else:
-        p_mat = pressure_matrix_object.p_mat_numerical
-        a = np.zeros([pressure_matrix_object.nz, pressure_matrix_object.ntheta])
-        b = np.zeros([pressure_matrix_object.nz, pressure_matrix_object.ntheta])
-        g1 = np.zeros(pressure_matrix_object.nz)
-        g2 = np.zeros(pressure_matrix_object.nz)
-        for i in range(pressure_matrix_object.nz):
-            for j in range(pressure_matrix_object.ntheta):
-                a[i][j] = p_mat[i][j] * np.cos(pressure_matrix_object.gama[i][j])
-                b[i][j] = p_mat[i][j] * np.sin(pressure_matrix_object.gama[i][j])
-
-        for i in range(pressure_matrix_object.nz):
-            g1[i] = integrate.simps(a[i][:], pressure_matrix_object.gama[0])
-            g2[i] = integrate.simps(b[i][:], pressure_matrix_object.gama[0])
-
-        integral1 = integrate.simps(g1, pressure_matrix_object.z_list)
-        integral2 = integrate.simps(g2, pressure_matrix_object.z_list)
-
-        normal_force = - pressure_matrix_object.radius_rotor * integral1
-        tangential_force = pressure_matrix_object.radius_rotor * integral2
-    return normal_force, tangential_force
-
-
-
-
-
-
-
-
-
-
+        p_mat = fluid_flow_object.p_mat_numerical
+        a = np.zeros([fluid_flow_object.nz, fluid_flow_object.ntheta])
+        b = np.zeros([fluid_flow_object.nz, fluid_flow_object.ntheta])
+        g1 = np.zeros(fluid_flow_object.nz)
+        g2 = np.zeros(fluid_flow_object.nz)
+        base_vector = np.array([fluid_flow_object.xre[0][0] - fluid_flow_object.xi,
+                                fluid_flow_object.yre[0][0] - fluid_flow_object.yi])
+        for i in range(fluid_flow_object.nz):
+            for j in range(int(fluid_flow_object.ntheta/2)):
+                vector_from_rotor = np.array([fluid_flow_object.xre[i][j] - fluid_flow_object.xi,
+                                              fluid_flow_object.yre[i][j] - fluid_flow_object.yi])
+                angle_between_vectors = np.arccos(np.dot(base_vector, vector_from_rotor) /
+                                                  (np.linalg.norm(base_vector) * np.linalg.norm(vector_from_rotor)))
+                if isnan(angle_between_vectors):
+                    angle_between_vectors = 0
+                if angle_between_vectors != 0 and j*fluid_flow_object.dtheta > np.pi:
+                    angle_between_vectors += np.pi
+                a[i][j] = p_mat[i][j] * np.cos(angle_between_vectors)
+                b[i][j] = p_mat[i][j] * np.sin(angle_between_vectors)
+
+        for i in range(fluid_flow_object.nz):
+            g1[i] = integrate.simps(a[i][:], fluid_flow_object.gama[0])
+            g2[i] = integrate.simps(b[i][:], fluid_flow_object.gama[0])
+
+        integral1 = integrate.simps(g1, fluid_flow_object.z_list)
+        integral2 = integrate.simps(g2, fluid_flow_object.z_list)
+
+        radial_force = - fluid_flow_object.radius_rotor * integral1
+        tangential_force = fluid_flow_object.radius_rotor * integral2
+    force_x = - radial_force * np.sin(fluid_flow_object.attitude_angle) \
+        + tangential_force * np.cos(fluid_flow_object.attitude_angle)
+    force_y = radial_force * np.cos(fluid_flow_object.attitude_angle) \
+        + tangential_force * np.sin(fluid_flow_object.attitude_angle)
+    return radial_force, tangential_force, force_x, force_y
diff --git a/ross/fluid_flow/fluid_flow_geometry.py b/ross/fluid_flow/fluid_flow_geometry.py
index 0eaf7707..0aa3ad65 100644
--- a/ross/fluid_flow/fluid_flow_geometry.py
+++ b/ross/fluid_flow/fluid_flow_geometry.py
@@ -14,16 +14,18 @@ def calculate_attitude_angle(eccentricity_ratio):
         Attitude angle
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> calculate_attitude_angle(my_fluid_flow.eccentricity_ratio) # doctest: +ELLIPSIS
     1.5...
     """
-    return np.arctan(np.pi * (1 - eccentricity_ratio ** 2) /
-                     (4 * eccentricity_ratio))
+    return np.arctan(
+                    (np.pi * (1 - eccentricity_ratio ** 2)**(1/2)) /
+                    (4 * eccentricity_ratio)
+                     )
 
 
-def internal_radius_function(gama, beta, radius_rotor, eccentricity):
+def internal_radius_function(gama, attitude_angle, radius_rotor, eccentricity):
     """This function calculates the x and y of the internal radius of the rotor,
     as well as its distance from the origin, given the distance in the theta-axis,
     the attitude angle, the radius of the rotor and the eccentricity.
@@ -31,7 +33,7 @@ def internal_radius_function(gama, beta, radius_rotor, eccentricity):
     ----------
     gama: float
         Gama is the distance in the theta-axis. It should range from 0 to 2*np.pi.
-    beta: float
+    attitude_angle: float
         Attitude angle. Angle between the origin and the eccentricity (rad).
     radius_rotor: float
         The radius of the journal.
@@ -47,19 +49,19 @@ def internal_radius_function(gama, beta, radius_rotor, eccentricity):
         The position y of the returned internal radius.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
-    >>> beta = my_fluid_flow.beta
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
+    >>> attitude_angle = my_fluid_flow.attitude_angle
     >>> radius_rotor = my_fluid_flow.radius_rotor
     >>> eccentricity = my_fluid_flow.eccentricity
-    >>> radius_internal, xri, yri = internal_radius_function(0, beta, radius_rotor, eccentricity)
-    >>> radius_internal
-    0.079
+    >>> radius_internal, xri, yri = internal_radius_function(0, attitude_angle, radius_rotor, eccentricity)
+    >>> radius_internal # doctest: +ELLIPSIS
+    0.079...
     """
-    if (np.pi - beta) < gama < (2 * np.pi - beta):
-        alpha = np.absolute(2 * np.pi - gama - beta)
+    if (np.pi / 2 + attitude_angle) < gama < (3 * np.pi / 2 + attitude_angle):
+        alpha = np.absolute(3 * np.pi / 2 - gama + attitude_angle)
     else:
-        alpha = beta + gama
+        alpha = gama + np.pi / 2 - attitude_angle
     radius_internal = np.sqrt(radius_rotor ** 2 - (eccentricity * np.sin(alpha)) ** 2) + eccentricity * np.cos(alpha)
     xri = radius_internal * np.cos(gama)
     yri = radius_internal * np.sin(gama)
@@ -86,8 +88,8 @@ def external_radius_function(gama, radius_stator):
         The position y of the returned external radius.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> radius_stator = my_fluid_flow.radius_stator
     >>> radius_external, xre, yre = external_radius_function(0, radius_stator)
     >>> radius_external
@@ -123,8 +125,8 @@ def modified_sommerfeld_number(radius_stator, omega, viscosity, length, load, ra
         The modified sommerfeld number.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> radius_stator = my_fluid_flow.radius_stator
     >>> omega = my_fluid_flow.omega
     >>> viscosity = my_fluid_flow.viscosity
@@ -156,8 +158,8 @@ def sommerfeld_number(modified_s, radius_stator, length):
         The sommerfeld number.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> radius_stator = my_fluid_flow.radius_stator
     >>> omega = my_fluid_flow.omega
     >>> viscosity = my_fluid_flow.viscosity
@@ -227,8 +229,8 @@ def calculate_rotor_load(radius_stator, omega, viscosity, length, radial_clearan
         Load applied to the rotor.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> radius_stator = my_fluid_flow.radius_stator
     >>> omega = my_fluid_flow.omega
     >>> viscosity = my_fluid_flow.viscosity
diff --git a/ross/fluid_flow/fluid_flow_graphics.py b/ross/fluid_flow/fluid_flow_graphics.py
index e77d6548..76fe0d6b 100644
--- a/ross/fluid_flow/fluid_flow_graphics.py
+++ b/ross/fluid_flow/fluid_flow_graphics.py
@@ -1,14 +1,14 @@
 import matplotlib.pyplot as plt
-import numpy as np
 from bokeh.plotting import figure
+import numpy as np
 
 
-def plot_eccentricity(pressure_matrix_object, z=0):
+def plot_eccentricity(fluid_flow_object, 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.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     z: int, optional
         The distance in z where to cut and plot.
     Returns
@@ -17,39 +17,33 @@ def plot_eccentricity(pressure_matrix_object, z=0):
         An object containing the plot.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> fig = plot_eccentricity(my_fluid_flow, z=int(my_fluid_flow.nz/2))
     >>> # to show the plots you can use:
     >>> # show(fig)
     """
     p = figure(
-        title="Cut in plane Z=" + str(z), x_axis_label="X axis", y_axis_label="Y axis"
+        title="Cut in plane Z=" + str(z),
+        x_axis_label="X axis",
+        y_axis_label="Y axis",
     )
-    for j in range(0, pressure_matrix_object.ntheta):
-        p.circle(
-            pressure_matrix_object.xre[z][j],
-            pressure_matrix_object.yre[z][j],
-            color="red",
-        )
-        p.circle(
-            pressure_matrix_object.xri[z][j],
-            pressure_matrix_object.yri[z][j],
-            color="blue",
-        )
+    for j in range(0, fluid_flow_object.ntheta):
+        p.circle(fluid_flow_object.xre[z][j], fluid_flow_object.yre[z][j], color="red")
+        p.circle(fluid_flow_object.xri[z][j], fluid_flow_object.yri[z][j], color="blue")
         p.circle(0, 0, color="blue")
-        p.circle(pressure_matrix_object.xi, pressure_matrix_object.yi, color="red")
+        p.circle(fluid_flow_object.xi, fluid_flow_object.yi, color="red")
     p.circle(0, 0, color="black")
     return p
 
 
-def plot_pressure_z(pressure_matrix_object, theta=0):
+def plot_pressure_z(fluid_flow_object, theta=0):
     """This function assembles pressure graphic along the z-axis for one or both the
     numerically (blue) and analytically (red) calculated pressure matrices, depending on if
     one or both were calculated.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     theta: int, optional
         The theta to be considered.
     Returns
@@ -58,8 +52,8 @@ def plot_pressure_z(pressure_matrix_object, theta=0):
         An object containing the plot.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
     array([[...
     >>> fig = plot_pressure_z(my_fluid_flow, theta=int(my_fluid_flow.ntheta/2))
@@ -67,8 +61,8 @@ def plot_pressure_z(pressure_matrix_object, theta=0):
     >>> # show(fig)
     """
     if (
-        not pressure_matrix_object.numerical_pressure_matrix_available
-        and not pressure_matrix_object.analytical_pressure_matrix_available
+            not fluid_flow_object.numerical_pressure_matrix_available
+            and not fluid_flow_object.analytical_pressure_matrix_available
     ):
         raise ValueError(
             "Must calculate the pressure matrix. "
@@ -76,37 +70,26 @@ def plot_pressure_z(pressure_matrix_object, theta=0):
         )
     y_n = []
     y_a = []
-    for i in range(0, pressure_matrix_object.nz):
-        y_n.append(pressure_matrix_object.p_mat_numerical[i][theta])
-        y_a.append(pressure_matrix_object.p_mat_analytical[i][theta])
+    for i in range(0, fluid_flow_object.nz):
+        y_n.append(fluid_flow_object.p_mat_numerical[i][theta])
+        y_a.append(fluid_flow_object.p_mat_analytical[i][theta])
     p = figure(
-        title="Pressure along the Z direction (direction of flow); Theta=" + str(theta),
+        title="Pressure along the Z direction (direction of flow); Theta="
+              + str(theta),
         x_axis_label="Points along Z",
     )
-    if pressure_matrix_object.numerical_pressure_matrix_available:
-        p.line(
-            pressure_matrix_object.z_list,
-            y_n,
-            legend="Numerical pressure",
-            color="blue",
-            line_width=2,
-        )
-    if pressure_matrix_object.analytical_pressure_matrix_available:
-        p.line(
-            pressure_matrix_object.z_list,
-            y_a,
-            legend="Analytical pressure",
-            color="red",
-            line_width=2,
-        )
+    if fluid_flow_object.numerical_pressure_matrix_available:
+        p.line(fluid_flow_object.z_list, y_n, legend="Numerical pressure", color="blue", line_width=2)
+    if fluid_flow_object.analytical_pressure_matrix_available:
+        p.line(fluid_flow_object.z_list, y_a, legend="Analytical pressure", color="red", line_width=2)
     return p
 
 
-def plot_shape(pressure_matrix_object, theta=0):
+def plot_shape(fluid_flow_object, theta=0):
     """This function assembles a graphic representing the geometry of the rotor.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     theta: int, optional
         The theta to be considered.
     Returns
@@ -115,34 +98,34 @@ def plot_shape(pressure_matrix_object, theta=0):
         An object containing the plot.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> fig = plot_shape(my_fluid_flow, theta=int(my_fluid_flow.ntheta/2))
     >>> # to show the plots you can use:
     >>> # show(fig)
     """
-    y_re = np.zeros(pressure_matrix_object.nz)
-    y_ri = np.zeros(pressure_matrix_object.nz)
-    for i in range(0, pressure_matrix_object.nz):
-        y_re[i] = pressure_matrix_object.re[i][theta]
-        y_ri[i] = pressure_matrix_object.ri[i][theta]
+    y_re = np.zeros(fluid_flow_object.nz)
+    y_ri = np.zeros(fluid_flow_object.nz)
+    for i in range(0, fluid_flow_object.nz):
+        y_re[i] = fluid_flow_object.re[i][theta]
+        y_ri[i] = fluid_flow_object.ri[i][theta]
     p = figure(
         title="Shapes of stator and rotor along Z; Theta=" + str(theta),
         x_axis_label="Points along Z",
         y_axis_label="Radial direction",
     )
-    p.line(pressure_matrix_object.z_list, y_re, line_width=2, color="red")
-    p.line(pressure_matrix_object.z_list, y_ri, line_width=2, color="blue")
+    p.line(fluid_flow_object.z_list, y_re, line_width=2, color="red")
+    p.line(fluid_flow_object.z_list, y_ri, line_width=2, color="blue")
     return p
 
 
-def plot_pressure_theta(pressure_matrix_object, z=0):
+def plot_pressure_theta(fluid_flow_object, z=0):
     """This function assembles pressure graphic in the theta direction for a given z
     for one or both the numerically (blue) and analytically (red) calculated pressure matrices,
     depending on if one or both were calculated.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     z: int, optional
         The distance along z-axis to be considered.
     Returns
@@ -151,8 +134,8 @@ def plot_pressure_theta(pressure_matrix_object, z=0):
         An object containing the plot.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
     array([[...
     >>> fig = plot_pressure_theta(my_fluid_flow, z=int(my_fluid_flow.nz/2))
@@ -160,8 +143,8 @@ def plot_pressure_theta(pressure_matrix_object, z=0):
     >>> # show(fig)
     """
     if (
-        not pressure_matrix_object.numerical_pressure_matrix_available
-        and not pressure_matrix_object.analytical_pressure_matrix_available
+            not fluid_flow_object.numerical_pressure_matrix_available
+            and not fluid_flow_object.analytical_pressure_matrix_available
     ):
         raise ValueError(
             "Must calculate the pressure matrix. "
@@ -172,18 +155,18 @@ def plot_pressure_theta(pressure_matrix_object, z=0):
         x_axis_label="Points along Theta",
         y_axis_label="Pressure",
     )
-    if pressure_matrix_object.numerical_pressure_matrix_available:
+    if fluid_flow_object.numerical_pressure_matrix_available:
         p.line(
-            pressure_matrix_object.gama[z],
-            pressure_matrix_object.p_mat_numerical[z],
+            fluid_flow_object.gama[z],
+            fluid_flow_object.p_mat_numerical[z],
             legend="Numerical pressure",
             line_width=2,
             color="blue",
         )
-    elif pressure_matrix_object.analytical_pressure_matrix_available:
+    elif fluid_flow_object.analytical_pressure_matrix_available:
         p.line(
-            pressure_matrix_object.gama[z],
-            pressure_matrix_object.p_mat_analytical[z],
+            fluid_flow_object.gama[z],
+            fluid_flow_object.p_mat_analytical[z],
             legend="Analytical pressure",
             line_width=2,
             color="red",
@@ -191,12 +174,12 @@ def plot_pressure_theta(pressure_matrix_object, z=0):
     return p
 
 
-def matplot_eccentricity(pressure_matrix_object, z=0, ax=None):
+def matplot_eccentricity(fluid_flow_object, z=0, ax=None):
     """This function assembles pressure graphic along the z-axis using matplotlib.
     The first few plots are of a different color to indicate where theta begins.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     z: int, optional
         The distance in z where to cut and plot.
     ax : matplotlib axes, optional
@@ -207,9 +190,8 @@ def matplot_eccentricity(pressure_matrix_object, z=0, ax=None):
         Returns the axes object with the plot.
     Examples
     --------
-    >>> import matplotlib.pyplot as plt
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> fig, ax = plt.subplots()
     >>> ax = matplot_eccentricity(my_fluid_flow, z=int(my_fluid_flow.nz/2), ax=ax)
     >>> # to show the plots you can use:
@@ -217,15 +199,19 @@ def matplot_eccentricity(pressure_matrix_object, z=0, ax=None):
     """
     if ax is None:
         ax = plt.gca()
-    for j in range(0, pressure_matrix_object.ntheta):
-        ax.plot(
-            pressure_matrix_object.xre[z][j], pressure_matrix_object.yre[z][j], "r."
-        )
-        ax.plot(
-            pressure_matrix_object.xri[z][j], pressure_matrix_object.yri[z][j], "b."
-        )
-        ax.plot(0, 0, "r*")
-        ax.plot(pressure_matrix_object.xi, pressure_matrix_object.yi, "b*")
+    x_r = []
+    x_b = []
+    y_r = []
+    y_b = []
+    for j in range(0, fluid_flow_object.ntheta):
+        x_r.append(fluid_flow_object.xre[z][j])
+        y_r.append(fluid_flow_object.yre[z][j])
+        x_b.append(fluid_flow_object.xri[z][j])
+        y_b.append(fluid_flow_object.yri[z][j])
+    ax.plot(x_r, y_r, "r.")
+    ax.plot(x_b, y_b, "b.")
+    ax.plot(0, 0, "r*")
+    ax.plot(fluid_flow_object.xi, fluid_flow_object.yi, "b*")
     ax.set_title("Cut in plane Z=" + str(z))
     ax.set_xlabel("X axis")
     ax.set_ylabel("Y axis")
@@ -233,13 +219,13 @@ def matplot_eccentricity(pressure_matrix_object, z=0, ax=None):
     return ax
 
 
-def matplot_pressure_z(pressure_matrix_object, theta=0, ax=None):
+def matplot_pressure_z(fluid_flow_object, theta=0, ax=None):
     """This function assembles pressure graphic along the z-axis using matplotlib
     for one or both the numerically (blue) and analytically (red) calculated pressure matrices,
     depending on if one or both were calculated.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     theta: int, optional
         The distance in theta where to cut and plot.
     ax : matplotlib axes, optional
@@ -250,8 +236,8 @@ def matplot_pressure_z(pressure_matrix_object, theta=0, ax=None):
         Returns the axes object with the plot.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
     array([[...
     >>> ax = matplot_pressure_z(my_fluid_flow, theta=int(my_fluid_flow.ntheta/2))
@@ -259,8 +245,8 @@ def matplot_pressure_z(pressure_matrix_object, theta=0, ax=None):
     >>> # plt.show()
     """
     if (
-        not pressure_matrix_object.numerical_pressure_matrix_available
-        and not pressure_matrix_object.analytical_pressure_matrix_available
+            not fluid_flow_object.numerical_pressure_matrix_available
+            and not fluid_flow_object.analytical_pressure_matrix_available
     ):
         raise ValueError(
             "Must calculate the pressure matrix. "
@@ -268,15 +254,15 @@ def matplot_pressure_z(pressure_matrix_object, theta=0, ax=None):
         )
     if ax is None:
         ax = plt.gca()
-    y_n = np.zeros(pressure_matrix_object.nz)
-    y_a = np.zeros(pressure_matrix_object.nz)
-    for i in range(0, pressure_matrix_object.nz):
-        y_n[i] = pressure_matrix_object.p_mat_numerical[i][theta]
-        y_a[i] = pressure_matrix_object.p_mat_analytical[i][theta]
-    if pressure_matrix_object.numerical_pressure_matrix_available:
-        ax.plot(pressure_matrix_object.z_list, y_n, "b", label="Numerical pressure")
-    if pressure_matrix_object.analytical_pressure_matrix_available:
-        ax.plot(pressure_matrix_object.z_list, y_a, "r", label="Analytical pressure")
+    y_n = np.zeros(fluid_flow_object.nz)
+    y_a = np.zeros(fluid_flow_object.nz)
+    for i in range(0, fluid_flow_object.nz):
+        y_n[i] = fluid_flow_object.p_mat_numerical[i][theta]
+        y_a[i] = fluid_flow_object.p_mat_analytical[i][theta]
+    if fluid_flow_object.numerical_pressure_matrix_available:
+        ax.plot(fluid_flow_object.z_list, y_n, "b", label="Numerical pressure")
+    if fluid_flow_object.analytical_pressure_matrix_available:
+        ax.plot(fluid_flow_object.z_list, y_a, "r", label="Analytical pressure")
     ax.set_title(
         "Pressure along the Z direction (direction of flow); Theta=" + str(theta)
     )
@@ -285,11 +271,11 @@ def matplot_pressure_z(pressure_matrix_object, theta=0, ax=None):
     return ax
 
 
-def matplot_shape(pressure_matrix_object, theta=0, ax=None):
+def matplot_shape(fluid_flow_object, theta=0, ax=None):
     """This function assembles a graphic representing the geometry of the rotor using matplotlib.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     theta: int, optional
         The theta to be considered.
     ax : matplotlib axes, optional
@@ -300,35 +286,33 @@ def matplot_shape(pressure_matrix_object, theta=0, ax=None):
         Returns the axes object with the plot.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> ax = matplot_shape(my_fluid_flow, theta=int(my_fluid_flow.ntheta/2))
     >>> # to show the plots you can use:
     >>> # plt.show()
     """
     if ax is None:
         ax = plt.gca()
-    y_ext = np.zeros(pressure_matrix_object.nz)
-    y_int = np.zeros(pressure_matrix_object.nz)
-    for i in range(0, pressure_matrix_object.nz):
-        y_ext[i] = pressure_matrix_object.re[i][theta]
-        y_int[i] = pressure_matrix_object.ri[i][theta]
-    ax.plot(pressure_matrix_object.z_list, y_ext, "r")
-    ax.plot(pressure_matrix_object.z_list, y_int, "b")
+    y_ext = np.zeros(fluid_flow_object.nz)
+    y_int = np.zeros(fluid_flow_object.nz)
+    for i in range(0, fluid_flow_object.nz):
+        y_ext[i] = fluid_flow_object.re[i][theta]
+        y_int[i] = fluid_flow_object.ri[i][theta]
+    ax.plot(fluid_flow_object.z_list, y_ext, "r")
+    ax.plot(fluid_flow_object.z_list, y_int, "b")
     ax.set_title("Shapes of stator and rotor along Z; Theta=" + str(theta))
     ax.set_xlabel("Points along Z")
     ax.set_ylabel("Radial direction")
     return ax
 
 
-def matplot_pressure_theta_cylindrical(
-    pressure_matrix_object, z=0, from_numerical=True, ax=None
-):
+def matplot_pressure_theta_cylindrical(fluid_flow_object, z=0, from_numerical=True, ax=None):
     """This function assembles cylindrical pressure graphic in the theta direction for a given z,
     using matplotlib.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     z: int, optional
         The distance along z-axis to be considered.
     from_numerical: bool, optional
@@ -344,8 +328,8 @@ def matplot_pressure_theta_cylindrical(
         Returns the axes object with the plot.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
     array([[...
     >>> ax = matplot_pressure_theta_cylindrical(my_fluid_flow, z=int(my_fluid_flow.nz/2))
@@ -353,51 +337,45 @@ def matplot_pressure_theta_cylindrical(
     >>> # plt.show()
     """
     if (
-        not pressure_matrix_object.numerical_pressure_matrix_available
-        and not pressure_matrix_object.analytical_pressure_matrix_available
+            not fluid_flow_object.numerical_pressure_matrix_available
+            and not fluid_flow_object.analytical_pressure_matrix_available
     ):
         raise ValueError(
             "Must calculate the pressure matrix. "
             "Try calling calculate_pressure_matrix_numerical() or calculate_pressure_matrix_analytical() first."
         )
     if from_numerical:
-        if pressure_matrix_object.numerical_pressure_matrix_available:
-            p_mat = pressure_matrix_object.p_mat_numerical
+        if fluid_flow_object.numerical_pressure_matrix_available:
+            p_mat = fluid_flow_object.p_mat_numerical
         else:
-            p_mat = pressure_matrix_object.p_mat_analytical
+            p_mat = fluid_flow_object.p_mat_analytical
     else:
-        if pressure_matrix_object.analytical_pressure_matrix_available:
-            p_mat = pressure_matrix_object.p_mat_analytical
+        if fluid_flow_object.analytical_pressure_matrix_available:
+            p_mat = fluid_flow_object.p_mat_analytical
         else:
-            p_mat = pressure_matrix_object.p_mat_numerical
+            p_mat = fluid_flow_object.p_mat_numerical
     if ax is None:
         fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
     r = np.arange(
         0,
-        pressure_matrix_object.radius_stator + 0.0001,
-        (pressure_matrix_object.radius_stator - pressure_matrix_object.radius_rotor)
-        / pressure_matrix_object.nradius,
-    )
-    theta = np.arange(
-        -np.pi * 0.25,
-        1.75 * np.pi + pressure_matrix_object.dtheta / 2,
-        pressure_matrix_object.dtheta,
+        fluid_flow_object.radius_stator + 0.0001,
+        (fluid_flow_object.radius_stator - fluid_flow_object.radius_rotor) / fluid_flow_object.nradius,
     )
+    theta = np.arange(-np.pi * 0.25, 1.75 * np.pi + fluid_flow_object.dtheta / 2, fluid_flow_object.dtheta)
 
-    pressure_along_theta = np.zeros(pressure_matrix_object.ntheta)
-    for i in range(0, pressure_matrix_object.ntheta):
+    pressure_along_theta = np.zeros(fluid_flow_object.ntheta)
+    for i in range(0, fluid_flow_object.ntheta):
         pressure_along_theta[i] = p_mat[0][i]
 
     min_pressure = np.amin(pressure_along_theta)
 
     r_matrix, theta_matrix = np.meshgrid(r, theta)
     z_matrix = np.zeros((theta.size, r.size))
-    inner_radius_list = np.zeros(pressure_matrix_object.ntheta)
+    inner_radius_list = np.zeros(fluid_flow_object.ntheta)
     pressure_list = np.zeros((theta.size, r.size))
     for i in range(0, theta.size):
         inner_radius = np.sqrt(
-            pressure_matrix_object.xri[z][i] * pressure_matrix_object.xri[z][i]
-            + pressure_matrix_object.yri[z][i] * pressure_matrix_object.yri[z][i]
+            fluid_flow_object.xri[z][i] * fluid_flow_object.xri[z][i] + fluid_flow_object.yri[z][i] * fluid_flow_object.yri[z][i]
         )
         inner_radius_list[i] = inner_radius
         for j in range(0, r.size):
@@ -410,12 +388,12 @@ def matplot_pressure_theta_cylindrical(
     return ax
 
 
-def matplot_pressure_theta(pressure_matrix_object, z=0, ax=None):
+def matplot_pressure_theta(fluid_flow_object, z=0, ax=None):
     """This function assembles pressure graphic in the theta direction for a given z,
     using matplotlib.
     Parameters
     ----------
-    pressure_matrix_object: a PressureMatrix object
+    fluid_flow_object: a FluidFlow object
     z: int, optional
         The distance along z-axis to be considered.
     ax : matplotlib axes, optional
@@ -426,8 +404,8 @@ def matplot_pressure_theta(pressure_matrix_object, z=0, ax=None):
         Returns the axes object with the plot.
     Examples
     --------
-    >>> from ross.fluid_flow.fluid_flow import pressure_matrix_example
-    >>> my_fluid_flow = pressure_matrix_example()
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example
+    >>> my_fluid_flow = fluid_flow_example()
     >>> my_fluid_flow.calculate_pressure_matrix_numerical() # doctest: +ELLIPSIS
     array([[...
     >>> ax = matplot_pressure_theta(my_fluid_flow, z=int(my_fluid_flow.nz/2))
@@ -435,8 +413,8 @@ def matplot_pressure_theta(pressure_matrix_object, z=0, ax=None):
     >>> # plt.show()
     """
     if (
-        not pressure_matrix_object.numerical_pressure_matrix_available
-        and not pressure_matrix_object.analytical_pressure_matrix_available
+            not fluid_flow_object.numerical_pressure_matrix_available
+            and not fluid_flow_object.analytical_pressure_matrix_available
     ):
         raise ValueError(
             "Must calculate the pressure matrix. "
@@ -444,17 +422,17 @@ def matplot_pressure_theta(pressure_matrix_object, z=0, ax=None):
         )
     if ax is None:
         ax = plt.gca()
-    if pressure_matrix_object.numerical_pressure_matrix_available:
+    if fluid_flow_object.numerical_pressure_matrix_available:
         ax.plot(
-            pressure_matrix_object.gama[z],
-            pressure_matrix_object.p_mat_numerical[z],
+            fluid_flow_object.gama[z],
+            fluid_flow_object.p_mat_numerical[z],
             "b",
-            label="Numerical pressure",
+            label="Numerical pressure"
         )
-    if pressure_matrix_object.analytical_pressure_matrix_available:
+    if fluid_flow_object.analytical_pressure_matrix_available:
         ax.plot(
-            pressure_matrix_object.gama[z],
-            pressure_matrix_object.p_mat_analytical[z],
+            fluid_flow_object.gama[z],
+            fluid_flow_object.p_mat_analytical[z],
             "r",
             label="Analytical pressure",
         )
@@ -462,3 +440,5 @@ def matplot_pressure_theta(pressure_matrix_object, z=0, ax=None):
     ax.set_xlabel("Points along Theta")
     ax.set_ylabel("Pressure")
     return ax
+
+


Test patch:
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
index 61c6ad94..db62ac99 100644
--- a/ross/tests/test_fluid_flow.py
+++ b/ross/tests/test_fluid_flow.py
@@ -24,9 +24,9 @@ def fluid_flow_short_eccentricity():
     eccentricity = (radius_stator - radius_rotor)*0.2663
     visc = 0.1
     rho = 860.
-    return flow.PressureMatrix(nz, ntheta, nradius, length, omega, p_in,
-                               p_out, radius_rotor, radius_stator,
-                               visc, rho, eccentricity=eccentricity)
+    return flow.FluidFlow(nz, ntheta, nradius, length, omega, p_in,
+                          p_out, radius_rotor, radius_stator,
+                          visc, rho, eccentricity=eccentricity)
 
 
 def fluid_flow_short_load():
@@ -42,9 +42,9 @@ def fluid_flow_short_load():
     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)
+    return flow.FluidFlow(nz, ntheta, nradius, length, omega, p_in,
+                          p_out, radius_rotor, radius_stator,
+                          visc, rho, load=load)
 
 
 def test_sommerfeld_number():
@@ -112,11 +112,11 @@ def fluid_flow_short_numerical():
     radius_stator = 0.1
     visc = 0.015
     rho = 860.
-    beta = np.pi
+    attitude_angle = None
     eccentricity = 0.001
-    return flow.PressureMatrix(nz, ntheta, nradius, length,
-                               omega, p_in, p_out, radius_rotor,
-                               radius_stator, visc, rho, beta=beta, eccentricity=eccentricity)
+    return flow.FluidFlow(nz, ntheta, nradius, length,
+                          omega, p_in, p_out, radius_rotor,
+                          radius_stator, visc, rho, attitude_angle=attitude_angle, eccentricity=eccentricity)
 
 
 def fluid_flow_long_numerical():
@@ -129,14 +129,13 @@ def fluid_flow_long_numerical():
     radius_rotor = 1
     h = 0.000194564
     radius_stator = radius_rotor + h
-    length = 8 * radius_stator
+    length = 8 * 2 * radius_stator
     visc = 0.015
     rho = 860.
-    beta = np.pi
     eccentricity = 0.0001
-    return flow.PressureMatrix(nz, ntheta, nradius, length,
-                               omega, p_in, p_out, radius_rotor,
-                               radius_stator, visc, rho, beta=beta, eccentricity=eccentricity)
+    return flow.FluidFlow(nz, ntheta, nradius, length,
+                          omega, p_in, p_out, radius_rotor,
+                          radius_stator, visc, rho, eccentricity=eccentricity)
 
 
 def test_numerical_abs_error():
@@ -163,29 +162,29 @@ def test_long_bearing():
     bearing.calculate_pressure_matrix_numerical()
     error = (max(bearing.p_mat_analytical[int(bearing.nz / 2)]) -
              max(bearing.p_mat_numerical[int(bearing.nz / 2)])) / max(bearing.p_mat_numerical[int(bearing.nz / 2)])
-    assert math.isclose(error, 0, abs_tol=0.02)
+    assert math.isclose(error, 0, abs_tol=0.007)
 
 
 def test_oil_film_force_short():
     bearing = fluid_flow_short_numerical()
     bearing.calculate_pressure_matrix_numerical()
-    n, t = calculate_oil_film_force(bearing)
-    n_numerical, t_numerical = calculate_oil_film_force(bearing, force_type='numerical')
-    error_n = (n - n_numerical) / n_numerical
-    error_t = (t - t_numerical) / t_numerical
-    assert_allclose(error_n, 0, atol=0.009)
-    assert_allclose(error_t, 0, atol=0.7)
+    n, t, force_x, force_y = calculate_oil_film_force(bearing)
+    n_numerical, t_numerical, force_x_numerical, force_y_numerical = \
+        calculate_oil_film_force(bearing, force_type='numerical')
+    assert_allclose(n, n_numerical, rtol=0.08)
+    assert_allclose(t, t_numerical, rtol=0.7)
+    assert_allclose(force_x_numerical, 0, atol=1e-06)
+    assert_allclose(force_y_numerical, bearing.load, atol=1e-05)
 
 
 def test_oil_film_force_long():
     bearing = fluid_flow_long_numerical()
     bearing.calculate_pressure_matrix_numerical()
-    n, t = calculate_oil_film_force(bearing)
-    n_numerical, t_numerical = calculate_oil_film_force(bearing, force_type='numerical')
-    error_n = (n - n_numerical) / n_numerical
-    error_t = (t - t_numerical) / t_numerical
-    assert_allclose(error_n, 0, atol=0.2)
-    assert_allclose(error_t, 0, atol=0.4)
+    n, t, force_x, force_y = calculate_oil_film_force(bearing)
+    n_numerical, t_numerical, force_x_numerical, force_y_numerical = \
+        calculate_oil_film_force(bearing, force_type='numerical')
+    assert_allclose(n, n_numerical, rtol=0.2)
+    assert_allclose(t, t_numerical, rtol=0.4)
 
 
 def test_bokeh_plots():
@@ -206,3 +205,4 @@ def test_matplotlib_plots():
     assert isinstance(matplot_pressure_theta(bearing), ax_type)
     assert isinstance(matplot_pressure_z(bearing), ax_type)
     assert isinstance(matplot_shape(bearing), ax_type)
+