W Leclerc
email: willy.leclerc@u-picardie.fr

A Ammar

D Moukadiri

T Dridi

M Guessasma

Halo approach to model cracks initiation and propagation in 3D Discrete Element Method simulation of homogeneous and heterogeneous materials

Keywords: Discrete Element Method, Numerical simulation, Crack, Composite materials, Stress field

The present contribution deals with an investigation of the halo approach to simulate cracks initiation and propagation in homogeneous and heterogeneous materials by Discrete Element Method using the hybrid lattice-particle model. Recent developments highlighted that the halo, which can be described as a spherical interaction domain of controlled size, enables to efficiently reduce the intrinsic dispersion of the stress in discrete modelling. This paper aims at studying the benefits of such a concept to provide suitable cracks initiation and propagation through two applications: the diametral compression test of a cylindrical tablet composed of a brittle marble and the failure process of a flax fibre/bio-based composite material. For that purpose, based on an original process ensuring the connexity within each halo during cracks propagation, numerical simulations are led to determine the influence of the halo size on the validity of stress field before damaging, and ensuing failure process with respect to theoretical expectations and finite element data. Results exhibit that a minimum halo diameter about 15 times the particles size is required to obtain proper stress fields and crack patterns, but care must be taken to avoid too large halos which could lead to excessive stress smoothing and unrealistic cracks.

Introduction

Composite materials are a promising class of engineering materials enabling a wide range of industrial applications due to their multiplicity of shape, size and properties. The main idea consists in combining two or more components to produce a material with superior properties than their individual constituents. Thus, composite materials are widely used in various industries as automotive, aerospace or civil engineering sectors to take benefit of their advantageous characteristics in terms of weight, toughness, acoustic behaviour, cost, ... Nevertheless, composite structures require a continuous optimisation to respond to current economic and environmental issues. As a result, efforts are placed on manufacturing processes, but also on numerical simulation to extend their lifetime. Indeed, due to severe operating conditions including heavy loads and high temperatures, the durability of composite materials remains poorly controlled while it is involved in major security problems. In such a context, the main scientific challenge lies in setting up reliable and robust numerical approaches to predict the mechanical behaviour and the failure of composites by reproducing real conditions of use. In the literature, several numerical methods give the possibility to simulate cracks onset and their propagation in homogeneous and heterogeneous media. Most of them are based on a mesh discretisation as the Finite Element Method (FEM) which is very flexible and enables to model cracks in several manners. Thus, a first methodology consists in introducing pre-existing cracks path between elements using the Virtual Crack Closure Technique (VCCT) [START_REF] Krüger | The virtual crack closure technique : History, approach and applications[END_REF] or the Cohesive Zone Model (CZM) [START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF][START_REF] Zhang | Cohesive zone modeling of dynamic failure in homogeneous and functionally graded materials[END_REF]. However, such approaches exhibit some limitations related to the mesh updating during the failure process. That is why, a second class of methods taking benefit of an enriched formulation can be privileged. The most famous example is the eXtended Finite Element Method (XFEM) [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF] which is well-adapted to describe simple cracks onset and propagation without bifurcation and nucleation. Moreover, other techniques designated as meshless as the cracking particle method [START_REF] Rabczuk | A three-dimensional large deformation meshfree method for arbitrary evolving cracks[END_REF] does not lie on a mesh discretisation or topology, and are de facto more adapted to treat simultaneous cracks opening and closure arise. This type of methods is clearly advantageous to reproduce complex crack phenomena but turns out to be also more difficult to handle.

Distinct approaches and more particularly the Discrete Element Method (DEM) constitute a third class of methods based on a discrete representation of a material. Typically, the DEM is used to treat physical systems that are difficult to handle using a classical continuous approach. Initially introduced by Cundall and Strack [START_REF] Cundall | Discrete numerical model for granular assemblies[END_REF] in the field of rock mechanics, this was later used for different applications in various domains such as the granular flow via the study of multi-contact frictional problems as the silo discharge [START_REF] Nicot | On the definition of the stress tensor in granular media[END_REF], the tribology [START_REF] Cao | Coupling continuous and discontinuous descriptions to model first body deformation in third body flows[END_REF][START_REF] Haddad | A DEM-FEM coupling based approach simulating thermomechanical behaviour of frictional bodies with interface layer[END_REF], the mechanics of soils [START_REF] Wang | Effect of soil particle size on soil-subsoiler interactions using the discrete element method simulations[END_REF][START_REF] Sun | A p-y model for predicting the lateral nonlinear interaction between pile and soil-rock mixture material based on discrete element modeling[END_REF], electrical [START_REF] Bourbatache | DEM ball bearing model and defect diagnosis by electrical measurement[END_REF] and thermal [START_REF] Nguyen | Discrete modeling of granular flow with thermal transfer : Application to the discharge of silos[END_REF] transfer processes in granular media, fluid-structure interaction phenomena as fluidisation [START_REF] Arkawazi | Modeling the hydrodynamic forces between fluid-granular medium by coupling DEM-CFD[END_REF] or even additive manufacturing [START_REF] Steuben | Discrete element modeling of particle-based additive manufacturing processes[END_REF][START_REF] Xin | Discrete element simulations of powder-bed sintering-based additive manufacturing[END_REF]. Moreover, other studies have been carried out to model continuous media using DEM. In this context, Potyondy and Cundall [START_REF] Potyondy | A bonded-particle model for rock[END_REF] introduced the Bonded Particle Model (BPM) to simulate a rock via a dense particulate system composed of polydisperse spherical particles that are bonded together at their contact points by spring elements. Furthermore, based on the works of Schlangen et al. [START_REF] Schlangen | Fracture simulations of concrete using lattice models : Computational aspects[END_REF], André et al. [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF] developed a hybrid lattice-particle model in which a set of spherical particles are bonded using cohesive Euler-Bernoulli beam elements. Leclerc et al. [START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF] later investigated a similar paradigm which is distinguishable by coupling terms between bending and tangential terms. Such an approach has been applied to simulate the mechanical [START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF][START_REF] Haddad | Application of DEM to predict the elastic behaviour of particulate composite materials[END_REF] and thermo-mechanical [START_REF] André | A discrete element thermo-mechanical modelling of diffuse damage induced by thermal expansion mismatch of two-phase materials[END_REF][START_REF] Leclerc | On a discrete element method to simulate thermal-induced damage in 2D composite materials[END_REF] behaviour of composite materials. Besides, André et al. [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF] introduced the Removed Discrete Element Failure (RDEF) process which is based on the deletion of a set of particles when a given failure criterion is reached. In the case of indentation test, the authors exhibited that such a concept provides more realistic crack patterns than the classical Rankine criterion which lies on the bond deletion. This paradigm was later exploited to simulate cracks onset and propagation in homogeneous [START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF] and heterogeneous materials [START_REF] Le | Modeling damages and cracks growth in composite with a 3D discrete element method[END_REF][START_REF] Leclerc | On the suitability of a Discrete Element Method to simulate cracks initiation and propagation in heterogeneous media[END_REF][START_REF] Alhajj-Hassan | On the suitability of a 3D discrete element method to model the composite damage induced by thermal expansion mismatch[END_REF]. In addition, other approaches based on CZM have also been set up to reproduce interfacial debonding in composite materials using DEM [START_REF] Leclerc | On a discrete element method to simulate thermal-induced damage in 2D composite materials[END_REF][START_REF] Le | Modeling damages and cracks growth in composite with a 3D discrete element method[END_REF].

