Repository: ross-rotordynamics/ross
Pull request #899 (2022-09-01 13:20:17)

Cylindrical bearing starvation
Updating cylindrical bearing code to be able to run analysis with starvation conditions.


Patch:
diff --git a/ross/fluid_flow/cylindrical.py b/ross/fluid_flow/cylindrical.py
index cf69448d..21c703b2 100644
--- a/ross/fluid_flow/cylindrical.py
+++ b/ross/fluid_flow/cylindrical.py
@@ -2,10 +2,9 @@ import time
 
 import numpy as np
 from numpy.linalg import pinv
-from scipy.optimize import curve_fit, minimize
-
 from ross.bearing_seal_element import BearingElement
 from ross.units import Q_, check_units
+from scipy.optimize import curve_fit, minimize
 
 
 class THDCylindrical(BearingElement):
@@ -54,6 +53,10 @@ class THDCylindrical(BearingElement):
         Load in X direction. The unit is newton.
     load_y_direction : Float
         Load in Y direction. The unit is newton.
+    operating_type : string
+        Choose the operating condition that bearing is operating.
+        - 'flooded'
+        - 'starvation'
 
     Fluid propierties
     ^^^^^^^^^^^^^^^^^
@@ -64,12 +67,15 @@ class THDCylindrical(BearingElement):
         - 'ISOVG46'
     reference_temperature : float
         Oil reference temperature. The unit is celsius.
-    reference_viscosity : float
-        Oil viscosity at reference temperature. Unit is Pa.s.
     groove_factor : list, numpy array, tuple or float
         Ratio of oil in reservoir temperature that mixes with the circulating oil.
         Is required one factor per segment.
-
+    oil_flow: float
+        Suply oil flow to bearing. Only used when operating type 'starvation' is
+        selected. Unit is Litre per minute (l/min)
+    injection_pressure: float
+        Suply oil pressure that bearing receives at groove regions. Only used
+        when operating type 'starvation' is selected. Unit is Pascal (Pa).
 
     Turbulence Model
     ^^^^^^^^^^^^^^^^
@@ -125,7 +131,7 @@ class THDCylindrical(BearingElement):
     >>> from ross.fluid_flow.cylindrical import cylindrical_bearing_example
     >>> bearing = cylindrical_bearing_example()
     >>> bearing.equilibrium_pos
-    array([ 0.60678516, -0.73288691])
+    array([ 0.6873316 , -0.79393636])
     """
 
     @check_units
@@ -136,11 +142,9 @@ class THDCylindrical(BearingElement):
         radial_clearance,
         elements_circumferential,
         elements_axial,
-        n_y,
         n_pad,
         pad_arc_length,
         reference_temperature,
-        reference_viscosity,
         speed,
         load_x_direction,
         load_y_direction,
@@ -150,6 +154,9 @@ class THDCylindrical(BearingElement):
         sommerfeld_type=2,
         initial_guess=[0.1, -0.1],
         method="perturbation",
+        operating_type="flooded",
+        injection_pressure=None,
+        oil_flow=None,
         show_coef=False,
         print_result=False,
         print_progress=False,
@@ -161,10 +168,8 @@ class THDCylindrical(BearingElement):
         self.radial_clearance = radial_clearance
         self.elements_circumferential = elements_circumferential
         self.elements_axial = elements_axial
-        self.n_y = n_y
         self.n_pad = n_pad
         self.reference_temperature = reference_temperature
-        self.reference_viscosity = reference_viscosity
         self.load_x_direction = load_x_direction
         self.load_y_direction = load_y_direction
         self.lubricant = lubricant
@@ -173,14 +178,14 @@ class THDCylindrical(BearingElement):
         self.sommerfeld_type = sommerfeld_type
         self.initial_guess = initial_guess
         self.method = method
+        self.operating_type = operating_type
+        self.injection_pressure = injection_pressure
+        self.oil_flow = oil_flow
         self.show_coef = show_coef
         self.print_result = print_result
         self.print_progress = print_progress
         self.print_time = print_time
 
-        if self.n_y == None:
-            self.n_y = self.elements_circumferential
-
         self.betha_s_dg = pad_arc_length
         self.betha_s = pad_arc_length * np.pi / 180
 
@@ -191,7 +196,7 @@ class THDCylindrical(BearingElement):
         ##
         # Dimensionless discretization variables
 
-        self.dY = 1 / self.n_y
+        self.dY = 1 / self.elements_circumferential
         self.dZ = 1 / self.elements_axial
 
         # Z-axis direction
@@ -214,6 +219,8 @@ class THDCylindrical(BearingElement):
 
         self.Zdim = self.Z * self.axial_length
 
+        self.oil_flow = self.oil_flow / 60000
+
         self.lubricant_dict = {
             "ISOVG32": {
                 "viscosity1": Q_(4.05640e-06, "reyn").to_base_units().m,
@@ -234,13 +241,13 @@ class THDCylindrical(BearingElement):
                 "lube_conduct": Q_(0.15, "W/(m*degC)").to_base_units().m,
             },
             "TEST": {
-                "viscosity1": Q_(2.758e-6, "reyn").to_base_units().m,
-                "temp1": Q_(121.7, "degF").to_base_units().m,
-                "viscosity2": Q_(1.119e-6, "reyn").to_base_units().m,
-                "temp2": Q_(175.7, "degF").to_base_units().m,
-                "lube_density": Q_(8.0e-5, "lbf*s²/in⁴").to_base_units().m,
-                "lube_cp": Q_(1.79959e2, "BTU*in/(lbf*s²*degF)").to_base_units().m,
-                "lube_conduct": Q_(2.00621e-6, "BTU/(in*s*degF)").to_base_units().m,
+                "viscosity1": Q_(0.04, "Pa*s").to_base_units().m,
+                "temp1": Q_(40, "degC").to_base_units().m,
+                "viscosity2": Q_(0.01, "Pa*s").to_base_units().m,
+                "temp2": Q_(100, "degC").to_base_units().m,
+                "lube_density": Q_(863.61302696, "kg/m³").to_base_units().m,
+                "lube_cp": Q_(1951.88616, "J/(kg*degK)").to_base_units().m,
+                "lube_conduct": Q_(0.15, "W/(m*degC)").to_base_units().m,
             },
         }
 
@@ -256,6 +263,8 @@ class THDCylindrical(BearingElement):
         # Interpolation coefficients
         self.a, self.b = self._interpol(T_muI, T_muF, mu_I, mu_F)
 
+        self.reference_viscosity = self.a * (self.reference_temperature**self.b)
+
         number_of_freq = np.shape(speed)[0]
 
         kxx = np.zeros(number_of_freq)
@@ -292,6 +301,692 @@ class THDCylindrical(BearingElement):
 
         super().__init__(node, kxx, cxx, kyy, kxy, kyx, cyy, cxy, cyx, speed)
 
+    def _flooded(self, n_p, Mat_coef, b_P, mu):
+        """Provides an analysis in which the bearing always receive sufficient oil feed to operate.
+
+        Parameters
+        ----------
+        n_p : integer,
+           current pad in analysis.
+        Mat_coef : np.array
+            Coeficient matrix.
+        b_P: np.array
+            Coefficients to pressure independent terms.
+        mu : np.array
+            Viscosity matrix.
+
+        Returns
+        -------
+        self.P : np.array
+            Pressure distribution in current pad vector.
+        """
+
+        ki = 0
+        kj = 0
+        k = 0
+        for ii in np.arange((self.Z_I + 0.5 * self.dZ), self.Z_F, self.dZ):
+            for jj in np.arange(
+                self.thetaI[n_p] + (self.dtheta / 2),
+                self.thetaF[n_p],
+                self.dtheta,
+            ):
+
+                hP = 1 - self.X * np.cos(jj) - self.Y * np.sin(jj)
+                he = (
+                    1
+                    - self.X * np.cos(jj + 0.5 * self.dtheta)
+                    - self.Y * np.sin(jj + 0.5 * self.dtheta)
+                )
+                hw = (
+                    1
+                    - self.X * np.cos(jj - 0.5 * self.dtheta)
+                    - self.Y * np.sin(jj - 0.5 * self.dtheta)
+                )
+                hn = hP
+                hs = hn
+
+                if kj == 0 and ki == 0:
+                    MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                    MU_w = mu[ki, kj]
+                    MU_s = mu[ki, kj]
+                    MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                if kj == 0 and ki > 0 and ki < self.elements_axial - 1:
+                    MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                    MU_w = mu[ki, kj]
+                    MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                    MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                if kj == 0 and ki == self.elements_axial - 1:
+                    MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                    MU_w = mu[ki, kj]
+                    MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                    MU_n = mu[ki, kj]
+
+                if ki == 0 and kj > 0 and kj < self.elements_circumferential - 1:
+                    MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                    MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                    MU_s = mu[ki, kj]
+                    MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                if (
+                    kj > 0
+                    and kj < self.elements_circumferential - 1
+                    and ki > 0
+                    and ki < self.elements_axial - 1
+                ):
+                    MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                    MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                    MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                    MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                if (
+                    ki == self.elements_axial - 1
+                    and kj > 0
+                    and kj < self.elements_circumferential - 1
+                ):
+                    MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                    MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                    MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                    MU_n = mu[ki, kj]
+
+                if ki == 0 and kj == self.elements_circumferential - 1:
+                    MU_e = mu[ki, kj]
+                    MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                    MU_s = mu[ki, kj]
+                    MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                if (
+                    kj == self.elements_circumferential - 1
+                    and ki > 0
+                    and ki < self.elements_axial - 1
+                ):
+                    MU_e = mu[ki, kj]
+                    MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                    MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                    MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                if (
+                    kj == self.elements_circumferential - 1
+                    and ki == self.elements_axial - 1
+                ):
+                    MU_e = mu[ki, kj]
+                    MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                    MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                    MU_n = mu[ki, kj]
+
+                CE = (self.dZ * he**3) / (
+                    12 * MU_e[n_p] * self.dY * self.betha_s**2
+                )
+                CW = (self.dZ * hw**3) / (
+                    12 * MU_w[n_p] * self.dY * self.betha_s**2
+                )
+                CN = (self.dY * (self.journal_radius**2) * hn**3) / (
+                    12 * MU_n[n_p] * self.dZ * self.axial_length**2
+                )
+                CS = (self.dY * (self.journal_radius**2) * hs**3) / (
+                    12 * MU_s[n_p] * self.dZ * self.axial_length**2
+                )
+                CP = -(CE + CW + CN + CS)
+
+                B = (self.dZ / (2 * self.betha_s)) * (he - hw) - (
+                    (self.Xpt * np.cos(jj) + self.Ypt * np.sin(jj)) * self.dy * self.dZ
+                )
+
+                k = k + 1
+                b_P[k - 1, 0] = B
+
+                if ki == 0 and kj == 0:
+                    Mat_coef[k - 1, k - 1] = CP - CS - CW
+                    Mat_coef[k - 1, k] = CE
+                    Mat_coef[k - 1, k + self.elements_circumferential - 1] = CN
+
+                elif kj == 0 and ki > 0 and ki < self.elements_axial - 1:
+                    Mat_coef[k - 1, k - 1] = CP - CW
+                    Mat_coef[k - 1, k] = CE
+                    Mat_coef[k - 1, k - self.elements_circumferential - 1] = CS
+                    Mat_coef[k - 1, k + self.elements_circumferential - 1] = CN
+
+                elif kj == 0 and ki == self.elements_axial - 1:
+                    Mat_coef[k - 1, k - 1] = CP - CN - CW
+                    Mat_coef[k - 1, k] = CE
+                    Mat_coef[k - 1, k - self.elements_circumferential - 1] = CS
+
+                elif ki == 0 and kj > 0 and kj < self.elements_circumferential - 1:
+                    Mat_coef[k - 1, k - 1] = CP - CS
+                    Mat_coef[k - 1, k] = CE
+                    Mat_coef[k - 1, k - 2] = CW
+                    Mat_coef[k - 1, k + self.elements_circumferential - 1] = CN
+
+                elif (
+                    ki > 0
+                    and ki < self.elements_axial - 1
+                    and kj > 0
+                    and kj < self.elements_circumferential - 1
+                ):
+                    Mat_coef[k - 1, k - 1] = CP
+                    Mat_coef[k - 1, k - 2] = CW
+                    Mat_coef[k - 1, k - self.elements_circumferential - 1] = CS
+                    Mat_coef[k - 1, k + self.elements_circumferential - 1] = CN
+                    Mat_coef[k - 1, k] = CE
+
+                elif (
+                    ki == self.elements_axial - 1
+                    and kj > 0
+                    and kj < self.elements_circumferential - 1
+                ):
+                    Mat_coef[k - 1, k - 1] = CP - CN
+                    Mat_coef[k - 1, k] = CE
+                    Mat_coef[k - 1, k - 2] = CW
+                    Mat_coef[k - 1, k - self.elements_circumferential - 1] = CS
+
+                elif ki == 0 and kj == self.elements_circumferential - 1:
+                    Mat_coef[k - 1, k - 1] = CP - CE - CS
+                    Mat_coef[k - 1, k - 2] = CW
+                    Mat_coef[k - 1, k + self.elements_circumferential - 1] = CN
+
+                elif (
+                    kj == self.elements_circumferential - 1
+                    and ki > 0
+                    and ki < self.elements_axial - 1
+                ):
+                    Mat_coef[k - 1, k - 1] = CP - CE
+                    Mat_coef[k - 1, k - 2] = CW
+                    Mat_coef[k - 1, k - self.elements_circumferential - 1] = CS
+                    Mat_coef[k - 1, k + self.elements_circumferential - 1] = CN
+
+                elif (
+                    ki == self.elements_axial - 1
+                    and kj == self.elements_circumferential - 1
+                ):
+                    Mat_coef[k - 1, k - 1] = CP - CE - CN
+                    Mat_coef[k - 1, k - 2] = CW
+                    Mat_coef[k - 1, k - self.elements_circumferential - 1] = CS
+
+                kj = kj + 1
+
+            kj = 0
+            ki = ki + 1
+
+        # Solution of pressure field end
+
+        p = np.linalg.solve(Mat_coef, b_P)
+        cont = 0
+
+        for i in np.arange(self.elements_axial):
+            for j in np.arange(self.elements_circumferential):
+
+                self.P[i, j, n_p] = p[cont]
+                cont = cont + 1
+
+                if self.P[i, j, n_p] < 0:
+                    self.P[i, j, n_p] = 0
+
+        # Dimensional pressure fied
+
+        self.Pdim = (
+            self.P * self.reference_viscosity * self.speed * (self.journal_radius**2)
+        ) / (self.radial_clearance**2)
+
+        return self.P
+
+    def _starvation(self, n_p, Mat_coef_st, mu, p_old, p, B, B_theta, nk):
+        """Provides an analysis in which the bearing may receive insufficient oil feed.
+
+        Parameters
+        ----------
+        n_p : integer,
+           current pad in analysis.
+        Mat_coef_st : np.array
+            Coeficient matrix.
+        mu : np.array
+            Viscosity matrix.
+        p_old : np.array
+            Past pressure matrix.
+        p : np.array
+            Current pressure matrix.
+        B: np.array
+            Coefficients to independent terms.
+        B: np.array
+            Coefficients to volumetric fraction independent terms.
+        nk: integer
+            counter.
+
+
+        Returns
+        -------
+        self.P : np.array
+            Pressure distribution in current pad vector.
+        """
+
+        while self.erro >= 0.01:
+
+            p_old = np.array(p)
+
+            theta_vol_old = np.array(self.theta_vol)
+
+            # cont = cont+1
+
+            k = 0
+            ki = 0
+            kj = 0
+
+            for ii in np.arange((self.Z_I + 0.5 * self.dZ), self.Z_F, self.dZ):
+                for jj in np.arange(
+                    self.thetaI[n_p] + (self.dtheta / 2),
+                    self.thetaF[n_p],
+                    self.dtheta,
+                ):
+
+                    hP = 1 - self.X * np.cos(jj) - self.Y * np.sin(jj)
+                    he = (
+                        1
+                        - self.X * np.cos(jj + 0.5 * self.dtheta)
+                        - self.Y * np.sin(jj + 0.5 * self.dtheta)
+                    )
+                    hw = (
+                        1
+                        - self.X * np.cos(jj - 0.5 * self.dtheta)
+                        - self.Y * np.sin(jj - 0.5 * self.dtheta)
+                    )
+                    hn = hP
+                    hs = hn
+
+                    hpt = -self.Xpt * np.cos(jj) - self.Ypt * np.sin(jj)
+
+                    if kj == 0 and ki == 0:
+                        MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                        MU_w = mu[ki, kj]
+                        MU_s = mu[ki, kj]
+                        MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                    if kj == 0 and ki > 0 and ki < self.elements_axial - 1:
+                        MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                        MU_w = mu[ki, kj]
+                        MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                        MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                    if kj == 0 and ki == self.elements_axial - 1:
+                        MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                        MU_w = mu[ki, kj]
+                        MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                        MU_n = mu[ki, kj]
+
+                    if ki == 0 and kj > 0 and kj < self.elements_circumferential - 1:
+                        MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                        MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                        MU_s = mu[ki, kj]
+                        MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                    if (
+                        kj > 0
+                        and kj < self.elements_circumferential - 1
+                        and ki > 0
+                        and ki < self.elements_axial - 1
+                    ):
+                        MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                        MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                        MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                        MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                    if (
+                        ki == self.elements_axial - 1
+                        and kj > 0
+                        and kj < self.elements_circumferential - 1
+                    ):
+                        MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
+                        MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                        MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                        MU_n = mu[ki, kj]
+
+                    if ki == 0 and kj == self.elements_circumferential - 1:
+                        MU_e = mu[ki, kj]
+                        MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                        MU_s = mu[ki, kj]
+                        MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                    if (
+                        kj == self.elements_circumferential - 1
+                        and ki > 0
+                        and ki < self.elements_axial - 1
+                    ):
+                        MU_e = mu[ki, kj]
+                        MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                        MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                        MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+
+                    if (
+                        kj == self.elements_circumferential - 1
+                        and ki == self.elements_axial - 1
+                    ):
+                        MU_e = mu[ki, kj]
+                        MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
+                        MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
+                        MU_n = mu[ki, kj]
+
+                    CE = (self.dZ * he**3) / (
+                        12 * MU_e[n_p] * self.dY * self.betha_s**2
+                    )
+                    CW = (self.dZ * hw**3) / (
+                        12 * MU_w[n_p] * self.dY * self.betha_s**2
+                    )
+                    CN = (self.dY * (self.journal_radius**2) * hn**3) / (
+                        12 * MU_n[n_p] * self.dZ * self.axial_length**2
+                    )
+                    CS = (self.dY * (self.journal_radius**2) * hs**3) / (
+                        12 * MU_s[n_p] * self.dZ * self.axial_length**2
+                    )
+                    CP = -(CE + CW + CN + CS)
+
+                    # Termo Fonte
+                    KP1 = -(self.dZ / (2 * self.betha_s)) * he
+
+                    KP2 = -hpt * self.dY * self.dZ
+
+                    KP = KP1 + KP2
+
+                    KW = (self.dZ / (2 * self.betha_s)) * hw
+
+                    if (
+                        kj > 0
+                        and kj < self.elements_circumferential - 1
+                        and ki > 0
+                        and ki < self.elements_axial - 1
+                    ):  # Center region
+
+                        Mat_coef_st[k, k] = CP
+                        Mat_coef_st[k, k + 1] = CE
+                        Mat_coef_st[k, k - 1] = CW
+                        Mat_coef_st[k, k + self.elements_circumferential] = CN
+                        Mat_coef_st[k, k - self.elements_circumferential] = CS
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = -KP * self.theta_vol[k] - KW * self.theta_vol[k - 1]
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - KW * self.theta_vol[k - 1]
+                            ) / KP
+
+                    elif (
+                        kj > 0 and kj < self.elements_circumferential - 1 and ki == 0
+                    ):  # Inferior edge without corners
+
+                        Mat_coef_st[k, k] = CP - CS
+                        Mat_coef_st[k, k + 1] = CE
+                        Mat_coef_st[k, k - 1] = CW
+                        Mat_coef_st[k, k + self.elements_circumferential] = CN
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = -KP * self.theta_vol[k] - KW * self.theta_vol[k - 1]
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - KW * self.theta_vol[k - 1]
+                            ) / KP
+
+                    elif (
+                        kj > 0
+                        and kj < self.elements_circumferential - 1
+                        and ki == self.elements_axial - 1
+                    ):  # Superior edge without corners
+
+                        Mat_coef_st[k, k] = CP - CN
+                        Mat_coef_st[k, k + 1] = CE
+                        Mat_coef_st[k, k - 1] = CW
+                        Mat_coef_st[k, k - self.elements_circumferential] = CS
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = -KP * self.theta_vol[k] - KW * self.theta_vol[k - 1]
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - KW * self.theta_vol[k - 1]
+                            ) / KP
+
+                    elif (
+                        kj == 0 and ki > 0 and ki < self.elements_axial - 1
+                    ):  # Left edge without corners
+
+                        Mat_coef_st[k, k] = CP - CW
+                        Mat_coef_st[k, k + 1] = CE
+                        Mat_coef_st[k, k - self.elements_circumferential] = CS
+                        Mat_coef_st[k, k + self.elements_circumferential] = CN
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = (
+                                -KP * self.theta_vol[k]
+                                - self.theta_vol_groove[n_p] * KW
+                                - 2 * CW * self.injection_pressure
+                            )
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - self.theta_vol_groove[n_p] * KW * 1
+                            ) / KP
+
+                    elif (
+                        kj == self.elements_circumferential - 1
+                        and ki > 0
+                        and ki < self.elements_axial - 1
+                    ):  # Right edge without corners
+
+                        Mat_coef_st[k, k] = CP - CE
+
+                        Mat_coef_st[k, k - 1] = CW
+                        Mat_coef_st[k, k + self.elements_circumferential] = CN
+                        Mat_coef_st[k, k - self.elements_circumferential] = CS
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = -KP * self.theta_vol[k] - KW * self.theta_vol[k - 1]
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - KW * self.theta_vol[k - 1]
+                            ) / KP
+
+                    elif kj == 0 and ki == 0:  # Corner inferior left
+
+                        Mat_coef_st[k, k] = CP - CS - CW
+                        Mat_coef_st[k, k + 1] = CE
+                        Mat_coef_st[k, k + self.elements_circumferential] = CN
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = (
+                                -KP * self.theta_vol[k]
+                                - self.theta_vol_groove[n_p] * KW
+                                - 2 * CW * self.injection_pressure
+                            )
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - self.theta_vol_groove[n_p] * KW * 1
+                            ) / KP
+
+                    elif (
+                        kj == self.elements_circumferential - 1 and ki == 0
+                    ):  # Corner inferior right
+
+                        Mat_coef_st[k, k] = CP - CS - CE
+                        Mat_coef_st[k, k - 1] = CW
+                        Mat_coef_st[k, k + self.elements_circumferential] = CN
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = -KP * self.theta_vol[k] - KW * self.theta_vol[k - 1]
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - KW * self.theta_vol[k - 1]
+                            ) / KP
+
+                    elif (
+                        kj == 0 and ki == self.elements_axial - 1
+                    ):  # Corner superior left
+
+                        Mat_coef_st[k, k] = CP - CN - CW
+                        Mat_coef_st[k, k + 1] = CE
+                        Mat_coef_st[k, k - self.elements_circumferential] = CS
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = (
+                                -KP * self.theta_vol[k]
+                                - self.theta_vol_groove[n_p] * KW
+                                - 2 * CW * self.injection_pressure
+                            )
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - self.theta_vol_groove[n_p] * KW * 1
+                            ) / KP
+
+                    elif (
+                        kj == self.elements_circumferential - 1
+                        and ki == self.elements_axial - 1
+                    ):  # Corner superior right
+
+                        Mat_coef_st[k, k] = CP - CN - CE
+                        Mat_coef_st[k, k - 1] = CW
+                        Mat_coef_st[k, k - self.elements_circumferential] = CS
+
+                        if p[k] > 0:
+
+                            self.theta_vol[k] = 1
+                            B[k] = -KP * self.theta_vol[k] - KW * self.theta_vol[k - 1]
+                            pp = np.zeros((nk - 1, 1))
+                            pp = np.delete(p, k)
+                            C = np.zeros((1, nk))
+                            C = Mat_coef_st[k, :]
+                            C = np.delete(C, k)
+                            p[k] = (B[k] - np.matmul(C, pp)) / Mat_coef_st[k, k]
+
+                        else:
+
+                            p[k] = 0
+                            B_theta[k] = -np.matmul(Mat_coef_st[k, :], p)
+                            self.theta_vol[k] = (
+                                B_theta[k] - KW * self.theta_vol[k - 1]
+                            ) / KP
+
+                    k = k + 1
+                    kj = kj + 1
+
+                kj = 0
+                ki = ki + 1
+
+            self.erro = np.linalg.norm(p - p_old) + np.linalg.norm(
+                self.theta_vol - theta_vol_old
+            )
+
+        cont = 0
+
+        for i in np.arange(self.elements_axial):
+            for j in np.arange(self.elements_circumferential):
+
+                self.P[i, j, n_p] = p[cont]
+
+                self.Theta_vol[i, j, n_p] = self.theta_vol[cont]
+
+                cont = cont + 1
+
+                if self.P[i, j, n_p] < 0:
+                    self.P[i, j, n_p] = 0
+
+        # Dimensional pressure fied
+
+        self.Pdim = (
+            self.P * self.reference_viscosity * self.speed * (self.journal_radius**2)
+        ) / (self.radial_clearance**2)
+
+        return self.P
+
     def _forces(self, initial_guess, y0, xpt0, ypt0):
         """Calculates the forces in Y and X direction.
 
