Damien André
email: damien.andre@unilim.fr

Jérémie Girardot

Cédric Hubert

A novel DEM approach for modeling brittle elastic media based on distinct lattice spring model

Keywords: Discrete Element Method, Voronoi, Elastic, Brittle fracture, Calibration, Continuum

The Discrete Element Method (DEM), also known as Distinct Element Method (DEM), is extensively used to study divided media such as granular materials. When brittle failure occurs in continuum such as concrete or ceramics, the considered media can be viewed as divided. In such cases, DEM offers an interesting way to study and simulate complex fracture phenomena such as crack branching, crack extension, crack deviation under coupled mode or crack lip closure with friction.

The fundamental difficulty with DEM is the inability of the method to deal directly with the constitutive equations of continuum mechanics. DEM uses forcedisplacement interaction laws between particles instead of stress-strain relationships. Generally, this difficulty is bypassed by using inverse methods, also known as calibration processes, able to translate macroscopic stress-strain relationships into local force-displacement interaction laws compatible within DEM frameworks. However, this calibration process may be fastidious and really hard to manage.

The presented work proposes to improve the Distinct Lattice Spring Model in order to deal with non-regular domains, by using Voronoi cells, which allows to completely fill the volume space of discrete domains. With this approach, the rotational effects must be included in the contact formulation, which enables the management of large rigid body rotations. This work also introduces a simple

Introduction, motivation of study

Since the earlier work of Cundall and Strack (1), the Discrete Element Method (DEM) has been widely used for modeling dynamic behaviour of granular media. The aim of this method, also called Distinct Element Method (DEM), is for solving a class of problems that can not be treated by continuous or analytical approaches such as packing problems or hourglass flows of granular media.

Based on the pioneering works of Hans J. Herrmann and al. concerning fracture simulation [START_REF] Pompe | Statistical models for the fracture of disordered media[END_REF], the wide range of methods grouped under the "DEM" appellation was enriched for investigating crack propagation problems in brittle media. Several methods were developed such as the Lattice Element Methods (LEM) [START_REF] Schlangen | New method for simulating fracture using an elastically uniform random geometry lattice[END_REF][START_REF] Mora | A lattice solid model for the nonlinear dynamics of earthquakes[END_REF] and the Bonded Particle Element Methods (BPEM) [START_REF] Potyondy | A bonded-particle model for rock[END_REF]. Both LEM and BPEM implement interaction laws between elements (or particle) responsible for representing the apparent mechanical behaviour. Failure criteria are introduced in order to break these interactions and simulate crack propagation. Mainly, LEM implements material points as elements while BPEM implements material bodies as elements. Indeed, this distinction is not so clear, so both models are considered here as 'DEM'.

Today, the most popular model appears to be the BPEM implemented in the PFC software (Itasca Inc.) that uses the parallel bond model [START_REF] Potyondy | A bonded-particle model for rock[END_REF]. Currently, more than 1,800 citations of this article are numbered. More recently, the flat joint model [START_REF] Potyondy | Parallel-bond refinements to match macroproperties of hard rock[END_REF] was developed in order to overcome limitations of parallel bonds concerning the compressive-to-tensile strength ratio of rocks. Indeed, many variations of BPEM and LEM exist. These variations often concern the interaction range, the geometrical arrangement of packing, the rheological model of the interaction and the failure criterion. Most of the interaction models are based on a pair of tangential and linear dash-pot springs. In such models, many efforts were done in the past years for enriching the failure criterion [START_REF] Oñate | A local constitutive model for the discrete element method. application to geomaterials and concrete[END_REF][START_REF] Scholtès | Modelling progressive failure in fractured rock masses using a 3d discrete element method[END_REF][START_REF] Potapov | Mixed dem/fem modeling of advanced damage in reinforced concrete structures[END_REF]. Other variations exist such as the cohesive beam model [START_REF] Lilliu | 3d lattice type fracture model for concrete[END_REF][START_REF] Carmona | Fragmentation processes in impact of spheres[END_REF][START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF][START_REF] Halász | Effect of disorder on shrinkage-induced fragmentation of a thin brittle layer[END_REF] where discrete elements are bonded by Euler-Bernouilli beams able to work in tensile, compressive, twisted and bending modes. In such cases, fracture criteria are mainly based on maximal bond stresses [START_REF] Carmona | Fragmentation processes in impact of spheres[END_REF] or strains [START_REF] Tavarez | Discrete element method for modelling solid and particulate materials[END_REF][START_REF] Shiu | Discrete element modelling of missile impacts on a reinforced concrete target[END_REF]. Again, many variations exist such as fracture criterion based on virial stress computation [START_REF] André | Using the discrete element method to simulate brittle fracture in the indentation of a silica glass with a blunt indenter[END_REF][START_REF] André | A discrete element thermo-mechanical modelling of diffuse damage induced by thermal expan-sion mismatch of two-phase materials[END_REF]. A review of the main used models can be found in [START_REF] Lisjak | A review of discrete modeling techniques for fracturing processes in discontinuous rock masses[END_REF].

Among these models, challenging problems often concern the ability of DEM to reproduce complex fracture phenomena such as coupled mode fracture, internal friction, compressive hardening or fracture softening. Indeed, such DEM models are not able to model easily simple elastic behaviour [START_REF] Zhao | On the linear elastic responses of the 2d bonded discrete element model[END_REF][START_REF] Celigueta | Accurate modelling of the elastic behavior of a continuum with the discrete element method[END_REF] defined by only two physical (apparent) parameters : a Young's modulus and a Poisson's ratio. The number of local parameters to manage is high regarding the few number of apparent parameters. In addition, influences of local parameters on apparent behaviours are difficult to predict because of possible coupled effects. For example, with the parallel bond model, the apparent Poisson's ratio can be overcome by adjusting the normal to shear stiffness ratio. However, apparent Poisson's ratio value depends also on the local coordination number. So, apparent properties depends on both mechanical and structural properties. These coupled effects between structural and mechanical properties is complex to address especially for random networks [START_REF] Ostoja-Starzewski | Lattice models in micromechanics[END_REF]. Calibration processes are currently used in many simulation works to overcome this difficulty, for example in [START_REF] Leclerc | On the suitability of a discrete element method to simulate cracks initiation and propagation in heterogeneous media[END_REF][START_REF] Delaplace | Discrete modeling of cracking of brittle materials in large relative motion and localization problem[END_REF][START_REF] Yang | Modeling progressive delamination of laminated composites by discrete element method[END_REF].

In such cases, DEM is difficult to manage even for simple behaviours, such as the purely elastic materials, in comparison with the Finite Element Method (FEM). Because of coupled effects with structural parameters, DEM often involves calibration processes in order to fit local parameters. These calibrations are non-normalized processes, time consuming and really difficult to manage. In addition, many couple of solutions may exist for addressing one set of apparent parameters. In such case, how to choose the best one ? Are these set of local parameters really equivalent ? In the opinion of the authors, these difficulties and questions limit the usage of DEM for solving fracture mechanic problems. Mainly, the mechanical engineering community is most familiar with classical continuum mechanics that requires no calibration processes.

A key point for DEM improvement regarding continuous mechanics is the interaction models between discrete elements. In the DEM framework, these interactions are widely expressed in term of force-displacement and not directly in term of stress-strain laws. The presented work addresses this difficulty by combining and improving several techniques for modeling brittle elastic materials without calibration process. This work was motivated by the recent study of Celigueta and al. in (20) where authors propose a modified contact law that "takes into account the surroundings of the contact, not just the contact itself. This brings in the missing terms that provide an accurate approximation to an elastic continuum, and avoids calibration of the DEM parameters for the purely linear elastic range". Here, the cited authors highlight the main problem of discrete to continuum transition. In most of DEM models, computing contact reaction forces involves only information from the two discrete elements in contact such as their relative positions. However, in the continuum elastic framework, these informations are not sufficient for computing reaction forces considering Poisson's and shear effects in the orthogonal directions of the contact. The solution proposed by Celigueta and al. consists in enriching the contact information thanks to virial stress tensors [START_REF] Zhou | A new look at the atomic level virial stress: on continuummolecular system equivalence[END_REF]. However, the authors could not run stable computations with the proposed solution with explicit time scheme.

A different approach was developed in [START_REF] Mariotti | Lamb's problem with the lattice model mka3d[END_REF][START_REF] Monasse | An energy-preserving discrete element method for elastodynamics[END_REF][START_REF] Christian | A least-squares coupling method between a finite element code and a discrete element code[END_REF] by proposing an energetic homogeneization of a discrete domain. This is using a Cosserat equivalent media giving cohesive forces but also torques thanks to the quaternions. This approach appears to be the single one reaching the bonding between the element scale and the continuum one. Nevertheless, the use of quaternions and the non-description of the full strain tensor involves limitation concerning complex fracture criteria [START_REF] Mariotti | Modeling of the fragmentation by discrete element method[END_REF].

Another approach, named Distinct Lattice Spring Model (DLSM), was proposed by Zhao and al. in (30). In this work, Cauchy's strain tensors are deduced from the displacement of surrounding elements thanks to least square regressions. Contact reactions forces are deduced from classical continuum mechanics using this strain tensors. In this study, authors recommend the usage of regular arrangement. Indeed, regular arrangement conduct to non-desired apparent orthotropy. Regular arrangement promotes non realistic crack paths as it was highlighted by Schlangen and al. in (3). In addition, DLSM does not take into account the local rotation of elements and some questions remain on the ability of the method to deal with large transformations and, more specifically, rigid body rotations.

The model proposed in the current study tackles these limitations. The DSLM is reformulated here to deal accurately with non-regular domains. While the original DLSM model involves only material points (with a radius), no rotational effects are taken into account. The proposed model, inspired from the work of Mariotti, involves 3D polyhedral elements given by Voronoi tessellation that completely fills a volume space [START_REF] Halász | Effect of disorder on shrinkage-induced fragmentation of a thin brittle layer[END_REF][START_REF] Mariotti | Modeling of the fragmentation by discrete element method[END_REF]. With such a description, rotational effects and moments of polyhedral element must be taken into account. Based on the geometry of elements, a natural method is proposed here for computing moments and rotation of elements. A simple benchmark, a deformable sphere that bounces on an inclined plane, is carried out in order to highlight the ability of the method to deal with large transformations.

The last part of this study is dedicated to the fracture criterion for purely elastic brittle materials. A simple criterion is introduced able to reproduce qualitatively and quantitatively the failure in complex tests such as the indirect tensile test, named Brazilian test, very popular in rock mechanics and civil engineering [START_REF] Andreev | A review of the brazilian test for rock tensile strength determination. part i: calculation formula[END_REF].

Displacement field computation

In this section, the formulation proposed by [START_REF] Zhao | A 3d distinct lattice spring model for elasticity and dynamic failure[END_REF] is reformulated for isotropic elastic media. Strain tensors and displacement fields are computed through local frames compatible with large transformations. This formulation is thus able to deal with large transformations as illustrated in section 8.1.

Problem definition

Supposing a polyhedral discrete element i surrounded by n discrete elements j 1 , j 2 , ..., j n as shown in 2D on Figure 1. The contact surface between i and j is denoted S ij . The points O i and O j are related to the geometric centers of discrete elements i and j. Because of the Voronoi tessellation of the domain (see section 7), the (O i , O j ) line is orthogonal to its related contact surface S ij . The infinite plane P ij that contains S ij cut the (O i , O j ) line at its midway. In addition, three different frames are managed in this study:

1. the global Galilean frame F g . This frame is unique and fixed. 2. local frame F i attached to discrete element i. The local frames have no relative movement relatively to their attached discrete elements. The origin of F i is located at the center of the related discrete elements i. At the initial state t 0 , the axes of F i are aligned with the axes of the global frame F g . 3. local contact frame F ij attached to a contact between two discrete elements i and j. The origin of F ij is located at the intersection between the contact plane P ij and the (O i , O j ) line. Its first axis x ij is normal to the contact surface and oriented from i to j.

Displacement field computation with least square method

The main objective of this theoretical part is to compute a Cauchy's strain tensor associated to each discrete element i. The Voigt's hypothesis is applied. It supposes a linear displacement field. With this assumption, the local displacement field d i|F i related to a discrete element i and expressed in its local frame F i is :

d i|F i (x, y, z) =   u(x, y, z) v(x, y, z) w(x, y, z)   =   a 1 x + b 1 y + c 1 z a 2 x + b 2 y + c 2 z a 3 x + b 3 y + c 3 z   F i = a • x + b • y + c • z (1) O i O j y ij x ij O ij xg yg Og x i y i local frame F i Galilean frame Fg contact frame F ij discrete element (i) discrete element (j) (O i , O j ) contact surface S ij contact plane P ij geometric center
a =   a 1 a 2 a 3   b =   b 1 b 2 b 3   c =   c 1 c 2 c 3   (2)
In equation 1, the special case where x=y=z=0 corresponds to the origin of F i , e.g, the O i point. At these coordinates, the related displacement d i|F i is null :

d i|F i (0, 0, 0) = d i|F i (O i ) = 0 (3)
This property is expected because the origin of the frame F i is attached to O i . In addition, the displacement field d i|F i shall respect the relative displacements of the discrete elements j in contact with the regarded discrete element i as :

d i|F i (O j ) = d j|F i (4)
where d j|F i is the displacement of the discrete element j computed at its center O j and expressed in the local frame F i . Finally, the displacement d j|F i of j can be deduced from its position as follow :

d j|F i = p j|F i -p 0j|F i (5)
where :

• p j|F i is the current position vector expressed in F i of the discrete element j and

• p 0j|F i is the initial position vector expressed in F i of the discrete element j at the initial iteration.

Thus, local displacement field must satisfy the following equations for each discrete element i in contact with j:

d i|F i (x j , y j , z j ) = p j|F i -p 0j|F i (6)
where (x j , y j , z j ) are the initial coordinates of O j expressed in F i :

x j = p 0j|F i • x i y j = p 0j|F i • y i z j = p 0j|F i • z i
So, it gives for the n discrete elements j in contact with the regarded discrete i the n following equations that form a system :

           a • x 1 + b • y 1 + c • z 1 = d 1|F i a • x 2 + b • y 2 + c • z 2 = d 2|F i . . . a • x n + b • y n + c • z n = d n|F i (7)
where :

• a, b and c are the three vector unknowns to be determined.

• (x j , y j , z j ) are the initial coordinates of O j expressed in the local frame F i . These coordinate values can be deduced from the initial positions of the discrete elements i and j.

• d j|F i are the displacements of a discrete element j in contact with i. These displacement values can be deduced from the current positions and orientations of the discrete elements i and j.

To determine a, b and c, three surrounding discrete elements j in contact with the regarded discrete element i are required.

In most cases, more than three discrete elements are in contact with a given discrete element. For discrete domains built thanks to Voronoi tessellation (see section 7), the average coordination number is about 11. So, the system of equation is overdetermined. The solution proposed here is to estimate a, b and c thanks to least square method (the least square method is a really common method and it will not be detailed here). Indeed, three separated cases must be treated :

1. less than three surrounding contacts. In this case, the computation of the displacement field is considered as null. 2. three surrounding contacts. In this case, the linear system 7 can be exactly solved. 3. more than three surrounding contacts. In this case, the linear system 7 is overdetermined. The solution is approximated thanks to least square method.

The proposed method allows us to determine the three vector unknownq a, b and c thanks to a least square method. So, at this stage, the local displacement field d i|F i (x, y, z) of a discrete element i expressed in F i can be computed with the proposed method. The next section will treat about the computation of strains.

Cauchy's strain tensor computation for a discrete element i

The infinitesimal strain theory postulates that the Cauchy's strain tensor ε is the gradient of the continuous displacement field d as :

ε = ∇ d (8)
Considering :

d(x, y, z) =   u(x, y, z) v(x, y, z) w(x, y, z)   (9)
and :

ε =   ε xx ε xy ε xz ε yx ε yy ε yz ε zx ε zy ε zz   (10)
It gives :

ε xx = ∂u ∂x and ε xy = 1 2 ∂v ∂x + ∂u ∂y (11)
ε yy = ∂v ∂y and ε xz = 1 2 ∂w ∂x + ∂u ∂z (12)
ε zz = ∂w ∂z and ε yz = 1 2 ∂w ∂y + ∂v ∂z (13)
Now, considering the local displacement field d i|F i computed with the least square method which is detailed in the previous section, the local strain tensor ε i|F i for the discrete element i expressed in the local frame F i can be deduced as :

ε i|F i = ∇ d i|F i (14)
Considering equation 1, 11, 12 and 13, the local strain tensor ε i|F i can be expressed as follows :

ε i|F i = 1 2   2a 1 a 2 + b 1 a 3 + c 1 a 2 + b 1 2b 2 b 3 + c 2 a 3 + c 1 b 3 + c 2 2c 3   F i (15)
In conclusion, the local strain tensor ε i|F i associated to discrete element i can be determined from the a(a 1 , a 2 , a 3 ), b(b 1 , b 2 , b 3 ) and c(c 1 , c 2 , c 3 ) coefficients previously determined (see section 2.2). Here, it is important to note that the Voigt's hypothesis induces a uniform and constant strain inside a discrete element i. However, in order to compute the contact reaction forces, strains should be determined at the contact interface between a discrete element i and j. The next section will be dedicated to this computation.

Cauchy's strain tensor at contact interface between a discrete 1 and 2

The previous section has shown how to compute the strain tensor related to each discrete element i by taking into account its neighbors j. Indeed, for computing contact reaction forces, these strain tensors must be interpolated at their related contact interfaces. Considering two discrete elements (noted 1 and 2) in contact as illustrated on Figure 2.

The Cauchy's strain tensors ε 1|F 1 and ε 2|F 2 are perfectly known thanks to the previous section. The strain tensor related to the contact at the point interface O 12 is simply considered as the average of these two strain tensors. However, these two tensors are expressed in two different frames. So, a preliminary calculation is needed to express both tensors in the local contact frame F 12 (see section 2.1) to give ε 1|F 12 and ε 2|F 12 . So, the strain tensor ε 12|F 12 at the contact interface can be computed as :

O 1 O 2 y 12 x 12 O 12 contact frame F 12 discrete element 1 discrete element 2 (O 1 , O 2 ) contact surface S 12 contact plane P 12
ε 12|F 12 = 1 2 ε 1|F 12 + ε 2|F 12 (16)
This simple interpolation method comes from the geometrical configuration of contacts where the contact surface is located at midway of O 1 and O 2 . This property is ensured by Voronoi diagrams used for building the discrete domains in this study (see section 7). Indeed, the obtained strain tensor ε 12|F 12 should be considered as non-local regarding the contact itself because it combines information from all the discrete elements in contact with 1 and 2. In order to force information from the local contact itself, the normal component of the normal strain is extracted as highlighted here : [START_REF] André | A discrete element thermo-mechanical modelling of diffuse damage induced by thermal expan-sion mismatch of two-phase materials[END_REF] and this value is forced to be equal to the engineering strain :

