Effect of member initial curvature on a flexible mechanism response

EFFECT OF MEMBER INITIAL CURVATURE ON

A FLEXIBLE MECHANISM RESPONSE

K. M. H  R. T. Y

Department of Mechanical Engineering, National Chiao Tung University, Hsinchu, Taiwan, Republic of China

(Received 16 February 1994, and in final form 15 February 1995)

A co-rotational finite element formulation of slender curved beam element is presented to investigate the effect of member initial curvature on the dynamic behaviour of planar flexible mechanisms. The Euler–Bernoulli hypothesis and the initial curvature are properly considered for the kinematics of curved beam. The nodal co-ordinates, incremental displacements and rotations, velocities, accelerations and the equations of motion of the system are defined in terms of a fixed global co-ordinate system, while the total strains in the beam element are measured in element co-ordinates which are constructed at the current configuration of the beam element. The element equations are constructed first in the element co-ordinate system and then transformed to the global co-ordinate system by using standard procedures. Both the deformation nodal forces and the inertia nodal forces of the beam element are systematically derived by consistent linearization of the non-linear beam theory in the element co-ordinates. An incremental-iterative method based on the Newmark direct integration method and the Newton–Raphson method is employed here for the solution of the non-linear dynamic equilibrium equations. Numerical examples are presented to demonstrate the effectiveness of the proposed element and to investigate the effect of the initial curvature on the dynamic response of the flexible mechanisms.

7 1996 Academic Press Limited

1. INTRODUCTION

The dynamic behaviour of flexible mechanisms has been the subject of considerable research [1–13]. Currently, the most popular approach for this analysis is to develop finite element models. Hsiao and Jang [13] proposed a co-rotational formulation for the dynamic analysis of planar flexible mechanisms, which was proven to be very effective by numerous numerical examples studied. However, only a straight beam element was considered.

In order to investigate the effects of member initial curvature on the dynamic behaviour of planar flexible mechanisms, a consistent co-rotational formulation of a slender curved beam with large rigid body motion and small elastic deformations is proposed in what follows here. It is assumed that the product of the initial curvature and the depth of the slender curved beam is much smaller than unity. As was done by Hsiao and Jang [10, 12], the nodal co-ordinates, incremental displacements and rotations, velocities, accelerations and the equations of motion of the system are defined in terms of a fixed global co-ordinate system, while the total strains in the beam element are measured in element co-ordinates which are constructed at the current configuration of the beam element. The element equations are constructed first in the element co-ordinate system and then transformed to the global co-ordinate system by using standard procedures. The element inertia nodal forces and deformation nodal forces are derived by using the d’Alembert principle and the virtual work principle. In order to capture all inertia effects and coupling between

177

of the motion is always small relative to the local element axes; thus, in conjunction with the co-rotational formulation, the higher order terms of the nodal parameters in the element deformation nodal forces and inertia nodal forces may be neglected by consistent linearization [13].

An incremental-iterative method based on the Newmark direct integration method and the Newton–Raphson method is employed here for the solution of the non-linear dynamic equilibrium equations. Numerical examples are presented to demonstrate the effectiveness of the proposed beam element and to investigate the effect of the initial curvature on the dynamic response of the flexible mechanisms.

2. NON-LINEAR FORMULATION

2.1.  

The following assumptions are made in the derivation of the non-linear behaviour: (1) the Euler–Bernoulli hypothesis is valid; (2) the unit extension of the centroid axis of the beam element is uniform; (3) the deflections of the beam element measured in the element co-ordinates are small; (4) the product of the initial curvature and the depth of the beam is much smaller than unity; (5) the strain of the beam element is small.

The third assumption can always be satisfied if the element size is properly chosen. Due to the assumption of small strain, the engineering strain and stress are used for the measure of the strain and stress. For convenience, the engineering strain is obtained from the corresponding Green strain in this study.

2.2. - 

