Computational Modelling of RC Slabs Cracking with an Embedded Discontinuity Formulation 1

The study of collapse in concrete elements is an interest topic in engineering; particularly, the determination of the first crack load and the crack pattern. It is well known that concrete strength under compression is from 10 to 20 times greater than the strength under tension, as shown the experimental results reported by Kupfer and Gerstle (1973). In the way to collapse of reinforced concrete elements, their behaviour at the beginning is approximately linear elastic; next, cracking occurs. Then, crushing appears and finally, plasticity in reinforcing steel initiates, although strucJuárez-Luna G. a, 1

tural collapse may occur before yielding of steel bars. Particularly, in clamped reinforced concrete slabs, cracking initiates on the top surface, then at the centre on the bottom surface, growing as the load increases; whereas in simple supported slabs, cracking initiates at the centre of the span on the bottom surface, growing to the edges (Juárez-Luna and Caballero-Garatachea 2014).
Laboratory tests have been perform to obtain the cracking paths and the moment coefficients of the rectangular slabs such as Bach and Graf (1915), who tested 52 simple supported slabs on their edges and 35 strips supported as beams, which were loaded until failure occurs. In these experimental tests, the displacements at some points of the slabs and the slopes at the centre of the edges were measured; also, the propagation of the cracks was registered as reported by Westergaard and Slater (1921). Afterwards, other experimental tests were performed as those reported by Casadei et al. (2005), Foster et al. (2004), Galati et al. (2008), Gamble et al. (1961), Girolami et al. (1970), Hatcher et al. (1960), Hatcher et al. (1961), Jirsa et al. (1962), Mayes et al. (1959), Vanderbilt (1961), among others. It is interested to say that these references are the basis of the research and applications of the current analysis and design of the rectangular slabs.
In the modelling of the reinforced concrete slabs, de Borst and Nauta (1985) applied the smeared crack model to study an axisymmetric slab under shear penetration, showing that cracking initiated at the bottom face of the slab and the corresponding cracking paths. Then, Kwak and Filippou (1990) modelled a square slab supported on its corners with a concentrated load at the centre of the span, obtained the load vs. displacement curve, which was congruent with experimental results reported by Jofriet and McNeice (1971) and Mcneice (1967); in the reported results by Kwak and Filippou (1990), neither the first crack load nor the cracking pattern was given. There were other proposals for modelling reinforced concrete slabs such as Gilbert and Warner (1978), Hand et al. (1973), Hinton et al. (1981), Lin and Scordelis (1975), Wang et al. (2013) among others, most of them used the smeared crack model.
There are some commercial software for modelling reinforced concrete elements such as ABAQUS (ABAQUS 2011), ANSYS (ANSYS 2010), DIANA (DIANA 2008), ATENA (Kabele et al. 2010), NLFEAS (Smadi and Belakhdar 2007), among others. These software mainly use the finite element method with the smeared crack model for the behaviour of the concrete, equipped with a failure surface with different threshold value in tension and compression, necessary to determine the first crack load and crack propagation. However, the smeared crack model may have numerical problems of stress locking and spurious kinematic modes (Rots 1988), which may be overcome with heuristic shear retention factors.
In this paper, finite elements with embedded discontinuities (FEED) were used for studying reinforced concrete slabs, computing their load-displacement capacity curves and their cracking patterns. The advantages of FEED are the capability for representing highly localized strains by improving the kinematic, the possibility to statically condense out the displacement jump and the nearly mesh-independent. Concrete was discretized with hexahedral FEED and steel reinforcement was discretized with 3D bar elements, both kinds of elements have three degrees of freedom per node.
The outline of this paper is as follows. Section 2 presents the details of FEED formulation. Section 3 provides the constitutive models to describe the behaviour of the materials, a discrete damage model equipped with softening for concrete and a plasticity model for the steel reinforcement. Numerical examples of reinforce concrete slabs which validate the proposed formulation are presented in Section 4. Finally, in Section 5, conclusions derived from this work are given.

Variational formulation
The FEED are formulated from an energy functional which has the displacement, u, and the displacement jump [|u|] as independent variable (Alfaiate et al. 2003, Juárez and Ayala 2009, Lotfi and Shing 1995, Wells and Sluys 2001. This functional is given by: where the free energy density, ( ), u Ψ ε ε ε ε depends on the continuous strain field u ε ε ε ε , and the free discrete energy density, u ( ), S   Ψ   depends on the displacement jump. These energy densities are respectively given by: where the elastic stresses, σ σ σ σ , are defined by: and T S is the traction vector at the discontinuity.

Finite Element Approximation
It is not possible to prescribe the boundary conditions, u * , in only one of the displacement fields, i.e., u or [|u|], a difficulty overcame, according to Oliver (1996), defining the displacement as in Eq. (5), shown in Figures 1a and b: Then, the strain field is defined by: where û is the regular displacement field and x ( ) S M is a function given by: is a continuous function such that: The function S M , has two properties:  The continuous displacement field is defined as: In the continuous part of the solid, which may be linear elastic, the continuous strain field, ε ε ε ε is given by: Substituting Eq. (9) into Eq. (10), If the displacement jump is constant in Eq. (11), the continuous strain field may be rewritten as:

Approximation of the Displacement and Strain Fields
The regular displacement field is approximated by: ˆ= u Nd (13) where N is the standard vector of shape functions of the element and, d, is the nodal displacement vector. The function, M S (x), is defined in the finite element approximation as: where e φ is constructed by: The continuous strain field in Eq. (12) is approximated as: where B , is the standard strain interpolation matrix, containing the derivatives of the standard shape functions Nd The equilibrium equations corresponding to this formulation are obtained by substituting Eqs.
(17) and (18) In Eqs. (19) and ( where R has the direction cosines, R 1 and R 2 are defined as: To reduce the size of the system given in Eq. (21), the additional degrees of freedom, Δ[|u|], may be condensed. In Eq. (22), R 1 means the equilibrium between the external and the internal forces in the domain \ S Ω , whereas R 2 , in Eq. (23), the equilibrium between the forces in the domain \ S Ω and forces in the discontinuity S Γ .
Tractions at the discontinuity are: which expressed in the local system becomes The definitions of the traction vector in Eqs. (24) and (25) are dependent on the discontinuity area, A d , and the direction cosines to the normal vector n, as shown in Figure 2. The FEED, given in eq. (21), were implemented in the finite element analysis program (FEAP), developed by Taylor (2008). These FEED capture a discontinuity surface at their geometric centre, which is placed perpendicular to the major principal stress direction. The discontinuity sur- faces of the surrounding elements are not aligned as Sancho et al. (2007) did for 2D problems. Although there are formulations which consider linear displacement jumps inside the element such as Alfaiate at al. (2003) and Contrafatto et al. (2012), the displacement jump is constant into these FEED. It is important to say that these elements do not have problems of spurious shear deformations and they satisfied the following requirements: (1) equilibrium, traction continuity across the discontinuity interface and (2) kinematics, free relative rigid body motions of the two portions of an element split up by a discontinuity (Juárez-Luna and Ayala 2014).

Concrete
The concrete behaviour was modelled with a discrete damage model, which has different threshold values under tension and compression, as shown in Figure 3a. This model is equipped with softening after reaching the ultimate tensile strength, T ut , or the ultimate compressive strength, T uc , shown in Figure 3b. This model is defined by the following equations: Hardening rule ; 0 Loading-unloading 0; 0; 0; 0 (consistency) conditions whereφ is the discrete free energy density, T is the traction vector. The damage variableω is defined in terms of the hardening/softening variable _ q ɺ , which is dependent on the hardening/softening parameter. The damage multiplier λ determines the loading-unloading conditions, the function ( f T , ) q , bounds the elastic domain defining the damage surface in the tractions space. The tangent constitutive equation, in terms of rates from the model in Eq. (26), is: where C T d is the tangent constitutive operator, relating the traction and the displacement jump of the nonlinear loading interval, which is defined by This constitutive model considers a material fictitious interpenetration, [|u|] c , when T uc is overcome. Because of the energy dissipation in compression takes place over a surface, rather than within a volume, the use of a fictitious interpenetration instead of a strain is justified (Carpinteri et al. 2010). In this model, it is assumed that both post-peak regimes, tension and compression, have decreasing functions with the same slope, EH, as ca be seen in Figure 3b.
The area under the traction versus jump displacement curve represents the fracture energy density, G f , in the first quadrant of Figure 3b, i.e.: where [|u|] t cr is the critical value of the crack opening, which is given by: Substituting Eq. (31) into Eq. (30) and solving for EH, the slope of the decreasing function under tension is given by: Considering that the area under the compression versus the fictitious jump interpenetration displacement curve represents the crushing energy density, G c , in the third quadrant of Figure 3b, then The fictitious critical value of interpenetration can be computed as: Substituting Eq. (34) into Eq. (33), then: As the slope of the decreasing function under tension is the same under compression, Eq. (32) is substituted into Eq. (35): This equation provides a relationship between the crushing and the fracture energy densities.

Steel
A 1D rate independent plasticity model with isotropic hardening was used for modelling the steel reinforcement. This model has the same threshold value in tension and compression, as shown in Figure 4a; the hardening of the steel reinforcement, after reaching the yield stress σ y , was considered with an idealized bilinear function as shown in Figure 4b. The plasticity model is defined by the following equations: Loading-unloading 0; 0; 0; 0 (consistency) conditions where ψ is the free energy density, σ is the stress tensor. The plastic variable α is defined in terms of the hardening variable H. The plastic multiplierλ determines the loading-unloading con- ditions, the function f(σ,σ y ), bounds the elastic domain defining the plastic surface in the stress space. The tangent constitutive equation, in terms of rates from the model in Eq. (38), is: where C T is the tangent constitutive operator, relating the stresses and the strain of the nonlinear loading interval, which is defined by

NUMERICAL EXAMPLES
In the presented examples, the reinforcement was meshed with 3D linear finite elements with two nodes, which have three degrees of freedom each. The constitutive behaviour of the steel reinforcement was modelled with a plasticity model with hardening. The steel elements were placed on the edges of the solid elements, coupling the degree of freedom of both kinds of elements. Perfect bond between steel bars and concrete was assumed, as the failure of this type of slabs occurs mainly on flexure without evidence of debonding.

Square Slabs Supported at the Corners
The FEED and the constitutive models were validated by the numerical modelling of the experimental results reported by Girolami et al. (1970). The test specimen, shown in Figure 5a, is a square slab of sides 1.829 m long and a thickness 0.044 m, which was simple supported at its cor- ners. The top and bottom reinforcement used in the slab consisted of 3.66 mm diameter steel bars cut from No.7 gage wire, which were spaced 10.954 cm in both orthogonal directions, as shown in Figure 6a. Also, the stirrups were bent from the No.7 gage steel wire as shown in Figure 6b. The vertical loads were applied at the top of the slab using four jacks with loading trees distributed to 16 load plates, as shown in Figure 5b. The mechanical properties of the con-

Latin American Journal of Solids and Structures
Only one quarter of the slab was modelled, considering the two axes of symmetry of its geometry. To describe the symmetry conditions at the boundary constraints, no translation was imposed on the degree of freedoms perpendicular to the corresponding axes of symmet ment mesh is shown in Figure 7a, whereas the concrete The load vs displacement at the centre of the span curves are shown in curve with FEED shows numerical results ported by Girolami et al. (1970). In fact, it is observed that both curves are similar at the begi ning; however, there is a backward of the displacement. This movement may be attributed to a sliding of the measurement devices, cons move upwards as the load is applied downwards ments such as beams, where snapback behaviour may occur show that the numerical and the experimental curves are greater than the horizontal curve puted with the yield line theory. The yield line theory is a method to estimate the ultimate load of a slab system by postulating a collapse mechanism compatible with the bound The moments at the plastic hinge lines are the ultimate moments of resistance of the sections, and the ultimate load is determined using librium (Park and Gamble 2000). quarter of the slab was modelled, considering the two axes of symmetry of its geometry. To describe the symmetry conditions at the boundary constraints, no translation was imposed on the degree of freedoms perpendicular to the corresponding axes of symmetry. The steel reinforc a, whereas the concrete mesh of the slab and beams is shown in b) : Meshes: a) steel reinforcement in the slab and b) concrete.
displacement at the centre of the span curves are shown in Figure 8. The computed numerical results, which are congruent with the experimental curve . In fact, it is observed that both curves are similar at the begi ning; however, there is a backward of the displacement. This movement may be attributed to a , considering that the displacement at the centre would never the load is applied downwards; this effect only happens at plain concrete el where snapback behaviour may occur (Carpinteri 1988). It is important to numerical and the experimental curves are greater than the horizontal curve co The yield line theory is a method to estimate the ultimate load of a slab system by postulating a collapse mechanism compatible with the boundary conditions. The moments at the plastic hinge lines are the ultimate moments of resistance of the sections, ultimate load is determined using the principle of virtual work or the equations of equ : Comparison between experimental and numerical results. 3.0

d (cm)
Embedded discontinuity Girolami et al. (1970) 1.0 1.5 2.0 2.5 Yield line theory quarter of the slab was modelled, considering the two axes of symmetry of its geometry. To describe the symmetry conditions at the boundary constraints, no translation was imposed on ry. The steel reinforcethe slab and beams is shown in computed curve re-. In fact, it is observed that both curves are similar at the beginning; however, there is a backward of the displacement. This movement may be attributed to a idering that the displacement at the centre would never at plain concrete ele-It is important to com-The yield line theory is a method to estimate the ultimate load ary conditions. The moments at the plastic hinge lines are the ultimate moments of resistance of the sections, the principle of virtual work or the equations of equi-Cracking initiated at the corners on the top surface, growing along the edges to the centre as shown in Figure 9a; on the bottom surface, cracking occurs at the centre, growing to the corner as shown in Figure 9b. The crack patterns at both, top and bottom surfaces, are congruent with the experimental results reported by Girolami et al. (1970).

Simple Supported and Clamped Slabs
Other two slabs were also analysed, a 4m square slab and a 2m width by 4m length rectangular slab respectively shown in Figure 10; both of them were 10 cm thick and subjected to an increasing uniform distributed load. Two boundary conditions were considered along their edges: simply supported and clamped. The reinforcement used in the slabs consisted of 3/8 in diameter steel bars, which were spaced 20 cm in both orthogonal directions, as shown respectively in Figures 11  and 12 area A=0.71 cm 2 and hardening modulus H=2.871 GPa. Like the previous example, the two axes or symmetry were considered, so only a quarter of the slab was modelled, as shown in Figure 13, which reduces the time consuming of the computational calculation.  The distributed load vs. displacement curves at the centre of the spans of both slabs are tively shown in Figure 14. Here, it may be observed that clamped slabs are approximately five times greater than at the displacements of 2.5 and 0.5 cm, respectively supported slabs, the ultimate load computed with the yield line theory is greater than the load computed with FEED; on the contrary, for clamped slabs, the ultimate load computed with the yield line theory is lower than the load computed with FEED. effect is attributed to the steel reinforcement placed stress before the steel reinforcement placed at the is assumed that every bar instantly reaches the yield at the centre of the spans of both slabs are respec-. Here, it may be observed that the intensity of the distributed load on clamped slabs are approximately five times greater than the intensity on simple supported slabs s of 2.5 and 0.5 cm, respectively. Additionally, it is observed that, for simple ltimate load computed with the yield line theory is greater than the load computed with FEED; on the contrary, for clamped slabs, the ultimate load computed with the yield line theory is lower than the load computed with FEED. In simple supported slabs, this placed at the centre of the span reaches the yield stress before the steel reinforcement placed at the edges does; whereas in the yield line theory, it is assumed that every bar instantly reaches the yield stress. clamped square slab, cracking simultaneously initiated along the edges on the top surface, growing to the centre until a ring is formed, as shown in Figure 15a, whereas on the bottom surgrowing to the corners and to the edges, as shown in b. In the simple supported condition, cracking initiated at the corners and at the centre Figure 16b, whereas on top surface, it initiated at the Figure 16a. In the clamped rectangular slab, cracking simultaneously initiated along its long edges on the top surface. Then, cracking appeared along its short edges, growing to the centre until a ring is formed as shown in Figure 17a. On the bottom surface, cracking initiated along a strip at the centre of the span, parallel to the long edges, growing to the corners and to the edges, as shown in Figure 17b. In the simple supported rectangular slab, cracking initiated along a strip at the centre on the bottom surface, parallel to the long edges as shown in Figure 18b; some cracks appeared at the corners on the top surface, as shown in Figure 18a. In general, the cracking paths on the bottom surfaces of the corresponding studied slabs were similar to the paths proposed by the classical yield line theory as well as the experimental results reported by Bach and Graf (1915). To know if the steel reinforcement yields, the yield strain, ε y = σ y /E s =0.0021, was compared with the strains of the steel bars placed perpendicular to the central and edge strips of the slabs, i.e., perpendicular to strips 0-A, 0-a, 0-B, A-C, a-c and b-c, respectively, as shown in Figure 10. The strains of the steel bars were only computed in half of the strips as they are symmetric.
In the clamped square slab, only the central reinforcing top steel bars placed perpendicular to the edge strip yielded as shown in Figure 19b, while the other bars remained linear elastic as shown in Figure 19a. This result is congruent with the occurrence of cracking, which simultaneously initiated along the edges on the top surface. On the contrary, in the simple supported square slab, only the central reinforcing bottom steel bars placed perpendicular to the central strip yielded as shown in Figure 20a, while the other steel bars remained linear elastic as shown in Figure 20b. This result is also congruent with the occurrence of cracking, which initiated at the centre on the bottom surface.
In the clamped rectangular slab, only the central reinforcing top steel bars placed perpendicular to the long edge strip yielded, strip b-c, as shown in Figure 21c, while the other bars remained linear elastic as shown in Figure 21a, b and d. This result is congruent with the occurrence of cracking, which simultaneously initiated along the long edges on the top surface. On the contrary, in the simple supported rectangular slab, only the central reinforcing bottom steel bars placed perpendicular to the long central strip yielded, strip 0-a as shown in Figure 22a, while the other bars remained linear elastic as shown in Figure 22b, c and d. This result is also congruent with the occurrence of cracking, which initiated at the centre on the bottom surface.

CONCLUSIONS
The FEED with a discrete damage model were used for modelling concrete slabs under vertical load, involving their crack pattern and the load-displacement capacity curve. The discrete 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 1D rate independent plasticity model was assigned to the reinforcement, which take into account the hardening of the steel as isotropic.
The coupling of solid FEED and bar elements was validated with experimental results reported in the literature. Although the experimental curve shows a recovery of the displacement, it is considered that it was an effect of a slipping of the measurement equipment, since under an increasing vertical load it is not possible a recovery of the displacement at the centre of the span.
In clamped slabs, cracking initiated on the top surface at the centre of its edges, growing to the corners until a kind of ring is formed. On the bottom surface, cracking occurred at the centre of the span, propagating to the corners until a cross is formed. Whereas, in simple supported a c slabs, cracking initiated at the centre of the span on the bottom surface, propagating to the corners until a cross was formed. On the top surface, incipient cracking occurred at the corners.
The ultimate load computed for simple supported slabs with the yield line theory is greater than the load computed with FEED; on the contrary, for clamped slabs, the ultimate load computed with the yield line theory is lower than the load computed with FEED. In simple supported slabs analysed with FEED, this effect is attributed that the steel reinforcement at the centre of the span reached the yield stress before the steel reinforcement placed at the edges did; whereas with the yield line theory, it is assumed that every bar instantly reaches the yield stress.
In the clamped square slab, the central reinforcing steel bars placed perpendicular to the edges yielded, while in the simple supported square slab the central reinforcing steel bars placed perpendicular to central strip yielded. On the other hand, in the clamped rectangular slab, the reinforcing steel bars placed perpendicular to the long edges yielded, while in the simple supported rectangular slab, the reinforcing steel bars placed perpendicular to long central strip yielded. These results were congruent with the places of the slabs where cracking initiated, corresponding to the places where grater tension stresses occurred.