# Complete modeling of hydrodynamic bearings with a boundary parameterization
approach

J.A. Mota  D.J.G. Maldonado  J.V. Valério  T.G. Ritto

###### Abstract

The present work aims to revisit the simplifications made in the Navier-Stokes
equations for the flow between two cylinders with a small thickness of
lubricating oil film. Through a dimensionless analysis, the terms of these
equations are mapped and ordered by importance for the hydrodynamic bearing
application. An effective parameterization of the geometry is proposed,
enabling a more detailed description of the problem and its adaptation to
other contexts. At the end, an elliptical partial differential equation is
reached and solved by the centered finite difference method, whose solution is
the pressure field between the cylinders. To illustrate the effectiveness of
the proposed approach, the model is applied to hydrodynamic bearings, where
the pressure field and some parameters resulting from it, such as stiffness
and damping coefficients, are computed. Based on the facilities offered by the
parameterization of the geometry, two different configurations are presented:
(1) elliptical and (2) worn bearings. Their responses are evaluated and a
comparative analysis is performed. The modeling exposed in this text, as well
as all its simulations were developed to integrate Ross-Rotordynamics, an open
library in Python, available on the GitHub platform.

###### keywords:

Lubrication Theory , Geometry parameterization , Numerical Simulation ,
Hydrodynamic Bearings , Dynamic Coefficients

††journal: ???

[inst1]organization=Institute of Computing,addressline=UFRJ, city=Rio de
Janeiro, country=Brazil

[inst2]organization=Department of Mechanical Engineering,addressline=UFRJ,
city=Rio de Janeiro, country=Brazil

##  1 Introduction

The dynamics of rotating machinery has been actively studied since the
beginning of the last century. Many industries rely on the performance of such
machines. Any failure or malfunction affects the productivity of the process.
Thus, it is of great importance to know in advance the behavior of the machine
under different operating conditions. In that study, special attention has
been drawn to elements with fluid-induced forces. This is the case of
hydrodynamic bearings and annular seals. These components are composed of a
stator, a rotor and a fluid between them that reacts to external loads applied
to the rotor. Depending on the geometry and the operating conditions, those
forces may induce high vibration on the rotor. In general, the prediction of
the rotor vibration under a certain operating condition is the main concern of
the investigations.

Journal bearings play an important role in rotating machines since they
support the rotor and the external loads applied to it. Since the reaction
forces of the fluid may increase the vibration of the machine, it is necessary
to understand the bearings forces in order to know the behavior of the machine
under different conditions. Based on this information, predictions about
malfunctions or different observed phenomena can be understand and corrected
in advance.

One of the main concerns of the rotordynamics are the high vibrations due to
instabilities and resonances. Rankine [1], Dunkerley [2], Föppl [3] and
Jeffcott [4] were the pioneers that developed investigations on the critical
speeds of a spinning rotor supported by bearings. The objective of their
studies was the stability of the rotor’s motion beyond the first critical
speed. Since then, stability has been one of the most important
characteristics of rotating machines; specially today, where there is a need
of high performance machines with reliable operating conditions. A more
detailed review of bibliography can be found in the books of Childs [5],
Muszynska [6], Tiwari [7], Vance [8] and Ishida and Yamamoto [9].

Regarding hydrodynamic bearings, the research began in the 1880s, with the
experimental investigations of Beauchamp Tower (1845-1904) and Nicolai P.
Petrov (1836-1920). However, the theoretical basis for their experimental
observation was only presented by Osborne Reynolds (1842- 1912). The three of
them can be considered the founding fathers of the hydrodynamic lubrication
and their findings are still used these days. Sommerfeld [10] found an
analytical expression for infinitely long bearings, but only Swift [11] and
Stieber [12] used realistic boundary conditions to the Reynolds equation.
Ocvirk [13] developed a detailed solution to infinitely short bearings, which
is still widely used today. For the finite bearing, a strong effort was given
by Cameron and Wood [14] with numerical calculations without the aid of
electronic devices. After the arrival of computers, Pinkus [15] was the first
researcher to use the computer for the numerical calculations of the solutions
from circular bearings, elliptical bearings and a three lobes bearing. In
1981, Singhal [16] applied finite differences to the Reynolds equation.

In addition to the lubrication theory, the studies carried out by Pina and
Carvalho [17], Andrade [18] and Queiroz [19] are relevant because their work
has inspired the parameterization used in the description of the geometry. The
well-known collection of lecture notes recorded by San Andres [20], especially
the $1$ to $5$ notes, were essential for understanding the context of
hydrodynamic bearings. The works of Frene et al [21], Hamrock [22] and Ishida
and Yamamoto [23] stand out, for this text, in what concerns the
approximations of the Reynolds equation. The three authors present basic
concepts of hydrodynamic bearing geometry and use the approximate equation for
short and infinitely long bearings. Finally, the research of Machado and
Cavalca [24] presents solutions for bearings of different geometries. Machado
[25] analyzes dynamic and operational characteristics of different
configurations of hydrodynamic radial bearings: a cylindrical, an elliptical
and a trilobular. Already Machado and Cavalca [24] expose a numerical model
that characterizes the bearing with wear on its wall, analyzing its influence
on the dynamic response of the rotor.

In this article, a study of a fluid film flow between two surfaces is carried
out in a novel way. A new and effective parameterization of the geometry is
presented to describe the journal bearing, allowing a more detailed
description of the problem. Instead of describing the geometry in terms of a
gap function, this work describes the geometry of both the rotor and stator.
Thus, different type of bearings can be modeled by simply defining the desired
shape for each surface. Once the domain is defined, the pressure distribution
of the fluid is obtained by means of the lubrication theory. First, the
governing equations of a differential fluid element are determined by the
Navier-Stokes equations and they are simplified to obtained the Reynolds
equation, whose solution is the pressure field of the fluid flow. This
equation is discretized by means of the finite differences method. The
pressure field is then obtained and verified with results available in the
literature. Afterwards, the fluid forces are linearized around the equilibrium
point of the rotor and the dynamic coefficients (stiffness and damping) are
computed by the least squares method. This methodology is tested for three
different geometries: (1) a cylindrical bearing, (2) a lemon bearing and a (c)
bearing with a worn stator. More details about the modeling presented in this
text can be found in Mota [26].

This article is organized as follows. In Section 2, the description of the
fluid flow problem is stated and detailed. In Section 3, the model of the
fluid flow is described and the governing equation relating the pressure field
and the geometrical characteristics is obtained. In Section 4, the pressure
field is obtained by using the finite differences method. Then, in Section 5
the reaction forces of the fluid film are computed by a numerical integration
of the pressure field in the whole domain. Section 6 shows the stiffness and
damping coefficients obtained with the least-square method. In Section 7, the
model is explored and three geometries are analyzed: a short bearing, an
elliptical bearing and a worn bearing. Finally, in Section 8, the conclusions
are thrown.

##  2 Problem description

In this article, the dynamics of the fluid flow in the annular space between
two cylinders is studied. This is a representation of a hydrodynamic bearing.
The system is shown in Fig. 1. The external cylinder, with radius $R_{o}$ is
called stator and the internal cylinder, with radius $R_{i}$, is the rotor
which spins with a speed $\omega$.

Figure 1: Schematic drawing of an eccentricity journal bearing

With the rotation speed, the rotor moves to an eccentric position. Any point
on the rotor’s surface can be described as (see Fig. 1)

$\displaystyle
R_{\theta}=\sqrt{R_{i}^{2}-e^{2}\sin^{2}{\alpha}}+e\cos{\alpha}\,,$ (1)

with

$\displaystyle\alpha=\begin{cases}\dfrac{3\pi}{2}-\theta+\beta\text{,}&\mbox{if
}\dfrac{\pi}{2}+\beta\leq\theta<\dfrac{3\pi}{2}+\beta\quad\textbf{({I})}\\\
\\\ -\left(\dfrac{3\pi}{2}-\theta+\beta\right)\text{,}&\mbox{if
}\dfrac{3\pi}{2}+\beta\leq\theta<\dfrac{5\pi}{2}+\beta\quad\textbf{({II})}\end{cases},$

where $e$ is the eccentricity, $\beta$ is the attitude angle, $R_{o}$ is the
stator’s radius and $\theta$ is the angle used to define $R_{\theta}$. Note
that $\theta$ is measured from an axis passing through the points of maximum
and minimum fluid width. Also, it is important to note that the difference
between the radius of the stator and the rotor is very small compared to the
radius of the stator.

##  3 Theoretical Modeling

The dynamic of the fluid between the rotor and stator is modeled using the
Navier-Stokes equations and the continuity equation. For this problem, the
equations are

$\displaystyle\rho\Bigr{(}\frac{\partial\mathbf{v}}{\partial
t}+\mathbf{v}\cdot\nabla\mathbf{v}\Bigl{)}$ $\displaystyle=\nabla\cdot\sigma,$
(2a) $\displaystyle\frac{\partial\rho}{\partial
t}+\nabla\cdot(\rho\mathbf{v})$ $\displaystyle=0,$ (2b)

where $\rho$ is the specific mass of the fluid; $\mathbf{v}$ is the velocity
field vector in cylindrical coordinates with components $u$ (axial), $v$
(radial) and $w$ (tangential); and $\sigma=-p\mathbf{I}+\tau$ is Cauchy’s
tensor in terms of the pressure field $p$, the tension tensor $\tau$ and the
identity tensor $\mathbf{I}$.

Assuming an incompressible Newtonian fluid in a steady-state laminar flow,
Eqs. (2) can be simplified [27]. First, $\rho$ can be considered constant.
Also, the shear stress is proportional to its deformation rate, i.e.,
${\mathbf{\tau}}={\mathbf{\mu}}(\nabla\mathbf{v})$. Moreover, the fluid is
considered in a steady state, where the quantities do not vary over time.
Thus, Eq. (2) can be rewritten as:

$\displaystyle\rho(\mathbf{v}\cdot\nabla\mathbf{v})$ $\displaystyle=-\nabla
p+\mu\nabla^{2}\mathbf{v},$ (3a) $\displaystyle\nabla$
$\displaystyle\cdot\mathbf{v}=0\,,$ (3b)

which are expressed in cylindrical coordinates to preserve the curvature
effects, following [17, 18]. In most references it is assumed that the radius
of curvature of the boundaries is large compared to the thickness of the
lubricating film, and these effects are neglected. Developing Eq. (3) yields
the following equations.

Axial direction,  z:
$\displaystyle{\rho\left(u\dfrac{\partial{u}}{\partial{z}}+v\dfrac{\partial{u}}{\partial{r}}+\dfrac{w}{r}\dfrac{\partial{u}}{\partial{\theta}}\right)}={-\dfrac{\partial{p}}{\partial{z}}+\mu\left(\dfrac{1}{r}\dfrac{\partial{}}{\partial{r}}\left[r\dfrac{\partial{u}}{\partial{r}}\right]+\dfrac{1}{r^{2}}\dfrac{\partial^{2}{u}}{\partial{\theta^{2}}}+\dfrac{\partial^{2}{u}}{\partial{z^{2}}}\right)}\text{.}$
(4a) Radial direction,  r:
$\displaystyle{\rho\left(v\dfrac{\partial{v}}{\partial{r}}+u\dfrac{\partial{v}}{\partial{z}}+\dfrac{w}{r}\dfrac{\partial{v}}{\partial{\theta}}-\dfrac{w^{2}}{r}\right)}$
$\displaystyle=$
$\displaystyle{{-\dfrac{\partial{p}}{\partial{r}}+\mu\left(\dfrac{\partial}{\partial
r}\left[\dfrac{1}{r}\dfrac{\partial(rv)}{\partial
r}\right]+\dfrac{1}{r^{2}}\dfrac{\partial^{2}{v}}{\partial{\theta^{2}}}-\dfrac{2}{r^{2}}\dfrac{\partial{w}}{\partial{\theta}}+\dfrac{\partial^{2}{v}}{\partial{z^{2}}}\right)}}\text{.}$
(4b) Tangent direction, $\theta$:
$\displaystyle{\rho\left(v\dfrac{\partial{w}}{\partial{r}}+u\dfrac{\partial{w}}{\partial{z}}+\dfrac{w}{r}\dfrac{\partial{w}}{\partial{\theta}}+\dfrac{vw}{r}\right)}$
$\displaystyle=$
$\displaystyle{{-\dfrac{1}{r}\dfrac{\partial{p}}{\partial{\theta}}+\mu\left(\dfrac{\partial}{\partial
r}\left[\dfrac{1}{r}\dfrac{\partial(rw)}{\partial
r}\right]+\dfrac{1}{r^{2}}\dfrac{\partial^{2}{w}}{\partial{\theta^{2}}}+\dfrac{2}{r^{2}}\dfrac{\partial{v}}{\partial{\theta}}+\dfrac{\partial^{2}{w}}{\partial{z^{2}}}\right)}}\text{.}$
(4c) Continuity:
$\displaystyle\dfrac{1}{r}\dfrac{\partial{\left(rv\right)}}{\partial{r}}+\dfrac{1}{r}\dfrac{\partial{w}}{\partial{\theta}}+\dfrac{\partial{u}}{\partial{z}}=0\text{.}$
(4d)

This system of partial differential equations is nonlinear, and the velocity
and pressure fields are unknown. After applying the proper boundary
conditions, it can be approximated numerically. However, additional
simplifications can be applied when a dimensional analysis is performed.

###  3.1 Dimensionless Navier-Stokes equations

The fact that the radial clearance $F=R_{o}-R_{i}$ is very small compared to
the radii $R_{o}$ and $R_{i}$, and the length $L$ is used to decide which
terms can be neglected. First, consider $U$ a typical speed and $L$ a
characteristic length with the same order of magnitude as $R_{o}$ and $R_{i}$.
Then, the following dimensionless variables are introduced:

$\displaystyle u=U\hat{u},\quad v=U\hat{v},\quad w=U\hat{w},\quad
p=P\hat{p},\quad P=\dfrac{\mu UL}{F^{2}},$ $\displaystyle\quad\Delta
z=L\Delta\hat{z},\quad\Delta
r=F\Delta\hat{r},\quad\Delta(r\theta)=L\Delta(\hat{r}{\theta})\,,$ (5)

where the hat symbol means dimensionless variable. Note that the expression
for $P$ was used since the pressure exhibits a higher variation near the
smallest fluid width. Now, using the dimensionless variables (3.1) in (4a)
yields:

$\displaystyle\begin{split}&\rho\left(U\hat{u}\partialderivative{U\hat{u}}{L\hat{z}}+U\hat{v}\partialderivative{U\hat{u}}{F\hat{r}}+\frac{U\hat{w}}{L\hat{r}}\partialderivative{U\hat{u}}{\theta}\right)=\\\
&-\partialderivative{P\hat{p}}{L\hat{z}}+\mu\left(\frac{1}{L\hat{r}}\partialderivative{F\hat{r}}\left[L\hat{r}\partialderivative{U\hat{u}}{F\hat{r}}\right]+\frac{1}{L^{2}\hat{r}^{2}}\partialderivative[2]{U\hat{u}}{\theta}+\partialderivative[2]{U\hat{u}}{L^{2}\hat{z}}\right)\,.\end{split}$
(6)

After some algebraic manipulations one obtains

$\displaystyle\begin{split}&Re\left(\left(\frac{F^{2}}{L^{2}}\right)\hat{u}\partialderivative{\hat{u}}{\hat{z}}+\left(\frac{F}{L}\right)\hat{v}\partialderivative{\hat{u}}{\hat{r}}+\left(\frac{F^{2}}{L^{2}}\right)\frac{\hat{w}}{\hat{r}}\partialderivative{\hat{u}}{\theta}\right)=\\\
&-\partialderivative{\hat{p}}{\hat{z}}+\left(\frac{1}{\hat{r}}\partialderivative{\hat{r}}\left[\hat{r}\partialderivative{\hat{u}}{\hat{r}}\right]+\left(\frac{F^{2}}{L^{2}}\right)\frac{1}{\hat{r}^{2}}\partialderivative[2]{\hat{u}}{\theta}+\left(\frac{F^{2}}{L^{2}}\right)\partialderivative[2]{\hat{u}}{\hat{z}}\right),\end{split}$
(7)

where the Reynolds number

$\displaystyle Re=\frac{\rho UL}{\mu}$ (8)

is used. Following the traditional theory of lubrication, the terms multiplied
by $(F/L)^{2}$ and $(F/L)$ are neglected and Eqs. (4a), (4) and (4) become:

$\displaystyle-\partialderivative{p}{z}+\frac{1}{r}\partialderivative{r}\left(r\partialderivative{u}{r}\right)$
$\displaystyle=0,$ (9a) $\displaystyle-\partialderivative{p}{r}$
$\displaystyle=0,$ (9b)
$\displaystyle\partialderivative{r}\left(\frac{1}{r}\partialderivative{(rw)}{r}\right)-\frac{1}{r}\partialderivative{p}{\theta}$
$\displaystyle=0.$ (9c)

After this process, the number of terms was reduced considerably. For the
equation along the radial direction (9b), the radial velocity disappeared and
the pressure gradient along this direction is zero.

###  3.2 Velocities

Equations (9a) and (9c) can now be integrated in order to compute the
velocities $u$ and $w$:

$u=\left(\partialderivative{p}{z}\right)\frac{1}{4}r^{2}+c_{1}\ln r+c_{2}\,,$
(10a) $w=\frac{1}{2}\partialderivative{p}{\theta}r\left(\ln
r-\frac{1}{2}\right)+c_{3}r+\frac{c_{4}}{r}\,,$ (10b)

where $c_{1}$, $c_{2}$, $c_{3}$ and $c_{4}$ are constants. Since the fluid
speed is zero at the stator and equal to the tangential speed $W$ at the
rotor’s surface, the boundary conditions are: $u(R_{o})=0$, $u(R_{\theta})=0$,
$w(R_{o})=0$ and $w(R_{\theta})=\omega R_{i}=W$. Thus, the description of the
velocities $u$ and $w$ is obtained as a function of the pressure gradient:

$\displaystyle u=$
$\displaystyle\partialderivative{p}{z}\frac{R^{2}_{\theta}}{4}\left[\left(\frac{r}{R_{\theta}}\right)^{2}-\frac{(R^{2}_{o}-R^{2}_{\theta})}{R^{2}_{i}\ln(R_{o}/R_{\theta})}\ln\left(\frac{r}{R_{\theta}}\right)-1\right],$
(11a) $\displaystyle w$ $\displaystyle=\dfrac{1}{2}\dfrac{\partial
p}{\partial\theta}\left[r\left(\ln
r-\dfrac{1}{2}\right)+Kr-\dfrac{R_{o}^{2}}{r}\left(\ln
R_{o}+K-\dfrac{1}{2}\right)\right]$
$\displaystyle+\dfrac{WR_{\theta}}{\left(R_{\theta}^{2}-R_{o}^{2}\right)}\left(r-\dfrac{R_{o}^{2}}{r}\right),$
(11b)

where

$K=\frac{1}{R^{2}_{o}-R^{2}_{\theta}}\left[R^{2}_{\theta}\left(\ln
R_{\theta}-\frac{1}{2}\right)-R^{2}_{o}\left(\ln
R_{o}-\frac{1}{2}\right)\right].$ (12)

###  3.3 Continuity Equation

In order to obtain the pressure field, the continuity equation (4d) is
integrated in the annular region of interest,