In order to describe the system, following Hsiao and Jang [10, 12], we define two sets of co-ordinate systems (see Figure 1).

(1) The first is a fixed global set of co-ordinates, X1, X2; the nodal co-ordinates,

incremental displacements and rotations, velocities, accelerations and the equations of motion of the system are defined in this co-ordinate system.

(2) The second set is the element co-ordinates, x1, x2; a set of element co-ordinates

associated with each element, which is constructed at the current configuration of the beam element. The origin of this co-ordinate system is located at node 1, and the x1axis is chosen

to pass through the two end nodes of the element. Note that this co-ordinate system is a local co-ordinate system not a moving co-ordinate system. However, this co-ordinate system is updated at each iteration. The element equations are constructed first in the element co-ordinate system and then transformed to the global co-ordinate system for assemblage using standard procedure.

Figure 1. Co-ordinate systems.

2.3.    

The geometry of the beam element is described in the current element co-ordinate system. In this study, the symbol{} denotes a column matrix, and the symbol ( ) denotes the variable ( ) in the undeformed state. Let P (see Figure 2) be an arbitrary point in the beam element, and Q be the point corresponding to P on the centroid axis. The position vector of point P in the undeformed and deformed configurations can be expressed as

r¯ =

6

y¯p

7

=

6

v¯(s¯)+y cosu

7

Figure 3. Nodal displacements u1and u2. and r=

6

7

6

7

where xc(s) and v(s) are the x1 and x2co-ordinates of point Q, respectively, s is the arc

length of the deformed centroid axis measured from node 1 to point Q, y is the distance between points P and Q, andu is the angle measured from the x1-axis to the tangent of

the centroid axis. The relationship among xc(s) and v(s) and s may be given as

xc(s)=u1+

S 2

g

j

−1

cosu dj, cosu=(1−v'2₎1/2_{, v'=}dv(s)

ds =sinu,

j=1+2s_S, (3–6)

where u1 is the displacement of node 1 in the x1-direction, S is the current arc length of

the centroid axis of the beam element. Note that due to the definition of the element co-ordinate system, the value of u1 is equal to zero. However, the variation and time

derivatives of u1 are not zero. Making use of equation (3), one obtains

S=2l

>g

1

−1

cosu dj=2l/b and l=xc(S)−xc(0)=L−u1+u2, (7, 8)

in which l is the current chord length of the beam axis, L is the chord length of the undeformed beam axis, and u2is the displacement of node 2 in the x1-direction as shown

in Figure 3.

If s¯ and y in equation (1) are regarded as the Lagrange co-ordinates, the Green strains eij (i=1, 2; j=1, 2) are given by [14]

eij=12(G T iGj−gTigj), (9) where G1=1r_1s¯=1r_1s1s_1s¯=(1+e0)(1−ky)

6

7

6

7

6

7

6

7

in which e0is the unit extension of the centroid axis, c=1/cosu, v0=d2v/ds2, andk is an

exact expression for the physical curvature of the deformed beam centroid axis. An equivalent expression fork has been given by Hodges [15]. Making use of the assumption of uniform unit extension, one may rewrite the unit extension e0 in equation (11) as

e0=(S/S)−1. (13)

Due to the use of the Euler–Bernoulli hypothesis, as expected,e11is the only non-zero

component of eij, and is given by

e11=12[(1+e0)2(1−ky)2−(1−k¯y)2]. (14)

The engineering strain corresponding to e11 is given by [14]

e=

0

1

Note thate in equation (15) is an exact expression for the engineering strain for the curved Euler beam. When k¯yW1, the engineering strain in equation (15) may be approximated by e1e0−(k−k¯)y.

Here, the lateral deflection of the centroid axis, v(s) is assumed to be represented by cubic Hermitian polynomials in s and is given by

v(s)={N1, N2, N3, N4}T{v1, v'1, v2, v'2}=NTbub, (16)

where vj and v'j ( j=1, 2) are the nodal values of v and v' at nodes, j, respectively. Note

