Phase-field fracture of a notched plate | PDDO-MultiPhysics Examples

1The Problem

A square plate with a horizontal notch running halfway in from the left is pulled at its top edge while its bottom edge is held. The phase field d, 0 in intact material and 1 in a crack, smears the crack over a regularization length; the elastic stress is degraded where d grows, and an irreversible history field makes sure the crack never heals.

$$\nabla\cdot\boldsymbol{\sigma} = \mathbf{0}, \qquad \boldsymbol{\sigma} = \big((1-d)^2 + k\big)\,\boldsymbol{\sigma}_0(\boldsymbol{\varepsilon})$$

$$G_c\Big(\frac{d}{\ell} - \ell\,\nabla^2 d\Big) = 2(1-d)\,\mathcal{H}, \qquad \mathcal{H} = \max_{s \le t} \psi^+(s), \quad \psi^+ = \tfrac12\,\boldsymbol{\sigma}_0:\boldsymbol{\varepsilon}$$

2The Problem File

The history is a field with an update rule, H: {type: scalar, update: "max(H, psi_plus)"}, applied after every converged step. The notch is not cut out of the geometry: it is a large initial history along y = 0.5, x < 0.5, so the damage there starts at one. The top displacement grows linearly with time, and Newton uses a line search.

The complete problem file, examples/fracture_phase_field.yaml

# Phase-field fracture (AT2) - notched plate in tension
# (plane-strain slab: thin, periodic in z). Irreversibility through the history field
#   H = max(H, psi_plus)        (update rule after every step)
# Damage d (stationary):   Gc (d/l - l lap(d)) = 2 (1 - d) H
# Elasticity:              div(sigma) = 0
# Isotropic degradation: sigma = ((1 - d)^2 + k) sigma_0, psi_plus = sigma_0 : eps / 2.
# A volumetric-deviatoric split (Amor: only tension and shear degrade) can be written with
# relations, see tests/test_internal_variables.py; on this notched plate its Newton solve
# still stalls near the crack faces (contact-like switch), even smoothed (C1), with a line
# search and adaptive load steps - an alternate-minimization (staggered u / d) solve is the
# usual remedy (TODO).
# The initial notch (x < 0.5 at y = 0.5) is given as a large initial history.
"k">name: fracture_phase_field
"k">domain:
  "k">size: [1.0, 1.0, 0.125]
  "k">n: [41, 41, 5]
  "k">periodic: [false, false, true]
"k">parameters:
  "k">E: 210.0
  "k">nu: 0.3
  "k">Gc: 2.7e-3
  "k">l0: 0.05                         # regularization length (2 nodes)
  "k">kres: 1.0e-4
  "k">K: "E/(3*(1 - 2*nu))"
  "k">G: "E/(2*(1 + nu))"
  "k">ubar: "1.0e-3*t"                 # top displacement
"k">fields:
  "k">u: vector
  "k">d: scalar
  "k">H: {type: scalar, update: "max(H, psi_plus)"}
"k">relations:
  "k">eps: "sym(grad(u))"
  "k">sigma0: "K*tr(eps)*I + 2*G*dev(eps)"
  "k">psi_plus: "1/2*ddot(sigma0, eps)"
  "k">sigma: "((1 - d)**2 + kres)*sigma0"
"k">equations:
  "k">u: "div(sigma) = 0"
  "k">d: "Gc*(d/l0 - l0*lap(d)) = 2*(1 - d)*H"
"k">initial_conditions:
  "k">H: "1.0e3*Gc/l0*Heaviside(0.5 - x)*exp(-((y - 0.5)/(0.5*l0))**2)"
"k">boundary_conditions:
  - "k">faces: [y-]
    "k">u: {dirichlet: [0, 0, 0]}
  - "k">faces: [y+]
    "k">u: {dirichlet: [0, "ubar", 0]}
  - "k">faces: all
    "k">u: {equation: "dot(sigma, n) = [0, 0, 0]"}
    "k">d: {neumann: 0}
"k">time:
  "k">end: 8.0
  "k">dt: 0.1
  "k">scheme: backward-euler
"k">solver:
  "k">newton: {line_search: true}
"k">output:
  "k">every: 10
  "k">fields: [u, d, H]
"k">monitors:
  "k">d_max: "max(d)"
  "k">cracked: "mean(Heaviside(d - 0.9))"
  "k">sigma_yy: "mean(sigma[1, 1])"

**Figure 1.** The model in the GUI: the Model Builder tree on the left, the settings of the selected node in the middle, the graphics window on the right.

**Figure 2.** Strain, undamaged stress, tensile energy and degraded stress.

**Figure 3.** Equilibrium and the damage equation.

**Figure 4.** The fields: displacement, damage, and the history H with its update rule.

3Results

Nodes: 41×41×5 = 8,405
Periodic: z
PDDO order: 2, horizon 2.015 Δx
Time scheme: backward-euler
Steps: 80 × Δt = 0.1
Newton iterations: 160
Factorizations: 159
Run time: 14 min

At first the plate responds elastically with the damage confined to the notch (6.3 % of the plate cracked). The mean stress peaks at 0.399 at t = 6.0; then the crack runs from the notch tip, the cracked fraction grows to 8.6 %, and the stress drops by 85 % to 0.060.

**Figure 5.** Damage at the final time, viewed along z: the crack has grown from the notch tip across the plate.

**Figure 6.** Vertical displacement on the deformed plate: the two halves separate along the crack.

Mean stress σyy against time: elastic rise, peak, and the drop as the crack grows. — **Figure 7.** Mean stress σ_yy against time: elastic rise, peak, and the drop as the crack grows.

**Figure 8.** Fraction of the plate with d > 0.9.

4Run It Yourself

The problem file ships with PDDO-MultiPhysics. Open it in the graphical interface and press Compute, or run it from the command line:

python -m pddo_mp gui examples/fracture_phase_field.yaml
python -m pddo_mp run examples/fracture_phase_field.yaml

Results go to results/fracture_phase_field/ next to the problem file: a ParaView time series (.pvd), the monitors as CSV, and the final state. Any entry can be changed on the command line, for example --set time.end=4 for a run to half the end time. The user manual describes every part of the problem file and of the interface.

1The Problem

2The Problem File

3Results

4Run It Yourself

Have a coupled problem of your own?