Error using ode45

Hi All,

I have recently written the following script, and after realising a couple of mistakes I get an error surrounding ode45. The error message I get is:

Error in ode45 (line 434)

f3d = cat(3,f3d,zeros(neq,7,chunk,dataType));

Error in file_name (line 80)

sol = ode45(Mi,tspan,0);

The only changes I made were in the definition of 'Lambda' which has gone from this:

Lambda = (6*pi^2*a*Dbi^3*L^5*n.^2*Mu_0^2)/(rho*cl*Co^2*ps)

to this:

Lambda = (6*pi^2*a*Dbi^3*L^2*n.^2*Mu_0^2)/(rho*cl*Co^5*ps)

If anyone can shed some light on why I am getting this error, I would be really grateful.

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^2*n.^2*Mu_0^2)/(rho*cl*Co^5*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;
e=0.1;

f_QHS_ecc=@(e) (2*pi*Dbi*C^3/(12*Mu_0)*(ps/L))*(1+1.5*e^2);
f_QHD_ecc=@(e) (2*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 differential 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 10];

%%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.1;

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)')                                 

[Torsten]

And did you check whether your code works now with the old definition of Lambda ?

[Max]

I have swapped the definition of Lambda back to its original format and it runs, but just the values cannot be correct due to the equation being wrong.

[Torsten]

Could you include your ODE function by inserting

matlabFunction(V,'Vars',{'t','Y'},'File','myfile');

after the line

Mi = matlabFunction(V,'Vars',{'t','Y'});

and attaching "myfile" here ?

[Max]

Thanks for your help so far Torsten, really appreciate it.

Please see below the output from "my file". Please note the code stopped after expressing the same error message originally stated.

function V = myfile(t,Y)
%MYFILE
%    V = MYFILE(T,Y)

%    This function was generated by the Symbolic Math Toolbox version 8.3.
%    02-Jul-2019 11:40:56

t2 = Y(1);
V = [t2.*(t2.*(-2.383503026633638e+1)+(t2-1.0).^2.*1.304316168028935e-5+2.383503026633638e+1).*(-3.734896195706258e-3)+1.768347707341557e+16];

[Torsten]

Then your code is equivalent to

function main
  %%time span 
  tspan = [0 10];  
  %%initial condition
  delta0 = 0;
  %%solve equation with ode45
  [T delta] = ode45(@myfile,tspan,delta0);
  plot(T,delta)
end

function V = myfile(t,Y)
  t2 = Y(1);
  V = [t2.*(t2.*(-2.383503026633638e+1)+(t2-1.0).^2.*1.304316168028935e-5+2.383503026633638e+1).*(-3.734896195706258e-3)+1.768347707341557e+16];
end

The constant 1.76...e+16 looks quite big.

[Max]

hmmm it does look a little large. I will revisit the equations and input data to see if theres anything glaring that I have missed.

Thanks again.