General sliding-beam formulation: A non-material description for analysis of sliding structures and axially moving beams

Journal of Sound and Vibration 480 (2020) 115341

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

General sliding-beam formulation: A non-material description for analysis of sliding structures and axially moving beams Alexander Humer a, ∗ , Ivo Steinbrecher b , Loc Vu-Quoc c a

Institute of Technical Mechanics, Johannes Kepler University, Altenberger Str. 69, 4040, Linz, Austria Institute for Mathematics and Computer-Based Simulation, University of the Bundeswehr Munich, Werner-Heisenberg-Weg 39, 85577, Neubiberg, Germany c Aerospace Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA b

article info

abstract

Article history: Received 27 October 2019 Revised 4 March 2020 Accepted 18 March 2020 Available online 26 March 2020 Handling Editor: Erasmo Carrera

We present a general formulation for problems of sliding structures and axially moving beams that undergo large deformations. The formulation relies on a coordinate transformation that facilitates the analysis of beams characterized by a large sliding motion, for which conventional approaches typically become inefficient. The transformation maps variable domains in the material coordinate, which result, e.g., from the beam’s sliding motion relative to its supports and external loads, onto fixed domains with respect to the new stretched coordinate. We do not only consider supports and loads prescribed at variable points and domains, but their current position relative to the material points of the structure may additionally depend on the current state of deformation. Hamilton’s principle and the geometrically exact theory for shear-deformable beams serve as basis for the derivation of the equations of motion in the stretched coordinate. We introduce a generalized notion of variation that includes the boundaries of variable domains as unknowns and we discuss the implications on the governing equations. Upon a spatial semi-discretization, symmetric mass and tangent stiffness matrices are obtained from the variational formulation of the equations of motion along with non-linear velocity and stiffness-convection terms. Several numerical examples demonstrate both the range of applications and the advantages of the proposed formulation in problems of sliding structures and axially moving beams. © 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

Keywords: Axially moving beams Sliding structures Large deformations Geometrically exact beam theory

1. Introduction In the present paper, we develop a new formulation for the dynamics of highly deformable sliding beams with multiple supports (prismatic or pointwise), which exhibit translational motion in axial direction. The formulation starts from the energy of the system (Hamilton’s principle) leading to symmetric mass and (tangent) stiffness matrices. The classic “spaghetti problem” [1,2] is a particular case of a sliding beam with one prismatic support, motivated by embarrassing real-life experience of having the vibrating tail end of the spaghetti scatter red sauce all over a white shirt when sucking spaghetti into one’s mouth a bit too fast, see Fig. 1a. The rich dynamic behavior of sliding beams is also of great importance in many engineering applications. Early contributions to the sliding-beam problem were motivated by deployable structures in space like spacecraft antennas [3] and tethered satellites [4]. More down-to-earth problems with flexible structures sub-

∗ Corresponding author. E-mail addresses: [email protected] (A. Humer), [email protected] (I. Steinbrecher), [email protected] (L. Vu-Quoc).

https://doi.org/10.1016/j.jsv.2020.115341 0022-460X/© 2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

2

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

Fig. 1. Typical problems of sliding beams: In the so-called “sliding spaghetti problem”, a beam that is retrieved into a guide through a prismatic joint (a). In problems of axially moving continua, we consider structures that undergo a more or less continued sliding motion. The focus lies on the behavior of those parts of a structure, which currently occupy some spatial region (control volume), e.g., the material points of a beam instantaneously located in between two supports (b).

jected to sliding motion include, for instance, belt drives [5–8], drilling [9], band saws [10], rolling mills [11] and brake discs [12,13]. Investigations in these fields of application have coined the term “axially moving continua”, which is frequently used for structures undergoing a continued sliding motion. In problems of axially moving continua, typically not the entire structure, but only the part currently located within some spatial domain is of interest, see Fig. 1b. In fluid dynamics terminology, such non-material (spatial) domain represents a control volume, which the material points of a body flow through. Sliding structures and axially moving continua are characterized by a sliding motion into (one of) their preferred dimensions. To simplify the modeling in problems of this kind, the sliding motion is often decomposed into a translational rigid body motion and a superimposed flexible deformation. As opposed to conventional problems in structural mechanics, support reactions and external loads are not prescribed at fixed material points of the structure under consideration, but their position relative to the structure generally changes in the course of motion. The contributions on sliding structures and axially moving continua have become so numerous that a comprehensive overview is beyond the scope of the present paper. For a detailed exposition, we refer to the review papers by Chen [14] and, more recently, by Marynowski and Kapitaniak [15]. Early investigations on traveling strings and the influence of inertia forces date back to the late 19th century [16,17]. We further mention contributions on the dynamics of axially moving continua by Sack [18], Swope and Ames [19] and Mote [20]. Buffinton and Kane [21] investigated the dynamics of a finite-length beam that is moving over its supports. Wickert and Mote performed classical vibration analysis of axially moving beams [22], which was later extended to non-linear problems, in which the von Kármán theory was adopted [23]. The problem of a beam whose length changes with time was studied by Tabarrok et al. [3], who presented derivations of the equations of motion of an Euler-Bernoulli beam from global balance equations of mass and momentum on the one hand, and based on Hamilton’s principle on the other hand. Hamilton’s principle also was the foundation for the formulation derived by Behdinan et al. [24,25], who discussed the specific nature of a sliding beam as an open system of changing mass, for which McIver’s extensions to Hamilton’s principle [26] need to be accounted for. The seminal paper by Vu-Quoc and Li [2] on the dynamics of sliding beams that are being retrieved into (and rejected from) a prismatic guide provides much of the inspiration for the present paper. Unlike previous contributions, no approximations were introduced in the kinematics of the beam’s deformation, since their formulation is based on the geometrically exact theory for shear-deformable beams [27–30]. Vu-Quoc and Li proposed two equivalent approaches to the sliding-beam problem, i.e., a “full Lagrangian” formulation and an “Eulerian-Lagrangian” formulation. In the full Lagrangian formulation, the governing equations are first formulated in material coordinates before a stretched coordinate is introduced to map the variable material domain onto a fixed domain. The Eulerian-Lagrangian formulation makes use of a spatially fixed intermediate configuration, relative to which the flexible deformation is described. This idea shares similarities with what is referred to as Arbitrary Lagrangian-Eulerian (ALE) formulation in today’s computational mechanics, where suitable non-material control volumes are introduced to facilitate the numerical modeling, e.g., in fluid-structure interaction and material forming processes. Therefore, the term ALE has been adopted in problems of sliding structures and axially moving continua [5,8] as well. Pechstein and Gerstmayr [5] rely on the extension of Lagrange’s equations for open systems with a variable mass present by Irschik and Holl [31]. In a recent contribution, Vetyukov [32] presented a mixed Eulerian-Lagrangian model for axially moving strings and beams, in which a spatial coordinate is used as axial coordinate of the moving structure. Lowe and Cooley [33] presented an Eulerian formulation for translating and rotating continua based on the global balance laws and Reynold’s transport theorem. They referred to their formulation as “Newtonian” in order to distinguish their approach from energy-based formulations. An Eulerian description for the modeling of belt drives was adopted in Refs. [6,34]. Structures subjected to moving loads can be regarded as a dual problem to sliding beams and axially moving continua [2]. Exemplarily, we mention the contributions by Olsson [35], Vu-Quoc and Olsson [36,37] and Cojocaru et al. [38,39] on vehicle-structure interaction. Despite the sliding motion and possibly large deformation, a common assumption in problems of sliding beams and axially moving continua is that the material points momentarily located at supports of the structure are a priori known. Either because the axial deformation at the supports is neglected, or because kinematic constraints are deliberately imposed. Humer [40] extended the original spaghetti problem of Vu-Quoc and Li [2] by also considering that part of the beam which has already been retrieved into a rigid guide, see Fig. 3. The sliding motion of the spaghetti was no longer prescribed kinematically at the “mouth”, i.e., the opening of the guide, which allows the vibrating outside part of the beam and the already retrieved part inside the guide to interact dynamically. As a consequence, the material domains inside and outside the guide are not only variable, but also depend on the current state of deformation. The sliding-beam formulation, which relies on a deformation-dependent coordinate transformation, was proposed in Ref. [40] to map these variable domains onto fixed domains with respect to a stretched

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

3

axial coordinate. The extended sliding spaghetti problem was studied in more detail by Steinbrecher et al. [41]. While slender structures, for which the Euler-Bernoulli hypothesis holds, were considered in Ref. [40], shear-deformation was included in the modeling and a comparison of the sliding-beam formulation to fully Lagrangian models using conventional beam elements and 2-D solid elements was presented. The sliding-beam formulation proved more effective than the Lagrangian approaches, which require carefully tuned contact models to accurately describe the interaction of the structure with the rigid guide. A drawback of the sliding-beam formulation, which is derived from the principle of virtual work, lies in the non-symmetry of the mass matrix and the tangent stiffness matrix. The idea of a support at a variable position that depends on the state of deformation was also adopted in the variable length beam element of Hong and Ren [42] and the extension by Yang et al. [43]. To model sliding joints, they used length-dependent shape-functions and an ad hoc notion of variation in their formulation, which is less general than the proposed approach. The formulation proposed in the present paper is based on the original sliding-beam formulation by Vu-Quoc and Li [2] and generalizations in Refs. [40,41]. We consider the planar motion of shear-deformable beams, which may be subjected to large deformations. Unlike our previous contributions, which focused on the specifics of the sliding spaghetti problem with one [2] or two [40,41] non-material domains, we introduce a coordinate transformation that maps an arbitrary number of variable domains in the material coordinate of the structure onto fixed domains with respect to a new stretched coordinate. We derive the equations of motion from a continuous formulation of Hamilton’s principle, where we introduce a generalized notion of variation as a key ingredient of the proposed formulation. Analytical solutions for large deformation problems of beams are only available for few selected static problems, see, e.g., Refs. [44–46]. For the transient analysis, we therefore adopt the finite element formulation by Simo and Vu-Quoc [29,30], which is applied to the transformed set of the equations of motion formulated in the stretched coordinate. The augmented variational formulation of the equations of motion results in symmetric mass and tangent stiffness matrices upon spatial semi-discretization, by which numerical instabilities can be avoided. The present paper is organized as follows: We first outline the fundamentals of the geometrically exact theory for the planar deformation of shear-deformable beams. We discuss the kinematics of the sliding motion before a deformation-dependent transformation that maps variable domains in the material coordinate onto fixed domains is introduced. Subsequently, the equations of motion are derived from Hamilton’s principle. We do not only take the variation with respect to the displacement field but also include the variable boundaries of the non-material domains. In Section 3, we introduce the spatial semi-discretization of the variational formulation and discuss the symmetry properties of the individual terms that follow from the coordinate transformation. Section 4 presents three numerical examples that show the potential of the proposed formulation in the transient analysis of both sliding structures and axially moving beams. 2. Generalized sliding-beam formulation In the present paper, we focus on the sliding motion of beams whose deformation is restricted to a single plane. Without loss of generality, we assume and coincides with the x1 -axis of { that the } beam’s axis is straight in the undeformed configuration 1 2 a Cartesian inertial frame e1 , e2 , e3 ; the deformation of the beam takes place in the (x , x )-plane of the inertial frame. 2.1. Geometrically exact beam theory For our investigations, we adopt the geometrically exact theory for beams subjected to bending, extensional and shear deformation. The geometrically exact beam theory, which goes back to Reissner [27] and Antman [28], can be regarded as a generalization of Timoshenko’s linear relations [47] to the case of large deformations, which requires appropriate strain measures capable of accounting for large rotations. Later, Simo and Vu-Quoc [29,30] addressed the dynamic response and proposed the well-known finite element formulation for the numerical analysis of geometrically exact beams, which is also used in the present investigations. According to Timoshenko’s kinematic hypothesis, cross-section initially plane and perpendicular to the beam’s axis remain plane but not necessarily perpendicular to the axis in the deformed configuration, {see Fig. } 2 for an illustration. The orientation of the cross-section in the deformed configuration is described by the local frame t 1 , t 2 , which is rotated by the angle 𝜃 about the x3 -axis relative to the inertial frame, t 1 = cos 𝜃 e1 + sin 𝜃 e2 ,

t 2 = − sin 𝜃 e1 + cos 𝜃 e2 .

(1)

Owing to Timoshenko’s kinematic assumption, a beam’s deformed configuration is completely determined from the current position of its axis r0 and the rotation 𝜃 of the local frame. The current position r of a material point initially located at X = X1 e1 + X2 e2 + X3 e3 consequently reads

(

)

(

)

r = r 0 + X 2 t2 + X 3 e3 = X 1 + u1 − X 2 sin 𝜃 e1 + u2 + X 2 cos 𝜃 e2 + X 3 e3 ,

(2)

where the components of the axis’ displacement vector u = r0 − X1 e1 = u1 e1 + u2 e2 have been introduced. The derivative of the axis’ position with respect to the axial material coordinate X1 gives a vector tangent to the deformed axis,

( ) 𝜕 r0 = Λ cos 𝜂 t1 + sin 𝜂 t2 , 1 𝜕X

(3)

whose length represents the axial stretch, i.e., the length ratio of a material line element in the deformed and the undeformed configuration. The tangent to the beam’s axis and the normal to the cross-section enclose the shear angle 𝜂 .

4

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

{

}

Fig. 2. Kinematics of shear-deformable beams: The orientation of cross-sections in the deformed configuration is described by the local frame t 1 , t2 . The normal force N and the shear force Q are introduced in the directions of the local frame.

Fig. 3. Initial and deformed configurations of the (extended) sliding spaghetti problem. The material point currently located at the opening of the guide X1 = l1 (t) changes in the course of motion. In the classic sliding beam problem [2,3], the displacement u1 (X1 = l1 (t), t) at the opening of the guide is prescribed kinematically, which also determines the current length l2 (t) of the outside part that remains to be retrieved. As a consequence, the (already retrieved) part inside the guide and the outside part are decoupled. In the extended problem discussed in Refs. [40,41], no kinematic constraint for the axial motion at the opening is introduced. Therefore, the inside and outside parts, whose material domain depends on the current state of deformation, dynamically interact in the course of the retrieval.

The geometrically exact theory relies on appropriate strain measures that are invariant with respect to arbitrary rigid body motions. Reissner [27] identified the axial force strain 𝜀, the shear strain 𝛾 and the bending strain 𝜅 from variational considerations as

𝜀 = t1 ·

𝜕 r0 − 1 = Λ cos 𝜂 − 1, 𝜕X1

(4)

𝛾 = t2 ·

𝜕 r0 = Λ sin 𝜂, 𝜕X1

(5)

𝜅=

𝜕𝜃 . 𝜕X1

(6)

Not least for notational convenience, we combine the strain measures by introducing a generalized strain vector 𝜞 as

𝜞 =

𝜕 r0 − t1 + 𝜅 e3 , 𝜕X1

(7)

which is related to Reissner’s strain measures (4)–(6) by means of inner products with the basis vectors that span the local frame:

𝜀 = 𝜞 · t1 ,

𝛾 = 𝜞 · t2 ,

𝜅 = 𝜞 · e3 .

(8)

The generalized strain vector 𝜞 is the first of several “generalized vectors” to be introduced in what follows, where the physical dimensions of the individual components are not consistent. As the respective third components always refer to quantities related to bending about the out-of-plane x3 -axis, we sacrifice dimensional consistency in favor of compactness and clarity of our exposition. The components of the cross-sectional stress resultant in the local frame, f = Nt1 + Qt 2 , are referred to as normal force and shear force, respectively, also see Fig. 2. The bending moment M about the x3 -axis is the only non-trivial component of the

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

5

resulting moment with respect to the beam’s axis. In what follows, we employ the conventional linear constitutive equations, which are typically adopted in structural mechanics when physical non-linearities are neglected and the cross-sectional shape allows a decoupling in the strains: N = EA 𝜀,

Q = kGA 𝛾,

M = EI 𝜅.

(9)

In the above relations, EA is referred to as extensional stiffness, kGA as shear stiffness, and EI is the bending stiffness. Generally, these material constants are to be understood as effective stiffnesses. For simple cross-sections, the effective stiffnesses can be related to the elastic moduli from elasticity, where E denotes Young’s modules, G is the shear stiffness and k is the shear correction factor that accounts for the non-uniform distribution of shear stresses over the cross-section. The cross-sectional area and the second moment of inertia about the x3 -axis are denoted by A and I, respectively. In analogy to the vector of generalized strain measures (7), we introduce the vector of stress resultants as F = f + M e3 ,

(10)

and a second-order tensor of (generalized) cross-sectional stiffnesses, D = EA t 1 ⊗ t1 + kGA t 2 ⊗ t 2 + EI e3 ⊗ e3 .

(11)

With these definitions made, we can express the constitutive relations (9) in short as F = D ·𝜞.

(12)

In the geometrically exact beam theory [27,28], the local balance laws of linear momentum and angular momentum are given by

𝜌A ü = 𝜌I 𝜃̈ =

𝜕f + f ext , 𝜕X1 𝜕M + 𝜕X1

(

𝜕 r0 ×f 𝜕X1

(13)

)

· e3 + mext ,

(14)

where f ext and mext represent the external forces and moments per unit length in the reference configuration; 𝜌 denotes the beam’s referential density. Note that material time derivatives are indicated by superimposed dots throughout the present paper. The notion of material time derivatives will subsequently play a crucial role when we introduce a non-material stretched coordinate through a deformation-dependent coordinate transformation. 2.2. Kinematics of sliding motion A key property that distinguishes sliding beams and axially moving continua from conventional problems in structural mechanics is the motion relative to their supports. Kinematic constraints imposed on the deformation by supports are therefore not prescribed at fixed material points of the structure. To illustrate problems addressed by the proposed approach, we exemplarily introduce the (extended) sliding spaghetti problem, which will be generalized subsequently. As spaghetti problem, one refers to a sliding beam that is retrieved into (or ejected from) a guide through a prismatic joint as depicted in Fig. 3. The schematic representation of the guide already reveals an important feature which distinguishes the problem from previous investigations on the classical spaghetti problem [2,3], which only considered that part of the structure which still is to be retrieved (or ejected). To fully decouple the part inside the guide from the outside part, the motion of the beam is required to be completely prescribed at the prismatic joint, i.e., the opening of the guide. In addition to the transverse deflection and the rotation of the cross-section, also the sliding motion in axial direction of the material point instantaneously located at the prismatic joint needs to be constrained kinematically. Following the ideas of our earlier work [40,41], we abandon the strong requirement of a prescribed sliding motion in order to generalize the proposed approach for a larger range of problems. As the inside and outside part then mutually influence each other in general, the entire structure has to be considered in the modeling. For this reason, we subsequently distinguish two domains of the material coordinate  1 and  2 , i.e., the part inside the guide X 1 ∈  1 = [0, l1 (t )] and the outside part X 1 ∈  2 = [l1 (t ), L], respectively. With l2 (t) = L − l1 (t), the initial (t = 0) material length inside the prismatic joint is l1 (0); the initial length of the outside part is denoted by l2 (0), see Fig. 3. Obviously, the domains change in the course of motion, which is also reflected in the kinematic constraints imposed by the support of the structure. Without further specifying the nature of the prismatic joint, let’s assume that the beam’s transverse deflection is prohibited inside the guide. For all times t, the solution has to obey the algebraic equation u2 (X 1 , t ) = 0,

∀X 1 ∈  1 = [0, l1 (t )].

(15)

Due to the sliding motion, the constraint equation is not formulated at fixed material points of the structure. Constraints prescribed at varying material points represent an intrinsic property of sliding structures and axially moving continua and have been referred to as non-material constraints or boundary conditions [40]. Moreover, if the sliding motion is not a priori prescribed kinematically, the material points currently located at the opening and within the guide do not only change over time but also depend on the structure’s current state of deformation. From a finite element perspective, problems of this type involve several complications that may adversely affect both efficiency and accuracy. In a straight-forward approach, the rigid guide of the sliding spaghetti problem could be realized by means

6

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

of a contact formulation. Regardless of the actual implementation, however, contact problems typically imply significant computational expense due to their discontinuous nature. Contact forces not only need to be evaluated, but general purpose contact formulations also require contact search algorithms. In view of the large sliding motion, a comparatively fine discretization is required for all parts of the structures that may get in contact with the guide. Besides computational aspects, the continuity of the contact lines or surfaces is crucial for both efficiency and accuracy in problems involving a large sliding motion. To avoid spurious oscillations resulting from contact forces transitioning inter-element boundaries, higher-order finite elements or surface smoothing techniques are necessary. In Ref. [41], a detailed guidance how to realize the (extended) sliding spaghetti problems within the commercial finite element code ABAQUS [48] by means of conventional Lagrangian beam elements and 2-D solid elements is provided. The challenges of constraints and loadings applied to varying parts of structures have motivated a key idea of the present approach, i.e., the introduction of a simple coordinate transformation that maps variable domains with respect to the material coordinate X1 onto fixed domains. For this purpose, we leave the introductory spaghetti problem behind and continue our derivations from a more general perspective subsequently. 2.3. Coordinate transformation Points of interest that move relative to the structure, regardless of whether kinematic constraints are instantaneously prescribed at these points or moving loads are passing by, partition the structure into one or more domains that comprise a variable set of material points currently located within. Let Xi1 (t ) and Xi1+1 (t ) denote the material points instantaneously located at the boundaries of the i-th domain, we can regard the part of the beam  that is being considered, which not necessarily coincides with the entire physical body , as the union of n subdomains  i

=

n ⋃

 i ⊆ ,

 i = [Xi1 (t ), Xi1+1 (t )].

(16)

i=1

We propose a transformation that maps each variable domain with respect to the axial material coordinate X1 onto a fixed domain with respect to a new coordinate 𝜉 1 . Without loss of generality, we assume that a constant domain in the transformed coordinate coincides with the corresponding material domain in some referential configuration of the considered structure, say, the initial configuration at time t = 0. The simplest possible mapping from X 1 ∈ [Xi1 (t ), Xi1+1 (t )] onto 𝜉 1 ∈ [Xi1 (0), Xi1+1 (0)] is a piecewise linear coordinate transformation 𝜒 ,

