A El

Ouazani Tuhami

S Feld-Payet

S Quilici

N Osipov

J Besson

A two characteristic length nonlocal GTN model: application to cup-cone and slant fracture

Keywords: Nonlocal GTN model, Characteristic length, Ductile failure, Cup cone fracture, Slant fracture

Ductile failure prediction is essential to avoid loss of structural integrity or to control crack propagation during forming processes. One of the main difficulties is the prediction of complex crack paths. This study proposes a nonlocal extension of a local Gurson-Tvergaard-Needleman (GTN) ductile damage model at finite strains able to capture cup-cone or slant fracture. The proposed model is based on an implicit gradient formulation which enables to solve the problem of spurious strain and damage localization. The model integrates two different material characteristic lengths, which are used to separately regularize damage by void growth and damage by damage nucleation. After parameter fitting using data for a line pipe steel, conditions to obtain converged solutions are studied, which can be used to select the mesh size in the localization band. With the appropriate mesh design, it was possible to study the effect of the characteristic lengths on the formation of cup-cone fracture and slant fracture. It is observed that, for a given specimen size, larger characteristic lengths favor flat crack advance. Similarly, for given material lengths, size effects can be predicted with small specimens being more prone to flat fracture. This paves the way to a more direct determination of material lengths by using homothetic specimens so as to obtain different crack paths.

Introduction

Predicting ductile failure is a major challenge when trying to prevent catastrophic failure of structures or to control crack propagation in the case of metal forming. One of the main challenges is the prediction of complex crack paths in metals such as cup-cone or slant fracture. To achieve this goal, models able to robustly represent the local degradation phenomena (damage nucleation, growth and coalescence) are needed. Among such models is the micromechanics based Gurson-Tvergaard-Needleman (or GTN) model [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]. Based on this seminal framework, numerous extensions have been proposed, which are reviewed e.g. in [START_REF] Tvergaard | Material failure by void growth to coalescence[END_REF][START_REF] Besson | Damage of ductile materials deforming under multiple plastic or viscoplastic mechanisms[END_REF][START_REF] Benzerga | Ductile Fracture by Void Growth to Coalescence[END_REF][START_REF] Pineau | Failure of metals I: Brittle and ductile fracture[END_REF]. They are able to account for void shape, void orientation, plastic anisotropy, strain rate sensitivity, coalescence by internal necking. . . All these models must lead to softening up to a point where the material entirely loses its load-carrying capacity so as to be able to model fracture. These softening models were first developed within the usual local framework and thus suffered from the well-known numerical issues such as mesh size and mesh orientation dependence [START_REF] Rousselier | Non-Linear Fracture Mechanics[END_REF][START_REF] Liu | Mesh-dependence and stress singularity in finite element analysis of creep crack growth by continuum damage mechanics approach[END_REF][START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF]. Due to the loss of ellipticity of the equations in statics, the solution to the mechanical problem is no longer unique. Overcoming these issues is fundamental for the reliability of fracture simulations. Using a finite strain framework is also required as ductile fracture is always accompanied by large deformations. Several solutions have been proposed to overcome mesh dependency. A first solution is based on a local enrichment by embedding a finite thickness band [START_REF] Huespe | A finite thickness band method for ductile fracture analysis[END_REF]. In that case, the band is introduced when loss of ellipticity is detected within one element. The normal to the band is determined following Rice's bifurcation analysis [START_REF] Rice | The localisation of plastic deformation[END_REF]. Another solution is based on implicit gradient methods [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF][START_REF] Geers | Strain-based transient-gradient damage model for failure analyses[END_REF] which facilitate the use of integral methods as originally proposed in (Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Bazant | Non local continuum damage. localization, instability and convergence[END_REF]). These methods have been used to model ductile fracture for metals within a finite strain framework [START_REF] Enakoutsa | Numerical implementation and assessment of a phenomenological nonlocal model of ductile rupture[END_REF][START_REF] Mediavilla | A nonlocal triaxiality-dependent ductile damage model for finite strain plasticity[END_REF][START_REF] Linse | Simulation of crack propagation using a gradient-enriched ductile damage model based on dilatational strain[END_REF][START_REF] Hütter | Simulation of ductile crack initiation and propagation by means of a non-local Gurson-model[END_REF][START_REF] Javani | Three-dimensional finite element modeling of ductile crack initiation and propagation[END_REF][START_REF] Seupel | On the identification and uniqueness of constitutive parameters for a non-local GTN-model[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] although they were initially developed for quasi-brittle failure. Using an implicit gradient formulation, the local constitutive equations are preserved at the cost of small adaptations. This makes this formulation particularly attractive in the case of the GTN model, as its initial micro-mechanical foundations are then preserved. Multiple internal lengths were introduced in [START_REF] Nguyen | A nonlocal approach of ductile failure incorporating void growth, internal necking, and shear dominated coalescence mechanisms[END_REF]. This model uses an implicit gradient formulation written in the initial configuration. An internal length was associated with each damage evolution mechanism: (i) void growth based on the classical GTN model, (ii) internal necking governed by a heuristic extension of the Thomason model (Thomason, 1985b,a) based on the maximum principal stress, (iii) shear-dominated coalescence mechanisms controlled by the maximum shear stress [START_REF] Torki | On void coalescence under combined tension and shear[END_REF]. Although the model proposes different internal length scales, the same characteristic length was assigned to the three nonlocal variables in the presented simulations. Micromorphic models [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF] can also be used to solve problems related to the loss of ellipticity. They use local state variables and their "micromorphic" counterparts. The free energy of the system depends on the gradient of the micromorphic variables and on coupling terms. Ductile fracture was modeled using this class of models in [START_REF] Brepols | Gradient-extended two-surface damage-plasticity: Micromorphic formulation and numerical aspects[END_REF][START_REF] Diamantopoulou | Micromorphic constitutive equations with damage applied to metal forming[END_REF]. In particular, models based on a microdilatational theory [START_REF] Huetter | A micromechanical gradient extension of gurson's model of ductile damage within the theory of microdilatational media[END_REF] appear to be well suited to represent ductile damage by void growth.

