Investigation of Couple Stress Fluid and Surface Roughness Effects in the Elastohydrodynamic Lubrication Problems using Wavelet-Based Decoupled Method

Abstract: The standard decoupled method with Newton-generalized minimum residual procedure performs poorly or may break down when used to solve certain elastohydrodynamic lubrication problems. The authors Zargari et al., 2007 presented decoupled and coupled methods in which the limitations of the decoupled method for some set of physical parameters and slight variation in these values (physical parameters) give a non-convergence solution. To overcome this aspect the wavelet-based preconditioners technique is used in this paper to solve the elastohydrodynamic lubrication problem. The effect of coupled stress fluids on elastohydrodynamic lubrication behavior in smooth, as well as rough, contact is investigated using the proposed method, numerically, in a low-speed-high load condition. The elastohydrodynamic lubrication characteristics computed for couple stress fluids are found to have strong dependence on the couple stress parameter, which shows the versatility and applicability of the proposed method.


Introduction
The elastohydrodynamic lubrication (EHL) problem is modeled by two main groups of equations, the first concerned with the physical model of the lubricant and the second one governing the lubrication problem, itself.The lubricant model specifies the dependence of the fluid viscosity pηq and density pρq on the pressur pPq.These are empirical relations and are highly nonlinear, since most lubricants are non-Newtonian fluids [1,2].The governing equations consist of the following: the Reynolds equation, the film thickness equation, and the force balance equation.The system of equations is derived from a dimensionless, thin-film approximation of the Stokes equation, coupled with a linear elastic model of the contacting surfaces.The inlet boundary and the free cavitation point boundary define the region where the Reynolds equation is valid.The main objective of the present work is to detect the location of the cavitation point using the proposed wavelet-based decoupled method.Subsequently, the present study addresses the compressible EHL problem involving the effect of couple stress fluid in smooth as well as rough line contacts.The constitutive relationship for couple stress fluids is presented by Stokes [3] in the analysis to derive the modified Reynolds equation in terms of a dimensionless couple stress parameter which represents the molecular length of the additives.The surface roughness is assumed to be single-sided, transverse, and sinusoidal [4][5][6][7].
The ill-conditioned matrices in the system of non-linear algebraic equations arising in EHL confine the application of even more powerful schemes of the Krylov subspace method.The suitable remedy is the preconditioning for such matrices.Classical preconditioners are incomplete LU decomposition (ILU), the approximate inverse preconditioner, etc. [8][9][10].However, there are large classes of matrices occurring in the modeling of problems which are not amenable to these latest non-stationary iterative schemes with classical preconditioners for their solutions.Zargari et al. [11] used versatile MatLab code from KINSOL in analyzing EHL line contact problems.They observed that the decoupled scheme does not give a convergent solution for some set of physical parameters.Wavelet-based preconditioners provide some remedy in such challenging cases.Chen [12] highlighted the construction and implementation of various wavelet-based preconditioners for the solution of large classes of ill-conditioned systems very effectively.In EHL, an appeal to such schemes has resulted in the construction of suitable wavelet-based preconditioners [13].The dense and diagonal singular structure of the Jacobian demands more sophisticated iteration schemes from Krylov subspaces.The linear system is A x " b, if A is a symmetric and positive definite then the conjugate gradient (CG) method works well, whereas for non-symmetric and other systems, the generalized minimum residual method (GMRES) can be used [14][15][16].As the Jacobian matrix is full in EHL problems, the convergence of the iteration is not guaranteed and takes a large number of iterations to converge.It is essential to find effective preconditioners.The Jacobian matrices in EHL are dense with a non-smooth diagonal and smooth away from the diagonal band.This smoothness of the matrix transforms into smallness in the wavelet transform and facilitates in the design or construction of efficient preconditioners using discrete wavelet transform (DWT).DWT with permutation (DWTPer) are designed and implemented by [12,13].In this paper, we developed the wavelet-based decoupled method for the investigation of surface roughness effects in EHL line contact problems using a couple stress fluid.
This paper is arranged as follows: Section 2 contains the discrete wavelet transform.The wavelet-based decoupled method of the solution is discussed in Section 3. Section 4 presents the method of implementation and numerical findings.In Section 5, results and discussions are given.Finally, the conclusion of the proposed work is discussed in Section 6.

Discrete Wavelet Transform (DWT)
A major problem in the growth of wavelets during the 1980s was the search for a multiresolution analysis where the scaling function was compactly supported and continuous.As we already know, the Haar multiresolution analysis is generated by a compactly-supported scaling function, but it is not continuous.The B-splines are continuous and compactly supported but fail to form an orthonormal basis.A family of multiresolution analyses was generated by scaling functions, which are both compactly supported and continuous.These multiresolution analyses were first constructed by Daubechies [17,18], which created great eagerness among mathematicians and scientists conducting research in the area of wavelets.