ε 12|F 12 =   ε xx ε xy ε xz ε yx ε yy ε yz ε zx ε zy ε zz   F 12
ε xx = L -L 0 L 0 = ∆L L 0 ( 18
)
where L is the current interaction length and L 0 is the initial length. The lengths L and L 0 are simply defined as the norm of the ---→ O 1 O 2 vector as :

L = ---→ O 1 O 2 (19)
To summarize, the computation of the contact strain tensor involves the three following steps.

1. Expressing the local strain tensors ε 1|F 1 and ε 2|F 2 computed with the least square method (see section 2.2) into the contact local frame F 12 . This step involves only changing frame operation with rotation matrix that can be deduced from quaternion :

ε 1|F 1 → ε 1|F 12 ε 2|F 2 → ε 2|F 12 2.
Interpolating the strain tensor at the contact interface by averaging both strain tensor :

ε 12|F 12 = 1 2 ε 1|F 12 + ε 2|F 12
3. Substituting the first diagonal component of ε 12|F 12 by the local engineering strain :

ε 12|F 12 =   ε xx ε xy ε xz ε yx ε yy ε yz ε zx ε zy ε zz   F 12 → ε xx = ∆L L 0
Thanks to these three main steps, the contact strain tensor ε 12|F 12 expressed in the local contact frame F 12 is deduced. It combines information from local contact (the ε xx component) and non-local information from surrounding discrete elements. Such tensor should be considered as both local and non-local. The next section will highlight the method for computing reaction force and torque from ε 12|F 12 .

Elastic brittle contact model

In this section, the mechanical contact model is described. The developed contact model is able to simulate cohesive forces between discrete elements compatible with a brittle elastic continuum description.

Cohesive force computation

The stress tensor σ 12|F 12 at the contact (see Figure 2) interface can be deduced from the strain tensor ε 12|F 12 thanks to the Hook's law as follows :

σ 12|F 12 = 2µ ε 12|F 12 + λ tr(ε 12|F 12 ) • I ( 20
)
where µ and λ are the Lamé parameters that come from elastic properties of the considered material. The contact reaction force f 12|F 12 , expressed in the local contact frame F 12 is deduced from the stress tensor σ 12|F 12 as follows :

f 12|F 12 = S 12 σ 12|F 12 • x 12 (21)
where S 12 is the contact surface between the discrete element 1 and 2. By considering the σ 12|F 12 as :

σ 12|F 12 =   σ xx τ xy τ xz τ xy σ yy τ yz τ xz τ yz σ zz   F 12 (22)
Expanding equation 21 gives :

f 12|F 12 = S 12   σ xx τ xy τ xz τ xy σ yy τ yz τ xz τ yz σ zz   F 12 •   1 0 0   F 12 = S ij   σ xx τ xy τ xz   F 12 (23)
If the related Cauchy's strain tensor ε 12|F 12 expressed in F 12 is :

ε 12|F 12 =   ε xx ε xy ε xz ε xy ε yy ε yz ε xz ε yz ε zz   F 12 (24)
By combining the last equations, it gives :

f 12|F 12 = S ij • E 1 + ν   ε xx + ν 1-2ν (ε xx + ε yy + ε zz ) ε xy ε xz   F 12 (25)
where E and ν are respectively the Young's modulus and Poisson's ratio of the considered elastic material. They can be deduced from the Lamé parameters with :

E = µ(3λ + 2µ) λ + µ (26) ν = λ 2(λ + µ) (27)
Thanks to equation 25 the contact reaction force are deduced from the Hook's law which is compatible within a perfect elastic continuum.

The last part of the model is related to the consideration of local moments, which is necessary to deal with large transformations. This is detailed in the next section.

Including moments

The cohesive contact force f 12 is considered to act at the contact interface point O 12 . The induced moment m 12 of this force acting on the discrete element 1 at its center O 1 can be simply deduced with the classical formula :

m 12 = ----→ O 1 O 12 ∧ f 12 (28)
and the dual moment m 21 , acting on the discrete element 2 at its center O 2 is :

m 21 = - ----→ O 2 O 12 ∧ f 12 (29)
Considering moments induces possible rotations of discrete elements. So, a rotational stiffness K 12 must be introduced to ensure the cohesion in rotation of the whole assembly. Inspired from micro-polar elasticity, this rotational stiffness is defined as :

K 12 = E × L 0 ( 30
)
where E is the Young's modulus and L 0 is the initial interaction length between discrete elements 1 and 2. From micro-polar point of view, L 12 can be considered as the internal length. Now, supposing the relative rotational vector θ 12 between the discrete element 1 and 2. This vector represents the relative rotation through the axis-angle formulation : the direction of θ is related to the rotational axis and the magnitude of θ is related to the angle in radian. This rotation vector can be deduced from quaternions associated to discrete elements able to describe rigid body rotations [START_REF] Diebel | Representing attitude : Euler angles, unit quaternions, and rotation vectors[END_REF].So, additional moments m 12 and m 21 that comes from this rotational stiffness K 12 must be added as :

m 12 =K 12 θ 12 (31)
m 21 = -K 12 θ 12 (32)
Finally, the total moments m 12 and m 21 that act on discrete elements 1 and 2 due to their cohesive elastic contact are :

m 12 = ----→ O 1 O 12 ∧ f 12 + K 12 θ 12 (33)
m 21 = - ----→ O 2 O 12 ∧ f 12 -K 12 θ 12 (34)

Crack management and brittle failure criterion

In comparison with FEM, a great advantage of DEM is the management of discontinuities. Here, a local criterion is proposed for modeling brittle fracture of elastic media. Fracture management with DEM is extensively used to study fragmentation processes ( 9), [START_REF] Watson | Simulating hypervelocity impact phenomena with discrete elements[END_REF], machining [START_REF] Tan | Simulation of ceramic grinding mechanism based on discrete element method[END_REF], microstructures [START_REF] André | A discrete element thermo-mechanical modelling of diffuse damage induced by thermal expan-sion mismatch of two-phase materials[END_REF] or damage processes [START_REF] Coré | Dynamic energy release rate evaluation of rapid crack propagation in discrete element analysis[END_REF]. The proposed criterion is based on the assumption that crack opening remains only for mode I. This assumption is relevant for brittle media. In this case, a crack surface corresponds to a breaking contact if the cohesive normal contact stress σ xx reaches a given value σ f . Considering a contact between the discrete elements 1 and 2, the contact normal stress σ xx can be computed as :

