A C1 Beam Element Based on Overhauser Interpolation

A new C1 element is proposed to model Euler-Bernoulli beams in one and two-dimensional problems. The proposed formulation assures C1 continuity requirement without the use of rotational degrees of freedom, used in traditional elements, through the use of an Overhauser interpolation scheme for bending displacements. The principle of virtual displacements is used to determine the equilibrium equations and boundary conditions for one and twodimensional Euler-Bernoulli beams. The Overhauser interpolation is introduced and the new bending interpolation functions are defined. Finally, beam and frame problems are solved with the new formulation and the results are compared to the traditional Euler-Bernoulli element and exact solutions.


INTRODUCTION
In the classical Euler-Bernoulli beam theory the main assumption is that the cross sections remain plane and orthogonal to the beam axis after deformation. This makes the beam shear undeformable and the shear stresses remain indeterminate. This theory presents good results for isotropic homogeneous thin beams.
In order to account for shear, the assumption of orthogonality to the beam axis was removed by the Timoshenko beam theory, but the cross sections remain plane. Then, the beam is made shear deformable, but the predicted shear strain and stress proved to be constant through the beam thickness. Which is in contradiction with the principle that shear stress is zero at the beam surface. Therefore, a shear correction factor was introduced to account the effect of this deficiency of the Timoshenko beam theory.

Latin American Journal of Solids and Structures 14 (2017) 92-112
The applications of those two theories in Finite Element Method are already standardized, as in Reddy (1993).
Other beam theories were formulated to include shear effects. Among those theories is the Reddy third-order theory where the application for beams was developed as a specific case of plates. This theory was used by Heyliger and Reddy (1988) to develop a beam finite element and study bending and vibrations of isotropic beams. Reddy (1997) later defined locking-free elements based on the simplified third-order theory, having the Euler-Bernoulli and Timoshenko solution as special cases. Such elements can be used to obtain exact nodal displacements and forces (from the equilibrium equations of the element) in a frame structural analysis using one finite element per structural member. Which is not achieved with the traditional equal interpolation, reduced integration element because the elements exhibit shear locking for thick members unless two or more elements per structural member are used.
Polizotto (2015) presented shear deformable beams of increasing order where a Reddy beam without shear coincides with the Euler-Bernoulli beam, and a Reddy beam shear deformable of infinite order can be identified with the Timoshenko beam.
Nevertheless, in all presented theories and their derived elements, rotational degrees of freedom are present in order to guarantee C 1 continuity at the expense of increasing the number of degrees of freedom (DOF) necessary to model beams and frames. The objective of this work is to eliminate the rotational degrees of freedom related to bending and yet still satisfy the C 1 continuity requirement necessary for the convergence of the finite element method applied to beams, reducing the computational burden involved in the solution. In order to achieve this goal an Overhauser interpolation scheme will be used to generate a cubic interpolation functions.
The new formulation proposed in the present work is developed based on the Euler-Bernoulli assumptions, where the rotational degrees of freedom are related to the transverse displacements derivatives only. Overhauser (1968) proposed an interpolation scheme for curves and surfaces where, for curves, a linear blending of implicit parabolas defines a cubic interpolation between each pair of adjacent points, as an exclusive function of the point's positions. Applying this concept to finite elements method, the beam and frame elements formulated in the present work possess cubic interpolation functions to represent bending based only in nodal displacements. But as it will be explained further, the application of the Overhauser interpolation is not straight forward and some solutions were proposed in order to incorporate such interpolation in finite elements formulation.
The Overhauser interpolation has already been explored to solve numerical methods problems to guarantee C 1 continuity. In Boundary Elements Method (BEM), Overhauser quadrilateral elements were first applied by Hibbs (1987, 1988). Later, Walters and Gipson (1994) demonstrated that Overhauser elements give more accurate results then the standard quadratic and cubic elements based on Lagrange interpolation in BEM and Hadavinia et al. (2000) developed generalized parabolic blending C 1 elements using Overhauser concept enhancing the accuracy of the solution over non-uniform meshing in BEM. Archer (2006) proposed continuous solutions from the Green element method using Overhauser elements, as an effective alternative to the Hermitian approach presented by Taigbenu (1998). And in Finite Element method Ulaga et al. (1999) applied the Overhauser interpolation in contact problems of gears, in this case, the contact area was ap-proximated with Overhauser spline interpolation functions over boundary nodes but the geometry was completely decoupled from the finite element polynomial shape functions.
In the sequence of this work, the finite element formulation based on the classical beam theory for the beams and frames will be constructed using the principle of virtual displacements (PVD). The Euler-Bernoulli beam and frame interpolation functions will be briefly introduced. Then, Overhauser interpolation will be presented and the new bending interpolation functions using Overhauser scheme will be developed and applied to one and two-dimensional formulations. Typical beam and frame problems will be solved using the new elements and results will be compared to the analytical solution, for the one dimensional problems, and to the classical Euler-Bernoulli beam solution obtained using an academic license of the commercial Finite Element Analysis (FEA) software ABAQUS from Dassault Systèmes.

THE EULER-BERNOULLI BEAM FORMULATION
The governing equations of a three-dimensional beam can be obtained applying the Principle of Virtual Displacements (PVD). Then, one-dimensional and two-dimensional equations can be generated with the correct simplifications.
The displacement field for a three-dimensional beam is represented by: The strain-displacement relations are: In the Euler-Bernoulli (classical) beam theory, it is assumed that the plane cross sections perpendicular to the axis of the beam remain plane and perpendicular to the beam axis after deformation. Therefore θz = v,x, θy = -w,x and as a result the terms of the virtual work related to shear vanish. Then, the strain-displacement relations become: Another assumption of this theory is that no normal stresses are present in the y and z directions. Along with the non-deformation of the cross-section, we can verify that σyy = σzz = τyz = 0.
The principle of virtual displacement states: Where δWi is the virtual work of the internal strains and δWe is the external efforts virtual work. The first term is defined as: Where {σ} is a column matrix in the Voigt notation of the Cauchy tensor components and {δe} is a column matrix in the Voigt notation of the tensor [δe], this tensor is a quadratic function of the displacement gradients in Lagrangian coordinates.
Assuming the displacement gradients are small in modulus, compared to the unity, it is valid to approximate the displacement gradient with the components of the Green Strain Tensor. Then, the virtual work of the internal stresses can be written as Applying the assumptions of the classical beam theory for an isotropic material δWi reduces to: The virtual work of external loads for a three dimensional case considering distributed loads, concentrated forces and moments is written as: As the main objective is to propose a formulation without rotational degrees of freedom for bending we can conclude that this formulation does not support, at first, prescribed moments as boundary conditions. But later, we will see that this is not entirely true. Also, three-dimensional applications, where torsion is involved, are impractical. Because in three-dimensional problems bending couples with torsion at the junctions of structural members with different orientations. Therefore, we cannot eliminate the degrees of freedom only related to bending. Reducing the applicability of a three-dimensional element to structures composed by only one member. Then, no advantage arises from the use of the proposed element compared to the traditional Euler-Bernoulli element in 3D problems.
However, in one and two-dimensional applications, this limitation does not exist. Because only one rotation degree of freedom is present and it is exclusively related to bending.

