Re-entry simulator - Part 1

Ages ago, I've had to make a launch simulator for a multi-staged rocket. The goal was to find the efficiency of multi-stageg rocket compared to single staged ones. That simulator was pretty simplistic; 1D (altitude only), no easy way to change initial conditions, lift was neglected, and, as it was only 2D, there was no way to change orientation or go into orbit.

Now, for my space sciences class (Phy-2100), I have a homework that requires me to calculate the re-entry window of an Apollo capsule so that astronauts don't feel more than 6,5G. Now, I had 2 ways to do it: either calculate a few values by hand and easilly find the maxima and minimum, or do it the Kerbal way; Trying every possibilities and see which sticks. So I took my old Simulink program, and used that as an excuse to fix the problems with it.

First was making it simpler to change initial conditions. That was pretty straightforward; I simply changed some values to be inputs from the workspace. The initial conditions I added are:

• Initial speed (Absolute value and angle)
• Initial weight of propelant
• Initial altitude

Next was turning my 1D sim in 2D...now that was an headache. I'm not great at setting up motion equations, so the following took a few textbooks; along with wikipedia, to set up and make sure I understood. Movements, especially relative motion, and in polar coordinates to boot, are a pain.

First was setting up my hypotheses:

• Coordinates are geocentric.
• The earth surface is a round shell at mean earth radius inside of which relative speed is forced to 0. That shell rotates at 1 turn per 23 hours and 56 minutes.

The equations of movement for polar coordinates are the following (From Wikipedia, because it's easier to copy paste than trying to do rdotdot in blogger):

 

And as such, our forces equilibrium is the following (Again, from wikipedia):

Now, in order to make my things simpler, I need to isolate the acceleration terms.

Now comes the question: what are my forces? and what information do I need to calculate them?

The forces I've chosen to model are the following:

• Drag (Dependent on speed, atmospheric conditions, and geometry)
• Lift (Dependent on the same things as drag)
• Thrust (Depends on ISP, Chamber pressure, angle and atmospheric conditions)
• Gravity (Depends on altitude)

Just for the sake of completeness, here are the things I neglected/assumed:

• The pilot is not Jebediah Kerman
• For this particular case, Drag and lift are hardcoded, and are only dependent on altitude. They are based on this paper:
• Trans-sonic effects (Additionnal drag between mach 0.8 and 1.2) are neglected
• Capsule dynamics (Roll, pitch, yaw and such) are neglected. Capsule is assumed stable at a 25 degrees flight angle.
• Flow is supposed incompressible below mach 1. Past mach 1, the drag coefficient depends solely on altitude.
• No control whatsoever (although this could easilly be added later, along with the capsule/rocket yaw pitch and roll)
• The atmosphere is the standard atmosphere, but it also has a speed, which depends on the altitude. It is 0 at 100km, and the same speed as the earth's shell at sea level. The distribution is assumed linear for now. (Hint: it's not, but that's easier to estimate)

Allright, now that everything is set-up, time to code my life away!

For now I got the orbit part mostly functionning, except for a weird "bug" where my orbit precesses. It shouldn't. I'm guessing it might have to do with rounding errors...time to debug!

EDIT on may 19th: Since it's been less than a day from my previous post, I thought I'd just edit it. I found my bug. Gravity was acting weird, because I had re-used my gravity calculator from my first launch simulator, I had added an "arbitrary" Earth Radius constant. As such, my gravity force was getting lower and lower and lower, due to it adding an earth radius every time! Which explains the precession.

Here's how it looks so far

It works, it's a circle tracer!