Calcul des caractéristiques aérodynamiques d’un rotor

en mouvement de translation rectiligne et rotation

uniformes

Jean Fourcade <autogire@[Link]>

18 janvier 2016
Table des matières

1 Modélisation théorique 7
1.1 Définition des repères . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Calcul des forces aérodynamiques . . . . . . . . . . . . . . . . . . 11
1.3 Calcul des moments de battement et de trainée . . . . . . . . . . 17
1.4 Calcul du moment de torsion . . . . . . . . . . . . . . . . . . . . 18
1.4.1 Calcul du moment de torsion aérodynamique . . . . . . . 18
1.4.2 Calcul du moment de torsion dû à la force centrifuge . . . 19
1.4.3 Calcul du moment dû à l’inertie de torsion . . . . . . . . 21
1.5 Tableau récapitulatif des forces et moments . . . . . . . . . . . . 22
1.6 Intégration des forces et moments . . . . . . . . . . . . . . . . . . 22
1.7 Etude du battement . . . . . . . . . . . . . . . . . . . . . . . . . 23
1.8 Equation de déformation de torsion . . . . . . . . . . . . . . . . . 27
1.9 Resultats de calculs . . . . . . . . . . . . . . . . . . . . . . . . . 28
1.10 Expressions des forces par rapport au disque rotor . . . . . . . . 30
1.11 Calcul de la polaire rotorique . . . . . . . . . . . . . . . . . . . . 31
1.12 Calcul de la vitesse induite . . . . . . . . . . . . . . . . . . . . . 31
1.13 Codage des équations . . . . . . . . . . . . . . . . . . . . . . . . 37
1.14 Revue des hypothèses . . . . . . . . . . . . . . . . . . . . . . . . 38

2 Caractéristiques remarquables d’un rotor en autorotation 40

2.1 Variation de la vitesse à incidence constante . . . . . . . . . . . . 40
2.2 Etude du coefficient de perméabilité . . . . . . . . . . . . . . . . 41
2.3 Fonctionnement en mode autogire . . . . . . . . . . . . . . . . . 42

3 Applications numériques 43
3.1 Test de validation des fonctions analytiques . . . . . . . . . . . . 43
3.2 Test de validation en mode autogire . . . . . . . . . . . . . . . . 44
3.3 Test de validation en mode hélicoptère . . . . . . . . . . . . . . . 53
3.4 Etude des coefficients de portance et de trainée . . . . . . . . . . 56
3.5 Evolution des forces et moments en fonction de l’azimut . . . . . 59

A Résultats Maple 63

1
Introduction
Ce rapport a pour objet le calcul des paramètres aérodynamiques d’un rotor
dont le mouvement se caractérise par une vitesse rectiligne uniforme et une
rotation uniforme. Ce document comporte trois parties et une annexe :
• la première partie traite des équations de la mécanique qui conduisent
aux expressions analytiques des forces, moments et coefficients de bat-
tement du rotor ainsi qu’aux équations de calcul de la déformation de
torsion des pales ;
• la deuxième partie expose les propriétés remarquables d’un rotor en ré-
gime établi d’autorotation ;
• la troisième partie donnes quelques exemples d’applications numériques
dont un cas de validation sur des mesures réelles d’un vol.
Les intégrales des équations conduisant aux formules analytiques exposées
dans ce document ont été calculées avec le logiciel de calcul formel Maple dont
les résultats sont fournis en annexe.
Les développements mathématiques sont principalement basés sur les tra-
vaux de H. Glauert (référence [1]) et John B. Wheatley (références [2], [3] et
[4]).
Les équations développées dans cette note ont été codées dans la biblio-
thèque Java GyroKit et ont conduit à la réalisation du programme GyroRotor.

2
Nomenclature

αB incidence de la pale à la distance r du pied de pale

αS incidence de la vitesse par rapport au plan d’entrainement
λ̄ vitesse de descente normalisée
µ̄ vitesse d’avancement normalisée
v̄ vitesse induite normalisée
β angle de battement
δ coefficient de traı̂née du profil de la pale
η rapport définissant le mode de fonctionnement du rotor
γ nombre de Lock
p̂ vitesse de roulis normalisée (p̂ = p/ψ̇)
q̂ vitesse de tangage normalisée (q̂ = q/ψ̇)
λ facteur de perméabilité
λD coefficient de perméabilité dans le repère TD
P pale du rotor considéré dans son ensemble
TB repère lié à la pale
TD repère constitué à partir du plan formé par le saumon des pales
TF repère en translation et rotation uniforme
TI repère inertiel de référence
TP repère pale définit en un point P de celle-ci
TR repère rotor tournant
TS repère rotor
µ facteur d’avancement
µD coefficient d’avancement dans le repère TD
ν(x) déformation de torsion de la pale
ψ angle d’azimut de la pale
ρ densité atmosphérique

