Multibody dynamics 1

N-bodies

In this session we will further expand the physics simulation to incorporate multiple bodies interacting with one another.

Forces as Relationships

The core update of the physics logic is in the force analysis, where we introduce the general form of Newton’s gravity relating two bodies [>] i.e. F = GmM/d2, where “G” is a gravitational scaling constant, “m” is the mass of one body, “M” is the mass of another body and “d” is the distance between them.

Unlike the form F = mg used for projectile motion, the new equation expresses relationships between two different objects. Hence we may consider the former as a “unary” type, while the later as a “binary” operation between particles; in the same sense negation and multiplication operations work for numbers.

Moreover, we are interested to understand how multiple objects, say “n” number of objects, interact with one another in an all-pairs sense. The implementation below follows an “n-ary” force application, or colloquially everything-to-everything.

The computation work like this: First we create vectors from every particle to every other particle. This is achieved by subtracting particle positions as usual. Having the second operand grafted causes the creation of all pair-wise combinations. In the figure below we can see that supplying four positions produces a table of four by four vectors. Each vector represents the line connecting two bodies.

The equivalent symbolic tabular form is seen below, where the “A” input parameter is represented by row and “B” by columns. Using the Cartesian product behaviour of grafting produces all pair-wise combinations conveniently and for free.

${\bf A - B}$	${\bf p_{0}}$	${\bf p_{1}}$	${\bf p_{2}}$	${\bf p_{3}}$
${\bf p_{0}}$	$p_{0} - p_{0}$	$p_{0} - p_{1}$	$p_{0} - p_{2}$	$p_{0} - p_{3}$
${\bf p_{1}}$	$p_{1} - p_{0}$	$p_{1} - p_{1}$	$p_{1} - p_{2}$	$p_{1} - p_{3}$
${\bf p_{2}}$	$p_{2} - p_{0}$	$p_{2} - p_{1}$	$p_{2} - p_{2}$	$p_{2} - p_{3}$
${\bf p_{3}}$	$p_{3} - p_{0}$	$p_{3} - p_{1}$	$p_{3} - p_{2}$	$p_{3} - p_{3}$

With a bit of compaction i.e. replacing point differences with vectors $u_{ij}$, we can see that the table is symmetric and the diagonal is zero. This makes perfect sense as it represents the idea of one body seen from the perspective of another; even by itself.

${\bf A - B}$	${\bf p_{0}}$	${\bf p_{1}}$	${\bf p_{2}}$	${\bf p_{3}}$
${\bf p_{0}}$	$0$	$-u_{10}$	$-u_{20}$	$-u_{30}$
${\bf p_{1}}$	$u_{10}$	$0$	$-u_{21}$	$-u_{31}$
${\bf p_{2}}$	$u_{20}$	$u_{21}$	$0$	$-u_{32}$
${\bf p_{3}}$	$u_{30}$	$u_{31}$	$u_{32}$	$0$

Next we compute the square of the distances between particles by taking the dot product between all vector pairs. This works because the dot product of a vector by itself is $u · u = |u| * |u| * cos( a )$ and with the angle being zero it reduces down to the desired $|u|^2 = u · u$ . Again we note that the diagonal elements are zero because the distance of a particle from itself is zero.

We compute the inverse square distance using a python component named $f(x) = x * e^{2 pi i \xi x}$1 / A. The expression of this component is B = 0 if A == 0 else 1 / A. It represents “safe” division where the result is zero if the denominator is zero otherwise it returns one over the input value.

This may appear as a bizarre hack but it actually makes sense: we want to compute the force magnitudes between each pair of particles. But doing so for a particle with itself does not make sense. Zeroing the force thus expresses the idea of ignoring.

We continue computing the pair-wise mass products using the Cartesian product approach used for force directions. This produces a table of of scalars.

We conclude the force magnitude calculation by factoring the inverse distances with the mass products and finally the gravitational constant as seen below.

The force vectors are expressed as the product of the normalized direction between particles times the gravity force magnitude. Finally, we need to convert the table of partial pair-wise forces into a list of sum of force vectors per particle. This is trivially performed using the summation component which causes the table’s row or column values to collapse into the desired outcome.

${\bf A - B}$	${\bf p_{0}}$	${\bf p_{1}}$	${\bf p_{2}}$	${\bf p_{3}}$
${\bf p_{0}}$	$F_{00} = 0$	$F_{01}$	$F_{02}$	$F_{03}$
${\bf p_{1}}$	$F_{10}$	$F_{11} = 0$	$F_{12}$	$F_{13}$
${\bf p_{2}}$	$F_{20}$	$F_{21}$	$F_{22} = 0$	$F_{23}$
${\bf p_{3}}$	$F_{30}$	$F_{31}$	$F_{32}$	$F_{33} = 0$
	$\sum F_{i0}$	$\sum F_{i1}$	$\sum F_{i2}$	$\sum F_{i3}$

This concludes the major modification of the graph’s logic. Overall, we used a table to capture the idea of computing all pairs of gravity forces. Interactions are expressed as the intersections of columns and rows of all-combinations table. There is nothing special about the use of matrices other than they concisely allow us to perform the computations.

Practice

Tasks

• The difference between Newton’s gravity law and Coulomb’s electromagnetic force law [>] is small. You can re-interpret the particle’s masses as charges i.e. allow them to take both positive and negative values to produce both attraction as well as repulsion forces!
• As an advanced challenge you may implement what is known as swarming or flocking behaviour [>] using the same techniques. The article [>] explains how to express the relevant forces.

Hints

2-Body Setup

The first demonstration of the new force setup involves two bodies representing the earth revolving around the sun. Despite its simplicity it is quite amazing to visually verify that just a magical combination of mass ratios, distances and velocities can keep a planet in orbit producing circular motion out of nowhere.

The sun is positioned at the world’s origin, has no initial velocity and 10,000 units of mass. Earth is situated 100 units along the X-direction, has 10 units/sec velocity along the Y-direction and 10 units of mass. The numbers are relative after all.

Practice

Tasks

• Try to modify the initial parameters of this setup such as the masses, velocities and positions and observe the catastrophic effects.

• Draw a circle from the origin with radius 100 units and for each time step visualize the deviation of the earth’s position from the ideal orbit.

• This is also a good way to appreciate the stability of the integration algorithm over time ie. whether we can observe drift after say 100 revolutions.

Hints

58KB

5B_Sun-and-Earth.gh