Modeling a gravity turn - could you help me out with differential equations? ;)

[S1gmoid]

Hey,

I noticed that trying to do proper gravity turns requires some heavy trial and error, so I decided to create a 2-dimensional simulation of a gravity turn, assuming no active steering, except that the thrust is always applied in prograde direction.

For now I'm assuming a vacuum, because it gets complex enough without accounting for changing thrust and drag. This is the model I wrote up:

d[SUP]2[/SUP]s[SUB]x[/SUB](t) / dt[SUP]2[/SUP] = |a[SUB]thrust[/SUB](t)|*(v[SUB]x[/SUB](t)/sqrt(v[SUB]x[/SUB](t)[SUP]2[/SUP]+v[SUB]y[/SUB](t)[SUP]2[/SUP])) = (F[SUB]thrust[/SUB]/(m[SUB]dry[/SUB]+m[SUB]fuel[/SUB]-t*(F[SUB]thrust[/SUB]/I[SUB]sp[/SUB]))) * ((ds[SUB]x[/SUB](t)/dt)/sqrt((ds[SUB]x[/SUB](t)/dt)[SUP]2[/SUP]+(ds[SUB]y[/SUB](t)/dt)[SUP]2[/SUP]))

d[SUP]2[/SUP]s[SUB]y[/SUB](t) / dt[SUP]2[/SUP] = |a[SUB]thrust[/SUB](t)|*(v[SUB]y[/SUB](t)/sqrt(v[SUB]x[/SUB](t)[SUP]2[/SUP]+v[SUB]y[/SUB](t)[SUP]2[/SUP])) - g(t) = (F[SUB]thrust[/SUB]/(m[SUB]dry[/SUB]+m[SUB]fuel[/SUB]-t*(F[SUB]thrust[/SUB]/I[SUB]sp[/SUB]))) * ((ds[SUB]y[/SUB](t)/dt)/sqrt((ds[SUB]x[/SUB](t)/dt)[SUP]2[/SUP]+(ds[SUB]y[/SUB](t)/dt)[SUP]2[/SUP])) - g[SUB]o[/SUB]*(R/(R+s[SUB]y[/SUB](t)))

Now this is a bit more involved than we ever got during calculus class, and if I remember correctly, there probably doesn't exist a closed form expression for s[SUB]x[/SUB](t) and s[SUB]y[/SUB](t). Can someone point me in the right direction for programming a reasonably accurate numerical evaluation of this? Edited November 21, 2015 by S1gmoid

[SomeGuy12]

I'm rather unclear why it needs to be this complicated. You have a velocity and a sum of acceleration vectors. You can trivially approximate this using various numerical integration methods - the simplest one being euler's method. They get a little more complex but you can use a library and you build your programming code in layers - define the function that computes the derivative as a separation function that gets called by the function doing the numerical integration.

At every moment in time, you make the vector component for engine thrust caused acceleration (the other component is for gravity from the planet) the same direction as your vector for velocity. (prograde like you say)

Just break the problem into tiny, logical pieces instead of that horrific mess you have above. Adding changing thrust becomes trivial - call a function to compute the actual thrust magnitude, and in a separation function, subtract fuel burned, recompute vehicle mass, and calculate the new effective acceleration. Drag is similarly easy, assuming you already have drag coefficients.

You don't even need calculus...(well, the numerical integration is technically calculus but Euler's method is braindead simple. The more complex 2nd order + methods are more complex, but in programming terms are still simple)

[PB666]

[quote name='S1gmoid']Hey,

I noticed that trying to do proper gravity turns requires some heavy trial and error, so I decided to create a 2-dimensional simulation of a gravity turn, assuming no active steering, except that the thrust is always applied in prograde direction.

For now I'm assuming a vacuum, because it gets complex enough without accounting for changing thrust and drag. This is the model I wrote up:

d[SUP]2[/SUP]s[SUB]x[/SUB](t) / dt[SUP]2[/SUP] = |a[SUB]thrust[/SUB](t)|*(v[SUB]x[/SUB](t)/sqrt(v[SUB]x[/SUB](t)[SUP]2[/SUP]+v[SUB]y[/SUB](t)[SUP]2[/SUP])) = (F[SUB]thrust[/SUB]/(m[SUB]dry[/SUB]+m[SUB]fuel[/SUB]-t*(F[SUB]thrust[/SUB]/I[SUB]sp[/SUB]))) * ((ds[SUB]x[/SUB](t)/dt)/sqrt((ds[SUB]x[/SUB](t)/dt)[SUP]2[/SUP]+(ds[SUB]y[/SUB](t)/dt)[SUP]2[/SUP]))