$\underset{\mathrm{I}}{\underbrace{\int^{R_{o}}_{R_{\theta}}\left(\frac{\partial{(rv)}}{\partial{r}}\right)\,dr}}+\underset{\mathrm{II}}{\underbrace{\int^{R_{o}}_{R_{\theta}}\left(\frac{\partial{w}}{\partial{\theta}}\right)\,dr}}+\underset{\mathrm{III}}{\underbrace{\int^{R_{o}}_{R_{\theta}}\left(\frac{\partial{(ru)}}{\partial{z}}\right)\,dr}}=0\,.$
(13)

Each term is detailed as follows:

$\displaystyle(\mathrm{I})\quad\int^{R_{o}}_{R_{\theta}}\left(\frac{\partial{(rv)}}{\partial{r}}\right)\,dr=R_{o}v(R_{o})-R_{\theta}v(R_{\theta})\,,$
$\displaystyle(\mathrm{II})\quad\int^{R_{o}}_{R_{\theta}}\left(\frac{\partial{w}}{\partial{\theta}}\right)\,dr=\dfrac{\partial}{\partial\theta}\int^{R_{o}}_{R_{\theta}}w\,dr-\left[w(R_{o})\dfrac{\partial
R_{o}}{\partial\theta}-w(R_{\theta})\dfrac{\partial
R_{\theta}}{\partial\theta}\right]\,,$
$\displaystyle(\mathrm{III})\quad\int^{R_{o}}_{R_{\theta}}\left(\frac{\partial{(ru)}}{\partial{z}}\right)\,dr=\dfrac{\partial}{\partial
z}\int^{R_{o}}_{R_{\theta}}(ru)\,dr-\left[u(R_{o})\dfrac{\partial
R_{o}}{\partial z}-u(R_{\theta})\dfrac{\partial R_{i}}{\partial z}\right]\,.$

Some remarks about the boundary conditions should be registered. The radial
speed at the stator $v(R_{o})$ is zero; although the speed at the rotor’s
surface is tangential to it, it has a projection along the radial direction in
the coordinate system with origin at the stator’s center (see Fig. 2); since
the external radius does not vary along the tangential direction,
$\dfrac{\partial R_{o}}{\partial\theta}$ is zero; due to the eccentricity,
$\dfrac{\partial R_{\theta}}{\partial\theta}$ is not zero; the speeds
$u(R_{o})$ and $u(R_{\theta})$ are zero.

After applying the aforementioned boundary conditions, Eq. (13) is rewritten
as:

$\displaystyle\dfrac{\partial}{\partial\theta}\int^{R_{o}}_{R_{\theta}}w\,dr+w(R_{\theta})\dfrac{\partial
R_{\theta}}{\partial\theta}+\dfrac{\partial}{\partial
z}\int^{R_{o}}_{R_{\theta}}ru\,dr-R_{\theta}v(R_{\theta})=0$ (14)

The remaining boundary conditions $w(R_{\theta})$ and $v(R_{\theta})$ are
obtained by kinematics. First, consider a point $A$ at the surface of the
rotor, as shown in Fig. 2. The poition vector of this point, with respect to a
cylindrical coordinate system with origin at the stator’s center, can be
described as

$\displaystyle\mathbf{p}_{A}=R_{\theta}\mathbf{e}_{\theta},$ (15)

and its time derivative as

$\displaystyle\mathbf{v}_{A}=\dot{R_{\theta}}\,\mathbf{e}_{r}+R_{\theta}\,\dot{\mathbf{e}}_{r}=\underbrace{\omega\dfrac{\partial
R_{\theta}}{\partial\theta}}_{v_{\text{rad}}}\,\mathbf{e}_{r}+\underbrace{\omega
R_{\theta}}_{v_{\text{tan}}}\,\mathbf{e}_{\theta},$ (16)

where

$\displaystyle v_{\text{rad}}$
$\displaystyle=u(R_{\theta})=\omega\dfrac{\partial
R_{\theta}}{\partial\theta},$ (17a) $\displaystyle v_{\text{tan}}$
$\displaystyle=w(R_{\theta})=\omega R_{\theta}.$ (17b)

Figure 2: Graphical description of the speed of point A.

The boundary conditions in Eq. (17) are substituted in Eq. (14), which yields:

$\displaystyle\dfrac{\partial}{\partial\theta}\int^{R_{o}}_{R_{\theta}}w\,dr+\dfrac{\partial}{\partial
z}\int^{R_{o}}_{R_{\theta}}ru\,dr=0$ (18)

Finally, the speeds $u$ and $w$ in Eqs. (11a) and (11b) are substituted in Eq.
(18) and integrated:

$\partialderivative{\theta}\left(C_{1}\partialderivative{p}{\theta}\right)+\partialderivative{z}\left(C_{2}\partialderivative{p}{z}\right)=\partialderivative{\theta}C_{0}$
(19)

where:

$C_{0}=-WR_{\theta}\left[\ln\left(\frac{R_{o}}{R_{\theta}}\right)\left(1+\frac{R^{2}_{\theta}}{(R^{2}_{o}-R^{2}_{\theta})}\right)-\frac{1}{2}\right]$
(20)
$\displaystyle\begin{split}C_{1}=&\frac{R_{\theta}}{2}\left[\frac{1}{2R_{\theta}}[R^{2}_{o}\ln
R_{o}-R^{2}_{\theta}\ln
R_{\theta}-\left(R^{2}_{o}-R^{2}_{\theta}\right)\left(1+K\right)]\right]-\\\
&\frac{R^{2}_{\theta}}{2\mu}\left[\left(\ln
R_{\theta}-\frac{1}{2}+K\right)\ln\left(\frac{R_{o}}{R_{\theta}}\right)\right]\end{split}$
(21)
$\displaystyle\begin{split}C_{2}=&\frac{-R^{2}_{\theta}}{8}\left[R^{2}_{o}-R^{2}_{\theta}-\frac{(R^{4}_{o}-R^{4}_{\theta})}{2R^{2}_{\theta}}\right]+\\\
&\left(\frac{R^{2}_{o}-R^{2}_{\theta}}{R^{2}_{r}\ln(R_{o}/R_{\theta})}\right)\left[R^{2}_{o}\ln\left(\frac{R_{o}}{R_{\theta}}\right)-\frac{(R^{2}_{o}-R^{2}_{\theta})}{2}\right]\end{split}$
(22)

The original Navier-Stokes equation, Eq. (4), was reduced to an elliptical
differential equation, Eq. (19), which is known as the Reynolds equation that
can be approximated numerically.

##  4 Numerical Approximation for the pressure field

In this section, an approximate solution to Eq. (19) is obtained for the
pressure field $p$. A central finite difference approximation method is used
(see Ref. [28] for more details).

###  4.1 Discretization of the domain

The cylindrical region of interest is opened up and the discretization is
performed in a rectangular region with $N_{z}\times N_{\theta}$ nodes. Thus,
the fluid film has $N_{z}$ divisions along the axial direction and $N_{z}$
division along the tangential direction, both equally spaced. Fig. 3 shows
this procedure.

Figure 3: Rectangular mesh used for discretization of the equation and graphic
translation of the centered finite difference method; adapted from [18].

###  4.2 Central finite differences formula

The partial derivatives of Eq. (19) are dicretized according to the following
formulas:

$\displaystyle\partialderivative{f(i,j)}{\theta}=\frac{f(i,j-1)-f(i,j+1)}{2\Delta\theta}\quad,\quad\partialderivative{f(i,j)}{z}=\frac{f(i-1,j)-f(i+1,j)}{2\Delta
z}$ (23)

After the substitution of Eq. (23) in Eq. (19), the following algebraic
equation is obtained:

$\displaystyle p_{i-1,j}\frac{(C_{2(i-1,j)})}{\Delta
z^{2}}+p_{i,j-1}\frac{(C_{1(i,j-1)})}{\Delta\theta^{2}}-p_{i,j}\left(\frac{(C_{1(i,j)}+C_{1(i,j-1)})}{\Delta\theta^{2}}+\frac{(C_{1(i,j)}+C_{1(i,j-1)})}{\Delta
z^{2}}\right)+$ $\displaystyle
p_{i,j+1}\frac{(C_{1(i,j)})}{\Delta\theta^{2}}+p_{i+1,j}\frac{(C_{2(i,j)})}{\Delta
z^{2}}=\frac{1}{\Delta\theta}\left[C_{0W(i,j)}-C_{0W(i,j-1)}\right],$ (24)

where the pressure at the node $(i,j)$ is written as $p(i,j)$. The boundary
conditions of Eq. (4.2) are:

$\displaystyle p(z=0)=p(1,j)=P_{in},$ $\displaystyle
p(z=L)=p(N_{Z},j)=P_{out},$ $\displaystyle
p(\theta=0)=p(\theta=2\pi)=p(i,1)=p(i,N_{\theta}).$

Eq. (4.2) can be expressed in a matrix form as

$\displaystyle\mathbf{M}\mathbf{p}=\mathbf{f}\,,$ (25)

where

$\displaystyle\mathbf{p}=\begin{bmatrix}p_{0,0}&p_{0,1}&\dots&p_{N_{Z},N_{\theta}}\end{bmatrix}^{T},$
(26)
$\displaystyle\mathbf{f}=\begin{bmatrix}f_{0,0}&f_{0,1}&\dots&f_{N_{Z},N_{\theta}}\end{bmatrix}^{T},$
(27)

and $\mathbf{M}$ is a sparse matrix. Eq. (25) is assembled in Python and
solved with the
SymPy111https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html
library.

##  5 Fluid Film Forces

Integrating the pressure field one obtains the fluid film force that is
applied to the rotor. The decomposition of this force in the orthogonal axis
N-T shown in Fig. 1 are given by [23]:

$\displaystyle f_{N}$
$\displaystyle=-R_{i}\int_{0}^{L}\int_{0}^{\pi}p(\theta,z)\cos{\theta}d\theta
dz$ (28a) $\displaystyle f_{T}$
$\displaystyle=R_{i}\int_{0}^{L}\int_{0}^{\pi}p(\theta,z)\sin{\theta}d\theta
dz$ (28b)

or in a X-Y axis by:

$\displaystyle f_{x}=f_{T}\cos{(\beta)}-f_{N}\sin{(\beta)}$ (29a)
$\displaystyle f_{y}=f_{T}\sin{(\beta)}+f_{N}\cos{(\beta)}$ (29b)

Since the pressure field is discrete, the double integrals are computed
numerically using the Simpson’s method.

###  5.1 Equilibrium position

During the normal operation of a rotating machine supported by hydrodynamic
bearings, external loads are applied to the rotor and the fluid film forces
oppose to them. If only the weight of the rotor is considered, the center of
the spinning rotor stays in an equilibrium position where the weight is
balanced with the fluid film forces. This position can be described by the
eccentric displacement $e$ and the attitude angle $\beta$, as shown in Fig. 1.
At this position, the following condition must be satisfied:

$\displaystyle f_{x_{0}}=$ $\displaystyle 0,$ (30a) $\displaystyle f_{y_{0}}=$
$\displaystyle W,$ (30b)

where $W$ is the weight of the rotor.

An iterative method is used to find the position of the rotor that meets the
condition of Eq. (30). Starting with an initial (guess) position, the
difference between the current forces and the forces shown in Eq. (30) is
calculated. Next, the rotor center is slightly moved, and the forces are
recalculated. When a tolerance condition is met, the method stops. More
details can be found in the Python library optimize.least_squares

The static performance characteristics of the bearing are listed in a
nondimensional number known as Sommerfeld Number ($S$) [29]:

$\displaystyle S=\frac{\mu\omega LR_{o}^{3}}{\pi WF^{2}}\text{,}$ (31)

where $F=R_{o}-R_{i}$ is the radial clearance.

This nondimensional number is directly linked to the equilibrium position.
Fig. 4 illustrates this relation.

Figure 4: Path taken by the center of the rotor related to the Sommerfeld
Number; adapted from [30].

##  6 Dynamic Coefficients

In the previous section, the fluid film forces calculated are valid only at
the static equilibrium of the rotor’s center. If the rotor performs small
movements around such position, the fluid film forces vary. One way to
characterize the relation between the forces and small displacements around
the equilibrium position is to determine the stiffness and damping
coefficients.

###  6.1 Radial and tangential speeds

Here, it is considered that the rotor is not only spinning but also moving
along the horizontal / vertical directions. Thus, the speeds $w(R_{\theta})$
and $v(R_{\theta})$ in the continuity Eq. (14) must be determined with this
new dynamic state.

Let us consider again a point $A$ belonging to the surface of the rotor. If
the rotor is spinning and the velocity of its center is being perturbed, its
total velocity is

$\displaystyle\mathbf{v}_{A}=\mathbf{v}_{r}+\mathbf{v}_{p}$ (32)

where $\mathbf{v}_{r}$ is the velocity of $A$ due to the rotation speed and
$\mathbf{v}_{p}$ is the velocity acquired by the perturbation. In order to
derive an equation for $\mathbf{v}_{p}$, the small displacements around the
equilibrium position $(x_{0},y_{0})$ are assumed of the following form:

$\displaystyle x=x_{0}+x_{p}\sin(\omega_{p}t),$ (33a) $\displaystyle
y=y_{0}+y_{p}\sin(\omega_{p}t),$ (33b)

where $x_{p}$ and $y_{p}$ are the maximum amplitudes of the perturbations and
$\omega_{p}$ is the perturbation frequency along each axis. The time
derivative of Eq. (33) gives the two terms of $\mathbf{v}_{p}$:

$\displaystyle(\mathbf{v}_{p})_{x}=\omega_{p}x_{p}\cos(\omega_{p}t)\mathbf{e}_{x},$
(34a)
$\displaystyle(\mathbf{v}_{p})_{y}=\omega_{p}y_{p}\cos(\omega_{p}t)\mathbf{e}_{y},$
(34b)

where $(\mathbf{v}_{p})_{x}$ and $(\mathbf{v}_{p})_{y}$ stand for the speeds
when perturbations along the horizontal and vertical directions are applied,
respectively. In cylindrical coordinates, and using Eq. (34), one arrives to:

$\displaystyle(\mathbf{v}_{A})_{x}=\underset{(v_{\text{rad}})_{x}}{\underbrace{\left(\omega\dfrac{\partial
R_{\theta}}{\partial\theta}+\omega_{p}x_{p}\cos{(\omega_{p}t)\cos{(\theta)}}\right)}}\mathbf{e}_{r}+\underset{(v_{\text{tan}})_{x}}{\underbrace{\left(\omega
R_{\theta}-\omega_{p}x_{p}\cos{(\omega_{p}t)}\sin{(\theta)}\right)}}\mathbf{e}_{\theta}$
$\displaystyle(\mathbf{v}_{A})_{y}=\underset{(v_{\text{rad}})_{y}}{\underbrace{\left(\omega\dfrac{\partial
R_{\theta}}{\partial\theta}+\omega_{p}y_{p}\cos{(\omega_{p}t)\sin{(\theta)}}\right)}}\mathbf{e}_{r}+\underset{(v_{\text{tan}})_{y}}{\underbrace{\left(\omega
R_{\theta}-\omega_{p}y_{p}\cos{(\omega_{p}t)}\cos{(\theta)}\right)}}\mathbf{e}_{\theta}$
(35)

where

$\displaystyle(v_{\text{rad}})_{x}$
$\displaystyle=u(R_{\theta})_{x}=\omega\dfrac{\partial
R_{\theta}}{\partial\theta}+\omega_{p}x_{p}\cos{(\omega_{p}t)\cos{(\theta)}},$
(36a) $\displaystyle(v_{\text{tan}})_{x}$
$\displaystyle=w(R_{\theta})_{x}=\omega
R_{\theta}-\omega_{p}x_{p}\cos{(\omega_{p}t)}\sin{(\theta)},$ (36b)
$\displaystyle(v_{\text{rad}})_{y}$
$\displaystyle=u(R_{\theta})_{y}=\omega\dfrac{\partial
R_{\theta}}{\partial\theta}+\omega_{p}y_{p}\cos{(\omega_{p}t)\sin{(\theta)}},$
(36c) $\displaystyle(v_{\text{tan}})_{y}$
$\displaystyle=w(R_{\theta})_{y}=\omega
R_{\theta}-\omega_{p}y_{p}\cos{(\omega_{p}t)}\cos{(\theta)}.$ (36d)

###  6.2 Additional terms in the continuity equation

Equations (36)a-d are used as the new boundary conditions in the continuity
equation, Eq. (14). This consideration affects only the coefficient $C_{0}$ of
Eq. (14), which is

$\displaystyle(C_{0})_{x}$
$\displaystyle=C_{0}+R_{\theta}\omega_{p}x_{p}\cos{(\omega_{p}t)}\sin{(\theta)},$
(37a) $\displaystyle(C_{0})_{y}$
$\displaystyle=C_{0}+R_{\theta}\omega_{p}y_{p}\cos{(\omega_{p}t)}\cos{(\theta)}.$
(37b)

Thus, for an excitation along the $x$ direction, the fluid film forces are
computed according to the equation:

$\displaystyle(f_{x}(t))_{x}=$ $\displaystyle
R_{i}\\!\\!\\!\int_{0}^{L}\\!\\!\\!\\!\int_{0}^{\pi}\\!\\!\\!\\!p(\theta,z,t;x_{p},\omega_{p})\sin{\theta}\cos{\beta}d\theta
dz$ $\displaystyle\hskip
85.35826pt+R_{i}\\!\\!\\!\int_{0}^{L}\\!\\!\\!\\!\int_{0}^{\pi}\\!\\!\\!\\!p(\theta,z,t;x_{p},\omega_{p})\cos{\theta}\sin{\beta}d\theta
dz$ (38a) $\displaystyle(f_{y}(t))_{x}=$ $\displaystyle
R_{i}\\!\\!\\!\int_{0}^{L}\\!\\!\\!\\!\int_{0}^{\pi}\\!\\!\\!\\!p(\theta,z,t;x_{p},\omega_{p})\sin{\theta}\sin{\beta}d\theta
dz$ $\displaystyle\hskip 85.35826pt-
R_{i}\\!\\!\\!\int_{0}^{L}\\!\\!\\!\\!\int_{0}^{\pi}\\!\\!\\!\\!p(\theta,z,t;x_{p},\omega_{p})\cos{\theta}\cos{\beta}d\theta
dz,$ (38b)

where $p(\theta,z,t;x_{p},\omega_{p})$ is the pressure field computed at time
$t$ with perturbation amplitude and frequency $x_{p}$ and $\omega_{p}$,
respectively. Similar equations are valid for a perturbation along the $y$
direction.

###  6.3 Stiffness and damping coefficients

The forces of the fluid film can be considered general functions of the
displacements and speeds of the center of the rotor. Thus, if small
displacements around the center are applied, the expansions of the fluid
forces in Taylor series yield

$\displaystyle f_{x}$ $\displaystyle=f_{x_{0}}+K_{xx}\Delta x+K_{xy}\Delta
y+C_{xx}\Delta\dot{x}+C_{xy}\Delta\dot{y},$ (39a) $\displaystyle f_{y}$
$\displaystyle=f_{y_{0}}+K_{yx}\Delta x+K_{yy}\Delta
y+C_{yx}\Delta\dot{x}+C_{yy}\Delta\dot{y}.$ (39b)

where

$\displaystyle
K_{xx}=\partialderivative{f_{x}}{x},K_{xy}=\partialderivative{f_{x}}{y},$
$\displaystyle
K_{yx}=\partialderivative{f_{y}}{x},K_{yy}=\partialderivative{f_{y}}{y},$
$\displaystyle
C_{xx}=\partialderivative{f_{x}}{\dot{x}},C_{xy}=\partialderivative{f_{x}}{\dot{y}},$
$\displaystyle
C_{yx}=\partialderivative{f_{y}}{\dot{x}},C_{yy}=\partialderivative{f_{y}}{\dot{y}}.$

The terms $K_{xx}$, $K_{yy}$, $C_{xx}$ and $C_{yy}$ are called the direct
stifness coefficients and direct damping coefficients, respectively. $K_{xy}$,
$K_{yx}$, $C_{xy}$ and $C_{yx}$ are called the cross-coupled stiffness
coefficients and cross-coupled damping coefficients, respectively.

By including the effects of the perturbation in the simplified equation (19),
it is possible to rewrite the Taylor Series expansion of the forces in Eq.
(39) as a function of the time $t$. Thus, for each direction of disturbance
the fluid film forces are:

$\displaystyle\left(f_{x}(t)\right)_{x}$ $\displaystyle=f_{x_{0}}+k_{xx}\Delta
x(t)+c_{xx}\Delta\dot{x}(t)\text{,}$ $\displaystyle\left(f_{y}(t)\right)_{x}$
$\displaystyle=f_{y_{0}}+k_{yx}\Delta x(t)+c_{yx}\Delta\dot{x}(t)\text{,}$
$\displaystyle\left(f_{x}(t)\right)_{y}$ $\displaystyle=f_{x_{0}}+k_{xy}\Delta
y(t)+c_{xy}\Delta\dot{y}(t)\text{,}$ $\displaystyle\left(f_{y}(t)\right)_{y}$
$\displaystyle=f_{y_{0}}+k_{yy}\Delta y(t)+c_{yy}\Delta\dot{y}(t).$ (40)

Note that the $\Delta x(t)$, $\Delta\dot{x}(t)$, $\Delta y(t)$ and
$\Delta\dot{y}(t)$ are known by Eq. (33); $f_{x0}$ and $f_{y0}$ are the forces
at the equilibrium position from Eq. (30); and $(f_{x}(t))_{x}$,
$(f_{y}(t))_{x}$, $(f_{x}(t))_{y}$ and $(f_{y}(t))_{y}$ are obtained from Eq.
(38).

Equation (6.3) can be written in matrix form as:

  * 1.

Towards $x$:

$\displaystyle\underset{\mathbf{A}_{x}}{\underbrace{\begin{bmatrix}\Delta
x[0]&\Delta\dot{x}[0]\\\ \Delta x[1]&\Delta\dot{x}[1]\\\ \vdots&\vdots\\\
\Delta x[N-1]&\Delta\dot{x}[N-1]\end{bmatrix}}}\begin{bmatrix}k_{xx}&k_{yx}\\\
c_{xx}&c_{yx}\end{bmatrix}=\underset{\mathbf{F}_{x}}{\underbrace{\begin{bmatrix}(f_{x}[0])_{x}-f_{x_{0}}&(f_{y}[0])_{x}-f_{y_{0}}\\\
(f_{x}[1])_{x}-f_{x_{0}}&(f_{y}[1])_{x}-f_{y_{0}}\\\ \vdots&\vdots\\\
(f_{x}[N-1])_{x}-f_{x_{0}}&(f_{y}[N-1])_{x}-f_{y_{0}}\end{bmatrix}}}\text{.}$
(41)

  * 1.

Towards $y$:

$\displaystyle\underset{\mathbf{A}_{y}}{\underbrace{\begin{bmatrix}\Delta
y[0]&\Delta\dot{y}[0]\\\ \Delta y[1]&\Delta\dot{y}[1]\\\ \vdots&\vdots\\\
\Delta y[N-1]&\Delta\dot{y}[N-1]\end{bmatrix}}}\begin{bmatrix}k_{xy}&k_{yy}\\\
c_{xy}&c_{yy}\end{bmatrix}=\underset{\mathbf{F}_{y}}{\underbrace{\begin{bmatrix}(f_{x}[0])_{y}-f_{x_{0}}&(f_{y}[0])_{y}-f_{y_{0}}\\\
(f_{x}[1])_{y}-f_{x_{0}}&(f_{y}[1])_{y}-f_{y_{0}}\\\ \vdots&\vdots\\\
(f_{x}[N-1])_{y}-f_{x_{0}}&(f_{y}[N-1])_{y}-f_{y_{0}}\end{bmatrix}}}\text{.}$
(42)

Or, in order to simplify the notation:

$\displaystyle\mathbf{A}_{x}\begin{bmatrix}k_{xx}&k_{yx}\\\
c_{xx}&c_{yx}\end{bmatrix}=\mathbf{F}_{x}\quad\quad\text{e}\quad\quad\mathbf{A}_{y}\begin{bmatrix}k_{xy}&k_{yy}\\\
c_{xy}&c_{yy}\end{bmatrix}=\mathbf{F}_{y}\text{.}$ (43)

In order to calculate the values of the coefficients, the Least Squares method
is employed. For that, one uses the Moore-Penrose pseudo inverse of the
matrices $\mathbf{A}_{x}$ and $\mathbf{A}_{y}$.

$\displaystyle\mathbf{A}_{x}^{+}=\left(\mathbf{A}_{x}^{T}\mathbf{A}_{x}\right)^{-1}\mathbf{A}_{x}^{T}\quad\quad\text{e}\quad\quad\mathbf{A}_{y}^{+}=\left(\mathbf{A}_{y}^{T}\mathbf{A}_{y}\right)^{-1}\mathbf{A}_{y}^{T}\text{.}$
(44)

And finally, the stiffness and damping coefficients are obtained:

$\displaystyle\begin{bmatrix}k_{xx}&k_{yx}\\\
c_{xx}&c_{yx}\end{bmatrix}=\mathbf{A}_{x}^{+}\mathbf{F}_{x}\quad\quad\text{e}\quad\quad\begin{bmatrix}k_{xy}&k_{yy}\\\
c_{xy}&c_{yy}\end{bmatrix}=\mathbf{A}_{y}^{+}\mathbf{F}_{y}\text{.}$ (45)

##  7 Results and Discussion

The choice of parameterization of the model presented in this work allows easy
adaptation to different geometries and operational conditions of rotating
machine components. This section analyzes three different cases: (1) short
journal bearing, (2) elliptical bearing, and (3) worn bearing.

The algorithm used was developed in Python, and is available on the GitHub
platform, through the software Ross-Rotordynamics222https://github.com/ross-
rotordynamics/ [31], an open source library for rotodynamic analysis.

###  7.1 Short journal bearing

In the literature, it is common to use the Reynolds equation approximated for
bearings with short or infinitely long length. In this way, there an
analytical resolution is possible.

In order to verify the computational model developed in the present work,
comparisons will be made with the results obtained by the current literature
approximations. According to [21], a bearing can be considered short if the
ratio $L/D$ ($D=2R_{o}$ is the diameter) is less than or equal to $1/8$ and
long if it is greater than $4$. However, according to [20], the model for the
short bearing provides accurate results (in the case of simple cylindrical
bearings with small to moderate values of the bearing eccentricity, $e\leq
0.75F$) already with the ratio $L/D\leq 1/2$.

To analyze the short journal bearing, the pressure variation in the $z$
direction is considered much larger than in the $\theta$ direction, i.e.
$\partial p/\partial\theta\ll\partial p/\partial z$. Thus, the first term in
the Eq. (11a) is neglected, and one obtains [23]:

$p=\dfrac{-3\mu\epsilon\omega\sin{\theta}}{\left(R_{\theta}-R_{i}\right)^{2}\left(1+\epsilon\cos{\theta}\right)^{3}}\left[\left(z-\dfrac{L}{2}\right)^{2}-\dfrac{L^{2}}{4}\right]\,,$
(46)

where $\epsilon=\dfrac{e}{R_{o}-R_{i}}$ is the eccentricity ratio and $\omega$
the rotor speed.

Element | Value
---|---
Rotor radius $(R_{i})$ | $0.2\left[m\right]$
Gap $(F)$ | $0.0001\left[m\right]$
Load $(W)$ | $50\left[N\right]$
Ratio $L/D$ | $1/8$
Viscosity $(\mu)$ | $0.015\left[Pa.s\right]$
Rotation speed $(\omega)$ | $10.472\left[rad/s\right]$
Table 1: Parameters used for the simulation on a short cylindrical bearing.

The numerical simulation, using data shown in Table 1, generates a relative
error $E\approx 10^{-3}$. Figure 5 shows the pressure field computed together
with the reference result [23]. This result indicates that the present model
generate consistent results.

Figure 5: Pressure field obtained for a short journal bearing; Table 1.
Comparison of the current results with the one found in [23].

The forces can be obtained integrating the pressure field given in Eq. (46)
using Eq. (28). The relative errors, between the forces computed numerically
and analytically (pressure fields shown in Fig. 5), are $E_{N}\approx 0.02$
for the radial force and $E_{T}\approx 0.005$ for the tangential force.
Concerning the equilibrium position, the are were also small:
$f_{x}\approx|f_{y}-W|\approx 10^{-6}[N]$.

The errors in the coefficients are obtained comparing the numerical results
with the ones presented in [30]. The following relative errors are obtained:

$\displaystyle\begin{bmatrix}E_{k_{xx}}&E_{k_{xy}}\\\
E_{k_{yx}}&E_{k_{yy}}\end{bmatrix}\approx\begin{bmatrix}0.017&0.003\\\
0.001&0.04\end{bmatrix}\quad\text{e}\quad\begin{bmatrix}E_{c_{xx}}&E_{c_{xy}}\\\
E_{c_{yx}}&E_{c_{yy}}\end{bmatrix}\approx\begin{bmatrix}0.005&0.04\\\
0.06&0.002\end{bmatrix}$ (47)

Now that the proposed computational model was verified with the short bearing
approximation found in the literature, the new geometry parameterization,
developed in Section 2, can be explored in other bearing configurations, or
other machine elements.

###  7.2 Elliptical bearing

Cylindrical geometry bearings, at high rotation speeds and with small external
loads applied, are prone to instability or high vibrations [25]. This is
because, in these situations, the machine tends to operate centered, reducing
the variation in the thickness of the lubricating film and, consequently, the
ability to generate hydrodynamic pressure from the bearing. Several
configurations of bearings were developed in order to circumvent this problem,
producing preloads that make the shaft work outside the new established
centers, thus creating the desired wedge effect.

