# A probabilistic model for fast-to-evaluate 2D crack path prediction in heterogeneous materials

Kathleen Pele, Jean Baccou, Loïc Daridon, Jacques Liandrat, Thibaut Le Gouic, Yann Monerie, Frédéric Péralès

► To cite this version:

Kathleen Pele, Jean Baccou, Loïc Daridon, Jacques Liandrat, Thibaut Le Gouic, et al. A probabilistic model for fast-to-evaluate 2D crack path prediction in heterogeneous materials. Journal of Theoretical, Computational and Applied Mechanics, 2023, Volume XV, pp.1-21. 10.46298/jtcam.8322. hal-03295926v4

HAL Id: hal-03295926

https://hal.science/hal-03295926v4

Submitted on 18 Jan 2023

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.

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

Distributed under a Creative Commons Attribution 4.0 International License
**Identifiers**

DOI 10.46298/jtcam.8322

OAI hal-03295926v4

**History**

Received Jul 30, 2021

Accepted Mar 23, 2022

Published Jan 12, 2023

**Associate Editor**

Anna PANDOLFI

**Reviewers**

J.M. CHANDRA KISHEN

Anonymous

**Open Review**

OAI hal-03739347

**Licence**

CC BY 4.0

©The Authors

# A probabilistic model for fast-to-evaluate 2D crack path prediction in heterogeneous materials

Kathleen PELE¹,⁴, Jean BACCOU³,⁴, Loïc DARIDON²,⁴, Jacques LIANDRAT¹,

Thibaut LE GOUIC¹, Yann MONERIE²,⁴, and Frédéric PÉRALÈS³,⁴

¹ Aix Marseille Université, CNRS, Centrale Marseille, I2M, UMR7353, Marseille, France

² LMGC, Université de Montpellier, CNRS, Montpellier, France

³ Institut de Radioprotection et de Sureté Nucléaire (IRSN), Cadarache, France

⁴ MIST, Université de Montpellier, CNRS, IRSN, France

This paper is devoted to the construction of a new fast-to-evaluate model for the prediction of 2D crack paths in concrete-like microstructures. The model generates piecewise linear cracks paths with segmentation points selected using a Markov chain model. The Markov chain kernel involves local indicators of mechanical interest and its parameters are learned from numerical full-field 2D simulations of cracking using a cohesive-volumetric finite element solver called XPER. This model does not include any mechanical elements. It is the database, derived from the XPER crack, that contains the mechanical information and optimizes the probabilistic model. The resulting model exhibits a drastic improvement of CPU time in comparison to simulations from XPER.

**Keywords:** Markov chain, concrete, cracking prediction, machine learning

## 1 Introduction

Aging of materials is a major issue in industrial applications. It is particularly the case in nuclear framework where the aging of cementitious materials can cause safety issues. Among various degradations due to aging, a special attention is devoted to concrete cracking induced by Internal Sulfate Attack (Safiuddin et al. 2018). These degradations generally lead to the development of a network of cracks. They considerably influence the strength of structures, reduce their tightness. Knowing the characteristics of cracking is essential in the study of the life extension of nuclear power plants.

At mesoscopic scale, concrete, used in design of nuclear power plant containments, can be viewed as a two-phase composite material with mortar matrix and aggregates inclusions. The granular particles are randomly distributed in a matrix of mortar. The heterogeneities of this kind of material (Gangnant 2015) and the different associated scales (Le 2011; Le et al. 2009) increase the complexity of its study.

In the context of research on the safety of Pressurized Water Reactors, IRSN, in collaboration with the LMGC through the joint MIST laboratory, has developed a micro-mechanical approach for the analysis of the behavior of materials during hypothetical accidental transients (Péralès 2005). This approach is based on the concepts of Cohesive Friction Zone Model (MZCF) (Monerie and Acary 2001) associated with numerical modeling methods for multi-body systems based on the Non Smooth Contact Dynamics (NSCD) approach (Jean et al. 2001). The developed parallelized numerical platform, XPER (eXtended cohesive zone models and PERiodic homogenization) (Péralès et al. 2010; Péralès et al. 2008), allows to simulate the initiation and two-dimensional propagation of multi-cracks in heterogeneous materials. However, each simulation is very expensive in terms of computing time, since it can involve several days of computation on a few dozen processors. This CPU time can therefore become prohibitive in the context of probabilistic safety studies.
In order to reduce this CPU time, this paper proposes the construction of a probabilistic model that allows to quickly predict a set of crack paths from the discretization of the microstructure associated with a local law of probability based on a Markov chain model. The transition core of the introduced model depends on two geometrical indicators. The parameters of the model are estimated from a set of training crack paths obtained numerically using full-field cohesive-volumetric finite element analysis with XPER. All the mechanics are contained exclusively in the training of the model.

Several recent works are dealing with the development of statistical and machine learning tools for cracking data analysis. They study for example crack classification or detection of type of cracks (Das et al. 2019; Kim et al. 2018). The crack path prediction has been also addressed in (Willot 2015; Le Mire et al. 2021) or (Bayar and Bilir 2019) but none of them leads to a surrogate model providing local information on the crack path which is the originality of our development.

This paper is organized as follows. Section 2.1 is devoted to the reminder of some elements related to cracking of materials that are exploited in the construction of the model and allows to define working hypotheses. Section 2.2 describes the test case considered in this document and the computer code used for the simulations. Sections 3 and 4 deal with the construction of the prediction model. After the introduction of the discretization of the problem, the two geometric indicators chosen to capture the local configuration of the aggregates are defined and an efficient procedure is proposed to evaluate them. On the basis of these indicators, a Markov chain model is then developed to perform the prediction of the crack path. Finally, the performance of the new model is studied in Section 5 for different shapes of aggregates.

# 2 Concrete crack

## 2.1 Phenomenology of cracking

The objective here is to describe the behaviour of cracks in concrete in order to retain general assumptions that characterize the crack path. The crack path is directly related to the heterogeneous composition of the concrete. In general, cracking studies are performed at the mesoscopic scale where concrete can be considered as a bi-material composed of a matrix (mortar) and inclusions (aggregates). In this paper, rectangular microstructures with 25 % uniformly distributed aggregates of different shapes are studied. Concrete is often considered as *quasi-brittle material* (Gangnant 2015). It is clearly established that aggregates have a strong influence (Chen et al. 2012) on the fracture crack path. This influence is due, on the one hand, to a *high fracture resistance of the aggregates* and, on the other hand, to the *relative weakness of the aggregate/matrix interface*. The interface properties depend on a zone surrounding the aggregates called the Interfacial Transition Zone (ITZ) (Elices and Rocco 2008; Pope and Jennings 1992). This area is very porous, reducing significantly its strength. The cracks thus preferentially propagate along the aggregate/matrix interfaces, see Figure 1 (Chen et al. 2012; Scholer 1966; Akçaoğlu et al. 2004).

<figure><img src="image_2.png" /><figcaption>**Figure 1** Example of a crack that propagates across the matrix and along the aggregate-matrix interface</figcaption></figure>

It is established that aggregates and their boundaries have a significant influence on the
mechanical behavior of concrete. Indeed, according to Husem (2003), the compressive strength of a concrete depends for about 20% on cement paste, for more than 10% on aggregates and about 70% on interphase. To summarize, the presence of aggregates leads to (Walker and Bloem 1960):

* heterogeneous local mechanical fields,

* weak zones (aggregate/matrix interfaces) along which the cracks propagate.

In addition, since the overall loading remains constant in orientation and just increases in intensity, the local stress state implies *a main direction of propagation of the crack*. The crack propagation is restricted to a half-space oriented ahead of the crack tip. In 2D, this half-space is an half-plane: the main crack can not move backwards but goes ahead of the crack tip where the material is still sound.

This analysis of the behavior of cracking in concrete leads to four working hypotheses:

(H1) the matrix is elastic brittle,

(H2) the aggregates are stiff and unbreakable,

(H3) the matrix aggregate interface is the weakest zone in the composite material,

(H4) the crack propagates mainly in mode I and thus in 2D in a half-plane oriented in the local direction of propagation, in front of the crack tip perpendicularly to the local tensile stress.

## 2.2 Numerical simulation of cracking

In this section, the numerical representations of the concrete and the XPER software, that allows to numerically simulate the crack, will be presented. Numerical samples of concrete are constructed from morphological descriptors. XPER implementation involves a multi-contact modeling strategy based on the Non-Smooth Contact Dynamics (NSCD) (Jean et al. 2001) method where cohesive models are introduced as mixed boundary conditions between each volumetric finite element. Each element or group of elements of the mesh can be considered as an independent body and the interface between bodies follows a frictional CZM (Nordmann et al. 2019) with no regularization nor penalization (Péralès et al. 2010). These descriptors are statistical and geometric information identified on real concrete, like n-point moments and spatial covariance (Jeulin 2015) allowing to characterize the spatial arrangements of different phases. The numerical microstructure are qualified as statistically similar to real concrete in the sense of these descriptors. An example of such microstructure is given in Figure 2.

<figure><img src="image_3.png" /><figcaption>**Figure 2** Example of microstructure with square inclusions randomly distributed in space and in orientation according to a hard core point process. [Left] Pre-notched test case submitted to uniaxial traction. [Right] Example of finite element mesh, zoom on the light gray zone on the bottom left corner of the left figure; volumetric finite elements are connected to each other using a cohesive zone model.</figcaption></figure>

The test case deals with a rectangular domain ($l \times L = 0.225 \text{ m} \times 0.6 \text{ m}$) submitted to an uniaxial traction, see Figure 2. The microstructure is composed of an isotropic matrix containing about 25% of polygonal aggregates (square shapes in Figure 2) randomly distributed in space and in orientation: the heterogeneous medium is thus considered as isotropic. Remark 2 provides comments on the choice of this volume fraction. Cracks can propagate through the matrix and not through aggregates, see Assumptions (H1) and (H2).

All full-field crack simulations are performed with XPER code. The CPU time for a complete microstructure cracking simulation like Figure 2 (from left to right) is about 44 hours on an Intel(R) Xeon(R) Gold 6126 CPU @ 2.60GHz connected in Infiniband - 256Go RAM - Nodes of 24
processors. The simulation parameters are listed in Table 1 to Table 3. For more information on the parameters, see (Blal et al. 2012; Bichet 2017; Bichet et al. 2016; Socié 2019; Socié et al. 2021).

<table><caption>**Table 1** Mechanical parameters of bulk materials in the matrix/aggregates composite simulations.</caption><thead><tr><th></th><th>Matrix</th><th>Aggregates</th></tr></thead><tbody><tr><td>Mass density ρ (kg/m³)</td><td>2900</td><td>2900</td></tr><tr><td>Young's modulus E (Pa)</td><td>12 × 10<sup>9</sup></td><td>60 × 10<sup>9</sup></td></tr><tr><td>Poisson's ratio ν</td><td>0.2</td><td>0.2</td></tr></tbody></table>

<table><caption>**Table 2** Mechanical parameters of cohesive models at the interfaces: C<sub>n</sub> (resp. C<sub>t</sub>) normal (resp. tangent) initial cohesive stiffness, σ<sub>0</sub> peak stress at onset of fracture, w toughness.</caption><thead><tr><th></th><th>Matrix-Matrix</th><th>Aggregate-Aggregate</th><th>Matrix-Aggregate</th></tr></thead><tbody><tr><td>C<sub>n</sub> = C<sub>t</sub> (Pa/m)</td><td>1 × 10<sup>17</sup></td><td>1 × 10<sup>17</sup></td><td>1 × 10<sup>17</sup></td></tr><tr><td>σ<sub>0</sub> (Pa)</td><td>4.6 × 10<sup>7</sup></td><td>2.4 × 10<sup>9</sup></td><td>1.4 × 10<sup>7</sup></td></tr><tr><td>w (J/m<sup>2</sup>)</td><td>20</td><td>∞</td><td>20</td></tr></tbody></table>

<table><caption>**Table 3** Numerical parameters of the simulations. See Figure 2. Right for an example of corresponding finite element mesh.</caption><tr><td>Area of aggregates</td><td>0.6 m × 0.225 m</td></tr><tr><td>Degrees of freedom</td><td>48 666</td></tr><tr><td>Number of processors</td><td>24</td></tr><tr><td>Computation time</td><td>44 h</td></tr></table>

Figure 3 shows an example of simulation with square-shape aggregates. The crack path corresponds to the colour discontinuity: the colour stands for vertical displacement, in red a material point goes up, in blue down, in-between the crack occurs. The main direction of crack path is therefore almost perpendicular to the direction of the load, referred as main mode I (perpendicular to the uniaxial tension).

<figure><img src="image_4.png" /><figcaption>**Figure 3** Crack path (colour discontinuity) obtained with XPER for the case described in Figure 2. Colours stand for the vertical displacement: in red, mater goes up; in blue, maters goes down; in-between crack occurs.</figcaption></figure>

It turns out that two numerical concretes sharing the same morphological descriptors (covariogram, percentage of aggregate) can exhibit different crack paths. To illustrate this point, we propose an example in Figure 4 with the covariogram (Jeulin 2015) under stationary and isotropic assumptions. The covariogram allows to study the spatial distribution and the relative

<figure><img src="image_5.png" /><figcaption>**Figure 4** Two different microstructures with same volume fraction of inclusion, same geometric covariogram C(h) and submitted to same load, while exhibiting two different crack paths.</figcaption></figure>

organisation of the phases of a random medium. It represents the probability that for a point x
located in a given phase A, here an aggregate, $x + h$ where $h$ is a translation vector is in the same phase:

$$C(x, x + h) = P(x \in A, x + h \in A). \qquad (1)$$

In the case of a stationary and isotropic process, the covariogram only depends on the norm $||h||$ and is noted for simplicity $C(||h||)$. As a result, crack path is a complex local process strongly influenced by the local configuration of the microstructure even when overall microscopic descriptors are well known. To different microstructures with same, at least, second moment descriptors can lead to different values for mechanical quantities of interest such as the tortuosity, see Figure 5. The latter is strongly correlated to permeability of cracked concrete which is of primary interest in the context of nuclear safety. The definition of tortuosity is recalled.

**Definition 1**
[Tortuosity] The tortuosity $\tau$ of a crack is defined in 2D as the ratio between the length $\ell$ of the crack path to the length $L_{\text{mean}}$ of the shortest path joining the first and last points of the crack: $\tau = \ell/L_{\text{mean}}$.

In view of the cracked microstructure of Figure 3, it appears that the *local* configuration of the aggregates strongly influences the path of the crack. This information is central in our construction of the new crack indicators.

<figure><img src="image_6.png" /><figcaption>**Figure 5** Tortuosity $\tau$ defined as the ratio of the crack length $\ell$ (—) over the length of the shortest path $L_{\text{mean}}$ (---) between two ending points at the scale of a few aggregates. The considered crack path is a piecewise linear approximation of the real crack: the edges of the aggregates are discretized and thus the crack path follows this discretization, the real crack path in the matrix is summarized as a straight path joining two aggregates. See also Section 3.1.</figcaption></figure>

# 3 Definition and evaluation of local crack indicators

In this section, several key notions are presented in order to satisfy Assumptions (H1) to (H4) describing the local behavior of a crack. The definition and the evaluation of the local crack indicators are described in Sections 3.2 and 3.3.

## 3.1 Discretization of the microstructure

