A MATLAB finite-element implementation of an atomistically-informed phase-field fracture model for graphene under uniaxial tension. The phase-field model is calibrated against molecular dynamics (MD) simulations performed with LAMMPS using the AIREBO interatomic potential.
The figure below compares the force–displacement response from the phase-field (PF) model and the MD simulation for a 50 × 50 nm pristine graphene sheet containing a central crack of 7.5 nm under uniaxial tension along the y-axis.
Both models predict:
- Linear elastic regime up to ~1.0 nm displacement with good agreement in stiffness
- Peak force of ~265–285 nN at ~1.0 nm displacement
- Brittle fracture: the MD simulation shows an instantaneous load drop; the phase-field model reproduces this sharp transition followed by a gradual post-peak tail arising from the AT2 regularisation
The slight stiffness overestimate in the PF model pre-peak is consistent with the continuum assumption not capturing atomic-scale bond non-linearity near the crack tip.
The model couples two sub-problems solved with a staggered Newton–Raphson scheme:
| Sub-problem | Governing equation |
|---|---|
| Elasticity | ∇·σ + b = 0, σ = g(φ)·C:ε |
| Phase field (AT2) | (G_c/l₀)φ − G_c l₀ Δφ = −g′(φ) H |
where H = maxτ∈[0,t] ψ⁺(ε(τ)) is the crack-driving history field and g(φ) = (1−φ)² is the AT2 degradation function.
Three strain-energy decompositions are implemented:
| Key | Model | Reference |
|---|---|---|
"NoSplit" |
Total energy (no split) | Default for mode-I graphene tension |
"Amor" |
Volumetric-deviatoric split | Amor et al. (2009) |
"Miehe" |
Spectral (principal-strain) split | Miehe et al. (2010) |
| Parameter | Value | Unit | Source |
|---|---|---|---|
| Young's modulus E | 1000 | GPa | Lee et al., Science 2008 |
| Poisson's ratio ν | 0.16 | — | DFT consensus |
| Fracture toughness K_Ic | 4.0 | MPa√m | Zhang et al., Nat. Commun. 2014 |
| Critical energy release rate G_c | 9.0 | J/m² (= nN/nm) | Calibrated from MD |
| Regularisation length l₀ | 0.2 | nm | Calibrated from MD |
| Graphene thickness h | 0.34 | nm | Interlayer spacing |
All quantities in the code use:
length → nm, force → nN, stress → nN/nm² = GPa, G_c → nN/nm = J/m²
The MD simulations used to calibrate and validate the phase-field model were performed following the methodology of Dewapriya (2022).
- Simulator: LAMMPS (Thompson et al., Comput. Phys. Commun. 2022)
- Potential: AIREBO (Stuart et al., J. Chem. Phys. 2000) with the REBO cut-off distance reduced to 2 Å to suppress non-physical strain hardening (Dilrukshi et al. 2015)
| Parameter | Value |
|---|---|
| Sheet dimensions | 50 nm × 50 nm |
| Central crack length (2a) | 7.5 nm |
| Crack type | Row of missing hexagons, horizontal |
| Loading direction | y-axis (armchair or zigzag, depending on orientation) |
The sheet dimensions satisfy the rule-of-thumb that length and width must exceed ten times the half crack length to avoid finite-size effects on K_I (Cleri et al. 1998).
| Step | Details |
|---|---|
| Equilibration | 25 ps, time step 0.5 fs |
| Ensemble | NPT (isothermal–isobaric), Nosé–Hoover thermostat |
| Temperature | 300 K |
| Boundary conditions | Periodic along x and y |
| Initial perturbation | ±0.001 nm atomic displacements along x, y, z |
| Loading | Uniaxial strain along y at strain rate 0.001 ps⁻¹ until fracture |
The virial stress tensor was used (Thompson et al., J. Chem. Phys. 2009):
σᵢⱼ = (1/V) Σ_α [ ½ Σ_β (R_i^β − R_i^α) F_j^αβ − mα vᵢ^α vⱼ^α ]
The graphene thickness was taken as h = 3.4 Å = 0.34 nm for volume calculations. Fracture was identified by a sudden drop in the far-field stress.
The MD force–displacement curve was used to calibrate the two phase-field fracture parameters G_c and l₀:
-
G_c was estimated from the mode-I stress intensity factor at crack initiation: K_I = σ_f √(πa), G_c = K_I² / E where σ_f is the far-field stress at the onset of crack propagation and a is the half crack length.
-
l₀ was chosen so that the phase-field crack bandwidth (≈ 4l₀) resolves the MD crack-tip process zone while remaining computationally tractable. The value l₀ = 0.2 nm corresponds to approximately one graphene unit cell width (~0.25 nm).
-
Both parameters were refined iteratively by matching the peak force and the displacement at peak from the MD curve.
graphene_phasefield/
│
├── main_graphene.m ← Entry point — run this
├── inputdata_graphene.m ← Simulation parameters & decomposition selector
├── mat_params_graphene.m ← Atomistically calibrated material constants
├── initiate_fun_graphene.m ← Variable and DOF initialisation
│
├── mesh/
│ ├── mesh_graphene.geo ← Gmsh geometry (half-model, nested refinement)
│ ├── mesh_graphene.msh ← Pre-generated Q4 mesh
│ └── gmsh_read_graphene.m ← Gmsh v2 ASCII mesh reader
│
├── src/
│ ├── fem/
│ │ ├── basis.m ← Shape functions & reference Jacobians
│ │ ├── basis_j.m ← Physical-space B-matrices (Bu, Be)
│ │ ├── def_grad.m ← Deformation gradient / Jacobian
│ │ ├── element_elastic.m ← Displacement element (small strain, Q4)
│ │ ├── elementphase_elastic.m← Phase-field element (isotropic AT2)
│ │ ├── gaussfun.m ← Gauss quadrature points and weights
│ │ ├── indexfun.m ← Global DOF index arrays
│ │ ├── phasefield_elastic.m ← Global assembly wrapper
│ │ └── shapefun.m ← Isoparametric shape functions
│ │
│ ├── mat_model/
│ │ ├── cons_model_elastic.m ← Plane-stress Hooke's law + energy splits
│ │ └── gfun.m ← Degradation function g(φ) and derivatives
│ │
│ ├── postproc/
│ │ ├── postproc_graphene.m ← Phase-field contour + F–u curve plots
│ │ └── store_data_graphene.m ← Result storage per load step
│ │
│ ├── solver/
│ │ ├── boundary.m ← Apply homogeneous Dirichlet BCs
│ │ ├── newton_raphson_disp_elastic.m ← NR solver, displacement sub-problem
│ │ ├── newton_raphson_phase_elastic.m ← NR solver, phase-field (active-set + line search)
│ │ ├── predofs.m ← Prescribed DOF treatment (stiffness only)
│ │ └── predofs1.m ← Prescribed DOF treatment (stiffness + RHS)
│ │
│ └── utilities/
│ ├── dev.m ← Deviatoric part of a tensor
│ ├── signfun_p.m ← Macaulay bracket ⟨x⟩₊ = (x+|x|)/2
│ └── tr.m ← Trace of a 3×3 tensor
│
├── docs/
│ └── figures/
│ └── phasefield_vs_MD.png ← PF model vs. MD simulation comparison
│
└── output/ ← Created automatically at runtime
├── phase_field_step*.png ← Crack propagation contours
├── force_disp_step*.png ← Force-displacement curves
└── force_displacement_final.txt
The mesh models the top half of a 50 × 50 nm graphene sheet (symmetry along the crack plane). A 7.5 nm central crack is embedded at the mid-height boundary.
# Requires Gmsh ≥ 4.x
gmsh mesh/mesh_graphene.geo -2 -o mesh/mesh_graphene.msh -format msh2Mesh sizing is controlled by nested Box fields in the .geo file:
| Zone | Element size | Strip half-width |
|---|---|---|
| Innermost (crack band) | l₀/3 ≈ 0.067 nm | l₀/2 = 0.1 nm |
| Zone 2 | 2l₀/3 ≈ 0.133 nm | l₀ = 0.2 nm |
| Zone 3 | 4l₀/3 ≈ 0.267 nm | 2l₀ = 0.4 nm |
| Zone 4 | 8l₀/3 ≈ 0.533 nm | 4l₀ = 0.8 nm |
| Background | 3.2 nm | — |
Open MATLAB, navigate to the repository root, and run:
main_grapheneKey parameters to adjust in inputdata_graphene.m:
input.solver.decom = "Amor"; % "NoSplit" | "Amor" | "Miehe"
input.load_inc = 1e-2; % nm per step
input.disp_max = 1.8; % nm — stop after this displacement| File | Description |
|---|---|
output/phase_field_step*.png |
φ contour at each saved step |
output/force_disp_step*.png |
Force–displacement curve |
output/force_displacement_final.txt |
Tabulated F–u data (nm, nN) |
output/checkpoint_step*.mat |
Full MATLAB workspace snapshots |
The displacement and phase-field sub-problems are solved alternately within each load step. Under-relaxation of φ is applied at the stagger level (pre-peak ω = 0.7, post-peak ω = 0.4) to suppress stagger oscillations independently of the inner NR.
Post-peak convergence of the naive clamped NR breaks down because the irreversibility constraint (φ ≥ φ_old) is applied after the linear solve, making K and R inconsistent. The solver in newton_raphson_phase_elastic.m fixes this with:
- Active-set irreversibility: constrained DOFs are identified before the solve and treated as Dirichlet BCs on K and R, keeping the system consistent at every iteration.
- Armijo line search: prevents overshooting past φ = 1 in the crack band.
- Post-peak fallback: if the strict tolerance is not met after 10 load-step halvings but displacement has converged, the best iterate is accepted (irreversibility is guaranteed by the active-set regardless).
The load increment is automatically halved (up to 6× pre-peak, 10× post-peak) when Newton fails to converge, and refined by a fixed factor at first crack onset (phi_max > 0.6).
- MATLAB R2021a or later (uses
parfor,sparse, string literals) - Parallel Computing Toolbox (for
parpool; disable by removingparpoolcall and replacingparforwithfor) - Gmsh 4.x (only needed to re-generate the mesh)
- Zakavati S. & Arash B. (2023). Atomistically-informed phase-field fracture modeling of defective graphenes. Oslo Metropolitan University.
- Dewapriya M.A.N. (2022). Modeling fracture of graphene: a molecular dynamics tutorial. MODULUS, 32(2), 10–12.
- Thompson A.P. et al. (2022). LAMMPS — a flexible simulation tool for particle-based materials modeling. Computer Physics Communications, 271, 108171.
- Stuart S.J., Tutein A.B. & Harrison J.A. (2000). A reactive potential for hydrocarbons with intermolecular interactions. Journal of Chemical Physics, 112, 6472–6486.
- Dilrukshi K.G.S., Dewapriya M.A.N. & Puswewala U.G.A. (2015). Size dependency and potential field influence on deriving mechanical properties of carbon nanotubes using molecular dynamics. Theoretical and Applied Mechanics Letters, 5, 167–172.
- Cleri F. et al. (1998). Atomistic simulations of materials fracture. Journal of the American Ceramic Society, 81, 501–516.
- Thompson A.P., Plimpton S.J. & Mattson W. (2009). General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions. Journal of Chemical Physics, 131, 154107.
- Lee C. et al. (2008). Measurement of the elastic properties and intrinsic strength of monolayer graphene. Science, 321, 385–388.
- Zhang P. et al. (2014). Fracture toughness of graphene. Nature Communications, 5, 3782.
- Miehe C., Hofacker M. & Welschinger F. (2010). A phase field model for rate-independent crack propagation. Computer Methods in Applied Mechanics and Engineering, 199, 2765–2778.
- Amor H., Marigo J.J. & Maurini C. (2009). Regularized formulation of the variational brittle fracture. Journal of the Mechanics and Physics of Solids, 57, 1209–1229.
- Bourdin B., Francfort G.A. & Marigo J.J. (2000). Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids, 48, 797–826.
MIT License — see LICENSE for details.
If you use this code in your research, please cite:
@software{graphene_phasefield,
author = {Zakavati, Shadab and Arash, Behrouz},
title = {Atomistically-Informed Phase-Field Fracture of Graphene},
year = {2025},
url = {https://github.com/YOUR_USERNAME/graphene_phasefield}
}