An energy-based computational scheme for the analysis of reinforced concrete structures with geometric and material nonlinearities

ABSTRACT This paper proposes a new energy-based numerical technique for the nonlinear analysis of reinforced concrete beam elements. The mathematical concept and equations are presented in detail for a general nonlinear quasi-static problem. This method can create a physical sense of the types of energy in a system for an analyst to easily consider any material type and geometric nonlinearity effects from the analysis of reinforced concrete members. The governing equations of this computational scheme have a quadratic form, which leads to proving a complementary relationship for a stability check. In other words, a distinguishing feature of the proposed method is paving the path for selecting the optimum size of displacement increments during a nonlinear numerical analysis. In the proposed method, the steps have been made easier for analysts. Unlike several existing methods, the formulation has also considered the shear and normal deformations, making the analysis results closer to experimental realities. Moreover, the results obtained from solving, evaluating, and comparing several examples demonstrate the acceptable accuracy and convergence of this energy-based approach for reinforced concrete beams. In addition, the quadratic nature of its formulation provides an analyst with information about the accuracy and stability of the obtained structural responses.


Introduction
Today, reinforced concrete (RC) structures constitute a significant part of engineering structures due to their durability and resistance against corrosion and fire, despite their heavier weight than steel structures. In these systems, the significant compressive strength of the concrete, along with the high tensile strength of the steel bars, provide a resistant composite material. While steel and concrete are well suited to one another, the presence of two distinct materials complicates finite element (FE) analysis of such structures, especially in cases of near-failure exhibiting extreme nonlinear behavior. Many researchers have worked on modeling and analyzing the nonlinear behavior of RC structures. Jakobsen and Grenacher were among the first researchers who investigated slender RC frames (Aas-Jakobsen and Grenacher 1974). Kayal et al. conducted a study into the nonlinear interaction of RC frame-wall systems (Kayal 1986). Gunnin et al. proposed a computational method involving material nonlinearity, axial shortening, and P-Δ effects (Gunnin, Rad, and Furlong 1977). Subsequently, with advances in computer science, modeling and structural analysis methods also made significant progress, and in-depth research was pursued specifically on the nonlinear FE analysis of RC structures. For example, Marni and Spacone (2006) focused on the shear effects on the analysis of RC elements, which led to presenting a force-based 2D element based on the Timoshenko beam theory. Wackerfuß propped an FE formulation to analyze localized failure in beam structures (Wackerfuß 2008).
The fracture process in this study was based upon the cohesive crack concept using the traction-separation law. In another study, Lezgy-Nazargah et al. provided a refined global-local FE model for bending analysis of laminated composite rectangular beams (Lezgy-Nazargah, Beheshti-Aval, and Shariyat 2011). The presented element was free of shear locking and satisfied the continuity conditions of the displacements and transverse shear stresses at the interfaces and the non-homogenous shear tractions on the upper and lower surfaces of the beam. Overall, comparing the results obtained by the study with previous studies and the 3D theory of elasticity has shown that this formulation is accurate for thin and relatively thick beams.
Regarding the numerical complications in the FE analysis of brittle systems, Slobbe et al. introduced an alternative FE method and compared the obtained results with nonlinear analyses using an incrementaliterative solution procedure (Slobbe, Hendriks, and Rots 2012). The influence of the adopted amount of material fracture energy, the number of damage increments in the discretized stress-strain curve, the shear retention after cracking, mesh refinement and alignment, and the possibility of non-symmetric failure modes were also determined by performing a series of sensitivity analyses. In another study, Stramandinoli and La Rovere (2012) developed an FE model based on the Timoshenko Beam Theory (TBT), together with the Modified Compression Field Theory (MCFT). The numerical results showed good agreement with experimental data, and it was claimed that the presented 1D model was more stable and economical than 2D models. Bui et al. also proposed an enriched Timoshenko beam finite element method for modeling the bending and shear failure modes in RC frames (Bui et al. 2014), in which the FE enhancement was based on the Strong Discontinuity Approach (SDA). The numerical simulations carried out in this study, and the comparison of the results with those of the experiments generally indicated the appropriate performance of the proposed methodology to compute the ultimate load of RC frames. Šćulac, Jelenić, and Škec (2014) introduced a novel FE technique with embedded transversal cracking to analyze RC beams. The main focus of their research was on the definition and FE execution of kinematics of a layered beam setup.
The nonlinear finite element (NFE) modeling of Reinforced Concrete Haunched Beams (RCHBs) designed to develop a shear failure was another relevant subject investigated by Godínez-Domínguez, Tena-Colunga, and Juárez-Luna (2015). They assessed the ability and limitation of their NFE model to predict the shear failure considering eight simply supported RCHBs. There were two types of modeling in this study. In the first one, the effects of longitudinal bars and stirrups were indirectly simulated using SAP2000 software, while in the second one, ANSYS software was employed to model the longitudinal and transverse bars used in experimental models. Hawileh (2015) also developed a 3D FE model including material nonlinearity and bond-slip behavior between the reinforcement and adjacent concrete in the presence of aramid-fiber-reinforced polymer (AFRP) rebars (Hawileh 2015). Experimental results validated the presented model, and the accurate prediction of the loadcarrying capacity of these structural components was reported. In a similar work, Nistico et al. numerically investigated the behavior of RC beams strengthened against shear using carbon-fiber-reinforced polymer (CFRP) (Nisticò, Ožbolt, and Polimanti 2016). A 3D FE analysis was performed using the MASA code based on the microplane modeling. They found that it was consistent with the experimental data in predicting the damage propagation in CFRP. Dere and Koroglu (2017) studied the NFE modeling of plain RC structures through concrete damage plasticity models available within the commercial ABAQUS software package. Based on the Refined High-order Global-local Beam (RHGB) theory, Lezgy-Nazargah (2017) developed an effective and time-consuming 1D FE model to assess the progressive failure of laminated composite beams. They adopted an iterative method to predict the first-ply failure load. Afterward, the material properties of the failed elements were modified using a stiffness degradation factor. Finally, the applied load was increased step-by-step for each load increment until the ultimate strength of the laminate was achieved. In another study, a materially nonlinear model was formulated for RC beams.
Contrary to the kinematic assumptions of 1D beam elements generally developed based on the Euler-Bernoulli or Timoshenko theory, a Layered Global-Local (LGL) theory was used to consider the effects of shear deformation and normal flexibility. The nonlinear behavior of cracked concrete and reinforcing bars was also simulated using fixed smeared crack and elastoplastic models. In summary, a comparison of the results revealed that the proposed method yielded accurate results at a low computational cost (Lezgy-Nazargah 2018).
Regarding the critical role of beam-column connections' behavior in RC frames' resistance, ductility, and stability, Najafgholipour et al. developed a 3D FE modeling for these components governing joint shear failure mode (Najafgholipour et al. 2017). The joint shear capacity, deformations, and cracking pattern were highlighted in their investigation. The comparison between experimental and numerical results demonstrates that their proposed FE model could properly simulate the behavior of the joint while capturing the shear failure. Similarly, Pan et al. studied the effect of local joint deformations on the global frame response (Pan, Guner, and Vecchio 2017). They modified a previously developed global frame analysis process. They found that the predicted to observed peak loads ratio had mean values of 1.25 and 1.05 before and after the modification. Kakavand et al. (2018) proposed a 3D continuum FE model to evaluate RC frames' nonlinear behavior and failure modes in pushover analysis. To this end, four RC frames were analyzed through the nonlinear static procedure. The comparison of numerical and experimental results indicated the capability of the refined FE models to capture the lateral loadcarrying capacity and the location and evolution of concrete damage.
Mbewe and van Zijl (2018) performed a simplified nonlinear analysis of RC frames with masonry infill subjected to seismic loading. They proposed a novel truss analogy analysis method and an analytical solution procedure validated by experimental data. Li et al. also proposed a modeling method for collapse analysis of RC frames with infill walls (Li et al. 2019). Through the work, a three-strut model was proposed to model the overall response of infill walls. Finally, it was found that the presented model was able to predict the maximum resistance force accurately. Vafaei, Alih, and Fallah (2019) investigated the accuracy of the lumped plasticity model for predicting the nonlinear behavior of RC frames under gradually increasing vertical loads. For this purpose, two full-scale RC frames with different shear spans were modeled in SAP2000, and the results were compared to those obtained from the experiments. The effects of different plastic hinge lengths, initial effective stiffnesses, and locations of plastic hinges were examined entirely. Eventually, Mohemmi, Broujerdian, and Rajaeian (2020) proposed an equivalent method for rebar slip simulation of RC frames based on pull-out test results. The FE results obtained from the ABAQUS environment suggested a steel strength reduction factor equal to 0.6 for practical analyses.
As it can be concluded from the literature review, various studies have been carried out on the nonlinear behavior of reinforced concrete structures. However, more research is still required regarding the complexity of modeling such components, especially in simulating their nonlinear behavior and determining their structural response near failure.To this end, this study proposes a novel methodology for nonlinear FE analysis of RC beam elements based on the energy balance equation instead of using the force equilibrium equations. The main idea of this research is based on the modified energy method (MEM) proposed by Abad et al. in nonlinear structural dynamics (Abad et al., 2017). In fact, this method is being combined with LGL theory to solve the set of nonlinear equations involved in the response analysis of RC beams. It is noteworthy that the numerical efficiency of this method is investigated to determine the postbuckling response of framed structures with bifurcation points studied by Razaghi, Asgari Marnani, and Rohanimanesh (2019). The efficiency of this method in the nonlinear FE analysis of RC beams is studied considering both nonlinear material and nonlinear geometrical behavior.
The structure of the research is as follows. Section 2 provides the nonlinear FE formulation of RC beams. Section 3 describes the presented methodology along with its computer implementation. Section 4 presents the results of some illustrative numerical examples and comparisons with experimental and other related works. Finally, Section 5 provides the major outcomes of the study based on the interpretation of the obtained data.

Finite element formulation
The following section first introduces the concepts of equivalent multilayered cross-sections and other modeling assumptions. The equilibrium equations in the analysis of the structural components are then expressed at the end of the section.

Equivalent multilayer cross-section for reinforced concrete members
As mentioned earlier, one of the most significant complexities in the analysis of RC members is the presence of two different materials (steel and concrete) simultaneously inside the cross-section, having entirely different mechanical behaviors against applied loads. On the other hand, the two materials generally interact with each other, which adds to the analysis's complexity.
The present study uses an equivalent multilayered section to simulate the response of such members by assuming a complete bonding between the surface of rebars and the adjacent concrete (i.e., no separation occurs between the two materials during the analysis). Figure 1a depicts a concrete rectangular section (width = b and height = h), where the coordinate system is assumed to be centered. In this section, two layers of longitudinal rebars are used, including compression rebars with an area of A′ s located at a distance of d′ from the upper edge and tensile rebars with an area of A s and an effective depth of d. According to Lezgy-Nazargah, Beheshti-Aval, and Shariyat (2011), this section can be considered an equivalent n-layered section, as shown in Figure 1c, by defining equivalent layers for steel rebars (the equivalent thickness is calculated by dividing the area of the rebars by the width of the section).

Displacement and strain fields
Herein, a two-dimensional displacement field in the form of U = [u, w] T is assumed for the FE formulation. In general, the Bernoulli or Timoshenko hypotheses can be used to express the components of transitional displacements u and w, respectively, along the x 1 -and x 3 -directions. However, the beam theory of Lezgy-Nazargah, Beheshti-Aval, and Shariyat (2011) has been used in this research, in which high-order displacement functions, considering shear deformations, are expressed as follows: According to Figure 2, in these expressions, u 0 (x 1 ) and w 0 (x 1 ) respectively stand for the mid-translational displacements in x 1 -and x 3 -directions and ψ(x 1 ) is the mid-plane rotation about x 2 -axis. The comma refers to the differentiation concerning global coordinates. Moreover, z (x 1 ) and x (x 1 ) denote the prescribed shear and normal tractions, in which the + and ̶ superscripts denote the upper and lower planes, respectively. � u 1 1 (x 1 ) defines an additional unknown quantity associated with the refinement per layer, and k is the layer number. Ψ, Ω, Λ, and ϒ are some functions defined in terms of shear modulus and layer coordinates (please refer to Lezgy-Nazargah (2018)).
Subsequently, the strain field components via assuming the plane-stress state can be expressed as ε = [ε 11 ,ε 13 ,ε 33 ] T , where the components of the Lagrangian strain tensor, ε, are defined as follows by employing indicial notation: Note that the last term on the right-hand side of Eq. (3) involves nonlinear geometric terms. Hence, the strain field can be stated as follows for the problems considered in this study.
Consequently, by combining Eq. (4-6), the displacement-strain relationship can be rewritten as a matrix: In which L is a 3 × 2 tensor that relates strains to displacements. As mentioned earlier, nonlinear terms exist in the present problem. Therefore, in practice, it would be appropriate to write this matrix as a sum of linear and nonlinear terms denoted by L l and L Nl , respectively.

Constitutive laws
The models used to define nonlinear behaviors of each material, i.e., steel and concrete, are discussed in this section, which involve the constitutive equations for longitudinal steel rebars, uncracked and cracked concrete, and transverse reinforcement.

Longitudinal steel rebars
The behavior of steel reinforcing rebars is usually assumed to be elastoplastic in both tension and compression. Accordingly, its stress-strain relationships in elastic and plastic states are controlled by Eq. (8. a) and Eq. (8. b), respectively.
In these equations, σ and ε are stress and strain tensors, while dσ and dε represent their differential changes, respectively. For the elastic behavior in planestress state, generalized Hooke's matrix, C s , is defined using steel's modulus of elasticity, E s , together with Poisson's ratio, υ s ; that is Furthermore, concerning plastic behavior, the constitutive tangent matrix of the steel can be expressed as follows; In the above expression, f(σ) indicates the yield function, which is considered based on the Von-Mises criterion (Kubosek and Vaskova 2017) herein, using the uniaxial yield strength, σ y , according to Eq. (9. b).

Concrete
Prior to cracking and elastic behavior, the linear stressstrain relationship for this material can be expressed as follows: By assuming a homogeneous and isotropic behavior, the stiffness tensor in the linear state of concrete, C c , is also defined according to Eq. (11). In this matrix equation, E s and υ s represent the modulus of elasticity and Poisson's ratio of uncracked concrete, respectively.
If the maximum tensile stress inside an element exceeds the concrete rupture modulus (i.e., f t > f r ), as shown in Figure 3, with the appearance of cracks, Eq.
(10) can no longer be used. In this case, a twodimensional coordinate system is considered whose x-axis is parallel with crack direction, and the y-axis is perpendicular to the direction of the crack (the same direction of the principal tensile stress). The incremental constitutive equation for the cracked condition is expressed as follows: In Eq. (12), β is the shear retention factor, which depends on the aggregate interlocking and the arching phenomenon that determines the conditions of the crack face. It is also used to surmount the difficulties with a singular matrix and varies from 0 to 1. By transforming the tangential stiffness of concrete C t c in the global coordinate system (x 1 ,x 2 ,x 3 ), the governing equation takes the form below: where the transformation matrix,T is defined by Eq. (14).
The nonlinear compressive behavior of concrete is also considered by the stress-strain relationship provided by Hognestad (1951), as illustrated in Figure 4. By using fiber elements, the modeling accuracy of the contribution of concrete between cracks can be improved through models like Vecchio and Collins according to Figure 5 (1986). The models employed to consider the tension of concrete between cracks are generally expressed based on the rebar ratio of the section (Stramandinoli and La 2008). For more details, please refer to references (Edifícios and Concreto 2002; Filipe and Soares 2010; Stramandinoli et al. 2007). Various models have also been provided for the concrete under tension based on the experimental works (Carreira and Chu 1986). Before the concrete cracking, the material is assumed to be in the linear elastic mode, and after it, the effect of strain hardening is expressed using Eq. (15) (Stramandinoli and La 2008).
In Eq. (15), f ct denotes the tensile strength of the concrete, and ε cr is its corresponding strain. Moreover, ε y is the strain corresponding to the yield of rebars, after which the σ ct stress becomes zero. In addition, α is the exponential reduction parameter, which is generally a function of the rebar percentage of the section (ρ) and the elastic modulus ratio of steel to concrete (n), which is given by Eq. (16):  When using the equation above to calculate ρ for reinforced concrete sections, the effective area in the tensile part of the section should be defined using Eq. (17) (C. MC90 1993): In this equation, b is the section width, and h and d denote the total height and effective height of the section, respectively. Moreover, k x denotes the height of the neutral axis in the section. In practice, Eq. (18) is generally approximated as the following:

Transverse reinforcement
By assuming that stirrups are distributed in the concrete layer (smeared transverse reinforcement), the stress-strain relationship for transverse rebars with linear behavior can be expressed by Eq. (19).
Herein, the linear constitutive matrix, C sv , is defined using modulus of elasticity, E sv , cross-sectional area, A v , and the longitudinal spacing, s, of transverse rebars, and width of the section, b (as shown in Figure 6) as follows; It should be noted that in the case of incremental stress-strain relationship, the constitutive tangent matrix of transverse reinforcements is precisely equal to Eq. (20), i.e., C t sv ¼ C sv . As a rule, by assuming the embedded transverse reinforcement, the stress-strain relationships for concrete are modified by Eq. (21).

Finite element formulation
A three-node beam element is used to derive the governing equilibrium equations in local coordinates, as shown in Figure 7, where the middle node is considered to increase the accuracy of the interpolation procedure. In this figure, variable ξ is a dimensionless coordinate parameter whose relation to the longitudinal coordinate of the element, x, is expressed as � ¼ 2x=l.
If the nodal displacement of each element is expressed by vector δ, the relationship between the displacement field inside the element and nodal variables can be obtained using the shape functions, N.
It is worth mentioning that the cubic Hermite and quadratic Lagrangian shape functions are employed to interpolate the components w and u of the displacement field. By substituting Eq. (22) into Eq. (7) the strain-displacement matrix, B, can be introduced as follows: The total potential energy (TPE), П, stored in each element shall be expressed to derive the equilibrium equations.   where V e and S e represent the volume and boundary surface of each element, while traction and body force vectors are denoted by F and f, respectively. Subsequently, minimizing the TPE leads to the following equilibrium equation for each beam element.
In Eq. (25), element stiffness matrix, k e , is given by (26) where n c ,n s , and n sc are the number of concrete layers, longitudinal steel layers, and concrete layers with embedded transverse reinforcement, respectively. The load vector, p e , is also defined by the following integral expression: Eventually, through the assembling procedure, the system's set of force equilibrium equations can be obtained by Eq. (28), in which n e denotes the number of elements.

