On the geometrically exact low-order modelling of a flexible beam: formulation and numerical tests.

This paper proposes a low-order geometrically exact flexible beam formulation based on the utilization of generic beam shape functions to approximate distributed kinematic properties of the deformed structure. The proposed nonlinear beam shapes approach is in contrast to the majority of geometrically nonlinear treatments in the literature in which element-based-and hence high-order-discretizations are adopted. The kinematic quantities approximated specifically pertain to shear and extensional gradients as well as local orientation parameters based on an arbitrary set of globally referenced attitude parameters. In developing the dynamic equations of motion, an Euler angle parametrization is selected as it is found to yield fast computational performance. The resulting dynamic formulation is closed using an example shape function set satisfying the single generic kinematic constraint. The formulation is demonstrated via its application to the modelling of a series of static and dynamic test cases of both simple and non-prismatic structures; the simulated results are verified using MSC Nastran and an element-based intrinsic beam formulation. Through these examples, it is shown that the nonlinear beam shapes approach is able to accurately capture the beam behaviour with a very minimal number of system states.

Introduction

Within the field of flexible structural modelling, one-dimensional beam models display a surprisingly wide applicability to the representation of numerous structural applications such as the modelling of aircraft wings, flexible satellites, turbine blades, communication towers, energy harvesters, atomic force microscopes and DNA molecules, to name but a few. Such treatments are typically suited to relatively slender structures with a single characteristic large dimension and, together with an appropriate constitutive material law, a one-dimensional representation may capture many aspects of the full three-dimensional problem. The same slender properties that make many physical structures amenable to beam representation often result in significant structural flexibility and, particularly for systems with free boundary conditions, this can lead to large deflections and subsequent geometric nonlinearity. Therefore, many studies treat such nonlinearity in kinematic frameworks permitting arbitrarily large rotation of either the structure directly, or to element frames in which structural sub-domains are treated. Of particular interest to the present study is the class of geometrically exact modelling methods first proposed in the seminal work of Reissner [1]. These formulations are centred around the concept of a local body fixed coordinate system in which the material deformation of the beam is cast. The kinematic description of the beam follows from knowledge of this local material deformation, as well as the span-varying orientation of the local body system, which may comprise a large rotation from the initial undeformed configuration; consequently, a global formulation of the flexible beam is attained in which geometric nonlinearity is inherently captured.

In a majority of cases, such geometrically nonlinear formulations are implemented using an element-based discretization of the structure. In [2,3], Simo and Vu-Quoc present a geometrically exact parametrization based upon the vector components of the global rotation and displacement from the initial undeformed configuration. This approach yields a geometrically nonlinear finite-element formulation of the flexible structure in terms of rotation and displacement states. Bauchau et al. [4] compare this formulation with the absolute nodal coordinate formulation of Shabana et al. [5] that uses gradient information rather than a rotational field in expressing the deflected geometry of the system. With a sufficiently fine discretization, linear strain–displacement laws may be used within each element frame. The class of co-rotational formalisms (first proposed in [6-8]) takes advantage of this by combining small local material deformation with arbitrarily large total rotations of the structure via the treatment of an element fixed (co-rotational) reference frame. Cesnik & Brown [9] allow for larger local deformations by accounting for the non-small linear variation of curvature within each element frame; the resulting formulation is cast directly in terms of curvature rather than displacement states. Consequently, the larger deformation allowed within each element frame affords the use of fewer elements in capturing geometric nonlinearity at the cost of increased complexity and bandwidth of the system matrices. Santos et al. [10] employ stress-resultant and displacement states in formulating a complimentary energy-based method for the representation of large static deformations. Together with additional displacement constraints, this allows for the treatment of nonlinearly deformed frame structures, free from the shear-locking phenomenon [11] and allowing for the accurate recovery of distributed stress. However, this comes at the cost of 24 degrees of freedom per element (see [12] for an overview of past research efforts and ongoing challenges in developing complementary energy principles for large deformation problems). Hodges [13] introduces the intrinsic beam formulation wherein constituent relations and the subsequent equations of motion are cast purely in a span-varying body-referenced intrinsic coordinate system. Patil et al. [14] use an element-based implementation of this formulation in studying the static aeroelastic response of a simple high-aspect ratio aircraft. Here, the intrinsic reference system effectively becomes the co-rotational system for each element and global geometric information is recovered via quaternion integration from root to tip. The review papers of Wasfy & Noor [15] and Shabana [16] further capture the large body of flexible multibody dynamic modelling literature.

A key feature of the above element-based descriptions is the generality with which structural configurations may be represented; however, this comes at the expense of a large number of system states. To reduce the problem size, one approach is to use model reduction methods; here one projects the full system onto a subset of modes using a regression analysis (for example, via implicit condensation [17,18] or enforced displacement [19] tests; see the review of Mignolet et al. [20] for further details) to fit the coefficients of higher-order terms, thus capturing nonlinear effects. However, the degree of nonlinear truncation required to render the regression problem tractable, as well as the inherent errors introduced when approximating the transient dynamics of the system, limits the domain of validity over which reduction techniques suitably approximate the full system.

For the current study, a different approach is adopted. Rather than reducing a high-order (large number of states) element-based model, the focus is cast instead on the low-order (minimal state) approximation of the full geometrically nonlinear dynamic problem. To this end, a nonlinear beam shapes approach is employed which approximates the kinematic quantities distributed along the full length of the flexible structure without subdivision into smaller elements. Compared to a typical element-based representation the reduction in system states can be significant. Such an approach is typical for the treatment of linear systems via, for example, Rayleigh Ritz, Galerkin and Modal type approximations (see standard texts on vibration modelling, e.g. [21]). In a geometrically nonlinear context, Patil & Althoff [22] consider higher-order shape functions as applied to the intrinsic beam formulation of [13]; here a set of shifted Legendre functions are used in approximating local angular velocity, translational velocity, force and moment quantities, expressed in the intrinsic reference frame. The resulting formulation is applied to the modelling of a flexible beam; however, a number of simplifying assumptions are required (i.e. treatment of a prismatic, constant property beam) in rendering the problem tractable to analysis. In the current paper, a geometrically exact flexible beam formulation is presented using a shape-based approximation of the kinematic beam quantities to yield a low-order continuous formulation of the problem. By using a global referencing attitude representation of the beam orientation at any spanwise location, analyticity is retained up to and including the integrands of positional quantities and their derivatives, thus lessening the computational requirements on the numerical integration of terms. The resulting dynamic equations are suitable for the representation of beam problems displaying generic span-varying geometric and material properties that typically exist in real structures; the underlying shape functions are kept general with candidate sets requiring only the satisfaction of a single kinematic constraint.

This paper is divided into two parts detailing firstly the modelling approach, followed by examples of its application. In §2, the mathematical formulation of the dynamic equations is presented, beginning with the variational statement of the geometrically exact problem and then introducing the chosen parametrization. The generic requirements and example application of an approximating shape basis is detailed in §3. In §4, a variety of validating test cases, which exercise various aspects integral to the treatment of a typical flexible beam problem, are presented. Throughout these test cases, additional results are generated as a baseline for comparison using both MSC Nastran (a well-established finite-element package capable of treating geometric nonlinearity) and the intrinsic beam formulation of [13] implemented as a finite-element code [23,24]. The latter modelling approach forms a particularly suitable base of comparison as, in addition to constituting a well-documented and validated nonlinear treatment, the resulting equations of motion take the form of a high-order system with fairly sparse system matrices and maximum order 2 nonlinearity. This serves as an interesting counterpart to the proposed formulation in this study which admits a greatly reduced number of states with fully populated form and a high degree of nonlinearity. Following the demonstration of the proposed nonlinear beam shapes approach in treating the discussed examples, conclusions are drawn in §5.

Presentation of dynamic equations

This section presents the derivation of the equations of motion describing a generic flexible nonlinear beam. This derivation proceeds with a geometrically exact beam treatment [2,3,25] wherein consideration of a local intrinsic coordinate system is used to capture the full geometric nonlinearity inherent in this class of flexible beam problem. By relating generic shape function sets to a spanwise-varying kinematic vector—consisting of three attitude, one extension and two shear parameters—a minimal state representation of the system is developed here which yields fast computational performance, easy truncation of shear and extensional deformations and admits representation as a nonlinear system of ordinary differential equations.

In §2a,b, a general statement of the geometrically exact variational problem is given wherein virtual work contributions are presented in terms of incremental displacement, shear, extensional and curvature deformations of the beam. In §2c, an Euler angle attitude representation is used to link these variational terms to a minimal set of shear, extension and rotational parameters. Following this, in §2d, a shape function discretization is used to split these parameters into spatial and temporal components yielding a minimal state representation of the problem in terms of the set of time-dependent generalized states from which the equations of motion follow.

Kinematic description

Consider first the flexible beam of length L depicted in figure 1, deemed to be a suitable representation of some slender three-dimensional structure, reduced about an arbitrary material reference line that forms the beam axis; the curvilinear coordinate s∈[0, L] is used to parametrize the distance along this axis. The axis itself is defined by the position vector (s) taken in some arbitrary reference system denoted A (here considered with origin placed at the beam root). This system is in turn embedded in the inertial reference system G with positional offset . Now take the intrinsic body-fixed system I (depicted at an example location s along the beam) represented by the orthonormal vector triad {, , }. In the undeformed configuration, lies parallel to the beam axis and and are mutually perpendicular directions in the cross-sectional plane normal to . As the beam deflects the orientation of and (in the absence of cross-sectional warping) follow the evolution of their respective material lines while remains normal to the cross-sectional plane they span. In figure 2, example orientations of , and ′(s) (where the superscript •′ denotes differentiation with respect to s) are shown for the cases of planar pure bending (a) and pure shear (b). The difference in cross-sectional orientation for each sub-case illustrates the important point that, in the presence of shear deformation, ≠′. Therefore, in general, the reference line ′(s) comprises a spanwise integration along (including axial extension) in addition to shear deformations in the and directions. [(s) is thus given the form: (cf. Reissner [1]), where ε(s), τ(s) and τ(s) denote axial extension, shear and shear deformations of the beam at each spanwise location, respectively. Each subscript in square brackets denotes the reference system (G, A or I depicted in figure 1) in which each vector quantity is cast. These coordinate systems are related by the rotation matrices [] and [] such that for any vector

Figure 1.

Reference coordinate systems and position vectors used in characterizing the three-dimensional beam geometry.

Figure 2.

Illustration of the distinction between (a) bending and (b) shear-type structural deformation.

Virtual work terms

Following on from the kinematic description, the equations of motion describing the geometrically exact beam are now developed in weak form; individual terms are formulated via the principle of virtual work done by internal and external forces acting over incremental displacements, rotations and beam strains. These work terms are detailed in the following subsections and in each case, the variational quantities are related back to the intrinsic coordinate system depicted in figure 1.

Material stress

In treating the stress terms of the system, strain deformations of the structural material are considered as lumped translational and rotational gradients along the parametric beam line (s). Hence the deformation vector (s) is introduced here which consists of the components

Here, κ(s), κ(s) and κ(s) give the , and components of the beam curvature ′′(s), respectively. The curvature may be written as[1] often expressed in the compact form where denotes the skew symmetric matrix formed from the vector  = (κ, κ, κ).

With these curvatures defined, the strain-induced internal force and moment distribution may be written as the six-component force/moment vector where denotes some minimal stress pre-curvature of the beam and the generic stiffness matrix is written as For the examples of this study, the linear stress–strain relationship is assumed wherein diag() = (GA, EA, GA, EI, GJ, EI) and all other entries are zero. Note however that, more generally, one may simply consider to be some arbitrary constitutive law relating the internal stress distribution to deformations Δ about the beam reference axis (s) [25]. For structures with simple geometries and known material properties, this may be obtained analytically by relating the beam deformation vector to the material strain tensor (e.g. [26]); for more complex structures, the constitutive relation may be obtained from physical or modelled tests of the full structure followed by the subsequent one-dimensional reduction about an arbitrary reference line [27]. In this general case of equation (2.8), the virtual work done over the beam, given the incremental strain deformation , is

Structural damping

Formulation of the virtual work done by structural dissipation follows in much the same form as the preceding strain term. Considering again the variation , one may specify the function relating the internal structural damping force/moment to the strain rates . Thus, the virtual work for this incremental deformation takes the form

Applied loads

External forces. The work done by an applied span-varying force vector [(s) with components in the global reference system is expressed with respect to the infinitesimal displacement of the reference line along which it is applied. This reference line in the global system is given by Recall from figure 1 and define a and c as the offset of from the beam reference line in the and directions, respectively. This reference line follows the deformation of the beam structure and, in the absence of cross-sectional warping, a and c remain independent of the beam strain. Thus, the virtual work due to the applied load [(s) takes the form where the variational change in reference line is and Note that [ is general, in that its magnitude and orientation may vary arbitrarily in time and/or depend upon the states of the system, e.g. as in the case of a time-varying follower force.

External moment. Applied moments are denoted by the distributed vector [(s) for which the work done is considered with respect to the rotational variation (s). This variation is written as and thus

Kinetic terms

Consider the reference line (s) intersecting the centre of mass of all infinitesimal cross sections along the length of the beam. The inertial terms are furnished by considering the work performed when imparting a change in the momentum of the beam over incremental displacements (s) and rotations (s) of this reference line. The virtual work is thus written Expressions for (s) and its variation δ(s) take the identical form to equations (2.11) and (2.13) with the substitution . m(s) gives the mass per unit length along (s) and (s) the sectional inertia matrix. Additional point masses may be incorporated via substitution of the appropriate Dirac delta representations onto the reference line (s) and mass distribution m(s). The remaining accelerations and are detailed in the appendices.

Attitude representation

To further develop the describing equations of motion, one must relate the variational quantities , and to a consistent set of kinematic parameters over which the problem will be solved. Together these variations depend upon the deformation of the flexible beam as well as the translational and rotational motion of the reference system A.

For the deformation of the flexible structure, one may consider the components of the vector  = (τ(s), ε(s), τ(s), κ(s), κ(s), κ(s))T in the material frame as a natural description of this deformation. However, although the above variations are uniquely determined by these components, there is no closed-form solution directly relating the curvature in the material frame to [ and [ for all but the simplest κ(s), κ(s), κ(s) distributions (note that Cesnik & Brown circumvent this issue by assuming simplified κ approximations within individual element frames, building up complexity via discretization [9]).

For the continuous formulation developed here, this lack of a closed-form relation complicates the treatment of these quantities and their derivatives. To address this issue, one notes from equations (2.1), (2.15), (2.3) and (2.4) that these variations (in the reference system A) are dependent on the intrinsic coordinate system e[ such that: where denotes the set of attitude parameters uniquely defining e[. There are numerous potential choices for such a parametrization including, for example, quaternions, rotation vector, Rodrigues parameters; these and other such parametrizations are detailed in [28]. For the remainder of this study, an example set of asymmetric 3-1-2 Euler angles is employed; like other three-parameter descriptions this has the advantage of a minimal state representation and in the context of the current formulation is found to yield fast computational performance. The attitude states—here denoted  = (θ, ψ, ϕ)—relating to this parametrization are depicted in figure 3. The indicated angular deflections are labelled in a consistent manner with the particular positive convention employed in equation (2.19). Given these angular definitions, and taking the reference system A, where the shorthand convention s and c for the trigonometric operations ‘sin()’ and ‘cos()’ is used; E denotes the optional transform specifying the mapping from the Euler axis (E, E, E) to the intrinsic system (, , ). Via the choice of E one may place the θ =  ± 90° Euler singularities of the problem at generic opposing orientations of the intrinsic reference frame A. For example, when E takes the value of the 3 × 3 identity matrix, the Euler singularities fall at intrinsic frame orientations satisfying the equality  = (0, ± 1, 0)T.

Figure 3.

The set of Euler angles used to parametrize the intrinsic system orientation.

Note that an additional attitude and position parametrization can be assigned to the motion of system A with respect to the global reference G for problems where the root kinematics are not prescribed. This is useful in treating free problems, e.g. the modelling of a flexible aircraft in free flight [9,14,24,29]. For brevity, these terms are not furnished in this derivation but are easily appended to the problem; therefore, from here and [](t) take the form of prescribed time-varying functions.

Given the aforementioned attitude parametrization, the following quantities may now be related directly to the Euler angles θ, ψ and ϕ.

Curvature: Spin: and In summary, the variational terms of the problem have now effectively been cast in the form parametrized by the components of the kinematic vector where

Generalized coordinates

In line with d'Alembert's principle, one may state that the dynamic evolution of the system will proceed such that the virtual work contributions of equations (2.9), (2.10), (2.12), (2.16) and (2.17) will sum to zero for all admissible variations; thus It has been shown that all of these virtual work contributions may be expressed in terms of the kinematic vector . In order to construct the temporal ordinary differential equations capturing the system dynamics, the components of are written in a typical s, t separable form as a summation of shape functions where is the set of shape functions approximating the deflected state of the flexible beam continuum and ζ(s) denotes some optional reference distribution for each kinematic parameter from which the shape functions deform the system. q are the time-dependent generalized coordinates that will form the state vector of this beam formulation. Note that the tensor notation is adopted throughout the remainder of this paper wherein one is to sum over all values of any index that appears only in multiplicative pairs. Hence, equation (2.26) is written equivalently as

Equations of motion

The partial derivatives of (2.25) are now taken with respect to each of the generalized coordinates q: These terms are broken down into their constituent parts via chain differentiation in appendix Aa.

Writing the above equation (2.28) as a matrix system yields where the differential matrix operator To facilitate solution of this system, the acceleration terms in the above formulation are rearranged, factoring out any expressions that include the second-order state derivative . This rearrangement concerns the terms and and is detailed in appendix Af. Consequently, the kinetic virtual work ∂ may be split into the following parts (example matrix dimensions are also provided where Q is equal to the number of states in ). where Here, is the angular velocity of the coordinate system A.

Given this rearrangement of the dynamic terms, the final equations of motion for the nonlinear beam take the form of the ODE system: Here and are given by equation (2.31); the remaining terms , , and are as defined by equations (A 1), (A 2), (A 3) and (A 4), respectively, with constituent relations feeding into these terms as detailed in appendix A(b–e).

Shape function basis

In order to implement the formulation detailed in the preceding section, one must choose a suitable set (or sets) of shape functions from which approximations of the components of can be constructed. The choice of such a set is quite general, requiring only satisfaction of the essential kinematic boundary conditions of the problem. In the context of this formulation, these take the form of prescribed boundary values applied to a subset of components at one or more spanwise locations along the beam. Considering equation (2.27), one may write this condition where the value is constant. This condition holds for any constrained component of . It is noted here that the reference system A, in which this condition is cast, is a non-inertial system and thus (3.1) still admits arbitrary translation and rotation of the beam in the global coordinate system G (see, for example, the spinning beam test case of §4c).

A suitable set of shape functions satisfying (3.1) is depicted in figure 4a; these functions are based on shifted and reversed Chebyshev polynomials of the first kind. In panel (b), the same set is depicted following multiplication by a scaling function E(s) to bring the left side to zero (similar such scaling may be applied to the right-hand side for tip boundary conditions). This scaling procedure is applied throughout the examples of this study to constrain the applicable components of . A recursive generating relation for this set of functions y(x) is given by Note that the use of the scaling function E(s) retains orthogonality of the underlying Chebyshev set; however, this is not a required property, nor must each component of draw from the same shape function basis. Furthermore, because of the use of this scaling function, it becomes trivial to satisfy the boundary condition (3.1) for other arbitrary shape sets. However, the relative merits of different candidate shape bases are not pursued further in this study; rather the intention is to show that for the general set (3.2) satisfying the minimal essential boundary condition of the problem (3.1), the various test cases detailed in §4 can be efficiently treated (note also that throughout the cases of this paper, ζ(s) = 0). For similar reasons, the properties of orthogonality and natural load-dependent constraints are not discussed further here.

Figure 4.

The first eight functions of the candidate shape set based on Chebyshev polynomials of the first kind. Panel (a) depicts the pre-multiplied set; (b) depicts the scaled set following application of the envelope function E(s).

Before proceeding, it is noted that no condition on the smoothness or continuity of these shape functions has thus far been imposed; indeed any discontinuities of these sets may be treated via appropriate placement of the integration bounds of the weak formulation (2.35). Thus, discontinuous features of the physical problem may be treated by embedding the same discontinuities in the approximating shape basis. To illustrate this, consider the modelling of a straight horizontal cantilever beam with the following characteristics: firstly, the inner half running between the root and midpoint is assigned a much greater stiffness than the outer half from midpoint to tip; secondly, a point mass is fixed to the midpoint of the beam at the boundary of these two sections. One (basic) way of embedding these discontinuities into the underlying shape function set is to furnish the existing smooth functions with two additional non-smooth functions with appropriate mid-span C1 and C2 discontinuities; an example Chebyshev derived set, extended in this fashion, is depicted in figure 5a with the additional C1 and C2 discontinuous functions given by the dashed line and thicker line labelled 1 and 2, respectively (the remaining Chebyshev functions are indicated in increasing order by the grey curves labelled 3–9). In panel (b), the treatment of an example static test is shown; here a spring is connected between the beam tip and point ‘P’. The combined effects of self-weight and tip force induced by the spring cause the beam to assume the indicated shape with a clearly visible larger deformation of the more flexible outer section. The contributions of each of the shape functions to this static solution are given by the inset bar diagram in panel (b); the index of each bar refers to the corresponding shape function labelled in panel (a). Note that a more sophisticated tailoring of the underlying shape set may be employed in the efficient treatment these or other such localized characteristics; however, such a discussion is not pursued further in this study.

Figure 5.

An example non-smooth beam problem. Panel (a) gives the modified shape function set used to embed the stiffness and mass discontinuities in the numerical problem and (b) depicts the predicted static deflection obtained as a consequence of the mass distribution and tip spring.

It is pointed out here that the treatment of such discrete features also extends to geometric discontinuities in the structure such as for a kinked beam or framed structure (see examples in [3,10,30-32]). In such cases, the discontinuity must also be represented in the kinematic description of the system. There are multiple entry points in the above formulation where this discontinuous information may be embedded. For example, one possibility is to define a generic mapping from the underlying Euler-referenced triad system e[ to the structure-aligned intrinsic triad system e[ used in calculating the global reference line [, i.e. Another possibility is to directly incorporate appropriate Dirac-delta functions into the initial curvature distribution 0 to build up the discontinuous global geometry. These potential extensions of the formulation and their relative merits in the treatment of frame structures and other discontinuous geometries lies outside of the scope of this study.

Finally, it is worth mentioning that, alternatively, for the above examples one may of course treat smooth sub-domains of the problem breaking the discontinuities across designated element boundaries. Combined with a suitable choice of finite-element type basis set one approaches the more typical treatment of other element-based geometrically exact methods. However, in doing so one invariably moves away from the minimal state representation sought in this study; furthermore one notes that although smooth problems are exclusively treated in the examples of this study, the arbitrary smooth variation of physical parameters as displayed in the example of §4e is successfully treated in a geometrically exact manner without subdivision of the problem; this feature is novel to this nonlinear beam shapes formulation.

A note on numerical integration

Consider the numerical integration of the nonlinear beam shapes formulation. Recall the previously derived equation of motion Since the choice of shape functions is arbitrary, the constituent integral terms of and cannot typically be treated analytically. Specifically, the terms , , ∂, , , ∂, ∂, ∂ and ∂ are integrated numerically over s∈[O, L] when evaluating the equation of motion at a particular time step. Here, numerical evaluation of these terms is accomplished using quadratic interpolation. Noting that the full system is stiff, one may compute a dynamic trajectory in time using any stiff first- or second-order solver capable of treating systems of the form (3.3). Since this method is implemented in Matlab, the in-built variable-step, variable-order stiff ODE solver ‘ODE15s’ [33] is applied to the first-order form of (3.3).

Numerical test cases

A variety of test cases, comprising benchmark tests from the literature as well as some further examples, is now presented. Each case targets particular characteristics of the geometrically nonlinear flexible beam problem. In addition to published results, the responses generated by MSC Nastran and the intrinsic beam formulation of Hodges [13,23] are in places used to validate the nonlinear beam shapes approach proposed in this paper. For all the cases, the Chebyshev shapes detailed in §3 are used in approximating τ(s), ε(s), τ(s), θ(s), ψ(s) and ϕ(s).

Large static deformation

The first test case considered is well documented in the literature and treats the large deformation of a pre-curved beam. The initial shape of the beam forms one-eighth of a circle of radius 100 m in the horizontal (x,y) plane (see the thin curves in figure 6). The beam has a 1 m2 cross section, a Young's modulus of E = 107 N m−2 and a Poisson's ratio of zero (thus EI = EI = 25/3 Nmˆ2, GJ = 7.02885 × 105 Nm rad−1); gravitational acceleration is not considered. The system is solved statically using the Euler mapping moving the singular points away from vertical (all other test cases in this study are treated with E equal to the 3 × 3 identity matrix). Two load cases are shown in figure 6, each involving the application of a 600 N tip load. In panel (a), the tip load is applied in the vertical z-direction; in panel (b) the same load is applied as a follower force (i.e. remaining parallel to (L)). In table 1, the static tip deflection in the global coordinate system is stated as predicted by equivalent tests from a number of published sources. Results using Nastran and the intrinsic beam formulation are also displayed and a good agreement is observed.

Figure 6.

Beam deformation given (a) a 600 N vertical tip load and (b) a 600 N follower tip load.

Table 1.

Comparison of predicted tip positions for the curved beam static test case.

source	tip position (m) Fvertical	tip position (m) Ffollower
Bathe & Bolourchi [34]	(15.9, 47.2, 53.4)	—
Simo & Vu-Quoc [3]	(15.79, 47.23, 53.37)	—
Géradin & Cardona [35]	(15.55, 47.04, 53.50)	—
Li & Vu-Quoc [36]	(15.54, 46.85, 53.64)	—
Crisfield [30]	(15.61, 46.84, 53.71)	—
MSC Nastran	(15.56, 46.89, 53.60)	(−10.92, 24.53, 59.42)
intrinsic beam	(15.56, 46.90, 53.60)	(−10.93, 24.55, 59.41)
nonlinear beam shapes	(15.55, 46.90, 53.60)	(−10.95, 24.54, 59.41)

Problem size

To assess the number of states required for the various methods, consider the following convergence criteria on the required number of states n*: i.e. doubling the number of states alters the predicted tip deflection by less than 1 mm. For the nonlinear beam shapes method, this doubling is applied to the discretization for each component; for the special case of zero component shapes, one shape is added rather than doubling the set size.

Using equation (4.2), the total number of states required by the nonlinear beam shapes method, the intrinsic beam formulation and MSC Nastran is shown in table 2; results are given for both the tip vertical (figure 6a) and tip follower (figure 6b) load cases. For the current formulation, the number of shape functions assigned to τ(s), ϵ(s), τ(s), θ(s), ψ(s) and ϕ(s) is indicated in brackets, respectively. In addition, table 2 also provides a rough comparison of the required computational effort in treating these systems. To achieve this, the follower and vertical load cases are run as dynamic simulations over a period of 60 s. A large stiffness proportional damping is used ( =  − , ) such that the responses are settled to their respective equilibria by the end of the 60 s simulation. The nonlinear beam shapes and intrinsic beam methods are run in Matlab using the stiff ODE solver ‘ODE15s’ [33]; no Jacobian information is fed to the solver for the former method. Nastran is run using its in-built SOL400 solution routine with adaptive time stepping. Owing to the large deformations involved in this test case, it proved difficult to obtain an accurate Nastran transient response for the follower 600 N tip load. Hence only the Nastran benchmark time based upon the vertical load case is provided in table 2. Because of the differences in the numerical integration schemes only a crude comparison of computational speed can be inferred here. Looking to the results of table 2, the nonlinear beam shapes formulation required an order of magnitude fewer states, accompanied by a significantly reduced simulation time when compared with the other methods. Of particular note is that no shear or extensional states were required by the nonlinear beam shapes formulation in meeting the convergence criterion in this example. Indeed, this is the case for a majority of the example problems of this study wherein capturing large nonlinear deformations results in the treatment of flexible and slender structures where shear effects are not significant. One advantage of this formulation is that the shear and extensional compliances can easily be removed by simply not assigning shapes to τ(s), ϵ(s) and τ(s). Note however that the shallow arch example of the following section introduces a test case for which these additional shear/extensional flexibilities do prove influential on the resulting structural response.

Table 2.

Static degrees of freedom required in satisfying the convergence criterion (4.2), and computational benchmarks.

method and load case	problem discretization	static doFs	dynamic solution time (s)
Nastran (vertical)	22 elements	132	40.1
Nastran (follower)	60 elements	360	—
intrinsic beam (vertical)	31 elements	186	132
intrinsic beam (follower)	69 elements	414	760
nonlinear beam shapes (vertical)	11 shapes (0, 0, 0, 4, 4, 3)	11	3.26
nonlinear beam shapes (follower)	15 shapes (0, 0, 0, 5, 5, 5)	15	4.25

Shallow arch example

This test case focuses on the response of a pin-jointed shallow arch subject to a distributed radial load. The motivating features of this test case are two-fold. Firstly, the treatment of non-cantilever boundary conditions provides a further example of the construction of equivalent kinematic boundary conditions in the context of this formulation. Secondly, the loaded arch provides a system for which the shear/extension states prove influential on the observed response.

The shallow arch is based upon the pre-curved beam of the previous example. The cross-sectional properties remain the same; however, the radius of curvature is now reduced to 20 m—consequently reducing the beam length over the 45° arc. This curved beam is supported on two pin joints to form an un-stressed arch to which a radial load is applied. This problem is planar and oriented in the (y,z)–plane; thus ψ(s) = ϕ(s) = τ(s) = 0, a condition achieved simply by assigning no shape functions to the 1st, 5th and 6th components of . The remaining θ(s), ε(s) and τ(s) kinematic parameters are unconstrained at both boundaries of the arch, thus the unscaled shape set of figure 4a is selected for these components. [(0) = 0 and [(L) =  correspond to the displacement constraints at each end of the arch. [(0) = 0 is implicitly satisfied by equation (2.1) where ; to constrain the right-hand boundary equation (2.28) is furnished with the additional energy term where the multiplier Υ is given a sufficiently large value to constrain [(L).

Application of a distributed radial load will cause this arch to buckle. Despite the simplicity of the configuration, there are a number of distinct buckling behaviours that may be observed for this system. In [37], Pi et al. classify these as snap-through and bifurcation buckling modes, the latter of which is typically asymmetric in nature. The specific buckling characteristic observed is related to the non-dimensional ‘shallowness’ parameter λ. For a rigid-pinned circular arch such as that treated in this example this parameter takes the form where A is the cross-sectional area of the beam, I1 is the second moment of area and R is the radius of curvature. This yields a lambda value of 10.68 for this test case.

Looking now to figure 7, the shape of the shallow arch is shown following application of a radial load, and snapshots of the deformation of the arch as it buckles are depicted. In panel (a), a distributed load of 7500 N m−1 is applied; no shear or extensional compliance is permitted and the arch is observed to buckle under an asymmetric bifurcation. This agrees with analytic results for a rigidly supported arch with λ above the critical value of approximately 9.2. The effect of elastically compliant supports are also considered in [37]. It is shown that as the flexibility of these supports is increased, the effective λ is reduced, allowing for symmetric snap-through behaviour. Rather than directly modifying the support flexibility in this example, the addition of shear and extensional states is used to produce an alternative source of radial and axial compliance to the system. In figure 7b, five shapes are assigned to each of the extensional ε(s) and shear τ(s) parameters; the applied load is reduced to 5500 N m−1 in light of the increased flexibility of the arch. Looking at the snapshots of the arch deformation, one observes that the addition of shear and extensional states is indeed sufficient to capture the existence of a symmetric snap-though mode.

Figure 7.

Arch buckling following application of a distributed load. Panel (a) excludes shear and extensional states; (b) includes shear/extension. In both cases, five shapes are used for each active component of (s).

Rotating pre-curved beam

The following two examples are based upon the test beam treated in Pai [31]. The parameters are summarized in table 3. For the first test case, the beam is oriented in the global y-direction ((0) = (0, 1, 0)T) and assigned a constant pre curvature κ =  − (3π)/(18L) m−1.

Table 3.

Pai beam parameters from [31].

parameter	value	units
length	0.479	m
width	5.08 × 10−2	m
height	4.5 × 10−4	m
EIxx	4.899184 × 10−2	Nm2
EIzz	6.243471	Nm2
GJ	7.164459 × 10−2	Nm rad−1
mass per length	1.012698 × 10−1	kg m−1
torsional inertia per length	2.178012 × 10−1	kg m

This pre-curved geometry is shown in figure 8a by the curve labelled ‘undeformed’. The additional results are obtained by first allowing the beam to deflect under its own weight (0 Hz curve). Then the root is spun in the global z-direction at a number of angular velocities. The settled deformation as predicted by the shape function formulation is depicted by the labelled red lines for rotation speeds of 1, 2, 3, 4, 6 and 8 Hz. Corresponding data points from [31] are overlaid onto these curves illustrating a very close agreement. For the second phase of this test, the root of the beam is fed a 4.5 Hz harmonic excitation in the vertical direction of amplitude ± 2 cm in addition to the applied angular rotation. To help with visualization, the smaller plots at the bottom of figure 8 (panel (f)) illustrate the oscillating and spinning motion of this system in the global coordinate system (X[, Y[, Z[). A stiffness proportional, viscous damping law is assumed at 1% of the structural stiffness; i.e. in equation (2.10). The blue solid line in panel (b) shows the tip response for 10 s of this test case. At first, the rotational velocity of the system is zero and thus the initial beam response is due purely to the root vertical oscillation. The settled response to this excitation is shown for the first 2 s of the time series. Then for 2 ≤ t ≤ 10 the rotational velocity |[| is ramped linearly from 0 to 8 Hz. Over this time, the beam tip continues to oscillate and remains centred about the dashed red line which corresponds to the static solution without the root vertical oscillation (compare with panel (a)). At the onset of the rotational sweep (t = 2) one also notes both an out-of-plane and torsional component of the beam response (depicted in panels (d) and (e) by the non-zero ψ(L) and ϕ(L) components of the response). This illustrates the expected gyroscopic effects resulting from large local rotational motion of the structure in the rotating beam reference frame.

Figure 8.

Predicted responses for the pre-curved spinning beam test case.

Looking at these responses, one notes that the amplitude of oscillation grows as the root angular velocity is increased reaching a maximum before dropping off again. Physically, one is observing here the increasing centrifugal dynamic component having a stiffening effect on the flexible system and consequently increasing the frequencies at which the flexible modes reside (see, for example, the modelling studies of [38-40] for treatments of the centrifugal stiffening effect in the context of simple rotating systems). This stiffening is illustrated in panel (c) where the relationship between the root angular velocity and the first modal linear frequency is depicted. For a rotational velocity of 3.91 Hz—occurring at t = 5.91 s in the simulated test case and indicated by the vertical line in panel (b)—this modal frequency is equal to 4.5 Hz and therefore equal to the frequency of root excitation of the system. Given a small amplitude excitation and sufficiently slow sweep, the maximum response amplitude would be observed close to this angular velocity. However, for the level of harmonic excitation applied in this example the amplitude of the system response is sufficient to push into the nonlinear regime of the system dynamics and consequently the resonant frequency at the observed amplitudes will not coincide with this linear frequency; coupled with the relatively quick traversal through this resonant region one observes a delayed maximal response at a larger rotational frequency. To verify this case, the same test was performed using the intrinsic beam model; the resulting tip deflection is also depicted in the background of panel (b) by the thicker background line. Here, one sees that the obtained envelope of oscillation is almost identical between the two methods and that the time series match well, deviating only very slightly towards the end of the simulation.

Harmonically excited vertical cantilever

This next test case also concerns the beam detailed in the previous example. A harmonic root excitation is again fed to this system, this time in isolation. The beam configuration matches the corresponding experimental/numerical study conducted in [31]. Specifically, the beam is oriented pointing vertically upwards and subject to its own weight. The transverse harmonic oscillation applied to the beam root acts in the ± (0) direction and has an amplitude of ±1 mm. Structural damping is again stiffness proportional and set to a very light value of 0.005% of the beam stiffness. Two slow frequency sweeps of the root harmonic excitation are performed over the range 7–11 Hz in both the forward and reverse directions. The amplitude of tip oscillation relative to the base position, as this excitation frequency is varied, is plotted by the thicker curves in figure 9. These curves form approximations of the upper and lower stable periodic branches of the frequency response function for this level of root excitation and damping (the curves are clipped off at points where a loss of stability is observed). The colour at any location along these curves corresponds to the observed phase shift between the tip response and base excitation as indicated by the bar to the right of the figure.

Figure 9.

Example frequency response diagram approximated via numerical simulation of the system for a 0.005% stiffness proportional damping and ±1 mm base excitation. Estimates of the underlying nonlinear normal mode, as predicted by the current formulation and intrinsic beam model, are indicated by the thin red and thin grey curves, respectively. Snapshots of the beam oscillations are provided at points A (8 Hz), B (9.64 Hz) and C (10 Hz) along the upper periodic branch.

The phase of the periodic branch formed by the two stable parts and unstable part (not indicated) traverses a large portion of the 0–180° phase envelope as the excitation frequency is increased over the range illustrated. Three sub-plots are also provided, illustrating snapshots of the beam deflection at three sample locations along the upper branch labelled A–C. In the 8 Hz and 10 Hz sub-plots, the oscillatory response of the flexible beam is observed to correlate closely with the second linear mode shape of the system with a clearly defined nodal point observable. At the 9.64 Hz sample location, this vibration mode and characteristic nodal point is not cleanly observed. Looking to the upper branch topology at this sample location, one notes a localized ‘kink’ in the periodic curve and accompanying variation in phase lag. This feature corresponds to a 3 : 1 resonance between the second and third nonlinear flexible modes (for reference, the first three linear modal frequencies for the vertically oriented beam are 1.44, 10.42 and 29.55 Hz). Note that the precise topology of the periodic structure around this resonance point is not known, as the illustrated curve simply reflects the transient response as the base excitation sweep traverses this region; however, clearly the current formulation is capable of capturing these higher-order resonances. For further discussion of such resonant branch phenomena see [41] wherein the out-of-plane resonance of a similar square cross-sectioned flexible beam is detailed.

Following the upper periodic curve towards the left, one notes the phase lag approaches (and passes through) 90°. Nonlinear vibration theory dictates that (providing the strain rate damping law is linear) this point of quadrature at 90° coincides with a point on the nonlinear normal mode (or backbone curve) of the flexible system. By allowing the unforced system to decay from this point of quadrature, one may trace a portion of this backbone curve by sampling the transient frequency and amplitude of the decaying response. The resulting backbone section is shown in the figure by the thinner red curve, beginning at the quadrature point and decaying in amplitude until the linear frequency at 10.42 Hz is met (note that the exact backbone curve for this system extends indefinitely beyond this quadrature point). The overall curve has a clear left tilt in the figure and thus indicates a softening of the second nonlinear normal mode as the amplitude of oscillation is increased (also observed in [31]). One notes that the resonant ‘kink’ previously observed in the upper branch also helps shape the transient decay and subsequent backbone prediction in the region 9.5–10 Hz.

As in the previous examples validation is sought via comparison with the intrinsic beam formulation; performing an identical decay from the same initial state yields a second backbone approximation predicted by the latter method and indicated by the grey curve. One notes close agreement between each backbone result, both tracing a similar shape and indicating the nonlinear softening effect. Note that there is some slight difference between the two results in the 9.5–10 Hz range; however, due to the speed of decay it is unlikely that either method closely captures the localized backbone structure in this region. For the precise tracking of this backbone curve, one could employ more sophisticated methods based upon shooting or pseudo-arclength continuation [42-44], however, this is not pursued here.

Wind turbine blade

The final test case details the static deformation of a wind turbine blade based upon the ‘NREL offshore 5-MW baseline wind turbine’ detailed in the technical report of Jonkman et al. [45] with supplementary parameters taken from the pre-design study of Kooijman et al. [46] (part of the Dutch Offshore Wind Energy Converter (DOWEC) project). The blade has a span of 61.5 m from its hub attachment point to the tip; the unloaded structural geometry exhibits variable material properties, pre curvature and twist from its root to tip. For the specific turbine blade of this test, arbitrary parametrized functions are fitted to the data and used to concisely express these span-dependent quantities; these functions are detailed in table 4. A depiction of the undeformed blade geometry is provided at the top of figure 10. For this test, the turbine blade is oriented such that E is equal to the identity matrix. A distributed load is applied to this structure along the reference line (s). The load has an elliptical profile and acts normal to the blade chord line, i.e. Two load cases corresponding to values of and in expression (4.6) are applied to the blade in panels (a) and (b) of figure 10, respectively. Each panel shows the static deformation achieved via dynamic simulation of the system, allowing sufficient time for transient motion to settle; a conservative choice of eight Chebyshev polynomials per Euler angle was used for the treatment of this test. For the first load case of panel (a), the applied load is just sufficient to remove the pre-curvature from the turbine blade. For the second load case of panel (b), the load produces an extreme deflection of the structure far in excess of any physically realizable load; this latter case is used to demonstrate the ability to solve the geometrically exact problem for large deformations of the generic non-prismatic structure. To verify these deformations both the net applied moment (a consequence of the applied load distribution ) and the internal moment (calculated from the Euler-Bernoulli constituent relation) are compared at each spanwise location in the range s∈[0, L]. In panel (c), the out-of-plane () component of the external (solid) and internal (dashed) moment distributions are plotted for the first load case; panels (d) and (e) depict the twist () and in-plane () components, respectively. Panels (f–h) provide the same three plots for the second, larger load case. Correct solution of the static problem may be verified by noting the close agreement indicated in these sub plots.

Table 4.

Wind turbine blade parameters.

span-dependent function	coefficients
length: L = 61.5 (m)
EI_xx=a_1+a_2x^1/a_3+a_4x+a_5x^2+a_6x^3+a_7x^4(Nm^2)x=1−s/L	a_1=12.04a_2=7.4296a_3=5.9294a_4=0.6791a_5=6.3115a_6=−2.6386a_7=−0.9242
GJ=a_1+a_2x^1/a_3+a_4x+a_5x^2+a_6x^3+a_7x^4(Nm rad)^−1x=1−s/L	a_1=12.15a_2=6.8864a_3=3.5466a_4=−1.2809a_5=2.0892a_6=2.8858a_7=0.0322
EI_zz=a_1+a_2x^1/a_3+a_4x+a_5x^2+a_6x^3+a_7x^4(Nm^2)x=1−s/L	a_1=15.43a_2=7.2698a_3=4.0351a_4=0.1553a_5=−0.2459a_6=0.3072a_7=0.4667
m=a_1x^1/2(1+a_2x^2+a_3x^3+a_4x^4+a_5x^5+a_6x^6)(kg m)^−1x=1−s/(L+0.2)	a_1=250a_2=−1.4119a_3=17.8617a_4=−21.5880a_5=−9.7330a_6=17.0449
θ_0=−tan^−1⁡((∑n=06a_nx^n)(tan^−1⁡(s−14)/π+0.5))(rad)x=(s−30.75)/18.02	a_0=0.029954a_1=0.039748a_2=−0.0031189a_3=−0.0074728a_4=0.0080149a_5=0.00047782a_6=−0.0021779
ψ0 = 0 (rad)	—
ϕ_0=(13.5−∑n=09a_nx^n) (tan^−1⁡(s−10)/π+0.5)−13.5(rad)x=(s−30.75)/18.02	a_0=6.5963a_1=−5.3307a_2=−0.71825a_3=0.73322a_4=2.6739a_5=−1.2031a_6=−1.7149a_7=0.80577a_8=0.30529a_9=−0.14539

Figure 10.

Example deflections for the wind turbine test case.

Conclusion

This paper has detailed the development and application of a nonlinear beam shapes approach capable of the low-order geometrically exact representation of a flexible beam; the formulation is based upon the shape-based discretization of attitude and incremental shear kinematic quantities distributed along the flexible structure. The development of the describing equations of motion was demonstrated for an example Euler angle parametrization and coupled to a set of Chebyshev polynomials of the first kind, modified to admit the generalized kinematic condition. The formulation was subsequently applied to a variety of test cases including the modelling of planar and non-planar transient dynamics, the prediction of static equilibria and the treatment of simple beam, pre-curved and non-prismatic structures, all within a geometrically nonlinear context. The shape functions used are generic and may be tailored to the characteristics of the modelled physical system. Complementary Nastran and intrinsic beam calculations were used to verify specific results throughout these tests; altogether the formulation of this study demonstrated the accurate treatment of the considered test cases, consistently achieved using a small set of states, representing an order of magnitude reduction in problem size over its element-based counterparts. This low-order formulation admits the efficient treatment of flexible beam analyses for which the size of system generated by traditional element-based representations would incur a significant computational penalty; such examples include structural optimization problems over large parametric design spaces, sensitivity and uncertainty quantification analyses, and numerical continuation methods exploring limiting dynamic behaviours.