3
θ0 pas de la pale au pied de pale
θB angle de pas
θT W vrillage de la pale
~ B /TI )
Ω(T vecteur rotation instantané du repère TB par rapport au repère TI
~ S /TI )
Ω(T ~ F /TI )
rotation instantanée du repère rotor = Ω(T
~σ (S, P/TI ) moment cinétique de la pale P calculé au point S
~
D traı̂née du rotor
~
L portance du rotor
~
R résultante des forces aérodynamiques générée par le rotor
~ (P/TI )
U vitesse du point P de la pale
~ (S/TI )
U vitesse de translation du rotor
~
V vitesse de l’air par rapport à la pale
A surface du disque rotorique
a pente du coefficient de portance du profil de la pale (1/rad)
a0 angle de conicité du rotor
a1 angle de basculement longitudinal du rotor
a2 angle du deuxième ordre de basculement longitudinal du rotor
B facteur de perte de portance en bout de pale
b nombre de pales
b1 angle de basculement latéral du rotor
b2 angle du deuxième ordre de basculement latéral du rotor
c corde du profil de la pale
Cd coefficient de trainée du profil de la pale
Cl coefficient de portance du profil de la pale
CM coefficient de moment aérodynamique du profil exprimé au centre
de gravité
Cm coefficient de moment aérodynamique du profil exprimé au foyer
CQ coefficient du couple robotique total
CT coefficient de la force axiale rotorique
CHD coefficient de la force latérale totale exprimée dans le plan rotor
CHi coefficient de la force arrière induite
CHp coefficient de la force arrière de profil
CQi coefficient du couple rotorique induit

4
CQp coefficient du couple rotorique de profil
CY i coefficient de la force latérale induite
H résultante des forces aérodynamique suivant l’axe Y~S ; compté po-
sitivement quand elle est de sens opposée à la vitesse
h ~ P (m)
coordonnée du foyer aérodynamique suivant l’axe Z
HD force arrière dans le repère TD
I(S, P) tenseur d’inertie de la pale P calculé au point S
Iy moment quadratique de la section de pale
Iz moment quadratique de la section de pale
Ixx inertie de pas de la pale
Iyy inertie de battement de la pale
Izz inertie de traı̂née de la pale
KF facteur de normalisation des forces
KQ facteur de normalisation des moments
l ~P (m)
coordonnée du foyer aérodynamique suivant l’axe Y
M moment rotorique en battement
n facteur de charge
P un point de la pale
p vitesse angulaire de roulis
Q moment rotorique en azimut
q vitesse angulaire de tangage
R longueur de la pale
r abscisse du point P de la pale
r vitesse angulaire de lacet
s coefficient de plénitude du rotor
T ~ S ; composante
résultante des forces aérodynamique suivant l’axe Z
axiale
TD force axiale dans le repère TD
UP composante perpendiculaire de la vitesse air par rapport à la pale
UR composante radiale de la vitesse air par rapport à la pale
UT composante tangentielle de la vitesse air par rapport à la pale
v ~H
composante de la vitesse suivant Y
vi vitesse induite par les forces rotoriques

5
W poids de l’autogire
w ~H
composante de la vitesse suivant Z
x abscisse normalisée du point P (x = r/R)
xac position du foyer du profil par rapport au bord d’attaque ramené
à la corde
xcg position du centre de gravité du profil par rapport au bord d’at-
taque ramené à la corde
Y ~ S ; composante
résultante des forces aérodynamique suivant l’axe X
latérale compté positivement vers la pale avançante
YD force latérale dans le repère TD

6
Chapitre 1

Modélisation théorique

1.1 Définition des repères

Notons TI = (O, X ~I, Y
~I , Z
~ I ) le repère inertiel dont l’axe Z~ I constitue la
verticale locale descendante. Le point O est le centre de ce repère.
Notons TF = (S, X ~F, Y
~F , Z~ F ) un repère animé d’un mouvement de transla-
tion et de rotation uniformes par rapport au repère inertiel TI . Le point S est
l’intersection entre le moyeu et l’axe de battement du rotor. Le plan (X ~F, Y
~F )
matérialise le plan d’entrainement du rotor (plan perpendiculaire au moyeu).
L’axe X ~ F est orienté vers l’avant et l’axe Z ~ F vers le bas (voir figure 1.1). Ces
axes sont tels que la vitesse du point S appartient au plan (X ~F, Z~ F ).
On note Ω(T~ F /TI ) le vecteur vitesse angulaire de rotation du repère TF
par rapport au repère TI . Les composantes de ce vecteur dans le repère TF
s’expriment par :  
p
~
Ω(TF /TI ) =  q  (1.1)
r T
F

~ S, Y
Soit TS = (S, X ~S , Z
~ S ) le repère rotor proprement dit. L’axe Z ~ S constitue
le moyeu et est dirigé vers le haut. Le plan (X ~ S, Y
~S ) est le plan d’entrainement
du rotor. On passe du repère TF au repère TS par une transformation définie
par la matrice :  
0 1 0
MSF =  1 0 0  (1.2)
0 0 −1
~ F /TI ) =
Le repère TS étant fixe par rapport au repère TF , on en déduit que Ω(T
~ S /TI ). Les composantes de ce vecteur dans le repère TS s’expriment par :
Ω(T
 
q
~ S /TI ) =  p 
Ω(T (1.3)
−r T
F

Notons αS l’angle d’incidence de la vitesse de telle manière que cet angle

soit positif quand sa composante w est négative (voir la figure 1.1). On déduit :
w
αS = arctan(− ) (1.4)
u

7
 Figure 1.1 – Forces aérodynamiques

~ (S/TI ) s’écrit :
Le vecteur vitesse U
 
0
~ (S/TI ) =  U cos αS 
U (1.5)
−U sin αS T
I

On passe du repère rotor TS au repère rotor tournant TR = (S, X ~ R, Y

~R , Z
~R)
~
par une rotation d’angle ψ autour de l’axe ZS (voir la figure 1.2). L’angle ψ
constitue l’azimut de la pale. L’origine est telle que ψ vaut 180 degrés quand le
~ R est confondu avec Y
vecteur X ~S (la pale se situe en avant de l’aéronef).

Figure 1.2 – Repère tournant

Ainsi quand :
ψ ∈ [0 . . . π] la pale se trouve dans la région avançante (la vitesse de trans-
lation du rotor s’additionne à la vitesse propre de la pale) ;
ψ ∈ [π . . . 2π] la pale se trouve dans la région reculante (la vitesse de trans-
lation du rotor se soustrait à la vitesse propre de la pale).

8
 Cette rotation définit le changement de coordonnées suivant :
    
x sin ψ − cos ψ 0 x
 y  =  cos ψ sin ψ 0   y  (1.6)
z T 0 0 1 z T
R S

On passe du repère rotor tournant TR au repère TB = (S, X ~ B, Y

~B , Z
~ B ) ma-
~
térialisant la pale par une rotation d’angle β autour de l’axe −YR (voir la figure
1.3).

Figure 1.3 – Repères pales

L’angle β constitue l’angle de battement. Il est orienté par l’axe −Y~R afin
d’être positif quand la pale située en avant de l’aéronef se lève. Nous ferons
l’hypotèse que l’angle de battement est un petit angle et nous limiterons son
developpement à l’ordre un.
La rotation de battement définit la transformation de coordonnées suivante :
    
x 1 0 β x
 y  =  0 1 0  y  (1.7)
z T −β 0 1 z T
B R

La composition des rotations en azimuth et en battement s’écrit alors :

    
x sin ψ − cos ψ β x
 y  =  cos ψ sin ψ 0  y  (1.8)
z T −β sin ψ β cos ψ 1 z T
B S

On déduit l’expression du vecteur vitesse dans TB :

 
−U cos αS cos ψ − U β sin αS
~ =
U U cos αS sin ψ  (1.9)
U β cos αS cos ψ − U sin αS T B

9
 Soit P un point de la pale situé sur l’axe X ~ B . On passe du repère TB au
~ ~ ~
repère local de la pale TP = (P, XP , YP , ZP ) par une rotation d’angle θB autour
de l’axe X~ B (voir la figure 1.3).
L’angle θB est l’angle de pas et X ~ B est l’axe de pas. L’angle de pas varie
le long de la pale et dépend donc de la coordonnée du point P . L’axe Y ~P est
par définition la ligne de portance nulle du profil. La rotation de pas définit la
transformation de coordonnées suivante :
    
x 1 0 0 x
 y  =  0 cos θB sin θB   y  (1.10)
z T 0 − sin θB cos θB z T
P B

On note I(S, P) le tenseur d’inertie de la pale P calculé au point S. On

suppose que le centre de gravité de chaque section de pale est situé sur l’axe
X~ B.
Considérons dans un premier temps que le pas est nul en tout point de la
pale ce qui revient à confondre les repères TP et TB . Nous supposerons que
le repère TP est le repère principal d’inertie de la pale. Le tenseur d’inertie
s’exprime alors sous la forme :
 
Ixx 0 0
I(S, P) =  0 Iyy 0  (1.11)
0 0 Izz T
P

Par construction :
• Ixx définit l’inertie de pas de la pale ;
• Iyy définit l’inertie de battement de la pale ;
• Izz définit l’inertie de traı̂née de la pale.
Des expressions :
Z
Ixx = (y 2 + z 2 )dm (1.12)
Z P

Iyy = (z 2 + x2 )dm (1.13)

P
On déduit : Z
Ixx + Iyy = (y 2 + 2z 2 + x2 )dm (1.14)
P
Etant donnée que la coordonnée z du profil de la pale est toujours petite
devant les coordonnées x et y, on peut écrire :
Z
Ixx + Iyy = (y 2 + x2 )dm = Izz (1.15)
P
On déduit finalement la forme du tenseur d’inertie :
 
Izz − Iyy 0 0
I(S, P) =  0 Iyy 0  (1.16)
0 0 Izz T
P

Nous considérerons que cette forme est encore valable quand le pas n’est
pas nul et que I(S, P) s’exprime de cette manière aussi bien dans TB que dans
TP

10
1.2 Calcul des forces aérodynamiques
On suppose le rotor composé de b pales de longueur R et de corde c. On
note A = πR2 la surface du disque rotorique et ρ la densité atmosphérique.
On note L ~ et D
~ respectivement la portance et la traı̂née générées par le
rotor.
~ générée par le rotor se décompose
La résultante des forces aérodynamiques R
traditionnellement en trois forces T, H et Y (voir la figure 1.1) telles que :
 
Y
~ =  −H 
R (1.17)
T T S

• T est la résultante suivant l’axe Z~ S (composante axiale) ;

• H est la résultante suivant l’axe Y~S (composante arrière compté positi-
vement dans le sens opposé à la vitesse) ;
• Y est la résultante suivant l’axe X ~ S (composante latérale compté posi-
tivement vers la pale avançante).
Exprimons le vecteur de rotation instantané de la pale. Il s’agit de la rotation
du repère TB par rapport au repère TI . Cette rotation est formée des vitesses
de rotation suivantes :
• q autour de X ~ S (mouvement de tangage) ;
• p autour de Y ~S (mouvement de roulis) ;
• −r autour de Z ~ S (mouvement de lacet) ;
• ψ̇ autour de Z ~ S (rotation de la pale) ;
• β̇ autour de −Y ~R (battement de la pale).
Nous pouvons donc écrire :
~ B /TI ) = q X
Ω(T ~ S + pY
~S + (ψ̇ − r)Z
~ S − β̇ Y
~R (1.18)

Exprimons les composantes de ce vecteur dans le repère pale. Nous obte-

nons :
    
sin ψ − cos ψ β q 0
~ B /TI ) =  cos ψ
Ω(T sin ψ 0  p  +  −β̇  (1.19)
−β sin ψ β cos ψ 1 ψ̇ − r T 0 T
S B

En considérant que les vitesses de rotation p, q, r sont petites devant la

vitesse de rotation du rotor ψ̇, on obtient :
 
β ψ̇ + q sin ψ − p cos ψ
~
Ω(TB /TI ) =  −β̇ + q cos ψ + p sin ψ
 
 (1.20)
ψ̇ TB

La vitesse de lacet n’intervient donc plus dans les équations. A partir de

maintenant nous noterons r l’abscisse du point P de la pale dans le repère TB :
−→ ~B
SP = rX (1.21)

11
 On déduit la vitesse du point P par rapport au repère TI :
 
0
~ (P/TI ) = Ω(T
U ~ B /TI ) ∧ rX~B =  rψ̇

 (1.22)
rβ̇ − rq cos ψ − rp sin ψ
TB

Nous pouvons à présent calculer la vitesse de l’air par rapport à la pale noté
~ afin de calculer les forces aérodynamiques. Celle-ci s’exprime par :
V
~ = −U
V ~ (S/TI ) − U
~ (P/TI ) − vi Z
~S (1.23)
vi est la vitesse induite communiquée aux molécules du fait de l’apparition
des forces rotoriques. Celle-ci, supposée constante sur toute la surface du disque
rotor, sera calculée au paragraphe 1.12.
On obtient en projetant dans le repère pale :
 
U cos αS cos ψ + U β sin αS
~ (P/TI ) = 
V  −U cos αS sin ψ − rψ̇

 (1.24)
−U β cos αS cos ψ + U sin αS − rβ̇ + rq cos ψ + rp sin ψ − vi
TB

Notons UR , UT , UP les trois composantes de ce vecteur dans le repère TB .

 
UR
~ =  −UT 
V (1.25)
UP T B

UT est la composante tangentielle (direction principale de la vitesse) compté

positivement quand le flux d’air souffle du bord d’attaque vers le bord de fuite,
UR est la composante radiale (composante le long de la pale) et UP la compo-
sante perpendiculaire au plan que décrit la pale. On déduit :
UP = −U β cos αS cos ψ + U sin αS − rβ̇ + rq cos ψ + rp sin ψ − vi (1.26)
UT = U cos αS sin ψ + rψ̇ (1.27)
On définit les facteurs suivants :
• le facteur d’avancement
U cos αS
µ= (1.28)
Rψ̇
• le facteur de perméabilité :
U sin αS − vi
λ= (1.29)
Rψ̇
Le facteur d’avancement caractérise le flux d’air perpendiculaire au mat
rotor ramené à la vitesse du bout de pale. Le facteur de perméabilité caractérise
le flux d’air qui traverse le rotor ramené à la vitesse du bout de pale.
Les deux équations ci-dessus permettent de calculer l’incidence du disque
rotor. On obtient :
λ vi
tan αS = + (1.30)
µ µRψ̇
On définit les quantités suivantes :

12
 r
• x= abscisse normalisé du point P ;
R
p
• p̂ = vitesse normalisée de roulis ;
ψ̇
q
• q̂ = vitesse normalisée de tangage.
ψ̇
En introduisant les coefficients µ et λ dans les expressions de UP et UR et
d d
en effectuant le changement de variable tel que = ψ̇ on obtient :
dt dψ

dβ
UP = Rψ̇(−µβ cos ψ + λ − x + xq̂ cos ψ + xp̂ sin ψ) (1.31)
dψ
UT = Rψ̇(µ sin ψ + x) (1.32)

L’expression de la vitesse tangentielle UT montre que quand µ sin ψ + x est

négatif, le fux d’air alimente la pale par le bord de fuite. Cette équation définit
une région dont la forme est un disque de rayon µR (voir la figure 1.4).

Figure 1.4 – Cercle d’inversion

Calculons à présent la portance dL~ et la traı̂née dD

~ élémentaires que déve-
loppe un segment de pale de longueur dr. Notons Cl et Cd respectivement le
coefficient de portance et de trainée du profil. Nous avons :

1
dL = ρcCl V 2 dr (1.33)
2
1
dD = ρcCd V 2 dr (1.34)
2

13
 Nous supposerons que le pas θB de la pale varie linéairement le long de la
pale et nous noterons ν(x) la déformation de torsion de la pale due aux forces
qui s’appliquent sur la pale. Le pas s’écrit donc :

θB = θ0 + xθT W + ν(x) (1.35)

L’angle d’incidence se calcule à partir du pas et de l’angle que prend le

vecteur vitesse par rapport à l’axe Y ~B . Nous négligerons la composante radiale
UR de la vitesse pour le calcul de cet angle.
Pour calculer l’incidence, il nous faut considérer deux cas : le cas nominal
où le flux d’air alimente le profil par le bord d’attaque (UT > 0) et le cas où le
flux d’air alimente le profil par le bord de fuite (UT < 0).
Le premier cas est représenté par la figure 1.5.

Figure 1.5 – Incidence de la pale : flux normal

L’incidence de la pale est l’angle αB = θB + φ avec :

UP
φ = arctan( ) (1.36)
UT
Etant donné que φ est un petit angle, on obtient au premier ordre :

UP
φ= (1.37)
UT
Soit :
UP
αB = θB + (1.38)
UT
La figure 1.6 montre le cas où la vitesse relative attaque le profil par le bord
de fuite. Dans ce cas l’angle φ de la figure s’exprime par :

UP
φ=− (1.39)
UT

14
 Figure 1.6 – Incidence de la pale : flux inversé

L’incidence s’exprime par αB = φ − θB , d’où l’on déduit :

UP
αB = −θB − (1.40)
UT
Notons δ le coefficient de traı̂née de la pale supposée constant. Notons a la
pente du coefficient de portance du profil par rapport à l’incidence également
supposé constant. On déduit l’expression des coefficients aérodynamiques du
profil :
Cl = aαB (1.41)
Cd = δ (1.42)
Soit :
1
dL = ρcaαB V 2 dr (1.43)
2
1
dD = ρcδV 2 dr (1.44)
2
Ayant ignoré UR et étant donné que UP est généralement faible devant UT ,
le module de la vitesse s’écrit :
V = UT (1.45)
Nous considérerons que le profil possède les mêmes coefficients aérodyna-
miques (a, δ) que celui-ci soit alimenté par le bord d’attaque ou le bord de
fuite. On déduit l’expression des forces aérodynamiques élémentaires :

1
 ρac(θB UT2 + UP UT )dr

 si UT > 0
dL = 2 (1.46)
1
− ρac(θB U 2 + UP UT )dr si UT < 0


T
2
1
dD= ρδcUT2 dr (1.47)
2

15
 En introduisant la coordonnée normalisée x du point P , on déduit :

1
 ρacR(θB UT2 + UP UT )dx

 si UT > 0
dL = 2 (1.48)
1
− ρacR(θB U 2 + UP UT )dx si UT < 0


T
2
1
dD= ρδcRUT2 dx (1.49)
2
Les forces élémentaires de portance et de traı̂née du rotor dans le cas de la
figure 1.5 (flux d’air abordant la pale par le bord d’attaque) s’expriment dans
TB en négligeant φdD devant dL par :
   
0 0
dR~ =  −dD + φdL  =  −dD + φdL  (1.50)
dL + φdD T
dL T
B B

Ces forces s’expriment dans le repère TS par :

 
−dD cos ψ + dL(−β sin ψ + φ cos ψ)
dR~ =  −dD sin ψ + dL(β cos ψ + φ sin ψ)  (1.51)
dL T S

En ramenant ces composantes à la définition des forces générées par le rotor

T, H et Y on obtient :

dT =dL (1.52)
dH=dD sin ψ − dL(β cos ψ + φ sin ψ) (1.53)
dY = − dD cos ψ + dL(−β sin ψ + φ cos ψ) (1.54)

Pour simplifier les calculs nous ne prendrons pas en compte l’effet d’inversion
du flux pour le calcul de dH et dY et considérerons que les équations ci-dessus
sont valables quel que soit le signe de UT .
On remarque que les équations qui fournissent les expressions des forces
H, Y font apparaı̂tre deux termes : un terme qui dépend de la traı̂née dD et un
autre qui dépend de la portance dL.
Nous distinguerons ces termes dans les calculs. Les termes dépendant de la
traı̂née seront indicés de la lettre P indiquant que leur origine provient direc-
tement du profil et les termes dépendant de la portance seront indicés par la
lettre I indiquant qu’ils sont induits par cette dernière.
En remplaçant dL et dD par leurs valeurs dans l’expression des forces ci-

16
dessus on obtient finalement :

1
 ρacR(θB UT2 + UP UT )dx

 si UT > 0
dT = 2 (1.55)
1
− ρacR(θB U 2 + UP UT )dx si UT < 0


T
2
1
dHP = ρδcRUT2 sin ψdx (1.56)
2
1
dHI = − ρacR((θB UT2 + UP UT )β cos ψ + (θB UP UT + UP2 ) sin ψ)dx (1.57)
2
1
dYP = − ρδcRUT2 cos ψdx (1.58)
2
1
dYI = ρacR(−(θB UT2 + UP UT )β sin ψ + (θB UP UT + UP2 ) cos ψ)dx (1.59)
2

1.3 Calcul des moments de battement et de trainée

Il s’agit de calculer les moments par rapport aux axes d’articulation du
rotor. Considérons dans un premier temps le cas où le flux d’air alimente le
profil par le bord d’attaque.
Le torseur des forces aérodynamique au point P se compose d’une résultante
~
dR et d’un moment de torsion suivant l’axe X ~ P . Intéressons nous pour l’instant
~
au moment dû uniquement à la résultante dR. Celle-ci s’écrit :
 
0
~ =  −dD + UP dL 
 
dR (1.60)
 U T

dL TB

Le moment de la force aérodynamique est alors :

 
  0  
x 0
−→ −→  UP 
 −dM 
SP ∧ dF = R  0  ∧   −dD + U dL  =
 (1.61)
0 T T dQ
B dL TB
T B

~B égal à :
Ce calcul conduit à un moment suivant l’axe Y

dM = RxdL (1.62)
~ R égal à :
Et un moment suivant l’axe Z

UP
dQ = Rx(−dD + dL) (1.63)
UT
Le premier moment est le moment de battement (il régit la variation de
l’angle β) et le deuxième moment celui de rotation proprement dit du rotor (il
régit la variation de l’angle ψ).

17
 Dans le cas où le flux d’air alimente le rotor par le bord de fuite, la résultante
~
dR s’exprime par :  
0
~ =  dD + UP dL 
 
dR (1.64)
 UT 
dL T B

On en déduit les moments :

dM = RxdL (1.65)
UP
dQ = Rx(dD + dL) (1.66)
UT
On obtient finalement en distinguant les termes de profil des termes induits :

1
− ρδcR2 xUT2 dx si UT > 0


dQP = 2 (1.67)
1
 ρδcR2 xU 2 dx

si UT < 0

T
2

1
 ρacR2 x(θB UP UT + UP2 )dx

 si UT > 0
dQI = 2 (1.68)
 1 2 2
− ρacR x(θB UP UT + U )dx si UT < 0

P
2

1
 ρacR2 x(θB UT2 + UP UT )dx

 si UT > 0
dM = 2 (1.69)
1
− ρacR2 x(θB U 2 + UP UT )dx si UT < 0


T
2

1.4 Calcul du moment de torsion

Calculons à présent le moment de torsion de la pale (moment suivant l’axe
X~ P ). La pale est soumise à trois moments de torsion qui sont :
• le moment dû aux forces aérodynamiques ;
• le moment dû à la force centrifuge ;
• le moment dû à l’inertie de rotation.
Pour simplifier, nous négligerons l’effet du cercle d’inversion. L’incidence
s’exprime donc quelque soit le signe de UT par

UP
αB = θB + (1.70)
UT

1.4.1 Calcul du moment de torsion aérodynamique

Le schéma de la figure 1.7 représente la position du foyer (AC) et du centre
de gravité (CG) du profil de la pale.
On note l et h les coordonnées du foyer dans le repère TP .

18
 Figure 1.7 – Position du foyer par rapport au centre de gravité

Notons Cm le coefficient de moment aérodynamique du profil exprimé au

foyer. Calculons le coefficient CM de ce même moment exprimé au centre de
gravité. Nous avons :

l h
CM = Cm + (Cl cos αB + Cd sin αB ) − (Cl sin αB − Cd cos αB ) (1.71)
c c
En remplacant Cl et Cd par leur valeur, en negligeant δ devant a et en
négligeant les terme du second ordre de l’incidence, on obtient :

la
CM = Cm + αB (1.72)
c
Le moment aérodynamique élémentaire s’écrit alors :

1
dN = ρc2 UT2 CM dr (1.73)
2
On déduit :
1 la la
dN = ρc2 [UT2 (Cm + θB ) + UP UT )]dr (1.74)
2 c c
Soit en se ramenant à la variable x :
1 la la
dN = ρc2 R[UT2 (Cm + θB ) + UP UT )]dx (1.75)
2 c c

1.4.2 Calcul du moment de torsion dû à la force centrifuge

La force centrifuge crée un moment de torsion qui a pour effet de s’opposer
à toute variation de pas : c’est le rappel à plat. La composante principale de
cette force étant due à la rotation du rotor (ψ̇), nous considérerons pour son
calcul que le battement est nul ainsi que les vitesse de roulis et de tangage.
Soit un point M de la pale défini dans le repère TB par ses coordonnées
r, y, z.
La figure 1.8 représente le point M dans le plan de rotation de la pale et
dans le plan du profil.

19
 Figure 1.8 – Force centrifuge sur un élement de pale

~ que subit se point se décompose en deux

La force centrifuge élémentaire dF
~ ~
forces élementaire dF1 et dF2 avec dF ~ 1 perpendiculaire à l’axe de pas.
Nous avons :
dF1 = dF sin ξ (1.76)
En notant dm la masse de cet élément de pale, on obtient :

dF1 = ψ̇ 2 r sin ξdm = ψ̇ 2 ydm (1.77)

Le moment de cette force s’écrit :

d2 P = −dF1 z = −ψ̇ 2 yzdm (1.78)

Le moment s’integre sur toute la section du profil. Il est préférable de changer

de repère et se placer dans le repère TP qui constitue les axes propres du tenseur
d’inertie.
En notant y 0 et z 0 les coordonnées du point dans le repère TP , on obtient :

y = y 0 cos θB − z 0 sin θB (1.79)

0 0
z = y sin θB + z cos θB (1.80)

Soit :

yz = (cos2 θB − sin2 θB )y 0 z 0 + sin θB cos θB (y 02 − z 02 ) (1.81)

En notant ρ la densité du matéiaux consituant la pale, la masse élementaire

s’exprime par :
dm = ρdrdy 0 dz 0 (1.82)

20
 L’intégrale du produit y 0 z 0 est nulle puisque le repère TP est un repère prin-
cipal d’inertie. En se limitant au premier ordre en θB , et en intégrant sur la
surface du profil on obtient :
Z
dP = −ψ̇ 2 θB dr (y 02 − z 02 )ρdy 0 dz 0 (1.83)
S

Notons Iz et Iy les moments quadratiques de la section de la pale. Nous

obtenons :
dP = −ψ̇ 2 θB (Iz − Iy )ρdr (1.84)
Nous avons :
Z
Iz = y 02 dy 0 dz 0 (1.85)
ZS
Iy = z 02 dy 0 dz 0 (1.86)
S

D’autre part les moments d’inerties s’écrivent :

Z
Izz = (y 02 + r2 )ρdrdy 0 dz 0 (1.87)
ZP
Iyy = (y 02 + r2 )ρdrdy 0 dz 0 (1.88)
P

On déduit :
Z Z R
02 02 0 0
Izz − Iyy = (y − z )ρdrdy dz = (Iz − Iy )ρdr (1.89)
P 0

Soit :
Izz − Iyy = (Iz − Iy )ρR (1.90)
En introduisant la variable x, on obtient finalement l’expression du moment :

dP = −ψ̇ 2 (Izz − Iyy )θB dx (1.91)

1.4.3 Calcul du moment dû à l’inertie de torsion

Le moment d’inertie de pas étant Izz − Iyy (voir équation (1.16)), le moment
d’inertie s’exprime par :
d2 θB
− (Izz − Iyy ) 2 (1.92)
dt
Nous verrons que ce moment est négligeable devant les deux précédents. Il
ne sera donc pas pris en compte.

21
1.5 Tableau récapitulatif des forces et moments
On obtient finalement le tableau récapitulatif suivant :

θB =θ0 + xθT W + ν(x) (1.93)


1
 ρacR(θB UT2 + UP UT )dx

 si UT > 0
dL = 2 (1.94)
1
− ρacR(θB U 2 + UP UT )dx si UT < 0


T
2
1
dHP = ρδcRUT2 sin ψdx (1.95)
2
1
dHI = − ρacR((θB UT2 + UP UT )β cos ψ + (θB UP UT + UP2 ) sin ψ)dx (1.96)
2
1
dYP = − ρδcRUT2 cos ψdx (1.97)
2
1
dYI = ρacR(−(θB UT2 + UP UT )β sin ψ + (θB UP UT + UP2 ) cos ψ)dx (1.98)
2
1
− ρδcR2 xUT2 dx si UT > 0


dQP = 2 (1.99)
1
 ρδcR2 xU 2 dx

si UT < 0

T
2

1
 ρacR2 x(θB UP UT + UP2 )dx

 si UT > 0
dQI = 2 (1.100)
1
− ρacR2 x(θB UP UT + U 2 )dx si UT < 0


P
2

1
 ρacR2 x(θB UT2 + UP UT )dx

 si UT > 0