A feature commonly found in high-speed machines is the elliptical bearing or
”lemon bearing”, as it is also known [32]. This is a variation of the
cylindrical bearing with axial groove and reduced clearance in one direction.
Ostayen and Beek [32] claim that this bearing configuration is relatively
cheap and easy to manufacture, as it is made from adaptations in the
cylindrical bearing itself, being composed of two circular arcs (cuts from the
original structure), with aligned centers, generating two horizontal slits, as
shown in Fig. 6.

Figure 6: Schematic drawing of an elliptical bearing.

To consider this new geometry, stator radius is is no longer constant in
$\theta$. As seen in Fig. 6, the new stator is composed of the arc $C_{1}$,
with center in $O_{1}$, joined to the arc $C_{2}$, centered in $O_{2}$, both
with radius $R_{o}$. In this new configuration, the centers are at a distance
$\epsilon$ from the origin, called ellipticity. Now, the stator radius is
described from the origin. This new distance will be called $R_{o}^{*}$ and it
varies along the angular position.

Using the cosine law, the following relation is obtained:

$\displaystyle
R_{o}^{*}=\sqrt{R_{o}^{2}-\epsilon^{2}\sin^{2}{\alpha}}+\epsilon\cos{\alpha},$
(48)

where $\alpha=$
$\begin{cases}\pi/2+\theta\text{,}&\mbox{se}\quad\theta\in\text{1st
quadrant}\\\ 3\pi/2+\theta\text{,}&\mbox{se}\quad\theta\in\text{2nd
quadrant}\\\ \theta-\pi/2\text{,}&\mbox{se}\quad\theta\in\text{3rd
quadrant}\\\ 5\pi/2-\theta\text{,}&\mbox{se}\quad\theta\in\text{4th
quadrant}\end{cases}.$

Another important parameter to be defined is the preload $m$ which, in this
text, will be established as:

$\displaystyle m=\dfrac{\epsilon}{F}$ (49)

where $\epsilon$ is the ellipticity and $F=R_{o}-R_{i}$ is the radial
clearance.

For $m=0$, the bearing becomes cylindrical, while when $m$ tends to $1$ the
stator arcs tend to touch the axis. Garner et al [33] suggest that the maximum
preload adopted should be 0.6. According to the authors, values above 0.7 can
be used in an effort to improve stability, however, greater lateral clearances
are generated, causing the horizontal dynamic coefficients to decrease, in
addition to requiring a greater flow of oil.

Element | Value
---|---
Stator radius $(R_{o})$ | $0.015\left[m\right]$
Gap $(F)$ | $9\cdot 10^{-5}\left[m\right]$
Length $(L)$ | $0.02\left[m\right]$
Load $(W)$ | $100\left[N\right]$
Viscosity $(\mu)$ | $5.449\cdot 10^{-2}\left[Pa.s\right]$
Rotation speed $(\omega)$ | $261.8\left[rad/s\right]$
Table 2: Numerical simulation for the elliptical bearing. Fonte: Adapted from
[34]

Given the data shown in Table 2, the results for the pressure on the
elliptical bearing, using the modeling presented in this work, are compatible
with those exposed in [34]. Figures 7 (a-b) illustrate the pressure field for
a bearing described as elliptical, but with $m=0$, that is, cylindrical. In
contrast, the graphs shown in Figs. 8 (a-b) and 9 (a-b), indicate the pressure
behavior for preloads equal to $0.4$ and $0.6$, respectively.

(a) Pressure distribution.

(b) Side view.

Figure 7: Cylindrical bearing $(m=0)$.

(a) Pressure distribution.

(b) Side view.

Figure 8: Elliptical bearing with $(m=0.4)$.

(a) Pressure distribution.

(b) Side view.

Figure 9: Elliptical bearing with $(m=0.6)$.

It can be seen that Figs. 8 and 9 show two pressure peaks: one for the upper
bearing arch ($0$ to $\pi$) and one for the lower arc ($\pi$ to $2\pi$). This
phenomenon can be explained due to the configuration of the elliptical
bearing. Even if the shaft operates centered, there will be variations in the
thickness of the fluid, causing the oil film to be compressed, both in the
upper and lower arc.

Analyzing also the graph of Fig. 10, it is understood that the reaction force
increases when the preload increases, and, hence the pressure. For $m=0.8$,
the pressures generated are considerably higher than for other preloads, being
an option for cases in which a greater demand is sought to rotor stability.

Figure 10: Pressure along $\theta$ in the middle of the elliptical bearing for
different values of the preload $m$.

With regard to the equilibrium position, Fig. 11(a) illustrates a different
behavior comparing to the cylindrical bearing. That is, the Sommerfeld Number
used in the cylindrical bearing is not valid to analyze the static performance
characteristics of the elliptical bearing. For the elliptical bearing,
increasing the preload, first the rotor moves vertically and, after $m=0.3$,
it walks in the horizontal direction, approaching the center of the stator,
always reducing the eccentricity (Fig. 11(b)).

(a) Rotor center.

(b) Variation of eccentricity.

Figure 11: Variation of the equilibrium position on the elliptical bearing for
different preloads. Figure 12: Behavior of radial and tangential forces on the
elliptical bearing for different preloads.

The forces of the fluid film can be seen in Fig. 12, in which it is evident
that the tangential force is always much higher than the radial force for any
preload applied to the elliptical bearing.

It is also noticeable, when analyzing Figs. 13 (a-b), that the dynamic
coefficients of the elliptical bearing suffer great variation between their
terms. For the preload $m$ up to $0.2$, the terms $k_{xx}$, $k_{xy}$ and
$k_{yy}$ are very close to each other. From $m=0.3$ it is observed that the
difference between these terms starts to increase. The cross term $k_{xy}$ is
increasingly larger than the direct term $k_{xx}$, which might indicate
instability, as pointed out by [34]. Concerning the damping coefficients, all
terms increase as the preload increases. It is worth remembering that the main
function of damping is to mitigate the effects of vibration on the system, and
that high values of these coefficients are desired.

(a) Stiffness coefficients.

(b) Damping coefficients.

Figure 13: Dynamic coefficients of the elliptical bearing for different
preloads.

The difference between the results of the cylindrical and elliptical bearing
is clear. In the first configuration, all the factors that reduce the
eccentricity imply a decrease in the pressure field, forces and variation of
the dynamic coefficients. For the elliptical bearing configuration, the exact
opposite occurs. By increasing the preload, rotor and stator become closer to
concentricity and, this yields an increase in all the results just mentioned.
However, the increase of the stiffness cross coefficient indicates that the
system risks to be instable.

###  7.3 Wear bearing

Although lubrication reduces the friction between the metallic surfaces of the
bearing, these structures usually suffer wear after a long operating period or
due to a certain number of starting cycles. One of the most common causes of
interruptions in rotating systems is the occurrence of failures related to
hydrodynamic bearings [24]. As the wear suffered by these components affects
the radial clearance, they certainly influence the results of the pressure
field and, consequently, the reaction forces and dynamic coefficients. For
this reason, it is relevant to investigate the effect of wear, under certain
operating conditions, on the dynamic response of the system.

The wear geometry that will be used in this text was adapted from [24]. The
description was initially proposed by [35] and validated experimentally [36].

Figure 14: Schematic drawing of a worn bearing.

Dufrane et al [35] consider the wear thickness uniform in the axial direction,
and located at the bottom of the bearing, symmetrically around the $Y$ axis.
On the other hand, Machado and Cavalca [24] states that, for practical cases,
most wear and tear occurs in the fourth quadrant, the location of the axis
balance position for an anti-clockwise rotation. For this reason, taking into
account that the spent region is symmetrical in relation to the greatest depth
$d_{0}$, an angular displacement $\gamma$ will be considered from the vertical
axis to the deepest part, as shown in Fig. 14.

To include wear in the geometry, it is necessary to make some adaptations to
the stator radius. Considering that the fault starts at the angular position
$\theta=\theta_{s}$ and ends at $\theta=\theta_{f}$, the stator description
from the origin is defined as

$\displaystyle R_{o}^{*}=R_{o}+d_{\theta}$ (50)

where $d_{\theta}=$ $\begin{cases}0\text{,}&\mbox{se}\quad
0\leq\theta\leq\theta_{s}\text{,}\quad\theta_{f}\leq\theta\leq 2\pi\\\
d_{0}-F\left(1+\cos{\left(\theta-\pi/2\right)}\right)\text{,}&\mbox{se}\quad\theta_{s}<\theta<\theta_{f}\end{cases}$.

In $\theta_{s}$ and $\theta_{f}$, the wear depth is zero, so the location of
the edges can be defined as follows:

$\displaystyle\theta_{s}=\pi/2+\cos^{-1}{\left(d_{0}/F-1\right)}+\gamma$
$\displaystyle\theta_{f}=\pi/2-\cos^{-1}{\left(d_{0}/F-1\right)}+\gamma$ (51)

Element | Value
---|---
Stator radius $(R_{o})$ | $0.015\left[m\right]$
Gap $(F)$ | $9\cdot 10^{-6}\left[m\right]$
Length $(L)$ | $0.02\left[m\right]$
Load $(W)$ | $18.9\left[N\right]$
Viscosity $(\mu)$ | $0.1044\left[Pa.s\right]$
Rotation speed $(\omega)$ | $104.72\left[rad/s\right]$
Table 3: Numerical simulation for the worn bearing. Fonte: Adapted from [24]

With the geometry described as above, and data shown in Table 3, the results
for the pressure field in the bearing with wear are verified with [24]. For a
cylindrical bearing without wear, the pressure field is represented by Figs.
15 (a-b). Figs. 16 (a-b) illustrate the pressure field for the bearing with
maximum depth wear $d_{0}=10[\mu m]$ and an offset $\gamma=10^{\circ}\approx
0.17[rad]$. In Figs. 17 (a-b) $d_{0}=50[\mu m]$ is considered for the same
displacement and, finally, Figs. 18 (a-b) show the pressure for the same wear
with displacement $\gamma=30^{\circ}\approx 0.52[rad]$.

