Merge pull request #1211 from gsabinoo/fix/cylindrical

Correct formulation when using Lund method

Author: raphaeltimbo

docs/references/api.rst

@@ -23,7 +23,7 @@ Elements
     SealElement
     BallBearingElement
     RollerBearingElement
-    THDCylindrical
+    PlainJournal
     BearingFluidFlow
     MagneticBearingElement
     GearElement

docs/user_guide/tutorial_part_1.ipynb

@@ -1215,7 +1215,7 @@
     "- 4. `RollerBearingElement`: A bearing element for roller bearings.\n",
     "- 5. `MagneticBearingElement`: A bearing element for magnetic bearings.\n",
     "    - 5.1. `param_to_coef`: A bearing element for magnetic bearings from electromagnetic parameters\n",
-    "- 6. `THDCylindrical`: A cylindrical bearing element based on the pressure and temperature field in oil film.\n",
+    "- 6. `PlainJournal`: A cylindrical bearing element based on the pressure and temperature field in oil film.\n",
     "\n",
     "The classes from item 2 to 6 inherits from `BearingElement` class. It means, you can use the same methods and commands, set up to `BearingElement`, in the other classes."
    ]
@@ -1870,7 +1870,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### 4.6 THDCylindrical Class\n",
+    "### 4.6 PlainJournal Class\n",
     "\n",
     "This class computes the pressure and temperature fields within the oil film of a cylindrical bearing while also determining the stiffness and damping coefficients. These calculations are based on a wide range of inputs, including bearing geometry, operating conditions, fluid properties, turbulence modeling, and mesh discretization.\n",
     "\n",
