Nonlinear Finite Element Analysis Formulation for Shear in Reinforced Concrete Beams

: This study suggests a novel beam-column element formulation that utilizes an equilibrium-driven shear stress function. The beam shear is obtained from the bi-axial states of micro-planes, through matrix condensation and zero vertical traction assumptions. This properly remedies the shear sti ﬀ ening of a one-dimensional beam-column element, keeping its degrees of freedom to a minimum. For veriﬁcation of the proposed method, a total of seven shear test results of reinforced concrete (RC) beams were collected from the literature, in which the key variables were the reinforcement ratio, the presence of shear reinforcement, and section shape. The advantages are clearly shown in the shear stresses distributions being accurately described and the global load-displacement relations being successfully obtained and matching well with various test results. The proposed model shows satisfactory descriptions of the monotonic load-displacement response of the RC beams failing in multiple modes that vary from diagonal-tension to ﬂexural-compression. In addition, more accurate and reliable information of sectional responses including sectional shear deformation and stresses is collected, leading to better prediction of a potential shear failure mode. Finally, the advantages of the proposed model are demonstrated by comparing the analysis results of an RCT-beam by using the di ﬀ erent shear assumptions that include the constant and parabolic shear strains, constant shear ﬂow, and the proposed shear stress function.


Introduction
The modern seismic design provisions for reinforced concrete (RC) structures enforce ductility requirements to meet target limit states, avoiding brittle failure [1]. As a typical shear failure of an RC member is often very brittle, more accurate descriptions of shear behavior are imperative in achieving proper safety goals. Therefore, the shear performance of RC members should be clearly investigated and adequately modeled. In general, the global responses of reinforced concrete (RC) structures are often obtained by means of a frame analysis model, where the members can be idealized by a one-dimensional beam-column element. The conventional method in considering beam shear utilizes the Timoshenko beam with a shear coefficient. Most shear coefficients known today are strain-driven functions, such as parabola shear strains [2][3][4][5]. Such strain-driven shear functions do not properly describe nonlinear shear evolution. As a result, the conventional method rarely shows the shear strength degradation and shear failure that are often observed in loaded RC beams and columns. Many shear models have been proposed to improve such shear simulations, as follows.
Petrangeli et al. [2,3] developed a flexibility-based element to model shear behavior and its interaction with axial force and bending moment in RC beams and columns. Based on the fiber section discretization, the state of a section was determined by numerically integrating biaxial constitutive relations over the cross section. While the longitudinal strain field was determined from the plane section assumption, the shear distributions over the cross section were defined by strain-driven shear functions, such as constant shear strain, or parabolic shear strain functions. Despite the use of the micro-plane where the two-dimensional realization of shear was possible, the stain-driven shear functions gave little hint of shear degradation. A similar development was made by Taylor et al. [4] but based on the mixed finite element method using the Hu-Washizu variational form, and later extended by Filippou and Saritas [5] for RC members. As in Petrangeli's model [2], the shear strain distribution over the cross section was assumed to be parabolic.
The equilibrium-driven stress function was first introduced in the beam formulation by Vecchio and Collins [6] called the dual-section analysis. This method well predicted the shear stress distributions over the cross section of an RC beam. The shear flow of the cross section was determined from the normal stress difference between two cross sections. Bentz [7] proposed the longitudinal stiffness method, where the shear flow or the normal force gradients were directly related to the sectional properties and the sectional shear force. It was a challenging task to write a consistent element formulation that utilized equilibrium-driven stress functions. The method worked best in the sectional analysis. However, both Vecchio [6] and Bentz [7] formulations have not yet been studied in the conventional finite element method.
Recently, finite element analysis software programs [8][9][10] have been developed based on the finite element formulations proposed by the afore-mentioned researchers, and they were utilized by many researchers for their numerical simulations on shear behavior of RC members. Ju et al. [11] have conducted finite element analyses to investigate the shear stress concentration phenomenon near a tapered cross section, and Chen et al. [12] carried out a hybrid simulation testing on a precast concrete shear wall, which combines the structural test and numerical simulation. In addition, Roudane et al. [13] performed numerical simulations on masonry infilled RC buildings considering construction stages by using ABAQUS software.
Many shear models available today use strain-driven functions. Unfortunately, the strain-driven shear functions cause redundant constraints resulting in an over-estimation of shear stiffness and strength. To overcome such a shortcoming, this paper aims to propose a force-based solution adopting the equilibrium-driven shear stress function. The longitudinal stiffness method is rewritten in the context of the conventional finite element method, and a consistent beam formulation is derived based on the Timoshenko beam theory.

