Tidal forces simulator

Interface

Launching the tool will bring you to the rappture interface, where you will be able to set the simulation parameters. An overview of using the rappture interface can be found at the using rappture page. After you have set the parameters to the values desired, click “simulate” to run the simulation. A progress bar will show up at the bottom. Below is a list of parameters and their description.

• Apogee (ratio)
  This is the furthest away the cloud will get away from the planet expressed as a ratio of the Roche limit radius. For example: a value of ‘3’ means that the cloud’s maximum distance from the planet would be three times the Roche limit.
• Perigee (ratio)
  This is the closest the cloud will get to the planet expressed as a ratio of the Roche limit radius. A value of ‘0.5’, for example, means that the cloud’s closest approach will be half of the Roche limit.
• Mass ratio
  This parameter set the ratio of the planet (central) mass to the cloud mass. For example, if you want to simulate the Earth-Moon system, you would set the mass ratio to 81.3 because the Earth is 81.3 times heavier than the Moon.
• Number of particles
  This parameters sets the (approximate) number of particles in the cloud that will be simulated. The actual particle count in the simulation will deviate from the value entered because of the way the particle cloud is generated. Details on the cloud generation method can be read in the cloud generation section. Computation time will increase with the square of particle count.
• Seed
  This parameter sets the seed used for the random number generator. Using the same seed from a previous simulation (and keeping other parameters identical) will generate the exact same result. Changing the seed will change the random particle velocities when the cloud is generated.
• Number of orbits
  This parameter sets number of orbits the cloud will make around the planet. Increasing the number of orbits will increase the computation time proportionally.

Simulation outline

After the parameters are set in the interface, the simulation goes through several steps: initializing the cloud, calculating the timestep, running the n-body simulation, and displaying the results. Click the following links to read more about the theory and steps involved.

• Background
• Cloud generation
• Timestep calculation
• N-body simulation
• Visualization

Background

The Roche limit the distance between two bodies where the gravitational force becomes zero on the surface of the lighter body. The Roche limit can be calculated directly from Newton’s law of gravitation, resulting in

\(r_{roche}=r\left(\frac{2M}{m}\right)^{\frac{1}{3}}\)

where \(m\) and \(r\) are the mass and radius of the lighter body, and \(M\) is the mass of the heavier body.
As you approach the limit, the tidal forces become significant compared to gravitational force on the surface. As a particle cloud moves closer to the Roche limit, we expect the cloud to spread in its orbit more quickly.

Cloud generation

A particle cloud is used to visualize the tidal forces because simulating a solid body, like an asteroid, is very hard to do. A constant-density particle cloud is required for the cloud to be stable. The cloud can be generated by creating a cube in cartesian coordinates and carving a sphere out of it, but such a cloud lacks spherical symmetry. A spherically symmetric cloud is needed. Consider the equation for the mass of sphere:

\( M=\int_0^r\int_0^\theta\int_0^{\phi}\rho \left(r,\theta,\phi\right)r^2 \sin\left(\theta\right)\,dr\, d\theta\, d\phi \)

Assuming constant density:

\( M=-\rho\left(\frac{1}{3}\right)r^3\cos\left(\theta\right)\phi \)

If particles were placed equally spaced along the radial coordinate, the density would decrease from the center outward as \( r^{-3} \). To have constant density, the distance from the center of the cloud to a particle should increase in intervals proportional to \( r^{-3} \). A similar argument holds for the altitude coordinate (theta). If simple constant angle spacings are used, then there will be a larger density closer to the vertical axis. The density would increase according to cos(theta). Thus, the angle should increase with constant intervals of cos(theta). The azimuthal coordinate remains unchanged, as seen in the latter equation. We now know how to construct a sphere of average constant density.

This particle cloud is in virial equilibrium: the gravitational potential energy is twice the kinetic energy, or

\( 2\langle T\rangle=\langle U\rangle \)

where \( \langle x\rangle \) denotes the time average of the quantity \(x\).
The virial theorem gives us an order-of-magnitude approximation for the velocities inside the cloud. Each (cartesian) component of the velocity is chosen randomly from a normal distribution. The mean of the distribution is zero -- on average the same number of particles will be going in the positive coordinate direction as the negative direction. The variance of the distribution is determined by the virial theorem and an empirical constant. The variance is

\( \sigma^2 = C\frac{G M}{R} \)

where \( C = 1/5 \). \(C\) was chosen to be \(\frac{1}{5}\) because that is what gave the best results after testing. Some particles have speeds greater than the speed needed to escape the cloud (the escape velocity). The speed of these particles are reduced to half the escape velocity so that all particles are initially confined to the cloud.
The constant average density, particle velocity distribution, and fast particle culling method lead to a confined and well-behaved cloud of particles.
Now that the cloud is created, we must insert the cloud into an orbit around a central mass. The orbital velocity is given by

\( v_{orbital}=\left(G M_{central} \frac{b^2 A_p}{a}\right)^{1/2} \)

where \(a\) is the semimajor axis, \(b\) is the semiminor axis, and \(A_p\) is the apoapsis. The positions of the particles in the cloud are set to their current position plus the initial orbital position. The velocities are set to their current velocities plus the orbital velocity. This process inserts the particle cloud into an orbit around the central mass (the planet).

Timestep calculation

The period of the orbit can be found by the equation

\( \tau=2\pi\left(\frac{a^3}{G M_{central}}\right)^{\frac{1}{2}} \)

The timestep is then

\( \Delta t = \frac{N_{\frac{steps}{orbit}}}{\tau} \)

The timestep defaults to a value of 1600 steps per orbit at a distance of three times the Roche limit. The Roche limit distance is arbitrary, but it does not affect the simulation output because all orbital and cloud parameters depend on this value. This timestep is fixed for all simulations regardless of the number of orbits or major axis so that accurate comparisons can be made between any set of orbits.

N-body simulation

Newton’s law of gravity

\( F=G\frac{M m}{r^2} \)

is used to calculate the force on each particle on every other particle. This simulation is an “N-body” simulation because there are N bodies interacting gravitationally with one another. The number of force calculations increases as \( N^2 \).
As \(r \rightarrow 0\), the force between the masses approaches infinity; the particles will scatter and escape the system if they are too close. A force softening length, \(r_0\), must be introduced to eliminate this behavior. Newton’s law of gravity is modified so that

\( F=G\frac{M m}{\left(r+r_0\right)^2} \)

When \(r \rightarrow 0\), the force on the particles will be

\( F=G\frac{M m}{r_0^2} \)

At particle distances much greater than the force softening length, \(r \gg r_0\), the softening length term is negligible. The force-softening length is calculated by finding the inter-particle distance

\( \left( \frac{4}{3}\pi \frac{r_{cloud}^3}{N_{particles}}\right)^\frac{1}{3} \)

and dividing it by the softening length ratio of 1000. The softening length depends on the cloud radius so that very small clouds can be simulated realistically.
The acceleration of the particle depends on the position:

\( a_i = \ddot{x_i} = \sum_{j=1}^N \frac{G M}{\left(x_j - x_i\right)^2} \)

To find the positions of the particles as a function of time, the force equation must be twice integrated through time.
Many methods of numerically integrating the equations of motion. The method chosen for this simulation is the velocity verlet algorithm because it is fast and conserves the total energy of the system. The algorithm is shown below:

\( x_i\left( t + \Delta t \right) = x_i\left( t\right) + v_i\left( t\right) \Delta t + \frac{1}{2} a_i\left( t \right) \Delta t^2 \)
\( a_i\left( t + \Delta t\right) = \sum_{j=1}^N \frac{G M}{\left(x_j\left(t + \Delta t\right) - x_i\left(t + \Delta t\right)\right)^2} \)
\( v_i\left(t + \Delta t\right) = v_i\left( t\right) + \frac{1}{2} \left(a_i\left( t\right)+a_i\left(t + \Delta t\right)\right) \Delta t \)

In pseudocode, the algorithm is implemented as

   
x += v*∆t + 1/2*a*∆t^2
a_new = calculate forces(x)
v += ∆t * (a + a_new)/2
a = a_new

Each repetition is one timestep (which is a constant time interval).

Visualization

The orbit is drawn by analytically calculating the orbit. The Roche limit is a circle drawn with a radius of the Roche limit distance. The cloud particles and central mass particle are drawn as points.
In the lower part of the window, histograms of the angular displacement, the radial displacement, and vertical displacement are shown. These histograms make the spreading of the cloud in each dimension easier to see. The data for the histogram are calculated by comparing each particle to a test particle in the same orbit.
SFML is used to display text and create the OpenGL context. The scene is rendered using OpenGL. Screenshots and movies can be created and are saved in a zip file after the simulation visualization is closed.