Lubrication theory for Bingham plastics

Abstract The rheological behavior of many lubricants used in hydrodynamic bearings can reasonably be modeled using the Bingham plastic material model. This behavior is characterized by a strong discontinuity, from a pure solid state to a viscous fluid state depending on the local shear stress. In literature three methods have been presented to model a Bingham plastic lubricated film. A full CFD and thus expensive, numerical simulation has been used. A general Reynolds equation based simulation has been used, however with a less accurate numerical regularization of the material discontinuity. Or a general Reynolds equation based simulation has been used, but with a severe reduction of the geometric complexity. In this paper, an ’exact’ thin film lubrication simulation for a Bingham plastic fluid is presented. The model is said to be exact in the sense that it requires no additional approximations to those used in the derivation of the general Reynolds equation, and requires no numerical regularization of the Bingham plastic fluid model and can still be used on any lubricating film geometry. Simulations on both infinite and finite journal bearings shows that the results of this new method are in good accordance with literature, demonstrating the validity of the method.


Introduction
The hydrodynamic bearing is a machine component used to transmit a load between two parts of a mechanical system moving with respect to each other. In this type of bearing a lubricating film separates the moving surfaces. The fluid flow in this relatively thin lubricating film can in general be modeled using full 3D fluid flow relations or preferably, by using a numerically much more efficient 2D thin film approach. This thin film approach is based on a list of approximations first derived by Reynolds (1886), with the main observations being that the flow is dominated by viscosity and that the pressure in the film is constant in the thickness direction of the film. These observations and approximations result in the so-called general Reynolds equation, a 2D partial differential equation describing the pressure in the film [6].
In many bearing systems lubricants are used of which the rheological behavior is reasonably modeled using the Bingham plastic material model, which combines a solid behavior at shear stress values below a critical yield stress, with a constant viscosity above that yield stress. Some frequently used applications are, for example, the grease in roller bearings or low speed journal bearings. Some less common applications are, for example, magnetorheological and electrorheological bearing systems. The rheological properties of these fluids are a function of the magnetic and electric field. Using this property it is possible to realize an active bearing by controlling the rheological properties of the lubricant [12,13,21,25]. Proper design of these bearing systems is only possible when an efficient and accurate numerical simulation method is used to predict the behavior of a certain design. The discontinuity in the viscosity for the Bingham plastic fluid model caused by the yield stress has been shown to be a major challenge here.
Wada et al. [23] was the first to present a significant lubrication theory for Bingham plastics in 2D lubricating films. The major drawback of this model is that it is an approximation in the sense that it requires additional approximations to accommodate for the strong non-linearities inherent in the Bingham material model. The theory was validated for the case of a step bearing [22] and a journal bearing [24]. The method of Wada et al. was later implemented in Refs. [4,14]. In Ref. [4] the calculations are limited to situations where no plug flow occurs, significantly reducing the numerical complexity but also the practical usefulness.
In the work of Tichy [19] a standard lubrication model has been used in combination with the Bingham plastic fluid model, however only applied to 1D lubricating films. The downside of the method is that it is only valid for 1D lubricating films which makes it only accurate when applied to for example infinitely long journal bearings, or infinitely wide slider bearings. The benefit is that the method can be said to be exact since it does not use any additional approximations to accommodate for the non-linearities due to the yield stress. The method of Tichy [19] was later implemented in Refs. [7,20].
The seminal work of Dorier and Tichy [5] presents a regularization strategy for the lubrication theory for Bingham plastic fluids applicable to a 2D lubricating film. The disadvantage of the method is that it approximates the Bingham plastic fluid model by using a model that is similar to the Bingham-Papanastasiou model [15]. This approximation converges to the conventional Bingham model for an increasing value of the regularization parameter but simultaneously gets more computationally demanding. The advantage of the method is that it removes the necessity of explicitly tracking the location of the plug in the film as is the case in work of Tichy [19] and Wada et al. [23]. The method of Dorrier and Tichy [5] was later implemented by for example Sharma and Khatri [18], Jang and Khonsari [9] and Khlifi et al. [10].
The work of Gertzos et al. [8] presents a full CFD method to model the behavior of a Bingham Plastic in a lubricating film. It uses the Bingham-Papanastasiou model to approximate the Bingham behavior. The drawback of this theory is that it can be said to be very computationally demanding since it desires the solution of all fluid flow variables (velocity vector and pressure) on a full 3D domain and it includes a regularization technique. The advantage of the method is that it is also applicable when the thin film approximation is not valid. The method of Gertzos et al. [8] was later implemented by Peng and Zhu [16,17]. In which the authors furthermore extended the method with a cavitation algorithm and a varying yield stress over the lubricating film.
The works of Bompos et al. [2,3] and Kim et al. [11] demonstratethe use of a bi-viscous material model implemented in a CFD model to mimic the behavior of the Bingham plastic material in a lubricating film. This model uses an artificial and relatively high viscosity to model the solid phase behavior and the physically accurate low viscosity to model the fluid phase behavior.
In literature three methods to model a lubricating film with a Bingham plastic fluid have been published. The first method is to use the full 3D flow relations to model the flow in the full 3D domain of the lubricating filmThe second method is to use the approximated thin film model presented by Wada et al. [23]. The third method is to use the regularization technique used by Dorier and Tichy [5]. An additional fourth method exists when modelling only a 1D lubricating film. The work of Tichy [19] presents an exact lubrication theory that is applicable to a 1D lubricating film. The overview demonstrates that all presented methods that are applicable to 2D lubricating films are either computationally demanding or provide only an approximation to the exact but strongly discontinuous Bingham behavior. Therefore, this paper presents an exact lubrication theory for a Bingham plastic fluid that is more computationally efficient. The theory is said to be exact in the sense that it requires no additional approximations other than the ones already used in the derivation of the generalized Reynolds equation [6]. Compared to other published methods, the method is less computationally demanding due to the lower number of degrees of freedom needed and absence of a regularization method. The method is tested against results presented in literature. Simulation results both for infinite and finite journal bearings are obtained.

Method
The method description consists of two parts. In the first part the lubrication theory for a Bingham plastic fluid is derived. In the second part the numerical method used to solve the set of equations derived in the first part is discussed. Data from literature (Wada et al. [24], Tichy [19] and Gertzos et al. [8]) is used to compare the modelling results for both finite and infinite journal bearings. The data was gathered from these sources using a data extraction utility [1].

Bingham material model
The flow of a Bingham plastic fluid is described using a model where the material behaves as a solid when the shear stress magnitude in the material j τ ! j is lower than a given yield stress τ 0 . The material behaves as a liquid when the stress in the material is higher than that yield stress. Eq. (1) presents the formula that describes the material behaviour for this paper, it shows that the material is solid below the yield stress τ 0 and flows above the yield stress.
In Eq. (2) we introduce the proportionality between shear rate vector and the shear stress vector which includes flow factor f and viscosity η.
Eq. (3) describes this flow factor which for all values of the shear stress and yield stress has a value between 0 and 1.

Stress distribution
The theory considers only lubricating films and is valid for the approximations laid out in Ref. [6]. The Cauchy momentum equation is thus reduced to the form presented in Eq. (4), where p stands for the local film pressure and z presents the z-coordinate in the film (see also Integration of Eq. (4) over the film height results in Eq. (5) that describes the stress as a function of the film height and the stress τ ! c at z ¼ 0. Fig. 1 provides a graphical presentation of the distribution of the xcomponent of the stress across the film height.
The magnitude of the shear stress vector is given by Eq. (6). Fig. 2 presents a graphical representation of the shear stress magnitude j τ ! j.
The blue part of the graph presents the amount of stress that causes flow since this is the only part above the yield stress. The green part represents the part of the shear stress that does not cause flow since it is below the yield stress. The red part represents the amount of additional shear stress that would be required in order to cause flow. � � τ x-component of the shear stress vector.

Plug analysis
Between the two surfaces part of the material behaves as a solid and part as a fluid, depending on the shear stress distribution. Due to the nature of this shear stress distribution there are two possible transition points (z bottom and z top ) where the lubricant changes from solid to fluid. The lubricant behaves as a solid between these two transition points and as a fluid outside these two transition points. In order to locate these points, first the z-coordinate of the minimum shear stress z m and the corresponding minimum stress value τ m are derived from Eq. (6) and presented in respectively Eq. (7) and Eq. (8).
Solving Eq. (9) for z when the stress is equal to the yield stress τ 0 , leads to Eq. (10). The parameter Δ presents the distance from the point of minimum stress z m to the point where the stress is equal to the yield stress, see Fig. 2 for a graphical presentation of this.
Eq. (10) now finally leads to the location of the two transition points z bottom and z top as given by (11) and (12) The max and min functions in Eq. (11) and Eq. (12) ensure that the resulting value is between h 2 and h 2 , and therefore somewhere inside the extent of the film.

Flow field
The lubricant velocity across the film height can be found by integrating Eq. (2) over the height of the channel, resulting in Eq. (13). Fig. 3 gives a graphical presentation of a typical fluid profile. The green and blue parts represent the lubricant above and below the plug where viscous flow occurs, respectively. The red part represents the solid plug.
For convenience we define the flow factor integrals Eq. (14) such that the velocity difference between the two surfaces is described using Eq.
Substituting Eq. (16) into Eq. (13) results in Eq. (18) which describes the flow velocity across the height of the channel as a function of the surface velocities.

Flow rate
Integration of Eq. (13) across the height of the channel results in Eq. (18) that describes the flow rate through the system.
Using integration by parts to simplify this result leads to Eq. (19), in  which the flow rate is described as a function of the as yet unknown minimal shear stress τ ! c . In this step this derivation deviates from the derivation of the general Reynolds equation in Ref. [6].
More traditionally, by substituting Eq. (16) into Eq. (19), we obtain Eq. (20), in which the flow rate is expressed using the surface velocities.
Both Eq. (19) and Eq. (20) describe the flow rate through the system and can as such be used to solve for the one remaining unknown, the film pressure p, using the conservation of mass in Eq. (21). However, this does require a solution for the flow factor integrals F n . This is described in the next section.

Flow factor integrals
Goal of this section is to get rid of the integral in the flow factor integrals F n defined in Eq. (14). Eliminating these integrals removes the necessity to do a numerical integration step as required for instance in the regularization method used by Tichy [19]. To achieve this, Eq. (3) and the results of Eq. (11) and Eq. (12) are substituted in Eq. (14) resulting in Eq. (22). where: and: or, after introducing the following notation short-cuts: this integral is equal to: Z B A ðZ þ bÞ n ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi a 2 þ ðcZÞ 2 A Z k ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi a 2 þ ðcZÞ 2 where the integrals H k can be determined analytically: � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi a 2 þ ðcBÞ 2 q ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi a 2 þ ðcAÞ 2 q � (28) ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi a 2 þ ðcBÞ 2 q A ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi

Numerical model
The system of differential equations that describes the behavior of a Bingham plastic fluid in a lubrication film is solved using the commercial finite element package COMSOL Multiphysics® 5.4.
The lubricating film in a journal bearing with one axial supply line is considered. The supply line is located at the top of the bearing ðθ ¼ 0Þ with a supply pressure equal to ambient pressure. The height of the film is described by Eq. (30) where c is radial clearance, e the eccentricity and θ the radial bearing position. The pressure is constrained at all four sides to ambient pressure by using a Dirichlet boundary condition.
The model uses a triangular mesh with a maximum element size of 5% of the bearing length. The elements use a Lagrange shape function and are quadratic order for the PDE and quartic order for the weak contribution.
Two solver strategies are tested within this work: Strategy A uses Eq. (19) to solve the system as a partial differential equation for the pressure p while implementing Eq. (15) as a weak contribution to force the correct surface velocities.
Strategy B uses Eq. (20) to solve the system as a partial differential equation for the pressure p. The values of F 0 , F 1 and F 2 are updated every iteration based on the current pressure distribution. The process can be summarized with the following solver sequence: 0. Initial guess: [1.], iterate, until convergence criteria are met Both solver strategies use a parametric solver to guarantee convergence of the non-linear system. The solver starts with a zero yield stress and iterates towards the desired yield stress. The step size is generated a z n dz ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ðZ þ z m Þ n dZ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi � automatically by the software package. The solver is assumed to be converged when a relative tolerance smaller than 1e-6 is achieved. Computation is performed on an Intel Xeon CPU E5-1620 V3 @ 3.50 GHz with 32 GB of Ram.

Results
Fig . 4 shows the normalized pressure distribution p � ¼ h 2 0 p ηuL in a finite journal bearing for a Newtonian fluid for different eccentricities e. The figure shows the results published by Wada et al. [24], Gertzos et al. [8] and the present work. Fig. 5 shows the pressure distribution in an infinite journal bearing both for a Newtonian fluid and for a Bingham plastic fluid, where the normalized yield stress is given by relation T 0 ¼ τ0c uη . The figure shows the results of Tichy [19] compared to the current study.
Figs. 6 and 7 present the pressure distribution in a finite journal bearing lubricated with a Bingham plastic for different values of T 0 and different bearing eccentricities. The figure shows the work published by Wada et al. [24], Gertzos et al. [8] and the present work of this paper. Fig. 8 and Fig. 9 present the plug configuration derived with the method presented in this paper for a normalized yield stress of T The models typically solves within a minute for one bearing configuration. The solution process tends to be relatively slow for large plug configurations and relatively quick for no or small plug configurations. No significant difference between method A and method B concerning speed or accuracy were found during this research. Fig. 4 shows that for a Newtonian fluid the method of Wada et al. [23] corresponds exactly with the method presented in this work. This good fit is expected since in this case the two methods use exactly the same approximations; here the two methods are fundamentally the   Fig. 6. Pressure distribution for a finite journal bearing for the case.T 0 ¼ 0:4 Fig. 7. Pressure distribution for a finite journal bearing for the case.T 0 ¼ 0:8 same. The figure shows furthermore a slight offset between the method of Gertzos et al. [8] and the method presented in this paper. This is explained by the fact that a full 3D CFD simulation is compared with a 2D thin film approximation. Fig. 5 shows that the method of Tichy [19] is in good accordance with the method presented in this paper. Extracting the data from the original journal paper caused a slight offset between the two lines. It must be noted here that this analysis is only an infinite journal bearing analysis. This means that the shape of the fluid profile is more basic and there is no cross-coupling between the flow in x-and y-direction.

Discussion
Figs. 6 and 7 show a good fit between the methods of Wada et al. [24], Gertzos et al. [8] and the one presented in this paper. The minor offset between the CFD analyses of Gertzos et al. [8] and the present method is also there for the Bingham plastic fluid. Interesting to note is that the pressure profile obtained using the present method always seems to exceed the pressure profile compared to Ref. [8]. The difference is again explained by the nature of the simulation, one is full turbulent 3D CFD while the other is a laminar thin film approximation.
The close fit between the method of Wada et al. [24] and the one introduced in this paper is not there any more for a Bingham plastic fluid. This can be explained by the additional approximations on the fluid flow required for the derivation of the method of Wada et al. [24]; the two methods are therefore not fundamentally the same any more. Since the method in the present work is based on only the classic thin film approximations combined with the Bingham plastic fluid model, it is considered to be the more exact method. All three lines are in good accordance with the measurements performed in the work of Wada et al. [24].
Figs. 8 and 9 demonstrate that the resulting pressure field and plug configuration have a smooth and symmetrical solution. Note that cavitation is not included in this model.
The computation power needed for the present method is comparable to that of the method published by Wada et al. [23]. Both methods are of a similar form and require a similar solver strategy. The present method is less computationally demanding than the method of Dorier and Tichy [5] since it does not require an numerical integration step in the solver sequence to get the proper viscous behavior. The present method is also significantly less computationally demanding than the method of Gertzos et al. [8] as this is a full 3D turbulent CFD calculation. To the contrary, the method proposed in this paper uses the thin film approximations and obtains a solution for a laminar pressure distribution in 2D, resulting in a huge reduction in degrees of freedom.

Conclusion
This paper presents an efficient numerical simulation method to describe the pressure in a lubrication film for Bingham plastic fluids. The derivation uses the 'exact' Bingham plastic fluid rheology model in combination with only the approximations required for the general Reynolds equation. A further approximation of the Bingham model or regularization method is not required. Also, the method is computationally efficient since no regularization technique is needed. The results are shown to be in good agreement with literature.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.