Sage Journals: Discover world-class research

Abstract

Hydropower rotors and pumps have the specificity to be oriented vertically, meaning that the bearing forces have to be evaluated at each time-step depending on the position of the rotor for dynamical analyses. If the bearing forces cannot be evaluated analytically, a suitable numerical method should be used to calculate the pressure distribution over the bearing domain. This process can be computationally expensive as it should be performed for each discrete time-step. As a result, a comparison between the spectral method, the finite difference method, and the finite element method is performed to investigate which method is more adapted to dynamical analysis of the bearing. It is observed that the spectral method has the advantage of having a reasonable simulation time for any eccentricity magnitude with a moderate number of interpolation points. However, this method should be restricted to simple bearing models such as plain bearings or multilobe bearings due to the advantage of finding a global numerical solution directly on the entire bearing/pad domain.

Keywords

Vertical machine multilobe bearing convergence

Introduction

In horizontal turbomachines, the stiffness and damping properties of the bearings are calculated around a stationary position due to dead weight. On the contrary, vertical machines such as pumps and water turbines do not operate around a static equilibrium position but wobbles around the geometric center of the bearing due to unbalance forces. As a result, the bearing forces have to be calculated by solving the Reynolds equation at each time-step. Since rotordynamical simulations are usually performed for more than a hundred revolutions to reach steady state or for analyses of run-ups and shut-downs, there is a need to reduce the computational time in order to evaluate the system dynamics as function of design parameters such as oil viscosity, bearing type, or unbalance load in the case of vertical machines. As few complete rotordynamic codes are available on the market, it is important to develop functions to evaluate dynamical properties of bearings and rotors that could be of use for stand-alone computers in order to obtain a reasonable computational time.

Reynolds equation can be solved using different numerical methods. The oldest and most common way is to use the finite difference method (FDM). Cardinali et al.¹ used this method in a 160-MW Francis type hydro-unit to describe the nonlinear force model of tilting-pad bearing models. Another way to solve the Reynolds equation is to use the finite element method (FEM). Ma and Zhang² adopted this method to evaluate the nonlinear oil film forces of the tilting-pad guide bearing in order to simulate the dynamics of a hydro-unit supported by a generator and a turbine guide bearing. The main concern with FDM and FEM is the computational cost required to obtain accurate results. An alternative choice to these two methods is the spectral method described by Boyd³ and Trefethen.⁴ The spectral method is widely used in fluid dynamics, but its use to solve Reynolds equation in bearings has been limited to the work of Raad and Karageorghis⁵ and Gantasala et al.⁶ Fundamentally, the FEM and spectral method are based on the same assumption: the solution of the differential equation is approximated by a set of trial functions satisfying the boundary conditions, and the residuals are minimized to obtain a solution. In the FEM case, the domain is split in smaller elements where the solution is approximated by local low-order trial functions. On the contrary, the spectral method uses global higher order trial functions to approximate the solution over the entire computational domain. The advantage of using the spectral method is the high accuracy of the solution due to the choice of global basis functions.

In this study, multilobe bearings with three and four lobes will be investigated since they are of common use in vertical machines. Transient simulations using multilobe bearings with different numerical schemes on rigid rotors have been widely investigated in the 1980s by Allaire et al.⁷ and Li et al.⁸ and more recently by Rao et al.,⁹ but not by using the spectral method. As a result, a comparison of the convergence of the FDM, FEM, and spectral method is performed with different mesh sizes. A focus will be made on the required computational time to perform one calculation of the fluid-film forces. Moreover, an application to a dynamic problem of a vertical machine will be qualitatively investigated. The numerical method used in this article can then be extended to other bearing types such as arc bearings or tilting-pad journal bearings.

Model description

Bearing modeling

The bearing used in this study is a three-lobe bearing with grooves. The geometry of the bearing is available in Figure 1 with the corresponding parameters in Table 1. The fluid-film pressure in the bearing is found from the solution of the Reynolds equation

\frac{1}{R^{2}} \frac{\partial}{\partial θ} (h^{3} \frac{\partial p}{\partial θ}) + \frac{\partial}{\partial z} (h^{3} \frac{\partial p}{\partial z}) = 6 μ Ω \frac{\partial h}{\partial θ} + 12 μ \frac{\partial h}{\partial t}

(1)

Figure 1.

Description of the multilobe bearing.

Table 1.

Geometric properties of the multilobe bearing.

Symbol	Item	Value	Unit
$μ$	Dynamic viscosity	$0.89 \times 10^{- 4}$	Pa s
$R$	Journal radius	$1.905 \times 10^{- 2}$	m
$L$	Bearing width	$3.81 \times 10^{- 2}$	m
$α$	Offset	$0.6$	–
$C_{b}$	Bearing clearance	$1.676 \times 10^{- 4}$	M
$m_{p}$	Preload factor	$0.476$	–
$θ_{p}$	Lobe angle	$110$	°
$N_{lobes}$	Number of lobes	$3$	–

where $θ$ represents the circumferential coordinate and $z$ the axial coordinate. The oil film thickness $h$ of the multilobe bearing can be derived from geometrical relationships in Figure 1. Using Lund’s convention, it can be found that

\begin{matrix} h (θ) = C_{p} + e_{x} \cos (θ) + e_{y} \sin (θ) \\ - (C_{p} - C_{b}) \cos (θ - θ_{0} - α θ_{p}) \end{matrix}

where $C_{p} = R_{p} - R$ is the pad clearance; $C_{b} = R_{b} - R$ the bearing clearance; $e_{x}$ and $e_{y}$ the positions of the journal bearing center $O_{j}$ in the X and Y direction, respectively; $θ_{0}$ the leading edge angle of the lobe; and $θ_{p}$ the arc angle. The lobe is positioned so that the minimum clearance occurs at the angle $θ_{c}$ , meaning that the offset is defined by $α = (θ_{c} - θ_{0}) / θ_{p}$ . The film thickness is not dependent on $z$ since there is no misalignment in the model. Equation (1) is integrated numerically to obtain the pressure distribution on each lobe. In this study, the Gumbel boundary condition is used, where the negative pressures are replaced by the atmospheric pressure $P_{a} \approx 0$ . Moreover, the pressure at the leading edge and trailing edge of each lobe should also be equal to the atmospheric pressure as well as the pressure at each end of the bearing in the axial coordinate z, that is

{\begin{matrix} \forall z p (θ_{0}, z) = p (θ_{0} + θ_{p}, z) = 0 \\ \forall θ p (θ, - L / 2) = p (θ, + L / 2) = 0 \end{matrix}

(3)

By integrating along the circumferential and axial directions, the forces can be found in the fixed coordinate system using

{\begin{matrix} F_{x} = \sum_{k = 1}^{N_{lobes}} \int_{θ_{0}^{k}}^{θ_{p}^{k}} \int_{- L / 2}^{L / 2} p (θ, z) R \cos (θ) d θ dz \\ F_{y} = \sum_{k = 1}^{N_{lobes}} \int_{θ_{0}^{k}}^{θ_{p}^{k}} \int_{- L / 2}^{L / 2} p (θ, z) R \sin (θ) d θ dz \end{matrix}

(4)

Numerical methods

In this section, a short description of the three numerical methods used to solve the pressure distribution is available, with more details given for the spectral method since it is the numerical scheme of interest.

FDM

One of the simplest methods to solve the Reynolds equation is the FDM, which is based on a finite difference approximation of the derivatives. By applying it to equation (1), we obtain the following discrete form

\begin{matrix} \frac{1}{R^{2}} [\frac{h_{i + 1 / 2, j}^{3} (p_{i + 1, j} - p_{i, j}) - h_{i - 1 / 2, j}^{3} (p_{i, j} - p_{i - 1, j})}{Δ Θ^{2}}] \\ + [\frac{h_{i, j + 1 / 2}^{3} (p_{i, j + 1} - p_{i, j}) - h_{i, j - 1 / 2}^{3} (p_{i, j} - p_{i, j - 1})}{Δ Z^{2}}] = f_{i} \end{matrix}

(5)

The pressure can be found

e p_{i, j} = a p_{i + 1, j} + b p_{i - 1, j} + c p_{i, j + 1} + d p_{i, j - 1} - f_{i}

(6)

where

{\begin{matrix} a = 1 / (R Δ Θ)^{2} {[0.5 (h_{i + 1, j} + h_{i, j})]}^{3} \\ b = 1 / (R Δ Θ)^{2} {[0.5 (h_{i - 1, j} + h_{i, j})]}^{3} \\ c = 1 / (Δ Z)^{2} {[0.5 (h_{i, j + 1} + h_{i, j})]}^{3} \\ d = 1 / (Δ Z)^{2} {[0.5 (h_{i, j - 1} + h_{i, j})]}^{3} \end{matrix}

(7)

with $e = a + b + c + d$ and

\begin{array}{l} f_{i} = 6 μ [Ω - e_{x} \sin (i Δ Θ) + e_{y} \cos (i Δ Θ) + (C_{p} - C_{p}) \sin (i Δ Θ - θ_{c})] \\ + 12 μ [{\dot{e}}_{x} \cos (i Δ Θ) + {\dot{e}}_{y} \sin (i Δ Θ)] \end{array}

(6)

represents the source term on the right side of Reynolds equation. Since there is no misalignment in the axial direction, one can notice that $h_{i, j + 1} = h_{i, j} = h_{i, j - 1}$ . The pressure distribution is solved using the Gauss Seidel iterative scheme

\begin{matrix} p_{i, j}^{k + 1} = (1 - α) p_{i, j}^{k} \\ + α (a p_{i + 1, j} + b p_{i - 1, j} + c p_{i, j + 1} + d p_{i, j - 1}) / e \end{matrix}

where $α$ is the over-relaxation coefficient. In this case, the optimum value of the over-relaxation coefficient follows the description of Singhal,¹⁰ where its value is function of the number of mesh points in the axial and circumferential directions as well as the radius to length ratio

{\begin{matrix} α = 2 (\frac{\sqrt{1 - (1 - μ^{2})}}{μ^{2}}) \\ μ = 1 - \frac{π^{2}}{2} \frac{4 + {(π D / L)}^{2}}{N^{2} + {(M π D / L)}^{2}} \end{matrix}

(10)

The pressure distribution is iterated until a convergence condition is reached

\sum \sum ∥ p_{i, j}^{k + 1} - p_{i, j}^{k} ∥ / \sum \sum p_{i, j}^{k + 1} \leq ϵ

(11)

The integration of the pressure distribution to find the bearing forces is done using the Simpson numerical integration.

FEM

In the FEM, the main concept is to split the domain in small elements where the solution is approximated by local low-order trial functions. Some details about this method applied to bearings can be found in Allaire et al.¹¹ as well as Booker and Huebner.¹² In our case, the fluid domain is discretized with 9-node quadrilateral elements to obtain the global matrices corresponding to the pressure fluidity component $K_{p}$ , the shear fluidity component $K_{Ux}$ , and the squeeze component $K_{\overset{\cdot}{h}}$ . The isoparametric formulation is used to obtain the stiffness matrix for each element. Due to the symmetry of the bearing along its length, only half of the bearing is considered as the computational domain. The pressure distribution is then obtained by solving the following equation

K_{p} {p} = [K_{U_{x}} {U_{x}} + K_{\overset{\cdot}{h}} {\overset{\cdot}{h}}]

(12)

and applying the appropriate boundary conditions.

Pseudo-spectral method

The idea of the spectral method is to assume the solution as a linear combination of orthogonal polynomials. The differential equation can be written in a simplified way

L u (x) = s (x) for x \in V

(13)

with the boundary condition $B u (y) = 0$ for $y \in \partial V$ , where $L$ and $B$ are differential operators. This boundary condition justifies the use of the spectral method to bearings since the atmospheric pressure is usually assumed $P_{a} \approx 0$ on all the edges in comparison with the pressure magnitude within the computational domain. The solution is written in terms of orthogonal polynomials as

u (x) \approx \hat{u} (x) = \sum_{k = 0}^{N} u_{k} ϕ_{k} (x)

(14)

where $ϕ_{k} (x)$ can be either Fourier functions or Chebyshev polynomials. The residual is obtained by substitution of equation (14) in equation (13) and is given by

R_{N} (x) = L \hat{u} (x) - s (x)

(15)

The solution of the differential equation will be a function $\hat{u} (x)$ which should satisfy the boundary conditions and make the residual as small as possible. The main method used to identify the solution with minimum residual is the weighted residuals method. In this method, the inner product of the residual and the weight function are made equal to zero, that is

〈 R_{N} (x), w (x) 〉 = \int (L \hat{u} (x) - s (x)) w (x) dx = 0

(16)

\Leftrightarrow \int L \hat{u} (x) w (x) dx = \int s (x) w (x) dx

(17)

The collocation method is implemented where the residual is made equal to zero at specific points $x_{n}$ . In the pseudo-spectral method (PSM), the weight functions are chosen as Dirac delta functions at the specific points, that is, $w (x) = δ (x - x_{n})$ . Equation (17) finally becomes

\sum_{k = 0}^{N} \int L (u_{k} ϕ_{k} (x)) δ (x - x_{n}) dx = \int s (x) δ (x - x_{n}) dx

(18)

\Leftrightarrow \sum_{k = 0}^{N} L (ϕ_{k} (x_{n})) u_{k} = s (x_{n})

(19)

Since there is $N$ unknowns in the function $\hat{u} (x)$ , we get N algebraic equations corresponding to $N$ weight functions at the specified points ${x_{1}, x_{2}, . . ., x_{n}}$ . Equation (19) is usually written in the form

Du = s

(20)

where $D$ is referred as the differentiation matrix. Since the bearing geometry is split in several lobes, the periodicity in the circumferential direction is lost. As a result, the Chebyshev spectral method is used for both the circumferential and axial direction. In this method, the points are not equally spaced but defined by the Gauss–Lobatto points $x_{i} = \cos (π i / N)$ to avoid oscillation at the edge due to Runge’s phenomenon. The derivation of the differentiation matrix $D$ for the Chebyshev spectral method can be found in Trefethen.⁴ Since the Chebyshev spectral method is applied to a domain defined in the interval $[- 1, 1]$ , the following coordinate transformation $z = (L / 4) (\tilde{z} + 1)$ is used for the axial coordinate to obtain the interval $[0, L / 2]$ (as only half the length of the bearing is considered due to symmetry).

As the Chebyshev polynomials are used both in the circumferential and axial directions, the pressure is approximated with $N + 1$ and $M + 1$ non-uniformly spaced grid points as

p (θ, z) = \sum_{n = 0}^{N} \sum_{m = 0}^{M} p_{nm} T_{n} (θ) T_{m} (z)

(21)

Reynolds equation (1) can be written in a tensor form by substitution of the differentiation matrices as described in Trefethen⁴

\begin{matrix} (D_{θ} \otimes I_{M + 1}) (A \otimes I_{M + 1}) (D_{θ} \otimes I_{M + 1}) \\ + (I_{N + 1} \otimes D_{z}) (B \otimes I_{M + 1}) (I_{N + 1} \otimes D_{z}) {p} = {f} \end{matrix}

(22)

where $I_{M + 1}$ and $I_{N + 1}$ are the identity matrices of size $M + 1$ and $N + 1$ , respectively, and

{\begin{matrix} A_{jl} = {(\frac{h^{3}}{R^{2}})}_{θ = θ_{j}} δ_{jl} \\ B_{jl} = {(\frac{16 h^{3}}{L^{2}})}_{θ = θ_{j}} δ_{jl} \end{matrix}

(23)

Results

Simulation procedure

Comparison study

To study the accuracy of the spectral method, the forces resulting from the pressure distribution will be compared with the FEM with quadrilateral elements and the FDM. The comparison is done for several mesh sizes, and the evaluation of the convergence is performed as well as the simulation time depending on the mesh size. An overview of the different types of mesh associated with each numerical method is available in Figure 2.

Figure 2.

Example of mesh grids of size $8 \times 4$ for the different numerical methods: (a) grid for the finite difference method, (b) mesh for the finite element method with 9-node quadrilateral elements, and (c) grid for the spectral method with Gauss–Lobatto points.

Dynamic analysis

A dynamic analysis of an unbalanced rigid rotor set in the horizontal and vertical position is also performed. The model is very simplified since it corresponds to a rigid point mass moving in a bearing. However, this type of model is often used to evaluate the nonlinear dynamic behavior of rotors induced by the bearings using analytical methods for short bearings as in Adiletta et al.¹³ or numerical methods in Meybodi et al.¹⁴ Fundamentally, the way to include the bearing forces in a rigid rotor and a full turbine is similar. The main difference is that the rotor elastic equations have to be coupled with the Reynolds equation, as in De Castro et al.¹⁵ and Jing et al.¹⁶ who evaluated oil whirl and whip on a vertical rotor and continuous model, respectively.

Numerical comparison

The parameters of the bearings used for simulating Reynolds equation are given in Table 1. The eccentricities at which the calculation are performed correspond to a percentage of the bearing clearance $C_{b}$ . For instance, for an eccentricity of 0.1 at an angle of 45°, one has the following displacements $e_{x} = 0.1 C_{b} \cos (π / 4)$ and $e_{y} = 0.1 C_{b} \cos (π / 4)$ , while the velocities are fixed with the following values ${\overset{\cdot}{e}}_{x} = 0.1 Ω e_{x}$ and ${\overset{\cdot}{e}}_{y} = 0.1 Ω e_{y}$ . The convergence of the three numerical schemes is verified for different eccentricities $ε = [0.1; 0.3; 0.5; 0.7; 0.9]$ at the same attitude angle. The mesh sizes are specified by $N \times M$ , where $N$ is the number of nodes in the circumferential direction and $M$ the number of nodes in the axial direction. The average simulation time is also given for one mesh size, which will of course depend on the capacities of the computer. In our case, we use a computer equipped with an Intel Core i72620M processor with a frequency of 2.70 GHz, while the simulations are performed using the MATLAB R2015b software. A code color has been used to make Tables 2 and 3 more visible. The red color indicates that the method has not fully converged, while the blue color corresponds to a value of a force which is less than 1% of the final converged value. The 1% convergence threshold is arbitrary and could be modified to obtain a better (or worse) accuracy, while the converged value corresponds to the force calculated using FEM for a large number of elements since it is the most accurate method as observed afterwards.

Table 2.

Comparison of fluid-film forces for different eccentricities and mesh sizes using the finite difference method.

FDM		Mesh size
Journal position		$8 \times 6$		$12 \times 5$		$18 \times 12$		$24 \times 15$		$36 \times 20$
Eccentricity	Attitude angle (°)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)
0.1	45	−1.39	−8.97	−1.54	−8.77	−1.57	−9.40	−1.57	−9.44	−1.57	−9.49
0.3	45	−1.12	−32.1	−1.17	31.4	−1.22	−33.4	−1.22	−33.6	−1.24	−33.8
0.5	45	11.1	−75.1	10.6	−73.7	11.2	−78.9	11.3	−79.4	11.4	−79.8
0.7	45	64.5	−193	63.0	−191	67.6	−206.5	68.0	−208.1	68.3	−209.1
0.9	45	526.7	−930	535.9	−959.1	607.0	−1079.6	616.1	−1097	626.2	−1113
Average time (s)		0.0093		0.0095		0.0286		0.0470		0.1132

Table 3.

Comparison of fluid-film forces for different eccentricities and mesh sizes using the finite element method.

FEM quad. elements		Mesh size
Journal position		$8 \times 6$		$12 \times 5$		$18 \times 12$		$25 \times 15$		$36 \times 20$
Eccentricity	Attitude angle (°)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)
0.1	45	−1.58	−9.51	−1.56	−9.48	−1.57	−9.52	−1.57	−9.52	−1.57	−9.53
0.3	45	−1.24	−33.8	−1.24	−33.8	−1.23	−33.9	−1.24	−33.9	−1.24	−33.9
0.5	45	11.5	−79.9	11.4	−79.9	11.4	−80.2	11.4	−80.2	11.4	−80.2
0.7	45	69.0	−210.1	69.0	−209.9	69.0	−210.8	69.0	−210.8	69.0	−210.9
0.9	45	627.7	−1116	630.0	−1119	635.2	−1126	635.5	−1127.5	635.7	−1128.0
Average time (s)		0.129		0.161		0.560		0.952		1.985

The simulation for FDM in Table 2 shows that the convergence is reached for a low number of nodes $(18 \times 12)$ for an eccentricity of 0.1 corresponding to a computational time of 0.0286. For the FEM in Table 3, the convergence is even reached for a lower number of elements for the same eccentricity of 0.1, but the simulation time is five times greater than FDM. Concerning the PSM in Table 4, it shows that the convergence is also attained for a mesh size of $18 \times 12$ , for a corresponding time of 0.040 which is still slower than the FDM. However, when the eccentricity is increased, the FDM takes a longer time to converge and needs a large increase of number to reach convergence. As observed for the $36 \times 20$ mesh size, the convergence is not still reached as the final force in the x-direction 626.2 N instead of a value of 635 N for the FEM and PSM.

Table 4.

Comparison of fluid-film forces for different eccentricities and mesh sizes for the pseudo-spectral method.

Pseudo-spectral		Mesh size
Journal position		$8 \times 6$		$12 \times 5$		$18 \times 12$		$25 \times 15$		$36 \times 20$
Eccentricity	Attitude angle (°)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)	$F_{x}$ (N)	$F_{y}$ (N)
0.1	45	−1.53	−9.39	−1.49	−9.35	−1.56	−9.49	−1.56	−9.50	−1.57	−9.51
0.3	45	−1.20	−33.45	−1.08	−33.42	−1.21	−33.8	−1.22	−33.8	−1.22	−33.9
0.5	45	10.9	−78.9	11.5	−79.4	11.4	−80.0	11.4	−80.1	11.4	−80.1
0.7	45	60.9	−199.0	68.5	−209.0	68.9	−210.5	68.9	−210	68.9	−210.6
0.9	45	293.4	−660.8	584.7	−1035	635.0	−1125	635.3	−1126	635.2	−1126
Average time (s)		0.0120		0.0125		0.040		0.106		0.530

As seen in Table 3, the FEM results show that the convergence of the force is obtained even for a low mesh size. Indeed, for a $8 \times 6$ size, the force in the x-direction is 69 N for an eccentricity of 0.7. An appropriate size mesh for one lobe using the FEM is $18 \times 12$ for any eccentricity. However, the simulation time is around $0.56 s$ for this particular size, so that the total computational time for a dynamic analysis will be too long since the forces should be calculated at each time-step. Concerning the PSM, convergence is not achieved for a low number of interpolation points, but the convergence is increased in a fast way by adding a few interpolation points. The main advantage is that the simulation time of $0.0286 s$ is reasonable at all eccentricities for a mesh size of $18 \times 12$ .

As a result, the FDM is more appropriate for low eccentricities due to the fast calculation speed, but it is less useful for large eccentricities as the number of nodes has to be high to fully represent the pressure distribution, resulting in an increase of simulation time. As seen in Figure 3(a), the convergence time is always lower for FDM when compared with the same number of elements. However, the convergence for an eccentricity of 0.9 is not attained for a mesh size $40 \times 20$ , so that the number of elements has to be increased significantly. In Figure 3(b), one can observe that the convergence is obtained rather fast for the PSM with a mesh size of $16 \times 8$ , which is similar to the FEM. However, the simulation time of FEM is far greater than PSM for this mesh size. The FEM achieves a good convergence even for a low number of elements, but the simulation time is substantially larger than for the FDM due to the assembly process. However, the FEM should be more suitable when the geometry of the bearing is more complex in comparison with the FDM. For the PSM, the convergence is not acceptable for a low number of elements, but the accuracy increases rapidly with regard to the number of interpolation points, so that a good compromise can be found between speed and accuracy for all eccentricities using this method. A qualitative summary of these results is given in Table 5 by indicating if the method is suitable depending on the computational time, the convergence of the forces, and the complexity of the geometry and implementation. Since the concept of PSM is to find a global solution on one single domain, it is more efficient to use on bearings with simple geometries such as plain bearings or multilobe bearings. When a bearing is split in numerous lobes, recess domains, and pressurized sections, the use of PSM is not relevant any longer as it loses the fundamental property of finding the solution on the entire computational domain. In that case, the use of FEM is more appropriate instead.

Figure 3.

(a) Simulation time as function of the number of elements and (b) evaluation of the convergence of the bearing force in the x-direction for the different methods. The simulations are performed for an eccentricity of 0.9.

Table 5.

Qualitative summary of the advantages and inconveniences of the three different numerical methods.

Mesh	Comp. time			Convergence			Geometry
	Coarse	Medium	Fine	Coarse	Medium	Fine	Simple	Complex
FDM	++	+	o	−	o	+	++	−−
FEM	o	−−	−−	+	++	++	o	++
PSM	++	+	−−	−−	+	+	+	−

FDM: finite difference method; FEM: finite element method; PSM: pseudo-spectral method; ++: very good; +: good; o: acceptable; −: poor; −−: unsatisfactory.

Dynamic application

The equation of motion of a rigid rotor supported by fluid bearings is simply given by

{\begin{matrix} m \overset{\cdot\cdot}{x} = F_{x} (x, y, \overset{\cdot}{x}, \overset{\cdot}{y}) + me Ω^{2} \cos (Ω t) \\ m \overset{\cdot\cdot}{y} = F_{y} (x, y, \overset{\cdot}{x}, \overset{\cdot}{y}) + me Ω^{2} \sin (Ω t) \end{matrix}

(24)

where $m$ is the rotor mass; $e$ is the eccentricity, that is, the distance between the center of mass and the geometric center of the rotor; $Ω$ is the rotating speed; and $F_{x}$ and $F_{y}$ are the bearing forces calculated with equation (1) using the spectral method. If the rotor is set horizontally, the static weight $- mg$ of the rotor should be included on the right-hand side of the second equation in equation (24). For this dynamic case, the equation of motion is written in the state-space form and solved using a fifth-order Runge–Kutta integration with adaptive time-step solver. Figures 4 and 5 display the orbit simulation at different speeds and eccentricities for the bearing described in Table 1. Figure 4(a) and (b) shows the first 10 periods of simulation for a speed of 7000 r/min with an unbalance eccentricity of 0.2 and 0.8, respectively. Figure 4(a) displays a transient behavior with a period two motion similar to Falkenhagen et al.¹⁷ for the same parameters, while Figure 4(b) shows a typical triangular-shaped orbit due to the variation of stiffness depending on the position of the rotor with regard to the three lobes. According to Tables 2 –4, it should be more appropriate to use FDM for the unbalance eccentricity of 0.2, while the PSM should be more suitable for the unbalance eccentricity of 0.8. However, using a mesh size of $12 \times 8$ for both methods, it is observed that the simulation time is 9.42 s for FDM and 7.75 s for PSM for an eccentricity $e = 0.2$ , while the simulation time is 52.40 s for FDM and 42.76 s for PSM for an eccentricity $e = 1$ . Even though a few elements have been used for the FDM for the higher unbalance case, the orbits of the rotor are very similar in shape and magnitude for the two methods. Hence, it is possible to reduce the number of elements for PSM and FDM to perform faster simulations and to still obtain reasonable results depending on the accuracy one wants to achieve in rotordynamics. Figure 5(a) shows also the orbit plots of the rotor for an unbalance eccentricity $e = [0.2; 0.4; 0.6; 0.8; 1]$ , but at a speed of 2000 r/min. Even for that case, the nonlinear response is similar to Figure 4(b) but with a smaller amplitude since the unbalance load is scaled with the square of the speed as seen in equation (24). For the smallest eccentricity $e = 0.2$ , the orbit is reasonably circular since the difference of stiffness on the lobe and between the lobe is small for small eccentricities. By increasing the eccentricity, the orbit shape becomes more triangular, but the maximum displacement is not proportional to the increase of unbalance eccentricity due to the stiffening of the bearing. Figure 5(b) shows the orbit response for the same conditions as Figure 5(a) except that the weight of the rotor is taken into consideration, meaning that the rotor rotates around a static equilibrium position when compared with the vertical case. Figure 6 displays a similar simulation, but with a four-multilobe bearing with a lobe angle of 72°, with the other bearing parameters kept the same. The simulation shown in Figure 6(a) is performed for 200 periods, while plotting the steady-state response after 100 periods. It shows a typical quasiperiodic motion that can be observed for some parameter value combination. The orbit in Figure 6(b) has the same type of behavior as for Figure 4(b) due to the variation of load along the bearing, but with a square shape due to the presence of four lobes instead of three. In this case, the simulation time for 10 periods of revolution is 84.22 s for the PSM with a mesh size of $12 \times 8$ . It is approximately two times slower with only one more pad, meaning that the dynamic simulation needs more steps to reach dynamic equilibrium which increase the total computational due to the repetition of bearing force calculations. The user could decrease the relative tolerance of the numerical integration of the equations of motion to obtain faster results, but it may affect the final solution accuracy.

Figure 4.

Orbit plots of the rigid rotor for dynamic simulations: (a) vertical rotor at $Ω = 7000 r / \min$ and $e = 0.2$ and (b) vertical rotor at $Ω = 7000 r / \min$ and $e = 0.8$ .

Figure 5.

Orbit plots of the rigid rotor for dynamic simulations: (a) vertical rotor at $Ω = 2000 r / \min$ and $e = [0.2; 0.4; 0.6; 0.8; 1]$ and (b) horizontal rotor at $Ω = 2000 r / \min$ and $e = [0.2; 0.4; 0.6; 0.8; 1]$ .

Figure 6.

Orbit plots of the rigid rotor with four-lobe bearings: (a) vertical rotor at $Ω = 7000 r / \min$ and $e = 0.2$ and (b) vertical rotor at $Ω = 7000 r / \min$ and $e = 0.8$ .

Conclusion

The aim of this article was to investigate numerical methods for solving the Reynolds equation in journal bearings. Using the geometry of a multilobe bearing, it is observed that the choice of a specific numerical method should be done with regard to the geometry of the bearing, the eccentricity, and the computational time. For instance, the FDM is appropriate for small eccentricities since the Central Processing Unit time is faster than the other methods due to a reasonable number of elements to reach convergence. However, the accuracy of this method decreases when the eccentricity becomes larger. On the contrary, the FEM has a good convergence for a low number of elements at any eccentricity, but the drawback is the simulation time which is expensive due to the assembly process. However, this method is more suitable for complex bearing geometry. The PSM also has the advantage of having a good convergence for a moderate number of elements at all eccentricities, with a simulation about 10 times faster than FEM. This method is suitable for simple journal geometries with a few number of lobes, which is frequently the case for vertical machines. Hence, the dynamic simulations of vertical machines can be performed in a faster and more accurate manner by applying the PSM to calculate the bearing forces at each time-step. However, the FDM is still an appropriate numerical method to use due to its programming simplicity and computational speed for most of the eccentricities, which may be the reason it is still widely used in rotordynamics software. Depending on the type of dynamic problem investigated, the user can choose the most appropriate method between FDM, FEM, and PSM if one wants to decrease computational time or increase the accuracy of the results.

Footnotes

Appendix 1

Academic Editor: Aditya Sharma

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The research presented was carried out as a part of “Swedish Hydropower Center – SVC”. SVC has been established by the Swedish Energy Agency, Elforsk and Svenska Krafnät together with Luleå University of Technology, The Royal Institute of Technology, Chalmers University of Technology and Uppsala University ().

References

Cardinali

Nordmann

Sperber

. Dynamic simulation of non-linear models of hydroelectric machinery. Mech Syst Signal Pr 1993; 7: 29–44.

L-F

Zhang

X-Z

. Numerical simulation of nonlinear oil film forces of tilting-pad guide bearing in large hydro-unit. Int J Rotat Mach 2000; 6: 345–353.

Boyd

. Chebyshev and Fourier spectral methods. Mineola, NY: Dover Publications, 2001.

Trefethen

. Spectral methods in MATLAB. Philadelphia, PA: Society for Industrial and Applied Mathematics, 2000.

Raad

Karageorghis

. A Chebyshev spectral collocation method for the solution of the Reynolds equation of lubrication. J Comput Phys 1993; 106: 42–51.

Gantasala

Praveen Krishna

Sekhar

. Dynamic analysis of rotors supported on journal bearings by solving Reynolds equation using pseudospectral method. Cham: Springer, pp.1009–1019.

Allaire

Choy

. Transient unbalance response of four multilobe journal bearings. J Lubr Technol 1980; 102: 300–306.

Choy

Allaire

. Stability and transient characteristics of four multilobe journal bearing configurations. J Lubr Technol 1980; 102: 291–298.

Rao

Biswas

Athre

. A methodology for dynamic coefficients and nonlinear response of multi-lobe journal bearings. Tribol T 2001; 44: 111–117.

10.

Singhal

. Computation methods for hydrodynamic problems (Reynold’s equation). Comput Aided Design 1981; 13: 151–154.

11.

Allaire

Nicholas

Gunter

. Systems of finite elements for finite bearings. J Lubr Technol 1977; 99: 187–194.

12.

Booker

Huebner

. Application of finite element methods to lubrication: an engineering approach. J Lubr Technol 1972; 94: 313–323.

13.

Adiletta

Guido

Rossi

. Chaotic motions of a rigid rotor in short journal bearings. Nonlinear Dynam 1996; 10: 251–269.

14.

Meybodi

Mohammadi

Bakhtiari-Nejad

. Numerical analysis of a rigid rotor supported by aerodynamic four-lobe journal bearing system with mass unbalance. Commun Nonlinear Sci 2012; 17: 454–471.

15.

De Castro

Cavalca

Nordmann

. Whirl and whip instabilities in rotor-bearing system considering a nonlinear force model. J Sound Vib 2008; 317: 273–293.

16.

Jing

Meng

Sun

. On the oil-whipping of a rotorbearing system by a continuum model. Appl Math Model 2005; 29: 461–475.

17.

Falkenhagen

Gunter

Schuller

. Stability and transient motion of a vertical three-lobe bearing system. J Eng Ind 1972; 94: 665–676.

Numerical evaluation of multilobe bearings using the spectral method

Abstract

Keywords

Introduction

Model description

Bearing modeling

Numerical methods

FDM

FEM

Pseudo-spectral method

Results

Simulation procedure

Comparison study

Dynamic analysis

Numerical comparison

Dynamic application

Conclusion

Footnotes

Appendix 1

Declaration of conflicting interests

Funding

References