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

Question

Kellie Matthews 2017-7-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/350667-using-ode45-for-high-altitude-balloon-trajectory-need-some-variables-to-update-every-iteration-and

编辑： Torsten 2017-7-31

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;

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Ari 2017-7-28

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/350667-using-ode45-for-high-altitude-balloon-trajectory-need-some-variables-to-update-every-iteration-and#answer_275933

编辑：Ari 2017-7-28

在 MATLAB Online 中打开

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Kellie Matthews 2017-7-31

在 MATLAB Online 中打开

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 2017-7-31

编辑：Torsten 2017-7-31

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.

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论