Gravity turn solver only works for a gamma of 90 degrees

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 :)

Answer 1

[William Rose]

@Jordan Booth

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.

[William Rose]

@Jordan Booth,

The calculation of has the same issue: multiply the right hand side of dpdt(4) by 180/pi to get in degrees and therefore to get gamma in degrees.

Or have both theta and gamma and their derivatives in radians. That would be my personal preference. In that case you would replace cosd() and sind() with cos() and sin(). Then your eqations for and in function dpdt are correct without units conversions.

[Jordan Booth]

dpdt = [p(3)*sind(p(4)); % r
        (180/pi)*(p(3)*cosd(p(4)))/p(1); % theta in degrees
        T/m-A*Cd*aD*p(3)^2/(2*m)-u*sind(p(4))/p(1)^2; % velocity
        (180/pi)*((-u*cosd(p(4)))/(p(3)*p(1)^2)+(p(3)*cosd(p(4)))/p(1))]; % gamma in degrees

I have tried the above and tried using sin/cos instead of sind/cosd but still cant get the result I'm looking for. I understand if you are too busy but have you tried to get it to run on your end? It makes sense that theta and gamma need to be converted to degrees but dont know how to implement it correctly.

[William Rose]

@Jordan Booth, I have not tried it. I will.

[William Rose]

@Jordan Booth,

[edited: attached code with a fix]

A modified version of your code is attached. Angle calculations are in degrees. I have modified the figures a bit. When gamma0=90, it performs as you have already seen: rocket goes straight up, then down.

When gamma=89.999, the x,y figure seems to show practically no motion (left hand figure below). Zooming in (right figure) shows that the rocket goes "straight up" (89.999 degrees), reaches a max altitude of 2.5 mm, then impacts the Earth surface about 16 mm downrange. (I changed rE to 6.371e6, exactly, to make it easier to interpret the x,y plot.)

The result above suggests to me that something is wrong with the equation for d(gamma)/dt. I recommend that you check the equation and the parameters. I see an interesting Stack Exchange discussion here. The notation at S.E. differs somewhat from yours, and the equations there includes terms not in your model. These terms are lift (L), which is zero in your model, and , the angle between the thrust vector and the velocity vector, which is assumed to be zero in your model. Those simplifications in your model are very reasonable, in my opinion.

[William Rose]

@Jordan Booth,

Zoom in on the plots of h, v, gamma versus time as shown below, showing the first 0.1 seconds.

The plot of gamma versus t shows that the rocket tips over very quickly: by t=0.05 s, the rocket is flying horizontally (gamma=0), and after that time it is pointing down (gamma<0). The plot of height versus time shows that it hits the ground (h=0) around t=0.065 s.

I wrote a version of the code in which all angles are in radians instead of degrees. The results are the same when the initial angle is the radian equivalent of 89.999 degrees. This gives me confidence that the surprising results are not due to a radians versus degrees conversion issue.

Consider doing the following to figure out what is going on: Turn the thrust down to zero. This should cause the mass to be constant, since no fuel is leaving. Make the initial velocity large. Make the air density zero, to eliminate drag. Then the trajectory should be a simple parabola. Is the trajectory what you expect? This could lead to some insight.

[Jordan Booth]

After setting aD and T to 0, I set the initial velocity to 3000 m/s and gamma to 85 degrees. The result is a nice parabolic curve however the maximum height is 489706 m whereas following standard projectile motion equations you get 455387 m. Im assuming this is due to the ode solver not being perfect although I dont have enough experience to know if this is the case.

With an initial velocity of 50, gamma 89.3 and 0kg of payload you get:

At first glance it looks promising but at 700km and 12km/s it should be an elliptical orbit not suborbital. I am thinking of modeling the rocket going vertical until it gains speed then pitching over but I dont think that will solve the problem.

[William Rose]

@Jordan Booth,

I suspect the difference in max altitude (490 km with the ODEs, 455 km from standard parabolic equation) is due to the decrease in gravitational force with altitude, which is accounted for in the ODEs but is not included in the parabolic equation.

In the two figures you posted above, with v0=50 m/s and gamma0=89.3, are the thrust and air density back to their normal values? It appears that they are. Why did you make v0=50 m/s?

[Jordan Booth]

Yes thrust and air density where set back to their normal values. I was seeing what the response was to a higher initial velocity because it worked for the test without thrust and drag. Im certain the thrust to weight is greater than 1 so there shouldnt be a problem with using a low initial velocity which means I'm missing something.

I have tried to add the following into the ode function:

flag = 0;
if t >= 10 && flag == 0
    p(4) = 89.5;
    flag = 1;
end

Using the initial conditions:

p0 = [rE; 0; 0.001; 90];

With the intention of building velocity before pitching over but Im still trying to get that to work.

[William Rose]

@Jordan Booth,

It is the rapid change in gamma at the start, when v is approximately zero, that results in the unexpected rapid return to Earth of the rocket. This makes me wonder if the differential equation for gamma is unreliable when v~=0. Note that term 1 in the differential equation for gamma includes v in the denominator, so approaches infinity as v approaches zero.

Therefore I rewrote the differential equations in Cartesian coordiantes. In this version of the model, the state vector is p = [x, y, vx, vy]. This version is called rocketLaunch2stagexy.m (attached).