σ xx = σ 12|F 12 • x 12 • x 12 (35)
In fact, σ xx is the first diagonal component of the stress tensor σ 12|F 12 expressed in the F 12 frame. The failure criterion is simply :

σ xx > σ f ( 36
)
where σ f is the failure strength in uniaxial tension of the material. If a contact reaches this criterion, it vanishes. The related interaction is definitively broken and the contact should be considered as a classical one (non cohesive). In addition, the related broken contact should be excluded from the computation of local strain tensors because the crack is considered as a discontinuity in the displacement field.

Explicit numerical scheme

The numerical solution is based on an explicit integration scheme that is well suited to massive DEM simulations [START_REF] Rougier | Numerical comparison of some explicit time integration schemes used in dem, fem/dem and molecular dynamics[END_REF] and high velocity phenomena, such as fracturing or impact. The velocity Verlet scheme is chosen for its simplicity and efficiency. The discrete element orientations are described by quaternions, which provide an efficient way of computing the discrete element rotations (37, §2.5). The quaternions must be normalized at every time step to prevent numerical drift.

However, one observes that the proposed algorithm exhibits instabilities after thousands of iteration steps. During a computation, high frequency vibrations appear. Their amplitude increases, and they eventually lead to a divergence of the scheme and the simulation fails. To avoid this problem, a small numerical damping factor is introduced into the integration scheme as described by Tchamwa and Wielgosz [START_REF] Mahéo | Damping efficiency of the tchamwa-wielgosz explicit dissipative scheme under instantaneous loading conditions[END_REF]. This is a decentered explicit integration scheme that allows high frequencies dissipation. This scheme is very similar to the Velocity Verlet algorithm. The dissipation is controlled with a single parameter ϕ. Only the second time derivative equality is modified. The related Velocity Verlet equation becomes:

˙ p(t + ∆t) = ˙ p(t) + ϕ∆t ¨ p(t) + ¨ p(t + ∆t) ( 37
)
where p is the position vector, t is the time variable and ∆t is the time step. The case where ϕ = 0.5 is strictly equivalent to the Velocity Verlet integrator scheme.

For the whole study, values between 1 ϕ 1.3 are chosen. This range of ϕ ensures a good stability of the numerical scheme. The same approach is applied to quaternions :

q(t + ∆t) = q(t) + ϕ∆t (q(t) + q(t + ∆t)) (38)
where q is related to a quaternion. The same ϕ value is applied for both linear position vectors and angular quaternions. In addition, the explicit time step ∆t is chosen as three times lower than the critical time step :

∆t = min(∆t ij ) = min( 1 3
m ij L ij ES ij ) (39)
where:

• i and j are the two discrete elements in contact,

• ∆t ij is the critical time step for the ij contact,

• m ij is the equivalent mass of the contact :

m ij = m i × m j m i + m j (40)
• S ij is the contact surface and

• L ij is the initial length.

Virtual sample preparation

In order to simulate continuum with DEM, without voids, spherical discrete elements are proscribed because they prevent to fill completely and continuously a volume space. This limitation is tackled here by using polyhedral discrete elements. Polyhedral shapes allow to fill completely a volume. The most simple 3D pavement involves cubes or rectangular parallelepipeds as illustrated on figure 3. In these cases of regular grids, continuum problems can be approximated with finite difference methods. Indeed, these ordered arrangements lead to configurations where contact surfaces between elements are always oriented in same directions. As it was shown in (3), this kind of ordered arrangement conducts to geometrical contact anisotropy and promotes non-realistic crack paths. For this reason, ordered configurations are proscribed. Here, two stochastic methods of virtual sample preparation are tested.

The noisy grid method

The method used here was first proposed by (39) and used for example in [START_REF] Delaplace | Discrete modeling of cracking of brittle materials in large relative motion and localization problem[END_REF]. A regular 3D cartesian grid is built. Points are inserted at the center of each elementary cell of this grid. Then, a controlled noise is applied to the coordinates of these points. Finally, these points are used as germs for a Voronoi diagram that completely fills the space. Here, the Voro++1 library is used to implement Voronoi tessellations.

