Clustering analysis in a mesoscopic model of a concrete-like material and its impact on dynamic fracture

– In this paper, a numerical tool to evaluate the quality of a meso-structure is presented. This tool is then modiﬁed to produce a local measurement of the element density in a matrix (local clustering). Finally, from this local quantity, a method is proposed to predict the crack nucleation and propagation in concrete-like material under tension loading. The predictions that can be obtained with this method are compared to results from numerical simulations of dynamic traction using cohesive elements in Akantu Finite Element code on various meso-structures.


Introduction
Failure of concrete-like materials, is a complex and non-linear phenomenon related to their heterogeneities, composition, geometry and loading conditions. From a local point of view, fracture of these materials involves at the beginning of the loading the opening of micro-cracks in the matrix which propagate and may coalesce in interaction with the aggregates. Explicitly representing the meso-structure allows then to study its influence on the mechanical properties of the different phases, and their interfaces.
A meso-scale level of observation (as proposed in [1]) allows an explicit representation of some concrete constituents, which enables reducing the number of model parameters and to describe the interactions between matrix and inclusions. This also allows to describe concrete as a biphasic material: aggregates embedded in a mortar paste matrix, to better understand their role in fracture mechanisms. It is indeed known that representing granular meso-structures allowed some authors to represent certain macroscopical effects, such as in [2,3].
Several assumptions are usually made on the role of inclusions in the dynamic fracture mechanisms: creation of weak zones in the structure and local stresses increase, changes in crack paths, important for representing scattering and size effects.
Being able to draw such conclusions is highly dependent on the quality of the numerical representation of concrete a Corresponding author: pierre@lmt.ens-cachan.fr meso-structure. The aim of this paper is thus to present the method of representation used, and some analytical tools that were developed to study the meso-structure impact on mechanical properties. The numerical simulations used to represent the fracture will be detailed and the resulting cracking patterns will be compared to predictions derived from the proposed meso-structure analysis tools.
The paper is composed as follows. Section 2 describes the meso-scale approach and choices made. In Section 3 the tools for a geometrical analysis of a meso-structures are presented. With its material parameters is presented the chosen finite element framework for representing compressive cracking at the interfaces. Comparison of these geometrical tools with numerical simulation of tensile tests is reported in Section 4. Finally, concluding remarks are stated in Section 5.

Concrete representation
Concrete is a heterogeneous (quasi) brittle material made of various components, which are present in different proportions. The failure of this kind of materials is a complex and nonlinear phenomenon, which dissipates energy, according to its (meso-structural) composition, geometry and loading conditions. Fracture of these materials involves the opening of local micro-cracks, which may propagate, coalesce and subsequently opposing crack surfaces enter into contact influencing the nonlinear failure process. Nevertheless, it is commonly assumed that cracks tend to initiate at the interface between the cement paste and the aggregates, which can be seen more or less as hard Article published by EDP Sciences S. Pierre and F. Gatuingt: Mechanics & Industry 18, 306 (2017) inclusions. To explain this fracture behavior, one can note that: -The Young's modulus of cement paste and aggregates are different, which leads to strain incompatibilities, a porous layer with a high water content develops on the aggregate surface. This thin layer is called Interfacial Transition Zones (ITZ) (see [4]) and its lower tensile strength explains local weakness of the composite material, the inclusions shapes and their spatial arrangement induce local stress concentrations.
Moreover, when the first crack initiates, inclusions distribution can affect its coalescence and propagation. In dynamics, an important aspect to take into account is that the inclusions introduce heterogeneity in the stress field. This, then, may induce a transition from one main crack in static tension to multiple cracks if the strain rates are high enough.

