Analytical Solution of Reynolds Equation: Grubin’s Approximation

Analytical Solution of Reynolds Equation: Grubin's Approximation

Derivation of famous analytical solution of EHL equations is presented here. In 1949, Grubin obtained a solution for so called elasto-hydrodynamic lubrication (EHL) line contact problem with simplifications known as Grubin approximation. He was the first to combine both elastic deformation and lubricant hydrodynamic flow, including pressure-viscosity dependence. Although his solution did not satisfy strictly hydrodynamic equations of EHL, his analysis was recognized as very useful under highly loaded conditions, not only in stationary case, but also in transient conditions. Corresponding MATLAB code can be found <a href="http://www.tribonet.org/cmdownloads/semi-analytical-solution-of-1d-transient-reynolds-equationgrubins-approximation/">here</a>. <lite-youtube videoid="M7g7MecZoHU"></lite-youtube> Here, the solution approach is described for the case of stationary line contact (1D Reynolds equation). The base for the validity of Grubin's approximation is the fact that at highly loaded contacts, the hydrodynamic pressure is almost the same as Hertz contact pressure obtained for dry contacts (cylinder on plane). On top of it, the film thickness is almost constant within Hertz contact area, as shown in Fig. 1. [caption id="attachment_1884" align="aligncenter" width="536"]<img class=" wp-image-1884" src="http://www.tribonet.org/wp-content/uploads/2017/02/Hertz-Vs.-EHL-Pressure.jpg" alt="Hertz Vs. EHL Pressure" width="536" height="518" data-wp-pid="1884" /> Fig. 1. Hertz Vs. EHL[/caption] Following system of equations is considered: [math] \begin{eqnarray} \label{complete_sys1} \frac{d}{x} (\frac{h^3}{12\mu} \frac{dp}{dx}) = U_m\frac{dh}{dx} \\ \label{complete_sys2} \mu= \mu_0 e^{\alpha p} \\ \label{complete_sys3} p(-\infty) = 0 \end{eqnarray} [/math] where [math] h, p, \mu, \alpha [/math] are hydrodynamic film thickness, pressure, viscosity, and pressure-viscosity coefficient correspondingly. Integration of Reynolds equation gives: <p style="text-align: center">[math] \frac{h^3}{12\mu} \frac{dp}{dx} = h + C [/math]</p> where C is an integration constant. Since the pressure drops to zero in the infinity, there exist a point [math] x' [/math], so that [math] \frac{dp}{dx}=0 [/math]. Hence,[math] C=-h' [/math], where [math] h'=h(x') [/math]. Following expression can be obtained from the integrated Reynolds equation: [math] \begin{eqnarray} \label{complete_sys1} \frac{dp}{dx} = \frac{h-h'}{h^3} 12\mu U_m\\ \end{eqnarray} [/math] Further, a following variable is introduced: [math] p_{\mu} = \frac{1-e^{-\alpha p}}{\alpha} [/math]. Then, the equation is true: [math] \begin{eqnarray} \label{complete_sys1} \frac{dp_{\mu}}{dx} = e^{-\alpha p} \frac {dp}{dx} \\ \label{complete_sys2} =&gt; \frac{dp_{\mu}}{dx} = 12\mu_0 U_m  \frac{h-h'}{h^3} \\ \end{eqnarray} [/math] When pressure is very high, [math] \ p_{\mu} -&gt; \frac{1}{\alpha} \equiv const [/math] and [math] =&gt; \frac{dp_{\mu}}{dx} = 0 [/math]. From this it follows that [math] h'=h_0 [/math] and film is constant within the Hertz contact area. It is clear from Fig. 1 that outside the Hertz area, the hydrodynamic pressures are small compared to those within the contact and therefore it can be assumed that all the elastic defection is due to Herztian pressure. Therefore, the film thickness can be written in the form: [math] \begin{eqnarray} \label{complete_sys1} h = h_0 + h_{elastic} = h_0 + \frac{W_L \delta}{E_L} \\ \end{eqnarray} [/math] where [math] h_{elastic} [/math] is the combined elastic deflection of the surfaces (given by Hertz solution), [math] W_L [/math] is the applied load per unit length of the cylinder, [math] \frac{1}{E_L} = \frac{1}{\pi}(\frac{1-\nu_{1}^{2}}{E_1} + \frac{1-\nu_2^{2}}{E_2}) [/math] is the combined elastic modulus and [math] \delta [/math] is given by the following expression: [math] \begin{eqnarray} \label{complete_sys1} \delta = 2 \{\frac{x}{a} \sqrt{\frac{x^2}{a^2}-1} - ln{\frac{x}{a} + \sqrt{\frac{x^2}{a^2} - 1}} \} \\ \end{eqnarray} [/math] From the expression of the film thickness [math] h-h' = h-h_0 =  \frac{W_L \delta}{E_L}  [/math] and hence (equation 2): [math] \begin{eqnarray} \label{complete_sys1} =&gt; \frac{dp_{\mu}}{dx} = 12\mu_0 U_m  \frac{\frac{W_L \delta}{E_L}}{(h_0 + \frac{W_L \delta}{E_L})^3} \\ =&gt; \frac{dp_{\mu}}{dx} = \frac{12\mu_0 U_m}{(\frac{W_L }{E_L})^2}\frac{\delta}{H^3} \\ =&gt; \frac{(\frac{W_L }{E_L})^2}{12\mu_0 U_m}\frac{dp_{\mu}}{dx} = \frac{\delta}{H^3}\\ \end{eqnarray} [/math] By introducing dimensionless variables [math] p_{\mu}^{'} =p_{\mu} \frac{(\frac{W_L }{E_L})^2}{12\mu_0 U_m a} [/math] and [math] x'=\frac{x}{a} [/math], the equation reduces to the following form: [math] \begin{eqnarray} \label{complete_sys1} \frac{dp_{\mu}^{'}}{dx'} = \frac{\delta}{H^3} \end{eqnarray} [/math] Integration of this equation from [math] -\infty [/math] to [math] 1 [/math] gives [math] p'(-1) = \int_{-\infty}^{-1}{\frac{\delta}{H^3}} dx' [/math]. This equation was integrated by Grubin numerically for various values of [math] H_0 [/math] and further fitted with the following simple formula: [math] \begin{eqnarray} \label{complete_sys1} p'(-1) = 0.0986H_0^{-\frac{11}{8}} \\ \end{eqnarray} [/math] Now, for the line contact, [math] a=2\sqrt{\frac{W_L R}{E_L}} [/math] with [math] R [/math] being equivalent radius of the bodies. At [math] x'=-1 [/math] in dimensional coordinates [math] x=-a [/math]. Hence, [math] \begin{eqnarray} \label{complete_sys1} p_{\mu}(x=-a) = \frac{1}{\alpha} = \frac{12\mu_0 U_m a}{\frac{W_L }{E_L})^2} p'(-1) \\ \end{eqnarray} [/math] Substituting the expression for dimensionless pressure one can obtain following: [math] \begin{eqnarray} \label{complete_sys1} \frac{1}{\alpha} = \frac{12\mu_0 U_m a}{\frac{W_L }{E_L})^2} 0.0986(h_0 W_L/E_L)^{-\frac{11}{8}}\\ \end{eqnarray} [/math] And finally, the film thickness can be expressed as [math] \begin{eqnarray} \label{complete_sys1} h_0^{\frac{11}{8}} = 2.3664 \frac{\mu_0 U_m \alpha \sqrt{R}}{ (W_L/E_L)^{\frac{1}{8}}}\\ \end{eqnarray} [/math]