I've been trying to figure out how to build a rocket launch simulation that allows the user to specify the desired altitude at which a payload will be inserted into a stable circular orbit, but have been having a little trouble. It seems that a fairly straightforward, albeit outdated, approach is to use a so-called "linear tangent steering law" (http://www.dtic.mil/dtic/tr/fulltext/u2/643209.pdf) that controls the rocket's thrust angle and is given by:


Where $\theta_{0}$, $\theta_{f}$, $t$ and $t_{f}$ are the initial thrust pitch angle, final thrust pitch angle, current time and final time respectively. I applied the linear tangent steering law to the $\theta(t)$ terms in the system of differential equations below:

$$r'(t)=v(t)\sin\gamma(t),$$ $$\phi'(t)=\frac{v(t)}{r(t)}\cos\gamma(t),$$ $$v'(t)=\frac{F(t)}{m(t)}\cos\theta(t)-\frac{GM}{r(t)^{2}}\sin\gamma(t),$$ $$v(t)\gamma'(t)=\frac{F(t)}{m(t)}\sin\theta(t)-\frac{GM}{r(t)^{2}}\cos\gamma(t)+\frac{v(t)^{2}}{r(t)}\cos\gamma(t),$$ $$m'(t)=-\frac{F(t)}{g I_{sp}},$$

given by the following Mathematica code, which models a 2-stage Falcon 9 launch with a 10000kg payload:

G = 6.672*10^-11;(*Gravitational constant*)
M = 5.97219*10^24 ;(*Mass of Earth*)
R = 6.378 *10^6;(*Radius of Earth*)
g = 9.81; (*Gravitational acceleration*)

mprop1 = 409500;(*First stage propellant mass*)
mstruct1 = 22200; (*First stage structural mass*)
mfirststage = mprop1 + mstruct1;(*Total initial mass of first stage*)

T1 = 9*845000 ;(*Thrust from 9 first stage engines*)
Isp1 = 282; (*First stage specific impulse*)

mprop2 = 107500;(*Second stage propellant mass*)
mstruct2 = 4000; (*Second stage structural mass*)
msecondstage = mprop2 + mstruct2;(*Total initial mass of second stage*)

T2 = 934000;(*Thrust from second stage engine*)
Isp2 = 345; (*Second stage specific impulse*)

mpayload = 10000; (*Payload mass*)
mtotal = mfirststage + msecondstage + mpayload; (*Total mass of first stage, second stage and payload*)
mMECO = 0.15 mprop1 + mstruct1 + mprop2 + mstruct2 + mpayload; (*Mass of rocket at main engine cut off with 15% of first stage propellant left over*)
mSECO = mstruct2 + mpayload; (*Mass of rocket at second stage engine cutoff*)
tmax = 517; (*Total burn time*)

F = If[m[t] > mMECO, T1, If[mSECO < m[t] <= mMECO, T2, 0]]; (*Varying thrust for first and second stages*)
Isp = If[m[t] > mMECO, Isp1, If[mSECO < m[t] <= mMECO, Isp2, Isp2]]; (*Varying specific impulse for first and second stages*)

vorbital = Sqrt[(G M)/(R + 250000)] (*Orbital velocity for a desired orbit height of 250km*)

θ0 = 0 Degree;
θf = 0.1085 Degree;
θ[t_] := ArcTan[Tan[θ0] - (Tan[θ0] - Tan[θf]) t/tmax] (*Linear tangent steering law*)

(*Thrust control angle vs. time*)
Plot[θ[t], {t, 0, tmax}, AxesLabel -> {t, θ}, PlotRange -> Full]

