Abstract
A three-dimensional fiber-based frame element accounting for multiaxial stress conditions in reinforced concrete structures is presented. The 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-of-plane displacements perpendicular to the cross section due to warping effects. The formulation was implemented into a force-based 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.
1. Introduction
The 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 one-dimensional 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]. This creates multidirectional shear flow in the spandrel beam coupled with normal stresses arising from bending. The 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–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. The 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. The constitutive model was based on the MCFT. The 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. The beam kinematics was based on the Timoshenko beam theory in conjunction with the Coulomb torsion, that is, no warping was considered.
The abovementioned work was based on the classical Timoshenko beam theory with either uniform or parabolic shear profiles. This 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–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. The 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. The Saint-Venant theory is used for modeling torsion, which is suitable for noncircular solid cross sections with warping effects. The 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.
2. 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. The kinematics of any point of the beam is completely described by three independent displacements and rotations of the beam axis as follows (Figure 2):
The total displacement at any given point of the section is related to the beam axis displacements aswhere , , and are the axial and transverse displacements of a generic point of the section with coordinates (), is the warping function, and is the warping parameter equal to ∂θx/∂x according to the Saint-Venant theory. Assuming small displacements, axial and shear strains at a given section fiber are given aswhere is the average axial strain, and the bending curvatures, and the section shear deformations, and is the derivative of the torsional rotation θx. The above relationships can be expressed in matrix as ε = Bεe, where Bε is the strain-displacement compatibility matrix, ε = (εx, γxy, γxz), and e = (εo, γoy, γoz).
The section forces consistent with the section strains are obtained from integration of the axial and shear stress distributions over the cross section as follows:
The 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 is the condensed (3 × 3) tangent material stiffness matrix. Note that and are different; hence, will not be symmetric in the general case, with the exception of the circular cross section for which the warping function is zero. The matrix is obtained from condensation of the (6 × 6) material stiffness matrix assuming that the total stresses , , and are zero in the beam element as follows:where and are the stresses in concrete and steel, respectively, and (i = y, z) is the reinforcement ratio at each individual section fiber. This process is performed iteratively at each section fiber until convergence, which requires implementing an additional loop within the fiber state-determination procedure.
The 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 represent the cross-sectional domain and its boundary, respectively, and and 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.
(a)
(b)
3. 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]. The model is formulated in the material directions (Figure 4), which are coincident with the crack directions determined at first cracking aswhere is the plastic compressive strain (positive) and is the principal strain (negative for tension). can be obtained upon solving the eigenvalue problem for the strain tensor in the xyz reference system, . 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, , are obtained from aswhere is the standard (6 × 6) transformation matrix. The corresponding concrete stresses, , 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 aswhere and are the total tensile strains in the orthogonal directions and .
(a)
(b)
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 different 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)). Thus, average concrete stress in the ith material direction can be obtained aswhere , , and are concrete stresses in the (b), (sl), and (cr) regions, respectively. , , and are the average spacing, unbonded length parameter, and tensile strain in the ith direction, respectively. This approach allows full characterization of the tension stiffening, crack-closing, and crack-opening regimes based on principles of mechanics. Further details can be found elsewhere [23, 25].
Shear stresses, , arising from deviation between principal directions and crack directions are related to the shear strains, , aswhere , and , 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 is the maximum aggregate size and is the average crack width (mm). Assuming smeared cracks, the average crack width is given aswhere is the maximum (in absolute value) tensile strain normal to the crack. The average crack spacing, , 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 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 modified isotropic hardening [27]. For the ith reinforcement group, which orientation is defined by the unitary vector , the total axial strain can be obtained upon transforming to the corresponding direction. Then, steel stresses are obtained using the abovementioned model and transformed back to the xyz reference system aswhere is the uniaxial steel stress, is the reinforcement ratio for the ith reinforcement group, and is the corresponding transformation matrix between reference systems. stresses are added to the concrete stresses, , in order to obtain the total stresses, , which are used for the calculation of the vector of internal element forces. Similarly, the ith reinforcement contribution to the material stiffness matrix, , is added to that provided by concrete, , such that the total material stiffness matrix is given aswherewhere and are the tangent moduli for concrete and steel and is the shear stiffness.
4. Frame Element Implementation
The 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. This is performed in the basic element configuration upon removal of the rigid body modes as follows (Figure 6):where is the increment of element displacements in the basic configuration, is the increment of global displacements, at a given load step, and is the rigid body transformation matrix given as
(a)
(b)
The element forces, , can be obtained from the element flexibility matrix, , which is known from the previous converged load step as
Equation (25) is only an approximation since 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 aswhere is the force interpolation matrix at the kth integration point with coordinate , which is given as
An estimate of the section deformations at the kth integration point, , can be made based on the section flexibility matrix, , 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, , and an updated section flexibility matrix. The unbalanced section forces, , given by the difference between and , are used to update the section deformations as , which are further used to correct the element forces by the following quantity:where are the integration weights at the kth integration point and is the updated element flexibility matrix given as
The procedure is repeated with the new element forces until satisfying element convergence, that is, when < 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.
5. Verification Examples
5.1. Biaxial Bending Test
Bousias et al. [29] tested flexure-dominated 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. The shear span-to-depth ratio was 6. The 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.
(a)
(b)
The columns were modeled with one single FB element with 5 IPs. The cross section was discretized into 10 × 10 fibers, with the longitudinal and transverse reinforcement smeared as shown in Figure 9. The 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. The main effect of the transverse reinforcement is confining of concrete, which was accounted for with an estimated confinement factor of 1.1.
The hysteretic lateral force-displacement response is shown in Figure 10. The 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. This 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.
(a)
(b)
(a)
(b)
5.2. Torsion and Bending
The series of RC beams tested by Onsongo [30] in combined torsion and flexure are analyzed herein. The 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. The test setup consisted of adjustable steel level arms at the ends of the specimen, where 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. The 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. The 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. The 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. The moment-curvature response is better captured for the over-reinforced beam (TBO3) since it essentially behaves elastically in bending, whereas some differences can be observed for the under-reinforced beam (TBU3).
(a)
(b)
(a)
(b)
Results for the hoop strains at middepth are shown in Figure 15 for both beams. These 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 cross-sectional fibers.
(a)
(b)
From inspection of local strains in the longitudinal reinforcement, it was found that this did not yield for beam TBO3 but did yield for beam TBU3, with the values of 1.7‰ for TBO3 and 2.5‰ for TBU3 at the bottom flange. Concrete crushing occurred for the beam TBO3, reaching principal compressive strains of about 2.5‰ and corresponding stresses of 18 MPa, slightly less than due to concrete softening at orthogonal tensile strains. For TBU3, principal compressive strains remained below 2‰.
5.3. Bending, Shear, and Torsion
The test of the floor subassembly by Collins and Lampert [1] presented in Figure 1 is considered here for verification. This test setup produces torsion, bending, and shear in the spandrel beam, and bending and shear in the floor beam. The deflection of the floor beam depends on the postcracked torsional and bending stiffness of the spandrel beam. Both beams were reinforced with longitudinal and transverse reinforcements, with the material properties summarized in Table 3.
The 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 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. The 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 load-deflection 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. The 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 longitudinal reinforcement also occurred in the spandrel beam.
(a)
(b)
6. Conclusions
Results from a 3D fiber-based frame element with multiaxial stress interaction have been presented. The 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. The proposed frame element formulation considers full coupling between bending, shear, and torsion in a consistent manner based on the fiber-section 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. The 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.
Data Availability
The authors declare that the data that support the findings of this study are available from the corresponding author on request. Furthermore, the data related to the verification examples are available within the article and in the papers cited in the references.
Conflicts of Interest
The authors declare that there are not conflicts of interest regarding the publication of this paper.
Acknowledgments
The Ph.D. scholarship support provided by the University School for Advanced Studies IUSS Pavia, Pavia, is greatly acknowledged.