air_data.m

导弹控制系统的鲁棒控制设计实例

function [p,ro,M,Tm] = air_data(V,H)
% Air parameters as function of altitude
%
% H                     % m       altitude above sea level
g0  = 9.80665;          % m/s^2   acceleration of gravity at sea level
T0  = 288.15;           % K       air temperature at sea level
p0  = 101325;           % N/m^2   air pressure at sea level
ro0 = 1.225;            % kg/m^3  air density at sea level
a0  = 340.29;           % m/s     speed of sound
R   = 287.039;          % J/K.kg  specific gas constant
Re  = 6371210;          % m       radius of Earth
%
% Acceleration of gravity
g = g0*(Re/(Re+H))^2;
%
% Geopotential altitudes
f = H*Re/(Re+H);
f11  =  11000*Re/(Re+11000);
f25  =  25000*Re/(Re+25000);
f46  =  46000*Re/(Re+46000);
f54  =  54000*Re/(Re+54000);
f80  =  80000*Re/(Re+80000);
f95  =  95000*Re/(Re+95000);
f110 = 110000*Re/(Re+110000);
f120 = 120000*Re/(Re+120000);
f150 = 150000*Re/(Re+150000);
f160 = 160000*Re/(Re+160000);
f200 = 200000*Re/(Re+200000);
%
% Temperature gradients
a11  = -0.00651122;
a46  =  0.00276098;
a80  = -0.00349544;
a110 = 0.005;
a120 = 0.00801741;
a150 = 0.02346357;
a160 = 0.01987408;
a200 = 0.00308461;
%
% Molecular temperature and air pressure
T11  = T0 + a11*f11;
p11  = exp(log(p0) - log(T11/T0)*g0/(a11*R));
T25  = 216.66;
p25  = exp(log(p11) - g0*(f25 - f11)/(R*T25));
T46  = T25 + a46*(f46 - f25); 
p46  = exp(log(p25) - log(T46/T25)*g0/(a46*R));
T54  = 274.0;
p54  = exp(log(p46) - g0*(f54 - f46)/(R*T54));
T80  = T54 + a80*(f80 - f54);
p80  = exp(log(p54) - log(T80/T54)*g0/(a80*R));
T95  = 185.0;
p95  = exp(log(p80) - g0*(f95 - f80)/(R*T95));
T110 = T95 + a110*(f110 - f95);
p110 = exp(log(p95) - log(T110/T95)*g0/(a110*R));
T120 = 257.64 + a120*(f120 - f110);
p120 = exp(log(p110) - log(T120/T110)*g0/(a120*R));
T150 = 335.0 + a150*(f150 - f120);
p150 = exp(log(p120) - log(T150/T120)*g0/(a150*R));
T160 = 1010.0 + a160*(f160 - f150);
p160 = exp(log(p150) - log(T160/T150)*g0/(a160*R));
T200 = T160 + a200*(f200 - f160);
p200 = exp(log(p160) - log(T200/T160)*g0/(a200*R));
%
if H <= 11000
   Tm = T0 + a11*f;
   p = exp(log(p0) - log(Tm/T0)*g0/(a11*R));
elseif H <= 25000
   Tm = T25;
   p = exp(log(p11) - g0*(f - f11)/(R*Tm));
elseif H <= 46000
   Tm = T25 + a46*(f - f25); 
   p = exp(log(p25) - log(Tm/T25)*g0/(a46*R));
elseif H <= 54000
   Tm = T54;
   p = exp(log(p46) - g0*(f - f46)/(R*Tm));
elseif H <= 80000
   Tm = T54 + a80*(f - f54);
   p = exp(log(p54) - log(Tm/T54)*g0/(a80*R));
elseif H <= 95000
   Tm = T95;
   p = exp(log(p80) - g0*(f - f80)/(R*Tm));
elseif H <= 110000
   Tm = T95 + a110*(f - f95);
   p = exp(log(p95) - log(Tm/T95)*g0/(a110*R));
elseif H <= 120000
   Tm = T110 + a120*(f - f110);
   p = exp(log(p110) - log(Tm/T110)*g0/(a120*R));
elseif H <= 150000
   Tm = T120 + a150*(f - f120);
   p = exp(log(p120) - log(Tm/T120)*g0/(a150*R));
elseif H <= 160000
   Tm = T150 + a160*(f - f150);
   p = exp(log(p150) - log(Tm/T150)*g0/(a160*R));
elseif H <= 200000
   Tm = T160 + a200*(f - f160);
   p = exp(log(p160) - log(Tm/T160)*g0/(a200*R));
elseif H > 200000
   Tm = T200; 
   p = p200; 
end 
%
% Air density
ro = p/(R*Tm);                                  % kg/m^3   
%
% Air speed
a = sqrt(1.4*R*Tm);                             % m/s   
%
% Mach number
M = V/a;
%
% Coefficient of dynamic viscosity
mu_dyn = 1.458*10^(-6)*Tm^(3/2)/(Tm + 110.4);   % N.s/m^2
%
% Coefficient of kinematic viscosity
nu = mu_dyn/ro;                                 % m^2/s