Simulating velocity behavior w/ IF conditionals and ODE45

[halleluhwah]

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:

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!

 

[Filip Larsen]

Welcome to PF!

You should in general not expect numerical integration of discontinuous fields to produce accurate solutions, and since your field is discontinuous I would not be surprised if that is why you see weird results (this does of course not exclude the possibility that you get weird results because of some other issues in your code).

If you want to integrate a "multi-mode" field in a more exact way using RK45 or similar, you should split the field up into multiple fields such that each field is continuous and sufficiently smooth. In your case that would be two fields, one for acceleration and one for braking. Starting with the field for a particular mode (deceleration in your case) you then integrate a solution (using only the field for that mode) until the trajectory reach a boundary to another mode (in your case until speed equals a certain speed) and then continue integrating the solution with the field for this new mode and so forth.

If your mode boundary conditions had been determined by time alone (the independent variable) then you could just have integrated to that time in order to get a solution point on the boundary. However, in your case the boundary conditions between your two modes are determined by speed alone (a dependent variable) so you will have to use interpolation, bisection or similar to find a solution point that is close enough to the boundary value (the initial solution samples from the RK45 will in general not be close enough by themselves since the sample distance is controlled by error measures and not by any effort to "hit" the boundary).

 

[halleluhwah]

Thanks for this Filip. I did not realize this was possible using ODE45. This seems like the solution to my problem, and while I understand generally the process you're describing, I'm not sure I understand how I would actually go about splitting a multi-mode field into two fields within the algorithm of ODE45 - I've done a search for an example of how to do this, but I'm not seeing any documentation. Could you point me in the right direction?

Thanks again for your help.

 

[Filip Larsen]

I'm not sure what would be the best way to solve such systems in general using matlab. When I was involved with such stuff (20+ years ago) we coded a RKF78 in C++ with the capability to provide high order interpolation with only a few extra field evaluations extra on top of those required by the solver (cannot at the moment recall what this technique was called). The interpolations was used in the phase when the solver had advanced the solution from a point that was within the state boundaries for the current mode to a point that lay outside it, where bisection or a Newton-Ralphson like method applied to some kind of boundary distance function was used to find a point "on the boundary".

However, for your particular problem it will probably be possible to skip the interpolation I mentioned above by noting that the speed is changed by constant acceleration and thus it is possible to calculate in advance precisely when the speed in either mode will reach the boundary to the next mode. So one method for you should be to start at some initial value in some mode (like decelerating), calculate how long this mode will last (time until speed is at limit), solve with RK45 to this time, change to the other mode, calculate new time interval using last point of previous solution as initial condition, and so forth.