Thermo-Hydrodynamic Analysis of a Plain Journal Bearing on the Basis of a New Mass Conserving Cavitation Algorithm

Accurate prediction of cavitation is an important feature in hydrodynamic bearing modeling. Especially for thermo-hydrodynamic modeling, it is crucial to use a mass-conservative cavitation algorithm. This paper introduces a new mass-conserving Reynolds cavitation algorithm, which provides fast convergence and easy implementation in finite element models. The proposed algorithm is based on a variable transformation for both the pressure and mass fraction, which is presented in the form of a complementary condition. Stabilization in the streamline and crosswind direction is provided by artificial diffusion. The model is completed by including a simple and efficient thermal model and is validated using the numerical values of a reference plain journal bearing experiment under steady-state conditions. In addition, a transient analysis is performed of a journal bearing subjected to a harmonic load. It is shown that the proposed cavitation algorithm results are in good agreement with the reference measurement results. Moreover, the algorithm proves to be stable and requires only a small number of iterations to convergence in the Reynolds-based finite element model.


Introduction
Hydrodynamic bearings are machine elements used for separating moving components that have a relative velocity. In these bearings, the relative sliding velocity of the surfaces that border the fluid film causes a shear flow in the fluid. When this flow enters a converging volume, the pressure increases, which explains the carrying capacity of the bearing. Accurate predictions of this carrying capacity are crucial for an accurate design and, thus, for preventing any unwanted contact that can result in noise, wear, excessive friction or even failure of the bearing [1].
The commonly-used mathematical description of the pressure distribution in a hydrodynamic bearing is the classical Reynolds equation. Solutions of the Reynolds equation can contain negative pressures for flow through diverging sections or for the normal motion of walls moving apart. Fluids, however, can only hold a very limited amount of tensile stress before film rupture occurs [2]. This cavitation thus needs to be included in the bearing modeling. Numerous methods have been proposed for this purpose, as summarized in the review papers of Dowson and Braun [2,3].
In many of the early approaches, such as the Sommerfeld, the Gümbel "half-Sommerfeld" and the Swift and Stieber Reynolds approach, conservation of mass throughout the fluid and cavitation domain is not considered [3]. This shortcoming might sometimes only slightly affect the bearing reaction forces; conservation of mass is, however, an important condition for models that should include thermal modeling [4] or for models that treat surface texture, especially if accurate friction losses are to be calculated [5,6].
A popular cavitation algorithm nowadays was proposed by Elrod [7], who developed a method in which a switch function is used, automatically satisfying the boundary conditions at the cavitation interfaces, and which is based on the mass-conservative boundary conditions formulated by Floberg, Jakobson and Olsson [8]. Brewe [9,10] found excellent agreement with Jakobson's experimental data for dynamically loaded bearings using Elrod's algorithm.
Cavitation algorithms based on a finite element formulation are compared in Bayada et al. [11]. Many of these algorithms are presented in the form of a complimentary condition, i.e., in the cavitated regions, the mass fraction distribution is unknown, and the pressure distribution is known, whereas in the full film region, the pressure distribution is unknown, and the mass fraction distribution is known.
In the proposed algorithm in this paper, the same logic has been implemented, based on a variable transformation where both the mass fraction f and the pressure p are replaced by known functions of a new variable ξ. Thus, the modified Reynolds equation can be solved for this variable, and the pressure and mass fraction distribution can be obtained from the inverse transformation functions. The proposed mass-conserving cavitation algorithm is developed and used in a general purpose finite element code [12], which allows for the formulation of problems as a set of strongly coupled non-linear partial differential equations (PDEs). Moreover, it can be used in any finite element code that allows the arbitrary definition of the PDE coefficients and the solution of non-linear systems of equations using a simple fixed-point iteration.
In this paper, the cavitation algorithm is first applied on a 1D test problem, the infinitely long journal bearing, in order to demonstrate the streamline stabilization technique. Subsequently, the 1D stabilized model is expanded to a 2D finite journal bearing model in order to demonstrate the crosswind stabilization technique from which the complete stabilized 2D Reynolds equation is obtained. This model is then expanded to include a classic adiabatic thermal model, i.e., assuming adiabatic boundary conditions at the fluid-solid interfaces and assuming an effective temperature field across the film thickness (see, for instance, Pinkus [13]). These simplifications reduce the 3D heat transport model to a simple and efficient 2D model for the thin oil film. Subsequently, the 2D THD model is validated by a numerical case study of a finite length journal bearing under steady load conditions [14]. Finally, for the sake of completeness, a time transient analysis is performed using a dynamic load on the journal, after which the results are discussed.

Development of a New Mass-Conserving Reynolds Cavitation Algorithm
The distribution of pressure p and mass fraction f in a 2D lubricating film with film thickness h can be described using the Reynolds equation: where η, ρ and U are the viscosity, density and surface velocity, respectively. xand y-directions are in the direction of and perpendicular to the surface sliding velocity, respectively. It is assumed that only the shaft has velocity U in the x-direction; all other surfaces are stationary. Moreover, all surfaces are assumed to be rigid. As schematically depicted in Figure 1a, the attitude angle X ε can be calculated as follows: where ε is the eccentricity ratio (ε = e h 0 ).x andȳ are the global coordinates with the origin in the center of the bearing. Assuming perfect alignment of the journal and taking into account the attitude angle X ε , the expression for the film thickness may be written as: where h 0 is the nominal radial clearance. For the sake of numerical robustness, the problem is scaled using the following dimensionless variables: which results in the following dimensionless Reynolds equation: where the normalized film thickness is given by: Note that the circumferential and axial flows, q X and q Y , respectively, are represented by the terms between brackets ∂ ∂X (q X ) and ∂ ∂Y (q Y ) in Equation 5 and can thus simply be computed as follows: In Figure 1a,b, the kinematics of the plain journal bearing and the computational fluid domain with the imposed boundary conditions are presented, respectively. The supply groove is positioned at X = −π, X = π, over the full width of the bearing. The specific boundary conditions (BCs) for the 1D and 2D case will be defined in each section. In Figure 1, the X-and Y -directions are in the direction of and perpendicular to the surface sliding velocity, respectively. For most practical applications where the hydrodynamic pressure will increase to values much higher than the cavitation pressure, it may be assumed that the cavitation pressure is equal to zero. With this simplification, the lubricating film can be divided into two parts: • the full film region, where P is unknown (but larger than zero) and f is known (f = 1), and • the cavitated region, where P is known (P = 0) and f is unknown (but between zero and one). This logic has been captured by the proposed cavitation algorithm, based on a variable transformation where both the mass fraction f and the pressure p are replaced by known functions of a new variable ξ: where (ξ < 0) and (ξ ≥ 0) are Boolean expressions, equal to one if true and zero if not true. This transformation resembles Elrod's algorithm [7], although the algorithm proposed here has some different properties. Firstly, the Boolean functions are present both in the pressure and in the mass fraction functions, instead of just in the pressure function. Furthermore, the pressure and mass fraction functions are independent of each other. Thus, both functions can be optimally scaled for a stable numerical solution, which is an advantage compared to Elrod's transformation, i.e., the transformation constant c f can be equal to any desired positive value. For Elrod's transformation to be numerically stable, the bulk modulus has to be 10-to 100-times smaller than the value considered by Elrod, which at least reduces the physical meaning of the transformation [15,16]. Moreover, the accuracy of Elrod's transformation is not always optimal, as it will result in mass fractions slightly above one in full film regions.

