Interactive Lecture
Newton-Raphson, Load Stepping, and Advanced Solution Methods
by Dr Arnaldo Delli Carri
Nonlinear problems arise from three main sources, often occurring simultaneously:
1. Material Nonlinearity: Stress-strain relationship is nonlinear
\[\sigma = f(\varepsilon) \quad \text{(e.g., plasticity, hyperelasticity)}\]
2. Geometric Nonlinearity: Large displacements change equilibrium equations
\[\text{Equilibrium on deformed configuration}\]
3. Contact/Boundary Nonlinearity: Changing boundary conditions
\[\text{Contact pressure active only when } g(u) \leq 0\]
Nonlinear problems cannot be solved directly like linear systems. Instead, we linearise about a current state and iterate to convergence. The choice of solution algorithm — Newton-Raphson, modified Newton, arc-length — depends on the specific nonlinearity type.
The standard approach for solving nonlinear FEM equations:
Residual Form: Find \(\{u\}\) such that
\[\{R\}(\{u\}) = \{F_{ext}\} - \{F_{int}\}(\{u\}) = \{0\}\]
where \(\{F_{int}\}\) are internal forces depending on displacement.
Newton-Raphson Iteration:
\[[K_T]^i \{\Delta u\}^{i+1} = \{R\}^i\]
\[\{u\}^{i+1} = \{u\}^i + \{\Delta u\}^{i+1}\]
where \([K_T]\) is the tangent stiffness matrix:
\[[K_T] = \frac{\partial \{R\}}{\partial \{u\}} = -\frac{\partial \{F_{int}\}}{\partial \{u\}}\]
The tangent stiffness represents the instantaneous stiffness at the current state:
For material nonlinearity:
\[[K_T] = \int_\Omega [B]^T[D_{tan}][B]d\Omega\]
where \([D_{tan}] = \frac{d\sigma}{d\varepsilon}\) is the tangent modulus
For geometric nonlinearity:
\[[K_T] = [K_L] + [K_\sigma]\]
\([K_L]\) = linear stiffness, \([K_\sigma]\) = stress stiffness (geometric stiffness)
Multiple criteria ensure solution accuracy:
1. Residual Norm (Force Equilibrium):
\[\frac{\|\{R\}^i\|}{\|\{F_{ext}\}\|} < \epsilon_F \quad (\text{typically } 10^{-3} \text{ to } 10^{-6})\]
2. Displacement Norm (Solution Change):
\[\frac{\|\{\Delta u\}^i\|}{\|\{u\}^i\|} < \epsilon_u \quad (\text{typically } 10^{-4} \text{ to } 10^{-6})\]
3. Energy Norm:
\[\frac{\{\Delta u\}^T \{R\}}{\{u\}^T \{F_{ext}\}} < \epsilon_E\]
For problems with complex nonlinear response:
Incremental-Iterative Method:
Adaptive Load Stepping:
For problems with snap-through or snap-back behaviour:
Constraint Equation:
\[\{\Delta u\}^T \{\Delta u\} + \psi^2 \Delta\lambda^2 \|\{F_{ref}\}\|^2 = \Delta s^2\]
where:
This allows the solution to follow the equilibrium path through limit points where load-displacement curves have vertical or negative tangents.
Implicit Dynamic Analysis: Newmark method with Newton-Raphson at each time step
\[[M]\{\ddot{u}\}_{n+1} + [C]\{\dot{u}\}_{n+1} + \{F_{int}\}(\{u\}_{n+1}) = \{F_{ext}\}_{n+1}\]
Predictor-corrector scheme:
Watch the Newton-Raphson method solve a nonlinear spring problem. The simulation shows how the algorithm iterates to find equilibrium by linearising about the current state.
Nonlinear FEM solvers show up wherever a structure leaves the linear-elastic, small-displacement world:
| Concept | Equation | Description |
|---|---|---|
| Residual | \(\{R\} = \{F_{ext}\} - \{F_{int}\}(\{u\})\) | Out-of-balance force |
| Newton-Raphson Update | \([K_T]^i \{\Delta u\}^{i+1} = \{R\}^i\) | Linearised equilibrium |
| Tangent Stiffness | \([K_T] = -\frac{\partial \{F_{int}\}}{\partial \{u\}}\) | Instantaneous stiffness |
| Residual Convergence | \(\frac{\|\{R\}\|}{\|\{F_{ext}\}\|} < \epsilon_F\) | Force equilibrium criterion |
| Displacement Convergence | \(\frac{\|\{\Delta u\}\|}{\|\{u\}\|} < \epsilon_u\) | Solution change criterion |
| Arc-Length Constraint | \(\{\Delta u\}^T\{\Delta u\} + \psi^2\Delta\lambda^2\|\{F_{ref}\}\|^2 = \Delta s^2\) | Path-following method |
| Effective Stiffness (Dynamics) | \([K_T^{eff}] = [K_T] + a_0[M] + a_1[C]\) | Includes inertia and damping |
A nonlinear spring has force-displacement relation F = 1000u + 200u³ (N, m). An external force F = 150 N is applied. Perform three Newton-Raphson iterations starting from u⁰ = 0. Calculate the residual and tangent stiffness at each iteration.
Show Solution
Iteration 0: \(u^0 = 0\)
\[F_{int} = 1000(0) + 200(0)^3 = 0\]
\[R^0 = 150 - 0 = 150\,\text{N}\]
\[K_T = \frac{dF}{du} = 1000 + 600u^2 = 1000\,\text{N/m}\]
\[\Delta u^1 = \frac{R^0}{K_T} = \frac{150}{1000} = 0.15\,\text{m}\]
\[u^1 = 0 + 0.15 = 0.15\,\text{m}\]
Iteration 1: \(u^1 = 0.15\)
\[F_{int} = 1000(0.15) + 200(0.15)^3 = 150 + 0.675 = 150.675\,\text{N}\]
\[R^1 = 150 - 150.675 = -0.675\,\text{N}\]
\[K_T = 1000 + 600(0.15)^2 = 1000 + 13.5 = 1013.5\,\text{N/m}\]
\[\Delta u^2 = \frac{-0.675}{1013.5} = -0.000666\,\text{m}\]
\[u^2 = 0.15 - 0.000666 = 0.149334\,\text{m}\]
Iteration 2: \(u^2 = 0.149334\)
\[F_{int} = 1000(0.149334) + 200(0.149334)^3 = 149.334 + 0.666 = 150.0\,\text{N}\]
\[R^2 = 150 - 150.0 \approx 0\,\text{N}\]
After just three iterations \(|R| < 0.001\,\text{N}\): the residual falls from 150 to 0.675 to ~0, the signature of Newton-Raphson's quadratic convergence.
Answer: u ≈ 0.1493 m (converged after 3 iterations).
For the problem above, if the full load F = 150 N requires 3 iterations to converge, estimate how many total iterations would be needed using: (a) 1 load step (apply full load), (b) 3 equal load steps, (c) 5 equal load steps. Assume convergence requires 3-5 iterations per step.
Show Solution
(a) 1 load step (full load):
Iterations: 3 (as shown in Problem 1)
(b) 3 equal load steps (\(F = 50, 100, 150\,\text{N}\)):
Each smaller load increment typically converges faster (2-3 iterations).
\[\text{Estimated total} = 3 \times 3 = 9 \text{ iterations}\]
(c) 5 equal load steps (\(F = 30, 60, 90, 120, 150\,\text{N}\)):
Very small increments converge quickly (2 iterations each).
\[\text{Estimated total} = 5 \times 2 = 10 \text{ iterations}\]
Analysis:
For this mildly nonlinear problem, 1 load step is the most efficient (3 iterations). For severely nonlinear problems (plasticity, contact), multiple load steps still pay off:
Rule of thumb: use 5-10 load steps for plasticity, 20-50 for contact problems.
Answer: (a) 3 iterations, (b) ≈9 iterations, (c) ≈10 iterations — 1 load step wins here, but that reverses for strongly nonlinear problems.
For a nonlinear spring F = ku + αu³, compare the tangent stiffness K_T = dF/du with the secant stiffness K_S = F/u. For k = 1000, α = 100, and u = 0.2 m, calculate both values. Which would converge faster in Newton-Raphson iterations?
Show Solution
Calculate force at u = 0.2 m:
\[F = 1000(0.2) + 100(0.2)^3 = 200 + 0.8 = 200.8 \text{ N}\]
Tangent stiffness:
\[K_T = \frac{dF}{du} = k + 3\alpha u^2 = 1000 + 3(100)(0.2)^2\]
\[K_T = 1000 + 3(100)(0.04) = 1000 + 12 = 1012 \text{ N/m}\]
Secant stiffness:
\[K_S = \frac{F}{u} = \frac{200.8}{0.2} = 1004 \text{ N/m}\]
Comparison:
• Tangent stiffness (1012 N/m) is 0.8% higher than secant (1004 N/m)
• Tangent stiffness represents the local slope at current displacement
• Secant stiffness represents the average slope from origin
Convergence:
Newton-Raphson with tangent stiffness converges quadratically (very fast). Using secant stiffness (Modified Newton or Secant method) converges linearly (slower but avoids recalculating K_T every iteration).
Trade-off: Tangent requires computing K_T each iteration (expensive) but needs fewer iterations. For expensive K_T calculations (large FEM models), Modified Newton (reuse K_T for several iterations) is often more efficient overall.
Answer: K_T = 1012 N/m, K_S = 1004 N/m (0.8% higher); tangent stiffness converges faster (quadratically).