domain: spi phys chim

text:
Valentin Davaze

Nicolas Vallino

Bertrand Langrand

Jacques Besson

Sylvia Feld-Payet
email: sylvia.feld-payet@onera.fr

Sylvia Feld

A non-local damage approach

Keywords: Mesh-dependency, non-local damage, dynamic explicit simulation, crash, parallel computing

come   L'archive ouverte pluridisciplinaire

Introduction

Numerical crack prediction in car crash simulations has been a major issue for automotive industry over the last decades. Indeed, structures are becoming increasingly thinner due to ecological concern and thus can be subjected to tearing in crash conditions. For safety and structural integrity reasons, tearing must be prevented. To do so finite element simulations are used for dimensioning during the car design process. However, ductile failure is not yet reliably predicted by the numerical computations. Then, when cracking is not numerically predicted but detected by the certification tests, parts and tools have to be re-designed, which increases the development cost and entails delays.

A previous work [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] has proposed to improve the numerical prediction of crack initiation and propagation with the use of a fully coupled damage model instead of fracture criteria.

The computed damage variable causes a decrease of the local element resistance, offering a better modelling of strain localization. In the same work, simulations of several geometries with different stress triaxialities tested at different loading rates have validated the performance of the proposed Gurson type model, in terms of both simulated load-displacement curves and crack path.

However, stress softening associated with coupled damage models is responsible for the localization of the damage variable on a row of integration points. The numerical prediction of cracking is thus strongly dependent on mesh size, mesh orientation and element type. This well-known difficulty requires to always use the same mesh size, then considered as constitutive model parameter, and sometimes to know and prescribe beforehand the crack path [START_REF] Besson | Modeling flat to slant fracture transition using the computational cell methodology[END_REF] so as to obtain the experimentally observed macroscopic response. The prediction capability of the model is therefore limited. Applying this strategy is not satisfactory for an industrial application since car parts may present very complex geometries. As the corresponding meshes are automatically generated, it becomes difficult to obtain a constant mesh size and to prescribe mesh orientation. In addition these meshes may use different types of elements (e.g. hexahedrons, prisms and tetrahedrons).

To obtain mesh-independent results, several regularization methods are proposed in the literature. For example, [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF][START_REF] Bažant | Crack band theory for fracture of concrete[END_REF] suggest to maintain a constant fracture energy G c , which is a material parameter, by adapting the elementary fracture energy according on the elements' size (L e ). However, this well-known Crack-Band method is more suitable for brittle materials, since fracture energy is difficult to estimate for ductile materials. Moreover, damage remains localized and results may still depend on the mesh orientation.

Other papers [START_REF] Ladevèze | On a damage mechanics approach[END_REF][START_REF] Allix | Delayed-Damage Modelling for Fracture Prediction of Laminated Composites under Dynamic Loading[END_REF][START_REF] Suffis | Damage model with delay effect: Analytical and numerical studies of the evolution of the characteristic damage length[END_REF] have proposed, from a mathematical analysis of localization in dynamic simulations, a delayed damage approach. This method suggests computing a new damage variable, delayed with respect to the original one, with a rate bounded [START_REF] Allix | The bounded rate concept: A framework to deal with objective failure predictions in dynamic within a local constitutive model[END_REF] by a given value to identify. Despite its easy implementation, since it only requires a rather simple modification of the constitutive model, this approach strongly influences the material behavior. In addition, the regularized band width depends on the loading rate. Thus, there is no direct control on the dissipated energy and material ductility may be overestimated at very high velocity if damage is over-delayed.

The last type of regularization approach groups the so called non-local methods. They consist in computing non-local damage variables at a given material point from the values in its neighborhood. These methods introduce a material characteristic length (referred to as ℓ in this paper) which affects the regularization band width (L r ) and enables the control of the dissipated energy independently of the loading rate. Results are thus mesh-independent as long as the mesh size is smaller than the regularization band width (3L e L r see [START_REF] Zhang | Ductile damage modelling with locking-free regularised gtn model[END_REF]).

A first non-local approach [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF] consists in replacing the usual damage variable with its weighted spatial average over a representative volume (Ω c ) of the same size as the material characteristic length (ℓ). The weighting function is often a Gaussian function but other forms can be used [START_REF] Pijaudier-Cabot | Nonlocal Damage Model with Evolving Internal Length[END_REF]. However, these non-local integral methods necessitate establishing a dialogue between Gauss points considering neighboring elements. These operations may compromise parallel computing, and necessitates a specific treatment at boundaries to account for the fact that the integration volume is truncated. They are believed to be very detrimental in the case of explicit computations with respect to the computation time as they do not follow the "element by element" treatment required to achieve efficiency.

To solve the numerical and algorithmic problems associated with the previous integral non-local method, it was proposed to use a Taylor expansion of the integral used to computed the non-local variables [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF][START_REF] Engelen | Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour[END_REF]. This leads to an equation introducing the second order derivative (Laplacian ∇ 2 ) of the non-local variable. The explicit form [START_REF] Mühlhaus | A variational principle for gradient plasticity[END_REF][START_REF] Borst | Gradient-dependent plasticity: Formulation and algorithmic aspects[END_REF] uses the Laplacian of the local variable to define its non-local counterpart, which limits the non-locality to an infinitesimal volume and is not suited to realistically describe dynamic problems and fracture problems according to Peerlings et al. [START_REF] Peerlings | A critical comparison of nonlocal and gradientenhanced softening continua[END_REF]. In contrast, the implicit form [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF] involves the gradient of the non-local variable and is thus truly non-local. It has been found to be equivalent to the integral non-local method [START_REF] Peerlings | A critical comparison of nonlocal and gradientenhanced softening continua[END_REF], but with the advantage of preserving the "element by element" treatment, so that algorithms for parallel computing remain unchanged.

Other solutions consider the enrichment of the free energy with the gradient of a variable which plays a non-local role. The internal variable gradient method [START_REF] Zhang | Ductile damage modelling with locking-free regularised gtn model[END_REF] considers one (or several) variable(s) of the local model, while the micromorphic approach [START_REF] Forest | Micromorphic approach for gradient elasticity, viscoplasticity, and damage[END_REF][START_REF] Saanouni | Micromorphic approach for finite gradient-elastoplasticity fully coupled with ductile damage: Formulation and computational aspects[END_REF] introduces one (or several) new variable(s) related to an already existing variable through an additional term in the free energy. If this additional term plays a penalization role to impose equality between an internal variable and its non-local counterpart, then the micromorphic model tends to be equivalent to the internal variable gradient method. In both cases, differentiation of the free energy leads to the same partial differential equation as the implicit second gradient methods for the non-local variable. Besides, for both methods, the adaption of the constitutive behavior ensues directly from the differentiation of the free energy, whereas the adaption is user-dependent for the implicit gradient method.

