Meso-scale size effects of material heterogeneities on crack propagation in brittle solids: Perspectives from phase-field simulations

TL;DR

This work investigates how mesoscale heterogeneities influence dynamic fracture in brittle solids using a variational phase-field approach. By modeling inclusions as an array with controlled size, spacing, and material contrasts, the authors quantify crack speed and fracture energy dissipation under Mode-I loading and identify a K-dominant zone, $D_K$, as a critical design scale. The key finding is that matching inclusion size to $D_K$ yields optimal toughening; for a fixed area fraction, many smaller inclusions near $D_K$ enhance fracture resistance, while for fixed inclusion size, larger, well-spaced inclusions approaching $D_K$ can maximize toughness gains. These insights enable a physics-guided pathway to design mesoscale heterogeneous materials and meta-materials that resist crack propagation in extreme environments.

Abstract

Brittle solids are often toughened by adding a second-phase material. This practice often results in composites with material heterogeneities on the meso scale: large compared to the scale of the process zone but small compared to that of the application. The specific configuration (both geometrical and mechanical) of this mesoscale heterogeneity is generally recognized as important in determining crack propagation and, subsequently, the (effective) toughness of the composite. Here, we systematically investigate how dynamic crack propagation is affected by mesoscale heterogeneities taking the form of an array of inclusions. Using a variational phase-field approach, we compute the apparent crack speed and fracture energy dissipation rate to compare crack propagation under Mode-I loading across different configurations of these inclusions. If fixing the volume fraction of inclusions, matching the inclusion size to the K-dominance zone size gives rise to the best toughening outcome. Conversely, if varying the volume fraction of inclusions, a lower volume fraction configuration can lead to a better toughening outcome if and only if the inclusion size approaches from above the size of the K-dominance zone. Since the size of the K-dominance zone can be estimated \textit{a priori} given an understanding of the application scenario and material availability, we can, in principle, exploit this estimation to design a material's mesoscale heterogeneity that optimally balances the tradeoff between strength and toughness. This paves the way for realizing functional (meta-)materials against crack propagation in extreme environments.

Figures (11)

• Figure 1: (a) The setup of our numerical model: a single-notched three-point bending beam subjected to a constant indentation speed $v_\text{load}$. The beam has a span $L$, a height $H$, and a notch with length $a$. (b) The geometric configuration of a single array of mono-sized ($d$) and equally spaced ($h$) inclusions embedded along a line of length $L_\text{in}$ that starts at a distance of $D_\text{in}$ (as a buffer zone) from the notch tip. The width $W = 10L_\text{in}$ considered for computing the area fraction is not drawn to scale. In this particular image, $N = 3$ using $c_0 = 0.2, N_0 = 5$. (c) Three examples of the geometrical configuration produced with $c_0 = 0.2, N_0 =5$ for $N = 5$, $N = 12$, and $N = 32$.
• Figure 2: (a) Evolution of normalized crack length $l/a$ as a function of normalized indentation displacement $v_\text{load}t/a$ for a homogeneous case with no inclusion (yellow curve) and a heterogeneous case with five inclusions (blue curve, $d = 0.2a$ resulting from setting $c_0 = 0.2, N_0 = 5$). The inclusion material has an elastic contrast $\alpha = 0.4$ and a toughness contrast $\beta = 2.4$. The two alternating shaded areas indicate two different types of crack tip locations for the heterogeneous case. The dark gray area indicates that the crack tip is outside of an inclusion, whereas the light gray area indicates that the crack tip reaches the boundary of or is penetrating through an inclusion. The (almost) constant crack length within each light gray area indicates that the crack is almost arrested as reaching the boundary of and penetrating through an inclusion. (b) A plot similar to (a) but showing the variation of the normalized instantaneous crack tip speed $V/v_R$, where $v_R$ is the Rayleigh wave speed of the base material. (c) A visualization of the final crack trajectory based on $\phi$ showing the interaction between the crack and the inclusions for the heterogeneous case. The white cross indicates the identified crack tip, and the white dashed rectangle indicates the region within which $[\langle V\rangle_l]_t$ and $[\langle G \rangle_l]_t$ are calculated. (d) A similar visualization as (c) but showing the result of the homogeneous case.
• Figure 3: Variations of the normalized stress $(\sigma_{xx}+\sigma_{yy})/E_0$ (blue squares) and $\phi$ (red circles) as a function of $r/a$ as the crack approaches (a), about to enter (b), penetrates through (c), nucleates at the interface (d), and merges and leaves (e) an inclusion of size $d = 0.2a$. The grey shaded area in each sub-figure shows the portion of that inclusion ahead of the crack tip (the black cross shown in all sub-figures). The time is set to zero at (a), and the value is relative, used to indicate the rate of crack propagation at each stage. In this particular case, crack entering the inclusion ($t_3-t_2$) takes more time than that for crack nucleating at the compliant-to-stiff interface ($t_4-t_3$) after entering the inclusion.
• Figure 4: (a) A schematic showing the variation of stress $\sigma$ as a function of the distance $r$ to the crack tip, based on either LEFM (blue curve) or the expected response of a brittle material (red curve). Three regions can be identified. Very close to the crack tip is the process zone (colored in grey) where both dissipation and nonlinear elasticity prevail chen2017instability, and LEFM breaks down. Following the process zone is an annulus called the K-dominant zone (colored in yellow) where LEFM holds. At the junction of these two zones, the crack nucleation stress $\sigma_c$ is achieved. After the K-dominant zone is where boundary effects play a role and where LEFM breaks down again. (b) A schematic resembling (a) but showing the expected outcome from a typical phase-field simulation (red curve). Immediately after the numerically identified crack tip ($r=0$) is the PF-regularized zone (colored in green and with a size of about $\ell$) representing mechanisms active in the process zone janssen2004fracture. It ends where the stress peaks at the numerical crack nucleation stress ($\sigma_c^\text{num}$) whose value is related to $\ell$ (see Eqn. \ref{['criticalstress']}) and $\delta$. Following the PF-regularized zone is the proposed "effective K-zone" (colored in yellow and with a size denoted as $D_K$) that ends where the PF simulation result deviates from the LEFM solution. The solid line-shaded area after the PF-regularized zone indicates stress deviations resulting from FE discretization errors and non-zero (but small) $\phi$ values. (c) A visualization showing the spatial distribution of $(\sigma_{xx}+\sigma_{yy})/E_0$ at the moment of crack initiation for the homogeneous case obtained using $\ell = \ell_0$. We indicate the crack tip location using the black cross and denote the vertical distance away from the crack tip as $r$.
• Figure 5: Variations of $(\sigma_{xx}+\sigma_{yy})/E_0$ (left axis) and $\phi$ (right axis) as a function of $r$ from the configuration shown in Fig. \ref{['kzonediscussion']}(c), obtained with fixed $\delta$. Symbols (hollow squares and circles) are from PF simulations, while solid lines are from LEFM calculations. For all simulations, we fix $\delta \simeq \ell_0/26$ near the notch tip area. The PF-regularized zone and the "effective K-zone" discussed in Fig. \ref{['kzonediscussion']}(b) are highlighted using the same color scheme. Results from the left to the right are from different $\ell$ used in simulations: $\ell = \ell_0/6$ for (a) with a zoomed-in plot shown in (d); $\ell = \ell/3$ for (b) with a zoomed-in plot shown in (e), and $\ell = \ell_0$ for (c) with a zoomed-in plot shown in (f).