Damage and plasticity evolution of reinforced concrete beams using laminated Euler-Bernoulli finite elements

abstract: This paper analyses the behaviour of reinforced concrete beams subjected to progressive loads. Although this topic is addressed in recent works through tri-dimensional finite elements, the present study adopts a simpler alternative, using Euler-Bernoulli’s beam finite element. The evolution of the generalized cracking is addressed by Damage Mechanics theory, and the nonlinear constitutive models of concrete and steel are considered, including plasticity and post-peak ‘. The laminated approach permits not only to ascertain the damage evolution but, in practical aspects, makes easy to define vertically symmetrical transversal sections and to distribute the reinforcement as desired, even if the steel rebars have different lengths. It is highlighted the changing of the neutral axis during the damage evolution process, which is updated automatically by the developed numeric-computational code. To validate the proposed methodology, two examples are detailed, both based on experimental tests found in the literature. The results obtained are close to the experimental ones and confirm the applicability of the proposed approach.


INTRODUCTION
For civil engineering applications, the adoption of one-dimensional finite elements in structural analysis is frequent due to its low computational cost.Over the last years, however, the necessity of more complex analyses has been growing, to increase the structure reliability.The pursuance of conducting more advanced studies, however, usually demands more computational resources.Several research take advantage of one-dimensional finite elements, making necessary adjustments, according to the conditions.In reinforced concrete, subjected to cracking and plasticity, or dynamic loadings, many methodologies can be implemented to consider these effects in simulations.
Álvares [1] proposes a simplified methodology to consider the damage and plasticity, which are applied to the element's nodes.Biondi and Caddemi [2] presented a model which uses distributions such as unit step and Dirac's delta functions to simulate discontinuities in the curvature and in the slope functions of Euler-Bernoulli beams.Ghrib et al. [3] presented two computational procedures to recalculate the stiffness due to damage in Euler-Bernoulli beams, identified through static deflection measurements.One of the formulations is based on the principle of the equilibrium gap along with a finite-element discretization, obtaining a solution by minimizing a regularized functional using a Tikhonov total variation scheme.The second formulation refers to a minimization of a data-discrepancy functional between measured and model-based deflections.In order to account for multi-axial coupling of axial force, shear, flexure and torsion in a beam element, Corvec [4] proposes a model with warping degrees of freedom at arbitrary points of the cross-section, in order to accommodate higher order strain kinematics.Regarding dynamic effects, Xie et al. [5] compared results obtained from beam and plane state finite elements in rail-wheel interaction, obtaining similar outcomes on both models.
To consider the physical non-linearity, some research adopt the beam element's discretization into fibers or layers.Spacone et al. [6] present a fibre beam-column element in which the equilibrium state is obtained through an iterative algorithm that satisfies the element constitutive relation.The method is suitable for the simulating reinforced concrete columns under varying axial load, subjected to non-linear hysteretic behaviour of softening members.Mazars et al. [7] presented a multifiber beam element capable of reproducing shear from Timoshenko's theory or shear due to torsion.Oliveira et al. [8] studied the bond-slip effects on the steel-concrete interface of reinforced concrete beams.For such, is presented a beam-layered model, which accounts the materials interaction.Lezgy-Nazargah et al. [9] proposed a refined global-local laminate theory, based on the superposition hypothesis for bending and vibration analyses of the laminated/sandwich composite beams.
By dividing a transversal section into layers, the following advantages may be highlighted: a) each concrete layer may have its own stress state and may be damaged differently; b) any type of vertically symmetrical transversal section can be adopted; c) the reinforcement bars can be distributed as desired in any position of cross-section; and d) the verification of the neutral axis can be updated for each step.
The laminated beam theory for Euler-Bernoulli model is adopted in the present work.Details about this methodology, based on the work of Abeche et al. [10], are addressed.The concrete is taken as bi-modular material, with different behaviour for tension, compression and both combined with shear.The damage criterion in this work is the Mazars' damage model [11] and it is applied for each layer in all Gaussian points of the elements.If the damage occurs, the neutral axis is updated since the cross-section's stiffness is changing.An elastic perfectly plastic stressstrain constitutive model is adopted for the steel reinforcement, with same behaviour for tension and compression.The shear stress and the consideration of stirrups are also addressed in this work.The bonding between concrete and steel rebars is assumed perfect without slipping.
To validate this approach, two examples are presented.The first one is a single supported beam which was modeled and tested by Mazars and Grange [12].The second problem is related to the famous work of Bresler and Scordelis [13], who tested many reinforced concrete beams with different cross-sections and for many kinds of reinforcements.Many models were developed [14], [15], [16] aiming to reproduce the experimental results of Bresler and Scordelis [13].
The present work developed a finite element code in C++ programming language wrapped with MATLAB.It can solve linear and nonlinear static models.For progressive loads, a multi-step scheme is adopted with the Newton-Raphson iterative technique for the nonlinear problems.
This work is divided into seven more items following this Introduction.A brief review of Euler-Bernoulli beam theory is presented in the subsequent section, which also brings the laminated beam theory.The following section addresses the constitutive models for concrete and steel.The damage model is discussed, highlighting the Mazars' model.Then, some aspects of the variation of neutral axis are pointed out, followed by an item that presents the scheme for nonlinear static equilibrium.The next section details the two examples analysed in this work, and the results are presented.The last item is reserved for discussion and conclusions.

