// Hydrodynamic Journal Bearing, iterative solution with Scilab
// 28th January 2015 - David J Grieve
printf(" \n");
printf("Analysis of a journal bearing with hydrodynamic lubrication, \n\n");
// Set up frame for input data
// The sample data is a test example for comparison see related web page
labels=["Journal diam. mm";"Journal length, mm";"Radial clearance, mm";"Load, N";"Shaft RPM";..
"Lubricant viscosity, Pas";"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,JDiam,JLength,radialClear,JLoad,shaftRPM,lubeVisc,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),..
["40";"40";"0.04";"2500";"1800";"0.0275";"180";"16";"0.9";"50";"0.59";"850";"1800";"100000";"4"]);
printf("Input data: \n\n");
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 Bearing characteristic number or Sommerfeld number;
somNum=(JDiam/(radialClear*2.0))*(JDiam/(radialClear*2.0))*lubeVisc*shaftRPM*JDiam*JLength/(60.0*JLoad*1000000.0);
printf(" \n");
printf("Calculation results: \n\n");
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*lubeVisc*JLength*JDiam*dTheta/(h(i:i,1:1)*2000000.0),
shearForce=shearForce+fric,
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*lubeVisc*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;
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, 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),..
cT(i,j)=3.0*dx*dy2*vel*lubeVisc*(h(i+1,j)-h(i-1,j))/(dx2dy2*hi3),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=3:circInts-2,
for j=2:widthInts,
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);end;end;
//replace p0 with p1 using sor factortimes difference from above
for i=3:circInts-2,
for j=2:widthInts,
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');