- Email: [email protected]
Nonlinear dynamics of elastic rods using the Cosserat theory: Modelling and simulation
Available online at www.sciencedirect.com
International Journal of Solids and Structures 45 (2008) 460–477 www.elsevier.com/locate/ijsolstr
Nonlinear dynamics of elastic rods using the Cosserat theory: Modelling and simulation D.Q. Cao a
a,*,1
, Robin W. Tucker
b
School of Astronautics, Harbin Institute of Technology, P.O. Box 137, Harbin 150001, China b Department of Physics, Lancaster University, Lancaster LA1 4YB, UK Received 20 March 2006; received in revised form 9 July 2007 Available online 6 September 2007
Abstract The method of Cosserat dynamics is employed to explore the nonplanar nonlinear dynamics of elastic rods. The rod, which is assumed to undergo flexure about two principal axes, extension, shear and torsion, are described by a general geometrically exact theory. Based on the Cosserat theory, a set of governing partial differential equations of motion with arbitrary boundary conditions is formulated in terms of the displacements and angular variables, thus the dynamical analysis of elastic rods can be carried out rather simply. The case of doubly symmetric cross-section of the rod is considered and the Kirchoff constitutive relations are adopted to provide an adequate description of elastic properties in terms of a few elastic moduli. A cantilever is given as a simple example to demonstrate the use of the formulation developed. The nonlinear dynamic model with the corresponding boundary and initial conditions are numerically solved using the Femlab/Matlab software packages. The corresponding nonlinear dynamical responses of the cantilever under external harmonic excitations are presented through numerical simulations. 2007 Elsevier Ltd. All rights reserved. Keywords: Cosserat model; Nonlinear dynamics; Slender rod; Modelling and simulation; Femlab; Matlab
1. Introduction Nonlinear dynamic analysis of slender elastic structures under external forces and torques and parametric excitations remains an active area of study. Such a study can find application in accelerating missiles and space crafts, components of high-speed machinery, manipulator arm, microelectronic mechanical structures (MEMS), components of bridges (such as towers and cables) and other structural elements. Considerable attention has been devoted to the study of nonlinear dynamics of rods or beams subject to both external and parametric excitations (see, for example, Nayfeh and D.T.Mook, 1979; Saito and Koizumi, 1982; Kar and Dwivedy, 1999 and the references cited therein). While attention so far has mainly been devoted to the *
1
Corresponding author. Tel.: +86 (0) 451 86414479; fax: +86 (0) 451 86402822. E-mail addresses: [email protected] (D.Q. Cao), [email protected] (R.W. Tucker). On leave from the Department of Physics, Lancaster University, Lancaster LA1 4YB, UK.
0020-7683/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2007.08.016
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
461
study of planar, nonlinear dynamic analysis of beams, research has been done concerning nonplanar, nonlinear motions of beams. Bolotin (1964) addressed such motions in beams, but restricted himself to consideration of nonlinear inertia terms and stability of the planar response. Crespo da Silva and Glynn (1978a) formulated the equations of motion describing the nonplanar, nonlinear dynamics of an inextensional beam. The nonplanar, nonlinear forced oscillations of a cantilever are then analyzed in Crespo da Silva and Glynn (1978b) using the perturbation method. Cartmell (1990) and Forehand and Cartmell (2001) derived the nonlinear equation of motion for the in-plane and out-of-plane forced vibration of cantilever beams with a lumped mass. The studies mentioned above are restricted to systems with no extension of the beam’s neutral axis and no warping or shear deformation. Crespo da Silva (1988a,b) investigated the problem of nonlinear dynamics of the nonplanar flexural-torsional-extensional beams, but the effects of rotary inertia and shear deformation were neglected. Thus, cross-sectional dimensions of the beam were assumed to be small enough in comparison to the beam length. However, the shear deformation may be of considerable importance and can not be negligible for studying the vibration of high frequencies when a comparative short rod is investigated. In such a case, the effect of shear deformation should be taken into account for. For the planar problem, although such effects can be included by using the Timoshenko beam theory, most of the studies are limited to the determination of natural frequencies and eigenfunctions (Huang, 1961; Bishop and Price, 1977; Grant, 1978; Bruch and Mitchell, 1987). For the three-dimensional problem, we refer to Simo and Vu-Quoc (1988, 1991) for a geometrically exact rod model incorporating shear and torsion-warping deformation. Further studies on the dynamic formulation of sandwich beams have been presented in Vu-Quoc and Deng (1995) and Vu-Quoc and Ebcioglu (1995) based on the geometrically-exact description of the kinematics of deformation. Moreover, Esmailzadeh and Jalili (1998) investigated the parametric response of cantilever Timoshenko beams with lumped mass, but restricted themselves to consideration of nonlinear inertia terms. Based on the geometrically-exact model of sliding beams, parametric resonance has been presented in Vu-Quoc and Li (1995), where the beams can undergo large deformation with shear deformation accounted for. In the case involving the full dynamic response, the strong nonlinearity and fully coupling introduce a challenge for solving the partial differential equations of motion of rods. We refer to Rubin (2001) for a formulation of a numerical solution procedure for three-dimensional dynamic analysis of rods by modelling the rod as a set of connected Cosserat points, also Rubin and Tufekci (2005) and Rubin (2000) for three-dimensional dynamics of a circular arch and shells using the theory of a Cosserat point, respectively. A Galerkin projection has been applied to discretize the governing partial differential equations of sliding beams in Vu-Quoc and Li (1995) and sandwich beams in Vu-Quoc and Deng (1997). Recently, a new modelling strategy has been proposed in Cao et al. (2006) to discretize the rod and to derive the ordinary differential equations of motion with third order nonlinear generic nodal displacements. This modelling strategy has been successfully used to investigate the nonlinear dynamics of typical MEMS device that comprises a resonator mass supported by four flexible beams (Cao et al., 2005). In this paper we explore the nonplanar, nonlinear dynamics of an extensional shearable rod by using the simple Cosserat model. The method of Cosserat dynamics for elastic structures is employed since it can accommodate to a good approximation the nonlinear behavior of complex elastic structures composed of materials with different constitutive properties, variable geometry and damping characteristics (Antman, 1991; Antman et al., 1998; Tucker and Wang, 1999; Cull et al., 2000; Gratus and Tucker, 2003). With arbitrary boundary conditions, the Cosserat theory is used to formulate a set of governing partial differential equations of motion in terms of the displacements and angular variables, describing the nonplanar, nonlinear dynamics of an extensional rod. Bending about two principal axes, extension, shear and torsion are considered, and care is taken into account for all the nonlinear terms in the resulting equations. As an example, a simple cantilever is given to demonstrate the use of the formulation developed. The nonlinear dynamic model with the corresponding boundary and initial conditions are numerically solved using the commercially available software packages Femlab and Matlab. Corresponding nonlinear dynamical responses of the cantilever under external harmonic excitations are presented and discussed through numerical simulations. The following conventions and nomenclature will be used through out this paper. Vectors, which are elements of Euclidean 3-space R3 , are denoted by lowercase, bold-face symbols, e.g., u, v; vector-valued functions are denoted by lowercase, italic, bold-face symbols, e.g., u, v; tensors are denoted by upper-case, bold-face symbols, e.g., I, J; matrices are denoted by upper-case, italic, bold-face symbols, e.g., M, K. The symbols
462
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
ot and os denote differentiation with respect to time t and arc-length parameter s, respectively. The symbols ð : Þ and ( 0 ) denote differentiation with respect to dimensionless time parameter s and dimensionless length parameter r, respectively. 2. Background on the simple Cosserat model In this section the basic concepts of the Cosserat theory for an elastic rod of unstressed length ‘ are summarized. The three vectors {e1, e2, e3} are assumed to form a fixed right-handed orthogonal basis. Elements of the rod are labelled in terms of the Lagrangian coordinate 06 s 6 ‘ at time t. The motion of an elastic segment may be described in terms of the motion in space of a vector r that locates the line of centroids (shown dashdotted) of the cross-sections as shown in Fig. 1. Specifying a unit vector d3 (which may be identified with the normal to the cross-section) at each point along this line enables the state of shear to be related to the angle between this vector and the tangent osr to the centroid space-curve. Specifying a second vector d2 orthogonal to the first vector (thereby placing it in the plane of the cross-section) can be used to encode the state of bending and twist along the element. In such a way, the current orientation of each cross-section at s 2 [0, ‘] is defined by specifying the orientation of a moving basis d1, d2, d3 = d1 · d2, see Simo and Vu-Quoc (1988) and Vu-Quoc and Deng (1995) for details. Elastic deformations about the line of centroids are then coded into the rates of change of r and the triad d1, d2, d3 of the cross-section at s. Thus a time dependent field of two mutually orthogonal unit vectors along the segment provides three continuous dynamical degrees of freedom that, together with the three continuous degrees of freedom describing the centroid space-curve relative to some arbitrary origin in space (with fixed inertial frame e1, e2, e3), define a simple Cosserat rod model. Supplemented with appropriate constitutive relations and boundary conditions such a model can fully accommodate the modes of vibration that are traditionally associated with the motion of slender rods in the engineering literature: namely axial motion along its length, torsional or rotational motion and transverse or lateral motion. The tangent to the space-curve, v = v(s, t), is given by vðs; tÞ ¼ os rðs; tÞ;
ð1Þ
the rotation of the directions along the space curve and the local angular velocity vector of the director triad are measured by u = u(s, t), w = w(s, t) according to os di ðs; tÞ ¼ uðs; tÞ di ðs; tÞ; ð2Þ ot di ðs; tÞ ¼ wðs; tÞ di ðs; tÞ; respectively. It follows from the first equation of (2) that
s=b
s
r
d3
d2
e3 r
s d1
o
e2
e1 Fig. 1. A simple Cosserat model.
s=a
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477 3 3 3 X X X ðdi os di Þ ¼ ðdi ðu di ÞÞ ¼ ðuðdi di Þ di ðdi uÞÞ ¼ 2u: i¼1
i¼1
463
ð3Þ
i¼1
Similarly, we have 3 3 3 X X X ðdi ot di Þ ¼ ðdi ðw di ÞÞ ¼ ðwðdi di Þ di ðdi wÞÞ ¼ 2w: i¼1
i¼1
ð4Þ
i¼1
Since the basis {d1, d2, d3} is natural for the intrinsic description of deformation, we decompose relevant vector valued functions with respect to it: 8 3 P > > > vðs; tÞ ¼ vi ðs; tÞdi ðs; tÞ; > > > i¼1 > > < 3 P ð5Þ uðs; tÞ ¼ ui ðs; tÞdi ðs; tÞ; > i¼1 > > > > 3 > P > > wi ðs; tÞdi ðs; tÞ: : wðs; tÞ ¼ i¼1
The contact forces n(s, t) and contact torques m(s, t) are related to the extension and shear strains v(s, t), and flexure and torsion strains u(s, t) by constitutive relations, respectively. The angular momentum h(s, t) is related to the rotary inertia and the director angular velocity w(s, t). The dynamical evolution of the rod with density, q(s), and cross-section area, A(s) is governed by the Newton’s dynamical laws: qðsÞAðsÞott r ¼ os nðs; tÞ þ fðs; tÞ; ð6Þ ot hðs; tÞ ¼ os mðs; tÞ þ vðs; tÞ nðs; tÞ þ lðs; tÞ; where f(s, t) and l(s, t) denote external force and torque densities, respectively, the contact force vector, contact torque vector, and the angular momentum vector can be written as 8 3 P > > > ni ðs; tÞdi ðs; tÞ; nðs; tÞ ¼ < i¼1 ð7Þ 3 > P > > mi ðs; tÞdi ðs; tÞ; : mðs; tÞ ¼ i¼1
and hðs; tÞ ¼
3 X
hi ðs; tÞdi ðs; tÞ;
ð8Þ
i¼1
respectively. For more details about Eqs. (6)–(8), we refer to the governing equations (i)–(iii) derived in terms of the noncommutative Lie group of proper orthogonal transformations in Simo and Vu-Quoc (1988), see also in Vu-Quoc and Deng (1995) and Vu-Quoc and Ebcioglu (1995). It will be assumed that the Young’s modulus E, the shear modulus G, and the material density q along the rod are functions of s only, and the mass center coincides with the area centroid of the cross-section at s. In this case, the simplest constitutive model for the linear elastic material is based on the Kirchhoff constitutive relations with shear deformation, see the discussions in Simo and Vu-Quoc (1988, 1991) for details of constitutive relations. Such a model provide an adequate description of elastic properties in terms of a few elastic moduli. The presence of arbitrary rotations relating the local director frame to the global inertial frame renders the equations of motion inherently nonlinear. The contact forces, contact torques and the angular momentum are then given as n ¼ Kðv d3 Þ;
m ¼ Ju;
h ¼ Iw;
where, for a symmetric cross-section of a rod, the tensors K, J and I are described as
ð9Þ
464
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
8 3 P > > > Kðs; tÞ ¼ K ii ðs; tÞðdi ðs; tÞ di ðs; tÞÞ; > > > i¼1 > > < 3 P Jðs; tÞ ¼ J ii ðs; tÞðdi ðs; tÞ di ðs; tÞÞ; > i¼1 > > > > 3 > P > > Iðs; tÞ ¼ I ii ðs; tÞðdi ðs; tÞ di ðs; tÞÞ: :
ð10Þ
i¼1
Here the symbol is used to denote the tensorial product u v of two vectors u and v that assigns to each vector w the vector (u v)w = (v Æ w)u. The corresponding terms in (10) can be written as 8 K 11 ¼ K 22 ¼ jðsÞGðsÞAðsÞ; K 33 ¼ EðsÞAðsÞ; > > R R > 2 > > J 11 ¼ AðsÞ EðsÞg dA; J 22 ¼ AðsÞ EðsÞn2 dA; > > > R R < I 11 ¼ AðsÞ qðsÞg2 dA; I 22 ¼ AðsÞ qðsÞn2 dA; ð11Þ > R > 2 oW oW 2 > > J 33 ¼ AðsÞ GðsÞ n þ g þ n og þ g on dA; > > > > R : I 33 ¼ AðsÞ qðsÞðn2 þ g2 Þ dA; where W is the warping function which is the same for all cross-sections according to the Saint-Venant Theory, and j(s) is a numerical factor depending on the shape of the cross-section at s. Making use of (5), (9) and (10), we have 8 > < n ¼ K 11 v1 d1 þ K 22 v2 d2 þ K 33 ðv3 1Þd3 ; m ¼ J 11 u1 d1 þ J 22 u2 d2 þ J 33 u3 d3 ; ð12Þ > : h ¼ I 11 w1 d1 þ I 22 w2 d2 þ I 33 w3 d3 : Comparing with (7), we get 8 m1 ¼ J 11 u1 ; > < n1 ¼ K 11 v1 ; n2 ¼ K 22 v2 ; m2 ¼ J 22 u2 ; > : n3 ¼ K 33 ðv3 1Þ; m3 ¼ J 33 u3 ;
h1 ¼ I 11 w1 ; h2 ¼ I 22 w2 ;
ð13Þ
h3 ¼ I 33 w3 :
3. Specifications for the deformed configuration space Consider a uniform and initially straight rod of constant length ‘, supported in an arbitrary manner at each end as shown in Fig. 2. It is assumed that the static equilibrium of the rod corresponds to the situation where the directions of d3 and e3 coincide, d1 and d2 are parallel to e1 and e2, respectively. The components of the elastic displacement vector of the centroid at an arbitrary location s is assumed to be x(s, t), y(s, t), z(s, t), i.e. d1
e1
d3
s
r O
d2 A
e3 e2 Fig. 2. Coordinate systems used in the development of governing equations.
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
rðs; tÞ ¼ xðs; tÞe1 þ yðs; tÞe2 þ ðs þ zðs; tÞÞe3 :
465
ð14Þ
The orientation of the directors d1, d2, d3 at location s, relative to the inertial basis {e1, e2, e3}, can be described by three successive rotations, called Euler angles. Alternatively, the moving frame can be described by the rotation group SO(3) as the formulations in Simo and Vu-Quoc (1988, 1991), Vu-Quoc and Deng (1995) and Vu-Quoc and Ebcioglu (1995). The three successive rotations starts by rotating the basis {e1, e2, e3} an angle w(s, t) about the axis aligned with the director e2 as shown in Fig. 3a. Next, we rotate the basis f^e1 ; e2 ; ^e3 g an angle h(s, t) about the axis aligned with the new director ^e1 as shown in Fig. 3b. Finally, we rotate the basis f^e1 ; ^e2 ; d3 g an angle /(s, t) about the axis aligned with the new director d3 as shown in Fig. 3c. The representative matrices for the three successive rotations can be written as 2 3 cos wðs; tÞ 0 sin wðs; tÞ 6 7 Dw ðs; tÞ ¼ 4 0 1 0 ð15Þ 5; sin wðs; tÞ
2
1 6 Dh ðs; tÞ ¼ 4 0 0 and
0
0 cos hðs; tÞ sin hðs; tÞ
cos wðs; tÞ
3 0 7 sin hðs; tÞ 5;
ð16Þ
cos hðs; tÞ
2
cos /ðs; tÞ sin /ðs; tÞ 6 D/ ðs; tÞ ¼ 4 sin /ðs; tÞ cos /ðs; tÞ 0
0
3 0 7 0 5;
ð17Þ
1
respectively. Denote dij = di Æ ej. Then, using Eqs. (15)–(17), one can write the transformation matrix between the inertial basis {e1, e2, e3} and the basis {d1, d2, d3} as follows: Dðs; tÞ ¼ ½d ij ðs; tÞ ¼ D/ ðs; tÞDh ðs; tÞDw ðs; tÞ 2 cos / cos w þ sin / sin h sin w sin / cos h 6 ¼ 4 sin / cos w þ cos / sin h sin w cos / cos h
3 cos / sin w þ sin / sin h cos w 7 sin / sin w þ cos / sin h cos w 5:
sin h
cos h sin w
cos h cos w
Since D is an orthogonal matrix, from the definition of dij = di Æ ej, we have
e1 e1
e1
ψ
e1
e3 ψ
d1 o (e 2 )
φ ψ
o
θ
φ
e3
d3 θ ψ
e2
(a)
e3
d3
e3
(b)
θ
e2
e3
θ
o(e1 ) d2
φ
e2
e2
(c) d1
e2
φ
d2 o (d 3)
e1
Fig. 3. The successive angle rotations of a differential element of arc length ds.
ð18Þ
466
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
di ¼
3 X
ej ¼
d ij ej ;
j¼1
ð19Þ
d kj dk :
k¼1
Thus os di ¼ os
3 X
3 X
! d ij ej
j¼1
¼
3 X
os d ij ej ¼
3 X 3 X
j¼1
j¼1
ðos d ij Þd kj dk :
Substituting (20) into (3) yields 3 3 1X 1X di os di ¼ di u¼ 2 i¼1 2 i¼1
¼
3 X 3 X j¼1
ð20Þ
k¼1
! ðos d ij Þd kj dk
k¼1
3 3 1X 1X ½ðos d 1j Þd 2j ðos d 2j Þd 1j d3 þ ½ðos d 2j Þd 3j ðos d 3j Þd 2j d1 2 j¼1 2 j¼1
þ
3 1X ½ðos d 3j Þd 1j ðos d 1j Þd 3j d2 : 2 j¼1
ð21Þ
Taking notice of 3 X
d ij ðs; tÞd kj ðs; tÞ ¼ dik ;
ð22Þ
j¼1
where d is the Kronecker d-function, we have 3 X
½ðos d ij ðs; tÞÞd kj ðs; tÞ þ d ij ðs; tÞðos d kj ðs; tÞÞ ¼ 0:
ð23Þ
j¼1
Hence, (21) can be written as u¼
3 X ½ðos d 2j Þd 3j d1 þ ðos d 3j Þd 1j d2 þ ðos d 1j Þd 2j d3 :
ð24Þ
j¼1
Therefore, 8 > < u1 u2 > : u3
from (18) and (24) we can obtain ¼ os h cos / þ os w cos h sin /; ð25Þ
¼ os h sin / þ os w cos h cos /; ¼ os / os w sin h:
Similarly, replacing the time derivative by the space derivative, we can easily obtain 8 > < w1 ¼ ot h cos / þ ot w cos h sin /; w2 ¼ ot h sin / þ ot w cos h cos /; > : w3 ¼ ot / ot w sin h:
ð26Þ
In order to express the extension and shear strains as functions of displacements x(s, t), y(s, t), z(s, t) and angles h(s, t), w(s, t), /(s, t), we define re1 ðs; tÞ ¼ xðs; tÞ;
re2 ðs; tÞ ¼ yðs; tÞ;
re3 ðs; tÞ ¼ s þ zðs; tÞ:
ð27Þ
Then, the strain vector v(s, t) can be written as v ¼ os r ¼
3 3 X X ðos rei Þei ¼ ðos rei Þd ji dj : i¼1
i¼1
Substituting (18) and (27) into (28) yields
ð28Þ
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
8 v1 ¼ os xðcos / cos w þ sin / sin h sin wÞ þ os y sin / cos h > > > > > þð1 þ os zÞð cos / sin w þ sin / sin h cos wÞ; > < v2 ¼ os xð sin / cos w þ cos / sin h sin wÞ þ os y cos / cos h > > > þð1 þ os zÞðsin / sin w þ cos / sin h cos wÞ; > > > : v3 ¼ os x cos h sin w os y sin h þ ð1 þ os zÞ cos h cos w:
467
ð29Þ
4. Governing equations of motion To obtain the partial differential equations of motion in terms of the displacements x(s, t), y(s, t), z(s, t) and angles h(s, t), w(s, t), /(s, t), we decompose the contact force n(s, t) with respect to the fixed basis {e1, e2, e3} as nðs; tÞ ¼
3 X
nei ðs; tÞei :
ð30Þ
i¼1
Then, making use of (9) and (19), we have n¼
3 X
ni di ¼
i¼1
3 X 3 X i¼1
3 X
ni d ij ej ¼
j¼1
nej ej :
ð31Þ
i¼1
It follows from (13) that 8 3 P > > > ne1 ¼ d i1 ni ¼ d 11 K 11 v1 þ d 21 K 22 v2 þ d 31 K 33 ðv3 1Þ; > > > i¼1 > > < 3 P ne2 ¼ d i2 ni ¼ d 12 K 11 v1 þ d 22 K 22 v2 þ d 32 K 33 ðv3 1Þ; > i¼1 > > > > 3 > P > > d i3 ni ¼ d 13 K 11 v1 þ d 23 K 22 v2 þ d 33 K 33 ðv3 1Þ: : ne3 ¼
ð32Þ
i¼1
To proceed, using the relations between the strain vector u(s, t) and the differentiation of the directors d1(s, t), d2(s, t) and d3(s, t) with respect to s, and the relations between the angular velocity vector w(s, t) and the differentiation of the directors d1(s, t), d2(s, t) and d3(s, t) with respect to t, we can easily obtain the following expressions: os m ¼
3 X
3 X
os ðmi di Þ ¼
i¼1
ðos mi di þ mi u di Þ
i¼1
¼ ðos m1 m2 u3 þ m3 u2 Þd1 þ ðos m2 m3 u1 þ m1 u3 Þd2 þ ðos m3 m1 u2 þ m2 u1 Þd3 ; ot h ¼
3 X
ð33Þ
3 X ot ðhi di Þ ¼ ðot hi di þ hi w di Þ
i¼1
i¼1
¼ ðot h1 h2 w3 þ h3 w2 Þd1 þ ðot h2 h3 w1 þ h1 w3 Þd2 þ ðot h3 h1 w2 þ h2 w1 Þd3 : In addition, we record ! 3 X vn¼ v i di i¼1
3 X j¼1
ð34Þ
! nj dj
¼ ðv2 n3 v3 n2 Þd1 þ ðv3 n1 v1 n3 Þd2 þ ðv1 n2 v2 n1 Þd3 :
ð35Þ
468
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
Now, let us assume that the external force and torque are distributed loadings with fixed direction and prescribed intensity. Consequently, the external force and torque densities can be written as: fðs; tÞ ¼ fx ðs; tÞe1 þ fy ðs; tÞe2 þ fz ðs; tÞe3 ; ð36Þ lðs; tÞ ¼ lx ðs; tÞe1 þ ly ðs; tÞe2 þ lz ðs; tÞe3 : Decomposing the torque density with respect to the moving basis {d1, d2, d3}, we have lðs; tÞ ¼
3 X
li ðs; tÞdi ðs; tÞ;
ð37Þ
i¼1
where, li(s,t) (i = 1, 2, 3) can be expressed by using the transformation matrix (18) and the relation (19) as 8 l1 ¼ lx ðcos / cos w þ sin / sin h sin wÞ þ ly sin / cos h > > > > > þlz ð cos / sin w þ sin / sin h cos wÞ; > < l2 ¼ lx ð sin / cos w þ cos / sin h sin wÞ þ ly cos / cos h ð38Þ > > > þlz ðsin / sin w þ cos / sin h cos wÞ; > > > : l3 ¼ lx cos h sin w ly sin h þ lz cos h cos w: Substituting the expressions (32)–(37) into the dynamic equation (6), we can obtain 8 qAott x ¼ os ðd 11 K 11 v1 þ d 21 K 22 v2 þ d 31 K 33 ðv3 1ÞÞ þ fx ; > > > > > qAo tt y ¼ os ðd 12 K 11 v1 þ d 22 K 22 v2 þ d 32 K 33 ðv3 1ÞÞ þ fy ; > > > < qAo z ¼ o ðd K v þ d K v þ d K ðv 1ÞÞ þ f ; tt s 13 11 1 23 22 2 33 33 3 z > ot h1 h2 w3 þ h3 w2 ¼ os m1 m2 u3 þ m3 u2 þ v2 n3 v3 n2 þ l1 ; > > > > > ot h2 h3 w1 þ h1 w3 ¼ os m2 m3 u1 þ m1 u3 þ v3 n1 v1 n3 þ l2 ; > > : ot h3 h1 w2 þ h2 w1 ¼ os m3 m1 u2 þ m2 u1 þ v1 n2 v2 n1 þ l3 ;
ð39Þ
where dij is defined by (18), ni, mi and hi are given by (13), while ui, wi, vi and li are given by (25), (26), (29) and (38), respectively. The boundary conditions at the ends of the rod should be determined in each particular case. For example, at a fixed end the deflections and the rotations are equal to zero. In this case the boundary conditions are x¼y¼z¼0
and
h ¼ w ¼ / ¼ 0:
ð40Þ
At a free end the boundary conditions are am n þ f b ðtÞ ¼ 0
and
am m þ lb ðtÞ ¼ 0;
where am = +1 at the right end of the rod, am = 1 at the left end, and ( b f ðtÞ ¼ fxb ðtÞe1 þ fyb ðtÞe2 þ fzb ðtÞe3 ; lb ðtÞ ¼ lbx ðtÞe1 þ lby ðtÞe2 þ lbz ðtÞe3
ð41Þ
ð42Þ
are the external force and torque applied to the corresponding end of the rod. As the external distributed torque l(s, t), the torque vector lb(t) can be expressed in the moving frame {d1, d2, d3} and the corresponding components lbi ðtÞ ði ¼ 1; 2; 3Þ can be determined through the relation (19) and the transformation matrix (18) as follows: 8 b l1 ¼ lbx ðcos / cos w þ sin / sin h sin wÞ þ lby sin / cos h > > > > > > þlbz ð cos / sin w þ sin / sin h cos wÞ; > < b l2 ¼ lbx ð sin / cos w þ cos / sin h sin wÞ þ lby cos / cos h ð43Þ > > > b > þlz ðsin / sin w þ cos / sin h cos wÞ; > > > : lb ¼ lb cos h sin w lb sin h þ lb cos h cos w; 3
x
y
z
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
469
with w = w(s*, t), h = h(s*, t) and / = /(s*, t) at the boundary s = s* (s* = 0 or ‘). Substituting the relations (32) and (42) into the first equation of (41), the force boundary conditions can be rewritten as 8 b > < am ½d 11 K 11 v1 þ d 21 K 22 v2 þ d 31 K 33 ðv3 1Þ þ fx ðtÞ ¼ 0; am ½d 12 K 11 v1 þ d 22 K 22 v2 þ d 32 K 33 ðv3 1Þ þ fyb ðtÞ ¼ 0; ð44Þ > : b am ½d 13 K 11 v1 þ d 23 K 22 v2 þ d 33 K 33 ðv3 1Þ þ fz ðtÞ ¼ 0: Similarly, the torque boundary conditions can be obtained as am J 11 u1 þ lb1 ðtÞ ¼ 0;
am J 22 u2 þ lb2 ðtÞ ¼ 0;
am J 33 u3 þ lb3 ðtÞ ¼ 0:
ð45Þ
For simplicity, it will be assumed that the axes along the directors d1, d2 and d3 are chosen to be the principal axes of inertia of the cross-section at s, and centered at the cross-section’s center of mass. Then, we have 8 > < K 11 ¼ jGA; K 22 ¼ jGA; K 33 ¼ EA; J 11 ¼ Ec11 ; J 22 ¼ Ec22 ; J 33 ¼ Gc33 ; ð46Þ > : I 11 ¼ qc11 ; I 22 ¼ qc22 ; I 33 ¼ qðc11 þ c22 Þ; where c11 and c22 are the moments of inertia of the rod cross-section, and c33 is the torsional moment of inertia. For a rectangular cross-section, the torsional moment of inertia is as follows (Rubin, 2000): " # 1 WH 3 192H X 1 pð2n 1ÞW c33 ¼ 1 5 tanh for W P H ; p W n¼1 ð2n 1Þ5 2H 3 where W and H are the dimensions of the cross-section. Moreover, the approximate formula of the torsional moment of inertia can be written as WH 3 H H4 1 1 0:63 c33 ¼ for W P H : W 3 12W 4 As a prelude to expanding the nonlinear partial differential equations (39) to a form suitable for a perturbation analysis of the motion, it is useful to introduce some natural scales to obtain a dimensionless equation of motion. Introduce the dimensionless variables s r¼ ; ‘
x p1 ¼ ; ‘
y p2 ¼ ; ‘
z p3 ¼ ; ‘
p4 ¼ h;
p5 ¼ w;
p6 ¼ /;
s ¼ x0 t;
where x0 is the reference natural frequency yet to be determined. Moreover, let 8 0 0 g0 ¼ qAx1 2 ‘ ; n01 ¼ jGAg ; n02 ¼ EAg ; > ‘ ‘ > > 0 > > > > g3 ¼ qc 1 x2 ; g1 ¼ qc 1 x2 ; g2 ¼ qc 1 x2 ; > > 11 0 22 0 33 0 > < n11 ¼ q‘2Ex2 ; n12 ¼ Gcq‘332 cEcx222 ; n13 ¼ jGAg1 ; n14 ¼ EAg1 ; 11 0 0 > > > Gc33 Ec11 E > n21 ¼ q‘2 x2 ; n22 ¼ q‘2 c x2 ; n23 ¼ jGAg2 ; n24 ¼ EAg2 ; > > 22 0 > 0 > > > : n31 ¼ G2 2 ; n32 ¼ Eðc222 c112 Þ ; n33 ¼ c22 c11 : c33 q‘ x q‘ c x 0
33
ð47Þ
ð48Þ
0
Then, substituting the relations (18), (13), (25), (26), (29), (38) and (46) into the nonlinear equation (39), we obtain the governing differential equation of motion in terms of dimensionless displacements pi(s, t) (i = 1, 2, . . . , 6) as follows:
470
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
€ p1 ¼ fðn01 n02 Þðp01 cos p4 sin p5 þ p02 sin p4 ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 sin p5 0
þ n01 p01 n02 cos p4 sin p5 g þ g0 fx ; € p2 ¼ fðn01 n02 Þðp1 0 sin p4 sin p5 þ þ
n02 ðp02
p02
ð49Þ
cos p4 þ ð1 þ
p03 Þ sin p4
cos p5 Þ cos p4
0
þ sin p4 Þg þ g0 fy ;
ð50Þ
€ p3 ¼ fðn01 n02 Þðp01 cos p4 sin p5 þ p02 sin p4 þ n01 ð1 þ p03 Þ n02 cos p4 cos p5 g0 þ g0 fz ;
ð1 þ
p03 Þ cos p4
cos p5 Þ cos p4 cos p5 ð51Þ
:
ðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ ðp_ 6 p_ 5 sin p4 Þðp_ 4 sin p6 p_ 5 cos p4 cos p6 Þ 0
¼ n11 ðp04 cos p6 þ p05 sin p6 cos p4 Þ n12 ðp06 p05 sin p4 Þðp04 sin p6 p05 cos p6 cos p4 Þ ððn13 n14 Þðp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ þ n14 Þ fp01 ð sin p6 cos p5 þ cos p6 sin p4 sin p5 Þ þ p02 cos p6 cos p4 þ ð1 þ p03 Þðsin p6 sin p5 þ cos p6 sin p4 cos p5 Þg þ g1 lx ðcos p6 cos p5 þ sin p6 sin p4 sin p5 Þ þ g1 ly sin p6 cos p4 þ g1 lz ð cos p6 sin p5 þ sin p6 sin p4 cos p5 Þ; : ðp_ 4 sin p6 þ p_ 5 cos p6 cos p4 Þ ðp_ 6 p_ 5 sin p4 Þðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ
ð52Þ
0
¼ n21 ðp05 cos p6 cos p4 p04 sin p6 Þ n22 ðp06 p05 sin p4 Þðp04 cos p6 þ p05 sin p6 cos p4 Þ þ ððn23 n24 Þðp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ þ n24 Þ fp01 ðcos p6 cos p5 þ sin p6 sin p4 sin p5 Þ þ p02 sin p6 cos p4 þ ð1 þ p03 Þð cos p6 sin p5 þ sin p6 sin p4 cos p5 Þg þ g2 lx ð sin p6 cos p5 þ cos p6 sin p4 sin p5 Þ þ g2 ly cos p6 cos p4 þ g2 lz ðsin p6 sin p5 þ cos p6 sin p4 cos p5 Þ; ðp_ 6 p_ 5 sin p4 Þ: þ n33 ðp_ 4 sin p6 þ p_ 5 cos p4 cos p6 Þðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ
ð53Þ
0
¼ n31 ðp06 p05 sin p4 Þ þ n32 ðp04 sin p6 þ p05 cos p6 cos p4 Þðp04 cos p6 þ p05 sin p6 cos p4 Þ þ g3 ðlx cos p4 sin p5 ly sin p4 þ lz cos p4 cos p5 Þ;
ð54Þ
where the symbols (.) and ( 0 ) denote differentiation with respect to dimensionless time parameter s and dimensionless length parameter r, respectively. The boundary conditions at a fixed end are p1 ¼ p2 ¼ p3 ¼ p4 ¼ p5 ¼ p6 ¼ 0:
ð55Þ
At a free end, the force boundary conditions (44) can be rewritten as am fðn01 n02 Þðp01 cos p4 sin p5 þ p02 sin p4 ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 sin p5 g þ n01 p01 n02 cos p4 sin p5 g þ 0 fxb ðsÞ ¼ 0; ‘ am fðn01 n02 Þðp01 sin p4 sin p5 þ p02 cos p4 þ ð1 þ p03 Þ sin p4 cos p5 Þ cos p4 g þ n02 ðp02 þ sin p4 Þg þ 0 fyb ðsÞ ¼ 0; ‘ am fðn01 n02 Þðp01 cos p4 sin p5 þ p02 sin p4 ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 cos p5 g þ n01 ð1 þ p03 Þ n02 cos p4 cos p5 g þ 0 fzb ðsÞ ¼ 0; ‘ while the torque boundary conditions (45) can be rewritten as
ð56Þ
ð57Þ
ð58Þ
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
am n11 ðp04 cos p6 þ p05 sin p6 cos p4 Þ þ
471
g1 b l ðcos p6 cos p5 þ sin p6 sin p4 sin p5 Þ ‘ x
g1 b g l sin p6 cos p4 þ 1 lbz ð cos p6 sin p5 þ sin p6 sin p4 cos p5 Þ ¼ 0; ‘ y ‘ g 0 0 am n21 ðp5 cos p6 cos p4 p4 sin p6 Þ þ 2 lbx ð sin p6 cos p5 þ cos p6 sin p4 sin p5 Þ ‘ g g þ 2 lby cos p6 cos p4 þ 2 lbz ðsin p6 sin p5 þ cos p6 sin p4 cos p5 Þ ¼ 0; ‘ ‘ g3 b 0 0 am n31 ðp6 p5 sin p4 Þ þ ðlx cos p4 sin p5 lby sin p4 þ lbz cos p4 cos p5 Þ ¼ 0: ‘ þ
ð59Þ
ð60Þ ð61Þ
5. Femlab implementation Femlab is an open flexible system used to model all types of scientific and engineering problems based on partial differential equations (Femlab homepage). This gives us freedoms to implement numerical solution procedure of the nonlinear partial differential equations of motion with the corresponding boundary conditions and initial conditions. We introduce the solution procedure of the nonlinear partial differential equations (49)–(54) with the boundary conditions (55) and/or (56)–(61). Let 2 b 3 2 3 2 3 2 3 2 3 fx ðsÞ p_ 1 b1 g1 p1 6 b 7 6p 7 6 p_ 7 6b 7 6g 7 6 f ðsÞ 7 6 27 6 27 6 27 6 27 7 6 y 6 7 6 7 6 7 6 7 6 6 p3 7 6 p_ 3 7 6 b3 7 6 g3 7 g0 6 fzb ðsÞ 7 7 b 7 7 7 7 6 6 6 6 ð62Þ p ¼ 6 7; q ¼ p_ ¼ 6 7; b ¼ 6 7; g ¼ 6 7; g ¼ 6 b 7: ‘ 6 l1 ðsÞ 7 6 p4 7 6 p_ 4 7 6 b4 7 6 g4 7 7 6 6 7 6 7 6 7 6 7 6 lb ðsÞ 7 4 p5 5 4 p_ 5 5 4 b5 5 4 g5 5 5 4 2 b p6 b6 g6 p_ 6 l3 ðsÞ where 8 b1 ¼ ðn01 n02 Þðþp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 sin p5 > > > > > n01 p01 þ n02 cos p4 sin p5 ; > > > > > b2 ¼ ðn02 n01 Þðp01 sin p4 sin p5 þ p02 cos p4 þ ð1 þ p03 Þ sin p4 cos p5 Þ cos p4 > > > < n02 ðp02 þ sin p4 Þ; > b3 ¼ ðn01 n02 Þðp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 cos p5 > > > > > n01 ð1 þ p03 Þ þ n02 cos p4 cos p5 ; > > > > > b4 ¼ n11 ðp04 cos p6 þ p05 sin p6 cos p4 Þ; > > : b5 ¼ n21 ðp04 sin p6 p05 cos p6 cos p4 Þ; b6 ¼ n31 ðp06 þ p05 sin p4 Þ: Taking notice of 8 ðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ: ðp_ 6 p_ 5 sin p4 Þðp_ 4 sin p6 p_ 5 cos p4 cos p6 Þ > > > > > ¼€ p4 cos p6 þ € p5 sin p6 cos p4 2p_ 4 p_ 6 sin p6 > > > > > _ _ þ2 p cos p4 cos p6 p_ 25 sin p4 cos p4 cos p6 ; p > 5 6 > > > : > > < ðp_ 4 sin p6 þ p_ 5 cos p6 cos p4 Þ ðp_ 6 p_ 5 sin p4 Þðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ ¼ € p4 sin p6 þ € p5 cos p6 cos p4 2p_ 4 p_ 6 cos p6 > > > 2p_ 5 p_ 6 sin p6 cos p4 þ p_ 25 sin p4 cos p4 sin p6 ; > > > > : > > ðp_ 6 p_ 5 sin p4 Þ þ n33 ðp_ 4 sin p6 þ p_ 5 cos p4 cos p6 Þðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ > > > > > p6 p€5 sin p4 Þ p_ 5 p_ 4 cos p4 > ¼ ð€ > : þn33 ðp_ 24 sin p6 cos p6 þ p_ 4 p_ 5 cos 2p6 cos p4 þ p_ 25 sin p6 cos p6 cos2 p4 Þ; the nonlinear partial differential equations (49)–(54) can be rewritten as
ð63Þ
ð64Þ
472
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
p_
þ
M a q_
0
b0
q ¼ ; g
where Ma ¼
I 33 Mm
;
ð65Þ 2
cos p6
sin p6 cos p4
0
cos p6 cos p4 sin p4
6 M m ¼ 4 sin p6
0
3
7 0 5; 1
ð66Þ
and 8 g1 > > > > > g4 > > > > > > > > > > > > > > > > > > > > > < g5 > > > > > > > > > > > > > > > > > > > > > > > > g6 > > :
¼ g0 fx ; ¼
g 2 ¼ g0 f y ;
n12 ðp06
p05
g3 ¼ g0 fz ;
sin p4 Þðp04 sin p6 p05 cos p6 cos p4 Þ
ððn13 n14 Þðp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ þ n14 Þ fp01 ð sin p6 cos p5 þ cos p6 sin p4 sin p5 Þ þ p02 cos p6 cos p4 þð1 þ p03 Þðsin p6 sin p5 þ cos p6 sin p4 cos p5 Þg þ g1 l1 þ2p_ 4 p_ 6 sin p6 2p_ 5 p_ 6 cos p4 cos p6 þ p_ 25 sin p4 cos p4 cos p6 ; ¼ n22 ðp06 p05 sin p4 Þðp04 cos p6 þ p05 sin p6 cos p4 Þ n24 Þðp01
p02
ð67Þ p03 Þ cos p4
þððn23 cos p4 sin p5 sin p4 þ ð1 þ cos p5 Þ þ n24 Þ fp01 ðcos p6 cos p5 þ sin p6 sin p4 sin p5 Þ þ p02 sin p6 cos p4 þð1 þ p03 Þð cos p6 sin p5 þ sin p6 sin p4 cos p5 Þg þ g2 l2 þ2p_ 4 p_ 6 cos p6 þ 2p_ 5 p_ 6 sin p6 cos p4 p_ 25 sin p4 cos p4 sin p6 ; ¼ n32 ðp04 sin p6 þ p05 cos p6 cos p4 Þðp04 cos p6 þ p05 sin p6 cos p4 Þ þ g3 l3 þ p_ 5 p_ 4 cos p4 þn33 ðp_ 24 sin p6 cos p6 p_ 4 p_ 5 cos 2p6 cos p4 p_ 25 sin p6 cos p6 cos2 p4 Þ:
Using the general form for nonlinear PDEs with 12 dependent variables, Eq. (65) can be implemented in Femlab in terms of the Matlab code (for details, see Comsol, 2004). The Dirichlet boundary conditions should be adopted at a fixed end, i.e., p = 0. At a free end, the Neumann boundary conditions should be adopted, i.e., amb = gb, where b and gb are defined by (62). For the formulated nonlinear model of a Cosserat rod, the numerical analysis using the Femlab programming language at the Matlab command line consists of the following steps: 1. 2. 3. 4.
Formulation of the governing equations; Discretization of the equations; Solution of the equations; Interpretation of the results.
Given the one-dimensional geometry data, an initial finite element mesh is easily generated by meshinit. The meshextend command is used for discretization of the PDE problem and can be modified to improve accuracy. The geometry, nonlinear PDEs and boundary conditions are defined by a set of fields in the Matlab codes. For solving purposes Femlab contains specific solvers (like static, dynamic, linear, nonlinear solvers) for specific PDE problems. Here, we use femtime to solve the dynamic problem of the system. 6. Simulation results and discussion for a simple cantilever A cantilever (fixed at the left hand and free at the right hand), as shown in Fig. 2, is now presented as a simple example to demonstrate the numerical solution procedure of the proposed Cosserat rod model. Numerical calculations based on (65) are carried out for the uniform horizontal cantilever of length ‘ = 1.0 m, of constant cross-section with width b = 0.06 m and thickness a = 0.04 m. The mass density, the Young’s modulus and the Poisson’s ratio are assumed to be q = 2.730 · 103 kg/m3, E = 7.1 · 1010 Pa and m = 0.32, respectively. The shear correction factor for rectangular cross-section is taken to be j = 0.833. The reference frequency is
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
473
chosen to be the lowest natural frequency of the linearized system of a cantilever without considering the effect of shearing deformations, i.e., sffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffi 1:8752 Ec22 1:8752 Ea2 x0 ¼ ¼ ¼ 207:0236 rad s1 : qA 12q ‘2 ‘2 Based on the derived nonlinear system (65), numerical simulations are performed to investigate the dynamic responses of the cantilever under harmonic excitations. The differential equations of motion are full coupled by the nonlinear terms and could exhibit internal resonance introduced by the nonlinearities. They also exhibit external resonances when the external excitation is periodic and the frequency of a component of its Fourier series is near one of the natural frequencies of the system, or near a multiple of a natural frequencies. The detailed analysis of complex dynamic behavior, such as bifurcation (Vu-Quoc and Li, 1995) and chaos, of the system is not the main focus of this paper. We only present here the simulation results of the dynamic responses of the system under external harmonic excitations. The dynamical responses of the free end of the cantilever (with zero initial displacements and velocities) under external distributing loads fx ðtÞ ¼ 2 sinðpsÞ sinð8x0 tÞ kN m1 ;
f y ðtÞ ¼ 2 sinðpsÞ sinð8x0 tÞ kN m1 ;
are obtained using the Femlab/Matlab, and shown in Fig. 4. It is observed that the transverse displacements x(1, t) and y(1, t) consist of more than one harmonic motions with different frequencies. The longitudinal displacement z(1, t) and the torsional oscillation /(1, t) are obviously not periodic motions, which is not surprising, since various dynamic phenomena are expected for such a strongly nonlinear and fully coupling system. Moreover, it can be observed that the longitudinal displacement z(1, t) is a higher order small term comparing with other variables like x(1, t) and y(1, t). The transverse displacement time histories x(s, t) and y(s, t) of the nonlinear coupling model under external harmonic excitations are shown in Figs. 5 and 6, respectively. Dividing the cantilever into 10 Cosserat Rod Elements (CRE) of equal length, see Cao et al. (2006) for the formulation of CREs, one can establish the nonlinear ordinary differential equations of motion for solving the free displacements. Numerical simulations for the responses of the model under the same external harmonic −3
−3
x 10
θ(1,t) [rad]
x(1,t) [m]
2 0 −2
0
0.02
0.04
0.06
5 0 −5
−3
0 −2
0
0.02
0.04
0.06
5
φ(1,t) [rad]
z(1,t) [m]
0 −4
0
0.04
0.06
0.02
0.04
0.06
0.04
0.06
x 10
−5
0 −5
x 10
−2 −6
0.02
0
−6
2
0 −3
x 10
ψ(1,t) [rad]
y(1,t) [m]
2
x 10
0.02
0.04
t [sec]
0.06
2
x 10
0 −2 −4
0
0.02
t [sec]
Fig. 4. Dynamical responses of the cantilever at the free end under harmonic excitations.
474
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477 −3
x 10 1.5 1
x(s,t) [m]
0.5 0 −0.5 −1 −1.5 0.06 0.04
t [sec]
0.02 0
0.2
0
0.4
0.6
0.8
1
s [m]
Fig. 5. Displacement time histories of the cantilever x(s, t) [m] under harmonic excitations.
−3 x 10 1
y(s,t) [m]
0.5 0 −0.5 −1 −1.5 0.06 1
0.04
t [sec]
0.8 0.6
0.02
0.4 0
0.2
s [m]
0
Fig. 6. Displacement time histories of the cantilever y(s, t) [m] under harmonic excitations.
excitations are performed with Matlab. The corresponding dynamical responses of the CRE model are depicted in Figs. 7–9 for a comparison of the results in Figs. 4–6 obtained in this paper. 7. Conclusions The nonlinear partial differential equations of motion for slender elastic rods have been formulated using the Cosserat theory. Unlike other formulations presented in the literature, the nonlinear dynamic model formulated here are explicitly in terms of the displacements and angular variables. Flexure along two principal axes, extension, shear and torsion are considered, and care is taken to account for all the nonlinearities include contributions from both the curvature expression and inertia terms.
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
θ(1,t) [rad]
5
0
0.02
x 10
0.04
−3
5
0 −2
2
0
0.02
x 10
0.04
−6
2
0 −2 −4 −6
0
0.02
0.04
0 x 10
0.02
0.04
0.06
0.02
0.04
0.06
0.04
0.06
−3
0 x 10
−5
0 −2 −4
0.06
−3
0 −5
0.06
x 10
0 −5
0.06
ψ(1,t) [rad]
2
y(1,t) [m]
−3
0 −2
z(1,t) [m]
x 10
φ(1,t) [rad]
x(1,t) [m]
2
475
0
0.02
t [sec]
t [sec]
Fig. 7. Dynamical responses of the cantilever at the free end under harmonic excitations: Cosserat Rod Element Approach.
−3
x 10 1.5
x(s,t) [m]
1 0.5 0 −0.5 −1 −1.5 0.06 0.04
t [sec]
0.02 0
0
0.2
0.4
0.6
0.8
1
s [m]
Fig. 8. Displacement time histories of the cantilever x(s, t) [m] under harmonic excitations: Cosserat Rod Element Approach.
A numerical solution procedure in terms of the Femlab/Matlab interfaces has been presented for numerically solving the formulated nonlinear dynamical model with the corresponding boundary and initial conditions. A simple cantilever under external harmonic excitations has been presented to demonstrate the use of the formulation developed, and to illustrate the numerical solution procedure. Based on the case study of applying the proposed modelling strategy for nonlinear dynamics of elastic rods, the following conclusions have been drawn.
476
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
−3
x 10 1
y(s,t) [m]
0.5 0 −0.5 −1 −1.5 0.06 1
0.04
t [sec]
0.8 0.6
0.02
0.4 0
0.2 0
s [m]
Fig. 9. Displacement time histories of the cantilever y(s, t) [m] under harmonic excitations: Cosserat Rod Element Approach.
1. The explicit nature of the construction permits us to directly express the nonlinear dynamic model of the elastic rod and the corresponding boundary conditions in terms of the displacements and angular variables, thus the dynamical analysis of elastic rods can be carried out rather simply. 2. The governing partial differential equations of motion can be easily expanded to contain nonlinearities up to order three to render them amenable to the study of moderately large amplitude flexural-torsional-shearable oscillations by perturbation techniques. Moreover, the approach on dimension reduction of dynamical systems, such as the nonlinear Galerkin method (Rega and Troger, 2005) and the Karhunen–Loeve method (which in the mathematical literature is also called Proper Orthogonal Decomposition (POD) method (Georgiou, 2005)), could be used to develop reduced-order dynamic models of the coupled nonlinear partial differential systems derived here. Therefore, the full coupled nonlinear model formulated here could be very useful for the detailed analysis of complex dynamic behavior, such as bifurcation and chaos, of the slender rods. 3. The nonlinear dynamic model could be numerically solved in terms of the Femlab/Matlab interfaces. The Femlab environment offers a powerful simulation tool for analyzing the responses of nonlinear dynamic models (of slender material elements such as those found in an MEMS device) presented here.
Acknowledgements The authors are grateful to the EPSRC (Computational Engineering Mathematics Programme) for financial support in this study. References Antman, S.S., 1991. Nonlinear Problems of Elasticity. In: Applied Mathematical Sciences. Springer-Verlag, New York. Antman, S.S., Marlow, R.S., Vlahacos, C.P., 1998. The complicated dynamics of heavy rigid bodies attached to deformable rods. Quarterly of Applied Mathematics 86, 431–460. Bishop, R.E.D., Price, W.G., 1977. Coupled bending and twisting of a Timoshenko beam. Journal of Sound and Vibration 50, 469–477. Bolotin, V.V., 1964. The Dynamic Stability of Elastic Systems. Holden-Day, San Francisco. Bruch Jr., J.C., Mitchell, T.P., 1987. Vibrations of mass-loaded clamped-free Timoshenko beam. Mitchell Journal of Sound and Vibration 114 (2), 341–345. Cao, D.Q., Liu, D., Wang, C.H.-T., 2005. Nonlinear dynamic modelling for MEMS components via the Cosserat rod element approach. Journal of Micromechanics and Microengineering 15, 1334–1343.
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
477
Cao, D.Q., Liu, D., Wang, C.H.-T., 2006. Three-dimensional nonlinear dynamics of slender structures: Cosserat rod element approach. International Journal of Solids and Structures 43, 760–783. Cartmell, M.P., 1990. The equations of motion for a parametrically excited cantilever beam. Journal of Sound and Vibration 143 (3), 395– 406. Comsol, AB., 2004. Femlab 3.0: Matlab Interface Guide. COMSOL Ltd, Oxford. Crespo da Silva, M.R.M., 1988a. Nonlinear flexural-flexural-torsional-extensional dynamics of beams – I. Formulation. International Journal of Solids and Structures 24 (12), 1225–1234. Crespo da Silva, M.R.M., 1988b. Nonlinear flexural-flexural-torsional-extensional dynamics of beams – II. Response analysis. International Journal of Solids and Structures 24 (12), 1235–1242. Crespo da Silva, M.R.M., Glynn, C.C., 1978a. Nonlinear flexural-flexural-torsional dynamics of inextensional beams, I: Equations of motion. Journal of Structural Mechanics 6 (4), 437–448. Crespo da Silva, M.R.M., Glynn, C.C., 1978b. Nonlinear flexural-flexural-torsional dynamics of inextensional beams, II: Forced motions. Journal of Structural Mechanics 6 (4), 449–461. Cull, S.J., Tucker, R.W., Tung, R.S., Hartley, D.H., 2000. On parametrically excited flexural motion of an extensible and shearable rod with a heavy attachment. Technische Mechanik 20, 147–158. Esmailzadeh, E., Jalili, N., 1998. Parametric response of cantilever Timoshenko beams with tip mass under harmonic support motion. International Journal of Non-Linear Mechanics 33, 765–781. Femlab homepage:. Forehand, D.I.M., Cartmell, M.P., 2001. On the derivation of the equations of motion for a parametrically excited cantilever beam. Journal of Sound and Vibration 245 (1), 165–177. Georgiou, I., 2005. Advanced proper orthobonal decomposition tools: using reduced order models to identify normal modes of vibration and slow invariant manifolds in the dynamics of planar nonlinear rods. Nonlinear Dynamics 41, 69–110. Grant, D.A., 1978. The effect of rotary inertia and shear deformation on the frequency and normal mode equations of uniform beams carrying a concentrated mass. Journal of Sound and Vibration 57 (3), 357–365. Gratus, J., Tucker, R.W., 2003. The dynamics of Cosserat nets. Journal of Applied Mathematics (4), 187–226. Huang, T.C., 1961. The Effect of rotatory inertia and of shear deformation on the frequency and normal mode equations of uniform beams with simple end conditions. Transactions ASME Journal of Applied Mechanics 28, 579–584. Kar, R.C., Dwivedy, S.K., 1999. Nonlinear dynamics of a slender beam carrying a lumped mass with principal parametric and internal resonances. International Journal of Non-Linear Mechanics 34, 515–529. Nayfeh, A.H., Mook, D.T., 1979. Nonlinear Oscillations. Wiley-Interscience, New York. Rega, G., Troger, H., 2005. Dimension reduction of dynamical systems: methods, models, applications. Nonlinear Dynamics 41, 1–15. Rubin, M.B., 2000. Cosserat Theories: Shells, Rods and Points. Kluwer Academic Publishers, New York. Rubin, M.B., 2001. Numerical solution procedures for nonlinear elastic rods using the theory of a Cosserat point. International Journal of Solids and Structures 38, 4395–4437. Rubin, M.B., Tufekci, E., 2005. Three-dimensional free vibration of a circular arch using the theory of a Cosserat point. Journal of Sound and Vibration 286, 799–816. Saito, H., Koizumi, N., 1982. Parametric vibrations of a horizontal beam with a concentrated mass at one end. International Journal of Mechanical Sciences 24, 755–761. Simo, J.C., Vu-Quoc, L., 1988. On the dynamics in space of rods undergoing large motions – a geometrically exact approach. Computer Methods in Applied Mechanics and Engineering 66, 125–161. Simo, J.C., Vu-Quoc, L., 1991. A geometrically-exact beam model incorporating shear and torsion warping deformation. International Journal of Solids and Structures 27, 371–393. Tucker, R.W., Wang, C., 1999. An integrated model for drill-string dynamics. Journal of Sound and Vibration 224, 123–165. Vu-Quoc, L., Deng, H., 1995. Galerkin projection for geometrically-exact sandwich beams allowing for ply drop-off. ASME Journal of Applied Mechanics 62, 479–488. Vu-Quoc, L., Deng, H., 1997. Dynamics of geometrically-exact sandwich beams: computational aspects. Computer Methods in Applied Mechanics and Engineering 146, 135–172. Vu-Quoc, L., Ebcioglu, I.K., 1995. Dynamic formulation for geometrically-exact sandwich beams and 1-D plates. ASME Journal of Applied Mechanics 62, 756–763. Vu-Quoc, L., Li, S., 1995. Dynamics of sliding geometrically-exact beams: large angle maneuver and parametric resonance. Computer Methods in Applied Mechanics and Engineering 120, 65–118.
International Journal of Solids and Structures 45 (2008) 460–477 www.elsevier.com/locate/ijsolstr
Nonlinear dynamics of elastic rods using the Cosserat theory: Modelling and simulation D.Q. Cao a
a,*,1
, Robin W. Tucker
b
School of Astronautics, Harbin Institute of Technology, P.O. Box 137, Harbin 150001, China b Department of Physics, Lancaster University, Lancaster LA1 4YB, UK Received 20 March 2006; received in revised form 9 July 2007 Available online 6 September 2007
Abstract The method of Cosserat dynamics is employed to explore the nonplanar nonlinear dynamics of elastic rods. The rod, which is assumed to undergo flexure about two principal axes, extension, shear and torsion, are described by a general geometrically exact theory. Based on the Cosserat theory, a set of governing partial differential equations of motion with arbitrary boundary conditions is formulated in terms of the displacements and angular variables, thus the dynamical analysis of elastic rods can be carried out rather simply. The case of doubly symmetric cross-section of the rod is considered and the Kirchoff constitutive relations are adopted to provide an adequate description of elastic properties in terms of a few elastic moduli. A cantilever is given as a simple example to demonstrate the use of the formulation developed. The nonlinear dynamic model with the corresponding boundary and initial conditions are numerically solved using the Femlab/Matlab software packages. The corresponding nonlinear dynamical responses of the cantilever under external harmonic excitations are presented through numerical simulations. 2007 Elsevier Ltd. All rights reserved. Keywords: Cosserat model; Nonlinear dynamics; Slender rod; Modelling and simulation; Femlab; Matlab
1. Introduction Nonlinear dynamic analysis of slender elastic structures under external forces and torques and parametric excitations remains an active area of study. Such a study can find application in accelerating missiles and space crafts, components of high-speed machinery, manipulator arm, microelectronic mechanical structures (MEMS), components of bridges (such as towers and cables) and other structural elements. Considerable attention has been devoted to the study of nonlinear dynamics of rods or beams subject to both external and parametric excitations (see, for example, Nayfeh and D.T.Mook, 1979; Saito and Koizumi, 1982; Kar and Dwivedy, 1999 and the references cited therein). While attention so far has mainly been devoted to the *
1
Corresponding author. Tel.: +86 (0) 451 86414479; fax: +86 (0) 451 86402822. E-mail addresses: [email protected] (D.Q. Cao), [email protected] (R.W. Tucker). On leave from the Department of Physics, Lancaster University, Lancaster LA1 4YB, UK.
0020-7683/$ - see front matter 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2007.08.016
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
461
study of planar, nonlinear dynamic analysis of beams, research has been done concerning nonplanar, nonlinear motions of beams. Bolotin (1964) addressed such motions in beams, but restricted himself to consideration of nonlinear inertia terms and stability of the planar response. Crespo da Silva and Glynn (1978a) formulated the equations of motion describing the nonplanar, nonlinear dynamics of an inextensional beam. The nonplanar, nonlinear forced oscillations of a cantilever are then analyzed in Crespo da Silva and Glynn (1978b) using the perturbation method. Cartmell (1990) and Forehand and Cartmell (2001) derived the nonlinear equation of motion for the in-plane and out-of-plane forced vibration of cantilever beams with a lumped mass. The studies mentioned above are restricted to systems with no extension of the beam’s neutral axis and no warping or shear deformation. Crespo da Silva (1988a,b) investigated the problem of nonlinear dynamics of the nonplanar flexural-torsional-extensional beams, but the effects of rotary inertia and shear deformation were neglected. Thus, cross-sectional dimensions of the beam were assumed to be small enough in comparison to the beam length. However, the shear deformation may be of considerable importance and can not be negligible for studying the vibration of high frequencies when a comparative short rod is investigated. In such a case, the effect of shear deformation should be taken into account for. For the planar problem, although such effects can be included by using the Timoshenko beam theory, most of the studies are limited to the determination of natural frequencies and eigenfunctions (Huang, 1961; Bishop and Price, 1977; Grant, 1978; Bruch and Mitchell, 1987). For the three-dimensional problem, we refer to Simo and Vu-Quoc (1988, 1991) for a geometrically exact rod model incorporating shear and torsion-warping deformation. Further studies on the dynamic formulation of sandwich beams have been presented in Vu-Quoc and Deng (1995) and Vu-Quoc and Ebcioglu (1995) based on the geometrically-exact description of the kinematics of deformation. Moreover, Esmailzadeh and Jalili (1998) investigated the parametric response of cantilever Timoshenko beams with lumped mass, but restricted themselves to consideration of nonlinear inertia terms. Based on the geometrically-exact model of sliding beams, parametric resonance has been presented in Vu-Quoc and Li (1995), where the beams can undergo large deformation with shear deformation accounted for. In the case involving the full dynamic response, the strong nonlinearity and fully coupling introduce a challenge for solving the partial differential equations of motion of rods. We refer to Rubin (2001) for a formulation of a numerical solution procedure for three-dimensional dynamic analysis of rods by modelling the rod as a set of connected Cosserat points, also Rubin and Tufekci (2005) and Rubin (2000) for three-dimensional dynamics of a circular arch and shells using the theory of a Cosserat point, respectively. A Galerkin projection has been applied to discretize the governing partial differential equations of sliding beams in Vu-Quoc and Li (1995) and sandwich beams in Vu-Quoc and Deng (1997). Recently, a new modelling strategy has been proposed in Cao et al. (2006) to discretize the rod and to derive the ordinary differential equations of motion with third order nonlinear generic nodal displacements. This modelling strategy has been successfully used to investigate the nonlinear dynamics of typical MEMS device that comprises a resonator mass supported by four flexible beams (Cao et al., 2005). In this paper we explore the nonplanar, nonlinear dynamics of an extensional shearable rod by using the simple Cosserat model. The method of Cosserat dynamics for elastic structures is employed since it can accommodate to a good approximation the nonlinear behavior of complex elastic structures composed of materials with different constitutive properties, variable geometry and damping characteristics (Antman, 1991; Antman et al., 1998; Tucker and Wang, 1999; Cull et al., 2000; Gratus and Tucker, 2003). With arbitrary boundary conditions, the Cosserat theory is used to formulate a set of governing partial differential equations of motion in terms of the displacements and angular variables, describing the nonplanar, nonlinear dynamics of an extensional rod. Bending about two principal axes, extension, shear and torsion are considered, and care is taken into account for all the nonlinear terms in the resulting equations. As an example, a simple cantilever is given to demonstrate the use of the formulation developed. The nonlinear dynamic model with the corresponding boundary and initial conditions are numerically solved using the commercially available software packages Femlab and Matlab. Corresponding nonlinear dynamical responses of the cantilever under external harmonic excitations are presented and discussed through numerical simulations. The following conventions and nomenclature will be used through out this paper. Vectors, which are elements of Euclidean 3-space R3 , are denoted by lowercase, bold-face symbols, e.g., u, v; vector-valued functions are denoted by lowercase, italic, bold-face symbols, e.g., u, v; tensors are denoted by upper-case, bold-face symbols, e.g., I, J; matrices are denoted by upper-case, italic, bold-face symbols, e.g., M, K. The symbols
462
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
ot and os denote differentiation with respect to time t and arc-length parameter s, respectively. The symbols ð : Þ and ( 0 ) denote differentiation with respect to dimensionless time parameter s and dimensionless length parameter r, respectively. 2. Background on the simple Cosserat model In this section the basic concepts of the Cosserat theory for an elastic rod of unstressed length ‘ are summarized. The three vectors {e1, e2, e3} are assumed to form a fixed right-handed orthogonal basis. Elements of the rod are labelled in terms of the Lagrangian coordinate 06 s 6 ‘ at time t. The motion of an elastic segment may be described in terms of the motion in space of a vector r that locates the line of centroids (shown dashdotted) of the cross-sections as shown in Fig. 1. Specifying a unit vector d3 (which may be identified with the normal to the cross-section) at each point along this line enables the state of shear to be related to the angle between this vector and the tangent osr to the centroid space-curve. Specifying a second vector d2 orthogonal to the first vector (thereby placing it in the plane of the cross-section) can be used to encode the state of bending and twist along the element. In such a way, the current orientation of each cross-section at s 2 [0, ‘] is defined by specifying the orientation of a moving basis d1, d2, d3 = d1 · d2, see Simo and Vu-Quoc (1988) and Vu-Quoc and Deng (1995) for details. Elastic deformations about the line of centroids are then coded into the rates of change of r and the triad d1, d2, d3 of the cross-section at s. Thus a time dependent field of two mutually orthogonal unit vectors along the segment provides three continuous dynamical degrees of freedom that, together with the three continuous degrees of freedom describing the centroid space-curve relative to some arbitrary origin in space (with fixed inertial frame e1, e2, e3), define a simple Cosserat rod model. Supplemented with appropriate constitutive relations and boundary conditions such a model can fully accommodate the modes of vibration that are traditionally associated with the motion of slender rods in the engineering literature: namely axial motion along its length, torsional or rotational motion and transverse or lateral motion. The tangent to the space-curve, v = v(s, t), is given by vðs; tÞ ¼ os rðs; tÞ;
ð1Þ
the rotation of the directions along the space curve and the local angular velocity vector of the director triad are measured by u = u(s, t), w = w(s, t) according to os di ðs; tÞ ¼ uðs; tÞ di ðs; tÞ; ð2Þ ot di ðs; tÞ ¼ wðs; tÞ di ðs; tÞ; respectively. It follows from the first equation of (2) that
s=b
s
r
d3
d2
e3 r
s d1
o
e2
e1 Fig. 1. A simple Cosserat model.
s=a
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477 3 3 3 X X X ðdi os di Þ ¼ ðdi ðu di ÞÞ ¼ ðuðdi di Þ di ðdi uÞÞ ¼ 2u: i¼1
i¼1
463
ð3Þ
i¼1
Similarly, we have 3 3 3 X X X ðdi ot di Þ ¼ ðdi ðw di ÞÞ ¼ ðwðdi di Þ di ðdi wÞÞ ¼ 2w: i¼1
i¼1
ð4Þ
i¼1
Since the basis {d1, d2, d3} is natural for the intrinsic description of deformation, we decompose relevant vector valued functions with respect to it: 8 3 P > > > vðs; tÞ ¼ vi ðs; tÞdi ðs; tÞ; > > > i¼1 > > < 3 P ð5Þ uðs; tÞ ¼ ui ðs; tÞdi ðs; tÞ; > i¼1 > > > > 3 > P > > wi ðs; tÞdi ðs; tÞ: : wðs; tÞ ¼ i¼1
The contact forces n(s, t) and contact torques m(s, t) are related to the extension and shear strains v(s, t), and flexure and torsion strains u(s, t) by constitutive relations, respectively. The angular momentum h(s, t) is related to the rotary inertia and the director angular velocity w(s, t). The dynamical evolution of the rod with density, q(s), and cross-section area, A(s) is governed by the Newton’s dynamical laws: qðsÞAðsÞott r ¼ os nðs; tÞ þ fðs; tÞ; ð6Þ ot hðs; tÞ ¼ os mðs; tÞ þ vðs; tÞ nðs; tÞ þ lðs; tÞ; where f(s, t) and l(s, t) denote external force and torque densities, respectively, the contact force vector, contact torque vector, and the angular momentum vector can be written as 8 3 P > > > ni ðs; tÞdi ðs; tÞ; nðs; tÞ ¼ < i¼1 ð7Þ 3 > P > > mi ðs; tÞdi ðs; tÞ; : mðs; tÞ ¼ i¼1
and hðs; tÞ ¼
3 X
hi ðs; tÞdi ðs; tÞ;
ð8Þ
i¼1
respectively. For more details about Eqs. (6)–(8), we refer to the governing equations (i)–(iii) derived in terms of the noncommutative Lie group of proper orthogonal transformations in Simo and Vu-Quoc (1988), see also in Vu-Quoc and Deng (1995) and Vu-Quoc and Ebcioglu (1995). It will be assumed that the Young’s modulus E, the shear modulus G, and the material density q along the rod are functions of s only, and the mass center coincides with the area centroid of the cross-section at s. In this case, the simplest constitutive model for the linear elastic material is based on the Kirchhoff constitutive relations with shear deformation, see the discussions in Simo and Vu-Quoc (1988, 1991) for details of constitutive relations. Such a model provide an adequate description of elastic properties in terms of a few elastic moduli. The presence of arbitrary rotations relating the local director frame to the global inertial frame renders the equations of motion inherently nonlinear. The contact forces, contact torques and the angular momentum are then given as n ¼ Kðv d3 Þ;
m ¼ Ju;
h ¼ Iw;
where, for a symmetric cross-section of a rod, the tensors K, J and I are described as
ð9Þ
464
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
8 3 P > > > Kðs; tÞ ¼ K ii ðs; tÞðdi ðs; tÞ di ðs; tÞÞ; > > > i¼1 > > < 3 P Jðs; tÞ ¼ J ii ðs; tÞðdi ðs; tÞ di ðs; tÞÞ; > i¼1 > > > > 3 > P > > Iðs; tÞ ¼ I ii ðs; tÞðdi ðs; tÞ di ðs; tÞÞ: :
ð10Þ
i¼1
Here the symbol is used to denote the tensorial product u v of two vectors u and v that assigns to each vector w the vector (u v)w = (v Æ w)u. The corresponding terms in (10) can be written as 8 K 11 ¼ K 22 ¼ jðsÞGðsÞAðsÞ; K 33 ¼ EðsÞAðsÞ; > > R R > 2 > > J 11 ¼ AðsÞ EðsÞg dA; J 22 ¼ AðsÞ EðsÞn2 dA; > > > R R < I 11 ¼ AðsÞ qðsÞg2 dA; I 22 ¼ AðsÞ qðsÞn2 dA; ð11Þ > R > 2 oW oW 2 > > J 33 ¼ AðsÞ GðsÞ n þ g þ n og þ g on dA; > > > > R : I 33 ¼ AðsÞ qðsÞðn2 þ g2 Þ dA; where W is the warping function which is the same for all cross-sections according to the Saint-Venant Theory, and j(s) is a numerical factor depending on the shape of the cross-section at s. Making use of (5), (9) and (10), we have 8 > < n ¼ K 11 v1 d1 þ K 22 v2 d2 þ K 33 ðv3 1Þd3 ; m ¼ J 11 u1 d1 þ J 22 u2 d2 þ J 33 u3 d3 ; ð12Þ > : h ¼ I 11 w1 d1 þ I 22 w2 d2 þ I 33 w3 d3 : Comparing with (7), we get 8 m1 ¼ J 11 u1 ; > < n1 ¼ K 11 v1 ; n2 ¼ K 22 v2 ; m2 ¼ J 22 u2 ; > : n3 ¼ K 33 ðv3 1Þ; m3 ¼ J 33 u3 ;
h1 ¼ I 11 w1 ; h2 ¼ I 22 w2 ;
ð13Þ
h3 ¼ I 33 w3 :
3. Specifications for the deformed configuration space Consider a uniform and initially straight rod of constant length ‘, supported in an arbitrary manner at each end as shown in Fig. 2. It is assumed that the static equilibrium of the rod corresponds to the situation where the directions of d3 and e3 coincide, d1 and d2 are parallel to e1 and e2, respectively. The components of the elastic displacement vector of the centroid at an arbitrary location s is assumed to be x(s, t), y(s, t), z(s, t), i.e. d1
e1
d3
s
r O
d2 A
e3 e2 Fig. 2. Coordinate systems used in the development of governing equations.
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
rðs; tÞ ¼ xðs; tÞe1 þ yðs; tÞe2 þ ðs þ zðs; tÞÞe3 :
465
ð14Þ
The orientation of the directors d1, d2, d3 at location s, relative to the inertial basis {e1, e2, e3}, can be described by three successive rotations, called Euler angles. Alternatively, the moving frame can be described by the rotation group SO(3) as the formulations in Simo and Vu-Quoc (1988, 1991), Vu-Quoc and Deng (1995) and Vu-Quoc and Ebcioglu (1995). The three successive rotations starts by rotating the basis {e1, e2, e3} an angle w(s, t) about the axis aligned with the director e2 as shown in Fig. 3a. Next, we rotate the basis f^e1 ; e2 ; ^e3 g an angle h(s, t) about the axis aligned with the new director ^e1 as shown in Fig. 3b. Finally, we rotate the basis f^e1 ; ^e2 ; d3 g an angle /(s, t) about the axis aligned with the new director d3 as shown in Fig. 3c. The representative matrices for the three successive rotations can be written as 2 3 cos wðs; tÞ 0 sin wðs; tÞ 6 7 Dw ðs; tÞ ¼ 4 0 1 0 ð15Þ 5; sin wðs; tÞ
2
1 6 Dh ðs; tÞ ¼ 4 0 0 and
0
0 cos hðs; tÞ sin hðs; tÞ
cos wðs; tÞ
3 0 7 sin hðs; tÞ 5;
ð16Þ
cos hðs; tÞ
2
cos /ðs; tÞ sin /ðs; tÞ 6 D/ ðs; tÞ ¼ 4 sin /ðs; tÞ cos /ðs; tÞ 0
0
3 0 7 0 5;
ð17Þ
1
respectively. Denote dij = di Æ ej. Then, using Eqs. (15)–(17), one can write the transformation matrix between the inertial basis {e1, e2, e3} and the basis {d1, d2, d3} as follows: Dðs; tÞ ¼ ½d ij ðs; tÞ ¼ D/ ðs; tÞDh ðs; tÞDw ðs; tÞ 2 cos / cos w þ sin / sin h sin w sin / cos h 6 ¼ 4 sin / cos w þ cos / sin h sin w cos / cos h
3 cos / sin w þ sin / sin h cos w 7 sin / sin w þ cos / sin h cos w 5:
sin h
cos h sin w
cos h cos w
Since D is an orthogonal matrix, from the definition of dij = di Æ ej, we have
e1 e1
e1
ψ
e1
e3 ψ
d1 o (e 2 )
φ ψ
o
θ
φ
e3
d3 θ ψ
e2
(a)
e3
d3
e3
(b)
θ
e2
e3
θ
o(e1 ) d2
φ
e2
e2
(c) d1
e2
φ
d2 o (d 3)
e1
Fig. 3. The successive angle rotations of a differential element of arc length ds.
ð18Þ
466
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
di ¼
3 X
ej ¼
d ij ej ;
j¼1
ð19Þ
d kj dk :
k¼1
Thus os di ¼ os
3 X
3 X
! d ij ej
j¼1
¼
3 X
os d ij ej ¼
3 X 3 X
j¼1
j¼1
ðos d ij Þd kj dk :
Substituting (20) into (3) yields 3 3 1X 1X di os di ¼ di u¼ 2 i¼1 2 i¼1
¼
3 X 3 X j¼1
ð20Þ
k¼1
! ðos d ij Þd kj dk
k¼1
3 3 1X 1X ½ðos d 1j Þd 2j ðos d 2j Þd 1j d3 þ ½ðos d 2j Þd 3j ðos d 3j Þd 2j d1 2 j¼1 2 j¼1
þ
3 1X ½ðos d 3j Þd 1j ðos d 1j Þd 3j d2 : 2 j¼1
ð21Þ
Taking notice of 3 X
d ij ðs; tÞd kj ðs; tÞ ¼ dik ;
ð22Þ
j¼1
where d is the Kronecker d-function, we have 3 X
½ðos d ij ðs; tÞÞd kj ðs; tÞ þ d ij ðs; tÞðos d kj ðs; tÞÞ ¼ 0:
ð23Þ
j¼1
Hence, (21) can be written as u¼
3 X ½ðos d 2j Þd 3j d1 þ ðos d 3j Þd 1j d2 þ ðos d 1j Þd 2j d3 :
ð24Þ
j¼1
Therefore, 8 > < u1 u2 > : u3
from (18) and (24) we can obtain ¼ os h cos / þ os w cos h sin /; ð25Þ
¼ os h sin / þ os w cos h cos /; ¼ os / os w sin h:
Similarly, replacing the time derivative by the space derivative, we can easily obtain 8 > < w1 ¼ ot h cos / þ ot w cos h sin /; w2 ¼ ot h sin / þ ot w cos h cos /; > : w3 ¼ ot / ot w sin h:
ð26Þ
In order to express the extension and shear strains as functions of displacements x(s, t), y(s, t), z(s, t) and angles h(s, t), w(s, t), /(s, t), we define re1 ðs; tÞ ¼ xðs; tÞ;
re2 ðs; tÞ ¼ yðs; tÞ;
re3 ðs; tÞ ¼ s þ zðs; tÞ:
ð27Þ
Then, the strain vector v(s, t) can be written as v ¼ os r ¼
3 3 X X ðos rei Þei ¼ ðos rei Þd ji dj : i¼1
i¼1
Substituting (18) and (27) into (28) yields
ð28Þ
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
8 v1 ¼ os xðcos / cos w þ sin / sin h sin wÞ þ os y sin / cos h > > > > > þð1 þ os zÞð cos / sin w þ sin / sin h cos wÞ; > < v2 ¼ os xð sin / cos w þ cos / sin h sin wÞ þ os y cos / cos h > > > þð1 þ os zÞðsin / sin w þ cos / sin h cos wÞ; > > > : v3 ¼ os x cos h sin w os y sin h þ ð1 þ os zÞ cos h cos w:
467
ð29Þ
4. Governing equations of motion To obtain the partial differential equations of motion in terms of the displacements x(s, t), y(s, t), z(s, t) and angles h(s, t), w(s, t), /(s, t), we decompose the contact force n(s, t) with respect to the fixed basis {e1, e2, e3} as nðs; tÞ ¼
3 X
nei ðs; tÞei :
ð30Þ
i¼1
Then, making use of (9) and (19), we have n¼
3 X
ni di ¼
i¼1
3 X 3 X i¼1
3 X
ni d ij ej ¼
j¼1
nej ej :
ð31Þ
i¼1
It follows from (13) that 8 3 P > > > ne1 ¼ d i1 ni ¼ d 11 K 11 v1 þ d 21 K 22 v2 þ d 31 K 33 ðv3 1Þ; > > > i¼1 > > < 3 P ne2 ¼ d i2 ni ¼ d 12 K 11 v1 þ d 22 K 22 v2 þ d 32 K 33 ðv3 1Þ; > i¼1 > > > > 3 > P > > d i3 ni ¼ d 13 K 11 v1 þ d 23 K 22 v2 þ d 33 K 33 ðv3 1Þ: : ne3 ¼
ð32Þ
i¼1
To proceed, using the relations between the strain vector u(s, t) and the differentiation of the directors d1(s, t), d2(s, t) and d3(s, t) with respect to s, and the relations between the angular velocity vector w(s, t) and the differentiation of the directors d1(s, t), d2(s, t) and d3(s, t) with respect to t, we can easily obtain the following expressions: os m ¼
3 X
3 X
os ðmi di Þ ¼
i¼1
ðos mi di þ mi u di Þ
i¼1
¼ ðos m1 m2 u3 þ m3 u2 Þd1 þ ðos m2 m3 u1 þ m1 u3 Þd2 þ ðos m3 m1 u2 þ m2 u1 Þd3 ; ot h ¼
3 X
ð33Þ
3 X ot ðhi di Þ ¼ ðot hi di þ hi w di Þ
i¼1
i¼1
¼ ðot h1 h2 w3 þ h3 w2 Þd1 þ ðot h2 h3 w1 þ h1 w3 Þd2 þ ðot h3 h1 w2 þ h2 w1 Þd3 : In addition, we record ! 3 X vn¼ v i di i¼1
3 X j¼1
ð34Þ
! nj dj
¼ ðv2 n3 v3 n2 Þd1 þ ðv3 n1 v1 n3 Þd2 þ ðv1 n2 v2 n1 Þd3 :
ð35Þ
468
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
Now, let us assume that the external force and torque are distributed loadings with fixed direction and prescribed intensity. Consequently, the external force and torque densities can be written as: fðs; tÞ ¼ fx ðs; tÞe1 þ fy ðs; tÞe2 þ fz ðs; tÞe3 ; ð36Þ lðs; tÞ ¼ lx ðs; tÞe1 þ ly ðs; tÞe2 þ lz ðs; tÞe3 : Decomposing the torque density with respect to the moving basis {d1, d2, d3}, we have lðs; tÞ ¼
3 X
li ðs; tÞdi ðs; tÞ;
ð37Þ
i¼1
where, li(s,t) (i = 1, 2, 3) can be expressed by using the transformation matrix (18) and the relation (19) as 8 l1 ¼ lx ðcos / cos w þ sin / sin h sin wÞ þ ly sin / cos h > > > > > þlz ð cos / sin w þ sin / sin h cos wÞ; > < l2 ¼ lx ð sin / cos w þ cos / sin h sin wÞ þ ly cos / cos h ð38Þ > > > þlz ðsin / sin w þ cos / sin h cos wÞ; > > > : l3 ¼ lx cos h sin w ly sin h þ lz cos h cos w: Substituting the expressions (32)–(37) into the dynamic equation (6), we can obtain 8 qAott x ¼ os ðd 11 K 11 v1 þ d 21 K 22 v2 þ d 31 K 33 ðv3 1ÞÞ þ fx ; > > > > > qAo tt y ¼ os ðd 12 K 11 v1 þ d 22 K 22 v2 þ d 32 K 33 ðv3 1ÞÞ þ fy ; > > > < qAo z ¼ o ðd K v þ d K v þ d K ðv 1ÞÞ þ f ; tt s 13 11 1 23 22 2 33 33 3 z > ot h1 h2 w3 þ h3 w2 ¼ os m1 m2 u3 þ m3 u2 þ v2 n3 v3 n2 þ l1 ; > > > > > ot h2 h3 w1 þ h1 w3 ¼ os m2 m3 u1 þ m1 u3 þ v3 n1 v1 n3 þ l2 ; > > : ot h3 h1 w2 þ h2 w1 ¼ os m3 m1 u2 þ m2 u1 þ v1 n2 v2 n1 þ l3 ;
ð39Þ
where dij is defined by (18), ni, mi and hi are given by (13), while ui, wi, vi and li are given by (25), (26), (29) and (38), respectively. The boundary conditions at the ends of the rod should be determined in each particular case. For example, at a fixed end the deflections and the rotations are equal to zero. In this case the boundary conditions are x¼y¼z¼0
and
h ¼ w ¼ / ¼ 0:
ð40Þ
At a free end the boundary conditions are am n þ f b ðtÞ ¼ 0
and
am m þ lb ðtÞ ¼ 0;
where am = +1 at the right end of the rod, am = 1 at the left end, and ( b f ðtÞ ¼ fxb ðtÞe1 þ fyb ðtÞe2 þ fzb ðtÞe3 ; lb ðtÞ ¼ lbx ðtÞe1 þ lby ðtÞe2 þ lbz ðtÞe3
ð41Þ
ð42Þ
are the external force and torque applied to the corresponding end of the rod. As the external distributed torque l(s, t), the torque vector lb(t) can be expressed in the moving frame {d1, d2, d3} and the corresponding components lbi ðtÞ ði ¼ 1; 2; 3Þ can be determined through the relation (19) and the transformation matrix (18) as follows: 8 b l1 ¼ lbx ðcos / cos w þ sin / sin h sin wÞ þ lby sin / cos h > > > > > > þlbz ð cos / sin w þ sin / sin h cos wÞ; > < b l2 ¼ lbx ð sin / cos w þ cos / sin h sin wÞ þ lby cos / cos h ð43Þ > > > b > þlz ðsin / sin w þ cos / sin h cos wÞ; > > > : lb ¼ lb cos h sin w lb sin h þ lb cos h cos w; 3
x
y
z
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
469
with w = w(s*, t), h = h(s*, t) and / = /(s*, t) at the boundary s = s* (s* = 0 or ‘). Substituting the relations (32) and (42) into the first equation of (41), the force boundary conditions can be rewritten as 8 b > < am ½d 11 K 11 v1 þ d 21 K 22 v2 þ d 31 K 33 ðv3 1Þ þ fx ðtÞ ¼ 0; am ½d 12 K 11 v1 þ d 22 K 22 v2 þ d 32 K 33 ðv3 1Þ þ fyb ðtÞ ¼ 0; ð44Þ > : b am ½d 13 K 11 v1 þ d 23 K 22 v2 þ d 33 K 33 ðv3 1Þ þ fz ðtÞ ¼ 0: Similarly, the torque boundary conditions can be obtained as am J 11 u1 þ lb1 ðtÞ ¼ 0;
am J 22 u2 þ lb2 ðtÞ ¼ 0;
am J 33 u3 þ lb3 ðtÞ ¼ 0:
ð45Þ
For simplicity, it will be assumed that the axes along the directors d1, d2 and d3 are chosen to be the principal axes of inertia of the cross-section at s, and centered at the cross-section’s center of mass. Then, we have 8 > < K 11 ¼ jGA; K 22 ¼ jGA; K 33 ¼ EA; J 11 ¼ Ec11 ; J 22 ¼ Ec22 ; J 33 ¼ Gc33 ; ð46Þ > : I 11 ¼ qc11 ; I 22 ¼ qc22 ; I 33 ¼ qðc11 þ c22 Þ; where c11 and c22 are the moments of inertia of the rod cross-section, and c33 is the torsional moment of inertia. For a rectangular cross-section, the torsional moment of inertia is as follows (Rubin, 2000): " # 1 WH 3 192H X 1 pð2n 1ÞW c33 ¼ 1 5 tanh for W P H ; p W n¼1 ð2n 1Þ5 2H 3 where W and H are the dimensions of the cross-section. Moreover, the approximate formula of the torsional moment of inertia can be written as WH 3 H H4 1 1 0:63 c33 ¼ for W P H : W 3 12W 4 As a prelude to expanding the nonlinear partial differential equations (39) to a form suitable for a perturbation analysis of the motion, it is useful to introduce some natural scales to obtain a dimensionless equation of motion. Introduce the dimensionless variables s r¼ ; ‘
x p1 ¼ ; ‘
y p2 ¼ ; ‘
z p3 ¼ ; ‘
p4 ¼ h;
p5 ¼ w;
p6 ¼ /;
s ¼ x0 t;
where x0 is the reference natural frequency yet to be determined. Moreover, let 8 0 0 g0 ¼ qAx1 2 ‘ ; n01 ¼ jGAg ; n02 ¼ EAg ; > ‘ ‘ > > 0 > > > > g3 ¼ qc 1 x2 ; g1 ¼ qc 1 x2 ; g2 ¼ qc 1 x2 ; > > 11 0 22 0 33 0 > < n11 ¼ q‘2Ex2 ; n12 ¼ Gcq‘332 cEcx222 ; n13 ¼ jGAg1 ; n14 ¼ EAg1 ; 11 0 0 > > > Gc33 Ec11 E > n21 ¼ q‘2 x2 ; n22 ¼ q‘2 c x2 ; n23 ¼ jGAg2 ; n24 ¼ EAg2 ; > > 22 0 > 0 > > > : n31 ¼ G2 2 ; n32 ¼ Eðc222 c112 Þ ; n33 ¼ c22 c11 : c33 q‘ x q‘ c x 0
33
ð47Þ
ð48Þ
0
Then, substituting the relations (18), (13), (25), (26), (29), (38) and (46) into the nonlinear equation (39), we obtain the governing differential equation of motion in terms of dimensionless displacements pi(s, t) (i = 1, 2, . . . , 6) as follows:
470
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
€ p1 ¼ fðn01 n02 Þðp01 cos p4 sin p5 þ p02 sin p4 ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 sin p5 0
þ n01 p01 n02 cos p4 sin p5 g þ g0 fx ; € p2 ¼ fðn01 n02 Þðp1 0 sin p4 sin p5 þ þ
n02 ðp02
p02
ð49Þ
cos p4 þ ð1 þ
p03 Þ sin p4
cos p5 Þ cos p4
0
þ sin p4 Þg þ g0 fy ;
ð50Þ
€ p3 ¼ fðn01 n02 Þðp01 cos p4 sin p5 þ p02 sin p4 þ n01 ð1 þ p03 Þ n02 cos p4 cos p5 g0 þ g0 fz ;
ð1 þ
p03 Þ cos p4
cos p5 Þ cos p4 cos p5 ð51Þ
:
ðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ ðp_ 6 p_ 5 sin p4 Þðp_ 4 sin p6 p_ 5 cos p4 cos p6 Þ 0
¼ n11 ðp04 cos p6 þ p05 sin p6 cos p4 Þ n12 ðp06 p05 sin p4 Þðp04 sin p6 p05 cos p6 cos p4 Þ ððn13 n14 Þðp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ þ n14 Þ fp01 ð sin p6 cos p5 þ cos p6 sin p4 sin p5 Þ þ p02 cos p6 cos p4 þ ð1 þ p03 Þðsin p6 sin p5 þ cos p6 sin p4 cos p5 Þg þ g1 lx ðcos p6 cos p5 þ sin p6 sin p4 sin p5 Þ þ g1 ly sin p6 cos p4 þ g1 lz ð cos p6 sin p5 þ sin p6 sin p4 cos p5 Þ; : ðp_ 4 sin p6 þ p_ 5 cos p6 cos p4 Þ ðp_ 6 p_ 5 sin p4 Þðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ
ð52Þ
0
¼ n21 ðp05 cos p6 cos p4 p04 sin p6 Þ n22 ðp06 p05 sin p4 Þðp04 cos p6 þ p05 sin p6 cos p4 Þ þ ððn23 n24 Þðp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ þ n24 Þ fp01 ðcos p6 cos p5 þ sin p6 sin p4 sin p5 Þ þ p02 sin p6 cos p4 þ ð1 þ p03 Þð cos p6 sin p5 þ sin p6 sin p4 cos p5 Þg þ g2 lx ð sin p6 cos p5 þ cos p6 sin p4 sin p5 Þ þ g2 ly cos p6 cos p4 þ g2 lz ðsin p6 sin p5 þ cos p6 sin p4 cos p5 Þ; ðp_ 6 p_ 5 sin p4 Þ: þ n33 ðp_ 4 sin p6 þ p_ 5 cos p4 cos p6 Þðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ
ð53Þ
0
¼ n31 ðp06 p05 sin p4 Þ þ n32 ðp04 sin p6 þ p05 cos p6 cos p4 Þðp04 cos p6 þ p05 sin p6 cos p4 Þ þ g3 ðlx cos p4 sin p5 ly sin p4 þ lz cos p4 cos p5 Þ;
ð54Þ
where the symbols (.) and ( 0 ) denote differentiation with respect to dimensionless time parameter s and dimensionless length parameter r, respectively. The boundary conditions at a fixed end are p1 ¼ p2 ¼ p3 ¼ p4 ¼ p5 ¼ p6 ¼ 0:
ð55Þ
At a free end, the force boundary conditions (44) can be rewritten as am fðn01 n02 Þðp01 cos p4 sin p5 þ p02 sin p4 ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 sin p5 g þ n01 p01 n02 cos p4 sin p5 g þ 0 fxb ðsÞ ¼ 0; ‘ am fðn01 n02 Þðp01 sin p4 sin p5 þ p02 cos p4 þ ð1 þ p03 Þ sin p4 cos p5 Þ cos p4 g þ n02 ðp02 þ sin p4 Þg þ 0 fyb ðsÞ ¼ 0; ‘ am fðn01 n02 Þðp01 cos p4 sin p5 þ p02 sin p4 ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 cos p5 g þ n01 ð1 þ p03 Þ n02 cos p4 cos p5 g þ 0 fzb ðsÞ ¼ 0; ‘ while the torque boundary conditions (45) can be rewritten as
ð56Þ
ð57Þ
ð58Þ
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
am n11 ðp04 cos p6 þ p05 sin p6 cos p4 Þ þ
471
g1 b l ðcos p6 cos p5 þ sin p6 sin p4 sin p5 Þ ‘ x
g1 b g l sin p6 cos p4 þ 1 lbz ð cos p6 sin p5 þ sin p6 sin p4 cos p5 Þ ¼ 0; ‘ y ‘ g 0 0 am n21 ðp5 cos p6 cos p4 p4 sin p6 Þ þ 2 lbx ð sin p6 cos p5 þ cos p6 sin p4 sin p5 Þ ‘ g g þ 2 lby cos p6 cos p4 þ 2 lbz ðsin p6 sin p5 þ cos p6 sin p4 cos p5 Þ ¼ 0; ‘ ‘ g3 b 0 0 am n31 ðp6 p5 sin p4 Þ þ ðlx cos p4 sin p5 lby sin p4 þ lbz cos p4 cos p5 Þ ¼ 0: ‘ þ
ð59Þ
ð60Þ ð61Þ
5. Femlab implementation Femlab is an open flexible system used to model all types of scientific and engineering problems based on partial differential equations (Femlab homepage). This gives us freedoms to implement numerical solution procedure of the nonlinear partial differential equations of motion with the corresponding boundary conditions and initial conditions. We introduce the solution procedure of the nonlinear partial differential equations (49)–(54) with the boundary conditions (55) and/or (56)–(61). Let 2 b 3 2 3 2 3 2 3 2 3 fx ðsÞ p_ 1 b1 g1 p1 6 b 7 6p 7 6 p_ 7 6b 7 6g 7 6 f ðsÞ 7 6 27 6 27 6 27 6 27 7 6 y 6 7 6 7 6 7 6 7 6 6 p3 7 6 p_ 3 7 6 b3 7 6 g3 7 g0 6 fzb ðsÞ 7 7 b 7 7 7 7 6 6 6 6 ð62Þ p ¼ 6 7; q ¼ p_ ¼ 6 7; b ¼ 6 7; g ¼ 6 7; g ¼ 6 b 7: ‘ 6 l1 ðsÞ 7 6 p4 7 6 p_ 4 7 6 b4 7 6 g4 7 7 6 6 7 6 7 6 7 6 7 6 lb ðsÞ 7 4 p5 5 4 p_ 5 5 4 b5 5 4 g5 5 5 4 2 b p6 b6 g6 p_ 6 l3 ðsÞ where 8 b1 ¼ ðn01 n02 Þðþp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 sin p5 > > > > > n01 p01 þ n02 cos p4 sin p5 ; > > > > > b2 ¼ ðn02 n01 Þðp01 sin p4 sin p5 þ p02 cos p4 þ ð1 þ p03 Þ sin p4 cos p5 Þ cos p4 > > > < n02 ðp02 þ sin p4 Þ; > b3 ¼ ðn01 n02 Þðp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ cos p4 cos p5 > > > > > n01 ð1 þ p03 Þ þ n02 cos p4 cos p5 ; > > > > > b4 ¼ n11 ðp04 cos p6 þ p05 sin p6 cos p4 Þ; > > : b5 ¼ n21 ðp04 sin p6 p05 cos p6 cos p4 Þ; b6 ¼ n31 ðp06 þ p05 sin p4 Þ: Taking notice of 8 ðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ: ðp_ 6 p_ 5 sin p4 Þðp_ 4 sin p6 p_ 5 cos p4 cos p6 Þ > > > > > ¼€ p4 cos p6 þ € p5 sin p6 cos p4 2p_ 4 p_ 6 sin p6 > > > > > _ _ þ2 p cos p4 cos p6 p_ 25 sin p4 cos p4 cos p6 ; p > 5 6 > > > : > > < ðp_ 4 sin p6 þ p_ 5 cos p6 cos p4 Þ ðp_ 6 p_ 5 sin p4 Þðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ ¼ € p4 sin p6 þ € p5 cos p6 cos p4 2p_ 4 p_ 6 cos p6 > > > 2p_ 5 p_ 6 sin p6 cos p4 þ p_ 25 sin p4 cos p4 sin p6 ; > > > > : > > ðp_ 6 p_ 5 sin p4 Þ þ n33 ðp_ 4 sin p6 þ p_ 5 cos p4 cos p6 Þðp_ 4 cos p6 þ p_ 5 sin p6 cos p4 Þ > > > > > p6 p€5 sin p4 Þ p_ 5 p_ 4 cos p4 > ¼ ð€ > : þn33 ðp_ 24 sin p6 cos p6 þ p_ 4 p_ 5 cos 2p6 cos p4 þ p_ 25 sin p6 cos p6 cos2 p4 Þ; the nonlinear partial differential equations (49)–(54) can be rewritten as
ð63Þ
ð64Þ
472
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
p_
þ
M a q_
0
b0
q ¼ ; g
where Ma ¼
I 33 Mm
;
ð65Þ 2
cos p6
sin p6 cos p4
0
cos p6 cos p4 sin p4
6 M m ¼ 4 sin p6
0
3
7 0 5; 1
ð66Þ
and 8 g1 > > > > > g4 > > > > > > > > > > > > > > > > > > > > > < g5 > > > > > > > > > > > > > > > > > > > > > > > > g6 > > :
¼ g0 fx ; ¼
g 2 ¼ g0 f y ;
n12 ðp06
p05
g3 ¼ g0 fz ;
sin p4 Þðp04 sin p6 p05 cos p6 cos p4 Þ
ððn13 n14 Þðp01 cos p4 sin p5 p02 sin p4 þ ð1 þ p03 Þ cos p4 cos p5 Þ þ n14 Þ fp01 ð sin p6 cos p5 þ cos p6 sin p4 sin p5 Þ þ p02 cos p6 cos p4 þð1 þ p03 Þðsin p6 sin p5 þ cos p6 sin p4 cos p5 Þg þ g1 l1 þ2p_ 4 p_ 6 sin p6 2p_ 5 p_ 6 cos p4 cos p6 þ p_ 25 sin p4 cos p4 cos p6 ; ¼ n22 ðp06 p05 sin p4 Þðp04 cos p6 þ p05 sin p6 cos p4 Þ n24 Þðp01
p02
ð67Þ p03 Þ cos p4
þððn23 cos p4 sin p5 sin p4 þ ð1 þ cos p5 Þ þ n24 Þ fp01 ðcos p6 cos p5 þ sin p6 sin p4 sin p5 Þ þ p02 sin p6 cos p4 þð1 þ p03 Þð cos p6 sin p5 þ sin p6 sin p4 cos p5 Þg þ g2 l2 þ2p_ 4 p_ 6 cos p6 þ 2p_ 5 p_ 6 sin p6 cos p4 p_ 25 sin p4 cos p4 sin p6 ; ¼ n32 ðp04 sin p6 þ p05 cos p6 cos p4 Þðp04 cos p6 þ p05 sin p6 cos p4 Þ þ g3 l3 þ p_ 5 p_ 4 cos p4 þn33 ðp_ 24 sin p6 cos p6 p_ 4 p_ 5 cos 2p6 cos p4 p_ 25 sin p6 cos p6 cos2 p4 Þ:
Using the general form for nonlinear PDEs with 12 dependent variables, Eq. (65) can be implemented in Femlab in terms of the Matlab code (for details, see Comsol, 2004). The Dirichlet boundary conditions should be adopted at a fixed end, i.e., p = 0. At a free end, the Neumann boundary conditions should be adopted, i.e., amb = gb, where b and gb are defined by (62). For the formulated nonlinear model of a Cosserat rod, the numerical analysis using the Femlab programming language at the Matlab command line consists of the following steps: 1. 2. 3. 4.
Formulation of the governing equations; Discretization of the equations; Solution of the equations; Interpretation of the results.
Given the one-dimensional geometry data, an initial finite element mesh is easily generated by meshinit. The meshextend command is used for discretization of the PDE problem and can be modified to improve accuracy. The geometry, nonlinear PDEs and boundary conditions are defined by a set of fields in the Matlab codes. For solving purposes Femlab contains specific solvers (like static, dynamic, linear, nonlinear solvers) for specific PDE problems. Here, we use femtime to solve the dynamic problem of the system. 6. Simulation results and discussion for a simple cantilever A cantilever (fixed at the left hand and free at the right hand), as shown in Fig. 2, is now presented as a simple example to demonstrate the numerical solution procedure of the proposed Cosserat rod model. Numerical calculations based on (65) are carried out for the uniform horizontal cantilever of length ‘ = 1.0 m, of constant cross-section with width b = 0.06 m and thickness a = 0.04 m. The mass density, the Young’s modulus and the Poisson’s ratio are assumed to be q = 2.730 · 103 kg/m3, E = 7.1 · 1010 Pa and m = 0.32, respectively. The shear correction factor for rectangular cross-section is taken to be j = 0.833. The reference frequency is
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
473
chosen to be the lowest natural frequency of the linearized system of a cantilever without considering the effect of shearing deformations, i.e., sffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffi 1:8752 Ec22 1:8752 Ea2 x0 ¼ ¼ ¼ 207:0236 rad s1 : qA 12q ‘2 ‘2 Based on the derived nonlinear system (65), numerical simulations are performed to investigate the dynamic responses of the cantilever under harmonic excitations. The differential equations of motion are full coupled by the nonlinear terms and could exhibit internal resonance introduced by the nonlinearities. They also exhibit external resonances when the external excitation is periodic and the frequency of a component of its Fourier series is near one of the natural frequencies of the system, or near a multiple of a natural frequencies. The detailed analysis of complex dynamic behavior, such as bifurcation (Vu-Quoc and Li, 1995) and chaos, of the system is not the main focus of this paper. We only present here the simulation results of the dynamic responses of the system under external harmonic excitations. The dynamical responses of the free end of the cantilever (with zero initial displacements and velocities) under external distributing loads fx ðtÞ ¼ 2 sinðpsÞ sinð8x0 tÞ kN m1 ;
f y ðtÞ ¼ 2 sinðpsÞ sinð8x0 tÞ kN m1 ;
are obtained using the Femlab/Matlab, and shown in Fig. 4. It is observed that the transverse displacements x(1, t) and y(1, t) consist of more than one harmonic motions with different frequencies. The longitudinal displacement z(1, t) and the torsional oscillation /(1, t) are obviously not periodic motions, which is not surprising, since various dynamic phenomena are expected for such a strongly nonlinear and fully coupling system. Moreover, it can be observed that the longitudinal displacement z(1, t) is a higher order small term comparing with other variables like x(1, t) and y(1, t). The transverse displacement time histories x(s, t) and y(s, t) of the nonlinear coupling model under external harmonic excitations are shown in Figs. 5 and 6, respectively. Dividing the cantilever into 10 Cosserat Rod Elements (CRE) of equal length, see Cao et al. (2006) for the formulation of CREs, one can establish the nonlinear ordinary differential equations of motion for solving the free displacements. Numerical simulations for the responses of the model under the same external harmonic −3
−3
x 10
θ(1,t) [rad]
x(1,t) [m]
2 0 −2
0
0.02
0.04
0.06
5 0 −5
−3
0 −2
0
0.02
0.04
0.06
5
φ(1,t) [rad]
z(1,t) [m]
0 −4
0
0.04
0.06
0.02
0.04
0.06
0.04
0.06
x 10
−5
0 −5
x 10
−2 −6
0.02
0
−6
2
0 −3
x 10
ψ(1,t) [rad]
y(1,t) [m]
2
x 10
0.02
0.04
t [sec]
0.06
2
x 10
0 −2 −4
0
0.02
t [sec]
Fig. 4. Dynamical responses of the cantilever at the free end under harmonic excitations.
474
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477 −3
x 10 1.5 1
x(s,t) [m]
0.5 0 −0.5 −1 −1.5 0.06 0.04
t [sec]
0.02 0
0.2
0
0.4
0.6
0.8
1
s [m]
Fig. 5. Displacement time histories of the cantilever x(s, t) [m] under harmonic excitations.
−3 x 10 1
y(s,t) [m]
0.5 0 −0.5 −1 −1.5 0.06 1
0.04
t [sec]
0.8 0.6
0.02
0.4 0
0.2
s [m]
0
Fig. 6. Displacement time histories of the cantilever y(s, t) [m] under harmonic excitations.
excitations are performed with Matlab. The corresponding dynamical responses of the CRE model are depicted in Figs. 7–9 for a comparison of the results in Figs. 4–6 obtained in this paper. 7. Conclusions The nonlinear partial differential equations of motion for slender elastic rods have been formulated using the Cosserat theory. Unlike other formulations presented in the literature, the nonlinear dynamic model formulated here are explicitly in terms of the displacements and angular variables. Flexure along two principal axes, extension, shear and torsion are considered, and care is taken to account for all the nonlinearities include contributions from both the curvature expression and inertia terms.
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
θ(1,t) [rad]
5
0
0.02
x 10
0.04
−3
5
0 −2
2
0
0.02
x 10
0.04
−6
2
0 −2 −4 −6
0
0.02
0.04
0 x 10
0.02
0.04
0.06
0.02
0.04
0.06
0.04
0.06
−3
0 x 10
−5
0 −2 −4
0.06
−3
0 −5
0.06
x 10
0 −5
0.06
ψ(1,t) [rad]
2
y(1,t) [m]
−3
0 −2
z(1,t) [m]
x 10
φ(1,t) [rad]
x(1,t) [m]
2
475
0
0.02
t [sec]
t [sec]
Fig. 7. Dynamical responses of the cantilever at the free end under harmonic excitations: Cosserat Rod Element Approach.
−3
x 10 1.5
x(s,t) [m]
1 0.5 0 −0.5 −1 −1.5 0.06 0.04
t [sec]
0.02 0
0
0.2
0.4
0.6
0.8
1
s [m]
Fig. 8. Displacement time histories of the cantilever x(s, t) [m] under harmonic excitations: Cosserat Rod Element Approach.
A numerical solution procedure in terms of the Femlab/Matlab interfaces has been presented for numerically solving the formulated nonlinear dynamical model with the corresponding boundary and initial conditions. A simple cantilever under external harmonic excitations has been presented to demonstrate the use of the formulation developed, and to illustrate the numerical solution procedure. Based on the case study of applying the proposed modelling strategy for nonlinear dynamics of elastic rods, the following conclusions have been drawn.
476
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
−3
x 10 1
y(s,t) [m]
0.5 0 −0.5 −1 −1.5 0.06 1
0.04
t [sec]
0.8 0.6
0.02
0.4 0
0.2 0
s [m]
Fig. 9. Displacement time histories of the cantilever y(s, t) [m] under harmonic excitations: Cosserat Rod Element Approach.
1. The explicit nature of the construction permits us to directly express the nonlinear dynamic model of the elastic rod and the corresponding boundary conditions in terms of the displacements and angular variables, thus the dynamical analysis of elastic rods can be carried out rather simply. 2. The governing partial differential equations of motion can be easily expanded to contain nonlinearities up to order three to render them amenable to the study of moderately large amplitude flexural-torsional-shearable oscillations by perturbation techniques. Moreover, the approach on dimension reduction of dynamical systems, such as the nonlinear Galerkin method (Rega and Troger, 2005) and the Karhunen–Loeve method (which in the mathematical literature is also called Proper Orthogonal Decomposition (POD) method (Georgiou, 2005)), could be used to develop reduced-order dynamic models of the coupled nonlinear partial differential systems derived here. Therefore, the full coupled nonlinear model formulated here could be very useful for the detailed analysis of complex dynamic behavior, such as bifurcation and chaos, of the slender rods. 3. The nonlinear dynamic model could be numerically solved in terms of the Femlab/Matlab interfaces. The Femlab environment offers a powerful simulation tool for analyzing the responses of nonlinear dynamic models (of slender material elements such as those found in an MEMS device) presented here.
Acknowledgements The authors are grateful to the EPSRC (Computational Engineering Mathematics Programme) for financial support in this study. References Antman, S.S., 1991. Nonlinear Problems of Elasticity. In: Applied Mathematical Sciences. Springer-Verlag, New York. Antman, S.S., Marlow, R.S., Vlahacos, C.P., 1998. The complicated dynamics of heavy rigid bodies attached to deformable rods. Quarterly of Applied Mathematics 86, 431–460. Bishop, R.E.D., Price, W.G., 1977. Coupled bending and twisting of a Timoshenko beam. Journal of Sound and Vibration 50, 469–477. Bolotin, V.V., 1964. The Dynamic Stability of Elastic Systems. Holden-Day, San Francisco. Bruch Jr., J.C., Mitchell, T.P., 1987. Vibrations of mass-loaded clamped-free Timoshenko beam. Mitchell Journal of Sound and Vibration 114 (2), 341–345. Cao, D.Q., Liu, D., Wang, C.H.-T., 2005. Nonlinear dynamic modelling for MEMS components via the Cosserat rod element approach. Journal of Micromechanics and Microengineering 15, 1334–1343.
D.Q. Cao, R.W. Tucker / International Journal of Solids and Structures 45 (2008) 460–477
477
Cao, D.Q., Liu, D., Wang, C.H.-T., 2006. Three-dimensional nonlinear dynamics of slender structures: Cosserat rod element approach. International Journal of Solids and Structures 43, 760–783. Cartmell, M.P., 1990. The equations of motion for a parametrically excited cantilever beam. Journal of Sound and Vibration 143 (3), 395– 406. Comsol, AB., 2004. Femlab 3.0: Matlab Interface Guide. COMSOL Ltd, Oxford. Crespo da Silva, M.R.M., 1988a. Nonlinear flexural-flexural-torsional-extensional dynamics of beams – I. Formulation. International Journal of Solids and Structures 24 (12), 1225–1234. Crespo da Silva, M.R.M., 1988b. Nonlinear flexural-flexural-torsional-extensional dynamics of beams – II. Response analysis. International Journal of Solids and Structures 24 (12), 1235–1242. Crespo da Silva, M.R.M., Glynn, C.C., 1978a. Nonlinear flexural-flexural-torsional dynamics of inextensional beams, I: Equations of motion. Journal of Structural Mechanics 6 (4), 437–448. Crespo da Silva, M.R.M., Glynn, C.C., 1978b. Nonlinear flexural-flexural-torsional dynamics of inextensional beams, II: Forced motions. Journal of Structural Mechanics 6 (4), 449–461. Cull, S.J., Tucker, R.W., Tung, R.S., Hartley, D.H., 2000. On parametrically excited flexural motion of an extensible and shearable rod with a heavy attachment. Technische Mechanik 20, 147–158. Esmailzadeh, E., Jalili, N., 1998. Parametric response of cantilever Timoshenko beams with tip mass under harmonic support motion. International Journal of Non-Linear Mechanics 33, 765–781. Femlab homepage: