Generalized Formulation for the Behavior of Geometrically Curved and Twisted Three-Dimensional Timoshenko Beams and Its Isogeometric Analysis Implementation

Curved and twisted beams are commonly used in many applications in civil, mechanical, and aerospace engineering due to their esthetics and unique load-bearing properties. Tall buildings with curved and twisted columns have been designed and constructed in many parts of the world in recent years [1–3]. This type of columns not only leads to stunning building façades but they are also efficient in resisting both gravity and lateral loads. In contrast, straight columns are in most cases designed to only resist gravity loads. Wind turbine blades and helicopter blades, which are commonly found in the energy industry and aerospace engineering, can be also modeled as beam-like structures [4–6]. Straight beam models have been used in the past in many of the dynamic and stability analyses of blades. However, the continuous effort on design optimization of the aerodynamic and structural performances of blades makes the beam geometry more complex; hence, analytical methods for curved and twisted beams have become increasingly prevalent. Geometrically curved and twisted smart beams that can sense and respond to stimuli also gained attention recently [7,8]. The need for analytical capabilities for smart beams with curved and twisted geometry has inspired many studies, including piezoelectric and multiphysical behavior of smart beams [8,9].

Analysis methods for beams with increasing geometric complexities have been extensively studied by several authors in the past. Reissner [10] presented a variational analysis of small deformations of pretwisted elastic beams. Sandhu et al. [11] and Crisfield [12] developed co-rotation formulations for a curved and twisted beam element. Simo and Vu-Quoc [13] developed a geometrically exact beam model including shear and torsion warping deformations. The limitation of all the published formulations is that the kinematic model that relates the strains at one point of a beam cross section with the beam axis elastic deformation and elastic curvature is assumed and only valid for slender geometries, as opposed to rigorously derived from the continuum definition of strains. In addition, these formulations assume the axis of the beam to coincide with the centroid of the cross section and the local system of reference to be the principal axes of inertia. This is convenient for analytical hand calculations, but it is instead cumbersome in computational analysis because the cross-sectional geometrical properties need to be calculated before defining the beam axis. This is not convenient in complex cases.

The classical finite-element formulation of beam theories uses straight beam elements, in which the axial behavior is decoupled from the transverse behavior. However, by using straight finite elements to approximate a curved beam, locking issues arise from the interplay of shear and membrane behaviors. This leads to a spurious stiffer response and an overestimation of shear stresses. The fundamental underlying issue is that the axial and transverse behaviors are not decoupled in the actual curved beam [14–17]. A solution to this issue is to exactly describe the beam geometry via isogeometric analysis (IGA).

Starting from the pioneering work of several researchers, e.g., Kagan et al. [18], Rogers [19], and Hughes et al. [20], isogeometric analysis (IGA) uses nonuniform rational B-splines (NURBS) basis functions to represent both the geometry and the field variables. Among the studies of IGA in structural mechanics, shell element and rod element formulations are frequently discussed. These include the work of Kiendl et al. [21], Benson et al. [22], Echter et al. [23], Auricchio et al. [24], Hu et al. [25], and Weeger et al. [16]. The structural analysis of beams, especially those with complex geometries, can be accurately performed with the help of IGA, while the computational cost is significantly reduced compared to IGA with solid elements. The isogeometric beam element formulation of curved beams has been presented for both two-dimensional and three-dimensional cases and for both Euler–Bernoulli beam and Timoshenko beam in Refs. [15,26,27]. Locking issues as well as the locking free formulations of curved beams are also discussed in the literature and can be found in Refs. [14,28,29]. Nonlinear analysis of isogeometric curved beams gain more attention nowadays and are discussed in Refs. [30,31], among others.

The underlying assumptions for the new beam formulation are the same as those made in classical Timoshenko beam theory: (1) the beam axis is orthogonal to the beam cross sections before the deformation; (2) the cross sections remain planar and preserve their shape and size during deformation; and (3) displacements and rotations are small compared to the beam size (first-order theory). The warping effects of the section planes are neglected in this work. The authors recognize that warping effects might be important particularly for open thin-walled cross sections, but they leave this additional complexity to future work.

The geometry of a curved and twisted beam can be represented by the mathematical description of the beam axis and its cross sections. The generic position, r(s), of a point on the beam axis can be expressed as a function of the arc-length s, where $s∈[0,L]→R3$ and L denotes the initial length of the curve.

The vector r(s) allows calculating the Frenet–Serret local basis as follows:

The position of any generic point P on a given cross section centered at r(s) is calculated as x(s, p_n, p_b) = r(s) + p = r(s) + p_nn + p_bb. The out-of-plane component of p is zero, p_t = 0, because the cross section is orthogonal to the beam axis in the undeformed configuration (Fig. 1).

According to the beam assumptions, the displacement of a point at a generic cross section can be calculated as u = u₀ + θ × p, where u₀(s) = [u_0t, u_0n, u_0b]^T is the cross-sectional translation, θ(s) = [θ_t, θ_n, θ_b]^T is the cross-sectional rotation with reference to point O corresponding to the intersection between the axis and the cross section (Fig. 1). Point O is any point in the cross section and it does not need to be the cross section centroid.

while ɛ_nn = n^T · ϵ · n = 0, ɛ_bb = b^T · ϵ · b = 0, and γ_nb = b^T · ϵ · n + n^T · ϵ · b = 0.

Equation (8) differs from the strain definition in classical Timoshenko beam formulations, which do not have the multiplier term 1/J = 1/(1 − κp_n).

One has J = 1 for a straight beam (κ = 0) and J ≈ 1 if κh ≪ 1, where h is the characteristic size of the cross section. However, the effect of J on the local strains cannot be neglected for large values of κh, which occurs in the case of stocky geometries. The definition of κh as the curviness of the beam was first introduced by Borković et al. [34]. Equation (8) leads to cross-sectional strain profiles that are nonlinear. From a physical point of view, this is due to the fact that material fibers away from the geometrical center of curvature are longer than materials fibers closer to the radius of curvature in their undeformed configuration. For a circular beam of radius R with a rectangular cross section of depth h, the error in the strain calculation without the curvature effect is 50h/R %, that is, for example, 5% for h/R = 0.1 and 50% for h/R = 1.

In the linear elastic regime, one can write the stresses as σ_tt = Eɛ_tt, τ_tn = Gγ_tn, and τ_tb = Gγ_tb, where E is the elastic modulus, G = E/(2 + 2ν) is the elastic shear modulus, and ν is Poisson’s ratio.

Following Hughes et al. [20], this study employs NURBS (nonuniform rational B-spline) as shape functions to interpolate both the beam geometry and the unknown fields. This technique is known in the literature as isogeometric analysis (IGA). The main advantage of IGA is the accurate and sometimes exact representation of the geometry: This is a critical aspect for the simulation of spatially curved and twisted beams. Furthermore, a unique advantage of IGA compared to the classical finite-element (FE) method is the possibility of global regularity refinement, which enables high-order interpolation of unknown fields without significantly increasing the computational cost [20,36,37].

A NURBS basis function on the parametric domain $Ω^=[ξ1,ξm]⊂R$ can be defined by specifying a knot vector with a nondecreasing order $Ξ={ξ1,ξ2,…,ξm}$⁠, an associated set of B-spline basis functions $NIp$ and a set of NURBS weights {w_I}, where I is the Ith knot, n is the number of basis functions, and p is the polynomial order. In IGA, the relation m = n + p + 1 always holds. The B-spline basis function $NIp$ can be constructed starting from p = 0 with $NI0(ξ)=1$⁠, if ξ ∈ [ξ_I, ξ_I+1[, otherwise $NI0(ξ)=0$⁠.

One then defines the nonzero entries in the knot vector $Ξ$ to span the parametric domain, $Ω^=[0,1]$ if normalized. The element after spatial discretization in the parametric domain now can be defined as a span of the unique entries of the knot vector $Ω^e=[ξI,ξI+1](ξI≠ξI+1,I=p+1,p+2,…,ns)$⁠, where n_s is the number of unique knots.

Another domain that is commonly used for numerical quadrature is referred to as the parent domain $Ω~=[−1,1]$⁠. It is worth mentioning that the parent domain in IGA is always referred to as the parametric domain in conventional FE formulations, and the parametric domain used in IGA is absent in the FE context. The parametric domain is essentially an additional domain in IGA, and hence, an additional mapping is needed. Figure 2 illustrates the spatial mapping from the parent domain to the physical domain via the parametric domain. The mapping from the parent domain $Ω~$ to the elemental parametric domain $Ω^e$⁠, $φ^e:Ω~→Ω^e$ and the mapping from the parametric domain $Ω^$ to the physical domain Ω, $φ:Ω^→Ω$ are assumed to be sufficiently smooth and invertible [39].

As mentioned earlier, considering a spatially curved beam in the physical domain $Ω⊂R3$⁠, IGA requires a set of control points P_I, the corresponding weights of the control points w_I, a knot vector $Ξ=[ξ1,ξ2,…,ξI+p+1](I=1,2,…,n)$⁠, the number of control points n, and the polynomial order p. This information is commonly found in most CAD software applications and packages and must be imported before the analysis.

It is worth noting that the smoothness condition for the classical Galerkin approach used in this study requires shape functions with only C⁰ continuity; this is typical of Timoshenko beam numerical implementations. However, the smoothness for the Frenet–Serret local basis requires C² continuity. Since the NURBS basis function $RIp(s)$ is C^p−k continuous, at least p = 2 degree shape functions are needed to exactly capture the geometry of the beam.

Finally, by using the weak form of the equilibrium equations, one can compute the element stiffness matrix and nodal load vectors as customarily done in Galerkin FE implementations [27,39].

To verify the proposed beam formulation, numerical examples of 3D beams with various geometrical complexities are presented in this section. Three different geometries are included: (1) a curved cantilever arch, (2) a circular balcony, and (3) a helical rod. They all represent respective complexities in terms of geometry and boundary conditions. One additional numerical example of a curved cantilever arch with a cruciform cross section is provided as well, to investigate the capability of using the new beam formulation for beam problems with irregular cross-sectional shapes.

The first example is a cantilever quarter circle arch subjected to an in-plane tip load. The geometry of the quarter circle arch axis can be categorized as an in-plane curve with a constant curvature κ and zero torsion τ = 0 along the arc-length. The quarter circle arch of curvature radius R has a rectangular cross section with the dimensions of h × w. The curved arch is clamped at one end and loaded at the other end with a concentrated force F pointing toward its curvature center (see Fig. 3(a)).

A representative convergence study for the classical beam formulation (1/J = 1) with a slenderness ratio h/R = 0.1 is first performed, to investigate the convergence properties of IGA beam simulations using both the standard h- (mesh size) and p- (degree of basis functions) refinements. The L²-norm relative errors of nodal displacements u₁, u₂, and nodal rotation θ₃ versus the mesh size with quadratic and cubic NURBS basis functions are reported in Figs. 3(b)–3(d), respectively. The L²-norm relative error can be calculated as follows: $‖■−■h‖/‖■‖$⁠, where $■h$ denotes the numerical values and $■$ denotes the reference values reported in Cazzani et al. [15]. It can be observed that higher degrees of the basis functions lead to higher convergence rates, as well as more accurate results.

The influence of the multiplier term 1/J in the new beam formulation is then investigated by comparing the beam simulations of the new beam formulation with those of the classical beam formulation (1/J = 1) and those of 3D solid finite elements. Beams with slenderness ratios h/R ranging from 0.1 to 1.0 were simulated. Figures 3(e) and 3(f) report the normalized, dimensionless x₁-displacements $u1A=u1,oriA⋅[Ewh3/(FR3)]$ and x₂-displacements $u2A=u2,oriA⋅[Ewh3/(FR3)]$ at point A on the edge center of the tip cross section (see Fig. 3(a)), respectively, where $u1,oriA$⁠, $u2,oriA$⁠, E, w, h, F, and R are the original x₁ displacement; x₂ displacement at point A; beam elastic modulus; cross-sectional width, height, and magnitude of applied load; and curvature radius, respectively. The new beam formulation and classical beam formulation results were obtained with 16 IGA beam elements with cubic NURBS basis functions; the 3D finite-element solution was calculated by using 1024 × 16 × 16 solid finite elements. It is worth noting that the results of the new beam formulation are relatively close to the reference 3D FE results, while the classical beam formulation with 1/J = 1 is inadequate to accurately simulate the beam deflections. This is particularly true for slenderness ratios h/R > 0.5 (thick beams).

The difference between the beam solutions and the 3D finite-element solution is due to two limitations of the Timoshenko beam theory: (1) the higher the slenderness ratio is, the harder the shape of the beam sections can be approximated by a plane, and the planar integration used in sectional stress calculations is not accurate anymore; (2) the change in the reference length in strain calculations is more significant for higher slenderness ratio cases. While the new beam formulation adopts the multiplier term 1/J to resolve the second issue, the classical beam formulation basically has no mitigation for any of the issues mentioned earlier.

Another set of simulations was conducted for beams with the same geometry but with arbitrary positions of the beam axis. Figures 4(a)–4(c) show curved arches that are simulated with the beam axis located at the center, top, and bottom of the cross section, respectively; a diagram of all the locations of the beam axis considered in this comparison is shown in Fig. 4(d); and the local coordinates of the generic point (denoted “X” in Fig. 4(d)) are: [−0.25h, 0.25w]. It is worth noting that for the cases with the beam axis that does not intersect the cross section at the centroid location, a bending moment M = d × F must be applied, where d is the axis eccentricity. This ensures that the solved physical problem is always the same, i.e., the cantilever circular arch with radial load applied at the centroid of the cross section. Figures 4(e) and 4(f) report the normalized tip x₁ displacement $u1tip$ and tip x₂ displacement $u2tip$ versus slenderness ratio with different positions of the beam axis. Similar to the results shown in Figs. 3(e) and 3(f), the dimensionless, normalized displacements are calculated as follows: $u1tip=u1,oritip⋅[Ewh3/(FR3)]$ and $u2tip=u2,oritip⋅[Ewh3/(FR3)]$⁠, where $u1,oritip$ and $u2,oritip$ are the original x₁ displacement and x₂ displacement at the centroid at the free-tip section, respectively. The overlapped results of variously positioned beam axes in Figs. 4(e) and 4(f) show that the new beam formulation can account for the effect of changing positions (and hence changing reference length in strain calculations) of the beam axis on the beam computations, which can be considered as one of the advantages of the new beam formulation over the classical beam formulation as the location of the beam axis can be arbitrarily selected within the cross section.

The second example is a semi-circular balcony subjected to an out-of-plane distributed load. The geometry of the circular balcony can be described by the expression x₁(s) = Rcos(s/R), x₂(s) = Rsin(s/R), where R is the radius of the curvature and s is the arc-length. The dimensions of the circular balcony are selected to be consistent with the dimensions R = 3 m, h = 0.3 m, and w = 0.3 m of a numerical example in the study by Zhang et al. [27]. The semi-circular structure was clamped at both ends; a uniformly distributed load q = 5 kN/m was applied in the negative x₂ direction (Fig. 5(a)). After a convergence study, a mesh of 32 elements with the cubic basis functions was selected. The calculated local displacement u_b, the local rotation about t-axis θ_t, and the local rotation about n-axis θ_n versus the arch length s with the aforementioned mesh are compared with the values in the study by Zhang et al. [27] and are reported in Figs. 5(b)–5(d), respectively. Because the slenderness ratio of the curved arch (h/R = 0.1 for this example) is small, the differences between the results calculated by the new beam formulation and those calculated by the classical beam formulation are negligible. An excellent overall agreement shows that the current formulation has high accuracy with a relatively few number of elements and low degrees of the basis functions.

The next example is a helical rod subjected to a tip load. The helical rod has the expression x₁(s) = acos(s/c), x₂(s) = asin(s/c), x₃(s) = b s/c, where a = 2, b = 3/2π and $c=a2+b2=2.06$⁠, the beam axis has a curvature radius of 2 m, a total height of H = 3 m, and can be categorized as a 3D structure with a constant curvature κ and torsion τ along the arc-length. The cross section is circular with a diameter d, which is constant along the arc-length. Varying diameters d were selected to make the slenderness ratios equal to d/H = 0.33, 0.1, 0.05, 0.033, and 0.01, respectively. The curved beam is fixed at one end and loaded at the other end with concentrated force F = 10 kN in the negative x₂ direction (Fig. 6(a)). The global vertical displacement u₂ and rotation about x₂-axis θ₂ versus the arch length s of the beam axis for slenderness ratio d/H = 0.05 are as shown in Figs. 6(b) and 6(c), respectively. The comparison of the tip displacement and rotation for the beams with arbitrarily positioned beam axis is shown in Figs. 6(d)–6(f). Figure 6(d) shows the positions of the beam axis in this comparison, the local n − b coordinates of the generic point (denoted “X” in Fig. 6(d)) is: [0.25d, 0.25d].

Figures 6(e) and 6(f) report the normalized, dimensionless tip displacement $u2tip=u2,oritip⋅[Ed4/(FH3)]$ and rotation $θ2tip=θ2,oritip⋅[Ed4/(FH2)]$ versus slenderness ratio with various beam axes, where $u2,oritip$⁠, $θ2,oritip$⁠, E, d, F, and H are the original x₂ displacement, rotation around x₂-axis at the centroid at the free-tip section, beam elastic modulus, cross-sectional diameter, magnitude of applied load, and total height of the beam, respectively. Again, the overlapping results of the helical rod show that the new beam formulation can accurately simulate beam deflections with arbitrarily selected positions of the beam axis.

One additional numerical example is provided to demonstrate the possibility of simulating beams with irregular cross sections with the generalized beam formulation. A cantilever quarter circle arch with a “tri-webs” cross section subjected to an in-plane tip load was simulated (see Fig. 7(a)). A clamped-free boundary condition was used, and a tip concentrated force F acted at the free end toward the curvature center. The shape of the cross section can be approximately seen as an assembly of three rectangles with dimensions w_i × h_i (i = 1, 2, 3), the beam axis passes through the mid-point of the bottom edge of each rectangle, each rectangle rotates counter-clockwise around the beam axis with angle θ_i within the local coordinate system n − b, and the overlapped area can be neglected if one assumes w_i ≪ h_i, as shown in Fig. 7(b). The sectional properties of the “tri-webs” cross section can be calculated by taking the superposition of those properties of each web, i.e., $A*=∑i=13Ai*$⁠, $Sn*∑i=13Sni*$⁠, $Ibb*∑i=13Ibbi*$⁠, etc. The shear coefficient has no general estimation for the irregular cross sections, but it can always be evaluated by the ratio of the average shear strain on a section to the shear strain at the shear center. After calculation, approximate shear coefficients α_n = 0.4 and α_b = 0.35 are used.

The beam dimensions in this numerical example are: radius of curvature R = 5 m, web dimensions h₁ = 0.8, h₂ = 0.5, h₃ = 0.3 m, w₁ = 0.08, w₂ = 0.05, w₃ = 0.03 m, rotation angles θ₁ = 1π/3, θ₂ = 7π/8, θ₃ = 13π/8. The material properties used are as follows: elastic modulus E = 200 GPa and Poisson’s ratio ν = 0.3. The applied tip load was F = 10 kN. Because of the absence of the reference solutions, the results of the IGA beam simulation with the finest mesh (1024 elements) and the highest degree of the basis functions (6th degree) are used as the reference solution. The initial and deformed shapes of the circular arch corresponding to the reference solution are shown in Fig. 7(a), and the deformation is multiplied with the scale factor 100. It can be observed that the in-plane load F leads to not only the in-plane bending of the beam but also in the out-of-plane bending and the torsion around the beam axis, and this reflects the fully coupled behaviors of the beam with an irregular cross section. The displacement u₂ and rotation around x₂-axis θ₂ along the arch length s of the beam axis are shown in Figs. 7(c) and 7(d), respectively. With the reference solutions, the convergence studies of the L²-norm relative errors of nodal displacements u₂, and nodal rotation θ₂ versus the mesh size are reported in Figs. 7(e) and 7(f).

In this study, a new generalized Timoshenko beam formulation was developed to accurately capture the deformation of geometrically curved and twisted beams. The proposed beam formulation employs a parameterization of the beam axis with its arc-length and a local system of reference described by the Frenet–Serret basis. Furthermore, a beam kinematic model, more accurate than the ones currently available in the literature, is derived rigorously imposing the kinematic constraints dictated by the Timoshenko beam assumptions. Compared to existing formulations, the derived kinematic model features the effect of the initial curvature of the beam via a multiplicative term and leads to a nonlinear distribution of strains over the cross section. The resulting theory was implemented using isogeometric analysis and was used to solve four examples with various degrees of complexity.

From the obtained results one may draw the following conclusions.

1. The generalized Timoshenko beam formulation presented in this article allows the seamless analysis of spatially curved and twisted beam geometry.

2. The beam geometry can be directly imported and used from CAD software packages without the need of any preprocessing including precalculation of cross section centroids and/or principal axis of inertia.

3. The axis of the beam can intersect the cross section at any generic point of the cross-sectional plane. This simplifies the analysis of beams with complex cross sections.

4. The IGA implementation of the proposed formulation leads to optimal convergence.

5. The numerical results are free of any stress locking issue.

6. The obtained results are more accurate than the ones obtained with classical Timoshenko beam for a wide range of slenderness ratios.

The presented work was supported in part by the US National Science Foundation through grant CMMI–1762757. Any opinions, findings, conclusions, and recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the sponsors.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.