@@ -1913,7 +1913,7 @@
     {
      "data": {
       "text/plain": [
-       "THDCylindrical(n=3, n_link=None,\n",
+       "PlainJournal(n=3, n_link=None,\n",
        " kxx=[2.49767453e+09], kxy=[-3.14056282e+09],\n",
        " kyx=[1.00992315e+08], kyy=[7.8393767e+08],\n",
        " kzz=[0], cxx=[36950674.61976394],\n",
@@ -1931,7 +1931,7 @@
     }
    ],
    "source": [
-    "cybearing = rs.THDCylindrical(\n",
+    "cybearing = rs.PlainJournal(\n",
     "    axial_length=0.263144,\n",
     "    journal_radius=0.2,\n",
     "    radial_clearance=1.95e-4,\n",

ross/__init__.py

@@ -19,7 +19,7 @@
 from .units import Q_
 from .utils import get_data_from_figure, visualize_matrix
 from ross.bearings.lubricants import lubricants_dict
-from ross.bearings.cylindrical import *
+from ross.bearings.plain_journal import *
 from ross.seals.labyrinth_seal import *
 from ross.bearings.thrust_pad import *
 from ross.bearings.tilting_pad import *

ross/bearings/cylindrical.py renamed to ross/bearings/plain_journal.py

@@ -12,9 +12,9 @@
 from ross.bearings.lubricants import lubricants_dict
 
 
-class THDCylindrical(BearingElement):
+class PlainJournal(BearingElement):
     """This class calculates the pressure and temperature field in oil film of
-    a cylindrical bearing. It is also possible to obtain the stiffness and
+    a plain journal bearing. It is also possible to obtain the stiffness and
     damping coefficients.
     The basic references for the code is found in :cite:t:`barbosa2018`, :cite:t:`daniel2012` and :cite:t:`nicoletti1999`.
 
@@ -46,6 +46,9 @@ class THDCylindrical(BearingElement):
         Choose the method to calculate the dynamics coefficients. Options are:
         - 'lund'
         - 'perturbation'
+    model_type : str, optional
+        Type of model to be used. Options:
+        - 'thermo_hydro_dynamic': Thermo-Hydro-Dynamic model
     print_progress : bool
         Set it True to print the score and forces on each iteration.
         False by default.
@@ -111,7 +114,7 @@ class THDCylindrical(BearingElement):
 
     Returns
     -------
-    A THDCylindrical object.
+    A PlainJournal object.
 
     References
     ----------
@@ -138,8 +141,8 @@ class THDCylindrical(BearingElement):
 
     Example
     --------
-    >>> from ross.bearings.cylindrical import THDCylindrical
-    >>> bearing = THDCylindrical(
+    >>> from ross.bearings.plain_journal import PlainJournal
+    >>> bearing = PlainJournal(
     ...    n=3,
     ...    axial_length=0.263144,
     ...    journal_radius=0.2,
@@ -193,6 +196,7 @@ def __init__(
         sommerfeld_type=2,
         initial_guess=[0.1, -0.1],
         method="perturbation",
+        model_type="thermo_hydro_dynamic",
         operating_type="flooded",
         oil_supply_pressure=None,
         oil_flow_v=None,
@@ -220,6 +224,7 @@ def __init__(
         self.initial_guess = initial_guess
         self.method = method
         self.operating_type = operating_type
+        self.model_type = model_type
         self.oil_supply_pressure = oil_supply_pressure
         self.oil_flow_v = oil_flow_v
         self.show_coeffs = show_coeffs
@@ -234,14 +239,13 @@ def __init__(
         self.thetaF = self.betha_s
         self.dtheta = (self.thetaF - self.thetaI) / (self.elements_circumferential)
 
-        if self.method == "perturbation":
-            pad_ct = np.arange(0, 360, int(360 / self.n_pad))
-            self.thetaI = np.radians(pad_ct + 180 / self.n_pad - self.betha_s_dg / 2)
-            self.thetaF = np.radians(pad_ct + 180 / self.n_pad + self.betha_s_dg / 2)
-            self.theta_range = [
-                np.arange(start_rad + (self.dtheta / 2), end_rad, self.dtheta)
-                for start_rad, end_rad in zip(self.thetaI, self.thetaF)
-            ]
+        pad_ct = np.arange(0, 360, int(360 / self.n_pad))
+        self.thetaI = np.radians(pad_ct + 180 / self.n_pad - self.betha_s_dg / 2)
+        self.thetaF = np.radians(pad_ct + 180 / self.n_pad + self.betha_s_dg / 2)
+        self.theta_range = [
+            np.arange(start_rad + (self.dtheta / 2), end_rad, self.dtheta)
+            for start_rad, end_rad in zip(self.thetaI, self.thetaF)
+        ]
 
         # Dimensionless discretization variables
         self.dY = 1 / self.elements_circumferential
@@ -289,7 +293,7 @@ def __init__(
         )  # Interpolation function
         self.reference_viscosity = self.interpolate(self.reference_temperature)
 
-        # if self.geometry == "lobe":
+        # Pivot angle for lobe geometry
         self.theta_pivot = np.array([90, 270]) * np.pi / 180
 
         n_freq = np.shape(frequency)[0]
@@ -307,7 +311,8 @@ def __init__(
         for i in range(n_freq):
             speed = frequency[i]
 
-            self.run(speed)
+            if self.model_type == "thermo_hydro_dynamic":
+                self.run_thermo_hydro_dynamic(speed)
 
             coeffs = self.coefficients(speed)
             kxx[i], kxy[i], kyx[i], kyy[i] = coeffs[0]
@@ -608,12 +613,12 @@ def _forces(self, initial_guess, speed, y0=None, xpt0=None, ypt0=None):
                     if self.theta_vol_groove[n_p] > 1:
                         self.theta_vol_groove[n_p] = 1
 
-        self.P = Pdim
-        self.T = Tdim
+        # self.P must receive the adimensional pressure field
+        self.P = Pdim * (self.radial_clearance**2) / (self.reference_viscosity * speed * (self.journal_radius**2))
         self.Theta_vol = Theta_vol
 
-        PPlot = self.P.reshape(self.elements_axial, -1, order="F")
-        TPlot = self.T.reshape(self.elements_axial, -1, order="F")
+        # Reshape dimensional pressure field from (axial, circumferential, pads) to (axial, all_other_dims)
+        PPlot = Pdim.reshape(self.elements_axial, -1, order="F")
 
         Ytheta = np.sort(
             np.linspace(
@@ -632,7 +637,7 @@ def _forces(self, initial_guess, speed, y0=None, xpt0=None, ypt0=None):
 
         return Fhx, Fhy
 
-    def run(self, speed):
+    def run_thermo_hydro_dynamic(self, speed):
         """This method runs the optimization to find the equilibrium position of
         the rotor's center.
         """
@@ -682,7 +687,7 @@ def _get_interp_coeffs(self, T_muI, T_muF, mu_I, mu_F):
         def viscosity(x, a, b):
             return a * (x**b)
 
-        xdata = [T_muI, T_muF]  # changed boundary conditions to avoid division by ]
+        xdata = [T_muI, T_muF]  # changed boundary conditions to avoid division by zero
         ydata = [mu_I, mu_F]
 
         popt, pcov = curve_fit(viscosity, xdata, ydata, p0=(6.0, 1.0))
@@ -713,7 +718,7 @@ def coefficients(self, speed):
         """
 
         if self.equilibrium_pos is None:
-            self.run(speed)
+            self.run_thermo_hydro_dynamic(speed)
             self.coefficients(speed)
         else:
             if self.method == "lund":
@@ -1961,7 +1966,7 @@ def plot_pressure_distribution(self, axial_element_index=None, fig=None, **kwarg
         )
 
         fig.update_layout(
-            title="Cylindrical Bearing",
+            title="Plain Journal Bearing",
             xaxis=dict(
                 showgrid=False,
                 zeroline=False,