>> rocketLaunch2stage  %polar version
Max altitude=0.0028 m.
Earth impact time=0.071 s.
>> rocketLaunch2stagexy  %Cartesian version
Max altitude=0.0029 m.
Earth impact time=0.071 s.

The results from the polar coordinates model and the Cartesian coordinates model are are basically the same, as shown in the plots below and in the console output above. In other words, the change to Cartesian did not fix the rapid turnover at the start. In both models shown below, v0=1e-3 and gamma0=89.999. I think this shows us that a rocket is unstable.

left: polar coordinates model; right: Cartesian. v0=1e-3, gamma0=89.999 in both

Why does this problem not happen in real rockets? Due to moment of inertia (not included in our model) and because large rockets have closed loop attitude control: they adjust the direction of the thrust vector dynamically to acheive a desired trajectory and orientation. Small model rockets do not have dynamic control, and they do not launch at exactly 90, and yet they do not immediately turn over. Why not? Because of rapid acceleration, fins on the tail, and a rigid launch guide wire about a meter long. The guide wire keeps the rocket oriented, i.e. prevents tip-over, for the first meter of flight. This distance is enough for a model rocket to reach v>7 m/s, fast enough for the fins to provide stability. (Model rocket acceleration=27-32 m/s^2, different web sites.)

[Jordan Booth]

Thank you for your continued help!

After checking some parameters A should be equal to 10.75 m^2 for the chosen rocket although this doesnt change the issues in the model.

p0 = [rE+5e5; 0; 7616; 0];

starting at 500km and 7616 m/s should be a circular orbit but you get this:

Do you think this is due to the same issue of stability and if I found a value for the moment of inertia is this something that could be implemented?

[William Rose]

What is gamma0 for the launch above?

I think the cure is to use thrust vectoring, i.e. thrust will have to be allowed to be something other than straight ahead. Among the possible approaches are :

• Implement a schedule where the thrust vector angle changes to specified vaues at specified times.
• Implement closed loop control. This is more complicated, and more like a real large rocket. For this, you must have a target, which could be a profile of downrange distance versus time, or gamma versus time, etc. Then you adjust the thrust angle in the direction that will help correct for the current deviation. This can cause instability if the gain of the control system is too high.

[William Rose]

@Jordan Booth, I htink there is still an error in the diffrential equaiton for gamma. The code says

dpdt = [p(3)*sind(p(4));                            % r
    (180/pi)*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
    -(180/pi)*u*cosd(p(4))/(p(3)*p(1)^2)+p(3)*cosd(p(4))/p(1)]; % gamma

but it should be

dpdt = [p(3)*sind(p(4));                            % r
    (180/pi)*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
    -(180/pi)*(u*cosd(p(4))/(p(3)*p(1)^2)+p(3)*cosd(p(4))/p(1))]; % gamma

where I have added a set of (...) in the gamma equation. This is a radians-to-degrees error, again. I thought I had found all the radians-to-degrees errors, but I guess not. It is yet another example of why I do not like to do calculations in degrees: too many chances for mistakes like this. I suggest you make this correction and see if it makes a difference.

[William Rose]

@Jordan Booth,

Please send me email by clicking on the envolope at the upper right of the pop-up window that appears when you click on my WR icon.

[Jordan Booth]

I have added the expected solution in green

gamma0 = 0;
p0 = [rE+10e5; 0; 8000; gamma0]; % initial r, theta, velocity, gamma
[t,p] = ode45(@ode,[0 tfinal],p0); % ode solver

u = 3.986004418e14; % standard gravitational parameter
SemiMajor = (u*p0(1))/(2*u-p0(1)*p0(3)^2); % semi major axis
foci = 2*(SemiMajor-p0(1)); % distance between foci
semiMinor = sqrt(-(foci^2-SemiMajor^2)); % semi minor axis
C = [-(SemiMajor-p0(1)) 0]; % center of orbit
th = linspace(0,2*pi);
xe = C(1)+SemiMajor*cos(th);
ye = C(2)+semiMinor*sin(th);

[William Rose]

@Jordan Booth,

I found another mistake in the equation for dgamma/dt: an incorrect sign for the second term. WIth this error fixed, we get the expected result for a zero-thrust flight starting at alt=500 km, v=7616: a circle with an orbit time of slightly over 90 minutes. See plots, console output, and corrected script.

>> rocketLaunch2stage
Max and final altitudes=500000.0000, 499932.13 m.
No Earth impact.

This looks good. This adjustment to the differential equaiton for gamma does not change the result for a launch with normal thrust, from altitude=0, with v0=1e-3. In that case, the rocket still tips over very quickly, and impacts 16-17 mm downrange after 0.07 seconds. The reason that the change to the code did not affect the solution for such a launch is that the term whose signed was changed is insignificantly small for such a launch. Therefore I think you will need to add a thrust vector angle to the model. The current model assumes the thrust vector is 0 (ie. aligned with the velocity vector) at all times.

[Jordan Booth]

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.

[William Rose]

@Jordan Booth, well done! You're welcome.

[Rohit]

@Jordan Booth Sir can you please attach updated script i'm also looking solution for the same.