Lauren Thevin
email: lauren.thevin@inria.fr

Anke M Brock
email: anke.brock@enac.fr

Augmented Reality for People with Visual Impairments: Designing and Creating Audio-Tactile Content from Existing Objects

Tactile maps and diagrams are widely used as accessible graphical media for people with visual impairments, in particular in the context of education. They can be made interactive by augmenting them with audio feedback. It is however complicated to create audiotactile graphics that have rich and realistic tactile textures. To overcome these limitations, we propose a new augmented reality approach allowing novices to easily and quickly augment real objects with audio feedback. In our user study, six teachers created their own audio-augmentation of objects, such as a botanical atlas, within 30 minutes or less. Teachers found the tool easy to use and were confident about re-using it. The resulting augmented objects allow two modes: exploration mode provides feedback on demand about an element, while quiz mode provides questions and answers. We evaluated the resulting audio-tactile material with five visually impaired children. Participants found the resulting interactive graphics exciting to use independently of their mental imagery skills.

Introduction

In special education environments, tactile maps and graphics are widely used, for example in O&M (Orientation & Mobility) classes, or for mathematics and biology. There are two design methods for augmenting tactile media with audio feedback. First, tactile graphics can be hand-made, e.g. using paper, tissue or real objects [START_REF] Edman | Tactile graphics. American Fnd. for the Blind[END_REF] so that they possess a rich tactile dimension (e.g. texture, thermal properties). Yet, it is difficult to make those graphics interactive. Second, the tactile mediums can be digitally designed and edited, then printed and easily made interactive using Text-to-Speech synthesis (TTS) or Braille displays [START_REF] Anke | Interactivity Improves Usability of Geographic Maps for Visually Impaired People[END_REF]. But such objects are limited, for instance to swellpaper or 3D print. In specialized schools, many objects for visually impaired (VI) people are hand-made, such as small-scale models or students' works made in class. At the same time, people who want to augment tactile objects with audio feedback (e.g. teachers), do not necessarily master 2D and 3D modeling techniques. It is likely that they have knowledge about office software applications, but no specialized computer skills. scanning it and using the resulting image as basis for (c) creating the SVG file with interactive zones, using PapARt (Fig. 2, left) to project this image on a table, (d) superimposing the projected interactive zones and the real objects, and (e) using the object as an audio-tactile interactive medium.

We propose a first approach of creating augmented content for Augmented Reality (AR) with real objects, thus enabling both rich tactile feedback and interactivity. The audio-content is associated with a SVG file (vector format often used for designing raised-line maps [START_REF] Anke | Interactivity Improves Usability of Geographic Maps for Visually Impaired People[END_REF]). We present the necessary steps in Figure 1. We illustrate our work with two use cases: a botanical atlas and a geographic map. This approach allows teachers to design their own pedagogical audio-tactile content from existing objects (leaves, branches, small-scale models, etc.). For the resulting audio-tactile medium, two modes of use are proposed : exploration mode with audio feedback on demand, and quiz mode. The augmented prototype was designed and evaluated in collaboration with IRSA (special education school in Bordeaux). We tested the creation of audio-tactile content with teachers, and the usability of the resulting audio-tactile graphics with VI students.

State of the Art

Tactile and Audio Modalities in Interactive Graphics

Tactile graphics serve the same use cases as pictures for sighted people: they are small-scale representations used e.g. for education. Audio feedback allows to create interactive tactile graphics, using for instance TTS to describe content and events [START_REF] Ducasse | Accessible interactive maps for visually impaired users[END_REF]. The tactile modality encodes various types of information, and supports the creation of a rich mental imagery also for VI people [START_REF] Cattaneo | Blind vision: the neuroscience of visual impairment[END_REF]. Tactile maps, interactive or not, can be created in relief with a single material (e.g. using laser cutting, swell-paper [START_REF] Tatham | The design of tactile maps: theoretical and practical considerations[END_REF], 3D printing [START_REF] Gotzelmann | LucentMaps 3D Printed Audiovisual Tactile Maps for Blind and Visually Impaired People[END_REF]), or composed from objects with various textures [START_REF] Giraud | Map learning with a 3d printed interactive small-scale model: Improvement of space and text memorization in visually impaired students[END_REF]. Tactile graphics can support several haptic sensations: graining, softness/stiffness, shapes (circle/square, full/hollow), relief (cavity/bump), thermal behavior (cold/warm), etc. [START_REF] Klatzky | Touch. Handbook of psychology[END_REF]. Using various textures helps to distinguish tactile elements. Using realistic textures can help to recognize elements in real situations. In our use case of a botanical atlas, some plants have a particular texture (such as succulent plants) that makes them easy to recognize by touch.

Braille text normally accompanies tactile graphics [START_REF] Tatham | The design of tactile maps: theoretical and practical considerations[END_REF]. As an alternative, audio feedback allows to augment objects [START_REF] Ducasse | Accessible interactive maps for visually impaired users[END_REF]. Unlike with Braille text, there is no need for moving the hands between the object and the associated caption [START_REF] Anke | Interactivity Improves Usability of Geographic Maps for Visually Impaired People[END_REF]. The exploration and use of audio-tactile maps has proved more effective for complex content [START_REF] Giraud | Map learning with a 3d printed interactive small-scale model: Improvement of space and text memorization in visually impaired students[END_REF] than the use of Braille maps, and more efficient and satisfying for simple content [START_REF] Anke | Interactivity Improves Usability of Geographic Maps for Visually Impaired People[END_REF]. Moreover, interactivity makes it possible to provide adaptive feedback depending on the context for the same tactile support: generic captions for superficial exploration, detailed descriptions for precise exploration, and an interactive answer-question system for testing knowledge [START_REF] Albouys-Perrois | Towards a Multisensory Augmented Reality Map for Blind and Low Vision People: a Participatory Design Approach[END_REF].

Augmentation Process

We are interested in augmenting objects and making them interactive using AR. AR is not limited to visual modalities, but is defined as the combination of real objects with virtual objects in real time within a real environment independently of the modalities (vision, audition, touch,...) [START_REF] Azuma | Recent Advances in Augmented Reality[END_REF]. In our case, physical objects are combined with virtual audio feedback. We observe two main processes used to augment real objects with digital content. The first approach is to augment the object itself with electrical components [START_REF] Giraud | Map learning with a 3d printed interactive small-scale model: Improvement of space and text memorization in visually impaired students[END_REF] or to integrate components directly inside the object material [START_REF] Roquet | 3D Folded PrintGami: Transforming Passive 3D Printed Objects to Interactive by Inserted Paper Origami Circuits[END_REF]. But this approach offers little flexibility for the augmented content as it is integrated in the hardware. The second approach is to integrate augmented content in the digital representation of the object, e.g. vector images printed later as a raised-line map [START_REF] Albouys-Perrois | Towards a Multisensory Augmented Reality Map for Blind and Low Vision People: a Participatory Design Approach[END_REF]. This object is then made physical, e.g. with 3D printing. Modeling objects with the associated digital content allows to attach several annotations to the same object (e.g. content with increasing detail). The physical object can then be linked with the digital object and its augmented content. A tracking system can follow the object and the users action [START_REF] Shi | Markit and Talkit: A Low-Barrier Toolkit to Augment 3D Printed Models with Audio Annotations[END_REF]. However in this approach, the augmented object is either artificially created or digitally modeled based on an existing object.

To sum up, creating audio-tactile media is possible with real objects as a starting point or respectively by modeling and annotating digital objects. Yet, these two approaches require advanced skills in computer sciences and for building specific hardware. As far as we know, even for the tools made and tested in participatory design approaches, such as [START_REF] Albouys-Perrois | Towards a Multisensory Augmented Reality Map for Blind and Low Vision People: a Participatory Design Approach[END_REF], the audio-tactile supports were developed by the researchers and not the target users (e.g. teachers).

Proposed System

Our audio-tactile AR system should be compatible with the tactile contents usually used by teachers (swellpaper, 3D printing, scratch, small-scale models), and their design should be possible for people with no specific computer skills.

Hardware

The technical system is composed of the existing PapARt framework (software and hardware) [START_REF] Laviole | PapARt: Interactive 3D graphics and multi-touch augmented paper for artistic creation[END_REF], a PC for the SVG file creation, a scanner, and a sheet of paper with tactile elements placed on it (see Fig. 2). PapARt uses a depth camera to trigger the interactive zones usually on a table surface based on the hand position (the arm is detected, then the hand and fingers position). The contact of the finger with the tactile map is detected using depth analysis.

Creation of Interactive Content

Defining the Content. To create the tactile content, real objects are glued to an A3/A4 paper sheet which serves as support (see Fig. 1).

Defining the Reference Frame to Superimpose Annotations. One of the main issues in AR is to link an object with its augmented content. We propose to scan a A3/A4 paper with the tactile content. The scanned image becomes a background in a vector graphics software (Inkscape) to draw shapes directly on top of it, respecting de facto the scale and alignment with the real graphic. the shapes in SVG format have a textual description rendered by TTS when using the interactive prototype and additional XML fields. The users can add the fields "step" and "question" to any shape to create a quiz (see below).

Usage. When PapARt hardware and software are launched, the user chooses the correct SVG image in the application. The image is then projected, and the projection overlaid on the physical support. Thus, the map or diagram becomes audio-tactile and usable in two modes: 1) in exploration mode an audiodescription is provided once the objects are pointed at as in [START_REF] Albouys-Perrois | Towards a Multisensory Augmented Reality Map for Blind and Low Vision People: a Participatory Design Approach[END_REF], and 2) in quiz mode the system reads questions which the user solves by pointing to the requested element. If the pointed object is not correct, the computer provides directional help to redirect the finger as described in [START_REF] Albouys-Perrois | Towards a Multisensory Augmented Reality Map for Blind and Low Vision People: a Participatory Design Approach[END_REF]. It is possible to load another SVG file for the same tactile content so the audio-content can evolve in the same session (e.g. quiz getting progressively more difficult).

Evaluation

Six O&M Instructors and teachers were asked to evaluate the creation process, and five VI students were asked to evaluate the resulting audio-tactile media. We User Test 1: Can people without Specific Computer Skills Annotate Existing Tactile Material with Audio Feedback? Three O&M instructors from IRSA were asked to augment an existing tactile map of the school. Three different teachers (2 from IRSA: biology, tutor; 1 teacher in a primary school) were asked to create and augment a tactile biological atlas. We assessed the workload using Nasa-TLX [START_REF] Hart | Development of NASA-TLX (Task Load Index): Results of Empirical and Theoretical Research[END_REF] and the usability using UEQ [START_REF] Laugwitz | Construction and Evaluation of a User Experience Questionnaire[END_REF]. To be sure that the system is compatible with the available resources and the usual activity of the instructors and teachers, we first asked which material is available at their workplace. We asked the users to spontaneously design a tactile support and to describe the audio feedback which the support should have. This allowed to verify that our audio augmentation process is consistent with the audio interaction which teachers expect. To be sure the users were able to create the audio augmentation autonomously, we provided a manual and asked the participants to follow the instructions.

User Test 2: Are the Created Audio-Tactile Supports Acceptable and Usable by Students with VI? We tested whether the created maps are usable by the end users (i.e. VI children). Five VI children (11 to 13 years old, 2F and 3M, 2 low vision and 3 blind) were asked to use the audio-tactile supports created in User Test 1 by one of the previously participating teachers (map and botanical atlas). The children's teacher was present during the experimentation. We assessed the usability with the UEQ [START_REF] Laugwitz | Construction and Evaluation of a User Experience Questionnaire[END_REF].

Proof of Concept: Are Existing Tactile Supports Easy to Annotate? We validated that the proposed design process provides the possibility to augment various existing tactile supports which are not specifically designed to be audio-augmented. Three existing tactile supports were augmented by the research team with audio feedback as requested by instructors of IRSA.

Results

User Test 1. All users succeeded in creating the audio-tactile content. The content for exploration mode was intuitive and easy for all, while the creation of quiz content was perceived as complicated by some users even if they succeeded. The average time for describing the tactile support was 11 min (min: 2, max: 23), 13 min (min: 8, max: 25) for creating the tactile support, 4 min (min: 1, max: 10) for describing the audio-augmentation and 27 min (min: 15, max: 51) for creating the audio augmentation. All participants needed one more iteration to adjust the audio feedback. The main errors in creating audio feedback were: confusion between questions that can be answered by pointing at an element (e.g. where is the library?) and questions that require verbal replies which is not foreseen in the system (e.g. which building is here?); and adjustment if the answering zone was first too small. Because a scanner was not always available, some of the augmentations were done using a photo, taken as vertical as possible. The results of Nasa-TLX and UEQ (Fig. 4) show that the system is usable without knowledge in programming and vector drawing. A participant noticed that "It is complicated the first time, because I never used it [Bezier curves, inkscape]. But now I feel I can do it again, and it will be ok. I need to manipulate only one software." One user was missing labels indicating which parts of the leaves or map elements were interactive. Such "interactive points" can be made with our system without electronic devices on the tactile maps. Some negative remarks concerned Inkscape (not usable with a picture that has a large filesize as it become slow, no feedback regarding actions).

User Test 2. In this user test we tested the audio-tactile document designed by the professors with VI students. The main result, noticed also by the teachers present during the experiment, was that the technology is interesting for chil-dren with good mental representation skills who are already attracted by tactile images (3 of 5 participants), as well as for the children for whom mental imagery tasks are more challenging and who do not necessarily appreciate tactile images (2 of 5 participants). Indeed, 4 children want to participate in future sessions with our prototype; and 3 wanted to continue exploring new tactile supports after the experiment ended. Even though the system is usable by both low vision and blind people, the feedback about the detection of finger positions is entirely visual and thus not accessible to blind people. Moreover, concepts as "no finger should hide other fingers" (to avoid occlusion for the image recognition) are hard to explain. Yet, such feedback can be useful as sometimes the finger detection for the younger participants was not very precise. We think this could be fixed with a better calibration of the system.

Proof of Concept The instructors from IRSA provided 3 existing tactile supports which we augmented. Students had built a tram station with magnets on a magnet board (Fig. 3). The time to augment was 6 min, and it was fully operational after 1 additional iteration (adapting the text read by the TTS). One usage was identified: by projecting a "photo" of the original magnets' positions, it is possible to re-create the work of the students later on. The Velcro board had similar results. The small-scale model of a house had cardboard as vertical walls and textures on the floor. The time for augmentation of this model was 15 min, the audio-content was fully operational after 1 more iteration (adapting the interactive zones with the finger detection due to the 3D structure). We also built a quiz mode where a question could have multiple solutions.

Discussion & Perspectives

In this paper we presented the use of the PapARt AR toolkit [START_REF] Laviole | PapARt: Interactive 3D graphics and multi-touch augmented paper for artistic creation[END_REF] and Inkscape to easily and quickly make tactile material interactive. We tested the creation process with teachers. Generally, teachers found the tool easy to use and were confident about re-using it. Then, we evaluated the created material with VI children. They found it exiting and were willing to use it again. However, the finger detection needs to be improved. The annotation of real objects presented in this paper opens new perspectives, such as creating tactile and tangible prototypes in creativity sessions based on manual activities with physical materials (lego, playdough etc.), which can become interactive within 30min. Some negative feedback concerned the usability of inkscape. Therefore we propose that in the future, the interactive content could be directly created using PapARt only.

Fig. 1 .

 1 Fig.1. Augmentation process: (a) designing the tactile content with real objects, (b) scanning it and using the resulting image as basis for (c) creating the SVG file with interactive zones, using PapARt (Fig.2, left) to project this image on a table, (d) superimposing the projected interactive zones and the real objects, and (e) using the object as an audio-tactile interactive medium.

Fig. 2 .

 2 Fig. 2. Left: PapARt Hardware. Middle: SVG file with drawing over the scanned image. An audio-description is associated with the drawn elements. Right: Overlaying process. The projection displays the interactive zones. Once the tactile medium is aligned with the projection, the audio-tactile content is ready to use (here projection and objects are not aligned for clarity of image).

Fig. 3 .

 3 Fig. 3. Left: User Test 1, augmentation of a map and botanical atlas. Middle: User Test 2, two botanical atlases visible: one with felt and baize, one with real leaves on the background. Right: proof of concept: augmentation of a magnetic board

Fig. 4 .

 4 Fig. 4. TOP. Left Results of UEQ for User Test 1: all criteria are evaluated as good or very good. Right: Results of Nasa-TLX for User Test 1: none of the average indicators is associated with important workload. BOTTOM. Left Results of UEQ for User Test 2: lower scores for efficiency and dependability can be explained by finger detection. Right A visually impaired child exploring the map.

Acknowledgment. We thank the participants (teachers and students) especially C. Briant; and J. Laviole (Reality Tech) and J. Albouys for participating in the design of this tool. Done with the support of the EU Erasmus+ Pr. no 2016-1-EL01-KA201-023731.


Emilien Azéma

Farhang Radjai

Robert Peyroux

Vincent Richefeu

Gilles Saussine

Short-time dynamics of a packing of polyhedral grains under horizontal vibrations

Keywords: 83, 80, Fg Granular solids -45, 70, Cc Static sandpiles; granular compaction

We analyze the dynamics of a 3D granular packing composed of particles of irregular polyhedral shape confined inside a rectangular box with a retaining wall subjected to horizontal harmonic forcing. The simulations are performed by means of the contact dynamics method for a broad set of loading parameters. We explore the vibrational dynamics of the packing, the evolution of solid fraction and the scaling of dynamics with the loading parameters. We show that the motion of the retaining wall is strongly anharmonic as a result of jamming and grain rearrangements. It is found that the mean particle displacement scales with inverse square of frequency, the inverse of the force amplitude and the square of gravity. The shorttime compaction rate grows in proportion to frequency up to a characteristic frequency, corresponding to collective particle rearrangements between equilibrium states, and then it declines in inverse proportion to frequency.

Introduction

The dynamics of dense granular materials subjected to vibrations involves collective phenomena resulting from kinematic constraints (steric exclusions, weak spatial dimensions, . . . ) and energy dissipation [1]. Well-known examples of the vibration-induced phenomena are compaction, convective flow, size segregation and standing wave patterns at the free surface [2,3,4,5,6,7,8]. Three different states can be distinguished depending on the intensity and frequency of vibrations: 1) Gas-like or fluidized state: The rate of energy input is such that there are no enduring contacts between particles and the material behaves as a dissipative gas [9,[START_REF] Brennen | Vertical oscillation of a bed of granular material in Powders and Grains 93[END_REF]17]. 2) Solid-like state: Vibrational energy propagates through the network of enduring contacts between particles and the material undergoes slow rearrangements and progressive compaction [START_REF] Weathers | Effects of Horizontal Vibration on Hopper Flows of Granular Material in Symposium of Mechanics of deformation and flow of particulate materials[END_REF]21,10,20].

3) Liquid-like state: Both particle migration and enduring contact networks are involved in the dynamics and various collective effects can be observed [3,18,[START_REF] Aoki | [END_REF]5,19].

