Guillaume Hervé
email: herve@lmt.ens-cachan.fr

Fabrice Gatuingt

Adnan Ibrahimbegovic

On numerical implementation of a coupled rate dependent damage-plasticity constitutive model for concrete in application to high-rate dynamics

Keywords: Concrete constitutive model, damage, plasticity, high-rate dynamic, impact

. We developed an explicit-implicit integration scheme, implemented it in LS-DYNA finite element code, with explicit part dealing with large global problem of equations of motion and implicit part treating local evolution equations. We also provide a detailed consideration of the numerical stability of this kind of scheme for ratedependent damage model. Several numerical tests, both simple ones and a complex problem of the large aircraft impact on a large concrete slab demonstrate the efficiency of the proposed numerical implementation. Comparison between simulations of impact of equivalent aircraft engine missiles and the tests carried out in Sandia laboratory [25] also demonstrates its efficiency.

Introduction

Transient high rate dynamic behavior of concrete is a very important to take into account for design concrete structures in the case of dynamic loading conditions, such as an impact on the structure. This impact loading can be due to explosion, mind blast or an accidental collision of cars, trains or airplanes with the structure. In particular for structures that involve public safety, they have to be design to resist not only the static loading but also the dynamic loading produced extreme conditions. In order to perform a 3D non-linear analysis with a wide variety of damage mechanism in concrete we need a model sufficiently robust and capable of providing a reliable representations for different loading path. But we also need to have a numerical implementation of this model firmly under control in order to ensure a robust computation. In fact, the complexity of any practically interesting model is such that a fully implicit schemes are excluded. In other words, for the numerical simulations of fast transient dynamic, such as an impact on concrete structures, one uses an explicit time integration scheme. The computer program architecture can thus be simplified accordingly, in order to provide a solution for very complex industrial application; with the code as LS-DYNA as one of the most prominent examples, which is also used herein.

The outline of the paper is as follow. In the next section we briefly review the salient feature on the coupled damage plasticity model for high rate dynamics of concrete, adapted to our cases of interest. Then the pertinent details of the chosen scheme for the numerical implementation is presented. In Section 4 we briefly discuss the proposed procedure for model parameters identification. Then the stability of the numerical scheme for a rate dependent damage model is discuss in Section 5. Section 6 provides several illustrative example to show the model performance in representing different inelastic modes of damage in dynamics. Closing remarks and conclusions are given in Section 6.

Coupled rate dependant damage-plasticity constitutive model for concrete

This model is first proposed by Gatuingt [START_REF] Gatuingt | Coupled damage and plasticity modelling in transient dynamic analysis of concrete[END_REF] to simulate explosion in contact or impact of hard projectiles at velocity less than 350 m/s. For this kind of loading three failure mechanisms have to be described. First one which is observed under the impact is a decrease of the material porosity. This phenomenon is represented with homogenization technics by considering concrete as a matrix (cement paste and aggregates) with pores. To model the penetration of a hard striker introducing a large deviatoric strain, a plastic model based on modified Gurson's yield function (to take into account porosity evolution) is used. The third and final mechanism is supplied to handle the case where the compressive wave can reflect on a free surface producing a traction state of stress which is represented with a damage model.

These mechanical effects are combined in the relationships which relate the stresses to the elastic strains:

σ ij = (1 -D) Kε e kk δ ij + 2G(ε e ij - 1 3 ε e kk δ ij ) (1)
where the shear G and bulk moduli K of the coupled model are defined by Mori-Tanaka's expressions:

K = 4K M G M (1 -f * ) 4G M + 3K M f * ; G = G M (1 -f * ) 1 + 6K M +12G M 9K M +8G M f * (2)
with K M and G M respectively the bulk and shear moduli of the matrix material without pores . In the case of a smooth impact (such as impact by cars or airplanes), we can consider that the porosity will not decrease enough under the projectile to induce a variation of the concrete moduli. For the study of these kind of problems we will consider in this paper that K and G remain constant.

The phenomenon of microcracking, in uniaxial tension and compression, is captured with a rate dependent damage model [START_REF] Dube | A Rate Dependent Damage Model for Concrete in Dynamics[END_REF]. Accounting for rate effects is necessary in order to represent the type of response one finds in dynamic experiments, mostly dynamic tensile tests [START_REF] Klepaczko | An experimental method for dynamic tensile testing of concrete by spalling[END_REF]). In addition, the added benefit of the rate dependency is to preserve the well-posedness of the equations of motion when strain softening occurs [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF], [START_REF] Sluys | Wave Propagation, Localisation and Dispersion in Softening Solids[END_REF]. The extensive program of experiments carried out within the French research network GEO showed that there was a marked dependence between the loading rate and the curve relating the volumetric strain to the hydrostatic stress [START_REF] Gatuingt | Prévision de la rupture des ouvrages en béton sollicités en dynamique rapide[END_REF]. The latter is thus captured by implementing a viscoplasticity model. Within the classical framework of small strain cinematics, we use the basic assumption of additive strain decomposition:

εij = εe ij + εvp ij ( 3
)
where εij is the total strain rate, εe ij , the elastic one and εvp ij the viscoplastic strain rate.

The viscoplastic strain rates are obtained following Perzyna's approach:

εvp ij = λ ∂F N T ∂σ ij (4)
The viscoplastic multiplier λ is defined with the power law which also takes into account the porosity f * :

λ = f * (1 -f * ) F N T m vp nvp (5)
where m vp and n vp are material parameters. The porosity evolution is controlled by the irreversible volumetric strain only according to:

df * = k(1 -f * )f * dε p kk (6)
where parameter k is introduced in equation ( 6) in order to be able to calibrate the velocity with which voids are closed.

In equation above, F N T is the modified Gurson's yield function proposed by Needleman and Tvergaard [START_REF] Needleman | An analysis of ductile rupture in notched bars[END_REF]:

F N T σ ij , σ M , f * := 3J 2 σ 2 M + 2q 1 f * cosh q 2 I 1 2σ M -1 + (q 3 f * ) 2 = 0 (7)
where I 1 = T r(σ), the first invariant of stress tensor, J 2 = dev(σ) , the second invariant of the deviatoric part, σ M is the stress in concrete matrix without voids and q 1 , q 2 , q 3 are scalars parameters. The model of this kind treats the concrete as a porous material. This porosity has a great importance on the material behavior when the hydrostatic stress contribution is not negligible. Indeed, this model improves upon the Drucker-Prager [START_REF] Drucker | Soil Mechanics and Plastic Analysis of Limit Design[END_REF] yield surface which is often used for concrete with consequences that the material remains elastic for triaxial compression, which is in contradiction with the compaction experimentally observed [START_REF] Burlion | Experimental analysis of compaction of concrete and mortar[END_REF]. The main interest in the modified Gurson's yield function is to be closed for a hydrostatic state of stress and to provide a kind of CAP model [START_REF] Gatuingt | Gurson's plasticity coupled to damage as a CAP Model for concrete compaction in dynamics[END_REF], [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] or [START_REF] Ibrahimbegovic | Constitutive model of coupled damage-plasticity and its finite element implementation[END_REF].

The constitutive response in tension and compression is controlled by the damage evolution law governed by a rate dependent model given by :

Ḋ =   < εe -ε D0 -1 a D 1-D (1/b) > m D   n D (8)
with the damage criterion of Mazars [START_REF] Mazars | A description of micro and macroscale damage of concrete structures[END_REF] f = εe -κ

where εe = i ε e i + 2 is the elastic equivalent strain for quasi-brittle materials [START_REF] Mazars | A description of micro and macroscale damage of concrete structures[END_REF]. In Equation ( 8) above m D , n D are material parameters which control the rate effect, a, b are material parameters which govern the growth of damage in quasistatic tension and compression and ε D0 is the initial value of damage threshold.

Figure 1 shows the coupled plasticity-damage model response for the hydrostatic stress state for a loading/unloading cycle. We can see that the model captures both strain hardening effect with concrete compaction and the irreversible plastic residual strain upon unloading. Given the complexity of the constitutive relations and the main application domain which pertains to concrete in fast transient dynamics, the model has been implemented in the explicit finite element code LS-DYNA [START_REF] Hallquist | LS-DYNA 3D -Theoretical Manual[END_REF]. The objective of the numerical implementation of the model is to be able to calculate the new state of stress at time t + ∆t knowing the increment of strain ∆ε and the state of strain at time t. In this computational process, the evolution of the damage and viscoplastic strain are totally decoupled in order to preserve the computational efficiency. We refer to Figure 2 for the flow-chart of the implemented computational procedure, which provides the summary of all different cases.

In a first time, we assume that all the strain increment are elastic and we compute the equivalent elastic strain:

ε e ij (t + ∆t) = ε ij (t + ∆t) -ε vp ij (t) εe (t + ∆t) = i ε e i (t + ∆t) + 2 (10)
Furthermore, we assume that all internal variables remain fixed with their rates ( ḟ * , σM and Ḋ) equal to zero, so that we can compute the corresponding stress rate in the elastic predictor:

ε ij (t+∆t) = ε ij (t) + dε ij ε ij (t+∆t) = ε ij (t+∆t) -ε ij (t) 0 ~> e ε Case 1 0 < NT F 0 ~D e ε ε < 0 < NT F vp e 0 ~D ε ε < e F.E. imput Elastic equivalent strain σ ij (t+∆t)= ∆t e Elastic predictor Case 2 Case 4 Case 3 Case 3 Case 1 YES YES NO NO NO NO YES YES YES NO σ ij (t) + σ ij . * Fig. 2. Numerical integration scheme εij (t + ∆t) = εe ij (t + ∆t) + εvp ij (t + ∆t) = εe ij (t + ∆t) σ = (1 -D(t)) K εkk (t + ∆t)δ ij + 2G( εij (t + ∆t) -1 3 εkk (t + ∆t)δ ij ) εvp ij (t + ∆t) = 0; ḟ * (t + ∆t) = 0; σM (t + ∆t) = 0 (11)
In order to verify if the correction is needed to obtain the real state of stress, we test the positivity of the elastic equivalent strain predicted. If this strain is negative, we deal with a loading path mainly in triaxial compression and the damage variable D will never evolve. We are thus concerned with either the case 1 (elasticity) or with the case 3 (only plasticity) as shown in Figure 2. In case 2 we only have a damage evolution with no plasticity. For the "coupled" case 4, we have evolutions of the damage variable and of the plastic strain. In this case, we first compute the new state of damage and we use it to compute the plastic correction without any subsequent iteration.

For the plastic correction we use the return mapping algorithm [START_REF] Ortiz | An analysis of a new class of integration algoriths for elastoplastic constitutive relations[END_REF] for the plastic part of the model assuming that the damage D is known at the state t + ∆t. The latter is computed first from the evolution of the damage during the incremental process, by using an explicit scheme due to the formulation of the damage growth (see equation 8). We can then write:

εij (t + ∆t) = εe ij (t + ∆t) + εvp ij (t + ∆t) = 0 =⇒ εe ij (t + ∆t) = -εvp ij (t + ∆t) σij = -(1 -D(t + ∆t)) K εvp kk δ ij + 2G( εvp ij - 1 3 εvp kk δ ij ) εvp ij = λ∂F N T ∂σ ij ḟ * = k(1 -f * )f * εvp kk = λk(1 -f * )f * ∂F N T ∂σ kk σM = λ E t (1 -f * )σ M ∂F N T ∂σ ij σ ij (12)
Substituting equation ( 12) 3 into equation ( 12) 2 we obtain:

σij = -λ(1 -D(t + ∆t)) K ∂F N T ∂σ kk δ ij + 2G( ∂F N T ∂σ ij - 1 3
∂F N T ∂σ kk δ ij ) ( 13
)
with λ = f * (1-f * ) F N T
mvp (see equation ( 2)) for linear viscoplasticity (n vp = 1) and where plastic loading (F N T > 0) is implicitly assumed.

The rate of change of the yield function during the relaxation process can be written as follows (see Figure 3):

ḞNT = ∂F N T ∂σ ij : σij + ∂F N T ∂f * . ḟ * + ∂F N T ∂σ M . σ M (14)
Substituting σij , ḟ * and σM with their evolutions in equation ( 12), defined as a function of the plastic multiplier λ, we obtain an ordinary differential equation to solve:

ḞNT = - F N T m vp f * 1 -f *    ∂F N T ∂σ ij :(1 -D) K ∂F N T ∂σ kk δ ij + 2G( ∂F N T ∂σ ij - 1 3
∂F N T ∂σ kk δ ij ) + ∂F N T ∂f * . k(1 -f * )f * ∂F N T ∂σ kk + ∂F N T ∂σ M . E t (1 -f * )σ M ∂F N T ∂σ ij σ ij    ( 15
)
This kind of problem has been solved by Ortiz and Simo [START_REF] Ortiz | An analysis of a new class of integration algoriths for elastoplastic constitutive relations[END_REF] introducing the instantaneous relaxation time t:

t =m vp 1 -f * f *    ∂F N T ∂σ :(1 -D) K ∂F N T ∂σ kk δ ij + 2G( ∂F N T ∂σ ij - 1 3
∂F N T ∂σ kk δ ij ) + ∂F N T ∂f * . k(1 -f * )f * ∂F N T ∂σ kk + ∂F N T ∂σ M . E t (1 -f * )σ M ∂F N T ∂σ ij σ ij    (16)
equation ( 15) reduces to :

ḞNT = - F N T t =⇒ ln(F N T ) = -t t ( 17
)
The following algorithm is then applied: the elastic predictor is computed. The return path can then reach a suitably updated yield surface by means of a sequence of straight segment ∆ i t (see figure 3) which is the instantaneous relaxation time of the linear differential equation.

∆ i t =m vp 1 -f * (t i ) f * (t i )    ∂F N T i ∂σ ij :(1 -D) K ∂F N T i ∂σ kk δ ij + 2G( ∂F N T i ∂σ ij - 1 3
∂F N T i ∂σ kk δ ij ) - ∂F N T i ∂f * • k(1 -f * (t i ))f * (t i ) ∂F N T i ∂σ kk - ∂F N T i ∂σ M • E t (t i ) (1 -f * (t i ))σ M (t i ) ∂F N T i ∂σ ij : σ ij (t i )    (18)
To update the variable σ ij , f * and σ M during the return mapping iterative process, we use:

σ ij (t i+1 ) = σ ij (t i ) + σij (t i )∆ i t σ M (t i+1 ) = σ M (t i ) + σM (t i )∆ i t f * (t i+1 ) = f * (t i ) + ḟ * (t i )∆ i t (19)
We can thus compute F N T i+1 as well as the new values of the internal variables. The total relaxation time is obtained with the following equation:

t i+1 = t i + ∆ i t.log F N T i F N T i+1 (20)
We obtain the exact plastic correction for F N T = 0 which is obtained when the total relaxing time is equal to the real time increment ∆t [START_REF] Ortiz | An analysis of a new class of integration algoriths for elastoplastic constitutive relations[END_REF].

During this internal variable computation and in particular for the elastic prediction, we can obtain numerical problem to compute the modified Gurson's yield function (see eq. ( 7)). This is due to the high numerical values obtained with the cosh function when the term q 2 I 1 2σ M becomes large enough. This phenomenon is represented in Figure 4. We follow Mahnken [START_REF] Mahnken | Aspects on the finite-element implementation of the Gurson model including parameter identification[END_REF] who proposed to modify the expression of the hyperbolic cosine with the power series development around a critical point. We choose for critical value q 2 I 1 2σ M = X c = 30 and we obtain the expression:

q 2 I 1 2σ M ≤ X c : F N T cosh q 2 I 1 2σ M q 2 I 1 2σ M > X c : F N T cosh(X c ) + sinh(X c )( q 2 I 1 2σ M -X c ) + 1 2 cosh(X c )( q 2 I 1 2σ M -X c ) 2 (21)
This method ensures continuity under the critical point and limits the growth up of this term when X > X c . Figure 5 shows the response of the model for a hydrostatic compression (where only the plastic criterion is activated) for different strain increment but with the same strain rate. We can notice that the model response is now totally independent of the strain increment. This is due to the implicit integration scheme used for the plasticity.

f (q 2 I 1 / 2 σ M ) q 2 I 1 / 2 σ M 0 2 4

Parameters identification

The main goal of this section is to give a systematic procedure in order to calibrate and choose the model parameters. We saw in Section 2 that the constitutive equations are governing with two inelastic threshold functions. The first one is based on Mazars's work and concerns the damage evolution. The second one is the modified Gurson's yield function for the evolution of the plastic strain. It is important to note that, the Mazars's threshold function is activated in traction before the Gurson's yield function, means that only the damage constitutive law is used see Figure 6. The representation of the yield functions was generated with MatSGen software developed at LMT-Cachan by François [START_REF] François | MatSGen : user guide[END_REF]. So if we want to calibrate the model damage parameters we have to perform a traction test in statics and in dynamics. In addition, for a triaxial compression load path, only the Gurson's yield function is activated which means that we will use only the viscoplastic constitutive law. To calibrate the model parameters for the viscoplastic strain evolution, we have to perform the hydrostatic compression test in statics and in dynamics. On the other hand, if we are in a load path in simple compression with low confinement, we activate the two thresholds and we will have a coupled response of the constitutive equation. Nevertheless, if we choose appropriate parameters, we can obtain a simple compression depending only on the damage law which will permit to calibrate the damage parameters in compression.

In table 1 we present a summary of the parameter values for MB50.

Stability of the numerical implementation

Rate dependant damage model

An homogeneous material element at a macroscopic point of view deforms in a homogeneous manner if a homogeneous stress is applied at its boundaries. But when the strain becomes larger, due to the loading for example, concentration can occur in one element which lead to localization of strain over more or less extended area and the deformation of the considered element ceased to be homogeneous.

q 1 1.5 - ǫ D 0 1.10 -4 - q 2 0.8 - a t 20000 - q 3 1 - b t 1.6 - σ M 0 60 MPa n Dt 5 - n 15 - m Dt 0.510 -4 - k 45 - a c 3000 - m vp 1.110 -2 - n Dc 20 - n vp 1 - m Dc 0.510 -3 -
A good way for describing localization in terms of continuum theory is the strain rate discontinuity [START_REF] Rudnicki | Conditions for the localization of deformations in pressure-sensitive dilatant materials[END_REF][START_REF] Rice | The localization of plastic deformation[END_REF][START_REF] Rice | A note on somes features of the theory of localization of deformations[END_REF]. The localization implies a non uniqueness in the incremental response of a homogeneous, homogeneously strained body and also implies a vanishing speed of acceleration waves [START_REF] Hadamar | Lecons sur la Propagation des Ondes et les Equations de l'Hydrodynamique[END_REF][START_REF] Hill | Acceleration waves in solids[END_REF].

The equation of the damage model is:

σ = (1 -D)Eε (22)
with the evolution law of the damage governed by a rate dependent model (see equation ( 8)) considering n D = 1. This is an evolution problem which is solved by a time integration scheme. This latter provides a discretized solution of this evolution problem, which can be written as :

div σ n+1 = 0 σ n+1 = (1 -D n+1 )Eε n+1 = G(ε n+1 , D n+1 ) ∆D = D n+1 -D n = ∆t Ḋn+1 = F (ε n+1 , D n+1 ) (23)
This is a non-linear problem which can be solved by the Newton scheme. We have to ensure that this remains a well-posed problem by verifying the positiveness of the consistent tangent operator.

Therefore, we linearize the non-linear evolution laws:

δσ n+1 = ∂G ∂ε n+1 : δε n+1 + ∂G ∂D n+1 δD n+1 δD n+1 = ∂F ∂ε n+1 : δε n+1 + ∂F ∂D n+1 δD n+1 (24)
With the equation ( 24) 2 we have:

δD n+1 = ∂F ∂ε n+1 : δε n+1 1 -∂F ∂D n+1 (25)
Using this expression of δD n+1 in equation ( 24) 1 we can write:

δσ n+1 = H : δε n+1 (26)
with

H = ∂G ∂ε n+1 + ∂G ∂D n+1 ⊗ ∂F ∂ε n+1 1 -∂F ∂D n+1 (27)
the tangent modulus of the non-linear evolution law. In our case, we have:

∂G ∂ε n+1 = (1 -D n+1 )E ∂G ∂D n+1 = -Eε n+1 ∂F ∂ε n+1 = ∆t m D < ε n+1 > + εn+1 ∂F ∂D n+1 = - ∆t m D . 1 a D n+1 1-D n+1 (1/b) bD n+1 (1 -D n+1 ) = -α ∆t m D (28)
Which leads to :

H = (1 -D n+1 )E + -Eε n+1 ⊗ ∆t m D <ε n+1 > + εn+1 1 + α ∆t m D (29)
In order to solve the problem, we have to verify the equilibrium for the linearized problem: div δσ n+1 = 0 (30) which is a well-posed problem under the condition:

det nH n = 0 (31)
With the tangent modulus obtained in equation ( 29) we obtain from the last expression :

det (1 -D n+1 ) nEn - ∆t m D .ε n+1 .(1 + α ∆t m D ) n(Eε n+1 ⊗ < ε n+1 > + )n = 0
(32) this is verified for [START_REF] Bigoni | Uniqueness and localization-I. Associative and non-associative elastoplasticity[END_REF]:

m D .ε n+1 .(1 + α ∆t m D ) ∆t = (< ε n+1 > + n){nEn(1 -D n+1 )} -1 (Eε n+1 ) (33)

Application to 1D bar in traction

For a 1D bar in traction we can obtain a simple form of the tangent modulus: For this simple case, the stability condition is H = 0. With this condition, we can obtain a stability time step ∆t:

H = (1 -D n+1 )E -Eε n+1 ∆t/m D 1 + α ∆t m D (34)
∆t = m d ε n+1 /(1 -D n+1 ) -α (35)
Figure 7 shows the reduction of the critical time step with the damage increase. This curve has been obtain for model parameters identify for concrete MB50, with the values of parameters given in table 1. We can observe that for these values, the critical time step is big enough to ensure stability of the numerical scheme even at the end of the damage evolution.

6 Numerical simulations

Patch Test

In order to see if the model was implemented correctly in the FE code LS-DYNA, we decided to test two different meshes (one regular and one distorted see Figure 8), under a homogeneous stress field -the classical path test (see [START_REF] Zienkiewicz | The Finite Element Method -Volume 1 : The Basis[END_REF] carried out in dynamics.

First we carried out traction tests on a regular mesh cube and a distorted mesh one in order to check if the damage phenomena exhibit mesh dependency, see The stress vs strain diagram clearly indicates no mesh dependency as shown in Figure 10, with two curves which remain in very good accordance.

Then we carried out the same kind of tests for an hydrostatic state of stress as illustrated in Figure 11, with the same purpose to verify the good accordance of the model response for different meshes. Figure 12 showing the hydrostatic pressure vs voluminal strain curves for regular and distorted mesh indeed confirms good accordance in both cases.

Reinforced concrete slab

Another test of the model capabilities was to carry out a large scale computation. As it would be dedicated to simulations of soft impact phenomena and especially aircraft impact, we carried out the computation of impact between 747 class Boeing like aircraft and a reinforced concrete slab.

The slab was a parallelepiped of 80 by 80 meters large and 80 centimeters The aircraft impacted the slab following a 20 degrees angle, at 252 m/s this means 900 kph the ultimate speed of this kind of aircraft in such a configuration, as it is shown in figure 16.

The finite element model of the aircraft was constructed by a mesh of Belytschko-Lin-Tsay shell elements, the slab is constituted of under-integrated (one Gauss point) brick elements . The automatic contact surface to surface LS-DYNA option for contact was used in this simulation.

The reinforcement was represented with beams and a perfect contact between concrete and steel bar was assumed, by merging nodes between beams elements and 3D elements. Reinforcement ratio was 0.4 % in each direction. The results presented in figure 17 show the evolutions of the damage on the impacted face and the opposite one. When damage reaches the maximum value of 1 the concrete is locally fully destroyed. case that a huge crater has occurred on the top face of the slab, yet we can assume that the scabbing on the rear face is limited to a 15x15 meters square area containing nearly 7 disseminated scabbed discs. According to the 0.25 damage value on the half of the opposite face, we can assume that half of rear face is considerably weakened with the presence of several cracks.

Sandia's Laboratory tests simulations

In order to check the relevance of the constitutive relations we carried out simulations of impact tests performed in Sandia Laboratory by Sugano et al. and presented in [START_REF] Sugano | Local damage to reinforced concrete structures caused by impact of aircraft engine missiles -Part 1. Test program, method and results[END_REF] and [START_REF] Sugano | Local damage to reinforced concrete structures caused by impact of aircraft engine missiles -Part 2. Evaluation of test results[END_REF]. Following the test program, we simulate impact of n aircraft engine equivalent missiles, considering several sizes. In particular we carried out simulations of LED, MED and SER missiles which means, respectively, Large size Equivalent and Deformable missile, Medium size Equivalent and Deformable missile and Small size Equivalent Rigid missile. The details for the selected three missiles are presented in Figure 18.

The missiles impacted three slabs, with a particular choice depending on the kind of missiles. The slabs are presented in figure 19, from left to right with the two slabs placed in a box have been used for small size and large size missile tests. Among several impact tests, carried out by Sugano et al., we chose to simulate impact of those exposed in Table 2.

In our simulations missiles where represented by using the finite element models constructed with Belytschko-Lin-Tsay shell elements, except for SER missile, whose mesh is built with 3D elements. The concrete part of the slab was represented with under-integrated 3D elements, the reinforcement was represented with truss-bars. In each case reinforcement ratio was 0.4% in each direction. The qualitative results of numerical simulations, the same as those obtained in the tests, are presented in Table 3; we note that penetration means that the missile penetrated the slab without having gone through it, perforation means that missile went through, scabbing means that the impact generated a scab on the rear face of the slab.

All computations where stopped when the velocity of the missile stopped decreasing, this means when the missile finally got stuck in the concrete, or when concrete did not bring any resistance to its penetration, since the slab was totally damaged. The damage fringes of the simulations at the end of the computations are exposed in Figures 13, 14 and 15, respectively, for impact tests S10, S28 and L5. In figure 13, we can observe that the damage is not enough continuous to let the missile goes through, the damage fringe on the rear face exposes that a scabbing is occurring.

In Figure 14 the damage fringe is continuous through the slab thickness under the impacted zone, thus we can conclude that a perforation has occurred and a scabbing too.

In Figure 15, we see that the damage is localized in the impact area, a large undamaged volume separates damage on the rear face and damage on the front face. We can here conclude on a penetration with maybe some cracks or a small scab.

Considering these observations and Table 3, we can conclude that obtained a very good accordance with the tests. The latter also confirms the validation of the constitutive equations that had been developed and also their implementation. We have presented a very careful consideration of the numerical implementation aspects for the constitutive model of coupled damage-plasticity, within the framework of explicit-implicit scheme used by computer code LS-DYNA. The proposed integration scheme can be considered as optimal one for the given problem, in which we combine an explicit solution of global equations of motion, with the implicit solution of local evolution equations, which ensures both the computational efficiency and admissible values of stress field during the damage computation. The model is especially suited for representing different damage mechanisms which develop in concrete under high rate loading, both in tension and compression. The former pertains to tension cut-off damage criterion, whereas the latter is represented by the plasticity model with Gurson-like criterion, with a particular feature of accounting for hardening due to concrete compaction.

Fig. 1 .

 1 Fig. 1. Hydrostatic response of the model used 3 Numerical implementation

Fig. 3 .

 3 Fig. 3. Principle of the return mapping

Fig. 4 . 3 Fig. 5 .

 435 Fig. 4. Modification of the cosh function to prevent numerical overflow

Fig. 6 .

 6 Fig. 6. Needleman-Tvergaard and Mazars yield function in stress space

Fig. 7 .

 7 Fig. 7. Critical time step versus damage

Fig. 8 .

 8 Fig. 8. Regular and distorted meshes used for the patch test u(t)

.10 - 4 10 - 5 Fig. 9 .

 459 Fig. 9. Patch test in traction

Figure 9 .

 9 Figure 9. The stress vs strain diagram clearly indicates no mesh dependency as shown in Figure 10, with two curves which remain in very good accordance.

Fig. 10 .

 10 Fig. 10. Traction response for the two different meshes

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Patch test in hydrostatic compression

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Simulations of L5 test -Top, bottom and side views

Fig. 17 .Fig. 18 .Fig. 19 .

 171819 Fig. 17. Damage evolution during the impact

Table 1

 1

	Model parameters for MB50
	Plasticity par. Values Units Damage par. Values Unity
	E 0	5.510 10	Pa	E 0	3.510 10	Pa
	ν	0.2	-	ν	0.2	-

  Thus we can conclude in this

	u(t)	u(t)
		u (m)
	u(t)
		3.10 -4
	u(t)	u(t)
		10 -5	t (s)
	u(t)
	u(t)

Table 2

 2 Simulated tests characteristics No. Missile type Missile velocity Slab thickness Reinf. ratio Slab type

			(m/s)	(m)
	S10	SER	141	0.15	0.4	Small # 1
	S28	SED	196	0.06	0.4	Small # 1
	L5	LED	214	1.60	0.4	Large # 3

Table 3

 3 Tests results for each impact

	No. Perforation	Scabbing	Penetration
	S10	No	Yes	Yes
	S28	Yes	Yes	Yes
	L5	No	Some cracks	Yes
	Fig. 13. Simulations of S10 test -Top, bottom and side views
	Fig. 14. Simulations of S28 test -Top, bottom and side views


E Parizet

F Lerme

Y Padayachi

O Collery

Loudness models in the very low frequency range: application to the prediction of airplane cockpit noise comfort

Aircraft interior noise is a major design stake with respect to airline requirements for passenger and crew comfort. A study dealing with airplane cockpit noise annoyance had shown that loudness computed according to ISO 532-B standard was a good gauge of qualitative evaluation. This was not true for loudness computed according to ANSI S3.4 standard. It seems that this was due to differences of the two models for very low frequency sounds, as these models were adjusted to fit to different iso-loudness curves. In order to check this hypothesis, direct loudness estimations will be presented in the low frequency range (from 35 Hz to 65 Hz). These estimations will be compared to predictions given by the two models.

INTRODUCTION

Loudness is a key factor of unpleasantness and annoyance. Two main models are available to predict loudness of sounds. They are both standardized and can be used by manufacturers for sound quality improvement. In this paper, the use of these models will be illustrated in the case of airplane cockpit noise. Some discrepancy has been observed between the models, especially for one sound with a very low frequency tone. In order to have a more clear understanding of this discrepancy, a second experiment was conducted, during which listeners had to estimate the loudness of very low frequency sounds (from 35 to 65 Hz).

UNPLEASANTNESS EVALUATION

Experimental setup

Six stimuli were provided by Airbus corresponding to in-flight recordings with a microphone located in the cockpit of six aircraft (fs = 48 kHz, 16 bits). One of the stimuli exhibited an abnormally high emerging frequency (63 Hz). An additional stimulus was created by filtering out this frequency; the spectra of the original and modified signals are presented in Figure 1.

Fig. 1 -Examples of inside cockpit spectra of stimuli used in the experiment.

All stimuli had a duration of 5 seconds, including a 150 ms cosine in and out fading. These stimuli were played back by a couple of loudspeakers (Mackie Tapco S8) placed approximately two meters from the subject. The loudspeakers and the listener were placed at the edges of an equilateral triangle. This set up was arranged in a sound proof booth (approx. 3.5 m long and 2.5 m large).

The experiment was a pair comparison one. Each of the 21 possible pairs of different stimuli was played to the listener, plus two for training at the beginning of the experiment. The presentation order was randomized and arranged according to the rules of Ross' series. The listener was asked to compare the unpleasantness of the two sounds and had to give his answer on a continuous scale, which indicated five labels (from "A is much more unpleasant than B" to the reverse, including an equality label). He could play again the pair as often as needed before giving his answer. The whole experiment was conducted by the "Jury Testing" software (developed by 01dB-Metravib and running on a PC outside the booth).

31 subjects participated to the experiment: 15 students and 16 older searchers of the laboratory, 12 women and 19 men. They did not relate any hearing impairment, but their hearing ability was not checked.

Results

First of all, analysis of circular triads showed that all individuals could be considered as giving reliable answers, so that all subjects were considered in the following analysis. Then, a hierarchical cluster analysis of individual answers revealed a good inter-subjects agreement. 10dB Therefore, the results of the 31 listeners were averaged, giving a set of 21 preference probabilities. These probabilities were analyzed using a BTL model, which gave the unpleasantness scores of sounds shown on figure 2. In this figure, the 95% confidence intervals of each score are also presented. They were estimated using a boot-strap technique (500 trials, 20 subjects, randomly selected with replacement). As can be seen in figure 2, filtering out the low-frequency component reduced the unpleantness of sound n°6 (n°7 being the modified signal). 4 groups of sounds can be seen: n° 4 and 5, n° 2, 3 and 7, n°1 and n°6, these groups being ordered according to an increasing unpleasantness.

Metrics

In most cases, unpleasantness is strongly related to loudness (some examples concern noise in trains [START_REF] Parizet | Noise assessment in a high-speed train[END_REF] or in cars 2 ). Therefore, loudness metrics were natural candidates for the description of subjective/qualitative results. Two models were used. They are described by two standards: ISO 532-B 3 and ANSI S3.4 -2005 [START_REF]Procedure for the computation of loudness of steady-state sounds[END_REF] . It appeared that ISO model represented results in an accurate way (figure 3). The left part of the figure shows loudness values of sounds (X-axis) and BTL scores (Y-axis), indicating a high correlation between the two sets of values (R = 0.86, p = 0.013). The right part of figure 3 shows the ratio of the loudness of sounds in each of the 21 pairs presented to listeners (X-axis) and the preference probabilities in pairs (Y-axis). The relation is still clear, the ISO-loudness can be used to predict unpleasantness of sounds in airplane cockpits. On the other hand, ANSI loudness failed to represent unpleasantness (figure 4). Namely, this is mainly due to sound 6, whose loudness is underestimated by this model. Indeed, the relation between sound pressure level and loudness of a 63 Hz pure tone is different in the two models (figure 5). Loudness computed according to the ISO model is approximately twice loudness computed according to the ANSI one. This is certainly due to the fact that these two models have been adjusted according to different equal-loudness level contours. This difference in loudness values can be observed in different versions of the ISO 226 curves [START_REF]Acoustics: normal equal-loudness level contours[END_REF][START_REF]Acoustics: normal equal-loudness-level contours[END_REF] .

LOUDNESS EXPERIMENT

The goal of this experiment was to evaluate the loudness of very low frequency sounds, in order to get a better idea of the validity of the two loudness models in that frequency range. Very few studies have been conducted to establish the relation between sound pressure level and loudness for frequencies lower than 100 Hz (which was the lowest frequency used by Hellman and Zwislocki 7 ). As perception of very low frequency sounds can be extra-auditive [START_REF] Berglund | Sources and effects of low-frequency noise[END_REF] , it was decided to use a loudspeaker, as in the previous experiment.

Experimental setup

Experiments took place in the same sound proof booth. The first resonance frequency of the room is approximately 43 Hz. A subwoofer was placed in a corner of the room. This subwoofer has been designed by Ph. Herzog (Laboratoire de Mécanique et d'Acoustique, Marseille) and built by the Centre de Transfert et de Technologie du Mans (CTTM). It is made of a 30 cm loudspeaker (PHL Audio W384) mounted in a vented wooden box. The first requirement was to place the listener at a position at which the level of the sound did not vary too much with the exact position of the subject. The opposite corner of the room was a good candidate position. A set of 27 microphones was placed in an arrangement made of 3 sets of 9 microphones located at three different heights (1 m, 1.1 m, 1.2 m), which described the range of the head of the sitting subject. In each plane, 8 microphones were arranged in a circle (with a 11 cm radius) and another one located at the center of the circle. This center was 33 cm far from the two walls defining the corner of the room.

A white noise generator was used for this experiment and the frequency response function was measured for each microphone. In figure 6, the median, minimum and maximum values of these functions are shown. As it can be seen in figure 6, the useful frequency range was limited to 30 -65 Hz. Below 30 Hz the level was too low and above 65 Hz the spatial variation was considered as too important.

It was decided to use pure tones as stimuli. In order to reduce the number of experiments, three frequencies were considered: 35, 45 and 65 Hz. In a second step, harmonic distorsion of the subwoofer was measured. For each frequency, the maximum level was selected so that the level of all higher harmonics was below normal hearing threshold. These maximum levels can be seen in table 1. For each frequency, 10 levels were selected, going from slightly above normal hearing threshold to maximum level. These levels were equally distributed between the minimum and maximum values of table 1. A set of 10 signals was prepared, with a duration of 3 seconds each, including a 100 ms in and out fading.

Procedure

The procedure consists in absolute magnitude estimation, as it has been shown by Boullet 9 that it gives reliable results with a limited duration of the experiment. When hearing a sound, the listener was asked to give a number representing loudness. Three experimental sessions were conducted (one for each frequency). In each session, each stimulus was presented four times to the subjects, so that the session was arranged in four blocks of 10 stimuli. In each block, the presentation order was randomized. The sessions order (there were 6 possible orders) was balanced within listeners and 30 subjects participated to the experiment.

Results

The experiment has just been conducted and its results will be presented at the conference. Subjective loudness evaluations will be compared to the predictions given by the two models.

CONCLUSIONS

These two experiments give a deeper knowledge of unpleasantness of low frequency sounds. More precisely, 2 questions shall be addressed: (1) was the unpleasantness of sound 6, as determined in experiment 1, related to loudness only, or did the tonal character of that sound contribute to unpleasantness ? (2) Which of the two standards should be recommended for the prediction of loudness of such sounds? We hope answering from the analysis of the related experiments.

Fig. 2 -

 2 Fig. 2 -Merit scores of the 7 stimuli.

Fig. 3 -

 3 Fig. 3 -Relation between ISO-loudness values of sounds and subjective unpleasantness. Left:BTL scores; right: preference probabilities.

Fig. 4 -

 4 Fig. 4 -Relation between ANSI-loudness values of sounds and unpleasantness.

Fig. 5 -

 5 Fig. 5 -Loudness of a 63 Hz pure tone, as predicted by the two standardized models.

Fig. 6 -

 6 Fig. 6 -Median (black curve), minimum et maximum (red curves) of the frequency responsefunctions measured in the room.

Table 1 -

 1 Sound pressure level range for each tone

	Frequency Max. level Max. level
	35 Hz	63 dB	85 dB
	45 Hz	60 dB	85 dB
	65 Hz	47 dB	83.5 dB


Adnan Ibrahimbegović

Damijan Markovič

Fabrice Gatuingt

Constitutive model of coupled damage-plasticity and its finite element implementation

Keywords: plasticité-endommagement couplé, implantation éléments finis coupled plasticity-damage, finite element implementation

Dans ce travail on présente un cadre théorique général du développement des modèles de couplage de deux types de comportement anélastique, la plasticité et l'endommagement. On introduit la nouveauté principale par rapport aux modèles précédents de ce type en utilisant un critère pour définir le domaine élastique valable aussi bien pour la plasticité que pour l'endommagement, qui peut être adopté pour une grande variété des matériaux, d'une part pour les metaux poreux et d'autre part pour le béton en compaction. L'implantation numérique est d'abord présentée pour un cas unidimensionnel simple et ensuite généralisée pour les critères 2D et 3D, pertinents à des métaux ou des bétons.

ABSTRACT. In this work we present a general theoretical framework for developing a constitutive model capable of coupling two basic types of inelastic behaviour, plasticity and damage. We elaborate upon the main novelty with respect to the previous models of this type, which pertains to a systematic use of criteria for defining the elastic domain, both for plasticity and damage, which can be adapted to a very wide variety of engineering materials, from metals with voids on one side to concrete compaction on the other side. The numerical implementation is first presented for a simple one-dimensional case, and subsequently extended to 2D and 3D criteria which are adequate for either metals or concrete.

Introduction

Two classes of constitutive models of inelastic behaviour, those of plasticity and damage, are most frequently employed for engineering materials. Plasticity model (e.g. see [HIL 50], [LUB 90] or [SIM 98] for a comprehensive account), with a clearly defined yield criterion for identifying the occurrence of plastic deformation which does not affect the elastic response which remains the same in elastic loading and unloading, is mostly applicable to simplified representation of metals and alloys. Similarly, the continuum damage model (e.g. see [LEM 88] or [KRA 96]), with its ability to account for the micro-cracking induced change of elastic response but without residual deformation upon loading, is only a simplified representation of cracking of concrete or soil-mechanics materials. It has been long recognized (e.g. see [LEM 85]) that the validity of these basic constitutive models can be significantly extended when combined into a single coupled model of damage and plasticity. The applications which can be tackled in such a manner are typically outside of the scope of each of the basic models when acting alone, ranging from the metals with voids (e.g. see [GUR 77], [NEE 84]) to compacting concrete (e.g. see [ GAT 99]). Typical attempts of merging these two constitutive models have been developed along a one-track approach (e.g. see [LEM 85], [BEN 88], [SIM 87] or [JU 89], among others) where the damage model is first applied to produce the equivalent stress measure for virgin material or so called effective stress, followed by the plasticity model defining the yield criterion in terms of this effective stress. Even nowadays, the same kind of developments are being carried out to provide a reliable representation of constitutive behavior of concrete (e.g. see [MES 98], [JOH 99] or [MAH 00], among others), each proposing a rather complex computational procedure to obtain stresses and tangent moduli.

In this work we depart from these previous developments to propose a new kind of coupled damage-plasticity model which has the following desirable features: (i) both plasticity and damage model are built around the corresponding criteria indicating the presence of inelastic deformation. For the former this amounts to constructing the classical yield criterion only in terms of total rather than the effective stress, whereas for the latter this implies using a somewhat less standard continuum damage formulation with a clearly defined damage criterion. (ii) the numerical computation on both plasticity and damage part of the model can be carried out in parallel (even making use of the presently available parallel computer architecture) and their final results can then be merged through a local iterative procedure which imposes the uniqueness of stresses.

The outline of the paper is as follows. In the next section we briefly present the main idea of coupling the plasticity and damage in a single constitutive model, starting with a simple one-dimensional setting and extending it to a more general case. Numerical implementation is discussed in Section 3 and several numerical examples are presented in Section 4. In Section 5 we state some closing remarks.

Theoretical formulation

For clarity, we first start with a simple one-dimensional setting allowing for only non-trivial component of the stress tensor and the strain tensor. It is mostly a matter of tensor calculus to develop the corresponding 2D or 3D version of any equation governing the model. However, the particular forms of the yield or damage criterion provide a very significant wealth of possible applications in higher dimensional case, making the model useful for both metals and soil-mechanics materials.

Basic hypothesis

The one-dimensional constitutive model of coupling the damage and plasticity can be built on 3 basic hypotheses: additive decomposition of the total strain field, the strain energy and finally the plasticity and damage criteria. The first of them implies that the total deformation can be additively decomposed into elastic part ǫ e , plastic part ǫ p and damage part ǫ d , leading to

ǫ = ǫ e + ǫ p + ǫ d . [1]
We note in passing that the damage models rarely make use of the notion of the damage strain, with few exceptions (e.g. see [CAR 97] or [YAZ 90]) dealing only with rate form of the equation (1) .

The second ingredient of the model governing elastic response is specified in terms of strain energy. Assuming the simplest quadratic form in terms of the corresponding state variables we can write the strain energy as the sum of elastic, damage and plastic parts:

Ψ(ǫ, ǫ d , D, ξ d , ǫ p , ξ p ) = Ψ e (ǫ e ) + Ψ d (ǫ d , D) + Ξ p (ξ p ) + Ξ d (ξ d ), [2]
where,

Ψ e (ǫ e ) = 1 2 ǫ e C e ǫ e [3]
In equations (2) and (3) above, the elastic modulus is denoted with C e , ǫ p and D are internal variables of plastic strain and damage compliance, ξ p and ξ d are strain like internal variables describing the hardening phenomenon for plasticity and damage, respectively, and Ξ p (ξ p ) and Ξ d (ξ d ) are the corresponding hardening functions. Instead of working with the damage strain energy Ψ d (•), it is convenient to appeal to the Legendre transformation (e.g. see [STR 86]) to define complementary energy, which is postulated as a quadratic form in stress with

χ d (σ, D) := σǫ d -Ψ d (ǫ d , D) = 1 2 σDσ , [4]
The final group of basic model ingredients is provided to specify the elastic domain, where no change of internal variables takes place, along with the yield criteria Φ p and damage criteria Φ d with

Φ p (σ, q p ) = |σ| -(σ p y -q p ) ≤ 0 Φ d (σ, q d ) = |σ| -(σ d f -q d ) ≤ 0 , [5]
where q p and q d are stress like variables describing the hardening phenomena, σ y is the yield limit and σ f is the fracture limit.

We next show that all the remaining ingredients can be derived from these three basic ones simply by the standard thermodynamics consideration and the principles of maximum plastic and maximum damage dissipation. To start, one can provide the local form of the second principle of thermodynamics (e.g. see [LUB 90], [LEM 88]), stating that the total inelastic dissipation is always non-negative. Subsequently, by making use of the results in (2) to (4) we can further write

0 ≤ Ḋ = σ ǫ - Ψ = (σ - ∂Ψ e ∂ǫ e ) ǫe + Ḋp σ ǫp - ∂Ξ p ∂ξ p ξp + σ( ∂χ d ∂σ -ǫ d ) + ∂χ d ∂D Ḋ - ∂Ξ d ∂ξ d ξd Ḋd , [6]
where D p and D d denote, respectively, plastic and damage dissipation.

In an elastic process with no change of plastic variables, with ǫp = 0 and ξp = 0, and no change of damage variables, with Ḋ = 0 and ξd = 0, which further implies that no dissipation takes place, with Ḋp = 0, Ḋd = 0 and Ḋ = 0, it follows from (6) above that the stress can be obtained from the elastic strain energy:

σ := ∂Ψ e ∂ǫ e = C e ǫ e [7]
and that the damage strain is defined through stress and the current value of he damage compliance:

ǫ d := ∂χ d ∂σ = Dσ. [8]
Assuming that the stress constitutive equation in (7) and damage strain definition in (8) remain the same in an inelastic process we can conclude from dissipation inequality in (6) that

0 < Ḋ = σ ǫp + q p ξp Ḋp + 1 2 σ Ḋσ + q d ξd Ḋd , [9]
where we introduced stress-like variables q p and q d which are as thermodynamically conjugate to hardening variables ξ p and ξ d and which can be computed from

q p = - ∂Ξ p ∂ξ p q d = - ∂Ξ d ∂ξ d .
[10]

Plasticity model

In an inelastic process where the plastic module is activated we can appeal to the principle of maximum plastic dissipation (e.g. see [HIL 50] or [LUB 90]) to conclude that among all plastically admissible fields of stress and stress-like hardening variable we can choose those which render the maximum of the plastic dissipation. This can be formally presented as a constrained minimization problem which can also be recast in a non-constrained form by making use of the Lagrange multiplier procedure (e.g. see [STR 86]) σ = ? q p = ? =⇒ min σ,q p :Φ p (σ,q p )=0 (-Ḋp (σ, q p )) min ∀σ,q p (L p (σ, q p , γp ), [11] where γp is the plastic multiplier and L p = -Ḋp (σ, q p ) + γp Φ p (σ, q p ).

[12]

The Kuhn-Tucker optimality conditions (e.g. see [STR 86]) for this problem can be readily obtained as:

0 = ∂L p (σ, q p , γp ) ∂σ = -ǫp + γp ∂Φ p ∂σ =⇒ ǫp = γp ∂Φ p ∂σ 0 = ∂L p (σ, q p , γp ) ∂q p = -ξp + γp ∂Φ p ∂q p =⇒ ξp = γp ∂Φ p ∂q p , [13]
which can be interpreted as the evolution equations of the internal variables of the plastic model with the Lagrange multiplier γp as the plastic multiplier. If we assume that the latter takes a zero value in an elastic process, the last Kuhn-Tucker equation associated with the Lagrangian in (11) leads to the standard form of the loading /unloading conditions which can be written as γp ≥ 0 ; Φ p ≤ 0 ; γp Φ p = 0.

[14]

The positive value of the plastic multiplier can be obtained from the plastic consistency condition, imposing for a state of plastic loading the plastic admissibility on the subsequent state, resulting with 0 = Φp = ∂Φ p ∂σ σ + ∂Φ p ∂q p qp .

[15] >From (10) and (13) we get:

qp = ∂q p ∂ξ p ξp = -γp ∂ 2 Ξ p ∂ξ p 2 ∂Φ p ∂q p , [16]
and from (7) and (13):

σ = C e ( ǫ -ǫd -ǫp ) = C e ( ǫ -ǫd ) -γp C e ∂Φ p ∂σ . [17]
Combining the equations (15), ( 16) and (17) we obtain:

γp = ∂Φ p
∂σ C e ( ǫǫd )

∂Φ p ∂σ C e ∂Φ p ∂σ + ∂Φ p ∂q p ∂ 2 Ξ p ∂ξ p 2 ∂Φ p ∂q p [18]
Finally, using (17) and (18) we determine the elasto-plastic consistent tangent modulus (with damage kept fixed) that

σ = [C e - C e ∂Φ p ∂σ ∂Φ p ∂σ C e ∂Φ p ∂σ C e ∂Φ p ∂σ + ∂Φ p ∂q p ∂ 2 Ξ p ∂ξ p 2 ∂Φ p ∂q p
]( ǫǫd ).

[19]

Damage model

The damage model can also be cast in an equivalent format as the one given for the plasticity, although that is not the usual manner of presenting it (e.g. see [LEM 88] or [KRA 96]). Namely, in the case where the damage model is activated while the plasticity remains inactive, we can appeal to the principle of maximum damage dissipation to select among all admissible value of stress and hardening damage variables those which maximize the damage dissipation. The latter can be formulated as an constrained minimization problem and further transformed into an unconstrained one by using the Lagrange multiplier method σ = ? q d = ?

=⇒ min σ,q d :Φ d (σ,q d )=0 (-Ḋd (σ, q d ))

min ∀σ,q d (L d (σ, q d ), γd ), [ 20
]
where γd is the Lagrange multiplier for damage and

L d = -Ḋd (σ, q d ) + γd Φ d (σ, q d ).
[21]

The Kuhn-Tucker optimality conditions corresponding to the damage Lagrangian can readily be obtained as

0 = ∂L d (σ, q d , γd ) ∂σ = -Ḋσ + γd ∂Φ d ∂σ =⇒ Ḋ = γp 1 σ ∂Φ d ∂σ 0 = ∂L d (σ, q d , γd ) ∂q d = -ξd + γd ∂Φ d ∂q d =⇒ ξd = γd ∂Φ d ∂q d [22]
which specify the evolution equations for the internal variables of the damage model.

To complete this description of the Kuhn-Tucker optimality condition we also provide the loading/unloading conditions with

γd ≥ 0 ; Φ d ≤ 0 ; γd Φ d = 0, [23]
where the zero value of the damage multiplier is introduced for an elastic process, where no change of damage internal variables would take place. The positive value of the damage multiplier, which occurs for the case of damage loading can be computed from the damage consistency condition, by imposing the admissibility on the subsequent state to get

0 = Φd = ∂Φ d ∂σ σ + ∂Φ d ∂q d
qd .

[24] >From (10) and ( 22) we have

qd = ∂q d ∂ξ d ξd = -γd ∂ 2 Ξ d ∂ξ d 2 ∂Φ d ∂q d , [25]
and from (8) and ( 22) we get

σ = (D -1 ǫ d ) = D -1 ǫd -D -1 Ḋ σ D -1 ǫ d = D -1 ǫd -γd D -1 ∂Φ d ∂σ , [ 26
]
where we used the relation 24),( 25) and ( 26) we obtain:

(D -1 ) = -D -1 ḊD -1 . Combining equations (
γd = ∂Φ d ∂σ D -1 ǫd ∂Φ d ∂σ D -1 ∂Φ d ∂σ + ∂Φ d ∂q d ∂ 2 Ξ d ∂ξ d 2 ∂Φ d ∂q d
.

[27]

Hence using ( 26) and ( 27) we get the damage consistent tangent modulus as

σ = [D -1 - D -1 ∂Φ d ∂σ ∂Φ d ∂σ D -1 ∂Φ d ∂σ D -1 ∂Φ d ∂σ + ∂Φ d ∂q d ∂ 2 Ξ d ∂ξ d 2 ∂Φ d ∂q d
] ǫd .

[28]

Elasto-plastic-damage coupling

Finally, for the case where both plasticity and damage models are active one can apply simultaneously the principle of maximum plastic and maximum damage dissipation to recover both sets of results as presented in (7), ( 13) and (15) for plasticity and (8), ( 22) and (24) for damage. This would imply in particular that we can write two different forms of the stress rate equation, one for plasticity and another for damage, according to

σ = C ep ( ǫ -ǫd ), [29]
where where

C ep =    C e ; γp = 0 C e - C e
C ed =    D -1 ; γd = 0 D -1 - D -1 ∂Φ d ∂σ ∂Φ d ∂σ D -1 ∂Φ d ∂σ D -1 ∂Φ d ∂σ + ∂Φ d ∂q d ∂ 2 Ξ d ∂ξ d 2 ∂Φ d ∂q d ; γd > 0 [32]
Therefore, from equality of stress rates in ( 29) and (31) we can obtain the damage strain rate according to

C ep ( ǫ -ǫd ) = C ed ǫd =⇒ ǫd = [C ep + C ed ] -1 C ep ǫ.
[33]

Replacing the last result in (31) we can rewrite the stress rate equation as

σ = C epd ǫ, [34]
where C epd is the elasto-plastic-damage consistent tangent modulus is given as

C epd = C ed [C ep + C ed ] -1 C ep . [35] ǫ σ C e C ep C epd σ y σ f C e ǫ p ǫ d ǫ e
Figure 1. The σǫ diagram for the one-dimensional elasto-plastic-damage constitutive model. We first reach the yield surface and then the fracture threshold surface. By unloading we can observe both phenomena, the appearance of plastic strains and decreasing of elastic properties. We also give the graphical representation of the strain decomposition, ǫ = ǫ e + ǫ p + ǫ d .

Numerical Implementation

In this section we address the numerical implementation issues of the proposed damage-plasticity model, within the framework of one-step, implicit backward Euler time-integration scheme. The central problem of computational plasticity can thus be presented as follows:

Central problem of computational damage-plasticity

Given: ǫ n = ǫ(•, t n ), ǫ p n = ǫ p n (•, t n ), ǫ d n = ǫ d n (•, t n ), ξ p n = ξ p n (•, t n ), D n = D(•, t n ), ξ d n = ξ d (•, t n ), ∆t = t n+1 -t n Find: ǫ n+1 , ǫ p n+1 , ǫ d n+1 , ξ p n+1 , D n+1 , ξ d n+1
Such that:

A n el el=1 [f el,int n+1 -f el,ext n+1 ] = 0, f el,int n+1 = le B el σ(ǫ n+1 , ǫ p n+1 , ξ p n+1 , D n+1 , ξ d n+1 , ǫ d n+1 )dx, f el,ext n+1 = le N el bdx + [N el t] ∂le , γp n+1 ≥ 0, Φ p n+1 (σ n+1 , q p n+1 ) ≤ 0, γp n+1 Φ p n+1 = 0 γd n+1 ≥ 0, Φ d n+1 (σ n+1 , q d n+1 ) ≤ 0, γd n+1 Φ d n+1 = 0. [36]
In (36) above the first equation expresses equilibrium of the complete system, whereas the other two are the equations imposing the plastic and damage admissibility of the computed stress state at time t n+1 . The symbol A n el el=1 denotes the standard finite element assembly procedure summing up all element contributions (e.g. see [BAT 90]), f el,int and f el,ext represent the internal and external forces of the element ′ el ′ , N el is element displacement variation and B el the corresponding strain variation. The essential role of the finite element method and the commonly used numerical integration scheme is to drastically simplify the central problem reducing it to corresponding numerical integration points which are chosen according to the given finite element mesh. Furthermore, in accordance to the usual operator split procedure (e.g. see [SIM 98]), we can separate the solution of the equilibrium equations from the one carried out to find the admissible stress state. The first group of equations, namely the equilibrium equations are solved iteratively in order to provide the current guess of the total strain field. If the Newton method is employed, one would need the tangent stiffness matrix, or more precisely the elasto-plastic-damage tangent modulus (e.g. see [BEN 88]). However, this is not essential and many other iterative schemes can be used for such a purpose. For any global solution procedure we are given the best iterative value for total strain field at each integration point, ǫ

σ n+1 = C e (ǫ (i) n+1 -ǫ p(k) n+1 -ǫ d(k) n+1
).

[37]

As already indicated, the computation of this kind is carried out independently on plasticity and damage part of the model, producing the stress values, σ p n+1 and σ d n+1 , respectively. At the final stage, we then impose that these two stresses coincide, which provides the converged values of the plastic and damage deformation.

Plasticity computation

The computation for the plasticity model is carried out by the procedure described in this section. By keeping the damage variable fixed and, in particular, the current value of the damage strain, ǫ d(k) n+1 , frozen, we first assume that the step remains elastic (with respect to the potential change of plasticity model internal variables), which allows us to compute the elastic trial state:

σ p,trial n+1 = C e (ǫ (i) n+1 -ǫ p n -ǫ d(k) n+1 ) [38]
q p,trial n+1 = q(ξ p n ).

[39]

If the trial elastic step is indeed plastically admissible in the sense that Φ p,trial n+1 = Φ p (σ p,trial n+1 , q p,trial n+1 ) ≤ 0, we conclude that the total step is the exact solution with γp n+1 = 0. In the opposite, we need to correct this result for the stress by finding the positive value of the plastic multiplier which will re-establish the plastic admissibility. To that end we first integrate by backward Euler scheme the evolution equation in (13) to obtain the correction of internal variables:

ǫ p n+1 = ǫ p n + γ p n+1 ∂Φ p n+1 ∂σ p n+1 ξ p n+1 = ξ p n + γ p n+1 ∂Φ p n+1 ∂q p n+1 q p n+1 = q p (ξ p n+1 ), [40]
where

γ p n+1 = γp n+1 ∆t [41]
and the corresponding value of the stress,

σ p n+1 = C e (ǫ (i) n+1 -ǫ p n+1 -ǫ d(k) n+1 ) = C e (ǫ (i) n+1 -ǫ p n -ǫ d(k) n+1 ) -C e (ǫ p n+1 -ǫ p n ) = σ p,trial n+1 -C e γ p n+1 ∂Φ p n+1 ∂σ n+1 , [ 42]
which should be tested for the plastic admissibility,

Φ p n+1 = Φ p (σ p n+1 , q p n+1 ) = 0 . [43]
Using the yield surface function defined in (5) and the equations ( 40) and ( 42), (43) can be rewritten as

Φ p n+1 = |σ p,trial n+1 | -C e γ p n+1 -(σ p y -q p (ξ p n+1 )) = 0 . [44]
For a convex hardening function (e.g. saturation hardening ) we can always solve the last equation using the Newton procedure to get converged value γ p n+1 from which we can calculate

ǫ p n+1 = ǫ p n + γ p n+1 sign(σ p,trial n+1 ) ξ p n+1 = ξ p n + γ p n+1 σ p n+1 = σ p,trial n+1
-C e γ p n+1 sign(σ p,trial n+1 ) .

[45]

Damage computation

We then turn towards the damage part of the model to compute the potential evolution of the damage internal variables. We start again this computation by assuming the elastic trial step and keeping the damage variables the same as in the previous step, which leads to

σ d,trial n+1 = D -1 n ǫ d(k) n+1 [46] q d,trial n+1 = q(ξ d n ).
[47]

We recall that the damage strain value at t n+1 , ǫ d n+1 , is determined by the iteration process, where for either plastic and damage calculation at kth iteration we keep the damage strain value frozen at ǫ

d(k) n+1 .
If the damage admissibility of the trial step is confirmed with

Φ d,trial n+1 = Φ d (σ d,trial n+1 , q d,trial n+1 ) ≤ 0 ,
the trial step is indeed the exact solution with γ d n+1 = 0. In the opposite, one needs to compute the corresponding positive value of the damage multiplier that will reestablish the damage admissibility. More precisely, integrating by backward Euler scheme the damage variables evolution equations in (20) we can write

σ d n+1 = D -1 n+1 ǫ d(k) n+1 = σ d,trial n+1 -D -1 n γ d n+1 ∂Φ d n+1 ∂σ d n+1 D n+1 = D n + γ d n+1 1 σ d n+1 ∂Φ d n+1 ∂σ d n+1 ξ d n+1 = ξ d n + γ d n+1 ∂Φ d n+1 ∂q d n+1 q d n+1 = q d (ξ d n+1
).

[48]

By using the damage function form introduced in (5) and the updates in (48) above, the damage admissibility condition can be written as

Φ d n+1 = |σ d,trial n+1 | -D -1 n γ d n+1 -(σ d f -q d (ξ d n+1 )) = 0. [49]
Solving this non-linear equation by the Newton method gives us the value of the damage multiplier γ d n+1 and hence the final values of stresses and the damage internal variables:

σ d n+1 = σ d,trial n+1 -D -1 n γ d n+1 sign(σ d,trial n+1 ) D n+1 = D n + γ d n+1 1 |σ d n+1 | ξ d n+1 = ξ d n + γ d n+1 q d n+1 = q d (ξ d n+1
).

[50]

In equation (46) the trial stress is calculated assuming that the initial value of the damage variable is non-zero, D 0 > 0, to avoid a potential problem of dividing by zero. By this assumption we have to change the elastic property values in the program to Ce , so that the real elastic properties of the model do not change:

C e = (D 0 + Ce -1 ) -1 . [51]
However, by taking the value of D non-zero, we slightly change the damage hardening variable evolution since D(ξ d = 0) > 0. Thus, D 0 becomes another material parameter for the hardening evolution. In our calculations we took D 0 = C e-1 , which means that within the present model calculations the elastic properties have to be Ce = 2C e .

Coupling

We note that the computation on the plastic model described in the previous section is carried out in a completely independent manner from the computation on the damage part of the model summarized in (46) to (50); In particular these two computational procedures can be advanced in parallel with no exchange of results, dedicating to each its particular processor on a parallel computer, to carry out the corresponding part of the work.

The final results produced by the two parts of the model should finally be compared against one another and any discrepancy, or residual, should be eliminated with r(ǫ

d(k) n+1 ) := σ p n+1 (ǫ d(k) n+1 ) -σ d n+1 (ǫ d(k) n+1 ) = 0. [52]
The last equation is then solved by an iterative procedure to which we supply the results available from the two parts to perform the next iterative step with

0 = r (k) n+1 -C ep n+1 ∆ǫ d(k) n+1 -C ed n+1 ∆ǫ d(k) n+1 [53] leading to ∆ǫ d(k) n+1 = r (k) n+1 C ep n+1 + C ed n+1 ǫ d(k+1) n+1 = ǫ d(k) n+1 + ∆ǫ d(k)
n+1 .

[54]

The computation then continues with this improved value of the damage strain until the convergence in ( 52) is finally achieved. The number of iterations depends significantly on the choice of the starting value of the damage strain, ǫ d (k=0) n+1

. As the first guess we should not take the last converged value, ǫ d n , since the damage strain evolves also in an elastic process (ǫ d = Dσ = const. even if D = const). So, we first assume that process remains elastic, which means that the damage variable, D, does not change and the trial stress equals:

σ trial n+1 = σ n + (D n + C e -1 ) -1 (ǫ n+1 -ǫ n ). [55]
By definition the corresponding damage strain is equal to:

ǫ d (k=0) n+1 = D n σ trial n+1 . [56]
The value above is then taken as the starting value of the damage strain iteration process to solve the non-linear case (equations ( 52) to (54)).

In (54) above we need the corresponding values of the elasto-plastic as well as the elasto-damage moduli, which can be easily computed from (45) for plasticity and from (50) for damage to obtain:

∂σ p n+1 ∂(ǫ n+1 -ǫ d n+1 ) =: C ep n+1 =    C e ; γ p n+1 = 0 C e ( ∂ 2 Ξ p ∂ξ p 2 )n+1 C e +( ∂ 2 Ξ p ∂ξ p 2 )n+1 ; γ p n+1 ≥ 0 ∂σ d n+1 ∂ǫ d n+1 =: C ed n+1 =      D -1 n ; γ d n+1 = 0 D -1 n ( ∂ 2 Ξ d ∂ξ d 2 )n+1 D -1 n +( ∂ 2 Ξ d ∂ξ d 2 )n+1 ; γ d n+1 ≥ 0 . [57]
The given values of tangent moduli for both parts of the model are also of direct use for constructing the tangent modulus of the coupled model, according to

C ep n+1 (dǫ n+1 -dǫ d n+1 ) = C ed n+1 dǫ d n+1 [58] =⇒ dǫ d n+1 = C ep n+1 C ep n+1 + C ed n+1 dǫ n+1 =⇒ dσ n+1 = C ep n+1 C ed n+1 C ep n+1 + C ed n+1 C epd dǫ n+1 .
[59]

The tangent elasto-plastic-damage modulus given in ( 59) is used to compute the corresponding element contribution to the global stiffness matrix for the Newton iterative scheme applied in solving global equations with

A n el el=1 [K el,(i) n+1 (d el,(i+1) n+1 -d el,(i) n+1 ) = f el,ext n+1 -f el,int n+1 ] K el,(i) n+1 = le B el C epd n+1 B el,T dx. [60]
The computational procedure presented in the foregoing can be generalized to 2D or 3D problems with no modification other than dealing with the second order tensors used for stress and strain, and the fourth order tensor for the tangent modulus. The outline of such a computational procedure is given in Appendix A.

Numerical simulations

In this section we present the results of a couple of numerical simulations. In order to illustrate the versatility of the model to represent large spectrum of different materials, one of the example draws from porous metal failure and another from compacting concrete .

Criteria for porous metals in tension

The porous metal coupled model was built along the lines of the pioneering work of Gurson ([GUR 77]), however with important difference regarding the present model, which has the ability to describe the closing of pores at unloading. Postulating that it is only spherical part of stress which determines the porosity, the damage criterion is given as

Φ d (σ, q d ) = tr(σ) -(σ d f -q d ), [61]
where tr(σ) denotes the trace of the tensor σ and < • > the Macauley brackets:

< x >= x ; x ≥ 0 0 ; x < 0 . [62]
Here we neglect the possibility that the material can be damaged in compression. To model the plasticity of metal matrix, we used the von Mises criterion:

Φ p (σ, q p ) = dev(σ) : dev(σ) -(σ p y -q p ), [63]
where dev(σ) denotes the deviatoric part of the tensor σ, dev(σ) ≡ σ -1 3 tr(σ). >From the choice of the criteria it follows that the evolution of damage variables depends only upon the spherical part of the stress tensor and the evolution of plastic variables upon its deviatoric part. Hence, the two nonlinear phenomena appear uncoupled in strain space. This is the direct consequence of the initial physical presumption that the opening of the micro-cracks is due to positive spherical part of the stress and sliding of crystal planes due to the stress deviator. The former corresponding to damage and the latter to plasticity.

Finally, we use an exponential hardening law for either phenomenon, plasticity and damage,

q p (ξ p ) = (σ p y -σ p ∞ )(1 -e -b p ξ p ) q d (ξ d ) = (σ d f -σ d ∞ )(1 -e -b d ξ d ), [64]
where σ p ∞ and σ d ∞ are saturation values of stress, whereas b p and b d are the material parameters governing the rate of saturation.

Results

The model is illustrated on an example of a rectangular plate with a circular hole in the middle, submitted to a simple tension test. By exploiting symmetry conditions, only one quarter of the model is used in the analysis; See Figure 2. The material properties taken in the calculation were the following; (i) for elasticity: Young's modulus, E = 240GP a and the shear modulus, µ = 92GP a; (ii) for plasticity: yield stress, σ y = 170M P a, hardening limit stress, σ p ∞ = 210M P a and saturation parameter, b p = 50; (iii) for damage: fracture stress, σ f = 170M P a, hardening limit stress, σ d ∞ = 210M P a and saturation parameter, b d = 50.

In Figures 3 and4 we show how the spreading of plastified and damaged regions will change with the other phenomenon being activated. Different stages of activation of either plasticity or damage models are illustrated by contours of hardening variables ξ p and ξ d , respectively. We observe the complete disappearance of shear band (Figure 3), a typical response of metals or alloys with von Mises criterion, when damage is also taken into account. Besides, we notice that in the case where both phenomena are activated either region is reduced to a smaller volume, but the differences between the maximum and minimum value of ξ p and ξ d is larger. With other words, the phenomena are, when activated simultaneously, more localized.

Convergence characteristics

Since the plastic and damage model state variables computations are carried out independently, the convergence rate of each one is not affected by another. In particular, the implementation of the von Mises model is completely the same as for the plasticity model and thus keeps its robustness in the coupled model as well. Moreover, since the damage model is formulated in a completely analogous way to the plasticity, that is using the fracture criterion and the principle of the maximum dissipation, we obtain comparable computational efficiency and similar robustness. In both cases, for saturation hardening models defined in (64), we need 2 -8 iterations to converge with the Newton procedure applied to each model.

The iterative procedure for the coupled plasticity-damage model, assuring the final uniqueness of stresses, requires in general only a few iterations to converge. The convergence criteria at k-th iteration can be written as:

3 i,j=1 (σ p (ǫ d(k) ) -σ d (ǫ d(k) )) 2 ≤ η min(σ y , σ f ), [ 65
]
where η is a chosen tolerance The solution uniqueness for the model of this kind is guaranteed if the plasticity and damage criteria do not concern the same inelastic mechanisms. Nevertheless, uniqueness of the solution does not ensure the easy convergence. We have remarked that for the plastic and damage criteria too similar, we need smaller time steps to converge, while the convergence is remaining quadratic.

Criteria for compacting concrete

It is quite interesting that the same kind of model can be used to represent compacting concrete behavior, although the latter can be seen as placed at the opposite end of the spectrum with respect to porous metal behavior described by Gurson's criterion. Namely, while the previous study deals with the increase of metal porosity and softening response under tension loading, in this example we model the decrease of concrete porosity under hydrostatic pressure loading and increasingly hardening response of compacting concrete [ GAT 99].

Since we do not necessarily seek to use the most elaborate model for concrete capable of providing a very good agreement between the numerical simulations and experimental results, we use rather a simple model for plasticity component in terms of the Drucker-Prager criterion with no hardening:

Φ p (σ) = dev(σ) : dev(σ) -tan(α) 1 3 tr(σ) -2/3 σ p f , [ 66]
where tan(α) is a (positive) material constant, roughly representing internal friction.

The damage model component is chosen to describe an essential mechanism of compacting concrete "damage" leading to an increase in resistance. Therefore, the damage criterion for concrete compaction concerns only the spherical part of the stress:

Φ d (σ, q d ) = tr(σ) -(σ d f -q d ), [67]
where σ d f is the damage elastic limit and q d is the damage hardening variable. By introducing a linear hardening law for the damage:

q d (ξ d ) = -K d ξ d [68]
with fairly large value of hardening modulus K d we hope to provide a reliable representation of concrete behavior in compaction.

Results

The computations are performed for the same model of a perforated plate (see Figure 2), but this time submitted to compression loading. The numerical values for the material parameters taken in the calculation are the following: (i) elasticity: Young's modulus E = 240GP a and the shear modulus µ = 92GP a; (ii) plasticity: yield stress, σ y = 170M P a, material parameter, tanα = 0.6; (iii) damage: fracture stress, σ f = 210M P a, K d = 200.

In Figure 5 we show how the damaged regions change with or without the plasticity component activated, by tracing the contours of the damage hardening variable ξ d . We observe that the damage hardening is much less present when both model components are activated.

Conclusions

The general theoretical framework for developing a coupled plasticity-damage constitutive model presented herein allows to accommodate a wide variety of materials, as illustrated on porous metal and compacting concrete. The main novelty of the model with respect to the standard coupled models of this kind pertains to treating each model component, plasticity and damage, in an independent manner attributing to each its own yield or damage criteria. This kind of feature may especially be advantageous when the behavior of one or both components of the coupled model is well under control. In the opposite case, in order to simulate a real material behaviour, we need to identify the material parameters, either in the phenomenological sense by comparing it with experimental results or by using a microscale-based representation to provide the material parameters for a particular component.

Another important novelty introduced herein concerns the numerical implementation, where the plastic and damage component state variables computations are carried out independently from each other. The latter allows that the classical return mapping algorithm providing the quadratic convergence rates is directly used for plasticity component computations. Moreover, since the damage model is formulated in a completely analogous way to the plasticity, that is using the damage criterion and the principle of the maximum dissipation, the same type of return mapping algorithm can be used to obtain computational efficiency and robustness comparable to those of plasticity component.

The chosen strategy of separating the computations for two components to reduce the model complexity is in general taxed mildly, since the additional computational cost concerns typically only a few iterations to obtain the correct value of damage strain and to converge the stress provided by two components to the same value.

  the superscript ′ (i) ′ denotes the iteration counter. The central problem thus reduces to computing the corresponding values of the internal variables which will provide an admissible stress field. This computation is accomplished by a local iterative scheme, with iteration counter denoted by ′ (k) ′ , where two ingredients of the model exploited independently produce the current values of plastic and damage strains, ǫ p(k) n+1 and ǫ d(k) n+1 , and the corresponding stress field according to the following expression:

Figure 2 .

 2 Figure 2. One quarter of the specimen, used in calculation and prescribed boundary conditions .

Figure 3 .Figure 4 .

 34 Figure 3. Contours of the plastic hardening variable ξ p , when a) only plasticity is activated and b) both plasticity and damage are activated.

Figure 5 .Figure 6 .

 56 Figure 5. The value of the damage hardening variable ξ d , when a) only damage is activated and b) both, plasticity and damage are activated.In Figure6we show contours of the largest stress component in terms of its absolute value, σ 11 , for three different cases where only the plasticity, only the damage or both model components are activated. It is for the latter case when the stress response appears least localized.

Table 1 .

 1 . A typical sequence of the intermediate results obtained in the iterative process is presented in Table 1. Iterative values for stress difference norm and damage strain computed by plasticity-damage model.

	iter. 1	σ p -σ d 2 2.400521 10 17 0.01681690 0.006922074 -0.002569930 ǫ d 11 ǫ d 22 ǫ d 12
	2	7.003845 10 12 0.01614148 0.007330000 -0.002210263
	3	3.758989 10 9	0.01615158 0.007325925 -0.002293375
	4	1.504535 10 2	0.01615126 0.007326051 -0.002292093
	5	9.815405 10 -5 0.01615126 0.007326049 -0.002292127

Acknowledgements:

This work was supported by the French Ministry of Research and by the Slovene Ministry of Education, Science and Sport. This support is gratefully acknowledged.

Appendix: Numerical algorithm of the elasto-plastic-damage model

given

≤ 0 :

)) = 0 where where


Marion Berton
email: marion.berton@insa-lyon.fr

Laurent Maxit

Daniel Juvé

C Audoly

PREDICTION OF FLOW-INDUCED SOUND AND VIBRATION : ON DIFFERENT METHODS FOR INTRODUCING THE TBL EXCITATION IN THE VIBROACOUSTIC MODEL

come     Prediction of flow-induced sound and vibration: On different methods for introducing the TBL excitation in the vibroacoustic model

The modeling of the vibroacoustic response of a heavy fluid loaded panel excited by a turbulent boundary layer is investigated. The pressure fluctuations due to the Turbulent Boundary Layer (TBL) at the surface of the structure are characterized by cross-spectra which can either be expressed in the physical space or in the wavenumber space. Various models have been proposed in the literature to obtain closed form expressions of these spectra starting from global TBL characteristics. A difficulty to estimate the vibroacoustic response of the panel excited by TBL relies on the coupling of the deterministic vibroacoustic model with the statistical wall pressure fluctuations model. In the present paper, we study different techniques to achieve this coupling: (a), in the physical space; (b), considering an uncorrelated wall plane waves field; (c), using a reciprocity technique; (d), using realizations of the stochastic uncorrelated plane waves field. These techniques require calculating different transfer functions with a vibroacoustic model. For that, we use the Patch Transfer Function formalism (PTF). This approach allows us to couple the structure problem with the heavy fluid problem. The added mass effect of the fluid on the panel (i.e. fluid reactance) is taken into account in the modal frequencies of the panel whereas the acoustic radiation from the panel into the fluid (i.e. fluid resistance) is modeled using the real part of the patch acoustic impedances of the fluid domain. For validation purpose, the results of these vibroacoustic calculations are compared with a result from the literature for a test case composed of a simply supported rectangular thin plate immersed in water. Afterwards, the response of the panel excited by the turbulent boundary layer is predicted in the low and medium frequency domains using the four techniques described above. The response is characterized here by the velocity spectrum at one point on the plate. Convergence of the four techniques is studied. The results with different values of calculation parameters are presented and analyzed in order to find the optimal ones. We observe that the four techniques give very similar results for these optimal parameters. Finally, we discuss the interest of each technique.

Introduction

In the context of the underwater missions, submarine vehicles should be designed to minimize noise sources and acoustic radiation. One of these noise sources is the turbulent flow induced by the vehicle movement. Numerical modeling can help to improve knowledge of interaction between Turbulent Boundary Layer (TBL) and heavy fluid loaded structures. The classical method consists in characterizing the wall pressure fluctuations by a spectrum depending on the TBL parameters, and then to inject it in a vibroacoustic model. The coupling between the statistical model used to describe the wall pressure fluctuations and the deterministic vibroacoustic model constitute a difficult point in the methodology. Several approaches to achieve this coupling are investigated in this paper.

Different methods for introducing TBL excitation in the vibroacoustic model

Four methods to evaluate the vibro-acoustic response to TBL excitation are briefly described in this section. Let us consider a fluid loaded baffled panel of surface Σ excited by a fully developed, stationary and homogeneous TBL of wall pressure field , at point as a function of time. The assumptions are made that the wall pressure fluctuations are not affected by the vibrations of the plate, and that the acoustic waves do not interfere with the flow. The spacefrequency spectrum of the plate velocity at point in response to a discrete spatial distribution of point pressure forces is defined as [START_REF] Strawderman | Wavevector-Frequency Analysis with Applications to Acoustics[END_REF][START_REF] Maxit | Prediction of flow induced sound and vibration of periodically stiffened plates[END_REF]:

, = ∑ ∑ * , , , , , , ∈ ∈ , (1)
where: -is the discretization of Σ defined with a criteria which will be discussed later; -, , is the space-frequency cross-spectrum of the wall pressure between points and ; -, , is the transfer function in velocity at point driven by a point force at point . It can be obtained using a vibro-acoustic model considering a driving point force at . The method to estimate the panel response with (1) and with , , evaluated from a vibro-acoustic model with a driving point force at is called the "spatial method".

The TBL excitation can also be described by the space Fourier transform of , , which gives the wavenumber-frequency spectrum Φ , . The discrete space Fourier transform of , , is defined by:

, , = ∑ , , ∈ . (2)
The combination of equations ( 1) and ( 2) and the use of the wavenumber-frequency spectrum Φ , lead to:

, = ∑ , , ! Φ , ∈" . (3)
# is a truncation and a sampling of the wavenumber domain. Criterion defining # will be discussed later.

, ,

can be interpreted as the velocity response at point of the plate, due to a generalized plane wave of wavevector = $ % , $ & , exciting the plate at points (see the definition of the generalized plane waves in [START_REF] Aucejo | Experimental simulation of turbulent boundary layer induced vibrations by using a synthetic array[END_REF]). The method to estimate the panel response with (3) and with , , evaluated from a vibro-acoustic model considering a generalized wall plane wave excitation is called the "uncorrelated plane waves method".

Another method [START_REF] Maxit | Prediction of flow induced sound and vibration of periodically stiffened plates[END_REF] is based on the reciprocity principle, which indicates that the ratio of the velocity at point over the normal force at point is equivalent to the ratio of the velocity at point over the normal force at point :

, , = , , . (4)
This expression can be introduced in (2), and , , can then be interpreted as the space Fourier transform of the response for a point force at . It results that only one single force should be considered in the vibro-acoustic model to deduce , , . The method to estimate the panel response with (3) and using the reciprocity principle to evaluate , , is called the "wavepoint reciprocity method".

The last method consists in considering several realizations of the stochastic uncorrelated plane waves. For each sample (i.e. realization), a random phase is assigned to each wall plane wave of amplitude 'Φ , /2*. The deterministic wall pressure of the sample is then:

= ∑ + , 'Φ , /2* ∈" , (5)
whereis the random phase associated to the wall plane wave of wavevector . The panel response to the TBL excitation is estimated from the average of the panel response over the samples of plane waves. This method is called the "sampling of uncorrelated plane waves method".

The four presented approaches are based on equations (1) or [START_REF] Aucejo | Experimental simulation of turbulent boundary layer induced vibrations by using a synthetic array[END_REF]. Criteria should be considered to define the discretization spaces and #. For (related to (1)), spatial discretization has to be fine enough to well describe the correlation decrease between wall pressure signals at points and . For # (related to (3)), the maximal wavenumber $ ./% has to be chosen to include the parts of the wall pressure spectrum that impact the response of the panel. Moreover the wavenumber resolution $ should be small enough to take properly into account the variations of , , and Φ , over . These criteria will be studied and discussed on numerical examples in Sec. 4.

Patch Transfer Functions (PTF) formalism validation

In this paper, the different transfer functions of the fluid loaded panel will be evaluated with the PTF (Patch Transfer Functions) approach ( [START_REF] Maxit | Improving the Patch Transfer Function approach for fluid-structure modeling in heavy fluid[END_REF]). It is a substructuring approach which allows us to decompose our problem into two parts (i.e. the panel and the semi-infinite fluid). The coupling interface is divided in patches which size depends on the considered wavelengths. Each part is characterized separately by Patch Transfer Functions (i.e. patch mobilities for the panel, patch impedances for the fluid). Writing the continuity conditions at the coupling interface enables to assemble the two parts. The particularity of the present model compared to [START_REF] Maxit | Improving the Patch Transfer Function approach for fluid-structure modeling in heavy fluid[END_REF][START_REF] Aucejo | Convergence acceleration using the residual shape technique when solving structure-acoustic coupling with the Patch Transfer Functions method[END_REF] is that the fluid added mass effect is taken into account through the "wet modal frequencies", (which are estimated by assuming the fluid incompressible) instead of using the imaginary part of the acoustic impedance of the fluid domain. This allows to overcome the convergence issue evoked in [START_REF] Aucejo | Convergence acceleration using the residual shape technique when solving structure-acoustic coupling with the Patch Transfer Functions method[END_REF] concerning the patch size criterion. Here, a patch size lower than half the wavelength gives converged results.

The test case used here to validate the PTF code corresponds to the case studied by Berry [START_REF] Berry | A new formulation for the vibrations and sound radiation of fluid-loaded plates with elastic boundary conditions[END_REF]. A simply supported plate immersed in water is excited by one point force at its center. Its dimensions are 455 mm x 375 mm x 1 mm and its structural damping is taken to be 0 = 0.01. The critical frequency is around 220 kHz. The obtained mean quadratic velocities over the plate surface are compared in Figure 1 between 10 Hz and 500Hz, with a comparison to the PTF calculation without taking the added mass effect into account. The results show that the modal frequencies calculated with the PTF method with added mass effect are satisfactory in comparison to the case without added mass effect. The frequency shift between the corresponding curves is well visible: for example the resonant peak at 121 Hz in the case without added mass effect is shifted to 43 Hz for both Berry's results and the PTF results.

Panel response calculations to TBL excitation

This section investigates the structural response of the plate described in Sec. 3 under a water flow with a convective velocity 4 5 = 57/8. The TBL wall pressure fluctuations which excite the panel are described here by the Corcos cross-spectrum model [START_REF] Corcos | Resolution of pressure in Turbulence[END_REF], either in the physical space or the wavenumber space. The observation point M of the plate velocity response is arbitrarily fixed at 9 = 5.4 <7 and = = 18.6 <7. The transfer functions , , and , , are calculated with PTF. We propose here to study the results of the different methods described in Sec. 2 in function of different calculation parameters. We indicate that the considered frequency band [10 Hz-500 Hz] is above the hydrodynamic coincidence frequency (which is around 2 Hz). For the "spatial method", the spatial resolution is set dependant to a characteristic wavelength A 5 = 4 5 /B in order to adapt the spatial discretization to the decrease of correlation of the wall pressure fluctuations. Results with different spatial resolutions are shown on Figure 2. A spatial resolution 9 = A 5 /3 is a just balance between time computation cost and satisfactory precision (computing time of 660 seconds when 9 = A 5 /3, against 1440 seconds when 9 = A 5 /4). However the computation time becomes prohibitive when the frequency increases because A 5 decreases, and the number of excitation points which should be considered becomes large.

The criterion related to the "uncorrelated plane waves method" concerns the wavenumber range defined with $ ./% and $: at the first sight, the analysis of the wavenumber-frequency spectrum Φ , leads us to consider a maximum wavenumber $ ./% greater than the convective wavenumber $ 5 = /4 5 . However, the dynamic behaviour of the plate can filter the excitation and decrease the impact of the convective region of the wall pressure spectrum on the plate response. If it is the case, $ ./% could be defined from the natural bending wavenumber of the structure $ D (which is well lower than $ 5 for the present case: at 500 Hz $ 5 ~630 7 + and $ D ~637 + ). Figure 3 shows the results for the two cases: $ ./% = 2$ D and $ ./% = $ 5 F G$ with G$ = 1007 + . The results show that the convective region of the TBL wall pressure spectrum does not influence the plate response in the considered frequency range. It is thus unnecessary to take into account the convective part of Φ , . The criterion $ ./% = 2$ D can be applied for this case, which saves computing time. We underline that this criterion can be questioned in regards to the TBL model. For example, as the Chase model presents a lower pressure level than the Corcos model in the low wavenumber region, the convective part could have a higher contribution on the plate response. The converged results for both "spatial" and "uncorrelated plane waves" methods are also compared in Figure 3, and show a good agreement.

The result of the "wave-point reciprocity method" is proposed in Figure 4 a. A good agreement with the result of the "uncorrelated plane waves method" is observed. Although one single point force excitation should be considered in the vibro-acoustic model, the time cost (due to the space Fourier transform) appears to be higher than for the "uncorrelated plane waves method" for this case (650 seconds for the reciprocity method against 460 seconds for the plane waves method).

Finally, the results obtained from the last method, the "sampling of uncorrelated plane waves", are presented in Figure 4 b. The results for each realization are plotted in grey lines whereas the mean values of the 10 realizations are plotted in black line. The figure shows a good agreement with the "uncorrelated plane waves method" for a nearly equivalent computation time. One advantage of this approach is that few load case (here 10) should be considered in the vibroacoustic model. If a commercial vibro-acoustic code is used to evaluate the transfer functions, the small number of load case allows reducing the amount of input/output data compared to the "uncorrelated plane waves method".

Conclusion

Four methods have been tested and compared with good agreement in order to introduce the TBL excitation in the calculation of the response of an elastic plate submitted to a water flow. The "spatial method" is the most time consuming. In order to increase its efficiency, investigations could be carried out to limit the excitation to an area around a reference point. The process based on the Cholesky decomposition of the space-frequency cross-spectrum of the wall pressure described in [START_REF]Installation, operations, theory and utilities. Chap. 20: Random excitation[END_REF] could be an alternative to reduce the computing time. The "uncorrelated plane waves method" has shown its efficiency when the cut-off wavenumber $ ./% can be fixed in regards to the panel characteristics (i.e. natural bending wavenumber) and not of the excitation characteristics (i.e. convective wavenumber). Additional investigations should be done to study the influence of the hydrodynamic part of other TBL models. The "wave-point reciprocity method" has also been validated. For the studied case, it has not allowed to reduce the computing time. This approach is more appropriate for stiffened structure ( [START_REF] Maxit | Prediction of flow induced sound and vibration of periodically stiffened plates[END_REF]) for which the point response in the wavenumber domain (i.e.

, , ) can easily be estimated. Finally the technique using samples of plane waves with random phase gives also a good accuracy with few realizations. It constitutes another way to introduce the TBL excitation with few load case in the vibro-acoustic model.

Figure 1 :

 1 Figure 1 : Mean quadratic velocity of a thin plate driven by a point force at its centre.

Figure 2 :

 2 Figure 2 : Velocity spectrum at point M for different spatial resolutions @ .

Figure 3 :

 3 Figure 3 : Velocity spectrum at point M for different cut-off wavenumbers, and comparison with the spatial method.

Figure 4 :

 4 Figure 4 : Velocity spectrum at point M. a.: Reciprocity method. b.: Sampling method.

Acknowledgements: This work was performed within the framework of the LabEx CeLyA of Université de Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ ANR-11-IDEX-0007). This work was performed with the financial support of DCNS.