A transfer matrix model of large deformations of curved rods

Computers and Structures 87 (2009) 467–484

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

A transfer matrix model of large deformations of curved rods Aviv Rosen *, Ohad Gur Technion – Israel Institute of Technology, Faculty of Aerospace Engineering, Technion City, Haifa 32000, Israel

a r t i c l e

i n f o

Article history: Received 27 April 2008 Accepted 29 December 2008 Available online 11 February 2009 Keywords: Transfer matrix Curved rods Pretwisted rods Large deformations

a b s t r a c t Transfer matrix models have been frequently used to analyze the structural behavior of rods, including curved and pretwisted ones. The advantages of these models include their relative simplicity, numerical efficiency and ease of implementation. Previous investigations did not include nonlinear analyses of curved rods that undergo large deformations. The present paper describes a nonlinear transfer matrix model of curved and pretwisted rods, which is capable of analyzing very large spatial deformations. The rod is divided into segments. A body system of coordinates is attached to each segment. This system translates and rotates with the segment during the deformation. If the segments are kept small enough, the local deformations of each segment, relative to its body system of coordinates, are small. The segments’ systems of coordinates rotate relative to their neighbors and if this rotation is dealt with properly, large rotations and displacement of the curved rod can be analyzed. In spite of its nonlinear nature, the model remains relatively simple and efficient. The new model is used to solve a few problems and the results show very good agreement with other analytical, finite-element, and experimental results. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction Curved rods are important structural elements that are commonly used for many applications. Because of their importance curved rods have been the subject of intensive research in the past. In many cases curved rods undergo large deformations that require a nonlinear analysis. The planar nonlinear behavior of curved rods was discussed in detail by Frisch-Fay [1]. Nevertheless, even the planar nonlinear behavior of curved rods continues to attract the interest of researchers, Refs. [2,3] are recent examples of such a research. The problem of the spatial deformations of a spatially curved and pretwisted rod is much more complicated than the planar one. Many researchers were concerned with the formulation of the nonlinear equations of equilibrium of spatially curved rods [4–7]. It was shown that in certain cases it is possible to obtain exact solutions for the out of plane problem of an arch [8]. Other researchers extended the equations to include composite curved rods [9,10]. Nonlinear models of curved rods were also used to investigate the stability of those structures [11–13]. Although various methods have been used to analyze large deformations of curved rod, for example Galerkin method [14], the common method of calculating the nonlinear behavior of curved rods involves the use of nonlinear finite elements. In * Corresponding author. Tel.: +972 4 8292711; fax: +972 4 8292030. E-mail addresses: [email protected] (A. Rosen), ohadg@aerodyne. technion.ac.il (O. Gur). 0045-7949/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2008.12.014

[15,16] the total Lagrangian formulation of a three-dimensional beam element was presented. Other formulations and discussions appear in [17–19]. A finite element implementation of the threedimensional finite-strain beam theory of Reissner was presented in [20]. A more recent investigation is described in [21] where thin-walled curved rectangular box beams, under torsion and out-of-plane bending, were considered. Other recent papers include finite-element formulations of geometrically exact threedimensional beam theories [22,23]. The past investigations included also dynamic analyses [24]. Of special interest and relevance to the present derivation is the co-rotational formulation that has been applied to the analysis of nonlinear behavior of structures [25,26]. In the case of finite element analysis this means the use of a frame of reference that is attached to the element and continuously rotates with it. For practical cases of small strains a regular linear formulation is used with respect to this rotating frame, while the rotation of this frame relative to a general fixed frame introduces the nonlinear effects. This method has found increasing use in the finite element nonlinear analysis of structures [27–30]. A weakness of the co-rotational formulation is its tendency to lead to a lack of objectivity of the strain measures [31,32]. This means a violation of the principle of invariance of strains to rigid body rotations. Special effort is required in order to maintain the property of objectivity. The purpose of the present paper is to present the application of the transfer matrix method for the nonlinear analysis of spatially curved rods. The transfer matrix method is one of the first numerical methods that was used to analyze structures [33]. It is beyond

468

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

Nomenclature a(i), b(i), c(i) components of the boundary displacements A(i), B(i), C(i), D(i) matrices associated with the ith segment the position vector of the ith nodal point before deford(i) mation 0 the position vector of the ith nodal point after deformadðiÞ tion dx(i), dy(i), dz(i) components of d(i) 0 0 0 0 dxðiÞ ; dyðiÞ ; dzðiÞ components of dðiÞ ex(i), ey(i), ez(i) a triad of orthogonal unit vectors in the direction, x(i), y(i), and z(i), respectively, before the deformation e0xðiÞ ; e0yðiÞ ; e0zðiÞ a triad of orthogonal unit vectors in the directions x(i), y(i) and z(i), respectively, after the deformation e00xðiÞ ; e00yðiÞ ; e00zðiÞ a triad of orthogonal unit vectors that are attached to the elastic axis and describe the elastic rotations along the segment matrix of elastic rotation along the ith segment E(i) (EIyy)(i), (EIzz)(i), (EIyz)(i) cross-sectional bending moments of inertia about the neutral axis of the ith segment F the magnitude of the applied force resultant cross-sectional force vector along the ith segF(i) ment F exðiÞ ; F eyðiÞ ; F ezðiÞ components of the resultant force that acts along the ith segment torsional stiffness of the rod, along the ith segment (GJ)(i) h(i), s(i), t(i) three parameters that describe the boundary rotation (i = 0, N + 1) i index of the ith segment kx, ky, kz linear spring constants boundary stiffness matrix (i = 0, N + 1) K(i) length vector of the ith segment l(i) coordinates of the circular rod l a, lb length of the ith segment L(i) resultant cross-sectional moment vector along the ith M(i) segment Mx(i), My(i), Mz(i) components of M(i) N number of nodes number of complex boundary conditions (i = 0, N + 1) N(i) distributed force per unit length, that acts along the ith p(i) segment px(i), py(i), pz(i) components of p(i) P(i) resultant cross-sectional tensile force along the ith segment

the scope of this paper to give a detailed survey of the use of this method in the past, only previous applications of the method to analyze curved rods will be discussed. Hayashi [34] represented a curved rod as a series of straight beam elements and used this representation to find the stiffness constants. He proposed a method to increase the accuracy of the transfer matrix calculations such that the numerical error would not increase as the number of partitions along the rod is increased. Irie et al. [35] used a transfer matrix approach to calculate the steady state in-plane response of a curved Timoshenko beam with internal damping, to a sinusoidally varying point force or moment. A method to analyze free and forced vibrations of one-dimensional structures was presented by Nitzsche [36]. His scheme linked the classical transfer matrix method and the integrating matrix technique. Haktanir and Kiral [37] investigated the static behavior of elastically and continuously supported helicoidal structures by the transfer and stiffness matrix methods. Refs. [35–37] presented investigations of various aspects of the structural behavior of planar curved rods by using transfer matrix

distributed moment vector, per unit length, that acts along the ith segment qx(i), qy(i), qz(i) components of q(i) Q exðiÞ ; Q eyðiÞ ; Q ezðiÞ components of the resultant moment that acts along the ith segment Position vector of the ith nodal point before deformation r(i) R radius Position vector of the ith nodal point after deformation R(i) rotation matrix between the systems (X, Y, Z) and S(i) (x(i), y(i), z(i)), before deformation rotation matrix between the systems (X, Y, Z) and S0ðiÞ (x(i), y(i), z(i)), after deformation rotation matrix, representing the initial curvature and T(i) pretwist at the ith node rotation matrix, representing the total curvature and T0ðiÞ pretwist of the deformed axis at the ith node u(i), v(i), w(i) components of d(i) utip, vtip, wtip components of the tip deflection of a rod relative to the system (X, Y, Z) Vy(i), Vz(i) components of the cross-sectional resultant shear force x(i), y(i), z(i) the body system of coordinates of the ith segment X, Y, Z a frame of reference that is fixed in space X coordinate of a cross-section at zero load X0 X coordinate of the rod tip at zero load X0-tip b Y; b Z b triad of unit vectors in the directions (X, Y, Z), respecX; tively cF angle of the applied force the local displacement vector of each point along the d(i) elastic axis of the ith segment hx(i), hy(i), hz(i) Elastic rotation angles along segment I rotation angle between the (i  1)th and the ith segf(i) ments, in the two dimensional case vector of boundary displacement and rotation k(i) l(i) vector of boundary conditions ~ðiÞ vector of complex boundary reactions l mx, my, mz components of angular spring constants /(i), u(i) elastic torsional rotation angle ~ tip projection of the tip rotation angle / v(i) state vector of the ith segment

q(i)

models. Kashimoto et al. [38] studied the dynamic stress concentration problem of an inhomogeneous rod of infinite length, consisting of two infinite straight portions and one finite curved portion of arbitrary curvature. Yildirim [39] investigated the inplane free vibrations of symmetric cross-ply laminated circular arches, using either Bernoulli–Euler or Timoshenko theories. Stephen and Ghosh [40] analyzed a curved repetitive pin-jointed structure using the state variable transfer matrix technique. The transfer matrix method is a very useful method to analyze curved rods because of its simplicity, numerical efficiency and ease of implementation. However, previous investigations were confined to small deformations. In what follows extension of the transfer matrix method for the case of large spatial deformations of curved rods will be presented. The rod will be described as a series of small segments. Due to the initial curvature and pretwist (of the undeformed rod), there are in general finite rotations between neighboring segments. Similar to the co-rotational formulation, the present derivation will also include the use of a body frame of reference that is attached to the two edges of each segment.

