International Journal of Differential Equations

International Journal of Differential Equations / 2012 / Article
Special Issue

Differential Equations with Applications to Industry

View this Special Issue

Research Article | Open Access

Volume 2012 |Article ID 718308 |

Thorsten Fütterer, Axel Klar, Raimund Wegener, "An Energy Conserving Numerical Scheme for the Dynamics of Hyperelastic Rods", International Journal of Differential Equations, vol. 2012, Article ID 718308, 15 pages, 2012.

An Energy Conserving Numerical Scheme for the Dynamics of Hyperelastic Rods

Academic Editor: Mapundi Banda
Received12 Jan 2012
Revised21 Jun 2012
Accepted05 Jul 2012
Published26 Aug 2012


A numerical method for special Cosserat rods based on Antman's description Antman, 2005 is developed for hyperelastic materials and potential forces. This method preserves the relevant properties of the underlying PDE system, namely, the orthonormality of the directors and the conservation of the energy.

1. Introduction

Elastic rods are considered in many fields of science and technology; see, for example, [13]. For the simulation of their dynamics, the correct description of the local and global physical properties and reasonable computation times are the essential requirements.

Over the past years many different approaches to the numerical simulation of elastic rods have been developed; see, for example, [49] for different approaches involving finite element methods, finite difference methods and discrete mechanics.

In this paper we use the description of the rod as a special one-dimensional Cosserat continuum following the formulation in Antman [10]. This is a geometrically correct one-dimensional description based on partial differential equations. This type of modeling fulfills the above requirements with respect to simplicity and correctness. Numerical schemes for the special case of Kirchhoff equations have been developed by Pai [1] (stationary) and Weber et al. [3] (instationary). Energy conservation and the directors’ orthonormality are not strictly fulfilled in these schemes. Schemes for structural mechanics problems conserving energy and further invariants of the equation are developed, investigated, and applied, for example, in [1113].

We formulate the kinematic and dynamic equations of the Cosserat rod theory in the so-called director-basis. The rod’s curve 𝗋 and the outer potential force 𝖿 are the only fields described in an external Cartesian basis. For technical reasons we use a representation of the rotational group by unit-quaternions. We note that the first use of quaternions in geometrically-exact rod models was in [14], see also [15]. In Section 2 we present an overview of the model resulting in the following system:𝜕𝑡𝗋=𝖣1(𝗊)𝗉,𝜕𝑡𝗊=𝛀(𝗐)𝗊,𝜕𝑡𝗏=𝜕𝑠𝗉+𝗎×𝗉+𝗏×𝗐,𝜕𝑡𝗎=𝜕𝑠𝗐+𝗎×𝗐,𝜕𝑡𝗉=1(𝜌𝐴)𝜕𝑠𝗇+1(𝜌𝐴)𝗎×𝗇+𝗉×𝗐+1(𝜌𝐴)𝖣(𝗊)𝖿,𝜕𝑡𝗐=(𝜌𝖩)1𝜕𝑠𝗆+(𝜌𝖩)1(𝗎×𝗆+𝗏×𝗇+((𝜌𝖩)𝗐)×𝗐).(1.1) Here, 𝑡 is the time and 𝑠 the parameter determining a material cross section of the rod. The unit-quaternions 𝗊 and the associated orthogonal matrix 𝖣(𝗊) describe the transformation between the fixed external basis and the director-basis. The vector fields 𝗏, 𝗎, 𝗉, and 𝗐 are the tangent of the curve, the generalized curvature, the velocity, and the angular velocity, respectively. Moreover, (𝜌𝐴) is the rod’s line density and the positive definite matrix (𝜌𝖩) is defined by the moments of inertia of its cross sections. In this paper the contact force 𝗇 and the contact couple 𝗆 are specified by a hyperelastic material law. In Section 3 we introduce the concept of energy as a constant of motion. In Section 4 we develop a straightforward finite difference scheme for the above equations with appropriate boundary conditions. The scheme conserves the energy and the orthonormality of the directors. In Section 5 the method is investigated for several examples using Timoshenko’s material law.

2. Model