Shear Solutions in One-Dimensional Beams
Nodal displacements of a typical beam-column element are obtained by solving differential equations expressed in terms of the longitudinal coordinate. Therefore, the typical beam-column is often referred to as a one-dimensional element, although the nodal degrees of freedom can be defined in two-or three-dimensional spaces. The bending of a slender beam is accurately solved using the plane section assumption in the realization of a one-dimensional element, a.k.a. the Bernoulli-Euler beam [14,15]. This simplification yields the stiffer behavior of a short beam, as it neglects shear deformation, which imposes additional rotation of the cross section. The Timoshenko beam suggests a solution to this problem. The equilibrium and strain-displacement equations of the Timoshenko beam can be derived as follows [16], where, n = axial force field; m = bending moment field; q and p = applied transverse and axial load per unit length, respectively; v = shear force field; κ = curvature field; γ = shear strain field; ε = axial strain field; w = transverse displacement field; θ = rotation field; and u = axial displacement field. A prime denotes differentiation with respect to ξ, which locates a position in the reference coordinate system, as depicted in Figure 1.
where, n = axial force field; m = bending moment field; q and p = applied transverse and axial load per unit length, respectively; v = shear force field; κ = curvature field; γ = shear strain field; ε = axial strain field; w = transverse displacement field; θ = rotation field; and u = axial displacement field. A prime denotes differentiation with respect to ξ , which locates a position in the reference coordinate system, as depicted in Figure 1. The Timoshenko beam can accurately describe the shear behavior of homogeneous and linear beams. The sectional displacements, including axial and transverse displacements and rotation, can define the normal deformations of each point or layer of the cross section. However, as the required constitutive model to calculate the corresponding stresses is uniaxial, ad-hoc functions, such as shear functions, are required to define shear profiles over the cross section, mostly with parabolic shapes or constant ones. Moreover, the beam shear is uncoupled from the flexural responses, and the description of nonlinear shear evolution of the composite beams that are composed of nonlinear materials is hardly achievable. To resolve this problem, this study suggests a new element formulation to determine shear stress distributions over the cross section.

Shear Formulations for Beam Section
The longitudinal stiffness method [7] determines the shear stress distributions based upon equilibrium. Therefore, the method is rigorous in solving nonlinear beams, such as RC beams. The longitudinal stiffness method is rewritten in the context of the finite element method, where more consistent element formulations can be derived.
The shear stresses acting on the surface of a layer of the cross section should be equal to the horizontal surface shear stresses acting between layers of the beam from equilibrium. This is the basis by which the shear flow is used to calculate the shear stress distributions. To demonstrate this, let us consider a free body diagram (distributed load omitted), as shown in Figure 2. The shear flow-the horizontal shear force per unit distance along the longitudinal axis of the beam-can be evaluated from equilibrium as: ( ) x df ζ denotes the difference of longitudinal stresses between cross-sections mn and m1n1, ( ) w b ζ is the width of the layer at the depth of ζ , and h is the overall height of the section. The horizontal shear stress at ζ is obtained by dividing the shear flow by the width of the layer at ζ ; i.e., The Timoshenko beam can accurately describe the shear behavior of homogeneous and linear beams. The sectional displacements, including axial and transverse displacements and rotation, can define the normal deformations of each point or layer of the cross section. However, as the required constitutive model to calculate the corresponding stresses is uniaxial, ad-hoc functions, such as shear functions, are required to define shear profiles over the cross section, mostly with parabolic shapes or constant ones. Moreover, the beam shear is uncoupled from the flexural responses, and the description of nonlinear shear evolution of the composite beams that are composed of nonlinear materials is hardly achievable. To resolve this problem, this study suggests a new element formulation to determine shear stress distributions over the cross section.

Shear Formulations for Beam Section
The longitudinal stiffness method [7] determines the shear stress distributions based upon equilibrium. Therefore, the method is rigorous in solving nonlinear beams, such as RC beams. The longitudinal stiffness method is rewritten in the context of the finite element method, where more consistent element formulations can be derived.
The shear stresses acting on the surface of a layer of the cross section should be equal to the horizontal surface shear stresses acting between layers of the beam from equilibrium. This is the basis by which the shear flow is used to calculate the shear stress distributions. To demonstrate this, let us consider a free body diagram (distributed load omitted), as shown in Figure 2. The shear flow-the horizontal shear force per unit distance along the longitudinal axis of the beam-can be evaluated from equilibrium as: where, − h 2 ≤ ζ < h 2 ; d f x (ζ) denotes the difference of longitudinal stresses between cross-sections mn and m 1 n 1 , b w (ζ) is the width of the layer at the depth of ζ, and h is the overall height of the section. The horizontal shear stress at ζ is obtained by dividing the shear flow by the width of the layer at ζ; i.e., Appl. Sci. 2019, 9,3503 4 of 20 The normal stress gradient or the stress differences between the two sections d f x can be defined as: where, the first term [10] is multiplied to obtain only the scalar normal stress, d is the 2 × 2 sectional stiffness matrix of the cross section, which will be defined later, and dε is the 2 × 1 sectional strain vector. The last term is the sectional strain changes due to the moment change ∆M as shown in Figure 2a, and related to the 3 × 1 full strains as, dε = ςde And where, ψ is the shear strain function. Depending on how this is defined results in a constant, parabolic, or other shear distribution. For example, if ψ(ζ k ) = 1, a constant shear distribution is obtained.

Two-Dimensional Realization of Beam Shear
As shear is at least two-dimensional, two-dimensional representation of shear is indispensable to obtain more realistic shear distributions. The micro-plane concept is the simplest and widely accepted method among beam shear models because it can be realized in a one-dimensional beamcolumn element. The micro-plane describes bi-axial states, where two correlated uniaxial states are defined in their principal directions. For an RC beam, the bi-axial states can be defined as follows.

Concrete Material
There are many constitutive models to simulate the concrete behavior [17][18][19], and those models are widely used in the finite element analysis. In all analyses presented in this paper, the uniaxial behavior of concrete in compression is obtained by applying the modified Kent and Park model [20] to account for the confining effects provided by stirrups. That uniaxial behavior is scaled by using the compression softening relation suggested by Vecchio and Collins [21]. In addition, other concrete models with confinement effect [16][17][18] can be interchangeably used. To describe tension stiffening, the shear cracking stress of concrete is calculated by using the ACI expression 0.17 c f ′ , where c f ′ is in MPa [22]. The post-cracking envelope is obtained by using the equation proposed by Collins and Mitchell [22]. Figure 3 illustrates the uniaxial behavior of concrete material. It is assumed that this Knowing that Equation (5) defines the sectional strain changes due to the moment change ∆M = V leads to: Herein, K s is the full stiffness matrix of the cross section, a.k.a. the longitudinal stiffness matrix. The discretization of Equation (3) and writing the shear stress at the k-th fiber leads to: where, ℘(ζ k ) = k j=1 h f ζ j · d f x ζ j · b w ζ j , and h f ζ j is the height of the j-th fiber. Finally, the shear stress function at the k-th fiber is obtained from: It should be noted that the shear distributions are determined from the sectional equilibrium. For equilibrium-driven equations, it is extremely difficult to derive a consistent element formulation in the conventional finite element method, which is based on the displacement-based formulation. Therefore, the authors seek the answer from the force-based formulation that strictly satisfies equilibrium.

Two-Dimensional Realization of Beam Shear
As shear is at least two-dimensional, two-dimensional representation of shear is indispensable to obtain more realistic shear distributions. The micro-plane concept is the simplest and widely accepted method among beam shear models because it can be realized in a one-dimensional beam-column element. The micro-plane describes bi-axial states, where two correlated uniaxial states are defined in their principal directions. For an RC beam, the bi-axial states can be defined as follows.

Concrete Material
There are many constitutive models to simulate the concrete behavior [17][18][19], and those models are widely used in the finite element analysis. In all analyses presented in this paper, the uniaxial behavior of concrete in compression is obtained by applying the modified Kent and Park model [20] to account for the confining effects provided by stirrups. That uniaxial behavior is scaled by using the compression softening relation suggested by Vecchio and Collins [21]. In addition, other concrete models with confinement effect [16][17][18] can be interchangeably used. To describe tension stiffening, the shear cracking stress of concrete is calculated by using the ACI expression 0.17 f c , where f c is in MPa [22]. The post-cracking envelope is obtained by using the equation proposed by Collins and Mitchell [22]. Figure 3 illustrates the uniaxial behavior of concrete material. It is assumed that this uniaxial behavior discussed thus far describes the orthotropic concrete material in the principal directions. The composition of material tangent stiffness for an orthotropic material is discussed as follows.

Force-Based Element Formulation of the Timoshenko Beam
For the realization of the aforementioned shear functions and two-dimensional beam shear, first, the element formulation is derived followed by the presentation of the solution strategy. The forcebased element formulation of the Timoshenko beam is derived by directly integrating the governing equations. For an orthotropic material such as concrete, by neglecting the Poisson's ratio, the material tangent in the principal directions can be written as: where, E 1 and E 2 are the Young's modulus in each principal direction, respectively, which can be calculated by using the ACI expression 4700 f c , and G is the shear modulus. The latter can be determined by: where, f 1 and f 2 are the principal stresses, and ε 1 and ε 2 are the principal strains. Equation (11) satisfies the condition of invariance under the principal axis rotation [23]. The orthotropic material stiffness should be transformed to the current local coordinate system of the element by using the following relation: where the transformation matrix is defined as: The failure criterion of concrete material can be found in the modified compression field theory (MCFT) formulation [7], where the two-dimensional plane reinforced concrete response was simplified with two uni-axial stress-strain relations in the two principal directions. In the principal directions, the compression failure occurs at the compression strain of −0.003. The steel was considered as smeared to the corresponding concrete layers, and the confining effects were also reflected, as described in Figure 3.

Steel Material
A bi-linear behavior with a certain strain hardening ratio is assumed for steel. The material tangent for steel can be written as:d where, the subscripts x and y represent the longitudinal and vertical directions, respectively. The former is the axis of the longitudinal bars, and the latter is that of the stirrups. In a beam-column element, the aforementioned two-dimensional states should be driven from the uniaxial strain of each layer of the cross section. This involves the matrix condensation with the assumption of zero traction in the vertical direction and the shear stress function that was presented before.

Force-Based Element Formulation of the Timoshenko Beam
For the realization of the aforementioned shear functions and two-dimensional beam shear, first, the element formulation is derived followed by the presentation of the solution strategy. The force-based element formulation of the Timoshenko beam is derived by directly integrating the governing equations.

Direct Integration of the Equations of Equilibrium
Define the following functions as integrals of the applied transverse and axial loads [16]: The integration of Equation (1c) gives: Applying the integration by parts, and integrating the equation above, leads to: By solving for the constants and rearranging the equation above, we obtain: where, ν 0 = ν(0). From Equation (19) and knowing that m = −ν, we can obtain: Along with Equation (15), substituting L into Equation (20) and solving for v 0 give: Rearranging Equation (19) leads to: The integration of Equation (1a) and rearrangement give: where, n 0 = n(0). Let us denote the force field and the discrete element force parameters as: The force field can be written as: where, b T (ξ) and the particular force field σ p (ξ) are defined as: The end forces can be defined as: as shown in Figure 4. This can be expressed in terms of the essential force parameters and the applied loads as [16]: where, the equivalent nodal load vector is denoted: and the matrix L is given by: Appl. Sci. 2019, 9, 3503 9 of 20 (a) Element force.
(b) Nodal displacement and element end force sign convention.

Direct Integration of the Differential Strain-Displacement Equations
Integrating the strain-displacement equations of Equations (1b), (1d), and (1f) gives: This can be written in matrix form as: Beam-Column Element

Direct Integration of the Differential Strain-Displacement Equations
Integrating the strain-displacement equations of Equations (1b), (1d), and (1f) gives: This can be written in matrix form as: The detailed derivation of the above equations can be found in reference [24], and it was developed by adopting the direct integration [25].

State Determination
In the direct stiffness method, the equations of motion, which are constructed by assembling the element state quantities, can be solved by using any Newton-type iteration method. This reduces the difference between the internal force and external force through iterations. This also holds for the nonlinear flexibility-based method, as the element stiffness matrix is obtained by inverting the element flexibility matrix. The latter is obtained from the integration of the quadratic product of the section flexibility and force interpolation function matrix. However, in contrast to the stiffness-based method, the flexibility-based method requires an additional iteration loop to satisfy compatibility conditions. This is discussed in the following.

Fiber State Determination
The section strain is discretized to obtain the strain at the k-th fiber by using the following equation: where, ε mik ≡ ε m (ξ ik ); ξ ik ≡ ξ i (ζ k ); ξ ik is the vertical coordinate, as shown in Figure 5; where ε x (ζ) and γ xy (ζ) are the normal and shear strains, respectively; ς mik is the base function as defined earlier. The fiber stress ϕ is then obtained from: where, ϕ ≡ f x (ζ), ν xy (ζ) T ; f x (ζ) and ν xy (ζ) are the normal and shear stresses, respectively. If a biaxial constitutive law is used for the material, additional strain and stress fields should be considered; those are ε y and f y for strain and stress, respectively. The complete biaxial state of an infinitesimal element is depicted in Figure 6, and the stress-strain relations are given in Equation (30) in the incremental form after linearization.
where,φ ≡ f x , f y , ν xy ;ε ≡ ε x , ε y , γ xy . A general approach to retrieve the additional quantities assumes f y = 0 at any location of a beam. Therefore, through static condensation, the biaxial state can be reduced to: Appl

Solution Strategy
The flexibility-based method works together with force interpolation functions, for which reason it is also known as the force-based method. This method allows the solution of the system of equations without shape functions. A similar principle is taken in the determination of the fiber state,

Solution Strategy
The flexibility-based method works together with force interpolation functions, for which reason it is also known as the force-based method. This method allows the solution of the system of equations without shape functions. A similar principle is taken in the determination of the fiber state,

Solution Strategy
The flexibility-based method works together with force interpolation functions, for which reason it is also known as the force-based method. This method allows the solution of the system of equations without shape functions. A similar principle is taken in the determination of the fiber state, where shear stress functions (or shear interpolation functions) are known. The predictor and corrector phases, which determine the section shear force V mi for the known section strain e mi , will be explained in the following.

Predictor Phase
Step 1: The trial section shear force is determined from the following relation: where, ∆V mi = 0 0 1 K ν−1 mi ∆e mi .
Step 2: The axial strain is evaluated by using the plane-section assumption but the vertical and shear strains are calculated through matrix condensation; i.e., in an incremental form: where, ∆ν xy = ν trial xy − ν v−1 xy . Since all the strain values are known, the corresponding stresses can be evaluated by using any biaxial constitutive model. Then, the corrector phase follows, because the shear stress obtained from a constitutive relation will differ from the trial shear stress.

Corrector Phase
Step 3: In the corrector phase, the shear strain is corrected in a way that reduces the shear stress residuals. However, in the flexibility method, strain residuals instead of stress are integrated to give the counterpart in the upper level. Therefore, the section shear strain residual is obtained as: and ω is the shear stress functions discussed later. The superscript R denotes residual.
Step 4: The section shear force is now corrected by the following relation; i.e., where, V R ≡ K v mi γ R (ξ mi ).
Step 5: The fiber state should also be corrected as given below: where, 1 ≡ [ 0 0 1 ]; ι ≡ diag 1 1 0 ; ϕ v mik ≡ ϕ(ε v mik ) is the fiber resisting stress obtained from a constitutive relation for a known ε v mik . Note that each step is indicated with a circular number at the corresponding location in Figure 7 for the illustration purpose.

Update of Section Stiffness and Resisting Force
The section stiffness is updated by using Equation (26) with the following updated base function: The numerator of the latter is the work done by the internal shear stress. Therefore, it satisfies the principle of virtual work. Finally, the resisting force vector has the expression:

Constitutive Relation
For the v-th iteration, the internal force vector for section i of element m can be obtained from the relation [25]: where, is the section tangent compliance from the previous step; and Δ denotes an increment operator. The tangent stiffness at the known strain e is evaluated for element m as:

Constitutive Relation
For the v-th iteration, the internal force vector for section i of element m can be obtained from the relation [25]: where, is the section tangent compliance from the previous step; and ∆ denotes an increment operator. The tangent stiffness at the known strain e is evaluated for element m as: Through numerical integration, the element flexibility matrix has the expression: where, ω i is the weights. Inverting this gives the element stiffness matrix; i.e., The section tangent stiffness can be calculated from the fiber model. The section tangent stiffness can be obtained by: where, A k is the area of the k-th fiber on the cross-section, and d mik ≡ ∂ϕ mik (ε)/∂ε is the material tangent as defined earlier.

Model Evaluation
The proposed model was evaluated on the RC beam test data reported by Bresler and Scordelis [26]. In this test, they obtained various failure modes, such as diagonal tension, shear-compression, and flexure-compression failures, by changing the shear span ratios. Another key variable was the existence of web reinforcement. Among those test specimens, the specimens OA-1 through OA-3, which did not contain stirrups, and A-1 through A-3, which contained stirrups, were selected in this presentation. The numbers 1, 2, and 3 denote the IDs that indicate the different shear span ratios, 4.0, 5.0, and 7.0, respectively. Figure 8 illustrates the sectional dimensions of the selected members, while Table 1 summarizes the test program and the results.