The classical EB beam element
The classical Euler-Bernoulli beam finite element is presented for modelling reinforced concrete beams.The main hypotheses are presented, where  1 and  2 are, respectively, the horizontal and vertical axes: 1. the existence of a neutral axis; 2. plane sections stay plane after the bending; 3. the material is linear elastic and homogeneous; 4. normal and shear stresses,   2 and   1  2 , are too small in relation of   1 and may be ignored; 5. the  1  2 is a principal plane.
The differential equation that represents the problem is (Equation 1) Where  is the Young's modulus,  is the beam's moment of inertia,  is the external applied force and  2 is the vertical displacement.It must be highlighted that these parameters depend on the position and load steps for physical nonlinearities, as presented in further sections.
In finite element approach, the approximated displacement in vertical direction,  � 2 may be obtained by interpolating the nodal responses,   , with the shape functions   , that is, By using the approximated variable, Equation 2, there are residues, , in the domain, , and in the contour, , such that (Equation 3), The weighted residual sentence is used to minimize them, Where   and  �  are, respectively, the weighting functions in  and .
In finite element method the Dirichlet boundary conditions are automatically satisfied prior to solving the system of equations, making null the second integral of Equation 4. In this sense, only the domain portion of the equation should be considered.Furthermore, by using the Petrov-Galerkin method, weighting functions, , are the shape functions,   , themselves.Thus, one can obtain (Equation 5) Taking the element domain as [−/2, +/2] and transferring the differential operator acting in the variable of interest, two times, to the shape function, through Green's divergence theorem, one obtains (Equation 6) or, in expanded way, in terms of nodal degrees of freedom,  ‾ 1 ,  ‾ 1 ,  ‾ 2 ,  ‾ 2 , which means the approximated vertical displacements and rotations at initial and final nodes (Equation 7) Rearranging the terms in Einstein notation, one has (Equation 8) Where   is the FEM stiffness matrix,   is the applied external forces vector and  ∂Ω  is the vector which contains the local components of bending moment and shear force that are canceled when assembling the global system.The stiffness matrix can also be defined through FEM  tensor, as further detailed ahead, by (Equation 9)

The laminated beam element
To consider the nonlinear behaviour of reinforced concrete by Damage Mechanics and Plasticity Theory in the Euler-Bernoulli FEM beam model, the cross-section is divided in layers, as indicated in Figure 1.During the progressive loading, the mechanical properties of each layer (elasticity modulus, damage parameter, etc.), and the effective area and moment of inertia of the whole section may be changed.Therefore, stresses and strains for each layer may be computed individually.

Mapping and equivalent stiffness concept
Despite the finite element of Euler-Bernoulli beam is a one-dimensional element with 2 degrees of freedom per node, it is possible to consider the other dimensions and physical parameters through mapping it to an auxiliary mesh, as represented by Figure 2.For this reason, it is necessary to determine an equivalent stiffness for each element.In the particular case of  2 symmetry in laminated beams of composite materials with constant width , the bending equivalent stiffness   can be determined by [17] (Equation 10) In which ( 2 )  is the  ℎ layer's bottom edge cartesian coordinate and ( 2 ) −1 is the ( − 1) ℎ layer's bottom edge cartesian coordinate.

Gaussian quadrature
The Gaussian quadrature is a procedure that numerically approximates the value of an integral exactly, enables practical use in domains whose integral functions are not standardized and minimizes computational processing time.Therefore, its use in association with FEM is very opportune.
Through this numerical integration method, it is only necessary to choose a quantity of points sufficiently to integrate a function exactly.As the Euler-Bernoulli beam element utilize Hermitian third order polynomial functions, 2 points are enough to approximate it with exact precision.These points have coordinates ±1/√3 and unitary weights.Those are the number of points and coordinates adopted in this study.
In addition to bringing a precision advantage to linear FEM analysis when calculating responses at Gaussian points instead of directly evaluating them at nodes, from the Damage Mechanics and Plasticity Theory perspectives it is also better to evaluate the nonlinearities effects at these points and then extrapolate the responses to the nodes, since numerical precision can affect the responses path due to the dependence of response history, especially in concrete damage.

The 𝑩𝑩 FEM tensor evaluated in Gaussian points
In bending beams theory, the rotation  is defined by (Equation 11) The strain in direction  1 of Euler-Bernoulli's beam is calculated by As can be seen, the strain defined by Equation 12 is defined in  1 coordinate.To calculate the strains in the Gaussian points, it is necessary to substitute the variable  1 to the variable  1 .Using the chain rule to transform the coordinate system, the   tensor is defined dependently of the Jacobian , as shown below (Equation 13), In this sense, the strain in each point of the cross-sections' layers is evaluated at Gaussian integration points of the elements,   , in the following way (Equation 14) Where  ‾  are the nodal degrees of freedom, or the approximated nodal responses of displacements and rotations.Since a variable substitution was performed, the stiffness matrix in Einstein notation, Equation 9 can be now calculated by (Equation 15) Where the Jacobian is defined by  = , and  is the length of the local finite element.

CONSTITUTIVE MODELS AND THE CONSIDERATION OF NONLINEARITIES
It is possible to introduce the effects of physical nonlinearities in a sub local region, or better, in a representative volume element, RVE.Thus, the third hypothesis presented of the classical Euler-Bernoulli beam element can now be disregarded.
Unlike a linear elastic regime, the stiffness of the problem now depends on the nodal degrees of freedom, displacements, and rotations, which in turn depend on the position vector of each layer.In this sense, we have (Equation 16) 3.1 Fundamentals of solid thermodynamics

Clausius-Duhem inequality
By combining the first and second laws of thermodynamics, so that a mechanical process if thermodynamically admissible, through Clausius-Duhem inequality [18], one has (Equation 17) Being (Equation 18) And (Equation 19) Where  is the mass density and  is the absolute temperature defined as a scalar field of positive values at each point in the domain  under consideration. is the specific entropy per mass unit,   is the stress tensor,  ˙ is the strain rate,  is the heat flow vector. ˙stored is the rate of change of the stored energy,  ˙dissipated is the rate of change of the dissipated energy,  ˙stored and  ˙dissipated are the rate of change in the density of internal energy, stored and dissipated, respectively, per mass unit.In a pure mechanical process, the Clausius-Duhem inequality is defined by

Thermodynamic formalism of Damage Mechanics
Part of  ˙stored remains stored in the system while another portion, recoverable, is a function of a reversible state variable, such as elastic strain.In this sense,  ˙stored is a function that depends on elastic strain and associated internal variables, that in rate terms is defined as [18]  ˙stored = Being  ˙  is the elastic strain rate for small strains consideration and  ˙ is the internal state variables rate of change.Combining the Equation 20and Equation 21, one can obtain [18] The Equation 22is extremely important for understanding of nonlinear physical processes, such as damage and plasticity.The second law of thermodynamics states that the rate of dissipated internal energy is always positive.In Equation 22, a portion of the recoverable stored energy rate and an unrecoverable one is observed.If the unrecoverable portion, third term, is null, remaining only the recoverable portion, there will be a fully reversible process.Therefore, the recoverable portion is always reversible.Purely elastic processes have no energy dissipation, being fully recoverable and reversible.In this sense, the energy dissipated can only arise from the unrecoverable portion of the energy rate stored in the system [18].
The unrecoverable portion is divided according to the material's state variables.In the most general case, the unrecoverable portion is divided into a reversible and an irreversible portion, which is the dissipated energy.The unrecoverable and not dissipated portion can be reversible, but only by an additional thermal process.In purely mechanical processes, this portion remains blocked in the volume.It can be unrecoverable, but it is not completely irreversible.Part of it is irreversible, dissipated, and another part, not dissipated, is reversible, but only by an additional thermal process [18].
In damaged elastic medium, as for the concrete, there are only recoverable energy and dissipated energy, that is, the entire portion of the unrecoverable energy is dissipated.In this sense, if the third term of Equation 22is null there will be no dissipated energy, that is, being a fully reversible process.If it is greater than zero, on the other hand, there will be dissipated energy.In elastoplastic responses with hardening, as for the reinforcement, otherwise, there is a portion of recoverable energy, reversible, an unrecoverable portion dissipated, irreversible, and an unrecoverable portion not dissipated, blocked in volume and reversible only by additional thermal process [18].

Damage Mechanics
Here, some concepts of Damage Mechanics will be briefly presented.

Elastic media with damage
In the case of elastic damaged materials, when unloading process occurs all strains can be recovered.Nevertheless, the Young's modulus is directly affected by the damage.In this way, the internal energy is determined by two state variables, the strain tensor   and the damage parameter , the latter being an internal variable.Thus,  = (  , ).Consequently [11], [19], In uniaxial models, strain energy can be defined as follows [19], By combining the Equation 23in Equation 24, one obtains [19] (Equation 25) And (Equation 26) Where (Equation 27) is the damage energy release rate associated with the damage variable.

The damage of each layer
In each layer of each cross-section of its respective Gaussian point, are attained a strain and a stress tensorial states,  and , dependents of the distance from the neutral axis, the stirrup distribution, and the moment of inertia variation.In these are calculated the concrete's damage and the steel's plasticity evolution processes.
The beam can suffer damage effect in a generalized way.The following hypotheses are considered: a) in a local scale, the damage effect occurs due to stretching; b) the damage is represented by a scalar variable in the sense that, when the variable reaches a certain value, a damage evolution happens; c) the damage is considered here as an isotropic variable; d) it is assumed that concrete is a nonlinear elastic medium with damage.
The damage evolution is based on the stretch that the material is subjected during the loading.For the equivalent strain of the Damage Mechanics, , it is used the definition proposed by Mazars and Pijaudier-Cabot [20], in order to adequate the density energy release rate for concrete's mechanical behaviour, defined as (Equation 28) With (Equation 29) Where 〈•〉 + are the Macauley brackets and   are the principal strains.
When  ˜ reaches certain value, the damage process begins.Taking  as the damage variable,  0 the maximum linear strain for the uniaxial tension test, and () a function that links the beginning damage strain with the damage variable, this process can be expressed by the damage loading function (Equation 30) So, the damage evolution can be measured by the rate of damage variable using the Karush-Kuhn-Tucker conditions, which is [21] (Equation 31), In which  ˙ is the damage rate, ⟨ ˜⟩+ is the positive portion of the equivalent strain rate, ( ˜) is a continuous and positive function of the equivalent strain  ˜ such that  ˙≥ 0 for any  ˜.This ensures that the damage is increasing and irreversible.In this sense, the function ( ˜) must be able to reproduce the experimental behaviour of one-dimensional, twodimensional, and three-dimensional experiments.According to Mazars [11], two rate damage variables can be deal with,  ˙ and  ˙, for the traction and compression, respectively (Equation 32), Where   () and   () are the () function for the tension and compression situations, respectively.Since the damage is a dissipative and accumulative process, the damage variable can raise in traction or in compression.This is an important aspect to be highlighted since all points in a beam can suffer tension or compressions due to loading and unloading or, for dynamic cases, due to dynamic responses and induced vibrations, as proposed by Abeche et al. [10].
The damage variable is split into two parts (Equation 33) Where   and   are the damage variables in tension and compression, respectively.They are combined with the weighting coefficients   and   , defined as functions of the principal values of the strains, due to positive and negative stresses [20] (Equation 34), And (Equation 35) In which 〈  〉 is the tensile strain vector, in other words, the strain associated with tensile stresses.In the same way, 〈  〉 is the compressive strain vector, associated with compressive stresses. is an adjustment relative to the response to the material's shear behaviour and, usually, is defined as  = 1.For the case of the material is subjected to uniaxial tension,   = 1 and   = 0. Similarly, if the material is under uniaxial compression,   = 0 and   = 1 [22].The damage variables   and   can be calculated through the Mazars' damage constituve model [11], as (Equation 36) And (Equation 37) Where   ,   ,   and   are characteristics parameters obtained from single tension and compression state tests.Mazars [11] initially suggested the following limits for these variables, to define the constitutive relation that associates the physical and mathematical formalisms of Damage Mechanics (Equation 38), Mazars later proposed changes in the limits of the calibration parameters [12].The stress-strain diagrams for the tension and compression cases considering the damage process are illustrated in Figure 3.

The plasticity in steel layers
The steel rebars of the reinforced concrete beam are reshaped into steel layers, maintaining its area and center of gravity properties.For the behaviour of the steel layers in Euler-Bernoulli beam finite element, a bilinear constitutive model with yielding and permanent strains is used, with symmetric behaviour in tension and compression, as illustrated in Figure 4.

Analytical variation of the neutral axis position
If the variation of the neutral axis position, ( 2 )  , is not contemplated, i.e, remaining unchanged in relation to the initial position, this prescribes that damages of the layers above the neutral axis, or above the center of gravity, would always occur by compression while the layers below the unchanged neutral axis would always suffer the phenomenon of traction.This can be valid if the neutral axis changes slightly.
On the other hand, if the position of the neutral axis varies significantly, the damage evolution itself can occur differently, and the consideration of a stagnant neutral axis can no longer be made.For example, a layer that would initially suffer compression phenomenon, in the case of ascension of the neutral axis, can now suffer the traction effect.The work of Abeche et al. [10] considers the variation of the neutral axis position iteratively.This paper considers the variation of the neutral axis' position analytically, inside each element, as follows (Equation 40) Where  is the area of the -th layer of the -th element.
This consideration allows a great reduction of the computational effort to find the equilibrium position, in comparison with iterative procedures.
After updating the position of the neutral axis ( 2 )   , the new coordinates are updated analytically according to the absolute coordinates of the layers and the new position of the neutral axis as follows (Equation 41) Where {( 2 ) rel }  is the vector that stores these updated  2 coordinates of each element inside each iteration.

Variation of the moments of inertia of the layers
In addition to considering the variation of the position of the neutral axis in an analytical way, this work adds as innovation the variation of the moments of inertia of each layer when the damage process occurs.
No studies were found that considered the variation in the moments of inertia of the layers when the damage process occurs.This physical consideration is plausible, since by varying the Young's modulus in the damage process the position of the neutral axis varies and, consequently, the moments of inertia of the layers, indeed, also vary.
The absolute global coordinates of the center of gravity of each layer are stored in the �  � vector.After updating the position of the neutral axis, according to Equation 40, the vector containing the moments of inertia of each layer of each element's Gaussian points is recalculated using Steiner's theorem, as follows (Equation 42) The Figure 5 illustrates the variation of the moment of inertia of the layers.

Shear stress and stirrups
Despite the classical Euler-Bernoulli beam theory neglects the shear deformation in its formulation, when the segmentation of finite elements in layers is applied, the shear strain and stresses can have a major influence on its Damage evolution.
To consider the shear stresses at the  ℎ layer of a Euler-Bernoulli beam, one has [23] (Equation 43) Where (( 2 )  ) is the  ℎ layer's first moment of inertia and  is the shear force acting on the element.When a layer is subjected to axial stresses due to bending and to shear stresses due to shearing forces, one must consider the principal stresses and strains to calculate the Damage according to Equations 36-37.
If the beam has stirrups, these must be accounted when calculating the shear stresses.When the concrete is still intact, that is,  = 0, the concrete is responsible for the shear resistance.The stirrup's contribution occurs when the first microcracks appear ( > 0), acting on the tensioned axis of the principal plane generated by the shear stress portion.This contribution, in terms of shear forces, is given by [24] (Equation 44) in which  is the effective depth of section,   the transversal reinforcement ratio, given by   () ⁄ , where   is the stirrup cross-section area and  is the spacing between stirrups.The stirrup axial stress,   , is obtained through (Equation 45) where   is the stirrup's Young Modulus,  the orientation of the principal plane 1, and  1 the principal strain associated with the axis 1, wherein To calculate the Damage, one must regard only the fraction supported by the concrete.Therefore, the effective shear force to be considered on the concrete,   , is given by (Equation 46)

NONLINEAR STATIC EQUILIBRIUM
To capture the effect of the damage evolution, the total load is divided into load steps, ∆, so that the model looks for the equilibrium damage state.
In each iteration, the displacements are obtained, the finite element matrix (  =      ) is calculated at the Gauss integration points of the beam elements.The strains of each layer are calculated at the same direction of Gaussian integration points, and, in these points, the damage status is checked.It is important to point that the neutral axis position is updated in each iteration.When damage occurs in any layer, the model updates the stiffness matrix, rechecks the balance setting, as the loss of stiffness due to the Damage Mechanics in any layer may result in some damage in another part of the structure.Then, the stress tensor is calculated considering the damage and, consequently, the global internal force vector is obtained with the contribution of degrees of freedom.

Deformed configuration, the second Piola-Kirchhoff stress tensor and its variation for varied materials
As shown, the nonlinear static equilibrium in this model is analysed in the undeformed configuration,  0 , that is, the total Lagrangian description is used.The stresses are calculated by the second Piola-Kirchhoff tensor,   .Although the Green-Lagrange strain tensor should be used in this case, in the context of small strains field, it is possible to use the infinitesimal strain tensor,   .It can be noted that the Cauchy stress tensor,   , could also be used through the updated Lagrangian description,   , as both descriptions are simply two ways for expressing the same problem [25].
The Damage Mechanics approach takes into consideration the concept of the effective stress tensor [18].The variation of the position of the neutral axis and the moments of inertia of the layers, due to damage, cause a change in the components of the stress tensor.Since the effective stress tensor depends on the damage variable, the stiffness loss increases the stress components.

Gain of moment of inertia of damaged layers due to the increase of distance from the neutral axis
When the damage process occurs in any layer there will be loss of resistance due to the reduction of Young's modulus.However, such degradation causes the position of the neutral axis to move away from the damaged layer, increasing their distance and resulting in a gain of the layer's moment of inertia.
One should note that this consideration could cause an equivocal gain of stiffness.Nonetheless, the gain of moment of inertia is never greater than the loss of Young's modulus, since the beam will never increase its original stiffness, as the damage is progressive.

Nonlinear static equation
The global equation considering physical nonlinearities and discretized through laminated Euler-Bernoulli beam finite element is (Equation 47) Although the verification occurs in the undeformed configuration, the damage is recalculated at each iteration within each loading step through Newton-Raphson iterative technique.
The convergence criterion here includes both terms of forces and displacements in relation to the tolerated equilibrium values, comparing the internal energy increment during each iteration, which is the amount of work done by the unbalanced load vector in displacement increments, with the initial internal energy increment [26] (Equation 48), Where   is the energy convergence tolerance that includes both the convergence tolerances in terms of force and displacement.

