<figure><img src="image_1.png" /></figure>

HAL
open science

# Crack nucleation and propagation in highly heterogeneous microstructure models based on X-ray CT images of real materials

Thanh Tung T.T. Nguyen, Julien Yvonnet, Qizhi Zhu, Michel Bornert, Camille Chateau

<figure><img src="image_2.png" /></figure>

## ► To cite this version:

Thanh Tung T.T. Nguyen, Julien Yvonnet, Qizhi Zhu, Michel Bornert, Camille Chateau. Crack nucleation and propagation in highly heterogeneous microstructure models based on X-ray CT images of real materials. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. 4p. hal-01514290

HAL Id: hal-01514290

https://hal.science/hal-01514290

Submitted on 26 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Crack nucleation and propagation in highly heterogeneous microstructure models based on X-ray CT images of real materials

TT. Nguyen <fnref n="1" />,<fnref n="2" />, J. Yvonnet <fnref n="1" />, Q.-Z. Zhu <fnref n="1" />, M. Bornert <fnref n="2" />, C. Chateau <fnref n="2" />

**Résumé** — In this work, crack nucleation and propagation in highly heterogeneous materials models, such as those obtained by micro-CT imagery of real materials, are investigated for the first time by means of the phase field method. The Phase field method is based on a variational formulation of brittle fracture with a regularized approximation of discontinuities. Its various advantages for voxel-based models of microstructures are discussed. A shifted strain split operator algorithm is proposed to handle unilateral contact within cracks in a very efficient manner.

**Mots clés** — Phase Field, Cracks, Nucleation, Microtomography, Voxel models, Heterogeneous materials.

## 1 Phase-Field method to simulate crack propagation

Starting from the pioneer works of Francfort and Marigo [1], numerical difficulties of the classical fracture framework can be overcome by a variational-based energy minimization framework for brittle fracture and by using a regularized description of the discontinuities related to the crack front. In this work, crack propagation in highly heterogeneous microstructures, such as segmented X-ray CT images of real material, which are used as direct input of the simulations, are investigated for the first time by means of the phase field method, which here follows the algorithmic framework proposed by Miehe et al. [2]. To handle the large voxel-based models, we propose a shifted algorithm to compute the strain tensor split used to handle unilateral contact which leads to a very simple and efficient algorithm. The advantages of such an approach are demonstrated for crack nucleation and propagation in microtomography models of cement-based materials.

### 1.1 Smeared approximation of cracks

Let $\Omega \subset \mathbb{R}^d$ an open domain describing a cracked solid. We introduce the time-dependent crack phase field $d \in [0, 1]$ characterizing for $d = 0$ the unbroken and for $d = 1$ the fully broken state of the material. A regularized crack phase field $d(x, t)$ is obtained from the boundary value problem [2] :

$$ \left\{ \begin{array}{ll} d - l^2 \Delta d = 0 & \text{in } \Omega, \\ d(\mathbf{x}) = 1 & \text{on } \Gamma, \\ \nabla d(\mathbf{x}) \cdot \mathbf{n} = 0 & \text{on } \partial \Omega, \end{array} \right. \qquad (1) $$

where $\Delta(.)$ is the Laplacian, $l$ is a regularization parameter describing the actual width of the smeared crack, and $\mathbf{n}$ the outward normal to $\partial\Omega$. It can be shown that (1) is the Euler-Lagrange equation associated with the variational problem :

$$ d(\mathbf{x}) = \operatorname{Arg} \left\{ \inf_{d \in S_d} \Gamma_l(d) \right\}, \qquad (2) $$

with $S_d = \{d|d(\mathbf{x}) = 1 \text{ on } \Gamma \forall \mathbf{x} \in \Gamma\}$ and $\Gamma_l(d) = \int_{\Omega} \gamma(d)d\Omega$, where $\Gamma_l(d)$ represents the total crack length, $\gamma(d)$ denotes the crack density function per unit area, defined by : $\gamma(d) = \frac{1}{2l}d^2 + \frac{l}{2}\nabla d \cdot \nabla d$ and which governs the regularization by the length scale parameter $l$ giving for $l \to 0$ the sharp crack topology.

<footnote n="1">¹ Université Paris-Est, Laboratoire Modélisation et Simulation Multi Echelle, MSME UMR 8208 CNRS, 5 Bd Descartes, 77454 Marne-la-Vallée Cedex 2, France {thanht-tung.nguyen, Julien.Yvonnet, qizhi.zhu}@univ-paris-est.fr</footnote>

<footnote n="2">² Université Paris-Est, Laboratoire Navier, CNRS UMR8205, ENPC, IFSTTAR, 6 et 8 avenue Blaise Pascal, 77455 Marne-la-Vallée Cedex, France., {michel.bornert, camille.chateau}@enpc.fr</footnote>
## 1.2 Brittle fracture modelling by using the phase-field method.

The variational approach to fracture mechanics provided by Francfort and Marigo [1] introduces the following energy functional for cracked body:

$$E(\mathbf{u}, \Gamma) = E_u(\mathbf{u}, \Gamma) + E_s(\Gamma) = \int_{\Omega \setminus \Gamma} W_u(\varepsilon(\mathbf{u})) d\Omega + g_c \mathcal{H}^{d-1}(\Gamma) \quad (3)$$

where $W_u$ is the elastic strain energy density function, $\varepsilon = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T)$, $\mathbf{u}$ is the displacement field, $g_c$ is the fracture toughness, and $\mathcal{H}^{d-1}$ is the Hausdorff surface measure giving the crack length ($d=2$) or surface ($d=3$). The term $E_u(\mathbf{u}, \Gamma)$ represents the stored elastic energy in the cracked body, and $E_s(\Gamma)$ is the energy required to create the crack according to the Griffith criterion. The a priori unknown crack surface is approximated by the crack density function. The work required to create a unit crack area can be rewritten as:

$$g_c \mathcal{H}^{d-1}(\Gamma) = g_c \int_{\Omega} \gamma(d) d\Omega. \quad (4)$$

Using an unilateral contact formulation [2] with an assumption of isotropic elastic behavior of the body accounting for damage induced by traction only, the elastic strain density function can be expressed as:

$$W_u(\mathbf{u}, d) = g(d)\Psi^+(\varepsilon(\mathbf{u})) + \Psi^-(\varepsilon(\mathbf{u})), \quad (5)$$

where $g(d)$ is a simple degradation function chosen as:

$$g(d) = (1 - d)^2 \quad (6)$$

and where $\Psi^+(\varepsilon) = \frac{\lambda}{2} (\langle Tr(\varepsilon) \rangle)_{\pm}^2 + \mu Tr\{(\varepsilon^\pm)^2\}$, with $\langle x \rangle_{\pm} = (x \pm |x|)/2$, and $\varepsilon^\pm$ and $\varepsilon^-$ are extensive and compressive parts of the strain tensor [2]. We introduce here [4] a shifted strain tensor split algorithms to avoid the nonlinearity of this decomposition:

$$\varepsilon_{n+1}^{+} \approx \tilde{P}_{n+1}^{+} : \varepsilon_{n+1}, \quad \varepsilon_{n+1}^{-} \approx \tilde{P}_{n+1}^{-} : \varepsilon_{n+1}, \quad (7)$$

where:

$$\mathbb{P}_{n+1}^{+} \approx \mathbb{P}_{n}^{+} + \beta(d_n) \Delta t_n \frac{\partial \mathbb{P}_{n}^{+}}{\partial t}. \quad (8)$$

where we have introduced a weight function $\beta(d)$ such that $\beta(d=0) = 0$, $\beta(d=1) = 0$, e.g. $\beta(d) = (1-d)^2$ to avoid unphysical values of $\varepsilon^\pm$ in the broken zone [4]. It yields the smeared form of the total potential energy for brittle fracture:

$$E(\mathbf{u}, d) = \int_{\Omega} \{g(d)\Psi^{+}(\varepsilon(\mathbf{u})) + \Psi^{-}(\varepsilon(\mathbf{u}))\} d\Omega + g_c \int_{\Omega} \gamma(d) d\Omega. \quad (9)$$

Minimization of $E(\mathbf{u}, d)$ and introduction of the history field $\mathcal{H}$ to enforce irreversibility [2]:

$$\mathcal{H}(t) = \max_{\tau \in [0, t]} \{\Psi^{+}(\mathbf{x}, \tau)\}, \quad (10)$$

and introducing a time-stepping, the problem is solved in a staggered procedure, where two linear problems are obtained at each iteration. At time $t^{n+1}$, given $\mathbf{u}^n(\mathbf{x})$ and $d^n(\mathbf{x})$, the weak form for the phase-field problem is given by:

$$\int_{\Omega} \left\{ \left( 2\mathcal{H}_n(\mathbf{u}^n) + \frac{g_c}{l} \right) d^{n+1} \delta d + g_c l \nabla d^{n+1} \cdot \nabla(\delta d) \right\} d\Omega = \int_{\Omega} 2\mathcal{H}_n(\mathbf{u}^n) \delta d d\Omega \quad (11)$$

and provides the phase field $d^{n+1}(\mathbf{x})$. The weak form for the elastic problem is obtained as:

$$\int_{\Omega} W_u(\mathbf{u}^{n+1}, d^{n+1}) : \varepsilon(\delta\mathbf{u}) d\Omega = \int_{\Omega} \mathbf{f} \cdot \delta\mathbf{u} d\Omega + \int_{\partial\Omega_F} \bar{\mathbf{F}} \cdot \delta\mathbf{u} d\Gamma \quad \forall \delta\mathbf{u} \in H_0^1(\Omega), \quad (12)$$

which provides the displacement field $\mathbf{u}^{n+1}(\mathbf{x})$.
## 2 Numerical example

A microtomography-based microstructure of a three-phases porous cementitious material is under consideration. The studied material is an EPS lightweight concrete [3], made from quartz sand and EPS beads embedded in a cement matrix. The grey level image was segmented in order to separate the three phases of the microstructure. The dimension of the sample is L = 1 mm (see figure 1) (b). As depicted in the same figure, the white, grey and black phases correspond to matrix, inclusions and pores, respectively.

<figure><img src="image_3.png" /><figcaption>FIGURE 1 – Compression test of a microtomography image-based model of cementitious material : (a) microstructure : the white, grey and black phases correspond to matrix, inclusions and pores, respectively ; (b) geometry of the domain and boundary conditions.</figcaption></figure>

On the lower end, the y-displacements are blocked while the x-displacements are free. On the upper end, the x-displacements are free, while the y-displacements are prescribed at value of U which increases during the simulation. Plane strain is assumed. The model consists of 550 × 550 pixels, each associated with a material property of matrix, inclusion or holes, according to the data obtained from the microtomography image segmentation. The voxel data are transferred into a regular grid of square domains associated with voxels, each divided into 2 3-nodes elements

The material parameters of each phase are : $E_i = 30$ Gpa, $v_i = 0.3$, and $E_m = 10$ Gpa, $v_m = 0.2$. The fracture toughness is $g_c = 2.5 \times 10^{-4}$ kN/mm. The simulation is performed with monotonic displacement increments of $\bar{U} = 10^{-4}$ mm during the first 110 time steps and $\bar{U} = 10^{-6}$ mm during the last 240 time steps which correspond to the propagation of the micro cracks. The length scale parameter is chosen as $l = 0.0075$ mm. In this example, the domain is initially not cracked, and the cracks first nucleate and then propagate with an increase of the compressive load.

The crack distribution evolution for different time steps is depicted in figure 2. We can observe that several cracks are nucleated from the pores and can propagate either in the matrix or in the inclusions, with complex paths. When the microcracks start nucleating, the materials strength quickly drops. Then, this example shows the potential of the method for describing microcracking, involving nucleation and complex crack patterns in real microstructures.

The load-displacement curve is provided in figure 3.

## Références

[1] G.A. Francfort and J.J. Marigo. *Revisiting brittle fracture as an energy minimization problem*, J. Mech. Phys. Solids, 46(8):1319-1342, 1998.

[2] C. Miehe, M. Hofacker, and F. Welschinger. *A phase field model for rateindependent crack propagation : Robust algorithmic implementation based on operator splits*, Comput. Meth. Appl. Mech. Eng., 199 :2776-2778, 2010.

[3] K. Miled, K. Sab, and R. Le Roy. *Particle size effect on eps lightweight concrete compressive strength : Experimental investigation and modelling*, Mech. Mater., 39 :222-240, 2007.
<figure><img src="image_4.png" /><figcaption>FIGURE 2 – Compression test of a microtomography image-based model of cementitious material : crack propagation for : (a) $\overline{U} = 20.5 \times 10^{-3}$ mm ; (b) $\overline{U} = 23.5 \times 10^{-3}$ mm ; (c) $\overline{U} = 24 \times 10^{-3}$ mm, and (d) $\overline{U} = 25 \times 10^{-3}$ mm.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 3 – Compression problem : load-deflection curve.</figcaption></figure>

[4] T.T. Nguyen, J. Yvonnet, M. Bornert, C. Chateau, Q.-Z. Zhu. *A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their microstructure*, submitted.