More recently, nonlocal gradient enhanced energy (GEE) models, first developed for quasi-brittle fracture [START_REF] Lorentz | A variational formulation for nonlocal damage models[END_REF] were also used to model ductile failure [START_REF] Zhang | Ductile damage modelling with lockingfree regularised gtn model[END_REF][START_REF] Chen | Properties of a nonlocal GTN model within the context of small-scale yielding[END_REF]. They only use the gradient of a local state variable. A decomposition-coordination technique is used to treat the non-locality. The variable of interest is duplicated: a first instance is used at the (global) scale of the structure while a second instance is used at the (local) constitutive law level. As both variables represent the same field, they should be equal. A Lagrange multiplier is introduced to ensure this equality weakly. Finally, the regularization of ill-posed problems can also be achieved by coupling the elastoplastic models with phase-field formulations. This type of formulation was initially introduced for brittle failure [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Hofacker | Continuum phase field modeling of dynamic fracture: variational principles and staggered FE implementation[END_REF][START_REF] Miehe | A phase field model for rateindependent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Heider | A phase-field modeling approach of hydraulic fracture in saturated porous media[END_REF]. However, many extensions to ductile fracture have been introduced for 2D simulations [START_REF] Ambati | Phase-field modeling of ductile fracture[END_REF][START_REF] Aldakheel | A modified gurson-type plasticity model at finite strains: formulation, numerical analysis and phase-field coupling[END_REF], for 3D simulations [START_REF] Aldakheel | A modified gurson-type plasticity model at finite strains: formulation, numerical analysis and phase-field coupling[END_REF][START_REF] Ambati | A phase-field model for ductile fracture at finite strains and its experimental verification[END_REF][START_REF] Borden | A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects[END_REF][START_REF] Miehe | Phase field modeling of fracture in multi-physics problems. Part II. Coupled brittle-toductile failure criteria and crack propagation in thermo-elastic-plastic solids[END_REF][START_REF] Miehe | Phase field modeling of ductile fracture at finite strains: A variational gradient-extended plasticity-damage theory[END_REF][START_REF] Hu | A variational phase-field model For ductile fracture with coalescence dissipation[END_REF] and also for 3D simulations with remeshing techniques [START_REF] Eldahshan | Phase field modeling of ductile fracture at large plastic strains using adaptive isotropic remeshing[END_REF]. More recently, a porous ductile model using a phase-field formulation was introduced in which the critical energy is decomposed into elastic and plastic contributions where the plastic part is described using a GTN model [START_REF] Dittmann | Phase-field modeling of porous-ductile fracture in non-linear thermoelasto-plastic solids[END_REF].

This work proposes a nonlocal extension of the GTN material model (based on [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]) using an implicit gradient framework incorporating two nonlocal characteristic lengths to capture complex ductile failure patterns, such as cup-cone or slant crack paths. The most common issue with nonlocal formulations is the choice of the characteristic lengths. In this paper, comparisons with experimental data on steel specimens first enable obtaining an order of magnitude for the model parameters. Then two length scales are distinguished in the model to account for the different spacings that exist between material defaults responsible for damage nucleation and void growth. In the case of ferritic steels used in modern pipelines, two damage mechanisms exist: (i) void growth from MnS and oxides and (ii) damage nucleation at iron carbides (Fe 3 C), which occurs at high strains. Mns and oxides early debond from the metallic matrix so that they can be considered as initial voids. Similar damage behavior is e.g. observed in aluminum alloys where coarse particles are Fe-rich inclusions and small particles are strengthening dispersoids [START_REF] Bron | Ductile rupture in thin sheets of two grades of 2024 aluminum alloy[END_REF]. In order to guide the choice of these characteristic lengths, an original study is conducted on their effect on the occurrence of cup-cone and slant fracture. This work is structured as follows. First, in section 2, the local version of the model is briefly recalled, and its nonlocal modification is introduced. The proposed finite strain formulation is also introduced within an updated Lagrangian framework. In section 3, a parameter fit is carried out for both the local and nonlocal versions of the GTN model based on existing literature data. In section 4, the independence of the results to both mesh size and mesh orientation is checked using the same characteristic length for both mechanisms. The relationship between the numerically obtained localization band width and the characteristic lengths is also investigated. This allows the selection of appropriate mesh sizes for given characteristic lengths. Finally, in section 6, an investigation is presented of the effect of different characteristic lengths for damage nucleation and void growth on the formation of cup-cone, and slant crack paths and the results are translated in terms of size effects.

Material models and their framework

In the first part of this section, the local version of the considered GTN model used throughout the paper is briefly recalled. From this basis, the proposed modifications to obtain a nonlocal model with two characteristic lengths are developed. A corotational finite strain formulation is used to deal with large strains (see section 2.2). This allows to use a simple small strain like formulation for the constitutive model based on an additive strain decomposition of the strain tensor.

Material models

The used local model is based on the GTN model presented by [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] but with a different nucleation function using the accumulated plastic strain κ as opposed to the original model which uses the growth porosity f g .

Local GTN model

One assumes the additive decomposition of the strain tensor (ε) into an elastic (ε e ) and a plastic (ε p ) part. The elastic strain tensor and the stress tensor are related using Hooke's law:

σ = E : ε e (1)
where E is the fourth-order elasticity tensor. Work hardening is assumed to be isotropic, and the flow stress of the sound material (R(κ)) is expressed as a function of the accumulated plastic strain κ assuming a power law:

R(κ) = K(e 0 + κ) n (2)
e 0 , K and n are parameters to be identified (any other function could indeed be used). The yield surface is then expressed following [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] as:

Φ = σ -R(κ)
where σ is an effective stress measure depending on both the stress tensor and the void volume fraction, which is the only damage parameter. In the case of the GTN model, the effective scalar stress is implicitly defined as:

G(σ, f, σ ) = σ 2 eq σ 2 + 2q 1 f cosh q 2 2 σ kk σ -1 -q 2 1 f 2 def = 0 (3)
where σ eq is the von Mises stress invariant1 of the Cauchy stress, σ kk designates its trace, q 1 and q 2 are two parameters and f is a function of porosity. Function f is defined as in [START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF]:

f = f if f ≤ f c f c + δ(f -f c ) otherwise (4)
where f c represents the porosity at the onset of void coalescence and δ ≥ 1 is a parameter which represents the increased deleterious effect of porosity above f c . The damage variable f can be itself decomposed into two parts which reflect the fact that damage is caused by damage nucleation (f n ) and void growth (f g ), so that the total damage is equal to f = f g + f n . The evolution of these damage variables is given by the following set of equations:

ḟg = (1 -f g )trace( εp ) mass conservation (5) ḟn = A n (κ) κ strain controlled nucleation (6)
where A n is a coefficient representing the damage nucleation rate. It is expressed as a function of κ so as to represent strain-controlled nucleation [START_REF] Chu | Void nucleation effects in biaxially stretched sheets[END_REF][START_REF] Zhang | A complete Gurson model approach for ductile fracture[END_REF]. Note that eq. 5 slightly differs from the original one where (1 -f g -f n ) is used and not (1 -f g ). As MnS inclusions early debond from the matrix, it is unnecessary to describe the nucleation of voids on these inclusions. In this work, nucleation corresponds to damage created on iron carbides which are much smaller than the MnS inclusions. It is assumed that this creates very small voids but high damage. For that reason, the equation for ḟg is modified as growth is only attributed to large voids originated at MnS inclusions. In practice, this modification has little effect on the overall behavior. The plastic strain rate tensor εp is obtained using the normality rule:

εp = (1 -f ) κ ∂Φ ∂σ = (1 -f ) κ∂σ ∂σ = (1 -f ) κn (7)
where n designates the normal to the yield surface. Using this expression, the following equivalence between the macroscopic (left handside) and the microscopic (right handside) plastic dissipations is obtained as:

εp : σ = (1 -f ) κσ (8)
as σ is an homogeneous function of order 1 of σ.

It is assumed that the considered material can exhibit a slight strain rate dependence, so that κ is expressed as:

κ = F(Φ) = ṗ0 σ -R σ 0 n (9)
where . is the positive part function and n, ṗ0 and σ 0 are material parameters.

Nonlocal GTN model

The proposed extension of this GTN model uses the implicit gradient methodology proposed in [START_REF] Peerlings | Some observations on localisation in non-local and gradient damage models[END_REF][START_REF] Geers | Strain-based transient-gradient damage model for failure analyses[END_REF][START_REF] Engelen | Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour[END_REF] to regularize two state variables: the volume variation ω = trace(ε p ) and the effective accumulated plastic strain κ. Their nonlocal counterparts are referred to as ω and κ. Their evolution within the considered material body Ω is governed by the following Helmholtz-type equations:

ω -l 2 ω ∆ω = ω in Ω (10) κ -l 2 κ ∆κ = κ in Ω (11)
with the following natural boundary conditions:

∇ω. n = 0 on ∂Ω (12) ∇κ. n = 0 on ∂Ω (13)
where ∂Ω and n respectively designate the boundary of the body Ω and its outer normal vector.

Let us note that the gradients are calculated on the current configuration and not on the initial configuration as considered by [START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF][START_REF] Nguyen | A nonlocal approach of ductile failure incorporating void growth, internal necking, and shear dominated coalescence mechanisms[END_REF]. Besides, this model introduces two characteristic lengths in equations 10 and 11, l ω and l κ , respectively associated to void growth and nucleation mechanisms.

The nonlocal variables ω and κ are then used to formulate the evolution of the damage variables f g and f n as:

ḟg = (1 -f g ) ω (14) ḟn = A n (κ) κ (15)
These two equations now replace equations 5 and 6. All other equations remain unchanged. In particular, the hardening law R(κ) remains identical, which would not be the case with micromorphic models [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF] or gradient enhanced energy models [START_REF] Chen | Properties of a nonlocal GTN model within the context of small-scale yielding[END_REF] as it is expressed as a function of κ in those cases. In the present case, the flow stress then remains unaffected by the gradient of κ.

Implicit local resolution

The state variables describing the material behavior are then the elastic strain tensor (ε e ), the accumulated plastic strain (κ), the porosity due to void growth (f g ), the porosity due to damage nucleation (f n ) and the volume variation ω. Their evolution laws can be expressed as functions of the rates of the input variables ε, ω, κ as:

εe = ε -(1 -f g ) κn (16) κ = F(Φ) (17) ḟg = (1 -f g ) ω (18) ḟn = A n ω (19) ω = trace ((1 -f g ) κn) = (1 -f g ) κtrace(n) (20)
The implicit resolution of this system of equations for each integration point, detailed in Appendix A, requires its integration over a finite time step ∆t. This leads to a system of equations relating the residuals of these equations to the increments of the state variables ∆V S = (∆ε e , ∆κ, ∆f g , ∆f n , ∆ω), for a given set of the increments of input variables ∆V IN = (∆ε, ∆ω, ∆κ).

The resolution of this system directly gives the increments of output variables [∆σ, ∆ω, ∆κ]. It then becomes possible to numerically compute the consistent tangent matrix by considering that any given small variation of the input variables δV IN leads to a modification of the state variables so that the residual equations remain null (see details in Appendix B). The consistent tangent matrix can be expressed as a block-matrix as:

∆V OUT =
K mat = ∂∆V OUT ∂∆V IN =         ∂∆σ ∂∆ε ∂∆σ ∂∆ω ∂∆σ ∂∆κ ∂∆ω ∂∆ε ∂∆ω ∂∆ω ∂∆ω ∂∆κ ∂∆κ ∂∆ε ∂∆κ ∂∆ω ∂∆κ ∂∆κ         (21)
The terms of this matrix are required to express the elementary stiffness matrix of the associated finite elements.

Finite strain framework

Both local and nonlocal GTN models must be formulated within a finite strain framework. This is done using a rotating frame concept.

Corotational framework. To do so, a corotational formulation [START_REF] Sidoroff | Some issues about anisotropic elastic-plastic models at finite strain[END_REF] is used. The corotational frame is defined using the rotation tensor Q such that:

Q = W .Q with Q(t = 0) = 1 (22)
where the W tensor designates the skew-symmetric part of the velocity gradient tensor L. The symetric part of the velocity gradient tensor is denoted D so that: L = D + W . This velocity gradient tensor is defined from the transformation tensor F = ∂ x/∂ X (where X is the position of a material point in the initial configuration and x its position in the current configuration) and its derivative as: L = Ḟ .F -1 . Knowing the rate of deformation D, one gets the material deformation rate through the following expression:

ε = Q T .D.Q (23)
Knowing the corotational Cauchy stress (i.e., the Cauchy stress expressed in the coordinate system that rotates with the material), the Cauchy stress in the unrotated frame is given by:

Σ = Q.σ.Q T (24)
Let us note that the rotated rate of σ is the Jaumann rate of Σ.

Finite element formulation.

Considering cases where the load is applied sufficiently slowly so that the inertial forces can be neglected, the weak form of the equilibrium equations then reads:

∀ v Ω Σ : L dΩ = Ω Σ : D dΩ = W ext , with L = ∂ v ∂ x ( 25
)
where v is a virtual velocity field, W ext represents the virtual power of external forces and D is the symmetric part of L .

In the case of a nonlocal formulation, based on equations 10, 11 and 12, the weak form corresponding to the nonlocal variables is expressed as:

∀ ẇ Ω (w -w) ẇ + l 2 w ∇w. ∇ ẇ dΩ = 0 ∀ ṗ Ω (p -p) ṗ + l 2 p ∇p. ∇ ṗ dΩ = 0 (26)
where ẇ and ṗ designates a virtual rate of the nonlocal field. Gradients are evaluated in the current configuration as

∇• = ∂•/∂ x.
These equations must then be discretized both in time and in space, as detailed in Appendix C. For the spatial discretization, quadratic shape functions are used for displacement degrees of freedom (DOFs), and linear shape functions are used for the nonlocal variables (The notation P2P1P1 is used in this paper to refer to this choice of interpolation order).

The discretization of the set of equations 25 and 26 leads to a global system built thanks to the assembly of elementary reactions. This system relates the unknown variables to the external forces through the assembly of elementary stiffness matrices. The calculation of the different terms of these elementary matrices requires the evaluation of the consistent tangent matrix K mat , as detailed in Appendix D. The resolution of this system is performed implicitly using a Newton-Raphson algorithm and leads to the simultaneous determination of the DOFs associated with the displacements and the nonlocal variables (as opposed to a staggered resolution scheme).

Fitting of the model parameters

Model parameters

The constitutive model of [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF] being relatively close to the local one presented in this paper, it is possible to use the same elasto-viscoplastic parameters and most of the GTN parameters (see table 1). The low initial porosity (f 0 ) corresponds to that of a modern line pipe steel. It corresponds to the MnS inclusion volume fraction as these particles easily debond from the matrix.

Table 1: Set of material parameters from [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF] Young modulus

E 210 GPa Poisson ratio ν 0.3 Isotropic hardening (power law) K 795 MPa e 0 0.002 n 0.13 Gurson criterion q 1 1.5 q 2 1. f 0 1.5 10 -4 Viscosity ṗ0 1s -1 σ 0 55 MPa n 5
The main difference between the constitutive model of [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF] and the local model in this paper is the damage nucleation rate function A n (κ) in equation 6. In the model described in [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF][START_REF] Besson | Modeling of plane strain ductile rupture[END_REF], A n (κ) is a function of the growth porosity f g , whereas in this paper, this function is chosen to be a constant (A n ) when the accumulated plastic strain reaches a threshold κ c :

A n (κ) = A n if κ > κ c 0 otherwise (27)
The model was changed as the initial version did not allow forming cupcone fracture using the nonlocal formulation. This leaves thus two material parameters, A n and κ c , to be fitted. Besides, the nonlocal model requires the identification of the two nonlocal characteristic lengths l ω and l κ .

Tests used for the fit of the nucleation law

The experimental database consists of an axisymmetric tensile test presented in [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] and a plane strain test described in [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]. Both tests were carried out on the same material, i.e., X70 HSLA (high strength low alloyed) ferritic-pearlitic steel, and are experimental results from [START_REF] Rivalin | Développement d'aciers pour gazoducs à haute limite d'élasticité et ténacité élevée : Mécanique et mécanismes de la rupture ductile à grande vitesse[END_REF]. They have been selected because they result in complex crack paths, with cup-cone fracture for the axisymmetric specimen and slant fracture for the plane strain specimen. In addition, the tests have been used in the literature [START_REF] Scheider | Simulation of cup-cone fracture using the cohesive model[END_REF][START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF] as a reference to identify and study ductile fracture models.

Finite element simulations are carried out for both tests in order to fit the parameters (trial and error procedure). Meshes, dimensions as well as boundary conditions are shown in fig. 1a and fig. 1b. Mesh design is similar in both cases. In the case of the plane strain specimen, the entire specimen is meshed to enable the formation of a slanted crack path. In both cases, elements have quadratic shape functions with eight nodes and reduced integration (4 integration points). In the following S 0 represents the initial minimal cross-section of each specimen; S 0 = π /4φ 2 0 for the axisymmetric specimen (with φ 0 the initial diameter at the center of the specimen); S 0 = e 0 w 0 for the plane strain specimen where e 0 represents the specimen thickness and w 0 the width of the specimen. The normalized force (F/S 0 ) is then plotted as a function of the minimum diameter variation (∆Φ/Φ 0 ) or thickness variation (∆e/e 0 ).

Parameter fit for the local model

In the case of local models, it is important to keep the same mesh size in areas where cracks propagate [START_REF] Liu | Mesh-dependence and stress singularity in finite element analysis of creep crack growth by continuum damage mechanics approach[END_REF][START_REF] Skallerud | Effects of finite element mesh on the numerical prediction of ductile tearing[END_REF]. Because the diameter of the tensile bar (10 mm) is twice the thickness of the plane stress specimen (5 mm) and because symmetry is not accounted for in the latter case, the same number of elements (N h ) is used to discretize the minimum cross-sections in order to keep the same mesh size. The initial aspect ratio of these elements (r h ) is set to 6:1 so that it leads to approximately square elements at the onset of fracture. Using N h = 60, the initial mesh size is consequently: 83µm × 14µm.

The fitted material parameters are:

A n = 0.2 and κ c = 1.2 (28)
Indeed, using these values, the global responses for both test cases show a good agreement with both reference results, as can be seen in fig. 2. In the case of the axisymmetric specimen (figure 2a), it can be noted that the fitted local model response is identical to the reference one before crack initiation, which corresponds to the sharp load drop. In this regime, the damage has minimal effect on the overall behavior. Crack initiation occurs at a slightly higher diameter reduction (-∆φ/φ 0 = 0.453) compared to the reference (-∆φ/φ 0 = 0.44).

Besides, as load drop is faster for the new fitted parameters, the total diameter reduction (-∆φ/φ 0 ) c = 0.467) is ultimately slightly smaller than in the case of the reference simulation (-∆φ/φ 0 ) c = 0.48. The difference is less than 3 % and is acceptable. As can be seen on the image of total porosity at the end of the simulation (fig. 2a), a cup-cone crack path is obtained.

In the case of the plane strain specimen (fig. 2b), the computed forces before crack initiation are slightly higher than the reference ones. This can be explained by the fact that in the reference simulation [START_REF] Besson | Modeling of plane strain ductile rupture[END_REF], a Hill anisotropic model was used to describe plasticity. Such behavior is not accounted for in the present material model. In this case too, the thickness reduction (-∆e/e 0 ) c = 0.59 predicted using the new model parameters is slightly different from the reference (-∆e/e 0 ) c = 0.566. As in the case of the tensile test, the difference is acceptable. As can be seen on the image of total porosity at the end of the simulation (fig. 2b), a V-shaped crack path is obtained. Note that both V-shaped and S-shaped crack paths can be experimentally obtained [START_REF] Besson | Modeling flat to slant fracture transition using the computational cell methodology[END_REF].

Parameter fit for the nonlocal model

Using the local model, the band width is about the element height in the deformed configuration. For the above-considered cases, this size is l local b ≈ . It is therefore assumed that the mesh size for the local model allows representing the material's internal length following the early work by [START_REF] Rousselier | Ductile fracture models and their potential in local approach of fracture[END_REF]. Having l nl b ≈ l local b is a simple way to establish a correspondence between both models. Is is shown below (section 5.3) that the band width is related to the nonlocal characteristic length by the following relation:

l nl b ≈ 1.5l c .
Using this relation, one gets l ω = l κ ≈ 40µm. Using these values, 90 elements are required to discretize the specimen radius or thickness with at least three elements in the band width. Using fewer elements results in simulations being mesh dependent. This results in an initial element height equal to 10µm. As for the damage nucleation rate parameter A n , since the nonlocal model will always lead to higher ductilities than the local model for a given set of material parameters, it was necessary to impose a larger value: A n = 0.4. Nucleation parameters for the local and nonlocal models are compared in table 2. Note that the number of elements used to mesh the minimum cross-section differs using the above hypotheses. Using N h = 60 for the nonlocal model would result in slightly mesh-dependent simulations.

This fit leads to a relatively good agreement of the global responses for both test cases with the reference results as well as results from the literature [START_REF] Huespe | A finite strain, finite band method for modeling ductile fracture[END_REF][START_REF] Leclerc | A micromechanics-based non-local damage to crack transition framework for porous elastoplastic solids[END_REF], as can be seen in figures 3 and 4. For the simulation of the smooth axisymmetric bar (see figure 3), the engineering stressdiameter reduction curve is particularly close to the reference from [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF], and the total failure diameter reduction is well predicted. As can be seen in the images of total porosity taken at different stages of the simulation in the same figure, a cup-cone crack path is also obtained. In contrast, the global response for the plane strain simulation (presented in figure 4) shows a slight overestimation of the strain at crack initiation compared to the reference and a higher total failure strain. This difference can be explained by the fact that the crack path remains flat (as can be seen on the images of the total porosity maps taken at different stages of propagation). This result is not in agreement with experimental data as well as the simulation using the local framework, which both exhibit slant fracture. The effect of the characteristic lengths on the crack path is discussed below in section 6 where it is shown that slant fracture may indeed be obtained using the nonlocal framework and a proper choice for the characteristic lengths.

Study of the convergence properties of the nonlocal model

In this section, the convergence of the solution with respect to mesh size and mesh orientation is first checked. This allows to determine the mesh size needed to obtain convergence as a function of the internal lengths, which are assumed to be equal. Using the appropriate mesh size, it is then possible to establish a relation between the internal lengths and the numerically obtained band widths.

Effect of mesh size

Several meshes with the same element initial aspect ratio r h = 6 and the same orientation, but different element sizes are used to study the effect of mesh size.

Local model. In the case of the local model, the considered set of number of elements in the width of the central section N h is: 45, 60, 75 or 90. The engineering stress (F/S 0 ) is plotted as a function of the variation of the minimum diameter -∆φ/φ 0 for different mesh sizes in figure 5. On the same graph, the distribution of total porosity f at total failure is shown for the different meshes. Let us note that all the simulations predict crack initiation for the same diameter reduction and that the diameters at full failure (F/S 0 = 0) only slightly differ. In this particular case, the well-known mesh size dependency (see e.g. [START_REF] Liu | Mesh-dependence and stress singularity in finite element analysis of creep crack growth by continuum damage mechanics approach[END_REF][START_REF] Rousselier | Non-Linear Fracture Mechanics[END_REF]) has no significant effect on the global response as the load drop is very sharp. Indeed, whatever the mesh size, the highly damaged area is localized within one row of elements. The similarity of the global responses does not enable to observe a dependence on the mesh size. The total porosity for each mesh size at total failure shows a cup-cone failure with a slight deviation near the axis of symmetry.

Nonlocal model. In the case of the nonlocal model, the considered set of the number of elements in the width of the central section N h is: 20, 45, 60, 75, 90, or 180. In figure 6, the engineering stress-diameter reduction curves and the total porosity f at total failure are shown for the different meshes. The global response appears to be converged, and the simulation leads to a cup cone crack path as soon as N h ≥ 45. Let us note that, for the coarsest mesh (N h = 20), the global response is very different than for the other meshes close to full failure. In that case, there is no crack bifurcation so that the minimum diameter always plastically deforms. In the case of cup-cone failure, the material undergoes elastic unloading at the notch root due to crack bifurcation at some stage of loading so that the minimum diameter no longer changes. Obviously, for N h = 20, the mesh is too coarse to properly capture the width of the highly damaged zone, and the simulation is, in that case, mesh size-dependent. In all cases where cup-cone fracture is obtained, there exist two possible symmetric crack paths when flat propagation ends. One path is selected, but the other path can still be visualized. In the case of the smallest mesh size (N h = 180) a zigzagging crack path is obtained. In that case, the crack first runs downwards (-45 • ) and then upwards (+45 • ). Such a behavior is often experimentally observed [START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF]. The global responses are similar for the mesh sizes N h > 20. The total porosity for each mesh size at total failure shows a cup-cone failure with a straight first horizontal part except for the coarsest mesh size N h = 20, which exhibits a flat crack path. The finest mesh size N h = 180 enables the crack path to bifurcate twice, which can also be experimentally observed. These results indicate mesh independence.

Effect of mesh orientation

In this section, three meshes having the same number of elements in the central part and the same initial ratio r h = 6 : 1 are used. However, the mesh is slightly tilted with respect to the horizontal axis. The tilt angles are: 0 • , 5 • and 10 • . Note that because of large strain, the final tilt angle is much larger than the initial one.