Following Antman [10], a special Cosserat rod in the three-dimensional Euclidian space 𝔼3 is geometrically characterized by three vector-valued functions 𝐫,𝐝1,𝐝2(𝑠𝑎,𝑠𝑏)×𝔼3. The parameter 𝑠(𝑠𝑎,𝑠𝑏) identifies a material cross section (material point) of the rod, 𝐫(𝑠,𝑡) characterizes the position of this cross section at time 𝑡. The derivatives of the curve 𝐫 with respect to 𝑡 and 𝑠, 𝐩=𝜕𝑡𝐫,𝐯=𝜕𝑠𝐫,(2.1) are the velocity and the tangent field. The orthonormal directors 𝐝1,𝐝2 characterize the orientation of the cross sections. Defining 𝐝3=𝐝1×𝐝2, we get a local orthonormal basis at all material points. Due to the orthonormality of the directors there exist vector-valued functions 𝐮 and 𝐰 such that 𝜕𝑡𝐝𝑘=𝐰×𝐝𝑘,𝜕𝑠𝐝𝑘=𝐮×𝐝𝑘(2.2) for 𝑘{1,2,3}. We call 𝐰 the angular velocity and 𝐮 the generalized curvature. The definitions of the vector fields 𝐩, 𝐯, 𝐰, and 𝐮 imply the compatibility conditions: 𝜕𝑡𝐯=𝜕𝑠𝐩,𝜕𝑡𝐮=𝜕𝑠𝐰+𝐰×𝐮,(2.3) that complete the kinematic equations of the the special Cosserat rod theory in an invariant notation.

To rewrite these kinematic equations in an appropriate basis, we decompose an arbitrary vector field 𝐱 of our rod theory in the director-basis 𝐝𝑘 as well as in a fixed external Cartesian basis 𝐞𝑘, that is, 𝐱=3𝑘=1𝑥𝑘𝐝𝑘=3𝑘=1𝑥𝑘𝐞𝑘𝔼3. The corresponding component triple, 𝗑=𝑥1,𝑥2,𝑥33,𝗑=𝑥1,𝑥2,𝑥33,(2.4) are strictly to distinguish from the vector field 𝐱𝔼3. The component triples of its derivatives with respect to 𝑡 and 𝑠 in the director-basis are 𝜕𝑡𝐱𝐝𝑘𝑘=1,2,3=𝜕𝑡𝗑+𝗐×𝗑,𝜕𝑠𝐱𝐝𝑘𝑘=1,2,3=𝜕𝑠𝗑+𝗎×𝗑.(2.5) The transformation between the director-basis and the external basis is given by an orthonormal matrix 𝖣 with components 𝐷𝑖𝑗=𝐝𝑖𝐞𝑗. We use a representation of 𝖣 in unit-quaternions 𝗊=(𝑞0,𝑞1,𝑞2,𝑞3)4, cf. [16], 𝖣=𝖣(𝗊)=𝑞20+𝑞21𝑞22𝑞232𝑞1𝑞2𝑞0𝑞32𝑞1𝑞3+𝑞0𝑞22𝑞1𝑞2+𝑞0𝑞3𝑞20𝑞21+𝑞22𝑞232𝑞2𝑞3𝑞0𝑞12𝑞1𝑞3𝑞0𝑞22𝑞2𝑞3+𝑞0𝑞1𝑞20𝑞21𝑞22+𝑞23.(2.6) Here, the orthonormality of the matrix 𝖣(𝗊) is guaranteed by 𝑞=1. For an arbitrary 𝐱𝔼3 with 𝐱=3𝑘=1𝑥𝑘𝐝𝑘=3𝑘=1𝑥𝑘𝐞𝑘 we obtain 𝗑=𝖣(𝗊)𝗑. Instead of the kinematic equation for the directors, we formulate an equivalent equation for the quaternions, cf. [16], 𝜕𝑡𝗊=120𝑤1𝑤2𝑤3𝑤10𝑤3𝑤2𝑤2𝑤30𝑤1𝑤3𝑤2𝑤10𝗊=𝛀(𝗐)𝗊.(2.7) Initializing (2.7) with the unit-quaternions, the skew symmetry of Ω(𝗐) guarantees the preservation of the norm of the quaternions in time, that is, the orthonormality of 𝖣(𝗊) and hence the orthonormality of the director-basis. Using the presented formalism we obtain our final version of the kinematic equations of the rod: 𝜕𝑡𝗋=𝖣(𝗊)1𝗉,𝜕𝑡𝗊=𝛀(𝗐)𝗊,𝜕𝑡𝗏=𝜕𝑠𝗉+𝗎×𝗉+𝗏×𝗐,𝜕𝑡𝗎=𝜕𝑠𝗐+𝗎×𝗐.(2.8)

The balance laws for momentum and angular momentum yield the dynamic equations of the Cosserat rod theory, cf. [10], (𝜌𝐴)𝜕𝑡𝐩=𝜕𝑠𝐧+𝐟,𝛼,𝛽=1,2𝜌J𝛼𝛽𝜕𝑡𝐝𝛼×𝜕𝑡𝐝𝛽=𝜕𝑠𝐦+𝐯×𝐧+𝐥.(2.9) Here, (𝜌𝐴) is the line-density and (𝜌J𝛼𝛽) are the moments of inertia. These quantities are time independent, since they are defined in the reference configuration as Lagrangian quantities. We assume that they are constant with respect to 𝑠, too. The body force line density 𝐟 has to be specified in the applications. We assume 𝐥=𝟎 for the corresponding body couple density. The contact force 𝐧 and contact couple 𝐦 have to be defined by material laws. Usually, this is done in the director-basis for reasons of objectivity. Thus, we decompose also the dynamical equations in the director-basis: 𝜕𝑡𝗉=1(𝜌𝐴)𝜕𝑠𝗇+1(𝜌𝐴)𝗎×𝗇+𝗉×𝗐+1(𝜌𝐴)𝖣𝖿,𝜕𝑡𝗐=(𝜌𝖩)1𝜕𝑠𝗆+(𝜌𝖩)1(𝗎×𝗆+𝗏×𝗇+((𝜌𝐽)𝗐)×𝗐).(2.10) The positive definite matrix (𝜌𝖩) is given by the moments of inertia: (𝜌𝖩)=𝜌𝐽22𝜌𝐽120𝜌𝐽21𝜌𝐽11000𝜌𝐽11+𝜌𝐽22.(2.11) The body force, for example, gravity, is obviously defined in the external basis. Thus, (2.10) are coupled to the kinematic equation for the quaternions.

Remark 2.1. The special Cosserat rod theory describes the angular momentum as a linear function of the angular velocity. The choice of the representation of the vector fields in the director-basis leads to the time independent matrix (𝜌𝖩) characterizing this linear dependence. Besides the proper formulation of the material laws, this time independence is one of the major advantages of the choice of the director-basis.

In this paper we restrict ourselves to hyperelastic materials. That means there exists a potential 6(𝗏,𝗎)𝜓(𝗏,𝗎), such that 𝗇=𝜕𝗏𝜓,𝗆=𝜕𝗎𝜓.(2.12) Moreover, we assume that only potential forces act on the rod. Thus there exists a function 3𝗋𝑉(𝗋), such that 𝖿=𝜕𝗋𝑉.(2.13)

Remark 2.2. The more general class of elastic materials are materials where 𝗇 and 𝗆 are functions of the so-called strain variables 𝗏 and 𝗎. These functions may also depend explicitly on 𝑠.

The kinematic equations (2.8) and the dynamic equations (2.10) together with the restriction to hyperelastic materials and potential forces constitute our rod theory, see also (1.1). We consider system (1.1) with two types of boundary conditions defining a fixed or a free end. For simplicity we restrict our description to a fixed end at 𝑠=𝑠𝑎𝗉𝑠𝑎,𝑡=0,𝗐𝑠𝑎,𝑡=0,(2.14) and a free end at 𝑠=𝑠𝑏𝗇𝑠𝑏,𝑡=0,𝗆𝑠𝑏,𝑡=0.(2.15)

The presentation of the energy conserving numerical algorithm in Section 4 deals with the above general class of rods. For the numerical examples in Section 5 we specify the rod’s geometry, a hyperelastic material law, and the potential forces. We consider a homogeneous cylinder with diameter 𝑑>0, cross section area 𝐴=𝜋/4𝑑2, and moment of inertia 𝐼=𝜋/64𝑑4. In this case, the matrix of inertia is (𝜌𝖩)=(𝜌𝐼)diag(1,1,2).(2.16) We use the material law of Timoshenko [17] for Poisson number 𝜇=1/2: 𝜓=12(𝐸𝐴)13𝑣21+13𝑣22+𝑣312+12(𝐸𝐼)𝑢21+𝑢22+23𝑢23,(2.17) where 𝐸 is Young’s modulus. Additionally, we restrict to gravitational forces, that is, 𝑉=(𝜌𝐴)𝑔𝖾𝑔𝗋,(2.18) where 𝑔 is the gravitational constant and 𝖾𝑔 is the direction of gravity in the external Cartesian basis.

Remark 2.3. For hyperelastic materials (1.1) is an inhomogeneous hyperbolic system. In the special case of Timoshenko for a homogeneous cylinder the hyperbolic part is linear with eigenvalues 0 (sevenfold), ±𝑐 (threefold), and ±𝑐/3 (threefold), where 𝑐=𝐸/𝜌 is the speed of sound. Computing the eigenvectors one can easily show that the fixed and free end boundary condition correctly handle the characteristic variables 𝗋 and 𝗊 corresponding to the eigenvalue 0. With respect to the remaining variables we do not prescribe characteristic variables, but the correct number of variables on both sides of the rod.

3. Energy as a Constant of Motion

The system (1.1) for the state vector 𝝓=(𝗋,𝗊,𝗏,𝗎,𝗉,𝗐)19 can be written in the general form of a conservation law as 𝜕𝑡𝝓=𝜕𝑠𝐟(𝝓)+𝐡(𝝓)(3.1) with flux-function 𝐟(𝝓)=(𝟢,𝟢,𝗉,𝗐,(1/(𝜌𝐴))𝗇,(𝜌𝖩)1𝗆) and source term 𝐡(𝝓) that is easy to identify from (1.1). We introduce the energy density: 𝜀(𝝓)=12(𝜌𝐴)𝗉2+12𝗐(𝜌𝖩)𝗐+𝜓(𝗏,𝗎)+𝑉𝗋.(3.2) and the symmetric function: 𝑎𝝓1,𝝓2=12𝗇1𝗉2+𝗆1𝗐2+𝗉1𝗇2+𝗐1𝗆2(3.3) for arbitrary states 𝝓1, 𝝓2. The derivative of the energy density with respect to the state vector 𝝓 is given by 𝜕𝝓𝜀=𝖿,𝟢,𝗇,𝗆,(𝜌𝐴)𝗉,(𝜌𝖩)𝗐,(3.4) that leads to the properties: 𝜕𝝓𝜀𝐡=0,12𝜕𝝓𝜀𝝓1𝐟𝝓2=𝑎𝝓1,𝝓2,𝜕𝝓𝜀𝜕𝑠𝐟=𝜕𝑠𝑎(𝝓,𝝓).(3.5) We conclude the local energy balance: 𝜕𝑡𝜀=𝜕𝑠𝑎(𝝓,𝝓),(3.6) that is, 𝑎(𝝓,𝝓) is the energy flux and there is no energy source term. For the presented fixed and free end boundary conditions (2.14), (2.15) we have a vanishing energy flux at the boundaries. Therefore, the total energy is a constant of motion: dd𝑡=dd𝑡𝑠𝑏𝑠𝑎𝜀𝑑𝑠=0.(3.7)

4. Discretization

For the spatial discretization we use a simple finite difference scheme. We note that similar finite difference schemes have been developed and their properties, in particular, the conservation of invariants, have been investigated in [8, 12, 18].

We discretize (𝑠𝑎,𝑠𝑏) with (𝑁+1) equidistant mesh points 𝑠𝑗 and denote the corresponding length of the cells by Δ𝑠. The boundary points are 𝑠1=𝑠𝑎 and 𝑠𝑁+1=𝑠𝑏. As usual, the numerical fluxes are denoted by 𝐅𝑗+1/2 for all 𝑗{1,,𝑁}. The state vector and the source terms are discretely given at the mesh points 𝑗{1,,𝑁+1}, that is, 𝝓𝑗𝝓(𝑠𝑗,𝑡) and 𝐇𝑗𝐡(𝝓(𝑠𝑗,𝑡)). For the inner points 𝑗{2,,𝑁} the spatial discretization results in dd𝑡𝝓𝑗=1Δ𝑠𝐅𝑗+1/2𝐅𝑗1/2+𝐇𝑗,𝐅𝑗+1/212𝐟𝝓𝑗+𝐟𝝓𝑗+1.(4.1) Thereby, the numerical flux function is approximated by the arithmetic average of the flux at neighboring mesh points. For the flux calculation at the boundaries we have to fulfill the Dirichlet boundary conditions at the fixed end 𝑠𝑎=𝑠1: 𝗉1=0,𝗐1=0,(4.2) and at the free end 𝑠𝑏=𝑠𝑁+1: 𝗇𝑁+1=0,𝗆𝑁+1=0.(4.3) Thus, for the remaining components our finite difference scheme reads at the fixed end 𝑠𝑎=𝑠1: dd𝑡𝗋1=0,dd𝑡𝗊1=0,dd𝑡𝗏1=1Δ𝑠𝗉2,dd𝑡𝗎1=1Δ𝑠𝗐2,(4.4) and at the free end 𝑠𝑏=𝑠𝑁+1: dd𝑡𝗋𝑁+1=𝖣1𝗊𝑁+1𝗉𝑁+1,dd𝑡𝗊𝑁+1=𝛀𝗐𝑁+1𝗊𝑁+1,dd𝑡𝗉𝑁+1=1Δ𝑠1𝜌𝐴𝗇𝑁+𝗉𝑁+1×𝗐𝑁+1+1(𝜌𝐴)𝖣𝗊𝑁+1𝖿𝑁+1,dd𝑡𝗐𝑁+1=1Δ𝑠(𝜌𝖩)1𝗆𝑁+(𝜌𝖩)1(𝜌𝖩)𝗐𝑁+1×𝗐𝑁+1.(4.5)

Now, we come up with our main point, the (semi)discrete energy conservation of this scheme. The energy density is approximated locally at the mesh points 𝑗{1,,𝑁+1}, that is, 𝜀𝑗𝜀(𝝓𝑗). We obtain for the inner points 𝑗{2,,𝑁}: dd𝑡𝜀𝑗=𝜕𝝓𝜀𝝓𝑗dd𝑡𝝓𝑗(19)=𝜕𝝓𝜀𝝓𝑗1Δ𝑠𝐅𝑗+1/2𝐅𝑗1/2+𝐇𝑗(17)=1Δ𝑠𝑎𝝓𝑗,𝝓𝑗+1𝑎𝝓𝑗1,𝝓𝑗.(4.6) For the boundaries we have to take into account the Dirichlet conditions: dd𝑡𝜀1=𝜕𝝓𝜀𝝓1dd𝑡𝝓1(16,20,22)=1Δ𝑠𝗇1𝗉2+𝗆1𝗐2(15,20)=2Δ𝑠𝑎𝝓1,𝝓2,dd𝑡𝜀𝑁+1=𝜕𝝓𝜀𝝓𝑁+1dd𝑡𝝓𝑁+1(16,20,23)=1Δ𝑠𝗇𝑁𝗉𝑁+1+𝗆𝑁𝗐𝑁+1(15,21)=2Δ𝑠𝑎𝝓𝑁,𝝓𝑁+1.(4.7) Applying the trapezoidal quadrature rule for the discrete energy, Disc=Δ𝑠2𝜀1+Δ𝑠𝑁𝑗=2𝜀𝑗+Δ𝑠2𝜀𝑁+1,(4.8) guarantees its conservation: dd𝑡disc=Δ𝑠2dd𝑡𝜀1+Δ𝑠𝑁𝑗=2dd𝑡𝜀𝑗+Δ𝑠2dd𝑡𝜀𝑁+1=0.(4.9) This means, the chosen semidiscretization in space ensures that the discrete energy (4.8) is a first integral of the ODE-system (4.1)–(4.5).

For the time discretization any energy conserving method can be used. We choose a Gauss method, that also guarantees the preservation of the norm of the quaternions. In the numerical realization, we make use of the second order Gauss method, that is, the midpoint rule, to obtain a temporal order that is consistent with the spatial one, at least at the inner points. For the discretization of space and time, and the use of the midpoint rule for the conservation of certain properties we refer also to the above mentioned papers [12, 18].

To solve the resulting nonlinear equations a Newton method is used. The strict conservation of energy and orthogonality are the main advantages of the straightforward finite difference scheme presented here.

Remark 4.1. The scheme described above concentrates on preserving the energy of the rod and the orthonormality of the directors. In the sense of numerical methods for hyperbolic systems the scheme is not able to handle shocks properly. It does not have the usual properties like being a TVD scheme or satisfying the entropy condition.

Remark 4.2. Higher order discretizations are also possible. For example, we could consider the following fourth order numerical flux function: 𝐅𝑗+1/212𝐟𝝓𝑗+𝐟𝝓𝑗+1+116𝐟𝝓𝑗+2+𝐟𝝓𝑗+1+𝐟𝝓𝑗𝐟𝝓𝑗1.(4.10) Then, the time derivative of the discrete energy density at the inner mesh point 𝑗 is given by dd𝑡𝜀𝑗=1Δ𝑠𝑎𝝓𝑗,𝝓𝑗+1=𝑂1𝑎𝝓𝑗1,𝝓𝑗=𝑈118𝑎𝝓𝑗,𝝓𝑗+2=𝑂3+14𝑎𝝓𝑗,𝝓𝑗+1=𝑂214𝑎𝝓𝑗1,𝝓𝑗=𝑈2+18𝑎𝝓𝑗2,𝝓𝑗=𝑈3.(4.11) The terms 𝑂1, 𝑂2, 𝑂3, 𝑈1, 𝑈2, and 𝑈3 are eliminated at the mesh points 𝑗+1, 𝑗+1, 𝑗+2, 𝑗1, 𝑗1, and 𝑗2, respectively. Neglecting the boundary points this yields again the conservation of the discrete energy. The discretization near the boundary points has to be considered separately.

5. Numerical Examples

In this section we present three numerical examples, restricting ourselves to Timoshenko’s material law for a homogeneous cylinder as discussed at the end of Section 2. Introducing a typical length, a typical time, and a typical mass: 𝑠typ=𝑠𝑏𝑠𝑎,𝑡typ=𝑠typ𝐸/𝜌,𝑚typ=𝜌𝜋4𝑑2𝑠typ,(5.1) the dimensionless parameters of the model are the slenderness ratio and the gravity number: 𝛿=𝑑𝑠typ,𝛾=𝜌𝑔𝑠typ𝐸.(5.2) In more details, we have (𝐸𝐴)=(𝜌𝐴)=1, (𝐸𝐼)=(𝜌𝐼)=𝜋/16𝛿2, and for the speed of sound 𝑐=1. In the dimensionless form Timoshenko’s material law reads 𝗇=13𝑣1,13𝑣2,𝑣31,𝗆=𝜋16𝛿2𝑢1,𝑢2,23𝑢3.(5.3) To simplify the formulation of the equations we define for 𝗑3 and 𝜆: [𝗑]𝜆=𝑥1,𝑥2,𝜆𝑥3,̂𝗑=𝑥2,𝑥1,0.(5.4) Then, the system (1.1) reads 𝜕𝑡𝗋=𝖣1(𝗊)𝗉,𝜕𝑡𝗊=𝛀(𝗐)𝗊,𝜕𝑡𝗏=𝜕𝑠𝗉+𝗎×𝗉+𝗏×𝗐,𝜕𝑡𝗎=𝜕𝑠𝗐+𝗎×𝗐,𝜕𝑡𝗉=13𝜕𝑠[𝗏]3+13𝗎×[𝗏]3+̂𝗎+𝗉×𝗐𝛾𝖣(𝗊)(0,1,0),𝜕𝑡𝗐=𝜕𝑠[𝗎]1/3+13𝑢3̂𝗎+16𝜋1𝛿2123𝑣3̂𝗏+𝑤3𝗐.(5.5) As mentioned, the rod is fixed at one side, 𝗉(0,𝑡)=(0,0,0),𝗐(0,𝑡)=(0,0,0),(5.6) and the other side is free, that is, for Timoshenkos’s material law: 𝗏(1,𝑡)=(0,0,1),𝗎(1,𝑡)=(0,0,0).(5.7)

The chosen initial configuration of a straight rod and direction of gravity 𝖾𝑔=(0,1,0) can be seen in Figure 1. In the following examples different initial torsions will be considered. More precisely, 𝐫(𝑠,0)=𝑠𝐝3,𝐝1(𝑠,0)=𝐞2cos(𝜇(𝑠))+𝐞3sin(𝜇(𝑠)),𝐝2(𝑠,0)=𝐞2sin(𝜇(𝑠))+𝐞3cos(𝜇(𝑠)),(5.8) where 𝜇(0,1) is a field of torsion angles that has to be defined. These conditions are equivalent to initial conditions for 𝗋(𝑠,0) and 𝗊(𝑠,0). Moreover, due to the definitions of 𝐯 and 𝐮 we have 𝗏(𝑠,0)=(0,0,1),𝗎(𝑠,0)=0,0,𝜕𝑠𝜇.(5.9) Finally we prescribe 𝗉(𝑠,0)=(0,0,0),𝗐(𝑠,0)=(0,0,0),(5.10) that is, the rod is initially at rest. Our initial conditions are compatible to the boundary conditions if 𝜇(0)=0 and 𝜕𝑠𝜇(1)=0.

In all simulations we choose the CFL-number equal to 1. This implies Δ𝑡=Δ𝑠 as for the speed of sound 𝑐=1 in the dimensionless form.

We remark that in all simulations energy is strictly conserved according to the above analysis.

Example 5.1 (Torsional Oscillation without Gravity). We choose 𝛾=0,𝛿=103, and a field of torsion angles fulfilling the compatibility condition: 𝜇(𝑠)=2𝜋sin𝜋2𝑠.(5.11) In this example, only the torsion 𝑢3 and the angular velocity 𝑤3 are involved, because of the vanishing gravity. The equations reduce to the homogenous wave-equation: 𝜕𝑡𝑢3𝜕𝑠𝑤3=0,𝜕𝑡𝑤313𝜕𝑠𝑢3=0.(5.12) Due to the chosen torsion angle we exactly initialize the fundamental mode of the wave equation with a frequency 𝜔theo=𝜋12.(5.13) We use this example as a benchmark for the convergence properties of our scheme comparing the computed frequencies with the analytical one for different grid sizes, see Figure 2. Comparing the identified simulation frequency with the fundamental frequency of the analytical solution 𝜔theo=0.9069 we note that, for example, for 𝑁=25 the relative error is of the order 103.

Example 5.2 (Transversal Oscillation). We choose 𝛾=107, 𝜇=0, and vary 𝛿 from 103 to 101. This example describes the classical case of transversal oscillation of an one-sided fixed rod without twist. Since the rod is strainless in the initial configuration and—as mentioned—at rest initially, gravity is used to excite the oscillation. The transversal oscillation frequency for the rod can be compared to the well-known result for vanishing gravitation given for example in Peterson [19] or Timoshenko [20]. It reads in the chosen dimensionless form 𝜔=𝜆24𝛿,(5.14) where the constant 𝜆 is given by 𝜆=1,875. Defining the temporal difference between two maxima of the potential energy as a period of oscillation, we identify the frequency 𝜔sim. Figure 3 shows for 𝑁=10 the expected linearity of the frequency 𝜔 with respect to the slenderness ratio 𝛿. Determining 𝜆 from a best fit line, we obtain 𝜆sim=1.870. Thus, the relative error is very small1—of the order 103—even for the very coarse discretization with 𝑁=10.
We note that the solution is very accurate although we have chosen a linear material law for the contact force 𝗇 instead of a Kirchhoff constraint 𝗏=(0,0,1) normally used.

Example 5.3 (Three-dimensional Problem Including all Strain Variables). We choose 𝛾=104, 𝛿=102, and a field of torsion angles: 𝜇=6𝜋𝑠(2𝑠).(5.15) In this case, see Figure 4, the rod is initially twisted three times and the gravitational number is comparatively high. That means all strain variables of the rod are excited during the evolution. The algorithm is able to deal with such a situation, in particular, from the point of view of conservation of energy. The results are reasonable as long as the linear material law is valid. One has to take into account that for large strain values Timoshenko’s material law does not guarantee that the deformation of the rod preserves the orientation. See Antman [10] for a precise definition of the preservation of orientation.

6. Conclusion

In this paper we use the description of a hyperelastic rod in the formulation of Antman [10]. For the resulting hyperbolic system we developed a numerical method, which conserves the energy of the rod as well as the orthonormality of the directors. However, the scheme is not able to handle shocks properly. It is neither a TVD scheme nor does it satisfy the entropy condition.

For the material law of Timoshenko [17] we illustrated the method using some numerical examples. For these examples, the conservation properties are strictly fulfilled and very good agreement of the numerical and analytical results can be observed. Finally, we mention that for realistic three-dimensional problems a nonlinear material law has to be used, which can be easily incorporated in the scheme.


This work has been supported by Deutsche Forschungsgemeinschaft (DFG), KL 1105/18-1 and WE 2003/3-1 and by Rheinland-Pfalz Excellence Center for Mathematical and Computational Modeling (CM)2.


  1. D. K. Pai, “Interactive simulation of thin solids using Cosserat models,” Computer Graphics Forum, vol. 21, no. 3, pp. 347–352, 2002. View at: Google Scholar
  2. J. Barbic and D. James, “Real-time subspace integration for St. Venant-Kirchhoff deformable models,” ACM Transactions on Graphics, vol. 24, no. 3, pp. 982–990, 2005. View at: Google Scholar
  3. A. Weber, T. Lay, and G. Sobottka, “Stable integration of the dynamic Cosserat equations with application to hair modeling,” Journal of WSCG, vol. 16, pp. 73–80, 2008. View at: Google Scholar
  4. J. C. Simo, “A finite beam formulation. The three dimensional dynamic problem. Part I,” Computer Methods in Applied Mechanics and Engineering, vol. 49, pp. 55–70, 1985. View at: Google Scholar
  5. J. C. Simo and L. Vu-Quoc, “A three dimensional finite-strain rod model. Part II: computational aspects,” Computer Methods in Applied Mechanics and Engineering, vol. 58, pp. 79–116, 1986. View at: Google Scholar
  6. J. C. Simo and L. Vu-Quoc, “On the dynamics in space of rods undergoing large motions—a geometrically exact approach,” Computer Methods in Applied Mechanics and Engineering, vol. 66, no. 2, pp. 125–161, 1988. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  7. J. C. Simo, J. E. Marsden, and P. S. Krishnaprasad, “The Hamiltonian structure of nonlinear elasticity: the material and convective representations of solids, rods, and plates,” Archive for Rational Mechanics and Analysis, vol. 104, no. 2, pp. 125–183, 1988. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  8. L. Vu-Quoc and X. Tan, “Optimal solid shells for nonlinear analyses of multilayer composites. part ii: dynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 192, no. 9-10, pp. 1017–1059, 2003. View at: Google Scholar
  9. E. Wittbrodt, I. Adamiec-Wojcik, and S. Wojciech, Dynamic of Flexible Multibody Systems, Springer, 2006.
  10. S. S. Antman, Nonlinear Problems of Elasticity, Springer, New York, NY, USA, 2005. View at: Zentralblatt MATH
  11. J. C. Simo, N. Tarnow, and K. K. Wong, “Exact energymomentum conserving algorithms and symplectic schemes for nonlinear dynamics,” Computer Methods in Applied Mechanics and Engineering, vol. 100, no. 1, pp. 63–116, 1992. View at: Google Scholar
  12. S. Li and L. Vu-Quoc, “Finite difference calculus invariant structure of a class of algorithms for the nonlinear Klein-Gordon equation,” SIAM Journal on Numerical Analysis, vol. 32, no. 6, pp. 1839–1875, 1995. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  13. L. Vu-Quoc and J. Simo, “On the dynamics of earth-orbiting flexible satellites with multibody components,” AIAA Journal of Guidance, Control, and Dynamics, vol. 10, no. 6, pp. 549–558, 1987. View at: Google Scholar
  14. L. Vu-Quoc, Dynamics of flexible structures performing large overall motions: a geometrically-nonlinear approach, technical report no. UCB/ERL M86/36 [Ph.D. thesis], University of California at Berkeley, 1986.
  15. F. A. McRobie and J. Lasenby, “Simo-Vu Quoc rods using Clifford algebra,” International Journal for Numerical Methods in Engineering, vol. 45, no. 4, pp. 377–398, 1999. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  16. L. Mahadevan and J. B. Keller, “Coiling of flexible ropes,” Proceedings of the Royal Society A, vol. 452, no. 1950, pp. 1679–1694, 1996. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  17. S. P. Timoshenko, “On the correction for shear of differential equations for transverse vibrations of prismatic bars,” Philosophical Magazine, vol. 6, no. 41, pp. 744–746, 1921. View at: Google Scholar
  18. L. Vu-Quoc and S. Li, “Invariant-conserving finite difference algorithms for the nonlinear Klein-Gordon equation,” Computer Methods in Applied Mechanics and Engineering, vol. 107, no. 3, pp. 341–391, 1993. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  19. C. Peterson, Dynamik der Baukonstruktionen, Vieweg, 1996.
  20. S. P. Timoshenko, Schwingungsproblem der Technik, Springer, 1932.

Copyright © 2012 Thorsten Fütterer et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

 PDF Download Citation Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

We are committed to sharing findings related to COVID-19 as quickly as possible. We will be providing unlimited waivers of publication charges for accepted research articles as well as case reports and case series related to COVID-19. Review articles are excluded from this waiver policy. Sign up here as a reviewer to help fast-track new submissions.