## On-line version ISSN 1983-4195

### Rev. IBRACON Estrut. Mater. vol.10 no.2 São Paulo March/Apr. 2017

#### http://dx.doi.org/10.1590/s1983-41952017000200008

Articles

Finite Element Model for Nonlinear Analysis of Reinforced Concrete Beams and Plane Frames

bDepartamento de Egenharia Civil, Universidade Federal de Santa Catarina, Florianópolis, SC, Brasil. henriettelarovere@gmail.com

Abstract

In this work, a two-dimensional finite element (FE) model for physical and geometric nonlinear analysis of reinforced concrete beams and plane frames, developed by the authors, is presented. The FE model is based on the Euler-Bernoulli Beam Theory, in which shear deformations are neglected. The bar elements have three nodes with a total of seven degrees of freedom. Three Gauss-points are utilized for the element integration, with the element section discretized into layers at each Gauss point (Fiber Model). It is assumed that concrete and reinforcing bars are perfectly bonded, and each section layer is assumed to be under a uniaxial stress-state. Nonlinear constitutive laws are utilized for both concrete and reinforcing steel layers, and a refined tension-stiffening model, developed by the authors, is included. The Total Lagrangean Formulation is adopted for geometric nonlinear consideration and several methods can be utilized to achieve equilibrium convergence of the nonlinear equations. The developed model is implemented into a computer program named ANEST/CA, which is validated by comparison with some tests on RC beams and plane frames, showing an excellent correlation between numerical and experimental results.

Keywords: nonlinear analysis, finite element; reinforced concrete; beams; plane frames

1. Introduction

Due to the advance in technology combined with the use of more resistant materials, more complex and slender structures are currently being designed, arising thus the necessity of more elaborated computational methods for structural analysis and design.

In structural analysis of reinforced concrete (RC) structures, the more refined methods should account for the structure nonlinear behavior, due to material nonlinearities (physical nonlinearity), as well as to changes in the deformed shape of the structure (geometric nonlinearity). Among the more refined methods, the Finite Element Method (FEM) stands out as one of the most utilized nowadays, which, for the case of plane structures, can be employed by using either bar or plane elements. Yet only a few computer programs are available for nonlinear finite element analysis of RC structures, and their cost is high in comparison with other programs. Due to this fact, refined finite element (FE) models have been mostly used by researchers, and many researchers opt in developing their own models and computer programs. Although several nonlinear FE models have already been developed, this is still an advanced topic in the scientific-technical community, in view of the difficulty on accurately model the reinforced concrete material, due to cracking of concrete, yielding of steel and the interaction between these materials. Hence, the development of models that combine computational efficiency and good accuracy needs to be more and more incentivized (Silva and Matos [1]). In this work, emphasis is given to the bar element model, as it yields a reduced number of degrees of freedom as compared to plane element models, making viable the analysis of large structures, which is the aim of the research project under way at the Federal University of Santa Catarina (UFSC), with participation of the authors.

The early FE bar models were developed in the 60's, and were limited to the analysis of small members and small structures, using simplified constitutive models. One of the first bar models was the one by Giberson [2], which consisted of a linear elastic element connected by nonlinear springs at its ends, in which predefined moment-rotation relations were utilized. Since then the bar models evolved, particularly with the introduction of the Fiber Model (Kaba and Mahin [3]), which subdivides the element section into overlaid concrete and reinforcement layers, by considering nonlinear constitutive laws for the materials in each layer. Another important contribution to the evolution of this kind of model was the introduction of an internal node at the element midpoint with only one axial degree of freedom. As demonstrated by Chan [4], the inclusion of this third node allows a proper representation of the element flexural stiffness with variation of the neutral axis position, caused by material nonlinearities. Not including this third node imposes a constraint on the element, making it artificially stiffer. Holzer et al [5] have used this kind of model for physical and geometric nonlinear analysis of RC beams/columns. Marí [6] has extended this model to a three-dimensional bar element, with a total of thirteen degrees of freedom (six in each external node and one axial degree of freedom in the internal node), by discretizing the element section into filaments, and by taking into account long term effects caused by creep and shrinkage. In the model two Gauss points are utilized to integrate the stiffness matrix and the internal forces vector of the element, with the constitutive matrix being evaluated only at the element midpoint, and by neglecting the tensile contribution of the intact concrete between cracks (tension-stiffening). The geometric nonlinearity is also considered, by means of the Updated Lagrangian formulation. More recently Marí [7] has improved this model, by including a tension-stiffening model developed by Carreira and Chu [8]. In Brazil, Schulz and Reis [9] have utilized a model similar to the one by Marí [6] to analyze three-dimensional RC frame structures, by considering physical nonlinearity by means of constitutive equations recommended by Design Codes (NBR-6118 and CEB 90), disregarding the tension-stiffening effect, and by considering the Total Lagrangian formulation for geometric nonlinearity. There is also another kind of FE bar model that uses a formulation in terms of forces instead of displacements, as the model by Taucer, Spacone and Filippou [10]. In Brazil this kind of model was utilized by Teixeira and de Souza [11] in the three-dimensional analysis of a reinforced concrete building, by using a computer program named OpenSees from the University of California in Berkeley. This model was validated by comparison with a model that uses a co-rotational formulation and with the method known as P-delta. The concept of co-rotational formulation, which allows that nodes undergo large displacements and rotations, as well as bars display large elongations and curvatures, has been presented by Pimenta [12]. Pimenta and Soler [13] applied this formulation to analyze one RC beam and two RC plane frames, by adopting a compressive constitutive law for concrete similar to the one given in NBR-6118, by neglecting the concrete tensile strength, and by considering the steel as an elastic-perfectly plastic material. A similar model, with a co-rotational coordinate system attached to the element, was utilized by Silva and Matos [1]. The authors utilized the Fiber Model and considered the contribution of concrete between cracks by means of the tension-stiffening model developed by Vecchio and Collins [14]. Pinto [15] and Carvalho [16] have also utilized a co-rotational formulation for physical and geometric nonlinear analysis of RC structures. Further details on the literature review of FE models for the analysis of RC structures can be consulted in Stramandinoli [17].

From this review one can conclude that the models based on the Force Method have presented excellent results, however their computer implementation becomes more difficult, especially in the usual FE programs that utilize a formulation in terms of displacements instead of forces. Regarding the geometric nonlinearity, despite of the co-rotational formulation being more complete, the large displacement with moderate rotations assumption is in general sufficient to represent the behavior of usual RC structures, such as beams and frames, since large displacements and rotations would be incompatible with the structure utilization. Hence, in this work, a two-dimensional bar element model based on the model developed by Marí [6] is presented. The assumption of large displacement with moderate rotations is considered and the Total Lagrangian formulation is utilized, which is easier to be computationally implemented as compared to the co-rotational formulation. The advantage of this model with respect to similar models described above, lies on the inclusion of a novel tension-stiffening model proposed by the authors (Stramandinoli and La Rovere [18]), which presents a tensile constitutive law for the concrete between cracks as a function of the reinforcement ratio of the bar element, whereas in the previously proposed tension-stiffening models - Vecchio and Collins [14], Carreira and Chu [8], etc., the same constitutive law is used independently of the reinforcement ratio. With that the proposed model can represent more realistically the nonlinear behavior of RC structures after cracking. Another advantage is that in the proposed model the constitutive matrix is evaluated along the element at the three integration points (Gauss points), which allows the use of coarser FE meshes to capture the spread of nonlinearities along the structure. At each Gauss point the element section is discretized into layers, and it is assumed that each layer is under a uniaxial stress-state. The model does not take creep and shrinkage effects into account and it is limited to RC structures with a dominant flexural behavior, where shear deformation is neglected. It should be pointed out that, in those structures where shear effects become important, bar models based on Timoshenko Beam Theory should be used, as for instance the one developed by the authors in Stramandinoli and La Rovere [19], and Stramandinoli [17], or, alternatively, plane finite element models, as, for example, the ones developed by d´Avilla [20], should be employed.

The FE model developed by the authors is presented in Section 2, in the following, and the constitutive equations of the materials utilized are described in Section 3. The model is implemented into a computer program named ANEST/CA, which is validated in comparison with an analytical model developed by another author for the case of geometric nonlinearity for large displacements, and in comparison with experimental tests on beams, by considering physical nonlinearity only, and on plane frames, by considering both nonlinearities, as presented in Section 4. At the end of the work, in Section 5, a few conclusions are extracted.

2. Finite Element Model

A nonlinear model based on the Finite Element Method with isoparametric formulation where the structure is discretized into bar finite elements, is developed. The Euler-Bernoulli assumption, in which shear deformation is disregarded, is adopted.

The bar finite element utilized has three nodes and a total of seven degrees of freedom (Figure [1]). The two external nodes have three degrees of freedom each: two translations (axial and transversal) and one rotation. The internal node at the element midpoint has only one axial degree of freedom, similar to the one utilized by Chan [16] and by Mari [2]. Upon inclusion of this node, the horizontal displacement field in the element becomes compatible (see equation 1 ahead, the first term shows a parabolic variation with x, the horizontal axis, as well as the second term due to bending, since v(x) varies cubically with x). This allows the axis x to have an arbitrary position, fixed during the analysis but not necessarily coincident with the line passing by the section centroids, in such a way that the element stiffness can be properly represented upon variation of the neutral axis position, caused by cracking and other material nonlinearities (further details can be found in Stramandinoli [17]).

The numerical integration of the element is performed by means of three Gauss points, and at each point the section is discretized into concrete layers overlaid to longitudinal reinforcement layers (Fiber Model). It is assumed that concrete and reinforcing steel are perfectly bonded, and that each layer is under a uniaxial stress-state. Nonlinear constitutive laws are utilized for the materials, as described in Section 3. Regarding the geometric nonlinearity, the Total Lagrangian Formulation, considering moderate rotations, is utilized.

The element formulation is described in the following, by considering initially linear-elastic material, and in the sequence including physical and geometric nonlinearities.

Formulation for linear-elastic material

By considering the Euler-Bernoulli Beam Theory, the displacement field along the element is given by:

(1)

where u is the longitudinal displacement, 𝑢 0 is the longitudinal displacement along the reference axis x, 𝑣 is the transversal displacement and 𝜃 is the rotation of the cross-section.

The longitudinal strain and stress in the element are:

(2)

(3)

Where 𝐸 is the longitudinal modulus of elasticity of the linear-elastic material.

By introducing the natural coordinate 𝜉= 𝑥 𝐿 , the displacement field in terms of nodal displacements can be written as:

(4)

(5)

(6)

(7)

