97 views (last 30 days)

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');

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.

Cris LaPierre
on 21 Mar 2020

Edited: Cris LaPierre
on 21 Mar 2020

James Tursa
on 31 Aug 2020

Edited: James Tursa
on 31 Aug 2020

See this link for an example of the parachute problem:

Opportunities for recent engineering grads.

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

Start Hunting!
## 5 Comments

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_809095

⋮## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_809095

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_809835

⋮## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_809835

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_813254

⋮## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_813254

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_817393

⋮## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_817393

## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_988976

⋮## Direct link to this comment

https://se.mathworks.com/matlabcentral/answers/385016-the-spacecraft-free-fall-math-model#comment_988976

Sign in to comment.