Geometrically exact Cosserat rods with Kelvin–Voigt type viscous damping

We present the derivation of a simple viscous damping model of Kelvin–Voigt type for geometrically exact Cosserat rods from three–dimensional continuum theory. Assuming a homogeneous and isotropic material, we obtain explicit formulas for the damping parameters of the model in terms of the well known stiffness parameters of the rod and the retardation time constants defined as the ratios of bulk and shear viscosities to the respective elastic moduli. We briefly discuss the range of validity of our damping model and illustrate its behaviour with a numerical example.


Introduction
Simulation models for computing the transient response of structural members to dynamic excitations should contain a good approach to account for dissipative effects in order to be useful in realistic applications.If the structure considered may be treated within the range of linear dynamics with small vibration amplitudes, there is a well established set of standard approaches, e.g.Rayleigh damping, or a more general modal damping ansatz, to add such effects on the level of discretized versions of linear elastic structural models (see e.g.Craig and Kurdila, 2006).In the case of geometrically exact structure models for rods and shells (Antman, 2005), such linear approaches are not applicable.Geometrically exact rods, in particular, have a wide range of applications in flexible multibody dynamics.We refer to the brief introduction given in ch.6 of Géradin and Cardona (2001) for a summary of the related work published before 2000, and to ch. 15 of Bauchau (2011) for a more recent account on this subject.Here the proper way to model viscous damping requires the inclusion of a frame-indifferent viscoelastic constitutive model into the continuum formulation of the structure model that is capable of dealing with large displacements and finite rotations (see Bauchau et al., 2008).

Viscous Kelvin-Voigt damping for Cosserat rods
In our recent work (Lang et al., 2011), we suggested the possibly simplest model of this kind to introduce viscous material damping in our quaternionic reformulation of Simo's dynamic continuum model for Cosserat rods (Simo, 1985).Following general considerations of Antman (2005) about the functional form of viscoelatic constitutive laws for Cosserat rods, we simply added viscous contributions, which we assumed to be proportional to the rates of the material strain measures U(s, t) and V(s, t) of the rod, to the material stress resultants F(s, t) and stress couples M(s, t), resulting in a constitutive model of Kelvin-Voigt type: A detailed presentation of the kinematical quantities and dynamic equilibrium equations of a Cosserat rod is given in Sect. 2 (see Figs. 1 and 2 for a compact summary).

3D momentum balance:
1 st Piola-Kirchhoff stress: with stiffness parameters given by the elastic moduli E and G and geometric parameters (area A, geometric moments I k ) of the cross section.In Lang et al. (2011) we assumed a similar structure for the matrices VF and VM , which determine the viscous response: The set of six effective viscosity parameters γ xx introduced in Eq. ( 3) represents the integrated cross-sectional viscous damping behaviour associated to the basic deformation modes (bending, twisting, transverse shearing and extension) of the rod, in the same way as the well known set of stiffness parameters given above determines the corresponding elastic response.

Effective damping parameter formulas
However, in Lang et al. (2011) the damping parameters γ xx remained undetermined w.r.t.their specific dependence on material and geometric properties.Considering the special case of homogeneous and isotropic material properties, they certainly cannot be independent, but rather should be mutually related in a similar way as the stiffness parameters of the rod in terms of two material parameters (E, G) and the geometrical quantities (A, I k ) associated to the cross section.
Assuming moderate curvature of the rod in its reference configuration, strains remaining small in its deformed configurations, strain rates that vary slowly compared to internal relaxation processes within the material, and a homogeneous and isotropic material, we will show that they are given by where ζ and η are the bulk and shear viscosities of a viscoelastic Kelvin-Voigt solid (Lemaitre and Chaboche, 1990) with elastic moduli G and E = 2G(1 + ν).While the viscous damping of the deformation modes of pure shear type is solely affected by shear viscosity η, extensional and bending deformations are both associated to normal stresses in the direction orthogonal to the cross section, which are damped by a specific combination of both bulk and shear viscosity that depends on the compressibilty of the material and may be interpreted as extensional viscosity parameter Introducing the retardation time constants τ S = η/G and τ B = ζ/K, which relate the viscosities η and ζ to the shear and bulk moduli G and 3K = E/(1 − 2ν), as well as the time constant relating extensional viscosity to Young's modulus, the formulas (4) may be rewritten equivalently as in terms of the stiffness parameters of the rod and the retardation time constants.Interesting special cases of Eq. ( 6) are the simplified expressions for completely compressible materials (ν = 0), and η E = 3η, τ E = τ S for incompressible materials (ν = 1 2 ).The relation η E /η = 3 between shear and extensional viscosity is well known as Trouton's ratio for incompressible Newtonian fluids (Trouton, 1906) and holds more generally for viscoelastic fluids in the limit of very small strain rates (Petrie, 2006)

Effective parameters modified by shear correction factors
It is well known that the stiffness parameters GA and GI 3 related to shearing type deformation modes systematically overestimate the actual stiffness of the structure for cross section geometries that display non-negligible warping.In the case of transverse shearing, this is accounted for via a modification of the corresponding stiffness parameter GA → GA α := GAκ α by introducing dimensionless shear correction factors κ α ≤ 1 depending on the cross section geometry (see Cowper, 1966;Gruttmann and Wagner, 2001).Likewise, the torsional rigidity C T = GJ T of a rod exactly equals GI 3 in the case of (annular) circular cross sections only, but is smaller than this value otherwise due to the presence of out-ofplane warping of cross sections.The replacement GI 3 → C T correcting this deficieny corresponds to the introduction of another dimensionless correction factor κ 3 = J T /I 3 ≤ 1 depending on the cross section geometry 1 which modifies the torsional stiffness according to the replacement rule GI 3 → GJ T = GI 3 κ 3 .Altogether the various shear corrections mentioned above yield the corrected set of stiffness parameter values 1 In the case of an elliptic cross section with half axes a and b, the area moments are given by I 1 = π 4 a 3 b and 1 in this case.Equality (κ 3 = 1) holds in the case of a circular cross section with a = b = r ⇒ I 1/2 = π 4 r 4 = 1 2 I 3 only.According to Nikolai's inequality C T ≤ 4GI 1 I 2 /(I 1 + I 2 ) the special case of an elliptic cross section maximes torsional rigidity among all asymmetric cross section geometries, and the value GI 3 = 2GI valid for circular cross sections provides the absolute maximum of torsional rigidity (Berdichevsky, 1981). 2 The stiffness parameters EA and EI α are not affected by shear warping effects.However, they already account for uniform lateral contraction, which is a simple specific type of in plane cross section warping.This topic is discussed further in Sect.3.4 below.
We argue that the analogously modified damping parameters associated to shearing type rod deformations likewise provide a corresponding improvement of the formulas (6), which accounts for the influence of cross section warping on effective viscous dissipation, such that the effective viscosity matrices VF and VM introduced in Eq. ( 3) may be rewritten as in terms of the effective stiffness matrices and retardation time constants given above.

