I’m surprised that this question has such few views and votes. Stacking is a fundamental pathological use case of physics engines that implement a simple contact solver model. It showcases how the principle behind simulating the dynamics on an individual per-contact basis breaks down when long chains or complex interactions of contact are at play.
In the typical situation of a general purpose engine and a stack of boxes, I’d analogize the situation with a chain of nodes that are connected. In both situations the forces computed between them are integrated one at a time only at their contact boundaries.
E.g. as far as the solver is concerned, the top box isn’t interacting with any other box than the second-to-top box. Therefore, in the course of one timestep, only the penetration resolution impulse (or positional correction, as necessary) from the second-to-top box can affect the top box.
So let’s imagine 5 nodes connected with constraints
T=0: * - * - * - * - *
They must be 1 unit of space between them (because for example the blocks are one unit thick). Represented here by the distance of 3 monospace characters.
A force (such as gravity, but let’s assume for the sake of easy analysis this stack is static in a weightless environment and we experience a single impulse) is experienced on the top of the stack pushing it down.
T=1: * * - * - * - *
T=2: * - * * - * - *
T=3: * * - * * - *
T=4: * - * - * * *
...Etc. Interactions become more diffused.
The point here is that the way the “speed of sound” has to travel through the material shows how it is limited by your timestep and number of iterations. This is where the realism breaks down. Indeed if you double the number of blocks stacked together (and halve their height), you have something that is more or less twice as unstable in the simulation,
In the real world, there is no such thing as a timestep or iterations, and yielding to a squishing force should not depend in any significant amount on the number of discrete structures involved (we do not see that books with a higher count of pages are proportionately more squishy, for example. Nor that books are the squishiest of all things.)
We understand that the correct interaction is that no squishing of the stack happens (provided we are modeling an incompressible block material, of course):
T=1: * - * - * - * - *
In order for it to be possible to have this “correct” computation, the engine has to realize that all of the blocks are connected to each other based on their geometric contact state, and perform a global constraint resolution on that connectivity. In our simplified node representation, node 0 and node 4 are connected through a chain of 3-character spacing constraints, so upon node 0’s receipt of an impulse, the solver must somehow determine that since node 4 will not “give” in that direction, the normal force impulse shall be delivered immediately, canceling out the squish motion.
I think that this is all quite possible to resolve in theory, and the TOI solver for the computation of the time of impact by sweeping in the box2d engine is an example of an engine taking some steps forward (though that one does not do much in the way of directly addressing massive contact interactions such as large stacks: it simply allows for detecting and adding “in-between” iterations at the time of impact, so that bullet-like fast moving objects will not freely warp through objects without interacting with them).
So this has been my approach to explaining the situation around poor stacking performance in typical rigid body simulations and how to start thinking about improving the situation.
The trouble with putting it into practice... Well, simple solvers that only resolve local contact constraints are very beneficial as they give you straightforward implementations that have nice computational scaling. You resolve your contact constraints for all of the points of contact that your collision detector computed, and you get on with the next frame or timestep.
I have not done much exhaustive brainstorming around this, though I have pondered this on and off over the years.
Consider that if you want to actually compute the correct result of resolving forces imposed on stacks or jumbles or piles of objects, it stands to reason that the entire connectivity graph based on contact for that pile of objects must be resolved each frame. This is not in and of itself prohibitive, given that the collision detection would already be required anyway for the basic contact solver, so you just additionally compute the connectivity graph based on the contacting objects.
Then, this graph (one for each group of contacting objects) must be used each step of the simulation to implement and solve a global constraint system (one for each group of contacting objects) to resolve the interacting forces.
It’s clear that the computational complexity scales rather dramatically with the dimension (given that the cardinality of this graph will likely scale with the size of the pile squared or cubed depending on if you’re doing 2D or 3D physics) and further modulated by whatever solver you are using to solve the contact constraint model that you have to have to keep things behaving as they should. Even though the runtime of such models might be subquadratic, remember that the model is based on the contact connectivity, which would be subject to almost arbitrary change each timestep, meaning that the time complexity of setting up the model in order to run it is also very much implicated here. Likely quadratic or worse, with worse yet constants.
In almost all game-type simulations, it is difficult to make a case for the importance of having a high fidelity model of just how a jumble or pile or stack of objects (which would themselves exhibit idealized material properties, e.g. that objects are totally rigid to begin with, reducing the resultant realism anyway) reacts to arbitrary disruption. The result is always highly chaotic, and the computation involved to evaluate the most correct response is always going to be rather prohibitive: A massive constraint system must be modeled out of the connectivity graph for that particular physics state, and then every force impulse that is applied to any object on that group of objects is accumulated for input to that model and evaluated. The good news is that the model can be evaluated one time per timestep, but the overhead of having to initialize the model for every single timestep (as the scene evolves) is most likely to be the performance killer.
So in conclusion, although I referred to the pairwise contact solver model as "basic", it is usually the best fit for the problem. You're really opening up a can of worms trying to move beyond it. But I do believe that under special circumstances given a particular type of scene, a more sophisticated and flexible solver can make great strides on something like stack stability.
But you'd be about as well served pursuing heuristical solutions (such as putting simulated objects to sleep etc) or design changes to address the issue.