An isotropic damage model to simulate collapse in reinforced concrete elements 1

The simulation of structures in their evolution to collapse is an interesting topic in engineering to know the remaining capacity of a structure after reaching the ultimate load, which is dependent on the constitutive behaviour of their construction materials. After concrete reaches its ultimate load, it presents a reduction of its load-carrying capacity, so the stresses decreases while the strains increases; this behaviour is typically called strains softening. The continuum damage mechanics can be used for modelling concrete failure, one of the pioneers of damage mechanics theory was Kachanov (1958), who introduced the fundamental concepts of the continuous damage mechanics, i.e., effective strain, effective stress and scalar damage variable. There are many local damage models based on the damage framework such as: isotropic damage models (Kachanov 1958, Krajcinovic 1984, 1985, Lemaitre 1986, Simo and Ju 1987a, Lubliner et al. 1989, Mazars and Pijaudier-Cabot 1989, Lubarda et al. 1994, Faria et al. 1998, Lee and Fenves Juárez-Luna G. a,* Méndez-Martínez H.a,1 Ruiz-Sandoval M.E.a,2


INTRODUCTION
The simulation of structures in their evolution to collapse is an interesting topic in engineering to know the remaining capacity of a structure after reaching the ultimate load, which is dependent on the constitutive behaviour of their construction materials.After concrete reaches its ultimate load, it presents a reduction of its load-carrying capacity, so the stresses decreases while the strains increases; this behaviour is typically called strains softening.
The continuum damage mechanics can be used for modelling concrete failure, one of the pioneers of damage mechanics theory was Kachanov (1958), who introduced the fundamental concepts of the continuous damage mechanics, i.e., effective strain, effective stress and scalar damage variable.There are many local damage models based on the damage framework such as: isotropic damage models (Kachanov 1958, Krajcinovic 1984, 1985, Lemaitre 1986, Simo and Ju 1987a, Lubliner et al. 1989, Mazars and Pijaudier-Cabot 1989, Lubarda et al. 1994, Faria et al. 1998 1998, Comi 2001, Jason et al. 2006, Contrafatto and Cuomo 2006, Wu et al. 2006and Gopinath et al. 2012); anisotropic damage models (Dragon and Mroz 1979, Cordebois and Sidoroff 1979, Krajcinovic and Fonseka 1981, Ortiz 1985, Simo and Ju 1987a, b, Ju 1989, 1990, Voyiadjis and Abu-Lebdeh 1994, Govindjee et al. 1995, Halm and Dragon 1996, Fichant et al. 1999, Lemaitre et al. 2000, Carol et al. 2001and Cicekli et al. 2007).It is well kwon that local continuum damage models lead to a physically unrealistic description of the strain localization phenomenon.When the tangent stiffness matrix is not positive definite, the governing differential equations may lose ellipticity, which renders the boundary value problem ill posed (Needleman 1988, Jirásek 1998).In finite elements computations spurious mesh sensitivity problems occur, strain localizes into a narrow band whose width depends on the element size, the corresponding load-displacement curve always shows snapback for a fine mesh and the dissipate fracture energy converges to zero.To overcome the shortcomings of local theories for modelling strain softening, some alternatives have been proposed such as: non-local modelling of the constitutive behaviour (Costa Mattos et al. 1992, Costa Mattos and Sampaio 1995, Da Costa Mattos et al. 2009), gradient material description (Triantafylidis andAifantis 1986, Peerlings et al. 1996), micropolar continuum theory (de Borst 1991, Steinmann 1995), viscous regularization (de Borst and Mühlhaus 1992) and local manipulations of the material properties depending on the element size (Bazant and Oh 1983, Crisfield 1986, Oliver 1989).Recently, Da Costa- Mattos et al. (2009) proposed a special gradient-enhanced damage theory in which the material is considered to possess a substructure or microstructure, where the free energy is supposed to depend not only on the strain and the damage variable but also on the damage gradient.
If the concrete is modelled with the finite element method and local damage models, the damage of a cracking surface is distributed into the material volume of a finite element, which is inconsistent with the fracture energy dissipated per unit of area and per unit of volume.This disagreement is overcome with a characteristic length such as the proposed one by Bazant and Oh (1983), who used a thickness of the cracking band, depending on the element area and on the propagation direction; Crisfield (1986) proposed a characteristic length defined by the Jacobian in each integration point; and Oliver (1989) developed a general method for computing the characteristic length depending on the element size and the elastic stress state.
The Software Ansys 12.0 uses a damage surface for concrete with different threshold values under tension and compression, coherent with the behaviour reported in laboratory test such as Kupfer and Gerstle (1973).However, this model has as a drawback that fracture energy has a fixed value, which is not possible to modify, so the results are dependent on that fixed energy fracture.Also, Ansys only considers linear softening curve, even though there are other curves with a better approximation such as bilinear, trilinear and exponential.Alternatively, the Software Diana 9.0 uses the smear cracking model for concrete (Rashid 1968, Rots 1988), in which is possible to assign the value of the fracture energy density.Nevertheless, this model has problems of mesh dependency, stress locking and spurious kinematic modes.
In the attempt to simulate material failure by mesh independent models, Hillerborg et al. (1976) developed the discrete crack model, where a previous stress analysis must be performed to a priori localize the strain concentration zones, where interfaces are placed on the boundaries of the elements, bringing about that the results are mesh dependent.Another alternative to simulate materi-Latin American Journal of Solids and Structures 11 (2014) 2444-2459 al failure is the strain localization model (Simo et al. 1993, Oliver 1996a, 1996b, Armero and Garikipati 1996, Juárez and Ayala 2009, Juárez-Luna and Ayala 2014), which considers that only one discontinuity (crack) may propagate into a finite element, so in the case that multi cracks appears, the computational cost of their propagation is high.
For the aforementioned reasons, this paper promotes the use of continuous damage models, which do not need computational costly tracking algorithms.A continuous damage model is developed and implemented to simulate the constitutive behaviour of the concrete in its evolution to collapse.This damage model, called different tension-compression (DTC), has a failure surface with the corresponding threshold values in tension and compression of the concrete.Additionally, in this model, it is possible to assign the density energy magnitude and it does not have the problems reported by Rots (1988).
The proposed DTC damage model uses the failure surface proposed by Oliver et al. (1990).Nevertheless, the tangent constitutive tensor is different to the one developed by Linero (2006).The DTC damage model was implemented into the subroutines of the Finite Element Analysis Program (FEAP) (Taylor 2008), then it was validated with the numerical reproduction of numerical examples and experimental reinforce concrete tests reported in the literature.It is important to mention that in this paper, the material failure is represented by a isotropic damage variable, which represents the degradation level, so the material degradation into the finite element analyses is presented as localization zones where important displacements occurs, but not as physical discontinuities (cracks).

Equal tension-compression
The isotropic continuum damage model by Simo and Ju (1987a) and Oliver et al. (2002) with a failure surface equal tension-compression (ETC) is shown in Figure 1 for 1D and 2D problems, where the ultimate tensile and compressive strength, σ ut , and, σ uc , respectively, have the same magnitude.This model, typically used to simulate the constitutive behaviour of some metals, is defined by the following equations: where Ψ is the free energy density, ε is the strain tensor, σ is the stress tensor and C is the elastic tensor.The damage variable, d, is defined in terms of hardening/softening variable q, which is dependent on the hardening/softening parameter, H.The damage multiplier γ determine the loadingunloading condition, the function, f(τ σ ,q), bounds the elastic domain defining the damage surface in the stress space.The value, r o , is the threshold that determines the limit of the initial elastic domain, i.e, q=r o .In this ETC damage surface, any stress state is transformed to a norm, bounded by r o , where every stress state outside of this circle is inelastic as shown in Figure 2.  The tangent constitutive equation, in terms of rates, from the model in Eq. ( 1) is:

 
Where C T is the continuum tangent constitutive operator, relating the stress and the strains rates, of the nonlinear loading interval, which is defined by ( )

Different tension-compression
The formulation of the damage model with failure surface DTC is based on the ETC damage model.The main difference of this model is the definition of the norm, which is modified for a parameter proposed by Lubliner et al. (1989).It is well known that concrete has a ratio between the compression and tension strength of 10 to 20, i.e., n=σ uc /σ ut , as shown in Figure 3.To have a failure surface with that characteristic, the norm given in Eq. ( 1)e, is modified by the parameter χ, then Also, the tangent constitutive tensor, C T , for the nonlinear interval is: The parameter  in Eqs. ( 5) where ϕ a is weight factor, depending on the principal stresses, σ i , given by: Latin American Journal of Solids and Structures 11 (2014) 2444-2459 Where the Maclauy 〈〉 operator and the symbol || consider, respectively, the positive and the absolute magnitudes of the principal stresses.The interval of ϕ is [0,1], bounded by 0 for a triaxial compression (σ 3 ≤ σ 2 ≤ σ 1 ≤ 0) and 1 for triaxial tension (0 ≤ σ 3 ≤ σ 2 ≤ σ 1 ).Consequently, the corresponding interval of χ is [1/n,1], bounded by 1/n for a triaxial compresion and 1 for triaxial tension.
The parameter χ scales down 1/n times the norm, as shown in Figure 4, in such a way it is compared with the elastic interval [0,r o ].The value of r o =σ ut / depends on the threshold value of σ ut and the Young's modulus E. Note that the initial elastic interval is the same for 1D, 2D and 3D problems, all with a limit point r o , and that for a 2D principal stress state, the parameter χ takes the value of 1 in the firts quadrant, 1/n in the third quadrant and the interval [1/n,1] in the second and fourth quadrants.The algorithm to implement the DTC damage model is as follows: a) If τ t+1 <r t , the elastic/unload/neutral load case, then ( ) 4. Output values: r t , C T and σ t+1 .Note that if the value of χ=1 in the norm and in the tangent constitutive tensor, the ETC continuous damage model is recovered.An equivalent way of evaluate the damage criterion is proposed, modifying the threshold value of r o in the limit point for the value of r om =r o /χ.With this modification, in a 2D principal stress state, the value of r o is scaled up 1 time in the first quadrant, n times in the third quadrant and 1 to n times in the second and fourth quadrants, as shown in Figure 5.Then, the algorithm for the ETC continuous damage needs only a modification in the threshold value of r o .Both ways of the DTC continuous damage model were implemented in the user routines of FEAP, having the same results.

NUMERICAL EXAMPLES
In the three presented examples, the reinforcement was meshed with 3D linear elements with two nodes equipped with three degrees of freedom.The constitutive behaviour of the steel reinforcement was modelled with a plasticity surface of von Mises.The steel elements were placed on the boundary of the solid elements, coupling the degree of freedom of both kinds of elements and then perfect bond was considered.This DTC damage model was validated for plain concrete elements in Méndez-Martínez and Juárez-Luna (2012).

Lightly reinforced concrete beam
In this example, a simple supported concrete beam is analysed, the dimensions of the reinforcement and boundary conditions are shown in Figure 6.A load is applied at the centre of the span by displacement control, which is gradually imposed downwards.The beam was meshed with hexahedral elements as shown in Figure 7a, where only a quarter of the beam was modelled because of its symmetry, reducing the time consuming of the computational cost.As can be seen, the failure started into the finite elements placed around the centre on the span, at the bottom of the beam as shown in Figure 7b, which is coherent with the experimental results reported by Burns and Seiss (1962).The load vs. displacement curves at the centre of the span are shown in Figure 8, where the curve computed with the DTC damage model is acceptable with the curve reported by Burns and Seiss (1962), showing the validity of the implemented DTC damage model to simulate reinforced concrete elements.

Reinforced concrete circular slab
In this example, a reinforced concrete slab with circular shape was studied, in which a distributed uniform load is gradually applied downwards.Two boundary conditions were considered at the border of the slab, simple supported and fixed.The geometry and reinforcement of the slab with thickness of 0.10 m are shown in Figure 9.As the previous example, the two axes or symmetry were considered, so only a quarter of the slab was modelled with hexahedral elements, as shown in Figure 10a, to reduce the time consuming of the computational calculation.The undeformed meshes are shown in Figure 10b and c, for the simple supported and fixed at the border of the slab, respectively.The load vs. displacements curves are shown in Figure 11a and b for simple supported and fixed conditions, respectively.These curves are compared with those reported by Caballero (2012), computed by a smear crack model in the Software ANSYS.As can be seen, both computed curves are coherent with those simulated with a smear crack model.

Reinforced concrete deep beam
In this example, a reinforced concreted deep beam (l/h=1) is analysed, where a distributed uniformed load, with magnitude q=0.0078125P, is gradually applied downwards on the top of the beam shown in Figure 12a.In this example, an elastoplastic constitutive model was assigned to the steel as Leonhardt and Walther (1966) did.Taking advantage of the symmetry, only a quarter of the beam was modelled to reduce the time consuming computational cost.Figure 13 shows the 3D mesh, discretized with hexahedral elements.Figure 14 shows a comparative between the computed load-displacement curve and the experimental curve reported by Leonhardt and Walther (1966).It is observed that the numerical results are congruent whit the experimental results.

CONCLUSIONS
A damage model to simulate concrete elements to collapse was developed, implemented and validated, based on the damage mechanics continuous theory.The damage model considers a constitutive behaviour equipped with softening for the concrete elements limited for a damage surface, whereas for the reinforced concreted elements, a plasticity model was assigned to the reinforcement, which take into account the hardening of the steel with a Von Mises surface.
A DTC damage surface for concrete is obtained by scaling down the norm and the tangent modulus with the parameter χ, but it is equivalent by scaling up the threshold value of r o .In both cases, the integration algorithm of the ETC damage model needs two and one modification, respectively, to have a surface with different threshold values in tension and compression.
The global behaviour of reinforced concrete elements is controlled by the hardening of the reinforcement steel.Their simulation with damage models for concrete and plasticity model for the steel provides the cracking zones, the crack paths and the ultimate load of the specimen.
The numerical results computed with the DTC damage model are coherent with the experimental and numerical results reported in the scientific literature, since it properly simulates the strain softening of the concrete and localizes the degradation zones.
To consider perfect bond between steel reinforcement and concrete does not show evident effects in the presented numerical examples.
It is useful over some commercial software that it is possible to assign the fracture energy density value in the implemented damage model, because the numerical results are dependent on this value.
The implemented damage model did not show the stress locking problem, a fact that guaranties a suitable simulation of the plain and reinforced concrete elements in their evolution to collapse.

Figure 2 :
Figure 2: Transformation of the stress state to a norm.

Figure 4 :
Figure 4: Transformation of the stresses to a norm.

Figure 5 :
Figure 5: Transformation of the stresses to a norm with a modified r 0 .

Figure 11 :
Figure 11: Load vs. displacement curves at the circular slab: a) simple supported and b) fixed.
The mechanical properties of the concrete are: Young's modulus E c =14710 MPa, Poisson ratio ν=0.2, ultimate tensile strength σ tu =2.45 MPa, ultimate compressive strength σ cu =29.42 M Pa fracture energy density G f =0.098 N/mm.The mechanical properties of the steel reinforcement are: Young's modulus E s =196.13MPa, yield stress σ y =411.88MPa, Poisson ratio ν=0.3, area A=0.71 cm 2 and hardening modulus H=2.871 GPa.
The mechanical properties of the concrete are: Young's modulus E c =31400 MPa, Poisson ratio ν=0.2, ultimate tensile strength σ tu =2.5 MPa, ultimate compressive strength σ tu =29.63 MPa and fracture energy density G f =0.098 N/mm.The mechanical properties of the steel reinforcement are: Young's modulus E s =206 GPa, yield stress σ y =536.6 MPa, Poisson ratio ν=0.3, area of the principal reinforcement bars A 1 =0.503 cm 2 and the secondary reinforcement bars A 2 =0.178 cm 2 and hardening modulus H=2.871 GPa.The distribution of the reinforcement is shown in Figure 12b.