Guilhem Mollon
email: guilhem.mollon@insa-lyon.fr

Adriana Quacquarelli

Edward Andò

Gioacchino Cinno

Adriana Quacquarelli1

• Edward Andò

Gioacchino Viggiani

Can friction replace roughness in the numerical simulation of granular materials?

Keywords: DEM, Friction, Roughness, Shearing, Localization

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Granular simulations based on various implementations of the Discrete Element Method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] are nowadays commonly used in the various fields of science and technology interested in divided solid matter. This includes civil engineering [START_REF] Donze | Advances in discrete element method applied to soil, rock and concrete mechanics[END_REF], geomechanics [START_REF] Darve | DEM modelling in geomechanics: some recent breakthroughs[END_REF], environmental sciences [START_REF] Richards | Discreteelement modelling: methods and applications in the environmental sciences[END_REF], geophysics [START_REF] Mollon | Discrete modelling of rock avalanches: sensitivity to block and slope geometries[END_REF], granular physics [START_REF] Azema | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF], pharmaceutical processes [START_REF] Kodam | Discrete element mothod modeling of bi-convex pharmaceutical tablets: contact detection algorithms and validation[END_REF] and tribology [START_REF] Fillot | A granular dynamic model for the degradation of material[END_REF][START_REF] Mollon | A numerical framework for discrete modelling of friction and wear using Voronoi polyhedrons[END_REF], among others. While early developments of DEM were only considering discs or spheres (depending on the 2D or 3D context of the simulation), there is now a growing interest in the introduction of more complex, and maybe more realistic shapes. A part of this effort is dedicated to a better characterization of the shapes of various populations of grains, and to the extraction of relevant shape features using various experimental and analytical procedures [START_REF] Matsushima | 3D shape characterization and image-based DEM simulation of the Lunar soil simulant FJS-1[END_REF][START_REF] Blott | Particle shape: a review and new methods of characterization and classification[END_REF][START_REF] Wei | Generation of realistic sand particles with fractal nature using an improved spherical harmonic analysis[END_REF]. Another part, closely related to the first one, is dedicated to the generation of granular shapes and samples in order to reproduce morphological features of target grains populations [START_REF] Mollon | Fourier-Voronoi-based generation of realistic samples for discrete modelling of granular materials[END_REF][START_REF] Mollon | 3D generation of realistic granular samples based on random fields theory and Fourier shape descriptors[END_REF][START_REF] Ouhbi | 3D particle shape modelling and optimazation through proper orthogonal decomposition[END_REF]. Finally, the last part of this process is dedicated to the development of numerical techniques making it possible to introduce and simulate such generated shapes in DEM codes [START_REF] Azema | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF][START_REF] Lin | A three-dimensional discrete element model using arrays of ellipsoids[END_REF][START_REF] Richefeu | Dissipative contacts and realistic block shapes for modelling rock avalanches[END_REF][START_REF] Houlsby | Potential particles: a method for modelling non-circular particles in DEM[END_REF][START_REF] Kawamoto | All you need is shape: predicting shear banding in sand with LS-DEM[END_REF].

Meanwhile, an important question remains unanswered, which is related to the relevance of the geometrical features that are successively extracted from the real grains, used to generate synthetic grains, and simulated in a DEM framework. By relevance, we mean here relevance with respect to a certain simulation purpose. Indeed, each simulation is performed for a specific reason and with specific expectations regarding the output quantities to be extracted from it, be it in an academic or a practical framework. This question arises because, to date, a very large number of morphological descriptors were proposed based on real grains shapes analyses, but it is still very difficult for the practitioner to determine what should be the key morphological parameters to introduce in his/her simulation in order to reproduce some desired property of a given target granular sample (either mechanical, or thermal, electrical, etc.). Besides, each class of morphological properties comes with an associated numerical cost. For example, it is much less costly to simulate ellipsoids [START_REF] Lin | A three-dimensional discrete element model using arrays of ellipsoids[END_REF] than spheropolyhedra [START_REF] Richefeu | Dissipative contacts and realistic block shapes for modelling rock avalanches[END_REF], while both approaches are able to introduce an elongation of the grains. Hence, if the only relevant physical content needed for the success of a given simulation lies in the elongation, ellipsoids should be preferred. But how to be sure that this single property is sufficient in a given framework?

In the present paper, we only address a small portion of this important question. More specifically, we focus our attention on the notions of roughness and of friction, with respect to strain localization in sheared granular materials [START_REF] Desrues | Strain localization in sand: an overview of the experimental results obtained in Grenoble using stereophotogrammetry[END_REF]. We wish to evaluate to what extent the roughness of the surface of the grains (which requires a fine discretization of their frontier and thus implies a large computational cost) may be replaced by an artificial increase of their friction coefficient (which comes at no additional cost). Two concepts need to be clearly defined in the context of this study: friction and roughness.

Friction, on one hand, is a very complicated topic, and tribology is not a unified science yet. Establishing a direct link between roughness and friction has been tried many times (see for example [START_REF] Greenwood | Contact of nominally flat surfaces[END_REF][START_REF] Pöschel | Static friction phenomena in granular materials: coulomb law versus particle geometry[END_REF][START_REF] Pöschel | A simple geometrical model for solid friction[END_REF][START_REF] Ivkovich | The influence of the contact surface roughness on the static friction coefficient[END_REF]) and proved useful to enhance our understanding of contact mechanics. But the current state of the art is that there is no straightforward causality between these two concepts. Friction is now regarded, not as an intrinsic property of the bodies into contact, but as an emerging phenomenon [START_REF] Godet | The third-body approach: a mechanical view of wear[END_REF][START_REF] Colas | Describing third body flows to solve dry lubrication issue-MoS 2 case study under ultrahigh vacuum[END_REF] related to several physical properties and phenomena (surface adhesion, heat production, wear, acoustic wave emissions, local phase change, microstructural evolutions, interactions with the gaseous environment, etc.) occurring at several scales (the scale of the contacting bodies, the scale of the wear debris filling the contact interface, the scale of surface topography, the scale of surface chemistry, the scale of the mechanical system driving the motion of the contacting bodies, etc.). Roughness only plays a small part in this plot. With this complexity in mind, it is important to distinguish two concepts: (i) the complex and emerging physical phenomenon that we call friction, and (ii) the phenomenological quantity that is introduced in a granular simulation in order to prescribe satisfactory values to the force resisting to sliding motion. We clearly place our work in the context of this second definition of the friction. Since the physical complexity of a mechanical contact is out of reach in a DEM simulation, we mostly consider our coefficient of friction as a contact parameter leading to satisfactory granular behaviors. As such it is not different at all from those commonly used in DEM, where this shortcut is always tacitly used.

Roughness, on the other hand, seems to be rather illdefined in the granular mechanics literature. There are several occurrences [START_REF] Jensen | DEM simulation of granular media-structure interface: effects of surface roughness and particle shape[END_REF][START_REF] Jensen | Effect of particle shape on interface behavior of DEM-simulated granular materials[END_REF][START_REF] Kozicki | Numerical simulations of sand behavior using DEM with two different descriptions of grain roughness[END_REF][START_REF] Kozicki | Effect of grain roughness on strength, volume changes, elastic and dissipated energies during quasi-static homogeneoue triaxial compression using DEM[END_REF][START_REF] Hare | The influence of aspect ratio and roughness on flowability[END_REF][START_REF] Mede | A medial axis based method for irregular grain shape representation in DEM simulations[END_REF] into which this word is used to describe the whole shape of the grain, at every scales, especially in DEM simulations using clumps, i.e. rigid clusters of discs. In some other papers, roughness is restricted to the topography of grains at very small scales, typically below the micrometer [START_REF] Santamarina | Effect of surface roughness on wave propagation parameters[END_REF][START_REF] Alshibli | Characterizing surface roughness and shape of sands using digital microscopy[END_REF][START_REF] Otsubo | Experimental and DEM assessment of the stress-dependency of surface roughness effects on shearmodulus[END_REF][START_REF] Nadimi | Effect of particle roughness on the bulk deformation using coupled boundary element and discrete element methods[END_REF]. Such a usage is common in scientific communities close to tribology or contact mechanics, where the very notion of "shape" is irrelevant at the scale of the contact. In this case, contacting solids are seen as planes with a superimposed irregularity, and roughness is used to describe this irregularity. Between these two extremes, it is often needed to consider in a single study the irregularities of the contacting objects at several scales: the scale of the object and the scale of its surface [START_REF] Gallier | Rheology of suspensions of rough frictional particles[END_REF][START_REF] Zheng | Traditional soil particle sphericity, roundness and surface roughness by computational geometry[END_REF][START_REF] Zheng | A corner preserving algorithm for realistic DEM soil particle generation[END_REF][START_REF] Wilson | The influence of surface roughness and adhesion on particle rolling[END_REF][START_REF] Hsu | Roughness-dependent tribology effects on discontinuous shear thicknening[END_REF]. A common approach for sand grains [START_REF] Blott | Particle shape: a review and new methods of characterization and classification[END_REF][START_REF] Jensen | Effect of particle shape on interface behavior of DEM-simulated granular materials[END_REF][START_REF] Zheng | Traditional soil particle sphericity, roundness and surface roughness by computational geometry[END_REF][START_REF] Zheng | A corner preserving algorithm for realistic DEM soil particle generation[END_REF][START_REF] Alshibli | Characterizing surface roughness and shape of sands using digital microscopy[END_REF][START_REF] Ehrlich | An exact method for characterization of grain shape[END_REF][START_REF] Bowman | Particle shape characterization using Fourier descriptor analysis[END_REF][START_REF] Andrade | Granular element method for computational particle mechanics[END_REF][START_REF] Chen | Numerical study of particle morphology effect on the angle of repose for coarse assemblies using DEM[END_REF] is to consider three independent scales of shape description: circularity/sphericity (main deviation of the object shape from a perfect disc/sphere), angularity/roundness (related to the sharpness of the corners and the variability of local curvatures around the solid frontiers), and roughness (related to irregularities in the surface). There is no general consensus, however, on the exact scale limits that separate these three descriptors, especially because the second one is ill-defined (what should be called a "corner" was never clearly stated on a geometrical basis). From a more mathematical point of view, self-affinity of surface roughness has been demonstrated in many studies [43, 46-49, 53, 54], for any kind of natural surface of polycrystalline materials (it often breaks down for glassy materials or for artificial machined surfaces, for different reasons). Hence, a rigorous way to define the scale where irregularities should be called "roughness" instead of "shape" could be: the scale beyond which the spectrum of the amplitudes of these irregularities becomes self-affine (i.e. follows a decreasing power-law). This is the definition that was used in several studies [START_REF] Mollon | Fourier-Voronoi-based generation of realistic samples for discrete modelling of granular materials[END_REF][START_REF] Mollon | 3D generation of realistic granular samples based on random fields theory and Fourier shape descriptors[END_REF][START_REF] Ehrlich | An exact method for characterization of grain shape[END_REF][START_REF] Bowman | Particle shape characterization using Fourier descriptor analysis[END_REF][START_REF] Meloy | Fast Fourier transforms applied to shape analysis of particle silhouettes to obtain morphological data[END_REF][START_REF] Das | Modeling three-dimensional shape of sand grains using discrete element method[END_REF][START_REF] Mollon | Generating realistic 3D sand particles using Fourier descriptors[END_REF], and that we use here.

To investigate whether roughness can be replaced by friction in simulations of granular materials, we perform 2D granular simulations reproducing some published experimental results [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF], using an appropriate numerical framework into which it is possible to introduce complex shapes and realistic boundary conditions and to observe the appearance of strain localization. Two numerical campaigns are performed, the first one paying attention to the influence of surface roughness on the localization phenomenon, and the second one focusing on the consequences of an artificial replacement of the surface roughness by an increase in the friction coefficient.

Methodology

Description of the reference experiments

The simulations presented in this paper are inspired by an experimental campaign performed some years ago at Laboratoire 3SR in Grenoble, and described in details in [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF][START_REF] Hall | Discrete and continuum analysis of localised deformation in sand using X-ray μCT and volumetric digital image correlation[END_REF][START_REF] Andò | Grainscale experimental investigation of localised deformation in sand: a discrete particle tracking approach[END_REF] and in Fig. 1. In these experiments, sand samples underwent triaxial compression tests within a lab-size µCT scanner in order to follow at the grain scale the deformation and failure process in situ. The granular sample was a cylinder of diameter 11 mm and height 22 mm, containing many tens of thousands of sand grains. An appropriate amount of sand was first carefully deposited by dry pluviation inside a flexible membrane stretched within a mold and between two porous stones, and suction was then applied inside the membrane in order to create a confining pressure and to give the ability to the sample to hold steady. The mold was then removed, the geometric quality of the sample was checked, and the sample was positioned in the triaxial cell (which had to attenuate few x-rays) inside the scanner and submitted to a confining pressure applied on the membrane by pressurized water (vacuum in the sample was meanwhile released).

The loading was applied from below by the hemispherical extremity of a ram, in a strain-controlled manner with a prescribed loading rate of 21 µm/min (Fig. 1). This loading was regularly stopped in order to perform a full 3D scan of the sample. From these experimental developments and the associated postprocessing of the 3D digital images, a large number of statements were made on the micromechanics of strain localization at the grain scale, and the reader is redirected to the large body of published work for more details on the technical challenges that had to be addressed and on the mechanical interpretations that could be derived [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF][START_REF] Hall | Discrete and continuum analysis of localised deformation in sand using X-ray μCT and volumetric digital image correlation[END_REF][START_REF] Andò | Grainscale experimental investigation of localised deformation in sand: a discrete particle tracking approach[END_REF].

Samples characterization and generation

During the experimental campaign, several types of sands were tested and analyzed, but in the present paper we will only focus on one specific case, namely Ottawa sand. This is a sand of intermediate angularity, coming from sedimentary deposits, made up of quartz grains, with a D 50 of the order of 250 µm and a narrow grain size distribution (sieved between 210 and 300 µm). An SEM view of typical Ottawa sand grains is provided in Fig. 2a. In order to extract geometrical features, six Ottawa sand grains were scanned in a higher resolution machine at 1 µm/px. The voxels corresponding to the external surface of each grain were identified (their number was around 200,000 with the chosen scanning resolution), and these surfaces were then interpolated between these voxels on a set of 10,242 radial directions forming a geodesic structure, as described in more details in [START_REF] Mollon | Generating realistic 3D sand particles using Fourier descriptors[END_REF].

The resulting discretised surfaces of the six grains are shown in Fig. 2b, and are now defined as a discrete set of 10,242 radii R i i , i . Since the angular direction of each of these radii is well-known, it is possible to compute the angle they form with respect to each other, and thus to compute the angular autocorrelation function of the discrete random field R i . As explained in [START_REF] Mollon | 3D generation of realistic granular samples based on random fields theory and Fourier shape descriptors[END_REF], the Fourier spectrum of this field can then be computed in a univocal way, for each grain. The six spectra are plotted in Fig. 2 (their first modes in a Cartesian frame, and the distribution tails in a log-log frame), along with the mean spectrum obtained by averaging the six cases. There seems to be a good repeatability between the six spectra, indicating that their average bears some statistical significance despite the very limited size of the sample. The log-log chart indicates a power-law decrease of the spectra amplitude with their mode (up to mode 48, beyond which the noise created by the interpolation on the geodesic structure disturbs the high-frequency measurements), which is very consistent with the classical self-affine decay of roughness observed on sand grains [START_REF] Das | Modeling three-dimensional shape of sand grains using discrete element method[END_REF] and more generally on any natural surface at any scale [START_REF] Candela | Roughness of fault surfaces over nine decades of length scales[END_REF].

To improve the quality of the samples to be generated with this spectrum, the modes from 8 to 64 were replaced by a smooth power-law decay fitted on the modes 8 to 48 of the measured average spectrum, as proposed in [START_REF] Mollon | Fourier-Voronoi-based generation of realistic samples for discrete modelling of granular materials[END_REF]. Modes beyond 64 (i.e., corresponding to an angular period lower than 5.6°, and to a spatial distance smaller than about 12 µm for the average grain size considered here) were ignored, although it may theoretically be possible to extrapolate the power-law roughness decay to much smaller scales. 2D samples can then be generated based on this spectrum, using the Packing2D code described in [START_REF] Mollon | Fourier-Voronoi-based generation of realistic samples for discrete modelling of granular materials[END_REF], see Sect. 2.4.

Numerical framework

The simulations described in the next sections of this paper were performed with the code MELODY (described in [START_REF] Mollon | A multibody meshfree strategy for the simulation of highly deformable granular materials[END_REF][START_REF] Mollon | A unified numerical framework for rigid and compliant granular materials[END_REF] and available for download at http://guilh em.mollo n.free.fr). This code is based on a multibody meshfree approach, which is an extension of the classical DEM to the case of highly deformable grains. This framework is able to deal in a unified manner with a large number of rigid and deformable bodies, with arbitrary geometries. Both kinds of bodies are represented by polygons (i.e., by a set of boundary nodes linked by segments), and any kind of contact interaction between such two bodies (or self-contact of a deformable body) is handled by a robust penalty-based two-ways node-to-segment contact algorithm. The bulk of deformable bodies is discretized by a certain number of field nodes which carry the degrees-of-freedom in displacement, and the displacement field is interpolated between these nodes using Moving-Least-Square meshfree shape functions. The weak form so-obtained is then integrated in space (based on a Gaussian quadrature on a set of triangles discretizing the surface of each body), and in time (using an explicit adaptive time-stepping scheme). Newtonian dynamics of rigid bodies are also solved in time using the same solver, in a manner similar to classical DEM. The code is written in C++ and parallelized in an OpenMP framework, with a very good performance scaling up to 32 CPUs. A time-step close to 10 -7 s is reached by the adaptive solver for all the simulations presented hereafter.

Numerical simulations, preprocessing and load phasing

Since MELODY is currently restricted to plane-strain kinematics, we present here 2D simulations which aim to reproduce some features of the real 3D triaxial tests described in Sect. 2.1. Obviously, a large part of the physical reality is lost when dropping one dimension. A 2D sample is much more constrained from a kinematic point of view. However, 2D DEM is still used a lot in current scientific literature because some major aspects of the granular mechanical response (peak-plateau behavior, dilatancy, interlocking, etc.) appear to be well-reproduced in 2D. However, apart from implementation complexity and computational cost, nothing seems to prevent future extension of this work to 3D. The samples are composed of about 9600 grains, with a width of 22 mm and a height of about 45 mm after compaction. These samples are thus twice as large as the experimental one, but contain many fewer grains because they are in 2D. In the reference case, the full spectrum of Ottawa sand, measured and computed in Sect. 2.2, is used up to mode 64 in order to generate 2D grains with a contour discretized by 128 nodes with equal angular spacing. Although this is only an approximation, we will consider in the remainder of this paper that this case contains the full roughness of the grains (Fig. 3). However, since we are primarily interested in the effect of the roughness, we also generate several other granular samples with different properties. First, we introduce a smoothing of the surfaces of the grains, by cutting off the high angular frequencies of the Fourier spectrum used for their generation. This is done by considering only the N m first modes of the spectrum, with N m being equal respec- tively to 64 (reference case), 32, 16, and 8. Examples of such grains are provided in the first row of Fig. 3, where the smoothing effect of this operation is evident. This change in the spectrum, however, does not decrease the computational cost of the simulation because the number of nodes N n on the contour of each grain remains equal to 128. Hence, the number of proximity and contact detections and the number of contact resolutions at each time step do not change much.

In order to investigate the effect of this parameter, several samples are generated with smaller numbers N n of nodes on the contour (namely 128 in the reference case, 64, 32, and 16). The number of nodes is kept larger than half the number of modes in order to avoid artificial disappearance of some amount of roughness by removing nodes while it is still included in the considered spectrum. The list of the 10 generated samples is provided in Fig. 3. As clearly appears in this figure, decreasing N m reduces the roughness of the surface, but decreasing N n seems to introduce a somewhat artificial angularity because of the coarseness of the description of the grains contours.

A simulation is divided into a number of successive stages (Fig. 4) which allow a controlled and reproducible loading of the granular sample. At stage 0, the granular sample is generated using a Fourier-Voronoi technique, in such a way that the grains do not touch each other [START_REF] Mollon | Fourier-Voronoi-based generation of realistic samples for discrete modelling of granular materials[END_REF]. This sample is bounded by two porous stones (rigid rectangular bodies located above and below the sample) and the two walls of the mold (rigid rectangular bodies located on the left-and right-hand side of the sample). These four bodies are fixed in all directions, except for the lower porous stone which is allowed to move vertically, and which is in contact with the semi-circular extremity of the ram on its lower side.

During stage 1, an upward vertical force is applied to the ram, corresponding to a pressure of 100 kPa when distributed on the width of the lower porous stone. During this uniaxial compression stage (which mimics the pluviation process but with a better control on the obtained solid fraction), the coefficient of friction between the grains (and also between the rigid walls and the grains) is set to a certain value c , which controls the compacted solid fraction (the larger this friction coefficient, the looser the compacted sample). For all the simulations presented in this paper, this friction value is calibrated by trial and error in order to obtain a constant solid fraction in the final sample (i.e., before stage 4). Stage 1 lasts until the sample is compacted and stabilized.

At the beginning of stage 2, the flexible membranes are activated in the simulation. These membranes are slender rectangles (0.3 mm thick), initially located slightly remote from the granular sample. They are deformable, and follow a hyperelastic Neo-Hookean constitutive law with a Young's modulus of 1 MPa and a Poisson coefficient of 0.49 (in agreement with the properties of the latex used experimentally [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF]). In order for them to transmit the load applied by a virtual fluid located in the pressurized chamber around them, they are submitted on their outer edge to a pressure of 100 kPa. This pressure is applied in a fluidlike manner, meaning that, at each border node of the outer edge of a given membrane, the outwards normal vector is computed at each time step, and the fluid pressure is applied inwards in that specific direction. Since the membranes are very flexible, they quickly come into contact with the grains Fig. 3 Illustrative 2D Ottawa sand grain submitted to smoothing by a reduction of the number of Fourier modes and of the number of contour nodes (no contact detection is activated between the mold and the membrane), and transmit the external pressure to the sample in a rather realistic way (as was done in [START_REF] Kawamoto | All you need is shape: predicting shear banding in sand with LS-DEM[END_REF][START_REF] De Bono | DEM of triaxial tests on crushable sand[END_REF][START_REF] Yang | Particle-scale mechanisms in undrained triaxial compression of biocemented sands: insights from 3D DEM simulations with flexible boundary[END_REF]). The coefficient of friction between the membranes and the grains is set at c during this stage.

Stage 3 corresponds to the removal of the mold. The lower porous stone is completely freed from its kinematic constraints, the coefficient of friction in the whole sample is set to a new value d (which will be the one used for the deviatoric loading). The two lateral walls representing the mold are removed by giving them an outwards horizontal velocity. This means that the load from the grains is progressively taken up by the flexible membranes. At the end of this stage the sample is in an equilibrium state, the only fixed bodies being the upper porous stone (in all directions) and the loading ram (in all directions but the vertical one, where it is still submitted to an upwards force).

Stage 4 corresponds to the deviatoric loading. In this stage, the loading ram becomes displacement-controlled, and moves upwards with a velocity of 0.01 m/s. This velocity is much higher than the experimental one, but this is necessary in order to be able to run the simulation in a reasonable amount of time. Careful calibration was performed in order to ensure that this loading rate did not introduce any inertial effects in the simulations (by checking that the norms of the total vertical contact forces transmitted by the lower and upper boundaries are very close along time, indicating quasi-static loading). The coefficient of friction between the loading ram and the lower porous stone is kept low (0.05) during the whole stage 4 in order to allow slip at this interface, since this feature was observed experimentally. Stage 4 is performed until a strain level of about 15% is reached, similar to the actual experiment. During this stage, all the forces applied to the different parts of the apparatus, and all the kinematics of the grains are recorded. We observe that, during simulations, the membrane is often strongly folded (Fig. 4), requiring a self-contact detection in order to avoid breakdown of the simulation.

Postprocessing techniques

In order to analyze the mechanical behavior of the sample under deviatoric loading, a specific postprocessing technique was implemented. It consists in mapping a given state of the sample (i.e., the geometric positions of the 9600 grains) on a pixel-based "image", with a size of 3500 * 4500 pixels. A typical grain has a dimension of about 30 pixels in such image. Geometric detection techniques are used to attribute Fig. 4 Description of the different stages of a simulation, from sample preparation to deviatoric loading the value "1" to each pixel with a center located inside a sand grain, and "0" otherwise. This forms the image A (Fig. 5). Then, an image B is created by "closing" image A (i.e. by first dilating the image A by 20 pixels and then eroding the result by 20 pixels). This process (dilation-erosion) is very common in image processing in order to get rid of fortuitous holes into homogeneous regions [START_REF] Russ | The Image Processing Handbook[END_REF]. The resulting image B provides the spatial extension of the whole sample without any remaining porosity inside it.

An image C is created by convoluting image A with a Gaussian filter with a radius equal to 5R ( R being the aver- age sand grain radius). The consequence of this operation is to smooth the spatial of the matter in the image, with a spatial smoothing scale of the order of a representative elementary volume (REV) of about 15-25 grains. As shown in Fig. 5, the obtained image allows the observation of some local variations in the sample density, but is very inaccurate on the boundaries of the sample. More specifically, all the areas located inside the sample but within a distance of 5R from its boundary are affected by the fact that a portion of the Gaussian filter covers a region where there are no grains at all (and they are thus attributed a density smaller than it should be). To address this issue, an image D is created by applying the same Gaussian filtering to image B. The final image E is obtained after correcting image C by image D (i.e. dividing the value of each pixel of C by its value in D) and multiplying the result by image B (in order to cancel the density in the pixels located outside of the sample). The result, provided in Fig. 5, gives a very clear picture of the field of granular density in the sample.

From image E (at the final state of a simulation), it becomes rather straightforward to assess quantitatively (and in a reproducible way) the features of the shear band. First, its location and orientation are obtained by labelling all the pixels belonging to it (i.e., all the pixels with a local density value lower than the average density minus one standard deviation) and computing the center of mass and the principal axes of inertia of the obtained area (blue area in the insert of Fig. 5). A geometric zone of interest is then defined (dotted-lines square in the insert of Fig. 5) in order to focus on the central area of the shear band and to ignore membrane-related boundary effects. Longitudinal slices (dottedlines thin rectangle in the insert of Fig. 5) are then defined in order to average relevant quantities along the shear band and to extract quantitative profiles across it and through time.

Results

Reference simulation

As a reference case, the 9600 grains are introduced in a simulation with a complete Ottawa spectrum (i.e., N m = 64 ) and a maximum number of nodes on their contours (i.e., N n = 128 ). The coefficient of friction during compaction is set to c = 0 , leading to an initial solid fraction close to 0.86, and during the deviatoric loading it is set to d = 0.5 . The curve of the stress ratio (ratio of the vertical stress to the horizontal stress v ∕ c ) as a function of the axial shortening of the sample (equal to h inih ∕h ini , where h ini and h are the initial and current height of the sample) is provided in the right-hand part of Fig. 6. The term "axial shortening" is preferred to the more classical "vertical strain" because strain is heterogeneous in the sample. This curve was smoothed in time in order to ease its reading. It shows the typical behavior of a dense granular medium with a distinct stress peak (stress ratio close to 5.40) followed by a decrease of the stress ratio, until its stabilization at a value close to 2.41. For comparison, the plot also shows an experimental curve obtained with Ottawa sand in [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF]. It appears that a contact friction coefficient d = 0.5 , combined with the chosen ini- tial solid fraction and sample grain-size and morphological properties, leads to a peak stress ratio which is very similar to the experimental one. However, it appears that the axial shortening at which this peak occurs is much smaller in the numerical response (0.7%) than in the experimental one (4.5%), that the peak is much sharper in the simulation, and that the numerical sample is much stiffer in the pre-peak Fig. 6 Left: six successive stages of the reference simulation; Right: Stress ratio measured during simulation and obtained in the reference experiment in [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF], and experimental field of porosity obtained for a similar sample [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF] zone than the experimental one. Actually, the final level of axial shortening is not sufficient to reach a stress plateau in the experimental case. The explanation for this discrepancy could be related to the 2D character of the simulation, but another possible explanation lies in the number of grains involved in the sample (i.e., 9600 in the simulation and about 50,000 in the experiment). Indeed, it was already pointed out in [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF] that there is a very important effect of the sample size on the axial shortening values at which the peak and the plateau do occur. A comparison with experiments from [START_REF] Combe | Comportement du sable d'Hostun S28 au triaxial axisymetrique[END_REF], conducted with the same sand and under similar conditions but with a much large sample (height of 100.2 mm instead of 11 mm) showed that the magnitude of the peak was not modified by the sample size, but that it was much more delayed for the large sample. Likewise, the pre-peak stress response was much softer in the case of the larger sample. The result obtained with the numerical model could thus be a direct consequence of these experimental observations. Figure 6 also shows the fields of solid fraction in the sample at six different stages of the simulation. In the prepeak zone, it clearly appears that the sample is largely undisturbed, but a broad zone of increased porosity appears in the middle of the sample when the peak is reached. During softening, a cross-like localization pattern appears, with two shear bands. When the stress plateau is reached, it appears that the lower porous stone has rotated and that eventually only one of them is active (the left-lateral one). This behavior, which is consistent with the experimental observations from [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF], is made possible by the lateral (leftward) motion of the lower porous stone. Indeed, a simulation where porous stone rotation was prevented (thus restraining this lateral motion) was performed and showed that, in that case, the cross-like pattern persisted until the end of the loading (Fig. 7). This analysis thus shows that the conditions that are applied in this study at the boundaries of the sample (i.e., two flexible membranes transmitting the fluid pressure on the lateral faces and a lower porous stone with an unconstrained motion but loaded by a contact with a hemispherical cap) make it possible to observe the appearance of a single band. This result is in good agreement with experiments, as shown by the insert of Fig. 6, which provides a cross-section of the porosity field in a sheared sample of Ottawa sand as obtained by in situ µCT [START_REF] Andò | Experimental investigation of microstructural changes in deforming granular media using x-ray tomography[END_REF], and exhibits very similar patterns.

Figure 8 provides several quantitative results for the reference simulation described in Fig. 6. Figure 8a shows the stress ratio as a function of the axial shortening (raw data and smoothed curve, as presented in Fig. 6). In Fig. 8c, the "slice-averaged" solid fraction across the shear band (see Sect. 2.5) is plotted as a function of the axial shortening. This plot clearly shows the different stages of the development of a region of higher porosity. For low values of the axial shortening (i.e., before 0.7%, which corresponds to the stress peak), the dilation is homogeneous in the whole around the (not yet existing) shear band, and the solid fraction evolves from 0.86 to 0.84. Localization takes place 1% and 5% of axial shortening, i.e., during the stage of post-peak softening. During this period, the regions of the sample located at more than about 8 mm from the central axis of the shear band remain undisturbed, while the central part dilates more and more until a solid fraction close to 0.775 is reached. Afterwards, the profile of the dilation zone remains constant in time, showing that the shear band is now well-established and that the sample has reached its stress plateau.

Figure 8b shows the evolution of the average solid fraction in the shear band and in the whole sample, and confirms a volumetric stabilization at a shortening close to 5%, with stabilized solid fractions of 0.775 and 0.830 respectively. Figure 8d shows the final "slice-averaged" profile of solid fraction across the shear band, and a similar profile for the shear rate. A comparison between these two profiles is instructive, because it shows that the dilated zone is much larger than the area where shear strain actually takes place. The total width of the dilated area is indeed close to 17.5 mm (i.e. about 70D 50 ), while the total width of the sheared zone is close to 7.9 mm (i.e. about 32D 50 ). This is further confirmed by Fig. 8e, which provides similar "slice-averaged" profiles for grain rotations across the shear band. Two measures of the grain rotations are performed. The first one ("Net rotations") is based on the difference between the initial and final orientations of the grains, and disregards all the history of the grain rotations between these two times. The second one ("Absolute rotations") is based on an integration of the absolute rotational motions in time, meaning that each grain rotation that occurs at a given time step is counted as positive and cumulated with the previous rotations of this grain. As expected, the magnitude of the absolute rotation is much larger than that of the net rotation (meaning that the grains are submitted to a complex rotation history, with a lot of back-and-forth small rotations in order to accommodate the shear). We also observe that the width of significant absolute rotations (10.8 mm, i.e., 43D 50 ) is much larger than that of the net rotations (6.8 mm, i.e., 27D 50 ). It draws the picture of a shear band where actual shear and significant (i.e., net) grain rotations occur in a narrow band of about 30 grains diameters, but which necessitates important dilation and grain relative motions in a broader band of about 60 grains diameters.

Influence of a decrease in roughness

In this subsection, we report the result of a campaign of 10 simulations performed with a varying number of Fourier modes (N m = 64, 32, 16, and 8) and a varying numof contour nodes (N n = 128, 64, 32, and 16). The 10 cases are summarized in Fig. 3, and among them the case with N m = 64 and N n = 128 is the reference case which was described in the previous subsection. Each sample was compacted with an appropriate value of the coefficient of friction c (determined by trial and error) in order to reach a constant and homogeneous initial solid fraction of 0.86 in all samples. During the deviatoric loading, the coefficient of friction was set to d = 0.5 in all cases in order to isolate the influence of the surface roughness of the grains. Quantitative results are summarized in Fig. 9.

Figure 9a shows the stress ratios obtained at the peak and on the plateau for each simulation. The residual strength does not seem to be much influenced by the surface roughness of the grains, although some deviations are observed for N m = 8. In contrast, the peak strength shows very clear trends, and decreases in a regular way when the grains get smoother (from 5.40 in the reference case to 3.79 in the smoothest case). It appears however that the peak strength is much more influenced by the number of modes used in the Fourier Spectrum (i.e., by the frequency content of the surface roughness) than by the number of contour nodes (i.e., by the spatial discretization of the grains). Rather clear trends are also observed in Fig. which shows the evolution of solid fraction in the whole sample and in the shear band. In both cases, solid fraction increases steadily when the grains get smoother, especially within the shear band where it goes from 0.775 in the reference case to 0.801 when the number of Fourier modes is reduced to 8.

In Fig. 9c the absolute and net rotations of the grains in the core of the shear band are plotted. While the net rotation does not seem to be strongly affected by the smoothing of the grains, the results are more scattered for the absolute rotation. It seems, however, that this kind of motion decreases when the grains are smooth, especially if they are discretized by a low number of contour nodes (407° in the reference case, and only 218° for N m = 8 and N n = 16). Figure 9d shows results in terms of shear band width, as measured in terms of dilation (as described in Fig. and net rotation (as described in Fig. 8f). Those widths can thus be considered to represent respectively the broad area of disturbed material and the narrower core of the shear band where strain localization takes place. Both widths seem to be more influenced by the quality of the discretization, since a decrease of N n leads to a wider disturbance and shear zone, especially at low numbers of Fourier modes.

Overall, these results indicate that a smoothing of the grains leads to a lower peak strength, and to shear bands that are wider, denser, and with less absolute rotations of the grains. In the next subsection, we explore the possibility that an increase of the contact friction coefficient between the grains may counteract these effects and mimic roughness.

Combined influence of a decreased roughness and an increased contact friction

In order to study the combined influence of roughness and friction, another simulation campaign was performed with calibrated values of the friction coefficient. Namely, the friction coefficient used for compaction c was chosen in order for all the samples to have the same initial solid fraction as the reference case (as in previous subsection) and the friction coefficient used for deviatoric loading d was chosen in order for all the samples to reach the same stress ratio at peak as the reference case (i.e., about 5.40). This calibration was done by trial and error, and the resulting values of d are provided in Fig. 10. It is interesting to notice that, for the case N m = 8 and N n = 32 it was necessary to use a non-physical contact friction coefficient d = 1.08 instead of 0.5 in order to reach the same peak strength. However, in the smoothest case (i.e., for N m = 8 and N n = 16), it was not possible to reach the target peak stress ratio of 5.40, even with a contact friction coefficient as high as d = 2.

The results of this series of simulations are summarized in Fig. 11. Figure 11a shows that the peak strength is indeed constant in the 10 samples, except in the smoothest case where it saturates at 4.19 whatever the value of the contact friction coefficient. The residual strength remains rather constant in all cases. This means that, if one only focuses on the mechanical response of the sample, friction can actually replace roughness up to a certain point.

Figure 11b-d show however that, although friction was artificially increased in order to reach a target peak strength, this cannot fully compensate a decrease in roughness in terms of solid fraction, either in the shear band or in the whole sample. Hence, despite this calibration of d , the trends observed in Fig. 10 are still valid: a decrease in roughness still increases the density and the width of the shear band, and decreases the absolute grains rotations. This observation, however, is valid for very low number of modes (i.e., for Nm = 8).

Discussion and conclusion

Based on the results presented in the previous section, it that friction can only replace roughness in granular simulations under certain conditions. If one focuses only on the macroscopic mechanical response of the sample, it is possible to compensate a certain amount of smoothing of the grains by an increase in contact friction, and still to achieve peak and residual strength that are identical to the original sample. However, if the grains are too smoothed, this may require non-physical values of the contact friction. In the extreme case where only 8 Fourier modes are used (i.e., when only the general shapes of the grains are kept, but all significant roughness disappears), the fact of discretizing the grains contours too poorly (i.e., with only 16 nodes) prevents the peak strength reaching the target value. In this case, the peak strength saturates and cannot increase anymore, whatever the increase of the contact friction. This might be related to the fact that such poorly discretized grains tend to form some planar facets, which facilitates the intergranular motion without requiring a significant dilation.

From a micromechanical point of view, a decrease of the grains roughness has some consequences on the properties of the shear band. The smoother the grains, the denser and the thicker the shear band, and this trend cannot be strongly counteracted by an artificial increase of the contact friction. However, a moderate smoothing only has moderate consequences. Since such simulations are typically very costly, it is of primary importance to reduce the number of nodes on the contour of each grain, since this will reduce the number of proximities and contact detections to be performed by the contact algorithm. Meanwhile, for a given number of nodes on the contour, it is worth accounting for as much roughness as possible, since this comes at no additional cost. For this reason, we present in Fig. the localization mechanisms for four particular simulations, which correspond to respectively 128, 64, 32 and 16 contour nodes. For each of these simulations, the number of Fourier modes was the maximum that can be reasonably represented with such numbers of nodes (i.e., respectively 64, 32, 16, and 8), and the friction coefficient was calibrated to achieve the same peak strength, apart for the last one which could only achieve a peak stress ratio of 4.19.

the shear band could occur in either direction depending on the simulations, some of them were flipped horizontally in Fig. 12 in order to ease the comparison. First, it is clear that, from one simulation to another, the patterns may differ slightly. For example, the simulations with 128 and 32 nodes started to develop a cross-like pattern before one particular shear band developed, and these patterns are still visible in the final stages of the simulations. In contrast, the simulation with 64 nodes developed a single shear band from the very beginning toward a lateral sliding, and its shear band is thus clearer and less disturbed. This observation indicates that such simulations are just as prone to variability as actual experiments, and future studies should be dedicated to a deeper analysis of the repeatability of such simulations. Because of their cost, however, a reduction of the number of nodes per grain contour is necessary. Quantitative results detailed in the previous section, as well as a visual analysis of Fig. 12, indicate that a number of 32 nodes on the grain contours and of 16 Fourier modes is largely acceptable and can be used with a good level of confidence, provided that the contact friction coefficient is properly calibrated to compensate this smoothing. Indeed, such a sample can still reach the target peak strength with Fig. 12 Localization patterns for simulations with 128, 64, 32 and 16 contour nodes per using the maximum number of Fourier modes and a calibrated contact friction a reasonable calibrated value of d = 0.715 , and the attrib- utes of the shear band (thickness, solid fraction, amount of rotations, etc.) remain not too far from the reference ones (Fig. 11). This smoothing is identified as the "best compromise", as illustrated in Fig. 3. Considering fewer nodes or fewer Fourier descriptors, however, does not make it possible to reproduce both qualitatively and quantitatively the behavior of the reference sample, even for very high values of d . The expected computational cost reduction is very implementation-dependent, especially when parallelization comes into play, but in our case the simulation duration was roughly proportional to the number of nodes on the grains contour. Hence, passing from 128 nodes (reference case) to 32 nodes (best compromise) leads to a reduction of the cost by a factor 4. It could be argued that this improvement does not include the cost of the calibration of the appropriate contact friction coefficient needed to counteract the smoothing. However, it should be kept in mind that this calibration is always necessary when running DEM simulations: direct measurements of the contact friction between two grains of sand are very seldom in the literature (see for example [START_REF] Senetakis | Experimental study of sand grain behavior at their contacts with force-and displacement-controlled sliding tests[END_REF]), and can hardly be generalized to different grains of the same sample or to grains of other sands: contact friction in DEM is always back-calibrated. The roughness correction that is proposed in this paper, limited to an extent where the physics of strain localization are not strongly modified, hence comes at no additional calibration cost.

Based on this knowledge, future work will be dedicated to larger scale simulations of similar systems, and will make it possible to explore the influence of other granulometric (D50, uniformity coefficient) and morphological (First Fourier modes) parameters on the qualitative and quantitative behavior of shear localization in granular materials.

Compliance with ethical standards

Conflict of interest

The authors acknowledge that this study contains original material, as a result of a purely academic study without any kind of private funding or conflict of interest. Its publication has been approved tacitly by the responsible authorities at the institutes where the work has been carried out.

Fig. 1

 1 Fig. 1 Description of the reference experiment (taken from [50]); a Membrane fixed on the lower porous stone; b Pluviation of the sand inside the mould; c Removal of the mold; d close view of the loading system; e. Raw tomographic image of the sample before shearing

Fig. 2

 2 Fig. 2 Sample characterization; a SEM views of Ottawa sand grains; b Reconstruction of 6 grains from high resolution x-ray scans and representation as geodesic structures; c First modes of the measured

Fig. 5

 5 Fig. 5 Illustration of the postprocessing techniques used for the study of the shear bands

Fig. 7 Fig. 8

 78 Fig. 7 Cross-like localization pattern obtained in a simulation identical to the reference one but restraining the rotation of the lower porous stone

Fig. 9

 9 Fig. 9 Influence of a reduction in roughness; a Stress ratios; b Solid fractions; c Rotations; d Thickness of the shear band

Fig. 10

 10 Fig. 10 Contact friction coefficient µ d calibrated in order to reach the target peak stress ratio

Fig. 11

 11 Fig. 11 Influence of a reduction in roughness compensated by an increase in contact friction; a Stress ratios; b Solid fractions; c Rotations; d Thickness of the shear band


L Illoul

P Lorong

F Chinesta

Study of High Speed Blanking Process with the C-NEM

Keywords: C-NEM, Finite Transformations, Simulation of Forming Processes

diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The C-NEM has still a number of common points with the finite elements approach: the interpolation is defined on a set of nodes, degrees of freedom (DOF) are associated to the nodes, a localised in space shape function is associated to each node, the interpolation is nodal and is exact for linear fields. The C-NEM interpolation uses the constrained Voronoi diagram (dual of the constrained Delaunay's tesselation) associated to the set of nodes and the boundaries of the domain. The quality of the produced interpolation depends primarily on the distribution of the nodes. The presence of flat tetrahedrons in the dual mesh of Delaunay does not affect the interpolation quality. In next paragraphs we initially reconsider the NEM interpolation. We then focuse on the non-convex domains for which the construction of a constrained Voronoi diagram is needed. The aspects related to the partitioning of the domain are then evoked. Finally we give the example of a 3D simulation of blanking.

THE C-NEM INTERPOLANT

Headlines of the natural element method

The NEM (Natural Element Method) interpolant is based on the Sibson's natural neighbor coordinates (shape functions) [START_REF] Sambridge | Geophysical parameterization and interpolation of irregular data using natural neighbours[END_REF][START_REF] Sukumar | The natural elements method in solid mechanics[END_REF] and is constructed on the basis of the Voronoi diagram. For a set of nodes S = {n 1 ,n 2 ,...,n N } in ℜ dim ,(dim ∈{ 2, 3}), the Voronoi diagram is the subdivision of ℜ dim into regions T i (Voronoi cells) defined by:

T i = {x ∈ℜ dim : d(x, x i ) <d(x, x j ), ∀j = i} (1)
Sibson's coordinates of x with respect to a natural neighbor n i (see figure 1) are defined as the ratio of the overlap area (volume in 3D) of their Voronoi cells to the total area (volume in 3D) of the Voronoi cell related to point x: If the point x coincides with the node n i , i.e. (x = x i ), φ i (x i )=1 , and all other shape functions are zero, i.e. φ i (x j )=δ ij (δ ij being the Kronecker delta). The properties of positivity, interpolation, partition of unity [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF] and local coordinate property [START_REF] Sibson | A vector identity for the dirichlet tesselations[END_REF] are verified :

φ i (x)= Area(af ghe) Area(abcde) (2)
⎧ ⎪ ⎨ ⎪ ⎩ 0 ≤ φ i (x) ≤ 1 φ i (x j )=δ ij n i=1 φ i (x)=1 x = n i=1 φ i (x)x i (3)
It turns out that the support of φ i (x) is the union of the n circles (spheres in 3D) passing through the vertices of the n Delaunay triangles (tetrahedrons) connecting the node n i (in this case n is the number of natural neighbors of node n i ). Another important property of the NEM interpolant is the ability to reproduce linear functions over the boundary of convex domains. The proof can be found in Sukumar et al. [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF]. This is not true in the case of non convex boundaries, and the next section focuses on an approach to circumvent this difficulty.

The Constrained natural element method

In its original form [START_REF] Sukumar | The natural elements method in solid mechanics[END_REF], the NEM can only be applied to strictly convex domains. For strongly non-convex domains (cracks, auto-contact...) some spurious influences between nodes of the boundaries appear [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (c-nem)[END_REF]. In order to avoid these drawbacks, we have proposed in a previous paper [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (c-nem)[END_REF] an extension of the NEM in which a visibility criterion is introduced in order to restrict influent nodes among natural neighbors. The computation of the shape functions is done on the basis of the so-called Constrained (or extended) Voronoi diagram (CVD), which is the strict dual to the constrained Delaunay triangulation introduced by Seidel in [START_REF] Seidel | Constrained delaunay triangulations and voronoi diagrams with obstacles[END_REF] (see [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (c-nem)[END_REF] for further details). The new cells T C i , called constrained Voronoi cells, are defined formally by:

T C i = {x ∈ℜ n : d(x, x i ) <d(x, x j ), ∀j = i, S x→n i ∩ Γ=∅,S x→n j ∩ Γ=∅} ( 4
)
where Γ is the domain boundary, composed by a set of segments in 2D (triangular plane facets in 3D) and S a→b denotes the segment between the points a and b.

In this framework, a point located inside a cell T C i is closer to the node n i than to any other visible node n j . For both the trial and the test functions the new approximation is then :

u h (x)= i∈V φ C i (x)u i ( 5
)
where V is the set of natural neighbors which are visible from point x and φ C i is the constrained natural neighbor shape function related to the i-th node at point x. It was shown in [START_REF] Yvonnet | A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (c-nem)[END_REF] that the use of the constrained Voronoi diagram does not affect the properties of the NEM interpolation.

Constrained Delaunay tetraedrisation

The constrained Delaunay triangulation does not always exist in 3D without adding new nodes [START_REF] Schönhardt | Uber die zerlegung von dreieckspolyedern in tetraeder[END_REF]. Some techniques for constructing 3D constrained Delaunay tessellations are available by adding Steiner points [START_REF] Shewchuck | Tetrahedral mesh generation by delaunay refinement[END_REF][START_REF] Shewchuck | Sweep algorithms for constructing higherdimensional constrained delaunay triangulations[END_REF]. These points are always set on triangles belonging to the boundary of the domain. It is possible to connect each of these points to at most 3 initial nodes (those defining triangles containing the new points). The new points are treated as Slave nodes, the initial nodes being the Master nodes. For each Slave node n S i , we set the following kinematic linear constraint :

u S i = 3 j=1 η ij . u M (S i ) j (6)
where u M (S i ) j

are the nodal values at the 3 Master nodes of n S i , and η ij the reduced coordinates of the Slave node on its supporting triangle. To keep the NEM interpolation properties, the interpolation must be built using all the Master and Slave nodes values. With this partition, the interpolation becomes :

u h (x)= i∈M i φ C i (x)u M i + i∈S i φ C i (x)u S i (7)
M i and S i being respectively the set of Master and Slave natural neighbors which are visible from x.

The software TETGEN (Quality Tetrahedral Mesh Generator and Three-Dimensional Delaunay Triangulator, [START_REF] Si | Meshing piecewise linear complexes by constrained delaunay tetrahedralizations[END_REF]) is currently used to build the constrained Delaunay tetrahedrisation.

EXPLICIT LAGRANGIAN PROCEDURE

With the principle of virtual work as a basis of kinematically based C-NEM solution scheme, the corresponding continuum incremental boundary value problem is formulated as follows:

Ω t ρ(t)ü • ηdΩ t + Ω t σ t : ∇ x ηdΩ t = Ω t ρ(t)b • ηdΩ t + ∂Ω t σ τ • ηdΓ t ∀η ∈ ϑ ( 8
)
where ρ is the density, b and τ represent the body forces and applied tractions respectively, and ϑ is the space of virtual displacements. The properties dΩ t = J t dΩ 0 and ρ 0 dΩ 0 = ρ(t)dΩ t are used, which leads to:

Ω 0 ρ 0 ü • ηdΩ 0 + Ω 0 P t : ∇ X ηdΩ 0 = Ω 0 ρ 0 b • ηdΩ 0 + ∂Ω t σ τ • ηdΓ t ∀η ∈ ϑ (9)
where P denotes the first Piola-Kirschhoff stress tensor related to σ by P = JF -1 σ. The C-NEM discretization (5) of the variational form ( 9) results in the discrete set of algebraic time dependent equations, in matrix form, as:

Mü n+1 (t)=F ext n (t) -F int n (u n ,t) ( 10
)
where t is the time, M denotes the mass matrix, F int n (u,t) the internal force vector, while F ext n (t) is the external force vector. As shown in the next section, the use of the SCNI quadrature [START_REF] Chen | Non-linear version of stabilized conforming nodal integration for galerkin meshfree methods[END_REF] results in a M matrix diagonal, whose diagonal terms are given by m i = ρ 0 Ω i , with Ω i the volume of the Voronoi cell related to node n i . The velocity v = u and acceleration ü = v are approximated by using central differences with variable time steps:

v n+1/2 = v n-1/2 + ∆t 1 +∆t 2 2 ün ( 11
)
u n+1 = u n +∆t 2 v n+1/2 (12) 4 NUMERICAL INTEGRATION
In the context of finite elements, a natural partition of the domain is carried out by the elements. With C-NEM the concept of element does not exist. A natural partitioning is given by the cells of the Voronoi diagram. When cells do intersect the boundary of the domain, only their part inside the domain is kept. However, for strongly non-convex domains in 3D, some boundary cells may have negative volume inside the domain (some faces may intersect themselves inside the domain). It is thus necessary to use another partitioning. For reasons of simplicity we have chosen a partitioning based on the Delaunay tetrahedrons : each elementary volume Ω i is defined by the quarters of all tetrahedrons connected to the node n i .

We then use the stabilized conforming nodal integration (SCNI) proposed by Chen et al. [START_REF] Chen | Non-linear version of stabilized conforming nodal integration for galerkin meshfree methods[END_REF] to define a deformation gradient at node n i :

∇f h i = 1 |Ω i | Ω i ∇f h (x)dΩ (13)
In order to avoid projections between two successive actualisations of the reference configuration all variables are computed with this gradient at the nodes.

NUMERICAL SIMULATION OF BLANKING

All the simulations given below are dealing with the same example. It is a simulation of blanking where the gap e between the punch and the die is equal to 0.6mm and the radius r of the cutting edges are equal to 0.2mm (Figure 2). The sample has a thikness (in the direction of cutting) equal to 5mm, and its behavior is modeled via a Johnson-Cook law [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF] (Eq.14) with the following coefficients (SI): A =2 5 3 . 10 6 , B =6 8 5 . 10 6 , C =0 .097, m =2 .044, n =0 .3128, ǫp 0 =1., T 0 = 296. These coefficients correspond to a 306L stainless steel. ǫp the rate of plastic strain and T the current temperature. The punch and the die are considered as rigid and the contact with the sample is without friction. The two symmetries of the experiment are taken into account and only a quarter of the sample is studied. The elastic behavior take a major place in the localization of the shearing band. The elastic waves take 0.91 0 -6 s to pass across the sample and figure 4 shows that the shearing band is clearly localize after only 2.010 -6 s when elastic waves are still governing the stress distribution in the sample. On figure 4 we compare the equivalent plastic strain for 2 different punch velocities.

t =0.25μs t =0.50μs t =1.00μs t =2.00μs This figure exhibits the asymmetry of the shearing band which is more important close to the die, this effect being more important when the velocity of the punch goes from 10 to 30m.s -1 .

Punch Velocity 10ms -1 Punch Velocity 30ms -1 On both previous simulations the punch and the die are taken as rigid although their young modulus and mass density are close to those of the sample. This increases the effect of plastic strain near the contact and inside the sample. This aspect will be take into account for future simulations.

CONCLUSIONS

With the possibility of dealing with three-dimensional problems involving domains of any type (convex or not), C-NEM is from now a general tool. Its development is currently oriented towards the adaptation of the node distribution (refinement in particular) specially in zones where high gradients occur, and also to the adaptation of the boundary description when high deformations have occurred. For instance, a correct of the domain boundary is necessary to update the constrained Voronoi diagram in the context of an updated Lagrangian approach. First results on blanking are nevertheless interesting.

Figure 1 :

 1 Figure 1: Sibson shape functions construction in 2D

σ

  y =[A + B.(ǭ p ) n ] 1+C. ln ǭp represents the equivalent plastic strain,

Figure 2 :

 2 Figure 2: The blanking device

Figure 3 :

 3 Figure 3: Stress field in blanking -Punch velocity : 10ms -1

Figure 4 :

 4 Figure 4: Plastic strain localisation after 2.00µs

ACKNOWLEDGEMENT his work has been done under the technical and financial support of CETIM Senlis (Centre Technique des Industries Mécanique).