𝜉 1 = 𝜒 (X 1 , t ) =

) Xi1+1 (0) − Xi1 (0) ( 1 X − Xi1 (t ) + Xi1 (0), 1 1 Xi+1 (t ) − Xi (t )

X 1 ∈ [Xi1 (t ), Xi1+1 (t )].

(17)

Evaluating the transformation at the boundaries in the material domain, we can immediately convince ourselves that the corresponding boundaries are constant with respect to the new coordinate, which suggests the introduction of

𝜉i1 = 𝜒 (Xi1 (t), t) = Xi1 (0),

𝜉i1+1 = 𝜒 (Xi1+1 (t), t) = Xi1+1 (0).

(18)

We further introduce the length of the i-th domain with respect to the material coordinate and the stretched coordinate as li (t ) = Xi1+1 (t ) − Xi1 (t ),

li (0) = Xi1+1 (0) − Xi1 (0).

(19)

In terms of the domain’s length, the coordinate transformation (17) can be expressed in a more compact form as

𝜉1 =

(

)

li (0) X 1 − Xi1 (t ) + 𝜉i1 , li (t )

X 1 ∈ [Xi1 (t ), Xi1+1 (t )].

(20)

Note that the above transformation is a bijective function, for which the inverse transformation is given by X 1 = 𝜒 −1 ( 𝜉 1 , t ) =

(

)

li (t ) 𝜉 1 − 𝜉i1 + Xi1 (t), li (0)

[

]

𝜉 ∈ 𝜉i1 , 𝜉i1+1 .

(21)

Before we proceed, we illustrate the coordinate transformation by means of three representative examples that are wellestablished in the literature. Sliding spaghetti problem. The problem of a sliding beam that is being retrieved through a prismatic joint has already been introduced above. In the classic problem [2,3], only that part of the beam which still remains to be retrieved is considered in the modeling. We therefore deal with a single variable domain in the material coordinate, for which the left boundary moves relative to the structure. The right boundary coincides with the beam’s outside tip and therefore represents a fixed material point. Assuming that the material point X1 = 0 is initially located at the opening of the guide and the outside tip is X1 = L, i.e., the beam’s initial length is L, the coordinate transformation (20) mapping the variable domain X 1 ∈ [X11 (t ) = 𝓁 (t ), X21 (t ) = L] onto 𝜉 1 ∈ [𝜉11 = 0, 𝜉21 = L] becomes

𝜉1 =

) L ( 1 X − 𝓁 (t ) , l1 (t )

X 1 ∈ [𝓁 (t ), L].

(22)

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

7

The length of the part already retrieved into the guide is denoted by 𝓁 (t) = L − l1 (t) and l1 (t) is the current length of the outside part that remains to be retrieved into the guide. Extended sliding spaghetti problem. In the extended sliding spaghetti problem [40,41], the beam is retracted through a prismatic joint into a straight guide that allows an axial motion but prevents any transversal deflection, see Figs. 1a and 3. Here, the opening of the guide divides the structure into two parts, whose (material) length varies in the course of motion. For this reason, we have two variable domains  1 = [X11 (t ) = 0, X21 (t ) = l1 (t )],  2 = [X21 (t ) = l1 (t ), X31 (t ) = L] in the material coordinate, which are mapped onto constant domains in the stretched coordinate 𝜉 1 ∈ [𝜉11 = 0, 𝜉21 = l1 (0)] and 𝜉 1 ∈ [𝜉21 = l1 (0), 𝜉31 = L] by means of the piecewise linear transformation

⎧ l1 (0) X 1 , ⎪ 𝜉 (X , t) = ⎨ ll1((0t)) ( ) X 1 − l1 (t ) + l1 (0), ⎪2 ⎩ l2 (t ) 1

X 1 ∈ [0, l1 (t )]

1

X 1 ∈ [l1 (t ), L]

.

(23)

Axially moving beam. As a third example, we consider axially moving beams, i.e., structures showing a more or less continued motion in axial direction, where one is typically interested in the deformation relative to the underlying sliding motion, see Fig. 1b. In problems of axially moving beams, we usually study the behavior of those parts of the structure which momentarily occupy some particular spatial region of interest, which can be regarded as a control volume. In the simplest case, we therefore deal with one variable domain that has two variable boundaries, where we adopt the typical assumption of both boundaries showing the same motion, which we denote by s(t) in what follows. The coordinate transformation needs to provide a mapping from X 1 ∈ [X11 (t ) = −s(t ), X21 (t ) = L − s(t )] onto 𝜉 1 ∈ [𝜉11 = 0, 𝜉21 = L], which follows from (17) as

𝜉 1 = X 1 + s(t),

X 1 ∈ [−s(t ), L − s(t )].

(24)

With the same motion prescribed at both boundaries, the same amount of material is contained within the domain of interest. As a consequence, the transformation (24) does not introduce any stretching but merely represents a rigid body translation in our example. The above examples show that the general structure of the coordinate transformation (17) allows us to describe problems with one or more variable domains, in which, e.g., boundary conditions are not prescribed at fixed material points. Introducing f (𝜉 1 , t ), a coordinate transformation 𝜒 , every function f(X1 , t) is converted into a function ̃

̃f (𝜉 1 , t) = ̃f (𝜒 (X 1 , t), t) = f (X 1 , t),

(25)

where a tilde is used to clearly indicate the change in arguments. According to the chain rule of differential calculus, the derivative of f with respect to the material coordinate X1 can be expressed in terms of ̃ f as

𝜕f 𝜕̃f 𝜕𝜉 1 𝜕̃f = = J, 𝜕 X 1 𝜕𝜉 1 𝜕 X 1 𝜕𝜉 1 i

(26)

where the Jacobian of the coordinate transformation Ji = ∂𝜉 1 ∕∂X1 has been introduced. In each variable domain, the Jacobian can be expressed in terms of the ratio between the initial and the current length of the domain with respect to the material coordinate X1 , Ji =

𝜉i1+1 − 𝜉i1 l (0) 𝜕𝜉 1 = i , = 𝜕 X 1 Xi1+1 (t) − Xi1 (t) li (t)

X 1 ∈ [Xi1 (t ), Xi1+1 (t )].

(27)

In case of several domains  i , the coordinate transformation (17) represents a piecewise constant stretching of the material coordinate, i.e., Ji is (spatially) constant within each domain but generally discontinuous at the boundaries Xi1 (t ) and Xi1+1 (t ). f ∕𝜕𝜉 1 is generally Given a continuous derivative ∂f∕∂X1 , the corresponding derivative with respect to the stretched coordinate 𝜕̃ discontinuous at the domain boundaries. The relation among the second spatial derivatives in the i-th domain reads

𝜕2f

(

𝜕 (X 1 )2

=

𝜕 2̃f 𝜕𝜉 1 1 2 𝜕 (𝜉 ) 𝜕 X 1

)2

=

𝜕 2̃f 2 J . 𝜕 (𝜉 1 )2 i

(28)

Clearly, the material time derivative plays a crucial role when introducing a non-material coordinate. Computing the temporal change of some material quantity f expressed in terms of the stretched coordinate 𝜉 1 for a fixed material point X1 , we need to account for the change of 𝜉 1 over time. While the spatial derivative only adds a scalar factor, the material time derivative therefore entails an additional convective term, ḟ =

𝜕 f (X 1 , t) 𝜕̃f (𝜉 1 , t) 𝜕̃f (𝜉 1 , t) 𝜕𝜉 1 = + . 𝜕t 𝜕t 𝜕𝜉 1 𝜕 t

(29)

The time derivative of the stretched coordinate vi = ∂𝜉 1 ∕∂t follows from the transformation (20) as vi =

1 ) l (0) Ẋ 1 (t ) ( ) 𝜕𝜉 1 li (0) Ẋ i+1 (t) ( 1 i 1 i = X ( t ) − X − Xi1+1 (t ) − X 1 , i 2 2 𝜕t li (t ) li (t )

X 1 ∈ [Xi1 (t ), Xi1+1 (t )].

(30)

8

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

Upon substitution of the inverse transformation (21), we can rewrite the above relation in terms of 𝜉 1 as

̃vi =

Ẋ 1i+1 (t ) (

)

𝜉i1 − 𝜉 1 −

li (t )

Ẋ 1i (t ) (

)

[

𝜉i1+1 − 𝜉 1 ,

li (t )

]

𝜉 1 ∈ 𝜉i1 , 𝜉i1+1 ,

(31)

f in 𝜉 1 as which allows us to express the material time derivative of a function ̃

{

𝜕̃f 𝜕̃f 𝜕̃f 𝜕̃f ḟ = + 1̃ vi = + 1 𝜕 t 𝜕𝜉 𝜕 t 𝜕𝜉

Ẋ 1i+1 (

𝜉i1

li

−𝜉

1

)

−

Ẋ 1i (

𝜉i1+1

li

−𝜉

1

)

}

,

[

]

𝜉 1 ∈ 𝜉i1 , 𝜉i1+1 .

(32)

For notational convenience, we have deliberately omitted function arguments in the above relation. The first term represents the local change of ̃ f at a fixed 𝜉 1 , whereas the convective term describes the change of ̃ f due to the motion of the material point instantaneously located at 𝜉 1 . The second material time derivative of ̃ f reads

(

𝜕2f 𝜕 2̃f 𝜕 2̃f 𝜕𝜉 1 𝜕̃f 𝜕 2 𝜉 1 𝜕 2̃f 𝜕𝜉 1 f̈ = 2 = 2 + 2 + 1 + ( )2 1 2 𝜕t 𝜕t 𝜕𝜉 𝜕 t 𝜕 t 𝜕𝜉 𝜕 t 𝜕t 𝜕 𝜉1

)2

,

(33)

where the second time derivative of the transformation (20) in the i-th domain X 1 ∈ [Xi1 (t ), Xi1+1 (t )] gives 1 ) l (0) Ẍ 1 ( ) ) ) li (0) l̇ i Ẋ 1i+1 ( 1 li (0) l̇ i Ẋ 1i ( 1 𝜕 2 𝜉 1 𝜕 vi li (0) Ẍ i+1 ( 1 i 1 1 1 i = = X − X − X − X − 2 Xi − X 1 − 2 Xi+1 − X 1 . 2 2 3 3 i i + 1 2 𝜕t 𝜕t li li li li

(34)

To determine the according relation in terms of the stretched coordinate 𝜉 1 , we can either substitute the inverse transformation (21) in the above relation, or, alternatively, we can make use of (29), by which we obtain

𝜕 2 𝜉 1 𝜕̃vi 𝜕̃vi 𝜕𝜉 1 𝜕̃vi 𝜕̃vi ̃v = + = + 𝜕 t2 𝜕 t 𝜕𝜉 1 𝜕 t 𝜕 t 𝜕𝜉 1 i =

Ẍ 1i+1 (

𝜉i1

li

−𝜉

1

)

Ẍ 1i (

−

𝜉i1+1

li

−𝜉

1

)

) ) Ẋ 1 l̇ ( Ẋ 1 l̇ ( − 2 i+1 i 𝜉i1 − 𝜉 1 + 2 i i 𝜉i1+1 − 𝜉 1 , li li li li

𝜉∈

[

𝜉i1 , 𝜉i1+1

]

(35)

.

In terms of the stretched coordinate, the second material time derivative of ̃ f therefore follows as

𝜕 2̃f 𝜕 2̃f f̈ = 2 + 2 𝜕t 𝜕𝜉 1 𝜕 t 𝜕̃f + 1 𝜕𝜉

{

𝜕 2̃f + 𝜕 (𝜉 1 )2

{

{

)

Ẋ 1i (

−𝜉 −

𝜉i1+1

li

−𝜉

} )

}

) Ẍ 1 ( ) ) Ẋ 1 l̇ ( Ẋ 1 l̇ ( − 𝜉 − 2 i+1 i 𝜉i1 − 𝜉 − i 𝜉i1+1 − 𝜉 + 2 i i 𝜉i1+1 − 𝜉 li li li li li )

𝜉i1

Ẋ 1i+1 (

𝜉i1

li

𝜉i1

li

Ẍ 1i+1 ( li

Ẋ 1i+1 (

)

Ẋ 1i (

−𝜉 −

𝜉i1+1

li

−𝜉

} ) 2

(36)

.

In the above relation, we can pull out the spatial derivative in the last term, which yields the following relation for the second material time derivative of ̃ f in the i-th domain:

𝜕 2̃f 𝜕 2̃f f̈ = 2 + 2 1 𝜕t 𝜕𝜉 𝜕 t +

𝜕̃f 𝜕𝜉 1

{

Ẍ 1i+1 ( li

{

Ẋ 1i+1 (

𝜉i1

li

)

𝜉i1 − 𝜉 1 −

−𝜉

Ẍ 1i ( li

1

)

−

Ẋ 1i (

𝜉i1+1

li

𝜉i1+1 − 𝜉 1

)

−𝜉

}

+

1

)

}

⎛ 𝜕 ⎜ 𝜕̃f 1 𝜕𝜉 ⎜ 𝜕𝜉 1

{

Ẋ 1i+1 ( li

⎝

)

𝜉i1 − 𝜉 1 −

Ẋ 1i ( li

𝜉i1+1 − 𝜉 1

)

}2

⎞ ⎟, ⎟ ⎠

[

]

𝜉 ∈ 𝜉i1 , 𝜉i1+1 , (37)

which boils down to f̈ =

𝜕 2̃f 𝜕 2̃f 𝜕̃f +2 1 ̃ v + 2 𝜕t 𝜕𝜉 𝜕 t i 𝜕𝜉 1

(

)

𝜕̃vi 𝜕̃vi 𝜕 ̃v + 1 − 𝜕 t 𝜕𝜉 1 i 𝜕𝜉

(

𝜕̃f 2 ̃v 𝜕𝜉 1 i

) (38)

vi (31) and its derivatives. The relations for the time derivatives of material quantities (32) and if we make use of the definition of ̃ (37) already reveal that the advantages of the proposed approach do not come free of charge. On the one hand, the algorithmic efforts involved, e.g., in kinematic constraints prescribed at variable material points, are reduced; on the other hand, however, the complexity of the inertia terms of the equations of motion increases significantly. Instead of just a single expression that is linear in the acceleration, we have to deal with several inertia terms that generally depend on the current configuration of the structure.

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

9

2.4. Equations of motion in the stretched coordinate As in conventional problems of structural mechanics, we have several options for deriving the equations of motion in the stretched coordinate 𝜉 1 . With the local balance laws being readily available in the material formulation (13)–(14), we could immediately introduce the stretched coordinate (20) and transform the derivatives of the material quantities accordingly. In their Lagrangian approach, Vu-Quoc and Li [2] follow this procedure to subsequently derive the variational formulation of the equations of motion in the stretched coordinate 𝜉 1 . In Refs. [40,41], we started from d’Alembert’s principle (principle of virtual work), i.e., the variational formulation was first formulated in the material coordinate, before the coordinate transformation was applied to introduce the stretched coordinate. The approach proposed in the present paper aims at computationally efficient analysis of sliding beams and axially moving continua by means of finite elements, which is why we are interested in the variational formulation of the equations of motion. Unlike formulations based on Lagrange’s equations, see, e.g. Refs. [8,43], which anticipate an appropriate discretization and a finite-dimensional set of generalized coordinates, we prefer to remain in a continuous setting to gain a clearer understanding for the physical meaning and all implications of the proposed formulation. For this reason, our considerations start from a continuous form of Hamilton’s principle written in the material representation. In Hamilton’s principle, the variation of the action integral S, which is formed as time integral of the Lagrangian  over the timespan of interest, say, t ∈ [t0 , t1 ], plays a key role. The Lagrangian of conservative systems, in turn, represents a functional that is defined as the difference of the kinetic energy T and the potential energy V, i.e., t1

S=

∫t

 dt,

 = T − V.

(39)

0

In case of (elastically) deformable bodies, the potential energy comprises the work of both internal and external forces, V = W int − W ext .

(40)

The kinetic energy of the structure, or, rather the domain of interest  , is given by T=

( ) 1 1 𝜌A u̇ · u̇ + 𝜌I 𝜃̇ 2 dX 1 = U̇ · I𝜌 · U̇ dX 1 , 2 ∫ 2 ∫

(41)

where the generalized displacement vector U augments the displacement field of the beam’s axis u by the rotation of the crosssection 𝜃 , i.e., U = u + 𝜃 e3 = u1 e1 + u2 e2 + 𝜃 e3 .

(42)

The generalized displacement field is also referred to as pseudo-displacement in the literature [2]. For a compact notation, we have further combined the mass properties of the beam’s cross-sections in a second-order inertia tensor,

(

)

I𝜌 = 𝜌A e1 ⊗ e1 + e2 ⊗ e2 + 𝜌I e3 ⊗ e3 .

(43)

As a first step towards our formulation, we decompose the domain of interest  into n domains  i , which need not be fixed relative to the material points of structure. The kinetic energy therefore follows as the sum over the contributions from each variable domain  i , T=

n ∑

X 1 (t)

Ti ,

i+ 1 1 2 ∫X1 (t)

Ti =

i=1

U̇ · I𝜌 · U̇ dX 1 ,

(44)

i

and so does the potential energy, V=

n ∑

Vi ,

Vi = Wiint − Wiext .

(45)

i=1

For the constitutive behavior assumed in the present paper (12), the work performed by the internal forces is quadratic in the generalized strains, X 1 (t)

Wiint =

i+ 1 1 2 ∫X1 (t) i

(

)

EA 𝜀2 + kGA 𝛾 2 + EI 𝜅 2 dX 1 =

X 1 (t)

i+ 1 1 2 ∫X1 (t) i

𝜞 · D · 𝜞 dX 1 .

(46)

In analogy to the vector of stress resultants (10), we combine the externally applied forces f ext and the external torque mext in the vector of generalized external forces, F ext = f ext + mext e3 ,

(47)

which allows us to write the work of the external loads as Wiext =

X 1 (t) i+ 1

∫X1 (t) i

(

)

f ext · u + mext 𝜃 dX 1 =

X 1 (t) i+ 1

∫X1 (t) i

F ext · U dX 1 .

(48)

10

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

The difference of the kinetic energy and the potential energy gives the following compact statement for the Lagrangian, n ∑

=

i ,

i =

i=1

{ ( 1

X 1 (t) i+ 1

∫X1 (t)

2

i

)

U̇ · I𝜌 · U̇ − 𝜞 · D · 𝜞 + F ext · U

}

dX 1 .

(49)

The action integral S can also be split into contributions from each variable domain Si , which, in turn, are integrals over the timespan of interest [t0 , t1 ], S=

n ∑

t1

Si ,

Si =

i=1

∫t

i dt =

0

∫t

0

{ ( 1

X 1 (t)

t1

)

U̇ · I𝜌 · U̇ − 𝜞 · D · 𝜞 + F ext · U 2

i+ 1

∫X1 (t) i

}

dX 1 dt.

(50)

From a mathematical point of view, the action integral can be regarded as a functional in the unknown fields and their rates, i.e., the (generalized) displacement and velocity fields, S = S(U , U̇ ).

(51)

Hamilton’s principle states that, among all admissible motions, the actual motion of a body is such that the action integral is stationary. Equivalently, we require that the first variation of the functional S vanishes identically,

𝛿 S = 0.

(52)

The variation of the action integral is separately evaluated in each variable domain,

𝛿S =

n ∑

𝛿 Si ,

𝛿 Si = 𝛿

i=1

t1

∫t

i dt =

0

(

t1

∫t

)

𝛿 Ti − 𝛿 Wiint + 𝛿 Wiext dt.

(53)

0

The variation itself is obtained as the directional derivative of the functional with respect to the generalized displacement field and the corresponding velocity field,

𝛿 Si (U , U̇ ) = lim

{

s→0

d S (U + s 𝛿 U , U̇ + s 𝛿 U̇ ) ds i

}

.

(54)

To recover d’Alembert’s principle and the local balance equations in strong form from Hamilton’s principle, further steps typically involve integration by parts on the variation of the generalized velocities followed by a localization of the spatial integrals. In our derivations, however, we leave the beaten path and reformulate Hamilton’s principle in terms of the stretched coordinate 𝜉 1 , which plays the pivotal rule in the present approach. For this purpose, we first consider the kinetic energy, 𝜉1

Ti =

i+ 1 ̇ 1 ̃ · I𝜌 · U ̃̇ J−1 d𝜉 1 , U i 2 ∫𝜉 1 i

(55)

where the inverse of the Jacobian of the coordinate transformation in the i-th domain (27) enters the integral. Accordingly, the work of the internal forces and the work of the external forces are given by 𝜉1

i+ 1 1 ̃ ·D ̃ J −1 d 𝜉 1 , ̃·𝜞 𝜞 i 2 ∫𝜉 1

Wiint =

(56)

i

𝜉i1+1

Wiext =

∫𝜉 1

̃ ̃ J −1 d 𝜉 1 . F ext · U

(57)

i

i

The Lagrangian of the i-th domain therefore reads

i =

𝜉i1+1

{ (

1 ̃̇ ̃ ·D ̃ ̃·𝜞 ̃̇ − 𝜞 U · I𝜌 · U 2

∫𝜉 1 i

)

} ̃ J−1 d𝜉 1 , +̃ F ext · U i

(58)

and, upon integration over time, we obtain the action integral as 𝜉i1+1

t1

Si =

∫t

0

∫𝜉 1 i

{ (

1 ̃̇ ̃ ·D ̃ ̃·𝜞 ̃̇ − 𝜞 U · I𝜌 · U 2

)

} ̃ J−1 d𝜉 1 dt. +̃ F ext · U i

(59)

Comparing (59) and (50), little has changed upon the introduction of the stretched coordinate (17) at first glance. Due to the nature of the transformation, the domain of integration becomes constant with respect to new independent spatial variable 𝜉 1 . Note that we refrain from putting a tilde over the integral quantities, since their values are not affected by the transformation. What does change, however, are two important aspects in the interpretation of the action integral as a functional. First, we no ̃. longer assume Si to be a functional of the material time derivatives (absolute velocities) of the generalized displacement field U

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

11

Instead, the material time derivatives are replaced by the local changes of the displacement for fixed points with respect to the stretched coordinate 𝜉 1 , i.e., the partial derivatives

̃ 𝜕̃ 𝜕U u 𝜕 𝜃̃ = + e . 𝜕t 𝜕t 𝜕t 3

(60)

The second—less intuitive, but all the more significant—difference is related to the length of the variable domains. Without further specifying the nature of variable domains for the moment, we treat the motion of the domain boundaries Xi1 (t ) and Xi1+1 (t ) as independent coordinates that enter the arguments of the functional,

̃∕𝜕 t, X 1 , Ẋ 1 , X 1 , Ẋ 1 ). ̃, 𝜕U Si = Si ( U i i i+1 i+1

(61)

Therefore, the variation of the integral does not only require the derivatives with respect to the axis’ displacement and its temporal change, but it also includes derivatives with respect to Xi1 and Xi1+1 and their rates, respectively,

{

𝛿 Si = lim s→0

d ̃, 𝜕U ̃ ∕𝜕 t + s 𝜕𝛿 U ̃∕𝜕 t, X 1 + s 𝛿 X 1 , Ẋ 1 + s 𝛿 Ẋ 1 , X 1 + s 𝛿 X 1 , Ẋ 1 + s 𝛿 Ẋ 1 ) ̃ + s 𝛿U S (U i i i i i+1 i+1 i+1 i+1 ds i

}

.

(62)

For notational convenience, we split the variation into two parts,

𝛿 S i = 𝛿u S i + 𝛿l S i ,

(63)

̃ with the domain  i kept constant, where the first part represents the variation in the generalized displacement U {

𝛿u Si = lim s→0

d ̃, 𝜕U ̃ ∕𝜕 t + s 𝜕𝛿 U ̃ ∕𝜕 t, X 1 , Ẋ 1 , X 1 , Ẋ 1 ) ̃ + s 𝛿U S (U i i i+1 i+1 ds i

}

.

(64)

The second part represents the contribution from a change of the variable domain  i ,

{

𝛿l Si = lim s→0

d ̃∕𝜕 t, X 1 + s 𝛿 X 1 , Ẋ 1 + s 𝛿 Ẋ 1 , X 1 + s 𝛿 X 1 , Ẋ 1 + s 𝛿 Ẋ 1 ) ̃, 𝜕U S (U i i i i i+1 i+1 i+1 i+1 ds i

}

,

(65)

where the displacement field is fixed. Accordingly, we have introduced subscripts in 𝛿 u and 𝛿 l to clearly distinguish the respective contributions to the total variation. In what follows, we first consider the variation of the action integral with respect to the (generalized) displacement field,

𝛿u S i =

t1

∫t0

(

)

𝛿u Ti − Wiint + Wiext dt.

(66)

Starting with the variation of the kinetic energy, we obtain

𝛿u T i =

𝜉i1+1

∫𝜉 1

̇

̇

̃ · I𝜌 · 𝛿u U ̃ J −1 d 𝜉 1 = U i

i

𝜉i1+1

∫𝜉 1

̇

̃ · I𝜌 · U

i

(

)

̃ 𝜕𝛿 U ̃ 𝜕𝛿 U ̃v J−1 d𝜉 1 , + 𝜕t 𝜕𝜉 1 i i

(67)

where we have introduced the variation of the material time derivative of the displacement field, which follows from (32) as

(

̃̇ = 𝛿u 𝛿u U

)

̃ 𝜕U ̃ ̃ 𝜕𝛿 U ̃ 𝜕U 𝜕𝛿 U ̃v = ̃v . + + 𝜕 t 𝜕𝜉 1 i 𝜕t 𝜕𝜉 1 i

(68)

̃ in the action integral. For this purpose, the variation of the kinetic energy is We want to eliminate the time derivatives of 𝛿 U split into two parts, 𝛿u T i =

𝜉i1+1

∫𝜉 1 i

̇

̃ · I𝜌 · U

𝜉i + 1 ̃ −1 1 ̃ 𝜕𝛿 U ̃̇ · I𝜌 · 𝜕𝛿 U ̃vi J−1 d𝜉 1 , U Ji d 𝜉 + ∫𝜉 1 𝜕t 𝜕𝜉 1 i 1

(69)

i

with integration by parts performed on the former only, i.e., 𝜉i1+1

t1

∫t

0

∫𝜉 1 i

[

𝜉i + 1 ̃ ̃ J −1 d 𝜉 1 ̃̇ · I𝜌 · 𝜕𝛿 U J−1 d𝜉 1 dt = ̃̇ · I𝜌 · 𝛿 U U U i i ∫ 1 𝜕t 𝜉 1

i

]t1 t0

𝜉i1+1

t1

−

∫t

0

∫𝜉 1 i

{

𝜕 ( ̃̇ ) ̃−U ̃ 𝜕 J i J −1 ̃̇ · I𝜌 · 𝛿 U U · I𝜌 · 𝛿 U 𝜕t 𝜕t i

} Ji−1 d𝜉 1 dt. (70)

̃ is assumed to vanish both at the beginning and the end of the timespan The boundary terms vanish since the variation 𝛿 U considered, ̃ (t0 ) = 𝛿 U ̃ ( t 1 ) = 0. 𝛿U

(71)

12

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

Inserting the material time derivative in the stretched coordinate (29),

𝜕 ( ̃̇ ) 𝜕 U = 𝜕t 𝜕t

(

)

̃ 𝜕U ̃ ̃ 𝜕̃v ̃ ̃ 𝜕U 𝜕2U 𝜕2U 𝜕U ̃v = 2 + 1 ̃vi + 1 i , + 𝜕 t 𝜕𝜉 1 i 𝜕t 𝜕𝜉 𝜕 t 𝜕𝜉 𝜕 t

(72)

along with the following elementary relation among the derivatives of the coordinate transformation,

𝜕 Ji −1 𝜕̃vi , J = 𝜕t i 𝜕𝜉 1

(73)

the time integral of the first term in (69) can be rewritten as 𝜉i1+1

t1

∫t

∫𝜉 1

0

i

𝜉1

1 ̃ i+ 1 ̃̇ · I𝜌 · 𝜕𝛿 U J−1 d𝜉 1 dt = U i ∫ 𝜕t t0 ∫𝜉 1 i

t

{

̃ ̃ ̃ 𝜕2U 𝜕2U 𝜕U + 1 ̃ vi + 2 𝜕t 𝜕𝜉 𝜕 t 𝜕𝜉 1

(

)

̃ 𝜕̃vi 𝜕̃vi 𝜕̃vi 𝜕U − 1̃ vi − 𝜕 t 𝜕𝜉 𝜕 t 𝜕𝜉 1

}

̃ J−1 d𝜉 1 dt. · I𝜌 · 𝛿 U i (74)

For the second term in the variation of the kinetic energy (69), we integrate by parts only on the local term leaving the convective term untouched, 𝜉i1+1

∫𝜉 1 i

𝜉1

̃ i+ 1 ̃̇ · I𝜌 · 𝜕𝛿 U ̃vi J−1 d𝜉 1 = U ∫𝜉 1 𝜕𝜉 1 i

(

i

[

̃ 𝜕U ̃ ̃vi J−1 = · I · 𝛿U i 𝜕t 𝜌

]𝜉i1+1

− 𝜉i1

)

̃ 𝜕U ̃ ̃ 𝜕U 𝜕𝛿 U ̃v · I𝜌 · 1 ̃vi Ji−1 d𝜉 1 + 𝜕 t 𝜕𝜉 1 i 𝜕𝜉

𝜉i1+1

(

∫𝜉 1 i

̃ 𝜕̃vi ̃ 𝜕2U 𝜕U ̃v + 𝜕𝜉 1 𝜕 t i 𝜕 t 𝜕𝜉 1

)

̃ J d𝜉 + · I𝜌 · 𝛿 U i −1

1

𝜉i1+1

∫𝜉 1 i

(75)

̃ ̃ 2 −1 1 𝜕U 𝜕𝛿 U ̃v J d𝜉 . ·I · 𝜕𝜉 1 𝜌 𝜕𝜉 1 i i

Combining the results of (74) and (75), the variation of the kinetic energy with respect to the generalized displacement field in the i-th domain becomes t1

∫t

𝛿u Ti dt =

0

[

t1

∫t

0

̃ 𝜕U ̃ ̃vi J−1 · I · 𝛿U i 𝜕t 𝜌 𝜉i1+1

t1

+

∫t0 ∫𝜉 1 i

]𝜉i1+1 𝜉i1

𝜉i1+1

t1

dt −

∫t

0

{

∫𝜉 1 i

̃ ̃ ̃ 𝜕2U 𝜕2U 𝜕U +2 1 ̃ v + 2 𝜕t 𝜕𝜉 𝜕 t i 𝜕𝜉 1

(

𝜕̃vi 𝜕̃vi ̃v − 𝜕 t 𝜕𝜉 1 i

)}

̃ J−1 d𝜉 1 dt · I𝜌 · 𝛿 U i

̃ ̃ 2 −1 1 𝜕U 𝜕𝛿 U ̃v J d𝜉 dt. · I𝜌 · 1 𝜕𝜉 𝜕𝜉 1 i i (76)

Next, we consider the variation of the work of the internal forces (56) with respect to the (generalized) displacement field. As the stiffnesses (11) are specified in the local frame spanned by the beam’s cross-section in the deformed configuration, we need to account for the change of the frame when taking the variation, i.e.,

𝛿ũt1 = 𝛿 𝜃̃̃t2 ,

𝛿ũt 2 = −𝛿 𝜃̃̃t1 .

(77)

The directional derivative of the work of the internal forces therefore comprises two terms

𝛿u Wiint =

𝜉i1+1

∫𝜉 1

(

) ̃ ·D ̃+1𝜞 ̃ · 𝛿u D ̃ J−1 d𝜉 1 , ̃ · 𝛿u 𝜞 ̃·𝜞 𝜞 i

(78)

2

i

where the variation of the vector of generalized strains follows from the definition (7) as

̃= 𝛿u 𝜞

̃ 𝜕𝛿 U ̃̃t2 . J − 𝛿𝜃 𝜕𝜉 1 i

(79)

We obtain the variation of the stiffness tensor with respect to the displacement field as

(

)

(

)

̃ = EA 𝛿𝜃 ̃t 2 ⊗ ̃t 1 + ̃t 1 ⊗ ̃t2 − kGA 𝛿𝜃 ̃t 1 ⊗ ̃t2 + ̃t2 ⊗ ̃t 1 . 𝛿u D

(80)

F , the variation of the work of the internal forces can be written as Recalling the definition (10) of the vector of stress resultants ̃

𝛿u W int =

𝜉i1+1

{

(

)

̃ 𝜕̃r0 ̃ 𝜕𝛿 U ̃ F· − 𝛿 𝜃̃ × f · e3 𝜕𝜉 1 𝜕𝜉 1

∫𝜉 1 i

} d𝜉 1 .

(81)

If we introduce the components of stress-resultants (9) and the variations of the strains (4)–(6), we recover Reissner’s [27] virtual work expression for the internal forces with the coordinate transformation (17) being applied,

𝛿u W int =

𝜉i1+1

∫𝜉 1 i

(

)

̃ 𝛿u 𝜀̃ + Q ̃ 𝛿u ̃ ̃ 𝛿u 𝜅 N 𝛾 +M ̃ Ji−1 d𝜉 1 .

(82)

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

13

The variation of the work performed by the externally applied loads remains to be considered. Under the assumption that F ext = 0, we obtain the variation 𝛿 u Wext from (57) as the loads are independent of the structure’s state of deformation, i.e., 𝛿u ̃

𝛿u W ext =

𝜉i1+1

∫𝜉 1

̃ J −1 d 𝜉 1 . ̃ F ext · 𝛿 U

(83)

i

i

Forming the variation of the action integral with respect to the displacement in the i-th domain (66) using the results (76), (81) and (83), we obtain the relation

𝛿u S i =

t1

∫t

(

𝛿u T i −

∫t

𝜉i1+1

{

∫𝜉 1

0

i

t1

−

+

𝛿u Wiext

)

t1

dt =

0

t1

−

𝛿u Wiint

𝜉i1+1

∫t

̃ ̃ ̃ 𝜕2U 𝜕2U 𝜕U +2 1 ̃ v + 2 𝜕t 𝜕𝜉 𝜕 t i 𝜕𝜉 1

{

(

0

(

[

̃ 𝜕U ̃ ̃vi J−1 · I · 𝛿U i 𝜕t 𝜌

𝜕̃vi 𝜕̃vi ̃v − 𝜕 t 𝜕𝜉 1 i

)}

)

̃ 𝜕̃r0 ̃ 𝜕𝛿 U ̃ J −1 ̃ F· − 𝛿 𝜃̃ × f · e3 − ̃ F ext · 𝛿 U i 𝜕𝜉 1 𝜕𝜉 1

∫t0 ∫𝜉 1 i

]𝜉i1+1 dt 𝜉i1

̃ J−1 d𝜉 1 dt + · I𝜌 · 𝛿 U i

𝜉i1+1

t1

∫t

0

∫𝜉 1 i

̃ ̃ 2 −1 1 𝜕U 𝜕𝛿 U ̃v J d𝜉 dt ·I · 𝜕𝜉 1 𝜌 𝜕𝜉 1 i i

} d𝜉 1 dt. (84)

The above relation—or rather the integrand of the time integral—corresponds to the variational formulation of the governing equations in sliding-beam problems that has been used in our previous contributions [2,40,41]. The local balance equations of momentum and angular momentum follow from the above variational formulation upon localization once we integrate the ̃ ∕𝜕𝜉 1 by parts and require the virtual displacements 𝛿 U ̃ to be arbitrary. The first term on the right-hand side of term with 𝜕𝛿 U (84) is to be evaluated at the boundaries 𝜉i1 and 𝜉i1+1 of the i-th variable domain. Using (31) and inserting the inverse of the Jacobian of the coordinate transformation (27), we obtain the relation

[

̃ 𝜕U ̃ ̃vi J−1 · I · 𝛿U i 𝜕t 𝜌

]𝜉i1+1

[

=− 𝜉i1

̃ 𝜕U ̃ · I · 𝛿U 𝜕t 𝜌

] 𝜉 1 =𝜉i1+1

Ẋ 1i+1 (t ) +

[

̃ 𝜕U ̃ · I · 𝛿U 𝜕t 𝜌

] 𝜉 1 =𝜉i1

Ẋ 1i (t ).

(85)

Upon summation of all variable domains, these boundary terms cancel for neighboring domains  i and  i+1 in the interior of

̃ are sufficiently smooth. At the boundaries 𝜉 1 and 𝜉 1 of the (sub-)structure  under ̃ and 𝛿 U the structure provided that U 1 n+ 1 consideration, however, the boundary terms do not cancel because no neighboring domains exist. Still, the contributions from these terms vanish in the situations we have outlined before. Consider, for example, the case of a boundary 𝜉11 (or 𝜉n1+1 ) that is a fixed material point, i.e., X11 = const (or Xn1+1 = const). The boundary terms (85) must vanish then, since Ẋ 11 = 0 (or Ẋ 1n+1 = 0)

holds. The boundary terms also vanish if the motion of the beam is fully prescribed kinematically at 𝜉11 and 𝜉n1+1 . In this case, ̃ and variations 𝛿 X 1 , 𝛿 X 1 must be zero. which is typical of many problems of axially moving beams, the virtual displacement 𝛿 U i i+1

̃ is assumed to be (spatially) continuous across domain boundaries, the continuity of the local time derivative of While 𝛿 U ̃∕𝜕 t is not immediately clear, e.g., in the case of concentrated forces acting at the variable the generalized displacement 𝜕 U boundaries. To gain insight, we recall the kinematic compatibility condition at singular surfaces of continuum mechanics [49], U [[F]] + [[ẋ ]] ⊗ n = 0,

(86)

which relates the jump of the deformation gradient F to the jump in the material time derivative of some point x at a singular surface that has the normal velocity U and the outward normal n. The compatibility condition is equivalent to

[[ẋ ]] = −U [[F · n]] .

(87)

In the present problem of a one-dimensional continuum, the domain boundaries velocities Ẋ 1i . At some domain boundary, say Xk1 , the above relation translates into

[[

̃ 𝜕U ̃ 𝜕U + 1̃ v 𝜕 t 𝜕𝜉 i

]]

= −Ẋ 1k

[[

Xi1

can be regarded as singular surfaces with

]]

e1 +

̃ 𝜕U J . 𝜕𝜉 1 k

(88)

From the time derivative of the stretched coordinate (31) and the Jacobian of the transformation (27), we obtain the relation

[[

̃ 𝜕U ̃v 𝜕𝜉 1 i

]]

= −Ẋ 1k

[[

̃ 𝜕U J 𝜕𝜉 1 k

]]

.

(89)

Substituting this result in the kinematic jump condition (88), we immediately find that no jumps of the local time derivative ̃∕𝜕 t across domain boundaries occur: 𝜕U

14

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

[[

̃ 𝜕U 𝜕t

]]

[[ ]] = −Ẋ 1k e1 = 0.

(90)

We can therefore conclude that the boundary terms (85) of neighboring domains do cancel upon summation in general. Before we proceed with the variations of the action integral with respect to the boundaries of the variable domains, we also want to anticipate a few thoughts on the symmetry (or anti-symmetry) of the individual terms in the variational formulation, which plays an important role in our reasoning. As we will discuss in more detail below, the notion of symmetry refers to the ̃ and the structure of matrices we obtain from the variational formulation (84) upon a Galerkin projection of the displacement U ̃ onto finite-dimensional (approximation) spaces. If the motion of the domain boundaries is prescribed virtual displacement 𝛿 U ̃ (and ̃ and 𝛿 U kinematically, ̃ vi is some given function in time that renders the inertia terms non-stationary but bilinear in U their derivatives). If the motion of the domain boundaries is considered as additional degrees of freedom, however, we find a ̃ and variable boundaries X 1 (and their respective derivatives) throughout non-linear coupling among the displacement field U i the individual terms in (84). Including variable boundaries Xi1 as unknowns, however, appropriate equations that specify the nature of these boundaries need to be supplemented, see Section 2.5 below. Irrespective of the actual algebraic constraint for a variable boundary Xi1 and irrespective of whether Lagrange multipliers are introduced to enforce these constraints as in the

present paper or whether unknowns Xi1 are eliminated as in the original formulation [40], we cannot expect symmetric mass and tangent stiffness matrices upon a spatial semi-discretization. The key idea of the proposed approach is to recover a symmetric formulation by including variations of the action integral with respect to the boundaries of the variable domains,

𝛿l S i =

t1

∫t

(

)

𝛿l Ti − Wiint + Wiext dt.

(91)

0

Again, we start with the kinetic energy (55) of the i-th domain and take the directional derivative with respect to Xi1 , Xi1+1 and the corresponding rates,

𝛿l T i =

{

𝜉i1+1

̃ ̃̇ · I𝜌 · 𝜕 U 𝛿l̃vi − 1 U ̃̇ · I𝜌 · U ̃̇ J−1 𝛿l Ji U i 𝜕𝜉 1 2

∫𝜉 1 i

} Ji−1 d𝜉 1 .

(92)

Introducing the variations of the Jacobian of the transformation (27) and the time derivative of the stretched coordinate (31),

(

Ji li

𝛿l J i = −

)

we obtain the relation

𝛿l T i =

𝜉i1+1

̇

̃ · I𝜌 · U

∫𝜉 1 i

𝜉i1+1

+

𝛿l̃vi =

𝛿 Xi1+1 − 𝛿 Xi1 ,

̃ 𝜕U 𝜕𝜉 1

{

) 𝛿 Ẋ 1 ( ) ̃ ( ) 𝛿 Ẋ 1i+1 ( 1 v 𝜉i − 𝜉 1 − i 𝜉i1+1 − 𝜉 1 − i 𝛿 Xi1+1 − 𝛿 Xi1 , li

li

(93)

li

}

) 𝛿 Ẋ 1 ( ) ̃ ( ) 𝛿 Ẋ 1i+1 ( 1 v 𝜉i − 𝜉 1 − i 𝜉i1+1 − 𝜉 1 − i 𝛿 Xi1+1 − 𝛿 Xi1 li

li

li

Ji−1 d𝜉 1

𝛿 Xi1+1 − 𝛿 Xi1 −1 1 1 ̃̇ ̃̇ U · I𝜌 · U Ji d𝜉 . 2 li

∫𝜉 1 i

(94)

When forming the action integral, we want to eliminate all time derivatives of the variations. For this purpose, integration by parts is performed on the terms involving 𝛿 Ẋ 1i and 𝛿 Ẋ 1i+1 , which results in t1

∫t

[

𝛿l Ti dt =

𝜉i1+1

∫𝜉 1

0

i

𝜉i1+1

t1

−

∫t

̃ ̃̇ · I𝜌 · 𝜕 U 𝛿l 𝜒̃ J−1 d𝜉 1 U i 𝜕𝜉 1

∫𝜉 1

0

i

]t1 t0

∫t

0

(

𝜉i1+1

t1

−

∫𝜉 1 i

̃ 𝜕̃vi ̃ ̃ 𝜕2U 𝜕2U 𝜕U +2 1 ̃ v + 2 𝜕t 𝜕𝜉 𝜕 t i 𝜕𝜉 1 𝜕 t

t1 𝜉i + 1 ̃ ̃ 𝜕U 𝜕2U 1 · I𝜌 · 1 𝛿l 𝜒̃ Ji−1 d𝜉 1 dt + ∫t ∫𝜉 1 𝜕t 𝜕𝜉 𝜕 t 2 0 1

i

(

)

· I𝜌 ·

̃ 𝜕U 𝛿 𝜒̃ J−1 d𝜉 1 dt 𝜕𝜉 1 l i

̃ ̃ 𝜕U ̃ ̃ 2 𝜕U 𝜕U 𝜕U ̃v ·I · − ·I · 𝜕 t 𝜌 𝜕 t 𝜕𝜉 1 𝜌 𝜕𝜉 1 i

)

𝛿 Xi1+1 − 𝛿 Xi1 li

Ji−1 d𝜉 1 dt, (95)

where the abbreviation

𝛿l 𝜒̃ =

) 𝛿X1 ( ) 𝛿 Xi1+1 ( 1 𝜉i − 𝜉 1 − i 𝜉i1+1 − 𝜉 1 li

li

(96)

has been introduced. As for the displacement field, we assume all variations to vanish at the beginning and the end of the timespan,

𝛿 Xi1 (t0 ) = 𝛿 Xi1 (t1 ) = 𝛿 Xi1+1 (t0 ) = 𝛿 Xi1+1 (t1 ) = 0.

(97)

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

15

We proceed with the work of the internal forces, for which the variations with respect to the boundaries of a variable domain are obtained as 𝜉i1+1

𝛿l Wiint =

(

∫𝜉 1

) ̃ ·D ̃−1𝜞 ̃ ·D ̃ J−1 𝛿l Ji J−1 d𝜉 1 . ̃ · 𝛿l 𝜞 ̃·𝜞 𝜞 i i

(98)

2

i

With the local frame (1) being independent of changes of the domain boundaries Xi1 and Xi1+1 , the corresponding variations of the local basis vectors vanish identically, as does the variation of the tensor of cross-sectional stiffnesses (11),

̃ = 0. 𝛿l D

𝛿l̃t 1 = 𝛿l̃t2 = 0,

(99)

Using the variation of the Jacobian of the coordinate transformation (93), we can express the variation of the generalized strain vector (7) as

̃ ̃ 𝛿 Xi1+1 − 𝛿 Xi1 𝜕U 𝜕U 𝛿 J = − J . 𝜕𝜉 1 l i 𝜕𝜉 1 i li

̃= 𝛿l 𝜞

(100)

Substituting the above relation along with (93) into (98) gives the variations of the work of the internal forces with respect to the domain boundaries as 𝜉i1+1

𝛿l Wiint =

(

∫𝜉 1

̃ ·D ̃· −𝜞

i

̃ 𝜕U 1 ̃ ̃ ̃ J + 𝜞 ·D·𝜞 𝜕𝜉 1 i 2

)

𝛿 Xi1+1 − 𝛿 Xi1 li

Ji−1 d𝜉 1 .

(101)

The work of the external loadings (57) naturally involves the Jacobian of the coordinate transformation (27) when being formulated in the stretched coordinate. Its variation with respect to the variable boundaries of the i-th domain follows from the directional derivative as

𝛿l Wiext

𝜉i1+1

=

∫𝜉 1

̃ext F

̃ ·U

i

𝛿 Xi1+1 − 𝛿 Xi1 li

−1

Ji

d𝜉 + 1

𝜉i1+1

(

∫𝜉 1 i

𝜕̃ F ext 𝜕̃ F ext 𝛿 Xi1 + 1 𝛿 Xi1+1 1 𝜕 Xi 𝜕 Xi+1

)

̃ J −1 d 𝜉 1 . ·U i

(102)

While the first term follows more or less immediately from the variation of the Jacobian of the coordinate transformation, the second integral is less intuitive from a first glance. A dependence of external loads on the boundaries of the variable domains is introduced upon the coordinate transformation F ext (X 1 ) → ̃ F ext (𝜉 1 , Xi1 , Xi1+1 ), unless loads are uniformly distributed over the structure, i.e., F ext = const. The derivatives are taken keeping the stretched coordinate constant, i.e.,

𝜕 F ext (X 1 (𝜉 1 , Xi1 , Xi1+1 )) 𝜕 X 1 𝜕̃ 𝜕̃ F ext F ext 𝜕𝜉 1 𝜕 X 1 𝜕̃ F ext li (0) = = = 𝜕X1 𝜕𝜉 1 𝜕 X 1 𝜕 Xi1 𝜕𝜉 1 li 𝜕 Xi1 𝜕 Xi1

{ 1−

(

1 𝜉 1 − 𝜉i1 li (0)

)}

,

(103)

where (21) and (27) have been used. Combining the above relation and the corresponding result for the directional derivative with respect to Xi1+1 , the variations of the external loads with respect the domain boundaries can be expressed as

( ) 𝛿X1 ( ) 𝛿X1 𝜕̃ F ext 𝜕̃ F ext 𝜕̃ F ext 𝜕̃ F ext 𝜕̃ F ext i+1 1 1 1 1 1 1 i 𝛿 X + 𝛿 X = 𝜉 − 𝜉 − 𝜉 − 𝜉 = − 1 𝛿l 𝜒. ̃ i i+1 i+1 i 1 1 1 1 𝜕𝜉 l 𝜕𝜉 l 𝜕𝜉 𝜕 Xi 𝜕 Xi+1 i i

(104)

In view of their relevance for practical applications, we want to address two important special cases of external forces. First, we consider a load that is uniformly distributed over the i-th non-material domain of the structure  i , i.e., the load is discontinuous at the domain’s boundaries, which may move relative to the material points:

{

F

ext

=

F 0 = const, 0,

X 1 ∈  i = [Xi1 (t ), Xi1+1 (t )] otherwise

.

(105)

For a more compact notation, we can introduce the Heaviside function H, which allows us two express the work performed by the load as W ext =

(

∫

)

̃ J−1 d𝜉 1 , H(𝜉 1 − 𝜉i1 ) − H(𝜉 1 − 𝜉i1+1 ) F 0 · U i

(106)

F ext no longer needs to be restricted to the i-th domain. According to (102) and (104), the variation of the work performed where ̃ by external load follows as

16

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

𝛿l W ext = −

+

(

)

̃ H(𝜉 1 − 𝜉i1 ) − H(𝜉 1 − 𝜉i1+1 ) F 0 · U

∫

𝛿 Xi1+1 − 𝛿 Xi1 li

Ji−1 d𝜉 1

( )( ) 𝛿X1 ̃ i J−1 d𝜉 1 δ(𝜉 1 − 𝜉i1 ) − δ(𝜉 1 − 𝜉i1+1 ) 𝜉 1 − 𝜉i1+1 F 0 · U i

∫

(107)

li

( )( ) 𝛿X1 ̃ i+1 J−1 d𝜉 1 . δ(𝜉 1 − 𝜉i1 ) − δ(𝜉 1 − 𝜉i1+1 ) 𝜉 1 − 𝜉i1 F 0 · U i

∫

li

where the defining property of the Dirac delta function δ and its role as derivative of the Heaviside function have been utilized. The Heaviside function under the first integral can be eliminated by restricting the integral to the variable domain onto which the load is imposed, which allows us to express the variation of the work of the external load as

𝛿l W ext =

𝛿 Xi1+1 − 𝛿 Xi1 li (0)

(

𝜉i1+1

)

̃ d𝜉 1 + F 0 · U ̃ (𝜉 1 ) 𝛿 X 1 − U ̃ (𝜉 1 ) 𝛿 X 1 . F0 · U i i i+1 i+1

∫𝜉 1

(108)

i

Note that the first term on the right-hand side of the above results corresponds to the expression we obtain for a uniformly distributed load on the entire structure. The additional terms, which can be regarded as concentrated forces that act on the boundaries of the variable domain originate from the discontinuities of the load. If a domain boundary Xi1 is fixed with respect

to the material points or if its motion is prescribed kinematically, the corresponding additional term vanishes, since 𝛿 Xi1 = 0 holds in this case. The second important case is that of a concentrated load which is applied at a variable domain boundary, say Xk1 . The work performed by the concentrated load, which we also denote by F 0 , can be expressed in terms of the Dirac delta function as W ext = F 0 · U (Xk1 ) =

∫

δ(X 1 − Xk1 ) F 0 · U dX 1 .

(109)

According to (102), taking the variations with respect to the variable domain boundaries gives

𝛿l W

X1 ext

=

k

∫X1

k− 1

X1

+

𝜕 δ(X 1 − Xk1 ) 𝜕X1

k+ 1

∫X 1 k

(

𝜕 δ(X 1 − Xk1 ) 𝜕X1

)

X 𝛿 Xk1 − 𝛿 Xk1−1 1 k 𝜕X1 𝜕X1 1 1 1 1 1 𝛿 X + 𝛿 X F · U dX + δ( X − X ) F · U dX 0 0 k ∫X 1 lk−1 𝜕 Xk1−1 k−1 𝜕 Xk1 k k− 1

(

1

)

X1 𝛿 Xk1+1 − 𝛿 Xk1 1 k+ 1 𝜕X1 1 𝜕X1 1 1 1 1 𝛿 X + 𝛿 X F · U dX + δ( X − X ) F · U dX . 0 0 k ∫X 1 lk 𝜕 Xk1 k 𝜕 Xk1+1 k+1

(110)

k

Clearly, only those domains adjacent to the boundary at which the load is applied contribute to the variation. Note that we remain in the material domain for the sake of simplicity, but still, we have to include the dependence of X1 on li for constant 𝜉 1 , i.e., 𝛿l dX 1 = (𝛿 Xi1+1 − 𝛿 Xi1 )∕li dX 1 , which is reflected in the second integrals of each domain. These terms, however, cancel upon integration by parts on the first integrals. Substituting the (inverse) coordinate transformation (20)–(21), we obtain the relation

𝛿l W ext = − 𝛿 Xk1 F 0 ·

X1

k+ 1

∫X 1

δ(X 1 − Xk1 )

k− 1

𝜕U dX 1 . 𝜕X1

(111)

Within Reissner’s beam theory, concentrated loads result in (generalized) displacement fields U that, in general, are not continuously differentiable at the those points at which these loads are applied. Therefore, the integral yields the arithmetic mean of the one-sided limits of U rather than its derivative, i.e.,

(

𝛿l W

ext

1 = − 𝛿 Xk1 F 0 · 2

𝜕 U || 𝜕 U || + | 1 1 𝜕 X |X1 =(X )− 𝜕 X 1 ||X1 =(X1 )+ k

k

)

⎛

=−

⎝

⎞

|

̃ 1 𝜕U 𝛿 Xk1 F 0 · ⎜ Jk−1 1 || ⎜ 2 𝜕𝜉 ||

𝜉 1 =(𝜉k1 )−

+ Jk

̃ || 𝜕U ⎟, | 𝜕𝜉 1 ||𝜉 1 =(𝜉 1 )+ ⎟ k

(112)

⎠

see, e.g., Refs. [50,51]. Combining the results for the kinetic energy (95), the work of the internal forces (101) and the work of the external loadings (102), we eventually obtain the variations of the action integral with respect to the boundaries of the i-th variable domain (91) from the respective time integrals,

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

𝛿l S i = −

𝜉i1+1

t1

∫t

∫𝜉 1

0

i

𝜉i1+1

t1

+

1 2

∫t0 ∫𝜉 1 i

∫t

i

𝜉i1+1

t1

+

∫t

̃ ̃ F ext · U

∫𝜉 1

0

)

t1 𝜉i + 1 ̃ ̃ ̃ 𝜕U 𝜕U 𝜕2U 𝛿 𝜒̃ J−1 d𝜉 1 dt − · I𝜌 · 1 𝛿l 𝜒̃ Ji−1 d𝜉 1 dt ∫t ∫𝜉 1 𝜕𝜉 1 l i 𝜕 t 𝜕𝜉 𝜕t 0 1

· I𝜌 ·

i

̃ ̃ 𝜕U ̃ ̃ 2 𝜕U 𝜕U 𝜕U ̃v ·I · − ·I · 𝜕 t 𝜌 𝜕 t 𝜕𝜉 1 𝜌 𝜕𝜉 1 i 𝜕𝜉

∫𝜉 1

0

(

̃ 𝜕̃vi ̃ ̃ 𝜕2U 𝜕2U 𝜕U +2 1 ̃ v + 𝜕 t2 𝜕𝜉 𝜕 t i 𝜕𝜉 1 𝜕 t

)

𝛿 Xi1+1

− li

𝛿 Xi1 −1 1 Ji d𝜉 dt

( ) 1 𝛿 Xi+1 − 𝛿 Xi1 −1 1 ̃ ̃ ·D ̃ ·D ̃ ̃ · 𝜕 U Ji + 1 𝜞 ̃·𝜞 −𝜞 Ji d𝜉 dt 1

𝜉i1+1

t1

−

(

17

2

𝛿 Xi1+1 − 𝛿 Xi1 li

i

li

Ji−1 d𝜉 1 dt −

𝜉i1+1

t1

∫t

∫𝜉 1

0

i

𝜕̃ F ext ̃ · U 𝛿l 𝜒̃ Ji−1 d𝜉 1 dt. 𝜕𝜉 1 (113)

The sum of the variations with respect to both the (generalized) displacement field and the variable boundaries gives the total variation of the action integral 𝛿 Si in the i-th domain of the structure. After collecting the respective terms related to inertia forces, internal forces and external forces, we arrive at the expression

𝛿 Si =

t1

∫t

[

̃ 𝜕U ̃ ̃vi J−1 · I · 𝛿U i 𝜕t 𝜌

0

𝜉i1+1

t1

+

−

∫t

−

i

𝜉i1+1

−

+

{

𝜉i1+1

∫𝜉 1

0

)

( ) ̃ ̃ + 𝜕 U 𝛿l 𝜒̃ J−1 d𝜉 1 dt · I𝜌 · 𝛿 U i 1

𝜕𝜉

(

1

t1

̃ 𝜕U ̃ 𝛿 Xi1+1 − 𝛿 Xi1 𝜕𝛿 U − 𝜕𝜉 1 𝜕𝜉 1 li

i

)

̃ Ji Ji − 𝛿 𝜃

(

̃ ̃ 𝜕U ̃ ̃ 2 𝜕U 𝜕U 𝜕U ̃v ·I · − ·I · 𝜕 t 𝜌 𝜕 t 𝜕𝜉 1 𝜌 𝜕𝜉 1 i )

𝜕̃r0 ̃ × f · e3 𝜕𝜉 1

)

𝛿 Xi1+1 − 𝛿 Xi1 li

Ji−1 d𝜉 1 dt

} Ji−1 d𝜉 1 dt

1 1 1 ̃ ̃ 𝛿 Xi+1 − 𝛿 Xi −1 1 F·𝜞 Ji d𝜉 dt 2 li

∫t0 ∫𝜉 1 i ∫t

(

̃ F·

𝜉i1+1

t1

i

𝜉i + 1 ̃ 𝜕U 𝜕 U 1 · I𝜌 · 1 𝛿l 𝜒̃ Ji−1 d𝜉 1 dt + ∫t ∫𝜉 1 𝜕t 𝜕𝜉 𝜕 t 2 0

i

t1

∫𝜉 1

0

2̃

∫t ∫𝜉 1 0

∫t

̃ 𝜕̃vi ̃ ̃ 𝜕2U 𝜕2U 𝜕U +2 1 ̃ vi + 2 𝜕t 𝜕𝜉 𝜕 t 𝜕𝜉 1 𝜕 t

i

𝜉i1+1

t1

𝜉i1

(

1

∫𝜉 1

0

𝜉i1+1

t1

dt −

t1 𝜉i + 1 ̃ ̃ ̃ 2 −1 1 𝜕U 𝜕U 𝜕𝛿 U ̃ 𝜕̃vi ̃vi J−1 d𝜉 1 dt + ̃v J d𝜉 dt · I · 𝛿 U ·I · 𝜌 ∫t0 ∫𝜉 1 𝜕𝜉 1 𝜕𝜉 1 i 𝜕𝜉 1 𝜌 𝜕𝜉 1 i i

∫t0 ∫𝜉 1 i t1

]𝜉i1+1

(

̃ext F

·

̃+U ̃ 𝛿U

𝛿 Xi1+1 − 𝛿 Xi1

)

li

i

Ji−1 d𝜉 1 dt −

𝜉i1+1

t1

∫t

0

∫𝜉 1 i

𝜕̃ F ext ̃ · U 𝛿l 𝜒̃ Ji−1 d𝜉 1 dt. 𝜕𝜉 1 (114)

The variational formulation of the equations of motion, which serves as basis for the subsequent finite-element discretization, follows from Hamilton’s principle once we require the integrand of the time integral to vanish. Summing up the contributions from each variable domain, we then obtain the relation n ∑

{

−

i=1

𝜉i1+1

(

∫𝜉 1 i

+

𝜉i1+1

∫𝜉 1 i

−

∫𝜉 1 𝜉i1+1

∫𝜉 1 i

+

𝜉i1+1

∫𝜉 1 i

)

( ) ̃ ̃ + 𝜕 U 𝛿l 𝜒̃ J−1 d𝜉 1 · I𝜌 · 𝛿 U i 1

𝜕𝜉

𝜉i + 1 ̃ ̃ ̃ 2 −1 1 𝜕U 𝜕U 𝜕𝛿 U ̃ 𝜕̃vi ̃vi J−1 d𝜉 1 + ̃v J d𝜉 · I𝜌 · 𝛿 U ·I · i 1 1 ∫𝜉 1 𝜕𝜉 𝜕𝜉 𝜕𝜉 1 𝜌 𝜕𝜉 1 i i 1

i

𝜉i1+1 i

−

̃ 𝜕̃vi ̃ ̃ 𝜕2U 𝜕2U 𝜕U +2 1 ̃ v + 𝜕 t2 𝜕𝜉 𝜕 t i 𝜕𝜉 1 𝜕 t

𝜉i + 1 ̃ ̃ 𝜕U 𝜕2U 1 · I𝜌 · 1 𝛿l 𝜒̃ Ji−1 d𝜉 1 + ∫ 𝜕t 𝜕𝜉 𝜕 t 2 𝜉1 1

{

(

̃ F·

i

̃ 𝜕U ̃ 𝜕𝛿 U − 𝜕𝜉 1 𝜕𝜉 1

(

̃+U ̃ ̃ F ext · 𝛿 U

𝛿 Xi1+1

−

𝛿 Xi1

li

𝛿 Xi1+1 − 𝛿 Xi1 li

(

̃ ̃ 𝜕U ̃ ̃ 2 𝜕U 𝜕U 𝜕U · I𝜌 · − 1 · I𝜌 · 1 ̃ v 𝜕t 𝜕 t 𝜕𝜉 𝜕𝜉 i

)

̃ Ji Ji − 𝛿 𝜃

) Ji−1 d𝜉 1 −

(

𝜉i1+1

∫𝜉 1 i

)

𝜕̃r0 ̃ × f · e3 𝜕𝜉 1 𝜕̃ F ext 𝜕𝜉 1

)

𝛿 Xi1+1 − 𝛿 Xi1 li

} Ji−1 d𝜉 1 −

}

̃ 𝛿l 𝜒̃ J−1 d𝜉 1 ·U i

𝜉i1+1

∫𝜉 1 i

Ji−1 d𝜉 1

1 1 1 ̃ ̃ 𝛿 Xi+1 − 𝛿 Xi −1 1 F·𝜞 Ji d𝜉 2 li

= 0. (115)

18

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

In the following section, we will show that the proposed approach provides the intended symmetry properties of the weak formulation of the governing equation (115) as compared to the original sliding-beam formulation [40,41]. On the downside, however, the variation with respect to variable domain boundaries 𝛿 Xi1 adds a significant number of terms, whose physical meaning has not been addressed so far. First of all, we note that there has been no need so far to further specify the nature of the variable domains in the derivation of the augmented variational formulation of the governing equation. In fact, we will discuss the constraint equations for the variable domain boundaries in some physically meaningful and practical way in what follows. For the moment, however, we can regard the coordinate transformation (20), which introduces the variable domain boundaries Xi1 as additional unknowns, as a mathematical exercise. From such mathematical point of view, we expect the (local) balance laws that govern the dynamics of the structure to remain preserved upon the coordinate transformation. And indeed, the variation with respect to the (generalized) displacement (84) gives the local balance laws of linear momentum (13) and angular momentum (14) also in the stretched coordinate 𝜉 1 . For this purpose, integration by parts eliminates the spatial derivatives of the (generalized) virtual displacements in (84) such that integrals can be localized following standard arguments as in the original formulation in the material coordinate X1 . With the dynamics of sliding beams being fully described by the local balance laws, we face the question what the additional terms introduced by including variations with respect to the boundaries of the variable domains in Hamilton’s principle imply. In the Appendix to this paper, we show in detail that the additional terms do not modify the fundamental balance relations. As a matter of fact, these terms cancel exactly once we re-introduce the material coordinate X1 , which, from a mathematical point of view, is a result that does not come unexpected. As soon as we approximate the continuous relations by means of a finite-dimensional spatial discretization, however, the additional terms no longer vanish identically. Hence, we cannot expect identical results from the original formulation presented in Refs. [40,41] and the approach proposed in the present paper, though both formulations should converge to the exact solution and therefore also to one another. 2.5. Constraint equations for domain boundaries In view of the extensive derivations presented in the previous sections, one might easily overlook the fact that the variational formulation of the equations of motion (115) has not been fully completed yet. The nature of the variable domains, i.e., the motion of each domain’s boundaries relative to the material points of the structure, still remains to be further specified. In what follows, we discuss typical cases of variable domain boundaries and the corresponding algebraic relations that enter the system of equations. As a first example, we return to the extended sliding spaghetti problem (see Fig. 3) that has motivated the generalization of the original sliding-beam formulation in the first place. We recall that a beam is being retrieved (“sucked”) into a guide through an opening, which can be regarded as a prismatic joint. At each moment in time, we find two variable domains: the portion of the beam inside the guide,  1 = [X11 (t ) = 0, X21 (t ) = l1 (t )], and the outside part of the beam,  2 = [X21 (t ) = l1 (t ), X31 (t ) = L]. The prismatic joint, whose location is assumed to be spatially fixed, prohibits any deflection of those points (and possibly a rotation of cross-sections) which are momentarily located at the opening, whereas it does not influence the relative axial motion of the beam and guide. In other words, the position of the point currently located at the opening and the opening itself coincide, i.e., r0 (X 1 = X21 (t ), t ) = popening = const

⇔

r0 (X 1 = X21 (t ), t ) − popening = 0.

(116)

Adopting the sliding-beam formulation, we introduce the coordinate transformation (20) such that the opening of the guide is always located at a fixed point with respect to the stretched coordinate 𝜉 1 = 𝜉21 = l1 (0), i.e., the following condition must hold:

̃r0 (𝜉 1 = 𝜉21 , t) − popening = l1 (t) e1 + ̃ u(𝜉 1 = 𝜉21 , t ) = 0.

(117)

Let the opening of the guide coincide with the origin of the global frame, popening = 0, then we obtain the following scalar-valued algebraic constraint equations C1 , C2 from the above relation: u1 (𝜉 1 = 𝜉21 , t ) = 0, C1 = l1 (t ) + ̃

C2 = ̃ u2 (𝜉 1 = 𝜉21 , t ) = 0.

(118)

The second equation requires the transverse deflection at the opening of guide to vanish, which is similar to conventional beam formulations except that the deflection is prohibited not at a fixed material point but at a fixed point in the stretched coordiu1 , however, is not present in a conventional Lagrangian nate 𝜉 1 = l1 (0). The first equation involving the axial displacement ̃ approach. Indeed, the algebraic equation establishes an explicit relationship between the motion of the variable domain boundary X21 (t ) = l1 (t ) and the structure’s current state of deformation, X21 (t ) = −̃ u1 (𝜉 1 = 𝜉21 , t ),

(119)

i.e., the coordinate transformation is deformation-dependent. Adding the relation C1 = 0 (along with the constraint condition for the transverse deflection within the guide), completes the equations of motion for the extended sliding spaghetti problem. We can further generalize the situation described above by requiring that the spatial position p, at which a domain boundary X 1 (t ) is supposed to be located, also depends on the current state of deformation, i.e., the displacement vector of the material k

u(𝜉k1 , t ). For this situation, we obtain a generic formulation of constraint equations C that may occur point currently located there ̃ in the sliding-beam formulation as u(𝜉k1 , t ), Xk1 (t ), t ) = Xk1 (t ) e1 + ̃ u(𝜉k1 , t ) − p(̃ u(𝜉k1 , t ), t ) = 0. C (̃

(120)

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

19

If the motion of the domain boundary is prescribed kinematically as, e.g., in the original sliding-beam problem [2,3] or typical problems of axially moving continua, the generic relation (120) no longer contains the current state of deformation and consequently reduces to Xk1 (t ) − p(t ) = 0.

(121)

The trivial case of a domain boundary being a fixed material point can also be regarded as a kinematically prescribed motion, in which the boundary remains at the original material point throughout the course of motion, Xk1 (t ) = Xk1 (0).

(122)

If the motion of a domain boundary Xk1 does not depend on the current state of deformation but is fully described kinemat-

ically, the corresponding variation vanishes, i.e., 𝛿 Xk1 = 0. Rather than augmenting the system of equations by an additional

algebraic relation, Xk1 can easily be eliminated from the remaining equations by substituting Xk1 = p(t ) and omitting the equation

associated with 𝛿 Xk1 . In the numerical examples of the present paper, we use the method of Lagrange multipliers to enforce the constraint equations that determine the motion of variable boundaries. The work performed by constraint forces is given by the inner product of the constraint equation (120) and the vector of Lagrange multipliers 𝝀, Wconstr = C (̃ u(𝜉k1 , t ), Xk1 (t ), t ) · 𝝀.

(123)

Taking the variations in the displacements, the Lagrange multipliers and the boundaries of the variable domains, we obtain the relation

(

𝛿 Wconstr = 𝛿u Wconstr + 𝛿l Wconstr + 𝛿𝜆 Wconstr =

𝜕C 𝜕C · 𝛿̃ u+ 𝛿X1 𝜕̃ u 𝜕 Xk1 k

)

· 𝝀 + C · 𝛿𝝀,

(124)

by which the variational formulation of the equations of motion (115) needs to be augmented. If one decides to use finite elements, as we do in our numerical examples, nodal points of elements are naturally placed at boundaries of the non-material domains, which are located at fixed places with respect to the stretched coordinate 𝜉 1 = 𝜉k1 . 3. Algorithmic treatment and semi-discretization In what follows, we introduce the spatial discretization in the variational formulation of the equations of motion (115). The (generalized) displacement field in the i-th variable domain is approximated by a set of interpolation functions represented by the matrix Ni (𝜉 1 ) and the corresponding vector of k generalized coordinates qi (t),

̃ (𝜉 1 , t) = Ni (𝜉 1 ) qi (t), U

Ni ∈ ℝ3×k ,

qi ∈ ℝk .

(125)

At this point, there is no need to further particularize the nature of the spatial discretization and the coordinates. One could employ finite elements, i.e., interpolation functions defined on a local support, or, alternatively, global functions defined on the entire domain  i . However, we do require that the interpolation functions are defined in the stretched coordinate 𝜉 1 separately on each variable domain. Pursuing a Galerkin approach, the same interpolation functions Ni are used for the variations, i.e., the virtual displacements, as for the displacement field itself,

̃ (𝜉 1 , t) = Ni (𝜉 1 ) 𝛿 qi (t), 𝛿U

𝛿 q i ∈ ℝk .

(126)

For the sake of clarity, we proceed term-by-term in our discussion of the semi-discretization of the variational formulation (115). As in our derivations of the continuous problem, we start with the contribution from the inertia forces. Upon the spatial discretization, a first term of the variation of the kinetic energy with respect to the generalized displacement is written as

𝛿u Ti,1 =

𝜉i1+1

∫𝜉 1 i

̃ 𝜕2U ̃ J−1 d𝜉 1 = 𝛿 qT li Mi q̈ i , · I · 𝛿U i i l (0) 𝜕 t2 𝜌 i

Mi =

𝜉i1+1

∫𝜉 1

NTi I𝜌 Ni d𝜉 1 ,

Mi ∈ ℝk×k ,

(127)

i

where the constant symmetric k × k matrix Mi = MTi has been introduced. Similar to a conventional approach based on the material coordinate X1 , Mi represents the mass matrix. Unlike conventional Lagrangian formulations, however, the deformationdependent ratio of the domain’s current and initial lengths li ∕li (0) enters the above mass term as a coefficient in the sliding-beam

20

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

formulation. The second Coriolis-type (gyroscopic) inertia term in (115) is first decomposed into two equal terms with integration by parts being subsequently performed on only one of them. Substituting the definition (31), we obtain 𝜉i1+1

𝛿u Ti,2 =

∫𝜉 1 i

[

𝜉i + 1 ̃ ̃ 𝜕2U 𝜕2U ̃ ̃vi J−1 d𝜉 1 + ̃ ̃vi J−1 d𝜉 1 · I𝜌 · 𝛿 U · I · 𝛿U i i 1 ∫𝜉 1 𝜕𝜉 𝜕 t 𝜕𝜉 1 𝜕 t 𝜌 1

i

̃ 𝜕U ̃ J −1 ̃v · I · 𝛿 U i 𝜕t i 𝜌

=

𝜉i1+1

−

=

+ 𝜉i1

𝜉i1+1

𝜉i + 1 ̃ ̃ −1 1 ̃ 𝜕2U 𝜕U 𝜕𝛿 U ̃ J −1 d 𝜉 1 − ̃ ̃v · I · v · I · 𝛿 U J d𝜉 i 𝜌 i ∫𝜉 1 𝜕𝜉 1 𝜕 t 𝜕 t i 𝜌 𝜕𝜉 1 i 1

∫𝜉 1 i

i

(128)

̃ 𝜕̃vi 𝜕U ̃ J −1 d 𝜉 1 · I · 𝛿U i 𝜕 t 𝜕𝜉 1 𝜌

∫𝜉 1 i

[

]𝜉i1+1

̃ 𝜕U ̃ J −1 ̃v · I · 𝛿 U i 𝜕t i 𝜌

]𝜉i1+1 𝜉i1

{

+ 𝛿 qTi

) ) Ẋ 1i+1 ( Ẋ 1 ( l̇ Vi − VTi − i Wi − WTi + i Mi li (0) li (0) li (0)

} q̇ i ,

with the k × k matrices Vi and Wi being introduced as Vi =

𝜉i1+1

NTi I𝜌

∫𝜉 1 i

) 𝜕 Ni ( 1 𝜉i − 𝜉 1 d𝜉 1 , 1 𝜕𝜉

Wi =

𝜉i1+1

∫𝜉 1

NTi I𝜌

i

) 𝜕 Ni ( 1 𝜉i+1 − 𝜉 1 d𝜉 1 . 1 𝜕𝜉

(129)

Contrary to our intuition from elementary kinematics of relative motion, we note that the semi-discrete Coriolis-type term is not entirely anti-symmetric. If the motion of two domain boundaries does not coincide, i.e., Ẋ 1i ≠ Ẋ 1i+1 ⇔ l̇ i ≠ 0, a symmetric part expressed in terms of Mi enters the convective terms in the equations of motion. In other words, the symmetric term only vanishes if no stretching is introduced instantaneously by the coordinate transformation. We have omitted a discretization of the boundary terms that enter the above relation upon integration by parts. These terms cancel with corresponding terms of ̃∕𝜕 t and the variation 𝛿 U ̃ , see the result (90). At the neighboring domains due to the continuity of the local time-derivative 𝜕 U outer boundaries 𝜉 1 , 𝜉 n+1 , we again assume the motion of the boundaries relative to the material points and/or the displacement

̃ = 0 hold. of the beam to be prescribed kinematically, i.e., 𝛿 X11 = 𝛿 Xn1+1 = 0 and/or 𝛿 U For the third term, we find the semi-discretized representation 𝜉i1+1

𝛿u Ti,3 =

∫𝜉 1 i

̃ 𝜕U 𝜕𝜉 1

(

)

𝜕̃vi 𝜕̃vi ̃ J−1 d𝜉 1 = 𝛿 qT ̃v · I𝜌 · 𝛿 U − i i 𝜕 t 𝜕𝜉 1 i

(

Ẍ 1i+1 li (0)

Vi −

Ẍ 1i li (0)

)

Wi

qi ,

(130)

which contains the second time derivatives of the variable domain boundaries. Substituting the approximations of the generalized displacement and the virtual displacement (125)–(126), a fourth term follows from the variation of the kinetic energy with respect to the generalized displacement as 𝜉i1+1

𝛿u Ti,4 =

∫𝜉 1 i

( )2 ⎫ ⎧ ( ̇ 1 )2 1 ̇1 X Ẋ 1i ̇ X X ̃ ̃ 2 −1 1 ⎪ i+1 𝜕U 𝜕𝛿 U i+1 i T ⎪ ̃ · I · v J d 𝜉 = 𝛿 q R − 2 S + Ti ⎬ qi . i i i ⎨ l (0) l 𝜕𝜉 1 𝜌 𝜕𝜉 1 i i l ( 0 ) l l ( 0 ) l i i i i i i ⎪ ⎪ ⎩ ⎭

(131)

Here, we have introduced three symmetric k × k matrices Ri , Si , Ti that are defined as Ri =

𝜉i1+1

(

∫𝜉 1 i

Ti =

𝜉i1+1

(

∫𝜉 1 i

𝜕 Ni 𝜕𝜉 1 𝜕 Ni 𝜕𝜉 1

)T

I𝜌

)T I𝜌

)2 𝜕 Ni ( 1 𝜉 − 𝜉 d𝜉 1 , 𝜕𝜉 1 i

𝜕 Ni (

𝜕𝜉 1

𝜉i1+1

−𝜉

)2

Si =

𝜉i1+1

(

∫𝜉 1 i

𝜕 Ni 𝜕𝜉 1

)T

I𝜌

)( ) 𝜕 Ni ( 1 1 𝜉 − 𝜉 𝜉 − 𝜉 d𝜉 1 , i+1 𝜕𝜉 1 i (132)

d𝜉 . 1

The semi-discrete expressions (127), (128), (130), (131) are generalizations of the virtual work of the inertia forces in the formulation presented in Ref. [40] with respect to two aspects: While Euler-Bernoulli beams were considered in Ref. [40], shear deformation is included in the proposed approach. More importantly, arbitrarily many variable domains are regarded in the present formulation, whereas the original sliding spaghetti problem [2,3] and the extended sliding spaghetti problem in Refs. [40,41] only exhibit one respectively two variable domains. The key idea of the present formulation is to include variations with respect to the domain boundaries in the variational formulation. The spatial discretization of the additional terms, by which the variation of the kinetic energy has been augmented, is considered next. For the first term in 𝛿 l Ti , we obtain the semi-discretized expression

𝛿l Ti,1 =

𝜉i1+1

∫𝜉 1 i

̃ ̃ 𝜕2U 𝜕U ·I · 𝛿 𝜒̃ J−1 d𝜉 1 = qTi 𝜕 t2 𝜌 𝜕𝜉 1 l i

(

𝛿 Xi1+1 li (0)

VTi −

𝛿 Xi1 li (0)

)

WTi

q̈ i ,

(133)

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

21

where the matrices Vi and Wi re-appear once we pull out the vector of generalized coordinates qi . ̃ and 𝛿 U ̃ along with the definition of ̃vi (31) in the second term, we recover the symmetric Inserting the approximations for U matrices Ri , Si and Ti :

𝛿l Ti,2 = 2

𝜉i1+1

∫𝜉 1 i

̃ ̃ 𝜕2U 𝜕U ̃v · I · 𝛿 𝜒̃J−1 d𝜉 1 = 2 q̇ Ti 𝜕𝜉 1 𝜕 t i 𝜌 𝜕𝜉 1 l i

(

𝛿 Xi1+1 Ẋ 1i+1 li (0) li

Ri −

𝛿 Xi1+1 Ẋ 1i + 𝛿 Xi1 Ẋ 1i+1 li (0) li

Si +

𝛿 Xi1 Ẋ 1i li (0) li

)

Ti

qi .

(134)

From the third term resulting from the variation with respect to the variable boundaries, we obtain the symmetric relation

𝛿l Ti,3 =

𝜉i1+1

̃ 𝜕̃vi ̃ 𝜕U 𝜕U ·I · 𝛿 𝜒̃ J−1 d𝜉 1 𝜕𝜉 1 𝜕 t 𝜌 𝜕𝜉 1 l i

∫𝜉 1 i

(

= qTi

Ẍ 1i+1 𝛿 Xi1+1 li (0) li

(

̇

l qTi i li

−

Ri −

Ẋ 1i+1 𝛿 Xi1+1 li (0) li

Ẍ 1i+1 𝛿 Xi1 li (0) li

Si −

Ẋ 1i+1 𝛿 Xi1

Ri −

li (0) li

Ẍ 1i 𝛿 Xi1+1 li (0) li

Si −

Si +

Ẋ 1i 𝛿 Xi1+1 li (0) li

)

Ẍ 1i 𝛿 Xi1

Ti

li (0) li

Ẋ 1i 𝛿 Xi1

Si +

li (0) li

(135)

qi

) qi .

Ti

The fourth term apparently contains an odd number of spatial derivatives. Therefore, we introduce the additive splitting and selective integration by parts once again, by which we obtain the following semi-discrete representation:

𝛿l Ti,4 =

𝜉1

𝜉1

̃ ̃ ̃ ̃ i+ 1 𝜕 U i+ 1 𝜕 U 1 𝜕2U 1 𝜕2U ·I · 𝛿 𝜒̃ J−1 d𝜉 1 + ·I · 𝛿 𝜒̃ J−1 d𝜉 1 2 ∫𝜉 1 𝜕 t 𝜌 𝜕𝜉 1 𝜕 t l i 2 ∫𝜉 1 𝜕 t 𝜌 𝜕𝜉 1 𝜕 t l i i i [

̃ ̃ 𝜕U 𝜕U = · I𝜌 · 𝜕t 𝜕t

{

}]𝜉 1

) 𝛿X1 ( ) 𝛿 Xi1+1 ( 1 𝜉i − 𝜉 1 − i 𝜉i1+1 − 𝜉 1 li (0)

li (0)

(

)

i+ 1

𝜉i1

1 𝛿 Xi+1 − 𝛿 Xi T + q̇ i Mi q̇ i . 2 li (0) 1

(136)

1

To obtain the above relation, we have used the property vT A − AT v = 0, which holds for arbitrary vectors v and matrices A.

̃∕𝜕 t, the boundary terms in the above relation cancel upon summation of Due to the continuity of the local time derivative 𝜕 U adjacent domains. At the outer boundaries 𝜉11 and 𝜉n1+1 , in turn, the corresponding variations are required to vanish identically,

i.e, 𝛿 X11 and 𝛿 Xn1+1 = 0. The first of the remaining two terms turns out quadratic in the generalized velocities q̇ i , 𝜉i + 1 ̃ ̃ 𝛿 li −1 1 1 𝛿 Xi1+1 − 𝛿 Xi1 T 1 𝜕U 𝜕U · I𝜌 · J d𝜉 = q̇ i Mi q̇ i , 2 ∫𝜉 1 𝜕t 𝜕 t li i 2 li (0) 1

𝛿l Ti,5 =

(137)

i

whereas the second term is quadratic in the coordinates qi , 𝜉i + 1 ̃ ̃ 2 𝛿 Xi1+1 − 𝛿 Xi1 −1 1 1 𝜕U 𝜕U · I𝜌 · 1 ̃ v Ji d 𝜉 1 ∫ 2 𝜉1 𝜕𝜉 𝜕𝜉 i li 1

𝛿l Ti,6 =

i

⎧(

1 1 1 𝛿 Xi+1 − 𝛿 Xi T ⎪ = qi ⎨ 2 li (0) ⎪

⎩

Ẋ 1i+1 li

)2 Ri − 2

Ẋ 1i+1 Ẋ 1i li 2

( Si +

(138)

)2 ⎫ ⎪ Ti ⎬ qi . li ⎪ ⎭

Ẋ 1i

In a next step, we collect all terms that contain second time derivatives, by which we obtain the mass terms of the total variation of the kinetic energy:

𝛿 Tim = +

1 li (0)

{

(

)

(

) }

𝛿 qTi li Mi q̈ i + 𝛿 qTi Ẍ 1i+1 Vi − Ẍ 1i Wi qi + qTi 𝛿 Xi1+1 VTi − 𝛿 Xi1 WTi q̈ i (

)

1 qT Ẍ 1i+1 𝛿 Xi1+1 Ri − Ẍ 1i+1 𝛿 Xi1 Si − Ẍ 1i 𝛿 Xi1+1 Si + Ẍ 1i 𝛿 Xi1 Ti qi . li (0) li i

(139)

For a more compact representation, we can augment the vector of generalized coordinates qi by the variable domain boundaries Xi1 , Xi1+1 ,

(

Q i = qTi Xi1+1 Xi1

)T

,

(

𝛿 Q i = 𝛿 qTi 𝛿 Xi1+1 𝛿 Xi1

)T

,

(140)

22

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341 aug

which suggests the introduction of an augmented symmetric (k + 2) × (k + 2) mass matrix Mi aug Mi

⎛ li 2 Mi 1 ⎜ T T = ⎜ l q V li (0) li ⎜ i i i T T ⎝−li qi Wi

−li Wi qi ⎞ ⎟ −qTi Si qi ⎟ , ⎟ qTi Ti qi ⎠

li Vi qi qTi Ri qi −qTi Si qi

as

(141)

such that equation (139) simplifies to

𝛿 Tim = 𝛿 Q Ti Maug Q̈ i. i

(142) aug

If the variations 𝛿 Xi1 , 𝛿 Xi1+1 are not included, only the first row in the augmented matrix Mi

remains. In view of their algebraic

̃ (see Section 2.5) do not add enter nature, constraints that specify domain boundaries Xi1 , Xi1+1 in terms of the displacement field U

̃ are considered in Hamilton’s the mass terms. For this reason, the mass matrix is generally non-symmetric if only variations 𝛿 U principle as in our original approach [40,41], which was based on the principle of virtual work. Following the reasoning of Vu-Quoc and Li [2], terms containing the (generalized) velocities are referred to as velocityconvection, which can be collected in the following expression: 𝛿 Tivc =

1 𝛿 qT li (0) i

+

2 q̇ T li (0) li i

−

+

l̇ i li (0) li

2

{

(

Ẋ 1i+1 Vi − VTi

{

qTi

)

( ) } − Ẋ 1i Wi − WTi + l̇ i Mi q̇ i

(

)

}

𝛿 Xi1+1 Ẋ 1i+1 Ri − 𝛿 Xi1+1 Ẋ 1i + 𝛿 Xi1 Ẋ 1i+1 Si + 𝛿 Xi1 Ẋ 1i Ti qi

{

(

)

}

𝛿 Xi1+1 Ẋ 1i+1 Ri − 𝛿 Xi1+1 Ẋ 1i + 𝛿 Xi1 Ẋ 1i+1 Si + 𝛿 Xi1 Ẋ 1i Ti qi

1 1 1 𝛿 Xi+1 − 𝛿 Xi T qi 2 li (0) li 2

{(

Ẋ 1i+1

)2

(

Ri − 2 Ẋ 1i+1 Ẋ 1i Si + Ẋ 1i

(143)

)2 } Ti

qi

( ) aug,asym aug,sym = 𝛿 Q Ti Vi + Vi Q̇ i . By introducing the augmented vector of generalized coordinates Q i , we can represent the velocity-convection in terms of an anti-symmetric matrix and a symmetric matrix, which are defined as aug,asym Vi

( ) ( ) ⎛Ẋ 1 Vi − VTi − Ẋ 1i Wi − WTi i+1 1 ⎜ = ⎜ 0 li (0) ⎜ 0 ⎝

0

0⎞

0

0⎟

0

⎟

(144)

⎟ 0⎠

and 0 0 ⎛l̇ i li Mi ⎞ {( ) } 1 1 ⎜ ⎟ ̇ ̇ ̇ X X l ⎜ 0 − i + Ẋ 1i Ri − Ẋ 1i Si qi −2li q̇ Ti Si qi + l̇ i qTi Si qi − i+1 qTi Ri qi + i qTi Ti qi ⎟ . 2li q̇ Ti Ri qi + qTi ⎜ ⎟ 2 2 2 {( ) } ⎟ ⎜ ̇1 ̇1 ̇ X X l ⎜ 0 −2li q̇ T Si qi + l̇ i qT Si qi − i+1 qT Ri qi + i qT Ti qi 2li q̇ T Ti qi + qT − i − Ẋ 1i+1 Ti + Ẋ 1i+1 Si qi ⎟⎠ i i i i i i ⎝ 2

aug,sym

Vi

=

1 li (0) li 2

2

2

2

(145) The anti-symmetric matrix does not contain terms related to the variation of the variable domain-length 𝛿 li , which, in turn, are strongly represented in the symmetric part of the velocity-convection. The symmetric part further contains a multiple of the mass-type matrix Mi in the upper left, which vanishes identically if the (material) length of the variable domain is preserved, i.e., l̇ i = 0. The remaining terms form the stiffness-convection part of the variation of the kinetic energy,

𝛿 Tisc = −

1 𝛿 qTi li (0) li

{(

Ẋ 1i+1

)2

(

Ri − 2Ẋ 1i+1 Ẋ 1i Si + Ẋ 1i

)2 } Ti

aug

which can be expressed in terms of a symmetric matrix Si

aug

Si

( )2 ⎛( ̇ 1 )2 1 1 1 ⎜ X i+1 Ri − 2Ẋ i+1 Ẋ i Si + Ẋ i Ti 1 ⎜ =− 0 li (0) li ⎜ ⎜ 0 ⎝

aug

qi = 𝛿 Q Ti Si

Q i,

(146)

̃, , which only contains the variation 𝛿 U

0 0 0

⎞

0⎟

⎟

0⎟ . 0⎟⎠

(147)

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

23

The attribution of individual terms to velocity-convection and stiffness-convection, respectively, may appear somewhat ambiguous, and, indeed, there is little physical underpinning as far as variations in the variable domain boundaries 𝛿 Xi1 , 𝛿 Xi1+1 are concerned. For terms related to the variation of the kinetic energy with respect to the generalized displacement, the respec̃ and 𝛿 U ̃ determine whether a term belongs to velocity or stiffness convection following the reasoning in tive derivatives of U Ref. [2]. We choose to attribute terms related to 𝛿 Xi1 , 𝛿 Xi1+1 either to mass terms or velocity convection terms according to the time derivatives involved. We proceed with the spatial discretization of the variation of the work performed by the internal forces. The variation with respect to the generalized displacement (81) gives

𝛿u Wiint =

𝜉i1+1

{

(

)

̃ 𝜕̃r0 ̃ 𝜕𝛿 U ̃ F· − 𝛿 𝜃̃ × f · e3 𝜕𝜉 1 𝜕𝜉 1

∫𝜉 1 i

} d𝜉 1 = 𝛿 qTi Fui ,

Fui ∈ ℝk ,

(148)

where Fui represents the corresponding k-dimensional discrete vector of generalized forces. The components of Fui are computed

just as in a conventional material formulation except that the interpolation functions are functions of 𝜉 1 rather than X1 . In the present formulation, additional terms emerge from the variation with respect to the boundaries of variable domains (101). We and 𝛿l Wiint , where the first one can be expressed in terms of a force vector Fli and the split these terms into two parts 𝛿l Wiint ,1 ,2 generalized coordinates qi as

𝛿l Wiint ,1 = −

𝛿 Xi1+1 − 𝛿 Xi1 li

𝜉i1+1

∫𝜉 1 i

𝛿 X − 𝛿 Xi T l ̃ 1 𝜕U ̃ F· d𝜉 = − i+1 qi Fi , 𝜕𝜉 1 li 1

𝜉i1+1

1

Fli =

(

∫𝜉 1 i

𝜕 Ni 𝜕𝜉 1

)T

̃ F d𝜉 1 .

(149)

The second part can be represented as a multiple of the work performed by internal forces in the i-th domain Wiint ,

𝛿l Wiint = ,2

𝜉i1+1

𝛿 Xi1+1 − 𝛿 Xi1 1

2 ∫𝜉 1 i

li

̃ J −1 d 𝜉 1 = ̃ F·𝜞 i

𝛿 Xi1+1 − 𝛿 Xi1 li

Wiint .

(150)

The augmented vector (140) of generalized coordinates Q i suggests the introduction of an augmented vector of internal forces aug Fi such that the total variation of the work of the internal forces reads

(

𝛿 Wiint = 𝛿 Q Ti Faug , i

aug

Fi

=

(Fui )T

−

1 T l 1 int q F + W li i i li i

1 T l 1 int q F − W li i i li i

)T

.

(151)

When using Newton’s method to iteratively solve the non-linear algebraic problems that result from the spatial and possibly temporal discretization, the tangent stiffness matrix has to be computed by linearization of the (augmented) vector of internal aug forces Fi with respect to the generalized coordinates Q i . If the variations with respect to the variable domain boundaries are not included, as it was the case in our original formulation [40,41], the tangent stiffness matrix turns out to be non-symmetric in analogy to what we have observed for the mass matrix above. For the sake of brevity, however, we omit the linearization of the governing equations of the proposed formulation. To conclude the discussion of the semi-discretization, we consider the variation of the work performed by external loads. ext,u The variation in the (generalized) displacement (83) gives the vector Fi

𝛿u W ext =

𝜉i1+1

∫𝜉 1

̃ J−1 d𝜉 1 = 𝛿 qT Fext,u , ̃ F ext · 𝛿 U i i i

ext,u

Fi

𝜉1

=

i

i+ 1 li NTi ̃ F ext d𝜉 1 . li (0) ∫𝜉 1

(152)

i

ext,u

A first term, which results from the variation in the variable boundaries (102), can be written in terms of Fi

𝛿l Wiext ,1 =

𝛿 Xi1+1 − 𝛿 Xi1 li

𝜉i1+1

∫𝜉 1

̃ ̃ J −1 d 𝜉 1 = F ext · U i

𝛿 Xi1+1 − 𝛿 Xi1 li

i

ext,u

qTi Fi

.

as (153)

F ext ∕𝜕𝜉 1 , we obtain the semi-discrete representation as For the second term, which involves the derivative 𝜕̃

𝛿l Wiext ,2 = −

𝜉i1+1

∫𝜉 1 i

𝜕̃ F ext ̃ ,l ,l · U 𝛿l 𝜒̃ Ji−1 d𝜉 1 = 𝛿 Xi1+1 qTi Fext + 𝛿 Xi1 qTi Fext . i,1 i,2 𝜕𝜉 1 ext,l

(154)

ext,l

In the above relation, the generalized force vectors Fi,1 and Fi,2 are introduced as ext,l

Fi,1 = −

𝜉1

i+ 1 li 𝜕̃ F ext NTi li (0) ∫𝜉 1 𝜕𝜉 1 i

(

)

𝜉i1 − 𝜉 1 d𝜉 1 ,

ext,l

Fi,2 =

𝜉1

i+ 1 li 𝜕̃ F ext NTi li (0) ∫𝜉 1 𝜕𝜉 1 i

(

)

𝜉i1+1 − 𝜉 1 d𝜉 1 .

(155)

24

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

Fig. 4. Undeformed reference configuration of a rod under a uniformly distributed tensile load f01 and a compressive tip-force F01 . The structure is partitioned into three variable domains of initial lengths l1 (0) = L∕2 and l2 (0) = l3 (0) = L∕4. The interface points between the variable domains remain spatially fixed in the course of motion.

In terms of the augmented vector of generalized coordinates Q i , the total variation of the work of the external forces follows from the combination of the three parts as

𝛿u W ext = 𝛿 Q Ti Fiext,aug ,

ext,aug

Fi

(

=

(Fiext,u )T

1 T ext,u ,l q F + qTi Fiext ,1 li i i

−

1 T ext,u ,l q F + qTi Fiext ,2 li i i

)T

,

(156)

ext,aug

where Fi denotes the augmented vector of external forces. Eventually, we can put together all terms discussed in the present section to obtain the (semi-)discretized representation of the variational formulation (115) of the proposed formulation: n ∑

{

(

)

}

aug,sym aug,asym aug aug ext,aug 𝛿 Q Ti Maug Q̈ i + Vi + Vi Q̇ i + Si Q i + Fi − Fi + 𝛿 Wconstr = 0. i

(157)

i=1

In view of the algebraic nature of the constraint equations 𝛿 Wconstr , a finite-dimensional system of differential algebraic equations (DAEs) is obtained from the above scalar-valued relation by requiring 𝛿 Q i to be compatible with kinematic boundary conditions but arbitrary otherwise. Note that all coefficient matrices in (157) and the vectors of internal and external forces depend on both the current state of deformation and the motion of the variable boundaries Xi1 . In transient problems, we need to make use of some appropriate time-integration scheme to solve the system of DAEs. Stiff problems and systems of DAEs typically demand implicit schemes, see, e.g., Hairer and Wanner [52]. Using implicit schemes, the non-linearity of the governing equations in the proposed formulation requires a non-linear system of algebraic equations to be solved in every time step, e.g., by Newton’s method. 4. Example problems 4.1. Implementation of the proposed approach The proposed extension of the sliding-beam formulation is implemented in the open-source simulation software HOTINT [53,54], in which the original approach [40,41] has also been realized. For this purpose, an implementation of the geometrically exact beam theory of Simo and Vu-Quoc [29,30] is modified to provide the additional terms required in the present formulation. We use two-noded finite elements with linear shape functions for the interpolation of the components of the generalized dis̃. Note ̃ , i.e., the axial displacement ̃ u1 , the transversal deflection ̃ u2 and the cross-sectional angle of rotation 𝜃 placement field U that the stretched coordinate 𝜉 1 is discretized in the sliding-beam formulation rather than the material coordinate X1 , which is used in conventional problems. To eliminate shear locking, reduced integration with a single integration point is applied to terms related to the work of internal forces. Time-integration is performed by means of the trapezoidal rule [52] using constant step sizes, where the modified Newton method [55] is used to solve the resulting non-linear system of equations, i.e., the tangent stiffness matrix is not updated in every iteration but according to an appropriate metric. We refer the interested reader to the documentation of HOTINT (www.hotint.org) for further details. 4.2. Axial vibrations of a rod As a first exemplary problem, we study axial vibrations of a rod subjected to a uniformly distributed tensile load f01 and a

compressive tip-force F01 , which is applied at its right end X1 = L, see Fig. 4. The left end of the rod is fixed, i.e., u1 (X 1 = 0, t ) = 0.

(158)

We consider the entire structure of length L, which is partitioned into three domains  =  =  1 ∪  2 ∪  3 . The length of the first domain from the left-hand end of the structure in the undeformed reference configuration is half of the rod’s total

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

25

Table 1 Geometry, material properties and external loads.

total length initial length domain 1 initial length domain 2 initial length domain 3

L l1 (0) l2 (0) l3 (0)

2m 1m 0.5 m 0.5 m

extensional stiffness translatory inertia tip force distributed load

EA

100 N 2 kgm−1 −5 N 10 Nm−1

𝜌A

F01 = −EA∕20 f01 = EA∕(5l)

length l1 (0) = L∕2. The remaining part is again divided into two equally long parts, i.e., the domains have the dimensions l2 (0) = l3 (0) = L∕4. All geometric and material parameters { } along with the external loadings are listed in Table 1. Unlike conventional approaches, the three domains  i , i = 1, 2, 3 are not fixed to the material points X1 of the structure. Only the domain boundaries at the rod’s ends are fixed, i.e., X11 (t ) = X11 (0) = 0, X41 (t ) = X41 (0) = L, whereas we assume the inter-

face points between two adjacent variable domains to move relative to the material points of the structure, i.e., X21 (t ) = l1 (t ),

X31 (t ) = l1 (t ) + l2 (t ). By means of the piecewise linear coordinate transformation (17), we introduce the stretched coordinate 𝜉 1

̃ i , i = 1, 2, 3, which coincide with the such that the generally variable domains  i , i = 1, 2, 3 are mapped onto fixed domains  variable domains in the initial undeformed configuration: [

X 1 ∈  1 = X11 (t ) = 0, X21 (t ) = l1 (t )

[

X 1 ∈  2 = X21 (t ) = l1 (t ), X31 (t ) = l1 (t ) + l2 (t )

[

X 1 ∈  3 = X31 (t ) = l1 (t ) + l2 (t ), X41 (t ) = L

] ]

]

[

]

↔

̃ 1 = 𝜉 1 = 0, 𝜉 1 = l1 (0) , 𝜉1 ∈  1 2

↔

̃ 2 = 𝜉 1 = l1 (0), 𝜉 1 = l1 (0) + l2 (0) , 𝜉1 ∈  2 3

↔

̃ 3 = 𝜉 1 = l1 (0) + l2 (0), 𝜉 1 = L . 𝜉1 ∈  3 4

[ [

]

(159)

]

Besides the two terminal points, the nature of the variable domains still remains to be specified. In the present problem, we require the variable domain boundaries to remain at spatially fixed points throughout the course of motion. Mathematically, this idea translates into four constraint equations Ci = 0 for the four domain boundaries Xi1 , i = 1, … , 4: C1 = X11 (t ) = 0,

C2 = X21 (t ) + u1 (X 1 = X21 (t ), t ) − l1 (0) = 0,

C3 = X31 (t ) + u1 (X 1 = X31 (t ), t ) − l1 (0) − l2 (0) = 0,

C4 = X41 (t ) − L = 0.

(160)

Once we introduce the coordinate transformation (17), the above constraint equations no longer involve the displacements at variable points, but at constant domain boundaries in the stretched coordinate 𝜉i1 , i = 1, … , 4: C1 = X11 (t ) = 0,

C2 = X21 (t ) + ̃ u1 (𝜉 1 = 𝜉21 , t ) − l1 (0) = 0,

C3 = X31 (t ) + ̃ u1 (𝜉 1 = 𝜉31 , t ) − l1 (0) − l2 (0) = 0,

C4 = X41 (t ) − L = 0.

(161)

We deliberately impose a distributed load that is constant and a concentrated load that acts at a domain boundary, which is fixed to a material point, i.e., the rod’s right-hand tip. For this choice, no dependency on the stretched coordinate is introduced in u1 (𝜉 1 = L, t) the external loadings. Therefore, we can immediately compare the results for the tip-trajectories u1 (X1 = L, t) = ̃ of the proposed formulation to trajectories computed with the conventional geometrically exact beam theory (GEBT) of Simo and Vu-Quoc [29,30]. In Fig. 5, tip-trajectories u1 (L, t) obtained from both approaches are illustrated for a coarse discretization (5a) and a refined mesh (5b). A constant time step size of Δt = 1 × 10−3 s has proven to be sufficiently small in the present problem. For both the coarse and the fine mesh, we adopt the same uniform discretization with respect to the material coordinate X1 and the stretched coordinate 𝜉 1 , i.e., the nodal points in the undeformed reference configuration coincide in both approaches. In view of the comparatively large deformation, the tip-trajectories obtained with coarse meshes of eight finite elements each (Fig. 5a) already show a good agreement. Using 256 finite elements, the tip-trajectories coincide up to the specified numerical precision of the Newton iteration, see Fig. 5b. The proposed formulation allows us to study the deformation at spatially fixed points by placing domain boundaries at these ̃ 1 and points. By means of the constraints (161), we have introduced the variable domains such that the boundaries between  1 1 ̃ ̃ ̃  2 as well as  2 and  3 remain located at x = l1 (0) = L∕2 and x = l1 (0) + l2 (0) = 3L∕4, respectively, throughout the u1 of those material points which are entire course of motion. In other words, we compute the evolution of the displacement ̃ momentarily located at the domain boundaries. Note that the nodal points within each domain are neither spatially fixed nor u1 (𝜉 1 = L∕2, t) and ̃ u1 (𝜉 1 = 3L∕4, t) are illustrated do they correspond to fixed material points in general. In Fig. 6, trajectories ̃ u1 (L, t). The coarse discretization (Fig. 6a) is already able to reproduce the qualitative along with the tip-displacement u1 (L, t) = ̃ behavior. For the fine mesh, we again obtain converged results, see Fig. 6b. We note that a comparison to the conventional geometrically exact beam theory is not immediately feasible for the variable domain boundaries. For a comparison, one would have to determine the respective finite elements and local element coordinates of the material points momentarily located at the spatially fixed points corresponding to 𝜉 1 = L∕2 and 𝜉 1 = 3L∕4 in every time step.

26

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

Fig. 5. Comparison of the tip-displacement u1 (L,t) obtained with the present sliding-beam formulation (SBF) and the geometrically exact beam theory (GEBT): (a) coarse discretization, (b) fine discretization.

Fig. 6. Displacement ̃ u1 (𝜉 1 , t) at spatially fixed domain boundaries 𝜉 1 = L∕2 and 𝜉 1 = 3L∕4 as well as the beam’s free tip 𝜉 1 = L: (a) coarse discretization, (b) fine discretization.

4.3. Sliding spaghetti problem In the original sliding-beam problem [2,3], the (transverse) vibrations of a beam that is being retrieved into (or deployed from) a guide are studied. The resemblance to a spaghetti sucked into one’s mouth has coined the term “sliding spaghetti problem”. The part of the beam yet to be retrieved can be viewed as a cantilever whose length decreases in the course of the retrieval process. At the same time, it is observed that the vibrations of the beam’s outside part gradually increase in both frequency and amplitude. The original problem and the extended sliding spaghetti problem [40,41] share the same ideas and setup, but they distinguish themselves in an equally subtle as crucial detail: In the original problem, the motion of the beam is completely prescribed kinematically at the opening of the guide. The kinematic constraints at the opening completely decouple the already retrieved part and the outside part kinetically. By prescribing the axial motion, the material point that is momentarily located at the opening is a priori determined throughout the entire retrieval process. The already retrieved part of the beam can therefore be disregarded when considering the vibrations of the outside part. Irrespective of the actual application, however, kinematic constraints that completely prescribe the motion at the opening the guide appear hardly viable from a technological point of view. This limitation motivated our investigations on the so-called extended sliding spaghetti problem [40,41], in which the part of the structure that has already been retrieved into the guide and the outside part still to be retrieved do influence each other. In the extended sliding spaghetti problem, the axial sliding motion was driven either by a force [40] or kinematically, by prescribing the motion of the beam’s already retrieved tip inside the guide [41]. As a consequence, the material point that is momentarily located at the opening of the guide is not a priori known but depends on the beam’s current state of deformation. The variable domains inside and outside the guide are no longer decoupled, but they dynamically interact with one another. In what follows, we adopt the problem setup of [41] to validate the proposed formulation against the original sliding-beam formulation [40], which does not include the variations with respect to the variable domain boundaries. Fig. 7 shows the undeformed reference configuration of the beam as well as a deformed configuration at time t. We assume an initially straight beam

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

27

Fig. 7. Undeformed reference configuration (top) and deformed current configuration (bottom) of the sliding beam: Initially the beam is retrieved by l1 (0) into a guide, which constrains the transversal deflection. The part outside the guide (l2 (0) = L − l1 (0)) is retrieved into the guide by prescribing the motion (displacement) of the left-hand tip u1 (X1 = 0, t). At time t, the material point X 1 = X21 (t) = l1 (t) is currently located at the opening of the guide, i.e., the variable domain inside the guide contains the material points X 1 ∈  1 = [0, l1 (t)], the outside part is formed by the material points X 1 ∈  2 = [l1 (t), L].

Table 2 Geometry, material properties and problem parameters.

total length initial length domain 1 initial length domain 2 height cross-section width cross-section Young’s modulus Poisson’s ratio mass density

L l1 (0) l2 (0) h b E

6m 1m 5m 0.2 m 1m 1 × 107 Nm−2 0 1 × 103 kgm−3

𝜈 𝜌

extensional stiffness shear correction coefficient shear stiffness bending stiffness translatory inertia rotatory inertia distributed load max. acceleration acceleration time

EA k kGA EI 𝜌A 𝜌I f02 a0 T

2 × 106 N 5/6 8.333 × 105 N 6.667 × 103 Nm2 200 kgm−1 6.667 × 10−1 kgm −200 Nm−1 1 ms−2 2s

of length L, which is partially retrieved by l1 (0) into the guide in its undeformed reference configuration, i.e., the part of the beam outside the guide initially has a length of l2 (0) = L − l1 (0). The geometry of the problem, the external loading and the material parameters of the structure are listed in Table 2. The nature of the sliding-beam formulation allows us to describe the interaction of the beam with the guide without a contact model. For this purpose, we introduce two variable domains: The first variable domain  1 comprises all material points that have already been retrieved into the guide. The second domain  2 represents the remaining part of the structure that is currently located outside the guide, where it can undergo large deflections. By means of the coordinate transformation (17), we map the domains that are variable with respect to the material coordinate X1 onto fixed domain in the stretched coordinate 𝜉 1 :

[

]

[

]

X 1 ∈  1 = X11 (t ) = 0, X21 (t ) = l1 (t ) X 1 ∈  2 = X21 (t ) = l1 (t ), X31 (t ) = L

[

]

[

]

↔

̃ 1 = 𝜉 1 = 0, 𝜉 1 = l1 (0) , 𝜉1 ∈  1 2

↔

̃ 2 = 𝜉 1 = l1 (0), 𝜉 1 = L . 𝜉1 ∈  2 3

(162)

The corresponding constraint relations for the variable boundaries of the two domains read C1 = X11 (t ) = 0,

C2 = X21 (t ) + u1 (X 1 = X21 (t ), t ) − l1 (0) = 0,

C3 = X31 (t ) − L = 0.

(163)

Rewritten in terms of the stretched coordinate 𝜉 1 , C1 = X11 (t ) = 0,

C2 = X21 (t ) + ̃ u1 (𝜉 1 = 𝜉21 , t ) − l1 (0) = 0,

C3 = X31 (t ) − L = 0,

(164)

𝜉21

= l1 (0). the constraints only involve the displacement at the opening of the guide, which is always located at the fixed point The guide prohibits the transverse deflection of the part of the structure which has already been retrieved into the guide. Due to the coordinate transformation, the kinematic constraints imposed by the guide on the beam’s transverse deflection can be formulated at a fixed domain in terms of the stretched coordinate 𝜉 1 ,

̃ u2 (𝜉 1 , t ) = 0,

[ ] ∀ 𝜉 1 ∈ 0, l1 (0) .

(165)

The investigations in Ref. [41] have shown that the kinematic constraints at the opening of the guide play a crucial role. We do not impose a constraint on the cross-sections’ angle of rotation at the opening in the present investigations. As a consequence,

28

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

the part of the beam inside the guide is not only subjected to axial strain, but it can also display shear deformation. In a first step, the beam is statically deflected under a uniformly distributed load in transverse direction f02 . Subsequently, the sliding motion is initiated from the static equilibrium configuration. We prescribe the sliding motion at the left-hand tip of the beam rather than at the opening of the guide as in the original spaghetti problem [2]. To set the structure smoothly into motion, a piecewise linear acceleration over the interval t ∈ [0, T] is chosen:

⎧ t −2a0 , ⎪ T ( ) ⎪ ü 1 (X 1 = 0, t ) = ⎨ −2a 1 − t , 0 T ⎪ 0, ⎪ ⎩

T 0≤t< 2 T ≤t
(166)

Already a first static analysis to determine the deflected starting configuration, from which the beam is retrieved into the guide, has revealed that the constraints (165) and the methodology to enforce them demand our attention. According to ̃ vanishes identid’Alembert’s principle, the work performed by constraint forces under admissible virtual displacements 𝛿 U cally. This property is exactly preserved in our algorithmic realization of the constraint equation (165), which is enforced at the respective nodal points of the discretized beam by means of Lagrange multipliers. Similar as for external loadings, we need to consider the variation of the constraint forces with respect to the variable domain boundary, cf. the corresponding relation (102). If we regard the constraint forces as a distributed external load, we deal with a non-uniform load that is only imposed ̃ 1 inside the guide. Under this assumption, neither the general relation (102) nor the particular result onto the variable domain  for distributed loads imposed on parts of the structure (108) let us expect a contribution from the constraint forces, since the conjugate transverse deflection vanishes exactly within the guide. Still, comparisons between the proposed approach and the original sliding-beam formulation [41] have shown non-negligible differences in the solution for the tip-displacement even in static computations, which have not reduced upon a refinement of the discretization, see the last two rows of Table 3. Due to the localized nature of the constraint realization at nodal points, we must treat the constraint forces as concentrated loads. For all nodes fully located within the guide, we do not expect a contribution to the variation of the action integral, since not only the transverse deflection but also its spatial derivative vanish there. At the opening of the guide, however, shear-deformable beams generally exhibit a discontinuity in the spatial derivative of the transverse deflection. For this reason, we need to apply the result for concentrated forces at domain boundaries (112) to the constraint force at the opening. Assume that the k-th node is located u2 (𝜉 1 = l1 (0), t) = 0, at the opening of the guide and let 𝜆k denote the Lagrange multiplier that is used to enforce the constraint ̃ the variation with respect to the domain boundary X21 (t ) follows from (112) as

(

1 𝛿l W = − 𝛿 X21 𝜆k 2 𝜆

𝜕̃ u2 || 𝜕̃ u2 | J1 + J2 1 || | 1 𝜕𝜉 |𝜉 1 =l1 (0)− 𝜕𝜉 |𝜉 1 =l1 (0)+

)

1 2

= − 𝛿 X21 𝜆k

l2 (0) 𝜕̃ u2 || . l2 (t ) 𝜕𝜉 1 ||𝜉 1 =l1 (0)+

(167)

u2 ∕∂𝜉 1 vanishes idenRecall that the transverse deflection is prohibited inside the guide, which is why the left-sided limit of ∂̃ tically. Once the above term is included, the results of the proposed approach and the original sliding-beam formulation [41] converge to one another in both static and dynamic computations, see Table 3 and Fig. 8, respectively. Table 3 also shows the convergence of the tip-displacement in the initial equilibrium configuration against the number of finite beam-elements. Fig. 8 shows the comparison of the axial displacement u1 (L, t) and the transverse deflection u2 (L, t) of the beam’s outside tip in the course of the retrieval process. Using a time step Δt = 5 × 10−4 s and a uniform discretization into 600 elements, i.e., 100 elẽ 1 and 500 elements for the outside part  ̃ 2 , the results obtained from present approach ments for domain inside the guide  and the original formulation almost coincide exactly. Starting from the static equilibrium configuration, the retrieval process induces vibrations whose amplitude—as compared to the decreasing length of the outside part l2 (t)—and frequency gradually increase in the course of motion. In Fig. 9a, the axial and transversal displacements are combined to illustrate the spatial trajectory r0 (L, t) = (L + u1 (L, t))e1 + u2 (L, t)e2 of the outside tip of the beam. The dynamic behavior of the sliding beam is visualized in Fig. 9b, which shows a series of configurations the beam occupies in the course of the retrieval process. The snapshots, which are colored from blue to red for better visibility, represent configurations at an interval of 50 ms starting from the initial state of static equilibrium. Note that the part inside guide, whose opening is located at x1 = 1 m, is not fully displayed for the sake of compactness. 4.4. Axially moving beam As a third example, we consider a problem of an axially moving beam undergoing large deformation. The problem is inspired by one of the examples discussed by Vetyukov [32] in his work on a mixed Eulerian-Lagrangian description for axially moving strings and beams. Though axially moving continua are often considered a separate class of problems, the example demonstrates that the proposed formulation for sliding beams can readily be applied to it. In fact, we show that our formulation unifies the treatment of sliding structures and axially moving continua, which distinguish themselves from the former by a continued nature of the underlying sliding motion. For this purpose, we study the problem illustrated in Fig. 10, for which all parameters are listed in Table 4. Consider a beam, or rather, a part of it, which enters the domain of interest—one could call it control volume—at a spatially fixed support located at the origin of the inertial frame. The beam leaves the domain at a second support,

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

29

Table 3 Axial displacement and transverse deflection of the beam’s axis at the right-hand tip in the static equilibrium configuration.

present

Steinbrecher et al. [41] present, without 𝛿 l W𝜆 (167)

elements

u1 (X1 = L)∕m

u2 (X1 = L)∕m

6 12 24 48 96 192 600 600 600

−0.652448 −0.554604 −0.548869 −0.548806 −0.548790 −0.548785 −0.548784 −0.548784 −0.548704

−2.353257 −2.157029 −2.144838 −2.144753 −2.144742 −2.144724 −2.144722 −2.144722 −2.145017

Fig. 8. Comparison of the tip-displacement obtained with the proposed approach and the original sliding-beam formulation [41]: (a) axial displacement u1 (X1 = L, t), (b) transverse tip-deflection u2 (X1 = L, t).

which is located at a constant horizontal distance D from the former. The vertical position of the right-hand support, however, is not fixed in space. Instead, it is lifted from x2 = 0 to x2 = h(t), with the vertical position being kinematically prescribed as

{1 h(t ) =

(

h0 1 − cos

2 h0 ,

𝜋 t) T

,

0≤t
.

(168)

Initially (t = 0), the structure is in a state of stationary motion, in which it is sliding in its axial direction. Following the example in Ref. [32], we further impose a uniform initial tensile strain 𝜀0 . As we are interested in the mid-span deflection, we

Fig. 9. Comparison of the tip-trajectory obtained with the proposed approach and the original sliding-beam formulation [41] (a); series of configurations occupied by the structure during the retrieval (b).

30

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

Fig. 10. Axially moving beam with a vertically moving support: The structure enters the domain of interest at a spatially fixed support with the velocity v0 . The material points of the beam leave the domain with the same velocity at the right-hand support, which is lifted in the course of motion.

Table 4 Geometry, material properties and problem parameters.

horizontal support distance static tensile strain material length domain 1 material length domain 2 support displacement

D

𝜀0 l 1 ( 0) l2 (0) = D∕(1 + 𝜀0 ) − l1 (0) h0

1m 0.1 0.5 m 0.4091 m 0.4 m

extensional stiffness shear stiffness bending stiffness translatory inertia rotatory inertia

EA kGA EI 𝜌A 𝜌I

1N 5∕ 6 × 1∕ 2 N 0.052 ∕12 Nm2 1 kgm−1 0.052 ∕12 kgm

introduce two variable domains  1 ,  2 with their interface being located at the spatially fixed point x1 = D∕2. The initial strain requires us to distinguish between the undeformed reference configuration and the initial configuration of the transient analysis, which starts from the strained state. To stay consistent with our notation, we denote the quantities in the initial configuration with a subscript ‘0’. The material length of the (part of the) structure initially located within the domain of interest is related to the spatial distance D between the supports and the strain 𝜀0 through L=

D . 1+ 𝜀0

(169)

Initially, the material point X11,0 = 0 is assumed to be located at the left-hand support and X31,0 = L is situated at the right-hand support. The displacement field that corresponds to a uniform tensile strain is given by u10 (X 1 ) = (D − L)

X1 . L

(170)

To determine the material point X21,0 that is located at mid-span, we consider the equation X1 D = X21,0 + u1 (X 1 = X21,0 ) = X21,0 + (D − L) 2,0 . 2 L

(171)

Solving the above equation for X21,0 , we obtain the (somewhat unsurprising) result for the positions of boundaries of the variable domains in the initial configuration: X11,0 = 0,

X21,0 =

L , 2

X31,0 = L.

(172)

The stretched coordinate 𝜉 1 , however, refers to the undeformed reference configuration. By means of the coordinate transformation (17), the two variable domains with respect to the material coordinate are mapped onto corresponding fixed domains

[

]

̃ 1 = 𝜉 1 = 0, 𝜉 1 = D∕2 , 𝜉1 ∈  1 2

[

]

̃ 2 = 𝜉 1 = D∕2, 𝜉 1 = L . 𝜉1 ∈  2 3

(173)

̃ 1 consequently has an undeformed material length of l1 (0) = 𝜉 1 − 𝜉 1 = D∕2, whereas the right domain  ̃2 The left domain  2 1

is l2 (0) = 𝜉31 − 𝜉21 = L − D∕2 long. Upon the coordinate transformation (20), the initial displacement field (170) from the undeformed reference configuration to the stretched initial configuration is obtained as

( ) ⎧ 1 − L 𝜉1, ⎪ D ̃ u10 (𝜉 1 ) = ⎨ ( ) ⎪ D − L D − L − 𝜉1 , ⎩ D − 2L

̃1 𝜉1 ∈  ̃2 𝜉1 ∈ 

.

(174)

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

31

Fig. 11. Mid-span displacement for v0 = 0, 0.1 and 0.2 ms−1 of the axial motion: (a) axial displacement relative to the sliding motion ̃ u1 (𝜉 1 = D∕2, t) − v0 t, (b) transverse u2 (𝜉 1 = D∕2, t). deflection ̃

In the transient analysis, we prescribe kinematic boundary conditions at both supports. The axial displacements of the material points currently located at the supports are assumed as u1 (𝜉 1 = 0, t ) = v0 t , u1 (X 1 = X11 (t ), t ) = ̃

u1 (X 1 = X31 (t ), t ) = ̃ u1 (𝜉 1 = L, t ) = D − L + v0 t .

(175)

With the supports being fixed in space horizontally, the material points currently located at the supports are determined by the constraint equations C1 = X11 (t ) + u1 (X 1 = X11 (t ), t ) = X11 (t ) + v0 t = 0, C3 = X31 (t ) + u1 (X 1 = X31 (t ), t ) − D = X31 (t ) + v0 t − L = 0.

(176)

The position of the supports relative to the material points in the reference configuration therefore follows as X11 (t ) = −v0 t ,

X31 (t ) = L − v0 t .

(177)

Note that v0 is not the absolute velocity of the material points currently located at the supports, but it rather represents a supply rate of material entering and leaving the domain. Due to the axial deformation, supply rates and velocities do not coincide in general, see the transformation of the material time derivative (29). We refer to Vetyukov [32] for a discussion on different boundary conditions and their implications in problems of axially moving continua. Unlike the points at the supports, the material point currently located in the middle is not known in advance but depends on the current state of deformation. The algebraic constraint that defines the point X21 (t ) follows from the requirement r0 (X 1 = X21 (t ), t ) · e1 = D∕2, i.e.,

C2 = X21 (t ) + u1 (X 1 = X21 (t ), t ) −

D D = X21 (t ) + ̃ u1 (𝜉 1 = 𝜉21 , t ) − = 0. 2 2

(178)

The transverse deflection is prohibited at the left support; on the right-hand side, the deflection is prescribed according to the vertical motion of the support given by (168): u2 (𝜉 1 = 0, t ) = 0, u2 (X 1 = X11 (t ), t ) = ̃

u2 (X 1 = X31 (t ), t ) = ̃ u2 (𝜉 1 = L, t ) = h(t ).

(179)

The cross-sectional angle of rotation must vanish at both supports, i.e.,

𝜃 (X 1 = X11 (t), t) = 𝜃̃(𝜉 1 = 0, t) = 0,

𝜃 (X 1 = X31 (t), t) = 𝜃̃(𝜉 1 = L, t) = 0.

(180)

In the numerical analysis, we use 64 beam elements to discretize each of the two variable domains. A time step size of

Δt = 1 × 10−3 s has proven sufficiently small in the transient analysis. Fig. 11 shows the components of the displacement u(𝜉 1 = D∕2, t ) for values v0 = 0, 0.1 and 0.2 ms−1 . For the sake of clarity, the axial at the interface of the variable domains ̃ u1 (𝜉 1 = D∕2, t) − v0 t, is depicted in Fig. 11a rather than the total displacement. displacement relative to the sliding motion, i.e., ̃ u2 (𝜉 1 = D∕2, t) (Fig. 11b) increases with v0 . We further note that the The amplitude of the mid-span transverse deflection ̃ axial sliding motion delays the onset of the vibrations at the midpoint of the domain of interest. Our results (v0 = 0.1 ms−1 ) agree very well with those presented by Vetyukov [32], in which slightly different boundary conditions were used and shear deformation was neglected. For the case v0 = 0.2 ms−1 , a series of configurations the structure occupies in the course motion is illustrated in Fig. 12. As with the spaghetti problem, the color is blended from blue to red to indicate the temporal sequence of the snapshots, which are taken at an interval of 0.5s.

32

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

Fig. 12. Series of configurations occupied by the axially moving beam (v0 = 0.2 ms−1 ) illustrated at an interval of 0.5s.

5. Discussion and conclusion In the present paper, we have proposed a general formulation for the dynamic analysis of sliding beams and axially moving beams. At its heart, our formulation relies on a piece-wise linear transformation that maps variable domains with respect to the beam’s material coordinate X1 onto fixed domains with respect to a new stretched coordinate 𝜉 1 . Reformulating the governing equations in terms of the stretched coordinate facilitates the analysis of sliding beams and axially moving beams, which are characterized by the fact that supports and loads are typically prescribed at variable positions with respect to their material points. In the proposed sliding-beam formulation, it makes no difference whether kinematic constraints and loads are prescribed at variable or fixed material points and domains, whereas a conventional Lagrangian (material) description requires, e.g., a computationally expensive contact approach. The formulation extends our previous approaches [2,40,41] with respect to several important aspects. First of all, we have generalized the coordinate transformation to arbitrarily many variable domains and we have unified the description of sliding beams and axially moving beams. We have considered large (planar) deformations of shear-deformable beams, which are described by the geometrically exact beam theory of Reissner [27] and Antman [28]. In view of possibly large deformations, we have abandoned the common but somewhat artificial assumption that the positions of supports and loads relative to the material points of the structure are prescribed kinematically. As a second key aspect of the proposed formulation, we have generalized Hamilton’s principle by including variable domain boundaries in the variation of the action integral, from which the variational formulation of the equations of motion in the stretched coordinate is derived. The additional terms introduced by variations in the boundaries of variable domains symmetrize the mass matrix and the tangent stiffness matrix as compared to our original formulation [40,41], in which the coordinate transformation was directly applied to the principle of virtual work. Besides the inertia terms and contributions from internal forces, the generalized notion of variation needs to be applied to both the work performed by external loadings and kinematic constraints. To correctly account for loads imposed on variable parts of the structure and concentrated loads at variable points, we have used representations of such loads in terms of the Heaviside function and Dirac’s delta function. We have outlined the spatial semi-discretization of the (generalized) variational formulation in the stretched coordinate 𝜉 1 , which demonstrates the symmetry of the mass matrix. Symmetric mass and tangent stiffness matrices are preferred from a numerical point of view, since both efficiency and convergence are improved and instabilities are less likely to occur. Contrary to the common notion in the kinematics of relative motion, the discretization has also revealed that the velocity-convection term is not anti-symmetric, but it contains a symmetric part if the coordinate transformation changes the length of the variable domains, i.e., a stretching is introduced. Three representative examples show the capabilities and advantages of the proposed formulation. As a first example, we have studied large axial vibrations of a rod under distributed and concentrated forces. The example has enabled us to validate the proposed generalization of the sliding-beam formulation against the conventional geometrically exact beam theory. We have introduced three variable domains to evaluate the displacement at two spatially fixed points along the rod’s axis and as well as its tip. The (extended) sliding spaghetti problem has been our second numerical example, which allows us to compare the proposed formulation to the original approach [40,41], in which variations with respect to the boundaries of variable domains were not included. A pre-tensioned axially moving beam has been studied as a third example, in which vibrations are induced by a vertical motion of one the supports. To evaluate the mid-span deflection, we have partitioned the region between the supports into two spatially fixed variable domains. The variable domains can be regarded as control volumes, which the material points of the structure flow through. The notion of control volumes is familiar from fluid dynamics, in which one typically adopts an Eulerian (spatial) description of the governing equations. Unlike Eulerian formulations for axially moving beams, a Lagrangian (material) description of Hamilton’s principle provides the basis for our formulation. Compared to Eulerian approaches as, e.g., by Vetyukov [32], the

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

33

proposed formulation allows arbitrarily deformed shapes including the formation of loops, in which the uniqueness of the mapping between material and spatial coordinates is generally lost. Vu-Quoc and Li [2] characterized their approach, which can be regarded as a special case of the proposed formulation, as “full Lagrangian”, whereas more recent contributions as, e.g. Refs. [8,9], have used the term “Arbitrary Lagrangian-Eulerian” in related problems. We deliberately evade a discussion and leave the classification to the reader. Though the proposed formulation is advantageous to a conventional Lagrangian description in problems of sliding beams and axially moving beams, since it naturally accounts for the underlying sliding motion, we do mention situations that still require some attention. The coordinate transformation (20) becomes singular if a variable domain collapses to zero material length. In the extended sliding spaghetti problem, for instance, this situation naturally occurs as soon as the beam is fully retrieved into the guide. Suitable approaches to cope with situations of this kind will be subject of future work. CRediT authorship contribution statement Alexander Humer: Conceptualization, Methodology, Software, Validation, Writing - original draft, Investigation, Visualization. Ivo Steinbrecher: Methodology, Writing - review & editing. Loc Vu-Quoc: Methodology, Writing - review & editing, Supervision. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements Alexander Humer and Loc Vu-Quoc acknowledge the support by the COMET-K2 Center of the Linz Center of Mechatronics (LCM) funded by the Austrian federal government and the federal state of Upper Austria. Ivo Steinbrecher contributed to this paper during his time at the Institute of Technical Mechanics at the Johannes Kepler University under the support of the international project “Eulerian Mechanics of Belts” (project no. I 2093) of the Austrian Science Fund (FWF). Appendix A The variation with respect to the variable domain boundaries Xi1 is a crucial ingredient of the proposed formulation, which has neither been present in the original sliding-beam formulation [40,41] nor does it occur in conventional beam formulations in the material coordinate. From a mathematical point of view, we have shown that the additional terms (113) introduced by including variations 𝛿 Xi1 establish the symmetry properties of the inertia terms and the tangent stiffness matrix. To find out about their physical meaning, we first consider the variations of the kinetic energy of the i-th domain  i with respect to its variable boundaries (94). As physical quantities are associated with the material points that form the structure, we revert to the material coordinate X1 using the transformation relation (20):

𝛿l T i =

X1

𝜕U U̇ · I𝜌 · 𝜕X1

i+ 1

∫X 1 i

X1

−

i+ 1

∫X 1 i

𝜕U U̇ · I𝜌 · 𝜕X1

X1

+

i+ 1

∫X 1 i

{

𝛿 Ẋ 1i+1 (

Xi1

li

{

−X

1

)

−

𝛿 Ẋ 1i (

Xi1+1

li

−X

1

)

}

) Ẋ 1 ( ) Ẋ 1i+1 ( 1 Xi − X 1 − i Xi1+1 − X 1 li li

}

dX 1

𝛿 Xi1+1 − 𝛿 Xi1 li

dX 1

(A.1)

𝛿 Xi1+1 − 𝛿 Xi1 1 1 ̇ U · I𝜌 · U̇ dX . 2 li

To form the action integral, we integrate the above relation over time, by which we obtain t1

∫t0

𝛿l Ti dt =

t1

∫t0 t1

−

[

𝜕U U̇ · I𝜌 · 𝜕X1 X1 i+ 1

∫t0 ∫X 1 i

t1

+

∫t

0

X1

i+ 1

∫X 1 i

(

] X 1 =X 1

Ẋ 1i+1 𝛿 Xi1+1 dt −

i+ 1

̇ ̈ · I𝜌 · 𝜕 U + U̇ · I𝜌 · 𝜕 U U 1 𝜕X 𝜕X1

t1

∫t0

){

[

𝜕U U̇ · I𝜌 · 𝜕X1

𝛿 Xi1+1 li

] X 1 =X 1

(

Xi1

−X

Ẋ 1i 𝛿 Xi1 dt

i

1

)

−

𝛿 Xi1 ( li

Xi1+1

−X

1

)

} dX 1 dt

(A.2)

𝛿 Xi1+1 − 𝛿 Xi1 1 1 ̇ U · I𝜌 · U̇ dX dt, 2 li

where we have already performed integration by parts on the first integral of (A.1) to eliminate the time derivatives of the variations of the domain boundaries. Assuming variations to vanish at the beginning and the end of the time interval considered,

34

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

i.e., 𝛿 Xi1 (t0 ) = 𝛿 Xi1 (t1 ) = 𝛿 Xi1+1 (t0 ) = 𝛿 Xi1+1 (t1 ) = 0, all (temporal) boundary terms vanish. Note that the limits of the integrals are no longer constant since we have returned to the material coordinate X1 . Therefore, we use Leibniz’ rule when pulling the time derivative in front, by which spatial boundary terms enter the above relation. As a next step, integration by parts is performed in the material coordinate X1 , i.e., X1

{

𝜕 U̇ U̇ · I𝜌 · 𝜕X1

i+ 1

∫X 1 i

𝛿 Xi1+1 (

)

Xi1 − X 1 −

li

𝛿 Xi1 ( li

Xi1+1 − X 1

)

}

dX 1 (A.3)

X 𝛿 Xi1+1 − 𝛿 Xi1 1 ] ] i+ 1 1[ 1[ 1 = − U̇ · I𝜌 · U̇ X1 =X1 𝛿 Xi1+1 + U̇ · I𝜌 · U̇ X1 =X1 𝛿 Xi1 + U̇ · I𝜌 · U̇ dX . 2 2 2 ∫X 1 li i+ 1 i 1

i

The above result allows us to rewrite the time integral of the variation of the kinetic energy with respect to the variable domain boundaries (A.2) as t1

∫t

([

t1

𝛿l Ti dt =

∫t

0

0

t1

1 2 ∫t0

+

∫t

( [

U̇ · I𝜌 · U̇

X1

t1

−

i+ 1

∫X 1

0

]

𝜕U U̇ · I𝜌 · 𝜕X1

i

[

X 1 =X 1

i+ 1

]

[

𝛿 Xi1+1 − U̇ · I𝜌 · U̇

X 1 =X 1

i+ 1

̈ · I𝜌 · 𝜕 U U 𝜕X1

{

𝛿 Xi1+1 (

X 1 =X 1

X 1 =X 1 i

)

Ẋ 1i 𝛿 Xi1 dt

i

)

]

Xi1 − X 1 −

li

)

]

𝜕U Ẋ 1i+1 𝛿 Xi1+1 − U̇ · I𝜌 · 𝜕X1

𝛿 Xi1 dt

(A.4)

𝛿 Xi1 (

Xi1+1 − X 1

li

)

} dX 1 dt.

Next, we substitute the local balance of linear momentum (13) and angular momentum (14), which, in terms of the generalized displacement U, can be written as

̈ = I𝜌 · U

𝜕r 𝜕F + 0 × f + F ext . 𝜕X1 𝜕X1

(A.5)

The contribution of the kinetic energy to the variation of the action integral 𝛿 l Si can then be expressed as t1

∫t0

t1

𝛿l Ti dt =

([

∫t0 t

1 1 2 ∫t0

+

t1

−

∫t

( [

U̇ · I𝜌 · U̇

X1

i+ 1

∫X 1

0

]

𝜕U U̇ · I𝜌 · 𝜕X1

i

{

[

X 1 =X 1

i+ 1

]

𝜕U Ẋ 1i+1 𝛿 Xi1+1 − U̇ · I𝜌 · 𝜕X1 [

𝛿 Xi1+1 − U̇ · I𝜌 · U̇

X 1 =X 1 i+ 1

𝜕r 𝜕F + 0 × f + F ext 𝜕X1 𝜕X1

}

𝜕U · 1 𝜕X

)

] X 1 =X 1

X 1 =X 1 i

{

dt

i

)

]

Ẋ 1i 𝛿 Xi1

𝛿 Xi1 dt

𝛿 Xi1+1 ( li

Xi1

(A.6)

−X

1

)

−

𝛿 Xi1 (

Xi1+1

li

−X

1

)

} dX 1 dt.

Eventually, we perform integration by parts on the term containing the external forces, X1

−

i+ 1

∫X1

F

ext

i

𝜕U · 1 𝜕X

X1

−

i+ 1

∫X 1 i

F

ext

·U

{

𝛿 Xi1+1 ( li

𝛿 Xi1+1 − 𝛿 Xi1 li

by which we obtain the relation

Xi1

−X

1

)

− X1

1

dX +

𝛿 Xi1 (

i+ 1

∫X1 i

li

Xi1+1

𝜕 F ext ·U 𝜕X1

−X

1

} )

[

dX 1 = F ext · U

]

[

X 1 =X 1

i+ 1

{

𝛿 Xi1+1 ( li

Xi1

−X

1

)

−

𝛿 Xi1 ( li

𝛿 Xi1+1 − F ext · U

Xi1+1

−X

1

)

] X 1 =X 1 i

𝛿 Xi1

} dX 1 ,

(A.7)

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341 t1

∫t

([

t1

𝛿l Ti dt =

∫t

0

0

t1

1 2 ∫t0

+

(

[

U̇ · I𝜌 · U̇

X1

t1

−

i+ 1

(

∫t0 ∫X1 i

∫t

i+ 1

∫X1

0

F

[

𝜕U Ẋ 1i+1 𝛿 Xi1+1 − U̇ · I𝜌 · 𝜕X1

X 1 =X 1

i+ 1

]

[

X 1 =X 1 i+ 1

𝛿 Xi1+1 − U̇ · I𝜌 · U̇

𝜕r 𝜕F + 0 ×f 𝜕X1 𝜕X1

X1

t1

−

]

𝜕U U̇ · I𝜌 · 𝜕X1

ext

·U

)

𝛿 Xi1+1 − 𝛿 Xi1 li

i

𝜕U 𝜕X1

·

{

X 1 =X 1 i

∫t

Ẋ 1i 𝛿 Xi1 dt

i

)

X1

∫X 1

0

i

[

∫t0

]

F ext · U

} )

Xi1+1 − X 1

li

{

𝛿 Xi1+1 (

Xi1

li

[

X 1 =X 1 i+ 1

𝛿 Xi1 (

𝜕 F ext ·U 𝜕X1

i+ 1

(

t1

𝛿 Xi1 dt +

Xi1 − X 1 −

t1

dX dt +

X 1 =X 1

𝛿 Xi1+1 ( li

1

)

] )

]

35

−X

1

𝛿 Xi1+1 − F ext · U

)

] X 1 =X 1 i

𝛿 Xi1 dt

dX 1 dt

)

−

𝛿 Xi1 (

Xi1+1

li

−X

1

} )

dX 1 dt. (A.8)

Next, we recall the variations of the work performed by the external forces with respect to the domain boundaries (102) and we use the definition (104) along with the coordinate transformation (20) to rewrite it in terms of the material coordinate X1 as 𝜉i1+1

𝛿l Wiext =

∫𝜉 1

𝛿 Xi1+1 − 𝛿 Xi1

̃ ̃ F ext · U

li

i

X1 i+ 1

=

∫X 1

F

ext

𝛿 Xi1+1

·U

∫𝜉 1 i

−

𝛿 Xi1

1

dX −

li

i

𝜉i1+1

Ji−1 d𝜉 1 +

X1 i+ 1

𝜕̃ F ext ̃ · U 𝛿 li Ji−1 d𝜉 1 𝜕 li

𝜕 F ext ·U 𝜕X1

∫X1 i

{

𝛿 Xi1+1 (

Xi1

li

−X

1

)

+

𝛿 Xi1 (

Xi1+1

li

−X

1

} )

(A.9) dX . 1

Note that we do not consider discontinuities and concentrated loads in what follows for the sake of brevity. The time integral of the above relation enters the variation of the action 𝛿 l Si , where we can already note that the (spatial) integrals involving the external forces F ext do cancel upon summation:

(

t1

)

𝛿l Ti + 𝛿l Wiext dt =

∫t0

t

1 1 2 ∫t0

+

∫t

0

[

U̇ · I𝜌 · U̇

X1

t1

−

(

i+ 1

(

∫X 1 i

t1

([

𝜕U U̇ · I𝜌 · 𝜕X1

∫t0 ]

]

[

X 1 =X 1

i+ 1

[

X 1 =X 1 i+ 1

𝛿 Xi1+1 − U̇ · I𝜌 · U̇

𝜕r 𝜕F + 0 ×f 𝜕X1 𝜕X1

)

{

𝜕U 𝜕X1

·

𝜕U Ẋ 1i+1 𝛿 Xi1+1 − U̇ · I𝜌 · 𝜕X1 )

] X 1 =X 1 i

𝛿 Xi1 dt +

𝛿 Xi1+1 (

)

Xi1 − X 1 −

li

t1

( [

∫t0

)

] X 1 =X 1

Ẋ 1i 𝛿 Xi1

F ext · U

]

[

X 1 =X 1

i+ 1

𝛿 Xi1 (

Xi1+1 − X 1

li

dt

i

)

𝛿 Xi1+1 − F ext · U

)

] X 1 =X 1 i

𝛿 Xi1 dt.

} dX 1 dt . (A.10)

As with the external forces, we try to recover the variation of the work performed by the internal forces in what follows. For this purpose, we consider the contribution from the internal forces in the above relation, for which integration by parts on the first term gives X1

−

i+ 1

∫X 1 i

𝜕F 𝜕U · 𝜕X1 𝜕X1 (

X1

i+ 1

+

∫X 1

F·

i

{

𝛿 Xi1+1 (

Xi1

li

𝜕t 𝜕𝜞 + 1 𝜕X1 𝜕X1

){

−X

𝛿 Xi1+1 li

1

)

−

𝛿 Xi1 (

Xi1+1

li

(

)

Xi1 − X 1 −

−X

𝛿 Xi1

1

} )

[ dX 1 = F ·

(

li

Xi1+1 − X 1

)

𝜕U 𝜕X1

}

]

[

X 1 =X 1

i+ 1

X1

dX 1 −

𝜕U 𝜕X1

𝛿 Xi1+1 − F ·

i+ 1

∫X1

F·

i

𝜕U 𝜕X1

𝛿 Xi1+1

− li

]

𝛿 Xi1

X 1 =X 1

𝛿 Xi1

i

dX 1 .

(A.11)

In the above relation, we have used the definition of the generalized strain (7) to substitute the derivative of the generalized displacement by

𝜕t 𝜕U = 𝜞 − e1 + 11 . 𝜕X1 𝜕X

(A.12)

We perform another integration by parts on the first (spatial) integral on the right-hand side of (A.11), by which we obtain the relation X1

i+ 1

∫X 1 i

{

𝜕𝜞 F· 𝜕X1

𝛿 Xi1+1 ( li

X1

Xi1

−X

i+ 1 1 𝜕D − 𝜞 · 1 ·𝜞 2 ∫X 1 𝜕X i

1

)

{

−

𝛿 Xi1 (

Xi1+1

li

𝛿 Xi1+1 ( li

−X )

1

Xi1 − X 1 −

} )

dX 1 = −

𝛿 Xi1 ( li

] ] 1[ 1[ F · 𝜞 X1 =X1 𝛿 Xi1+1 + F · 𝜞 X1 =X1 𝛿 Xi1 2 2 i+ 1 i

Xi1+1 − X 1

} )

X 𝛿 Xi1+1 − 𝛿 Xi1 1 i+ 1 1 dX 1 + F·𝜞 dX . 2 ∫X 1 li 1

i

(A.13)

36

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

We further consider the spatial derivative of the stiffness tensor, for which we obtain a relation analogous to its variation (80):

) ) 𝜕D 𝜕𝜃 ( 𝜕𝜃 ( = EA 1 t2 ⊗ t 1 + t 1 ⊗ t2 − kGA 1 t 1 ⊗ t2 + t 2 ⊗ t1 , 𝜕X1 𝜕X 𝜕X

(A.14)

and, therefore, 1 2

− 𝜞·

𝜕t 𝜕D 𝜕𝜃 𝜕 r 𝜕𝜃 𝜕 r 𝜕𝜃 𝜕𝜃 · 𝜞 + F · 11 = −N 1 01 · t 2 + Q 1 01 · t 1 − Q +Q = 𝜕X1 𝜕X 𝜕X 𝜕X 𝜕X 𝜕X 𝜕X1 𝜕X1

(

𝜕 r0 ×f 𝜕X1

)

·

𝜕U . 𝜕X1

(A.15)

Substituting the results (A.11), (A.13) and (A.15) in relation (A.10) gives t1

∫t

(

)

𝛿l Ti + 𝛿l Wiext dt =

0

(

t

+

1 1 2 ∫t0

t1

+

𝜕U 𝜕X1

F·

∫t0 ∫t

U̇ · I𝜌 · U̇

([

X1

t1

−

[

i+ 1

∫X 1

0

F·

i

([

t1

∫t

𝜕U U̇ · I𝜌 · 𝜕X1

0

]

]

[

X 1 =X 1

i+ 1

[

X 1 =X 1 i+ 1

𝛿 Xi1+1 − U̇ · I𝜌 · U̇

]

[

𝛿 Xi1+1 − F ·

X 1 =X 1

i+ 1

𝜕U 𝜕X1

𝜕U Ẋ 1i+1 𝛿 Xi1+1 − U̇ · I𝜌 · 𝜕X1 )

] X 1 =X 1 i

𝛿 Xi1 dt + )

] X 1 =X 1

𝛿 Xi1 dt +

i

(

t1

∫t

[

)

] X 1 =X 1

F ext · U

(

1 2 ∫t0

dt

i

]

[

X 1 =X 1 i+ 1

0

t1

Ẋ 1i 𝛿 Xi1

𝛿 Xi1+1 − F ext · U

)

] X 1 =X 1 i

𝛿 Xi1 dt

) [ ] [ ] − F · 𝜞 X1 =X1 𝛿 Xi1+1 + F · 𝜞 X1 =X1 𝛿 Xi1 dt i+ 1

i

1 1 t1 X1 𝛿 Xi1+1 − 𝛿 Xi1 1 i+ 1 𝜕 U 𝛿 Xi+1 − 𝛿 Xi 1 1 dX dt + F · 𝜞 dX dt. 𝜕X1 li 2 ∫t0 ∫X1 li i

(A.16) From (101), we obtain the variations of the work of the internal forces with respect to the domain boundaries formulated in the material coordinate X1 as X1

𝛿l Wiint

i+ 1

=

(

∫X 1 i

𝜕U 1 −F · 1 + F · 𝜞 𝜕X 2

)

𝛿 Xi1+1 − 𝛿 Xi1 li

dX 1 .

(A.17)

Once we form the variation of the action Si of the i-th domain  i with respect to its variable boundaries from (A.16) and (A.17), we obtain the result

𝛿l S i =

(

t1

𝛿l T i −

∫t0 t1

+

1 2 ∫t0 t1

+

∫t

(

[

𝛿l Wiint

U̇ · I𝜌 · U̇

([ F·

0

𝜕U 𝜕X1

+

𝛿l Wiext

]

)

t1

dt =

([

𝜕U U̇ · I𝜌 · 𝜕X1

∫t0

[

X 1 =X 1 i+ 1

]

𝛿 Xi1+1 − U̇ · I𝜌 · U̇ [

X 1 =X 1

i+ 1

𝛿 Xi1+1 − F ·

𝜕U 𝜕X1

] X 1 =X 1

)

] X 1 =X 1 i

X 1 =X 1 i

i+ 1

𝛿 Xi1 dt + )

]

[

𝛿 Xi1 dt +

𝜕U Ẋ 1i+1 𝛿 Xi1+1 − U̇ · I𝜌 · 𝜕X1 t1

( [

∫t0 t

1 1 ∫ 2 t0

(

F ext · U

]

)

] X 1 =X 1 i

[

X 1 =X 1 i+ 1

Ẋ 1i 𝛿 Xi1 dt

𝛿 Xi1+1 − F ext · U

)

] X 1 =X 1 i

𝛿 Xi1 dt

) [ ] [ ] − F · 𝜞 X1 =X1 𝛿 Xi1+1 + F · 𝜞 X1 =X1 𝛿 Xi1 dt. i+ 1

i

(A.18) Note that the above result does not contain any integrals over the spatial domain  i any longer. For this reason, contributions from neighboring domains  i and  i+1 do cancel upon summation. At the terminal points of the (sub-)structure, either the material points are fixed with respect to the stretched coordinate, or their motion is kinematically prescribed, i.e., 𝛿 X11 = 𝛿 Xn1+1 = 0. As a consequence, variations of the action integral with respect to the boundaries of the variable domains do vanish, i.e.,

𝛿l S =

n ∑

𝛿l S i = 0 .

(A.19)

i=1

We can therefore conclude that the terms related to the coordinate transformation (20) and the generalization of the variation (62) leave the physics of the problem unaltered, as we would expect when we consider the change of coordinates as a purely mathematical exercise. References [1] G.F. Carrier, The spaghetti problem, Am. Math. Mon. 56 (10) (1949) 669, https://doi.org/10.2307/2305560, http://www.jstor.org/stable/2305560? origincrossref. [2] L. Vu-Quoc, S. Li, Dynamics of sliding geometrically-exact beams: large angle maneuver and parametric resonance, Comput. Methods Appl. Mech. Eng. 120 (12) (1995) 65–118, https://doi.org/10.1016/0045-7825(94)00051-N.

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

37

[3] B. Tabarrok, C. Leech, Y. Kim, On the dynamics of an axially moving beam, J. Franklin Inst. 297 (3) (1974) 201–220, https://doi.org/10.1016/00160032(74)90104-5. [4] M. Krupa, W. Poth, M. Schagerl, A. Steindl, W. Steiner, H. Troger, G. Wiedermann, Modelling, dynamics and control of tethered satellite systems, Nonlinear Dynam. 43 (12) (2006) 73–96, https://doi.org/10.1007/s11071-006-0752-z. [5] A. Pechstein, J. Gerstmayr, A Lagrange–Eulerian formulation of an axially moving beam based on the absolute nodal coordinate formulation, Multibody Syst. Dyn. 30 (3) (2013) 343–358, https://doi.org/10.1007/s11044-013-9350-2. [6] E. Oborin, Y. Vetyukov, I. Steinbrecher, Eulerian description of non-stationary motion of an idealized belt-pulley system with dry friction, Int. J. Solid Struct. 147 (2018) 40–51, https://doi.org/10.1016/j.ijsolstr.2018.04.007. [7] E. Oborin, Y. Vetyukov, Steady state motion of a shear deformable beam in contact with a traveling surface, Acta Mech. (2019), https://doi.org/10.1007/ s00707-019-02476-x. [8] K. Grundl, T. Schindler, H. Ulbrich, D.J. Rixen, ALE beam using reference dynamics, Multibody Syst. Dyn. 46 (2) (2019) 127–146, https://doi.org/10.1007/ s11044-019-09671-7. [9] J.-P. Liu, Z.-B. Cheng, G.-X. Ren, An Arbitrary Lagrangian–Eulerian formulation of a geometrically exact Timoshenko beam running through a tube, Acta Mech. 229 (8) (2018) 3161–3188, https://doi.org/10.1007/s00707-018-2161-z. [10] C.D. Mote Jr., A study of band saw vibrations, J. Franklin Inst. 279 (6) (1965) 430–444, https://doi.org/10.1016/0016-0032(65)90273-5. [11] Y. Vetyukov, P.G. Gruber, M. Krommer, J. Gerstmayr, I. Gafur, G. Winter, Mixed Eulerian-Lagrangian description in materials processing: deformation of a metal sheet in a rolling mill, Int. J. Numer. Methods Eng. 109 (10) (2017) 1371–1390, https://doi.org/10.1002/nme.5314. [12] D. Hochlenert, G. Spelsberg-Korspeter, P. Hagedorn, Friction induced vibrations in moving continua and their application to brake squeal, J. Appl. Mech. 74 (3) (2007) 542–549, https://doi.org/10.1115/1.2424239, https://asmedigitalcollection.asme.org/appliedmechanics/article/74/3/542/476555/FrictionInduced-Vibrations-in-Moving-Continua-and. [13] G. Spelsberg-Korspeter, O.N. Kirillov, P. Hagedorn, Modeling and stability analysis of an axially moving beam with frictional contact, J. Appl. Mech. 75 (3) (2012), https://doi.org/10.1115/1.2755166, https://asmedigitalcollection.asme.org/appliedmechanics/article/doi/10.1115/1.2755166/465668/Modelingand-Stability-Analysis-of-an-Axially. [14] L.-Q. Chen, Analysis and control of transverse vibrations of axially moving strings, Appl. Mech. Rev. 58 (2) (2005) 91–116, https://doi.org/10.1115/1. 1849169, https://asmedigitalcollection.asme.org/appliedmechanicsreviews/article/58/2/91/446339/Analysis-and-Control-of-Transverse-Vibrations-of. [15] K. Marynowski, T. Kapitaniak, Dynamics of axially moving continua, Int. J. Mech. Sci. 81 (2014) 26–41, https://doi.org/10.1016/j.ijmecsci.2014.01.017. [16] J. Aitken, XIV. An account of some experiments on rigidity produced by centrifugal force, Lond. Edinb. Dublin Philos. Mag. J. Sci. 5 (29) (1878) 81–105, https://doi.org/10.1080/14786447808639394. [17] R. Skutsch, Ueber die Bewegung eines gespannten Fadens, welcher gezwungen ist, durch zwei feste Punkte mit einer constanten Geschwindigkeit zu gehen, und zwischen denselben in Transversalschwingungen von geringer Amplitude versetzt wird (On the motion and small amplitude vibrations of a tensioned string, which is forced to move through two fixed points with a constant velocity), Ann. Phys. Chem. 297 (5) (1897) 190–195, https://doi.org/10. 1002/andp.18972970510. [18] R.A. Sack, Transverse oscillations in travelling strings, Br. J. Appl. Phys. 5 (6) (1954) 224–226, https://doi.org/10.1088/0508-3443/5/6/307, http://stacks.iop. org/0508-3443/5/i6/a307?keycrossref.4ec16af690677c8e6299165b6beddf6d. [19] R.D. Swope, W. Ames, Vibrations of a moving threadline, J. Franklin Inst. 275 (1) (1963) 36–55, https://doi.org/10.1016/0016-0032(63)90619-7. [20] C.D. Mote Jr., Dynamic stability of an axially moving band, J. Franklin Inst. 285 (5) (1968) 329–346, https://doi.org/10.1016/0016-0032(68)90482-1. [21] K.W. Buffinton, T.R. Kane, Dynamics of a beam moving over supports, Int. J. Solid Struct. 21 (7) (1985) 617–643, https://doi.org/10.1016/00207683(85)90069-1. [22] J.A. Wickert, C.D. Mote Jr., Classical vibration analysis of axially moving continua, J. Appl. Mech. Trans. ASME 57 (3) (1990) 738–744, https://doi.org/10. 1115/1.2897085. [23] J.A. Wickert, Non-linear vibration of a traveling tensioned beam, Int. J. Non Lin. Mech. 27 (3) (1992) 503–517, https://doi.org/10.1016/00207462(92)90016-Z. [24] K. Behdinan, M.C. Stylianou, B. Tabarrok, Dynamics of flexible sliding beams - non-linear analysis part I: Formulation, J. Sound Vib. 208 (4) (1997) 517–539, https://doi.org/10.1006/jsvi.1997.1167. [25] K. Behdinan, B. Tabarrok, Dynamics of flexible sliding beams - non-linear analysis part II: transient response, J. Sound Vib. 208 (4) (1997) 541–565, https:// doi.org/10.1006/jsvi.1997.1168. [26] D.B. McIver, Hamilton’s principle for systems of changing mass, J. Eng. Math. 7 (3) (1973) 249–261, https://doi.org/10.1007/BF01535286. [27] E. Reissner, On one-dimensional finite-strain beam theory: the plane problem, Z. Angew. Math. Phys. (1972), https://doi.org/10.1007/BF01602645. [28] S.S. Antman, The theory of rods, in: C. Truesdell (Ed.), Linear Theories of Elasticity and Thermoelasticity, VIa/2, Springer Berlin Heidelberg, Berlin, Heidelberg, 1973, pp. 641–703, https://doi.org/10.1007/978-3-662-39776-3_6. [29] J.C. Simo, L. Vu-Quoc, On the dynamics of flexible beams under large overall motionsthe plane case: Part I, J. Appl. Mech. 53 (4) (1986) 849, https://doi.org/ 10.1115/1.3171870, http://appliedmechanics.asmedigitalcollection.asme.org/article.aspx?articleid1408677. [30] J.C. Simo, L. Vu-Quoc, On the dynamics of flexible beams under large overall motionsthe plane case: Part II, J. Appl. Mech. 53 (4) (1986) 855, https://doi. org/10.1115/1.3171871, http://appliedmechanics.asmedigitalcollection.asme.org/article.aspx?articleid1408679. [31] H. Irschik, H.J. Holl, The equations of Lagrange written for a non-material volume, Acta Mech. 153 (34) (2002) 231–248, https://doi.org/10.1007/ BF01177454. [32] Y. Vetyukov, Non-material finite element modelling of large vibrations of axially moving strings and beams, J. Sound Vib. 414 (2018) 299–317, https://doi. org/10.1016/j.jsv.2017.11.010. [33] R.L. Lowe, C.G. Cooley, A Newtonian mechanics formulation for the vibration of translating and rotating elastic continua, J. Vib. Contr. 25 (10) (2019) 1639–1652, https://doi.org/10.1177/1077546319825675. [34] Y. Vetyukov, E. Oborin, M. Krommer, V. Eliseev, Transient modelling of flexible belt drive dynamics using the equations of a deformable string with discontinuities, Math. Comput. Model. Dyn. Syst. 23 (1) (2017) 40–54, https://doi.org/10.1080/13873954.2016.1232281. [35] M. Olsson, On the fundamental moving load problem, J. Sound Vib. (1991), https://doi.org/10.1016/0022-460X(91)90593-9. [36] L. Vu-Quoc, M. Olsson, A computational procedure for interaction of high-speed vehicles on flexible structures without assuming known vehicle nominal motion, Comput. Methods Appl. Mech. Eng. 76 (3) (1989) 207–244, https://doi.org/10.1016/0045-7825(89)90058-3. [37] L. Vu-Quoc, M. Olsson, New predictor/corrector algorithms with improved energy balance for a recent formulation of dynamic vehicle/structure interaction, Int. J. Numer. Methods Eng. 32 (2) (1991) 223–253, https://doi.org/10.1002/nme.1620320202. [38] E.C. Cojocaru, H. Irschik, K. Schlacher, Concentrations of pressure between an elastically supported beam and a moving Timoshenko-beam, J. Eng. Mech. 129 (9) (2003) 1076–1082, https://doi.org/10.1061/(ASCE)0733-9399(2003)129:9(1076). [39] E.C. Cojocaru, H. Irschik, H. Gattringer, Dynamic response of an elastic bridge due to a moving elastic beam, Comput. Struct. 82 (1112) (2004) 931–943, https://doi.org/10.1016/j.compstruc.2004.02.001. [40] A. Humer, Dynamic modeling of beams with non-material, deformation-dependent boundary conditions, J. Sound Vib. 332 (3) (2013) 622–641, https:// doi.org/10.1016/j.jsv.2012.08.026. [41] I. Steinbrecher, A. Humer, L. Vu-Quoc, On the numerical modeling of sliding beams: a comparison of different approaches, J. Sound Vib. 408 (2017) 270–290, https://doi.org/10.1016/j.jsv.2017.07.010. [42] D. Hong, G. Ren, A modeling of sliding joint on one-dimensional flexible medium, Multibody Syst. Dyn. 26 (1) (2011) 91–106, https://doi.org/10.1007/ s11044-010-9242-7. [43] S. Yang, Z. Deng, J. Sun, Y. Zhao, S. Jiang, A variable-length beam element incorporating the effect of spinning, Lat. Am. J. Solid. Struct. 14 (8) (2017) 1506–1528, https://doi.org/10.1590/1679-78253894.

38

A. Humer et al. / Journal of Sound and Vibration 480 (2020) 115341

[44] A. Humer, Elliptic integral solution of the extensible elastica with a variable length under a concentrated force, Acta Mech. 222 (34) (2011) 209–223, https://doi.org/10.1007/s00707-011-0520-0. [45] A. Humer, Exact solutions for the buckling and postbuckling of shear-deformable beams, Acta Mech. 224 (7) (2013) 1493–1525, https://doi.org/10.1007/ s00707-013-0818-1. [46] A. Humer, A.S. Pechstein, Exact solutions for the buckling and postbuckling of a shear-deformable cantilever subjected to a follower force, Acta Mech. (2019), https://doi.org/10.1007/s00707-019-02472-1. [47] S. Timoshenko, LXVI. On the correction for shear of the differential equation for transverse vibrations of prismatic bars, Lond. Edinb. Dublin Philos. Mag. J. Sci. 41 (245) (1921) 744–746, https://doi.org/10.1080/14786442108636264. [48] Dassault Systmes Simulia Corp, Abaqus 6.14 Documentation, 2015. Providence, RI. [49] I.-S. Liu, Continuum Mechanics, Advanced Texts in Physics, Springer Berlin Heidelberg, Berlin, Heidelberg, 2002, https://doi.org/10.1007/978-3-66205056-9. [50] D. Griffiths, S. Walborn, Dirac deltas and discontinuous functions, Am. J. Phys. 67 (5) (1999) 446–447, https://doi.org/10.1119/1.19283. [51] F. Coutinho, Y. Nogami, F. Toyama, Unusual situations that arise with the Dirac delta function and its derivative, Rev. Bras. Ensino Fsica 31 (4) (2009) 4302–4308, https://doi.org/10.1590/S1806-11172009000400004. [52] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II, Vol. 14 of Springer Series in Computational Mathematics, Springer Berlin Heidelberg, Berlin, Heidelberg, 1996https://doi.org/10.1007/978-3-642-05221-7. [53] J. Gerstmayr, A. Dorninger, R. Eder, P. Gruber, D. Reischl, M. Saxinger, M. Schörgenhumer, A. Humer, K. Nachbagauer, A. Pechstein, Y. Vetyukov, HOTINT: a script language based framework for the simulation of multibody dynamics systems, in: 9th International Conference on Multibody Systems, Nonlinear Dynamics, and Control, 7B, ASME, 2013, p. 10, https://doi.org/10.1115/DETC2013-12299. [54] A. Humer, G. Jungmayr, W. Koppelstätter, M. Schörgenhumer, S. Silber, G. Weidenholzer, S. Weitzhofer, Multi-objective optimization of complex multibody systems by coupling HOTINT with MagOpt, in: 12th International Conference on Multibody Systems, Nonlinear Dynamics, and Control, vol. 6, ASME, 2016, p. 8, https://doi.org/10.1115/DETC2016-60204. [55] P. Wriggers, Nonlinear Finite Element Methods, Springer Berlin Heidelberg, Berlin, Heidelberg, 2008, https://doi.org/10.1007/978-3-540-71001-1.