A decoupled approach to enhance the Elrod algorithm

Abstract

The Elrod algorithm is widely used for the study of bearing cavitation problems. However, the potential instability of the algorithm makes it difficult to use in practical applications. A simplified Elrod algorithm based on a decoupled technique was developed. The new algorithm is well-regularized because the tough switch function used for the shear flow term is removed. The most significant advantage of the algorithm is the stability, unlike the traditional algorithm where the accuracy of results irregularly depends on the choice of grids. The numerical experiments, including a journal bearing and a surface-dimpled thrust bearing, show that the simplified Elrod algorithm converges consistently for all cases studied, while the traditional Elrod algorithm fails to converge for cases with a special number of grids. The simplified Elrod algorithm is considered to be a robust alternative to the traditional Elrod algorithm.

Introduction

The Elrod algorithm is a famous cavitation algorithm that can efficiently solve the bearing cavitation problem.^1–8 The algorithm was proposed by Elrod and Adams ⁹ in 1974 and later improved by Elrod ¹⁰ in 1981. Since its inception, the algorithm has attracted wide attention for its mass-conserving approach to the study of bearing cavitation. The most important feature of this algorithm was the proposed switch function, which can automatically capture the free boundaries to reduce numerical difficulties. The algorithm was considered to be the first practical approach to realizing the well-known Jakobsson and Floberg ¹¹ and Olsson ¹² (JFO) theory.

Two widely used versions of the Elrod algorithm are clarified here. ¹³ The first version refers to the original algorithm presented in work ¹⁰ where the shear flow term was innovatively discretized by a switching central and backward difference scheme, with the purpose of accuracy improvement. Another version refers to the modified algorithm developed by Vijayaraghavan and Keith. ¹⁴ In this version, the discrete scheme was improved by incorporating the idea of artificial viscosity. Thus, the compressibility effect was taken into account by preserving variable θ in the full film region. A detailed analysis of the differences between the two versions can be found in work 14.

The Elrod algorithm has the potential instability regardless of the version used. Elrod ¹⁰ initially reported that the alternating-direction implicit (ADI) approach occasionally experienced instability and suggested using a smaller time step to address the issue. Brewe ¹⁵ found that the instability appeared in cases of high journal bearing eccentricity, and the calculation had to be repeated to update the switch function immediately. Qiu and Khonsari ¹⁶ reported that numerical instability appeared in the areas near the grooves of a thrust bearing, which was circumvented by using an approximate pressure of 99999.999 Pa instead of 1 × 10⁵ Pa. Miraskari et al. ¹⁷ found that the instability strongly depended on the choice of the lubricant's bulk modulus.

The instability has caused convergence problems. ¹⁸ Therefore, many improvements to the Elrod algorithm have been proposed.^17–33 Vijayaraghavan and Keith ¹⁹ developed a combined algorithm instead of the ADI method. The algorithm, based on Newton iteration and approximate factorization (AF) techniques, showed superiority in efficiency and robustness. Ausas et al.^20,21 presented a relaxation-type finite-volume-based algorithm that eliminated the switch function by using the if-else-statement, and it was shown to be robust to complex surface-textured bearings. ²⁷ More recently, an optimized discrete scheme of the Ausas algorithm was introduced in work. ²² Giacopini et al. ²³ proposed a novel variable r to represent the void fraction (1 − θ) and then translated the cavitation problem into a linear complementarity problem (LCP). Standard LCP solvers can be used to assist in the solution. The LCP solvers reported so far have been robust. Many improvements have been made to the LCP.^{24–26,28–33} Fesanghary and Khonsari ¹⁸ proposed a modification of the switch function in which the step characteristic of the switch function was mitigated by a proposed gFactor. The initial crash problem of the Elrod algorithm was overcome. Miraskari et al. ¹⁷ developed an alternative solution scheme based on the finite volume method. The algorithm eliminated the scheme's dependence on the choice of lubricant bulk modulus and showed a significant improvement in robustness compared to the AF scheme.

The popular cavitation algorithms have recently been reviewed in work ³² and it can be seen that the balance between complexity and efficiency of these algorithms remains a challenge. In this study, we present a simplified algorithm that enhances the Elrod algorithm. The new algorithm is based on the authors’ previous decoupled technique. ²² Compared to the traditional Elrod algorithm, the new algorithm has the following advantages: (1) The discrete scheme is straightforward. The tough switch function used for the shear flow term is removed. This leads to a well-regularized algorithm with reduced arithmetic cost. (2) The algorithm is more stable. For all studied cases, including a journal and a thrust bearing, the algorithm converges consistently. (3) The algorithm is easy to understand and implement. The decoupled technique simplifies the Elrod algorithm without the need to introduce complex theory.

The outline of this article is as follows. First, the JFO theory is reviewed in part 2. Then, the simplified Elrod algorithm is introduced in part 3. Next, the characteristics of the algorithm are studied in part 4. Finally, the conclusion is summarized in part 5.

JFO theory

The compressible Reynolds equation for a Newtonian lubricant with constant viscosity is expressed as

$${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\rho h^3{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}\right)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}\left(\rho h^3{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }y}}\right)=6\mu U{\displaystyle \frac{\mathrm{\partial }(\rho h)}{\mathrm{\partial }x},}}}$$

(1)

where ρ is the density, h the film thickness, p the pressure, μ the dynamic viscosity, and U the sliding speed. The variable pressure is governed by equation (1). Generally, in the convergent clearance region, the oil is squeezed by the drive of sliding velocity, creating positive pressure. In the divergent clearance region, the oil is forced to expand; however, oil cannot bear large tension, causing oil film rupture. The air/vapor escapes from the oil as cavity bubbles, forming the striations like the finger pattern. ³⁴ Due to the presence of gas, the pressure in the cavitated region is almost constant. Therefore, the pressure flow terms in equation (1) are ignored. Equation (1) reduces to

$${\displaystyle \frac{\mathrm{\partial }(\rho h)}{\mathrm{\partial }x}=0,}$$

(2)

where ρ is the mixed density that represents the mixture of oil and gas. Equation (2) is the cavitation equation.

Equation (1) is a mass continuity equation and can be decomposed into

$${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}(\rho q_x)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}(\rho q_y)=0,}}$$

(\rm 3a)

$$q_x={\displaystyle \frac{Uh}{2}-{\displaystyle \frac{h^3}{12\mu }{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x},}}}$$

(\rm 3b)

$$q_y=-{\displaystyle \frac{h^3}{12\mu }{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }y},}}$$

(\rm 3c)

where q_x and q_y are the volume flows in the x and y directions, respectively. The mass flow is conserved by the x-direction control surface between the liquid and cavitated regions. Thus,

$$({\rho }_c{q}_x{)}_{\mathrm{liquid}}=(\rho q_x{)}_{\mathrm{cavitation}},$$

(4)

where ρ_c is the oil density at the cavitation pressure. Therefore,

$${\rho }_c\left({\displaystyle \frac{Uh}{2}-{\displaystyle \frac{h^3}{12\mu }{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}}}\right)=\rho \left({\displaystyle \frac{Uh}{2}}\right).$$

(5)

Equation(5) is rewritten as

$${\displaystyle \frac{h^3}{12\mu }{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}={\displaystyle \frac{Uh}{2}\left(1-{\displaystyle \frac{\rho }{{\rho }_c}}\right).}}}$$

(6)

Equation (6) is known as the JFO boundary condition. At the location of oil film reformation, mixture density ρ is unknown. Equation (6) indicates the oil film reformation boundary condition. At the location of oil film rupture, mixture density ρ equals oil density ρ_c. Equation (6) reduces to the Reynolds boundary condition:

$${\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}=0.}$$

(7)

Simplified Elrod algorithm

The cavitation problem can be modeled as follows. Find p and θ such that:

$$\{\begin{array}{cc}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\theta h^3{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}\right)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}\left(\theta h^3{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }y}}\right)=6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x}}}}& p>p_c\mathrm{and}\theta \ge 1\\ 6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x}=0}& p=p_c\mathrm{and}\theta <1\end{array},$$

(8)

where θ is the density ratio of ρ/ρ_c. Variable θ has two meanings that represent the density ratio of the film in the full film region and the fractional film content in the cavitated region. The term $1-\theta$ represents the void fraction. Because θ is close to 1 for common lubricants.^10,14,16,18 Equation (8) can be simplified by dropping θ in the full film region, and is integrated into a single equation as

$$\mathrm{Find}p\mathrm{and}\theta {\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(g{h}^3{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}\right)+{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}\left(g{h}^3{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }y}}\right)=6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x},}}}$$

