(B9) - Gravity Between Spheres

I just made a little 2D gravity toy using the javascript Paper library. Most similar n-body simulators on the web handle collision by simply checking if the circles are intersecting each frame. If they are intersecting, then the masses combine into a bigger circle and the velocity is updated so as to conserve momentum. This works well enough for most purposes, but oftentimes two bodies have so large a relative speed that they jump past one another from one frame to the next and the collision test fails. For this reason (among various others), I decided to turn collision off entirely. Indeed I think the simulation is more fun to play with when there is no collision anyway.

It is important to note that even though the simulation is 2-dimensional, the gravity masses are still modeled as spheres. In fact, mass doesn't really make sense in a 2D world. If you did introduce some sort of 2D mass-density, the force of "gravity" actually becomes proportional to $r^{-1}$ rather than the inverse-square law we are used to. Indeed, if you try to calculate the force of gravity on the plane of a uniform disk in accordance to the inverse-square law, you actually get an infinite force on the edge of the disk. But I digress.

When we turn collision off, we are faced with the problem of calculating the force of gravity between intersecting spheres. Suppose we have two uniform spheres $A,B$ with radii $R_a, R_b$ and masses $M_a, M_b$. Without loss of generality we will assume $R_a \geq R_b$. We denote the vector direction between the sphere centers by $\boldsymbol{\hat{r}}$ and the distance by $r$. By symmetry, we know that the force will point along $\boldsymbol{\hat{r}}$, so we can write it as $\boldsymbol{\vec{f}} = f(r) \boldsymbol{\hat{r}}$. Considering Newton's 3rd law, we need only concern ourselves with calculating the strength of gravity $f(r)$ that $A$ exerts on $B$. We break up $f$ into three intervals