that, due to the definition of the element co-ordinates, the value of vjis zero. However,

the variations and time derivatives of vjare not zero. Ni(i=1–4) are shape functions and

are given by

N1=14(1−j)2(2+j), N2=(S/8)(1−j2)(1−j), N3=14(1+j)2(2−j),

N4=(S/8)(−1+j2)(1+j), (17)

where S is the current arc length of the centroid axis, and j is the non-dimensional co-ordinate defined in equation (6).

The initial lateral deflection of the centroid axis, v¯(s¯) is required for the calculation of the initial curvaturek¯ of the beam axis. For consistency, v¯(s¯), should have the same form as v(s). Here v¯(s¯) is obtained by adding an overbar for each variable in equations (6), (16) and (17).

Note that because the value of vj is zero, from equations (16) and (17), v' in equation

(5) can be determined by v'j and j only. Thus, the value of S in equation (7) can be

calculated by using equations (7), (8), (16) and (17), and then xc(s) in equation (3) can

be determined by using v' and S.

2.4.     

The global nodal parameters for the system of equations associated with the individual elements are chosen to be the incremental translationsDUjandDVj( j=1, 2) in the X1and

X2directions, respectively, and the incremental counterclockwise rotationsDuj( j=1, 2) at

nodes j (see Figures 1 and 4). The nodal forces corresponding to the nodal parametersDUj,

DVjandDuj( j=1, 2) at nodes j are forces Fijin Xi(i=1, 2) directions and counterclockwise

moments mj ( j=1, 2) at nodes j (see Figure 4). The element employed here has three

