## Abstract

This article presents a novel derivation for the governing equations of geometrically curved and twisted three-dimensional Timoshenko beams. The kinematic model of the beam was derived rigorously by adopting a parametric description of the axis of the beam, using the local Frenet–Serret reference system, and introducing the constraint of the beam cross ection planarity into the classical, first-order strain versus displacement relations for Cauchy’s continua. The resulting beam kinematic model includes a multiplicative term consisting of the inverse of the Jacobian of the beam axis curve. This term is not included in classical beam formulations available in the literature; its contribution vanishes exactly for straight beams and is negligible only for curved and twisted beams with slender geometry. Furthermore, to simplify the description of complex beam geometries, the governing equations were derived with reference to a generic position of the beam axis within the beam cross section. Finally, this study pursued the numerical implementation of the curved beam formulation within the conceptual framework of isogeometric analysis, which allows the exact description of the beam geometry. This avoids stress locking issues and the corresponding convergence problems encountered when classical straight beam finite elements are used to discretize the geometry of curved and twisted beams. Finally, this article presents the solution of several numerical examples to demonstrate the accuracy and effectiveness of the proposed theoretical formulation and numerical implementation.

## 1 Introduction

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.

## 2 Generalized Beam Formulation

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.

### 2.1 Geometry.

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\u2208[0,L]\u2192R3$ and *L* denotes the initial length of the curve.

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

**t**(

*s*) is the unit vector tangent to the beam axis and orthogonal to the cross section,

**n**(

*s*) is the normal unit vector, and

**b**(

*s*) is the binormal unit vector. These mutually orthogonal unit vectors form a local orthonormal basis $Q(s)={t,n,b}\u2208R3\xd73$, which is also assumed to provide the orientation of the cross section. At any given location of the beam axis, the cross section is identical in the local system of reference.

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*_{n}**n** + *p*_{b}**b**. 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).

**t**,

**n**,

**b**can be obtained as follows:

*τ*(

*s*) =

*d*

**n**(

*s*)/

*ds*·

**b**is the torsion of the curve.

### 2.2 Kinematics.

According to the beam assumptions, the displacement of a point at a generic cross section can be calculated as **u** = **u**_{0} + ** θ** ×

**p**, where

**u**

_{0}(

*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.

**J**is the Jacobian of the local to global transformation. According to Strang [33], one has

*J*is the Jacobian determinant,

*J*= 1 −

*κp*

_{n}. By virtue of Eq. (3), the small strain tensor in the global system of reference reads

*D*=

*τp*

_{b}, and

*E*= −

*τp*

_{n}.

while *ɛ*_{nn} = **n**^{T} · * ϵ* ·

**n**= 0,

*ɛ*

_{bb}=

**b**

^{T}·

**·**

*ϵ***b**= 0, and

*γ*

_{nb}=

**b**

^{T}·

**·**

*ϵ***n**+

**n**

^{T}·

**·**

*ϵ***b**= 0.

*in the local system of reference can be then contracted in a 3 × 1 vector*

**ϵ***with nonzero components as follows:*

**ε**

*ε*_{0}=

*d*

**u**

_{0}/

*d s*−

**×**

*θ***t**is the generalized beam strain vector and

**is the beam torsional/flexural curvature vector. Note that the derivation of Eq. (8) used the condition**

*χ**d p*

_{t}=

*ds*.

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 50*h*/*R* %, that is, for example, 5% for *h*/*R* = 0.1 and 50% for *h*/*R* = 1.

### 2.3 Equilibrium.

_{h}is the boundary with prescribed tractions. Since the variation of the external work has the form $\delta Wext=\u222bl(qt\delta u0t+qn\delta u0n+qb\delta u0b+mt\delta \theta t+mn\delta \theta n+mb\delta \theta b)ds$, the equilibrium equations at any given cross section can be written as follows:

### 2.4 Elastic Behavior.

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.

**f**=

**E**.

*η***f**= [

*N*,

*Q*

_{n},

*Q*

_{b},

*M*

_{t},

*M*

_{n},

*M*

_{b}]

^{T}is the stress resultant vector,

**= [**

*η**ɛ*

_{0tt},

*γ*

_{0tn},

*γ*

_{0tb},

*χ*

_{t},

*χ*

_{n},

*χ*

_{b}]

^{T}is the generalized strain vector, and

**E**is the sectional stiffness matrix, which reads

*α*

_{n}and

*α*

_{b}are the shear correction factors in the

**n**and

**b**local directions, respectively [35]. They take into account that the actual shear stress distribution cannot be uniform over the cross section, and they depend on the shape of the cross sections. The definitions in Eq. (14) are generalized versions of the cross-sectional properties (area, first-order, and second order area moments), which take into account, again, the effect of the local to global transformation via the term

*J*= 1 −

*κp*

_{n}. Finally, the beam stiffness matrix in Eq. (14) is not diagonal. Indeed, the equivalent first-order area moments

*S**

_{n}and

*S**

_{b}are not zero because the beam axis intersects the cross section in a point that, in general, is not its centroid. In addition, the equivalent mixed moment of inertia

*I**

_{nb}is nonzero because the two local axes

**n**and

**b**are not, in general, principal axes of inertia.

## 3 Isogeometric Implementation

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 $\Omega ^=[\xi 1,\xi m]\u2282R$ can be defined by specifying a knot vector with a nondecreasing order $\Xi ={\xi 1,\xi 2,\u2026,\xi m}$, an associated set of B-spline basis functions $NIp$ and a set of NURBS weights {*w*_{I}}, where *I* is the *I*th 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(\xi )=1$, if *ξ* ∈ [*ξ*_{I}, *ξ*_{I+1}[, otherwise $NI0(\xi )=0$.

*p*≥ 1, it can be defined recursively using the Cox-de Boor formula:

*p*= 0,

*N*

_{I,0}(

*ξ*) are piecewise constant functions; when

*p*= 1,

*N*

_{I,0}(

*ξ*) are the same basis functions of classical constant-strain finite elements. B-spline basis functions are linearly independent and have a partition of unity property, and their support is compact. However, they, in general, do not satisfy the Kronecker delta property [38].

*w*

_{I}are selected depending on the type of curve to be represented exactly. Note that when all weights

*w*

_{I}are equal to 1, the NURBS basis function reduces to the B-spline basis function, which can be seen as a subset of the NURBS basis function.

One then defines the nonzero entries in the knot vector $\Xi $ to span the parametric domain, $\Omega ^=[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 $\Omega ^e=[\xi I,\xi I+1](\xi I\u2260\xi I+1,I=p+1,p+2,\u2026,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 $\Omega ~=[\u22121,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 $\Omega ~$ to the elemental parametric domain $\Omega ^e$, $\phi ^e:\Omega ~\u2192\Omega ^e$ and the mapping from the parametric domain $\Omega ^$ to the physical domain Ω, $\phi :\Omega ^\u2192\Omega $ are assumed to be sufficiently smooth and invertible [39].

As mentioned earlier, considering a spatially curved beam in the physical domain $\Omega \u2282R3$, IGA requires a set of control points **P**_{I}, the corresponding weights of the control points *w*_{I}, a knot vector $\Xi =[\xi 1,\xi 2,\u2026,\xi I+p+1](I=1,2,\u2026,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.

^{e}∈ [

*s*

_{I},

*s*

_{I+1}], displacements and rotations read

It is worth noting that the smoothness condition for the classical Galerkin approach used in this study requires shape functions with only *C*^{0} continuity; this is typical of Timoshenko beam numerical implementations. However, the smoothness for the Frenet–Serret local basis requires *C*^{2} 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.

## 4 Numerical Examples

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.

### 4.1 Curved Cantilever Arch.

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*^{2}-norm relative errors of nodal displacements *u*_{1}, *u*_{2}, and nodal rotation *θ*_{3} versus the mesh size with quadratic and cubic NURBS basis functions are reported in Figs. 3(b)–3(d), respectively. The *L*^{2}-norm relative error can be calculated as follows: $\Vert \u25a0\u2212\u25a0h\Vert /\Vert \u25a0\Vert $, where $\u25a0h$ denotes the numerical values and $\u25a0$ 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*_{1}-displacements $u1A=u1,oriA\u22c5[Ewh3/(FR3)]$ and *x*_{2}-displacements $u2A=u2,oriA\u22c5[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*_{1} displacement; *x*_{2} 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.25*h*, 0.25*w*]. 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*_{1} displacement $u1tip$ and tip *x*_{2} 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\u22c5[Ewh3/(FR3)]$ and $u2tip=u2,oritip\u22c5[Ewh3/(FR3)]$, where $u1,oritip$ and $u2,oritip$ are the original *x*_{1} displacement and *x*_{2} 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.

### 4.2 Circular Balcony.

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*_{1}(*s*) = *R*cos(*s*/*R*), *x*_{2}(*s*) = *R*sin(*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*_{2} 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.

### 4.3 Helical Rod.

The next example is a helical rod subjected to a tip load. The helical rod has the expression *x*_{1}(*s*) = *a*cos(*s*/*c*), *x*_{2}(*s*) = *a*sin(*s*/*c*), *x*_{3}(*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*_{2} direction (Fig. 6(a)). The global vertical displacement *u*_{2} and rotation about *x*_{2}-axis *θ*_{2} 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.25*d*, 0.25*d*].

Figures 6(e) and 6(f) report the normalized, dimensionless tip displacement $u2tip=u2,oritip\u22c5[Ed4/(FH3)]$ and rotation $\theta 2tip=\theta 2,oritip\u22c5[Ed4/(FH2)]$ versus slenderness ratio with various beam axes, where $u2,oritip$, $\theta 2,oritip$, *E*, *d*, *F*, and *H* are the original *x*_{2} displacement, rotation around *x*_{2}-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.

### 4.4 Beams With Arbitrarily Positioned Beam Axis and Irregular Cross Sections.

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*=\u2211i=13Ai*$, $Sn*\u2211i=13Sni*$, $Ibb*\u2211i=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*_{1} = 0.8, *h*_{2} = 0.5, *h*_{3} = 0.3 m, *w*_{1} = 0.08, *w*_{2} = 0.05, *w*_{3} = 0.03 m, rotation angles *θ*_{1} = 1*π*/3, *θ*_{2} = 7*π*/8, *θ*_{3} = 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*_{2} and rotation around *x*_{2}-axis *θ*_{2} 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*^{2}-norm relative errors of nodal displacements *u*_{2}, and nodal rotation *θ*_{2} versus the mesh size are reported in Figs. 7(e) and 7(f).

## 5 Conclusion

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.

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

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.

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.

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

The numerical results are free of any stress locking issue.

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

## Acknowledgment

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.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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