Staggered Solution Strategy for Coupled PDEs

Learning objectives

By the end of this lesson, you should be able to:

  • Explain why displacement and damage are coupled.
  • Distinguish monolithic and staggered solution strategies.
  • Write the staggered algorithm for phase-field fracture.
  • Understand load stepping, history-field updates, and irreversibility enforcement.

1. Why the system is coupled

Phase-field fracture solves for two fields:

\[ u,\qquad d. \]

The coupling is two-way:

u → strain → stress/energy → damage growth

d → stiffness degradation → stress redistribution

The displacement equation depends on \(d\) through the degraded stress, while the damage equation depends on \(u\) through the crack-driving energy.

2. Coupled residual form

The mechanical residual is

\[ R_u(u,d;v)= \int_\Omega \sigma(u,d):\varepsilon(v)\,d\Omega - \int_\Omega b\cdot v\,d\Omega - \int_{\Gamma_N}\bar t\cdot v\,dS. \]

The damage residual is

\[ R_d(u,d;q)= \int_\Omega g'(d)\psi_e^+(\varepsilon(u))q\,d\Omega + \int_\Omega G_c\left(\frac{d}{\ell}q+ \ell\nabla d\cdot\nabla q\right)d\Omega. \]

The coupled problem is

\[ R_u=0,\qquad R_d=0, \]

with irreversibility

\[ d_n\ge d_{n-1}. \]

3. Monolithic method

A monolithic method solves \(u\) and \(d\) together. Define

\[ w=\begin{bmatrix}u\\d\end{bmatrix}. \]

Then solve

\[ R(w)=0. \]

Newton's method gives

\[ J(w^k)\Delta w^k=-R(w^k), \qquad w^{k+1}=w^k+\Delta w^k. \]

This is strongly coupled and often efficient near fracture, but it is harder to implement and requires solving larger nonlinear systems.

4. Staggered method

A staggered method alternates between the two subproblems:

fix damage d → solve displacement u
fix displacement u → solve damage d
repeat until convergence

Advantages:

  • simpler implementation,
  • smaller linear systems,
  • easier debugging,
  • convenient for FEniCS-style codes.

Disadvantages:

  • slower convergence near crack growth,
  • sensitivity to load step size,
  • possible damage jumps if not carefully controlled.

5. Basic staggered algorithm

At load step \(n\), given \(d_{n-1}\):

\[ d_n^0=d_{n-1}. \]

For staggered iteration \(k=0,1,2,\dots\):

Mechanics solve

Given \(d_n^k\), find \(u_n^{k+1}\):

\[ R_u(u_n^{k+1},d_n^k;v)=0. \]

Damage solve

Given \(u_n^{k+1}\), find \(d_n^{k+1}\):

\[ R_d(u_n^{k+1},d_n^{k+1};q)=0. \]

Convergence check

Stop when

\[ \|u_n^{k+1}-u_n^k\|<\text{tol}_u, \qquad \|d_n^{k+1}-d_n^k\|<\text{tol}_d. \]

6. Load stepping

Fracture is usually solved incrementally:

\[ 0=\lambda_0<\lambda_1<\cdots<\lambda_N. \]

For displacement control,

\[ \bar u_n=\lambda_n\bar u_{\max}. \]

Small load steps are important near crack initiation and rapid crack propagation. If the load step is too large, the simulation may miss the peak load or create an unrealistic damage jump.

7. History-field update

A common way to approximate irreversibility is the history field:

\[ \mathcal{H}(x,t)=\max_{\tau\le t}\psi_e^+(\varepsilon(u(x,\tau))). \]

At load step \(n\), update

\[ \mathcal{H}_n(x)= \max\left(\mathcal{H}_{n-1}(x),\psi_e^+(\varepsilon(u_n(x)))\right). \]

The damage equation then uses \(\mathcal{H}\) instead of the current tensile energy.

8. Linear damage equation with history field

For \(g(d)=(1-d)^2+k\), the history-field damage equation becomes

\[ \left(\frac{G_c}{\ell}+2\mathcal{H}\right)d -G_c\ell\Delta d =2\mathcal{H}. \]

The weak form is

\[ \int_\Omega \left(\frac{G_c}{\ell}+2\mathcal{H}\right)dq\,d\Omega + \int_\Omega G_c\ell\nabla d\cdot\nabla q\,d\Omega = \int_\Omega 2\mathcal{H}q\,d\Omega. \]

If \(\mathcal{H}\) is fixed, this is a linear scalar elliptic problem.

9. Enforcing irreversibility

The discrete damage field should satisfy

\[ d_{n-1}\le d_n\le 1. \]

Common approaches are:

Method Main idea Comment
History field Prevents crack-driving force from decreasing Simple but not exact
Nodal clipping Set \(d_n=\max(d_n,d_{n-1})\) Easy but not variational
Penalty method Penalize healing Approximate
Bound-constrained solve Enforce \(d_{n-1}\le d_n\le1\) directly Most rigorous

For robust research simulations, bound-constrained minimization is preferred.

10. Compact pseudocode

for load step n:

    apply load λ_n
    d = d_old
    H = H_old

    for staggered iteration k:

        solve mechanics with fixed d

        H = max(H_old, tensile_energy(u))

        solve damage with fixed H

        enforce d_old <= d <= 1

        if ||d - d_previous|| < tolerance:
            break

    save u, d, H
    d_old = d
    H_old = H

The staggered strategy is the standard beginner-friendly route from weak forms to working phase-field fracture simulations.