dM = 2 (1.101)
1
− ρacR2 x(θB U 2 + UP UT )dx si UT < 0


T
2
1 la la
dN = ρc2 R[UT2 (Cm + θB ) + UP UT )]dx (1.102)
2 c c
dP = − ψ̇ 2 (Izz − Iyy )θB dx (1.103)

1.6 Intégration des forces et moments

L’expression analytique des forces et des moments des b pales s’obtient en
intégrant les expressions élémentaires des forces et moments (équations ci-dessus
(1.94) à (1.103)) une première fois le long de la pale et une deuxième fois en
azimut de 0 à 2π.
On introduit un paramètre noté B qui modélise les pertes de portance de
bout de pale du fait que la pale est de longueur finie. La valeur de B suggéré
par Prandtl est donné par la relation :

c
B =1− (1.104)
2R

22
 CT est le coefficient normalisé de la résultante rotorique axiale (voir para-
graphe 1.9). Dans la plupart des analyses ont prend pour simplifier la valeur
B = 0.97.
Pour calculer les forces et moments, on distingue les composantes de profil
et les composantes induites par la portance. On intègre les composantes induites
par la portance en prenant en compte le facteur de perte selon l’équation :
Z 2π Z B
b dFI
FI = dxdψ (1.105)
2π 0 0 dx
On intègre les composantes de traı̂nées en considérant toute l’étendue de la
pale car la traı̂née persiste même quand la portance disparait :
Z 2π Z 1
b dFP
FP = dxdψ (1.106)
2π 0 0 dx
Pour prendre en compte le changement de signe des équations sur l’étendue
du cercle d’inversion, on distingue dans l’intégration en azimuth la zone avan-
çante de la zone reculante. En prenant par exemple la résultante rotorique, on
obtient :
Z πZ B Z 2π Z −µ sin ψ Z B
b dT dT dT
T = ( dxdψ+ (− dx+ dx)dψ) (1.107)
2π 0 0 dx π 0 dx −µ sin ψ dx

Cette expression peut également s’écrire :

Z 2π Z B Z 2π Z −µ sin ψ
b dT dT
T = ( dxdψ − 2 dxdψ) (1.108)
2π 0 0 dx π 0 dx
Les expressions analytiques des forces et moments après intégration s’ex-
priment sous la forme de polynomes en µ. On montre (voir [2]) qu’une evalua-
tion à l’ordre 4 en fonction de µ suffit si a0 , a2 et b2 sont exprimés à l’ordre 2 et
a1 et b1 à l’ordre 3. La référence [3] montre que la même précision est obtenue
pour les forces et moments si u0 est évalué à l’ordre 4, u1 et v1 à l’ordre 3 et
u2 et v2 à l’ordre 2. Les forces et moments seront donc calculés en limitant le
developpement en fonction de µ à l’ordre 4 (voir annexe).

1.7 Etude du battement

Le calcul du battement consiste à exprimer l’équation différentielle qui régit
le mouvement en β . Pour cela calculons le moment cinétique au point S de la
pale P par rapport au repère TI . Nous avons :
~ B /TI )
~σ (S, P/TI ) = I(S, P)Ω(T (1.109)

Soit :
 
(Izz − Ixx )(β ψ̇ + q sin ψ − p cos ψ)
~σ (S, P/TI ) =  Iyy (−β̇ + q cos ψ + p sin ψ) (1.110)
 

Izz ψ̇ T B

23
 Calculons à présent le moment dynamique en S de la pale P par rapport au
repère TS . Celui-ci est égal à la dérivée du moment cinétique par rapport à ce
même repère. Le théorème de Varignon permet d’écrire :
 
d  d 
~σ (S, P/TI ) = ~σ (S, P/TI ) + Ω(T
  ~ B /TI ) ∧ ~σ (S, P/TI ) (1.111)
dt  dt 
TI TB

La coordonnée Y du premier terme du second membre de cette equation

donne :
− Iyy (β̈ + q ψ̇ sin ψ − pψ̇ cos ψ) (1.112)
La coordonnée Y du deuxième terme donne :
− Iyy (β ψ̇ 2 + q ψ̇ sin ψ − pψ̇ cos ψ) (1.113)
On obtient finalement la coordonnée Y du moment dynamique :
− Iyy (β̈ + β ψ̇ 2 + 2ψ̇q sin ψ − 2ψ̇p cos ψ) (1.114)
Le second principe fondamental de la dynamique stipule que le moment
dynamique est égal au moment des forces extérieures. En reliant cette dernière
équation à l’expression du moment aérodynamique précédemment calculés sur
~B et en l’intégrant le long de la pale nous obtenons l’équation différentielle
l’axe Y
du mouvement en battement :
Z B
Iyy (β̈ + β ψ̇ 2 + 2ψ̇q sin ψ − 2ψ̇p cos ψ) = dM (1.115)
0

Remarquons que le terme Iyy β̈ est la force d’inertie de la pale et le terme

Iyy β ψ̇ 2 la force de rappel centrifuge.
Pour intégrer β on utilise comme variable libre l’azimut au lieu du temps.
On obtient donc au final :
d2 β
Z B
1
+ β + 2q̂ sin ψ − 2p̂ cos ψ = dM (1.116)
dψ 2 Iyy ψ̇ 2 0
On introduit le nombre de Lock défini par :
ρacR4
γ= (1.117)
I
Ce nombre sans dimension caractérise le rapport des forces aérodynamiques
aux forces d’inerties. En remplaçant les vitesses UT , UP par leurs expressions et
en intégrant l’équation du moment dM le long de la pale on obtient l’équation
différentielle du battement.
L’expression de dM dépend comme nous l’avons vu du signe de UT . Nous
avons :
Z B Z B
1
dM = ρacR2 x(θB UT2 + UP UT )dx−
0 0 2
Z −µ sin ψ #2π
1
2 ρacR2 x(θB UT2 + UP UT )dx (1.118)
0 2
π

24
 L’expression ]2π π signifiant que ce terme n’est pris en compte que lorsque ψ
varie de π à 2π. Dans le but d’obtenir une seule expression du moment dans
l’intervalle [0 . . . 2π] il nous faut calculer les coefficient de Fourrier d’une fonction
définie sur deux intervalles, le premiere de [0 . . . π] et la deuxième de [π . . . 2π].
Soit une serie de Fourier définie de [0 . . . 2π] par :

y = A0 + A1 cos ψ + B1 sin ψ + A2 cos 2ψ + B2 sin ψ (1.119)

On additionne à cette série sur l’intervalle [π . . . 2π] une série de la forme :

∆y = ∆A0 + dA1 cos ψ + ∆B1 sin ψ + ∆A2 cos 2ψ + ∆B2 sin ψ (1.120)

L’analyse harmonique de la somme de ces deux séries conduit à une série

dont les coefficients sont (voir annexe) :

y + ∆y = C0 + C1 cos ψ + D1 sin ψ + C2 cos 2ψ + D2 sin ψ (1.121)

Avec :

∆a2 = − ∆a0 (1.122)

∆a1 =0 (1.123)
∆B1 1
C0 =A0 − + ∆A0 (1.124)
π 2
4 ∆B2
C1 =A1 − (1.125)
3 π
8 ∆A0 1
D1 =B1 − + ∆B1 (1.126)
3 π 2
1 2 ∆B1
C2 =A2 − ∆A0 + (1.127)
2 2 π
1
D2 =B2 + ∆B2 (1.128)
2
Les deux premières équations résultent de la continuité en π des deux séries.
L’expression analytique du moment de battement étant relativement longue,
elle ne sera pas détaillée dans ce chapitre. On la trouvera dans l’annexe.
Intéressons nous à présent à l’équation différentielle homogène du battement
(donc sans second membre). Celle-ci s’exprime dans le cas B = 1 :

d2 β γ 4 dβ γ 4
+ (1 + µ sin ψ) + [1 + ( µ cos ψ + µ2 sin 2ψ))]β = 0 (1.129)
dψ 2 8 3 dψ 8 3
On remarque qu’en plus des forces d’inerties, les forces aérodynamiques font
apparaı̂tre :
γ 4
• un terme d’amortissement : (1 + µ sin ψ) ;
8 3
γ 4
• une force de rappel : ( µ cos ψ + µ2 sin 2ψ)
8 3

25
 Cette équation est fondamentale car elle régit le mouvement de battement
et une attention toute particulière a été porté à sa stabilité. Cette équation
n’admet pas de solution générale mais les études montrent que des instabilités
apparaissent uniquement pour des valeurs de µ élevées.
Dans les domaines de vol courant le coefficient µ étant petit devant 1, il
peut être négligé. L’équation du battement se résume alors à :

d2 β γ dβ
+ +β =0 (1.130)
dψ 2 8 dψ
En revenant au temps comme variable libre au lieu de l’azimut, cette équa-
tion devient :
1 d2 β γ dβ
2
+ +β =0 (1.131)
ψ̇ dt
2 8ψ̇ dt
Il s’agit d’une équation différentielle du deuxième ordre aux coefficients
constants qui représente un oscillateur harmonique. Le facteur d’amortissement
γ γ
est ξ = et le temps de relaxation τ =
16 16ψ̇
La valeur de ce dernier dépend exclusivement du nombre de Lock et de
la vitesse de rotation. Les calculs donnent des valeurs de τ comprises entre
quelques centièmes de secondes pour des pales légères à un peu plus d’un dixième
de seconde pour des pales lourdes. Le mouvement de battement est donc un
mouvement fortement amorti.
Comme nous venons de le voir l’équation différentielle du battement n’admet
pas de solution générale. Pour calculer le battement nous allons le décomposer
en séries de Fourier en fonction de l’azimut. On pose :
n
X
β = a0 − (ai cos(iψ) − bi sin(iψ)) (1.132)
i=1

En se fixant un rang n et en introduisant l’expression de β ci-dessus dans

l’équation du battement, on obtient une expression discrète comportant 2n + 1
termes (n termes en cos(iψ) , n termes en sin(iψ) et le terme constant). Cette
équation devant être nulle quelque soit ψ , on en déduit 2n + 1 équations à
2n + 1 inconnus (les coefficients du battement).
Les calculs montrent que l’amplitude des coefficients diminue d’un facteur
10 par ordre de calcul et que les coefficients d’ordre 2 sont de l’ordre du degré.
Nous nous limiterons donc à l’ordre 2 pour développer les calculs dans cette
étude.
Les coefficients fondamentaux sont a0 , a1 , b1 dont les significations physiques
sont les suivantes :
• a0 est la conicité ;
• a1 est ele basculement longitudinal (positif quand le rotor bascule vers
l’arrière) ;
• b1 est le basculement latéral (positif quand le rotor bascule vers la pale
avançante).

26
 Nous avons vu que le moment de battement s’exprime à partir de l’équation
(1.115). En se limitant au premier ordre pour le calcul de l’angle β, soit :
β = a0 − a1 cos ψ − b1 sin ψ (1.133)
On obtient :
Z B
dM = I ψ̇ 2 (a0 + 2q̂ sin ψ − 2p̂ cos ψ) (1.134)
0
En l’absence de mouvement d’attitude (p̂ et q̂ nuls), on conclue que le bat-
tement au premier ordre est constant et vaut I ψ̇ 2 a0 .

1.8 Equation de déformation de torsion

Le calcul de la déformation s’obtient à partir de l’équation différentielle qui
régit la déformation de torsion :
dν M
= (1.135)
dr GJ
Dans cette equation ν est l’angle de torsion, GJ la rigidité de torsion de la
poutre et M la somme des moments des forces extérieures de r à R (convention
des efforts à droite).
En se ramenant à la variable x, on obtient pour la pale :
dν(x) RM
= (1.136)
dx GJ
On déduit l’expression de la déformation sachant que le moment M est la
somme des moments de torsion des forces aérodynamiques et inertielles (voir
paragraphe 1.4) :
Z x Z B Z 1
ν(x) = ( dN + dP )dx (1.137)
0 x x
Pour résoudre cette équation, on procède comme pour la résolution des
coefficients de battement en décomposant ν(x) en série de Fourier. On se limite
à l’ordre 2 pour être homogène avec l’équation de battement et on pose donc :
ν(x) = u0 (x) + u1 (x) cos ψ + v1 (x) sin ψ + u2 (x) cos ψ + v2 (x) sin ψ (1.138)
On exprime enfin les coefficients ci-dessus sous la forme de polynome en x
en se limitant à l’ordre 5, suffisamment précis en pratique (voir le paragraphe
1.8) :
u0 (x) = u01 x + u02 x2 + u03 x3 + u04 x4 + u05 x5 (1.139)
2 3 4 5
u1 (x) = u11 x + u12 x + u13 x + u14 x + u15 x (1.140)
2 3 4 5
v1 (x) = v11 x + v12 x + v13 x + v14 x + v15 x (1.141)
2 3 4 5
u2 (x) = u21 x + u22 x + u23 x + u24 x + u25 x (1.142)
v2 (x) = v21 x + v22 x2 + v23 x3 + v24 x4 + v25 x5 (1.143)
L’identification des membres de gauche et de droite de l’équation (1.137)
conduit à 25 équations dont les inconnues sont les coefficients de déformation.

27
1.9 Resultats de calculs
On introduit le coefficient de plénitude (ou solidité) défini par :

bcR
s= (1.144)
A
Ce coefficient exprime le rapport entre la surface des pales et la surface du
disque rotorique. On définit le coefficient de normalisation des forces par :

KF = ρsAR2 ψ̇ 2 = ρbcR3 ψ̇ 2 (1.145)

et celui des moments par :

KQ = ρsAR3 ψ̇ 2 = ρbcR4 ψ̇ 2 (1.146)

Les résultats d’intégration sont fournit en annexe. On ne précisera pas ici

les expressions dépendant de la déformation qui sont trop longues a développer.
La résolution analytique conduit au coefficients de battement suivant :

γ 1 5 µ4 1 1 µ3 1
a0 = {( µB 3 − )p̂ + ( B 3 + )λ + µ2 B 2 b2
2 6 48 π 3 4π 8
1 4 2 1 4 1 4 1 2 3 1 5
+ ( B µ − µ + B )θ0 + ( µ B + B )θT W } (1.147)
4 32 4 6 5
2µ B 4 5 8q̂ 1
a1 = {( + µ3 )p̂ − + (B 2 − µ2 )λ
1 2µ 48 γµ 4
B 4 − µ2 B 2
2
1 3 4 µ3
− B b2 + ( B 3 + )θ0 + B 4 θT W } (1.148)
3 3 3π
4µ 4p̂ 1 B4 1 1 1 µ3
b1 = {− − q̂ + B 3 a2 + ( B 3 + )a0 } (1.149)
1 γµ 4 µ 6 3 9 π
B 4 + µ2 B 2
2
γµ2 20 p̂ 128 1 7 q̂
a2 = 2 8 {− B 3 + (− − B γ)
γ B + 144 3 µ Bγ 3 µ
7 9 2 46 7 2 10
+ (16B + B γ )λ + ( B 2 + γ B )θ0
108 3 144
7 11 2
+( B γ + 12B 3 )θT W } (1.150)
180
γ 2 µ2 128 1 p̂ 20 B 3 q̂
b2 = − 2 8 {( 2 + B 7 ) −
γ B + 144 Bγ 3 µ 3 γ µ
5 5 25 6 8 7
+ B λ + B θ0 + B θT W } (1.151)
9 36 15
A partir des composantes des forces, on définit les coefficients suivant (on

28
remarquera que Yp étant nul, il n’y a pas de coefficient CY p ) :
T =KF CT (1.152)
Hp =KF CHp (1.153)
Hi =KF CHi (1.154)
Yi =KF CY i (1.155)
Qp =KQ CQp (1.156)
Qi =KQ CQi (1.157)
On obtient (voir annexe) :
a 1 1 1 1 1 1
CT = {(− µ3 + B 2 µ)p̂ + ( µ2 + B 2 )λ + µ3 a1 + Bµ2 b2
2 16 4 4 2 8 4
4 µ3 1 3 1 2 1 4 1 4 1 2 2
+ (− + B + Bµ )θ0 + ( B − µ + B µ )θT W } (1.158)
9π 3 2 4 32 4
1
CHp = δµ (1.159)
4
a 1 1 1 5 1
CHi = {(− B 2 λ − B 3 θ0 − B 4 θT W + B 3 b2 − B 2 µa1 )p̂
2 2 6 8 12 16
1 3 5 3 1 2
+ (− B a0 + B a2 − B µb1 )q̂
6 12 16
3 2 1 1 1
+ ( B a1 − Bθ0 µ − Bµb2 − B 2 µθT W )λ
4 2 4 4
1 3 1 1
+ ( B 3 a1 + B 2 µb2 )θ0 + ( B 4 a1 + B 3 µb2 )θT W
3 8 4 4
1 1
+ (− B 3 b1 − B 2 µa2 )a0
6 2
1 2 1 1 1
+ B µa02 + B 2 µa21 − B 3 b2 a1 + B 3 b1 a2 } (1.160)
4 4 4 4
a 1 3 5 3 5 2
CY i = {(− B a0 − B a2 + B µb1 )p̂
2 6 12 16
1 1 1 5 7
+ ( B 3 θ0 + B 4 θT W + B 2 λ + B 3 b2 + B 2 µa1 )q̂
6 8 2 12 16
3 2 1 3
+ ( B b1 + Bµa2 − Bµa0 )λ
4 4 2
1 3 1 2 3 3
+ ( B b1 + Bµ b1 − B 2 µa2 − B 2 µa0 )θ0
3 2 8 4
1 1 1 1
+ ( B 4 b1 + B 2 µ2 b1 − B 3 µa2 − B 3 µa0 )θT W
4 4 4 2
1 3 1
+ ( B a1 − Bµ2 a1 − B 2 µb2 )a0
6 2
1 3 1 2 1
+ ( B a2 + B µb1 )a1 + B 3 b1 b2 } (1.161)
4 4 4
1
CQp = δ{−8 − 8µ2 + µ4 } (1.162)
64

