Aerobraking and reentry

Let's stay with a simple 1 dimentional example for now. And let's assume gravity is the only force, no rockets or drag.

Sigurd,

Does F represent force? If so, I have to agree with nihiladrem that it is not right to say F = V + g.

Humm, It's not force.. "how" the force is achieve in that code is not right.. and how more I think about it how more terrible it is it's only good with that programs input.. atleast visual, and it starts with a velocity of as example 1000, but to make this work for the xls excel file, we need to start with 0 as velocity, add a force (engine), and gravity that changes depending on distance. So the stuff I posted before won't work for a rocket for 100% sure (not saying that the new one more below is right )

I would say you need something like this;
But that's just logical thinking (possible flawed)... not based on observations of how it really is.. so don't ask me if it's right or wrong..

a = (F / m) + g
k1 = a * dt
a = (F / m) + g + (k1/2)
k2 = a * dt
a = (F / m) + g + (k2/2)
k3 = a * dt
a = (F / m) + g + k3
k4 = a * dt
Vnew = V + (k1 + 2*k2 + 2*k3 + k4)/6
Snew = S + Vnew * dt

F as force, of the rocket engine
m as mass, g as gravity on that altitude, don't know if I had to put +g on that place..

Just a guess.. I can only try to search, cause I don't have the answers..

Bearing in mind the chances I have made no serious errors here are somewhat remote.

Say you are integrating using RK4 the path of a point particle in 3 dimensional space moving with an acceleration dependant on position and velocity. If the forcing does not explicitly depend on time, the state at various points in the RK4 iteration can be thought of as being six dimensional quantities labelled y1, y2, y3 and y4. Additionally, six dimensional vectors k1, k2, k3 and k4 would exist each containing a position like quantity and a velocity-like quantity.

Just as y1 would contain a position y1p and velocity y1v, so k1 would contain k1p, a position like part and k1v, a velocity like part. Positions, velocities and accelerations would in the three dimensional case contain three elements each.
e.g.
k1p=(k1px,k1py,k1pz)
k1v=(k1vx,k1vy,k1vz)
y1p=(y1px,y1py,y1pz)
y1v=(y1vx,y1vy,y1vz)

//y1=initial state
//y1=(y1p,y1v)
//k1=(k1p,k1v)
k1p=h*y1v;k1v=h*acceleration(y1p,y1v);

//y2=y1+k1/2
//y2=(y2p,y2v)
y2p=y1p+k1p/2;y2v=y1v+k1v/2;
//k2=(k2p,k2v)
k2p=h*y2v;k2v=h*acceleration(y2p,y2v);

//y3=y1+k2/2
//y3=(y3p,y3v)
y3p=y1p+k2p/2;
y3v=y1v+k2v/2;

//k3=(k3p,k3v)
k3p=h*y3v;k2v=h*acceleration(y3p,y3v);

//y4=y1+k3
//y4=(y4p,y4v)
y4p=y1p+k3p;
y4v=y1v+k3v;

//k4=(k4p,k4v)
k4p=h*y4v;k4v=h*acceleration(y4p,y4v);

//yn4=final_state=yn+k1/6+k2/3+k3/3+k4/6
//y4=(y4p,y4v)
y4p=y1p+k1p/6+k2p/3+k3p/3+k4p/6;
y4v=y1v+k1v/6+k2v/3+k3v/3+k4v/6;

I'll expand the last bit to give clarity(?!)...

y4px=y1px+k1px/6+k2px/3+k3px/3+k4px/6;
y4py=y1py+k1py/6+k2py/3+k3py/3+k4py/6;
y4pz=y1pz+k1pz/6+k2pz/3+k3pz/3+k4pz/6;

y4vx=y1vx+k1vx/6+k2vx/3+k3vx/3+k4vx/6;
y4vy=y1vy+k1vy/6+k2vy/3+k3vy/3+k4vy/6;
y4vz=y1vz+k1vz/6+k2vz/3+k3vz/3+k4vz/6;

With just one (vertical) dimension and gravity alone this would give rise to

k1p=h*y1v;k1v=h*g;
y2p=y1p+k1p/2;
y2v=y1v+k1v/2;
k2p=h*y2v;k2v=h*g;
y3p=y1p+k2p/2;
y3v=y1v+k2v/2;
k3p=h*y3v;k2v=h*g;
y4p=y1p+k3p;
y4v=y1v+k3v;
k4p=h*y4v;k4v=h*g;
y4p=y1p+k1p/6+k2p/3+k3p/3+k4p/6;
y4v=y1v+k1v/6+k2v/3+k3v/3+k4v/6;

Which *should* simplify down to
y4v=y1v+h*g;
y4p=y1p+h*y1v+h*h*g/2;

and with one second time increments just

y4v=y1v+g;
y4p=y1p+h*y1v+g/2;

Nice, I'm sure campbelp2002 should use stuff like this.. and not my code.. since the rocket may go up or down with my code... the results may be really wrong

Edit:

Reading it more in detail... it looks a lot more logical, to keep the "accuracy", it's diffrent than I had RK4 in mind, a lot larger..

(EDIT) Oh man! Two posts came in while I was typing this!

Sigurd,
For simplicity lest us assume no rocket force. Just consider a falling rock. Also, assume dt is 1 second. So F=1 and dt=1. and your equations simplify to this:
a = g
k1 = a
a = g + (k1/2)
k2 = a
a = g + (k2/2)
k3 = a
a = g + k3
k4 = a
Vnew = V + (k1 + 2*k2 + 2*k3 + k4)/6
Snew = S + Vnew

Now find all the k's in terms of g only.

k1 = g
k2 = g + (k1/2) = g + g/2 = 1.5g
k3 = g + (k2/2) = g + 1.5g/2 = 1.75g
k4 = g + k3 = g + 1.75g = 2.75g

To simplify more, assume g is constant. That is approximately true for a rock dropped from a few meters high. In that case, the expression V + (k1 + 2*k2 + 2*k3 + k4)/6 should simplify to V=V+g. Unfortunately:

(k1 + 2*k2 + 2*k3 + k4)/6 =
(g + 2*1.5g + 2*1.75G + 2.75g)/6 =
g*(1 + 2*1.5 + 2*1.75 + 2.75)/6 =
g*(1 + 3 + 3.5 + 2.75)/6 =
g*10.25/6 = about 1.7g

So something seems wrong.

I'm afraid it's going to take some time to fix the mistakes in my own effort. This is exactly why you should keep the integrator independent from what it is integrating at least till it is working...

Optimistically it's just that I've managed to try to define elements of y4 twice when in the second instance I should have been defining y5 - every element should be defined only once. Oh, and forget about the brief mention of yn, that's incompatible with my ordering scheme and I meant to abandon it. Also another similar mistake
k3p=h*y3v;k2v=h*acceleration(y3p,y3v);
becomes
k3p=h*y3v;k3v=h*acceleration(y3p,y3v);
again, this one propogates into the one dimensional incarnation.

It seems to work:

http://hhboard12.free.fr/SIG/RK4.xls

@nihiladrem, can you test if I used it right ?

The step can be excluded from repeating.. but I'm laizy
And gravity will have to change, in the future.

Yup, seems right. I've also done the same test in C. For the general version where acceleration isn't fixed a sufficient test is to make it recreate a few well known standard integrals. If it does so, that version works too.
BTW, thanks for writing this Sigurd!

Nothing better than testing it and playing with it, to make sure you know how it works

I've made a changing gravity version based on:
http://www.glenbrook.k12.il.us/gbssci/p ... u6l3e.html

http://hhboard12.free.fr/SIG/RK4Grav.xls

BUT, I'm wondering if it's right, the formula is really not right when you go below sea level. But maybe I just need to wait, until campbelp2002 add it to it's version with all nice numbers

Good work! Changing the starting y1p to 6378 km (the radius of Earth), the starting y1v to 11.2 km/s (escape velocity) and the acceleration to -398600/y1p^2 (for acceleration of gravity in km/s ^2 varying by the inverse square rule) and running the simulation out to 400,000 km gives about two and a half days to the Moon. It seems about right. For positions below sea level acceleration of gravity varies linearly from 9.8 m/s^2 at the surface to zero at the Earth's center, but that is of no interest for space flight.