Numerical integration
These values are found for the elements after formulating the stiffness matrix and load vector using the Gaussian integration. Then, in each Gaussian point, the structural matrix and stress vector is calculated using the sum of the values related to the concrete and rebar layers. After calculating the stiffness matrix and load vector in the elements, the middle degree of freedom is removed through a static condensation procedure, and the stiffness matrix and load vector are then assembled in the global coordinates. The numerical integration in each element is performed using three Gaussian points, in each of which the section is discretized to the concrete layers covered with longitudinal rebar layers (fiber model). Here, it is assumed that no segregation occurs between rebars and concrete, the rebars are properly restrained inside the concrete, and each layer is in a uniaxial stress state.

Methodology
In this section, after reviewing the conventional forcebased solution for analysis of a set of nonlinear FE equations, an energy-based solution is presented. The section focuses primarily on the fundamental concepts and equations underlying this method. Finally, the section concludes with a flowchart illustrating the proposed method's computer implementation.

Force-and energy-based solutions
These step-by-step approaches work on Eq. (25) that states the force equilibrium of the FE assemblage. By assuming that the current state of the system at the ith step is given by i D, the state of the system in the next step, i+1 D, may be computed by where the incremental displacement, i+1 ΔD, within each step is calculated by employing an iterative scheme that is expressed by Eq. (30).
In this equation, m refers to the number of iterations, and i+1 ΔF stands for the incremental load at the next step. This procedure is illustrated in Figure 8.
As an alternative approach, the energy equilibrium of the system can be considered instead of using the force equilibrium equation. The method presented in this research is primarily based on the modified energy method (MEM), which was first introduced in references (Abad et al., 2017; Jalili Sadr Abad, Mahmoudi, and Dowell 2018; Jalili Sadr Abad and Mahmoudi 2018) for the dynamic analysis of structures with nonlinearities. The current study proposes a numerical  technique to determine nonlinear problem responses in solid mechanics by developing the abovementioned method. In fact, the performance of this numerical scheme, which uses the system's energy equilibrium relationships, is compared with other common numerical methods, such as the Newton-Raphson method. The real analysis of structural problems involves a dynamic procedure in which the structural behavior changes over time. However, not to deal with analytical complexities, a quasi-static process is used when external loads are gradually applied to the structure, as assumed in this study. In quasi-static analyses, the load is applied to the system as incremental steps. In this case, t represents a pseudo-time variable that expresses the intensity of the loading in a given load step and should not be confused with the time variable used in nonlinear dynamic analyses.
In this study, instead of Eq. (28), which indicates the force balance of the system, an energy-based formulation is used as the governing equation for a nonlinear structural problem. In general, various types of energies in a static system include (I) the energy due to resistant forces of the structure, F int , and (II) the energy or work from the external loading applied to the structure, F ext , which are shown in this study, by E int and E ext , respectively. According to the law of conservation of energy, Eq. (29), these two values of energy must be equal at any moment.
In order to obtain an incremental form, if the above relationship is written in two consecutive time steps (t i , t i+1 ), we have The internal energy variation resulting from the work of the internal forces, ΔE int , is expressed for a nonlinear system as integral below.
This equation can be discretized using the trapezoidal rule as follows, in which superscript i indicates the time step number.
Expressing the internal force in the current step, i+1 F int , using the internal force in the previous step, i+1 F int , and the tangential stiffness matrix, t K, Eq. (33) takes the form below.
On the other hand, the external energy variation resulting from applied loads, ΔE ext , may also be defined by Eq. (36).
Using the definition of the system's velocity, V(t), the Eq. (36) can be rewritten as follows: with the following discretized form: Now by substituting Eq. (35) and (38) into the energy equilibrium equation, Eq. (32), it can be written: (39) On the other hand, the displacement changes, ΔD, can also be written using the Euler formula as Eq. (40), where β is the coefficient of velocity participation in determining the system's displacement, which controls the type of integration.
By substituting Eq. (40) into Eq. (39) and performing some mathematical simplification, a set of quadratic equations, Eq.(41), can be obtained using the current velocity of the system, i.e., i+1 V: The coefficients can be expressed as follows: Given that the existence of imaginary velocities in a structural system has no physical meaning, the discriminant of the quadratic equation, Δ, must always be positive. Considering this value is generally related to system properties, the nature of loading, and discretization parameters, an additional condition for the stability of the obtained numerical solutions can be obtained. It is worth noting that this is an advantage of this formulation that is not available in other commonly used force-based methods. Given that the value of the discriminant depends on parameter β, it can be proved that the formulation provided for β = 0 never creates imaginary velocities during the analysis. However, for other values of this parameter (i.e., 0 < β ≤ 1), the sign of this quantity needs to be checked during the analysis, and if the value is negative, the size of the time interval, Δt, should be reduced. By applying Eq. (41), it is clear that as the time step approaches zero (Δt → 0), the delta value of the quadratic equation always becomes positive.
The fundamental question that arises for readers here is that given the quadratic nature of the velocity equation and the creation of two velocities at any given time instance, how the real velocity of the system can be chosen from these two roots. For this purpose, according to the reference (Jalili Sadr Abad, Mahmoudi, and Dowell 2018), based on the continuity of velocity changes, the actual velocity at any time can be selected as the one that is closer to the velocity in the previous time step. Mathematically, if the two velocities calculated in the current step are called v 1 and v 2 , in this case, it can be written:

Convergence criteria
A convergence criterion is used to stop iterations in each step once a level of accuracy is obtained. In fact, the nonlinear analysis consists of several linear solutions in each iteration and is based on the residual force. In incremental-iterative methods, the convergence criteria are used after each iteration. The two criteria, CRI1 and CRI2, were used in this study to evaluate the convergence. 1. The first convergence criterion, CRI1: This criterion compares the normal residual forces with the normal forces applied at each increment. According to Eq. (44), we have F i and G i denote the applied components and the equivalent force of internal stresses in each degree of freedom in node i. In addition, n dof is the total number of a structure's degrees of freedom. 2. The second criterion, CRI2: This criterion compares the work of residual forces with the work of applied forces. According to Eq. (45), we have R ¼ W RES W F < Toler F i , G i , and n dof are defined for the first criterion. ΔU i is partial displacements in a given degree of freedom for node i, and ΔF i is the incremental force. The tolerance (Toler), which is defined by the user, is usually within the range of 0:001,0:01 (Frantzeskakis and Theillout 1989). Both CRI1 and CRI2 are appropriate for investigating convergence. The CRI2 criterion is mainly used in reinforced concrete structures. The convergence tolerance should also be chosen with caution. A high convergence rate would be uneconomical by increasing the computation time, while a low convergence rate leads to a considerable error in results. Finally, according to Figure 9 the computer implementation of the proposed method for analyzing a nonlinear system can be performed using the steps provided in the following flowchart.

Results and discussion
Three examples of RC beams previously studied in other works are presented in the following section, and the method presented in this research is applied to them. The first example is an RC beam with only tensile longitudinal reinforcements, but in the second and third examples, the effect of compressive and transverse rebars are also considered, respectively. The results obtained from the proposed method are validated with those of nonlinear force-displacement analyses in other studies.

A simply supported beam with only tensile reinforcement
The first example is an RC beam with two pinned supports subjected to a mid-span concentrated load, as shown in Figure 10a. The reinforcement of this structural member consists of only three longitudinal rebars with a diameter of 20 mm, which are placed at the bottom of the section (according to Figure 10b) and act under tensile stress. This structure has already been studied experimentally and numerically by Sarkhosh et al. (2012) and Slobbe et al. (2012), respectively. The assumed values for material properties of both concrete and steel in this example are given in Table 1.
In the first step in analyzing the problem, given the existing symmetry, only half of the structure is considered, as shown in Figure 11. Here, it should be noted that the mid-span displacement of the beam, i.e., Δ, is monitored under the load increments as the main degree of freedom.
As the shear span-to-depth ratio (a/d) is equal to 2.9, tension and shear-compression failure modes would be anticipated for the structure. According to table 2 Some of the assumptions made in FE analysis of the current problem can be expressed as the following; Based on the assumptions made and the formation of nonlinear equations governing the problem, these equations are solved to find the system's response. Accordingly, as stated in the flowchart presented in the previous section, the force variation against displacement at the middle of the beam is plotted in Figure 12.
As shown in Figure 12, the results obtained from the proposed numerical method have a good agreement with the experimental results, and they have converged to the actual response of the structure with    acceptable error. This is also evident in terms of the method's superior performance when compared to the method used by Slobbe et al. In addition, the numerical instability points (jumps in the curve) have been minimized appropriately using the present numerical method because, as mentioned before, an analyst can control the stability of the solution by monitoring the delta sign in the method. According to Figure 12, as we know, when an element with a larger size is used, the structure's stiffness increases to some extent, showing a relatively more ductile mode. Meanwhile, with the reduction in the size of elements, the ultimate load is generally reduced, resulting in a less ductile structure. In other words, the displacement at the ultimate load decreases with the element size. The ultimate load declines for a beam with a few rebars as the number of elements grows (Kheyroddin and Mortezaei 2011).
In the used FE method, after reaching the maximum capacity, a large number of elements crack, and their stiffness becomes zero. On the other hand, since the stiffness of the structure, which is updated in each step, is significantly reduced, the FE model becomes unstable, resulting in a sudden change in the slope.

A pinned-support beam reinforced with tensile and compressive rebars
As shown in Figure 13a, an RC beam with a span length of 2.4 m subjected to two symmetrical point loads creating a pure bending state along the middle onethird of the span length is considered another example. The cross-section of this member is a rectangle with a width of 120 mm and a height of 300 mm, which has been reinforced by 2φ5 and 3φ10, respectively, at its top and bottom as displayed in Figure 13b.
In this case, the mechanical properties of the materials, in this case, shown in Table 3, are selected according to an experimental work performed by Alvares (1993). In addition, the structural behavior of a beam quite similar to this example has been numerically examined by Oliveira, Ramalho, and Correa (2008), which is used here to validate the results.
As with the previous example, given the symmetric geometry and loads of the structure, only half of it is modeled, as depicted in Figure 14. It is noteworthy that in this case, the displacement in the middle of the beam is considered the master degree freedom of the FE assemblage.
The following are some of the assumptions made for the FE formulation according to table 4;   (Alvares 1993;Oliveira, Ramalho, and Correa 2008 Given the prior assumptions and the nonlinear set of equations derived from Eq. (28), the proposed numerical technique is used to determine the structure's structural behavior in terms of the load-displacement diagram. As a result of following the steps outlined in Section 3, Figure 15 illustrates the applied load versus mid-span deflection of the structure. According to the figure, the curve with fewer elements has provided a solution with acceptable accuracy. In fact, the results demonstrate that the effect of element size on the ultimate load is not significant, i.e., the beam's response is not substantially dependent on the number of elements.
From Figure 15, it is observed that the numerical response obtained from the present method for the displacement of up to 6 mm is almost consistent with the two curves related to the test data (Alvares 1993) and those reported by Oliveira, Ramalho, and Correa (2008). At greater displacements, both curves of the numerical analyses of the problem are significantly separated from the experimental curve. As can be seen, the maximum value of the diagram (failure load) is estimated with acceptable accuracy by the proposed method.

A simply supported beam with longitudinal and transverse reinforcements under a four-point bending test
The third example is a pinned 2.3 m long RC beam constrained by two pinned supports subjected to two point loads, as shown in Figure 16a. As can be seen in the figure, the width and height of the rectangular cross-section are 200 and 150 mm, respectively, where reinforcing rebars are selected as below (Figure 16b): • Top longitudinal rebars: 2ϕ8 • Bottom longitudinal rebars: 2ϕ16 • Stirrups: ϕ6@75 mm The material properties of the concrete and rebars are presented in Table 5. These properties are obtained from the experiments carried out by Rahimi and Hutchinson (2001). It should be noted that Lezgy-Nazargah (2018) also used this example in order to find the nonlinear force-deflection diagram of the system numerically.
Similar to the previous two examples, to reduce the computational effort and the number of unknowns of the problem (reducing the dimensions of the global   stiffness matrix), only half of the structure is considered, as shown in Figure 17, due to symmetry in the load and geometry conditions of the system. Other assumptions made in the formation of FE equations governing the problem can be summarized in the following items: According to the assumptions above, along with the nonlinear FE equations, these equations are solved to obtain the structural behavior of the system using the presented computational technique. The corresponding force-displacement curve is shown in Figure 18.
As shown in the previous examples, Figure 18 shows that the proposed numerical method is more accurate than the numerical data obtained from the study. The main advantage of the energy method, in this case, is that the analyst can choose displacement increments as the analysis progresses by using the delta sign of the energy method's characteristic equation, which allows for the avoidance of divergent responses.
When larger elements are used, the ultimate load can be estimated within an acceptable range (Kheyroddin and Mortezaei 2011). However, given the numerical method introduced in the article and the evaluation of the load-displacement curves in Figure  18, it can be concluded that the accuracy of the numerical method increases with the number of elements, the results obtained from it agree with the experimental results, and they converge with acceptable errors until reaching the real response of the structure. The suggested FE model can be applied for beam columns. It can also analyze reinforced concrete beams and   columns. Moreover, this method can be developed for structures like reinforced concrete frames and those with frame-wall systems.
It should be noted that the reason behind the high slope in the elastic area is the elastic modulus of the concrete. The equation suggested by the code does not always agree with reality and results in the overestimation of the elastic stiffness. Another reason is the ignoring of the sliding between rebars and concrete. The article has assumed the full connection between them; no segregation occurs between the rebars and concrete, and the rebars are properly restrained inside the concrete. Moreover, the software ignored the mechanical interaction between rebars and their surrounding concrete (bond-slip) compared to the experimental sample (Slobbe, Hendriks, and Rots 2012).

Conclusion
This study proposed an energy-based computational technique to solve nonlinear equations concerning the finite element analysis of reinforced concrete beams. It was demonstrated that the presented methodology could easily be applied to various types of material nonlinearities involving concrete and steel rebars and analysis of large deformations (nonlinear geometric effects). Unlike the force equilibrium equations, the discrete equations in this approach were a quadratic form. Therefore, by analyzing the delta sign of the quadratic expression, an analyst can ensure the stability and accuracy of the obtained response during nonlinear analysis, a unique feature in finite element analysis of nonlinear systems that can distinguish this method from other existing methods. Regarding the accuracy of this numerical scheme, by comparing the results obtained from this method for various types of reinforced concrete beams (with different types of longitudinal and transverse reinforcements), it was observed that the forcedisplacement curves obtained from this method were in good agreement with experimental and numerical data of other studies. Finally, the authors propose that future works extend the proposed energy method to nonlinear analysis of other structural engineering applications.

Disclosure statement
No potential conflict of interest was reported by the author (s).