Une méthode d'identification des coefficients aérodynamiques d'une flèche à partir de données prises en vol est proposée. Elle est basée sur l'optimisation d'une fonctionnelle où interviennent les mesures (vitesses, rotations, ), ainsi que les coefficients aérodynamiques. Le vol du projectile est représenté par un modèle à six d.d.l. où interviennent 8 coefficients aérodynamiques. L'optimisation de la fonctionnelle conduit à un système non linéaire qui fournit après résolution les coefficients aérodynamiques du projectile ainsi qu'une trajectoire recalculée s'approchant au mieux des mesures

BAILLY, Patrice

DELVARE, Franck

DEMAILLY H.

GRIGNON, Christophe

HEDDADJ, Settie

I-Revues

20ème Congrès Français de Mécanique                                                                  Besançon, 29 août au 2 septembre 2011 
  1 
Identification inverse des coefficients aérodynamiques d’un 
engin hypersonique à partir de données issues de vol 
H. DEMAILLYab, F. DELVAREb, S. HEDDADJa, C. GRIGNONc, P. BAILLYb 
a. Nexter Munitions, 7 route de Guerry, 18023 Bourges Cédex (France) 
b. Institut PRISME, ENSI de Bourges, 88 Boulevard Lahitolle, 18020 Bourges Cédex (France) 
c. Direction Générale de l’Armement (DGA) Techniques terrestres, Rocade Est – Echangeur de Guerry, 
18000 Bourges (FRANCE) 
Résumé : 
Une méthode d'identification des coefficients aérodynamiques d'un projectile flèche à partir de mesures 
prises en vol est proposée. Le vol du projectile est représenté par un modèle à six degrés de liberté intégrant 
les équations de la mécanique du vol, où interviennent 8 coefficients aérodynamiques. La technique 
développée est basée sur l'optimisation d'une fonctionnelle, au sens des moindres carrés, où interviennent les 
mesures (vitesses, rotations,…), ainsi que les coefficients aérodynamiques. L’intégration des équations du 
mouvement et l'optimisation de la fonctionnelle conduisent à un système non linéaire qui fournit, après 
résolution au moyen d’une technique de Newton-Raphson, les coefficients aérodynamiques du projectile 
ainsi qu'une trajectoire recalculée s'approchant au mieux des mesures. 
Abstract : 
An identification technique for the aerodynamic coefficients of a Kinetic Energy projectile from flight data is 
proposed. The flight of the projectile is represented by a six degrees of freedom model where eight 
aerodynamic coefficients are involved. The identification technique is based on the optimization of a 
functional which involves the measurements (velocities, rotations ...) and the aerodynamic coefficients. The 
numerical integration of the equations of motion and the minimization of the functional lead to a non linear 
system that, after resolution, provides the aerodynamic coefficients of the projectile and a recomputed flight 
which fits the best the measures. 
Mots clefs : Identification de paramètres, optimisation non linéaire, balistique extérieure 
1 Introduction  
La connaissance des coefficients aérodynamiques d’un projectile est indispensable pour comprendre, 
maîtriser et prévoir son attitude en vol, en termes de précision et de stabilité. L’objectif est donc de 
développer une méthode fournissant le jeu de coefficients aérodynamiques ayant généré le vol du projectile à 
partir de mesures prises en vol. Les mesures à notre disposition sont, par exemple, des vitesses de 
déplacement, des vitesses de rotations, certains angles pris par le projectile au cours de sa trajectoire, et la 
position du centre de gravité. Ces données peuvent être obtenues au moyen d’un radar Doppler, de panneaux 
de tir placés le long de la trajectoire du projectile, ou de magnétomètres [1] [2]. 
Différentes techniques d’identification des coefficients aérodynamiques d’un projectile ont été 
précédemment développées. La plupart sont basées sur l'optimisation au sens des moindres carrés [3], [4], 
[5], [6] de sorte à minimiser l’écart entre des mesures et une trajectoire calculée. A titre d’exemple, Chapman 
et Al. [3], pionniers dans l'identification de paramètres aérobalistiques à partir de mesures, associent 
l'optimisation au sens des moindres carrés à une technique de corrections différentielles dans un processus 
itératif. Ces travaux sont uniquement valables dans le cadre de petites perturbations, c’est-à-dire lors de petits 
angles d’attaque. Afin d’éviter une intégration numérique et les problèmes de convergence associés, Weiss et 
Al. [5] ont proposé d'associer la technique des moindres carrés à une résolution analytique des équations 
différentielles du mouvement à six degrés de liberté. Cette technique fournit des résultats probants à 
conditions de rester dans le cas où les effets aérodynamiques non linéaires sont ignorés (effet Magnus). 
20ème Congrès Français de Mécanique                                                                  Besançon, 29 août au 2 septembre 2011 
  2 
