Sinusoidal bearing in 1D: analytical solution

Analytical solution of Reynolds equation with mass-conserving cavitation model (JFO)

Supplementary material to the paper:
	J. Lengiewicz, M. Wichrowski and S. Stupkiewicz
	Mixed formulation and finite element treatment of the mass-conserving cavitation model
	Tribology International, vol. 72, pp. 143-155 (2014)
	http://dx.doi.org/10.1016/j.triboint.2013.12.012

Corresponding author:
	Stanislaw Stupkiewicz
	Institute of Fundamental Technological Research (IPPT PAN), Warsaw, Poland
	sstupkie@ippt.pan.pl
	http://www.ippt.pan.pl/~sstupkie

The file can be downloaded from:
	http://www.ippt.pan.pl/~sstupkie/files/bearing.html

version 1.3 (2014-02-01)

Brief description of the problem

This Mathematica notebook presents an analytical solution of a 1D hydrodynamic lubrication problem with cavitation for a sine-shaped lubricant film thickness. The mass-conserving  JFO cavitation model is used. The problem is defined in our paper, where all the details can be found.

The solution consists of two full-film regions separated by a cavitation region. The one-dimensional steady-state Reynolds equation holds in the full-film regions
(1)	(d/(dx))(u h-((h^3)/(12η))((d p)/(dx)))=0
where
(2)	h(x)=(1/2)((Subscript[h, 1])+(Subscript[h, 2]))-(1/2)((Subscript[h, 1])-(Subscript[h, 2]))cos(((π x)/L)),		-L≤x≤L
Equation (1) is integrated to yield
(3)	((d p)/(dx))=12η u ((h-(SuperStar[h]) )/(h^3))
where (SuperStar[h]) is an unknown integration constant.

The solution in the cavitated region is governed by the following mass-balance equation
(4)	(d/(dx))((1-λ)u h)=0
which after integration gives an algebraic equation for the unknown void fraction λ
(5)	(1-λ)u h=u (SuperStar[h])
Here, the unknown integration constant (SuperStar[h]) is equal to that in Eq. (3) so that the mass flux is conserved.

Equation (3) is integrated symbolically using Mathematica. The position of the film rupture boundary is then obtained numerically by requiring that the pressure and its derivative are equal to zero at this boundary. Finally, the position of the reformation boundary is obtained numerically by requiring that the pressure is equal to zero at this boundary and that the mass flux is preserved in the full film region. The void fraction λ is obtained from the mass-conservation equation (5) in the
cavitated region.

Note that direct integration of the Reynolds equation leads to a non-physical discontinuous solution due to the trigonometric functions involved. For this reason special parameterization by p0 (pressure at x=0) and (SuperStar[h]) (integration constant) is adopted, and appropriate limits for x->-L and x->L are exploited.

The notebook comes as it is. It gives the analytical solution for the input parameters used in the original paper. The solution procedure has not been tested for other input parameters. But it should work also for a modified input, possibly after the initial guess for the FindRoot[ ] is also modified.

Feel free to use and modify this notebook for your work provided the source, i.e. our paper, is adequately cited.

Symbolic integration of Reynolds equation

Clear variables

Clear[h,h1,h2,L,ηU,x,p]

Set the numerical values of the input parameters of the problem

Input parameters:
	h1p - maximum film thickness
	h2p - minimum film thickness
	Lp - half-length of the domain
	ηUp - product of viscosity η and entrainment speed U
	pBC - prescribed pressure at x=-L and x=L (boundary conditions)
[Attention! There is a misprint in the paper: the total length of the bearing should be l=125 mm and not l=12.5 mm.]

h1p=0.025;
h2p=0.015;
Lp=125./2;
ηUp=0.015*(10^(-6))*(4000/2);
pBC=1;

Set the numerical values for the initial guess for the root finding function

Rupture boundary:

xcini=25.;
p0cini=3.;