1D Stationary Situation: Streamline Numerical Stabilization
The characteristics of the proposed transformation will first be tested on a simple 1D problem: the statically-loaded, infinitely long plain journal bearing. In the remainder of this section, it is assumed that the supply groove is infinitely small in the circumferential direction. Furthermore, we assume incompressible (ρ = 1) lubricant properties. Considering this simplification and substituting Equation (8a) and (8b) in Equation (5), the following 1D Reynolds equation is obtained: In the 1D-domain, the dimensionless pressure is set equal to one at the supply (X = −π, X = π) (BC1: P = 1). It is worth mentioning that in this analysis, it is assumed that the derivative of the Boolean expression is zero and not equal to the Dirac-delta function. It is difficult to solve Equation (9) using the standard Galerkin finite element approach because in the cavitated region (ξ < 0), the equation is purely convection dominated and, thus, will exhibit oscillations when solved using a standard Galerkin FEM. This can be observed in Figure 2a, where the dependent variable ξ and the dimensionless flow (q X ) have been plotted against the dimensionless length. Plotting ξ provides direct insight into both the pressure and mass fraction distribution.
In the literature [17][18][19], one may find many approaches to numerically stabilize the FEM solution of the convection-dominated part in convection-diffusion equations. Some of these methods are the streamline upwind Petrov-Galerkin method (SUPG), the SU method and the artificial diffusion method (AD). The SUPG method adds an additional term to the test function to suppress the unwanted oscillations, whereas the SU method only adds the additional term to the convection-dominated part of the equation [19]. In the AD method, the diffusion part, of the convection-diffusion equation, is slightly modified with an additional artificial diffusion term in the convection-dominated part of the film.
Other techniques in the literature that deal with the convection-dominated aspect include semi-Lagrangian (or characteristic) methods [20,21]. The semi-Lagrangian technique has the advantage that it does not require the tuning parameters associated with the AD-method. However, this method requires substantial numerical effort and is not easily implemented in any general purpose software package, as it requires fundamental rewriting of the finite element code. Advanced general purpose FEM software is commonly used in the industry, as the lubrication problems usually encountered involve a multidisciplinary field, including thermal, structural, thermo-elastic and hydrodynamic simulations.
Hence, the chosen stabilization technique here is the artificial diffusion method (AD), as it can easily be implemented and only requires a slight modification to the coefficients in the general form PDE. In this method, the diffusion part in the convection-diffusion equation is slightly modified with an additional artificial diffusion term k AD in the convection-dominated part of the film: Note that the solution obtained using this method should be examined carefully as a slightly modified PDE is solved instead of the original one. Moreover, k AD should not be so large as to introduce excessive damping, but should be minimal just to suppress the generated oscillations. Hence, an optimum value for k AD needs to be chosen. The artificial diffusion coefficient k AD is chosen as follows: where h e is a typical element size in the streamline direction and α 0 is a constant that depends on the order of the elements that are used with α 0 = 0.5 for linear elements and α 0 = 0.25 for quadratic elements [22]. In the present analysis, linear elements are used. u and k are respectively equivalent to the velocity and diffusion coefficient in a typical convection-diffusion equation and can be calculated as follows: α AD is an upwind function and can (after making an asymptotic approximation [22]) be computed as follows: where P e AD is the cell Peclet number (P e AD = uhe 2k ). Due to the Boolean expressions, P e AD becomes zero in the full film region and infinite in the cavitated region. Consequently, the upwind function α AD reduces to one in the cavitated region and zero in the full film region. Substituting Equation (11) into Equation (10) results in the following Reynolds equation: This equation is strongly non-linear due to the included Boolean expressions. After discretization using finite elements, the system of non-linear equations is solved using a simple fixed-point iteration where the Boolean expressions are evaluated by means of the solution of the previous iteration, which then results in a linear system of equations.
The dependent variable ξ and flow q X are plotted in Figure 2b after solving Equation (14). The results are compared with another streamline stabilization technique (SUPG) in Figure 2b from which it can be concluded that the results obtained using AD and SUPG are in close agreement. Some discontinuity of flow at the reformation boundary remains, however, which should be continuous from the physical perspective. Minimization of this discontinuity can be obtained by choosing an optimum value for the transformation constant c f . An attempt to resolve this challenge is by equating the flows at this boundary and working out using numerical differences; or: where q X,cav and q X,full stand for the flow in the cavitated and full film region, respectively. Substituting the proposed transformation (Equation (8a) and (8b)) into Equation (15b) and working out in terms of numerical differences gives: where H r and ξ r are the film thickness and solution at this interface point. The partial derivative is approximated using the numerical difference between the values of point r and downstream neighboring point r + 1. In the cavitated region, the pressure is less than zero (ξ r < 0), whereas in the full film region, the pressure is equal to or larger than zero (ξ r+1 ≥ 0). According to this, a critical transformation constant (c * f ) may be introduced, which complies with the following expression: For values c f close to c * f , the discontinuity at the reformation boundary was minimal, and faster convergence was achieved; see Figure 3a. As can be observed, c * f depends on the local height at the reformation interface, which is initially unknown. However, after performing some numerical experiments, it was observed that in most cases, H r is found to be very close to the maximum film height of the oil film. Hence, the maximum film height was used in the calculation of c f .
With the dependence of c * f on H r , Equation (14) can now be simplified by suggesting Note the absence of any Boolean functions in the pressure term. This optimized transformation results in fast convergence, as it is slightly less non-linear and, thus, more stable. In Figure 3b, the mass fraction, film thickness, flow and pressure are plotted after solving Equation (18). It is clear that the flow discontinuity is minimal now. Note that due to this chosen value of c f = H 2 3ηhe , which is dependent on the local element size h e , the flow discontinuity is effectively minimized by locally refining the mesh size near the reformation boundary. This 1D cavitation algorithm has been verified for the full range of shaft eccentricity values and for a large range of mesh densities, for low to high numbers of degrees of freedom N dof . Figure 4a gives a representation of the required number of iterations N iter. and N dof , in order to obtain relative errors of 10 −4 − 10 −5 , as a function of ε. As can be observed, the maximum number of iterations required is seven, which contributes as evidence for fast convergence. Typically, the computation time is within a few seconds on a quad core desktop computer.
It should be noted that if one only considers the pressure and/or the mass fraction profiles, it may seem that the solution is smooth enough. The streamline flow q X , however, may exhibit oscillations due to the existence of steep pressure gradients, in particular at higher values for ε. Of course, this means that the mesh size has to be chosen smaller. This problem is clearly exposed in Figure 4b. Figure 4c,d shows that the results obtained using our AD-stabilized cavitation algorithm are identical to "exact" results calculated by SUPG-stabilization, even for relatively high values of ε. This comparison shows the numerical robustness of the cavitation algorithm for high ε values. (c) Comparison between pressure distributions obtained using the AD method with the SUPG method for relatively high eccentricity ratios. The results strictly follow each other. (d) Comparison between mass fraction distributions obtained using the AD method with the SUPG method. Even for high ε, minimal discontinuity is observed near the reformation boundary. Furthermore, excellent agreement exists between the results obtained using both methods.
To further validate our algorithm, the cavitation algorithm was used in several case studies for which there are reliable benchmark solutions in the literature. Giacopini et al. [23] developed a mass-conserving cavitation algorithm based on the concept of complementarity and successfully verified their model using several case studies with (semi)-analytical and numerical formulations. In the present paper, the proposed cavitation model will be adapted to three different case studies in order to test its numerical properties and finally verify the results with those of Giacopini et al. [23].
The first situation represents a simple journal bearing with an infinite width (see Figure 5a). The dimensions and relevant operating conditions are given in Figure 5. As can be observed, the film rupture and reformation boundary are accurately calculated. Excellent agreement is found between the predicted and reference pressure profiles.
The second situation is somewhat different, as the film thickness profile has a divergent-convergent shape (see Figure 5b). As can be extracted from Figure 5b, the film cavitates in the diverging region. Again, there is excellent accordance among both the predicted and reference results. Note that the diverging-converging shape of the film thickness can also be considered to be an idealization for situations involving surface texture. Hence, the last situation is that of a non-convergent/parallel bearing with a single microtexture pocket near the inlet (see Figure 6a). The reference dimensions and operating conditions are given in Figure 6. For this case, Olver et al. [24] derived an analytical solution, which pointed out the existence of a region in which the film cavitates. Comparisons of the predicted pressure profile with those of the reference authors are provided in Figure 6b in which the agreement is excellent. Figure 6c compares the predicted fractional film content distribution with the reference results in which fair accordance is obtained. Note that the fractional film content in the work of Giacopini et al. [23] is modeled using compressibility and is thus modeled as the non-dimensional density, which is why we have used a different notation for the fractional film content for the results of the reference authors. However, in the cavitated region, the pressure is constant, and thus, ρ/ρ 0 stands for the fractional film thickness. Giacopini et al. [23]. There is fair accordance between the results.

2D Stationary Situation: Crosswind Numerical Stabilization
In the previous section, the 1D cavitation algorithm was shown to be stable and robust and, moreover, produces results identical to the results found in the literature. Therefore, the algorithm will now be extended to a 2D cavitation algorithm. The extension of the stabilized 1D Reynolds equation (Equation (18)) to a 2D finite length journal bearing reads: For this case, the supply is over the full bearing length (BC1: P = 1), and the pressure on the edges of the bearing is set equal to zero (BC2: P = 0). Furthermore, a relatively coarse grid is used to show the numerical stability of the stabilization techniques used.
When Equation (19) is solved without further stabilization, the distribution of the mass fraction is not smooth, and spikes occur at the edges of the bearing, as can be seen in Figure 7a. This means that only applying streamline stabilization is not sufficient. Numerical stabilization is required in the crosswind direction (Y -direction), as well as in order to smooth the solution of the flow in the cavitated region. In Hughes et al. [19], Codina [22] and Tezduyar and Park [25], for instance, methods have been studied to choose the crosswind artificial diffusion coefficient. In this study, the diffusion term in the Y -component of the Reynolds equation is modified by adding an artificial diffusion term k AD,Y : In this study, the crosswind diffusion coefficient is chosen similar to k AD , but scaled relative to the angle between the flow solution gradient and flow direction. The crosswind artificial diffusion coefficient k AD,Y is chosen as follows: where δ is a very small number (order of magnitude 1e-10) to prevent numerical problems in areas where the gradients ∂ξ ∂X and ∂ξ ∂Y are zero. Substituting this value for k AD,Y in Equation (20) results in the stationary cavitated-modified Reynolds equation for a 2D lubricating film: A convergence analysis was performed from which it was concluded that an optimum combination of accuracy and the number of elements is obtained at N dof = 4100 (approximately 8050 linear elements). Relative errors of approximately 10 −4 − 10 −5 are usually found within 10-15 iterations. Typically, the computation time is within 10 s on a quad core desktop computer. Solving Equation (22) using a typical PDE routine gives a reasonable solution, as the overshoot for the mass fraction distribution is almost suppressed and the sharp transition at the edge of the bearing is retained; see Figure 7b. Several other non-linear crosswind numerical stabilization techniques have also been tested, but the one represented by Equation (21) has been shown to be most stable and the fastest to converge.
The artificial diffusion term for the crosswind stabilization slightly changes the original Reynolds PDE. In order to see if the results after this modification are still accurate, our cavitation algorithm was compared with those from Yu et al. [26] for a fully submerged bearing. Yu transformed Elrod's universal differential equation into a boundary integral equation by means of boundary element theory and compared their results with those obtained using a finite difference approach in which they found excellent accordance. Typical Dirichlet conditions (P = 0) are imposed on all edges of the computational domain. The pressure and mass fraction distributions are compared in Figure 8a,b, respectively, from which it is clear that the results closely follow each other. As mentioned earlier, the fractional film content in the full film region of the reference results is modeled using compressibility, which is the reason why the non-dimensional density slightly increases in that region.
Based on this comparison, it can safely be stated that the artificial diffusion terms have little influence on the overall solution accuracy, and thus, the AD-method may be used.

