BEAM ELEMENT UNDER FINITE ROTATIONS

The present work focuses on the 2-D formulation of a nonlinear beam model for slender structures that can exhibit large rotations of the cross sections while remaining in the small-strain regime. Bernoulli-Euler hypothesis that plane sections remain plane and perpendicular to the deformed beam centerline is combined with a linear elastic stress-strain law. The formulation is based on the integrated form of equilibrium equations and leads to a set of three first-order differential equations for the displacements and rotation, which are numerically integrated using a special version of the shooting method. The element has been implemented into an open-source finite element code to ease computations involving more complex structures. Numerical examples show a favorable comparison with standard beam elements formulated in the finite-strain framework and with analytical solutions.


Introduction
Highly slender fiber-or rod-like components represent essential constituents of mechanical systems in many fields of application such as civil, mechanical and biomedical engineering. It is widely recognized that slender bodies can be efficiently modeled applying a beam theory instead of a 3D continuum mechanics theory. Kirchoff proposed the first beam formulation which includes large three-dimensional deformations [1,2], and Reissner completed the theory for both two-dimensional [3] and three-dimensional cases [4] with two additional deformation measures representing the shear distortion of beam segments.
Reissner's finite-strain beam theory is one of the most important geometrically nonlinear models, subsequently extended and used by many other authors for 2-D and 3-D analysis of static as well as dynamic problems. Simo developed a dynamic formulation for Reissner's beam [5] and together with Vu-Quoc [6] initiated the finite element implementation. He also introduced the useful concept of a geometrically exact beam, based on recasting Reissner's theory in a form which is valid for any magnitude of displacements and rotations.
In this paper, a geometrically nonlinear beam model is formulated for the 2-D case. This model is applicable when the strains remain in a limited range while the rotations of beam sections can become arbitrarily large, and it properly accounts for the effect of curvature on the change of distance between end sections measured along the chord. Therefore the model considers large rotations and displacements while strains are treated as small and simple Hooke's law is used (an extension to a nonlinear material law would be relatively straightforward). Cross sections are still assumed to remain planar and perpendicular to the deformed beam axis. The formulation exploits the equilibrium equations in their strong form. They are combined with the kinematic equations and generalized material equations (linking internal forces to the curvature and centerline extension) and then numerically integrated.

Basic kinematic assumptions
In this section, the kinematics of the element is presented, leading to the strain-displacement equations. Each beam is analyzed in its local coordinate system, with the x-axis passing through the end joints in the undeformed configuration and the z-axis rotated by 90 • clockwise, while the rotations are positive counterclockwise (see Figure 1). We use u and w for horizontal and vertical displacements at an arbitrary point, while the displacements at the centerline are denoted by u s and w s .
Let us consider an infinitesimal segment of the centerline of initial length dx, which is after deformation transformed into a segment of length dx = ! (dx + du s ) 2 + dw 2 s . The corresponding centerline stretch can be expressed as where the prime denotes differentiation with respect to x. The assumption of planarity of the deformed section and of its perpendicularity to the deformed centerline implies that the rotation angle is linked to the centerline displacements by the relation and that the displacements of a general point with coordinate z can be expressed as In contrast to the geometrically linear beam model with u = u s + zϕ and w = w s , here the vertical displacement varies along the height of the section. However, the change of distance from the centerline due to deformation perpendicular to the beam axis is still neglected. Differentiating (3)-(4) with respect to x, we obtain and we can express the longitudinal stretch at an arbitrary location It is worth noting that the stretch is a linear function of the spatial coordinate z. Consequently, if we use a linear stress-strain law formulated in terms of the Biot strain defined as ε = λ − 1, the resulting work-conjugate stress will be a linear function of z, too. The choice of working with linear stress-strain relation between the Biot strain and stress can be justified in a regime where the strains are assumed to remain small. Making use of (7), we can express the Biot strain as where is the strain at the centerline and is the curvature.

Stress resultants and equilibrium
Based on (8) combined with Hooke's law σ = E ε, quantities ε s and κ that characterize the deformation of an infinitesimal beam segment can be linked to the normal force N and bending moment M by the standard linear relations in which EA is the sectional axial stiffness and EI is the sectional bending stiffness. As usual, E denotes Young's modulus, A is the sectional area and I is the sectional moment of inertia. The structure is assumed to be loaded only at its joints, i.e., no external loads are applied on individual beams, and so the internal forces in each beam can be expressed in terms of the end forces and moments, based on equilibrium conditions written for the deformed configuration. From equilibrium of the finite beam segment between the left end section and a generic section x, one can, according to Fig. 2, derive expressions in which X ab , and Z ab are the left-end forces (positive to the right and downwards) and M ab is the left-end moment (positive counterclockwise). For the sake of brevity, we will consider the moment as one of the (generalized) forces, and, in a similar spirit, the rotation as one of the (generalized) displacements.
Since the internal forces are involved only algebraically, they can be easily eliminated and equations (2) and (9)-(14) can be converted into the following set of first-order differential equations for unknown functions u s , w s and ϕ: Beam element under finite rotations Figure 2. Schematic representation of beam internal forces at a general cross section x and left-end forces and moment.
These equations are accurate even for very large rotations. The only limitation is that the strain remains small, which is true if ε s ≪ 1 and hκ ≪ 1 where h is the depth of the cross section.
In the context of a finite element analysis, the joint displacements are global degrees of freedom that are determined by equilibrium iterations on the structural level. In each iteration, one needs to determine, for each beam, the end forces that correspond to the given joint displacements and rotations. This means that constants X ab , Z ab and M ab in (15)-(17) are not known in advance. On the other hand, the displacements at both end sections are prescribed. At the left end, conditions can be considered as initial conditions for first-order differential equations (15)-(17). If the values of X ab , Z ab and M ab are somehow estimated, these equations can be integrated, at least numerically, and the values of ϕ(L), u s (L) and w s (L) can be determined. The values of forces X ab , Z ab and M ab then need to be adjusted such that the resulting kinematic quantities at the right end satisfy the remaining boundary conditions For a given set of end displacements, the initial estimate of the left-end forces can be constructed based on linear beam theory, or on the values at the end of the previous step if the calculation is done in the context of an incremental iterative structural analysis.
In fact, the suggested approach is a special version of the shooting method. Normally, the shooting method iterates on some of the initial values. Actually, the left-end forces would appear in static boundary conditions at the left end, which make it possible to integrate the differential equations of equilibrium and proceed from external forces to internal forces. Here, the equilibrium equations were directly presented in their integrated form (13)-(14), and the unknown variables play the role of integration constants instead of boundary values.

Numerical solution
The suggested numerical approach is based on the full form of equations (15)-(17). The right-end displacements u b , w b and ϕ b can be evaluated if the left-end displacements u a , w a and ϕ a and the left-end forces X ab , Z ab and M ab are given. To this end, the derivatives in (15)-(17) can be replaced by finite differences and the interval [0, L] is divided into N equally sized numerical segments.
Mapping of the right-end displacements, rotation, forces and moment based on the left-end displacements and rotation can formally be described by where and g is a column matrix of three functions of six variables. For given values of u a and u b , equation (24) represents a set of three nonlinear algebraic equations for unknowns f ab . The solution is found by the Newton-Raphson method, using the recursive formula where is the Jacobi matrix of mapping g with respect to variables f ab . The entries of the Jacobi matrix are evaluated numerically using the linearized version of the computational scheme. Once the left-end forces have converged, the right-end forces are easily evaluated from equilibrium conditions written for the whole beam. The foregoing equations were formulated in a local coordinate system aligned with the undeformed M Figure 3. Pure bending of a cantilever beam subjected to end moment.
beam. Since the deformed shape of the beam and the distribution of internal forces must be invariant with respect to rigid-body translations and rotations, one can also use a rotated coordinate systemx −z with the origin located at the left end of the beam in the deformed configuration (pointã in Fig. 2) and with thex axis in the direction of the tangent to the deformed centerline at the left end. With respect to this coordinate system, the left-end displacements are always zero, and so u a can be omitted from the list of arguments of function g. Of course, the rightend displacements u b must be computed from the global components of the joint displacements using an appropriate transformation rule, and the end forces obtained by the iterative process (26)-(27) must be transformed back into the global coordinate system. Consistent linearization of the relation between the end forces and end displacements leads to the element tangent stiffness matrix. Details will be presented in a full journal paper.

Numerical examples
A nonlinear beam element based on the proposed approach has been implemented into OOFEM [7,8], an object-oriented finite element code. To verify the implementation and demonstrate the potential of the suggested approach, two examples involving simple structures are shown next.

