A 3 D finite element model for reinforced concrete structures analysis

loads. The proposed model for concrete is orthotropic and uses the equivalent uniaxial strain concept. The equivalent uniaxial stress-strain relation is generalized to take into account the triaxial stress conditions. The parameters used in the equivalent uniaxial stress-strain curve are determined from the failure surface defined in the principal stress space. The implementation in finite elements is based on the consideration of smeared cracks with cracks rotating according to the directions of the principal stresses. Also, an embedded reinforcement model was implemented to represent existent reinforcing bars. Finally, some results are compared with experimental data from the literature to demonstrate the validity of the numerical model developed.


Introduction
The implementation of finite element models to investigate reinforced concrete structures has been the subject of many studies in the last years ( [8]).In several of these studies, a fundamental matter has been how to precisely define concrete's behavior when a multiaxial state of stress is in place.This task has been difficult because of the behavioral complexity of the material due to a series of factors, among which, these can be pointed out: concrete presents a nonlinear stress-strain relation; cracks tend to form when concrete is subjected to tensile stresses; concrete is subjected to creep and shrinkage; concrete's behavior is dependent of its load history; and concrete presents a volumetric expansion when near collapse.In the last decades, numerous constitutive models were developed for the analysis of concrete structures.However, most of them, in general, can describe only certain aspects of concrete's behavior, severely limiting their application.Among the most important characteristics a model for concrete must present, the following can be highlighted: to allow a general and precise representation of its real mechanical behavior; to employ a formulation that is clear to understand; and to make viable the implementation of a robust algorithm, stable enough for a nonlinear behavior determination.In this context, the objective of this work is to present a computational model to analyze the behavior of reinforced concrete elements subjected to multiaxial stress states.The constitutive models adopted for concrete and steel are herein described in details, and have been implemented in a computational code, developed by Bono [9], to analyze reinforced concrete structures through the Finite Element Method.Linear and quadratic isoparametric hexahedral elements can be chosen to model reinforced concrete structures, with the embedded reinforcing model proposed by Elwi and Hrudey [10] extended to the 3D case being used to represent the existent reinforcing bars.

Constitutive model for the concrete
In this work, an orthotropic nonlinear elastic constitutive model is used to represent the concrete's behavior.Among the models that have been already developed in this category, there is the one presented by Darwin and Pecknold [11] for plane stress state analyses.Elwi and Murray [12] continued the development of this model, extending its applicability to the 3D case.Afterwards, the model was further developed by Balan et al. [13], Balan et al. [14], Kwon [15], and Kwon and Spacone [16].The orthotropic model proposed is based on the work by Kwon [15], being capable of simulating the response of concrete under a multiaxial stress state.Nevertheless, some modifications have been implemented in the original model in order to further improve it, which are appropriately commented later along with the model's formulation.
The shear moduli, G ij , can be expressed by: (4) ( ) ( ) If the orthotropic axes are considered parallel to the directions of the updated principal stresses, then the constitutive relation can be reduced to: (5)

A 3D finite element model for reinforced concrete structures analysis
The equivalent uniaxial strains are used for determination of concrete's properties, i.e., the secant moduli of elasticity and the Poisson's ratios, which are used in Equation ( 5).Nevertheless, it can be noticed from Equation (6) that, to obtain the equivalent uniaxial strains for a nonlinear material, the secant moduli of elasticity, E i , are needed in the directions of the three principal stresses, being necessary the use of an iterative process for determination of these variables (see section 2.8).

Equivalent uniaxial stress-strain curves
As mentioned before, the adopted model for concrete follows the work of Darwin and Pecknold [11], admitting that the 3D constitutive law can be separated in three uniaxial relations, with actual stresses being functions of equivalent uniaxial strains.Elwi and Murray [12] used an expression by Saenz [17] to describe equivalent uniaxial stress-strain relations, which was both adopted for concrete's compressive and tensile response (see Figure 1): where In relation (5), the values of the secant moduli of elasticity, E i , and the Poisson's ratios, ν ij , have to be determined.These values are obtained from the uniaxial stress-strain curves for the concrete using the concept of equivalent uniaxial strain, which is detailed next.

Equivalent uniaxial strain
In a multiaxial stress state, the actual strain in a given direction is not only function of the stress in that direction, but also of the stresses acting in the other orthogonal directions.This makes the determination of the 3D response of concrete a complex task to perform.In order to make it an easier task, the concept of equivalent uniaxial strain can be used, a concept originally introduced by Darwin and Pecknold [11] for 2D analyses.Darwin and Pecknold [11] considered this procedure to be an ingenious scheme to separate the 2D response of concrete in two uniaxial curves, making it easier to determine the material's behavior.The technique provided an alternative to separate the effects of the Poisson's ratios of each strain, allowing a good approximation of experimental data.The equivalent uniaxial strains can be determined by: (6) where σ i -is the current principal stress in the orthotropic direction i; E i -is the secant modulus of elasticity in the orthotropic direction i, with i = 1, 2, 3 for 3D analyses.Kwon [15] also used the concept of equivalent uniaxial strain for 3D analyses in concrete, rewritting relation (5) in the following form: (7) and the equivalent uniaxial strains are then given by: (8) The principal stresses σ 1 , σ 2 , and σ 3 can be obtained from the uniaxial constitutive relations [Equation (6)], using the equivalent uniaxial strains ε u1 , ε u2 , and ε u3 .Therefore, the introduction of these equivalent uniaxial strains allows a representation of the triaxial behavior of the concrete through three independent stress-strain curves.
It is worth mentioning that these equivalent uniaxial strains are not actual strains, but fictitious strains defined in the updated directions of the principal stresses, being accumulated in these directions.Thus, ε ui does not give a strain history in a fixed direction, but rather in an always changing direction with the updated state of principal stresses (Darwin and Pecknold [11]). -- The other variables were defined in Equations ( 9)- (10).The use of the Popovics-Saenz curve (Equation ( 11)) as an equivalent uniaxial relation gives excellent results for plain concrete.Nevertheless, for reinforced concrete elements it has not shown to be adequate to take into account the collaboration of concrete between cracks (tension-stiffening). Therefore, for the proposed model herein, the Popovics-Saenz curve is used to describe only the response in compression of concrete under monotonic loads, as illustrated in Figure 2, and, for the response in tension, the formulation described next is used.When ε ui ≤ε ci : (12) E o -is the initial modulus of elasticity; ε ui -is the equivalent uniaxial strain in the orthotropic direction i; f ci -is the concrete strength or the peak stress; ε ci -is the strain that corresponds to f ci , i.e., it is the peak strain.f fi , ε fi -are the stresses and strains at the control point in the curve's descending branch.
The proposed curve by Saenz [17] has been considerably used as a stress-strain relation for concrete, with only one expression representing both ascending and descending branches.However, this curve works well only when , when the secant modulus of elasticity at the peak point, E ci , is not larger than half of the initial modulus of elasticity, E o .When this condition is not satisfied, the curve will present two curvatures between the origin and the peak point.This problem can be partially corrected by fixing E o /E ci to a value of 2, independently of the actual relation between the two moduli (Balan et al. [14]).
Popovics [18] proposed another curve to define concrete's stressstrain relation, which is expressed by: (10) ( ) where R i = K i / (K i -1) and the other variables have already been defined in Equation ( 9).This curve approximates well the actual initial stiffness of concrete in the ascending branch.However, it does not adjust well to experimental data when descending.
To avoid any limitation in the definition of the equivalent uniaxial curve, Kwon [15] proposed a modification in the curve used by Elwi and Murray [12]: two curves are used to describe the responses in tension and in compression of concrete.The curve proposed by Popovics [18] should describe the ascending branch until the peak point, while the curve proposed by Saenz [17] should be used in the response of the descending branch.The combination of the models by Popovics [18] and Saenz [17], named as Popovics-Saenz curve, from Equations (9) to (10), can then be expressed by the following relation: (1 ) 0.01 where α t -is a coefficient for reduction of cracking stresses, and the other variables were defined in Equation (9).The variable ε ctu , introduced in Figure 2, indicates the limiting strain for which the collaboration of concrete between cracks should not be considered any longer.In this study, in order to allow a better agreement with experimental data, a value of 0.01 has been adopted for the strain ε ctu , as indicated in Equation ( 13), while the value of parameter α t depended on the element analyzed.
In Equations ( 11) to ( 13), the expression proposed in CEB-FIP Model Code 1990 [19] was adopted for the initial modulus of elasticity, E o .Through those three equations, the secant moduli of elasticity, used in Equation ( 6), could be defined by considering

Failure surface of concrete
The determination of the variables K i , K εi , K σi , used in Equation (11), which are functions of the peak stresses and strains, f ci and ε ci , is necessary for definition of the three equivalent uniaxial curves.This peak stresses and strains, for a determined current stress state, are calculated from failure surfaces in the space of principal stresses (σ 1 , σ 2 , σ 3 ).Numerous formulations have been proposed as failure surfaces for concrete.A series of failure criteria was presented by Chen and Han [20] and by Menétrey and Willam [21], classifying these criteria according to the number of parameters that appear in the respective expressions.Among the more used failure surfaces for description of con-crete's triaxial strength, there are the four-parameter surface by Ottosen [22] and the five-parameter surface by Willam and Warnke [23].These surfaces have been highly used because the failure points they determine are quite close to available experimental data, as can be seen in Chen and Han [20].Therefore, for the present model, the failure surfaces by Ottosen [22] and Willam and Warnke [23] were implemented for concrete.
The equations used for these two surfaces are presented in details in the next sections.

The failure surface by Willam and Warnke
This failure surface can be expressed by the following equation: ( where 3 oct r τ = -is the stress component perpendicular to the hydrostatic axis; ρ f (σ m , q) -defines the failure curve on deviatoric planes, being determined by: , 2 cos The five-parameter surface by Willam and Warnke [23] presents parabolic curves for tensile and compression meridians (Figure 3) expressed by: (16) where σ m -is the mean normal stress; ρ t , ρ c -are stress components perpendicular to the hydrostatic axis for q = 0 o and q = 60 o , respectively; From Equation ( 16), the stress components perpendicular to the hydrostatic axis, ρ c and ρ t , are determined by: ( ) for the tensile meridian.
( ) for the compression meridian.
Since the two meridians should cross the hydrostatic axis at the same point, results that a 0 = b 0 .Based on biaxial tests by Kupfer et al. [24] and other triaxial tests, the five parameters of the failure surface by Willam and Warnke [23] can be determined from the following values on the failure surface (Chen and Han [20]): f c -uniaxial compression strength (q = 60 o ); In the present study, however, another approach, using the expressions for f t and f cc recommended by the CEB-FIP Model Code 1990 [19] have been adopted.In this approach, the uniaxial tensile strength, f t , is determined by:

The Failure Surface by Ottosen
The failure criterion proposed by Ottosen [22] defines that the failure surface for concrete subjected to multiaxial stress states is given by: (20) for cos 3q < 0; -is the second invariant of the deviatoric stress tensor; The four parameters, a, b, c 1 , and c 2 , used in this failure surface can be determined from the following conditions: f c -uniaxial compression strength (q = 60 o ); f t -uniaxial tensile strength (q = 0 o ); f cc @ 1.16 f c -biaxial compression strength (q = 0 o ); (x, r) = (-5f c ; 4f c ) -triaxial state in the compression meridian (q = 60 o ).

Determination of peak stresses (f c1 , f c2 , f c3 )
From the failure surface by Willam and Warnke [23] or from the one by Ottosen [22], the peak stresses, f c1 , f c2 , and f c3 , regarding the three directions for a determined current principal stress state (σ 1 , σ 2 , and σ 3 ), can be found.To calculate the values of these peak stresses, the following procedure has been used: 1.A straight line from the origin of the principal stress system to the point of current stresses M c (σ 1 , σ 2 , σ 3 ) is determined; 2. Next, this straight line is extended until the failure surface is reached, identifying the point M r (f c1 , f c2 , f c3 ), as can be seen in Figure 5.The equation of the straight line that goes from the origin to the updated principal stress state (σ 1 , σ 2 , σ 3 ) can be expressed as a function of the octahedral stresses σ oct and τ oct : Numerically, the procedure is carried out by solving a system of equations that corresponds to the determination of the intersection of the line, Equation (22), with the failure surface by Willam and Warnke [23], Equation ( 14), or with the failure surface by Ottosen [22], Equation (20).These calculations, in this case, were carried out by the subroutine DNEQNF, available in IMSL FORTRAN 90 MP Library [25] to solve the nonlinear system of equations.By solving the system of equations, the octahedral stresses ( r oct σ , r oct τ ) at the intersection point of the line with the specified failure surface is obtained, and with these octahedral stresses, the peak stresses are determined by the following expressions: (23) In Figure 5, it can be noticed that the failure surface has its origin in a failure point of hydrostatic tension and opens up in the negative direction of the hydrostatic axis.Therefore, a compressive hydrostatic load can not cause failure (Chen and Han [20]).With increasing lateral confinement, the updated stress state (σ 1 , σ 2 , σ 3 ) gets closer to the hydrostatic axis.In this situation, depending on the triaxial state of compression stresses, Equation ( 22) used to define peak stresses could not intercept the failure surface, not allowing the determination of values for f c1 , f c2 , and f c3 .In this case, a second alternative is used for definition of the ultimate stresses, where the octahedral normal stress, σ oct , is considered constant, while the octahedral shear stress, τ oct , should vary until reaching the failure surface.This procedure is similar to the previous one, only replacing Equation ( 22) by: (24)

Control Point of Popovics-Saenz curve
It is known that concrete exhibits a softening in its stress-strain curve after peak stresses are reached.This behavior is a phenomenon that depends on specimens' dimensions and shape, as well as on boundary conditions existent in the tests conducted.Moreover, the post-peak branch of the curve presents considerable variability, depending on test characteristics (Balan et al. [14]).Therefore, to perform a realistic analysis, the definition of an appropriate model to describe the material's post-peak behavior is essential.With Popovics-Saenz curve, an adequate estimative of concrete's compression response can be reached by using a control point located in the postpeak branch of the equivalent uniaxial stress-strain curve.According to Balan et al. [13], from experimental tests, the following values are proposed for the stresses and strains that define this control point: (28) f fi = 0.85 f ci e fi = 1.41 e ci (29) f fi = 0.25 f ci e fi = 4.0 e ci where f ci -is the material strength in the i direction, with i = 1, 2, 3; ε ci -is the strain that corresponds to f ci .
Increases in confinement stresses in triaxial compression tests were verified to affect strength, ductility, dilatancy, and failure modes of concrete.In addition, increases in lateral confinements make the failure mode of concrete change from brittle to ductile (Balan et al. [14]).However, the established values in ( 28) and ( 29) define the softening surface of the material as a surface that has its shape proportionally reduced in relation to the failure surface, independently of loading conditions.Thus, for a better consideration of concrete's experimental behavior when subjected to compression, the stress values f fi of the control point are defined as function of the peak stresses, which are dependent on the confinement stresses applied to the structure (Balan et al. [13]).
From the experimental results by Smith et al. [27], Kwon [15] defined the following coordinates for the control point: (30) where, f c is concrete's uniaxial compression strength.As can be noticed, in this work, the above relations have also been adopted for the control point used on the definition of the equivalent uniaxial stress-strain curve of the concrete subjected to compression.

Cracked and Crushed Concrete
In this work, a smeared crack model has been implemented, with cracks rotating according to the directions of updated principal stresses.No special treatment is given to post-peak behavior of concrete after crushing or cracking, since the equivalent uniaxial laws entirely govern the response of the material.The criterion adopted to verify cracking or crushing occurrence was to check the equivalent uniaxial strain, in a specific direction, to see if it is greater than the corresponding strain to the peak stress, also considering the signal of the ongoing stress.

Poisson's ratio
The values of the Poisson's ratios, ν ij , have to be determined to complete the definition of the constitutive law represented by Equation (5).It is known that the volume of concrete in a specimen under A 3D finite element model for reinforced concrete structures analysis compression initially decreases and, later on, increases when failure is imminent.This dilatancy occurs because of the opening (or expansion) of existent microcracks in the material.,Several works from the literature, like the ones published by Elwi and Murray [12], Balan et al. [13] and Kwon [15], suggest the use of an increasing function for the Poisson's ratio to simulate that volumetric variation.In this work, the expression proposed by Kwon [15] is the one used to define Poisson's ratios: (31) where ν ui is the Poisson's ratio for the equivalent uniaxial direction i.The symmetry of the constitutive matrix, D o , which was expressed in Equation (3), is guaranteed when Equation (31) is considered.For definition of ν ui , the following expression is used: (32) where ν o -is the initial Poisson's ratio;  32) is similar to the cubic expression defined by Elwi and Murray [12], which later on was used by Balan et al. [13].Kwon [15] proposed the addition of the coefficient K n to that cubic function to avoid Ω c , which appears in Equation (3), to become close to zero or a negative number, since a zeroed value for Ω c would create numerical problems in the finite element formulation of the constitutive law.

Implementation of the Numerical Model for the Concrete
In this section, the determination of the stresses in the concrete is detailed.Initially, an isotropic behavior is admitted for the concrete, considering the secant moduli of elasticity equal to the initial modulus of elasticity, i.e., . With these secant moduli of elasticity and the current strains for a determined load step, the total stresses are calculated, as presented in Table 1.With these total stresses related to the global system of reference, the current principal stresses, σ 1 , σ 2 , σ 3 , are determined.From these principal stresses, the equivalent uniaxial strains for the three principal directions are determined.With the principal stresses and the failure surface chosen for the concrete, the peak values of stresses and strains are also determined, as explained in 2.4.3 and 2.4.4.Next, the stresses and strains of the control point used in the descending branch of the Popovics-Saenz curve [Equation (11)] are also obtained.The new equivalent uniaxial stress-strain curves are defined for the three directions of principal stresses with the determination of these variables.From these three equivalent uniaxial curves, the new values of the current principal stresses are calculated.Finally, from the principal stresses, the total stresses related to the global system of reference are updated.In an iterative process, all the described steps are repeated until convergence, as indicated in Table 1.

Constitutive model for the steel
Since reinforcing bars resist to forces that primarily develop along their axial direction, i.e., perpendicular forces to rebars' axes can be neglected, it is enough to know rebars' properties relatively to an uniaxial stress state.The shape of the stress-strain diagram for steel can be represented by four branches: an elastic branch; a yielding branch; a hardening branch; and a softening branch.In the literature, however, there are several simplified representations for that uniaxial stress-strain behavior: elastic-perfectly plastic model; elastic-linear work hardening model; trilinear approximation; and a complete stress-strain curve.In this work, a bilinear diagram with work hardening was adopted for the monotonic stress-strain curve, as defined in Figure 6 and by: (33)  A 3D finite element model for reinforced concrete structures analysis

Numerical Application
In this section, the results obtained with the computational model developed are compared with values determined experimentally for reinforced concrete beams.These experimental results were presented by Bresler and Scordelis [28], concerning twelve beams in reinforced concrete primarily subjected to shear.The twelve beams were arranged in four series of three units, with each series concerning a different quantity of longitudinal reinforcement, quantity of stirrups, span length, cross-sectional dimensions, and concrete's strength.All beams had rectangular cross-sections and their details can be seen in Figure 7, with additional information presented in Tables 2 and 3.All bottom longitudinal bars were 28.7 mm in diameter, while bars of 12.7 mm were used at the top of the beams, and, when present, stirrups were made of 6.4 mm bars.One of the series had no stirrups at all, and was labeled as OA series.All beams were subjected to monotonic concentrated loads applied at midspan, as illustrated in Figure 8.
For the computational analysis, a finite element discretization of ten quadratic hexahedral elements was adopted to model the beams.Taking advantage of the problem's symmetry, the discretization and boundary conditions were the ones shown in Figure 9. Load-displacement curves, as presented in Figures 10 to 12, were used for validation of the numerical analysis.In Table 4, the failure loads obtained with the numerical program are presented and can be compared with the experimental results also indicated.In a general manner, an excellent correlation was obtained between the numerical response obtained and the experimental results given by Bresler and Scordelis [28].The failure loads obtained with the numerical model are quite close to the experimental responses for most of the beams analyzed.

Conclusions
This work presents a general formulation for nonlinear analyses through finite elements of reinforced concrete structures.The constitutive law used for concrete is orthotropic with its axes parallel to the updated principal stresses' directions.The model also uses the concept of equivalent uniaxial strain; a concept originally pre-sented by Darwin and Pecknold [11], for determination of ultimate stresses, and allows the use of two 3D failure criteria, namely the criterion by Willam and Warnke [23] and the one by Ottosen [22].The model describes well the response of concrete subjected to several loading types, being capable of representing crushing and cracking of concrete.The concept of smeared cracking was used for consideration of cracked concrete, while the modified Popovics-Saenz curve was used.

Figure 1 -
Figure 1 -Illustration of the model proposed by Saenz 1964

Figure 2 -
Figure 2 -Stress-strain curves for concrete subjected to monotonic loads: (a) response in compression; (b) response in tension
95 £ a £ 1.85 (MPa), while the biaxial compres- sion strength is given by f cc = 1.20 f c .Since, according to the Model Code, the tensile strength of concrete is more variable than its compression strength and can be substantially reduced by environmental factors, and after comparisons with experimental data, then the maximum value of a (1.85 MPa) has been considered.From the rupture data previously mentioned, expressions for a 0 , a 1 , a 2 , b 0 , b 1 , b 2 , all functions of α u = f t / f c , have then been determined.These expressions are given by: first invariant of the stress tensor.

Figure 4 -
Figure 4 -The model by Ottosen 1977: tensile and compression meridians

Figure 6 -
Figure 6 -Bilinear constitutive model for steel

Figure 7 -
Figure 7 -Beams' cross sections for the tests by Bresler and Scordelis 1963

Figure 8 -Figure 10 -
Figure 8 -Schematics of the beams tested by Bresler and Scordelis 1963