2D Transient Situation: Equations of Motion
For a time-dependent problem, the squeeze term ( ∂(f H) ∂t ) needs to be added to Equation (22). The stabilized 2D cavitated-Reynolds equation then becomes: with the dimensionless circumferential and axial flows, respectively: An inverse bearing model is used here: it is assumed that the applied load W a and load angle X W are known, and the corresponding eccentricity ε and attitude angle X ε have to be calculated. Note that the load vectors are dependent on the pressure and therefore on the eccentricity and attitude angle.
In case an external force (W a ) is acting on the shaft, this force and the load carrying force of the lubricating film combine, and the resulting equations of motion become: where M represents the reduced mass of the system. M , Wx, Wȳ and X W are given by: Equations (25) and (26) are simultaneously solved using a modified Newton-Raphson iteration scheme from which the correct load angle and eccentricity are obtained.

Expansion to a Thermo-Hydrodynamic Model
In the previous section, a novel and fast mass-conserving Reynolds cavitation algorithm was developed, from which the mass fraction and the pressure distribution throughout the fluid film domain can be calculated. In general, modeling a plain journal bearing using the assumption of an isoviscous fluid leads to unrealistic results [4]. Hence, in this section, a simple and efficient thermal model is developed from which the viscosity distribution can be obtained. Furthermore, in the previous section, it was assumed that the supply groove was infinitely small in the circumferential direction. In this section, the shape of the supply groove and its influence on the temperature and pressure distribution is taken into account.
In addition to the variables of Equation (4), the following dimensionless variables are introduced:

Fluid Rheology
In this model, the viscosity is a function of temperatureT based on the following exponential expression:η = e −βT (28) whereβ stands for the dimensionless temperature viscosity coefficient. The shear rate and pressure dependencies of viscosity are not taken into account in this analysis, but could straightforwardly be taken into account in the methods presented in this section.

Energy Equation and Oil Supply
The heat transport calculations are based on the assumption of adiabatic conditions, i.e., heat transport to and from the bearing housing and the shaft are neglected. This follows from the assumption that convective heat transport is dominant, i.e., the lubricant leaving the bearing carries away the major part of heat generated in the film. Furthermore, the use of an effective temperature across the film thickness is made. This reduces the 3D energy model to a simple and efficient 2D model. From these simplifications, the adiabatic energy equation reads [27]: where q X and q Y are the bulk flow velocities in the X and Y -directions, respectively, obtained from Equations (24a) and (24b). The viscous dissipation termψ and friction lossesφ can be calculated as follows:ψ where Ω stands for the bearing surface area. The second term on the RHS of Equation (29) denotes the inflow (heat flux) of oil from the supply groove(s) to the bearing. The calculation of this term follows later on in this section. Equation (29) assumes that finger-type cavitation occurs, which means that the flow in the cavitation area consists of fluid streamers separated by pockets of gas. This means that in the cavitated region, heat conduction in the Y -direction can safely be ignored, as heat conduction through gas is far less than through oil. Hence, the Boolean expression in front of the heat conduction term is in the Y -direction. Heat generation in the cavitated regions is assumed to be solely generated due to shear losses in the streamers. Furthermore, the heat capacity k and heat conduction coefficient C p are assumed to be constant throughout the calculations, as their temperature and pressure dependency is much less compared to that of the viscosity.
Similar to the Reynolds equation, the thermal Equation (29) is strongly convection dominated and, hence, also requires numerical stabilization. To this end, we have implemented the SUPG method in order to stabilize the solution of Equation (29) in the cavitated region.
The temperature distribution in the supply grooves is calculated based on perfect mixing conditions (see, for instance, van Ostayen and van Beek [27]). The oil inflow from the supply grooves to the fluid film is calculated, taking into account the shape of the groove, by: where Ω m stands for the surface area of the groove. In some cases, the hydrodynamically-generated pressure in the supply groove is larger than the supply pressure. In those situations, hot oil is forced back into the system (Q m < 0), and the Boolean function in Equation (31) is equal to zero, as no heat transport takes place. Weak Dirichlet pressure boundary conditions (P = P s ) are imposed on the edges of the supply grooves. The volumetric flow in the groove can then be obtained by integrating the Lagrange-multiplier associated with the pressure, over the edges of the supply groove: where ξ lm represents the Lagrange-multiplier associated with the pressure P . The assumption has been made that the oil inflow from the supply grooves to the fluid film is distributed across the grooves relative to the local groove depth according to the following relation [27]:

Boundary Conditions
At X = −π and X = π (BC1), periodicity is imposed to ensure flow continuity. On the edges of the supply groove, a typical Dirichlet boundary condition has been imposed with supply pressure P s . On the edges of the bearing (BC2), non-submerged boundary conditions are imposed, and the temperature gradient normal to the edges is assumed to be zero. The term "non-submerged" mathematically implies that in the full film areas, the pressure on the edges is set to zero (ξ = 0), as there is outflow of lubricant, whereas in the cavitated/inflow regions, zero flow is enforced (no boundary condition is set, 0 = 0). The non-submerged boundary condition can thus realistically be simulated using the following constraint: where n Y is the normal vector of the edges.

Results
The developed model in the previous sections is solved using a general purpose finite element analysis software [12]. Notice the strong non-linear coupling between Equations (2), (6) and (26), where W and X W are a function of P and, therefore, a function of ε and X ε . Moreover, the variable viscosity (η) couples the Reynolds equation with the energy equation.
Summarizing, the set of equations and constraints (load balances, the Reynolds and energy equation with its respective boundary conditions) is solved using this finite element formulation, where the resulting system of non-linear equations is solved using a monolithic approach, where all of the dependent variables ξ,T , ε, X ε are collected in one vector of unknowns and simultaneously solved using a modified Newton-Raphson iteration scheme coupled with a time stepping algorithm. For a detailed description of the used time stepping algorithm and integration scheme, see [12].

Validation of the 2D THD Steady State Model
Ferron et al. [14] performed experiments on a single axial grooved bearing and obtained satisfactory accordance between simulations and measurements. Their measurements are regularly used as a benchmark for verification of THD models (for example [28][29][30]). In order to validate the proposed model, the bearing geometry and operating conditions used by Ferron et al. [14] are used for a numerical study (see Table 1). In the present work, the assumption of adiabatic bounding surfaces has been made after consideration of the operating conditions under which the experiments were performed. A convergence analysis was performed for the reference operating conditions, from which it was concluded that an optimum combination of accuracy and the number of elements is obtained at N dof = 4100 (approximately 8050 linear elements), and hence, this mesh distribution is used further throughout the analysis. Relative errors of approximately 10 −4 − 10 −5 are usually found within 15-25 iterations. Typically, the computation time is 20-25 s on a quad core desktop computer.
It is worth mentioning that in the paper of Ferron et al., the viscosity-pressure dependence was not taken into account, although this dependence may compensate in part for the viscosity reduction due to the temperature rise. The pressure viscosity coefficient for a typical mineral oil is in the range of 10 −8 Pa −1 . Thus, it is safe to ignore the viscosity-pressure dependence (as compared to the viscosity-temperature dependence) considering the order of magnitude of the pressures in Figure 9b.
In Figure 9a, the dependent variable ξ is plotted, in which the hydrodynamic pressure peak can clearly be observed. It can also be observed from Figure 9a that due to the imposed non-submerged boundary condition at the edges of the bearing (see Equation (34)), ξ = 0 in the full film region, and no flow is enforced in the cavitated region.
The theoretical predictions are compared with the measurements of Ferron et al., in Figure 9b,c. Concerning the pressure profiles, excellent agreement is found. Note the jump in the pressure profile at both ends of the interval as a consequence of the imposed pressure constraint P s at these ends (the supply groove has a finite width, which is defined by the groove angle (see Table 1)). There are, however, some discrepancies between the predicted and measured temperature profiles. As can be seen, the magnitude of the predicted maximum temperature is in accordance with the measured maximum temperature. The location of the maximum temperature is not accurately predicted, though. A probable explanation follows from the treatment of heat transport in the cavitated region. In this model, the assumption was made that the cavitated area consists of fluid streamers (finger-type cavitation). However, the actual behavior of fluid in the cavitated region may be different. Moreover, the oil inflow profile, which was assumed, Equation (33), might be inaccurate. It is known that the inflow of oil can have significant influence on the calculated temperature field and, hence, should be investigated in more detail.
In Figure 10a, the eccentricity ratio is plotted against the varied applied load: fair agreement is found between the measurements and the predictions. Satisfactory agreement is also found between the outflow of the lubricant as a function of the eccentricity ratio; see Figure 10b. A similar statement can be made about the maximum pressures as a function of the eccentricity ratio at different rotation speeds; see Figure 10c,d.
In Figures 10e,f, the maximum temperatures are plotted against the eccentricity ratio at different rotational speeds. As can be noticed, there are some discrepancies between the predictions and theory of the maximum temperatures. This might be due to the aforementioned assumptions of the cavitation type and the oil inlet profile. In addition, changes in clearance of the bearing due to thermal growth could also have caused discrepancies, which are not taken into account in our model. In general, it can be stated that satisfactory agreement between the measurements and the predictions is found.
(e) (f) Figure 10. Comparison of the modeling results with the measurements of Ferron et al. [14].