(9)

$$\begin{array}{rl}& \mathrm{subject}\mathrm{to}g=1,\theta =1,\mathrm{for}p>p_c\\ & g=0,p=p_c,\mathrm{for}p\le p_c,\end{array}$$

(10)

where g is the switch function proposed by Elrod and Adams. ⁹

Equation (9) is known to have mixed behavior: elliptic for the full film region and hyperbolic (first-order differential equation) for the cavitated region. The shear flow term is needed to design a special discrete scheme that can be switched with the change of regions. The frequently used discrete schemes are as follows:

The Elrod version^10,14,15:

$$6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x}=-{\displaystyle \frac{6\mu U}{\mathrm{\Delta }x}\left[\begin{array}{c}(1-g_{i-1,j})h_{i-1,j}{\theta }_{i-1,j}-(1-g_{i,j})h_{i,j}{\theta }_{i,j}+{\displaystyle \frac{(2-g_{i,j})g_{i-1,j}}{2}{h}_{i-1,j}}\\ +{\displaystyle \frac{(g_{i-1,j}-2+g_{i+1,j})g_{i,j}}{2}{h}_{i,j}-{\displaystyle \frac{g_{i+1,j}{g}_{i,j}}{2}{h}_{i+1,j}}}\end{array}\right].}}$$

(11)

The Vijayaraghavan and Keith version ¹⁴ :

$$6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x}={\displaystyle \frac{6\mu U}{2\mathrm{\Delta }x}\left[\begin{array}{c}\left({\displaystyle \frac{g_{i,j}+g_{i+1,j}}{2}}\right)h_{i+1,j}{\theta }_{i+1,j}+\left(2-{\displaystyle \frac{g_{i-1,j}+2g_{i,j}+g_{i+1,j}}{2}}\right)h_{i,j}{\theta }_{i,j}\\ +\left(-2+{\displaystyle \frac{g_{i-1,j}+g_{i,j}}{2}}\right)h_{i-1,j}{\theta }_{i-1,j}\end{array}\right].}}$$

(12)

The discrete schemes of Equations (11) and (12) are complex, lengthy and dependent on the switch function. To make matters worse, they are believed to be the main source of potential instability. The numerical experiments performed later show that the accuracy of the results strangely depends on the choice of the number of grids. The results will be inaccurate even if the grid density is high.

Inspired by the discrete scheme, ²² the shear flow term can naturally be designed to satisfy the required special discrete scheme. The only thing needs to do is to decouple θ and h. Since θ is close to 1 in the full film region, variable θ can be written as $\theta =\theta -1+1$ . The shear flow term is then rearranged as

$$6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x}=6\mu U{\displaystyle \frac{\mathrm{\partial }((\theta -1+1)h)}{\mathrm{\partial }x}.}}$$

(13)

Thus,

$$6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x}=6\mu U{\displaystyle \frac{\mathrm{\partial }((\theta -1)h)}{\mathrm{\partial }x}+6\mu U{\displaystyle \frac{\mathrm{\partial }h}{\mathrm{\partial }x}.}}}$$

(14)

The new discrete scheme is designed as

$$6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x}=6\mu U{\displaystyle \frac{{(\theta -1)}_{i,j}{h}_{i,j}-{(\theta -1)}_{i-1,j}{h}_{i-1,j}}{\mathrm{\Delta }x}+6\mu U{\displaystyle \frac{h_{i+1,j}-h_{i-1,j}}{2\mathrm{\Delta }x}.}}}$$

(15)

Equation (15) achieves the same function as Equations (11) and (12), but switch function g is removed. If all the grid points are within the full film region, $\theta =1$ , the shear flow term produces the desired central difference scheme for h. If all the grid points are within the cavitated region, the shear flow term produces the desired upwind scheme for θ.

Table 1 compares the discrete schemes of the shear flow term for the three versions of the Elrod algorithm. The subtle differences between them can affect the accuracy and stability of the algorithm accordingly. In general bearing analysis, the accuracy differences can be safely ignored, as demonstrated in the numerical validation section in 4.1. However, for high pressure problems, such as the elastohydrodynamic lubrication problems,^35,36 the accuracy difference is worthy of study. Here, the concern is focused on the stability difference. The new discrete scheme uses θ to achieve the same function as g, resulting in two key differences: (1) g is strongly stepwise, while θ varies relatively continuously. (2) g is updated after a loop, while θ is updated within a loop. Variable θ changes more smoothly than g, which benefits iterative computations, especially stability.

Table 1.

Comparison of discrete schemes for the shear flow term.

Algorithm	Full film region	Cavitated region	Use of switch function
Original Elrod algorithm 10	$6\mu U\frac{h_{i+1,j}-h_{i-1,j}}{2\mathrm{\Delta }x}$	$6\mu U{\displaystyle \frac{{\theta }_{i,j}{h}_{i,j}-{\theta }_{i-1,j}{h}_{i-1,j}}{\mathrm{\Delta }x}}$	Yes
Modified Elrod algorithm 14	$6\mu U\frac{{\theta }_{i+1,j}{h}_{i+1,j}-{\theta }_{i-1,j}{h}_{i-1,j}}{2\mathrm{\Delta }x}$	$6\mu U{\displaystyle \frac{{\theta }_{i,j}{h}_{i,j}-{\theta }_{i-1,j}{h}_{i-1,j}}{\mathrm{\Delta }x}}$	Yes
Simplified Elrod algorithm	$6\mu U\frac{h_{i+1,j}-h_{i-1,j}}{2\mathrm{\Delta }x}$	$\begin{array}{rl}& 6\mu U{\displaystyle \frac{{(\theta -1)}_{i,j}{h}_{i,j}-{(\theta -1)}_{i-1,j}{h}_{i-1,j}}{\mathrm{\Delta }x}}\\ & +6\mu U{\displaystyle \frac{h_{i+1,j}-h_{i-1,j}}{2\mathrm{\Delta }x}}\end{array}$	No

Observing Equation (15), the term $(\theta -1)$ appears. This term happens to appear in the approximate pressure-density relation ¹⁰ :

$$p=p_c+g\beta (\theta -1).$$

(16)

Equation (16) is a key point in making a universal partial differential equation, because the discrete scheme of the pressure flow term is determined by Equation (16). The next intention is to simplify $(\theta -1)$ by defining a new variable $\phi =\theta -1$ .

The shear flow term is thus discretized as

$$6\mu U{\displaystyle \frac{\mathrm{\partial }(\theta h)}{\mathrm{\partial }x}=6\mu U{\displaystyle \frac{{\phi }_{i,j}{h}_{i,j}-{\phi }_{i-1,j}{h}_{i-1,j}}{\mathrm{\Delta }x}+6\mu U{\displaystyle \frac{h_{i+1,j}-h_{i-1,j}}{2\mathrm{\Delta }x}.}}}$$

(17)

The discrete scheme of Equation (17) is simpler than that of the traditional versions in Equations (11) and (12) and is independent of the switch function.

The pressure flow terms are approximated by substituting p with Equation (16):

$${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(h^3{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }x}}\right)={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\beta h^3{\displaystyle \frac{\mathrm{\partial }(g(\theta -1))}{\mathrm{\partial }x}}\right)={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(\beta h^3{\displaystyle \frac{\mathrm{\partial }(g\phi )}{\mathrm{\partial }x}}\right),}}}$$

(18)

$${\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}\left(h^3{\displaystyle \frac{\mathrm{\partial }p}{\mathrm{\partial }y}}\right)={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}\left(\beta h^3{\displaystyle \frac{\mathrm{\partial }(g(\theta -1))}{\mathrm{\partial }y}}\right)={\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}\left(\beta h^3{\displaystyle \frac{\mathrm{\partial }(g\phi )}{\mathrm{\partial }y}}\right).}}}$$

(19)

The discrete scheme is

$$\begin{array}{rl}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }x}}& \left(\beta h^3{\displaystyle \frac{\mathrm{\partial }(g\phi )}{\mathrm{\partial }x}}\right)={\displaystyle \frac{\beta }{{(\mathrm{\Delta }x)}^2}[h_{i+1/2,j}^3{g}_{i+1,j}{\phi }_{i+1,j}}\\ & -(h_{i+1/2,j}^3+h_{i-1/2,j}^3)g_{i,j}{\phi }_{i,j}+h_{i-1/2,j}^3{g}_{i-1,j}{\phi }_{i-1,j}],\end{array}$$

(20)

$$\begin{array}{r}{\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }y}\left(\beta h^3{\displaystyle \frac{\mathrm{\partial }(g\phi )}{\mathrm{\partial }y}}\right)={\displaystyle \frac{\beta }{{(\mathrm{\Delta }y)}^2}[h_{i,j+1/2}^3{g}_{i,j+1}{\phi }_{i,j+1}}}\\ -(h_{i,j+1/2}^3+h_{i,j-1/2}^3)g_{i,j}{\phi }_{i,j}+h_{i,j-1/2}^3{g}_{i,j-1}{\phi }_{i,j-1}].\end{array}$$

(21)

The switch function g has been moved to the derivative term of φ.

Equations (17), (20), and (21) are the core of the simplified Elrod algorithm. Variable φ represents the film increase in the full film region and the film decrease (negative void) in the cavitated region. Explicit and implicit methods are both applicable to the simplified Elrod algorithm. The basic iterative scheme for solving φ is

$${\phi }_{i,j}={\displaystyle \frac{\begin{array}{c}{\displaystyle \frac{\beta }{{(\mathrm{\Delta }x)}^2}(h_{i+1/2,j}^3{g}_{i+1,j}{\phi }_{i+1,j}+h_{i-1/2,j}^3{g}_{i-1,j}{\phi }_{i-1,j})+{\displaystyle \frac{\beta }{{(\mathrm{\Delta }y)}^2}(h_{i,j+1/2}^3{g}_{i,j+1}{\phi }_{i,j+1}+h_{i,j-1/2}^3{g}_{i,j-1}{\phi }_{i,j-1})}}\\ +{\displaystyle \frac{6\mu U}{\mathrm{\Delta }x}\left({\phi }_{i-1,j}{h}_{i-1,j}-{\displaystyle \frac{h_{i+1,j}-h_{i-1,j}}{2}}\right)}\end{array}}{{\displaystyle \frac{\beta }{{(\mathrm{\Delta }x)}^2}(h_{i+1/2,j}^3+h_{i-1/2,j}^3)g_{i,j}+{\displaystyle \frac{\beta }{{(\mathrm{\Delta }y)}^2}(h_{i,j+1/2}^3+h_{i,j-1/2}^3)g_{i,j}+6\mu U{\displaystyle \frac{h_{i,j}}{\mathrm{\Delta }x}}}}}.}$$

(22)

The arithmetic operation of equation (22) is less than that of the original Elrod algorithms presented in work. ¹⁸ The interested reader can reproduce the simplified Elrod algorithm with equation (22).

The solution scheme is shown in Table 2. The successive over-relaxation method can be used to accelerate convergence. The modified switch function developed by Fesanghary and Khonsari ¹⁸ is used to avoid the initial crash problem. The convergence of the solution is determined by monitoring the change in the 1-norm of the vector φ. The convergence criterion is defined as $\Vert {\phi }^{\mathrm{new}}-{\phi }^{\mathrm{old}}{\Vert }_1/\Vert {\phi }^{\mathrm{old}}{\Vert }_1<\xi$ , where ξ is the user-defined tolerance. The convergence criterion can be adjusted according to requirements.

Table 2.

Solution scheme.

while residual > tolerance   ${\phi }^{\mathrm{old}}=\phi$  for $i=1:\mathrm{end}$   for $j=1:\mathrm{end}$     ${\phi }^{new}\leftarrow$ Eq. (22)     ${\phi }_{i,j}=({\phi }^{new}-{\phi }_{i,j})\omega +{\phi }_{i,j}$   end  end  for $i=1:\mathrm{end}$   for $j=1:\mathrm{end}$    if ${\phi }_{i,j}\ge 0$ ; $g_{i,j}=g_{i,j}/gF$ ;else; $g_{i,j}=g_{i,j}\cdot gF$ ;end    if $g_{i,j}>1$ ; $g_{i,j}=1$ ;end   end  end  residual =  $\Vert (\phi -{\phi }^{\mathrm{old}}){\Vert }_1/\Vert {\phi }^{\mathrm{old}}{\Vert }_1$ end

Algorithm analysis

Numerical validation

Three different types of bearings were used to analyze the present algorithm. The first is the one-dimensional hydrodynamic bearing with a convergent–divergent sinusoidal profile. The pressure distribution of the bearing was validated by the different cavitation algorithms performed by Giacopini et al. ²³ and Miraskari et al. ¹⁷ The second is the submerged journal bearing studied experimentally by Jakobosson and Floberg. ¹¹ The pressure distribution measured was consistent to the numerical results obtained by Brewe ¹⁵ and Cupillard et al. ³⁷ The third is the dimple-enhanced seal-like thrust bearing studied by Qiu and Khonsari. ¹⁶ The pressure distribution based on the Elrod algorithm was provided by Qiu and Khonsari. ¹⁶ These bearing parameters are listed in Table 3. The periodic boundary condition is applied to cases 2 and 3.

Table 3.

Parameters for the studied bearings.

Case	Dimensional	Type	Textured	Periodic	Geometric parameter	Operation condition
Case 1	1d	Plane	No	No	hmax = 0.025 mm, hmin = 0.015 mm, a = 125 mm	μ = 0.015 Pa·s, U = 4 m/s, Pa = 1 MPa, Pc = 0 MPa, β = 10 GPa
Case 2	2d	Journal	No	Yes	D = 100 mm, L = 133 mm, ε = 0.60	n = 459 r/min, pa = 0.1 MPa, pc = 28 kPa, μ = 0.0127 Pa·s, β = 1.72 GPa
Case 3	2d	Thrust	Yes	Yes	r0 = 750 μm, hg = 10 μm, h0 = 4 μm, L = 3 mm, h0 = 4 μm	μ = 0.0035 Pa·s, U = 1.45 m/s, pa = 0.1 MPa, pc = 0.09 MPa, β = 1.62 GPa

Figure 1 shows the pressure distributions for the simplified Elrod algorithm and other methods. It demonstrates good agreement between the resulting pressures obtained by the simplified Elrod algorithm and those obtained by the other methods for all three cases, indicating the reliability of the simplified Elrod algorithm. In addition, to facilitate a convenient evaluation of the plane bearing, the one-dimensional bearing is considered as a two-dimensional bearing by setting the y-direction length of the bearing to ten times the x-direction length (i.e., L = 10a).

Figure 1.

Pressure distribution for simplified Elrod algorithm and other methods. (a) Plane bearing with gird 64 × 10, (b) journal bearing gird of 64 × 64, and (c) thrust bearing with gird 64 × 64.

Results and discussion

1. Algorithm instability

The instability of the Elrod algorithm has two meanings. The first meaning refers to the program crashing during iteration, where the asymptotic convergent behavior is broken. The crash typically occurs at the start of the iteration, resulting from the step change of the switch function, where the difference in the variable field between two neighboring iterations exceeds the tolerance of the system. However, the crash has already been overcome by using the modified switch function with the help of a constant gFactor. The choice of gFactor is wide and does not affect the accuracy of the results. The value of 0.8 for gFactor suggested by Fesanghary and Khonsari ¹⁸ makes all the calculations proceed normally and is adopted here. The initial crash problem is not the focus of this study.

The second meaning refers to the nonconvergent results caused by using a special number of grids. Some special grids, which are unknown a priori, can lead to nonconvergence, where the residual cannot fall sufficiently. In the study, the residual larger than 1 × 10⁻¹³ is considered nonconvergent. Figure 2 shows an example of the residual variations with iterative progress for the journal bearing case. The residuals were analyzed by the three versions of the Elrod algorithm mentioned previous. For the original Elrod algorithm, the residual for the grid with 64 × 63 nodes falls below 1 × 10⁻¹³, while the residual for the grid with 64 × 64 nodes stops falling at 1.9 × 10⁻⁶. For the modified Elrod algorithm by Vijayaraghavan and Keith, the residual for the grid with 64 × 64 nodes falls below 1 × 10⁻¹³, while the residual for the grid with 64 × 68 nodes stops falling at 4.7 × 10⁻⁶. The grids 64 × 64 and 64 × 68 are considered as the special grids for the two traditional Elrod algorithms. When the special grids are implemented in the simplified Elrod algorithm; however, the expected nonconvergence does not occur. Therefore, it would be interesting to investigate the stability of the simplified Elrod algorithm.

Figure 2.

Residual variations with iterative progress for the journal bearing case. (a) Original Elrod algorithm; (b) modified Elrod algorithm; and (c) simplified Elrod algorithm.