in which: `= 𝑑 𝑑𝜉

and where 𝛼 1 is related to the displacement of the internal node, 𝑢 3 :

(8)

and the interpolation functions 𝑁(𝜉) are given in the Annex.

From the displacement field, the strain in the element can be determined (see details in the Annex). The strain and stress vectors are given by:

(9)

where 𝑈 ~ is the nodal displacement vector; 𝐵 ~ 𝐿 is the linear matrix that relates strain to nodal displacements, which is obtained from equations (2) and (4) to (7):

(10)

and 𝐷 ~ =𝐸 is the constitutive matrix for linear-elastic material, in the uniaxial case.

By applying the Principle of Virtual Work and after some algebraic manipulations (Stramandinoli [17], Cook [21]), the equation that defines the stiffness matrix can be found:

(11)

and also the internal forces vector in the element:

(12)

The integrals in equations (11) and (12) can be obtained by means of numerical integration, and, in this work, the Gauss integration rules with three integration points are utilized.

Next, a static-condensation procedure is applied to condense out the seventh degree-of-freedom in the element stiffness matrix and in the element force vector, from which the global stiffness matrix and the global internal forces vector of the structure are assembled.

2.2 Formulation including physical nonlinearity only

In order to include physical nonlinearity, besides the Euler-Bernoulli assumption, it is assumed that: the element undergoes small strain and displacements; concrete and steel are homogeneous materials and there is a perfect bond between them; the element cross-section is discretized into layers, by considering that each layer is under a uniaxial stress-state; the tension-stiffening effect occurs in the concrete after cracking; the reinforcing steel is an elastic-plastic material with strain-hardening; the total internal forces at each section are obtained upon superposition of the forces derived from the stresses in the concrete layers to those arising from the stresses in the reinforcement layers.

Upon inclusion of physical nonlinearity, equations (11) and (12) need to be modified, since the stresses and the constitutive matrix vary along the element, 𝜎 ~ = 𝜎 ~ 𝜀 and 𝐷 ~ = 𝐷 ~ (𝜀), therefore an iterative procedure becomes necessary to achieve force equilibrium at each load step. Solution of the nonlinear equilibrium equations can be obtained by means of either the Newton-Raphson methods or the Arc-length method (Stramandinoli [17]). Loads are applied incrementally, and for each load step the internal forces vector and the tangent stiffness matrix of the element are calculated at each iteration of the iterative procedure:

(13)

It can be demonstrated that the tangent stiffness matrix 𝑘 ~ 𝑡 is given by (Stramandinoli [2]):

(14)

where 𝐷 ~ 𝑡 and 𝐷 ~ are the tangent and the secant constitutive matrix of the material, respectively.

For comparison with the well-known stiffness matrix (6 × 6) of a linear-elastic plane frame element, the secant stiffness matrix (7 × 7) of the element considering only physical nonlinearity, which contains coupling terms between axial and bending stiffness, is presented in the Annex. In the initial elastic stiffness matrix, these terms become null for the case that the longitudinal axis x of the element passes by the centroid of the section (the first moment of area, S, is zero)

2.3 Formulation including both physical and geometric nonlinearities

Regarding the geometric nonlinearity, the Total Lagrangian Formulation, with the simplification for moderate rotations, is utilized. By considering moderate rotations, equations (1) remain valid, however the strain becomes:

(15)

This strain can be separated into two parts, a linear one, 𝜀 𝐿 , equivalent to the one given by equation (2), and another one nonlinear, 𝜀 𝑁𝐿 . Rewriting the strain in matrix notation, it follows:

(16)

where:

which is equation (2) presented previously and

(17)

By rewriting equation (6) using matrix notation, it follows:

(18)

and its derivative with relation to x becomes:

where

equation (17) can be written in the following way:

(19)

The incremental form of strain is given by:

(20)

By replacing the expressions that define 𝜀 𝐿 , equation (9), and 𝜀 𝑁𝐿 , equation (19), into the equation (20) above, yields:

(21)

or else:

(22)

and where 𝐵 ~ 𝐿 is the matrix that relates strain with displacements (linear), equivalent to equation (10).

For this case, with both nonlinearities considered in the model, the internal forces vector becomes:

(23)

The tangent stiffness matrix is still given by the expression:

(24)

But in this case, upon derivation of the terms in equation (24), it leads to:

(25)

Hence, it follows that the tangent stiffness matrix 𝑘 ~ 𝑡 is composed of three matrices:

- One matrix 𝑘 ~ 0 obtained considering physical nonlinearity, expressed as:

(26)

- Another matrix, usually defined as geometric matrix, 𝑘 ~ 𝑔 , given by:

(27)

- And one matrix 𝑘 ~ 𝑢 caused by initial displacements by considering the structured deformed, which is given by the addition of three components:

(28)

If instead of the Total Lagrangian Formulation, the Updated Lagrangian Formulation were utilized, this last matrix would not enter in the formulation, however it would become necessary to update the coordinates in all steps.

For both this case, in which both nonlinearities are considered, and the case of physical nonlinearity only (section 2.2), the stiffness matrix and the internal forces vector are obtained by means of numerical integration, using the Gauss Rules with three integration points along the longitudinal axis. At each Gauss point the constitutive matrix and the stress vector are obtained upon addition of values arising from the concrete and from the reinforcement layers, which in turn are obtained from the nonlinear constitutive equations for uniaxial state, as described in the following Section 3.

After generation of the element stiffness matrix and internal forces vector, a static-condensation procedure is applied to condense out the seventh degree-of-freedom of the element, and the global stiffness matrix and the global internal forces vector of the structure are assembled. Solution of the nonlinear equilibrium equations can be obtained by means of either the Newton-Raphson methods (Cook [21]) - either the tangent or the initial stiffness method, or else the Arc-length method (Riks[22] and Wempner[23]), when the equilibrium path presents a limit point or exhibits snap-through phenomenon. This nonlinear FE model using bar elements was implemented into a computer program called ANEST/CA, previously named ANALEST (Stramandinoli [17]), developed in FORTRAN 90 language.

3. Uniaxial constitutive equations

3.1 Concrete under compression

For concrete subjected to compression, both the modified Hognestad model and the CEB model can be utilized in the ANEST/CA program, but in this work only the modified Hognestad model, described in the following, was utilized.

The Hognestad [24] model has already been utilized by several authors, showing good results in comparison with experimental tests. In this work this model is modified, using a parabola to describe the compressive stress-strain curve for both the ascending and the descending branch, after the peak:

(29)

where

/is the compressive strength of concrete and 𝜀 0 is the corresponding strain.

3.2 Concrete under tension

For concrete under tension, the model proposed by the authors in Stramandinoli and La Rovere [18] is utilized in this work for all examples. In this model, the material is assumed to behave linear-elastically until the tensile strength of concrete is reached, and, beyond cracking, the tension stiffening effect, described by the following equation, is considered:

(30)

(31)

where:

/ is the tensile strength of concrete and 𝜀 𝑐𝑟 is the corresponding strain; 𝜀 𝑦 is the strain value corresponding to yielding of the reinforcement, after which the stress 𝜎 𝑐𝑡 drops abruptly to zero; 𝛼 is an exponential decay parameter, which is a function of the reinforcement ratio (𝜌) and of the steel-to-concrete modular ratio (𝑛), defined by the equation:

(32)

This equation was derived for bars, where the entire member is subjected to tension, and it was validated by comparison with experimental results from reinforced concrete bars subjected to direct tension, using different reinforcement ratios, showing an excellent correlation (Stramandinoli and La Rovere [18]). In order to apply this model to beams, the effective area that corresponds to the tensile zone in the member section needs to be obtained, which can be estimated using the equation given in the CEB-FIP Model Code 1990 [25]:

(33)

where

h is the nominal depth of the beam, 𝑑 is the effective depth, and 𝑘 𝑥 is the neutral axis depth.

Recalling that in RC beams under bending, the relationship between nominal and effective depth is usually given by ℎ−𝑑≅0,1ℎ, the effective area can then be expressed, approximately, as:

(34)

This approximate equation will be utilized in this work to calculate the reinforcement ratio in all examples of beams and plane frames, by considering the tension-stiffening effect only in this effective area. In the numerical analyses where the tension-stiffening effect is not considered, the stress drops abruptly to zero in the stress-strain curve, after the tensile strength limit, 𝑓 𝑐𝑡 , is reached.

It should be noted that the tension-stiffening effect is more accentuated for smaller values of the 𝛼 parameter, which means smaller values of reinforcement ratios, since the 𝛼 parameter increases with increasing values of 𝜌 and/or 𝑛. In order to illustrate this effect and show the influence of the reinforcement ratio on the tension-stiffening model proposed by the authors, a graph of the proposed model upon variation of 𝑛𝜌 is shown in Figure 2b, in comparison with other tension-stiffening models - the one by Vecchio and Collins [14] and the bilinear model from Figueiras [26], both displayed in Figure 2a, which are independent of the reinforcement ratio.

3.3 Reinforcing steel

It is assumed that the reinforcing steel, under tension and compression, is an elastic-plastic material, modeled by a bilinear stress-strain curve. In order to avoid convergence problems and oscillations in the iterative process, a parabolic curve is fitted between the elastic and plastic branches of the bilinear stress-strain curve, between 0.8 and 1.2 𝜀 𝑦 (La Rovere [27]). Strain hardening of the reinforcing steel may or may not be considered, through the use of a coefficient 𝑠ℎ, which is the plastic-to-elastic modular ratio (𝑠ℎ = 0 for horizontal threshold, perfectly plastic steel). The ultimate strain is called 𝜀 𝑢 and the corresponding stress, 𝑓 𝑢 .

4. Comparison of results given by the proposed model with analytical and experimental results obtained by other authors

In order to validate the proposed nonlinear model, its numerical results were compared to several analytical and experimental results from examples of reinforced concrete plane structures, available in the literature, by Stramandinoli [17]. Among these, a few examples of beams and plane frames were selected to be presented here, giving emphasis to structures with dominant flexural behavior. Only one example of plane frame subjected to large displacements was chosen, by considering only geometric nonlinearity, and it will be compared to another analytical model in section 4.3. The other examples - RC beams considering only physical nonlinearity, and RC plane frames considering both nonlinearities, physical and geometric, will be compared to experimental tests.

Regarding the choice of the FE mesh, parametric studies conducted by Stramandinoli [17] in simply-supported beams have shown that, under 4-point bending, in which the central span is subjected to pure bending, the solution obtained using a mesh of 4 elements basically coincided with the one obtained using 24 elements, showing that the model is objective, without mesh dependency. As for the case of beams under a concentrated load applied at mid-span, convergence of the solution was obtained for meshes with at least 10 elements. On the other hand, for a plane frame of one span and one story, fully fixed at the base, a certain mesh dependency was observed, though the load-displacement curves were basically coincident until the peak, the values of ultimate load and corresponding displacement varied a little by refining the mesh up to 20 elements. The author thus recommends that a finer mesh be utilized for plane frames, but imposing a restraint of not employing elements of length less than its cross-section depth, as recommended by Bazant et al. [28]. Stramandinoli [17] has also investigated, for several examples, the effect of number of layers used to discretize the element section, concluding that for 10 layers or higher there was no change in the numerical solution, hence in this work 20 concrete layers were adopted in all examples.

In all numerical analyses, the modified Hognestad model was utilized for concrete under compression and the Newton-Raphson method was applied for solution of the nonlinear equations, except for the example in section 4.3, where the Arc-length method was employed.

In the examples of beams and frames tested experimentally, the material properties for concrete and steel measured experimentally were utilized in the numerical model, whenever available (shown in highlight in Tables 1 and 2), otherwise they were estimated using the values and equations recommended by the Brazilian Code for concrete structures, NBR-6118.

Table 1 Material properties used in numerical analysis of beams VT1/VT2/VB6/VC3 (experimental values in bold, calculated or estimated values without bold).

 Beam Concrete Tension - stiffening Reinforcing Steel f cm (MPa) f tm (MPa) ε o n ρ eff nρ eff α (5 layers) φ f y (MPa) E s (GPa) sh VT1 / VT2 33.58 2.62 0.0020 6.39 1.50 % 0.096 0.040 6 mm 738 214.83 0.016 10 mm 565 214.83 0.000 VB6 37.9 2.90 0.0020 5.15 4.60 % 0.240 0.072 3 mm 192 174 0.001 8 mm 497 195 0.0042 VC3 20.7 1.60 0.0020 9.00 3.80 % 0.342 0.093 12.5 mm 507 184.6 0.0014

Table 2 Material properties used in numerical analysis of frames (experimental values in bold, calculated or estimated values without bold).

 Frame Concrete Reinforcing Steel f cm (kN/m 2 ) f tm (kN/m 2 ) nρ eff α (5 layers) ε o f y (kN/m 2 ) E s (kN/m 2 ) ε u sh P2 36500 2814 0.464 0.113 0.0023 293000 200000000 0.01 0.01 A40 29096 2303 0.443 0.110 0.002 353000 189791000 0.015 0.029 A60 38955 2974 0.317 0.088 0.002 425406 181797000 0.007 0.063

4.1 Simply-supported beams tested by Beber

Among the simply-supported beams tested under 4-point bending by Beber [29], two of them, VT1 and VT2, are initially utilized for comparison with the results generated by ANEST/CA program, considering only physical nonlinearity. The two beams are identical, and their geometry and reinforcement, the applied loading, and the mesh utilized in the finite element analysis are all displayed in Figure 3, with dimensions given in cm. The material properties are presented in Table 3. In order to demonstrate the importance of the tension-stiffening effect, a numerical analysis without considering this effect was also performed.

Comparison between numerical and experimental results is presented in Figure 4 in terms of a graph "total applied load versus vertical displacement at mid-span". It can be observed that the FE model could capture very well the ascending branch of the experimental curve obtained in the experimental tests, showing a good approximation even after the onset of cracking when the tension-stiffening effect is considered. When such effect is not considered, the numerical model displays a much more flexible behavior than the experimental one, as expected, since this effect is more accentuated in beams with low reinforcement ratio. The onset of yielding of reinforcement was accurately captured by the numerical model, corresponding to a total applied load of 44 kN. However, beyond that, the post-yielding response of the specimens could not be measured experimentally since the instruments have been removed to avoid damage. The ultimate total load measured experimentally was 47 kN, while the numerical value obtained at failure by program ANEST/CA was a bit lower, 46 kN.

4.2 Simply-supported beams tested by Juvandes

Among the beams tested in the experimental program conducted by Juvandes [30], two simply-supported beams under 4-point bending (VB6 and VC3) were selected to be analyzed with ANEST/CA program, by considering only the physical nonlinearity. Beam geometry and reinforcement, loading and FE mesh used in the analyses are shown in Figure 5 (dimensions given in cm), while the material properties are given in Table 1. The numerical analyses were performed with and without considering the tension-stiffening effect, to illustrate the importance of such effect.

Comparison between numerical and experimental results is presented in terms of "total applied load versus vertical displacement at mid-span" graphs in Figures 6 and 7. For beam VB6, it can be observed from Figure 6 an excellent approximation of the numerical model, considering tension-stiffening, with respect to the experimental test, with the curve obtained numerically basically coinciding with the one obtained experimentally, in the post-cracking and post-yielding ranges. The numerical model was only a little stiffer at the beginning of the elastic range, before the onset of cracking, region more susceptible to instrumentation imprecision, and it predicted very well the ultimate load, with a slightly smaller corresponding displacement.

A more flexible behavior of the numerical model is again observed when the tension-stiffening effect is not considered in the analysis (NO T.S.). For the other beam, VC3, the numerical model considering tension-stiffening reproduced very well the experimental model in the elastic range but, beyond that, for a total load higher than 25 kN, it became slightly stiffer than the experimental model (see Figure 7). The numerical model also predicted well the ultimate load, but showed a corresponding displacement smaller as compared to the experiment. For this particular beam VC3, the differences in the analyses with and without tension-stiffening consideration are smaller because, besides its higher reinforcement ratio, the elastic modulus of concrete is low, resulting in a high value for the modular ratio 𝑛, which consequently results in a high value for the 𝛼 parameter (see Table 1); therefore there is less tension-stiffening effect as compared to the beams analyzed before.

4.3 Plane frame studied by Williams

In order to verify in this work the geometric nonlinearity formulation of the proposed model, and in view of the lack of examples of experimental testing on RC structures subjected to large displacements (since such tests are of difficult execution), an example of a plane frame composed of two bars slightly inclined and made of hypothetical material is analyzed. This frame was studied by Williams apud Peterson and Petersson [31], and is referred in several works related to this topic. The frame geometry is displayed in Figure 8 and the element properties are: modulus of elasticity E = 70.735 GPa; cross-sectional area A = 1.18 cm²; and moment of inertia I = 0.0374 cm4.

In the numerical analysis the structure was discretized into 20 elements of equal length and the material was assumed to be linear-elastic. Firstly the model was applied using the matrices 𝑘 𝑔 and 𝑘 𝑢 , and secondly just using the 𝑘 𝑔 matrix. For solving the nonlinear equations, the Arc-length method was utilized, since the structure presents a critical point. Comparison of both numerical analyses with the one obtained from the analytical model used by Petersson and Petersson [31], expressed in terms of a graph - "load versus vertical displacement at the center of the frame", is shown in Figure 9.

The model with both matrices 𝑘 𝑔 and 𝑘 𝑢 could capture the entire response of the structure, including the post-critical range, capturing also the effect known as snap-through, with the curve obtained by the numerical model coinciding with the one from the analytical model of Petersson and Petersson [31]. However, when the matrix 𝑘 𝑢 was not included in the formulation, converge problems arose in the analysis close to the response peak, and it was not possible to capture the response for the complete load history. It should be pointed out, however, that, in many examples of RC plane frames, the model with just the geometric matrix, 𝑘 ?? , is sufficient to obtain good results (Stramandinoli [17]).

4.4 Plane frame tested by Cranston

Cranston, apud Bazant et al. [28], has conducted several experimental testing on RC frames hinged at the base, from which one, P2, was selected here for comparison with ANEST/CA program. Several researches have analyzed this frame; besides Bazant et al. [17], it can be quoted: Lazaro and Richards [32], Sun et al. [33], and Bratina et al. [34].

Figure 10 shows the frame geometry, the cross-section of the members, the load application points, the supports, and the mesh utilized to discretize the structure (18 elements) in the numerical analysis. The material properties used in the numerical analysis are given in Table 2.

Figure 11 illustrates the graph of "total load (kN) versus vertical displacement at mid-span (mm)" for frame P2. From the graph one can realize that the curve obtained using the numerical model basically coincided with the one obtained experimentally until the load peak was reached; the numerical model just showed to be a little stiffer in the post-yielding range of the reinforcement. After reaching the peak resistance of the frame, the experimental curve showed a descending branch, effect known as softening. In the numerical analysis (considering tension stiffening), the onset of yielding of the reinforcing steel occurs at mid-span of the horizontal bar for a load P = 15.7 kN, whereas, for P=20.4 kN, yielding of the rebars initiates at the ends of the horizontal bar and also at the top of the columns. Close to the ultimate load (between 22.2 and 22.3 kN), started the convergence problems, and it was not possible to continue the analysis and capture the post-peak branch, by either using the Arc-length method or the Newton-Raphson method under displacement control. Such convergence problems are commonly found in analyses considering physical nonlinearity when one or more coefficients in the diagonal of the global stiffness matrix of the structure are close to zero.

4.5 Plane frames tested by Ernst et al.

Ernst et al. [35] perform a study on the behavior of RC plane frames, testing several frames hinged at the base, of one span and one story. Among the tested frames, two of them (A40 and A60) were selected for comparison with the numerical model. The frame geometry, the cross-section of the members, and the load application points are displayed in Figure 12. A total of 36 bar elements were used in the structure discretization for the numerical analysis, and the material properties and parameters utilized are shown in Table 2.

Figure 13 illustrates the graph of "total load (kN) versus vertical displacement at mid-span (mm)" for frame A40 and Figure 14 for frame A60. It can be observed from both graphs that the numerical model could capture satisfactorily the behavior of the frames observed in the tests. The numerical model showed to be a little stiffer since the initial elastic range, and this difference increased a little in the post-cracking range. For both frames, the ultimate load obtained numerically was larger than the one obtained experimentally, and for frame A40 the ultimate displacement was somewhat smaller in comparison with the experimental value.

5. Conclusions

A 2D finite element model, using bar elements with seven degrees of freedom for physical and geometric nonlinear analysis of reinforced concrete structures, was presented in this work. The model was implemented into a computer program named ANEST/CA and it was verified in comparison with analytical and experimental results obtained by other authors.

Regarding the geometric nonlinearity, the model could capture very well the behavior of linear-elastic structures under large displacements, by using the formulation with matrices 𝑘 𝑔 and 𝑘 𝑢 . As far as physical nonlinearity is concerned, comparison between the numerical model and experimental tests for structures with dominant flexural behavior, in terms of load-displacement curves, showed very good agreement for the case of beams and good agreement for the case of plane frames. This work also demonstrated the importance of considering the tension-stiffening effect, especially in beams. Besides the numerical results presented here, the ANEST/CA program can also provide other important information in structural analysis, such as the evolution of cracking and of stress in the concrete and in the reinforcement layers, along the load history (Stramandinoli [17]).

When shear effects become important, in the presence of inclined cracks, another model that takes that into consideration should be employed, as, for instance, the bar model developed by Stramandinoli and La Rovere [19], based on Timoshenko Beam Theory, or, alternatively, plane finite element models using bi-axial constitutive models, as, for example, the ones developed by d´Avilla [20].

The ANEST/CA program was also used for comparison with simplified methods by Junges [36] for the case of beams, by Gelatti [37] for plane frames, and by Junges and La Rovere [38] for the case of continuous beams. The numerical model described here is currently being extended and implemented into ANEST/CA program to allow for the analysis of three-dimensional reinforced concrete structures, by including also the effect of confinement in concrete provided by the stirrups.

6. References

[1] SILVA, R. M.; MATOS, E. F. Análise não-linear de pórticos planos de edifícios altos em concreto armado considerando a contribuição do concreto tracionado. In: Jornadas Sudamericanas de Ingeniería Estructural, Punta del Este, p. 1-15 Anais CD-ROM, 2000. [ Links ]

[2] GIBERSON, M. F. The response of nonlinear multi-story structures subjected to earthquake excitation. PhD Thesis. 1967. 232 p. California Institute of technology, Pasadena, California, USA, 1967. [ Links ]

[3] KABA, S.; MAHIN, S. A. Refined modeling of reinforced concrete columns for seismic analysis. EERC Report 84-03. Earthquake Engineering Research Center University of California, Berkeley, USA, 1984. [ Links ]

[4] CHAN E. C. Nonlinear geometric, material and time dependent analysis of reinforced concrete shells with edge beams. 1982. 361 p. Ph.D. Dissertation (Structural Engineering and Structural Mechanics). University of California, Berkeley, USA, 1982. [ Links ]

[5] HOLZER, S.M.; SOMERS, A.E.; BRADSHAW, J.C. Finite response of inelastic RC structures. Journal of the Structural Division (ASCE), v. 105, n. ST1, p. 17-33, 1979. [ Links ]

[6] MARÍ A.R. Nonlinear geometric, material and time dependent analysis of three dimensional reinforced and prestressed concrete frames. Report, no. UCB/SESM - 84/12. Department of Civil Engineering University of California, Berkeley, USA, 1984. [ Links ]

[7] MARÍ A.R. Numerical simulations of the segmental construction of three dimensional concrete frames. Engineering Structures, v. 22, p. 585-596, Ed. Elsevier, 2000. [ Links ]

[8] CARREIRA, D. J.; CHU, K.H. Stress-strain relationship for reinforced concrete in tension. ACI Journal, v. 83, n.3, p. 21-28, 1986. [ Links ]

[9] SCHULZ, M.; REIS, F.J.C. Estabilidade das estruturas de concreto para solicitações combinadas. In: V Simpósio EPUSP sobre Estruturas de Concreto. Anais CD-ROM, p. 1-18, São Paulo, 2003. [ Links ]

[10] TAUCER, F.F.; SPACONE, E.; FILIPPOU, F.C. A fiber beam-column element for seismic analysis of reinforced concrete structures. EERC Report 91/17 - Earthquake Engineering Research Center University of California, Berkeley, USA, 1991. [ Links ]

[11] TEIXEIRA, M. R.; de SOUZA, R.M. Análise não linear física e geométrica de um edifício de múltiplos andares em concreto armado utilizando-se a plataforma OpenSees. In: V Simpósio EPUSP sobre Estruturas de Concreto, Anais CD-ROM, p. 1-19, São Paulo, 2003. [ Links ]

[12] PIMENTA, P. M. Análise não linear de pórticos planos. In: Anais EPUSP, v. 1, n. 1a, p. 563-582, São Paulo, 1988. [ Links ]

[13] PIMENTA, P.M.; SOLER, J.G.M. Estabilidade de pórticos planos de concreto armado. In: Anais do Simpósio EPUSP sobre Estruturas de Concreto, 1, v.2., p. 501-527, São Paulo, 1989. [ Links ]

[14] VECCHIO, F. J.; COLLINS, M. P. The modified compression field theory for reinforced concrete elements subjected to shear. ACI Journal , v. 83, n. 2, p. 219-231, 1986. [ Links ]

[15] PINTO, R. S. Análise não-linear das estruturas de contraventamento de edifícios em concreto armado. 2002. 155 p. Tese (Doutorado) - Curso de Engenharia Civil, Escola de Engenharia de São Carlos - USP, São Carlos, 2002. [ Links ]

[16] CARVALHO, M. F. M. S. Formulação corrotacional para análise de vigas com elementos finitos. 2010. 71 p. Dissertação (Mestrado) - Curso de Engenharia Mecânica, Universidade Nova de Lisboa, Lisboa, 2010. [ Links ]

[17] STRAMANDINOLI, R.S.B. Modelos de elementos finitos para análise não linear física e geométrica de vigas e pórticos planos de concreto armado. 2007. 189 p. Tese (doutorado) - Universidade Federal de Santa Catarina, Florianópolis, Brasil, 2007. [ Links ]

[18] STRAMANDINOLI, R.S.B. e LA ROVERE, H.L. An Efficient Tension-Stiffening Model for Nonlinear Analysis of Reinforced Concrete Members. Engineering Structures v.30, n.7, p.2069-80, Ed. Elsevier, 2008. [ Links ]

[19] STRAMANDINOLI, R.S.B. e LA ROVERE, H.L. FE model for nonlinear analysis of reinforced concrete beams considering shear deformation. Engineering Structures v. 35 p. 244-253, Ed. Elsevier, 2012. [ Links ]

[20] D'AVILA, V. M. R. Estudo sobre modelos de fissuração de peças de concreto armado via método dos elementos finitos. 2003. 259 p. Tese (Doutorado) - Universidade Federal do Rio Grande do Sul, Porto Alegre, 2003. [ Links ]

[21] COOK, R. D.; MALKUS, D. S.; PLESHA, M.E. Concepts and Applications of Finite Element Analysis. 3rd. ed. Ed. Jonh Wiley & Sons, Inc., 1989. [ Links ]

[22] RIKS, E. The application of Newton's method to the problem of elastic stability. Journal of Applied Mechanics, v. 3, p.1060-1065, 1972. [ Links ]

[23] WEMPNER, G. A. Discrete approximation related to nonlinear theories of solids. International Journal of Solids and Structures, v. 7, p.1581-1599, 1971. [ Links ]

[24] HOGNESTAD, E. A study of combined bending and axial load in reinforced concrete members. Bulletin Series, 399:128 - University of Illinois, Urbana, USA, 1951. [ Links ]

[25] COMITÉ EURO-INTERNATIONAL DU BÉTON. CEB-FIP Model Code 1990. London, Thomas Telford, 1993. [ Links ]

[26] FIGUEIRAS, J. A. Practical approach for modelling the nonlinear response of RC shells. Computational Modeling of Reinforced Concrete Structures, p. 217-253, 1986. [ Links ]

[27] LA ROVERE, H. L. Nonlinear analysis of reinforced concrete masonry walls under simulated seismic loadings. 1990. 200p. Ph.D. Dissertation (Structural Engineering) - University of California, San Diego, USA, 1990. [ Links ]

[28] BAZANT, Z. P.; PAN, J.; CABOT, G. P. Softening in reinforced concrete beams and frames. Journal of Structural Engineering (ASCE), v. 113, n. 12, p. 2333-2347,1987. [ Links ]

[29] BEBER, A. J. Avaliação do desempenho de vigas de concreto armado reforçadas com lâminas de fibra de carbono. 1999. 108p. Dissertação (Mestrado em Engenharia Civil) - Universidade Federal do Rio Grande do Sul, Porto Alegre, Brasil, 1999. [ Links ]

[30] JUVANDES, L. F. P. Reforço e reabilitação de estruturas de betão usando materiais compósitos de "CFRP". 1999. 400p. Tese (Doutorado) - Faculdade de Engenharia, Universidade do Porto, Portugal, 1999. [ Links ]

[31] PETERSON, A.; PETERSSON, H. On finite element analysis of geometrically nonlinear problems. Computer Methods in Applied mechanics and Engineering, v. 51, p. 277-286, 1985. [ Links ]

[32] LAZARO, A.L. e RICHARDS JR, R. Full-range analysis of concrete frames. Journal of the Structural Division (ASCE) , v. 99, n. 8, p. 1761-1783, 1973. [ Links ]

[33] SUN, C. H.; BRADFORD, M. A.; GILBERT, R. I. A reliable numerical method for simulating the post-failure behaviour of concrete frame structures. Computers & Structures, v. 53, n. 3, p. 579-589, 1994. [ Links ]

[34] BRATINA, S.; SAJE, M.; PLANINC, I. On materially and geometrically non-linear analysis of reinforced concrete planar frames. International Journal of Solids and Structures , v. 41, n. 24-25, p. 7181-7207, 2004. [ Links ]

[35] ERNST, G. C., et al.. Basic reinforced concrete frame performance under vertical and lateral loads. ACI Journal , v. 70, n. 4, p. 261-269, 1973. [ Links ]

[36] JUNGES, E. Estudo comparativo entre métodos simplificados e modelos de elementos finitos não lineares para o cálculo de flecha imediata em vigas de concreto armado. 2011. 360p. Dissertação (Mestrado) - Universidade Federal de Santa Catarina, Florianópolis, Brasil, 2011. [ Links ]

[37] GELATTI, F. Análise não linear física e geométrica de pórticos planos de concreto armado: modelagem por elementos finitos de barra. 2012. 226p. Dissertação (Mestrado) - Universidade Federal de Santa Catarina, Florianópolis, Brasil, 2012. [ Links ]

[38] JUNGES, E.; LA ROVERE, H. L. Comparison between simplified and FE models for short-term deflection in continuous RC beams. Ibracon Structures and Materials Journal. (to be published) [ Links ]

7. ANNEX

The interpolation functions 𝑁(𝜉), as used in section 2.1, are given by:

From the displacement field, the strain in the element can be determined, by considering equation (2):

where:

and:

in which: ``= 𝑑 𝑑𝜉 and 𝜑 𝜉 is the curvature.

Hence the curvature can be written as follows:

And thus equation (2) can be rewritten as:

or:

The secant stiffness matrix of the element, according to section 2.2, is given by:

Received: December 09, 2015; Accepted: June 20, 2016