// Hydrodynamic Journal Bearing, iterative solution with Scilab
// Allowing for variation of oil viscosity with temperature
// 10th February 2015 - David J Grieve
printf(" \n");
printf("Analysis of a journal bearing with hydrodynamic lubrication, \n\n");
printf("To provide data about the variation of viscosity (Pa.s) with temperature, deg. C following is needed:,\n");
printf("Programme calculates viscosity in Pa.s using this formula: mU = (mu0*exp(b/(TC(i)+pluscon)))/1000.0 \n");
printf("Data is for straight SAE grades 10 30 or 60 \n");
printf("For grade 10: mu0=0.1975 b=468 pluscon=53 \n");
printf("For grade 30: mu0=0.1974 b=550 pluscon=53 \n");
printf("For grade 60: mu0=0.23375 b=660 pluscon=53 \n");// Set up frame for input data
// The sample data is a test example for comparison see related web page
labels=["Oil in Temp. deg. C";"mu0";"b";"pluscon";"Journal diam. mm";"Journal length, mm";"Radial clearance, mm";"Load, N";"Shaft RPM";..
"Number of intervals round journal";"Number of intervals along bearing length";..
"SoR factor to be used, less than 1";"Number of iterations";"Eccentricity ratio";"Oil density, kg/m3, about 850";..
"Oil specific heat, J/kg.deg.C, about 1800";"Oil supply pressure, Pa";"Oil supply hole diameter, mm"];
[ok,startTemp,mu0,b,pluscon,JDiam,JLength,radialClear,JLoad,shaftRPM,circInts,widthInts,sOR,iTerations,eccen,oilDen,oilSpecHeat,oilSupPress,oilHoleDia]=..
getvalue("Enter the fifteen items of data listed below",labels,list("vec",1,"vec",1,"vec",1,"vec",1,"vec",1,..
"vec",1,"vec",1,"vec",1,"vec",1,"vec",1,"vec",1,"vec",1,"vec",1,"vec",1,"vec",1,"vec",1,"vec",1,"vec",1),..
["60";"0.1974";"550";"53";"40";"40";"0.04";"2500";"1800";"180";"16";"0.9";"50";"0.59";"850";"1800";"100000";"4"]);
printf("Input data: \n\n");
printf("Oil in temperature deg C: %f\n",startTemp);
printf("mu0: %f\n",mu0);
printf("Constant b: %f\n",b);
printf("Constaaant pluscon: %f\n",pluscon);
printf("Journal diameter, mm: %f\n",JDiam);
printf("Journal Length, mm: %f\n",JLength);
printf("Radial clearance, mm: %f\n",radialClear);
printf("Load, N: %f\n",JLoad);
printf("Shaft RPM: %f\n",shaftRPM);
// printf("Lubricant viscosity, Pas: %f\n",lubeVisc);
printf("Number of intervals round full circumference of journal: %f\n",circInts);
printf("Number of intervals along length of journal: %f\n",widthInts);
printf("SoR factor to be used: (usually less than 1): %f\n",sOR);
printf("Number of iterations to be used: %f\n",iTerations);
printf("Eccentricity ratio: %f\n",eccen); //less than 1
printf("Oil density, kg/m3,(about 850): %f\n",oilDen);
printf("Specific heat of oil, J/kg.degC, (about 1800): %f\n",oilSpecHeat);
printf("Oil supply pressure, Pa: %f\n",oilSupPress);
printf("Oil hole diameter, mm: %f\n",oilHoleDia);
// Next - Start calculations;
// Calculate start viscosity from values input
mU=1:1:circInts;
mU(1)=(mu0*exp(b/(startTemp+pluscon)))/1000.0; // Using modified eqn from: A S Seireg and S Dandage, ref 1.
mU(2)=mU(1);
// Calculate Bearing characteristic number or Sommerfeld number;
somNum=(JDiam/(radialClear*2.0))*(JDiam/(radialClear*2.0))*mU(1)*shaftRPM*JDiam*JLength/(60.0*JLoad*1000000.0);
printf(" \n");
printf("Calculation results: \n\n");
printf("Oil viscosity at input temp. deg. C: %f\n",mU(1));
printf("Bearing characteristic number or Sommerfeld number: %f\n",somNum);
avProjPress=JLoad/(JDiam*JLength);
printf("Average projected pressure MPa: %f\n",avProjPress);
dx=JDiam*3.14159/(1000.0*circInts);
dx2=dx*dx;
dy=JLength/(1000.0*widthInts);
dy2=dy*dy;
dxdy=dx*dy;
dx2dy2=2*(dx*dx)+2*(dy*dy);
dTheta=2*3.14159/circInts; //dTheta is in radians
vel=3.14159*JDiam*shaftRPM/(60.0*1000.0);
// Next - Theoretical attitude angle
// theoAt=180.0*atan(sqrt(1-eccen*eccen)/eccen)/4.0;
// printf("Theoretical attitude angle, degrees: %f\n",theoAt);
// Next calc film thicknesses
x1=1;
for i=1:circInts,h(i:i,1:widthInts+1)=radialClear*(1+eccen*cos(i*dTheta))/1000.0, end;
//for i=1:circInts, Specifies a point at 'top' of an element in a column, round bearing
//for j=1:widthInts, Specifies point at the 'left' of an element in a row along the length of the bearing
shearForce=0;
for i=1:circInts,
fric=vel*mU(1)*JLength*JDiam*dTheta/(h(i:i,1:1)*2000000.0),
shearForce=shearForce+fric,
end;
// Use h values to determine the mass of oil in each longitudinal slice along journal
mSum=0;
for i=1:circInts-1,
m(i)=0.5*(h(i,1)+h(i+1,1))*JLength*dx*oilDen/1000.0, mSum=mSum+m(i),
end;
// Half m(1) is pulled in each slice advance, dx
mInTotPerRev=circInts*0.5*m(1); // kg
mInTotPerSec=mInTotPerRev*shaftRPM/60.0; // kg
volInTotPerSec=1000.0*mInTotPerSec/oilDen //litres
format(15);
disp(mInTotPerRev);
// disp(mInTotPerSec);
// disp(volInTotPerSec);
// Calculate force required to shear oil in each longitudinal slice and work done
// Also calculate temp rise caused by shear energy
startTemp=60.0;
pluscon=53; // factor for equation
temp(1)=startTemp; //in degrees C
TC=0; //1:1:circInts;
TC(1)=startTemp;
i=1;
// mu0=14*0.0141;
// b=550; // For SAE grade 30 oil
// mU=1:1:circInts;
// mU(1)=(mu0*exp(b/(TC(i)+pluscon)))/1000.0; // Using modified eqn from: A S Seireg and S Dandage, ref 1.
// mU(2)=mU(1);
f(1)=2*mU(1)*vel*dx*JLength/(1000.0*(h(1,1)+h(2,1)));
workDOS(1)=dx*f(1);
dTemp1(1)=workDOS(1)/(m(1)*oilSpecHeat);
tempSum=0.0;
for i=2:floor((circInts-1)/2),f(i)=2*mU(i)*vel*dx*JLength/(1000.0*(h(i,1)+h(i+1,1))),
workDOS(i)=dx*f(i), // Work done one slice per 1 slice advance
dTemp(i)=workDOS(i)/(m(i)*oilSpecHeat),
tempSum=tempSum+dTemp(i), // Temperature rise for slice i
// Add half temp rise to i+1 slice
TC(i)=TC(i-1)+((dTemp(i-1))/2.0)+dTemp(i),
// After each row recalculate viscosity
mU(i)=(mu0*exp(b/(TC(i)+pluscon)))/1000.0;
end;
// Do loop for half of journal circumference where h is increasing, assume constant mass
for i=ceil((circInts-1)/2):circInts-1,
f(i)=2*mU(i)*vel*dx*JLength/(1000.0*(h(i,1)+h(i+1,1))),
workDOS(i)=dx*f(i),
dTemp(i)=workDOS(i)/(m(ceil((circInts-1)/2))*oilSpecHeat),
TC(i)=TC(i-1)+((dTemp(i-1))/2.0)+dTemp(i),
end;
frictionTorque=shearForce*JDiam/2000.0;
fricPower=shearForce*vel; // Oil flow rate, 2 terms
oilQ1=vel*radialClear*JLength*eccen/1000000.0;
oilQ2=(1.2+11.0*oilHoleDia/JLength)*(h(1:1,1:1))*(h(1:1,1:1))*(h(1:1,1:1))*oilSupPress*oilHoleDia/(12.0*mU(1)*JLength);
oilQ=oilQ1+oilQ2;
oilTempRise=fricPower/((oilDen*oilQ)*oilSpecHeat);
printf("Friction force, N, is: %f\n",shearForce);
printf("Friction power, w, is: %f\n",fricPower);
printf("Oil flow, litres/s is: %f\n",oilQ*1000.0);
printf("Oil temperature rise, deg. C, is: %f\n",oilTempRise);
for i=1:circInts,
for j=1:widthInts+1, // zero matrices to store iterations. Start pressure distribution calcs
p1(i,j)=0.0;
p0(i,j)=0.0; // Calc coefficients cN, cS, cE, cW, cT - north, south, east, west and cT
end, end;
oilVisc=0.0;
for i=1:circInts, for j=1:widthInts+1,cN(i,j)=dx2/dx2dy2,cS(i,j)=cN(i,j), end,end;
for i=2:circInts-1, mU(i)=(mu0*exp(b/(TC(i)+pluscon)))/1000.0, oilVisc=mU(i),
for j=1:widthInts+1,hi3=h(i,j)*h(i,j)*h(i,j),p25=0.25*dy2/hi3,..
hi3p=h(i+1,j)*h(i+1,j)*h(i+1,j),..
hi3m=h(i-1,j)*h(i-1,j)*h(i-1,j),..
cE(i,j)=(dy2+p25*(hi3p-hi3m))/(dx2dy2),..
cW(i,j)=(dy2+p25*(-hi3p+hi3m))/(dx2dy2),..
end; end;
// Next solve matrix of eqns, k is iteration step no.
// To obtain a solution under relaxation must be used and thousands of iterations are normally needed
for k=1:iTerations,
for i=2:circInts-1,
f(i)=2*mU(i)*vel*dx*JLength/(1000.0*(h(i,1)+h(i+1,1))),
workDOS(i)=dx*f(i), // Work done one slice per 1 slice advance
dTemp(i)=workDOS(i)/(m(i)*oilSpecHeat),
tempSum=tempSum+dTemp(i), // Temperature rise for slice i
// Add half temp rise to i+1 slice
TC(i)=TC(i-1)+(dTemp(i-1)/2.0)+dTemp(i),
mU(i)=(mu0*exp(b/(TC(i)+pluscon)))/1000.0;
for j=2:widthInts,
hi3=h(i,j)*h(i,j)*h(i,j),
cT(i,j)=3.0*dx*dy2*vel*mU(i)*(h(i+1,j)-h(i-1,j))/(dx2dy2*hi3),
cT(i,1)=3.0*dx*dy2*vel*mU(i)*(h(i+1,j)-h(i-1,j))/(dx2dy2*hi3),
cT(i,widthInts+1)=3.0*dx*dy2*vel*mU(i)*(h(i+1,j)-h(i-1,j))/(dx2dy2*hi3),
p1(i,j)=(p0(i,j-1)+p0(i,j+1))*(cN(i,j))+p0(i+1,j)*cE(i,j)+p0(i-1,j)*cW(i,j)-cT(i,j), ///
//replace p0 with p1 using sor factortimes difference from above
if(p1(i,j))<0.0 p0(i,j)=0.0; end;
// SoR factor normally needs to be less than 1 - start with 0.9
pDiff=p0(i,j)-p1(i,j),
p0(i,j)=p0(i,j)-sOR*pDiff;
if(p1(i,j))<0.0 p1(i,j)=0.0; end;
end;end;end;
// Calculate element pressures as average of 4 corner node pressures
eSum=0.0,
for i=1:circInts,
for j=1:widthInts,
pE(i,j)=0.0; end, end;
for i=1:circInts-1,
for j=1:widthInts,
pE(i,j)=0.25*(p1(i,j)+p1(i+1,j)+p1(i,j+1)+p1(i+1,j+1));
if(pE(i,j))>0.0 eSum=eSum+pE(i,j)*dxdy;end;end;end;
printf("Sum of radial lubricant forces, (N) is: %f\n", eSum);
// Calculate element forces so load parallel and perpendicular to min. film thickness can be calculated
// Use negative sign
loadPerp=0.0; loadParallel=0.0;
for i=1:circInts-1,
for j=1:widthInts,
eLoad=pE(i,j)*dxdy;
loadParallel=loadParallel-eLoad*cos(i*dTheta);
loadPerp=loadPerp+eLoad*sin(i*dTheta);
end; end;
printf("Load parallel to minimum film thickness, (N) is: %f\n", loadParallel);
printf("Load perpendicular to minimum film thickness, (N) is: %f\n", loadPerp);
resLoad=sqrt((loadPerp*loadPerp)+(loadParallel*loadParallel));
printf("Resolved load, (N) is: %f\n", resLoad);
// Calculated angle of attitude is in degrees after the downwards vertical - which is the direction of the load
calcAt=180.0*atan(loadPerp/loadParallel)/3.14159;
calcAtRad=atan(loadPerp/loadParallel);
calcAtFloor=floor(calcAt);
printf("Calculated angle of attitude (degrees) is: %f\n", calcAt);
// Calculate max and av pressures
maxPress=0.0; maxPresPosni=0; maxPresPosnj=0; ;
for i=2:circInts-1,
for j=1:widthInts+1,
if(p1(i,j))>maxPress, maxPress=p1(i,j), maxPresPosni=i, maxPresPosnj=j; end;
end;end;
maxPressMPa=maxPress/1000000.0;
printf("Maximum pressure is: (MPa) %f\n", maxPressMPa );
ratPAvPMax=avProjPress/maxPressMPa;
printf("Ratio average pressure/max. pressure (load/projected area*max. pressure) is: %f\n", ratPAvPMax);
// Calculate angular position of max pressure past downward vertical
maxPressAng=calcAt-180.0+maxPresPosni*dTheta*180.0/3.14159;
printf("Angle of maximum pressure past downwards vertical %f\n\n", maxPressAng );
printf("For diagnostics, type variable name at prompt \n");
printf("eg: for list of pressures on each element at prompt type pE then RTN \n");
// Prepare plot of pressure (vert axis) round cente line (straight, horiz axis)
cT1=0;
cI1=round(360/circInts); //value in degrees of 1 circular interval
xx=[round(2*cT1/circInts):round(cT1/circInts):round(cT1*(circInts-1)/circInts)];
//xx=[round(2*%pi/circInts):round(%pi/circInts):round((circInts-1)*%pi/circInts)];
// xxxx=[4*%pi/180:2*%pi/180:358*%pi/180];
for i=2:circInts-2, yyA(1,i)=pE(i,maxPresPosnj)/1000000.0; end;
figure(2);
plot2d(xx,yyA);
title('Plot of element pressures in MPa along centre line of journal - radians - starting at maximum film thickness','fontsize',3);
xtitle('For angular position past max. oil film thickness multiply by 360/number of circular intervals');
// Polar plot of base circle
// xxbase=[round(2*%pi/circInts):round(%pi/circInts):round((circInts-1)*%pi/circInts)];
xxbase=[4*%pi/180:2*%pi/180:358*%pi/180];
// xxbase=[4:2:358];
yybase=[3.0004:0.0002:3.0358];
// Polar plot of pressure - first determine centre element width number
// Shaft rotates in clockwise direction
widthIntsMax=0;
if(widthInts-fix(widthInts/2)*2==0),widthIntsMax=widthInts/2,else widthIntsMax=ceil(widthInts/2),end;
// disp(widthIntsMax)
figure(3);
for i=2:circInts-2, yyB(1,i)=3+pE(i,widthIntsMax)/1000000.0; end;
ii=%pi*(90-calcAt)/180;
for i=2:circInts-2,xxx(i)=ii-2*%pi/circInts,ii=ii-2*%pi/circInts,end;
// for i=2:circInts-2,xxx(i)=ii-2*%pi/180,ii=ii-2*%pi/180,end;
polarplot(xxbase,yybase),polarplot(xxx,yyB),;
title('Polar plot of pressure in MPa on centreline of journal. Zero pressure is on 3. Load on shaft is down','fontsize',3);
figure(4);
x=[1:widthInts];
y=[1:circInts];
z=pE(y,x);
contour(y,x,z,12);
title('Contour plot of oil pressures on elements MPa','fontsize',3);
xtitle('For angular position past max. oil film thickness multiply by 360/number of circular intervals');
// Plot temperatures round journal
figure(5);
for i=2:circInts-2, yyATemp(i)=TC(i),end;
plot2d(xx,yyATemp);
title('Plot of oil temperature round journal journal - degrees C - starting at maximum film thickness','fontsize',3);
// xtitle('For angular position past max. oil film thickness multiply by 360/number of circular intervals');
mU(circInts)=mU(circInts-1);
figure(6);
plot2d(y,mU);
title('Plot of oil Viscosity, starting at max. oil film thickness','fontsize',3);