We discuss a solution algorithm for quasi-static elastoplastic problems with hardening. Such problems can be described by a time dependent variational inequality, where the displacement and the plastic strain fields serve as primal variables. After discretization in time, one variational inequality of the second kind is obtained per time step and can be reformulated as a minimization problem with a convex energy functional which depends smoothly on the displacement and nonsmoothly on the plastic strain. There exists an explicit formula for minimizing the energy functional with respect to the plastic strain for a given displacement. By substitution, the energy functional can be written as a functional depending only on the displacement. The theorem of Moreau from convex analysis states that the energy functional is differentiable with an explicitly computable first derivative. The second derivative of the energy functional does not exist, due to the lack of smoothness of the plastic strain across the elastoplastic interface, which separates the continuum in elastically and plastically deformed parts. A Newton-like method exploiting slanting functions of the energy functional's first derivative instead of the nonexistent second derivative is applied. Such a method is called a slant Newton method for short. The local superlinear convergence of the algorithm in the discrete case is shown, and sufficient regularity assumptions are formulated, which would guarantee the local superlinear convergence also in the continuous case.

The quasistatic, Prandtl--Reuss perfect plasticity at small strains is combined with a gradient, reversible (i.e., admitting healing) damage which influences both the elastic moduli and the yield stress. Existence of weak solutions of the resulting system of variational inequalities is proved by a suitable fractional-step discretization in time with guaranteed numerical stability and convergence. After finite-element approximation, this scheme is computationally implemented and illustrative two-dimensional simulations are performed. The model allows, e.g., for application in geophysical modeling of reoccurring rupture of lithospheric faults. Resulting incremental problems are solved in MATLAB by quasi-Newton method to resolve the elastoplasticity component of the solution, while the damage component is obtained by solving a quadratic programming problem.