I implemented an iterative constraint solver using impulses as described in Bender & Schmitt 2006:

http://www.researchgate.net/publication/228934234_Constraint-based_collision_and_contact_handling_using_impulses

Solving for the normal (penetration) constraint or the tangent constraint in isolation converges as expected, about 4-5 iterations before the penetration or tangent position is within the 'target' range (<0.0001). But running both constraints simultaneously causes 150+ iterations before they are both within target which sounds like too much for a single constraint, but maybe it is expected for a tight constraint?

The target range should be something that is numerical significant, but visually not. Try e.g. 0.005f if you use meters.

Oops, actually the target range was 0.0001 (edited the original post to fix this) which is 1/10 of the collision tolerance (0.001). Seems I can get away with 0.0004 (7 iterations) but after that it starts going up rather quickly (0.0003 -> 28 iterations, 0.0002 -> 70+ iterations). Hopefully that will be stable enough for stacking.