An Improved Formulation and Analysis of Reddy Beam Model for Framed Structures

Abstract

A structural analysis of framed structures using the finite element method considering both the Bernoulli- Euler and the Timoshenko beam theories can be performed adopting cubic interpolation functions that yield analytical solutions for the displacements. However, these theories may not provide stress results with sufficient accuracy. In such cases, it is necessary to employ higher order beam formulations, that may require a high level of discretization. Therefore, this study proposes an enhanced Reddy beam element, obtained by considering interpolation functions calculated directly from the solution of the differential equation system. This solution minimizes the impact of structural discretization on the analysis, and framed structures can be effectively modeled considering the minimum number of elements required to describe the geometry. The results obtained by the proposed formulation were compared against classical beam theories and the Reddy beam model adopting conventional shape functions, showing the efficacy of the proposed element in simulating the elastic behavior of framed structures in a FEM-like procedure.

Keywords: Reddy stiffness matrix; Analytical interpolation functions; Reduced discretization; High-order beams

Graphical Abstract

1 INTRODUCTION

The quality of information obtained from a numerical analysis of framed structures depends on the structural theory considered. Among these, the most classical theory is the Bernoulli- Euler beam theory (EBBT). EBBT assumes that beam cross-sections remain plane and orthogonal to the beam axis after deformation, thereby disregarding shear deformation. As a consequence, this formulation proves inadequate to predict the behavior of structures with a moderate slenderness ratio or with a small shear-to-bending ratio (Rodrigues et al. 2021a).

As an alternative, the Timoshenko bending theory (TBT) considers that cross sections remain plane, but not necessarily orthogonal to the deformed longitudinal axis. It introduces shear deformation as an additional rotation of the cross-section (Timoshenko 1921), resulting in an element where the transversal displacement and rotation are independent variables. This yields accurate structural displacements, regardless of the slenderness of the beam.

While both beams theories provide a reliable and often comparable displacement field, they exhibit different results regarding shear stress, for both not consistent with the correct shear stress distribution in the section, given that the shear stress is undetermined when considering the EBBT, and the TBT leads to a constant shear stress distribution (Polizzotto 2015) dependent of the shear correction factors (Dong et al. 2010), in contrast to the requirement that the upper and lower surfaces of the beam be stress free.

Usually, the dominant criterion in the design of a beam regarding its resistance is the maximum value of normal stress, related to the bending moment. However, shear stresses can be important, particularly in thin-walled bars and short, thick beams. According with Chakrabarti et al. (2012), despite of the Bernoulli-Euler beam theory does not consider the transverse shear deformation, its effect is very important for beams having low values of span to depth ratio, low shear rigidity or continuous spans. In Timoshenko’s beam theory (TBT), transverse shear strain distribution is assumed to be constant through the beam thickness and thus requires problem dependent shear correction factor, whereas the actual variation of shear stress is parabolic becoming zero at the beam top and bottom surfaces. While the shear correction factor helps to predict the global response such as deflection or vibration frequency well, it is not sufficient for an accurate prediction of the local response such as the shear stress distributions within structures like composite beams (Uddin et al. 2017a).

To improve shear stress prediction, several high-order bending theories were developed in the literature (Levinson 1981; Bickford 1982; Reddy 1984a; 1984b; Heyliger and Reddy 1988; Petrolito 1995), where cross sections deform according to polynomial, trigonometric, or hyperbolic functions (Shi and Voyiadjis 2011; Almeida et al. 2011; Nguyen et al. 2022; Vinh and Son 2022; Neves et al. 2011; Vidal and Polit 2008; El Meiche et al. 2011; Mahi et al. 2015).

High-order beam models typically have the capability to accurately predict shear stress distribution (Reddy 2022), making them widely applicable in structural engineering to simulate the behavior of steel and steel-concrete composite elements (Chakrabarti et al. 2012; Meghare and Jadhao 2015; Uddin et al. 2017a; 2017b; 2018; 2020). They are also relevant in structural analyses of laminated composite and sandwich beams, whether associated or not with refined theories such as the zig zag theory (Karama et al. 2003; Vidal and Polit 2011; Vo and Thai 2012; Szenkrényes 2014; Iurlaro et al. 2015; Karttunen et al. 2017; Zhen et al. 2019; Leite and Rocha 2023). Sayyad and Ghugal (2017) provide a comprehensive literature review on this topic.

These theories are also employed in the analysis of micro-beams and nonlocal theories (Reddy 2007; Wanji et al. 2012; Mohammad-Abadi and Daneshmehr 2014; Sahmani et al. 2014; Nazemnezhad and Zare 2016; Augello et al. 2019; Faroughi et al. 2020; Golbakhshi et al. 2022; Karamanli et al. 2023; Zheng et al. 2023), nonhomogeneous and functionally graded beams (Simsek and Reddy 2013; Khorshidi et al. 2016; Darijani and Mohammadabadi 2014; Karamanli and Vo 2020; Ruocco et al. 2020; Ruocco and Reddy 2021). Higher order beam have also shown satisfactory results in addressing buckling and vibration of beams and nanobeams (He and Yang 2014; Nazemnezhad and Zare 2016; Mittelstedt and Mittelstedt 2020; AkhavanAlavi et al. 2019; Lin et al. 2020; Muc 2020) as well as in analysis of beams resting on elastic foundations (Stojanović and Petković 2016), and porous beams (Keleshteri and Jelovica 2022).

The Reddy beam model (Reddy 1984b; 2022; Heyliger and Reddy 1988) represents the most classical theory among these higher-order beam models. In this model, a modification of the Timoshenko beam model is proposed, assuming that the cross section is capable of warping according to a suitably specified warping mode, with the profile of the warped cross section crosses the upper part and lower surfaces orthogonally (Polizzotto 2015). This effect is achieved by expanding the displacement field as cubic functions of the coordinate along the thickness of the cross section (Ruocco and Reddy 2023). In Reddy et al. (1997) is presented relationships between bending solutions of classical and shear deformation beam theories.

Nevertheless, the augmentation of the stress field is not free of challenges. In contrast to more straightforward beam models, higher-order beams typically lack analytical solutions and necessitate a numerical approach for resolution.

In a FEM context, the application of interpolation functions calculated from the solution of the differential equations governing the problem at hand - for both Bernoulli- Euler and Timoshenko beam theories cubic Hermitian interpolation functions (Rodrigues et al. 2019; Molina-Villegas and Ortega 2023; Reddy 1997) - leads to numerical results coincident to the analytical solution and, as a consequence, a structural response independent on the level of discretization (Rodrigues et al. 2021a). Structural analysis software usually employs these approaches (McGuire et al. 2000; Martha and Rangel 2022; Rangel and Martha 2019), returning the exact solution of structures with complex geometries using only one element per member (Yang and Leu 1994; Yang and Kuo 1994; Bathe 1996; Reddy 2019).

However, for the Reddy beam model, this is not the case. Adopting, as usual in the literature, Hermitian cubic functions to represent the out-of-plane displacements and linear functions to represent both axial displacements and rotations, even for a simple problem such as a clamped beam, a high level of discretization is necessary to converge to the exact response of the problem. Notably, in this manner the rotational components are governed by a function independent of transversal displacements.

Recently an analytical solution of the equations governing the Reddy beam has been presented in (Ruocco and Reddy 2023), proving effective in both linear and nonlinear scenarios. Nevertheless, this solution incorporates exponential terms that compromise its efficiency within the framework of a finite element numerical process (Rodrigues et al. 2021b). These terms still require a dense discretization to achieve solution convergence, thereby diminishing the advantages derived from the use of an analytical solution as shape functions.

In this paper, a new shape function and its corresponding stiffness matrix are introduced to better approximate the solution of a Reddy beam. This is achieved by developing a series expansion of the exponential terms from the exact solution, thereby regularizing their behavior and obtaining a polynomial kinematics that can be easily implemented in an automated computational code. The derivation from the analytical solution, along with the coupling between displacements and rotations stemming from it, yields a solution that is simultaneously stable and converging to the exact solution with a significantly reduced discretization if compared to conventional shape functions. This reduction in discretization does not compromise the quality of the solution, both in terms of displacement and the associated stress field, as demonstrated by the examples provided.

2 REDDY BEAM MODEL

In the Reddy beam theory, as described by Polizzotto (2015), the Timoshenko model is modified to permit the warping of cross-sections. This modification ensures that the profile of the warped section intersects the upper and lower surfaces orthogonally, as illustrated in Figure 1.

Figure 1
Bernoulli- Euler, Timoshenko and Reddy bending kinematics - Adapted fromRodrigues et al. (2021a).

This effect is achieved by proposing a cubic displacement field along the thickness of the cross section (Ruocco and Reddy, 2023). Assuming zero shear stress on the upper and lower faces of the section, the displacement field:

$u\left(x,y\right)=u_0\left(x\right)-y\theta \left(x\right)+\alpha y^3\left(\theta \left(x\right)-\frac{d{v}_0\left(x\right)}{dx}\right)$

$v\left(x,y\right)=v_0\left(x\right)$ (1)

is finally obtained. In equation (1), $u_0\left(x\right)$ and $v_0\left(x\right)$ represent the axial and transverse displacement, respectively, of the section centroidal axis, $\theta \left(x\right)$ indicates the rotation of the cross section, allowing to obtain the longitudinal and transversal displacements $u\left(x,y\right)$ and $v\left(x,y\right)$ of any point in the section, as illustrated in Figure 2. For rectangular sections, the Reddy constant $\alpha =4/3h^2$ is given in terms of the cross section height $h$.

Figure 2
Beam displacement field - Adapted from Rodrigues et al. (2019) .

The corresponding non-zero components of the linear part of the Green-Lagrange strain tensor are:

${\epsilon }_x=\frac{\partial u\left(x,y\right)}{\partial x}=\frac{d{u}_0}{dx}-y\frac{d\theta }{dx}+\alpha y^3\left(\frac{d\theta \left(x\right)}{dx}-\frac{d^2{v}_0\left(x\right)}{d{x}^2}\right)$ (2)