In parallel to these contributions, a great effort has been made to better define the notion of stress tensor in homogeneous and heterogeneous media described by a particulate system. Thus, Jebahi et al. [START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF][START_REF] Jebahi | Virial stress-based model to simulate the silica glass densification with the discrete element method[END_REF] proposed a symmetrised expression of the stress tensor given by the classical Zhou's formulation [START_REF] Zhou | A new look at the atomic level virial stress : on continuum-molecular system equivalence[END_REF] which takes into account the neighbourhood region of the investigated particle. Besides, based on specific tensile tests, Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF] studied the level of dispersion of the stress obtained at the scale of the particle which is always heterogeneous, even if it is theoretically homogeneous. Their results exhibited that this strongly depends on microscopic parameters describing the cohesive beam element and can reach a Coefficient of Variation (CoV), i.e. a ratio between the standard deviation and the mean, higher than 100% in specific cases. Furthermore, the same authors developed an approach called halo to relate the level of stress dispersion obtained using the symmetrised expression of Zhou's formulation to the span of the region of interest which is described by a sphere, qualified as "halo", of controlled size and same center than the investigated Discrete Element (DE). Moreover, several authors [START_REF] Kumar | Effect of packing characteristics on the discrete element simulation of elasticity and buckling[END_REF][START_REF] Leclerc | Effect of packing characteristics on the simulation of elasticity and brittle fracture by the Cohesive Discrete Element Method[END_REF] highlighted the great sensibility of BPM and cohesive beam model to intrinsic features of particulate system as the volume fraction of particles and the mean cardinal number. It turns out that cracks initiation and propagation are not only dependent on the packing characteristics as verified by Leclerc [33] but potentially also to the level of stress dispersion. Indeed, the high CoV of stress field induces local stress peaks which could lead to premature cracks onset and an improper failure process. Such effects could be also emphasised by the characteristics of particulate system. That is why, the main objective of the present contribution is to evaluate the influence of level of stress dispersion on the failure process by taking benefit of the halo approach. In order to limit undesirable effects due to packing characteristics, the cardinal number, i.e. the average number of particles in contact with a given DE, is set to 7.5 according to Leclerc's conclusions. In a first step, since Moukadiri et al. limited their investigations to the case of a lower cardinal number set to 6.2, we aim at evaluating the level of stress dispersion in the targeted configuration and determining the influence of the halo size as well as the cardinal number on the CoV of stress field. Based on the obtained results, a diametral compression test is then simulated in the context of a brittle marble and the failure process is investigated for several halo sizes to highlight the influence of the level of stress dispersion on cracks onset and propagation. Finally, in a third step, the case of a flax fibre/bio-based composite material is studied. Two models are investigated, namely a unidirectional composite composed of a single fibre strand and a representative volume element constituted of 8 independent aligned strand fibres generated from digi-tal images. Quasi-static transverse tensile tests are led and the failure process is examined as function of the halo size.

Following the present introduction, the paper is outlined as follows. Section 2 describes the hybrid lattice-particle model used to simulate a continuous medium by a twofold lattice and particulate description. Section 3 is dedicated to the determination of the level of stress dispersion in light of the halo approach. Section 4 discusses the influence of the halo size on stress fields and failure process during a diametral compression test. Finally, section 5 focuses on the case of a flax fibre/bio-based composite material for which single and multi-fibre configurations are examined in the context of quasi-static transverse tensile tests.

Hybrid lattice-particle model

The present section is dedicated to a brief introduction to the hybrid lattice-particle model which allows for simulating a continuous medium using a mixed lattice and particle description. In such a paradigm, a given continuous medium is modelled by a particulate system composed of rigid spheres in point contact in which all binary contacts are replaced by cohesive elements with targeted properties.

Particulate system

Particulate systems composed of spherical particles in point contact are randomly generated using a specific enlargingdisplacement process described in a previous contribution [START_REF] Leclerc | Effect of packing characteristics on the simulation of elasticity and brittle fracture by the Cohesive Discrete Element Method[END_REF]. In order to properly represent a continuous medium with given elastic properties and behaviour, a great care has to be taken to accurately control several intrinsic characteristics, namely the volume fraction of particles φ, the cardinal number Z, and the isotropy of the particulate system. In the present context of random packings, a maximum volume fraction is targeted according to the paradigm of Random Close Packing (RCP). Experiments and numerical simulations exhibited that an RCP is characterised by a volume fraction φ=0.64 which signifies that about 64% of a given container is occupied by spheres [START_REF] Donev | Improving the density of jammed disordered packings using ellipsoids[END_REF]. The cardinal number Z is often set to 6.2 according to previous works [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF][START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF]. However, as highlighted by Leclerc [33] in the context of indentation tests, this value leads to undesirable crack initiation and propagation due to very low local cardinal numbers (less than 3). As a result, as recommended by the author, a higher value of 7.5 is here handled to avoid such effects. Nevertheless, one has to keep in mind that the generation process leads to a Z lower than 7.5. Consequently, a sensitivity parameter ε is introduced to increase the range of admissible contacts and include binary interactions between two non-contacting particles as depicted on Figure 1. Thus, let l i j be the center to center interparticle distance between two particles i and j of radii R i and R j respectively, interaction occurs under the following condition : The isotropy of the particulate system is a crucial point which directly affects the representativeness of the continuous medium. This is only verified using random systems with a sufficient density of particles. Typically, as exhibited in the literature [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF], a minimum of 10,000 particles is required to ensure a quasi-randomised distribution of contact angles. However, such a result greatly depends of the tolerated level of geometric anisotropy on one hand, and could be largely underestimated in another hand to ensure isotropy from a mechanical standpoint as suggested by Leclerc [START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF]. Furthermore, there is a consensus that monodisperse particles lead to undesirable directional effects. A solution consists in introducing a small polydispersity controlled by a small CoV describing the variability of particle dimensions, defined as the ratio of the standard deviation to the mean particle radius. In the present work, this value is set to 0.3 according to previous studies [START_REF] Leclerc | Effect of packing characteristics on the simulation of elasticity and brittle fracture by the Cohesive Discrete Element Method[END_REF][START_REF] Leclerc | An efficient numerical model for investigating the effects of anisotropy on the effective thermal conductivity of alumina/Al composites[END_REF].

l i j ≤ R i + R j (1 + ε) (1)

Cohesive beam model

We consider a cohesive beam model based on Euler-Bernoulli formulation [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF][START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF] which proved to provide more realistic crack patterns than the BPM [START_REF] Cundall | Discrete numerical model for granular assemblies[END_REF][START_REF] Potyondy | A bonded-particle model for rock[END_REF] for specific tests as the indentation test [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF]. In such an approach, each beam element is described by a set of elastic coefficients, the Young's modulus E µ and the shear modulus G µ and geometrical parameters, namely the beam length L µ which corresponds to center to center interparticle distance l i j , a circular cross-section S µ of radius a µ and axial, I µ,t and I µ,b , and polar J µ quadratic moments (Figure 2). However, in the present context, it is more convenient to introduce a sole planar quadratic moment I µ so that

I µ =I µ,t =I µ,b =J µ /2.
Besides, for simplification purposes, we introduce a coefficient r µ > 0 which can be described as a dimensionless beam radius and reads:

r µ = 2a µ R i + R j ( 2
)
where R i and R j are the radii of both particles i and j in contact. Thus, S µ and I µ can be written as function of r µ as follows : As a result, knowing that L µ = (R i + R j ) and G µ has a little impact on the macroscopic behaviour, the beam model can be reduced to two parameters E µ and r µ if we consider dense enough particulate systems to avoid particle size effects. From a practical standpoint, cohesion is introduced by attractive internal forces described by a six-component vector of generalised forces. This is composed of normal and tangential components allowing to counteract the relative displacement, and moment ones allowing to resist to bending and torsion effects. In fact, the vector of generalised forces is derived from a 6-by-6 system the components of which depend on E µ and r µ parameters and the expression of which is given in previous works [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF][START_REF] Leclerc | Effect of packing characteristics on the simulation of elasticity and brittle fracture by the Cohesive Discrete Element Method[END_REF]. An expression of Newton's second law of motion, which derives from internal cohesive forces, is solved at each particle center by considering the contributions of each neighbouring DE. The numerical resolution is carried out using an explicit time integration based on the velocity Verlet scheme which is well-suited to large-scale calculations. Please keep in mind that a small Rayleigh damping is also added to internal forces to reduce dynamic effects over the time.

S µ = πa 2 µ = π r 2 µ 4 R i + R j 2 (3)
I µ = π a 4 µ 4 = π r 4 µ 64 R i + R j 4 (4)

Macroscopic elastic properties

Microscopic parameters E µ and r µ have to be accurately determined to target macroscopic elastic properties of the modelled continuous medium. Typically, the micro-macro relationship between E µ and r µ and a pair of independent elastic coefficients as the macroscopic Young's modulus E M and the Poisson's ratio ν M have to be determined using a calibration procedure. For that purpose, cubic discrete domains are generated and submitted to specific tensile tests as described in [START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF]. Macroscopic elastic coefficients (E M and ν M ) are extracted from numerical simulations and related to input microscopic parameters (E µ and r µ ). This process can be repeated for a large range of configurations so that it is possible to extract relationships taking the following expressions [START_REF] Jebahi | Virial stress-based model to simulate the silica glass densification with the discrete element method[END_REF]:

r µ = f (ν M ) (5)
E µ = g(r µ )E M ( 6
)
where f and g are two unknown functions which can be approximated by polynomial regression. However, one must keep in mind that f and g greatly depend on the characteristics of the particulate system Z and φ so that the micromacro relationship is only valid for given coordination and volume fraction. Besides, such a relation is also affected by the density of particles. That is why, dense enough particulate systems must be handled to minimise the influence of the number of particles on the macroscopic elastic properties. The required minimum number of particles n de varies in the literature from 10,000 [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF] to 700,000 [START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF] spherical particles without any real consensus. In fact, it turns out that the convergence is potentially sensitive to several parameters among them we can identify the characteristics of the particulate system (Z,φ), r µ , the domain geometry, the test configuration and the tolerated level of convergence which is quite subjective. In the present work, as discussed in subsection 2.1, packing characteristics Z and φ are respectively set to 7.5 and 0.64. Thus, for given domain geometry and configuration, the convergence could be mainly affected by r µ . Figures 3a and 3b GPa and ν M = 0.18 which correspond to microscopic parameters E µ = 68.9 GPa and r µ = 0.6. Thus, based on the fact that E µ has strictly no effect on macroscopic elastic properties due to f and g dependencies, numerical outputs provide a direct comparison of convergence curves obtained with a low (r µ = 0.1) and a high (r µ = 0.6) dimensionless radius. Please notice that each data dot represents the average value given by a set of 4 random particulate systems of same density. Results exhibit a two-step convergence in both cases. For a density less than about 200,000 particles, curves are characterised by a fast decrease/increase of properties. Conversely, for a density higher than about 200,000 DE, results evolve very slowly. Thus, r µ has a poor influence on convergence speed and n de can be reasonably set to 200,000 particles in the present context. Nevertheless, a higher density is required to get a strict convergence which is probably close to the number of 700,000 particles proposed by Leclerc in a previous contribution [START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF].

Stress field characterisation

Recent works [START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF][START_REF] Jebahi | Virial stress-based model to simulate the silica glass densification with the discrete element method[END_REF] proposed a symmetrised expression of the equivalent stress tensor associated to a given DE i:

σi = 1 2Ω i Z i j=1 1 2 f i j ⊗ l i j + f i j ⊗ l i j (7)
where Ω i is the effective volume of DE i which can be defined as the ratio between the volume of the particle V i and the volume fraction φ of the particulate system to ensure the volume conservation of the continuous domain [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF]. Z i is the local cardinal number associated to particle i, f i j is the interparticle force vector exerted by DE j on DE i, l i j is the relative position vector between particles i and j. ⊗ refers to tensor product. Typically, a stress field determined using Equation 7is quite heterogeneous, even if this is theoretically homogeneous, due to the limited number of contacts associated to a given particle. Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF] showed that the stress follows Gaussian distribution with a given standard deviation and a mean equal to the expected macroscopic value. Thus, in order to characterise this dispersion, they considered a coefficient of variation associated to the stress, CoV σ , which is defined as the ratio between the standard deviation and the mean of a given stress component σ. Based on a specific tensile test for which σ is theoretically homogeneous, the authors set up a methodology to determine CoV σ which potentially depends on geometrical parameters (Z and φ) and r µ . Figure 4 illustrates the evolution of Cov σ as function of r µ for Z = 6.2 and Z = 7.5, and φ set to 0.64. Please notice that all results are obtained using a homogeneous cubic domain composed of at least 400,000 particles which is dense enough to properly describe a continuous medium without any effect related to particle size according to the discussion of subsection 2.3. Furthermore, σ is set to be the maximal principal stress. Numerical outputs exhibit that a higher cardinal number leads to a lower stress dispersion whatever r µ . Such a result is consistent with Equation 7which defines the stress as an average quantity over the local cardinal number. Besides, CoV σ is greatly sensitive to r µ which signifies that the material also affects the level of stress dispersion. Thus, CoV σ respectively decreases from 120% to 43% and from 67% to 29% for Z = 6.2 and Z = 7.5 as r µ increases from 0.1 to 0.9.

Halo approach to control stress field dispersion

The present section describes the halo approach recently introduced by Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF] to better control the stress field dispersion in DEM simulation of a continuous medium. Numerical investigations are led in the context of a homogeneous cubic discrete domain. Results exhibit the influence of geometrical and microscopic parameters on its ability to reduce the level of stress dispersion characterised by CoV σ .

Halo approach

Stress field dispersion inherent to discrete simulations affects the modelling and simulation of damage and failure of materials in which cracks initiation and propagation arise. Indeed, this potentially leads to an acceleration of these phenomena due to local undesirable stress peaks. That is why, the halo approach was proposed to better control the level of dispersion and its influence on the failure process. This latter is based on a symmetrised Zhou's formulation of stress tensor [START_REF] Zhou | A new look at the atomic level virial stress : on continuum-molecular system equivalence[END_REF]:

σ h I = 1 2V h I i⊂V h I Z i j=1 1 2 f i j ⊗ l i j + f i j ⊗ l i j (8)
where V h I is a volume surrounding a given particle I. Z i is the local cardinal number associated to a particle i, f i j is the interparticle force vector exerted by DE j on DE i, l i j is the relative position vector between particles i and j. In halo approach, a spherical geometry of radius R h I is associated to each particle I of radius R h I (R h I ≤ R h I ) belonging to particulate system so as, in the most general case where the particle I is located far from the domain boundaries,

V h I = 4 3 πR h I 3
(in the opposite case, only a part of the sphere could be considered). Besides, we assume that each particle center coincides with the center of its associated halo as depicted on Figure 5. The main paradigm of halo approach is to control the level of stress dispersion by playing on the halo size which is defined as the ratio between R h I and R h I . Thus, the determination of the relationship between the halo size and CoV σ enables a direct control of the dispersion. The difficulty lies on the fact that such a relation is a priori unknown and depends on the characteristics of the particulate system (Z,φ) and the material elastic properties via r µ parameter. Nevertheless, for a given configuration of particulate system with set values of Z and φ, this is only sensitive to r µ and can be easily determined by calibration process and formulated using polynomial regression.

Stress field dispersion

Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF] exhibited that the stress determined using Equation 8 follows a quasi-Gaussian distribution so as the level of stress dispersion can be also characterised by CoV σ in the context of halo approach. Based on the configuration discussed in subsection 2.4 in terms of discrete domain and density of particles and the methodology introduced by Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF], the relationship between the halo size and CoV σ is investigated for particulate systems characterised by Z = 7.5 and φ = 0.64. Figure 6a illustrates the influence of the halo size (R h I /R h I ) on CoV σ in the context of the bio-based epoxy resin (r µ =0.1). Two threshold levels of 5% and 3% are also depicted. Results exhibit that CoV σ varies from 67.9% to 2.47% as the halo size varies from 1 to 15. This highlights the ability of the halo approach to reduce the stress level in the present context. Please notice that the present curve can be reversely used to determine the suitable halo size to reach an expected stress dispersion. Thus, for example, a minimum halo size of 10 is required to reach a CoV σ of 5%. Figure 6b shows the evolution of CoV σ as function of r µ for several halo sizes in the 1-15 interval. Please notice that results obtained for a halo size equal to 1, i.e. a halo radius equal to DE radius, are equal to results depicted on Figure 4 in the context of the stress determined using Equation 7. This exhibits the consistency of the definition of halo size as regards to the stress tensor expression at the scale of the particle. Two main conclusions can be drawn. First, whatever r µ , the greater the halo size is the lower the stress dispersion is and the decrease is fast for low halo sizes. Thus, for r µ set to 0.6 (case of the brittle marble), CoV σ varies from 35,6% to 1.08% as the halo size varies from 1 to 15. Nonetheless, CoV σ is already reduced to 4.18% for a halo size of 5. Second, the greater the halo size is, the lower r µ impacts on the level of dispersion. Thus, for a halo size equal to 15, CoV σ varies from 2.21% to 1.28% as r µ varies from 0.1 to 0.9.

Diametral compression test

The present section deals with the diametral compression test also called Brazilian test which allows for indirectly determining the tensile strength of brittle materials [START_REF] Hondros | The evaluation of Poisson's ratio and modulus of materials of a low tensile resistance by the Brazilian (indirect tensile) test with particular reference to concrete[END_REF][START_REF] Fairhurst | On the validity of the "Brazilian" test for brittle materials[END_REF]. In such a method, a compressive load is applied to a cylindrical tablet so that tensile stresses arise at the center of the specimen and lead to mode I crack initiation and propagation along the compression line. This study aims at investigating the influence of the halo approach on the validity of the numerical failure process which is typically thwarted by stress concentrations located near the loading area [START_REF] Wang | The flattened Brazilian disc specimen used for testing elastic modulus, tensile strength and fracture toughness of brittle rocks: analytical and numerical results[END_REF].

Study configuration

We consider a cylindrical tablet submitted to a compression loading along a diametral line as depicted on Figure 7. Please notice that compressive forces are directly applied at contact points between the sample and the experimental device and not over a strip of the circumference as usually observed in the literature to artificially reduce local stress concentrations. Thus, numerical simulations are closer to theoretical configuration on which is based the Hertz solution. In present investigations, we consider a cylindrical sample composed of a brittle marble which was briefly discussed in subsection 2.3. This choice of material is mainly motivated by the desire to avoid prohibitive calculations in dynamic mode. Indeed, the brittle marble has a low Ultimate Tensile Strength (UTS), set to 5.9 MPa in present works, in comparison with other materials such as ceramics. As a reminder, the material is supposed as homogeneous, elastic and isotropic with a Young's modulus E M = 17.3 GPa and a Poisson's ratio ν M = 0.18 to which are associated microscopic parameters E µ = 68.9 GPa and r µ = 0.6. Furthermore, the friction between the device and the specimen is not considered and we assumed that the failure follows the Griffith criterion without any effect of the intermediate principal stress. As regards the cylindrical tablet dimensions, the radius is set to 20mm and the height is set to 16mm.

The compression loading is imposed in the x direction using a constant velocity V equal to 1mm/min. All numerical simulations are performed using particulate systems of characteristics Z = 7.5 and φ = 0.64 and composed of 200,000 particles according to the discussion of subsection 2.3.

(a) R h I /R h I = 1 (b) R h I /R h I = 5 (c) R h I /R h I = 10 (d) R h I /R h I = 15
/R h I = 1 (b) R h I /R h I = 5 (c) R h I /R h I = 10 (d) R h I /R h I = 15

Stress fields

In the present context of a homogeneous elastic material, Hertz solution provides theoretical formulations of σ xx and σ yy stress fields:

σ xx (x, y) = -2F πh       (R -x) 3 β 4 1 + (R + x) 3 β 4 2 - 1 2R       (9)
σ yy (x, y) = -2F πh       y 2 (R -x) β 4 1 + y 2 (R + x) β 4 2 - 1 2R       ( 10
)
where

β 2 1 = (R -x) 2 + y 2 and β 2 2 = (R + x) 2 + y 2 ,
F is the loading applied to the tablet, and h and R are respectively the height and the radius of the cylindrical sample. Interesting indicative values can be extracted at the center of coordinates (x = 0, y = 0) where the shear stress is theoretically equal to 0, and σ xx and σ yy are given by the following expressions:

σ xx (0, 0) = -3F πRh (11)
σ yy (0, 0) = F πRh (12)
Thus, the center of the cylindrical sample is associated to normalised compression (σ * xx = σ xx πRh/F) and tensile (σ * yy = σ yy πRh/F) stresses respectively equal to -3 and +1. This signifies that, at the center of the sample, the absolute stress value in the compression direction is three times the value in the tensile direction. Based on these criteria, detailed studies allowing for evaluating the validity of DEM calculations in the present context are given in previous contributions [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF][START_REF] Leclerc | Effect of packing characteristics on the simulation of elasticity and brittle fracture by the Cohesive Discrete Element Method[END_REF]. In the present work, our main objective is to investigate the ability of the DEM to reproduce the failure process as function of halo size. That is why, our efforts focus on σ yy stress field. Indeed, crack initiation theoretically occurs at the center of the sample when σ yy is equal to the tensile strength of the material. Figures 8a, b, c andd illustrate the influence of the halo size on the normalised tensile stress for a halo size varying from 1 to 15. For comparison purposes, numerical estimates obtained by the FEM using a mesh composed of about 100,000 tetrahedral elements are also provided. Numerical outputs exhibit that a halo size equal to 1 leads to very disrupted values which oscillate around the targeted indicator of 1. However, except for undesirable boundary effects, the greater the halo size is the smoother the stress field becomes. Thus, a kind of convergence is observed for which halo sizes of 10 and 15 give the best results in comparison with theoretical expectations and numerical estimates obtained by the FEM. Please notice that this finding is in good agreement with the conclusions of Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF]. Nevertheless, as suggested by the same authors, one has to take care to avoid too large halo sizes, typically higher than 20, which could lead to a loss of information where local gradients occur if the density of particles is not defined accordingly.

4.3.

Failure process without halo 4.3.1. Griffith's failure criterion Griffith developed a theory of brittle fracture based on the pre-existence of a population of random sub-microscopic cracks and flaws in the material. In such a paradigm, stress concentrations at the tips of microscopic cracks lead to their growth over the time which produces an increase in surface energy related to the formation of new free surfaces. The source of the increase in energy surface is the elastic energy which is released during crack opening. Thus, a crack will propagate as long as the energy balance between increased surface energy and released elastic energy is verified, which is characterised by a given critical fracture stress equal to the ultimate tensile strength for a pure mode I propagation. In the present context, Griffith's failure criterion [START_REF] Fairhurst | On the validity of the "Brazilian" test for brittle materials[END_REF] reads:

         σ I = UTS if (3σ I + σ III ) ≥ 0 (σ I -σ III ) 2 + 8UTS(σ I + σ III ) = 0 if (3σ I + σ III ) < 0 ( 13
)
where UTS is the ultimate tensile strength, and σ I and σ III are major and minor principle stresses respectively.

Theoretical and numerical failure pattern

From a theoretical standpoint, crack should initiate at the center of the cylindrical sample, where shear stress is equal to 0 and the tensile stress is maximum, and propagate along the compression line. However, experimental observations [START_REF] Wang | The flattened Brazilian disc specimen used for testing elastic modulus, tensile strength and fracture toughness of brittle rocks: analytical and numerical results[END_REF] exhibited that stress concentrations located near the contact area between the loading device can affect the failure process. Thus, it is not uncommon to observe crack initiation close to contact area. Typical solutions consist in enlarging the surface application with or without flattening of the sample, or drilling a small hole at the center of the sample. From a numerical standpoint, boundary conditions affect the process in a similar way so that calculations led by both continuous and discrete approaches require analogous solutions. Thus, in a previous contribution, Leclerc [START_REF] Leclerc | Effect of packing characteristics on the simulation of elasticity and brittle fracture by the Cohesive Discrete Element Method[END_REF] performed discrete simulations to simulate the diametral compression of a ceramic tablet using the Removed Discrete Element Failure (RDEF) process developed by André et al. [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF][START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF] the concept of which is to remove a DE from the particulate system when a given stress criterion is reached. Based on the Griffith criterion discussed previously, Leclerc investigated the failure process for two configurations, a perfect cylinder and a flattened one, with particulate systems of characteristics Z = 7.5 and φ = 0.64. In the first case, results exhibited that high stress concentrations close to contact areas lead to improper cracks initiation and propagation (Figure 9a). In the second case, stress fields were much more diffuse and the crack onset occurred at the center of the sample as expected. On the other hand, the crack pattern characterised by triangular-like shapes at the ends was clearly influenced by the flattening part (Figure 9b). Thus, these investigations highlighted some difficulties to reproduce a proper Brazilian test using DEM. Besides, as discussed by the author, DEM results are very sensitive to packing characteristics which directly impact on the level of dispersion of the stress field. Stress dispersion affects the simulation by DEM of mechanical problems in which failure process is taken into account. Indeed, local stress peaks could lead to undesirable and premature cracks initiation with an acceleration of the macroscopic failure. One solution could consist in largely increasing the targeted cardinal number by enlarging the range of interaction controlled by parameter. However, preliminary tests exhibited that very large cardinal numbers would be required while a small increase in cardinal number is sufficient to cause sizeable additional computational cost due to supplementary number of contacts. That is why, it is preferred to consider a second solution which aims the stress determination processing based on Zhou's formulation [START_REF] Zhou | A new look at the atomic level virial stress : on continuum-molecular system equivalence[END_REF] which intrinsically introduces a range of interaction. In the present work, we expect to consider the symmetrised formulation given in Equation 8and control the domain of interaction using the halo approach. The main issue is that cracks propagation affects the connectivity within the discrete domain in the sense that a number of particles belonging to a halo could be separated from the central particle of the same halo by a crack (Figures 10a andb). Thus, the set of particles associated to each halo must be redefined over the time to take into account such effects and ensure that a physical connectivity exists between the central particle and each element of the corresponding halo. In this contribution, we propose to consider the well-known Depth-First Search (DFS) algorithm to respond to this objective. The idea is that, at a given time and knowing the position of removed particles at the previous step, a specific version of DFS algorithm is handled to verify the connectivity within each involved halo by searching the contact network. Thus, in case of local fracture, several connex components could be distinguished within a same halo (Figure 10c). A specific treatment is then applied to reduce the investigated halo to the connex component to which belongs the central particle (Figure 10d).

Halo size effects

We now aim at investigating the influence of the halo size on the failure process. For that purpose, we simulate cracks initiation and propagation within the tablet using the RDEF process and Griffith's criterion which were discussed in subsection 4.3. Please note that the RDEF process is applied at the scale of the particle and not at the scale of the halo. In a first step, on the basis of findings of studies led in subsection 4.2, we consider the case of a halo size set to 15 and we examine the failure process over the time. However, we must precise that, due to the poor efficiency of For information purposes, all removed particles are depicted in blue at their initial position in the particulate system to better highlight the crack propagation but these have strictly no effect on the failure process after deletion. Results exhibit that crack onset occurs at the center of the sample for a strain e = 0.001491. Cracks then propagate in the compression direction until reaching contact areas at both opposite side of the cylindrical sample for e = 0.00153. Triangular-like crack ends are still evolving until the final pattern is obtained for e = 0.001633. Thus, for a qualitative standpoint, the numerical failure process meets our expectations. Besides, the final crack pattern is much less diffuse than ones obtained without halo and illustrated on Figures 9a andb. However, the crack width remains quite large with about 9-11 particles in its width, which is less than the absolute halo size but denotes nevertheless some scale transition effects.

Figure 12 illustrates the evolution of the loading curve for several halo sizes in the 1-20 interval. For comparison purposes, based on the assumption of a brittle material, the analytical solution is also depicted in red dashed line.

Results exhibit a kind of convergence to theoretical reference which is only reached for a halo size of 15 while a too small halo size leads to premature failure of the sample. Thus, for example, a halo size of 5 is related to a maximum loading of 2452 N and a corresponding strain of 6.81×10 -4 which is very far from the analytical solution characterised by a maximum loading of 5931 N and a strain of 1.52×10 -3 . Conversely, a too large halo size of 20 provides a slightly overestimated maximum loading of 6116 N for a strain of 1.59 × 10 -3 . Thus, the sole case of a halo size of 15 gives suitable outcomes with a maximum loading of 5652 N for a corresponding strain of 1.49 × 10 -3 . At this stage, it could be reasonable to think that a simple shift suffices to meet our expectations using small halos. However, as shown by inset figures illustrating the final failure pattern for each investigated configuration, a halo size less than 10 provides improper failure process and crack pattern. As a result, the halo approach turns out to be quite relevant to model the failure process during a diametral compression test. Nevertheless, it requires a certain care in the determination of the halo size. Indeed, the smallest admissible halo size has to be targeted to avoid prohibitive additional calculation cost mainly due to stress evaluation at each time step (the specific DFS algorithm based treatment is in fact very lightweight). Moreover, Figure 13 shows the evolution of the rate of removed particles, as function of the applied strain for a halo size equal to 15. Results exhibit a sharp increase around a strain of 1.5 × 10 -3 which is in good agreement with the numerical loading curve and the theory. Please notice that the set of removed particles is also distinguished between tensile, shear and compression modes of failure. Thus, numerical outputs show that shear and tensile modes are predominant with a similar weight. More in-depth investigations enable to observe that shear failure prevails close to contact areas while the tensile failure is more prone to arise near the center of the specimen.

