再入式飞行器轨道建模与仿真 （Modeling of trajectory）

再入式飞行器轨道建模与仿真，建立各种模型，获得飞行器的动力学模型，用龙格库塔函数求解该微分方程初值问题。

%%atmosphere.m大气参数函数文件

function Y = atmosphere（x CL）
r=x（1）;vel=x（4）;lat=x（2）;

R = 287; %海平面气体常数 （J/kg.K）
% 重力加速度 应该 随高度、纬度而变化，忽略经度影响
%%g0 = 9.806; 海平面重力加速度 （m/s^2）
go=gravity（rlat）;
Na = 6.0220978e23; %Avogadro常数
sigma = 3.65e-10; %空气分子运动平均自由程
S = 110.4; %Sutherland温度 （K）
Mo = 28.964; %海平面空气分子质量 （g/mole）
To = 288.15; %海平面开氏温度 （K）
Po = 1.01325e5; %海平面大气压 （N/m^2）

re = 6378.14e3; %地球平均半径 （m）

Beta = 1.458e-6; %Sutherland常数 （kg/m.s.K^0.5）
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;

h=r-re;
% 若是h<0， 这已不在Z的范围之内，故以h=0时的值代替
% h>max（Z）时，以max（Z）的值代替
if h<0
    h=0;
elseif h>max（Z）
    h=max（Z）;
end

% 找出 h 夹在 Z的那两个元素之间，这里里k表示小的那个序号
% 即： Z（k） < h <= Z（k+1）
k=find（Z<=h 1 ‘last‘ ）;
% 既然找到了 h 处于第几层， 就没有必要计算全部的层，只需要计算需要的参数即可。
for i =1:k
    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
% 同上，直接算
i=k;
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; %Static Pressure （N/m^2）
    C7 = C2 + 1;
    rhoE = C5*rho（i）*C3^（-C7）; %Density （kg/m^3）
else
    TM = T（i）;
    C8 = -go*（h - Z（i））*（1 - （h + Z（i））*B/2）/（R*T（i））;
    PR = P（i）*exp（C8）; %Static Pressure （N/m^2）
    rhoE = rho（i）*exp（C8）; %Density （kg/m^3）
end
MOL = M（i） + （ M（i+1）-M（i） ）*（ h - Z（i） ）/（ Z（i+1） - Z（i） ）;
TM = MOL*TM/Mo; %Kinetic Temperature
asound = sqrt（gamma*R*TM）; % Speed of Sound （m/s）
MU = Beta*TM^1.5/（TM + S）; % Dynamic Viscosity Coeff. （N.s/m^2）
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; % Mean free-path （m）
Mach = vel/asound; % Mach Number
T0 = TM*（1 + （gamma - 1）*Mach^2/2）;

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----
     文件        3553  2020-11-27 06:24  ╘┘╚δ╩╜╖╔╨╨╞≈╣∞╡└╜¿─ú╙δ╖┬╒µ\atmosphere.m
     文件        1427  2020-11-27 06:24  ╘┘╚δ╩╜╖╔╨╨╞≈╣∞╡└╜¿─ú╙δ╖┬╒µ\get_drag.m
     文件         723  2020-11-27 06:24  ╘┘╚δ╩╜╖╔╨╨╞≈╣∞╡└╜¿─ú╙δ╖┬╒µ\gravity.m
     文件        1750  2020-11-27 06:24  ╘┘╚δ╩╜╖╔╨╨╞≈╣∞╡└╜¿─ú╙δ╖┬╒µ\kinetic_equation.m
     文件        2528  2020-11-27 06:24  ╘┘╚δ╩╜╖╔╨╨╞≈╣∞╡└╜¿─ú╙δ╖┬╒µ\main_function.asv
     文件        2547  2020-11-27 06:24  ╘┘╚δ╩╜╖╔╨╨╞≈╣∞╡└╜¿─ú╙δ╖┬╒µ\main_function.m