degrees of freedom per node defined in the current element co-ordinate system (see Figures

Figure 4. Nodal parameters and nodal forces.

in the x1and x2directions, respectively, and the counterclockwise rotationsuj( j=1, 2) at

nodes j. Note that the initial rotations are included in the measurement of uj. The nodal

forces corresponding to the nodal parameters ujand vjanduj( j=1, 2) at nodes j are forces

fij in xi (i=1, 2) directions and counterclockwise moments mj ( j=1, 2) at nodes j (see

Figure 4).

Due to the definition of the element co-ordinate system, the non-zero element nodal parameters are u2 anduj ( j=1, 2) only. LetuIj, X1jI and XI2j ( j=1, 2), and xiI (i=1, 2) be

the total nodal rotations, global nodal co-ordinates in the X1and X2axes at nodes j, and

the axes of the element co-ordinate system, respectively, for a single element at the equilibrium configuration of the Ith increment. The procedure used to determine the current u2 and uj ( j=1, 2) for individual elements corresponding to the incremental

translations DUjand DVj ( j=1, 2), and the incremental nodal rotationsDuj ( j=1, 2) at

nodes j (see Figure 1) is the same as that given in reference [12]. However, it is repeated here for completeness.

X1j and X2j ( j=1, 2), the current global nodal co-ordinates of the element at nodes j,

are obtained by addingDUjandDVj to X1jI and X2jI, respectively, and then the xi(i=1, 2)

axes, the axes of the current element co-ordinate system, can be constructed. The current chord length of the element l can be calculated by l=[(X12−X11)2+(X22−X21)2]1/2, and then

the value of u2can be determined by u2=l−L, where L is the initial chord length of the

element. Leta (Figure 1) denote the angle of rigid body rotation measured from the xI

1-axis

to the current x1-axis; then the current value of uj ( j=1, 2) is given by uj=uIj+Duj−a.

2.5.     

The element internal nodal forces are obtained from the d’Alembert principle and the virtual work principle in the current element co-ordinate system. The virtual work principle requires that duT afa+duuTbfb=

g

in which fd

j and f i

j ( j=a, b) are the deformation nodal force vector and the inertia nodal

force vector, respectively.de is the variation of e given in equation (15). s=Ee is the normal stress, where E is the Young’s modulus. r is the density, dr is the variation of r given in equation (2) with respect to the nodal parameters, and r¨=d2_r/dt2_{. In this paper, the symbol}

(  ) denotes differentiation with respect to time t. V is the volume of the undeformed beam.

For a curved beam the differential volume dV may be expressed as

dV=(1−k¯y) dA ds¯=(1−k¯y) dA ds/(1+e0), where dA is the differential cross section area.

Note that the element co-ordinate system is just a local co-ordinate system, which is updated at each iteration, not a moving co-ordinate system. Thus, r¨ is the absolute acceleration and fVrdrTr¨ dV comprises all the virtual inertia forces.

The exact expressions for faand fbmay be obtained by substituting the exact expressions

forde, e, dr and r¨ into equation (18). However, if the element size is properly chosen, the values of the nodal parameters (displacements and rotations) of the element defined in the current element co-ordinate system, which are the total deformational displacements and rotations, may always be much smaller than unity. Thus only the first order terms of the nodal parameters are retained in fd

aand fdb. However, in order to include the effect of axial

force on the lateral forces, a second order term of the nodal parameters is retained in fd b.

Because the values of the nodal parameters of the element may always be much smaller than unity, it is reasonable to assume that the coupling between the nodal parameters and their time derivatives are negligible. Thus only zeroth order terms of nodal parameters are retained in fi

a and fib. The approximations, 1−k¯y11, 1+e011, v'1u, v˙'1u, v¨'1u and

cosu11 are used in the derivation of fa and fb. In order to avoid improper omission in

the derivation of faand fb, these approximations are applied to the exact expressions for

de, e, dr and r¨. Note that because the shape functions of v' and v0 are functions of S, which is given in equation (7), the variations and differentiation of the shape functions are considered here.

From equations (7), (8), (12), (13) and (15), de may be expressed as

de=[1/(1−k¯y)][(1−ky)de0−(1+e0)ydk], (21) where dk=cdv0+c3_v'v0dv', de 0=dS/S=(2dl/bS)−(2ldb/b2S), (22, 23) dl={−1, 1}T_{du 1,du2}=GTadua, db=−

g

dv'=(dN'T b)ub+N'bTdub=−[de0/(1−e0)]NTcub+N'bTdub, (26) and dv0=(dN0T b )ub+N0bTdub=−[de0/(1+e0)]NTdub+N0bTdub, (27) where dub={dv1, cosu1du1,dv2, cosu2du2}, Nc={N'1, 0, N'3, 0}, Nd={2N01, N20, 2N03, N40}. (28–30)

Note that due to the definition of the element co-ordinate system, the values of v1 and v2

are zero. Thus NT

g

From equations (2)–(8), (16) and (17), dr can be expressed as

dr=

8

g

g

9

b)ub+NTbdub=[de0/(1+e0)]NTeub+NTbdub and Ne={0, N2, 0, N4}.

(35) Substituting equations (23)–(26) and (34) into equation (33), using the approximations 1+e011, v'1u, and cos u11, and retaining only zeroth order terms of nodal parameters,

one obtains dr=

6

7

6

7

From equations (2)–(8), (16) and (17), the exact expression for r¨ may be obtained. Using the approximations 1+e011, v'1u, v˙1u, v¨'1u and cos u11, and retaining only zeroth

order terms of nodal parameters in the exact expression of r¨, one obtains

r¨ =

g

G

G

F

f

g

g

h

G

G

J

j

u˙a={u˙1, u˙2}, u¨a={u¨1, u¨2}, and u˙b={v˙1,u1, v˙2,u2),

u¨b={v¨1,u1, v¨2,u2). (39, 40)

Note that u˙jand u¨j( j=a, b) are the absolute velocity and acceleration vectors of an element

referred to the element co-ordinate system which are obtained from the transformation of the corresponding global velocity and acceleration vectors extracted from the equations of motion of the system by using standard procedures [9].

Substituting equations (15), (31), (36) and (38) into equation (18), using the approximations 1−k¯y11, 1+e011, v'1u, and cos u11, dropping higher order terms

of nodal parameters, and equating the terms in both sides of equation (18) corresponding to virtual displacement vectors dua and duub, respectively, one obtains

fd a=AEe0Ga, fdb=kgub+kb(ub−u¯b), (41) kb= EIS 2

g

g

g

$

g

g

%

g

g

g

g

g

where A is the cross-section area, I=fAy2dA. e0 in equations (41) and (42) is given in

equation (13). Kband kgin equation (42) are bending and geometric stiffness matrices of

conventional beam element [12]. main equations (43) is the consistent mass matrix of bar

element [12]. mtand mrin equations (44) are the consistent mass matrices of conventional

beam element for lateral translation and rotation, respectively [12]. The underlined terms in equations (43) and (44) are called velocity coupling terms in this study.

Note that even equations (42)–(44) are functions of S, because the approximation 1+e011 is used in the derivation of equations (42)–(44); S is approximated by S in

equations (42)–(44). 2.6.  

The element stiffness matrices and mass matrices may be obtained by differentiating the element nodal force vectors in equations (41)–(44) with respect to nodal parameters, and time derivatives of nodal parameters. However, element matrices are used only to obtain predictors and correctors for incremental solutions of the non-linear equations in this study. Approximate element matrices can meet these requirements. Thus, the conventional stiffness matrix of a bar element

ka= AE S

$

%

and conventional element matrices given in equations (41)–(44) are adopted here for simplicity.

2.7.   

The non-linear equations of motion may be expressed as [13]

ew=>RTDU>/N>RTiDUi>Eewtol, (47–49)

where the subscript i of a quantity is used to denote that the quantity corresponds to the ith iteration; N is number of the equations of the system;DU is the displacement correction; eftol, edtoland ewtolare prescribed values of error tolerance. The error tolerance is set to 10−10,

and i is chosen to be 1 in this paper.

3. NUMERICAL ALGORITHM

An incremental iterative method based on the Newmark direct integration method [12, 16] and the Newton–Raphson method is employed. For clearness, the numerical procedure used here is given as follows.

Assume that the dynamic equilibrium configuration at time tnis known. Let Qnand Qn

denote the velocity and acceleration vectors of the discretized system at time tn,

respectively. DQ, Qn+1and Qn+1, the incremental displacement, velocity and acceleration

vectors at time tn+1=tn+Dt, respectively, may be obtained by the following

incremental-iterative procedure.

The initial incremental displacement vector DQ for the next time step may be chosen to be DQ=0. Then, from the Newmark method, Qn+1 and Qn+1 can be expressed as

Q n+1=_aDtDQ2−

Q n

aDt−

0

1

The incremental translations DUj and DVj, and the incremental nodal rotations Duj

( j=1, 2) of individual elements at element nodes j can be extracted from DQ. By using the method described in section 2.4, the current configuration and deformation of individual elements can be calculated. Then, by using equations (41) and (42), the element deformation nodal force vectors can be calculated in the current element co-ordinate system.

One can extract the global velocity and acceleration vectors of individual elements at element nodes j ( j=1, 2) from Qn+1 and Qn+1. Then, the absolute velocity and

acceleration vectors (equations (39) and (40)), of individual elements, referred to the element co-ordinate system, can be obtained from the transformation of the corresponding global velocity and acceleration vectors to the current element co-ordinate system by using standard procedures. Then, by using equations (43) and (44), the element inertia nodal force vectors can be calculated in the current element co-ordinate system.

Figure 5. Slider–crank mechanism.

From the assemblage of the element nodal force vectors and equation (46), the unbalanced force vector R can be calculated. If the convergence criteria in equations (47)–(49) are not satisfied, a displacement correctionDU is added to the previous DQ to obtain a new incremental displacement for the next iteration. The displacement correction DU may be determined by using the Newton–Raphson method as

DU=−K−1_{R (52)}

where K is the co-called effective stiffness matrix and may be expressed by

K =(1/aDt2_{)M+K, (53)}

where a is the parameter of the Newmark method. M and K are the mass matrix and the stiffness matrix of the system of equations, which are assembled from element matrices. This procedure is repeated until the convergence criteria are satisfied.

Figure 6. Midpoint deflection of coupler for slider-crank mechanism forv=124·8 rad/s, h/l=0 and ms=0.

sidered for the mass of the sliding block ms: (a) ms=0, (b) ms=0.03781 kg

(2.159×10−4_{lb s}2_{/in), (c) m}

s=0.07562 kg (4.318×10−4lb s2/in). The initial centroid axis of

the coupler is represented by a one-half sine function written as v¯=h sin (px/l), where

Figure 7. Midpoint deflection of coupler for slider–crank mechanism for v=124·8 rad/s. (a) ms=0; (b)

h is the amplitude of the sine function and l is the initial chord length of the couple. The values of h/l are chosen to be 0, 0·01 and 0·02 for this example.

The initial crank angle is set to zero and the initial elastic deformations are assumed to be zero. Three different crank speed,v=124.8, 250, 375 rad/s are considered here. The initial velocity and acceleration of the mechanism are calculated by using the kinematics of rigid mechanisms. In the numerical study, the rigid crank OA is simulated by a flexible member with the same cross-section and density as the member AB of the Young’s modulus 1010_{times larger than that of member AB. The member OA is discretized by one}

element for all cases and AB is discretized by using different numbers of element for different cases. LetDt denote the time step size and N denote the number of elements used for member AB. Unless stated otherwise, Dt is chosen to be 10−4_{, 5×10}−5 _and

Figure 8. Midpoint deflection of coupler for slider–crack mechanism for v=250 rad/s. (a) ms=0; (b)

curvature on the midpoint deflection of the coupler is not negligible for all the cases

Figure 9. Midpoint deflection of coupler for slider–crank mechanism for v=375 rad/s. (a) ms=0; (b)

Figure 10. Four-bar mechanism.

studied. The midpoint deflection of the coupler increases with increasing h/l. The discrepancies among the solutions for different h/l increase with increasing msand/orv.

4.2. - 

A four-bar mechanism shown in Figure 10 is considered. The crank OA is considered to be rigid, and the initially curved coupler AB and rocker BC are considered to be flexible. ABand BC are made of the same material and have the same cross section. The geometry, inertia properties and the material properties of the four-bar mechanism are as follows: length of OA=0·3048 m (12 in), chord length of AB=0· 9144 m (36 in), length of

Figure 11. Midpoint deflection of coupler for four-bar mechanism forv=10p rad/s. – – –, Dt=10−3_{s, N=1;} — - —,Dt=10−3_{s, N=2;}----, Dt=5×10−4_{s, N=4.}

discretized by one element for all cases, and AB and BC are discretized by the same number of elements. Let N denote the number of elements used for member AB (BC). Unless it is stated otherwise, Dt is chosen to be 5×10−4_{, 2·5×10}−4_{, 2·5×10}−4 _{and 10}−4_{s, and N is}

chosen to be 4, 4, 4 and 12 for crank speedsv=10p, 20p, 30p and 60p rad/s, respectively. The transverse deflection of the midpoint of the coupler vd, being measured relative to

its initial deflection, was computed for different h/l and v. For the case h/l=0 and v=10p rad/s, different Dt and N were used for convergence and accuracy studies. Figure 11 shows the time histories of vd/l obtained when using Dt=10−3s, N=1,

Figure 12. Midpoint deflection of coupler for four-bar mechanism. (a) 10p rad/s; (b) 20p rad/s; (c) 30p rad/s; (d) 60p rad/s. – – –, h/l=0; —-—, h/l=0·01; ----, h/l=0·2.