Local model. In the case of the local model, the number of elements in the central part is N h = 60. In figure 7, the global responses and the total porosity f at total failure are shown for the different meshes. Although the global responses are very similar, it is interesting to note that each mesh leads to a different crack path. Indeed, when there is no tilt, a classical cup-cone fracture path is obtained, as shown above. For a tilt angle of 5 • , the crack first runs upwards from the center of the specimen and then runs downwards after reaching the mid-radius. For a tilt angle of 10 • , the crack runs upwards following one row of elements. In both cases, cup-cone fracture is not observed, and the predicted crack path is only formed by cones. The local model thus clearly suffers from a strong dependence on mesh orientation. This implies that the crack path can hardly be predicted using local models. They can be useful if the predicted crack path corresponds to the expected one. In that case, the local stress and strain history can be investigated to analyze/understand crack formation. This can be done by prescribing the crack path using the computational cell methodology [START_REF] Xia | Ductile crack growth -I. A numerical study using computational cells with microstructurally-based length scales[END_REF][START_REF] Besson | Modeling flat to slant fracture transition using the computational cell methodology[END_REF]. The result also indicates that automatic meshing can hardly be used to generate the mesh because its design cannot be fully controlled. Simulations with automatic remeshing during the calculation are also not possible.

Nonlocal model. In the case of the nonlocal model, the number of elements in the central part is N h = 90. This number of elements was chosen so as to approximately obtain the same band width as in the case of the local model (N h = 60) following the discussion in section 3.4. Note, however, that results obtained with the local model for N h = 90 also lead to the disappearance of cup-cone fracture for a tilt angle equal to 5 and 10 • . In figure 8, the engineering stressdiameter reduction curves and the total porosity f at total failure are shown for the different meshes. All the simulations give the same global response and lead to a similar cup-cone crack path. In all cases, the band width of the damage localization band is the same. Therefore, it can be concluded that results are converged and insensitive to mesh orientation. This paves the way for the use of automatic meshing techniques of complex parts and remeshing.

Analysis of the damage localization bands

In this section, the relation between the characteristic length and the width of the localization band is investigated. To do so, it is first necessary to define how to measure the width of a localization band. The finite element discretization needed to obtain a converged band width is then investigated. Based on these results, the relation between the characteristic length and the band width is studied for both axisymmetric and plane strain specimens, as this relationship may depend on the stress state.

Definition of the band width

In the literature, there is no agreement, to the authors' knowledge, on a definition of a measure of the width of the localization bands. [START_REF] Maziere | Numerical aspects in the finite element simulation of the Portevin-Le Chatelier effect[END_REF] study localized bands in a nickel-based superalloy which exhibits a Portevin-Le Chatelier effect. Their proposed method is based on the analysis of the local plastic strain rate. It is applied to FE simulations but is designed to mimic actual experimental techniques (see references in the paper). In [START_REF] Labergere | Strain rate distribution and localization band width evolution during tensile test[END_REF], the authors study a copper sheet and also use the strain rate profile to determine the band width. The profile is fitted by a pseudo-Voigt function for which the band width is a fitting parameter. [START_REF] Gorodetskyi | Detecting precursors of localization by strain-field analysis[END_REF] propose a method able to detect discontinuities such as displacement jumps which can be applied to both FE simulations and experimental DIC results. In [START_REF] Ait-Amokhtar | Kinematics fields and spatial activity of Portevin-Le Chatelier bands using the digital image correlation method[END_REF], plastic strain localization is experimentally studied as a function of the applied strain rate in the case of an Al-Mg alloy exhibiting a PLC effect. The bandwidth is defined "as the width at the middle height of the localized strain band". In this study, it is proposed to define the band width (referred to as l v b for variable v = ω or v = κ) as the width over which the considered variable reaches half its maximum value. In that sense, this definition corresponds to that of [START_REF] Ait-Amokhtar | Kinematics fields and spatial activity of Portevin-Le Chatelier bands using the digital image correlation method[END_REF]. In the following, the case of the axisymmetric specimen will only be described as similar results are obtained for the plane strain specimen. Profiles in the deformed configuration along the symmetry axis for both ω and ω are exemplified in fig. 9 where the band width for ω is shown by the blue arrows and the band width for ω by the red arrows. The profiles for ω can easily be plotted as this variable is defined at nodes and can be interpolated. Plotting the profiles for the local variable ω first requires extrapolating values at Gauss points to nodes (red curve in figure 9). In all cases, the band width for the local variable is about the size of one single element in the deformed configuration. In the following, the width of the localization bands is determined when the maximum of ω reaches 0.5 along the profile. In that case, the material points are close to full failure. Applying this definition to measure the band width for the local and nonlocal volume variations along the vertical symmetry axis for the axisymmetric specimen leads to consistent relative measures. Indeed, as shown in fig. 9, the band width measured on ω is larger than the band width measured on its local counterpart ω, even if the maximum of the nonlocal variable is lower than the maximum of the local variable, as expected using an implicit gradient nonlocal formulation.

Conditions to obtain a converged band width

In this part, the evolution of the band width with element size is studied. As the present model makes use of gradients of ω and κ computed in the current configuration, it is important to evaluate the localization band width in the deformed configuration as proposed above. The element height after deformation and when fracture starts is referred to as H r1 in the following (see fig. 10). The element height along the element row just above is referred to as H r2 . Due to different loading histories, the element elongation at fracture can differ greatly from one simulation to another. It is therefore important that the initial element height (h min ) is small enough to assure that the element height at fracture (H r1 ) is still able to represent the localization band. The convergence of the band width should therefore be checked with respect to H r1 . Different values for the number of element along the thickness of the specimens are used: N h = 540, 270, 180, 90, 60, 45, 30, 20. This study is performed for two values of the characteristic length, taken here equal for both nonlocal variables (l κ = l ω = 80µm and l κ = l ω = 40µm). Both axisymmetric and plane strain specimens are considered. The considered variable is the nonlocal volume variation but similar results can be obtained if the nonlocal effective accumulated plastic strain is used to measure the localization band width.

Figure 11 shows the evolution of the band width size as a function of the element height in the current configuration (H r1 ), as defined in fig. 10. Each set of points can be approximated by a line, which means that it can be assumed that the band width is an affine function of the mesh size. Using this approximation, the band width for H r1 → 0 can be extrapolated to a non-zero value which is an estimated value of the band width l ∞ b free from discretization error. In practical applications, this error is always present.

Figure 11: Evolution of the measured band width for both l κ = l ω = 40µm and l κ = l ω = 80µm on both the axisymmetric and plane strain tensile tests as a function of the element height in the current configuration H r1 . There appears to be a linear relation between the band width and the element height. This means that it is possible to estimate a value of the band width l ∞ b free from any discretization error for H r1 = 0.

Relation between the band width and the internal lengths l κ and l ω

Let us note that in fig. 11, for a given geometry, the line corresponding to the largest characteristic length (i.e. l κ = l ω = 80µm) has approximately the same slope as the line corresponding to a smallest characteristic length (i.e. l κ = l ω = 40µm) but leads to a higher band width estimate l ∞ b . In order to validate this observation for other characteristic lengths and to establish a simple relation between the characteristic length and the band width l ∞ b , the same series of simulations was performed for additional values for of l ω = l κ : i.e. 20µm, 60µm, and 100µm. Plotting the evolution of the measured band width l ∞ b versus the imposed characteristic length (see fig. 12), it can be observed that for a given geometry, the measured points approximately lie on a line. Using this regression, it is observed, as expected, that a null band width is obtained when the characteristic length is also null. Let us note that the slope in this linear relationship is different for each geometry. However, the difference between the relation obtained for plane strain and axisymmetric elements remains small: l ∞ b ≈ 1.54l κ,ω in the first case and l ∞ b ≈ 1.69l κ,ω in the second case. This simple relation presents a major advantage as it enables to estimate, for a given geometry, and based on the choice of the characteristic length, the width of the localization band. It is then possible to choose the number and size of the elements in the localization band (respectively denoted a and H r1 ) in order to properly capture its gradients with: a × H r1 < l ∞ b . In all the nonlocal simulations of this paper (except for fig. 6), the factor a is at least equal to 3, i.e. there are at least three P2P1P1 elements in the band width at the onset of failure.

Figure 12: Evolution of the estimated band width free from discretization error for both the axisymmetric and plane strain tensile tests as a function of the characteristic length (l κ = l ω ). There appears to be a proportional relation between the band width and the characteristic length, with slightly different coefficients for each geometry. This means that it is possible to estimate the size the band width for a given characteristic length and a given geometry.

Using two characteristic lengths to model the cup-cone and slant fracture

In this section, two different characteristic lengths (l κ ≤ l ω ) are used to model the occurrence of cup-cone and slant fracture.

Procedure

Several simulations were performed for different characteristic lengths varying between 20µm and 100µm and different ratios l ω /l κ (i.e. 1, 2, 3, 4 and 5). Let us note that larger ratios were not considered in order to avoid too large computational costs due to mesh refinement associated with the requirement to have 3 elements in the band width at the onset of fracture (3

× H r1 < l ∞ b ).
The relative values of the two characteristic lengths are chosen to be consistent with the physical degradation of metallic materials, and especially steels. Indeed, in steels, two damage mechanisms exist: (i) void growth from sulfides/oxides and (ii) damage nucleation at iron carbides Fe 3 C, which occurs at high strains [START_REF] Tanguy | Plastic and damage behavior of a high strength X100 pipeline steel: experiments and modelling[END_REF]. As the spacing between sulfides/oxides is larger than between carbides, it makes sense to only consider l ω ≥ l κ .

Figure 13: Contours of total porosity at total fracture for different characteristic lengths l κ and l ω for the axisymmetric specimen. Cup-cone crack paths are observed only for l κ < 80µm. The characteristic length l κ should then be relatively small compared to the specimen size to observe a cup-cone fracture. The larger the characteristic lengths, the wider the localization band (which is linked to the displacement at failure) and the further from the symmetry axis the bifurcation..

Cup-cone fracture

For the axisymmetric specimen, the contours of total porosity at full fracture for different characteristic lengths are displayed in fig. 13. It appears that large values for l κ prevent the occurrence of crack bifurcation for the considered specimen. As seen in fig. 13, simulations with l κ ≥ 80µm predict a flat crack while all simulations with l κ < 80µm predict the cup-cone crack path. For a given value of l κ , it can be noticed that the larger l ω , the wider the localization band. Similarly, for a given value of l ω , the band width increases with increasing l κ . This is indeed expected since a large characteristic length leads to a larger band. In addition larger bands tend to increase the length of the flat part of the crack path. The flat central crack then bifurcates to form the slanted crack path at +45 • or -45 • . In one case (l ω = 60µm and l κ = 20µm) a zigzagging crack path is obtained, but it was impossible to determine why such a crack path was formed, although this phenomenon is often experimentally observed.

Slant fracture

For the plane strain specimen, the contours of total porosity at total fracture for different characteristic lengths are displayed in fig. 14. For l κ ≥ 40µm flat fracture is always obtained so that only one crack path is shown (l ω = l κ = 40µm). For l κ = 20µm slant fracture is always obtained (V-shape). In that case, the effect of l ω is small. The band width only slightly increases with increasing l ω .

Figure 14: Contours of total porosity at total fracture for different characteristic lengths l κ and l ω for the plane strain case. Slant fracture is observed only for l κ = 20µm. The effect of l ω appears to be limited in this case.

Size effect

It is interesting to revisit the above observations from a more physical point of view and to consider, for a given set of material lengths, how specimen size impacts the formation of the cup-cone and slant fracture. Indeed increasing the material lengths for a given specimen size is equivalent to decreasing the size of the specimen for given material lengths. Similar crack paths are obtained as long as the ratios l ω /L and l κ /L are constant where L is a characteristic dimension of the simulated structure (diameter or thickness in the present case).

From the results corresponding to the axisymmetric specimen (fig. 13), it can be concluded that it is easier to form a cup-cone crack path for large specimens (case corresponding to l κ < 80µm for a given geometry) than for smaller specimen sizes. From the results corresponding to the plane strain specimen (fig. 14), it can be concluded that only sufficiently large plane strain specimens break following the slant crack path. This opens the way to the determination of material lengths by using homothetic specimens so as to obtain different crack paths. Lengths could be fitted to represent the different fracture patterns. This will obviously require the use of very small specimens in the case of metals where lengths are expected to be of the order of the mean spacing between damage initiation sites.

Conclusions

In this work, a nonlocal GTN model for ductile fracture was proposed which uses two material lengths. The model is based on an implicit gradient formulation applied to the plastic volume variation (material length l ω ) and the accumulated plastic strain (material length l κ ). The model allows to regularize void growth and strain-controlled nucleation. The model parameters are fitted to reproduce the global response for tests carried out on a pipe line steel using a tensile bar and a plane strain specimen. These experiments have already been used in the literature to test several models and numerical techniques. Each test is representative of a characteristic crack path: cup-cone fracture (tensile test) and slant fracture (plane strain test).

The material parameters were first fitted assuming that both lengths are equal (l ω = l κ ) in order to reproduce the macroscopic behavior of both specimens. It was checked that mesh independence is obtained if a sufficiently fine mesh is used. In particular, it was shown that a local model is strongly dependent on mesh design (mesh orientation in the present study) whereas the nonlocal model is not. This implies that remeshing cannot be applied if a local model is used but with a nonlocal model, it can be a solution to reduce computational cost in cases where long crack propagation needs to be simulated.

A novel methodology was proposed to estimate, for a given geometry, the band width associated with a given characteristic length free from discretization error. By applying this methodology for different characteristic lengths, it was possible to establish a linear relation between the estimated band width and the characteristic length. This piece of information was then used to select the proper discretization of the localization bands for given values of the characteristic lengths in order to obtain converged solutions. In practice, the element size at failure is recommended to be three times smaller than the characteristic lengths.

Finally, two distinct characteristic lengths were used to simulate both test cases. It was first concluded that the characteristic length controlling strain controlled damage nucleation plays an important role on the occurrence of crack bifurcation: it must be small enough relative to the specimen size in order to observe a cup-cone crack path or a slanted crack path. It was also noted that increasing the characteristic lengths leads to an increase of both the band width and the radius at which the cup-cone crack path starts to bifurcate and propagate either at +45 • or -45 • . Larger values favor flat crack advance. This study also enabled to conclude that crack path should be affected by the size of the specimens. Pure flat fracture should be favored when using (very) small specimens. It could theoretically be possible to identify them separately. This paves the way for the more direct determination of material lengths by using homothetic specimens to obtain different crack paths. The fitted lengths should then be able to represent both the macroscopic behavior and the crack path for all sample sizes.

Appendix A. Implementation of the nonlocal GTN model

Integrating the set of equations 16 to 20 over a finite time step ∆t using a fully implicit scheme is equivalent to solving the following set of non linear equations with respect to the increments of the state variables ∆V S = (∆ε e , ∆κ, ∆f g , ∆f n , ∆ω) for a given increment of the input variables and .24) and

∆V IN = • Derivatives of R g : ∂R g ∂∆f g = 1 + ∆ω (A.15) ∂R g ∂∆f n = ∆ω (A.16) • Derivatives of R n : ∂R n ∂∆f n = 1 (A.17) • Derivatives of R ω : ∂R ω ∂∆ε e = -(1 -f )∆κI : N : E (A.18) ∂R ω ∂∆κ = -(1 -f )trace(n) (A.19) ∂R ω ∂∆f g = -(1 -f )∆κ trace(n f ) (A.
N = ∂n ∂σ = -h ∂ 2 G ∂σ 2 -h 3 ∂ 2 G ∂σ 2 ν ⊗ ν + h 2 ∂ 2 G ∂σ∂σ ⊗ ν + ν ⊗ ∂ 2 G ∂σ∂σ (A
n f = ∂n ∂f = -h ∂ 2 G ∂σ∂f + h 2 ∂G ∂f ∂ 2 G ∂σ∂σ + h 2 ∂ 2 G ∂σ ∂f -h ∂ 2 G ∂σ 2 ∂G ∂f ν (A.25)

Appendix B. Consistent tangent matrix

Once the solution of the system of equations A.1 to A.5 is found, it becomes possible to numerically compute the consistent tangent matrix. For any given small variation of the input variables δV IN , the state variables will be modified so that R remains null. One therefore has: The strain ∆ε is updated as:

δT = δ0 = ∂R
∆ε = Q T 1 2 .∆D.Q 1 2 (C.6)
The constitutive equations are then integrated (see Appendix A) using this strain increment. The resulting stress tensor is rotated back as Σ = Q 1 so that: The term K κu is the most complex of the three terms to calculate. Indeed it requires to calculate the derivative of Ḟ κ with respect to the displacements, which involves three terms: -J κN , Jl 2 l ( ĠT .G + G T . Ġ).κ e and J (κ -κ)N + l 2 κ G T .G.κ e def = JT k . Considering the first term, one has:

K uκ = Ω e
-J κu N = -J ∂∆κ ∂∆ε .B D . ue
To compute the second term, one first notices that G = F -T .G 0 (note that in this case, the Voigt notation is not used). Therefore

  (a) Axisymmetric mesh with (cax8r) element (b) Plane strain specimen mesh with (c2d8r) element

Figure 1 :

 1 Figure 1: Dimensions and boundary conditions for the tensile test simulations of the axisymmetric specimen 1a and the plane strain specimen 1b.

  (a) Axisymmetric specimen (b) Plane strain specimen

Figure 2 :

 2 Figure 2: Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for the axisymmetric specimen (left) and the thickness reduction -∆e/e 0 for the plane strain specimen (right). The comparison with global curves from the literature (respectively from[START_REF] Besson | Modeling of crack growth in round bars and plane strain specimens[END_REF] and[START_REF] Besson | Modeling of plane strain ductile rupture[END_REF]) indicates a correct parameter fit. Total porosity is displayed at total failure for both simulations: a cup-cone failure (left) and a slant failure (right) can be respectively observed.

Figure 3 :

 3 Figure 3: Evolution, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for the axisymmetric specimen with a central spatial discretization of N h = 90. The confrontation with global curves from the litterature (Besson et al., 2001; Huespe et al., 2012; Leclerc et al., 2020) and the local model indicates a correct parameter fit. Total porosity is displayed at various times of the simulation to illustrate the formation of the expected cup-cone crack path.

Figure 4 :

 4 Figure 4: Evolution, for the nonlocal model, of the engineering stress as a function of the thickness reduction -∆e/e 0 for the plane strain specimen with a central spatial discretization of N h = 90. The confrontation with global curves from the literature(Besson et al., 2003; Huespe et al., 2012; Leclerc et al., 2020) and the local model indicates a correct parameter fit. Total porosity is displayed at various steps of the simulation to illustrate the propagation of a flat crack instead of a slanted one. This illustrates the need to properly calibrate the nonlocal characteristic lengths.

Figure 5 :

 5 Figure 5: Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for different spatial discretization in the case of the axisymmetric specimen.The similarity of the global responses does not enable to observe a dependence on the mesh size. The total porosity for each mesh size at total failure shows a cup-cone failure with a slight deviation near the axis of symmetry.

Figure 6 :

 6 Figure 6: Evolution, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for different spatial discretization in the case of the axisymmetric specimen.The global responses are similar for the mesh sizes N h > 20. The total porosity for each mesh size at total failure shows a cup-cone failure with a straight first horizontal part except for the coarsest mesh size N h = 20, which exhibits a flat crack path. The finest mesh size N h = 180 enables the crack path to bifurcate twice, which can also be experimentally observed. These results indicate mesh independence.

Figure 7 :

 7 Figure7: Evolution, for the local model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for various mesh orientations in the case of the axisymmetric specimen. The global responses are similar for all the orientations. However, the total porosity for each mesh size at total failure shows a different crack path for each mesh orientation, which indicates mesh dependence.

Figure 8 :

 8 Figure 8: Evolution, for the nonlocal model, of the engineering stress as a function of the diameter reduction -∆φ/φ 0 for various mesh orientations in the case of the axisymmetric specimen. The global responses are similar for all orientations. Besides, the total porosity for each mesh size at total failure shows an identical cup-cone crack path for each mesh orientation. These observations indicate mesh independence.

Figure 9 :

 9 Figure 9: Band width measure for the nonlocal volume variation ω when ω max = 0.5. Evolution of the local and nonlocal volume variations along the symmetry axis. The double arrows show the measure of the corresponding band widths l ω b and l ω b .

  Figure 10: (a) Definition of H r1 and H r2 . (b) Contours of the nonlocal volume variation ω when ω max = 0.5. Arrows here define current element heights H r1 and H r2 .

  of σ are obtained considering the stationarity of the G function (eq. 3) : G = 0 and δG = 0. One then has:

  matrices ∂R/∂∆V IN and ∂∆V OUT /∆∂V S are computed as blockmatrices. The calculation is straightforward with: the case of κ and ω, the output and state variables are similar. All other terms are null.

  bloc K uω is computed in a similar way by replacing κ by ω from equation D.11 to equation D.14.Calculation of K κu , K κκ , K κωCalculation of the second bloc line can be obtained by writing the variation ofF κ : -κ)N + Jl 2 κ ĠT .G + G T . Ġ .κ e + Jl 2 κ G T .G. κe + J (κ -κ)N + l 2 κ G T .G.κ e dΩ 0 (D.15)and considering the derivative κ as the sum of the partial derivative of κ regarding the variables u e , κ e and ω e : κ = κu + κκ + κω (D.16)

Ġ

  = Ḟ -T .G 0 remembering that Ḟ -T = -L T .F -T Ġ = -L T .F -T .G 0 = -L T .G and ĠT = -G T .LLet V be the operator linking the Voigt representation of a tensor to the matrix representation such thata = V.{a} or a ij = V ijk a k and V * such that a T = V * .{a}. Indeed V * ijk = V jik . Using this notation, G T . Ġ.κ e = -G T .L T .G.κ e = -G T . (V * .(B L . ue )) .G.κ e = -G T . (V * .(B L . ue )) . ∇κ G T ij Ġjn κ e n = -G T ij V * jkl B L lm ue m ∇ k κ = -G T ij V kjl B L lm ue m ∇ k κ = -G T ij ∇ k κV kjl B L lm ue m G T . Ġ.κ e = -G T .( ∇κ.V.B L ).u e Similarly ĠT .G.κ e = -G T .( ∇κ.V * .B L ). ue m and ĠT .G.κ e + G T . Ġ.κ e = -G T .( ∇κ.(V + V * ).B L ). ue = -2G T .( ∇κ.V.B D ). ue Finally, one now considers the last term JT k .JT k = Jtrace(L)T k = J({1}.B D . ue )T κ = JT κ ⊗ ({1}.B D ). ue (D.17)In the end, one finally gets:κ G T .( ∇κ.V) + T κ ⊗ {1} .B D dΩ (D.18)where ∂∆κ ∂∆ε is computed using a sub block matrix of the consistent tangent matrix K mat (eq. B.4) as:Calculation of K ωu , K ωκ , K ωω Calculation of the last bloc line can be obtained exactly in the same way as the second bloc line by writing the variation of F ω . It leads to the same terms only replacing κ by ω.

Table 2 :

 2 Material parameters for the local and nonlocal models.

		local model nonlocal model
	κ c	1.2	1.2
	A n	0.2	0.4
	l κ = l ω	-	40µm
	70µm. The characteristic length for the nonlocal model was then selected so that
	the resulting band width (l nl b ) is approximatively equal to l local b

  The outputs variable are expressed as function of the state variables only so that:

	∂∆V IN	.δV IN +	∂R ∂∆V S	.δV S =	∂R ∂∆V IN	.δV IN + J .δV S (B.1)
	so that		δV S = -J -1 .	∂R ∂∆V IN	.δV IN	(B.2)
	δV OUT =	∂∆V OUT ∂∆V S	.δV S = -	∂∆V OUT ∂∆V S	.J -1 .	∂R ∂∆V IN	.δV IN	(B.3)

so that the consistent tangent matrix is equal to:

Note that the model can be easily extended to plastically anisotropic materials by replacing the von Mises by any stress measure accounting for anisotropy(Benzerga and Besson,

2001;[START_REF] Tanguy | Plastic and damage behavior of a high strength X100 pipeline steel: experiments and modelling[END_REF][START_REF] Shinohara | Anisotropic ductile failure of a highstrength line pipe steel[END_REF].

 (∆ε, ∆ω, ∆κ):

Using a fully implicit scheme, all variables in the previous system are evaluated at the end of the time increment. Solving the system requires the evaluation of its Jabobian matrix which is formally expressed as:

The Jacobian matrix can be computed blockwise as follows (the zero terms are omitted).

• Derivatives of R e :

)

Elements have nodal DOFs corresponding to the displacements and nonlocal variables κ and ω. They are represented as vector u e , κ e and ω e . Linear shape functions are used to interpolate the nonlocal nodal variables in the elements (i.e. κ e and ω e ), so that:

where N is a matrix formed with the linear shape functions. The gradients of κ and ω are computed as:

where the matrix G is formed with the derivatives of the linear shape functions with respect to the final configuration. The gradient with respect to the initial configuration is computed as ∇ 0 κ = G 0 .κ e where G 0 is formed with the derivatives of the linear shape functions with respect to the initial configuration.

Standard quadratic shape functions are used to interpolate the nodal variables associated with the displacements in the elements (i.e. u e ). The transformation tensor F and the velocity gradient tensor L are computed as:

where B F is a matrix formed with the derivatives of the quadratic shape functions with respect to the initial configuration whereas B L uses the derivatives of the same shape functions with respect to the final configuration. The deformation rate D is expressed as B D . ue .

Appendix C.2. Temporal discretization

The finite strain formulation is obtained using a mid-point integration scheme over the time step [t, t+∆t]. The transformation tensor at t+ 1 2 ∆t is first computed as

2 ∆F using the standard polar decomposition. The tranformation increment is then computed as:

where ∆D and ∆W are respectively the symmetric and skew-symmetric parts of ∆L. The rotation Q is updated as (based on the mid-point integration of eq. 22)

Elementary reactions

From the discretized form of the weak formulation (25 and 26), elementary reactions associated unknowns are given as:

where notation {Σ} indicates that the tensor Σ is expressed using Voigt notations.

Let us underline that the integrals are taken over the current configuration. Integration is performed using standard Gauss integration. Let us introduce here the function M R such that the product c of two second order tensors a and b, usually written c = a.b, can be expressed using Voigt notations as:

The matrix M R (b) thus depends linearly on b.

Appendix D.2. Elementary stiffness matrix

It is then necessary to evaluate the elementary stiffness matrix which is computed as a block matrix.

Calculation of K uu , K uκ , K uω Calculation of the first bloc line of the elementary stiffness matrix can be obtained through the derivative of F u with respect to u e , κ e and ω e :

Derivation of Ḟ u leads to:

In the central term, the derivative tensor Σ can be separated into three parts corresponding to variations relative to the different unknowns:

The calculation of K uu follows usual derivations for finite strain formulation and leads to:

where K M uu is the part related to the material non-linearity and K G uu is the geometrical non-linearity.

Then, to compute the K uκ term, let us focus on the part of the derivative of Σ that depends on κ:

The corresponding variation of F u is:

∂Σ ∂κ (N . κe ) dΩ 0 (D.12)


H Frånberg-Delahaye

Ganil

THE NEW GANIL BEAMS: COMMISSIONING OF SPIRAL 2 ACCELERATOR AND RESENT DEVELOPMENTS

The GANIL installation at Caen in France has been operating with warm temperatures cyclotrons for heavy ion beam physics since 1983. The accelerated stables beams widely ranges from Carbon to Uranium beams. Low energy and post accelerated radioactive ion beams are also being provided.

The GANIL laboratory has newly increased their different ion beams and energies available with the installation and commissioning of a superconducting linear accelerator -SPIRAL2 and its experimental areas. The construction of SPIRAL2 started in 2011, the first beam was extracted at low energy in late 2014 with pre-acceleration in 2017 and since 2021 the new installation delivers beam for nuclear physics experiments. This paper will cover the commissioning and power ramp up of the SPIRAL2 installation at GANIL with its superconducting LINAC -but also the latest development of stable and radioactive ion beams at the cyclotron facility of GANIL.

GANIL COMPLEX

GANIL will in November 2022 celebrate the 40 years from the first extracted beam from the two separated sector cyclotrons (SSC). The first accelerated and extracted beam was a 40 Ar 16+ beam, accelerated to 44 MeV/A. The same beam was used for an experiment 2 months later in January 1983.

Figure 1: The Running GANIL facilities in 2022. On the right is the cyclotron facility and on the left is the SC-LINAC facility.

In the 2015 version of HIAT conference O. Kamalou presented an article of the GANIL Operation Status and New Range of Post-Accelerated Exotic Beams [START_REF] Kamalou | GANIL Operation Status and New Range of Post-Accelerated Exotic Beams[END_REF]. In the same conference J. M. Lagniel presented the Advances of the SPIRAL2 PROJECT [START_REF] Lagniel | Advances of the Spiral 2 Project[END_REF] which had just commissioned the ion sources and the low energy beam transfer lines (LEBT), the SC-LINAC was still under installation.

Today GANIL is running two separated facilities, not yet connected as can be seen in Fig. 1. On the right hand side of the Fig. 1 is the cyclotron complex. The ion beams are produced on the far right hand side before accelerated in one or two steps in the two SS cyclotrons up to 95 MeV/A. On the left in the cyclotron buildings are the RIB factory SPIRAL1. The experimental halls, were the users from different fields explore the stable and radioactive beams, are on the upper part of the figure. On the left hand side in Fig. 1 is the SPIRAL2 installation with a more detailed description in Fig. 2. The four injectors used at GANIL are opposite to Figure 2: A schematic view over the SPIRAL2 installation, the two injectors on the right hand side. each other in the two facilities. In the near future the two facilities will be connected and beam from either of the two facilities will be used in a new experimental area as seen in the section "Next steps for GANIL".

From here we will refer to the cyclotron facility and the SC-LINAC facility (Fig. 3). In total GANIL is running seven ECR ion sources, one FEBIAD ion source, one ECR charge breeder, four injectors, five cyclotrons, one RFQ and one superconducting LINAC, delivering stable ion beams (SIB) and radioactive Ion beams (RIB) to nine experimental areas. The facility is operated 24h 7 days a week for 8 to 9 months per year. Since 2019 the cyclotrons and the SC-LINAC facilities are sharing the operation time, meaning that at GANIL there is a cyclotron season and a super conducting LINAC season. This arrangement divide individually the two facilities into separate 4-6 months uptime periods with 6-8 months maintenance and upgrading periods each year. This gives an opportunity to work on upgrades and regular improvements of ion beam transport and new beams while providing beams 8-9 months a year.

While the Cyclotron facility celebrate 40 years of operation the SC-LINAC is still in the first years of regular operation. In here we will present the resent upgrades at the cyclotron facility presented in details in several articles [3] and a short resume of the SC-LINAC commissioning as requested [4].

SC-LINAC -SPIRAL2 PROJECT

The SPIRAL2 project, launched in 2005 was designed to produce high-power proton, deuteron and heavy-ion beams to be accelerated in a superconducting linear accelerator. The beams were to be delivered in three different experimental areas for nuclear physics research. Two of the experimental areas are constructed and the third one was put on hold in 2013 [START_REF] Dolegieviez | Status of the SPIRAL2 project[END_REF]. The experimental areas and the production of RIBs were a part of the initial project. The idea was to produce RIB in one InFlight (S3) and one ISOL facility (Phase 2) to complement the neutron rich side of the table of isotopes and the production of super heavy ion for fundamental applications.

The construction was limited to the injectors and the highpower superconducting LINAC with two experimental areas called "Neutrons for Science" (NFS) [6] and "Super Separator Spectrometer" (S3) [START_REF]Toward the drip lines and the superheavy island of stability with the Super Separator Spectrometer S 3[END_REF] (see Fig. 3).

The new facility was constructed with two injectors, the proton/deuteron source is a 2.45 GHz ECR source (SILHI like [START_REF] Tuske | Light ion source for proton/deuteron production at CEA Saclay for the Spiral2 project[END_REF] it is designed to produce 5 mA proton and deuteron beams in CW or pulsed modes or pulsed modes.

The heavy-ion source was initially an 18 GHz ECR source PHOENIX V2 using three normal conducting coils and a large permanent magnet. It was designed to produce heavy ions with A/q ≤ 3 with a total extracted beam current up to 15 mA CW at 60 kV. The design allows the installation of a dedicated oven reaching 1600 °C for the production of metallic ion beams [START_REF] Barue | Metallic beam developments for the SPIRAL 2 project[END_REF]. It was later on upgraded to a Phoenix V3 [START_REF] Thuillier | First Plasma of the PHOENIX V3 ECR Ion Source[END_REF] also operating at 18 GHz. This modification was performed to be able to provide the intensities of the beams needed for the S3 experiments.

The installation and the commissioning of the accelerator complex was programmed in phases taking into account the availability of the personnel while continuing the operation of the cyclotron facility. The start-up and tests of all equipment have been performed with staff from GANIL and the external collaborating facilities for SPIRAL2. The ion sources and their low energy beam transport (LEBT) lines was commissioned while parts of the SC LINAC was still under installation.

The design goals of the installation are to be seen in Table 1.

The first SPIRAL2 beam at GANIL (protons) was produced by the proton/deuteron source on December 19, 2014. Later on, in January 2020, the first beam of 40 Ar 14+ was The building of SPIRAL2, is actually constructed by free standing blocks due to seismic reasons. The ion sources and LEBT are in one block and the RFQ, Medium energy beam transport line (MEBT), SC LINAC and high energy beam transport section (HEBT) in a second block and the two experimental halls in a third and fourth block. This was taken advantages of during the installation and commissioning of the equipment in the separated blocks. The RFQ and the MEBT was done in parallel of the commissioning of the ion sources from 2014. The low level RF tuning operations of the RFQ, including the voltage law bead-pull measurements and the 40 plunger adjustments ended in March 2015. From these primary tests the TOUTATIS simulations indicated that the theoretical transmission of beam should be 99.7%. The resonance frequency of the RFQ is 88.0159 MHz, and the beam is accelerate up to 0.75 MeV/A. The RFQ and the MEBT was commissioned using a specific diagnostic box before the SC LINAC, physically separating it from the LINAC but allowing to calibrate the beam and adjust the RFQ operation completely before the commissioning of the LINAC. The RFQ was commissioned at 114 kV with a pulsed 4.8 mA proton beam with 100% transmission in end of 2015 [START_REF] Giacomo | RF Commissioning of the SPIRAL2 RFQ in CW Mode and Beyond Nominal Field[END_REF].

In the years of 2015-2018 the commissioning the RFQ and the MEBT was performed with beams of H + , 4 He 2+ , 18 O 6+ and 40 Ca 14+ . The diagnostic box installed [START_REF] Jamet | SPIRAL2 Diagnostic Qualifications with RFQ beams[END_REF] before the SC-LINAC confirmed 100% transmission [START_REF] Ferdinand | Final Results of the SPIRAL2 Injector Commissioning[END_REF] showing no beam losses during this first accelerator stage of the installation (Fig. 4).

The diagnostic plate included temporary diagnostic tools for commissioning of the RFQ. The diagnostics used were ACCT and DCCT [START_REF] Leloir | Measurement Uncertainty Assessments of the Spiral2 ACCT/DCCT[END_REF] for intensity monitoring, wire profilers for the beam size and position, Alison type scanners for emittance measurements in H and V direction, TOF (Time of Flight) monitor for energy measurements and phase coupled with two BPMs [START_REF] Vignet | The Beam Profile Monitor for Spiral2[END_REF] (Beam profile monitor), a BEM (Beam extension monitors) [START_REF] Revenko | Bunch Extension Monitor for LINAC of SPIRAL2 Project[END_REF] for the longitudinal emittance measurements [START_REF] Normand | SPIRAL2 RFQ bunch lengths and longitudinal emittance measurements[END_REF] and FFC (fast faraday cups).

The validation of the RFQ line was done in November 2018 and the D-plate was removed, the full MEBT was In the full MEBT line the three re-bunchers and a beam bunch selector (SBS) [START_REF] Giacomo | Proton beam commissioning of the SPIRAL2 Single Bunch Selector[END_REF] was installed. The SBS allows an operation of 1 bunch over 100 up to 1/10000, needed primarily for experiments requiring a precise timing between the bunches to insure that there is no overlap of the bunches for precise energy measurements.

In the meantime the SC-LINAC installation was ongoing and in 2017 the cryogenic system for the cryo-modules were commissioned with a first complete cool down in November 2017 [START_REF] Ghribi | Cryogenic Tests of the SPIRAL2 LINAC Systems[END_REF]. The RF tests in the LINAC were performed in parallel to the operation of the cyclotrons during the first semester in 2019.

The GANIL cyclotron and SPIRAL2 installations are constructed and operated under French nuclear authority regulations. Before the authorisation of construction a file explaining the full operation and safety measurements were sent in. The different systems for the safety controls were then constructed and tested. The final authorisation to start the facility with acceleration of the beams in the LINAC was followed four months later after having complete testing and documentation of the different safety systems. Only after this the full authorisation was given and the beam could finally be sent to the LINAC.

The SC LINAC is designed with 26 cavities, 12 A-type cryo-modules each comprising one low-β cavity and 7 B -type cryo-modules with two high-β cavities [START_REF] Olry | Spiral2 Cryomodules B Tests Results[END_REF][START_REF] Marchand | Performances of SPIRAL 2 Low and High Beta Cryomodules[END_REF].

The Power Up of the Beams

In an early stage of the project the strategy for the beam commissioning of the LINAC was to go in four phases as shown in Fig. 5. The proton beam, were to be sent to the SAFARI beam dump during this commissioning.

1. Beam transmission at low energy as given by the RFQ. The first beam, a 1 W proton beam, without energy increase and at low intensity was transported through the SC-LINAC tuned in a rebuncher mode in the middle of November 2019. During the following two weeks a beam energy increase up to 33 MeV with a pencil beam of 200 µA and a full 6.4 W beam was sent to the SAFARI beam dump. Before the winter shut down in December, this beam was sent to the NFS experimental area for full identification of the energy and composition of the beam.

From here on, the LINAC operation is performed with the GANIL teams, also working on the operation of the Cyclotrons. The beam commissioning of the LINAC are therefore limited to 4-5 months every year at the moment.

The phase three continued in the second semester of 2020, after the cyclotron run programmed for the first semester of 2020. After restarting the full complex, retuning the beam and do all the initial tests as in the end of 2019, the proton beam intensity was increased to 4.8 mA and the in the following week, the beam power were then increased to 16 kW (Fig. 6) using a 10% duty cycle, with a beam transmission near of 100% into the LINAC and high energy lines. For each step of the power up the beam the transmission was monitored throughout the installation from the MEBT section, all through the LINAC and in the HEBT section [4]. The results was compared to theoretical TRACEWIN calculations for confirmation.

The next step in the commissioning were to proof the capacity of accelerating heavier ions, starting with deuterons. Although, accelerating the deuterons, an eventual problem of activation should be minimized. So, in the last weeks of operation in 2020 a 4He beam were accelerated to simulate the deuteron beam transport.

After the run in 2019 and 2020, 2021 was the start of the physics experiments and real challenges of the accelerator. The different systems (security, RF, LLRF, diagnostics) was pushed to their limits to validate the operation. While in 2020 the technical and scientific goals were to proof the operability and the capacity of the different systems of the installation, showing the capabilities of controlling and maintaining an accelerated beam of high intensity. The 2021 program was to proof different beams for physics could be produced, accelerated and delivered for physics. The commissioning is no longer driven by pushing the limits of the installation as designed but is in concurrence with the need for physics. It should be remembered that the SC LINAC installation was designed for a 200 kW deuteron beam (40 MeV / 5 mA) in continuous mode and the NFS experiments are requesting proton and deuteron beams of very low repetition rate. The constructed S3 experimental area is to be commissioned with beam within two years. Their needs are heavy ion beams from the A/q ≤ 3 ion source in a first time.

The first beam of 2021 was a deuteron beam of 20 MeV/A. Starting with the same procedure as in Fig. 5, the duty cycle went from 1 to 100 Hz using the SBS. A first experiment of 11 MeV/A deuteron beam on a Be and Li converter was performed in the middle of September in NFS followed by a 20 MeV/A beam on the rotational converter for the neutron beam production at NFS. Beam intensities of 47 µA was used during these tests.

During the shutdown periods, technical modifications have been performed and corrections identified during the run have been made. During the 2021 run one of the cavities went out of operation, to continue the program foreseen, the accelerator was tuned without this cavity. The tuning was performed and this mode validated before the correction of the breakdown during the winter break.

Experimental Areas for LINAC

The S3 installation will be ready for stable beams as of end of 2023. During the second semesters of 2022 and 2023 the facility shall proof the possibility to accelerate the ion beams of 48Ni, 48,50Ti, 50,54Cr with intensities up to 2 pμA for A/Q around 3. A new injector increasing the capacities for heavy ion production with a goal of A/q ≤ 7 are under study for GANIL. A superconducting ion source, a second RFQ for pre-acceleration and bunching should be installed. This project is under development, fully founded and should be installed and commissioned in 2026 as of the plans today.

CYCLOTRONS

The Cyclotron facility combined with the SPIRAL1 facility provides Stable Ion Beams (SIB) since 1983 and Radioactive Ion Beams (RIB) since 2001 for Physics experiments.

Two ECR4 ion sources, commissioned in 1992, are coupled to two cyclotrons (C01) and two separated section cyclotrons (SSC) where the beams can be successfully accelerated after striping and bunching in between.

At the cyclotron facility the SIB can be produced at three positions, any of the two injectors consisting of ECR4 (14.5 GHz) ion sources, or by the ECR booster directly.

The beam can therefore be delivered to five experiments in parallel if requested from three points in the production of the SIB. A low energy experimental area directly after the injectors (0.3 to 1 MeV/A), a medium energy experimental area (3.7-13.7 MeV/A) take the beams after the first SSC and finally after the full acceleration (24-95 MeV/A). The beam can then be shared by two experimental areas at the same time. In parallel, a SIB can be produced directly at the SPIRAL 1 installation, either from the ECR1+ ion source or the ECR booster, and sent to a dedicated experimental room.

If RIBs are requested, it can be delivered either as 1 + low energy or as post-accelerated beams from SPIRAL1, with the ISOL method. In one of the experimental areas there is the possibility to use a rotating or fixed target for beams up to 2 kW for fragmentation reactions with Inflight method to be used directly or on a secondary target.

The eight experimental areas at the cyclotron facility as seen on the right hand of the cyclotrons in Fig. 1 can take any of the SIB for experiments of fundamental nuclear physics, applied physics, industrial applications or any application needed SIB.

There are continuously ongoing upgrades and developments for the stable beams at the cyclotrons, driven by request from the physics.

The most recent developments are presented in [3] as the Tellurium beams and a 232 Th 30+ beam have been produced (results should be published in a short time), using LCO heating [START_REF] Lemagnen | New metallic stable beams for GANIL[END_REF]. Developments on adapting the MIVOC technique to the ECR4 sources are still ongoing to enlarge isotopes available to produce beams. It is the same group that work on the operation and development of the 4 injectors, p/d source, Phoenix V3 and ECR4 ion sources. The development of new beams and more stable operation are therefore mutualised.

Since the 2015 paper [START_REF] Kamalou | GANIL Operation Status and New Range of Post-Accelerated Exotic Beams[END_REF] an upgrade has been made of the SPIRAL1 facility shown in Fig. 7. This upgrade was motivated by the capacity to produce a larger span of radioactive ion beams at GANIL [START_REF] Dolegieviez | New exotic beams from the SPIRAL 1 upgrade[END_REF].

At the cyclotron facility there are two possibilities for production of RIBs, either through Fragmentation at the inflight facility LISE or by the ISOL method at the SPIRAL1 facility. The RIBs produced today [START_REF]Chart of beams at GANIL, stable and radioactive[END_REF] at the SPIRAL1 facility are seen in Fig. 8, this list is being updated regularly due to the ongoing developments.

The SPIRAL1 Installation

The SPIRAL1 installation, running in 2001-2014, underwent an upgrade to increase the capacity to deliver RIBs for the physics. This upgrade lead to changing totally the design of the target, ion source and the complete beam lines and install a Charge breeder to continue to permit the post accelerated beams for physics. At the SPIRAL1 installation the target is an ISOL target where the stable beams of 12 C at 95 MeV/A up to 238 U at 8 MeV/A from the SSC are stopped in a thick target (targets of Z < 41) or thin targets (no Z restriction). The elements are ionised either by an ECR or FEBIAD ion source. A charge breeder [START_REF] Delahaye | Charge breeding of light metallic ions: Prospects for SPIRAL[END_REF] allows to send all the RIBs to a post-accelerator CIME cyclotrons and use this beam for any of the eight experimental areas.

Since the upgrade a dense program for production of new elements is prepared every year and a new exciting exploration of these beams for physics are in the starting block opening up new scientific possibilities for the users of GANIL. Since 2019 the beams are produced with intensities of 106 pps. A 47 K 10+ @7 MeV/A beam was used for physics in 2021. In complementary to the noble gazes already produced earlier at GANIL beams of Na, Al, Mg, Cl, K, Rb have successfully been produced, the production of Br was recently produced but still under analysis while writing this article and are ready for further explorations for the needs for physics. The use of molecular beams, broken up in the charge breeder are under investigation and opens up further possibilities of RIBs for the users [START_REF] Maunoury | Charge Breeder at GANIL: metal chargebred elements[END_REF]. The charge breeder have been upgraded with variable 8-18 GHz RF injection at a power of 200 W. This allows to optimise further the charge state efficiencies and then the RIBs to be post-accelerated with CIME.

NEXT STEPS FOR GANIL

The installation of the S3 installation is ongoing and the installation should be ready for the first beams in end of 2023. The installation will in the first years use the beams from the A/q ≤ 3 ion source before going to the beams from the A/q ≤ 7 ion source (2026-2027).

In the end of the S3 installation is a low energy branch [START_REF] Ferrer | In gas laser ionization and spectroscopy experiments at the Superconducting Separator Spectrometer (S 3 ): Conceptual studies and preliminary design[END_REF] consisting of a gas cell setup. This part of S3 installation will be installed before the commissioning of the S3.

A seen in Fig. 9 the DESIR facility is to be installed between the current cyclotron and SC Linac facilities. The construction is foreseen to start in 2023, this installation consist of 150 m of transport beam lines from the SPIRAL1 and the S3-LEB installation towards the low energy experimental hall of DESIR. Exploration of the radioactive ion beams from either of the facilities will therefore be possible in this new experimental complexes as from 2026 after the full installation.

CONCLUSIONS

The last 5 years have been exciting for GANIL, upgrades of the current 40 year old cyclotron facility and the new installation of SPIRAL2 have largely increased the capacity of beams and the capacity of interesting physics to be performed at the installation. The next 5 years will follow by further installations and further commissioning of the installations already in place. As in the span of 5-10 years even more beams, stable and radioactive will come available for the users. For what is the program from 10-30 years from now, new plans are already discussed.

Figure 3 :

 3 Figure 3: SC LINAC layout. At the left hand side the two injectors coming together in the LEBT beam line. The RFQ and the MEBT beam line and then the 26 accelerating cavities and their cryo-modules, each one separated by a warm quadrupole section and a beam diagnostic.

Table 1 :

 1 The Technical Design Goals of the SPIRAL2 Installation an intensity of 100 μAe, and since then 58 Ni at 40 kV and 40 eµA have been produced, as requested by the S3 experimental facility. A strong and intense development program is ongoing for the injector, the focus of the next beams to be produced; extracted and accelerated are 48 Ni, 48,50 Ti, 50,54 Cr for the S3 experiments, with intensities up to 2 pμA for A/Q around 3 [3].

Figure 4 :

 4 Figure 4: Transmission in the RFQ of the three reference Q/A. installed and the line was physically connected to the LINAC.In the full MEBT line the three re-bunchers and a beam bunch selector (SBS)[START_REF] Giacomo | Proton beam commissioning of the SPIRAL2 Single Bunch Selector[END_REF] was installed. The SBS allows an operation of 1 bunch over 100 up to 1/10000, needed primarily for experiments requiring a precise timing between the bunches to insure that there is no overlap of the bunches for precise energy measurements.In the meantime the SC-LINAC installation was ongoing and in 2017 the cryogenic system for the cryo-modules were commissioned with a first complete cool down in November 2017[START_REF] Ghribi | Cryogenic Tests of the SPIRAL2 LINAC Systems[END_REF]. The RF tests in the LINAC were performed in parallel to the operation of the cyclotrons during the first semester in 2019.The GANIL cyclotron and SPIRAL2 installations are constructed and operated under French nuclear authority regulations. Before the authorisation of construction a file explaining the full operation and safety measurements were sent in. The different systems for the safety controls were then constructed and tested. The final authorisation to start the facility with acceleration of the beams in the LINAC was followed four months later after having complete testing and documentation of the different safety systems. Only after this the full authorisation was given and the beam could finally be sent to the LINAC.The SC LINAC is designed with 26 cavities, 12 A-type cryo-modules each comprising one low-β cavity and 7 B -type cryo-modules with two high-β cavities[START_REF] Olry | Spiral2 Cryomodules B Tests Results[END_REF][START_REF] Marchand | Performances of SPIRAL 2 Low and High Beta Cryomodules[END_REF].

2 .

 2 Increase of the Beam energy through the LINAC 3. Increase of the Beam intensity while opening the slits cutting the beam in the LEBT. 4. Increase of the duty cycle. During the phases 3 and 4 the beam losses were controlled in all the LINAC by the beam loss monitors (BLM) installed

Figure 5 :

 5 Figure 5: Transmission in the RFQ of the three reference Q/A. around the LINAC, the intensity (ACCTs) and the pressure variations all along the LINAC.The first beam, a 1 W proton beam, without energy increase and at low intensity was transported through the SC-LINAC tuned in a rebuncher mode in the middle of November 2019. During the following two weeks a beam energy increase up to 33 MeV with a pencil beam of 200 µA and a full 6.4 W beam was sent to the SAFARI beam dump. Before the winter shut down in December, this beam was sent to the NFS experimental area for full identification of the energy and composition of the beam.From here on, the LINAC operation is performed with the GANIL teams, also working on the operation of the Cyclotrons. The beam commissioning of the LINAC are therefore limited to 4-5 months every year at the moment.The phase three continued in the second semester of 2020, after the cyclotron run programmed for the first semester of 2020. After restarting the full complex, retuning the beam and do all the initial tests as in the end of 2019, the proton beam intensity was increased to 4.8 mA and the in the following week, the beam power were then increased to 16 kW (Fig.6) using a 10% duty cycle, with a beam transmission near of 100% into the LINAC and high energy lines.

Figure 6 :

 6 Figure 6: The beam power up to 16 kW proton beam.

15th

  Int. Conf. on Heavy Ion Acc. Technology HIAT2022, Darmstadt, Germany JACoW Publishing ISBN: 978-3-95450-240-0 ISSN: 2673-5547 doi:10.18429/JACoW-HIAT2022-WE2I1 WE2I1 Content from this work may be used under the terms of the CC BY 4.0 licence (© 2022). Any distribution of this work must maintain attribution to the author(s), title of the work, publisher, and DOI

Figure 7 :

 7 Figure 7: The SPIRAL1 for production of 1 + radioactive ion beams with the possibility of charge breeding in the ECR booster before post acceleration in the CIME cyclotron.

Figure 8 :

 8 Figure 8: List of beams from the SPIRAL1 facility.

Figure 9 :

 9 Figure 9: The full GANIL facility as from 2026.

15th Int. Conf. on Heavy Ion Acc. Technology HIAT2022, Darmstadt, Germany JACoW Publishing

Content from this work may be used under the terms of the CC BY 4.0 licence (© 2022). Any distribution of this work must maintain attribution to the author(s), title of the work, publisher, and DOI

Radioactive Ion Beam Facilities and other Facilities

ACKNOWLEDGEMENTS

I like to thank A.K. Orduz, G. Normand and M. Dubois for the help in writing this article.