# Using ode45 for high altitude balloon trajectory: need some variables to update every iteration and need to plot altitude vs time.

12 views (last 30 days)
Kellie Matthews on 28 Jul 2017
Edited: Torsten on 31 Jul 2017
I'm trying to predict the trajectory of a high altitude balloon. This is the second order ODE I am using: (mass+cb*rho*vol)*z"=g(rho*vol-mass)-.5*rho*Cd*z'*|z'|*Ca
cb, Cd, g, and mass are constants. density[rho],volume, (approximate) surface area of a circle [Ca] all change with altitude.
I used this link to help me set up my second order ODE as a function: http://www.math.purdue.edu/~shen/cs614/projects/ode45.pdf
function xp=F(t,x) %xp = x prime or x'
xp=zeros(2,1); %output must be column vector
xp(1)=x(2);
xp(2)=(g*(RhoA*vol-mass)-.5*RhoA*realCD*x(2)*abs(x(2))*Ca)/(mass+cb*RhoA*vol);
end
This is what I need for my atmospheric properties that rely on altitude:
if (z <= 11000) %Meters (Troposhpere)
temp = 15.04 - 0.00649*z;
tempK = temp + 273.15;
p = 101.29*((temp+273.1)/(288.08)).^5.256; %kPa
elseif (z > 11000 && z < 25000) %Meters (Lower Stratosphere)
temp = -56.46;
tempK = temp + 273.15;
p = 22.65*exp(1.73-0.000157*z); %kPa
else %Upper Stratosphere
temp = -131.21 + 0.00299*z;
tempK = temp + 273.15;
p = 2.488 * ((temp+273.1)/216.6).^-11.388; %kPa
end
dTempK = abs(tempK - oldTempK);
RhoA = (p/(.2869*tempK));
Wg = Mb.*(1000*p).*vol/(r.*tempK);
old_z = z;
[t,x]=ode45('F',[0,tf],[0,0]);
hold on
plot(t,x(:,1))
z=x(i,1);
dz = z - old_z; %this is the change in altitude from the last second
dVol = (r/(p*Mb))*(Wg*dTempK/dt)*dt + (RhoA/p)*(vol)*dz;
vol = vol + dVol;

Ari on 28 Jul 2017
Edited: Ari on 28 Jul 2017
or variables that change with time or state you should put their calculations inside the function xp. In your case, your states x seem to be [z;z']. Set z = x(1) in the beginning of the function and calculate the variable parameters before you calculate xp(2). It seems you will run into a problem trying to access the z value of the previous timestep (old_z) unless you use a persistent variable. You can try the following.
function xp = F(t,x)
persistent old_z;
z = x(1);
dz = z - old_z;
old_z = z; % set the value of old_z for next timestep
% calculate variable parameters
...
xp = zeros(2,1);
xp(1) = x(2);
xp(2) = (g*(RhoA*vol-mass)-.5*RhoA*realCD*x(2)*abs(x(2))*Ca)/(mass+cb*RhoA*vol);
end
The persisent variable will remain in memory between calls to the function.

Show 1 older comment
Ari on 28 Jul 2017
Inside the function:
if t >= 10800
% code to decrease volume here
end
Make sure you change the time span of simulation to a larger time span while calling the ode45 function.
Kellie Matthews on 31 Jul 2017
I'm sorry, but this didn't work for me. After xp(2)=%ascent eq.% I have
dz = z - old_z; %this is the change in altitude from the last second
dVol = (r/(p*Mb))*(Wg*dTempK/dt)*dt + (RhoA/p)*(vol)*dz;
vol = vol + dVol;
However, dVol and vol are underlined in red and matlab says "The value assigned to variable '..' might be unused" Does this mean that the volume is not updated throughout the code? I tried the if statement you suggested and my graph came out exactly the same as before, with no apparent decrease in altitude caused by a decrease in volume.
This is what I tried with an if statement:
dz = z - old_z; %this is the change in altitude from the last second
if t <= 10800
dVol = (r/(p*Mb))*(Wg*dTempK/dt)*dt + (RhoA/p)*(vol)*dz;
vol = vol + dVol;
else
vol = vol - 10;
end
The variables dVol and vol are still underlined in red when I try this.
Torsten on 31 Jul 2017
I wonder why you don't solve an additional (third) ODE for "vol" together with the two ODEs for height and velocity:
dVol/dt = (r/(p*Mb))*(Wg*dTempK/dt) + (RhoA/p)*(Vol)*dz/dt ;
This way, you overcome all the problems from above.
Best wishes
Torsten.