3D Fiber-Based Frame Element with Multiaxial Stress Interaction for RC Structures

A three-dimensional fiber-based frame element accounting for multiaxial stress conditions in reinforced concrete structures is presented. )e element formulation relies on the classical Timoshenko beam theory combined with sectional fiber discretization and a triaxial constitutive model for reinforced concrete consisting of an orthotropic, smeared crack material model based on the fixed crack assumption. Torsional effects are included through the Saint-Venant theory of torsion, which accounts for out-ofplane displacements perpendicular to the cross section due to warping effects. )e formulation was implemented into a forcebased beam-column element and verified against monotonic and cyclic tests of reinforced concrete columns in biaxial bending, beams in combined flexure-torsion, and flexure-torsion-shear.


Introduction
e increasing interest towards performance-based engineering requires robust and accurate numerical models for strength and ductility assessment of reinforced concrete (RC) structures.A great majority of RC structures, such as bridges and frame buildings, can be idealized with onedimensional frame elements for global nonlinear analysis.However, these elements usually rely on the classical Euler-Bernoulli beam assumption, that is, plane sections remain plane and orthogonal to the element axis, and uniaxial stress conditions, which is only valid in the so-called B-regions.As an example, Figure 1 shows the results from such an element applied to the modeling of a floor subassembly subjected to a vertical concentrated load in the middle of the floor beam [1]. is creates multidirectional shear flow in the spandrel beam coupled with normal stresses arising from bending.e frame element is unable to reproduce the cracked torsional stiffness of the spandrel beam, thus significantly overestimating the overall cracked stiffness and strength of the test.In the same figure, the results from 3D nonlinear finite element analysis (FEA) have been included, demonstrating much better agreement [2].
From the above, it is clear that standard flexural frame elements, commonly used in commercial software, present limitations under multiaxial loading conditions, which have motivated the development of more advanced frame elements for global nonlinear static and dynamic analysis.Several formulations have been proposed in 2D addressing the problem of axial-flexure-shear interaction.Some authors used an equivalent truss model in shear combined with the Euler-Bernoulli beam assumption and uniaxial stress response in flexure [3].Further improvements were obtained when shear deformations were included in the element kinematics (Timoshenko beam theory) [4], which also allowed using two-dimensional constitutive models at each individual fiber based on either fixed [5][6][7] or variable shear strain profile [8].Different constitutive models have been used, ranging from orthotropic smeared crack models (e.g., modified compression field theory (MCFT) [9] and fixed angle softened truss models (FA-STM) [10]) to softened plasticity damage models [11,12], and microplane models [13], among others.A detailed review can be found in [14].
Research devoted to RC 3D frame elements has been more limited mostly due to complexities in the development of robust and accurate triaxial material models for concrete and their implementation into fiber frame elements.Mazars et al. [15] developed a 3D displacement-based Timoshenko beam element with a 3D constitutive model for concrete based on continuum damage mechanics.e Saint-Venant torsion was introduced by means of the warping function derived from linear elastic analysis under pure torsion.Gregori et al. [16] used a similar element to that used by Mazars et al. but with a parabolic shear strain profile and Coulomb theory for torsion.
e constitutive model was based on the MCFT.e cross section was subdivided into three regions based on the dominant stress conditions: 1D, 2D, and 3D, where the corresponding constitutive models were applied.Mullapudi and Ayoub [17] formulated a 3D flexibility-based fiber element using the FA-STM.e beam kinematics was based on the Timoshenko beam theory in conjunction with the Coulomb torsion, that is, no warping was considered.
e abovementioned work was based on the classical Timoshenko beam theory with either uniform or parabolic shear profiles.is approach does not satisfy interfiber sectional equilibrium; hence, more advanced (higher-order) beam theories have been recently investigated which satisfy interfiber sectional equilibrium in an average sense [18][19][20][21].However, only a limited number have been successfully applied to 3D analysis of RC frame structures [22].
In the present work, a nonlinear frame element based on the Timoshenko beam assumption is developed.e cross section is discretized into several fibers, which allows direct consideration of multiaxial stress interaction at the material level.A triaxial smeared crack constitutive model is implemented for concrete, whereas reinforcement is treated as smeared within the concrete with uniaxial stress-strain behavior.
e Saint-Venant theory is used for modeling torsion, which is suitable for noncircular solid cross sections with warping effects.
e element formulation has been implemented into a force-based frame element and applied to monotonic and cyclic analysis of RC columns and beams.Comparison against experimental results shows similar levels of accuracy to those obtained from detailed 3D FEA with solid elements.

Frame Element Formulation
In the Timoshenko beam theory (TBT), cross sections are not constrained to remain orthogonal to the beam axis, as in the Euler-Bernoulli beam theory, due to shear distortions of an infinitesimal beam element.e kinematics of any point of the beam is completely described by three independent displacements and rotations of the beam axis as follows (Figure 2): e total displacement at any given point of the section is related to the beam axis displacements as where u, v, and w are the axial and transverse displacements of a generic point of the section with coordinates (x, y, z), ω is the warping function, and α is the warping parameter equal to zθ x /zx according to the Saint-Venant theory.Assuming small displacements, axial and shear strains at a given section fiber are given as e section forces consistent with the section strains are obtained from integration of the axial and shear stress distributions over the cross section as follows: e relation between section forces and section strains is given by the section stiffness matrix.In the nonlinear case, incremental section forces and strains are related by the tangent stiffness matrix of the section as follows: where D is the condensed (3 × 3) tangent material stiffness matrix.Note that B σ and B ε are different; hence, K s will not be symmetric in the general case, with the exception of the circular cross section for which the warping function is zero.e D matrix is obtained from condensation of the (6 × 6) material stiffness matrix assuming that the total stresses σ y , σ z , and τ yz are zero in the beam element as follows: where σ c,i and σ s,i are the stresses in concrete and steel, respectively, and ρ si (i � y, z) is the reinforcement ratio at each individual section fiber. is process is performed iteratively at each section fiber until convergence, which requires implementing an additional loop within the fiber state-determination procedure.e warping function in ( 2) is obtained from the solution of an elastic prismatic beam under pure torsion based on the Saint-Venant theory of torsion, which is defined by the following boundary value problem: where Ω and zΩ represent the cross-sectional domain and its boundary, respectively, and n y and n z are the components of the vector normal to the boundary.Equations ( 9) and ( 10) represent a boundary value problem that can be solved using the finite element technique.In the present case, four-node plane isoparametric elements were used to discretize the cross section and calculate the elastic warping function ω.
Figure 3 shows the calculated warping functions for the rectangular and hollow cross sections used for verification in Section 5.

Constitutive Model
A 3D orthotropic smeared crack material model is formulated for cracked reinforced concrete, which represents an extension of the 2D model presented in [23] for membrane elements.A similar approach was also adopted in [17,24].e model is formulated in the material directions (Figure 4), which are coincident with the crack directions determined at first cracking as where ε cpi is the plastic compressive strain (positive) and ε i is the principal strain (negative for tension).ε i can be obtained upon solving the eigenvalue problem for the strain tensor in the xyz reference system, ε xyz .After the first cracking, the material directions are kept constant throughout the analysis, assuming that secondary cracks are orthogonal to the initial ones.Total strains in the material directions, ε 123 , are obtained from ε xyz as where T is the standard (6 × 6) transformation matrix.e corresponding concrete stresses, σ 123 , are determined from equivalent uniaxial constitutive laws for axial and shear stresses (Figure 5(a)).In compression, the monotonic envelope is defined by the Popovics curve scaled by the lateral confining factor and a softening coefficient due to orthogonal tensile strains [9] given as where ε j and ε k are the total tensile strains in the orthogonal directions j and k.

Advances in Civil Engineering
In tension, a physically based approach is used such that average concrete stresses are determined using equilibrium, compatibility, and constitutive relations applied on a cracked concrete element with three di erent regions as follows: an elastic, fully bonded region (b), a crack region (cr), and a bond-slip region (sl) where the slip between concrete and reinforcement occurs (Figure 5(b)).us, average concrete stress in the ith material direction can be obtained as where σ b ci , σ sl ci , and σ cr ci are concrete stresses in the (b), (sl), and (cr) regions, respectively.s i , λ i , and ε i are the average spacing, unbonded length parameter, and tensile strain in the ith direction, respectively.is approach allows full characterization of the tension sti ening, crack-closing, and crack-opening regimes based on principles of mechanics.Further details can be found elsewhere [23,25].
Shear stresses, τ ij , arising from deviation between principal directions and crack directions are related to the shear strains, c ij , as where σ i , σ j and ε i , ε j are the normal stresses and strains in the ith and jth crack directions, respectively.Note that, in this approach, dowel action is not directly accounted for.Shear failure due to aggregate interlock is included by setting an upper limit on the maximum shear stress transferred along cracks [9] as follows: where a g is the maximum aggregate size and w ij is the average crack width (mm).Assuming smeared cracks, the average crack width is given as where ε n is the maximum (in absolute value) tensile strain normal to the crack.e average crack spacing, s ij , is taken as the average between i and j directions as follows: where the average crack spacing in the ith direction is determined from those in the x, y, and z directions as follows: where n ji are the direction cosines between xyz and the crack directions.
Finally, steel stresses are determined based on the uniaxial Menegotto-Pinto hysteretic model [26] with modi ed isotropic hardening [27].For the ith reinforcement group, which orientation is de ned by the unitary vector n si , the total axial strain can be obtained upon transforming ε xyz to the corresponding direction.en, steel stresses are obtained using the abovementioned model and transformed back to the xyz reference system as where σ si is the uniaxial steel stress, ρ si is the reinforcement ratio for the ith reinforcement group, and T si is the corresponding transformation matrix between reference systems.σ si,xyz stresses are added to the concrete stresses, σ c,xyz , in order to obtain the total stresses, σ xyz , which are used for the calculation of the vector of internal element forces.
Similarly, the ith reinforcement contribution to the material  Advances in Civil Engineering stiffness matrix, D si,xyz , is added to that provided by concrete, D c,xyz , such that the total material stiffness matrix is given as where where E ci and E si are the tangent moduli for concrete and steel and G ij is the shear stiffness.

Frame Element Implementation
e proposed frame model was implemented using Matlab © into a three-dimensional, force-based fiber element based on the state-determination formulation described by Spacone et al. [28].In this formulation, once the global displacements are known, an iterative procedure is invoked at the element level in order to obtain the corresponding section strains at each integration point.is is performed in the basic element configuration upon removal of the rigid body modes as follows (Figure 6): where Δq is the increment of element displacements in the basic configuration, Δd E is the increment of global displacements, at a given load step, and R is the rigid body transformation matrix given as e element forces, ΔQ, can be obtained from the element flexibility matrix, F e , which is known from the previous converged load step as Equation ( 25) is only an approximation since F e is changing within each load step due to material nonlinearities.Section forces, ΔS, are interpolated exactly from element forces.At the kth integration point, this is given as where b k is the force interpolation matrix at the kth integration point with coordinate x k , which is given as An estimate of the section deformations at the kth integration point, Δe k , can be made based on the section flexibility matrix, F s,k , which is obtained from the inverse of the section stiffness matrix given in (7): At this point, the total strains at each individual fiber can be determined, and the corresponding stresses can be obtained using the constitutive models presented in Section 3. Integration of the total stresses over the cross section yields the section forces, ΔS int,k , and an updated section flexibility Advances in Civil Engineering matrix.e unbalanced section forces, ΔS u,k , given by the difference between ΔS k and ΔS int,k , are used to update the section deformations as Δe k � Δe k + F s,k ΔS u,k , which are further used to correct the element forces by the following quantity: where w k are the integration weights at the kth integration point and F e is the updated element flexibility matrix given as e procedure is repeated with the new element forces until satisfying element convergence, that is, when ‖ΔS k − ΔS int,k ‖ < tol at each kth integration point.Once element convergence is achieved, the element forces obtained in the basic configuration are transformed back to the global configuration and assembled within the vector of internal forces of the structure.Global equilibrium between the vector of applied external forces and internal forces is checked at this point.If this is not satisfied, the residuum between external and internal forces and the tangent stiffness matrix of the structure are used in a subsequent iteration to compute a new increment of nodal displacements.A flow chart summarizing the numerical procedure is shown in Figure 7.

Verification Examples
5.1.Biaxial Bending Test.Bousias et al. [29] tested flexuredominated RC cantilever columns in biaxial cyclic bending.Different paths of imposed lateral displacement were investigated, ranging from uncoupled bidirectional to circular patterns.All columns were 250 mm × 250 mm × 1490 mm rectangular specimens fixed to a stiff base foundation.e shear span-to-depth ratio was 6.
e columns were well reinforced in shear with a double hoop pattern of 8 mm diameter stirrups at 70 mm spacing.Longitudinal reinforcement consisted of 16 mm bars uniformly distributed around the perimeter.Material properties are summarized in Table 1.A constant axial load was applied at the top of each column.Two tests are reproduced here, S1 and S7, with the displacement paths shown in Figure 8.
e columns were modeled with one single FB element with 5 IPs.e cross section was discretized into 10 × 10 fibers, with the longitudinal and transverse reinforcement smeared as shown in Figure 9.
e contribution of the inclined hoops was decomposed in the x and z local directions and added to the hoops aligned with the axes.In any case, the influence of transverse reinforcement in the response of the member is rather small since the behavior is dominated by flexure.
e main effect of the transverse reinforcement is confining of concrete, which was accounted for with an estimated confinement factor of 1.1.
e hysteretic lateral force-displacement response is shown in Figure 10.e numerical results agree well with the experimental data for both S1 and S7 specimens.Wide hysteretic loops can be observed for specimen S7 due to the nature of the imposed displacement.
is aspect is well captured by the analytical model.For specimen S1, analytical results slightly underestimate the strength of the member, which might be related to material modeling of steel, although good agreement is observed in terms of hysteretic pinching and residual displacements.Figure 11 also reports the measured horizontal reactions in each direction.It is noted that the horizontal reactions are not totally decoupled for the S1 column, as it happens with displacements, presumably due to biaxial yielding.It shall be mentioned that testing and measuring the response of biaxially loaded columns involves several difficulties, such as misalignment of the cross-sectional principal axes with those of the actuator, and inclusion of P-delta effects, among others.

Torsion and Bending.
e series of RC beams tested by Onsongo [30] in combined torsion and flexure are analyzed herein.
e test campaign consisted of two series of five beams each, being the series TBO over-reinforced and designed to fail in concrete crushing without yielding of neither longitudinal nor transverse reinforcement, whereas the series TBU was under-reinforced, with either longitudinal or transverse reinforcement yielding prior to failure.Within each series, the ratio of applied torque to applied moment was varied from 0 to 5.059, being the rest of the test variables essentially constant.
e test setup consisted of adjustable steel level arms at the ends of the specimen, where

6
Advances in Civil Engineering an eccentric jack load was applied, resulting in constant bending moment and torque along the specimen.Beams TBO3 and TBU3 were selected for analysis, both having the same cross section.Material and loading properties are summarized in Table 2.
Half of the beam was modeled using one single element with 5 IPs (Figure 12).Rigid elements were connected to the end of the beam in order to reproduce the support conditions and the eccentric load application.e torsional and bending rotations around z-axis were restrained at node 1.At node 4, the beam was free to rotate, and only the vertical displacement was restrained.Additionally, out-of-plane degrees of freedom were restrained in order to create an isostatic system and avoid stiffness matrix singularities.e Table 1: Material properties of the columns [29].Advances in Civil Engineering cross section was discretized into 64 fibers, with both longitudinal and transverse reinforcements smeared within the web and flange regions as shown in Figure 12. e warping function of the hollow section is already presented in Figure 3.It was calculated as the difference between the warping functions of two solid rectangular sections, one with the outer dimensions and the other one with the dimensions of the inner core.
Figures 13 and 14 compare analytical and experimental results in terms of applied bending moment-curvature and applied torque-twist for the two beams.In general, good agreement can be found with the experiment.e moment-curvature response is better captured for the over-reinforced beam (TBO3) since it essentially behaves elastically in bending, whereas some differences can be for the under-reinforced beam (TBU3).
Results for the hoop strains at middepth are shown in Figure 15 for both beams.ese are somewhat overestimated in the postcracking range, especially for beam TBO3, due to overestimation of the twist.In both beams, yielding of hoops occurred in both y and z reinforcements at several crosssectional fibers.
From inspection of local strains in the longitudinal reinforcement, it was found that this did not yield for beam  Advances in Civil Engineering TBO3 but did yield for beam TBU3, with the values of 1.7‰ for TBO3 and 2.5‰ for TBU3 at the bottom ange.Concrete crushing occurred for the beam TBO3, reaching principal compressive strains of about 2.5‰ and corresponding stresses of 18 MPa, slightly less than f ′ c due to concrete softening at orthogonal tensile strains.For TBU3, principal compressive strains remained below 2‰.

Bending, Shear, and Torsion.
e test of the oor subassembly by Collins and Lampert [1] presented in Figure 1 is considered here for veri cation.is test setup produces torsion, bending, and shear in the spandrel beam, and bending and shear in the oor beam.e deection of the oor beam depends on the postcracked torsional and bending sti ness of the spandrel beam.Both Table 2: Material properties of the beam in bending and torsion [30].Advances in Civil Engineering beams were reinforced with longitudinal and transverse reinforcements, with the material properties summarized in Table 3. e T specimen was modeled with 4 elements, with 5 IPs each (Figure 16).A mesh of 10 × 10 fibers was used to discretize both cross sections.Longitudinal reinforcement 10 Advances in Civil Engineering was smeared in the top three layers and in the four bottom layers for the floor beam.For the spandrel beam, it was smeared in the top and bottom three layers.Transverse reinforcement was smeared over the entire section domain for both beams.e torsional rotation was fixed at nodes 1 and 3, being the rest of rotations free.At node 5, the beam was simply supported.A vertical concentrated load was applied at node 4 under displacement control with 2 mm displacement increments.Figure 17 shows the results for the applied loaddeflection at node 4 and for the load-twist response of the spandrel beam.Furthermore, the results from 3D FEA presented by Valipour and Foster [2] are shown.e authors used brick elements with a 3D nonlinear cementitious material model, which is a fracture-plastic constitutive law for cracked concrete, and 1D elements for the reinforcement with a linear-hardening stress-strain relationship.Good correlation can be observed with the experiments in both cases.At 80 mm vertical deflection, the twist predicted by the model is 0.010 rad/m, whereas in the experiment, FEA is about 0.014 rad/m.In the frame element analysis, yielding of longitudinal top and bottom reinforcement occurred in the floor beam and only of the transverse reinforcement in the spandrel beam.In the experiment, yielding of the bottom Floor beam Spandrel beam ρ x = 2.7% ρ y = 0.13% ρ z = 0.08% ρ x = 2.4% ρ y = 0.13% ρ z = 0.08% ρ x = 0% ρ y = 0.13% ρ z = 0.08% ρ x = 1.8% ρ y = 0.24% ρ z = 0.24% ρ x = 0.7% ρ y = 0.24% ρ z = 0.24% ρ x = 0% ρ y = 0.24% ρ z = 0.24% Advances in Civil Engineering longitudinal reinforcement also occurred in the spandrel beam.

Conclusions
Results from a 3D fiber-based frame element with multiaxial stress interaction have been presented.e proposed element is intended for nonlinear analysis of 3D frame structures such as those in the field of earthquake engineering, when consideration of multiaxial effects due to the presence of disturbed regions is important.Established frame analysis theories only consider the interaction between bending and normal forces with some "adjustments" for flexure-shear interaction.e proposed frame element formulation considers full coupling between bending, shear, and torsion in a consistent manner based on the fibersection discretization approach and 3D constitutive models.Despite the conceptual simplicity of the implemented constitutive models and beam kinematics, the element accuracy was comparable to that from refined finite element analysis, at least for the investigated case study involving shear, bending, and torsion coupling.Further verification studies, for instance, regarding cyclic shear and torsion, are expected in the future depending on the availability of experimental data.Potential enhancement of the constitutive models will be the object of future research as well.e possibility of replacing traditional flexural frame elements, based on the Euler-Bernoulli beam assumption, for global nonlinear analysis of framed structures should be further addressed.An important aspect, not considered in the present work, is the computational time, which obviously lies in between flexural frame elements and solid finite elements.

Figure 3 :Figure 4 :
Figure 3: Calculated warping functions for rectangular and hollow-core section.

Figure 5 :
Figure 5: Cyclic response of concrete (a) and cracked concrete element for tensile model (b).

Figure 6 :
Figure 6: Frame element degrees of freedom in the global (a) and basic (b) configurations.

Figure 8 :
Figure 8: Imposed top displacement on the columns.

Figure 9 :Figure 10 :
Figure 9: Numerical model of the Bousias column and section discretization.

Figure 16 :Figure 17 :
Figure 16: Numerical model and section discretization of the floor subassembly.
, where ε o is the average axial strain, χ y and χ z the bending curvatures, c oy and c oz the section shear deformations, and α is the derivative of the torsional rotation θ x .e above