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

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);

      radius = ((3/(4*pi))*vol).^(1/3);
      Ca = pi*radius.^2;
      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;

Answer 1

[Ari]

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.

[Kellie Matthews]

Thanks, I'm not getting any errors now and the graph looks good!

[t,x]=ode45('F',[0,10800],[1,0]);
plot(t,x(:,1))

Now I just need to figure out how to model the balloon's descent. I think I want to decrease the volume by a little bit after these initial 10800 seconds, not sure how to fit that in or if I need another function.

[Ari]

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]

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]

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.