Pure bending of a cantilever beam
The first test is a cantilever of length L = 8 and bending stiffness EI = 1 loaded by a concentrated end moment M on its right-end. The exact solution to this problem is a circular arc with radius R = EI/M . An applied end moment, M = π/4, will force the rod to deform into a full closed circle. In this example, the moment is applied in six load steps, making the rod wind around itself at the end of the sixth step. The deformed shape of the beam at the end of each step is depicted in Fig. 3. The solution of the presented nonlinear model is compared with the one obtained by employing the geometrically exact finite beam element by Simo and Vu-Quoc [6] with a mesh of eight elements. The exact solution is reported as well. The overall agreement is good, and the simulation based on the present model, which uses only one two-noded element, is closer to the analytical solution.
The example demonstrates that the present model allows for a dramatic reduction of the number of global degrees of freedom, but of course the number of segments N used for numerical integration of the governing equations (15)-(17) must be chosen high enough to provide a good approximation. The presented results have been obtained using 100 segments. A close-up view of the sixth step circle is showed in Fig. 5 for calculations in which 6, 10, 20 and 50 numerical segments were employed. To ease the interpretation of the results we plot straight segments in between the integration points even though the curvature is constant along the beam. It turns out that 20 segments provide an acceptable solution, which is further improved and gets close to the analytical one if 50 segments are used. It is also noted that even with just 10 segments the method provides an accurate estimation of displacements at the integration points. To better quantify the previous comparisons, we calculated the sum of the distances between the the analytical and the computed solutions, normalized with respect to the number of sample points. The results for the sixth load step are reported in Table 1. Already for 8 integration segments, the error is smaller that 5%. The results indicate that comparable accuracy can be obtained if the present method replaces traditional finite elements by the same number of integration segments located within one single finite element.

Clamped-hinged circular arch
subject to point load  problem was investigated by many authors [6,9] and the exact solution based on the Kirchhoff-Love theory was given by DaDeppo and Schmidt [10]. The arch is circular with one boundary clamped and the other boundary hinged, and is loaded by a point load, as shown in Fig. 5. The cross section of the beam is a rectangle of depth h = 2.289 and width t = 1.0. The elastic modulus is set to E = 1.0 × 10 6 and the Poisson ratio is zero. Geometrical and constitutive parameters are dimensionless meaning that any length or force scale can be used for these values. Simo and Vu-Quoc [6] performed a FE analysis with a forty-element mesh, leading to a buckling load of 9.0528 EI/R 2 , while Wood and Zienkiewicz [9] calculated a buckling load of 9.24 EI/R 2 . The exact value of the buckling load reported by DaDeppo and Schmidt [10] is 8.97 EI/R 2 . In our simulations, we used forty elements as well, each with four integration segments, and we obtained a buckling load of 8.9874 EI/R 2 , which is very close to the analytical solution and exhibits a relative error less than 0.02%. The full comparison with the analytical solution is reported in Fig. 6 where we plot the load-displacement and the load-rotation curves relative to the top of the arch (loaded point). Again, the analytical solution is satisfactorily reproduced. The deformed configurations for the four load levels numbered in the curves are also shown together with the bending moments indicated by colors. The calculation was completed in 612 load steps using the arc-length method with a It is noted that Simo and Vu-Quoc [6] and Wood and Zienkiewicz [9] show a longer snap-through behavior together with the associated deformed shapes even though the response is not reasonable because of non-physical penetration of the left support into the structure. This aspect was investigated by Simo et al. [11], who showed that a contact constraint condition on the left support needs to be introduced to obtain a more realistic solution. For the aim of this paper we did not consider contact activation and the subsequent stiffening effect in the structure even though this extension could be added.

Conclusions
We present an efficient formulation for a geometrically nonlinear beam model that accounts for arbitrarily large rotations of the beam sections and the effect of curvature on the change of distance between end sections measured along the chord. In the first stage, the model is restricted to the small-strain framework and cross sections are still assumed to remain planar and perpendicular to the deformed beam axis. Equilibrium equations in their strong form are combined with the kinematic equations and generalized material equations and then numerically integrated using a special version of the shooting method.
Two examples show the capabilities of the element to replicate analytical results in a more accurate way with respect to existing beam finite elements that account for geometric nonlinearities. Moreover, the h-refinement of finite elements, which increases the number of global unknowns and is thus computationally demanding, can be substituted by employing a proper number of segments in the integration scheme, while keeping the number of elements limited. The element will be extended into a formulation that takes into account an initial curvature and also accounts for a possible contact mechanism activation.