Reformation boundary:

xrini=55.;
p0rini=-2.;

Define film thickness h(x) for -L≤x≤L

h=(h1+h2)/2-(h1-h2)/2Cos[ π x/L ]

Block[ {h1=h1p,h2=h2p,L=Lp},Plot[ h,{x,-L,L},PlotRange->{0,h1},AxesLabel->{"x","h(x)"} ] ]

Integrate Reynolds equation in the full-film region

p0 - unknown pressure at x=0
hs - unknown integration constant
h1, h2, L, ηU - parameters of the problem

pSolution[x_,h1_,h2_,L_,ηU_,hs_,p0_]=p[x]/.DSolve[ {p'[x]==12ηU/h^3( h-hs),p[0]==p0},p[x],x ][[1]]//Simplify

Determine pressure at x=-L and x=L as a function of p0 and hs

p1=Limit[ pSolution[x,h1,h2,L,ηU,hs,p0],x->-L,Direction->-1 ]//Simplify[ #,L>0&&h1>0&&h2>0&&μU>0 ]&

p2=Limit[ pSolution[x,h1,h2,L,ηU,hs,p0],x->L,Direction->1 ]//Simplify[ #,L>0&&h1>0&&h2>0&&μU>0 ]&

Compute hs corresponding to p1=pBC (denoted by hs1) and to p2=pBC (denoted by hs2)

hs1=hs/.Solve[ p1==pBC,hs ][[1]]

hs2=hs/.Solve[ p2==pBC,hs ][[1]]

Pressure gradient

DpSolution[x_,h1_,h2_,L_,ηU_,hs_,p0_]=D[pSolution[x,h1,h2,L,ηU,hs,p0],x]//Simplify

Find numerically the position xc of rupture boundary such that the pressure and its gradient at xc are equal to 0

The corresponding pressure p0 is also computed, and subsequently the integration constant hs.

Block[ {h1=h1p,h2=h2p,L=Lp,ηU=ηUp},
{xc,p0c}={x,p0}/.FindRoot[
{pSolution[x,h1,h2,L,ηU,hs1,p0],DpSolution[x,h1,h2,L,ηU,hs1,p0]},
{{x,xcini,xcini+1},{p0,p0cini,p0cini+0.1}}
];
hsc=hs1/.p0->p0c;
{xc,p0c,hsc}
]

Find numerically the position xr of the reformation boundary such that the pressure is equal to 0

The second equation (hs2=hsc) implies that the flux is identical in both full-film regions.

Block[ {h1=h1p,h2=h2p,L=Lp,ηU=ηUp},
{xr,p0r}={x,p0}/.FindRoot[
{pSolution[x,h1,h2,L,ηU,hs2,p0],hs2-hsc},
{{x,xrini,xrini+1},{p0,p0rini,p0rini+0.1}}
];
{xr,p0r}
]

Construct the solution

The void fraction λ in the cavitation zone is computed from the mass balance equation: (1-λ)h U=(SuperStar[h]) U

Block[ {h1=h1p,h2=h2p,L=Lp,ηU=ηUp},
pAnalytical[x_]=Piecewise[
{{pSolution[x,h1,h2,L,ηU,hsc,p0c],x≤xc},
{0,xc<x≤xr},
{pSolution[x,h1,h2,L,ηU,hsc,p0r],xr<x}}];
λAnalytical[x_]=Piecewise[
{{0.,x≤xc},
{λ/.Solve[(1-λ)h==hsc,λ][[1]],xc<x≤xr},
{0.,xr<x}}];
]

Plot the solution

Plot[pAnalytical[x],{x,-Lp,Lp},Axes->False,Frame->True,FrameLabel->{"x","p"},PlotRange->{-0.4,7.4}]
Plot[λAnalytical[x],{x,-Lp,Lp},Axes->False,Frame->True,FrameLabel->{"x","λ"},PlotRange->{-0.02,0.32}]