# Finite element analysis of the geometric stiffening effect. Part 1_ a correction in the floating frame of reference formulation英文

Finite element analysis of the geometric stiffening effect. Part 1: a correction in the fl oating frame of reference formulation D Garc? ′a-Vallejo1, H Sugiyama2, and A A Shabana2? 1Department of Mechanical and Materials Engineering, University of Seville, Seville, Spain 2Department of Mechanical and Industrial Engineering, University of Illinois at Chicago, Chicago, Illinois, USA The manuscript was received on 13 May 2004 and was accepted after revision for publication on 10 November 2004. DOI: 10.1243/146441905X10041 Abstract:The fact that incorrect unstable solutions are obtained for linearly elastic models motivates the analytical study presented in this paper. The increase in the number of fi nite elements only leads to an increase in the critical speed. Crucial in the analysis presented in this paper is the fact that the mass matrix and the form of the elastic forces obtained using the absolute nodal coordinate formulation remain the same under orthogonal coordinate trans- formation. The absolute nodal coordinate formulation, in contrast to conventional fi nite element formulations, does account for the effect of the coupling between bending and exten- sion. Based on the analytical results obtained using the absolute nodal coordinate formulation, a new correction is proposed for the fi nite element fl oating frame of reference formulation in order to introduce coupling between the axial and bending displacements. In this two-part paper, two- and three-dimensional fi nite element models are used to study the problem of rotating beams. The models are developed using the absolute nodal coordinate formulation that allows for accurate representation of the axial strain, thereby avoiding the ill-conditioning problem that arises when classical displacement-based fi nite element formulations are used. In the fi rst part of the paper, the case of linear elasticity is considered and assumptions used in the fi nite element fl oating frame of reference formulation are investigated. In the second part of the paper, non-linear elasticity is considered. A rotating helicopter blade is simulated, and the complexity of the motion suggests the inclusion of rotary inertia, shear deformation, and non-linear elastic forces in order to obtain an accurate solution that does not suffer from the instability problem regardless of the number of fi nite elements used. Keywords: geometric stiffening, absolute nodal coordinate formulation, fl oating frame of reference formulation, rotating beams, stability 1INTRODUCTION Since the late 1950s, special attention has been given to the problem of rotating beams. One of the early investigations of rotating beams can be attributed to Schilhans [1], who developed the partial differen- tial equations for fl exural vibration of a rotating beam under steady state conditions. For more than 25 years, most of the research was concerned with planar one-dimensional beams undergoing steady state motion. In these models neither the Coriolis effect nor the coupling between extensional and fl ex- ural vibrations was considered. Several later investi- gations presented more detailed analysis that made provisions for including the centrifugal stiffening and Coriolis effects. Some of these investigations presented a more complete theory for dealing with the motion of a rotating beam attached to a base and subjected to a prescribed motion [2]. In some of these investigations, an elongation variable is introduced such that the bending displacement of a ?Corresponding author: Department of Mechanical Engineering (M/C 251), University of Illinois at Chicago, 2031 Engineering Research Facility, 842 W. Taylor street, Chicago, Illinois 60607, USA. 187 K01204# #IMechE 2005 Proc. IMechE Vol. 219 Part K: J. Multi-body Dynamics at PENNSYLVANIA STATE UNIV on September 18, 2016pik.sagepub.comDownloaded from point on a beam gives rise to an axial displacement, thereby producing the stiffening force necessary to eliminate the source of instability. It was also demonstrated that classical modal expansion formu- lations that fail to account for the stiffening effect lead to incorrect solution. The study of the dynamics of rotating beams has been a subject of interest because of its importance in many engineering applications such as helicopter blades. The theory of helicopters early documented the need for considering the coupling between fl ap or lead-lag motions with extensional motion of blades [3]. Recently, Ruzicka and Hodges [4] used a mixed fi nite element method in which axial forces are introduced as an independent variable in the dynamic equations of blades, and this method has been shown to account for the foreshortening effect. Simo and Vu-Quoc [5] presented a study in which they showed that the use of the linear elasticity theory leads to a softening term that makes the stiff- ness matrix singular at a certain angular velocity, regardless of the stiffness of the material. Using a power series expansion in terms of a small para- meter, they concluded that, at least, quadratic terms of the strain–displacement relationship must be considered in order to avoid the singularity. Boutaghou and Erdman [6, 7] presented a review and comparison between the existing non-linear rod theories and concluded that Simo and Vu-Quoc’s approach [5] gives accurate results. Substructuring techniques were used to model the coupling between the axial and bending dis- placements [8]. Within each substructure, normal vibration and correction modes were used with the assumptions of the linear elasticity theory. No quan- titative analysis, however, was given to determine the relationship between the size of the substructures and the critical speed at which instability may occur. Geometric stiffening effects have also been taken into account in multi-body system formu- lations by considering the reference confi guration of the fl exible body as a prestressed state, and hence including centrifugal effects not in the kinetic but in the potential energy [9, 10]. Other authors discussed the importance of the non-linear strain– displacement relationship when the fl oating frame of reference formulation is used [11–14]. In the fi eld of robotics, the interest in real-time control has motivated the study of simple alternatives in order to eliminate the instability [15]. Some of the approaches, however, failed to account for the geo- metric stiffening effect as a result of premature line- arizations of the dynamic equations [2, 15]. The absolute nodal coordinate formulation is a fi nite element procedure for the study of fl exible bodies that experience large rotations and defor- mations [16]. In this formulation, the elastic forces canbeformulatedusingdifferentapproaches including the linear theory of elasticity and non- linear continuum mechanics. When the fi nite ele- ment forces are formulated using the linear theory of elasticity, the resulting equations are conceptually equivalent to the equations obtained using substruc- turingtechniquesthatemploylinearelasticity theory. Because the elastic forces are defi ned using a local element frame, the use of the assumptions of linear theory can be justifi ed. In a recent investi- gation [17] it was shown that, for rotating beams, there always exists a critical value of the angular velocity at which the fi nite element solution based on the linear theory of elasticity leads to incorrect unstable solution regardless of the number of fi nite elements used in the model. In contrast, the use of anon-linearcontinuummechanicsapproach always leads to a stable solution regardless of the value of the angular velocity [17]. It is one of the goals of the present paper to explain the reason for the existence of a critical angular velocity if the linear theory of elasticity is used within each element and explain why increasing the number of fi nite elements amounts to a shift in the critical speed without eliminating it. Crucial in the analysis pre- sented in this paper is the fact that the mass matrix and the expression of the elastic forces remain the same under an orthogonal coordinate transform- ation. The use of this important property is discussed in this paper. This paper is organized as follows. In the following section, the absolute nodal coordinate formulation is briefl y reviewed. Section 3 presents an analysis of a rotating beam modelled using one element based on linear elasticity. Sections 4 and 5 discuss the results presented in reference [17], which revealed that an increase in the number of fi nite elements does not always lead to the correct stable solution of rotating beams. Sections 6 and 7 present a com- parison between the absolute nodal coordinate formulation and the fi nite element fl oating frame of reference formulation. The results of this compari- son suggest a method to correct for the element orientation in order to obtain a stable solution using the fl oating frame of reference formulation. A summary and the conclusions drawn from this study are presented in section 8. 2FINITE ELEMENT FORMULATION The formulation used in this work differs from classi- cal fi nite element formulations in the sense that the nodal coordinates contain all the information about deformation and rigid body motion [16]. Owing to this representation of motion it is possible to formu- late isoparametric beam and plate/shell elements 188D Garc? ′a-Vallejo, H Sugiyama, and A A Shabana Proc. IMechE Vol. 219 Part K: J. Multi-body DynamicsK01204# #IMechE 2005 at PENNSYLVANIA STATE UNIV on September 18, 2016pik.sagepub.comDownloaded from [16, 18–22]. The nodal coordinate vector consists of Cartesian coordinates of the position vector of the node in a global frame and partial derivatives of the position vector with respect to the element or body spatial parameters. Thus, the global position of an arbitrary point on element j of fl exible body i is interpolated as rij? Sijeij(1) where Sijis the element shape function, and eijis the vector of element nodal coordinates. Using this kin- ematic description, the same interpolation is used for the velocity of the arbitrary point, and can be written as follows ˙ r ij ? Sij ˙ e ij (2) Thevelocityvectorcanbeusedtoformulatethekinetic energy of the elements. The system kinetic energy can be obtained as follows T ? 1 2 X nb i?1 X nie j?1 e Vij rij ˙ r ijT ˙ r ij dV ? 1 2 ˙ e TM˙ e (3) whererijand Vijare the material density and the volume of element ij respectively, nbis the number of fl exible bodies, nieis the number of elements of body i, M is the constant mass matrix of the system, and e is the vector of system nodal coordinates. The displacement fi eld used in the absolute nodal coordinateformulationcanaccountforslope discontinuity without the need for using non-linear constraints to represent the connectivity between the fi nite elements. To this end, a different body parametrization is used [23]. The elastic forces can be formulated using a linear approach [16] by intro- ducing a local frame that follows the element rigid body motion. The formulation of the elastic forces for two-dimensional Euler–Bernoulli beam elements using this linear approach is described in Appendix 2. Non-linear formulations of elastic forces in the absolute nodal coordinate formulation have been also presented in the literature [19–22, 24]. Usingthekineticandstrainenergies,the equations of motion can be obtained using the Lagrangian approach and can be written as follows M ¨ e ? Qe? Fe(e)(4) where Qeis the vector of generalized external forces, and Fe(e) is the elastic force vector, which is non- linear even in the case of the linear theory of elasticity [19–22, 24]. 3ONE-ELEMENT MODEL BASED ON LINEAR ELASTICITY In order to explain the results obtained by Berzeri and Shabana [17], the problem of a rotating beam is now analysed. In this section, two-dimensional Euler–Bernoulli beam elements formulated using linear elasticity theory are used. Firstly, the equations of motion of a beam element are expressed in a rotat- ing frame (Fig. 1). As shown in Appendix 2, the expression of the stiffness matrix of the element is the same in both the global and the rotating frames. Nevertheless, as shown in what follows, the equationsintherotatingframeincludetwo additional terms: one is expressed in terms of a skew-symmetric matrix and represents the Coriolis effect, and the other is expressed in terms of a sym- metric matrix and represents the centrifugal inertia forces. Since the analysis presented in this section is focused on one element, superscripts i and j referring to the body and the element numbers are omitted for simplicity. 3.1Equations of motion Assume at this point that the beam is modelled using only one beam element (Fig. 1). Assuming that no external forces are applied except for the driving con- straint, the equations of motion of the fi nite element of the rotating beam in the absolute nodal coordinate formulation can be written as M ¨ e ? Qr? Fe(e)(5) where M is the element mass matrix, Qris the gener- alized force vector due to the driving constraint, and Feis the generalized elastic force vector. In order to write the equations of motion in the rotating frame, the following coordinate transformation is introduced r1 r2 ?? ? Ar r1 r2 ?? (6) Fig. 1One-element model of the rotating beam Finite element analysis of the geometric stiffening effect189 K01204# #IMechE 2005 Proc. IMechE Vol. 219 Part K: J. Multi-body Dynamics at PENNSYLVANIA STATE UNIV on September 18, 2016pik.sagepub.comDownloaded from where r1and r2, as shown in Fig. 1, are the coordi- nates in the global frame,r1andr2are coordinates in the rotating frame, and Ar? cosu?sinu sinucosu ?? (7) whereuis the angle of rotation of the beam with respect to the global frame. Equations (6) and (1) allow for writing the relation- ship between the two sets of nodal coordinates. The fi rst set is the coordinates ek of node k defi ned in the inertial fi xed frame, while the second set is the vector qk defi ned in the rotating frame. The transformation is written as follows ek 1 ek 2 ek 3 ek 4 2 6 6 6 4 3 7 7 7 5 ? Ar0 0Ar ?? qk 1 qk 2 qk 3 qk 4 2 6 6 6 6 4 3 7 7 7 7 5 (8) where qk 1 ?rk 1 qk 2 ?rk 2 qk 3 ? @rk 1 @x qk 4 ? @rk 2 @x (9) are the coordinates of node k (k ? 1, 2) of the element expressed in the rotating frame. It is clear that the relationship between the two sets of element coordinates is as follows e ? Aq(10) where e ? e1e2. e8 ??T and q ? q1q2. q8 ??T, in which eit4(k?1)? ek i and qit4(k?1)? qk i with k ? 1, 2 and i ? 1,.,4. The transformation matrix is written as follows A ? Ar000 0Ar00 00Ar0 000Ar 2 6 6 4 3 7 7 5 (11) Differentiating equation (10) twice with respect to time, substituting into the equations of motion of the element (equation (5)), and then premultiplying by AT, it follows that M¨ q t C˙ q ?v2Mq ? ? Qr? Fe(q)(12) wherevis the angular velocity of the beam reference, C ? 2ATM˙A is a skew-symmetric matrix, and ? Qr? Px