The Good, the Bad, and the IMPL-EX
A practical guide to the integration scheme that saved my simulation
1. The 3 AM crash: why your contact simulation diverges
You've set up your model. Meshed the geometry, calibrated the materials, written the input file. You launch the analysis, go to sleep, and at 3 AM the solver quits: "Newton–Raphson failed to converge."
If you work in earthquake engineering, you've been there. The culprit is almost always the same: contact, softening, or both at once. The tangent stiffness matrix goes singular, the iterative solver chokes, and your 20-second earthquake record dies at second 8.
This post is about a scheme that fixes this. It's called IMPL-EX (IMPLicit-EXplicit), it was introduced by Oliver et al. (2008), and it deserves to be more widely known among earthquake engineers. Let me walk you through the math.
2. The Bad: implicit integration and Newton–Raphson
In the standard implicit approach, we solve a nonlinear system of equations at every load (or time) step. Given external forces \(\mathbf{f}_{\mathrm{ext}}\) and a displacement vector \(\mathbf{u}\), equilibrium demands:
\[\mathbf{R}(\mathbf{u}) \;=\; \mathbf{f}_{\mathrm{int}}(\mathbf{u}) - \mathbf{f}_{\mathrm{ext}} \;=\; \mathbf{0}\]We solve this with Newton–Raphson iteration. At each iteration \(k\):
\[\mathbf{K}_T^{(k)} \, \Delta\mathbf{u}^{(k)} = -\mathbf{R}^{(k)} \qquad\text{where}\quad \mathbf{K}_T = \frac{\partial \mathbf{f}_{\mathrm{int}}}{\partial \mathbf{u}}\]Update \(\mathbf{u}^{(k+1)} = \mathbf{u}^{(k)} + \Delta\mathbf{u}^{(k)}\), repeat until \(\|\mathbf{R}\| < \text{tol}\).
This works beautifully—as long as \(\mathbf{K}_T\) stays well-conditioned. But consider what happens at a contact interface: when two surfaces collide, the penalty stiffness jumps from zero to \(10^8\) in one step. Or in a softening material (masonry, concrete), the tangent goes negative. In either case, the condition number of \(\mathbf{K}_T\) collapses, and Newton–Raphson diverges.
The cruel irony: the implicit method gives up precisely at the moments that matter most in earthquake engineering—impact, cracking, crushing.
Why does the tangent go singular at contact? (expand)
Consider a penalty-based normal contact law: \(f_N = k_N \, g_N\) when \(g_N < 0\) (penetration), and \(f_N = 0\) when \(g_N \geq 0\) (gap open). The tangent \(\partial f_N / \partial g_N\) jumps discontinuously between \(k_N\) and \(0\) at the contact boundary. This discontinuity means the linearization used in Newton–Raphson is not valid in a neighborhood of the transition—the solver oscillates between open and closed states without converging.
Worse: if contact opens and closes repeatedly within a single step (chatter), the tangent matrix changes at every iteration, and the convergence radius shrinks to zero.
3. The Ugly: explicit integration
The alternative is explicit time integration (e.g., central differences). No system of equations to solve—just a forward march:
\[\mathbf{u}_{n+1} = 2\mathbf{u}_n - \mathbf{u}_{n-1} + \Delta t^2 \, \mathbf{M}^{-1} \left(\mathbf{f}_{\mathrm{ext},n} - \mathbf{f}_{\mathrm{int},n}\right)\]No tangent matrix, no iterations, no convergence failures. So what's the problem?
Conditional stability. The critical time step is:
\[\Delta t_{\mathrm{crit}} = \frac{2}{\omega_{\max}} \approx \frac{h_{\min}}{c}\]where \(h_{\min}\) is the smallest element size and \(c\) is the wave speed in the material. For a masonry building with 20 mm mesh elements and \(c \approx 2000\) m/s, this gives \(\Delta t_{\mathrm{crit}} \approx 10^{-5}\) s. A 20-second earthquake requires 2 million steps. Each step is cheap, but two million of them on a 3D building model? Impractical.
| Property | Implicit (NR) | Explicit |
|---|---|---|
| System solve per step | Yes (iterative) | No (lumped mass) |
| Stability | Unconditional | Conditional (\(\Delta t < \Delta t_{\mathrm{crit}}\)) |
| Typical \(\Delta t\) | \(10^{-3}\)–\(10^{-2}\) s | \(10^{-6}\)–\(10^{-5}\) s |
| Contact / softening | Diverges | Handles naturally |
| Large 3D models | Feasible (few steps) | Impractical (millions of steps) |
We're stuck. Implicit can't handle the physics, explicit can't handle the cost. Is there a third option?
4. The Good: IMPL-EX — the idea
The key insight of Oliver et al. (2008) is simple once you see it:
The IMPL-EX idea: Don't solve the nonlinear constitutive equations during the global equilibrium iteration. Instead, extrapolate the internal variables from previous converged steps, solve a linear system, and correct the constitutive response after convergence.
Let's make this precise. In a standard implicit scheme, the stress at a Gauss point depends on the current strain through the constitutive law:
\[\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}\!\left(\boldsymbol{\varepsilon}_{n+1},\, \boldsymbol{\alpha}_{n+1}\right)\]where \(\boldsymbol{\alpha}\) represents internal variables (damage, plastic strain, hardening variables, etc.). The problem is that \(\boldsymbol{\alpha}_{n+1}\) depends on \(\boldsymbol{\varepsilon}_{n+1}\), which depends on \(\mathbf{u}_{n+1}\), which is what we're trying to find. This circular dependency is what makes the problem nonlinear.
IMPL-EX breaks the circle. Instead of solving for \(\boldsymbol{\alpha}_{n+1}\) simultaneously, it extrapolates:
\[\widetilde{\boldsymbol{\alpha}}_{n+1} = \boldsymbol{\alpha}_n + \frac{\Delta t_{n+1}}{\Delta t_n} \left(\boldsymbol{\alpha}_n - \boldsymbol{\alpha}_{n-1}\right)\]For equal time steps, this simplifies to \(\widetilde{\boldsymbol{\alpha}}_{n+1} = 2\boldsymbol{\alpha}_n - \boldsymbol{\alpha}_{n-1}\).
With \(\widetilde{\boldsymbol{\alpha}}_{n+1}\) fixed (known from previous steps), the stress becomes a linear function of strain. The tangent matrix is constant, symmetric, and positive semi-definite. Newton–Raphson converges in exactly one iteration.
After the linear solve gives us \(\mathbf{u}_{n+1}\), we go back and evaluate the true constitutive response to get the corrected \(\boldsymbol{\alpha}_{n+1}\) and \(\boldsymbol{\sigma}_{n+1}\). These corrected values sit on the real constitutive curve and are stored for the next extrapolation.
5. The algorithm, step by step
Algorithm: IMPL-EX integration at step \(n+1\)
Given: converged states at steps \(n\) and \(n-1\) with internal variables \(\boldsymbol{\alpha}_n\), \(\boldsymbol{\alpha}_{n-1}\) and secant stiffness \(\kappa_n\).
- Extrapolate. Compute the predicted internal variables: \[\widetilde{\boldsymbol{\alpha}}_{n+1} = 2\boldsymbol{\alpha}_n - \boldsymbol{\alpha}_{n-1}\] This defines an extrapolated tangent (or secant) stiffness \(\widetilde{\mathbf{K}}_{n+1}\) that is constant—it does not depend on the unknown \(\mathbf{u}_{n+1}\).
- Linear solve. Solve the linearized equilibrium: \[\widetilde{\mathbf{K}}_{n+1} \, \Delta\mathbf{u} = \mathbf{f}_{\mathrm{ext},n+1} - \mathbf{f}_{\mathrm{int}}\!\left(\mathbf{u}_n\right)\] This is a single linear system—one factorization, one back-substitution, done. No iterations.
- Correct. With \(\mathbf{u}_{n+1} = \mathbf{u}_n + \Delta\mathbf{u}\) now known, evaluate the true constitutive response: \[\boldsymbol{\sigma}_{n+1} = \boldsymbol{\sigma}\!\left(\boldsymbol{\varepsilon}_{n+1},\, \boldsymbol{\alpha}_{n+1}^{\text{true}}\right)\] Store \(\boldsymbol{\alpha}_{n+1}^{\text{true}}\) for the next extrapolation. The corrected stress sits on the real constitutive curve.
The beauty is in step 2. Because the tangent is constant:
- It is always symmetric positive semi-definite (no negative tangents from softening).
- Newton–Raphson converges in one iteration by construction.
- There is no chatter, no divergence, no bisection of load steps.
The price is that equilibrium is only approximately satisfied at each step—there is a residual of order \(O(\Delta t)\). But the corrected internal variables are exact (on the constitutive curve), and the error decreases as \(\Delta t \to 0\).
Accuracy order and error bound (expand)
Oliver et al. (2008) prove that the IMPL-EX scheme is first-order accurate in \(\Delta t\), the same order as backward Euler. The global error in displacement satisfies:
\[\|\mathbf{u}^{\text{impl-ex}} - \mathbf{u}^{\text{exact}}\| \leq C \, \Delta t\]where \(C\) depends on the smoothness of the loading and the constitutive response. In practice, the error is often smaller than this bound because the extrapolation is good when the internal variables evolve smoothly—which they do for most of the analysis.
The scheme is not unconditionally stable. There is a critical \(\Delta t\) beyond which the extrapolation overshoots and the scheme diverges. However, this critical step is typically orders of magnitude larger than the explicit \(\Delta t_{\mathrm{crit}}\), because it depends on the rate of change of internal variables rather than the wave speed.
6. Toy model: seeing it in action
Theory is one thing. Let's watch it work.
Consider a system of two springs in series, fixed at one end, with a prescribed displacement \(u_{\mathrm{app}}\) at the other. Spring 1 has a nonlinear constitutive law with softening:
\[\sigma_1(\varepsilon) = 150\,\varepsilon\, e^{-0.5\varepsilon^2} + 60\, e^{-3(\varepsilon - 3)^2}\]Spring 2 is linear: \(\sigma_2 = K_2 \, (u_{\mathrm{app}} - u_1)\), with \(K_2 = 100\). Equilibrium at the internal node requires \(\sigma_1(u_1) = K_2(u_{\mathrm{app}} - u_1)\).
Watch what happens at each step:
- Implicit: Newton–Raphson iterates toward equilibrium. In the softening zone the tangent \(d\sigma_1/d\varepsilon\) becomes small, so convergence slows—sometimes requiring 5+ iterations per step. It gets there, but expensively.
- Explicit: one step using the tangent from the previous converged state. Fast, but the solution drifts off the equilibrium path because there's no correction. The residual accumulates.
- IMPL-EX: extrapolates the secant stiffness \(\tilde{\kappa}\), solves a single linear equation, then corrects back onto the constitutive curve. One "iteration" per step, and the solution tracks the exact path closely.
The IMPL-EX math for the toy model (expand)
In the toy model, the "internal variable" is the secant stiffness:
\[\kappa_n = \frac{\sigma_1(u_{1,n})}{u_{1,n}}\]Step 1 — Extrapolate:
\[\widetilde{\kappa}_{n+1} = 2\kappa_n - \kappa_{n-1}\]Step 2 — Linear solve: equilibrium with the extrapolated stiffness is:
\[\widetilde{\kappa}_{n+1} \cdot u_1 = K_2 (u_{\mathrm{app}} - u_1) \quad\Longrightarrow\quad u_1 = \frac{K_2 \, u_{\mathrm{app}}}{\widetilde{\kappa}_{n+1} + K_2}\]Step 3 — Correct:
\[\sigma_{1,n+1} = \sigma_1(u_1) \qquad \kappa_{n+1} = \frac{\sigma_{1,n+1}}{u_1}\]The corrected point \((u_1, \sigma_{1,n+1})\) lies exactly on the constitutive curve. The displacement \(u_1\) is approximate (linear solve), but the stress is exact for that displacement.
7. IMPL-EX for contact
The toy model demonstrates IMPL-EX for material nonlinearity. But what about contact—the other major source of convergence headaches?
In a penalty-based contact formulation, the contact traction depends on the gap function \(g_N\) and a contact projector \(p_c\) that indicates whether the surfaces are in contact:
\[f_N = \begin{cases} k_N \, |g_N| & \text{if } g_N < 0 \;\;(\text{penetration}) \\ 0 & \text{if } g_N \geq 0 \;\;(\text{open}) \end{cases}\]The projector \(p_c\) is the binary indicator: \(p_c = 1\) if in contact, \(p_c = 0\) if open. In the standard implicit approach, \(p_c\) is evaluated at every iteration, causing the tangent matrix to change structure mid-iteration—this is the origin of contact chatter.
The IMPL-EX treatment
The internal variables for contact are the gap \(g_N\) (and tangential slip \(g_T\) for friction). We extrapolate these:
\[\widetilde{g}_{N,n+1} = 2\,g_{N,n} - g_{N,n-1}\]But there's a subtlety: the contact projector \(p_c\) is not extrapolated. It is binary-valued (0 or 1), and linear extrapolation of a discrete indicator is meaningless. Instead, \(p_c\) is frozen from the last converged step:
\[\widetilde{p}_{c,n+1} = p_{c,n}\]Why freeze \(p_c\)? If contact was active at step \(n\), the IMPL-EX scheme assumes it remains active during the linear solve at step \(n+1\). The correction step then evaluates the true gap and updates \(p_c\) for the next step. This eliminates chatter entirely: the contact state can only change between steps, not within iterations.
With the frozen projector and extrapolated gap, the contact stiffness contribution is constant during the solve. The global tangent matrix remains symmetric and well-conditioned. After the solve, the true gap is computed, the true contact forces are evaluated, and the corrected state is stored.
Friction: Coulomb with IMPL-EX (expand)
For Coulomb friction, the tangential traction is:
\[f_T = \begin{cases} k_T \, g_T & \text{if } |f_T| \leq \mu \, f_N \;\;(\text{stick}) \\ \mu \, f_N \, \mathrm{sign}(g_T) & \text{if } |f_T| > \mu \, f_N \;\;(\text{slip}) \end{cases}\]The stick/slip transition is another discrete state change that causes chatter in implicit schemes. In IMPL-EX, the slip state is frozen from the previous converged step, and the tangential slip \(g_T\) is extrapolated. The friction contribution to the tangent is therefore constant during the solve.
8. When NOT to use IMPL-EX
IMPL-EX is not a silver bullet. Intellectual honesty requires listing its limitations:
- First-order accuracy. The global error is \(O(\Delta t)\). If you need high-accuracy results (energy balance, long-duration creep), you may need very small time steps, partially defeating the purpose.
- Conditional stability. There is a critical \(\Delta t\). It's much larger than the explicit one, but it exists and requires a sensitivity study. If the internal variables evolve rapidly (sudden brittle failure), the extrapolation can overshoot.
- Energy drift. Because equilibrium is only approximately satisfied, there is a small energy error at each step. Over many cycles (e.g., long earthquake records), this can accumulate. Monitor the energy balance.
- Not needed when implicit works. If your problem converges fine with Newton–Raphson (e.g., no contact, mild nonlinearity), standard implicit is more accurate at the same \(\Delta t\). Don't fix what isn't broken.
- Not a replacement for explicit in wave propagation. If you genuinely need to resolve stress waves (blast, impact at microsecond scale), explicit is the right tool.
| Property | Implicit (NR) | Explicit | IMPL-EX |
|---|---|---|---|
| Iterations per step | 3–20+ | 0 | 1 |
| Tangent matrix | Recomputed each iter | Not needed (lumped \(\mathbf{M}\)) | Constant per step |
| Stability | Unconditional | Conditional (wave speed) | Conditional (internal var. rate) |
| Accuracy | Exact equilibrium | 1st order | 1st order |
| Contact / softening | Diverges | OK | OK |
| Typical \(\Delta t\) | \(10^{-3}\) s | \(10^{-6}\) s | \(10^{-3}\) s |
| Best for | Smooth nonlinearity | Wave propagation, small models | Contact, softening, large models |
9. A real example: seismic pounding in masonry
To give you a sense of what IMPL-EX enables in practice: we used it to simulate the SERA-AIMS shaking table test—two half-scale rubble stone masonry building units separated by a dry joint, subjected to the Montenegro 1979 earthquake record.
The model includes a damage-plasticity material (ASDConcrete3D) for the masonry and penalty-based contact elements at the joint. Both the material and the contact use IMPL-EX integration.
The result: 20 seconds of earthquake simulation in ~90 minutes, with \(\Delta t = 0.00125\) s, a single iteration per step, and zero convergence failures. A standard implicit scheme failed to converge beyond 8 seconds with the same time step.
10. Where to go from here
If you want to try IMPL-EX yourself:
- The formulation is implemented in OpenSees via the ASDConcrete3D material and a custom contact element. It's open source.
- The original paper by Oliver et al. (2008) is the definitive reference. Read Section 3 for the algorithm and Section 5 for the convergence proof.
- For contact applications, see Oliver et al. (2009) which extends IMPL-EX to frictional contact.
The scheme is simple to implement—if you already have a standard implicit code, the changes are localized to the constitutive update routine. You keep your global solver, your element library, your post-processing. You just change how the stress is computed at each Gauss point.
For earthquake engineers: this is a practical tool. Not a theoretical curiosity, not a niche method for one specific problem. It works for contact, for softening, for cyclic loading, for dynamics. It let me sleep through the night while my simulation ran to completion. Maybe it can do the same for you.
- Oliver, J., Huespe, A.E., and Cante, J.C. (2008). "An implicit/explicit integration scheme to increase computability of non-linear material and contact/friction problems." Computer Methods in Applied Mechanics and Engineering, 197(21–24), 1865–1889. doi:10.1016/j.cma.2007.09.027
- Oliver, J., Huespe, A.E., and Cante, J.C. (2009). "Implicit/explicit integration of coupled thermo-mechanical contact problems." International Journal for Numerical Methods in Engineering. doi:10.1002/nme.2806