Frédéric Boyer

Guillaume De Nayer

Alban Leroyer

Michel Visonneau

Keywords: nonlinear Kirchhoff beam, geometrically exact, finite elements, cable, FSI, computational fluid dynamics

Geometrically Exact Kirchhoff Beam Theory: Application to Cable Dynamics

In this article, the finite element simulation of cables is investigated for future applica-tions to robotics and hydrodynamics. The solution is based on the geometrically exact approach of Cosserat beams in finite transformations, as initiated by Simo in the 1980s. However, the internal basic kinematics of the beam theory is not those of Reissner-Timoshenko but rather those of Kirchhoff. Based on these kinematics, the dynamic model adopted is a nonlinear extension of the so-called linear model of twisted and stretched Euler-Bernoulli beams. In agreement with the investigated applications, one or both of the ends of the cable are submitted to predefined motions. This model is also implemented into a computational fluid dynamics code, which solves the Reynolds-averaged Navier-Stokes equations. Regarding this last point, an implicit/iterative algorithm including a conservative load transfer for the variable hydrodynamic forces exerted all along the beam length has been used to reach a stable coupling. The relevance of the approach is tested through three advanced examples. The first is related to the prediction of cable motion in robotics. Then, the two last illustrations deal with fluidstructure interaction (FSI). A 2D classical benchmark in FSI is first investigated, and, at last, a computation illustrates the procedure in a 3D case.

Introduction

Cables and other flexible pipes are found in many applications of civil, aerospace, and mechanical engineering. They are also commonly used in robotics, to bring energy to the actuators of a robotic arm. In oceanic engineering, risers are used, for example, to convey fluids from the bottom of the sea up to the surface. In these applications, engineers need numerical models of nonlinear thin deformable structures in order to predict their motions and the internal stress they undergo for a given application. Hence, the need for numerical nonlinear beam models is still an intense subject of research. Among the most advanced finite element methods proposed for nonlinear beams, the geometrically exact approach initiated by Simo is probably one of the most efficient. This approach is based on the Cosserat assumption of rigid cross sections ͓1͔ and on Timoshenko's internal kinematics ͓2͔, as extended to finite transformations by Reissner ͓3͔. Following the pioneering work of Simo ͓4͔, Simo and Vu-Quoc ͓5,6͔, and Cardona and Gérardin ͓7͔, this approach was further improved through the work of Cardona et al. ͓8͔,Ibrahimbegovic et al. ͓9͔, Ibrahimbegovic ͓10͔, Ibrahimbegovic and Al Mikdad ͓11͔, Ibrahimbegovic and Mamouri ͓12͔, Chrisfield and Jelenic ͓13͔, Kapania and Li ͓14͔, and Zupan and Saje ͓15͔. As regards numerical performances, the good results of this approach are essentially due to the fact that it requires no approximation of the finite rotations of the cross sections except the unavoidable discretization of the structure domain into finite elements and of the time axis through the integration schemas ͓16͔.

However, when the "length/thickness" ratio of the beam becomes very high ͑i.e., when the structure is more of a rope or a cable than a beam͒, it becomes relevant to use the alternative internal kinematics of Kirchhoff and Clebsch ͓17͔, as is done in oceanic engineering for the simulation of very long undersea cables ͓18-21͔.

Based on the nonlinear Kirchhoff beam theory, much progress has been accomplished by the oceanic engineering community. However, most of the results related to undersea cables have been based to date on the finite difference methods ͑see Refs. ͓18-21͔͒ with a very few based on the usual ͑not geometrically exact͒ finite element method ͓22͔.

Hence, developing geometrically exact finite element methods for cables is still a relevant challenge in nonlinear structural dynamics. A first step toward this goal was made in Ref. ͓23͔, where a geometrically exact finite element method based on Kirchhoff kinematics is proposed. In this reference, two models are tested and compared with the Reissner-Simo model. Prolonging the linear beam theory, these models correspond to those of Euler-Bernoulli and Rayleigh ͓24͔. In comparison to the geometrically exact theory of Reissner-Timoshenko beams, these models are not defined on the same configuration space since, in the case of very slender beams, the orthonormal frames of the Reissner theory reduce to the field of the unit vectors tangent to the line of centroids, i.e., in the language of Cosserat media, a continuous set of directors ͓25,26͔. This difference requires replacing the Lie group SO͑3͒ of the configuration space of Reissner theory by the fibration S2 ϫ SO͑2͓͒ 27͔, where S 2 = SO͑3͒ / SO͑2͒ stands for the orientation of the cross section normals, which is parameterized by the position field of the centroid line while SO͑2͒ stands for the cross section roll around the centroid line. In the current article, the geometrical approach is applied to this definition of the configuration space in order to address the problem of the modeling and simulation of cables with a view to future application to the detection of collisions in robotics and to the installation of risers in oceanic engineering. In particular, this article proposes a model ͑and its finite element implementation͒ of a nonlinear Euler-Bernoulli beam whose ends are possibly controlled by imposed motion and submitted to external forces ͑gravity and fluid forces͒. The constraints are modeled through an augmented Lagrangian approach with penalty, consistent with the Kirchhoff parameterization. Furthermore, this beam solver is coupled to a computational fluid dynamics ͑CFD͒ code, enabling complex vortex induced vibration ͑VIV͒ configurations to be handled.

This article is structured as follows. We first start by establishing the basic kinematics of the cable and the parameterization of its configuration space ͑Sec. 2͒. Once the geometric model has been achieved, we derive, in Sec. 3, the weak forms of the dynamics of a cable constrained at its two ends to follow some imposed motions. Then, the numerical algorithm is explained in Sec. 4. It is based on a Newmark implicit integrator coupled to a Newton resolution algorithm of the resulting algebraic formulation. Finally, the algorithm must be completed with a finite element discretization ͑Sec. 5͒ and some updated formulae of the weak form and its linearization. This is the purpose of Sec. 6. In Sec. 7, the approach is illustrated by some simulation results obtained on three examples. The first is related to predicting the motion of a cable linked to a moving robot arm. The second one concerns the fluid-structure interaction ͑FSI͒ application field: The main steps of the flow/motion coupling with a CFD solver are first described, followed by a classical benchmark of a 2D filament clamped behind a square rigid body and then submitted to a VIV phenomenon. As the previous test-case, the third one is a FSI illustration: It was done to show the capabilities of the FSI code in 3D.

Basic Kinematics and Parameterization of the Cable Geometry

