Getting Your Head Around Large Deformation Plasticity
I spent about three months debugging a simulation that kept diverging at strains above 20 percent before I realized I was using small-strain plasticity logic on a problem that clearly wasn't small. The code was giving physically reasonable results for everything up to about 5 percent strain, then suddenly the Newton iterations started oscillating and the energy norm went off the rails. That's usually when you know your kinematic assumptions have outlived their usefulness. The core issue in finite strain elasto-plasticity isn't really the plasticity itself. It's figuring out what "strain" even means when the body has rotated and stretched enough that your reference configuration looks nothing like the current one. In small strain theory, you just take the symmetric part of the displacement gradient and you're done. In finite strain, you have to decide which measure of deformation you're actually using, and that decision ripples through every single equation in your formulation.
Introduction To Finite Strain Theory For Continuum Elasto Plasticity
The standard starting point is the multiplicative decomposition of the deformation gradient, which is F = Fe · Fp. This was introduced by Lee in 1969 and it's still the workhorse. The idea is clean enough: you split the total deformation into an elastic part that stores energy and a plastic part that represents permanent rearrangement. The elastic part goes into a hyperelastic stress-strain law, and the plastic part evolves according to a yield condition and flow rule. But the devil is in how you define Fe and Fp when the rotations get large. If you just multiply matrices and call it a day, you'll get frame-indifference violations that will quietly corrupt your results. The elastic strain measure you feed into your hyperelastic potential needs to be defined in a way that doesn't depend on how the body has rotated. That's why people typically use something like the right Cauchy-Green elastic tensor Ce = Fe^T · Fe, or equivalently the elastic right Cauchy-Green deformation based on the polar decomposition. Here's where most implementations trip up. You need to update Fp incrementally. At each time step, you compute a trial elastic deformation, check whether the yield function is exceeded, and if it is, you perform a plastic corrector. The corrector step requires solving a nonlinear equation because the flow rule is typically nonlinear in the stress space. In practice, people use a return-mapping algorithm, which is essentially a local Newton iteration on the stress update. For isotropic J2 plasticity, this can be solved semi-analytically, but for anisotropic yield functions or more complex hardening laws, you're doing full numerical integration.
I ran into a specific problem a few years ago with a rubber-metal adhesive joint simulation. The rubber was modeled with a hyperelastic Neo-Hookean material and the metal with finite strain plasticity. The issue was that the plastic spin, which is the skew-symmetric part of the plastic velocity gradient Lp = p · Fp^(-1), was accumulating rotation in a way that didn't match the physical behavior. The metal was essentially rotating itself unnecessarily because of how I had formulated the intermediate configuration. The workaround was to switch from the rate-form formulation to an updated Lagrangian approach where I tracked the plastic deformation relative to the last converged configuration rather than the initial reference. It added some complexity to the Jacobian but eliminated the spurious rotation drift entirely. Another thing nobody warns you about: the choice of stress measure matters enormously for convergence. The 2nd Piola-Kirchhoff stress is energy-conjugate to the Green-Lagrange strain, which makes it natural for a total Lagrangian formulation. But it can develop unrealistically large components when rotations are big, even if the actual physical stresses are modest. The Kirchhoff stress, which is = J · where is the true Cauchy stress, tends to behave better numerically in large rotation scenarios. I usually convert to Kirchhoff for the yield function evaluation and then back-convert for the global residual assembly. The consistent tangent operator is non-negotiable if you want quadratic convergence in the global Newton iteration. A lot of people skip deriving it properly and just use an approximate tangent, which might still converge but at a glacial pace. For finite strain plasticity with associative flow, the consistent tangent has additional terms coming from the plastic correction that don't appear in the small strain version. These terms involve the derivative of the plastic multiplier with respect to the trial stress, and they're essential for maintaining the symmetry of the global tangent in many formulations.
Get the Full Details

If you're implementing this from scratch, I'd recommend starting with a single element patch test in 2D plane strain before you go anywhere near a full structure. Use a simple bilinear quadrant element with a known analytical solution for large strain uniaxial tension with plasticity. If your element can't reproduce that, nothing else you build on top of it will work reliably either. I've seen people spend weeks debugging full models only to find the root cause was a sign error in the plastic flow direction that would have shown up in five minutes with a proper patch test. There are also commercial codes that handle this well if you don't want to write your own. Abaqus has a robust finite strain plasticity framework with several material models built in. Ansys and COMSOL have similar capabilities. The trade-off is that you lose some control over the constitutive integration, which matters if you're doing something unusual like strain-rate dependent plasticity at extreme strains or coupled damage-plasticity models. The main limitation of the multiplicative decomposition approach is that it assumes a clear separation between elastic and plastic deformation mechanisms. When you get into very large strains where the material microstructure changes significantly, or when damage and plasticity are tightly coupled, the concept of an intermediate configuration starts to break down. In those cases, some researchers have proposed additive decompositions of strain measures or other alternatives, but none of them have achieved the same level of maturity as Lee's framework. If you're working in a regime where the multiplicative split is questionable, be honest about it and document your assumptions rather than blindly trusting the output.
One practical tip that might save you some headaches: always check your trace conservation. In a pure plastic deformation with no volume change, the determinant of Fp should remain exactly 1. If you're seeing drift in det(Fp) over many load steps, your plastic update algorithm has a numerical issue, usually related to the implicit integration of the plastic flow. Reducing the increment size might mask it temporarily, but the real fix is usually in the consistency condition of your return map.