newtonian mechanics - Calculating distribution of force in a N-body system of balls

With three balls placed in a pyramid on a static ground and under the influence of gravity, how will the force of gravity be distributed? The balls have different radius and mass. The initial velocity is zero and all contact points are known. The goal is to calculate the acceleration and the new velocity for a small timestep so that the new positions do not cause intersections in previous contacts.

I'm trying to implement this for a custom 2D physics engine but I'm really struggling. It seemed so simple first but no matter how I try I can't get it right. Maybe it's a lot harder than it looks? This is of course just an example, I can have any setup and number of balls (but just balls for now). I'm not considering friction but if it's not that more complicated I would like to add that too. And rotation would be awesome of course! Anyway, I'm satisfied with some basic explanation on how to think about this and in particular with respect to Euler integration and generalized software algorithms.

An approximate solution is enough as long as it looks close to natural. Also, if contacts are broken as expected I will instead calculate new collisions in another way. What I'm looking for here is just how it works with multiple contacts when applying forces or gravity.

Note: When having multiple collisions, for instance one ball colliding with two balls at the exact same moment (perhaps unlikely) one can resort to calculate one collision at a time and get away with that pretty good. But when like in this case having a stack of balls pushing against each other in multiple contact points I guess it's not possible to calculate each contact separately anymore, or is it? Perhaps it's just a variant of the same problem? Or is there a better and quicker solution here that will produce an accurate result no matter the setup?

Answer

I presume that you are working in 2D. When you have an algorithm for that it should not be difficult to extend it to 3D. So I shall call the shapes circles.

Your question is similar to Reaction forces in pyramid stacking of steel coils which I answered by claiming that the structure is statically indeterminate. The reactions between the rigid bodies cannot be found because there are not enough independent constraints. The elastic properties and deformations of the circles must be considered to obtain more constraints. The Direct Stiffness Method is an efficient way of doing this, and is particularly suitable for computer programming.

In your case the centres of the circles become the nodes of a truss, and the reactions between them are the compression forces in the struts between nodes. You will need to include 2 equal and opposite lateral forces to ensure that the structure is initially in equilibrium. When you have determined the static reactions between the circles you then remove the lateral forces and calculate how the circles respond to the unbalanced forces on them.

Making the contacts frictionless will avoid having to consider the rotational dynamics of the circles.