The ambient space is given by a fixed Cartesian frame ͑O , e 1 , e 2 , e 3 ͒. Throughout this article, the tensor product will be denoted by , the contracted product of tensors by a point, and for any V R3 , V ˆ= V ∧ R 3 R 3 will be the skew symmetric tensor such that ∀x R 3 : V ϫ x = V ˆ• x. We consider a cable with a constant circular cross section and suppose it to be straight when at rest. We denote by X, the material abscissa along the cable and by ͑t i ͑t , X͒͒ i=1,2,3 , the moving orthonormal frame rigidly attached to the X cross section ͑centered on its center of gravity͒ and such that t 1 supports the material axis of the cable ͑see Fig. 1͒.A tt = 0, the cable is at rest, and such that its material axis ͑centroidal line͒ is along ͑O , e 1 ͒, i.e., ͑t i ͑0,X͒ = e i ͒ i=1,2,3 . Adopting the Euler-Bernoulli model of thin beams, the cross sections remain orthogonal to the centroidal line ͑Kirchhoff assumption͒. Consequently, any transformation undergone by the cable can be parameterized by the composition of a translation changing the position field X ۋ r͑X͒ into X ۋ r ء ͑X͒ = r͑X͒ + d͑X͒͑ where d is a displacement field͒, and a rotation changing the orientation field X ۋ R͑X͒ with respect to the configuration at rest into ͑see Fig. 2͒ R

ء ͑X͒ = exp͑t ˆ1 ء ͒⌳ t 1 ͑t 1 ء ͒R͑X͒͑ 1͒
where we introduce the field of unit tangent to the reference line after the transformation X ۋ t 1 ء ͑X͒ = ץ X r ء / ʈץ X r ء ʈ and the following rotation compatible with the Kirchhoff assumption:1

⌳ t 1 ͑t 1 ء ͒ = ͑t 1 • t 1 ء ͒I 3 + 1 1+t 1 • t 1 ء ͑t 1 ϫ t 1 ء ͒ ͑t 1 ϫ t 1 ء ͒ + ͑t 1 ϫ t 1 ء ͒ ͑2͒
which is the unique rotation applying t 1 S 2 onto t 1 ء S 2 while leaving the unit vector t 1 ϫ t 1 ء / ʈt 1 ϫ t 1 ء ʈ invariant ͓28,29͔. Finally, exp͑t ˆ1 ء ͒, where "exp" denotes the exponential of matrices, 2 is the rotation about the vector t 1 ء of an angle , which stands for a roll angle about the centroidal line. Once such a set of transformations is defined, we can use it to parameterize the configuration space of the cable. For that, we just distinguish a configuration ⌽ = ͑r , R͒ of the cable and define any other configuration ⌽ ء = ͑r ء , R ء ͒ as the action onto this arbitrary configuration ͑said to be "of reference"͒ of any transformation of the form ͑d , exp͑t ˆ1 ء ͒⌳ t 1 ͑t 1 ء ͒͒ R 3 ϫ SO͑3͒. In the numerical treatment that follows, the reference configuration is not fixed once and for all but updated at each time step of the simulation. Such a choice corresponds to the so-called updated Lagrangian approach of nonlinear elasticity and has to be distinguished from the "total Lagrangian" approach, where the reference configuration is unique and usually coincides with the initial one. This choice is adopted in order to avoid the singularity of Eq. ͑2͒, which occurs when t 1 =-t ˆ1 ء . As a matter of fact, this singularity corresponds to an antipodal ͑related on S 2 ͒ rotation of the cross section, easily avoidable by updating the reference configuration sufficiently often. Finally, in this approach, the current rotation field R͑t ,.͒ : X ۋ R͑t , X͒ will be parameterized as ͑v / ʈvʈ͒+ sin͑ʈvʈ͒͑v / ʈvʈ͒ ˆ. R͑t,X͒ = exp͑͑t,X͒t ˆ1͑t,X͒͒⌳ t 1n ͑t 1 ͑t,X͒͒R n ͑X͒͑ 3͒

where t 1n ͑X͒ = ץ X r n / ʈץ X r n ʈ and r n stands for the position field of the centroidal line at the last time step t n , where the equilibrium has been numerically reached by the algorithm. In the same manner, R n ͑X͒ is the orientation tensor of the X cross section at time t n measured with respect to the initial straight configuration.

Governing Balance Equations: Weak Form

In what follows, we denote by a prime symbol and a dot symbol, the operators ץ X and ץ t , respectively, and suppose the material abscissa of the beam normalized ͑X ͓0,1͔͒.

Strain and Stress Fields.

In the framework of finite displacements and small strains and under the Kirchhoff assumption, the strain state of the cable is defined by the three following fields ͑a vector with the subscript Ќ ͑respectively, with a ʈ͒ is normal ͑respectively, parallel͒ to the beam axis.

1. The vector field of curvature k Ќ :

X ۋ t 1 Ј=͑t 1 Ј•t 2 ͒t 2 + ͑t 1 Ј•t 3 ͒t 3 . 2. The scalar field of torsion K : X ۋ K = t 2 Ј•t 3 =-t 3 Ј•t 2 .
3. The scalar field of stretching : X ۋ = ʈrЈʈ -1, where we note that we have the following subsidiary relation:

RЈR T = ͑k Ќ + Kt 1 ͒ ˆ= ͑k Ќ + k ʈ ͒
ând that these strain measurements verify the material objectivity ͓30͔ required by the modeling of finite displacements.

On the dual side, we have the following internal force and torque densities per unit of beam length ͑the convention on repeated subscript ␣ =2,3 is adopted͒:

m = Mt 1 + M ␣ t ␣ = m ʈ + m Ќ , n = Nt 1 + N ␣ t ␣ = n ʈ + n Ќ
Finally, occurring in the framework of small strains and material elasticity, the internal stress and strains are related by the reduced Hook's law ͓4͔

m ʈ = GI p Kt 1 , m Ќ = EI a t 1 Ј, n ʈ = EAt 1
with EI a , GI p , and EA as the bending stiffness, twisting stiffness, and stretching stiffness per unit of beam length, respectively.

Statement of the Variational Formulation.

For the moment, we only consider the case where the external conservative forces are due to gravity, and we will consider more general nonconservative forces at the end of the section. Under the Euler-Bernoulli assumption, the rotational kinetic energy of cross sections is neglected and the Lagrangian of the free cable ͑without constraint͒ is simply

L f = 1 2 ͵ 0 1 Ar ˙2dX - 1 2 ͵ 0 1 EI a t 1 Ј 2 dX - 1 2 ͵ 0 1 GI p K 2 dX - 1 2 ͵ 0 1 EA 2 dX - 1 2 ͵ 0 1 Ag • rdX ͑4͒
where, from left to right, we find the kinetic energy of the cable, its bending energy, its torsion and stretching energies and, finally, the potential energy due to gravity ͑g and A are the gravity acceleration field and the mass density per unit of beam length, respectively͒.

