12 views (last 30 days)

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;

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.

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.

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.

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 0 Comments

Sign in to comment.