How to stabilize nonlinear solvers for rate-independent plasticity problems?

Plastic strain localization in a rate-independent limit is a physical phenomenon that demands robust and reliable nonlinear solves to be modeled in a computer. It is quite common on practice that standard iterative algorithms (e.g. Newton-Raphson), being applied to this demanding problem, lead to  convergence issues or even fail. From the physical viewpoint the problem can be attributed to the difference between the dilatation and friction angles for the pressure-sensitive materials, which gets worse as this difference gets larger.

Common remedies include deriving consistent Jacobian matrix, or switching to the equivalent rate-dependent visco-plastic formulation. Both methods have their specific side effects. A very high condition number of a heavily unsymmetrical Jacobian matrix renders it nearly useless in the context of an iterative linear solver, such as multigrid. Hence, the use of the Jacobian matrix is essentially limited to a 2D formulation, for which a direct solver is practical. The visco-plastic formulation is confusing form the conceptual viewpoint. It strives to achieve the convergence by modifying the physics of the problem. Hence, the stabilization viscosity is not a pure numerical parameter that can be freely selected, but it is a physical parameter that must be determined in the laboratory. The advantages of the visco-plastic formulation vanish, if rate-independent limit is considered, or if affordable grid size is (much) larger than the intrinsic localization length-scale. The latter condition is a dominant limiting factor for a 3D model. 

In this work we share a few recipes, that can potentially improve the convergence of the rate-independent plasticity problems, without relying on the availability of a direct solver, or perturbing the physics.