A Ammar

W Leclerc

M Guessasma

N Haddar

Discrete element approach to simulate debonding process in 3D short glass fibre composite materials: application to PA6/GF30

Keywords: Discrete Element Method, Short Fibre Reinforced Composites, Numerical approach, Interfacial debonding, Delamination

This paper discusses the development of a numerical approach based on the Discrete Element Method (DEM) to simulate the interfacial debonding in Short Fibre Reinforced Composites (SFRC). For that purpose, we consider a hybrid lattice-particle approach based on a cohesive beam model. Already used to study composite materials, we aim to extend this concept to investigate the damage mechanism of SFRC in particular the PA6/GF30. To do so, a specific technique to estimate the Representative Elementary Volume (REV) of PA6/GF30 is proposed and discussed. The suggested discrete modelling is then validated by comparison with micromechanical approaches and experimental data for elastic behaviour predictions. At this time, the development of a DEM-based framework to model SFRC composite represents one of the novelties of this work. A mixed-mode discrete damage model based on an energetic formulation (which constitutes another originality of this paper) is subsequently discussed. To do so, the cohesive beam model is first implemented and validated through standard delamination problems. Thereafter, the case of SFRC is explored, focusing on the case of UniDirectional (UD) single fibre and aligned short-fibre composites under bending test. Simulation results in terms of onset and propagation of the debonding of fibre/matrix interface confirm the potential of the discrete modelling presented.

Introduction

The transport sector, particularly the automotive industry, is going through a deep change as regards the environmental and social challenges it poses. The concept of sustainable mobility has become a major issue for everyone, car manufacturers as well as economic actors. Recent Conference Of Parties (COP21) legislations restrict pollutant emissions of light vehicles [START_REF] Wang | Tweeting the united nations climate change conference in paris (cop21): An analysis of a social network and factors determining the network influence[END_REF]. As a result, working on vehicle weight reduction has become necessary in order to meet the new COP21 standards. Polyamide 6 reinforced with 30% of short glass fibre (PA6/GF30) is an example of Short Fibre Reinforced Composites (SFRC) used in the automotive sector. It represents an alternative to metal structures as part of the engine compartment and interior equipment of vehicles thanks to its duality: mechanical strength-lightness. However, this composite is subjected to several environmental conditions that affect its mechanical behaviour and durability [START_REF] Ksouri | Effects of processing steps and hygrothermal ageing on mechanical performance of PA6GF30 composite: Interfacial shear strength[END_REF]. For example, under-the-hood automotive parts such as the air filter, the water coolant reservoir or the air intake manifold are subjected to permanent vibrations and are therefore exposed to cyclic loading. This phenomenon can influence the fatigue life, including mechanical properties at the microstructural level up to the damage of the composite material. Moreover, on European high speed railway lines, some parts manufactured using PA6/GF30 are affected by high humidity conditions which influence the fracture behaviour of this material [START_REF] Ferreno | Characterisation by means of a finite element model of the influence of moisture content on the mechanical and fracture properties of the polyamide 6 reinforced with short glass fibre[END_REF].

In order to obtain enough information and understand different phenomena related to SFRC PA6/GF30, several experimental tests have been conducted. Fatigue failure mechanisms in PA6/GF30 were characterised using experimental methods in [START_REF] Horst | Mechanisms of fatigue in short glass fiber reinforced polyamide 6[END_REF][START_REF] Horst | Fatigue fracture mechanisms and fractography of short-glass-fibre-reinforced polyamide 6[END_REF]. Meraghni et al. [START_REF] Meraghni | Parameters identification of fatigue damage model for short glass fiber reinforced polyamide (PA6-GF30) using digital image correlation[END_REF] presented the strategy for identifying damage parameters using fatigue tests. Detailed investigations on mechanical properties and fracture mechanisms of PA6/GF35 were carried out using tensile tests and three point bending tests on Single Edge Notch Bend (SNEB) samples [START_REF] Ferreno | Characterisation by means of a finite element model of the influence of moisture content on the mechanical and fracture properties of the polyamide 6 reinforced with short glass fibre[END_REF]. Nciri et al. [START_REF] Nciri | Modelling and characterisation of dynamic behaviour of short-fibre-reinforced composites[END_REF] characterised the microstructure of short glass fibre reinforced polypropylene composites using micro-computed microtomography. Ksouri et al. [START_REF] Ksouri | Effects of processing steps and hygrothermal ageing on mechanical performance of PA6GF30 composite: Interfacial shear strength[END_REF] focused on the hygrothermal ageing of the PA6/GF30 and its effects on mechanical performances. For that purpose, they identified the InterFacial Shear Strength (IFSS) and its consequences on the PA6/GF30 frictional effects, using the Kelly-Tyson model, under hygrothermal ageing. Based on IR spectroscopy and using experimental techniques, the same authors pointed out that long term ageing affects the fibre/matrix interface, reducing adhesion between fibre and matrix [START_REF] Ksouri | Long term ageing of polyamide 6 and polyamide 6 reinforced with 30% of glass fibers: physicochemical, mechanical and morphological characterization[END_REF]. In fact, swelling phenomena may occur in relation to moisture effects which favour the interfacial debonding damage [START_REF] Arif | In situ damage mechanisms investigation of PA66/GF30 composite: Effect of relative humidity[END_REF]. In another context, during fatigue loading [START_REF] Arif | Multiscale fatigue damage characterization in short glass fiber reinforced polyamide-66[END_REF] or bending tests [START_REF] Sato | Microfailure behaviour of randomly dispersed short fibre reinforced thermoplastic composites obtained by direct SEM observation[END_REF] damage in SFRC occurs at the microscopic scale precisely at the fibre/matrix interface. Moreover, observations show that under the transverse tensile load of the glass-polymer composite material [START_REF] Vajari | Effect of fiber positioning on mixed-mode fracture of interfacial debonding in composites[END_REF], cracks usually initiate and propagate at the fibre/matrix interface. The authors related this phenomenon to the weaker interface strength compared with the matrix and the fibre. Agrawal et al. [START_REF] Karlsson | Obtaining mode mixity for a bimaterial interface crack using the virtual crack closure technique[END_REF] highlighted the fact that interface cracking in bimaterials always involves mixed mode conditions even if macroscopically pure mode I or pure mode II failure is observed. Hence, experimental methods may be limited in identifying different parameters associated with mixed mode. Thus, a multi-scale characterisation is required to understand specific phenomena such as the material fracture behaviour occurring at the local scale. Resort to micromechanical models and numerical simulations is consequently necessary to better understand such mechanisms and their influence at the macroscopic scale.

