2D Re-entry Model Figure errors

18 visualizaciones (últimos 30 días)
Michael Fitzpatrick
Michael Fitzpatrick el 29 de Nov. de 2019
Comentada: Michael Fitzpatrick el 29 de Nov. de 2019
The figures for my dynamic model are showing nothing could anyone explain why? Here's the code.
altitudeerror.jpg
clc
clear
%Constants
mu = 3.986e14; % Gravitational Parameter of Earth
Re = 6378e3; % Radius of Earth, m
m = 5000; % Mass of Vehicle, kg
Sref = 100; % Surface Area of Spacecraft, m^2
Rn = 5; % Nose Radius, m
%Initial Conditions
v0 = 10766; % Re-entry Velocity, m/s
r0 = 6878e3; % Re-entry Altitude, m
gamma0 = 9*pi/180; %Re-entry Flight Path Angle, rad
alpha = 25*pi/180;
x0 = [r0; 0; -v0*sin(gamma0); v0*cos(gamma0)/r0];
t0=0; tf = 350;
t = [t0,tf];
options = odeset('reltol',1e-7,'abstol',1e-7);
[T,X] = ode45(@(t,x) dynamic_EoM(t,x,m,Sref,Rn,alpha,mu,Re), [t0,tf], x0,options);
r = X(:,1);
theta = X(:,2);
vr = X(:,3);
vt = X(:,4);
%% Figures
% Altitude vs Time
figure(1) = figure;
plot(T,X(:,1)-Re);
title ('Altitude');
ylabel('Altitude (m)');
xlabel ('Time (seconds)');
Dynamic Model Script
function dx = dynamic_EoM(t,x,m,Sref,Rn,alpha,mu,Re)
% State Vector
r = x(1);
theta = x(2);
vr = x(3);
vt = x(4);
h = r - Re; % Altitude
gamma = tan(vr./vt);
v = sqrt(vr^2 + vt^2); % Velocity Vector
% Atmospheric Properies
[T, P, rho] = standard_atm(h);
[L, D] = Lift_Drag_Components(h, alpha, v, rho, Sref);
% Rate of Heat Transfer
q = heat_flux(rho, v, Rn);
dx(1) = vr;
dx(2) = vt;
dx(3) = -((mu)/r^2) + vt^2/r + 1/m*(-D*sin(gamma)+L*cos(gamma));
dx(4) = -(vr*vt)/r + 1/m*(-D*cos(gamma)-L*sin(gamma));
dx = [dx(1); dx(2); dx(3); dx(4)];
end
Atmosphere, Lift and Drag
function [T, P, rho] = standard_atm(h)
%{
Function implementing the International Standard Atmospheric model, which
relies on look-up tables.
Currently, the function looks up the closest altitude in the table to the
one asked for in the input (alt), and returns that value. A better method
would be to assume a linear or polynomial interpolation for intermediate
values.
INPUTS
alt = altitude above Earth surface, m
OUTPUTS
T = atmospheric temperature (K)
P = atmospheric pressure (Pa)
rho = atmospheric density (kg/m3)
%}
% make sure you are working in the correct units! The table below uses km,
% not m for altitude
% alt sigma delta theta temp press dens a visc k.visc
% km K N/sq.m kg/cu.m m/s kg/m-s sq.m/s
satm = [-2 1.2067E+0 1.2611E+0 1.0451 301.2 1.278E+5 1.478E+0 347.9 18.51 1.25E-5;...
0 1.0000E+0 1.0000E+0 1.0000 288.1 1.013E+5 1.225E+0 340.3 17.89 1.46E-5;...
2 8.2168E-1 7.8462E-1 0.9549 275.2 7.950E+4 1.007E+0 332.5 17.26 1.71E-5;...
4 6.6885E-1 6.0854E-1 0.9098 262.2 6.166E+4 8.193E-1 324.6 16.61 2.03E-5;...
6 5.3887E-1 4.6600E-1 0.8648 249.2 4.722E+4 6.601E-1 316.5 15.95 2.42E-5;...
8 4.2921E-1 3.5185E-1 0.8198 236.2 3.565E+4 5.258E-1 308.1 15.27 2.90E-5;...
10 3.3756E-1 2.6153E-1 0.7748 223.3 2.650E+4 4.135E-1 299.5 14.58 3.53E-5;...
12 2.5464E-1 1.9146E-1 0.7519 216.6 1.940E+4 3.119E-1 295.1 14.22 4.56E-5;...
14 1.8600E-1 1.3985E-1 0.7519 216.6 1.417E+4 2.279E-1 295.1 14.22 6.24E-5;...
16 1.3589E-1 1.0217E-1 0.7519 216.6 1.035E+4 1.665E-1 295.1 14.22 8.54E-5;...
18 9.9302E-2 7.4662E-2 0.7519 216.6 7.565E+3 1.216E-1 295.1 14.22 1.17E-4;...
20 7.2578E-2 5.4569E-2 0.7519 216.6 5.529E+3 8.891E-2 295.1 14.22 1.60E-4;...
22 5.2660E-2 3.9945E-2 0.7585 218.6 4.047E+3 6.451E-2 296.4 14.32 2.22E-4;...
24 3.8316E-2 2.9328E-2 0.7654 220.6 2.972E+3 4.694E-2 297.7 14.43 3.07E-4;...
26 2.7964E-2 2.1597E-2 0.7723 222.5 2.188E+3 3.426E-2 299.1 14.54 4.24E-4;...
28 2.0470E-2 1.5950E-2 0.7792 224.5 1.616E+3 2.508E-2 300.4 14.65 5.84E-4;...
30 1.5028E-2 1.1813E-2 0.7861 226.5 1.197E+3 1.841E-2 301.7 14.75 8.01E-4;...
32 1.1065E-2 8.7740E-3 0.7930 228.5 8.890E+2 1.355E-2 303.0 14.86 1.10E-3;...
34 8.0709E-3 6.5470E-3 0.8112 233.7 6.634E+2 9.887E-3 306.5 15.14 1.53E-3;...
36 5.9245E-3 4.9198E-3 0.8304 239.3 4.985E+2 7.257E-3 310.1 15.43 2.13E-3;...
38 4.3806E-3 3.7218E-3 0.8496 244.8 3.771E+2 5.366E-3 313.7 15.72 2.93E-3;...
40 3.2615E-3 2.8337E-3 0.8688 250.4 2.871E+2 3.995E-3 317.2 16.01 4.01E-3;...
42 2.4445E-3 2.1708E-3 0.8880 255.9 2.200E+2 2.995E-3 320.7 16.29 5.44E-3;...
44 1.8438E-3 1.6727E-3 0.9072 261.4 1.695E+2 2.259E-3 324.1 16.57 7.34E-3;...
46 1.3992E-3 1.2961E-3 0.9263 266.9 1.313E+2 1.714E-3 327.5 16.85 9.83E-3;...
48 1.0748E-3 1.0095E-3 0.9393 270.6 1.023E+2 1.317E-3 329.8 17.04 1.29E-2;...
50 8.3819E-4 7.8728E-4 0.9393 270.6 7.977E+1 1.027E-3 329.8 17.04 1.66E-2;...
52 6.5759E-4 6.1395E-4 0.9336 269.0 6.221E+1 8.055E-4 328.8 16.96 2.10E-2;...
54 5.2158E-4 4.7700E-4 0.9145 263.5 4.833E+1 6.389E-4 325.4 16.68 2.61E-2;...
56 4.1175E-4 3.6869E-4 0.8954 258.0 3.736E+1 5.044E-4 322.0 16.40 3.25E-2;...
58 3.2344E-4 2.8344E-4 0.8763 252.5 2.872E+1 3.962E-4 318.6 16.12 4.07E-2;...
60 2.5276E-4 2.1668E-4 0.8573 247.0 2.196E+1 3.096E-4 315.1 15.84 5.11E-2;...
62 1.9647E-4 1.6468E-4 0.8382 241.5 1.669E+1 2.407E-4 311.5 15.55 6.46E-2;...
64 1.5185E-4 1.2439E-4 0.8191 236.0 1.260E+1 1.860E-4 308.0 15.26 8.20E-2;...
66 1.1668E-4 9.3354E-5 0.8001 230.5 9.459E+0 1.429E-4 304.4 14.97 1.05E-1;...
68 8.9101E-5 6.9593E-5 0.7811 225.1 7.051E+0 1.091E-4 300.7 14.67 1.34E-1;...
70 6.7601E-5 5.1515E-5 0.7620 219.6 5.220E+0 8.281E-5 297.1 14.38 1.74E-1;...
72 5.0905E-5 3.7852E-5 0.7436 214.3 3.835E+0 6.236E-5 293.4 14.08 2.26E-1;...
74 3.7856E-5 2.7635E-5 0.7300 210.3 2.800E+0 4.637E-5 290.7 13.87 2.99E-1;...
76 2.8001E-5 2.0061E-5 0.7164 206.4 2.033E+0 3.430E-5 288.0 13.65 3.98E-1;...
78 2.0597E-5 1.4477E-5 0.7029 202.5 1.467E+0 2.523E-5 285.3 13.43 5.32E-1;...
80 1.5063E-5 1.0384E-5 0.6893 198.6 1.052E+0 1.845E-5 282.5 13.21 7.16E-1;...
82 1.0950E-5 7.4002E-6 0.6758 194.7 7.498E-1 1.341E-5 279.7 12.98 9.68E-1;...
84 7.9106E-6 5.2391E-6 0.6623 190.8 5.308E-1 9.690E-6 276.9 12.76 1.32E+0;...
86 5.6777E-6 3.6835E-6 0.6488 186.9 3.732E-1 6.955E-6 274.1 12.53 1.80E+0];
k = h/1000;
T=interp1(satm(:,1),satm(:,5), k,'makima', 0);
P=interp1(satm(:,1),satm(:,6), k,'makima');
rho=interp1(satm(:,1),satm(:,7), k,'makima');
end
function [L, D] = Lift_Drag_Components(h, alpha, v, rho, Sref)
cL = get_cl(h, v, alpha);
cD = get_cd(h, v, alpha);
L=0.5*rho*v^2*Sref*cL;
D=0.5*rho*v^2*Sref*cD;
end
function cD = get_cd(h, v, alpha)
%{
INPUTS
h -- altitude (m)
v -- velocity (m/s)
alpha -- angle of attack (rad)
OUTPUT
cD -- Interpolated coefficient of drag
%}
% Computation of Mach number
k=1.4;
R=287.058; % J/(kg*K)
[T, P, rho] = standard_atm(h); % atmospheric parameters, including T temperature
c=(k*R*T)^(1/2); % speed of sound, m/s
Ma=v/c;
% Loading cdfile.txt as a matrix
cdfile =[-4 -2 0 2 4 6 8 10 12 14 16 18 20 22.5 25 30 35 40;...
0.3 0.0385 0.0325 0.0315 0.0365 0.0481 0.0665 0.0932 0.1293 0.1776 0.2422 0.4052 0.7744 1.3828 1.7971 2.8139 4.0548 5.468;...
0.6 0.0376 0.0325 0.0319 0.0365 0.0475 0.0656 0.0911 0.1256 0.173 0.2363 0.39 0.7176 1.3207 1.7197 2.7047 3.9168 5.3035;...
0.9 0.0401 0.0364 0.0364 0.0404 0.0503 0.0668 0.0894 0.1203 0.1636 0.2213 0.3524 0.6182 1.1713 1.5304 2.4271 3.5466 4.8951;...
1.2 0.1192 0.1173 0.1192 0.1262 0.1379 0.1555 0.1696 0.1872 0.2068 0.2296 0.2551 0.2862 0.3306 0.3775 0.4809 0.6001 0.7246;...
1.5 0.1212 0.119 0.1186 0.1222 0.129 0.14 0.1552 0.1743 0.1949 0.2179 0.2431 0.2741 0.3189 0.3666 0.4722 0.5936 0.7201;...
2 0.1104 0.1064 0.105 0.1073 0.1131 0.1225 0.1357 0.1533 0.1729 0.1956 0.2212 0.2523 0.2968 0.3438 0.4476 0.5671 0.6914;...
3 0.0946 0.0893 0.0877 0.09 0.0961 0.1059 0.1193 0.1367 0.156 0.1787 0.2041 0.2348 0.2783 0.3241 0.4253 0.5423 0.6641;...
5 0.0814 0.0754 0.0741 0.0773 0.0849 0.0958 0.1096 0.1269 0.146 0.1683 0.1931 0.2233 0.2659 0.3107 0.4093 0.5225 0.6398;...
7.5 0.0761 0.07 0.069 0.0732 0.0813 0.0919 0.1055 0.1225 0.1412 0.1634 0.1879 0.2173 0.2589 0.3026 0.3993 0.5117 0.6288;...
10 0.074 0.0679 0.0671 0.0718 0.0798 0.0904 0.1038 0.1206 0.1394 0.1614 0.1857 0.2149 0.256 0.2992 0.3949 0.5068 0.6236;...
15 0.0726 0.0664 0.0657 0.0708 0.0786 0.0892 0.1025 0.1192 0.1379 0.1597 0.1838 0.2127 0.2535 0.2964 0.3917 0.5034 0.6201;...
20 0.0719 0.0655 0.0651 0.0704 0.078 0.0885 0.1019 0.1185 0.137 0.1586 0.1825 0.2112 0.2518 0.2945 0.3894 0.501 0.6178;...
40 0.0719 0.0655 0.0651 0.0704 0.078 0.0885 0.1019 0.1185 0.137 0.1586 0.1825 0.2112 0.2518 0.2945 0.3894 0.501 0.6178];...
% the list of alpha (angle of attack) values are in the first row
% convert from deg -> rad
alpha_list = cdfile(1,2:end).*(pi/180);
% the list of Mach values are in the first column, in both cases, you need to skip
% the first entry (NaN)
Ma_list = cdfile(2:end,1);
cD_matrix = cdfile(2:end, 2:end);
% look up griddata in the help file of matlab to learn how it works
cD = griddata(alpha_list(:), Ma_list(:), cD_matrix, alpha, Ma);
end
function cL = get_cl(h, v, alpha)
%{
INTERPOL_CL loads the lift coefficient cL dataset, and interpolates a
single value of cL for given discrete values of altitude(h), velocity (v)
and angle of attack (alpha)
INPUTS
h -- altitude (m)
v -- velocity (m/s)
alpha -- angle of attack (rad)
OUTPUT
cL -- Interpolated coefficient of lift
%}
% Computation of Mach number
k=1.4;
R=287.058; % J/(kg*K)
[T, P, rho] = standard_atm(h); % atmospheric parameters, including T temperature
c=(k*R*T)^(1/2); % speed of sound, m/s
Ma=v/c;
% You can load in a text file, or copy & paste the data as a matrix in here
clfile = [-4 -2 0 2 4 6 8 10 12 14 16 18 20 22.5 25 30 35 40;...
0.3 -0.0982 -0.0222 0.0539 0.1299 0.2059 0.2854 0.3719 0.4649 0.564 0.6691 0.7787 0.8905 1.0336 1.1843 1.4933 1.8025 2.1016;...
0.6 -0.0971 -0.0209 0.0552 0.1313 0.2074 0.287 0.3735 0.4666 0.5658 0.6708 0.7803 0.8928 1.0366 1.188 1.4993 1.8119 2.1143;...
0.9 -0.0922 -0.017 0.0582 0.1333 0.2085 0.2872 0.3729 0.4652 0.5638 0.6684 0.7779 0.8912 1.0367 1.1898 1.5058 1.825 2.1465;...
1.2 -0.0711 -0.0093 0.0504 0.114 0.1856 0.2566 0.3014 0.3403 0.3787 0.4177 0.4578 0.4943 0.5351 0.5746 0.6437 0.693 0.7352;...
1.5 -0.0715 -0.024 0.0219 0.0698 0.1226 0.1773 0.2318 0.2835 0.3332 0.3779 0.4178 0.4541 0.4981 0.5413 0.6178 0.6712 0.7166;...
2 -0.0618 -0.0217 0.0176 0.0581 0.1016 0.1463 0.191 0.2368 0.282 0.3235 0.3619 0.3976 0.4411 0.4841 0.5616 0.6193 0.6702;...
3 -0.0549 -0.0201 0.014 0.0496 0.0864 0.1244 0.1627 0.2018 0.2408 0.2784 0.3148 0.3498 0.3928 0.4357 0.5141 0.575 0.6299;...
5 -0.0498 -0.0192 0.0116 0.0432 0.0761 0.1096 0.1435 0.1777 0.212 0.2466 0.2816 0.316 0.3584 0.4009 0.4793 0.5418 0.5984;...
7.5 -0.0471 -0.0187 0.0104 0.0403 0.0712 0.1026 0.135 0.1683 0.2017 0.2358 0.2706 0.3047 0.3469 0.3892 0.4677 0.5308 0.5883;...
10 -0.0457 -0.0185 0.0098 0.039 0.0686 0.0994 0.1314 0.1645 0.1981 0.2322 0.267 0.3011 0.3433 0.3856 0.4641 0.5273 0.5848;...
15 -0.0448 -0.0184 0.0092 0.0379 0.0666 0.0973 0.1293 0.1623 0.1958 0.2298 0.2644 0.2984 0.3405 0.3829 0.4615 0.5248 0.5825;...
20 -0.0441 -0.0183 0.0087 0.037 0.0651 0.0956 0.1279 0.1608 0.1941 0.2278 0.2622 0.2961 0.3383 0.3806 0.4593 0.5228 0.5807;...
40 -0.0441 -0.0183 0.0087 0.037 0.0651 0.0956 0.1279 0.1608 0.1941 0.2278 0.2622 0.2961 0.3383 0.3806 0.4593 0.5228 0.5807];...
% the list of alpha (angle of attack) values are in the first row
% convert from deg -> rad
alpha_list = clfile(1,2:end).*(pi/180);
% the list of Mach values are in the first column, in both cases, you need to skip
% the first entry (NaN)
Ma_list = clfile(2:end,1);
cL_matrix = clfile(2:end, 2:end);
% look up griddata in the help file of matlab to learn how it works
cL = griddata(alpha_list(:), Ma_list(:), cL_matrix, alpha, Ma);
end
  2 comentarios
Walter Roberson
Walter Roberson el 29 de Nov. de 2019
heat_flux is missing.
Michael Fitzpatrick
Michael Fitzpatrick el 29 de Nov. de 2019
function q = heat_flux(rho, v, Rn)
kh = 1.83e-4 %
q = kh*sqrt(rho/Rn)*v^3

Iniciar sesión para comentar.

Iniciar sesión para responder a esta pregunta.

Respuestas (0)

Categorías

Más información sobre Aerospace Applications en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by