Sunday, 19 July 2015

Missing something basic about simple orbital mechanics


I seem to be missing something basic.



I've been trying to get a simple orbital simulation working, and my two objects are Earth around the Sun.


My problem is this. I placed the Earth at 93M miles away from the Sun, or 155M km.


As I understand it, the orbital velocity of something at 155M km from the Sun is:


$$v = \sqrt{\frac{GM}{r}}$$


Plugging in the numbers for the Sun, I get a velocity of:


$$29261 \frac{\mathrm{m}}{\mathrm{s}}$$


However, if I want to get the acceleration that the Sun has upon the earth, I use:


$$g = \frac{GM}{r^2}$$


For the Sun, and 155M km, I get an acceleration of:


$$0.0055\frac{\mathrm{m}}{\mathrm{s}^2}$$



Now, I start with a simple body at the proper radius out along the X axis and give it a simply vector of 29261 m/s along the Y axis, then I start applying the 0.0055 m/s^2 acceleration to it. And the acceleration of the Sun is simply not enough to hold the Earth.


If the Earth starts with a vector of (0, 29261 m/s), and after and I add the acceleration vector of (-0.0055 m/s, 0) to it, you can see that after a single second, it doesn't move a whole lot. If I chunk things to days, 86400 seconds, then the acceleration vector is only, roughly, -477 m/day, but the velocity vector is:


$$2,325,974,400 \frac{\mathrm{m}}{\mathrm{day}} = 29,261 \frac{\mathrm{m}}{\mathrm{s}} \times 86,400 \frac{\mathrm{s}}{\mathrm{day}}$$


As you can imagine, the -477 isn't going to move that much towards the Sun.


I understand that better simulations use better techniques than simply adding basic vectors together, but that's not what this is. I seem to be missing something fundamental. I had assumed that given the correct velocity, that the pull of the Sun should keep the Earth in orbit, but the "pull" that I'm using doesn't seem to be having the desired effect.


So, I'm curious what basic "D'oh" thing I'm missing here.


Edit for Luboš Motl answer.


Perhaps there's something more fundamental I'm missing here. I understand you point, but .0055 m/s * 86,400 is -477. I was doing that math fine.


Simply, I have an object with a velocity vector. Then I apply an acceleration at a right angle. I do that for N seconds to come up with a new, right angle velocity vector. I then add that to the original vector to come up with the objects new vector. I then take that vector, apply to the current position of the object, and arrive at a new position.


Clearly there is a granularity issue which makes some amount of seconds a better choice for a model than others, but this is high school level simple mechanics, so there's going to be some stepping. I chose one day so that my little dot of a planet on my screen would move. If I update every 1/10th of a second "real time", and each update is a day, then I should get a rough orbit that's really a 365ish polygon in a little over 30s real time.



If I choose a step size of 1 second, then my acceleration (0.0055 m/s^2) * 1 s = a right angle velocity vector that's -0.0055 in magnitude. That vector is added to the original vector of 29261 (at right angles), giving me a new vector of (-0.0055, 29261). That's after one second. That's not much of a bump. It's barely a blip. If I apply one days full of acceleration, "all at once", I am obligated to not only multiply the acceleration by 86,400, but also the original vector (since it's 29261 m/s, and we have 86,400 s), thus giving me, proportionally, the same vector, just longer. And it's still just a bump.


So, I'm mis-applying something somewhere here, as I think the numbers are fine. I'm simply "doing it wrong". Trying to figure out what that wrong part is.


Edit 2, responding to Platypus Lover


Thank you very much for the simple code you posted. It showed me my error.


My confusion was conflating updating of the vector with the calculation of the velocity vector. I felt that I had to multiply both the original vector AND the acceleration amount by the time step, which would give me the silly results.


It was just confused in my head.



Answer



I suspect there is something wrong in the way you are adding the acceleration vector to the velocity vector after the first timestep. The simplest way to check you are doing things correctly is to write down your scheme in cartesian coordinates.


To point you in the right direction, I wrote a sample orbit integrator for you here: Simplest orbit integrator, which should be at a level appropriate for a high school student. This is probably simplest integrator you can possibly write. With your code, you should get a nice circle for an x-y plot, and a nice sinusoid for the position and velocity components:


x-y plot t-x plot



Notice that it uses Euler's method:


$\vec{v}(t+\Delta t) = \vec{v}(t) + \Delta t \ \vec{a}(t)$


$\vec{x}(t+\Delta t) = \vec{x}(t) + \Delta t \ \vec{v}(t)$


which is probably what you have been doing so far without realizing. This is the most inaccurate method to use when integrating, and once you introduce elliptical orbits, this method will give you wrong results after a few orbits. There are many simple recommendations I can give you on how to further improve your code once it is fixed (normalizing units, a better integration scheme, etc.).


No comments:

Post a Comment

Understanding Stagnation point in pitot fluid

What is stagnation point in fluid mechanics. At the open end of the pitot tube the velocity of the fluid becomes zero.But that should result...