Application to the failure of heterogeneous media

This section presents an application of the halo approach to the simulation of the failure process of heterogeneous media. We consider the context of a UniDirectional (UD) fibre composite composed of a bio-based epoxy matrix reinforced with flax fibres. Two configurations are studied. In the first one, a simple model constituted of a single fibre embedded in the matrix is investigated. Then, in a second step, image-based techniques are used to generate a Representative Volume Element (RVE) of the material containing 8 UD fibres. In both cases, transverse tensile tests are set up and the numerical failure process is studied, under the assumption of a brittle matrix, as function of the halo size.

Flax/bio-based epoxy composite

We consider the context of a laminate composite constituted of a set of layers of fibrous material. In the present work, each individual layer is a UD composite composed of a bio-based epoxy as matrix and aligned flax fibres as reinforcement. Thus, each layer has the same orientation and composition. For ecological purposes, we have selected the InfuGreen 810 epoxy provided by Sicomin, which has a lower environmental impact than standard epoxy resins. Unidirectional flax fibre yarn bundles commercialised as FLAXPLY UD 180 by LINEO and specifically treated for high performance composite applications have been used as reinforcement. Figure 14 shows two SEM images of the composite material in longitudinal and transverse planes. The macroscopic longitudinal Young's modulus of an impregnated flax fibre strand is evaluated from the macroscopic longitudinal Young's modulus of the composite, which has preliminary been determined with the help of tensile tests and is equal to 19.7 GPa ± 1.8 GPa, using the rule of mixture. Under the simplifying hypotheses discussed above, we assimilate the anisotropic fibre strand to a homogeneous isotropic material of Young's modulus equal to its longitudinal Young's modulus so that E f = 39.7 GPa. The Poisson's ratio ν f is a priori unknown and set to 0.28. The corresponding microscopic properties are E µ, f = 4387.9 GPa and r µ, f = 0.15.

Please notice that the volume fraction of flax fibre strands has been measured using SEM images to be 45%.

Single-fibre model 5.2.1. Configuration

In a first approach, we model the flax fibre/bio-based epoxy composite using a cubic representative pattern of length L= 234 µm. This pattern, which can be seen as an RVE of the heterogeneous material, is composed of one single fibre described by a perfectly aligned cylinder of diameter D= L/1.3 = 180 µm the axis of which is parallel to x direction. Thus, the volume fraction of flax strands is equal to 45% in agreement with experimental measurements. DEM simulations are carried out using particulate systems composed of 250,000 particles according to the discussion of subsection 2.3 and cohesive beam elements the formulation of which is based on microscopic parameters associated to each material. Numerical tensile tests are performed in the transverse direction of fibres by applying an imposed velocity v= 2.34 × 10 -5 µm/s in the y direction. Due to heterogeneous nature of the medium, particles could be associated to two possible phases, fibre or matrix, according to the position of their centroid. Interfacial parameters have to be also introduced to describe cohesive beam elements linking two particles located in different materials. Thus, in order to ensure a smooth transition between each material, we define the microscopic interfacial Young's modulus E µ,int and the interfacial dimensionless radius r µ,int as the arithmetic average of their counterparts :

E µ,int = E µ,m + E µ, f 2 (14)
r µ,int = r µ,m + r µ, f 2 (15)

Stress fields

We aim at investigating the influence of the halo size on the stress and its level of dispersion in the case of a heterogeneous medium. Based on tensile tests previously described, stress fields are extracted from DEM calculations at a given applied strain e=0.0005 for several halo sizes in the 1-20 interval. Please notice that in the present context of a composite material, the halo definition takes into account the heterogeneity of the medium so that all particles related to a given halo belongs to the same phase. Thus, the geometrical shape of the halo could only be a part of a sphere near the fibre-matrix interface. Figures 15a andb illustrate σ zz stress fields in the plane x=L/2 obtained for halo sizes respectively equal to 5 and 15. For information purposes, the case of a halo size equal to 1 has not been depicted since it leads to an unsuitable noise-like stress field, which does not allow for relevant comparisons. Besides, FE estimates determined in the same configuration using a mesh composed of 150,000 tetrahedral elements are provided (Figure 15c). From a qualitative standpoint, results show that the DEM with a halo size of 15 and the FEM produce very similar stress fields. Conversely, the case of a small halo size equal to 5 exhibits a quite disturbed stress field with improper local stress peaks. From a quantitative standpoint, indicative values are extracted at A(0.5L,0.942L,0.5L), B(0.5L,0.5L,0.5L), C(0.5L,0.942L,0.942L) and D(0.5L,0.5L,0.942L) positions. Globally, σ zz stress values obtained using a halo size of 15 are close to FE estimates while a halo size of 5 leads to quite inaccurate results due to local dispersion. Maximum and minimum values given by DEM calculations using a halo size of 15 are respectively 3.04 and -1.87 MPa. FEM calculations provide a quite close minimum value equal to -2.04 MPa and a higher maximum stress equal to 5.02 MPa. One could conclude that a halo size of 15 suffices to unsuitably smooth out some relevant stress gradients. Nevertheless, this finding must be tempered since extreme stress values given by FEM turn out to be very sensitive to mesh refinement in the present case. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF] which conclude that too large halos (larger than 25) could be associated to an excessive smoothing and a loss of information, we deliberately limit our investigations to the interval [START_REF] Krüger | The virtual crack closure technique : History, approach and applications[END_REF][START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF] of halo sizes. For comparison purposes, FEM estimates based on a mesh composed of 150,000 tetrahedral elements are also provided. Besides, each data dot represents the mean of a set of four numerical predictions obtained using four different particulate systems of same characteristics (number of particles, Z, φ). Results exhibit that, whatever the position, σ yy and σ zz stresses converge to FE data with an acceptable margin of absolute error less than 0.5 MPa. Thus, a too small halo (smaller than 5) is related to unsuitable values while large halos (strictly larger than 5) lead to proper ones, and the larger the halo is more accurate the results are. Based on this scope of data, we can conclude that a halo size in the 10-15 interval is probably a good compromise between accuracy and computational cost. These findings are in good agreement with those of subsection 4.2 in the context of the diametral compression test.

Failure process

DEM simulations are now carried out by taking into account cracks initiation and propagation in the matrix. Please notice that, for simplifying purposes, flax fibres are supposed unbreakable and the bio-based epoxy resin is supposed perfectly brittle. Based on these assumptions, material failure is modelled using the RDEF process [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF][START_REF] Jebahi | Simulation of Vickers indentation of silica glass[END_REF] discussed in subsection 4.3 and considering the classical Rankine criterion (or maximum stress criterion) for brittle materials: where UTS stands for the Ultimate Tensile Strength of the bio-based epoxy resin which is 60 MPa. σ I is the maximum normal stress which is determined using the halo approach coupled to the specific halo treatment discussed in subsection 4.4. On the basis of previous results related to stress field determination, we first study the failure process for a given halo size set to 15. Figures 17a,b,c,d, e and f illustrate the evolution of the crack pattern as function of the applied strain e. Two cracks initiate at both opposite inner sides of the matrix, where maximum stresses are localised, at a given applied strain equal to 0.00248. Then, the cracks progressively propagate symmetrically along the z direction until they reach outer surfaces at an applied strain of 0.01028. These results exhibit a suitable cracks propagation consistent with a mode I opening. However, as observed in the case of the diametral compression test and despite the fact that the failure criterion is applied at the scale of the particle and not the halo, the crack thickness is quite large and also irregular with from 7 to 12 particles in its width. This is less than the absolute size of the halo which is here set to 15 but nevertheless denotes some scale effects. A too small halo of size less than 5 leads to unsuitable cracks initiation and propagation near the loaded surfaces, which is related to the presence of undesirable stress gradients at these specific areas. A proper crack onset is observed in the case of a halo size of 10 but the propagation only occurs at one side of the inclusion so that the final pattern is asymmetrical. Figure 19 illustrates the evolution of the rate of removed particles as function of the applied strain for the same halo size of 1, 5, 10 and 15. Results show that, except for the case of a halo size of 1 which highlights an accelerated process, failure starts to initiate for a strain close to 0.0023. Furthermore, the larger the halo is the higher the number of removed particle is. We hypothesise that this could be related to the width of cracks which is sensibly larger for bigger halos. Nevertheless, the failure process strongly depends on the halo size so that the number and the length of cracks could be also affected. As a result we can conclude that, in accordance with the findings of the previous investigation of stress field and those of subsection 4.4 in the context of the diametral compression test, a halo size less than 10 leads to improper cracks onset and propagation. In the present case of heterogeneous medium, the sole configuration of a halo size equal to 15 provides a suitable failure pattern so that a minimal admissible size could be found in the 10-15 interval. For that purpose, a 4-step specific process is followed. First, the methodology consists in choosing the number of inclusions. SEM images allow to extract representative patterns composed of a maximum of 16 fibres. However, due to the limitations in terms of computational cost, we prefer to only handle 8 inclusions which is a good compromise between accuracy and efficiency. Figure 20a illustrates the sample extracted from a given SEM image of the flax/biobased resin composite. In a second step, based on a 8-bit grayscale image which contains 256 colors, a two-level thresholding technique is set up and results in a binary image (Figure 20b). This latter can not be directly exploited due to the presence of flaws and irregularities. Besides, each strand is composed of a set of flax fibres which are only taken into account as a whole in our modelling. That is why, a Gaussian filter is applied to reduce the image noise and undesirable details, and both operations (thresholding and Gaussian filter) are repeated until a satisfactory result is obtained. Then, as depicted on Figure 20c, an edge detector algorithm is carried out to well separate each body. This could be the last step before DEM simulations. Nevertheless, for practical reasons, we opt for a circular representation of strands with respect to their centroid and surface. Thus, Figure 20d illustrates the final multi-inclusion RVE the volume fraction of which is 47.4%, what is close to the initial percentage of 45%. DEM simulations are performed using the same assumptions as previously seen in the case of the single-inclusion model, except that the number of particles is increased to 1 million DE to ensure a proprer representation of each inclusion. For benchmarking purposes, FE predictions obtained in the same configuration using the above-mentioned mesh composed of 350,000 hexahedra are also provided. Results show that except for too small halo sizes less than 5 which lead to unsuitable stress predictions, the larger the halo is, more accurate σ yy and σ zz stress values become. Whatever the position, a near convergence is observed and the absolute difference between FE and DE results is always less than 0.4 MPa for a halo size of 15. Thus, for example, at A position σ yy and σ zz stress values vary respectively from 27.9 to 14.5 MPa and from -14.6 to -2.57 MPa as the halo size varies from 1 to 20. Corresponding FE values are respectively 14.2 MPa and -2.57 MPa. This highlights the ability of the halo approach to provide suitable stress fields in the present configuration. However, as already discussed in subsection 5.2, we remind that a halo size bigger than 25 could lead to an excessive smoothing and a loss of information, which signifies that the notion of convergence must be tempered beyond a halo size of 20. Besides, for large halos containing a consequent number of particles, the step of stress determination requires a sizeable additional calculation cost. Thus, the suitable halo size must be neither underestimated nor overestimated to ensure a good compromise between accuracy and calculation cost, what is verified for a halo size comprised between 10 and 15. This finding is in accordance with the conclusions of the diametral compression test and the previous study of a single fibre composite. In such an approach, the maximum stress criterion is applied at the scale of the particle using the stress determined at the scale of the halo. We consider the same configuration as previously investigated in the context of the single fibre composite. The bio-based epoxy resin is modelled as a perfectly brittle material with a UTS equal to 60 MPa and flax fibre strands are supposed unbreakable. Figures 23a, b, c andd show the influence of the halo size on the final crack pattern of the multi-fibre composite. Please notice that the failure pattern is produced at different values of the applied strain varying from 0.00179 to 0.00698 as the halo size varies from 1 to 15. The case of a halo size of 1 leads to an unsuitable crack pattern which is clearly affected by stress peaks located near the upper face. In all other cases, results exhibit a failure process consistent with a mode I opening and in which cracks initiate at the fibre/matrix interface and propagate along the z direction until reaching the edges of the RVE. Some variations can nevertheless be noticed. Halo sizes of 5 and 10 provide a very similar crack pattern but the one of a halo size of 15 is slight different what highlights the predominant role of the stress dispersion on cracks onset and nucleation phenomena. Besides, the crack spacing strongly depends on the halo size so that a halo size of 5 leads to cracks with from 2 to 5 particles in the width and a halo size of 15 produces cracks with from 6 to 12 particles in the width. Thus, in the present case, intermediate halo sizes of 5 and 10 exhibit more realistic cracks than a larger halo size of 15. Figure 24 illustrates the influence of the halo size on the loading curve. As observed in the case of the single inclusion problem, a halo size of 1 is related to an accelerated failure due to the high level of stress dispersion which induces undesirable stress gradients. Maximum loadings are respectively 5.353 N, 13.48 N 14.49 N and 15.16 N for halo sizes of 1, 5, 10 and 15 and corresponding strains of 8.975 × 10 -4 , 2.308 × 10 -3 , 3.783 × 10 -3 and 2.188 × 10 -3 . Thus, a near convergence is again observed except that the halo size of 15 exhibits a premature peak in comparison with halo sizes of 5 and 10. This could be due to the difference of failure process discussed previously but more in-depth investigations are required to support this hypothesis. Figure 25 illustrates the evolution of the rate of removed particles with respect to the total number of DE as function of the applied strain for halo sizes in the 1-15 interval. Numerical outcomes show that, with the exception of halo size of 1 which leads to an unsuitable failure process, cracks initiate at a strain of about 0.002 and the larger the halo is the higher the rate of removed particle is. This is probably connected to cracks width issue observed for larger halos. As a result, we conclude that in the present context of a multi-fibre composite a halo size in the 5-10 interval suffices to provide a suitable failure pattern. Nevertheless, in light of previous investigations of the diametral compression test and the single-fibre composite, this finding has to be tempered since a halo size less than 10 leads to an improper loading curve.

σ I = UTS (16)

Conclusions and prospects