Although the main crack propagates in the mode I direction, the local propagation of a crack depends strongly on the microstructure configuration (the aggregates) near the crack tip. Therefore, we propose a discretization of the microstructure, at this scale, allowing to estimate the step by step crack path. In this work, each aggregate is approximated by a polygon where each side is discretized by five points. This allows to simulate the propagation of the crack along the side, see Figure 6. The crack path can leave or join an aggregate at several positions along the side of the aggregate and not only at the corners of the polygon.

The matrix is not discretized which implies that either the crack crosses the matrix following a straight line between two points of the discretization belonging to two distinct aggregates, see Assumption (H1), or it follows the side of an aggregate, see Assumption (H3). For a microstructure E, the set of discretization points denoted by $Y_E = \{y_i^E\}_{i=1,...,N_E}$ constitutes the
<figure><img src="image_7.png" /><figcaption>**Figure 6** [Left] Real local microstructure and crack path. [Right] Discrete local representation of the microstructure and of the crack path, the blue dots stand for the discretization points ($y_i^E$) of the aggregate boundaries and the purple path represents a piecewise straight crack joining the points ($x_i^E$).</figcaption></figure>

discrete granular microstructure, see Figure 6. The crack path is then approximated from a subset of $\mathcal{Y}_E$ denoted by $\mathcal{X}_E = \{x_i^E\}_{i=1,...,m_E}$. We can define the discretized crack as follows.

**Definition 2** [Discretized crack] In a Cartesian coordinate system with the origin at the bottom left corner of the microstructure, the crack path can be parameterized by ($x^1, h_E(x^1)$) with

$$h_E(x^1) = \frac{x_{i+1}^{E,2} - x_i^{E,2}}{x_{i+1}^{E,1} - x_i^{E,1}} (x^1 - x_i^{E,1}) + x_i^{E,2}, \quad \forall x^1 \in [x_i^{E,1}, x_{i+1}^{E,1}] \quad (2)$$

where ($x_i^{E,1}, x_i^{E,2}$) denote the coordinates of $x_i^E$ in the Cartesian system.

Assumption (H4) implies that, starting from a point of the discretized crack path, only a subset of discretization points can be reached. They are located in a field of view defined below. An example is exhibited in Figure 7.

<figure><img src="image_8.png" /><figcaption>**Figure 7** [Left] Example of a field of view (blue area) ahead of the crack tip (crack in purple, crack tip is the white circle) according to the local direction of propagation. The crack can only propagates in the half space defined by the crack tip and the local direction of crack propagation. [Right] Screening effect, see Procedure 1: the three black aggregates screen all blue dots within shaded areas. This screening effect takes also place for further aggregates within the field of view.</figcaption></figure>

**Definition 3** [Field of view] If $E$ is a microstructure, $x \in \mathcal{X}_E$ the crack tip and $\vec{u}_x$ is the local direction of crack propagation at $x$, the field of view of the crack path at $x$ is defined by the area containing the set of $I$ discretization points $\{y_i^E\}_{i \in I} \subset \mathcal{Y}_E$ such that:

$$\langle \vec{x} y_i^E, \vec{u}_x \rangle \ge 0, \quad \forall i \in [1, \dots, N_E] \quad (3)$$

where $\langle \cdot, \cdot \rangle$ denotes the Euclidean scalar product in $\mathbb{R}^2$.

**Remark 1** It is possible to take into account a non constant local direction and locally orientate the field of view. In the numerical studies of this paper, $\vec{u}_x$ is a constant vector orthogonal to the loading
direction. We refer to (Pele 2021, ch.4) for an example where $\vec{u}_x$ is evaluated through a mechanical simulation. For simplicity, $\vec{u}_x$ corresponds here to the main mode I propagation.

Depending on the configuration associated of the tip of the crack, the field of view may not contain point in the same aggregate. This strongly influences the crack path and should be taken into account in the prediction model developped in this paper. Consequently, it is important to distinguish the two following configurations of the field of view:

**Configuration F1** field of view including possibilities to follow the aggregates sides or to cross the matrix, Figure 8 [Left].

**Configuration F2** field of view only including possibilities to cross the matrix, Figure 8 [Right].

<figure><img src="image_9.png" /><figcaption>**Figure 8** Two field of view configurations. White dots stand for the crack tip. Red dots stand for points that can be reached by the crack crossing the matrix. Yellow dots stand for points that can be reached by the crack and belonging to the crack tip aggregate. [Left] Configuration F1 where the crack can propagate along the same aggregate boundary. [Right] Configuration F2 where the crack should propagate within the matrix before reaching a new aggregate. The reachable points are subjected to the screening effect illustrated in Figure 7 [Right] and detailed in Section 3.3.1.</figcaption></figure>

## 3.2 Definition of the two indicators

This section is devoted to the definition of two geometrical indicators illustrated in Figure 9 and used to capture the local behavior of a crack. These indicators will characterize all the couples $(x, y_i^E)$ where $x$ is the tip of the crack and $y_i^E$ a point of the field of view $E$. The first indicator is the angle that each vector $\vec{x}y_i^E$ makes with the local direction of propagation, $\vec{u}_x$, see Assumption (H4). The second indicator is the Euclidean norm of $\vec{x}y_i^E$. These indicators are evaluated at each increment of the crack propagation.

**Definition 4** [Crack indicators] Let $E$ be a microstructure, $x \in \mathcal{X}_E$ be the tip of the crack and $\vec{u}_x$ be the local direction of propagation at $x$. If $\langle \cdot, \cdot \rangle$, resp. $||.||$, is the Euclidean scalar product, resp. norm, in $\mathbb{R}^2$, for any $y$ in the field of view of $x$, the two local indicators are defined by

$$d_x(y) = \|\vec{x}\vec{y}\|, \quad \theta_x(y) = \arccos\left(\frac{\|\langle\vec{x}\vec{y}, \vec{u}_x\rangle\|}{\|\vec{x}\vec{y}\| \cdot \|\vec{u}_x\|}\right). \qquad (4)$$

For any $x$, normalized distance and angle are also introduced as

$$\tilde{d}_x(y) = \frac{d_x(y) - d_{\min}}{d_{\max} - d_{\min}}, \quad \tilde{\theta}_x(y) = \frac{\theta_x(y) - \theta_{\min}}{\theta_{\max} - \theta_{\min}}, \qquad (5)$$

where $d_{\max}$, $d_{\min}$ and $\theta_{\max}$, $\theta_{\min}$ are the maximum and minimum values of the indicators associated to the candidate points in the field of view for a given location of the crack tip.

## 3.3 Evaluation of the indicators

This section is devoted to the evaluation of the two previous indicators. In Section 3.3.1, an efficient procedure to select the points actually reachable by the tip of the crack in the field of
<figure><img src="image_10.png" /><figcaption>**Figure 9** Indicators of the candidate point y for the crack tip x in configuration F2. The red points correspond to candidate points on neighbouring aggregates. In gray, distance $d_x(y)$ of candidate y to the crack tip. In blue, angle $\theta_x(y)$ between the potential crack direction and local direction of propagation.</figcaption></figure>

view is developed. In the second Section 3.3.2, a numerical study of these indicators is performed by analyzing a training set of microstructures cracked with XPER code.

### 3.3.1 Procedure to select the candidate points

For each point x of the crack, a straightforward strategy for indicator evaluation would require the computation of the distance and the angle, $d_x(y_i^E)$ and $\theta_x(y_i^E)$, for all the discretization points located in the field of view of x. However, following Assumption (H2), some points cannot be reached by the crack due to the presence of aggregates and should be deleted from the candidate points. For a given aggregate, the deletion procedure is based on the detection of an area not visible from x because it is hidden by the aggregate. This area is called in the sequel the shadow area. Procedure 1 identifies this shadow area and deletes the corresponding discretization points. To prevent a dependence to the local crack propagation orientation, the position and the shape of the aggregate, it relies on computational geometry (Boissonnat 2018).

**Procedure 1** [Deletion of the discretization points in the shadow area of an aggregate]