Multiresolution Analysis
Wavelets are functions generated from one single function called the mother wavelet by the simple operations of dilation and translation.A mother wavelet gives rise to a decomposition of the Hilbert space L 2 pRq, into a direct sum of closed subspaces W j , j P Z [17].
Let ψ j, k pxq " 2 j{2 ψp2 j x ´kq and: Then every f P L 2 pRq has a unique decomposition: where s j P W j for all j P Z, it is: Using this decomposition of L 2 pRq, a nested sequence of closed subspaces V j , j P Z of L 2 pRq can be obtained, defined by: These closed subspaces V j , j P Z ( of L 2 pRq, form a "multiresolution analysis" with the following properties: Let φ P V 0 , the so-called "scaling function" that generates the multiresolution analysis V j ( jPZ of L 2 pRq.Then: is a basis of V 0 , and by setting: it follows that, for each j P Z, the family: is also a basis of V j .Then, since φ P V 0 is in V 1 and since φ 1, k : k P Z ( is a basis of V 1 , there exists a unique sequence a k that describes the following "two-scale relation": of the scaling function φ [19].

Daubechies's Wavelets
Different choices for φ may yield different multiresolution analyses and the most useful scaling functions are those that have compact support.As an example of multiresolution analysis, a family of orthogonal Daubechies wavelets with compact support has been constructed by Daubechies [18].
A wavelet basis is orthonormal if any two translated or dilated wavelets satisfy the condition: where δ is the Kronecker Delta function.
Each wavelet family is governed by a set of 2 J (an even integer) coefficients, a k : k " 0, 1, . . ., 2 J ´1 through the two-scale relation: Based on the scaling function φ 2 J pxq, the mother wavelet can be written as: Since the wavelets are orthonormal to the scaling basis, the coefficients of the scaling function and the mother wavelet for the two-scale equation are related by: Daubechies [17] found and exploited the link between vanishing moments of the wavelet ψ and regularity of the wavelet and scaling functions, ψ and φ.The wavelet function ψ has K vanishing moments if: and a necessary and sufficient condition for this to hold is that the integer translates of the scaling function φ exactly interpolate polynomials of degree up to K.That is, for each k, there exists constants c l such that: Daubechies introduced scaling functions satisfying this property and distinguished by having the shortest possible support.The scaling function φ 2 J has suppor " 0, 2 J ´1‰ , while the corresponding wavelet function ψ 2 J has support in the interval " 1 ´2J {2, 2 J {2 ‰ and has p2 J {2 ´1q vanishing wavelet moments.Thus, according to Equation (8) Daubechies scaling functions of order 2 J can exactly represent any polynomial of order up to, but not greater than 2 J {2 ´1 [18].For example, Daubechies family of wavelets when N " 4, we have low-pass filter coefficients [19]

Discrete Wavelet Transform Matrix
A DWT matrix is a linear transformation that transforms vectors from the standard basis to a wavelet basis.Certain classes of linear operators that correspond to dense matrices in the standard basis can be approximated by sparse matrices in a suitably-chosen wavelet basis.We use this property to design preconditioners for linear systems based on representing the system in a wavelet basis and forming a sparse approximation to the transformed matrix [20].The process of applying a DWT matrix to a dense matrix in order to obtain a sparse approximation is known as wavelet compression.The DWT of a n ˆn matrix, we form the orthogonal matrices pW 1 , W 2 , . . ., W k q are as follows:

Discrete Wavelet Transform with Permutation (DWTPer) Matrix
The standard DWT is not adapted to compressing matrices with non-smooth (NS) diagonal bands for preconditioning purposes because the pattern of non-zero elements in the compressed matrix tends to lead to a large amount of fill-in during LU factorization.One way of avoiding this problem is to use NS-forms, which we consider briefly here.We then present the non-standard DWT with permutation (DWTPer) and look in detail at its effect on dense matrices with non-smooth features both on and off the diagonal.Ford [21] proposed the DWTPer as an alternative way of avoiding the creation of finger pattern matrices.This is implemented by performing a standard DWT followed by a permutation of the rows and columns of the matrix to center the fingers about the leading diagonal.The DWTPer matrices pW 1 , W 2 , . . ., W k q are as follows: Here, each I is an identity matrix of dimension 2 J´1 ´1 and the Φ's are block zero matrices of the appropriate sizes.For J = 1, both I and Φ are of dimension 0.