29
 a 5 1
CQi = {(− µ4 + B 4 )p̂2
2 64 8
4 µ4 θ0 1 4 1 1 1 1
+ (− + B µθT W + B 3 µθ0 + B 3 µb2 − B 4 a1 + µ3 λ)p̂
45 π 8 6 6 4 4
1 1 4
8 µ a0 1 3 1 1
+ (− µ4 + B 4 )q̂ 2 + ( + B µa2 − B 3 µa0 + B 4 b1 )q̂
64 8 45 π 6 3 4
1 2 1 2 2
+ ( B − µ )λ
2 4
1 4 1 1 1 2 µ3 θ0 3 3
+ ( µ θT W + B 2 µa1 + B 4 θT W + B 3 θ0 + − µ a1 )λ
32 2 4 3 9 π 8
1 1 1
+ (− B 2 µ2 a0 − B 3 µb1 + B 4 a2 )a2
4 6 2
1 2 1 1 1
+ ( B θ0 µ2 + B 3 µ2 θT W + B 3 µa1 + B 4 b2 )b2
8 12 6 2
1 2 2 1 4 2 1 4 3 2 2 2 1 1
+ ( B µ − µ )a0 + ( B + B µ )a1 + ( B 4 + B 2 µ2 )b21
4 16 8 16 8 16
1
− B 3 µa0 b1 } (1.163)
3

1.10 Expressions des forces par rapport au disque

rotor
~S
Les forces ont été exprimées juqu’à présent dans le repère TS dont l’axe Z
constitue le moyeu et le plan (X ~ S, Y
~S ) le plan d’entraı̂nement. Nous avons vu
que sous l’effet des forces aérodynamiques et inertielles le battement des pales
décrit un mouvement dont l’équation au premier ordre est :

β = a0 − a1 cos ψ − b1 sin ψ (1.164)

On déduit de ce mouvement que le rotor forme un cône d’angle a0 , bascule

longitudinalement en arrière d’un angle a1 en tournant autour de X ~ S et bascule
latéralement vers la pale avançante d’un angle b1 en tournant autour de Y ~S .
Il est alors intéressant de considérer un nouveau repère noté TD obtenu à
partir de TS et des rotations a1 et b1 . L’axe Z~ D de ce nouveau repère représente
l’axe du cône que forme le rotor et le plan (X ~ D, Y
~D ) est défini par le saumon
des pales. Ce plan est appelé plan du disque rotor.
La transformation de TS à TD s’écrit au premier ordre :
    
x 1 0 −b1 x
 y  = 0 1 a1   y  (1.165)
z T b1 −a1 1 z T
D S

Cette matrice permet de calculer l’expression des forces dans le repère TD :

    
YD 1 0 −b1 Y
 −HD  =  0 1 a1   −H  (1.166)
TD b1 −a1 1 T

30
 Soit :

TD = T (1.167)
HD = H − a1 T (1.168)
YD = Y − b1 T (1.169)

Les calculs montrent que les composantes HD , YD sont très faibles comparées
à T et que l’on peut donc considérer avec une très bonne approximation que la
résultante des forces aérodynamiques générée par le rotor est alignée avec l’axe
du cône formé par les pales.
L’incidence du disque rotorique s’écrit :

αD = αS + a1 (1.170)

Calculons à présent les coefficients d’avancement et de perméabilité dans ce

nouveau repère. Nous avons :

U cos(αS + a1 ) U cos αS
µD = = (1.171)
Rψ̇ Rψ̇
U cos(αS + a1 ) U sin αS + U a1 cos αS
λD = − vi = − vi (1.172)
Rψ̇ Rψ̇
Soit :

µD = µ (1.173)
λD = λ + µa1 (1.174)

1.11 Calcul de la polaire rotorique

Pour calculer les coefficients CD et CL du rotor il nous faut calculer la
portance L et la traı̂née D. Ces deux forces s’expriment directement en fonction
de T et H . Nous avons :

L = T cos αS − H sin αS (1.175)

D = T sin αS + H cos αS (1.176)

On déduit alors :
L D
CL = , CD =
1 1
ρAU 2 ρAU 2
2 2

1.12 Calcul de la vitesse induite

La vitesse induite est la vitesse communiquée aux molécules d’air du fait
de l’apparition des forces rotoriques. Pour calculer cette vitesse, nous ferons les
hypothèses suivantes :

31
 • la force générée par le rotor est uniformément distribuée à travers le
disque rotor. Cette distribution uniforme revient à considérer un nombre
infini de pale (la solidité vaut un) ;
• il y a continuité du flux d’air qui traverse le rotor ce qui conduit à
la continuité des vitesses des molécules d’air mais à l’apparition d’une
discontinuité de pression nécessaire à la génération de la force rotorique.
Nous considérerons donc que :
• la pression varie de p0 à p dans la partie aval (elle augmente) et de p+∆p
à p0 dans la partie amont (elle diminue) ;
• la vitesse de l’air par rapport au rotor varie de V ~ à V
~2 entre l’infini aval
~
et l’infini amont et vaut V1 au niveau du disque rotor.
Considérons le mouvement général de translation représenté par la figure
1.9.

Figure 1.9 – Flux en translation horizontale

On définit la vitesse induite par la relation :

~1 = V
V ~ + ~vi (1.177)

Appliquons le théorème de la quantité de mouvement. Nous avons :

~ = T~ + ṁV
ṁV ~2 (1.178)

Soit :
T~ = −ṁ(V
~2 − V
~) (1.179)
Calculons le bilan de puissance :

1 1
ṁV 2 = T~ .V
~1 + ṁV 2 (1.180)
2 2 2
Soit :
1
− T~ .(V
~ + ~vi ) = ṁ(V22 − V 2 ) (1.181)
2
En reportant l’expression de la force calculée précédemment nous obtenons :

32
 1
~2 − V
ṁ(V ~ ).(V
~ + ~vi ) = ṁ(V
~2 − V
~ ).(V
~2 + V
~) (1.182)
2
Soit :
~2 − V
V ~ = 2~vi (1.183)
La force rotorique T s’exprime alors par :

T~ = −2ṁ~vi (1.184)

La vitesse induite est alignée avec la force générée par le rotor et de sens
opposée. Il ne reste plus qu’à calculer le débit d’air. Glauert (voir note [1])
a suggéré qu’il ne faut pas considérer le flux d’air qui traverse le rotor mais
l’expression suivante :
ṁ = ρA|V~1 | = ρA|V~ + ~vi | (1.185)
Il n’y a égalité que dans le cas d’un mouvement vertical. On déduit finale-
ment les résultats : connaissant la force T~ générée par le rotor et la vitesse V
~ à
l’infini aval on a les relations :

T~ = −2ρA|V ~ + ~vi |~vi (1.186)

~1 = V
V ~ − ~vi (1.187)
~2 = V
V ~ − 2~vi (1.188)

La vitesse induite ainsi calculée est appelée vitesse induite de Froude.

~ S . Notons T tel que T~ = T Z
Notons ~vi tel que ~vi = −vi Z ~ S . La relation vec-
torielle donnant l’expression de la force se transforme en l’expression algébrique
suivante :
T = 2ρA|V ~ + ~vi |vi (1.189)
Soit VZ la composante de la vitesse de l’air suivant Z ~ S provenant du mou-
vement de translation de l’aéronef. La vitesse axiale des filet d’air au niveau du
disque rotor est :
VZ1 = VZ − vi (1.190)
A l’infini aval, nous avons :

VZ2 = VZ − 2vi (1.191)

Ces deux expressions montrent que la vitesse axiale varie de VZ à l’infini

aval à VZ − 2vi à l’infini amont en passant par VZ − vi au niveau du disque
rotor. On peut alors prédire que l’écoulement ne sera pas le même suivant le
sens et le module de VZ .
On considère quatre cas de fonctionnement :
• La descente rapide quand VZ > 2vi . Les filets d’air sont tous dirigées dans
le sens de l’axe Z ~ S . Ce régime correspond à celui du moulinet frein ;
• La descente modérée quand vi < VZ < 2vi . Il y a inversion du sens des
vitesses à l’amont du disque ce qui fait apparaı̂tre un sillage de sorte
que les filets au centre dirigés vers le haut se rabattent vers le bas à
l’extrémité des pales. Ce mode de fonctionnement porte le nom de régime
turbulent ;

33
 • La descente faible quand VZ < vi . Le sillage apparaı̂t cette fois à la
partie inférieure du rotor. Les filets d’air au centre du rotor se rabattent
vers le bas ce qui conduit les filets centraux supérieurs à créer une zone
tourbillonnaire d’où le nom d’état d’anneaux tourbillonnaire donné à ce
régime ;
• La montée quand VZ < 0. Il correspond au fonctionnement d’une hélice.
En autorotation, le rotor fonctionne soit en moulinet frein, soit en régime
turbulent. Nous montrerons au paragraphe 3.4 que le passe d’un régime à est
indépendant de la vitesse et ne dépend que de l’incidence.
Notons enfin que l’on définit le régime d’autorotation idéale celui dans lequel
VZ = vi ce qui conduit à VZ1 = 0.
Utilisons la vitesse de translation du rotor U~ qui, par définition, est l’opposé
~
de la vitesse V . Nous avons :

U ~S
~S − U sin αS Z
~ = U cos αS Y (1.192)

avec αS l’incidence du disque rotor positive dans le mode autogire. On pose :

Uy = U cos αS (1.193)
Uz = −U sin αS (1.194)

Posons vi0 la vitesse induite définie par :