d[SUP]2[/SUP]s[SUB]y[/SUB](t) / dt[SUP]2[/SUP] = |a[SUB]thrust[/SUB](t)|*(v[SUB]y[/SUB](t)/sqrt(v[SUB]x[/SUB](t)[SUP]2[/SUP]+v[SUB]y[/SUB](t)[SUP]2[/SUP])) - g(t) = (F[SUB]thrust[/SUB]/(m[SUB]dry[/SUB]+m[SUB]fuel[/SUB]-t*(F[SUB]thrust[/SUB]/I[SUB]sp[/SUB]))) * ((ds[SUB]y[/SUB](t)/dt)/sqrt((ds[SUB]x[/SUB](t)/dt)[SUP]2[/SUP]+(ds[SUB]y[/SUB](t)/dt)[SUP]2[/SUP])) - g[SUB]o[/SUB]*(R/(R+s[SUB]y[/SUB](t)))

Now this is a bit more involved than we ever got during calculus class, and if I remember correctly, there probably doesn't exist a closed form expression for s[SUB]x[/SUB](t) and s[SUB]y[/SUB](t). Can someone point me in the right direction for programming a reasonably accurate numerical evaluation of this?[/QUOTE]


If you are talking about post lauch turn to horizontal you have to include drag.
When i launch from the moon i turn to 45 as quickly as i have verticle ascent, and seconds later to horizontal, which is ussually going to zero thrust as the apo approaches the desired orbital radius.

A key component of vector change has to be launch pressure and atmospheric curve. So on kerbal you want the equilibrium velocity to be 100 m/s at msl. since the atmosphere curve at this altitude is 5000 feet it means that atm is 1/2.6 at 5000 than it is at msl. Given this two bits of information can be added.
1. assumming that you have leading edges at 45' or less relative to the axis of travel, you can disregard mach pressure effects to 280. Drag comperable to msl at 5000 is ~160
2. If this motion was sideways and according to a ground motion of 170 you have 330 m/s for omega squared r which is trivial negation of gravity, you have more loss of gravity due to altitude. Thus it is more productive to invest in altitude.
3. However, pitchin 1-5 degrees down offers some horizontal with essentially no cost of dr/dt

To 10k one can accelerate to 250, but the problem now becomes a drag versus thrust issue, if you have alot of thrust or you coefficient of drag is low you can significantly turn, if you have allot of dv and you thrust is lower or drag is high, keep the vertical ascent past 15k. Essentially hoovering waste dv and proximach drag eats dv at alts below 15k

At 15k the push past mach can begin but significantly turns become harder as momentum increases and a properly tuned ascent does not have the ship glowing red at 36k alt and shuting engines for circularization shortly thereafter. Its a good idea to have a stage event here or throw off some SRBs. A couple reasons, the space engines work well from here, and you want to control a rapid delta degree/delta time turn to 45 ~24k alt, this keeps vel below 1000. Here is were you want full thrust going because we want to stretch apo to 70k, but because r has not drppred g much yet, and because omega is still negating a fraction of gravity angle is up to keep a good climb, by 36 you want 30', 42 alt 20', and above zero until you establishe a desired apo. If you happen to be below 1000 m/s at 36k alt you can apply all that dragless theory, otherwise you are basically keeping things avoiding drag past 42k alt and you dont care at all past 55k.

[mardlamock]

Just do a 4 DoF simulation of the rocket launch with orbital mechanics and drag taken into account, subdivide gravity turn into say 100 steps and use a solver to get an estimate. Its not all that hard. When doing the actual numerical simulation you will have to split the diff equation as many times as needed in order to be able to even do it. Pop an excel sheet to do it first and it shouldnt take you longer than a couple of hours.

[S1gmoid]

Hey, thanks for the insight. I have since found the primer on numerical integrators by the guy behind Principia (the experimental n-body mod), and got a few answers from there.

Yea, by the way overlooking drag is quite a huge limitation from a practical point of view. What first sparked my interest was that after building a pretty large lifter, with relatively low Isp, high thrust and a load of fuel, as it gobbled up its fuel mass, even turned somewhat sideways, the vertical component of thrust overwhelmed gravity, and the gravity turn reversed, leading to wasting all the delta-V in a near-vertical launch (and due to the huge bulk of the thing, it was nearly impossible to steer in the upper atmosphere). :)