Trying to plot viscosity and kinematic viscosity vs altitude

close all
clear
clc

temp_sl = 288; %temperature at sea level 
rho_sl = 1.225; %density at sea level

alt = 1:1:47000;
temp_alt = zeros(1,47000);
d = zeros(1,47000);
d_viscosity = zeros(1,47000);
k_viscosity = zeros(1,47000);

%temperature in Kelvin
for i = 1:47000
    if alt(i) < 11000
        temp_alt(i) = temp_sl - 0.0065*alt(i);
    elseif alt(i) > 10999 && alt(i) < 25000
        temp_alt(i) = 216.54;
    else
        temp_alt(i) = 191.5 + 0.001*alt(i);
    end
end

%kinematic viscosity 
for i = 1:47000
    d_viscosity(i) = 2.287973*10^(-6) + 6.259793*(10^(-8))*temp_alt(i)-3.131956*(10^(-11))*temp_alt(i)^(2) + 8.15038*(10^(-11))*temp_alt(i)^(3);
    d(i) = rho_sl*(temp_alt(i)/temp_sl).^(4.2588);
    k_viscosity(i) = d_viscosity(i)/d(i);
end
%%%%%%%%%%%%%%%
plot(d_viscosity, alt, "DisplayName","Dynamic Viscosity")
hold on
plot(k_viscosity, alt, "DisplayName","Kinematic Viscosity")
hold off
       

[Torsten]

temp_alt can't be a function handle

temp_alt = @(T_sl,h) (T_sl - 0.0065*h);

and an array

temp_alt(i) = temp_sl - 0.0065*m;

at the same time.

T_sl in the line

temp_alt(i) = 0.0028*m+((0.001*32000+(T_sl - 0.0065*11000-0.001*20000))-0.0028*32000);

is undefined.

[Torsten]

temp_alt = zeros(1,47000);

instead of

temp_alt = zeros(1:47000);

Same for the other variables.

T undefined in

d_viscosity(i) = (2.287973*10^-6) + (6.259793*10^-8)*T(i) - (3.131956*10^11)*T(i)^2 + (8.15038*10^-11)*T(i)^3;

temp undefined in

d(i) = rho_sl*(temp(i)/temp_sl).^(4.2588);

And it's ok now ?

[Roxanne Williamson]

I'm just struggling to space the graph nicely but the values look right. Your assistance is greatly appreciated!

[William Rose]

You're welcome, @Roxanne Williamson.

Answer 1

[William Rose]

@Roxanne Williamson,

When I run your code I get the plot below. What do you not like about it? Th only revommendation I have is to add a legend, which I have done. I changed the plotting order so that Altitude is on the horizontal axis and viscosity is on the vertical axis.

%viscosity_RW.m
clear
temp_sl = 288; %temperature at sea level 
rho_sl = 1.225; %density at sea level
alt = 1:1:47000;
temp_alt = zeros(1,47000);
d = zeros(1,47000);
d_viscosity = zeros(1,47000);
k_viscosity = zeros(1,47000);
%temperature in Kelvin
for i = 1:47000
    if alt(i) < 11000
        temp_alt(i) = temp_sl - 0.0065*alt(i);
    elseif alt(i) > 10999 && alt(i) < 25000
        temp_alt(i) = 216.54;
    else
        temp_alt(i) = 191.5 + 0.001*alt(i);
    end
end
%kinematic viscosity 
for i = 1:47000
    d_viscosity(i) = 2.287973*10^(-6) + 6.259793*(10^(-8))*temp_alt(i)-3.131956*(10^(-11))*temp_alt(i)^(2) + 8.15038*(10^(-11))*temp_alt(i)^(3);
    d(i) = rho_sl*(temp_alt(i)/temp_sl).^(4.2588);  %denisty
    k_viscosity(i) = d_viscosity(i)/d(i);           %kinematic = dynamic/density
end
%%%%%%%%%%%%%%%
plot(alt,d_viscosity,'-r*',alt,k_viscosity,'-bs');
title('Viscosity');
legend('Dynamic','Kinematic');
xlabel('Altitude')

See above.

[William Rose]

@Roxanne Williamson

I uploaded that answer before I was done. Th escript in my answer generates the plot below. I would add a y-axis label, and units for each axis, and maybe a grid.

I realize now that your way of plotting shows altitude going up and down, which is nice.