MATLAB Answers

The spacecraft free-fall math model

97 views (last 30 days)
Sanzhar Januzakov
Sanzhar Januzakov on 26 Feb 2018
Edited: James Tursa on 31 Aug 2020
Please, help to figure out my problem! This code very roughly describes the free-fall movement of spacecraft in 3 external conditions: in outer atmosphere, in the atmosphere and in the atmosphere with a parachute. I'm trying to find the minimum parachute deployment height at which I will obtain the minimum velocity of contact with the earth but I'm always getting the same velocity whether it was opened at 100km height or at 1km. Guessing, something wrong with air density or with velocity in 3rd while loop. Please, could you help with this? Thank you in advance.
% Input data
As = 5; %effective area of the spacecraft (m^2)
Ap = 150; %effective area of the parachute(m^2)
M = 850; %the spacecraft mass (kg)
Hs = 150; %the initial height (km)
Hp=3; %the parachute deployment height (km)
U=0; % initial velocity (m/s)
%Hp=3; % parachute deployment height (km)
Ha=100; % height of the atmosphere above the earth (km)
C=0.7; % drag coefficient
Tstep=0.01; % time step (s)
% Increment initialisation:
i=1; H(i)=Hs; v(i)=U; s(i)=0; t(i)=0; a(i)=(40*10^7)/(6371+H(i))^2;
% Freefall without atmosphere
while (H(i)>Ha)
FD(i)=0;
g(i)=(40*10^7)/(6371+H(i))^2; % Gravitational acceleration
a(i+1)=g(i); % Spacecraft acceleration
s(i+1)=s(i)+0.001*(v(i)*Tstep+0.5*a(i+1)*Tstep^2); % Displacement
v(i+1)=v(i)+a(i+1)*Tstep; % Velocity
t(i+1)=t(i)+Tstep; % Next time moment
H(i+1)=Hs-s(i+1); % Height
i=i+1;
continue
end
%Freefall in the atmosphere
A=As;
while (H(i)>Hp)
p(i)=1.4-H(i)/71; % Density
FD(i)=0.5*C*p(i)*A*v(i)*v(i); % Resistant force
g(i)=(40*10^7)/(6371+H(i))^2; % Gravitational acceleration
a(i+1)=g(i)-FD(i)/M; % Spacecraft acceleration
s(i+1)=s(i)+0.001*(v(i)*Tstep+0.5*a(i+1)*Tstep^2); % Displacement
v(i+1)=v(i)+a(i+1)*Tstep; % Velocity
t(i+1)=t(i)+Tstep; % Next time moment
H(i+1)=Hs-s(i+1); % Height
i=i+1;
continue
end
% %Freefall in the atmosphere with parachute
A=Ap;
while (H(i)>0)
p(i)=1.4-H(i)/71; % Density
FD(i)=0.5*C*p(i)*A*v(i)*v(i); % Resistant force
g(i)=(40*10^7)/(6371+H(i))^2; % Gravitational acceleration
a(i+1)=g(i)-FD(i)/M; % Spacecraft acceleration
s(i+1)=s(i)+0.001*(v(i)*Tstep+0.5*a(i+1)*Tstep^2); % Displacement
v(i+1)=v(i)+a(i+1)*Tstep; % Velocity
t(i+1)=t(i)+Tstep; % Next time moment
H(i+1)=Hs-s(i+1); % Height
i=i+1;
continue
end
% Display graphs
% In one Diagram
subplot(2,2,1); plot(t,s); % 1st graph - Displacement
grid; axis([0 1500 0 150]);
title('The dependence of displacement on time s(t), km'),
ylabel('s(t), km'),
xlabel('t, s')
subplot(2,2,2); plot(t,H); % 2nd graph - Height
grid; axis([0 1500 0 150]);
title('The dependence of height on time H(t), km'),
ylabel('H(t), km'),
xlabel('t, s')
subplot(2,2,3); plot(t,v); % 3rd graph - Velocity
grid; axis([0 1500 0 1000]);
title('The dependence of velocity on time v(t), m/s'),
ylabel('v(t) ), m/s'),
xlabel('t, s')
subplot(2,2,4); plot(t,a); % 4th graph - Acceleration
grid; axis([0 1500 -300 10]);
title('The dependence of acceleration on time a(t), m/s^2'),
ylabel('a(t), m/s^2'),
xlabel('t, s');

  5 Comments

Show 2 older comments
Luke Todd
Luke Todd on 21 Mar 2020
The equation for air density just comes from an equation sheet provided in the coursework brief
Vigneshwaran Sankar
Vigneshwaran Sankar on 29 Mar 2020
Hi @Sanzhar Januzakov:
Did u find the answer for this problem: "....I'm always getting the same velocity whether it was opened at 100km height or at 1km....".
I am also having same pbm right now?
Can u help me with this now?
Regards,
VS.
Arqam  Zaheer
Arqam Zaheer on 28 Aug 2020
@Sanzhar Januzakov Did you get this problem solved? Can you share your results?
Regards

Sign in to comment.

Answers (3)

David Goodmanson
David Goodmanson on 29 Mar 2020
Hi Sanzhar & Luke,
This problem keeps bouncing around on other posts, so it's time to check out what is going on.
Although the plotting has some issues, your answer is basically correct. The reason your ground-strike velocity never changes is that when the parachute is released, the capsule attains terminal velocity very quickly. So unless you release the parachute at a recklessly low elevation, the earth-strike velocity is always the same.
Temporarily neglecting the factors of (1/2) and C, and using density = rho instead of p, terminal velocity occurs when
M*g = rho*vterm^2*A, so vterm = sqrt(M*g/(rho*A))
The thing is, this occurs very quickly in terms of the capsule flight time. To find the characteristic time tau to get to terminal velocity, since g is the only factor around with units of time, tau is given by
g*tau = vterm so tau = sqrt(M/(g*rho*A))
which if you throw the neglected constants back in turns out to be about 1 sec with the larger parachute. Since tau is just a characteristic time and terminal velocity is never quite reached exactly, let's say 5 sec for practical purposes. That's still very much less than the flight time in the lower atmosphere, which is about 300 sec.
The code used units of meters for v and a, but units of kilometers for H. That's not a very good idea, but it seems to have been incorporated correctly.
As for the linear dropoff of air density with height, I don't think that
"an equation sheet provided in the coursework brief"
is good justification at all, whether you are a student or not. It might be all right if someone was pointing out that it is laughably unrealistic. Air density falls off roughly as exp(-H/H0) where the scale height H0 is about 8.5 km. That exponential curve is not even remotely fit with a straight line. Besides, 1.4- H/71 gives negative density at 100 km, which is hard to achieve. And it gives 1.4 at ground level, wheras the real number at ground level is 1.22.

  0 Comments

Sign in to comment.


Cris LaPierre
Cris LaPierre on 21 Mar 2020
Edited: Cris LaPierre on 21 Mar 2020
Another student in the class has asked several questions about this problem as well. Perhaps take a look at the comments in these posts to see if they answer your questions.
The most recent question can be found here.

  0 Comments

Sign in to comment.


Categories

Community Treasure Hunt

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

Start Hunting!