Dt=10−3_{s, N=2, and} Dt=5×10−4_{s, N=4. The solutions corresponding to}

Dt=5×10−4_{s, N=4 are nearly identical to those obtained when using} Dt=2×10−4_s,

N=4 (not shown in Figure 11). Thus, the solutions corresponding toDt=5×10−4_{s, N=4}

may be taken to be the converged solutions. Note that the solutions corresponding to Dt=10−3_{, N=1 are nearly identical to those reported in reference [9] (not shown in}

Figure 11), which were obtained by using the same Dt and N. As can be seen from Figure 11, the solutions obtained by using Dt=10−3_{, N=1 are not accurate enough.}

Figure 12 shows the present results for different h/l and v. It can be seen that the effect of the initial curvature on the midpoint deflection of the coupler is not negligible for all the cases studied. The midpoint deflection of the coupler increases with increasingh/l. The discrepancies among the solutions for different h/l increase with increasing v.

5. CONCLUSIONS

A co-rotational finite element formulation of a slender curved beam element to investigate the effect of member initial curvature on the dynamic behaviour of planar flexible mechanisms has been presented. The Euler–Bernoulli hypothesis and the initial curvature are properly considered for the kinematics of curved beams. Both the inertia nodal forces and deformation nodal forces are systematically derived by consistent linearization of the fully geometrically non-linear beam theory by using the d’Alembert principle and the virtual work principle. In conjunction with the co-rotational formulation, the higher order terms of the nodal parameters in the element nodal forces are consistently neglected. From the numerical examples studied, it is found that the effects of the initial curvature on the dynamic response of flexible mechanisms is not negligible for all the cases studied. The midpoint deflection of the coupler increases with increasing h/l. The discrepancies among the solutions for different h/l increase with increasing v.

It is believed that the consistent co-rotational formulation for a curved beam element presented here may represent a valuable engineering tool for the dynamic analysis of planar flexible mechanisms with initially curved members.

ACKNOWLEDGMENTS

The research was sponsored by the National Science Council, Republic of China, under the contract NSC79-0401-E009-15.

REFERENCES

1. S. C. C and K. C. P 1975 ASME Journal of Engineering for Industry 97, 542–550. Dynamic response of a high-speed slider–crank mechanism with an elastic connecting rod.

2. J. O. S and E. J. H 1980 Computer Methods in Applied Mechanics and Engineering 24, 359–381. Dynamic analysis of planar flexible mechanisms.

3. A. K. N and N. F. K 1980 Computer Methods in Applied Mechanics and Engineering 23, 225–251. Nonlinear dynamic analysis of curved beams.

4. M. B and A. M 1982 ASME Journal of Mechanical Design 104, 159–167. Member initial curvature effects on the elastic slider–crank mechanism response.

5. J. C. S and V-Q 1987 Journal of Sound and Vibration 119, 487–508. The role of non-linear theories in transient dynamic analysis of flexible structures.

6. T. R. K, R. R. R and A. K. B 1987 Journal of Guidance, Control and Dynamics 10, 139–151. Dynamics of a cantilever beam attached to a moving base.

7. D. F. B and A. A. S 1989 Mechanism and Machine Theory 24, 411–429. Effect of beam initial curvature on the dynamics of deformable multibody systems.

1–14. Dynamic analysis of planar flexible mechanisms by co-rotational formulation.

13. K. M. H, R. T. Y and A. C. L 1994 International Journal for Numerical Methods in Engineering 37,75–89. A consistent finite element formulation for nonlinear dynamic analysis of planar beam.

14. T. J. C 1988 Continuum Mechanics. Englewood Cliffs, NJ: Prentice Hall.

15. D. H. H 1988 AIAA Journal 22, 1825–1827. Proper definition of curvature in nonlinear beam kinematics.

16. K. J. B 1992 Finite Element Procedures in Engineering Analysis. Englewood Cliffs, NJ: Prentice Hall.