Iterating final temperature minus initial temperature in density calculation

Hello all, Basically i have a matrix of temperatures, ranging from 223.15K to 3273.15K, and i am trying to calculate the final density of air. The formula i have for the change in density is: ρ1 = ρ0 / (1 + β (t1 - t0)) where ρ1 = final density (kg/m3) ρ0 = initial density (kg/m3) β = volumetric temperature expansion coefficient (m3/m3 K) t1 = final temperature t0 = initial temperature

So what i want to do is input the second element for my temperature matrix into t1, and the first element into t0. Then the next calculation i want to use the third element minus the second element and so on until the last element. Also, use the final density from the first calculation and input it into the initial density into the second equation and so on. Ive been trying to set it up with for and while loops and have gotten nowhere. Any help with the code would greatly be appreciated.

[John D'Errico]

Be clear in your question. Is ρ0 the density at temperature t0? Then ρ1 is the density at temperature t1?

[Doug]

Yes, for the first calculation the intial density is at t0, the next density is at t1. For the second calculation, ρ2 = ρ1 / (1 + β (t2 - t1)), where ρ2 will be calculated from ρ1 from the first calculation, and t2 would be the third element from the temperature matrix.

Answer 1

[Star Strider]

There are likely more efficient ways to do this, but if you want to use a loop, this would be my approach to your problem:

beta = 0.5;                                         % Beta (Obviously)
T = linspace(223.15, 3273.15);                      % Temperature Vector
DT = T(2:end) - T(1:end-1);                         % Temperature Difference Vector
p(1) = 1.5;                                         % Initial Air Density (‘rho’) At T(1)
for k1 = 2:length(DT)-1
    p(k1) = p(k1-1)/(1 +  beta*DT(k1));
end

figure(1)
semilogy(T(2:end-1), p)
grid

I’m obviously guessing at the data you didn’t provide, and making some assumptions.

[Doug]

Beta is a function of temperature where beta=1/T. Initial density is just some constant value dependent upon the fluid. Since beta is a vector i dont think this code will work. I think im just going to use a polynomial curve fit of order 3 that just takes a temperature input. Thank you for your help though.

[Star Strider]

My pleasure.

If beta = 1/T, that still makes the code straightforward:

T = linspace(223.15, 3273.15);                      % Temperature Vector
DT = T(2:end) - T(1:end-1);                         % Temperature Difference Vector
p(1) = 1.5;                                         % Initial Air Density (‘rho’) At T(1)
for k1 = 2:length(DT)-1
    p(k1) = p(k1-1)/(1 +  DT(k1)/T(k1));
end

figure(1)
plot(T(2:end-1), p)
grid

This does work, although you may want to subtract 1 from the ‘T’ subscript instead. I’m not familiar with the physics of your simulation (I know about density altitude, but that doesn’t seem to apply at these temperatures).

For example, this might be more like what you want:

p(k1) = p(k1-1)/(1 +  DT(k1)/T(k1-1));

Also, consider taking the mean of [T0,T1] or using the interp1 function if you want intermediate temperature values for the denominator.