Vibro-compaction is the main feature of the solid-like state. Most investigated systems are unconfined granular beds (with a free surface) subjected to vertical vibrations. The vibrations behave as a source of randomness allowing the system to explore metastable configurations and to reduce its potential energy. There is, however, another mechanism which prevails in the case of confined granular materials. It is well-known that under cyclic straining, Send offprint requests to: azema@lmgc.univ-montp2.fr a granular material accumulates plastic deformation and the solid fraction tends to a maximum value depending on the material. This phenomenon is sometimes called "granular ratcheting" due to the irreversible character of compaction under cyclic loading [14,[START_REF] Luding | A minimal model for slow dynamics: Compaction of granular media under vibration or shear in Compaction of Soils, Granulates and Powders[END_REF][START_REF] Vandewalle | [END_REF]. In both confined and unconfined geometries, the solid fraction evolves as a logarithmic function of the number of cycles. In most work reported on vibrated granular media, the collective dynamics of the particles and the influence of various parameters related to the material or the driving system have not been investigated in all details. Moreover, in nearly all studies, spherical or nearly spherical particles in 3D or disks or polygons in 2D have been used [14,[START_REF] Azéma | [END_REF][START_REF] Vandewalle | [END_REF].

In this paper, we present a numerical investigation of the dynamics and short-time compaction of a system of irregular polyhedral particles confined inside a rectangular box with a retaining wall subjected to horizontal harmonic loading. This system is different from nearly all experimental systems investigated under horizontal vibrations since the packing remains confined inside a box, so that the gravity plays little role during the inward motion of the retaining wall. This geometry is similar to that used in various industrial applications such as the casting of fresh concrete where efficient vibro-compaction of dry and wet granular materials represents a crucial issue [22,23]. The tamping operation on railway ballast is another interesting case where vibrating tamping bars are used to restore the initial geometry of the track distorted as a result of ballast settlement [START_REF] Saussine | Contribution à la modélisation de granulats tridimensionnels : application au ballast[END_REF][START_REF] Oviedo-Marlot | Etude du comportement du ballast par un modèle micromécanique[END_REF]. With the increase of commercial speed, a better understanding of the physics of compaction is important for long time stability of ballast.

We are interested here in the evolution of the packing in the course of harmonic loading, the short-time compaction (during the first cycles) and the scaling of the dynamics with loading parameters. We used discrete-element numerical simulations by means of the contact dynamics method in 3D with rigid irregular polyhedral particles [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics : application to granular materials[END_REF][START_REF] Radjai | Multicontact dynamics of granular systems[END_REF][START_REF] Dubois | LMGC90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF]29]. The system is explored for a broad set of loading parameters including the frequency and amplitude of the harmonic driving force. A similar study was recently performed in 2D with polygonal particles [START_REF] Azéma | [END_REF]. The passage from 2D to 3D from a purely numerical point of view involves numerical handling of particles of polyhedral shape and a higher numerical efficiency making 3D simulations over many cycles and for varying parameters possible. In this paper, we will revisit the same phenomenology as in 2D for irregular polyhedral particles. This shows that in granular materials the dynamics is not sensitive to space dimensionality although the influence of particle shape and the details of structural organization can only be appreciated in a 3D geometry.

We begin with an introduction to 3D contact dynamics method as applied to polyhedral particle shapes and the numerical procedures. Then, we present in three sections the dynamics of the packing, the evolution of solid fraction and scaling with loading parameters.

Numerical method

The simulations were carried out by means of the contact dynamics (CD) method with irregular polyhedral particles. In this section we present the properties of this numerical method and compare it to a more classical numerical approach molecular dynamics (MD). Then we presents the principle of the contact detection and numerical parameters which have been used for this study.

Contact dynamics

The CD method is based on implicit time integration and nonsmooth formulation of mutual exclusion and dry friction between particles in case of contact [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics : application to granular materials[END_REF][START_REF] Radjai | Multicontact dynamics of granular systems[END_REF][START_REF] Dubois | LMGC90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF]. The equations of motion for each particle are formulated as differential inclusions in which velocity jumps replace accelerations [START_REF] Moreau | Some numerical methods in multibody dynamics : application to granular materials[END_REF]. The unilateral contact interactions and Coulomb friction law are represented as set-valued force laws according to convex analysis. The implementation of the time-stepping scheme requires the geometrical description of each potential contact in terms of contact position and its unit normal vector.

At a given step of evolution, all kinematic constraints implied by enduring contacts and the possible rolling of some particles over others are simultaneously taken into account, together with the equations of dynamics, in order to determine all velocities and contact forces in the system. This problem is solved by an iterative process (nonlinear Gauss-Seidel method) which consists of solving a single contact problem, with other contact forces being treated as known, and iteratively updating the forces until a given convergence criterion is fulfilled. The method is thus able to deal properly with the nonlocal character of the momemtum transfers -resulting from the perfect rigidity of particles in contact. The CD method makes no difference between smooth evolution of a system of rigid particles during one time step and nonsmooth evolutions in time due to collisions or dry friction effects.

The MD-like methods are based on regularization schemes where impenetrability is approximated by a steep repulsive potential and Coulomb's law by a viscous-or elasticregularized friction law, to which smooth computation methods can be applied. In this case, the choice of viscous or elastic parameters has to be specified depending on the particle shapes. This regularization implies small time steps in order to ensure numerical stability whereas the implicit time integration method, inherent in the CD method, is unconditionally stable. The uniqueness is not guaranteed by CD approach for perfectly rigid particles in absolute terms. However, by initializing each step of calculation with the forces calculated in the preceding step, the set of admissible solutions shrinks to fluctuations which are basically below the numerical resolution. In MD-based simulations, this "force history" is encoded by construction in the particle positions.

For our simulations, we used the LMGC90 which is a multipurpose software developed in Montpellier, capable of modeling a collection of deformable or undeformable particles of various shapes (spherical, polyhedral, or polygonal) by different algorithms [START_REF] Dubois | LMGC90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF].

Simulation of polyhedral particles

The treatment of a contact interaction between two particles requires the identification of the contact zone and a plane (in 3D), the so-called "common plane". Obviously, for rigid particles it is possible to define this contact zone by a finite set of points. Before applying the contact detection algorithm between a pair of particles of irregular shapes (see Fig. 2), a "bounding box" method is used to compute a list of particle pairs potentially in contact. Then, for each pair, the first step is to determine if an overlap exists through a 3D extension of the "shadow overlap method" [START_REF] Saussine | Contribution à la modélisation de granulats tridimensionnels : application au ballast[END_REF][START_REF] Dubois | LMGC90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF]. Several algorithms exist for overlap determination between convex polyhedra [START_REF] Nezami | [END_REF][START_REF] Dubois | LMGC90 une plateforme de développement dédiée à la modélisation des problèmes d'intéraction[END_REF]36]. In the case of an overlap, the contact plane is determined by computing the intersection between the two particles.

The contacts between polyhedral particles belong to different categories, namely face-face, edge-face, vertexface, edge-edge, vertex-vertex, vertex-edge; see Fig. 1. The vertex-vertex and vertex-edge contacts are rare, i.e. in practice of "zero measure". For each case we determine one, two or three contact points which provide a good description of the contact zone. In this paper, the vertexedge and edge-edge contacts are referred to as "simple" contacts. On the other hand, the edge-face and face-face contacts are treated as "double" and "triple" contacts, respectively, as their representation involves 2 and 3 distinct points on the common plane. The detection procedure is fairly rapid and allows us to simulate large samples composed of polyhedral particles.

Numerical samples

Our numerical samples are composed of rigid polyhedral particles with shapes and sizes that represent those of ballast grains (Fig. 2). Each particle has at most 70 faces and 37 vertices and at least 12 faces and 8 vertices. A sample contains nearly 1200 particles. The particle size is characterized as the largest distance between the barycenter and the vertices of the particle, to which we will refer as "diameter" below. We used the following size distribution: 50% of diameter d min = 2.5 cm, 34% of diameter 3.75 cm, 16% of diameter d max = 5 cm. The particles are initially placed in a rectangular box and compressed by downward motion of the upper wall at zero gravity; see Fig. 3. Then, the gravity is set to g and the upper wall is raised 1 cm and fixed. The coefficient of friction between the particles and with the horizontal walls was fixed to 0.4, but it was 0 at the vertical walls. The coefficient of restitution between particles was fixed to zero because of the high solid fraction of the samples. One of the walls is allowed to move horizontally (x direction in Fig. 3) and subjected to a harmonic driving force. All other walls are immobile. For all simulations the time step was 2.10 -4 s.

Vibrational dynamics

The free wall is subjected to a harmonic force as a function of time,

f (t) = (f max + f min ) 2 - (f max -f min ) 2 sin ωt, (1)
where f max and f min are the largest and lowest compressive (positive) forces acting on the wall. The first term represents the mean confining force modulated by the second term. At t = 0, the external force f = (f max + f min )/2 causes the inward motion (contraction) of the free retaining wall. Jamming occurs when the gap left between the upper wall and the free surface of the packing is filled. If f min is above the (gravitational) force exerted by the particles on the free wall, f will be large enough to prevent the wall from backward motion (extension) during the whole cycle. Then, the granular material is in "passive state" and the major principal stress direction is horizontal [START_REF] Nedderman | Statics and Kinematics of Granular Materials[END_REF].

On the other hand, if f max is below the force exerted by the particles, f will never be large enough to prevent the extension of the packing. This corresponds to the "active state" where the major principal stress direction remains vertical. In all other cases, both contraction and extension occur during each period, and the displacement ∆x of the free wall will be controlled by f min and f max . Without loss of generality, we set f min = 0. This ensures the largest possible displacement of the wall in the active state. Four different values of f max were tested, ranging from 2.10 3 N to 10 4 N.

We first consider the trajectory x(t) of the free wall which reflects the dynamics of the particles in the box in response to harmonic forcing. Figure 4 shows x(t) for frequency ν = 5 Hz over a time interval ∆t = 1 s. We can observe a fast initial contraction (t < 0.1 s) followed by slow contraction over four periods. The initial contraction is a consequence of the gap left between the free surface of the packing and the upper wall. The subsequent periodic motion of the wall takes place around this jammed state and will be studied below. A zoom on a single period is shown in Fig. 5. It begins at the jamming position x = x j corresponding to the jamming position reached at the end of the preceding period. The motion of the wall begins (point a in Fig. 5) only when the applied force f declines near to its minimum f min = 0. The maximum displacement ∆x max occurs at a later time δt r (point b). From a to b, the force exerted by the packing on the free wall is above the applied force, so that the wall moves backward (extension). In this phase, the packing is in an active state. The inverse situation prevails from b to c where the particles are pushed towards the box (contraction). Then, the packing is in a passive state. The new jamming position x ′ j is below the jamming position x j reached at the end of the preceding period. The difference x j -x ′ j represents the net compaction of the packing over one period. The particle velocity field is not a simple oscillation around an average position. The particles undergo a clockwise convective motion in the cell as shown in Fig. 6. Figure 7 shows the horizontal force f g exerted by the packing on the wall as a function of x over four periods. In the active phase, f g grows slightly with x. In the passive phase, it grows faster and almost linearly as x decreases. The vertical line corresponds to the jammed state where f g decreases with f at x = x j . We also observe two transients : 1) unjamming and the onset of the active state, 2) jamming from the passive state. Inside the packing, the contact forces evolve between a fully jammed state, where horizontal force chains dominate (Fig. 8(a)), and the active state, where vertical gravity-induced chains can be observed (Fig. 8(b)).

Granular ratcheting

In our system, the solid fraction ρ of the packing increases due to horizontal vibrations. This accumulation of plastic strain under oscillatory loading is sometimes called "granular ratcheting" [14,[START_REF] Luding | A minimal model for slow dynamics: Compaction of granular media under vibration or shear in Compaction of Soils, Granulates and Powders[END_REF]. To evaluate ρ, we consider a control volume inside the box. The initial value of the solid fraction is 0.50. Figure 9 shows the evolution of the variation ∆ρ of solid fraction for several periods. An initial compaction of 2% is followed by oscillations with a small increase of ∆ρ in each period. The initial compaction should be attributed to the initial state of the packing is not yet fully confined. We use ρ 0 = 0.51, reached after a time laps of 0.2 s, as the reference value for the evolution of solid fraction. The relative compaction of the packing is given by ∆ρ/ρ 0 . The compaction rate η over several periods and for a total time interval ∆t is

η ≡ 1 ρ 0 ∆ρ ∆t . ( 2
)
The compaction of the packing slows down logarithmically at long times [18]. But, the short-time compaction can well be approximated by a linear function with a constant compaction per period ∆ρ 1 . Then, we have

η = ∆ρ 1 ρ 0 ν. (3)
For ν = 5 Hz and f max = 6 10 3 N, we have η ≃ 0.025 s -1 . This rate is faster in 3D compared to 2D simulations for the same frequency [START_REF] Azéma | [END_REF].

It is important to note that compaction occurs in the active state, i.e. during the extension of the packing. This is shown in Fig. 10, where the variation ∆ρ of the solid fraction is plotted as a function of x. The solid fraction increases during extension (increasing x) and decreases during contraction (decreasing x).

Granular ratcheting has been investigated by numerical simulations showing that the anisotropy of critical contacts, where the friction force is fully mobilized, plays an important role [14,[START_REF] Luding | A minimal model for slow dynamics: Compaction of granular media under vibration or shear in Compaction of Soils, Granulates and Powders[END_REF]. Quasi-static cyclic shearing also leads to cumulative compaction of a granular material at low strain amplitudes [START_REF] Radjai | Roux Contact Dynamics Study of 2D Granular Media : Critical States and Relevant Internal Variables in The Physics of Granular Media[END_REF][START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. At large amplitudes, the compaction is followed by decompaction (dilation) and no net compaction can be observed over a full cycle. In our system, compaction is a consequence of unjamming and it is pursued during the whole active state. Decompaction takes place in the passive state, but it is cut short by fast jamming. The outcome of a full cycle is thus a net compaction of the packing.

Influence of loading parameters

We performed a series of simulations for frequencies ν ranging from 1 Hz to 60 Hz and for a total time of 1 s. All simulations yield similar results both for dynamics and compaction. Moreover, a simple dimensional analysis leads to the collapse of the data on a single plot. Indeed, the frequency sets the time scale τ = ν -1 . Force scales are set by the largest driving force f max in the passive state and the particle weights mg as well as the smallest driving force f min in the active state. Hence, dimensionally, for fixed values of mg, f min and f max , all displacements are expected to scale with ν -2 and all velocities with ν -1 . This scaling is shown in Fig. 11 where the phase space trajectory is shown for ν = 5 Hz and ν = 10 Hz without scaling and after scaling the displacements ∆x by ν -2 and the velocities v by ν -1 . We see that the data from both simulations collapse nicely on the same trajectory after scaling. Figure 12 shows the maximum displacement ∆x max in the active state and the maximum velocity v max in the passive state as a function of ν. The corresponding fits by ν -2 and ν -1 are excellent.

The role of force parameters mg, f min and f max is less evident. Since we have f min = 0, we expect ∆x max to be dependent on the ratio mg/f max representing the relative importance of the gravitational to driving forces. Indeed, our data show that ∆x max varies as f -1 max ; Fig. 13. On the other hand, the mass ratio m w /m, where m w and m are the mass of the free wall and the total mass of the particles, must control the inertia and thus the maximum displacement of the wall. Our simulations with different values of m w show that ∆x max varies as m/(m + m w ).

Hence, we propose the following expression for the scaling of displacements with loading parameters:

∆x max = C m m + m w mg f max g ν 2 , ( 4
)
where C is a dimensionless prefactor. Fig. 15 shows ∆x max as a function of (mg) simulations with different values of ν, f max , g and m w . We see that the data are in excellent agreement with Eq. 4. The prefactor is C ≃ 0.01. This scaling is the same as in 2D simulations with a material constant C ≃ 0.05 for polygonal particles [START_REF] Azéma | [END_REF]. Let us also remark that Eq. 4 predicts that ∆x max varies as g 2 . This prediction agrees well with our simulation data shown in Fig. 14 for four different values of g.

Compaction rates

We now come back to granular ratcheting and we would like to evaluate short-time compaction rates as a function of frequency. According to Eq. 3, the compaction rate η varies linearly with the frequency ν if the total compaction per period ∆ρ 1 is independent of ν. Fig. 16 shows η as a function of ν. We see that only at low frequencies, η increases linearly with ν. represents the characteristic time for the relaxation of the packing. In the active state, the packing needs a finite rearrangement time to achieve a higher level of solid fraction. As long as the vibration period τ = ν -1 is longer than τ c , the packing has enough time to relax fully to a more compact equilibrium state. But, when the period τ is below τ c , the relaxation is incomplete so that ∆ρ 1 < ∆ρ max , where ∆ρ max is the largest compaction between two periods.

It is expected that ∆ρ 1 should follow the same scaling with the frequency as the displacement of the retaining wall, i.e. ∆ρ 1 ∝ ∆ρ max ν -2 . This is because the volume change ∆V is proportional to ∆x. Hence, from Eq. 3 and imposing the continuity at ν = ν c , we get

η = ∆ρmax ρ0 ν ν < ν c , ∆ρmax ρ0 ν 2 c ν -1 ν > ν c . (5)
Fig. 16 shows this prediction together with the data points. We see that, although ν c is the only fitting parameter, the compaction rate η is well adjusted by Eq. 5. The prefactor ∆ρ max /ρ 0 is ≃ 0.005, corresponding to ∆ρ max ≃ 0.0025. In spite of the sharp transition at ν = ν c , it is convenient to construct a single expression containing the correct behavior both at low and high frequencies. As in 2D for polygon packings [START_REF] Azéma | [END_REF], the following fitting form provides a good approximation as shown also in Fig. 16 (fitting form 2):

η = ∆ρ max ρ 0 1 + e -( ν νc -1) 2 1 + ν νc 2 ν. (6)
The characteristic time τ c = 0.1 s is of the same order of magnitude as the time required for one particle to fall down a distance equal to its diameter. Obviously, the above findings concern only short-time compaction (∆t < 1 s). At longer times, η declines with time, but the scaling with frequency according to Eq. 5 is expected to hold at each instant of evolution of the packing.

Conclusion

In this paper, the contact dynamics method was employed to simulate and analyze the dynamics of a system of polyhedral particles subjected to horizontal harmonic forcing of a retaining wall. Our system is devoid of elastic elements and, hence, the behavior is fully governed by particle rearrangements. Moreover, it involves a jammed state separating passive (loading) and active (unloading) states. Dimensional analysis was used to scale the displacements with the frequency of oscillations. It was shown that the data collapse by scaling the displacements by the inverse square of frequency. We also studied the scaling with confining force and particle weights.

Granular ratcheting under horizontal vibrations was investigated. During each vibration period a small compaction of the system occurs during unloading, i.e. upon sample extension, followed by decompaction upon contraction. The compaction rate increases linearly with frequency up to a characteristic frequency and then it declines in inverse proportion to frequency. The characteristic frequency was interpreted in terms of relaxation time of the packing under its own weight during the unloading phase.

The similarity of the phenomenology of vibrational dynamics and compaction at short times in the 3D system of polyhedral particles with that of a 2D system of polygonal particles suggests that space dimensionality plays a minor role in granular dynamics. The characteristic times and compaction rates are slightly different but the scaling behavior and the functional dependence of the compaction rate with frequency are the same. A comparison with spherical particles would be interesting in order to highlight the effect of particle shapes on these parameters. On the other hand, the characteristic time appears as a crucial parameter for the compaction rate and it merits further investigation as a function of various control parameters of the system.

The authors would like to thank specially F. Dubois for interesting discussions and help with the software LMGC90. This work was supported by a grant from the Région Languedoc-Roussillon and the french railway company SNCF.

Fig. 1 .

 1 Fig. 1. Different types of contacts between two polyhedra.

Fig. 2 .

 2 Fig. 2. Examples of polyhedral shapes used in the simulations.

Fig. 3 .

 3 Fig. 3. A snapshot of the packing inside a box with a free wall over which the driving force f (t) is applied along the x direction.

Fig. 4 .Fig. 5 .

 45 Fig. 4. The evolution of the displacement x of the free wall (up) in response to harmonic loading (down).

Fig. 6 .Fig. 7 .

 67 Fig. 6. Instantaneous particle velocity field in the passive state, i.e. during inward motion of the free wall.

Fig. 8 .

 8 Fig. 8. Normal forces in the passive (a) and active (b) states in a section of the packing parallel to the xz plane. The segments connect particle centers with a thickness proportional to the normal force.

Fig. 9 .

 9 Fig. 9. Evolution of the solid fraction ∆ρ from the initial state as a function of time over several periods.

Fig. 10 .

 10 Fig. 10. Variation ∆ρ of the solid fraction from the initial state as a function of the displacement x of the free wall.

Fig. 11 .

 11 Fig. 11. Phase space trajectories for two frequencies without scaling (a) and with scaling (b) of the displacements and velocities with respect to the frequency.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. Maximum displacement ∆xmax (a) and the maximum velocity vmax (b) as a function of frequency ν.

Fig. 14 .Fig. 15 .

 1415 Fig. 14. Scaling of the maximum displacement ∆xmax with gravity g.

2 .Fig. 16 .

 216 Fig. 16. The compaction rate η as a function of the frequency (circles) fitted by two different functions; see text.


L Diogo

Philippe R B Cecílio

Sonia M Devloo

Erick S R Gomes

Santos

Diogo L Cecílio

Philippe R B Devloo

Sônia M Gomes

Erick S R Santos

An adaptive finite element

Keywords: Geomechanics, Plasticity, Finite elements, Mesh Adaptation, Return-Mapping Algorithm, Welbore stability

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

During the drilling operation of a well, the in-situ stress state, acquired during millions of years, changes abruptly. Well-bore stability analysis refers to the control of stress changes during the well excavation process, which includes the movement of rock material by the drill bit, and its replacement by the drilling fluid. The confining pressure becomes balanced by the hydrostatic pressure generated by the column of drilling fluids, as illustrated in Figure 1.

The pressure of the drilling fluid is a very important parameter. If the hydrostatic pressure exerted by the column of drilling fluid is excessively high, the fluid may cause the the rock to fracture, causing significant loss of fluid and damage to the reservoir formation (hydraulic fracturing). If the pressure is insufficient, there may be the invasion of the fluid contained in the reservoir. This invasion can cause fluid separation, resulting in expansion of gas contained in the oil (e.g. Deep Horizon disaster). Even when the pressure of the drilling fluid is between the two thresholds mentioned, the stress state change can lead to plastic deformation of the material of the formation and to the eventual collapse of the well. This altered state of stress can be partially controlled by the adjusting the pressure of the drilling fluid by changing its density. Increase of reliability of wells can be obtained by a more detailed study of the stress state around the well, including the non-linear material behavior due to plastic deformation and by including numerical models of breakout of material. In this context, the main purpose of the present article is to study the stress state around the well-bore using computational elastoplasticity [START_REF] Chen | Plasticity for Structural Engineers[END_REF][START_REF] Neto | Computational Methods for Plasticity[END_REF][START_REF] Borja | Plasticity: Modeling & Computation[END_REF][START_REF] Lubliner | Plasticity theory[END_REF]. A finite element formulation of the stress state of excavated wells is developed. Additionally, the second invariant of the plastic deformation tensor is adopted as a damage criterion to simulate material collapse (breakout), which is represented numerically by adjusting the geometry of the domain.

The following advanced numerical finite element algorithms have been combined to construct the proposed numerical model of breakout:

• An improved numerical integration algorithm proposed in [START_REF] Cecílio | An improved numerical integration algorithm for elastoplastic constitutive equations[END_REF], which holds for associative models.

• hp-Refinement around the well-bore. In order to enhance the precision of the elastoplastic simulation, curved meshes directional refined towards the well-bore are used, and the polynomial order is increased in the plastification area.

• A specific transfer procedure for elastoplastic deformation history. Considering the fact that in elastoplastic simulations the stress state is path dependent, adaptive simulations require the communication of the deformation history from one mesh to another.

• An automatic numerical method adapting the geometry of the computational domain by removing the region where the material loss occurs, which is associated to an given value of plastic shear.

Concerning the constitutive plastic model, the focus on this work has been on cases that could test the versatility of the plastic integration algorithm and still resemble the major mechanical response of a rock under high burial stresses. In order to capture the true mechanical behavior of the rock, the shear and traction failure modes must be reproduced. Rocks are materials with intrinsic nonlinear mechanical response. Although all the nonlinearities, such as those in the elastic behavior, cap envelop with hardening/softening and non-associativeness of plasticity in the shear envelope could possibly be modeled, it is challenging to devise numerical properties to fulfill the most sophisticated constitutive models. Material availability and test facility capabilities often reduce the amount of real material information at hand when determining material parameters. The real rock behaviour also includes pore collapse and excessive dilatation.

The Dimaggio-Sandler model [START_REF] Dimaggio | Material Model for Granular Soils[END_REF][START_REF] Sandler | An algorithm and a modular subroutine for the cap model[END_REF] is an accepted model to represent the behavior of soil materials, and can reproduce all the phenomena described above. Originally, this model has been conceived for a specific sand material, but the perfectly-plastic shear and strain hardening compression cap potentials have shown good overall match of the laboratory observed rock behavior [START_REF] Santos | Simulador de meios porosos saturados elastoplásticos[END_REF][START_REF] Santos | Mechanical Behavior of a Brazilian Off-Shore Carbonate Reservoir[END_REF]. Based on these facts, the Dimaggio-Sandler model [START_REF] Dimaggio | Material Model for Granular Soils[END_REF][START_REF] Sandler | An algorithm and a modular subroutine for the cap model[END_REF] has been chosen to test the versatility of the plastic integration algorithm adopted in the currently proposed numerical finite model.

It could be argued that the use of associative models to simulate the behavior of brittle and dilative materials may be limited. But this behavior happens at low confining pressures. At high confining pressures the cap plasticity is dominant and the material compacts. This is what happens in the reservoir and in the nearwellbore domains, except at the wellbore walls where shear and shear dilation causes breakouts. However, this shear dilation occurs at small hydrostatic stresses and thus small shear strength, presenting a low impact in the overall cavity stresses. However, the "controled dilation" of the Dimaggio-Sandler model, caused by volume expansion, cap softening and transition from shear envelope to cap envelop plasticity, can be used to reduce this effect. These reasons also justify the choice of this model in the present numerical study.

Note also that other models have also been tested in the research that supported this paper. For instance, the Lade-Kin model, in its complete version with the 12 parameters, had been implemented and tested, but it was early abandoned due to its large amount of parameters, lack of clear mechanical meaning of each individual parameter. This mechanical meaning of each parameter is important in the material characterization so that the modelers are prone to detect misleading extrapolation along non-tested load paths. Other important aspect is that this model comprehends two distinct envelopes with a discontinuity in the derivatives in their interception. The code has been devised to support sub-differentiation in a numerical manner and this capability needed to be tested. A more conventional Cam-Clay constitutive model could have been used, but would not test this functionality.

Numerical analyses for well stability and breakout simulations have been previously considered in the literature. For instance, the study in [START_REF] Aadnoy | Stability of Highly Inclined Boreholes[END_REF] shows how to determine the collapse and fracture-initiation pressure in inclined boreholes using an ideal linear elastic and isotropic model. The conclusion is that the the borehole collapse at low pressures is caused by shear failure in combination with tensile failure. Using von Mises yield condition, the increase in borehole angle does not increase the sensitivity toward collapse. The method presented in [START_REF] Somerville | The Prediction of Well Stability Using the Yield Zone Concept[END_REF] determines the failure criteria from limited unconsolidated samples. The simplicity of the approach allows realistic appreciation of the mechanical characteristics derived in laboratory. A conclusion is that unconsolidated formation appears to react mechanically in similar manner to consolidated formation. Well stability analysis for horizontal or highly inclined wells in weak formations is considered in [START_REF] Morita | Core-Based Horizontal or Highly Inclined Well Stability Analysis for Unconsolidated Formations[END_REF], showing that local failure initiation does not necessarily lead to a global borehole instability for some anisotropic rocks. The work in [START_REF] Shaohua | On the Mechanical Stability of Inclinate Wellbores[END_REF] concerns with the development of a shear failure criterion, which is a combination of the effective strength concept and the Drucker-Prager criterion, objecting to quantify the stresses at which borehole breakout will occur. Furthermore, it is shown that the mechanical stability of inclined wells can be improved by adopting optimum deviation angle and drilling direction. The work in [START_REF] Morales | Borehole Failure: Safe Drawdown Pressures and Wellbore Damage Radius[END_REF] presents numerical examples and laboratory tests to show the onset of failure and the growth of the plastic zone development and propagation during production. A method is proposed in [START_REF] Wanga | A coupled conductive-convective thermo-poroelastic solution and implications for wellbore stability[END_REF] for the calculation of linear thermo-poroelastic tangential stresses around a borehole, in a medium with isotropic permeability for the case of constant borehole temperature and pore pressure change. In [START_REF] Papamichos | Breakout stability experiments and analysis[END_REF] the authors consider the stability of hollow cylinders with holes of various geometries, which is analyzed experimentally and numerically. They showed that the elliptical breakouts are a more stable structure in terms of hole stability. The numerical simulations showed that the stability of the hole increases with increasing breakouts depth, i.e. the breakout propagation is a stable process where higher external stress is required to drive further the breakout. Laboratory tests simulating a breakout in a vertical well are presented in [START_REF] Villarroel | Breakouts: Physical and Numerical Models[END_REF] using a synthetic sandstone. The same problem is solved using the commercial finite elements software Ansys and using analytical calculations, with good correlation results.

Adaptive numerical treatment of problems by high order finite element schemes has proved to be very effective [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF], by incorporating local mesh refinement (h-adaptivity) and/or different orders of approximating polynomials on separate elements in the same mesh (p-adaptivity). These methods are usually based on a posteriori error control which has been well established for linear problems. For elastoplasticity, and using a representation in terms of variational inequalities, this analysis is more recent [START_REF] Carstensen | An optimal adaptive finite element method for elastoplasticity[END_REF].

Our adaptive approach is more empirical, and the purpose is to guarantee high quality resolution in the part of the domain where the material is subjected to plastic deformation. A refinement criterion is introduced to define regions where the material has disintegrated, which is based on the magnitude of plastic deformation second invariant. Then, a numerical procedure is constructed to automatically adapt the geometry of the computational domain by removing the region where the values of the refinement indicator is larger than a given threshold.

Implementation of hp-adaptive strategies is generally more complex than for standard finite element schemes, since the differences in polynomial orders and the presence of hanging sides introduce additional difficulties in the enforcement of continuity constraints [START_REF] Calle | Implementation of continuous h p -adaptive finite element spaces without limitations on hanging sides and distribution of approximation orders[END_REF]. The mesh refinement process for elastoplastic problems is even more challenging. The tension at an integration point in the domain depends not only on the current deformation state, but also on its loading history. Therefore, a general purpose procedure had to be developed to transfer the elastoplastic deformation history to newly created integration points.

The text is organized in the following way. In Section 2 the elastoplastic mathematical model is described, and its discretized version is summarized in Section 3. The verification of the numerical model is presented in Section 4. The developed method for the breakout simulation is documented in Section 5, and the obtained results are discussed. Some important aspects of the adopted constitutive plastic model are presented in A.

The proposed numerical model has been implemented using the NeoPZ 1 library [START_REF] Devloo | An object oriented environment for scientific programming[END_REF][START_REF] Devloo | Object Oriented Tools for Scientific Computing[END_REF], which is an object-oriented programming environment providing a framework for the development of finite element simulations. Specific classes have been developed and integrated into NeoPZ for the implementation of elastoplastic constitutive models.

Mathematical Elastoplastic Model

The mathematical elastoplastic model is composed of two parts: the equilibrium equation coupled with a constitutive initial value problem, as described in [START_REF] Neto | Computational Methods for Plasticity[END_REF].

Finite Element Formulation for the Mechanical Equilibrium Problem

Given the reservoir occupying the region Ω, under the action of body force b, traction τ applied at the boundary part Γ N ⊂ ∂Ω, and assuming no displacement elsewhere at Γ D = ∂Ω \ Γ N , the infinitesimal quasi-static mechanical problem consists in finding the displacement field u ∈ V such that the weak statement of the equilibrium equation is satisfied

ˆΩ σ : ∇ v dω - ˆΩ b. v dω = ˆΓN τ . v ds ∀ v ∈ V, (1)
where

V = { q ∈ H 1 (Ω)| q = 0 on Γ D }, σ = σ(ε)
is the solution of the constitutive initial value problem, as described in the next section, and the deformation tensor is

ε = ∇ u + ∇ u T 2 .
1 http://github.com/labmec/neopz

Constitutive elastoplastic model

According to [START_REF] Neto | Computational Methods for Plasticity[END_REF], an elastic-plastic constitutive model defines the history of an irreversible dissipative process to which the material has been submitted, which is defined by three fundamental axioms: the plasticity criterion, the flow rule, and the hardening law. They are defined in terms of a plasticity function Φ and a plastic potential function Ψ.

The history of the infinitesimal deformation tensor ε(t), t ∈ [t 0 , T ] allows to compute the evolution of the plastic deformation tensor ε p (t), the internal damage variable α(t) and a plastic multiplier . γ(t) through a system of ordinary differential equations, as described in A

In the present paper, we shall adopt the associative DiMaggio-Sandler plastic model proposed in [START_REF] Dimaggio | Material Model for Granular Soils[END_REF], whose definition is included in A.1. It was conceived for granular soils, and is currently widely applied in the oil industry to represent the behavior of rocks subject to large confinement stresses. As described in [START_REF] Santos | Mechanical Behavior of a Brazilian Off-Shore Carbonate Reservoir[END_REF], this model fits well the plasticity envelopes of a weak carbonate reservoir, which exhibits a consistent shear envelop throughout the triaxial tests under different hydrostatic loadings. It fits well even in those cases where the rock is compacted prior to shear breaking, meaning that the cap plasticity with hardening took place prior to the achievement of a perfectly-plastic shear envelope. More information can found in [START_REF] Santos | Mechanical Behavior of a Brazilian Off-Shore Carbonate Reservoir[END_REF] 3 Numerical Elastoplastic Model

The numerical model for the elastoplastic problem uses a H 1 -conforming finite element discretization of the weak statement (1) of the mechanical equilibrium boundary value problem. The computation of the constitutive law associated with each integration point uses an efficient numerical integration method for the initial value problem corresponding to associative plastic models. These two numerical schemes are summarized in the next sections.

Numerical Integration of the plastic model

An algorithm for the integration of general elasto-plastic constitutive initial value problems is summarized in A. It is composed of two major steps: the elastic trial step, and the plastic correction step. For the numerical integration of the DiMaggio-Sandler model an improved algorithm is adopted, which has been proposed in [START_REF] Cecílio | An improved numerical integration algorithm for elastoplastic constitutive equations[END_REF] for associative models. It is based on the following aspects:

• The elastoplastic model is represented in terms of the principal stresses.

• For the plastic correction step, a distance function is defined, which has to be minimized in order to determine the projection of the trial stress onto the yield surface (closest point projection). This projection is made possible thanks to the Haigh-Westergaard representation. This geometric interpretation of the plastic corrector step not only improves the numerical efficiency, but also the stability of the computations.

Finite Element Model for the Mechanical Equilibrium Problem

The finite element discretization of the equilibrium problem starts by the construction of a partition T = {K} of the computational domain, and of a finite dimensional approximation space V T ⊂ V, formed by functions which are piecewise defined in terms of polynomials based on T . Then, the approximate solution u T ∈ T is searched such that the variational formulation (1) is verified for all v T ∈ V.

Due to the vertical extension of the problem, we shall consider a two dimensional plane strain model, and restrict the computational region Ω to one quadrant. The computational domain representing the reservoir, the boundary conditions, and the geometrical mesh are illustrated in Figure 2. Because of the symmetry properties, no penetration is imposed at Γ D = {(x, 0), (0, y)| 0 < x, y < 1}. At the well-bore, the normal stress corresponds to the hydrostatic stress of the drilling fluid, and at the far field the normal stress is equal to the in-situ stress. In our simulation the gravity effects are neglected, so that b = 0.

The stress state is sensitive to the accuracy of geometry representation. Element removing and/or non smooth boundary produce concentration of non physical stresses. The use of smooth curved elements adapted to the well-bore geometry is crucial. Thus, in order to enhance the precision of the elastoplastic simulation, elements with quadratic geometry are used, and they are directional refined towards the well-bore, as shown in Figure 2.

Constant polynomial degree k = 2 is adopted for approximations in the absence of material damage. However, for simulations involving breakout, the elements that intersect the plastic region are divided and their polynomial orders are increased, as described in Section 5. The construction of shape functions and assembly of the adopted H 1 -conforming finite element approximation spaces are described in [START_REF] Devloo | Systematic and generic construction of shape functions for padaptive meshes of multidimensional finite elements[END_REF][START_REF] Lucci | Descric ¸ão Matemática de Geometrias Curvas por Interpolac ¸ão Transfinita[END_REF][START_REF] Calle | Implementation of continuous h p -adaptive finite element spaces without limitations on hanging sides and distribution of approximation orders[END_REF] for h, p and hp-adaptive curved meshes.

Verification Tests

In order to verify the implementation of the proposed finite element numerical model, we simulate the stress state around a well-bore and compare the computed solutions with the analytic Airy solution for the elastic stress state [START_REF] Cecílio | Modelagem Elastoplástica Aplicada à Simulac ¸ão Numérica da Estabilidade de Poc ¸os de Petróleo[END_REF]. The stress state is given by the DiMaggio-Sandler model, whose parameters are defined in Table 1, corresponding to a carbonate reservoir material of the Brazil Offshore Santos Basin. Its characterization has been presented in paper [START_REF] Santos | Mechanical Behavior of a Brazilian Off-Shore Carbonate Reservoir[END_REF]. The well radius is r wb = 0.10795 m, the reservoir radius is r e = 3 m, and the effective in situ tension components are maximum horizontal stress in the y-direction S H = -62.1 M P a, minimum horizontal stress in the x-direction S h = -45.9 M P a, and vertical stress S v = -48.2 M P a.

The mesh configuration for this test is obtained by 40 subdivisions in the radial direction, and 80 circular ones, and the polynomial degree is k = 2. The results are documented in Figure 3 by the graphs of the values of the plasticity function Φ along the well-bore circular curve. Three cases are represented, corresponding to different values of well-bore fluid pressures P f . The corresponding values of the yield function Φ(I 1 , √ J 2 , L), the plastic deformation second invariant J 2 (ε p ), the first invariant of the stress tensor I 1 (σ) and the second invariant of the deviatoric stress tensor J 2 (σ) are graphically represented in Figure 4-Figure 6 Parameter Case 1: P f = 28.9 M P a In Figure 3 (top left side), the plot for the yield function along the well-bore circular curve is such that Φ < 0, which means that the problem is purely elastic. The values of Φ(I 1 , √ J 2 , L), J 2 (ε p ), I 1 (σ) and J 2 (σ) are graphically represented in Figure 4. As expected, Φ < 0 and J 2 ( p ) = 0.

Case 2: P f = 23.4 M P a The results for this case are graphically represented in Figure 5. The lower fluid pressure causes the occurrence of a region with plastic deformation around the well-bore, where Φ = 0. The extent of the region corresponds to approximately 1.3% of the area of the reservoir model. The comparison with the analytic Airy solution in Figure 3 (top right side) shows that the plastification extends to an angle of 46.5 • . Case 3: P f = 19.5 M P a Decreasing the fluid pressure even more, the region with plastific deformation grows, as shown in Figure 6, occupying around 4.28% of the modeled reservoir area, and extending to a 72.12 • angle at the well-bore curve, as illustrated in Figure 3 (bottom side). As suggested by [START_REF] Zoback | Wellbore breakouts and in situ stress[END_REF], a safe threshold is 60 • . Thus, according to this criterion, the stress state of Case 3 corresponds to a dangerous scenario for well stability. In order to evaluate the effect of discretization resolutions on the numerical solution, the method is applied to the reservoir configuration of Case 3 using three different mesh refinements, and polynomial degrees k = 1, 2 and 3. The first mesh is obtained by 20 subdivisions in the radial direction, and 40 circular ones (Mesh 1). The second (Mesh 2) and third (Mesh 3) finer meshes consists of subsequent uniform subdivisions of the first one. Defining the plasticized region where J 2 (ε p ) > 0, its area for different discretizations is represented graphically in Figure 8. Figure 7 displays the value of J 2 ( p ) for the solutions on Mesh 2 and Mesh 3, with k = 2. It can be observed that the plastified area converges as the mesh is refined either in h or in polynomial degree.

Breakout Simulation

The breakout phenomenon is characterized by the well-bore mechanical failure and loss of material in the regions where the shear strains are critical. The elastoplastic simulation allows to identify the regions where this is most likely to happen.

In the well drilling area of oil and gas industry, the information of damage extension is at a premium. Well design and stability monitoring while drilling hence rely on generally poor, interpolated information, mainly obtained from correlation from welllogs. On the other hand, rock is heterogeneous, brittle, anisotropic, sometimes is load-rate-dependent, etc. There must be a trade off between model sophistication and applicability. This is the main reason for requiring models based on "mechanically meaningful" parameters, so that correlation from physical properties inferred from the welllogs may be devised, as well as any misjudgements may be detected at early stages (parameter-definition stages). Having that in mind, it should become clear that the numerical solution on that backgroung must not be taken literally into consideration. The plastified zone does not necessarily represent a damaged zone which remains attached to the formation, nor continuously. There is insufficient amount of information to deploy real damage models that determine the volume of fragmented rock debris. There is also insufficient information to state whether the excavation could still be geometrically stable.

The plastified zone or opening angle should be regarded as imperfect quantitative indications of a mechanical damage zone. From the drilling experience, global excavation instability occurs rarely and is always preceded by local wall instabilities (breakout). The formation of breakouts poses more problems to the well cleaning than to the cavity stability, as the breakouts may occur over a large extent along the wellbore axis. The volume of rock debris produced by breakout may easily surpass the amount of cuttings from the drillbit. Pressure fluctuations in the excavation associated to the drillstring vibration ease the freeing of damaged material. The settlement of large amounts of debris around the drillstring may cause it to get stuck and prevent further drilling with some other serious consequences.

Thus, each drilling company rely on ultimate plastified area or on opening angle indications to drive their wellbore stability analysis. However, they usually rely on the simple Mohr-Coulomb shear-envelope-basedplasticity models to evaluate those indicators. Although it captures the shearing processes at the wellbore wall, they are unable to reproduce the amount of compaction experienced by the rock mass under depletion that drives the evolution of far-field stresses.

Our proposal is to use the plastic deformation second invariant J 2 (ε p ) to define a criterion for the regions where the material has disintegrated. As described in the next sections, the numerical method automatically adapts the geometry of the computational domain by removing the region where the values of J 2 (ε p ) are larger than a given threshold. Then, a moving-mesh algorithm is adopted to simulate how the cavity failure is supposed to evolve based on the complete detachment of the plastified zone material.

Automatic Procedure for Adjusting the Geometry of the Domain

An attempt to adapt the geometry simply removing elements with excessive plastic deformation is not a viable solution: the removal of the elements creates nonphysical re-entrant corners causing high stress concentrations that are non physical.

A different approach is considered based on the hypothesis that the isolines of the proposed breakout criterion can be closely matched by ellipses whose center is the well-bore center. This same geometry can also be observed in experiments conducted by [START_REF] Zoback | Wellbore breakouts and in situ stress[END_REF]. Based on this principle, a numerical method has been developed [START_REF] Cecílio | Modelagem Elastoplástica Aplicada à Simulac ¸ão Numérica da Estabilidade de Poc ¸os de Petróleo[END_REF] to adjust the axes of an ellipse to points with constant values of the breakout criterion J 2 (ε p ).

A motivation for this choice of criterion is that the plastic deformation is characterized by two components, one hydrostatic and the other deviatoric. One can conjecture that the hydrostatic deformation does cause the material to fail. The shear deformation, on the other hand, causes the sliding between the microscopic plains possibly leading to rupture. Hence the value of J 2 (ε p ) has been chosen as a criterion to quantify the tendency of the material to disintegrate.

The computation of the axes of the ellipse corresponding to a given value of J 2 (ε p ) is performed using a standard least square algorithm. Precisely, a list with coordinates (x i , y i ) of points with constant value of

J 2 (ε p ) are computed. The error function E(a, b) = i x i -a 1 - y 2 i b 2
is defined to measure the deviation of the points (x i , y i ) from a general ellipse with focus at (0, 0). Then, the least squares method is used for the determination of the values a and b minimizing this error function.

Once the ellipse equation is found, the mesh nodes on the contour of the well-bore are moved to fit the position of the ellipse. The remaining nodes of the mesh that are moved to maintain the smoothness of the mesh.

Refinement Process on the Plastic Zone

The mesh refinement process for elastoplastic problems is challenging. The tension at an integration point in the domain depends not only on the current deformation state, but also on its loading history. For meshes without adaptation, this history can be obtained by storing the of the nodal displacements. In the case of mesh movement, the displacements themselves do not allow to reconstruct the history of deformation. During the refinement process, integration points are created and others are deleted. The adopted technique to refine meshes arbitrarily is to implement a general procedure that allows to compute the solution at an arbitrary point in a mesh. At each refinement step a new mesh object is created and the deformation history is obtained by computing the history of deformation at a given coordinate by computing the incremental deformation at all previous meshes (e.g. solutions). For more details see [START_REF] Cecílio | Modelagem Elastoplástica Aplicada à Simulac ¸ão Numérica da Estabilidade de Poc ¸os de Petróleo[END_REF]. The implementation has two fundamental steps:

1. Identify the element and its coordinates in the master element domain corresponding to an (x, y) coordinate.

2. Submit the newly created material point to the deformation tensor evolution according to the solution of the previous meshes.

Refinements h and p in the plasticized zone are illustrated in Figure 10.

Results

The geometry adjustment and the hp mesh refinement are illustrate by numerical results for the confinement state configuration in Case 3, described in Section 4. This choice is justified because, among the three considered cases, this is the worst case scenario, since the breakout damage has a breakout angle greater than 60 • , as detailed in the reference [START_REF] Zoback | Wellbore breakouts and in situ stress[END_REF]. Two parameters are defined: Ref HP and AdjBK. In a first step, elements that have J 2 (ε p ) > Ref HP are searched in the mesh, were an hp mesh refinement is performed. Secondly a geometric adjustment criteria is specified to search a line with constant values of J 2 (ε p ) = AdjBK, which is adjusted by an ellipse. Therefore, four geometry adjustment stages are performed. Precisely, 1. Given an initial well-bore geometry, the problem is initially solved with p = 2, as described in Section 4.

2.

A mesh adjustment procedure is performed, subdividing each element that has J 2 (ε p ) ≥ Ref HP , and increasing it's polynomial order to p = 3.

3. The geometry is adjusted based on the criterion AdjBK and the problem is solved in the original mesh that is restricted to the new geometry.

4. Return to step 1.

The results in Figure 11 are for an initial mesh with 40 × 80 elements (Test 1). The values for J 2 (ε p ) are shown at four instances of the breakout simulation. The thresholds AdjBK = Ref HP = 0.0001 were chosen such that the relevant parts of the domain where plastification occurs are refined. For comparison, the curves corresponding to J 2 (ε p ) = 0.0001 obtained at the four stages are presented in all figures, and plotted in white color. At the first instance, shown in Figure 11 (a), the results are for the original geometry region, and the hp adapted mesh is shown. At this initial instant, a maximum second invariant plastic deformation value J 2 (ε p ) = 0.000829 is observed. In a next simulation instant, represented by Figure 11 (b), the geometry adjustment is performed, by the removal of the plasticized area delimited by the ellipse that has been adjusted to the points where J 2 (ε p ) = AdjBK (first white isoline seen in Figure 11 (a)). At this instant, a maximum value J 2 (ε p ) = 0.001669 is observed. Again, a white isoline is used to mark the points where J 2 (ε p ) = AdjBK. Next, the plasticized region delimited by this isoline is removed and the simulation returns to stage 1. Figure 11 (c) indicates the next result, after the hp refinement and running the problem in a new mesh. It is observed a maximum value of J 2 (ε p ) = 0.003159, which gives an evidence for a plastic deviatoric deformation increasingly higher. Figure 11 (d) shows the result of the last procedure instance. The last maximum value J 2 (ε p ) = 0.006338 is more than twice the value observed in the previous step and the plasticized area is considerably reduced.

Two tendencies have been noted. One is the tendency of the plasticized zone to move towards the inferior symmetry axis, with decreasing angle of plastic deformation extent. The other evidence is that the plasticized area decreases. The data shown in Table 2 for removed and plasticized areas, damage extension and breakout opening angles at the four breakout simulation instants described before confirm these tendencies. Note that in well-bore projects the angle of plastic deformation extent is an important parameter, and it is usually computed based on linear elastic stress states, and limited to a predefined value.

Step removed area (×10 -3 m 2 ) plasticized area (×10 -3 m 2 ) damage area (×10 2: Breakout evolution using different mesh configurations: with hp adapted meshes with initial uniform partition of 40x80 and polynomial order p = 2 (Test 1) and 80x160 elements, and polynomial degree p = 2 (Test 2), and without adaptativity, using the fine mesh of 80x160 elements and polynomial degree p = 3 (Test 3).

In Table 2 simulation results are presented for two other discretization configurations. Test 1 refers to an adaptive simulation starting from an initial uniform mesh with 40 × 80 elements and k = 2. Test 2 refers to an The results of the removed, plasticized and damage areas are in close agreement in all the three strategies. This table shows that the extent of plastic deformation and/or breakout is stable with respect to the choice of initial mesh and polynomial order of approximation.

Conclusions

The DiMaggio-Sandler elastoplastic model has been applied to the study of well-bore stability, applying hpadaptive finite element approximations combined with mesh movement. The robustness of this elastoplastic model is greatly improved by parametrizing the plastic yield surface and performing proper tensor analysis. A failure criterion is proposed assuming that the material will disintegrate due to excessive shear strain. Based on this criterion, the quality of the finite element simulation is guaranteed by applying h and p refinements in the region were plastic deformation exceeds a given old. For the newly created adaptive meshes, the complete deformation history is applied to each integration point. Stacks of meshes are used to keep track of this deformation history. After this adaptive simulation, an automatic elliptic geometry fitting is applied to the region of excessive deformation, and the corresponding material is removed, by using mesh movement. One of the results of our work is that using the automatic successive application of the breakout damage criterion, and after removing the area with excessive plastic deformation, realistically looking geometries are obtained, as compared with experimental results reported in the literature, the area of plastic deformation migrating towards the horizontal line, with decreasing angle of plastic deformation extent.

Figure 1 :

 1 Figure 1: Illustration of an off-shore well drilling process: fluid is pumped and circulated inside the perforation column for lubrication, transport of debrits and stress equilibrium purposes.

Figure 2 :

 2 Figure 2: Computation domain and boundary conditions.

Figure 3 :

 3 Figure 3: Plasticity function Φ (M P a) computed at the circular well-bore curve by the Airy solution (red lines) and by the numerical finite element solution (dotted lines) for the Case 1, Case 2, and Case 3.

Figure 4 :

 4 Figure 4: Case 1: Φ(I 1 , √ J 2 , L), and J 2 (ε p ) (top side); I 1 (σ), and J 2 (σ) (bottom side).

Figure 5 :

 5 Figure 5: Case 2: Φ(I 1 , √ J 2 , L), and J 2 (ε p ) (top side); I 1 (σ), and J 2 (σ) (bottom side).

Figure 6 :

 6 Figure 6: Case 3: Φ(I 1 , √ J 2 , L), and J 2 (ε p ) (top side); I 1 (σ), and J 2 (σ) (bottom side).

  (a) 40 x 80 elements: plasticized area = 1.12 × 10 -3 m 2 . (b) 80 x 160 elements: plasticized area = 1.13 × 10 -3 m 2

Figure 7 :

 7 Figure 7: Value of J 2 ( p ) associated to the solutions obtained with mesh resolutions and k = 2.

Figure 8 :

 8 Figure 8: The effect of the mesh resolution and polynomial degree on the plasticized area.

Figure 9 :

 9 Figure 9: Boundary adjustment using an ellipse shape.

  (a) p-refinement. (b) h-refinement

Figure 10 :

 10 Figure 10: Illustration of h and p mesh refinements on the plasticized zone.

Figure 11 :

 11 Figure 11: Breakout simulation: automatic geometric fitting to the plastic zone.

Table 1 :

 1 Parameters for DiMaggio-Sandler elastoplastic model.

		Value
	E (elasticity modulus)	29269 M P a
	ν (Poisson coefficient)	0.203
	A	152.54 M P a
	B	0.0015489 M P a -1
	C	146.29 M P a
	D	0.018768 M P a -1
	R	0.91969
	W	0.06605

Acknowledgments

The authors thankfully acknowledges financial support from ANP-Brazilian National Agency of Petroleum, Natural Gas and Biofuels (grant SAP4600333146). P. R. B. Devloo and S. M. Gomes also thank financial support from CNPq -the Brazilian Research Council ((305425/2013-7, and 304029/2013-0, respectively). The authors also thank O. Durán, T. Quinelato, and L. Marques for the help in setting some NeoPZ configurations for the current simulations.

A Constitutive elastoplastic model

As studied in [START_REF] Neto | Computational Methods for Plasticity[END_REF], the total deformation tensor ε is divided into two components: ε = ε e + ε p , an elastic strain ε e and a plastic strain ε p . The plastic strain is related to the history of irreversible dissipative processes to which the material was submitted based on three fundamental axioms: the plasticity criterion, the flow rule, and hardening law.

• Yield Criterion. Describes the transition between the elastic and plastic domains through a plasticity function Φ, assuming non-positive values if the stress state is elastic and null values if the stress state reaches its plastic limit.

• The plastic flow rule. The plastic flow rule is defined by a plastic potential function Ψ, which specifies how the plastic deformation tensor ε p evolves in a plasticity process

.

γN, in which N = ∂Ψ/∂σ is the flow direction.

• Hardening law. Specifies how the internal damage variable evolves γH, where H is defined as the hardening modulus.

Following [START_REF] Neto | Computational Methods for Plasticity[END_REF], the integration procedure to go from a time instant t n to the next one t n+1 is usually performed by two major steps. The first step is a purely elastic predictor process (elastic trial step), with ∆γ = 0. The trial elastic strain and internal variables are defined as ε e trial = ε e,n + ∆ε and α trial = α n . σ trial is calculated according to ε e trial , and the corresponding Φ(σ trial , A) is computed. If Φ(σ trial , A) ≤ 0, a valid solution to the system is reached, and the variables are updated by the trial values. Otherwise, a plastic corrector step (also known as plastic return-mapping scheme) is performed reformulating the incremental problem solving for ε e,n+1 , α n+1 and ∆γ satisfying

Next, the plastic strain is updated ε p,n+1 = ε p,n + ∆ε -∆ε e .

A.1 DiMaggio-Sandler plastic model

The plastic model adopted in the present study is the DiMaggio-Sandler model [START_REF] Dimaggio | Material Model for Granular Soils[END_REF]. Its plasticity function Φ is piecewise defined by the yield function F 1 (I 1 , √ J 2 ), and by the cap elliptical function F 2 (I 1 , √ J 2 , L), where I 1 and J 2 are the first invariant of the stress tensor and the second invariant of the deviatoric stress tensor, respectively, and L is a parameter whose value is generally negative. Precisely,

where

and F f (s) = A -C exp(Bs). A typical 2D profile of the yield surface is plotted in Figure 12. Note that |L| specifies the distance of the origin to the center of the ellipse representing the cap of the yield surface. The relation between the value of L and the plastic volumetric strain p v is defined by:

where

The constants A, B, C, D, R and W are material properties, which need to be obtained from laboratory tests.

The DiMaggio-Sandler model is associative such that Ψ = Φ. The corresponding hardening modulus H, and the flow direction N are related by H = tr(N).

Using this property, the equation for the damage variable at the plastic correction step can be reformulated as In terms of its Haigh-Westergaard stress representation σ * = [σ * 1 , σ * 2 , σ * 3 ] T , the yield surface is σ * 1 -axially oriented and is formed by two parts: one perfectly plastic surface S * 1 , where F 1 = 0, and a cap S * 2 (L), corresponding to F 2 = 0, as represented graphically in Figure 13. The plastic return-mapping considers these two cases, and the return-mapping scheme is reduced to the resolution of a single non-linear equation.