FINITE-ELEMENT MULTIBODY PROCEDURES
Carlo L. Bottasso, LorenzoTrainelli
DipartimentodiIngegneria Aerospaziale,PolitecnicodiMilano, Milano, Italy
Pierre Abdel-Nour, GianlucaLabo
Agusta S.p.A., Cascina Costadi Samarate, Varese, Italy
Abstract
Thepresentworkfocusesontheapplicationofthe
nite-elementmultibodytechniquetothedynamic
analysisoftiltrotors. Adetailedmultibodymodel
is rstdescribed andvalidated againstwell
estab-lished solution procedures. The non-linear model
developedherein includes afull description ofthe
control linkages and three possible realizations of
thegimbalmount. Thee ects ofthethree gimbal
joint design solutions on the stability of the
sys-temareanalyzedwiththehelpofnumerical
simu-lations.
Introduction
New rotorcraft con gurations can pose stringent
modeling requirements, noteasily met with
stan-dard simulation procedures. Forexample, tilt
ro-tors presenta numberof uniquefeatures that set
themapartfrom classicalhelicopterrotors: highly
twisted blades, large pitch excursions across the
ightenvelope,gimbalmountsthatallowthe
ap-ping of the rotor while at the same time try to
provideaconstanttransmissionoftheangular
ve-locity. Lackofdetailinsomeofthesecritical
com-ponents,orquestionable\equivalent"orlinearized
models, can seriouslyundermine the reliability of
the analysis and its nal e ectiveness in the
de-sign process. It is clear that the overall accuracy
delivered by an aeroelastic model will depend on
theaccuracyof itscomponents,such asthe
struc-turaldynamicsmodel,theaerodynamicmodel,the
possiblecontrolsmodel,thecouplingstrategy,etc.
Onthestructuraldynamicssideoftheproblem,it
is nowpossibleto usemorekinematically and
dy-namicallyconsistentnon-linearmodelsthanitwas
usuallydoneinthepast.
In fact, multibody formulations can deal with
complex exible mechanisms of arbitrary
topolo-gies. Using this approach, a given mechanism is
modeledthroughanidealizationprocessthat
iden-ti es the mechanism components from within a
largelibraryofelementsimplementedinasoftware
code. Each element provides a basic functional
buildingblock,forexamplearigidor exible
mem-ber, a hinge, a motor, etc. After assembling the
various elements,one canconstructavirtual
pro-totypeofthemechanismwiththerequiredlevelof
accuracy. Therefore,usingthistechnologyonecan
model acomplex systemsuch asa tilt rotorwith
ahigherlevelofdetailthanitispresentlypossible
usingconventionalindustrialtools.
Inref. [1],amultibodysimulationprocedurewas
proposed that is applicable to rotorcraft systems
and that provides a comprehensive, modular and
expandable simulation software. The same code
wasdemonstratedontheaeroelasticanalysisof
en-gageanddisengageoperationsin highwind
condi-tions for a detailed model of an articulated rotor
in ref. [4] and ref. [2]. The multibody dynamics
analysisiscastwithin theframeworkofnon-linear
niteelementmethods,andtheelementlibrary
in-cludesrigidand deformablebodiesaswellasjoint
elements. Deformablebodiesaremodeledwiththe
niteelementmethod,incontrastwiththeclassical
approachthatpredominantlyreliesonrigidbodies
orintroduces exibilityby means ofa modal
rep-resentation.
details of the control linkages and of the gimbal
mount. The systemhereconsideredis theAgusta
Ericaconcept,anon-goingprojectforaninnovative
tiltrotorinthe10ton,20seatclass. At rst,a
val-idationiscarefullycarriedoutbycomparisonwith
establishedindustrial codes,includingNASTRAN
[6]andCAMRAD/JA[5]. Thevalidatione ort
in-cludes modal analyses at di erent rotor speeds in
aircraft and helicopter mode, and comparisons of
therotoraerodynamicloadsin di erent ight
con-ditions.
Next,weusethehigh delitykinematicand
dy-namic modeling capability of the multibody
ap-proach for studying the in uence of the possible
e ects of the geometric non-linearities related to
the hub design on the whirl- utter stability. We
discussthreepossiblemodelsofthegimbalmount.
First, we consider asimple mount representedby
astandarduniversal (Hooke) joint, which is nota
constant speed joint. Next, we study a more
in-teresting joint that guarantees an exact constant
speedtransmissiononlyincertainoperating
condi-tions. This joint isagain modeled asamultibody
system, using avariety of joints and rigid bodies.
Finally,weinvestigatean\ideal"jointthatrealizes
an exact constant speed transmission in all ight
conditions. Thee ectsofthesevariousrealizations
ofthegimbalmountson utterlimitsarediscussed
indetail.
Weconcludethepaperwithaplan offuture
ac-tivities, that include in particular the analysis of
the drive train torsional loads and of the overall
vibratorylevel.
Overview of the Finite Element
Multibody Code
Thebasicfeaturesof the niteelementmultibody
codeused in this work are brie y reviewedin the
following.
Theelementlibraryincludesthebasicstructural
elements such as rigid bodies, composite capable
beamsand shells, and joint models. All elements
are referred to a single inertial frame, and hence,
arbitrarilylargedisplacementsand niterotations
aretreatedexactly. Inthisformulation,nomodal
reductionisperformed, andthe full nite element
basedcross-sectionalanalysisprocedureisfully
in-tegrated with the multibody dynamics code and
ensures the ability to model components made of
laminatedcompositematerials.
Joints are modeled through the use of
appro-priate holonomicor non-holonomicconstraints
en-forcedbymeansofLagrangemultipliers. Alljoints
are formulated with the explicit de nition of the
relative joint motion asadditional unknown
vari-ables. This allows the introduction of generic
spring and/or damper elements in the joints, as
usually required for the modeling of realistic
con- gurations. Furthermore,thetimehistoriesofjoint
relativemotionscanbedrivenaccordingtosuitably
speci edtimefunctions.Alljointscanbeequipped
withbacklash,freeplayandfrictionmodels.
The code implements implicit integration
pro-cedures that are non-linearly unconditionally
sta-ble. Theproofofnon-linearunconditionalstability
stems from two physical characteristics of
multi-body systems that are re ected in the numerical
scheme at the discrete level: the preservation of
the total mechanical energy and the vanishing of
the work performed by constraintforces.
Numer-ical dissipationis obtainedby letting the solution
drift from the constantenergy manifold in a
con-trolledmannerinsuchawaythatateachtimestep,
energycanbedissipatedbutnotcreated. More
de-tails on these non-linearly stable schemes can be
foundin thebibliographyofref. [1].
Onceamultibodyrepresentationofasystemhas
beende ned,severaltypesof analysescan be
per-formed on the virtual prototype. A static
analy-sis solvesthe static equations of theproblem,
ob-tainedbysetting alltime derivativesto zero. The
deformedcon gurationofthesystemunderthe
ap-plied static loads is then computed. The static
loads can be of various types, such as prescribed
static loads, steady aerodynamic loads, orthe
in-ertial loads associated with prescribed rigid body
motions.
Oncethestaticsolutionhasbeenfound,the
dy-namic behavior of small amplitude perturbations
about this equilibrium con guration can be
stud-ied. This is doneby rst linearizing the dynamic
equations of motion,then extracting the
eigenval-uesandeigenvectorsoftheresultinglinearsystem.
Finally, static analysis is also useful for providing
anal-A dynamic analysis solves the non-linear
equa-tions of motion for the complete nite element
multibodysystem. Theinitialconditionsaretaken
tobeatrest,orthosecorrespondingtoapreviously
determined static or dynamic equilibrium
con g-uration. Automated time step size adaptivity is
available toincreasetheeÆciencyandaccuracyof
thesimulation.
An importantaspectof theaeroelasticresponse
ofrotorcraftvehiclesisthepotentialpresenceof
in-stabilitieswhichcanoccurbothonthegroundand
in ight. Themultibody code implementsthe
im-plicitFloquetanalysismethod[3],whichevaluates
the dominant eigenvalues of thetransition matrix
using the Arnoldi algorithm, without the explicit
computation of the same matrix, which is
poten-tiallyveryexpensiveforadirect niteelement
ap-proach. This method is ideallysuited for systems
involvingalargenumberofdegreesoffreedom.
Finally,variousvisualizationandpost-processing
procedures,includinganimationsand timehistory
plots,areusedtohelptheanalystduringthemodel
preparationphaseandfortheinterpretationofthe
computedresults.
Tilt Rotor Multibody Model
A topological viewof themultibody model ofthe
aircraft is symbolically given in g. 1. The
vari-ous mechanical components of the system are
as-sociated withthe elements foundin the libraryof
the code. The model consist of the wing, engine
andnacelle andofadetailedrepresentationofthe
rotor. The rotoritself includes shaft, swashplate,
pitch-links, constant speed joint, ex-beam, cu
andblade.
Theswashplateisrepresentedbytworigid
bod-ies,rotatingandnon-rotating,connectedbya
rev-olute joint. The non-rotating lower swashplate is
connected to the nacelle by means of a universal
joint followedby aprismatic joint. The collective
inputsignalisprovidedbyprescribingtherelative
displacement of the prismatic joint, whereas the
cyclicinputsignalsareprovidedbyprescribingthe
tworelativerotationsoftheuniversaljoint. The
ro-tatingupperswashplateisconnectedto the
pitch-links by means of universal joints. In turn, the
pitch-linksareattached tothe pitch-horn through
Figure 1: Topologicalrepresentationofthetilt
ro-tormultibody model. Onesinglebladeisdepicted
forclarity.
sphericaljoints andtransferthe inputcontrol
sig-nalsfromtheupperswashplateto theblades. The
rotationoftheupperswashplateisenforcedby
scis-sorsthat connectto therotatingshaft.
For this project, it was deemed necessary to
model the exibility of the control linkages in
or-der to include this e ect in the global aeroelastic
model. At rst, the exibility measured at the
head of the pitch-links was computed through a
detailed nite element model of thewhole control
system. Ascommonlyfoundin thiscase,the
com-puted equivalent sti ness at thislocation depends
on the type of load condition considered, namely
cyclic, collective or reactionless. This e ect was
reproduced by introducing a special arrangement
of additional joints with suitable sti ness
charac-teristics. In particular, a prismatic joint is
intro-ducedin thepitch-links,whose sti nessrepresents
the reactionless sti ness of the whole system. An
additional prismatic joint followed by a universal
joint is introduced immediately below the
swash-plate. The spring in this prismatic joint is
com-putedsothatitssti ness,inserieswiththesprings
ofthepitch-links,yieldsatotalsti nessofthe
sti -computedinordertoyield,againinserieswiththe
springs associated with the pitch-links, the cyclic
sti ness of thesystem. This way, the load
depen-dentsti nesscharacteristicsofthecontrollinkages
canbemodeledusingjustafewadditionaldegrees
of freedom. Furthermore, during the preliminary
design phasethe exactstructural sti ness
charac-teristicsofeachpartofthecontrolsystemare
usu-allynot known. This approachallowsthen to
ac-countforthecontrolsystemcompliance usingjust
afewmacro-parameters,ascommonlydoneinthis
phaseofthedesign.
Continuing in thedescription of themain
com-ponentsofthevirtualmodel,wenotethattheshaft
isattachedto thehub,modeledas amassiverigid
body,thatinturnconnectstoa ex-beam,followed
bythecu and nallybythebladeitself. Cu and
ex-beamareconnectedattwopointsthrough
elas-tomeric joints and can rotate one with respect to
theothertoallowthebladepitchsetting.
The wing model includes an inboard section
which connects to the fuselage, and an outboard
sectionthatcanbetiltedbyprescribingtherelative
rotation in a revolute joint that connects the two
parts. Rigid bodies ofappropriateinertial
charac-teristicsmodeltheengine,nacelleandgear-box.
The aerodynamic model is based on lifting
lines associated with each blade and parts of the
cu . The lifting lines are based on classical
two-dimensionalstriptheoryanduselocalpro le
aero-dynamiccharacteristics,accountingforthe
aerody-namiccentero set,twist,sweep,andunsteady
cor-rections. Thetwo-dimensionalaerodynamicmodel
iscorrectedthroughtheuseof thedynamic in ow
modeldevelopedbyPeters[7]. However,giventhe
factthatweareheremainlyinterestedinthe
whirl- utter speed boundaries in airplanemode, the
ef-fectsofanin owmodelonthecomputedresultsis
negligible.
The nal multibody nite element model
in-cludes about 7,000 degrees of freedom, a clearly
prohibitivesizeforclassicalFloquetanalysis.
Constant-Speed Joint Models
The constant speed joint that connects shaft and
hub wasrealized in three di erent versions. The
rstsolutionsimplyadoptsauniversaljoint,which
null rotor apping. The other two joints
stud-ied here represent an approximate and an ideal
constant-speed joint which are describedin detail
in thenexttwosections.
The \Artichoke" Constant-Speed
Joint Model
Figure 2: Topological model of the approximate
constantspeedjoint.
Theapproximateconstantspeedjointis
symbol-icallydepictedin g. 2. Thedrivingshaftandhub
are connectedby asphericaljoint(A) that allows
thehubcompleterotationalfreedom. Transmission
of motionbetweenshaft andhubis provided bya
numberofscissors. Eachscissorissymmetricabout
theABplane,whichisnormaltotheshaftfornull
apping. A sphericaljointisusedatB,while
con-nections to shaft and hub are realized with
revo-lute joints. Therotor bladesare connectedto the
drivenhubat theABplane. Clearly,from a
kine-matical point of view one single scissor would be
enough,andredundancyishereintroducedonlyto
distributetheresultingloadsamongalarger
num-berofstructuralelements. Sincethereexistsan
al-ternativeversionofthisjointusingcomposite-made
torque transfer petals that somehow resemble an
artichoke, this joint will becalled the "artichoke"
in thefollowing,forthesakeofbrevity.
This implementation of the constant speed link
betweenshaftandhubismodeledusingrigid
bod-ies,revoluteandsphericaljoints. Therefore itwas
easily introducedin the globalmodeldescribedin
the previous section, and used for the aeroelastic
simulations.
Asecondmultibodymodelofthejointisdepicted
equiv-plementationallowstotestthejointbehavior
inde-pendently of the rest of therotor. This was done
in order to tryto give apreliminary
characteriza-tion of the joint itself by means of numerical
ex-periments. Itiseasily shownanalyticallythatthis
joint will provide an exact constant speed
trans-missionatthehubforarbitraryconstanttiltingof
the hub normal with respect to a xed reference
frame, i.e. for constanttip plathplane; thisisthe
same result that could be obtained by using two
universal jointsin series. However, the analytical
characterizationofthejointbecomesquitecomplex
for generalmotions of thehub normal, and afew
numericalexperimentscan behelpful in clarifying
itsbasicproperties.
Thesphericaljointthatallowsthetiltingofthe
hubinthe rstimplementationisherereplacedby
a universal and a revolute joint in series, which
are then connected to the ground. Testing of the
model can then beconducted by explicitlytilting
thehubthroughtheprescribedrotationsinthe
uni-versal jointwhile spinning the shaft at some
con-stant speed, and measuring the resulting angular
speedsatthehub.
Figure 3: Topological model of the approximate
constantspeedjointusedforitscharacterization.
Wehaveconducted twotests with twodi erent
waysoftiltingthehub:inthe rsttest,thenormal
to the huboscillates in a plane that containsthe
shaft axis, while in the second test thenormal to
the hub describes a cone. The shaft is driven by
imposing atimehistoryto therelativerotationat
the revolute joint in E. Tilting of the hub is
ob-tainedbyprescribingthetimehistoriesofthe
rela-tiverotationsintheuniversaljoint. Notethatthe
points labeledA, Cand D in the gureare in
re-alitycoincidentintheactualmodel. Rotationsare
measuredandplottedat thedrivingrevolutejoint
Angular velocities are measured in body attached
axesandplotted forthedrivingshaftandthehub.
Example1: Oscillationsin a Plane
Forthis rst case,oneangle in theuniversal joint
isdescribedbyacosinefunction whiletheotheris
held xedandnull. Theoscillationamplitudeis20
deg. with nullmeanvalue, while its speedis four
timestheangularspeedofthedrivingshaft.
Fig. 4givesthe time history of the shaft
rota-tions, of the universal joint rotations, and of the
relative rotation in the revolute joint at point A.
Fig. 5givesthetime historyof thebodyattached
componentsofangularvelocityalongtheshaftaxis
and thehubnormal. It isseenthat thejoint
pro-videsanexactconstantspeedtransmissioninthis
case.
0
1
2
3
4
5
6
7
8
9
10
−100
0
100
200
300
400
500
600
Time [sec]
Relative rotations [deg]
phi driver
phi hub
phi unj 1
phi unj 2
Figure 4: Oscillations in a plane: shaft rotations
(solidlines), universaljointrotations(dash-dotted
lines).
Example2: ConeMotions
In this second case, afully three-dimensional
mo-tion of the hub is considered. In particular, the
normal to the hub plane exactly describes acone
of semi-aperture = 20 deg. In order to achieve
this result, one angle (
1
) in the universal joint
is described by a sine function while the other is
computed as
2
=acos(cos()cos(
1
)). Since the
number of blades for the Ericadesign is four, we
usehereaprecessionspeedwhichisfourtimesthe
0
1
2
3
4
5
6
7
8
9
10
−0.2
0
0.2
0.4
0.6
0.8
1
1.2
Time [sec]
Angular speeds about shafts [rad/sec]
omega
3
driver
omega
3
hub
Figure 5: Oscillations in a plane: driving shaft
angular speed ! d
3
(solid line), slaveshaft angular
speed! s
3
(dash-dotted line).
Fig. 6gives the time historyof theshaft
rota-tions, of the universal joint rotations, and of the
relativerotationin therevolute jointat A. Fig. 7
givesthetimehistoryofthebodyattached
compo-nents ofangular velocityalong the shaft axisand
thehubnormal. Itisseeninthiscasethatthejoint
doesnotprovideanexactconstantspeed
transmis-sion: thehubrotationalspeedpresentsan
oscillat-ingbehavior,anditsmeanvalueislowerthanthe
angularspeedofthedrivingshaft.
0
1
2
3
4
5
6
7
8
9
10
−100
0
100
200
300
400
500
600
Time [sec]
Relative rotations [deg]
phi driver
phi hub
phi unj 1
phi unj 2
Figure 6: Cone motions (20 deg. semi-aperture):
shaft rotations (solid lines), universal joint
rota-tions(dash-dottedlines).
Thetransmissionbecomesclosertotheconstant
speedcasewhentheconesemi-apertureisreduced.
This is shownif g. 8, that reports the time
his-toryoftheangularvelocitycomponentsforacone
0
1
2
3
4
5
6
7
8
9
10
0.7
0.75
0.8
0.85
0.9
0.95
1
1.05
Time [sec]
Angular speeds about shafts [rad/sec]
omega
3
driver
omega
3
hub
Figure 7: Cone motions (20 deg. semi-aperture):
driving shaft angular speed ! d
3
(solid line), slave
shaftangularspeed! s
3
(dash-dottedline).
semi-apertureof5deg.,avaluethatismore
repre-sentativeofthe apping motionsof arotor. Note
thatthehubspeedoscillationsaregreatlyreduced
in this case, and that the velocity mean value is
nowclosertothedrivingangularspeed.
0
1
2
3
4
5
6
7
8
9
10
0.7
0.75
0.8
0.85
0.9
0.95
1
1.05
Time [sec]
Angular speeds about shafts [rad/sec]
omega
3
driver
omega
3
hub
Figure 8: Cone motions (5 deg. semi-aperture):
driving shaft angular speed ! d
3
(solid line), slave
shaftangularspeed! s
3
(dash-dottedline).
Ideal Constant-Speed Joint Model
Thethirdand nalgimbalmountconsideredinthis
work is an ideal joint, that guarantees an exact
transmissionof thedriving shaft angularspeedto
thehubangularspeedaboutitslocalnormal. The
other twocomponents ofthe hub angularvelocity
areclearlynotconstrainedbythejoint,sothatthe
beenforcedasanon-holonomicconstraint,
equat-ingthetwocomponentsofangularvelocity. In
re-ality, the same e ect can also be obtained with a
simplemodi cationof thesystemtopology, as
de-picted in g. 9. In fact we can takeadvantageof
thefact thattheshaftisrigidandthatthedriving
angular speedis constant and known a priori. In
thismodi edmodel,theshaftnowonlyservesthe
purposeofdriving thescissorsthatconnecttothe
upperpartoftheswashplate. Thehubisconnected
tothenacellethroughauniversaljointthatallows
its apping motions. A revolute joint is now
in-troducedbetweenthehubandthisuniversaljoint.
The relative rotationin the revolute joint can be
drivenwithalinearintimefunctionsoasto
guar-anteethedesiredvalueofangularvelocity.
Figure 9: Topological model of theideal constant
speedjoint.
Model Validation
Before conducting theaeroelastic simulations,the
modelwasvalidatedagainstalternativesimulation
procedures.
Structural Validation
The structural validation was determined rst in
terms ofnaturalfrequencies ofthe critical system
components. Thereference code usedforthe
vali-dationisNASTRAN.TheNASTRANmodelofthe
multibody code, and includes a detailed
descrip-tions ofthe ex-beamand cu with multiple load
paths. Nastran Present 2.74 2.75 4.00 4.00 7.29 7.28 19.08 19.05 25.90 25.93 45.80 45.50
Table1: Wingeigenfrequencies[Hz].
Table 1 gives the rst eigenfrequencies of the
wing. Very close agreement is observed for all
modes considered. Table 2 gives the rst
non-rotatingeigenfrequenciesoftherotor. Eveninthis
case,verycloseagreementisobservedforallmodes
considered. Table3 reports the rotating
eigenfre-quenciesinavacuumfortherotorata
representa-tiveairplanemode ightcondition. Thecorrelation
between the twocodes is slightly less satisfactory
in thiscase,inparticularforthetorsionalmodes.
Nastran Present Cyclic gimbal 0.66 0.66 1chord 10.33 10.30 1tors 28.80 28.84 2tors 42.68 42.63 1beam 77.59 77.24 Collective 1beam 4.25 4.24 1chord 10.35 10.31 1tors 33.11 33.16 2tors 49.54 49.41 2beam 88.82 89.33
Table2: Non-rotatingeigenfrequenciesin airplane
mode[Hz].
Forthe validation e ort, aCAMRAD model of
the tilt rotorwas also prepared. This model uses
an equivalentbeam for modeling the whole blade
system, since multiple load paths are notallowed
Nastran Present Cyclic gimbal 6.23 6.20 1chord 11.08 11.04 1tors 28.99 25.90 2tors 45.36 44.56 2beam 79.77 79.26 Collective 1beam 7.79 7.71 1chord 11.10 11.05 1tors 33.10 30.60 2tors 51.90 50.91 2beam 89.95 89.72
Table 3: Rotating eigenfrequencies in airplane
mode[Hz].
propertiesintheareaofthe ex-beamandcu are
adjusted in order to yield eigenfrequenciessimilar
tothoseofthedetailed model.
Various testsat di erent rotationalspeedswere
conductedin ane orttotryto comparethe
mod-elingofthetennisracquete ectinCAMRADand
the presentcode. In the latter case, the
geomet-ricallyexactbeamtheory adopted doesnot
intro-ducemodeling approximationsforthis potentially
importante ect.
Fig. 10 shows the pitch distribution along the
blade span computed using the two codes. Only
thebladeis deformablein this case,and in nitely
rigid control linkages are used. A somewhat
dif-ferent blade twisting is observed. Fig. 11 shows
the results obtained by considering an in nitely
sti blade,butcompliantcontrollinkages. This
ef-fectis obtainedin CAMRADintroducing suitable
torsional springs in the pitch hinge, while forthe
multibody code it is given by the special
mecha-nism described earlier. Even in this casea
di er-entbehaviorisobserved,demonstratingadi erent
modeling of the tennis racquet e ect in the two
codes.
Theapproximationsin the tennis racquete ect
modeling in CAMRAD might imply a somewhat
di erentaerodynamicbehavioroftherotor,for
ex-ample in termsof thrustversuspower,becauseof
the di erent pitch settings of the blade sections.
This was indeed observed for a number of ight
conditionsatsealevel. Thee ect ishowevermuch
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
45
50
55
60
65
70
75
80
85
r/R
Collective [deg]
Camrad
Present
Figure10: Tennisracquete ect,rigidcontrol
link-ages and exible blade: pitch distribution along
span.
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
45
50
55
60
65
70
75
80
85
r/R
Collective [deg]
Camrad
Present
Figure 11: Tennis racquet e ect, exible control
linkages and sti blade: pitch distribution along
span.
lesspronouncedforthe ightconditionsat7500m
usedfortheaeroelasticstudiesreportedhere.
Aerodynamic Validation
The aerodynamic validation was conducted by
comparisonwiththeCAMRADcode. At rst,
var-ious trim conditions in airplane mode were
com-puted using CAMRAD for the deformable blades
and control linkages model. The nacelle was
clampedto the groundduring these tests, so that
the rotor operates in axial ow. Next, the same
pitchcontrolsettingscomputedbyCAMRADwere
appliedtothemultibodymodel. Thesepitchvalues
swash-plate along the shaft. Finally, the inertial loads
correspondingto therotationabouttheshaftaxis
wereappliedtothesystem,togetherwiththe
corre-spondingaerodynamicloads. Thede ected
con g-urationof thesystem was consequentlycomputed
usingthestaticsolutionprocedure.
Fig. 12 shows the values of thrust and shaft
power computed by the two codes at 7500 m of
altitude. Very close agreement is observed at all
ightconditionsconsidered,exceptforthelastone
where the di erences in the modeling of the
ten-nis racquete ectinduce slightchangesin angleof
attackofthebladesections. Athighvaluesof
dy-namic pressure, small changes in angle of attack
can produce changes in global aerodynamic loads
ontherotorwhich arenottotallyneglectable.
0
2
4
6
8
10
12
14
16
18
x 10
5
1000
2000
3000
4000
5000
6000
7000
8000
Power [W]
Thrust [N]
Camrad
Present
Figure 12: Airplane mode: thrust vs. power at
di erenttrimpoints.
Tocheckthebehaviorofthetwocodesfor ight
conditionswhicharenotinaxial ow,anadditional
test wasconducted. In thiscase the rotorsystem
isfullyrigidand appingatthegimbaljointis
pre-vented. The rotor is then actuated at constant
speed at various nacelle tilt values. The global
forces and moments on the rotor were then
mea-suredandcomparedbetweenthetwocodes. Good
correlation was observed even in this case for all
measured loads at all nacelle tilt angles. For
ex-ample, g. 13 givestheH force component
(aero-dynamic resultant component in the plane of the
rotor)versusnacelletilt.
70
75
80
85
90
95
100
105
110
−5
−4
−3
−2
−1
0
1
2
3
4
5
x 10
4
Tilt angle [deg]
H Force [N]
Camrad
Present
Figure 13: Quasi-static conversion: rotor H force
vs. nacelletiltangle.
Whirl-Flutter Analysis
Analysis Procedure
The utteranalysiswasconductedaccordingtothe
followingprocess.
At rst, astatic solution procedure applies the
pitch control settings for the ight condition
un-derconsiderationbytranslatingtheswashplateas
previously discussed. Next,the procedure applies
the inertial loads correspondingto the rigid
rota-tionabouttheshaftandtheresultingaerodynamic
loads. Duringthis sequenceof operations, the
ro-tor is allowed to deform under the applied loads;
however, the wing tip is partially clamped to the
ground. Moreprecisely, chord and twistrotations
areprevented,sincethesewouldcausetherotorto
leavetheaxial owcondition.
The equilibrium con gurationcomputed in this
wayisnowusedastheinitialconditionfora
tran-sient dynamic analysis. At the beginning of this
process,thewingtipclampisremoved. Underthe
rotor thrust, the wingthen de ects and the rotor
leavesthe axial owcondition. The simulationis
thenadvancedintimeforafewcompleterotor
rev-olutions.
Ifthe ightconditionanalyzedisstable,the
sys-tem willquicklyreachaperiodicsolution,
charac-terizedbymoderate appingmotionsoftherotor.
If onthe other hand the ight conditionis
unsta-ble,thesystemwillquicklydivergefromtheinitial
condition. Inthis sense,theinitialcondition
simplewayofexcitingthesystemtotestits
stabil-ity.
Analternativewayoftestingthestabilityofthe
system is by using Floquet analysis. Using this
approach,theperiodicequilibriumsolutionshould
rst be computed. One has then to evaluate the
dominant eigenvaluesof the thetransition matrix
that maps initial perturbations of the various
de-grees of freedom of the systeminto perturbations
afteronerotorrevolution,perturbationscomputed
starting from the periodic solution. The Arnoldi
algorithmprovidesawaytocomputethis
informa-tion even for systems denoted by a considerably
largenumberofdegreesof freedom.
Whenthesystemisstable,theperiodicsolution
computedaspreviouslydiscussedprovidesthe
ini-tialcondition for theFloquet analysis. Whenthe
systemisunstable,weusethefollowingwayof
com-puting the periodic solution. First,the axial ow
conditionwith partially clampedwing tip is
com-puted as previously explained. Next, a transient
analysisisstartedfromtheresultssoobtained,
re-movingthewingtip clamp,asbefore. However,in
this case a certain amount of structural damping
is introduced in the wing. The damping is
pro-gressivelyincreaseduntil onecanreach aperiodic
solution. Theintroductionofdampingisclearlyan
arti cialdevicethathastheonlypurposeof
stabi-lizingthesystem. However,itse ectsonthe
peri-odicsolutionobtainedisminimal,sinceitwillonly
altertheperiodicde ectionsofthewing,whichare
howeververysmallattrim. Notethatnodamping
isintroducedintherotor,sincethismightnotably
changetheperiodiccondition.Oncethisarti cially
stabilizedperiodicconditionhasbeenobtained,the
structuraldampinginthewingisremovedandthe
Floquetanalysis isconducted.
Results
The whirl- utteranalysis processdescribedso far
wasapplied totheEricavirtualmultibody model.
A typicalsolutionobtainedwiththe proposed
ap-proach is shown in g. 14. The plot shows the
spectral radius of the dominant eigenvalue versus
ight speed at 7500 m of altitude. The systemis
unstable if the spectral radius is larger then one.
Fromtheplot,itappearsquiteclearlythatthe
de-tailsofthetransmissionofmotionformtheshaftto
thehubhaveanoticeablee ectonthestabilityof
250
300
350
400
0.8
1
1.2
1.4
1.6
1.8
2
Speed [kts]
Spectral radius
Universal joint
The "artichoke"
Ideal constant speed joint
Figure 14: Whirl-Flutter: spectral radius of the
transition matrix vs. ight speed for the three
joints.
the aircraft. In particular, the universal joint has
the lowest utter speed, while the ideal joint has
the highest. The \artichoke" gimbalmount, that
approximates quite well an ideal constant speed
transmissionfor smallrotor apping asshown
be-fore,behavesbetterthantheuniversaljoint,indeed
quitesimilarlytotheidealone.
Togain a better insight into the mechanism of
utter, the eigenvectorassociatedwith the
unsta-bleeigenmodecanbeanalyzed. Fig. 15showstwo
snap-shotsfrom ananimationof themode. It
ap-pears that the instability is caused by a coupling
ofthewing rstbendingandtorsionalmodeswith
laggingand appingmotionsoftherotor.
Inordertotrytocorroboratetheseresults,
tran-sient analyses were conducted starting from the
axial ow, partially clamped conditions, as
de-scribedabove. Fig. 16showsafewsnap-shots
ob-tainedfromthetransientsimulationinanunstable
regime. Theanimationshowsonceagainthe
qual-itativenatureoftheinstability,withlargebending
and torsion of the wing coupled with lag motions
ofthebladesandpronounced appingofthewhole
rotor.
Fig. 17showsthecomputedtimehistoriesofthe
wing tip displacements for the ight condition at
350 kts. Even usingthis approach to analyze the
system stability, it appears that the three joints
behavein markedly di erent ways. In particular,
theuniversaljointisinamoreunstableregimethan
then the \artichoke" joint, which in turn is more
Figure 15: Whirl-Flutter: animationof the
domi-nanteigenmodeofthetransitionmatrix( ight
con-ditionV=350kts) fortheuniversaljoint.
joint appears to be very near the stability limit,
sincethewingtipamplitudes growveryslowly.
Comparing with g. 14, it should be noticed
that the Floquet analysis had predicted a stable
solution for the ideal joint at 350 kts, prediction
whichis notherecon rmed. There areafew
pos-sible explanations for this fact. First, when
per-forming the transientanalysis, the excitation was
appliedtothesystemnotatitsperiodicsolutionas
in theFloquetcase,but startingfrom the
approx-imate trim in axial ow. Second, this excitation
is not guaranteed to be \small", so some system
non-linearitiesmightbeexcitedduringtheprocess.
Third, eventhe implicitFloquetmethod needs to
perturbthesysteminordertoextractthedominant
eigenvalues. Iftheperturbationistoosmall,its
ef-fectscouldbemaskedduringtheintegrationofthe
equationsofmotionthroughoutacomplete
revolu-tion. In fact, the integration is clearly conducted
in nite precision arithmetic and requires the
so-lution ofnon-linearsystemsofalgebraicequations
ateachtimestep,whichnecessarilyimpliestheuse
ofsomeconvergencetolerance. Ontheotherhand,
sientresponseintheunstableregime( ight
condi-tionV=350kts) fortheuniversaljointsolution.
if the perturbation is too large, it can excite the
non-linearities present in the system. While due
attention was paid to all these issues during the
analysis,itisclearthatthesetopicsneedadditional
investigationinthefuture.
The transientsimulationswere also analyzedin
order to try to understand thewhirl- utter
insta-bility herefound. Thegloballoads applied bythe
rotorat the wingtip provide someinsightonthis
problem. Fig. 18 shows the time history of the
wingtip forces. Itappearsthat themain
contrib-utorto theinstability istheloadpromotingbeam
de ections of the wing. Thechord loadoscillates
aboutameanvaluethatisincloseagreementwith
thenominalthrustatthis ightcondition( g. 12).
Althoughthegimbalmounthasasmallamountof
apping sti nessinthisrotor,themomentsatthe
wing tips are almost essentially due to the
trans-port of the rotor forces. Even though the three
jointsstudiedareclearlydi erentasfarasthe
sta-bility boundariesareconcerned,themechanismof
utter in the three cases does not seem to di er
100
100.5
101
101.5
102
102.5
−0.1
−0.08
−0.06
−0.04
−0.02
0
0.02
0.04
0.06
0.08
0.1
Time [sec]
Wing tip displacements [m]
Axial
Chord
Beam
100
100.5
101
101.5
102
102.5
−0.1
−0.08
−0.06
−0.04
−0.02
0
0.02
0.04
0.06
0.08
0.1
Time [sec]
Wing tip displacements [m]
Axial
Chord
Beam
100
100.5
101
101.5
102
102.5
−0.1
−0.08
−0.06
−0.04
−0.02
0
0.02
0.04
0.06
0.08
0.1
Time [sec]
Wing tip displacements [m]
Axial
Chord
Beam
Figure 17: Whirl-Flutter: wingtip de ections vs.
timeintheunstableregime( ightconditionV=350
kts)forthethree joints(top tobottom: universal,
artichoke,ideal).
100
100.5
101
101.5
102
102.5
−3
−2
−1
0
1
2
x 10
4
Time [sec]
Wing tip forces [N]
Axial
Chord
Beam
100
100.5
101
101.5
102
102.5
−3
−2
−1
0
1
2
x 10
4
Time [sec]
Wing tip forces [N]
Axial
Chord
Beam
100
100.5
101
101.5
102
102.5
−3
−2
−1
0
1
2
x 10
4
Time [sec]
Wing tip forces [N]
Axial
Chord
Beam
Figure18: Whirl-Flutter: wingtiptotalinternal
re-actionforcesvs. timeintheunstableregime( ight
condition V=350 kts) for the three joints (top to
Wehavestudiedthree alternativedesignsolutions
forthegimbalmount. The rstusesasimple
uni-versal joint, and is representative of the level of
detail allowed by other non-multibody based
sim-ulationprocedures. Thesecondis anapproximate
constantspeedjoint,andisrepresentativeofa
pos-sible actual hardware implementation of this
crit-ical component. The third is an ideal joint that
guarantees perfect constant speed transmission in
all ightconditions,butwouldclearlybenoteasy
torealizeandimplementin arealaircraft.
Using these three gimbal mounts, wehave
con-ducted a preliminary study of the whirl- utter
boundaries. Stability itself was analyzed in two
di erent ways, by means of transient simulations
and by using the implicit Floquet method. Both
approachesseemto indicateaprogressiveincrease
inthe utterspeedassociatedwithamoreaccurate
transmission ofthe constant speed of theshaft to
thehub. Thebestresultswerealwaysachievedby
theidealjoint,eventhoughttherealizable
approx-imateconstantspeed joint,heretermedthe
\arti-choke",wasseentoimply onlymodestreductions
onthe utterboundaries.
We plan to use the tilt rotor model developed
herein for a number of further studies. First, a
betterinsightonthestabilityofthesystemis
nec-essary. We intend to study the system response
toperturbationsofvariableamplitude,in orderto
assess the stability \robustness". This study will
beconducted perturbing thesystemstartingfrom
the periodic trim conditions, by applying suitable
forceexcitationsat thewingtip, similarlytowhat
is donein actualexperimental settings. Second, a
partfrom considerationsonwhirl- utter stability,
amoreconstanttransmissionoftheangularspeed
fromtheshafttothehubcouldalsohavean
impor-tantimpactonthedrivetrainloads andvibratory
levelsoftheaircraft. Weintendtoinvestigatethese
e ects by conducting transient simulationsfor
ro-torsoperatingat someangleofattack,in orderto
excitethefullrotordynamics,andbythen
analyz-ingthesystemresponse inthefrequencydomain.
Theresultsthatwereherereportedaretobe
con-sidered preliminary,and further investigationsare
surelynecessarybeforeassessingthetrueimpactof
thedesigndetailsoftheconstantspeedjointonthe
aeroelasticcharacteristicsofatiltrotor.
Nonethe-inarystudythatmodelingassumptionsand
simpli- cationsintheanalysisofthesemachinesincertain
ightconditionsmightseverelyunderminethe
ac-curacyofthecomputedresults. Inthissense,the
-niteelementmultibodyapproachseemstoo erthe
potentialforenhancedmodeling ofcomplex ying
machines,bysimplyprovidingthetoolsforadirect
numericalsimulationofthesystemcomponents.
Acknowledgments
The rst author acknowledges the help of Prof.
OlivierA. Bauchau, GeorgiaInstitute of T
echnol-ogy, for the many years of fruitful collaboration
onmultibody dynamics. The rst twoauthors
ac-knowledgethesupportofAgustathroughcontract
01/138/VILwiththePolitecnicodi Milano.
References
[1] O.A. Bauchau, C.L. Bottasso and Y.G.
Nik-ishkov. Modeling rotorcraft dynamics with
nite element multibody procedures. Math.
Comput.Modeling, 33:1113{1137,2001.
[2] O.A. Bauchau, J. Rodriguez and C.L.
Bot-tasso. Modeling of unilateral contact
condi-tionswithapplicationtoaerospacesystems
in-volvingbacklash,freeplay andfriction. Mech.
Res. Comm., 28:571{599,2001.
[3] O.A. Bauchau and Y.G. Nikishkov. An
im-plicit Floquet analysis for rotorcraft stability
evaluation.J. A.H.S.,46:200{209,2001.
[4] C.L. Bottasso and O.A. Bauchau. Multibody
modeling of engage and disengage operations
of helicopter rotors.J. A. H.S., 46:290{300,
2001.
[5] W.Johnson.CAMRAD/JA:Acomprehensive
analytical model of rotorcraft aerodynamics
and dynamics. Johnson Aeronautics Version.
VolumeI:Theorymanual.1988.
[6] MSC/NASTRANVersion70.5User'sGuide.
[7] D.A.PetersandC.J.He.Finitestateinduced
owmodels.PartII:Three-dimensionalrotor
