Brittle Fracture Simulation

• Brittle fracture simulation is the predictive computational modeling of crack initiation, growth, and coalescence in materials that fail without plastic deformation.
• It utilizes meshfree peridynamic methods combined with probabilistic collocation to capture material heterogeneity and energy release accurately.
• The framework demonstrates O(δ²) convergence and validates glass–ceramic fracture toughness through efficient damage evolution and uncertainty quantification.

Brittle fracture simulation refers to the predictive computational modeling of failure in materials that fracture without significant plastic deformation, typically under applied stresses that exceed a critical threshold. The central challenge lies in accurately capturing the nucleation, growth, and coalescence of cracks in materials whose response is governed by intrinsic atomic or microstructural mechanisms, geometric singularities, and stochastic material heterogeneity.

1. Fundamental Mathematical and Physical Frameworks

Brittle fracture mechanics is governed by the interplay of linear elasticity, energy release, and micromechanical failure criteria. The atomistic, continuum, and statistical-nonlocal regimes are bridged by tailored mathematical formalisms:

• Peridynamics: The Linear Peridynamic Solid (LPS) model generalizes classical elasticity by replacing partial derivatives with nonlocal integral operators over a horizon of radius δ. For a displacement field $u(x)$, the nonlocal dilatation $\theta(x)$ is employed, and the governing (state-based) peridynamic operator is

$$L_{H\delta} u(x) = - \int_{B_\delta(x)} (\lambda(x,y) - \mu(x,y)) K(|y-x|) (y-x)[\theta(x)+\theta(y)]\,dy -8 \int_{B_\delta(x)} \mu(x,y) K(|y-x|)\frac{(y-x)\otimes(y-x)}{|y-x|^2} [u(y)-u(x)]\,dy = f(x)$$

with $K$ a radially symmetric kernel, and the two-point averages of Lamé moduli $\lambda$, $\mu$ set by harmonic mean for heterogeneity (Fan et al., 2022). As $\delta\to 0$, $L_{H\delta}$ converges to the classical Navier operator for elasticity (Fan et al., 2022).

• Stochastic Representation: Micromechanical parameters $\lambda(x)$, $\mu(x)$ modeling material heterogeneity are represented as finite-dimensional random fields, such as through Karhunen–Loève or PCA expansions:

$$\lambda(x, \xi) = \lambda_0(x) + \sum_{n=1}^N \sqrt{\tau_n} \phi_n(x) \xi_n$$

with stochastic variables $\xi_n$ and $\rho(\xi)$ their probability measure (Fan et al., 2022).

• Fracture Nucleation and Growth: Damage evolution and crack tracking are handled locally through bond-based Boolean criteria: a bond $(i, j)$ is irreversibly broken if its stretch $s_{i,j}$ exceeds a critical value $s_0(x_i,x_j)$; for a Griffith-type fracture energy $G(x_i,x_j)$,