(a) Pressure distribution.

(b) Side view.

Figure 15: Bearing without wear.

(a) Pressure distribution.

(b) Side view.

Figure 16: Wear bearing ($d_{0}=10[\mu m]$, $\gamma=10^{\circ}$).

(a) Pressure distribution.

(b) Side view.

Figure 17: Wear bearing ($d_{0}=50[\mu m]$, $\gamma=10^{\circ}$).

(a) Pressure distribution.

(b) Side view.

Figure 18: Wear bearing ($d_{0}=50[\mu m]$, $\gamma=30^{\circ}$).

It is possible to see that, even for minor wears, the pressure field is
affected drastically. In addition to presenting two peaks, as occurred in the
elliptical bearing, the magnitude of the pressure is also changed. Fig. 19
shows that, as the depth $d_{0}$ increases, the shape of the pressure field
changes considerably.

Figure 19: Pressure along $\theta$ in the middle of the bearing with wear to
different depths $d_{0}$ ($\gamma=10^{\circ}$).

The depth of wear on the structure also influences the equilibrium position of
the rotor, as seen in Figs. 20. As the depth $d_{0}$ increases, the center of
the rotor moves away from the center of the stator, increasing the
eccentricity. In addition to the direction opposite to the center of the
stator, the path taken also differs in behavior when compared to the
cylindrical bearing without wear. As the wear depth increases, the center of
the rotor moves further down, describing a course close to a straight line.

(a) Rotor center.

(b) Variation of eccentricity.

Figure 20: Variation of the balance position in the bearing with wear for
different depths $d_{0}$ ($\gamma=10^{\circ}$).

The radial and tangential forces for the worn bearing develop a different
behavior, comparing to the elliptical bearing. Figure 21 indicates that the
radial component of the force increases rapidly at greater depths of wear,
while the value of the tangential component decreases.

Figure 21: Behavior of radial and tangential forces on the bearing with wear
to different depths $d_{0}$ ($\gamma=10^{\circ}$).

(a) Stiffness coefficients.

(b) Damping coefficients.

Figure 22: Dynamic bearing coefficients with wear for different depths $d_{0}$
($\gamma=10^{\circ}$).

Finally, with regard to the dynamic coefficients, Figs. 22 (a-b) describe the
behavior of stiffness and damping for different wear depths. There is an
oscillation in the stiffness coefficients, with all the terms distinct from
each other and only negative $k_{yx}$ values. For maximum depth wear up to
$d_{0}=30[\mu m]$, the cross term $k_{xy}$ is greater than the direct term
$k_{xx}$, indicating possibility of instability. From that value, $k_{xx}$
grows, distancing itself from the cross term, but then decreases again shortly
thereafter. The coefficients related to the vertical forces do not show very
relevant variations, the direct term being always much higher than the cross
term. The damping coefficients develop positive values in their direct terms,
and negative and similar values in their crossed coefficients.

Thus, it can be concluded that wear on the cylindrical bearing are variations
in geometry that directly affect the dynamic response of the system. The
stiffness coefficients show intervals in which the risk of instabilities is
more critical, which can serve to indicate operational conditions that should
be avoided.

###  7.4 Comparison between geometries

Nowadays, many rotating machines operate at high rotation speeds. This speed
influences the entire dynamic response of the bearing, from the equilibrium
position to the stiffness and damping coefficients. For this reason, it is
essential to understand how the model’s responses behave for different
rotation speeds.

Figure 23: Schematic drawing listing the three bearing configurations shown. Element | Value
---|---
Stator radius $(R_{o})$ | $0.015\left[m\right]$
Gap $(F)$ | $9\cdot 10^{-5}\left[m\right]$
Length $(L)$ | $0.02\left[m\right]$
Load $(W)$ | $100\left[N\right]$
Viscosity $(\mu)$ | $0.05449\left[Pa.s\right]$
Rotation speed $(\omega)$ | $50\leq\omega\leq 550\left[rad/s\right]$
Table 4: Numerical simulation for the comparison of bearing geometries.

In order to compare the results, the data shown in Table 4 is used, adapting
only the specifics of each geometry and varying the rotation speed in the
range of $50\leq\omega\leq 550\left[rad/s\right]$. For the elliptical bearing
the preload $m=0.4$ was considered and, for the bearing with wear, the maximum
depth $d_{0}=50[\mu m]$ with an angular displacement $\gamma=10^{\circ}$. The
schematic drawing that relates the three bearing configurations is illustrated
in Fig. 23, in which the geometries are indicated by the letters $C$, $E$ and
$D$, representing cylindrical, elliptical and worn bearings, respectively.

(a) Rotor center.

(b) Variation of eccentricity.

Figure 24: Balance position for different bearing configurations, varying the
rotation speed.

Figures. 24 (a-b) illustrate the variation of the equilibrium position for the
cylindrical, elliptical and worn bearing. In the three cases, the eccentricity
decreases with the increase of the rotation speed. However, it is possible to
observe the different paths. While the cylindrical bearing moves in a
semicircle-like path, the elliptical bearing generates a steeper curve, in
which most points are aligned with the horizontal axis. On the other hand, the
worn bearing has much higher eccentricities for all speeds. It is interesting
to highlight the difference between the bearing balance positions with and
without wear. Despite having very similar geometry, wear alters the position
of the rotor center and the way it behaves at different speeds.

Figure 25: Pressures along $\theta$ in the middle of the bearing for different
geometry configurations, varying the rotation speed.

As already shown, the pressure behavior is changed for each geometry. Fig. 25
demonstrates these changes in order to facilitate comparison. It is noticed
that, in addition to the new peak, the maximum pressure value is also
modified. The cylindrical bearing has the lowest pressure value, while the
elliptical reaches the greatest magnitude, even though the configuration with
the least eccentricity. It is also possible to notice that the wear causes an
increase in the pressure comparing to the cylindrical bearing, also modifying
the place where the pressure reaches its highest value.

Figure 26: Forces for different bearing configurations, varying the rotation
speed. The continuous line represents the radial component and the dashed line
represents the tangential component.

The radial and tangential components of the forces of the fluid film (Fig.
26), as well as the pressure fields, are different for the three cases
presented. In cylindrical and elliptical bearings, the radial force $N$
presents a clear decrease when the rotation speed increases, even reaching
negative values, for the case of the elliptical bearing. For wear, the
decrease of this component is moderate, with values much higher than the
tangential force $T$, which is also growing smoothly. The tangential force $T$
for the elliptical bearing is little influenced by the rotation speed. From
150 rad/s up, the this force is almost constant for the rotation speed range
analyzed. On the other hand, for cylindrical geometry, varying the rotation
speed has a higher impact on the forces.

(a) $k_{xx}$

(b) $k_{xy}$

(c) $k_{yx}$

(d) $k_{yy}$

Figure 27: Stiffness coefficients for different bearing configurations,
varying the rotation speed.

Finally, Figs. 27 (a-d) and 27 (a-d) show the variations of each term of the
stiffness and damping coefficients for the different geometries in the
rotation speed range adopted. Some characteristics stand out in these results
and are commented below.

Regarding the stiffness coefficients, it is noticed that, in the elliptical
bearing, the terms $k_{xy}$ and $k_{yy}$ have rapid growing up to
$\omega\approx 300[rad/s]$ and, then there is a drastic reduction in their
values, which does not occur in the other two geometries. The term $k_{xx}$
shows little variation with the rotation speed for the cylindrical bearing. It
decreases for the elliptical bearing and develops much higher values for the
wear geometry. The terms $k_{xy}$, $k_{yx}$ and $k_{yy}$ show similar growth
and decrease their values in cylindrical and worn bearings, with the greatest
difference between them in the crossover term $k_{yx}$.

(a) Direction $x$.

(b) Direction $y$.

Figure 28: Difference between the direct and cross terms of stiffness in each
direction.

In order to analyze the stability conditions in the different bearings, it is
necessary to compare the direct and crossed stiffness coefficients for the
same direction. Such comparison becomes more understandable in Figs. 28 (a-b),
which illustrate the differences $(k_{xx}-k_{xy})$ and $(k_{yy}-k_{yx})$ in
the applied rotation speed range. If these differences are negative, they mean
that the cross term has a higher value than the direct term, indicating
possibility of rotor instability. It is desired that these values are high and
positive. In this way, it is possible to notice that the differences, for the
three bearing configurations in the $x$ direction, tend to decrease. However,
the elliptical bearing presents a different behavior from the others, with the
cross term greater than the direct one for rotation speeds in the range
$100<\omega<300\left[rad/s\right]$. The difference in the terms in the $y$
direction is positive for the entire variation, for the three geometries, also
showing much less uniformity in the elliptical bearing.

(a) $c_{xx}$

(b) $c_{xy}$

(c) $c_{yx}$

(d) $c_{yy}$

Figure 29: Damping coefficients for different bearing configurations, varying
the speed.

Unlike what happens in stiffness, damping is not as sensitive to the presence
of faults. All coefficients develop close values for cylindrical and worn
bearings, except for the direct terms $c_{xx}$, which is different for lower
speeds, but end up getting closer as the speed increases. The elliptical
geometry results in damping coefficients that follow the same growth and
decrease patterns as the others, but always with higher values. The cross
terms $c_{xy}$ and $c_{yx}$ show increasing and negative values for all cases.
However, with higher rotation speeds, they become positive only on the
elliptical bearing. Direct terms decrease in all cases without reaching
negative results.

##  8 Conclusions

Considering the flow of a small thickness of lubricating oil film between two
cylinders, this work presents a review of the whole Lubrication Theory,
starting from the Navier-Stokes Eq. (2a) until reaching an elliptical partial
differential Eq. (19) that describes the pressure field in the annular space,
corresponding to the well-known Reynolds equation normally used in the
literature on the subject.