Wavelet-Based Decoupled Method for the Numerical Solution of EHL Problem
The multigrid method is commonly used to solve the Reynolds equation efficiently [22,23].Since EHL problems are nonlinear, the wavelet-based decoupled method is appropriate.(i) We begin with an initial guess for P, H 0 , and the cavitation point X c ; (ii) evaluate H from the film thickness equation; (iii) solve the Reynolds equation for P; (iv) update H 0 using the force balance equation; (v) move the cavitation point based upon the value of dP dX at the cavitation point; (vi) while not converged, go to (ii); (vii) the overall method is such that the cavitation point is only updated on the finest grid and H 0 is only updated on the coarsest level so as to ensure smooth convergence of the whole method; and (viii) numerical under-relaxation parameters C 1 and C 2 are used for pressure and film thickness equations, respectively, at (iii) and (iv).
A detailed method of solution of the proposed method is as follows.
To initiate the wavelet-based decoupled method for the solution of EHL problem with line contacts using couple stress fluid and surface roughness effects, an initial guess is made for pressure pPq distribution and the offset film thickness H 0 [24].These values use the calculation of film thicknes pHq, density pρq [25], and viscosity pηq [26], respectively.The Reynolds equation is discretized over a uniform mesh of different mesh sizes using finite difference approximation and solved along with the film thickness equation and force balance equation using the wavelet-based decoupled method as follows, Consider the Reynolds equation of the form: with: PpX in q " PpX out q " dPpX c q dX " 0 Taking the initial values of Hertzian pressure Ppxq " , then evaluates H from the film thickness equation.
Finite difference discretization of Equation (15) gives the system of nonlinear algebraic equations: where n " 2 J , n is the number of grid points and J is the level of resolution.
The nonlinear systems of Equation ( 17) can be solved using Newton's method as follows: where C 1 " 0.4 p0 ă C 1 ă 1q is the under-relaxation parameter: Let, Equation ( 19) reduces to the linear system of equations: where B " rF 1 pP n qs nˆn , x " rS n s nˆ1 , f " rFpP n qs nˆ1 .Applying the DWT-based decoupled method (DWT-DM) or DWTPer-based decoupled method (DWTPer-DM) up to Jth level to Equation (20) as follows: First level: Equation ( 20) reduced to: Second level: Equation ( 21) becomes: Jth level: At the coarsest level we have: Consider the transformed linear system at the Jth level (coarsest level) in Equation ( 23), from which we can write B J as B J " D J `CJ [7] based on operator (matrix) splitting.Then, we need to choose an efficient preconditioner matrix M " D J .Multiplying M ´1 on both sides of Equation ( 23), we get: Solving Equation ( 24) iteratively, we obtain x J .Now applying successively the inverse discrete wavelet transform-based decoupled method (IDWT-DM) or the inverse discrete wavelet transform with permutation-based decoupled method (IDWTPer-DM) to x J at the level J, which reduces to ‰ nˆn rx 1 s nˆ1 .Substitute x p" S n q in Equation ( 18) we get the required pressure P.
Fix the cavitation point X c using dP dX " 0 at X " X c , which we can write as P i ´Pi`1 ∆x " 0 at X " X c .If the pressure gradient is zero at the initial guess X c is fine.Otherwise, move X c one interval to the left or right depending on whether ˆPi ´Pi`1 ∆x ˙Xc ą 0 or ˆPi ´Pi`1 ∆x ˙Xc ă 0, respectively.
Approximate a new cavitation point X c1 and solve the EHL problem.If the pressure gradient is zero at X c1 then we get a converged solution; otherwise check for the condition whether: We obtain the residual (Res) from the following force balance Equation ( 26): and then update H 0 in the following Equation ( 26): where C 2 " 0.03 p0 ă C 2 ă 1q is the under-relaxation parameter for pressure.
If condition Equation ( 25) is not met, repeat the procedure.

Numerical Experiment
Consider the governing modified Reynolds equation [4,27] in its dimensionless form: where, k " 3Uπ 2 4W 2 , ξ " and Ω " 8W H π Lm with the boundary conditions, inlet boundary P " 0 at X " X in , and outlet boundary P " BP BX " 0 at X " X o u t .The film thickness equation incorporating the sinusoidal roughness term in dimensionless form is given as: where ν " νR b 2 is the dimensionless surface displacement given by: The dimensionless force balance equation is given by: The dimensionless form of viscosity ηpPq [25] and density ρpPq [26], are given by: η " exp " plogη 0 `9.67q ˆ ´1 `p1 `5.1E ´09 P p h q z (‰ (32) and: ρ " ˆ1 `0.6E ´09 P p h 1 `1.7E ´09 P p h
The spatial domain X P rX in , X o u t s is discretized with a uniform grid of n points X i p1 ď i ď nq.
Consider the cavitation point X c to be located at an unknown internal poin X j p2 ď j ď n ´1q.
The Equations ( 28), (29), and (31) are discretized using second-order finite differences with uniform grid points x i , 1 ď i ď n, and the domain of interest is rX in , X o u t s " r´2, 1.5s, X c is the cavitation point to be determined in the solution process.The discretized form of modified Reynolds equation [27]: The film thickness equation approximated at x i on the regular grid is given by: where: for i " 0, 1, 2, . . ., n and j " 0, 1, 2, . . ., n, and the force balance equation in discrete from is: Use boundary conditions in Equations ( 34), (35), and (37) by taking the initial values of P, X c , and H 0 , then solving Equation (35), we get H. Next, we take Equation (34) in the form of FpP i q " 0 to solve this system using procedure explained in Section 3. The numerical results are presented in figures and tables, which are shown in next section.

Results and Discussion
The EHL characteristics are presented for one set of physical parameters, i.e., low speed-high load pU " 2E ´11, W " 4E ´04q, maximum Hertzian pressure pp h " 5.8E `08q.The results are obtained using the values of input parameters given in Table 1.The method in [11] is not converged for some set of physical parameters, whereas the proposed method overcomes the above limitations, which shown in Table 2 and Figures 1-3 for different grid points.Figure 4a,b presents the numerical solutions of pressure profile and film shapes respectively, using wavelet-based decoupled method for pure base oil pLm " 0q with three values of couple stress parameter Lm " 1E ´05, 2E ´05, 5E ´05 for a low speed-high load pU " 2E ´11, W " 4E ´04q condition.It is clear from Figure 4a,b that, couple stress fluids cause negligible change in the pressure profile and film thickness, respectively.Figure 5a,b represents the detailed solution presented in Figure 4a,b.On closer examination, it is observed that pressure increases marginally with an increase in Lm from 0 to 1E ´05, whereas a marginal decrease in pressure as Lm is increased further to 5E ´05.Therefore, it implies that the initial rise in pressure due to increased lubricant viscosity is compensated by the pressure reduction caused by fluid film thickening.Figure 6a,b compares the pressure and film thickness solutions, respectively, for the case of rough surface pA " 0.2, λ " 0.12q obtained using pure base oil pLm " 0q with couple stress fluid at Lm " 1E ´05, 2E ´05 and 5E ´05 for low speed-high load pU " 2E ´11, W " 4E ´04q condition.It can be seen from Figure 6a that the local pressure peaks at asperity tips for the case of couple stress fluids are marginally lower than those for Newtonian fluid.This is attributed to a reduction in the flow resistance offered by surface roughness due to fluid film thickening, as apparent from Figure 6b, in the case of couple stress fluids.Similarly, Figure 7a,b shows the detailed observations presented in Figure 6a,b.

Conclusions
The new decoupled method presented in this paper, called the wavelet-based decoupled method has been shown to be effective and versatile.The standard and other decoupled methods fail to converge for EHL problems discussed in [11].The wavelet-based decoupled method does ensure such convergence.The properties of wavelets should permit the wavelet-based decoupled

Conclusions
The new decoupled method presented in this paper, called the wavelet-based decoupled method has been shown to be effective and versatile.The standard and other decoupled methods fail to converge for EHL problems discussed in [11].The wavelet-based decoupled method does ensure such convergence.The properties of wavelets should permit the wavelet-based decoupled method to be efficiently applied in many cases through use of compression.The results shown in this paper have demonstrated that, it is worthwhile to further explore the usefulness of this new decoupled method.The pressure distribution solution varies only marginally when couple stress fluid is employed instead of a Newtonian fluid at low speed-high load.The maximum value of the couple stress parameter used in the analysis.Additionally, the reduction in minimum film thickness due to surface roughness may be compensated by film thickening caused by couple stress fluid and the desired film thickness is obtained by selecting an appropriate value of the couple stress parameter.

Figure 7 .
Figure 7. Detailed numerical solution of the EHL problem with couple stress and surface roughness around the Petrusevich spike using wavelet based decoupled method pn " 256, U " 2E ´11, W " 4E ´05 and G " 5E `03q.(a) Pressure solution; and (b) film thickness solution.

Table 2 .
The sensitivity of convergence of three schemes to the physical parameters with couple stress pLm " 1E ´05q and surface roughness effects pA " 0.2, λ " 0.12q.