- #1
halleluhwah
- 2
- 0
Hello, I've made a script to solve a non-linear ordinary differential equation for the temperature (lumped-mass) of a brake rotor due to some constant deceleration of a race car (Qbrake in script). The problem I'm having is properly simulating the velocity of the car. What I'd like for it to do is start from 60mph (Vmax), decelerate (@ Gdec=-1.6 g's) to 30mph (Vmin), and accelerate (@ Gacc=0.5 g's) back to 60 mph, and repeat this process for the time span (a few minutes) of ODE45 which should drive the temperature to a roughly steady value as this velocity cycles between Vmin and Vmax. I've attempted this w/ IF conditionals and a persistent "gees" (the acceleration value changing based on velocity conditions), but running the script results in erratic velocity values over the time span. That is, it does not run smoothly between the velocity limits.
Would would be a better way to simulate this behavior within the ODE45 algorithm?
My script:
Thanks!
Would would be a better way to simulate this behavior within the ODE45 algorithm?
My script:
Code:
function thermalrotor
% Initial conditions
Vo=60*1.67; % mph --> ft/s
To=75; % degF, ambient
% ODE45
tspan=0:.1:120;
[t,y]=ode45(@thermalt,tspan,[Vo To])
Vc=y(:,1)/1.67
Tr=y(:,2);
Tm=2795; % melting temperature of ductile iron
plot(t,Tr,t,Tm,'g')
title('Temperature as a function of time, ductile iron rotor');
ylabel('Temperature, F');
xlabel('Time, s');
legend('Rotor Temp.','Melting Temp.');
clear
%% Subfunction
function dx=thermalt(t,x)
% Constants:
persistent gees
W=575; % lbf
g=32.18; % ft/s
Vmax=60*1.67; % ft/s
Vmin=30*1.67; % ft/s
Gacc=.5;
Gdec=-1.6;
As1=48.027/144; % ft^2
As2=38.799/144; % ft^2
Ah=4*0.2281/144; % ft^2
Lh=0.31846/12; % ft
M=0.9335; % lbm
Cp=0.128; % Btu/(lbm*degF)
k=46/3600; % Btu/(s*ft*degF)
h=25/3600; % Btu/(s*ft^2*degF)
fv=1.0;
fe=0.9;
sigma=1.714*10^-9/3600; % Btu/(ft^2*s*R^4)
Tamb=75; % degF
rad=fe*fv*sigma*As2;
if x(1)>=Vmax
gees=Gdec;
elseif x(1)<=Vmin
gees=Gacc;
end
% Energy absorbed by rotor from deceleration (Q=Fv):
Qbrake=-W/3*x(1)*gees/778;
% Energy radiated to air
Qrad=rad*((Tamb+460)^4-(x(2)+460)^4);
% x(1) = Vc = car velocity (linear)
% x(2) = Tr = lumped mass rotor temperature
dx=zeros(2,1);
dx(1)=gees*g;
dx(2)=(Qbrake+(h*As2+k*Ah/Lh)*(Tamb-x(2))+Qrad)/(M*Cp);
Thanks!