Philae lander simulation off by factor of ~3

Question

I'm trying to simulate the Philae landing by writing a program to compute the position of the lander vs time. According to various mission websites, the orbiter will match its orbit to the rotation of the comet, then move toward the comet to impart an initial velocity on the lander as it releases it. The lander will "free fall" toward the comet, landing approximately 7 hours later at a velocity no greater than 1.0 m/s.

My simulation will assume a straight-line trajectory from the orbiter to the comet.

I gathered some figures from mission websites:

Mass of comet (M_C): 1.0 x 10¹³ kg (Wikipedia)
Mass of lander (M_L): 100 kg (ESA)
Mean diameter of comet: 4 km (Wikipedia)
Initial velocity of lander (v₀): 0.187 m/s (ESA Blog)
Height of release: 22.5 km (ESA Blog)

I can't use a constant gravitational acceleration because the gravitational pull will be increasing as the lander approaches the comet. So I (with the help of a physics prof.) derived the following formulas. Since the forces will be changing as the simulation progress, we used an iterative set of formulas in which the next set of values is computed from the previous set. In the following, $n$ is the current iteration and $n-1$ is the previous iteration.

$$ a_n = \frac{GM_c}{r^2} $$

$$ v_n = v_{n-1} + a_{n-1}\Delta t $$

$$ r_n = r_{n-1} - v_{n-1} \Delta t - \tfrac{1}{2} a_{n-1} \Delta t^2 $$

(The signs are chosen so that $a_n$, $v_n$ and $r_n$ are non-negative quantities.) We can set up a table of computed values, using $\Delta t = 60$ seconds:

| n | a         | v        | r        |
----------------------------------------
| 0 | 1.1112E-6 | 0.187    | 24500.0  |
| 1 | 1.1122E-6 | 0.18706  | 24488.78 |
| 2 | 1.1132E-6 | 0.18713  | 24477.55 |
etc...

As you can see, the velocity and acceleration are both increasing and the radius (height) is decreasing. That's good.

I stop the simulation when $r$ drops below 2000 meters (the mean radius of the comet).

The problem is the simulation predicts the descent will take ~21 hours. The real descent will take ~7 hours. I'm off by a factor of 3.

I tried changing $\Delta t$ to 1 second, but the numbers don't change much.

Will these formulas work for the simulation? And am I using the correct initial values?

Answer 1

Like mentioned in the comments, you got the initial velocity wrong. You source says:

Further to this, the selected separation strategy foresees a fixed separation velocity of approximately 0.187 m/s. This imposes limitations on the capability of directing Philae towards the comet, i.e. it restricts the domain of possible positions and velocities for Rosetta at the time of separation.

So this 0.187 m/s is not the speed at which Philae is traveling towards the comet, but the relative velocity to Rosetta after separation.

John Rennie did found a source which does mention the actual relative velocity between Philae and 67P at separation. I tried to find another one myself, preferably one from the site of ESA, however I was not able to find one, so I will assume $2\ 1/2\ km/h$ is correct, which is $0.69\ m/s$.

Although your method should yields a good approximation, it is possible to calculate is much easier using the Kepler equation, however the accuracy of this results is up to discussion, since the comet's gravity is far from having spherical symmetry. For a Kepler orbit we do not have enough information, namely the direction of the velocities (radial and lateral components). You assumed complete vertical motion, however I think it would be more likely that at touchdown the lateral velocity would match that of the surface. For this you need to know the position of the landing site and the speed and direction of the rotation of the comet. The landing site is located at the "head" of the smaller of the two lobes of the comet, also called site J and later doped Agilkia:

I was not able to find the actual altitude relative to the center of gravity of the comet, so this will have to be estimated. For this I will be using this video showing the center of the comet and the following image to get the scale:

My best estimation from these two sources gives an altitude of about 2.8 km. This landing site is fairly close to the equator, so the surface velocity would be just the altitude times the angular velocity. The angular velocity can be calculated from the rotation period of 12.4043 ± 0.0007 hours and therefore the surface velocity would be equal to about $0.39\ m/s$. From this it can also be derived that the specific angular momentum, $h$, would be equal to $1.1 \times 10^3\ m^2/s$. To be clear the mentioned height of release is measured from the center of the comet. Combining these properties yields an orbital eccentricity of 1.46 and a semi-major axis of -1.6 km. The impact velocity calculated from these orbital elements would be 0.86 m/s (only the radial orbital velocity component, since the lateral component was chosen to be zero relative to the surface). The time between separation and impact would take about 27000 seconds or about 7.4 hours. If I would use a purely vertical descent it would take 7.3 hours. These results are a lot closer to the actual 7 hours (not sure how exact this number is). When considering that the comet is far from spherical and that the release velocity might as well be close to 2.4 or 2.6 km/h, which already changes in travel time to 7.7 or 7.2 hours respectively, I think that this is a reasonable good approximation.