Let E be a discrete granular microstructure, $x \in \mathcal{X}_E$ be a crack tip and $\{y_j^E\}_{j=1,...,n_E}$ be the $n_E$ points discretizing the sides of an aggregate A in the field of view of x. For any $(i, k) \in \{1, ..., n_E\}^2$, we denote by $\angle(\vec{x}y_i^E, \vec{x}y_k^E)$ the angle between two vectors $\vec{x}y_i^E$ and $\vec{x}y_k^E$ in $[0, \pi]$. The deletion is based on the identification of the points lying in the shadow area of the aggregate. We first find the two points $y_1^*$ and $y_2^*$ such that $\angle(xy_1^*xy_2^*) = \max_{(i,k) \in \{1,...,n_E\}^2} \angle(xy_i^*xy_k^*)$. This angle defines a search cone as illustrated in Figure 10. A point z in the field of view belongs to the cone if:

$$ \det(\vec{xy}_2^*, \vec{x}z) \det(\vec{xy}_1^*, \vec{x}z) < 0. \qquad (6) $$

A point z in the cone belongs to the shadow area and is deleted if, see Figure 10:

$$ \det(\vec{xy}_1, \vec{y}_1^* z) \det(\vec{y}_1^* z, \vec{y}_1 \vec{y}_2^*) < 0. \qquad (7) $$

This procedure is repeated for all aggregates, allowing the capture of the candidate points of the field of view satisfying Assumption (H2).

<figure><img src="image_11.png" /><figcaption>**Figure 10** Illustration of Procedure 1: construction of the search cone (dark lines) associated to any aggregate A located within the field of view and screening effect (gray area).</figcaption></figure>
### 3.3.2 Analysis of the relevance of the indicators

In this section, we analyze the relevance of the two indicators angle and distance to capture the local characteristics of the crack path. This analysis is also exploited to define the different terms of the prediction model developed in Section 4.

The two indicators are evaluated for each discretization point located on the crack path using Procedure 1 for a set of 35 microstructures whose cracking is simulated by the XPER code following Section 2.2. In this study the local direction of the field of view is always orthogonal to the loading (mode I). The indicator evaluation is performed for the two configurations of the field of view: F1 in Figure 8 [Left] and F2 in Figure 8 [Right]. Figure 11 displays the results of the set of indicator pairs selected by the cracks of the 35 microstructures.

<figure><img src="image_12.png" /><figcaption>**Figure 11** Indicator values in configurations F1 [Left] and F2 [Right] associated to points on the crack for 35 microstructures with square aggregates. Each point represents a couple of not normalized values ($d_x, \theta_x$). The red dots stand for candidate points present on neighbouring aggregates (cross the matrix) and the green square dots stand for the candidate points belonging to the same aggregate.</figcaption></figure>

**Configuration F1, Figure 8 [Left]** Figure 8 [Left] confirms that the crack is likely to follow the aggregate side. This situation corresponds to the large number of green points with the shortest distance and a variation range between 0 and $\pi/2$ for the angle. For these points on the same aggregate, the distance indicator is not influential. It is interesting to observe that the green points are all located on a vertical line since the discretization of each aggregate is uniform. Therefore, when the crack is propagating along the aggregate side, the points reached by the cracks are equally spaced leading to the same value for the distance indicator.

It is important to note that the attraction of the point on the same aggregate may be reduced in some cases. Indeed, the observation of our training base highlights that the crack sometimes chooses to leave the aggregate side and to cross the matrix even if it has the possibility to follow the aggregate (red points in Figure 11). This change of behaviour in the crack propagation is due to two situations that simultaneously occur: following the same aggregate involves a large angle while crossing the matrix is associated to the combination of small angle and short distance.

**Configuration F2, Figure 8 [Right]** In this configuration, Figure 8 [Right], since there is no candidate point on the aggregate sides, the crack can only propagate into the matrix. The indicators values of the candidates points reached by the crack are in the bottom left hand corner of the graph (i.e short distance and small angle). This reflects the fact that the crack is at the same time constrained by the direction of the stress and attracted by the closest aggregates.

**Remark 2** The previous numerical tests show that the crack can reach an aggregate further away from the crack tip but located in the direction of propagation. This situation does not occur for high volume fraction of aggregates. A volume fraction of 25% was thus chosen for the analysis to consider all possible situations that could be encountered in mechanical applications. It is expected that the analysis and prediction model developed in Section 4 remain valid for a wide range of volume fraction including [40%, 70%] corresponding to concrete.
The two local indicators are therefore relevant to characterize and discriminate the local behavior of a crack. They are integrated into a probabilistic prediction model that is described in the next section.

# 4 Markov chain model for crack prediction

This section is devoted to the construction of the fast-to-evaluate surrogate model for crack path prediction. The crack is modelled by a Markov chain with a set of parameters. For given crack tip and direction of propagation, several points can be reached and one of them will be selected by the model as the next crack tip. The Markov chain model proposed in this section allows to associate a probability to each of these points. The future crack tip is then chosen by a random draw based on these probabilities. The key point of this type of model is the construction of a transition kernel that precisely defines the probability for the crack tip at position x to reach a new point y.

In this section, we fully describe the development of the model with a specific transition kernel involving the two indicators presented in Section 3.2. The parameters are estimated in Section 5.1. Finally, Section 4.2 fully describes the crack path prediction using the proposed model.

## 4.1 Construction of the model

The prediction model is based on a Markov chain model. This type of model relies on the property that the prediction of future states of a system only depends on its present state. It is therefore particularly well suited to the modeling of the crack whose propagation only depends on local information at the crack tip and not on the whole path before the tip. Let us recall the definition of a Markov chain. For more details on Markov chains, see (Norris 1997, ch.1).

**Definition 5** [A Markov chain] A sequence of random variables $(X_n)_{n \in \mathbb{N}}$ with values in a set $\mathcal{E}$ is a Markov chain of state space $\mathcal{E}$ if for all $k \in \mathbb{N}$, $\forall (x_0, \dots, x_{k+1}) \in \mathcal{E}^{k+2}$ such that $P(X_0 = x_0, \dots, X_k = x_k) > 0$,

$$P(X_{k+1} = x_{k+1} | X_0 = x_0, \dots, X_k = x_k) = P(X_{k+1} = x_{k+1} | X_k = x_k). \quad (8)$$

Here, $P(X_{k+1} = y | X_k = x)$ denotes the conditional probability of $X_{k+1} = y$ given $X_k = x$. We set $v_0(x_0) = P(X_0 = x_0)$, so that for all $(x_0, \dots, x_n) \in \mathcal{E}^{n+1}$,

$$P(X_0 = x_0, \dots, X_n = x_n) = v_0(x_0) \prod_{k=0}^{n-1} P(X_{k+1} = x_{k+1} | X_k = x_k). \quad (9)$$

Moreover, the chain is said to be homogeneous if for all $k \in \mathbb{N}$ and for all $(x, y) \in \mathcal{E}^2$, $P(X_{k+1} = y | X_k = x) = P(X_1 = y | X_0 = x)$.

The discretization points of the crack are assumed to be a realization of homogeneous Markov chain. The discussion of Section 3 leads to propose that the probability of transition is evaluated from a transition kernel depending on the two local geometrical indicators introduced in Definition 4. The following modelling defines the general structure of transition probability between two points in the crack.

**Modeling 1** [Crack transition probability] For any microstructure $E$, we suppose that $(X_i^E)_{i=1,\dots,m_E}$ is a sequence of random variables that constitutes a homogeneous Markov chain. $(x_i^E)_{i=1,\dots,m_E}$ defines a sequence of realizations of the Markov chain and the transition kernel of the chain defining the probability of transition from $x$ to $y$ is given by

$$P_E(X_{i+1}^E = y | X_i^E = x) = c_{x,\Lambda} f_{\Lambda}(\tilde{d}_x(y), \tilde{\theta}_x(y)), \quad (10)$$

where $P_E$ denotes the conditional probability given $E$, $c_{x,\Lambda}$ is a normalization constant in order to define a probability and $f_{\Lambda}$ is a function parameterized by a set $\Lambda$ of parameters to be determined.
We set the first realization $X_0 = x_0$: it is the starting point of the crack and is chosen in the microstructure. In particular we have,

$$P_E(X_0^E = x_0, \dots, X_{m_E}^E = x_{m_E}^E) = \prod_{i=0}^{m_E-1} P_E(X_{i+1}^E = x_{i+1}^E | X_i^E = x_i^E). \quad (11)$$

The aim of this model is to reproduce the crack behaviour by weighting the probability of propagation of the crack with the importance of the indicators according to the local configuration at the tip of the crack in the microstructure. Since the presence of candidate points on the aggregate sides influences the local direction of propagation, the expression of $f_\Lambda$ depends on the configuration (F1 or F2). We use a decreasing exponential function in order to penalize high values of the indicators. Based on Modeling 1 and our assumptions, the following model presents the final transition probability by specifying $c_{x,\Lambda}$ and $f_\Lambda$.

**Modeling 2**
[Probabilistic cracking model]

Let $x$ be a discretization point reached by the crack and belonging to an aggregate A. If $\{y_k^E\}_{k=1,...,K_E}$ is the set of candidate points in the field of view of $x$, then the transition kernel is given by:

**Configuration F1, Figure 8 [Left]**

Let $\Lambda_{F1} = (\mu_1, \mu_2, \mu_3, \mu_4, \mu_5, \mu_6)$ be the parameters associated to this configuration. If $\{y_k^E\}_{k=r_E+1,...,K_E}$ defines the set of candidate points on A and $\{y_k^E\}_{k=r_E+1,...,K_E}$ the candidate points on the other aggregates, then $c_{x,\Lambda} = c_{x,\Lambda_{F1}}$ and $f_\Lambda = f_{\Lambda_{F1}}(\tilde{d}_x(y), \tilde{\theta}_x(y))$ with

$$f_{\Lambda_{F1}}(\tilde{d}_x(y), \tilde{\theta}_x(y)) = \begin{cases} e^{-\mu_1(\tilde{\theta}_x(y))^{\mu_2}} & \text{if } y \in A \\ e^{-\mu_3(\tilde{d}_x(y)\tilde{\theta}_x(y))^{\mu_6} - \mu_4(\tilde{d}_x(y))^{\mu_5}} & \text{else} \end{cases} \quad (12)$$

and the normalizing constant is

$$c_{x,\Lambda_{F1}} = \frac{1}{\sum_{k=1}^{K_E} f_{\Lambda_{F1}}(\tilde{d}_x(y), \tilde{\theta}_x(y))} \quad (13)$$

**Configuration F2, Figure 8 [Right]**

Let $\Lambda_{F2} = (\lambda_1, \lambda_2, \lambda_3, \lambda_4, \lambda_5, \lambda_6)$ be the parameters associated to this configuration, then $c_{x,\Lambda} = c_{x,\Lambda_{F2}}$ and $f_\Lambda = f_{\Lambda_{F2}}(\tilde{d}_x(y), \tilde{\theta}_x(y))$ with

$$f_{\Lambda_{F2}}(\tilde{d}_x(y), \tilde{\theta}_x(y)) = e^{-\lambda_1(\tilde{d}_x(y)\tilde{\theta}_x(y))^{\lambda_5} - \lambda_2(\tilde{d}_x(y))^{\lambda_6} - \lambda_3(\tilde{\theta}_x(y))^{\lambda_4}} \quad (14)$$

and the normalizing constant is

$$c_{x,\Lambda_{F2}} = \frac{1}{\sum_{k=1}^{K_E} f_{\Lambda_{F2}}(\tilde{d}_x(y), \tilde{\theta}_x(y))} \quad (15)$$

The different variables integrated in this modeling come from the analysis performed in Section 3.3.2. The mixed term allows to take values of distance and angle of the same order of magnitude. The distance term allows to take into account the neighboring attraction of the aggregates. The angle term allows to take into account the direction of propagation. In configuration $F_1$, the distinction between the possibility of following the aggregate and crossing the matrix is taken into account in the expression of $f_{\Lambda_{F1}}$ that depends on the two previous situations. Moreover, as observed in Section 3.3.2 for candidate points on the same aggregate, the distance has no influence and therefore, only the angle indicator is included in the model. For candidate points that require crossing the matrix in configuration $F_1$, experimentally, it is the presence of an aggregate close to the crack tip that has the most important influence on the change of crack behavior. That is why the angle indicator is only taken into account in the interaction term contrarily to the transition kernel in configuration $F_2$. It allows reducing the number of parameters in the model and the variability of their estimate.

More generally, in both configurations $F_1$ and $F_2$, the interaction term between distance and angle allows integrating the effect of the combination of small values of the two indicators that was observed in Section 3.3.2. It constrains the crack to follow a propagation direction while keeping the attraction of the closest aggregates
The model parameters are estimated from a training set of microstructures whose cracking has been simulated with XPER code. Under the independence assumption of the microstructures, the estimate is performed by maximization of the likelihood. More precisely, if {$E_i$}$_{i=1,...,Q}$ denotes the training set of microstructures and $(x_j^E)_{j \in I_{E_i}^{F_1}}$ and $(x_j^E)_{j \in I_{E_i}^{F_2}}$ are the sequences of I points selected by the crack on the microstructure $E_i$ in both configurations, $\Lambda_{F_1}^*$ and $\Lambda_{F_2}^*$ satisfy:

* Configuration F1:

$$ \Lambda_{F1}^* = \arg\max_{(\mu_1, ..., \mu_6)} \left\{ \prod_{i=1}^{Q} \prod_{j \in I_{E_i}^{F_1}} P_{E_i}(X_{j+1}^{E_i} = x_{j+1}^{E_i} | X_j^{E_i} = x_j^{E_i}) \right\} \quad (16) $$

* Configuration F2:

$$ \Lambda_{F2}^* = \arg\max_{(\lambda_1, ..., \lambda_6)} \left\{ \prod_{i=1}^{Q} \prod_{j \in I_{E_i}^{F_2}} P_{E_i}(X_{j+1}^{E_i} = x_{j+1}^{E_i} | X_j^{E_i} = x_j^{E_i}) \right\} \quad (17) $$

where $P_{E_i}$ is given by Modeling 2.

## 4.2 Cracking prediction with the probabilistic model

This section describes the crack path prediction for a given microstructure using the Markov chain model introduced in Modeling 2. The parameters of the transition kernel are assumed to be known, their estimate is studied in Section 5. This prediction relies on a procedure that starting from the crack tip, provides the next point of the field of view that is reached by the crack.

More precisely, the local indicators introduced in Definition 4 are first evaluated integrating Procedure 1 to reduce the set of candidate points in the field of view. Then, the Markov chain model can evaluate the probability of each point of the set to be the next point reached by the crack. Finally, a realization is retained according to the evaluated probability to select the next point of the crack. This procedure is described below:

**Procedure 2** [Local prediction of the next point of the crack] Given $x_i^E$ the crack tip in a microstructure E and $\overrightarrow{u_{x_i^E}}$ the local propagation direction, the prediction includes three steps:

**Step 1** Construct the field of view of $x_i^E$, identify the candidate points following Procedure 1 ($\{y_k^E\}_{k=1,...,K_E}$ is the set of remaining points) and evaluate the normalized indicators $\tilde{d}_{x_i^E}(y_k^E)$ and $\tilde{\theta}_{x_i^E}(y_k^E)$ $\forall k$.

**Step 2** According to the configuration, compute $\{P_E(X_{i+1}^E = y_k^E | X_i^E = x_i^E)\}_{k=1,...,K_E}$ using Modeling 2.