In this paper, non-local methods are preferred to regularize the constitutive model presented in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF]. Indeed, with these methods, the damage variable is not localized on a single row of elements which ensures independence on mesh orientation, unlike the Crack Band approach. Also, they should offer a better control of the dissipated energy via the internal length ℓ, independently to the loading rate, contrary to the delayed damage approaches.

Among the non-local methods, both the implicit gradient approach and the micromorphic approach enable to preserve the parallel computing, which is a crucial point for an industrial application. Extension of the micromorphic formulation to the explicit dynamics framework has recently been proposed by Diamantopoulou et al. [START_REF] Diamantopoulou | Micromorphic constitutive equations with damage applied to metal forming[END_REF]. Although these two non-local methods share some similarities regarding the discretized balance equations, it is worth noting that they lead to significant differences regarding the constitutive model. Comparison of the numerical responses with experimental results [START_REF] Davaze | Modélisation numérique de l'amorçage et la propagation des fissures dans les tôles métalliques ductiles pour les simulations de crash[END_REF] has eventually led the authors to focus this paper on the implicit gradient approach for which parameter identification is more straightforward, at least for the considered constitutive model.

Besides, contrary to the work of Diamantopoulou et al. [START_REF] Diamantopoulou | Micromorphic constitutive equations with damage applied to metal forming[END_REF] which deals with slow loadings problems (metal forming), this work considers faster loadings (crash simulations). For this loading rate range, the authors propose to improve the robustness of the considered regularization methods by introducing, in the non-local balance equation, a damping term.

The implicit gradient method has then been implemented by enriching classical linear underintegrated continuum elements of the Europlexus dynamic explicit solver. Since considering crash simulations, the updated lagrangian version of these Europlexus elements has been used. This enables to express the constitutive relations within the small strain framework.

In this paper, 3D solid elements (i.e. hexahedrons, prisms, and tetrahedrons) are considered, whereas Diamantopoulou et al. [START_REF] Diamantopoulou | Micromorphic constitutive equations with damage applied to metal forming[END_REF] used 2D plane strain elements.

In section 2 the local constitutive model obtained for the DP450 steel in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] is presented and used to illustrate the mesh-dependency problem on two practical examples. Then, the extended implicit second gradient method and its numerical implementation are described in section 3. The introduction of the additional non-local variable in the constitutive model is detailed in section 4. Finally, the implicit gradient formulation is used in section 5 to successfully simulate experimental results for five different specimens. Limits of the model are also discussed based on the comparison with a sixth case.

Note that all the simulations presented in this work use multi-domain parallel computing.

Damage constitutive model and mesh dependency

In this section, the constitutive model proposed in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF], which will be regularized later in this paper, is briefly described. The mesh dependency problem is then illustrated with several simulations realized with different mesh sizes and mesh orientations on specimens taken from the experimental campaign presented in the same paper.

Damage constitutive model

In a previous paper [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF], a Gurson type damage model was proposed to describe the behavior of dual-phase DP450 steel sheets. It uses the Gurson-Tvergaard-Needleman (GTN) yield criterion [START_REF] Gurson | Continuum Theory of Ductile Rupture by Void Nucleation and Growth: Part I-Yield Criteria and Flow Rules for Porous Ductile Media[END_REF][START_REF] Tvergaard | Influence of voids on shear band instabilities under plane strain conditions[END_REF][START_REF] Tvergaard | Analysis of the cup-cone fracture in a round tensile bar[END_REF] accounting for the effect of damage and pressure which is expressed as:

φ( σ, σ F , f * ) = σ σ F 2 -1 Plasticity + 2q 1 f * cosh q 2 2 Tr( σ) σ F -(q 1 f * ) 2 Effect of damage (1)
This yield criterion can be seen as composed of two parts. The plasticity part corresponds to an elasto-plastic criterion that compares an equivalent stress σ to the flow stress σ F . The Drucker equivalent stress [START_REF] Drucker | The Relation of Experiments to Mathematical Theories of Plasticity[END_REF] is used in this case, and its expression is given below:

σ( σ) = 6 1 27 -c 4 729 -1 (J 3 2 -cJ 2 3 ) (2)
where J 2 is the second deviatoric invariant and J 3 the third deviatoric invariant of the stress tensor. The expression of the flow stress describes the effect of work-hardening [START_REF] Maoût | Aluminum alloy damage evolution for different strain pathsapplication to hemming process[END_REF], strainrate dependency [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF] and thermal softening [START_REF] Zhao | A constitutive model for metals over a large range of strain rates identification for mild-steel and aluminium sheets[END_REF] due to quasi-adiabatic self-heating at high loading rates. It is expressed as:

σ F (p, ṗ, T ) = (R e + Hp + Q(1 -exp(-bp))) Work hardening 1 + C ln ṗ ε0 + Viscosity (1 -µ(T -T ref ))
Thermal softening

Strain-rate effects [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] where p designates the effective plastic strain, R e is the initial yield stress, H the linear hardening, Q and b the two Voce [START_REF] Voce | The Relationship between Stress and Strain for Homogeneous Deformation[END_REF] parameters, C the Johnson-Cook [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF] parameter, ε0 the inviscid limit strain-rate, µ the Zhao [START_REF] Zhao | A constitutive model for metals over a large range of strain rates identification for mild-steel and aluminium sheets[END_REF] thermal softening parameter, T the temperature, and T ref the reference temperature chosen as the room temperature (20 • C). The temperature evolution [START_REF] Roth | Effect of strain rate on ductile fracture initiation in advanced high strength steel sheets: Experiments and modeling[END_REF] is governed by a quasi-adiabatic conversion [START_REF] Taylor | The latent energy remaining in a metal after cold working[END_REF] of the plastic work into heat:

Ṫ = ω( ṗ) η k ρC p σ : εp (4)
The weight function ω is used to obtain a continuous transition between, on the one hand, isothermal conditions at low strain rates and, on the other hand, quasi-adiabatic conditions at high strain rates with purely adiabatic heat generation. Its expression is:

ω( ṗ) =            0 if ṗ < εis ( ṗ -εis ) 2 (3 εad -2 ṗ -εis ) ( εad -εis ) 3 if εis ≤ ṗ ≤ εad 1 if ṗ > εad (5)
The second part of the GTN yield criterion represents the effects of damage on the material behavior. The evolution of the damage variable f t is given by: ḟt = AH(p -ε d ) ṗ

Nucleation ḟn

+ (1 -f t )Tr( εp ) Growth ḟg + k w f t [1 -(cos(3θ)) 2 ] s : εp σ "Shear damage" ḟsh (6)
It incorporates a nucleation term using the Heaviside function H. Nucleation starts at a plastic strain equal to ε d . The nucleation rate is then constant (A). The void growth term is directly obtained from the plastic flow rule (mass conservation). Finally the model also uses the Nahshon and Hutchinson "shear damage" term [START_REF] Nahshon | Modification of the Gurson Model for shear failure[END_REF], enabling failure prediction under near pure shear conditions. The effective void volume fraction f * in the GTN criterion (see eq. ( 1)) is defined from f t and represents the coalescence of cavities that starts when f t exceeds a critical value f c :

f * =      f t if f t < f c f c + 1/q 1 -f c f r -f c (f t -f c ) if f t ≥ f c ( 7
)
Material failure is obtained when f t = f r which is the void volume fraction at failure.

Numerical implementation

The presented constitutive model is implemented in the Europlexus dynamic explicit finite element code using the NICE 1 "Next Increment Correct Error" explicit numerical scheme for plasticity proposed by Halilovič et al. [START_REF] Halilovič | NICE h : a higher-order explicit numerical scheme for integration of constitutive models in plasticity[END_REF].

The constitutive law is called to compute the new stress tensor σn from the strain increment ∆ εn . At the beginning of the time step, the strain increment is supposed to be purely elastic, and the trial stress is computed as follows:

σtrial = σn-1 + C : ∆ εn ( 8
)
Where C is the assumed isotropic stiffness tensor defined with the Young modulus E and the Poisson ratio ν. The corresponding value of the yield function is then evaluated:

φ pred = φ( σtrial n , σ n-1 F , f * n-1
). If φ pred ≤ 0, the material behavior is linear and elastic. Otherwise, i.e. φ pred > 0, a plastic correction must be applied to bring the yield function back to a value close to zero by updating the stresses and the internal variables. In this case, the NICE 1 scheme proposes to explicitly compute the plastic multiplier ∆λ by verifying an enhanced consistency condition:

φ n-1 + ∆φ n = 0 (9)
with:

∆φ n = ∂φ ∂ σ n-1 : ∆ σn + ∂φ ∂p n-1 ∆p n + ∂φ ∂f * n-1 ∆f * n + ∂φ ∂T n-1 ∆T n (10)
It has to be noted that this modified consistency condition introduces the previous value of φ to take into account the remaining error after correction created at the previous time step ∆t n-1 . The plastic multiplier at the current time step ∆t n will thus be computed by "correcting" this residual error. This self-correcting property helps to prevent the drift of the solution usually encountered with explicit plasticity algorithms. This drift is often corrected by performing costly iterations to ensure that the yield function remains close to zero. As drifting is prevented using the NICE scheme, no iterations are needed which helps to reduce the computation time.

Applying the normality rule, one has: ∆ εp = ∆λ ∂φ ∂ σ = ∆λ n where n is the normal to the yield surface. Following the GTN model, the plastic work is expressed as:

(1 -f t )σ F ∆p = σ : ∆ εp (11)
The effective plastic strain increment is therefore related to the plastic multiplier by a linear function involving a scalar parameter denoted a p :

∆p = σ : n (1 -f t )σ F ∆λ ≡ a p ∆λ (12)
All the internal variable increments depend either on the plastic strain increment tensor ∆ εp or on the effective plastic strain increment ∆p and can thus be rewritten introducing the plastic multiplier:

∆ σ = C : (∆ ε -∆λ n), ∆p = a p ∆λ, ∆f * = a f ∆λ, ∆T = a T ∆λ (13)
Where the scalars a f and a T can be derived from eq. ( 6) and eq. ( 7) on the one hand, and eq. ( 4) on the other hand, i.e.:

a f = ∂∆f * ∂∆f t AH(p -ε d )a p + (1 -f t )Tr( n) + k w f t [1 -(cos(3θ)) 2 ] s : n σ ( 14
)
and

a T = ω( ṗ) η k ρC p σ : ṅ (15)
Using eq. ( 9), the expression of the plastic multiplier can thus be obtained as:

∆λ = φ n-1 + n : C : ∆ ε n : C : n -a p ∂φ ∂p -a f ∂φ ∂f * -a T ∂φ ∂T (16)
Once ∆λ is calculated, the stress tensor and the internal variables are updated by computing the corresponding increments for time step ∆t n using eq. ( 13).

Due to mesh dependency, the implemented constitutive model was identified in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] for a constant mesh size of L e = 0.2 mm. The corresponding parameters are given in tab. 1.

E [GPa] ν R e [MPa] H [MPa] Q [MPa] b c 192.0 0.3 283.0 587.0 208.0 23.9 1.45 ε0 [s -1 ] C µ [ • C -1 ] η k ρ [kg.m -3 ] C p [J.kg -1 . • C -1 ] εis [s -1 ]
0.00361 0.0236 0.00134 0.9 7850.0 420.0 0.002

εad [s -1 ] ε d A k w f c f r f u 0.04 0.3 0.11 2.65 0.16 0.2 1/q 1 q 1 q 2 1.5 1.0
Table 1: Identified parameters for a fixed mesh size of L e = 0.2 mm.

Mesh-dependency phenomenon 2.3.1. Theory

If the simulations are performed with a mesh size different from the one used for identification (i.e. L e = 0.2 mm), a poorer numerical prediction might be observed. This phenomenon is due to the softening nature (due to damage or temperature increase) of the constitutive model which leads to the localization of strains in one row of elements or integration points. From a mathematical point of view, this localization can be explained, in the case of dynamic explicit simulations, by a loss of hyperbolicity of the discretized equilibrium equations [START_REF] Bažant | Wave propagation in a strain softening bar: Exact solution[END_REF][START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] De Borst | Fundamental issues in finite element analyses of localization of deformation[END_REF]. To illustrate it, a one-dimensional case problem is considered. The discretized equilibrium equation is in this case:

∂∆σ ∂x = ρ ∂ 2 ∆u ∂t 2 with ∆σ = κ∆ε ( 17
)
where κ is the tangent modulus (κ = E for a purely linear elastic behavior). Considering that ∆ε = ∂∆u/∂x, this equation can be rewritten as follows:

κ ρ ∂ 2 ∆u ∂x 2 = ∂ 2 ∆u ∂t 2 (18)
If the tangent modulus κ is positive (elasticity, work-hardening...), the equation above is hyperbolic and its solution is a wave that propagates at the speed κ/ρ. When the tangent modulus becomes negative (softening due to damage, self-heating ...), the equation becomes elliptic and its solution gets stationary. The propagation speed vanishes and the deformation wave is then trapped [START_REF] Wu | Deformation trapping due to thermoplastic instability in one-dimensional wave propagation[END_REF] in the first weakening element. This wave is no longer transmitted to neighboring elements. The damage variable will only evolve in the weakening element and damage localization is observed. As the elementary fracture energy depends on the elements size, the smaller the mesh size, the less "ductile" the material macroscopic response.

Therefore, spurious results might be obtained for different mesh sizes or mesh orientations.

Note that some works [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] De Borst | Fundamental issues in finite element analyses of localization of deformation[END_REF] have shown that introducing a plastic strain-rate dependency in the constitutive model helps to regularize the solution. The examples presented below will illustrate that the rate dependency introduced by the Johnson-Cook model [START_REF] Johnson | A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[END_REF] is not sufficient in this case.

Examples of mesh dependency

In the presented constitutive model, two softening variables may introduce mesh dependency: the damage variable f t and the temperature T at high velocities. To illustrate this phenomenon on practical examples, two specimens taken from the experimental campaign presented in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] are used: FN and V45 (sketches of the tests specimens are shown in section Appendix A). The FN specimen is relevant to see the mesh-dependence on both crack initiation and crack bifurcation due to its specific "S-shape" crack path. The V45 specimen is more interesting to observe the mesh-dependency on crack propagation. These specimens were tested with a high speed tensile machine at loading rates of respectively 1.0 m.s -1 and 3.0 m.s -1 , at room temperature. Simulations are performed using the dynamic explicit software Europlexus with the identified parameters presented in tab. 1 on three different mixed (hexahedrons and prisms) meshes for both specimen (see fig. 1).

The first mesh for each specimen (see fig. 1a and fig. 1d) is a regular mesh with a constant mesh size of L e = 0.2 mm corresponding to the identification size. For these meshes, the results are supposed to give the best predictions. The second mesh (see fig. 1b and fig. 1e) uses a smaller mesh size of L e = 0.15 mm and are disoriented on purpose (30 • for FN specimen and 50 • for V45 specimen) to observe a potential dependency to mesh orientation.

Finally, the third mesh (see fig. 1c and fig. 1f) is regular and uses a very small mesh size equal to L e = 0.1 mm. Usual symmetry conditions are used so that only half the thickness of the specimen is represented. For the three mesh designs, this half thickness is respectively divided in 3, 4 and 6 elements. Simulations are performed using under-integrated linear hexahedrons (CUBE) with viscous hourglass control and prisms (PRIS). and fig. 2e). As expected the total crack length is larger for the 0.15 mm mesh compared to the 0.2 mm mesh for the FN specimen. Regarding the FN disoriented mesh (see fig. 2b), the crack path shape seems to be in good agreement with the experiment, which leads to think that the local results are not influenced by the mesh orientation in this particular case. It is important to note, for the rest of the study, that this specimen has 3 crack initiation spots (see fig. 2a), which is also observed experimentally: two on the edges of the specimen, and one at the center. This gives a specific "S-shape" crack path. The fact that this crack path does not have the same orientation for the disoriented mesh cannot be attributed to mesh dependency. This is only due to small defects in the mesh that induce the crack to choose one orientation or the other. The scatter of crack path orientation with meshes could be fixed by using exact boundary conditions extracted from DIC pictures (as presented in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF]), that take into account a possible dissymmetry of the loading. In this work, nominal imposed displacement is used as DIC measurements are not available for these high velocities tests.

Regarding the V45 disoriented mesh (see fig. 2e), even if the crack path seems to be globally straight as expected, it shows local small deviations (zig-zag between element rows).

Moreover, the crack progression (i.e. its horizontal length) is surprisingly smaller for this disoriented mesh with L e = 0.15 mm than for the coarsest mesh with L e = 0.2 mm (see fig. 2d). A higher dependence to the mesh-orientation is therefore observed on the V45 specimen.

The load-displacement curves (see fig. 3) also show different results at failure for both specimens, particularly visible for the V45 specimen. The finest mesh presents an earlier crack initiation for the FN specimen than the two others (see fig. 3a). The displacement at failure (and thus the dissipated energy) tends to decrease with the increasing mesh fineness.

The same observation is made on the V45 curves (see fig. 3b). However, for this specimen, the disoriented 0.15 mm mesh shows a higher dissipated energy than observed for the 0.2 mm mesh, despite a quasi-identical final notch opening. This indicates that more energy is dissipated for the disoriented 0.15 mm mesh due to the zig-zag path (see fig. 2e).

Implicit second-gradient non-local approach

Mesh dependent results have been observed for both FN and V45 specimens. In order to obtain mesh-independent results, it is proposed in this section to use the recently proposed idea of [START_REF] Diamantopoulou | Micromorphic constitutive equations with damage applied to metal forming[END_REF] to obtain a modified version of the implicit second gradient equation suitable for dynamic explicit computations. Since the damage variable f t and the temperature T both depend on the plastic strain, it is chosen to regularize the effective plastic strain p.

The theory of the presented implicit second-gradient non-local approach is detailed. Its numerical implementation for continuum elements is then presented and illustrated.

Theory

In this paper, the following non-local implicit second gradient formulation is proposed:

   ℓ 2 ∇ 2 p nl -γ ṗnl + (p -p nl ) = ζ pnl on Ω ∇p nl . n = 0 on Γ ( 19
)
where p nl is the non-local value of the effective plastic strain and n is the external normal unit vector. The factor ζ is a non-local "density" (unit s 2 ). This term is not in the initial formulation of Peerlings et al. [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF] which only describes a "quasi-static" evolution on a nonlocal variable. Its addition follows the proposition of Diamantopoulou et al. [START_REF] Diamantopoulou | Micromorphic constitutive equations with damage applied to metal forming[END_REF] to take into account inertia effects for a micromorphic model. It makes the formulation compatible with the Newmark central differences scheme since it enables to evaluate the acceleration of the non-local degrees of freedom p nl . This formulation implies that the solution p nl is a wave (as the displacements) that might oscillate in the case of sudden variations of the local variable p, especially in the case of crash simulations. In order to avoid parasite oscillations in the p nl field, it is proposed in this work to introduce a non-local damping parameter γ (unit s) as a factor of the non-local variable velocity.