Notta-Cuvier et al. [START_REF] Notta-Cuvier | Modelling of progressive fibre/matrix debonding in short-fibre reinforced composites up to failure[END_REF] developed an approach to simulate progressive interfacial debonding in SFRC. The proposed failure criterion takes into account the critical amount of voids as well as the influence of environmental conditions at the fibre/matrix interface. A micromechanical approach, based on a generalised incremental Mori-Tanaka approach, was developed to model progressive damage at the interface of PA66/GF30 [START_REF] Despringre | Micromechanical modeling of damage and load transfer in particulate composites with partially debonded interface[END_REF]. Moreover, a hybrid model was proposed by Laribi et al. [START_REF] Laribi | Sheet Molding Compound Automotive Component Reliability Using a Micromechanical Damage Approach[END_REF] for Sheet Molding Compound (SMC) material. The latter combined a micromechanical approach, based on the Mori-Tanaka scheme, with a phenomenological model in order to investigate fibre/matrix interface decohesion under fatigue loading. Nevertheless, approaches presented in [START_REF] Notta-Cuvier | Modelling of progressive fibre/matrix debonding in short-fibre reinforced composites up to failure[END_REF] and [START_REF] Laribi | Sheet Molding Compound Automotive Component Reliability Using a Micromechanical Damage Approach[END_REF] require a large experimental database to provide more realistic predictions, representing the main drawback of these models. Besides, some limitations, such as sample preparation, testing procedure and measurement at the fibre scale, related to the use of the reverse engineering method to identify parameters at the fibre/matrix interface, reduce the scope of applications of the model presented in [START_REF] Despringre | Micromechanical modeling of damage and load transfer in particulate composites with partially debonded interface[END_REF]. Hence, current micromechanical models are not suited to dealing with complex debonding phenomena.

The literature provides a set of numerical approaches to simulate cracks and to deal with delamination or debonding problems. Among them, Virtual Crack Closure Technique (VCCT) is one of the most widely used procedures to predict delamination growth [START_REF] Krueger | The virtual crack closure technique for modeling interlaminar failure and delamination in advanced composite materials[END_REF]. This technique provides information regarding the onset and the propagation of delamination in composites [START_REF] Jimenez | Application of the finite-element method to predict the onset of delamination growth[END_REF]. However, a pre-crack must be defined for analysis, which makes this technique inefficient in predicting delamination initiation. Furthermore, interface delamination or debonding propagation requires an updating mesh which reduces its ability to treat complex geometries and some load cases [START_REF] Heidari-Rarani | Finite element modeling strategies for 2D and 3D delamination propagation in composite DCB specimens using VCCT, CZM and XFEM approaches[END_REF]. The Cohesive Zone Model (CZM) is used in numerical studies as a strong tool to model fracture of bimaterial interfaces particularly when Linear Elastic Fracture Mechanics (LEFM) can not be used. Unlike VCCT, CZM is able to predict the onset as well as the propagation of interface delamination without pre-crack [START_REF] Garg | Application of scaled boundary finite element method for delamination analysis of composite laminates using cohesive zone modelling[END_REF]. Thus, delamination growth was investigated in composite Double Cantilever Beam (DCB) specimens [START_REF] Heidari-Rarani | Finite element modeling of mode I delamination growth in laminated DCB specimens with R-curve effects[END_REF][START_REF] Fan | Cohesive zone with continuum damage properties for simulation of delamination development in fibre composites and failure of adhesive joints[END_REF]. Moreover, CZM was adapted to predict the debonding damage at the fibre/matrix interface under quasi-static transverse loading [START_REF] Hu | Stable numerical simulations of propagations of complex damages in composite structures under transverse loads[END_REF]. Among other applications, a unit cell-model taking into account CZM formulation was introduced by Tvergaard in order to simulate fibre/matrix decohesion in aligned short SiC fibre composite [START_REF] Tvergaard | Debonding of short fibres among particulates in a metal matrix composite[END_REF]. The fibre/matrix interface debonding in carbon nanotube reinforced composites was also analysed using CZM by Matveeva et al. [START_REF] Matveeva | Debonding at the fiber/matrix interface in carbon nanotube reinforced composites: Modelling investigation[END_REF]. However, CZM appears to be mesh dependent according to Harper and Hallett [START_REF] Harper | Cohesive zone length in numerical simulations of composite delamination[END_REF]. It requires a refined mesh at the interface to provide better accuracy of interface delamination phenomenon. Furthermore, the eXtended Finite Element Method (XFEM) [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] has been developed to investigate failure of composite materials [START_REF] Motamedi | 3D nonlinear XFEM simulation of delamination in unidirectional composite laminates: a sensitivity analysis of modeling parameters[END_REF]. However, XFEM requires a pre-existing crack which limits the application of this method. Besides, dealing with multiple fractures with joining and bifurcation cracks makes this method costly to track damage propagation [START_REF] Leclerc | On the suitability of a Discrete Element Method to simulate cracks initiation and propagation in heterogeneous media[END_REF].

As an alternative, the Discrete Element Method (DEM) gives an opportunity to solve mechanical problems where several scales and discontinuities are involved. Initially developed by Cundall and Strack to simulate geomechanics problems [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], the DEM was adapted to model continuous media using cohesive spring [START_REF] Potyondy | A bonded-particle model for rock[END_REF] or beam elements [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 behavior of 3D heterogeneous continuous media[END_REF] between particles in contact. Such an approach was used by Haddad et al. [START_REF] Haddad | Application of DEM to predict the elastic behavior of particulate composite materials[END_REF] to develop a 2D-DEM model to simulate the effective elastic properties of BaTiO 3 material. It was later extended by Leclerc [33] in 3D to model the elastic behaviour of particulate alumina/Al composite. Furthermore, 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] took advantage of the developed 3D-DEM approach and introduced a concept called Halo to control the level of stress dispersion at local scale in the context of flax/bio based epoxy composite material. The DEM was then implemented to evaluate thermal-induced damage [START_REF] Hassan | On the suitability of a 3D discrete element method to model the composite damage induced by thermal expansion mismatch[END_REF]. Moreover, it has also been used to simulate the initiation and propagation of cracks of homogeneous and heterogeneous materials. Jebahi et al. [START_REF] Jebahi | Simulation of vickers indentation of silica glass[END_REF] studied the cracking response of silica glass under Vickers indentation test. Damage mechanisms were also investigated in a carbon fibre/epoxy composite cell by Maheo et al. [START_REF] Maheo | A promising way to model cracks in composite using Discrete Element Method[END_REF]. Recently, Leclerc et al. [START_REF] Leclerc | Halo approach to model cracks initiation and propagation in 3D Discrete Element Method simulation of homogeneous and heterogeneous materials[END_REF] adapted the Halo approach to simulate damage process and to evaluate cracking patterns for composite materials.

In the context of delamination problems, several studies have been carrid out via DEM to simulate mode I, II as well as mixed mode delamination growth. In this case, mode I fracture was investigated in [START_REF] Yang | Modeling progressive delamination of laminated composites by discrete element method[END_REF][START_REF] Chen | A Particle-Based Cohesive Crack Model for Brittle Fracture Problems[END_REF] via DCB tests. End Loaded Split (ELS) [START_REF] Yang | Modeling progressive delamination of laminated composites by discrete element method[END_REF] and End Notched Flexural (ENF) [START_REF] Chen | A Particle-Based Cohesive Crack Model for Brittle Fracture Problems[END_REF] were used to reproduce mode II fracture. Also, Fixed-Ratio Mixed-Mode (FRMM) [START_REF] Yang | Modeling progressive delamination of laminated composites by discrete element method[END_REF] and Mixed Mode Bending (MMB) [START_REF] Chen | A Particle-Based Cohesive Crack Model for Brittle Fracture Problems[END_REF] were exploited to simulate mixed mode in composite laminate materials. Delamination tests are performed by Yang et al. [START_REF] Yang | Modeling progressive delamination of laminated composites by discrete element method[END_REF] where hexagonal configurations are considered. As well, Chen et al. [START_REF] Chen | A Particle-Based Cohesive Crack Model for Brittle Fracture Problems[END_REF] resorted to regular configurations of spheres to define the continuous domain. However, induced privileged directions could lead to undesirable effects in relation to the geometric anisotropy. For that purpose, these configurations may present a potential misunderstanding regarding the isotropy of the studied system in comparison to random particulate systems. The concept of cohesive zone model was exploited and adapted to deal with fibre/matrix interfacial decohesion in composite materials. In this context, a debonding phenomenon for single-and multi-fibre composites under transverse tension was investigated in [START_REF] Le | Modeling damages and cracks growth in composite with a 3D discrete element method[END_REF]. In this study, Breakable Bond Failure (BBF) and Removed Discrete Element Failure (RDEF) processes were applied in the case of a 3D UniDirectional (UD) fibre-reinforced composites to detect the interface degradation. In addition, Ismail et al. [START_REF] Ismail | Discrete element modelling of unidirectional fibre-reinforced polymers under transverse tension[END_REF] examined the progressive interfacial decohesion of a unidirectionally reinforced composite under transverse loading and via 2D-hexagonally packed particles configuration. The onset and propagation of interfacial debonding were detected through the microbond test of fibre reinforced composite [START_REF] Yang | Discrete element modeling of the microbond test of fiber reinforced composite[END_REF] and in composite laminate [START_REF] Ismail | A dem model for visualising damage evolution and predicting failure envelope of composite laminae under biaxial loads[END_REF]. Both studies used a 2D-Discrete Element (DE) model which may limit the description of the post-failure process.

The main objective of this paper is to develop a 3D numerical model based on DEM to predict the mechanical behaviour and the damage of SFRC PA6/GF30 under bending load. The present work is a part of a methodology using 3D DEM to simulate fractures in composite materials by taking into account multi-scale aspects and multi-physical coupling. Previous works have dealt with debonding effects induced by thermal expansion phenomena [START_REF] Hassan | On the suitability of a 3D discrete element method to model the composite damage induced by thermal expansion mismatch[END_REF], and matrix cracking using RDEF process [START_REF] Leclerc | Halo approach to model cracks initiation and propagation in 3D Discrete Element Method simulation of homogeneous and heterogeneous materials[END_REF]. Through the present work, we aim to set up a specific DEM-approach based on the Representative Elementary Volume (REV) concept in order to model SFRC PA6/GF30. To do so, we consider DEM based on a hybrid particulate-lattice model in which cohesive elements are used to model the continuous and heterogeneous medium [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. For information purposes, all DEM simulations are carried out using MULTICOR 3D++ code based on C++ programming and developed in LTI laboratory [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]. The main benefit of such a model compared to a full lattice approach is that it also allows binary collisions between particles to be handled which occur during interfacial debonding. In addition, we consider cohesive beam elements which provide more realistic crack patterns than spring ones under several mechanical loadings [START_REF] Schlangen | Fracture simulations of concrete using lattice models: computational aspects[END_REF]. Furthermore, we handle random particulate systems which, combined with a fine discretization, allow us to address anisotropy issues. Using this paradigm, an equivalent continuous domain is established for which mechanical properties related to the beam elements are determined using a calibration process [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF][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] Hassan | On the suitability of a 3D discrete element method to model the composite damage induced by thermal expansion mismatch[END_REF]. It should be noted that the model proposed in this work is applicable for aligned short-fibre reinforced composites, with the desired geometric parameters and fibre volume fraction. Effective elastic properties of PA6/GF30 are evaluated and compared to analytical, numerical predictions and experimental data as a validation means. The interfacial debonding of SFRC during bending tests and under mixed mode is well simulated. Thus, tensile and shear behaviours on failure mechanisms of PA6/GF30 are highlighted. Results show that the proposed DEM model allows a realistic modelling of interfacial debonding and provides useful information, in particular in terms of energy and force/displacement curve, at the fibre/matrix interface scale (Figure 1). Following these guidelines, the present work is organised as follows. In section 2, we characterise experimentally the macroscopic elastic behaviour of PA6/GF30 while taking into account the influence of different geometric parameters such as fibre dispersion and shape factor. In addition, damage mechanisms are highlighted by Scanning Electron Microscope (SEM) observations. In section 3, a short description of the cohesive beam model is given. Then, a suitable 3D DEM model of SFRC is set up, based on the determination of REV. The later is validated by the estimation of the effective elastic properties of PA6/GF30 in comparison to those determined by micromechanical approaches, numerical and experimental methods. In section 4, a discrete CZM is explored through DCB, ELS and MMB delamination tests. Comparisons between numerical and analytical solutions are performed for validation purposes. In section 5, the case of PA6/GF30 material is investigated. Two models are evaluated, the case of a composite made of a UD single fibre and the case of a multi-fibre composite described by the REV discussed in section 3. Thus, three point bending tests are established and the debonding process is investigated in the final step.

Experimental procedure and theoretical basis

The present section is dedicated to experimentally characterising the PA6/GF30. Firstly, we present the manufacturing process. Then, we determine some morphological and mechanical properties. For this purpose, Digital Image Correlation (DIC) and monotonic tensile tests are used to determine the Poisson's ratio and the Young's modulus of PA6/GF30 respectively. Thereafter, via those 2 tests we characterise elastic properties of PA6/GF30 which is subsequently used to validate the proposed numerical model of PA6/GF30. Finally, we study the damage mechanisms of PA6/GF30 at the microscopic scale using SEM observations of fracture surfaces after 3-point bending test.

Manufacturing process

A composite reinforced with 30% by mass of glass fibre (PA6/GF30) is considered in this study. Commercially, this material is known as ALYAMID C2122 GF30. It is made of a polyamide 6 provided by RADICI Chemicals, reinforced with E-type glass fibre from CIPC using an Amino Silane coupling agent (Table 1). The injection process is performed in accordance to ISO 527-2 type 1A. Before testing, samples were dried at 70

Materials characterization

In order to model composite material, we need to identify the volume fraction of fibres instead of their mass fraction. This latter is determined by pyrolysis technical which remains the most practicable one. Three samples was considered. The average mass fraction is about 29.35% which is in adequacy with the one provided by the technical data sheet (30%). Based on [START_REF] Ksouri | Effects of processing steps and hygrothermal ageing on mechanical performance of PA6GF30 composite: Interfacial shear strength[END_REF], fibre volume fraction V f is set to 15.7%. Glass fibre length is characterised by extraction of fibres from the matrix. The methodology consists in putting fibres on glass slide and making subsequent optical micro observations. Based on image analysis of a set of 500 fibres, fibre length distribution is established [START_REF] Ksouri | Effects of processing steps and hygrothermal ageing on mechanical performance of PA6GF30 composite: Interfacial shear strength[END_REF] as depicted in Figure 2. The number-average fibre length L n is determined to be equal to 183.4 micrometers. Through this preliminary characterization, it is possible to identify the Shape Factor (SF) based on L n with respect to fibre diameter D as:

SF = L n D (1)
Experimental measurements show that SF is close to 18. Moreover, fibre orientation was investigated in [START_REF] Ksouri | Effects of processing steps and hygrothermal ageing on mechanical performance of PA6GF30 composite: Interfacial shear strength[END_REF] using optical microscopy observation and performed on a cross section, cut perpendicular to the injection flow direction, of PA6/GF30 samples. Microstructures show that the majority of fibres present a quasi-circular shape which means that fibres are preferentially oriented in the injection direction. Indeed, short fibre orientation is difficult to quantify and requires advanced techniques to control it. Based on this assumption, we consider that all fibres are aligned.

Mechanical properties

In order to evaluate the longitudinal Young's modulus, 8 specimens were subjected to a tensile test with a monotonic crosshead spead of 5mm/min using a testing machine of type WDW-E (Figure 3). The mean Young's modulus was measured to be 7.363 GPa with a standard deviation of 0.284 GPa. In addition, the stress at rupture was also measured, averaging at 132.687 MPa with a standard deviation of 1.573 MPa. Poisson's ratio of PA6/GF30 is evaluated by Digital Image Correlation (DIC). Hence, it is possible to quantify the displacement and deformation over the entire two-dimensional surface of the specimen. Four specimens were tested according to the ISO 527 procedure at a constant velocity of 5mm/min. Experimental outcomes provide a mean value of Poisson's ratio of 0.358 with a standard deviation of 0.005. broken and the fibre is completely detached and the matrix cannot withstand the stresses transferred by the damage fibre. As a result, it cracks and the material is broken macroscopically. Fibre failure is hardly ever observed for this material. These conclusions lead us to neglect this phenomenon in the modelling. Hence, we focus only on the damage of the fibre-matrix interface.

Cohesive beam model

In this section we present a brief description of the cohesive beam model. In this approach, an Equivalent Continuous Domain (ECD) [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF] is modelled by a particulate system composed of spheres in point contact in which cohesion is introduced using beam elements at the scale of elementary contact. A calibration process is then set up to relate targeted macroscopic elastic properties to microscopic ones associated with beam elements. In adition, a methodology is discussed to investigate the case of a heterogeneous material through the case of PA6/GF30.

ECD

The concept of ECD is based on the discretization of a continuous medium by a particulate system composed of 3D spheres in contact. In the present contribution, particulate systems are generated using the process described by Leclerc [47] which takes advantage of a specific "enlargement-displacement" procedure. It should be noted that the present approach allows for a precise control of intrinsic parameters of the domain, in particular the compacity defined as the volume fraction filled by spheres in a given domain, the coordination number which is the average number of links per particle and the potential polydispersity of spheres. Under the hypothesis that the particulate system respects the assumptions of "Random Close Packing" (RCP) [START_REF] Donev | Improving the density of jammed disordered packings using ellipsoids[END_REF], we consider the ECD as a particulate system representative of the continuous medium. By default, the compacity of the domain is set to 64% and the coordination number to 6.2 for monodisperse spheres. However, to avoid undesirable effects as suggested by Leclerc [47] a sensitivity parameter ε is introduced to increase the area of interaction between DE and meet a higher coordination number of 8.2. Thus, let L ij be the interparticle distance between two particles i and j of radii R i and R j respectively, we assume that contact occurs when the following relationship is verified:

L ij ≤ (1 + ε)(R i + R j ) (2)
Furthermore, to avoid undesirable directional effects and to ensure the isotropy of the system, a slight polydispersity based on a Gaussian model with a coefficient of variation of 0.3 is also introduced [START_REF] Leclerc | Effect of packing characteristics on the simulation of elasticity and brittle fracture by the Cohesive Discrete Element Method[END_REF].

Cohesive beam model

The numerical modelling by DEM is based on the cohesive beam model introduced by André et al. [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. In such a paradigm, the cohesive bond between two particles is modelled by a Euler-Bernoulli beam element [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF] which is described by a set of intrinsic parameters, namely the beam length L µ = (R i + R j ), the microscopic Young's modulus E µ , a circular cross-section of planar quadratic moment I µ and a radius a µ (Figure 6). Nevertheless, for the sake of simplification, a dimensionless beam radius r µ is introduced as a function of R i and R j radii of the particles i and j in contact. a µ and I µ depend on r µ and are expressed as follows:

a µ = r µ R i + R j 2 (3)
I µ = π a 4 µ 4 = π r 4 µ 64 (R i + R j ) 4 (4)
From a practical point of view, beam reaction forces are controlled by an internal 6-component generalised force vector. This includes normal and tangential components able to counteract the relative displacement of particles, and moment ones resisting bending and torsion effects. For information purposes, the expression of every component depends on E µ and r µ microscopic parameters and derives from a 6-by-6 matrix system given in previous work [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] Hassan | On the suitability of a 3D discrete element method to model the composite damage induced by thermal expansion mismatch[END_REF]. Furthermore, the beam length L µ has no effect on the macroscopic behaviour for a particulate system dense enough to meet the RCP assumptions [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]. Internal cohesive forces between particles i, j are given by the classic Euler-Bernoulli theory. Newton's 2 nd law is applied to each particle in both translation and rotation. The motion of the particle is thus computed using an explicit temporal scheme of Velocity-Verlet type [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]. This choice is mainly motivated by the simplicity of the scheme and the applicability of the explicit scheme for large-scale DEM calculations. It should be noted that a Rayleigh-type damping is also added to internal forces to reduce dynamic effects over time.

Macroscopic elastic properties

The microscopic Young's modulus E µ does not correspond to its macroscopic equivalent E M and the Poisson's ratio ν M is a priori unknown. Therefore, a calibration process must be set up in order to identify local parameters leading to the expected macroscopic properties E M and ν M . For this purpose, specific tensile tests are carried out for a wide range of configurations as described in [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]. Polynomial functions are then established using regression analysis to correlate microscopic parameters E µ and r µ to their corresponding macroscopic properties E M and ν M as follows:

r µ = P (ν M ) (5)
E µ = Q(r µ )E M (6)
Thus, knowing ν M of a given material, we can obtain its corresponding r µ , and E µ can then be deduced from E M and r µ . Macroscopic properties of PA6 and GF30 are given in Table 1. Based on the calibration process, the corresponding microscopic parameters are E P A6 µ = 1738 GPa, r P A6 µ = 0.05 for the PA6 case and E GF 30 µ = 725 GPa, r GF 30 µ = 0.39 for the GF30 case.

DEM model of PA6/GF30

A suitable DEM model of PA6/GF30 is now set up using a specific approach. This is based on the generation of cubic REV of the targeted material the definition of which depends on some assumptions in size established from Fast Fourier Transform (FFT) calculations. The use of this approach is justified by its efficiency in estimating effective properties of heterogeneous media and provide information on the required REV size. In this work, we benefit from a FFT approach based on the Eyre-Milton scheme for the resolution [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF]. It should be noted that, according to the experimental characterization discussed in Section 2, we consider randomly distributed fibres described by cylinders with SF=18 and V f =15.7%. The proposed procedure is as follows (Figure 7). In a first preliminary step, the optimal number of voxels to be used in FFT calculations is investigated for a set of scale ratio, which is the ratio between the edge length of the REV and the length of the fibre. Thus, we study the influence of the number of voxels along each direction on the evolution of the effective longitudinal Young's modulus while considering a maximum scaling ratio of 8. For information purposes, the number of voxels is varied along each direction from 50 to 700 voxels and 3 different elementary volumes are handled for each case to ensure the reproducibility of the results. From Figure 8(a), we can conclude that a 500×500×500 resolution is accurate enough to ensure suitable calculations. As a second step, we aim to determine the optimal scale ratio (SR). Therefore, we evaluate the effect of this parameter on the effective longitudinal Young's modulus. Predictions are carried out using a 500×500×500 resolution, as previously discussed, and for several SR between 1.2 to 8. Results are illustrated in Figure 8(b) and from a quantitative standpoint, we can deduce that an optimal SR equal to 1.5 suffices to appropriately predict the effective Young's modulus of PA6/GF30. In a third and last step, the density of particles to be used in DEM calculations is estimated. For this purpose, based on a set of 4 cubic REV with a SR set to 1.5, V f is numerically determined for a large range of number of DE in the interval [8 × 10 5 , 1.6 × 10 7 ] (Figure 9). Results show a two-step convergence. For a number of particles less than 8×10 6 DE, V f increases gradually. However, for a number of particles higher than 8×10 6 DE, V f begins to stabilise around 15.7%. It should be borne in mind that due to the random nature of fibre generation, the convergence curve is not sufficiently smooth. In addition, the relative difference between the experimental and numerical volume fraction is less than 1% which exhibits a quite suitable representativity of the heterogeneous continuous domain. As a result, we can consider that 8 million DE is a reasonable choice to be used in the present context. Thus, Figure 10 (a) depicts an example of REV generated using a SR of 1.5, composed of 424 aligned fibres, SF=18 and V f =15.7%. Since all particles are randomly distributed with respect to assumptions of the RCP, each DE is associated with a given phase according to its distance "d" between the axis of the fibre and the center of the particle (Figure 10 (b)). Thus, this procedure allows to generate cylindrical fibres with particle-scale irregularities at the interface according to defined geometric parameters. Microscopic parameters E µ and r µ for both materials PA6 and GF30 are provided in the previous subsection. Since we are dealing with a heterogeneous medium, the transition from one phase to another should be as smooth as possible. Thus, microscopic interfacial parameters, corresponding to the beam elements linking every pair of particles at the interface (E int µ and r int µ ), were defined as the arithmetic average of their counterparts: Results provided in Table 2 show that the prediction given by FFT is quite similar to those estimated by TW or MT approaches with relative differences less than 2% for the Young's modulus, and 3% for the Poisson's ratio whatever the theoretical reference. This highlights the relevance of the proposed definition of REV in the case of SFRC PA6/GF30. Nevertheless, relative differences between DE results and FFT outcomes are estimated at 6.3% and 12.1% for Young's modulus and Poisson's ratio respectively. 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] showed that increasing the coordination number from 4.5 to 7.5, for a compacity set to 64%, leads to a decrease in the maximum admissible Poisson's ratio. In our case, a coordination number of 8.2 was used corresponding to a maximum Poisson's ratio of 0.28, which is lower than the targeted value of PA6/GF30 and could, at least partially, explain the difference in Poisson's ratio between DEM and FFT approaches. Furthermore, differences between FFT and DE results could also be related to an underestimation of the required number of particles to ensure the convergence of DEM predictions to FFT ones. Moreover, relative differences between DE results and experimental outcomes are 19.7% and 22% for Young's modulus and Poisson's ratio respectively. Discrepancies between DE predictions and experimental measurements can be explained by the level of detail of the PA6/GF30's microstructure and the fact that the proposed numerical model does not take into account some parameters such as interfacial effects, fibre orientation distribution or even local defects. However, DE calculations show promising results in terms of elastic properties in comparison to experimental measures, FFT and analytical approaches. Thus, we can conclude that the proposed DEM model allows for a suitable representation of PA6/GF30 in terms of effective elastic properties.

E int µ = E P A6 µ + E GF 30 µ 2 (7)
r int µ = r P A6 µ + r GF 30 µ 2 (8)

Approach

DEM modelling of interfacial debonding

One of the crucial micromechanical damage mechanisms in composites is fibre/matrix interfacial debonding.

Owing to the mismatch of elastic properties at the interface between two different materials, authors confirm that the mixed mode will always persist [START_REF] Karlsson | Obtaining mode mixity for a bimaterial interface crack using the virtual crack closure technique[END_REF]. Modelling this phenomenon requires a robust and efficient approach. DEM coupled with energy-based criterion may be a solution to more properly describe the interfacial degradation mechanism. Thus, interfacial debonding is first studied through standard delamination tests for DEM model validation purpose. Double Cantilever Beam (DCB), End Loaded Split (ELS) and Mixed Mode Bending (MMB) tests are investigated by DEM to evaluate the ability of the interfacial model to reproduce the delamination process in mode I, mode II and mixed mode (I+II), respectively.

DCB test

In a first step, we consider the standard DCB test to reproduce the mode I of delamination of a laminate beam with following dimensions [START_REF] Le | Modeling damages and cracks growth in composite with a 3D discrete element method[END_REF]: a length l = 45 mm, a width w = 6 mm and a thickness 2e = 3 mm (Figure 11(a)). All links between particles are cohesive beam elements except that at the interface located at the mid-thickness. In this specific domain, a pre-crack of length a 0 = 13.5 mm is defined by deleting interfacial contacts linking two particles belonging to different halves of the beam in this area. The remaining interfacial links are replaced by elements wih normal stiffness K n . For information purposes, DE calculations are performed using an adjustable time step which is determined as a function of critical time step ∆t c crit associated to each beam element as follows :

∆t = C t min c∈ξ (∆t c crit ) (9)
where ξ designates the set of contacts within the granular packing and C t is a security ratio. In the community, the latter is chosen between 1 2π [52] and 0.5 [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]. In our case C t is set to 0.3. In addition, Rayleigh damping which depends on a constant β coefficient is introduced in DEM simulations to reduce dynamic effects and ensure the convergence of numerical outcomes [START_REF] Leclerc | Discrete element method to simulate the elastic behavior of 3D heterogeneous continuous media[END_REF]. Thus, the time step is finally expressed as follows:

∆t =    ∆t β β > 1 ∆t β ≤ 1 (10)
Nevertheless, it should be remembered that the introduction of damping imposes the use of a rather small time step in discrete simulations. As a result, a compromise has therefore to be found to reduce those dynamic effects while avoiding a prohibitive calculation cost. Thus, the time step of DCB test is set to 2.4×10 -6 s with β = 15. In order to ensure the propagation of the crack in mode I, a bilinear softening model is defined at the level of interfacial links. This model presents three phases (Figure 11(b)). A first loading phase corresponds to a constant initial normal stiffness K 0 n up to a critical normal displacement U c n which reads as follows:

U c n = N K 0 n ( 11
)
where N is the critical normal force. Then, the interface starts to degrade until the normal displacement U n reaches a maximum value U m n which is defined as follows:

U m n = 2AG Ic N ( 12
)
where A is the average cross-section of cohesive beam elements which is also associated with interfacial links in the present study, and G Ic is the mode-I critical energy release rate of the investigated material. Finally, at this step, a crack initiates and propagates along the interface with a mode I propagation. Damage is introduced at the level of the normal stiffness related to interfacial links via a damage parameter D, the value of which depends on U n . The normal stiffness is then defined according to the following equation:

K n = (1 -D)K 0 n ( 13
)
Three phases can be identified: in the first phase, the interfacial behaviour is assumed to be elastic without damage and D is zero. In a second phase, in which U n >U c n , a linear softening is introduced via D which is then strictly positive and evolves according to U n . Finally, in a third phase, D is equal to 1 and the cohesive link is assumed to be broken. The expression of the damage parameter D is given below as a function of

D =          0 U n ≤ U c n U m n (Un-U c n ) Un(U m n -U c n ) U c n < U n < U m n 1 U n ≥ U m n ( 14
)
DEM calculations are performed using a constant velocity V = 5 mm/min at both sides of the left edge of the sample. The mechanical properties used for this study are [START_REF] Le | Modeling damages and cracks growth in composite with a 3D discrete element method[END_REF]: a Young's modulus E = 135 GPa, a Poisson's ratio ν 22 = 0.24, G Ic = 0.56 N/mm and an interfacial strength defined as the ratio of normal critical force (N ) to the average cross-section of cohesive beam elements (A) and set to σ max = 3 × 10 3

MPa. Nevertheless, one has to be borne in mind that the interfacial strength is defined at the contact level and has no macroscopic meaning. The particulate system is composed of 400 000 DE which corresponds to a discretisation level, defined as the ratio of the smallest dimension of the modelled domain to the average diameter of particles, close to 22 DE. For validation purposes, a comparison is carried out with an analytical solution which is derived from elastic beam theory [START_REF] Le | Modeling damages and cracks growth in composite with a 3D discrete element method[END_REF]. From Figures 12(a) and 12(b) it can be seen that macroscopic responses in terms of force/displacement and crack length-displacement curves are close to analytical expectations. As a first step, the displacement linearly evolves without any damage until ∆ = 0.23 mm. Then, in a second step, the crack initiates and propagates along the interface between the two layers leading to a decrease of the force (Figure 13). The fluctuations occuring during the failure step are related to explicit dynamic simulations and require a high damping and fine time steps to be deleted.

ELS test

ELS test is now carried out using a time step of 6.6×10 -5 s with β = 10 and a constant velocity V = 5 mm/min applied on the left edge of a sample composed of 198 000 particles (Figure 14(a)). This corresponds to a discretisation level of 10 DE along the thickness of the specimen. One has to keep in mind that the level of discretisation is different from one test to another in order to ensure a good compromise between the accuracy of numerical outcomes and the computational cost. The present configuration is characterised by the following dimensions [START_REF] De Moura | Equivalent crack based analyses of ENF and ELS tests[END_REF]: a length l = 105 mm, a width w = 20 mm, a thickness 2e = 3.6 mm and a pre-crack a 0 = 60 mm. The mechanical properties are: a Young's modulus E = 150 GPa, a shear modulus G 12 = 6 GPa, an interfacial strength τ max = 2.5 × 10 3 MPa and a mode-II critical energy release rate G IIc = 1.5 N/mm. All DE are connected using cohesive beam elements except for the interface. In this specific area, spring elements with tangential stiffness K s are introduced with the exception of a pre-crack area. Spring elements are governed by a shear/displacement curve using a bilinear softening model (Figure

U c s = S K 0 s ( 15
)
where S is the tangential critical force and K 0 s is the initial tangential stiffness. Then, in a second stage, the degradation of the interface at the cohesive beam element level occurs until the tangential displacement U s reaches the maximum displacement U m s which reads:

U m s = 2AG IIc S (16)
Lastly, interfacial bonds are assumed to be broken so that a crack initiates and propagates along the interface.

A damage parameter D is then introduced as a function of U s such as:

D =          0 U s ≤ U c s U m s (Us-U c s ) Us(U m s -U c s ) U c s < U s < U m s 1 U s ≥ U m s ( 17
)
The tangential stiffness K s is then computed according to the following expression:

K s = (1 -D)K 0 s ( 18
)
To verify the validity of this study, DE results are compared to an analytical solution given in [START_REF] Ilyas | Damage modeling of carbon epoxy laminated composites submitted to impact loading[END_REF]. The force as function of displacement and crack length with respect to the displacement curves are plotted in Figures 15(a) and 15(b) respectively. Results show that the force as function of displacement curve includes two stages: a linear step where the force increases without any damage until a peak is reached which corresponds to ∆ = 13.2mm. In the second part, the force gradually drops due to a mode II delamination effect. The crack propagates between the two layers along the interface until the failure of the material (Figure 16). The presence of slight differences between DE results and the analytical solution should be noted, in which some fluctuations appear during the damage process (Figure 15(a)). As discussed before, this can be explained by dynamic effects occuring during delamination. In the case of the crack length, a mismatch between the two plotted curves can be observed during the failure process (Figure 15(b)). We hypothesize that this discrepancy could be partially explained by the length of the pre-crack which is quite long with respect to the length of the pattern [START_REF] Blackman | Mode ii fracture testing of composites: a new look at an old problem[END_REF]. Moreover, the value of the Rayleigh damping coefficient used in DEM simulations could affect the onset and propagation of crack growth. Note that the ELS test is performed with a reduced level of discretisation compared to the DCB test, which may also influence the accuracy of numerical results.

MMB test

We now consider the mixed-mode fracture through MMB test the configuration of which is described in Figure 17(a). The following dimensions and material properties are used in this study [START_REF] Xie | Discrete cohesive zone model for mixed-mode fracture using finite element analysis[END_REF]: a length l = 50 mm, a width w = 10 mm, a thickness 2e = 3 mm, a parameter to adjust the ratio between the two applied forces (F 1 and F 2 ) c = 50 mm, a pre-crack length a 0 = 30 mm, the Young's modulus E = 135 GPa, the Poisson's ratio ν = 0.24, and mode-I and mode-II critical energy release rates G Ic = G IIc = 4 N/mm. The test can be conducted by applying a force F by a rigid loading lever with adjustable arm. We derive the expressions of resulting F 1 and F 2 loads from the static balance of force and moment as follows:

F 1 = F c l ( 19
)
F 2 = F c + l l ( 20
)
To perform the DE analysis, we consider a sample composed of 237 000 DE for a discretisation level of 12 approximately. We consider the geometry and boundary conditions presented in Figure 17(a) in which displacement control is used for the DE calculations. In order to respect the ratio between forces F2 F1 = c+l l , velocity control is defined as follows V2 V1 = 0.027 where V 1 = 5 mm/min. We also set the time step of DE calculations to 9.1×10 -6 s with β = 15. Spring elements in normal and tangential directions are inserted in the interfacial zone except in the pre-crack to detect the debonding phenomenon. We consider a damage model based on a bilinear softening law taking into account mixed-mode shear and tensile/compressive loadings [START_REF] Harper | Cohesive zone length in numerical simulations of composite delamination[END_REF]. The latter is studied as a function of a set of physical parameters such as G Ic and G IIc , and tensile N and shear S strength limits associated with the cohesive bond (Figure 17(b)). In this model, normal (U n ) and tangential (U s ) displacements are input parameters obtained numerically at a previous time step. Mode-I and mode-II bilinear models depend on critical displacements, respectively U c n and U c s , which are expressed by Eq.11 and Eq.15. Both bilinear models also depend on maximum displacements, respectively U m n and U m s in mode I and mode II as mentioned in Eq.12 and Eq.16. Let the total effective displacement U e be defined as a function of U n and U s as follows:

U e = U 2 n + U 2 s ( 21
)
It is possible to define a critical effective displacement U c e beyond which the cohesive bond is assumed to be damaged. In order to define such a parameter as a function of U c n and U c s , we have to introduce a coupling parameter η according to the following expression: The following relationship between critical displacements U c n and U c s , and displacements U n and U s is also assumed:

η =    Us Un U n >0 0 U n ≤ 0 (22)
U n U c n 2 + U s U c s 2 = 1 (23)
By injecting Eq.23 in Eq.21 and taking into account the coupling parameter η, we obtain the following expression of U c e :

U c e =    U c n U c s 1+η 2 (U c s ) 2 +(ηU c n ) 2 U n >0 U c s U n ≤ 0 (24)
In order to define a maximum effective displacement U m e , the strain energy release rates in mode I (G I ) and in mode II (G II ) are introduced according to the following expressions:

     G I = K 0 n U m e U c e 2A(1+η 2 ) G II = K 0 s η 2 U m e U c e 2A(1+η 2 ) (25)
The most common law to describe the phenomenon of mixed-mode delamination is the power law [START_REF] Harper | Cohesive zone length in numerical simulations of composite delamination[END_REF], which relates G I and G II to G Ic and G IIc :

G I G Ic α + G II G IIc α = 1 (26)
The energy criterion can follow either a linear (α=1) or a quadratic (α=2) criterion which are the most common. For information purposes, in the present work, α is set to 1 as typically used in the literature [START_REF] De Moura | Equivalent crack based analyses of ENF and ELS tests[END_REF][START_REF] De Lorenzis | Coupled mixed-mode cohesive zone modeling of interfacial debonding in simply supported plated beams[END_REF][START_REF] Chaves | Fracture mechanics tests in adhesively bonded joints: a literature review[END_REF]. By injecting Eq.26 into Eq.25, we then obtain the expression of U m e as a function of

U m n , U m s , U c e , K 0 n , K 0 s , G Ic , G IIc , η and α: U m e =      2A(1+η 2 ) U c e K 0 n G Ic α + η 2 K 0 s G IIc α -1 α U n >0 U m s U n ≤ 0 (27)
In addition, the constitutive equations for mixed-mode fracture depend on a damage parameter D which reads as follows:

D =          0 U e ≤ U c e U m e (Ue-U c e ) Ue(U m e -U c e ) U c e < U e < U m e 1 U e ≥ U m e ( 28
)
and expressions of K n and K s as a function of D, K 0 n and K 0 s are given in Eq.14 and Eq.18 respectively. For DE calculations, based on preliminary studies, tensile σ max and shear τ max interfacial strengths are respectively set to 1.4 × 10 4 MPa and 3.5 × 10 3 MPa. Figure 18 shows DEM predictions in terms of force/displacement and crack length-displacement curves compared to a reference analytical solution. Details of the analytical solution are given in [START_REF] Mi | Progressive delamination using interface elements[END_REF]. From Figure 18(a) one can observe that DEM and analytical This highlights the ability of DEM to predict the mixed-mode fracture during a delamination process (Figure 19).

Application to PA6/GF30

Based on SEM observations, tensile and shear stresses develop at the fibre/matrix interface of PA6/GF30 during a 3-point bending test. This mechanism is characterised by a mixed-mode debonding, which potentially reduces or even breaks the adhesion between the fibre and the matrix. The present section investigates these phenomena through an application to PA6/GF30 of the mixed-mode interfacial debonding model discussed in section 4. Two configurations are studied: the case of a UD fibre composite and the case of an behaviour was observed by Garcia et al. [START_REF] García | Fiber-size effects on the onset of fiber-matrix debonding under transverse tension: a comparison between cohesive zone and finite fracture mechanics models[END_REF] in the case of fibre/matrix debonding under transverse tension.

The authors related the post-peak nonlinearity to the evolution of interface decohesion. They also mentioned that fibre size in terms of fibre radius as well as fibre volume fraction affect the nonlinear force/displacement response. Figure 21(b) illustrates the evolution of the rate of broken bonds at the interface as function of the displacement in Z direction as well as debonding area for a cutting plane at X = 0.5 mm. Results show that damage appears for a displacement equal to 0.201 mm which is quite close to the value associated to the local force peak at the end of the elastic phase. Furthermore, the larger the displacement, the higher the number of broken bonds at the interface. Thus, the failure process evolves until 58.482% of interfacial bonds are broken, corresponding to a displacement of 0.800 mm. DE calculations are limited to 80% of the displacement in relation to the length of the studied system. The choice of this criterion is mainly justified by the necessity to remain within a real scale of displacement. Note that to ensure the complete detachment of the fibre from the matrix, we will be confronted with dramatically dynamic effects in relation to contact management. Furthermore, simulation should be conducted for a larger displacement which does not fit with the model used in the present work.

Multi-fibre composite

To improve the representativeness of the numerical model of PA6/GF30, we consider the case of an aligned short-fibre reinforced composite. As depicted on Figure 22 and based on assumptions discussed in subsection 3.4, a cubic REV (1mm×1mm×1mm) is generated from which a parallelepipedic elementary sub-volume (0.1mm×1mm×0.1mm) is extracted. Nevertheless, one must keep in mind that the SR related to DE model has been evolved from 1.5 to 3 to increase the number of fibres along their longitudinal axis in the elementary sub-volume and meet the hypotheses of representativity. A three point bending test is carried out along Z Fig. 22: Configuration and boundary conditions associated to multi-fibre composite direction using a sample composed of 400 000 DE and boundary conditions defined on Figure 22. For information purposes, the number of DE is purposely set to 400 000 particles so as to ensure a good compromise between the accuracy of numerical outcomes and the computational cost. DE calculations are performed using the same material properties as previously defined in the case of the single-fibre model, except that τ max and σ max are increased to 4.992 × 10 4 MPa and 1.248 × 10 4 MPa respectively to detect the different stages of damage through a wide range of displacement. Figure 23(a) illustrates the evolution of the force as function of the displacement in Z direction. Numerical outcomes reveal that as observed in the case of the single-fibre model, the force develops linearly with some premature debonding effects on the edges of the sample until it reaches 0.502 N for a displacement corresponding to 0.045 mm. At that point, the rate of broken bonds at interface starts to increase (Figure 23(b)) and the larger the displacement the higher the rate of broken bonds. The interfacial debonding evolves progressively from both ends of the sample until it reaches its middle plane for 49.89% of broken bonds at the interface. At a displacement of 0.055 mm which corresponds to a force of 0.528 N, a stage of saturation arises where the matrix is no longer able to transmit strengths to the reinforcement. This phenomenon is due to damage induced by the multiplication of fibre detachment. Hence, contrary to the single-fibre case, the force continues to decrease without regaining strength until material failure. Furthermore, one can notice that the level of force detected in the single-fibre case is different from what we obtain in the multi-fibre case. We assume that this difference is probably connected to the variation in the studied configurations as well as dimensions of samples used in each test. In addition, This difference can be explained by the influence of the adhesion mechanisms at the fibre/matrix interface; a phenomenon of mechanical attachment at the interface is more dominant in the case of single fibre composite. In this case the fibre/matrix interface is initially distributed over the entire length of the sample. This gives a potential to withstand to the interfacial debonding process. Thus, intact areas of the interface continue to transmit force from the matrix to the reinforcement during decohesion. This may justify the force recovery observed in the case of single fibre composite. On the other hand, interfacial zones in the case of multi-fibre composite are linked to lengths of short fibres, i.e. they are limited in size compared to the length of the sample. This reduces the ability of the intact interfacial links during the damage process to ensure the transfer of force between matrix and fibres, which clearly demonstrates the absence of rising force in the case of multi-fibre composite. Based on these results, we can conclude that the developed DEM model is quite suitable to simulate the debonding process at the fibre/matrix interface of PA6/GF30 during a 3-point bending test and opens a wide variety of perspectives which are suggested later.

Conclusions and prospects

The present contribution was devoted to highlighting the potential of the DEM to deal with interfacial problems in heterogeneous continuous domain especially with regard to SFRC PA6/GF30. In a first step, elastic properties of PA6/GF30 were determined using experimental campaign. Based on SEM observations on fracture surfaces after three point bending test, it turned out that interfacial debonding constitutes the main mechanism of failure in PA6/GF30. In a second step, a DE model was introduced in the context of short glass fibres. Numerical results, in terms of effective elastic properties, are in good agreement with micromechanical models as well as experimental data. Thus, we can conclude that the 3D DEM model is quite suitable for modelling SFRC and simulating their elastic behaviour, which constitues one of the originalities of this contribution. In a third step, an interfacial 3D DEM model was discussed and validated in the case of composite delamination tests in mode I, mode II and mixed mode. Suitable results were obtained by comparison with analytical solutions. Hence, DEM revealed its capacity to deal correctly with interfacial problems. Moreover, results show that unlike analytical approaches which are too perfect, DEM is able to detect some effects related to unstable propagation aspects of delamination at the interface scale. In a final step, debonding process in PA6/GF30 under bending load was investigated through two configurations: a single UD glass fibre embedded in PA6 matrix and an aligned short-fibre model. Some differences in term of the force/displacement curve between the two configurations were highlighted and discussed. In the near future, we aim to integrate hygrothermal parameters into the developed interfacial model to simulate the behaviour and the damage of the PA6/GF30 under a wide range of environmental conditions.

Fig. 1 :

 1 Fig. 1: Flowchart of the proposed 3D DEM approach to simulate interfacial debonding in 3D SFRC

Fig. 2 :

 2 Fig. 2: Histogram of glass fibre for PA6/GF30 after injection molding process

Fig. 3 :Fig. 4 :Fig. 5 :

 345 Fig. 3: (a) Stress-strain curves of 8 specimens of PA6/GF30, (b) Tensile test conditioning

Fig. 6 :

 6 Fig. 6: Cohesive beam model

Fig. 7 :

 7 Fig. 7: Methodology proposed to model PA6/GF30 using DEM and REV concept

Fig. 8 :Fig. 9 :

 89 Fig. 8: Effect of (a) voxel number along each direction and (b) SR on the effective longitudinal Young's modulus

Fig. 10 :

 10 Fig. 10: (a) REV of PA6/GF30 short fibre composite (SF=18 and SR=1.5), (b) Magnified view to highlight how each DE is associated to a given phase: PA6 or GF30

Fig. 11 :

 11 Fig. 11: (a) DCB test , (b) Mode-I bilinear softening law

Fig. 12 :Fig. 13 :Fig. 14 :

 121314 Fig. 12: (a) Force/displacement curves , (b) Crack length-displacement curves in the DCB test

Fig. 15 :Fig. 16 :

 1516 Fig. 15: (a) Force/displacement curves, (b) Crack length/displacement curves in the ELS test

Fig. 17 :

 17 Fig. 17: (a) Mixed-mode MMB test , (b) Mixed-mode bilinear softening law

Fig. 18 :

 18 Fig. 18: (a) Force/displacement curves , (b) Crack length/displacement curves in the MMB testapproaches lead to very close force/displacement curves except for some fluctuations in DE calculations which are related to dynamic effects. It can be seen in Figure18(b) that, in comparison with the analytical solution, interfacial bonds start to fail for a slightly lower displacement and the damage propagates a little faster until the final stage of the test but results are globally in quite good agreement. This highlights the ability of DEM to predict the mixed-mode fracture during a delamination process (Figure19).

Fig. 19 :

 19 Fig. 19: Crack propagation at (a) ∆ = 5 mm and (b) ∆ = 7 mm

Fig. 20 :

 20 Fig. 20: (a) Configuration and boundary conditions of the single-fibre composite, (b) Corresponding DE modelconsider microscopic parameters E µ and r µ defined in subsection 3.3 for both materials PA6 and GF30. The mixed-mode cohesive zone model discussed in subsection 4.3 is implemented in the heterogeneous medium between matrix and fibre to simulate the interfacial debonding during a 3-point bending test. Note that in our case, no pre-crack is considered and the fracture is only constrained to the fibre-matrix interface. The bending test is carried out using a symmetrical configuration where only one half of the beam is modelled as shown in Figure20(a). Boundary conditions are as follows: the sample is supported at the right lower edge and the test is conducted by applying an imposed constant velocity equal to 5mm/min along Z direction on the upper edge. Based on preliminary studies, interfacial strengths are set to τ max = 2.496 × 10 4 MPa and σ max = 6.24 × 10 3 MPa. Because of a lack of experimental data, G Ic and G IIc values come from a fibrereinforced epoxy composite and are set to G Ic = G IIc = 0.5 N/mm[START_REF] Le | Modeling damages and cracks growth in composite with a 3D discrete element method[END_REF]. Figure21(a) represents the force as function of the displacement in Z direction. Results exhibit that the force evolves in a linear way reflecting the elastic behaviour of the material until it reaches a value of 7.684 N for a displacement corresponding to 0.238 mm. From this value, a softening of the behaviour appears, indicating the beginning of the damage process. This stage is characterised by the appearance of debonding at the top and bottom of the fibre on the right edge of the sample. Afterwards, the debonding evolves longitudinally along the contact area between the fibre and the matrix. From this stage, there is a change in slope for a force of 8 N which corresponds to a displacement of 0.419 mm. The transmission of force from the matrix to the fibre is then limited to a smaller area governed by intact interfacial bonds which require more energy to break. Thus, the force rises again up to a value of 10.910 N for which the fibre debonding process almost reaches the left end of the sample. This

Fig. 21 :

 21 Fig. 21: (a) Force/displacement curve and debonding patterns of the single-fibre composite, (b) Crack length/displacement curve and debonding area of the single-fibre composite

Fig. 23 :

 23 Fig. 23: (a) Force/displacement curve and debonding patterns of the multi-fibre composite, (b) Crack length/displacement curve and debonding area (cutting plane Y = 1 mm) of the multi-fibre composite

  • C for 24 hours and immediately stored in vacuum desiccators.

	Material	Specifications	Supplier
	PA6	Density=1.13	RADICCI
		Young's modulus=2.03 GPa
		Poisson's ratio=0.33
	E-type glass fibre	Density=2.6	CEPIC
		D=10 µm -L=4.5mm
		Young's modulus=72.5 Pa
		Poisson's ratio=0.22

Tab. 1: Material properties

Acknowledgements

The authors would like to thank the UPJV Plateforme Microscopie Électronique team especially M. Arash Jamili for his support and the help provided for this research work.

Compliance with ethical standards

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.