(*Solving system of differential equations*)
Solution = NDSolve[{
   r'[t] == v[t] Sin[γ[t]],
   ϕ'[t] == v[t]/r[t] Cos[γ[t]],
   v'[t] == F /m[t] Cos[θ[t]] - (G M)/r[t]^2 Sin[γ[t]],
   v[t] γ'[t] == F /m[t] Sin[θ[t]] - (G M)/r[t]^2 Cos[γ[t]] + (v[t])^2/r[t] Cos[γ[t]],
   m'[t] == -(F/(g Isp )),
   x[t] == r[t] Cos[ϕ[t]],
   y[t] == r[t] Sin[ϕ[t]],
   v[0] == 1, γ[0] == 90 Degree, 
   r[0] == R, ϕ[0] == 90 Degree , m[0] == mtotal, 
   WhenEvent[m[t] == mMECO + 0.1, m[t] -> (msecondstage + mpayload)]}, {v[t], v'[t], γ[t], 
   r[t], r'[t], ϕ[t], x[t], y[t], x'[t], y'[t], m'[t], m[t]}, {t,
    0, tmax}]

(*Plotting rocket's trajectory*)
Show[ParametricPlot[{{{x[t], y[t]} /. Solution}}, {t, 0, tmax}, AxesLabel -> {x, y}, ImageSize -> Large, PlotRange -> Full], Graphics[{Green, Disk[{0, 0}, R]}]]

(*Mass vs. time*)
Plot[m[t] /. Solution, {t, 0, tmax}, ImageSize -> Large, PlotRange -> Full, AxesLabel -> {t, m}]
vtable = Table[m[t] /. Solution, {t, 0, tmax}]

(*Flight-path angle vs. time*)
Plot[180 - γ[t]*180/π /. Solution, {t, 0, tmax}, ImageSize -> Large, PlotRange -> Full, AxesLabel -> {t, γ}]
γtable = Table[180 - γ[t]*180/π /. Solution, {t, 0, tmax}]

(*Altitude vs. time*)
Plot[r[t] - R /. Solution, {t, 0, tmax}, ImageSize -> Large, PlotRange -> Full, AxesLabel -> {t, r}]
rtable = Table[r[t] - R /. Solution, {t, 0, tmax}]

(*Velocity vs. time*)
Plot[v[t] /. Solution, {t, 0, tmax}, ImageSize -> Large, PlotRange -> Full, AxesLabel -> {t, v}]
vtable = Table[v[t] /. Solution, {t, 0, tmax}]

With the above values for $\theta_{0}=0^{\circ}$ and $\theta_{f}=0.1085^{\circ}$ I was able to get a flight path angle of roughly $0^{\circ}$ at engine burnout, however I soon realised that I had no clue how to simultaneously satisfy the 3 constraints (orbital velocity, altitude and flight-path angle) that are required for a stable circular orbit.

For example, if I would like to have my payload reach a stable circular orbit of $250km$, I would also require the rocket's flight-path angle to be $\gamma=0^{\circ}$ (i.e. it's velocity vector is parallel to the local horizon) and would need its orbital velocity to be $v_{orbital}=\frac{\sqrt{GM}}{R+250000}=7753.6 m/s$.

As such, my questions are:

1) Am I even using the linear tangent steering law correctly? (I couldn't manage to find any papers showing its implementation)

2) Can a simple linear tangent steering law be used to simultaneously satisfy the altitude, velocity and flight-path angle requirements needed for a stable circular orbit?

3) As an aside, relative to optimal control approaches, how "optimal" is a linear tangent steering law in terms of fuel use?

Any help would be great!

Below are images of the output for the rocket's (1) trajectory, (2) altitude vs. time, (3) flight-path angle vs. time, (4) mass vs. time, (5) velocity vs. time and (6) thrust control angle theta vs. time

Rocket trajectory

Altitude vs. time

Flight-path angle vs. time

Mass vs. time

Velocity vs. time

Theta vs. time

  • 1
    $\begingroup$ Linear tangent steering is near-optimal to within a few lbs of payload weight in most cases, especially in a model like yours which does not have aerodynamic forces or an uneven gravity field. $\endgroup$ – Nate Berkopec Oct 25 '17 at 12:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.