469

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

Thus for the case of small strains linear formulations for the deformations will be used. The nonlinearity will be introduced through the finite elastic rotations between neighboring segments. Due to the special formulation of the present derivations, objectivity is retained throughout the entire derivation. Although the present derivation will be confined to isotropic Bernoulli–Euler rods, the extension of the method to include other models (Timoshenko, thin walled, etc.) is straight forward. In spite of the extension of the method to include very large deformations, the relative simplicity and efficiency of the original method are retained. The implementation is straight forward, intuitive and simple. The solution includes mainly multiplication of small matrices the dimensions of which do not exceed 7 in the present examples. The dimension of the matrix that is inverted during each iteration is of same order. The accuracy of the results for the cross-sectional forces and moments is excellent even at large nonlinear deformations. Comparisons with the results of other methods exhibit excellent agreement. 2. Description of the rod geometry and loads Consider a spatially curved rod. The rod geometry is defined by its axis. r is the position vector of each point along the rod axis, which is given relative to a frame of reference (X, Y, Z) that is fixed in space. The rod’s cross-sections are perpendicular to the rod axis. In general the cross-sectional geometry and structural properties vary along the rod. The rod axis is divided into segments. Each segment is defined by two nodal points: origin and end points. The origin of segment i is defined by the vector r(i). Segment i ends at point r(i+1), which is also the origin of segment (i + 1), see Fig. 1. For each segment a local Cartesian system of coordinates is defined: x(i), y(i), z(i). x(i) is a coordinate line along the segment, connecting the origin and end points. y(i) and z(i) are normal to x(i). At the origin of the segment x(i) = 0, at the end of the segment x(i) = L(i). ex(i), ey(i) and ez(i) form a triad of unit vectors in the directions of the coordinate lines: x(i), y(i) and z(i), respectively. Since a spatially curved and pretwisted rod is considered, there are in general rotations between consequent coordinates systems, thus:

eðiÞ ¼ TðiÞ  eði1Þ ;

ð1Þ

eTðiÞ

ð2Þ

¼ ½exðiÞ ; eyðiÞ ; ezðiÞ :

T(i) is the rotation matrix between the (i) and (i  1) systems of coordinates, and it is a function of the elements’ curvature and pretwist.

y (i )

z (i )

It is possible now to continue the derivation, taking into account the exact curvature and pretwist along each segment. This introduces a significant complexity into the derivations, as compared to a straight segment. Since the solution procedure includes a discretization scheme where the: geometry, loads and structural properties – are approximated along each segment, it looks appropriate to replace each curved rod segment by a straight segment, connecting points r(i) and r(i+1). The geometric and structural properties of the straight segment are identical to those of the original curved one. Thus, the influences of curvature and pretwist are solely introduced by the rotation matrices, T(i). The decision to replace the curved segment by a straight one is also based on the goal of the present derivation, to develop a model that is: simple, efficient and easy to implement. It is easier and more efficient to increase the number of segments along the rod, than to introduce the local curvature into the segment’s derivations. Furthermore, in many cases curved rods are assembled from straight segments. In other cases the rod geometry is described by the coordinates of discrete points along its axis, without any exact information about the geometry between those points. The results in what follows will justify the approximation of straight segments. The consequence of adopting this approximation is that the segments should be chosen short enough, such that the maximum distance of the curved elastic axis from the coordinate xi is very small compared to L(i). There are N segments along the rod. The matrices T(i) and scalars L(i) define the geometry of the elastic axis of the undeformed rod, as explained in Appendix A. The rod deforms due to the loads that act along it. After deformation the position vector of each nodal point along the rod axis becomes R(i). The (x(i), y(i), z(i)) system of coordinates is rigidly translated and rotated such that its origin is located at the point R(i) and the axis x(i) passes through the point R(i+1). e0xðiÞ ; e0yðiÞ and e0zðiÞ is a triad of orthogonal unit vectors in the directions x(i), y(i) and z(i), respectively, after the deformation. The distributed force and moment per unit length, that act along segment (i), p(i)(x(i)) and q(i)(x(i)), are described by their components, as follows:

pðiÞ ðxðiÞ Þ ¼ pxðiÞ ðxðiÞ Þ  e0xðiÞ þ pyðiÞ ðxðiÞ Þ  e0yðiÞ þ pzðiÞ ðxðiÞ Þ  e0zðiÞ ;

ð3Þ

qðiÞ ðxðiÞ Þ ¼ qxðiÞ ðxðiÞ Þ  e0xðiÞ þ qyðiÞ ðxðiÞ Þ  e0yðiÞ þ qzðiÞ ðxðiÞ Þ  e0zðiÞ :

ð4Þ

The displacement vector of each point along the elastic axis of element i, relative to the system (x(i), y(i), z(i)), is d(i)(x(i)):

dðiÞ ðxðiÞ Þ ¼ uðiÞ ðxðiÞ Þ  e0xðiÞ þ v i ðxðiÞ Þ  e0yðiÞ þ wðiÞ ðxðiÞ Þ  e0zðiÞ :

It should be noted that d(i)(x(i)) is not the global displacement of the points along the axis of segment (i). In order to obtain the global

Nodal point ( i + 1) Segment i

x( i )

Nodal point ( i )

z( i ) Nodal point ( i )

The axis of the deformed rod

R (i )

ð5Þ

y (i ) L (i )

R ( i +1) r( i )

r ( i +1)

Nodal point ( i + 1) x( i )

The axis of the undeformed rod

Z Y

X Fig. 1. The rod geometry before and after deformation.

470

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

displacement, the translations and rotations of the systems (x(i), y(i), z(i)) should be taken into account, as it is shown in Appendix A. Based on the definition of the system of coordinates of the segment after the deformation:

uðiÞ ð0Þ ¼ v ðiÞ ð0Þ ¼ wðiÞ ð0Þ ¼ v ðiÞ ðLðiÞ Þ ¼ wi ðLðiÞ Þ ¼ 0:

ð6Þ

Since most of the engineering materials are confined to small strains, it is also assumed that:

LðiÞ þ uðiÞ ðLðiÞ Þ ffi LðiÞ :

ð7Þ

In addition to displacements, there is also an elastic torsional rotation about the axis of the segment. This elastic torsional rotation is described by the angle u(i)(x(i)). Without a loss of generality it is convenient to choose:

uðiÞ ð0Þ ¼ 0:

ð8Þ

In the present investigation the segments’ lateral displacement components and elastic torsional rotation, are approximated as follows:

    dv dv dv   1 v ðiÞ xðiÞ ¼ ðiÞ ð0Þ  xðiÞ þ  2 ðiÞ ð0Þ  ðiÞ LðiÞ  x2ðiÞ LðiÞ dxðiÞ dxðiÞ dxðiÞ   dv ðiÞ dv ðiÞ   1 þ 2  ð0Þ þ LðiÞ  x3ðiÞ ; dxðiÞ LðiÞ dxðiÞ     dwðiÞ dwðiÞ dwðiÞ   1 ð0Þ  xðiÞ þ 2 ð0Þ  LðiÞ  x2ðiÞ wðiÞ xðiÞ ¼ LðiÞ dxðiÞ dxðiÞ dxðiÞ   dwðiÞ   1 dwðiÞ þ 2 ð0Þ þ LðiÞ  x3ðiÞ ; dxðiÞ LðiÞ dxðiÞ 



uðiÞ xðiÞ ¼

duðiÞ dxðiÞ

ð9aÞ

ð9bÞ

" # duðiÞ   duðiÞ 1 ð0Þ  xðiÞ þ  LðiÞ  ð0Þ  x2ðiÞ : 2  LðiÞ dxðiÞ dxðiÞ

segment i, to the axis x(i), is very small compared with L(i). In the same manner, although the global deformations are large, it is assumed that the segments are small enough such that:

   v ðiÞ xðiÞ  ;    L ðiÞ

ð10Þ

The local elastic rotations are also assumed to be small (although the global elastic rotations can be very large):

  dv ðiÞ     dx ; ðiÞ

  dwðiÞ     dx ; ðiÞ

    uðiÞ   1:

ð11Þ

Consider now a triad of orthogonal unit vectors e00xðiÞ ðxðiÞ Þ; e00yðiÞ ðxðiÞ Þ, and e00zðiÞ ðxðiÞ Þ that is attached to the elastic axis and describes the elastic rotations along the segment:

      e00ðiÞ xðiÞ ¼ EðiÞ xðiÞ  e0ðiÞ xðiÞ ;

ð12Þ

h   iT   ¼ ½e00xðiÞ xðiÞ ; e00yðiÞ ðxðiÞ Þ; e00zðiÞ ðxðiÞ Þ; e00ðiÞ xðiÞ h   iT h      i ¼ e0xðiÞ xðiÞ ; e0yðiÞ xðiÞ ; e0zðiÞ xðiÞ : e0ðiÞ xðiÞ

2

Eqs. (9a–c) satisfy Eqs. (6) and (8). As indicated above, the segments are chosen small enough, such that the distance from the ‘‘real” curved rod axis of the undeformed

6 ¼4

TðiÞ

cos fðiÞ

0  sin fðiÞ

0

1

0

sin fðiÞ

0

cos fðiÞ

z( i −1)

The elastic axis before deformation

7 5:

ð14Þ

(

)

z( i ) z( i )

x(i )

Before deformation

ð13bÞ

3

The elastic axis after deformation

θ (i −1) L(i −1)

ζ (i )

ð13aÞ

E(i)(x(i)) is the matrix of elastic rotation along segment i. The expression for E(i)(x(i)) is given in Appendix B. As indicated above, the axis of the initially curved and pretwisted undeformed rod has been replaced by a chain of straight rod segments. The initial curvature and pretwist are presented by the rotation matrices, T(i), between segments (i) and (i  1). In a similar manner, the elastic curvature and twist are also presented by elastic rotations between the straight segments, which are superimposed on the initial rotations. In order to understand this description, the two-dimensional case is presented in Fig. 2. In the figure, segments (i  1) and (i), before and after the deformation, are shown. Before the deformation, because of the initial curvature, there is an angle f(i) between coordinate lines x(i1) and x(i). Thus the rotation matrix, T(i) (see Eq. (1)) is:

ð9cÞ

z(i −1)

   wðiÞ xðiÞ    1:    L ðiÞ

x(i −1)

−θ (i ) ( 0 )

x( i )

x( i −1)

(

)

ζ (i ) + θ (i −1) L(i −1) − θ (i ) ( 0 ) After deformation

Fig. 2. The rotation between the systems of coordinates of the (i  1) and (i) segments – a two dimensional case.

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

Because of the deformation there is an elastic rotation angle h(i1)(x(i1)) along segment (i  1), and an elastic rotation angle h(i)(x(i)) along segment i, thus:

   3 cos hðiÞ ðxðiÞ Þ 0  sin hðiÞ ðxðiÞ Þ 6 7 EðiÞ ðxðiÞ Þ ¼ 4 0 1 0 5:     sin hðiÞ ðxðiÞ Þ 0 cos hðiÞ ðxðiÞ Þ 2

ð15Þ

The elastic rotations are defined in more details in Appendix B. Similar to Eq. (1), it is possible to write:

e0ðiÞ

¼

T0ðiÞ



e0ði1Þ ;

where the rotation matrix after the deformation,

T0ðiÞ

F exðiÞ ¼ F eyðiÞ ¼ F ezðiÞ ¼ Q exðiÞ ¼

ð16Þ T0ðiÞ ,

Q eyðiÞ ¼

is:

¼ EðiÞ ð0Þ  TðiÞ  Eði1Þ ðLði1Þ Þ:

ð17Þ

Eqs. (16) and (17) are also applicable in the three-dimensional case, where the elastic rotation matrices and the rotation matrices representing the initial curvature and pretwist, are the appropriate threedimensional matrices according to Eq. (1) and Appendix B. The matrices T0ðiÞ define the geometry of the elastic axis after deformation, as explained in Appendix A. It should be noted that due to the present formulation rigid body rotations or translations cannot result in any strains. Thus by definition the present formulation fulfills the objectivity criterion.

Q ezðiÞ ¼

Z

LðiÞ

0

Z Z

LðiÞ

0 LðiÞ

471

  pxðiÞ xðiÞ  dxðiÞ ;

ð22aÞ

  pyðiÞ xðiÞ  dxðiÞ ;

ð22bÞ

  pzðiÞ xðiÞ  dxðiÞ ;

ð22cÞ

  qxðiÞ xðiÞ  dxðiÞ ;

ð23aÞ

  qyðiÞ xðiÞ  dxðiÞ ;

ð23bÞ

  qzðiÞ xðiÞ  dxðiÞ :

ð23cÞ

0

Z

LðiÞ

0

Z

LðiÞ

0

Z 0

LðiÞ

Integration of Eqs. (19a–c) along the segment, lead to the following algebraic equations:

  PðiÞ LðiÞ  PðiÞ ð0Þ þ F exðiÞ ¼ 0;   V yðiÞ LðiÞ  V yðiÞ ð0Þ þ F eyðiÞ ¼ 0;   V zðiÞ LðiÞ  V zðiÞ ð0Þ þ F ezðiÞ ¼ 0:

ð24aÞ ð24bÞ ð24cÞ

3. The equilibrium equations of each segment

Integration of Eqs. (21a–c), assuming that Vy(i)(x(i)) and Vz(i)(x(i)) vary in a linear manner along the segment, result in the following algebraic equations:

As indicated above by Eqs. (10) and (11), the local deformations of the segment, namely the displacements and rotations relative to the local body system or coordinates, are very small. Thus the local elastic equations can be the linear equations of rods. The nonlinearity will be presented by the matrices T0ðiÞ . The resultant cross-sectional force along segment i, is:

  MxðIÞ LðiÞ  M xðiÞ ð0Þ þ Q exðiÞ ¼ 0;     V zðiÞ LðiÞ þ V zðiÞ ð0Þ MyðiÞ LðiÞ  M yðiÞ ð0Þ   LðiÞ þ Q eyðiÞ ¼ 0; 2     V yðiÞ LðiÞ þ V yðiÞ ð0Þ MzðiÞ LðiÞ  M zðiÞ ð0Þ þ  LðiÞ þ Q ezðiÞ ¼ 0: 2

FðiÞ ¼ P ðiÞ ðxðiÞ Þ  e0xðiÞ þ V yðiÞ ðxðiÞ Þ  e0yðiÞ þ V zðiÞ ðxðiÞ Þ  e0zðiÞ :

ð18Þ

P(i) is the resultant cross-sectional tensile force, while Vy(i) and Vz(i) are the components of the cross-sectional resultant shear force. Thus the three scalar equations of equilibrium of forces along segment (i) become (see Eq. (3)):

  dPðiÞ   xðiÞ þ pxðiÞ xðiÞ ¼ 0; dxðiÞ

ð19aÞ

  dV yðiÞ   xðiÞ þ pyðiÞ xðiÞ ¼ 0; dxðiÞ

ð19bÞ

  dV zðiÞ   xðiÞ þ pzðiÞ xðiÞ ¼ 0: dxðiÞ

ð19cÞ

The resultant cross-sectional moment about the elastic axis, along segment i, is:

      MðiÞ ¼ M xðiÞ xðiÞ  e0xðiÞ þ M yðiÞ xðiÞ  e0yðiÞ þ M zðiÞ xðiÞ  e0zðiÞ :

ð20Þ