THE OVERHAUSER INTERPOLATION
Overhauser (1968) proposed that for a set of points, implicit parabolas are defined using three consecutive points in a way that each adjacent pair of points belongs to two different parabolas. Then a linear blend of the two parabolas results in a C 1 continuous cubic function.
In order to apply this idea in finite element formulation, the points coordinates are replaced by the nodal transverse displacements and the implicit parabolas are converted into parabolic Lagrange interpolations of those displacements, as illustrated in Figure 1. The si independent variable originally related to the arc-length is naturally linked to the length of the elements. And the two parabolas vA and vB are defined as: Where coefficients a0, a1, a2, b0, b1, b2 are functions of the nodal displacements vi-1, vi, vi+1 and vi+2.
The linear blending functions A and B are defined by: Finally, the Overhauser curve, defined in the same interval is: Evaluating the properties of such curve we notice that, as A and B are linear functions of si and vA and vB are quadratic functions of si thus, the interpolation v(si) in Eq. 11 is a cubic function of si. And if we evaluate the values of v at point i, the independent variable assumes zero value and vA(si = 0) = vB(si = 0) = v(si = 0) = vi. Also, at point i+1, si = li and vA(si Latin American Journal of Solids and Structures 14 (2017) 92-112 The first derivative of v(si) is given by.
Then the first derivatives of v with respect to si at nodes i and i+1 are Therefore, it is seen that dv/dsi varies continuously in between each pair of adjacent points. So, the curve v(si), can adequately represent an inflection that should occur within an interval. It should be apparent that the cubic function defined in the interval between nodes i and i+1 does not pass through points i1 and i+2, this is an important difference in comparison to simple fit by cubic polynomials. The manner of construction (as a blend of two parabolas) guarantees that spurious wiggles will not be introduced, as frequently happens when simple cubics are forced to pass through four points of a curve.
However for the extremities, Overhauser (1968) uses the parabolic interpolation and the cubic interpolation derived applies only in the elements that share two parabolas, what we will call central elements (CE).
For those elements transverse displacement can be written as: Where ϕ i-1 , ϕ i , ϕ i+1 , and ϕ i+2 are the interpolation functions for the central element (Eq. 28 in Appendix A). Those functions are defined in the interval 0 ≤ si ≤ li. Assuming all elements possess the same length, li-1 = li = li+1 = l, we can plot the interpolation functions for a better understanding of its behavior. In Figure 2, we see clearly a local support character in the interpolation functions, and that the influence of the displacements of the nodes outside the element (i-1 and i+2) is much smaller than the ones of the element nodes (i and i+1). Now, in order to construct cubic interpolation functions for the elements on the extremities, some considerations must be made. First, their interpolations must satisfy C1 continuity with the Overhauser interpolation already formulated for the central elements. Second, as it was said before, no rotational degrees of freedom are present, as a consequence, no essential nor natural boundary conditions of this character can be applied. Nevertheless, for many beam and frame problems this information is primordial.
Considering those points, a solution proposed is to build an element with a rotational degree of freedom only at the node not connected to central elements (first or last node of the structural member; i.e. node A in Figure 3). And then eliminate this degree of freedom via static condensation already considering the essential boundary condition. This strategy allows applying the different boundary conditions and accounting the effect of those boundary conditions in the element stiffness matrix. For the left corner element, the new curve vC is a cubic defined in the interval li-1 ≤ si ≤ 0 by the following equation.
where coefficients c0, c1, c2 and c3 are obtained imposing the conditions, with a1 as the same coefficient from the parabola vA used in the blending of the central element to obtain a cubic for the interval from i to i+1, and θi-1 is the rotation in the point i-1.
Then the curve vC can be written as a function of the nodal displacements and the rotation θi-1 as: Where ηθ, ηi-1, ηi, and ηi+1 are the interpolation functions for the left corner element (Eq. 29 in Appendix A). Again assuming all elements possess the same length, li-1 = li = li+1 = l, we can plot the interpolation functions to understand its behavior. Again a local support character is observed in Figure 4, and we see that rotation influence is small, as also the displacement of the node outside the element (node C in Figure 3). But, as stated earlier, the main propose of this work is to eliminate the dependency of the rotational degrees of freedom. This can be accomplished by performing a static condensation of the degree of freedom θi-1. The bending problem for the element can be written as:  Considering the essential boundary condition of the node that possesses θi-1, there are two possibilities that imply two different stiffness matrices and load vectors. If the node is clamped θi-1 = 0, then we just need to eliminate the first row and the first column of the matrix But if this node is simply supported or free θi-1 is unknown, then we have to isolate θi-1 in the first equation of the system presented in Eq. 18.
Therefore, it is also possible to apply a concentrated moment at the degree of freedom θi-1 when we have this type of essential boundary condition. Then, the stiffness and load components related to θi-1 are redistributed in an equivalent problem. The elimination procedure in element level is only possible because the rotational degree of freedom is not shared with other elements. And as a result, the essential boundary condition is already considered in the element level, being simply assembled in the global problem. For the right corner element an analogous procedure is adopted to develop cubic curve vD and the interpolation functions are defined in the interval 0 ≤ x ≤ li+1 as presented in Figure 5, resulting: The interpolation functions ψθ, ψi, ψi+1, and ψi+2 presented in Eq. 30 of Appendix A and plotted in Figure 6 (considering equal length for all elements) also guarantee C 1 continuity with the central elements and originate a stiffness matrix and load vector that consider the essential boundary condition at the unshared node.

PLANE FRAME FORMULATION
For one-dimensional problems, the developed functions can be directly applied. But for twodimensional analyses, some considerations must be made. First, the axial displacements interpolation functions are the same linear interpolation functions used in traditional Euler-Bernoulli elements.
Consequently, as in traditional elements, only the degrees of freedom that belong to the element are taken in consideration. Contrary to the new bending formulation that uses information of nodes outside the element.
Second, the different orientations of the elements must be accounted and coherent transformation matrices must be elaborated in order to correctly transform the stiffness and load components to the global coordinate system. As they were formulated the Overhauser interpolation functions for bending assume the same orientation for all nodes present in the cubic formulation.
Therefore, the transformation matrices must act also on the components related to the nodes external to the elements for the bending effects.
In order to illustrate this concept we can observe Figure 7 and see how the orientation of the nodes changes from one element to another. If we imagine that elements e 2 and e 3 are defined as central elements, we note that for element e 2 (Figure 7.a) the nodes 1, 2 and 3 have a different orientation compared to the same nodes when they associate with the element e 3 (Figure 7 Then, in element e 2 (Figure 7.a), for node 3 a stiffness related to the transverse displacement y on the element (local) coordinate system would be transformed into a axial stiffness in the global coordinate system (GCS). This is not physically coherent nor numerically, inducing bad results. Because as stated before, axial behavior is not affected by nodes outside the element.
In the light of the above we propose that, when we have elements with different orientations sharing a node (i.e.: elements e 2 an e 3 sharing node 2 in Figure 7) we will treat each element that shares this node as left corner elements (LCE). With all matrices being generated including the rotational degree of freedom. Then, they will be transformed to the GCS and assembled in a superelement. Finally, the rotational degree of freedom can be eliminated by a static condensation procedure in element level analogous to the one presented before.
This solution numerically assures C 1 continuity, avoids incoherent coupling of stiffness and allows representing structures such as the one presented in Figure 8, where elements e 1 , e 3 and e 5 share node 0. And e 3 and e 5 are rotated of α 2 and α 3 respectively, in relation to the GCS (element e 1 LCS coincides with the GCS). Once the stiffness matrices and load vectors of these three elements have been transformed to the GCS they can be assembled as presented in Figure 9. Then, the elimination of the rotational degree of freedom θ 0 can be performed with a static condensation, resulting the final stiffness matrices and load vector of the superelement, which can finally be assembled in the global stiffness matrix and load vector of the problem.
It is important to remark that the transformations to the global coordinates system should have been applied before the elimination and the assembly, if not, the stiffness related to the θ0 will be redistributed in a wrong way.

NUMERICAL RESULTS
Here we study the behavior of the Overhauser elements discussed on this paper. First, five onedimensional analyses of different load conditions and essential boundary conditions are tested and compared to the exact solution according Euler-Bernoulli beam theory. Then, two frame problems are solved. For the convergence analysis, we adopted the criterion of 1% error in the displacements of specific points (not the whole structure) compared to the exact solutions for the beam problems, and compared to the Euler-Bernoulli solutions using ABAQUS for the frame problems.
The beam and frame configurations tested are presented in Figure 10.

Beam Problems
For the beam problems, the exact solutions using the classical beam theory according to Hibbler (2011) are: Where v is the transverse displacement, L is the beam length and x is the independent variable related to the length. And P represents the concentrated force, where in all cases a force P = 1000 N was applied, except for the distributed load case where a value of w = 1600 N/m was chosen in order to result the same maximum displacement of the force P. The maximum displacements were analyzed for the convergence tests, except for problem 5 where the displacement at x = 3L/4 was monitored.  We see in all problems that the Overhauser solution is more rigid, as expected of a numerical approximation. And despite being cubic interpolation functions they do not represent as well the solution as the Hermite cubic functions. As a consequence, more elements are required and the computational burden is bigger (more degrees of freedom) than using Euler-Bernoulli traditional elements.
It is also noted that the essential boundary conditions influence the convergence rate. Comparing problems 2 and 3, where the natural boundary conditions are the same, but in 2 the extremities are simply supported and in 3 they are clamped, we see that the convergence rate reduces by half for problem 3, even with both solutions being cubic. This behavior is not observed with traditional Euler-Bernoulli elements for this situation. But if we compare problems 2, 4 and 5 with the same essential boundary conditions we see that the convergence rates are similar.
For problem 4, where not only the stiffness is integrated using the interpolation functions, but also the load, we see that results are coherent, indicating good behavior of the interpolation functions in this aspect.
Finally, post-processing was performed in order to extract information about nodal rotation v'(x) and maximum bending stress σmax(x) (at the nodes and element centroid) for all onedimensional problems.
C 1 continuity was attested, confirming the mathematical formulation and also good behavior was attested concerning the critical points that were correctly captured. But as the displacement presented errors to the exact solution, those errors were passed to v'(x) and σmax(x) as consequence.
The errors to the exact solution as function of the beam length for v(x), v'(x) and σmax(x) are presented in Figure 11 to Figure 13 respectively, on the domains where the exact solutions are defined (Eqs. 23 to 27).  In Figure 12 and Figure 13 the arrows indicate that the errors extrapolate. And for this reason, we preferred to indicate its tendency.
It is interesting to note that the essential boundary conditions interfere not only in the convergence rate but at the errors magnitude. As it was briefly explained before, the clamped and simplysupported/free elements have different stiffness matrices. Also, we see that for clamped situations, the errors are bigger at the corner elements decreasing in the other elements. While for the cases where simply-supported condition exists, the errors are bigger next to the points were v' = 0 (critical points). Where, for all problems in which a critical point occurs in central elements, the error in displacement tends to 1%, as presented in Figure 11.
The same tendency is observed in the nodal rotation errors presented in Figure 12. And in problem 5 the error in v'(x) next to the critical point rises until 23% closer to where the first derivative changes sign. But as it was said before, the critical points are correctly captured.
Finally, concerning stresses, the corner elements are again the ones with the bigger errors, being the simply-supported / free corner elements the ones that exhibit the biggest errors, in the order of 73% for the problems 1, 2 and 5 where we have a concentrated force and 62% for problem 4 where we have a distributed load.
Therefore, we can conclude that both central and corner elements present errors, but it is clear that the proposed solution for the corner elements is more critical, especially when it comes to the second derivative of the interpolation functions.

Plane Frame Problems
For the frame problems the global coordinate system XY was defined in Figure 10 and the structural members with length L 1 are vertical (parallel to Y) and the ones with length L 2 are horizontal (parallel to X).
The L-frame of problem 6 and the U-frame of problem number 7 both have clamped bases and elements with equal length were used.
Both structures were first submitted to concentrated forces at the junctions of the structural members that produced simple traction to the vertical segments in order to validate the superelement solution according to the principles of the classical beam theory. According to this theory, for such loading, only axial deformation of the vertical segments and simple translation of the horizontal segments should occur. What indeed happened, validating not only the solution, but the implementation of the axial behavior as well, obtaining the same results of Euler-Bernoulli traditional elements within ABAQUS.
The same does not happen when we use central elements to model the edge. And a transverse deflection appears in both segments.
As an example, the displacements of L-frame structure with the force F 1 = 10 4 N at the junction of the structural members (see Figure 14) are presented in Table 2. The structure was modeled with six elements (three per structural member) using the superelement (SE) solution and also using central elements (CE) at the junction. As stated before, modeling the junctions using central elements induces coupling between axial and transverse behavior that is not in accordance with the classical beam theory. This can be seen in the different nodal displacements on Y direction of the nodes of the horizontal member. But the superelement solution proved to be consistent with the problem physics.

Node
With SE With CE X (10 -6 m) Y (10 -6 m) X (10 -6 m) Y (10 -6 m) 0  For the L-frame other two loading situations were tested, as shown in Figure 15. The convergence criterion was 1% error in both X and Y displacements of the free end. First a force F 2 = 1000 N was applied. This problem resembles the cantilever beam problem analyzed before. Convergence was indeed achieved using the same nine elements in the vertical segment (18 elements for the whole structure, see Figure 17) as in the cantilever beam and no axial displacement was exhibited by the nodes in the vertical member, as the classical theory dictates. Also, the displacements in the nodes of the horizontal member (nodes 9 to 18 for the converged solution) in the X direction increase linearly, as shown in Figure 16. Also in accordance with the classical beam theory.  Second, the horizontal force F 2 was replaced by a vertical force F 3 = 100 N at the same point.
Inducing bending in both members and axial deformation in the vertical member. The mesh convergence test in Figure 18 shows that the solution achieved the convergence criterion with already eight elements (four per structural member). Also it reveals that the X-displacement of the endnode is the same that the one obtained using Euler-Bernoulli traditional elements even with six elements. Figure 18: Mesh convergence test for L-frame with load F3 at the free end.
Comparing the two load cases and their respective convergence analyses ( Figure 17 and Figure  18), we see that, curiously, the displacements in the direction of the concentrated force are the ones that present error. And strain energy percent error of the structure coincides with the percent error in the displacement in that direction. This indicates that bending is the dominant phenomena in terms of energy.
Also, the differences in the convergence rates indicate that, the solution is somehow sensible to the direction of the force and the presence of bending in one or more members.
Finally, the U-frame under multiple forces according to Figure 19 was analyzed, where F 1 = -500 N and F 2 = 1000 N are in the Y direction and F 3 = 1000 N is in the X direction. In this case, all members are under bending and axial solicitations.
The convergence criterion adopted was 1% error to the Euler-Bernoulli traditional element in both X and Y displacements of the junction where F 2 and F 3 were applied.
Solution converged using 60 Overhauser elements (20 per structural member), as illustrated in Figure 20, much more than both situations of the L-frame. Also, we can note that the error associated to the X-displacement is bigger, and this is the direction where bending displacements are more expressive, as it happened in the L-frame tests. Concerning the strain energy, again we attest that bending is the dominant phenomena, therefore the percent error in energy coincides the percent error in X-displacement (principal bending direction).
Finally, for all frame tests that involved bending, the same behavior of bigger errors in the nodal transverse displacements closer to the clamped edges is exhibited, as it happened in the clamped beam problems. Revealing that despite the cubic interpolation functions, the Overhauser solution is less accurate, especially for the clamped corner elements.

CONCLUSIONS
In this paper, we formulated C 1 beam finite element based on Overhauser interpolation and the Euler-Bernoulli beam theory that goes without the use of rotational degrees of freedom to represent bending.
The main objective of such interpolation scheme was to reduce the computational burden associated with the presence of the rotational degrees of freedom, but as observed in all the numerical tests performed, the Overhauser solution demands more computational effort because despite having only displacement degrees of freedom, more elements are required to achieve an acceptable result.
Also we have seen that solution is much more influenced by the essential and natural boundary conditions compared to the traditional Euler-Bernoulli elements, especially the corner elements. It is not a trivial question to be answered why this influence of the boundary condition occurs and why the cubic interpolation does not represent the solution as well as the Hermite cubics. But the article brings this important breakthrough in understanding such interpolation as this is still a difficulty to be overcome by the proposed interpolation class. A future work suggestion is to investigate the use of different element sizes for the extremities, to evaluate the influence of such elements in the behavior of the whole structure.
Finally, the application of such formulation is restricted to one and two-dimensional problems because torsion cannot be represented in the classical beam theory, without the use of a rotational degree of freedom, that in complex frame structures, couples with bending rotations.