@@ -316,6 +1011,7 @@ class THDCylindrical(BearingElement):
         Fhy : float
             Force in Y direction. The unit is newton.
         """
+
         if y0 is None and xpt0 is None and ypt0 is None:
             self.initial_guess = initial_guess
 
@@ -360,27 +1056,20 @@ class THDCylindrical(BearingElement):
         )
 
         Ytheta = [
-            np.linspace(t1, t2, self.elements_circumferential)
+            np.linspace(
+                t1 + self.dtheta / 2,
+                t2 - self.dtheta / 2,
+                self.elements_circumferential,
+            )
             for t1, t2 in zip(self.thetaI, self.thetaF)
         ]
 
-        self.pad_ct = [ang for ang in range(0, 360, int(360 / self.n_pad))]
-
-        self.thetaI = np.radians(
-            [pad + (180 / self.n_pad) - (self.betha_s_dg / 2) for pad in self.pad_ct]
-        )
-
-        self.thetaF = np.radians(
-            [pad + (180 / self.n_pad) + (self.betha_s_dg / 2) for pad in self.pad_ct]
-        )
-
-        Ytheta = [
-            np.linspace(t1, t2, self.elements_circumferential)
-            for t1, t2 in zip(self.thetaI, self.thetaF)
-        ]
+        self.theta_vol_groove = 0.8 * np.ones(self.n_pad)
 
-        while (T_mist[0] - T_conv) >= 1e-2:
+        T_end = np.ones(self.n_pad)
 
+        while (T_mist[0] - T_conv) >= 0.5:
+            nk = (self.elements_axial) * (self.elements_circumferential)
             self.P = np.zeros(
                 (self.elements_axial, self.elements_circumferential, self.n_pad)
             )
@@ -400,7 +1089,9 @@ class THDCylindrical(BearingElement):
                 * 1.2
             )
 
-            T_conv = T_mist[0]
+            self.Theta_vol = np.zeros(
+                (self.elements_axial, self.elements_circumferential, self.n_pad)
+            )
 
             mu_new = 1.1 * np.ones(
                 (self.elements_axial, self.elements_circumferential, self.n_pad)
@@ -409,273 +1100,69 @@ class THDCylindrical(BearingElement):
                 (self.elements_axial, self.elements_circumferential, self.n_pad)
             )
 
-            PP = np.zeros(((self.elements_axial), (2 * self.elements_circumferential)))
-
-            nk = (self.elements_axial) * (self.elements_circumferential)
-
-            Mat_coef = np.zeros((nk, nk))
-            Mat_coef_T = np.zeros((nk, nk))
-            b = np.zeros((nk, 1))
-            b_T = np.zeros((nk, 1))
-
-            for n_p in np.arange(self.n_pad):
-
-                T_ref = T_mist[n_p - 1]
-
-                # Temperature convergence while
-
-                while (
-                    np.linalg.norm(T_new[:, :, n_p] - T[:, :, n_p])
-                    / np.linalg.norm(T[:, :, n_p])
-                    >= 1e-3
-                ):
-
-                    T_ref = T_mist[n_p - 1]
-
-                    mu = mu_new
-                    self.mu_l = mu_new
-
-                    T[:, :, n_p] = T_new[:, :, n_p]
-
-                    ki = 0
-                    kj = 0
-                    k = 0
-
-                    # Solution of pressure field initialization
-
-                    for ii in np.arange((self.Z_I + 0.5 * self.dZ), self.Z_F, self.dZ):
-                        for jj in np.arange(
-                            self.thetaI[n_p] + (self.dtheta / 2),
-                            self.thetaF[n_p],
-                            self.dtheta,
-                        ):
-
-                            hP = 1 - self.X * np.cos(jj) - self.Y * np.sin(jj)
-                            he = (
-                                1
-                                - self.X * np.cos(jj + 0.5 * self.dtheta)
-                                - self.Y * np.sin(jj + 0.5 * self.dtheta)
-                            )
-                            hw = (
-                                1
-                                - self.X * np.cos(jj - 0.5 * self.dtheta)
-                                - self.Y * np.sin(jj - 0.5 * self.dtheta)
-                            )
-                            hn = hP
-                            hs = hn
-
-                            if kj == 0 and ki == 0:
-                                MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
-                                MU_w = mu[ki, kj]
-                                MU_s = mu[ki, kj]
-                                MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
-
-                            if kj == 0 and ki > 0 and ki < self.elements_axial - 1:
-                                MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
-                                MU_w = mu[ki, kj]
-                                MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
-                                MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
-
-                            if kj == 0 and ki == self.elements_axial - 1:
-                                MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
-                                MU_w = mu[ki, kj]
-                                MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
-                                MU_n = mu[ki, kj]
-
-                            if (
-                                ki == 0
-                                and kj > 0
-                                and kj < self.elements_circumferential - 1
-                            ):
-                                MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
-                                MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
-                                MU_s = mu[ki, kj]
-                                MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+            T_conv = T_mist[0]
 
-                            if (
-                                kj > 0
-                                and kj < self.elements_circumferential - 1
-                                and ki > 0
-                                and ki < self.elements_axial - 1
-                            ):
-                                MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
-                                MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
-                                MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
-                                MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+            self.H = np.ones((self.elements_circumferential, self.n_pad))
 
-                            if (
-                                ki == self.elements_axial - 1
-                                and kj > 0
-                                and kj < self.elements_circumferential - 1
-                            ):
-                                MU_e = 0.5 * (mu[ki, kj] + mu[ki, kj + 1])
-                                MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
-                                MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
-                                MU_n = mu[ki, kj]
+            U = 0.5 * np.ones(
+                (self.elements_axial, self.elements_circumferential, self.n_pad)
+            )
 
-                            if ki == 0 and kj == self.elements_circumferential - 1:
-                                MU_e = mu[ki, kj]
-                                MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
-                                MU_s = mu[ki, kj]
-                                MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+            self.V = np.zeros(
+                (self.elements_axial, self.elements_circumferential, self.n_pad)
+            )
 
-                            if (
-                                kj == self.elements_circumferential - 1
-                                and ki > 0
-                                and ki < self.elements_axial - 1
-                            ):
-                                MU_e = mu[ki, kj]
-                                MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
-                                MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
-                                MU_n = 0.5 * (mu[ki, kj] + mu[ki + 1, kj])
+            self.Qedim = np.ones(self.n_pad)
 
-                            if (
-                                kj == self.elements_circumferential - 1
-                                and ki == self.elements_axial - 1
-                            ):
-                                MU_e = mu[ki, kj]
-                                MU_w = 0.5 * (mu[ki, kj] + mu[ki, kj - 1])
-                                MU_s = 0.5 * (mu[ki, kj] + mu[ki - 1, kj])
-                                MU_n = mu[ki, kj]
+            self.Qsdim = np.ones(self.n_pad)
 
-                            CE = (self.dZ * he**3) / (
-                                12 * MU_e[n_p] * self.dY * self.betha_s**2
-                            )
-                            CW = (self.dZ * hw**3) / (
-                                12 * MU_w[n_p] * self.dY * self.betha_s**2
-                            )
-                            CN = (self.dY * (self.journal_radius**2) * hn**3) / (
-                                12 * MU_n[n_p] * self.dZ * self.axial_length**2
-                            )
-                            CS = (self.dY * (self.journal_radius**2) * hs**3) / (
-                                12 * MU_s[n_p] * self.dZ * self.axial_length**2
-                            )
-                            CP = -(CE + CW + CN + CS)
+            self.Qldim = np.ones(self.n_pad)
 
-                            B = (self.dZ / (2 * self.betha_s)) * (he - hw) - (
-                                (self.Ypt * np.cos(jj) + self.Xpt * np.sin(jj))
-                                * self.dy
-                                * self.dZ
-                            )
+            b_T = np.zeros((nk, 1))
 
-                            k = k + 1
-                            b[k - 1, 0] = B
+            b_P = np.zeros((nk, 1))
 
-                            if ki == 0 and kj == 0:
-                                Mat_coef[k - 1, k - 1] = CP - CS - CW
-                                Mat_coef[k - 1, k] = CE
-                                Mat_coef[
-                                    k - 1, k + self.elements_circumferential - 1
-                                ] = CN
+            Mat_coef = np.zeros((nk, nk))  # Coeficients matrix
 
-                            elif kj == 0 and ki > 0 and ki < self.elements_axial - 1:
-                                Mat_coef[k - 1, k - 1] = CP - CW
-                                Mat_coef[k - 1, k] = CE
-                                Mat_coef[
-                                    k - 1, k - self.elements_circumferential - 1
-                                ] = CS
-                                Mat_coef[
-                                    k - 1, k + self.elements_circumferential - 1
-                                ] = CN
+            B = np.zeros((nk, 1))  # Termo fonte for pressure
 
-                            elif kj == 0 and ki == self.elements_axial - 1:
-                                Mat_coef[k - 1, k - 1] = CP - CN - CW
-                                Mat_coef[k - 1, k] = CE
-                                Mat_coef[
-                                    k - 1, k - self.elements_circumferential - 1
-                                ] = CS
+            for n_p in np.arange(self.n_pad):
 
-                            elif ki == 0 and kj > 0 and kj < self.n_y - 1:
-                                Mat_coef[k - 1, k - 1] = CP - CS
-                                Mat_coef[k - 1, k] = CE
-                                Mat_coef[k - 1, k - 2] = CW
-                                Mat_coef[
-                                    k - 1, k + self.elements_circumferential - 1
-                                ] = CN
+                T_ref = T_mist[n_p]
 
-                            elif (
-                                ki > 0
-                                and ki < self.elements_axial - 1
-                                and kj > 0
-                                and kj < self.n_y - 1
-                            ):
-                                Mat_coef[k - 1, k - 1] = CP
-                                Mat_coef[k - 1, k - 2] = CW
-                                Mat_coef[
-                                    k - 1, k - self.elements_circumferential - 1
-                                ] = CS
-                                Mat_coef[
-                                    k - 1, k + self.elements_circumferential - 1
-                                ] = CN
-                                Mat_coef[k - 1, k] = CE
+                while (
+                    np.linalg.norm(T_new[:, :, n_p] - T[:, :, n_p])
+                    / np.linalg.norm(T[:, :, n_p])
+                    >= 0.01
+                ):
 
-                            elif (
-                                ki == self.elements_axial - 1
-                                and kj > 0
-                                and kj < self.n_y - 1
-                            ):
-                                Mat_coef[k - 1, k - 1] = CP - CN
-                                Mat_coef[k - 1, k] = CE
-                                Mat_coef[k - 1, k - 2] = CW
-                                Mat_coef[
-                                    k - 1, k - self.elements_circumferential - 1
-                                ] = CS
+                    T_ref = T_mist[n_p]
 
-                            elif ki == 0 and kj == self.n_y - 1:
-                                Mat_coef[k - 1, k - 1] = CP - CE - CS
-                                Mat_coef[k - 1, k - 2] = CW
-                                Mat_coef[
-                                    k - 1, k + self.elements_circumferential - 1
-                                ] = CN
+                    mu = mu_new
 
-                            elif (
-                                kj == self.n_y - 1
-                                and ki > 0
-                                and ki < self.elements_axial - 1
-                            ):
-                                Mat_coef[k - 1, k - 1] = CP - CE
-                                Mat_coef[k - 1, k - 2] = CW
-                                Mat_coef[
-                                    k - 1, k - self.elements_circumferential - 1
-                                ] = CS
-                                Mat_coef[
-                                    k - 1, k + self.elements_circumferential - 1
-                                ] = CN
+                    self.mu_l = mu_new
 
-                            elif ki == self.elements_axial - 1 and kj == self.n_y - 1:
-                                Mat_coef[k - 1, k - 1] = CP - CE - CN
-                                Mat_coef[k - 1, k - 2] = CW
-                                Mat_coef[
-                                    k - 1, k - self.elements_circumferential - 1
-                                ] = CS
+                    T[:, :, n_p] = T_new[:, :, n_p]
 
-                            kj = kj + 1
+                    self.erro = 1
 
-                        kj = 0
-                        ki = ki + 1
+                    p_old = np.zeros((nk, 1))
 
-                    # Solution of pressure field end
+                    self.theta_vol = np.zeros((nk, 1))  # Theta volumetric vector
 
-                    p = np.linalg.solve(Mat_coef, b)
-                    cont = 0
+                    Mat_coef_st = np.zeros((nk, nk))  # Coeficients matrix
 
-                    for i in np.arange(self.elements_axial):
-                        for j in np.arange(self.elements_circumferential):
+                    Mat_coef_T = np.zeros((nk, nk))
 
-                            self.P[i, j, n_p] = p[cont]
-                            cont = cont + 1
+                    p = np.ones((nk, 1))  # Pressure vector
 
-                            if self.P[i, j, n_p] < 0:
-                                self.P[i, j, n_p] = 0
+                    B_theta = np.zeros((nk, 1))
 
-                    # Dimensional pressure fied
+                    if self.operating_type == "flooded":
+                        self._flooded(n_p, Mat_coef, b_P, mu)
 
-                    self.Pdim = (
-                        self.P
-                        * self.reference_viscosity
-                        * self.speed
-                        * (self.journal_radius**2)
-                    ) / (self.radial_clearance**2)
+                    elif self.operating_type == "starvation":
+                        self._starvation(n_p, Mat_coef_st, mu, p_old, p, B, B_theta, nk)
 
                     ki = 0
                     kj = 0
@@ -787,24 +1274,37 @@ class THDCylindrical(BearingElement):
                                 )
 
                             HP = 1 - self.X * np.cos(jj) - self.Y * np.sin(jj)
-                            hpt = -self.Ypt * np.cos(jj) + self.Xpt * np.sin(jj)
+                            hpt = -self.Xpt * np.cos(jj) - self.Ypt * np.sin(jj)
+                            self.H[kj, n_p] = HP
 
                             mu_p = mu[ki, kj, n_p]
 
-                            Reyn[ki, kj, n_p] = (
-                                self.rho
-                                * self.speed
-                                * self.journal_radius
-                                * (HP / self.axial_length)
-                                * self.radial_clearance
-                                / (self.reference_viscosity)
-                            )
+                            if self.operating_type == "starvation":
+                                Reyn[ki, kj, n_p] = (
+                                    self.Theta_vol[ki, kj, n_p]
+                                    * self.rho
+                                    * self.speed
+                                    * self.journal_radius
+                                    * (HP / self.axial_length)
+                                    * self.radial_clearance
+                                    / (self.reference_viscosity * mu_p)
+                                )
+
+                            else:
+                                Reyn[ki, kj, n_p] = (
+                                    self.rho
+                                    * self.speed
+                                    * self.journal_radius
+                                    * (HP / self.axial_length)
+                                    * self.radial_clearance
+                                    / (self.reference_viscosity * mu_p)
+                                )
 
                             if Reyn[ki, kj, n_p] <= 500:
 
                                 self.delta_turb = 0
 
-                            elif Reyn[ki, kj, n_p] > 400 and Reyn[ki, kj, n_p] <= 1000:
+                            elif Reyn[ki, kj, n_p] > 500 and Reyn[ki, kj, n_p] <= 1000:
 
                                 self.delta_turb = 1 - (
                                     (1000 - Reyn[ki, kj, n_p]) / 500
@@ -839,6 +1339,13 @@ class THDCylindrical(BearingElement):
 
                             mi_t = mu_turb[ki, kj, n_p]
 
+                            U[ki, kj, n_p] = (
+                                -(HP**2)
+                                / (12 * mi_t * self.betha_s)
+                                * dPdy[ki, kj, n_p]
+                                + 1 / 2
+                            )
+
                             AE = -(self.k_t * HP * self.dZ) / (
                                 self.rho
                                 * self.Cp
@@ -1015,7 +1522,11 @@ class THDCylindrical(BearingElement):
                                     T_ref / self.reference_temperature
                                 )
 
-                            elif ki == 0 and kj > 0 and kj < self.n_y - 1:
+                            elif (
+                                ki == 0
+                                and kj > 0
+                                and kj < self.elements_circumferential - 1
+                            ):
                                 Mat_coef_T[k - 1, k - 1] = AP + AS
                                 Mat_coef_T[k - 1, k] = AE
                                 Mat_coef_T[k - 1, k - 2] = AW
@@ -1027,7 +1538,7 @@ class THDCylindrical(BearingElement):
                                 ki > 0
                                 and ki < self.elements_axial - 1
                                 and kj > 0
-                                and kj < self.n_y - 1
+                                and kj < self.elements_circumferential - 1
                             ):
                                 Mat_coef_T[k - 1, k - 1] = AP
                                 Mat_coef_T[k - 1, k - 2] = AW
@@ -1042,7 +1553,7 @@ class THDCylindrical(BearingElement):
                             elif (
                                 ki == self.elements_axial - 1
                                 and kj > 0
-                                and kj < self.n_y - 1
+                                and kj < self.elements_circumferential - 1
                             ):
                                 Mat_coef_T[k - 1, k - 1] = AP + AN
                                 Mat_coef_T[k - 1, k] = AE
@@ -1051,7 +1562,7 @@ class THDCylindrical(BearingElement):
                                     k - 1, k - self.elements_circumferential - 1
                                 ] = AS
 
-                            elif ki == 0 and kj == self.n_y - 1:
+                            elif ki == 0 and kj == self.elements_circumferential - 1:
                                 Mat_coef_T[k - 1, k - 1] = AP + AE + AS
                                 Mat_coef_T[k - 1, k - 2] = AW
                                 Mat_coef_T[
@@ -1059,7 +1570,7 @@ class THDCylindrical(BearingElement):
                                 ] = AN
 
                             elif (
-                                kj == self.n_y - 1
+                                kj == self.elements_circumferential - 1
                                 and ki > 0
                                 and ki < self.elements_axial - 1
                             ):
@@ -1072,7 +1583,10 @@ class THDCylindrical(BearingElement):
                                     k - 1, k + self.elements_circumferential - 1
                                 ] = AN
 
-                            elif ki == self.elements_axial - 1 and kj == self.n_y - 1:
+                            elif (
+                                ki == self.elements_axial - 1
+                                and kj == self.elements_circumferential - 1
+                            ):
                                 Mat_coef_T[k - 1, k - 1] = AP + AE + AN
                                 Mat_coef_T[k - 1, k - 2] = AW
                                 Mat_coef_T[
@@ -1097,28 +1611,76 @@ class THDCylindrical(BearingElement):
 
                     Tdim = T_new * self.reference_temperature
 
-                    T_end = np.sum(Tdim[:, -1, n_p]) / self.elements_axial
+                    T_end[n_p] = np.sum(Tdim[:, -1, n_p]) / self.elements_axial
 
-                    T_mist[n_p] = (
-                        self.fat_mixt[n_p] * self.reference_temperature
-                        + (1 - self.fat_mixt[n_p]) * T_end
+                    if self.operating_type == "flooded":
+
+                        T_mist[n_p - 1] = (
+                            self.fat_mixt[n_p] * self.reference_temperature
+                            + (1 - self.fat_mixt[n_p]) * T_end[n_p]
+                        )
+
+                    mu_new[:, :, n_p] = (
+                        self.a * (Tdim[:, :, n_p]) ** self.b
+                    ) / self.reference_viscosity
+
+            if self.operating_type == "starvation":
+                for n_p in np.arange(self.n_pad):
+
+                    self.Qedim[n_p] = (
+                        self.radial_clearance
+                        * self.H[0, n_p]
+                        * self.speed
+                        * self.journal_radius
+                        * self.axial_length
+                        * self.Theta_vol[0, 0, n_p]
+                        * (np.mean(U[:, 0, n_p]))
                     )
 
-                    for i in np.arange(self.elements_axial):
-                        for j in np.arange(self.elements_circumferential):
+                    self.Qsdim[n_p] = (
+                        self.radial_clearance
+                        * self.H[-1, n_p]
+                        * self.speed
+                        * self.journal_radius
+                        * self.axial_length
+                        * self.Theta_vol[0, -1, n_p]
+                        * (np.mean(U[:, -1, n_p]))
+                    )
+
+                geometry_factor = np.ones(self.n_pad)
+
+                for n_p in np.arange(self.n_pad):
+                    geometry_factor[n_p] = (self.Qedim[n_p] + self.Qsdim[n_p - 1]) / (
+                        np.sum(self.Qedim) + np.sum(self.Qsdim)
+                    )
+
+                for n_p in np.arange(self.n_pad):
+
+                    T_mist[n_p] = (
+                        (self.Qsdim[n_p - 1] * T_end[n_p - 1])
+                        + (
+                            self.reference_temperature
+                            * geometry_factor[n_p]
+                            * self.oil_flow
+                        )
+                    ) / (geometry_factor[n_p] * self.oil_flow + self.Qsdim[n_p - 1])
+
+                    self.theta_vol_groove[n_p] = (
+                        0.8
+                        * (geometry_factor[n_p] * self.oil_flow + self.Qsdim[n_p - 1])
+                        / self.Qedim[n_p]
+                    )
 
-                            mu_new[i, j, n_p] = (
-                                self.a * (Tdim[i, j, n_p]) ** self.b
-                            ) / self.reference_viscosity
+                    if self.theta_vol_groove[n_p] > 1:
+                        self.theta_vol_groove[n_p] = 1
 
-        PP = np.zeros(
-            ((self.elements_axial), (self.n_pad * self.elements_circumferential))
+        PPlot = np.zeros(
+            (self.elements_axial, self.elements_circumferential * self.n_pad)
         )
 
-        i = 0
         for i in range(self.elements_axial):
 
-            PP[i] = self.Pdim[i, :, :].ravel("F")
+            PPlot[i] = self.Pdim[i, :, :].ravel("F")
 
         Ytheta = np.array(Ytheta)
         Ytheta = Ytheta.flatten()
@@ -1128,21 +1690,17 @@ class THDCylindrical(BearingElement):
         auxF[0, :] = np.cos(Ytheta)
         auxF[1, :] = np.sin(Ytheta)
 
-        dA = self.dy * self.dz
-
-        auxP = PP * dA
-
-        vector_auxF_x = auxF[0, :]
-        vector_auxF_y = auxF[1, :]
+        fx1 = np.trapz(PPlot * auxF[0, :], self.journal_radius * Ytheta)
+        Fhx = -np.trapz(fx1, self.axial_length * self.Z[1 : self.elements_axial + 1])
 
-        auxFx = auxP * vector_auxF_x
-        auxFy = auxP * vector_auxF_y
+        fy1 = np.trapz(PPlot * auxF[1, :], self.journal_radius * Ytheta)
+        Fhy = -np.trapz(fy1, self.axial_length * self.Z[1 : self.elements_axial + 1])
+        F1 = Fhx
+        F2 = Fhy
 
-        fxj = -np.sum(auxFx)
-        fyj = -np.sum(auxFy)
+        Fhx = F1
+        Fhy = F2
 
-        Fhx = fxj
-        Fhy = fyj
         self.Fhx = Fhx
         self.Fhy = Fhy
         return Fhx, Fhy
@@ -1151,8 +1709,8 @@ class THDCylindrical(BearingElement):
         """This method runs the optimization to find the equilibrium position of
         the rotor's center.
 