Model Evaluation
The proposed model was evaluated on the RC beam test data reported by Bresler and Scordelis [26]. In this test, they obtained various failure modes, such as diagonal tension, shear-compression, and flexure-compression failures, by changing the shear span ratios. Another key variable was the existence of web reinforcement. Among those test specimens, the specimens OA-1 through OA-3, which did not contain stirrups, and A-1 through A-3, which contained stirrups, were selected in this presentation. The numbers 1, 2, and 3 denote the IDs that indicate the different shear span ratios, 4.0, 5.0, and 7.0, respectively. Figure 8 illustrates the sectional dimensions of the selected members, while Table 1 summarizes the test program and the results.   Each member was analyzed by using two different elements: A Bernoulli-Euler beam element with cubic shape functions, and the proposed element. The former, driven from a stiffness-based formulation, was used to illustrate the flexural response of a beam, and the latter to include the shear effects. In the analysis model, only a half of each beam was considered, to take advantage of the symmetry. It was modeled with four elements, each of which had four Gauss-Lobatto integration points. Each member was analyzed by using two different elements: A Bernoulli-Euler beam element with cubic shape functions, and the proposed element. The former, driven from a stiffness-based formulation, was used to illustrate the flexural response of a beam, and the latter to include the shear effects. In the analysis model, only a half of each beam was considered, to take advantage of the symmetry. It was modeled with four elements, each of which had four Gauss-Lobatto integration points.
An arc-length method was used to capture a post-peak response, if any. The numerical response was terminated when numerical problems arose, usually with the satisfaction of the transverse equilibrium equations for all shear elements. For the Bernoulli-Euler beam element, the termination point was selected at a few steps after the post-peak response was obtained.  Figure 9 compares the predicted and observed curves of load versus midspan deflection relations. The prediction by using the Bernoulli-Euler beam element is denoted as predicted without shear, and the one by using the proposed element as predicted with shear. In the test results, it can be observed that when no stirrup is present, the strength drops. The strength reduction was significant, especially in the members with a small shear span ratio.  The predictions without shear present the flexural response only; therefore, there was no effect of the stirrups. This explains why the predictions without shear in series OA were similar to the counterparts in series A. As shown in the figure, the predictions without shear showed a higher member strength and stiffness and higher ductility than those with shear. The discrepancy between the two predictions, which represents the shear contribution to the deflection, increased proportionally as the applied load increased. Generally speaking, the curves of the test results The predictions without shear present the flexural response only; therefore, there was no effect of the stirrups. This explains why the predictions without shear in series OA were similar to the counterparts in series A. As shown in the figure, the predictions without shear showed a higher member strength and stiffness and higher ductility than those with shear. The discrepancy between the two predictions, which represents the shear contribution to the deflection, increased proportionally as the applied load increased. Generally speaking, the curves of the test results corresponded well to those without shear up to the point of half of the ultimate strength but became closer to them with shear afterwards. For most cases, the proposed model predicted the strength of the actual beams satisfactorily, while the model without shear consideration predicted twice as much as the actual shear strength, depending on the problem. For specimen A-1 and 2, the analysis with the proposed model was terminated shortly before it reached the actual strength. However, as can be seen in the figure, the termination points were closely related to the points where the slope of the curve of test changes abruptly.

Results
Comparison with Other Shear Models 1. Rectangular beam: The analysis results obtained for specimen OA-3 were compared with those by other shear models, such as (1) constant, or (2) parabolic shear strain model, and (3) constant shear flow model. The constant and parabolic shear strain models assumed the predefined shear strain functions over the cross-section first and then estimated the corresponding shear stresses. The advantage of the shear strain models was that the calculation procedure was simple, but those models can result in redundant constraints. On the other hand, the constant shear flow model, in which the shear force acting on the section considered is given first to obtain the shear strain of the section, can provide a more accurate shear response than the shear strain models.
The constant and parabolic shear strain distributions were obtained by using the following relations for ψ(ζ k ) in Equation (6): The constant shear flow model by using the secant stiffness method was discussed by Vecchio and Collins [6]. The shear stresses were determined by dividing a given section shear force by the height of each fiber (or layer). In the proposed model, the constant shear flow model was applied by replacing Equation (2) with the following relation: Figure 10 shows the comparison of the load-deflection relations obtained by using four different shear models. All four models showed similar initial stiffness up to the first cracking point, where the curves started to diverge. The parabolic and constant shear strain models, which had predefined shear strain distribution shape, showed a very similar response. In comparison with these models, the proposed model and constant shear flow models, which were a force-based method, showed a relatively softer behavior. This discrepancy between the constant/parabolic shear strain methods and the force-based methods comes from the fact that shear strains over the cross-section were assumed as the predefined shapes, such as a constant or parabolic, in the constant/parabolic shear strain methods, while they are not in the force-based method, because the predefined shear strains inevitably create redundant constraints in the deformability of the cross-section. In the force-based method, the shear stresses over the cross-section are calculated from equilibrium, and afterwards, the corresponding shear strains were determined, which does not cause such redundant constraints. The curve obtained by using the proposed model was closely related to the actual behavior with regards to both stiffness and strength. It is clear from the figure that the proposed model predicted the shear strength most accurately, provided that the termination point caused by the vertical equilibrium instability defined the shear failure point. 2. T-beam: The advantage that the proposed model possesses is that it describes the shear strain distributions over the cross section very accurately, and, therefore, the state of the section as well. This can best be illustrated by using T-beams, the section of which has larger shear strength than a rectangular section with the same depth and height. For this problem, the predefined shear strain methods, such as constant or parabolic shear strain method, would overestimate the shear strength by drawing more shear stress in the flange area, and keeping the web area stressed in lower intensity.
A T-beam (beam no.: T2) was selected from the extensive test conducted by Placas and Regan [27]. The sectional dimensions are: The web width is 152 mm; the total depth is 305 mm; the depth to the center of the steel is 272 mm;  Figure 11 shows the applied shear force versus the mid-span deflection curves obtained by using the four different shear models. The force-based model, including the proposed model and the constant shear flow model, terminated the analysis around the ultimate shear strength of the beam observed in the actual test. In contrast, the predefined shear strain models, including the parabolic shear strain model and the constant shear strain model, led to a shear strength more than twice as much as the actual shear strength; the analysis with the latter continued beyond the end point shown in the figure. The predefined shear strain models gave much stiffer response than the force-based models. For example, at the shear force of 54.7 kN, the constant shear strain model gave the mid-span deflection value of half as much as that obtained by the proposed model.  2. T-beam: The advantage that the proposed model possesses is that it describes the shear strain distributions over the cross section very accurately, and, therefore, the state of the section as well. This can best be illustrated by using T-beams, the section of which has larger shear strength than a rectangular section with the same depth and height. For this problem, the predefined shear strain methods, such as constant or parabolic shear strain method, would overestimate the shear strength by drawing more shear stress in the flange area, and keeping the web area stressed in lower intensity.
A T-beam (beam no.: T2) was selected from the extensive test conducted by Placas and Regan [27]. The sectional dimensions are: The web width is 152 mm; the total depth is 305 mm; the depth to the center of the steel is 272 mm; A s /b d = 1.46%; the flange breadth is 76 mm; and the flange width is 610 mm. The material properties are: The compressive concrete strength is 28.1 MPa, and the yield strength of the longitudinal bar is 620 MPa. The shear span ratio (a/d) is 3.36, and the beam was subjected to a central-point loading. The beam failed at the shear force of 54.7 kN by diagonal tension cracking. Figure 11 shows the applied shear force versus the mid-span deflection curves obtained by using the four different shear models. The force-based model, including the proposed model and the constant shear flow model, terminated the analysis around the ultimate shear strength of the beam observed in the actual test. In contrast, the predefined shear strain models, including the parabolic shear strain model and the constant shear strain model, led to a shear strength more than twice as much as the actual shear strength; the analysis with the latter continued beyond the end point shown in the figure. The predefined shear strain models gave much stiffer response than the force-based models. For example, at the shear force of 54.7 kN, the constant shear strain model gave the mid-span deflection value of half as much as that obtained by the proposed model.
To explain these discrepancies, the local behaviors were compared at the critical section (the effective depth away from the support) at the shear force of 54.7 kN, as shown in Figures 12-14. Figure 12 shows that the moment curvature values were analogous to the previous findings. The proposed model showed the largest curvature-more than twice as much as that of the constant shear strain model-and the smallest neutral surface depth. This suggests that the different shear model significantly affects the flexural behavior for the beam having a non-rectangular section. Figure 13 shows the shear strain distributions over the cross section at the same level of load. The proposed model described a parabolic shear distribution at the uncracked region (located top one-third portion), and dramatically increasing strains towards the bottom. Although a large discrepancy existed in the shear strain distributions, the contribution to the mid-span deflection of this seemed negligible. For example, the maximum shear strain estimated by the proposed model and the constant shear flow model were approximately 0.00135 and 0.00091, respectively, while the mid-span deflections were 2.54 mm and 2.79 mm, respectively. Their difference in maximum shear strain is 67% and that in mid-span deflection is 9%. Figure 14 shows the shear stress distributions. The proposed model gave a similar shape of the shear stress distributions to that of the shear strain distribution; a parabolic shape in the uncracked region, and increasing stress towards the bottom. It is interesting to see the distributions obtained from the constant shear flow were similar to those from the proposed model. Together with the previous observations, this suggests that the constant shear flow model is a good approximation for the shear behavior of a reinforced concrete beam. In contrast, with the predefined shear strain models, a large portion of the shear stress residues are expected in the flange area. This may have reduced the shear stress in the web leading to the stiffer response.
A T-beam (beam no.: T2) was selected from the extensive test conducted by Placas and Regan [27]. The sectional dimensions are: The web width is 152 mm; the total depth is 305 mm; the depth to the center of the steel is 272 mm; / ' s A b d = 1.46%; the flange breadth is 76 mm; and the flange width is 610 mm. The material properties are: The compressive concrete strength is 28.1 MPa, and the yield strength of the longitudinal bar is 620 MPa. The shear span ratio (a/d) is 3.36, and the beam was subjected to a central-point loading. The beam failed at the shear force of 54.7 kN by diagonal tension cracking. Figure 11 shows the applied shear force versus the mid-span deflection curves obtained by using the four different shear models. The force-based model, including the proposed model and the constant shear flow model, terminated the analysis around the ultimate shear strength of the beam observed in the actual test. In contrast, the predefined shear strain models, including the parabolic shear strain model and the constant shear strain model, led to a shear strength more than twice as much as the actual shear strength; the analysis with the latter continued beyond the end point shown in the figure. The predefined shear strain models gave much stiffer response than the force-based models. For example, at the shear force of 54.7 kN, the constant shear strain model gave the mid-span deflection value of half as much as that obtained by the proposed model.  To explain these discrepancies, the local behaviors were compared at the critical section (the effective depth away from the support) at the shear force of 54.7 kN, as shown in Figures 12-14. Figure 12 shows that the moment curvature values were analogous to the previous findings. The proposed model showed the largest curvature-more than twice as much as that of the constant shear strain model-and the smallest neutral surface depth. This suggests that the different shear model significantly affects the flexural behavior for the beam having a non-rectangular section.   Figure 13 shows the shear strain distributions over the cross section at the same level of load. The proposed model described a parabolic shear distribution at the uncracked region (located top one-third portion), and dramatically increasing strains towards the bottom. Although a large discrepancy existed in the shear strain distributions, the contribution to the mid-span deflection of this seemed negligible. For example, the maximum shear strain estimated by the proposed model and the constant shear flow model were approximately 0.00135 and 0.00091, respectively, while the mid-span deflections were 2.54 mm and 2.79 mm, respectively. Their difference in maximum shear strain is 67% and that in mid-span deflection is 9%. one-third portion), and dramatically increasing strains towards the bottom. Although a large discrepancy existed in the shear strain distributions, the contribution to the mid-span deflection of this seemed negligible. For example, the maximum shear strain estimated by the proposed model and the constant shear flow model were approximately 0.00135 and 0.00091, respectively, while the mid-span deflections were 2.54 mm and 2.79 mm, respectively. Their difference in maximum shear strain is 67% and that in mid-span deflection is 9%.  Figure 14 shows the shear stress distributions. The proposed model gave a similar shape of the shear stress distributions to that of the shear strain distribution; a parabolic shape in the uncracked region, and increasing stress towards the bottom. It is interesting to see the distributions obtained from the constant shear flow were similar to those from the proposed model. Together with the previous observations, this suggests that the constant shear flow model is a good approximation for the shear behavior of a reinforced concrete beam. In contrast, with the predefined shear strain models, a large portion of the shear stress residues are expected in the flange area. This may have reduced the shear stress in the web leading to the stiffer response.

Conclusions
In this study, the element formulation with the stress-driven shear force functions was derived and realized in the context of the finite element method. The proposed method showed excellent agreement with various RC beam test results. The two-dimensional representation of beam shear on the micro-plane successfully described the sectional shear stresses and strains, hence, allowing shear degradation and shear failure prediction in the global response. In addition, the coupled bending and shear behavior was accurately described, and the authors believe that in future studies, the axial load effects on shear behavior should be easily demonstrated, without any modification of the element formulation. In addition, since the formulation proposed in this study is very consistent with a strong applicability for nonlinear finite element analysis, it is expected that the proposed method can be applied to the existing finite element analysis software without interchangeability problems.
However, it is not clear at this point whether the incomplete post-peak response resulted from the inherent inconsistency in the element formulation. The inherent inconsistency may be ascribed to the use of the plane section assumption along with the stress-driven functions. Theoretically, the error due to this should disappear at convergence.

Conclusions
In this study, the element formulation with the stress-driven shear force functions was derived and realized in the context of the finite element method. The proposed method showed excellent agreement with various RC beam test results. The two-dimensional representation of beam shear on the micro-plane successfully described the sectional shear stresses and strains, hence, allowing shear degradation and shear failure prediction in the global response. In addition, the coupled bending and shear behavior was accurately described, and the authors believe that in future studies, the axial load effects on shear behavior should be easily demonstrated, without any modification of the element formulation. In addition, since the formulation proposed in this study is very consistent with a strong applicability for nonlinear finite element analysis, it is expected that the proposed method can be applied to the existing finite element analysis software without interchangeability problems.
However, it is not clear at this point whether the incomplete post-peak response resulted from the inherent inconsistency in the element formulation. The inherent inconsistency may be ascribed to the use of the plane section assumption along with the stress-driven functions. Theoretically, the error due to this should disappear at convergence.