Meso-structures generation techniques
Experimental observation on real concrete morphologies may be summarized as follows: the aggregates in the cement paste never are in contact, unlike what might be found for granular media. This is due to hydrates that appear rapidly on solid inclusion surfaces during hydration [4], the concrete microstructure is complex and a thorough representation of every phase (sand, anhydrous cement particles, small aggregates) would be extremely costly from a numerical point of view and would also drastically extend the number of input parameters.
To preserve the general aspect of this numerical study on a controlled material, it was decided to only represent a biphasic material: aggregates embedded in a mortar paste matrix. In our model only the largest aggregates are represented explicitly. Whereas, the smallest aggregates and other components are assumed to be mixed up with the cement paste establishing the matrix phase. In this case, the choices made with regards to the generation and placement of the inclusions are important when it comes to representing the macroscopic behavior of concrete. A large range of inclusions shapes might be found in the literature to numerically represent the aggregates in the mesoscopic representation of concrete. Among those, one may cite the most simple spherical inclusions [2,5,6], the enriched ellipsoidal inclusions [7], or the most pertinent (but complex and requiring much data) Fourier descriptors-generated inclusions [8].
Due to the contact between particles experimentally observed in granular media, it appears unjustified to use algorithms coming forth from this research field. In concrete meso-structures case, the most common solutions for inclusions placement are to use the Random Sequential Adsorption (RSA) algorithm [5,15] or based on Voronoi's polygons [16].
In our study, for the sake of simplicity, the choice was made to use circular inclusions placed with the RSA algorithm in 2D, detailed hereafter.

Circular inclusions generation
The most widely used approach to generate aggregate shapes in a meso-scopic description of concrete is to simplify them as much as possible. Many authors chose then to represent elliptical or spherical inclusions [2,5,6].
Inclusions radius are often computed to fit a gradation curve, either by estimating it with a probability law, or by discretizing it in several classes of inclusions of identical dimensions as in [2,17].
The main benefits of this type of approach are that: the numerical implementation is fairly simple and allows for fast results, this simple shape for the inclusions has been chosen for some controlled materials for which experimental results exist.
However, there are some drawbacks: stress concentrations are underestimated, involving probably a delay in the fracture process, the shape allows for large rotations when a total decohesion of a rigid inclusion is obtained, the inclusions does not look like real aggregates in concrete, raising issues concerning the representation of the results therefore obtained.
Nevertheless, one can notice that this simple approach is interesting in regard with many classical homogenization techniques which are well-adapted to these kind of shapes -e.g. Mori-Tanaka [18] homogenization.

Placing inclusions with the RSA algorithm
The Random Sequential Adsorption (RSA) method to place inclusions in a matrix is the most commonly used in literature. This can be explained by its simplicity and the little amount of information needed on inclusions real repartition. The RSA method is also easy to set up numerically.
The basic idea behind this algorithm is to place each inclusion sequentially and randomly. Placement of an inclusion is validated only if it does not overlap a previously placed. As this method is widely discussed in the literature, the morphological characteristics of resulting dispersions are well known [19][20][21]. This simple placement procedure yields to quick results and has met with success, even in recent studies [ Though this method is widely used for circular inclusions generations, it can be also used to place another type of shape [7,23]. In this case, it should be noted however that one loses the main advantage of this algorithm: its simplicity. Indeed, it may be very numerically costly to detect overlapping with complex shapes for the inclusions.

Mesoscopic representation choices
One of the aim of this work is to study the influence of the inclusions placement on the mechanical properties and fracture mechanisms in concrete-like material.
As presented earlier, we decide to use circular inclusions in our mesoscopic model. Considering circular shapes for the aggregates introduces a known bias in the macroscopic mechanical behavior due to the fact that this model will underestimate stress concentrations. Nevertheless, with this simplified approach, one can easily focus on the quantitative analysis of the influence of the aggregate placement on this response.
As discussed in the previous section, the RSA algorithm is disconnected from the physical mechanisms that occur during concrete casting and hydration. This is an acceptable compromise mainly due to the number of additional parameters otherwise needed and to the large number of authors using it for concrete mesoscopic representation.
The kind of morphology one might expect with the RSA method presented here is shown in Figure 1. This aggregate distribution was obtained using the realistic granulometry given in Table 1 for a 15 × 15 cm sample .
For further mesh generation, several aggregate classes of radius lesser than 4 mm (and accounting for 32% of the aggregate surface fraction) was not represented here because we consider that small inclusions will have a less important impact on dynamic fracture. The importance of this assumption should probably be studied more in detail in a future work if we consider the recent work of [24] in quasi-statics.

Geometrical analysis of mesostructures
The main disadvantage of random generation of mesostructures is that it might be difficult to study the resulting distributions. As previously mentioned, it might be useful to check the quality of the meso-structures generated. Moreover, being able to link informations obtained from the spatial arrangement of the inclusions to the results obtained by numerical simulations of a tensile test (e.g. the cracking mechanisms) is very interesting and challenging.

A measure of the aggregates distribution quality
In order to determine the generated meso-structures quality, a mathematical tool is proposed: Ripley's function.
This function introduced by Ripley [25] is a spatial indicator close to the functions nearest-neighbor distribution [26] or pair correlation [27]. It is expressed as: where: -N (t) refers to the number of events in a radius t from a randomly chosen event (illustrated in Fig. 2); λ is the density of events.
This measure introduces an influence radius, t. If the interval between two events is greater than t, it is assumed that they will not influence each other.
In a context where the events (inclusions in our case study) distribution is perfectly isotropic and follows a uniform distribution law, it can be proved [28] that: Any deviation from this theoretical value can be seen as a flaw in the quality of the meso-structure analyzed. To numerically assess to this quantity, one may use an esti-matorK(t) of K(t) that can be expressed in a population of elements e as: where: d i,j , is the distance between events e i and e j (inclusions/aggregates in our case); -I (condition), is an indicator equals to 1 if the condition is verified, 0 otherwise; t, is the influence radius; -N elements , is the number of mesostructure elements e.
In order to simplify the analysis of the distribution quality, a functionL(t) was introduced by Dixon [28] to easily detect any deviation from the ideal value K(t) = πt 2 : Thus,L(t) − t should be equal to 0 for an ideal Poisson's process. Plotting it indicates the deviation from the uniform and isotropic case. This method is very efficient to detect correlation lengths [29]. In order to test it, two distributions were generated in a 1 × 1 m square sample: in the first one, elements have been placed two by two with a fixed distance of 0.01 m. This generates a meso-structure with one correlation length in the spatial arrangement of the aggregates. in the second, elements have been placed three by three, forming an isosceles triangle which sides are 0.01 m and 0.015 m respectively. In this particular case, two correlation lengths are introduced.
In both cases, 10 000 elements were placed in the sample tested in order to obtain results independent of the drawing process. Indeed, a important variability of the measure is expected and observed for a small number of elements.
In Figure 3 theL(t) − t function is plotted for the two placements generated. One can observe in this figure the effect of the introduction of a correlation length in the meso-structures. Each correlation length produces a peak in the plot linked to how often this length can be detected in the distribution analyzed.
Due to the introduction of a radius t in the analysis, a major drawback is that the Ripley's functionK(t) is affected by edge effects. To overcome this problem, Dixon [28] proposed to introduce a weigh coefficient w in this function: where the weight w added in the original expression of K(t) can be used to minimize the edge influence on the result. In our work we proposed to use a weight function defined as: where x ei is the spatial position of particle e i ; -S real i (t) is the surface of the circle of radius t centered on particle e i that is in the surface of study. The efficiency of this correction was tested for a uniform distribution of 10 000 events in a square surface. The results obtained are presented in Figure 4. It appears that using w edge eliminated the progressive deviation of L(t) − t indicator from the expected theoretical value. The 0 value, expected for an infinite and uniform distribution, is reached with a good approximation in our case with only 10 000 points. We can then conclude that the expression proposed in Equation (5) with our weight function can be a good measure of the aggregates distribution quality in a meso-structure.

A clustering indicator
As shown in the previous part, K(t) and its estimator only provide global informations on the distribution quality. They do not account for inclusions size or type and do not consider what happen locally in specific places of the meso-structure.
The main interest in this work is to better understand phenomena such as crack nucleation and propagation in heterogeneous material. Assessing meso-structure's impact on these mechanisms is yet difficult to a priori quantify. Indeed, inclusions spatial arrangement has a huge influence on the stress concentration field and in concrete-like material the interface between the aggregates and the matrix is a weak zone. All these points should have an impact on the fracture of such materials.
To try to assess this point a local clustering indicator is proposed, deriving from the Ripley's function. Local term of theK(t) expression is then used: However, with this definition, the function is only defined and measured on the set of inclusions e i , which provides a quantity that is non-representative and difficult to handle.
To address this representation issue, another definition is provided for the clustering indicator, that is more general and defined at any point of the meso-structure: Several mappings of this clustering indicator are presented in Figure 5. They correspond to the same placement of inclusions in a square surface. They were computed using the same distribution of particles obtained with the RSA algorithm, which should ensure isotropy as said in part 2.1.2. In this figure, the higher the local value of C(x, t) is, the darker the color on the map. For each mapping, a different value of the influence radius t has been used for the computation. As we can see in Figure 5 the indicator proposed is not smooth and highly depend on the radius of influence which does not seem very appropriate. Meaning that if the influence radius chosen is not good enough, no correlation between the apparent cluster positions and the values of C(x, t) can be found. If the influence radius chosen is too small compared to inclusions size, no cluster would be detected as shown in Figure 5a. On the contrary, if the radius is too large compared to the cluster size (Fig. 5c), the indicator would detect important clusters that does not exist between the neighboring clusters.
To provide continuity and smoothness to this indicator, a regularization is then introduced. The form chosen is a classic Gaussian weight [30] used in nonlocal analysis: Several mappings of this regularized clustering indicator are presented in Figures 6 and 7. In this new analyse, the meso-structures are a 10 cm square where 1 cm diameter inclusions cover about 40% of the surface for two aggregates distributions. The influence of the radius t is also observed in theses figures.
We can observe that this regularized indicator is also dependent on the value of the chosen influence radius t as in Figure 5 for the non-regularized indicator. Nevertheless, even if the clusters seem to be larger with the increase influence radius t, the zones where the value of the regularized indicator is high (dark red in Figs. 6 and 7) remain the same. In this regularized case, the results are clearly less dependent on the influence radius t. From our experience, a radius equal to the size of the aggregates seems to be the better choice. It might be noted that the measurements presented here are not really impacted by the edge effects because the particles at the edge that are not entirely inside the square shape of our specimen are taken into account for the calculation of C(x, t).

Comparison with numerical simulations
To verify the relevance of the proposed indicator (Eq. (9)), predictions drawn from the clustering indicator are compared to the cracking patterns obtained through numerical simulations of dynamic tensile test on a concrete-like controlled material.

Numerical framework
To simulate the meso-structure's behavior, the Finite Elements Method was used. Fracture is modeled by means of cohesive elements and the extrinsic approach developed by [31,32] that allows to circumvent the mechanical properties alteration observed when interface elements are allowed to deform before fracture [33]. Interface elements, representing the cohesive separation, are introduced during the simulation when a criterion is verified, here in tension. The Akantu code [34] in which this method is implemented was used for the simulations.
The traction separation law used is the most classical one developed by Camacho and Ortiz [35] and is presented in Figure 8. It links the loss of strength as the damage progress to the effective distance (scalar) δ between the crack edges. This numerical method is presented in details in [2] or [16].  [35] (source [16]). For all following simulation results, the same material parameters has been used and are presented in Table 2. They define the linear behavior of the mesh elements.
The parameters controlling the traction separation laws are presented in Table 3. The role of each parameter is given in [2], but one can simply explain that they are responsible for the evolution of the damage and the behavior of the interface elements that are inserted during simulation to represent the cohesive separation of the materials.

Comparison with numerical simulations
Predictions made using the analytical tool developed in part 3.2 has been compared to the cracking patterns obtained through numerical simulations presented in the previous part.
The meso-structure that is used for the simulation is a model material with only one inclusions class of d = 10 mm diameter and accounting for 40% of the surface of a 10×10 cm square as used in Section 3 for the geometrical analysis. The sample is loaded under displacement control with a constant strain rate (ε = 1 s −1 ) .
In tension all the nodes of the finite element mesh which are located on the upper (respectively lower) boundary are forced to move at a constant velocity V 0y = V 0 (respectively V 0y = −V 0 ) as illustrated in Figure 9. To avoid important stress wave propagation and an early fracture near the boundaries [36], an initial velocity is prescribed to all nodes of the finite-element mesh as illustrated in Figure 9.
The result of the geometrical analysis of the mesostructure is shown in Figure 10a. The clustering indicator x y Fig. 9. Initial gradient and boundary conditions. mapping shows high values on the edges, such as what appears in the middle of the left side of the specimen. This zone displays the highest values of C(x, t) and is close to an edge, where macro-cracks are expected to initiate during tensile tests. Apart from this one, three high values zones may also be observed: top left hand corner, top right hand corner and at the bottom of the specimen. Even if these zones are lower in values than the previous one, they show the presence of spatial heterogeneities that might affect the cracking behaviour of the meso-structure. Figure 10b shows the fracture patterns of the sample that were produced under tensile test. The damage variable is associated to the gradual decrease of strength of the cohesive elements. It is defined as D = min(1, δmax δc ) (see Fig. 8).
If one considers cracking along the x-direction, which is expected for this type of loading, a band of high values of the C(x, t) indicator can be observed in Figure 10a at middle height in the specimen, where the main crack appears in the numerical simulation (Fig. 10b). Moreover, on the opposite edge of the specimen in the upper part, where relatively high values of the indicator where close to the edge, a crack started to open before being shielded from the main one. This can still be seen in the slight difference of displacement where elements where damaged.
For the two others zones of higher values (on the top left hand corner and at the bottom of the specimen), no crack initiated even if our indicator did detect clusters of aggregates. This is due to the vertical tensile loading that promote the macro-crack initiation at the left and right edges. These zones will be probably active in horizontal tensile loading that, in this case, promote the macro-crack initiation at the upper and lower edges.
For this first example, crack initiation seems to be closely correlated to our clustering indicator. In order to confirm this, another meso-structure was tested and Figure 11 shows the results of the indicator mapping and of the dynamic traction simulation.
In this second example, the zones with the highest values of our clustering indicator are at the top in the middle of the specimen, in the center of the specimen and on each edge at the bottom.
Once more, crack initiation appeared in areas where the element density was high: two fully formed cracks appeared in the middle of the specimen and on one edge on the bottom. Moreover, some crack initiated at the top left hand corner and on the bottom right hand corner where our indicator displayed high values, but these cracks were soon shielded by the two macro-cracks and stopped propagating. Once again due to the direction of the tension, the high value of the indicator on the upper side does not induce a crack.
In our simulations, the strain rate imposed (ε = 1 s −1 ) is small enough to produce localized macro-cracks. In [2] it was shown that a local inertial effect induces a more diffuse micro-cracking crack pattern for higher strain rates. In this more diffuse crack pattern, our indicator will lose its sense.

Conclusion
In this paper, we proposed a geometrical tool to evaluate the quality of an aggregate distribution in mesoscopic model of concrete-like materials. This model is based Ripley's function. Based on this global tool we proposed a new local indicator to provide a measure of the aggregates clustering in a meso-structure. This indicator depends on a single parameter t which is a king of radius of influence. This indicator had the disadvantage to give clustering map that is extremely dependent on the value of t. In order to overcome this issue, a regularized indicator has bee then developed. The results obtained are then less dependent on the value of t and we can now assume that this new indicator is more objective.
The second main interest of the paper is to compare the geometrical clustering map to the crack pattern obtained from a direct tensile test. The results obtained are hopeful. In this simple dynamic loading test (with small effect of local inertia), one can easily link the crack path to the position of the aggregate clusters in a mesoscopic representation of a heterogenous quasi-brittle material. Nevertheless, some work remains in order to capitalize the results obtained for complex loadings. This indicator should probably be coupled to the local stress states to succeed in a more accurate prediction of cracking, in particular for rotating paths in the loading state.
One might also keep in mind that weight functions have not been exploited in the calculation of the local indicator so far, and that many other weighting possibilities exist. It might be interesting to modulate the influence of different types of particles on one another, depending on their nature, position or size.