# I. Hydrodynamic Model for the Journal Bearing he journal bearings are sliding bearings which use the flow of a lubricating fluid around a pair of non concentric circular surfaces given by the shaft and the bearing that generate a pressure field. Under high rotations this pressure field supports the shaft, eliminating its contact with the surface of the bearing. The pressure field arises due to the variation of the clearance between these surfaces, as shown in Fig. 1: The most common mathematical model to approach this type of problem, usually proposed in the design of machine elements, is deduced from the Navier-Stokes Equation assuming several simplifying hypotheses. The following features are neglected: the curvature of bearing, the inertia of the lubricant, the external field forces, the gradient of pressure and the velocity of the fluid in the radial direction and their variation in this direction as well. Newtonian lubricant with constant viscosity is also assumed, as well is supposed laminar and incompressible flow. Thus, under conditions of one-dimensionality, the following differential equilibrium equation is given by (Shigley et. al, 2003): dx dh(x) V 6 dx ) x ( dp ? ) x ( h 3 dx ) x ( p d ? ) x ( ?h ] x d p(x) d ? ?h(x) [ x d d 2 2 2 3 3 ? = ? + = (1) In Eq. (1) p(x) is the overpressure with respect to the initial oil injection pressure, V is a mean flow velocity and ? is the viscosity. The clearance h(x) Author ? ? ?: Federal University of Espírito Santo -Technological Centre, Av. Fernando Ferrari, 540 -Bairro Goiabeiras, CEP 29075 -910, Vitória, ES -Brazil. e-mails: carlosloeffler@bol.com.br, depends on both the curvature and the eccentricity between the surfaces of the components, but it is known, see figure 1. Essential conditions are imposed in terms of the pressure p(x) and natural conditions in terms of the flow q(x), as follows: ?h(x)v(x) q(x) = (2) T Global Journal of Researches in Engineering ( ) In the previous equation, ? is the density, assumed to be constant. The lubricant velocity v(x) is given by: dx dp(x) dx dh(x) ? h(x) V v(x) 2 ? = (3) Equation ( 1) is a particular case of the Reynolds Equation in fluid motion. Despite in many journal bearings, the fluid flow in the z direction does not exist, the pressure can vary axially as a consequence oil output system. In this case, a source term W(z) can be assumed, resulting in a more general governing equation (Panday et. al, 2012): ) z ( W dx d[h(x)] V ] z z) p(x, ? ) x ( h [ z ] x z) p(x, ? ) x ( h [ x 3 3 ? ? = ? ? ? ? + ? ? ? ? (4) Thus, in this assumed model the lubricant can have a similar magnitude in both directions x and z. Concerning the source terms in the right hand side of Eq. ( 4), both result due to the simplifications imposed in the mathematical model, especially related to the integration along y-direction. Another kind of sources or sinks can be idealized, as those that concentrate in a specific region. It must be highlighted that the height h(x) does not change axially, what allows an easier numerical approach. Equation ( 4) it is a partial scalar differential equation in which the physical property continuously varies only the position x in the domain ?. It can be verified that this governing equation is similar to that which describe the following cases: torsion in bars with uniformly variable sections; potential flow in tanks with variable depth; and heat transfer in heterogeneous media. # II. # Integral Governing Equation For convenience, henceforth the indicial notation is used for easiness in subsequent mathematical manipulations. Therefore, the governing equation (Eq. ( 4)) can be adequately written in indicial notation, as shown below: s(X) ], [K(X)p(X), i i =(5) In this last expression it was considered the following definitions: ) z , x ( X X ; W(z) dx d[h(x)] V ) X ( s ; ? ) x ( h ) x ( K 3 = ? ? = = (6) The starting integral form (Brebbia, 1978) is obtained by integrating Eq. ( 5) over the physical domain ?(X), using as an auxiliary function the fundamental solution u*(?; X). This function is the solution of a Poisson´s problem which is assumed a homogenous infinite medium, submitted to a unitary concentrate source applied to an arbitrary point ?, that is (Brebbia et al., 1980): X) ; ( X) ; *( u, ii ? ? ? = ? (7) Thus, the integral equation associated with the Eq. ( 5) is given by: ? ? = ? ? ? i i ? X)d? ; ( * s(X)u X)d? ; ( * u ], p(X), [K(X)(8) The inverse integral form of the left hand side of Eq. ( 8) can be deduced by performing two integrations by parts, as shown below, where the Divergence Theorem has been properly applied: d? ], X), ; ( * p(X)[K(X)u [ d ]n X), ; ( * p(X)K(X)u [ d X)]n ; ( * u p(X), [K(X) X)d? ; ( * u ], (X)p(X), [K i i ? i i Î?" i i Î?" i i ? ? ? + Î?" ? ? ? Î?" ? ? = ? ? (9) Developing the product of two functions in the last term to the right hand side of the previous equation can be rewritten: ]d? X), ; ( * p(X)K(X)u [ ]d? X), ; ( * (X)u p(X)K, [ d? ], X), ; ( * p(X)[K(X)u [ ii ? i i ? i i ? ? ? + ? ? = ? ? (10) Using the properties of the Dirac Delta function (Raisinghania, 2011), one has: ) )p( )K( c( ]d? X), ; ( * (X)u p(X)K, [ X)]dÎ?" ; ( * [p(X)K(X)q X)] ; ( * q(X)u {[K(X) X)d? ; ( * u ], p(X), [K(X) i i ? Î?" Î?" i i ? ? ? ? + ? + ? ? ? = ? ? ? ? ? (11) In the previous expression, q*(?; X) is the normal derivative of the fundamental solution and q(X) is the normal derivative of the pressure p(X). The coefficient c(?) depends on the positioning of the source points ? with respect to the physical domain ? (X) and, in the case of being located on the boundary Î?"(X), also of the smoothness of this one (Brebbia and Walker, 1980). It is verified the presence of a domain integral in Eq. ( 11), which is treated in the following item, like the approach of the term source s(X), referring to the right side of Eq. ( 6). It should also be noted that in this case K(X) is treated as a nodal quantity, linearly interpolated on each boundary element. # III. # Quasi-Dual Procedure Focusing the elimination of the domain integrals that still persist in the right hand side of Eq. ( 11) the direct integration procedure (DIBEM) could be used. However, in this case, the Quasi-dual Reciprocity model (Loeffler and Mansur, 2003) is the most suitable and accurate procedure. Quasi-dual (QDR) is a technique similar to Dual Reciprocity (Partridge et al., 1992), but it was developed to approximate first order derivatives in diffusive-advective problems. The QDR approximation uses linear combinations of primitive radial basis functions ? j , which are multiplied by coefficients ? j in the following form: X) ; (X ?, ? (X) p(X)K, j j i j i ? (12) In Eq. ( 12) X j is the coordinates of interpolation basis points. Considering two-dimensional cases the QDR approach is just suitable for potential fields. However, since the height h does not vary in the z direction of journal bearings (see eq. ( 6)), the differential given by the left hand side of Eq. ( 12) is always exact and can be treated accurately using the QDR. Thus, replacing Eq. ( 12) in the domain integral that persists in Eq. ( 11), one has: ) ; X ( ? ? ) ( c X)]dÎ?" ; ( * X)q ; (X ? [ ? ]d? X), ; ( * X)u ; (X ? [ ? d? ], X), ; ( * X)u ; (X ? [ ? ]d? X), ; ( * u X) ; (X ?, [ ? ]d? X), ; ( * (X)u p(X)K, [ j j j j j Î?" j ii j j ? j i i j j ? j i j j i ? j i i ? ? ? ? ? ? = ? ? ? ? ? = ? ? ? ? ? (13) As a last remark, it is important to highlight that unlike other models that use radial basis functions, the inclusion of poles did not improve the results in QDR approach; in fact, the inclusion of poles only increased the round of errors. IV. # Dibem Procedure When the body force term is given by the simple mathematical function, the Galerkin Tensor is the most effective way to deal with its domain integral, transforming it in a boundary integral. Considering more elaborate functions, a solution of Poisson's problems can be more effectively by global functions in the Goldberg's sense (Goldberg,1994) or then, using the Radial Integration Method (RIM) (Gao, 2002), despite the huge computational time of this later. Here, the DIBEM procedure is applied in association with the QDR for convenience. DIBEM is also a similar technique to the Dual Reciprocity but the entire kernel of the domain integral is interpolated, aiming to transform it in a boundary integral. The technique already successfully applied to scalar problems governed by the Poisson equation and Helmholtz . DIBEM substitute advantageously standard procedures as the domain integration by cells and the DRBEM. Comparatively, it shows superior performance and mathematical suitableness, since it is more similar to a simple interpolation procedure. First DIBEM tests for performance are done solving Poisson's Equation, where domain integral is just comprised by known functions. Thus, different coordinates can be used to distinguish the interpolation points to the field points on the boundary, avoiding the singularity in the fundamental solution. Concerning the internal interpolation points, these points do not appear in the final matrix system. A different situation occurs in other more elaborate problems as Helmholtz and diffusive-Advective, in which the regularization procedure is required to avoid singularity since the Provided internal interpolation points are not necessary for the QDR, the application of DIBEM is easier. Thus, the complete kernel of the domain integral is interpolated directly, according to the following expression: ) X ; X ( F ) X ; ( u ) X ( s j j j ? = ? ? * (14) Similarly to the DRBEM, the proposed method also uses a primitive function ?, such as: ) X ( z d ) X ( d )) X ( n ) X ( ( d )) X ( ( d )) X ( F ( d ) X ; ( z j j j j i j i , j j ii , j j j ? = Î?" ? ? ? = Î?" ? ? ? = ? ? ? ? = ? ? ? = ? ? ? ? Î?" ? Î?" ? ? ? ? ? ?( # Discretization Using radial basis functions, for each source point ?, the interpolation given by Eq. ( 6) corresponds to scanning all points X j in relation to domain points X. For the QDR, a similar procedure to the DRBEM is followed, in which the matrix H already calculated by the discretization of the integrals related to the Laplacian is used (Partridge, 1992). Double points X j located in the corners should be departed to avoid singularity in the inversion of the interpolation matrix. Thus, one can write: p ] H [ ]p , [K ] , [ ] ][ H [ ] ][ H [ q ] G [ p ] H [ i -1 i = ? ? = ? ? + ? (17) In the last equation the ? vector was eliminated based on Eq. ( 12), that is: ]p , [K ] , [ ? i -1 i ? = (18) The governing equation is a scalar one, but the source term K, i in Eq. ( 18) taken separately is vectorial. So, it is necessary also to put the interpolation function ? in the dyadic form. Among other options, one such a class of functions is given by: ip 3 p i j i , p R R RR 3 ? + = ? (19) In Eq. (18), R=R(X;X j ) is the Euclidian distance between the interpolation point X j and the field point X, ? ip is the Kronecker Delta operator, and: R p = [x p (X j )-x p (X)] (20) It is easily demonstrated that: p 3 j p R R = ? (21) Considering now the source term analyzed in Eq. ( 13) by the DIBEM approach, the complete governing matrix equation take the following form: b z ] A [ p ] H [ q ] G [ p ] H [ = = + ? (22) In Eq. ( 8), the lines of matrix A are comprised of vectors ? ?, which may be obtained from following the basic interpolation equation: ?][s] [ [F] ?][F] [ [F] ] [ ? 1 - ? 1 - ? = ? = ? (23) The radial basis function used in the DIBEM approach is the well know thin plate function, given by: ] R [ln R F # a) First example: one-dimensional flow Aiming to evaluate the robustness of the proposed model, in the following tests the function that characterizes the variation of the clearance h(x) between the bearing and shaft surfaces is changed. Actually, this distance is defined by the gap between two circular surfaces; here, it will be successively approximated by polynomials with crescent order, that is, by linear, quadratic and cubic functions. The clearance h(x), the geometry and the boundary conditions are shown in Figure 2. In these simulations, the pressure values are prescribed zero at the input and output, since internally the fluid flow imposes the overpressure values, while the null values of the flow along the x direction impose the mathematical one-dimensionality of the model. The average velocity V (see Eq. ( 1)) is assumed to be equal 1/12 and the density ? and viscosity ? are unitary. The difference between analytical and numerical values, divided by the highest analytical value, was chosen as a measure of errors. For the one dimensional cases, the analytical solution is available for comparison. Three meshes with 40, 80, 160 and 320 linear boundary elements with double nodes at the corners were used to simulate the pressure field and velocities, while the thin-plate radial function is employed to the DIBEM interpolation. Improving the approximation of the source term according required by the DIBEM approach, a different number of internal interpolation points also is used. The quantity of these points is indicated in each simulation. # i. Linear variation of the h(x) This simplest and hypothetical case appears as an example solved in most of the classic books dealing with the hydrodynamics of the bearings, serving to show how the height differential -in this case, a linear variation of clearance h(x) -implies a pressure value that throughout the domain examined. Figure 3 shows the profile obtained by the MEC for the two meshes, with 44 and 84 nodes, and respectively 49 and 81 interpolating internal points, in comparison with the corresponding analytical value. The results presented an excellent concordance, as can be observed. # Numerical Simulations For better detailing of the effect of internal interpolation points and boundary mesh refinement, a convergence curve is presented in Fig. 3, in which the average percentage error is expressed as a function of the increasing number of nodal points and interpolation points. It is verified that the values of the percentage errors in each mesh are very small and are mainly reduced with the boundary refinement. It can also be seen that only the introduction of many interpolating points without the proper refinement of the boundary is not very effective. It occurs because only the DIBEM procedure requires such points. ii. Quadratic variation of the h(x) Now the effect of height variation amplifies exponentially comparatively to the previous example so that a significant reduction in the precision of the numerical model employed is expected. Indeed, numerical errors have grown, but the results of this simulation continued with very good accuracy, as shown in Fig. 5. As shown in the previous case, Fig. 6 shows an error curve as a function of the boundary mesh refinement and insertion of internal interpolation points. # iii. Cubic variation of the h(x) Considering a cubic variation for h(x), the pressure gradients at the inlet are strongly accentuated and the less refined mesh already presents errors above 2%. It must be highlighted that the function that describes the height h(x) appears to the third power. However, the results are shown in Fig. 7 are still reasonable and reach a very satisfactory precision if more refined meshes are used, accompanied by a regular number of poles to represent the term proactive source or term of the governing equation, which is composed of the mean inlet velocity versus the derivative of the function h(x). # b) Second example: two-dimensional flow Many journal bearings have an oil feedback system, as well as many other functional upgrades that cannot be described here. The representation of a squeeze lubricant flow in the axial direction is done hypothetically through a source W(z) (see Eq. ( 4)). This source should be located in a restricted region but for simplicity, it will be assumed distributed throughout the domain. The purpose here is only to show that the model can solve suitably two-dimensional cases pertinent to the hydrodynamic theory of the journal bearings. In order to compare the results, a model generated from the Finite Element Method (FEM) (Reddy, 2005) using triangular elements with 20000 nodal points was taken as reference solution, since no analytical solution is available. The sources adopted have the following form: It must be highlighted that the insertion of the source W(z) not only alters the profile of velocities such as also change the values of pressure. For clarity, in Fig. 9 one three-dimensional view of the pressures on the domain is shown. ) z z ( 10 ) z ( W ; 1 V 2 ? = = (25) # Conclusions The boundary element model developed here was successfully implemented to study the hydrodynamic of the journal bearing problem, a case of great industrial interest. The rotating shaft and the slider bearing are given by two non-concentric circumferences whose clearance defines the lubricant fluid flow and the pressure field. Regarding the numerical model, this variable distance can be computed directly at the nodal level using the simplicity of the BEM discretization. Mathematically, this problem is expressed in terms of a non-homogeneous scalar field equation, composed of three terms with different physical meaning: the variable diffusivity, the advective effect, and the body force. The diffusive term has been well represented, although the cube of the function h(x) is approximated by linear boundary elements. The advective term was suitably approached by the Quasidual model through radial functions, as well as the body force term, related here to the source. The Quasi-dual solves accurately onedimensional cases in general or then two-dimensional cases that can be expressed by a potential function, which is the case of the hydrodynamic bearing. It does not require the internal inclusion of interpolating points, which are necessary only to the source representation. Regarding this term, in this case, the application of the DIBEM procedure did not require numerous poles, due to its relative mathematical simplicity. Unlike to the DRBEM, using the DIBEM approach the insertion of an excessive number of interpolating points does not produce disturbance effects in the numerical solution, commonly reported in the literature as due to ill conditioning matrix problems. Despite the necessary matrix inversion, the computational cost of this model is comparatively lower than that spent using alternative formulations such as DRBEM and RIM. The successful association between the two techniques based on the approach with radial basis functions opens new options for the BEM application, due to the similarity of the problem addressed and the modeling of other cases in which the properties of the constitutive medium vary gradually along the domain, common in geophysics and soil mechanics problems. 1![Figure 1: A schematic model for the journal bearing highlighting the meaning of the clearance h(x).](image-2.png "Figure 1 :") ![Boundary Element Model Applied to the Simulation of Journal Bearings interpolation points also are taken as source points(Loeffler and Mansur, 2017).](image-3.png "A") ![15)A Boundary Element Model Applied to the Simulation of Journal Bearings](image-4.png "") 2![Figure 2: Geometry of the domain with the applied boundary conditions.](image-5.png "Figure 2 :") 3![Figure 3: Numerical and analytical solutions of pressure along the fluid flow for a linear variation of h(x).](image-6.png "Figure 3 :") ![Boundary Element Model Applied to the Simulation of Journal Bearings](image-7.png "A") 4![Figure 4: Average percentage error as a function of mesh refinement and insertion of internal interpolating points for a linear variation of h(x).](image-8.png "Figure 4 :") 5![Figure 5: Numerical and analytical solutions of pressure along the fluid flow for quadratic variation of h(x).](image-9.png "Figure 5 :") 6![Figure 6: Average percentage error as a function of mesh refinement and pole insertion for quadratic variation of h(x).](image-10.png "Figure 6 :") © 2019 Global Journals j = (24) © 2019 Global Journals * The Boundary Element Method for Engineering CABrebbia 1978 Pentech Press London * CABrebbia JC FTelles LCWrobel Boundary Element Techniques Berlin Springer-Verlag 1984 * Boundary Element Techniques in Engineering CABrebbia SWalker 1980 Newnes-Butterworths; London * The radial integration method for evaluation of domain integrals with boundary-only discretization XWGao Engineering Analysis with Boundary Elements 26 2002 * MAGolberg CSChen The theory of radial basis functions applied to the BEM for inhomogeneous partial differential equations BE Communication 1994 5 * Quasi-Dual Reciprocity Boundary Element Method for Incompressible Flow: Application to the Diffusive-Advective Equation CFLoeffler WJMansur International Journal for Numerical Methods in Engineering 58 2003 * Direct Use of Radial Basis Interpolation Functions for Modelling Source Terms with the Boundary Element Method CFLoeffler ALCruz ABulcão Engineering Analysis with Boundary Elements 50 2015 * Solving Helmholtz problems with the boundary element method using direct radial basis function interpolation CFLoeffler HMBarcelos WJMansur ABulcão Engineering Analysis with Boundary Elements 61 2015 * A Regularization scheme applied to the direct interpolation boundary element technique with radial basis functions for solving the eigenvalue problem CFLoeffler WJMansur Engineering Analysis with Boundary Elements 74 2017 * Numerical Unsteady Analysis of Thin Film Lubricated Journal Bearing KMPanday PLChoudhury NPKumar International Journal of Engineering and Technology 2012 v.4, 2 * The Dual Reciprocity, Boundary Element Method, Computational Mechanics Publications and PWPartridge CABrebbia LCWrobel 1992 Elsevier London * Integral Equations and Boundary Value Problems MDRaisinghania . S. Chand Ed., New Dehli 2011 * An Introduction to the Finite Element Method JReddy 2005 McGraw-Hill 3rd Ed * Mechanical Engineering Design JEShigley CMischke RBudynas 2003 7rd Ed * Mcgraw-Hill