The noise applied to the point coordinates is set through uniform distribution laws. The range of the distribution, noted d, is chosen between 0% and 100%. As example, if 0% is chosen, each point is perfectly located at the center of its related cell. If 100% is chosen, each point can be everywhere in its related cell. Figure 4 illustrates this process for different distribution ranges in two dimensions. For d = 0%, the Voronoi tessellation is perfectly ordered and aligned with the original grid and all polyhedrons are cubic. Higher spatial distribution ranges give non cubic polyhedrons. Therefore, a question remains for choosing the right value of this dispersion range in 3D. The following paragraph will discuss briefly about this point. To study this problem, cumulative orientation histograms were built using the procedure described in [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. These histograms highlight the level of geometrical isotropy of packed domains. Generally, these histograms are given in 2D. However these 2D graphs give only information for one given plane and not in the whole 3D space. Here, the histograms are built in the 3D space. If a 3D histogram related to a given packing system is closed to a spherical shape, it means that the related contact orientation distribution is closed to a uniform distribution. In this case, there is no privileged contact direction and the discrete domain can be considered as geometrically isotropic. This is an important property for further simulations where realistic crack propagation paths are expected. Figure 5 shows 3D histograms for several dispersion ranges from 0% to 100%. One considers that domains with spatial dispersion higher than 90% give correct level of geometrical isotropy. An example of cubic domain built with this method for d = 90% is shown on Figure 6.

The packing method

Another strategy for building virtual samples involves classical sphere packing. It promotes isotropy of contact direction distribution especially for random distribution of sphere radii higher than 15% [START_REF] André | Discrete element method to simulate continuous material by using the cohesive beam model[END_REF]. The first step of the virtual sample preparation consists in filling a given volume with spherical discrete elements. The packing method used here comes from the GranOO discrete element free code 2 . It allows to control some packing parameters and ensures acceptable level of contact isotropy. More information about the involved packing algorithms and related 3D contact orientation histograms can be found in [START_REF] Jebahi | Discrete Element Method to Model 3D Continuous Materials[END_REF]. Finally, the obtained packed domain is processed with a Voronoi tesselation. The centers of discrete element are used as germs for the Voronoi tesselation. Again, the Voro++ library is used to implement these tessellations. The Voro++ library was chosen for its ability to deal with surface boundaries. Boundary primitives such as cylinders, spheres or boxes, have been already implemented in Voro++. The non-planar surfaces are approximated by plane sets able to cut Voronoi cells located at the boundary of the domain. In addition, the Voro++ library offers simple C++ interface for implementing custom boundary surfaces. Figure 7 shows the result of a sample preparation for a cylindrical domain.

Discussion about these two methods

The first sample preparation method involves only geometrical algorithms. This is a great advantage because it is very fast in comparison with the second one that involves a granular DEM computation. However, as it was shown in the previous section, the level of dispersion must be high, more than 90% to produce a uniform distribution of contact direction. Indeed, this high dispersion value conducts to configurations where two discrete elements can be very close together. These configurations promotes instability of the numerical scheme and must be avoided.

This situation is overcome and controlled with the second method that involves DEM packing processes. The minimal distances between discrete elements are related to their radii. It makes the distance distribution between discrete elements more uniform and avoid the situation described above. Indeed, some improvements may be done for the noisy grid method to avoid geometrical configurations that promotes numerical instabilities. In the following study, the packing method is preferred for building virtual samples.

Validation through virtual tests

This section proposes to validate the implemented model through numerical quasi-static, dynamic and failure tests. Indeed, the current implemented code is a prototype. Considerations about performance and optimization are not regarded here.

Qualitative test : the rolling ball test

The following test consists in dropping a deformable ball on an inclined planar ground. Due to gravity, the ball bounces and rotates along the inclined ground as shown on Figure 8. This figure shows the configuration at three temporal states : the initial state, an intermediate state and the final state after 50,000 iteration steps. Two trajectory paths are drawn. They correspond to the trajectory of the centers of two discrete elements (1) and ( 2) diametrically opposed.

The main parameters of this simulation are given in Table 1. At the initial state t = 0 ms, the ball is relaxed and does not touch the plane. An initial vertical velocity v (along the gravity axis) is imposed to the ball. The contact detection between the inclined plane and the polyhedral discrete elements involves the EPA-GJK algorithm [START_REF] Bergen | A fast and robust gjk implementation for collision detection of convex objects[END_REF]. After 50,000 iteration steps, the ball has bounced 4 times, has translated about 4 times its diameter and has rotated about one time. (2)

(2)

(1)

(1)

(2) This simulation is used here as a qualitative benchmark. In opinion of authors this test is of high level of interest. Behind its apparent simplicity, this test is complex because it involves high transformations: non negligible rigid body motion in both translation and rotation and large deformation at the contact interface (about 200,000 µdef). In such cases, the simulation gives qualitatively acceptable results only if :

1. local rotations of discrete element are taken into account, 2. a rotational stiffness K is introduced (see section 5.2), and 3. strains tensor are deduced and expressed through local frames.

This test highlights the main improvements of the proposed method. If one of the above ingredients is missing from the model, the related simulation quickly diverges or gives non-realistic results. In this new LSDM formulation, usage of local frames for computing local strain tensors enable simulations where large transformations are expected. Such transformation may be present in failure, impact and fragmentation phenomena which represent a class of problems where DEM models are strongly involved. In addition, the proposed model is compatible within classical DEM frameworks because local rotations are taken into account. These local rotations are of high level of importance especially for non spherical discrete elements. Again, if local rotations are not taken into account, the simulation of the rolling ball quickly diverges.

Quantitative tests

This section focuses on the accuracy of the implemented model regarding theoretical expected results. This is a common way in DEM community for validating models because DEM contact models act only at local scale. In this study, equality between local and apparent global behaviours are expected.

The material used as reference in the following quantitative tests is virtual. The main material constants of this virtual material are given in the table 2. Virtual quasi-static tests are implemented to monitor the accuracy of the implemented model in terms of apparent Young's modulus < E > and Poisson's ratio < ν >. Virtual tests were performed on cubic samples of one meter length discretized with around 11,000 polyhedral discrete elements (see figure 9) and around 72,000 contacts. Following the relationship given in equation 39, the time step is ∆t = 2.81×10 -6 seconds. The tensile test involves 6,000 time increments, corresponding to a total relative elongation of 1.2%.

Young's modulus E 20 GP a Poisson's ratio ν [-1, 0.5] - Density ρ 1000 kg.m -3 Failure strength σ f 10 M P a
In order to simulate a tensile test, a progressive elongation at constant velocity is applied to both opposite faces of the cubic sample (called normal faces on figure 9). The apparent normal stress < σ n > is computed by taking the mean of the total internal forces F x-and F x-that act on these two faces as follows :

< σ n > = 1 2S ( F x-+ F x+ ) (41)
where S is the section of the cube. The apparent Young's modulus < E > is given as the ratio of the apparent normal stress < σ n > over the imposed elongation. In addition, average displacements of the orthogonal faces of the cubic sample are monitored. It allows us to deduce the apparent contraction of the sample and, so on, to deduce the apparent Poisson's ratio < ν >.

An example of result for ν = 0.4 is shown on Figure 10. Numerical results show a quick convergence of the apparent parameters towards the local ones with less than 0.5% of difference. This good agreement is also observed for different values of Poisson's ratio. Table 3 summarizes the results of computations for Poisson's ratio values in the [0.5, -1] range. For extremum values of ν close to 0.5 and -1, the numerical scheme becomes unstable. The reason comes from the related physical elastic constant K for bulk modulus and G for shear modulus. These constants can be computed from E and ν as :

K = E 3(1 -2ν) (42)
G = E 2(1 + ν) (43)
One observes that the bulk modulus K tends to infinity when Poisson's ratio is close to 0.5 and the shear modulus G tends to infinity when Poisson's ratio is close to -1. Another observation concerns the accuracy of the method. A small loss of accuracy is observed for extreme Poisson's ratio values. Indeed, these values of Poisson's ratio are not really common for the wide range of engineering materials.

In addition, several samples are tested with different level of discretization from 2,000 to 20,000 discrete elements. All the simulations carried out with these samples give good agreements between apparent parameters and local input parameter values. In conclusion, the proposed model is able to simulate accurately an elastic continuum in quasi static configurations.

Dynamic tests

The aim of this section is to validate the model from the dynamic point of view. Dynamics plays an important role during cracking phenomena because of the sudden elastic energy released during crack extension. So, the proposed model must be relevant regarding dynamic effects. As described in section 6, the implemented numerical scheme dissipates high frequency vibrations. Regarding stability this is a good solution, indeed this damping must have negligible effects on the propagation of mechanical waves.

To study this question, this section implements two tests : a dynamic tensile test able to validate longitudinal waves propagation and a dynamic torsion test able to validate transverse mechanical waves propagation. The numerical sample used

Input parameters

Output parameters Bernoulli beams. The virtual sample has 1 meter length and 0.1 meter radius. It is discretized with 22,200 discrete elements and contains 146,500 contacts.

E (GPa) ν < E > (GPa) < ν > 20
For both tests the S x-face is clamped. It means that the discrete element of the S x-face are not able to move. At the initial state, the sample is placed in a loaded configuration that corresponds to the static solution of a tensile or a torsion test. In order to generate first mode oscillations, the sample is suddenly released (except the clamped face).

For the tensile test, the applied initial displacement field is : where :

u = Ax v = 0 w = 0 x y S x+ S x-
• u, v and w are the components of the displacement field along x, y and z,

• x, y and z are the related Cartesian coordinates and

• A is related to the imposed elongation along x.

For the torsion test, the initial applied displacement field is :

u = 0 v = y cos(α × x) -z sin(α × x) -y w = y sin(α × x) + z sin(α × x) -z
where α is the twist unit angle (in rad.m -1 ) related to the imposed shear strain. The material strength theory gives the theoretical frequency of the first mode of natural vibration for a clamped-free sample in a tensile configuration [START_REF] Shabana | Vibration of Discrete and Continuous Systems[END_REF] :

f tension = 1 4L E ρ ( 44
)
where L is the length of the sample. In addition, the frequency of the first mode of natural vibration for a clamped-free sample in a twisted configuration is given by [START_REF] Shabana | Vibration of Discrete and Continuous Systems[END_REF] :

f twist = 1 4L G ρ ( 45
)
where G is the shear modulus (see equation 43). In conclusion, the presented dynamic tests give results in very good accordance with theory. The developed model is able to accurately simulate propagation of mechanical waves in both longitudinal (compression) and transverse (torsion) modes. The numerical damping used to improve the stability of the numerical scheme seems to have negligible influence on the velocity of mechanical waves. However, as it is shown on figure 12, it has a significant effect on the amplitude of oscillations. This decreasing means a possible loss of energy. Indeed, for short time simulations such as impact tests that involve few oscillations, the observed dissipation is not so high and remains acceptable.

Tensile failure test

In order to validate brittle behaviour, tensile failure tests are carried out. These tests are very similar to those described in section 8.2.1. The difference remains only in the introduction of a local tensile failure strength equal to σ f = 10 MPa as described in section 5.3.

Figure 13 plots a typical stress-strain evolution of a tensile failure (see red line curve on this graph). The sudden decreasing in the apparent normal stress < σ n > marks the apparition of a macroscopic crack. Here, the macroscopic crack appears close to the introduced σ f value. Difference between the input value and the apparent failure strength is about 4%. In addition, the crack path, highlighted on Figure 14, is correctly oriented, orthogonal to the main normal stress direction. The next paragraph will discuss about this remainging difference of 4%. From the earlier work of Griffith [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF], it is well established in literature that brittle fracture occurs on local defects [START_REF] Lawn | Fracture of brittle solids -second edition[END_REF]. So, macroscopic cracks are governed at the microscopic scale by the local defects inside materials. Here, the proposed DEM model induces variability that comes from sample preparation method (see section 7). These uncertainties conduct to non-homogeneous stress field and failure occurs on the most stressed contact interface. This intrinsic uncertainty conducts to the observed difference between the local failure criterion and the apparent observed failure value.

In order to confirm this scenario, additional tensile failure tests are carried out on a perfect sample with cubic discrete elements ordered on a regular grid. Figure 15 shows the ε xx strain distribution during a tensile test for both randomized domain and regular domain. For the ordered case, the strain distribution ε xx is perfectly homogeneous. For the randomized case, the strain distribution exhibits dispersion. This dispersion has no significant effect on the apparent elastic properties such as Young's modulus, Poisson's ratio, or mechanical wave propagation, but have significant effect on apparent failure strength. This is confirmed by the blue dashed curve in Figure 13 that highlights that the perfect sample breaks at 10 MPa exactly. So, in the ordered case, no significant difference is observed between local fracture criterion and apparent failure strength. Indeed, no local-ization occurs for the perfect case and the observed crack path is not realistic. In conclusion, random configurations are preferred even if the local and apparent failure strength does not match perfectly. The indirect tensile failure test, also named Brazilian test, is widely used for characterization of brittle materials such as ceramics, concretes or rocks. The Brazilian test is performed by applying a vertical compressive load across the diameter of a disk sample. In the fracture phase, a crack appears along the vertical diameter of the disk due to tensile stresses induced horizontally by the geometry of sample [START_REF] Li | The brazilian disc test for rock mechanics applications: review and new insights[END_REF].

Figure 16 shows the numerical experiment. A displacement imposed at constant velocity is applied to the discrete elements highlighted on figure 16a diametrically and vertically opposed. The main parameters of the simulation are given on Table 5.

The blue solid line on Figure 17 shows the evolution of the applied force versus iteration number. Again, a sudden decrease marks the apparent failure of the specimen. Theoretically, the failure strength force P f can be deduced from the analytic expression (45):

P f = π 2 × σ f × D × t (46)
With the values extracted from Table 5, the theoretical specimen strength is P = 64.2 kN . The difference with the numerical results is less than 1.5%. However, as shown on Figure 16b, failure does not start at the center of the disk. According to experimental and numerical studies, crack initiation at the center of the disk sample is considered crucial for the test validity [START_REF] Wang | The flattened brazilian disc specimen used for testing elastic modulus, tensile strength and fracture toughness of brittle rocks: analytical and numerical results[END_REF][START_REF] Fairhurst | On the validity of the 'brazilian' test for brittle materials[END_REF][START_REF] Erarslan | Experimental, numerical and analytical studies on tensile strength of rocks[END_REF]. Fairhurst stated that "failure may occur away from the center of the disk for small angles of loading contact area (...). In such cases the tensile strength as usually calculated from test results, is lower than the true value" [START_REF] Fairhurst | On the validity of the 'brazilian' test for brittle materials[END_REF]. According to Nazife Erarslan et al. [START_REF] Erarslan | Experimental, numerical and analytical studies on tensile strength of rocks[END_REF], "a 20 • loading arc gives the best estimate of the indirect tensile strength".

According to these observations, a second simulation with a loading arc of 2α = 20 • is carried out (see Figure 18a, the loading arc is highlighted in yellow). The red dashed curve on Figure 17 shows the evolution of the applied force versus iteration number. As expected, the strength of the specimen is higher. When the loading arc becomes non-negligible, the theoretical strength P is modified by a factor k (49, 50): k = (2 cos 3 α + cos α + sin α/α) 2 8(cos α + sin α/α) × α sin α [START_REF] Fairhurst | On the validity of the 'brazilian' test for brittle materials[END_REF] where α is the half contact angle. The modified strength P is given by :

P = 1 k × π 2 × σ f × D × t (48)
With the values extracted from Table 5, the specimen strength is P = 66.8 kN for a loading arc of 2α = 20 • . In this case, the difference with theory is about 3%. Another observation concerns the location of fracture. As expected, the failure occurs at the center of the disk (see Figure 18b). This last observation is relevant with the previously theoretical and experimental cited studies.

(a) initial state (b) cracks at it=18,400 (c) cracks at it=18,800

Perspectives and conclusions

The presented paper describes a novel DEM numerical method for modeling brittle elastic media. The proposed model is directly described in terms of continuum mechanics. This is achieved by proposing a method for computing local strain tensor based on a local linear displacement field assumption and a least square method. The second part of the presented study shows that the model is relevant with material strength theory in quasi-static and dynamic configurations. Indeed, some issues were observed regarding numerical stability of the proposed numerical scheme. A solution was proposed by introducing a numerical damping in the integration scheme able to stabilize simulations. In addition, dynamic tests indicate that the introduced damping factor has negligible effects regarding mechanical wave velocities. A solution, for fixing this issue, may be the introduction of a mechanical damping factor which can be related to physical material properties. Finally, a simple failure criterion was introduced and was able to simulate crack initiation and propagation of brittle media.

To conclude, this contribution may enhance the dissemination of DEM in the applied mechanics and engineering community for studying brittle elastic media such as rocks, glasses or ceramics. Mainly, material behaviours are expressed in terms of continuum mechanics, i.e, stress-strain behaviour laws. A common issue with DEM for modeling a continuum is the difficulty to introduce directly physical quantities into DEM contact laws. This limitation was tackled here for brittle elastic media. The proposed contact model comes directly from Hooke's law and model the cohesion of an elastic media. Further developments may concern more complex behaviour such as plasticity or anisotropy. This approach may also be used for introducing more refined failure criterion such as the Drucker-Prager yield criterion commonly used in rock mechanics and civil engineering.

Bibliography

Figure 1 :

 1 Figure 1: Configuration of problem in 2D

Figure 2 :

 2 Figure 2: Configuration of the contact problem problem in 2D

Figure 3 :

 3 Figure 3: 3D discrete domain with parallelepipedic discrete elements

  (a) result with d = 0% (b) result with d = 50% (c) result with d = 100%

Figure 4 :

 4 Figure 4: 2D Voronoi tessellations for different values of spatial distribution range d

Figure 5 :

 5 Figure 5: 3D histograms of surface normal for different values of spatial distribution range d

Figure 6 :

 6 Figure 6: 3D Voronoi discrete domain with d = 90%

Figure 7 :

 7 Figure 7: Example of a cylindrical virtual sample

20 t

 20

Figure 8 :

 8 Figure 8: The rolling ball test

Figure 9 :

 9 Figure 9: Configuration of virtual tensile test

Figure 10 :

 10 Figure 10: Apparent Poisson's ratio and Young's modulus deduced from virtual tensile test

Figure 11 :

 11 Figure 11: Cylindrical sample

Figure 12 :

 12 Figure 12: Dynamic response of the free face for the tensile configuration

Figure 13 :

 13 Figure 13: Evolution of apparent normal stress versus relative elongation strain

Figure 14 :

 14 Figure 14: Crack path generated during a tensile test (displacement ×50)

3 1

 3 × 10 -3 (b) Perfectly ordered domain

Figure 15 :

 15 Figure 15: ε xx strain distribution in a tensile test

Figure 17 :

 17 Figure 17: Evolution of force versus iteration

  (a) initial state (b) cracks at it=22,800 (c) cracks at it=23,300

Figure 18 :

 18 Figure 18: Indirect tensile test for 2α = 20 •

[ 1 ]

 1 P. A. Cundall, O. D. L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1979) 47-65. doi:10.1680/geot.1979. 29.1.47.

