# Nonlinear Finite Element Analysis Formulation for Shear in Reinforced Concrete Beams
## 3 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:  
\[\varphi(\zeta)=\int_{\zeta}^{\frac{h}{2}}df_{x}(\zeta)b_{w}(\zeta)d\zeta \tag{2}\]  
where, \(-\frac{h}{2}\leq\zeta<\frac{h}{2}\); \(df_{x}(\zeta)\) denotes the difference of longitudinal stresses between cross-sections \(mn\) and \(m_{1}n_{1}\), \(b_{w}(\zeta)\) is the width of the layer at the depth of \(\zeta\), and \(h\) is the overall height of the section. The horizontal shear stress at \(\zeta\) is obtained by dividing the shear flow by the width of the layer at \(\zeta\); i.e.,  
\[v_{xy}(\zeta)=\frac{\varphi(\zeta)}{b_{w}(\zeta)} \tag{3}\]  
Figure 1: Timoshenko beam.  
The normal stress gradient or the stress differences between the two sections \(df_{x}\) can be defined as:  
\[df_{x}(\zeta)=[1\quad 0]\mathbf{d}d\varepsilon \tag{4}\]  
where, the first term [10] is multiplied to obtain only the scalar normal stress, \(\mathbf{d}\) is the \(2\times 2\) sectional stiffness matrix of the cross section, which will be defined later, and \(d\varepsilon\) is the \(2\times 1\) sectional strain vector. The last term is the sectional strain changes due to the moment change \(\Delta M\) as shown in Figure 1(a), and related to the \(3\times 1\) full strains as,  
\[d\varepsilon=\varepsilon d\mathbf{e} \tag{5}\]  
And  
\[\varsigma\equiv\begin{bmatrix}1&-\zeta&0\\ 0&0&\psi(\zeta)\end{bmatrix} \tag{6}\]  
where, \(\psi\) is the shear strain function. Depending on how this is defined results in a constant, parabolic, or other shear distribution. For example, if \(\psi(\zeta_{k})=1\), a constant shear distribution is obtained.  
Knowing that Equation (5) defines the sectional strain changes due to the moment change \(\Delta M=V\) leads to:  
\[d\mathbf{e}=[K]^{-1}\begin{bmatrix}0\\ V\\ 0\end{bmatrix} \tag{7}\]  
Herein, \(\mathbf{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:  
\[v_{xy}(\zeta_{k})=\varphi(\zeta_{k})/b_{w}(\zeta_{k}) \tag{8}\]  
Figure 2: Equilibrium in a beam subjected to non-uniform bending moment.  
where, \(\varphi(\zeta_{k})=\sum\limits_{j=1}^{k}h_{f}\big{(}\zeta_{j}\big{)}\cdot df_{x} \big{(}\zeta_{j}\big{)}\cdot b_{w}\big{(}\zeta_{j}\big{)}\), and \(h_{f}\big{(}\zeta_{j}\big{)}\) is the height of the \(j\)-th fiber. Finally, the shear stress function at the \(k\)-th fiber is obtained from:  
\[\overline{\omega}(\zeta_{k})^{v}=v_{xy}(\zeta_{k})/V_{mi} \tag{9}\]  
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.