Nonlinear Earthquake Analysis of Reinforced Concrete Frames with Fiber and Bernoulli-Euler Beam-Column Element

A beam-column element based on the Euler-Bernoulli beam theory is researched for nonlinear dynamic analysis of reinforced concrete (RC) structural element. Stiffness matrix of this element is obtained by using rigidity method. A solution technique that included nonlinear dynamic substructure procedure is developed for dynamic analyses of RC frames. A predicted-corrected form of the Bossak-α method is applied for dynamic integration scheme. A comparison of experimental data of a RC column element with numerical results, obtained from proposed solution technique, is studied for verification the numerical solutions. Furthermore, nonlinear cyclic analysis results of a portal reinforced concrete frame are achieved for comparing the proposed solution technique with Fibre element, based on flexibility method. However, seismic damage analyses of an 8-story RC frame structure with soft-story are investigated for cases of lumped/distributed mass and load. Damage region, propagation, and intensities according to both approaches are researched.


Introduction
The realistic modeling of the nonlinear static or dynamic behavior of RC structures is a more sophisticated problem due to the inelastic behavior of concrete, plasticity of reinforcement, interface debonding of these materials, and so on. In the last thirty years, researchers have made more effort for the numerical modeling of RC structures and most of the state-of-the-art on this problem deals with two main approaches: lumped plasticity modeling [1][2][3] and distributed-inelasticity modeling (i.e., the so-called fibre beam-column elements, FBCE) [4][5][6][7][8]. In the first approach, nonlinear springs based on moment-rotation and force-displacement curves are used and nonlinear volume is assumed to be lumped in the specific location of the element. Mohr et al. [1] were used a series of polynomial shape functions for strain distributions of the vertical and shear on cross-section. Li et al. [2] was proposed a simplified lumped hinge for determination of failure mode of statically indeterminate structure. Reshotkina and Lau [3] was used a concentrated plasticity approach for axial-flexure-shear interactions on the inelastic behavior of reinforced concrete members under the seismic loads. In the second approach, the nonlinear behaviors of concrete and reinforcement materials are separately calculated. This method divided into two as approaches was based on flexibility and rigidity. Taucer et al. [9] were developed a flexibility method for nonlinear dynamic analysis of reinforced concrete element. This method based on Bernoulli-Euler hypothesis included biaxial bending and axial force conditions. Nonlinear behavior of concrete and steel were computed by using uniaxial stressstrain relationships. Force-displacement interpolation functions are used for the obtaining of element flexibility matrix. Element stiffness matrix is achieved by inverse of flexibility matrix. Furthermore, Ceresa et al. [10] were developed a beam-column element based on flexibility matrix by using to the Timoshenko beam theory under cyclic loading. Damages in the element are taken into account with coaxial rotating crack model. Lu et al. [4] were used multilayered shell element and fibre beam-column element based on flexibility method for modeling of the different frame and frame with shear wall. However, an important research report about methods based on rigidity was performed by Kawano et al. [6]. In the report, nonlinear analyses of reinforced concrete space frames with multistory and multibay are obtained under the earthquake 2 The Scientific World Journal loads. Cubic Hermitian polynomials for transversal displacements and the linear shape functions for axial displacements are used for obtaining of stiffness matrix. However, the Gauss-Lobatto integration rule was used for obtaining element mass and stiffness matrices. Element damage location was simply presented and element was not divided into subelement called as "segment. " Some calculation problems in heavily damage location was shown and they are reported to arise from flexural deformation. Furthermore, strain rate effects in response to reinforced concrete frames was investigated by Iribarren [7]. Used formulations are based on rigidity method and a strain rate dependent material formulation is developed for both the concrete and steel constitutive response. Brum [8] was obtained nonlinear dynamic analyses of frame and masonry structures by using a method which was based on stiffness matrix. In the research, a uniaxial constitutive model for concrete and masonry is proposed for compression and tension regions of the materials under the cyclic loadings.
In this study, a beam-column element based on the Euler-Bernoulli beam theory is presented for the fiber RC element. Element stiffness matrices are obtained by using rigidity method. The beam or column element is divided into subelements called "segment" [11]. Furthermore, the internal freedoms of this segment are dynamically condensed to external freedoms at the end of the element. Thus, nonlinear dynamic analysis of high RC building can be obtained within a short time. This procedure requires that nested loops of the element and structure are obtained at the same time (the so-called nonlinear dynamic substructure). However, this condensation procedure is not used in the modeling of the fibre element approach (FEA) of the RC element. In addition, uniform or trapezoidal loads of the segment are assumed to be zero or condensed to the external freedoms at the end of the element in the FEA [11]. This case is not preferred due to the required redistribution of loading. However, in this study, external loading of the segment is taken into account by considering the damage occurred in the element.
The present research is organized as follows: (i) presenting of fiber Bernoulli-Euler beam-column element based on rigidity method, (ii) developing of nonlinear dynamic substructure technique of RC frames, (iii) verification of the constitutive model with respect to experimental result of a column structural element and dynamic analysis results of fibre RC element of a portal frame, (iv) obtaining of seismic damage analyses of an 8-story RC frame with soft-story for distributed/lumped mass and load case, and (v) results.

Fiber and Bernoulli-Euler Approach (FBEA) for Reinforced Concrete Beam Column Element
In this section, theory of Fiber and Bernoulli-Euler element which is based on stiffness matrix is firstly presented for reinforced concrete section. Elements are divided into subelements, called segment, and they are assumed as substructures of element. Furthermore, cross-section of the each segment is also subdivided into a number of fibers/layers. In the next section, nonlinear dynamic substructures method is also mentioned to be applied to freedoms in element ends.

Obtaining Segment Stiffness and Mass
Matrices. The strain distribution on a cross-section of beam is assumed to be uniform due to axial forces and linear due to only bending according to the Bernoulli-Euler approach. Furthermore, if a plane section before bending is plane after bending and if the strain is assumed to be small and shear stresses are omitted, the strain of a point on a cross-section in the axial direction can be written as where and V are the displacement of the axial and vertical directions of element, respectively ( Figure 1). However, if cross-section of element is divided into fibers/layers, this equation for each fiber/layer in local axis can be rewritten as, where , is incremental axial strain of a fiber/layer on a segment in local axis direction. Therefore, if the cubic Hermitian polynomials for transversal displacements and the linear shape functions for axial displacements are used [12], (2) can be written as where { } is displacement vector, including displacement and rotation of a segment. [ ] is strain-displacement matrix. Thus, incremental stress in each fiber/layer can obtained as where , , determined by using uniaxial stress-strain relationship of using materials, is tangent elasticity modulus of each fiber. Therefore, total strain energy of a segment subelement is obtained as ] ] ] The Scientific World Journal where Seg and Seg are expressed to be length and crosssection area of segment, respectively.
is also area of a fiber/layer. Thus, element stiffness matrix can be obtained by using minimum potential energy principle. Stiffness matrix of a segment can be written as The segment stiffness matrix is obtained by using areas, coordinates, and tangent elasticity modulus of fiber/layer for reinforced concrete sections. In this study, linear superposition rule is used for different material properties of the fiber/layer in the section ( Figure 1). Furthermore, a fiber concrete or reinforced bar on the cross-section may be cracked or damaged due to external loading. For this reason, a relationship between the damaged and undamaged cases must be obtained for the solutions. Thus, damage in each concrete/reinforcement fiber can be written as [13] where , and , are undamaged and damaged/tangent elasticity module of the th fiber, respectively. , damage intensities are obtained separately under tensile and compressive stress for each fiber. However, total kinetic energy of particle velocities on the cross-section of an element throughout its neutral axis can be written as where is mass density and V is particle velocity. If fiber/layer elements are used for (8), the mass matrices of an element on the local axis are obtained as where [ ( )] is the element shape functions matrix on the local axis and is also mass density of the th fiber. Stiffness and mass matrices of a segment, obtained by local axis, are transformed to global axes by using transformation matrix.

Obtaining with Nonlinear Dynamic Substructures Technique of Element Stiffness and Mass Matrices.
In FBEA, cubic Hermitian shape functions for rotation and shear strain and linear shape functions for the axial deformation are used, respectively. The Gauss-Lobatto integration rule is generally used for obtaining element mass and stiffness matrices along with local axis [6][7][8]. However, substructure procedures are generally preferred for solutions of high buildings. However, this technique in the previous study is not used due to numerical difficultly. Element stiffness/flexibility matrix is obtained by using cross-section properties on the integration points. Therefore, effects of shape functions on the solution are very important. In this study, an element is divided into subelements which are called "segment" to remove disadvantage effects of the shape functions. The freedoms of segment are condensed to the end freedoms at the ends of the element. Thus, if freedoms at ends and internal regions of element are called as external and internal freedoms, respectively ( Figure 2), these external and internal freedoms for stiffness matrices can be written as, where [ ] and [ ] are stiffness matrices which include external and internal freedoms, respectively [11]. { } and { } are the displacement vectors referred by these freedom; { } and { } are also the external loads of these freedoms. When applying the substructure procedure to (10a) and (10b), the external load and stiffness matrices of a frame element can be rewritten as However, if the Rayleigh method may be used to obtain element damping matrices, this matrix is defined as where dam and dam are Rayleigh damping coefficients with respect to mass and stiffness matrices, respectively [14].
However, if substructure procedure in (11a), (11b), and (11c) is applied to mass and damping matrices [15], mass and damping matrices can be obtained for the external freedom. Global stiffness, damping, and mass matrices may be achieved by using the stiffness and mass matrices belonging to external freedom. Thus, global stiffness, mass matrices, and external load vector can be calculated as Therefore, the dynamic equilibrium equations of the structure may be given as where subscripts gr and stat indicate that the quantity is related to ground acceleration and static loads.

Bossak-Form of the Equation of
Motion. The equation of motion for the RC frame is given in (15). In this study, the Bossak-integration method, presented by Wood et al. [16], is used for the solution of the equation in the time domain. Integration scheme of the method retains the Newmark method. Furthermore, the Bossak-integration method is required to be modified to (15) in the time domain as follows: where is the Bossak parameter, used for controlling the numerical dissipation. The Bossak parameter should be chosen as in (17) for unconditional stability and second-order accuracy: The Scientific World Journal In this study, is selected as −0.10. To solve the nonlinear dynamic equation of motion for the RC frame, the Newton-Raphson method is used in conjunction with the predictorcorrector technique. The Bossak-time integration algorithm is given by Wood et al. [16]. Predicted displacement and velocity vectors for the time step + Δ are obtained by using displacement and velocity vectors at the time step which is known. Thus, they can be calculated as where and are Newmark's coefficients. The displacement and velocity vectors [17] of the RC frame can be written in terms of the predicted vectors shown in (18a) and (18b). The vectors may be defined as These relations can be substituted into (16) and a time marching algorithm can be applied to this equation as given in the Appendix.  [18]. The experimental setup is shown in Figure 3. Exponential decreasing functions for the softening region of tensile and compressive strengths of the concrete for the constitutive model of FBEA are used. These functions are shown in Figure 4(a). The bilinear kinematic hardening rule is used for nonlinear behavior of the steel for the two approaches (Figure 4(b)). Static loads are converted to masses which are condensed to the element ends for all solutions and displacement values obtained due to the loads being considered as the initial condition. Acceleration data, a sinus wave shape, and time varying are used for all dynamic solutions. This dynamic load is applied to the horizontal direction at 1280 mm height of RC column. Tangent stiffness matrix is used for the solution and damping matrix is also assumed to be proportional to stiffness matrix. Load-displacement curves of experimental and numerical results are given in Figure 5. Maximum and minimum values of cyclic displacement obtained from numerical results are approximately between −0.055 and 0.055 m. Numerical response of RC element which is obtained with FBEA is shown as similar to envelope curve of experimental results. All values of the horizontal displacements are on the envelope curve of the experimental result. It is said that this solution technique and material models of concrete and steel can be used for the solution of the RC structural element under the dynamic loading.

Nonlinear Dynamic Analyses of a Portal RC Frame
Structure. In this section, nonlinear dynamic analyses of a portal RC frame are obtained for the comparing of proposed solution technique with model based on flexibility. Seismo-Structure program [19] is used for the flexibility based solutions (with fibre element method). Finite element meshes are shown in Figure 6 for both approaches. Finite element mesh of proposed method is divided into segment for comparing with analysis results of Seismo-Structure Program. ACI 318-02 [20] code is used for the material properties of concrete  and steel. Cross-section and material properties of the beam and column of selected portal RC frame are given in Table 1.
The tensile softening region in the fibre element method is not taken into account while the compressive behavior of the concrete is used. In Seismo-Structure program, if the tensile stresses obtained from the solutions reach the tensile strength of the concrete, tensile strength suddenly drops to zero. In this case, a sudden collapse of the structure occurs. In this study, exponential decreasing functions in the softening region of tensile and compressive strengths of the concrete for constitutive model of FBEA are used. These functions are shown in Figure 4(a). The bilinear kinematic hardening rule is used for nonlinear behavior of the steel for the two approaches (Figure 4(b)). Static loads are converted to masses which are condensed to the element ends for all solutions and displacement values obtained due to the loads being considered as the initial condition. Acceleration data and a sinus wave shape are used for all dynamic solutions (Figure 7). This dynamic load is applied to the horizontal direction.
Displacement time history graphs of node 3 obtained from FBEA and FEA are shown in Figure 8 obtained from the FEA are bigger than those obtained from the FBEA. This case arises from not taking into account the concrete tensile softening region for the fibre element method. However, failure of structure is not seen for two approaches in all times. Accumulated tensile damage cases obtained for both approaches are given in Figure 9.

Seismic Damage Analyses of an 8-Story RC Frame Structure with Soft-Story.
In this numerical application, nonlinear dynamic analyses of an 8-story RC frame structure with softstory are investigated for cases of lumped/distributed mass and load. ACI 318-02 [20] code is used for the material properties of concrete and steel. Cross-section and material properties of the beam and column of selected RC frame are given in Table 2. Finite element mesh and gravity loading case are shown in Figure 10. Uniaxial stress-stain relationships of the concrete and steel are seen in Figures 4(a) and 4(b), respectively. Static loads and displacements are considered as initial conditions in all solutions. Spectrum acceleration The Scientific World Journal 9  curve, given by Z1 type soil in the Turkish Regulation Code on Building in Disaster Areas [20], is selected as the target spectrum curve. Synthetic earthquake acceleration data for maximum amplitude, 0.3 g, are produced. This synthetic acceleration-time graph is shown in Figure 11 and it is affected on the horizontal direction of the RC frame structure. Tangent stiffness matrix is used for the all solution and damping matrix is also assumed to be proportional to stiffness matrix. Displacement time history graphs of node 9 obtained from nonlinear dynamic analyses results for lumped/ distributed mass and load cases are shown in Figure 12. Displacement amplitude values are approximately similar until time of 1.08 sec and, after this time, displacement amplitude values obtained from the distributed mass and load case are bigger than those obtained from the lumped case. Absolute maximum horizontal displacements for the lumped and distributed cases are obtained as 22.6 and 41.4 mm, respectively. Thus, absolute maximum displacement value according to distributed case is approximately 83% of ratio bigger than that of lumped case. This result arises from the redistribution of loading in the internal regions of the elements for the distributed approach. However, the redistribution of loading in these internal regions is only obtained at the end regions of elements for the lumped approach. The redistribution procedure for loading in these internal regions requires more iterative steps for a distributed approach. However, numerical dissipation is not shown for both solutions and Newton-Raphson procedures for all element solutions which are converged for the two approaches. The Bossak-dynamic integration algorithm is applied successfully to the nonlinear dynamic solutions.
Accumulated tensile damage regions obtained for both mass/load case are given in Figures 13 and 14. Damage zones in the lumped case are obtained at the end of the beam of all floors and at the whole of the beam cross-section for the time = 1.00 sec. Damage zones occurred at upper parts of the beams at the beam-column join region and at lower parts of the mid region of the beams until this time for the distributed case. Furthermore, damage zones seen at lower parts of first-floor columns occurred for both approaches, but intensities and regions of these damages according to distributed approach are bigger than those of lumped mass/load case. The damage intensities are some increased and damage propagations are remained at the same regions for both approaches until = 3.00 sec; after this time, damage intensities according to two approaches are increased at between = 3.00 and 8.00 sec, and no change for the damage zones is obtained. However, increases in the damage intensities achieved a minimal level between = 8.00 and 10.00 sec. It is said that the damage zones obtained according to the distributed/lumped case had important differences.

Conclusions
In this study, a beam-column element based on the Euler-Bernoulli beam theory is researched for nonlinear dynamic analysis of reinforced concrete (RC) structural element. Stiffness matrix of this element is obtained by using rigidity method. A solution technique that included nonlinear dynamic substructure procedure is developed for dynamic analyses of RC frames. A predicted-corrected form of the Bossak-method is applied to dynamic integration scheme. A comparison of experimental data of a RC column element with numerical results, obtained from proposed solution technique, is studied for verification of the numerical solutions. Furthermore, nonlinear cyclic analysis results of a portal reinforced concrete frame are achieved for comparing of the proposed solution technique with fibre element, based on flexibility method. However, seismic damage analyses of an 8story RC frame structure with soft-story are investigated for cases of lumped/distributed mass and load. Damage region, propagation, and intensities according to both approaches are researched. Results obtained from this study are presented to be itemized as follow. (iv) When the results obtained from lumped and distributed approaches are investigated, absolute maximum displacement values for distributed approach are seen to be bigger than those for lumped approach. This case arises from the redistribution of loads in the internal regions of the elements for the distributed approach. However, the redistribution of loads in these internal regions is not used for the lumped approach. The redistribution procedure for loading in these internal regions requires more iterative steps.
(v) Obtained damage zones for the lumped mass/load case are seen at the end of the beam of all floors and at the whole of the beam cross-section, but obtained damage regions for the distributed mass/load case showed up at upper parts of the beams at the beamcolumn join region and at lower parts of the mid region of the beams. Furthermore, additional damage zones at lower parts of first-floor columns occurred for both approaches, but intensities and regions of these damages according to distributed approach are bigger than those of lumped approach.
(vi) However, numerical dissipation is not shown for all solutions and Newton-Raphson procedures for all element solutions are converged for the two approaches. The Bossak-time integration algorithm and proposed solution technique are successfully applied to the nonlinear dynamic solutions.

Time Marching Algorithm
The algorithm applied to (16) is as follows.
(9) Compute the effective dynamic stiffness matrix and the vector of unbalanced forces, where [ ] +Δ and {Δ ,res } +1 +Δ are the tangent stiffness matrix and the incremental restoring force vector of external joints on the global axis of the structure, respectively, which are assembled from the element contributions. (11) Update displacement, velocity, and acceleration vectors, Strain of a point on a cross-section in the axial direction of element , : Incremental axial strain of the th fiber/layer in local axis direction of a segment : Second Newmark's coefficient : Vertical coordinates of the th fiber/layer in local axis direction of a segment : Mass density Π: Total strain energy of a segment.