${\gamma }_{xy}=\frac{\partial v\left(x,y\right)}{\partial x}+\frac{\partial u\left(x,y\right)}{\partial y}=\left(\theta -\frac{d{v}_0}{dx}\right)\left(3\alpha y^2-1\right)$ (3)

The second Piola-Kirchhoff stress tensor and the Green-Lagrange strain tensor are energetically conjugated and satisfy the condition of being expressible in terms of the know configuration (McGuire et al. 2000). The stresses are obtained from the constitutive law of the material, for elastic behavior, and assuming homogenous and isotropic material, the stress-strain tensor of the know configuration is constant and leads to the stress:

${\sigma }_x=E{\epsilon }_x\tau =G\gamma$ (4)

This development is consistent with the nonlinear kinematics of the cross-sectional Reddy beam (Reddy 1984a; 1984b; Heyliger and Reddy 1988, Reddy 2022; Ruocco and Reddy 2023). Thus, indicating with $q$ the transversal load, and $m$ the distributed moment along the element $x$ axis, the virtual work principle can be written as:

${\int }_V^{}{\sigma }_x\delta \left[-y\frac{d\theta }{dx}+\alpha y^3\left(\frac{d\theta }{dx}-\frac{d^2{v}_0}{d{x}^2}\right)\right]dV+{\int }_V^{}\tau \delta \left[\left(\theta -\frac{d{v}_0}{dx}\right)\left(3\alpha y^2-1\right)\right]dV-{\int }_0^L\left(q\delta v_0+m\delta \theta \right)dx=0$ (5)

that is,

${\int }_0^L\left(M-\alpha P\right)\delta \frac{d\theta }{dx}dx+{\int }_0^L\alpha P\delta \frac{d^2{v}_0}{d{x}^2}dx+{\int }_0^L\left(Q-3\alpha V\right)\left(\delta \theta -\delta \frac{d{v}_0}{dx}\right)dx-{\int }_0^L\left(q\delta v_0+m\delta \theta \right)dx=0$ (6)

having indicate with:

$Q=-{\int }_A^{}\tau dA;\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}M=-{\int }_A^{}y{\sigma }_{x}dA;\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}V=-{\int }_A^{}{y}^2\tau dA;\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}P=-{\int }_A^{}{y}^3{\sigma }_{x}dA\mathrm{}$ (7)

the generalized stress components.

Integrating equation (6) by parts leads to the system of differential equations governing the Reddy beam model:

$\delta v_0:\mathrm{}\frac{d^2\left(\alpha P\right)}{d{x}^2}+\frac{d\left(Q-3\alpha V\right)}{dx}+q=0\to \frac{dT}{dx}+q=0$

$\delta \theta :\mathrm{}\left(Q-3\alpha V\right)-\frac{d\left(M-\alpha P\right)}{dx}+m=0\to T-\frac{dM}{dx}+m=0$ (8)

wherein:

$T=\left(Q-3\alpha V\right)+\frac{d\left(\alpha P\right)}{dx}\mathrm{}$ (9)

is the effective shear force.

The equation system (8) can be rewritten in terms of ($v_0$, $\theta$), returning:

$\left\{G\left(A-3\alpha I\right)\left(\frac{d^2{v}_0}{d{x}^2}-\frac{d\theta }{dx}\right)-3\alpha \left[GI\left(\frac{d^2{v}_0}{d{x}^2}-\frac{d\theta }{dx}\right)-3\alpha GF\left(\frac{d^2{v}_0}{d{x}^2}-\frac{d\theta }{dx}\right)\right]\right\}$

$+E\alpha F\frac{d^3\theta }{d{x}^3}-E{\alpha }^{2}H\left(\frac{d^3\theta }{d{x}^3}-\frac{d^4{v}_0}{d{x}^4}\right)=q$ (10)

$\left\{G\left(A-3\alpha I\right)\left(\frac{d{v}_0}{dx}-\theta \right)-3\alpha \left[GI\left(\frac{d{v}_0}{dx}-\theta \right)-3\alpha GF\left(\frac{d{v}_0}{dx}-\theta -\right)\right]\right\}$

$+E\alpha F\frac{d^2\theta }{d{x}^2}-E{\alpha }^{2}H\left(\frac{d^2\theta }{d{x}^2}-\frac{d^3{v}_0}{d{x}^3}\right)-\left[EI\frac{d^2\theta }{d{x}^2}-E\alpha F\left(\frac{d^2\theta }{d{x}^2}-\frac{d^3{v}_0\left(x\right)}{d{x}^3}\right)\right]=m$ (11)

where, $\left(I,F,H\right)={\int }_A^{}\left(y^2,y^4,y^6\right)dA$ correspond to the higher-order moment of inertia.

3 SOLUTIONS OF THE DIFFERENTIAL EQUATIONS

The Reddy beam theory’s differential equations can be rewritten into a homogeneous form by setting equations (10) and (11) to zero. Upon derivation of equation (11), and substituting expression (10), the resulting differential relation is given by:

$\left[EI\frac{d^3\theta }{d{x}^3}-E\alpha F\left(\frac{d^3\theta }{d{x}^3}-\frac{d^4{v}_0\left(x\right)}{d{x}^4}\right)\right]=0\mathrm{}$ (12)

The third derivative of the rotation is then calculated as:

$\frac{d^3\theta }{d{x}^3}=\frac{-E\alpha F\left(\frac{d^4{v}_0\left(x\right)}{d{x}^4}\right)}{\left(EI-E\alpha F\right)}$ (13)

integrating equation (12) twice yields expressions for the cross-sectional rotation and its derivative:

$\frac{d\theta }{dx}=\frac{-E\alpha F\left(\frac{d^2{v}_0\left(x\right)}{d{x}^2}\right)-\mathrm{}{c}_{5}x-c_6}{\left(EI-E.\alpha .F\right)}$ (14)

$\theta =\frac{-E\alpha F\left(\frac{d{v}_0\left(x\right)}{dx}\right)-\mathrm{}{c}_5\frac{x^2}{2}-c_{6}x}{\left(EI-E\alpha F\right)}+c_7$ (15)

Applying relations (13) to (15) to the equilibrium equation (10) results in a fourth-order differential equation in terms of the out-of-plane displacement$v_0$:

$\frac{EI}{105}\frac{d^4{v}_0}{d{x}^4}-\frac{4EI\frac{d^2{v}_0}{d{x}^2}+4\mathrm{}{c}_{5}x+c_6}{h^2\left(1+\nu \right)}=0$ (16)

returning:

$v_0=-\frac{\mathrm{}{c}_1{h}^2\left(1+\nu \right)}{420EI}x+c_2+c_3{e}^{-\mu x}+c_4{e}^{\mu x}-\frac{\mathrm{}{c}_5{x}^3}{6EI}-\frac{\mathrm{}{c}_6{x}^2}{2EI}\mathrm{}$ (17)

where $\mu =\frac{2\sqrt{105}}{h\sqrt{1+\nu }}$ is a constant depending on the cross-section height $h$ and on the Poisson coefficient $\nu$, and $c_i$ are constants of integration to be determined by posing the proper boundary conditions.

The cross-sectional rotation, determined from equation (15), takes the form:

$\theta =\frac{\mathrm{}{c}_1{h}^2\left(1+\nu \right)}{168EI}+\frac{\mu }{4}{c}_{3}3e^{-\mu x}-\frac{\mu }{4}{c}_4{e}^{\mu x}-\frac{\mathrm{}{c}_5{x}^2}{2EI}-\frac{\mathrm{}{c}_{6}x}{EI}\mathrm{}$ (18)

According with Rodrigues et al. (2019; 2021a), this homogeneous solution $v_h\left(x\right)$ corresponds to the final solution in the direct stiffness method when no transverse force q is present within the element.

Given that the interpolation functions are derived from the homogenous solution of the differential equation system ensures that the displacements and nodal rotations in the discrete model are exact, independent of the level of the discretization.

4 INTERPOLATION FUNCTIONS

The elastic behavior of a Reddy beam can be represented as a function of its nodal displacements (Figure 3) using the approximation provided by equations (19) to (22).

Figure 3
Deformed configuration of an isolated Reddy element - Adapted from Rodrigues et al. (2021a) .

$u_0\left(x\right)=N_1^u\left(x\right)d_1^{\mathrm{\text{'}}}+N_4^u\left(x\right)d_4^{\mathrm{\text{'}}}\to u_0\left(x\right)=\left\{N_u\left(x\right)\right\}\left\{u\right\}$ (19)

$v_0\left(x\right)=N_2^v\left(x\right)d_2^{\mathrm{\text{'}}}+N_3^v\left(x\right)d_{3dv}^{\mathrm{\text{'}}}+N_5^v\left(x\right)d_5^{\mathrm{\text{'}}}+N_6^v\left(x\right)d_{6v}^{\mathrm{\text{'}}}\to v_0\left(x\right)=\left\{N_v\left(x\right)\right\}\left\{v\right\}$ (20)

$\frac{{dv}_0\left(x\right)}{dx}=\frac{N_2^v\left(x\right)}{dx}{d}_2^{\mathrm{\text{'}}}+\frac{N_3^v\left(x\right)}{dx}{d}_{3dv}^{\mathrm{\text{'}}}+\frac{N_5^v\left(x\right)}{dx}{d}_5^{\mathrm{\text{'}}}+\frac{N_6^v\left(x\right)}{dx}{d}_{6dv}^{\mathrm{\text{'}}}\to \frac{{dv}_0\left(x\right)}{dx}=\left\{\frac{{dN}_v\left(x\right)}{dx}\right\}\left\{v\right\}$ (21)

$\theta \left(x\right)=N_1^{\theta }\left(x\right)d_{3\theta }^{\mathrm{\text{'}}}+N_4^{\theta }\left(x\right)d_{6\theta }^{\mathrm{\text{'}}}\to \mathrm{}\theta \left(x\right)=\left\{N_{\theta }\left(x\right)\right\}\left\{v\right\}$ (22)

Typically, numerical analyses are performed by employing Hermitian interpolating functions for the out-of-plane component $v_0\left(x\right)$, and linear interpolation for both the axial displacement $u_0\left(x\right)$ and rotation $\theta \left(x\right)$.

However this study, making use of the analytical solutions presented in equations (17) and (18), proposes improved interpolation functions to better represents the out-of-plane displacement field. Rewriting the obtained displacement field as:

$\left\{\begin{array}{c}{v}_0\left(x\right)\\ \frac{{dv}_0\left(x\right)}{dx}\\ \theta \left(x\right)\end{array}\right\}=\left[\begin{array}{c}\begin{array}{cc}-\frac{\mathrm{}{h}^2\left(1+\nu \right)x}{420EI}& \mathrm{}\mathrm{}\begin{array}{cc}1& \begin{array}{cc}\mathrm{}{\mathrm{}\mathrm{}\mathrm{}e}^{-\mu x}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{cc}{e}^{\mu x}& \begin{array}{cc}-\frac{\mathrm{}{x}^3}{6EI}& -\frac{\mathrm{}{x}^2}{2EI}\end{array}\end{array}\end{array}\end{array}\end{array}\\ \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}-\frac{\mathrm{}{h}^2\left(1+\nu \right)}{420EI}& \begin{array}{cc}0& \begin{array}{cc}-\mu e^{-\mu x}& \begin{array}{cc}\mu e^{\mu x}& \begin{array}{cc}-\frac{x^2}{2EI}& -\frac{x}{EI}\end{array}\end{array}\end{array}\end{array}\end{array}\\ \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{cc}\frac{\mathrm{}{h}^2\left(1+\nu \right)}{168EI}& \begin{array}{cc}0& \begin{array}{cc}\frac{\mu e^{-\mu x}}{4}& \begin{array}{cc}\frac{-\mu e^{\mu x}}{4}& \begin{array}{cc}-\frac{x^2}{2EI}& -\frac{x}{EI}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]\mathrm{}\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{c}_1\\ \begin{array}{c}{c}_2\\ \begin{array}{c}{c}_3\\ c_4\end{array}\end{array}\end{array}\\ c_5\end{array}\\ c_6\end{array}\right\}=\left[X\right]\left\{C\right\}$ (23)

it is possible to evaluate the nodal values represented in Figure 4 as:

Figure 4
Clamped beam with concentrated load.

$\left\{d^{\mathrm{\text{'}}}\right\}=\left\{\begin{array}{c}{d}_2^{\mathrm{\text{'}}}\\ d_{3dv}^{\mathrm{\text{'}}}\\ \begin{array}{c}{d}_{3\theta }^{\mathrm{\text{'}}}\\ \begin{array}{c}{d}_5^{\mathrm{\text{'}}}\\ \begin{array}{c}{d}_{6v}^{\mathrm{\text{'}}}\\ d_{6\theta }^{\mathrm{\text{'}}}\end{array}\end{array}\end{array}\end{array}\right\}=\left\{\begin{array}{c}{v}_0\left(0\right)\\ \frac{{dv}_0\left(0\right)}{dx}\\ \begin{array}{c}\theta \left(0\right)\\ \begin{array}{c}{v}_0\left(L\right)\\ \begin{array}{c}\frac{{dv}_0\left(L\right)}{dx}\\ \theta \left(L\right)\end{array}\end{array}\end{array}\end{array}\right\}=\left[\begin{array}{c}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}1\end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}1\end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}1\end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}\\ \begin{array}{c}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\mathrm{}-\frac{h^2(1+\nu )}{420EI}& 0\end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}-\mu \end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mu \end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}\\ \mathrm{}\mathrm{}\mathrm{}\begin{array}{c}\mathrm{}\mathrm{}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\frac{h^2(1+\nu )}{1680EI}& \mathrm{}0\end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}-\mu \end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mu \end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}\\ \begin{array}{c}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}-\frac{L{h}^2(1+\nu )}{420EI}& 1\end{array}& {\mathrm{}\mathrm{}\mathrm{}e}^{-L\mu }\end{array}& {\mathrm{}\mathrm{}\mathrm{}e}^{L\mu }\end{array}& \mathrm{}\mathrm{}\frac{-L^3}{6EI}\end{array}& \mathrm{}\mathrm{}\frac{-L^2}{2EI}\end{array}\\ \begin{array}{c}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}-\frac{h^2(1+\nu )}{420EI}& 0\end{array}& -\mu e^{-L\mu }\end{array}& \mu e^{L\mu }\end{array}& \frac{-L^3}{6EI}\end{array}& \frac{-L^2}{2EI}\end{array}\\ \begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\frac{h^2(1+\nu )}{1680EI}& \mathrm{}\mathrm{}0\end{array}& \mathrm{}\frac{\mu e^{-L\mu }}{4}\end{array}& \frac{-\mu e^{L\mu }}{4}\end{array}& \frac{-L^2}{2EI}\end{array}& \frac{-L}{EI}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}\begin{array}{c}{c}_1\\ \begin{array}{c}{c}_2\\ \begin{array}{c}{c}_3\\ c_4\end{array}\end{array}\end{array}\\ c_5\end{array}\\ c_6\end{array}\right\}\to \left\{d^{\mathrm{\text{'}}}\right\}=\left[H\right]\left\{C\right\}$ (24)

the interpolation functions can be then obtained as:

${\left\{\begin{array}{ccc}{v}_0\left(x\right)& \frac{{dv}_0\left(x\right)}{dx}& \theta \left(x\right)\end{array}\right\}}^T=\left[N\right]\left\{d^{\mathrm{\text{'}}}\right\}$ (25)

having indicate with:

$\left[N\right]=\left[X\right]{\left[H\right]}^{-1}=\left[\begin{array}{c}\mathrm{}\begin{array}{cc}{\mathrm{}N}_2^v\left(x\right)& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{cc}{N}_3^v\left(x\right)& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0& \begin{array}{cc}{\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}N}_5^v\left(x\right)& \begin{array}{cc}{\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}N}_6^v\left(x\right)\mathrm{}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}\end{array}\end{array}\end{array}\end{array}\\ \begin{array}{c}\begin{array}{cc}\frac{d{N}_2^v\left(x\right)}{dx}& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\frac{d{N}_3^v\left(x\right)}{dx}& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{cc}\frac{d{N}_5^v\left(x\right)}{dx}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{cc}\frac{d{N}_6^v\left(x\right)}{dx}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}\end{array}\end{array}\end{array}\end{array}\\ \begin{array}{cc}{\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}N}_2^{\theta }\left(x\right)& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{cc}{N}_3^{\theta }\left(x\right)& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{N}_1^{\theta }\left(x\right)& \begin{array}{cc}{\mathrm{}\mathrm{}N}_5^{\theta }\left(x\right)& \begin{array}{cc}{\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}N}_6^{\theta }\left(x\right)\mathrm{}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}{N}_4^{\theta }\left(x\right)\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]$ (26)

the interpolation functions.

From equation (23), it becomes evident that the determined shape functions contain exponential terms. The presence of these terms introduces numerical instabilities, offsetting the advantages derived from using an analytical solution. Unlike beam models characterized by polynomial solutions, exponential terms necessitate a dense discretization to converge to a stable solution. To eliminate these numerical instabilities, it is possible to replace the exponential terms with their expansion truncated to the second order, resulting in:

$\left[N\right]=\left[\begin{array}{c}\begin{array}{c}\begin{array}{cc}2\frac{x^3}{L^3}-3\frac{x^2}{L^2}+1\mathrm{}\mathrm{}\mathrm{}& \begin{array}{cc}\frac{x^3}{L^2}-2\frac{x^2}{L}+x& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}3\frac{x^2}{L^2}-2\frac{x^3}{L^3}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\begin{array}{cc}\frac{x^3}{L^2}-\frac{x^2}{L}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}0\end{array}\end{array}\end{array}\end{array}\mathrm{}\end{array}\\ \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}6\frac{x^2}{L^3}-6\frac{x}{L^2}& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}3\frac{x^2}{L^2}-4\frac{x}{L}+1& \mathrm{}\mathrm{}\begin{array}{cc}\mathrm{}\mathrm{}0& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}6\frac{x}{L^2}-6\frac{x^2}{L^3}& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}3\frac{x^2}{L^2}-2\frac{x}{L}& \mathrm{}\mathrm{}0\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\\ \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}6\frac{x^2}{L^3}-6\frac{x}{L^2}& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}3\frac{x^2}{L^2}-3\frac{x}{L}& \begin{array}{cc}\mathrm{}\mathrm{}\mathrm{}1-\frac{x}{L}& \mathrm{}\begin{array}{cc}6\frac{x}{L^2}-6\frac{x^2}{L^3}& \begin{array}{cc}\mathrm{}\mathrm{}3\frac{x^2}{L^2}-3\frac{x}{L}& \mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\frac{x}{L}\end{array}\end{array}\end{array}\end{array}\end{array}\end{array}\right]$ (27)

Finally, the total cross-sectional becomes:

$\theta \left(x\right)=N_2^{\theta }\left(x\right)d_2^{\mathrm{\text{'}}}+N_3^{\theta }\left(x\right)d_{3dv}^{\mathrm{\text{'}}}+N_1^{\theta }\left(x\right)d_{3\theta }^{\mathrm{\text{'}}}+N_5^{\theta }\left(x\right)d_5^{\mathrm{\text{'}}}+N_6^{\theta }\left(x\right)d_{6v}^{\mathrm{\text{'}}}+N_4^{\theta }\left(x\right)d_{6\theta }^{\mathrm{\text{'}}}$ (28)

To show the potentiality of the proposed interpolation functions, the problem of a clamped beam was analyzed (Figure 4). The proposed beam have a length of L = 1 m, Young’s modulus of E = 10⁷ kN/m², Poisson’s ratio $\nu =0.3$.

According with Ruocco and Reddy (2023), the analytical solution of the problem considering the Reddy beam model is given by:

$v_0\left(x\right)=\frac{P\left(3L-x\right)x^2}{6EI}+\frac{12P(1+\nu )x}{5EA}+\frac{12P(1+\nu )sech\left(\mu L\right)\left[sinh\left(\mu \left(L-x\right)\right)-sinh\left(\mu L\right)\right]}{5EA\mu }$ (29)

$\frac{{dv}_0\left(x\right)}{dx}=-\frac{P\left(2L-x\right)x}{2EI}+\frac{12P\left(1+\nu \right)}{5EA}-\frac{12P(1+\nu )sech\left(\mu L\right)cosh\left(\mu \left(L-x\right)\right)}{5EA}$ (30)

$\theta \left(x\right)=\frac{P\left(2L-x\right)x}{2EI}+\frac{3P\left(1+\nu \right)}{5EA}-\frac{3P(1+\nu )sech\left(\mu L\right)cosh\left(\mu \left(L-x\right)\right)}{5EA}$ (31)

Figures 5 to 7 shows the displacements obtained by comparing the analytical solution with numerical results, denoted as RBT, where Hermitian and linear shape functions are adopted for displacement and rotations, respectively. Additionally, are reported the results obtained with the proposed modified shape functions, and indicated with MRBT, considering a discretization with a single element. The solutions are presented for two different slenderness ratio $\lambda$($L/h$) = 10, and $\lambda =40$. The concentrated load is P = 1030 kN for the slenderness ratio of 10, and 16 kN for the slenderness ratio of 40. Both solutions, analytical and numerical are performed considering just a linear elastic analysis for comparison purposes only.

Figure 5
Displacement ( $v$ ) comparison between solutions - (a) $\lambda =10$ and (b) $\lambda =40$ .

Figure 6
Displacement ( $dv/dx$ ) comparison between solutions - (a) $\lambda =10$ and (b) $\lambda =40$ .

Figure 7
Cross-sectional rotation ( $\theta$ ) comparison between solutions - (a) $\lambda =10$ and (b) $\lambda =40$ .

The results demonstrate the effectiveness of the proposed interpolation in predicting the displacement field of the beam. Moreover, for the highest slenderness ratio $\lambda =40$, the analytical solution faces challenges in providing a solution for the transversal displacement due to numerical instabilities in the exponential terms.

5 LOCAL STIFFNESS MATRIX

The stiffness matrix can be calculated as:

$K^{RBT}={\left\{\delta u\right\}}^T{\int }_0^{L}EA\left\{N_u\mathrm{\text{'}}\right\}{\left\{N_u\mathrm{\text{'}}\right\}}^{T}dx\left\{u\right\}+{\left\{\delta \theta \right\}}^T{\int }_0^{L}E{I}_z\left\{N_{\theta }\mathrm{\text{'}}\right\}{\left\{N_{\theta }\mathrm{\text{'}}\right\}}^{T}dx\left\{\theta \right\}$

$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}-\frac{1}{5}{\left\{\delta \theta \right\}}^T{\int }_0^{L}E{I}_z\left\{N_{\theta }\mathrm{\text{'}}\right\}{\left\{N_{\theta }\mathrm{\text{'}}\right\}}^{T}dx\left\{\theta \right\}+\frac{1}{5}{\left\{\delta \theta \right\}}^T{\int }_0^{L}E{I}_z\left\{N_{\theta }\mathrm{\text{'}}\right\}{\left\{N_{\theta }\mathrm{\text{'}\text{'}}\right\}}^{T}dx\left\{\theta \right\}$

$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}-\frac{1}{5}{\left\{\delta \theta \right\}}^T{\int }_0^{L}E{I}_z\left\{N_{\theta }\mathrm{\text{'}}\right\}{\left\{N_{\theta }\mathrm{\text{'}}\right\}}^{T}dx\left\{\theta \right\}+\frac{1}{21}{\left\{\delta \theta \right\}}^T{\int }_0^{L}E{I}_z\left\{N_{\theta }\mathrm{\text{'}}\right\}{\left\{N_{\theta }\mathrm{\text{'}}\right\}}^{T}dx\left\{\theta \right\}$ 32

$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}-\frac{1}{21}{\left\{\delta v\right\}}^T{\int }_0^{L}E{I}_z\left\{N_{\theta }\mathrm{\text{'}}\right\}{\left\{N_v\mathrm{\text{'}}\mathrm{\text{'}}\right\}}^{T}dx\left\{\theta \right\}+\frac{1}{5}{\left\{\delta v\right\}}^T{\int }_0^{L}E{I}_z\left\{N_v\mathrm{\text{'}}\mathrm{\text{'}}\right\}{\left\{N_0\mathrm{\text{'}}\right\}}^{T}dx\left\{v\right\}$

$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}-\frac{1}{21}{\left\{\delta v\right\}}^T{\int }_0^{L}E{I}_z\left\{N_v\mathrm{\text{'}}\mathrm{\text{'}}\right\}{\left\{N_{\theta }\mathrm{\text{'}}\right\}}^{T}dx\left\{\theta \right\}+\frac{1}{21}{\left\{\delta v\right\}}^T{\int }_0^{L}E{I}_z\left\{N_v\mathrm{\text{'}}\mathrm{\text{'}}\right\}{\left\{N_v\mathrm{\text{'}}\mathrm{\text{'}}\right\}}^{T}dx\left\{v\right\}$

$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}+\frac{8}{15}GA\left({\left\{\delta \theta \right\}}^T{\int }_0^L\left\{N_v\mathrm{\text{'}}\right\}{\left\{N_v\mathrm{\text{'}}\right\}}^{T}dx\left\{v\right\}+{\left\{\delta v\right\}}^T{\int }_0^L\left\{N_{\theta }\mathrm{\text{'}}\right\}{\left\{N_{\theta }\mathrm{\text{'}}\right\}}^{T}dx\left\{\theta \right\}\right)$

$\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}-\frac{8}{15}GA\left({\left\{\delta \theta \right\}}^T{\int }_0^L\left\{N_{\theta }\mathrm{\text{'}}\right\}{\left\{N_v\mathrm{\text{'}}\right\}}^{T}dx\left\{v\right\}+{\left\{\delta v\right\}}^T{\int }_0^L\left\{N_v\mathrm{\text{'}}\right\}{\left\{N_{\theta }\mathrm{\text{'}}\right\}}^{T}dx\left\{\theta \right\}\right)$

Using cubic shape functions, equation (32) returns$K^{RBT}$:

$\left[\begin{array}{cccccccc}\frac{EA}{L}& 0& 0& 0& -\frac{EA}{L}& 0& 0& 0\\ 0& \left(\frac{4EI}{7L^3}+\frac{16GA}{25L}\right)& \left(\frac{2EI}{7L^2}+\frac{4GA}{75}\right)& \frac{4GA}{15}& 0& \left(-\frac{4EI}{7L^3}-\frac{16GA}{25L}\right)& \left(\frac{2EI}{7L^2}+\frac{4GA}{75}\right)& \frac{4GA}{15}\\ 0& \left(\frac{2EI}{7L^2}+\frac{4GA}{75}\right)& \left(\frac{4EI}{21L}+\frac{16GAL}{225}\right)& \left(\frac{16EI}{105L}-\frac{2GAL}{45}\right)& 0& \left(-\frac{2EI}{7L^2}-\frac{4GA}{75}\right)& \left(\frac{2EI}{21L}-\frac{4GAL}{225}\right)& \left(-\frac{16EI}{105L}+\frac{2GAL}{45}\right)\\ 0& \frac{4GA}{15}& \left(\frac{16EI}{105L}-\frac{2GAL}{45}\right)& \left(\frac{68EI}{105L}+\frac{8GAL}{45}\right)& 0& -\frac{4GA}{15}& \left(-\frac{16EI}{105L}+\frac{2GAL}{45}\right)& \left(-\frac{68EI}{105L}+\frac{4GAL}{45}\right)\\ -\frac{EA}{L}& 0& 0& 0& \frac{EA}{L}& 0& 0& 0\\ 0& \left(-\frac{4EI}{7L^3}-\frac{16GA}{25L}\right)& \left(-\frac{2EI}{7L^2}-\frac{4GA}{75}\right)& -\frac{4GA}{15}& 0& \left(\frac{4EI}{7L^3}+\frac{16GA}{25L}\right)& \left(-\frac{2EI}{7L^2}-\frac{4GA}{75}\right)& -\frac{4GA}{15}\\ 0& \left(\frac{2EI}{7L^2}+\frac{4GA}{75}\right)& \left(\frac{2EI}{21L}-\frac{4GAL}{225}\right)& \left(-\frac{16EI}{105L}+\frac{2GAL}{45}\right)& 0& \left(-\frac{2EI}{7L^2}-\frac{4GA}{75}\right)& \left(\frac{4EI}{21L}+\frac{16GAL}{225}\right)& \left(\frac{16EI}{105L}-\frac{2GAL}{45}\right)\\ 0& \frac{4GA}{15}& \left(-\frac{16EI}{105L}+\frac{2GAL}{45}\right)& \left(-\frac{68EI}{105L}+\frac{4GAL}{45}\right)& 0& -\frac{4GA}{15}& \left(\frac{16EI}{105L}-\frac{2GAL}{45}\right)& \left(\frac{68EI}{105L}+\frac{8GAL}{45}\right)\end{array}\right]$ 33

However, employing the proposed approximation for displacements given by equations (19) to (21) and (28), and using interpolation functions shown in equations (27), the corresponding stiffness matrix is $K^{MRBT}$:

$\left[\begin{array}{cccccccc}\frac{EA}{L}& 0& 0& 0& -\frac{EA}{L}& 0& 0& 0\\ 0& \frac{12EI}{L^3}& \frac{6EI}{L^2}& 0& 0& -\frac{12EI}{L^3}& \frac{6EI}{L^2}& 0\\ 0& \frac{6EI}{L^2}& \left(\frac{64EI}{21L}+\frac{8GAL}{45}\right)& \left(\frac{16EI}{105L}-\frac{8GAL}{45}\right)& 0& -\frac{6EI}{L^2}& \left(\frac{62EI}{21L}+\frac{4GAL}{45}\right)& \left(-\frac{16EI}{105L}-\frac{4GAL}{45}\right)\\ 0& 0& \left(\frac{16EI}{105L}-\frac{8GAL}{45}\right)& \left(\frac{68EI}{105L}+\frac{8GAL}{45}\right)& 0& 0& \left(-\frac{16EI}{105L}-\frac{4GAL}{45}\right)& \left(-\frac{68EI}{105L}+\frac{4GAL}{45}\right)\\ -\frac{EA}{L}& 0& 0& 0& \frac{EA}{L}& 0& 0& 0\\ 0& -\frac{12EI}{L^3}& -\frac{6EI}{L^2}& 0& 0& \frac{12EI}{L^3}& -\frac{6EI}{L^2}& 0\\ 0& \frac{6EI}{L^2}& \left(\frac{62EI}{21L}+\frac{4GAL}{45}\right)& \left(-\frac{16EI}{105L}-\frac{4GAL}{45}\right)& 0& -\frac{6EI}{L^2}& \left(\frac{64EI}{21L}+\frac{8GAL}{45}\right)& \left(\frac{16EI}{105L}-\frac{8GAL}{45}\right)\\ 0& 0& \left(-\frac{16EI}{105L}-\frac{4GAL}{45}\right)& \left(-\frac{68EI}{105L}+\frac{4GAL}{45}\right)& 0& 0& \left(\frac{16EI}{105L}-\frac{8GAL}{45}\right)& \left(\frac{68EI}{105L}+\frac{8GAL}{45}\right)\end{array}\right]$ (34)

For instance, considering the simple clamped beam modeled with just one element (Figure 4), $\lambda =10$, numerically, the stiffness matrix for both beams theories are written as,

$K^{RBT}=\left[\begin{array}{cccccccc}5& 0& 0& 0& -5& 0& 0& 0\\ 0& 1.2332& 0.1038& 0.5128& 0& -1.2332& 0.1038& 0.5128\\ 0& 0.1038& 0.1375& -0.0848& 0& -0.1038& -0.0338& 0.0848\\ 0& 0.5128& -0.0848& 0.3446& 0& -0.5128& 0.0848& 0.1682\\ -5& 0& 0& 0& 5& 0& 0& 0\\ 0& -1.2332& -0.1038& -0.5128& 0& 1.2332& -0.1038& -0.5128\\ 0& 0.1038& -0.0338& 0.0848& 0& -0.1038& 0.1375& -0.0848\\ 0& 0.5128& 0.0848& 0.1682& 0& -0.5128& -0.0848& 0.3446\end{array}\right]10^5$ (35)

$K^{MRBT}=\left[\begin{array}{cccccccc}5& 0& 0& 0& -5& 0& 0& 0\\ 0& 0.0500& 0.0250& 0& 0& -0.0500& 0.0250& 0\\ 0& 0.0250& 0.3546& -0.3412& 0& -0.0250& 0.1832& -0.1716\\ 0& 0& -0.3412& 0.3446& 0& 0& -0.1716& 0.1682\\ -5& 0& 0& 0& 5& 0& 0& 0\\ 0& -0.0500& -0.0250& -0.5128& 0& 0.0500& -0.0250& 0\\ 0& 0.0250& 0.1832& -0.1716& 0& -0.0250& 0.3546& -0.3412\\ 0& 0& -0.1716& 0.1682& 0& 0& -0.3412& 0.3446\end{array}\right]10^5$ (36)

It can be observed that the classical Reddy stiffness matrix represents a major rigidity for the transversal displacement of the element provided by the terms and$16GA/25L$. However, the MRBT theory increases the rigidity related to the derivative of transversal displacement and their relation with the total cross-sectional rotation. These terms that are responsible to adjust the MRBT matrix to consider the shear deformation, otherwise, the displacement result would correspond to EBBT theory.

According with equations (29) to (31), the analytical solution for displacements at the free node of the clamped beam is obtained as:

$d^{ANALYTICAL}=\left\{\begin{array}{ccc}-0.8304& -1.2296& -1.2376\end{array}\right\}$ (37)

meanwhile, the Reddy beam models, using the stiffness matrix presented in equations (35) and (36), yields:

$d^{RBT}=\left\{\begin{array}{ccc}\text{-0.626373}& \text{-1.234984}& \text{-1.236254}\end{array}\right\}{d}^{MRBT}=\left\{\begin{array}{ccc}\text{-0.828820}& \text{-1.245640}& \text{-1.233590}\end{array}\right\}$ (38)

In the case of a slenderness ratio $\lambda$= 40, the analytical solution diverges encountering difficulties in computing the transversal displacement:

$d^{ANALYTICAL}=\left\{\begin{array}{ccc}~& -1.2330& -1.2335\end{array}\right\}$ (39)

whereas the numerical approximations return:

$d^{RBT}=\left\{\begin{array}{ccc}\text{-0.617229}& \text{-1.233345}& \text{-1.233428}\end{array}\right\}{d}^{MRBT}=\left\{\begin{array}{ccc}\text{-0.822575}& \text{-1.234012}& \text{-1.233261}\end{array}\right\}$ (40)

6 EQUIVALENT NODAL LOADS

In usual structures, the most general loading case considered to obtain nodal reactions is a linearly distributed force $q\left(x\right)=q_0+q_{1}x/L$, Figure 8.

Figure 8
Equivalent nodal loads.

This linear load can be written by composing a uniformly distributed load $q_0=q_i$ with a linear (triangular) load, $q_1=q_f-q_i$. The nodal equivalent loads can be calculated by using the interpolation functions related to the transversal displacement ($N_2^v,N_3^v,N_5^v,N_6^v$).

$f_i^{\mathrm{\text{'}}}=-{\int }_0^l{N}_i^v\left(x\right)\mathrm{}q\left(x\right)\mathrm{}dx=\mathrm{}{f}_{0i}+f_{xi}=-q_0{\int }_0^l{N}_i^v\left(x\right)dx-\frac{q_1}{L}{\int }_0^l{N}_i^v\left(x\right).xdx$ (41)

Therefore, the equivalent nodal loads considering the proposed Reddy beam formulation correspond to the Bernoulli-Euler beam theory. Since that the interpolation functions are the same (first line of Equation 27), leading to the usual equivalent nodal loads.

$f_{02}=\mathrm{}{f}_{05}=\mathrm{}\frac{-q_{0}L}{2};\mathrm{}{f}_{03}=\frac{-q_0{L}^2}{12};\mathrm{}{f}_{06}=\frac{q_0{L}^2}{12};\mathrm{}{f}_{x2}=\mathrm{}\mathrm{}\frac{-3q_{1}L}{20};\mathrm{}{f}_{x5}=\frac{-q_1{L}^2}{30};\mathrm{}{f}_{x3}=\frac{-7q_{1}L}{20};\mathrm{}{f}_{x6}=\frac{q_1{L}^2}{20}$ (42)

7 NUMERICAL APPLICATIONS

In this section, numerical examples were developed to verify the accuracy of the proposed finite element. Framed structures were modeled with both a high and a reduced level of discretization. The following descriptions identify elements in the examples:

• EBBT – Classic Bernoulli- Euler beam theory using Hermitian functions.

• TBT – Classic Timoshenko beam theory using cubic functions [Martha (2018)].

• RBT – Classic Reddy beam theory using Hermitian and linear interpolation functions for displacement and rotation, respectively, matrix (33).

• MRBT – Modified Reddy beam theory, proposed in this study. The stiffness matrix is given by the presented in matrix (34).

For consistency, the proposed formulation is compared with a converged result according to the usual Reddy model with a high level of discretization.

7.1 Clamped Beam

To validate the developed element (MRBT), the clamped beam depicted in Figure 9 has been analyzed. A reduced slenderness ratio $\lambda =L/h=10$ was employed to assess the element’s behavior. The beam have a length of L = 1 m, Young’s modulus of E = 10⁷ kN/m², section form factor of $\chi =5/6$ (in case of TBT element), and Poisson’s ratio $\nu =0.3$. The applied load was P = 1030 kN, and M = 1.03 kNm.

Figure 9
Clamped beam.

For EBBT and TBT the results are consistent, irrespective of structural discretization. On the contrary, the RBT model returns results depending on the adopted discretization, as shown in the convergence curve reported in Figure 10, where is compared the vertical displacement (v) at the free node for the RBT and the proposed model MRBT. It is evident that considering a single element the displacement resulting traditional RBT does not align with the correct result.

Figure 10
Displacement (v) convergence in a clamped beam (λ = 10) - RBT and MRBT model.

The modified model also yields results that vary with discretization. Nevertheless, the influence is significantly reduced, and even with just one element it is possible to obtain numerical results close to the analytical ones, as reported in Table 1.

Table 1
Transversal displacement (v) in a clamped beam (λ = 10) - RBT and MRBT model.

elements	dof	RBT		MRBT
		v	Diff	v	Diff
1	4	0.62761	24.53%	0.83005	0.18%
2	8	0.78207	5.96%	0.83074	0.10%
4	16	0.82048	1.34%	0.8312	0.05%
10	40	0.83051	0.13%	0.83142	0.02%
20	80	0.83146	0.02%	0.83156	0.00%
40	160	0.83159	0.00%	0.83159	0.00%

The displacements found using all beam theories (EBBT, TBT, RBT and MRBT) and a single element are summarized in Figure 11.

Figure 11
Displacement in a clamped beam (λ = 10) for distinct theories using one element.

Another crucial aspect to assess is related to shear stress. As is well-known, for the EBBT the shear stress distribution is undetermined, and the TBT leads to a constant shear stress distribution. The Reddy beam model improves upon these theories by providing a parabolic shear distribution. Thus, the MRBT needs to be evaluated in this context.

Figure 12 illustrates the convergence study for shear stress using both Reddy stiffness matrices, the classical RBT and the MRBT proposed in this work. The shear stress was assessed at the middle of the cross-section (maximum shear stress) at the free node of the clamped beam. The numerical values are presented in Table 2.

Figure 12
Shear stress (τ) convergence (clamped beam), λ = 10 - RBT and MRBT model.

Table 2
Shear stress (τ) in a clamped beam (λ = 10) - RBT and MRBT model.

elements	dof	RBT		MRBT
		τ	Diff	τ	Diff
1	4	4946	83.71%	-46320	52.57%
2	8	5821	80.83%	-26440	12.91%
4	16	1918	93.68%	-30480	0.40%
10	40	-12550	58.66%	-30640	0.92%
20	80	-23370	23.02%	-30400	0.13%
40	160	-27890	8.14%
80	320	-30330	0.10%
160	640	-30360	0.00%

From Figure 12, it is evident that achieving a convergence value with the standard Reddy beam stiffness matrix requires more than 300 degrees of freedom. As a consequence, the beam must be modelled with 80 or 160 elements, implying a high-level discretization to attain a correct evaluation of the shear stress. For example, a discretization with 40 RBT elements returns a difference from the converged result of approximately 8.14%, and it rises 93.68% considering 4 elements.

However, with the modified Reddy beam theory, better approximations are achieved even with reduced discretization. These considerably reduces the dimension of the system to be solved. Convergence is obtained with only 16 degrees of freedom, meaning that employing 4 elements is sufficient to reach the expected results. In such case, the difference between the shear stress obtained with 160 elements using the RBT element is only 0.40%.

Figure 13 illustrates the results from Table 2, showing the shear stress distribution in the cross-section of the free node, and how it changes depending on the discretization and the considered element.

Figure 13
Cross-section shear stress convergence (clamped beam), λ = 10 – (a) RBT and (b) MRBT.

The shear stress distribution along the beam is illustrated in Figure 14 (a) for the RBT element’s converged result. Figure 14 (b) demonstrated that the same result can be represented using the proposed MRBT element with a reduced number of elements.

Figure 14
Shear stress (τ) clamped beam (λ = 10) - (a) RBT with 160 elements and (b) MRBT model with 20 elements.

A similar study was conducted for a clamped beam considering a slenderness ratio , corresponding to h = 0.25 m. The applied load was P = 16050 kN. The convergence curve for the transversal displacement of the free node is presented in Figure 15, and numerical results are shown in Table 3.

Figure 15
Displacement (v) convergence in a clamped beam (λ = 4) - RBT and MRBT model.

Table 3
Transversal displacement (v) in a clamped beam (λ = 4) - RBT and MRBT model.

elements	dof	RBT		MRBT
		v	Diff	v	Diff
1	4	0.66920	22.41%	0.853058	1.09%
2	8	0.821443	4.76%	0.857346	0.59%
4	16	0.855768	0.78%	0.860152	0.27%
10	40	0.86183	0.07%	0.861812	0.08%
20	80	0.862317	0.02%	0.862293	0.02%
40	160	0.862473	0.00%	0.862468	0.00%

Another noteworthy point is that, despite the similarity of the MRBT stiffness matrix with the EBBT, the displacement is more aligned with the results obtained from the Timoshenko beam model when a small slenderness is considered, and this alignment is achieved without the need of the section form factor. This observation is illustrated in Figure 16 with a single element.

Figure 16
Displacement in a clamped beam (λ = 4) for distinct theories using one element.

The modified model also yields results that vary with discretization. Nevertheless, the influence is significantly reduced, and even with just one element it is possible to obtain numerical results close to the analytical ones, as reported in Table 3.

The convergence of shear stress is shown in Figure 17, and numerical results are presented in Table 4. Additionally, Figure18 illustrates the parabolic shear stress distribution along the cross-section for various theories and the adopted discretization.

Figure 17
Shear stress (τ) convergence in clamped beam (λ = 4) - RBT and MRBT model.

Table 4
Shear stress (τ) in a clamped beam (λ = 4) - RBT and MRBT model.

elements	dof	RBT		MRBT
		τ	Diff	τ	Diff
1	4	22420	88.44%	-288600	48.84%
2	8	-5081	97.38%	-165200	14.80%
4	16	-79980	58.75%	-190200	1.91%
10	40	-161500	16.71%	-191100	1.44%
20	80	-182000	6.14%	-192900	0.52%
40	160	-187400	3.35%
80	320	-193400	0.26%
160	640	-193900	0.00%

Figure 18
Cross-section shear stress convergence (clamped beam), λ = 4 – (a) RBT and (b) MRBT.

The results indicate that employing the MRBT leads to shear stress with accurate approximation using 4 elements, with a difference of 1.91% from the converged solution (160 RBT elements). In contrast, using the RBT formulation, it is necessary to use 40 or more elements to achieve a similar difference. The shear stress distribution for the RBT along the clamped beam is reported in Figure 19 (a), and in Figure 19 (b) for the MRBT element.

Figure 19
Shear stress (τ) clamped beam (λ = 4) - (a) RBT and (b) MRBT with 160 elements.

Finally, the last example evaluated for the clamped beam considers a slenderness ratio of $\lambda =100$, thus the height of the cross section corresponds to h = 0.01 m. The convergence curve for the transversal displacement of the free node is presented in Figure 20 and in Table 5.

Figure 20
Displacement (v) convergence in a clamped beam (λ = 100) - RBT and MRBT model.

Table 5
Transversal displacement (v) in a clamped beam (λ = 100) - RBT and MRBT model.

elements	dof	RBT		MRBT
		v	Diff	v	Diff
1	4	0.619319	24.95%	0.825284	0.01%
2	8	0.773819	6.23%	0.825290	0.01%
4	16	0.812445	1.54%	0.825295	0.01%
10	40	0.823259	0.23%	0.825298	0.01%
20	80	0.824804	0.05%	0.825299	0.01%
40	160	0.825188	0.00%	0.825299	0.01%

Also, the displacement for all beam theories (EBBT, TBT, RBT, MRBT) cited in this work is shown in Figure 21.

Figure 21
Displacement in a clamped beam (λ = 100) for distinct theories using one element.

Despite of this similar behavior observed for the displacements in the clamped beam examples previously described, in this specific case achieving convergence for shear stress requires 640-1280 elements (2560-5120 degrees of freedom) when using the usual RBT model. However, when the proposed MRBT element is employed, accurate results are obtained with just 4 elements, mirroring the outcomes seen in the other examples, as illustrated in Figure 22 and summarized in Table 6.

Figure 22
Shear stress (τ) convergence in a clamped beam (λ = 100) - RBT and MRBT model.

Table 6
Shear stress (τ) in a clamped beam (λ = 100) - RBT and MRBT model.

elements	dof	RBT		MRBT
		τ	Diff	τ	Diff
1	4	50.86	83.43%	-463.3	50.91%
2	8	71.23	76.80%	-264.3	13.91%
4	16	73.78	75.97%	-304.7	0.75%
10	40	66.3	78.40%	-306.3	0.23%
20	80	47.26	84.61%
40	160	16.7	94.56%
80	320	-87.58	71.47%
160	640	-208.8	31.99%
320	1280	-279	9.12%
640	2560	-301	1.95%
1280	5120	-307	0.00%

The shear stress distribution along the cross section presented in Figure 23 also evidences the better convergence of the proposed element.

Figure 23
Cross-section shear stress convergence (clamped beam), λ = 100 – (a) RBT and (b) MRBT.

The shear distribution along the clamped beam is shown in Figure 24. However, in this case, using just 160 elements, it is not possible for the classical RBT theory to represent the converged result.

Figure 24
Shear stress (τ) clamped beam (λ = 100) - (a) RBT and (b) MRBT with 160 elements.

7.2 Portal Frame

This example assesses the proposed formulation in a more complex system, the portal frame depicted in Figure 25. The frame was modeled with a length of L = 1 m, and slenderness ratio of $\lambda =10$. To model the frame using the Timoshenko beam theory, a section form factor of $\chi =5/6$ is considered. The material parameters remain consistent with those used in the clamped beam example. The applied load P is 1000 kN.

Figure 25
Portal frame.

Following the same methodology as in the initial example, a convergence curve was elaborated for the displacement (u) at the top of the left column, considering both the RBT and the proposed MRBT. The results of the convergence study are presented in Figure 26. Additionally, the numerical results obtained for both theories are summarized in Table 7.

Figure 26
Displacement (u) convergence in a portal frame (λ = 10) - RBT and MRBT model.

Table 7
Transversal displacement (u) in a portal frame (λ = 10) - RBT and MRBT model.

elements	dof	RBT		MRBT
		v	Diff	v	Diff
1	4	0.007549	94.91%	0.144108	2.75%
2	8	0.115235	22.23%	0.147283	0.61%
4	16	0.140934	4.89%	0.147774	0.27%
10	40	0.147571	0.41%	0.148094	0.06%
20	80	0.148180	0.00%	0.148212	0.02%

From this analysis, it is evident that the RBT converges to a value of 0.14818 with 20 elements (80 degrees of freedom). However, when employing just one element, this theory yields to a result with a 94.91% error concerning the converged displacement. This discrepancy diminishes with an increased level of discretization, reaching acceptable results with 4-10 elements, with differences of 4.89% and 0.41%, respectively.

Using the proposed model in this work, it is possible to model the structure with just one element for both columns and beams, akin to the behavior of the EBBT model. In this scenario, the difference from the converged displacement is approximately 2.75%. When utilizing 2 elements, this difference further decreases to less than 1%. Figure 27 illustrates the displacement for the distinct beam theories obtained with just one element in each member.

Figure 27
Displacement in a portal frame (λ = 10) for distinct theories using one element.

To evaluate the shear stress, the middle of the beam is analyzed. The convergence study and numerical result for the maximum shear stress (which occurs in cross section centroidal axis) are shown in Figure 28 and Table 8, respectively.

Figure 28
Shear stress (τ) convergence in a portal frame (λ = 10) - RBT and MRBT model.

Table 8
Shear stress (τ) in a portal frame (λ = 10) - RBT and MRBT model.

elements	dof	RBT		MRBT
		τ	Diff	τ	Diff
2	20	-1702	113.39%	19670	54.76%
4	44	-668.2	105.26%	10810	14.95%
10	116	5384	57.64%	12800	0.71%
20	236	9938	21.81%	12770	0.47%
40	476	11950	5.98%
80	956	12550	1.26%
160	1916	12710	0.00%

This example shows that using 10 MRBT elements leads to a difference of the order of 0.71% related to the result adopted. While for the RBT element, if just 10 elements are used, this difference is around of 57.64%. When using 4 elements, the difference for MRBT is of 14.95%, while for the RBT this difference is 105.26%. This result of the discretization influence on the cross-section shear distribution is illustrated in Figure 29.

Figure 29
Cross-section shear stress convergence (beam of portal frame) – (a) RBT and (b) MRBT.

The shear stress distribution along the beam of the portal frame is shown in Figure 30. Using a high level of discretization (160 elements), both theories are capable of simulating the shear stress distribution, giving the same results.

Figure 30
Shear stress (τ) in the beam of the frame - (a) RBT and (b) MRBT with 160 elements.

7.3 Frame with Distributed Loads

The main objective of this example is to evaluate the proposed model when is acting distributed loads and with simply supported conditions. It is proposed an frame as shown in Figure 31. As in the other examples, the frame was modeled with a length of L = 1 m, slenderness ratio of $\lambda =10$, and a section form factor of $\chi =5/6$. The material parameters remain consistent with previous examples. The distributed load q was adopted as 10000 kN/m.

Figure 31
Frame distributed loads.

A convergence study presented in Figure 32 was performed to verify the structural response with the discretization variation. For this example, the rotation of the node 3 (support of the beam) was verified. The numerical results for both Reddy beams theories (RBT and MRBT) are given in Table 9.

Figure 32
Rotation ($\theta$) convergence in a frame with distributed load ($\lambda =10$) - RBT and MRBT model.

Table 9
Rotation ($\theta$) in a frame with distributed load ($\lambda =10$) - RBT and MRBT model.

elements	dof	RBT		MRBT
		θ	Diff	θ	Diff
1	8	0.29398	54.35%	0.653515	1.48%
2	16	0.58176	9.66%	0.640769	0.50%
4	32	0.635934	1.25%	0.642391	0.25%
10	80	0.646368	0.37%	0.643111	0.14%
20	160	0.644985	0.16%	0.643383	0.09%
40	320	0.643981	0.00%

The analysis shows that to reach a converged result for the cross-section rotation, the RBT element requires a discretization of 20-40 elements (160-320 dof). Using just one RBT element, the difference from the converged rotation is of 54.35%. In fact, is necessary at least 4 elements to reach acceptable results.

However, it can be verified that even considering distributed loads and simply supported conditions, the MRBT presents very reduced discretization influence. If the proposed element (MRBT) is considered with just 1 element, the difference from the converged result is about 1.48%, reducing to 0.50% with 2 elements.

Figure 33 shows the result for the rotation considering the distinct beams theories using one element per member. It is evident that employ the MRBT element permit a structural analysis for displacements similar to theories of EBBT and TBT, i.e, with a minimum number of elements required to describe the geometry.

Figure 33
Rotation in a frame with distributed load (λ = 10) for distinct theories using one element.

A convergence study for the shear stress was also performed (Figure 34). The shear stress considered for this study occurs at the node 3 (support of the beam), at the cross section centroidal axis (representing the maximum shear stress in this point). The obtained numerical values are presented in Table 10.

Figure 34
Shear stress (τ) convergence in a frame with distributed loads (λ = 10) - RBT and MRBT model.

Table 10
Shear stress (τ) in a portal frame with distributed loads (λ = 10) - RBT and MRBT model.

elements	dof	RBT		MRBT
		τ	Diff	τ	Diff
2	16	-30600	125.48%	90370	24.75%
4	32	-16260	113.54%	96750	19.44%
10	80	40950	65.90%	111500	7.16%
20	160	88960	25.93%	116800	2.75%
40	320	111300	7.33%	119300	0.67%
80	640	118200	1.58%
160	1280	120100	0.00%

The results evidence that to reach the shear stress using the usual RBT model is necessary a high level of discretization, around of 80 to 160 elements (640-1280 dof). If a reduced discretization is considered the difference from the response is around 125.48% (2 elements), or 113.54% (4 elements), and also the sign is inverted. Acceptable results are presented with a discretization of 40 elements.

However, modeling the structure with the proposed element (MRBT), this reduced discretization of 40 elements leads to a difference only of 0.67% from the response. Use 2 or four elements will lead to relevant differences, but this difference is only achieved with more than 20 elements in the classical RBT model. Acceptable results are reached with a discretization of 10 elements.

The influence of the discretization on the shear stress distribution is shown in Figure 35. The results demonstrate that the convergence is achieved more easily when the proposed element is applied. On the other hand, the usual RBT element demands a high level of discretization to reach reasonable values.

Figure 35
Cross-section shear stress convergence (beam of frame with distributed load) – (a) RBT and (b) MRBT.

Finally, the shear stress distribution is predicted with similar behavior if a hard number of elements are used. Figure 36 present this distribution considering 160 elements for both theories (RBT and MRBT).

Figure 36
Shear stress (τ) in the beam of the frame with distributed load - (a) RBT and (b) MRBT with 160 elements.

8 CONCLUSION

In this study, an enhanced shape function for the Reddy beam model is presented. It allows framed structural analysis using just one element in each member, a practice commonly employed with both the Bernoulli- Euler and Timoshenko beam models. This feature facilitates easy integration into any FEM code. The proposed displacement interpolation also proves to be effective in predicting structural behavior when analytical solution fall short.

The formulation underwent comprehensive evaluation through various numerical examples. Initially, we compared the proposed interpolation functions with the analytical solution for the Reddy beam model, demonstrating the element capability to reproduce this solution without the stability problems deriving by the exponential terms contained in the closed-form solution. Subsequently, a detailed analysis of a clamped beam was conducted, examining displacements and shear stress distribution for different structural slenderness ratios. Finally, was performed a study on a framed structure, evaluating displacements and shear stress distribution along the entire beam and cross-section.

The results demonstrate the efficiency of the proposed element, able to predict accurately both the displacements and the stress field with just one element for each member, a capability lacking in conventional theories (Bernoulli- Euler or Timoshenko) implemented for structural analysis of frames. Notably, the shape functions usually adopted for the Reddy beam model (Reddy 2022) necessitates higher-order discretization for correct displacement and shear stress determination, a challenge effectively mitigated with the proposed approach.

Acknowledgement

This work has been supported by University of Campania “Luigi Vanvitelli” [D.R. 869/2022], FAPES [Grant 380/2023 P 2023-6KPV7], Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) [Finance Code 001], the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) [Grant 308884/2021-3], and FAPERJ [process number E-26/201.224/2022].

References

• AkhavanAlavi, S., Mohammadimehr, M. and Edjtahed, S. (2019). Active control of micro Reddy beam integrated with functionally graded nanocomposite sensor and actuator based on linear quadratic regulator method. Eur. J. Mech. A Solids, 74:449–461.
• Almeida, C. A., et al. (2011). Geometric nonlinear analyses of functionally graded beams using a tailored Lagrangian formulation. Mech. Res. Commun., 38 8:553–559.
• Augello, R., Carrera, E. and Pagani A. (2019). Unified theory of structures based on micropolar elasticity. Meccanica, 54:1785–1800.
• Bathe, K. J. (1996). Finite Element Procedures. Prentice-Hall, Englewoods Cliffs, New Jersey.
• Bickford, W. B. (1982). A consistent higher order beam theory. Develop. Theor. and Appl. Mech., 11:137-150.
• Chakrabarti, A., et al. (2012). Analysis of composite beams with longitudinal and transverse partial interactions using higher order beam theory. Int. J. Mech. Sci., 59:115–125.
• Darijani, H. and Mohammadabadi, H. (2014). A new deformation beam theory for static and dynamic analysis of microbeams. Int. J. Mech. Sci., 89:31–39.
• Dong, S. B., Alpdogan, C., and Taciroglu, E. (2010). Much ado about shear correction factors in Timoshenko beam theory. Int. J. Solids Struct., 47:1651-1665.
• El Meiche, N., et al. (2011). A new hyperbolic shear deformation theory for buckling and vibration of functionally graded sandwich plate. Int. J. Mech. Sci., 53 4:237–247.
• Faroughi, S., Rahmani, A. and Friswell, M. (2020). On wave propagation in two-dimensional functionally graded porous rotating nano-beams using a general nonlocal higherorder beam model. Appl. Math. Model, 80:169–190.
• Golbakhshi, H., Dashtbayazi, M. R. and Saidi, A. R. (2022). A New Theoretical Framework for Couple Stress Analysis of Reddy-Levinson Micro-Beams. Int. J. Appl. Mech., 14 7:1-24.
• He, G. and Yang, X. (2014). Finite element analysis for buckling of two-layer composite beams using Reddy’s higher order beam theory. Finite Elem. Anal. Des., 83:49-57.
• Heyliger, H. and Reddy, J. N. (1988). A higher order beam finite element for bending and vibration problems. J. Sound Vib., 126 2:309–326.
• Iurlaro, L., Gherlone, M. and Di Sciuva, M. (2015). The (3,2)-Mixed Refined Zigzag Theory for generally laminated beams: Theoretical development and C0 finite element formulation. Int. J. Solids Struct., 73-74:1-19.
• Karama, M., Afaq, K. S. and Mistou S. (2003). Mechanical behaviour of laminated composite beam by the new multi-layered laminated composite structures model with transverse shear stress continuity. Int. J. Solids Struct., 40:1525–1546.
• Karamanli, A. and Vo, T. P. (2020). Size-dependent behaviour of functionally graded sandwich microbeams based on the modified strain gradient theory. Compos. Struct., 246:112401.
• Karamanli, A., Vo, T. P. and Civalek, O. (2023). Finite element formulation of metal foam microbeams via modified strain gradient theory. Eng. Comput., 39:751–772.
• Karttunen, A. T., et al. (2017). Bridging plate theories and elasticity solutions. Int. J. Solids Struct., 106-107: 251-263.
• Keleshteri, M. and Jelovica, J. (2022). Analytical assessment of nonlinear forced vibration of functionally graded porous higher order hinged beams. Compos. Struct., 298:115994.
• Khorshidi, M. A., Shariati, M. and Emam, S. A. (2016). Postbuckling of functionally graded nanobeams based on modified couple stress theory under general beam theory. Int. J. Mech. Sci., 110:160–169.
• Leite, L. F. P. and Rocha, F. C. (2023). A Novel Higher-Order Zigzag Function Applied to Refined Unified Beam Theory for the Analysis of Composite Laminated Materials. Periodica Polytechnica Civil Engineering, 67 3:867–874.
• Levinson, M. A. (1981). A new rectangular beam theory. J. Sound. Vib., 74 1: 81-87.
• Lin, F., et al. (2020). Assessment of first and third order shear deformation beam theories for the buckling and vibration analysis of nanobeams incorporating surface stress effects. Int. J. Mech. Sci., 186:105873.
• Mahi, A., Adda Bedia, E. A. and Tounsi, A. (2015). A new hyperbolic shear deformation theory for bending and free vibration analysis of isotropic, functionally graded, sandwich and laminated composite plates. Appl. Math. Model., 39 9:2489–2508.
• Martha, L. F. (2018). Matrix structural analysis with object oriented programming (in Portuguese). Elsevier, Rio de Janeiro.
• Martha, L.F. and Rangel, R.L. (2022). "FTOOL: Three decades of success as an educational program for structural analysis”. Proceedings of the XLIII Ibero-Latin-American Congress on Computational Methods in Engineering, ABMEC/UNILA, Foz do Iguaçú, PR.
• McGuire, W., Gallagher, R. H. and Ziemian, R. D. (2000). Matrix Structural Analysis. John Wiley & Sons Inc, New York.
• Meghare, T. K. and Jadhao, P. D. (2015). A simple higher order theory for bending analysis of steel beams. Int. J. Struct. Civ. Eng., 2: 31-38.
• Mittelstedt, S. and Mittelstedt, C. (2020). Mixed-mode buckling of shear-deformable composite laminated I-beams. Int. J. Mech. Sci., 169:105332.
• Mohammad-Abadi, M. and Daneshmehr, A. R. (2014). Size dependent buckling analysis of microbeams based on modified couple stress theory with high order theories and general boundary conditions. Int. J. Eng. Sci., 74:1–14.
• Molina-Villegas, J. C. and Ortega, J. E. B. (2023). Closed-form solution of Timoshenko frames using the Green’s Function Stiffness Method. Int. J. Solids Struct., 269:112180.
• Muc, A. (2020). Non-local approach to free vibrations and buckling problems for cylindrical nano-structures. Compos. Struct., 250:112541.
• Nazemnezhad, R. and Zare, M. (2016). Nonlocal Reddy beam model for free vibration analysis of multilayer nanoribbons incorporating interlayer shear effect. Eur. J. Mech. A Solids., 55:234–242.
• Neves, A., et al. (2011). Bending of FGM plates by a sinusoidal plate formulation and collocation with radial basis functions. Mech. Res. Commun., 38 5:368–371.
• Nguyen, N. D., et al. (2022). A new two-variable shear deformation theory for bending, free vibration and buckling analysis of functionally graded porous beams. Compos. Struct., 282: 115095.
• Petrolito, J. (1995). Stifness analysis of beams using a higher-order theory. Comput. Struct., 55:33-39.
• Polizzotto, C. (2015). From the Euler–Bernoulli beam to the Timoshenko one through a sequence of Reddy-type shear deformable beam models of increasing order. Eur. J. Mech. A Solids., 53: 62–74.
• Rangel, R. L. and Martha, L. F. (2019). LESM—An object-oriented MATLAB program for structural analysis of linear element models. Comput. Appl. Eng. Educ., 27 3:553-571
• Reddy, J. N. (1984a). A refined nonlinear theory of plates with transverse shear deformation. Int. J. Solids Struct., 20:881-896.
• Reddy, J. N. (1984b). A simple higher order theory for laminated composite plates. J. Appl. Mech., 51 4:745–752.
• Reddy, J.N. (1997). On locking-free shear deformable beam finite elements. Comput. Methods Appl. Mech. Eng., 149:113-132.
• Reddy, J. N. (2007). Nonlocal theories for bending, buckling and vibration of beams. Int. J. Eng. Sci., 45:288–307.
• Reddy, J. N. (2019). An Introduction to the Finite Element Method, 4th. McGraw-Hill, New York.
• Reddy, J.N. (2022). Theories and Analyses of Beams and Axisymmetric Circular Plates. CRC Press, Boca Raton, FL.
• Reddy, J. N., Wang, C. M. and Lee, K. H. (1997). Relationships between bending solutions of classical and shear deformation beam theories. Int. J. Solids Struct., 34-26:3373-3384.
• Rodrigues, M. A. C., Burgos, R. B. and Martha, L. F. (2019). A unified approach to the Timoshenko geometric stiffness matrix considering higher-order terms in the strain tensor. Lat. Am. J. Solids Struct., 16:1-22.
• Rodrigues, M. A. C., Burgos, R. B. and Martha, L. F. (2021a). A unified approach to the Timoshenko 3D beam-column element tangent stiffness matrix considering higher-order terms in the strain tensor and large rotations. Int. J. Solids Struct., 222–223 111003:1-22.
• Rodrigues, M. A. C., Burgos, R. B. and Martha, L. F. (2021b). Complete tangent stiffness matrix considering higher-order terms in the strain tensor and large rotations for a Euler Bernoulli – Timoshenko space beam-column element. MethodsX, 8 101591:1-37.
• Ruocco, E. and Reddy, J. N. (2021). A discrete differential geometry-based approach to buckling and vibration analyses of inhomogeneous Reddy plates. Appl. Math. Model., 100:342–364.
• Ruocco, E. and Reddy, J. N. (2023). Analytical solutions of Reddy, Timoshenko and Bernoulli beam models: A comparative analysis. Eur. J. Mech. A Solids, 99:104953.
• Ruocco, E., Reddy, J. N. and Wang C. (2020). An enhanced hencky bar-chain model for bending, buckling and vibration analyses of Reddy beams. Eng. Struct., 221:111056
• Sahmani, S., Bahrami, M. and Ansari, R. (2014). Nonlinear free vibration analysis of functionally graded third-order shear deformable microbeams based on the modified strain gradient elasticity theory. Compos. Struct., 110:219–230.
• Sayyad, A. S. and Ghugal, Y. M. (2017) Bending, buckling and free vibration of laminated composite and sandwich beams: A critical review of literature. Compos. Struct., 171:486–504.
• Shi, H. and Voyiadjis, G. (2011). A Sixth-order theory of shear deformable beams with variational consistent boundary conditions. J. Appl. Mech., 78 2:1–11.
• Simsek, M. and Reddy, J. N. (2013). A unified higher order beam theory for buckling of a functionally graded microbeam embedded in elastic medium using modified couple stress theory. Compos. Struct., 101:47–58.
• Stojanović, V. and Petković, M. D. (2016). Nonlinear dynamic analysis of damaged Reddy-Bickford beams supported on an elastic Pasternak foundation. J. Sound Vib., 385:239-266.
• Szenkrényes, A. (2014). Bending solution of third-order orthotropic Reddy plates with asymmetric interfacial crack. Int. J. Solids Struct., 51:2598-2619.
• Timoshenko, S. P. (1921). On the correction factor for shear of the differential equation for transverse vibrations of bars of uniform cross-section. Philos. Mag., 744.
• Uddin, Md. A., et al. (2017a). A higher order model for inelastic response of composite beams with interfacial slip using a dissipation based arc-length method. Eng. Struct., 139:120–134.
• Uddin, Md. A., et al. (2017b). Large deformation analysis of two layered composite beams with partial shear interaction using a higher order beam theory. Int. J. Mech. Sci., 122:331–340.
• Uddin, Md. A., et al. (2018).Geometrically nonlinear inelastic analysis of steel–concrete composite beams with partial interaction using a higher-order beam theory. Int. J. Non. Linear Mech., 100:34–47.
• Uddin, Md. A., et al. (2020). Convergence studies of finite element model for analysis of steel-concrete composite beam using a higher-order beam theory. Struct., 27:2025–2033.
• Vidal, P. and Polit, M. (2008). A family of sinus finite elements for the analysis of rectangular laminated beams. Compos. Struct., 84:56–72.
• Vidal, P. and Polit, O. (2011). A sine finite element using a zig-zag function for the analysis of laminated composite beams. Composites: Part B., 42:1671–1682.
• Vinh, P. V. and Son, L. T. (2022). A new first-order mixed beam element for static bending analysis of functionally graded graphene oxide powder-reinforced composite beams. Struct., 36:463–472.
• Vo, T. P. and Thai, H.-T. (2012). Static behavior of composite beams using various refined shear deformation theories. Compos. Struct., 94:2513–2522.
• Wanji, C., Chen, W. and Sze, K. (2012). A model of composite laminated Reddy beam based on a modified couple-stress theory. Compos. Struct., 94 8:2599–2609.
• Yang, Y. B. and Leu, L. J. (1994). Non-linear stiffnesses in analysis of planar frames. Comput. Methods Appl. Mech. Eng., 117;233-247.
• Yang, Y. B. and Kuo, S. R. (1994). Theory & Analysis of Nonlinear Framed Structures. Prentice Hall, Simon & Schuster (Asia) Pte ltd, Singapura.
• Zhen, W., et al. (2019). Stability of laminated composite and sandwich beams by a Reddy-type higher-order zig-zag theory. Mech. Adv. Mat. Struct., 26 19:1622-1635.
• Zheng, Y.-F., et al. (2023). Nonlinear postbuckling analysis of magneto-electro-thermo-elastic laminated microbeams based on modified couple stress theory. Appl. Math. Model., 118:89–106.

Edited by

Editor: Marcílio Alves

Publication Dates

• Publication in this collection
  05 July 2024
• Date of issue
  2024

History

• Received
  14 Mar 2024
• Reviewed
  27 Apr 2024
• Accepted
  02 May 2024
• Published
  07 May 2024