Lunar lander: why is my descent module crashing

Question

I've been reading this awesome simulation of the actual Apollo 17 descent: http://www.braeunig.us/apollo/LM-descent.htm

I've been trying to make my own simulation using Euler integration, with a time step of 100ms. The problem is my descent module is crashing during the first two minutes with a vertical speed of -220 m/s.

During the initial stage of the descent (braking), the lander is tilted 90º and almost all the thrust is spent in the horizontal component, with no force to counter gravity. In my simulation the vertical speed quickly builds up while the pitch is 80-90º, and then gravity does it job. Here's a trace:

    t: 0 s;  x: 0;   y: 17234.6112;  vx: 1697.1264 m/s;  vy: -20.4216 m/s;   ax: 1.2995855485404475 m/s2;    ay: -1.6121564579144938 m/s2 
        pitch: 90;   throttle: 0.1297;   fuel: 8874.2Kg 

    t: 10.015 s;     x: 17062.974621379497;  y: 16948.14950433809;   vx: 1710.2669538608811 m/s;     vy: -36.64853435841286 m/s;     ax: 1.3247309786564347 m/s2;    ay: -1.628484797082333 m/s2 
        pitch: 90;   throttle: 0.1297;   fuel: 8855.501190599978Kg 

    t: 20.019 s;     x: 34239.89169813059;   y: 16498.912567960208;  vx: 1723.6514116543285 m/s;     vy: -53.026163503416676 m/s;    ax: 1.3509707026458686 m/s2;    ay: -1.6456786784409965 m/s2 
        pitch: 90;   throttle: 0.1297;   fuel: 8836.82291908347Kg   

    t: 30.026 s;     x: 51527.03405575163;   y: 15884.74450896862;   vx: 1725.5675762461628 m/s;     vy: -69.5834660661016 m/s;  ax: -0.9584923539027193 m/s2;   ay: -1.663174744204344 m/s2 
        pitch: 90;   throttle: 1;    fuel: 8756.011226046036Kg 

    t: 40.041 s;     x: 68759.3217349438;    y: 15103.315677459259;  vx: 1715.7571868391299 m/s;     vy: -86.32821281208818 m/s;     ax: -1.0003892381595199 m/s2;   ay: -1.6806337716945954 m/s2 
        pitch: 90;   throttle: 1;    fuel: 8611.841531365935Kg 

    t: 50.502 s;     x: 86651.75776002374;   y: 14107.047450672597;  vx: 1705.0564069549016 m/s;     vy: -104.00697912916615 m/s;    ax: -1.0452064683102156 m/s2;   ay: -1.6991809596193168 m/s2 
        pitch: 90;   throttle: 1;    fuel: 8461.251498809865Kg 

    t: 134.865 s;    x: 226365.92349700694;  y: -5.915270639694082;  vx: 1602.384872970266 m/s;  vy: -225.9485740506923 m/s;     ax: -1.382111636427434 m/s2;    ay: -1.224439728661491 m/s2 
        pitch: 78;   throttle: 1;    fuel: 7246.814359289593Kg 

    FINISHED at t=2.24775 minutes

However, in the original simulation, vertical speed never exceeds -20 m/s. I don't know how could this be possible, as there's no drag due to atmosphere. And I've coded initial acceleration and velocity values exactly as in this very same page. You can check the original trace here:

http://www.braeunig.us/apollo/LM-descent.pdf

Here we can see the Moon's gravity is aprox. -1.5 m/s which heavily influences ay in my simulation. However in the real descent trace it was kept to lower values (+-0.1) most of the time, even during the first minute when there was no vertical thrust due to the engine.

I must be making very stupid mistakes but right now I can't see where.

BTW this is how I update my vars in each step:

    //Calculate acceleration as a function of mass, velocity, thrust and gravity
    ax = vx*vx/(y + MOON_MEAN_RADIUS_M) + (thrust_x/mass);
    ay = (vx*vy/(y + MOON_MEAN_RADIUS_M)) + g + (thrust_y/mass);

    // This would be the acceleration for a flat moon. Don't work either.
    //  ax =  thrust_x/mass;
    //  ay = g + (thrust_y/mass);

    //Update velocity as a function of acceleration and time.
    vx = vx + ax*dtSecs;
    vy = vy + ay*dtSecs;

    // Update position
    x = x + vx*dtSecs;
    y = y + vy*dtSecs;

Any hint would be appreciated.

Thanks in advance.

Answer 1

You have switched the two acceleration formulas around, and you forgot a negation, so:

ay = vx*vx/(y + MOON_MEAN_RADIUS_M) + (thrust_x/mass);
ax = -(vx*vy/(y + MOON_MEAN_RADIUS_M)) + g + (thrust_y/mass);

Basically this is computation of the centrifugal force.

But the centrifugal force isn't really a force, it is just an artefact of using a rotating coordinate system.

So a better way to explain it may be that it is a confusing way of introducing coordinate system rotation. As the rest of the simulation presume that the coordinate system is aligned according to local gravity, the coordinate system has to rotate with the landers movement around the Moon.

Calculating new coordinates when rotating a coordinate system by an angle of α counter-clockwise:

x' = cos(α)x + sin(α)y
y' = -sin(α)x + cos(α)y

But we would only ever turn the coordinate system by a very small amount, so small that it is acceptable to use limits to simplify the formula:

for α → 0
sin(α) → α
cos(α) → 1

Now we need to know what α is, that depends on the landers horizontal movement and the distance to the centre of the Moon:

α = dx / [distance to Moon]

Where dx is the amount x changes.

In total in code with clockwise rotation in the positive direction we get:

vx = vx - vx * timestep / moondist * vy
vy = vx * timestep / moondist * vx + vy

That can be shortened to:

vx += -vx * timestep / moondist * vy
vy += vx * timestep / moondist * vx

And as the job of acceleration in code is to be multiplied by timestep and added to velocity it is easy to integrate this into acceleration.

About the task

I find two things wrong about this task. First of all giving you a "magic" formula with no explanation doesn't help you to learn the relevant mathematics. Second, who the hell use rotating coordinate systems in simulation? Stick to one reference frame and resign to the fact that up and down are not universal directions. It is so much easier in the long run.