D’autres méthodes ont été proposées, elles consistent par exemple à employer un filtre de Kalman [7] [8] 
dans le but de limiter l’imprécision des estimations dues au bruit relatif à l’instrumentation. 
La méthode inverse proposée consiste à minimiser au sens des moindres carrés une fonctionnelle contenant 
deux termes. Le premier terme caractérise un écart entre  les données mesurées et des quantités qui vont être 
calculées. Cela permet, de relaxer les mesures, qui peuvent être bruitées, et de recalculer les paramètres 
d’état du système. Le second terme fait intervenir l'intégration des équations de la mécanique du vol et 
impose ainsi que ces équations soient vérifiées au mieux. La minimisation de cette fonctionnelle conduit à un 
système matriciel non linéaire. Dans une première étape, des données simulées (pouvant être bruitées) par un 
modèle direct à six degrés de liberté sont utilisées en entrée de la technique d'identification. A l’issue du 
calcul, la technique développée fournit les coefficients aérodynamiques et permet également de reconstruire 
la trajectoire suivie par le projectile (vitesses de déplacement et de rotation). 
Dans une première partie, la mise en équation du vol du projectile par un modèle dynamique à six degrés de 
liberté est présentée. Ensuite la démarche inverse utilisée pour l'identification des paramètres est détaillée. 
Dans une dernière partie, les premiers résultats numériques obtenus sont présentés et critiqués.  
2 Mise en équation du vol par un modèle à six degrés de liberté 
Un modèle à six degrés de liberté peut être utilisé pour représenter le comportement du projectile sur sa 
trajectoire. Les équations du mouvement, pour un projectile à symétrie de révolution non balourdé, s'écrivent 
alors sous la forme suivante : 
 
( )GCFFF
m
1
v KLD ++++=ɺ  (1) 
 
( )
( )
1 c
2 a c b 2 1 er ar tan a tan mag
2 b c a 1 2
I
I I I M M M M M
I I I
 ω
 
ω − ω ω − = + + + + 
 
ω − ω ω −  
ɺ
ɺ
ɺ
 (2) 
 
0 c 1 a 2 b 3
c 0 a 3 b 21
c 3 a 0 b 12
c 2 a 1 b 03
1
2
 λ
−ω ⋅λ − ω ⋅λ − ω ⋅λ 
   +ω ⋅λ − ω ⋅λ + ω ⋅λλ   
=   +ω ⋅λ + ω ⋅λ − ω ⋅λλ   
 
−ω ⋅λ + ω ⋅λ + ω ⋅λ  λ 
ɺ
ɺ
ɺ
ɺ
 (3) 
où : 
 vɺ  est le vecteur des accélérations du projectile, cωɺ , aωɺ , bωɺ  sont les dérivées des vitesses de 
rotation, et 0λɺ , 1λɺ , 2λɺ , 3λɺ  sont les dérivées des éléments du quaternion ; 
 I1 et I2 sont les inerties longitudinale et transversale du projectile ; 
 m est la masse du projectile ; 
 DF , LF  et KF  sont les forces aérodynamiques s’exerçant sur le projectile. A chacune d’elles est 
associé un coefficient aérodynamique ; 
 C  et G  sont les forces de Coriolis et de pesanteur ; 
 erM , arM , tanM , tanaM , magM  sont les moments aérodynamiques s’exerçant sur le projectile 
auxquels sont associés différents coefficients aérodynamiques. 
Les forces sont exprimées dans un repère fixe lié au canon alors que les moments sont exprimés dans un 
repère mobile lié au projectile. Plus particulièrement, les expressions des forces et moments aérodynamiques 
sont les suivantes : 
Force de traînée : t C S q F xD −=  (4) 
Force de portance : ( )( )L zF  q S C  t c tα= ∧ ∧  (5) 
20ème Congrès Français de Mécanique                                                                  Besançon, 29 août au 2 septembre 2011 
  3 
Force de Magnus ( )cK yp
r
D 
F  q S C  c t
v
α
ω
= ∧  (6) 
Moment d’entrainement en roulis : cCD S qM 0ler =  (7) 
Moment d’amortissement en 
roulis : 
2 c
ar lp
r
M q S D  C  c
v
ω
= −  (8) 
Moment de tangage : ( )tan mM q S D C  c tα= − ∧  (9) 
Moment d’amortissement en 
tangage : ( )cc C v
D
 S qM mq
r
2
tana ɺ∧=  (10) 
Moment de Magnus : ( )( )2cmag np
r
 D
M q S  C  c c t
v
α
ω
= − ∧ ∧  (11) 
où : 
 vr est la norme de la vitesse relative du projectile (c’est-à-dire la vitesse par rapport au vent) ; 
 Le vecteur 
r
r
v
v
t =  porte le vecteur vitesse du projectile, et le vecteur c  porte l’axe longitudinal du 
projectile ; 
 
2
rv2
1q ρ=  où ρ est la masse volumique de l’air ; 
 
4
πD
=S
2
est le maître couple du projectile, où D représente le calibre du projectile. 
Le modèle à six degrés de liberté peut alors s’écrire, de manière simplifiée, sous la forme d’un système 
d’équations différentielles non linéaires correspondant aux 10 équations provenant de (1), (2) et (3) : 
 ( )( )t,C,tUf=Uɺ  (12) 
Où Ti j k c a b 0 1 2 3U = v , v , v , ω , ω , ω , λ , λ , λ , λ    est le vecteur des paramètres d’état du système, 
[ ]Tmqnpαmαl0lpypαzαx C,C,C,C,C,C,C,C=C  est le vecteur des coefficients aérodynamiques 
recherchés. Il est également possible de prendre en compte les équations des positions. 
3 Formulation du problème d’optimisation 
Lors d'un vol, des mesures discrètes de vitesses de déplacement du projectile et de vitesses de rotation sont 
enregistrées. L'idée de la technique d'identification est de chercher l'ensemble des coefficients 
aérodynamiques qui permet de s’approcher au mieux des mesures tout en vérifiant à chaque instant les 
équations de la mécanique du vol. Pour cela un problème d'optimisation non linéaire en dimension finie est 
introduit : 
 
Trouver ( )1 n= U ,..., U ,Cφ  ∈  10n 8+ℝ  tel que : 
( ) ( )J     Jφ ≤ ψ  10n 8+∀ψ ∈ℝ  avec 
( ) ( ) 22d
M E
J    R  ψ = ψ − φ + α ψ  
(13) 
Où dφ  est un vecteur composé de l’ensemble des mesures discrètes, Tm m m m mi j 0 3U v , v ,..., ,..., = λ λ   est un 
vecteur composé des valeurs des paramètres d’état à l’instant tm et α est un paramètre de pénalisation. 
Le premier terme de la fonctionnelle J caractérise un écart entre les données mesurées et des quantités qui 
vont être calculées. Cela permet de relaxer les mesures, qui peuvent être bruitées, et de recalculer l'état du 
système. Le second est un terme de pénalisation qui prend en compte l'intégration numérique des équations 
20ème Congrès Français de Mécanique                                                                  Besançon, 29 août au 2 septembre 2011 
  4 
de la mécanique du vol (cf. équation 12). Il permet ainsi d'imposer aux paramètres d'état de vérifier aux 
mieux ces équations. Par exemple lorsque que l'on utilise le schéma d'Euler explicite : 
 ( )
2n 12
m 1 m m m 1 m
E
m 1
 R  U U f (U ,C) (t t )
−
+ +
=
 ψ = − − −
 ∑  (14) 
D'autres schémas d’intégration peuvent être utilisés comme par exemple les méthodes multipas d’Adams 
Bashforth explicites et ou d’Adams Moulton implicites ou des méthodes de type B.D.F. (Backward 
Differentiation Formulas) implicites. 
Le calcul du gradient de la fonctionnelle J par rapport à chacune des composantes du vecteur ψ  fournit un 
système d’équations non linéaires que l'on peut écrire sous la forme : 
 
( ) 0C,U,...,UJ
j
n1
=
ψ∂
∂
 ∀ j = 1 , (10n+8) (15) 
Nous utilisons ensuite une technique de Newton-Raphson pour résoudre ce système d'équations non linéaires 
et obtenir un élément optimal φ . 
4 Résultats numériques 
En premier lieu, pour valider le code de calcul, nous avons utilisé des données simulées générées par un code 
direct à six degrés de liberté. Une restriction de cette solution est ensuite utilisée pour générer des données 
d’entrée qui peuvent être entachées de bruit. Dans la mesure où la méthode est en cours de développement, 
les coefficients aérodynamiques sont progressivement activés. Nous considérons aussi que l’outil développé 
est destiné à être utilisé par des personnes expertes, c’est-à-dire ayant une connaissance approximative des 
coefficients aérodynamiques recherchés. C’est pourquoi, pour valider les codes de calcul, les valeurs initiales 
des paramètres C  sont fixées à deux fois leurs valeurs théoriques. 
Le premier cas de test présente les résultats obtenus suite à l’activation de 4 coefficients aérodynamiques 
T
x lp l0 mqC = C , C , C , C   , en présence de bruit. Différents schémas d’intégration sont employés et 
confrontés. 
Le bruit appliqué sur chaque paramètre d’état est exprimé en pourcentage de la valeur maximale prise 
pendant le vol. Le bruit sur les données est fixé à 1% pour vi, et 5% pour chacune des vitesses de rotations 
ωc, ωa, ωb. Le tableau suivant présente le pourcentage d’erreur relative sur les coefficients aérodynamiques 
identifiés par rapport à leurs valeurs exactes. Le pas de temps est fixé à 7.10-4 secondes. 
Erreur relative sur les coefficients aérodynamiques identifiés (%) Schéma 
d’intégration 
Cx Clp Cl0 Cmq 
B.D.F.1 0.1 3.2 2.6 867.6 
B.D.F.2 2.2 2.1 1.5 49.6 
B.D.F.3 0.7 0.7 1.2 38.9 
B.D.F.4 0.28 3.2 2.4 10.0 
B.D.F.5 1.0 11.2 4.8 2.1 
B.D.F.6 2.2 0.8 0.4 0.3 
La méthode B.D.F. à l’ordre 6 permet d’obtenir le meilleur compromis sur les erreurs relatives obtenues pour 
l’ensemble des coefficients aérodynamiques identifiés. Toutefois, un ordre inférieur à 6 est suffisant si l’on 
souhaite identifier individuellement certains coefficients. Les courbes suivantes présentent la reconstruction 
des données de vol en utilisant le schéma à l’ordre 6. Les paramètres d'état recalculés sont correctement 
reconstruits et permettent d’accéder à la trajectoire du projectile au cours du vol. Les courbes sont 
adimensionnées et certaines échelles sont absentes pour des raisons de confidentialité. 
 
20ème Congrès Français de Mécanique                                                                  Besançon, 29 août au 2 septembre 2011 
  5 
 
FIG. 1 – Reconstruction des paramètres d’état du problème 
Le second cas présente les résultats obtenus après l’activation de 6 coefficients aérodynamiques 
T
x lp l0 m np mqC C , C , C , C , C , Cα α =   , en utilisant le schéma d’intégration explicite d’Euler. Le pas de 
temps est fixé à 10-5 secondes. Le tableau ci-dessous présente les pourcentages d’erreur relative sur les 
coefficients aérodynamiques identifiés. 
Erreur relative sur les coefficients aérodynamiques 
identifiés (%) Bruit sur vi 
(%) 
Bruit 
sur ωc 
(%) 
Bruit 
sur ωa 
(%) 
Bruit 
sur ωb 
(%) Cx Clp Cl0 Cmα Cnpα Cmq 
Nombre 
d’itérations 
0 0 0 0 0.000052 0.0066 0.0066 0.032 25.8 0.44 13 
0.1 0 0 0 2.5 0.3 0.013 0.064 21.8 0.5 13 
0 0.1 0 0 0.00051 0.13 0.028 0.028 25.3 0.53 13 
0 0 0.1 0 0.000052 0.0066 0.0066 0.24 152.9 8.2 15 
0 0 0 0.1 0.000052 0.0066 0.0066 0.57 32.1 0.3 15 
0.1 0.1 0.1 0.1 2.5 0.3 0.0057 2.9 204.4 11.3 14 
Il ressort de cette simulation que les coefficients Cx, Clp, Cl0, Cmα et Cmq sont identifiés avec une précision 
satisfaisante en présence de bruit. En revanche, le coefficient Cnpα est plus délicat à identifier. Ce coefficient 
est couplé à la fois à la vitesse de roulis ωc et aux éléments du quaternion λ0, λ1, λ2 et λ3 (cf. équation 11). 
Ces grandeurs évoluent très rapidement au cours du vol. De plus nous n'avons pas de mesure des 
quaternions. D’autre part, en fixant le pas de temps à 10-5 secondes, le nombre de données du vol est 
tellement important qu’elles ne peuvent pas être toutes prises en compte. Seul le début du vol a donc été 
traité. Or, sur cette fenêtre de vol, la vitesse de roulis ωc n’a pas encore atteint sa valeur maximale (constante 
20ème Congrès Français de Mécanique                                                                  Besançon, 29 août au 2 septembre 2011 
  6 
en fin de vol, cf. FIG. 1). Une prise en compte de données situées en fin de vol permettrait peut être d’avoir 
de meilleurs résultats sur le coefficient Cnpα. 
5 Conclusion 
Une technique d'identification des coefficients aérodynamiques d'un projectile flèche à partir de mesures 
enregistrées en vol a été proposé. A l'heure actuelle, la technique a été développée avec six coefficients 
activés et permet d'identifier avec une précision certaine, même lorsque les données sont bruitées, cinq 
d'entre eux. 
Les prochains développements vont tout d'abord consister à prendre en compte les deux équations de la 
mécanique du vol restantes par la méthode d’Euler explicite afin de déterminer les derniers coefficients 
aérodynamiques (Czα et Cypα). Puis sur le système complet, différents schémas d'intégration numérique seront 
testés afin de pouvoir gérer des pas de temps plus grands. Les différents schémas numériques seront alors 
comparés en termes de précision et de stabilité vis à vis de l'initialisation ou de données bruitées. Finalement, 
la méthode devra être validée avec les données réellement enregistrées lors d'un vol. 
Références 
[1] Berner C., Fleck V., Sommer E., Lentz A., Dupuis A., Comparison of simultaneous free flight 
measurements between 3D magnetic sensors and aeroballistic range data, French German Research Institute 
(ISL), PU 627/2004, 2004 
[2] Fleck V., Sommer E., Berner C., Dupuis A., Aerodynamic characteristics of a long range spinning 
artillery shell obtained from 3D magnetic sensors, 22nd International Symposium on Ballistics, Vancouver 
BC, Canada, 2005 
[3] Chapman G. T, Kirk D.B., A method for extracting aerodynamic coefficients from free-flight Data, 
AIAA Journal 0001-1452, 8, 753-758, 1970 
[4] Lieske R. F., Mackenzie A. M., Determination of aerodynamic drag from radar data, Memorandum rept, 
rapport n° 0465057, Aberdeen Proving Ground Tech. Rep., 1972 
[5] Weiss S., Doherr K-F., Schiling H., Analytical solution and parameter estimation of projectile dynamics, 
Journal of spacecraft and rockets 32, 67-74, 1995 
[6] Mahmutyazicioğlu G., Platin B. E., Aerodynamic data identification using linear aeroballistic theory, 
RTO SCI Symposium, 1998 
[7] Quanwei J., Qiongkang C., Dynamic model for real-time estimation of aerodynamic characteristics, 
Journal of Aircraft 1989 0021-8669, 26, 315-321, 1986 
[8] Yamada H., Miyata T., Ichikawa H., Measurement of aerodynamic coefficients by system identification 
methods, Journal of Wind Engineering and Industrial Aerodynamics, 42, 1255-1263, 1992 


Identification inverse des coefficients aérodynamiques d'un engin hypersonique à partir de données issues de vol

http://documents.irevues.inist.fr/bitstream/2042/46313/1/cfm2011_487.pdf

Identification inverse des coefficients aérodynamiques d'un engin hypersonique à partir de données issues de vol

Abstract

Similar works

Full text

Available Versions

I-Revues