-
         """
+
         args = self.print_progress
         t1 = time.time()
         res = minimize(
@@ -1160,8 +1718,8 @@ class THDCylindrical(BearingElement):
             self.initial_guess,
             args,
             method="Nelder-Mead",
-            tol=10e-3,
-            options={"maxiter": 1000},
+            tol=10e-2,
+            options={"maxiter": 1e4},
         )
         self.equilibrium_pos = res.x
         t2 = time.time()
@@ -1200,7 +1758,7 @@ class THDCylindrical(BearingElement):
         xdata = [T_muI, T_muF]  # changed boundary conditions to avoid division by ]
         ydata = [mu_I, mu_F]
 
-        popt, pcov = curve_fit(viscosity, xdata, ydata, p0=(6.0, -1.0))
+        popt, pcov = curve_fit(viscosity, xdata, ydata, p0=(6.0, 1.0))
         a, b = popt
 
         return a, b
@@ -1296,20 +1854,20 @@ class THDCylindrical(BearingElement):
             (Auinitial_guess3[1] - Auinitial_guess4[1]) / (epiy / self.radial_clearance)
         )
 
-        Cxx = -self.sommerfeld(Auinitial_guess5[0], Auinitial_guess6[0]) * (
-            (Auinitial_guess6[0] - Auinitial_guess5[0])
+        Cxx = -self.sommerfeld(Auinitial_guess5[0], Auinitial_guess6[1]) * (
+            (Auinitial_guess5[0] - Auinitial_guess6[0])
             / (epixpt / self.radial_clearance / self.speed)
         )
-        Cxy = -self.sommerfeld(Auinitial_guess7[0], Auinitial_guess8[0]) * (
-            (Auinitial_guess8[0] - Auinitial_guess7[0])
+        Cxy = -self.sommerfeld(Auinitial_guess7[0], Auinitial_guess8[1]) * (
+            (Auinitial_guess7[0] - Auinitial_guess8[0])
             / (epiypt / self.radial_clearance / self.speed)
         )
-        Cyx = -self.sommerfeld(Auinitial_guess5[1], Auinitial_guess6[1]) * (
-            (Auinitial_guess6[1] - Auinitial_guess5[1])
+        Cyx = -self.sommerfeld(Auinitial_guess5[0], Auinitial_guess6[1]) * (
+            (Auinitial_guess5[1] - Auinitial_guess6[1])
             / (epixpt / self.radial_clearance / self.speed)
         )
-        Cyy = -self.sommerfeld(Auinitial_guess7[1], Auinitial_guess8[1]) * (
-            (Auinitial_guess8[1] - Auinitial_guess7[1])
+        Cyy = -self.sommerfeld(Auinitial_guess7[0], Auinitial_guess8[1]) * (
+            (Auinitial_guess7[1] - Auinitial_guess8[1])
             / (epiypt / self.radial_clearance / self.speed)
         )
 
@@ -1354,8 +1912,6 @@ class THDCylindrical(BearingElement):
         order perturbation solution) is aplied. The four stiffness coefficients,
         and the four damping coefficients is obtained by integration of the pressure
         field.
-
-
         """
 
         p = self.P
