Tutorial - Thermo-Hydro-Dynamic (THD) Bearings#
This tutorial covers the thermo-hydro-dynamic (THD) bearing classes available in ROSS. These classes compute dynamic coefficients by solving the Reynolds equation (pressure field) and the energy equation (temperature field) for the lubricant film. They model the physics of the oil film in detail and are the preferred choice for industrial-grade analysis.
Class |
Bearing type |
Model |
|---|---|---|
|
Fixed-geometry journal bearing |
THD, multiple geometries |
|
Tilting-pad journal bearing |
THD, adiabatic or full thermal |
|
Tilting-pad thrust bearing |
THD, axial load capacity |
|
Squeeze film damper |
Analytical short-bearing, 3 geometries |
All four classes inherit from BearingElement and can be placed in a rs.Rotor assembly exactly like the simpler classes described in the general bearings tutorial.
Section 1: PlainJournal#
This tutorial demonstrates how to use the PlainJournal class from ROSS to simulate the thermohydrodynamic behavior of fixed-geometry journal bearings (cylindrical, lemon-bore, elliptical, and lobed) and to determine their dynamic coefficients for different operation conditions.
1.1: Example of application#
First, to run the test case provided as an example in the documentation of the PlainJournal class, use:
from ross.bearings.plain_journal import PlainJournal
from ross.units import Q_
import plotly.io as pio
pio.renderers.default = "notebook"
plainjournal = PlainJournal(
n=3, # Nodal location of the bearing in the FE model of the rotor
axial_length=0.263144, # Axial length of the bearing
journal_radius=0.2, # Journal radius
radial_clearance=1.95e-4, # Bearing radial clearance
elements_circumferential=21, # Number of elements in the circumferential direction per pad
elements_axial=11, # Number of elements in the axial direction per pad
n_pad=2, # Total number of pads
pad_arc_length=176, # Angular length of the pads
preload=0, # Preload factor
geometry="circular", # Bearing geometry type
reference_temperature=50, # Oil supply temperature
frequency=Q_([900, 1200], "RPM"), # Journal rotational speed
fxs_load=0, # External load applied to the shaft in the horizontal direction
fys_load=-112814.91, # External load applied to the shaft in the vertical direction
lubricant="ISOVG32", # Lubricant type
sommerfeld_type=2, # Equation used to calculate the Sommerfeld number
initial_guess=[
0.1,
-0.1,
], # Initial guess to find the equilibrium position of the shaft using an optimization process
method="perturbation", # Method used to calculate the bearing dynamic coefficients
operating_type="flooded", # Lubrication condition in the bearing
groove_factor=[
0.52,
0.48,
], # Mixing factor between the inlet oil and recirculating oil
oil_supply_pressure=0, # Oil supply pressure
)
/home/gsabino/Documents/dev/rosstest/fernandarossi/venvross/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
/home/gsabino/Documents/dev/rosstest/fernandarossi/ross/ross/bearings/plain_journal.py:739: OptimizeWarning: Covariance of the parameters could not be estimated
popt, pcov = curve_fit(viscosity, xdata, ydata, p0=(6.0, 1.0))
As shown in the example above, when method="perturbation" is selected, small perturbations are applied to the rotor around its equilibrium position, with the stiffness and damping coefficients calculated using finite-difference approximations of the derivatives of the bearing forces with respect to rotor displacements and velocities.
On the other hand, when method="lund" is selected, the solution of the Reynolds equation is perturbed, with the stiffness and damping coefficients determined by integrating the dynamic pressure fields over the bearing surface.
In addition, to help the user set up a simulation case using the PlainJournal class, the table below specifies the required parameter for each operating_type option:
Option |
Required parameter |
|---|---|
|
|
|
|
1.2: Bearing representation#
To display a schematic representation of the bearing geometry (number of pads, groove locations), including the direction of the external load applied on the shaft, use:
fig = plainjournal.plot_bearing_representation()
fig.show()
Note: The way the PlainJournal class was built, the first groove is always located at the positive horizontal axis.
1.3: Optimization convergence#
To display in table and graph formats the magnitude of the residual force vector acting on the shaft throughout the iterations of the optimization process for each analyzed speed, use:
plainjournal.show_optimization_convergence(show_plots=True)
=========================================================
OPTIMIZATION CONVERGENCE - 900.0 RPM
=========================================================
+---------------------------+---------------------------+
| Iteration | Residual [N] |
+---------------------------+---------------------------+
| 0 | 1.0381e+05 |
| 1 | 1.0284e+05 |
| 2 | 1.0189e+05 |
| 3 | 9.9408e+04 |
| 4 | 9.6102e+04 |
| 5 | 8.8939e+04 |
| 6 | 7.6378e+04 |
| 7 | 5.3391e+04 |
| 8 | 5.3391e+04 |
| 9 | 5.3391e+04 |
| 10 | 5.2609e+04 |
| 11 | 5.2609e+04 |
| 12 | 5.2609e+04 |
| 13 | 5.1748e+04 |
| 14 | 5.1748e+04 |
| 15 | 5.1675e+04 |
| 16 | 5.0488e+04 |
| 17 | 5.0338e+04 |
| 18 | 4.6676e+04 |
| 19 | 4.6676e+04 |
| 20 | 3.8705e+04 |
| 21 | 3.8705e+04 |
| 22 | 2.7423e+04 |
| 23 | 9.9313e+03 |
| 24 | 9.9313e+03 |
| 25 | 9.9313e+03 |
| 26 | 9.9313e+03 |
| 27 | 9.9287e+03 |
| 28 | 4.5320e+03 |
| 29 | 3.7950e+03 |
| 30 | 3.7950e+03 |
| 31 | 1.1300e+03 |
| 32 | 1.1300e+03 |
| 33 | 1.1300e+03 |
| 34 | 8.1374e+02 |
| 35 | 7.4933e+02 |
| 36 | 3.2156e+02 |
| 37 | 3.2156e+02 |
| 38 | 3.2156e+02 |
| 39 | 1.2052e+02 |
| 40 | 1.2052e+02 |
| 41 | 1.2052e+02 |
| 42 | 5.8065e+01 |
| 43 | 5.8065e+01 |
| 44 | 5.8065e+01 |
| 45 | 5.4652e+01 |
| 46 | 2.7175e+01 |
| 47 | 1.0612e+01 |
| 48 | 1.0612e+01 |
| 49 | 9.6573e+00 |
| 50 | 6.3791e+00 |
| 51 | 2.5117e+00 |
| 52 | 2.5117e+00 |
| 53 | 2.5117e+00 |
| 54 | 1.8404e+00 |
| 55 | 1.4085e+00 |
| 56 | 1.0492e+00 |
+---------------------------+---------------------------+
=========================================================
=========================================================
OPTIMIZATION CONVERGENCE - 1200.0 RPM
=========================================================
+---------------------------+---------------------------+
| Iteration | Residual [N] |
+---------------------------+---------------------------+
| 0 | 9.7770e+03 |
| 1 | 3.5396e+03 |
| 2 | 3.5396e+03 |
| 3 | 3.5396e+03 |
| 4 | 3.5396e+03 |
| 5 | 3.5396e+03 |
| 6 | 2.7468e+02 |
| 7 | 2.7468e+02 |
| 8 | 2.7468e+02 |
| 9 | 2.7468e+02 |
| 10 | 2.7468e+02 |
| 11 | 2.7468e+02 |
| 12 | 2.7468e+02 |
| 13 | 2.7468e+02 |
| 14 | 2.7439e+02 |
| 15 | 9.0084e+01 |
| 16 | 9.0084e+01 |
| 17 | 9.0084e+01 |
| 18 | 6.2551e+01 |
| 19 | 4.2383e+01 |
| 20 | 2.7320e+01 |
| 21 | 1.9681e+01 |
| 22 | 1.3802e+01 |
| 23 | 1.0657e+01 |
| 24 | 8.5002e+00 |
| 25 | 5.2910e+00 |
| 26 | 3.2309e+00 |
| 27 | 2.2083e+00 |
| 28 | 1.8081e+00 |
| 29 | 1.1152e+00 |
| 30 | 1.0161e+00 |
+---------------------------+---------------------------+
=========================================================
1.4: Execution time#
To display the total time spent during the execution of the optimization process to find the equilibrium position of the shaft for all speeds specified as input in the PlainJournal class, use:
plainjournal.show_execution_time()
Execution time: 6.27 seconds
1.5: Pressure and temperature distributions#
To display the pressure and temperature distributions (2D contours and 3D profiles) for the last analyzed speed, use:
plainjournal.plot_results(show_plots=True);
1.6: Circumferential pressure profile#
To display the circumferential pressure profile at the middle plane of the bearing (or in the plane closest to it, depending on the number of elements in the axial direction) for the last analyzed speed, use:
plainjournal.plot_pressure_distribution()
Optionally, to display the circumferential pressure profile at a specific axial plane of the bearing, pass the axial element index as an argument to the plot_pressure_distribution method:
plainjournal.plot_pressure_distribution(2)
1.7: Coefficients comparison#
To display in table format the stiffness and damping coefficients of the bearing for all analyzed speeds, use:
plainjournal.show_coefficients_comparison()
=============================================================================================================================================================================================================================================================
DYNAMIC COEFFICIENTS COMPARISON TABLE
=============================================================================================================================================================================================================================================================
+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+
| Frequency [RPM] | kxx [N/m] | kxy [N/m] | kyx [N/m] | kyy [N/m] | cxx [N*s/m] | cxy [N*s/m] | cyx [N*s/m] | cyy [N*s/m] |
+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+
| 900.0 | 2.3587e+09 | 7.2577e+08 | -3.0873e+09 | 2.5047e+09 | 3.3502e+07 | -3.1875e+07 | -3.8124e+07 | 9.5236e+07 |
| 1200.0 | 2.2310e+09 | 7.5598e+08 | -2.9689e+09 | 2.3502e+09 | 2.3617e+07 | -2.0366e+07 | -2.4895e+07 | 6.7264e+07 |
+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+---------------------------+
=============================================================================================================================================================================================================================================================
1.8: Show results#
To display in table format the optimization results and the bearing dynamic coefficients for each analyzed speed, use:
plainjournal.show_results()
=====================================================================================
PLAIN JOURNAL RESULTS - 900.0 RPM
=====================================================================================
+---------------------------+---------------------------+---------------------------+
| Parameter | Value | Unit |
+---------------------------+---------------------------+---------------------------+
| Operating Speed | 900.0 | RPM |
| Eccentricity Ratio | 6.5259e-01 | - |
| Attitude Angle | -4.3884e+01 | deg |
| Load Fx | 0.0000e+00 | N |
| Load Fy | -1.1281e+05 | N |
| kxx (Stiffness) | 2.3587e+09 | N/m |
| kxy (Stiffness) | 7.2577e+08 | N/m |
| kyx (Stiffness) | -3.0873e+09 | N/m |
| kyy (Stiffness) | 2.5047e+09 | N/m |
| cxx (Damping) | 3.3502e+07 | N*s/m |
| cxy (Damping) | -3.1875e+07 | N*s/m |
| cyx (Damping) | -3.8124e+07 | N*s/m |
| cyy (Damping) | 9.5236e+07 | N*s/m |
| Optimization Success | True | - |
| Function Value | 1.0492e+00 | - |
| Iterations | 58 | - |
| Evaluations | 109 | - |
| Execution Time | 1.4276e+01 | s |
+---------------------------+---------------------------+---------------------------+
=====================================================================================
=====================================================================================
PLAIN JOURNAL RESULTS - 1200.0 RPM
=====================================================================================
+---------------------------+---------------------------+---------------------------+
| Parameter | Value | Unit |
+---------------------------+---------------------------+---------------------------+
| Operating Speed | 1200.0 | RPM |
| Eccentricity Ratio | 6.0528e-01 | - |
| Attitude Angle | -4.1085e+01 | deg |
| Load Fx | 0.0000e+00 | N |
| Load Fy | -1.1281e+05 | N |
| kxx (Stiffness) | 2.2310e+09 | N/m |
| kxy (Stiffness) | 7.5598e+08 | N/m |
| kyx (Stiffness) | -2.9689e+09 | N/m |
| kyy (Stiffness) | 2.3502e+09 | N/m |
| cxx (Damping) | 2.3617e+07 | N*s/m |
| cxy (Damping) | -2.0366e+07 | N*s/m |
| cyx (Damping) | -2.4895e+07 | N*s/m |
| cyy (Damping) | 6.7264e+07 | N*s/m |
| Optimization Success | True | - |
| Function Value | 1.0161e+00 | - |
| Iterations | 32 | - |
| Evaluations | 63 | - |
| Execution Time | 6.2727e+00 | s |
+---------------------------+---------------------------+---------------------------+
=====================================================================================
Optionally, to display the results only for a given speed (e.g., 900 rpm), use:
plainjournal._print_single_frequency_results(plainjournal.frequency[0])
=====================================================================================
PLAIN JOURNAL RESULTS - 900.0 RPM
=====================================================================================
+---------------------------+---------------------------+---------------------------+
| Parameter | Value | Unit |
+---------------------------+---------------------------+---------------------------+
| Operating Speed | 900.0 | RPM |
| Eccentricity Ratio | 6.5259e-01 | - |
| Attitude Angle | -4.3884e+01 | deg |
| Load Fx | 0.0000e+00 | N |
| Load Fy | -1.1281e+05 | N |
| kxx (Stiffness) | 2.3587e+09 | N/m |
| kxy (Stiffness) | 7.2577e+08 | N/m |
| kyx (Stiffness) | -3.0873e+09 | N/m |
| kyy (Stiffness) | 2.5047e+09 | N/m |
| cxx (Damping) | 3.3502e+07 | N*s/m |
| cxy (Damping) | -3.1875e+07 | N*s/m |
| cyx (Damping) | -3.8124e+07 | N*s/m |
| cyy (Damping) | 9.5236e+07 | N*s/m |
| Optimization Success | True | - |
| Function Value | 1.0492e+00 | - |
| Iterations | 58 | - |
| Evaluations | 109 | - |
| Execution Time | 1.4276e+01 | s |
+---------------------------+---------------------------+---------------------------+
=====================================================================================
1.9: Final remarks#
To perform quick analyses of short bearings, ROSS provides the CylindricalBearing class, which models the fluid film using a closed-form solution to the Reynolds equation. An example of how this class is used is shown below:
from ross import CylindricalBearing
from ross.units import Q_
cylindrical = CylindricalBearing(
n=0,
speed=Q_([1500, 2000], "RPM"),
weight=525,
bearing_length=Q_(30, "mm"),
journal_diameter=Q_(100, "mm"),
radial_clearance=Q_(0.1, "mm"),
oil_viscosity=0.1,
)
To display the value of a specific bearing dynamic coefficient (e.g., kxx) for all the analyzed speeds, use:
cylindrical.kxx
array([12807959.57740873, 13005276.62769462])
To summarize the inputs (geometric parameters and operational conditions) and outputs (equilibrium position and dynamic coefficients) of the CylindricalBraring class, use:
cylindrical._hover_info()
([0],
'Cylindrical Bearing at Node: 0<br>Length: 0.0300 m<br>Journal Diameter: 0.1000 m<br>Radial Clearance: 1.0000e-04 m<br>Oil Viscosity: 0.1000 Pa·s<br>Load: 525.00 N<br>Speed: 1500.00 ... 2000.00 RPM<br>Eccentricity: 0.2663 ... 0.2126<br>Attitude Angle: 0.1989 ... 0.1617 rad<br>Sommerfeld: 3.571e+00 ... 4.762e+00<br>Kxx: 1.281e+07 ... 1.301e+07 N/m<br>Kyy: 8.815e+06 ... 8.021e+06 N/m<br>Cxx: 2.329e+05 ... 2.249e+05 N·s/m<br>Cyy: 2.949e+05 ... 2.620e+05 N·s/m<br>')
Section 2: TiltingPad#
This tutorial provides a guide on how to simulate the thermo-hydrodynamic behavior of tilting pad journal bearings and compute their dynamic coefficients using the ROSS library.
Tilting pad bearings differ from fixed-geometry bearings because each pad can rotate independently around a pivot, allowing better stability and reducing the risk of oil whirl and oil whip instabilities. The analysis typically involves solving the Reynolds equation coupled with thermal effects.
Tilting pad bearings can be analyzed using two different thermal modeling approaches available in the ROSS library:
Adiabatic model (thermal_type=”adiabatic”): Only the oil film temperature is computed. Heat transfer within the pad is neglected, resulting in a faster solution suitable for preliminary analyses.
Full thermo-hydrodynamic model (thermal_type=”full”): Includes thermal coupling between the oil film and the pad material. This approach accounts for heat conduction inside the pads and requires iterative convergence, providing more accurate results for high-load or high-speed applications.
2.1: Adiabatic example of application#
First, to run a test case using an adiabatic thermal model (where only the oil film temperature is computed), use the following code:
from ross.bearings.tilting_pad import TiltingPad
from ross.units import Q_
import plotly.io as pio
pio.renderers.default = "notebook"
tilting_pad_adiabatic = TiltingPad(
n=1, # Nodal location of the bearing in the FE model of the rotor
frequency=Q_([3000, 5000], "RPM"), # Journal rotational speed
load=[884.05, -2670.4], # External load applied to the shaft [Fx, Fy]
journal_diameter=101.6e-3, # Journal diameter
radial_clearance=74.9e-6, # Radial clearance
pad_thickness=12.7e-3, # Pad thickness
pivot_angle=Q_([18, 90, 162, 234, 306], "deg"), # Pivot angular positions
pad_arc=Q_([60, 60, 60, 60, 60], "deg"), # Angular extent of each pad
pad_axial_length=Q_([50.8e-3] * 5, "m"), # Axial length of each pad
pre_load=[0.5] * 5, # Preload factor for each pad
offset=[0.5] * 5, # Pivot offset factor for each pad
lubricant="ISOVG32", # Lubricant type
oil_supply_temperature=Q_(40, "degC"), # Oil supply temperature
nx=30, # Number of elements in circumferential direction
nz=30, # Number of elements in axial direction
thermal_type="adiabatic", # Only oil film temperature is computed
equilibrium_type="match_eccentricity", # Method to determine equilibrium
eccentricity=0.35, # Prescribed eccentricity
attitude_angle=Q_(287.5, "deg"), # Attitude angle
initial_pads_angles=[
1.0747e-03,
7.2199e-04,
2.9340e-04,
3.4885e-04,
8.1554e-04,
], # Initial guesses to find the equilibrium position of the pads using an optimization process
solver_options={"xtol": 1e-2, "ftol": 1e-2, "maxiter": 1000}, # Solver options
)
2.2: Equilibrium calculation options#
In addition to the thermal model, the TiltingPad class requires defining how the equilibrium position of the rotor is determined. This is controlled by the equilibrium_type parameter.
Two options are available:
“match_eccentricity”: Uses the prescribed eccentricity and attitude angle, optimizing only the pad rotation angles to satisfy moment equilibrium.
“determine_eccentricity”: Computes the full equilibrium condition by simultaneously determining the eccentricity, attitude angle, and pad rotation angles required to balance the applied loads.
Using prescribed eccentricity (match_eccentricity)
In this approach, the eccentricity and attitude angle are prescribed, and an optimization process is performed to determine the pad rotation angles that satisfy the moment equilibrium condition.
tilting_pad_match = TiltingPad(
n=1,
frequency=Q_([3000], "RPM"),
load=[884.05, -2670.4],
journal_diameter=101.6e-3,
radial_clearance=74.9e-6,
pad_thickness=12.7e-3,
pivot_angle=Q_([18, 90, 162, 234, 306], "deg"),
pad_arc=Q_([60] * 5, "deg"),
pad_axial_length=Q_([50.8e-3] * 5, "m"),
pre_load=[0.5] * 5,
offset=[0.5] * 5,
lubricant="ISOVG32",
oil_supply_temperature=Q_(40, "degC"),
nx=30,
nz=30,
thermal_type="adiabatic",
equilibrium_type="match_eccentricity",
eccentricity=0.35,
attitude_angle=Q_(287.5, "deg"),
)
The convergence of the optimization process can be visualized using:
tilting_pad_match.show_optimization_convergence(show_plots=True)
=====================================================================================
OPTIMIZATION CONVERGENCE - 3000.0 RPM
=====================================================================================
+---------------------------+---------------------------+---------------------------+
| Pad | Iterations | Final Residual [N] |
+---------------------------+---------------------------+---------------------------+
| 1 | 20 | 0.000803 |
| 2 | 20 | 0.034964 |
| 3 | 20 | 0.010224 |
| 4 | 20 | 0.012591 |
| 5 | 24 | 0.056801 |
+---------------------------+---------------------------+---------------------------+
=====================================================================================