**Step 3** Select a realization from the discrete law $\sum_{k=1}^{K_E} P_E(X_{i+1}^E = y_k^E | X_i^E = x_i^E) \delta_{y_k^E}$.

In this procedure, the objective is to randomly drawn the next point of the crack among the points with the highest probabilities and not to select the candidate point with the highest probability. It is considered that several points can be good candidates for to be the next point of crack. Starting from the crack tip (initial position), this procedure is successively applied to any new point on the crack until the boundary of the domain is reached.

Thus the prediction model is stochastic. Starting from the same initial position, it can be applied to obtain several realizations of the crack path. Therefore, it allows to quantify the uncertainty associated to the prediction and to the mechanical quantities of interest. In Section 5, tools to use this set of realizations to determine the most optimal crack path will be presented. Algorithm 1 summarizes the three steps of the prediction procedure and Figure 12 provides an illustration of the procedure for configuration F2.

## 5 Numerical application

Several types of numerical applications are considered in this section. They first concern the estimate of the parameters (Section 5.1) of the transition kernel. Then, we focus on the ability of the model to correctly predict the crack path (Section 5.2) and the associated mechanical quantities of interest.
1 E is the discrete granular microstructure

2 $x^E$ is the point of crack tip

3 DirectionProp is the local direction in mode I

4 MarkovChainModel is Modeling 2 (gives the probabilities for the candidate points)

5 RandWeighted performs a weighted draw for select the next point of the crack

6 ParamModel are the estimated parameters of the prediction model

7 function CRACKING(E, x₀^E, ParamModel, DirectionProp)

8 FieldviewAgg ← set of aggregates in the field of view of point $x^E$ associated to DirectionProp

9 FieldviewPoint ← set of points located on the aggregates of FieldviewGran

10 for G in FieldviewAgg do

11 ShadowG ← shadow zone of aggregate G

12 for point in FieldviewPoint do

13 if point in ShadowG then

14 Delete point in FieldviewPoint

15 end if

16 end for

17 end for

18 ProbaPoint ← MarkovChainModel(FieldviewPoint,ParamModel)

19 PositionPointChosen ← RandWeighted(ProbaPoint)

20 return PositionPointChosen

21 end function

<figure><img src="image_13.png" /><figcaption>**Figure 12** Prediction step in Procedure 1 for configuration F2 with a zoom on the region of interest. The field of view is orientated with respect to a propagation direction orthogonal to the stress effort. The orange point is the crack tip, the red points are candidate points in the field of view where the probability to be reached is evaluated, and the green points are the candidate points selected by the model.</figcaption></figure>

## 5.1 Parameter estimate

The parameter estimate exploits a training set of 35 microstructures numerically cracked with the XPER code from the test case described in Section 2.2 with square shape aggregates. The estimation of the parameters of the Markov chain are obtained by maximising the likelihood, see Equations (16) and (17). The maximisation is performed with an optimized version of the BFGS algorithm (Xiao et al. 2008). Table 4 provides the estimated values for $\Lambda_{F1}$ and $\Lambda_{F2}$. Let us recall

<table><caption>**Table 4** Estimated parameters for $\Lambda_{F1}$ and $\Lambda_{F2}$, $\Lambda_{F1} = (\mu_1, \dots, \mu_6)$ and $\Lambda_{F2} = (\lambda_1, \dots, \lambda_6)$</caption><thead><tr><th>μ<sub>1</sub></th><th>μ<sub>2</sub></th><th>μ<sub>3</sub></th><th>μ<sub>4</sub></th><th>μ<sub>5</sub></th><th>μ<sub>6</sub></th><th>λ<sub>1</sub></th><th>λ<sub>2</sub></th><th>λ<sub>3</sub></th><th>λ<sub>4</sub></th><th>λ<sub>5</sub></th><th>λ<sub>6</sub></th></tr></thead><tbody><tr><td>7.06</td><td>4.1</td><td>30.2</td><td>8.9</td><td>0.2</td><td>0.85</td><td>34.2</td><td>9.2</td><td>13.16</td><td>1.79</td><td>1.08</td><td>0.42</td></tr></tbody></table>

that in Modeling 2, the parameters $\mu_1, \mu_3, \mu_4, \lambda_1, \lambda_2$, and $\lambda_3$ correspond to multiplicative factors, the others are exponents.

In configuration F1, the two highest values of the multiplicative factors ($\mu_3$ and $\mu_4$) and the two lowest values of the exponents ($\mu_6$ and $\mu_5$) correspond respectively to the parameters of the distance and the interaction terms. These two quantities are associated to candidate point located
on another aggregate. As a result, the crack will tend to favour a path on the same aggregate unless the angle to stay on the aggregate is high and the crossing of the matrix is associated to low distances and angles. This is consistent with the analysis of Section 3.3.2.

It is also interesting to study the variability of the estimate with respect to the size of the training set. More precisely, the parameter estimate is performed from a training set including an increasing number of microstructures. Figures 13 and 14 display the estimate of each parameter.

<figure><img src="image_14.png" /><figcaption>**Figure 13** Parameter estimate with respect to the size of the training set in configuration F1. From top to bottom, parameters associated to $\tilde{\theta}_x$ (points on the same aggregate), to the interaction term (points requiring crossing the matrix) and to $\tilde{d}_x$ (points requiring crossing the matrix)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14** Parameter estimate with respect to the size of the training set in configuration F2. From top to bottom, parameters associated to $\tilde{d}_x$, to the interaction term and to $\tilde{\theta}_x$.</figcaption></figure>

It is possible to observe a first stabilization for the majority of the parameters from about ten microstructures. In all cases, beyond 25 microstructures, the estimate can be considered to be stabilized for all parameters.
## 5.2 Cracking prediction

We evaluate the performances of the prediction model with the estimated parameters provided by Table 4. Since the model is stochastic, for a given microstructure, its output is a random crack path. When run M times, we obtain a set of M cracks.

In this section, the model is applied to the training set for verification and to a test set (i.e. not used for the parameter estimate) of 30 microstructures with multiform aggregates (aggregates represented by regular polygons with a number of sides varying from 3 to 8) for validation. It is important to keep in mind that the change of aggregate shapes between the training and the test sets does not require new XPER simulations to construct the prediction model, and therefore there is no extra significant computational time for the crack prediction.

We first focus on the analysis of the crack paths then consider the prediction of the tortuosity.

### 5.2.1 Crack path

For a given microstructure E, we denote by $\mathcal{X}_E^k = (x_{i,k}^E)_{i=1,...,m_E^k}$ the sequence associated to the $k^{th}$ ($k = 1, ..., M$) crack realization.

**Median path** When the objective is to derive a unique crack path that approximates the deterministic computer code simulation at a low cost, we use a median path from the set of points $\mathcal{X}_E^*$ satisfying

$$ \mathcal{X}_E^* = \arg\min_{\{\mathcal{X}_E^k\}_{k=1,...,M}} (\sum_{j=1, j \neq k}^M D(\mathcal{X}_E^k, \mathcal{X}_E^j)) \quad (18) $$

where D is the Fréchet distance (Driemel et al. 2012; Jekel et al. 2018) defined as follows.

**Definition 6** Let $X_1, X_2: [0, 1] \to \mathbb{R}^2$ be parametrizations of two crack paths. The Fréchet distance between $X_1$ and $X_2$ is

$$ D(X_1, X_2) = \inf_{\alpha, \beta \in [0, 1]} \max_{t \in [0, 1]} \|X_1(\alpha(t)) - X_2(\beta(t))\| \quad (19) $$

where $\alpha, \beta: [0, 1] \to [0, 1]$.

The median crack is considered as the most optimal crack for a given microstructure.

