Implementation of a damage model in a finite element program for computation of structures under dynamic loading

This work is a numerical simulation of nonlinear problems of the damage process and fracture of quasi-brittle materials especially concrete. In this study, we model the macroscopic behavior of concrete material, taking into account the phenomenon of damage. J. Mazars model whose principle is based on damage mechanics has been implemented in a finite element program written Fortran 90, it takes into account the dissymmetry of concrete behavior in tension and in compression, this model takes into account the cracking tensile and rupture in compression. It is a model that is commonly used for static and pseudo-static systems, but in this work, it was used in the dynamic case.


INTRODUCTION
he numerical program is employed in analyzing Koyna gravity dam and comparing the results with available results of many researchers find in literature.The choice of Koyna concrete gravity dam is due to the fact that many studies have been done on this structure which experienced a devastating earthquake in 1967.Since Koyna dam in the 1967 Koyna earthquake [1][2][3][4] is one of the few examples of cracking in a concrete dam, it has been selected to investigate its earthquake response.The earthquake records are in the form of accelerograms with horizontal and vertical components which are used as dynamic loads.The effect of the damping coefficient is obviously taken into account.This work proposes the implementation of a damage model based on the damage mechanics as part of the isotropic formulation proposed by J. Lemaitre in a finite element program.the J. Mazars model [5] whose principle is based on damage mechanics which is a theory describing the progressive reduction of the mechanical properties of a material due to initiation, growth and coalescence of microscopic cracks.These internal changes lead to the degradation of mechanical properties of the material.The model takes into account the asymmetry of concrete behavior and it considers the cracking in tension and compression failure.The method is based on the notion of equivalence of the strains.
T MODEL he influence of microcracking due to external loads is introduced via a single scalar damage variable d ranging from 0 for the undamaged material to 1 for completely damaged material.The stress-strain relation reads [6]: 0 E and 0  are the Young's modulus and the Poisson's ratio of the undamaged material; ij ε and ij σ are the strain and stress components, and ij δ is the Kronecker symbol.The elastic (i.e., free) energy per unit mass of material is Where 0 ijkl C is the stiffness of the undamaged material.This energy is assumed to be the state potential.The damage energy release rate is With the energy of dissipated energy: Since the dissipation of energy ought to be positive or zero, the damage rate is constrained to the same inequality because the damage energy release rate is always positive.

DAMAGE EVOLUTION
he evolution of damage is based on the amount of extension that the material is experiencing during the mechanical loading.An equivalent strain is defined as Where . is the Macauley bracket and i  are the principal strains.The loading function of damage is Where  is the threshold of damage growth.Initially, its value is 0  , which can be related to the peak stress t f of the material in uniaxial tension: In the course of loading  assumes the maximum value of the equivalent strain ever reached during the loading history.
The function   h  is detailed as follows: in order to capture the differences of mechanical responses of the material in tension and in compression, the damage variable is split into two parts: Where Note that in these expressions, strains labeled with a single indicia are principal strains.In uniaxial tension 1 The evolution of damage is provided in an integrated form, as a function of the variable  : Figure 1.Evolution of two parts of damage t d and c d [5].
A direct tensile test or three point bend test can provide the parameters which are related to damage in tension ( 0 ).Note that Eq. 7 provides a first approximation of the initial threshold of damage, and the tensile strength of the material can be deduced from the compressive strength according to standard code formulas.The parameters ( c A , c B ) are fitted from the response of the material to uniaxial compression.

DYNAMIC EQUATIONS OF MOTION
ynamic analysis of structures exhibiting non linear behavior is performed by using direct integration, to trace the response in the time domain.The nonlinear dynamic equilibrium equation can be written as where M and C are the global mass and damping matrices respectively, n p is the global vector of internal resisting nodal forces, n f is the vector of consistent nodal forces for the applied body and surfaces traction forces grouped together, the body force term ( g MIu   ) due to seismic excitation, is included in the body forces which are taken into account in n f , n u  is the global vector of nodal accelerations and n u  is the global vector of nodal velocities [7].When the structure is subjected to seismic excitation, the external applied body forces is Where g u  is ground acceleration and I is a vector indicating the direction of the earthquake excitation.

D
The geometry of a typical non-overflow monolith of the Koyna dam is illustrated in Fig. 4. The monolith is 103 m high and 71 m wide at its base.The upstream wall of the monolith is assumed to be straight and vertical, which is slightly different from the real configuration.The depth of the reservoir at the time of the earthquake is h w = 91.75 m.Following the work of other investigators, we consider a two-dimensional analysis of the non-overflow monolith assuming plane stress conditions.The finite element mesh used for the analysis is shown in Fig. 4. It consists of 760 first-order, reducedintegration, plane stress elements (CPS4R).Nodal definitions are referred to a global rectangular coordinate system centered at the lower left corner of the dam, with the vertical y-axis pointing in the upward direction and the horizontal xaxis pointing in the downstream direction.The transverse and vertical components of the ground accelerations recorded during the Koyna earthquake are shown in Fig. 3. (units of g = 9.81 m sec -2 ).Prior to the earthquake excitation, the dam is subjected to gravity loading due to its self-weight and to the hydrostatic pressure of the reservoir on the upstream wall [8].

RESULTS
e note that the displacements are relatively low during the first two seconds because of low the amplitudes of the excitations.The displacements reach their maximum at 3.7 s and 7.5 s, 30 mm was recorded at 3.8 s, the maximum displacement value does not correspond to the maximum amplitude of the excitement that is recorded at 3.65 s.The nodal displacements decrease after 7.5 s.W Figure 5: Horizontal displacements at the crest of the dam We note that the displacements are relatively low during the first two seconds because of low the amplitudes of the excitations.The displacements reach their maximum at 3.7 s and 7.5 s, 30 mm was recorded at 3.8 s, the maximum displacement value does not correspond to the maximum amplitude of the excitement that is recorded at 3.65 s.The nodal displacements decrease after 7.5 s.Fig. 6.Shows the evolution of t D and c D , the tensile and compressive part of the damage variable D .We note that traction damage is more significant and more predominant than in compression, which leads to a brittle fracture of the material by traction.The total rupture of the material can be achieved by traction, in compression the damage is less significant, so the rupture will not occur under compression, Fig. 1b presents some ductility of the material.It is also seen from this figure (Fig. 6) that the tensile damage starts firstly and evolves quickly to a value close to 1 ( 1 t D  ), then followed by the compression damage, but the temporary difference is only of the order of a thousandth of a second, so nearly or completely insignificant difference.In counterparty, compressive damage is less important, it increases until a value 0.15 c D  a value that does not allow rupture in compression of the element.Damage were observed (see Fig. 7) in the dam after 3.53 s at the integration point (815) of the element 204, then from 3.54 s at integration point (811) of the element 203 and finally, from 3.54 s at Gauss point (807) of the element 202.It may be noted that the evolution of the damage is mainly concentrated in the time interval where the maximum values of positive and negative displacements occur.CONCLUSION e have taken the example of a concrete gravity dam subjected to seismic excitation in the form of accelerograms.The results are, first the response of gravity dam subjected to seismic loading as displacements history, secondly the history of the evolution of integration points damage at the cracked or damaged zone.The dynamic equilibrium equation is solved, an integration algorithm by the method of central differences for nonlinear systems, we have used a very small time step (0.0001 s).The choice of weight Koyna dam in India is the fact that many studies have been conducted on this structure to a known destructive earthquake in 1967.The records of the earthquake accelerograms as with horizontal and vertical components were used as dynamic loads.
Hence, t d and t d can be obtained separately from uniaxial tests.

Figure 2 :
Figure 2: Uniaxial response of the model in (a) tension and (b) compression[5]

Figure 6 :
Figure 6: Damage evolution, t D and c D .