Code not working, velocity comes back the same each time
Show older comments
ive finally got to the point where my code actually runs without bringning back error messages, however now whatever values i plug in my answer comes out the same and the answer seems way too high for the impact velocity of a spacecraft hitting earth. if anyone could please correct this that would be great, thank you :) my code is below
% M = 850;
% As = 5;
% Ap = 150;
% Hs = 150000;
% Hp = 3000;
% U = 0;
promptA = 'What is the mass of the Spacecraftt in Kg? ';
M = input(promptA);
promptB = 'What is the Cross Sectional Area of the Spacecraft in m^2? ';
As = input(promptB);
promptC = 'What is the Cross Sectional Area of the Parachute in m^2? ';
Ap = input(promptC);
promptD = 'What is the Initial Height of the Spacecraft in m? ';
Hs = input(promptD);
promptE = 'What is the parachute deployment height in m? ';
Hp = input(promptE);
promptF = 'What is the initial velocity of the Spacecraft in m/s? ';
U = input(promptF);
%Constants
C=0.7; % Drag Coefficient
Ha=100000; % Height of the Atmosphere in m
Tstep=1; % Time step in seconds
% Initialisation of Increments
i=1;
H(i)=Hs;
V(i)=0;
S(i)=0;
t(i)=0;
% Freefall above the atmosphere
while (i <= length(H) && H(i) > Ha)
Tstep=1;
Fd=0;
g(i)=(40*(10^7))/(6371+H(i));
a(i)=g(i);
V(i+1)=V(i)+(a(i)*Tstep); % Velocity at next i
S(i+1)=S(i)+(V(i)*Tstep+0.5*a(i)*Tstep^2); % Displacement at next i
t(i+1)=t(i)+Tstep; % Time at next incriment
i=i+1;
continue
end
% Freefall within the Atmosphere (No parachute)
while (i <= length(H) && H(i) > Hp )
Tstep=1;
p(i)=(H(i)/71)+1.4; % Drag Coefficient for Fd
Fd(i)=(0.5*C*p(i)*As*(V(i)^2)); % Air Resistance
a(i)=(g(i)-(Fd(i)/M)); % New Acceleration formula inc. Fd
V(i+1)=V(i)+(a(i)*Tstep);
S(i+1)=S(i)+(V(i)*Tstep+0.5*a(i)*Tstep^2);
t(i+1)=t(i)+Tstep;
i=i+1;
continue
end
% Freefall within the Atmosphere with a parachute
while (i <= length(H) && H(i)> 0)
Tstep=1;
p(i)=(H(i)/71)+1.4; % Drag Coefficient for Fd
Fd(i)=(0.5*C*p(i)*Ap*(V(i)^2)); % Air Resistance
a(i)=(g(i)-(Fd(i)/M)); % New Acceleration formula inc. Fd
V(i+1)=V(i)+(a(i)*Tstep);
S(i+1)=S(i)+(V(i)*Tstep+0.5*a(i)*Tstep^2);
t(i+1)=t(i)+Tstep;
i=i+1;
continue
end
display (V(i))
6 Comments
darova
on 8 Apr 2020
Please edit your code according following instructions
No one will input those values
Oliver Stericker-Taylor
on 8 Apr 2020
Edited: Oliver Stericker-Taylor
on 8 Apr 2020
darova
on 8 Apr 2020
Yes. It will not change anything. It's for me
Oliver Stericker-Taylor
on 8 Apr 2020
Oliver Stericker-Taylor
on 8 Apr 2020
darova
on 8 Apr 2020
Thank you. Very generous of you!
Answers (3)
James Tursa
on 8 Apr 2020
Edited: James Tursa
on 8 Apr 2020
Lines like this are wrong:
S(i+1)=S(i)+(V(i)*Tstep+0.5*a(i)*Tstep^2); % Displacement at next i
You are double booking the effect of gravity. Once for the V(i)*Tstep (which already includes the gravity effect) and then again for the 0.5*a(i)*Tstep^2. You only need the V(i) part in all three of your loops:
S(i+1)=S(i)+V(i)*Tstep; % Displacement at next i
Also, you never update H(i+1) in your loops. You need to update the height.
Gravity needs to be negative since it is an attractive force. So
g(i) = -(40*(10^7))/(6371+H(i));
And this g(i) calculation needs to be part of every loop. Gravity doesn't stop acting just because you are in the atmosphere.
Finally, I don't see the purpose of testing i <= length(H). Why would you exit the loop if i is greater than length(H)? Get rid of that part of the test and only test on height. If you want to limit the number of iterations, then put that test within the loop and generate an error if it fails. E.g., you could put something like this in each loop
if( i > maxi )
error('Too many iterations, i > %d',maxi);
end
with an appropriate value of maxi set earlier in your code.
15 Comments
Oliver Stericker-Taylor
on 8 Apr 2020
Edited: Oliver Stericker-Taylor
on 8 Apr 2020
James Tursa
on 8 Apr 2020
H(i+1) = H(i) + V(i)*Tstep;
James Tursa
on 8 Apr 2020
Also, the drag force always acts opposite of the velocity direction. I prefer to use abs(v)*v instead of v^2 so that the equation always works regardless of the sign of v. E.g., change this
Fd(i)=(0.5*C*p(i)*As*(V(i)^2)); % Air Resistance
a(i)=(g(i)-(Fd(i)/M));
to this
Fd(i) = -0.5 * C * p(i) * As * abs(V(i))*V(i); % Air Resistance
a(i) = g(i) + Fd(i)/M;
Same for the last loop.
Oliver Stericker-Taylor
on 8 Apr 2020
Oliver Stericker-Taylor
on 8 Apr 2020
James Tursa
on 8 Apr 2020
Well, should the height change inside of each loop? Yes, of course. So this line needs to go inside each loop.
Oliver Stericker-Taylor
on 8 Apr 2020
James Tursa
on 8 Apr 2020
Please post your current code.
Oliver Stericker-Taylor
on 8 Apr 2020
James Tursa
on 8 Apr 2020
Edited: James Tursa
on 8 Apr 2020
As I stated before, get rid of the i <= length(H) part of the test. You should have
while (H(i) > Ha)
:
while (H(i) > Hp)
:
while (H(i)> 0)
And you have too many lines here:
a(i) = g(i) + Fd(i)/M;
a(i)=(g(i)-(Fd(i)/M)); % New Acceleration formula inc. Fd
Get rid of that second a(i) = etc. line.
And in your parachute loop:
Fd(i) = -0.5 * C * p(i) * As * abs(V(i))*V(i); % Air Resistance
I think that As should be Ap.
Oliver Stericker-Taylor
on 8 Apr 2020
James Tursa
on 8 Apr 2020
Edited: James Tursa
on 8 Apr 2020
You are missing your p(i) code in the last two loops. Also look at the fixes suggested by David.
Oliver Stericker-Taylor
on 9 Apr 2020
James Tursa
on 9 Apr 2020
I'm still looking at this. There's issues with the gravity and drag forces. I'll let you know what I find out later today ...
Oliver Stericker-Taylor
on 9 Apr 2020
David Goodmanson
on 8 Apr 2020
Hello Oliver,
could you say where you got this problem? I am interesed in finding out since the problem keeps reoccurring on this website. See the link below, where I believe the calculation to be numerically correct.
This is a very badly posed problem, not in terms of what is being asked but in terms of units and dimensions. The expression for g should be
g(i)=(40*(10^7))/(6371+H(i))^2 not g(i)=(40*(10^7))/(6371+H(i)) as you have it.
However, in this expression, g is in m/sec^2 but H must be in km. This is not a good idea, rife with the possibility for error.
The expression for p should be
p(i) = -(H(i)/71) +1.4 not p(i) = (H(i)/71) +1.4 as you have it
since I think you'll agree that air density increasing with height doesn't make sense. Here p is in kg/m^3 as it should be, but again H must be in km. [ The expression is ludicrously inaccurate as a model for air density vs. height in the atmosphere. But for problem purposes it's all right as long as people are aware of it. ]
In the link below, you can see how the signs are worked out and the height checks are done. But that calculation is done with velocity in m's, acceleration in m/sec^2 and height in km. You will see a couple of factors of .001 to make the conversion. That approach is tricky, and it's really better to go with meters everywhere, in which case H is in meters and
g(i) = (40e7)/(6371 +H(i)/1e3)^2 or equivalently (40e13)/(6371e3 +H(i))^2
and
p(i) = -(H(i)/71e3) +1.4
Then just go with
V(i+1) = V(i) +a(i)*Tstep;
S(i+1) = S(i) +V(i)*Tstep
at each stage, no conversion factors.
2 Comments
Oliver Stericker-Taylor
on 8 Apr 2020
Edited: Oliver Stericker-Taylor
on 8 Apr 2020
brilliant thank you ill have a go at changing these values over. im a first year mechanichal engineer and this is our maths coursework, ill show you the questions we have as id imagine were all doing the exact same ones. weve been asked to do this with no prior knowledge of matlab and we have no resources to use to learn as were all stuck at home so im guessing everyones struggling with it. im gonna retry with these new values now although i think ill probably put them in wrong again haha, thank you for your help
Oliver Stericker-Taylor
on 9 Apr 2020
Edited: Oliver Stericker-Taylor
on 9 Apr 2020
James Tursa
on 9 Apr 2020
Edited: James Tursa
on 10 Apr 2020
So here is a version of your code with the following changes:
Gravity model replaced with the model found here (also note that the sign of gravity in the code is negative):
Atmosphere density model replaced with the simple isothermal model found here (using the exp(etc) model instead of truncated Taylor series):
Time step dropped significantly when the parachute opens to account for the large spike in acceleration due to the high velocity. This turns out to be extremely important to keep your simulation stable.
So, only a few changes made to your code, but extremely important ones. With these slight changes, the code seems to work OK and you should be able to manipulate it to answer your homework questions. At least now you can work on getting answers to the specific questions instead of spinning your wheels just getting things to work.
The fundamental problems I found in your latest code were:
1) Incorrect sign on the gravity
2) Incorrect density formula (that 71e3 simply isn't correct even for the truncated Taylor Series version)
3) Way too big of a stepsize, particularly when the parachute opens, to keep the simulation stable.
clear all
M = 850; % mass (kg)
As = 5; % s/c area (m^2)
Ap = 150; % parachute area (m^2)
Hs = 150 * 1000; % initial height (m)
Hp = 3 * 1000; % height at parachute open (m)
U = 0; % initial velocity (m/s)
%Constants
C = 0.7; % Drag Coefficient (dimensionless)
Ha = 100 * 1000; % Height of the Atmosphere (m)
Tstep = 1; % Time step in seconds
% Earth constants
g0 = 9.80665; % standard gravity (m/s^2)
Re = 6371 * 1000; % Earth radius (m)
% Initialisation of Increments
i = 1;
H(i) = Hs;
V(i) = U;
t(i) = 0;
% All units below are (m,kg,s)
% Freefall above the atmosphere
disp('Freefall outside atmosphere')
while (H(i) > Ha)
Fd(i) = 0;
% g(i) = -(40e7)/(6371+ H(i)/1e3)^2;
g(i) = -g0 * (Re / (Re + H(i)))^2;
a(i) = g(i) + Fd(i)/M;
V(i+1) = V(i) + a(i)*Tstep; % Update velocity
H(i+1) = H(i) + V(i)*Tstep; % Update height
t(i+1) = t(i) + Tstep; % Update time
i = i + 1;
continue
end
% Freefall within the Atmosphere (No parachute)
disp('Inside atmosphere')
while (H(i) > Hp)
% p(i) = -(H(i)/71e3)+1.4; % Drag Coefficient for Fd
p(i) = 1.3 * exp(-H(i)/7000);
Fd(i) = -0.5 * C * p(i) * As * abs(V(i))*V(i); % Air Resistance
% g(i) = -(40e7)/(6371+ H(i)/1e3)^2;
g(i) = -g0 * (Re / (Re + H(i)))^2;
a(i) = g(i) + Fd(i)/M;
V(i+1) = V(i) + a(i)*Tstep; % Update velocity
H(i+1) = H(i) + V(i)*Tstep; % Update height
t(i+1) = t(i) + Tstep; % Update time
i = i + 1;
continue
end
% Freefall within the Atmosphere with a parachute
disp('deploy parachute, changing stepsize');
Tstep = 0.05;
while (H(i)> 0)
% p(i) = -(H(i)/71e3)+1.4; % Drag Coefficient for Fd
p(i) = 1.3 * exp(-H(i)/7000);
Fd(i) = -0.5 * C * p(i) * Ap * abs(V(i))*V(i); % Air Resistance
% g(i) = -(40e7)/(6371+ H(i)/1e3)^2;
g(i) = -g0 * (Re / (Re + H(i)))^2;
a(i) = g(i) + Fd(i)/M;
V(i+1) = V(i) + a(i)*Tstep; % Update velocity
H(i+1) = H(i) + V(i)*Tstep; % Update height
t(i+1) = t(i) + Tstep; % Update time
i = i + 1;
continue
end
disp(V(i))
n = numel(t)-1;
m = 1;
figure;
plot(t(m:n),H(m:n),'.-');
grid on
title('Height')
xlabel('Time (s)');
ylabel('Height (m)');
figure;
plot(t(m:n),V(m:n),'.-');
grid on
title('Velocity')
xlabel('Time (s)');
ylabel('Velocity (m/s)');
figure;
plot(t(m:n),a(m:n),'.-');
grid on
title('Acceleration')
xlabel('Time (s)');
ylabel('Acc (m/s^2)');
3 Comments
Oliver Stericker-Taylor
on 9 Apr 2020
Arqam Zaheer
on 28 Aug 2020
Hi Oliver Stericker-Taylor. Did you find the solution to this problem?
Can you please share it?
Regards
James Tursa
on 28 Aug 2020
A working code is already posted and shared. Can you tell us what more you need?
Categories
Find more on General Applications in Help Center and File Exchange
on 8 Apr 2020
on 28 Aug 2020
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