The present contribution was dedicated to the simulation of cracks initiation and propagation in 3D DEM simulations of continuous media by using the halo approach introduced by Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF]. In a first part, after a brief reminder of the hybrid lattice-particle model which allows to simulate a continuous medium using DEM, the level of stress dispersion was investigated as function of the cardinal number and the halo size. Results exhibit two main findings. Firstly, as expected, a higher coordination leads to a lower stress dispersion characterised by a small CoV σ . Thus, CoV σ varies respectively from 120% to 43% and from 67% to 29% for Z = 6.2 and Z = 7.5 as r µ varies from 0.1 to 0.9. Secondly, whatever r µ , the halo approach enables to efficiently reduce the level of dispersion. For example, in the case r µ = 0.6, CoV σ decreases from 35,6% to 1.08% as the halo size increases from 1 to 15. In a second step, the diametral compression of a cylindrical tablet composed of a brittle marble was studied. As concluded by Moukadiri et al. [START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF] in the same context, numerical simulations showed that a halo size between 10 and 15 provides accurate stress values with respect to FEM calculations. Based on the RDEF process and a specific halo treatment aiming at ensuring the connexity within each halo during a failure process, cracks onset and propagation were simulated as function of the halo size. In agreement with the results obtained for stress fields, both cracks process and loading curve turned out to be unsuitable for a halo size less or equal than 10. Nevertheless, even though a halo size of 15 exhibits proper outputs, the crack width is also affected with about 9-11 particles in the width which denotes some scale effects. Finally, the case of a flax fibre/bio-based composite material was explored through two configurations, a UD fibre composite composed of a single inclusion and an RVE composed of 8 inclusions which was generated from a SEM image. In both cases, numerical simulations of quasi-static transverse tensile tests confirmed that a halo size between 10 and 15 produces a proper stress field. However, examination of brittle failure process based on Rankine As observed in cases of Brazilian test and single-fibre composite, the halo approach turns out to be an interesting method to ensure suitable DEM simulations of failure process in homogeneous and heterogeneous materials. However, one must take care of avoiding too large halo sizes which can lead to unrealistic crack width, expecially when cracks nucleation occurs as observed for the multi-fibre composite. In fact, we hypothesize that scale effects are related to the stress smoothing that leads to a local chain reaction around a given crack onset. This could potentially be better controlled using a high Rayleigh damping and a smaller time step but at the direct expense of computational cost. In a next future, we expect to take benefit of the halo approach to simulate non-linear stress-strain behaviours and model the brittle ductile failure of homogeneous and heterogeneous materials using DEM. Besides, we are also currently working on the development of anisotropic fibre models.

Figure 1 :

 1 Figure 1: Influence of ε sensitivity parameter on the local cardinal number Z (a) case ε = 0 (b) case ε > 0

Figure 2 :

 2 Figure 2: Cohesive beam model

  depict the evolution of elastic coefficients E M and ν M as function of the number of particles for two materials, a bio-based epoxy resin and a brittle marble respectively, and a number of particles in the interval [10,000;1,500,000]. Macroscopic elastic properties of the bio-based epoxy resin are E M = 3.0 GPa and ν M = 0.3 which correspond to microscopic parameters E µ = 883.4 GPa and r µ = 0.1. For the brittle marble, E M = 17.3

Figure 3 :

 3 Figure 3: Evolution of macroscopic elastic properties E M and ν M as function of the number of particles

Figure 4 :Figure 5 :

 45 Figure 4: Evolution of CoV σ as function of r µ for Z = 6.2 and Z = 7.5, and φ set to 0.64

Figure 6 :

 6 Figure 6: Evolution of CoV σ as function of (a) the halo size and (b) r µ

Figure 7 :

 7 Figure 7: Numerical configuration of the diametral compression test

Figure 8 :

 8 Figure 8: Influence of the halo size on normalised σ yy stress field in the compression line (a) R h I /R h I = 1 (b) R h I /R h I = 5 (c) R h I /R h I = 10 (d)

Figure 9 :

 9 Figure 9: Crack pattern obtained by numerical simulation of a diametral compression test (a) Case of a perfect sample (b) Case of a flattened sample (Figures taken from [33])

Figure 10 :

 10 Figure 10: DFS algorithm based halo treatment

Figure 11 :

 11 Figure 11: Evolution of the crack pattern as function of the applied strain e

Figure 12 :Figure 13 :

 1213 Figure 12: Influence of the halo size on the loading curve and the crack pattern

Figure 14 :

 14 Figure 14: SEM images of the flax/bio-based epoxy composite material in longitudinal and transverse planes

Figure 15 :

 15 Figure 15: Influence of the halo size on σ zz stress field (expressed in MPa) in the plane x=L/2 and comparison with FE calculations Figures 16a and b illustrate the evolution of σ yy and σ zz stress values as function of the halo size. Based on the findings of Moukadiri et al.[START_REF] Moukadiri | Halo approach to evaluate the stress distribution in 3D Discrete Element Method simulation : Validation and application to flax/bio based epoxy composite[END_REF] which conclude that too large halos (larger than 25) could be associated to an excessive smoothing and a loss of information, we deliberately limit our investigations to the interval[START_REF] Krüger | The virtual crack closure technique : History, approach and applications[END_REF][START_REF] Leclerc | Discrete Element Method to simulate the elastic behaviour of 3D heterogeneous continuous media[END_REF] of halo sizes. For comparison purposes, FEM estimates based on a mesh composed of 150,000 tetrahedral elements are also provided. Besides, each data dot represents the mean of a set of four numerical predictions obtained using four different particulate systems of same characteristics (number of particles, Z, φ). Results exhibit that, whatever the position, σ yy and σ zz stresses converge to FE data with an acceptable margin of absolute error less than 0.5 MPa. Thus, a too small halo (smaller than 5) is related to unsuitable values while large halos (strictly larger than 5) lead to proper ones, and the larger the halo is more accurate the results are. Based on this scope of data, we can conclude that a halo size in the 10-15 interval is probably a good compromise between accuracy and computational cost. These findings are in good agreement with those of subsection 4.2 in the context of the diametral compression test.

Figure 16 :

 16 Figure 16: Evolution of σ yy and σ zz stress values (expressed in MPa) as function of the halo size

Figure 17 :

 17 Figure 17: Evolution of the crack pattern of the single-fibre model as function of the applied strain e

Figure 18 :Figure 19 :

 1819 Figure 18: Influence of the halo size on the loading curve and the crack pattern of the single-fibre model

Figure 20 :

 20 Figure 20: Image-based generation of an RVE of the flax/bio-based epoxy resin composite

Figure 21 :

 21 Figure 21: Influence of the halo size on σ zz stress field (expressed in MPa) in the plane x=L/2 and comparison with FE calculations (case of multi-fibre model) Figures 22a and b illustrate the influence of the halo size on σ yy and σ zz stress values at A, B, C and D positions.For benchmarking purposes, FE predictions obtained in the same configuration using the above-mentioned mesh composed of 350,000 hexahedra are also provided. Results show that except for too small halo sizes less than 5 which lead to unsuitable stress predictions, the larger the halo is, more accurate σ yy and σ zz stress values become. Whatever the position, a near convergence is observed and the absolute difference between FE and DE results is always less than 0.4 MPa for a halo size of 15. Thus, for example, at A position σ yy and σ zz stress values vary respectively from 27.9 to 14.5 MPa and from -14.6 to -2.57 MPa as the halo size varies from 1 to 20. Corresponding FE values are respectively 14.2 MPa and -2.57 MPa. This highlights the ability of the halo approach to provide suitable stress fields

Figure 22 :

 22 Figure 22: Evolution of σ yy and σ zz stress values (expressed in MPa) as function of the halo size (case of multi-fibre model)

Figure 23 :

 23 Figure 23: Influence of the halo size on the crack pattern of the multi-fibre composite

Figure 24 :Figure 25 :

 2425 Figure 24: Influence of the halo size on the loading curve of the multifibre model

Data availability statement

The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations.

Conflicts of interest

None declared.