Related work on viscoelastic rods
While there is a rather large number of articles considering various kinds of damping terms (also of Kelvin-Voigt type) added to linear Euler-Bernoulli or Timoshenko beam models (usually assumed to have a straight reference geometry), one hardly finds any work on viscous damping models for geometrically nonlinear beams or rods in the literature.One notable exception is Antman's work (2003), where a damping model as given by Eq. ( 1) with positive, but otherwise undetermined parameters (3) is suggested from a completely different, mathematically motivated viewpoint, namely: as a simple possibility to introduce dissipative terms (denoted as artificial viscosity) into the dynamic balance equations of a Cosserat rod, which constitute a nonlinear coupled hyperbolic system of PDEs (see also Weiss, 2002a), and thereby achieve a regularization effect in view of the possible formation of shock waves that might appear in the undamped hyperbolic equations.
The recent article of Abdel-Nasser and Shabana (2011) is another relevant work for our topic.By inserting a 3-D Kelvin-Voigt model into a geometrically nonlinear beam given in absolute nodal coordinate formulation (ANCF), the authors obtain a viscous damping model for such ANCF beams which (by construction) is closely related, but conceptually quite different from our approach proposed for Cosserat rods.Later we briefly discuss the relation of both damping models (see Sect. 4.3).We refer othwise to the article of Romero (2008) for a comparison of the geometrically exact and ANCF approaches to nonlinear rods.Mata et al. (2008) model the inelastic constitutive behaviour of composite beam structures under dynamic loading, using a Cosserat model as kinematical basis.However, they evaluate inelastic stresses by numerical integration of 3-D Piola-Kirchhoff stresses over 2-D discretizations of the local cross sections to obtain the stress resultants and couples of Simo's model.This differs from our approach aiming at a direct formulation of frame-indifferent inelastic constitutive laws in terms of F and M, as achieved e.g. by Simo et al. (1984) for viscoplastic rods.The viscous model proposed in Sect.3.2 of their paper is likewise of Kelvin-Voigt (KV) type, but formulated in terms of a vectorial strain measure related to the Biot strain (see also Sect.A2) and defined pointwise within the cross section.Moreover, they set up their model using only a single viscosity parameter.
Although there seems to be no further work on viscoelastic Cosserat rods made from solid material, viscoelastic flow in domains with rod-like geometries has been discussed in a number of articles.In his work on the coiling of viscous jets, Ribe (2004) presents a reduction of the three-dimensional Navier-Stokes equations to the dynamic equilibrium equations of a Kirchhoff/Love rod, endowed with Maxwell type constitutive equations for the viscous forces and moments which govern the finite resistance of the jet axis to stretching, bending and twisting.Although the derivation approach is different from ours, it represents its fluid-mechanical counterpart, as it likewise provides effective damping parameters3 as given in Eq. ( 4), in the special case of an incompressible viscous fluid (ν = 1 2 ) with extensional viscosity given by Trouton's relation η E = 3η, which in turn confirms our derivation of this special result.
A systematic derivation and mathematical investigation of viscous string and rod models in the context of Ribe's work is given by Panda et al. (2008) and Marheineke and Wegener (2009).Klar et al. (2009) and Arne et al. (2011) likewise use Ribe's Maxwell type constitutive law in their related work on the simulation of viscous fibers aiming at applications in the area of textile and nonwoven production.Lorenz et al. (2012) extend constitutive modelling for viscous strings by deriving an upper convected Maxwell model using mathematical methods of asymptotic analysis.
In the same context we finally mention the discrete modelling approach for viscous threads presented by Bergou et al. (2010), which extends earlier work of Bergou et al. (2008) on discrete elastic rods that, similar to our own approach as briefly presented in Linn et al. (2008) (see also Jung et al., 2011), relies on geometrically exact rod kinematics based on the discrete differential geometry of framed curves.

Overview of the remaining sections of the paper
After collecting a few basics of Cosserat rod theory in the following Sect.2, we proceed with our derivation of the formulas (4) in of a two-step procedure: in Sect. 3 we start with the derivation of the elastic (stored) energy function of a Cosserat rod, which is a quadratic functional of the terms ∆U(s, t) = U(s, t) − U 0 (s) and ∆V(s, t) = V(s, t) − V 0 (s) measuring the change of the strain measures w.r.t their reference values, from three-dimensional continuum theory.This sets the notational and conceptional framework for the subsequent derivation of the viscous part of our damping model given in Sect. 4 by an analogous procedure, which yields the dissipation function of a Cosserat rod introduced4 in Lang et al. (2011).The dissipation function (11), deduced from the three-dimensional (volumetric) continuum version of the dissipation function of a Kelvin-Voigt solid (Landau and Lifshitz, 1986;Lemaitre and Chaboche, 1990), corresponds to one half of the volumeintegrated viscous stress power of a rod-shaped Kelvin-Voigt solid, such that 2D v yields the rate at which the rod dissipates mechanical energy.
Having completed our derivation of the Kelvin-Voigt model, we proceed by a discussion of a seemingly straightforward, but, as it turns out, erroneous approach to derive the viscous parts of the forces and moments as given by Eq. (1) as resultants in analogy to the elastic counterparts.This shows that our energy-based approach to derive viscous damping is the proper one.After that, we briefly comment on the relation of our continuum model to the Kelvin-Voigt type model recently proposed by Abdel-Nasser and Shabana (2011) within their alternative ANCF approach to geometrically nonlinear rods, and conclude Sect. 4 by a short discussion of the validity of the Kelvin-Voigt model w.r.t. a more general viscoelastic model of generalized Maxwell type.
In Sect.5, we illustrate the behaviour of our viscous damping model (1) by some simple numerical experiments with a clamped cantilever beam subject to bending with large deflections.We conclude our article with a short summary.

Basic Cosserat rod theory
The configuration variables of a Cosserat rod (see Antman, 2005) are its centerline curve ϕ(s, t) = ϕ k (s, t) e k with cartesian component functions ϕ k (s, t) w.r.t. the fixed global ONB {e 1 , e 2 , e 3 } of Euclidian space and "moving frame" R(s, t) = a (k) (s, t) ⊗ e k ∈ SO(3) of orthonormal director vectors, both smooth functions of the curve parameter s and the time t, with the pair {a (1) , a (2) } of directors spanning the local cross sections with normals a (3) along the rod (see Fig. 1).

Material strain measures
The material strain measures associated to the configuration variables are given by (i) the components V k = a (k) •∂ s ϕ of the tangent vector in the local frame (i.e.: V = RT • ∂ s ϕ = V k e k ), with V 1 , V 2 measuring transverse shear deformation and V 3 measuring extensional dilatation, and (ii) the material Darboux vector U = RT •u = U k e k , obtained from its spatial counterpart u = U k a (k) governing the Frénet equations ∂ s a (k) = u×a (k) of the frame directors, with U 1 , U 2 measuring bending curvature w.r.t. the director axes {a (1) , a (2) }, and U 3 measuring torsional twist around the cross section normal.
In general, the reference configuration of the rod, given by its centerline ϕ 0 (s) and frame R0 (s) = a (k)  0 (s) ⊗ e k , may have non-zero curvature and twist (i.e.: U 0 0).However we may assume zero initial shear (V 01 = V 02 = 0), such that all cross sections of the reference configuration are orthogonal to the centerline tangent vector, which coincides with the cross section normal (i.e.: ∂ s ϕ 0 = a (3)  0 ⇒ V 03 = 1) if we choose the arc-length of the reference centerline as curve parameter s.

Dynamic equilibrium equations
The constitutive equations (1) -or more general ones of viscoelastic type (see ch. 8.2 in Antman, 2005) -are required to close the system of dynamic equilibrium equations (see Fig. 2) which has to be satisfied by the spatial stress resultants f = R • F and stress couples m = R • M with appropriate boundary conditions (see Simo, 1985).The inertial terms appearing on the r.h.s. of the equations of the balance of forces (linear momentum) ( 12) and the balance of moments (angular momentum) (13) depend parametrically on the local mass density ρ 0 (s) along the rod as well as on geometrical parameters of the local cross section (area A(s) and area moment tensor Ĵ(s, t) = R• Ĵ0 (s)• RT ) and contain the accelerations of the centerline positions ∂ 2 t ϕ(s, t) as well as the angular velocity vector ω(s, t), which is implicitely defined by the the temporal evolution equations ∂ t a (k) = ω × a (k) of the frame in close analogy to the Darboux vector, and its time derivative ∂ t ω(s, t) as dynamical variables (see Simo, 1985;Antman, 2005;Lang et al., 2011 for details).
Although we implemented Kelvin-Voigt type viscous damping given by Eq. (1) for our discrete5 Cosserat model formulated with unit quaternions as explained in detail by Lang et al. (2011) and investigated further in Lang and Arnold (2012) w.r.t.numerical aspects, we do not make use of this particular formulation here, as it is more practical to work with the directors associated to SO(3) frames for the vector-algebraic calculations which we have to carry out within our derivations of one-dimensional rod functionals from three-dimensional continuum formulation.

Spatial configurations of a Cosserat rod
Introducing cartesian coordinates (ξ 1 , ξ 2 ) w.r.t. the director basis {a (1)  0 (s), a (2) 0 (s)} of the cross section located at the centerline point ϕ 0 (s), the spatial positions of material points in the reference configuration of the rod are given by6 The positions of the same material points in the current (deformed) configuration are then given by in terms of the deformed centerline curve ϕ(s, t), the rotated orthonormal cross section basis vectors {a (1) (s, t), a (2) (s, t)}, the same pair of cartesian cross section coordinates (ξ 1 , ξ 2 ), and an additional displacement vector field w(ξ 1 , ξ 2 , s, t), which by definition describes the (in-plane and out-of-plane) warping deformations of the cross sections along the deformed rod.
The kinematic assumption that the cross sections of a rod remain plane and rigid in a configuration is equivalent to the assumption that the displacement field w vanishes identically.Although we will initially adhere to this very common assumption for rod models, we will later admit some specific form of in-plane deformation of cross sections -namely: a uniform lateral contraction -to correct a deficiency w.r.t.artificial in-plane normal stresses caused by the excessively rigid kinematical ansatz (15) with w ≡ 0.
For simplicity we assume the rod to be prismatic, such that all cross sections along the rod are identical, and the domain of the cartesian coordinates (ξ 1 , ξ 2 ) coincides with one fixed domain A ⊂ R 2 .As usual we choose the geometrical center of the domain A to coincide with the origin of R 2 recent collaboration with Schulze et al. (2012).We refer to the article of Zupan et al. (2009) for fundamental aspects of Cosserat rods with rotational d.o.f.represented by unit quaternions, as well as to the recent work (2012B) of the same authors discussing the undamped dynamics of quaternionic Cosserat rods with various time integration approaches.Appendix B contains additional remarks related to alternative discretization approaches and model variants.
such that ξ α A = 0 holds, where we introduced the shorthand notation f A := A f (ξ 1 , ξ 2 ) dξ 1 dξ 2 for the cross section integral of functions.In addition we choose the orientation of the orthonormal director pairs {a (1)  0 (s), a (2) 0 (s)} as well as {a (1) (s, t), a (2) (s, t)} to coincide with the principle geometrical axes of A, such that ξ 1 ξ 2 A = 0 holds.The quantities that characterize the geometric properties of the cross section in the Cosserat rod model are the cross section area A = 1 A , the two area moments With these definitions we obtain the centerline of the reference configuration as the average position ϕ 0 (s) = X A /A of all material points of the cross section located at fixed s.The same relation ϕ(s, t) = x A /A holds for deformed configurations provided that the warping field w(ξ 1 , ξ 2 , s, t) satisfies w A = 0.

The stored energy function of a Cosserat rod
In order to set the notational and conceptional framework for the derivation of the viscous part of our damping model, we first give a brief account of the derivation of its elastic part, i.e.: the stored energy function (10) of a Cosserat rod.Within this derivation we will encounter a variety of smallness assumptions w.r.t. the curvatures describing the reference geometry of the rod as well as the local strains occuring in its deformed configurations.In our subsequent derivation of the viscous dissipation function (11) we will use the same assumptions and thereby remain consistent with the derivation of the elastic part.

Three-dimensional strain measures
In the first step we compute the deformation gradient F = g k ⊗ G k , the right Cauchy-Green tensor Ĉ = FT • F and the Green-Lagrange strain tensor Ê = 1 2 ( Ĉ − Î) from the basis vectors G k = ∂ k X and g k = ∂ k x associated to the curvilinear coordinates of the rod configurations given by Eqs. ( 14) and ( 15), with The dual basis vectors G j and g j are defined by the relations G i • G j = δ i j and g i • g j = δ i j , respectively.Proceeding in this way we obtain the basis vectors of the reference configuration ( 14) as G α = a (α)  0 (s) and G 3 = a (3) 0 (s)+ξ α U 0α (s)a (α) 0 (s).Their duals may be computed from the general formula and the initial twist U 03 (s) of the reference configuration ( 14) influence the deviation of the dual vectors G k from the frame directors a (k)  0 (s) within the cross section.Both vectors coincide if the reference configuration of the rod is straight and untwisted (i.e.: U 0 = 0).We have approximate coincidence G k ≈ a (k)  0 (s) if curvature and twist of the reference configuration are suffi-ciently weak, in the sense that for the curvature radii given by R k = 1/|U 0k | the estimates |ξ α |/R 3 1 and |ξ α |/R β 1 ⇒ J 0 ≈ 1 hold throughout each cross section along the rod, such that all initial curvature radii R α are large compared to the cross section diameter.The geometric approximation J 0 (s) ≈ 1 will occur repeatedly and therefore play an important role in the derivation of the elastic energy and dissipation function of a Cosserat rod.To compute the deformation gradient we also need the basis vectors g α = a (α) (s, t) and g 3 = a (3) (s, t) + ξ α U α (s, t)a (α) (s, t) of the deformed configuration (15) with vanishing gradient of the warping vector field (∂ k w = 0).For the dual vectors g k one obtaines analogous expressions as those for the dual vectors G k given above, which we omit here.
For the special kinematical relations of a Cosserat rod, the deformation gradient F = g k ⊗ G k may be expressed in terms of a pseudo-polar decomposition (see Géradin and Cardona, 2001) by a factorization of the relative rotation Rrel (s, t) := R(s, t) • RT 0 (s) = a (k) (s, t) ⊗ a (k) 0 (s) connecting the moving frames of the reference and deformed configurations of the rod.The resulting formula depends on the absolute values of the curvatures of the reference configuration ( 14) through J 0 (s), and on the change of the strain measures of the Cosserat rod given by the difference vectors U(s, t) − U 0 (s) and V(s, t) which can be written more compactly7 in the form of a cartesian vector RT 0 Computing the right Cauchy-Green tensor Ĉ = FT • F with the deformation gradient given by Eq. ( 16) results in the following kinematically exact expression for the Green-Lagrange strain tensor: The approximate expression8 may be obtained from Eq. ( 18) by the geometric approximation J 0 ≈ 1 assumed to hold for the reference geometry and the additional assumption H 1 of a small material strain vector.Later we will make use of the approximate strain tensor ( 19), which is linear in the vector field H and therefore also in the change of the strain measures of the rod, to obtain the stored energy function (10), which then becomes a quadratic form in the change of the strain measures.Likewise we will use Eq. ( 19) to obtain an approximation of the strain rate ∂ t E in terms of the rate ∂ t H of the strain vector.

Validity of the small strain approximation
For deformed configurations of a slender rod one observes large displacements and rotations, but local strains remain small.To estimate the size of the strain tensor it is useful to compute its components E i j = a (i)  0 • ( Ê • a ( j) 0 ) w.r.t. the tensor basis a (i)  0 ⊗ a ( j)  0 obtained from the directors of the reference frame R0 (s).From Eqs. ( 18) and ( 19) we obtain identically vanishing in-plane components (E αβ = E βα ≡ 0), as well as the exact and approximate expressions of the components related to out-of-plane deformations of the local cross section.Introducing the the quantity |ξ| max := max (ξ 1 ,ξ 2 )∈A (|ξ 1 |, |ξ 2 |) to estimate the maximal linear extension of the cross section A, one may estimate the deviation of the determinant J 0 (s) from unity by |J 0 (s) as a coarse check of the validity of the approximation J 0 ≈ 1. Otherwise the smallness of the components of Ê is implied by the smallness of the components H k of the strain vector.According to Eq. ( 17) these components in turn become small if the change of the strain measures of the Cosserat rod is small, i.e. if the estimates For slender rods with moderately curved undeformed geometry these estimates are obviously easily satisfiable, except for extreme deformations of the rod that produce large curvatures or twists of the order of the inverse cross section diameter.In this case, the assumption of small strains obviously would be invalid.