$$f\,(r) = \left\{ \begin{array}{rl} f_{\mathtt{in}}\,(r) \quad \quad & r \leq R_a - R_b \\[1em] f_{\,\cap}\,(r) \quad \quad & r \in (R_a - R_b, R_a + R_b) \\[1em] f_{ext}\,(r) \quad \quad & r \geq R_a + R_b \end{array} \right.$$

Let's first review the external case where the spheres are far enough apart that they do not touch. By the shell theorem we can represent each sphere by a point-mass at its center and we immediately find the strength of gravity $f_{\mathtt{ext}}\,(r) = \frac{M_a M_b} {r^2}$ (where we ignore the gravitational constant).

Now let's suppose that $B$ is completely contained within $A$, that is $r \leq R_a - R_b$. Consider a point mass $p$ inside of $A$. By the shell theorem we know there is zero gravity within a spherical shell, so we can ignore the part of $A$ that lies outside of $p$. If we denote the mass density of $A$ by $\rho_a = \frac{3 M_a}{4 \pi R_a^3}$, then the effective mass that acts on $p$ is $\frac{4}{3} \pi \rho_a r^3$ and the strength of gravity becomes $\frac{M_a M_p}{R_a^3} r$ pointing from $p$ towards the center of $A$.

Without loss of generality we can choose our coordinate system so that $A$ lies at the origin and $B$ lies on the $x$-axis at $(r,0,0)$. By symmetry, we notice that the $y$-forces will cancel out and that the $x$-forces will average out to the center of $B$. The strength of gravity is thus $f_{\mathtt{in}}\,(r)=\frac{M_a M_b}{R_a^3} r$. Indeed, since the force is linear in $r$, any (non-uniform) rigid body of mass $m$ contained inside of a uniform sphere of mass $M$ and radius $R$ will experience a net gravitational force $\frac{mM}{R^3}r$ in the direction from the body's center of mass towards the center of the sphere.

Now let's consider the interesting case where $A$ and $B$ are intersecting. We again choose our coordinate system so that $A$ lies at the origin and $B$ lies at $(r,0,0)$.

The intersection $A \cap B$ consists of two spherical caps pushed together. On the left side of $q$ we have a cap of $B$ and on right side we have a cap of $A$. We notice that as we move along the $x$-axis we get cross sections of closed disks. With this in mind we can solve for $q$ by combining the surface equations of the two spheres.

$$\begin{array}{rcl} x^2 + y^2 + z^2 & = & R_a^2 \\ (x-r)^2 + y^2 +z^2 & = & R_b^2 \\ \end{array} \quad \Longrightarrow \quad q = \frac{r^2 - R_a^2 + R_b^2}{2r}$$

We can breakdown the gravity calculation as follows.

To find $f_{\cap_\mathtt{in}}$ we simply need to find the total mass $M_{A \cap B}$ and the center of mass $(C_{A \cap B}, 0, 0)$ of the intersection. Let $(C_1,0,0), \, (C_2,0,0)$ be the center of masses and $V_1, \, V_2$ be the volumes of the left and right caps, respectively. The equations for such can be found on Wolfram Mathworld. We then find

$$\begin{array}{rcl} M_{A \cap B} & = & \rho_b (V_1 + V_2) \\[1em] C_{A \cap B} & = & \displaystyle \frac{C_1 V_1 + C_2 V_2} {V_1 + V_2} \\ \\ f_{\,\cap_\mathtt{in}}\,(r) & = & \displaystyle \frac{M_a M_{A \cap B}}{R_a^3} C_{A \cap B} \\[1em] & = & \displaystyle \frac{M_a M_b}{32 R_a^3 R_b^3} \frac{1}{r^2} \left((r - R_b)^2 - R_a^2\right)^2 \, \left(r^2 + 4R_br - R_a^2 + R_b^2\right) \end{array}$$

Now on to $f_{\, \cap_\mathtt{out}} = f_{\mathtt{ext}} - f_{\mathtt{cap_1}} - f_{\mathtt{cap_2}}$. We already know $f_{\mathtt{ext}} = \frac{M_a M_b} {r^2}$. To find $f_{\mathtt{cap_2}}$ we first find a general formula for the gravity between a ring and a point-mass on the axis thru its center. We then use that result to find the force between a point-mass and a disk, and use that result to calculate the gravity between a point-mass and a spherical cap (where the point-mass lies at the center of the sphere).

Now we simply plug and chug $[m, \rho, R, h] \mapsto [M_a, \rho_b, R_a, R_a-q]$ into $f_{\mathtt{cap}}$ to find $f_{\mathtt{cap_2}}$

$$f_{\mathtt{cap_2}} = \frac{3 M_a M_b}{16 R_a R_b^3} \frac{1}{r^2} \left((r-R_a)^2 - R_b^2\right)^2$$

Finding $f_{\mathtt{cap_1}}$ is a little more difficult, but the concept is the same. We integrate a disk element over $x$.

We can then plug-in the corresponding variables $[m, \rho, R, h, d] \mapsto [M_a, \rho_b, R_b, R_b+q-r, r]$ to get

$$f_{\mathtt{cap_1}}\,(r) = \frac{M_a M_b}{4 R_b^3 r^2} \left( r^3 \, - \, 3 R_a \, r^2  \, + \, 3 (R_a^2 - R_b^2) \, r \, - \, (R_a-2R_b) (R_a + R_b)^2 \right)$$

Finally, putting everything together and simplifying using a symbolic manipulator (I used Sympy and Mathematica) we find

$$\begin{array} { lcl }     f_{\mathtt{in}}\,(r) & = & \frac{M_a M_b}{R_a^3} r \\     f_{\,\cap}\,(r)            & = & \frac{\;\; M_a M_b} {32 \, R_a^3 R_b^3} \frac{1}{r^2} \; \left(         \begin{array}{rcrcl}             r^6 & - & 9 \, (R_a^2 + R_b^2) \; r^4 & + & 16 \, (R_a^3 + R_b^3) \; r^3 \\             & - & 9 \, (R_a^2 - R_b^2)^2 \; r^2 & + & (R_a-R_b)^4 \, (R_a^2 + 4R_aR_b + R_b^2)         \end{array} \right) \\     f_{\mathtt{ext}}\,(r) & = & \frac{M_a M_b} {r^2} \\ \end{array}$$

To check our work it is always a good idea to plug in a few concrete values. Here is the graph for when $M_a = M_b = 4$ and $R_a=2$ and $R_b=1$.

We notice that this curve appears smooth at the seams. Indeed, $f$ is actually $C^2$ continuous, meaning both the first and second derivative are also continuous. In hindsight, we might have been able to guess that $f_{\,\cap}$ takes the form $\psi(r)/r^2$ where $\psi$ is the unique $\mathtt{6}$th order polynomial that matches $f_{\mathtt{in}}$ and $f_{\mathtt{ext}}$ and the first two derivatives at the seams. Interpolating $\psi$ perhaps might have been easier than taking the integrals above ¯\_(ツ)_/¯