Is it the code or my laptop?

Hi All,

I am trying to run the code I have created on a work computer, however after some changes ot the code I now wish to run it on my academic version of Matlab 2017a.

My laptop is useless at the best of times, 4GB of RAM, so I was wonderring whether I just need to be running the code on a more powerful machine.

After leaving it over night, I have managed to atleast get it to run but I got the following error message:

Error using cat

Requested 1x7x74350200 (3.9GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a

long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.

Error in ode45 (line 380)

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

Error in TG (line 86)

sol = ode45(Mi,tspan,0);

I'm not sure what I can do to help 'slim down' the script, but if it is possibl for someone with a more capable machine to just try running it for me, I would really appreciate it.

FYI matlab online just times out when I try and run it on there.

Script supplied below:

All the best.

clc
clear all
close all


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

n=3600;  %Rotational speed of the shaft in rpm 
%n=n*0.107; %Rotational speed of the shaft in rad/s
Co=0.000125; %initial radial clearance (m)
T_supply = 40; % 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+273)+B_oil)))+0.7; %kinematic viscosity at the supply temp
rho = rho15*(1-(0.7*(((T_supply+273)-289)/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
B=80;      %Width of one bearing land
Dbi=159.725;   %Inner diameter of the bushing
a=50;      %linear thermal expansion coefficient of the shaft
Dbo = 179.725; % Outer Diameter of the bearing
Ds = Dbo-Co;   %Shaft diameter
Lb = B;     %length of the bearing


Vs  = ((pi*Ds^2)/4)*Lb;  %Volume of the shaft
Vb = ((pi*(Dbo^2 - Dbi^2))/4)*Lb; % 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 = 250; % 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*B^5*n.^2*Mu_0^2)/(rho*cl*Co^2*ps);%coefficient related to friction power clearance losses
Gamma = (12*n*B^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*B*Mu_0)); %Dimensionless time
 

C=Co;
e=0.1;

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

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

%%check for load 
deltas=sol.y;
Cp=(1-deltas(end))*fcp(e);
load=(B^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,ds)

xlabel('time/60')
ylabel('dimensionless loss of clearance')

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


figure(3)
plot(time,T)
xlabel('time/60')
ylabel('Temperature (kelvin)')

[KALYAN ACHARJYA]

Please do read.