@@ -1857,6 +2413,7 @@ class THDCylindrical(BearingElement):
         Score coefficient.
 
         """
+
         Fhx, Fhy = self._forces(x, None, None, None)
         score = np.sqrt(
             ((self.load_x_direction + Fhx) ** 2) + ((self.load_y_direction + Fhy) ** 2)
@@ -1887,6 +2444,7 @@ class THDCylindrical(BearingElement):
         Ss : float
             Sommerfeld number.
         """
+
         if self.sommerfeld_type == 1:
             S = (
                 self.reference_viscosity
@@ -1928,16 +2486,14 @@ def cylindrical_bearing_example():
     """
 
     bearing = THDCylindrical(
-        axial_length=0.263144,
-        journal_radius=0.2,
-        radial_clearance=1.95e-4,
+        axial_length=(0.263144),
+        journal_radius=(0.2),
+        radial_clearance=(1.95e-4),
         elements_circumferential=11,
         elements_axial=3,
-        n_y=None,
         n_pad=2,
         pad_arc_length=176,
         reference_temperature=50,
-        reference_viscosity=0.02,
         speed=Q_([900], "RPM"),
         load_x_direction=0,
         load_y_direction=-112814.91,
@@ -1947,6 +2503,9 @@ def cylindrical_bearing_example():
         sommerfeld_type=2,
         initial_guess=[0.1, -0.1],
         method="perturbation",
+        operating_type="flooded",
+        injection_pressure=0,
+        oil_flow=37.86,
         show_coef=False,
         print_result=False,
         print_progress=False,


Test patch:
diff --git a/ross/tests/test_cylindrical.py b/ross/tests/test_cylindrical.py
index b30d74b3..49da74b2 100644
--- a/ross/tests/test_cylindrical.py
+++ b/ross/tests/test_cylindrical.py
@@ -18,13 +18,11 @@ def cylindrical():
         axial_length=L,
         journal_radius=0.2,
         radial_clearance=1.95e-4,
-        elements_circumferential=4,
+        elements_circumferential=11,
         elements_axial=3,
-        n_y=None,
         n_pad=2,
         pad_arc_length=176,
         reference_temperature=50,
-        reference_viscosity=0.02,
         speed=speed,
         load_x_direction=0,
         load_y_direction=-112814.91,
@@ -34,8 +32,11 @@ def cylindrical():
         sommerfeld_type=2,
         initial_guess=[0.1, -0.1],
         method="perturbation",
+        operating_type="flooded",
+        injection_pressure=0,
+        oil_flow=37.86,
         show_coef=False,
-        print_result=True,
+        print_result=False,
         print_progress=False,
         print_time=False,
     )
@@ -52,12 +53,8 @@ def test_cylindrical_parameters(cylindrical):
 
 
 def test_cylindrical_equilibrium_pos(cylindrical):
-    assert math.isclose(
-        cylindrical.equilibrium_pos[0], 0.6063220830888834, rel_tol=0.01
-    )
-    assert math.isclose(
-        cylindrical.equilibrium_pos[1], -0.8080463672910241, rel_tol=0.01
-    )
+    assert math.isclose(cylindrical.equilibrium_pos[0], 0.6873316, rel_tol=0.01)
+    assert math.isclose(cylindrical.equilibrium_pos[1], -0.79393636, rel_tol=0.01)
 
 
 def test_cylindrical_coefficients(cylindrical):
@@ -71,11 +68,11 @@ def test_cylindrical_coefficients(cylindrical):
     cyx = coefs[1][2]
     cyy = coefs[1][3]
 
-    assert math.isclose(kxx, 1154092585.2914207, rel_tol=0.0001)
-    assert math.isclose(kxy, -78688969.0436872, rel_tol=0.0001)
-    assert math.isclose(kyx, -1399698008.045516, rel_tol=0.0001)
-    assert math.isclose(kyy, 1333930910.0158768, rel_tol=0.0001)
-    assert math.isclose(cxx, 5853947611.807541, rel_tol=0.0001)
-    assert math.isclose(cxy, -6125302956.996735, rel_tol=0.0001)
-    assert math.isclose(cyx, -16878755.32663667, rel_tol=0.0001)
-    assert math.isclose(cyy, 4366663.510970946, rel_tol=0.0001)
+    assert math.isclose(kxx, 1080948741.8512235, rel_tol=0.0001)
+    assert math.isclose(kxy, 339258572.34310913, rel_tol=0.0001)
+    assert math.isclose(kyx, -1359170639.0567815, rel_tol=0.0001)
+    assert math.isclose(kyy, 1108970752.2456105, rel_tol=0.0001)
+    assert math.isclose(cxx, 9815899.503793057, rel_tol=0.0001)
+    assert math.isclose(cxy, -9963602.922357056, rel_tol=0.0001)
+    assert math.isclose(cyx, -11312462.69772395, rel_tol=0.0001)
+    assert math.isclose(cyy, 27194995.506247465, rel_tol=0.0001)