Gravity turn solver only works for a gamma of 90 degrees

52 visualizaciones (últimos 30 días)
I am trying to plot a gravtiy turn trajectory for a rocket given an inital pitch over angle using the following equations of motion:
When the intial gamma is set to 90 there will be no pitch over and the solution seems correct so the v_dot equation should be correct. As soon as some angle less than 90 is used it doesnt seem to work.
main:
clc, clear
rE = 6371008.8; % radius of earth in m
t0 = 0;
tfinal = 1500; % total time of model
% starting values
p0 = [rE; 0; 0.001; 90]; % r, theta, velocity, gamma
[t,p] = ode45(@ode,[t0 tfinal],p0); % ode solver
fig1 = figure(1); % velocity plot
plot(t,p(:,3))
xlabel('time (s)')
ylabel('velocity (m/s)')
fig2 = figure(2); % altitude plot
plot(t,p(:,1)-rE) % height above surface of earth
xlabel('time (s)')
ylabel('height (m)')
axis([0 1500 0 3e6])
fig3 = figure(3); % ground relative plot
x = p(:,1).*cosd(p(:,2)); % polar to cartesian conversion
y = p(:,1).*sind(p(:,2));
plot(x,y)
hold on;
viscircles([0 0],rE);
axis([0 10e6 0 10e6])
axis square
hold off
ode function:
function dpdt = ode(t,p)
rE = 6371008.8; % radius of earth in m
dSL = 1.225; % sea level density in kg/m^3
h0 = 10400; % scale alt in m
A = 43.0084; % cross-sectional area in m^2
Cd = 0.237; % drag co-efficent
u = 3.986004418e14; % standard gravitational parameter
gA = 9.80665; % gravtiational acceleration at sea level in m/s^2
pm = 20000; % payload mass
% second stage
smp = 92670; % propellant mass in kg
sIsp = 348; % specific impulse in secs
sF = 981000; % thrust
sm0 = smp+3900+pm; % total mass
stf = (smp*gA*sIsp)/sF; % second stage burn time in secs
% first stage
fmp = 395700; % propellant mass in kg
fIsp = 283; % specific impulse in secs
fF = 7607000; % thrust
fm0 = fmp+25600+sm0; % total mass
ftf = (fmp*gA*fIsp)/fF; % first stage burn time in secs
aD = dSL*exp(-(p(1)-rE)/h0); % air density function
T = thrust(t,ftf,fF,stf,sF); % thrust function
m = mass(t,ftf,stf,sm0,smp,gA,sIsp,fm0,fIsp,T); % mass function
dpdt = [p(3)*sind(p(4)); % r
(p(3)*cosd(p(4)))/p(1); % theta
T/m-A*Cd*aD*p(3)^2/(2*m)-u*sind(p(4))/p(1)^2; % velocity
(-u*cosd(p(4)))/(p(3)*p(1)^2)+(p(3)*cosd(p(4)))/p(1)]; % gamma
end
outputs for inital gamma of 90 degrees:
outputs for initial gamma of 89 degrees:
The relative graph should show some pitch over if the function is working correctly
thrust:
function T = thrust(t,ftf,fF,stf,sF)
if t < ftf
T = fF;
elseif t < ftf + 12
T = 0;
elseif t < stf + 12 + ftf
T = sF;
else
T = 0;
end
end
mass:
function m = mass(t,ftf,stf,sm0,smp,gA,sIsp,fm0,fIsp,T)
if t >= stf + 12 + ftf
m = sm0-smp;
elseif t > ftf + 12
m = sm0-((T*(t-ftf-12))/(gA*sIsp));
elseif t > ftf
m = sm0;
else
m = fm0-((T*t)/(gA*fIsp));
end
end
The thrust and mass functions perform as expected:
I have been trying for a while to figure it out and have looked online for solutions but cant find anything similar that works. I would be grateful for any help, thanks :)

Respuesta aceptada

William Rose
William Rose el 2 de Feb. de 2023
This looks like a nice model!
Your equation for = d(theta)/dt = dp(2)/dt is
dpdt = [...
(p(3)*cosd(p(4)))/p(1); % theta...
I think you want to do
dpdt = [...
(180/pi)*p(3)*cosd(p(4))/p(1); % theta in degrees...
because you compute x(t),y(t) with cosd(p(:,2)) and sind(p(:,2)), i.e. you assume theta is in degrees.
Alternatively, you could compute x and y with cos(p(:,2)) and sin(p(:,2)), but then you would have theta in radians and gamma in degrees, which would be kind of wierd.
  17 comentarios
Jordan Booth
Jordan Booth el 7 de Feb. de 2023
Thank you so much!!!! This is exactly what I was looking for.
With some tweaking I managed to get a stable orbit with a payload of 16000kg. I had to use an initial velocity of 50m/s which if I make the rocket fly vertical until it reaches this speed then pitch over it will negate the issue. I think this is a good estimate and will allow future models to use thrust vectoring to remain on the correct trajectory.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

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

Productos


Versión

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by