atmosphere.m

大气层内航天器飞行动力学matlab仿真模型

% Copyright Ashish Tewari (c) 2006
function Y = atmosphere(h, vel, CL)
R = 287; 
go = 9.806; 
Na = 6.0220978e23; 
sigma = 3.65e-10; 
S = 110.4;  
Mo = 28.964; 
To = 288.15; 
Po = 1.01325e5; 
re = 6378.14e3; 
Beta = 1.458e-6; 
gamma = 1.405; 
B = 2/re;
layers = 21;
Z = 1e3*[0.00; 11.0191; 20.0631; 32.1619; 47.3501; 51.4125; 
    71.8020; 86.00; 100.00; 110.00; 120.00; 150.00; 160.00; 170.00; 190.00;
    230.00; 300.00; 400.00; 500.00; 600.00; 700.00; 2000.00];
T = [To; 216.65; 216.65; 228.65; 270.65; 270.65; 214.65; 186.946;
    210.65; 260.65; 360.65; 960.65; 1110.60; 1210.65; 1350.65; 1550.65;
    1830.65; 2160.65; 2420.65; 2590.65; 2700.00; 2700.0];
M = [Mo; 28.964; 28.964; 28.964; 28.964; 28.964; 28.962; 28.962; 28.880;
    28.560; 28.070; 26.920; 26.660; 26.500; 25.850; 24.690;
    22.660; 19.940; 17.940; 16.840; 16.170; 16.17];
LR = [-6.5e-3; 0; 1e-3; 2.8e-3; 0; -2.8e-3; -2e-3;
    1.693e-3; 5.00e-3; 1e-2; 2e-2; 1.5e-2; 1e-2; 7e-3; 5e-3; 4e-3;
    3.3e-3; 2.6e-3; 1.7e-3; 1.1e-3; 0];
rho0 = Po/(R*To);
P(1) = Po;
T(1) = To;
rho(1) = rho0;
for i = 1:layers
    if ~(LR(i) == 0)
        C1 = 1 + B*( T(i)/LR(i) - Z(i) );
        C2 = C1*go/(R*LR(i));
        C3 = T(i+1)/T(i);
        C4 = C3^(-C2);
        C5 = exp( go*B*(Z(i+1)-Z(i))/(R*LR(i)) );
        P(i + 1) = P(i)*C4*C5;
        C7 = C2 + 1;
        rho(i + 1) = rho(i)*C5*C3^(-C7);
    else
        C8 = -go*(Z(i+1)-Z(i))*(1 - B*(Z(i + 1) + Z(i))/2)/(R*T(i));
        P(i+1) = P(i)*exp(C8);
        rho(i+1) = rho(i)*exp(C8);
    end
end

for i = 1:21
    if h < Z(i+1)
        if ~(LR(i)== 0)
            C1 = 1 + B*( T(i)/LR(i) - Z(i) );
            TM = T(i) + LR(i)*(h - Z(i));
            C2 = C1*go/(R*LR(i));
            C3 = TM/T(i);
            C4 = C3^(-C2);
            C5 = exp( B*go*(h - Z(i))/(R*LR(i)) );
            PR = P(i)*C4*C5;
            C7 = C2 + 1;
            rhoE = C5*rho(i)*C3^(-C7);
        else
            TM = T(i);
            C8 = -go*(h - Z(i))*(1 - (h + Z(i))*B/2)/(R*T(i));
            PR = P(i)*exp(C8);
            rhoE = rho(i)*exp(C8);
        end      
        MOL = M(i) + ( M(i+1)-M(i) )*( h - Z(i) )/( Z(i+1) - Z(i) );
        TM = MOL*TM/Mo;
        asound = sqrt(gamma*R*TM); 
        MU = Beta*TM^1.5/(TM + S); 
        KT = 2.64638e-3*TM^1.5/(TM + 245.4*10^(-12/TM));
        Vm = sqrt(8*R*TM/pi);
        m = MOL*1e-3/Na;
        n = rhoE/m;
        F = sqrt(2)*pi*n*sigma^2*Vm;
        L = Vm/F; 
        Mach = vel/asound; 
        T0 = TM*(1 + (gamma - 1)*Mach^2/2);
        MU0 = Beta*T0^1.5/(T0 + S);
        RE0 = rhoE*vel*CL/MU0;
        RE = rhoE*vel*CL/MU; 
        Kn = L/CL; 
        Kno = 1.25*sqrt(gamma)*Mach/RE0;
        if Kn >= 10
            d = 1;
        elseif Kn <= 0.01 
            d = 2;
        else         
            d = 3;
        end
        Y = [TM; rhoE; Mach; Kn; asound; d; PR; MU; RE];
        return;
    end
end