A geometrically-exact momentum-based non-linear theory applicable to beams in non-inertial frames

International Journal of Non-Linear Mechanics 113 (2019) 158–170

Contents lists available at ScienceDirect

International Journal of Non-Linear Mechanics journal homepage: www.elsevier.com/locate/nlm

A geometrically-exact momentum-based non-linear theory applicable to beams in non-inertial frames Shanran Tang, Bert Sweetman ∗ Department of Ocean Engineering, Texas A&M University, College Station, TX 77843, United States

ARTICLE

INFO

Keywords: Geometrically-exact theory Non-linear beam dynamics Separation of displacements Finite-volume beam analysis

ABSTRACT A geometrically-exact non-linear beam model is developed based on conservation of momentum for application to arbitrarily-shaped beams having large deformations and large overall motions. Coordinate transformations are used to derive the non-linear inertial forces and moments and the non-linear relationships between displacements and strains, enabling rigorous consideration of kinematic and geometric nonlinearities. General non-linear equations of motion are first derived in a differential form, and then are simplified using practical engineering assumptions, including developments of a linearized model for small angles and small strains and a discretized model using finite volumes. A separation of displacements technique is proposed, which enables non-linear beam dynamics to be rigorously reduced to a series of piecewise-linear models. The proposed model is unique amongst existing geometrically-exact beam formulations in that it is formulated using dynamic quantities relative to intermediate non-inertial coordinates, allowing general Lagrangian formulations to be established in floating frames with inertial coupling effects, which is compatible with multibody models. The practical value of these theoretical developments is demonstrated through numerical implementation and convergence analyses. The piecewise-linear dynamic solver in which displacements are computed from an evolving configuration is shown to capture non-linear dynamic behaviors using linear solutions without iteration at each time-step.

1. Introduction Beam theories are used to reduce three-dimensional (3D) slender structures to one-dimensional (1D) beam models for analysis. The Euler–Bernoulli beam theory was developed in the mid 18th century and continues to be widely used today. Timoshenko improves the Euler–Bernoulli beam theory by including the effects of shear deformation and rotational inertia. These classical beam theories are inadequate to properly address non-linear behaviors of some beam structures. Several types of nonlinearities are inherently associated with beam dynamics: (a) geometric nonlinearities are non-linear relationships between displacement and strain, which are introduced by large deformations and induce geometric stiffening effects [1,2]; (b) constitutive nonlinearities are non-linear relationships between stress and strain resulting from material properties, or non-linear relationships between stress and resultant forces and moments resulting from geometric deformations; (c) kinematic nonlinearities are non-linear relationships among displacements, velocities, and accelerations, typically associated with large angular motions; (d) inertial nonlinearities occur when the mass moment of inertia and mass per unit length depend on geometric deformations; (e) various non-linear damping effects also exist. Nonlinear strain measures are used with classical beam theories in order to

address geometric nonlinearities, yet this approach relies on the same fundamental assumptions made in the Euler or Timoshenko theories. Modern beam theories are developed to address beam nonlinearities. Geometrically-exact beam theory (GEBT) is a class of non-linear beam theories [3,4], which allows for non-linear displacement–strain relationships and non-linear beam kinematics. Reissner [5] derives the non-linear strain measures using the principle of virtual work for static two-dimensional (2D) beams. Reissner [6], Reissner [7] extends the derivation from a 2D plane to a 3D space to develop a largedisplacement static theory for space-curved beams. Simo [8],Simo and Vu-Quoc [9] derive non-linear kinematics and strain measures for beams in a 3D space and propose a geometrically-exact dynamic beam theory. Simo and Vu-Quoc [10] further develop the theory to include shear and torsion-warping deformations. Bauchau and Hong [11,12] develop a small-strain non-linear beam model in which cross-sections do not deform in their planes while St. Venant warping functions are used for out-of-plane deformations. Danielson and Hodges [13] derive a non-linear displacement–strain relationship for curved and twisted beams with large displacements and rotations, and significantly simplify the strain measures based on small-strain and small-localrotation assumptions. Danielson and Hodges [14] further develop this

∗ Corresponding author. E-mail address: [email protected] (B. Sweetman).

https://doi.org/10.1016/j.ijnonlinmec.2019.03.007 Received 24 August 2018; Received in revised form 25 February 2019; Accepted 18 March 2019 Available online 26 March 2019 0020-7462/© 2019 Published by Elsevier Ltd.

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

frame method, (b) the incremental FE approach, (c) the finite segment method, (d) the large rotation vector method, and (e) the absolute nodal coordinate approach. The floating frame method includes two sets of coordinates: one set is the intermediate floating frame that describes the rigid-body motion of a beam; the other set describes the local deformations in the floating frame. Local deformations are usually solved using classical beam theories because the formulation of modern non-linear beam theories in an intermediate non-inertial frame is not commonly available [36,37]. The incremental FE approach is based on convected coordinate systems in which individual beam elements are considered as deformable multibodies. The finite segment method differs from the incremental FE approach in that beam segments are assumed as rigid multibodies connected by springs and/or dampers. Both the incremental FE method and the finite segment method make use of multibody dynamic theories to compute the body motion of each beam element [51,52]. The large rotation vector method is equivalent to application of the FE method to non-linear beam theories in inertial frames, which requires special interpolation of finite rotations [42,43, 53]. The absolute nodal coordinate approach differs from the large rotation vector method in that the rotational variables are eliminated from the formulation: higher-order global shape functions are used to interpolate absolute translational displacements only, such that rotations are included implicitly [54–56]. The new theory presented here was originally motivated by dynamic response of blades on highly-compliant floating offshore wind turbines, in which the blades undergo not only planar rotation but also irregular large overall motions [57]. The objective of this work is to develop a new geometrically-exact model for arbitrarily-shaped beams having large deformations and large overall motions in non-inertial floating frames, and to demonstrate the new theoretical developments in an effective numerical implementation. The proposed momentumbased beam theory (MBBT) can be used to rigorously address geometric and kinematic nonlinearities and the inertial coupling introduced by floating frames. The new theory is different from other existing GEBT formulations. First, it is formulated in a floating frame with dynamic quantities relative to non-inertial frames, which enables the total Lagrangian formulation [58] to be established including inertial coupling effects. The resulting formulation explicitly includes overall rigid-body motion and no special treatment is required to remove or restore overall rotation during numerical implementation, which makes it readily compatible with multibody dynamic models. Second, the derivation is based on conservation of momentum in vector form and use of coordinate transformation matrices. The resulting equations of motion (EOM’s) retain intuitive physical significance, which allows rational linearization using practical engineering assumptions. The value of the coordinate transformation and of rational linearization is demonstrated by the subsequent development of a separation of displacements technique and derivation of a linearized updated Lagrangian formulation. The new linearized formulation can be solved using conventional discretization methods without special treatment, and has been shown to effectively represent non-linear beam dynamics using a series of piecewise-linear solutions. This paper includes the complete derivation of MBBT and its further theoretical developments, as well as numerical implementation of these developments. MBBT is derived in Section 2, and is linearized, simplified, and discretized in Section 3. Separation of displacements and the resulting piecewise-linear MBBT are presented in Section 4, including linearization of the EOM’s about a deformed configuration. A static and a dynamic MBBT numerical solver are each developed in Section 5, where convergence analyses and several numerical examples are presented.

work into a static non-linear beam theory for moderate local rotation. Hodges [15] derives a mixed-form intrinsic formulation of GEBT by extending the static theory to dynamic cases using a variational calculus approach to apply Hamilton’s principle. Other GEBT formulations also exist: Bauchau and Hong [12] establish a displacement-based formulation; Zupan and Saje [16] and Su and Cesnik [17] establish a strain-based formulation. Constitutive relations play an important role in non-linear beam theories. Deriving constitutive equations from strain energy functions is a widely used approach in the early works [8,10,15]. Using a separate cross-sectional analysis to reduce dimensions (3D body to 1D beam) and to derive constitutive equations are more popular in recent studies, especially for composite beams. Variational asymptotic beam sectional analysis (VABS) is proposed in which the variational asymptotic method is used to evaluate the constitutive properties of cross-sections [18–22]. Cross-sections are modeled using 2D finite elements (FE’s) in VABS such that a 3D geometrically non-linear analysis for slender structures is split into a 2D cross-sectional analysis and a 1D non-linear beam analysis. Recent studies on composite wind turbine blades [23,24] show that VABS is sufficient to provide crosssectional constitutive matrices including effects such as cross-section warping, irregular shape, and varying properties of composite materials. Additional assumptions may be applied in cross-sectional analyses depending on how the constitutive equations are derived: different material models can be assumed with different strain energy functions, and certain cross-sectional strains or warping measures often are assumed to be small. Flexible beam structures having large overall motions are found in many engineering applications, such as wind turbine blades, satellite antennas, and robotic arms. Many recent studies on beam dynamics address coupling effects between overall beam motion and local beam deformation. Vibration of simply rotating beams has been extensively investigated, including application of classical beam theories and non-linear strains in many studies [25–31]. Simo and Vu-Quoc [32,33] apply GEBT to planar beams having large deformations and simple rotations. Kim et al. [34] and Zhao and Wu [35] apply similar geometrically-exact theories to 3D beams having simple rotations. Dynamics of 3D beams undergoing general overall motions also has been studied, including application of classical beam theories [36–39]; these formulations generally neglect some of the effects of structure coupling and large deformations. Application of GEBT to 3D beams having large deformations and general overall motions is found in a few studies, e.g., [15,40,41]. Hodges [41] presents an inertial frame formulation similar to the one proposed by Simo and Vu-Quoc [40]: dynamic quantities are relative to inertial frames in both GEBT formulations. Implementation of modern non-linear beam theories relies on numerical methods, most commonly the FE method. The conventional FE method interpolates a continuous quantity between nodes as a sum of shape functions, which is not directly applicable to 3D finite rotations because most rotation parameterizations are non-additive quantities. Interpolation of rotations requires special treatment to address the non-linear character of 3D finite rotations [42–45]. Strain-based and intrinsic formulations of GEBT allow application of FE interpolation to additive quantities [46,47], in which rotations are not interpolated as primary variables but are reconstructed through other quantities. Numerical implementation becomes more complicated for beams undergoing large general overall motions. The overall rigid-body rotation of an undeformed beam and the local rigid-body rotation of each deformed element are not directly additive, and both need to be considered carefully to ensure strain is invariant with rigid-body motion [48]. Removal of rigid-body rotation prior to numerical interpolation is proposed by Wang et al. [49]; the rigid-body rotation is then restored after interpolation to establish the governing equations in a inertial frame. Several numerical methods are developed for moving beams in the context of flexible multibody dynamics [50]: (a) the floating

2. Governing equations MBBT is derived using two fundamental assumptions: (a) crosssectional deformations are assumed to be small, such that cross-sectional 159

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

Translational and rotational displacements can be expressed as 3D first[ ]𝑇 [ ]𝑇 order tensors in vector form: 𝒖𝑙 = 𝑢1 , 𝑢2 , 𝑢3 ; 𝜽 = 𝜃1 , 𝜃2 , 𝜃3 . Bold font is used to indicate vectors or matrices throughout this paper. Coordinate transformation matrices are used to transfer vectors from one coordinate system to another. The subscripts of a transformation matrix indicate the post- and pre-transformation coordinates, and the subscript of a vector indicate in which coordinate the vector is resolved. For example, transformation matrix 𝑻 𝑑𝑙 enables a translational displacement vector resolved in local undeformed coordinate 𝑙 to be resolved in deformed coordinate 𝑑 as: 𝒖𝑑 = 𝑻 𝑑𝑙 𝒖𝑙 . Transformation matrices represent the change of directions between coordinates and are equivalent to rotation tensors that describe the rotation between basis vectors. The resulting transformation matrices are orthogonal: 𝑻 𝑙𝑑 = 𝑻 𝑇𝑑𝑙 = 𝑻 −1 𝑑𝑙 and 𝒖𝑙 = 𝑻 𝑙𝑑 𝒖𝑑 . No Einstein summation convention is applied in this paper. Transformation matrices are functions with respect to 𝑠 and 𝑡, and can be computed using rotational displacements. Angular velocities and accelerations are expressed using transformation matrices: ̃ 𝑑 = 𝑻 𝑑𝑙 𝑻̇ 𝑙𝑑 𝝎

Fig. 1. Coordinate systems.

(1)

̃ 𝑓 = 𝑻 𝑓 𝑔 𝑻̇ 𝑔𝑓 (2) 𝜴 𝜕 ̇ ̇ ̇ ̈ (𝑻 𝑻 ) = 𝑻 𝑑𝑙 𝑻 𝑙𝑑 + 𝑻 𝑑𝑙 𝑻 𝑙𝑑 (3) 𝝎 ̃̇ 𝑑 = 𝜕𝑡 𝑑𝑙 𝑙𝑑 𝜕 ̃̇ 𝑓 = 𝜴 (𝑻 𝑻̇ ) = 𝑻̇ 𝑓 𝑔 𝑻̇ 𝑔𝑓 + 𝑻 𝑓 𝑔 𝑻̈ 𝑔𝑓 (4) 𝜕𝑡 𝑓 𝑔 𝑔𝑓 in which vector 𝝎𝑑 is the angular velocity of local deformation at a cross-section measured in coordinate 𝑑; vector 𝜴𝑓 is the angular velocity of overall rigid-body motion measured in coordinate 𝑓 ; ( ̇ ) and (̈) denote the first and second time derivatives, and (̃) denotes the skew-symmetric matrix of a vector, which can be used to represent cross products as matrix multiplications. The absolute translational velocity of a differential cross-sectional element on a beam is measured in inertial coordinate 𝑔. Taking the 𝑊𝑃 𝑂𝑊 time derivative of position vector 𝑹𝑂𝐺 +𝑻 𝑔𝑓 𝑹𝑓 0 +𝑻 𝑔𝑓 𝑻 𝑓 𝑙 𝒖𝑙 + 𝑔 = 𝑹𝑔 𝑃𝐺 𝑻 𝑔𝑓 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑹𝑑 yields the absolute translational velocity:

mass properties, mass per unit length and mass moment of inertia, are constant in time, the result of which is that inertial nonlinearities are not considered; (b) the cross-sectional internal forces and moments are functions of strain measures and/or time derivatives of strain measures, such that the constitutive relations do not introduce new unknowns to the system of governing equations. The derivation includes rigorous consideration of geometric and kinematic nonlinearities under these assumptions. 2.1. Coordinate systems The coordinate systems used in MBBT are presented in Fig. 1. An inertial Cartesian coordinate system, 𝑿𝒀 𝒁, with origin 𝑂 fixed to the earth, is denoted as global coordinate 𝑔. A non-inertial Cartesian coordinate, 𝒙𝒚𝒛, is fixed to the undeformed beam at point 𝑊 , and is defined as floating coordinate 𝑓 . Translations and rotations of 𝑓 relative to 𝑔 represent the overall rigid-body motion of an undeformed beam. Curvilinear coordinate 𝑠 is attached to the elastic axis of an undeformed, unstressed beam; another curvilinear coordinate 𝑆 is similarly attached to the elastic axis of the deformed beam. Two sets of local coordinates are defined at each point along 𝑠 and 𝑆. The non-inertial local undeformed coordinate 𝑙 is fixed to elastic center 𝑃0 , and includes Cartesian coordinates 𝒂𝒃𝒄 with 𝒃 and 𝒄 on the cross-sectional plane and 𝒂 tangent to the undeformed elastic axis 𝑠 at 𝑃0 . Correspondingly, deformed local coordinate 𝑑 includes Cartesian coordinates 𝑨𝑩𝑪, with origin 𝑃 fixed to the elastic center of the displaced cross-section. Deformed local coordinate 𝑑 rotates with the cross-section such that 𝑩 and 𝑪 remain on the cross-sectional plane but 𝑨 is not generally tangent to the deformed elastic axis 𝑆.

𝑊 𝑃0 𝑊 ̇ 𝒗𝐺 + 𝑻̇ 𝑔𝑓 𝑻 𝑓 𝑙 𝒖𝑙 + 𝑻 𝑔𝑓 𝑻 𝑓 𝑙 𝒖̇ 𝑙 𝑔 = 𝒗𝑔 + 𝑻 𝑔𝑓 𝑹 𝑓

+ 𝑻̇ 𝑔𝑓 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆𝑑 + 𝑻 𝑔𝑓 𝑻 𝑓 𝑙 𝑻̇ 𝑙𝑑 𝒆𝑑

(5)

𝑹𝑃𝑑 𝐺

in which 𝒆𝑑 = is the eccentricity measured in coordinate 𝑑, and ̇ 𝑂𝑊 . Time derivatives of 𝑹𝑊 𝑃0 and 𝑻 𝑓 𝑙 are equal to zero be𝒗𝑊 = 𝑹 𝑔 𝑔 𝑓 cause there is no relative motion between the local undeformed coordinates and floating coordinate 𝑓 . Cross-sectional eccentricity is assumed to be constant in time, 𝒆̇ 𝑑 = 𝟎, because cross-sectional deformations are assumed to be small per assumption (a). 2.3. Momentum balance Momentum balance is applicable to the absolute momentum measured in an inertial coordinate system. Any change in momentum results from the sum of forces and moments acting on each differential cross-sectional element: ∑ 𝑖𝑛𝑡 𝑳̇ 𝑔 = 𝑭 𝑔 = 𝑭 𝑒𝑥𝑡 (6) 𝑔 + 𝑭𝑔 ∑ 𝐺 𝐺 𝑒𝑥𝑡,𝐺 𝑖𝑛𝑡,𝐺 ̇ 𝑯𝑔 = 𝑴𝑔 = 𝑴𝑔 + 𝑴𝑔 (7)

2.2. Displacements and kinematics

in which 𝑳𝑔 = 𝑚𝒗𝐺 is the absolute linear momentum; 𝑔 𝐺 𝐺 𝑯 𝑔 = 𝑻 𝑔𝑓 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑰 𝑑 𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴𝑓 +𝑻 𝑔𝑓 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑰 𝐺 𝑑 𝝎𝑑 is the absolute angular momentum about the center of mass 𝐺; 𝑚 is the mass per unit length 𝑒𝑥𝑡 𝑒𝑥𝑡,𝐺 and 𝑰 𝐺 are the 𝑑 is the mass moment of inertia about 𝐺; 𝑭 𝑔 and 𝑴 𝑔 total forces and moments about 𝐺 resulting from external loads; 𝑭 𝑖𝑛𝑡 𝑔 and 𝑴 𝑖𝑛𝑡,𝐺 are the total internal forces and moments about 𝐺 resulting 𝑔 from beam deformation. The left-hand side of Eq. (6) results from taking the time derivative of the absolute translational velocity expressed as Eq. (5):

The absolute displacement measured in global coordinate 𝑔 of any point along the beam can be expressed as a sum of two parts: the displacement resulting from a prescribed overall rigid-body motion of the undeformed beam, plus the unknown displacements resulting from the local beam deformation. The unknown displacements are measured in local undeformed coordinate 𝑙: translational displacements 𝑢𝑖 (𝑠, 𝑡) are the components of a vector from 𝑃0 to 𝑃 resolved onto the 𝒂𝒃𝒄 coordinates; rotational displacements 𝜃𝑖 (𝑠, 𝑡) are defined using three rotational variables that describe the 3D finite rotations from coordinate 𝑙 to 𝑑, in which 𝑠 is the location of 𝑃0 ; 𝑡 is time, and 𝑖 = 1, 2, 3.

𝑊 𝑃0 ̈ 𝑳̇ 𝑔 = 𝑚(𝒂𝑊 + 𝑻̈ 𝑔𝑓 𝑻 𝑓 𝑙 𝒖𝑙 + 2𝑻̇ 𝑔𝑓 𝑻 𝑓 𝑙 𝒖̇ 𝑙 + 𝑻 𝑔𝑓 𝑻 𝑓 𝑙 𝒖̈ 𝑙 𝑔 + 𝑻 𝑔𝑓 𝑹 𝑓

160

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

+ 𝑻̈ 𝑔𝑓 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆𝑑 + 2𝑻̇ 𝑔𝑓 𝑻 𝑓 𝑙 𝑻̇ 𝑙𝑑 𝒆𝑑 + 𝑻 𝑔𝑓 𝑻 𝑓 𝑙 𝑻̈ 𝑙𝑑 𝒆𝑑 ) 𝒂𝑊 𝑔

(8)

𝒗̇ 𝑊 𝑔

in which = is the absolute acceleration of point 𝑊 . Multiplying the combination of Eqs. (6) and (8) by 𝑻 𝑓 𝑔 and substituting Eqs. (1)–(4) into the result yields a physically-intuitive formula that represents the change in linear momentum in non-inertial floating coordinate 𝑓 : 𝑊 𝑃0 𝑖𝑛𝑡 𝑊 ̃̇ 𝑭 𝑒𝑥𝑡 + 𝑻 𝑓 𝑙 𝒖𝑙 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆𝑑 ) 𝑓 + 𝑭 𝑓 = 𝑚[𝒂𝑓 + 𝜴𝑓 (𝑹 𝑓 𝑊 𝑃0

̃ 𝑓𝜴 ̃ 𝑓 (𝑹 +𝜴 𝑓

+ 𝑻 𝑓 𝑙 𝒖𝑙 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆𝑑 )

̃ 𝑑𝝎 ̃ 𝑑 𝒆𝑑 ) + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 (𝝎 ̃̇ 𝑑 𝒆𝑑 + 𝝎 ̃ 𝑓 𝑻 𝑓 𝑙 (𝒖̇ 𝑙 + 𝑻 𝑙𝑑 𝝎 ̃ 𝑑 𝒆𝑑 ) + 𝑻 𝑓 𝑙 𝒖̈ 𝑙 ] + 2𝜴

(9) Fig. 2. Linear strain.

The derivation for angular momentum balance in floating coordinate 𝑓 follows the same logic: computing the time derivative of 𝑯 𝐺 𝑔, 𝐺 multiplying the combination of Eq. (7) and 𝑯̇ by 𝑻 𝑓 𝑔 , and substituting 𝑔

The elastic force and moment of a cross-section at 𝑠 are assumed to be functions of strain measures and/or the time derivatives of strain measures:

Eqs. (1)–(4) into the result yields: ̃ 𝑓 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑰 𝐺 (𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴𝑓 + 𝝎𝑑 ) 𝑴 𝑒𝑥𝑡,𝐺 + 𝑴 𝑖𝑛𝑡,𝐺 = 𝜴 𝑑 𝑓 𝑓 ̃ 𝑑𝑰𝐺 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝝎 𝑑 (𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴𝑓 + 𝝎𝑑 ) + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑰 𝐺 (𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴̇ 𝑓 + 𝝎̇ 𝑑 ) 𝑑 ̃ − 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑰 𝐺 𝑑 𝝎𝑑 𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴𝑓

(10)

The relationships between unknown displacements, 𝒖𝑙 and 𝜽, along the beam and the associated internal forces and moments are needed for solution of Eqs. (9) and (10). Angular strains are computed from curvature vectors; linear strains are computed using differential position vectors. The curvature vector 𝒌𝑙 of the undeformed curvilinear coordinate 𝑠, and 𝒌𝑑 of the deformed curvilinear coordinate 𝑆, each describe the change in direction of a curvilinear coordinate per unit length, measured in the local Cartesian coordinates 𝑙 and 𝑑 respectively. The skew-symmetric matrix of the local curvature vector at any point can be computed as 𝒌̃ 𝑙 = 𝑻 𝑙𝑓 𝑻 ′𝑓 𝑙 or 𝒌̃ 𝑑 = 𝑻 𝑑𝑓 𝑻 ′𝑓 𝑑 , in which 𝑻 𝑑𝑓 = 𝑻 𝑑𝑙 𝑻 𝑙𝑓 and ( )′ denotes spatial derivatives with respect to 𝑠. The local angular strain, 𝜿, can be computed along the beam as the difference between these curvature vectors:

𝜸=

− 𝒓′𝑙

=

𝑃 ′ 𝑻 𝑑𝑓 (𝑹𝑊 𝑓 )

𝑊𝑃 − 𝑻 𝑙𝑓 (𝑹𝑓 0 )′

=

𝑻 𝑑𝑙 (𝒆1 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) − 𝒆1

(14)

′ ′ ′ 𝑭 𝑖𝑛𝑡 𝑓 = (𝑻 𝑓 𝑙 𝑻 𝑙𝑑 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 )𝑭 𝑑 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑭 𝑑

(15)

The total of the internal moments exerted at each end of the element expressed about the mass center, including elastic forces and eccentricities is: 𝑴 𝑖𝑛𝑡,𝐺 = (𝑻 ′𝑓 𝑙 𝑻 𝑙𝑑 + 𝑻 𝑓 𝑙 𝑻 ′𝑙𝑑 )𝑴 𝑃𝑑 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 (𝑴 𝑃𝑑 )′ 𝑓

(11)

− 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆̃ 𝑑 𝒌̃ 𝑑 𝑭 𝑑 − 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆̃ 𝑑 𝑭 ′𝑑 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 (𝜸̃ + 𝒆̃ 1 )𝑭 𝑑

where axial( ) denotes the axial vector computed from a skew-symmetric matrix. Derivation of linear strain follows a similar logic applied to differential position vectors instead of curvature vectors. Differential position vectors 𝒓′𝑙 of the undeformed curvilinear coordinate 𝑠, and 𝒓′𝑑 of the deformed curvilinear coordinate 𝑆, each describe the change in shape of a differential element along the beam axis, measured in the local Cartesian coordinates 𝑙 and 𝑑 respectively. Differential position vectors 𝑊𝑃 𝑃 ′ at a point can be computed as 𝒓′𝑙 = 𝑻 𝑙𝑓 (𝑹𝑓 0 )′ and 𝒓′𝑑 = 𝑻 𝑑𝑓 (𝑹𝑊 𝑓 ) . Transformation matrices 𝑻 𝑙𝑓 and 𝑻 𝑑𝑓 effectively rotate the spatial derivatives of position vectors to the common coordinate system 𝑓 (Fig. 2), such that these vectors can be subtracted. The resulting linear strain, 𝜸, includes axial and shear strains and excludes rigid element rotations: 𝒓′𝑑

(13)

in which 𝑭 𝑑 is the cross-sectional elastic force measured in deformed local coordinate 𝑑, and 𝑴 𝑃𝑑 is the cross-sectional elastic moment about elastic center in coordinate 𝑑. Eqs. (13) and (14) represent the constitutive relations for a material, and are applicable for linear elastic, hyperelasic, or viscoelastic models. Momentum balance equations, (9) and (10), are expressed about the mass center of the differential element. The total internal force exerted on the differential element can be computed directly as the sum of cross-sectional elastic forces, but the total internal moment must be transferred from the elastic center to the mass center including the effects of elastic forces and eccentricities. Internal forces are exerted at each end of a differential element (𝑠− and 𝑠+ ). The total internal force exerted on an element is:

2.4. Strain measures and internal forces and moments

𝜿 = 𝒌𝑑 − 𝒌𝑙 = axial(𝑻 𝑑𝑙 𝒌̃ 𝑙 𝑻 𝑙𝑑 + 𝑻 𝑑𝑙 𝑻 ′𝑙𝑑 − 𝒌̃ 𝑙 )

̇ 𝜿) ̇ 𝑭 𝑑 (𝑠, 𝑡) = 𝑭 𝑑 (𝜸, 𝜿, 𝜸, ̇ 𝜿) ̇ 𝑴 𝑃𝑑 (𝑠, 𝑡) = 𝑴 𝑃𝑑 (𝜸, 𝜿, 𝜸,

(16)

2.5. Equations of motion Translational and rotational EOM’s are derived based on conservation of momentum in non-inertial coordinate 𝑓 by substituting Eqs. (15) and (16) into Eqs. (9) and (10), respectively: 𝑭 𝑒𝑥𝑡 𝑓 =

𝑊 𝑃0 ̃̇ 𝑚[𝒂𝑊 + 𝑻 𝑓 𝑙 𝒖𝑙 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆𝑑 ) 𝑓 + 𝜴𝑓 (𝑹 𝑓 𝑊 𝑃0

̃ 𝑓𝜴 ̃ 𝑓 (𝑹 +𝜴 𝑓

+ 𝑻 𝑓 𝑙 𝒖𝑙 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆𝑑 )

̃ 𝑑𝝎 ̃ 𝑑 𝒆𝑑 ) + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 (𝝎 ̃̇ 𝑑 𝒆𝑑 + 𝝎 ̃ ̃ 𝑑 𝒆𝑑 ) + 𝑻 𝑓 𝑙 𝒖̈ 𝑙 ] + 2𝜴𝑓 𝑻 𝑓 𝑙 (𝒖̇ 𝑙 + 𝑻 𝑙𝑑 𝝎 − (𝑻 ′𝑓 𝑙 𝑻 𝑙𝑑 + 𝑻 𝑓 𝑙 𝑻 ′𝑙𝑑 )𝑭 𝑑 − 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑭 ′𝑑 𝑴 𝑒𝑥𝑡,𝐺 𝑓

(12)

=

(17)

̃ 𝑓 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑰 𝐺 (𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴𝑓 + 𝝎𝑑 ) 𝜴 𝑑 ̃ 𝑑𝑰𝐺 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝝎 𝑑 (𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴𝑓 + 𝝎𝑑 ) + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑰 𝐺 (𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴̇ 𝑓 + 𝝎̇ 𝑑 )

𝑊𝑃 𝑊𝑃 in which 𝒓′𝑙 = 𝒌̃𝑙 𝑹𝑙 0 + (𝑹𝑙 0 )′ = 𝒆1 and 𝒆1 = [1, 0, 0]𝑇 is a unit vector, ′ because 𝒓𝑙 yields a tangent unit vector that is always aligned with the 𝒂direction of local undeformed coordinate 𝑙. Eqs. (11) and (12) represent the strain measures of a beam subject to large angular and translational deflections, and include full consideration of geometrical nonlinearities. These strain measures are invariant with rigid-body motion.

𝑑

̃ − 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝑰 𝐺 𝑑 𝝎𝑑 𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝜴𝑓 − (𝑻 ′𝑓 𝑙 𝑻 𝑙𝑑 + 𝑻 𝑓 𝑙 𝑻 ′𝑙𝑑 )𝑴 𝑃𝑑 − 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 (𝑴 𝑃𝑑 )′ + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆̃ 𝑑 𝒌̃ 𝑑 𝑭 𝑑 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 𝒆̃ 𝑑 𝑭 ′𝑑 − 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 (𝜸̃ + 𝒆̃ 1 )𝑭 𝑑 161

(18)

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

Multiplying the time derivative of Eq. (20) by 𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝑻 𝑓 𝑔 and combining Eqs. (6), (7), (15), and (16) with the resulting equation yields:

Eqs. (17) and (18) are the non-linear 6-DOF EOM’s expressed using transformation matrices in the floating frame. These EOM’s are established as a part of a total Lagrangian formulation; the primary unknown variables are displacements 𝒖𝑙 and 𝜽 relative to non-inertial coordinate 𝑙, which represents the undeformed configuration. A complete, displacement-based formulation of MBBT is represented by the EOM’s of Eqs. (17)–(18), constitutive relations of Eqs. (13)–(14), strain measures of Eqs. (11)–(12), angular motions of Eqs. (1)–(4), and the rotation formulations for transformation matrices. Any three-variable finite rotation formulation can be used to compute transformation matrices 𝑻 𝑙𝑑 and 𝑻 𝑑𝑙 as functions of 𝜽, such as Euler angles, Rodrigues’ formula, or Lie group. Inertial forces and moments associated with the overall rigid-body motion of an undeformed beam are explicitly included in Eqs. (17) and (18). Significant inertial forces and moments generally result from translational acceleration, 𝒂𝑊 , rotational velocity, 𝑓 𝜴𝑓 , and rotational acceleration, 𝜴̇ 𝑓 . The inertial term that includes ̃ 𝑓𝜴 ̃ 𝑓 in Eq. (17) represents the centrifugal force introduced by rigid𝜴 ̃ 𝑓 in Eq. (17) represents body motion. The inertial term that includes 2𝜴 the Coriolis effect. The first term in the right-hand-side of Eq. (18) implies that rigid-body motion alone (𝜴𝑓 ≠ 0, 𝝎𝑑 = 0) can cause an inertial moment in the floating frame. The explicit representation of rigid-body motion in MBBT enables the beam theory to be readily coupled with multibody dynamics in which rigid-body motion and local deformation can be numerically solved either by separate codes or by a coupled routine. The floating coordinate system represents the arbitrary absolute motion of an undeformed beam described by vectors 𝒂𝑊 , 𝜴𝑓 , and 𝑓 𝜴̇ 𝑓 in Eqs. (17)–(18). MBBT enables computation of local dynamics relative to the floating coordinate system. Moreover, Eqs. (17)–(18) can be transferred to global inertial coordinate 𝑔 or to local deformed (material) coordinate 𝑑 through multiplication of the EOM’s by transformation matrices 𝑻 𝑓 𝑔 or 𝑻 𝑓 𝑑 , respectively. No re-derivation of the EOM’s is required because the original EOM’s include full consideration of relative motions among these coordinate systems.

𝑃 ̃ ∗𝑯𝑃 𝑯̇ 𝑑 + 𝜴 𝑑 𝑑 𝐺

= 𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝑻 𝑓 𝑔 (𝒆̃̇ 𝑔 𝑳𝑔 + 𝒆̃ 𝑔 𝑳̇ 𝑔 + 𝑯̇ 𝑔 ) ̃ ∗ 𝒆̃ 𝑑 𝑳𝑑 + 𝒆̃ 𝑑 [𝑭 𝑒𝑥𝑡 + 𝒌̃ 𝑑 𝑭 𝑑 + (𝑭 𝑑 )′ ] + 𝑴 𝑒𝑥𝑡,𝐺 + 𝒌̃ 𝑑 𝑴 𝑃 =𝜴 𝑑 𝑑 𝑑 𝑑 + (𝑴 𝑃 )′ − 𝒆̃ 𝑑 [𝒌̃ 𝑑 𝑭 𝑑 + (𝑭 𝑑 )′ ] + (𝜸̃ + 𝒆̃ 1 )𝑭 𝑑 𝑑

= −𝒗̃ 𝑃𝑑 𝑳𝑑 + 𝑴 𝑒𝑥𝑡,𝑃 + 𝒌̃ 𝑑 𝑴 𝑃𝑑 + (𝑴 𝑃𝑑 )′ + (𝜸̃ + 𝒆̃ 1 )𝑭 𝑑 𝑑 ̃ ∗ 𝒆̃ 𝑑 𝑳𝑑 𝜴 𝑑

The governing equations of the fully non-linear MBBT are linearized based on physical engineering assumptions. The linearized MBBT is then further simplified for application to three practical engineering cases. The linearized formulation is also re-derived in a discretized integral format using the finite-volume method. 3.1. Linearized MBBT Three additional assumptions are introduced for the linearization of MBBT: (a) rotational deformations along the beam are assumed to be small; (b) strain measures are assumed to be small, and (c) constitutive equations are assumed to be linear functions of strain measures. Transformation matrices 𝑻 𝑙𝑑 and 𝑻 𝑑𝑙 and derivatives of these matrices can be simplified considerably if the magnitudes of rotational displacė 𝑡)|, |𝜽(𝑠, ̈ 𝑡)|, and ments and their derivatives are small: |𝜽(𝑠, 𝑡)|, |𝜽(𝑠, |𝜽′ (𝑠, 𝑡)| ≪ 1. The 1-2-3 sequenced Euler angle formulation is selected here to compute coordinate transformation matrices [59]. This Euler angle formulation is straightforward and no singularity issues arise in this application, which is limited to small angles. Alternate rotation formulations could be applied. Assuming sin(𝜃) ≈ 𝜃 and cos(𝜃) ≈ 1 in Euler angle formulation yields the simplified transformation matrices 𝑻 𝑙𝑑 and 𝑻 𝑑𝑙 for small angles:

′ ′ ′ 𝑭 𝑒𝑥𝑡 𝑑 + 𝑻 𝑑𝑙 𝑻 𝑙𝑓 (𝑻 𝑓 𝑙 𝑻 𝑙𝑑 + 𝑻 𝑓 𝑙 𝑻 𝑙𝑑 )𝑭 𝑑 + (𝑭 𝑑 ) 𝑒𝑥𝑡 ′ = 𝑭 + 𝒌̃ 𝑑 𝑭 𝑑 + (𝑭 𝑑 )

𝑻 𝑑𝑙 =

𝑑

(19)

where 𝜴∗𝑑 = 𝜴𝑑 +𝝎𝑑 = 𝑻 𝑑𝑓 𝜴𝑓 +𝝎𝑑 is the absolute angular velocity, and ̃ ∗ 𝒆𝑑 ) is the absolute linear momentum measured in 𝑳𝑑 = 𝑚𝒗𝐺 = 𝑚(𝒗𝑃𝑑 +𝜴 𝑑 𝑑 coordinate 𝑑. The absolute angular momentum about the elastic center can be computed in inertial coordinate 𝑔: ∗ 𝐺 ̃ 𝑯 𝑃𝑔 = 𝑻 𝑔𝑑 𝑯 𝑃𝑑 = 𝑻 𝑔𝑑 (̃𝒆𝑑 𝑳𝑑 + 𝑰 𝐺 𝑑 𝜴𝑑 ) = 𝒆𝑔 𝑳𝑔 + 𝑯 𝑔

𝑻 𝑙𝑑 ≈ 𝑰 + 𝜽̃

(22)

𝑻 𝑇𝑙𝑑

(23)

≈ 𝑰 − 𝜽̃

in which 𝑰 is identity matrix. Linearized equations for the angular velocity and acceleration are obtained by substituting Eqs. (22)–(23) into Eqs. (1) and (3) and neglecting all terms of second-order and higher in 𝜃𝑖 , 𝜃̇ 𝑖 , and 𝜃̈𝑖 : [ ]𝑇 𝝎𝑑 ≈ 𝜃̇ 1 , 𝜃̇ 2 , 𝜃̇ 3 = 𝜽̇ (24) [ ]𝑇 ̈ ̈ ̈ ̈ 𝝎̇ 𝑑 ≈ 𝜃1 , 𝜃2 , 𝜃3 = 𝜽 (25)

𝜕 (𝑻 𝑻 𝑻 𝒗𝐺 ) 𝜕𝑡 𝑔𝑓 𝑓 𝑙 𝑙𝑑 𝑑

∗

(21)

−𝒗̃ 𝑃𝑑 𝑳𝑑

3. Linearization, simplification, and discretization

MBBT is unique amongst other GEBT formulations in that dynamic quantities are relative to non-inertial frames and coordinate transformation matrices are used throughout the derivation. The vector-based formulation and the overall framework of the derivation are both useful for subsequent theoretical developments. The vector-based formulation using the relative dynamic quantities retains physical significance for each term and enables rigorous application of practical engineering assumptions, which leads to rational linearization in Section 3.1 and subsequent simplification in Section 3.2. The overall derivation framework using transformation matrices is essential for the separation of displacements technique developed in Section 4.1. The final non-linear result fo MBBT can be shown to be mathematically equivalent to the existing inertial frame formulation of GEBT, as would be expected because each derivation is based on the same physical assumptions. An inertial frame formulation of MBBT can be developed following the non-inertial frame formulation. Multiplying Eq. (6) by 𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝑻 𝑓 𝑔 and substituting Eq. (15) into the resulting equation yields:

̃ 𝑳𝑑 + 𝑳̇ 𝑑 =𝜴 𝑑

𝒗̃ 𝑃𝑑 )𝑚𝒗𝐺 𝑑

where = − = is the inertial-coupling term that represents the moment about the elastic center introduced by the inertial force acting at the center of mass, which does not exist if the EOM is established about the mass center. Eqs. (19) and (21) are the inertial frame formulation, and are identical to the EOM’s presented as Eqs. (1a) and (1b) in [41], in which dynamic quantities (𝜴∗𝑑 and 𝒗𝑃𝑑 ) and the underlying displacements are relative to an inertial frame. Hodges’ formulation or Eqs. (19) and (21) can also be used as the basis for an alternate derivation of MBBT in a non-inertial frame. This alternate derivation is approximately as complicated as the proposed MBBT derivation, and does not offer the useful intermediate results that makes MBBT valuable. The alternate derivation includes several steps: (a) transfer the reference coordinate into the non-inertial frame using variable transformation matrices; (b) transfer the reference point of angular momentum from the elastic center to the center of cross-sectional mass using Eq. (20); (c) Algebraically re-formulate the resulting dynamic quantities, and then finally (d) derive the noninertial frame formulation using the identity regarding antisymmetric ̃ ̃ operations: 𝜶̃ 𝜷̃ = 𝜷̃ 𝜶̃ + (𝜶𝜷), in which 𝜶 and 𝜷 are arbitrary 3D vectors.

2.6. Relation to existing GEBT

= 𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝑻 𝑓 𝑔 𝑚

(𝒗̃ 𝐺 𝑑

The spatial derivatives of translational displacements are small under small-strain assumption: |𝜿|, |𝜸|, and |𝒖′𝑙 (𝑠, 𝑡)| ≪ 1. The linearized

(20) 162

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

− 𝑻 ′𝑓 𝑙 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )]

angular and linear strain measures are found by substituting Eqs. (22)– (23) into Eqs. (11)–(12) and neglecting all terms of second-order and higher: ′

𝜿 ≈ axial(𝒌̃ 𝑙 𝜽̃ − 𝜽̃ 𝒌̃ 𝑙 + 𝜽̃ ) = 𝒌̃ 𝑙 𝜽 + 𝜽′

(26)

𝜸 ≈ 𝒆̃ 1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙

(27)

′ − 𝑻 𝑓 𝑙 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )]

𝑴 𝑒𝑥𝑡,𝐺 𝑓

− 𝑻 ′𝑓 𝑙 [𝑲 21 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 22 (𝒌̃ 𝑙 𝜽 + 𝜽′ )] ′ − 𝑻 𝑓 𝑙 [𝑲 21 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 22 (𝒌̃ 𝑙 𝜽 + 𝜽′ )] + 𝑻 𝑓 𝑙 𝒆̃ 𝑑 𝒌̃ 𝑙 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′ ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )]

Internal forces and moments are computed using the strain measures and constitutive matrices assuming linear elasticity, such that Eqs. (13) and (14) can be expressed as: [ ] [ ] [ ][ ] 𝑭 𝑑 (𝑠, 𝑡) 𝜸(𝑠, 𝑡) 𝑲 11 𝑲 12 𝜸(𝑠, 𝑡) =𝑲 = (28) 𝑴 𝑑 (𝑠, 𝑡) 𝜿(𝑠, 𝑡) 𝑲 21 𝑲 22 𝜿(𝑠, 𝑡)

𝑙

′

+ 𝑻 𝑓 𝑙 𝒆̃ 𝑑 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )] − 𝑻 𝑓 𝑙 𝒆̃ 1 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )]

′

𝑭 𝑒𝑥𝑡 ≈ 𝑚𝒖̈ 𝑙 − 𝑚̃𝒆𝑑 𝜽̈ − [𝑲 11 (̃𝒆1 𝜽 + 𝒖′𝑙 ) + 𝑲 12 𝜽′ ] 𝑙 ′ ̈ 𝒆1 𝜽 + 𝒖′𝑙 ) + 𝑲 22 𝜽′ ] 𝑴 𝑒𝑥𝑡,𝐺 ≈ 𝑰𝐺 𝑑 𝜽 − [𝑲 21 (̃ 𝑙

𝑊 𝑃0 𝑊 ̃ ̃ ̃̇ 𝑭 𝑒𝑥𝑡 + 𝑻 𝑓 𝑙 𝒆𝑑 ) 𝑓 ≈ 𝑚[𝒂𝑓 + (𝜴𝑓 + 𝜴𝑓 𝜴𝑓 )(𝑹 𝑓

𝑴 𝑒𝑥𝑡,𝐺 𝑓

′

≈

− 𝒆̃ 1 [𝑲 11 (̃𝒆1 𝜽 + 𝒖′𝑙 ) + 𝑲 12 𝜽′ ]

̇ + 𝑻 𝑓𝑙𝑰𝐺 𝑑 𝑻 𝑙𝑓 𝜴𝑓

(29)

𝑭 𝑒𝑥𝑡 ≈ 𝑚𝒖̈ 𝑙 − 𝑲 ′11 (̃𝒆1 𝜽 + 𝒖′𝑙 ) − 𝑲 11 (̃𝒆1 𝜽 + 𝒖′𝑙 )′ 𝑙 𝑴 𝑒𝑥𝑡 𝑙

̈ + 𝑻 𝑓𝑙𝑰𝐺 𝑑𝜽

𝐺̃ ̃ 𝑓 𝑻 𝑓 𝑙 𝑰 𝐺 + 𝑻 𝑓 𝑙 𝑰 𝐺 (𝑻̃ ̇ + [𝜴 𝑙𝑓 𝜴𝑓 ) − 𝑻 𝑓 𝑙 (𝑰 𝑑 𝑻 𝑙𝑓 𝜴𝑓 )]𝜽 𝑑 𝑑

𝐺̃̇ ̃ ̇ + [𝑻 𝑓 𝑙 𝑰 𝐺 𝑑 (𝑻 𝑙𝑓 𝜴𝑓 ) − 𝑻 𝑓 𝑙 (𝑰 𝑑 𝑻 𝑙𝑓 𝜴𝑓 )]𝜽 − 𝑻 ′ [𝑲 21 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′ ) + 𝑲 22 (𝒌̃ 𝑙 𝜽 + 𝜽′ )] 𝑙

′ − 𝑻 𝑓 𝑙 [𝑲 21 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 22 (𝒌̃ 𝑙 𝜽 + 𝜽′ )] + 𝑻 𝑓 𝑙 𝒆̃ 𝑑 𝒌̃ 𝑙 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′ ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )]

𝑚

𝑙

′ + 𝑻 𝑓 𝑙 𝒆̃ 𝑑 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )] − 𝑻 𝑓 𝑙 𝒆̃ 1 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )]

≈

′ ′ ̈ 𝑰𝐺 𝑑 𝜽 − 𝑲 22 𝜽

′′

− 𝑲 22 𝜽

− 𝒆̃ 1 𝑲 11 (̃𝒆1 𝜽 + 𝒖′𝑙 )

(35) (36)

3.2.3. 2D beams having negligible axial deformation The EOM’s represented by Eqs. (35) and (36) can be simplified further if axial force and deformation are negligible. The resulting equations are applicable to straight 2D beams with loading limited to shear and bending in 𝑥𝑦-plane. Eqs. (37) and (38) are commonly used to represent Timoshenko beam theory:

𝐺̃ ̃ 𝑓 𝑻 𝑓 𝑙 𝑰 𝐺 (𝑻̃ ̃ + [𝜴 𝑙𝑓 𝜴𝑓 ) − 𝜴𝑓 𝑻 𝑓 𝑙 (𝑰 𝑑 𝑻 𝑙𝑓 𝜴𝑓 )]𝜽 𝑑

𝑓𝑙

(34)

Eqs. (33) and (34) can be simplified further if the beam is symmetric such that eccentricity is zero (𝒆𝑑 = 𝟎):

𝑙

̃ 𝑓 𝑻 𝑓 𝑙 𝑰 𝐺 𝑻 𝑙𝑓 𝜴𝑓 𝜴 𝑑

(33)

+ 𝒆̃ 𝑑 [𝑲 11 (̃𝒆1 𝜽 + 𝒖′𝑙 ) + 𝑲 12 𝜽′ ]

̃ 𝑓𝜴 ̃ 𝑓 )𝑻 𝑓 𝑙 𝒖𝑙 + 2𝜴 ̃ 𝑓 𝑻 𝑓 𝑙 𝒖̇ 𝑙 + 𝑻 𝑓 𝑙 𝒖̈ 𝑙 ̃̇ 𝑓 + 𝜴 + (𝜴 ̇ ̃ 𝑓𝜴 ̃ 𝑓 )𝑻 𝑓 𝑙 𝒆̃ 𝑑 𝜽 − 2𝜴 ̃ 𝑓 𝑻 𝑓 𝑙 𝒆̃ 𝑑 𝜽̇ − 𝑻 𝑓 𝑙 𝒆̃ 𝑑 𝜽] ̈ ̃ +𝜴 − (𝜴 ′ ′ ′ ̃ ̃ − 𝑻 [𝑲 11 (̃𝒆1 𝜽 + 𝒌𝑙 𝒖𝑙 + 𝒖 ) + 𝑲 12 (𝒌𝑙 𝜽 + 𝜽 )] ′ − 𝑻 𝑓 𝑙 [𝑲 11 (̃𝒆1 𝜽 + 𝒌̃ 𝑙 𝒖𝑙 + 𝒖′𝑙 ) + 𝑲 12 (𝒌̃ 𝑙 𝜽 + 𝜽′ )]

(32)

3.2.2. Straight beams without initial curve or twist The EOM’s for beams in an inertial coordinate system can be further simplified if the beam is straight, such that local undeformed coordinate 𝑙 is aligned with Cartesian coordinate system 𝑥𝑦𝑧. Substituting 𝑻 𝑓 𝑙 = 𝑰 and 𝒌𝑙 = 𝟎 into Eqs. (31) and (32) yields:

where 𝑲 is the cross-sectional constitutive matrix, which is constant in time. The constitutive matrix can be obtained using one of several existing cross-sectional analysis methods, some of which enable inclusion of small cross-sectional deformations [19–21]. The linearized EOM’s of MBBT are derived by substituting Eqs. (22)– (28) into the non-linear Eqs. (17)–(18) and neglecting the terms of second-order and higher in those quantities assumed to be small:

𝑓𝑙

(31)

̈ ≈ 𝑻 𝑓𝑙𝑰𝐺 𝑑𝜽

𝜕 2 𝑢𝑦

𝜕𝑡2 𝜕 2 𝜃𝑧

𝜕𝑢𝑦 𝜕 [𝐾 ( − 𝜃𝑧 )] 𝜕𝑥 1 𝜕𝑥 𝜕𝑢𝑦 𝜕𝜃𝑧 𝜕 (𝐾 ) + 𝐾1 ( − 𝜃𝑧 ) − 𝑀𝑧𝑒𝑥𝑡 (𝑥, 𝑡) ≈ 𝜕𝑥 2 𝜕𝑥 𝜕𝑥 − 𝐹𝑦𝑒𝑥𝑡 (𝑥, 𝑡) ≈

(37)

(38) 𝜕𝑡2 in which 𝐼𝑧𝑧 is moment of inertia per unit length; 𝑥 is the beam axis; 𝐾1 = 𝜅𝐴𝐺 and 𝐾2 = 𝐸𝐼, where 𝜅 is Timoshenko shear coefficient, 𝐴 is the cross-section area, 𝐺 is the shear modulus, 𝐸 is the elastic modulus, and 𝐼 is the second moment of area. 𝐼𝑧𝑧

(30)

Eqs. (29) and (30) are the linearized translational and rotational EOM’s of MBBT resolved in non-inertial coordinate 𝑓 . No additional equations are required to solve for the unknown displacements 𝒖𝑙 and 𝜽, because the linearized rotation formulations and strain measures have been substituted into the EOM’s. Significant engineering insights into moving beam structures can be inferred from Eqs. (29) and (30): structural coupling between translational and rotational responses is captured by strain measures in the elastic terms, and dynamic coupling between translation and rotation is captured by eccentricity in the inertial terms.

3.3. Spatial discretization A spatially discretized model of linearized MBBT is derived to demonstrate that conventional discretization methods can be applied directly to the linearized formulation without special treatment. The finite-volume method is selected because the discretized formulation can be established easily for a control finite volume, as MBBT is based on a conservation law. The discretization could alternatively be done using the FE method, which is more commonly used for structural analyses, but the discretization would be slightly more complicated. A finite control volume attached to a deformed beam segment is shown in Fig. 3. The 𝑗th and (𝑗 + 1)th cross-sections are end boundaries of the 𝑛th segment. The balance of absolute linear momentum is derived for the 𝑛th finite volume in integral form:

3.2. Simplifications for practical engineering cases Three reduced formulations of linearized MBBT are derived for practical engineering applications. Linearized MBBT is shown to be equivalent to Timoshenko beam theory for straight, non-twisted beams in an inertial coordinate system. Linearized MBBT also can be recognized as an extension of Timoshenko beam theory for more complicated beam structures that have eccentricities, irregular shapes, fully-populated stiffness matrices, or overall rigid-body motion.

𝑠𝑗+1

𝑑 𝑑𝑡 ∫𝑠𝑗

3.2.1. Beams in inertial frames Coordinate 𝑓 becomes an inertial system if 𝒂𝑊 = 𝟎, 𝜴𝑓 = 𝟎, and 𝑓 ̇𝜴𝑓 = 𝟎. The EOM’s for beams in an inertial coordinate system can be simplified from Eqs. (29) and (30):

𝑛,𝑒𝑥𝑡 𝑚(𝑠)𝒗𝐺 + 𝑭 𝑛,𝑖𝑛𝑡 𝑔 𝑔 (𝑠, 𝑡)𝑑𝑠 = 𝑭 𝑔

(39)

and 𝑭 𝑛,𝑖𝑛𝑡 are the total external and internal forces in which 𝑭 𝑛,𝑒𝑥𝑡 𝑔 𝑔 exerted on the 𝑛th segment, respectively. The total internal force, 𝑭 𝑛,𝑖𝑛𝑡 𝑔 , is the sum of elastic forces applied at the 𝑗th and (𝑗 +1)th cross-sections, which can be computed using Eq. (28). There is no momentum flux

̈ ̃ ̈ 𝑭 𝑒𝑥𝑡 𝑓 ≈ 𝑚𝑻 𝑓 𝑙 𝒖𝑙 − 𝑚𝑻 𝑓 𝑙 𝒆𝑑 𝜽 163

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

are the cross-sectional constitutive matrices of the 𝑗th cross-section, and vectors 𝜸 𝑗 and 𝜿 𝑗 are the linear and angular strains at 𝑠𝑗 . Crosssectional constitutive matrices, 𝑀 𝑛 , 𝑹𝑛𝑓 , 𝑻 𝑛𝑓 𝑙 , 𝑻 𝑛𝑗 , 𝑻 𝑛𝑗+1 , and 𝒆𝑛 are 𝑙 𝑙 known constants related only to the undeformed beam properties. A deformed local coordinate, 𝑑 𝑛 , is introduced with origin fixed to 𝑃 𝑛 and Cartesian coordinates 𝑨𝑛 𝑩 𝑛 𝑪 𝑛 (Fig. 3). The translation and rotation of coordinate 𝑑 𝑛 are specified by the averaged displacements, 𝒖𝑛𝑙 and 𝜽𝑛 . The mass moment of inertia of the 𝑛th finite control volume, 𝑰 𝑛𝑑 𝑛 , is evaluated in coordinate 𝑑 𝑛 about the center of mass, 𝐺𝑛 . The mass moment 𝑰 𝑛𝑑 𝑛 is generally a time-dependent parameter that has a non-linear relationship with the instantaneous deformation of the 𝑛th segment, but which can reasonably be assumed constant for most structural applications having relatively small finite segments and small strains. The balance of absolute angular momentum for the 𝑛th segment similarly is derived in a linear, discretized form by setting 𝑰 𝑛𝑑 𝑛 to a constant value of 𝑰 𝑛 : ̃ 𝑓 𝑻 𝑛 𝑰 𝑛 𝑻 𝑛 𝜴𝑓 + 𝑻 𝑛 𝑰 𝑛 𝑻 𝑛 𝜴̇ 𝑓 + 𝑻 𝑛 𝑰 𝑛 𝜽̈ 𝑛 𝜴 𝑓𝑙 𝑓𝑙 𝑙𝑓 𝑓𝑙 𝑙𝑓

Fig. 3. The 𝑛th beam finite control volume.

𝑛 𝑛 𝑛̃ 𝑛 ̃ 𝑓 𝑻 𝑛 𝑰 𝑛 + 𝑻 𝑛 𝑰 𝑛 (𝑻̃ ̇𝑛 + [𝜴 𝑓𝑙 𝑓𝑙 𝑙𝑓 𝜴𝑓 ) − 𝑻 𝑓 𝑙 (𝑰 𝑻 𝑙𝑓 𝜴𝑓 )]𝜽 𝑛 𝑛 𝑛 𝑛̃ 𝑛 ̃ 𝑓 𝑻 𝑛 𝑰 𝑛 (𝑻̃ ̃ + [𝜴 𝑓𝑙 𝑙𝑓 𝜴𝑓 ) − 𝜴𝑓 𝑻 𝑓 𝑙 (𝑰 𝑻 𝑙𝑓 𝜴𝑓 )]𝜽

term in the right-hand side of Eq. (39) because the control volume itself moves and deforms with the beam segment. Substituting Eq. (5) into Eq. (39), multiplying the resulting equation by 𝑻 𝑓 𝑔 , and then applying the small-angle and small-strain assumptions (Eqs. (22)–(27)) yields:

𝑛 𝑛 𝑛 ̇ 𝑛̃ 𝑛 ̇ + [𝑻 𝑛𝑓 𝑙 𝑰 𝑛 (𝑻̃ 𝑙𝑓 𝜴𝑓 ) − 𝑻 𝑓 𝑙 (𝑰 𝑻 𝑙𝑓 𝜴𝑓 )]𝜽 𝑃𝑗𝑃𝑛

̃ ≈ 𝑴 𝑛,𝑒𝑥𝑡 + 𝑻 𝑛𝑓 𝑙 (̃𝒆𝑛 + 𝑹 𝑙 𝑓 𝑃 𝑛 𝑃 𝑗+1

𝑛 𝑛 𝑛 ̃ ̃ ̃̇ 𝑀 𝑛 [𝒂𝑊 𝑓 + (𝜴𝑓 + 𝜴𝑓 𝜴𝑓 )(𝑹 𝑓 + 𝑻 𝑓 𝑙 𝒆 )

̃ − 𝑻 𝑛𝑓 𝑙 (̃𝒆𝑛 − 𝑹 𝑙

̃ 𝑓𝜴 ̃ 𝑓 )𝑻 𝑛 𝒖𝑛 + 2𝜴 ̃ 𝑓 𝑻 𝑛 𝒖̇ 𝑛 + 𝑻 𝑛 𝒖̈ 𝑛 ̃̇ 𝑓 + 𝜴 + (𝜴 𝑓𝑙 𝑙 𝑓𝑙 𝑙 𝑓𝑙 𝑙

ment about the center of mass; 𝑲 𝑗21 and 𝑲 𝑗22 are the cross-sectional 𝑗

𝑠𝑗+1

4. Piecewise-linear MBBT

𝑠𝑗+1 ∫𝑠𝑗

𝑚𝑹

𝑊 𝑃0

A special formulation of MBBT is developed that allows application of linearized dynamic theory to large non-linear structural displacements. Two critical advantages of the new MBBT over exiting GEBT formulations are the capacity for rational linearization of the EOM’s and for expression of the EOM’s relative to arbitrary non-inertial frames. The combined strength of these two capabilities enables development of a piecewise-linear version of MBBT using a separation of displacements technique. The piecewise-linear MBBT can be implemented in the timedomain when the incremental rotation within each time-step conforms to small-angle assumptions.

𝑑𝑠

𝑀𝑛

𝑻 𝑛𝑓 𝑙 is the transformation matrix between coordinates 𝑓 and 𝑙 averaged over the 𝑛th segment: 𝑠𝑗+1

𝑻 𝑛𝑓 𝑙 =

∫𝑠𝑗

𝑚𝑻 𝑓 𝑙 𝑑𝑠

𝑀𝑛 𝒆𝑛 is the eccentricity averaged over cross-sections along the undeformed 𝑛th segment: 𝑠𝑗+1

𝑛

𝒆 =

∫𝑠𝑗

4.1. Theory 𝑚𝑻 𝑓 𝑙 𝒆𝑑 𝑑𝑠

𝑠𝑗+1

∫𝑠𝑗

The theoretical basis of separation of displacements is to express the total deformation at any time as a sum of a quasi-static base configuration plus dynamic displacements relative to that base. This separation of displacements is a purely mathematical manipulation: an arbitrary beam configuration without any particular physical significance can be selected as the base configuration and considered to be quasi-static; the time derivatives of the dynamic displacements represent the total velocities and accelerations resulting from deformation. Selection of a base configuration such that dynamic rotations are sufficiently small for application of small-angle assumptions allows linearized MBBT to be used to solve for dynamic displacements, while nonlinearities associated with large displacements can be preserved for the base configuration. The total translational and rotational displacements are each separated into a quasi-static displacement plus a dynamic displacement:

𝑚𝑻 𝑓 𝑙 𝑑𝑠

𝒖𝑛𝑙

is the unknown translational displacement averaged over the 𝑛th segment: 𝑠𝑗+1

𝒖𝑛𝑙 =

∫𝑠𝑗

𝑚𝑻 𝑓 𝑙 𝒖𝑙 𝑑𝑠

𝑠𝑗+1 ∫𝑠𝑗

𝑚𝑻 𝑓 𝑙 𝑑𝑠

𝑛

𝜽 is the unknown rotational displacement averaged over the 𝑛th segment: 𝑠𝑗+1

𝜽𝑛 = 𝑭 𝑛,𝑒𝑥𝑡 𝑓

𝑛 𝑗+1

𝑚𝑑𝑠

∫𝑠𝑗

𝑹𝑛𝑓 is the position vector pointing from 𝑊 to the elastic center 𝑃0𝑛 averaged over the 𝑛th segment in the undeformed configuration: 𝑹𝑛𝑓 =

𝑛

constitutive matrices of the 𝑗th cross-section, and 𝑹𝑙𝑃 𝑃 and 𝑹𝑃𝑙 𝑃 are vectors pointing from elastic centers 𝑃 𝑗 to 𝑃 𝑛 and from 𝑃 𝑛 to 𝑃 𝑗+1 , respectively. Eqs. (40) and (41) are discretized finite-volume formulations of the EOM’s of linearized MBBT.

(40)

in which 𝑀 𝑛 is the total mass of the 𝑛th segment: 𝑀𝑛 =

(41)

where 𝑴 𝑛,𝑒𝑥𝑡 is the sum of external moments acting on the 𝑛th seg𝑓

≈ 𝑭 𝑛,𝑒𝑥𝑡 − 𝑻 𝑛𝑓 𝑙 𝑻 𝑛𝑗 (𝑲 𝑗11 𝜸 𝑗 + 𝑲 𝑗12 𝜿 𝑗 ) 𝑓 𝑙 + 𝑲 𝑗+1 𝜿 𝑗+1 ) 12

)𝑻 𝑛𝑗+1 (𝑲 𝑗+1 𝜸 𝑗+1 + 𝑲 𝑗+1 𝜿 𝑗+1 ) 𝑙 11 12

− 𝑻 𝑛𝑓 𝑙 𝑻 𝑛𝑗 (𝑲 𝑗21 𝜸 𝑗 + 𝑲 𝑗22 𝜿 𝑗 ) + 𝑻 𝑛𝑓 𝑙 𝑻 𝑛𝑗+1 (𝑲 𝑗+1 𝜸 𝑗+1 + 𝑲 𝑗+1 𝜿 𝑗+1 ) 𝑙 𝑙 21 22

̃ 𝑓𝜴 ̃ 𝑓 )𝑻 𝑛 𝒆̃ 𝑛 𝜽𝑛 − 2𝜴 ̃ 𝑓 𝑻 𝑛 𝒆̃ 𝑛 𝜽̇ 𝑛 − 𝑻 𝑛 𝒆̃ 𝑛 𝜽̈ 𝑛 ] ̃̇ + 𝜴 − (𝜴 𝑓𝑙 𝑓𝑙 𝑓𝑙 +𝑻 𝑛𝑓 𝑙 𝑻 𝑛𝑗+1 (𝑲 𝑗+1 𝜸 𝑗+1 𝑙 11

)𝑻 𝑛𝑗 (𝑲 𝑗11 𝜸 𝑗 + 𝑲 𝑗12 𝜿 𝑗 ) 𝑙

∫𝑠𝑗

𝑠𝑗+1 ∫𝑠𝑗

𝑚𝑻 𝑓 𝑙 𝜽𝑑𝑠 𝑚𝑻 𝑓 𝑙 𝑑𝑠

is the sum of external forces exerted on the 𝑛th segment measured

in coordinate 𝑓 ; 𝑻 𝑛𝑗 and 𝑻 𝑛𝑗+1 are the transformation matrices among 𝑙 𝑙 different undeformed local coordinates at 𝑠𝑗 , 𝑠𝑛 , and 𝑠𝑗+1 ; 𝑲 𝑗11 and 𝑲 𝑗12

𝒖𝑙 (𝑠, 𝑡) = 𝑼 𝑙 (𝑠) + 𝒖∗𝑙 (𝑠, 𝑡) 164

(42)

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

where 𝜿 𝑞 = 𝒌𝑞 − 𝒌𝑙 is the angular strain caused by deflection of base configuration. The linear strain is treated in a similar fashion: substituting Eq. (44) into Eq. (12) assuming small dynamic angles. Linearizing the resulting equation yields an expression for strain linearized about the base configuration: ′ 𝜸 ≈ 𝜸 𝑞 + (𝜸̃ 𝑞 + 𝒆̃ 1 )𝜽∗ + 𝒌̃ 𝑞 𝒖∗𝑞 + 𝒖∗𝑞

(49) 𝑼 ′𝑙 )