s s
|T | sR2 ψ̇ 2 CT
vi0 = = (1.195)
2ρA 2
~ = 0) qui génère la force T .
vi0 est la vitesse induite du vol stationnaire (V
On introduit alors les coefficients suivants :
vi
v̄ = (1.196)
vi0
Uy
µ̄ = (1.197)
vi0
Uz
λ̄ = (1.198)
vi0
avec :
• v̄ la vitesse induite normalisée (positive vers le bas) ;
• µ̄ la vitesse d’avancement normalisée ;
• λ̄ la vitesse de descente normalisée (positive vers le haut).
λ̄ par définition est négatif quand Uz est négatif c’est-à-dire quand l’inci-
dence αS est positive (mode autogire). A l’opposé λ̄ est positif quand l’incidence
est négative (mode hélicoptère).
On pose :
λ̄ Uz
η= = (1.199)
v̄ vi
Ce paramètre détermine le mode de fonctionnement du rotor que nous avons
précédemment défini. En effet nous avons :

34
 Figure 1.10 – Mode autogire et helicoptère

• η̄ < −2 soit Vz > 2vi régime moulinet frein ;

• −2 < η̄ < −1 soit vi < Vz < 2vi régime turbulant ;
• −1 < η̄ < 0 soit 0 < Vz < vi état d’anneaux tourbillonnaires ;
• 0 < η̄ soit Vz < 0 mode hélice.
Le régime d’autorotation idéale est déterminé par η̄ = −1 soit la droite
λ̄ = −v̄. Le passage du régime turbulent au régime moulinet frein est défini par
la droite λ̄ = −2v̄.
Nous ne nous intéresserons dans la suite qu’au cas T ≥ 0 correspondant à
vi ≥ 0 soit v̄ ≥ 0. En effet la solution du cas T < 0 est −vi avec vi solution de
−T .
Exprimons le module du vecteur V ~ + ~vi . On obtient :
~ + ~vi |2 = vi0
|V 2
(µ̄2 + (λ̄ + v̄)2 ) (1.200)
L’expression de la force T de la théorie de Froude conduit à la relation :
q
1 = v̄ µ̄2 + (λ̄ + v̄)2 (1.201)
Cette expression ne donne pas de solution explicite. On peut cependant
tracer v̄ à µ̄ fixé en fonction de λ̄ en calculant λ̄ en fonction de v̄ ce qui conduit
à une équation du second degré. On obtient deux familles de solution :
s
1
λ̄ = −v̄ ± − µ̄2 (1.202)
v̄ 2
Ces courbes sont représentées sur le schéma ci-dessous pour trois valeurs
différentes de µ̄ . Le schéma ci-dessous fait également apparaı̂tre les zones de
fonctionnement en moulinet frein (intersection entre la valeur de v̄ et la pente
λ̄ = −2v̄) et en mode hélice (λ̄ > 0).
On remarque que pour µ̄ = 0 la solution de Froude correspond à deux
branches d’hyperboles admettant chacune la droite v̄ = −λ̄ comme asymptote.
Cette droite, qui correspond à la vitesse idéale d’autorotation, joue également un
rôle important quand µ̄ 6= 0 car les maxima des courbes v̄ = f (λ̄) se produisent
sur cette droite.
Dans le cas d’une descente verticale (µ̄ = 0) les résultats de la théorie de
Froude comparés aux mesures expérimentales permettent de déterminer trois
zones :

35
 Figure 1.11 – Vitesse induite de Froude

• lorsque λ̄ ≥ −1 la branche de l’asymptote de droite est assez en accord

avec les différentes expériences tout en minorants celle-ci ;
• lorsque −2 ≤ λ̄ ≤ −1 la vitesse induite diminue fortement jusqu’à re-
joindre l’asymptote de gauche ;
• lorsque λ̄ ≤ −2 les résultats expérimentaux suivent la branche de l’hy-
perbole en majorant celle-ci.
La théorie de Froude dans le cas d’un autogire ne peut donc plus s’appli-
quer dès que le fonctionnement du rotor s’approche du régime turbulent. Nous
montrerons que c’est le cas dès que l’incidence dépasse 30 degrés.
Des modèles plus réalistes ont été élaborés dont notamment le modèle de
Shaydakov basé sur la théorie des anneaux tourbillonnaires.
Ce modèle est abordé dans la note Ref. [4]. Il conduit à deux jeux d’équa-
tions : s
v̄
• quand η̄ > 1 alors 1 = v̄ µ̄2 + (λ̄ + )2 ;
2 s
p λ̄2
• lorsque η̄ < 1 alors 1 = (λ̄ + v̄) µ̄2 + (λ̄ + v̄)2 − λ̄ µ̄ +
4
On remarque encore une fois l’importance de la droite définissant le régime
d’autorotation idéale qui définit la transition entre les deux jeux d’équations.
Dans le cas η̄ > 1 l’équation n’admet pas de solution explicite et on procède
comme dans le cas de la théorie de Froudre pour tracer v̄ en fonction de λ̄ à µ̄

36
 Figure 1.12 – Vitesse induite de Shaydakov

fixé. Dans le cas η̄ > 1 on obtient :

v v
u u s
u
2
u − µ̄2 + µ̄4 + 4(1 + λ̄ µ̄2 + λ̄ )
u u
t
4
u
t
v̄ = −λ̄ + (1.203)
2
Le tracé de cette courbe montre que la valeur de v̄ obtenue pour η̄ = 1 est
la valeur minimale. Elle s’obtient (voir l’équation ci-dessus) pour :
s
λ̄2
1 + λ̄ µ̄2 + =0 (1.204)
4
Soit : q p
λ̄ = − 2( µ̄4 + 1 − µ̄2 (1.205)
Cette dernière équation définit la transition de manière explicite entre les
deux familles de courbes.
Le schéma de la figure 1.12 représente les modèles de Froude (en rouge) et
de Shaydakov (en bleu). Le modèle de Shaydakov permet de distinguer la zone
du régime turbulent de la zone d’anneaux tourbillonnaires.

1.13 Codage des équations

Les équations analytiques exposées dans ce document ont été codées dans
la bibliothèque GyroKit et donne lieu au logiciel GyroRotor qui calcule les
caractéristiques aérodynamiques d’un rotor dans lequel la vitesse de rotation
est imposée ou fonctionnant en auto-rotation.

37
 L’expression de la vitesse induite vi fait apparaı̂tre un couplage avec le co-
efficient de la force rotorique axiale CT . Pour une vitesse de rotation donnée du
rotor, on résout ce couplage en utilisant un algorithme itératif de type Newton
(recherche du zéro d’une fonction).
Une difficulté apparaı̂t lorsqu’on utilise le modèle de Froude pour de petits
facteurs d’avancement. Dans l’intervale λ̄ = −2 à λ̄ = −1, outre le fait que le
modèle de Froude est incorrect, il apparaı̂t une discontinuité dans la solution
quand on passe d’une branche à l’autre de l’hyperbole. Cette discontinuité peut
conduire à une non convergence de la force générée par le rotor. Dans cette
région, il est nécessaire d’utiliser le modèle de Shaydakov.
Dans la programmation des équations, la définition de vi0 est légèrement
différente de celle du paragraphe 1.7. Afin de toujours avoir v̄ ≥ 0, nous avons
pris : s
|T |
vi0 = signe(T ) 2ρ (1.206)
A
Dans le cas du fonctionnement du rotor en auto-rotation, il faut déterminer
de plus sa vitesse de rotation. Celle-ci est telle que le couple rotorique total
CQp + CQi est nul. On utilise pareillement un algorithme de type Newton pour
calculer cette vitesse de rotation.
On peut résumer ainsi l’algorithme de calcul dans le cas d’un fonctionnement
en auto-rotation :
• itération sur la vitesse de rotation du rotor jusqu’à annulation du couple
rotorique CQp + CQi ;
• itération sur la vitesse induite jusqu’à convergence de la force axiale CT ;
• calcul des coefficients de déformation ;
• calcul des coefficients de battement ;
• calcul de CT ;
• calcul de CQp + CQi .
Les coefficients de déformations sont calculés à partir de l’équation (1.137).
Cette équation se ramène à un système linéaire de 25 équations à 25 inconnues
qui est résolu numériquement.

1.14 Revue des hypothèses

Les principales hypothèses prises en compte lors du développement des équa-
tions et qui limitent la porté des calculs sont :
• la vitesse induite est supposée constante sur toute la surface du disque
rotor ;
• la zone de décrochage de la pale reculante n’est pas prise en compte ;
• le profil des pales est caractérisé par deux nombres constants : la traı̂née
δ et la pente du coefficient de portance a.
La première hypothèse a comme conséquence que les équation ne modé-
lisent pas correctement le comportement du rotor dans la zone des anneaux
tourbillonaires.

38
 La deuxième hypothèse concerne le calcul de la force de portane générée par
l’élément de pale qui est supposé proportionnelle à la pente du CL du profil.
Ce calcul n’est valide que tant que l’angle d’incidence est inférieur à l’angle de
décrochage du profil. En conséquence, le domaine de validité des calcul se réduit
aux faibles valeurs du facteur d’avancement.
La dernière hypothèse concerne les paramètres modélisant la traı̂née et la
pente de portance du profil. Ces paramètres doivent être considérés comme
des paramètres d’ajustement. On pourra ainsi les modifier afin d’obtenir une
portance et une vitesse de rotation données. En augmentant a on augmente
les forces aérodynamiques générées par le rotor en changeant peu la vitesse de
rotation. En augmentant δ on diminue la vitesse de rotation et l’amplitude des
forces aérodynamiques.

39
Chapitre 2

Caractéristiques remarquables
d’un rotor en autorotation

Dans ce chapitre nous négligerons la déformation de torsion des pales. Consi-

dérons un rotor aux caractéristiques géométriques données définies par les pa-
ramètres R, c, θ0 , θT W , B, Ib , Ic , a, δ.

2.1 Variation de la vitesse à incidence constante

Supposons que pour pour une vitesse U et une incidence du disque αS
données, les caractéristiques aérodynamiques de ce rotor soient :
ψ̇, vi , µ, λ, a0 , a1 , b1 , a2 , b2 , CT , CH , CY , T, H, Y (2.1)
Que deviennent ces valeurs pour une nouvelle vitesse à la même incidence ?
Notons kU cette nouvelle vitesse avec k un coefficient strictement positif. Mon-
trons que la nouvelle solution s’écrit :
k ψ̇, kvi , µ, λ, a0 , a1 , b1 , a2 , b2 , CT , CH , CY , k 2 T, k 2 H, k 2 Y (2.2)
Supposons que la vitesse de rotation du rotor soit multipliée par k.
L’expression de la vitesse induite vi0 (équation (1.195)) montre que celle-ci
est multipliée par k quand la vitesse de rotation du rotor est multipliée par k.
On en déduit à travers les équations (1.197) et (1.198) que les coefficients µ̄ et
λ̄ sont inchangés. Pour que l’équation (1.201) soit respectée, il faut donc que v̄
soit inchangé. On en déduit à partir de l’équation (1.201) que la vitesse induite
vi est multiplié par k.
Les expressions des paramètres µ, λ (equations (1.28) et (1.29)) montrent
que ceux-ci sont inchangés quand on multiplie les vitesses U, ψ̇, vi par un même
coefficient.
µ et λ étant inchangés, les equations (1.147), (1.148), (1.149), (1.150) et
(1.151) montrent que les coefficients de battements a0 , a1 , b1 , a2 et b2 restent
inchangés.
Il en est de même du couple rotor total CQ (équations (1.162) et (1.163)).
Celui-ci reste donc nul et ainsi la nouvelle solution respecte bien la condition
d’équilibre du régime d’autorotation.

40
 Enfin la vitesse de rotation étant en facteur au carré dans l’expression des
forces T, H, Y (voir équation (1.145)) on déduit que les valeurs de ces forces
sont multipliées par le coefficient k 2 .
Ainsi quand la vitesse de translation d’un rotor en autorotation
est multipliée par un coefficient k alors que l’incidence du disque est
inchangée :
• les coefficients µ, λ, a0 , a1 , b1 , a2 , b2 sont inchangés ;
• la vitesse du rotor ψ̇ et la vitesse induite vi sont multipliées par
k;
• les forces aérodynamique T, H, Y sont multipliées par k 2 .
Cette dernière propriété justifie la définition des coefficients CD , CL de la
polaire rotorique.
Si l’on examine à présent le rapport η introduit lors du calcul de la vitesse
induite :
U sin αS
η=− (2.3)
vi
on remarque que ce rapport est également indépendant de la vitesse.
Autrement dit, le mode de fonctionnement d’un rotor en auto-
rotation (moulinet frein ou régime turbulent) est indépendant de la
vitesse. Il ne dépend que de l’incidence.

2.2 Etude du coefficient de perméabilité

Considérons un rotor aux pales non vrillées, infiniment rigides et en transla-
tion pure. Pour étudier les caractéristiques du coefficient de perméabilité, nous
allons exprimer le couple rotorique sous une forme différente de celle de l’équa-
tion (1.163). Nous avons (voir la référence [6] page 102) :
Qi = −(T λ + Hi µ)R (2.4)
Cette equation qui pouvait être attendue montre que le couple rotorique
induit est le produit scalaire entre la résultant rotorique et le flux d’air qui
traverse le rotor.
On déduit de cette expression l’expression du coefficient total du couple
rotorique :
δ µ4
CQ = (1 + µ2 + ) − λCT − µCHi (2.5)
8 8
A partir de l’expression du coefficient CHp et en se plaçant dans le plan du
rotor, on déduit :
δ µ4
CQ = (1 + 3µ2 + ) − λD CT − µCHD (2.6)
8 8
Les termes µ2 et µ4 petit devant 1 peut être négligé. Le calcul montre que
CHD est négligeable devant CT étant donné que la force rotorique est prati-
quement alignée avec l’axe du cône rotor. La condition d’autorotation s’écrit
donc :
δ
− λD CT = 0; (2.7)
8

41
 A partir de l’expression de CT donné par l’équation (1.158), on obtient en
négligeant les puissance de µ devant 1 et se plaçant dans le plan du rotor :
a θ0 λD
CT = ( + ) (2.8)
2 3 2
La condition d’autorotation s’écrit finalement :
s
1 9δ
λD = ( θ02 + − θ0 ) (2.9)
3 2a
Le coefficient de perméabilité exprimé dans l’axe du cône d’un
rotor au profil des pales donnée est constant. Il ne dépend ni de la
vitesse, ni de l’incidence

2.3 Fonctionnement en mode autogire

L’expression de la portance du rotor est :
L = T cos αS (2.10)
En considérant que l’angle d’incidences est suffisamment petit de telle ma-
nière que cos αS = 1, on obtient :
L=T (2.11)
Nous avons :
T = ρbcR3 ψ̇ 2 CT (2.12)
Notons W le poids de l’autogire et n facteur de charge qu’il subit. Nous
avons par définition du facteur de charge :
T
n= (2.13)
W
A partir de la condition d’autorotation (equation (2.7)), on déduit :
nW
ψ̇ 2 = 8 λD (2.14)
ρR3 bcδ
Soit : s
8 nW 9δ
ψ̇ 2 = 3
[ θ02 + − θ0 ] (2.15)
3 ρR bcδ 2a
Ainsi pour de faibles valeurs d’incidence et au premier ordre, on déduit les
propriétés suivantes :
Pour un facteur de charge donné, le carré de la vitesse de rotation
du rotor est :
• indépendante de la vitesse et de l’incidence ;
• inversement proportionnelle à la densité, la corde, le nombre
de pale et le cube du rayon du disque.
La vitesse de rotation du rotor varie comme la racine carrée du fac-
teur de charge.

42
Chapitre 3

Applications numériques

3.1 Test de validation des fonctions analytiques

Un premier test de validation a consisté à vérifier le codage correct des
équations en comparant le résultat de la bibliothèque GyroKit et les calculs du
logiciel Maple. On utilise pour cela un rotor dont les caractéristiques sont les
suivantes :
• Longueur des pales, R 4m
• Masse de la pale, M 12 kg
• Corde du profil, c 0.2 m
• Centre de gravité du profil ramené à la corde, xac 0.3
• Foyer du profil ramené à la corde, xcg 0.278
• Nombre de pale, b 2
• Vrillage de la pale, θ1 2 deg
• Coefficient de moment aérodynamique du profil, Cm 0.005
• Rigidité de torsion de la pale, GJ 6350 N.m2 /rad
Le profil correspond à un NACA 8h12. La rigidité de torsion a été calculée
sur des pales AVERSO à partir des résultats de la mesure suivante : une force
de 100 N appliquée perpendiculairement à la pale à une distance d’un mètre à
3,81 m du pied de pale conduit à une flèche de 60 mm au point d’application
de la force. On déduit de ces données :
• Le moment de torsion 100 N.m
• L’angle de torsion 0.03 rad.
L’équation de torsion s’écrit :

dν(r) M
= (3.1)
dr GJ
Dans le cas d’un moment constant, celle-ci s’intègre en :

M
ν(r) = r (3.2)
GJ
d’où l’on déduit :
M
GJ = r (3.3)
ν(r)

43
 En remplaçant par les données de la mesure : r = 3, 81 m, M = 100 N.m,
ν(r) = 0.03 rad, on obtient la valeur de rigidité de torsion mentionnée.
A partir des données du rotor, on calcule les moments d’inerties en supposant
la masse linéique de la pale constante. On obtient :
R2
• Moment d’inertie de battement, Ib = m = 64 kg.m2
3
c2
• Moment d’inertie de traı̂née, Ic = Ib + m = 64.04 kg.m2
12
On simule alors le fonctionnement du rotor en auto-rotation dans les condi-
tions suivantes :
• Vitesse de translation, U 90 km/h
• Incidence, αS 7 deg
• Pas au pied de palle, θ0 2 deg
• Vitesse de tangage, q 2 deg/s
• Vitesse de roulis, q -3 deg/s
On obtient après convergence :
• Vitesse de rotation du rotor, ψ̇ 353.0027 rpm
• Facteur d’avancement, µ 0.1678124092350715
• Facteur de perméabilité λ, 0.013486848994919694
La saisie dans le logiciel Maple des caractéristiques du rotor ansi que les
données de son point de fonctionnement (ψ̇, µ, λ) permet de comparer les ré-
sultats entre le logiciel Maple et la bibliothèque GyroKit et vérifier ainsi la
cohérence entre le développement des équations et leurs codages (voir le code
[Link]).
On peut ensuite procéder à plusieurs tracés dans le logiciel Maple (voir les
figures à partir de la ligne 113) afin de confirmer certaines hypothèses de calcul.
Le premier tracé correspond aux moments de torsion dus aux forces aérody-
namiques, à la force centrifuge et à l’inertie de torsion pour différents azimuth
de la pale. On vérifie ainsi que l’inertie de torsion est négligeable devant les
deux premiers.
Les tracés qui suivent permettent de comparer les déformations calculées à
partir des coefficients de déformation et à celles calculées à partir du moment
de torsion ce qui permet de vérifier que l’ordre 5 du polynôme de déformation
est suffisant pour bien identifier la forme de la déformée.

3.2 Test de validation en mode autogire

La validation à partir de mesures en vol a été réalisé à partir du rapport
NACA No 600 ([4]) qui présente des mesures effectuées sur un autogire Kel-
lett KD-1. Ce rapport fournit sous forme de courbes les coefficients de bat-
tement a0, a1, b1, a2 et b2 ainsi que les coefficients de déformation de torsion
u0, u1, v1, u2 et v2 (calculés au 3/4 du rayon) en fonction du paramètre d’avan-
cement µ. Les valeurs mesurées ont été digitalisée avec le logiciel Engauge V4.1
(voir [5])
Les caractéristiques du KD-1 sont les suivantes :

44
 • Poids, W 2100 pounds
• Longueur des pales, R 20 feet
• Moment de battement de la pale, Ib , Ic 175 slug.feet2
• Corde du profil, c 1 foot
• Position du centre de gravité du profil 0.280 foot
• Position du foyer du profil 0.242 foot
• Nombre de pale, b 3
• Pas constant, θ0 0.096 radian
• Coefficient de moment aérodynamique du profil, Cm -0.056
• Rigidité de torsion de la pale, G 1700 [Link]
A partir des coefficients de conversion suivants :
• 1 slug = 14.5939029 kg
• 1 foot = 0.3048 m
• 1 pound = 0.45359237 kg
et en prenant 9.81 m/s2 pour valeur de l’accélération de la pesanteur, on
obtient les caractéristiques du KD-1 en unités SI :
• Poids, W 934.44 daN
• Longueur des pales, R 6.096 m
• Moment de battement de la pale, Ib 237.27 kg.m2
• Corde du profil, c 0.3048 m
• Centre de gravité du profil ramené à la corde, xcg 0.280
• Foyer du profil ramené à la corde, xac 0.242
• Nombre de pale, b 3
• Pas constant, θ0 0.096 radian
• Coefficient de moment aérodynamique du profil, Cm -0.056
• Rigidité de torsion de la pale, G 2305.67 N.m
La définition de la rigidité de torsion G dans la note [4] est différente de
celle du paragraphe 1.8 de ce document. Dans la note [4], G est le moment de
torsion qui s’applique au pied de pale, distribué le long de la pale suivant la loi
(B 2 − x2 ) et tel que la déformation au 3/4 du rayon vaille 0.8 radian. Calculons
ce moment. Les conditions aux limites conduisent à écrire le moment de torsion
sous la forme :
x2
M(x) = G(1 − 2 ) (3.4)
B
A partie de l’équation de torsion :
dν(x) RM(x)
= (3.5)
dx GJ
on déduit en intégrant :
GR x3
ν(x) = (x − ) (3.6)
GJ 3B 2
En prenant :
x = 0.75 (3.7)
ν(x) = 0.8 (3.8)
B = 0.97 (3.9)

45
 On obtient :
GJ = 0.751GR (3.10)
Soit : GJ = 10551 N.m2 /rad.
L’équation donnée dans la note [4] est :
585
GJ = GR (3.11)
64
avec G = 1700 [Link], R = 240 inches et GJ en pounds.inches2 . Le cal-
cul conduit à 10706 N.m2 /rad, une valeur légèrement supérieure mais cohérente
avec la précédente. Nous prendrons cette dernière valeur pour les calculs.
Le profil des pales du KD-1 est un Göttingeh 606 dont il nous faut disposer
des coefficients a et δ afin de simuler le comportant de l’autogire. On peut
calculer la valeur de a à partir du nombre de Lock donné dans le rapport [3].
Celui-ci, évalué au niveau de la mer, vaut : 12.74. On déduit :

γI
a= (3.12)
cρR4

En prenant ρ = 1.225 kg/m2 , on obtient : a = 5.86. Pour déterminer la

trainée de profil des pales, on procède comme suit : à partir du diagramme de la
vitesse de rotation du rotor en fonction de la vitesse de translation de l’autogire
donné dans le rapport [4] on calcule l’incidence et la traı̂née de profil de telle
manière à ce que la portance égale le poids et que le moment total de traı̂née
du rotor à la vitesse de rotation donnée soit nul (condition d’autorotation, voir
le programme [Link]).
On en déduit alors les coefficients de battement et les coefficients de la
déformés des pales évalués au trois quart du rayon que l’on peut comparer avec
les données mesurées en vols.

46
Figure 3.1 – Vitesse de rotation du rotor et coefficient de résultante rotorique

La figure 3.1 représente la vitesse de rotation du rotor (en rouge) et les

coefficient CT de la résultante rotorique (en bleu et vert). La courbe rouge
continue est obtenue par un ajustement polynomial à partir des points extraits
de la note [4]. Cette courbe interpolée est utilisé pour calculer l’incidence comme
indiqué ci-dessus. La courbe bleu est le coefficient CT calculé comme dans la
note [4]. Celle-ci passe parfaitement dans les points extraits de cette note.
La courbe verte est le coefficient CT calculé en prenant en compte l’inci-
dence. On remarque un écart entre la courbe bleue et verte pour les faibles
facteurs d’avancements où l’incidence du disque rotor est importante.

47
 Figure 3.2 – Incidence du rotor et traı̂née de profil de la pale

La figure 3.2 représente l’incidence du disque rotor et le coefficient de traı̂née

de profil de la pale. Le coefficient de traı̂née moyen du profil Gottingen 606 tel
que monté sur le KD1 est d’environ 0.0162.

48
 Figure 3.3 – Coefficients de déformation

La figure 3.3 représente les coefficients de deformation calculés au trois quart

du rayon de la pale. On remarque que les coefficients u0 et v1 qui sont les plus
importants sont surestimés. Dans la note [4] dont la déformée est un modèle
linéaire, le coefficient v1 est très bien estimé et le coefficient u0 est sous estimé.
Une augmentant de la rigidité de torsion conduit a diminuer v1. On peut
alors chercher la valeur qui permet un bon ajustement de ce paramètre. On
obtient : GJ = 12700 N.m2 /rad.

49
 Figure 3.4 – Coefficients de déformation avec GJ = 12700 N.m2 /rad

Le résultats est donné figure 3.4. On constate qu’avec cette valeur l’ensemble
des coefficients est parfaitement bien estimé.

50
 Figure 3.5 – Coefficients de battement

La figure 3.5 représente les coefficients de battement du rotor. L’ensemble

des coefficients exceptés b1 sont bien estimés. L’origine de erreur concernant le
coefficient b1 est connue. Elle provient du fait que la vitesse induite est supposée
constante dans le présent calcul or une variation de celle-ci augmentant de
l’avant à l’arrière du rotor accroı̂t le coefficient b1.
Il est intéressant de faire remarquer que la prise en compte de la torsion des
pales a peu d’effet sur les coefficients b1, b2, a2, une influence légère sur a0 et
une influence importante sur a1 (voir la note [3]).

51
 Figure 3.6 – Coefficients u0 et v1 en fonction du rayon pour mu=0.322

La figure 3.6 représente les deux coefficients de déformation les plus impor-
tants, u0 et v1 en fonction de rayon normalisé x de la pale pour µ = 0.32.
Il est intéressant pour terminer de tracer la distribution des incidences sur le
disque rotor. La figure 3.7 représente cette distribution obtenue avec le logiciel
GyroRotor pour µ=0.32. On y remarque l’importance de la zone de décrochage
des pales (incidence supérieure à environ 12 degrés), phénomène qui n’est pas
modélisée dans cette étude comme cela a été indiqué lors de la revue des hypo-
thèses (voir le chapitre 1.14).
Bien que les résultats soient conformes aux mesures, les simulations sont
conduites ici au-delà des limites de la modélisation.

52
 Figure 3.7 – Distribution des incidences des pales pour µ=0.32

3.3 Test de validation en mode hélicoptère

La validation s’est porté sur le robinson R22 et a consisté à comparer le
diagramme de puissance calculé a celui du manuel de vol (voir le programme
[Link]).
Les pales sont supposées infiniment rigides. Les caractéristiques du R22 sont
les suivantes :
• Poids, W 608 daN
• Longueur des pales, R 3.835 m
• Moment de battement de la pale, Ib 63.73 kg.m2
• Corde du profil, c 0.183 m
• Nombre de pale, b 2
• Vrillage des pales, θT W -7 degrés
• Pente du coefficient de portance, a 5.7
• Coefficient de traı̂née du fuselage, SCD 0.49
• Vitesse de rotation du rotor, ψ̇ 520 rpm
Le coefficient de traı̂née du fuselage a été calculé de manière à ajuster la
puissance de traı̂née du fuselage à celle donnée par le diagramme du manuel de

53
vol.
Le programme a consisté à calculer en fonction de la vitesse : l’incidence du
disque rotor, le pas au pied de pale et le coefficient de traı̂née de profil de la pale
pour que la portance égale le poids, que la traı̂née du fuselage soit compensée
par la force motrice du rotor et enfin que la puissance de profil soit celle donnée
par le diagramme du manuel de vol. A partir de ces données, on peut calculer
la puissance de profil, la puissance induite, la puissance de traı̂née du fuselage
et la puissance totale. Pour prendre en compte la puissance débité par le rotor
de queue, 8 % de la puissance que fournit le moteur au rotor a été ajouté à la
puissance induite.
La courbe 3.8 donne le coefficient de traı̂née de profil des pales en fonction
de la vitesse.

Figure 3.8 – Coefficients de traı̂née de profil des pales

La courbe de la figure 3.9 donne le diagramme de puissance avec les points

issus de celui du manuel de vol. Les résultats sont tés proches.

54
 Figure 3.9 – Diagramme de puissance

Enfin les deux courbes suivantes fournissent le pas au pieds de pale et l’in-
cidence du disque rotor toujours en fonction de la vitesse.

55
 Figure 3.10 – Pas au pied de pale

Figure 3.11 – Incidence du rotor

3.4 Etude des coefficients de portance et de trainée

Considérons maintenant un rotor aux caractéristiques suivantes : b = 2,
R = 4 m, c = 0.2 m, θT W = 0, θ0 = 3 deg, Ib = 64 kg.m2 doté d’un profil
tel que a = 5.7 et δ = 0.011. On ne considère pas de déformation de torsion.
Traçons la polaire rotorique de ce rotor en utilisant le modèle de vitesse induite

56
de Shaydakov (voir le programme [Link]). On obtient la
courbe de la figure 3.12. Le coefficient de portance maximum vaut 1 pour une
traı̂née équivalente.

Figure 3.12 – Polaire rotorique

Les deux figures suivantes comparent les valeurs des coefficients CD et CL

calculé avec le modéle de Froude (courbe bleue) et Shaydakov (courbe rouge).

Figure 3.13 – Coefficient de portance

57
 On remarque que le modèle de Froude sous estime la valeur des deux coef-
ficients dès 30 degrés avec une erreur très importante du coefficient de traı̂née
pour les forts angles d’incidence. Le modèle de Froude ne peut être utilisé pour
calculer la vitesse de descente verticale d’un autogire où l’incidence est de 90
degrés.
La valeur maximum du CL voisine de 1 est obtenue vers 45 degrés d’inci-
dence.

Figure 3.14 – Coefficient de trainée

Enfin la courbe 3.15 représente l’évolution du paramètre η en fonction de

l’incidence. On remarque que le rotor passe du mode moulinet frein au mode
turbulent (η=-2) autour de 35 degrés d’incidence

58
 Figure 3.15 – Etat de fonctionnement du rotor

3.5 Evolution des forces et moments en fonction de

l’azimut
La courbe de la figure 3.16 montre l’évolution de la force robotique axiale
d’une seule pale en fonction de l’azimut (voir le programme BladeForcesOve-
[Link]).
On remarque que la force est minimale au niveau de la pale reculante (ψ =
270 degrés) et que cette force est supérieure à la moyenne dans toute la région
avançante et inférieure à la moyenne dans toute la région reculante. Toutefois
à cause des variations à deux fois la fréquence de rotation, la résultante n’est
pas maximale en position avançante (ψ = 90 degrés) mais à un azimut proche
de 180 degrés. Les variations autour de la moyenne sont d’environ +- 10 %.

59
 Figure 3.16 – Evolution de la force rotorique en fonction de l’azimut

La courbe 3.17 représente la force dans le plan du rotor des deux pales.
L’axe x est la force latérale Y et l’axe y et la force longitudinale H. La figure
de cette courbe est parcourue deux fois par tour. Ces forces sont à l’origine des
vibrations ressenties dans le manche.

Figure 3.17 – Evolution de la force dans le plan du rotor

Enfin la courbe de la figure suivante montre l’évolution du couple rotorique

d’une seule pale en fonction de l’azimut (voir le programme BladeTorqueOve-

60
[Link]).

Figure 3.18 – Evolution du couple rotorique en fonction de l’azimut

La variation est quasiment sinusoı̈dale. Le couple maximum est obtenu à

l’azimut 270 degrés (pale reculante) et minimum à l’azimut 90 degrés (pale
avançante).

61
Bibliographie

[1] H. Glauert : A general Theory of the Autogiro. Reports and Memoranda

111. November 1926.
[2] John B. Wheatley : An Aerodynamic Analysis of the Autogiro Rotor with a
Comparison between Calculated and Experimental Results. Report NACA
No. 487, January 1934.
[3] John B. Wheatley : An Analytical and Experimental Study of the Effect
of Periodic Blade Twist on the Thrust, Torque, and Flapping Motion of an
Autogiro Rotor. Report NACA No. 591, January 1937.
[4] John B. Wheatley : An Analysis of the Factors that Determine the Periodic
Twist of an Autogyro Rotor Blade, with a Comparison of Predicted and
Measured Resuls. Report NACA No. 600, April 1937.
[5] V. I. Shaydakov : Aerodynamic calculation for helicopter lifting rotors in
steep descent. Army Foriegn Science and Technology Center. October 1971.
[6] Bramwell’s Helicopter Dynamic
[7] Engauge Digitizer : [Link]

62
Annexe A

Résultats Maple

63
Combined Fourier series
restart ;
Fourier series over pi .. 2*pi
dy := dA0+dB1*sin(psi)-dA0*cos(2*psi)+dB2*sin(2*psi);
(1)
Fourier coefficient of combined series
C0 := expand(1/2/Pi*int(dy,psi=Pi..2*Pi));
(2)

C1 := expand(1/Pi*int(dy*cos(psi),psi=Pi..2*Pi));
(3)

D1 := expand(1/Pi*int(dy*sin(psi),psi=Pi..2*Pi));
(4)

C2 := expand(1/Pi*int(dy*cos(2*psi),psi=Pi..2*Pi));
(5)

D2 := expand(1/Pi*int(dy*sin(2*psi),psi=Pi..2*Pi));
(6)

Forces and moments computation

restart ;
with(CodeGeneration); with(linalg);

(7)

Elementary forces and moments

Flapping moment [Eq ]

(8)
 dM := 1/2*rho*a*c*R^2*x*(theta[R]*UT^2+UP*UT) ;
(8)

Thrust [eq ]
d2L := 1/2*rho*a*c*R*(theta[R]*UT^2+UP*UT) ;
(9)

Rear profile force [Eq ]

d2Hp := 1/2*rho*delta*c*R*UT^2*sin(psi) ;
(10)

Rear induced force [Eq ]

d2Hi := -1/2*rho*a*c*R*((theta[R]*UT^2+UP*UT)*beta(psi)*cos(psi)
+(theta[R]*UP*UT+UP^2)*sin(psi)) ;
(11)

Lateral profile force [Eq ]

d2Yp := -1/2*rho*delta*c*R*UT^2*cos(psi) ;
(12)

Lateral induced force [Eq ]

d2Yi := 1/2*rho*a*c*R*(-(theta[R]*UT^2+UP*UT)*beta(psi)*sin(psi)
+(theta[R]*UP*UT+UP^2)*cos(psi)) ;
(13)

Rotor profile torque [Eq ]

d2Qp := -R*x*1/2*rho*delta*c*R*UT^2 ;
(14)

Rotor induced torque [Eq ]

d2Qi := 1/2*R^2*x*rho*a*c*(theta[R]*UP*UT+UP^2) ;
(15)

Aerodynamic torsional moment [Eq ]

d2N := 1/2*rho*R*c^2*(Cm*UT^2+a*l/c*(theta[R]*UT^2+UP*UT));

(16)

Propeller moment [Eq ]

d2P := -Omega^2*Icb*theta[R];
(17)
Velocity component UP [eq ]
UP := R*Omega*(-mu*beta(psi)*cos(psi)+lambda-x*diff(beta(psi),

(18)
 psi)+x*q*cos(psi)+x*p*sin(psi));
(18)

Velocity component UT [eq ]

UT := R*Omega*(mu*sin(psi)+x) ;
(19)
Pitch deformation [Eq ]
nu := u01*x+u02*x^2+u03*x^3+u04*x^4+u05*x^5+(u11*x+u12*x^2+u13*
x^3+u14*x^4+u15*x^5)*cos(psi)+(v11*x+v12*x^2+v13*x^3+v14*x^4+
v15*x^5)*sin(psi)+(u21*x+u22*x^2+u23*x^3+u24*x^4+u25*x^5)*cos(2*
psi)+(v21*x+v22*x^2+v23*x^3+v24*x^4+v25*x^5)*sin(2*psi);
(20)

Total blade pitch [Eq ]

theta[R] := theta[0]+x*theta[TW]+nu;
(21)

Integration along the blade from x to B or x to 1

dMR := factor(int(dM,x=0..B)):
dMRInv := factor(int(dM,x=0..-mu*sin(psi))):
dL := factor(int(d2L,x=0..B)):
dLInv := factor(int(d2L,x=0..-mu*sin(psi))):
dHp := factor(int(d2Hp,x=0..1)):
dHi := factor(int(d2Hi,x=0..B)):
dYp := factor(int(d2Yp,x=0..1)):
dYi := factor(int(d2Yi,x=0..B)):
dQp := factor(int(d2Qp,x=0..1)):
dQpInv := factor(int(d2Qp,x=0..-mu*sin(psi))):
dQi := factor(int(d2Qi,x=0..B)):
dQiInv := factor(int(d2Qi,x=0..-mu*sin(psi))):
dN := factor(int(d2N,x=x..B)):
dP := factor(int(d2P,x=x..1)):
Flapping equation [eq ]
beta(psi) := a0-a1*cos(psi)-b1*sin(psi)-a2*cos(2*psi)-b2*sin(2*

(22)
 psi) ;
(22)
Force normalisation coefficient [Eq 14]
Kf := rho*b*c*R^3*Omega^2 ;
(23)
Torque normalisation coefficient [Eq 16]
Kq := rho*b*c*R^4*Omega^2 ;
(24)
Evaluation of blade velocity components
UP := factor(expand(UP));
(25)

UT := factor(expand(UT));
(26)
Code generation : UT,UP
Java([ut=subs(sin(psi)=sinPsi,cos(psi)=cosPsi,UT),up=subs(sin
(psi)=sinPsi,cos(psi)=cosPsi,UP)],coercetypes=false,optimize=
true,output="GyroRotorFiles/[Link]");
Angle of attack distribution over rotor disk
Case of no blade twist
eq := aoa = theta[0]+simplify(UP/UT) ;
(27)

X := solve(eq,x) ;
(28)

Code generation : aoa in case of no blade twist

Java([x=subs(theta[0]=theta0,theta[TW]=thetaTW,cos(psi)=cosPsi,
sin(psi)=sinPsi,X)],coercetypes=false,output=
"GyroRotorFiles/[Link]");
Case of linear twist (theta[1]=theta[TW]+u00+u10+v10+u20+v20 ; only first order of pitch
deformation is taking into account)
eq := simplify((aoa*UT-(theta[0]+x*theta[1])*UT - UP)/R/Omega);
(29)

des := discrim(eq,x) ;
(30)

r1 := coeff(eq,x^2) ;
(31)

r2 := coeff(eq,x) ;
(32)

X1 := (-r2 + sqrt(delta))/(2*r1);
(33)
 (33)

X2 := (-r2 - sqrt(delta))/(2*r1);
(34)

Code generation : aoa with linear twist

Java([delta=subs(theta[0]=theta0,theta[1]=theta1,cos(psi)=
cosPsi,sin(psi)=sinPsi,des),x1=subs(theta[0]=theta0,theta[1]=
theta1,cos(psi)=cosPsi,sin(psi)=sinPsi,X1),x2=subs(theta[0]=
theta0,theta[1]=theta1,cos(psi)=cosPsi,sin(psi)=sinPsi,X2)],
coercetypes=false,optimize=true,output="GyroRotorFiles/AOA2.
java");
Thrust
Mean value over one rev, b blades
T := simplify(b/(2*Pi)*(int(dL,psi=0..2*Pi)-2*int(dLInv,psi=Pi.
.2*Pi))) ;
(35)

Define truncate function of power of mu

truncate := proc(input,order)
local res;
res := subs(a1=a1*mu,b1=b1*mu,a2=a2*mu^2,b2=b2*mu^2,simplify
(input)):
res := subs(u10=u10*mu,u11=u11*mu,u12=u12*mu,u13=u13*mu,u14=
u14*mu,res):
res := subs(v10=v10*mu,v11=v11*mu,v12=v12*mu,v13=v13*mu,v14=
v14*mu,res):
res := subs(u20=u20*mu^2,u21=u21*mu^2,u22=u22*mu^2,u23=u23*
mu^2,u24=u24*mu^2,res):
res := subs(v20=v20*mu^2,v21=v21*mu^2,v22=v22*mu^2,v23=v23*
mu^2,v24=v24*mu^2,res):
res := expand(convert(taylor(res,mu,order),polynom)):
res := subs(a1=a1/mu,b1=b1/mu,a2=a2/mu^2,b2=b2/mu^2,res):
res := subs(u10=u10/mu,u11=u11/mu,u12=u12/mu,u13=u13/mu,u14=
u14/mu,res):
res := subs(v10=v10/mu,v11=v11/mu,v12=v12/mu,v13=v13/mu,v14=
v14/mu,res):
res := subs(u20=u20/mu^2,u21=u21/mu^2,u22=u22/mu^2,u23=
u23/mu^2,u24=u24/mu^2,res):
res := subs(v20=v20/mu^2,v21=v21/mu^2,v22=v22/mu^2,v23=
v23/mu^2,v24=v24/mu^2,res):
res ;
end proc;
(36)
Thrust coefficient
CT := truncate(T/Kf,5); T := Kf*CT:
(37)

Functions to get terms ordered

ex1 := proc(r)
global input;
local res;
res:=expand(coeff(input,r))*r;
input:=simplify(input-res);
res;
end proc;
(38)

ex2 := proc(r,rank)
local i,res;
global input;
res := 0;
for i from 1 to rank do
res := res + expand(coeff(input,r||i))*r||i;
od;
input:=simplify(input-res);

(39)
 res;
end proc;
(39)

CT Terms
input := 2*CT/a:
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);

0 (40)
ex1(p);ex1(q);ex1(lambda);ex1(theta[0]);ex1(theta[TW]);ex1(a0);
ex1(a1);ex1(b1);ex1(a2);ex1(b2);

0
0
(41)

input;
(42)
 0 (42)
Periodic component of thrust
TP := b*combine(dL/Kf,trig):
TP1S := coeff(TP,sin(psi)):
TP1C := coeff(TP,cos(psi)):
TP2S := coeff(TP,sin(2*psi)):
TP2C := coeff(TP,cos(2*psi)):
Rear force
Profile mean value over one rev, b blades
Hp := b/(2*Pi)*int(dHp,psi=0..2*Pi) ;
(43)

Profile coefficient
CHp := simplify(Hp/Kf);
(44)

Induced mean value over one rev, b blades

Hi := simplify(b/(2*Pi)*int(dHi,psi=0..2*Pi)) :
Induced coefficient
CHi := truncate(Hi/Kf,5): Hi := Kf*CHi:
CHi Terms
input := 2*CHi/a:
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);
 (45)

ex1(p);ex1(q);ex1(lambda);ex1(theta[0]);ex1(theta[TW]);ex1(a0);
ex1(a1);ex1(b1);ex1(a0^2);ex1(a1^2);ex1(a2);ex1(b2);

0
0 (46)
input;
0 (47)
Periodic component of rear profile coefficient
HpP := b*combine(dHp/Kf,trig):
HpP1S := coeff(HpP,sin(psi)):
HpP2C := coeff(HpP,cos(2*psi)):
Periodic component of rear induced coefficient
 HiP := b*combine(dHi/Kf,trig) :
HiP1S := coeff(HiP,sin(psi)):
HiP1C := coeff(HiP,cos(psi)):
HiP2S := coeff(HiP,sin(2*psi)):
HiP2C := coeff(HiP,cos(2*psi)):
Lateral force
Profile mean value over one rev, b blades
Yp := b/(2*Pi)*int(dYp,psi=0..2*Pi) ;
(48)
Induced mean value over one rev, b blades
Yi := simplify(b/(2*Pi)*int(dYi,psi=0..2*Pi)) :
Induced coefficient
CYi := truncate(Yi/Kf,5): Yi := Kf*CYi:
CYi Terms
input := 2*CYi/a:
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);
 (49)

