Computational Dynamics of a 3D Elastic String Pendulum Attached to a Rigid Body and an Inertially Fixed Reel Mechanism
aa r X i v : . [ m a t h . D S ] S e p Nonlinear Dynamics manuscript No. (will be inserted by the editor)
Computational Dynamics of a 3D Elastic String PendulumAttached to a Rigid Body and an Inertially Fixed ReelMechanism
Taeyoung Lee · Melvin Leok · N. Harris McClamroch
September 10, 2009
Abstract
A high fidelity model is developed for anelastic string pendulum, one end of which is attachedto a rigid body while the other end is attached to aninertially fixed reel mechanism which allows the un-stretched length of the string to be dynamically varied.The string is assumed to have distributed mass andelasticity that permits axial deformations. The rigidbody is attached to the string at an arbitrary point, andthe resulting string pendulum system exhibits nontriv-ial coupling between the elastic wave propagation in thestring and the rigid body dynamics. Variational meth-ods are used to develop coupled ordinary and partialdifferential equations of motion. Computational meth-ods, referred to as Lie group variational integrators, arethen developed, based on a finite element approxima-tion and the use of variational methods in a discrete-time setting to obtain discrete-time equations of mo-tion. This approach preserves the geometry of the con-figurations, and leads to accurate and efficient algo-rithms that have guaranteed accuracy properties thatmake them suitable for many dynamic simulations, es-pecially over long simulation times. Numerical resultsare presented for typical examples involving a constantlength string, string deployment, and string retrieval.These demonstrate the complicated dynamics that arisein a string pendulum from the interaction of the rigid
Taeyoung Lee, Assistant ProfessorDepartment of Mechanical and Aerospace Engineering, FloridaInstitute of Technology, Melbourne, FL 32901. E-mail: taey-oung@fit.eduMelvin Leok, Associate Professor,Department of Mathematics, University of California, San Diego,CA 92093. E-mail: [email protected]. Harris McClamroch, Professor,Department of Aerospace Engineering, University of Michigan,Ann Arbor, MI 48109. E-mail: [email protected] body motion, elastic wave dynamics in the string, andthe disturbances introduced by the reeling mechanism.Such interactions are dynamically important in manyengineering problems, but tend be obscured in lowerfidelity models.
Keywords
Lagrangian mechanics · geometric inte-grator · variational integrator · string pendulum · reelmechanism · rigid body The dynamics of a body connected to a string appearin several engineering problems such as cable cranes,towed underwater vehicles, and tethered spacecraft. Sev-eral types of analytical and numerical models have beendeveloped. Lumped mass models, where the string isspatially discretized into connected point masses, aredeveloped in [1,2,3]. Finite difference methods in boththe spatial domain and the time domain are appliedin [4,5]. Finite element discretizations of the weak formof the equations of motion are used in [5,6]. Variable-length string models also have been developed: a vari-able length string is modeled based on a continuousplastic impact assumption in [7,8], and a reel mecha-nism is considered in [9,10]. But, the reel mechanismsdeveloped in those papers are problematic. In [9], thedeployed portion of the string is assumed to move alonga fixed line. The dynamic model of reeling developedin [10] is erroneous (this will be discussed further inSection 2). Instead of a point mass, a rigid body modelis considered in [8], but this paper does not provide anycomputational results. Analytical and numerical mod-els of a rigid body connected to an elastic string appearin [11].
The goal of this paper is to develop an analyticalmodel and a numerical algorithm that can be used forsimulation of an an elastic string attached to a rigidbody and an inertially fixed reel mechanism, acting un-der a constant gravitational potential. The string hasdistributed mass; it can move in a three-dimensionalspace while deforming axially; the rigid body attachedto the string can translate and rotate. We assume thatthe point where the string is attached to the rigid bodyis displaced from the center of mass of the rigid body sothat there exist nonlinear coupling effects between thestring deformation dynamics and the rigid body dy-namics. The reel mechanism is an inertially-fixed sys-tem consisting of a cylindrical reel on which the stringwinds and unwinds and a guide way that acts as a pivotfor the deployed portion of the string. The portion ofthe string on the reel mechanism is assumed to be inex-tensible. The combined system of the string, the rigidbody and the reel mechanism provides a realistic andaccurate dynamic model of cable cranes and towing sys-tems.In this paper, we first show that the governing equa-tions of motion of the presented string pendulum canbe developed according to Hamilton’s variational prin-ciple. The configuration manifold of the string pendu-lum is expressed as the product of the real space R representing the configuration of the reel mechanism,the space of connected curve segments on R describingthe deployed portion of the string, and the special or-thogonal group SO ( ) defining the attitude of the rigidbody [12]. The variational principle is carefully appliedto respect the geometry of the Lie group configurationmanifold. We incorporate an additional modificationterm, referred to as the Carnot energy loss term [13],in the variational principle to take account of the factthat the portion of the string in the reel mechanismis inextensible. The resulting Euler-Lagrange equationsare expressed as coupled partial and ordinary differen-tial equations.The second part of this paper deals with a geomet-ric numerical integrator for the model we presented forthe string pendulum. Geometric numerical integrationis concerned with developing numerical integrators thatpreserve geometric features of a system, such as in-variants, symmetry, and reversibility [14]. The stringpendulum is a Lagrangian/Hamiltonian system evolv-ing on a Lie group. When numerically simulating suchsystems, it is critical to preserve both the symplecticproperty of Hamiltonian flows and the Lie group struc-ture for numerical accuracy and efficiency [15]. A geo-metric numerical integrator, referred to as a Lie groupvariational integrator, has been developed for a Hamil-tonian system on an arbitrary Lie group and it has been applied to several multibody systems ranging from bi-nary asteroids to articulated rigid bodies and magneticsystems in [16,17].This paper develops a Lie group variational inte-grator for the proposed string pendulum model. Thisextends the results presented in [17] by incorporatingdeformation of the string using a finite element modeland by including a discrete-time Carnot energy lossterm. The proposed geometric numerical integrator pre-serves the symplectic structure and momentum maps,and exhibits desirable energy conservation properties.It also respects the Lie group structure of the configu-ration manifold, and avoids the singularities and com-putational complexities associated with the use of lo-cal coordinates, explicit constraints or projection. As aresult, this computational approach can represent arbi-trary translations and rotations of the rigid body andlarge deformations of the string.In summary, this paper develops an analytical modeland a geometric numerical integrator for a string pen-dulum attached to a rigid body and a reel mechanism.These provide a realistic mathematical model for teth-ered systems and a reliable numerical simulation toolthat characterizes the nonlinear coupling between thestring dynamics, the rigid body dynamics, and the reelmechanism accurately. The proposed high-fidelity com-putational framework can be naturally extended to for-mulating and solving control problems associated withstring deployment, retrieval, and vibration suppressionas in [18].This paper is organized as follows. A string pendu-lum is described and the corresponding Euler-Lagrangeequations are presented in Section 2. A Lie group vari-ational integrator is derived in Section 3, followed bynumerical examples and conclusions in Section 4 and 5. is connected to an inertially-fixed reel mechanism com-posed of a drum and a guide way. The string is woundaround the drum at a constant radius, and the stringon the drum and in the guide way is assumed to be in-extensible. A control moment is applied at the rotatingdrum. This system of the string, the rigid body, and thereel mechanism, acting under a constant gravitationalpotential, is referred to as a string pendulum. This isillustrated in Figure 1.We choose an inertially fixed reference frame and abody-fixed frame. The origin of the body-fixed frameis located at the end of the string where the string isattached to the rigid body, and it is fixed to the rigidbody. Since the string is extensible, we need to distin-guish between the arc length for the stretched deformedconfiguration and the arc length for the unstretched ref-erence configuration. Define r d ∈ R the location of the origin of the axisof the drum d ∈ R the radius of the drum I d = κ d d ∈ R the rotational inertia of the drumfor κ d ∈ R b ∈ R the length from the drum to theguide way u ∈ R the control moment applied at drum L ∈ R the total unstretched length of thestring µ ∈ R the mass of the string per the unitunstretched length O the point at which the string is at-tached to the drum s ∈ [0 , L ] the unstretched arc length of thestring between the point O and amaterial point P on the string s ( s, t ) ∈ R the stretched arc length to a mate-rial point Ps p ( t ) ∈ [ b, L ] the arc length of the string betweenthe point O and the material pointon the string located at the guideway entrance r ( s, t ) ∈ R the deformed location of a materialpoint Pθ ( s ) ∈ R θ = (( s p − b ) − s ) /d for s ∈ [0 , s p − b ] M ∈ R the mass of the rigid body J ∈ R the inertia matrix of the rigid bodywith respect to the body fixed frame ρ c ∈ R the vector from the origin of thebody fixed frame to the center ofmass of the rigid body representedin the body fixed frame R ∈ SO ( ) the rotation matrix from the bodyfixed frame to the reference frame Ω ∈ R the angular velocity of the rigid bodyrepresented in the body fixed frameThe configuration of the string on the drum and inthe guide way is completely determined by the variable s p ( t ), since the string there is inextensible. The con-figurations of the deployed portion of the string andthe rigid body are described by the curve r ( s, t ) for s ∈ [ s p , L ], and the rotation matrix R ∈ SO ( ), respec-tively, where the special orthogonal group is SO ( ) = { R ∈ R × | R T R = I, det[ R ] = 1 } . Therefore, the con-figuration manifold of the string pendulum is the prod-uct of the real space R , the space of connected curveson R , and the special orthogonal group SO ( ).The attitude kinematics equation of the rigid bodyis given by˙ R = R ˆ Ω, (1)where the hat map ˆ · : R → so (3) is defined by thecondition that ˆ xy = x × y for any x, y ∈ R . Since ˆ x is a 3 × x T = − ˆ x .The inverse map of the hat map is referred to as the vee map : ( · ) ∨ : so (3) → R .2.2 LagrangianWe develop Euler-Lagrange equations for the string pen-dulum according to Hamilton’s variational principle.The Lagrangian of the string pendulum is derived, andthe corresponding action integral is defined. Due to theunique dynamic characteristics of the string pendulum,the variation of the action integral should be carefullydeveloped: (i) since the unstretched length of the de-ployed portion of the string is not fixed, when deriv-ing the variation of the corresponding part of the ac-tion integral, we need to apply Green’s theorem; (ii)since the attitude of the rigid body is represented inthe special orthogonal group, the variation of rotationmatrices are carefully expressed by using the exponen-tial map [17,16]; (iii) since the portion of the stringon the guide way and the drum is inextensible, the ve-locity of the string is not continuous at the guide wayentrance. To take account of the effect of this velocitydiscontinuity, an additional modification term, referredto as a Carnot energy loss term is incorporated [13].Then, Euler-Lagrange equations are derived accordingto Hamilton’s principle, and they are expressed as cou-pled ordinary and partial differential equations. Lagrangian
The total kinetic energy is composed of thekinetic energy of the portion of the string on the drumand the guide way T r , the kinetic energy of the deployed e e s P r d r p bρ c d O global reference framebody fixed frame(a) Reference configuration e e r ( s, t ) Pr ( s p ( t ) , t ) s = 0 O s = s p ( t ) s = L ρ c R ( t ) θ (0)(b) Deformed configuration Fig. 1
String Pendulum Model portion of the string T s , and the kinetic energy of therigid body T b . The kinetic energy T r can be written as T r = Z s p µ ˙ r ( s ) · ˙ r ( s ) ds + 12 I d ˙ θ (0) , where the dot represents the partial derivative with re-spect to time. Here, the dependency of variables on time t is omitted for simplicity, i.e. r ( s ) = r ( s, t ). The veloc-ity of the string in the reel mechanism is equal to ˙ s p asthe string is inextensible. From the definitions, we have˙ θ (0) = ˙ s p /d , and I d = κ d d . Then, the kinetic energy T r can be written as T r = 12 ( µs p + κ d ) ˙ s p . (2)The kinetic energy of the deployed portion of the stringis given by T s = Z Ls p µ ˙ r ( s ) · ˙ r ( s ) ds. (3)Let ρ ∈ R be the vector from the free end of the string r ( L ) to a mass element of the rigid body, expressed inthe body fixed frame. The location of the mass element in the reference frame is given by r ( L ) + Rρ . Then, thekinetic energy of the rigid body is given by T b = Z body k ˙ r ( L ) + R ˆ Ωρ k dm = 12 M ˙ r ( L ) · ˙ r ( L ) + M ˙ r ( L ) · R ˆ Ωρ c + 12 Ω · JΩ, (4)where J = R − ˆ ρ dm is the inertia matrix of the rigidbody in the body fixed frame.Now we obtain expressions for the potential energyof each part. The gravitational potential energy of theportion of the string on the drum and the guide way isgiven by V r = − Z s p − b µg ( r d · e − d sin θ ) ds = − µg (cid:8) ( s p − d ) r d · e + d (cos(( s p − b ) /d ) − (cid:9) . (5)The strain of the string at a material point located at r ( s ) is given by ǫ = lim ∆s → ∆s ( s ) − ∆s∆s = s ′ ( s ) − , where ( ) ′ denote the partial derivative with respect to s . The tangent vector at the material point is given by e t = ∂r ( s ) ∂s = ∂r ( s ) ∂s ∂s∂s ( s ) = r ′ ( s ) s ′ ( s ) . Since this tangent vector has unit length, we have s ′ ( s ) = k r ′ ( s ) k . Therefore, the strain of the string is given by ǫ = k r ′ ( s ) k −
1. The potential energy of the deployedportion of the string is composed of the elastic potentialand the gravitational potential energy: V s = Z Ls p EA ( k r ′ ( s ) k − − µgr ( s ) · e ds, (6)where E and A denote the Young’s modulus and thecross sectional area of the string, respectively. The grav-itational potential of the rigid body is given by V b = − M g ( r ( L ) + Rρ c ) · e . (7)In summary, the Lagrangian of the string pendulumis given by L = ( T r − V r ) + ( T s − V s ) + ( T b − V b ) = L r + L s + L b . (8)2.3 Variational Approach Action Integral
The action integral is defined by G = Z t f t L r + L s + L b dt = G r + G s + G b . (9)We find expressions for the variation of each term ofthe action integral. Variation of G r From (2) and (5), the variation of G r is given by δ G r = Z t f t (cid:26) − ( µs p + κ d )¨ s p − µ ˙ s p + µg ( r d · e ) − µgd sin(( s p − b ) /d ) (cid:27) δs p dt, (10)where we used integration by parts. Variation of G s From (3), (6), (9), the second term ofthe action integral G s is a double integral on ( t, s ) ∈ [ t , t f ] × [ s p ( t ) , L ]. Since the variable s p ( t ) is dependenton the time t , the variation of G s should take into ac-count the variation of s p ( t ): δ G s = Z t f t Z Ls p ( t ) µ ˙ r ( s ) · δ ˙ r ( s ) − EA k r ′ ( s ) k − k r ′ ( s ) k r ′ ( s ) · δr ′ ( s ) + µge · δr ( s ) ds dt − Z t f t (cid:26) µ ˙ r ( s + p ) · ˙ r ( s + p ) − EA ( (cid:13)(cid:13) r ′ ( s + p ) (cid:13)(cid:13) − + µgr ( s p ) · e (cid:27) δs p dt, (11)where r ( s + p ) represents the material point of the stringlocated just outside the guide way.Now we focus on the first term of (11). Here, wecannot apply integration by parts at time t , since theorder of the integrals in (11) cannot be interchanged dueto the time dependence in the variable s p ( t ). Instead,we use Green’s theorem, I B ˙ r ( s ) · δr ( s ) ds = Z t f t Z Ls p ( t ) ddt ( ˙ r ( s ) · δr ( s )) dsdt, (12)where H B represents the counterclockwise line integralon the boundary B of the region [ t , t f ] × [ s p ( t ) , L ]. Theboundary B is composed of four lines: ( t = t , s ∈ [ s p ( t ) , L ]), ( t = t f , s ∈ [ s p ( t f ) , L ]), ( t ∈ [ t , t f ] , s = s p ( t )), and ( t ∈ [ t , t f ] , s = L ). For the first two lines, δr ( s ) = 0 since t = t , t f . For the last line, ds = 0 since s is fixed. Thus, parameterizing the third line by t , weobtain I B ˙ r ( s ) · δr ( s ) ds = Z t f t ˙ r ( s p ( t )) · δr ( s p ( t )) ˙ s p ( t ) dt. Substituting this into (12) and rearranging, the firstterm of (11) is given by Z t f t Z Ls p ˙ r ( s ) · δ ˙ r ( s )= Z t f t "Z Ls p − ¨ r ( s ) · δr ( s ) ds + ˙ r ( s p ) · δr ( s p ) ˙ s p dt. (13)Substituting this into (11), and using integration byparts with respect to s for the second term of (11), the variation of G s can be written as δ G s = Z t f t Z Ls p {− µ ¨ r ( s ) + F ′ ( s ) + µge } · δr ( s ) ds dt + Z t f t (cid:26) − µ ˙ r ( s + p ) · ˙ r ( s + p ) + 12 EA ( (cid:13)(cid:13) r ′ ( s + p ) (cid:13)(cid:13) − − µgr ( s p ) · e (cid:27) δs p + µ ˙ r ( s + p ) · δr ( s + p ) ˙ s p − F ( L ) · δr ( L ) + F ( s p ) · δr ( s + p ) dt, where F ( s ) = EA k r ′ ( s ) k − k r ′ ( s ) k r ′ ( s ) represents the tensionof the string.We simplify this using the boundary condition atthe guide way. The location of the guide way entrance isgiven by r p = r ( s p ( t ) , t ). Since the location is inertiallyfixed, we have δr p = δr ( s + p ) + r ′ ( s + p ) δs p = 0, and ˙ r p =˙ r ( s + p ) + r ′ ( s + p ) ˙ s p = 0. Substituting these, we obtain δ G s = Z t f t Z Ls p ( t ) {− µ ¨ r ( s ) + F ′ ( s ) + µge } · δr ( s ) ds dt + Z t f t (cid:26) µ (cid:13)(cid:13) r ′ ( s + p ) (cid:13)(cid:13) ˙ s p + 12 EA ( (cid:13)(cid:13) r ′ ( s + p ) (cid:13)(cid:13) − − µgr ( s p ) · e (cid:27) δs p dt − Z t f t F ( s p ) · r ′ ( s + p ) δs p + F ( L ) · δr ( L ) dt. (14) Variation of G b From (4), (7), the variation of G b isgiven by δ G b = Z t f t { M ˙ r ( L ) + M R ˆ Ωρ c } · δ ˙ r ( L )+ { JΩ + M ˆ ρ c R T ˙ r ( L ) } · δΩ + M ˙ r ( L ) · δR ˆ Ωρ c + M ge · δr ( L ) + M ge · δRρ c dt. (15)The attitude of the rigid body is represented by therotation matrix R ∈ SO ( ). Therefore, the variationof the rotation matrix should be consistent with thegeometry of the special orthogonal group. In [17,16], itis expressed in terms of the exponential map as δR = ddǫ (cid:12)(cid:12)(cid:12)(cid:12) ǫ =0 R ǫ = ddǫ (cid:12)(cid:12)(cid:12)(cid:12) ǫ =0 R exp ǫ ˆ η = R ˆ η (16)for η ∈ R . The key idea is expressing the variationof a Lie group element in terms of a Lie algebra ele-ment. This is desirable since the Lie algebra so (3) ofthe special orthogonal group, represented by 3 × R using the vee map. As a result, the variation of thethree-dimenstional rotation matrix R is expressed in terms of a vector η ∈ R . We can directly show that(16) satisfies δ ( R T R ) = δR T R + R T δR = − ˆ η + ˆ η = 0.The corresponding variation of the angular velocity isobtained from the kinematics equation (1): δ ˆ Ω = ddǫ (cid:12)(cid:12)(cid:12)(cid:12) ǫ =0 ( R ǫ ) T ˙ R ǫ = ( ˙ η + Ω × η ) ∧ . (17)Substituting (16), (17) into (15), and using integra-tion by parts for η , we obtain δ G b = Z t f t M {− ¨ r ( L ) − R ˆ Ω ρ c − R ˆ˙ Ωρ c + ge } · r ( L )+ {− J ˙ Ω − M ˆ ρ c R T ¨ r ( L ) − ˆ ΩJΩ + M g ˆ ρ c R T e } · η dt. (18) Variation of G From (10), (14), (18), the variation ofthe action integral is given by δ G = δ G r + δ G s + δ G b . (19) Variational Principle with Discontinuity
Let r p = r ( s p ( t ) , t )be the location of the pivot in the reference frame. Sinceit is fixed, we have ˙ r p = ˙ r ( s p , t ) + r ′ ( s p , t ) ˙ s p = 0. Let r ( s − p ), and r ( s + p ) be the material point of the string justinside the guide way, and the material point just outsidethe guide way, respectively. Since the string is inextensi-ble inside the guide way, k r ′ ( s − p ) k = 1. Since the stringis extensible outside the guide way, k r ′ ( s + p ) k = 1 + ǫ + ,where ǫ + represents the strain of the string just outsidethe guide way. Using these, the speeds of the string atthose points are given by k ˙ r ( s − p ) k = k − r ′ ( s − p ) ˙ s p k = | ˙ s p |k ˙ r ( s + p ) k = k − r ′ ( s + p ) ˙ s p k = (1 + ǫ + ) | ˙ s p | . Therefore, the speed of the string changes instanta-neously by the amount ǫ + | ˙ s p | at the guide way.Due to this velocity and strain discontinuity, thevariation of the action integral is not equal to the neg-ative of the virtual work done by the external controlmoment u at the drum. In order to derive equations ofmotion using Hamilton’s principle, an additional term Q , referred to as Carnot energy loss term should be in-troduced [13,8]. The resulting variational principle isgiven by δ G + Z t f t ( Q + u/d ) δs p dt = 0 . (20)The corresponding time rate of change of the total en-ergy is given by ˙ E = ( Q + u/d ) ˙ s p , where the first term Q ˙ s p represents the energy dissipation rate due to thevelocity and strain discontinuity. Consider the infinitesimal mass element dm = µ ˙ s p dt located just outside the guide way. Without loss of gen-erality, we assume that ˙ s p > ǫ + ) ˙ s p can be considered as a plastic impact intothe the portion of the string on the guide way movingwith velocity ˙ s p . The corresponding energy dissipationrate is given by Q ˙ s p = − µ ( ǫ + ) ˙ s p − EA ( ǫ + ) ˙ s p . (21)(See [13,8]). Dividing both side by ˙ s p , we obtain theexpression for the Carnot energy loss term Q .2.4 Euler-Lagrange EquationsSubstituting (19), (21) into (20), we obtain Euler-Lagrangeequations for the string pendulum: − ( µs p + κ d )¨ s p + µg ( r d − r p ) · e − µgd sin(( s p − b ) /d ) − F ( s + p ) · r ′ ( s + p ) + µ ( (cid:13)(cid:13) r ′ ( s + p ) (cid:13)(cid:13) −
1) ˙ s p + ud = 0 , (22) − µ ¨ r ( s ) + F ′ ( s ) + µge = 0 , s ∈ [ s p , L ] , (23) − M ¨ r ( L ) − M R ˆ Ω ρ c − M R ˆ˙ Ωρ c + M ge − F ( L ) = 0 , (24) − J ˙ Ω − M ˆ ρ c R T ¨ r ( L ) − ˆ ΩJΩ + M g ˆ ρ c R T e = 0 , (25)where F ( s ) = EA k r ′ ( s ) k − k r ′ ( s ) k r ′ ( s ) is the tension of thestring. These are coupled ordinary and partial differen-tial equations. The motion of the reel mechanism andthe deployed portion of the string are described by (22)and (23), respectively. The translational and rotationaldynamics of the rigid body are determined by (24), (25).All of these equations are coupled.In (22), the fifth term, µǫ + ˙ s p , represents the effect ofthe velocity discontinuity. Note that this term vanishesif the deployed portion of the string is also inextensi-ble, i.e. k r ′ ( s ) k = 1 or ǫ + = 0. A similar expression isdeveloped in [10] from momentum balance, but theirexpression is erroneous. Special Cases
Suppose that the length of the string onthe reel mechanism is fixed, i.e. s p ( t ) ≡ s p ( t ) for any t > t . Then, the equations of motion (23)-(25) describethe dynamics of an elastic string pendulum model witha fixed unstretched length, which is studied in [11]. In this case, the total energy and the total angular mo-mentum about the gravity direction e are conserved: E = ( T r + V r ) + ( T s + V s ) + ( T b + V b ) ,π = e · (cid:20) Z Ls p µr ( s ) × ˙ r ( s ) ds + M r ( L ) × ( ˙ r ( L ) + R ˆ Ωρ c ) − M ˙ r ( L ) × Rρ c + JΩ (cid:21) . If we choose ρ c = 0, then the rotational dynamics ofthe rigid body (25) is decoupled from the other equa-tions. In this case, (22)-(24) describe the dynamics ofan elastic string attached to a point mass and a reelmechanism. Geometric numerical integration deals with numericalintegration methods that preserve geometric propertiesof a dynamic system, such as invariants, symmetries, re-versibility, or structure of the configuration manifold [14,19]. The geometric structure of a dynamic system de-termines its qualitative dynamical behavior, and there-fore, the geometric structure-preserving properties of ageometric numerical integrator play an important rolein the qualitatively accurate computation of long-termdynamics. The continuous-time Euler-Lagrange equa-tions developed in the previous section provide an ana-lytical model for a string pendulum. However, the pop-ular finite difference approximations or finite elementapproximations of those equations using a general pur-pose numerical integrator may not accurately preservethe geometric properties of the system [14].Variational integrators provide a systematic methodof developing geometric numerical integrators for La-grangian/Hamiltonian systems [20]. Discrete-time Euler-Lagrange equations, referred to as variational integra-tors, are constructed by discretizing Hamilton’s princi-ple, rather than discretizing the continuous-time Euler-Lagrange equations using finite difference approxima-tions. This is in contrast to the conventional viewpointthat a numerical integrator of a dynamic system is adiscrete approximation of its continuous-time equationsof motion. As it is derived from a discrete analogue ofHamilton’s principle, it preserves symplecticity and themomentum map, and it exhibits good total energy be-havior for an extremely long time period.On the other hand, Lie group methods conserve thestructure of a Lie group configuration manifold as it up-dates a group element using the group operation [21].As opposed to computational methods based on localcoordinates, projections, or constraints, this approachpreserves the group structure naturally without any singularities associated with local coordinates or theadditional computational overhead introduced by con-straints.These two methods have been unified to obtain a Liegroup variational integrator for Lagrangian/Hamiltoniansystems evolving on a Lie group [17]. This geometric in-tegrator preserves symplecticity and group structure ofthose systems concurrently. It has been shown that thisproperty is critical for accurate and efficient simulationsof rigid body dynamics [15]. This is particularly usefulfor dynamic simulation of a string pendulum that un-dergoes large displacements, deformation, and rotationsover an exponentially long time period.In this section, we develop a Lie group variationalintegrator for a string pendulum. We first construct adiscretized string pendulum model, and derive an ex-pression for a discrete Lagrangian, which is substitutedinto discrete-time Euler-Lagrange equations on a Liegroup.3.1 Discretized String Pendulum ModelLet h > t = t + kh is denoted by a subscript k . We discretizethe deployed portion of the string using N identical lineelements. Since the unstretched length of the deployedportion of the string is L − s p k , the unstretched length ofeach element is l k = L − s pk N . Let the subscript a denotethe variables related to the a -th element. The naturalcoordinate of the a -th element is defined by ζ k,a ( s ) = ( s − s p k ) − ( a − l k l k (26)for s ∈ [ s p k + ( a − l k , s p k + al k ]. This varies between 0and 1 for the a -th element. Let S , S be shape functionsgiven by S ( ζ ) = 1 − ζ , and S ( ζ ) = ζ . These shapefunctions are also referred to as tent functions . Define q k,a to be the relative location of a string element withrespect to the guide way entrance, i.e q k,a = r k,a − r p .Using this finite element model, the position vector r ( s, t ) of a material point in the a -th element is approx-imated as follows: r k ( s ) = S ( ζ k,a ) q k,a + S ( ζ k,a ) q k,a +1 + r p . (27)Therefore, a configuration of the presented discretizedstring pendulum at t = kh + t is described by g k =( s p k ; q k, , . . . , q k,N +1 ; R k ), and the corresponding con-figuration manifold is G = R × ( R ) N +1 × SO ( ). This isa Lie group where the group acts on itself by the diago-nal action [12]: the group action on s p k and q ,k · · · q N +1 ,k is addition, and the group action on R k is matrix mul-tiplication. We define a discrete-time kinematics equation usingthe group action as follows. Define f k = ( ∆s p k ; ∆q k, , . . . , ∆q k,N +1 ; F k ) ∈ G such that g k +1 = g k f k :( s p k +1 ; q k +1 , , . . . , q k +1 ,N +1 , R k +1 ) =( s p k + ∆s p k ; q k, + ∆q k, , . . . , q k,N +1 + ∆q k,N +1 ; R k F k ) . (28)Therefore, f k ∈ G represents the relative update be-tween two integration steps. This ensures that the struc-ture of the Lie group configuration manifold is numeri-cally preserved since g k is updated by f k using the rightLie group action of G on itself.3.2 Discrete LagrangianA discrete Lagrangian L d ( g k , f k ) : G × G → R is anapproximation of the Jacobi solution of the Hamilton–Jacobi equation, which is given by the integral of theLagrangian along the exact solution of the Euler-Lagrangeequations over a single time step: L d ( g k , f k ) ≈ Z h L (˜ g ( t ) , ˜ g − ( t ) ˙˜ g ( t )) dt, where ˜ g ( t ) : [0 , h ] → G satisfies Euler-Lagrange equa-tions with boundary conditions ˜ g (0) = g k , ˜ g ( h ) = g k f k .The resulting discrete-time Lagrangian system, referredto as a variational integrator, approximates the Euler-Lagrange equations to the same order of accuracy asthe discrete Lagrangian approximates the Jacobi solu-tion [17,20].We construct a discrete Lagrangian for the stringpendulum using the trapezoidal rule. We first find thecontributions of each component to the kinetic energy.From the given discretized string pendulum model and(2), the kinetic energy of the reel mechanism is approx-imated by T k,r = 12 h ( µs p k + κ d ) ∆s p k . (29)Using the chain rule, the partial derivative of r k ( s )given by (27) with respect to t is given by˙ r k ( s ) = 1 h (cid:26) S ( ζ k,a ) ∆q k,a + S ( ζ k,a ) ∆q k,a +1 + ( L − s )( L − s p k ) ( q k,a − q k,a +1 ) l k ∆s p k (cid:27) . Substituting this into (3), the contribution of the a -thstring element to the kinetic energy of the string is given by T k,a = 12 h M k ∆q k,a · ∆q k,a + 12 h M k ∆q k,a +1 · ∆q k,a +1 + 12 h M k,a ∆s p k + 1 h M k ∆q k,a · ∆q k,a +1 + 1 h M k,a ∆s p k · ∆q k,a +1 + 1 h M k,a ∆s p k · ∆q k,a , (30)where inertia matrices are defined in Appendix A.1.From the attitude kinetics equations (1), the angularvelocity is approximated byˆ Ω k ≈ h R Tk ( R k +1 − R k ) = 1 h ( F k − I ) . Define a non-standard inertia matrix J d = tr[ J ] I − J .Using the trace operation and the non-standard inertiamatrix, the last term of the kinetic energy of the rigidbody given by (2) can be written in terms of ˆ Ω as Ω · JΩ = tr[ ˆ ΩJ d ˆ Ω T ]. Then, the kinetic energy of the rigidbody is given by T k,b = 12 h M ∆q k,N +1 · ∆q k,N +1 + 1 h tr[( I − F k ) J d ]+ 1 h M ∆q k,N +1 · R k ( F k − I ) ρ c , (31)where we use properties of the trace operator: tr[ AB ] =tr[ BA ] = tr[ A T B T ] for any matrices A, B ∈ R × .From (29), (30), (31), the total kinetic energy of thediscretized string pendulum is given by T k = T k,r + N X a =1 T k,a + T k,b . (32)Similarly, the total potential energy for the given dis-cretized string pendulum can be written as V k = − µg (cid:8) ( s p k − d ) r d · e + d (cos(( s p − b ) /d ) − (cid:9) + N X a =1 − µgl k e · (2 r p + q k,a + q k,a +1 )+ 12 EAl k ( k q k,a +1 − q k,a k − l k ) − M ge · ( q k,N +1 + r p + R k ρ c ) . (33)This yields the discrete-Lagrangian of the discretizedstring pendulum L d k ( g k , f k ) = hT k ( g k , f k ) − h V k ( g k , f k ) − h V k +1 ( g k , f k ) . (34) 3.3 Lie Group Variational IntegratorFor a given discrete Lagrangian, the discrete action sumis given by G d = P k L d k . As the discrete Lagrangianapproximates the action integral over a single discretetime step, the action sum approximates the action in-tegral. According to the discrete Lagrange–d’Alembertprinciple, the sum of the variation of the action sumand the discrete virtual work done by external controlmoments and constraints is equal to zero. This yieldsdiscrete-time forced Euler-Lagrange equations referredto as variational integrators [20]. This procedure is fol-lowed for an arbitrary discrete Lagrangian defined on aLie group configuration manifold in [17] to obtain a Liegroup variational integrator: T ∗ e L f k − · D f k − L d k − − Ad ∗ f − k · ( T ∗ e L f k · D f k L d k )+ T ∗ e L g k · D g k L d k + u d k + Q d k = 0 , (35) g k +1 = g k f k , (36)where T ∗ L : G × T ∗ G → T ∗ G is the co-tangent lift ofthe left translation action, D f represents the derivativewith respect to f , and Ad ∗ : G × g ∗ → g ∗ is the co-Adjoint operator [12].The contribution of the external control momentand the Carnot energy loss term are denoted by u d k and Q d k . They are defined to approximate the addi-tional term in the variational principle (20) that arisesdue to a discontinuity: Z ( k +1) hkh ( Q + u/d ) δs p dt ≈ ( Q d k + u d k ) δs p k . From (21), these are chosen as Q d k = − h l k ( µ∆s p k /h + EA )( k q k, k − l k ) , (37) u d k = hu k /d. (38)We substitute the expressions for the discrete La-grangian (34), the Carnot energy loss term (37), and thecontrol moment (38) into (35) and (36) to obtain a Liegroup variational integrator for the discretized stringpendulum model. This involves deriving the derivativesof the discrete Lagrangian and their co-tangent lift. Thedetailed procedure and the resulting expressions for Liegroup variational integrators are summarized in the Ap-pendix. Computational Properties
The proposed Lie group vari-ational integrators have desirable computational prop-erties. The Lie group configuration manifold is oftenparameterized. But, the local parameterizations of thespecial orthogonal group, such as Euler angles or Ro-drigues parameters, have singularities. In a numerical simulation of large angle maneuvers of a string pen-dulum, these local parameters should be successivelyswitched from one type to another in order to avoidtheir singularities. They also lead to excessive complex-ity. Non-parametric representations such as quaternionsalso have associated difficulties: there is an ambiguityin representing an attitude since the group SU (2) ofquaternions double cover SO ( ). Furthermore, as theunit length of quaternions is not preserved in numeri-cal simulations, attitudes cannot be determined accu-rately. Sometimes, numerical solutions updated by anyone-step integration method are projected onto the Liegroup at each time step [22]. Such projections may de-stroy the desirable long-time behavior of one-step meth-ods, since the projection typically corrupts the numeri-cal results. Lie group variational integrators update thegroup elements by using a group operation. Therefore,the Lie group structure is naturally preserved at thelevel of machine precision, and they avoid any singular-ity and complexity associated with other approaches.Since Lie group variational integrators are constructedaccording to Hamilton’s principle, their numerical tra-jectories preserve a symplectic form and a momentummap associated with any symmetry. These ensure long-term structural stability and avoid artificial numericaldissipation. These properties are difficult to achieve inconventional approaches based on finite difference ap-proximation of continuous equations of motion.In summary, the proposed Lie group variational in-tegrators for a string pendulum will be particularly use-ful when studying nontrivial maneuvers that combinelarge elastic deformations and large rigid motions accu-rately over a long time period. We now numerically demonstrate the computationalproperties of the Lie group variational integrators de-veloped in the previous section. The properties of thereel mechanism are as follows: b = d = 0 . κ d = 1 kg.The material properties of the string are chosen to rep-resent a rubber string [5]: µ = 0 .
025 kg / m, EA = 40 N, L = 100 m. The rigid body is chosen as an elliptic cylin-der with a semimajor axis 0 . . . M = 0 . ρ c = [0 . , . , .
4] m, respectively.We consider three cases: (1) dynamics of a fixedlength string pendulum released from a horizontal con-figuration, (2) deployment dynamics due to gravity froma horizontal configuration, and (3) retrieval dynamicsusing a constant control moment. Initial conditions areas follows: s p (m) q ,a (m) u k (Nm)(1) 90 l e -(2) 99 l ( a − e l ( a −
1) (sin 15 ◦ e + cos 15 ◦ e ) 2.09 For all cases, we choose ˙ s p = 0 m / s, ˙ q ,a = 0 m / sfor a ∈ { , N } , ˙ q ,N +1 = 0 . e m / s, R = I , Ω =0 rad / sec. The deployed portion of the string is dis-cretized by N = 20 elements, and the time step is h = 0 . T = 10, T = 8,and T = 10 seconds for each case, respectively. Energy transfer
The following figures show the simu-lation results. The maneuver of the string pendulum isillustrated by snapshots, where the relative elastic po-tential distribution at each instant is denoted by colorshading (the corresponding animations are available athttp://my.fit.edu /˜taeyoung). As the point where thestring is attached to the rigid body is displaced from thecenter of mass of the rigid body, the rigid body dynam-ics are directly coupled to the elastic string dynamics.The illustrated maneuvers clearly show the nontrivialcoupling between the strain deformation, the rigid bodydynamics, and the reel mechanism.The energy exchange plots also show that there issignificant energy transfer between the kinetic energy,the gravitational potential energy, and the elastic po-tential energy. In Figure 2(b), there is an energy ex-change between the kinetic energy and the gravitationalpotential energy. But, when the string is mostly stretchedat t = 1 . t = 5 . t = 7 seconds. The elastic potential energytransfer along the string is observed in Figure 2(a). Forthe deployment case shown in Figure 3(b), the gravi-tational potential energy is generally transferred to thekinetic energy. As the length of the deployed portionof the string increases, the elastic potential energy in-creases and the rotational kinetic energy decreases. Forthe third retrieval case, both the gravitational potentialenergy and the total energy increase due to the constantcontrol moment. In addition, there is a smaller-scaleperiodic energy exchange between the elastic potentialand the gravitational potential energy with an approx-imate period of 1.2 seconds. The rigid body starts tum-bling at t = 3 seconds. Conservation Properties
The proposed Lie group vari-ational integrators exhibit excellent conservation prop-erties for these complicated maneuvers of the stringpendulum. For the fixed length string dynamics, thetotal energy and the total angular momentum aboutthe gravity direction should be preserved. The devia-tions of those quantities are shown in Figure 2(d), where the maximum deviation of the total energy is less than0 .
01% of the maximum kinetic energy, and the devia-tion of the angular momentum is less than 3 × − %of its initial value. For the second deployment case, thetotal energy dissipates only due to the velocity discon-tinuity. Figure 3(d) shows the difference between thecomputed total energy change and the energy dissipa-tion computed by the Carnot energy loss term (21).The difference is less than 0 . (cid:13)(cid:13) I − R T R (cid:13)(cid:13) , is less than 10 − . We have developed continuous-time equations of mo-tion and geometric numerical integrators, referred to asLie group variational integrators, for a 3D elastic stringpendulum attached to a rigid body and a reel mecha-nism. They are carefully derived while taking accountof the length change of the deployed portion of thestring, the Lie group configuration manifold of the rigidbody, and the velocity discontinuity at the guide wayentrance. The continuous-time equations of motion pro-vide an analytical model that is defined globally on theLie group configuration manifold. The Lie group varia-tional integrator preserves the geometric features of thesystem, thereby yielding a reliable numerical methodto compute the nonlinear coupling between the largestring deformation and the nontrivial rigid body dy-namics accurately over a long time period. In short,this paper provides high fidelity analytical and compu-tational models for a string pendulum.The numerical experiments suggest that accuratelymodeling the reeling mechanism is of critical impor-tance in order to capture the correct dynamics, due tothe disturbance that is introduced in the string at thepoint of contact with the reeling mechanism when thestring is deployed or retracted. One can observe thatthis disturbance propagates down the string at a veloc-ity that is determined by the elastic properties of thestring. Since the point of contact between the stringand the rigid body does not go through the center ofmass of the rigid body, the elastic disturbance excitesa rotational response in the rigid body. As such, ac-curately modeling the reel mechanism, elastic stringdynamics, rigid body motion, and their interactions,is critical for obtaining realistic predictions about howtowed underwater vehicles and tethered spacecraft be-have when performing aggressive maneuvers. The proposed string pendulum model and compu-tational approach can be extended in several ways. Forexample, different types of string models can be consid-ered, such as an inextensible string, nonlinear elasticity,and bending stiffness. The reel mechanism can be gen-eralized by assuming that the portion of the string onthe drum is also extensible. These results can be ex-tended to model tethered spacecraft in orbit, and theycan be used to study associated optimal control prob-lems by adopting the discrete mechanics and optimalcontrol approach [23]. t ∈ [0 , t E n e r gy ETT rot V gravity V elastic (b) Energy exchange ( E :solid, T :solid, T rot :dashed, V gravity :dash-dotted, V elastic :dotted) L e n g t h t Ω (c) Stretched length of the deployed portion of the string, and thesecond component of the angular velocity Ω
101 x 10 ∆ E
101 x 10 t ∆ π (d) Deviation of conserved quantities: total energy and the totalangular momentum about the gravity direction Fig. 2
Fixed length string pendulum3(a) Snapshots at each 0.4 second t ∈ [0 , t E n e r gy ET T rot V gravity V elastic (b) Energy exchange ( E :solid, T :solid, T rot :dashed, V gravity :dash-dotted, V elastic :dotted) L e n g t h t Ω (c) Length of the deployed portion of the string (stretched:solid,unstreched:dashed), and the second component of the angularvelocity Ω
202 x 10 ∆ E − Z Q ˙ s p d t t k I − R T R k (d) Deviation of conserved quantities: the difference between thecomputed total energy change and the energy dissipation due tothe velocity discontinuity, the orthogonality error of the rotationmatrix Fig. 3
Deployment due to gravity4 (a) Snapshots at each 0.5 second t ∈ [0 , t E n e r gy TET rot V gravity V elastic (b) Energy exchange ( E :solid, T :solid, T rot :dashed, V gravity :dash-dotted, V elastic :dotted) L e n g t h t Ω (c) Length of the deployed portion of the string (stretched:solid,unstreched:dashed), and the second component of the angularvelocity Ω t k I − R T R k (d) Orthogonality error of the rotation matrix Fig. 4
Retrieval using a constant control moment5
A Development of the Lie Group VariationalIntegrator for a String Pendulum
A.1 Inertia Matrices for the Discrete Lagrangian
The inertia matrices for the discrete Lagrangian are defined asfollows. M k = 13 µl k , M k = M k ,M k,a = 13 µl k (3 N + 3 N + 1 − Na − a + 3 a ) N ,M k = 16 µl k , M k,a = 16 µ (1 + 3 N − a ) N ( q k,a − q k,a +1 ) ,M k,a = 16 µ (2 + 3 N − a ) N ( q k,a − q k,a +1 ) . A.2 Derivatives of the discrete Lagrangian
The Lie group variational integrator given by (35) is expressed interms of the derivatives of the discrete Lagrangian and their co-tangent lift. Here, we describe how to compute the co-tangent liftand the co-Adjoint operator on the configuration manifold G = R × ( R ) N +1 × SO ( ) without introducing the formal definitionof those operators.The co-tangent lift of the left translation on a real space isthe identity map on that real space. Using the product struc-ture of the configuration manifold G = R × ( R ) N +1 × SO ( ),the derivative of the discrete Lagrangian with respect to f k =( ∆s p k ; ∆q k, , . . . , ∆q k,N +1 ; F k ) ∈ G is given by T ∗ e L f k · D f k L d k = ˆ D ∆s pk L d k ; D ∆q ,k L d k , · · · , D ∆q N +1 ,k L d k ; T ∗ I L F k · D F k L d k ˜ . (39)Deriving the derivatives of the discrete Lagrangian with re-spect to ∆s p k or ∆q k,a is straightforward. For example, from(30), (32), (34), the derivative of the discrete Lagrangian withrespect to ∆q k,a for any a ∈ { , . . . , N } is given by D ∆q k,a L d k = h D ∆q k,a T k,a − + h D ∆q k,a T k,a − − h D ∆q k,a V k +1 = 1 h M k ∆q k,a − + 2 h M k ∆q k,a + 1 h M k ∆q k,a +1 + 1 h ( M k,a + M k,a − ) ∆s p k − h D q k +1 ,a V k +1 , (40)where the derivative of the potential energy is given by D q k,a V k = − µgl k e + ∇ V ek,a − − ∇ V ek,a , (41) ∇ V ek,a = EAl k ‚‚ q k,a +1 − q k,a ‚‚ − l k ‚‚ q k,a +1 − q k,a ‚‚ ( q k,a +1 − q k,a ) . (42)Expressions for the other derivatives of the discrete Lagrangianwith respect to s p k , ∆s p k , q k,a are similarly developed and theyare summarized later.Now we find the derivative of the discrete Lagrangian withrespect to F k . From (31), (33), (34), we have D F k L d k · δF k = 1 h tr[ − δF k J d ] + Mh R Tk ∆q k,N +1 · δF k ρ c + h Mge · R k δF k ρ c . Similar to (16), the variation of the rotation matrix F k can bewritten as δF k = F k ˆ ζ for a vector ζ ∈ R . From the definition ofthe co-tangent lift of the left translation, we have( T ∗ I L F k · D F k L d k ) · ζ = 1 h tr h − F k ˆ ζ k J d i + Mh R Tk ∆q k,N +1 · F k ˆ ζ k ρ c + h Mge · R k F k ˆ ζ k ρ c . By repeatedly applying the following property of the trace op-erator, tr[ AB ] = tr[ BA ] = tr[ A T B T ] for any A, B ∈ R × ,the first term can be written as tr[ − F k ˆ ζ k J d ] = tr[ − ˆ ζ k J d F k ] =tr[ˆ ζ k F Tk J d ] = − tr[ˆ ζ k ( J d F k − F Tk J d )]. Using the property of thehat map, x T y = − tr[ˆ x ˆ y ] for any x, y ∈ R , this can be fur-ther written as (( J d F k − F Tk J d ) ∨ ) · ζ k . As y · ˆ xz = ˆ zy · x for any x, y, z ∈ R , the second term can be written as F Tk R Tk ∆q k,N +1 · ˆ ζ k ρ c = ˆ ρ c F Tk R Tk ∆q k,N +1 · ζ k . Using these, we obtain T ∗ I L F k · D F k L d k = 1 h ( J d F k − F Tk J d ) ∨ + Mh ˆ ρ c F Tk R Tk ∆q k,N +1 + h Mg ˆ ρ c F Tk R Tk e . (43)Expression for the derivatives of the discrete Lagrangian withrespect to R k is similarly developed.In summary, in addition to (40) and (43), derivatives of thediscrete Lagrangian are summarized as follows. D ∆s pk L d k = 1 h M k ∆s p k + 1 h N X a =2 ( M k,a + M k,a − ) · ∆q k,a + 1 h M k,N ∆q k,N +1 − h D s pk +1 V k +1 ,M k = µs p k + κ d + 13 µ ( L − s p k ) , D ∆q k,N +1 L d k = 1 h ( M k + M ) ∆q k,N +1 + 1 h M k ∆q k,N + 1 h M k,N ∆s p k + 1 h MR k ( F k − I ) ρ c − h D q k +1 ,N +1 V k +1 , D s pk L d k = µ h ∆s p k − µ Nh N X a =1 ( ∆q k,a · ∆q k,a + ∆q k,a +1 · ∆q k,a +1 + ∆q k,a · ∆q k,a +1 ) − h D s pk V k − h D s pk +1 V k +1 , D q k,a L d k = µ Nh (1 + 3 N − a ) ∆s p k ∆q k,a +1 − µ Nh ∆s p k ∆q k,a − µ Nh (5 + 3 N − a ) ∆s p k ∆q k,a − − h D q k,a V k − h D q k,a V k +1 , D q k,N +1 L d k = − µ Nh ∆s p k ∆q k,N +1 − µ Nh ∆s p k ∆q k,N − h D q k,N +1 V k − h D q k,N +1 V k +1 , T ∗ I L R k · D R k L d k = Mh (( F k − I ) ρ c ) ∧ R Tk ∆q k,N +1 + h Mg ˆ ρ c R Tk e + h MgF k ˆ ρ c F Tk R Tk e , D s pk V k = − µgr d · e + µgd sin(( s p k − b ) /d )+ 12 N N X a =1 µge · (2 r p + q k,a + q k,a +1 )6 + EAl k ( ‚‚ q k,a +1 − q k,a ‚‚ − l k ) , D q k,a V k = − µgl k e + ∇ V ek,a − − ∇ V ek,a , D q k,N +1 V k = − ( 12 µl k + M ) ge + ∇ V ek,N . The co-Adjoint map on a real space is the identity mapon that real space. The co-Adjoint map on SO ( ) is given byAd ∗ F − k p = F k p = ( F k ˆ pF Tk ) ∨ for any p ∈ ( R ) ∗ ≃ so (3) ∗ . Usingthe product structure of the configuration manifold, we haveAd ∗ f − k ( T ∗ e L f k · D f k L d k ) = ˆ D ∆s pk L d k ; D ∆q ,k L d k , · · · , D ∆q N +1 ,k L d k ; Ad ∗ F Tk ( T ∗ I L F k · D F k L d k ) ˜ , (44)whereAd ∗ F Tk ( T ∗ I L F k · D F k L d k ) = 1 h ( F k J d − J d F k ) ∨ + Mh F k ˆ ρ c F Tk R Tk ∆q k,N +1 + h MgF k ˆ ρ c F Tk R Tk e . (45) Discrete-time Euler-Lagrange Equations
Substituting thederivatives of the discrete Lagrangian given in (40), (43) and theappendix, the co-Adjoint map given by (44), and the contribu-tions of the external control moment and the Carnot energy lossterm (37), (38) into the Lie group variational integrator on an ar-bitrary Lie group given by (35), (36), we obtain the discrete-timeEuler-Lagrange equations of the string pendulum at (46)-(51).Equation (48) is satisfied for a ∈ { , . . . , N } , and (51) is satis-fied for a ∈ { , . . . , N } . For given g k = ( s p k ; q k, , . . . q k,N +1 ; R k ),we solve (46)-(50) for the relative update f k = ( ∆s p k ; ∆q k, , . . . , ∆q k,N +1 ; F k ). Then, the configuration at the next step g k +1 =( s p k +1 ; q k +1 , , . . . , q k +1 ,N +1 ; R k +1 ) can be obtained by (51). Thisyields a discrete-time Lagrangian flow map ( g k , f k ) → ( g k +1 , f k +1 ),and this is repeated. Special Cases
If we set ∆s p k ≡ k , then the discrete-time Euler-Lagrange equations and Hamilton’s equations providea geometric numerical integrator for a string pendulum modelwith a fixed unstretched string length, studied in [11]. If we chose ρ c = 0, then these equations describe the dynamics of an elasticstring attached to a point mass and a fixed pivot. Computational Approach
These Lie group variational inte-grators for a string pendulum are implicit: at each time step,we need to solve nonlinear implicit equations to find the relativeupdate f k ∈ G . Therefore, it is important to develop an effi-cient computational approach for these implicit equations. Thiscomputational method should preserve the group structure of f k , in particular, the orthogonal structure of the rotation matrix F k ∈ SO ( ). The key idea of the computational approach pro-posed in this paper is to express the rotation matrix F k in termsof a vector c k ∈ R using the Cayley transformation [14]: F k = ( I + ˆ c k )( I − ˆ c k ) − . (52)Since the rotation matrix F k represents the relative attitude up-date between two adjacent integration steps, it converges to theidentity matrix as the integration step h approaches zero. There-fore, this expression is valid for numerical simulations even thoughthe Cayley transformation is a local diffeomorphism between R and SO ( ).Our computational approach is as follows. The implicit equa-tions for F k given by (50) are rewritten in terms of a vector c k using (52), and a relative update expressed by a vector X k = [ ∆s p k ; ∆q k, , . . . ∆q k,N +1 ; c k ] ∈ R × ( R ) N +1 × R is solved byusing a Newton iteration. After the vector X k converges, the ro-tation matrix F k is obtained by (52).This computational approach is desirable, since the implicitequations are solved numerically using operations in a linear vec-tor space. The three-dimensional rotation matrix F k is computedby numerical iterations on R , and its orthogonal structure is au-tomatically preserved by (52). It has been shown that this compu-tational approach is so numerically efficient that the correspond-ing computational load is comparable to explicit integrators [15]. Acknowledgements
This research has been supported in partby National Science Foundation Grants DMS-0714223, DMS-0726263,DMS-0747659, ECS-0244977, CMS-0555797.
References
1. D. Chapman, “Towed cable behaviour during ship turningmaneuvers,”
Ocean Engineering , vol. 11, no. 4, pp. 327–361,1984.2. R. Driscoll, R. Lueck, and M. Nahon, “Development and val-idation of a lumped-mass dynamics model of a deep sea ROVsystem,”
Applied Ocean Research , vol. 22, no. 3, pp. 169–182,2000.3. T. Walton and H. Polacheck, “Calculation of transient mo-tion of submerged cables,”
Mathematics of Computation ,vol. 14, no. 69, pp. 27–46, 1960.4. C. Koh, Y. Zhang, and S. Quek, “Low-tension cable dynam-ics: Numerical and experimental studies,”
Journal of Engi-neering Mechanics , vol. 125, no. 3, pp. 347–354, 1999.5. A. Kuhn, W. Seiner, J. Zemann, D. Dinevski, and H. Troger,“A comparison of various mathematical formulations and nu-merical solution methods for the large amplitude oscillationsof a string pendulum,”
Applied Mathematics and Computa-tion , vol. 67, pp. 227–264, 1995.6. B. Buckham, F. Driscoll, and M. Nahon, “Development of afinite element cable model for use in low-tension dynamicssimulation,”
Journal of Applied Mechanics , vol. 71, pp. 476–485, 2004.7. W. Steiner, J. Zemann, A. Steindl, and H. Trooger, “Numer-ical study of large amplitude oscillations of a two-satellitecontinuous tether systems with a varying length,”
Acta As-tronautica , vol. 35, no. 9-11, pp. 607–621, 1995.8. M. Krupa, W. Poth, M. Schagerl, A. Steindl, W. Steiner,and H. Troger, “Modeling, dynamics and control of tetheredsatellite systems,”
Nonlinear Dynamics , vol. 43, pp. 73–96,2006.9. A. Banerjee and T. Kane, “Tether deployment dynamics,”
The Journal of the Astronautical Sciences , pp. 347–365,1982.10. K. Mankala and S. Agrawal, “Dynamic modeling and sim-ulation of satellite tethered systems,”
Transactions of theASME , vol. 127, pp. 144–156, 2005.11. T. Lee, M. Leok, and N. McClamroch, “Dynamics ofa 3D elastic string pendulum,” in
Proceedings of IEEEConference on Decision and Control , 2009, submitted.[Online]. Available: http://arxiv.org/abs/0903.033212. J. Marsden and T. Ratiu,
Introduction to Mechanics andSymmetry , 2nd ed., ser. Texts in Applied Mathematics.Springer-Verlag, 1999, vol. 17.13. E. Crellin, F. Janssens, D. Poelaert, W. Steiner, andH. Troger, “On balance and variational formulations of theequations of motion of a body deploying along a cable,”
Jour-nal of Applied Mechanics , vol. 64, pp. 369–374, 1997.71 h M k ∆s p k + 1 h N X a =2 ( M k,a + M k,a − ) · ∆q k,a + 1 h M k,N ∆q k,N +1 − h M k − ∆s p k − − h N X a =2 ( M k − ,a + M k − ,a − ) · ∆q k − ,a − h M k − ,N ∆q k − ,N +1 − µ h ∆s p k + µ Nh N X a =1 ( ∆q k,a · ∆q k,a + ∆q k,a +1 · ∆q k,a +1 + ∆q k,a · ∆q k,a +1 ) + h D s pk V k + hd u k = h l k ( µ∆s p k /h + EA )( ‚‚ q k, ‚‚ − l k ) , (46) q k, = 0 , (47)1 h M k ∆q k,a − + 2 h M k ∆q k,a + 1 h M k ∆q k,a +1 + 1 h ( M k,a + M k,a − ) ∆s p k − h M k − ∆q k − ,a − − h M k − ∆q k − ,a − h M k − ∆q k − ,a +1 − h ( M k − ,a + M k − ,a − ) ∆s p k − − µ Nh (1 + 3 N − a ) ∆s p k ∆q k,a +1 + µ Nh ∆s p k ∆q k,a + h D q k,a V k + µ Nh (5 + 3 N − a ) ∆s p k ∆q k,a − = 0 , (48)1 h ( M k + M ) ∆q k,N +1 + 1 h M k ∆q k,N + 1 h M k,N ∆s p k + 1 h MR k ( F k − I ) ρ c − h ( M k − + M ) ∆q k − ,N +1 − h M k − ∆q k − ,N − h M k − ,N ∆s p k − − h MR k − ( F k − − I ) ρ c + µ Nh ∆s p k ∆q k,N +1 + µ Nh ∆s p k ∆q k,N + h D q k,N +1 V k = 0 , (49)1 h ( F k J d − J d F Tk − J d F k − + F Tk − J d ) ∨ + Mh ˆ ρ c R Tk ( ∆q k,N +1 − ∆q k − ,N +1 ) − hMg ˆ ρ c R Tk e = 0 , (50) s p k +1 = s p k + ∆s p k , q k +1 ,a = q k,a + ∆q k,a , R k +1 = R k F k . (51)14. E. Hairer, C. Lubich, and G. Wanner, Geometric numericalintegration , ser. Springer Series in Computational Mechanics31. Springer, 2000.15. T. Lee, M. Leok, and N. H. McClamroch, “Lie group varia-tional integrators for the full body problem in orbital me-chanics,”
Celestial Mechanics and Dynamical Astronomy ,vol. 98, no. 2, pp. 121–144, June 2007.16. ——, “Lie group variational integrators for the full bodyproblem,”
Computer Methods in Applied Mechanics and En-gineering , vol. 196, pp. 2907–2924, May 2007.17. T. Lee, “Computational geometric mechanics and controlof rigid bodies,” Ph.D. dissertation, University of Michigan,2008.18. T. Lee, M. Leok, and N. McClamroch, “Optimal attitudecontrol of a rigid body using geometrically exact computa-tions on SO(3),”
Journal of Dynamical and Control Systems ,vol. 14, no. 4, pp. 465–487, 2008.19. B. Leimkuhler and S. Reich,
Simulating Hamiltonian Dy-namics , ser. Cambridge Monographs on Applied and Com-putational Mathematics. Cambridge University Press, 2004,vol. 14.20. J. Marsden and M. West, “Discrete mechanics and varia-tional integrators,” in
Acta Numerica . Cambridge Univer-sity Press, 2001, vol. 10, pp. 317–514.21. A. Iserles, H. Munthe-Kaas, S. Nørsett, and A. Zanna, “Lie-group methods,” in
Acta Numerica . Cambridge UniversityPress, 2000, vol. 9, pp. 215–365.22. L. Dieci, R. Russell, and E. van Vleck, “Unitary integra-tors and applications to continuous orthonormalization tech-niques,”
SIAM Journal on Numerical Analysis , vol. 31, no. 1,pp. 261–281, 1994.23. O. Junge, J. Marsden, and S. Ober-Bl¨obaum, “Discrete me-chanics and optimal control,” in