Physics, Units, and Formulation¶

This page summarizes the continuum model used by PhAST and maps the main symbols to the YAML and fluent Python interfaces. It is intended as the bridge between the paper formulation, the public examples, and the executable solver configuration.

For sparse-direct and autograd-enabled linear-solve details, see Sparse and Matrix-Free Solves.

Variational Fracture Model¶

Sharp Griffith fracture minimizes elastic energy plus crack-surface energy:

where \(\mathbf{u}\) is displacement, \(\Gamma\) is the crack set, \(G_c\) is the critical energy release rate, and \(\mathcal{H}^{n-1}\) measures crack surface. Phase-field fracture replaces the sharp crack by a scalar damage field \(d \in [0, 1]\), where \(d=0\) is intact and \(d=1\) is fully damaged. The damaged set can be interpreted as a regularized crack band whose width is controlled by \(\ell_0\); in practical meshes, the visible diffuse band is typically spread over a few elements around the crack centreline.

PhAST uses the regularized energy

The degradation function is

with a small residual stiffness \(\eta\) to avoid singular stiffness rows in fully damaged zones. The length scale \(\ell_0\) controls the width of the diffuse crack band.

The equivalent crack surface density is often written separately as

The displayed expression is the standard split internal energy. Complete potential energy also includes external work from prescribed tractions and body forces. The current linear AT1/AT2 damage operator retains conventional \(\eta\)-independent driving coefficients; the small \(\eta\) is therefore treated as a mechanics regularization rather than an exactly consistent term in that damage operator.

Governing Equations¶

For quasistatic brittle fracture, the displacement field satisfies mechanical equilibrium with degraded tensile energy:

with

The usual mechanical boundary conditions are

For AT2 with history-field irreversibility, the damage equation has the common linearized form

where

This is the default hard-history update. Optional smooth-history routes are separate differentiability approximations rather than the same discrete irreversibility mechanism.

The natural phase-field boundary condition is

For explicit dynamics, PhAST solves the momentum balance

with damage updated through the configured phase-field substep cadence.

AT1 and AT2 Crack-Density Terms¶

The regularization choice enters through \(w(d)\) and \(c_w\):

The AT1 threshold is often written as

This is the threshold implied by the same \(c_w=2/3\) normalization shown above: the implemented AT1 weak form has gradient coefficient \(3G_c\ell_0/4\) and source \(2\mathcal{H}-3G_c/(8\ell_0)\). The source coefficient is not an additional normalization convention.

In PhAST, choose the model with material.overrides.pf_model: AT1 or material.overrides.pf_model: AT2.

For AT1, the damage update is a bound-constrained variational problem. The lower bound must be enforced during the iterative solve; a post-solve clamp is not equivalent to the active-set solve.

Energy Splits¶

An unsplit degradation can allow compressive energy to drive damage and can remove compressive stiffness in damaged zones. PhAST therefore offers splits of the elastic energy into tensile and compressive parts:

and applies \(g(d)\) only to \(\psi^+\) for split formulations. The isotropic route is the intentional unsplit exception. The public configuration key is material.overrides.energy_split.

The distinction is visible in bending, impact, and branching examples, where compressive regions should not create artificial damage. The public dynamic gallery includes crack curvature and branching examples such as examples/dynamic/B2_kalthoff_winkler/ and examples/dynamic/B7_dynamic_crack_branching_comsol/.

For a picture-first explanation of these terms, see the visual glossary.

Mixed-mode impact	Dynamic branching

Plane Strain and Plane Stress¶

Two-dimensional examples use either plane strain or plane stress. The main difference is the effective elastic operator. In plane strain, \(\varepsilon_{zz}=0\) and the Lamé parameter is

In plane stress, \(\sigma_{zz}=0\) and the in-plane reduction uses

The public fracture examples choose the assumption through their material and workflow settings. Miehe-style spectral split studies are usually documented as plane strain, while thin PMMA dynamic examples are typically treated with plane-stress settings and an Amor-style split.

Staggered Solver Loop¶

The supported quasistatic fracture route uses alternate mechanics-damage iterations. Convergence ordinarily assesses both displacement and damage changes. Explicit dynamics performs one segregated pass per time step and solves damage at the configured cadence without a quasistatic inner stagger test. The joint monolithic \((\mathbf{u},d)\) route is experimental and does not use projected-CG active-set enforcement.

        flowchart TD
    A[Start load or time step] --> B[Freeze damage d]
    B --> C[Solve mechanics for displacement u]
    C --> D[Update tensile history field H]
    D --> E[Freeze displacement u]
    E --> F[Solve damage equation for d]
    F --> G[Apply irreversibility and bounds]
    G --> H{displacement and damage criteria satisfied?}
    H -- no --> B
    H -- yes --> I[Store fields, histories, and visuals]
    I --> J[Advance to next step]
    

One representative damage component of the convergence check is

where \(\varepsilon_{\mathrm{stag}}\) maps to the configured staggered tolerance.

Matrix-Free Operator View¶

PhAST’s public fracture path applies finite-element operators through gather-compute-scatter tensor kernels rather than assembling a global stiffness matrix for every operator application. The mechanics JVP used by the quasi-static Newton route is

and is evaluated in the implementation with an automatic-differentiation JVP. This is distinct from the legacy CG route, whose secant_matvec applies a perfectly linear operator formed from projections frozen at the current state. That frozen secant operator is not the consistent Jacobian for spectral splits because it omits eigenvector-rotation terms.

This is the matrix-free sense used throughout the repository: the solver materializes the action of the weak-form operator on a vector, while avoiding a persistent global stiffness matrix unless an optional sparse-direct pathway or preconditioner setup explicitly requires one.

Symbol to Configuration Map¶

Unit System¶

The public examples use a consistent mm-N-MPa-s convention unless a local README states otherwise:

Mesh resolution should be judged relative to \(\ell_0\). A common rule of thumb is

near the expected crack path, with local refinement around notches, holes, or branching regions.

This criterion is not a convergence guarantee. Refine \(h\) at fixed \(\ell_0\) when assessing discretization error, and study \(\ell_0\) separately because it changes the regularized fracture model.

References¶

• Bourdin, B., Francfort, G. A., and Marigo, J.-J. (2000). Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids.

• Amor, H., Marigo, J.-J., and Maurini, C. (2009). Regularized formulation of the variational brittle fracture with unilateral contact. International Journal for Numerical Methods in Engineering.

• Miehe, C., Welschinger, F., and Hofacker, M. (2010). Thermodynamically consistent phase-field models of fracture. International Journal for Numerical Methods in Engineering.

• Kumar, A., Francfort, G. A., and Lopez-Pamies, O. (2020). Revisiting nucleation in the phase-field approach to brittle fracture. Journal of the Mechanics and Physics of Solids.

• PhaseFieldX documentation: https://phasefieldx.readthedocs.io/