Transient 2D THD Results
Although the reference results of Ferron [14] only contain steady-state results, a time transient analysis is performed here for the sake of completeness. The transient analysis is performed by means of applying a synchronous harmonic load on the journal, to demonstrate the performance of the journal bearing model under dynamic loading. The applied load in its respectivex-andȳ-directions can be written as: In Figure 11a, the shaft orbit resulting from the transient analysis is plotted. It is observed in Figure 11a,b that after, the initial transient effects, both pressure p and film thickness h, reach a limit cycle and continue to move in that cycle. Moreover, the solution for the maximum pressure and temperature in the limit cycle is in the vicinity of the steady-state solution (see Figures 11b and 12). An interesting result that is provided is the behavior of the maximum temperature T max with time, which is that it reaches equilibrium with the surroundings and, hence, approximates a steady value after a certain period of time. The typical computational time for the time-dependent calculation is 2-4 minutes per revolution on a typical quad core desktop computer.
Careful study of Figure 11a reveals that it is clear that the shape of the profile (ε, X ε ,t) is not a perfect circle, which means that the eccentricity (ε) has a maximum and a minimum (in this case, near an altitude angle (X ε ) of 180 • ). This was found to be due to the hydrostatic pressure of the oil inlet channel. As can be observed in Figure 11b, a peak pressure occurs every time the hydrodynamic peak pressure coincides with the start and end of the oil inlet channel, explaining the double pressure peaks each revolution. During this passage, also the supply flow decreases, which logically implies that less cold lubricant is introduced into the film (see Equation (31)) and, hence, results in a slight increase in temperature, which, after the passage of the oil inlet channel, decreases again to its former value ( Figure 12).  and ω s 2 krpm.

Conclusions
In this paper, a simple and efficient thermo-hydrodynamic (THD) model has been developed and applied to a plain journal bearing. To this end, a new mass/conserving Reynolds cavitation algorithm has been developed for finite element analysis. The proposed pressure and mass fraction transformation, in combination with anisotropic artificial diffusion, provides a fast, stable and robust approach to solve full film lubrication problems. The complete THD model simultaneously solves the fundamental equations (Reynolds and energy equation) combined with a set of load balance equations and has proven to give results in accordance with published measurements for steady-state cases. To be more specific, design parameters, such as the eccentricity ratio, the maximum temperature and pressure, are satisfactorily similar to the measurement results. However, some discrepancies were found in matching of theoretical and predicted circumferential temperature profiles. A probable explanation was given by the treatment of heat transport in the cavitated region, the inflow profile assumed and the thermal growth of the components. Hence, these should be investigated in more detail. Finally, a transient analysis of the shaft under harmonic loading was performed, which revealed distinct knowledge about the dynamic and transient thermal behavior of the bearing.

Author Contributions
Shivam Alakhramsing and Ron van Ostayen developed the model. Shivam Alakhramsing performed the validation studies and wrote the manuscript. Rob Eling coordinated the present work. x, y, X, Y coordinate along resp. perpendicular to the sliding direction ((m), (-)) X ε , X W attitude, load angle (rad)