Geometrically exact physics-based modeling and computer animation of highly flexible 1D mechanical systems

Graphical Models 75 (2013) 56–68

Contents lists available at SciVerse ScienceDirect

Graphical Models journal homepage: www.elsevier.com/locate/gmod

Geometrically exact physics-based modeling and computer animation of highly flexible 1D mechanical systems q Ye Duan a,⇑, Dong Li a, P. Frank Pai b a b

Department of Computer Science, University of Missouri, Columbia, MO 65211, USA Department of Mechanical and Aerospace Engineering, University of Missouri, Columbia, MO 65211, USA

a r t i c l e

i n f o

Article history: Received 18 April 2012 Received in revised form 22 December 2012 Accepted 1 January 2013 Available online 4 February 2013 Keywords: Physics-based animation Geometrically exact beam theory Nonlinear finite element Flexible multibody systems

a b s t r a c t This paper presents a geometrically exact beam theory and a corresponding displacementbased finite-element model for modeling, analysis and natural-looking animation of highly flexible beam components of multibody systems undergoing huge static/dynamic rigidelastic deformations. The beam theory fully accounts for geometric nonlinearities and initial curvatures by using Jaumann strains, concepts of local displacements and orthogonal virtual rotations, and three Euler angles to exactly describe the coordinate transformation between the undeformed and deformed configurations. To demonstrate the accuracy and capability of this nonlinear beam element, nonlinear static and dynamic analysis of two highly flexible beams are performed, including the twisting a circular ring into three small rings and the spinup of a flexible helicopter rotor blade (Graphical abstract). These numerical results reveal that the proposed nonlinear beam element is accurate and versatile for modeling, analysis and 3D rendering and animation of multibody systems with highly flexible beam components. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Thin and flexible one-dimensional (1D) mechanical components are important parts for computer animation of dynamical systems related to human locomotion [1– 3]. Examples include human hairs, jump ropes, rope swing, rope dance, cowboy ropes, telephone cables, grape vines, archery bows, clad poles, sabre fencing, violin strings, fishing rods, and surgical sutures. Moreover, highly flexible beams are used in many mechanical, civil, aerospace, and architectural systems [4–7], such as helicopter rotor blades, wings of high-altitude long-endurance aircrafts, aviation propeller blades, wind-turbine blades, robot manipulators, slender space structures for buildings, armtype positioning mechanisms of magnetic disk drives, and flexible links for high-speed slider-crank mechanisms. q

This paper has been recommended for acceptance by Peter Lindstrom.

⇑ Corresponding author.

E-mail address: [email protected] (Y. Duan). 1524-0703/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.gmod.2013.01.001

Moreover, NASA’s various science missions have extensively used deployable/inflatable structures consisting of highly flexible beams in order to reduce the stowed volume and weight during launch, minimize extra vehicular activities in space, and/or reduce the operation time and cost [6–8]. Furthermore, although cables are beams with very small bending rigidity, study of loops and kinks of cables requires the use of a nonlinear beam model that can account for dramatic geometric nonlinearities. Accurate description of static/dynamic deformations of such systems is important for structural engineering as well as computer animation. However, accurate, natural-looking, and real-time 3D rendering of such highly nonlinear systems requires a geometrically exact beam theory with a high-fidelity discretization and solution method. Today’s surgical and medical treatments for many diseases depend on the use of highly flexible beam- or cable-like medical devices and tools. For example, natural orifice translumenal endoscopic surgery (NOTES) has recently emerged as a favorable surgical technology under

Y. Duan et al. / Graphical Models 75 (2013) 56–68

worldwide research and development [9,10]. NOTES uses flexible endoscopic devices to access body cavities through the mouth, vagina, or anus without skin incisions and perform common surgical procedures and retraction. Such surgical tools are required to be small in diameter and highly flexible in order to follow and pass through the gullet, intestine, natural orifices and/or endoscope channels easily and safely but do not undergo kinking, permanent deformation, or breaking during surgery. Also, accurate and smooth handling of sutures is important for a surgeon to have a successful surgery when s/he performs delicate operations while looking at a computer monitor. Unfortunately, a suture is a highly flexible beam that exhibits a variety of complex mechanical behaviors under a confined space. Hence, high-fidelity computer animation of static/ dynamic characteristics of sutures is an important part of any surgical simulation system. Even in molecular biology, an accurate nonlinear beam theory is needed [11–13]. The supercoiling (or writhing) of DNA is known to affect every physical, chemical, and biological property of a molecule. Because supercoiled DNA is an important functional state active in the processes of replication, transcription, and recombination, it is important to model and predict the structural properties of DNA in higher-order forms (supercoils, knots, catenanes, and protein-DNA complexes) in order to understand DNA’s fundamental functions, including strand unwinding (replication, transcription) and passage (knotting and catenation), looping, and slithering. For modeling and analysis, a double-helical DNA polymer can be modeled as a very thin beam with initial curvatures, and nonlinear buckling analysis is the main task for understanding DNA’s fundamental functions. This has stimulated many theoretical, computational, and experimental studies of mechanics and stability of thin beams, including the controlled buckling of elastic beams and the ingenious manipulations of DNA strands [13]. Hence, to advance theoretical structural mechanics for today’s science and engineering applications, it is important to derive a geometrically exact beam theory that can be used to investigate influences of highorder geometric nonlinearities on static and dynamic characteristics of highly flexible beams and cables. Beam theories are 1D mathematical models for structures with one dimension being much larger than the other two. Modeling of beams is more difficult than modeling of plates and shells because the former transforms a 3D problem in nature into a 1D problem whereas the latter transforms a 3D problem into a 2D problem. Beam theories can be divided into three groups of different complexity: (a) Euler–Bernoulli beam theory, (b) shear-deformable beam theories (i.e., Timoshenko’s theory, third-order shear theory, higher-order shear theory, layer-wise shear theory, etc.), and (c) 3D beam theories [7]. In the Euler–Bernoulli beam theory, only the axial stress r11 is considered and a plane cross section perpendicular to the reference axis before deformation is assumed to remain plane and perpendicular to the deformed reference axis after deformation. In shear-deformable beam theories, transverse shear stresses r12 and r13 and out-of-plane warpings are taken into account. In 3D beam theories, both out-of-plane and inplane warpings are taken into account. In-plane warpings

57

and hence transverse normal strains e22 and e33 and inplane shear strain e23 may result in significant stresses r22, r33, and r23. Hence, all six stresses are accounted for in a 3D beam theory. However, because in-plane and outof-plane warpings are relative displacements with respect to the deformed cross section and are much smaller than global displacements, inertial forces caused by warpings are negligible. But, because these warpings offer extra degrees of freedom for the cross-section to deform, they significantly affect the global stiffnesses of a beam. Hence, the geometrically nonlinear problem of elastic beams can be decoupled into a nonlinear 1D problem and a linear 2D sectional problem [14,15]. In other words, a 1D nonlinear beam model with global stiffnesses determined from a static 2D sectional analysis of warpings is a general and practical approach in solving nonlinear beam problems [15,16]. Mechanical systems are multibody systems. A multibody system consists of interconnected rigid and deformable components, and each component may undergo large translations and rotations [7,17,18]. Modeling and analysis of a flexible component that undergoes large rotations is very challenging because geometric nonlinearities exist and they couple with the nonlinear ordinary differential equations that govern the motions of rigid components. Hence, finite element techniques are often used in the modeling and analysis of flexible multibody systems. Even with the use of finite elements, however, many challenging problems still exist, and the most serious challenge is how to accurately describe the large rotations of flexible components. One way to release this complexity is to derive and use totalLagrangian structural theories. Because the strain–displacement relations of a total-Lagrangian beam theory account for rigid-body effects, there is no need of separate consideration of the rigid-body movement of the rigid component that the beam is connected to. However, the most serious difficulty is to derive objective strains using the absolute displacements and the not-well-defined torsional variable because they contain rigid-body effects [7]. In a formulation using a floating frame or many element-based reference frames, local displacements are defined with respect to a corotated local coordinate system that is assumed to account for rigid-body displacements due to high flexibility [19,20]. Hence, the local displacements are small and mainly caused by straining. If linear strain–displacement relations are used, the stiffness matrix is constant. Unfortunately, the origin of the local coordinate system needs three absolute displacements to describe its motion with respect to the fixed global coordinate system. These two sets of displacements are coupled through a deformation-dependent coordinate transformation matrix, and this nonlinear coupling and the so-induced gyroscopic effects (including Coriolis, centrifugal, and Euler forces) make the obtained mass matrix nonlinear and displacement-dependent. Moreover, when an element undergoes large relative rotations, nonlinear strain–displacement relations (either a fully nonlinear version or a truncated version) are needed and they cause a complicated displacement-dependent stiffness matrix [21]. If linear strain–displacement relations are used for such cases, the axial displacement cannot be accurately accounted for and artificial dynamic stiffening may exist

58

Y. Duan et al. / Graphical Models 75 (2013) 56–68

[19]. A typical problematic case is the spinup of a gravitydeflected helicopter rotor blade. The spin-induced centrifugal force raises up the blade with a significant axial displacement due to rigid-body rotations, instead of forced stretching. Hence, concepts of substructuring and static/ dynamic condensations, use of non-Cartesian variables and different orders or types of strain–displacement relations, and many other approaches have been proposed to overcome such drawbacks in using floating or elementbased reference frames [22–28]. The Cosserat rod theory is often used in the structural engineering and computer graphics communities for modeling of thin beams undergoing large deformations [2–5]. The curvatures (and hence strains) are spatial derivatives of the three unit vectors (also called directors) fixed on the deformed geometry. However, the challenge is how to relate the curvatures to the three displacement variables and one torsional rotation at each point on the reference line [29,30]. Four Euler parameters or three original or modified Rodrigues parameters are often used to relate the deformed unit vectors to the undeformed ones. Unfortunately, because these parameters are not explicit functions of displacement variables, the so-derived beam theory is not fully displacement-based [29,30]. Hence, stress- and strain-based theories are often adopted in order to use the Cosserat rod theory [1–3,16]. Unfortunately, such beam theories do not have fully nonlinear, explicit strain–displacement relations and cannot provide geometrically exact solutions because the displacements need to be obtained from spatial, approximate step-by-step integration [30]. This approach is essentially an updated-Lagrangian approach and the obtained solution cannot be claimed to be geometrically exact. The accumulation of numerical errors and other problems related to the use of updated Lagrangian formulations also exist in the deformed geometry obtained from this approach. Literature reviews indicate that totalLagrangian formulation is more favorable than updatedLagrangian formulation for highly flexible elastic structures. A geometrically exact total-Lagrangian displacementbased beam theory based on the concept of local relative displacements, Jaumann strains, exact coordinate transformations, and orthogonal virtual rotations exists in the literature [7,31]. Because this beam theory uses only two Euler angles, it is versatile for many applications since its governing equations are independent of the rotation sequence of Euler angles [7]. However, this beam theory’s torsional variable can be discontinuous and hence invalidates any continuity-based methods (including the finite element method) for solving some beam structures undergoing large rigid-elastic rotations, such as the spinup of rotor blades [32]. To overcome this problem, this work derives a geometrically exact displacement-based beam theory using a 3–2–1 rotation sequence for three Euler rotations to describe the transformation between the deformed and the undeformed coordinate systems. This beam theory uses only four dependent variables, including three displacement variables and one torsional angle. Moreover, the corresponding finite element formulation is derived, and the accuracy and capability of this new nonlinear beam theory is verified through two advanced numerical examples.

2. Geometrically exact modeling of initially curved beams We consider a flexible beam undergoing large overall displacements, as shown in Fig. 1a. Three coordinate systems are needed for describing the undeformed and deformed geometries of an initially curved beam undergoing large rigid-body movement and elastic deformation [7]. The abc is a fixed rectangular coordinate system used for reference; the xyz is an orthogonal curvilinear coordinate system whose x axis is the reference line formed by connecting the reference points of all cross sections of the undeformed beam; and the ngf is a local orthogonal curvilinear coordinate system whose n axis represents the deformed reference line and g and f represent the deformed axes y and z only if there is no cross-section warping. Moreover, s denotes the arc length along the reference line starting from the beam root, and u, v, and w represent the absolute displacements of the observed reference point O (see Fig. 1b) with respect to the axes x, y, and z, respectively. We let ia, ib, and ic be the unit vectors of the abc system; ix, iy, and iz be the unit vectors of the xyz system; and i1, i2, and i3 be the unit vectors of the ngf system. Fig. 1b shows the displacement components of an arbitrary point on the observed cross section. Because the undeformed beam geometry is assumed to be known, the position vector R = A(s)ia + B(s)ib + C(s)ic of the undeformed reference point O and the angles h21(s), h22(s), and h23(s) between the y axis and the a, b, and c axes are known. Hence, the transformation matrix [T0] that relates the coordinate systems abc and xyz is obtained as fixyz g ¼ ½T o fiabc g; fixyz g  fix ; iy ; iz gT ; fiabc g  fia ; ib ; ic gT 2 6 ½T  ¼ 4 o

ð1aÞ

3 A0 B0 C0 7 cos h21 cos h22 cos h23 5 0 0 0 0 0 0 B cos h23  C cos h22 C cos h21  A cos h23 A cos h22  B cos h21

ð1bÞ 0

where ( )  d( )/ds. Moreover, using the orthonormality of ix, iy, and iz we obtain

2 d fixyz g ¼ ½Kfixyz g; ds

0

6 ½K  4 k3 k2 o

k1 ¼

3 X diy dT 2i o  iz ¼ T ; ds ds 3i i¼1

k3 ¼

3 X dix dT 1i o  iy ¼ T ds ds 2i i¼1

k2 ¼ 

k3 0 k1

k2

3

7 k1 5 0

ð2aÞ

o 3 X dix dT 1i o  iz ¼  T ; ds ds 3i i¼1

o

ð2bÞ

where k1, k2, and k3 are the initial curvatures with respect to the axes x, y, and z, respectively. Here we use three successive Euler angles h3, h2, and / to describe the rotation of a cross section from its undeformed position to its deformed one, as shown in Fig. 2. In Fig. 1b, if OP ¼ ds, the deformed one is given by P  ¼ ð1 þ eÞdsi1 , where e is the axial strain on the n axis. O Hence, we have

59

Y. Duan et al. / Graphical Models 75 (2013) 56–68

(a)

(b)

Fig. 1. An elastic beam undergoing large overall motion: (a) three coordinate systems and (b) different displacements.

i1  T 11 ix þ T 12 iy þ T 13 iz ¼

P  O ð1 þ eÞds 0

0

0

ix þ u0 ix þ v 0 iy þ w0 iz þ uix þ v iy þ wiz 1þe 1 þ u0  v k3 þ wk2 v 0 þ uk3  wk1 w0  uk2 þ v k1 T 11 ¼ ; T 12 ¼ ; T 13 ¼ 1þe 1þe 1þe qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 0 0 0 e ¼ ð1 þ u  v k3 þ wk2 Þ þ ðv þ uk3  wk1 Þ þ ðw  uk2 þ v k1 Þ2  1 ¼

ð3Þ

where the expression of e is derived from the fact qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ji1 j ¼ T 211 þ T 212 þ T 213 ¼ 1. Eq. (3) shows that the n axis is fully related to the xyz system by displacements u, v, and w. Moreover, the coordinate systems ng1 and xyz are related 8 9 as 8 9

2 3 T 11 T 12 T 13 > > < i1 > < ix > = = 6 7 i2 ¼ ½T iy ; ½T  4 T 21 T 22 T 23 5 > > : > : > ; ; T 31 T 32 T 33 i3 iz 2 32 3 T 11 T 12 T 13 1 0 0 6 76 7 ¼ 4 0 cos / sin / 54 T 21 T 22 0 5 T 31

0  sin / cos /

T 32

Eq. (4) reveals that [T] is a function of T11, T12, T13, and /. However, because T 211 þ T 212 þ T 213 ¼ 1, [T] is a function of only three independent variables. Eqs. (4) and (3) also show that the large rotation between the coordinate systems xyz and ng 1 can be exactly described by u, v, w, u0 , v0 , w0 , and / (see Fig. 1a). Taking derivatives of Eq. (4) with respect to s, using Eq. (2a) and the identity [T]T[T] = [I], we obtain

@ fi123 g ¼ ½Kfi123 g @s 2 3 0 q3 q2 6 0 T T 0 q1 7 ½K  4 q3 5 ¼ ½T ½T þ ½T½K½T q2 q1 0

ð4Þ

q1  i02  i3 ¼

T 33

3 X  0  T 2i T 3i þ T 1i ki i¼1

T 11 T 11 T 13 ffi ; T 31 ¼ p ffiffiffiffiffiffiffiffiffi ffi; T 22 ¼ pffiffiffiffiffiffiffiffiffi 1T 213 1T 213 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T 13 ffi ; T 33 ¼ 1  T 213 ¼  pT 12ffiffiffiffiffiffiffiffiffi 2

¼ /0 þ

1T 13

1T 13

q2

T 11 T 012 T 13  T 011 T 12 T 13

þ T 11 k1 þ T 12 k2 þ T 13 k3 1  T 213 3 X  0  1 0 ½T 31 ðu0  v k3 þ wk2 Þ0  i1  i3 ¼ T 1i T 3i þ T 2i ki ¼ 1þe i¼1

η

z

ð5bÞ

where q1 is the deformed twisting curvature, and q2 and q3 are the deformed bending curvatures with respect to the axes g and f. Moreover, it follows from Eqs. (5b), (4) and (2a) that

12 ffi ffiffiffiffiffiffiffiffiffi T 21 ¼ pT ; 2

T 32

ð5aÞ

þ T 32 ðv 0 þ uk3  wk1 Þ0 þ T 33 ðw0  uk2 þ v k1 Þ0  þ

3 X T 2i ki

ð6Þ

i¼1

ζ

y

θ3

z

φ

y

ξ

q3  i01  i2 ¼

i¼1

-θ2

φ

-θ2

φ

3 X  0  T 1i T 2i þ T 3i ki

¼

θ3

T13

i1

-θ2

i¼1

x

T12

θ3

x

T11 Fig. 2. Kinematic relations between three Euler angles.

1 ½T 21 ðu0  v k3 þ wk2 Þ0 þ T 22 ðv 0 þ uk3  wk1 Þ0 1þe 3 X þ T 23 ðw0  uk2 þ v k1 Þ0  þ T 3i ki

Note that e, q1, q2, and q3 are global strains and they are invariant under rigid-body motions. It follows from Eq. (3) that the variation of the axial strain is

de ¼ T 11 dð1 þ u0  v k3 þ wk2 Þ þ T 12 dðv 0 þ uk3  wk1 Þ þ T 13 dðw0  uk2 þ v k1 Þ

ð7Þ

60

Y. Duan et al. / Graphical Models 75 (2013) 56–68

Variations of the unit vectors i1, i2, and i3 are due to virtual rotations of the axes n, g, and f and are given by

8 9 2 0 > < di1 > = 6 di2 ¼ 4 dh3 > > : ; dh2 di3

dh3 0 dh1

38 9 dh2 > < i1 > = 7 dh1 5 i2 > : > ; 0 i3

ð8Þ

where dh1, dh2, and dh3 are virtual rotations with respect to axes n, g, and f, respectively. By using Eqs. (8), (5a), (4) and (3), we obtain

dh1 ¼ di2  i3

characteristics of geometrically exact beam theories without those complex mathematics and notations related to shear-deformable and 3D beam theories [7]. For fully nonlinear strain–displacement relations and inertial terms that include the influences of U, the derivation steps are essentially the same and the reader is referred to [7] for details. 3. Finite element formulation The extended Hamilton principle is used to formulate the problem. The variation of elastic energy is given by

T 13  ½T 11 k3 du þ T 12 k3 dv ¼ d/ þ  1  T 213 ð1 þ eÞ  ðT 11 k1 þ T 12 k2 Þdw  T 12 du0 þ T 11 dv 0 dh2

dP ¼

2

T 33 k2  T 32 k3 T 31 k3  T 33 k1 du þ dv 1þe 1þe T 32 k1  T 31 k2 T 31 T 32 þ dw  du0  dv 0 1þe 1þe 1þe T 33  dw0 dh3 1þe ¼ di1  i2 ¼

T 22 k3  T 23 k2 T 23 k1  T 21 k3 du þ dv 1þe 1þe T 21 k2  T 22 k1 T 21 T 22 þ dw þ du0 þ dv 0 1þe 1þe 1þe T 23 þ dw0 1þe

ð9Þ

dq1 ¼ ðdh1 Þ0  q3 dh2 þ q2 dh3 ð10Þ

0

dq3 ¼ ðdh3 Þ  q2 dh1 þ q1 dh2 Eqs. (10) and (9) can be used to obtain variations of curvatures in terms of variations of u, v, w, /, and their spatial derivatives. In order to fully account for geometric nonlinearities, we use Jaumann strains because they are fully nonlinear, objective, and geometric strain measures and their directions are defined with respect to the deformed coordinate system, which implies co-rotation for every material point [7]. If the local displacement vector U (see Fig. 1b) caused by out-of-plane and in-plane warpings are neglect (i.e. the Euler–Bernoulli beam theory), Jaumann strains Bij can be derived to be

8 9 > < B11 > = feg ¼ ½Sfwg; feg ¼ 2B12 ; > > : ; 2B13 8 9 2 3 e > > > > 1 0 z y > = 1 6 7 1 ; ½S  4 0 z 0 0 5 fwg ¼ > q  k 2> > > 2 > > 0 y 0 0 : ; q3  k3

Z

fdegT frg dA ds;

frg 

A

E

0

0

0

0

G

b ¼ 6 ½Q 40 G

Taking variations of Eq. (6) and using the fact  0 d ik ¼ ðdik Þ0 and Eqs. (8), (5a) and (5b), we obtain

dq2 ¼ ðdh2 Þ0 þ q3 dh1  q1 dh3

L

0

¼ di1  i3

¼

Z

3

8 9 > < r11 > = > :

r12 ¼ ½ Qb feg; > r13 ;

7 05

ð12Þ

where A denotes the cross-sectional area, L is the beam b  is the length, {r} is the Jaumann stress vector, and ½ Q reduced material stiffness matrix. For an isotropic beam, b  is given in Eq. (12), where E is Young’s modulus its ½ Q and G is the shear modulus. For a laminated or composb  can be a full 3  3 matrix [7]. Note that ite beam, its ½ Q b  can be determined by experiments using engithe ½ Q neering strain and stress measures because Jaumann strains are nothing but co-rotated engineering strains without influences of rigid-body movement [7]. On the b other hand, if Green–Lagrange strains are used, the ½ Q needs to be determined by experiments using second Piola–Kirchhoff stresses and Green–Lagrange strains, which is not usually done in experiments and is nonlinear and hence computationally awkward. By using Eq. (12), we obtain

dP ¼

Z

L

0

Z

b feg dA ds ¼ fdegT ½ Q

Z

A

L

fdwgT ½Dfwg ds

ð13aÞ

0

where

2 ½D 

Z

EA

6 0 b ½S dA ¼ 6 ½ST ½ Q 6 4 0 A

fI1 ; I2 ; I3 g 

Z

0

3

0

0

0

GI1

0

0

EI2

0 7 7 7 0 5

0

0  2  ðy þ z2 Þc1 ; z2 ; y2 dA

EI3

A

fdwg ¼ ½WfdUg;

U ¼ fu; u0 ; u00 ; v ; v 0 ; v 00 ; w; w0 ; w00 ; /; /0 gT ð13bÞ

ð11Þ

Here we use the simple Euler–Bernoulli beam theory to clearly demonstrate the derivation process and the main

Here, c1 is a correction factor accounting for the decrease of torsional rigidity due to out-of-plane torsional warping and can be calculated using the theory of elasticity [33]. For example, c1 = 1 for a circular cross section, and c1 < 1 for non-circular cross sections. The entries Wij (@ wi/ @Uj) of [W] can be derived from Eqs. (7), (10) and (9). Next we use two-node beam elements to discretize the beam into ne elements. The absolute displacements u, v, w, and / of the jth element are discretized as

61

Y. Duan et al. / Graphical Models 75 (2013) 56–68

fdg  fu; v ; w; /gT ¼ ½NfqðjÞ g

ð14aÞ

dK e ¼ 

where [N] is a 4  14 matrix of shape functions and {q(j)} is the element displacement vector of the jth element defined as

Substituting Eq. (14a) into Eq. (13b), we get

fUg ¼ ½@NfqðjÞ g; ½@N  ½@½N

dP ¼

ne Z X LðjÞ

j¼1

fdqðjÞ gT ½@NT ½WT ½Dfwg ds

ne X b fqg ¼ fdqðjÞ gT ½K ðjÞ fqðjÞ g ¼ fdqgT ½ K

ð16aÞ

j¼1

ðjÞ

ðjÞ

j¼1

½K fq g ¼

LðjÞ

¼

Z

Z 0

L 0

Z

½@N ½W ½Dfwg ds

½M ðjÞ  

A

Z

0

0

ð18bÞ

0

L

ðq1 du þ q2 dv þ q3 dw þ q4 dh1 þ q5 dh2

0

j¼1 T

¼ fdqg fRg

ð19aÞ

where q1, q2, and q3 are distributed forces along the axes x, y, and z, respectively, and q4, q5, q6 are distributed twisting and bending moments along the axes n, g, and f, respectively. Moreover, {R} is the global nodal loading vector and {R(j)} is the elemental nodal loading vector given by

Z LðjÞ

e ½@NT f Rgds

ð19bÞ

The extended Hamilton principle states as

Z

t2

ðdK e  dP þ dW nc Þdt ¼ 0

ð20Þ

t1

Ah2  j2 €h2  ðj3  j1 Þh_ 3 h_ 1 ;

Substituting Eqs. (16a), (18a) and (19a) into Eq. (20) yields the following equation of motion:

Ah3  j3 €h3  ðj1  j2 Þh_ 1 h_ 2 Z Z Z m qdA; j1  qðy2 þ z2 ÞdA; j2  qz2 dA; j3 

~ ½NT ½m½N ds;

þ q6 dh3 Þ ds Z L ne X e ds ¼ fdUgT f Rg fdqðjÞ gT fRðjÞ g ¼

€ du þ mv€ dv þ mwdw € ðmu þ Ah1 dh1

A

LðjÞ

For a discrete rigid body, its inertial moments Ahi are the same as those in Eq. (17) except that the rotary inertias jk need to be obtained by integration with respect to its volume, instead of the cross-section area. The variation of non-conservative work due to six general external distributed loads is given by

L

A

Z

0

fRðjÞ g 

þ Ah2 dh2 þ Ah3 dh3 Þ ds

Z

3 m 0 0 0 6 0 m 0 07 7 ~ 6 ½m 7 6 4 0 0 m 05 2

ð16bÞ

€  dD dA ds qD

Ah1  j1 €h1  ðj2  j3 Þh_ 2 h_ 3 ;

ð18aÞ

0

T

ne is the total number of elements, L(j) is the length of jth b  is the global stiffness matrix, and {q} is the element, ½ K global displacement vector. We note that element stiffness matrix [K(j)] and element displacement vector {q(j)} are nonlinearly coupled into a vector and cannot be separated in this fully nonlinear formulation. If all cross-section warpings are neglected, the variation of kinetic energy dKe is given by [7]

dK e ¼ 

€ðjÞ g ds ~ fdqðjÞ gT ½NT ½m½Nf q

where [M] is the global mass matrix assembled from the elemental mass matrix [M(j)], which is defined as

dW nc ¼ T

LðjÞ

€g ¼ fdqgT ½Mfq

where

Z

€ du þ mv€ dv þ mwdwÞ € ðmu ds

¼

ð15Þ

where [@] is a 11  4 matrix of differential operators. Then we obtain the variation of elastic energy dP as

L

0 ne Z X

 T fqðjÞ g ¼ ui ; v i ; wi ; /i ; w0i ; v 0i ; u0i ; uk ; v k ; wk ; /k ; w0k ; v 0k ; u0k ð14bÞ

Z

b fqg ¼ fRg €g þ ½Cfqg _ þ ½K ½Mfq

A

2

qy dA

A

ð17Þ where ü  @ 2u/@t2, t is the time, q is the mass density, jk are rotary inertias per unit length of the beam, and D is the absolute displacement vector of an arbitrary point on the observed cross section, as shown in Fig. 1b. Eqs. (17) and (9) show that the inertial moments Ahi are nonlinear functions of u, v, w, and / and their spatial derivatives. However, because rotary inertias jk of highly flexible beams are negligibly small, the inertial moments Ahi are negligible, especially if the vibration frequency is not too high. If inertial moments are neglected, Eq. (17) reduces to

ð21Þ

where the damping matrix [C] is added. Popular methods for constructing a damping matrix use the concepts of modal damping and/or proportional damping from linear finite element modeling. Each mode’s modal damping ratio needs to be experimentally obtained in order to construct a damping matrix, or two parameters need to be determined in order to combine the mass and stiffness matrices into a proportional damping matrix [34]. For static problems, Eq. (21) is solved by using an incremental/iterative method based the modified Riks method [7]. For dynamic problems, Eq. (21) is solved by direct numerical integration using the Newmark-b method [7]. If q4 = q5 = q6 = 0 and/or the linear forms of dhi are adopted, the [M] and {R} are constant and only [K] is a function of

62

Y. Duan et al. / Graphical Models 75 (2013) 56–68

displacements. To derive the incremental form of Eq. (21) for static and dynamic analyses using any incremental/ iterative methods, first we define

g þ fDqðiÞ g; fqðiÞ g ¼ fq

fUg ¼ fUg þ fDUg

ð22aÞ (i)

g denotes an equilibrium state and {Dq } denotes where fq a displacement increment vector when the loads increase and/or time proceeds. Substituting Eq. (22a) into {w}, [W] and Eq. (16b) yields

 þ ½WfDUg; ½W ¼ ½W þ ½H; fwg ¼ fwg @ Wij @ 2 wi  DU k ¼ DU k @U k @U j @U k Z  ðiÞ  þ ½@NT ½HT ½Dfwg  ½k fqðiÞ g  ½@NT ½WT ½Dfwg Li  þ½@NT ½WT ½D½WfDUg ds

Hij 

ð22bÞ

By direct expansion, one can show that [7]

 ¼ ½CfDUg; ½HT ½Dfwg ¼

Cij 

@ 2 wm n Dmn w @U i @U j

@ 2 wm  n ¼ Cji Dmn w @U j @U i

^ðiÞ fDqðiÞ g; ½k fqðiÞ g ¼ ½k fqðiÞ gjfqg þ ½k Z   ^ðiÞ   ½@NT ½C þ ½WT ½D½W ½@Nds ½k ðiÞ

ð23aÞ

3.2. Collision detection

ðiÞ

ð23bÞ

Li

^ðiÞ  is the elemental tangential stiffness matrix and where ½k it is symmetric because [C] is symmetric. Since [C] requires second-order differentiations on {w}, it is more efficient to use finite difference (with, e.g., DUi = 105) in programming. The use of finite difference for [C] affects its quadratic convergence rate (because of the Newton– Raphson method) of the iteration solution process, but it does not affect the accuracy of converged solutions because this is a total-Lagrangian formulation instead of an updated-Lagrangian formulation. If nonlinear effects of dhi on the nodal loading vector {R(i)} and Ahi dhi on the mass matrix [m(i)] need to be investigated, they can be similarly treated as those shown in Eqs. (22a), (22b), (23a) and (23b). 3.1. Implementation of geometric constraints

v mðiÞ ¼ v ðjÞ n ;

ðjÞ wðiÞ m ¼ wn

Collision detection algorithms can be easily developed for this beam theory because this is a displacement-based theory and the displacements u, v, w (and hence the location) and the coordinate transformation matrix [T] (and hence the orientation) of each node is directly available from the finite element solution. For stress- or strain-based beam theories, the implementation of geometric constraints and collision detection is awkward, inaccurate, and computationally expensive. For high computational efficiency, one can compute and monitor each node’s position vector from the origin of the abc coordinate system and then check whether two elements are in contact only when the four relative position vectors between the two elements are smaller than the elemental lengths. Development of collision detection algorithm and implementation of contact constraints will be reported separately.

4. Multiple shooting formulation

Absolute constraints on displacements and rotations can be easily applied in this displacement-based beam theory. For example, if the mth node of the ith beam element is connected by a socket joint to the nth node of the jth element, we have ðiÞ um ¼ uðjÞ n ;

are needed for assembly of elements (see Fig. 3a). If relative rotations between two connected beam elements along different directions are restrained by a joint, the [T] matrix of each element needs to be calculated and related, which results in nonlinear multiple-point constraints in the finite element algorithm [30]. For two beam elements smoothly connected at a node, their u, v, w, u0 , v0 , w0 , and / are the same, and hence no coordinate transformation and multiple-point constraints are needed. Hence, one can use an initially curved small beam element to smoothly connect two oblique elements, as shown in Fig. 3c. By the way, the actual beam is more close to the one shown in Fig. 3c than the one shown in Fig. 3b. However, small elements with sizes similar to the curved connecting element need to be used around the node in order to have accurate solutions. A simpler approach is that the rigid angle between two connected beam elements can be enforced by using one or two massless rigid truss elements to connect them at locations very close to the node. Another approach similar to the use of rigid truss elements is to use multiple-point constraints on a few points on each of the two elements at locations close to the node to enforce a rigid relative angle between the two elements.

ð24Þ

which can be easily implemented in programming. On the other hand, it is difficult to implement such constraints when stress- or strain-based beam theories are used. Because u, v and w are defined w.r.t. the x, y and z axes, the exact relations between _ the_ u,_v, w, u0_, v0 , w0 , and / of _ _ _ element #1 and the u ; v ; w; u0 ; v 0 ; w0 , and / of element #2

Substituting Eqs. (9)–(12) and (17) into Eq. (20) and setting the coefficients of du, dv, dw, and dh1 to zero yields the following equations of motion [7]

8 91 8 9 8 9 8 9 €> > > < mu < F1 > < F1 > < q1 > = = = > = > C @ B ½TT F 2 A  ½k½TT F 2 þ q2 ¼ mv€ @s @ > > > > > > > : ; : ; : ; : €> ; F3 F3 q3 mw 0

M01 þ M3 q2  M 2 q3 þ q4 ¼ Ah1 ð25aÞ where F2 and F3 are inertia-induced transverse shear forces and

63

Y. Duan et al. / Graphical Models 75 (2013) 56–68

Fig. 3. Methods of connecting two beam elements: (a) an L-frame, (b) discrete connection, and (c) smooth connection.

  M03  M2 q1 þ M 1 q2 þ Ah3  q6

1 1þe  1  0 F3  M 2  M 3 q1 þ M 1 q3  Ah2 þ q5 1þe F2 

8 9 F1 > > > > Z > = 1  > M2 > > A > > > : ; M3

8 > > > <

ð25bÞ

9

r11 > > > r13 y  r12 z = dA ¼ ½Dfwg > r11 z > > > > > : ; r11 y

ð25cÞ

€  q1 Þ þ T 12 ðmv€  q2 Þ F 01 ¼ q3 F 2  q2 F 3 þ T 11 ðmu €  q3 Þ þ T 13 ðmw

The corresponding boundary conditions are:

€  q1 Þ þ T 22 ðmv€  q2 Þ F 02 ¼ q1 F 3  q3 F 1 þ T 21 ðmu €  q3 Þ þ T 23 ðmw

du ¼ 0 or F x ð F 1 T 11 þ F 2 T 21 þ F 3 T 31 Þ known dv ¼ 0 or F y ð F 1 T 12 þ F 2 T 22 þ F 3 T 32 Þ known dw ¼ 0 or F z ð F 1 T 13 þ F 2 T 23 þ F 3 T 33 Þ known dh1 ¼ 0 or M 1 known

s as the only independent variable. If these governing equations can be put into a group of nonlinear first-order ODEs, they can be solved for numerically exact solutions using a multiple shooting algorithm based on direct numerical integration using Runge–Kutta methods [7]. Because there are six boundary conditions at each end, this is a 12th-order system. The governing equations can be arranged into the following 13 first-order ODEs:

ð25dÞ

dh2 ¼ 0 or M 2 known dh3 ¼ 0 or M 3 known The governing Eqs. (25a) and (25b) can also be derived using a vector approach based on Newton’s second law and the free-body diagram of a differential beam element [7]. This shows that the energy formulation starting from the extended Hamilton principle (i.e., Eq. (20)) is fully correlated with the vector formulation, and governing equations obtained from these two different approaches are essentially the same. On the other hand, if Green–Lagrange strains and second Piola–Kirchhoff stresses are used in the extended Hamilton principle and the material stiffness matrix shown in Eq. (12) is used, one can never show that the so-obtained governing equations are the same as those from the vector formulation [7]. This geometrically exact beam theory is completely and explicitly described by u, v, w, and /. If the influences of transverse shear deformations c5 and c6 are to be included, one just needs to keep the local displacement vector U in the total displacement vector D and then follow the same derivation process [7]. However, because shear rotations are independent of bending rotations, two equations governing c5 and c6 will be added to Eqs. (25a) and (25b) [7]. For many actual applications, analysis of quasi-static deformation of flexible beams is very important. For nonlinear static problems or steady-state dynamic problems, this geometrically exact beam theory can be transformed into nonlinear ordinary differential equations (ODEs) with

€  q1 Þ þ T 32 ðmv€  q2 Þ F 03 ¼ q2 F 1  q1 F 2 þ T 31 ðmu €  q3 Þ þ T 33 ðmw M 01 ¼ q3 M 2  q2 M 3  q4 M 02 ¼ q1 M 3  q3 M 1 þ ð1 þ eÞF 3  q5 M 03 ¼ q2 M 1  q1 M 2  ð1 þ eÞF 2  q6 T 011 ¼ q3 T 21  q2 T 31 þ T 12 k3  T 13 k2 T 012 ¼ q3 T 22  q2 T 32 þ T 13 k1  T 11 k3 T 013 ¼ q3 T 23  q2 T 33 þ T 11 k2  T 12 k1

  /0 ¼ q1  ðq3 T 22  q2 T 32 þ T 13 k1  T 11 k3 ÞT 11 T 13 = 1  T 213   þ ðq3 T 21  q2 T 31 þ T 12 k3  T 13 k2 ÞT 12 T 13 = 1  T 213  T 11 k1  T 12 k2  T 13 k3 u0 ¼ 1 þ v k3  wk2 þ ð1 þ eÞT 11 v 0 ¼ wk1  uk3 þ ð1 þ eÞT 12 w0 ¼ uk2  v k1 þ ð1 þ eÞT 13 ð26aÞ

Because e and qi are functions of F1 and Mi as shown in Eqs. (25b), (25c) and (11), there are only 13 unknowns in Eq. (26a), i.e., F1, F2, F3, M1, M2, M3, T11, T12, T13, /, u, v, and w. However, T11 can be determined by using

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1  v k3 þ wk2 þ Du=Ds T 11 ¼  1  T 212  T 213  1þe

ð26bÞ

where the second expression obtained from Eq. (3) is proposed for determining the sign of T11 using spatial finite difference. Hence, only 12 equations with 12 unknowns need to be solved with 12 boundary conditions, indicating a 12th-order system.

64

Y. Duan et al. / Graphical Models 75 (2013) 56–68

The equations of T 011 ; T 012 , and T 013 are obtained from the alternative form [T]0 = [K][T]  [T][k] of Eq. (5b), and the equations for u0 , v0 , and w0 are obtained from Eq. (3). However, the use of the 13 equations in Eq. (26a) is easier for programming, assigning boundary conditions, and computation. 5. Numerical examples and discussions Two highly flexible beam structures are used to illustrate the accuracy and capability of the proposed geometrically exact beam theory and the nonlinear finite element formulation. The first example considered is the quasi-static twisting of a circular ring and the second one is the spin-up dynamics of a flexible helicopter rotor blade. 5.1. Twisting of a ring We consider a circular titanium alloy ring twisted by an angle of a3 at the two ends of a diameter, as shown in Fig. 4. The ring has 3

E ¼ 127 GPa; m ¼ 0:36; q ¼ 4430 kg=m b ¼ 3h ¼ 10 mm; Rp ¼ 800 mm

ð27Þ

The boundary conditions for the multiple shooting Eq. (26a) are given by

T 11 ¼ cos a3 ; T 12 ¼ sin a3 ; T 13 ¼ / ¼ u ¼ v ¼ w ¼ 0 at s ¼ 0 T 12 ¼ sin a3 ;

F 3 ¼ T 13 ¼ / ¼ u ¼ v ¼ 0 at s ¼ Rp

Fig. 4. Circular ring twisted by a3 at the two opposite ends of a diameter.

obtained from the exact nonlinear ordinary differential equations by direct numerical integrations. If boundary conditions are appropriately applied, the unity of T 211 þ T 212 þ T 213 ¼ 1 is automatically maintained without using the constraint shown in Eq. (26b) during the iteration process because exact coordinate transformation and nonlinear coupling terms are used in the theory. Hence, we highly recommend the multiple shooting formulation shown in Eq. (26a) because it is easy for programming, assigning boundary conditions, and computation.

ð28Þ Only one half of the ring needs to be analyzed because the deformation of the other half can be determined by using the symmetry of geometry and loading as

uð2p  hÞ ¼ uðhÞ; v ð2p  hÞ ¼ v ðhÞ; wð2p  hÞ ¼ wðhÞ; /ð2p  hÞ ¼ /ðhÞ

ð29Þ

Multiple shooting analysis is performed using the formulation shown in Eq. (26a) and 30 shooting points uniformly distributed along one half ring. Fig. 5a shows the internal forces Fi at h = p under different angles of twisting. Because the right end is free to move along the a axis, F3 = 0 at any angle of twisting. Fig. 5b shows the internal moments Mi at h = p under different angles of twisting. In Fig. 5, all the curves are highly nonlinear and have multiple turning points. When a3 = p, F1 = F2 = F3 = M1 = M3 = 0, indicating a self-locked deformed configuration. Although M2 – 0, dM2/da3 = 0 reveals that it is a self-balanced configuration. Because there is a lower half ring, the required b 3 shown in Fig. 4 should be external twisting moment M two times the internal bending moment M3 shown in Fig. 5b. The M3 curve agrees well with the numerical results obtained using a beam theory that uses two Euler angles [7] and with the experimental results obtained using a special experimental setup [7]. Fig. 6 shows different 3D deformed geometries and their 2D projections on the ab, ac, and bc planes. Note that the ring is twisted into three small rings when a3 = p. These solutions are numerically exact because they are

5.2. Spinup dynamics of a flexible rotor blade Here we consider the spinup of a 479  50.8  0.45 mm clamped-free flexible titanium rotor blade subject to a linearly increasing angular speed X = 2pt rad/s starting from its static equilibrium position, as shown in Fig. 7a with r = 0. Twenty of the beam element derived in Section 3 are used. The static vertical displacement of the rotor tip due to gravity is computed to be w = 12.6 cm, and the experimental value was w = 12.7 cm [7]. For transient analysis using direct numerical integration, the time step Dt is initially set at 0.002 s, and it is automatically adjusted by comparing the previous step’s number of iterations for convergence with a reference number of iterations (20 is used in this work). The rotor tip’s vertical displacement in Fig. 7b and the deformed configurations at different times in Fig. 7c–f shows that the spinning-induced centrifugal force raises up the rotor when the spinning speed X increases. Fig. 8 shows the time-varying frequencies of the first two intrinsic mode functions (IMFs) c1 and c2 extracted from the signal shown in Fig. 7b by using a conjugate-pair decomposition (CPD) method [35,36]. The first IMF is caused by the second and other higher-order transverse modal vibrations, and the second IMF is caused by the first transverse modal vibration [35]. The amplitude of c1 is about 0.3 mm and the amplitude of c2 is about 5 mm. The time-varying frequencies of IMFs (especially c1) are not smooth because there is no damping added in the

65

Y. Duan et al. / Graphical Models 75 (2013) 56–68

(b) 40

(a) 20 F3 0

30

F2

M i (N-m)

Fi (N)

-20

F1

-40

M3

-M 1

20

M2 10

-60 -80

0 0

50

100

150

α3 (degree)

0

50

100

150

α3 (degree)

Fig. 5. Internal forces and moments at: (a) internal forces and (b) internal moments.

Fig. 6. Different 3D deformed geometries under different twisting angles and their 2D projections on the wall. The degrees of the twisting angles are from top to bottom, then left to right: 0, 30, 60, 90, 110, 120, 130, 140, 150, 160, 170, 180, respectively.

numerical simulation and hence propagating high-frequency waves never disappear. However, Fig. 8 shows that the average frequency of c2 increases from about 1.5 Hz to 4.5 Hz and that of c1 increases from about 10 Hz to 15 Hz. These frequency ranges of c1 and c2 agree well with those numerically exact answers (the dotted lines) obtained from multiple shooting analysis with the assumption of steady-state vibrations at different constant spinning

speeds [7]. Note that this problem cannot be solved by using a geometrically exact beam derived using two Euler angles for coordinate transformation because the corresponding torsional variable is discontinuous [7,29]. Because the proposed total-Lagrangian beam element is based on a geometrically exact beam theory that accounts for any arbitrarily large rigid-elastic displacements, when it is used to model a beam attached to a moving base, there

66

Y. Duan et al. / Graphical Models 75 (2013) 56–68

z Ω

r

w (cm)

y θ η ζ ξ

x

Time (s)

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 7. Spinup of a helicopter blade with initial static deformation and gravity under a constant acceleration: (a) coordinate systems and (b) the timevarying vertical displacement of the tip end (node #21), and (c–f) deformed configurations during the first, second, fourth, and eighth cycles.

is no need of extra equation derivations in order to account for the base’s motion. Hence, this beam element is very useful for modeling and analysis of flexible multibody systems. Moreover, because the proposed beam element is fully nonlinear, it can be used to provide noise-free nonlinear response for testing, development, and validation of nonlinear system identification methods [37]. Of course, this element is very useful for natural-looking computer animation of static/dynamic deformation of 1D mechanical components.

the previous step’s solution as the initial guess of next step’s solution, the multiple shooting algorithm converges in about 10 iterations and the wall-clock time is about 25 s. For the rotor blade spinup with an average time step of 0.0012 s (Dt being automatically adjusted), the solution for each time step converges in about four iterations and the wall-clock time is about 1.5 s. The significant computational time in the finite-element analysis is expected because the code is written in MATLAB and the many nonlinear terms shown in Eqs. (3)-(11) need to be calculated for each element during each iteration.

5.3. Running time performance 5.4. Discussions The two numerical examples (Sections 5.1 and 5.2) are solved using a Dell Latitude-E6500 laptop computer with 2.0 GB RAM and a CPU speed of 2.26 GHz. For the ring twisting with an increment step of 2.5° for a3 and using

Because rotary inertias of highly flexible beams are very small, if one really wants to investigate their influences on dynamic response, it is not necessary to use a nonlinear

Y. Duan et al. / Graphical Models 75 (2013) 56–68

67

References

c1

c2

Fig. 8. Time–frequency analysis of the rotor tip’s vertical vibration using CPD.

mass matrix with those nonlinear terms shown in Eq. (17). One effective approach is to use the previous time step’s solution in Eq. (17) to derive a mass matrix and keep it as a constant matrix during the iterations for the next time step’s solution. Another way to build a database of dynamic behaviors of highly flexible 1D components is to use a camera-based motion analysis system to directly record the typical motions of such components and then put them in a database for use in computer animation [38]. This is how the game industry develops nature-looking computer games.

6. Conclusions This paper presents a geometrically exact displacement-based beam theory and its finite-element formulation. The beam theory uses three Euler angles to exactly account for the transformation between the undeformed and deformed configurations, and it solves the singularity problem of another geometrically exact beam theory based on the use of two Euler angles for coordinate transformation. Static and transient analyses of two beam structures with different boundary and loading conditions and undergoing huge overall motions are used to verify the accuracy and capabilities of the proposed beam element. Numerical results show that the proposed beam theory and finite element formulation is accurate and versatile for analysis and computer animation of multibody systems with highly flexible beams. Acknowledgments The authors thank for the support by the National Science Foundation under Grants CMS 0319853, CMMI0856206 and CMMI-1039433. The authors thank for the support by the National Science Foundation under Grants CMS 0319853, CMMI-0856206, CMMI-1039433, and CCNIE 1245795.

[1] R. Barzel, Faking dynamics of ropes and springs, IEEE Computer Graphics and Applications 17 (3) (1997) 31–39. [2] D.K. Pai, Interactive simulation of thin solids using Cosserat models, Computer Graphics Forum 21 (3) (2002) 347–352. [3] J. Kautz, T.Y. Lee, M.C. Lin, Stable robust elastic rod simulation using heterogeneous materials, Pacific Graphics 30 (7) (2011). [4] S.S. Antman, Nonlinear Problems of Elasticity, Springer-Verlag, 1995. [5] M.B. Rubin, Cosserat Theories: Shells, Rods and Points, Kluwer, 2000. [6] C.J. Gantes, Deployable Structures: Analysis and Design, WIT Press, Southampton, England, UK, 2001. [7] P.F. Pai, Highly Flexible Structures: Modeling, Computation, and Experimentation, AIAA, Reston, Virginia, 2007. [8] C.H.M. Jenkins (Ed.), Gossamer Spacecraft: Membrane and Inflatable Structures Technology for Space Applications, AIAA, Reston, Virginia, 2001. [9] E. Sporn, B.W. Miedema, S.L. Bachman, T.S. Loy, R.D. Calaluce, K. Thaler, Endoscopic colotomy closure after full thickness excision: comparison of T-fastener versus multi-clip applier, Endoscopy 40 (7) (2008) 589–594. [10] B.W. Miedema, J.A. Astudillo, E. Sporn, K. Thaler, NOTES techniques: present and future, Review, European Surgery 40 (3) (2008) 103– 110. [11] Y. Yang, I. Tobias, W.K. Olson, Finite element analysis of DNA supercoiling, Journal of Chemical Physics 98 (1993) 1673– 1686. [12] Y. Shi, J.E. Hearst, The Kirchhoff elastic rod, the nonlinear Schrodinger equation and DNA supercoiling, Journal of Chemical Physics 101 (1994) 5186–5200. [13] T. Schlick, Modeling superhelical DNA: recent analytical and dynamical approaches, Current Opinions in Structural Biology 5 (1995) 245–262. [14] D.F. Parker, The role of Saint Venant’s solutions in rod and beam theories, Journal of Applied Mechanics 46 (1979) 861–866. [15] M. Borri, T. Merlini, A large displacement formulation for anisotropic beam analysis, Meccanica 21 (1986) 30–37. [16] D.H. Hodges, Nonlinear Composite Beam Theory, AIAA, Reston, Virginia, 2006. [17] A.A. Shabana, Dynamics of Multibody Systems, third ed., Cambridge University Press, New York, 2005. [18] O.A. Bauchau, Flexible Multibody Dynamics, Springer Verlag, New York, 2010. [19] T.R. Kane, R.R. Ryan, A.K. Banerjee, Dynamics of a cantilever beam attached to a moving base, Journal of Guidance and Control 10 (1987) 139–151. [20] C.C. Rankin, F.A. Brogan, An element-independent corotational procedure for the treatment of large rotations, Journal of Pressure Vessel Technology 108 (1986) 165–174. [21] F. Deng, X. He, L. Li, J. Zhang, Dynamics modeling for a rigid-flexible coupling system with nonlinear deformation field, Multibody System Dynamics 18 (2007) 559–578. [22] S.C. Wu, E.J. Haug, Geometric nonlinear substructuring for dynamics of flexible mechanical systems, International Journal for Numerical Methods in Engineering 26 (1988) 2211–2226. [23] H.H. Yoo, R.R. Ryan, R.A. Scott, Dynamics of flexible beams undergoing large overall motions, Journal of Sound and Vibration 181 (2) (1995) 261–278. [24] I. Sharf, Nonlinear strain measures, shape functions and beam elements for dynamics of flexible beams, Multibody System Dynamics 3 (1999) 189–205. [25] J. Mayo, J. Dominguez, A.A. Shabana, Geometrically nonlinear formulations of beams in flexible multibody dynamics, Journal of Vibration and Acoustic 117 (1995) 501–509. [26] J.C. Simo, L. Vu-Quoc, On the dynamics of flexible beams under large overall motions – the plane case: Part I, Journal of Applied Mechanics 53 (1986) 849–854. [27] J.C. Simo, L. Vu-Quoc, On the dynamics of flexible beams under large overall motions – the plane case: Part II, Journal of Applied Mechanics 53 (1986) 855–863. ´ nguez, Study of the geometric [28] J.M. Mayo, D. Garcı´a-Vallejo, J. Domı stiffening effect: comparison of different formulations, Multibody System Dynamics 11 (2004) 321–341. [29] P.F. Pai, Three kinematic representations for modeling of highly flexible beams and their applications, International Journal of Solids and Structures 48 (2011) 2764–2777. [30] P.F. Pai, Geometrically exact beam theory without Euler angles, International Journal of Solids and Structures 48 (2011) 3075–3090.

68

Y. Duan et al. / Graphical Models 75 (2013) 56–68

[31] P.F. Pai, T.J. Anderson, E.A. Wheater, Large-deformation tests and total-lagrangian finite-element analyses of flexible beams, International Journal of Solids and Structures 37 (2000) 2951– 2980. [32] G. Wu, X. He, P.F. Pai, Geometrically exact total-Lagrangian modeling and finite element analysis of highly flexible multibody systems, in: 51st AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference, Orlando, Florida, April 12–15, 2010. [33] S.P. Timoshenko, J.N. Goodier, Theory of Elasticity, third ed., McGraw-Hill, New York, 1970. [34] D.J. Ewins, Modal Testing: Theory, Practice and Application, second ed., Research Studies Press, New York, 2000.

[35] N.E. Huang, N.O. Attoh-Okine (Eds.), The Hilbert–Huang Transform in Engineering, CRC Press, Boca Raton, FL, 2005. [36] P.F. Pai, Online tracking of instantaneous frequency and amplitude of dynamical system response, Mechanical Systems and Signal Processing 24 (4) (2010) 1007–1024. [37] P.F. Pai, A.N. Palazotto, Detection and identification of nonlinearities by amplitude and frequency modulation analysis, Journal of Mechanical Systems and Signal Processing 22 (2008) 1107–1132. [38] P.F. Pai, S. Ramanathan, J. Hu, D.K. Chernova, X. Qian, G. Wu, Camerabased noncontact metrology for static/dynamic testing of flexible multibody systems, Measurement Science and Technology 21 (8) (2010) 085302 (14p).