ODE to solve for density

7 views (last 30 days)
Christian Boaitey
Christian Boaitey on 2 Mar 2021
Hi,
I'm currently trying to calculate air density which involves solving an ODE for altitude. I have plotted a graph for it, I get the right shape but the axis and numbers are wrong. The first figure is what I get for my output, whereas the second is the correct one (using the line vector). Can anybody help me? Much appreciated
I get this for my output
[h,p]=ode45(@p_prime,[0 7.5],0);
plot(p,h)
ylabel('Altitude (km)')
xlabel('Air density (kg/m^3)')
function dp=p_prime(h,p)
p0=1.225
hscale=7.5
dp=(p0*exp(-h/hscale))
end

Accepted Answer

Mathieu NOE
Mathieu NOE on 2 Mar 2021
hello
as far as i can tell , the second plot is simply the derivative given by you and not the solution of the ODE
hx = linspace(0,7.5,100);
p0=1.225;
hscale=7.5;
dp=(p0*exp(-hx/hscale)); % derivative
figure(1), plot(dp,hx);
The curve altitude /density should have a negative slope , so at least there is a sign error in the derivative equation, drho/dh must be negative - beside that I don't see consistence with what I read here :
  4 Comments
Christian Boaitey
Christian Boaitey on 2 Mar 2021
Edited: Christian Boaitey on 2 Mar 2021
The equation we were given was po*e^(h/hscale), this is an approximation of density as a function of altitude. These come from the curved earth equations. Other than that, I'm still confused as to why the plot is reversed?
Bjorn Gustavsson
Bjorn Gustavsson on 2 Mar 2021
You get the plot you get because you interpret things wrong. That expression fo the pressure-variation are the solution to the ODE I gave you below. You can derive that ODE from the equation-of-motion for a slice of atmosphere under steady-state conditions (air-parcel at rest so, no acceleration) for which the difference bewteen the pressure-force upward from below and downward from above is balanced by the force of gravity on the air-parcel. That leads to the ODE below. That ode you can implement in a matlab-function, then with some initial condition, like ground-pressure, you can call ode45 and obtain a solution for the pressure-variation.
Where you go wrong is that you plug in the solution of the barometric ODE in a function for an ODE and integrate that one, starting with an initial condition for the ground-pressure that is zero (which in itself is wrong, we have approximately 1 atm pressure at ground-level...).
The dp/dz should be strictly negative, since the weight of each atmospheric layer is kept steady by a correspondinly larger pressure from below than from above. Your pressure-gradient is set to be larger than zero a all altitudes which is unphysical. This also does not come from any "curved earth equation" (other than that gravity in itself makes planetary-sized bodies rather curved...)

Sign in to comment.

More Answers (1)

Bjorn Gustavsson
Bjorn Gustavsson on 2 Mar 2021
Your ODE is not correct for the pressure-variation. What you have in the p_prime function is closer to the solution of the pressure-balance. If you go back to the book you'll see something like this:
That is what you should have in the p_prime-function that you then integrate with altitude.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!