Elastic constitutive behaviour of rods at small strains
If we assume the rod material to behave hyperelastically with a stored energy density function Ψ e ( Ê), a simple Taylor expansion argument9 shows that the behaviour of the energy density within the range of small strains may be well approximated by the quadratic function Ψ e ( Ê) ≈ 1 2 Ê : H : Ê, where H = ∂ 2 ÊΨ e ( 0) is the fourth order Hookean material tensor known from linear elasticity.This quadratic approximation yields a well defined frame-indifferent elastic energy density that is suitable for structure deformations at small local strains, but arbitrary large displacements and rotations, and therefore serves as a proper basis for the derivation of the stored energy function of a Cosserat rod.
The corresponding approximation of the stress-strain relation yields the 2nd Piola-Kirchhoff stress tensor Ŝ = ∂ ÊΨ e ( Ê) ≈ H : Ê for small strains.The 1st Piola-Kirchhoff stress tensor P, which is used to define the stress resultants and stress couples of the Cosserat rod model (see Simo, 1985, for details), is obtained by the transformation P = F • Ŝ using the deformation gradient, and the Cauchy stress tensor as the inverse Piola transformation σ = J −1 P• FT depending also on J = det( F).If we approximate the strain tensor Ê by Eq. ( 19) and consistently discard all terms that are of second order in H in accordance with our assumption of small strains, we have to use the approximation F ≈ Rrel (s) (which implies J ≈ 1) for the deformation gradient in all stress tensor transformations.This means that all pull back or push forward transformations are carried out approximately as simple relative rotations connecting corresponding frames R0 (s) and R(s, t) of the undeformed and deformed configurations of a Cosserat rod.Alltogether we obtain the approximate expressions10 for the various stress tensors, which are valid for the specific type of small strain assumptions encountered for Cosserat rods, as discussed above.
In the case of a homogeneous and isotropic material, the Hookean tensor acquires the special form of an isotropic fourth order tensor H SVK = λ Î ⊗ Î + 2µ I depending on two constant elastic moduli: the Lamé parameters λ and µ.Here Î and I are the second and fourth order identity tensors, which act on (symmetric) second order tensors Q by double contraction as where P = I − 1 3 Î ⊗ Î is the orthogonal projector on the subspace of traceless second order tensors, such that P : Ê = Ê − 1 3 Tr( Ê) Î yields the traceless (deviatoric) part of the strain tensor, and K = λ + 2 3 µ is the bulk modulus.

Modified strain tensor including lateral contraction
The stress-strain relation obtained from ( 22) is given by Inserting the approximate expressions ( 19) and ( 20) of the strain tensor and its components into Eq.( 23) yields the small strain approximation ŜSVK ≈ λH 3 Î + µ[H ⊗ a (3) 0 + a (3) 0 ⊗ H] of the stress tensor ŜSVK for Cosserat rods.The computation of the stress components w.r.t. the basis of R0 (s) directors yields normal stress components S αα ≈ λH 3 and S 33 ≈ (λ + 2µ)H 3 , and the shear stress components are given by S 12 = S 21 = 0 and S α3 = S 3α ≈ µH α , respectively.
As both elastic moduli λ = 2µν/(1−2ν) and λ+2µ = 2µ(1− ν)/(1 − 2ν) appearing in the expressions for the normal stress components, expressed in terms of the shear modulus µ = G and Poisson's ratio given by 2ν = λ/(λ + µ), diverge in the incompressible limit ν → 1 2 (just as the bulk modulus K = 2 3 1+ν 1−2ν G does), the normal stresses would become infinitely large whenever the normal strain E 33 ≈ H 3 becomes nonzero.This unphysical behaviour is a direct consequence of the kinematical assumption of plain and rigid cross section, which prevents any lateral contraction of the cross section in the case of a longitudinal extension.Therefore the assumption of a perfectly rigid cross section, as well as the expressions ( 18) and ( 19) derived under this assumption, are strictly compatible only with perfectly compressible materials (i.e.: in the special case ν = 0).
The standard procedure to fix this deficiency (see e.g.Weiss, 2002a) is based on the plausible requirement that all in-plane stress components S αβ (including the normal stresses S αα ), which for rods in practice are very small compared to the out of plain normal and shear stresses S α3 and S 33 , should vanish completely.This may be achieved by imposing a uniform lateral contraction with in-plane normal strain components E αα = −νE 33 upon the cross section.Although this procedure seems to be rather ad hoc, it may be justified by an asymptotic analysis 11 of the local strain field for rods, e.g. in the way as presented by Love (1927) in the paragraph §256 on the "Nature of the strain in a bent and 11 See Berdichevsky (1981) and ch.15 of Berdichevsky (2009) for a modern comprehensive analysis within Berdichevsky's variational asymptotic approach.
twisted rod" in ch.XVIII of his book.Following Love's analysis, we obtain the in-plane normal strains to leading order as E αα = ∂ α w α = −νE 33 with the additional requirement that E 12 = E 21 = ∂ 1 w 2 + ∂ 2 w 1 = 0, which determines the in-plane components w α of the the warping field w corresponding to the lateral contraction in terms of E 33 .
To obtain the modified value of E αα = −νE 33 one has to add an additional term −νE 33 a (α)  0 ⊗ a (α)  0 to the exact expression (18) of the strain tensor.Using the identity Î = a (k)  0 ⊗ a (k) 0 , we obtain the modified expression for the strain tensor, with E 33 ≈ H 3 as small strain approximation according to Eq. ( 19).Inserting the modified strain tensor (24) into the stress-strain equation of the Saint-Venant-Kirchhoff material with Tr( Ê ) = (1 − 2ν)E 33 ≈ (1 − 2ν)H 3 , and using the relation λ(1−2ν) = ν 1+ν E that relates the Lamé parameter λ to Young's modulus E, we obtain the following modified expression for the stress of a Cosserat rod: By construction, we now obtain vanishing in-plane stress components S 12 = S 21 = S αα ≡ 0, while the transverse shear stresses remain unaffected by the modification (i.e.: S α3 = S 3α ≈ GH α with G = µ).As 2G = E/(1 + ν), we likewise obtain the modified expression S 33 ≈ E H 3 for the normal stress component orthogonal to the cross section, which corresponds to the familiar expression from elementary linear beam theory, with Young's modulus E replacing λ + 2µ.

Elastic energy of a Cosserat rod
Next we demonstrate briefly that the modified expressions ( 24) and ( 25) immediately lead to the known stored energy function (10) mentioned in the introduction.
In the case of a hyperelastic material with an elastic (stored) energy density Ψ e the elastic potential energy of a body is given by the volume integral V 0 dV Ψ e of the energy density over the volume V 0 of the reference configuration of the body.In the case of a rod shaped body parametrized by the coordinates (ξ 1 , ξ 2 , s) of the reference configuration ( 14), the volume measure of V 0 is given by dV = J 0 dsdξ 1 dξ 2 , where J 0 is the Jacobian of the reference configuration (see Sect. 3.1).Using the geometric approximation J 0 ≈ 1, the stored energy function of a rod shaped body is obtained as the integral V 0 dV Ψ e ≈ L 0 ds Ψ e A of the density over the cross sections and along the centerline of the reference configuration of the rod.
In the special case of the energy density ( 22) this leads to the stored energy function W e = L 0 ds Ψ SVK ( Ê ) A , using the modified strain tensor Ê from Eq. ( 24).Applying our previously introduced approximations of small strains and small initial curvature, we obtain the approximate expression for the energy density.Its cross section integral Ψ SVK ( Ê ) A may be evaluated in terms of the integrals corresponding exactly to the stored energy function ( 10) with effective stiffness parameters given by Eq. ( 2).The subsequent introduction of shear correction factors (GA → GAκ α ) as well as the corresponding correction GI 3 → GJ T = GI 3 κ 3 of torsional rigidity 12 finally yields the stored energy function (10) with correspondingly modified effective stiffnesses as given by Eq. ( 7) (see also Sect.4.1 for a more detailed discussion of this point).

Kinetic energy and energy balance for Cosserat rods
In general, the kinetic energy of a body is given by the volume integral V 0 dV 1 2 ρ 0 v 2 , where ρ 0 (X) is the local mass 12 The correction of torsional rigidity accounts for the contribution of out-of-plane cross section warping in terms of a corresponding torsional stress function Φ(ξ 1 , ξ 2 ) and leads to an improved approximation of the strain and stress fields as well as the resulting elastic energy given by Eq. ( 10) compared to its 3-D volumetric counterpart.Similar arguments apply to an improved approximation of transverse shear strains and stresses as well as the associated part of the elastic energy density by accounting for additional contributions given by a corresponding pair of stress functions χ α (ξ 1 , ξ 2 ).The classical results obtained by St.-Venant are given in ch.XIV of Love's treatise (Love, 1927) (see also ch.II §16 in Landau and Lifshitz, 1986).They are contained as a special (and simplified) case within Berdichevsky's more comprehensive and modern treatment in terms of his method of variational asymptotic analysis applied to rods (see Berdichevsky, 1981, 1983and ch. 15 of Berdichevsky, 2009).Apart of Timoshenko's original treatment of shear correction factors, the article of Cowper (1966) is a classical reference on this subject, with correction factors obtained from pointwise (centroidal) and cross section averaged values of transverse shear stresses σ α3 (see also the discussions in ch.II, section 11 of Villagio, 1997 and section 2.1 of Simo et al., 1984).More recently an alternative approach based on energy balance as utilized e.g. in (Gruttmann and Wagner, 2001) and likewise fits to our considerations, is considered as standard due to superior results.However, the issue of correction factors for transverse shear in Timoshenko-type rod models is still subject of discussion and research activities (see e.g.Dong et al., 2010).density of the body in the reference volume, and v(X, t) = ∂ t x(X, t) is the velocity of the respective material point.Using the kinematic ansatz (15) with the geometric approximation J 0 ≈ 1, assuming a homogeneous mass density, and neglecting the contribution of cross section warping (w ≡ 0), we obtain the integral expression (α) ) 2 ] for the kinetic energy of the rod as a quadratic functional of the time derivatives of its kinematic variables.The rotatory part may be reformulated in terms of the material components Ω j = ω • a ( j) of the angular velocity vector ω = Ω j a ( j) of the rotating frame, which is implicitely defined by ∂ t a (k) = ω × a ( j) , by substituting k .This finally yields the familiar expression for the kinetic energy of a Cosserat rod as given in Lang et al. (2011) with Ω k expressed in quaternionic formulation.Altogether we obtain the approximation V 0 dV [ 1 2 ρ 0 v 2 + Ψ e ] ≈ W e + W k =: W m of the three-dimensional mechanical energy of a rod shaped body in terms of the corresponding sum of the kinetic and stored energy functions W k and W e of the Cosserat rod model as given above.In the absence of any dissipative effects, the mechanical energy must be conserved exactly in both the 3-D as well as the 1-D setting, such that the identities hold identically as a consequence of the respective balance equations for both the 3-D volumetric body and the 1-D rod.

Kelvin-Voigt damping for Cosserat rods
Now we have collected all technical prerequisites and approximate results that enable us to derive the dissipation function (11) of a Cosserat rod from a three-dimensional Kelvin-Voigt model in analogy to the derivation of the stored energy function (10) in a consistent way.
In Landau and Lifshitz (1986) (see ch. V §34) the dissipation function V dV 1 2 η i jkl εi j εkl is considered as an appropriate model of dissipative effects within a solid body near thermodynamic equilibrium, with constant fourth order tensor components η i jkl that are the viscous analogon of the components of the Hookean elasticity tensor.Transfering this ansatz to the formalism used in our paper, the dissipation function of Landau and Lifshitz (1986) becomes that of a Kelvin-Voigt solid as given in Lemaitre and Chaboche (1990) which is a quadratic form in the material strain rate ∂ t Ê defined as the time derivative of the Green-Lagrange strain tensor.The constant fourth order viscosity tensor V may be assumed to have the same symmetries as the Hookean tensor H, with its components depending on viscosity parameters in the same way as the components of H depend on elastic www.mech-sci.net/4/79/2013/Mech.Sci., 4, 79-96, 2013 moduli.The stress-strain relation of the Kelvin-Voigt model is given by Ŝ = H : Ê+V : ∂ t Ê, with the viscous stress 13 given by the term Ŝv := V : The dissipation function for a Cosserat rod results by inserting the rate ∂ t Ê of the modified strain tensor (24) into the dissipation density function Ψ KV of the Kelvin-Voigt model.We will compute this dissipation function explicitely in closed form for the special case of a homogeneous and isotropic material.In this special case, the viscosity tensor assumes the form depending on two constant parameters: bulk viscosity ζ and shear viscosity η.
To compute ∂ t Ê we use the expression (24) for the modified Green-Lagrange strain tensor of a Cosserat rod including the small strain approximation (19), with the result depending on the time derivative 0 (s) of the material strain vector with components written as a cartesian vector w.r.t. the global basis {e 1 , e 2 , e 3 }.
Inserting Eqs. ( 30) and (31) into the dissipation density function Ψ IKV (∂ t Ê ) = 1 2 ∂ t Ê : V IKV : ∂ t Ê of the isotropic Kelvin-Voigt model, analogous computational steps as those 13 Note that Ŝv : ∂ t Ê = 2Ψ KV (∂ t Ê) corresponds to the viscous stress power density, such that the integral P v (t) := 2 V dV Ψ KV (∂ t Ê) over the body volume yields the (time dependent) rate at which a Kelvin-Voigt solid dissipates mechanical energy under approximately isothermal conditions near thermodynamic equilibrium, (see ch. V §34 and §35) of Landau and Lifshitz, 1986).For a thorough discussion of the role of the dissipation function within the theory of small fluctuations near thermodynamic equilibrium from the viewpoint of statistical physics we refer to the the corresponding paragraphs in ch.XII in Landau and Lifshitz (1980) (in particular §121), as well as V. Berdichevsky's recent article 2003.In section VI of the latter, the author points out that a Kelvin-Voigt type constitutive relation holds also at finite strains, with the dissipative part governed by a fourth order viscosity tensor V[ Ê, ∂ t Ê] depending on the local strain and its rate.While a dependence of V on the invariants of ∂ t Ê in general prevents the existence of a dissipation function, the latter does indeed exist according to V.B.'s arguments if V = V[ Ê] is independent of the strain rate.This holds e.g. in the case of the Kelvin-Voigt limit of constitutive laws belonging to the class of finite linear viscoelasticity (Coleman and Noll, 1961) at sufficiently small strain rates (i.e.sufficiently slow deformations of a body).
done for the derivation of the stored energy Ψ SVK ( Ê ) in the previous subsection yield the expression with the extensional viscosity parameter η E as defined in Eq. ( 5) appearing as the prefactor14 of (∂ t H 3 ) 2 .The computation of the cross section integrals of the squared time derivatives (∂ t H k ) 2 yields the expressions from which we obtain the desired cross section integral of the dissipation density function: The dissipation function ( 11) of the Cosserat rod with diagonal damping coefficient matrices (3) and damping parameters ( 4) is then obtained as

Modification by shear correction factors
There is obviously a high degree of formal algebraic similarity in the derivations of the stored energy function (10) as presented in Sect.3.5 and the dissipation function (11) as presented above: both functionals result by inserting the specific strain tensor (24) of a Cosserat rod or respectively its rate (30) into a volume integral over the 3-D body domain of a density function defined as a quadratic form given by constant isotropic fourth order material tensors H and V, making use of the same geometric as well as "small strain" approximations implied by the specific kinematical ansatz (15) for the configurations of a Cosserat rod.The formal analogy in the derivation procedure leads to a dissipation density (32) that may be obtained from its elastic counterpart ( 27) by substituting viscosity parameters for corresponding elastic moduli (G → η, E → η E ) and strain rates for strain measures.
In the case of the stored energy function (10) the effective stiffness parameters (2) of the rod model are obtained from a derivation using a kinematical ansatz that completely neglects out-of-plane warping (i.e.: w 3 = 0 = ∂ k w 3 ) due to transverse shearing and twisting, but accounts for in-plane warping (i.e.: w α 0) in a simplified way by assuming a uniform lateral contraction (ULC) of the cross section according to the linear elastic theory (see Sect. 3.4).Softening effects due to out-of-plane warping are then accounted for by introducing shear correction factors 0 < κ j ≤ 1, which in the case of a homogeneous and isotropic material enter the model as multipliers A → A α = Aκ α and I 3 → J T = I 3 κ 3 of the area A and polar moment I 3 of the cross section and -according to the linear theory -depend solely on the cross section geometry.The modified stiffness constants (7) are obtained in combination with the elastic moduli G = µ and E, the latter appearing instead of λ + 2µ due to the enforcment of vanishing in-plane stresses by allowing for ULC according to Eq. ( 24).
Although the derivation of explicit formulas 15 for κ j is carried out for static boundary value problems, the same κ j , as well as the kinematic ansatz accounting for ULC, may be used for dynamic problems, due to the negligible influence of dynamic effects on the warping behaviour of cross sections, provided that the rod geometry is sufficiently slender.Therefore the geometric modifications A → A α = Aκ α and I 3 → J T = I 3 κ 3 , which have already been used to provide modified stiffness parameters ( 7) for an improved approximation of the 3-D (volumetric) elastic energy by the stored energy function (10) in the static as well as in the dynamic case, remain likewise valid to achieve a comparable improvement for the approximation of the 3-D integrated viscous stress power by the dissipation function ( 11), with modified damping parameters given by Eq. ( 8), leading to the modified expressions (9) for the effective viscosity matrices.
This completes our derivation of the Kelvin-Voigt type dissipation function of a Cosserat rod.Although the arguments given above would certainly benefit from a mathematical confirmation by rigorous (asymptotic) analysis, the latter is beyond the scope of this work.

An (erroneous) alternative derivation approach
The formulation of the Cosserat rod model given by Simo (1985) introduces spatial force and moment vectors f and m, usually denoted as stress resultants and stress couples, as the cross section integrals A of the traction forces of the 1st Piola-Kirchhoff stress tensor acting on the cross section area and the corresponding moments generated by the Piola-Kirchhoff tractions w.r.t. the cross section centroid, which are obtained by means of the "lever arm" vector ξ(s) = ξ α a (α) 0 (s).Both integrants may be expressed in terms of the 2nd Piola-Kirchhoff stress by means of the transformation P = F • Ŝ with the deformation gradient.In view of the small strain approximation P ≈ Rrel • Ŝ with Ŝ ≈ H : Ê discussed in Sect.3.3 we obtain the relations A 15 We refer to footnote 12 for a discussion of this issue.
connecting the spatial stress resultants f = R • F and stress couples m = R • M to their material counterparts rotated to the local reference frame R0 (s) = a k 0 (s) ⊗ e k .Expanding the material force and moment vectors w.r.t. the local ONB given by the reference frame R0 (s) as R0 (s) • F(s, t) = F k (s, t) a k 0 (s) and R0 (s) • M(s, t) = M k (s, t) a k 0 (s) yields their components in terms of the cross section integrals of the components of Ŝ w.r.t.this basis.To compute these components of the material force and moment vectors in closed form for the special case Ŝ = H SVK : Ê + V IKV : IKV with the approximate expressions ( 24) and ( 30) of the Green-Lagrange strain tensor and its rate and the constant isotropic material tensors H SVK = K Î ⊗ Î + 2G P and V IKV = ζ Î⊗ Î+2η P, we have to evaluate the cross section integrals with the stress components S α3 = GH α + η∂ t H α and Therefore ηE has to be interpreted as extensional viscosity, but obviously differs from the expression η E given in Eq. ( 5) and derived above by computing the dissipation function.Therefore the corresponding retardation time constant τE := ηE /E = 1 3 (τ B + 2τ S ), which is independent of the value of Poisson's ratio ν, likewise differs from the expression of the extensional retardation time τ E given in Eq. ( 6).Both expressions ηE and η E yield extensional viscosity as a combination of shear and bulk viscosity, but agree only in the special case ν = 0.The same assertion likewise holds for the corresponding retardation times, of course.However, only η E yields the correct incompressible limit η E → 3η for ν → 1 2 , while ηE tends to the smaller (and incorrect) value of 2η in this case.
The resulting expressions for the material force components are given by and the material moment components correspondingly by A comparison with the stiffness and damping parameters (2) and ( 6) entering the constitutive equations (1) shows that the derivation approach sketched above correctly yields all of the stiffness parameters as well as the damping parameters associated to transverse and torsional shear deformations.However, the damping parameters governed by normal stresses and extensional viscosity do not agree due to the appearance of τE instead of the correct time constant τ E .
The discrepancy between the results of both derivation approaches can be traced back to the fact that the integration www.mech-sci.net/4/79/2013/Mech.Sci., 4, 79-96, 2013 of the traction forces and their associated moments over the cross section fails to account for the non-vanishing contributions of the in-plane strain rates ∂ t E αα = −ν∂ t H 3 associated to uniform lateral contraction to the total energy dissipation of the rod.Paired with the corresponding viscous stress components to the dissipation function.As the cross section integrals given above involve only the stress components S α3 and S 33 , this additional source of damping is, by definition, not contained in the resulting formulas for the material force and moment components F j and M j obtained via this approach.However, this deficiency affects only the viscous part of the constitutive equations.The elastic part does not show any discrepancy, as the modified strain tensor (24) by construction provides vanishing in-plane elastic stress components (see Sect. 3.4), such that the stored energy function does not contain any contributions from non-vanishing in-plane elastic stresses to the elastic energy, and the cross section integrals of the traction forces and their moments yield all stiffness parameters correctly.
In summary, the considerations above suggest that, also in the case of more general viscoelastic constitutive laws, our approach to derive effective constitutive equations for Cosserat rods by computing the stored energy and dissipation functions is superior to the alternative approach based on a direct computation of the forces and moments as resultant cross section integrals of the traction forces and associated moments, as the latter yields an effective extensional viscosity which is systematically too small for partially compressible and incompressible solids (i.e.: 0 < ν ≤ 1 2 ).

ANCF beams with Kelvin-Voigt damping
In the recent article of Abdel-Nasser and Shabana (2011), a damping model for geometrically nonlinear beams given in the ANCF (absolute nodal coordinates) formulation has been proposed.The authors obtained their model by inserting the 3-D isotropic Kelvin-Voigt model as described above into their ANCF element ansatz.They used the Lamé parameters λ and µ as elastic moduli, and introduced corresponding viscosity parameters λ v and µ v , which they related to the elastic moduli by dissipation factors γ v1 and γ v2 .From the context it seems clear that in our notation γ v2 = τ S , such that µ v = Gτ S = η.Likewise we may identify , and the viscosities are related by the same relation as the elastic moduli (i.e.: If the ANCF ansatz chosen in Abdel-Nasser and Shabana (2011) handles lateral contraction effects correctly, both models should behave similar and yield similar simulation results.However, the appearance of the unmodified elastic moduli λ = 2µν/(1 − 2ν) and λ + 2µ = 2µ(1 − ν)/(1 − 2ν) in the element stiffness matrix (see Eq. 25 of the paper) indicates that the formulation chosen in Abdel-Nasser and Shabana (2011) may have problems in the case of incompressible materials (ν → 1 2 ).A clarifying investigation of this issue as well as a detailed comparison of both models remains to be done in future work.

Validity of the Kelvin-Voigt model
As remarked already in Landau and Lifshitz (1986), the modelling of viscous dissipation for solids by a dissipation function of Kelvin-Voigt type is valid only for relatively slow processes near thermodynamic equilibrium, which means that the temperature within the solid should be approximately constant, and the macroscopic velocities of the material particles of the solid should be sufficiently slow w.r.t. the time scale of all internal relaxation processes.
To illustrate and quantify this statement, we briefly discuss the one-dimensional example of a linear viscoelastic stress-strain relation σ Using a 1-D Kelvin-Voigt model σ KV (t) = Gε(t)+η ε(t) we obtain the simple expression σKV (ω) = [G + iηω] ε(ω), which approximates the generalized Maxwell model at sufficiently low frequencies with G = G ∞ and η = N j=1 G j τ j .The deviation between the generalized Maxwell model and its Kelvin-Voigt approximation may be estimated as This deviation may indeed become small, provided that the modulus | ε(ω)| of the strain spectrum, which appears as a weighting factor for the terms of the sum on the r.h.s., takes on non-vanishing values only at frequencies much smaller than those given by the discrete spectrum of the inverse relaxation times ω j = 1/τ j .The estimate given above also shows that in this case the Kelvin-Voigt model provides a low frequency approximation of second order accuracy.

Numerical examples
To illustrate the behaviour of our damping model, we show the results of numerical simulations of nonlinear vibrations of a cantilever beam in Fig. 3 obtained with the discrete Cosserat rod model presented in Lang et al. (2011).The parameters of the beam are: length L = 30 cm, quadratic cross-section area A = 1 × 1 cm 2 , mass density ρ = The beam is fully clamped at one end, the other end is initially pulled sideways by applying a force f L = Fe 1 of magnitude F = 0.05 N to the other end.The resulting initial deformation state in static equilibrium16 deviates far from the linear range of deformations governed by (infinitesimally) small displacements and rotations w.r.t. the reference configuration, while local strains are small in accordance with the constitutive assumptions.Starting from this initial equilibrium configuration, the beam is then released to vibrate transversally.The deformations of the beam shown in the inset of Fig. 3 are snapshots taken during the first half period of the oscillations which illustrate that in the initial phase of the oscillations substantial geometric nonlinearities are present.During the vibrations the beam remains in the plane of its initial deformation, such that all deformations are of plane bending type, and the extensional viscosity η E = Eτ E becomes the main influence for damping.
As expected, the plots of the transverse oscillation amplitude x(t) = e 1 • ϕ(L, t) recorded at the free end of the beam show an exponential dying out in the range of small amplitudes (linear regime).The deviations from the exponential envelope adapted to the linear regime that are observed during the initial phase clearly show the influence of geometric nonlinearity.The plots also suggest that damping becomes weaker in the nonlinear range.However, linear behaviour seems to start already with the fifth oscillation period, where the amplitude still has a large value of ≈ L/3.
This may be further analyzed by evaluating the logarithmic decrements δ k = ln(x(t k )/x(t k+1 )) recorded between succesive maxima x(t k ) of the amplitude as well as the corresponding damping ratios ζ k implicitely defined (see Craig and Kurdila, 2006, ch. 3.5, p. 75) k .The plots for the values of ζ k determined in this way are shown in the inset of Fig. 3.As expected, the ratios approach constant values in the linear regime, which scale as 1 : 2 : 4 proportional to the values of the time constant τ E used in the simulations.The simulations also show that the decrements become lower in the range of large amplitudes, which confirms the observation that the damping effect of our Kelvin-Voigt model is extenuated by the presence of geometrical nonlinearity.Nevertheless, ζ k still scales approximately proportional to τ E also in the nonlinear range.
To investigate the influence of a variation of the bending stiffness on the damping behaviour, an additional test with quadrupled Young's modulus E = 4 MPa was performed.In the corresponding amplitude plot shown in Fig. 3 the time axis of the plot with quadrupled E was streched twofold, such that the oscillations could be compared directly.After time stretching the (E = 4 MPa, τ E = 0.02 s) plot coincides with the (E = 1 MPa, τ E = 0.04 s) plot, surprisingly even throughout the whole nonlinear range.Since the oscillation period T of the four times stiffer (E = 4 MPa) beam is twice smaller than that of the softer (E = 1 MPa) beam, this suggests that the damping ratio varies proportional to the ratio τ E /T .Again this would be the expected behaviour in the linear regime, but is observed here in the nonlinear range as well.
For small amplitudes, the oscillation period may be estimated as T ≈ (2π/3.561)L 2 ρA/EI using the well known formula for the fundamental transverse vibration frequency of a cantilever beam obtained from Euler-Bernoulli theory (see Craig and Kurdila, 2006, ch. 13.2, Ex. 13.3, eq. 8).Inserting the parameters assumed above, we get T ≈ 1.81 s as an estimate, which correponds well to the time intervals of approximately 1.8 s between successive maxima shown in Fig. 3 that are also observed throughout the range of geometrically nonlinear deformations.For linear vibrations, damping ratio values ζ ≈ 1 correspond to a critical damping of the vibrating system, while values 0 < ζ 1 indicate a weak damping.According to that, the values ζ k observed in our experiments are in the range of weak to moderate damping, and are well approximated by the empirical formula ζ ≈ 1 π τ E /T .This provides a rough guideline for estimating the strenght of damping, or likewise an adjustment of the the retardation time τ E relative to the fundamental period T , if the Kelvin-Voigt model is utilized to provide artificial viscous damping in the sense of Antman (2003).According to this, a critical damping of transverse bending vibrations would be observed at a value of τ E ≈ πT .
Corresponding experiments for axial or torsional vibrations are limited to the range of small vibrations amplitudes, similar to the ones shown by Abdel-Nasser and Shabana (2011), as for large amplitudes one would inevitably induce buckling to bending deformations, such that all deformation modes would occur simultaneously, which greatly hampers a systematic investigation of different damping effects in the geometrically nonlinear range.Nevertheless, experiments at small amplitudes are helpful to determine the ranges of weak, moderate and critical damping for the respective deformation modes, quantifyable by explicit formulas similar to the one given above for the case of transverse vibrations.These could then be used e.g. to adjust damping of different deformation modes to experimental obervations.

Conclusions
In our paper we presented the derivation of a viscous Kelvin-Voigt type damping model for geometrically exact Cosserat rods.For homogeneous and isotropic materials, we obtained explicit formulas for the damping parameters given in terms of the stiffness parameters and retardation time constants, assuming moderate reference curvatures, small strains and sufficiently low strain rates.In numerical simulations of vibrations of a clamped cantilever beam we observed a slightly weakening influence of geometric nonlinearities on the damping of the oscillation amplitudes.We also found that the variation of retardation time and bending stiffness has a similar effect on the damping ratio as in the linear regime.In view of the limitations of the Kelvin-Voigt model w.r.t.higher frequencies it would be worthwile to develop more complex viscoelastic models (e.g. of generalized Maxwell type) for Cosserat rods.Our approach to derive Kelvin-Voigt damping for Cosserat rods may be helpful to obtain such models from three-dimensional continuum theory in an analogous way.

Appendix A Measuring 3-D strains and stresses for rods
From a mathematical point of view, the tensor Ĉ may be regarded as the fundamental quantity to decribe the shape of a body, as it corresponds to the metric which determines the shape up to rigid body motions, provided that certain integrability conditions (i.e.: the vanishing of the Riemann curvature tensor) are satisfied.Other strain measures may be obtained as invertible functions of Ĉ via its spectral decomposition.As a supplement to the brief discussion given in Sect.3.1, we mention a few alternatives to measure 3-D strains and stresses used elsewhere in connection with geometrically exact rod theory.

A1 The Biot strain and its approximation
In the case of small strains, the Biot strain tensor defined as ÊB := Û − Î, with the right stretch tensor Û given implicitely either by the polar decompostion F = Rpd • Û of the deformation gradient, or as Û = Ĉ1/2 in terms of the right Cauchy-Green tensor, is likewise an appropriate alternative choice of a frame-indifferent material strain measure.Due to the algebraic identity Ê = 1 2 ( Û2 − Î) = 1 2 ( Î+ Û)• ÊB the Biot and Green-Lagrange strains agree up to leading order for small strains, i.e.: Ê ≈ ÊB holds whenever Û ≈ Î.
One might argue that for small strains it is preferable to use ÊB as a strain measure, as it is linear in Û and therefore a first order quantity in terms of in the principal stretches, different from Ê, which is quadratic in Û.However, while (18) provides a kinematically exact expression for Î + 2 Ê = Ĉ = Û2 , a comparably simple closed form expression for Û itself is not available.In general the tensor Û has to be constructed via the spectral decomposition of Ĉ, which in 3-D cannot be expressed easily17 in closed form.
For special simplified problems, like the plane deformation of an extensible Kirchhoff rod as discussed by Irschik and Gerstmayr (2009) and Humer and Irschik (2011), it is possible to derive simple, kinematically exact closed form expressions18 for Û and Rpd by inspection of the deformation gradient.Also in the more general case of Ĉ given by Eq. ( 18) an analytical solution of the spectral problem is possible: by inspection N 3 := H×a (3)  0 /(H 2 1 +H 2 2 ) 1/2 is found to be one of its eigenvectors, with eigenvalue λ 2 3 = 1.The remaining 2-D spectral problem may then be solved analytically by a Jacobi rotation which diagonalizes the matrix representing Ĉ w.r.t. the ONB in the plane orthogonal to H × a (3)  0 given by a (3)  0 and the unit vector along the direction of the projection a (3)  0 × (H × a (3) 0 ) = H α a (α) 0 of the material strain vector H onto the local reference cross section.The resulting analytical formulas19 for the two eigenvalues λ 2 ± and orthonormal eigenvectors N 1/2 of which we present below without providing further details of their derivation, are given by: 0 , with H := H/J 0 , and the angle φ given implicitely by H2 1 + H2 2 cos(2φ) + ( H3 + H 2 /2) sin(2φ) = 0 .
They provide the spectral decomposition Ĉ = 3 k=1 λ 2 k N k ⊗N k of the right CG tensor (see Gurtin, 1981, ch. I and II), and the closed form expression ÊB = 3 k=1 (λ k −1) N k ⊗ N k of the Biot strain tensor, as Û = Ĉ1/2 .These considerations confirm that, although a kinematically exact closed form expression of ÊB for deformed configurations of a Cosserat rod (H 0) may be derived in this way, it consists of algebraically rather complicated expressions in terms of the vector H/J 0 and its components, compared to the relatively simple formula (18) for the Green-Lagrange strain.Otherwise, it is straightforward to show that Û ≈ Î+ 1 2J 0 H ⊗ a (3) 0 + a (3)  0 ⊗ H provides an approximate expression for the right stretch tensor of leading order in H/J 0 , as its square agrees with the exact expression for Ĉ up to terms of order O(H 2 /J 2 0 ).Therefore, we obtain ÊB ≈ 1 2J 0 H ⊗ a (3) 0 + a (3)  0 ⊗ H as an approximate expression for the Biot strain, which reduces to Eq. ( 19) for J 0 ≈ 1 and in this way provides an alternative interpretation of Eq. ( 19).Within the same order we may use Rpd (ξ 1 , ξ 2 , s, t) ≈ Rrel (s, t) to approximate the rotational part of the polar decomposition of F.
A2 Relation of the material strain vector to the Biot strain Following Kapania and Li (2003), Mata et al. (2007Mata et al. ( , 2008) ) use the spatial vector quantity (k)   with F given by a kinematically exact expression for the deformation gradient of a Cosserat rod equivalent to Eq. ( 16) to measure the strain at the individual points of a cross section.Its material counterpart J −1 0 RT 0 • H = J −1 0 H k e k as well as objective rates of both vector quantities are then used by these authors to formulate inelastic constitutive laws for their rod model on the 3-D level, which are required for a subsequent numerical evaluation of the spatial stress resultants and couples of the rod in its deformed configurations by numerical integration over the cross section areas.
Following our discussion of the Biot strain and its approximation given above, one recognizes that the strain measure used by Mata et al. (2008) likewise may be interpreted in terms of an approximation of the Biot strain via F − Rrel ≈ F − Rpd = Rpd • ÊB ≈ Rrel • ÊB .
Using F − Rpd as a strain measure is directly related to the geometric idea to quantify the strains caused by the deformation of a body by the deviation of a deformation mapping to a rigid body motion, as discussed by Chao et al. (2010).For a given deformation gradient F with positive determinant, this deviation may be measured by the distance of F to the group SO(3) of proper rotations defined as min R∈SO(3) F − R F , where • F denotes the Frobenius norm.It can be shown that the minimum is actually reached for the unique rotation R = Rpd provided implicitely by the polar decomposition of F, such that min R∈SO(3) F − R F = Rpd • ( Û − Î) F = ÊB F holds due to the invariance of the norm under rotations.Altogether these considerations, combined with the approximation Rpd ≈ Rrel , provide a geometric interpretation for the strain measure considered by Mata et al. (2008) and its relation to the Biot strain.

A3 The Biot stress and its approximation
In some works dealing with geometrically exact rods, e.g. in the articles of Irschik and Gerstmayr (2009) and Humer and Irschik (2011), 3-D stress distributions within cross sections are analyzed in terms of the (unsymmetric) Biot stress tensor TB := RT pd • P = Û • Ŝ, which is related to the (true) Cauchy stress σ via the co-rotational stress tensor RT pd • σ • Rpd = J −1 TB • Û.The Biot stress tensor TB as well as its symmetric part T(s) B := 1 2 ( TB + TT B ) are both work-conjugate stresses related to the Biot strain ÊB , as TB − T(s) B : δ ÊB = 0 holds, such that both yield identical virtual work expressions.
Small strain approximations of these stress quantities are obtained by substituting Û ≈ Î (implying F ≈ Rpd and J ≈ 1) into the various transformation identities for the stresses as given above.This yields the set of approximate relations TB ≈ RT pd • σ • Rpd ≈ T(s) B ≈ Ŝ, which are valid to leading order, analogous to the approximate relations ÊB ≈ Ê for the corresponding strain quantities.The approximate stress relations ( 21) are obtained by the additional approximation F • Û−1 = Rpd ≈ Rrel , likewise valid to the same order, which effectively amounts to applying the approximation F ≈ Rrel (implying J ≈ 1) within all transformations of stress tensors.
In summary, due to the assumption of small strains, the Biot and 2 nd Piola-Kichhoff stress tensors approximately coincide to leading order (i.e.: TB ≈ Ŝ), such that both stresses approximately correspond to the co-rotational stress tensor given by the components of the Cauchy stress (i.e.: TB ≈ RT rel • σ • Rrel ≈ Ŝ) w.r.t. the approximate material basis (i.e.: G k ≈ a (k)  0 ≈ G k ) given by the reference frames R0 (s).

Discretizations of the Kelvin-Voigt model
Our recent articles (Lang et al., 2011;Lang and Arnold, 2012) provide one concrete example of an implementation of a discrete version of our constitutive model (1), as an integral part of (and taylored to) our specific continuum formulation of the Cosserat rod model using unimodular quaternions, our specific spatial discretization approach -finite differences for the centerline, finite quotients for the quaternion field, both on a staggered grid -applied on the level of the stored energy (10), kinetic energy (see Sect. 3.6) and the dissipation function (11), the specific formulation of the resulting semidiscrete system as a first order DAE or ODE (depending on the kind of internal kinematical constraints and their treatment), and the class of time integration methods we choose to solve the semidiscrete equations for various initial-boundary value problems.
Our treatment differs substantially from other approaches as discussed e.g. in the textbooks (Géradin andCardona, 2001 andBauchau, 2011) or the article (Bauchau et al., 2008) already mentioned in the introduction, which mainly use finite elements (of first or higher order) for the spatial discretization, but again differ among each other in the treatment of rotational variables and the related interpolation strategy.
In addition, other model variants for geometrically nonlinear rods or beams exist, like the already mentioned ANCF approach used by Abdel-Nasser and Shabana (2011), or the recent approach of dynamic splines investigated by Theetten et al. (2008) and Valentini and Pennestri (2011), where geometrically exact extensible Kirchhoff rods, which require only a single angle variable to account for twisting, are desribed using computer-aided geometrical design functions, very similar to the usage of cubic Hermite splines on the element level as employed by Weiss (2002b).
In view of the great variety of discretization approaches applied to different geometrically exact rod models, a corresponding discussion of discrete versions of our Kelvin-Voigt model (1) for each variant is clearly beyond the scope of this article.In general, any implementation may be obtained most easily by a semidiscrete approach in terms of material strain quantities as used in Eq. (1).In this way one circumvents the technically rather complicated issue of constructing (and implementing) objective strain rates, which for the discretized material strain measures U h and V h of a Cosserat rod are given by simple partial time derivatives ∂ t U h and ∂ t V h .An adaption of Eq. ( 1) for the dynamic spline model mentioned above is obtained by setting the transverse shear strains and their rates to zero (V α = 0 = ∂ t V α ), such that V 3 = ∂ s ϕ(s, t) remains as the measure for elongational strain.

Figure 2 .
Figure 2. Dynamic equilibrium equations of a Cosserat rod.
e.: a Prony series) of a generalized Maxwell model.By Fourier transformation we obtain the relation σ(ω) = Ĝ(ω) ε(ω) in the frequency domain, where the real and imaginary parts of the complex modulus function Ĝ(ω) = G ∞ + N j=1 G j iτ j ω 1+iτ j ωmodel the frequency dependent stiffness and damping properties of the material.

Figure 3 .
Figure 3. Damped non-linear bending vibrations of a clamped cantilever beam (see text for further details).