# Simulating velocity behavior w/ IF conditionals and ODE45

1. Mar 31, 2012

### 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:

Code (Text):
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

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;
% x(1) = Vc = car velocity (linear)
% x(2) = Tr = lumped mass rotor temperature
dx=zeros(2,1);
dx(1)=gees*g;
Thanks!

2. Mar 31, 2012

### 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).

3. Apr 1, 2012

### 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?