4. Sensitivity analysis in nonlinear CAE

In linear structural optimization, the sensitivity of a response such as displacement or compliance with respect to a design variable can often be written in closed form because the stiffness matrix is independent of the displacement field. Nonlinear CAE breaks that convenience. Once the equilibrium state depends on the solution itself — through geometric nonlinearity, contact, material plasticity, or follower loads — the derivative of a response must account for the fact that the internal force vector, the tangent stiffness, and even the load vector may all change when a design variable changes. Sensitivity analysis is therefore not a post-processing step; it is a core solver capability that determines whether gradient-based optimization of nonlinear models is accurate, affordable, and robust.

What the derivative actually means in a nonlinear setting

Consider a discretized nonlinear equilibrium problem written as R(u, x) = F_ext(u, x) − F_int(u, x) = 0, where u is the converged displacement vector, x is a vector of design variables (shell thickness, material parameter, shape control point), and R is the residual. A scalar response g(u, x) — for example, peak von Mises stress, tip deflection, or buckling load factor — depends on x both directly and indirectly through u(x). The total derivative is

dg/dx = ∂g/∂x + (∂g/∂u) · (du/dx)

The term du/dx is the implicit part and is the source of all difficulty. Differentiating the equilibrium equation R(u(x), x) = 0 with respect to x gives the tangent relation

K_T · (du/dx) = −∂R/∂x ,   where K_T = ∂R/∂u

Here K_T is the converged tangent stiffness matrix at the final load step. Every sensitivity method is essentially a strategy for solving this equation without paying the cost of a full factorization for every design variable.

Direct differentiation

The direct method solves the tangent equation once per design variable. For each x_i, you form the pseudo-load vector −∂R/∂x_i, solve K_T · (du/dx_i) = −∂R/∂x_i, and then contract with ∂g/∂u to obtain dg/dx_i. The tangent matrix is factorized once and reused across all right-hand sides, so the marginal cost per variable is one back-substitution.

The pseudo-load is where nonlinearity bites. For a geometrically nonlinear model, ∂R/∂x_i includes derivatives of the internal force vector with respect to thickness or shape, which involve the strain-displacement matrix and the second Piola–Kirchhoff stress. For contact, it includes derivatives of the contact gap and normal direction. For plasticity, it includes derivatives of the consistent tangent and the yield condition. In practice these terms are assembled element-by-element, exactly as the residual is assembled, using the same quadrature points and the same converged state variables.

Direct differentiation is attractive when the number of design variables is small relative to the number of responses — for example, optimizing five thickness zones against hundreds of stress constraints. Its weakness is that the cost grows linearly with the number of design variables, and each additional variable requires a new pseudo-load assembly.

Adjoint methods

The adjoint method reverses the order of operations. Instead of solving for du/dx, it introduces an adjoint vector λ defined by

K_T^T · λ = (∂g/∂u)^T

Once λ is known, the sensitivity of a single response is

dg/dx = ∂g/∂x − λ^T · (∂R/∂x)

The key advantage is that λ depends only on the response g, not on the design variable. One adjoint solve therefore yields the complete gradient of one response with respect to all design variables. This is the method of choice when the design space is large — topology optimization with thousands of element densities, or shape optimization with hundreds of control points — and the number of responses is modest.

Two practical cautions apply. First, K_T is generally unsymmetric in nonlinear problems (contact, follower loads, non-associated plasticity), so the adjoint solve must use K_T^T, not K_T. Reusing a symmetric factorization silently produces wrong gradients. Second, the adjoint state must be computed at the converged load step; if the nonlinear solve is not fully converged, the adjoint gradient inherits that error and can point in the wrong direction.

Finite-difference verification

No sensitivity implementation should be trusted without verification. The standard check is a central finite difference:

dg/dx_i ≈ [ g(x + h·e_i) − g(x − h·e_i) ] / (2h)

where e_i is the unit vector for variable i and h is a perturbation size. The comparison is meaningful only if the perturbed analyses converge to the same load level with the same convergence tolerances and the same number of load steps. If the baseline and perturbed runs stop at different equilibrium states, the difference quotient measures solver inconsistency, not sensitivity.

Choosing h is a trade-off. Too large, and truncation error from the O(h²) term dominates; too small, and subtractive cancellation in floating-point arithmetic dominates. For thickness variables on the order of 1–10 mm, a relative perturbation of h = 10⁻³ · x_i is a reasonable starting point. For shape variables, perturb the control point by a small fraction of the local element size. A useful diagnostic is to sweep h over two or three decades and confirm that the relative error between the analytic and finite-difference gradient reaches a minimum and then rises — the classic V-shaped error curve.

Worked verification example

Consider a thin-walled cantilever modeled with shell elements, loaded by a transverse tip force large enough to produce moderate rotation. The design variables are the thicknesses of three zones: x = [2.0, 2.5, 3.0] mm. The response is tip deflection g = u_tip. The nonlinear solve uses a Newton–Raphson scheme with a residual tolerance of 1e-6 and five load steps.

Running the adjoint solver gives dg/dx = [−0.42, −0.31, −0.18] mm/mm. A central finite difference with h = 0.002 mm (0.1% of the nominal thickness) gives [−0.419, −0.309, −0.181]. The relative error is below 0.5% for all three variables, which is acceptable for gradient-based optimization. If the error exceeds a few percent, the first thing to inspect is the convergence of the perturbed runs — tightening the residual tolerance from 1e-6 to 1e-8 typically reduces the discrepancy by an order of magnitude.

Practical guidance for implementation

When implementing or selecting a sensitivity capability, keep the following points in mind. First, always verify against finite differences on a small subset of variables before trusting the full gradient. Second, monitor the adjoint residual; if it does not converge as tightly as the primal residual, the gradient will be inconsistent. Third, for contact-dominated problems, freeze the active set during the sensitivity computation and verify that the set does not change under the perturbation — if it does, the gradient is non-smooth and a subgradient or smoothed formulation is needed. Fourth, for path-dependent materials, the sensitivity must be computed at the final converged state with the consistent tangent, not the continuum tangent; using the wrong tangent is a frequent source of gradient error that finite differences will expose.

The payoff is direct: accurate nonlinear sensitivities enable gradient-based optimization to converge in far fewer design iterations than derivative-free methods, and they make it possible to handle design spaces with hundreds or thousands of variables that would be intractable with finite differences alone. The finite-difference check is not a formality — it is the contract that guarantees the optimizer is being driven by the physics, not by numerical noise.