{'fraction_non_alphanumeric': 0.024629106718658957, 'fraction_numerical': 0.005131997669311102, 'mean_word_length': 4.59816836030611, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Studying the influence of inclusion characteristics on the characteristic length involved in quasi-brittle materials using the lattice element method

Unlike nonlocal models, there is no need to introduce an internal length in the constitutive law for lattice model at the mesoscopic scale. Actually, the internal length is not explicitly introduced but rather governed by the mesostructure characteristics themselves. The influence of the mesostructure on the width of the fracture process zone which is assumed to be correlated to the characteristic length of the homogenized quasi-brittle material is studied. The influence of the ligament size (a structural parameter) is also investigated. This analysis provides recommendations/warnings when extracting an internal length required for nonlocal damage models from the material mesostructure


Introduction
Fracture of quasi-brittle materials is characterized by a zone with a finite size around and ahead the crack tip, in which damage occurs and causes the softening behavior of the materials. This is the fracture process zone (FPZ). For instance for concrete, the size (width) of the FPZ, denoted by FPZ hereafter, is believed to be proportional to the maximum aggregate size d max , see, e.g., [1,2]. Therefore, in nonlocal models (gradient or integral form [2][3][4]), the FPZ size which only depends on the internal length c introduced, depends on (is proportional to) the maximum aggregate size. Accordingly, neither loading nor structural effect is considered to affect the resulting size of the FPZ except in the latest integral nonlocal model proposed in [3]. In the latter, the internal length parameter evolves depending on the stress state during the damage process and also depends on the intrinsic (characteristic) length that can be correlated with aggregate size of the material. However, the correlation between the characteristic length and the aggregate size has not been explicitly calibrated yet.
The literature often reports a linear or affine relation between c and d max , see, e.g., [5,6]. But actu-ally, varying d max in experiments may lead to a number of changes in the aggregate structure characterized by other parameters such as the volume fraction of aggregate, their size distribution, their fabric or connectivity. Basic questions may be raised: what does affect the internal length of a nonlocal model? Is it only the maximum size of aggregates or some less obvious parameter(s)? Does the structure itself (size or ligament) play a role in the internal length?
To address these questions, numerical simulations of uniaxial tensile tests are carried out using the lattice model in which the geometry and mechanical properties of the material mesostructure are explicitly introduced. The output of the simulations is the FPZ size and the characteristic length of the material. The characteristic length is a priori regarded as the internal length that would be introduced in nonlocal models. The same notation c is thus used in the following. From the lattice simulations, the relationship between the two lengths FPZ and c , and some relevant characteristics of the material mesostructure is found out. The study is restricted to the case of two-dimensional analysis of a brittle elastic model material with circular inclusions and is also restricted to mode-I failure problems occurring with small deformations under quasi-static loading conditions.
It is important to stress before reading the following that the inclusions and matrix have a brittle elastic be-havior together with highly simplified geometry. As a consequence, our observations and conclusions must be translated with caution to the case of real concrete.
The lattice model used in our study is briefly recalled. The model is implemented in our self-writing code using C++ programing language. The method to assess the FPZ size and the characteristic length of the material will be next pointed out before performing numerous numerical experiments to study the influence of the material mesostructure and of the structural parameter (ligament size) on these lengths.

Numerical model
The lattice element method (LEM) is a convenient way to model the fracturing of quasi-brittle materials for the problems in which the discontinuities are dominant since it provides a discrete representation of material disorder and failure. By using the LEM, the microcracking, crack branching, crack tortuosity and bridging of quasi-brittle materials can easily be identified and captured. It allows the fracture process to be followed until complete failure. There exist two different types of lattice models. The first one is called classical lattice models in which the material is discretized as a network of discrete 1D-elements that can transfer forces and possibly moments [7][8][9][10]. The second type of lattice models, called particle lattice models, are classified as a discrete element method [11] in which the material is discretized as an assemblage of rigid particles interconnected along their boundaries through normal and shear springs [12]. The models in this category also include the rigid-body-spring networks [13], bonded-particle model [14], random particle models [15], beam-particle model [16,17], confinement-shear lattice model [18]. The main advantage of particle lattice models with respect to classical lattice models is that they account for the fact that crack surfaces may act on each other causing the repulsive force during the loading process. So the particle lattice models are more suitable for predicting the failure behavior in mode II or mode I under cyclic loadings whereas the classical ones are enough when the mode I failure prevails.
In this work, only the mode-I failure of the material submitted to monotonic mechanical loadings is considered. Moreover, for studying the influence of the material mesostructure on the FPZ which is related to the characteristic length of the material, a detailed description of tortuous crack patterns is important. Therefore, a lattice model, based on the classical lattice models, in which the normal and shear springs are introduced.
The constitutive laws of the 1D-elements are simple elastic relations in the normal and tangential directions defined by each element, see Figure 1a. Only small perturbations are considered, the positions of the lattice nodes are assumed fixed and unknown variables are the node displacements u. The axial direction n i j 0 and transverse direction t i j 0 associated with each element i j remain thus fixed. Length variations between the node i and j are defined by δ i j n = (u i − u j ) · n i j 0 and δ i j t = (u i − u j ) · t i j 0 for the normal and tangential directions, respectively. The forces are related to this length variations by f i j n = K i j n δ i j n and f i j t = K i j t δ i j t , where K i j n and K i j t are the normal and shear stiffnesses of the element, respectively. Figure 1: 1D-element with its local coordinate system (a) and its effective width A i j (b).
The approach consists in finding the set of node displacements [u] -among which some are imposed along the boundaries -that minimize the total elastic energy of the system: To proceed this minimization, the conjugate gradient method is used with the following definition of the gradient: where e α stands for the two directions of the global frame. The damage (in the form of diffuse or macroscopic cracks) of the whole lattice system is accounted for by removing each element that breaks according to a criterion ψ( f i j n , f i j t ) ≥ 0. The Mohr-Coulomb surface with a cut-off of the tensile strength [19] can be adopted. However, we chose to use another model that has the advantage of being more generic while it is expressed in a single function: where σ 0 n and σ 0 t are the ultimate stresses for pure normal and tangential loadings, respectively; n is a positive parameter that changes the yield surface from a linear form (n = 1) -corresponding to the classical Mohr-Coulomb criterion -to a non-linear form (n > 1). In this study, n = 5 is used.
Let us now consider a system of lattice elements where small displacements are imposed for some nodes on the boundary. A reference solution [u ref ], corresponding to the free displacements of the other nodes, can be found by minimizing U e as described above. Provided that the elements remain elastic and intact, any other elastic state is an uniform scaling of the reference solution: As a consequence, elastic forces can be scaled by the same factor and it becomes possible to find, for each element, a factor η i j so that ψ(η i j f n , η i j f t ) = 0. The state corresponding to the failure of the weakest element can thus be obtained by scaling the reference solution by the factor η min = min i j {η i j }, and then recorded. The next loading state will result from another reference solution beginning from a new configuration into which the broken element is removed. By repeating this procedure for each element failure, one by one, the loading course is controlled by these events, rather than a time-stepping which could involve more than one element removal within a single time step. This would results in nonphysical solutions that make the mechanical response dependent on the loading magnitude [20].
With the LEM, heterogeneities appear de facto at the mesh level. The required disorder in the mesh, introduces a variation in lengths i j and effective width A i j of the elements . It results in an unwanted parasitic heterogeneity in the stiffness properties that can be limited by accounting for the local geometry in the element behavior: whereK ϕ n andK ϕ t are the stiffnesses that can be set "uniformly" to a phase ϕ. The effective width A i j is the distance between centroids C i jk and C i jm of the triangles adjacent to the element i j, projected onto the local direction t i j 0 as proposed in [21], see Figure 1b. Since the state of plane stress or plane strain is not explicitly defined in LEM-based simulations, the quantity A i j can also be regarded as a surface by assuming an unit length in out-of-plane direction. In this picture,K ϕ n and K ϕ t have a dimension of material stiffness. As a consequence of the weighting of imposed stiffnesses (or modulii)K ϕ n andK ϕ n in a phase, actual stiffnesses of elements differ from each other.
The targeted Young's modulus E ϕ and Poisson's ratio ν ϕ of the phase ϕ can be used to determine the element stiffnesses by the following relations: These relations are derived from the equations given in [22] for a regular and triangular lattice, by replacing a factor √ 3 by 1 (found empirically from a number of single-phase simulations).
From there, heterogeneity intrinsic to the mesh geometry is limited as much as possible, and a structure of inclusions (grains) can be generated using the take-andplace processes [23,24]. After generating the inclusion structure, different material phases are defined and different local mechanical properties are assigned to the elements falling in each phase. At the mesoscale, three phases can be distinguished: inclusion, matrix and interfacial transition zone (ITZ) , see Figure 2. If both ends of an element are located in the same phase, then this element is assigned the same mechanical properties of the corresponding phase (inclusion or matrix), otherwise it is considered as interface or inclusion element depending on the location of its midpoint. If its midpoint is located within the grain, the element is classified as inclusion element, or else it will be ranked as ITZ element. The reason for this definition of ITZ element is that the resulting fraction of inclusions (the ratio between the number of inclusion elements and the total number of elements) is closer to desired fraction of inclusions in material than those developed by other authors [7,25,26]. In their models, all elements that connect two different zones of grain structures are considered as ITZ elements. inclusion ITZ matrix Figure 2: Distinction between inclusion, matrix and ITZ phases according to the location of a lattice element in the grain structure.

Assessment of characteristic length
To account for damage in continuous (and homogenized) modeling of concrete, a length parameters is required [2,4]. This length, denoted by c , called characteristic length is seen as an intrinsic property of the material, however it is not so simple to determine and to connect with the heterogeneities at lower scales. The method proposed in [6] is used here to assess this characteristic length for a material modeled by lattice elements. The basic idea is that the characteristic length of the material is approximated by the effective width h of the zone in which the fracture energy of the material is dissipated. This effective width is defined as the ratio of the fracture energy G f (energy per unit area of crack surface) dissipated by the cracking that localizes in a narrow band of the specimen in localized tensile test and the energy density W s dissipated by the cracking that is nearly homogeneously distributed in the whole volume of the specimen of the same material in distributed tensile test. Finally, the characteristic length is approximated by h which can be assessed by To evaluate c with LEM simulations, both numerical tensile tests (localized and distributed) have to be performed to determine G f and W s . G f is determined from the tensile test performed on a notched specimen so that the damage can be localized whereas W s has to be determined from the tensile test carried out on an unnotched specimen with specific design of loading such that the damage is as homogeneously as possible distributed in the specimen volume. To this end, the numerical simulations of tensile tests using the lattice model can be performed in which the tensile loading is indirectly applied to the notched and unnotched specimens by elongating the steel bars "glued" to the specimens as proposed in [6], see Figure 3. These two tests were performed on numerical specimens of the same size, the loading is applied by means of lateral bars that are "glued" to the specimen and set 10 times stiffer than the material tested. The main difference between the two types of tensile tests is that the steel bars are only glued to the ends of the notched specimen within a certain length while they are entirely glued to the unnotched specimen within the whole height of the specimen. For the following, the tensile tests performed on notched specimens, where the Location of Damage is forced, are referred to as LD-tests. The tensile tests performed on unnotched specimens, designed to identify Distributed Damage, are mentioned as the DD-tests. These tests are known as the PIED (Pour Identifier l'Endommagement Diffus) tests in the French community, as introduced in [27]. Note however that a diffuse damage is actually not achievable, that is why we prefer to talk about distributed rather than diffuse damage. In the lattice simulations, the steel bars with the width of 2 mm are also discretized in 2D by the lattice elements (2D mesh) but their stiffnesses are set 10 times greater than those of the material tested and they always have an elastic behavior. The steel bars are perfectly "glued" to the specimens via compatible nodes. Figure 3: Sketch of the specimens used to determine the characteristic length as proposed in [6]. The tensile test performed on the notched specimen to obtain the localized cracking process (a) and on the unnotched specimen to obtain the distributed cracking process (b).
In LD-tests (Figure 3a), a crack is initiated and then propagates until the specimen breaks. The fracture energy G f is simply the sum of all elastic energy dissipated by the rupture of each element i j divided by the total cracking surface: The DD-tests ( Figure 3b) aims to avoid any onset of crack so that the straining and damage are as uniform as possible. The energy density W s is thus given by the total elastic energy dissipated within the specimen volume V: Direct measurement of the effective width h, denoted by FPZ for the following, of the fracture process zone (FPZ) is another characteristic dimension. We also made this estimation from single tensile tests performed on notched specimen, by treating the fracture energy of each element similarly to acoustic emission [28,29]. A density map of the dissipated elastic energy can be drawn from broken elements. Based on this map, the size of the FPZ can be determined by analyzing the density distribution of dissipated energy around the macrocrack. This distribution, when represented as a probability density function (pdf), can be fitted by a Gaussian distribution in order to extract a width. Rather than that, we choose to rely on the cumulative density function (cdf) of the dissipated energy to determine the size of the FPZ since that curve can be more smoothly defined by sorting the dissipated energy along a direction. The direction chosen here is the one perpendicular to the mean direction of the final crack which may not be strictly perpendicular to the loading direction depending on the microstructure setting. A fit of the cumulated form by a "Gaussian bell" allows to assess FPZ as being four times larger than the standard deviation σ of the Gaussian curve. This choice corresponds to a width containing a bit more than 95% of energy dissipated (provided that only one process zone exists).
It is worth pointing out that the FPZ size and the characteristic length of the material determined by lattice simulations also result from the mesh size, i.e. the lattice element size. This means that the LEM introduces a characteristic length by its mesh. An analysis of the mesh-size influence on the FPZ size is performed. A series of LD tensile tests is performed in which the specimen is discretized with five different mean values of the mesh size l m . Furthermore, for each discretization, five independent meshes are generated by randomly moving the nodes within the radius of l min (the minimum mesh size) to take into consideration mesh orientation effect on FPZ . The dependence of the PFZ size on the mesh size is shown in Figure 4.l m is the average value of the discretization size . As expected, the FPZ size does statistically tend to "zero" upon mesh refinement. Note however that the intercept of the fit is not exactly zero, its value is 0.18 mm. This is probably due to the fact that there are only five discretizations were used and there was not any mesh finer than 1 mm to be generated for the sake of saving computational time. Once the influence of the mesh on the material internal length is known, it can be subtracted from the relationship between the internal length and the inclusion properties. The latter defines the aim of the present study.

Numerical experiments
To study the role played by coarse inclusions in the internal length, a number of simulations have been performed. In the modeling of the material mesostructure, inclusions are considered, which are embedded in the matrix separated by the interfacial transition zones (ITZ). The inclusions, matrix and ITZs are assumed to be linear brittle elastic. The inclusions are also assumed to be stiffer and more resistant than the matrix, whereas the ITZs are assumed to be less stiff and with a smaller strength than the matrix. In the following simulations, the stiffness and the strength of inclusions are 10 times larger than those of the matrix. In turn, the stiffness and the strength of the matrix is 2 times larger than those of ITZs. Elastic and strength parameters of the matrix are set to values listed in Table 1, and they are kept fixed for all simulations. The way coarse inclusions are structured -referred to as "grain structure" in the sequelwas restricted to two characteristics in this study: the mono-sized grain diameters d and their surface fraction P a . In the (P ad) parameter space, shown in Figure 5, three variation paths were considered:  In addition to the variation of grain structure, the presence of a weak interfacial transition zone between inclusion and matrix phases is analyzed. Without weak ITZ, only two phases (2ϕ) are modeled in the sense that the properties of the ITZs defined as in Figure 2 are those of the matrix. With weak ITZ, a less stiff phase with smaller strength is added in-between inclusions and matrix, bringing the number of phases to three (3ϕ).
Typical force-displacement and stress-strain curves obtained for the LD and DD tests, respectively, are shown in Figure 6. The corresponding crack patterns are also presented. It is seen that there is only one macro crack which traverses the notched specimen of the LD test while about fifteen macro cracks cross over In fact, the last one is characterized by a monotonic increase of the displacement. Therefore, in order to have a corresponding response, the "envelope" of the numerical curve should be taken. The envelope curve is obtained by the so-called smoothing procedure. The procedure is described as follows. By connecting from the first to the last point that describes the specimen state and as soon as a decrease of the displacement is observed, the decrease of the computed load is kept vertically until an intersection with the original curve is observed. The envelope curve then follows the original curve until the new decrease of the displacement is met again and the procedure is repeated. The zoom-in figure in Figure 6a shows the procedure. Note that envelope curves were also proposed in [30,31]. However, when using the envelope curve alone, some essential information may be lost such as a possible snap-back. Also the area under the envelope curve is overestimated. So, the values of G f and W s should not be taken from the corresponding areas under the envelopes of the forcedisplacement and stress-strain curves. Instead, G f and W s are directly computed from the stored elastic energies of the broken elements by Equations (7) and (8).
In all tests presented herein, the characteristic length 0 c intrinsic to the lattice mesh is determined by performing the LD and DD tests on several mesh configurations without inclusions. The intrinsic effective width of the FPZ 0 FPZ is also determined via direct measurements. These values are shown in the plots of lengths as if the inclusion diameter or the surface fraction is zero. It is noted that these fits are calculated only from the mean values of FPZ in the cases of inclusions are introduced, so the value of FPZ of the homogeneous material is not taken into consideration in the fits. Also, the displayed fitted lines do not necessary mean that an affine relation is enlightened. It must rather be seen as a tendency since the data presents significant variations. As a consequence, the intersection of the fitted line with the vertical axis has no particular meaning i.e., one could also say the fit is only valid between the limits studied.  Numerical result Envelope (c) (d) Figure 6: The force-displacement curve (a) and the corresponding crack pattern (b) of the localized tensile test on the notched specimen; the stress-strain curve (c) and the crack pattern (d) of the distributed tensile test on the unnotched specimen.
The main observation from the Figure 7 is that when the inclusions are introduced, they have a strong effect on the FPZ size in both two-and three-phase materials. First, the mean values of FPZ in the case of heterogeneous material are greater then the value of FPZ in the case of homogeneous one. Second, in the case of inclusions are introduced, the fitted slope of the mean values of FPZ of the three-phase material is greater than that of the two-phase material. This means that when the ITZ is taken into account, the inclusion size plays a stronger effect on the variation of FPZ than the case in which the ITZ is not taken into consideration. So, according to our model, the internal length does not only depend on the size of the inclusions but also their constituents and therefore the presence of ITZ. The second observation is probably explained by the increase of the ITZ fraction when increasing the inclusion size of the threephase material, see Figure 8. Here, the ITZ plays a role of attractive zones for the crack propagation because of their lower strengths and stiffnesses. Accordingly, the greater fraction of the ITZ results in the larger mean value of FPZ compared to the mean value of FPZ of the two-phase material (in which the ITZ fraction is zero). In the case of d = 4 mm, the mean value of FPZ of the three-phase material does not differ from that of the two-phase one. This is related to the fact that the matrix prevails in the mesostructure in the case of d = 4 mm, as shown in Figure 8, and thus few inclusions are found on the crack path.  Furthermore, bigger values of the standard deviations are observed when increasing the particle size of the three-phase material as well as of the two-phase one even though this is less clearly observed in the twophase material than the three-phase one. We believe that increasing the inclusion size results in the increase of the inclusion fraction and as a result, the spatial distribution of inclusions plays a stronger role in the resulting value of FPZ . In two-phase material, by comparison between d = 8 mm and d = 10 mm, the standard deviation of FPZ does not significantly change. This is due to the fact that from d = 8 mm, the inclusions get dense in the mesostructure of the material, in that a change of the position of the inclusions does not have a strong effect on the value of FPZ . However, in three-phase material, the spatial distribution of inclusions still make sense on the variation of FPZ even though the inclusions get dense. This is reflected by the greater value of the standard deviation of FPZ in the case of d = 10 mm than that in the case of d = 8 mm, of the three-phase material. There is no doubt that this is due to the effect of the ITZ.

Path (b). A second series of tests is performed by following the variation path (b), that is with fixed surface fraction of inclusions and varied inclusion diam-
eter. The reason of this choice relies on the fact that the fundamental role of inclusion size d must be checked while keeping other parameters unchanged to suppress their possible effect. Figure 9 shows the plot of the mean value of the FPZ size FPZ with respect to the size d of inclusions of the path (b) of variation. For the sake of comparison, the same plot of the path (a) is shown as well. Surprisingly, it exhibits that the mean value of the FPZ size does not depend on the inclusion size for the path (b) of variation. It means that the FPZ size developed in this type of model material (brittle elastic) may not always be related to the inclusion size itself as usually observed in the literature. The observation is in agreement with that of [32], in which the width of the FPZ was experimentally measured on the surface of concrete specimens using a Digital Image Correlation (DIC) technique. However, it is in contrast to the results of [5,33] for concrete material, in which the experiments were carried out with X-rays and three-dimensional Acoustic Emission techniques leading to the conclusion that the width of the FPZ increases with the increase of the maximum inclusion size. But for the path (a) of variation, as previously shown, it is seen that the mean values of FPZ increases with the increase of inclusion size d that also results in the increase of the inclusion surface fraction. Figure 10 shows the crack patterns and the corre- sponding value of FPZ obtained when changing the inclusion size d within the variation path (b). The figure shows for each inclusion size d only one random distribution of inclusions in the specimen. It seems that the position of inclusions around the notches have an essential role on the resulting FPZ size. In fact, the crack is always initiated at the weak ITZ between the inclusion and matrix phases. With regard to the position of a notch -that can also be seen as a "weak link" -the crack is then propagated via the development of microcracks and at the end the macrocrack is formed by connecting the notch(es) and the broken elements (mainly in the ITZs). However, sometimes an inclusion is found just in front of the notch(es) and it plays a role of an obstacle that prevents the rupture of elements in the vicinity of the notch(es) and consequently, prevents the macrocrack to reach the notch(es). In this case, the macrocrack is finally formed by mainly connecting the broken ITZ elements. Therefore, the spatial distribution (positions) of the inclusions actually have an important role on the FPZ size in conjunction with the size of the inclusions. Nevertheless, for the case where the reference surface fraction of inclusion is kept almost constant (path (b)), the spacing between the inclusions seems to be constant regardless of the size of the particles, and thus the spatial distribution of the particles prevails more and more on their size in the resulting FPZ size FPZ . Actually, as shown in Figure 11 in which four different sets of inclusion positions with the diameter being 6 mm, the values of FPZ are finally different depending on the spatial distribution of the inclusions with regard to the notch position. This explains why changing the size of the in- clusions according to the path (b) does not change the value of FPZ averaged over five random spatial distributions of inclusions. On the contrary, within the path (a) of variation, changing the inclusion size leads to a change in the inclusion surface fraction together with the spacing between inclusions, and the FPZ size is affected by not only the size of inclusions but also the other structuring parameters (position and surface fraction of inclusions in our case). Still, the smaller the inclusion particle size, the larger the spacing between the particles in the path (a). This leads to the weaker influence of the spatial distribution of inclusions observed on the FPZ size. It is revealed in Figure 9 by the value of the standard deviation that is increased with the inclusion size.
Path (c). A third series of tests were performed by following the variation path (c) in which the inclusion surface fraction is varied while keeping the inclusion size constant at 8 mm in order to evaluate the only influence of the inclusion surface fraction (or equivalently the inclusion spacing, since the latter is inversely proportional to the former) on the FPZ size FPZ . Figure 12 shows the variation of the mean value of FPZ with respect to the inclusion surface fraction P a . For the sake of comparison, the results of above studies for the variation paths (a) and (b) are plotted as well, but in the ( FPZ -P a ) space. The main observation is that the mean value of FPZ of the path (a) and the path (c) does increase with the increase of the inclusion surface fraction P a , whereas that of the path (b) does not change. This is simply explained by the fact that the spacing between inclusions decreases with the increase of the inclusion surface fraction within the variation paths (a) and (c), whereas it seems to be "constant" (or hardly changed) within the path (b). By comparison the path (c) with the path (a), it is observed, however, that the increase rate of FPZ with respect to P a , which is represented by the fitted slope, of the path (c) is smaller than that of the path (a). A suitable explanation for this observation is that within the path (a), the size of and the spacing between the inclusions do change (increase and decrease, respectively) at the same time with respect to the increase of P a whereas only the spacing of the inclusions does decrease with respect to the increase of P a within the path (c). Therefore, the observation could lead to the evidence that the FPZ size depends on both the inclusion spacing (which is just a consequence of the inclusion surface fraction) and inclusion size.
Larger specimen width compared to inclusion size. The above studies of FPZ are performed on slender specimens (small ligament size, i.e. in the order of 3 × d compared to inclusion size d). These slender specimens were used to ensure the damage distribution in DD-tests as homogeneous as possible for studies of the characteristic length c which is presented in Section 4.2. However, it was shown that the FPZ size FPZ performed on these slender specimens strongly depends on the position of inclusions. This can be considered as a drawback for an attempt to correlate the FPZ size and the inclusion size. Therefore, it would be better if the study of FPZ is performed on a larger specimen size compared to the inclusion size. To this end, tensile tests are performed on the specimen shown in Figure 13 with mono-sized inclusion structures and the inclusion size is varied by taking the value 4, 6, 8 and 10 mm. Three variation paths above are also considered here.  Figure 16 shows the variation of FPZ with respect to the surface fraction P a when keeping the inclusion size constant at 8 mm (path (c)). It exhibits that the variation of FPZ with respect to the position of inclusions is less important than the previous cases. Also, better fit is obtained with higher coefficients of correlation (0.99, 0.93 and 0.96 for paths (a), (b) and (c) respectively). Figures 14 and 15, it can be seen that, in contrast to the results tested on the slender specimens, a higher influence of the inclusion size d on the FPZ size FPZ of the variation path (b) compared to that of the variation path (a). Indeed, a higher value of the fit slope is obtained within the variation path (b). It is likely due to the fact that when analyzing on the larger specimen (compared to the inclusion size), the sensitiv- ity of FPZ with respect to the position of inclusions is less important than testing on the slender specimen, and thus the role of the inclusion size in the FPZ size prevails over the position. So, in conjunction with the influence of the inclusion surface fraction on the FPZ size (which can be observed in Figure 16), the higher influence of the inclusion size on the FPZ size is obtained within the path (b) than that within the path (a) because varying the inclusion size in the path (b) is combined with a higher surface fraction of inclusions than in the path (a). Therefore, a partial conclusion is that depending on the relative size between the macroscopic size (specimen size) and the mesoscopic size (inclusion size), the influence of the mesostructure on the FPZ size is different. When this relative size is small, the effect of the position of inclusions or the spacing between inclusions of the mesostructure prevails over the effect of the inclusion surface fraction. When the relative size is more important (≥ 5), an inverse influence is observed. In any case, the influence of the inclusion size on the FPZ size is always recognized.

Ligament size
The specimen geometry is shown in Figure 17 and the dimensions of the specimens used in the numerical test are given in Table 2, with the same size but with different notch lengths that results in different ligament lengths: 90, 80, 65 and 50 mm. They are labeled by L, M, S and XS respectively for convenience. Figure 17: Specimen geometry. In order to study the influence of the maximum inclusion size d max on the FPZ size, the tensile tests are performed on the specimens made of the material with a polydisperse inclusion structure with d max being 6.3, 8, 10, 12.5 and 16 mm and the reference inclusion surface fraction is kept constant at 45%. The minimum inclusion size d min is 3.15 mm. All the inclusion gradings are generated by the Fuller's curve which is an "ideal" grading curve [34]. In the study, for each inclusion grading up to d max , five inclusion structure realizations are generated with independently random distribution of inclusion positions. The specimens are loaded in tension by directly imposing the vertical displacement increment on the nodes of the top boundary of the specimens while vertically fixing the nodes of their bottom boundary. Figure 18 shows the relationship between the FPZ size FPZ with respect to the maximum inclusion size d max for the specimens corresponding to four ligament lengths L, M, S and XS. This figure also shows the FPZ size of L, M, S, XS specimens in which no inclusion structure has been introduced. It is seen that, when the inclusion structures are introduced, it always results in a larger FPZ size than the one computed with the homogeneous cases. For a given value of d max , the mean value of FPZ is systematically increased when the ligament size is increased. In addition, the increase rate of FPZ is also increased with d max . It results in an higher slope of variation of FPZ as a function of d max for a larger ligament size. It is also observed that the increase rate of the slope of variation of FPZ decreases with the increase of the ligament size from XS specimens to L specimens. So, a stabilized value of variation slope can be achieved for specimens with the ligament size being in order of specimen width. When the ligament size is half (and may be lower by extrapolation) specimen width (the XS specimens), the variation slope is negligible, which means that the inclusion size appears to have no influence on the mean value of FPZ . It may suggest that the FPZ has not enough time to develop completely within the specimens with "too short" ligament length. Between these limits, the slope variation evolves progressively, indicating that both the inclusion structure and the specimen dimension itself can play a role on the FPZ size. The maximum-inclusion-size independence of FPZ for the specimens with too short ligament length is in agreement with the previous study performed on the specimen which also has a short ligament length.  Figure 19 shows the crack patterns (selected among several realizations of inclusion positions) and the value of FPZ corresponding to the smallest inclusion sizes (d max = 6.3 mm) and the biggest ones (d max = 16 mm) for the two extreme ligament lengths (XS and L). In the case of XS specimens (Figures 19a and 19c), whatever the maximum inclusion size, a crack without bifurcation crosses the ligament by connecting ITZ elements with a path that seems to be the shortest. Whereas, in the case of L specimens (Figures 19b and 19d), even if only one crack finally crosses the ligament, a number of microcracks occur either side of the inclusions. As a consequence, the FPZ size is in direct proportion with the maximum inclusion size in the latter case. This should be simply because the microcracks have enough space to develop in specimens with large ligament. The black dots indicate broken elements with the largest opening -read as the macrocrack -while the blue dots stand for the remaining broken elements -read as the microcrack.

Material characteristic length versus FPZ size
The influence of the material mesostructure on the FPZ size has been studied. The aim is now to question whether the same influence can be observed on the characteristic length of the material. Although many simulations were carried out to answer this question, we only focus herein on two mesoscopic features that may influence the characteristic length c , the inclusion size with (i) fixed position and (ii) fixed surface fraction.
First, the lattice simulations are performed by varying the inclusion size while both the positions and the number of inclusions remain unchanged. This concerns the path (a) of mesostructure variation in which the monodisperse diameter of inclusions is changed by setting their values to 4, 6, 8, and then 10 mm. Figure 20 shows the relation between the characteristic length c and the inclusion size d. For a comparison with the FPZ size FPZ , the relation between FPZ , computed from the LD tests, and the inclusion size is shown as well. The main observation is that c and FPZ have the same order of magnitude and trend with respect to the inclusion size.
The increase in d results in an increase of standard deviations of c as previously observed in the variation of the FPZ size FPZ . This is explained by the same reasons mentioned above for the FPZ size. Now the question is whether the variation of the characteristic length c with respect to the inclusion size d still follows the variation of the FPZ size FPZ with respect to d if we only do vary the size d of inclusions while keeping the surface fraction of inclusions as constant as possible? For this end, the path (b) of the mesostructure variation is used to study the influence of the inclusion size d on the characteristic length c , in which the "reference" surface fraction of inclusions is kept at 45% when changing the inclusion size. Figure 21 shows the characteristic length of the material c as a function of the inclusion size d. The plot between the FPZ size FPZ and d is shown as well. It exhibits that increasing d does not lead to an increase of c , as previously observed in the case of FPZ . With a fixed value of the inclusion size, the resulting characteristic length of the material varies upon the spatial distribution of inclusions. However, the mean value of the characteristic length with respect to the spatial distribution of inclusions seems to be unchanged upon the increase of the inclusion size. The reason for this non-sensitivity may be related to the fact that the spacing between the inclusions, thus the spacing between the ITZs, seems to be insignificantly changed when the inclusion size is increased, as previously shown for the case of the FPZ size.

Conclusion
Two types of tensile tests were performed to study the key features that influence the FPZ size FPZ and the material characteristic length c . The assessment of FPZ is achieved via localized damage (LD) tests while c is measured via both LD tests and distributed damage (DD) tests. The numerical simulations are performed on the brittle elastic model material with inclusions. The material is then modeled as the three-phase material with the inclusion and matrix phases and the interfacial transition zone (ITZ) in-between them. Not only the mesostructure of the material but also the specimen geometry and the ligament size are varied in order to analyze their effect on the resulting FPZ size and material characteristic length. Five independent realizations of inclusion positions are generated for each case of the mesostructure so that the average values of FPZ and of c over that five realizations are used to analyze the effect of the mesostructure. The study points out the influences of: the inclusion size with fixed surface fraction, the inclusion size with the number and the position of inclusions unchanged, the inclusion surface fraction with fixed size, and the ligament size of the specimen on FPZ and c .
From that extensive study, the following conclusions can be drawn: • It appears that not basically the size, but other parameters that characterize the inclusion structure of the material such as the surface fraction, the fabric, the connectivity. . . strongly affects the size of the FPZ, and thus the characteristic length of the material.
• The measured value of the FPZ size is also dependent on the specimen geometry and the ligament size of specimens. Therefore, it is difficult to avoid the conclusion that the FPZ size is not an intrinsic property of the material as usually believed. However, it seems true that the FPZ size remains in the same order when the tested system is the same (mesostructure, global geometry and dimensions, loading conditions. . . ).
• The assessment of the characteristic length of the material is essential for using its value as the internal length in nonlocal models. However, just like the FPZ size, it is difficult to avoid structural effects in the method of measurement of the characteristic length.
This is a first step to study the influence of inclusion properties on the characteristic length and several interesting qualitative conclusions on the numerical model material have already been pointed out. For the future work, we plan to study the effect of the mechanical properties, especially the ratios of the different stiffnesses and strengths of the material phases, on the resulting FPZ size and material characteristic length. Future developments will also aim to develop the numerical model to be more representative of quasi-brittle materials, especially concrete.