Although the steps are so small here you won't easily notice it, unfortunately acceleration changing only once every full step is a no-no for RK4. k1v needs its own gravity calculation based on y1p, the gravity for k2v needs to come from y2p and so on... The same would also need to happen in drag calculations, except they are based on yp and yv.

Humm, true.. all 4 steps use the same gravity, and have a diffrent position.. I'll play more with it later tonight.. when we have more math available, I'm thinking to build it into the site (usable in the forum) some math tools, that can create graphs on the fly etc (but that will ask some CPU and Bandwidth usage, so I'll first need to fix (finish the new news, company profiles sections etc) some other site sections to get higher advertisement revenue... else I can't pay for it all...)

² ³ etc works now.. feww.. it was just a stupid ascii I used instead of UTF-8 , so forum should be easier to use now for math formula's

Hello all. I applaud the discussion of numerical integration techniques. It’s always good to have a means of performing your own analysis, and Euler’s method & RK integration are very good techniques. I have two technical points of my own to add, too.

First, the accuracy of numerical integration can be further improved by decreasing the time step used. For example, instead of calculating change over an entire second of time (or a minute, or whatever), you can run the computations using a time step of 0.1 or even less. This is very easy to do with Euler’s method, and even RK2, although higher order Runge-Kutta integration requires deriving special formulas for variable time steps. Reducing the time step is not elegant, but it can give increased accuracy without adding extra code.

Second, if you want to compare the accuracy of Euler’s method and RK integration, you should be aware that the system of equations in Peter’s spreadsheet has an exact solution.

The radial acceleration due to gravity without factoring in drag is:

ar = - G M / r²

The relationship between acceleration, velocity and position is that acceleration is the first order time derivative of velocity, which is the first order time derivative of position, making acceleration the second order time derivative of position.

r = r
vr = r’
ar = v’ = r’’

You can substitute this relation back into the formula for ar to make it a second order differential equation:

r’’ + G M / r² = 0

Now, the formula Peter used for acceleration had an extra term to account for drag:

r’’ + G M / r² + k (r’)² exp ( (RE – r) / hs ) = 0

If you let z = k hs² exp ( (RE – r) / hs ), then the extra term in Peter’s acceleration equation is just z’’. Substitute that back in, and you get:

r’’ + z’’ + G M / r² = 0

If you let u = r + z, then this can be further simplified to:

u’’ + G M / r² = 0

which is an equation of the same form as the equation for gravitational acceleration without drag, and you can solve it for u’ exactly as you would solve the dragless equation for r’.

If you substitute r and z back in, you can then use the relationship F = d/dr (Potential Energy) to solve the resulting equation for r’.

v = 2 G M ( 1 / r – 1 / ro ) / (1 – z / h )

and the time for the descent is then just the integral of dt = dr / v (which I didn’t bother to do… ).

Not every system of interest is lucky enough to have an exact equation. People don’t go through all the work of coding RK4 routines and tweaking time steps just because they’re too lazy to do the math. But it’s interesting to have an exact solution for comparison from time to time.

Thanks for the discussion.

Picking up here, if ro is the radial distance of the moon, so the atmospheric z terms there are effectively zero. Modulo errors, I get for a trajectory just stopping at the moon.

(r' + z')² / 2 - G M / r = - G M / ro

v² = 2 G M ( 1 / r – 1 / ro ) / (1 – 2z / hs +z² / hs²)




Ha! Good luck to anyone who tries to get an integral for that


[EDIT] - It's been over a day since my post, just in case anyone's interested, what I was going to get round to saying is that C M Edwards' equation doesn't seem to match my interpretation of what you'd expect if
But regardless of my interpretation of it, I think the whole premise for that is faulty. It was a brave attempt, but IMO there's no reason to believe there's a nice simple equation including drag in this situation. I suppose in a way this was what C M Edwards was talking about at the end of his post. There are relatively few systems that are fully integrable and even when you can integrate and get a formula, unless you are going to spend a lot of time trying to understand it, it might not be more useful than a numerical result.