c in which 𝜸 𝑞 = 𝑻 𝑞𝑙 (𝒆1 + 𝒌̃ 𝑙 𝑼 𝑙 + − 𝒆1 is the strain caused by base configuration deflections. Substituting Eqs. (44)–(46) into the non-linear EOM’s of Eqs. (17)– (18) yields the governing EOM’s including separation of displacements. The resulting EOM’s are linearized by applying small-angle assumptions, substituting (48)–(49), and then neglecting the second- and higher-order terms in dynamic displacements: 𝑊 𝑃0

𝐵 ̃ ̃ ̃̇ 𝑭 𝑒𝑥𝑡 𝑓 ≈ 𝑚[𝒂𝑓 + (𝜴 + 𝜴𝑓 𝜴𝑓 )(𝑹 𝑓

Fig. 4. Base configuration and coordinate systems.

𝝓𝑙 (𝑠, 𝑡) = 𝜱𝑙 (𝑠) + 𝝓∗𝑙 (𝑠, 𝑡)

̃ 𝑓 )𝑻 𝑓 𝑞 𝒖∗ + 2𝜴 ̃ 𝑓 𝑻 𝑓 𝑞 𝒖̇ ∗ + 𝑻 𝑓 𝑞 𝒖̈ ∗ ̃ 𝑓𝜴 ̃̇ 𝑓 + 𝜴 + (𝜴 𝑞 𝑞 𝑞 ∗ ̇ ̃ 𝑓𝜴 ̃ 𝑓 )𝑻 𝑓 𝑞 𝒆̃ 𝑑 𝜽 − 2𝜴 ̃ 𝑓 𝑻 𝑓 𝑞 𝒆̃ 𝑑 𝜽̇ ∗ − 𝑻 𝑓 𝑞 𝒆̃ 𝑑 𝜽̈ ∗ ] ̃𝑓 +𝜴 − (𝜴 ′ − 𝑻 ′𝑓 𝑞 [𝑲 11 (𝜸 𝑞 + 𝜸̃ 𝑞 𝜽∗ + 𝒆̃ 1 𝜽∗ + 𝒌̃ 𝑞 𝒖∗𝑞 + 𝒖∗𝑞 )

(43)

′

+ 𝑲 12 (𝒌𝑞 − 𝒌𝑙 + 𝒌̃ 𝑞 𝜽∗ + 𝜽∗ )]

in which 𝑼 𝑙 (𝑠) is the quasi-static portion of the translational displacement measured from undeformed coordinate 𝑙; 𝒖∗𝑙 (𝑠, 𝑡) is the dynamic part of the translational displacement; total rotation vector 𝝓𝑙 (𝑠, 𝑡) describes the rotation of a differential element from its undeformed to deformed configurations; rotation vector 𝜱𝑙 (𝑠) represents the rotation from the undeformed configuration to the base configuration, and vector 𝝓∗𝑙 (𝑠, 𝑡) is the dynamic rotation vector. The base configuration is represented by an additional curvilinear coordinate system denoted as coordinate 𝑄 (Fig. 4). Coordinate 𝑄 shares the same origin as curvilinear coordinates 𝑠 and 𝑆, but its curvilinear axis is through each cross-sectional elastic center for a prescribed base configuration. Coordinate 𝑄 is displaced from 𝑠 by constant displacement 𝑼 𝑙 (𝑠). A local element-fixed Cartesian coordinate system, 𝒂∗ 𝒃∗ 𝒄 ∗ , is introduced as coordinate 𝑞, with origin at the elastic center 𝑃 ∗ and with directions described by rotation vector 𝜱𝑙 (𝑠) (Fig. 4). The transformation matrices between coordinate 𝑞 and coordinates 𝑙 and 𝑑 are given by 𝑻 𝑙𝑑 = 𝑻 𝑙𝑞 𝑻 𝑞𝑑 and 𝑻 𝑑𝑙 = 𝑻 𝑑𝑞 𝑻 𝑞𝑙 . Eq. (42) is rewritten as: 𝒖𝑙 (𝑠, 𝑡) = 𝑼 𝑙 (𝑠) + 𝑻 𝑙𝑞 (𝑠)𝒖∗𝑞 (𝑠, 𝑡) 𝒖∗𝑞

′ − 𝑻 𝑓 𝑞 [𝑲 11 (𝜸 𝑞 + 𝜸̃ 𝑞 𝜽∗ + 𝒆̃ 1 𝜽∗ + 𝒌̃ 𝑞 𝒖∗𝑞 + 𝒖∗𝑞 ) ′ ′ +𝑲 12 (𝒌𝑞 − 𝒌𝑙 + 𝒌̃ 𝑞 𝜽∗ + 𝜽∗ )]

𝐺 ̃ ̃ 𝑓 𝑻 𝑓 𝑞 𝑰 𝐺 − 𝑻 𝑓 𝑞 (𝑰 𝐺̃ ̇∗ + [𝜴 𝑑 𝑑 𝑻 𝑞𝑓 𝜴𝑓 ) + 𝑻 𝑓 𝑞 𝑰 𝑑 (𝑻 𝑞𝑓 𝜴𝑓 )]𝜽 ∗ 𝐺̃ ̃ 𝑓 𝑻 𝑓 𝑞 𝑰 𝐺 (𝑻̃ ̃ + [𝜴 𝑞𝑓 𝜴𝑓 ) − 𝜴𝑓 𝑻 𝑓 𝑞 (𝑰 𝑑 𝑻 𝑞𝑓 𝜴𝑓 )]𝜽 𝑑 ∗ 𝐺̃̇ ̃ ̇ + [𝑻 𝑓 𝑞 𝑰 𝐺 𝑑 (𝑻 𝑞𝑓 𝜴𝑓 ) − 𝑻 𝑓 𝑞 (𝑰 𝑑 𝑻 𝑞𝑓 𝜴𝑓 )]𝜽 ′

− 𝑻 ′𝑓 𝑞 [𝑲 21 (𝜸 𝑞 + 𝜸̃ 𝑞 𝜽∗ + 𝒆̃ 1 𝜽∗ + 𝒌̃ 𝑞 𝒖∗𝑞 + 𝒖∗𝑞 ) ∗′

+ 𝑲 22 (𝒌𝑞 − 𝒌𝑙 + 𝒌̃ 𝑞 𝜽∗ + 𝜽 )] ′

− 𝑻 𝑓 𝑞 [𝑲 21 (𝜸 𝑞 + 𝜸̃ 𝑞 𝜽∗ + 𝒆̃ 1 𝜽∗ + 𝒌̃ 𝑞 𝒖∗𝑞 + 𝒖∗𝑞 ) ′

′ + 𝑻 𝑓 𝑞 𝒆̃ 𝑑 𝒌̃ 𝑞 [𝑲 11 (𝜸 𝑞 + 𝜸̃ 𝑞 𝜽∗ + 𝒆̃ 1 𝜽∗ + 𝒌̃ 𝑞 𝒖∗𝑞 + 𝒖∗𝑞 ) ′ + 𝑲 12 (𝒌𝑞 − 𝒌𝑙 + 𝒌̃ 𝑞 𝜽∗ + 𝜽∗ )]

is the dynamic translational displacement = in which measured in coordinate 𝑞. Taking time derivatives of Eq. (44) yields the translational velocity and acceleration:

′

+ 𝑻 𝑓 𝑞 𝒆̃ 𝑑 [𝑲 11 (𝜸 𝑞 + 𝜸̃ 𝑞 𝜽∗ + 𝒆̃ 1 𝜽∗ + 𝒌̃ 𝑞 𝒖∗𝑞 + 𝒖∗𝑞 ) ′

′ − 𝑻 𝑓 𝑞 (𝜸̃ 𝑞 + 𝒆̃ 1 )[𝑲 11 (𝜸 𝑞 + 𝜸̃ 𝑞 𝜽∗ + 𝒆̃ 1 𝜽∗ + 𝒌̃ 𝑞 𝒖∗𝑞 + 𝒖∗𝑞 )

𝜕 (𝑼 + 𝑻 𝑙𝑞 𝒖∗𝑞 ) = 𝑻 𝑙𝑞 𝒖̇ ∗𝑞 (45) 𝜕𝑡 𝑙 𝜕2 𝒂𝑙 = (𝑼 𝑙 + 𝑻 𝑙𝑞 𝒖∗𝑞 ) = 𝑻 𝑙𝑞 𝒖̈ ∗𝑞 (46) 𝜕𝑡2 where 𝑼̇ 𝑙 = 𝑼̈ 𝑙 = 𝟎, and 𝑻̇ 𝑙𝑞 = 𝑻̈ 𝑙𝑞 = 𝟎. Application of small-angle assumptions to Euler angles of the dynamic rotational displacements (𝜃1∗ , 𝜃2∗ , 𝜃3∗ ) measured in coordinate 𝑞 ∗ ∗ ∗ ∗ yields: 𝝓∗𝑞 ≈ 𝜽∗ ; 𝑻 𝑞𝑑 ≈ 𝑰 + 𝜽̃ ; 𝑻 𝑑𝑞 ≈ 𝑰 − 𝜽̃ ; 𝝎𝑑 ≈ 𝜽̇ , and 𝝎̇ 𝑑 ≈ 𝜽̈ . Computation of strain measures necessarily is based on the total deformation including both quasi-static and dynamic displacements. The strain measures associated with the base configuration are computed using non-linear strain and Euler angle formulation. The strain measures associated with dynamic displacements can be linearized about the base configuration by applying small-angle assumptions to the dynamic rotations. The total curvature vector is computed as:

′ + 𝑲 12 (𝒌𝑞 − 𝒌𝑙 + 𝒌̃ 𝑞 𝜽∗ + 𝜽∗ )]

(51)

in which 𝑻 𝑓 𝑞 = 𝑻 𝑓 𝑙 𝑻 𝑙𝑞 . Eqs. (50) and (51) are a general Lagrangian formulation of the EOM’s of MBBT that are linearized about a prescribed base configuration using separation of displacements. Geometric and kinematic nonlinearities are fully preserved for the base configuration. Inertial properties 𝑚 and 𝑰 𝐺 𝑑 can be evaluated either for the undeformed configuration, or evaluated for the base configuration such that inertial nonlinearities associated with base deformations are also preserved. 4.2. Discretized matrix formulation A discretized finite-volume representation of beam motions about a prescribed base configuration is derived from Eqs. (50) and (51) by applying the finite-volume method, following the same steps as in Section 3.3. The resulting linear discretized EOM’s are:

(47)

in which 𝒌̃ 𝑞 = 𝑻 𝑞𝑓 𝑻 ′𝑓 𝑞 = 𝑻 𝑞𝑙 𝑻 𝑙𝑓 𝑻 ′𝑓 𝑙 𝑻 𝑙𝑞 + 𝑻 𝑞𝑙 𝑻 ′𝑙𝑞 is the skew-symmetric matrix of the local curvature vector of the base configuration measured in coordinate 𝑞. Angular strain is linearized about the base configuration and is computed by substituting Eq. (47) into Eq. (11): ′

′

+𝑲 12 (𝒌𝑞 − 𝒌𝑙 + 𝒌̃ 𝑞 𝜽∗ + 𝜽∗ )]

𝒗𝑙 =

𝜿 = axial(𝒌̃ 𝑑 − 𝒌̃ 𝑙 ) ≈ 𝜿 𝑞 + 𝒌̃ 𝑞 𝜽∗ + 𝜽∗

′

+𝑲 22 (𝒌𝑞 − 𝒌𝑙 + 𝒌̃ 𝑞 𝜽∗ + 𝜽∗ )]

[𝑢∗1 , 𝑢∗2 , 𝑢∗3 ]𝑇

′

(50)

̃ 𝑓 𝑻 𝑓 𝑞 𝑰 𝐺 𝑻 𝑞𝑓 𝜴𝑓 + 𝑻 𝑓 𝑞 𝑰 𝐺 𝑻 𝑞𝑓 𝜴̇ 𝑓 + 𝑻 𝑓 𝑞 𝑰 𝐺 𝜽̈ ∗ 𝑴 𝑒𝑥𝑡,𝐺 ≈ 𝜴 𝑑 𝑑 𝑑 𝑓

(44)

∗ ∗ ∗ 𝒌̃ 𝑑 = 𝑻 𝑑𝑙 𝑻 𝑙𝑓 𝑻 ′𝑓 𝑙 𝑻 𝑙𝑑 + 𝑻 𝑑𝑙 𝑻 ′𝑙𝑑 ≈ 𝒌̃ 𝑞 − 𝜽̃ 𝒌̃ 𝑞 + 𝒌̃ 𝑞 𝜽̃ + 𝜽̃

+ 𝑻 𝑓 𝑙 𝑼 𝑙 + 𝑻 𝑓 𝑞 𝒆𝑑 )

𝑛 𝑛 𝑛 𝑛 𝑛 ̃ ̃ ̃̇ 𝑀 𝑛 [𝒂𝑊 𝑓 + (𝜴 + 𝜴𝑓 𝜴𝑓 )(𝑹 𝑓 + 𝑻 𝑓 𝑙 𝑼 𝑙 + 𝑻 𝑓 𝑞 𝒆 )

̃ 𝑓𝜴 ̃ 𝑓 )𝑻 𝑛 𝒖∗,𝑛 + 2𝜴 ̃ 𝑓 𝑻 𝑛 𝒖̇ ∗,𝑛 + 𝑻 𝑛 𝒖̈ ∗,𝑛 ̃̇ 𝑓 + 𝜴 + (𝜴 𝑓𝑞 𝑞 𝑓𝑞 𝑞 𝑓𝑞 𝑞 ̃ 𝑓𝜴 ̃ 𝑓 )𝑻 𝑛 𝒆̃ 𝑛 𝜽∗,𝑛 − 2𝜴 ̃ 𝑓 𝑻 𝑛 𝒆̃ 𝑛 𝜽̇ ∗,𝑛 − 𝑻 𝑛 𝒆̃ 𝑛 𝜽̈ ∗,𝑛 ] ̃̇ 𝑓 + 𝜴 − (𝜴 𝑓𝑞 𝑓𝑞 𝑓𝑞 𝑗 𝑗 𝑗 𝑗 ≈ 𝑭 𝑛,𝑒𝑥𝑡 − 𝑻 𝑛𝑓 𝑞 𝑻 𝑛𝑗 𝑞 (𝑲 11 𝜸 + 𝑲 12 𝜿 ) 𝑓

(48) 165

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

+ 𝑻 𝑛𝑓 𝑞 𝑻 𝑛𝑗+1 (𝑲 𝑗+1 𝜸 𝑗+1 + 𝑲 𝑗+1 𝜿 𝑗+1 ) 𝑞 11 12

(52)

The first type captures conservative inertial effects associated with the [ ] overall rigid-body motion of the beam, 𝑪 𝑛𝑖𝑛𝑒 , which can introduce either positive or negative damping. The second type is traditional [ ] nonconservative structural damping, 𝑪 𝑛𝑠𝑡𝑟 . Total effective damping can be assumed to be the sum of these two types, with the latter type proportional to the corresponding stiffness and mass matrices: [ ] ] [ [ ] ] [ ] [ ] [ (58) 𝑪 𝑛𝑒𝑞𝑢 = 𝑪 𝑛𝑖𝑛𝑒 + 𝑪 𝑛𝑠𝑡𝑟 = 𝑪 𝑛𝑖𝑛𝑒 + 𝜇𝐾 𝑲 𝑛𝑒𝑙𝑎 + 𝜇𝑀 𝑴 𝑛𝑒𝑞𝑢

̃ 𝑓 𝑻 𝑛 𝑰 𝑛 𝑻 𝑛 𝜴𝑓 + 𝑻 𝑛 𝑰 𝑛 𝑻 𝑛 𝜴̇ 𝑓 + 𝑻 𝑛 𝑰 𝑛 𝜽̈ ∗,𝑛 𝜴 𝑓𝑞 𝑓𝑞 𝑞𝑓 𝑓𝑞 𝑞𝑓 ∗,𝑛

𝑛 ̃ 𝑓 𝑻 𝑛 𝑰 𝑛 − 𝑻 𝑛 (𝑰 𝑛̃ ̇ + [𝜴 𝑻 𝑛𝑞𝑓 𝜴𝑓 ) + 𝑻 𝑛𝑓 𝑞 𝑰 𝑛 (𝑻̃ 𝑓𝑞 𝑓𝑞 𝑞𝑓 𝜴𝑓 )]𝜽 𝑛 ∗,𝑛 𝑛 𝑛̃ 𝑛 ̃ 𝑓 𝑻 𝑛 𝑰 𝑛 (𝑻̃ ̃ + [𝜴 𝑓𝑞 𝑞𝑓 𝜴𝑓 ) − 𝜴𝑓 𝑻 𝑓 𝑞 (𝑰 𝑻 𝑞𝑓 𝜴𝑓 )]𝜽 ∗,𝑛 𝑛 𝑛̃ 𝑛 ̇ 𝑛 ̇ + [𝑻 𝑛𝑓 𝑞 𝑰 𝑛 (𝑻̃ 𝑞𝑓 𝜴𝑓 ) − 𝑻 𝑓 𝑞 (𝑰 𝑻 𝑞𝑓 𝜴𝑓 )]𝜽 𝑃𝑗𝑃𝑛

̃ ≈ 𝑴 𝑛,𝑒𝑥𝑡 + 𝑻 𝑛𝑓 𝑞 (̃𝒆𝑛 + 𝑹 𝑞 𝑓 ̃𝑃 − 𝑻 𝑛𝑓 𝑞 (̃𝒆𝑛 − 𝑹 𝑞

𝑛 𝑃 𝑗+1

𝑗 𝑗 𝑗 𝑗 )𝑻 𝑛𝑗 𝑞 (𝑲 11 𝜸 + 𝑲 12 𝜿 )

in which 𝜇𝐾 and 𝜇𝑀 are stiffness and mass damping coefficients, and [ 𝑛 ] 𝑪 𝑖𝑛𝑒 is computed as: [ ] [ 𝑛 ] 𝟎 𝟎 𝑪 𝑛13 𝑪 𝑛14 𝟎 𝟎 (59) 𝑪 𝑖𝑛𝑒 = 𝑛 𝟎 𝟎 𝟎 𝑪 24 𝟎 𝟎

)𝑻 𝑛𝑗+1 (𝑲 𝑗+1 𝜸 𝑗+1 + 𝑲 𝑗+1 𝜿 𝑗+1 ) 𝑞 11 12

𝑗 𝑗 𝑗 𝑗 − 𝑻 𝑛𝑓 𝑞 𝑻 𝑛𝑗 𝑞 (𝑲 21 𝜸 + 𝑲 22 𝜿 )

+ 𝑻 𝑛𝑓 𝑞 𝑻 𝑛𝑗+1 (𝑲 𝑗+1 𝜸 𝑗+1 + 𝑲 𝑗+1 𝜿 𝑗+1 ) 𝑞 21 22 1 𝑀𝑛

(53)

where

𝑠𝑗+1 ∫𝑠𝑗

𝑚𝑻 𝑓 𝑙 𝑻 𝑙𝑞 𝑑𝑠 is the transformation matrix in which 𝑻 𝑛𝑓 𝑞 = between coordinates 𝑓 and 𝑞 averaged over the 𝑛th segment; 𝑼 𝑛𝑙 = 𝑠𝑗+1 ∫𝑠𝑗

𝑠𝑗+1 𝑚𝑻 𝑓 𝑙 𝑼 𝑙 𝑑𝑠∕ ∫𝑠𝑗

̃ 𝑓𝑻𝑛 𝑪 𝑛13 = 2𝑴 𝑛 𝜴 𝑓𝑞

𝑚𝑻 𝑓 𝑙 𝑑𝑠 is the base translational displacement 𝑠𝑗+1

̃ 𝑓 𝑻 𝑛 𝒆̃ 𝑛 𝑪 𝑛14 = −2𝑴 𝑛 𝜴 𝑓𝑞

𝑠𝑗+1

averaged over the 𝑛th segment; 𝒖∗,𝑛 = ∫𝑠𝑗 𝑚𝑻 𝑓 𝑞 𝒖∗𝑞 𝑑𝑠∕ ∫𝑠𝑗 𝑚𝑻 𝑓 𝑞 𝑑𝑠 𝑞 is the unknown averaged dynamic translational displacement of the 𝑠𝑗+1 𝑠𝑗+1 𝑛th segment, and 𝜽∗,𝑛 = ∫𝑠𝑗 𝑚𝑻 𝑓 𝑞 𝜽∗ 𝑑𝑠∕ ∫𝑠𝑗 𝑚𝑻 𝑓 𝑞 𝑑𝑠 is the unknown averaged dynamic rotational displacement of the 𝑛th segment. Eqs. (52) and (53) are equivalent to Eqs. (40) and (41) with local undeformed coordinate 𝑙 replaced by coordinate 𝑞 such that dynamic displacements are relative to a base configuration. The strain measures in Eqs. (52) and (53) include the effects associated with base deflections. These strain measures can be computed in a discretized form by applying linear interpolation to Eqs. (48) and (49). Eqs. (52) and (53) can be combined with strain measure interpolations and be converted into conventional matrix form: ⎡ 𝒖̈ ∗,𝑛−1 ⎤ ⎡ 𝒖̇ ∗,𝑛−1 𝑞 𝑞 ⎢ ̈ ∗,𝑛−1 ⎥ ⎢ ̇ ∗,𝑛−1 ⎢ 𝜽 ∗,𝑛 ⎥ ⎢ 𝜽 ∗,𝑛 [ ] ⎢ 𝒖̈ ] ⎢ 𝒖̇ ⎥ [ 𝑴 𝑛𝑒𝑞𝑢 ⎢ ̈ 𝑞∗,𝑛 ⎥ + 𝑪 𝑛𝑒𝑞𝑢 ⎢ ̇ 𝑞∗,𝑛 ⎢ 𝜽∗,𝑛+1 ⎥ ⎢ 𝜽∗,𝑛+1 ⎢ 𝒖̈ 𝑞 ⎥ ⎢ 𝒖̇ 𝑞 ⎢ ̈ ∗,𝑛+1 ⎥ ⎢ ̇ ∗,𝑛+1 ⎣ 𝜽 ⎦ ⎣ 𝜽 [ 𝑛,𝑖𝑛𝑒 ] [ 𝑛,𝑒𝑥𝑡 ] [ 𝑛,𝑒𝑙𝑎 ] ̂ ̂ ̂ = 𝑭𝑓 + 𝑭𝑓 + 𝑭𝑓

⎤ ⎡ ⎥ ⎢ ⎥ ⎢ ⎥ [ 𝑛 ]⎢ ⎥ + 𝑲 𝑒𝑞𝑢 ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎦ ⎣

𝒖∗,𝑛−1 𝑞 𝜽∗,𝑛−1 𝒖∗,𝑛 𝑞 𝜽∗,𝑛 𝒖∗,𝑛+1 𝑞 𝜽∗,𝑛+1

𝑛 𝑛 𝑛̃ 𝑛 ̃ 𝑓 𝑻 𝑛 𝑰 𝑛 + 𝑻 𝑛 𝑰 𝑛 (𝑻̃ 𝑪 𝑛24 = 𝜴 𝑓𝑞 𝑓𝑞 𝑞𝑓 𝜴𝑓 ) − 𝑻 𝑓 𝑞 (𝑰 𝑻 𝑞𝑓 𝜴𝑓 )

The internal forces and moments include linear contributions due to dynamic displacements plus non-linear contributions due to base [ ] deflections. The linear contributions are represented by 𝑲 𝑛𝑒𝑙𝑎 . The internal forces and moments associated with the base deflections are generally non-linear and are represented in the right-hand side of Eq. (54) [ 𝑛,𝑒𝑙𝑎 ] [ ]𝑇 by 𝑭̂ = 𝑭 𝑛,𝑒𝑙𝑎 , 𝑴 𝑛,𝑒𝑙𝑎 . Assuming linear elastic constitutive 𝑓

𝑓

𝑓

equations enables computation of 𝑭 𝑛,𝑒𝑙𝑎 and 𝑴 𝑛,𝑒𝑙𝑎 using cross-sectional 𝑓 𝑓 stiffness matrices with base configuration strain measures: 𝑭 𝑛,𝑒𝑙𝑎 = 𝑻 𝑛𝑓 𝑞 𝑻 𝑛𝑗+1 [𝑲 𝑗+1 𝜸 𝑗+1 + 𝑲 𝑗+1 (𝒌𝑗+1 − 𝒌𝑗+1 )] 𝑞 𝑞 𝑓 𝑙 11 𝑞 12

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

𝑗 𝑗 𝑗 𝑗 𝑗 − 𝑻 𝑛𝑓 𝑞 𝑻 𝑛𝑗 𝑞 [𝑲 11 𝜸 𝑞 + 𝑲 12 (𝒌𝑞 − 𝒌𝑙 )]

𝑴 𝑛,𝑒𝑙𝑎 𝑓

=

𝑻 𝑛𝑓 𝑞 (̃𝒆𝑛

𝑗 𝑛 ̃ 𝑃 𝑃 )𝑻 𝑛𝑗 [𝑲 𝑗 𝜸 𝑗 +𝑹 𝑞 𝑞 11 𝑞

𝑃 𝑛 𝑃 𝑗+1

̃ − 𝑻 𝑛𝑓 𝑞 (̃𝒆𝑛 − 𝑹 𝑞

+ 𝑲 𝑗12 (𝒌𝑗𝑞

(60) − 𝒌𝑗𝑙 )]

)𝑻 𝑛𝑗+1 [𝑲 𝑗+1 𝜸 𝑗+1 + 𝑲 𝑗+1 (𝒌𝑗+1 − 𝒌𝑗+1 )] 𝑞 𝑞 𝑙 11 𝑞 12

𝑗 𝑗 𝑗 𝑗 𝑗 − 𝑻 𝑛𝑓 𝑞 𝑻 𝑛𝑗 𝑞 [𝑲 21 𝜸 𝑞 + 𝑲 22 (𝒌𝑞 − 𝒌𝑙 )]

+ 𝑻 𝑛𝑓 𝑞 𝑻 𝑛𝑗+1 [𝑲 𝑗+1 𝜸 𝑗+1 + 𝑲 𝑗+1 (𝒌𝑗+1 − 𝒌𝑗+1 )] 𝑞 𝑞 𝑙 21 𝑞 22

(54)

The inertial properties of the 𝑛th beam segment can be represented as an equivalent mass matrix: [ ] [ ] 𝟎 𝟎 𝑴 𝑛 𝑻 𝑛𝑓 𝑞 −𝑴 𝑛 𝑻 𝑛𝑓 𝑞 𝒆̃ 𝑛 𝟎 𝟎 𝑴 𝑛𝑒𝑞𝑢 = (55) 𝟎 𝟎 𝟎 𝑻 𝑛𝑓 𝑞 𝑰 𝑛 𝟎 𝟎

(61)

The other two terms in the right-hand side of Eq. (54) represent [ 𝑛,𝑒𝑥𝑡 ] [ ]𝑇 two distinct quantities: 𝑭̂ 𝑓 = 𝑭 𝑛,𝑒𝑥𝑡 , 𝑴 𝑛,𝑒𝑥𝑡 represents the ex𝑓 𝑓 [ 𝑛,𝑖𝑛𝑒 ] ternally applied loading on the 𝑛th beam segment, and 𝑭̂ 𝑓 = [ ]𝑇 𝑛,𝑖𝑛𝑒 𝑛,𝑖𝑛𝑒 𝑭 𝑓 , 𝑴𝑓 represents the inertial loading vector of the 𝑛th beam segment. The inertial loading results from the overall rigid-body motion of the beam:

where 𝑴 𝑛 = Diag(𝑀 𝑛 , 𝑀 𝑛 , 𝑀 𝑛 ). The stiffness properties of the 𝑛th beam segment can be represented as an equivalent stiffness matrix: [ ] [ ] [ ] 𝑲 𝑛𝑒𝑞𝑢 = 𝑲 𝑛𝑒𝑙𝑎 + 𝑲 𝑛𝑖𝑛𝑒 (56) [ 𝑛 ] [ 𝑛 ] where 𝑲 𝑒𝑙𝑎 is the elastic stiffness and 𝑲 𝑖𝑛𝑒 is the inertial stiffness representing the inertial effects on stiffness due to overall rigid-body motion. The elastic stiffness matrix can be computed from the crosssectional constitutive matrices and the geometry of the undeformed configuration. The inertial stiffness matrix can be computed as: [ ] [ 𝑛 ] 𝟎 𝟎 𝑫 𝑛13 𝑫 𝑛14 𝟎 𝟎 𝑲 𝑖𝑛𝑒 = (57) 𝑛 𝟎 𝟎 𝟎 𝑫 24 𝟎 𝟎

𝑛 ̃̇ 𝑛 𝑛 𝑛 𝑛 ̃̇ 𝑛 𝑛 𝑭 𝑛,𝑖𝑛𝑒 = −𝑴 𝑛 𝒂𝑊 𝑓 − 𝑴 𝜴𝑓 (𝑹 𝑓 + 𝑻 𝑓 𝑙 𝑼 𝑙 ) − 𝑴 𝜴𝑓 𝑻 𝑓 𝑞 𝒆 𝑓

𝑴 𝑛,𝑖𝑛𝑒 𝑓

=

̃ 𝑓𝜴 ̃ 𝑓 (𝑹𝑛 + 𝑻 𝑛 𝑼 𝑛 ) − 𝑴 𝑛 𝜴 ̃ 𝑓𝜴 ̃ 𝑓 𝑻 𝑛 𝒆𝑛 − 𝑴 𝑛𝜴 𝑓 𝑓𝑙 𝑙 𝑓𝑞

(62)

̃ 𝑓 𝑻 𝑛 𝑰 𝑛 𝑻 𝑛 𝜴𝑓 −𝜴 𝑓𝑞 𝑞𝑓

(63)

− 𝑻 𝑛𝑓 𝑞 𝑰 𝑛 𝑻 𝑛𝑞𝑓 𝜴̇ 𝑓

A matrix formulation for the undeformed beam configuration similar to Eq. (54) can be similarly obtained through Eqs. (55)–(63) by replacing 𝑻 𝑛𝑓 𝑞 /𝑻 𝑛𝑞𝑓 with 𝑻 𝑛𝑓 𝑙 /𝑻 𝑛𝑙𝑓 , replacing (𝑹𝑛𝑓 + 𝑻 𝑛𝑓 𝑙 𝑼 𝑛𝑙 ) with 𝑹𝑛𝑓 , and letting 𝜿 𝑞 and 𝜸 𝑞 equal to zero.

where

4.3. Implementation

̃ 𝑓𝜴 ̃ 𝑓𝑻𝑛 ̃̇ 𝑓 𝑻 𝑛 + 𝑴 𝑛 𝜴 𝑫 𝑛13 = 𝑴 𝑛 𝜴 𝑓𝑞 𝑓𝑞 The piecewise-linear MBBT can be implemented in the time-domain as a linearized updated Lagrangian formulation for large-displacement non-linear beams using separation of displacements. A time-stepping scheme is developed to solve the assembled system of Eq. (54) in which the base configuration for any time-step is selected as the configuration at the end of the prior time-step. Incremental time-steps are sufficiently small such that incremental dynamic displacements can be computed using the linearized MBBT. Quantities that need to be updated for each

̃ 𝑓𝜴 ̃ 𝑓 𝑻 𝑛 𝒆̃ 𝑛 ̃̇ 𝑓 𝑻 𝑛 𝒆̃ 𝑛 − 𝑴 𝑛 𝜴 𝑫 𝑛14 = −𝑴 𝑛 𝜴 𝑓𝑞 𝑓𝑞 𝑛 ̃ 𝑓 𝑻 𝑛 𝑰 𝑛 (𝑻̃ ̃ 𝑓 𝑻 𝑛 (𝑰 𝑛̃ 𝑫 𝑛24 = −𝜴 𝑻 𝑛𝑞𝑓 𝜴𝑓 ) + 𝜴 𝑓𝑞 𝑓𝑞 𝑞𝑓 𝜴𝑓 ) 𝑛 ̇ − 𝑻 𝑛𝑓 𝑞 (𝑰 𝑛̃ 𝑻 𝑛𝑞𝑓 𝜴̇ 𝑓 ) + 𝑻 𝑛𝑓 𝑞 𝑰 𝑛 (𝑻̃ 𝑞𝑓 𝜴𝑓 )

The damping properties of the 𝑛th [segment can also be repre] sented as an equivalent damping matrix, 𝑪 𝑛𝑒𝑞𝑢 . Two types of velocityproportional forces are accounted for in the effective damping terms. 166

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

time-step are: the base translational displacement, 𝑚 𝑼 𝑛𝑙 ; the transformation matrix between floating coordinate and local base coordinates, 𝑚 𝑻 𝑛 ; the transformation matrix between local coordinates along the 𝑓𝑞 𝑚 𝑛𝑗+1 ; the strain associated with base base configuration, 𝑚 𝑻 𝑛𝑗 𝑞 and 𝑻 𝑞 𝑗 𝑚 deflections, 𝜸 𝑞 ; the curvature vector of the base configuration, 𝑚 𝒌𝑗𝑞 , 𝑗

𝑛

𝑛 𝑗+1

and vectors 𝑚 𝑹𝑃𝑞 𝑃 and 𝑚 𝑹𝑃𝑞 𝑃 . The left superscript 𝑚 denotes the quantity evaluated at the 𝑚th time-step; the right superscript indicates locations along the beam, e.g., 𝑚 𝒖∗,𝑛 and 𝑚 𝜽∗,𝑛 are the incremental 𝑞 displacements of the 𝑛th beam segment solved at the 𝑚th time-step based on the configuration updated from the (𝑚 − 1)th time-step. These quantities are updated at the end of the 𝑚th time-step: 𝑚

𝑻 𝑛𝑙𝑞 =

𝑚 𝑚

𝑚 𝑚

𝑻 𝑛𝑙𝑞

𝑚−1

𝑼 𝑛𝑙 +

𝑻 𝑛𝑓 𝑞 = 𝑻 𝑛𝑓 𝑙

𝑚 𝑚

𝑼 𝑛𝑙 =

𝑚−1

𝑻 𝑛𝑗 𝑞 =

𝑻 𝑛𝑗+1 = 𝑞 𝑃𝑗𝑃𝑛

𝑹𝑞

𝑃 𝑛 𝑃 𝑗+1

𝑹𝑞

𝑚 𝑗 𝜸𝑞

𝑚

𝑚

𝑻 𝑛𝑞𝑑 = 𝑚

𝑻 𝑛𝑙𝑞

𝑚−1

𝑻 𝑛𝑙𝑞 (𝑰 +

(64)

)

(65)

𝑻 𝑛𝑙𝑞

(66)

𝑚 𝑗 𝑚 𝑛 𝑻 𝑙𝑞 𝑻 𝑞𝑙 𝑻 𝑛𝑗 𝑙

(67)

𝑚 𝑗+1 𝑚 𝑛 𝑻 𝑙𝑞 𝑻 𝑞𝑙 𝑻 𝑛𝑗+1 𝑙

(68)

=

𝑚

𝑻 𝑛𝑞𝑙

=

𝑚

𝑻 𝑛𝑞𝑙

=

𝑚

𝑻 𝑗𝑞𝑙

=

𝜽

𝑚 ∗,𝑛 𝒖𝑞

𝑃𝑗𝑃𝑛

𝑚 𝑗 𝑼𝑙 ) (𝑹𝑙 + 𝑚 𝑼 𝑛𝑙 − 𝑻 𝑛𝑗 𝑙 𝑛 𝑗+1 𝑚 𝑗+1 𝑼𝑙 ) (𝑹𝑃𝑙 𝑃 − 𝑚 𝑼 𝑛𝑙 + 𝑻 𝑛𝑗+1 𝑙 ) ( 𝑚 𝑼 𝑛 − 𝑚 𝑼 𝑛−1 𝑗 𝑚 𝑗 𝑙 𝑙 ̃ − 𝒆1 𝒆1 + 𝒌𝑙 𝑼 𝑙 + 𝛥𝑠𝑗

(𝑚 𝑚 ̃𝑗 𝒌𝑞

𝑚 ̃ ∗,𝑛

𝑚

𝑗 𝑚

𝑻 𝑗𝑞𝑙 𝒌̃ 𝑙

𝑻 𝑗𝑙𝑞 +

𝑚

𝑻 𝑗𝑞𝑙

𝑻 𝑛𝑙𝑞 −

𝑚 𝑻 𝑛−1 𝑙𝑞 𝛥𝑠𝑗

(69) Fig. 5. Change of root-mean-square error with element size.

(70) (71)

linearized system within each time-step. The resulting incremental dynamic displacements are used to update the base configuration for the next time-step. The relative tolerance of ode45 is set to 10−6 in all examples. The fixed- and free-end boundary conditions are considered at cross-sections 𝑗 = 1 and 𝑗 = 𝑁 + 1, respectively. The fixed end is simulated using an imaginary zeroth beam segment (𝑛=0) at 𝑠 = −𝛥𝑠11 with displacements 𝜽0 = −𝜽1 and 𝒖0𝑙 = −𝒖1𝑙 , such that the fixedend condition is enforced at the first cross-section. The free end is represented by the last cross-section that has zero stiffness: 𝑲 𝑁+1 = 11 𝑁+1 𝑁+1 𝑲 𝑁+1 = 𝑲 = 𝑲 = 𝟎. A constant mass moment of inertia is 12 21 22 selected for each beam segment as the mass moment evaluated for the undeformed segment. Four numerical examples are presented to demonstrate the effectiveness of the proposed MBBT solvers and the underlying theories. The research codes used in these examples are available for academic purposes from Texas A&M University.

) (72)

in which 𝑚 𝑻 𝑗𝑙𝑞 , 𝑚 𝑻 𝑗𝑞𝑙 , and 𝑚 𝑼 𝑗𝑙 are evaluated at the 𝑗th cross-section by linearly interpolating along the 𝑠 axis: [ ] ∗,𝑛−1 ∗,𝑛 ) 𝛥𝑠𝑛−1𝑗 (𝑰 + 𝑚 𝜽̃ ) − 𝛥𝑠𝑗𝑛 (𝑰 + 𝑚 𝜽̃ 𝑚 𝑗 𝑻 𝑙𝑞 = (𝑚 𝑻 𝑗𝑞𝑙 )𝑇 = 𝑚−1 𝑻 𝑗𝑙𝑞 (73) 𝛥𝑠𝑗 ] [ 𝑗𝑛 𝑚 ∗,𝑛−1 ) 𝛥𝑠𝑛−1𝑗 (𝑚 𝒖∗,𝑛 𝑞 ) − 𝛥𝑠 ( 𝒖𝑞 𝑚 𝑗 (74) 𝑼 𝑙 = 𝑚−1 𝑼 𝑗𝑙 + 𝑚−1 𝑻 𝑗𝑙𝑞 𝛥𝑠𝑗 where 𝛥𝑠𝑗 , 𝛥𝑠𝑛−1𝑗 , and 𝛥𝑠𝑗𝑛 are the curvilinear distances: 𝛥𝑠𝑗 = 𝑠𝑛 -𝑠𝑛−1 ; 𝛥𝑠𝑛−1𝑗 = 𝑠𝑗 -𝑠𝑛−1 ; 𝛥𝑠𝑗𝑛 = 𝑠𝑛 -𝑠𝑗 (see Fig. 3). The evolving base configuration is updated from the prior timestep including all non-linear and large-angle effects. The coefficient matrices and loading vectors in Eq. (54) are renewed for the (𝑚+1)th time-step at the end of the 𝑚th time-step using the updated quantities from Eqs. (64)–(74). Initial conditions for Eq. (54) are also renewed at each time-step. Initial displacements are set to zero because they are relative to the new base configuration. Initial velocities are obtained by transferring the velocities computed at the prior step to the newly updated base coordinates: 𝑚+1 ∗,𝑛 𝒖̇ 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 𝑞 𝑚+1 ̇ ∗,𝑛 𝜽 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 𝑞

= =

𝑚 ∗,𝑛 𝒖̇ 𝑞 𝑚 𝑛 𝑚 ̇ ∗,𝑛 𝑻 𝑑𝑞 𝜽𝑞

𝑚

𝑻 𝑛𝑑𝑞

5.1. Convergence and order of accuracy The convergence and order of accuracy of the proposed piecewiselinear method are quantitatively verified. The benchmark problem used to verify the spatial convergence is pure bending of a cantilever beam subject to tip moment. This case is an established benchmark for which a non-linear analytic static solution exists [9]. A straight cantilever beam of length 1 and uniform bending stiffness 2 is subject to static tip moment of 4𝜋. The exact static solution predicts a closed circle with a diameter of 1∕2𝜋. The static MBBT solver is used to compute this benchmark problem with seven progressively finer discretization models. The number of equal-length elements in each these beam models is 5, 10, 20, 50, 100, 200, and 500. Starting from the initial straight configuration, the static solver computes element displacements and updates the configuration until the root mean square of the residual of element displacements is less than 10−6 . All models are found to reach their convergent solutions within 20 iterations. The convergent solution of each model is compared to the exact solution using the root-meansquare error (RMSE). The RMSE is computed by taking the root mean square (RMS) of the errors between the exact translational element displacements predicted by the analytical circle and the numerical results predicted by MBBT. Fig. 5 shows that RMSE decreases with element size with second-order accuracy in spatial discretization. The same benchmark problem is used to show convergence in temporal discretization. The cantilever beam now is provided with physical properties necessary to quantify dynamic effects: the beam is

(75) (76)

An explicit linear system of ordinary differential equations (ODE’s) for the (𝑚+1)th time-step is then established by assembling the renewed Eq. (54) for all beam segments, which can be solved explicitly for the renewed initial conditions using an existing ODE solver. 5. Numerical examples Two finite-volume solvers are developed to solve the system of governing equations. The first solver is a piecewise-linear dynamic solver that computes beam motions based on piecewise-linear MBBT using a continuously evolving base configuration. The second solver is a static simplification of the dynamic solver, in which the mass and damping matrices are set to be zero and only the elastic and forcing terms are updated with solution steps. The dynamic solver uses the MATLAB ‘‘ode45’’ numerical integrator to integrate the locally 167

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

Table 1 Convergence analysis for time-step size. 𝑅𝑀𝑆(𝒖 −𝒖𝛥𝑡∕2 )

𝛥𝑡 (s)

𝑅𝑀𝑆(𝒖𝛥𝑡 − 𝒖5𝑒−6 )

𝑅𝑀𝑆(𝒖𝛥𝑡 −𝒖5𝑒−6 ) 𝑅𝑀𝑆(𝒖5𝑒−6 )

𝑅𝑀𝑆(𝒖𝛥𝑡 − 𝒖𝛥𝑡∕2 )

log2 [ 𝑅𝑀𝑆(𝒖 𝛥𝑡

128e−5 64e−5 32e−5 16e−5 8e−5 4e−5 2e−5 1e−5

7.480e−2 7.948e−3 1.286e−3 1.014e−3 2.988e−4 6.982e−5 3.308e−5 1.359e−5

144% 15.3% 2.48% 1.96% 0.576% 0.135% 0.064% 0.026%

7.913e−2 7.418e−3 1.538e−3 1.277e−3 2.580e−4 5.778e−5 2.009e−5 1.359e−5

3.415 2.270 0.268 2.308 2.159 1.524 0.564 –

𝛥𝑡∕2 −𝒖𝛥𝑡∕4 )

]

Fig. 7. 45◦ -curved cantilever beam.

Table 2 Tip coordinates (𝑥, 𝑦, 𝑧) of the deformed curved beam.

Fig. 6. The cantilever beam subject to tip moment.

0.00331 m high, 0.00331 m wide, 1 m long, and is made of steel with Young’s modulus of 200 GPa, shear modulus of 79.3 GPa, and density of 7850 kg/m3 . The stiffness and mass damping coefficients are set to 10−6 and 0 in the solver, respectively. The beam is discretized into 200 elements, and its deformation from the initial straight configuration at 0.0064 s is computed. The exact dynamic solution is unknown, so the solution computed using time-step size 𝛥𝑡 = 5 × 10−6 s is treated as a reference solution for analyses of ratios of differences between solutions computed with different time-step sizes. Table 1 shows the time-step sizes used and the convergence analyses, in which vector 𝒖𝛥𝑡 includes all translational element displacements along the beam computed using time-step size of 𝛥𝑡. The results demonstrate the piecewise-linear scheme is consistently convergent in this case but large time-step size can lead to unacceptable solutions (𝛥𝑡 = 128×10−5 s). The order of accuracy (the last column of Table 1) is not constant because the specified time-step is for the base configuration update. The ode45 integrator computes an explicit linear ODE system within each timestep, using a six-stage, fifth-order, Runge–Kutta method with its own adaptive internal steps. Fig. 6 shows the static solutions computed using the 10-element and 200-element discretizations, and the dynamic solutions at 0.0064 s computed using the 200-element model with various time-step sizes. Convergence is verified for the finite-volume discretization and the piecewise-linear MBBT method. The finite-volume method is shown to have second-order accuracy and is used to demonstrate that conventional discretization methods are applicable without special treatment; other discretization methods could alternatively be used to improve accuracy or efficiency. The piecewise-linear time-stepping scheme is shown to be numerically stable when the time-steps are sufficiently small for incremental dynamic displacements to be small enough for linearization. Convergence should be tested for a particular problem to select a suitable time-step size. The time-step size needed for piecewiselinear models is expected to be smaller than needed for a non-linear iterative scheme, such as the Newton–Raphson method, though each step is expected to be less computationally intensive as no iteration is used.

Tip load

𝐹 = 300

Present (20 elements) Present (40 elements) Present (60 elements) ANSYS (20 beam elements) ANSYS (40 beam elements) ANSYS (60 beam elements) ANSYS (632 solid elements) Simo and Vu-Quoc [9] Cardona and Geradin [61] Zupan and Saje [16] Romero [62]

58.79, 58.78, 58.78, 58.56, 58.57, 58.57, 58.74, 58.84, 58.64, 58.78, 58.54,

40.23, 40.20, 40.20, 40.45, 40.44, 40.44, 40.22, 40.08, 40.35, 40.16, 40.48,

𝐹 = 600 22.26 22.25 22.25 22.13 22.12 22.12 22.23 22.33 22.14 22.28 22.12

47.15, 47.15, 47.15, 46.91, 46.91, 46.91, 46.93, 47.23, 47.04, 47.15, 46.89,

53.43, 53.47, 53.47, 53.59, 53.60, 53.60, 53.59, 53.37, 53.50, 53.43, 53.60,

15.69 15.69 15.69 15.58 15.57 15.57 15.59 15.79 15.55 15.74 15.56

It includes bending, shear, torsion, and extension deformations of a 3D beam and has been studied in many articles on non-linear beams with large deformations [9,16,61,62]. The 45◦ -curved cantilever has a radius of 100 with a unit square cross-section (Fig. 7). Young’s modulus is 107 , and shear modulus is 0.5 × 107 . A static tip load 𝐹 is applied in the positive 𝑦-direction with two different magnitudes: 𝐹 =300 and 𝐹 =600. The piecewise-linear static solver is used to compute the beam deflection using three different discretization models. The models differ only by the number of elements into which the beam is divided. Table 2 shows the spatial coordinates of the tip computed using various non-linear beam theories and by ANSYS v18.2. Displacements are computed at the mass center of each element in the finite-volume method. Therefore, the curved beam is modeled with its free-end extended by one half-element length and with the static load acting at the center node of that last element, such that the finite-volume results are directly comparable with the FE results at the same tip node. The 2-node Timoshenko beam element used in ANSYS has the same secondorder accuracy as the proposed finite-volume method, and the results from these two methods converge similarly. The piecewise-linear MBBT solutions agree well with equivalent results from other non-linear beam solvers, demonstrating capture of 3D large deformations of an initially-curved beam. 5.3. Rotating straight beam Direct application of MBBT in a non-inertial frame is demonstrated in this example. The piecewise-linear dynamic solver and a non-linear FE code are applied to simulate the response of a straight cantilever beam having an overall rotation of 10 rad/s about the 𝑦-axis at its fixedend (𝛺𝑦 =10 rad/s). The cantilever beam is 0.5 m high, 0.25 m wide, 10 m long, and is made of steel with the same material properties as in Section 5.1. The beam is subject to 5000 kN static tip force in negative 𝑦-direction plus a dynamic load distributed uniformly along the beam: 𝑤𝑒𝑥𝑡 𝑦 (𝑡) = 𝑊 sin(𝜔𝑡) kN/m. The dynamic loading frequency is selected as the natural frequency of the beam estimated using Timoshenko

5.2. Initially-curved beam A 45◦ -curved cantilever beam subject to static tip load is another classical problem that was first presented by Bathe and Bolourchi [60]. 168

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

Fig. 8. The cantilever beam subject to static and dynamic loads and rotation.

Fig. 9. Time histories of blade tip displacements.

beam theory: 𝜔 = 25.6 rad/s. The amplitude of the dynamic load is equivalent to application of five times the magnitude of gravity: 𝑊 = −48.13 kN/m. The beam is discretized into 30 equal-length elements and the time-step size is 0.005 s in both the MBBT and ANSYS solvers. Fig. 8 shows time histories of tip displacements. Good agreement between solvers is observed in the 𝑥- and 𝑦-directions but not in the 𝑧-direction. Results from the MBBT solver show the expected Coriolis effect on response in the 𝑧-direction. This inertial effect caused by the beam rotation about its 𝑦-axis couples velocity in the 𝑥-direction with Coriolis force in the 𝑧-direction, and velocity in the 𝑧-direction with Coriolis force in the 𝑥-direction. In this example, bending in the 𝑦-direction introduces axial motion in the 𝑥-direction due to foreshortening, and the resulting axial motion introduces Coriolis force in the 𝑧-direction. Coriolis force also exists in the 𝑥-direction, but the beam axial stiffness in this case is sufficiently large that the resulting deformation is minimal. The Coriolis effect demonstrated here is captured through the coupled inertial ‘‘damping term’’ in the MBBT formulation: [ 𝑛 ] 𝑛 ̃ 𝑓 𝑻 𝑛 𝒖̇ 𝑛 . This inertial term is proportional to velocity 𝑪 13 𝒖̇ 𝑙 = 2𝑴 𝑛 𝜴 𝑓𝑙 𝑙 such that the Coriolis-excited response in the 𝑧-direction is out of phase with displacements in the other directions (Fig. 8).

Fig. 10. Time histories of blade root bending moments.

Table 3 Comparisons of blade tip displacements.

5.4. Initially-twisted wind turbine blade

Displacement (m)

Edgewise

Flapwise

Axial

𝑅𝑀𝑆𝑀𝐵𝐵𝑇 𝑅𝑀𝑆𝐵𝑒𝑎𝑚𝐷𝑦𝑛 𝑅𝑀𝑆𝑑𝑖𝑓 𝑓 .

3.167e−1 3.130e−1 3.861e−3

5.730e−2 5.553e−2 3.308e−3

1.408e−3 1.388e−3 2.773e−5

1.23%

5.96%

2.00%

𝑅𝑀𝑆𝑑𝑖𝑓 𝑓 . 𝑅𝑀𝑆𝐵𝑒𝑎𝑚𝐷𝑦𝑛

Wind turbine blades are complex beam structures with varying cross-sections and twisted geometry. Blades used on floating offshore wind turbines also experience large irregular overall motions. The proposed beam theory is perfect for application to these complex beam structures in non-inertial frames. The theory is demonstrated on one blade of the National Renewable Energy Laboratory (NREL) 5-MW reference wind turbine [63]. The BeamDyn code is used as a benchmark. BeamDyn is a modular code included in the FAST package made available by the NREL. BeamDyn is based on non-linear FE formulation of GEBT [49,64,65]. The blade is 61.5 m long and is rotating about its root at 12.1 RPM in a vertical plane. Gravity is the only external load applied in this example (𝑔 = 9.80665 m∕s2 ). The blade is discretized into 23 elements in the MBBT solver, and is represented by a single fifth-order Legendre spectral finite element in BeamDyn. The stiffness damping coefficient is set to 0.003 and the time-step is set to 0.005 s in both solvers. Figs. 9 and 10 respectively show the blade tip displacements and root bending moments for the first minute of simulation. Table 3 lists the RMS of displacements of tip translation evaluated for the last 3 min of a 5 min simulation to neglect startup transients. The results from two very different geometrically-exact beam solvers show excellent agreement. BeamDyn is based on an existing GEBT formulation, in which rotational quantities are total rotations referred to an inertial frame. Special treatment is used for discretization of rotational quantities. The rigid-body rotations of the blade and of each blade element are removed from the total rotation. The remaining relative rotation within

each element is then interpolated using Lagrangian-interpolant shape functions. The total rotation of each element is restored afterwards to establish governing equations in an inertial frame [49]. The MBBT implementation is more straightforward. The overall rigid-body rotation is separated explicitly in the MBBT formulation and spatial discretization is applied directly to the linearized MBBT formulation in a floating non-inertial frame. 6. Conclusions A new momentum-based beam theory (MBBT) has been derived in a floating non-inertial frame for arbitrarily-shaped beams having large deformations and overall motions. MBBT is a geometrically-exact non-linear dynamic beam theory that fully addresses geometric and kinematic nonlinearities and inertial coupling effects. Expressing MBBT in an inertial frame has been shown to be mathematically equivalent to an existing formulation of GEBT. Specific advantages of the new theory have been demonstrated through linearization, simplification, and the development of a separation of displacements technique. Separation of displacements enables application of piecewise-linear MBBT to large non-linear dynamic displacements. The effectiveness of the theoretical developments has been demonstrated through numerical implementation. A finite-volume formulation of linearized MBBT has 169

S. Tang and B. Sweetman

International Journal of Non-Linear Mechanics 113 (2019) 158–170

been developed and implemented in two solvers: a piecewise-linear static solver, and a piecewise-linear dynamic solver that computes dynamic motion relative to a series of evolving base configurations in a non-inertial frame. The static MBBT solver has been shown to effectively capture geometric nonlinearities. The dynamic MBBT solver has been demonstrated to effectively capture non-linear geometric and kinematic behaviors and inertial coupling effects using only linear solutions without iteration at each time-step.

[33] J.C. Simo, L. Vu-Quoc, The role of non-linear theories in transient dynamic analysis of flexible structures, J. Sound Vib. 119 (1987) 487–508. [34] H. Kim, H.H. Yoo, J. Chung, Dynamic model for free vibration and response analysis of rotating beams, J. Sound Vib. 332 (2013) 5917–5928. [35] G. Zhao, Z. Wu, Coupling vibration analysis of rotating three-dimensional cantilever beam, Comput. Struct. 179 (2017) 64–74. [36] H.H. Yoo, R.R. Ryan, R.A. Scott, Dynamics of flexible beams undergoing overall motions, J. Sound Vib. 181 (1995) 261–278. [37] J.L. Meek, H. Liu, Nonlinear dynamics analysis of flexible beams under large overall motions and the flexible manipulator simulation, Comput. Struct. 56 (1995) 1–14. [38] J.Y. Liu, J.Z. Hong, Dynamics of three-dimensional beams undergoing large overall motion, Eur. J. Mech. A/Solids 23 (2004) 1051–1068. [39] J.Y. Liu, J.Z. Hong, Geometric stiffening effect on rigid-flexible coupling dynamics of an elastic beam, J. Sound Vib. 278 (2004) 1147–1162. [40] J.C. Simo, L. Vu-Quoc, On the dynamics in space of rods undergoing large motions - A geometrically exact approach, Comput. Methods Appl. Mech. Engrg. 66 (1988) 125–161. [41] D.H. Hodges, Geometrically exact, intrinsic theory for dynamics of curved and twisted anisotropic beams, AIAA J. 41 (2003) 1131–1137. [42] A. Ibrahimbegović, On finite element implementation of geometrically nonlinear Reissner’s beam theory: three-dimensional curved beam elements, Comput. Methods Appl. Mech. Engrg. 122 (1995) 11–26. [43] G. Jelenić, M.A. Crisfield, Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and dynamics, Comput. Methods Appl. Mech. Engrg. 171 (1999) 141–171. [44] E. Zupan, M. Saje, D. Zupan, The quaternion-based three-dimensional beam theory, Comput. Methods Appl. Mech. Engrg. 198 (2009) 3944–3956. [45] V. Sonneville, A. Cardona, O. Brüls, Geometrically exact beam finite element formulated on the special Euclidean group SE (3), Comput. Methods Appl. Mech. Engrg. 268 (2014) 451–474. [46] M.J. Patil, D.H. Hodges, Variable-order finite elements for nonlinear, fully intrinsic beam equations, J. Mech. Mater. Struct. 6 (2011) 479–493. [47] E. Zupan, M. Saje, D. Zupan, Dynamics of spatial beams in quaternion description based on the Newmark integration scheme, Comput. Mech. 51 (2013) 47–64. [48] M.A. Crisfield, G. Jelenić, Objectivity of strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation, in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, The Royal Society, 1999, pp. 1125–1147. [49] Q. Wang, M.A. Sprague, J. Jonkman, N. Johnson, B. Jonkman, Beamdyn: a highfidelity wind turbine blade solver in the FAST modular framework, Wind Energy 20 (2017) 1439–1462. [50] A.A. Shabana, Flexible multibody dynamics: review of past and recent developments, Multibody Syst. Dyn. 1 (1997) 189–222. [51] J.D. Downer, K.C. Park, J.C. Chiou, Dynamics of flexible beams for multibody systems: A computational procedure, Comput. Methods Appl. Mech. Engrg. 96 (1992) 373–408. [52] A.K. Banerjee, S. Nagarajan, Efficient simulation of large overall motion of beams undergoing large deflection, Multibody Syst. Dyn. 1 (1997) 113–126. [53] I. Romero, The interpolation of rotations and its application to finite element models of geometrically exact rods, Comput. Mech. 34 (2004) 121–133. [54] A.A. Shabana, R.Y. Yakoub, Three dimensional absolute nodal coordinate formulation for beam elements: theory, J. Mech. Des. 123 (2001) 606–613. [55] R.Y. Yakoub, A.A. Shabana, Three dimensional absolute nodal coordinate formulation for beam elements: implementation and applications, J. Mech. Des. 123 (2001) 614–621. [56] J. Gerstmayr, M.K. Matikainen, A.M. Mikkola, A geometrically exact beam element based on the absolute nodal coordinate formulation, Multibody Syst. Dyn. 20 (2008) 359–384. [57] J. Gao, B. Sweetman, Design optimization of hull size for spar-based floating offshore wind turbines, J. Ocean Eng. Mar. Energy 4 (2018) 217–229. [58] Y.B. Yang, J.D. Yau, L.J. Leu, Recent developments in geometrically nonlinear and postbuckling analysis of framed structures, Appl. Mech. Rev. 56 (2003) 431–449. [59] B. Sweetman, L. Wang, Momentum cloud method for dynamic simulation of rigid body systems, J. Eng. Mech. 140 (2014) 257–267. [60] K.J. Bathe, S. Bolourchi, Large displacement analysis of three-dimensional beam structures, Internat. J. Numer. Methods Engrg. 14 (1979) 961–986. [61] A. Cardona, M. Geradin, A beam finite element non-linear theory with finite rotations, Internat. J. Numer. Methods Engrg. 26 (1988) 2403–2438. [62] I. Romero, A comparison of finite elements for nonlinear beams: the absolute nodal coordinate and geometrically exact formulations, Multibody Syst. Dyn. 20 (2008) 51–68. [63] J. Jonkman, S. Butterfield, W. Musial, G. Scott, Definition of a 5-MW Reference Wind Turbine for Offshore System Development, Technical Report, National Renewable Energy Laboratory, Golden, CO, USA, 2009. [64] B. Jonkman, J. Jonkman, FAST v8. 16.00 a–bjj, Technical Report, National Renewable Energy Lab, Golden, CO, USA, 2016. [65] Q. Wang, J. Jonkman, M. Sprague, B. Jonkman, BeamDyn User’s Guide and Theory Manual, National Renewable Energy Laboratory, CO, USA, 2016.

References [1] T.R. Kane, R. Ryan, A.K. Banerjee, Dynamics of a cantilever beam attached to a moving base, J. Guid. Control Dyn. 10 (1987) 139–151. [2] I. Sharf, Geometric stiffening in multibody dynamics formulations, J. Guid. Control Dyn. 18 (1995) 882–890. [3] D.H. Hodges, Nonlinear composite beam theory, American Institute of Aeronautics and Astronautics, 2006. [4] O.A. Bauchau, Flexible Multibody Dynamics, Springer, 2011. [5] E. Reissner, On one-dimensional finite-strain beam theory: the plane problem, J. Appl. Math. Phys. 23 (1972) 795–804. [6] E. Reissner, On one-dimensional large-displacement finite-strain beam theory, Stud. Appl. Math. 52 (1973) 87–95. [7] E. Reissner, On finite deformations of space-curved beams, J. Appl. Math. Phys. 32 (1981) 734–744. [8] J.C. Simo, A finite strain beam formulation. The three-dimensional dynamic problem. Part I, Comput. Methods Appl. Mech. Engrg. 49 (1985) 55–70. [9] J.C. Simo, L. Vu-Quoc, A three-dimensional finite-strain rod model. Part II: Computational aspects, Comput. Methods Appl. Mech. Engrg. 58 (1986) 79–116. [10] J.C. Simo, L. Vu-Quoc, A geometrically-exact rod model incorporating shear and torsion-warping deformation, Int. J. Solids Struct. 27 (1991) 371–393. [11] O.A. Bauchau, C.H. Hong, Large displacement analysis of naturally curved and twisted composite beams, AIAA J. 25 (1987) 1469–1475. [12] O.A. Bauchau, C.H. Hong, Nonlinear composite beam theory, J. Appl. Mech. 55 (1988) 156–163. [13] D.A. Danielson, D.H. Hodges, Nonlinear beam kinematics by decomposition of the rotation tensor, J. Appl. Mech. 54 (1987) 258–262. [14] D.A. Danielson, D.H. Hodges, A beam theory for large global rotation, moderate local rotation, and small strain, J. Appl. Mech. 55 (1988) 179–184. [15] D.H. Hodges, A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams, Int. J. Solids Struct. 26 (1990) 1253–1273. [16] D. Zupan, M. Saje, Finite-element formulation of geometrically exact threedimensional beam theories based on interpolation of strain measures, Comput. Methods Appl. Mech. Engrg. 192 (2003) 5209–5248. [17] W. Su, C.E.S. Cesnik, Strain-based geometrically nonlinear beam formulation for modeling very flexible aircraft, Int. J. Solids Struct. 48 (2011) 2349–2360. [18] C.E.S. Cesnik, D.H. Hodges, Variational-asymptotical analysis of initially curved and twisted composite beams, Appl. Mech. Rev. 46 (1993) S211–S220. [19] C.E.S. Cesnik, D.H. Hodges, VABS: A new concept for composite rotor blade cross-sectional modeling, J. Am. Helicopter Soc. 42 (1997) 27–38. [20] W. Yu, D.H. Hodges, V. Volovoi, C.E.S. Cesnik, On Timoshenko-like modeling of initially curved and twisted composite beams, Int. J. Solids Struct. 39 (2002) 5101–5121. [21] W. Yu, V. Volovoi, D.H. Hodges, X. Hong, Validation of the variational asymptotic beam sectional analysis (VABS), AIAA J. 40 (2002) 2015–2113. [22] W. Yu, D.H. Hodges, Elasticity solutions versus asymptotic sectional analysis of homogeneous, isotropic, prismatic beams, J. Appl. Mech. 71 (2004) 15–23. [23] D.H. Hodges, W. Yu, A rigorous, engineer-friendly approach for modelling realistic, composite rotor blades, Wind Energy 10 (2007) 179–193. [24] Y. Jonnalagadda, J.D. Whitcomb, Calculation of effective section stiffness properties for wind turbine blades using homogenization, Wind Energy 17 (2014) 297–316. [25] H. Yang, J. Hong, Z. Yu, Dynamics modelling of a flexible hub–beam system with a tip mass, J. Sound Vib. 266 (2003) 759–774. [26] J.B. Yang, L.J. Jiang, D.C. Chen, Dynamic modelling and control of a rotating Euler–Bernoulli beam, J. Sound Vib. 274 (2004) 863–875. [27] S.K. Sinha, Non-linear dynamic response of a rotating radial Timoshenko beam with periodic pulse loading at the free-end, Int. J. Non-Linear Mech. 40 (2005) 113–149. [28] S.Y. Lee, J.J. Sheu, Free vibration of an extensible rotating inclined Timoshenko beam, J. Sound Vib. 304 (2007) 606–624. [29] C.L. Huang, W.Y. Lin, K.M. Hsiao, Free vibration analysis of rotating Euler beams at high angular velocity, Comput. Struct. 88 (2010) 991–1001. [30] M.H. Tsai, W.Y. Lin, Y.C. Zhou, K.M. Hsiao, Investigation on steady state deformation and free vibration of a rotating inclined Euler beam, Int. J. Mech. Sci. 53 (2011) 1050–1068. [31] Z. Zhao, C. Liu, W. Ma, Characteristics of steady vibration in a rotating hub-beam system, J. Sound Vib. 363 (2016) 571–583. [32] J.C. Simo, L. Vu-Quoc, On the dynamics of flexible beams under large overall motions - The plane case: Part I, J. Appl. Mech. 53 (1986) 849–854. 170