Bearing model - viscosity/clearance loop?

Hi All,

I was wondering if there is someone out there who can potentially help out. After many sleepless nights I have managed to get a rough grasp of matlab code to be able to produce the below, and although it runs, I believe that it's not running quite correctly as the temperature change is nowhere near as high as expected.

This code intends to simulate the temperature increase in a journal bearing given a rapid acceleration.

Please take a look at the code below. I believe the issue may lie because the viscosity of the lubricant is not taken into account with every loop. The viscosity is a function of clearance. Unfortunately I am unsure how to build this into the code. So as the clearance 'C' changes, the viscosity changes also.

Will I need to set up another loop in order to achieve this?

%%New Viscosity = Mu_0*(C/Co) 
% How do I build this into an iterative loop? The condiditons to break would be if C>Co

Sorry if this is vague, but I have been looking at this for so long now I cant see the wood for the trees.

And again apologies for what is probably an infantile question.

Thanks.

clc
clear all
close all

%%=========parameters===========

n=3600;  %Rotational speed of the shaft in rpm 
Co=0.000125; %initial radial clearance (m)
T_supply = 40+273; % Oil Supply temperature %initial temperature of the entire system
rho15 = 886.2;     %Density of supplied oil at 15degC
V_kin_40 = 71.83;  % Kinematic viscosity of oil at 40degC
V_kin_100 = 9.072; % Kinematic viscosity of oil at 100degC

A_oil=log10(log10(V_kin_40 + 0.7)/log10(V_kin_100 + 0.7))/(log10(313/373));
B_oil = log10(log10(V_kin_40 +0.7))-A_oil*(log10(313));
Vkin = 10^((10^(A_oil*log10(T_supply)+B_oil)))+0.7; %kinematic viscosity at the supply temp
rho = rho15*(1-(0.7*(((T_supply)-288)/rho15))); %density of lubricant at supply temp
Mu_0 = (Vkin*10^-6)*rho; %%Initial dynamic viscosity of oil at supply temp


W=1000;       %Applied load on the bearing
L=49.3e-3;      %Length of bearing
Dbi=80e-3;   %Inner diameter of the bushing
a=50e-6;      %linear thermal expansion coefficient of the shaft
Dbo = 110e-3; % Outer Diameter of the bearing
Ds = Dbo-Co;   %Shaft diameter

Vs  = ((pi*Ds^2)/4)*L;  %Volume of the shaft
Vb = ((pi*(Dbo^2 - Dbi^2))/4)*L; % Volume of the bearing

rho_s = 7800; % Density of the shaft material
rho_b = 7800; % Density of the bearing material
cl = 2000; % Specific heat capacity of lubricant (guess)
cs = 456; % Specific heat capacity of shaft material
cb = 444; % Specific heat capcity of the bearing material
cc = cb+cs; % Specific heat capacity of bearing system (in accordance with lumped model)

ps=2e5; %supply pressure of the lubricant (relative to the atmosphere) 
M=(rho_b*cb*Vb+rho_s*cs*Vs)/cc; % Equivalent mass of the bearing

%%functions to get temperature
Tfun=@(C) 2*(Co-C)/(Dbi*a)+T_supply;    

%%eccentricity's influence on load capacity
fcp=@(e) ((1-e^2)^2)/(pi*e*sqrt(pi^2*(1-e^2)+16*e^2));

%dimensionless terms
Lambda = (6*pi^2*a*Dbi^3*L^5*n.^2*Mu_0^2)/(rho*cl*Co^2*ps);%coefficient related to friction power clearance losses
Gamma = (12*n*L^2*Mu_0)/(Co^2*ps); %Coefficient related to hydrodynamic flow rate clearance losses
tau= ((pi*Co^3*ps*rho*cl*Dbi)/(6*cc*M*L*Mu_0)); %Dimensionless time

C=Co;%initial clearance
e=0.1; %initial eccentricity

f_QHS_ecc=@(e) (2*pi*Dbi*C^3*ps/(12*Mu_0*L))*(1+1.5*e^2);
f_QHD_ecc=@(e) (pi*L*Dbi*n*C)*e;    
f_Pf_ecc =@(Mu_0,e) (2*pi*Dbi*L*Mu_0*(((pi*Dbi*n).^2)/(2*C)))*((2+e)/(2*(1+e)*(1-e^2)^0.5)); 

while true
  
%%define symbolic variable for differential equation
syms delta(t)

%%define the differential equation
fun =diff(delta,t)== tau*(Lambda*f_Pf_ecc(Mu_0,e)- delta*((1-delta)^2*f_QHS_ecc(e)+(Gamma*(1-delta)*f_QHD_ecc(e))));

%define initial condition of the differentiaon equation
cond=[delta(0)==0];

%%convert the differential equation to system of first order for ode45
V=odeToVectorField(fun,cond);
Mi = matlabFunction(V,'Vars',{'t','Y'});

%%time span 
tspan = [0 3600];

%%solve equation with ode45
sol = ode45(Mi,tspan,0);

%%check for load 
deltas=sol.y;
Cp=(1-deltas(end))*fcp(e);
load=(L^3*Dbi*Mu_0*n)/(2*Cp*Co^2);

fprintf('eccentricity =%.2f \n',e)
fprintf('Load =%.2f N \n',load)

%if load condition is valid then break loop 
if load>W
   break 
end

e=e+0.01;
%%New Viscosity = Mu_0*(C/Co) 
% How do I build this into an iterative loop? The condiditons to break would be if C>Co


end

figure(1)

time=sol.x;
ds=sol.y;

C=Co*(1-ds);

T=Tfun(C);

figure(1)
plot(time/60,ds)

xlabel('time (mins)')
ylabel('dimensionless loss of clearance')

figure(2)
plot(time/60,C*1000000)
xlabel('time (mins)')
ylabel('Radial Clearance C (um)')

figure(3)
plot(time/60,T-273)
xlabel('time (mins)')
ylabel('Temperature (DegC)')                              

T=Tfun(C);

figure(1)
plot(time/60,ds)

xlabel('time (mins)')
ylabel('dimensionless loss of clearance')

figure(2)
plot(time/60,C*1000000)
xlabel('time (mins)')
ylabel('Radial Clearance C (um)')

figure(3)
plot(time/60,T-273)
xlabel('time (mins)')
ylabel('Temperature (DegC)')