There are two other interesting phenomena that can be seen from Figure 2. The first is that the convergence speed of the simplified algorithm is slower than that of the traditional algorithms. The speed difference, however, is very small and still within the same order of level. As these algorithms are explicit relaxation methods and only suitable for small numbers of grid cases, the calculated time difference cannot be reflected effectively. The second is that the residual oscillation of the simplified algorithm is severer than that of the traditional algorithms. This is caused by the discrete scheme difference. The oscillation does not seem to have a negative impact on the simplified algorithm.

2. Hazards of nonconvergence

This part is to analyze the reason why the residuals stop falling and what the hazard of nonconvergence is. Figure 3 shows the difference in θ between two neighboring iterations. High residuals can be seen at some nodes close to the film reformation boundary, as shown by the white lines. ³⁸ Although these residuals are limited to only the local 2 ∼ 6 nodes, the maximum value is as high as the level of 1 × 10⁻³. The accuracy of the resulting pressure must be reduced.

Figure 3.

Difference in θ between two neighboring iterations. (a) Grid 64 × 64 and (b) grid 64 × 68.

Local slight difference has little effect on the overall accuracy of the results. However, in some cases, unreliable pressure may appear. Figure 4 shows the results for the journal bearing case under special grids. The high residual near the film rupture boundary results in the singularity of pressure and density ratio. Notice that the correct maximum pressure only reaches up to 0.38 MPa, as shown in Figure 1(b). The predicted load-carrying capacity should become unreliable. The Elrod algorithm fails on the special grids. Nonconvergent grids should be avoided, even if they are unknown a priori. The practical approach to ensure accuracy is to try different grids until the residuals converge, but this is time consuming.

3. Stability analysis

Figure 4.

Results of the journal bearing case under special grids. (a) Grid 60 × 68, original Elrod algorithm. (b) Grid 65 × 63, modified Elrod algorithm.

The accuracy of the Elrod algorithm has been shown to be related to the number of grids. The primary attention is on the nonconvergence probability for the three versions of the Elrod algorithm. Three grid densities, representing coarse, medium, and fine grids, were studied, as listed in Table 4. Each case has 100 grids to calculate. Notice that relaxation factor ω is kept as 1 to exclude the instability induced by ω. Table 5 shows the statistics of nonconvergent cases for coarse, medium, and fine grids. The detailed nonconvergent distributions of x-nodes and y-nodes are given in the appendix. From Table 5, it can be concluded that: (1) Both the original Elrod algorithm and the modified Elrod algorithm by Vijayaraghavan and Keith have potential instability. The nonconvergence probability is very high for the journal bearing case. The probability can reach up to 82% for the original Elrod algorithm. The thrust bearing cases can almost converge, with only three nonconvergent cases appearing in the modified Elrod algorithm by Vijayaraghavan and Keith. (2) The nonconvergence probability cannot be reduced by increasing the grid density. Conversely, the probability increases as the grid density increases. It makes the Elrod algorithm inapplicable to complex bearings. (3) The simplified Elrod algorithm converges consistently for all cases studied. The switch function used for the shear flow term is considered to be one of the main reasons for nonconvergence. The removal of the switch function for the shear flow term makes the Elrod algorithm more stable.

Table 4.

Coarse, medium, and fine grids.

Grid density	x-nodes × y-nodes
Coarse	(30 ∼ 39) × (30 ∼ 39)
Medium	(60 ∼ 69) × (60 ∼ 69)
Fine	(100 ∼ 109) × (100 ∼ 109)

Table 5.

Statistics of nonconvergences for coarse, medium, and fine grids.

Algorithm	Number of nonconvergent grids
	Journal bearing	Thrust bearing
Original Elrod algorithm	71,65,82	0,0,0
Modified Elrod algorithm	44,43,56	0,1,2
Simplified Elrod algorithm	0,0,0	0,0,0

Conclusion

The study clarified the versions of the Elrod algorithm and highlighted the drawback of the potential instability. A simplified Elrod algorithm based on the decoupled technique was developed. The discrete scheme, arithmetic cost, and regularity of the algorithm are superior to those of the traditional Elrod algorithms. The most significant advantage of the algorithm is the stability, unlike the traditional algorithm where the accuracy of results depends irregularly on the choice of grids and the convergent grids cannot be known a priori. The new algorithm is considered to be a robust alternative to the traditional Elrod algorithm.