**Confidence region** To quantify the uncertainty of the median, we define a confidence region. The confidence region is constructed from point-value evaluations of percentiles. This type of construction requires a parameterization of each predicted path, see Definition 2. In this test case, the crack is initialized on the left hand side boundary of the microstructure and propagates until the failure of the microstructure. For any $x^1 \in [0, L]$, the uncertainty is quantified by estimating the 5% and 95% percentiles from the sample $(\tilde{h}_E^1(x^1), ..., \tilde{h}_E^M(x^1))$. Denoting $(\tilde{h}_E^1(x^1), ..., \tilde{h}_E^M(x^1))$ the ordered sample (increasing order) and using the notation $[.]$ for ceiling and $[.]$ for floor, the uncertainty interval at $x^1$ is defined by $[\tilde{h}_E^{[0.05M]}(x^1), \tilde{h}_E^{[0.95M]}(x^1)]$ which constitutes the upper curve and the lower curve. The region of confidence is finally taken as the convex hull and the uncertainty of the model prediction can be quantified by the diameter of this region computed as the Fréchet distance between the lower and upper curves defining this hull. The confidence region gives the area where the crack has the highest probability of passing.

Figures 15 and 16 provide two examples of median crack path and confidence region as well as a comparison with the XPER simulation for a microstructure of the training set (only square inclusions) and a microstructure of the test set (various shape inclusions). The number of paths calculated by the prediction model is $M = 100$. To complement the previous comparison, Figures 17 and 18 display the Fréchet-distance-based error between the XPER simulation and the median path as well as the uncertainty for all the microstructures of the training and test sets. For a better understanding of their values, the results of Figures 15 and 16 correspond to the microstructures 31 and 28 in these figures. It can be observed that the error between the median crack path and the XPER simulation is always smaller than the uncertainty. The accuracy of the median prediction can be reduced for some microstructures of the training and test sets. As an example, Figure 19 shows the case associated to the largest error on the test set.
<figure><img src="image_16.png" /><figcaption>**Figure 15** Example of path prediction for square aggregates: XPER simulation (dashed purple), median path (solid green), region of confidence (dotted red).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16** Example of path prediction for multiform aggregates: XPER simulation (dashed purple), median path (solid green), region of confidence (dotted red).</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 18** Prediction on the multiform aggregate test set: Fréchet-distance-based error (green square dots) between the median crack and XPER simulations; largest diameter of the confidence region (red dots).</figcaption></figure>

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

This lack of accuracy can be explained by a situation encountered by the crack where two candidate points of the field of view have similar values of their indicators. As a result, they exhibit the same probability to be the next point on the crack path. Since our model is stochastic, these two points are therefore reached by several crack realizations and it is not surprising that, for some microstructures, the median crack path does not correspond to the same local choice
than the XPER simulation. However, despite this local bifurcation, the rest of the crack path is in good agreement with the mode I propagation direction. Moreover, the XPER simulation lies in the region of confidence. Note that this result holds for both microstructures of the training and test sets.

<figure><img src="image_20.png" /><figcaption>**Figure 19** Example of path prediction associated to the largest error on the test set for multiform aggregates: XPER simulation (dashed purple), median path (solid green), region of confidence (dotted red).</figcaption></figure>

**Influence of the training set** We also study the influence of training set on the predictions for microstructures of the test set. Besides the training set with square aggregates, we consider two new training sets of 35 microstructures with pentagonal, resp. octagonal, aggregates numerically cracked with XPER code under the same conditions as those of Section 2.2. The parameters of the prediction model are estimated with the BFGS algorithm. Figure 20 gives the results associated to the construction of the median crack path for the three training sets. The set with square aggregates leads to the smallest error for the largest number of microstructures. However on the 30 microstructures studied, the three training bases show good performances on the test set with the same error magnitude. We can therefore conclude that there is no significant influence of the shape of the training set.

<figure><img src="image_21.png" /><figcaption>**Figure 20** Fréchet distances between the median crack path and the XPER simulations for different training sets: square set (■), pentagon set (●), and octogon set (●).</figcaption></figure>

**Evaluation of the computational time** We conclude this section by focusing in Table 5 on the CPU time associated to the crack prediction by our model for both training and test sets. By comparing to the XPER simulation, see Table 3, our model allows a drastic reduction of the computational cost (reduction by a factor 2575 for the median crack path when considering the same number of processors).

<table><caption>**Table 5** Predicted average computational time with the probabilistic model on a single processor.</caption><thead><tr><th></th><th>computational time [s]</th></tr></thead><tbody><tr><td>one iteration of Algorithm 1</td><td>0.13</td></tr><tr><td>one crack realization (M = 1)</td><td>13</td></tr><tr><td>median crack (M = 100)</td><td>1375</td></tr></tbody></table>

### 5.2.2 Analysis of a mechanical quantity: the tortuosity

For each microstructure, the tortuosities of the *M* predicted paths are computed and we focus on their median and uncertainty (confidence) intervals. This interval is obtained from percentile
estimate similarly to the previous section. Figure 21 displays the results when $M = 100$ for each microstructure of the test set and a comparison with the tortuosity coming from the XPER simulation is performed as well. This statistical quantity can be exploited to evaluate the probability for the tortuosity to belong to a given variation range. This type of information is relevant for computational cost reduction since it allows performing targeted XPER simulations i.e. simulations leading to a tortuosity in a variation range of interest. It is also important to keep in mind that the model provides more information than a median and a confidence interval. It allows to derive the tortuosity histogram of each microstructure, such as the one in Figure 22.

Finally, the performance of the model is studied in term of tortuosity density estimate taking into account the variability due to all microstructures. Figure 23 shows that the XPER tortuosity density can be accurately approximated.

<figure><img src="image_22.png" /><figcaption>**Figure 22** Tortuosity histogram for a microstructure of the test set.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 23** Tortuosity densities constructed from the XPER simulations (dashed) and the model predictions (solid) for microstructures with multiform aggregates. Density obtained by the kernel density estimate method. The mean tortuosity is respectively 1.181 for the model and 1.185 for XPER.</figcaption></figure>

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

# 6 Conclusion

We developed a probabilistic model for 2D fast crack prediction. This model leads to an algorithm that approximates at a low computational cost the local path of a crack.
The algorithm has been fully described and includes two steps. The first step is the computation of two local geometrical indicators (distance and angle) to capture the aggregate configurations in the vicinity of the crack tip. An efficient procedure exploiting a shadow cone has been proposed to evaluate them. The second step is the realization of a Markov chain model on the basis of the two indicators. The model parameters are estimated from a training set of numerically cracked microstructures. From any point of a crack in a microstructure, the model is able to evaluate the capability of any new point to be reached by the crack.

The numerical results show that, for a given microstructure with different aggregate shapes, the new model provides a set of predicted cracks that are in agreement with the computer code simulation. Moreover, it allows recovering the variability of mechanical quantities of interest such as the tortuosity.

The main advantage for practical issues is that the model is sufficiently flexible to be automatically adapted to any change in the local direction of the crack. A more complex application (a three-point bending beam test) is currently conducted on a mechanical test case where the crack direction is obtained through a mechanical simulation.

In this paper, it was decided to study microstructures of similar uniform density with different shapes of aggregates. Further investigations should concern different microstructures, in particular non-uniform densities (presence of clusters) and aggregates of different sizes. The first step will be to test the robustness of the current model on these new types of microstructures. According to the results, adaptation of the model should be considered, taking in mind that our goal is not to calculate a new training set for each new case.

## References

Akçaoğlu, T., M. Tokyay, and T. Çelik (2004). Effect of coarse aggregate size and matrix quality on ITZ and failure behavior of concrete under uniaxial compression. *Cement and Concrete Composites* 26(6):633-638. [DOI].

Bayar, G. and T. Bilir (2019). A novel study for the estimation of crack propagation in concrete using machine learning algorithms. *Construction and Building Materials* 215:670-685. [DOI].

Bichet, L. (2017). Mécanismes de transports dans la fissuration des matériaux hétérogènes : application à la durée de vie d'exploitation des centrales nucléaires. French. PhD thesis. France: Université de Montpellier. [HAL].

Bichet, L., F. Dubois, Y. Monerie, and F. Péralès (2016). Taking into account transport mechanisms in crack tracking problem of heterogeneous materials: application to the operating lifespan of nuclear power plant. *International Conference on Computational Methods - ICCM2016* (Berkeley, USA, Aug. 1, 2016-Aug. 4, 2016). [HAL].

Blal, N., L. Daridon, Y. Monerie, and S. Pagano (2012). Artificial compliance inherent to the intrinsic cohesive zone models: criteria and application to planar meshes. *International Journal of Fracture* 178(1-2):71-83. [DOI], [HAL].

Boissonnat, J.-D. (2018). *Géométrie algorithmique : des données géométriques à la géométrie des données*. Collège de France - OpenEdition Books. [DOI], [OA].

Chen, X., J. J. Yan, and H. Q. Yang (2012). Influence of aggregates on cracking sensitivity of concrete. *Applied Mechanics and Materials* 204-208:3299-3302. [DOI], [HAL].

Das, A. K., D. Suthar, and C. Leung (2019). Machine learning based crack mode classification from unlabeled acoustic emission waveform features. *Cement and Concrete Research* 121:42-57. [DOI].

Driemel, A., S. Har-Peled, and C. Wenk (2012). Approximating the Fréchet distance for realistic curves in near linear time. *Discrete & Computational Geometry* 48(1):94-127. [DOI], [OA].

Elices, M. and C. Rocco (2008). Effect of aggregate size on the fracture and mechanical properties of a simple concrete. *Engineering Fracture Mechanics* 75(13):3839-3851. [DOI].

Gangnant, A. (2015). Modélisation de la rupture quasi-fragile d'un béton à l'échelle mésoscopique. *34es Rencontres Universitaires de Génie Civil de l'AUGC* (Bayonne, France, May 26, 2015-May 29, 2015), pp 1-8. [HAL].

Husem, M. (2003). The effects of bond strengths between lightweight and ordinary aggregate-mortar, aggregate-cement paste on the mechanical properties of concrete. *Materials Science*
and Engineering: A 363(1-2):152-158. [DOI].

Jean, M., V. Acary, and Y. Monerie (2001). Non-smooth contact dynamics approach of cohesive materials. *Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences* 359(1789):2497-2518. [DOI], [HAL].

Jekel, C., G. Venter, M. Venter, N. Stander, and R. Haftka (2018). Similarity measures for identifying material parameters from hysteresis loops using inverse analysis. *International Journal of Material Forming* 12(3):355-378. [DOI].

Jeulin, D. (2015). Caractérisation morphologique et modèles de structures aléatoires. *Homogénéisation en mécanique des matériaux : Matériaux aléatoires élastiques et milieux périodiques*. Ed. by M. Bornert, T. Bretheau, and P. Gilormini. Vol. 1. Hermes Science. Chap. 4, pp 95-132. ISBN: 9782746201996.

Kim, H., E. Ahn, M. Shin, and S.-H. Sim (2018). Crack and noncrack classification from concrete surface images using machine learning. *Structural Health Monitoring* 18(3):725-738. [DOI].

Le, T. T. H. (2011). Étude multi-échelle du comportement thermo-hydro-mécanique des matériaux cimentaires : approche morphologique pour la prise en compte de la mésostructure. French. PhD thesis. France: Université Paris-Est. [HAL].

Le, T. T. H., F. Meftah, and H. Boussa (2009). A multi-scale thermo-hydral modelling of concrete behaviour at high temperatures. *Poro-Mechanics IV: 4th Biot Conference on Poromechanics* (Columbia University, New York, USA, June 8, 2009-June 10, 2009).

Le Mire, E., E. Burger, B. Iooss, and C. Mai (2021). Prediction of crack propagation kinetics through multipoint stochastic simulations of microscopic fields. *EPJ Nuclear Sciences & Technologies* 7(4):1-21. [DOI], [OA].

Monerie, Y. and V. Acary (2001). Formulation dynamique d'un modèle de zone cohésive tridimensionnel couplant endommagement et frottement. *Revue Européenne des Éléments Finis* 10(2-4):489-503. [DOI], [HAL].

Nordmann, J., K. Naumenko, and H. Altenbach (2019). Cohesive Zone Models-Theory, Numerics and Usage in High-Temperature Applications to Describe Cracking and Delamination. *Advances in Mechanics of High-Temperature Materials*. Springer, pp 131-168. [DOI].

Norris, J. R. (1997). *Markov Chains*. Cambridge University Press. [DOI].

Pele, K. (2021). Prédiction statistique de chemins de fissuration : application au matériau cimentaire. French. PhD thesis. France: École Centrale de Marseille. [HAL].

Pope, A. W. and H. M. Jennings (1992). The influence of mixing on the microstructure of the cement paste/aggregate interfacial zone and on the strength of mortar. *Journal of Materials Science* 27(23):6452-6462. [DOI].

Péralès, F. (2005). Fissuration des matériaux à gradient de propriétés : application au zircaloy hydrure. French. PhD thesis. France: Université de Montpellier. [HAL].

Péralès, F., S. Bourgeois, A. Chrysochoos, and Y. Monerie (2008). Two field multibody method for periodic homogenization in fracture mechanics of nonlinear heterogeneous materials. *Engineering Fracture Mechanics* 75(11):3378-3398. [DOI], [HAL].

Péralès, F., F. Dubois, Y. Monerie, B. Piar, and L. Stainier (2010). A NonSmooth Contact Dynamics-based multi-domain solver. Code coupling (Xper) and application to fracture. *European Journal of Computational Mechanics*, 389-417. [DOI], [HAL].

Safiuddin, M., A. Kaish, C.-O. Woon, and S. Raman (2018). Early-age cracking in concrete: Causes, consequences, remedial measures, and recommendations. *Applied Sciences* 8(10):1730. [DOI], [OA].

Scholer, C. (1966). The role of mortar-aggregate bond in the strength of concrete. *Symposium on concrete strength* (Washington DC, USA). [URL].

Socié, A. (2019). Modélisation chimio-mécanique de la fissuration de matériaux cimentaires : vieillissement et tenue des enceintes de confinement des centrales nucléaires. French. PhD thesis. France: Université de Montpellier. [HAL].

Socié, A., F. Dubois, and Y. Monerie (2021). Multibody approach for reactive transport modeling in discontinuous-heterogeneous porous media. *Computational Geosciences* 25:1473-1491. [DOI], [HAL].

Walker, S. and D. L. Bloem (1960). Effect of aggregate size on properties of concrete. *ACI Journal Proceedings* 57(9). [DOI].
Willot, F. (2015). The power laws of geodesics in some random sets with dilute concentration of inclusions. *Lecture Notes in Computer Science*. Springer, pp 535–546. [DOI], [HAL].

Xiao, Y., Z. Wei, and Z. Wang (2008). A limited memory BFGS-type method for large-scale unconstrained optimization. *Computers & Mathematics with Applications* 56(4):1001–1009. [DOI], [OA].

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the authors—the copyright holder. To view a copy of this license, visit creativecommons.org/licenses/by/4.0.

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

**Authors' contributions** KP: Conceptualization, Methodology, Software, Visualization, Writing - Original Draft, Review and Editing. JB: Conceptualization, Writing - Review and Editing. LD: Conceptualization, Writing - Review, Supervision. JL: Conceptualization, Writing - Review, Supervision. TLG: Conceptualization, Writing - Review. YM: Conceptualization, Writing - Review and Editing. FP: Conceptualization, Writing - Review

**Supplementary Material** None.

**Acknowledgements** None.

**Ethics approval and consent to participate** Not applicable.

**Consent for publication** Not applicable.

**Competing interests** The authors declare that they have no competing interests.

**Journal's Note** JTCAM remains neutral with regard to the content of the publication and institutional affiliations.