COMPUTATIONAL MODEL APPLICATIONS IN EXPERIMENTAL TESTS
A computational model was developed in C++ programming language with MATLAB wrapper to consider physical nonlinearities models by Damage Mechanics and Plasticity Theory in laminated Euler-Bernoulli beam finite elements.Secondly, to analyse the shear effects in reinforced concrete beams, Bresler and Scordelis [13] elaborated a set of 12 experimental beams, with different cross-section geometry, span and reinforcement ratio.These 12 beams are divided into 3 series of 4 beams with similar aspect ratios /ℎ.Due to the use of the Euler-Bernoulli beam element, one of these series will not be investigated in this work on the account of its short span.The regarded beams presented in this work are referred as BS-OA2, BS-A2, BS-B2, BS-C2, BS-OA3, BS-A3, BS-B3 and BS-C3.The relevance of these experimental tests, presented in Bresler and Scordelis [13], for the engineering structures field granted it to be investigated in several works [15], [16], [28].The test setup consists in three-point bending, as presented in Figure 6b-6c.Each cross-section geometry and reinforcement details are shown in Figure 7b-7i.The steel properties are described in the Table 1 and the concrete properties, obtained experimentally, are presented in the Table 2.The Damage parameters, presented in the Table 3, were calibrated to accord with the concrete resistances obtained by Bresler and Scordelis [13].

RESULTS
The maximum vertical displacement of the Mazars' beam is presented in the Figure 8.With the proposed model, the beam presented a flexure-tensile failure mode.As the load increased, the bottom concrete layers were submitted to damage due to tension, and neutral axis started to ascend.The lower steel rebars then reached the elastic threshold and the yielding amplified the damage process, leading the beam to large displacements.
The damage evolution of the Mazars' beam is presented in the Figure 9.It is noticeable, through the analysis of the neutral axis position and the damage configuration, that the damage occurred from bottom to top.The Figure 10 shows the maximum vertical displacement for the Bresler-Scordelis' beams.At the experimental tests, the BS-OA2 and BS-OA3 failed due to diagonal tension (D-T).The failure was reported as sudden and brittlelike [16].These beams are similar to the BS-A2 and BS-A3, respectively.The difference is the presence of stirrups in the latter, which granted a higher ultimate load strength.The failure mode of the BS-A2, BS-B2 and BS-C2 beams were reported as a combination between shear and compression stresses (V-C).
Lastly, the BS-A3, BS-B3 and BS-C3, due to the longer span, failed due to flexure-compression (F-C) [16].The top concrete was crushed because of the great amount of bottom steel reinforcement.The proposed model was able to reproduce this effect numerically in the simulation.On the numerical analysis with the proposed model, was noticeable that in these beams, the damage started at the bottom concrete layers.However, differently from the Mazars' beam, the amount of bottom rebars maintained the beam's stiffness and the neutral axis near its original position.The position of the neutral axis presents significant importance to the beginning of damage on the top concrete layers, once the distance between the neutral axis and a certain layer height is proportional to the strain of this layer.As the load increases, the damage also increases in both top and bottom layers, due to compression and tension, respectively.The final damage configurations of these beams are presented in the Figures 11-12.

CONCLUSIONS
The laminated Euler-Bernoulli beam theory was used here to analyse nonlinear reinforced concrete beams subjected to damage and plasticity due to progressive loading.The numerical results were compared with experimental values and numeric ones found in the literature.The nonlinear behaviour expressed in terms of load-displacement curves was accurately approximated by the proposed method, as can be seen in Figures 8 and 10.The damage process was well captured by the laminated beam model and the variation of neutral axis position can be highlighted in Figures 9,11   The present methodology was applied in experimental cases which the collapse was due to bending since the model adopted is based on the Euler-Bernoulli beam theory.
The relevance of considering the variation of the neutral axis in the damage and yielding procedures was highlighted, especially in comparison with bending tests on beams with a low rate of tensile reinforcement.
Furthermore, the previous consideration combined with the moment of inertia variation, the Damage Mechanics of concrete and Plasticity Theory of the steel and the shear stresses acting on the beam was able to numerically reproduce, with satisfactory approach, the results of the experimental tests both from the recent [12] work and the classic [13] research, when comparing force-displacement responses and damage configurations during the progressive loading.

Figure 1 .
Figure 1.Cross-section of a laminated Euler-Bernoulli beam element.

Figure 2 .
Figure 2. Consider of other dimensions, physical parameters, and local refinement in laminated Euler-Bernoulli beam element.

Figure 3 .
Figure 3. Concrete damage model for tension and compression stresses.

Figure 5 .
Figure 5. Layer properties and strain in each layer.

Figure 10 .
Figure 10.Maximum displacement in the BS beams.