Margaux Sage
email: margaux.sage@u-bordeaux.fr

Jérémie Girardot

Jean-Benoît Kopp

Stéphane Morel

A damaging beam-lattice model for quasi-brittle fracture

Keywords: Quasi-brittle, Discrete Element Method, Damage, mixed mode, energetic criterion

  L'archive ouverte pluridisciplinaire

Introduction

The modeling of quasi-brittle fracture has been widely studied. Of the many approaches considered, the finite element method is the most extensively used.

The most popular finite element-based approaches are eXtended-Finite Element Method (X-FEM), damage models, cohesive zone models, and more recently 5 the phase-field model. Extended finite element method, developed in [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF], [START_REF] Moës | A finite element method for crack 780 growth without remeshing[END_REF], [START_REF] Moës | X-fem, de nouvelles frontières pour les éléments finis[END_REF], simulates cracks within the finite element mesh using functions that enrich the degrees of freedom of the mesh nodes. This model performs well when simulating crack propagation for complex cracking paths. However, it struggles when modeling processes such as fragmentation or multi-cracking, due to a large num-10 ber of cracks present.

Another approach lies in damage models [START_REF] Ragueneau | Damage model for concrete-like 785 materials coupling cracking and friction, contribution towards structural damping: first uniaxial applications[END_REF], [START_REF] Bazant | Finite element modeling of crack band propagation[END_REF], [START_REF] Mazars | The unilateral behaviour of damaged concrete[END_REF], [START_REF] Fichant | Isotropic and anisotropic descriptions of damage in concrete structures[END_REF]. They allow the stiffness of the material to be degraded locally in the meshes where cracks are present.

Initially, faced with mesh dependency problems, work introducing regularisation methods, such as non-local damage models [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] and gradient damage methods 15 [START_REF] Frémond | Damage, gradient of damage and principle of vir-800 tual power[END_REF], made it possible to get rid of mesh dependency and snapback problems. Cohesive zone models are another approach that aims to model crack development zones ahead of cracks or defects [START_REF] Van Den Bosch | An improved description of the exponential xu and needleman cohesive zone law for mixed-mode decohesion[END_REF], [START_REF] Dourado | Comparison of fracture properties of two wood species through cohesive crack simulations[END_REF], [START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF]. This model is applied with the element method by adding cohesive elements between the meshes. Within these elements placed ahead of the crack understudy, an adhesion force is defined 20 that accounts for lower stresses ahead of cracks than those predicted by linear fracture mechanics. Different laws are used to define the forces of the cohesive used, with the result that the stiffness of the structure depends on the initial stiffness within the cohesive elements. More recently, work on the phase-field model [START_REF] Francfort | Revisiting brittle fracture as an energy mini-815 mization problem[END_REF], [START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field fe implementations[END_REF], [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF], [START_REF] Nguyen | A phase field method to simulate crack nucleation and propagation in strongly heterogeneous materials from direct imaging of their 835 microstructure[END_REF] has been carried out. This technique is based on the definition of a phase-field in the system. The presence of a crack is introduced as the interaction between 2 phases. The separation of the material is allowed 30 through the definition of the continuous phase field function which takes the values +1 and -1 on both sides of the crack. This model appears to be particularly effective in modeling phenomena such as branching and coalescence of cracks.

Although efficient in some cases, these models are confronted with the dis-35 advantages of the finite element method, such as the mesh dependency and the high computational costs in the presence of a large number of cracks. The disadvantage of mesh dependence has been partially solved for the damage and phase-field models by introducing an internal length, the physical meaning of which remains to be established. These difficulties in modeling fracture phenom-40 ena with the finite element method lead us to look at other numerical methods such as the discrete element method.

In general, a discrete disordered media is well known for its numerical ability to easily describe fracture mechanisms, and especially for elastic-brittle behaviors [START_REF] Herrmann | 5 -modelization of fracture in disordered sys-840 tems[END_REF][START_REF] Carmona | From fracture to fragmentation: 845 Discrete element modeling[END_REF]. The beam lattice approach in the Discrete Element Method 45 paradigm is one of the various examples to simulate complex crack patterns [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a Journal Pre-proof silica glass with a blunt indenter[END_REF].

The mainly used criterion for beam links consists in a cut-off method [START_REF] Schlangen | New method for simulating fracture using an elastically uniform random geometry lattice[END_REF], [START_REF] Šavija | Cracking of the concrete cover due to reinforcement corrosion: A two-dimensional lattice 860 model study[END_REF], [START_REF] Iordanoff | A discrete element 865 model to investigate sub-surface damage due to surface polishing[END_REF] in which the fracture of the beam is considered when a given threshold is reached.

Among cut-off criteria, some are based on threshold stresses [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a Journal Pre-proof silica glass with a blunt indenter[END_REF], strain threshold defined from axial and rotational deformations [START_REF] Schneider | Modeling of material failure by the 870 discrete element method[END_REF] or using a cohesive model 50 only in the longitudinal direction [START_REF] Mariotti | Modeling of the fragmentation by discrete element method[END_REF]. Improved cut-off criteria can also be used as in [START_REF] Tavarez | Discrete element method for modelling solid and particulate materials[END_REF] where fracture only occurs if a stress threshold is maintained for a while to reach the fracture energy of the material. Similarly, many works, such Among the damaging link models, the one consisting of replacing a beam link by damaging springs (in normal and tangential directions) when the threshold is reached leads to various studies in the literature [START_REF] Leclerc | On a discrete element method J o u r n a l P r e -p r o o f Journal Pre-proof ers Structures[END_REF]. However, the analysis in [START_REF] Schlangen | Fracture simulations of concrete using lattice models: Computational aspects[END_REF] shows that beam links possess a better accuracy on fracture patterns than spring links. An improvement such as the development of a new criterion for 70 beam links built on the spring link work would merely enhance the modeling of fracture behavior. Indeed, most of the papers in the literature, such as detailed thereafter, define damage variable on the normal spring and the differences between the models come from the damage definition.

75

In [START_REF] Gao | A discrete element model for damage and fatigue crack growth of quasi-brittle materials[END_REF], an equivalent deformation and principal stresses are evaluated. These quantities are deduced from mean strains and stresses in the normal and tangential directions expressed with normal and tangential forces and displacements.

On the other side of this mathematical approach, a simpler physical approach [START_REF] Gui | Modelling the dynamic failure of brittle rocks using a hybrid continuumdiscrete element method with a mixed-mode cohesive fracture model[END_REF] is based on cohesive zone models where the linear slope damage law is set by 80 cohesive energy in mode I, energy dissipated by the normal spring, and mixed maximum displacement. The normal and tangential stiffnesses are degraded by this damage. To take into account mixed-mode failure, a Mohr-Coulomb fracture criterion is defined from the mode I and mode II fracture energies.

J o u r n a l P r e -p r o o f

Journal Pre-proof Another work inspired by cohesive zone models suggests giving an exponen-85 tial evolution to the damage variable described thanks to displacement in the normal and shear directions [START_REF] Nguyen | A cohesive damageplasticity model for dem and its application for numerical investigation of soft rock fracture properties[END_REF]. A Mohr-Coulomb breakage criterion is also established.

The major drawbacks of this type of model are that they are only applied to normal and tangential forces respectively associated with the normal and 90 tangential spring. The damage does not influence the moments. Moreover, these models are only satisfied when one of the failure modes (I or II) is predominant.

When a mixed-mode is present, i. e. no mode of rupture is predominant, the energy required for fracture exceeds those related to pure modes that is physically inconsistent.

95

Nevertheless, the results obtained with the different cohesive zone models applied on spring links lead us to focus on this approach of fracture.

A cohesive zone model is used to represent phenomena happening in the process zone. This consists in the simulation of a fictitious crack along with a 100 cohesive interface for which a particular stress-displacement law is considered. This stress-displacement law consists generally of a first elastic regime up to the strength of the interface followed by a softening behavior describing the progressive damage of the stress and the stiffness leading to dissipated energy.

The complete separation of the interface is obtained when the cohesive energy 105 G f is dissipated (area under the stress-displacement law). The initial works on cohesive zone were defined for pure mode fracture with the fictitious crack concept [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]. Moreover, different damage evolution laws can be used depending on the modeled material [START_REF] Vandellos | Comparison between cohesive zone 930 models and a coupled criterion for prediction of edge debonding[END_REF]. In the literature, four main types of laws seem to emerge: a linear softening part [START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF], [START_REF] Snozzi | A meso-mechanical model for concrete under dynamic tensile and compressive loading[END_REF],a bi-linear softening slope [START_REF] Petersson | Crack growth and development of fracture zones in plain concrete and similar materials[END_REF], [START_REF] Dourado | Comparison of fracture properties of two wood species through cohesive crack simulations[END_REF], [START_REF] Phan | Influence of moisture content on mode i fracture process of pinus pinaster: evolution of micro-940 cracking and crack-bridging energies highlighted by bilinear softening in cohesive zone model[END_REF], 110 a tri-linear softening slope [START_REF] Park | Cohesive fracture model for functionally graded fiber reinforced concrete[END_REF], [START_REF] Evangelista | Three-dimensional cohesive zone model for fracture of cementitious materials based on the thermodynamics of irreversible processes[END_REF] and an exponential softening [START_REF] Xu | Void nucleation by inclusion debonding in a crystal matrix[END_REF], [START_REF] Zhu | Effect of characteristic 955 parameters of exponential cohesive zone model on mode i fracture of laminated composites[END_REF], [START_REF] Van Den Bosch | An improved description of the exponential xu and needleman cohesive zone law for mixed-mode decohesion[END_REF] (This last softening is illustrated in the figure 1). It can be noticed that only three parameters are required to describe an exponential softening, the initial stiffness K 0 , the maximum elastic displacement δ e , and the cohesive energy G f .

115

For mixed loading such as tensile and shear loading, the use of a cohesive zone model for one of the three pure modes can be suitable only in the case of a preferential direction of crack propagation. In other cases, the need for a coupling between pure mode models is essential, in line with the wide range of mixed-mode works carried out for example in [START_REF] Song | Cohesive zone simulation of mode i and 960 10[END_REF], [START_REF] Van Den Bosch | A new cohesive zone model for mixed-mode decohesion[END_REF] [START_REF] Venzal | Frictional cohesive zone model for quasi-brittle fracture: Mixed-mode and coupling between cohesive and frictional behaviors[END_REF]. Indeed, an uncoupled 120 model would lead to dissipated energy peak for mixed-mode rupture, i.e. more energy would be necessary for rupture happening than in pure modes, which is non-physical. In terms of energy, a coupled model ensures that energy dissipated for fracture in mixed mode is bounded by energies corresponding to pure modes and that energy tends to mode I energy (respectively mode II energy) when 125 mixed-mode fracture tends to pure mode I fracture (resp. mode II). Such a model is usually defined from pure modes. A parameter β monitoring mixedmode rate is expressed from displacements related to mode I and mode II such as β = δII δI . From this parameter and the ones provided by pure mode laws, the mixed-

J o u r n a l P r e -p r o o f

Journal Pre-proof

This cohesive zone model provides a useful approximation of fracture processes of materials and its installation at a mesoscopic scale allows reducing simulation time. Moreover, its ability to taking into account mixed mode at fracture appears as a powerful model applicable to the Discrete Element Method 135 (DEM). However, if the use of damaging link model through degraded spring is the subject of many works today, the solution consisting in using damage model in beam link is much less studied whereas DEM based on beam links is known to describe crack paths with better accuracy compared to those obtained from spring links [START_REF] Schlangen | Fracture simulations of concrete using lattice models: Computational aspects[END_REF].

140

The present work is inspired by the ones in [START_REF] Grassl | Meso-scale approach to modelling the fracture process zone of concrete subjected to uniaxial tension[END_REF][START_REF] Grassl | Meso-scale modelling of the size effect on the fracture process zone of concrete[END_REF] which propose damageable beam-lattice domain in 2D.

Here, an original damage model acting on tridimensional beam-lattice domain is proposed. Its capabilities to describe the general aspects of the quasi-145 brittle fracture and in particular the cracking paths are shown and discussed.

Nevertheless, in this first version of the model, the frictional aspects of the quasibrittle failure are not taken into account. Hence, the model will fail to describe the values of the residual displacements observed from an unloading following a damage in traction for example, [START_REF] Bažant | Analysis of work-of-fracture method for measuring fracture en-980 ergy of concrete[END_REF], [START_REF] Landis | Microstructure and fracture in three dimensions[END_REF]. Nevertheless, if the capacities of a 150 discrete element model to naturally describe the friction between particles are not exploited in this first version of the model(usually, a classical Hertz-Mindlin model is implemented for the broken particles that are in contact like in [START_REF] Hedjazi | Application of the discrete element method to crack propagation and crack branching in a vitreous dense biopolymer material[END_REF]), the contact between discrete elements is used in order to describe the unilateral effect of quasi-brittle materials and the capacities of the model to describe 155 quasi-brittle mixed mode I+II cracking is studied. In addition, the definition of a cohesive zone model at the virtual Voronoi interface shared by two neighboring discrete elements makes it possible to get rid of the problems of initializing the stiffness of the cohesive element observed with the finite element method. Indeed, the stiffness of the beams connecting the discrete elements is naturally 160 used for this purpose. Moreover, the introduction of this cohesive behavior, via a damage variable within the Euler-Bernoulli beams of the lattice, provides an [START_REF] Carpiuc-Prisacari | Experimental database of mixed-mode crack propagation tests performed on mortar specimens with a hexapod and full-field measurements. part ii: interactive loading[END_REF] are simulated. They consist of a mixed loading including tensile and shear forces as well as a bending moment applied to pre-cut specimens (single-and double-notched).

The Damaging Beam Model (DBM)

Beam lattice domain

Contrary to the finite element method, the continuum mechanics are not applicable for the discrete element method. Indeed, the principle of the discrete element method consists of discretizing a continuous media with particles. The discrete element method, created for granular materials, is based on the inter-190 actions between the particles. These interactions consist of forces and moments

J o u r n a l P r e -p r o o f

Journal Pre-proof applied on the discrete elements. In the case of cohesive media, the interactions between discrete elements are represented thanks to beam links. The forces and moments induced by the beam on the elements are defined from the theory of beams. Thus, for each beam, three forces, one normal and two tangential forces Newton law, accelerations are deduced. Thanks to an integration scheme, velocities and displacements are calculated from accelerations. In GranOO, an explicit scheme is used, the Verlet Velocity scheme [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a Journal Pre-proof silica glass with a blunt indenter[END_REF]. To ensure stability and convergence, a critical time step ∆t crit is determined and the calculation time step ∆t is taken with a ratio of ten percent i. e. ∆t = 0.1∆t crit .

225

The Damaging Beam Model (DBM) detailed thereafter consists in an improvement of this beam model inspired by cohesive zone models.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Mixed mode cohesive model

The damaging beam model consists of the definition of a damage variable D for each beam. The beam Young modulus is thus degraded inducing a decrease 230 in forces and torques :

E = (1 -D)E 0 , (1)
where E 0 and E correspond to the initial and current Young modulus. The evolution law of D is inspired from the mixed-mode cohesive zone model found in [START_REF] Venzal | Frictional cohesive zone model for quasi-brittle fracture: Mixed-mode and coupling between cohesive and frictional behaviors[END_REF], and in this work, all the three classical crack opening modes, I, II, and III are used to explicit the evolution of the damage.

235

The activation of the damage evolution is governed by a mixed initiation criterion expressed as :

k=I,II,III δ k δ e,k 2 = 1, (2)
where δ k is the current displacement for the mode k and δ e,k is the elastic limit also for each mode.

Assuming an exponential shape of the softening as shown in figure 1 for all 240 modes, the evolution of the damage variable can then be written as:

D = 1 - δ e δ exp      2 K 2 G f δ e -K δ e ( δ e -δ)      , (3)
where K, δ e and G f are the cohesive law parameters in the mixed space (i.e. a combination of cohesive parameters for each crack opening modes I, II, III. δ is the current displacement also projected in the mixed space. The projection in the mixed space will be described in the section 2.3.2.

245

Finally, a cut-off criterion is used to make the fracture of the beam effective when D tends to 1. It is assumed that total fracture, i.e. the beam is deleted, happens when the mixed stress in the beam reaches one percent of the mixed

J o u r n a l P r e -p r o o f

Journal Pre-proof peak stress after the damage activation. In the DBM, the criterion is translated into an ultimate displacement δ u :

250 δ u = δ e -      2 G f δ e -K δ e 2 K      ln 1 100 . (4)
2.3. Definition of the pseudo modes

Kinematic equivalences

The damage evolution law and the damage initiation criterion depend on four variables δ, δ e , K and G f that need to be evaluated at each iterations. These mixed law variables are defined from the variables of the three cohesive laws for 255 each pure mode. As these modes are connected to a specific opening kinematic that is not similar to a Bernoulli beam one, it is then required to propose an equivalence. The representation of the pure mode fractures is displayed in the table 1 as well as their assumed equivalences for the beam with respect to the surface A V oronoi (Figure 3). From these beam kinematics, a cohesive law is 260 defined for each of the pure modes. To determine a cohesive law, only three parameters are needed: a stiffness K, a maximum elastic displacement δ e and a cohesive energy G f . Note that cohesive laws relate stress to displacement, therefore the stiffness K must be in [N/m] and the displacement δ in [m].

Tension. First, the equivalence of mode I of failure appears obviously as a result 265 of the tensile effect. Thus, the displacement describing this mode of failure is the elongation of the beam. The stiffness defined for this mode corresponds to the tensile stiffness of the beam. Their respective expressions are given by:

         δ tension = L -L 0 K tension = E0 L0 (5)
J o u r n a l P r e -p r o o f

Journal Pre-proof

Shear and bending. The definition of mode II becomes less straightforward because of the two contributions, shear, and bending. Indeed, defining a displacement and a stiffness providing monitoring of shear contributions F y and F z as well as bending contributions T y and T z is more complicated. Nevertheless, it can be highlighted that both shear and bending loads induce an increase in beam length. Hence, the difference between the curvilinear length of the beam and the beam length obtained only by tensile loading could be a displacement monitoring shear and bending effects. The mode II stiffness is then determined by averaging these two stiffnesses. For Mode II stiffness, the shear and bending stiffnesses can be written as:

                       δ shear/bending = L curve -L K shear = E 0 I g L 2 0 R beam K bending = E 0 L 0 (6)
Torsion. For mode III failure, the parameters of the cohesive law are derived from torsional effects. Thus, the displacement is naturally expressed as that induced by the torsion angle. The mode III stiffness is also assigned the torsional stiffness. These displacement and stiffness are expressed as:

         δ torsion = R beam θ torsion K torsion = E 0 (1 + ν)L 0 (7)

Standard modes Bernoulli beam

Table 1: Conventional modes of fractures and the equivalent modes for a beam with respect to the surface A V oronoi

Pseudo modes definition

Three pseudo pure modes of fracture applied to the beam equivalent to the 270 conventional modes of fracture are defined. However, a mixed-mode combining three pure modes would be unnecessarily complicated. Indeed, as explained in the paper [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF], when the beam is deformed, the energy is principally composed of tension energy. In this way, tension effects can be considered as the major mode of fracture of the mixed-mode model and will be named as pseudo mode 275 I noted as mode I * . On the other hand, contributions from torsion loads are low and torsion and bending fractures appear on the same surface as shown on kinematic schemes in the table 1. Thus, a single pseudo mode II can be defined from both torsion and bending effects. More generally, all loading on the beam different from tension loading is considered belonging to pseudo mode II noted 280 as mode II * .

Moreover, the low contributions of bending and torsion loads ensure that the mode II * will not be predominant and thus, that the assumption of combining torsion and bending effects is suitable. By an energy analysis, the pseudo

                                   δ I * = L -L 0 K I * = E 0 L 0 δ e,I * = ε e,I L 0 G f,I * = A beam A voronoï G f,I (8)
Concerning the cohesive energy for mode I * , the choice is made to introduce directly the energy extracted from experimental tests in mode I failure.

Similarly, the maximum elastic displacement in mode I * is deduced from experimental maximum elastic strain noted ε e,I . However, the lattice character of the method implies a resetting step of the cohesive energy. In the aim of avoiding 295 another calibration step, the beam parameters are directly calculated from the material ones. However, it should be noted that energy conservation is not valid here. Indeed, the elastic calibration of the Young modulus induces a rigidity of the beam which is greater than the rigidity of the material. Thus, in the case of energy conservation, the beam fracture energy would be lower than the elastic 300 energy accumulated in the beam, which is not physically allowed. Hence, to estimate the fracture energy of the beam, it is assumed that the ratio of the fracture energy of the beam over the cross-section area of the beam is equal to the ratio of the fracture energy of the material over the surface area A V oronoi .

The area A V oronoi is the material surface modeled as defined in Section 2.1 305 and obtained from voronoï mesh. The energy release rate of each beam is set as shown in equation [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF].

Mode II * represents the second-order of elongation of the beam. The pa-

J o u r n a l P r e -p r o o f

Journal Pre-proof rameters of the pseudo mode II (mode II * ) are written as :

310                                δ II * = L curve -L + R beam θ torsion K II * = f shear K shear + f bending K bending + f torsion K torsion δ e,II * = αδ e,I * G f,II * = αG f,I * (9)
The mode II * displacement is expressed as the sum of the mode II and mode III displacements. To express a mode II * stiffness, a weighted average is obtained from the shear, bending, and torsion stiffnesses. The weights used to calculate the Mode II * stiffness are given by :

                                   f shear = F y + F z F y + F z + T y L + T z L δ shear/bending δ II * f bending = T y L + T z L F y + F z + T y L + T z L δ shear/bending δ II * f torsion = δ torsion δ II * ( 10
)
where f shear , f bending and f torsion correspond to the weight of the shear, 315 bending, and torsion stiffnesses in the estimation of the related to the pseudo mode II. Note that these weights are functions of the shear forces F y and F z and the bending torques T y and T z from the theory of beams.

Since the mode II * energy has no physical equivalent, the choice is made to 320 express this energy in proportion to the mode I * cohesive energy. For this pur-

J o u r n a l P r e -p r o o f

Journal Pre-proof pose, an α parameter is introduced and monitors the mixed-mode contribution and will need to be calibrated.

In summary, only one parameter needs a calibration step, the other parameters are set from experimental data.

325

Based on the pseudo pure modes previously defined, the parameters corresponding to a "pseudo" mixed-mode are expressed according to [START_REF] Camanho | Numerical simulation of mixed-mode progressive delamination in composite materials[END_REF] and [START_REF] Venzal | Frictional cohesive zone model for quasi-brittle fracture: Mixed-mode and coupling between cohesive and frictional behaviors[END_REF] such as:

                                                           δ = δ I * δ II * 1 + β 2 δ 2 II * + β 2 δ 2 I * K = K 2 I * + β 2 K 2 II * 1 + β 2
= G f,I * + G f,II * , (11)
where β corresponds to the pseudo mixed mode ratio with :

β = δ II * δ I * (12)
The mixed mode stiffness K, the maximum elastic displacement of the mixed 330 mode δ e and the mixed mode cohesive energy G f are functions of β. The last term G f can be written as a combination of two contributions of both pseudo mixed modes G f,I * and G f,II * . These decomposition will be used especially in section 2.6 for validation.

J o u r n a l P r e -p r o o f

Journal Pre-proof

The four terms in equation 11 are used to feed the initation and evolution 335 criterion defined previously in equations ( 2) and (3). At every step of the simulation, a different damage law is defined. Indeed, between two incremental times, displacements δ I * and δ II * evolve with the fluctuation of loads applied to the structure. An update of parameter β value is required, modifying all parameters from the damage evolution law. In this way, for a zero β value, the 340 pure mode I * law is restored while for an infinity β value, the mode II * law is found.

Closure model

A quasi-brittle material submitted to a tension loading exhibits a decrease of its stiffness in the tension direction due to the development of micro-cracks 345 oriented preferentially in the normal direction to one of the tension loadings.

Nevertheless, if this damaged material is then submitted to a compression loading applied in the direction of the preceding tension loading, the material exhibits a restoration of its initial stiffness due to the closure of micro-cracks. To reproduce this phenomenon called as the unilateral effect, closing forces and 350 torques are added. Thus, the sum of the forces and torques in the re-closed link is equal to the forces and torques of an undamaged beam loaded in compression,

i.e. such that :

   F T    =    F Damagebeam T Damagebeam    +    F closure T closure    (13)
In addition, micro-crack reclosure occurs within the entire representative elementary volume (REV) modeled by the reclosed beam. It is therefore necessary to the reclosure force F closure can be expressed as :

360 K closure = D E 0 A voronoï L 0 (14)
Note that the natural capacity of DEM to detect contact between discrete elements is used in order to activate closure forces (see equation 13). Nevertheless, if contact detection is activated, the frictional behavior between discrete elements is not considered in order to evaluate the capacities of the alone damage beam model to describe quasi-brittle failure. 365

Sequential Algorithm

The process carried out at each time step for each beam is detailed in figure 4 and includes the application of the damage model as well as the closure model. The analysis of the energy dissipated shows that the proposed damage model is consistent in terms of energy.

405

der to test the model in various fracture modes. Indeed, the main characteristics of the concrete failure under compressive and tensile loading are well-known and are here simulated to evaluate the relevance of the damage cohesive beam model.

The simulated responses in tension and compression will be simply qualitatively compared to the experimental ones while the ratio of compression and tensile 415 strengths obtained from simulation will provide a quantitative comparison with respect to the well-known ratio of 10 expected for quasi-brittle materials.

Prior to simulating tests, according to [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF], a calibration step is realized to determine the model parameters allowing to simulate concrete elastic behavior.

The elastic parameter values for macroscopic and model scales are displayed in 420 the table 3.

Macro scale

Model scale

E M = 39.4GP a E µ = 189.5GP a ν M = 0.2 ν µ = 0.3 r µ = 0.6
Table 3: Model and macroscopic elastic parameters for concrete, where rµ is the radius ratio i.

e. the ratio between each beam radius and the average radius of its linked discrete elements.

Tensile Test on a cylindrical specimen

A tensile test is performed on a cylindrical specimen. The dimensions of the specimen are 70 mm for the length and 35 mm for the radius as indicated in figure 6. in order to simulate a tensile strength of the specimen equal to f t = 4M P a.

430

Finally, the value for α is given arbitrarily at 8.

To remain in a quasi-static regime, the displacements imposed on the end cross-sections are such that a speed of 1.6 mm/s is ensured. The response stress-strain obtained from the simulation is plotted in figure 7. The trend displayed in the figure 7 corresponds to the expected quasi-brittle 435 behavior. Indeed, peak stress of 4.1 MPa is obtained which is consistent with the experimental values observed for this material [START_REF] Xenos | Calibration of nonlocal models for tensile fracture in quasi-brittle heterogeneous materials[END_REF], [START_REF] Ulfkjaer | Fracture energy of plain concrete beams at different rates of loading[END_REF], [START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF]. Furthermore, the sharp softening occurring after the peak stress is characteristic of the brittle character of fracture expected from a non-stabilized tensile loading. The localization of the damage occurs randomly at the height of the cylinder. Every beam of the network is allowed to damage, no artifices need to be added to obtain the localization.

The repeatability of this test on different initial beam networks is shown in figure 9 and allows to verify that the crack is randomly located in the cylinder but always around a horizontal plane. Finally, the influence of the number of discrete elements on the stress-strain response is studied, to quantify the number of elements ensuring a qualitative result. This convergence study is undertaken thereafter, from a sample consti-

465

tuted by a number of discrete elements ranged between 1 700 and 22 000 as shown in figure 10.

For any element numbers, the expected stress-strain shape is exhibited and the brittle tendency is recovered even for a coarser specimen. However, for the sake of accuracy, a minimum of 5700 elements is suitable to provide a deviation 470 under 10 % of finer mesh peak stress and energy. In other words, an elastically accurate specimen will display precisely the damaging behavior. On the left side, it can be noted that the fracture energy converges rapidly. Indeed, even for a coarse specimen, the average error on the fracture energy 480 obtained is close to zero. The case of a large error on the fracture energy is induced by the various initial beam networks available. For the finer specimens, this effect is noticeable. In addition, the impact of the variability of initial networks seems to become constant for specimens composed of more than 13000 discrete elements. This minimum influence of the variability can be quantified 485 as approximately 10%.

If the convergence of the energy is reached from 13 000 discrete elements, the convergence of the peak stress seems to require the use of 20 000 discrete

J o u r n a l P r e -p r o o f

Journal Pre-proof elements to be reached, as shown in figure 11b. Nevertheless, the shape of the stress-strain responses becomes acceptable from 5 700 discrete elements in so far 490 as all responses exhibit a post-peak consistent with the one expected in tension for a quasi-brittle material, i. e. marked from a first sudden decrease of the stress followed by a softening tail. As a consequence, for the peak stress as well as the fracture energy, the values obtained for the 5700 discrete elements remain fairly accurate with an error of 7.3% and 3% respectively. In such a case, the 495 tests on a 5700 sample are privileged.

Compression Tests

First, an interrupted tensile test followed by a compressing test is carried out on the cylindrical specimen. When no closure model is used for the simulation, the stiffness in the compressing reloading is equal to the stiffness observed during the unloading phase of the tensile test. Therefore, no unilateral effect 515 can be displayed without a closure model. The results for three different meshes using our closure model are illustrated in the figure 13. As for the convergence studied for the damaging behavior, a small number of discrete elements seems to ensure the accuracy of unilateral behavior. Indeed, for a sample composed

J o u r n a l P r e -p r o o f

Journal Pre-proof of 5700 elements, unilateral behavior is displayed. Similar results are obtained 520 for coarse and finer meshes as long as the number of elements in the sample is sufficient to model the elastic behavior. Our closure model is convenient and validated for unilateral behavior. However, due to the fact that frictional behavior is not activated in the model, this one fails to describe the value of the residual displacement at zero stress observed from experimental tests or simu-525 lated from other discrete element models considering frictional phenomena [START_REF] Sinaie | Application of the discrete element method for the simulation of size effects in concrete samples[END_REF], [START_REF] Nguyen | Discrete element modelling of fracture in quasi-brittle materials[END_REF]. Nevertheless, even if the friction phenomenon plays a major role for the quasi-brittle failure especially under compression loading and it is not taken into account in the present model, we choose to simulate a compression test in order to evaluate the capability of the alone damaging beam lattice model under this 530 kind of loading and because compression test involves shear failure required to calibrate the pseudo mode II parameters of the model. For such boundary conditions, a fracture of the plans parallel to the loading is expected. Figure 15 shows the evolution of the damage field with a top view of the cylinder. Numerous damaged beams are distinguishable through 545 the cylinder. Finally, multiple through-breaks of the cylinder appear in the cylinder, as expected. Nevertheless, damage does not localize in planes parallel to the compressive loading direction. Such a localization undoubtedly requires introducing some noise into the system as internal stresses which will be the object of a future study. Concerning the stress-strain curve presented in figure 16, the behavior of the specimen displays a smoother softening than in tension. In addition, the peak stress is 43 MPa. The compressing and tensile strength ratio is then approximately 10, as observed experimentally [START_REF] Snozzi | A meso-mechanical model for concrete under dynamic tensile and compressive loading[END_REF], [START_REF] Nouailletas | Ex-1030 perimental study of crack closure on heterogeneous quasi-brittle material[END_REF], [START_REF] Gangnant | Modeling of the quasibrittle fracture of concrete at meso-scale: Effect of classes of aggregates on global and local behavior[END_REF]. This first compressing testing is accurately simulated in terms of global mechanical behavior but leads to 555 diffuse damage compared to the experimental one observed for concrete. To test our model, particularly the crack path prediction, second compressing testing is carried out where a y-displacement is still applied but x and z displacements are blocked to simulate the frictional effect acting on the area of the boundary cross-sections frequently observed during experimental tests. From a mechanical behavior point of view, a similar stress-strain curve is obtained. A higher stress peak is displayed which is consistent with the larger surface of fracture happening in the specimen. In addition, the evolution of the damage field, illustrated in figure 17, reveals that damaged beams, represented in black in the figure, are more numerous than for the previous compressing test.

565

In order to observe the cracking path, the final damage field where the beams are fully damaged, i.e. for which D=1, is presented in Figure 18. For the sake of clarity, a slice of the specimen is made in order to analyse the cracking profile in one slice of the specimen. Damage is activated in many beams, including those located in cones, due to the embedded end conditions inducing high loading 570 in the specimen, especially in the material volume of the cones. However, a location of the heavily damaged beams as well as the cracking zone is noticeable along the conical surface, as experimentally. This type of break is provided by experimental testing due to the frictional effects between the sample and the plates of the device [START_REF] Simon | The effect of friction on the compressive strength of concrete specimens[END_REF], [START_REF] Caldarone | Importance of end surface preparation when 1070 testing high strength concrete cylinders[END_REF]. Consequently, our model is predictive for the 575 perfect compressive tests, with no frictional effect, and for endless frictional Let us remember that friction between discrete elements is not activated in the mode while frictional dissipation is well known to be a dominant character of the quasi-brittle failure [START_REF] Bažant | Analysis of work-of-fracture method for measuring fracture en-980 ergy of concrete[END_REF], [START_REF] Landis | Microstructure and fracture in three dimensions[END_REF] . This benchmark is an ideal candidate for testing our model due to its focus on the complex mixed-mode loading of mortar specimens. For these tests, a hexapod machine is used to apply complex loading. Indeed, the 6-axis testing machine, illustrated in figure 20 allows to apply up to 3 forces and 3 mo-610 ments simultaneously. The specimen is inserted between two plates subjected to these loads to better control the boundary conditions. For both tests simulated here, the applied loading can be split into three components: tension, shear, and bending, as shown in Figure 21.

J o u r n a l P r e -p r o o f

Journal Pre-proof Indeed, the positions of the randomly generated discrete elements do not coincide with those of the DIC measurements. Thus, an interpolation from experimental data is then performed. In a second stage, a calibration of the elastic and damage parameters is undertaken based on the experimental values of the mortar. The input parameters of the simulations are given in table 5. Indeed, for times greater than 0.1 s, a difference between the five numerical results and the experimental result is notable. The observation of the tensile and shear forces shows that the double notch test is faithfully reproduced, although the crack coalescence at the end of the simulation is dependent on the initial 700 beam network. However, the experimental result studied here corresponds to a specimen with a particular mesostructure. For a specimen with a different mesostructure, the evolution of the tensile and shear forces would not be identical to that presented here. Thus, due to these experimental and numerical variations of the tensile and shear forces, it can be assumed that the model As experimentally, a first crack is initiated at the right notch, observed in figure 27a. In a second stage, the specimen is unloaded to close the crack. A 710 new loading is then applied to the specimen to initiate a second crack without propagation of the first crack. The initiation of the second crack is displayed in figure 27b. In the final step, the loading is increased dramatically to reopen the two cracks and allow them to reorient and coalesce. This coalescence mechanism is highly unstable. The simulation results show a coalescence due to the 715 propagation of the second crack in contrast to the experimental results where the first crack propagates, illustrated in figure 27c. There are many reasons for these differences. First, the initiation part of the test has a great impact on the crack path. In addition, in this experimental work, a pre-cycle of elastic loading and unloading is applied to the specimen inducing a locally pre-stressed state, 720 different for each specimen tested. In addition, the mortar mesostructure influences the location of crack initiation. These pre-stress and mesostructure effects are not taken into account in the model and explain the differences observed during the unsteady propagation phase. stresses within the lattice. This will be tested using a tension-compression test [START_REF] Nouailletas | Ex-1030 perimental study of crack closure on heterogeneous quasi-brittle material[END_REF], which is more suitable for studying the dissipative mechanisms involved in crack reclosure.

Conclusion and Perspectives

765 Finally, the model was compared with the experimental results of the Carpiuc benchmark in which cracking is mostly induced in mixed mode I+II. Tests with complex mixed-mode loading on single and double notched specimens high-

J 1 δeFigure 1 :

 11 Figure 1: Stress-Displacement curve for a typical exponential damage law

J

  o u r n a l P r e -p r o o f Journal Pre-proof energetic control of the softening phase of the failure. The damage model is implemented within an Euler-Bernoulli beam-lattice rather than a Timoshenko one because if Timoshenko beams provide a more accurate description of the 165 contribution of shear, they would lead to new parameters to be calibrated. Their use would then require the development of a new calibration process. The reliability of the elastic calibration process for Euler-Bernoulli beams led us to choose these beams for our model. The model of the damaging beam and the closure management method are 170 detailed in the first part of this paper. Then, three main classical mechanical tests are carried out. As the last section of the paper is related to experimental comparisons, the authors analyze the results of these preliminary simulation tests to evaluate only the qualitative response to the model related to the calibration of the damage parameters. Typical failure mechanisms for quasi-brittle 175 materials are observed such as the unilateral effect or the ratio between compressing and tensile strengths. Tensile and compressive tests under monotonic and cycling loading are performed in this sense. Finally, to compare the model with experimental results, two tests from the recent Carpiuc benchmark proposed in the two previous ECCOMAS International Conferences on Computa-180 tional Modeling of Fracture and Failure of Materials and Structure (CFRAC)

195

  are calculated. Likewise, three moments, one torsion, and two bending moments are determined. These forces and moments and the induced kinematics for an Euler-Bernoulli beam can be found in the work of André et al. in [19]. These forces and torques are applied on each discrete element attached to 200 the beam. These elements can be of different shapes, such as spherical or polygonal shapes and may or may not be in contact with each other. Thereafter, two examples of discretization are shown in figure 2. The first one consists in spherical discrete elements which are not in contact while in the second one, polygonal discrete elements in contact are represented. In both cases, each 205 discrete element is connected to its neighbors thanks to links as beams.

  (a) Discretization with Spherical discrete elements linked by beams (b) Discretization with Voronoï discrete elements linked by beams

Figure 2 :

 2 Figure 2: Examples of discretization

Figure 3 :

 3 Figure 3: Scheme of material volumes associated to spherical discrete elements

  -of-plane shear Torsion

  285 mode II (mode II * ) can be considered as an energy E II * , the sum of energies of mode II E II and mode III E III of fracture which will be dissipated when the surface shared by Modes II and III is fractured. This new energy can be as well expressed as a function of a stiffness K II * and a displacement δ II * . 290 Consequently, two pseudo pure modes, mode I * and mode II * are defined from the equivalences of the 3 conventional pure modes detailed in the previous paragraph. Mode I * is similar to mode I, and corresponds to first order elongation of the beam. Thus, fracture parameters of the pseudo mode I (mode I * ) J o u r n a l P r e -p r o o f Journal Pre-proof can be expressed as :

δ 2 e* + β 2 δ 2 e,I * G f,I * + β 2 δ 2 e,I * δ 2 e,II * + β 2 δ 2 e

 22222 e = δ e,I * δ e,II * 1 + β 2 δ 2 e,II * + β 2 δ ,I * G f,II *

355

  to take into account the part of the REV represented by void. Thus, during the reclosure of two discrete elements, the forces F closure are transmitted along a surface that does not correspond to the surface of the beam but to the Voronoi surface common to both elements, illustrated in figure3. The stiffness related J o u r n a l P r e -p r o o f Journal Pre-proof

Figure 4 :

 4 Figure 4: Sequences for the application of the damage model and closure for each beams.

2. 6 . 4 α 5 Table 2 : 2 ]

 64522 Mixed mode Test on a single beamPhysically, the mixed-mode fracture energy of the beam has to be bounded 370 by the pure mode fracture energies. For various materials, the most energetic mode (i.e. the mode exhibiting the larger cohesive energy) corresponds to pseudo mode II while the least energetic corresponds to pseudo mode I. Thus, mixed-mode should require cohesive energy ranged between those of pure modes, and the evolution of this one with respect to the mixed-mode indicator β should 375 be monotonic. As our model is based on energy driving, it appears essential to verify this energy trend.The non-common mode definition used in the model leads us to take precautions and to achieve a deeper checking phase. For this purpose, a test is carried out on a single beam loaded with a mixed tensile-torsion boundary condition to 380 simulate various mixed-mode cases.Thus, a displacement in the longitudinal direction is imposed at a beam end. The other end is subjected to a torsion rotation around the beam axis.With such a loading, pseudo mode I displacement is provided by longitudinal displacement while pseudo mode II displacement is induced by torsion angle.385Hence, the different pseudo mixed mode ratio (β) values are generated in the beam with the variations of the displacement and rotation imposed at the beam ends. When β is equal to 0, a pure tension loading is applied, then is increasing with the rotation value on the beam.The input parameters for the damage law are given in the table 2. 390 G f,I 60 J.m -2 ε e,I 10 -Damage law parameters for the validation tests on a single beam. The evolution of the mixed-mode fracture energy as a function of the mixedmode indicator β is plotted in the figure 5. The mixed-mode fracture energy ( G f ) of the beam corresponds to the sum of the components of mode I * ( G f,I * ) J o u r n a l P r e -p r o o f Journal Pre-proof and of mode II * ( G f,II * ) of the energy and here can be either tensile energy or torsion energy. Thus, both mode I * and mode II * energy components are drawn. 395 It can be noted that each component (pseudo mode I energy or pseudo mode II energy) is lower than the corresponding pure mode energy value. Then, the mixed-mode fracture energy is bounded by pure mode I * fracture energy and mode II * fracture energy. Moreover, it can be observed that when β tends towards 0, the mixed-mode fracture energy tends to the one expected from 400 pure mode I * as well as its tensile component. Similarly, for a high value of β parameter, the mixed-mode fracture energy reaches the one associated with mode II * . Mixed-mode Fracture Energy G f Mode I * Fracture Energy G f,I * Mode II * Fracture Energy G f,II * G f,I * G f,II *

Figure 5 :

 5 Figure 5: Fracture Energies as function of the pseudo mixed mode ratio β on a single beam.

  425

Figure 6 :

 6 Figure 6: Voronoi mesh of the cylindrical specimen (values in mm)

Figure 7 :

 7 Figure 7: Stress-Strain response obtained from the DBM for a tensile test on a cylindrical specimen

Figure 8 :

 8 Figure 8: Evolution of damage fields during a tensile test

Figure 9 :

 9 Figure 9: Stress-strain response for five samples of different initial distribution under a tension loading.

Figure 10 :

 10 Figure 10: Study of the convergence on a cylinder under tension for an increasing number of elements.

Figure 11 :

 11 Figure 11: Study of the convergence on a cylinder under tension

Figure 12 :

 12 Figure 12: Stress-stain curve for tensile test with three different α values

  510

Figure 13 :

 13 Figure 13: Stress-Strain curve for tensile-compressing tests

Figure 14 :

 14 Figure 14: Boundary conditions applied during the friction-less compression test (on the left) and the clamped compression test (on the right)

Figure 15 :

 15 Figure 15: Evolution of damage fields during a compressing test with free extremities

Figure 16 :

 16 Figure 16: Stress-strain response for compression tests with no friction and infinite friction at the plates.

J

  o u r n a l P r e -p r o o f Journal Pre-proof effects.

Figure 17 :Figure 18 :

 1718 Figure 17: Evolution of damage fields during a compressing test with clamped extremities

Figure 19 :

 19 Figure 19: Influence of α value on frictionless compressing stress-strain curve

  . As a consequence, taking friction into account in a future version of the model should lead to consider a decrease of the value of the pseudo mode II fracture energy G f II * and/or of the parameter α because a 595 large amount of energy should be dissipated through friction under compression loading. Thus, even if the fracture parameters of pseudo mode II are certainly overestimated due to the absence of friction, the three parameters of our model are determined from compression and tensile tests and the model can be used for simulation of experimental tests in which most of the cracking is induced in 600 mixed-mode I+II. 4. Confrontation to experimental mixed-mode cracking : Benchmark Carpiuc In this section, the model is confronted with experimental work. For this purpose, particular attention is directed to two trials carried out in the frame-605 work of the Carpiuc benchmark [71] where data are available in the Zenodo website [72]

Figure 20 :

 20 Figure 20: Experimental set-up developped for the Carpiuc benchmark [71]

Figure 21 :

 21 Figure 21: The three elementary loading applied during the tests: tension, shear and torsion [71]A digital image correlation (DIC) is performed to measure the displacements

  (a) simple notch specimen (b) Double notch specimen

Figure 22 :

 22 Figure 22: 2D Voronoi Specimens used for the simulation of the benchmark tests

8 Table 5 :

 85 G f,I 114.6 J.m -2 ε e,I 2.2 ×10 -4 α Elastic and Damaging law parameters for the benchmark mortar To avoid specimen damage near the lower and upper zones of loading application, a correction factor is applied to the damage variables. In fact, the corrective factor F D follows a top hat law. As shown in figure 23 for a single notched specimen, the corrective factor varies between 0 at the edges of the specimen and 1 in the core of the specimen and thus allows the damage on the 635 edges of the specimen to be canceled out without modifying the damage values on the central part of the specimen.

Figure 23 :Figure 24 :Figure 25 :Figure 26 :

 23242526 Figure 23: Evolution of the corrective factor F D value within a single notched specimen

725A

  new damage beam-lattice model inspired by cohesive zone models has been presented. If this model is dedicated to discrete element modeling, it is implemented in this first version without resorting to the friction modeling in order to study the alone capabilities of the damage beam-lattice to describe quasi-brittle cracking mostly induced in mixed mode I+II. Only the contact de-730 tection between neighboring discrete elements is activated in order to describe unilateral behavior. The main advantage of this model is the reduced number of parameters which can be easily determined from two experimental tests, one tensile and one compression test. Five qualitative tests have been performed on concrete-like specimens. Among these tests, only two, monotonic tensile and 735 compressive tests are required to calibrate the damaging beam model. A second compressing test has been used to verify that the set of calibrated parameters is suitable and in accordance with the classical shape of mechanical response of quasi-brittle materials. The peak stresses and cracking paths expected are displayed. In particular, the compressing and tensile peak stress ratio is es-740 timated at 10 which is the typical value observed for concrete. Nevertheless, comparisons with data from experimental tests need to be performed to finally validated either the calibration protocol and the prediction capability of the proposed approach. In addition, it has been shown that a microscopic damaging behavior induces a macroscopic damaging behavior. The methodology 745 of energy calibration which has been presented allows displaying macroscopically energy dissipated equally to the cohesive energy input. Furthermore, it has demonstrated that convergence appears for specimens with few elements. The elasticity convergence, reached for approximately 5 000 discrete elements, ensures a good accuracy of the damaging behavior. Thus, the damaging beam 750 model could be a good candidate for structure simulation. The mechanism of crack reclosure has also been studied. The interrupted tensile test with a release followed by a compressing test has shown that the contact model developed allows taking into account the unilateral behavior of the material. However, as expected due to the fact that frictional behavior between 755 discrete elements is deactivated, the model fails to describe the value of residual displacement at zero stress observed from experiments or simulated by other discrete models in the litterature that already take into account the coulomb law for friction. Special interest should be focused on these frictional aspects during further work. Indeed, better modeling of the frictional contribution could result 760 in a larger residual displacement and hysteresis effects. The improvement of the present model is being developed to introduce friction effects as well as internal

Table 4 :

 4 Macroscopic Energy dissipated by the different samples

Preprint submitted to Journal of L A T E X TemplatesDecember 10, 2021

mode damage evolution law is updated at every step of the simulation.

Journal Pre-proof 3. Calibration and qualitative response of the model

Material

The model developed here, is built to be suitable for any quasi-brittle material simulation. Among quasi-brittle materials, concrete is subject of numerous studies available in literature [START_REF] Grégoire | Size and boundary effects during failure in quasi-brittle materials: Experimental and numerical investigations[END_REF], [START_REF] Vu | Revisiting statistical size effects on compressive failure of heterogeneous materials, with a special focus on concrete[END_REF], [START_REF] Fernández-Canteli | Determining fracture energy parameters of concrete from the modified compact tension test[END_REF], [START_REF] Lee | An experimental study on fracture energy of plain concrete[END_REF], [START_REF] Nouailletas | Ex-1030 perimental study of crack closure on heterogeneous quasi-brittle material[END_REF], [START_REF] Nguyen | Initiation and propagation of complex 3d networks of cracks in heterogeneous quasi-brittle materials: 1035 Direct comparison between in situ testing-microct experiments and phase field simulations[END_REF] and is chosen here in or-

J o u r n a l P r e -p r o o f

Journal Pre-proof . Firstly, during the simulation, the forces are measured at the extremities of the digital specimen corresponding to the area where the image correlation is performed. However, the experimental force values are measured at the hexapod level. Thus, due to these differences in the force measurement areas, and the use of image correlation data, noise in the numerical results are noticeable. At 655 the beginning of the test, for a time of less than 0.1 s, the experimentally and numerically measured forces are similar. The only difference lies in the peak force in tension. Indeed, due to the use of specimens composed of 8,000 discrete elements, it has been shown previously that convergence at the peak stress level was not achieved. As a result, an underestimation of the peak force at crack 660 initiation is observed numerically. Concerning the shear force, the experimental and numerical results are confounded for times lower than 0.1 s, i.e. during the crack initiation phase. During the crack bifurcation and branching phases, i. e.

for times greater than 0.1 s, differences are observed in the experimentally and numerically measured shear forces. These differences can be explained, on the

Appendix : curvilinear lenght computation

Curve length determination is the result of mathematical work. Indeed, the curve length is calculated from the space parametrization of the beam.

1080

Following the figure 28 presenting the framework for the cylindrical beam, the y and z expressions as functions of x are related to curve length with the equation: Beam displacements can be defined by a 3 rd polynomial function as following:

To determine the constants of these polynomials, it is necessary to write 4 equations for each position variable y and z. However, the only information known is on discrete element displacements and rotations. In other terms, the y and z positions of extremities x = 0 and x = L are given and provide two equations. The two last equations missing are obtained thanks to rota-1090 tions which are derivative values in x = 0 and x = L. Taking into account these boundary conditions allows system solving. From the resulting y and z expressions and as well their derivatives, curve length is obtained through numerical integration.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Declaratio if ioterettt

☒ The authors declare that they have no known competng fnancial interests or personal relatonships that could have appeared to infuence the work reported in this paper.

☐The authors declare the following fnancial interestsppersonal relatonships which may be considered as potental competng interests: