3. Weak forms and variational principles

The strong form of a nonlinear boundary value problem—equilibrium in the deformed configuration, traction boundary conditions, and a constitutive law—is rarely solved directly. Instead, nonlinear finite element analysis reformulates it into a weak form, where the governing equations are satisfied in an integral (weighted-residual) sense rather than pointwise. This reformulation lowers the continuity requirements on the displacement field, naturally incorporates traction and body-force boundary conditions, and provides the algebraic foundation for Newton-type solution schemes. This chapter develops the virtual work principle, the total potential energy functional, and the consistent linearization that produces the tangent stiffness matrix.

From Strong Form to Weak Form

Consider a body occupying the reference configuration Ω₀ with boundary Γ₀ = Γ_u ∪ Γ_t. The strong form of static equilibrium in the reference configuration, written in terms of the first Piola–Kirchhoff stress P and the displacement field u, is:

Div P + ρ₀ b = 0        in Ω₀
u = ū                    on Γ_u   (Dirichlet)
P · N = t̄                on Γ_t   (Neumann)

Here ρ₀ is the reference density, b is the body force per unit mass, N is the outward unit normal in the reference configuration, and t̄ is the prescribed traction. The divergence operator acts on the reference coordinates X.

To obtain the weak form, multiply the equilibrium equation by an admissible virtual displacement δu that vanishes on Γ_u, integrate over Ω₀, and apply the divergence theorem. The result is the principle of virtual work:

∫_Ω₀ P : δF dV − ∫_Ω₀ ρ₀ b · δu dV − ∫_Γt t̄ · δu dA = 0   ∀ δu ∈ V₀

where F = I + ∂u/∂X is the deformation gradient and δF = ∂δu/∂X is its variation. The test space V₀ contains all admissible variations satisfying the homogeneous Dirichlet condition. This single scalar equation encodes equilibrium in all directions and is the starting point for every nonlinear element formulation.

Energy Functionals and Hyperelasticity

When the material is hyperelastic, the first Piola–Kirchhoff stress derives from a stored-energy density W(F):

P = ∂W/∂F

In this case the weak form is equivalent to the stationarity condition of the total potential energy functional:

Π(u) = ∫_Ω₀ W(F) dV − ∫_Ω₀ ρ₀ b · u dV − ∫_Γt t̄ · u dA

Setting the first variation δΠ = 0 recovers the virtual work equation exactly. This equivalence is powerful: it guarantees that the tangent operator derived from the second variation is symmetric for hyperelastic materials, and it provides a natural energy-based convergence measure for Newton iterations.

For a compressible neo-Hookean model, a common choice is:

W = (μ/2)(I₁ − 3) − μ ln J + (λ/2)(ln J)²

where I₁ = tr(FᵀF), J = det F, and μ, λ are the Lamé parameters. The corresponding stress and tangent moduli follow from systematic differentiation—a process that must be carried out consistently to preserve quadratic convergence.

Consistent Linearization

The weak form is a nonlinear equation in u. Newton's method requires its linearization about a known state u in the direction of an increment Δu. The linearized virtual work is:

D(δW)[Δu] = ∫_Ω₀ δF : A : ΔF dV + ∫_Ω₀ P : Δ(δF) dV

where A = ∂P/∂F is the fourth-order material tangent modulus, and Δ(δF) captures the geometric (initial-stress) contribution. The first term is the material stiffness; the second is the geometric stiffness. Together they form the element tangent stiffness matrix Kₜ.

The term "consistent" means that A is obtained by exact differentiation of the stress update, not by numerical approximation. An inconsistent tangent still converges, but typically degrades Newton's quadratic rate to linear—often doubling or tripling the number of iterations.

Workflow: From Variational Statement to Solver

The following sequence summarizes how the variational formulation feeds a nonlinear solution procedure. Each step corresponds to a concrete implementation task in a finite element code.

1. Choose the kinematic measure. Total Lagrangian (reference configuration, F, P) or updated Lagrangian (current configuration, Cauchy stress, spatial tangent). The choice affects which stress and tangent objects are stored.
2. Define the stored-energy or stress update. For hyperelasticity, provide W(F); for plasticity, provide a return-mapping algorithm that updates stress and internal variables.
3. Derive the consistent tangent. Differentiate the stress update analytically. For complex models, use automatic differentiation or a verified numerical tangent.
4. Assemble the residual and tangent. Loop over elements, compute Rₑ = ∫ Bᵀ P dV and Kₜ,ₑ = ∫ Bᵀ A B dV + geometric term, then scatter into global arrays.
5. Solve and update. Solve Kₜ Δu = −R, update u ← u + Δu, and check the energy or residual norm against a tolerance.

For a representative 3D neo-Hookean patch test with μ = 80 MPa, λ = 120 MPa, and a prescribed 5% stretch, the first Newton iteration typically reduces the residual norm by several orders of magnitude. If the reduction is only a factor of 2–3, the tangent is likely inconsistent or the load step is too large.

Visual Overview of the Variational Chain

The diagram below traces the logical dependencies from the strong form through the weak form and energy functional to the linearized tangent and Newton solution. It highlights where consistency is enforced and where the material and geometric contributions enter.

Practical Verification of the Tangent

Before trusting a nonlinear element, verify the consistent tangent numerically. The procedure is straightforward and catches the majority of implementation errors:

1. Select a single element with a random but admissible deformation gradient F.
2. Compute the stress P(F) from the constitutive routine.
3. Perturb each component Fᵢⱼ by a small ε (e.g., 1×10⁻⁶) and recompute P.
4. Form the numerical tangent Aᵢⱼₖₗ ≈ (Pᵢⱼ(F + ε eₖₗ) − Pᵢⱼ(F)) / ε.
5. Compare against the analytical A returned by the material routine. The relative error should be below 1×10⁻⁵ for a correct implementation.

If the error is large, check the sign conventions, the ordering of tensor components, and whether the geometric term has been included in the element-level tangent. A frequent mistake is to omit the initial-stress contribution, which becomes dominant under large compressive stresses.

Key Takeaways

• The weak form replaces pointwise equilibrium with an integral statement, relaxing continuity requirements and embedding traction boundary conditions naturally.
• For hyperelastic materials, the weak form is the first variation of a total potential energy functional, guaranteeing symmetry of the tangent.
• Consistent linearization requires exact differentiation of the stress update; the resulting tangent splits into material and geometric contributions.
• Newton's quadratic convergence depends on tangent consistency—verify it numerically before production runs.
• The variational framework is independent of the specific constitutive model; only the stress update and its tangent change when switching materials.