$$s_0(x_i,x_j) = \sqrt{ \frac{ G(x_i, x_j) }{ 4(\lambda_{i,j}-\mu_{i,j})\,\beta' + 8\mu_{i,j} \beta } }$$

with geometric constants $\beta,\beta'$ set by the horizon $\delta$ (Fan et al., 2022).

2. Numerical Methods: Meshfree Quadrature and Uncertainty Quantification

• Optimization-Based Meshfree Quadrature: To discretize the peridynamic integral operators on a quasi-uniform point cloud $\{ x_i \}$, an optimization-based quadrature is constructed at each center $x_i$ via the solution of

$$\min_{\{ \omega_{j,i} \}} \sum_j \omega_{j,i}^2 \quad \text{subject to} \quad \sum_j q_k(x_i, x_j) \omega_{j,i} = \int_{B_\delta(x_i)} q_k(x_i, y)\,dy$$

for all basis functions $q_k \in \mathcal{V}_{h,x_i}$, e.g., polynomials of degree ≤5 divided by $|y-x_i|^3$ (Fan et al., 2022). This guarantees asymptotic compatibility: as $h\to 0$ and $\delta \to 0$ (with $\delta/h$ fixed), the discrete peridynamic model converges to classical elasticity at $O(\delta^2)$ (Fan et al., 2022).

• Probabilistic Collocation Method (PCM): Uncertainty quantification is integrated via PCM by evaluating the deterministic peridynamic model at specially selected collocation points $\{\xi^k\}$ in parameter space. For $N$ stochastic dimensions, tensor-product Gauss or Chebyshev nodes (for $N\leq 4$) or Smolyak sparse grids (for $N>4$) are used. The stochastic solution is interpolated as

$$u_\delta^Q(x, \xi) = \sum_{k=1}^{Q} u_\delta(x, \xi^k) \ell_k(\xi)$$

with $Q$ deterministic samples and multivariate Lagrange polynomials $\ell_k$ (Fan et al., 2022). The expected value and variance of the solution follow as weighted sums over the collocation nodes.

• Convergence Guarantees: For analytic dependence of the solution $u_\delta(x,\xi)$ on $\xi$, algebraic

$$\|u_\delta(x, \xi) - u_\delta^Q(x, \xi)\|_{S_{H\delta}} \leq C_1 Q^{-\beta_1}$$

and, for sufficiently high sparse-grid levels, sub-exponential convergence rates for PCM can be achieved (Fan et al., 2022). Spatial convergence remains $O(\delta^2)$ under meshfree quadrature as $h\to 0$.

3. Computational Workflow: Damage, Load Incrementation, and Post-Processing

A practical brittle fracture simulation using the meshfree peridynamic-PCM framework involves:

1. Spatial Discretization: Generate a quasi-uniform point cloud $\{ x_i \}$ over the domain $\Omega \cup B_\delta(\Omega)$, selecting fill distance $h$ and horizon $\delta$ such that $\delta/h$ is constant (Fan et al., 2022).
2. Kernel and Modulus Precomputation: Evaluate $K_{i,j}=K(\|x_j-x_i\|)$, harmonic means $\lambda_{i,j}, \mu_{i,j}$, and optimization-based quadrature weights $\{ \omega_{j, i} \}$ at each quadrature center $x_i$.
3. Stochastic Sampling: Choose PCM collocation points $\{ \xi^k \}$ in the reduced stochastic space.
4. Per Sample Problem: For each $\xi^k$,
   • Initialize all bond variables $\gamma_{i,j} = 1$ (intact).
   • For each load increment:
     • Solve the linear peridynamic system for $u^n$ with current damage/tension.
     • Update nonlocal dilatations.
     • Evaluate bond stretches $s_{i,j}^n$; break new bonds by setting $\gamma_{i,j}^n = 0$ if $s > s_0$.
     • If new bonds break, repeat until quiescence under current load.
   • Record final displacement and bond-damage fields.
5. Statistical Post-Processing: Aggregate metrics (e.g., average fracture toughness) using collocation quadrature.

4. Benchmarking, Materials Applications, and Quantitative Validation

• Materials Instantiation: In glass–ceramic microstructures, phase fields $R(x, \omega) \in \{0,1\}$ distinguish glass vs. crystal. PCA reduces real microstructures to $N$-dimensional random field $\xi$, and spatially-varying elastic and fracture parameters are mapped to $E(x,\xi), G(x,\xi)$ (Fan et al., 2022).
• Crack Propagation and Toughness Estimation: Under simulated loading (e.g., top/bottom displacement), bonds are broken according to $s_{i,j} > s_0$; the resulting crack paths are classified by the medium traversed. The effective energy release rate is constructed from crack lengths in glass ($L_1$), crystal ($L_2$), and interfaces ($L_i$):

$$G_{IC} = \frac{G_1 L_1 + G_2 L_2 + (G_1 + G_2)L_i/2}{W}$$

with $W$ the projected crack length. The effective fracture toughness is $$K_{IC} = \sqrt{ \frac{ E_{\mathrm{eff}} }{1-\nu^2} G_{IC} }$$ where $E_{\mathrm{eff}} = (1-f)E_1 + fE_2$ for crystal volume fraction $f$ (Fan et al., 2022).

• Quantitative Agreement: In lithium-disilicate glass–ceramic simulations (800μm × 400μm samples, microstructural crystal volume fractions $f = \{20\%, 40\%, 60\%, 80\%\}$, $Q=41$ Smolyak collocation points), simulated average $K_{IC}(f)$ increases nearly linearly with $f$, matching experimental trends with substantially fewer samples than fully random Monte Carlo (Fan et al., 2022).

5. Damage Criteria, Crack Tracking, and Physical Realism

• Irreversible Bond Breaking: Each spatial bond $(i,j)$ possesses a state variable $\gamma_{i,j}(t)\in\{0,1\}$, broken irreversibly when the bond stretch exceeds threshold. Cracks emerge and propagate naturally as clusters of broken bonds; no explicit crack tracking or mesh adjustment is needed (Fan et al., 2022).
• Critical Stretch and Material Energy: The threshold $s_0(x_i,x_j)$ is explicitly tied to the local fracture energy via a derived relationship, ensuring physical consistency with Griffith-type energy dissipation (Fan et al., 2022).
• Load Incrementation: Simulations are performed under quasi-static incremental boundary displacement, ensuring the equilibrium of the damage field at each step. Nonlinear loops resolve within-step bond breaking, preventing instability or missed failure events.

6. Theoretical Properties, Scalability, and Limitations

• Convergence and Scalability: The meshfree quadrature ensures asymptotic compatibility; spatial discretization converges at $O(\delta^2)$ as $h,\delta\to 0$, with $\delta/h$ fixed. In stochastic space, PCM achieves algebraic or sub-exponential convergence (dependent on analytic regularity in parameters and Smolyak level) (Fan et al., 2022).
• Parallelizability: Each deterministic sample at a PCM collocation node is independent, enabling trivial parallel scaling for large parameter studies.
• Limitations and Applicability: The methodology is restricted to quasi-static fracture; dynamic effects (inertia, rate-dependent damage) are not included. While the framework can represent arbitrary parameter heterogeneity via random fields, computational complexity increases with stochastic dimension, though the use of sparse grids mitigates combinatorial explosion. Bond-based peridynamics, as employed, imposes a fixed Poisson ratio; state-based generalizations are required for arbitrary isotropy.

7. Context and Extensions

• Statistical Uncertainty Quantification: The PCM framework enables efficient quantification of uncertainty in fracture toughness and crack path statistics, matching or exceeding accuracy of Monte Carlo with fewer samples.
• Broad Applicability: While glass–ceramic microstructures exemplify the capability for real-world validation, the methodology is extensible to other brittle heterogeneous solids, provided appropriate parameter fields and fracture energies are supplied (Fan et al., 2022).

This approach integrates meshfree peridynamics, stochastic collocation, and optimization-based quadrature to provide a robust, scalable framework for the simulation and statistical analysis of brittle fracture in both homogeneous and heterogeneously toughened materials, with rigorous spatial and stochastic convergence guarantees and physical consistency connected to classical fracture mechanics.