Table 1 :

 1 Main simulation parameters of the rolling ball test

	Discrete element number n 2,000	-
	Young's modulus E	0.01	GP a
	Poisson's ratio ν	0.1	-
	Density ρ	1000 kg.m -3
	Time step ∆t	60	µs
	Ground inclination α	10	•
	Contact stiffness Ks 10 6	N.m -1
	Contact friction coeff µ	0.4	-
	Initial velocity v	5	m.s -1

Table 2 :

 2 Main constant values of the reference model material

8.2.1. Quasi-static tensile test

Table 3 :

 3 Elastic tensile tests results for different values of Poisson's ratio

	y
		normal faces
	S x-	imposed elongation
		x
	orthogonal faces	S x+

Table 4

 4 summarizes the main values for these two dynamic studies. During the simulations, linear and angular displacements of the free S x+ face are monitored. For the tensile test, the result is shown on Figure12. It clearly shows that the free face S x+ exhibits oscillations. To deduce the related oscillation frequency, a Fast Fourier Transform (FFT) is applied to this temporal signal. The result of the FFT gives a peak at 1105 Hz. The theoretical value deduced from equation 44 gives 1118Hz. Difference between numerical and theoretical results is very close, less than 0.12%. The same procedure is applied for the twisted test. Again, theoretical (668 Hz) and numerical results (670 Hz) are very close, less than 0.1%.

	Initial tensile strain	K	10 -5	def
	Initial unit angle twist	α	10 -4	rad.m -1
	Sample length	L	1	m
	Sample radius	R	0.1	m
	Theoretical freq. in tension f tension	1118	Hz
	Theoretical freq. in torsion f twist	668	Hz
	Time step	∆t	6.96 10 -7	s
	Total iteration number	-	50 000	-
	Simulation duration	-	10 hours	h

Table 4 :

 4 The main parameter values of dynamic tensile and torsion tests

Table 5 :

 5 The main parameter values of indirect tensile test

	Figure 16: Indirect tensile test

http://math.lbl.gov/voro++/

www.granoo.org

LaTeX Source Files Click here to download LaTeX Source Files: article.zip