An 8-Nodes 3D Hexahedral Finite Element and an 1D 2-Nodes Structural Element for Timoshenko Beams, Both Based on Hermitian Intepolation, in Linear Range

The following article presents the elaboration and results obtained from a 3D finite element, of the 8-node hexahedron type with 6 degrees of freedom (DOF) per node (48 DOF per element) based on third degree Hermitian polynomials, and of a 2-node structural element, with 6 DOF per node (12 DOF per element), based on third degree Hermitian polynomials and the theory of Timoshenko for beams. This article has two purposes; the first one is the formulation of a finite element capable of capturing bending effects, and the second one is to verify whether it is possible to obtain the deformation of the beam’s cross section of a structural member of the beam type, based on the deformations of its axis. The results obtained showed that the 8-node hexahedron FE was able to reproduce satisfactory results by simulating some cases of beams with different contour and load conditions, obtaining errors between 1% and 4% compared to the ANSYS software, educational version. Regarding the structural element of the beam type, it reproduced results that were not as precise as the FE Hexa 8, presenting errors of between 6% and 7% with regard to the axis but with error rounding between 10% and 20%.


Introduction
Currently there are a large number of analytical and finite element models [1][2][3][4][5][6][7] capable of reproducing the behavior of beams with accuracy and precision. Generally, finite element models demand a considerable computational cost depending on the analysis that is proposed and the meshing established by the user.
Regarding recent developments employing Timoshenko's beam theory, ref. [8] derived a non-local formulation using Timoshenko's Theory and non-local elasticity, formulating the elasticity matrix using two parameters that combined the assumptions of Winkler and Pasternak. For their part, ref. [9] used the Theory of Functional Connections (TFC) to solve numerical boundary value problems. The authors solved equations that govern the response of beams using the Timoshenko-Ehrenfest Theory, which incorporates transverse shear strains.
To solve Timoshenko beams with nonlinear behavior at both geometric and material levels, the authors of [10] developed a formulation that directly applied three-dimensional nonlinear constitutive equations without zero-stress conditions. The researchers compared the results with those obtained by applying brick-type elements, managing to verify the precision and efficacy of the formulation. In another recent work, the authors of [11] incorporated interpolation functions in the Timoshenko 3D beam solution, which were obtained directly from the solution of the equilibrium differential equation of an infinitesimal deformed element, thus reducing the discretization of elements used. The interpolation functions were used in the definition of the tangential stiffness matrix in an updated Lagrangian formulation. The formulation was used to solve various three-dimensional examples with elements using different bending theories, including Timoshenko's, achieving accurate prediction of results with minimal discretization.
To simulate the rotation phenomenon of a Timoshenko beam with different slenderness ratios, the authors of [12] implemented a 3D finite element formulation. The novelty introduced in the finite element consisted of correctly considering the moment of inertia under small rotational deformations around the equilibrium position. The authors carried out various numerical simulations that allowed them to compare the effectiveness of the formulation against commercial finite element programs. The dynamic response of Timoshenko beams has also been studied by the authors of [13], who studied perforated beams under accelerated loads in a variable thermal environment. The equations of the system have been solved numerically by means of implicit integration in time. The results obtained show correspondence with respect to the results obtained by applying conventional methods.
Another important phenomenon is the one involving lateral-torsional buckling in steel beams subjected to high temperatures. To study it, the authors of [14] have used a formulation based on Timoshenko beam-type finite elements, thus avoiding using shelltype elements that transform the problem into a much more complex and computationally expensive one. They implemented two strategies: the first one corresponded to a fiberbased beam model (FBM) and the second one corresponded to a cruciform frame model (CFM). The results obtained from the numerical simulations show an adequate reproduction of the local buckling using a reduced computational time, which represents a promising projection towards the performance-based design of this type of beam when subjected to the action of fire.
In another recent study, this time focused on the dynamic response under moving loads of fiber-reinforced Timoshenko composite beams, the authors of [15] used a Lagrangian procedure to determine the constitutive equations of motion accompanied by the use of the Ritz Method. The results obtained with the numerical simulations demonstrated the effectiveness of the formation, which, however, depends on the orientation of the reinforcing fibers and the volume fraction on the dynamic response of the models. Continuing with the Timoshenko beams, the authors of [16] studied the effect of complex loads applied at multiple nodes. The authors used the inverse finite element method (iFEM) as a basic conceptual framework and Isogeometric Analysis (IGA) as a complementary strategy to update the shape of beams in two dimensions. Numerical simulations were compared against experimental results, indicating acceptable errors for deformations.
As for finite elements applied to structural members such as beams, there are no limits to modeling, including deformations in the cross section and in the longitudinal axis, while with analytical models, only beams that have specific formulations can be modeled at a low computational cost, limiting the types of analysis of non-deformable beams that can be carried out on their cross section, and their longitudinal axis can also be modeled at a low computational cost . Commonly, for buildings design, the cross section deformation is neglected, but when the analysis requires cross section deformation, or the steel reinforcement contribution for concrete beams, it is mandatory to evaluate it by using a finite element (FE) modeling or a fiber modeling.
Most of 3D FE are based on tethraedrons or hexaedrons, generally with 3 DOF each node, whose shape functions are generally linear (2 nodes per edge) or quadratic (3 nodes per edge); this means, for each polinomyal degree, an extra node is needed; with this formulation, the beam's cross section and steel reinforcement can be modeled [17][18][19][20]. For 1D FE, generally the cross section is considered through the geometrical properties within the formulation, and 6 DOF per node for full axial, shear, and bending effects. When the formulation is treated as 1D, the effects of the reinforcement in concrete beams can be represented by transforming the steel reinforcement into concrete, due to the limitations of the formulation.
When beams are being formulated, with a 3D FE or 1D FE, it is necessary a formulation capable to capture the bending and shear effects; for a 3D FE, the bending and shear effects are obtained by the three main displacements of each node, with a quadratic shape function (minimum), due to the complexity of the formulation; for a 1D FE, the formulation is quite easy, with a full 6 DOF each node (axial forces, shear forces, and bending moments). There are two main theories for beam formulation: Bernoulli's formulation, which neglects the shear deformation (valid for small beams in heigth), and Timoshenko's Theory, which considers the shear deformation (valid for all types of beams) [21][22][23][24][25][26].
In this article, an Hexa 3D 8-node (48 DOF per element) and a 1D 2-node (12 DOF per element) will be developed, both based on their weak formulation and the Hermitian interpolation. This formulation has the advantage of being able to capture bending effects directly with only 8 nodes in a regular hexahedron, in the isoparametric space [7]. The hexaedron element and its DOF, in the isporametric space (which domain is Ω = (−1, 1) for each edge), number of nodes and local axis, is shown in Figure 1. The hermitian polynomial for the interpolation was selected as third degree polynomial, as shown in Equation (1), which will be used to generate the shape functions for the three local axis (ξ, η, ζ).

Materials and Methods
Rotation at any point of P(ξ) will be ∂ ∂ξ P(ξ), and, if the boundary conditions are taken in 2D for each edge as θw 1 , θw 2 , w 1 , w 2 (Figure 2), the constants a, b, c, d can be calculated as shown in Equation (2).
w 1 and w 2 represent the displacement perpendicular to the edge axis for each node, while θw 1 and θw 2 represent the rotation of each node. Solving Equation (2) for constants a, b, c, d, with each one as a function of θw 1 , θw 2 , w 1 , w 2 , then grouping for θw 1 , θw 2 , w 1 , w 2 and expanding the solution for each local axe, the shape functions are obtained as shown in Equation (3). It has to be noticed that the same procedure can be applied for axis η and ζ, obtaining the same equations.
In order to construct the shape functions for each node, in the isoparametric space, a superposition of all them are proposed, as shown in Figure 3, where H 1 , H 3 represent displacements, and H 2 , H 4 represent rotations. By making the superposition of the corresponding shape functions for each node, the equations for displacements and rotations are shown from Equations (4)- (11). Notice that this element is 'complete', it means, if any equation for displacement is evaluated in any point with its coordinate, the result will be one (1), and for rotations, the result will be zero (0).
Once obtained, the shape functions in the isoparametric space, a transformation to the Cartesian space is required, through the Jacobian matrix. The shape functions can be expressed in Cartesian coordinates according to Equations (12)-(14) In matrix form, the Jacobian matrix can be expressed as shown in Equation (15), and the transformation to the Cartesian coordinates is shown in Equation (16).

Elasticity Parameters
The deformation tensor, based on the six DOF per node, is shown in Equation (17) (matrix B), and the displacements matrix is shown in Equation (18). Notice that Equation (18) is showing 6 rows and 6 columns, but it is only for 1 node; hence, for the 8 nodes, the dimension of the matrix will be by 6 × 48.
The properties of the materials are shown in the Equation (19).

Numerical Integration
The numerical integration used for calculations was based on Gaussian method, with 3 integration points for edge. The three main points used for the calculations were 3 , 0, 5 3 ], and w = 5 9 , 8 9 , 5 9 . The 27 integration points were obtained by the permutation of x and w.

Elemental Stifness and Forces Matrices
The strong formulation for the element stiffnes matrix is shown in Equation (20), while the strong formulation of the forces matrix is shown in Equation (21). The e subfix denotes the elemental matrix.
The strong formulation is not viable to evaluate, due to the complexity of the coordinate transformation between the isoparametric space and the Cartesian space; hence, a finite elements formulation based on the weak formulation, which implies numerical integration to the Equations (20) and (21), which become Equations (22) and (23) where: B is the deformation tensor (Equation (17)). Notice that this tensor depends directly of the shape functions (η, ξ, ζ), which asume the values of the Gaussian points. N is the displacements matrix (Equation (18)) D is the matrix which contains the materials properties (Equation (19)) q(x) is the Neumann condition or boundary conditions (external forces) n is the number of integration points The global matrices K and F are obtained from the assembly of K e and F e .

Formulation of the 1D 2-Nodes Structural Element Based on Timoshenko's Beams Theory
In this section, a formulation of a 1D 2-nodes structural element with 6 DOF per node is based on Timoshenko's beams theory starting from a variational formulation (this type of formulation is based on virtual work commonly), and also based on the main differential equation of displacement. Once this FE is formulated, the cross section will be deformated using geometric relationships in order to compare if, for simple beams, a deformed shape is based on the axis displacement.

2D Variational Formulation (Weak Formulation)
Similar to the deduction of the beam under Bernoulli's theory, the resolution of the Equation (24) leads to a varational formulation.
where I is the second moment of area of the cross section, E is the Young's module, v(x) is a displacement function, and q(x) is the applied external force. Following the general procedure, a solution is proposed by means of a test function, as shown in the Equation (25), which is integrated to reduce the order, as shown in Equation (26); this is w(x) C 4 (Ω, ), where Ω = (0, L) is the domain of the beam; v : Ω → is the displacement at any point of the beam; and q : Ω → is the Neumann condition or boundary conditions (external forces).
where: E, is the Young's module. Iz, is the second moment of areas around the axis submited to bending moments.
When solving the integral by parts, and knowing that the test function w(x) must have the same boundary condition (Dirichlet) as the function v(x), Equations (27) and (28) are obtained.
Equation (28) is solved again through integrating by parts, leading to Equation (29), where the values w(L), w(0) and their derivatives evaluated at 0 and L are inputs; so, the final Equation to solve is expressed as Equation (30). M(0) is the bending moment at the initial node, and M(L) is the bending moment at the final node.

Finite Elements Solution
A Hermitian polynomial is proposed for the solution of Equation (3), as shown in Equation (31). v Based on Timoshenko's beam theory, the additional shear deformation can be appended to the polynomial of Equation (31), posing a balance of shear forces, as shown in Equation (32), where the term Ks GA γ represents the shear force. Timoshenko's beam theory states that, unlike the Euler-Bernoulli's beam theory, where the cross section does not rotate after loads are applied to it, the shear distortion v is considered as an additional rotation of the cross section. Therefore, transverse rotation and transverse displacement are considered as independent variables [22].
where Iz represents the second moment of area around the bending axis, Ks is a shear correction factor depending of the cross section, G is the shear module, and A is the shear area.
In order to include this shear deformation, the value of γ is determined from Equation (32), obtaining Equation (33), and this value is substituted in the rotation expression, as shown in Equation (34).
To solve Equation (34), the boundary conditions shown in Figure 2 are set, obtaining the expressions of the integration constants, as shown in Equation (35).
Returning to Equation (30), the terms of stiffness and force are obtained, which are shown in Equations (40)-(42) Using Equation (41) for the indices i = 1, 2, 3, 4 and j = 1, 2, 3, 4, the stiffness matrix for bending and shear is obtained, shown in Equation (43), while the volume forces term is used Equation (42); it is shown in Equation (44). The ful 3D stiffness matrix is shown in Figure A2, and the schema for the ful 3D DOF is shown in Figure A3.
It can be seen that if shear effects are neglected due to the effects explained in Timoshenko's beam theory, that is, α = 0, the matrix K e is equal to the one obtained by applying the Euler-Bernoulli Theory. Since what is sought is to obtain a matrix of 6 degrees of freedom per node (12 × 12 matrix), the same process can be carried out by taking into account the rotation around the z axis, resulting in a 4 × 4 K e matrix, including Iy instead of Iz, and including axial and torsional effects.

Validation of the Hexahedral FE with 8 Nodes
Firstly, the finite element Hexa 8 was validated for various types of beams subjected to different types of load, comparing the results with the ANSYS ® software. The cross section is 20 cm in base and 40 cm in height (Concrete, fc = 24.52 MPa, E = 23,414 MPa), to which the minimum amount of steel reinforcement has been introduced and transformed into rectangular bars. The minimum reinforcing steel yielded 6 bars of 14 mm diameter, which yielded rectangular bars 12.5 mm on a side (fy = 411.9 MPa, E = 205,940 MPa). In Figures 4-7, the graphic at the right represents the modeling in ANSYS ® software [27], while the graphic at the left represents the output of the finite element Hexa 8, displayed in Paraview ® software [28] and calculated using Python ® software [29]. The mesh was constructed with a Python script, and there were three different types of mesh for the cross section ( Figure A1 and Table A1), using the medium mesh for the validation.  The error in each simulation is shown in Table 1, where the maximun error is about 3.21%, which means that the results for Hexa 8 FE are reliable.

Validation of the 1D 2-Nodes Structural Element for Timoshenko Beams
To determine if the deformed beam with the hypothesis is accurate in comparisson with the Hexahedral FE with 8 nodes, the evaluation of three different beams with different Dirichlet and Neuman conditions was proposed, as shown in the following paragraphs.

One Side Fixed Beam
This beam has the section already defined and has a length of 1 m. The results for this beam are shown in Figures 8-10. The left Figure represents the EF simulation, while the right represents the Timoshenko simulation. Notice that in these results, the error has been calculated for the Timoshenko beam formulation in comparison with the Hexahedral FE with 8 nodes. The error for each simulation, for the Timoshenko beam, is shown in Table 2. Notice that errors for Timoshenko beam formulation are not negible for the torsional case (case 3), and there are several reasons for this; one reason is that for Timoshenko beam formulation, the supports were added to a single point, allowing displacements for the nodes above and below the support nodes, while for the Hexahedral FE with 8 nodes, supports were added to a surface, restricting displacements for all the support nodes. These errors, from the point of view of structural analysis, should be not allowed, due to the fact that for structural analysis the main goal is to reproduce the behavior of a real structure as well as posible; on the other hand, for structural design, even when these error values are high, for the estimation of structural cross sections, steel reinforcement, or steel profiles, these errors probably do not affect the design. For instance, a common structural parameter for acceptation or rejection of beams is the deflection or maximun vertical displacement, which is expressed in terms of the total length of the beam (Length/480, Length/360, among others), and if the displacements are in the order of 10 −3 cm, with these errors, the beam still can be accepted.
Another reason could be that the deformed shape for Timoshenko beam formulation is based on linear geometric relationships, which cannot be accurate in order to estimate the deformed cross section, while with the Hexahedral FE with 8 nodes, the displacements of each node are calculated based on the FE formulation explained above. In addition, these errors were calculated as an average of the error of each point, which can be unfavorable if one point presents a huge error.

Double Pinned Beam
The results are presented for a beam with a length of 2 m, with pinned supports at its ends (Figures 11-13). The left Figure   The error for each simulation, for the Timoshenko beam, is shown in Table 3.  The error for each simulation, for the Timoshenko beam, is shown in Table 4.

Discussion
According to the results obtained, it is possible to affirm that although the results appear to be satisfactory for the structural element under Timoshenko's theory, the error calculated in relative terms is high in some cases. In this paper, two different FE were created and evaluated, with expected results: the Hexahedral FE with 8 nodes performed better that the 1D FE Timoshenko beam formulation. In comparison, the axis displacements were relatively close, but the deformed cross section were not, essencialy with torsional effects.
In the case of the beam embedded at one end, it is observed that in the first case (Table A2), the error is 24.34% on average, while for case 2, the error remains quite low, at 3.56%, and for the torsional case, the values again increase to 17.00% on average, which can mean important changes in the design of structures susceptible to displacement, by underestimating the maximum deflections; however, the behavior of the cross section in all cases is similar, and in addition to this, when observing the results just at the maximum length (L = 1 m), where the maximum deflection occurs, the error is low , being the lowest of the entire table in cases 1 to 3; therefore, from the point of view of maximum deflections, the error is low.
In the case of the beam articulated at both ends, it is observed that according to the Table A3, the error for case 1 is low, at 5.14%, while the error for case 2 increases to 9.07%, and in the case of torsion, it continues to increase to approximately 17%, again, this can cause important changes in structures susceptible to displacement. Moreover, the same behavior is presented, so that the error values at the maximum deflection point (L = 1 m) are the lowest; so, an estimate of deflections is quite close.
Something similar happens with the cantilever beam (Table A4), observing high average percentages in cases 1, 4, 5, and 6, but not in cases 2 and 3; in case 1, it is observed that in the section between L = 0 m and L = 2 m, the lowest error occurs in the zone of maximum deflection, and in the same way in L = 3 m, but the error is high in absolute terms. For cases 2 and 3, the behavior is the same, with a lower error value, so it can be considered acceptable. With respect to the other cases, the error has been calculated without taking into account the values of 100 because this comes from the way in which the degrees of freedom of both models are restricted, showing that for the FE model, there are some displacements that are translated between sections, while in the Timoshenko model, the displacements are restricted on the axis and no deformation is transmitted to the cross section.

Conclusions
As a conclusion, it can be established that the finite element 3D Hexa 8 created with Hermite polynomials accurately represents the results of the simulations when compared with a commercial software, with the disadvantage that when implemented in Python-Cython, the computation times are quite high. This finite element 3D Hexa 8 has the advantage to represent the rotations of the nodes with the translation and rotation component, which is the correct way to represent the real behavior of an element submitted to axial, bending and shear forces, in comparisson with some finite elements which the rotations are based only on the translations of each node.
The FE 3D formulation for the hexahedron was easy to formulate with third degree Hermitian polynomials due to the number of equations nedded for representing bending effects and adding the corresponding DOF for rotation, in comparison with common formulations. The proposed FE is based on hermitian polynomial interpolation, which considerably reduces the quantity of nodes in the hexahedron element and catches the bending effect more accurately by adding three additional DOF of rotation for each node. Rotations in common FE are caught by adding more nodes in the hexahedron and can include three or more nodes for each Edge, with 3 DOF for each one but only in translation, and rotations are calculated based on the translations.
It is also well known that the addition of rotations of DOF using hermitian polynomial interpolation is used for beam formulation (similar to the Timoshenko beam formulation presented in this paper) and for plates, which work with the axis of the structural element, but in this paper, the using hermitian polynomial interpolation is used directly in the FE formulation, in order to reduce the quantity of nodes and add three more DOF for rotations.
The advantage of this FE could lie in the processing time, but it depends because for a common FE formulation, there has to be a minumum of 16 nodes in order to correctly catch the bending effect (which is not very precise, because the polynomical order is barely 2), generating 48 equations per element, so the number of equations increases, but, using the Hexahedral FE with eight nodes, the number of equations is the same (48) per element but the order of the interpolation is 3 (used in this paper), not 2, so, with a common FE formulation (3 DOF per node) there will need 24 nodes, generating 72 equations per element, and it will always catch the rotations based on the translations. Based on the paragraph above, the advantages of using an FE lie in accuracy obtained by adding three additional DOF for rotations.
With respect to the Timoshenko simulations, it is observed that the results oscillate between high and low for the different beams in the different cases, so it is recommended to continue with the development until a better approximation is achieved in terms of deflections. One of the possible causes of these differences is the way in which the restrictions are carried out since by EF, these are carried out on a surface or on a line, while by the Timoshenko simulation, they are carried out at a point on the shaft, and this might make a difference.
In some cases, the Timoshenko simulations and the finite element 3D Hexa 8 matched almost perfectly, which implies that for academic purposes, it is possible to use the total section deformed shape based on the axis deformations, and by applying some trigonometric relations.
The accuracy of the finite element model must obviously be contrasted with the experimental results carried out on beams with geometry and mechanical characteristics similar to those used in this work.   Note: the matrix is symmetric, and the shear factor depends on the direction of the bending, that is, as a function of Iz e Iy.