In what follows, we consider the case of a cable constrained at one or both of its ends to follow some imposed motions. Thus, for a given cable configuration ⌽ = ͑r , R͒ with R in the form Eq. ͑3͒, the geometric constraints at the two ends take the form

⌼͑⌽,t͒ = ͫ ⌼ 0 ͑⌽,t͒ ⌼ 1 ͑⌽,t͒ ͬ=0, with ⌼ X ͑⌽,t͒ = ΄ r͑X͒ -r d ͑X,t͒ rЈ͑X͒ -r d Ј͑X,t͒ ͑X͒ -d ͑X,t͒ ΅ ͑5͒
The case with only one constrained end is a particular subcase, where ⌼ 0 ͑respectively, ⌼ 1 ͒ has to be canceled. In ⌼ X , the first ͑from top to bottom͒ set of constraints imposes a desired position of the X cross section while the next two ͑with defined by 3 for a field of arbitrary reference orientations R n ͒ completely specify its orientation. With these constraints, we define the augmented Lagrangian of the cable

L + = L f + • ⌼ ͑6͒
where we introduced the vector of Lagrange multipliers

= ͑ T ͑0͒ , T ͑l͒͒ T with ͑X͒ = ͑ r T ͑X͒ , r Ј T ͑X͒ , ͑X͒͒ T .
The Hamilton principle of constrained systems postulates that the motion of the cable makes the augmented action stationary

͓8͔, i.e., ␦ ͵ t 1 t 2 L + dt =0, ∀ ␦⌽, ∀ ␦ ͑7͒
where ␦⌽ = ͑␦r , ␦R͒ and ␦R are compatible with Kirchhoff kinematics; it has to verify the infinitesimal version 3 of Eq. ͑1͒

␦R = ͯ d d⑀ ͯ ⑀=0 ͑exp͑⑀␦t ˆ1͒⌳ t 1 ͑t 1,⑀ ͒R͒ = ͑␦t 1 + t 1 ϫ ␦t 1 ͒ ˆR ͑8͒
where we introduce the scalar field of variation ␦ = ␦t 2 • t 3 = -␦t 3 • t 2 , and where the fields subscripted by ⑀ are deduced from those nonsubscripted by substituting the field r by its perturbation r ⑀ = r + ⑀␦r into their definition. As the variation is achieved with fixed time and ends, Eq. ͑7͒ details as

͵ t 1 t 2 ␦L f + ␦⌼ • + ␦ • ⌼dt =0, ∀ ␦⌽, ∀ ␦ ͑9͒
Moreover, computing the variation of Eq. ͑3͒ with ⑀ = + ⑀␦, and identifying the result with Eq. ͑8͒, gives

␦R = ͑␦t 1 + t 1 ϫ ␦t 1 ͒ ˆR = ͯ d d⑀ ͯ ⑀=0 ͑exp͑ ⑀ t ˆ1,⑀ ͒⌳ t 1n ͑t 1,⑀ ͒R n ͒
After computation, we find the following scalar relation, which

links the iterative roll angle ␦ with the incremental one ␦:

␦ = ␦ + u n ͑X͒ • ␦rЈ ͑10͒
where we introduced the vector field

u n = t 1n ϫ t 1 ʈrЈʈ͑1+t 1n • t 1 ͒ ͑11͒
Then, inserting Eq. ͑10͒ into ␦⌼ • of Eq. ͑9͒ allows one to write

␦L + = ␦L f + ΄ ␦r͑0͒ ␦rЈ͑0͒ ␦͑0͒ ΅ • f c ͑0͒ + ΄ ␦r͑1͒ ␦rЈ͑1͒ ␦͑1͒ ΅ • f c ͑1͒ + ␦ • ⌼ =0

͑12͒

where we introduce the following generalized forces induced by the constraints:

f c = ΄ r r Ј + u n ΅ ͑13͒
Finally, the variation of K is defined by

␦K = ͯ d d⑀ ͯ ⑀=0 ͑t 2,⑀ Ј • t 3,⑀ ͒ =-ͯ d d⑀ ͯ ⑀=0 ͑t 3,⑀ Ј • t 2,⑀ ͒
and after computation, it can be detailed as

␦K = ␦Ј + ͩ t 1 ϫ t 1 Ј ʈrЈʈ ͪ•␦rЈ ͑14͒
Now, let us apply the extended Hamilton principle

␦ ͵ t 1 t 2 L + dt = ͵ t 1 t 2 ␦W ext dt, ∀ ␦⌽, ∀ ␦ ͑15͒
where ␦W ext is the virtual work produced by an eventual field of nonconservative external fluid forces ͑such as those applied on an underwater cable͒. Such a field of forces will be approximated through a wrench density of follower loads of the form

͓0,1͔ → R 3 ϫ SO͑3͒ ء X ۋ ͑n ¯,m ¯͒ = ͑N ¯it i ,M ¯it i ͒͑ 16͒
where N ¯i and M ¯i for i =1,2,3 are considered in the following as explicit time-dependent functions. Finally, by inserting Eq. ͑14͒ into Eq. ͑15͒ and achieving standard integrations by parts, the set of partial differential equations ͑PDEs͒ in the Cauchy form ͑i.e., the strong form of a nonlinear Euler-Bernoulli clamped at both ends͒ could be written. Since numerical applications are involved in this paper, only the weak form is under consideration in what follows.

Weak Form of Virtual

Works. Now, let us consider the following Lagrangian:

L + = L f + • ⌼ + p 2 ⌼ • ⌼ ͑17͒
where we add to Eq. ͑6͒ a penalty term ͑p is an arbitrary factor͒ in order to improve the convergence of the subsequent algorithm ͓8͔.

Then, repeating the same calculus as before ͑starting again from Eq. ͑15͒͒ but without achieving the spatial integration by parts,4 we obtain the weak form of virtual works that will be used later to solve the finite element formulation

͵ 0 1 ΄ ␦r ␦rЈ ␦ ΅ • ͑f dyn -f g -f ext ͒dX + ͵ 0 1 ΄ ␦rЈ ␦rЉ ␦Ј ΅ • f stat dX + ␦⌼ • ͑ + p⌼͒ + ␦ • ⌼ =0 ͑18͒ where ͑f dyn -f g -f ext ͒ T = ͩA͑r ¨-g͒ T -n ¯T,- ͑m ¯ϫ t 1 ͒ T ʈrЈʈ ,-t 1 • m ¯ͪ͑19͒
and

f stat = ΄ EAͩ1- 1 ʈrЈʈ ͪrЈ+EI a ͩ2 ͑rЈ • rЉ͒ 2 ʈrЈʈ 6 - ͑rЈ • rЉ͒rЉ + ʈrЉʈ 2 rЈ ʈrЈʈ 4 ͪ EI a ͩ rЉ rЈ 2 - ͑rЈ • rЉ͒ ʈrЈʈ 4 rЈͪ 0 ΅ + ΄ GI p Kͩ rЈ ϫ rЉ ʈrЈʈ 3 ͪ 0 3ϫ1 GI p K ΅ ͑20͒
Moreover, inserting Eq. ͑10͒ into ␦⌼ • ͑ + p⌼͒ allows one to write the weak form of virtual works ͑denoted G in what follows͒ as

G = ͵ 0 1 ΄ ␦r ␦rЈ ␦ ΅ • ͑f dyn -f g -f ext ͒ + ΄ ␦rЈ ␦rЉ ␦Ј ΅ • f stat dX + ΄ ␦r͑0͒ ␦rЈ͑0͒ ␦͑0͒ ΅ • f c ͑t,0͒ + ΄ ␦r͑1͒ ␦rЈ͑1͒ ␦͑1͒ ΅ • f c ͑t,1͒ + ␦ • ⌼ =0

͑21͒

where f c ͑t ,0͒ and f c ͑t ,1͒ are given by

f c ͑t,X͒ = ΄ r + p͑r -r d ͑t,X͒͒ r Ј + ͑ + p͑ -d ͑t,X͒͒͒u n ͑X͒ + p͑rЈ -r d Ј͑t,X͒͒ + p͑ -d ͑t,X͒͒ ΅ ͑22͒

Principle of the Numerical Algorithm

At this point, the dynamic problem turns into a variational problem, which consists in finding at each time t, the fields R, r, r ˙, r ¨, and , satisfying Eq. ͑21͒ for any ␦⌽ and ␦ imposed on the cable. The numerical algorithm proposed to solve this problem is based on a Newmark time integrator coupled to a Newton algorithm ͓31͔. Moreover, due to the Euler-Bernoulli assumption, only the position field r is concerned by the time integration, the torsion state being governed by the static equilibrium. This has two important consequences on the following: first, the "roll dynamics" can be interpreted as infinitely fast ͑in the sense of the "singular perturbation theory"͒ and, second, the roll dynamics will present some singularities when the boundary conditions at both ends of the cable are free. Fortunately, this case is not relevant for standard applications ͑see the examples at the end of this article͒. Whatever the case, if the cable is free at both ends, it is necessary to invoke a more refined model, such as that of Rayleigh ͓23͔ or Reissner ͓6͔͑ if the length/thickness ratio is not too high͒, which both take into account the rotational kinetic energy of cross section, and so removing the eventual roll singularities. However, these two models require the development of a time integrator on the curved manifolds SO͑3͓͒ 32͔ or S 2 ϫ SO͑2͓͒ 23͔ or, alternatively, the introduction of a vector parameterization of finite rotation ͓7͔ so compromising the simplicity of the geometrically exact approach. In our case, since the field r evolves in a vector space, the time integration has been readily achieved through the Newmark algorithm, which is quite standard for structural dynamics. In particular, no special work has been performed here to enhance energy and momentum conservation properties of the time integration.

In order to apply the Newton algorithm between two time steps, we have to compute the Jacobian of the resulting algebraic nonlinear problem. This is achieved, thanks to the linearization of the weak form Eq. ͑21͒ under the constraints imposed by the Newmark integration scheme. First of all, let us note that these constraints force the velocity and acceleration increment constructed on two steps of the Newton algorithm to verify

⌬r ˙͑t,.͒ = ␥ ␤⌬t ⌬r͑t,.͒, ⌬r ¨͑t,.͒ = 1 ␤⌬t 2 ⌬r͑t,.͒͑ 23͒
where, for any field f, we use the notation ⌬f͑t ,.͒ = f ͑k͒ ͑t ,.͒ f ͑k-1͒ ͑t ,.͒, with f ͑k͒ denoting the field computed at the kth Newton iteration. Then, the linearized weak form required by the Newton algorithm is computed by applying the operator "second variation" denoted by ⌬ to Eq. ͑21͒͑ see Ref. ͓33͔͒. Finally, we obtain

⌬G = ͵ 0 1 ΄ ␦r ␦rЈ ␦ ΅ • ͑⌬f dyn -⌬f g -⌬f ext ͒ + ΄ ␦rЈ ␦rЉ ␦Ј ΅ • ⌬f stat dX + ΄ ␦r͑0͒ ␦rЈ͑0͒ ␦͑0͒ ΅ • ⌬f c ͑t,0͒ + ΄ ␦r͑1͒ ␦rЈ͑1͒ ␦͑1͒ ΅ • ⌬f c ͑t,1͒ + ␦ • ⌬⌼ =0

͑24͒

where we make the second variations appear ͑produced by ⌬⌽ = ͑⌬r , ⌬R͒, with ⌬R = ͑⌬t 1 + t 1 ϫ⌬t 1 ͒ ˆR, i.e., by an infinitesimal transformation compatible with Kirchhoff kinematics͒

⌬f dyn = ΄ A⌬r 03ϫ1 0 ΅ , ⌬f g = 0 7ϫ1 , ⌬f ext = K ext ΄ ⌬r ⌬rЈ ⌬ ΅ , ⌬f stat = K stat ΄ ⌬rЈ ⌬rЉ ⌬Ј ΅ ⌬f c ͑t,X͒ = K c ͑t,X͒ ΄ ⌬r͑X͒ ⌬rЈ͑X͒ ⌬͑X͒ ΅ + ␦⌼ X • ⌬͑X͒͑ 25͒
Furthermore, we introduce in Eq. ͑25͒ the tangent stiffness operators of external, internal, and constraint forces denoted by K ext , K stat , and K c , respectively, and detailed in Appendices A and B. Then, identifying the second variations ⌬r ˙and ⌬r ¨with the Newton increment Eq. ͑23͒ allows one to rewrite the tangent dynamics Eq. ͑24͒ as

⌬G = ͵ 0 1 ΄ ␦r ␦rЈ ␦ ΅ • ͑K dyn -K ext ͒ ΄ ⌬r͑X͒ ⌬rЈ͑X͒ ⌬͑X͒ ΅ + ΄ ␦rЈ ␦rЉ ␦Ј ΅ • K stat ΄ ⌬rЈ͑X͒ ⌬rЉ͑X͒ ⌬Ј͑X͒ ΅ dX + ΄ ␦r͑0͒ ␦rЈ͑0͒ ␦͑0͒ ΅ • K c ͑t,0͒ ϫ ΄ ⌬r͑0͒ ⌬rЈ͑0͒ ⌬͑0͒ ΅ + ΄ ␦r͑1͒ ␦rЈ͑1͒ ␦͑1͒ ΅ • K c ͑t,1͒ ΄ ⌬r͑1͒ ⌬rЈ͑1͒ ⌬͑1͒ ΅ + ␦⌼ • ⌬ + ⌬⌼ • ␦ =0 ͑26͒
where K dyn is the dynamic tangent stiffness operator given in Appendix A, whereas the last two constraint terms will be detailed in Sec. 5.

Finite Element Method

We now have all the ingredients required by the finite element resolution. For that, let us discretize the weak form Eq. ͑21͒ and its linearization Eq. ͑26͒ as follows. The beam of initial length l = 1 is first of all subdivided into P elements subscripted e. Each element has a length l e and a first point X e ͑going along the beam from X =0 to l͒.

͑28͒

where ͑N I , M I ͒ and L I are the Hermit and Lagrange polynomials, respectively. Note that the summation convention on repeated subscripts ͑I =1,2͒ is adopted. Finally, considering the partition Eq. ͑27͒ and inserting Eq. ͑28͒ in Eqs. ͑21͒ and ͑26͒ give the discrete equilibrium and its discrete linearization ͑a tilde symbol denotes a discrete vector coming from a continuous field͒

G ˜= ͑␦q T ,␦ T ͒ ͫ f ˜dyn + f ˜stat -f ˜ext + f ˜c ⌼ ͬ = ͑␦q T ,␦ T ͒r ˜+ =0 ͑29͒ and ⌬G ˜= ͑␦q T ,␦ T ͒ ͫ K ˜dyn + K ˜stat -K ˜ext + K ˜c J ˜T J ˜0 ͬ ͫ ⌬q ⌬ ͬ = ͑␦q T ,␦ T ͒K ˜+ͫ ⌬q ⌬ ͬ=0 ͑30͒
where ␦q and ⌬q are the 7͑P +1͒ ϫ 1 vectors of nodal variations appearing in Eq. ͑28͒ while f ˜dyn , f ˜stat , and f ˜ext and K ˜dyn , K ˜stat , and K ˜ext are 7͑P +1͒ ϫ 1 force vectors and 7͑P +1͒ ϫ 7͑P +1͒ band matrices, respectively, obtained by assembling the corresponding elementary forces and stiffness matrices. On the other hand, the constraint forces and stiffness do not require any assembling procedure since we simply have f ˜c = ͑f c T ͑t ,0͒ , 0 1ϫ7͑P-1͒ , f c T ͑t ,1͒͒ T and

K ˜c = diag e=1,. . .,P+1 ͑K c,e ͒, with K c,e = 0 7ϫ7 , except K c,1 = K c ͑t ,0͒
and K c,P+1 = K c ͑t ,1͒. Finally, due to the discretization of ␦⌼⌬

+ ⌬⌼␦, we also introduce the 14ϫ 7͑P +1͒ matrix J ˜given by

J ˜= ͫ J͑0͒ 0 7ϫ7P 0 7ϫ7P J͑1͒ ͬ, with J͑X͒ = ΄ I 3 0 3ϫ3 0 3ϫ1 0 3ϫ3 I 3 0 3ϫ1 0 1ϫ3 u n T ͑X͒ 1 ΅ ͑31͒
as well as the augmented residue vector and stiffness matrix denoted by r ˜+ and K ˜+, respectively. Finally, once discretized, the linear expansion of the equilibrium at t + ⌬t gives the linear discrete system

K ˜+͑t+⌬t͒ ͑k͒ ⌬q ͑t+⌬t͒ ͑k+1͒ =-r ˜+͑t+⌬t͒ ͑k͒ ͑32͒
that we solve at each Newton step k with respect to the vector of nodal increments ⌬q t+⌬t ͑k+1͒ .

Updating Formulae

In order to solve the balance equilibrium, it is necessary to update the forces and the tangent stiffness operators at each Newton iteration k. This is achieved by updating the position, tangent, and curvature fields r, rЈ, and rЉ, the velocity and acceleration fields r ˙and r ¨, the torsion rate field K, the roll angles at the two boundaries ͑0͒ and ͑1͒, as well as the associated vectors of Lagrange multipliers ͑0͒ and ͑1͒. We give below the updating formulae of these fields, at any point X ͓X e , X e + l e ͔. To be used, these formulae first require Eq. ͑32͒ to be solved with respect to the vector of nodal increments ⌬q ͑t+⌬t͒ ͑k+1͒ and second, Eq. ͑28͒ to be invoked to get the current incremental fields that we simply note ⌬r, ⌬rЈ, and ⌬rЉ. These are the following. The following field is not explicitly required by the formulation; nevertheless, it is of practical interest since it intervenes as an intermediate variable in many calculations. 3. Updating of the unit tangent field t 1,͑t+⌬t͒ ͑k+1͒ ͑X͒ = r ͑t+⌬t͒ Ј ͑k+1͒ ͑X͒ ʈr ͑t+⌬t͒ Ј ͑k+1͒ ͑X͒ʈ 4. Updating of the orientation field. Replacing the two arbitrary configurations of Eq. ͑1͒, by two successive configurations of the Newton algorithm, the updating of the rotation field is done as follows:

R ͑t+⌬t͒ ͑k+1͒ = exp͑⌬͑t 1 + ⌬t 1 ͒ ˆ͒⌳ t 1 ͑t 1 + ⌬t 1 ͒R ͑t+⌬t͒ ͑k͒ ͑33͒ with t 1 = t 1,͑t+⌬t͒ ͑k͒ ͑X͒, ⌬t 1 = t 1,͑t+⌬t͒ ͑k+1͒ ͑X͒ -t 1,͑t+⌬t͒ ͑k͒ ͑X͒, ⌬ = ⌬ ˜͑t+⌬t͒ ͑k͒ ͑X͒
5. Updating of the torsion rate. The updating formula of K results from the answer to the following question: "What is the finite increment ⌬K generated by a finite variation of the configuration defined by ⌬r and ⌬?" To get this result, we can, for instance, apply the Frenet formulae K = t 2 Ј•t 3 = -t 3 Ј•t 2 to the updated configuration of the beam, and obtain after computations

K ͑t+⌬t͒ ͑k+1͒ = K ͑t+⌬t͒ ͑k͒ - 1 1+t 1,͑t+⌬t͒ ͑k͒ • t 1,͑t+⌬t͒ ͑k+1͒ ͓͑͑t 1,͑t+⌬t͒ ͑k+1͒ ϫ t 1,͑t+⌬t͒ Ј ͑k+1͒ ͒ • t 1,͑t+⌬t͒ ͑k͒ ͒ -͑t 1,͑t+⌬t͒ ͑k͒ ϫ t 1,͑t+⌬t͒ Ј ͑k͒ ͒ • t 1,͑t+⌬t͒ ͑k+1͒ ͔ + ⌬Ј 6.
Updating of the roll angle. After having updated the rotation field as in Eq. ͑33͒, we can extract from it the new incremental roll angle at any point of the beam. This is achieved by inverting the identity Eq. ͑3͒, where the current configuration here is the ͑k +1͒th configuration calculated by the Newton algorithm

͑t+⌬t͒ ͑k+1͒ = exp t ˆ͑t+⌬t͒ ͑k+1͒ -1 ͑⌳ t 1n T ͑t 1 ͑t+⌬t͒ ͑k+1͒ ͒R n T R ͑t+⌬t͒ ͑k+1͒ ͒
where we introduce the inverse map exp V -1 such that exp V -1 ͑exp͑␣V ˆ͒͒= ␣ ͓0,2͔. 7. Updating of the Lagrange multipliers. Because they belong to a linear space, the Lagrange multipliers are updated by the simple additive formula

͑t+⌬t͒ ͑k+1͒ = ͑t+⌬t͒ ͑k͒ + ⌬ ͑t+⌬t͒ ͑k+1͒ , with ͑t+⌬t͒ ͑0͒ =0

Numerical Illustrations

In what follows, the previous model is illustrated on three cases related to robotics and fluid-structure interaction fields. Due to weak instabilities introduced by the algebraic constraints in the first case and in order to reach a stable FSI coupling while keeping standard time step in the second, a modified version of the Newmark time scheme has to be used for the structural solver such as the Hilber-Hughes-Taylor ͑HHT͓͒34͔, the Bossak ͓35͔ method, or the ␣-generalized technique ͓36͔. These algorithms introduce a numerical dissipation of the higher frequencies without degrading the accuracy. In the following results, the Bossak scheme was used with an ␣ coefficient fixed to 0.1.

Robotics Application: Cable and Six-Axis Puma

Robot. In the robotics field, cables or flexible tubes are found everywhere since they bring electric or hydraulic power to the actuators. The main problems are how to guide them, how to avoid collisions, and how to minimize their wear. The present illustrations will show that the beam solver can predict the movement and positions of a complex robotics problem.

We consider a six degrees of freedom ͑DOFs͒ manipulator of the Puma type ͑cf. Fig. 3͒. The time evolution of the joints is defined in Table 1. The cable is clamped to the robot at the origin of the tool frame O 7 , and to the earth at the point C of the environment. Its characteristics are presented in Table 2. Figure 3 shows the robot and the cable in the initial configuration of the motion. The initial configuration of the cable is deduced from a vertical straight reference configuration by applying the constraint Eq. ͑5͒ to it with an imposed motion compatible with its manipulation by a person. This preliminary mounting phase first considers the cable embedded in C and lying under gravity. Then, the free end of the cable is brought along the circular ␣-parameterized path defined in Table 3 until it reaches the point O 7 while the robot is in the initial configuration of Fig. 3.

Once the mounting phase has been achieved, the tool begins to move at t = 0, according to the imposed motion of the robot ͑Table 1͒. The time-varying constraints ͑see Eq. ͑5͒͒ imposed on the moving end of the cable are deduced from the direct geometric model of the Puma type manipulator with Denavit-Hartenberg parameters ͓37͔ r 1 = d 4 = 1 m and r 5 = r 7 = 0.3 m. The spatial discretization uses ten finite elements. The mounting phase is first performed using a static version of the algorithm. Then, the simulation is achieved over a duration of4sw i t hat i m estep of ⌬t = 0.001 s. Figure 4, which shows snapshots at every 25 ms, highlights the dynamic effects such as the pendular oscillations due to gravity and the high velocity of the robot.

Table 4 shows the typical evolution of the residual vector norm during the iterations of the beam solver for one time step ͑the convergence for the others is quite similar͒.

Today, this solver is in progress toward virtual reality and haptic applications. In this case, it becomes necessary to address the problem of collision modeling, which is not considered here.

2D FSI Illustration:

Hübner Test-Case. Another application, where this model is relevant, concerns risers used in oceanic engineering. In this case, FSI is involved. Thus, the previous cable model has to be coupled with a fluid solver. As a first step toward this goal, we consider here the benchmark due to Hübner. This FSI application was originally carried out by Wall and Ramm ͓38͔ in order to demonstrate the ability of their FSI formulation to deal with complex flow-flexible structure interactions exhibiting large deformations. In 2004, Huebner et al. modified the inflow velocity and used this new problem to test his monolithic code

͓39͔.

The application consists of a thin elastic nonlinear beam attached to a fixed square rigid body, plunged in an incompressible fluid flow ͑cf. Fig. 5͒. The square rigid body produces a very unsteady phenomenon: Vortices separate alternatively from its corners. These vortices create an oscillating force on the filament, which deforms itself.

The fluid and structure characteristics are given in Table 5. The Reynolds number is fixed to 204; the computation is consequently done with a laminar flow. The boundary conditions are set to the following:

1. imposed velocity at the inflow, upper, and lower boundaries 2. imposed pressure at the outflow 3. adherence boundary on the filament Two meshes are generated using the commercial software GRIDGEN ™ . One is totally unstructured and coarse ͑Ӎ4000 cells and Ӎ2000 nodes͒. The second is finer ͑Ӎ33,000 cells and Ӎ27,000 nodes͒ and it is structured near the body and then fully unstructured.

Flows are computed using an in-house solver ISIS-CFD, developed by the Equipe Modélisation Numérique ͑EMN͒ of the Fluid Mechanics Laboratory of the Ecole Centrale Nantes. It solves the incompressible unsteady Reynolds-averaged Navier-Stokes ͑RANS͒ equations under isothermal conditions on generalized unstructured meshes, enabling complex geometries to be dealt with. The solver is based on a finite-volume method to build a secondorder accurate discretization of the RANS equations. An arbitrary Lagrangian Eulerian ͑ALE͒ approach is used to deal with moving bodies. Ad-hoc deformation techniques have thus been developed to keep a body-fitted mesh during its motion.

At first, the filament is supposed to be rigid and only a fluid problem is solved. The flow is then established after 8 s and a Fourier analysis is performed on the lift effort: The dominant frequency of the fluid field yields f 1 f = 3.74 Hz, which is in good agreement with the Hübner result ͓39͔. Now, the filament is considered as a deformable beam. In principle, the physics of the FSI problem is the following: The beam is

Imposed motions in X = 0 Imposed motions in X = L r d T ͑0͒ = ͑0,0,0͒ r d T ͑L͒ = ͑͑L -2͒ +2/ ␣ sin͑␣͒ ,-2/ ␣ +2/ ␣ cos͑␣͒ ,0͒, ␣ ͓0,͔ r d Ј T ͑0͒ = ͑1,0,0͒ r d Ј T ͑L͒ = ͑cos͑␣͒ , -sin͑␣͒ ,0͒, ␣ ͓0,͔ d ͑0͒ =0 d ͑L͒ =0
submitted to the hydrodynamic ͑viscous and pressure͒ forces exerted by the fluid through its boundaries. These forces make the beam move. As a feedback, the motion of the beam modifies the boundary conditions of the flow, so closing the chain of causalities ͑see Fig. 6͒.

To handle this problem numerically, the CFD solver is coupled with the geometrically exact beam solver presented here. More precisely, an implicit/iterative algorithm ͑see Fig. 7͒ is used to coordinate the data exchanges between the fluid and solid solvers and obtain a stable coupling ͓40͔. This algorithm is structured into three hierarchical loops. The highest level corresponds to the temporal loop. The second one, denoted here as the FSI loop, aims to solve the nonlinearities of the fluid problem. Included in it, the third loop calls the beam solver. Thus, at each new FSI iteration ͑where the time is fixed͒, the structure problem is solved until convergence on the geometric nonlinearities using the current fluid loads ͑beam loop͒. It can be done at each FSI iteration since the CPU time of the structure resolution is negligible compared with that of the fluid resolution. The computed configuration of the beam is then updated. In particular, an adjustment of the fluid mesh is carried out to fit with the updated configuration of the beam. At that time, the core of the fluid solver is called to obtain As an illustration, Table 6 shows the residual norm for the structure problem during the iterations of the beam solver ͑at a given fluid nonlinear iteration of a given time step͒.

X(m ) -3 -2 -1 0 Y(m ) -1 0 1 (a) X(m ) -3 -2 -1 0 Z(m ) -4 -3 -2 -1 0 (b) Y (m) -1 0 1 Z(m) -4 -3 -2 -1 0 X ( m ) -3 -2 -1 0 Y (m) -1 0 1 Z (m) -4 -3 -2 -1 0 X Y Z (d) (a)( b) (c)( d)
One of the difficulties in fluid-structure interaction is due to the fact that the fluid and solid meshes do not generally match ͑cf. Ref. ͓41͔͒. In the present case, the dimensions of the two media models are not even the same. In fact, from the structure point of view, the beam is represented by its discretized 1D centroidal line, whereas the fluid views the structure by a set of boundary faces around this line ͑see Fig. 8͑a͒͒. As a consequence, the 3D contact forces ͑pressure and viscous͒ exerted by the fluid on the boundaries of the beam have to be converted into a set of nodal wrenches compatible with Kirchhoff kinematics. The fluid forces acting on the beam are initially described by elementary forces applied on each face of the fluid mesh in contact with the beam boundary ͑Fig. 8͑b͒͒. Due to possible overlaps of fluid faces on beam elements, the faces are split on each element into a set of subfaces. Then, for each subface ͑which belongs to a single beam element͒, the corresponding part of the contact force is applied as a localized one at the geometric center of the subface ͑Fig. 8͑c͒͒. In the end, we have, on each beam element, a set of localized 3D loads exerted onto the wetted boundaries of some intermediate cross sections ͑between the two nodes of each element͒͑ see Fig. 9͑a͒͒. Then, invoking Kirchhoff kinematics, these forces are first changed into localized generalized forces f ext ͑Fig. 9͑b͒͒. Next, these forces are discretized and assembled into the vector f ˜ext ͑of the form given by Eq. ͑19͒͒, which is applied at the two end nodes of the beam by invoking the nodal interpolation ͑Fig. 9͑c͒͒. This process ensures the global conservation of the forces, which is an important criterion to fulfill. Finally, when computing K ˜ext , only the beam configuration dependence is taken into account in the linearization of f ext . Indeed, the complete computation of K ˜ext would require the numerical Jacobian of the CFD solver to be computed, which is prohibitive ͑in terms of complexity and CPU time͒. This choice explains why we assumed the hydrodynamic wrench density ͑see Eq. ͑16͒͒ to be of the follower type. Note that even if f ˜ext is only considered as an imposed force ͑without computation K ˜ext ͒, the FSI coupling is also successful.

In the Hübner problem with a deformable beam, the time step is set at ⌬t =10 -3 s. The inflow leads to the production of vortices and so the filament begins to oscillate in the second mode. Figure 10 represents the deflection of the free end of the beam y B in time. This curve can be compared with the following results of Valdés Vázquez ͓42͔ and Hübner.

1. The free end evolution is very similar ͑the maximum amplitude is reached at 8 ϫ 10 -3 m͒. 2. The first oscillating frequency obtained by FFT around 3 Hz The coupled problem frequency f 1 c Ϸ 3 Hz is similar to the dominant frequency of the fluid ͑f 1 f = 3.74 Hz͒ and the second natural frequency of the structural part ͑f 2 s = 4.09 Hz͒ but it is lower, and so we can think that the first frequency of the structure ͑f 1 s = 0.61 Hz͒ influences the coupled problem.

From the beam point of view, the results obtained with the coarse or fine fluid meshes are very similar. However, with the fine mesh, the details of the fluid flow are better described ͑cf. Fig. 11͒. The number of beam finite elements can be reduced from 100 to 20 without a noticeable degradation of the beam motion. Finally, we recover the good numerical performance of the geometrically exact approach since a limited number of elements leads to a grid-independent solution ͑as a comparison, 123 elements were used in Ref. ͓42͔͒.

3D FSI Illustration: A Towed Cable in Water.

A 3D IFS simple test-case has been done in order to show the capabilities of the FSI code in 3D and to prepare the program to work on riser/ VIV problems: A deformable cable is fixed at an extremity, free at the other, and is partly dived into water ͑see Fig. 12͒. The following two different configurations were realized using the same coupling procedure described in Sec. 7.2:

1. one with a fixed cable in a fluid with constant velocity ͑see Fig. 13͒ 2. one with a towed cable in still water ͑see Fig. 14͒ In both configurations, the cable deforms itself a lot ͑see Fig. 13͒ and periodically oscillates. It is a consequence of the VIV. The position of nodes is updated using the same grid deformation technique used for the Hübner test-case 7.2 to recover a body-fitted mesh each time it is needed. Near the body, the cells are not deformed and, so, are always orthogonal to the cable ͑see Fig. 15͒, which is good for the quality of the results.

The two configurations are physically identical. However, the second one ͑where the cable is towed in still water ͑see Fig. 14͒ enables to use complex imposed motion at the cable extremity, for example, the cable could be slowly accelerated to reach the nominal velocity, which avoids a violent start͒.

Conclusion

In this article, we propose a geometrically exact finite element approach in order to solve the problem of nonlinear cable simulations in finite transformation with small strains. Two potential applications are envisaged by this work. The first is related to the prediction of motions and forces undergone by cables or other flexible pipes, which are often used in robotics for carrying energy to the actuators. The second concerns the interaction with a fluid flow and especially the VIV phenomenon for undersea cables, such as risers, that are used in oceanic engineering or electric cables subjected to wind. As is usual when dealing with very elongated structures, such as ropes or cables, the model of the structure is that of Kirchhoff beams. However, the proposed finite element approach differs from the previous work on that topic by the fact that it is directly inspired by the geometrically exact approach initiated by Simo and his successors for dealing with Reissner-Timoshenko beams. In fact, like in this modeling framework, no approximation of the finite rotations is assumed in all the analytical computations required by the finite element method, except those related to the unavoidable time and space discretizations. Here, no advanced work has been investigated on this numerical part, and classical time integration has been used. Beyond these considerations, it will be useful to design specific time integrators able to preserve energy and momentum of the cable. This could be done through the growing results of geometric integrators ͓43-46͔, extending them to the configuration space S 2 ϫ SO͑2͒. Based on this nonlinear Euler-Bernoulli dynamic model, the problem of cable dynamics with imposed motions at the ends is investigated through an augmented Lagrangian approach. The resulting dynamic beam solver is also coupled with a CFD solver that deals with the Reynolds-averaged Navier-Stokes equations. The methodologies are presented here and validated through a 2D benchmark test-case of VIV. An illustration dealing with a 3D case is also shown. These results are currently being extended to more realistic 3D tests related to a riser, where experimental data are available. This integration opens the way to analyze more accurately FSI configurations with high Reynolds number turbulent flow around complex geometry. The challenge to simulate such configurations is vital in the industrial domain ͑such as off-shore͒. Finally, as regards its application to robotics, the solver is currently coupled with a numerical code devoted to the detection of collisions and the simulation of the contacts of the cable with itself and with its environment through nonsmooth dynamic techniques.

n i n i+1 m i m i+1 node i node i+1 (c)

Appendix A: Dynamic and Internal Tangent Stiffness

Because only the stiffness operators are concerned in this Appendix, in what follows, all the fields subscripted by ⑀ are deduced from those nonsubscripted by substituting the field r by its perturbation r ⑀ = r + ⑀⌬r and the roll field by ⑀ = + ⑀⌬.

Fig. 1

 1 Fig. 1 Basic kinematics of a beam Fig. 2 Parameterization of the rotation field

  ͓0,1͔ = ഫ e=1 P ͓X e ,X e + l e ͔͑ 27͒ The first and second variations of the position field r are interpolated using Hermit polynomials while the first and second variations of the ͑iterative͒ roll angle ͑␦ , ⌬͒ are interpolated using Lagrange ones, i.e., for any X ͓X e , X e + l e ͔ ␦r͑X͒ = N I,e ͑X͒␦r I,e + M I,e ͑X͒␦r I,e Ј , ⌬r͑X͒ = N I,e ͑X͒⌬r I,e + M I,e ͑X͒⌬r I,e Ј ␦͑X͒ = L I,e ͑X͒␦ I,e , ⌬͑X͒ = L I,e ͑X͒⌬ I,e

1.

  Updating of the position, tangent, and curvature fieldsr ͑t+⌬t͒ ͑k+1͒ = r ͑t+⌬t͒ ͑k͒ + ⌬r, r ͑t+⌬t͒ Ј ͑k+1͒ = r ͑t+⌬t͒ Ј ͑k͒ + ⌬rЈ, r ͑t+⌬t͒ Љ ͑k+1͒ = r ͑t+⌬t͒ Љ ͑k͒

Fig. 3

 3 Fig. 3 Initial configuration of the six-axis "Puma robot"

Fig. 4 A

 4 Fig.4A cable fixed at the end of a moving Puma type manipulator: "a… planar projection in "e 1 ,e 2 …, "b… planar projection in "e 3 ,e 1 …, "c… planar projection in "e 2 ,e 3 …, and "d… 3D visualization

Fig. 5

 5 Fig. 5 Elastic filament fixed to a square rigid body in an incompressible flow

  Fig. 8 Models of the structure and fluid loads

Fig. 9

 9 Fig. 9 Transfer of fluid contact forces to external wrenches compatible with Kirchhoff beam kinematics

KFig

  Fig. 10 Hübner test-case results with the Bossak method "⌬t =10 -3 s… on the coarse mesh

Table 1 Joint time evolutions

 1

	j	j	j
	1	1 ͑t͒ = / 6 sin͑ 1 t͒	1 =
	20	0
	30	0
	4 5	4 ͑t͒ =-/ 2+ / 6 sin͑ 4 ͒ 5 ͑t͒ = / 6 sin͑ 5 t͒	1 = 5 =
	60	0
	70	0
	Table 2 Cable characteristics "SI units…
	Length	L =5 m
	Diameter	d = 0.015 m
	Lineic mass	= 5 kg m -1
	Stretching stiffness	EA=5ϫ 10 7 N
	Twisting stiffness	GI p =9ϫ 10 2 N.m 2
	Bending stiffness	EI a =5ϫ 10 2 N.m 2

Table 3 Imposed motions during the mounting phase

 3

"SI units…

Table 6 Typical convergence history of the beam solver

 6

	Iteration	0	12345
	Residual	4.22	2.8e -2	2.1e -2	1.1e -4	1.0e -4	5.3e -7

Referring to Riemannian geometry, any ⌳ u ͑v͒ defines on the unit sphere S

,t h e action of the parallel transport of the Levi-Civita connection along the geodesic ͑great circle͒ linking u S 2 with v S 2 .2 Given in the case of SO͑3͒ by exp͑v ˆ͒ = cos͑ʈvʈ͒I

+ ͑1 -cos͑ʈvʈ͒͒͑v / ʈvʈ͒

d / d⑀ ͉ ⑀=0 denotes the directional ͑Gateau͒ derivative.

In fact, one of the two integrations by parts is done; the second will be "replaced" by the polynomial interpolation ͑Hermit͒ of the finite element method.

Appendix B: Tangent Stiffness of Follower Loads and Generalized Constraint Forces

Tangent stiffness of follower loads