Through a dimensionless analysis, it was possible to promote a mapping of the
terms of the equations, ordering them by relevance in the physical phenomenon
for this type of situation. By revisiting this formulation, it was also
possible to choose a new and effective parameterization of geometry, which
enables a more detailed description of the problem, allowing to adapt the
modeling of this work to other contexts.

Without changing the proposed model, simulations were carried out in different
configurations of hydrodynamic bearings. This parameterization also allows the
model to be adapted in the future for many other geometries, including changes
not only in the tangential direction, but also in the axial direction. In
addition, the equations were written in cylindrical coordinates in order to
preserve the effects of curvature, which may be interesting for future works.

The pressure field was obtained numerically through the centered finite
differences method. This result was verified from its application in
simulations for cylindrical hydrodynamic bearings, comparing the responses of
this model with those obtained through classic approaches in the literature.

A new methodology was presented, different from the one commonly used in the
literature, to compute the stiffness and damping coefficients. This
alternative, through a Least Squares solution obtained by the Moore-Penrose
pseudo inverse, making the process computationally effective.

Another relevant result is the construction of the algorithm in Python,
developed by the authors of this research and used in all the simulations
presented. Through Fluid Flow it is possible to choose geometric and
operational parameters according to analyzes of the desired pressure
performance, forces, equilibrium position, dynamic coefficients, etc. In
addition, the low computational cost of the implemented modeling allows to
generate different simulations, which enables its use in machine learning
algorithms, a function already explored in Ross, which contains a stochastic
module to provide this type of analysis. This module also explores
uncertainties about the model presented in this work, such as, for example,
the impact of the equilibrium position on the dynamic coefficients.

Two other configurations of hydrodynamic bearings were also presented: the
elliptical and wear bearing. These new structures, obtained through changes
only in the geometric description of the problem, without changing the
modeling, portray the gain resulting from the chosen parameterization. In each
situation, the impact of these changes on the dynamic response of the system
was analyzed and, finally, comparisons were made between the three geometries,
exploring the changes in the responses for different rotation speeds.

The results obtained make clear the impact of the variation of the bearing
type on the fluid model’s responses. The elliptical bearing, introduced as an
alternative to correct the instabilities generated in the cylindrical bearing,
showed the high stiffness coefficients, which might lead to rotor instability.
As for the bearing with wear, it is clear that even small changes compromise
the entire response of the model. The pressure field changes dramatically, as
do the forces and stiffness. On the other hand, the worn does not seem to
influence the system’s damping so much.

It is concluded that the description of geometry is fundamental for the
success of the modeling. It is worth noting that the parameterization chosen
for the description of the problem allows such changes to be made naturally,
without the need to modify the model. Furthermore the proposed model may
constitute a tool for choosing the best geometry, in order to optimize the
operation of rotating machines within their specific contexts.

##  9 _Acknowledgements_

The code used for the simulations presented in this paper is part of
ROSS333https://github.com/ross-rotordynamics/ross, an open source library
written in Python for rotordynamic analysis, by [31].

## References

  * [1] W. M. Rankine, On the centrifugal force of rotating shafts, Van Nostrand’s Eclectic Engineering Magazine (1869-1879) 1 (7) (1869) 598, publisher: Center for Research Libraries.
  * [2] S. Dunkerley, On the Whirling and Vibration of Shafts., Proceedings of the Royal Society of London 54 (1893) 365–370.
  * [3] A. Foppl, Das problem der laval’schen turbinenwelle, Civilingenieur 41 (1895) 248–250.
  * [4] H. H. Jeffcott, The lateral vibration of loaded shafts in the neighbourhood of a whirling speed.—The effect of want of balance, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 37 (219) (1919) 304–314.
  * [5] D. Childs, Turbomachinery Rotordynamics: Phenomena, Modeling, and Analysis, John Wiley & Sons, 1993.
  * [6] A. Muszynska, Rotordynamics, CRC Press, 2005.
  * [7] R. Tiwari, Rotor Systems: Analysis and Identification, CRC Press, 2017.
  * [8] J. M. Vance, Rotordynamics of Turbomachinery, John Wiley & Sons, 1988.
  * [9] Y. Ishida, T. Yamamoto, Linear and Nonlinear Rotordynamics: A Modern Treatment with Applications, John Wiley & Sons, 2013.
  * [10] A. Sommerfeld, Zur hydrodynamischen theorie der schmiermittelreibung, Z. Math. und Physik 50 (1904) 97.
  * [11] H. W. Swift, The stability of lubricating films in journal bearings, Proceedings of the Institution of Civil Engineers 233 (1932) 289–322.
  * [12] W. Stieber, Das schwimmlager, V.D.I. (1933).
  * [13] F. W. Ocvirk, Short bearing approximation for full journal bearings, NACA TN2808. (1952).
  * [14] A. Cameron, L. Wood, The full journal bearing, Proceedings of the Institution of Mechanical Engineers 161 (1949) 59.
  * [15] O. Pinkus, Analysis of elliptical bearings, Transactions ASME 78 (jul 1956).
  * [16] G. C. Singhal, Computation methods for hydrodynamic problems (reynold’s equation), Cumputer-Aided Design 13 (3) (1981) 151–154.
  * [17] E. P. Pina, M. S. Carvalho, Three-dimensional flow of a newtonian liquid through an annular space with axially varying eccentricity, Journal of Fluids Engineering 128 (2) (2006) 223–231.
  * [18] S. Andrade, Asymptotic model for monophasic flow through progressive cavities pump, Master’s thesis, Pontifícia Universidade Católica do Rio de Janeiro, Rio de Janeiro (2008).
  * [19] L. H. M. Queiroz, Análise de Escoamento de Fluido no Interior de Componentes de Máquinas Rotativas, projeto de Graduação — Curso de Engenharia Mecânica da Escola Politécnica, Universidade Federal do Rio de Janeiro, Rio de Janeiro, 2018. (2018).
  * [20] L. San Andrés, Modern lubrication theory: Static load performance of plain journal bearings (notes 4) (2010).
URL http://oaktrust.library.tamu.edu/handle/1969.1/93197

  * [21] J. Frêne, D. Nicolas, B. Degueurce, D. Berthe, M. Godet, Static characteristics of journal bearings, in: Hydrodynamic Lubrication - Bearings and Thrust Bearings, Vol. 33 of Tribology Series, Elsevier, Amsterdan, 1990, pp. 113–141.
  * [22] B. J. Hamrock, Hydrodynamic journal bearings - analytical solutions, in: Fundamentals of fluid lubrication, NASA, Washington, 1991, pp. 219–232.
  * [23] Y. Ishida, T. Yamamoto, Flow-induced vibrations, in: Linear and nonlinear rotordynamics, Wiley Online Library, Weinheim, 2012, pp. 235–261.
  * [24] T. H. Machado, K. L. Cavalca, Modeling of hydrodynamic bearing wear in rotor-bearing systems, Mechanics Research Communications - Elsevier 69 (2015) 15–23.  doi:10.1016/j.mechrescom.2015.05.008.
  * [25] T. H. Machado, Avaliação de mancais hidrodinâmicos com descontinuidades geométricas, dissertação de Mestrado — Faculdade de Engenharia Mecânica, Faculdade Estadual de Campinas, Campinas, 2011. (2011).
  * [26] J. de Araújo Mota, Estudo da teoria de lubrificação com parametrização diferenciada da geometria e aplicações em mancais hidrodinâmicos., Master’s thesis, PUniversidade Federal do Rio de Janeiro, Rio de Janeiro (2020).
  * [27] R. W. Fox, A. T. McDonald, J. W. Mitchell, Introduction to Fluid Mechanics, 10th Edition, Wiley, 2020.
  * [28] R. L. Burden, J. D. Faires, Numerical Analysis, 9th Edition, Cengage Learning, Boston, 2010.
  * [29] L. San Andres, Modern lubrication theory, “static load performance of plain journal bearings”, notes 4, http://oaktrust.library.tamu.edu/handle/1969.1/93199, acesso em 08-02-2020 (2010).
  * [30] M. I. Friswell, J. E. T. Penny, S. D. Garvey, L. A. W., Dynamics of Rotating Machines, Cambridge University Press, 2010.
  * [31] R. Timbó, R. Martins, G. Bachmann, F. Rangel, J. Mota, J. Valério, T. Ritto, Ross - rotordynamic open source software, Journal of Open Source Software 5 (48) (2020) 2120.  doi:10.21105/joss.02120.
URL https://doi.org/10.21105/joss.02120

  * [32] R. A. Ostayen, A. Beek, Thermal modelling of the lemon-bore hydrodynamic bearing, Tribology International - Elsevier 42 (2009) 23–32.  doi:10.1016/j.triboint.2008.05.01.
  * [33] D.R.Garner, C.S.Lee, F.A.Martin, Stability of profile bore bearings: influence of bearing type selection, Tribology International - Elsevier 13 (5) (1980) 204–210.  doi:10.1016/j.triboint.2008.05.01.
  * [34] L. H. J. Machado, Análise Experimental do Comportamento Dinâmico de Mancais Radiais Hidrodinâmicos de Geometria Fixa, dissertação de Mestrado — Programa de Pós-Graduação em Engenharia Mecânica, Universidade Federal de Minas Gerais, Minas Gerais, 2006. (2006).
  * [35] K. F. Dufrane, J. W. Kannel, T. H. McCloskey, Wear of steam turbine journal bearings at low operating speeds, J. of Lubrication Tech 105 (1983) 313–317.  doi:10.1115/1.3254599.
  * [36] H. Hashimoto, S. Wada, K. Nojima, Performance characteristics of worn journal bearings in both laminar and turbulent regimes. part i: Steady-state characteristics, ASLE Transactions 29 (2008) 565–571.  doi:10.1080/05698198608981721.