Non-local implicit second gradient methods are based on a rewriting [START_REF] Engelen | Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour[END_REF][START_REF] Peerlings | A critical comparison of nonlocal and gradientenhanced softening continua[END_REF] of the non-local weighted integral models proposed by Pijaudier-Cabot et al. [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Pijaudier-Cabot | Comparison of various models for strain-softening[END_REF]. The main advantage of gradient methods is that there is no longer a need to perform integrals over a patch of elements for each integration point. However new degrees of freedom are introduced which correspond to the "non-local" counterpart of the local variable. A standard finite element implementation is then possible as it preserves the "element by element" evaluation of internal forces. This approach also enables to preserve parallel computing which is an important issue for large size simulations performed in the automotive industry.

The implementation of this set of equations in an explicit finite element solver (Europlexus [START_REF] Team | Europlexus: a computer program for the finite element simulation of fluid-structure systems under transient dynamic loading[END_REF] in this case) is detailed in the following. Simulation are still performed with under-integrated linear continuum elements, but with additional degrees of freedom for the non-local variable p nl .

Numerical implementation for continuum elements

The weak non-local formulation is obtained by multiplying eq. ( 19) by a test function p *

The "non-local equilibrium" above shows the balance between the different non-local forces.

The non-local external forces contain the variable at the integration points to be regularized, i.e. p in this case. This external solicitation will induce the evolution of the non-local degree of freedom p nl .

Using the standard finite element procedures, the previous integrals can be expressed as a sum over all elements of the structure. For each element, the integral is evaluated using Gauss quadrature (1 Gauss point per element in the present case). Let {p e nl } be the vector of the nodal values for the non-local plastic strain. The values at the Gauss points of p nl and its gradient ∇p nl are evaluated as:

p nl = {N e }.{p e nl } ∇p nl = [B e ].{p e nl } (21)
where matrices {N e } and [B e ] contain respectively the shape functions and their spatial derivatives. As usual, an appropriate choice of the test functions among the shape functions in eq. ( 20) leads to the following elementary contributions: Standard assembly procedures are then used to obtain "forces" associated to each degree of freedom.

Description of the equilibrium procedure

The presented non-local equilibrium equation is used to compute the acceleration of the non-local variable, as shown in fig. 4 (and detailed below), which causes the non-local variable to evolve (as with any kinematic degree of freedom). Each node of the discretized structure now has a extra degree of freedom p nl in addition to to the nodal displacements.

Three steps are then necessary to calculate both sides of the equilibrium equation:

1. In the initialization of the computation (step 0. in the flow chart presented below in fig. 4), the global "mass matrix" is computed for the non-local variable considering elementary mass matrices. Considering eq. ( 22), the non-local elementary mass matrix is expressed as:

[M nl e ] = ζ Ωe {N e } T .{N e }dΩ (23)
Note that in dynamic explicit codes the mass matrix is often lumped to make it diagonal and so easier to inverse. The elementary non-local lumped mass matrix is defined as follows:

[M nl e ] L i,i = n j=1 [M nl e ] i,j (24)
where n is the number of nodes for the considered element. For the rest of the study, the mass matrix is supposed to be lumped, and the subscript L will be omitted. The non-local global mass matrix, which is lumped and so diagonal, will eventually be inverted to compute the acceleration of the additional degree of freedom in step 8.

The mass matrix associated to displacements (denoted [M ]

) is also computed and lumped.

2. At every time step, non-local forces are computed in each element (step 3. to 5. in fig. 4). As step 4. corresponds to the introduction of the non-local variable into the constitutive model, this step will be detailed in section 4. The local variable p is then supposed to be known at step 5. As the inertial forces shown in eq. ( 20) are not computed at the element level, the non-local elementary forces are defined so as to gather the internal, damping and external non-local forces:

{f nl e } = Ωe ℓ 2 [B e ] T .
When all the elementary contributions {f nl e } are computed, they are assembled in the same way as described for the mass matrix to obtain the global non-local forces vector {f nl } (see step 7.). Note that the dialogue between elements is realized at this step. Indeed, all elements with a common non-local degree of freedom will add their contribution to its evolution.

3. To finish, as the non-local variable p nl evolution is computed following the numerical scheme of Newmark central differences, it must respect a stability condition. A specific non-local critical time step must not be exceeded to preserve the numerical stability.

A stability analysis (see section Appendix B) on a one dimensional wire element leads to the following expression for the non-local critical time:

∆t nl = min e 2L e √ 3ζ 12ℓ 2 + L 2 e (26)
For the global computation stability, the minimum between the local mechanical equilibrium critical time step, referred to as ∆t meca , and the non-local critical time step ∆t nl is retained. In addition a safety factor (C stab < 1) is applied to this minimum time increment. In the following, C stab = 0.8 will be used.

∆t c = min(∆t meca , ∆t nl ) with ∆t meca = L e E/ρ (27)
It is important to note that the non-local critical time step depends on the value of the non-local parameters. From a physical point of view, the internal length ℓ generally remain of the order of magnitude of a few hundreds of micrometers and is thus not very critical. However, the parameter ζ must remain very low in order to prevent the "filtering" of the non-local variable evolution especially for very fast variations of p.

A very low value of ζ might increase the computation time. As the minimum between mechanical and non-local critical time steps is retained, it is chosen to fix ζ so that it does not penalize the computation time in comparison with a classical simulation without non-local regularization: This second gradient non-local approach is then implemented in the Europlexus solver by modifying some classical under-integrated linear continuum elements: hexahedral (CUBE), prisme (PRIS) and tetrahedral (TETR). The non-local version of these three element types will be referred to as CBNL, PRNL and TTNL (see fig. 5). Note that because of the under-integration, a viscous hourglass treatment is used for hexahedral elements.

∆t nl ≥ ∆t meca =⇒ ζ ≥ ρ(12ℓ 2 + L 2 e ) 12E ( 28

Non-local constitutive model

In this section, step number 4. of the algorithm presented in fig. 4 is detailed. For each element of the structure, the constitutive law (denoted F) is called to compute the stress tensor and the local effective plastic strain from the increment of total strain and the increment of non-local plastic strain:

( σ, p) = F(∆ ε, ∆p nl ) ( 29
)
Compared to a use with a classical element, the constitutive equations must then be modified to introduced the non-local plastic strain p nl . To do so, it is proposed to follow the Peerlings et al. [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF] approach. This method has also been used considering a non-local effective plastic strain by Engelen et al. [START_REF] Engelen | Nonlocal implicit gradient-enhanced elasto-plasticity for the modelling of softening behaviour[END_REF] and Feld-Payet et al. [START_REF] Feld-Payet | Finite element analysis of damage in ductile structures using a nonlocal model combined with a three-field formulation[END_REF], in combination with 0. Computation of the global lumped mass matrix:

[M] = [M e ] , [M nl ] = [M nl e ] Initialization of d.o.f.s: {ü} 0 = [M] -1 .{f ext } 0 , { u} 1/2 = { u} 0 + {ü} 0 ∆t 0 /2 {p nl } 0 = 0, { ṗnl } 1/2 = 0, {p nl } 0 = 0
For each element:

1. {∆u} n = { u} n-1/2 ∆t n {∆p nl } n = { ṗnl } n-1/2 ∆t n 2. New configuration: {u} n-1 + {∆u} n , {p nl } n-1 + {∆p nl } n
{f int e } n = V e [G e ] T . σn (⋆) {f nl e } n = V e ℓ 2 [
{f int } n = {f int e } n {f nl } n = {f nl e } n 8. {ü} n = [M] -1 .({f ext } n -{f int } n ) {p nl } n = -[M nl ] -1 .{f nl } n 9. { u} n+1/2 = { u} n-1/2 + {ü} n (∆t n+1 + ∆t n )/2 { ṗnl } n+1/2 = { ṗnl } n-1/2 + {p nl } n (∆t n+1 + ∆t n )/2 { ṗnl } pr ed n+1 = { ṗnl } n+1/2 + {p nl } n (∆t n+1 )/2 ∆t n+1 = C stab min(∆t e c )
Figure 4: Global algorithm used for non-local dynamic explicit simulations. Nodal data are in red and Gauss points data in blue. ( ⋆ : using under-integrated continuum elements, internal forces are computed using stresses at the single Gauss point and the element volume (V e ).) a solution to volumetric locking in the latter case. This practical method consists in the introduction of the non-local plastic strain p nl in the rate equations of the softening variables: i.e. damage f t and temperature T . Consequently eq. ( 6) and eq. ( 4) are modified to respectively compute the rates of the non-local damage and of the non-local temperature as follows:

d.o.fs              u x u y u z + p nl (a) CBNL element d.o.fs              ux uy uz + pnl (b) PRNL element d.o.fs              u x u y u z + p nl (c) TTNL element
ḟ nl t = AH(p nl -ε d ) ṗnl + (1 -f nl t )Tr( εnl p ) + k w f nl t w(θ) s : εnl p σ (30)
Ṫnl = ω( ṗnl ) η k ρC p σ : εnl p (31)
To calculate the non-local plastic strain rate tensor p from the non-local effective plastic strain ṗnl , a non-local plastic multiplier is computed using the Gurson energy equivalence:

(1 -f t )σ F ṗnl = λnl σ : n =⇒ λnl = (1 -f t )σ F ṗnl σ : n (32)
Note that this non-local plastic multiplier is only used to compute a non-local plastic strain tensor and remains different from the one computed with the NICE method for return mapping correction (see eq. ( 16)). The non-local plastic strain-rate tensor is then obtained using the normality rule:

εnl p = λnl n = (1 -f t )σ F ṗnl σ : n n (33)
The local effective damage f * and the local temperature T are then replaced by their nonlocal counterparts f * nl (obtained from f nl t using eq. ( 7)) and T nl in the yield criterion eq. ( 1)

and the flow stress eq. ( 3):

φ( σ, σ F , f * nl ) = σ σ F 2 -1 + 2q 1 f * nl cosh q 2 2 Tr( σ) σ F -(q 1 f * nl ) 2 (34)
σ F (p, ṗ, T nl ) = (R e + Hp + Q(1 -exp(-bp))) 1 + C ln ṗ ε0 + (1 -µ(T nl -T ref )) (35)
In the expression of the flow stress, the work-hardening and the rate-dependency, are not modified and still depend on the local plastic strain p. This implies that only the softening behavior due to heating or damage is modified. Plasticity and rate-dependency remain unchanged.

Using this approach, the numerical plastic algorithm scheme (presented in section 2.2)

needs to be adapted with the following two steps:

1. As the non-local plastic strain is computed outside the element, it is an input variable of the constitutive law and is considered as constant during material computation. The non-local damage and temperature are then updated at the beginning of the increment.

This operation is realized whether the material behavior is elastic or plastic.

f nl t = f nl t + ∆f nl t (∆p nl , ∆ εnl p ) and T nl = T nl + ∆T nl (∆p nl , ∆ εnl p , ṗnl ) (36)
2. As these two softening variables are not updated during the plastic correction but remain constant through the entire material computation process, they become weakly coupled with the element behavior. This is acceptable since dynamic computations use very small time steps. In case of plastic flow, the computation of the plastic multiplier is then modified since the only remaining internal variable is the local plastic strain p:

∆λ = φ n-1 + n : C : ∆ ε n : C : n -a p ∂φ ∂p (37)
The stress tensor σ and the local plastic strain p can be updated using eq. ( 13). This couple of outputs variables is then used for the computation of elementary forces (see step 5. of the flow chart given in fig. 4).

Simulation of experimental results

Parameters' choice

As explained in section 3.3, using the non-local approach, the width of the localization band should be greater than the largest considered element size within the investigated range (i.e. L e ≤ 0.2 mm). This guides the choice of the non-local parameter ℓ = 0.4 mm. The choice of the non-local parameter ζ is made so that it does not penalize the computation time: ζ = 1.0 × 10 -13 s 2 . Finally, to limit parasite oscillations, the non-local damping parameter γ is set to γ = 1.0 × 10 -6 s.

The parameters of the local constitutive model which gave results in agreement with the experiments with a mesh size of L e = 0.2 mm must be adapted to preserve this agreement.

In [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF], it has been mentioned that the main parameter that must be modified in the event of a change in localization bandwidth/element size is the nucleation slope A proportionally to the size of the localization band. The value A = 0.37 (instead of 0.11) is then retained.

The considered parameters of the entire non-local Gurson -implicit gradient constitutive model are then summarized in tab. 2.

In the end, the adaptation of the implicit gradient model parameters from the values identified for the local model remains straightforward, with only two different parameters that can be identified based on the load-displacement curves: the non-local internal length ℓ and the nucleation slope A, and two purely numerical parameters associated with the explicit framework: ζ and γ.

Simulation of experimental results for the FN and V45 specimens

To illustrate the performance of this approach, let us first consider the test cases of the FN and V45 specimens with different mesh sizes and orientations as presented above (see

2.3.2).

Similarly to the local formulation, when an element has reached complete failure, its mechanical internal forces are null. The explicit resolution with the Newmark algorithm then enables to naturally deal with the corresponding degrees of freedom for the displacements.

E [GPa] ν R e [MPa] H [MPa] Q [MPa] b c
192.0 0.3 283.0 587.0 208.0 23.9 1.45

ε0 [s -1 ] C µ [ • C -1 ] η k ρ [kg.m -3 ] C p [J.kg -1 . • C -1 ] εis [s -1 ]
0.00361 0.0236 0.00134 0.9 7850.0 420.0 0.002

εad [s -1 ] ε d A k w f c f r f u 0.04 0.3 0.37 2.65 0.16 0.2 1/q 1 q 1 q 2 ℓ [mm] ζ [s 2 ] γ [s] 1.5 1.0 0.4 1.0×10 -13 1.0×10 -6
Table 2: Identified parameters for a non-local constitutive model using implicit second gradient approach enabling convergence for L e ≤ 0.2 mm.

However, with the proposed non-local formulation, the evolution of the non-local degrees of freedom is still computed even for the failed elements.

The obtained damage fields for the three investigated mesh sizes are shown in fig. 6 at simulation time t = 6.95 ms for the FN specimen and t = 1.03 ms for the V45 specimen.

In the case of the FN specimen, the same damage field is obtained so that complete failure occurs at the same time for all mesh sizes and orientations. It must be noticed that the "S-shape" crack path of the FN specimen is preserved as well as the three crack initiation spots (two on the edges, one in the middle, see fig. 6a). In the case of the V45 specimen, the mesh orientation dependence appears to be corrected so that small deviations are no longer observed on the crack path (see fig. 6e). In both cases, the damage variable is much less localized than for the local model (see fig. 2). Maximal damage is reached on more than one row of integration points which is consistent with the non-local approach principle detailed in section 3.

Regarding the load-displacement curves fig. 7, numerical results with different mesh sizes are also quasi-identical for both crack initiation (see FN results in fig. 7a) and crack propagation V45 results in fig. 7b). The displacements at failure and dissipated energy are very close. The mesh size dependence as well as mesh orientation dependence (for V45 specimen) have thus been solved. Moreover, the simulated curves are in very good agreement with the experimental results.

Simulation of experimental results for additional specimens

To validate the proposed non-local modeling strategy, additional experimental results taken from [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] are simulated. Simulations are then compared with experimental results.

Sketches of these specimens can be found in section Appendix A. These specimens have been chosen to assess the performance of the non-local model for different stress triaxialities and strain-rates for both crack initiation and propagation. In all cases, simulations using mesh sizes equal to 0.1, 0.15 and 0.2 mm were investigated. The meshes having a minimum element size equal to 0.15 mm are disoriented similarly to what has been done for the FN and V45 specimens (see fig. 1).

The large flat tensile specimen tested at a loading rate of 1.0 m.s -1 is first investigated. This specimen has a large central area where strain remains homogeneous up to necking. A large portion of the specimen is therefore affected by adiabatic heating. This specimen is interesting to verify that the overall plastic behavior (which was fitted in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF] assuming no damage) is not modified using the non-local model. The damage fields given in fig. 8a, fig. 8b and fig. 8c show complete failure for the three meshes for the same simulation time, i.e. t = 9.0 ms. Maximal damage occurs on several rows of elements. The corresponding loaddisplacement curves (see fig. 10a) show a very good mesh convergence and good agreement with experiment.

The NT1 specimen tested at a loading rate equal to 3.0 m.s -1 is then considered. This specimen has the same cross section as the V45 specimen but has a larger notch radius equal to 1 mm. Its stress is therefore higher than the flat specimen triaxiality but lower than the V45 specimen triaxiality. The NT1 specimen is particularly interesting to evaluate both crack initiation and crack propagation. Moreover, the original local constitutive model was identified using this specimen [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF]. Damage fields (see fig. Finally, the M-shape specimen [START_REF] Roth | Ductile fracture experiments with locally proportional loading histories[END_REF] is simulated for a loading rate equal to 3.0 m.s -1 . This specimen enables testing the material at low triaxiality under shear. Damage fields at t = 1.0 ms (see fig. 9) and load-displacement curves (see fig. 10d) are almost identical for all three meshes. However a significant difference in terms of displacement at failure and dissipated energy is observed in comparison with experimental results. In the case of the M-shape specimens, the ligament between the notches is only 3 mm long. These notches also somehow constrain strain localization to occur in a band having a thickness of 0.5 mm (see fig. A.11) in which stress/strain gradients are very high even before damage is nucleated and starts to grow. The size of the highly deformed zone is of the same order of magnitude than the internal length ℓ used to regularize the problem, so that a size effect is observed.

The fact that ductility is overestimated indicates that a smaller internal length should have been chosen in order to closely reproduce the tests on the M-shape specimens. This is one of the drawbacks of regularized models which have difficulties capturing the physics of failure up to cracking which ultimately occurs in a very narrow band with regard to the internal length. Using a smaller internal length also implies using smaller elements (so as to obtain convergence) so that the size and duration of the calculations would be strongly increased.

A compromise needs to be found. The case of the M-shape specimen illustrates the fact that the failure of parts in which the failure zone size is constrained by the geometry will not be well described if their sizes is not large enough compared to the internal length. One possible way to determine ℓ would be to use homothetic specimens. As size effect should be observed for very small specimens and could be used to tune the internal length.

Conclusion

To improve the numerical prediction of metal sheets ductile tearing in crash simulations, the use of a coupled Gurson-type model was proposed in [START_REF] Davaze | Plastic and fracture behavior of a dual phase steel sheet under quasi-static and dynamic loadings[END_REF]. However, the strong mesh dependency (size, orientation) due to the softening behavior was an important shortcoming. This had to be solved to propose a reliable numerical crack prediction method that could be used for crash simulations in the automotive industry.

To obtain mesh-independent results, the implicit second gradient model was considered. This method requires an enrichment of under-integrated continuum elements with an additional degree of freedom, i.e. the non-local variable. For the considered model, it was proposed to compute a non-local effective plastic strain that can be used in the constitutive model to derive non-local counterparts of the two softening variables: damage and temperature. Compared to the original quasi-static formulation, two terms have been added in the non-local equation in order to comply with the dynamic explicit framework and deal with the possibly significant accelerations. With this method, parallel computing is automatically preserved, which is another important point for an industrial application.

Simulations with the implicit second gradient model were carried out on six geometries for which mesh independence (size or orientation) was indeed demonstrated. In addition, simulations were shown to be in good agreement with the corresponding experimental results, except for the M-shape shear specimen. In this particular case, discrepancy was attributed to the size of the highly deformed area in the M-shape specimen which was comparable to the internal length. This means that in the case of car crash simulations, which demand relatively coarse mesh sizes (since for explicit calculations the allowed time increment decreases as the mesh size gets smaller), the internal length ℓ may be adapted to the mesh size, but failure of zones relatively small compared to ℓ will no be properly described. This also suggests that a more physically based intrinsic internal length could be determined using homothetic specimens having very small sizes.

As thin metal sheets are often modeled using shell elements in car crash simulations, it would be interesting for further work to address the adaption of the non-local method so it is compatible with these elements.

To compute the mass matrix, these shape functions are computed at the two Gauss points The stability condition on time step (see eq. (B.4)) can be rewritten using the equation above. The following "non-local" critical time step is obtained: [START_REF] Borst | Gradient-dependent plasticity: Formulation and algorithmic aspects[END_REF]

∆t nl = 2 ω max = 2L e √ 3ζ 12ℓ 2 + L 2 e (B.

  40 mm (a) Regular -L e = 0.2 mm.

  Disoriented -L e = 0.15 mm. (c) Regular -L e = 0.1 mm. 24 mm (d) Regular -L e = 0.2 mm. 50 o (e) Disoriented -L e = 0.15 mm. (f) Regular -L e = 0.1 mm.

Figure 1 :

 1 Figure 1: Meshes of FN specimen (top) and V45 specimen (bottom) used to illustrate mesh dependency.

3Figure 2 :

 2 Figure 2: Damage fields obtained on FN specimen (1.0 m.s -1 at t = 6.9 ms) and V45 specimen (3.0 m.s -1 and t = 0.91 ms) with the local constitutive model (i.e. without regularization).

Figure 3 :

 3 Figure 3: Load-displacement curves obtained with different meshes on FN and V45 specimens illustrating the mesh dependency.

nl and integrating the

  result on the structure Ω. The result is analogous to the one obtained with the local mechanical equilibrium equation:Ω ℓ 2 ∇p nl . ∇p * nl + p nl p * nl dV

Ωe ℓ 2 [

 2 B e ] T .[B e ] + {N e } T .{N e } dΩ .{p e nl } + Ωe γ{N e } T .{N e }dΩ .{ ṗe nl } -Ωe {N e } T pdΩ = -Ωe ζ{N e } T .{N e }dΩ .{p e nl } (22)

)

  Considering a given value for ℓ and a possible range for L e (i.e. considering several meshes with different minimum mesh sizes), ζ must be computed using the largest value for L e if one wants a constant value for ζ. In the following, one uses ℓ = 0.4 mm and a maximum value for L e of 0.2 mm. This corresponds to a value for ζ equal to 6.67 × 10 -15 s 2 . However, in case of a large elongation due to loading, the element size L e might be increased. In this case, the critical value of ζ is also increased. That is why, as ζ is chosen as a constant, a conservative value of ζ = 1.0 × 10 -13 s 2 is retained for all the computations.No other modifications are required for a non-local parallel computation. The additional degrees of freedom are exchanged at the boundary of every domains like displacements. To summarize the implementation of the proposed method, fig.4shows the global algorithm including the modifications regarding a classic element necessary to introduce the implicit second gradient non-local approach. Standard procedures are used to time integrate displacements ({u}). The strain tensor increments are computed from the elementary nodal displacements thanks to the matrix [G e ] build from the shape functions' spatial derivatives as follows: ∆ ε = [G e ].{∆u e }.

3 .n 4 . 5 .

 345 Computation at Gauss point: ∆ εn = [G e ].{∆u e } n , ∆p n nl = {N e }.{∆p e nl } Calling constitutive law : ( σn , p n ) = F(∆ εn , ∆p n nl ) Computation of element forces:

Figure 5 :

 5 Figure 5: Non-local modified under-integrated elements.

Figure 6 :

 6 Figure 6: Damage fields obtained on FN specimen (1.0 m.s -1 and t = 6.95 ms) and V45 specimen (3.0 m.s -1 and t = 1.03 ms) with the non-local implicit second gradient approach.

Figure 7 :

 7 Figure 7: Load-displacement obtained with different meshes on FN and V45 specimens with the non-local implicit second gradient approach.

  8d, fig. 8e and fig. 8f) at t = 0.9 ms show quasi-identical crack lengths for the three meshes. The load-displacement curves for the different mesh sizes and orientations (see fig. 10b) are also very close and in relatively good agreement with experimental results. To evaluate the performance of the non-local model over large crack propagation, a Central Crack Panel (CCP) specimen is used. The loading rate is in that case equal to 7.0 m.s -1 . Crack lengths at t = 0.658 ms (see fig. 8g, fig. 8h and fig. 8i) are quasi-identicalfor the different mesh sizes and orientations. In particular, the results for the disoriented mesh (0.15 mm) are in good agreement with the other results. Displacement at failure and dissipated energy (see fig.10c) are very close illustrating a proper regularization in this case too. Relatively good agreement with experiment is also obtained.

29

 29

  (a) Lg. flat -L e = 0.2 mm. (b) Lg. flat -L e = 0.15 mm. (c) Lg. flat -L e = 0.1 mm. 20 mm (d) NT1 -L e = 0.2 mm. (e) NT1 -L e = 0.15 mm. (f) NT1 -L e = 0.1 mm. 20 mm (g) CCP -L e = 0.2 mm. (h) CCP -L e = 0.15 mm. (i) CCP -L e = 0.1 mm.

Figure 8 :

 8 Figure 8: Damage fields obtained on Large flat, NT1 and CCP specimen with the non-local constitutive model using the implicit second gradient approach.

Figure 9 :

 9 Figure 9: Damage fields obtained on M-shape specimen at the same simulation time t = 1.0 ms with the non-local constitutive model using the implicit second gradient approach.

  Exp. 3.0 m.s -1 0.20 mm regular 0.15 mm disoriented 0.10 mm regular (d) M-shape specimen.

Figure 10 :

 10 Figure 10: Load-displacement obtained with the non-local constitutive model using the implicit second gradient approach on several specimens.

which are located in ξ = - √ 3 /3 and ξ = √ 3 / 3 . 8 ) 11 ) 12 )

 33381112 Their integration weights are both equal to one. The non-local mass matrix [M nl e ] is then computed and lumped a follows:The matrix [B e ] contains the derivatives of the shape functions with respect to the globalThe characteristic equation (see eq. (B.5)) can be developed and the maximal eigenfrequency can be found by solving:By variable changing X = ω 2 n , this equation corresponds to a classical second order equation and the maximal solution is obtained for:

The shape functions of the reference element are: