M Chambart

R Desmorat

F Gatuingt

Intrinsic dissipation of a modular anisotropic damage model : application to concrete under impact

Keywords: Anisotropic damage, delay-effect, concrete, intrinsic dissipation, impact

come

Intrinsic dissipation of a modular anisotropic damage model : application to concrete under impact

Marion Chambart, Rodrigue Desmorat, Fabrice

Introduction

Damage mechanics is a powerful tool to handle micro-cracking and failure of quasi-brittle materials and structures. The corresponding constitutive models, even with induced anisotropy, are now most often implemented with localization limiters, i.e. either in a nonlocal form [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF][START_REF] Peerlings | Enhanced damage modelling for fracture and fatigue[END_REF] or as visco-damage models [START_REF] Needleman | Material rate dependance and mesh sensitivity in localization problems[END_REF][START_REF] Dub | Rate dependent damage model for concrete in dynamics[END_REF][START_REF] Allix | Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF] or both [START_REF] Voyiadjis | Thermodynamic framework for coupling of non-local viscoplasticity and non-local anisotropic viscodamage for dynamic localization problems using gradient theory[END_REF][START_REF] Gatuingt | Anisotropic 3D delay-damage model to stimulate concrete structures[END_REF]. This allows to gain some numerical robustness as mesh independency of the converged finite element solution.

Numerical robustness is nevertheless still difficult to obtain when induced anisotropic damage is considered, i.e. when the state of damage is represented by a tensorial variable. The modeling of damage anisotropy can be made by use of microplanes models [START_REF] Bazant | Microplane model for progressive fracture of concrete and rock[END_REF][START_REF] Bazant | Microplane model for brittle plastic material: I. Theory and II. Verification[END_REF][START_REF] Bazant | Nonlocal microplane model for fracture, damage, and size effect in structures[END_REF] or by use of a fourth [START_REF] Chaboche | Description thermodynamique et phénomenologique de la viscoplasticité cyclique avec endommagement[END_REF][START_REF] Leckie | Tensorial nature of damage measuring internal variables, Physical Non-Linearities in Structural Analysis[END_REF] or now more commonly of a second order damage tensor [START_REF] Cordebois | Endommagement anisotrope en élasticité et plasticité[END_REF][START_REF] Ladeveze | On an anisotropic damage theory[END_REF][START_REF] Murakami | Mechanical modeling of material damage[END_REF][START_REF] Mazars | The unilateral behavior of damage concrete[END_REF][START_REF] Chaboche | Development of continuum damage mechanics for elastic solids sustaining anisotropic and unilateral damage[END_REF][START_REF] Ladevèze | Modeling and simulation of the mechanical behavior of CMCs, High Temperature Ceramic-Matrix Composite[END_REF][START_REF] Dragon | An anisotropic model of damage and frictional sliding for brittle materials[END_REF][START_REF] Fichant | Isotropic and anisotropic descriptions of damage in concrete structures[END_REF][START_REF] Dragon | A model of anisotropic damage by mesocrack growth; unilateral effect[END_REF][START_REF] Carol | On the formulation of anisotropic elastic degradation. II. Generalized pseudo-Rankine model for tensile damage[END_REF][START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF][START_REF] Badel | Application of some anisotropic damage model to the prediction of the failure of some complex industrial concrete structure[END_REF][START_REF] Kitzig | Modeling of plain concrete structures based on an anisotropic damage formulation[END_REF]. The key point of a thermodynamics consistent modelling is to ensure the positivity of the intrinsic dissipation and in particular the positivity of the dissipation due to damage [START_REF] Lemaitre | Evaluation of dissipation, damage in metals submitted to dynamic loading[END_REF][START_REF] Jirasek | Mathematical analysis of strain localization[END_REF][START_REF] Lemaitre | Mécanique des matériaux solides[END_REF].

In case of isotropic (scalar) damage varable D, if Y stands for the (positive) thermodynamics force associated with D, the dissipation due to damage writes Ḋ = Y Ḋ so that one simply has to enforce Ḋ ≥ 0 from the damage evolution law [START_REF] Lemaitre | A course on damage mechanics[END_REF] in order to satisfy the positivity of the dissipation.

In case of anisotropic damage, two cases arise: a) the damage evolution law is written in the framework of standard generalized materials [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF], i.e. damage is proportional to the tensorial (positive) thermodynamics force Y or to J J J : Y with J J J a positive definite fourth order tensor [START_REF] Chaboche | Description thermodynamique et phénomenologique de la viscoplasticité cyclique avec endommagement[END_REF][START_REF] Leckie | Tensorial nature of damage measuring internal variables, Physical Non-Linearities in Structural Analysis[END_REF][START_REF] Kracinovic | Continuous damage mechanics revisited: basic concepts and definitions[END_REF][START_REF] Ju | On energy-based coupled elastoplastic damage theories: constitutive modeling and computational aspects[END_REF]. One has Ḋ = λ Y : J J J : Y ≥ 0, if λ stands for the (positive) damage multiplier.

The drawback of such a "standard" modelling is the difficulty to conceptualize the anisotropy represented in the different loading cases and the number of material parameters introduced as components J i jkl .

b) a practical possibility for quasi-brittle materials is to consider damage but also its anisotropy gouverned by the extensions [START_REF] Mazars | The unilateral behavior of damage concrete[END_REF][START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linaire et la rupture du bton de structure[END_REF] setting for example the damage rate proportional to the positive part of the strain tensor.

Email address: desmorat@lmt.ens-cachan.fr (R. Desmorat)

Preprint submitted to Engineering Fracture Mechanics July 1, 2014

Such a modelling is non standard, the positivity of the dissipation is not guaranteed for any models and has to be mathematically proven [START_REF] Desmorat | Positivity of intrinsic dissipation of a class of nonstandard anisotropic damage models[END_REF].

One considers here a modular anisotropic damage model, non standard (case b), for which the proof of the positivity of the dissipation can be mathematically given [START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF]. The proof, recalled in section 2.4, avoids the calculation of the tensorial thermodynamics force Y associated with second order damage D. Using this feature, one proposes in present work differents ways to determine within finite elements computations the dissipation for the considered anisotropic damage model. As pointed out in [START_REF] Jirasek | Mathematical analysis of strain localization[END_REF], the quantification of the energy dissipated in the rupture of structures made of quasi-brittle materials gives complementary information to the usual damage maps.

Thermodynamics of anisotropic damage

A state of micro-cracking at the Representative Volume Element (RVE) scale can be represented by a thermodynamics state variable, namely the damage variable. For instance, in concrete-like materials the micro-cracks are mainly orthogonal to the loading direction in tension and parallel to the load in compression. A second order damage tensor is then often considered to be able to represent such an induced anisotropic damage pattern [START_REF] Cordebois | Endommagement anisotrope en élasticité et plasticité[END_REF][START_REF] Ladeveze | On an anisotropic damage theory[END_REF][START_REF] Murakami | Mechanical modeling of material damage[END_REF][START_REF] Mazars | The unilateral behavior of damage concrete[END_REF][START_REF] Dragon | A model of anisotropic damage by mesocrack growth; unilateral effect[END_REF][START_REF] Papa | Anisotropic damage model for the multi-axial static and fatigue behavior of plain concrete[END_REF][START_REF] Fichant | Isotropic and anisotropic descriptions of damage in concrete structures[END_REF][START_REF] Lemaitre | Anisotropic damage law of evolution[END_REF][START_REF] Carol | On the formulation of anisotropic elastic degradation. II. Generalized pseudo-Rankine model for tensile damage[END_REF][START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF][START_REF] Badel | Application of some anisotropic damage model to the prediction of the failure of some complex industrial concrete structure[END_REF]. The choice to use of a second order tensor for damage -D of components D i j is made here. Such a choice is not a necessity (a fourth order tensor can be used for instance [START_REF] Chaboche | Description thermodynamique et phénomenologique de la viscoplasticité cyclique avec endommagement[END_REF][START_REF] Leckie | Tensorial nature of damage measuring internal variables[END_REF]) but it has proven very helpful to solve the incompatibility with thermodynamics of the consideration of two damage variables, one for tension and one for compression: damage represents the state of micro-cracking whatever the sign of the loading so that only one damage variable represents the microcracking pattern [START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF].

Hydrostatic or mean damage is then

D H = 1 3 tr D (1)
In order to simplify further expressions, a specific integrity tensor H is defined [START_REF] Cordebois | Endommagement anisotrope en élasticité et plasticité[END_REF][START_REF] Ladeveze | On an anisotropic damage theory[END_REF][START_REF] Papa | Anisotropic damage model for the multi-axial static and fatigue behavior of plain concrete[END_REF][START_REF] Lemaitre | Anisotropic damage law of evolution[END_REF]]

H = (1 1 1 -D) -1/2 (2)
where to take power α (here -1/2) of a symmetric tensor, one makes it first diagonal, one takes then the power α of the diagonal components and one finally turns back the tensor obtained in the initial working basis.

State potential and state laws

Following [START_REF] Desmorat | Modéle d'endommagement anisotrope avec forte dissymétrie traction/compression. (RF)2B -5es Journées Scientifiques du Regroupement Francophone pour la Recherche et la Formation sur le Béton[END_REF][START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF], the permanent strains are neglected and the free enthalpy density is written

ρψ ⋆ = 1 + ν 2E tr H σ D H σ D + 1 -2ν 6E tr σ 2 + 1 -ηD H + tr σ 2 - (3)
with E the Young modulus, ν the Poisson ratio, ρ the density, η the hydrostatic sensitivity parameter, and where (.) D = (.) -1 3 tr(.) 1 1 1 stands for deviatoric part of a tensor, . + and .

-respectively for positive and negative parts of a scalar, x + = max(x, 0) and x -= min(x, 0).

The elasticity law derives from state potential (3),

ǫ = ρ ∂ψ ⋆ ∂σ = 1 + ν E H σ D H D + 1 -2ν 3E tr σ + 1 -ηD H + tr σ -1 1 1 (4)
The thermodynamics force associated with damage D is

Y = ρ ∂ψ ⋆ ∂D ( 5
)
Its calculation is not an easy task (see the Appendix). Key-point of present work, one will avoid to have to calculate such a derivative with respect to D: one will only use the simpler derivative with respect to tensor H, with the property

∂ ∂H tr H σ D H σ D = 2 σ D H σ D (6)
which defines a symmetric second order tensor. The elasticity law coupled with anisotropic damage (4) can be inverted as

σ = E 1 + ν H -1 ǫ H -1 - (1 1 1 -D) : ǫ 3 -tr D (1 1 1 -D) + E 3(1 -2ν) (1 -ηD H ) tr ǫ + + tr ǫ -1 1 1 (7)
with then H -1 = (1 1 1 -D) 1/2 . Note that the term in factor of shear modulus E/(1 + ν) = 2G is deviatoric (traceless), the second term -in factor of bulk modulus K = E/3(1 -2ν) -being purely hydrostatic. Previous equation can also be expressed in terms of components as

σ i j = E 1 + ν H -1 ik ǫ kl H -1 l j - (δ kl -D kl )ǫ kl 3 -D kk δ i j -D i j + E 3(1 -2ν) (1 -ηD H ) ǫ kk + + ǫ kk -δ i j ( 8
)
If D = 0, equation [START_REF] Bazant | Microplane model for progressive fracture of concrete and rock[END_REF] gives

σ = E 1 + ν ǫ - tr ǫ 3 1 1 1 + E 3(1 -2ν) tr ǫ 1 1 1 (9)
which is the classical elasticity law without damage. The elasticity law ( 4) is the derivative with respect to stresses of the free enthalpy density (3). The Helmholtz free energy ρψ is defined either as the Legendre transform ρψ(ǫ, D) = sup σ σ σ σ : ǫρψ ⋆ (σ σ σ, D) , or in an equivalent manner as the primitive of elasticity law [START_REF] Bazant | Microplane model for progressive fracture of concrete and rock[END_REF] rewritten in terms of strains and damage. This defines :

ρψ = E 2(1 + ν) tr H -1 ǫ H -1 ǫ - [(1 1 1 -D) : ǫ] 2 3 -tr D + E 6(1 -2ν) (1 -ηD H ) tr ǫ 2 + + tr ǫ 2 - ( 10
)
so that one classically has

σ = ρ ∂ψ ∂ǫ Y = -ρ ∂ψ ∂D (11)

Criterion functions

For concrete-like materials, the elasticity domain can advantageously be defined in terms of strains [START_REF] Mazars | Application de la mécanique de l'endommagement au comportement non linaire et la rupture du bton de structure[END_REF]. Let us consider Mazars criterion function

f Mazars = ε -κ ε = ǫ + : ǫ + = I ε I 2 + ( 12
)
where . + stands for the positive part of a tensor in terms of principal values [START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF], ε I are the principal strains, ε I + are the (positive) extensions. The elasticity domain corresponds to f < 0. In the time independent case, damage growth occurs at the consistency condition f = 0 & ḟ = 0 and quasi-static concrete behavior is modeled by an adequate choice of consolidation function κ, either function of the trace of damage tensor (initial formulation [START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF]) or as a function of the active damage d act [START_REF] Souid | Pseudodynamic testing and nonlinear substructuring of damaging structures under earthquake loading[END_REF].

For dynamics and impact applications, the visco-damage framework is used [START_REF] Dub | Rate dependent damage model for concrete in dynamics[END_REF] as it can represent the observed strain rate effect on apparent ultimate stress. The criterion function will then in a equivalent manner be defined either as

f = g( ε) -tr D ( 13
)
in case of initial trace formulation, or as f = g( ε)d act [START_REF] Ladeveze | On an anisotropic damage theory[END_REF] in case of active damage formulation. In both cases, function g = κ -1 is the inverse of function κ.

Visco-damage framework states that the elasticity domain still corresponds to f < 0 but damage growth occurs at

f = D v ≥ 0 ( 15
)
D v is the viscous damage (it acts as the viscous stress in visco-plasticity). A visco-damage law has to be defined. Following [START_REF] Gatuingt | Anisotropic 3D delay-damage model to stimulate concrete structures[END_REF] a delay-damage law which bounds the damage rate can be used, expressed next either in terms of trace rate or in terms of active damage rate

D v = D v H(tr(ǫ)) tr Ḋ or D v = D v H(tr(ǫ)) ḋact ( 16
)
where Heaviside function of the trace of the strain tensor allows to model a material rate effect in tension and not in compression (H(x) = 1 if x ≥ 0, H(x) = 0 else). In compression one will assume that the rate effect experimentally observed is due to inertial effect that confine the sample tested [START_REF] Tang | Rate Effects in Uniaxial Dynamic Compression of Concrete[END_REF][START_REF] Gatuingt | Coupled damage and plasticity modeling in transient dynamic analysis of concrete[END_REF].

In both cases the same function D v (x) is used [START_REF] Allix | Delayed-damage modelling for fracture prediction of laminated composites under dynamic loading[END_REF],

D v (x) = - 1 b ln Ḋ∞ -x Ḋ∞ (17)
Ḋ∞ (equal to the invert of a characteristic time) and b are material parameters.

Active damage

The active damage is defined as the normed projection of the full damage tensor on the positive strain tensor ǫ + of components ǫ + i j ,

d act = D : ǫ + max I (ε I ) = D i j ( ǫ + ) i j max I (ε I ) (18)
To better understand this concept of active damage, we can illustrate it with a simple tension/compression test. In such a particular case the tensors are diagonal but note that Eq. ( 18) handles general 3D cases with non diagonal tensors. In a monotonic tension test performed in direction 1 we obtain the damage matrix

D =           D 1 0 0 0 0 0 0 0 0           d act = D 1 (in simple tension test) (19)
this allows to recovers d act = tr D as in the initial anisotropic damage model [START_REF] Desmorat | Modéle d'endommagement anisotrope avec forte dissymétrie traction/compression. (RF)2B -5es Journées Scientifiques du Regroupement Francophone pour la Recherche et la Formation sur le Béton[END_REF][START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF].

In monotonic compression test performed in direction 1 we obtain

D =           0 0 0 0 D 2 0 0 0 D 2           d act = 2D 2 (in simple compression test) (20)
d act = tr D is again recovered. But if we considered an alternate loading, for example a simple tension test in direction 1 follows by a simple compression test (still in direction 1), the damage state and the active damage become

D =           D 1 0 0 0 D 2 0 0 0 D 2           d act = 2D 2 (
in the compression phase after damaging in tension) [START_REF] Fichant | Isotropic and anisotropic descriptions of damage in concrete structures[END_REF] in this case d act tr D showing that damage D 1 due to tension is not active in compression anymore in alternate loading. The same feature stands for D 2 due to compression which becomes inactive in tension as then d act = D 1 if the sign of the loading path change one more time. This concept is useful to deal properly with alternated loading [START_REF] Souid | Pseudodynamic testing and nonlinear substructuring of damaging structures under earthquake loading[END_REF][START_REF] Desmorat | Non standard thermodynamics framework for robust computations with induced anisotropic damage[END_REF], as illustrated in Section 4.4.

Non standard damage evolution law and positivity of the dissipation

For standard generalized materials [START_REF] Halphen | Sur les matériaux standards généralisés[END_REF], the evolution laws for the internal variables are derived by generalized normality. Good mathematical properties of the corresponding evolution potential, such as convexity with respect to the thermodynamics forces, allows to naturally satisfy the positivity of the intrinsic dissipation Ḋ = σ : ǫρ ψ, which allows then to naturally satisfy the second principle of thermodynamics.

More modeling flexibility is obtained when the damage evolution law does not derive from an evolution potential. For instance the non standard evolution laws

Ḋ = λ ǫ + α (with then Ḋ ≥ 0) (22)
of damage anisotropy governed by the positive part of the strain tensor (usually at a power α = 1 or 2) have shown efficient to represent concrete micro-cracking pattern in different loading cases. The damage rate is gained as a positive tensor, λ is the (positive) damage multiplier determined either from the consistency conditions f = 0 & ḟ = 0 in time independent modeling or from the viscosity law f = D v in visco-damage modeling.

The corresponding damage model is non standard so that one has to prove the positivity of the intrinsic dissipation Ḋ. As only one internal variable is introduced it is expressed as Ḋ = Y : Ḋ. Such a proof has been given in [START_REF] Desmorat | Positivity of intrinsic dissipation of a class of nonstandard anisotropic damage models[END_REF][START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF] and part of it is the basis for the numerical computation of the dissipation proposed in next section. It uses the splitting of the state potential (4) into a shear (deviatoric) contribution and a hydrostatic (spherical) one. It also uses the fact that one does not need to express Y variable itself but only the contracted product Y : Ḋ which will later be recognized as

Ḋ = Y : Ḋ = 1 + ν E σ D H σ D : Ḣ + 1 -2ν 6E tr σ 2 + d dt g(D H ) g(D H ) = 1 1 -ηD H ( 23
)
with g a positive monotonic function of the hydrostatic damage D H = 1 3 tr D, so that dg dt ≥ 0. The right inside term of Eq. ( 23) is therefore positive. It uses the fact that due to Ḋ ≥ 0 and definition H = (1 1 1 -D) -1/2 the tensor Ḣ is also a positive tensor, so that if (σ D H σ D ) is found positive its contracted product with tensor Ḣ ≥ 0 (the right inside term of Eq. ( 23)) is positive, completing then the proof for positive dissipation Ḋ ≥ 0 for any loading, monotonic or not, uniaxial or 3D, proportional or non proportional.

The positivity of tensor (σ D H σ D ) is obtained from the positivity of its eigenvalues µ I (corresponding eigenvectors are denoted T I ) solution of

σ D H σ D T I = µ I T I (no sum) (24)
or after left multiplication by positive symmetric tensor H and then by transposed eigenvector T T I

T T I H σ D 2 T I = µ I T T I H T I (no sum) (25)
or

µ I = T T I H σ D 2 T I T T I H T I ≥ 0 (no sum) (26)
Eigenvalues µ I are ratios of positive scalars and are then positive, as announced.

Computation of the energy density dissipated by anisotropic damage

In the studied case of elasticity coupled with second order tensorial damage, the determination of the thermodynamics force associated with damage needs the derivative of (1 1 1 -D) -1/2 terms with respect to D. This can be handled in a closed form (see Appendix) but leads to a quite complex expression for intrinsic dissipation Ḋ.

One prefers to recognize in the derivative of the state potential with respect to the time two terms: the rate at constant damage of an elastic energy density and the intrinsic dissipation Ḋ,

ρ ψ⋆ = ρ ∂ψ ⋆ ∂σ : σ + ρ ∂ψ ⋆ ∂D : Ḋ = ǫ : σ + Y : Ḋ Ḋ ρ ψ = ρ ∂ψ ∂ǫ : ǫ + ρ ∂ψ ∂D : Ḋ = σ : ǫ -Y : Ḋ Ḋ ( 27
)
One can foresee then two ways to express the intrinsic dissipation :

• Directly, from the damage variable change as Ḋ = Y : Ḋ,

• Indirectly, as Ḋ = ρ ψ⋆ǫ : σ = σ : ǫρ ψ.

Expressions for the intrinsic dissipation

For the direct calculation, let us rewrite both the free enthalpy ρψ ⋆ and the elasticity law by making appear function g(D H ) defined in Eq. ( 23),

ρψ ⋆ = 1 + ν 2E tr H σ D H σ D + 1 -2ν 6E g (D H ) tr σ 2 + + tr σ 2 - ( 28
) ǫ = ρ ∂ψ ⋆ ∂σ = 1 + ν E H σ D H D + 1 -2ν 3E (g (D H ) tr σ + + tr σ -) 1 1 1 (29)
Take then the derivative of potential ρψ ⋆ with respect to time,

ρ dψ ⋆ dt = 1 + ν E H σ D H D : σ + 1 + 2ν 3E (g(D H ) tr σ + + tr σ -) tr σ + 1 + ν E (σ D H σ D ) : Ḣ + 1 -2ν 6E tr σ 2 + : d dt g(D H ) (30)
which enables to identify first term of previous equation as (see Eq. ( 27))

1 + ν E H σ D H D : σ + 1 + 2ν 3E (g(D H ) tr σ + + tr σ -) tr σ = ρ ∂ψ ⋆ ∂σ : σ (31)
so that Ḋ is second term of Eq. ( 30). It is obtained as Eq. ( 23) which is the analytical expression for the intrinsic dissipation by the direct method.

For the indirect method, one prefers to invert the elasticity law and to calculate the dissipation as

Ḋ = σ : ǫ -ρ ψ (32)
where the free energy ρψ given by Eq. ( 10) is the dual of free enthalpy ρψ ⋆ .

Choice of a discretization

The implementation in a finite element computer code requires a proper time discretization of the energy dissipated by damage per unit of volume D = Y : Ḋ dt for the direct method and of D = (σ : ǫρ ψ) dt for the indirect method. Over one time step ∆t = t n+1t n between times t n and t n+1 , one sets

D n+1 = D n + ∆D (33)
For the direct method the increment of energy dissipated over a time step is:

∆D = t n+1 t n Y : Ḋ dt. (34)
In the case of indirect method, one calculates the energy dissipated D n+1 = D(t n+1 ) as

D n+1 = ǫ ǫ ǫ n 0 σ : dǫ w e n + ǫ ǫ ǫ n+1 ǫ ǫ ǫ n σ : dǫ -ρψ(ǫ, D) n+1 0 with ρψ(0 0 0,0 0 0)=0 (35)
where w e = ǫ ǫ ǫ 0 σ : dǫ, w e n = w e (ǫ n ) is set. D n+1 is expressed as :

D n+1 = w e n + ǫ ǫ ǫ n+1 ǫ ǫ ǫ n σ : dǫ -ρψ n+1 ( 36
)
where ρψ n+1 = ρψ(ǫ n+1 , D n+1 ). For these two methods, direct and indirect, two simple numerical schemes for the time discretization are compared in Fig. 1: the rectangle rule and the trapezium rule.

Schemes for computation with direct method

From Eq. ( 23) and by discretization with the rectangles method (upperscript DER for Direct Exact Rectangles), one gets :

∆D = ∆D DER = 1 + ν E (σ D H σ D ) n+1 : ∆H + 1 -2ν 6E tr σ n+1 2 + g D n+1 H -g D n H ( 37
)
and with the trapezium rule (DET pour Direct Exact Trapezium):

∆D = ∆D DET = 1 + ν E (σ D H σ D ) n+1 + (σ D H σ D ) n 2 : ∆H + 1 -2ν 6E tr σ n+1 2 + tr σ n 2 2 g(D n+1 H ) -g(D n H ) (38)
with

∆H = H n+1 -H n = (1 1 1 -D n+1 ) -1 2 -(1 1 1 -D n ) -1 2 (39)
In both cases, a difference between nonlinear functions of damage D at time t n and time t n+1 appears. Numerically, this operation may cause important round up errors and one will prefer to approach these differences by developments at the first order in order to replace them by products and derivatives. For instance, the first order development for g(D n+1 H )g(D n H ) gives:

∆(g(D H )) = g(D n+1 H ) -g(D n H ) ≃ g ′ (D n H )∆D H = η (1 -ηD n H ) 2 (D n+1 H -D n H ) (40)
For tensor ∆H, one starts by expressing H n+1 as a function of damage D n and of damage increment ∆D :

H n+1 = (1 1 1 -D n+1 ) -1 2 = (1 1 1 -D n -∆D) -1 2 H n+1 = (1 1 1 -D n ) -1 4 1 1 1 -(1 1 1 -D n ) -1 2 ∆D (1 1 1 -D n ) -1 2 -1 2 (1 1 1 -D n ) -1 4 ( 41
)
H n+1 ≃ (1 1 1 -D n ) -1 4 1 1 1 + 1 2 (1 1 1 -D n ) -1 2 ∆D (1 1 1 -D n ) -1 2 (1 1 1 -D n ) -1 4
One now substitutes the new expression of H n+1 in ∆H:

∆H = H n+1 -H n ≃(1 1 1 -D n ) -1 4 1 2 (1 1 1 -D n ) -1 2 ∆D (1 1 1 -D n ) -1 2 (1 1 1 -D n ) -1 4 ≃ 1 2 (1 1 1 -D n ) -3 4 ∆D (1 1 1 -D n ) -3 4 ( 42
)
One finally gets two new expressions for the dissipation increment for the direct method. For the rectangle discretization (DAR for Direct Approached Rectangles):

∆D = ∆D DRA = 1 + ν E (σ D H σ D ) n+1 : 1 2 (1 1 1 -D n ) -3 4 ∆D (1 1 1 -D n ) -3 4 + 1 -2ν 6E tr σ n+1 2 + η (1 -ηD n H ) 2 (D n+1 H -D n H ) (43)
and for the trapezium one (DAT pour Direct Approached Trapezium) :

∆D = ∆D DT A = 1 + ν E (σ D H σ D ) n+1 + (σ D H σ D ) n 2 : 1 2 (1 1 1 -D n ) -3 4 ∆D (1 1 1 -D n ) -3 4 + 1 -2ν 6E tr σ n+1 2 + + tr σ n 2 + 2 η (1 -ηD n H ) 2 (D n+1 H -D n H ) (44)

Schemes for computation with indirect method

For the indirect method, the integral terms are also discretized by means of the rectangle and trapezium rules. This leads to two additional expressions for the dissipated energy density:

• for the rectangle rule (IR for Indirect Rectangles) one gets:

D n+1 = D IR n+1 = D n + σ n+1 : ∆ǫ -ρψ(ǫ n+1 , D n+1 ) (45)
• for trapezium rule (IT pour Indiret Trapezium):

D n+1 = D IT n+1 = D n + σ n+1 + σ n 2 : ∆ǫ -ρψ(ǫ n+1 , D n+1 ) (46)

Temperature increase

The complete heat transfer equation reads in the general case [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF] :

ρC Ṫ = k∇ 2 T + σ : ǫ p -A K : VK + T ∂σ ∂T : ǫe + ∂A K ∂T : VK (47)
where ∇ 2 is the Laplacian operator, C is the specific heat of the material, k is the thermal conductivity. The internal thermodynamics variables are noted V K and the associated variables (or thermodynamics forces) A K . In studied case of elasticity coupled with anisotropic damage without permanent strains (ǫ p = 0, ǫ e = ǫǫ p = ǫ), one has K = 1 and V 1 = D is the only internal variable, associated with A 1 = -Y. For low temperature changes (as it will be further shown to be the case), the simplification of no thermo-mechanical coupling ∂σ σ σ ∂T : ǫ + ∂A 1 ∂T : V1 = 0 can be made. The heat transfer equation becomes then:

ρC Ṫ ≃ k∇ 2 T + Y : Ḋ = k∇ 2 T + Ḋ (48)
Furthermore, under adiabatic conditions -which are justified in fast dynamics and impact applications -the temperature rate Ṫ and increase from room temperature ∆T become

Ṫ ≃ Y : Ḋ ρC = Ḋ ρC or ∆T ≃ D ρC (49)
The temperature increase is simply the energy density dissipated due to damage divided by material constant ρC.

Dissipated energy density evolution at RVE scale

The implementation in a Finite Element code of previous methods and schemes allows to compute first at the Representative Volume Element (RVE) scale, therefore at every Gauss point of a structure, the intrinsic dissipation Ḋ (as increment ∆D) and its evolution D(t) (as D n ) which is the evolution of the energy dissipated by damage per unit of volume. When there will be no risk of confusion (units precised for instance), the terminology "dissipation" will also stands for dissipated energy density D.

Anisotropic damage model -DYN and INI versions

Let us first recall the constitutive equations of the anisotropic delay-damage model used for concrete and which belongs to second order damage framework of Section 2. The value η = 3 is set for the hydrostatic sensitivity parameter [START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF].

• Elasticity ǫ = 1 + ν E (1 1 1 -D) -1/2 σ D (1 1 1 -D) -1/2 D + 1 -2ν 3E tr σ + 1 -tr D + tr σ -1 1 1 (50)
• Damage criterion

f = g( ε) -tr D g( ε) = aA atan ε a -atan κ 0 a ( 51
)
where ε is Mazars equivalent strain and where κ 0 is the damage threshold and A and a are the damage parameters.

• Visco-damage loading-unloading conditions,

f < 0 → elastic loading or unloading f = D v ≥ 0 → damage growth (52)
• Induced damage anisotropy governed by the positive extensions,

Ḋ = λ ǫ + 2 (53)
• The delay-damage evolution law (in terms of trace of damage tensor)

H(tr(ǫ)) tr Ḋ = Ḋ∞ 1 -exp (-b (g ( ε) -tr D)) ( 54
)
from which is determined the damage multiplier λ. Delay-damage law ( 54) is derived from visco-damage law [START_REF] Ladevèze | Modeling and simulation of the mechanical behavior of CMCs, High Temperature Ceramic-Matrix Composite[END_REF], it bounds the damage rate to Ḋ∞ and deactivates the delay (or strain rate) effect for negative hydrostatic strains. H(x) is Heaviside function.

When D v = 0 (or either Ḋ∞ → ∞ or b → ∞), the damage multiplier is determined from consistency condition f = 0 & ḟ = 0 instead of f = D v and the time independent initial model of [START_REF] Desmorat | Modéle d'endommagement anisotrope avec forte dissymétrie traction/compression. (RF)2B -5es Journées Scientifiques du Regroupement Francophone pour la Recherche et la Formation sur le Béton[END_REF][START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF]

is recovered (version designated next as INI).

The material parameters of the anisotropic delay-damage model for concrete are those of Table 2.

Add ρ = 2400 kg• m -3 , C = 880 J• kg -1 / • C

Monotonic tension and monotonic compression

This set of parameters leads to the response of Figure 2(a) in stress-strain plane for uniaxial monotonic quasi-static tension and compression. The delay-damage law considered with the viscosity parameters b and Ḋ∞ of Table 2 leads to a sensitive strain rate effect for strain rates larger than 1 s -1 (Fig. 2(b)). Below 1 s -1 , the stress-strain response in tension is found quasi independent from the strain rate. Recall that due to Heaviside function H(tr ǫ) in delay-damage law there is no material strain rate at all in compression (such a strain rate effect is assumed to be due to inertial confinement). As expected one observes that all the methods converge to the same result when the strain increment decreases, which validates the implementations. In both monotonic tension and compression cases, the DET method (Direct Exact Trapezium) is the one that presents the smallest difference between the total intrinsic dissipation at failure (maximum principal damage equal to 1, plateau reached) with either the large or the small strain increments. This method is the one retained for further structures analyses.

At failure, the energy dissipated by damage in tension is found 25 times lower than in compression (1000 J/m 3 in tension compared to 25000 J/m 3 in compression).

Strain rate effect on intrinsic dissipation

Let us focus now on the effect of the strain rate in tension on the dissipation (Fig. 2(b)). The strain rate effect introduced through the delay-damage law increases the apparent strength (or ultimate stress) in tension when the strain rate increases. The model rupture strain (strain at vanishing stress due to principal damage equal to 1) also increases. The intrinsic dissipation is therefore quite larger in dynamics tension than is quasi-static tension. The evolution of the dissipation is plotted for different strain rates in Figure 5. As for tensile strength, there is no substantial enhancement in the dissipation below the strain rate of 1 s -1 . When the strain rate increases, the dissipation evolution slows down (hence the denomination "delay-damage") but reaches at failure a much larger value (plateau of Fig. 5). At ε=10 s -1, the dissipation at failure worths 3100 J/m 3 and when it worths 71000 J/m 3 at 100 s -1. The energy dissipated is so multiplied by a factor of 3 at rate 10 s -1 (it versus a factor of 2 for the strength computed at the same strain rate) and it is multiplied by factor 71 at 100.s-1 (versus factor 10 for the strength at this strain rate).

The strain rate effect is more influent on the dissipation that on the peak stress. This phenomenon is also well captured with mesoscopic approaches [START_REF] Gatuingt | Numerical determination of the tensile response and the dissipated fracture energy of concrete: role of the meso-structure and influence of the loading rate[END_REF].

Micro-cracks closure effect on dissipation under alternated loading -INI-d act and DYN-d act model versions

Under alternated loading, such as those encountered in seismic or in impact applications, the micro-cracks opened in tension close in the compressive stages so that the stiffness is recovered [START_REF] Mazars | The unilateral behavior of damage concrete[END_REF][START_REF] Chaboche | Development of continuum damage mechanics for elastic solids sustaining anisotropic and unilateral damage[END_REF]. The anisotropic damage model of Section 4.1 represents a partial stiffness recovery thanks to the partition of the hydrostatic part into a positive part (affected by damage through parameter η = 3) and a negative one (not affected). To gain a complete stiffness recovery, a recovery not only of the bulk modulus but also of the shear modulus, it is necessary to do the same operation on the deviatoric part of the potential (to split this tensorial part into positive and negative parts and properly affect them differently by tensorial damage D). This additional splitting makes the numerical scheme much more complex (special positive parts have to be defined [START_REF] Ladeveze | On an anisotropic damage theory[END_REF][START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF]) for a relatively small difference on the result [START_REF] Lebon | Analyse de l'endommagement des structures de génie civil: techniques de sous-structuration hybrides couplées un modèle d'endommagement anisotrope[END_REF]. The choice to only keep a partial stiffness recovery is made here.

Consolidation function κ(tr D) of Eq. ( 12) acts as isotropic hardening in metal plasticity and increases by too much to yield stress in compression after a damaging stage in tension. To model correctly the behavior of concrete under alternated loading, a possibility is to deactivate in the criterion function f Mazars or in an equivalent manner in f the damage due to tension when going to compression (i.e. to consider to corresponding micro-cracks as closed). In that purpose, a new quantity called active damage, has been introduced and defined as Eq. [START_REF] Dragon | A model of anisotropic damage by mesocrack growth; unilateral effect[END_REF]. In an improved delay-damage model [START_REF] Desmorat | Non standard thermodynamics framework for robust computations with induced anisotropic damage[END_REF] it replaces tr D term in criterion function f as

f = g( ε) -d act d act = D : ǫ + max I (ε I ) ( 55
)
For monotonic loading, d act is equal to tr D (see Section 2.3), so that there is no need to re-identify the model material parameters (such as κ 0 , A and a) which remain those of Table 2.

The consequences of this first modification on the alternated stress-strain response are presented in Figure 6. This is the case of a compression loading (up to failure, loading in direction 1) after a damaging tensile stage (the damage due to tension reaches D 1 = D xx = 0.5 before unloading). The stress-strain response in monotonic compression is reported as a reference curve not to cross for previously damaged material A non physical response is once more encountered due to tr D term but this time present in initial delay-damage law [START_REF] Verpeaux | CASTEM 2000 : une approche moderne du calcul des structures Calcul des structures et intelligence artificielle[END_REF]. A too large overstress is obtained for the dot curve without damage deactivation of Fig. 7(a) where a dynamics tension stage is considered after damaging compression up to D 2 = D 3 = 0.5.

So in order to build such a version DYN-d act of anisotropic damage model able to deal with alternated loading in dynamics, the active damage concept also has to be used for delay-damage law, setting viscosity law as D v = D v H(tr(ǫ)) ḋact with D v (x) function defined by Eq. ( 17), so that

H(tr(ǫ)) ḋact = Ḋ∞ 1 -exp (-b (g (ε) -d act )) ( 56
)
with the fully satisfactory continuous black response "With damage deactivation" of Fig. 7(a). The response for the monotonic dynamics tension is reported as a reference in this figure, at the same strain rate in tension of 10 s -1 . The dissipation evolution with respect to strain is given in Fig. 7(b) for the different modeling choices. Note that due to the damaging stage in compression the (plateau) value for the dissipated energy density at rupture is larger than the one obtained at the same strain rate in the reference monotonic tension case (see Fig. 5).

- Again when the delay-damage law is not introduced (D v = 0), the version of the modular anisotropic damage model obtained is a time independent one (active damage acting if criterion function f , Eq. ( 55), damage multiplier determined from consistency condition f = 0 & ḟ = 0) and will be referred to as version INI-d act .

Dissipation in non proportional loading: example of Willam test

The schemes and numerical implementation for the calculation of dissipation D is incremental and 3D so that non-proportional loading can be handled.

Let us perform classical (numerical) Willam test which submits a RVE to tension/shear with significant rotation of the principal strain directions [START_REF] Willam | Fundamental issues of smeared crack models[END_REF]. This test is now widely used to verify that no stress-strain discontinuities are introduced by complex 3D constitutive damage models for concrete [START_REF] Rots | Computational modeling of concrete fracture[END_REF][START_REF] Oliver | Isotropic damage models and smeared crack analysis of concrete[END_REF][START_REF] Feenstra | The Rankine-plasticity model for concrete cracking[END_REF][START_REF] Meschke | A damage-plasticity model for concrete accounting for fracture-induced anisotropy[END_REF][START_REF] Carol | On the formulation of anisotropic elastic degradation. II. Generalized pseudo-Rankine model for tensile damage[END_REF][START_REF] Badel | Application of some anisotropic damage model to the prediction of the failure of some complex industrial concrete structure[END_REF][START_REF] Matallah | Inelasticity-damage-based model for numerical modeling of concrete cracking[END_REF]. It consists in a two stage applied loading: in a first stage, a simple tension loading is applied in direction z ≡ 3, bringing the material to the onset of softening. In a second stage, the increments of the strain components ε zz , ε xx and ε xz are prescribed in the fixed proportions 1/ 1.5/ 1. This represents an increment of positive (tensile) strain for both principal directions, accompanied by a rotation of principal axes.

The model response in terms of stresses is presented in Figure 8(a) and the dissipation evolution D(t) is plotted in Figure 8(b). This illustrated the abilities to handle complex loading of both the anisotropic delay-active-damage model and the proposed numerical schemes for the dissipation computation. This also illustrates in non-proportional loading that the intrinsic dissipation Ḋ remains positive (as D increases).

First structural application -3D dynamics tension test by spalling

The first 3D application example considered is the numerical simulation of a dynamics tensile test by spalling, corresponding to the experimental test performed by Brara and Klepasczko [START_REF] Klepaczko | An experimental method for dynamic tensile testing of concrete by spalling[END_REF]. This test consists in submitting a concrete specimen to an incident compression wave created by the impact of a projectile on the measuring bar (Fig. 9). The compression stresses generated by the incident wave are not sufficient to bring the specimen to failure. The reflection of this wave on the concrete specimen free face generates a tensile wave that breaks it into two parts. In this experimental apparatus the concrete specimen diameter is 40 mm and its length 120 mm.

Modeling parameters

The numerical simulation is carried out with CEA Finite Element code Cast3M [START_REF] Verpeaux | CASTEM 2000 : une approche moderne du calcul des structures Calcul des structures et intelligence artificielle[END_REF]. This code is based on a Newmark implicit integration scheme with a correction due to material non-linearities [START_REF] Verpeaux | Dynamique du solide : Modification du schéma de Newmark aux cas non linéaire, 6ème Colloque National en Calcul des Structures[END_REF]. In order to avoid numerical problem due to the damage model used, a specific treatment of the numerical control of rupture as in [START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF]. Only the specimen is mesh and the applied loading consists in prescribing a displacement history on one of the specimen faces (called the impacted face in the following), corresponding to the integration of experimentally measured velocity of this face. Using displacements for the applied load allows us to prevent convergence issues of the implicit Newmark integration scheme used. The mesh of the specimen is made up of 13632 prismatic elements (Fig. 10). The designation of the different versions are summarized in Table 3. The influence on dissipated energy of the delay-damagei.e. of strain rate effect -modeling and of micro-cracks closure by damage deactivation are studied next. The material parameters used for these simulation are the ones given in Table 2 ( Ḋ∞ and b delay-damage parameters are considered only for DYN model versions).

Model

Comparison between damage and dissipation maps

Damages maps (Fig. 11(a)) and dissipated energy density or dissipation maps (Fig. 11(b)) can be compared. The most dissipative areas correspond in all cases to the localization and failure zone. For the time independent damage models (no strain rate effect, INI and INI-d act model versions), the dissipation is localized in one or two rows of elements (the computation is pathologically mesh dependent). The other examples with DYN and DYN-d act model versions illustrate the role of delay-damage as a localization limiter. Nice feature, and as discussed in detail in [START_REF] Desmorat | Non standard thermodynamics framework for robust computations with induced anisotropic damage[END_REF], the computations regularization is obtained in dynamics with the value of the material (delay-damage) parameters identified for concrete from the real material strain rate effect (in other words it is not an unphysical viscosity put for regularization purposes).

Part of energy dissipated by damage

The total energy dissipated by damage or total dissipation D = Ḋ(t)dV is obtained by spatial integration of the dissipation density Ḋ over the whole mesh. This energy is compared to the injected energy E inj which corresponds to the integration over time of the product of the prescribed velocity v by the resulting force F on the impacted face

E inj = t 0 F(t) v(t) dt (57)
The total injected energy is calculated as E ∞ inj = 3.49 J at the end of the computation. Figure 12 shows for the different model versions the evolution at different times of the total dissipated energy. The ratio D/E inj (in percent) of the energy injected that is dissipated by damage is plotted in Figure 13. The evolutions superimpose perfectly for the different model versions until about 0.06 ms, i.e. until the time at which the maximum loading compressive stress reaches the free face. After, the effects of the modeling choices are clearly visible.

• The strain rate effect increases the value of the total energy dissipated by damage (up to almost 1.4 J with delay-damage DYN model versions, Fig. 12).

• In both INI and DYN cases, without and with delay-damage, i.e. without and with strain rate effect modeling, the damage deactivation leads to a decrease of the total dissipated energy, which is consistent with the fact that the tensile strength is smaller in this last case (see Fig. 7(a) and 7(b) and associated comments).

• The part of energy dissipated by damage fluctuates between 8 and 40 % of the injected energy. The later the failure occurs (the more distant of the free face), the larger the dissipated energy is.

• If one compares the total energy dissipated for the INI-d act model (D = 0.25 J) to the fracture energy per unit area for concrete G c , the values are of the same order of magnitude: with G c = 81 J/m 2 after [START_REF] Guinea | Measurement of the fracture energy using three-point bend tests: Part 1 Influence of experimental procedures[END_REF] and for our specimen of cross section S spec , the corresponding (total) fracture energy is S spec × G c = 0.11 J. Values which are found, due to dynamics effects, much lower than the 1.4 J obtained with delay-damage DYN model versions.

In order to validate the results but also to help for the choice of a model, the computed energies dissipated should be compared to experimental results. Unfortunately, in the tests used here [START_REF] Klepaczko | An experimental method for dynamic tensile testing of concrete by spalling[END_REF] the energy dissipated by the dynamic tensile tests was not measured. This kind of measure was made in [START_REF] Weerheijm | Research developments and experimental data on dynamic concrete behaviour[END_REF] and a dynamic increase factor of 12.5 was found for the higher loading rate (1685 GPa/s) which gives G c = 1505 J/m 2 in that case. However, in this kind of indirect measurement of the dynamic fracture energy, the results have to be used with precaution due to the inertia effects in dynamics that artificially increase the energy measured. That is why in this first numerical approach of the dissipated energy computed with the anisotropic damage model we only focus on the order of magnitude of the results obtained and on the the fact that the dissipated energy can be computed even in the complex case of loading induced anisotropic damage.

Second structural application -Quasi-statics and dynamics Brazilian tests

Brazilian test or tensile test by splitting (Fig. 14(a)) consists in applying a compression force on lateral height of a cylindric concrete specimen (or any material exhibiting tension/compression dissymmetry). This creates in the middle of the specimen a tensile stress state that leads to the complete rupture of the specimen. To perform such a test is a classical way to determine the tensile strength of concrete.

In quasi-statics, the compression force is usually applied by a cylindrical plate (Fig. 14(a)), while in dynamics the loading corresponds to the impact of a projectile dropped off from a certain height (Fig. 14(b)).

Quasi-statics brazilian test

The mesh, made of 8000 cubic elements, and the boundary conditions are represented in Figure 15. The (x, y, z) basis of Figure 15 is used for components plotting all along this example. The applied loading consists in a prescribed vertical displacement (parallel to z-axis) identical for the upper edge points of the specimen. The modular anisotropic damage model is used still with the material parameters of Table 2.

Figure 16(a) represents the damage map of component D yy , that is mainly micro-cracks orthogonal to y-direction. It exhibits strain and damage localization in a major tensile crack along the vertical loading axis as observed experimentally. In Figure 16(b), the isovalues of D zz -for micro-cracks mainly perpendicular to z-direction -show two small zones damaged in compression (micro-cracks parallel to the load), at the top right under the loaded edge and at the bottom along the support x = y = 0.

Figure 17 shows that the most dissipative areas correspond to the small zones damaged in compression (Figure 16(b)) and not to the middle zone (damaged in tension) which leads to the main crack and to failure. This is due to the fact, as seen in Section 4, that in quasi-statics it requires 25 times more dissipated energy to reach the failure in compression than in tension. A small compression damage involves a high value for the dissipation which can be larger than the total value for the dissipation needs to break the specimen in tension.

In the same Figure, the scale of temperature increase ∆T is reported (from Eq. 49). The maximum temperature elevation computed only reaches 0.01 • C, which is so small that it can not be experimentally measured for comparison.

Dynamics brazilian test

As already mentioned, the dynamics brazilian test is performed experimentally on a drop weight tower. The C.E.A. "ORION" drop tower enables to drop a 100 kg projectile from a given height h (max 8 m), which worths in the The simulation is carried out with C.E.A. explicit Finite Element code Europlexus [START_REF]Elements finis en dynamique rapide. Manuel de Rfrence Europlexus -CCR/CEA/EDF R&D[END_REF] dedicated to dynamics and impact applications. This computer code is able to handle the contact between two solids, one can then explicitly model the impact of the projectile on the specimen. Although in the Newmark explicit scheme used there is no convergence issues, we also used the numerical control of rupture in the integration of the anisotropic damage model as in Cast3M.

The projectile is meshed with 900 prismatic elements and its behavior is supposed linear elastic. The material parameters are the ones of steel, i.e. ρ = 7850 kg/m 3 , E = 210 GPa and ν =0.3. The 3D mesh of the specimen is made of 17600 cubic elements (Figure 18). The damage model used for concrete is again the anisotropic delay-damage model with the material parameters of Table 2. For both the projectile and the tested specimen the elements are fully integrated (i.e. with 6 Gauss point for prismatic elements and 8 for the cubic ones). One authorizes sliding without friction between the two solids. The time step is automatically adapted by the code to satisfy at each time step the Courant condition [START_REF] Courant | ber die partiellen Differenzengleichungen der mathematischen Physik[END_REF].

As in quasi-statics, the ruin of the specimen is caused by a major traversing tensile crack parallel to the vertical load direction z. The map of damage D yy (representing mainly micro-cracks parallel to z-direction) is presented in Figure 19. While in quasi-statics case the crack initiated from the center of the specimen, it can be seen that here, in dynamics, the major crack is the result of the coalescence of two cracks initiated from the edges. When focusing on the intrinsic dissipation map (Fig. 20), one notices that the most dissipative areas are once again the one loaded in compression. The maximum value for the energy dissipated by damage is found 4 times larger than in quasi-statics.

If one compares the experimentally measured displacement of the projectile to the one obtained by computation (Figure 21), one observes that the maximum (in absolute value) displacement is very well obtained but that the rebound is experimentally slower that numerically. The energy dissipated in the computation is then of the correct order of magnitude but it is too low by a factor of 3, feature which can be partly explained by the fact that the anisotropic damage model used in this paper does not take into account the permanent strains and friction that can occur along the cracks created. This may be important and represents an important source of energy dissipation. Furthermore in that simulation the contact between the two solids is considered as a sliding without friction contact, which does not correspond exactly to reality since the concrete right under the impact is crushed, which involves a non negligible friction.

Conclusion

This paper presents, from a theoretical and a numerical point of vue, the computation of intrinsic dissipation due to damage for an anisotropic 3D delay-damage model. Different methods and temporal discretizations are proposed in order first to verify the analytical expressions established but also to choose the most efficient numerically. Simulations on a single cubic elements allows to estimate the volumic intrinsic dissipation at failure, which worths 1000 J/m 3 in One has A A A : Ḣ = [H 2 ⊗ H 2 ] : Ḋ, rewritten A i jkl Ḣkl = H 2 ip Ḋpq H 2 q j , with fourth order tensor A A A defined as

A A A = [H ⊗ 1 1 1 + 1 1 1 ⊗ H]
or A i jkl = (1-ηD H ) 2 δ i j .

  (a) Rectangle rule (b) Trapezium rule

Figure 1 :

 1 Figure 1: Numerical schemes for the time discretization.

-Figure 2 :

 2 Figure 2: Uniaxial monotonic response in quasi-static and dynamics.

Figure 3 :

 3 Figure 3: Evolution of dissipation at one Gauss point in tension .

∆ǫ= 1 . 10 - 6 Figure 4 :

 11064 Figure 4: Evolution of dissipation at one Gauss point in compression .

Figure 5 :

 5 Figure 5: Evolution of the energy density dissipated in tension

  . T-C Initial criterion response corresponds to Section 4.1 model in trace of damage, satisfactory T-C Modified criterion response corresponds to change of tr D into active damage d act in criterion function f .

Figure 6 :

 6 Figure 6: Stress-strain response under an alternated loading tension-compression.

Figure 7 :

 7 Figure 7: Response under alternated loadings in dynamics.

Figure 8 :

 8 Figure 8: Willam test response with the anisotropic damage model.

Figure 9 :

 9 Figure 9: Principle of a dynamics tensile test by spalling, after [53].

Figure 10 :

 10 Figure 10: 3D mesh for computation of the dynamics tensile test.

Figure 11 :Figure 12 :Figure 13 :

 111213 Figure 11: Comparison between the different versions.

  (a) Quasi-statics Brazilian test ,after [51]. (b) Dynamics brazilian test on a drop weight tower.

Figure 14 :

 14 Figure 14: Experimental set-up for a brazilian test.

Figure 15 :

 15 Figure 15: 3D mesh for the quasi-static brazilian test simulation

Figure 17 :

 17 Figure 16: Quasi-static Brazilian test

Figure 18 :Figure 19 :

 1819 Figure 18: 3D Mesh DYN

Table 1

 1 summarizes the different methods and schemes and their designation.

	Method		Rectangle Trapezium
	Direct:	Exact	DER	DET
	Direct:	Approximate	DAR	DAT
	Indirect		IR	IT

Table 1 :

 1 Summary of designations for the different methods and schemes studied.

  for heat transfer.

	Parameter	Value
	E	42 GPa
	ν	0.2
	κ 0	5.10 -5
	A	5000
	a	2.93 10 -4
	b	1
	Ḋ∞	50000 s -1

Table 2 :

 2 Material parameters for modular damage anisotropic model.

Table 3 :

 3 Different versions of the anisotropic damage model .

		Designation
	Initial	INI
	Initial with damage deactivation	INI-d act
	Dynamics without damage deactivation	DYN
	Dynamics with damage deactivation	DYN-d act

  1 2 (H ik δ jl + H jl δ ik + H il δ jk + H jk δ il ). Finally ∂H ∂D = A A A -1 : [H 2 ⊗ H 2 ]

	so that Y i j = 1+ν E σ D kp H pq σ D ql A -1 klmn H 2 mi H 2 jn + η(1-2ν) 18E	tr σ 2 +

Acknowledgments

This work has been funded by the CEA. The authors would like to thank D. Guilbaud for hosting this research and his valuable comments.

tension and 25000 J/m 3 in compression in quasi-static. The role of strain rate effect introduced through the delaydamage law on the dissipation due to damage is very important from a strain rate of 100.s -1 since the dissipation at failure is multiplied by 71 in that case. Comparison with experimental results, which has not be done yet, is crucial in that case. This should bring the future work towards a better identification of the delay-damage law, and in particular, the post-peak response in dynamics. Indeed, no experiments enable to identify directly the post-peak part of the stress-strain response in dynamics, intrinsic dissipation would be an interesting alternative for identification.

We have also numerically checked that the model described here leads to a positive dissipation in the case of a nonproportional loading ( the case proposed by Willam et al). Simulations on 3D structures, in quasi-statics and dynamics, enable to check, once again, that dissipation remain positive at every Gauss points of the structure. Intrinsic dissipation maps can be drawn in order to see the most dissipative areas. We can also estimate the total energy dissipated and compare it to the injected energy.

The comparison with experimental results for a dynamic Brazilian test shows that the total energy dissipated is too weak, which can be explained by the fact that no permanent strain are modeled and that the contact conditions are non-realistic.

Appendix -Calculation of (tensorial) thermodynamics force Y

The associated variable with damage Y is a second order tensor in present framework for anisotropic damage. It is derived from the free enthalpy,

as ∂D H ∂D = 1 3 1 1 1 and where ∂H ∂D can be calculated as follows. First from the definition for H,

where tensorial product ⊗ for symmetric tensors a, b, means