The three scalar equations of equilibrium of moments, along segment (i), become (see Eq. (4):

dMxðiÞ   xðiÞ þ qxðiÞ ¼ 0; dxðiÞ

ð21aÞ

    dMyðiÞ   xðiÞ  V zðiÞ xðiÞ þ qyðiÞ xðiÞ ¼ 0; dxðiÞ

ð21bÞ

    dMzðiÞ   xðiÞ þ V yðiÞ xðiÞ þ qzðiÞ xðiÞ ¼ 0: dxðiÞ

ð21cÞ

Now, resultant applied force and moment components that act on each segment, will be defined:

ð25aÞ ð25bÞ ð25cÞ

4. Calculation of the deformation In order to calculate the deformation along each segment, a set of equations are defined describing the relations between the segment deformations and the resultant cross-sectional forces and moments. In the present case the following assumptions are adopted: a. The rod is made of isotropic elastic material that has identical material properties at any point of the rod segment. b. Bernoulli–Euler bending, namely: cross-sections that are normal to the rod axis before the deformation, remain normal to the same axis after deformation. c. Saint-Venant torsion. d. The rod axis connects the shear-centers of all the cross-section. This means that the resultant cross-sectional shear force does not result in a torsional moment about the rod axis. Based on the above assumptions, the following three equations are obtained:

duðiÞ     MxðiÞ xðiÞ ¼ ðGJÞðiÞ  xðiÞ ; dxðiÞ

ð26aÞ

      d2 v ðiÞ   xðiÞ MyðiÞ xðiÞ ¼ PðiÞ xðiÞ  zcðiÞ  EIyz ðiÞ  2 dxðiÞ   d2 wðiÞ   xðiÞ ;  EIyy ðiÞ  2 dxðiÞ

ð26bÞ

472

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

2     d v ðiÞ   M zðiÞ xðiÞ ¼ PðiÞ xðiÞ  ycðiÞ þ ðEIzz ÞðiÞ  xðiÞ 2 dxðiÞ



þ EIyz

 ðiÞ

2



d wðiÞ  2

dxðiÞ



xðiÞ :

From Eqs. (9a–c) the following relations are obtained:

ð26cÞ

(GJ)(i) is the torsional stiffness of the rod along segment i. In the case of a pretwisted rod the increase of the torsional stiffness due to pretwist [41] is introduced into (GJ)(i). The formulation can also be quite easily extended to model non-uniform torsion along the segment, including the influence of warping variations along it. yc(i) and zc(i) are the cross-sectional coordinates of the tension center (neutral axis). (EIyy)(i), (EIzz)(i) and (EIyz)(i) are the cross-sectional bending moments of inertia, about the tension center. It should be pointed out that Eqs. (26a–c) can easily be replaced by the appropriate equations of any other model of bending and torsion, including: Timoshenko bending, anisotropic materials, thin cross-section models, etc. The nonlinear, large deformation analysis, which is the subject matter of the present paper, will remain unchanged. Eq. (26a) can be written at the segment origin, x(i) = 0, and end point, x(i) = L(i), thus:

duðiÞ dxðiÞ

ð0Þ ¼

M xðiÞ ð0Þ ; ðGJÞðiÞ

duðiÞ  dxðiÞ

   M xðiÞ LðiÞ LðiÞ ¼ : ðGJÞðiÞ

ð27Þ

Eqs. (26b–c) are also written at the segment origin, x(i) = 0, and the segment end point, x(i) = L(i), thus the following relations are obtained:

2

d

v ðiÞ 2

dxðiÞ

2

d

2



LðiÞ ¼

2

d wðiÞ 2

dxðiÞ

2

dxðiÞ

2

ð0Þ ¼

ðLðiÞ Þ ¼

2

2

ð29cÞ

2

ð29dÞ

2

dwðiÞ LðiÞ d wðiÞ LðiÞ d wðiÞ ðLðiÞ Þ ¼ ð0Þ þ ðLðiÞ Þ:   dxðiÞ 6 dx2ðiÞ 3 dx2ðiÞ

ð29eÞ

Substituting Eqs. (27) and (28a–c) into Eqs. (29a–e) enables the dv calculation of the segment’s deformations, uðiÞ ðLðiÞ Þ; dxðiÞðiÞ ð0Þ; dv ðiÞ dwðiÞ dwðiÞ ðL Þ; ð0Þ; ðL Þ, as functions of the resultant forces and ðiÞ ðiÞ dxðiÞ dxðiÞ dxðiÞ moments at the origin and end points of the segment. As it will be described in the next section, the resultant forces and moments are obtained after applying the transfer matrix formulation and then the deformation calculations are carried out based on the transfer matrix solution.

ð28aÞ

ðEIyy ÞðiÞ  ðEIzz ÞðiÞ  ðEIyz Þ2ðiÞ

;

ð28bÞ

ð28cÞ

;

h i ðEIzz ÞðiÞ  ðM yðiÞ ðLðiÞ Þ þ PðiÞ ðLðiÞ Þ  zcðiÞ Þ  ðEIyz ÞðiÞ  MzðiÞ ðLðiÞ Þ þ PðiÞ ðLðiÞ Þ  ycðiÞ ðEIyy ÞðiÞ  ðEIzz ÞðiÞ  ðEIyz Þ2ðiÞ

ð29bÞ

2

dwðiÞ LðiÞ d wðiÞ LðiÞ d wðiÞ ð0Þ ¼  ð0Þ  ðLðiÞ Þ;   dxðiÞ 3 dx2ðiÞ 6 dx2ðiÞ

h i             ðEIyy ÞðiÞ  M zðiÞ LðiÞ þ PðiÞ LðiÞ  ycðiÞ  EIyz ðiÞ  M yðiÞ LðiÞ þ PðiÞ LðiÞ  zcðiÞ

ðEIyy ÞðiÞ  ðEIzz ÞðiÞ  ðEIyz Þ2ðiÞ

ð29aÞ

2

dv ðiÞ LðiÞ d v ðiÞ LðiÞ d v ðiÞ ðLðiÞ Þ ¼ ð0Þ þ ðLðiÞ Þ;   dxðiÞ 6 dx2ðiÞ 3 dx2ðiÞ

h i   ðEIzz ÞðiÞ  M yðiÞ ð0Þ þ PðiÞ ð0Þ  zcðiÞ  ðEIyz ÞðiÞ  M zðiÞ ð0Þ þ PðiÞ ð0Þ  ycðiÞ

2

d wðiÞ

2

dv ðiÞ LðiÞ d v ðiÞ LðiÞ d v ðiÞ ð0Þ ¼  ð0Þ  ðLðiÞ Þ;   dxðiÞ 3 dx2ðiÞ 6 dx2ðiÞ

h i     ðEIyy ÞðiÞ  MzðiÞ ð0Þ þ P ðiÞ ð0Þ  ycðiÞ  EIyz ðiÞ  M yðiÞ ð0Þ þ PðiÞ ð0Þ  zcðiÞ ð0Þ ¼ ;  2 ðEIyy ÞðiÞ  ðEIzz ÞðiÞ  EIyz ðiÞ

v ðiÞ 

dxðiÞ

uðiÞ ð0Þ ¼ 0;

" # duðiÞ LðiÞ duðiÞ uðiÞ ðLðiÞ Þ ¼  ð0Þ þ ðLðiÞ Þ ; 2 dxðiÞ dxðiÞ

:

ð28dÞ

473

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

F

Tip projection after deformation

F

Y

F

Y vtip F

t1

tip

X

Z

x

t2

Z

wtip Fig. 3. The spatial deformation of a straight cantilevered rod as a result of a tip force.

12

9

Test 8

[42]

γ F = 40 [ deg ]

Linear 10

Non Linear 7

γ F = 30 [ deg ] 8

6 5

wtip ["]

6

4

wtip ["] 4

3 2

vtip ["]

2

vtip ["]

1 0 0

1

2

3

0

4 F [lbf .] 5

0

1

2

4 F [lbf .] 5

3

Fig. 4. Nonlinear deformations of a cantilevered rod – calculations and test results [42] of the tip deflection in the case of a 1/200  1/800 cross-section.

4

7

Te st

γ F = 30 [ deg ]

3.5

[42]

Linear

6

Non Linear 3

γ F = 50 [ deg ]

5 2.5

wtip ["]

4

wtip ["]

2 3 1.5 2 1 1

0.5

vtip ["]

0 0

1

2

3

4

vtip ["] F [lbf .]

0 5

0

1

2

3

4

F [lbf .] 00

5 00

Fig. 5. Nonlinear deformations of a cantilevered rod – calculations and test results [42] of the tip deflection in the case of 1  1/8 cross-section.

474

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

n o A1ðiÞ ¼ 1; 0; 0; 0; 0; 0; F exðiÞ ;

5. The transfer matrix formulation

n o A2ðiÞ ¼ 0; 1; 0; 0; 0; 0; F eyðiÞ ;

A state vector of order seven, v(i)(x(i)), is defined as follows:

n o A3ðiÞ ¼ 0; 0; 1; 0; 0; 0; F ezðiÞ ;

vTðiÞ ðxðiÞ Þ ¼ ½PðiÞ ðxðiÞ Þ; V yðiÞ ðxðiÞ Þ; V zðiÞ ðxðiÞ Þ; M xðiÞ ðxðiÞ Þ; MyðiÞ ðxðiÞ Þ; M zðiÞ ðxðiÞ Þ; 1: ð30Þ According to Eqs. (24a–c) and (25a–c):

n o A4ðiÞ ¼ 0; 0; 0; 1; 0; 0; Q exðiÞ ; A5ðiÞ

vðiÞ ðLðiÞ Þ ¼ AðiÞ  vðiÞ ð0Þ:

ð31Þ

A(i) is square matrix of order seven. The rows of this matrix are defined below:

7

A6ðiÞ

ð32Þ



LðiÞ ¼ 0; 0; LðiÞ ; 0; 1; 0; F ezðiÞ   Q eyðiÞ ; 2

LðiÞ ¼ 0; LðiÞ ; 0; 0; 0; 1; F eyðiÞ   Q ezðiÞ ; 2

A7ðiÞ ¼ f0; 0; 0; 0; 0; 0; 1g:

10

φtip [ deg ]

Te st [42]

φtip [ deg ]

9

Non-Linear

6

Te st

[42]

Non-Linear

8

5

γ F = 15 [ deg ]

γ F = 30 [ deg ]

7 6

4

5

3

4 3

2

2

1 1

F [lbf .]

0 0

1

2

3

4

5

0

1.8

4

φtip [ deg ]

3

1

1.6

3

4

5

γ F = 60 [ deg ]

1.4

γ F = 45 [ deg ]

2

Te st [42] Non-Linear

φtip [ deg ]

Test [42] Non-Linear

3.5

F [lbf .]

0

1.2

2.5 1

2 0.8

1.5

0.6

1

0.4

0.5

0.2

F [lbf .]

F [lbf .]

0

0 0

1

2

3

4

0

0.5

1

1.5

2

2.5

Fig. 6. Nonlinear deformations of a cantilevered rod – calculations of the projection of the tip rotation onto the Y–Z plane, case of a 1/200  1/800 cross-section.

475

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

14

8 tip

[33]

deg

Te st 50 Nodes 25 Nodes 10 Nodes 5 Nodes

7

6

F

30 deg

tip

deg

Te st

[33]

F

30 deg

50 Nodes

12

25 Nodes 10 Nodes 10

5 Nodes

5

8 wtip "

4

6 3

4 2

vtip "

2

1

F lbf .

F lbf .

0

0 0

1

2

3

4

5

0

1

2

3

4

5

Fig. 7. Nonlinear deformations of a cantilevered rod – the influence of the number of nodes on the tip deflection in the case of a 1/200  1/800 cross-section and cF = 30°.

The relation between the state vectors on both sides of nodal point (i), namely the joint between the end of segment (i  1) and the origin of segment (i), is given by:

Eqs. (31) and (33) lead to the following relation between the state vectors at the tip of two consequent segments:

vðiÞ ð0Þ ¼ BðiÞ  vði1Þ ðLði1Þ Þ; 2 0 3 TðiÞ 6 7 BðiÞ ¼ 4 5: T0ðiÞ

where:

ð33Þ

0.14

i

i ¼ 0; ðN þ 1Þ:

ð37Þ

It is clear that:

0.1

, dx i

dv i

TðNþ1Þ ¼ Eð0Þ ðLð0Þ Þ ¼ EðNþ1Þ ð0Þ ¼ I33 ;

ð38aÞ

AðNþ1Þ ¼ I77 :

ð38bÞ

The relation between the vectors of boundary conditions is:

0.08

lðNþ1Þ ¼ D  lð0Þ ;

dx i

,

ð36Þ

Special attention should be given to the two edges of the rod, where the boundary conditions are applied. In order to do that two additional ‘‘virtual” segments of zero length are added to the rod at its two edges. The axes of these segments are tangent to the axis of the deformed rod, at those edges. These are segments (i = 0) and (i = N + 1). A vector of boundary reactions, l(i), is defined in a similar manner to v(i):

  lTðiÞ ¼ PðiÞ ; V yðiÞ ; V zðiÞ ; M xðiÞ ; M yðiÞ ; M zðiÞ ; 1 ;

5 Nodes 10 Nodes 25 Nodes

0.12

dw i

ð35Þ

CðiÞ ¼ AðiÞ  BðiÞ : ð34Þ

1

max

vðiÞ ðLðiÞ Þ ¼ CðiÞ  vði1Þ ðLði1Þ Þ;

ð39Þ

where:

0.06

D ¼ CðNþ1Þ  CðNÞ  CðN1Þ    Cð2Þ  Cð1Þ :

0.04

0.02

F lbf . 0 0

1

2

3

4

5

Fig. 8. Nonlinear deformations of a cantilevered rod – the maximum value of the ‘‘local” rotation components as a function of the load magnitude and number of nodes, in the case of a 1/200  1/800 cross-section and cF = 30°.

ð40Þ

There are 12 parameters that define the boundary conditions of forces and moments. Six of those parameters are assumed to be known, while the other six are obtained by solving Eq. (39). In the case of statically determinate structures, six of the force or moment reactions at the edges, are indeed known. Thus, for example, in the case of clamped/free boundary conditions, where the rod is loaded by a concentrate force at the free tip, (see examples (a) and (c) of the next section), all the force and moment components at the free edge are known. The approach in the case of statically indeterminate cases is discussed in Appendix C. Eq. (39) presents a system of seven equations, of which one is trivial (1 = 1). The solution procedure is as follows:

476

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

b. The resultant applied forces and moments that act on each segment are calculated according to Eqs. (22a–c) and (23a–c). In the case where these resultants depend on the

a. The undeformed rod is divided into N segments. The structural properties of each segment are defined and the matrices T(i) and segment lengths L(i) are obtained.

Z F

R 2F

2F

R

la lb

la

R

X

lb Fig. 9. Planar bending of a circular rod. The solid lines present the geometry of the undeformed rod while dashed lines exhibit the geometry after deformation.

10

F ⋅ R2 EI

8

lb R

la R

6 4 Analysis Ref.[44] 2

Linear Non-Linear

0 -2 -4

lb R

-6 la R

-8 -10 -0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

Fig. 10. Planar deformation of a circular rod – components of the normalized displacements at the middle of the rod as functions of the normalized applied force.

1.6

Z R

[44]

Analysis Ref.[44]

F ⋅ R2 = −3 EI

1.4

Non-Linear

1.2 F ⋅ R2 =0 EI

1

F ⋅ R2 =3 EI

0.8 0.6 0.4

F ⋅ R2 = −9 EI

F ⋅ R2 =9 EI

0.2

X R

0 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Fig. 11. Planar deformation of a circular rod – comparison with Ref. [44] of the geometry after the deformation.

477

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

0.25

5 Nodes 10 Nodes 25 Nodes 0.2

Y

, φ( i )

Z

X

0.15

dx( i )

i.

F

dx( i )

g. h.

It should be noted that the present formulation do not lead to the use of the tangent stiffness matrix that plays a major roll in other formulations (for example see [19,25,26]). In all the cases that have been analyzed by the present method, including the

dv( i )

f.

,

e.

dw( i )

d.

j. The segments’ deformations are calculated according to Section 4. k. Convergence of the solution is checked. If convergence has not been reached, the elastic rotations along the rod are updated and the iterative procedure starts again from (b). l. In case of convergence the geometry of the deformed rod is calculated according to the derivations of Appendix A and the appropriate boundary conditions.

0.1

max

c.

rod deformation, they are updated at the beginning of each iteration. The matrices A(i), E(i)(0) and E(i)(L(i)), for 1 6 i 6 N, are calculated using Eq. (32) and the appropriate expressions of Appendix B. The matrices T0ðiÞ , for 1 6 i 6 (N + 1), are calculated according to Eq. (17), using also Eq. (38a). Matrices B(i) for 1 6 i 6 (N + 1) are assembled according to Eq. (34). Matrices C(i) for 1 6 i 6 (N + 1) are calculated according to Eq. (36). Matrix D is calculated according to Eq. (40). The six known boundary reactions are substituted into Eq. (39) and the system of equations is solved in order to calculate the six unknown boundary parameters. Using the transfer matrix formulation, the resultant forces and moment along the rod are calculated.

π 4

R

0.05

F ⋅ R2 EI

0 0

4

6

[15]

FE Analysis Linear Non-Linear, 50 Nodes Non-Linear, 25 Nodes Non-Linear, 10 Nodes Non-Linear, 5 Nodes

0.5

0.4

vtip R

wtip

0.3

utip

utip

vtip

R 0.2

0.1

wtip R F ⋅ R2 EI

0 0

1

2

8

Fig. 14. Spatial bending of a circular rod – the maximum value of the ‘‘local” rotation components as a function of the load magnitude and number of nodes.

Fig. 12. Spatial bending of a circular rod.

0.6

2

3

4

5

6

7

8

Fig. 13. Spatial bending of a circular rod – the components of the tip displacement, as functions of the tip force and the number of nodes along the rod.

478

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

cases that are presented in the next section, convergence was very fast. In fact incremental loading was not necessary for convergence. Even at highly nonlinear cases at high loads, convergence was obtained within a few iterations while starting from the undeformed geometry and applying the actual full load from the start.

b. Planar bending of a circular rod that undergoes very large deformations. The results will be compared with analytical results from the literature. c. Spatial bending of a circular rod. The results will be compared with finite-element results from the literature.

6. Results 6.1. Spatial bending of a straight cantilevered rod The above described model will be used to calculate the deformations in the following cases: a. Spatial bending of a straight cantilevered rod. This example will help in verifying the capability of the model to analyze large deformations. The results will be compared with experimental results from the literature.

8

The purpose of this example is to present the model’s capability of calculating large nonlinear spatial deformations. Ref. [42] presents test results for the spatial bending of cantilevered rods under the action of a concentrated tip force. Two 2000 long rods, with rectangular cross-sections, were tested. The cross-sectional dimensions are 1/200  1/800 and 100  1/800 . The rod is made of Aluminum

8 Integration Transfer Matrix

P ⋅ R2 7 E⋅I

Integration Transfer Matrix

Vy ⋅ R2 E ⋅7I

F ⋅ R2 =8 EI

F ⋅ R2 =8 EI

6

6

5

6

6

5 4 4 3

4

4

3 2 2

2

X0 0.75 X 0−tip 1

1

1 0 0

0.25

0.5

2

0

0.25

0.5

X0 0.75 X 0−tip 1

2 Integration Transfer Matrix

Vz ⋅ R2 E1.5 ⋅I

1

0.5

0

2 4

-0.5 6

-1 F ⋅ R2 =8 EI

-1.5 0

0.25

0.5

0.75

X0 X 0−tip

1

Fig. 15. Spatial bending of a circular rod – distributions along the rod of the cross-sectional resultant force components.

479

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

7075. The force, F, which is applied at the free tip, is normal to the rod’s axis (before deformation) and forms an angle cF with the ‘‘long” cross-sectional axis of symmetry (see Fig. 3). When c – 0 the tip is deflected and rotates (Fig. 3). For a rectangular cross-section of dimensions t1 and t2, the torsional stiffness, (GJ)(i), is given by the following equation [43]:

ðGJÞðiÞ ¼ k1  G  t31  t2

for t 1 < t 2 ;

k1 ¼ 0:281 for t 2 =t1 ¼ 4;

ð41Þ

k1 ¼ 0:304 for t 2 =t 1 ¼ 6: The results that will be presented include: test results, results of the linear theory and converged results of the present nonlinear model. Fig. 4 presents the tip deflections of the rod having a 1/200  1/800 cross-section, for two different force angles, cF, as functions of the force magnitude, F. The linear model gives good result up to F = 2[lb.]. For higher loads the agreement deteriorates. The flatwise deflection, wtip, is over predicted by the linear model, while the edgewise deflection, vtip, is under predicted. The nonlinear model exhibits a very good agreement with the test results. Fig. 5 presents similar results for the second rod (with a crosssection of 100  1/800 ). In this case the differences between the linear and nonlinear models are smaller, because the rod is stiffer and the deflections are smaller. Yet, the trends are similar to those of Fig. 4. Fig. 6 presents a projection, onto the Y–Z plane, of the tip rotation. It should be noted that this angle is solely due to nonlinear effects and equals to zero in the case of the linear calculations. The comparison is shown for the rod having a cross-section 1/2  1/8, at four different force angles. In all the cases the agreement between the results of the nonlinear model and the test results is very good. In Fig. 7 the influence of the number of nodes on the convergence of the results is presented for the case of a 1/200  1/800 cross-section and a load angle of 30°. All the segments have the same length. The differences between the various number of nodes increase as the load is increased. It can be concluded that for all the loads 10 nodes give tip deflections that are within 2.7% of the converged results. In the case of the tip rotation, which represents a solely nonlinear effect, for a load of 5[lbf.], 10 nodes give errors of 14%. The error is reduced to 3.4% for 25 nodes. As indicated by Eq. (11), the derivation is based on the assumption that the local elastic rotations are small. In Fig. 8 the maximum local elastic rotation component along the rod is shown as a function of the tip load and number of nodes. As expected these values increase as the tip load is increased (and thus the deformations increase) or as the result of reducing the number of nodes. If the maximum elastic rotation component is smaller than 0.05, the error of the tip deflection calculation does not exceed 5%. In the case of the tip rotation, in order to obtain an error smaller than 0.1° it is necessary to have local elastic rotation components smaller than 0.02.

of the force (X direction). This is an indeterminate case that is solved by the procedure described in Appendix C. The tip coordinates, in the direction of the force and perpendicular to it, are la and lb, respectively (see Fig. 9). The problem was solved analytically by Bunce and Brown [44]. Fig. 10 presents a comparison between linear results, nonlinear results of the present transfer matrix model and analytical results from Bunce and Brown [44]. The normalized tip coordinates (dividing the dimensional coordinates by the rod radius R), are presented as functions of the normalized load, (F  R2/E  I). The linear model predicts only the X component of the displacement, while the normal displacement component is always equal to zero. The results of the linear model for the X component of the displacement exhibit a fair accuracy only for low loads (|F  R2/E  I|  1), while for higher loads the inaccuracy increases. The nonlinear model agrees very well with the analytical results, for both compression and tension, even for extremely high deformations. Fig. 11 presents the shape of the rod after deformation, for two positive (tension) load magnitudes and two negative (compression) ones. For reference, the undeformed rod axis is also shown (zero loading). The numerical results of the present model exhibit good agreement with the analytical results, even for extremely large deformations. 6.3. Spatial deformation of a circular rod Fig. 12 presents the geometry of an undeformed cantilevered circular rod, lying in the plane X–Z. The rod is clamped at its root and a concentrated force F, normal to the plane of the rod (in the Y direction), is applied at the free tip. The rod radius is R and its opening angle is equal to 45°. It has a square cross-section with bending stiffness (E  I) about its principle axes. This problem was analyzed by Bathe and Bolourchi [15] using a finite element approach.

8 Integration 7

Transfer Matrix Linear Vy

6

P ⋅ R2 E⋅I

5

4

Vy ⋅ R2 E⋅I 3

6.2. Planar bending of a circular rod

2 As a first step of examining the capability of the new model to analyze large deformations of curved rods, the planar bending of a circular rod is studied as shown in Fig. 9. A circular rod is acted upon by a tensile force of magnitude 2F (when F is negative, compression is considered). Instead of solving the original problem, an equivalent problem of a rod that before deformation has the shape of a quarter circle, is considered (see Fig. 9). This quarter circle rod is clamped at the root and pulled by a force F at the tip. As the quarter ring deforms, the direction of the elastic axis at the free tip, where the force acts, is forced to coincide with the direction

Vz ⋅ R2 E⋅I

1

0

0

1

2

3

4

5

6

F ⋅ R2 7 EI

8

Fig. 16. Spatial bending of a circular rod – components of the cross-sectional resultant force at X0/X0-tip = 0.9, as function of the applied tip force.

480

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

0.12

0.6 Mx ⋅ R E⋅I

Integration Transfer Matrix

F ⋅ R2 =8 EI

6

0.08

4

0.06

F ⋅ R2 =8 EI

Integration Transfer Matrix

My ⋅R E0.1 ⋅I

6

0.4

0.3 4

0.04 2

0.2

0.02

0.1

2

0

-0.02

0 0

0.25

0.5

0.75

0

X0 1 X 0−tip

0.25

0.5

0.75

X0 1 X 0−tip

2 Mz ⋅ R E⋅I

F ⋅R =8 EI 2

1.6

Integration Transfer Matrix

1.4 1.2

6

1 4

0.8 0.6 2

0.4 0.2 0 0

0.25

0.5

0.75

X0 1 X 0−tip

Fig. 17. Spatial bending of a circular rod – distributions along the rod of the cross-sectional resultant moment components.

Fig. 13 presents the components of the tip displacement in the directions (X, Y, Z), as functions of the normalized load, (F  R2/E  I). The linear analysis predicts a displacement in the direction of the applied force, while the other two components are equal to zero. The results of the converged present nonlinear model practically coincide with the results of Ref. [15]. In Fig. 13 the convergence of the results with the number of nodes is also shown. The influence of the number of nodes increases with the load. Twenty-five nodes give results that deviate by less than 3% from the converged results. Fig. 14 presents the maximum local elastic rotation component along the entire rod. Figs. 13 and 14 indicate that if this rotation is smaller than 0.05, then the results are within 3% of the converged ones. Rod deflections are typically easier to obtain accurately than forces and moments. One of the advantages of the present formu-

lation is that resultant cross-sectional forces and moments are obtained directly as part of the solution procedure and are very accurate. Fig. 15 presents distributions along the rod of the three components of the cross-sectional resultant force, P, Vy, and Vz, at increasing applied tip force, F. For each applied force the components are given as functions of the non-dimensional coordinate of that cross-section at zero load, X0/X0-tip. The direct results of the calculations are compared with results that are based on integrating the externally applied forces from the free tip up to the specific cross-section, and then projecting the resultant in the directions of the local system of coordinates (x, y, and z) after deformation. In the present case the integration is trivial because it includes only the applied tip force. Excellent agreement is shown between both results (practically the curves coincide). In the linear case only the component Vy is different from zero. Thus P and Vz are

481

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

2 Integration Transfer Matrix Linear 1.6

Mz ⋅ R E⋅I

1.2

Mx ⋅ R E⋅I

0.8

0.4

0

0

2

4

F ⋅ R2 EI 6

8

Fig. 18. Spatial bending of a circular rod – components of the cross-sectional resultant moment at X0/X0-tip = 0, as functions of the applied tip force.

solely due to nonlinear effects. Moreover, the nonlinear effects influencing Vy are also evident (in the linear case Vy is constant along the rod and equal to the applied force). Fig. 16 presents the components of the cross-sectional resultant force at X0/X0-tip = 0.9, as functions of the applied force, F. Nonlinear effects become very significant as the applied force increases, yet the agreement between the results of the present model and direct integration remains excellent. Fig. 17 presents distributions along the rod of the components of the cross-sectional resultant moment, at different applied forces. The direct results of the calculations are compared again with results that are based on integration of the contributions of the external applied forces to the resultant cross-sectional moment along the rod. These calculations include the influence of the deflections along the rods on the moment arms. The cross-sectional resultant moments are projected onto the directions of the local coordinates after the deformation. The agreement between both results is again excellent. Nonlinear effects are evident and significant. My is solely due to nonlinear effects. Fig. 18 presents the components of the resultant cross-sectional moment at the root, as functions of the magnitude of the applied tip force. The nonlinear effects are significant. The agreement between the results of the present model and direct integration is excellent along the entire range. 7. Conclusions The derivation of a nonlinear transfer matrix model to analyze the large spatial deformations of curved and pretwisted rods has been presented. The rod is divided into small segments. A local body Cartesian system of coordinates is attached to each segment. Due to the large deformations, these local systems of coordinates

experience large translations and rotations. Yet, if the segments are chosen small enough, the deformations of the segment relative to its local body system of coordinates, are small. The nonlinearity appears in the rotation matrix between neighboring systems of coordinates, which is a function of the initial curvature and pretwist, as well as the elastic deformations. The rotations due to the elastic deformations are superimposed on the initial rotations due to the initial curvature and pretwist. The rotations between the segments are treated as finite ones. In spite of the introduction of nonlinear effects, the formulation is relatively simple, the solution procedure is efficient and the implementation of the model in a code is straightforward. Thus, the main advantages of the transfer matrix method exist also after the introduction of high nonlinearity. Because of the present formulation that separates between the segment’s deformations and finite rotations between the segments, the strains are not affected by the finite rigid body rotations or translations. Thus the objectivity criterion is met. The model can deal with any combination of boundary conditions at the edges of the rod. The present derivation has been confined to an isotropic material. Bernoulli–Euler assumption has been applied for bending, while Saint-Venant torsion is assumed. Yet, the extension of the model to other cases that may include: anisotropic materials, Timoshenko beam, thin cross-sections etc. is straightforward. Only the linear relations between the local deformations and the resultant forces and moments, have to be adapted. An increase in the number of the unknown variables can be handled quite easily, while using the present approach of nonlinear analysis. The local elastic rotations can be used as a measure of the accuracy of the results. Keeping the maximum local elastic rotation small ensures good convergence. Thus it is possible to control the accuracy by increasing the number of segments at regions where the local maximum elastic rotations exceed predetermined limits. The accuracy of the model has been validated by comparing its results with analytical, experimental, and finite-element results from the literature. The agreement between the results is very good. The cross-sectional resultant forces and moments are integral part of the present formulation and solution procedure. They are very accurate even in the case of large nonlinearities, representing another advantage of the present method. Because of the simplicity and efficiency of the new model, it can easily be added, as a sub-model, to comprehensive codes for simulating complicated systems. The present formulation can very easily cope with complex boundary conditions, loads and dynamic phenomena. Appendix A. The geometry of the rod axis before and after deformation The geometry of the axis of the rod is described relative to the (X, Y, Z) system of coordinates. A.1. Before deformation Before deformation, Ti is the rotation matrix between the systems (x(i1), y(i1), z(i1)) and (x(i), y(i), z(i)), see Eq. (1). If S(i) is the rotation matrix between the systems (X, Y, Z) and (x(i), y(i), z(i)), then:

SðiÞ ¼ TðiÞ  Sði1Þ :

ða:1Þ

The position vector of the ith nodal point is:

b þ dyðiÞ  Y b þ dzðiÞ  Z b rðiÞ ¼ dxðiÞ  X

ða:2Þ

b Y; b ZÞ b is a triad of unit vectors in the directions ðX; Y; ZÞ, respecð X; tively.It is clear that:

482

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

dðiÞ ¼ STði1Þ  lði1Þ þ dði1Þ ;

ða:3Þ

hxðiÞ ðxðiÞ Þ ¼ uðiÞ ðxðiÞ Þ;

where: T

According to Fig. b.1 the following geometric relations exist:





dðiÞ ¼ dxðiÞ ; dyðiÞ ; dzðiÞ ;

ða:4aÞ

T lðiÞ

ða:4bÞ

¼ ½LðiÞ ; 0; 0;

d(0) and T(0) represent rigid body translation and rotation, respectively, of the entire rod, and are assumed to be known. A.2. After deformation

sin hyðiÞ ðxðiÞ Þ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi



ffi;

0

ða:5Þ 0

0

Appendix B. The elastic rotation matrix The elastic rotation matrices are used to calculate the rotation matrix after the deformation, T0ðiÞ , as shown by Eq. (17). The elastic rotation of any point along the axis of segment (i), relative to the system (x(i), y(i), z(i)) after deformation, is described by three Euler angles. Fig. b.1 describes an element of the rod axis that undergoes these three rotations. The sequence of the rotations is as follows: A rotation hz(i)(x(i)) of the system of coordinates about the axis z(i), followed by a rotation hy(i)(x(i)) about the new direction of the axis y(i) (after the rotation hz(i)(x(i))), and finally a rotation hx(i)(x(i)) about the new direction of the axis x(i). The elastic rotation matrix in this case is:

ðb:3Þ

dv ðx Þ 2

ðiÞ ðiÞ 1þ dxðiÞ cos hyðiÞ ðxðiÞ Þ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi



ffi;

dv ðiÞ ðxðiÞ Þ 2 dxðiÞ

þ

dwðiÞ ðxðiÞ Þ 2 dxðiÞ

dv ðiÞ ðxðiÞ Þ dxðiÞ

ðb:4Þ

sin hzðiÞ ðxðiÞ Þ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

;

ðb:5Þ

1 cos hzðiÞ ðxðiÞ Þ ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

:

ðb:6Þ

1þ

0

The calculation of the components dxðiÞ ; dyðiÞ and dzðiÞ is carried out in the same manner as for the undeformed rod, where T0ðiÞ (see Eq. 0 (17)) replaces T(i). Again dð0Þ and T0ð0Þ represent rigid body translation and rotation of the deformed rod, and are assumed to be known input to the present problem. The matrices S0ðiÞ are calculated in the same manner as the matrices S(i), by using Eq. (a.1). Often these matrices are also needed for calculating the components of the loads in the directions (x(i), y(i), z(i)) after deformation.

dwðiÞ ðxðiÞ Þ 2 dxðiÞ

þ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1þ

b þd Y b þ d  Z: b RðiÞ ¼ dxðiÞ  X yðiÞ zðiÞ 0

dv ðiÞ ðxðiÞ Þ 2 dxðiÞ

1þ

After deformation the position vector of the ith nodal point is: 0

ðb:2Þ dwðiÞ ðxðiÞ Þ  dx ðiÞ

1þ

dv ðiÞ ðxðiÞ Þ 2 dxðiÞ

dv ðiÞ ðxðiÞ Þ 2 dxðiÞ

Appendix C. The case of complex boundary conditions The present model can cope directly with cases of statically determinate structures. The application of the model to cases of statically indeterminate structures is more complicated. In these cases less than six components of the force or moment reactions at the edges are known a priori. Cases of statically indeterminate rods belong to the more general case of complex boundary conditions. In the case of complex boundary conditions the reactions at the rod edges are functions of the displacements and rotations there. The root and tip of the curved rod are nodes i = 0 and (N + 1), respectively. The three components of the displacements at those nodes are: a(i), b(i), and c(i)(i = 0, N + 1). In most cases it is convenient to give these components relative to the inertial system of coordinates: X, Y, Z. The rotations at the root and tip are also described by three parameters: h(i), s(i), and t(i) (i = 0, N + 1). It is also convenient in many cases to describe the rotations relative to the system (X, Y, Z). h(i), s(i), and t(i) can represent Euler angles,

          2 3 cos hyðiÞ ðxðiÞ Þ  cos hzðiÞ ðxðiÞ Þ cos hyðiÞ ðxðiÞ Þ  sin hzðiÞ ðxðiÞ Þ  sin hyðiÞ ðxðiÞ Þ                        7 6 EðiÞ ðxðiÞ Þ ¼ 4 sin hxðiÞ ðxðiÞ Þ  sin hyðiÞ ðxðiÞ Þ  cos hzðiÞ ðxðiÞ Þ  cos hxðiÞ ðxðiÞ Þ  sin hzðiÞ ðxðiÞ Þ sin hxðiÞ ðxðiÞ Þ  sin hyðiÞ ðxðiÞ Þ  sin hzðiÞ ðxðiÞ Þ þ cos hxðiÞ ðxðiÞ Þ  cos hzðiÞ ðxðiÞ Þ sin hxðiÞ ðxðiÞ Þ  cos hyðiÞ ðxðiÞ Þ 5                         cos hxðiÞ ðxðiÞ Þ  sin hyðiÞ ðxðiÞ Þ  cos hzðiÞ ðxðiÞ Þ þ sin hxðiÞ ðxðiÞ Þ  sin hzðiÞ ðxðiÞ Þ cos hxðiÞ ðxðiÞ Þ  sin hyðiÞ ðxðiÞ Þ  sin hzðiÞ ðxðiÞ Þ  sin hxðiÞ ðxðiÞ Þ  cos hzðiÞ ðxðiÞ Þ cos hxðiÞ ðxðiÞ Þ  cos hyðiÞ ðxðiÞ Þ ðb:1Þ

e z (i )

dx(i ) ⋅ 1 +

dv(i )

e y(i) dx(i ) −

2

dx(i )

dx(i ) ⋅ 1 +

dv(i ) dx(i )

2

+

dw(i )

2

dx(i )

⋅ dx(i )

ϕ(i)

θ y (i )

θ z(i)

dw(i )

dv(i ) dx(i )

⋅ dx(i )

dx(i ) Fig. b.1. Rotation of an element dx(i) by three Euler angles.

e x (i )

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

Euler parameters or any other set of parameters that define a rotation. A vector of boundary displacement and rotation is defined as follows:

  kTðiÞ ¼ aðiÞ ; bðiÞ ; cðiÞ ; hðiÞ ; sðiÞ ; t ðiÞ ;

i ¼ 0; ðN þ 1Þ

ðc:1Þ

In the case of complex boundary conditions:

~ðiÞ ¼ KðiÞ  kðiÞ ; l

i ¼ 0; ðN þ 1Þ

ðc:2Þ

~ðiÞ is a vector of order N(i), where i = 0, (N + 1), and: l

NðiÞ  6 for i ¼ 0; ðN þ 1Þ:

ðc:3Þ

~ðiÞ includes all, or part, of the components of l(i). The vector l KðiÞ is a N(i)  6 ‘‘Stiffness” matrix. In general K(i) is a nonlinear ~ðiÞ . function of k(i) and l As an example, consider the case where there is a linear spring at the root (i = 0), which is defined by linear spring constants kx, ky, and kz in the directions X; Y, and Z, respectively. In addition there is also an angular spring at the root, having spring constants mx, my, and mz, about the axes X, Y, and Z, respectively. a(0), b(0), and c(0), define the root displacement components in the directions X, Y, and Z, respectively. In addition there are small root rotations h(0), s(0), and t(0), about the axes X, Y, and Z, respectively. In this case:

Nð0Þ ¼ 6

ðc:4Þ

and:

2  Kð0Þ ¼

Tð0Þ

3

kx

6 6  6 6 6 Tð0Þ 6 6 6 4

ky

7 7 7 7 7: 7 7 7 5

0 kz

mx 0

my

ðc:5Þ

mz An iterative solution procedure is applied in the case of complex ~ðiÞ is calboundary conditions: Initial value of kðiÞ is assumed. Then l culated according to Eq. (c.2). Afterwards, the iterative procedure that is described Section 5 of the text is followed. At the end of each iteration cycle kðiÞ is calculated according to Appendix A, and then the procedure is repeated until convergence. If the rod is fully clamped at the root, then:

kx ; ky ; kz ; mx ; my ; mz ! 1:

ðc:6Þ

In this case an iterative procedure is applied where the root reactions for the case of zero displacement and zero rotation at the root, are calculated. References [1] Frisch-Fay R. Flexible bars. London: Butterworth; 1962. 220p. [2] González C, LLorca J. Stiffness of a curved beam subjected to axial load and large displacements. Int J Solids Struct 2005;42(5–6):1537–45. [3] Pi YL, Bradford MA, Tin-Loi F. Nonlinear analysis and buckling of elastically supported circular shallow arches. Int J Solids Struct 2007;44(7–8): 2401–25. [4] Libai A. Intrinsic equations for the nonlinear dynamics of space beams. AIAA J 1996;34(8):1657–63. [5] Atluri SN, Iura M, Vasudevan S. A consistent theory of finite stretches and finite rotations in space-curved beams of arbitrary cross-section. Comput Mech 2001;27(4):271–81. [6] Pi YL, Bradford MA, Uy B. Nonlinear analysis of members curved in space with warping and Wagner effects. Int J Solids Struct 2005;42(11–12):3147–69. [7] Rubin MB. Numerical solution procedures for nonlinear elastic curved rods using the theory of a Cosserat point. Math Mech Solids 2005;10:89–126.

483

[8] Tufekci E, Dogruer OY. Exact solution of out-of-plane problems of an arch with varying curvature and cross section. J Eng Mech 2006;132(6):600–9. [9] Bauchau OA, Hong CH. Large displacement analysis of naturally curved and twisted composite beams. AIAA J 1987;25(1):1469–75. [10] Yu W, Hodges DH, Volovoi V, Cesnik CES. Timoshenko-like modeling of initially curved and twisted composite beams with oblique cross sections. In: AIAA Paper 2000-1405, 41st AIAA/ASME/ASCE/AHS/ASC Structures, structural dynamics, and materials conference and exhibit, Atlanta, Georgia, 3rd–6th April, 2000. [11] Pai PF, Lee SY. Three-dimensional post-buckling analysis of highly flexible curved beams. In: AIAA Paper 2003–1933, 44th AIAA/ASME./ASCE/AHS structures, structural dynamics, and materials conference, Norfolk, Virginia, 7th–10th April, 2003. [12] Moon-Young K, Kim II N, Kim S-B. Spatial stability of shear deformable curved beams with non-symmetric thin-walled sections I: Stability formulation and closed form solutions. Comput Struct 2005;83(31–32):2525–41. [13] Moon-Young K, Kim II N, Kim S-B. Spatial stability of shear deformable nonsymmetric thin-walled curved beams: a centroid-shear center formulation. J Eng Mech 2006;132(12):1313–25. [14] Rosen A, Rand O. Numerical model of the nonlinear behavior of curved rods. Comput Struct 1986;22(5):785–99. [15] Bathe KJ, Bolourchi S. Large displacement analysis of three-dimensional beam structures. Int J Numer Methods Eng 1979;14(7):961–86. [16] Surana SK, Sorem MR. Geometrically non-linear formulation for three dimensional curved beam elements with large rotations. Int J Numer Methods Eng 1989;28(1):43–73. [17] Dutta A, White DW. Large displacement formulation of a three-dimensional beam element with cross-sectional warping. Comput Struct 1992;45(1): 9–24. [18] Lo SH. Geometrically nonlinear formulation of 3D finite strain beam element with large rotations. Comput Struct 1992;44(1/2):147–57. [19] Bathe KJ. Finite element procedures. New Jersey: Prentice-Hall; 1996. [20] Ibrahimbegovic´ A. On finite element implementation of geometrically nonlinear Reissner’s beam theory: three-dimensional curved beam elements. Comput Methods Appl Mech Eng 1995;122(1):11–26. [21] Yoon YK, Kim Y. A one-dimensional theory of thin-walled curved rectangular box beams under torsion and out-of-plane bending. Int J Numer Methods Eng 2002;53(7):1675–93. [22] Zupan D, Saje M. Finite-element formulation of geometrically exact threedimensional beam theories based on interpolation of strain measures. Comput Methods Appl Mech Eng 2003;192(49–50):5209–48. [23] Kapania RK, Li J. A formulation and implementation of geometrically exact curved beam elements incorporating finite strains and finite rotations. Comput Mech 2003;30(5–6):444–59. [24] Mukhopadhyay M, Sheikh AH. Large amplitude vibration of horizontally curved beams: a finite element approach. J Sound Vibration 1995;180(2):239–51. [25] Crisfield MA. A consistent co-rotational formulation for non-linear, three dimensional beam-element. Comput Methods Appl Mech Eng 1990;81: 131–50. [26] Shabana AA. Flexible multibody dynamics: review of past and recent developments. Multibody Syst Dyn 1977;1(2):189–222. [27] Belytschko T, Hsieh BJ. Nonlinear transient finite element analysis with convected coordinates. Int J Numer Methods Eng 1973;7:255–71. [28] Cardona A, Géradin M. Time integration of the equations of motion in mechanism analysis. Comput Struct 1989;33(3):801–20. [29] Bauchau OA. Computational schemes for flexible, nonlinear multi-body systems. Multibody Syst Dyn 1998;2(2):169–225. [30] Géradin M, Cardona A. Flexible multibody system: a finite element approach. New York: John Wiley & Sons; 2001. [31] Crisfield MA, Jelenic´ G. Objectivity of Strain measures in the geometrically exact three-dimensional beam theory and its finite-element implementation. In: Proceedings of the Royal Society, London: Mathematical, Physical and Engineering Sciences, vol. 455(1983), 1999. p. 1125–47. [32] Jelenic´ G, Crisfield MA. Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for static and dynamics. Comput Methods Appl Mech Eng 1999;171:141–71. [33] Leckie F, Pestel E. Transfer-matrix fundamentals. Int J Mech Sci 1960;2:137–67. [34] Hayashi N. On the explicit representation of transfer-matrices of curved beams and their applications. Bull JSME 1975;18(122):797–806. [35] Irie T, Yamada G, Takahashi IS. The steady state in-plane response of a curved Timoshenko beam with internal damping. Ing Arch 1980;49(1):41–9. [36] Nitzsche F. A revised version of the transfer matrix method to analyze onedimensional structures. In: AIAA Paper 83-0825, 24th structures, structural dynamics and materials conference, Lake Tahoe, Nevada, 2nd–4th May, 1983. [37] Haktanir V, Kiral E. Statical analysis of elastically and continuously supported helicoidal structures by the transfer and stiffness matrix methods. Comput Struct 1993;49(4):663–77. [38] Kashimoto H, Shiraishi A, Nagaya K. Dynamic stress concentration in an arbitrary curved and inhomogeneous rod joined with infinite straight rods and excited by a twisting wave. J Sound Vibration 1994;178(3):395–409. [39] Yildirim V. Rotary inertia, axial and shear deformation effects on the in-plane natural frequencies of symmetric cross-ply laminated circular arches. J Sound Vibration 1999;224(4):575–89.

484

A. Rosen, O. Gur / Computers and Structures 87 (2009) 467–484

[40] Stephen NG, Ghosh S. Eigenanalysis and continuum modelling of a curved repetitive beam-like structure. Int J Mech Sci 2005;47(12):1854–73. [41] Rosen A. Structural and dynamic behavior of pretwisted rods and beams. Appl Mech Rev 1991;44(12):483–515. [42] Dowell EH, Traybar J. An experimental study of the nonlinear stiffness of a rotor blade undergoing flap, lag and twist deformation. AMS report no. 1194, Ames Research Center, January 1975.

[43] Timoshenko S, Goodier JN. Theory of elasticity. 2nd ed. Berlin: McGraw-Hill Inc.; 1951. [44] Bunce JW, Brown EH. Nonlinear bending of thin ideally elastic rods. Int J Mech Sci 1976;18(September–October):435–41.