ex1(p);ex1(q);ex1(lambda);ex1(theta[0]);ex1(theta[TW]);ex1(a0);
ex1(a1);ex1(b1);ex1(a2);ex1(b2);

0
0 (50)
input;
0 (51)
Periodic component of lateral profile coefficient
YpP := b*combine(dYp/Kf,trig) :
YpP1C := coeff(YpP,cos(psi)):
YpP2S := coeff(YpP,sin(2*psi)):
Periodic component of lateral induced coefficient
YiP := b*combine(dYi/Kf,trig) :
YiP1S := coeff(YiP,sin(psi)):
 YiP1C := coeff(YiP,cos(psi)):
YiP2S := coeff(YiP,sin(2*psi)):
YiP2C := coeff(YiP,cos(2*psi)):
Rotor torque
Profile mean value over one rev, b blades
Qp := simplify(b/(2*Pi)*(int(dQp,psi=0..2*Pi)-2*int(dQpInv,psi=
Pi..2*Pi)));
(52)

Profile coefficient
CQp := simplify(Qp/Kq);
(53)

Induced mean value over one rev, b blades

Qi := simplify(b/(2*Pi)*(int(dQi,psi=0..2*Pi)-2*int(dQiInv,psi=
Pi..2*Pi))):
Induced torque coefficient. Force dependency of mu
CQi := truncate(Qi/Kq,5): Qi := Kq*CQi:
CQi terms
input := 2*CQi/a:
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);
 (54)

ex1(p^2);ex1(p);ex1(q^2);ex1(q);ex1(lambda^2);ex1(lambda);ex1
(a2^2);ex1(a2);ex1(a0^2);ex1(a0);ex1(a1^2);ex1(b1^2);ex1(b2^2);
ex1(b2);
 (55)

input;
0 (56)
Periodic component of rotor torque profile coefficient
QpP := b*combine(dQp/Kq,trig) :
QpP1S := coeff(QpP,sin(psi)):
QpP2C := coeff(QpP,cos(2*psi)):
Periodic component of rotor torque induced coefficient
QiP := b*combine(dQi/Kq,trig) :
QiP1S := coeff(QiP,sin(psi)):
QiP1C := coeff(QiP,cos(psi)):
QiP2S := coeff(QiP,sin(2*psi)):
QiP2C := coeff(QiP,cos(2*psi)):
Aerodynamic torsional moment
N := simplify(1/(2*Pi)*int(dN,psi=0..2*Pi)):
Propeller moment
P := simplify(1/(2*Pi)*int(dP,psi=0..2*Pi)):
Code generation : CQp, CQi
Java([cqp=subs(theta[0]=theta0,theta[TW]=thetaTW,CQp),cqi=subs
(theta[0]=theta0,theta[TW]=thetaTW,CQi)],coercetypes=false,
optimize=true,output="GyroRotorFiles/[Link]");
Code generation : CT
Java([ct=subs(theta[0]=theta0,theta[TW]=thetaTW,CT)],
coercetypes=false,optimize=true,output="GyroRotorFiles/[Link]")
;
Code generation : CHp, CHi, CYi
Java([chp=CHp,chi=subs(theta[0]=theta0,theta[TW]=thetaTW,CHi),
cyi=subs(theta[0]=theta0,theta[TW]=thetaTW,CYi)],coercetypes=
false,optimize=true,output="GyroRotorFiles/[Link]");
Code generation : N
Java([nx=subs(theta[0]=theta0,theta[TW]=thetaTW,dN)],
coercetypes=false,optimize=true,output="GyroRotorFiles/[Link]");
Code generation : P (Propeller moment)
Java([px=subs(theta[0]=theta0,theta[TW]=thetaTW,dP)],
coercetypes=false,optimize=true,output="GyroRotorFiles/[Link]");
Code generation : TP1S, TP1C, TP2S, TP2C
Java([tp1S=subs(theta[0]=theta0,theta[TW]=thetaTW,TP1S),tp1C=
subs(theta[0]=theta0,theta[TW]=thetaTW,TP1C),tp2S=subs(theta[0]=
theta0,theta[TW]=thetaTW,TP2S),tp2C=subs(theta[0]=theta0,theta
[TW]=thetaTW,TP2C)],coercetypes=false,optimize=true,output=
"GyroRotorFiles/[Link]");
Code generation : HpP1S, HpP2C
Java([hpP1S=subs(theta[0]=theta0,theta[TW]=thetaTW,HpP1S),hpP2C=
subs(theta[0]=theta0,theta[TW]=thetaTW,HpP2C)],coercetypes=
false,optimize=true,output="GyroRotorFiles/[Link]");
Code generation : HiP1S, HiP1C, HiP2S, HiP2C
Java([hiP1S=subs(theta[0]=theta0,theta[TW]=thetaTW,HiP1S),hiP1C=
subs(theta[0]=theta0,theta[TW]=thetaTW,HiP1C),hiP2S=subs(theta
[0]=theta0,theta[TW]=thetaTW,HiP2S),hiP2C=subs(theta[0]=theta0,
theta[TW]=thetaTW,HiP2C)],coercetypes=false,optimize=true,
output="GyroRotorFiles/[Link]");
Code generation : YpP1C, YpP2S
Java([ypP1C=subs(theta[0]=theta0,theta[TW]=thetaTW,YpP1C),ypP2S=
subs(theta[0]=theta0,theta[TW]=thetaTW,YpP2S)],coercetypes=
false,optimize=true,output="GyroRotorFiles/[Link]");
Code generation : YiP1S, YiP1C, YiP2S, YiP2C
Java([yiP1S=subs(theta[0]=theta0,theta[TW]=thetaTW,YiP1S),yiP1C=
subs(theta[0]=theta0,theta[TW]=thetaTW,YiP1C),yiP2S=subs(theta
[0]=theta0,theta[TW]=thetaTW,YiP2S),yiP2C=subs(theta[0]=theta0,
theta[TW]=thetaTW,YiP2C)],coercetypes=false,optimize=true,
output="GyroRotorFiles/[Link]");
Code generation : QpP1S, QpP2C
Java([qpP1S=subs(theta[0]=theta0,theta[TW]=thetaTW,QpP1S),qpP2C=
subs(theta[0]=theta0,theta[TW]=thetaTW,QpP2C)],coercetypes=
false,optimize=true,output="GyroRotorFiles/[Link]");
Code generation : QiP1S, QiP1C, QiP2S, QiP2C
Java([qiP1S=subs(theta[0]=theta0,theta[TW]=thetaTW,QiP1S),qiP1C=
subs(theta[0]=theta0,theta[TW]=thetaTW,QiP1C),qiP2S=subs(theta
[0]=theta0,theta[TW]=thetaTW,QiP2S),qiP2C=subs(theta[0]=theta0,
theta[TW]=thetaTW,QiP2C)],coercetypes=false,optimize=true,
output="GyroRotorFiles/[Link]");
Compute torsional deformation coefficient
Total integral torsionnal moment
NP := int(dN+dP,x=0..x):
NP := combine(NP,trig):
Deformation equation
DP := K*nu-NP:
DP := subs(cos(psi)=cpsi,sin(psi)=spsi,cos(2*psi)=c2psi,sin(2*
 psi)=s2psi,cos(3*psi)=0,sin(3*psi)=0,cos(4*psi)=0,sin(4*psi)=0,
DP):
DP := subs(x^6=0,x^7=0,x^8=0,x^9=0,DP):
Define matrix to solve linear system
m := Matrix(25,25);

(57)

n := Vector(25);

(58)

Torsional inertia force for plotting

S := Omega^2*Icb*((u11*x+u12*x^2+u13*x^3+u14*x^4+u15*x^5)*cos
(psi)+(v11*x+v12*x^2+v13*x^3+v14*x^4+v15*x^5)*sin(psi)+4*(u21*x+
u22*x^2+u23*x^3+u24*x^4+u25*x^5)*cos(2*psi)+4*(v21*x+v22*x^2+
v23*x^3+v24*x^4+v25*x^5)*sin(2*psi));
(59)

Flaping equation
unassign('beta(psi)');
A := diff(beta(psi),[psi$2])+beta(psi) + 2*q*sin(psi)- 2*p*cos
(psi) ;
(60)

Two orders Fourier series [Eq ]

beta(psi) := a0-a1*cos(psi)-b1*sin(psi)-a2*cos(2*psi)-b2*sin(2*
psi) ;
(61)
AP := subs(cos(psi)=cpsi,sin(psi)=spsi,cos(2*psi)=c2psi,sin(2*
psi)=s2psi,A);
(62)
Fourrier coefficient of nominal flow
dMR1 := combine(dMR,trig):
 dMR1 := subs(sin(3*psi)=0,cos(3*psi)=0,sin(4*psi)=0,cos(4*psi)=
0,sin(5*psi)=0,cos(5*psi)=0,sin(6*psi)=0,cos(6*psi)=0,sin(7*psi)
=0,cos(7*psi)=0,sin(8*psi)=0,cos(8*psi)=0,sin(9*psi)=0,cos(9*
psi)=0,sin(10*psi)=0,cos(10*psi)=0,sin(psi)=spsi,cos(psi)=cpsi,
sin(2*psi)=s2psi,cos(2*psi)=c2psi,dMR1):
A1 := coeff(dMR1,cpsi): B1 := coeff(dMR1,spsi) : A2 := coeff
(dMR1,c2psi) : B2 := coeff(dMR1,s2psi) :
A0 := simplify(dMR1-A1*cpsi-B1*spsi-A2*c2psi-B2*s2psi):
Fourrier coefficient of reverse flow
dMR2 := combine(-2*dMRInv,trig):
dMR2 := subs(sin(3*psi)=0,cos(3*psi)=0,sin(4*psi)=0,cos(4*psi)=
0,sin(5*psi)=0,cos(5*psi)=0,sin(6*psi)=0,cos(6*psi)=0,sin(7*psi)
=0,cos(7*psi)=0,sin(8*psi)=0,cos(8*psi)=0,sin(9*psi)=0,cos(9*
psi)=0,sin(10*psi)=0,cos(10*psi)=0,sin(11*psi)=0,cos(11*psi)=0,
sin(psi)=spsi,cos(psi)=cpsi,sin(2*psi)=s2psi,cos(2*psi)=c2psi,
dMR2):
dA1 := coeff(dMR2,cpsi) : dB1 := coeff(dMR2,spsi) : dA2 := coeff
(dMR2,c2psi) : dB2 := coeff(dMR2,s2psi) :
dA0 := simplify(dMR2-dA1*cpsi-dB1*spsi-dA2*c2psi-dB2*s2psi):
Combined fourier coefficient
C0 := A0-dB1/Pi+(1/2)*dA0:
C1 := A1-(4/3)*dB2/Pi:
D1 := B1+(1/2)*dB1-(8/3)*dA0/Pi:
C2 := A2-(1/2)*dA0+(2/3)*dB1/Pi:
D2 := B2+(1/2)*dB2:
dMRc := (C0+C1*cpsi+D1*spsi+C2*c2psi+D2*s2psi)/Ib/Omega^2 :
Gamma is the Lock number [Eq ]
Ib := rho*a*c*R^4/Gamma ;

(63)

dMRc := truncate(dMRc,5):
Code generation : dMRc
Java([Mx=subs(theta[0]=theta0,theta[TW]=thetaTW,dMRc)],
coercetypes=false,optimize=true,output="GyroRotorFiles/[Link]");
Flapping moment coefficient
input := simplify(coeff(2*dMRc/Gamma,cpsi)):
residual := input*cpsi:
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);
0

(64)
 0
0
(64)

ex1(p);ex1(q);ex1(lambda);ex1(a2);ex1(b2);ex1(theta[0]);ex1
(theta[TW]);ex1(a0);ex1(a1);ex1(b1);
0

0
0
0

0
(65)

input := simplify(coeff(2*dMRc/Gamma,spsi)):
residual := residual + input*spsi:
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);

0 (66)
ex1(p);ex1(q);ex1(lambda);ex1(a2);ex1(b2);ex1(theta[0]);ex1
(theta[TW]);ex1(a0);ex1(a1);ex1(b1);

0
 0

0 (67)
input := simplify(coeff(2*dMRc/Gamma,c2psi)):
residual := residual + input*c2psi:
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);

0 (68)
ex1(p);ex1(q);ex1(lambda);ex1(a2);ex1(b2);ex1(theta[0]);ex1
(theta[TW]);ex1(a0);ex1(a1);ex1(b1);

(69)
 0 (69)
input := simplify(coeff(2*dMRc/Gamma,s2psi)):
residual := residual + input*s2psi:
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);
0

0
0
(70)

ex1(p);ex1(q);ex1(lambda);ex1(a2);ex1(b2);ex1(theta[0]);ex1
(theta[TW]);ex1(a0);ex1(a1);ex1(b1);
0

0
0
0

0
(71)

input := simplify(simplify(2*dMRc/Gamma) - residual):

ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);

0 (72)
ex1(p);ex1(q);ex1(lambda);ex1(a2);ex1(b2);ex1(theta[0]);ex1
 (theta[TW]);ex1(a0);ex1(a1);ex1(b1);

0
0
0 (73)
Final flapping equation
AP := AP - dMRc:
Sine and cosine coefficient
Aspsi := coeff(AP,spsi): Acpsi := coeff(AP,cpsi): As2psi :=
coeff(AP,s2psi): Ac2psi := coeff(AP,c2psi):
Constant term
Acst := simplify(AP-Aspsi*spsi-Acpsi*cpsi-As2psi*s2psi-Ac2psi*
c2psi):
Flapping coefficient equations
eq1 := Acst = 0 : eq2 := Aspsi = 0 : eq3 := Acpsi = 0 : eq4 :=
As2psi = 0 : eq5 := Ac2psi = 0:
Solution
solve({eq1,eq2,eq3},{a0,a1,b1}) :
assign(%) ;
Code generation : a0, a1, b1
Java([ga0=subs(theta[0]=theta0,theta[TW]=thetaTW,a0),ga1=subs
(theta[0]=theta0,theta[TW]=thetaTW,a1),gb1=subs(theta[0]=theta0,
theta[TW]=thetaTW,b1)],coercetypes=false,optimize=true,output=
"GyroRotorFiles/[Link]");
Conning angle
expand(a0) ;
(74)
 (74)

input := simplify(2*a0/Gamma):
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);

0 (75)
ex1(p) ; ex1(q) ; ex1(a2) ; ex1(b2) ; ex1(lambda) ; ex1(theta[0]
); ex1(theta[TW]);

0
0

(76)

Longitudinal flapping
expand(a1);

(77)
 (77)

input := simplify(a1/(2*mu/(B^4-1/2*mu^2*B^2))):
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);

0 (78)
ex1(p) ; ex1(q) ; ex1(a2) ; ex1(b2) ; ex1(lambda) ; ex1(theta[0]
); ex1(theta[TW]);

(79)
 (79)
Lateral flapping
expand(b1):
b1r := simplify(b1-4*mu/(B^4+1/2*mu^2*B^2)*a0*(1/3*B^3+1/9*
mu^3/Pi)) ;
(80)

b1bis := 4*mu/(B^4+1/2*mu^2*B^2)*a0bis*(1/3*B^3+1/9*mu^3/Pi)+
b1r;

(81)

input := simplify(b1bis/(4*mu/(B^4+1/2*mu^2*B^2))):
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);
0

0
0
(82)

ex1(p) ; ex1(q) ; ex1(a0bis) ; ex1(a2) ; ex1(b2) ; ex1(lambda) ;

ex1(theta[0]); ex1(theta[TW]);
 0
0
0
0 (83)
Second order flapping coefficient solution
solve({eq4,eq5},{a2,b2}) :
assign(%) ;
Longitudinal flapping Taylor series
a2dev := series(a2,mu=0):
Constant term (u2, v2)
a20 := simplify(convert(series(a2,mu=0,1),polynom));
(84)

mu factor (u1, v1)

a21 := mu*simplify(subs(lambda=0,theta[0]=0,theta[TW]=0,coeff
(a2dev,mu)));
(85)

mu^2 factor (lambda, theta[0], theta[TW], u0)

a22 := mu^2*simplify(subs(u11=0,u12=0,u13=0,u14=0,u15=0,v11=0,
v12=0,v13=0,v14=0,v15=0,u21=0,u22=0,u23=0,u24=0,u25=0,v21=0,v22=
0,v23=0,v24=0,v25=0,coeff(a2dev,mu^2)));
(86)

a2:= simplify(a20+a21+a22);
(87)
 (87)

a2 terms
input := simplify(a2/(mu^2*Gamma/(Gamma^2*B^8+144))):
ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);

(88)

ex1(p) ; ex1(q) ; ex1(lambda) ; ex1(theta[0]); ex1(theta[TW]);

(89)
 (89)

Lateral flapping Taylor series

b2dev := series(b2,mu=0):
Constant term (u2, v2)
b20 := simplify(convert(series(b2,mu=0,1),polynom));
(90)

mu factor (u1, v1)

b21 := mu*simplify(subs(lambda=0,theta[0]=0,theta[TW]=0,coeff
(b2dev,mu)));
(91)

mu^2 factor (lambda, theta[0], theta[TW], u0)

b22 := mu^2*simplify(subs(u11=0,u12=0,u13=0,u14=0,u15=0,v11=0,
v12=0,v13=0,v14=0,v15=0,u21=0,u22=0,u23=0,u24=0,u25=0,v21=0,v22=
0,v23=0,v24=0,v25=0,coeff(b2dev,mu^2)));
(92)

b2 := simplify(b20+b21+b22);
(93)

b2 terms
input := simplify(b2/(-mu^2*Gamma^2/(Gamma^2*B^8+144))):
 ex2(u0,5);ex2(u1,5);ex2(v1,5);ex2(u2,5);ex2(v2,5);

(94)

ex1(p) ; ex1(q) ; ex1(lambda) ; ex1(theta[0]); ex1(theta[TW]);

(95)

Code generation : a2, b2

Java([ga2=subs(theta[0]=theta0,theta[TW]=thetaTW,a2),gb2=subs
(theta[0]=theta0,theta[TW]=thetaTW,b2)],coercetypes=false,
optimize=true,output="GyroRotorFiles/[Link]");
Compute torsional deformation coefficient
Define procedure to compute matrix coefficients
computeMatrix := proc(input,order)
local j,l,coef;
global k,m,n;
for j from 1 to order do;
coef := coeff(input,x^j);
m[k,1]:=coeff(coef,u01):m[k,2]:=coeff(coef,u02):m[k,3]:=coeff
(coef,u03):m[k,4]:=coeff(coef,u04):m[k,5]:=coeff(coef,u05):
m[k,6]:=coeff(coef,u11):m[k,7]:=coeff(coef,u12):m[k,8]:=coeff
(coef,u13):m[k,9]:=coeff(coef,u14):m[k,10]:=coeff(coef,u15):
m[k,11]:=coeff(coef,v11):m[k,12]:=coeff(coef,v12):m[k,13]:=
coeff(coef,v13):m[k,14]:=coeff(coef,v14):m[k,15]:=coeff(coef,
v15):
 m[k,16]:=coeff(coef,u21):m[k,17]:=coeff(coef,u22):m[k,18]:=
coeff(coef,u23):m[k,19]:=coeff(coef,u24):m[k,20]:=coeff(coef,
u25):
m[k,21]:=coeff(coef,v21):m[k,22]:=coeff(coef,v22):m[k,23]:=
coeff(coef,v23):m[k,24]:=coeff(coef,v24):m[k,25]:=coeff(coef,
v25):
n[k] := m[k,1]*u01+m[k,2]*u02+m[k,3]*u03+m[k,4]*u04+m[k,5]*
u05+m[k,6]*u11+m[k,7]*u12+m[k,8]*u13+m[k,9]*u14+m[k,10]*u15+m[k,
11]*v11+m[k,12]*v12+m[k,13]*v13+m[k,14]*v14+m[k,15]*v15:
n[k] := n[k]+m[k,16]*u21+m[k,17]*u22+m[k,18]*u23+m[k,19]*u24+m
[k,20]*u25+m[k,21]*v21+m[k,22]*v22+m[k,23]*v23+m[k,24]*v24+m[k,
25]*v25:
n[k] := simplify(n[k]-coef):
k:=k+1;
od;
end proc;
(96)
Compute matrix m and n
k:=1; K1 := coeff(DP,cpsi): K2 := coeff(DP,spsi): K3 := coeff
(DP,c2psi): K4 := coeff(DP,s2psi):
(97)
computeMatrix(K1,5):
computeMatrix(K2,5):
computeMatrix(K3,5):
computeMatrix(K4,5):
K5 := simplify(DP-K1*cpsi-K2*spsi-K3*c2psi-K4*s2psi):
computeMatrix(K5,5):
m:= map(simplify,m):
Code génération : m
Java(subs(theta[0]=theta0,theta[TW]=thetaTW,m),coercetypes=
false,optimize=true,output="GyroRotorFiles/[Link]");
Code génération : n
Java(subs(theta[0]=theta0,theta[TW]=thetaTW,n),coercetypes=
false,optimize=true,output="GyroRotorFiles/[Link]");
Extract sine and cosine components of total torsionnal moment for plotting
NP := subs(cos(psi)=cpsi,sin(psi)=spsi,cos(2*psi)=c2psi,sin(2*
psi)=s2psi,sin(3*psi)=0,cos(3*psi)=0,sin(4*psi)=0,cos(4*psi)=0,
NP):
Mu1 := coeff(NP,cpsi):
Mv1 := coeff(NP,spsi):
Mu2 := coeff(NP,c2psi):
Mv2 := coeff(NP,s2psi):
Mu0 := simplify(NP-Mu1*cpsi-Mv1*spsi-Mu2*c2psi-Mv2*s2psi):
Numerical application : validation test

(98)
 mu := 0.1678124092350715 ; lambda := 0.013486848994919694 ;
Omega := 36.96635711289101;

(98)
a := 5.7 ; delta := 0.011; Cm := 0.005 ; xac := 0.278 ; rho :=
1.225; b := 2 ; R := 4; c := 0.2 ; xcg := 0.3 ; theta[TW] :=
evalf[15](2*Pi/180) ; GJ := 6350 ; B := 0.97 ; Ib := 64.0 ; Ic
:= 64.04 ; Gamma := rho*a*c*R^4/Ib ; l:= (xcg-xac)*c ; Icb :=
Ic - Ib ; K := evalf(GJ/R) ; theta[0] := evalf[15](2*Pi/180) ;
q:=evalf[15](2/Omega*Pi/180) ; p:=evalf[15](-3/Omega*Pi/180) ;

(99)
m1:=map(evalf,m):n1:=map(evalf,n):
Pitch deformation coefficient
d := evalm(inverse(m1)&*n1);
(100)
 (100)

u01:=d[1]:u02:=d[2]:u03:=d[3]:u04:=d[4]:u05:=d[5]:u11:=d[6]
:u12:=d[7]:u13:=d[8]:u14:=d[9]:u15:=d[10]:v11:=d[11]:v12:=d[12]
:v13:=d[13]:v14:=d[14]:v15:=d[15]:u21:=d[16]:u22:=d[17]:u23:=d
[18]:u24:=d[19]:u25:=d[20]:v21:=d[21]:v22:=d[22]:v23:=d[23]
:v24:=d[24]:v25:=d[25]:
Pitch deformation at tip
evalf[15]((u01+u02+u03+u04+u05)*180/Pi) ;
0.211135519168619 (101)
evalf[15]((u11+u12+u13+u14+u15)*180/Pi) ;
0.00972186712599437 (102)
evalf[15]((v11+v12+v13+v14+v15)*180/Pi) ;
0.0450938074704633 (103)
evalf[15]((u21+u22+u23+u24+u25)*180/Pi) ;
0.00896996365230856 (104)
evalf[15]((v21+v22+v23+v24+v25)*180/Pi) ;
(105)
Flapping coefficient
evalf[15](a0*180/Pi); evalf[15](a1*180/Pi);evalf[15](b1*180/Pi);
evalf[15](a2*180/Pi); evalf[15](b2*180/Pi);
3.03409110926474
1.76342065635421
0.894100707588128
0.0627955850279689
(106)
Normalization coefficient for one blade
Kf1 := evalf[15](Kf/b) ; Kq1 := evalf[15](Kq/2) ;

(107)
Thrust
evalf[15](T);evalf[15](Kf1*TP1S);evalf[15](Kf1*TP1C);evalf[15]
(Kf1*TP2S);evalf[15](Kf1*TP2C);
3232.77629417085
143.732273638253
58.609531427432

105.851485983171 (108)
Rear force
evalf[15](Hp) ; evalf[15](Hi) ;
19.7763495505108
(109)
 105.385938824211 (109)
evalf[15](Kf1*HpP1S);evalf[15](Kf1*HpP2C);
41.7716899076188
(110)
evalf[15](Kf1*HiP1S);evalf[15](Kf1*HiP1C);evalf[15](Kf1*HiP2S);
evalf[15](Kf1*HiP2C);

6.92035704001788
(111)
Lateral force
evalf[15](Yi) ;
19.3624497239309 (112)
evalf[15](Kf1*YpP1C);evalf[15](Kf1*YpP2S);

(113)
evalf[15](Kf1*YiP1S);evalf[15](Kf1*YiP1C);evalf[15](Kf1*YiP2S);
evalf[15](Kf1*YiP2C);

53.0443573303871

(114)
Torque
evalf[15](Qp) ; evalf[15](Qi) ;

242.309982671682 (115)
evalf[15](Kq1*QpP1S);evalf[15](Kq1*QpP2C);

3.31871686394616 (116)
evalf[15](Kq1*QiP1S);evalf[15](Kq1*QiP1C);evalf[15](Kq1*QiP2S);
evalf[15](Kq1*QiP2C);

25.3555997328135

8.09729727480028 (117)
Plot aerodynamic torsional moment, propeller moment and 10 time inertial forces at given azimuth
plot ({subs(psi=0,dN),subs(psi=0,dP),10*subs(psi=0,S)},x = 0 ..
B) ;
 10

Plot pitch deformation at given azimuth and compare with valur computed with total torsional moment
plot({subs(psi=0,nu)*180/Pi,int(subs(psi=0,dN+dP),x=0..x)*
180/Pi/K},x=0..B);
 0
0
x
Check u1
plot({(u11*x+u12*x^2+u13*x^3+u14*x^4+u15*x^5)*180/Pi,evalf(Mu1*
180/Pi/K)},x=0..B);
 0
0
x
Check v1
plot({(v11*x+v12*x^2+v13*x^3+v14*x^4+v15*x^5)*180/Pi,evalf(Mv1*
180/Pi/K)},x=0..B);
 0
0
x
Check u2
plot({(u21*x+u22*x^2+u23*x^3+u24*x^4+u25*x^5)*180/Pi,evalf(Mu2*
180/Pi/K)},x=0..B);
 0
0
x
Check v2
plot({(v21*x+v22*x^2+v23*x^3+v24*x^4+v25*x^5)*180/Pi,evalf(Mv2*
180/Pi/K)},x=0..B);
 0

Check mean term

plot({(u01*x+u02*x^2+u03*x^3+u04*x^4+u05*x^5)*180/Pi,evalf(Mu0*
180/Pi/K)},x=0..B);
 0
0
x
Harmonic torque
QiP2 := evalf(Qi/Kq+QiP1S*sin(psi)+QiP1C*cos(psi)+QiP2S*sin(2*
psi)+QiP2C*cos(2*psi));
(118)

plot ({QiP,QiP2},psi = 0 .. 2*Pi);

 0
1 2 3 4 5 6

Harmonic thrust
TP2 := evalf(T/Kf+TP1S*sin(psi)+TP1C*cos(psi)+TP2S*sin(2*psi)+
TP2C*cos(2*psi));
(119)

plot ({TP,TP2},psi = 0 .. 2*Pi);

 0 1 2 3 4 5 6

Numerical application : KD1 Gyro

mu := 0.31388398812868695 ; lambda := 0.013511805577481526 ;
Omega := 22.412411639740625;

(120)
a:=5.862 ; delta := 0.013101 ; Cm := -0.056 ; xac := 0.242 ; rho
:= 1.191; b := 3 ; R := 6.096; c := 0.3048 ; xcg := 0.28 ;
theta[0] := 0.096; theta[TW] := 0 ; GJ := 11021.26 ; B := 0.97 ;
Ib := 237.268 ; Ic := 237.268 ; Gamma := rho*a*c*R^4/Ib ; l:=
(xcg-xac)*c ; Icb := Ic - Ib ; K := evalf(GJ/R) ; q:=0 ; p:=0 ;
 (121)
m1:=map(evalf,m):n1:=map(evalf,n):
Pitch deformation coefficient
d := evalm(inverse(m1)&*n1);
(122)

u01:=d[1]:u02:=d[2]:u03:=d[3]:u04:=d[4]:u05:=d[5]:u11:=d[6]
:u12:=d[7]:u13:=d[8]:u14:=d[9]:u15:=d[10]:v11:=d[11]:v12:=d[12]
:v13:=d[13]:v14:=d[14]:v15:=d[15]:u21:=d[16]:u22:=d[17]:u23:=d
[18]:u24:=d[19]:u25:=d[20]:v21:=d[21]:v22:=d[22]:v23:=d[23]
:v24:=d[24]:v25:=d[25]:
Check u1
plot({(u11*x+u12*x^2+u13*x^3+u14*x^4+u15*x^5)*180/Pi,evalf(Mu1*
180/Pi/K)},x=0..B);
 0

Check v1
plot({(v11*x+v12*x^2+v13*x^3+v14*x^4+v15*x^5)*180/Pi,evalf(Mv1*
180/Pi/K)},x=0..B);
 0

Check u2
plot({(u21*x+u22*x^2+u23*x^3+u24*x^4+u25*x^5)*180/Pi,evalf(Mu2*
180/Pi/K)},x=0..B);
 0
0
x
Check v2
plot({(v21*x+v22*x^2+v23*x^3+v24*x^4+v25*x^5)*180/Pi,evalf(Mv2*
180/Pi/K)},x=0..B);
 0

Check mean term

plot({(u01*x+u02*x^2+u03*x^3+u04*x^4+u05*x^5)*180/Pi,evalf(Mu0*
180/Pi/K)},x=0..B);
0