3D MODELLING OF CONCRETE FOR EARTHQUAKE ANALYSIS: DAMAGE MECHANICS AND PLASTICITY COUPLING

This contribution aims at presenting constitutive equations for concrete-like materials modeling in the framework of cyclic or seismic 3D analysis. A continuum damage based approach is considered coupled to plasticity theory allows accounting for localized damage in tension, smeared one in compression, permanent strains, hysteretic effects and unilateral features for load reversal conditions.


INTRODUCTION
Within the framework of earthquake analysis of critical civil engineering structures such as industrial sites or nuclear power-plants, one has to consider 3D finite elements to ensure the predictive computations of main local mechanisms such has crack opening and hysteretic energy dissipation.Simplified approaches, based on kinematics reduction, fail to represent these two levels of analysis.Quasi-brittle materials such as concrete deal with complex local mechanisms such as stiffness reduction due to damage, stiffness recovery for alternate loadings (case of earthquakes), damping increase due to frictional sliding and permanent displacement due to permanent strain.The expressions of material constitutive equations in a 3D framework still remain a challenge if numerical robustness for large scale computations is expected In this contribution, a new constitutive model for quasi-brittle materials subjected to earthquake loadings is presented.The main features are the asymmetry between tension and compression, permanent strains and unilateral effect, ensuring not only the relevancy but also the reliability of numerical predictions.The hysteretic effects are also included in the model since their importance for seismic applications has been emphasized in literature.Mainly driven by cracking, the tensile behaviour is governed by damage mechanics in a nonlocal framework avoiding localisation during strain softening.Considering smeared nonlinear behaviors, plasticity framework is adopted for compression state of stresses.To ensure numeri-cal robustness, hysteretic effects are only coupled to damage.The numerical implementation of the model is presented and local tests are shown for validation.
Recent experiments performed on CEA-Saclay shaking tables emphasized the 3D and torsional responses of a reinforced concrete mock up subject to strong earthquakes (SMART program).The comparisons between experiments and numerical computations help one to appreciate the relevancy of the approach regarding global responses (displacement, ductility) end local behaviors such as cracks pattern.

CONSTITUTIVE FORMULATION
This section aims at exposing a new constitutive model expressed within the framework of isotropic continuum damage mechanics and plasticity for seismic applications.It is based on the following observations: In tension, quasi-brittle materials exhibit localized cracking although in compression diffuse (or smeared) cracking appears.Therefore, isotropic damage is used in tension and plasticity is used in compression.The use of plasticity in compression is frequent in the field of soil mechanics when nonlinearities, dilatancy and permanent strains must be described [1].Tension and compression can be split into two distinct parts by considering the sign of the Cauchy stress tensor in 1D (or a related quantity in 3D).Since hysteretic phenomena are linked with frictional sliding, this mechanism is considered in tension, allowing a realistic description not only of hysteretic loops but also of permanent strains in tension.The unilateral effect [2] is also taken into account using a closure function which ensures the continuity stress/strain relation whatever the loading path.With these assumptions, one can state that the behavior of a representative elementary volume is accurately represented in tension (brittleness, hysteretic loops and permanent strains) and globally described in compression (nonlinearity and permanent strains).Indeed, close to the failure, the hysteretic effects are more due to steel plasticity than to cracking and frictional sliding.

State equations
The reversible part of the mechanical behavior can be defined thanks to the state equations.They are obtained by differentiating the state potential with respect to the state variables.The first state law allows to define the second-order Cauchy stress tensor: Where ρ is the material density, Ψ the Helmohltz free energy, ijkl C the fouth-order Hooke's tensor (computed from elastic material parameters) and ij ε is the second-order total strain ten- sor, d is the scalar damage variable, π ε kl is the second-order internal sliding tensor, η is the closure variable ranging from 1 (cracks opened) to 0 (cracks closed) that is used to cancel the permanent strains created in tension and p kl ε is the second-order plastic strain tensor.
The second state equation concerns the second order sliding tensor π σ ij allowing introducing hysteretic dissipation during crack opening and closing: At last, the energy released rate by damage mechanism may be expressed computing . To achieve the state equations definitions, the expression of consolidations or hardening functions should be performed.Those ones will be defined in the following sections.

Evolution laws and threshold functions
Three mains mechanisms are accounted for: damage, plasticity and hysteresis.These three irreversible processes require the expression of three threshold functions.
The damage mechanism and the isotropic hardenings are managed in a coupled way.Moreover, an associative low rule is assumed.The threshold surface d f is expressed as a function of the energy rate released by damage: Where Y is a part of the energy rate released by damage and is defined next, 0 Y is an initial threshold that must be overcome and Z is the thermodynamic force associated with isotropic hardening.Y is only a part of the damage energy release rate accounting only for the positive part of the elastic free energy : This definition allows keeping the explicit nature of the damage variable and therefore does not require local iterations, ensuring the numerical robustness of the integration of the flow rules associated with damage and isotropic hardening.Following the original works of Laborderie et al. [3][4], the damage evolution law is obtained as following, ( ) d A is a material parameter to be identified, drving the slope of the softening branch.
The internal sliding mechanism is assumed to be coupled with the kinematic hardening at the flow level.This approach aims at modelling the hysteretic effects for seismic engineering applications in a consistent way.The nonlinearities associated with the hysteretic effects are managed by a non-associative flow rule.The threshold surface is expressed in terms of sliding stress and back stress.One can note that there is no initial threshold in the expression of the threshold surface π f , allowing this mechanism to be activated only when damage is activated.

(
)( ) ij X is the back stress for kinematic hardening anf H(.) is the heaviside's function.The non- associative framework allows resolving the well-known problem of the kinematic hardening (if the rule had been associative) which is its linearity.Since the hysteretic effects are not liear, one needs to take their non-linearities into account.A pseudo potential of dissipation is introduced [5] ( )( ) a is a material parameter driving the hysteretic loop shape.The evolution of the internal sliding strain tensor is obtained thanks to consistency conditions introducing the sliding multiplier π λ & : The unilateral effect can be defined as an initial (undamaged) stiffness recovering when switching from tension to compression.This effect can be included in the model by various ways (changes of signs of strains, stresses or a coupling stress/strain expression).In this study, the choice to link the unilateral effect with the stresses has been made.As pointed out in literature, this choice allows describing the progressive closure of the cracks [3].A closure function is considered to gradually cancel the permanent strains created in tension when unloading.This function ranges from 1 (cracks opened) to 0 (cracks closed).In order to avoid local iterations to compute the stress from equation 5, a linear function is assumed.The following function is considered [3]: f σ is the closure stress, a material parameter to identify The plastic strain is managed thanks to a non-associative ow rule to ensure a satisfying description of the dilatancy.The threshold function has been choflsen as a Drucker-Prager criterion that is function of the Cauchy stress.

(
) where fp α , is a material parameter which manages the bi-compression response and c f an initial threshold that has to be overcome.One can notice that the stress state must lead to a global compression to activate plasticity.The consolidation function R has been chosen to describe the softening in compression and is expressed as: where Rp a and Rp b are material parameters that drive the shape of the softening.The pseudo potential of dissipation can be chosen as: p φ α is a material parameter that manages the dilatancy.The plastic strains flow rules are ob- tained by normality ryles :

Numerical implementations
The numerical implementation of a constitutive model lies in computing the stress increment knowing the thermodynamic state of the material at a given step.The proposed numerical scheme couples an explicit/implicit integration procedure.The flow rules associated with damage are explicitly integrated, allowing an exact computation of the damage variable.According to the kind of the stress state (tension or compression), the internal sliding and the plasticity mechanisms are treated thanks to a return mapping-based algorithm [6].

Local responses
A complex loading path is considered to emphasize numerical robustness of the proposed constitutive model as well as the relevancy of the local concrete behavior description.The loading is displacement controlled to capture the post-peak regime.In tension, the model exhibits a brittle behavior as well as hysteretic loops.Moreover, when switching from tension to compression, the cracks are progressively closed up to recover the initial stiffness.
In compression, permanent strains are taken into account as the only dissipative mechanism is plasticity.A linear behavior is recovered when unloading.The cracks are also re-opened progressively when a tension stress state appears after carrying out a compression (see figure 1).
The proposed constitutive model accounts for nonlinearities in compression thanks to a non associative plasticity criterion.One of the major drawbacks related to this approach is due to the plasticity theory.Indeed, the modulus of the unloading branch (in compression) is equal to the initial elastic stiffness.Nevertheless, in this study, a progressive crack closure/reopening in considered thanks to the closure function introduced previously.Therefore, the effective unloading modulus is not equal to the elastic one.To highlight this feature, the induced damage in compression is computed and plotted in figure 2. One can appreciate the ability of the model to account for stiffness decrease in compression, although using plasticity theory, depending on the level of pre-damage in tension.Another point of interest lies in the ability for the model to describe both the elasticity and the failure surfaces.Based on the results reported by Kupfer et al [7], the numerical elasticity and the failure surfaces are compared with the experimental ones.The numerical results are given in figure 3.One can observe that the effect of the hydrostatic pressure is well taken into account in the case of a biaxial compression loading path.On the contrary, the response in biaxial tension is not fully satisfying.Nevertheless, the kind of loading path is not common for application and gives a conservative property to the model.

Structural responses
It is now well-known that the use of continuum damage mechanics-based constitutive models leads to mesh-dependency issues.This numerical issue is mainly due to the loss of ellipticity of the local equilibrium equations.To avoid pathological mesh-dependencies issues, several approaches have been presented.The nonlocal approach proposed has been selected as it is available in the opened-source finite element software Cast3M [8].The local energy rate released by damage Y is averaged over given vicinity defined by a characteristic length that is usually chosen equal to 3 or 4 times the maximum aggregate size [9].
The aim of this section is to expose the ability of a complex damage constitutive model to be used in the framework of 3D seismic nonlinear analysis.The case study of the SMART mock-up is analyzed.Due to the design of SMART mock-up, it is possible to adopt different models for different parts of the structure.Indeed, the mock-up is designed to promote failure in the walls, and more accurately, design calculation predicts a rupture in the bottom of walls.Consequently, only the walls are modelled with the damage model using 8 Gauss point cubic elements.The other concrete parts of the structure are considered as elastic isotropic and linear.If the foundation and beams are discretized using 8 Gauss point cubic elements, the floors are made of DKT shell elements.This method aims at optimizing calculation cost.At last, steel reinforcement is modelled with elasto-plastic model and 1D bar elements.To check our model hypothesis, a modal analysis is performed in a first step.The complex shape of the mock-up makes occurring flexion-torsion coupled modes.Transient nonlinear analyses have been performed for the different sets of accelerograms.The 0.1g is below the design seismic level.It can be noted that the bottom of the wall is softly damaged, the mock-up do not remain elastic.But the maximum values of the spectrum correspond to initial modal eigen-values.At the 0.4g seismic level, the damage map shows that the mock up is really damaged.The bottom of the wall begins to be severely damaged and longer lintels present an important density of micro-cracks.The induced damage provokes a larger band of high intensity in the spectrum.Obviously the displacement is more important but the softening can be noted due to the increasing of the period.At last, the 0.8g seismic level is an over design intensity level.It is confirmed by the transient analysis: the spectrum becomes larger again, the wall is strongly damaged at the bottom and the lintels are crossed by damaged areas.

CONCLUSIONS
In this contribution, an original constitutive model devoted to seismic applications has been presented.Expressed within the theoretical framework of the 3D irreversible processes thermodynamics, an appropriate state coupling allows representing nonlinear effects such as cracking and hysteretic loops in a fairly satisfactory way.The compression is driven by nonassociative plasticity mechanism with isotropic hardening.Permanent strains are taken into account not only in tension but also in compression.The model accounts for a full and linear unilateral effect when switching from tension to compression, allowing a full recovery of the elastic (undamaged) stiffness.
The material parameters identification has been carried out and some local results (at the Gauss point level) have been exposed to show the main features of the constitutive law.The main problem related to the use of plasticity in compression concerning the stiffness decrease is solved fairly simply thanks to a progressive crack re-opening when switching from compression to tension.The ability of the model to represent the structural cyclic response of RC structure has been investigated thanks to the 3D SMART case-study.

Figure 2 .
Figure 2. Induced damage in compression -comparison with experimental results obtained for an initial damage equal to zero.

Figure 4 .
Figure 4. SMART mock mesh : concrete and reinforcement

Table 2 .
. Modal analysis of the SMART mock-up Damage prediction of the SMART mock-up of different level of earthquake loading PGA : 0.1 g PGA : 0.4 g PGA : 0.8g