Allumage des moteurs fusées cryotechniques
Allumage des moteurs fusées cryotechniques
THÈSE
présentée pour obtenir le titre de
DOCTEUR DE L’INSTITUT NATIONAL POLYTECHNIQUE DE TOULOUSE
spécialité : Dynamique des fluides
Antoine Dauptain
Réf. CERFACS/TH/CFD/06/85
A Guenaëlle
Sommaire
1 Introduction 5
1.1 Histoire des moteurs fusées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Marché actuel des lanceurs spatiaux . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3 Cadre de la thèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2 Etat de l’art 23
2.1 Description de l’allumage d’un moteur fusée . . . . . . . . . . . . . . . . . . . . . 23
2.2 La chimie de l’auto-allumage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.3 La dynamique des jets supersoniques sous-détendus . . . . . . . . . . . . . . . . . 32
2.4 La combustion supersonique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.5 Stratégie de la thèse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3 Equations 39
3.1 Equations d’écoulements réactifs . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.1 Equations de conservation . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1.2 Variables thermodynamiques . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.1.3 Equation d’état des gaz parfaits . . . . . . . . . . . . . . . . . . . . . . . 41
3.1.4 Conservation de la masse et vitesse de correction . . . . . . . . . . . . . . 42
3.1.5 Coefficients de transport . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.1.6 Flux de chaleur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.1.7 Cinétique chimique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.2 Equations SGE d’écoulements réactifs . . . . . . . . . . . . . . . . . . . . . . . . 45
3.2.1 La Simulation aux Grandes Echelles (SGE) . . . . . . . . . . . . . . . . . 45
3.2.2 Equations SGE pour les milieux réactifs . . . . . . . . . . . . . . . . . . . 46
3.2.3 Modèles pour le tenseur de sous-maille . . . . . . . . . . . . . . . . . . . 49
3.3 Présentation du code AVBP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.3.1 Discrétisation Cell-Vertex . . . . . . . . . . . . . . . . . . . . . . . . . . 51
1
2 SOMMAIRE
4 Auto-allumage 57
4.1 Configurations de références . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1.1 Auto-allumage homogène . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4.1.2 Auto-allumage en diffusion . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.1.3 Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.2 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7 Conclusion 179
7.1 Récapitulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179
7.2 Analyse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180
7.3 Perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 181
8 Remerciements 183
9 Annexe 1 185
9.1 Lois du mouvement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185
9.1.1 Mouvement dans le champ gravitationnel terrestre loin du sol . . . . . . . 185
9.1.2 Mouvement dans le champ gravitationnel terrestre près du sol . . . . . . . 187
9.2 Loi du vol d’une fusée . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 187
10 Annexe 2 191
10.1 Méthode des caractéristiques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191
10.2 Taille des bouteilles de Mach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
10.3 Analyse du champ lointain d’un jet sous-détendu . . . . . . . . . . . . . . . . . . 199
10.3.1 Diamètre virtuel du jet . . . . . . . . . . . . . . . . . . . . . . . . . . . . 199
10.3.2 Vitesses et concentrations axiales . . . . . . . . . . . . . . . . . . . . . . 201
11 Annexe 3 205
11.1 Contrainte du coût de calcul . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
11.2 Panorama des moyens de calculs . . . . . . . . . . . . . . . . . . . . . . . . . . . 205
11.3 Coût des codes de calculs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 206
11.4 Coût des SGE présentées . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 207
4 SOMMAIRE
Chapitre 1
Introduction
Cette introduction décrit les contraintes relatives à la construction des moteurs fusées cryotech-
niques, et identifie le besoin d’étudier en détail leur allumage. Pour commencer, quelques éléments
de l’histoire des fusées établissent les contraintes techniques incontournables. Ces contraintes sont
à replacer dans le panorama économique actuel pour mettre en évidence les besoins particuliers des
futurs lanceurs spatiaux. L’ensemble de ces élements forment le cadre général de cette thèse.
La première mise en évidence indiscutable du principe de propulsion par réaction date de 100
av. J.C (c.f. (134), Mediteranean Civilizations and Middle Ages). A cette époque Hero d’Alexandrie
invente l’aeolipile (litt.’boule à vent’), une sorte de roue à réaction illustrée par la Fig. 1.1. Un
chaudron clos produit de la vapeur s’évacuant par une sphère. Cette sphère peut tourner autour de
l’axe d’alimentation en vapeur. Les tuyaux d’évacuation de la sphère mobile étant coudés, l’éjection
de vapeur provoque un mouvement de rotation. Cette invention n’a cependant aucune application
pratique.
Origine Chinoise, de -200 à 1232
Dès 200 av. J.C., des chinois connaissent empiriquement les propriétés pyrotechniques des mélange
de poudre de charbon, de soufre et de salpêtre (c.f. Fig 1.8,
). 1 Sur cette base de chimie, ils
développènt l’art des feux d’artifices. Ces feux sont utilisés pour les rituels religieux, dans le
but de chasser les mauvais esprits. La poudre noire, ou plutôt l’utilisation du salpêtre en com-
bustion, se répand en occident au e siècle grâce au commerce. Les Byzantins se servent des
propriétés du salpêtre, en l’ajoutant à leurs feux grégeois pendant les batailles navales. Pour une
rétrospective exhaustive de l’histoire de la poudre noire, on peut se référrer aux travaux de M.
Mercier (105). Concernant les fusées proprement dites, on préférera l’analyse historique précise de
Jixing Pan (113). Pour résumer son propos, les historiens hésitent entre plusieurs dates de mise en
service de projectiles auto-propulsés.
Il y a tout d’abord en 1161 le p’i-li-pao (fusée Tonnerre, Fig. 1.2), un cylindre de papier conte-
nant deux compartiment de poudre delimités par trois couches de terre cuite. La première charge
5
6 Chapitre 1. Introduction
de propulsion est forée sur toute sa longueur. La deuxième contient un agent fumigène. Ces fusées
permettent de dresser des écrans de fumée à distance. Les témoignages écrits comme Hai-ch’iu fou-
hsü (Postface de l’Ode aux Navires de Guerre Anguilles) par Yang Wan-li en 1170, et les techniques
de l’époque, depuis la fabrication de la poudre jusqu’aux techniques de forages d’un pain de poudre
comprimé, attestent que cet engin peut avoir été utilisé en 1161. Un doute reste quant au rôle de
la couche de boue inférieure. Elle peut être une simple protection avant lancement, ou jouer le rôle
de tuyère rudimentaire, accélérant les gaz par effet Venturi. Seul Jixing Pan la décrit comme une
tuyère, c’est pourquoi cette conclusion reste à étayer.
On trouve ensuite en 1232 la célèbre Fei-huo-ch’iang (Lance de feu volante, Fig. 1.2). La traduc-
tion controversée de Chin-Shih (Histoire de la dynastie Chin), écrite par T’o T’o et Ouyang Hsüan
en 1345, laisse place à de nombreuses ambigüités. En effet, lors de la description de la victoire de
l’armée Chin conte le siège des troupes mongoles autour de la ville de K’ai-feng-fu en 1232, il est
question de lances de feu volantes. Ce sont des lances de grande taille (approximativement deux
mètres) auxquelles sont fixées une cartouche de poudre de ’deux pieds de long’ (approximativement
50 cm). Aucun système de lancement (arbalète, catapulte) n’étant décrit, ceci fait penser naturelle-
ment aux fusées. En revanche, le texte parle de ’répandre le feu sur une distance de dix pas’, une
portée trop courte pour une fusée de cette taille, à moins qu’il ne s’agisse du panache d’éjection.
Enfin, le texte ajoute que la cartouche contenant la poudre est intacte aprés combustion, une dif-
ficulté technique inutile pour un projectile explosif ou incendiaire, plaidant en faveur du caractère
propulsif de la poudre. La Fei-huo-ch’iang est donc probablement la première fusée au monde à
emporter une charge utile conséquente.
iron head
top
fuse
fuse
paper tube
power column
wooden stick
propellant
nozzle
fuse
L’Empire du Milieu est cependant limité par une barrière technique de l’ordre de la dizaine
de kilogrammes de charge utile : au delà d’une certaine quantité de poudre noire, la combustion
détruit son support. Inversement, un support suffisament solide possède un poids bien supérieur à la
8 Chapitre 1. Introduction
F. 1.3 – L’idéogramme japonais feu - flèche signifie ’fusée’, un dispositif régulièrement utilisé
autour de 1200 ap. J.C. à des fins religieuses et militaires. Dessin postérieur au e siècle d’un
soldat japonais lançant un fusée militaire.
1.1. Histoire des moteurs fusées 9
poussée disponible. Il est donc impossible de passer à des applications de grande échelle, comme
la construction d’armes de destruction plus massive ou la propulsion d’un homme (cf. Fig. 1.4). On
voit déjà se poser ici le problème du rapport de masse mis en évidence à l’aube du e siècle : la
charge utile et la distance à parcourir définissent la masse de carburant nécessaire.
De l’Orient à l’Occident, fin du e siècle
L’efficacité militaire des fusées étant prouvée lors du siège de K’ai-feng-fu en 1232, la technologie
des fusées de combat se propage rapidement dans le monde pendant le e siècle (c.f. Fig 1.8,
). 2
On sait que les Mongols attaquent le Japon par ce moyen en 1274 et 1281. Le Japon, le Java, la
Corée et l’Inde intègrent vite la nouvelle arme dans leur arsenal. A des milliers de kilomètres de
là, l’empire arabo-musulman collecte et enrichit depuis six siècles les sciences des civilisations voi-
sines, notamment des empires grec, égyptien, perse, indien et chinois (c.f. Fig 1.8,
). 3 Les liens
commerciaux soutenus, en particulier par la route de la soie, induisent des échanges culturels et
techniques rapides. La manne de savoir ainsi obtenue, fruit d’une volonté politique initiée par les
califes Al Malik (685-705) et Al Mamun (813-833), permet de déployer les fusées militaires arabo-
musulmanes moins de 50 ans après l’épisode du siège de Kai-Fung-Fu, comme l’attestent les écrits
de l’érudit syrien al-Hasan al-Rammah en 1285. Leurs alchimistes améliorent le processus de fa-
brication de la poudre noire en utilisant des composés purifiés. Le produit de leurs recherches est un
explosif rapide autorisant les premières applications ballistiques.
A l’aube du e siècle, l’Europe commence enfin à combler son retard sur l’Orient avec les
travaux sur la poudre noire de Roger Bacon, Albertus Magnus et Marcus Gaecus (c.f. Fig 1.8,
). 4
Un poête distrait nommé Wan Hu souhaitait par dessus tout explorer les cieux et rencontrer les étoiles. Aprés avoir
escaladé des montagnes et même essayé des catapultes, il construisit une chaise de bambou équipée de cerfs volants et
de 47 fusées. Lorsque la fumée se dissipa, il ne restait aucune trace de Wan-Hu, ni de son véhicule.
F. 1.4 – Un conte situe le premier homme à voler dans l’espace en l’an 1500. Illustration contem-
poraine d’un livre de contes pour enfants
1.1. Histoire des moteurs fusées 11
F. 1.5 – La civilisation arabo-musulmane joue un rôle essentiel dans la diffusion de la science
au moyen-âge. Le cas des fusées en est une illustration parfaite. Peinture d’influence européene
(perspective) d’une scène turque ou perse (Arche à voute circulaire bicolore).
12 Chapitre 1. Introduction
obtenues presque simultanément par ces trois hommes, est présentée en Annexe I.
science et guerre donnent lieu à bien des réflexions intéressantes à lire, comme celle de Edgarton
dans ”Science and War” (111).
F. 1.7 – Wernher von Braun, qui construisit le V2, puis contribua directement aux programmes
d’exploration spatiale et de défense des Etats Unis d’Amérique, fut au centre d’un débat éthique
toute sa vie.
La guerre froide
La conquête spatiale est souvent présentée comme une des plus belles pages de l’histoire humaine.
La légende bien lisse de ’l’étoffe des héros’, ces explorateurs intrépides repoussant les limites hu-
maines, ne doit cependant pas occulter que le lancement de Spoutnik ou les premiers vols des
hommes dans l’espace et vers la Lune sont autant de manières de montrer la maı̂trise de chaque
camp dans le domaine des lanceurs (c.f. Fig 1.8,
). 9 Le développement des théories de dissua-
sion issues de la crise de Cuba en 1962, induisant l’escalade dans l’équilibre de la terreur, exige
des lanceurs toujours plus rapides à mettre en œuvre. L’objectif militaire est de convaincre l’ad-
versaire qu’un missile balistique intercontinental peut décoller en quelques minutes (d’où le nom
du missile américain Minuteman1 ). C’est pourquoi les recherches sont développées en direction
des carburants dit ”stockables”. Ces carburants ont l’avantage de simplifier la préparation d’une
fusée, mais possèdent une très haute toxicité. Le risque technique (et environnemental) soulevé par
ces réactifs fut cependant accepté dans les deux camps, malgré l’opposition farouche de plusieurs
spécialistes comme S.P. Korolev, le célèbre ’père des fusées russes’. La technologie des carburants
stockables correspond aujourd’hui aux divers propulseurs à poudre encore largement utilisés. Bien
que ces réactifs soient particulièrment stables, leur synthèse fait toujours intervenir des gaz haute-
ment toxiques comme le Phosgen. La guerre froide se termine dans les années 80 avec la politique
d’ouverture de Michaël Gorbatchev, la perestroı̈ka. C’est alors l’avènement de l’économie spatiale,
décrite dans la prochaine section.
Vue d’ensemble
1
Lors de la guerre de sécession americaine, les minutemen étaient des volontaires prêts à rejoindre l’armée au premier
appel, en moins d’une minute.
14 Chapitre 1. Introduction
Tous ces événements possédent une organisation intéressante. En effet, dans une première période,
de -500 à 1900, plusieurs civilisations acquièrent successivement la technologie des fusées de com-
bat. A chaque fois, la technologie de la poudre apparaı̂t comme un prérequis. Par exemple, les Grecs
possédaient presque le principe du moteur fusée avec la vapeur, mais ne pouvaient le mettre en pra-
tique faute d’un carburant suffisamment efficace. Cela corrobore la prépondérance de la nature du
carburant, identifiée plus tard par K.E. Tsiolkovski dans se définition de l’impulsion spécifique. Par
la suite, un grand regroupement scientifique permet d’obtenir la technologie elle même, et le cycle
se termine par un document écrit prouvant et datant la maı̂trise de cette nouvelle technologie. On
pourra noter le rôle essentiel du Moyen Orient dans la propagation de la science en général et des
fusées en particulier.
La deuxième période, de 1900 à 1980 concerne les lanceurs de haute altitude. Les moyens fi-
nanciers et technologiques nécessaires limitent le nombre de nations capables de construire ces
engins. Les grandes guerres du e siècle conditionnent la recherche, effectuée sous forme de
compétition parallèle, principalement entre les blocs de l’Est et de l’Ouest de 1945 à 1980. Ce
contexte d’émulation et d’oppositions a mis en évidence les problèmes éthiques liés aux fusées, à
savoir les applications militaires toujours proches et la pollution colossale provoquée par les car-
burants. La troisième période, de 1980 à aujourd’hui, est caractérisée par l’ouverture mondiale
(c.f. Fig 1.8,
)
10 qui succéde à la guerre froide. Le nombre de pays lanceurs de fusées spatiales s’est
considérablement élargi et les montages internationaux se sont développés. Les transferts de tech-
nologies ne sont plus rares, bien qu’ils ne soient pas encouragés par les nations, comme l’a prouvé le
Congrés Américain en condamnant en 2003 Boeing et Hughes à verser 32 000 000$ d’amende suite
à des transferts illégaux de technologies vers la Chine dans les années 90. Boeing avait déjà dù payer
en 1998 une amende de 10 000 000$ pour plusieurs centaines de délits ”d’exportations de services
ou d’articles de défense sans autorisation” dans la création du montage international Sealaunch, de
1994 à 1998. Cette nouvelle ère complexe de partenariats est la conséquence de la priorité donnés à
la composante économique de l’industrie spatiale, qui va être abordée dans la partie suivante.
1.1. Histoire des moteurs fusées 15
temps
10
1981 Navette Spatiale
9 1946−1981
1975 50 ans
1969 Armstrong sur Lien scientifique
la lune Guerre froide
1962 Glenn en orbite 1961 Gagarine en orbite
1957 Spoutnik en orbite Ecrits
1950
Equipe Equipe Deuxieme 1939−1945 Equipe Equipe Poudre noire
americaine allemande
8 allemande russe
guerre mondiale
Regroupement
Verein fur 1927−
1925 1934 scientifique
Raumshiffart
1922 1923
6 Method for reaching 7 Die rakete su den Guerres
Extreme altitudes planetenraumen
1900
5 Travaux de
1900 K.E. Tsiolkovsky
1687
Siecle des Philosophiae Naturalis
Principa Mathematica
Lumieres
1500 Empire
1300
4 Poudre noire Arabo−musulman
Chine
Hassan 1285
Al Rammah Siege de 1232
2
Croisades Kai Fung Fu
1000 Initiatives des califes
Al Malik & Al Mamum 1045
685−705 813−833 Tseng Lung Kiang
673
Feux grecs
500
3
1
0 −200
Feux d’artifice
Science Science Science Science
Grecque Egyptienne Indienne Chinoise
−500
F. 1.8 – Chronologie et répartition mondiale des principaux événements constituant l’histoire des
moteurs fusées.
16 Chapitre 1. Introduction
une orbite géostationnaire (lancement GTO2 ) par une fusée chinoise Longue Marche 3B revient
à 60 M$, tandis que quatre lancements par des fusées chinoises Longue Marche 2C reviennent à
90 m$, soit 50% plus cher avec une technologie identique. Cela explique l’existence des lanceurs
’gros porteurs’ comme Delta IV Heavy (Boeing, 13 tonnes en orbite géostationnaire), Ariane 5 ECA
(Arianespace, 10 tonnes en orbite géostationnaire) illustrés en Fig. 1.11. Pour soulever de telles
massses, on utilise de préférence des moteurs cryotechniques Hydrogène liquide - Oxygène liquide,
comme les quatre moteurs de Delta IV Heavy . Le problème commercial qui en découle se situe dans
le carnet de commande : rares sont les missions qui permettent de tels regroupements de satellites. Si
le lanceur est capable de varier les orbites d’injections des charges utiles lors d’une même mission,
les possibilités de missions se multiplient. C’est la raison d’être du projet de moteur de deuxième
étage Vinci, réallumable cinq fois en cours de vol. Cela étant, l’allumage d’un moteur est l’une des
actions les plus critiques pour le lanceur, causant de nombreux accidents. Ainsi, la baisse du prix
d’accès à l’espace est étroitement liée à la maı̂trise des phénomènes physiques mis en jeu lors
de l’allumage des moteurs cryotechniques.
2
L’orbite géostationnaire est la plus éloignée de la terre, et donc la plus coûteuse des missions de lancement en
orbite terrestre. Le lanceur place la charge utile sur une trajectoire elliptique excentrée appelée Orbite de Transfert
Géostationnaire ou GTO. L’impulsion finale vers l’orbite circulaire est réalisée par le satellite.
18 Chapitre 1. Introduction
Nom Longue Dnepr Delta 2 Soyuz Atlas 2AS Longue Ariane 44L
Marche 2C Marche 2E
Origine China Russie USA Russie USA China Europe
Capacité LEO 3 200 Kg 4 400 Kg 5 144 Kg 7 000 Kg 8 618 Kg 9 200 Kg 10 200 Kg
Capacité GTO 1 000 Kg - 1 800 Kg 1 350 Kg 3 719 Kg 3 370 Kg 4 790 Kg
Prix lancement 22.5M$ 15M$ 55M$ 37.5M$ 97.5M$ 50M$ 112.5M$
Prix LEO 7 031$/Kg 3 409$/Kg 10 692$/Kg 5 357$/Kg 11 314$/Kg 5 435$/Kg 11 029$/Kg
Prix GTO 22 500$/Kg - 30 556$/Kg 27 778$/Kg 26 217$/Kg 14 837$/Kg 23 486$/Kg
T. 1.1 – Lanceurs légers, moyens et lourds, par ordre croissant de charge utile en orbite basse.
GTO et LEO désignent respectivement des missions vers des orbites géostationnaires (36 000 Km
d’altitude) et vers des orbites basses (400 km). Source : documents publics Futron Corp. Sept. 2002.
1.2. Marché actuel des lanceurs spatiaux 19
60000
50000
Lancements GTO
Autres lancements
Prix de lancement [$/Kg]
30000
20000
10000
0
1988 1990 1992 1994 1996 1998 2000 2002
Annee
F. 1.9 – Evolution annuelle du coût de lancement de satellites, évalué sur 264 lancements de sa-
tellites commerciaux, normalisé en dollars (valeur Etats-Unis en 2000). Avec l’aimable autorisation
de T. Freeland, Futron Corporation
20 Chapitre 1. Introduction
F. 1.10 – Les lanceurs transnationaux : Sea Launch est un partenariat Boeing (USA), Kvaerner
ASA (Norvège), RSC-Energia (Russie) et SDO Yuzhnoye/PO Yuzhmash (Ukraine). Starsem,
ou Soyuz Company, est un partenariat EADS (Europe), Arianespace (Europe), Russian Federal
Space Agency (Russie) et The Samara Space Center (Russie).
1.2. Marché actuel des lanceurs spatiaux 21
F. 1.11 – Les lanceurs gros porteurs Delta IV ’Heavy’ de Boeing (USA) et Ariane V ECA
d’Arianespace (Europe).
22 Chapitre 1. Introduction
Etat de l’art
Plusieurs disciplines techniques interviennent dans l’allumage d’un moteur fusée cryothech-
niques. La partie de la fusée qui nous intéresse se situe dans la chambre de combustion d’un mo-
teur à carburant liquide, depuis l’injection des ergols jusqu’au col précédent le divergent. Trois
mécanismes de base contrôlent le fonctionnement : l’auto-allumage, les jets supersoniques sous-
détendus, et la combustion supersonique. L’état de l’art dans la recherche sur chacun des ces phéno-
mènes constituent le point de départ de ce travail.
Oxygène
Allumeur
Hydrogène
Voyons tout d’abord le lieu sur lequel vont porter les recherches. Les informations techniques
qui suivent ont été fournies par les ingénieurs de la division ”grosse propulsion liquide” de la
Snecma. Le schéma de la Fig. 2.1 illustre la chambre de combustion d’un moteur cryotechnique,
23
24 Chapitre 2. Etat de l’art
depuis la plaque d’admission des réactifs, également appelés ergols, jusqu’au col de la tuyère. Le
diamètre de la chambre fait plusieurs dizaines de centimètres. Les injecteurs coaxiaux, au nombre
d’une centaine sur l’ensemble de la plaque, ont un diamètre extérieur de l’ordre du centimètre, et
l’allumeur un diamètre de quelques millimètres.
Le démarrage d’un moteur fusée se déroule dans l’ordre suivant : purge complète par un gaz
neutre, remplissage par le carburant sans combustion, mise en route de l’allumeur, et enfin admis-
sion de l’oxydant. On remarque que procédure d’allumage doit créer un mélange réducteur avant
l’allumage car un mélange oxydant attaque les parois du moteur (Fig. 2.3). Voici les quatre étapes
détaillées :
He
H2
Allumeur
O2
I II III IV Temps
He H2
He
O2
Allumeur
H2 H2
He He
F. 2.3 – Destruction par excés d’oxydant d’une chambre de combustion expérimentale. Avec l’ai-
mable autorisation de B. Berger, Swiss Propulsion Laboratory.
2.1. Description de l’allumage d’un moteur fusée 25
1. Purge : La mise en route d’un moteur fusée commence par la mise en place d’un état ’no-
minal’. C’est l’objectif de la purge (Fig. 2.2-I), qui consiste à balayer tous les systèmes par
de l’hélium (He), un gaz neutre. Lors d’un allumage au sol ce gaz chasse l’air présent dans
le moteur afin d’obtenir un milieu sans oxygène, tandis qu’en altitude cet hélium comble ra-
pidement le vide spatial, créant une faible pression avant l’arrivée des ergols. La purge fait
intervenir la mécanique des fluides et des gaz raréfiés.
2. Admission du carburant : Ensuite, le carburant, ici de l’hydrogène pur (H2 ), est injecté à
haute pression par une turbo-pompe (Fig. 2.2-II). La pression dans la chambre après purge
en altitude s’approchant de quelques centièmes de Bars, carburant liquide se détend forte-
ment, remplissant la chambre sous forme gazeuse. L’admission de carburant fait appel à
la mécanique des fluides multi-espèces appliquée à des jets supersoniques issus de fortes
détentes : les jets supersoniques sous-détendus.
La durée de l’admission du carburant, quelques dixièmes de secondes, est à comparer avec
le temps caractéristique d’injection de l’hydrogène, inférieur à quelques dixièmes de millise-
condes. Ainsi, le calcul complet d’un remplissage de carburant est excessivement long.
3. Allumage : L’allumeur se met alors en route (Fig. 2.2-III). Son rôle consiste à créer une
zone robuste de gaz chauds. Deux stratégies sont possibles. Les allumeurs pyrotechniques
(Fig. 2.5.a) sont des fusées d’ergols solides qui créent un jet chaud de produits de combus-
tions inertes. Cette technologie est bien maı̂trisée, mais s’adapte mal au réallumage. En effet,
il faudrait utiliser un barillet de plusieurs allumeurs, c’est-à-dire un système complexe et
lourd. L’alternative est un allumeur cryotechnique (Fig. 2.4 et 2.5.b), une torche créée par
la combustion des carburants ét oxydants liquides. Dans le cas d’un allumeur cryotechnique
F. 2.4 – Allumeur cryotechnique. Le diamètre de la canne d’injection est de l’ordre du centimètre.
Avec l’aimable autorisation de Stork / Aerospace Propulsion Products
(Fig. 2.5.b), l’oxygène réagit avec une partie de l’hydrogène (typiquement un tiers du débit
massique) pour créer un jet initial de gaz chauds riches en oxygène. L’autre partie de l’hy-
26 Chapitre 2. Etat de l’art
Réducteur
a) Pyrotechnique b) Cryotechnique
drogène (deux tiers du débit massique) permet de refroidir les parois de l’allumeur et brûle
dans la chambre au contact du jet initial oxydant chaud, ce qui provoque le panache d’allu-
mage. Dans son ensemble, l’allumeur est un moteur fusée stœchiométrique. La technologie
de l’allumeur utilise la chimie à cinétique finie, les écoulements supersoniques multi-espèces
réactifs couplés aux transferts thermiques et à la résistance des matériaux.
Lors de la mise en route de l’allumeur, sa pression statique en sortie est bien supérieure à
la pression chambre (NPR=50 dans les premiers instants). Cela produit un nouveau jet su-
personique fortement sous-détendu (c.f. Annexe II) qui contient un choc hydrodynamique
avec un rapport de pression statique aval/amont de plusieurs centaines. Le calcul d’une telle
discontinuité suppose l’emploi de techniques numériques spécifiques, c.f. Chap. 5.
4. Admission de l’oxydant : Enfin l’oxydant, ici de l’oxygène pur (O2 ), est admis de la même
manière que le carburant (Fig. 2.2-IV). La combustion augmente progressivement la tempéra-
ture et la pression. Le jet d’oxydant gazeux redevient liquide, puis supercritique. La présence
antérieure d’un mélange chaud doit assurer une combustion progressive, sans apparition de
poches froides de mélange carburant/oxydant. Dans le cas contraire, ces poches sont des zones
de combustion incontrôlée violente qui secouent le lanceur par des pics de pression brutaux :
c’est un allumage ’dur’.
Lorsque l’admission d’oxygène débute, la combustion principale commence. Cet ergol initia-
lement injecté à l’état liquide. Les gaz produits font augmenter alors la pression de la chambre
2.1. Description de l’allumage d’un moteur fusée 27
L’auto-allumage signifie l’établissement spontané d’une zone réactive entre deux réactifs mis en
contact. Pour la plupart des mélanges, la réaction s’effectue seulement au-dessus d’une température
critique. Généralement, les réactifs sont mélangés froids, puis le mélange est chauffé de différentes
manières pour obtenir une combustion : flamme pilote, poche de gaz chauds, diffusion et rayonne-
2.2. La chimie de l’auto-allumage 29
ment de la flamme principale, etc... Si la température à laquelle se fait le mélange est déjà au-dessus
de la température critique, alors la combustion est ”spontanée”, et on parle d’auto-allumage.
La chronologie des études menées sur l’auto-allumage, présentée en figure 2.7, fait apparaı̂tre
trois ensembles de recherches. Le premier ensemble concerne la modélisation de la cinétique chi-
mique. En 1889, Arrhé-nius (7) propose un modèle mathématique pour décrire un taux de réaction.
Les schémas cinétiques développés par la suite sont un ensemble d’étapes chimiques de type Arrhénius.
En 1980, Kee, Miller et Rupley (81) proposent un ensemble d’outils informatiques, CHEMKIN, qui
devient un standard de la discipline. Les schémas cinétiques plus ou moins complexes sont ajustés
pour simuler différentes applications. Certains schémas sont créés pour être cohérents dans plu-
sieurs configurations (flamme prémélangée, flamme de diffusion, autoallumage en tube à choc,
etc...). Concernant plus particulièrement l’hydrogène, la thermodynamique montre que les hautes
température atteintes lors de la combustion (2500 à 3000 ◦ K) autorisent la présence d’ions (H, O, OH)
en proportion non négligeable dans les produits. D’un point de vue cinétique, il apparaı̂t que les in-
termédiaires réactionnels (HO2 , H2 O2 ) prennent une part incontournable dans la combustion. Plu-
sieurs schémas globaux sont également disponibles, comme celui de Yetter (155). Les échelles de
temps ou d’espace associées aux réactions contenant les intermédiaires (HO2 , H2 O2 ) étant très pe-
tites, elles nécessitent une puissance de calcul très au-delà de ce qui est acceptable dans le cadre de
la simulation de systèmes complets. augmentent dramatiquement la puissance de calcul nécessaire.
On a donc souvent recours à des schémas simplifiés ne faisant pas intervenir ces intermédiaires ra-
pides, comme celui de Baurle (9) avec seulement 6 espèces (H2 , O2 , H2 O, OH, O, H) pour seulement
7 réactions.
Un deuxième ensemble de travaux regroupe de 1994 à 1997 l’étude de la technologie Diesel par
les motoristes (Baritaud - Institut Français du Pétrole (8), Kong - Engine Research Center (86) ). La
puissance des calculateurs augmentant, il est suivi d’une série de simulations numériques directes
d’auto-allumage (8; 86; 116; 151) mêlant cinétique chimique détaillée, turbulence résolue et pro-
cessus de diffusion complexes. Ces travaux ont montré que le temps d’apparition des premiers sites
d’auto-allumage est peu sensible à la turbulence, et qu’ils se situent sur une richesse correspondant à
celle du prémélange parfait le plus rapide à s’allumer. Sur cette iso-surface de richesse, ils débutent
là ou la dissipation scalaire est la plus basse. Enfin, ces premiers sites, apparus en temps t, transitent
vers une flamme non-prémélangée complètement turbulente pendant une durée allant de 0.8t à 1.8t.
30 Chapitre 2. Etat de l’art
T. 2.1 – Schéma cinétique de Yetter (155) pour la combustion de l’Hydrogène. Unités en cm3 −
mol − Kcal − K.
2.2. La chimie de l’auto-allumage 31
1977 Halstead(68)
Modélisation
de la chimie
1980
1980 CKEMKIN(81)
Etats de l’art
1982 Spadaccini(137)
1984 Cheng(25)
1985 Cox(36)
1990
1991 Yetter(155)
Etudes du Diesel
1994 Baritaud(8) 1994 Spadaccini(136)
Modélisation 1995 Kong(87) Simulations numériques
du mélange 1996 Pfahl(116) directes
1997 Mastorakos 1997 Wan(151) 1997 Mastorakos(102)
(103)
1998 Im(74)
2000
La Chronologie bibliographique ci-dessus et celles qui vont suivre mettent en évidence des
regroupements thématiques non exhaustifs de travaux en relation directe (traits pleins noirs) ou
indirecte (traits blancs) avec le sujet.
32 Chapitre 2. Etat de l’art
F. 2.8 – Photographie d’un micro-propulseur / allumeur en phase de test. La structure en ’queue
de tigre’ est caractéristique d’une détente dans un jet supersonique. (Photo : Stork / Aerospace
Propulsion Products)
thèse, ce sont les allumeurs de moteurs fusées qui génèrent ce type de jet : ils produisent des gaz
chauds à une pression supérieure à celle de la chambre de combustion (Fig. 2.8). Lorsque les rap-
ports de pression statique sont de l’ordre de 10 ou plus, les fluides subissent de fortes détentes, des
chocs, et des couches de cisaillement intenses. Ce type particulier de jet est illustré en Fig. 2.9
Les premiers éléments de la théorie des jets supersoniques sont issus de la guerre mondiale
1939-1945. La thèse de I. S. Chang propose dès 1945 une prédiction de la forme des jets sous-
détendus par la méthode des caractéristiques. On peut également citer A. Ferri (58; 59; 60; 61),
précurseur des études supersoniques à la NACA puis la NASA.
Les dispositifs expérimentaux dédiés spécifiquement aux jets sous-détendus commencent en
1959. Les caractéristiques non visqueuses (longueur et rayon de la cellule de détente) sont assez
vite connues. A. D. Birch (14; 15)établit ensuite le lien entre le panache d’un jet sous-détendu et
celui d’un jet supersonique établi1 équivalent. Ces connaissances sont détaillées dans l’Annexe ’Di-
mensionnement de jets sous-détendus’. Les couches de mélanges supersoniques courbées de ces
jets sous-détendus font encore aujourd’hui l’objet d’observations poussées. On notera par exemple
1
Un jet supersonique établi est un jet sortant à la même pression que le milieu ambient.
2.3. La dynamique des jets supersoniques sous-détendus 33
F. 2.9 – Photo de type Schlieren d’un jet sous-détendu d’air à 270K, rapport de pression statique :
32.45, Diametre de buse : 30 mm : exposition 1/30s. Avec l’aimable autorisation de V.I. Zapryagaev,
ITAM Novossibirsk.
l’image de B. Yip (156) , une coupe laser, montrant une structure tourbillonnaire proche de l’insta-
bilité Kelvin-Helmotz, ou la série de mesures de V. Zapryagaev (160; 159; 161) sur l’instabilité de
Görtler, qui s’avère déterminante dans les couches de mélange des jets sous-détendus réactifs.
Les simulations numériques de ce phénomène ont commencé tôt avec la méthode des ca-
ractéristiques. Cette méthode est toujours utilisée aujourd’hui lorsque l’on s’intéresse seulement
à l’intérieur de la cellule de détente (112). Les méthodes aux volumes finis sont devenues appli-
cables à ce type d’écoulements à partir de 1985 (Code SCIPVIS (43)) et ont été améliorées depuis
avec le développement de modèles RANS (2). Actuellement, le code utilisé par la NASA se nomme
VULCAIN, issu d’une collaboration NASA-Taitech (152).
Les simulations aux grandes échelles sont difficiles à réaliser sur des configurations contenant
des chocs. C’est pourquoi elles se limitent pour l’essentiel à l’acoustique de jets supersoniques
adaptés, comme dans le travail de Bogey (17). A l’heure actuelle, les chocs sont majoritairement
traités en LES par des schémas de type WENO, sur des maillages structurés. Ils ont été récemment
appliqués avec succès en deux dimensions sur un jet sous-détendus par T.S. Cheng (26).
34 Chapitre 2. Etat de l’art
Théorie supersonique
1945 Chang(21)
1970
1974 Antsupov(5)
Simulations numériques
1984 Birch(14) de jets supersoniques
1985 SCIPVIS(43)
1986 Yoshisawa(158)
1987 Birch(15)
1988 Zapryagaev(160)
Etudes liees aux 1988 Birch(13)
jets supersoniques
1990 1989 Yip(156) 1989 Abdol-Hamid(1; 2)
1991 Sandham(129)
Simulations numériques
de l’acoustique des
jets supersoniques
1993 Arnette(6)
1994 Lele(96) 1994 Mankbadi(101)
1995 Clemens(30) 1995 Cumber(38) 1995 Tam(142)
1997 Papamoschou(115)
1997 Lee(95)
1998 Palmer(112)
1998 Day(44) 1999 VULCAN(152) 1998 Shen(133)
1999 Panda(114) 1999 CE/SE(163)
2000
2000 Freund(63)
2001 Zapryagaev(159)
2002 Kourta(88)
2003 Bogey(17)
2004 Zapryagaev(161)
2005 Cheng(26)
Plusieurs revues complètes de la combustion supersonique sont disponibles : Waltrup (150), Ja-
chimovski (75), Billig (12), Spadaccini et Colket (136; 33) et enfin Curran (39; 40). On y apprend
essentiellement que le problème de combustion des super-stato-réacteurs est d’obtenir une com-
bustion suffisamment rapide pour s’effectuer à l’intérieur du moteur malgré des temps convectifs
extrêmement courts. Le dimensionnement du moteur et de la zone de réaction deviennent des pa-
ramètres critiques. C’est pourquoi de nombreuses études se focalisent sur le mélange, l’injection de
carburant, et la cinétique chimique. Concernant la simulation numérique, les avis sur le choix des
modèles de turbulence ou de combustion sont moins tranchés car l’augmentation directe de la puis-
sance de calcul (c.a.d. des maillages plus fins ou un plus grand nombre de particules stochastiques)
reste encore aujourd’hui le meilleur moyen d’améliorer les résultats. Comme pour les jets super-
soniques, la simulation des grandes échelles a essentiellement été réservée à des cas académiques
de flammes supersoniques bidimensionelles. Le modèle de combustion le plus adapté passe par des
fonctions de densité de probabilité, couplées à des approches RANS ou SGE.
2.4. La combustion supersonique 37
Etudes liées à la
combustion supersonique 1964 Ferri (60)
1970
1972 Suttrop (141)
1973 Burrows (19)
Flammes
supersoniques
1980
1980 Evans (55)
Simulations de
combustion supersonique Etats de l’art
1986 Drummond (50)
1986 Kumar (89) 1987 Waltrup (150)
1988 Birch (13) 1988 Jachimowski (75)
1990
Equations
Les équations résolues par l’approche SGE sont résumées dans ce chapitre. Le code utilisé
dans cette thèse pour résoudre ces équations se nomme AVBP. Il est développé conjointement par
l’Institut Français du Pétrole et le CERFACS.
Les équations de la combustion et de la dynamique des fluides, communes à toute approche
SGE, sont décrites en début de chapitre. Elles sont suivies d’une sélection des modèles existants
pour la dynamique et de la combustion, puis de l’implantation numérique spécifique utilisée dans le
code AVBP.
∂w
+∇·F = s (3.1)
∂t
où w = ( ρu, ρv, ρw, ρE, ρk )T est le vecteur des variables conservatives avec respectivement ρ, u, v,
w, E, ρk la densité, les trois composantes du vecteur vitesse V ~ = (u, v, w)T , l’énergie par unité de
masse et ρk = ρYk avec Yk la fraction massique de l’espèce k, k variant de 1 à N, N étant le nombre
d’espèces. Une décomposition classique en parties visqueuse et non visqueuse est utilisée pour le
tenseur des flux :
F = F(w)I + F(w, ∇w)V (3.2)
Les trois composantes du tenseur non visqueux F(w)I sont :
Termes non visqueux :
39
40 Chapitre 3. Equations
où la pression hydrostatique P est déterminée par l’équation d’état des gaz parfaits (Eq. 3.15).
Termes visqueux : Les trois composantes du tenseur visqueux F(w, ∇w)V sont :
−τ xx
−τ xy
f V = −τ xz
−(uτ xx + vτ xy + wτ xz ) + q x
J x,k
−τ xy
−τyy
g =
V
−τyz
(3.4)
−(uτ + vτ yy + wτyz ) + qy
xy
Jy,k
−τ xz
−τyz
hV = −τzz
−(uτ + vτ yz + wτzz ) + qz
xz
Jz,k
avec µ la viscosité dynamique explicitée en section 3.1.5. Ji,k est le flux diffusif de l’espèce k suivant
i défini en section 3.1.4 (Eq. 3.26) et qi est le flux de chaleur défini en section 3.1.6 (Eq. 3.32).
Les capacités calorifiques à pression constante C p,k et volume constant Cv,k sont supposées constantes
entre T i et T i+1 = T i + 100. Elles sont respectivement déterminées par les Eq. 3.10 et 3.11. L’énergie
sensible est définie par une interpolation linéaire via la température (Eq. 3.12). L’énergie sensible et
l’enthalpie sensible du mélange sont respectivement définies par les Eq. 3.13 et 3.14.
∂h s,k
C p,k = (3.10)
∂T
∂e s,k
Cv,k = (3.11)
∂T
ẽ s,k (T i+1 ) − ẽ s,k (T i )
e s,k (T ) = ẽ s,k (T i ) + (T − T i ) (3.12)
T i+1 − T i
N
X N
X
ρe s = ρk e s,k = ρ Yk e s,k (3.13)
k=1 k=1
XN XN
ρh s = ρk h s,k = ρ Yk h s,k (3.14)
k=1 k=1
R
P=ρ T (3.15)
W
N
1 X Yk
= (3.16)
W k=1
Wk
où W est la masse molaire du mélange. La constante du mélange variable en temps et en espace r et
les capacités calorifiques dépendent des fractions massiques :
N N
R X Yk X
r= = R= Y k rk (3.17)
W k=1
Wk k=1
N
X
Cp = Yk C p,k (3.18)
k=1
XN
Cv = Yk Cv,k (3.19)
k=1
où R = 8.3143 J/mol.K est la constante universelle des gaz parfaits. L’exposant polytropique du
mélange est donné par γ = C p /Cv . Donc, la constante du mélange, les capacités calorifiques et
l’exposant polytropique dépendent de la composition local du mélange définie par les fractions
massiques locales Yk (x, t) :
A partir de l’énergie sensible, on déduit la température en utilisant les Eq. 3.12 et 3.13. Enfin, les
conditions limites font intervenir la vitesse du son du mélange c définie par l’Eq. 3.21.
c2 = γrT (3.21)
µC p
λ= (3.28)
Pr
1
Les erreurs liées à cette hypothèse sont faibles par rapport aux propriétés thermodynamiques
3.1. Equations d’écoulements réactifs 43
La diffusivité moléculaire Dk est exprimée à l’aide des coefficients binaires Di j obtenus à l’aide de
la théorie cinétique des gaz (72). La diffusivité moléculaire Dk est définie par l’Eq. 3.29 (16).
1 − Yk
Dk = PN (3.29)
j,k X j /D jk
Les coefficients binaires Di j sont des fonctions complexes dépendant des intégrales de collision et
des variables thermodynamiques. Dans un code de simulation numérique directe utilisant une chi-
mie complexe, utiliser l’Eq. 3.29 a un sens. Cependant, dans la plupart des codes de simulation
numérique directe, un schéma simplifié est utilisé et la modélisation de la diffusivité n’a pas besoin
d’être aussi précise. En conséquence, une approche plus simplifiée est utilisée. En faisant l’hy-
pothèse que les nombres de Schmidt de chaque espèce S ck sont constants, la diffusivité moléculaire
Dk est définie par l’Eq. 3.30.
µ
Dk = (3.30)
ρ S ck
où S c,k est le nombre de Schmidt de l’espèce k supposé constant en temps et en espace. Pr et S ck
modélisent la diffusion laminaire thermique et moléculaire. Les valeurs sont obtenues à l’aide du
logiciel PREMIX, inclus dans le package CHEMKIN (81) en calculant leur valeur dans les gaz
brûlés pour une flamme de prémélange laminaire monodimensionnelle.
E(k)
grandes échelles résolues petites échelles modélisées
−5/3
k
fréquence
de coupure
Grâce à l’utilisation du filtre pour séparer petites et grandes échelles, la SGE permet une représen-
tation dynamique des tourbillons de plus grande taille dont les contributions jouent le premier rôle
en géométrie complexe. Ainsi la prédiction des écoulements turbulents est améliorée puisque les
phéno-mènes à grande échelle, comme la propagation des ondes acoustiques, sont intrinsèquement
présents dans les équations de conservation (119).
La SGE a donc un potentiel évident pour la prédiction des écoulements turbulents réactifs ren-
contrés dans les applications industrielles. Cependant, cette utilisation est limitée par les hypothèses
introduites pour construire les modèles de la SGE.
Les équations de conservation SGE sont décrites en section 3.2.2 et les modèles de sous-maille
en section 3.2.3.
Wk ∂Xk
!
en non réactif : Ji,k = −ρ Dk − Yk Vi c
(3.55)
W ∂xi
Wk ∂X
ek c
approximation : Ji,k ≈ −ρ Dk − Yk Vi , i = 1, 2, 3 (3.56)
e e
W ∂xi
le flux de chaleur qi
N
∂T X
qi = −λ + Ji,k h s,k (3.57)
∂xi k=1
N
∂Te X
approximation : qi ≈ −λ + h s,k ,
Ji,k e i = 1, 2, 3 (3.58)
∂xi k=1
48 Chapitre 3. Equations
Les variations spatiales des flux de diffusion moléculaire sont négligeables et un modèle de gradient
suffit pour les modéliser. Les équations de conservation font apparaı̂tre des termes de sous-maille
dont les modélisations sont présentées dans la section 3.2.2.
Termes de sous-maille
Les trois composantes du tenseur de sous-maille sont définies par l’Eq. 3.59.
t t t
−τ xx −τ xy −τ xz
−τ t −τ t −τ t
xy yy yz
t
t
t t t
f = −τ xz , gt = −τyz , h = −τzz
(3.59)
q t q t q t
x t y t z t
J x,k Jy,k Jz,k
Les différents termes s’écrivent :
le tenseur des contraintes τi j t
τi j t = −ρ(ugiu j − e uie
u j ), i, j = 1, 3 (3.60)
1
modèle : τi j t = 2 ρ νt (Sei j − δi j Sell ), i, j = 1, 3 (3.61)
3
La modélisation de νt est détaillée en section 3.2.3.
t
le tenseur de diffusion des espèces Ji,k
t
Ji,k = ρ (ug i Yk − e
ui Y
ek ), i = 1, 2, 3 (3.62)
t Wk ∂X
t ek c,t
modèle : Ji,k = −ρ Dk ei ,
− Yek V i = 1, 2, 3 (3.63)
W ∂x i
νt
avec : Dtk = (3.64)
S ctk
Le nombre de Schmidt turbulent S ctk est un paramètre délicat à déterminer. Des études ont
montré qu’il n’était pas constant, et des modèles plus complexes ont été développés. Cepen-
dant, dans une grande majorité d’études, S ctk est pris constant, proche de l’unité.
le flux de chaleur qi t
qi t = ρ(ug
iE − e
ui E),
e i = 1, 2, 3 (3.65)
N
∂Te X t e
modèle : qi t = −λt + Ji,k h s,k , i = 1, 2, 3 (3.66)
∂xi k=1
µt C p
avec : λt = et Prt = 0.9 (3.67)
Prt
Les vitesses de correction sont obtenues à partir de l’Eq. 3.68.
N
µ µt Wk ∂X
X !
c,t
ek
Vi + Vi =
e c e + , i = 1, 2, 3 (3.68)
k=1
ρS ck ρS ctk W ∂xi
3.2. Equations SGE d’écoulements réactifs 49
Le modèle de Smagorinsky
q
νt = (CS 4) 2 Sei j Sei j ,
2
i, j = 1, 3 (3.72)
avec : CS = (0.1 − 0.18) suivant l’écoulement (3.73)
Le modèle de Smagorinsky (135), développé dans les années 1960, fait l’objet de très nombreux
tests dans la littérature sur de multiples types d’écoulement. Il est très utilisé du fait de sa simplicité.
Cette fermeture a la particularité, dans le cas d’une THI, de fournir le bon niveau de dissipation de
l’énergie cinétique. Ce modèle est connu pour être trop dissipatif, spécialement près des murs et son
utilisation pour des régimes de transition vers la turbulence n’est pas recommandée (128).
Le modèle WALE
1 2 1 2
sdij = gi j + e
(e g2ji ) − e g δi j , i, j = 1, 3 (3.74)
2 3 kk
(sdij sdij )3/2
νt = (Cw 4) 2
, i, j = 1, 3 (3.75)
(Sei j Sei j )5/2 +(sdij sdij )5/4
avec : Cw = 0.4929 (3.76)
50 Chapitre 3. Equations
q
νt = (CS F 4)2 2 HP(Sei j ) HP(Sei j ), i, j = 1, 3 (3.77)
avec : CS F = 0.37 (3.78)
Dans l’Eq. 3.77, HP(Sei j ) est le tenseur des déformations résolu obtenu à partir d’un champ de vi-
tesse filtré au travers d’un filtre passe-haut. Ce modèle fut développé afin d’améliorer la représentation
de phénomènes locaux typiques des écoulements turbulents (51). Avec le modèle de Smagorinsky
filtré, la transition vers la turbulence est mieux prédite.
entrées/sorties parallèles de données et les méthodes itératives. AVBP est portable sur les principales
plateformes de calculs (PCs, stations de travail, super-calulateurs), restant efficace sur la majorité
des architectures parallèles.
où w est le vecteur des variables conservatives et F~ est le tenseur de flux correspondant. Ce tenseur
est décomposé en une partie visqueuse et une partie non visqueuse au travers de l’équation 3.80.
Les termes spatiaux sont approchés dans chaque volume de contrôle pour en déduire le résidu RΩ j
défini par l’équation 3.81.
où ∂Ω j est la projetée de Ω j le long de la normale ~n. L’approche Cell-Vertex est applicable à des
types de cellules quelconques et apparaı̂t donc très utile pour des maillages hybride. Le résidu RΩ j
défini par l’équation 3.81 est d’abord calculé pour chaque cellule en utilisant une loi d’intégration
simple appliquée aux faces. Pour des faces triangulaires, les composantes individuelles du flux sont
supposées varier linéairement. Pour des faces quadrilatères, où les nœuds peuvent ne pas être co-
planaires, afin d’assurer l’exactitude de l’intégration pour des éléments quelconques, chaque face
est divisée en triangles et intégrée sur chaque triangle. La valeur du flux est obtenue après moyen-
nage des quatre triangles ( deux divisions le long des deux diagonales). Cette propriété qui permet
de conserver la linéarité joue un rôle important car elle assure une bonne précision même sur des
52 Chapitre 3. Equations
maillages irréguliers. Il est utile d’écrire le résidu RΩ j défini par l’équation 3.81 au travers d’une
cellule quelconque. L’équation 3.81 devient alors l’équation 3.82.
1 X ~ ~
RΩ j = Fi · dS i , (3.82)
Nd VΩ j i∈Ω
j
où F~i est une approximation de F~ aux nœuds, Nd représente le nombre de dimensions de l’espace
et {i ∈ Ω j } sont les sommets de la cellule. Dans cette formulation, l’information géométrique a été
factorisée en termes dS~ i qui sont associés aux nœuds individuels de la cellule et non à ses faces. dS ~ i
est la moyenne des normales pondérées par les surfaces des faces triangulaires d’un nœud commun
i, i ∈ Ω j . Pour assurer la consistance, i∈Ω j dS
P ~ i = ~0. La linéarité de l’opérateur de divergence est
préservée si le volume VΩ j est défini par l’équation 3.83.
1 X ~ i,
VΩ j = ~xi · dS avec ∇ · ~x = Nd (3.83)
Nd2 i∈Ω
j
Une fois les résidus calculés, le schéma semi-discret est défini par l’équation 3.84.
dwk 1 X k
=− D VΩ RΩ , (3.84)
dt Vk j|k∈Ω Ω j j j
j
où DkΩ j est la matrice de distribution qui fait une projection pondérée du centre Ω j vers le nœud
k (scatter) et Vk est le volume de contrôle associé à chaque nœud k. La conservation est garantie
si k∈Ω j DkΩ j = I. Dans le contexte présent, l’équation 3.84) est résolue, en utilisant une méthode
P
explicite Euler ou une méthode Runge-Kutta à plusieurs étapes en temps, afin d’obtenir une solution
stationnaire.
La matrice de distribution est définie par l’équation 3.85.
1 δtΩ j
DkΩ j = (I + C A ~ k ),
~ Ω j · dS (3.85)
nn VΩ j
Le nombre de nœuds de Ω j est nn et A ~ est la matrice jacobienne du tenseur des flux. Le schéma
classique centré, obtenu en choisissant C = 0, est stable combiné avec des pas de temps Runge-
Kutta. Un schéma de type Lax-Wendroff est obtenu en choisissant la constante C dépendante du
nombre de dimensions du problème et du type de cellule. Une forme simple de C est définie par :
C = n2n /2Nd .
l’équation 3.87. Une approximation de ce gradient au nœud, définie par l’équation 3.88, est obtenue
en utilisant une moyenne pondérée par le volume de la cellule.
∂w
! Z Z
1
≈ w · ~n∂S (3.86)
∂x C VC ∂ΩC
~
1 X ~
∇w = wi dS i (3.87)
Ωj VΩ j i∈Ω
j
~
1 X ~ Ω
∇w = VΩ j (∇w) j (3.88)
k Vk j|k∈Ω
j
Senseurs
Un senseur ζΩ j est un paramètre compris entre 0 et 1 défini pour chaque cellule Ω j . Dans le
cas où la solution est bien résolue, le senseur est égal à 0 alors que dans le cas où la solution a de
fortes variations locales, le senseur est égal à 1 et la VA est appliquée. Ce senseur est obtenu en
comparant différentes évaluations du gradient d’un scalaire comme la pression, l’énergie totale ou
la fraction massique. Si les évaluations sont identiques, le senseur est fixé à 0. Par contre, si les
deux évaluations donnent des valeurs différentes, le senseur est déclenché. Le point crucial est de
trouver une fonction pour le senseur qui soit différente de zéro seulement dans les zones utiles. Deux
fonctions de senseur sont utilisées dans ces travaux : le senseur de Jameson ζΩJ j (76) et le senseur de
Colin ζΩC (31) qui est dérivé du senseur de Jameson.
j
Notation l’indice k désigne les variables liées à un sommet k de la cellule considérée et l’indice Ω j
désigne les variables liées à la cellule Ω j .
54 Chapitre 3. Equations
Senseur de Jameson
Le senseur de Jameson ζΩJ j lié à la cellule Ω j (défini par l’Eq. 3.89) est le maximum de tous les
senseurs ζkJ liés aux sommets k (définis par l’Eq. 3.90). S est le scalaire évalué par le senseur et (∆k1 ,
∆k2 ) sont des évaluations différentes du gradient définies par l’Eq. 3.91. ∆k1 mesure la variation de S
au sein de la cellule Ω j . ∆k2 est une estimation de la même grandeur en utilisant (∇S ~ )k , le gradient
de S au nœud k. Ce senseur varie proportionnellement à l’amplitude de la déviation par rapport à
l’évolution linéaire. Ce senseur a une évolution douce d’un point de vue numérique et s’applique
parfaitement pour des cas quasi-stationnaires.
|∆k1 − ∆k2 |
ζkJ = (3.90)
|∆k1 | + |∆k2 | + |S k |
∆k1 = S Ω j − S k ~ )k .(~xΩ − ~xk )
∆k2 = (∇S (3.91)
j
Senseur de Colin
– ζΩ
C est très petit lorsque ∆k et ∆k sont petits comparés à S . Ceci correspond à des erreurs
j 1 2 Ωj
numériques de faible amplitude (si ∆k1 et ∆k2 sont de signes opposés) ou à des faibles gradients
bien résolus par le schéma (si ∆k1 et ∆k2 sont de même signe).
– ζΩ
C est petit lorsque ∆k et ∆k sont du même signe et du même ordre de grandeur (même si
j 1 2
cet ordre de grandeur est grand). Ceci correspond à des gradients raides bien résolus par le
schéma.
– ζΩ
C est grand lorsque ∆k et ∆k sont de signes opposés et qu’un des deux est beaucoup plus
j 1 2
grand que l’autre. Ceci correspond à une oscillation numérique de grande amplitude.
– ζΩ
C est grand si ∆k ou ∆k est du même ordre de grandeur que S . Ceci correspond à une
j 1 2 Ωj
situation non physique résultant d’un problème numérique.
3.4. Les conditions limites 55
Il est à noter que les définitions de Ψ et k changent pour l’équation de conservation des fractions
massiques : la valeur de référence n’est plus S k mais 1, valeur maximum de la fraction massique.
Ψ − Ψ0
!! !!
1 1 −Ψ0
ζΩ
C
= 1 + tanh − 1 + tanh (3.92)
j 2 δ 2 δ
∆ k
!
avec : Ψ = max 0, k ζJ (3.93)
k∈Ω j |∆ | + 1 S k k
∆k = |∆k1 − ∆k2 | − k max |∆k1 |, |∆k2 | (3.94)
max |∆k1 |, |∆k2 |
k = 2 1 − 3 k (3.95)
|∆1 | + |∆k2 | + S k
Ψ0 = 2.10−2 δ = 1.10−2 1 = 1.10−2 2 = 0.95 3 = 0.5 (3.96)
Opérateurs
Les modèles de viscosité artificielle utilisent deux opérateurs dont les propriétés sont :
2ème ordre cet opérateur agit comme une viscosité classique. Il adoucit les gradients et introduit de
la dissipation artificielle. Il est associé à un senseur. Ainsi, le schéma numérique garde son
ordre de précision dans les zones à faible gradient et assure la stabilité et la robustesse du
schéma dans les zones critiques. A l’origine, il était utilisé pour les chocs mais peut en fait
adoucir n’importe quel gradient trop fort.
4ème ordre cet opérateur est utilisé pour diminuer les wiggles.
Les contributions à la cellule de l’opérateur du 2ème ordre (Eq. 3.98) et de l’opérateur du 4ème ordre
(Eq. 3.99) sont reportées sur les nœuds de cette cellule Ω j (Eq. 3.97).
X X
dwk = R2k∈Ω j + R4k∈Ω j (3.97)
j j
1 VΩ j
avec : R2k∈Ω j = − smu2 ζΩ j (wΩ j − wk ) (3.98)
Nv ∆tΩ j
1 VΩ j h
~ Ω · (~xΩ − ~xk ) − (wΩ − wk )
i
avec : R4k∈Ω j = smu4 (∇w) (3.99)
Nv ∆tΩ j j j j
Notation smu2 et smu4 sont des coefficients sans dimension fixés par l’utilisateur.
Connaissant la solution wn au temps t, la solution wn+1 au temps t + ∆t s’écrit pour chaque noeud i :
∆t
wn+1
i = wni − .dwni (3.100)
Vi
avec wni = w(t, ~xi ) et wn+1
i = w(t + ∆t, ~xi ), ~xi étant le vecteur des coordonnées, Vi le volume autour
n
du noeud i et dwi le résidu au noeud i calculé par le schéma numérique.
L’exposant n est utilisé pour rappeler que le code explicite ne se sert que des données au temps
n. L’Eq. 3.100 est utilisée à chaque noeud du domaine Ω. Afin d’imposer une valeur sur les limites
du domaine de calcul ∂Ω, l’Eq. 3.101 est définie.
En chaque noeud du bord, le résidu calculé par le schéma est remplacé par un résidu imposé par la
condition limite. En plus de cette méthode, l’utilisation d’un coefficient de relaxation peut être uti-
lisé pour tendre vers la valeur cible de manière plus ou moins douce. Ensuite, les conditions limites
sont classés en conditions non caractéristiques qui modifient directement les variables conserva-
tives et en conditions caractéristiques qui effectuent une décomposition en ondes pour modifier les
résidus (138). La méthode NSCBC développée par Poinsot et Lele (118) est utilisé dans ce second
cas.
Chapitre 4
Auto-allumage
Z Kg de carburant
1−Z Kg d’oxydant
Pression constante
Enthalphie Constante
57
58 Chapitre 4. Auto-allumage
Dans cette configuration, seule la cinétique chimique est à l’œuvre. Les mélanges oxydant/réducteur
sont étudiés au moyen d’une quantité scalaire normalisée, la fraction de mélange Z. Elle identifie la
composition du mélange homogène et non réactif de deux fluides. Un kilogramme de mélange à la
fraction Z correspond à ZKg de fluide réducteur additionné de 1 − ZKg de fluide oxydant. Ce sca-
laire est basé sur un bilan de masse des atomes mis en jeu, qu’ils soient dans les molécules initiales
(CH4 , O2 , H2 , ...) ou les molécules finales (H2 O, CO2 , ...). Ainsi, la quantité Z est indépendante de
la réaction. Par exemple, la fraction de mélange dans une flamme prémélangée est une constante.
La température initiale T (Z) se calcule avec l’équilibre adiabatique de deux fluides non réactifs.
Cette configuration est notée HMI (Homogeneous Mixing Ignition) dans la suite de cette thèse. Le
prémélange homogène, placé dans un réacteur adiabatique isobare (enthalpie et pression constante),
va réagir en un temps tHMI dépendant du mélange initial Z. Un tel dispositif se calcule rapide-
ment au moyen de logiciels de cinétique chimique courants, tels les modules PSR ou AURORA
de CHEMKIN. Dans cette configuration, on teste l’ensemble des mélanges possibles, et on évalue
Temperature
Time
THMI
0
ve etry
acti iom
1 t Re ech
Mos Sto
ZMR
ZST
Mixture Fraction
0
F. 4.2 – Evolution qualitative de la température dans des réacteurs remplis initialement de
différents mélanges auto-allumants.
leurs temps d’auto-allumage. Le mélange le plus rapide à l’auto-allumage est identifié par sa frac-
tion de mélange Z MR (Most Reactive). Le mélange stœchiométrique ZS T reste le plus exothermique,
sans être le plus rapide, voit Fig. 4.2. Les mélanges trop pauvres ou trop riches ne donnent que des
réactions faibles et lentes. Généralement, le mélange Z MR se situe du coté du réactif le plus chaud.
Ainsi, pour un oxydant chaud et un réducteur froid 0 ≈ Z MR < ZS T . Pour un réducteur chaud et
un oxydant froid, ZS T < Z MR ≈ 1. Les résultats obtenus pour les couples hydrogène/oxygène et
méthane/oxygène sont présentés dans l’article ”Comparison of computational methodologies for
ignition of diffusion mixing layers” (82) intégralement reproduit ci-après. On retrouve bien dans
la Fig. 4 de cette publication l’évolution du temps d’allumage décrite précédemment. Les calculs
faisant intervenir le méthane sont en accord avec les mesures expérimentales, validant la procédure
de calcul. Le comportement des deux carburants est qualitativement le même. On pourra cependant
4.1. Configurations de références 59
noter que la fraction de mélange la plus réactive du méthane est extrêmement pauvre (Z MR = 0.002)
ce qui induit une augmentation de température de seulement 25 degrés. La dénomination ’mélange
le plus réactif’ devrait donc se limiter à ’mélange le plus rapide’. La chaleur dégagée est trop faible
pour parler d’allumage. Ceci mène à l’étape suivante : comment les processus diffusifs propagent-ils
la faible combustion des premiers sites d’allumage à l’ensemble de la flamme de diffusion ?
Pression constante
Enthalphie Constante
T T T T
Z ZMR Z Z ZST Z
1
Fraction de melange
Flamme riche
Flamme de diffusion
Z
ST
Z
MR
Flamme pauvre
0
0 t Temps
LMI
nombre Lewis contrôle le temps d’allumage (cf. Fig.9a et 9b de l’article). Pour résumer, tLMI peut
majorer d’un ordre de grandeur tHMI . Un mélange HMI avec un Z MR suffisamment exothermique
(cas de l’hydrogène) contrôle toute la phase d’induction (Fig. 10a). Un Lewis égal à un implique
une diffusion qui retarde le temps d’allumage jusqu’a tLMI ≈ 7tHMI en Fig.9a. Les Lewis faibles
diminuent substantiellement le temps d’allumage. Dans un mélange HMI, un Z MR trop pauvre (cas
du méthane) ne contrôle plus l’induction, et augmente tLMI . Cette majoration ne peut être compensée
par un Lewis faible, comme l’indique la Fig. 9b.
4.1.3 Turbulence
L’impact de la turbulence sur l’auto-allumage nécessite des outils de simulation plus complexes,
tels que la simulation numérique directe (102; 74; 69; 53). Cette approche résoud les équations de
Navier-Stokes avec une précision suffisante pour se passer de toute modélisation. La puissance de
calcul nécessaire limite généralement ces simulations à des nombres de Reynolds faibles (< 200)
pour des milieux bidimensionnels. Les travaux sur l’auto-allumage (8; 86; 116; 151) s’intéressent à
l’impact d’une turbulence homogène isotrope sur une couche de mélange réactive. Sans entrer dans
les détails, on peut en retenir quelques points essentiels, illustrés par la fig. 4.5 :
1. Les premiers sites d’allumage se situent sur une iso-richesse correspondant à celle du pré-
mélange parfait le plus rapide à s’allumer Z MR . Sur cette iso-richesse, ils débutent là ou la
dissipation scalaire χ est la plus basse. Elle mesure un gradient d’espèce normalisé dans un
62 Chapitre 4. Auto-allumage
1 ∂Z ∂Z
" #
χ= (4.5)
S c ∂xi ∂xi
2. Le temps d’apparition des premiers sites d’auto-allumage tT RB1 est peu sensible à la turbu-
lence.
3. La turbulence réduit le temps d’auto-allumage de l’ensemble de la couche de mélange tT RB2 .
ZMR ZST
Dans notre recherche d’un estimateur du temps d’auto-allumage, la turbulence simplifie les choses
en faisant tendre le nombre Lewis global vers l’unité. En effet, le mélange par agitation turbulente
agit sur la chaleur comme sur les espèces. Lorsque que ce mélange est suffisamment intense, il
supplante la diffusion molóculaire et la diffusion thermique. Les effets de diffusion différentielle
sont alors réduits. Au final, l’auto-allumage d’un couple de réactifs séloigne de la configuration
HMI lorsque la diffusion différentielle est forte , c.a.d. le nombre de Lewis est loin de l’unité.
L’écart le plus net est obtenu avec la configuration LMI. Enfin, la turbulence diminue globalement
cet effet.
4.2 Simulations
L’auto-allumage intervient dans les brûleurs séquentiels de certaines turbines à gaz. C’est pour-
quoi l’étude présente à fait l’objet d’une publication conjointe parue dans le numéro 175 du journal
Combustion Science and Technology : Le Dr. R. Knikker s’est intéressé au couple méthane/oxygène
pour les applications de turbine à gaz, tandis que l’auteur s’est focalisé sur le couple hydrogène/oxy-
gène à destination des moteurs fusées. L’article reproduit ci-après présente les calculs HMI et LMI
effectués pour les deux couples, et interprète les différences de comportement. Il met en évidence
une forte dépendance de l’auto-allumage devant la diffusion différentielle. Par ailleurs, les infor-
mations obtenues vont permettre de créer une cinétique de combustion réduite adaptée à l’auto-
allumage (Chap.6).
4.2. Simulations 63
Pression constante
Enthalphie Constante
4.2. Simulations 65
) ) & /( =, . 0
) )
(- "
'
$
0
>0 / (-
)
& #&
) &
/
)
&- "
/ (
& +%
?63? ,>+ - 2@@.
& , . )
,
+
/( # ) . /
) )
- "
( #
/ & /
)
& # ) - 3
)
& 9
/
(
&- 0 & , .-
"
:13
/(
)
- :0
$ &
#& ( ) & - "
) %
' # 0&
#
/(
&- 0
)
( ) &
$
'
#&
'
&
(
& / --
& ) & &
/&
40-
"
( 0(
)
#
:13- "
>13
) #& &- 3
#&
/ ( ) /
#& /
%
' ( ,&- .- "
)
#
+
) #
- "
0
4
3
,A
- .-
"
>13 )
/
) %
( / 0
Æ 9 0
#
- 1 # & /
) # B? ) /
#& & ,:/
"
7 3 - 22@7 19 - 22!.- A
>13 ) / ) (
/
/$ #& (
%
'
&/ & /
-
"
:13
>13 )
(
&
*
"
)
:13
&- -
& #
/
/(
0
) ) (
/
(
&
#(& 7
& 0&
/%
)
(
)
-
0&
)
/ & 0
- "
#
+
) C 6
) D )
4.2. Simulations 67
#(&
!2 D ) & , C
C
!2.- "
)
(
&
& )
C 6- "
(
&
&
) (
C
C
! - "
, C .- "
# ) )
C
E- ( :/
"
,.
9
) F - ,22.
)
(
& - "
/
)(
&
, ) # 3 -
22@7 6+
>0 22E.- "
)
:13
0 &- -
"
&
(
0
& ) # )
)
# ) C
2 )
#
)- "
& &
(
C , .
C
- ?
) # )
&
(
( &( ) )
# ) / -
&
)
# )
&
(
0 #
/
$
&
&
(
(
)
# 0
( 0
- "
)
/
&
#, .- G&
%
6 0
) &
) )
- 3 &-
&
)
# )
(
#( & &
(-
"
)
4
& &
(
)
) &
)
(
& 4- "
#
+ &
&
C 6
&
(
- "
) # )
&
) , C
C
.
)
C 6- "
C -
# )
C
-
:& &
(
&
'
&- /- "
3 -
)- "
' 0
-
(
&
0
*
)
)
/-
/
0
%
/( & -
,22@. 0
/
)
#
)
3 -
-
3
&
0
'
)
/
)
0
& *
( '-
/0- (
&
(
# ,(/
&.
/(
6
9 ,2@. 0
$
C
# ,E 2 . H I H
I ,.
0
&
(
HI
H
I
,J. ) #(&
- ?
)
) &
)
(
& 4- "
&
68 Chapitre 4. Auto-allumage
&
( )
#
(
- 1
'
)
(
/ )
)
(
:
(
(
/
)
)
( /
%
&
%/
& ,8/9 .- ?
#
# ) )
- ,.-
?
&
( /
)
&- /-
8
&
%
# ) )
C
&
( ) C
2 - :0
& &
# ) (
E 6- "
(
- ?
0
# )
& &
(
) )
0
(
&
- #
)
#
-
(
;&<
:13 # / )
& - B)
' 9 -
4.2. Simulations 69
A
(
&
?%9
/
(
,? A . ) ,L
> 22. 0
/ #
% / ,A - 22 .- "
) /
( ( %
4& 0
#
&
/
(
) )
- ? #
( /
)& /
-
"
M
/(
/(
/( C
MM
E )
C MM
)
'
0
M
/0 0 &
9%
(- 3 0 90
> / )
#
( /
- "
)
/
)
40
&
( ' ( ) - "
#
( (
/ Æ
( / / - "
/
0
&
#
3 -
-
)
& )
- 5
%
0
&
/( )&
)
)
& /-
" (
) /0 )
&
( ( 9
)
&
- "
%
+
0
# )
/($
C = ), . ,.
0
)
0
)
%)
%
- "
,
(
).
&
0 /
)
# ) )0&
)
& - - ,.
/
( )
' /
)
) ( >0 / ,
'.
' Æ- 3
9
&
'
- " /
0 /
&
( / -- (
)
/
/-
"0
' &
0
&
&
)
:13 - "
)0&
)
>13 0 )
)0& $
4 )
' ) %( >0
/
/0
&
"
' )
>13 )
(
&
&- - A0 C
C
*
E
70 Chapitre 4. Auto-allumage
' 0
)
/
-
"
C
(
)
,-- ) ) .- "
&
/& &
4 & (
( 0
& /
& 0
- "
& )
4
0 &- / (
(
&
4 &
/
' 0 /
0
)0& (-
/
)
&
( )
%
&%
/ &-
/* / )
:/
:0
/(
&
4 0
& (
& 0
&- "
&
0
&
& - ,E. 0 /
&
& ( )
4 )-
0
&& &
(
&
( & /($
N
C ,!.
)
/
% 0
&&/
)0
/( # )
&
- "
0 (/
&
&
(
& - ,E.
,!.$
,.
#
)
,- !.
&
( /)
,.
&
# )
,- E.- "
0 (
)
&
(
& - ,!. / )
) &-
"
4 )
9 )
#& (
&
(
0 &- !- )
C &
# )
!
4.2. Simulations 71
.
&
& &
& , M
0
.-
0 /( &- ! ) /
(
& 4 ,.
4
,.
&
(
0
)
& - )%
4 )
&&/
) (
)
& /
- %
&
(
/( /
- ,E.
,!.
- #
) #
,- E. &
( )
&
(
)0&
,!. / /
0
) ) -
" &
4
&
( )
( >0 / %
)
)
(
&
&
+
& %( >0 / '- "
(
/
/( %) 0
0
/
0
) (
(
&
#(& ) 0 & /( C
C - "
& )
)
) ) / (
&
- 3
)0& %
(
&
#(& >0 / C
C
& &
(
)
)
(
- /
0
+
&- !- 3
/
)
&
( ) ) #(
)
( >0 / +- 5)
& (
& /0
)
(
/
&
)
) Æ(
) -
"
(
& ( &
&
:/
"
,.- 3
09 >13
&
>0 / )
)
&
(- "
(
/
&- 2 )
(
& )
,
&
. )
&- 2/- 3 &
&
,:/
"
7 3 - 22@7 19 - 22!.
& / 40
)
#& (
# 0
)
-
' )
(
& 0
&
& 0
/
)
#(&
)
%( >0 / ' )
#
&& 4-
&
( /
0
( >0 / )
& /(
%
(-
(
& 4
/ >0 /
&
( / C
0
)
(
& >0 /-
(
>0 / ) C - "
/ &&
&
'( )
(
& ) & / )
)
&
/(
/&
(
/(
#& - "
& 0
& &
(
)
>0 /- 5
>0 / )
)
- "
)
' /0
/
0
( >0 / &- 2/ - 3
( )
&
/( #&
& &
(
&
& &
(-
"
& / ) / ( 0
%&
(
4 /+
-
2
4.2. Simulations 73
# ) )
+
0
'
)
(
& %
( 0
&(
)
-
(
/
#& & 0
%
( ,19 - 22!. &&
&
0(
# )
&
/ &
(
(
>13
:13 - 3
)
&
(
/0
&
# ) ) #
(
&
(
&
(-
3
& (
) & &-
0
) )
# ) 0
# )
) #(& /
- 3 /
#
&(
) &
-
&-
)
(
&
&
4
$ #
4
&& 0
#
4 && 0
)
' 4 0
)
# ) - "
(
4 /
)
' 4 ,B&
O
22E7 6 - 22.- "
0 #
74 Chapitre 4. Auto-allumage
"
& ) )
0
(- "0
'
)
) &
(
- "
,:13.- 3 ( )
#& /0 )
+
& # 0
#
) - "
&
( )
# ) (
0
( / - "
,>13.
)
) %
#& (-
&
4.2. Simulations 75
& ) 0
:13 - "
)
) )
( ) >13-
' ( / (
(
0
:13
0
(
& / 0& /(
) &
/
( & '-
"
) L0 0+
&)( 90
&
-
0
/( 3B3
) ?-
76 Chapitre 4. Auto-allumage
A
:- - : - A9
?- L ,.- ) /
*
& & 4
-
EPE@-
A 1- ,22.-
- L
-
! P@@-
B& L-
>- O
,22E.- " 4
( #
/
%& ) %#
/ #- P
"
-
:/ -
B- "
,.- & ) / %#
4
&
&
-
! P!-
3 :- Q-
- >0 ,[email protected] 3& )
(
&% #& ( /
40-
# !P E-
6 - - (
Q- 1 ,[email protected]
9 33$ ) 9& )
(
) &%
9- $ %
& '
$(%)*+,,*- -
6 L- A- && 6- A(
- >R ,22.-
#&
($
4-
*. ,. !EP2-
6& - S- :
- + ,22 .- "
)
&
/
)
& - $(
*.,# -
6& - - 1 - +
1-
,.- 1
&
#
)
& / &
9 0
)
-$(
,,!+,!+!," -
6+ "-
- >0 ,22E.- 3& %#
40&
(
&
$
( 0
(-
!,/ !P
! -
6
6-
- 9 ,2@.- $ 0 / 2-
> - ,22.-
'
0
%9 -
!,1 EP-
>+ - - 6
Q- 1 ,2@@.- ?63?7 ) & )
&
& &
9 0
( (- "
@!%@@
? >/-
19 - "- A
"- L ,22!.- ? ) &%
/ #& 40-
!,* 2@P-
4.2. Simulations 77
2 P-
#
78 Chapitre 4. Auto-allumage
methane flame
nozzle
& $ 3
# & %&$
%
/- *
/
&
& )
/
/ ,
0.
&
/
0
&-
9& , . 9&
1# , .
)
& $
& )
)
#& & ,>13. /-
4.2. Simulations 79
1000 0.1
100
Yetter (1991)
10
Ignition delay (s)
0.1
1E-3
0.01
GRI Mech 3.0
GRI Mech 1.2
1E-3 12-step (Sung et al. 1998)
EXP (Krishman and Ravikumar 1981)
1E-4 1E-4
0.0 0.1 0.2 0.3 0.4 0.00 0.05 0.10 0.15 0.20
Mixture fraction Mixture fraction
E
80 Chapitre 4. Auto-allumage
3000
t=0 ms
t=0.33 ms
t=0.35 ms
t=0.37 ms
t=0.41 ms
2000
Temperature (K)
t=0.83 ms
1000
0
0.00 0.05 0.10 0.15 0.20
x (cm)
,.
3000
t=0 ms
t=0.82 ms
t=1.18 ms
t=1.28 ms
t=1.33 ms
2000
Temperature (K)
t=1.54 ms
1000
0
0.00 0.05 0.10 0.15 0.20
x (cm)
,/.
!
4.2. Simulations 81
1.0
0.8
Normalized temperature
0.6
Hydrogen
0.4
Methane
0.2
0.0
0 1 2 3 4 5
t / tHMI
,.
1.2
1.0
Heat release (x10 Jm s )
−3 −1
0.8
9
0.6
Hydrogen
Methane
0.4
0.2
0.0
0 1 2 3 4 5
t / tHMI
,/.
@
82 Chapitre 4. Auto-allumage
3.0
2.5
Methane
2.0
tign / tHMI
1.5
1.0
Hydrogen
0.5
0.0
0.00 0.05 0.10 0.15 0.20
d (mm)
& !$ > &
(
) )
0
)
#& (-
&
(
&
%
,E.
,!.- "
(
&
/(
&
&
%( >0 / +
/(
-
−2
10
−4
10
Mass fractions
−6
10
−8
10
YO (HMI)
YH (HMI)
−10 YO (LMI, x=xign)
10
YH (LMI, x=xign)
YO (LMI, max. heat release)
YH (LMI, max. heat release)
−12
10
0 0.5 1 1.5
t / tHMI
& @$ " )
)
5
: )
: J?%
- "
) )
:13
>13
)
& - (/
>13 )
) #
-
2
4.2. Simulations 83
10 10
Le=1
Le=1
tLMI / tHMI
tLMI / tHMI
1 1
Detailed transport Detailed transport
0.1 0.1
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
LeH2 LeCH4
,. ,/.
& 2$ 34 )
>0 / )
)
>13 ) ,.
(
& 4
,/.
4- &
/
/(
'- )
&
( )
'7
%
( >0 / ) -
0.8 0.3
0.6 zst
0.2
Mixture fraction
Mixture fraction
0.4 zst
0.1
0.2
zmr
zmr
0.0 0.0
0.0 0.5 1.0 1.5 2.0 0 1 2 3 4
t / tHMI t / tHMI
,. ,/.
& $ 1# ) # )
) #(&-
)
(
& ,.
,/. ) )
+
/(
&
(
&
,E.-
84 Chapitre 4. Auto-allumage
4.3 Conclusions
Ce chapitre a permis de mieux comprendre le mécanisme de l’auto-allumage dans le cas parti-
culier des moteurs-fusées. Des simulations de type HMI et LMI réalisées pour l’hydrogène mettent
en évidence l’impact de la diffusion différentielle. Dans le cas de l’hydrogène, le mélange HMI le
plus réactif contrôle l’induction. Le faible Lewis du carburant diminue l’effet retardant, ce qui mène
finalement à tHMI ≈ tLMI . Ainsi, tHMI est un bon estimateur du temps d’allumage. Cette propriété
sera utilisé lors de la création de schémas réduits capables de modéliser le temps d’auto-allumage de
l’hydrogène par un oxydant chaud. Cette conclusion fort arrangeante ne s’applique pas au méthane,
où tHMI ≈ 2.7tLMI Au final, chaque réducteur requiert une étude préliminaire de type LMI, car
la thermodynamique et le nombre de Lewis influencent substantiellement, et de plusieurs façons
contradictoires, la valeur de tLMI . Cette quantité peut dépasser d’un ordre de grandeur tHMI , ou bien
y être inférieure.
Chapitre 5
5.1 Objectif
Il s’agit d’adapter la simulation aux grandes échelles telle qu’elle est calculée dans AVBP aux
structures supersoniques transitoires. La technique de capture de choc est l’approche la plus utilisée
actuellement, et de nombreuse recherches sont en cours pour augmenter encore son applicabilité.
Cependant, l’objectif final est de calculer des jets réactifs. Dans ce cadre, il est possible de capturer
la structure globale d’un jet sous-détendu avec des schémas du type de ceux implantés dans AVBP.
85
86 Chapitre 5. Jets supersoniques sous-détendus
Zone de
mélange
Cisaillement
Frontière Frontière
réelle non Disque
visqueuse de Mach
Ondes
de choc
Coeur obliques
potentiel
w
U
t=t0 Solution initiale
t=t0+dt Dispersion
t=t0+dt Dissipation
La positivité, ou monotonicité est l’aptitude d’un schéma à interdire la création de valeurs non
physiques.
La dissipation, ou erreur de l’amortissement des hautes fréquences, apparaı̂t lorsque l’erreur
produite par le schéma est d’ordre pair, c.a.d. lorsque le schéma possède un ordre de précision
impair (premier, troisième, cinquième ordre). Il faut noter que l’opérateur de diffusion, qui calcule
l’effet de la viscosité du fluide, dissipe les hautes fréquences de la même manière. Ainsi, un schéma
dissipatif induit une viscosité ’artificielle’.
Le contrôle de ces erreurs peut se faire suivant de nombreuses approches, que l’on choisit de
réunir en deux grandes familles. L’approche de la viscosité artificielle non-linéaire vien de Von
Neumann et Richtmyer (148). Elle consiste à ajouter à un schéma un terme de viscosité artificielle
en volume. La diffusion supplémentaire qui en résulte doit être suffisante pour assurer la positivité
de la solution près des discontinuités, tout en s’annulant dans le reste de l’écoulement.
L’approche par les solveurs de Riemann s’inspire directement de la résolution d’un problème
hyperbolique par la méthode des caractéristiques. Par leur construction, ces outils doivent se res-
treindre strictement aux parties hyperboliques. Ces schémas sont précis. En effet une discontinuité
peut être calculée à l’intérieur d’une seule cellule, sans aucune erreur d’évaluation. Cela explique
l’appellation familière de ’schéma à capture de choc’. Par contre, ils contiennent implicitement l’hy-
pothèse d’un gaz parfait non réactif. Les principales sous-familles des solveurs de Riemann sont :
– Le schéma de Godunov (67). C’est le plus exact. Il produit des solutions d’un problème de
Riemann local à chaque cellule, au prix d’un temps de calcul prohibitif.
– Les ’solveurs de Riemann approchés’ qui s’appprochent d’un schéma de Godunov avec un
temps de calcul plus acceptable, comme les schémas de Roe (126).
88 Chapitre 5. Jets supersoniques sous-détendus
– La méthode MUSCL (145) (Monotonic Upwind Scheme for Conservative Laws) proposée
par Van Leer permet de porter la précision des solveurs de Riemann au second ordre par
reconstruction de la solution. Le second ordre implique la création d’oscillations, que l’on
peut détruire à leur source par l’utilisation de limiteurs.
– Les méthodes ENO, puis WENO (Weighted Essentially Non Oscillatory) qui poussent un peu
plus loin le concept MUSCL, en autorisant la reconstruction de la solution jusqu’au sixième
ou au huitième ordre (78). L’ensemble des travaux de Shu (78)donne un bon aperçu de cette
classe de méthodes.
Il est préférable de limiter l’utilisation des solveurs de Riemann aux régions avec discontinuités,
et de choisir des schémas moins dissipatifs et plus rapides dans le reste de la simulation : c’est une
approche hybride. Pullin et Hill ont proposé en 2004 un schéma centré adapté à la LES, le TCD
(Tuned Centered-Difference) associé au schéma WENO près des chocs (70). Selon D.I. Pullin, le
passage d’un schéma à l’autre est encore heuristique. Ajoutons à cela que les schéma WENO sont
encore rarement utilisés en maillages non structurés tridimensionnels (162), et leurs formulation en
gaz réel est encore récente (107). Dans le cadre de cette thèse, on utilisera donc une approche de
type Von Neumannn-Richtmyer.
∂u j
σi j = τi j − pδi j = µ ∂xi + ∂ui
∂x j + λ ∂u
∂xk δi j
k
− pδi j
(5.2)
Tenseur de contraintes Viscosité de cisaillement Viscosité de compression Pression
On introduit le paramètre de viscosité de volume κ = 2/3µ+λ. La théorie cinétique des gaz (72)
montre que ce dernier paramêtre est directement lié à l’indice de compressibilité Z par l’Eq, 5.3.
!
P
κ = 2/3µ + λ = µ1.002 − 1 = µ1.002 (Z − 1) (µ > 0), (5.3)
ρRT
5.4. Viscosité de Von Neumann-Richtmyer 89
Pour un gaz parfait , i.e. Z = 1, la viscosité de volume κ ainsi définie devient nulle, ce qui donne
l’hypothèse de Stokes :
κ = 2/3µ + λ = 0 (µ > 0) (5.4)
On remarque que la viscosité de compression est comparable à une pression. La divergence de
la vitesse est le terme source, tandis que λ dépend du fluide et peut être positif ou négatif. Dans le
cas d’un choc hydrodynamique, cette viscosité influence directement l’épaisseur du choc.
En 1949, Von Neumann et Richtmyer (148) proposent une viscosité en compression localement
prépondérante autour des chocs. D’un point de vue formel, un tenseur de contraintes vérifiant la
relation de Stokes, avec une augmentation λ0 de la viscosité en compression, est introduit (Eq. 5.5).
En interprétant cette augmentation comme un supplément local de pression p0 (Eq. 5.6), le tenseur
des contraintes peut se ré-écrire comme dans l’Eq. 5.7.
0 ∂uk ∂u j ∂ui
!
2
σi j modi f ied = (− µ + λ ) +µ + − pδi j (5.5)
3 ∂xk ∂xi ∂x j
∂uk
p0 = −λ0 (5.6)
∂xk
∂uk ∂u j ∂ui
!
σi j modi f ied =λ +µ + − (p + p0 )δi j (5.7)
∂xk ∂xi ∂x j
Il reste à determiner la valeur de λ0 , ou plutôt de p0 . Pour garantir une influence négligeable
dès que l’on s’éloigne du choc, la viscosité artificielle de compression dépend de la divergence de
vitesse au carré. L’ensemble est adimensionné par la taille de la cellule et la densité. On obtient une
forme telle qu’exprimée dans l’Eq. 5.8, ajustée grâce à une constante sans dimension CV NRV . Pour
l’ensemble de ces travaux, CV NRV = 1 s’est illustré comme étant la valeur optimale.
∂uk ∂uk
q
p0 = −λ0 = −CV NRV ρ∆x2 5 · ~u 2 (5.8)
∂xk ∂xk
Dans un code compressible à n dimensions, cette viscosité peut se calculer de plusieurs manières.
On retiendra l’Eq. 5.9 qui conserve la masse par construction. Cependant, l’expression est identique-
ment nulle dans le cas mono-dimensionnel. On retiendra alors également l’Eq. 5.10. Ces viscosités
seront désignées par l’acronyme ”VNRV” pour ”Von Neumann-Richtmyer Viscosity”.
0 ∂uk ∂u j ∂ρui
pi = −λi
0
= −CV NRV (Volcell ) max
2/3
(5.9)
∂xk j=1,n ∂x j ∂xi
0 ∂uk ∂u j ∂ui
pi = −λi
0
= −CV NRV (Volcell ) max
2/3
ρ (5.10)
∂xk j=1,n ∂x j ∂xi
Ce modèle simple agit sur les gradients de pression avec au opérateur du deuxième ordre. Les
modèles plus élaborés de l’hyperviscosité agissent sur les gradients de température et d’espèces,
avec des opérateurs d’ordre élevé, comme par exemple le modèle de Fiorina et Lele en 2006 (62)
sur maillage structuré.
90 Chapitre 5. Jets supersoniques sous-détendus
t=0 t = 1.93tc
Density Velocity Density Velocity
1.50 100 1.50 100
1.40 80 1.40 80
1.30 60 1.30 60
1.20 40 1.20 40
1.10 20 1.10 20
1.00 0 1.00 0
0.0 2.0 4.0 6.0 8.0 10.0 0.0 2.0 4.0 6.0 8.0 10.0 0.0 2.0 4.0 6.0 8.0 10.0 0.0 2.0 4.0 6.0 8.0 10.0
t = 3.86tc t = 19.32tc
F. 5.3 – L’onde de Riemann. Schéma au troisième ordre avec viscosité de Von-Neumann Richt-
myer. Solution analytique (trait plein), Simulation (trait interrompu)
5.4. Viscosité de Von Neumann-Richtmyer 91
La Fig. 5.3 montre l’évolution de l’onde de Riemann, d’après la solution analytique issue de
Landau & Lifchitz (91) jusqu’à presque 20 temps convectifs. Un schéma centré de troisième ordre
stabilisé par la viscosité de Von-Neumann Richtmyer suit bien la solution analytique. Un léger
tremblement est visible au voisinage gauche du choc à t = 1.93tc puis diminue jusqu’à disparition
complète à t = 19.32tc . On note un écart de position d’environ 40cm, après avoir parcouru 200m,
soit un décalage de 0.4%.
Afin de se convaincre de l’efficacité de la VNRV, le schéma à été appliqué seul à l’onde de
Riemann. A t = 19.32tc , le phénomène de Gibbs a totalement détruit la solution, comme le montre
la Fig. 5.4.a. Avec une viscosité artificielle plus classique (76) basée sur l’augmentation directe
de la viscosité cinématique, une solution stable est obtenue (cf. Fig. 5.4.b) mais des oscillations
sont encore visibles à gauche du choc. Ainsi, la viscosité artificielle classique est suffisante pour
stabiliser le schéma sur une solution continue, mais la qualité de la solution devient insuffisante en
cas de discontinuités.
1.00 0 1.00 0
0.0 2.0 4.0 6.0 8.0 10.0 0.0 2.0 4.0 6.0 8.0 10.0 0.0 2.0 4.0 6.0 8.0 10.0 0.0 2.0 4.0 6.0 8.0 10.0
F. 5.4 – L’onde de Riemann. Solution analytique (trait plein), Simulation (trait interrompu), à
t = 3.86tc
92 Chapitre 5. Jets supersoniques sous-détendus
80000.0
600.00
60000.0
400.00
40000.0
200.00
20000.0
0.0 0.00
0.0m 2.0m 4.0m 6.0m 8.0m 10.0m 0.0m 2.0m 4.0m 6.0m 8.0m 10.0m
3
T [K] ρ [Kg/m ]
800.0 1.0
t=0ms
0.8
600.0 t=1ms
t=2ms
0.6 t=3ms
400.0 t=4ms
0.4 t=5ms
200.0
0.2
0.0 0.0
0.0m 2.0m 4.0m 6.0m 8.0m 10.0m 0.0m 2.0m 4.0m 6.0m 8.0m 10.0m
F. 5.5 – Evolution d’un choc dans un tube à choc, solution analytique : Pamont /Paval = 100,
Température initiale constante
Si l’on se concentre sur le choc mobile issu du tube à choc, on voit sur la Fig. 5.6 que l’approche
de Von Neumann-Richtmyer est efficace, en maintenant une dizaine de points dans l’épaisseur du
choc. Le schéma centré de troisième ordre TTGC doté d’une viscosité artificielle classique com-
mence à osciller sérieusement autour du choc. Une fois que la viscosité VNRV remplace la viscosité
artificielle classique, on peut constater la positivité de la solution.
Ainsi, l’approche de Von Neumann-Richtmyer présentée dans cette thèse fonctionne sur des
chocs forts. Elle autorise l’utilisation de schémas centrés au troisième ordre sur des chocs mobiles
forts. De plus la viscosité est active seulement lorsque la divergence de la vitesse est forte.
5.4. Viscosité de Von Neumann-Richtmyer 93
P [Pa] u [m/s]
40000.0 600.0
TTGC
TTGC + VNRV
30000.0 400.0 Solution analytique
20000.0 200.0
10000.0 0.0
0.0 −200.0
8.0m 8.5m 9.0m 9.5m 10.0m 8.0m 8.5m 9.0m 9.5m 10.0m
3
T [K] ρ [Kg/m ]
500.0 0.3
450.0
0.2
400.0
0.1
350.0
300.0 0.0
8.0m 8.5m 9.0m 9.5m 10.0m 8.0m 8.5m 9.0m 9.5m 10.0m
a) Pamont /Paval = 10
P [Pa] u [m/s]
10000.0 1000.0
TTGC
TTGC + VNRV
800.0
8000.0 Solution analytique
600.0
6000.0
400.0
4000.0
200.0
2000.0
0.0
0.0 −200.0
8.50m 8.75m 9.00m 9.25m 9.50m 9.75m 10.00m 8.50m 8.75m 9.00m 9.25m 9.50m 9.75m 10.00m
3
T [K] ρ [Kg/m ]
800.0 0.04
700.0
0.03
600.0
500.0
0.02
400.0
300.0
0.01
200.0
100.0 0.00
8.50m 8.75m 9.00m 9.25m 9.50m 9.75m 10.00m 8.50m 8.75m 9.00m 9.25m 9.50m 9.75m 10.00m
F. 5.7 – Couche de mélange à fort gradient de densité : Pression 8 atm., Hélium à 10.6 m/s, Azote
à 4 m/s, ρN2 /ρHe = 7.
5.5. Outils d’analyse 95
2
Distance [m]
0
1 2 3 4 5 kHz
Frequence [Khz]
Ici, f est la fréquence. On obtient ainsi la matrice des interspectres ψ( f ) jl obtenus sur une ligne de
n sondes Xi=1,n placées en ( x~j , x~l ) dans la direction de développement de l’écoulement :
Le vecteur des nombres d’ondes k associés à chaque sonde et noté ~η(k) j = eikx j : On introduit
l’écart spatial entre les sondes en multipliant la matrice ψ( f ) à gauche et à droite par le vecteur
des nombres d’ondes ~η(k) . On notera en particulier que les sondes peuvent être irrégulièrement
réparties : X
FNO( f, k) = ~ηT (k) ψ( f ) ~η(k) = ψ( f ) j,l eikx j eikxl (5.13)
j,l
Plus localement, l’identification d’une même sinusoı̈de (fréquence f , amplitude β, vitesse ~u) de
la forme βcos(2π f (x/u − t)) dans les signaux issus d’un doublet de sondes séparés de la distance
~ conduit à l’ Eq. 5.14, qui fait intervenir u x , projection de ~u sur le support de ∆x
∆x, ~ dans la direction
orthogonale de propagation. Si l’onde de vitesse u se propage dans une direction faisant un angle α
par rapport à la ligne de sondes, alors u x = u/cos(α), cf. Fig. 5.9.
x + ∆x
!! ! !
x
βcos 2π f − t = βcos 2π f − t + ∆φ (5.14)
ux ux
Dans cette expression, ∆Φ représente le décalage en temps. Par définition ∆φ = 2π f ∆x/u x , soit :
2π f
ux = (5.15)
k
où k = ∆φ/∆x. Afin d’éviter le phénomène de repliement, le nombre d’onde est limité à 2π/2∆x.
La résolution en fréquence fres est imposée par le nombre d’enregistrements nrec et le pas de temps
∆t selon l’Eq. 5.16, et se ramène à la moitié de l’inverse de la durée d’enregistrement trec .
fNyquist 1
fres = = = 0.5/trec (5.16)
nrec 2 ∗ nrec ∆t
5.5. Outils d’analyse 97
α u
u/cos(α )
F. 5.9 – Onde sinusoı̈dale de vitesse u traversant une ligne de sondes avec un angle de déviation
α.
Ainsi, la vitesse minimale observable umin sur un domaine fréquentiel [ fres , N fres ] (N de l’ordre de
100) est donnée par l’Eq. 5.17
N∆x
umin = (5.17)
trec
La vitesse minimale observable impose le temps de simulation et l’espacement des sondes. Par
ailleurs, N fres doit être dans le domaine fréquentiel résolu, imposant ∆t < 1/(N fres ). Le choix de
ces deux paramètres doit se faire avant de réaliser la simulation.
La Fig. 5.10 donne un exemple de résultat typique de l’analyse fréquence-nombre d’onde. On y
voit des structures appartenant à la gamme de fréquences [0, 3000Hz] convectées parallèlement aux
sondes à la vitesse de 2π5000/314 = 100m/s. De la même manière, une onde à 1kHz et ses har-
moniques à 2, 3, 4, 5kHz se propagent vers l’amont et l’aval aux vitesses u ± c = 100 ± 360m/s.
Il est intéressant de noter que les structures de même fréquence (par ex. 1kHz) se propageant
en des sens opposés ont pu être distinguées. De plus, la présence de la structure advectée sur la
gamme [0, 3000Hz] ne cache pas les structures à 1, 2 et 3kHz. L’analyse fréquence nombre d’onde
décompose ainsi les signaux temporels par rapport à leur vitesse de convection. Autrement dit, cet
outil permet d’identifier les vitesses de propagation de différentes structures, même si elles corres-
pondent à une même fréquence. La direction de la ligne de sondes joue un rôle crucial.
98 Chapitre 5. Jets supersoniques sous-détendus
68 u+c
0
−120 u−c
1 2 3 4 5
Frequence [kHz]
5.5.3 Bicohérence
Les simulations aux grandes échelles peuvent reproduire des processus complexes d’interactions
non linéaires entre les phénomènes dynamiques de différentes fréquences. La bicohérence permet
d’observer en détail ces processus. D’un point de vue spectral, deux fréquences sources f1 et f2
peuvent se combiner non linéairement pour alimenter les fréquences cibles f1 + f2 et f1 − f2 . L’inter-
spectre du produit de deux sondes sources F̂( f1 )Ĝ( f2 ) et du signal d’une sonde cible Ĥ( f1 + f2 ) met
statistiquement en évidence ces interactions. L’Eq. 5.18 donne l’estimateur de bicohérence complet.
Le numérateur comporte l’interspectre entre le produit des deux TF sources F̂( f1 )Ĝ( f2 ) multiplié
par la TF cible Ĥ ∗ ( f1 + f2 ) conjugée. Cet interspectre est moyenné sur tous les triplets de sondes,
une opération notée h i. L’opération de moyenne conserve la norme du numérateur si son argument,
i.e. le déphasage source/cible, est fortement corrélé. On garde ensuite le module du numérateur,
opération notée | |. Enfin, on effectue une normalisation d’après la moyenne du signal source et du
signal cible, visibles au dénominateur.
D E 2
F̂( f1 )Ĝ( f2 )Ĥ ∗ ( f1 + f2 )
Bic2FGH ( f1 , f2 ) = (5.18)
F̂( f )Ĝ( f )2 Ĥ( f + f )2
1 2 1 2
L’estimateur de bicohérence est non nul si une interaction entre deux signaux sources est observée
sur le signal cible, avec un déphasage constant. La normalisation permet d’observer les interac-
tions entre des fréquences d’amplitude très variées. Par contre, l’estimateur dégénère en forme
indéterminée lorsque des fréquences observées f1 et f2 ont des amplitudes nulles.
La bicohérence se représente sur une carte à deux dimensions f1 et f2 . Une troisième dimension
implicite est présente : f1 + f2 . Les interactions peuvent concerner des additions (100 + 10 = 110Hz)
ou des soustractions (100 − 10 = 90Hz) dans l’espace fréquentiel. C’est pourquoi, sur un domaine
fréquentiel donné [0 − 1000Hz], le domaine d’évaluation de la bicohérence s’étend sur des valeurs
fréquentielles algébriques : [−1000Hz, 1000Hz] × [−1000Hz, 1000Hz]. On notera que la fréquence
nulle correspond à la composante stationnaire, c.a.d. l’écoulement moyen. Ce domaine contient
cependant deux fois les mêmes informations car Bic( f1 , f2 ) = Bic(− f1 , − f2 ). Les positions telles
que f1 + f2 < 0 sont donc ignorées. Par ailleurs, si les signaux sources f1 et f2 sont interchangeables,
on a Bic( f1 , f2 ) = Bic( f2 , f1 ), ce qui annule la pertinence du domaine f2 < f1 .
Trois cartes de bicohérences schématiques permettent d’illustrer le propos :
– La fig. 5.11.a montre plusieurs interactions (symbole •) très localisées. La première (0Hz +
1kHz → 1kHz) traduit l’alimentation d’un mode à 1kHz par l’écoulement moyen (0Hz).
Ce mode fondamental interagit avec lui même pour créer le premier mode harmonique à
2kHz (1kHz + 1kHz → 2kHz). Le deuxième mode harmonique est ensuite nourri par le
fondamental et son premier harmonique (1kHz+2kHz → 3kHz), et ainsi de suite. Finalement,
la bicohérence montre comment un mode interagissant non linéairement avec lui même peut
créer toute une famille d’harmoniques.
– La fig. 5.11.b correspond au cas opposé. Un fondamental et ses harmoniques (1,2,3,4,5 kHz)
sont alimentés par l’écoulement moyen. Comparé à cet apport d’énergie, il n’y a pas d’inter-
action non linéaire entre les modes.
– La fig. 5.11.c montre trois lignes d’interactions qui alimentent chacune un mode : 1, 3 et 9kHz.
Une ligne ininterrompue implique l’excitation du mode par l’ensemble du spectre. Dans ce
100 Chapitre 5. Jets supersoniques sous-détendus
cas particulier, trois modes particulièrement réceptifs recoivent de l’énergie par un signal au
spectre plein (bruit blanc, dirac).
Une carte de bicohérence révèle le scénario des interactions dans son intégralité. Elle affirme ou
exclut les liens entre deux phénomènes, et détermine le sens du transfert d’énergie. En un mot, elle
sert de ”test de parternité” parmi les modes présents.
10
4
3
9
F2 [kHz]
2
1 8
−4 −3 −2 −1 1 2 3 4 6
F1 [kHz]
5
4
3
2
F2+F1 [kHz]
1
10
4
3
9
F2 [kHz]
2
8
1
−4 −3 −2 −1 1 2 3 4 6
F1 [kHz]
5
4
3
2
F2+F1 [kHz]
1
F2 [kHz]
3
1
−9 −3 −1 1 3 9
F1 [kHz]
F2+F1 [kHz]
3
1
c) Excitation de trois modes particuliers par l’ensemble du spectre (par [Link] bruit blanc)
Le résultat de l’estimateur dépend avant tout des données qu’il traite. C’est pourquoi l’utilisation
de la bicohérence se fait par rapport à une question très précise. Par exemple : ”Dans cette chambre,
est-ce le mode tournant acoustique à 400Hz ou le tourbillon de jet PVC de même fréquence, ou
encore une combinaison des deux qui forment des fluctuations à 800Hz en sortie ?”. La question
impose la localisation des données sources. La Fig. 5.12 montre schématiquement où placer les
sondes avec les points F pour capter l’injection, G pour le mode tournant et H pour la sortie. Une
fois les signaux obtenus, on réalise une bicohérence entre le jet et lui même (FFH, Fig. 5.12.a), une
bicohérence entre le jet et le mode tournant (FGH, Fig. 5.12.b) et une bicohérence entre le mode
tournant et lui même (GGH, Fig. 5.12.c). On notera que dans le cas b), les sources F et G n’étant
pas intercheangeables, il est nécessaire de représenter F2 < F1.
L’interaction
1 du cas a) signifie que le jet alimente seul le mode à 800Hz en sortie. L’interac-
tion 2 du cas b) révèle que c’est la conjonction entre jet et mode tournant qui alimente la sortie.
L’interaction 3 du cas c) signifie enfin que le mode tournant sollicite seul la sortie. Par conséquent,
si l’analyse des trois bicohérences révèle que l’on a l’interaction
1 mais pas
2 et
,
3 alors on en
conclut que le mode de sortie est du aux interactions non-linéaires du jet d’injection avec lui-même.
Le mode tournant ne joue aucun rôle direct (
) 3 ou indirect (
).
2 Cet exemple montre comment il est
possible de prouver qu’un phémomène dynamique est à la source d’une instabilité à l’intérieur s’une
simulation numérique. La simulation des grandes échelles devient grâce à ce type d’exploitation un
outil d’investigation actif.
a)
400 1 400
H H
F
F
80
F H
0
H
40
0
−400 400
b) G
G 400 2
H H
F
F
80
F H
0
H
G −400
40
0
G
−400 400
c) G
G 400 3 400
H H
80
H
0
H
G
40
0
G
−400 400
7.5
2.5
5 10 15 [ms]
F. 5.13 – Trajectoire des tourbillons observés. Pression 7 atm., Hélium à 10.6 m/s, Azote à 4 m/s,
ρN2 /ρHe = 7.
5.5. Outils d’analyse 103
Commençons par la position des différentes structures. La Fig. 5.14 montre ce que donnerait une
évolution spatiale longitudinale de spectre pour la couche de mélange compressible, conformément
au scénario que nous avons retenu. Elle schématise le résultat obtenu d’après les signaux de sondes
placées longitudinalement dans la couche de mélange. Il est aisé d’en déduire les localisations spa-
tiales et fréquentielles de chaque mode :
– Mode I, Fréquence 1392Hz, Position [1cm, 5cm]
– Mode II, Fréquence 691Hz, Position [2.5cm, 5cm]
– Mode III, Fréquence 346Hz, Position [2.5cm, +∞]
Cela permet de délimiter la zone [2.5cm, 5cm], qui pourrait être le théatre d’échanges énergétiques
entre les différents modes. C’est pourquoi l’analyse par bicohérence sera alimentée par des signaux
provenant exclusivement de cette zone privilégiée.
10
7.5
Distance [cm]
2.5
0
346 691 1392
Frequence [Hz]
F. 5.14 – Shéma de l’évolution spatiale de spectre pour la couche mélange de Brown et Roshko
Poursuivons par la vitesse de déplacement des différentes structures. La Fig. 5.15 montre ce
que donnerait une analyse fréquence - nombre d’onde pour la couche de mélange compressible,
conformément au scénario que nous avons retenu. Cette analyse commence par le dimensionnement
de la ligne de sondes à placer dans l’écoulement. La vitesse de convection minimale est u = 4m/s.
Pour la trouver sur une fréquence de 2000Hz, il faut 3141, 6rad/s de nombre d’onde maximal. On
note au passage qu’à 2000Hz, 5.618m/s correspondent à 2236.8m/s, et 10.6m/s correspondent à
1185m/s. Le nombre d’onde maximal impose de disposer des sondes espacées de 1 10−3 mm. Le
résultat de l’analyse est schématisé sur la Fig. 5.15. Trois vitesses de convections se détachent
nettement. Tout d’abord, les trois classes de tourbillons sont convectées à la même vitesse, 5.6m/s.
Ensuite, deux vitesses de convection sont visibles sur toute la gamme des fréquences, à 4m/s et
10.6m/s, révélées par les inévitables perturbations que contiennent les flux d’azote et d’hélium.
L’analyse fréquence-nombre d’onde a mis en évidence la vitesse de convection très particuliére
de toutes les structures générées par la couche de mélange. Malgré leurs différences de taille, des
tourbillons issus du cisaillement ont la même vitesse de déplacement.
104 Chapitre 5. Jets supersoniques sous-détendus
+5.618m/s
1500
+10.6m/s
−1500
−3000
346 691 1392 2000
Frequence [Hz]
F. 5.15 – Shéma de l’analyse fréquence - nombre d’onde associée à la couche mélange de Brown
et Roshko
Enfin, l’analyse par bicohérence est illustrée. On retrouve dans la Fig. 5.16 les trois fréquences
346, 691, 1382Hz correspondant aux trois catégories identifiées par l’évolution spatiale de spectre
en Fig. 5.8. Pour plus de lisibilité, l’axe des fréquences cibles f1 + f2 est également tracé sur la figure.
Chaque lieu de bicohérence non nulle contribue au mode noté sur cet axe. Les structures de la classe
I présentent une interaction significative avec l’écoulement moyen :
1 ≡ 1392 + 0 → 1392Hz.
Ensuite la classe II interagit avec la classe I : 2 ≡ 1392 − 691 → 691Hz. Enfin la classe III interagit
avec la classe II : 3 ≡ 691 − 346 → 346Hz. Cet outil met en évidence la cascade energétique du
scénario retenu : écoulement moyen → Mode I → Mode II → Mode III. Le scénario de transfert
énergétique apparaı̂t sur l’analyse par bicohérence. Si il y avait eu une alimentation de tous les
tourbillons par l’écoulement moyen au lieu de la cascade, les interactions
2 et
3 auraient été
absentes. Par contre, on aurait trouvé les interactions (691 + 0 → 691Hz) et (346 + 0 → 346Hz). La
bicohérence est donc un outil puissant de compréhension des écoulements.
5.5. Outils d’analyse 105
1392 2 1
F2 [Hz]
691 3
346
1392
5.6 Simulations
Les écoulements supersoniques sont difficiles à aborder par la simulation des grandes échelles,
car on y rencontre des nombres de Reynolds très élevés. Dans ce cas, seule une petite partie du
spectre de la turbulence est calculée. Cependant, l’article soumis à Physics of Fluids et reproduit ici
montre que pour une simulation de cette courte gamme de Reynolds, une grande quantité d’infor-
mation est cependant disponible. Pour mieux illustrer la difficulté du problème, une l’étude dimen-
sionnelle en termes de nombres de Reynolds est présentée ci-dessous.
Sachant que l’échelle intégrale lt subit une dissipation identique, on a lt = u03 /. Il vient donc la
relation entre les échelles de Kolmogorov η et intégrale lt :
lt u03 /
= 3 1/4 = Re3/4 (5.20)
η (ν /) t
Les résultats de l’analyse dimensionnelle pour les cas des jets supersoniques traités dans ce chapitre
sont reportés dans le tableau 5.1. Les vitesses u0 et tailles lt caractéristiques sont estimées respec-
tivement à 5% U et 30% D. L’échelle de Kolmogorov est de l’ordre de 6µm pour les deux cas. Les
Reynolds Re en jeu sont très élevés, donc les grandes échelles de cet écoulement ne sont pas affectés
5.6. Simulations 107
par la dissipation visqueuse. La turbulence du jet de Seiner se développe sur une décade de plus que
celle du jet de Chuech. Les études quantitatives se feront de préférence sur ce dernier.
Concernant la compressibilité, le Mach convectif doit se calculer p en tenant compte de l’accéléra-
tion du fluide due à la pression d’injection, approchée par ∆U = P/ρ. On trouve ainsi les valeurs
0.679 et 1.39 pour les Mach convectifs respectifs des jet de Chuech et Seiner, ce qui signifie que les
effets de compressibilité sont présents dans la turbulence, tout particulièrement dans le jet de Seiner.
Ainsi, les deux configurations retenues sont des jets fortements turbulents et compressibles.
5.6.3 Publication
L’article reproduit ci-après est en cours de soumission au journal Physics of fluids. Il illustre
l’application de la viscosité artificielle de Von Neumann-Richtmyer à un jet supersonique en SGE
et les outils d’analyse spectrale introduits précedemment.
108 Chapitre 5. Jets supersoniques sous-détendus
Ret (r)
10000
1000
Seiner
Coupure
100
Chuech
10
r
−6 −5 −4 −3 −2 −1
10 10 10 10 10 10
F. 5.17 – Gammes de Reynolds modélisées et résolues lors de la simulation des jets de Chuech et
de Seiner.
5.6. Simulations 109
Abstract. Weakly underexpanded jets have been investigated numerically, with a compressible explicit Large Eddy
Simulation code. The study is mainly focused on the near injection zone. The LES approach gives a correct description
of the flow compared to experimental results. Three types of injections has been applied to the flow: laminar, uncorre-
lated white noise and correlated modes. Each injection type has an impact on the shockwaves, giving insight on their
acoustic activity. Spectral analysis enhances the frequency response of the flow to the pattern injected. Bicoherence
statistics proves the non-linear production of harmonics. Spectral tracking of coherent structure showed the strong
penetration of upstream oblique waves in the supersonic core of the jet. This work points out the critical impact of
turbulence injection in such jets, in particular concerning the unsteady behaviour of the shock structure.
I. INTRODUCTION
A. Context
An underexpanded jet is the flow resulting from the discharge of a high pressure tank. In this flow, the gases are
expanded to an intermediary low pressure state, then recompressed to the ambient state. This flow is encountered in
various applications. In fundamental physics, particle accelerators use strongly underexpanded jets. Indeed, a small
amount of particles must be injected with the same velocity, the same direction. These requirements are obtained
with the intermediary low pressure state of underexpanded jets, providing a rarefied gas with a good homogeneity of
velocities [22].
In aeronautical turbines, particular care is taken of the fatigue of engines. After years of operations, some cracks
appear in the walls of the engine. This concerns in particular high pressure compressor stages, combustion chamber,
and high pressure turbine. Because of the pressure gap betwen the two sides of the wall, possible flow leaks become
underexpanded jets. Like a high pressure water jet can cut plastic sheets, the resulting jets can destroy devices nearby.
It is a matter of safety, and the operation or the control of the engine can be altered [15].
In rocket launchers, cryotechnic igniters starts the combustion of second stage engines (e.g. Vinci) with the au-
toignition of two reactants. Fuel and oxidizer, injected in near vacuum conditions, are mixed in the interaction of two
underexpanded jets. These configurations are focused on the near-injection region, in the first diameters of the jet
development.
Underexpanded jets are generated by an injector fed with a gas at a high pressure compared to the ambiant
pressure. After injection, the fluid is expanded through a complex pattern of expansion/recompression zones, The
Nozzle Pressure Ratio (NPR), i.e. the static pressure at the nozzle divided by the ambiant pressure, and the Mach
number at the injection are both essential parameters. For this study, we will consider only sonic jets, i.e. without
diverging nozzle at the exit.
Fig. 1 shows a typical sonic underexpanded jet at NPR 4.7 [34], with a Schlieren visualization done at the Institute
of Theoretical and Applied Mechanics (ITAM) of Novossibirsk. A first cell of expansion is clearly visible, but the flow
is blurred drownstram by the growth of turbulent motion. This picture illustrates distinctly the struggle between
inviscid structures and turbulence.
In Case of sonic jets with a NPR slightly superior to unity, the jet is said weakly underexpanded (Fig. 2, N P R ≈ 1.2),
showing typical series of several expansion/compression cells often refered as ’diamond shock pattern’. For higher
pressure ratios, the cells become barrel shaped [5, 28] because of the reflection of the expansion fan on the slip line
(Fig. 2, N P R ≈ 3) [4]. These structures can be commonly observed on the exhaust of jet engines, especially in
afterburner regime. The strongest compression zones lie in the downstream tip of compression bundles. For weakly
110 Chapitre 5. Jets supersoniques sous-détendus
FIG. 1: 1/30 sec. Schlieren picture of a sonic jet, NPR=4.7 (Courtesy of ITAM, Novossibirsk [34])
underexpanded jet, these discontinuities are located at the end of each ’barrel’, and the extremum of the compression
occurs at the downstream rim of the barrel.
For pressure ratios higher than 4, the jet is said highly underexpanded (Fig. 2, N P R ≈ 10), and presents an intense
expansion zone limited by a barrel shaped weak shock, itself ended by a strong disc-shaped shock [6, 9, 16]. This
shock pattern, sometimes referred to the ’Mach Bottle’ is the inner side of the compression zone surrounding the
jet. A thin and intense mixing layer makes the outer side of the unique expansion cell. Turbulence is created in a
convective curved zone througg the development of longitudinal vortices by the Goertler instability [35]. For NPR
rising from 2 to 4 [7], the barrels become thicker, the disc shaped shocks appear at their bottom, and their number
decreases to one.
These flows contain strong discontinuities with huge pressure jumps in the shocks, and ’slip lines’, i.e. curved and
narrow mixing layers, between supersonic flow and fluid at rest. The balance between the mixing effect of turbulence
and the inviscid segregated flow controls the development of the plume. The inviscid compressible flow near the
injection becomes downstream a nearly incompressible viscous flow. Such jets generate strong acoustic noise near
the injection. The Screech phenomenon, first observed by Powell [24] corresponds to a coupling between structures
convected downstream by the jet and acoustics travelling upstream outside of the jet [25]. This particular feature,
critical for the acoustic fatigue in high speed aircraft, has been experimentally [21] and numerically [19] analysed.
Mixing Mixing
region region
Actual Shear
jet layer
boundary
Mach
Jet disk
boundary
Oblique Invicid Oblique
Core shock jet shock
waves boundary waves
Supersonic jets have been numerically investigated very early in the 1950s with the Method of Characteristics [4, 9].
This Euler-based method gave the inviscid parameters of plumes such as the length and diameter of the barrel shock.
Finite differences schemes with shock capturing techniques allowed to capture stationary viscous flow properties [1, 7,
8], with good matching of the plume development compared to experiments. Recently, Large Eddy Simulation gave
more quantitative results on noise generated by adapted and weakly under-expanded supersonic jets [2]. The impact
5.6. Simulations 111
of turbulence is difficult to take into account in RANS as no turbulent model is available for the curved supersonic
mixing layer.
The screech phenomenon is commonly simulated with Unsteady RANS approaches [29] associated to specific tur-
bulent models [32], bidimensional [12, 14, 18] and three dimensional LES [13] and finally DNS [20] in a simplified
framewrork.
LES is usually applied to supersonic flows through a Riemann Solver Approach. The numerical schemes involved
(MUSCL, WENO) are accurate but CPU-time consuming, and often limit simulations to small and/or structured
grids. The Flux Corrected Transport (FCT) approach uses compact centered schemes and artificial viscosity. This
last strategy is rarely encountered in the bibliography, but shows some advantages. Compact centered schemes are
compatible with unstructured grids, which is useful in complex geometries. Riemann solvers approach is built upon a
gas law, and assumes that the fluid is non reactive. FCT approach avoid these hypothesis, allowing easier developments
towards combustion and real gas applications. In this frame, this paper presents the results of a FCT approach on two
under-expanded jets. The analysis of steady and unsteady features resolved shows the informations available through
these computations.
The paper is organized as follows. [Link] presents the numerical method, clarifying the advantages and drawbacks
of the explicit fully compressible approach. The simulations are validated in Sec. III. Calculated profiles of mean
and fluctuating values calculated are compared to experimental profiles of Seiner [28] and Chuech [5] on Pressure
and Velocity field. Sec. IV shows observations of a fluctuating diamond-shock pattern. It starts with the temporal
evolution of the profile of pressure on jet axis. The localization of acoustic activity for three different excitations of
the jets follows, preparing the spectral part of the paper. Then, Sec. V gathers the acoustic analysis of the jet. It
begins with the spatial evolution of frequency spectrum with downstream distance. A bicoherence analysis follows,
proving any non linear coupling between main modes. In the end, a coherent structure tracking isolates the main
acoustic fluxes and their direction of propagation.
A. Governing equations
where G denotes the filter function and f (x, t) is the filtered value of the variable f (x, t). The filter function is defined
R +∞
as G(x0 , x) ≡ g(x0 − x) with properties g(x) = g(−x) and −∞ g(x)dx = 1. In the mathematical description of
compressible turbulent flows, the primary variables are the density ρ(x, t), the velocity vector ui (x, t)and the total
energy E(x, t). The application of the filtering operation to the instantaneous transport equations yields [23].
∂ ρ̄ ∂
a) ∂t + ∂xi (ρ̄ũi ) = 0
∂ ρ̄u˜j ∂τ
b) ∂t
∂
+ ∂x i
∂ p̄
(ρ̄ũi u˜j ) = − ∂x j
+ ∂xjkk ∂
− ∂x i
ρ̄T˜ij (2)
∂ q¯
c) ∂ ρ̄Ẽ
∂t
∂
+ ∂x i
∂
ρ̄ũi Ẽ = − ∂xjj + ∂x j
∂
[(τij − pδij )ui ] − ∂x j
∂
ρ̄Q̃j − ∂x j
ρ̄T˜ij ũi
In 2, one uses the Favre-filtered variable [10] f˜ = ρf /ρ̄. The fluid follows the ideal gas law p = ρRT and es = Cv T
where T stands for temperature. The tensors of viscosity and Von Neumann-Richtmyer artificial viscosity, and the
heat diffusion vector read respectively
∂u
τij = µ ∂x ∂ui
+ ∂xji + 32 ∂u
∂xk δij
k
j
∂T
(3)
qi = −λ ∂x i
In 3, µ is the fluid viscosity following Sutherland’s law and λ is the heat diffusion coefficient following Fourier’s
law. The objective of LES is to compute the largest structures of the flow (these structures are typically larger
112 Chapitre 5. Jets supersoniques sous-détendus
than the computational mesh size), whereas the effects of the smaller scales are modelled. This scale separation is
obtained through the filtering operation 1, and the unknowns T˜ij ,Q̃i correspond to the so-called subgrid scale (SGS)
(c.f. [11, 26]). The unkresolved SGS stress tensor T˜ij requires a sub-grid turbulence model. Introducing the concept
of SGS turbulent viscosity most models read [30]:
1 ˜
T˜ij = (ug ˜
i uj − ũi u˜j ) = −2νt Sij + Tll δij (4)
3
with
1 ∂ ũi ∂ u˜j 2 ∂ u˜k
S˜ij = + − δij (5)
2 ∂xj ∂xi 3 ∂xk
In equations 4 and 5, S˜ij is the resolved strain tensor and νt the SGS turbulent [Link] analysis yields
√ √
νt ∝ lSGS × qSGS , where lSGS is the length scale of the unresolved motion and qSGS its velocity scale. In the
Smagorinsky model [30] used here, the expression of νt follows:
∂ T̃
Q̃i = −κt (7)
∂xi
Note that T̃ is the modified filtered temperature and satisfies the modeled filtered state equation p̄ = ρ̄RT̃ .
The description does not intend to be comprehensive and for further information the reader is referred to [26]. The
numerical implementation of LES in the computer code AVBP is also given. Further information about AVBP is
found in [27].
The LES code (AVBP) solves the LES transport equations on structured, unstructured or hybrids grids (cf:
[Link] The numerical approach is based on finite volumes schemes using the cell-vertex
method [27] and offers third-order spatial and temporal accuracies. Variations in the filter sizes due to non-uniform
meshes are not directly accounted for in the LES models. Changes in the cell topologies and sizes are only accounted
1/3
for through the use of the local cell volume, i.e. ∆ = Vcell . Grid refinement needs therefore to be carefully controlled
for the LES model to operate efficiently. Such effects are beyond of the scope this work although great care has been
taken to minimize the consequences on the predictions.
Concerning artificial viscosity aspects, the Jameson formulation, with the second order and fourth order coefficients
set respectively to 1/40 and 1/256 is used, triggered on by sensors on strong gradients. The Von Neumann-Richtmyer
viscosity [33] is also used to handle shocks. The right hand side of the momentum transport equation receives the
term ∂pVij N R /∂xi . It performs a local perturbation of the pressure field triggered by the velocity divergence, which
thickens shocks. This tensor is expressed as
∂uk ∂ρui
pVij N R = −∆2 max
δkl δij (8)
k=1,n ∂xl ∂xj
A dimension analysis of this term shows a negligible effect on acoustic perturbations. Indeed, using subscript 0 for
reference state, subscript ac for quantities related to acoustics and F for the frequency of acoustic perturbations, the
Von Neumann-Richtmyer term related to acoustics sizes:
2
ρ0 ∆2 ∂u
V N R
pij ac
≈ ∂x iac
h 2
∆2 ∂p
≈ ρ0 c20 ∂x
2
ac (9)
∆Pac ∆2 Fac 2
≈ ρ0 c20 2 c20
(ρ0 c20 )
I II III
5.6. Simulations 113
Term I is the size of Euler terms in Eq. 2.b. For acoustic features, term II is negligible. Hence, the VNR term is
negligible as long as term III is lower than unity, i.e. F < cO /∆. This last condition is always verified, since the
maximal frequency reachable in an explicit code is c0 /∆.
C. Inlet conditions
The inlet is a Dirichlet supersonic inlet. Three versions are tested. Case 1 is a laminar supersonic inlet. In Case 2, a
random noise is added to the velocity at each node. The amplitude of the perturbation is about 5%. In Case 3, several
coherent perturbations are added on the velocity field, following the approach of Smirnov [31]. The strength and scale
are set assuming an homogeneous isotropic turbulence spectrum. The local amplitude of the coherent perturbations
follows approximately the fully developed RMS velocity profile of a turbulent duct. A constant pressure outlet is set
downstream with an NSCBC weak formulation allowing small inflows.
The serie of cases investigated begins with the jet described in Chuech [5](Mi = 1; N P R = 1.2). Cases 1,2 and 3 are
all related to this jet, differing by their respectives injection conditions. Case 1 injects a sonic under-expanded laminar,
i.e. without any temporal perturbation of the mean profile. In Case 2, the injected velocity receives a perturbation
~ x, t) of random magnitude and direction. The random step is estimated at each iteration for each node and shows
δv(~
a flat power spectral density. The perturbation is then considered as a white noise spatially and temporally non-
correlated. Case 3 uses the turbulence injection previously presented, with the RMS profile on velocity inspired from
fully developed flows (FDF) : 5% to 7% from the axis to the rim, with a circumferential layer, 1/10D thick , reaching
15% at the rim. Case 4 is the jet studied by Seiner and Norum [28](Mi = 1.99; N P R = 1.47), with the same turbulent
injection.
III. VALIDATIONS
A. Inviscid phenomenas
The diamond shock pattern is a repetition of expansion and compression zones. In the comparison of two similar
jets, the small differences on these compressions/expansions accumulate over the downstream distance. Hence the
longitudinal pressure profile is a rough test for the validation of simulations.
The experiment of Seiner gave the longitudinal pressure profile of a Mach 1.88 weakly underexpanded jet (N P R =
1.47). These conditions give a long diamond shock pattern with ten compression/expansion cells visible down to
x = 25D. The simulation of this jet is compared to this experimental data. To set up the plot of Fig. 3, a hundred
of probes are distributed on the jet axis, for 0 < x < 20D. No interpolation is used. The mean pressure profile is
displayed, with error bars giving the minimum and the maximum pressure recorded at each location.
The cell length observed is 2.5D in the simulation, to be compared with 2.6D in the experiment. The oscillations
vanishes at x = 20D in the simulation, and x = 25D in the experiment, while the difference of extremas is growing
downstream, reaching the amplitude of the mean signal oscillations around x = 15D. One can note that the oscillations
of minimas are slightly decayed downstream.
This LES proved a fair agreement between experiments and simulations concerning the inviscid structure. However,
the simulation overestimates dissipation. The fluctuating nature of the flow must now be investigated.
2.5
Experiment
Simulation
2
1.5
P/P0
0.5
0
0 10 20
x/D
FIG. 3: Seiner jet and Case 4 : Mach 1.99 and pressure ratio 1.47 (verifer si resultats meilleurs avec turb FDF)
B. Viscous phenomenas
Turbulence is visible through the development of the jet plume. The fluctuations can be measured at the first order
with averaged quantities of velocity ū, and at the second order with the Root Mean Squared velocity uRM S = u¯2 − ū2 .
Chuech held experimental measurements of this nature on a sonic weakly underexpanded jet (N P R = 1.2), both on
longitudinal and radial locations. The simulation of this jet is profiled the same way with a hundred of probes in
the jet axis for x < 20D and four orthogonal rakes of ten probes. The unstructured grid imposed slightly staggered
probes, but no interpolation was necessary. Mean and fluctuating quantities has been averaged according to the
Reynolds definition, from 8.25ms to 10.5ms. This time must be compared to the convective time, i.e 0.65ms for a
particle travelling the twenty first diameters.
5.6. Simulations 115
The comparison between simulation and to Chuech’s data is displayed in Fig. 4 for longitudinal results and in
Fig. 5 for radial results. Fig. 4.a points out the longitudinal segregation between the potential core (0 to 5D) and
the turbulent decay of the plume (7D and downstream). Simulated velocity profiles oscillate in the potential core,
as expected by theory. Experimental measurements are steady but give a spatially averaged profile where the global
increase of the velocity is visible. Concerning the different injections, the transition zone of Case 3 is located around
5D, like in the experiments, at least 3D upstream of the transition zones found for Cases 1 and 2. In Fig. 4.b, the
RMS values are normalised by the centerline velocity. Case 1 begins with a low fluctuating level around 4%, kept
until 10D. The fluctuations jumps up to 30%. Case 2 start with 15% of RMS velocity, but this is quickly damped
down to 3%. Then it rises up around 8D, reaching the level of 25%. Case 3 starts with 5% of RMS fluctuations, and
keep this level until 4 − 5D where the profile rises up to 20%. Longitudinal profiles shows then that the turbulence
injection develops a correct jet concerning the fluctuations level. The transition zone is slightly underestimated.
Fig. 5 shows the speading of the jet for the three cases. The first plane x/D = 0.2 confirm the accuracy of the
injected profile on mean velocity. Fluctuations are negligible for Case 1, and lower than 3% for Case 2, while the
imposed level was 15% 0.2D upstream. Case 3 shows the correct level of fluctuations in the core of the jet. High
magnitude of RMS in the turbulent layer are damped down to 10%. Profiles of mean velocity for the three planes
x/D = 5, 10, 20 follow the autosimilar classic behavior of free jets.
These results show again a good agreement between simulations and experiments concerning the viscous part of the
flow. The injection of coherent structure, i.e. Case 3, shows the best performances, especially at the beginning of the
jet. The two other injections underestimates the spreading of the jet. The differences fade out with the development
of the jet. The hierarchy of performances for the three injections is conform to intuition: Turbulent like injection is the
best choice, followed by the white noise injection. The inner processes leading of these differences may be investigated
with the analysis of the exhaustive data produced by these numerical simulations.
0.80
1.4 Exp. Exp.
Flat. Flat.
Noise. Noise.
1.2 Turb. Turb.
0.60
1.0
0.8
U/Ue
U/Uc
0.40
0.6
0.4 0.20
0.2
0.0 0.00
0.0 5.0 10.0 15.0 20.0 0.0 5.0 10.0 15.0 20.0
r/D r/D
a) Mean b) RMS
FIG. 4: Mean and RMS values of longitudinal velocity for Cases 1,2,and 3 compared to experiments
116 Chapitre 5. Jets supersoniques sous-détendus
0.40
1.00
Exp. Exp.
Flat. Flat.
Noise. Noise.
Turb. Turb.
U/Uc
U/Uc
0.20
0.50
0.00 0.00
0.0 1.0 0.00 1.00 2.00
r/D r/D
x/D = 0.2D
1.00 0.40
Exp. Exp.
Flat. Flat.
0.80 Noise. Noise.
Turb. 0.30 Turb.
0.60
U/Uc
U/Uc
0.20
0.40
0.10
0.20
0.00 0.00
0.00 0.10 0.20 0.30 0.40 0.00 0.10 0.20 0.30 0.40
r/D r/D
x/D = 5D
1.00 0.40
Exp. Exp.
Flat. Flat.
0.80 Noise. Noise.
Turb. 0.30 Turb.
0.60
U/Uc
U/Uc
0.20
0.40
0.10
0.20
0.00 0.00
0.00 0.10 0.20 0.30 0.40 0.00 0.10 0.20 0.30 0.40
r/D r/D
x/D = 10D
1.00 0.40
Exp. Exp.
Flat. Flat.
0.80 Noise. Noise.
Turb. 0.30 Turb.
0.60
U/Uc
U/Uc
0.20
0.40
0.10
0.20
0.00 0.00
0.00 0.10 0.20 0.30 0.40 0.00 0.10 0.20 0.30 0.40
r/D r/D
x/D = 20D
Mean RMS
FIG. 5: Mean and RMS values of velocity for Cases 1,2,and 3 compared to experiments
5.6. Simulations 117
IV. RESULTS
The diamond shock pattern is not steady, as observed in Fig. 3. This unsteadiness must be investigated from a
temporal point of view. Previous numerical investigation has been focused on shock motion, in order to study the
receptivity of a sreeching jet [17]. The simulated jets are not screeching, but the same mechanisms are present.
The jet of Seiner (Mach=2, NPR=1.47) is used here, with the pressure records of a hundred of probes on the jet
axis. The local pressure is mapped versus time in the range [9, 11.5ms] and centerline location for x < 20D on Fig. 6.
The averaged local velocity ū(x) and the speed of sound c are computed from the dataset, and are represented with
lines on the pressure map. Sound speed slopes lines (white) going upstream and downstream are placed on each
alignment of three pressure minima.
For x < 7D, the vertical stripes of the pressure map show three quasi-steady compression expansion cells. For
7D < x < 15, the cells are more and more staggered. Beyond x = 15D, coherent cells are not visible anymore. In the
staggered region, several peaks are visible both for both maximas and minimas. These peaks travel downstream at
the sonic speed, according to the reference slope +c. The locations of these peaks are also aligned with the reference
slope of upstream sonic speed −c. The reference slope u proves that these peaks are definitely not advected.
Fig. 6 shows that the peaks observed on pressure extremas are at the crossroads of structures travelling downstream
and upstream at the sound speed. The advected coherent signal does not affect directly the diamond shock structure,
a more complex process is in action.
10
The fluctuations observed on Fig. 6 were measured on the centerline, but strong unsteady processes also take place
in the mixing [Link] RMS value of the density gradient norm is used as an acoustic activity indicator on Fig. 7,
~ RM S on longitudinal planes for the three types in injection used in the simulation of the Seiner’s Jet
showing |∇ρ|
(Mach 1, N P R = 1.2). These planes are direct cuts, i.e. without azimuthal averaging.
The spreading of the jet increases from Case 1 to Case 3: '5◦ for laminar injection, '7◦ for white noise injection,
'9◦ for turbulent-like injection. The acoustic activity zone is concentrated closer to the injection in Case 3, compared
to Case 1 and 2. Close to the injection, Case 1 is very quiet, excepted on the diamond shock structure. Narrow zones
of activity are located on the end of each expansion/compression cell for all simulations. The zoomed views displayed
on the left hand side show that peaks are higher but lose their strenght quicker for Case 3, compared to Cases 1 and
2. A regular and intense activity is visible in the mixing layer of Case 3 (x < 7D), with peaks at the crossing of
compression zone with the mixing layer. On the opposite, the mixing layer of Case 1 and 2 are quiet with transient
high peaks on the same spots.
To summarize, each of the three cases simulated has its own scenario of development. Compared to Case 2, Case
1 starts with a quiet cell and develops slower. On the opposite, Case 3 simulation is quickly developped. Compared
to centerline shock, the mixing layer is almost inactive for Cases 1 and 2, while it is preponderant for Case 3.
Case 1
Case 2
Case 3
~
FIG. 7: Acoustic sources indicated by Density gradient fluctuations ∇ρ . Left: Near injection close up (3D × 3D), Right:
RM S
overall aspect (20D × 6D)
5.6. Simulations 119
11
In order to identify the three different scenarii, it is usefull to evaluate the energy distribution between the scales,
and its evolution over the jet development. A temporal Fourier transform gives this information. Velocity records of
forty probes in the mixing layer (x < 10D) are available for a temporal range long enough to perform a Fast Fourier
Transformation (16384 samples). The frequency resolution is 465Hz. The left hand side of Fig. 8 is a mapping of
the FFT Power Spectral Density (PSD) versus frequency (abcissa) and probe location (ordinates). The left hand side
plots are the FFT results for five equidistant positions x/D = 0, 2.5, 5, 7.5, 10. Scale is logarithmic for Frequencies
and PSD. A running average over five consecutive points has been applied to clean up the spectrums.
The spectral map of Case 1 shows the progessive filling of the spectrum from low frequencies. From an empty
energy spectra, the successive FFT grow up smoothly to a developed spectrum. The only exception is the brutal
growth of the 20000 Hz mode, at x/D = 6, followed by a broad 16000 Hz mode, at x/D = 8. The main part of the
growth occurs in the zone 0 < x < 2.5D.
The white noise of Case 2 is visible on the first profile (x/D = 0). High frequencies are quickly damped while low
frequencies are growing up. For x/D = 7, the flow reaches a developed spectrum. The map traduces the apparition and
disparition of several peaks, roughly corresponding to 15000, 21000, 28000 and 34000 Hz for the range 0.5D < x < 3D.
Then, a complex distribution of energy is rising beyond x > 3.5D for the frequency range [0,30000Hz]. The main part
of the evolution occurs again in the zone 0 < x < 2.5D.
The turbulent-like injection of Case 3 leads to a initial spectrum already developed at high frequencies. The energy
of low frequencies [0, 400000Hz] at x = 0D is as low as in Case 2, excepted for the modes 10000 and 14000 Hz. The
other visible modes are 57000, 120000, 180000 and 340000 Hz. The spectral map shows a continuous evolution from
the first spectrum to the developed one, the energy concentrating toward low frequencies. Developed spectrum is
almost reached at x/D = 4.
Fig. 8 is giving information about the development of each case. The simulation of Case 1 fills up the energy spectrum
with the lowest frequencies, i.e. the largest structures and the mean flow. Case 2 proved a selection of several modes
by the white noise in the very first diameters, followed by the growth of low frequency modes [0, 30000Hz]. Case 3
develops straightly the same low frequency growth [0, 30000Hz] closer to the injection. All cases lead to developed
spectra similar in shape and amplitude, but the three scenarii are very different from a spectral point of view. The
most discriminant behavior is probably visible in the first diameters, where the variations of energy distribution are
the highest for all cases.
120 Chapitre 5. Jets supersoniques sous-détendus
12
FIG. 8: Power Spectral Density along the mixing layer (FFT estimation)
5.6. Simulations 121
13
Fig. 8 indicated that a spectral analysis focused between x/D = 1 and x/D = 2 would allow to investigate further
the differences between the three screnarii of jet development. As this is a restriction to a small region, simple energy
transfers might be reachable. Bicoherence analysis is designed to identify the energy transfers between the frequencies.
The quantity Bic2F GH is defined in Eq. 10 where F̂ , Ĝ, Ĥ are the Fourier transforms of the velocity signal recorded
at the location F, G, H, a longitudinal triplet of probes separated by a tenth of diameter. In this equation ∗ is used
for the conjugate form, and hi is the averaging on 32 triplets spread in the mixing layer. Quantity Bic2F GH is then
mapped in the Fig. 9 for Case 1,2 and 3. A peak in the right half (F 1 > 0) means that modes F 1 and F 2 feeds the
mode F 1 + F 2. A peak in the left half (F 1 < 0) means that modes F 1 and F 2 feeds the mode F 2 − F 1. All peaks
over the thresold 0.35 are circled.
D E2
F̂ (f1 )Ĝ(f2 )Ĥ ∗ (f1 + f2 )
2
BicF GH (f1 , f2 ) = 2 2 (10)
F̂ (f1 )Ĝ(f2 ) Ĥ(f1 + f2 )
To ease the understanding of this tool, the initial x/D = 1 and final x/D = 2 spectra are displayed on the right hand
side of the figure.
Concerning Case 1, the vertical alignment of correlations (F 1 > 2000Hz, F 2 ' 0Hz) means that the mean flow
feeds all the range of frequencies. The diagonal (F 2 − F 1 ' 0) means that all the frequencies feeds the lowest one,
i.e. the extreme left of the spectrum.
Case 2 gives several diagonal alignments of correlations (F 1+F 2 ' 15000, 21000, 28000, 34000, 40000, etc...) meaning
that all F 1 frequencies feed several modes, as a white noise is expected to do. The same analysis on higher frequencies
(F 1, F 2 > 100000Hz) showed no bicoherence, indicating that the high frequencies injected are simply damped, not
redirected to low frequencies.
Case 3 contained an intense 10000Hz mode in the near injection region. The peak 0-10000Hz on the Fig. 9 proves
that this mode is still fed by low frequencies, one diameter after injection. A broad concentration of significant
10000-10000Hz interactions are also revealed by the bicoherence analysis, feeding the 20000 Hz mode. Then, a last
concentration of peaks in the region [−6000 < F 1 < 6000, 20000 < F 2 < 26000] explains the broad energy stripe
in the range [20000, 26000Hz] by a two-ways transfer from 20000Hz to 26000Hz (F 1 > 0 peaks), and back (F 1 < 0
peaks). No significative correlation can be seen beyond 30000Hz.
Fig. 9 enhances the contrast between the three scenarii. The laminar injection (Case 1) leads to a development
feeded by the low frequencies, i.e. the largest structures. The white noise injection excites several modes in the jet.
The evolution of spectrum, in Fig. 8 shows that these singular modes are quickly replaced by the developed spectrum.
The turbulent-like injection excites only the [0, 30000Hz] range, leading sooner than Case 2 to the developed spectrum.
122 Chapitre 5. Jets supersoniques sous-détendus
14
Case 1
7e+08
x/D=2
x/D=1
6e+08
5e+08
4e+08
PSD(U)
3e+08
2e+08
1e+08
0
0.0 10000.0 20000.0 30000.0 40000.0 50000.0
Frequency [Hz]
Case 2
5e+09
x/D=2
21000Hz
x/D=1
4e+09
15000Hz
3e+09
34000Hz
PSD(U)
28000Hz
40000Hz
2e+09
1e+09
0
0.0 10000.0 20000.0 30000.0 40000.0 50000.0
Frequency [Hz]
Case 3
1.4e+10 x/D=2
x/D=1
10000Hz
1.2e+10
1e+10
18000Hz
PSD(U)
8e+09
6e+09
20000Hz
26000Hz
4e+09
2e+09
0
0.0 10000.0 20000.0 30000.0 40000.0 50000.0
Frequency [Hz]
15
The previous analysis enlighted the development of the jet downstream. The temporal evolution of the pressure
waves along the axis (Fig. 6) revealed that both upstream and downstream sonic velocity carried the fluctuations.
The frequency-wavenumber analysis allows to track the advection velocity of coherent structure, giving the energy
balance between upstream and downstream directions. This is particularly interesting in the near injection region,
where the receptivity of the jet is important.
Nine rakes of ten probes are placed longitudinally at [D < x < 2D], 2r/D = 0; ±0.3, ; ±0.4; ±0.5; ±0.6 in the
simulations of Cases 1,2 and 3. Each rake gives the necessary data to perform a Frequency-Wavenumber Spectrum
Analysis (FWSA), using the interspectra between all the probes, as introduced by Capon [3]. The distance between
probes being ∆x, the phase decay ∆Φ(f ) observed between two signals of frequency f is interpreted in the wave-
number form λ(f ) = ∆φ(f )/∆x = 2π∆φ(f )/∆t. A peak on the FWSA map, located at the position (f, λ) indicates
that a coherent signal of frequency f is travelling trough the probes at the speed 2πf /λ. Fig. 10 gives on the right
hand side the results at the center line, referring to the central rake (r = 0), and on the left hand side the mixing
layer results, referring to the other rakes. Concerning the mixing layer, the data recorded at the outer side of the jet
gave the same FWSA results as the inner side data. Hence, all mixing layer results are averaged. As all the spectra
are fading out with the high frequencies, the FWSA coherence must be normalized. Fig. 10.a and .b gives the norm
nF W SA (f ) used for each case, taken as the coherence of frequency f averaged over all the wavenumbers.
The convective speed expected inside the sonic jet is +c. The related acoustic speeds are then 0 and +2c. How-
ever acoustic waves propagate also in the still air surrounding the jet, leading to the acoustic speeds −c and +c.
Consequently, FWSA may capture the velocities −c (external upstream acoustic), +c (external downstream acoustic,
internal convection) and +2c (internal downstream acoustic)
First, the norm nF W SA (f ) profiles confirm the coherence fading out for high frequencies in both centerline and
mixing layer region, in particular for Case 1 (Fig. 10.a,.b). The laminar injection of Case 1 has shown a spectrum
limited to low frequencies. FWSA finds out these frequencies travelling downstream (+c) in the mixing layer (f <
40000Hz). A large range of structures is also identified, that travels upstream (−c) in both regions. Due to the
normalization, the upstream of waves for f < 40000Hz is hidden by the downstream traveling part. White noise
injection of Case 2 and Turbulent noise of Case 3 show only downstream travelling structures. Their velocity matches
(+c) excepted in the mixing layer of Case 2, where the velocities are observed between +0.6c and +c.
The large structures created in Case 1, and the injected structures in Case 2 and 3 are all captured by FWSA and
convected by the jet at the velocity +c. However, the laminar injection being very quiet, some upstream travelling
acoustic structures are observed, and should be linked to the feedback process occuring in screeching jets.
124 Chapitre 5. Jets supersoniques sous-détendus
16
22 22
10 10
20 20
10 10
18 18
10 10
16 16
10 10
14 14
10 10
CSST Norm
CSST Norm
12 12
10 10
10 10
10 10
8 8
10 10
6 6
10 Case 1 10 Case 1
Case 2 Case 2
4 4
10 Case 3 10 Case 3
2 2
10 10
0 0
10 10
0.0 20000.0 40000.0 60000.0 80000.0 100000.0 0.0 20000.0 40000.0 60000.0 80000.0 100000.0
Frequency [Hz] Frequency [Hz]
FIG. 10: Coherent Structure Spectral Tracking applied in the region 1 < x/D < 2 at the centerline and at the mixing layer.
5.6. Simulations 125
17
[1] Abdol-Hamid, K. S., and R. G. Wilmoth, 1989, AIAA Journal 27(3), 315.
[2] Bogey, C., C. Bailly, and D. Juvé, 2003, Theoretical and Computational Fluid Dynamics 16(4), 273.
[3] Capon, J., 1969, in IEEE, volume 57, pp. 1408–1418.
[4] Chang, I.-S., 1945, Mach Reflection, Mach Disc, and the Associated Nozzle Free Jet Flows, Ph.D. thesis, University of
Illinois at Urbana-Champaign.
[5] Chuech, S. G., M. C. Lai, and G. M. Faeth, 1989, AIAA Journal 27, 549.
[6] Crist, S., D. R. Glass, and P. M. Sherman, 1966, AIAA Journal 4, 68.
[7] Cumber, P., M. Fairweather, S. Falle, and J. Giddings, 1995, Journal of Fluids Engineering 117, 599.
[8] Dash, S. M., and D. E. Wolf, 1985, AIAA Journal 23(4), 505.
[9] Eastman, D., and P. Radtke, 1963, AIAA Journal 1, 918.
[10] Favre, A., 1969, in Problems of hydrodynamics and continuum mechanics (SIAM, Philadelphia), pp. 231–266.
[11] Ferziger, J., 1997, in New tools in turbulence modelling, edited by O. Métais and J. Ferziger (Les Editions de Physique -
Springer Verlag), pp. 29 – 47.
[12] Imamoglu, B., and P. Balakumar, 2002, AIAA Paper 2002-2527.
[13] Imamoglu, B., and P. Balakumar, 2003, in 9th AIAA/CEAS Aeroacoustics Conference and Exhibit; Hilton Head, SC.
[14] Jorgenson, P., and C. Loh, 2002, AIAA Paper 2002-3889.
[15] Lehnasch, G., 2005, Contribution à l’étude numérique des jets supersoniques sous-détendus, Ph.D. thesis, Université de
Poitiers.
[16] Lewis, C., and D. Carlson, 1964, AIAA Journal 2(4), 776.
[17] Li, X., and J. Gao, 2005, Phys. Fluids 17(085105).
[18] Loh, C., H. S.C., and P. Jorgenson, 2001, AIAA Paper 2001-2252.
[19] Manning, T., 2000, A numerical investigation of sound generation in supersonic jet screech, Ph.D. thesis, Standford
University.
[20] Manning, T., and L. S. K., 1998, AIAA paper 1998-0282.
[21] Norum, T. D., 1982, in AIAA, Aerospace Sciences Meeting and Exhibit, 20th, Orlando.
[22] Palmer, J., and R. Hanson, 1998, AIAA Journal 36(2), 193.
[23] Poinsot, T., and D. Veynante, 2001, Theoretical and numerical combustion (R.T. Edwards).
[24] Powell, A., 1953, Aeronaut. Quart. 4(103-122).
[25] Raman, G., 1998, Prog. Aerosp. Sci. 34(1), 45.
[26] Sagaut, P., 2000, Large Eddy Simulation for incompressible flows, Scientific computation series (Springer-Verlag).
[27] Schönfeld, T., and M. Rudgyard, 1999, AIAA Journal 37(11), 1378.
[28] Seiner, J. M., and T. D. Norum, 1980, AIAA Paper 1980-965.
[29] Shen, H., and C. K. W. Tam, 1998, AIAA Journal 36(10), 1801.
[30] Smagorinsky, J., 1963, Mon. Weather Rev. 91, 99.
[31] Smirnov, A., S. Shi, and I. Celik, 2001, J. Fluid Eng. 123(2), 359.
[32] Thies, A., and C. K. W. Tam, 1996, AIAA Journal 34(2), 309.
[33] VonNeumann, J., and R. Richtmyer, 1950, Journal of Applied Physics 21(3), 232.
[34] Zapryagaev, V., V. Pickalov, N. Kiselev, and A. Nepomnyashchiy, 2004, Theor. Comput. Fluid Dynamics 18(2-4), 301.
[35] Zapryagaev, V. I., N. P. Kiselev, and A. A. Pavlov, 2004, J. of Applied Mechanics and Technical Physics(Russia) 45(3),
335.
126 Chapitre 5. Jets supersoniques sous-détendus
1
Laboratoire de Simulation Numérique en Mécanique des Fluides
2
Ecole Nationale Supérieure d’Art et Metiers
5.6. Simulations 127
T. 5.2 – Résultats d’une l’analyse de stabilité linéaire effectué sur un champ de base obtenu par
SGE du jet de S.G. Chuech (29).
5e+09
x/D=2
21000Hz
x/D=1
4e+09
15000Hz
3e+09
34000Hz
PSD(U)
28000Hz
40000Hz
2e+09
1e+09
0
0.0 10000.0 20000.0 30000.0 40000.0 50000.0
Frequency [Hz]
F. 5.18 – Spectres basses fréquences d’un jet sous-détendu excité par un bruit blanc. Coupes à un
et deux diamètres
128 Chapitre 5. Jets supersoniques sous-détendus
o
66
o
33
F. 5.19 – Vitesse longitudinale observée sur le cylindre r/D = 0.55, et angle des premières struc-
tures turbulente obtenues par analyse linéaire de l’instabilité.
5.7. Elements particuliers au calcul de jets supersoniques 129
Domaine de calcul
Zone Tampon
a) Ecoulement externe
Domaine de calcul
Zone Tampon
b) Ecoulement interne
Le passage du domaine de calcul à la zone tampon exige une variation faible et continue des
130 Chapitre 5. Jets supersoniques sous-détendus
mailles. Dans le cas contraire, la zone tampon crée des perturbations parasites dans la solution : elle
est ”bruyante”. Cette contrainte se traduit par un ratio maximal entre le volume de deux cellules
adjacentes. L’ordre de grandeur du ratio volumique est de 2 pour un maillage de bonne qualité
, soit un ratio de 1.25 sur les cotés des cellules. Les zones tampons utilisent conjointement une
augmentation de la taille des cellules avec une augmentation de la viscosité artificielle. On prendra
soin de placer le gradient de viscosité à l’aval dans une zone à pas de maillage constant, comme sur
la Fig. 5.21. En effet, l’addition des deux variations (maille et viscosité) peut rendre la simulation
bruyante localement. Avec une zone tampon qui respecte les régles énoncées précédemment, le
gain en robustesse est substantiel. Par contre, le faible ratio entre les cellules peut se révéler coûteux
en terme de maillage. Cette technique est utilisée dans tous les travaux présentés dans cette thèse
faisant intervenir un jet supersonique.
Periodique
Domaine de calcul
F. 5.22 – Domaine de calcul d’une chambre de combustion suivie d’un distributeur transsonique
Periodique
Domaine de calcul
F. 5.23 – Domaine de calcul corrigé d’une chambre de combustion suivie d’un distributeur trans-
sonique
132 Chapitre 5. Jets supersoniques sous-détendus
”Dear Antoine,
(...)About your question : ”In other words, does an adapted non parallel jet can develop diamond
shock structures ?” I think it is possible. Term ”adapted” is characteristics of one from some
gasdynamic parameters of supersonic jet in local region (nozzle exit only). Therefore if only one
from jet parameters is adapted to ambient condition (P jet = Pambient ) but angle of velocity vector
α is nonparallel of axis you have diamond structure. If you have P jet = Pambient and α = 0, the
jet (without viscosity) is a parallel flow without any structures. For your case, how I understand,
for Mexit = 1 any disturbance (P jet , Pambient or α , 0 i.e boundary is nonparallel) is a cause
for appearance of diamond structure. You can see a likewise phenomena - pseudoshock which
corresponds to transform supersonic flow to subsonic flow in channel. (...)
Best regards
Valery Zapryagaev”
Dr. Valery Ivanovich Zapryagaev, Lecturer, Russian Academy of Sciences Siberian Division ;
Institutskaya Str. 4/1, Novosibirsk 630090, Russia
Il faut donc noter qu’un jet supersonique dont la pression est adaptée sur la section de sortie peut
présenter des structures de compression/expansion si il n’est pas strictement parallèle à l’injection.
Il est fort probable que les expérimentateurs aient été les premiers spectateurs du paradoxe jet
adapté / écoulement en diamants. En effet, un jet fortement supersonique (Mach = 2 par exemple)
doit être accéléré au moyen d’un écoulement divergent. Cependant, il faut réaliser une augmenta-
tion de section tout en minimisant l’encombrement du dispositif. La comparaison de deux injecteurs
coaxiaux montre l’évolution de la solution technologique. La Fig.5.24.a montre une coupe de l’in-
jecteur utilisé par Cheng et Wehrmeyer (28) en 1994. Dans cette géométrie, le fluide subit une diver-
gence centrifuge. Au contraire, dans l’injecteur choisi par Cutler (41) en 2001 illustré en Fig. 5.24.b,
la divergence est centripète. Les lignes de courant extérieures sont rigoureusement parallèles à l’axe
du jet. La géométrie choisie par Cutler permet d’obtenir un jet fortement supersonique en concen-
trant les perturbations de parallélisme à l’intérieur. Les expérimentateurs ont adapté leur installations
pour réaliser un écoulement peu perturbé par les phénomènes supersoniques. Afin de suivre au plus
près les résultats expérimentaux, les simulations numériques doivent également prendre en compte
la composante radiale des vitesses imposées par les géométries de chaque injecteur.
5.7. Elements particuliers au calcul de jets supersoniques 133
a) 1994 - Jet coaxial de Cheng (28) b) 2001 - Jet coaxial de Cutler (41)
F. 5.24 – Deux types de divergents montés sur des jets coaxiaux
134 Chapitre 5. Jets supersoniques sous-détendus
Les bords d’un injecteur supersonique posent un problème numérique : la ligne sonique est
décalée du bord sur certaines simulations. La Fig. 5.25 montre la solution numérique donnée par
deux codes. Les deux utilisent des solveurs de Roe. Le premier est un code de recherche développé
lors de la thèse de D. Chargy (22). Le second est un code commercial bien connu, FLUENT 6.0 . Les
contours de Mach permettent d’indentifier clairement une structure de jet fortement sous-détendu,
avec un disque de Mach nettement marqué. On observe dans les deux cas un décalage très net de la
ligne sonique marquant les bords du jet sur le plan d’injection. Cette figure met donc en évidence un
comportement imprévu des simulations par différences, éléments, ou volume finis, comparé à une
résolution par la méthode des caractéristiques qui suppose une expansion centrée sur le coin. Il faut
donc analyser le décalage observé dans ces simulations.
F. 5.25 – Contournement obvervé dans des simulations de jets sous-détendus, Contours du nombre
de Mach
Cependant, il est étonnamment difficile de prouver théroriquement que cette ligne de glissement
doit correspondre au bord :
5.7. Elements particuliers au calcul de jets supersoniques 135
Dear Antoine, (...) Concerning your final point about the slip line, I do not see why it could not
be possible (to observe a decayed slip line), but we never looked there experimentally. Sincerely
Godfrey
Godfrey Mungal, Professor Mechanical Engineering Department Stanford University, Stanford,
CA 94305-3032
(...)Une chose est sure : il est numériquement trés difficile de calculer les écoulements d’arrière
corps. Il y a des effets numériques que l’on contrôle difficilement
Le modèle de Prandtl Meyer est effectivement un modèle idéal avec attachement du choc sur
l’arète. Et le second modèle(ligne décalée) ne me choque pas, Numériquement on observe des
vitesses qui glissent le long de la paroi : la recirculation ”n’arrive pas au coin” , ce qui justifie
la solution obtenue. Je ne sais pas dire si cela est physique ou non : La nature ne connait pas les
”coins” !
Merci de me tenir informer car le sujet est intéressant.
Dr. Didier Chargy INCKA-Simulog - Groupe Robinson Agence de Cannes
Lorsque l’on observe l’écoulement réel près d’une injection supersonique, il existe une échelle
en dessous de laquelle le bord de l’injection ne peut plus être considéré comme géométriquement
parfait, c’est à dire formant un angle droit. Prenons un écoulement de coin contenant une expansion
de Prandtl-Meyer centrée sur ce coin, situation illustrée par la Fig. 5.26.a. A une échelle beaucoup
plus petite, c.a.d. en dessous de la précision d’usinage de l’injecteur, ce coin ne forme plus un
angle droit parfait. Supposons que la forme de ce bord d’injection soit arrondie, comme illustré en
Fig. 5.26.b. L’expansion de Prandtl Meyer existe toujours, mais la ligne Mach = 1 est décalée vers
l’extérieur. Le fluide rampe le long de la paroi, s’accélère puis se recomprime rapidement, formant
une petite ligne de choc contre le plan d’injection, qui se transforme en ligne de glissement plus loin
en aval. Ainsi, la ligne sonique marquant la frontière du jet est légerement décalée vers l’extérieur
lorsque l’on regarde suffisamment près. Cette observation mène à la reflexion suivante : L’expansion
de Prandtl-Meyer centrée sur un coin orthogonal est une approximation qui permet de simplifier le
problème de l’écoulement de coin. Cependant, choisir un point géométrique présentant plusieurs
états physiques est une singularité mathématique.
Le décalage que nous avons décrit ne remet pas en cause la présence d’une expansion de Prandtl-
Meyer. C’est grâce à cette propriété que le Dr. Katanoda, professeur à l’université de Kyoto, établit
le test permettant de différencier un décalage physique d’une aberration numérique. Suite à une
conversation par courrier électronique sur ce sujet, il propose simplement de vérifier la conservation
de l’entropie à l’intérieur de l’expansion. Une expansion de Prandtl-Meyer étant rigoureusement
isentropique le long des lignes de courant, si une variation de l’entropie apparaı̂t, la solution est
non physique. Ce test est utilisable sur le champ d’entropie donné en Fig. 5.27 correspondant au
jet illustré en Fig. 5.25.a. La ligne de courant la plus proche de la paroi traverse trois contours
d’entropie au niveau du coin. Le coutournement observé est donc non-physique, d’après le critère
du Dr. Katanoda. Lorsqu’un décalage apparaı̂t dans une solution, le test de l’entropie permet ainsi
de vérifier la validité de la solution à cet endroit.
136 Chapitre 5. Jets supersoniques sous-détendus
M=0
M=0
Pamb
Pamb
M>M 0
M>M 0
Pamb Pamb
M0 M0
P > Pamb P > Pamb
”Dear Antoine
Thank you for an interesting inquiry. I have never seen underexpanded jets flowing ( with a
decayed slip line) in both of my experimental and computational work.
Solution similar to (the decayed configuration) can be seen at around the starting iterations in
computation if you put atmospheric condition as initial one. However, I have never seen ( a
decayed configuration ) in steady state. In steady state, when air-flow of Mach number 1 truns 90
degrees it has to be accelerated to nearly Mach 7.0. I wonder if your computational result is like
that or not. I recommend to check entropy distribution along stream-line around the nozzle lip.
Regards,
Hiroshi Katanoda
Dr. H. Katanoda, Lecturer, Department of Mechanical Systems and Environmental Engineering ;
The University of Kitakyushu
La méthode de l’injection décalée a permis de traiter avec succès un jet sonique présentant un
rapport de pression de 3. Cette simulation a été opérée sur un maillage de 500000 tetraèdres. Le
diamètre du jet vaut 10mm, le coté d’un tétraèdre 0.5mm. Le calcul commence par un jet sonique
adapté dont la pression augmente de 1 bar chaque milliseconde. Les contours de Mach sont illustrés
sur la Fig. 5.29. Il montre clairement un disque de mach intense, et la présence d’une deuxième
cellule de détente et compression. Ce calcul démontre l’efficacité de la stratégie d’injection décalée
pour la simulation de jets fortements sous-détendus.
Dans cette simulation, les variations des quantités dans la couche de mélange sont dues princi-
palement à du bruit numérique. Il est en effet impossible de distinguer une structure nette dans cette
zone, même au moyen d’un critère Q (73) par exemple. Ce défaut implique que cette simulation ne
peut être considérée comme une simulation aux grandes échelles. Malgré cela, le calcul montre un
phénomène intéressant. La figure Fig. 5.30 montre l’isosurface M = 1. Cette surface se plisse en al-
lant vers l’aval, et les plis ainsi formés sont instationaires. Un telle observation prouve que le calcul
ainsi posé développe une vorticité longitudinale, avec des vortex alternés et alignés dans la couche
de mélange parallèlement au jet. Il s’agit du premier pas vers l’instabilité de Görtler, présente dans
la couche de mélange courbe des jets sous-détendus.
En parallèle, des calculs effectués par le Dr. A. Beer en 2005 ont repoussé les limites des simu-
lations. Tout d’abord, un jet sonique fortement sous-détendu avec un rapport de pression NPR = 10
138 Chapitre 5. Jets supersoniques sous-détendus
F. 5.29 – Contours de Mach d’un jet sous-détendu NPR = 3, Mach = 1. Code AVBP.
F. 5.30 – Isosurfaces Mach 1 (blanc) et gradients de densité (gris) d’un jet sous-détendu NPR =
3, Mach = 1. Code AVBP.
est simulé avec un maillage de seulement 250000 nœuds. On peut observer l’allure caractéristique
de ce jet sur la Fig. 5.32, avec une seule cellule de détente terminée par un disque de Mach. Pour les
rapports de pression supérieurs, la couche de cisaillement exige une plus grande densité de points,
et donc un maillage plus lourd.
Ensuite, une configuration plus complexe est étudiée : deux jets fortement sous-détendus (NPR =
5, Mach = 1) entrent en collision avec un angle de 60◦ . La simulation montre une jonction des deux
jets en un seul, plus étroit dans le plan d’impact, qui contient encore des cellules de détente et
5.7. Elements particuliers au calcul de jets supersoniques 139
F. 5.32 – Contours de vitesse longitudinale et ligne sonique d’un jet fortement sous-détendu
NPR = 10, Mach = 1. Avec l’aimable autorisation d’A. Beer. Code AVBP.
compression (Cf. Fig. 5.33). Ce comportement est comparé à celui observé sur les strioscopies de
l’université de Kinki sur des jets presentant un rapport de pression et un angle légèrement différent,
respectivement NPR ≡ 2 et Angle = 30◦ . Le calcul n’ayant pas été réalisé dans ce but, il n’y a pas
de comparaison quantitative avec l’expérience de Kinki.
L’impact de deux jets est généralement réalisé dans le but d’augmenter le mélange entre les deux
jets. La SGE montre avec la Fig. 5.34 l’écoulement complexe instationnaire et le mélange efficace
mais inhomogène créé par cette interaction sur un maillage de 800000 nœuds. La morphologie
pour le moins ’originale’ de la zone d’interêt profite pleinement des avantages d’un maillage non-
structuré. Ainsi, les travaux de A. Beer prouvent que l’approche proposée autorise la simulation de
configurations complexes.
140 Chapitre 5. Jets supersoniques sous-détendus
F. 5.33 – Calcul de l’impact de deux jets fortement sous-détendus : NPR = 5, Mach = 1, Angle =
60◦ . En haut : contours de pression et ligne sonique (noir) de la simulation. En bas strioscopie de
deux jets sous-détendus réalisée à l’université de Kinki. Avec l’aimable autorisation d’A. Beer. Code
AVBP.
5.7. Elements particuliers au calcul de jets supersoniques 141
F. 5.34 – Calcul de l’impact de deux jets fortement sous-détendus : NPR = 5, Mach = 1, Angle =
60◦ . Isosurface de critère Q (73). Avec l’aimable autorisation d’A. Beer. Code AVBP.
142 Chapitre 5. Jets supersoniques sous-détendus
Chapitre 6
Combustion Supersonique
6.1 Objectif
La combustion utilise conjointement les processus de mélange hydrodynamique et de réaction
chimique. Les SGE réactives sont largement appliquées à des systèmes à combustion subsonique où
le mélange hydrodynamique est bien plus lent que la réaction. Le fort rapport d’échelles donne lieu à
différentes approches de la combustion, telles que la flamme épaissie, ou encore le suivi de flamme.
Pour une combustion supersonique, les mélanges hydrodynamiques sont aussi rapides que la com-
bustion. Le travail de simulation se focalise alors sur la compétition entre les deux phénomènes.
La simulation de la combustion supersonique a été relativement peu étudiée. La modélisation
de l’impact direct sur la chimie de la turbulence compressible d’écoulement supersonique n’est pas
encore très développée. La tendance actuelle consiste à augmenter la résolution des simulations tout
en gardant une chimie laminaire. Un article de Jeung (77) illustre bien cette tendance : il étudie la
géométrie de scramjets par une série de simulations Unsteady-RANS au deuxième ordre (Schema
MUSCL) sur des grilles structurées (384000 cellules en 2D) avec un schéma cinétique complet pour
l’hydrogène (8 espèces , 25 réactions).
L’approche proposée dans ce chapitre poursuit cette tendance. La simulation des grandes échelles
implique un maillage fin et tridimensionnel. Une chimie simplifiée compense le coût supplémentaire
en terme de calcul. L’outil final s’applique dans ce chapitre à un jet coaxial stationnaire en situation
d’autoallumage. Outre sa ressemblance avec un allumeur cryotechnique de moteur fusée, elle exige
une chimie simplifiée d’autoallumage, réalisée grâce aux résultats du Chapitre 4. Ce chapitre permet
d’éprouver l’approche actuellement observée dans la bibliographie pour des simulations au grandes
échelles, et identifie les difficultés et points délicats d’une SGE d’allumage d’un moteur fusée.
143
144 Chapitre 6. Combustion Supersonique
Un large jet coaxial chaud (1250◦ K, Mach 2) riche en oxygène apporte l’oxydant et la chaleur. Ce
gaz est chauffé par une combustion préalable pauvre (φ = 0.364). Pour information, en prenant
pour l’hydrogène un pouvoir calorifique inférieur de 120MJ/kg, le préchauffage dégage 207kW. La
flamme externe est en en régime pauvre (φ = 0.123), ce qui assure la consommation complète de
l’hydrogène, dégageant au passage une puissance de 43kW. La flamme extérieure est photographiée
dans le domaine visible sur la Fig. 6.2. Le dégagement de lumière de la flamme se fait au delà de
25 diamètres en aval de l’injection. Un faible dégagement de lumière transitoire se produit autour
de 13 diamètres. Le jet turbulent est visible sur la photographie de type schlieren de la Fig. 6.3.
L’agitation de l’écoulement et l’angle du flash cachent les structures supersoniques internes telles
que le cône de sillage de l’injecteur d’hydrogène.
La configuration retenue regroupe donc les thèmes du jet supersonique coaxial, et de la com-
bustion supersonique en auto-allumage. Driscoll a proposé un dimensionnement de flammes super-
soniques (48) qui donne dans le cas de Cheng une flamme atteignant 300 diamètres. Seule la partie
amont de cette flamme nous intéresse ici, dans les 50 premiers diamètres. La couche de mélange
interne hydogène-oxygène y joue un rôle central. La couche de mélange externe oxygène-extérieur
a un rôle secondaire, car elle peut autoriser un écoulement en diamant. Pour une première étude, sa
simulation détaillée n’est pas nécessaire.
On notera que la flamme a une apparence de flamme détachée. Cependant, le mode de stabili-
sation est l’auto-allumage, non la propagation. Dans ce mode, le contact entre les deux réactifs est
suivi obligatoirement d’une combustion. On ne parlera donc pas de flamme attachée ou détachée
mais de flamme auto-allumée.
Dimensions :
Diamètre int. de l’injecteur carburant 2.36 mm
Diamètre ext. de l’injecteur carburant 3.81 mm
Diamètre int. de l’injecteur comburant 17.78 mm
Injection de comburant chaud :
Alimentation en air (±2%) 0.0735 kg/s
Alimentation en hydrogène (±2%) 0.000173 kg/s
Alimentation en d’oxygène (±3%) 0.0211 kg/s
Pression 107 kPa
Température 1250 K
Mach 2.0
Vitesse 1420 m/s
Fraction molaire d’oxygène 0.201
Fraction molaire d’azote 0.544
Fraction molaire d’eau 0.255
Injection de carburant :
Alimentation hydrogène (±3%) 0.000362 kg/s
Pressure 112 kPa
Température 540 K
Mach 1.0
Vitesse 1780 m/s
Fraction molaire d’hydrogène 1.0
Les Reynolds Re observés sont élevés, donc les grandes échelles de cet écoulement ne sont pas
affectées par la dissipation visqueuse. La turbulence issue du jet d’hydrogène est bien plus faible
que l’oxydant. De plus, la viscosité y est largement surévaluée (+103%) bien qu’elle soit correcte
pour l’oxydant (−3%). Cependant, le débit massique de carburant représente 0.4% du débit total, ce
qui permet de négliger la turbulence émise par ce dernier.
La gamme de Reynolds [15 000, 150 000] est à la portée de la simulation des grandes échelles.
Il convient de préciser les echelles présentes, et de les comparer aux échelles de résolution pour
anticiper la qualité des calculs.
148 Chapitre 6. Combustion Supersonique
Ret (r)
1000
100
nt
ra
bu
ar
Coupure
nt
tc
10
da
Je
xy
to
Je
−6 −5 −4 −3 −2
10 10 10 10 10 r
F. 6.4 – Gamme des échelles couvertes par les jets de carburant et d’oxydant, et choix d’une
échelle de coupure à 1 10−4 m
6.3. Etude dimensionnelle 149
6.3.3 Compressibilité
Le brûleur de Cheng utilisant un jet coaxial supersonique, des effets de compressibilité sont pro-
bablement présents dans la dynamique, induisant une turbulence compressible. Cela pourrait exiger
l’utilisation de modèles de turbulence compressibles. L’importance des d’effets de compressibilité
peut être évaluée par une nouvelle analyse dimensionnelle.
Une couche de cisaillement entre deux fluides 1 et 2 produit de la turbulence compressible si le
cisaillement est suffisamment élevé. Pour l’évaluer, il faut d’abord avoir une vitesse de convection
Uc qui depend des caractéristiques de chaque fluide : vitesse U1,2 , vitesse du son C1,2 , rapport des
chaleurs massiques γ1,2 . Cette vitesse donne pour chaque fluide un Mach convectif Mc,1,2 = Uc /C1,2
indiquant l’impact de la couche limite sur le fluide. La ’pression convective’ Pc doit être équivalente
dans les deux fluides, ce qui mène à l’Eq. 6.1 :
γ1 γ2
γ1 − 1 γ2 − 1
! γ1 −1
! γ2 −1
1− 2
∗ Mc,1 = 1− 2
∗ Mc,2 (6.1)
2 2
Cette relation détermine Uc et permet d’en déduire Mc,1 et Mc,2 . Le Mach convectif indique le degré
de compressibilité pour le fluide. Pour Mc 0.3, la turbulence générée est incompressible. Pour
Mc 0.3, la turbulence générée est compressible.
Pour la couche de mélange carburant/oxydant, le calcul de Uc donne 1527 m/s. Cela implique un
Mach comvectif valant 0.149 pour l’oxydant et 0.143 pour le carburant. Pour la couche de mélange
oxydant/extérieur, le calcul de Uc donne 456.85 m/s. Cela implique un Mach comvective valant 1.35
pour l’oxydant et 1.27 pour l’extérieur.
Ainsi, le mélange entre hydrogène et oxydant chaud qui nous intéresse est soumis à une turbu-
lence faiblement compressible. La turbulence de la couche de mélange extérieure est très fortement
compressible, mais n’a pas d’influence majeure sur les premiers stades de la combustion.
Dans le cadre de notre étude, les modèles de combustion et de turbulence qui ne prennent pas
en compte de turbulence compressible sont donc suffisants.
Oxydant Oxydant
le Tableau 6.1. L’impact de la turbulence sur le temps de réaction se mesure avec le nombre de
Damköhler Da = tt /tc , et vaut ici Da = 2 10− 3. Par conséquent, la chimie est complètement pi-
lotée par la turbulence. Cela se traduit par une très grande hétérogénéı̈té dans l’avancement de la
réaction : statistiquement, on peut trouver en tout point de la couche de mélange, toutes les frac-
tions de mélanges associée à un avancement de réaction allant de l’absence de réaction à la réaction
complète.
L’utilisation de fonctions de densités de probabilité (PDF) est sans doute la meilleure méthode
actuelle pour étudier une combustion aussi hétérogène. Une fonction statistique peut représenter
la diversité des mélanges et des taux d’avancement de la réaction, au moyen d’hypothèses faites
sur la loi de micro-mélange. Les différentes manières d’utiliser cette technique ont été testées par
Mobus et al. (106), précisément sur la flamme de Cheng, en utilisant un code RANS pour résoudre
la dynamique. Notons au passage que la diffusion différentielle, essentielle pour l’autoallumage,
complique sérieusement le calcul de la combustion par PDF. Les techniques de PDF sont également
utilisables en LES, au détriment de la prédictivité : dans une approche LES idéale, l’hétérogénéité
d’un écoulement devrait être résolue, pas modélisée.
Dans le souci d’aborder progressivement ce point difficile, les simulations présentées ici se
limitent à une combustion ’laminaire’, sans prendre en compte l’impact de la turbulence ou de la
ségrégation dans un modèle de combustion de sous-maille. Ce choix permet d’observer la résolution
de l’hétérogénéité de la combustion supersonique par la simulation des grandes échelles. Il faut
s’attendre à une surestimation du mélange, et donc de la vitesse de réaction.
6.4. Création d’un schéma cinétique simplifié 151
pauvre (φ = 0.123), donc le mélange entre les gaz brûlés et l’oxydant en excés donne en aval à une
température de combustion globale très inférieure.
(0.201 moles) dilué dans de l’azote (0.544 moles) et de l’eau (0.255 moles). Cependant, l’oxygène
de l’air ambiant peut participer à la combustion. Pour prendre en compte ce deuxième oxydant, il
faut raisonner sur une réaction entre hydrogène et oxygène non dilués : 0.8903 moles d’hydrogène
réagissant avec 0.4451 moles d’oxygène , à 1200◦ K. La description du produit de cette réaction
non diluée est donnée dans le Tableau 6.4. On y remarque que ce mélange P contient toujours les
espèces minoritaires OH, O et H. Par ailleurs, son enthalpie de formation est inférieure de 34% aux
−238922K j/mol de l’eau pure. on verra par la suite que ce mélange P pourra être utilisé pour la
cinétique d’auto-allumage.
temperature
Tst
t all temps
Le jeu de paramètres retenus favorise largement la première étape. Ainsi, le fait que les deux
étapes soient succéssives et que seule la deuxième soit exothermique implique un emballement
thermique brutal :
!1 !1
20000 ρYH2 ρYO2
Q1 = 1 10 exp(−
15
) (6.6)
RT WH2 WO2
!1
27000 ρYI
Q2 = 1 109 exp(− ) (6.7)
RT WI
Les exposants des espèces sont tous fixés à l’unité. Dans la formulation Arrhénius, ils sont
égaux aux coefficients stœchiométriques, mais le rapport ν0H2 = 2 × νO 0 fixe une forte dépendance
2
en hydrogène. Cela entraı̂nerait une réaction trop brutale. La dépendance linéaire sur chaque réactif
entraı̂ne une réaction associée plus progressive, c.a.d. plus de robustesse pour le calcul.
Les temps d’auto-allumages obervés en configuration HMI sont tracés sur la Fig. 6.8 en fonction
de la fraction de mélange. Les cinétiques utilisées sont les schémas de Yetter (155) de Baurle (10)
et le schéma PI, avec une thermodynamique complète (JANAF (139)). On observe que le temps
d’allumage le plus court est d’environ tHMI = 5 10−5 sec. Le schéma PI ne reproduit naturellement
pas la dépendant en fraction de mélange. Au contraire, il effectue un allumage en bloc de l’ensemble
des mélanges.
Si l’on regarde la montée en température de la couche de mélange HMI cf. Fig. 6.9, on obsere
une bonne correspondance entre la cinétique complète de Yetter et la cinétique simplifiée du schéma
PI. D’après ce test d’autoallumage HMI, le schéma PI réalise un auto-allumage similaire au schéma
complet de Yetter du point de vue de l’évolution en température. La dépendance envers le mélange
n’est pas capturée par le schéma PI. Pour retrouver une dépendance correcte, on prendra le schéma
de Baurle.
6.5 Simulations
Ces simulations ont fait l’objet d’une d’une session pleinière lors de la conférence internationale
CY-LES, complex effects in large eddy simulation tenue à Limassol en septembre 2005. L’article qui
l’accompagne est destiné à terme un journal, et est intégralement reproduit ci-après :
6.5. Simulations 157
10000
1000 Yetter
Baurle
PI
100
10
Autoignition time [s]
0.1
0.01
0.001
0.0001
1e-05
0 0.1 0.2
Mixture Fraction Z [-]
F. 6.8 – Comparaison de l’autoallumage HMI réalisé par le schéma complet de Yetter , le schéma
semi réduit de Baurle et le schéma PI
3000.0
Maximum instantaneous Temperature [K]
2500.0
2000.0
Yetter
PI
1500.0
1000.0
0.00000 0.00020 0.00040 0.00060
Time [s]
F. 6.9 – Comparaison des montées en température dans la couche HMI entre le schéma complet
de Yetter et le schéma PI
158 Chapitre 6. Combustion Supersonique
1 Introduction
2 Experiment Description
The SSB sketched in Fig. 1 produces an axi-symmetric flame from a sonic
pure hydrogen cold jet (Diameter D = 2.36mm) surrounded by a largely
supersonic (Mach 2) jet of hot products generated by a lean combustor. It
gives an auto-ignited flame, i.e. fuel and oxidizer are reacting as soon as they
are mixed. However, photographs of the experiment [1] showed that there is
no emission of light for 0 < x/D < 25. The SSB global equivalence ratio is
close to 0.1. The Reynolds number is 100 000 and the convective mach number
is 0.12.
Cheng and Wehrmeyer measured temperature, oxygen and nitrogen con-
centrations and velocity with coherent anti-Stokes Raman spectroscopy and
laser Doppler velocimetry. They obtained simultaneous measurements of tem-
perature and concentration of major species (H2 ,O2 ,N2 ,H2 O) and OH radi-
cals with ultraviolet spontaneous vibrational Raman scaterring combined with
laser induced predissociative fluorescence. Measurements were made in the
radial direction at the downstream locations x/D = 0.85, 10.8, 21.5, 32.3 and
43.1, x being the downstream distance and D = 2.36mm the inside diameter
of the fuel jet.
Dimensions
Nozzle exit inner diameter 17.78 mm
Fuel injector inner diameter 2.36 mm
Fuel injector outer diameter 3.81 mm
Vitiated Air Exit Conditions
Pressure 107 kP a
Temperature 1250 K
Mach number 2.0
Velocity 1420 m/s
O2 mole fraction 0.201
N2 mole fraction 0.544
H2 O mole fraction 0.255
Fuel Exit Conditions
Pressure 112 kP a
Temperature 540 K
Mach number 1.0
Velocity 1780 m/s
H2 mole fraction 1.0
Table 1. Experimental parameters.
160 Chapitre 6. Combustion Supersonique
x Hot products
with excess O 2 H2
y
3 Computational Domain
The computed domain is a half sphere of radius 70D shown in Fig. 2. The
plane boundary is set as a tri coaxial inlet. At the center, a sonic fuel inlet
is specified with a turbulence injection. It is surrounded by an oxidizer inlet
with a second turbulence injection. The turbulence injected in the fuel has an
integral scale of 0.8mm and a RMS velocity of 180m/s, while in the oxidizer
the integral scale is 6mm and the RMS velocity is 125m/s. A co-flow of 20m/s
is finally imposed on the outer ring. The convex boundary on the sphere is a
freestream outlet for subsonic zones. All boundary conditions are set with the
NSCBC method [7].
The mesh contains 807 520 tetrahedras for 139 560 nodes, allowing a res-
olution of 25 to 30 nodes in the diameter of both fuel an oxidizer jets. It is
designed to resolve the main structures occuring around the cylindrical mixing
layer between oxidizer and fuel in the first ten diameters.
4 Models
Simulations were performed with the code AVBP developed at CERFACS.
The multi-species fluid dynamics solver involves a Third order Taylor Galerkin
Compact (TTGC) scheme for space direction, and a third order Runge-Kutta
scheme for explicit time advancement. A specificity of the present flow is the
necessity of handling chemistry and turbulence but also shocks: LES is not
well established for non-reacting flows with shocks and is still exploratory
6.5. Simulations 161
when combustion must also be taken into account as it is the case here. For
the present computation, the following strategy was used:
1. a penalty method on Euler fluxes based on the Von-Neumann Richtmyer
viscosity [8] was employed to capture the shocks.
2. local sensors were used to add artificial viscosity in strong gradient zones
(typically near the inlet conditions).
3. chemistry was modeled explicitly by using Arrhenius rates. The focus is
on the resolved impact of turbulence on combustion.
4. Subgrid stresses are modeled with a Smagorinsky model.
In a very first approach, the chemistry model is a two-step scheme (Table 2).
The produced species pool P results from the equilibrium mixture reached
after a stoichiometric combustion between Oxydant and Fuel. The induction
step is mimicked through the non-exothermic formation of an intermediate
species pool I .
Chemical kinetics are fitted on autoignition. According to the work of
R. Knikker et al. [9], the ignition time of a mixing layer between hot oxi-
dizer and cold fuel can be estimated from the homogeneous mixing ignition
time. Simulation performed with the CHEMKIN package and the 9 species -
19 reactions scheme of Yetter et al. [10] predict the minimum ignition time
tHM
ign
I
= 5.6 10−5 s. The simplified chemistry model detailed in Table 2 is
designed to give the same minimal ignition time for a one-dimensional auto-
igniting mixing layer.
Step reaction A β E0
Induction 0.68H2 + 0.34O2 ⇒ I 1.0 1015 0 20000
9
Run off I ⇒ 0.73P 1.0 10 0 27000
Table 2. The two-steps auto-ignition scheme based upon the ignition time
162 Chapitre 6. Combustion Supersonique
5 Results
The simulation leads to a lifted flame as in experiement. Fig. 3.a shows
the abrupt rise of temperature of 1500K around x/D = 15 on the instanta-
neous image. This rise is smoothed downstream over 5D on the time averaged
image. The lift-off length is then around x/D = 17, to be compared with
the experimental value of x/D = 25. The fluctuating nature of the reactive
mixing layer is illustrated on Fig. 3.b where the production of P is plotted.
Strong large structures comparable to Kelvin-Helmoltz eddies can be seen in
the lift-off region 0 < x/D < 15. The pressure field in Fig. 3.c presents two
diamond-shock patterns relative to Oxidizer and Fuel injection. The shock
pattern relative to the fuel moves in response to the turbulence injected so
it vanishes quickly on the average field before x/D = 2.5. The shock pattern
relative to the Oxidizer gets stronger from the injected static pressure ratio
of 1.05 up to 1.3. The pressure oscillation, which reaches 0.6Bars, is strongly
coupled with the reaction zones (See Driscoll et al. [11] ).
Along the longitudinal axis, the mean velocity, temperature, water concen-
tration and mixture fractions are compared to experimental data. In Fig. 4.a,
the axial velocity decreases from hydrogen to hot products bulk velocity (1780
to 1420 m/s) in the first ten diameters. An acceleration about 200m/s occur
because of the diamond-shock pattern. The mean inlet velocity is slightly
higher than the experimental one because of mesh requirements, to keep the
correct mass flow rate
The mean temperature and water concentration show similar behaviors
(Figs. 4.b and 4.c), reaching their maximum around x/D = 17, then decreas-
ing slowly to 50% of the maximum value at x/D ∼ 60. An induction zone
is visible up to x/D ∼ 10 then the equilibrium is reached as fast as the ex-
periment. This too fast combustion is due to the simplified kinetic model of
Table 2. Finally, the slow decrease is the result of diffusion that is too fast in
the range 20 < x/D < 70, where the Smagorinsky model is active on an un-
derresolved grid. This effect is also visible on the mixture fraction longitudinal
(Fig. 4.d) and transversal profiles (Fig. 5.d).
6 Conclusion
Acknowledgements
References
1. T.S. Cheng, J.A. Wehrmeyer, R.W. Pitz, O. Jarrett, G.B. Northam: Combus-
tion and Flame 99, 157–173 (1994)
2. P. Toone: (2002), “A computation fluid dynamic model of a supersonic axi-
symmetric jet using a beta probability function combustion model”, Tech. rep.,
Departement of Mechanical Engineering, Purdue University
3. H. Mobus, P. Gerlinger, D. Bruggemann: Comb. Flame 132(1), 3–24 (2003)
4. P. Colucci, F. Jaberi, P. Givi, S.B. Pope: Phys. Fluids (1998)
5. F. Jaberi, P. Colucci, S. James, S.B. Pope: J. Fluid Mech. 401, 85–121 (1999)
164 Chapitre 6. Combustion Supersonique
2000
2000
1000
1000
0 0
0 20 40 0 20 40 60
x/D x/D
a) b)
H2O mole fraction Mixture fraction
0.5
0.8
0.25
0.4
0 0
0 20 40 60 0 20 40 60
x/D x/D
c) d)
Fig. 4. Mean profiles on the longitudinal axis. Symbols: expt, line: LES.
2000
Velocity [m/s] Temperature [K]
3000
1500
1500 500
x/D=64.7
1000 0 0
2000 1500
3000
x/D=43.1
500 1500 x/D=43.1
0 1000 0
2000 1500
3000 x/D=32.3
1500 500
x/D=32.3
1000 0 0
2000 1500
3000
500 x/D=21.5
1500
x/D=21.5
0 1000 0
2000 1500
3000
500
x/D=0.85
x/D=0.85
0 0
-10 -5 0 5 10 -10 0 10
x/D x/D
a) b)
H2O Mole fraction 0.025
Mixture fraction
0.4
0.05
0.2
0.4 x/D=64.7 x/D=64.7
0 0
0.1
0.2
0.4 x/D=43.1
x/D=43.1
0 0
0.2 0.2
0.4
x/D=32.3
x/D=32.3
0 0
0.2 0.4
0.4
x/D=21.5 x/D=21.5
0 0
0.2
1.5
0.4
x/D=10.8
x/D=10.8
0 0
0.2
x/D=0.85 x/D=0.85
0 0
-10 -5 0 5 10 -10 -5 0 5 10
x/D x/D
c) d)
Fig. 5. Mean profiles for six transversal planes. Symbols: expt, line: LES.
6.5. Simulations 167
0 0
200 300
400 x/D=43.1 x/D=43.1
600
0 0
200 300
400 600
x/D=32.3
x/D=32.3
0 0
200 300
400 600
x/D=21.5 x/D=21.5
0 0
200 300
400 600
x/D=10.8
x/D=10.8
0 0
200 300
x/D=0.85 x/D=0.85
0 0
-10 -5 0 5 10 -10 -5 0 5 10
x/D x/D
a) b)
H2O RMS Mole fraction
0.1
x/D=64.7
0.1
0
x/D=43.1
0.1
0.1 x/D=32.3
0.1
x/D=21.5
0
0.1
x/D=10.8
0
x/D=0.85
0
-10 -5 0 5 10
x/D
c)
Fig. 6. Root Mean Square profiles for six transversal planes. Symbols: expt, line:
LES.
168 Chapitre 6. Combustion Supersonique
1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction Mixing w/o reaction
0.75 0.75
H2 Mole fraction
H2 Mole fraction
0.5 0.5
0.25 0.25
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
a) Hydrogène b) Oxygène
0.1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction
0.08
0.75
H2O Mole fraction
OH Mole fraction
0.06
0.5
0.04
0.25
0.02
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
2500
Temperature [K]
2000
1500
1000
500
0
0 0.05 0.1
Mixture fraction
e) Température
F. 6.10 – Répartition du mélange à x/D = 0.85 et r/D = 0.65.◦ expérience, • simulation.
170 Chapitre 6. Combustion Supersonique
Le second point de mesure, x/D = 10.8 et r/D = 0.65, se situe au milieu de la zone d’induction,
dans la couche de mélange instationnaire des deux jets coaxiaux. Les observations expérimentales
(◦) de la Fig. 6.11 , sont toujours sur la ligne de mélange sans réaction (trait interrompu). La dis-
persion des points de mesure a augmenté du fait du mélange. La dispersion de la température reste
entre 1000 et 1500 ◦ K.
Les points de la simulation numérique (•) montrent une réaction d’induction en plein essor. En
effet la masse molaire des radicaux OH atteint son maximum, 0.06 (Fig. 6.10.c). Les réactifs sont
en partie consommés, au vu des quelques points qui s’éloignent de la courbe mélange sans réaction
(trait interrompu) , comme le prouve la courbe de l’hydrogène (Fig. 6.10.a,.b). Enfin, quelques
points concernant la fraction molaire de l’eau ou la température se sont dispersés vers la courbe
d’équilibre adiabatique (trait plein, Fig. 6.10.d,.e), prouvant que quelques poches de fluides sont
passées par cet endroit en ayant complètement réagi.
La surestimation du mélange est encore en cause. Cependant, le calcul montre une grande
dispersion des points de mesure, illustrant la capacité de la SGE à reproduire une réaction très
hétérogène : des particules fluides à différents stades de la réaction (non-réctives, en phase d’induc-
tion, en phase d’emballement thermique) passent par le même endroit.
6.6. Analyse de l’hétérogénéité 171
1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction Mixing w/o reaction
0.75 0.75
O2 Mole fraction
H2 Mole fraction
0.5 0.5
0.25 0.25
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
a) Hydrogène b) Oxygène
0.1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction
0.08
0.75
H2O Mole fraction
OH Mole fraction
0.06
0.5
0.04
0.25
0.02
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
2500
Temperature [K]
2000
1500
1000
500
0
0 0.05 0.1
Mixture fraction
e) Température
F. 6.11 – Répartition du mélange à x/D = 10.8 et r/D = 0.65.◦ expérience, • simulation.
172 Chapitre 6. Combustion Supersonique
Le troisème point de mesure, x/D = 21.5 et r/D = 1.1, est placé en amont de la flamme visible
sur la photographie de l’expérience, Fig. 6.2. Les observations expérimentales (◦) de la Fig. 6.12,
sont dispersées de facon presque homogène entre les courbes d’équilibre adiabatique (trait plein)
et de mélange sans réaction (trait interrompu). Le nuage de points de la température est néanmoins
plus proche de l’état non-réagi (Fig. 6.12.c,.e), indiquant encore une prépondérence de la phase
d’induction. On note également l’accumulation des points vers les zones les plus pauvres (Z = 0)
causée par la combustion globalement pauvre.
Les points de la simulation numérique (•) sont tous regroupés prés l’équilibre adiabatique, ce qui
signifie que la réaction a atteint son terme. Le choix d’une thermodynamique simplifiée réglée sur la
stœchiométrie explique les décalages visibles entre l’équilibre atteint par le calcul et l’équilibre ob-
tenu avec une chimie détaillée, particulièrement visibles sur l’eau et l’ion hydroxyde (Fig. 6.12.c,.d)
dans les zones riches Z > 0.05. Enfin, le domaine de dispersion du mélange 0 < Z < 0.1 correspond
à la diversité observée dans l’expérience.
Ainsi, la comparaison calcul/expérience montre que la thermodynamique simplifiée s’approche
de facon satisfaisante de l’équilibre adiabatique. De plus, la SGE est en mesure de reproduire
l’hétérogénéité du mélange. Il faudra à l’avenir une analyse plus quantitative de la dispersion des
points pour préciser cette tendance. En revanche, la combustion hétérogène simulée atteint son terme
trop tôt, toujours en raison de la surestimation du mélange.
6.6. Analyse de l’hétérogénéité 173
1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction Mixing w/o reaction
0.75 0.75
O2 Mole fraction
H2 Mole fraction
0.5 0.5
0.25 0.25
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
a) Hydrogène b) Oxygène
0.1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction
0.08
0.75
H2O Mole fraction
OH Mole fraction
0.06
0.5
0.04
0.25
0.02
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
2500
Temperature [K]
2000
1500
1000
500
0
0 0.05 0.1
Mixture fraction
e) Température
F. 6.12 – Répartition du mélange à x/D = 21.5 et r/D = 1.1.◦ expérience, • simulation.
174 Chapitre 6. Combustion Supersonique
Le quatrième point de mesure, x/D = 32.3 et r/D = 1.1, se situe à l’intérieur de la flamme.
Les observations expérimentales (◦) de la Fig. 6.13 sont toujours très dispersées. La concentra-
tion d’ion hydroxyde atteint son maximum (Fig. 6.13.c), ce qui implique que les réactions in-
termédiaires d’induction sont toujours présentes. La concentration d’eau dépasse même l’équilibre
adiabatique (Fig. 6.13.d) du coté riche, tandis que les réactifs sont en dessous de leur coubre
d’équilibre (Fig. 6.13.a,.b). Enfin, la concentration en eau varie entre les courbes de mélange non-
réagi, et de mélange à l’équilibre (Fig. 6.13.d).
Les points de la simulation numérique (•) s’accumulent vers les zones pauvres. L’eau ne dépasse
pas la concentration d’équilibre (Fig. 6.13.d) conformément à la cinétique et la thermodynamique
retenue : l’espèce intermédiaire I est pauvre en eau, donc la concentration locale ne peut dépasser
celle de l’espèce des produits P.
Cette position est le lieu le plus hétérogène dans l’expérience, avec des concentrations en ions
hydroxyde et en eau qui dépassent largement l’équilibre adiabatique. Les choix thermodynamiques
et cinétiques présents pour la combustion ne reproduisent pas ces phénomènes. Le dépassement des
concentrations d’équilibre adiabatique dans l’expérience peut bien sûr être reproduit en utilisant une
thermodynamique et une cinétique chimique complète.
6.6. Analyse de l’hétérogénéité 175
1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction Mixing w/o reaction
0.75 0.75
O2 Mole fraction
H2 Mole fraction
0.5 0.5
0.25 0.25
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
a) Hydrogène b) Oxygène
0.1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction
0.08
0.75
H2O Mole fraction
OH Mole fraction
0.06
0.5
0.04
0.25
0.02
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
2500
Temperature [K]
2000
1500
1000
500
0
0 0.05 0.1
Mixture fraction
e) Température
F. 6.13 – Répartition du mélange à x/D = 32.3 et r/D = 1.1. ◦ expérience, • simulation.
176 Chapitre 6. Combustion Supersonique
Le cinquième et dernier point de mesure, x/D = 43.1 et r/D = 0, est plus en aval à l’intérieur
de la flamme. Les observations expérimentales (◦) de la Fig. 6.14 lorsqu’elles sont comparés aux
positions plus en amont (Fig. 6.13, Fig. 6.12 ) commencent à se concentrer vers l’équilibre. La
combustion reste très hétérogène, avec des concentrations en ions OH toujours élevées, entre 0.02
et 0.06. Par ailleurs, on trouve encore beaucoup de mélanges riches 0.025 < Z < 0.1 à plus de
quarante diamètres d’injection en aval, alors que la flamme globale est pauvre. Cela signifie que le
jet central ne s’est pas encore complétement mélangé à cette position
Au contraire, les points de la simulation numérique (•) sont entièrement limités aux zones
pauvres (Fig. 6.14.e). Cette observation signifie que la position x/D = 43.1 est hors de la zone
utile du calcul. La région centrale est insuffisamment résolue pour conserver des poches riches
convectées dans l’écoulement.
La dispersion du mélange est un bon indicateur de l’hétérogénéité locale. La Fig. 6.14 montre
que le domaine utile de la simulation est bien inférieur à quarante diamètre en aval, puisqu’aucune
poche de mélange riche n’apparaı̂t aussi loin, contrairement à l’expérience.
L’analyse présente montre que la SGE appliquée à la flamme de Cheng reproduit une longue
phase d’induction, suivie d’un emballement thermique brutal. Par ailleurs, les zones les plus résolues
autorisent la combustion hétérogène qui a lieu dans ce type de flamme. Une augmentation significa-
tive de la finesse du maillage permettrait de diminuer la combustion initiale, trop rapide par rapport
à l’expérience, et de conserver l’hétérogénéité de l’écoulement réel plus loin en aval. La thermody-
namique et la cinétique chimique simplifiée reproduisent le comportement global de la flamme au
premier ordre. Des progrés sont encore possibles sur les quatre paramètres du schéma cinétique
PI, actuellement trop brutal lors de la transition induction/emballement thermique. Concernant
les phénomènes de dépassement d’équilibre, il est préférable d’attendre la puissance de calcul
nécessaire pour un schéma complet (schéma de Yetter (155)), car les schémas simplifiés soulèvent
plus de questions qu’ils n’en règlent.
6.6. Analyse de l’hétérogénéité 177
1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction Mixing w/o reaction
0.75 0.75
O2 Mole fraction
H2 Mole fraction
0.5 0.5
0.25 0.25
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
a) Hydrogène b) Oxygène
0.1 1
Adiab. equil. Adiab. equil.
Mixing w/o reaction
0.08
0.75
H2O Mole fraction
OH Mole fraction
0.06
0.5
0.04
0.25
0.02
0 0
0 0.05 0.1 0 0.05 0.1
Mixture fraction Mixture fraction
2500
Temperature [K]
2000
1500
1000
500
0
0 0.05 0.1
Mixture fraction
e) Température
F. 6.14 – Répartition du mélange à x/D = 43.1 et r/D = 0.◦ expérience, • simulation.
178 Chapitre 6. Combustion Supersonique
Chapitre 7
Conclusion
7.1 Récapitulation
Pour commencer la conclusion, voici la récapitulation du contenu scientifique de cette thèse.
L’auto-allumage, les jets supersoniques sous-détendus et la combustion supersonique en sont les
trois thèmes principaux. L’état de l’art (Chap.2) a défini un angle de recherche pour chacun de ces
thèmes. Après introduction des équations résolues et de la méthode numérique employée (Chap.3),
les travaux présentés ci-avant ont mené aux résultats suivants :
Le phénomène d’auto-allumage a montré au chapitre 4 un comportement complexe en présence
de diffusion différentielle. Le couple méthane/oxygène s’auto-allume plus lentement lorsque les
réactifs sont initialement séparés (LMI), comparé au cas idéal d’une configuration prémélangée
(HMI), avec un retard de 100% sur le temps d’allumage. A l’opposé, le temps d’auto-allumage
du couple hydrogène/oxygène ne varie presque pas entre les deux configurations. Il s’agit d’un
cas particulier où le mecanisme de diffusion différentielle compense le délai de l’allumage dû à la
ségrégation des réactifs.
Le chapitre 5 montre comment un code de simulation des grandes échelles utilisant un schéma
centré de troisième ordre peut être adapté aux écoulements supersoniques. La méthode utilise prin-
cipalement la viscosité artificielle de Von Neumann- Richtmyer (148) pour calculer des chocs sans
modifier le schéma numérique et le maillage. Outre les cas test du tube à choc et de l’onde de
Riemann, cet approche a permis de traiter avec succés les chocs des jets sous-détendus. L’analyse
poussée de ces simulations valide cette approche, et montre la richesse des informations créées par
une SGE.
Le chapitre 6 décrit une simulation de flamme supersonique auto-allumée. Cette flamme est
un cas test de la NASA destinée à l’étude de la cinétique chimique intervenant dans les super-
statoréacteurs. La simulation utilise une cinétique chimique simplifiée grâce aux résultats du cha-
pitre 4, et reproduit clairement la longue zone d’induction endothermique caractéristique de cette
flamme.
Cette somme d’information ouvre de nouvelles possibilités. D’un point de vue industriel, les
configurations de géométrie complexe présentant des phénomènes transitoires sont désormais à
portée. D’un point de vue scientifique, l’outil développé peut directement prendre en compte une
thermodynamique de gaz réels, en vue de configurations supercritiques. En résumé, cette thèse a
179
180 Chapitre 7. Conclusion
produit des avancées notables sur le chemin de la simulation des moteurs-fusées cryotechniques par
la SGE, et les communications scientifiques qui assurent la légitimité de l’ensemble.
7.2 Analyse
Il est utile d’ajouter une analyse des options prises durant cette thèse. L’ensemble de ces tra-
vaux est destiné à s’approcher de la simulation de phénomènes supersoniques réactifs dans des
géométries complexes. Par la suite, ces simulations seront étendues à des phénomènes transitoires
et supercritiques. Ainsi, les techniques les plus évoluées observés dans l’état de l’art ne peuvent être
retenues :
En premier lieu, l’étude du chapitre 4 prouve que l’hypothèse très commune d’un nombre de
Lewis pris égal à un est définitivement inutilisable dans une combustion par auto-allumage. Les
résultats de l’auto-allumage du couple hydrogène/oxygène, comparés au couple méthane/oxygène,
ont fait l’objet d’une session pleinière au Quatrième colloque international sur la technologie des
lanceurs : propulsion liquide pour lanceurs spatiaux à Liège en décembre 2002, puis d’une pu-
blication dans Combustion Science and Technology en 2003 (82) en collaboration avec le Dr R.
Knikker.
En second lieu, l’état de l’art montre que les écoulements supersoniques ont toujours été abordé
dans la simulation des grandes échelles par des techniques de solveurs de Riemann. Le besoin d’un
maillage non structuré et d’un fluide réactif suivant une loi de gaz réels rend ce choix technique in-
approprié. Cette constatation a dirigé les travaux du chapitre 5. Un article destiné au journal Physics
of fluids présente la simulation de plusieurs jets sous détendus, puis l’exploitation intensive de leurs
résultats.
Enfin, l’étude de combustion supersonique du chapitre 6 montre la création d’une cinétique
chimique simplifiée destinée au processus d’auto-allumage, un point qui n’existait pas dans la
littérature. Les résultats d’une simulation utilisant cette cinétique ont fait l’objet d’une session
pleinière lors de la conférence internationale CY-LES, complex effects in large eddy simulation te-
nue à Limassol en septembre 2005. Un article, destiné à terme au journal Combustion and Flame,
reprend l’ensemble de ces travaux et termine ce chapitre.
Ainsi, par trois fois ces travaux ont privilégié des techniques simples comparées aux standards
actuels. Il est amusant de constater que le retour à la ’sobriété’ face à un problème plus complexe
suit l’un des obscurs préceptes du célèbre maı̂tre de stratègie chinois Sun Bin :”C’est par son
aspect fort qu’une chose perceptible l’emporte sur une autre. Néanmoins, il est impossible d’utiliser
l’aspect fort d’une chose pour venir à bout de tout. C’est pourquoi les raisons pour lesquelles une
chose l’emporte sur une autre sont toujours les mêmes, mais les choses qui peuvent l’emporter sur
d’autres diffèrent les unes des autres.” L’art de la Guerre, e siècle avant J.-C. Cette citation légère
mise à part, nous retiendrons que la stratégie utilisée nous permet de continuer le développement de
l’approche SGE vers des configurations plus réalistes des points de vue scientifique et industriel.
7.3. Perspectives 181
7.3 Perspectives
Pour terminer, voici les travaux actuellement en cours qui découlent directement de cette thèse.
Deux nouvelles thèses se déroulent au CERFACS depuis le mois de Novembre 2005, toujours dans
le cadre de l’allumage des moteurs fusées cryotechniques, soutenues par la SNECMA et le CNES.
La première est menée par Guilhem Lacaze. Son objectif final est de réaliser des simulations
d’écoulements réactifs transitoires puis de jets fortement sous-détendus comprenant de la combus-
tion supersonique et une injection diphasique. Le premier point est largement entamé aujourd’hui,
avec la simulation de l’allumage ”dur” d’une chambre de combustion expérimentale. Dans ce dis-
positif installé au DLR de Lampoldhausen, une chambre initialement pleine d’azote est progressive-
ment remplie d’hydrogène et d’oxygène. Un allumage par laser déclenche la combustion brutale du
prémélange, faisant augmenter la pression jusqu’à 8 bars. La présence d’un jet supersonique coaxial
en injection, d’une tuyère amorcée en sortie de chambre, et d’une phase de remplissage non réactive
extrêmement longues font de ce cas test un calcul particulièrement difficile. Des résultats concrets
sont déjà disponibles, et ont fait l’objet d’une session pleinière à l’atelier organisé par la SNECMA
”Third International Workshop on Rocket Combustion Modeling ” tenu à Versailles en Mars 2006.
La seconde thèse est conduite par Thomas Schmitt. Elle se focalise sur les phénomènes super-
critiques. En effet, lors de l’allumage les températures et pressions atteintes font que les gaz ne
suivent plus une loi de gaz parfaits. Il faut donc étudier les conséquences d’une telle perturbation du
système déquations habituellement étudié, puis adapter en conséquence le code de SGE voulu, ici
AVBP. Le travail initial consiste à choisir une nouvelle équation d’état, puis de générer la thermo-
dynamique qui lui est associée. Par la suite, les lois de transports d’espèces et de chaleur devraient
évoluer. Pour finir, les simulations d’écoulement supercritiques devront s’accommoder de gradients
de densités au comportement numérique aussi critique que des chocs hydrodynamiques.
Au terme de la thèse présente, ”Allumage des moteurs fusées cryotechniques”, l’effort de re-
cherche dans ce domaine continue en prise directe avec les résultats présentés dans ce document.
D’ici peu, de nouvelles simulations vont voir le jour, prenant en compte de nouveaux phénomènes
physiques et profitant de l’augmentation continue de la puissance de calcul. Peut être pourront-elles
conforter un peu plus l’idée que la simulation numérique est un outil d’investigation scientifique à
part entière, mais comme dirait R. Kipling, ”ceci est une autre histoire...”
182 Chapitre 7. Conclusion
Chapitre 8
Remerciements
Cette thèse est la somme de multiples contributions. En premier lieu, plusieurs points précis ont
bénéficé du travail de personnes extérieures à l’équipe initiale :
– La première partie de l’introduction, c.a.d. le résumé épistémologique de l’histoire des mo-
teurs fusées, doit beaucoup aux entretiens électroniques et téléphoniques avec Marie Cathe-
rine Rey (Musée Guimet, Conservateur section Chine bouddhique) et Gérard Emptoz (Uni-
versité de Nantes, Professeur émérite d’histoire et techniques).
– La réalité des problèmes techniques posés par les moteurs fusées m’a été révélée par l’équipe
FLTC de la Snecma Division Grosse Propulsion Liquide. En plus de l’apport fondamental de
Stéphan Zurbach, je remercie Didier Saucereau pour ses lumières sur les simulations RANS.
– Les conditions exceptionnelles des écoulements supersoniques ont soulevé de multiples ques-
tions ouvertes. Je remercie toux ceux qui ont pris le temps d’y répondre par voie électronique,
notamment le Dr. Jean Louis Estivalèzes (Toulouse, ONERA), le Pr. G. Mungal (Stanford,
Faculty of thermosciences group), le Pr. Olivier Chazot (Bruxelles, Von Karman Institute),
le Dr. Hiroshi Katanoda (Fukuoka, University of Kitakyushu) et le Pr. Valery Zapryagaev
(Novossibirsk, ITAM).
– L’emploi de la méthode de Von Neumann-Ricthmyer au lieu de méthodes plus classiques a
soulevé une discussion intéressante. Je remercie le Pr. Dale Pullin (California Institute of tech-
nology), le Dr. Benoı̂t Fiorina (Stanford, Center for Turbulence Research) et le Pr. Andrew
Cook (Lawrence Livermore National Laboratory) pour leurs réponses rapides et décisives sur
le sujet.
– L’analyse linéaire de stabilité est un domaine scientifique à part entière. La collaboration de
plusieurs mois avec le Pr. Jean Christophe Robinet (Paris, Ecole Nationale Supérieure des
Arts et Métiers, Lab. SINUMEF) a dégagé des résultats tout à fait passionnants. Le stage de
Master de Zaki Mazen fut entièmenent dédiée à l’étude du jet de Seiner d’après les résultats
de simulation présents. Pour cet apport inespéré, je souhaite signifier ma reconnaissance au
Pr. Robinet.
– Des travaux importants ont été effectués sur l’expérience d’allumage ”dur” installée par le
DLR à Lampoldshausen, en collaboration avec Guilhem Lacaze. Les résultats seront détaillés
dans le manuscrit de sa thèse. En attendant, je remercie le Dr. Michael Oschwald (Lampold-
shausen, DLR) pour ses réponses précises, et Guilhem Lacaze, pour sa ténacité face à un
183
184 Chapitre 8. Remerciements
Je remercie également l’aide indirecte mais essentielle de Gabriel Staffelbach, Matthieu Boi-
leau, Stéphane Pascaud, David Seigneuric, Vincent Hennebert, Amélie Baudon, Pierre Crovisier et,
tout particulièrement, Guenaëlle Mesthé durant ces trois années.
Pour finir, je remercie le Dr. Stéphan Zurbach, à l’origine de cette thèse, et le Dr. Bénédicte
Cuenot, qui a su m’épauler au long de cette recherche.
Chapitre 9
m × 4.002 1014
Fterre = (9.1)
d2
185
186 Chapitre 9. Annexe 1
F. 9.1 – Isaac Newton est à l’origine, entre autres, des lois régissant le vol spatial
Le travail T lib nécessaire pour emmener ce corps de sa distance initiale d0 à l’infini, c’est à dire de
le libérer du champ gravitationnel, s’écrit :
Z ∞ " #∞
Mm 1 K Mm
T lib = k 2 = kMm − = (9.2)
d0 d d d0 d0
On peut donc en déduire la vitesse à donner à ce corps au départ de la surface de la terre (i.e.
d0 = 6.37 106 m), pour que son énergie cinétique corresponde au travail de libération :
s
mVlib2
K Mm 2kM
= => Vlib = = 11 200m/s (9.3)
2 d0 d0
En poussant plus loin les développements des lois du mouvement de Newton, on peut retrouver
les trois lois de Kepler (1605) :
1. Le mobile décrit une orbite elliptique dont le point d’attraction est un foyer ;
2. La droite mobile-point d’attraction balaye des surfaces proportionnelles au temps qu’elle met
à les balayer ;
3. Le carré du temps de révolution est proportionnel au cube du demi-grand axe a de l’ellipse.
Le paramêtre principal d’un mobile dans un champ gravitationnel est la vitesse radiale Vr (r).
Comparée à la vitesse de libération locale Vlib (r), elle détermine si le mobile est captif, c.a.d. sa-
tellisé, ou s’il est libre. Par exemple, les comêtes peuvent être sur une orbite terrestre, solaire ou
extra-solaire. Par conséquent, l’incrément de vitesse qu’un lanceur peut donner à son projectile
détermine l’essentiel de la performance du lanceur.
9.2. Loi du vol d’une fusée 187
Vz (t) = z(t)
˙ = −gt + Vz (t0 ) ⇔ ∆Vz = g(t1 − t0 ) (9.5)
portionnels au carré de la vitesse et qui ŕeduisent la portée, ou encore la force de Coriolis qui crée
une déviation latérale,
Fequiv
[ln (m(t))]tt10 = V(t1 ) − V(t0 ) (9.9)
qequiv
1
Cette caractéristique rendue publique notamment par R. Goddard, fit un petit scandale au début du vingtième siècle,
quand le bon sens considérait que la présence d’un milieu matériel sur lequel s’appuyer était forcément préférable. Ce
contre-sens vient de la propulsion par hélice qui évoque plus facilement un tire-bouchon dans du liège qu’une éjection de
matière. Une hélice reste cependant un moteur à réaction, pouvant fonctionner efficacement dans un fluide non visqueux.
188 Chapitre 9. Annexe 1
!
Fequiv m(t1 )
ln = ∆V (9.10)
qequiv m(t0 )
L’équation finale du moteur idéal, Eq. 9.10, montre que les deux paramêtres essentiels d’un
moteur fusée sont la force équivalente d’une part et le rapport entre masses initiales et finales d’autre
part. Cette approche fut introduite en 1900 par le russe Tsiolkovsky. Il introduisit les quantité I sp
dite ’impulsion spécifique’, et le rapport de masse Rm , ce qui mène à l’Eq. 9.11.
Z t1 h F i
e ject (t) + F press (t)
I sp × g × ln(Rm ) = ∆V = ∂t (9.11)
t0 m(t)
Ces deux paramêtres sont les deux clefs des nouveaux moteurs fusées. Le rapport de masse devient
le paramêtre principal de construction. Il donne en premier lieu l’intérêt des fusées multi-étages.
L’application numérique montre par exemple que la fusée lunaire Saturn V dans une version mono-
étage verrait sa masse multipliée par 2.38, comme indiqué dans le Tableau 9.2. Un tel engin est
beaucoup plus difficile à mettre au point du point de vue de la résistance des matériaux. L’impulsion
Fequiv
spécifique est une quantité plus abstraite : qequiv q = I sp × g. Multipliée par l’accélération de la
pesanteur, elle donne une vitesse légèrment supérieure à la vitesse déjection des gaz propulsifs,
mais l’intérêt est ailleurs. L’impulsion spécifique est donnée au premier ordre par la nature des
carburants, et seulement au second ordre par l’optimisation technique du moteur. Ainsi, connaissant
la charge utile et l’incrément de vitesse à donner, on déduit des travaux de Tsiolkovsky la nature
des carburants et la masse du lanceur. Le nombre d’étages dépend alors d’une optimisation entre
efficacité de la propulsion et contraintes structurelles.
Pour vérifier la pertinence de ces outils, voici l’exemple exemple précis du V2. Le cahier des
charges se résume ainsi : ”une arme de type missile balistique capable d’envoyer une charge militaire
d’une tonne à 300 km de distance”. Le tir optimal étant réalisé à 45 degrés, on trouve une vitesse
initiale de 1715 m/s. Un calcul ballistique plus complet prenant en compte la résistance de l’air
augmente cette vitesse initiale de 25%, soit environ 2300m/s. Le couple Oxygène/Alcool est le
mélange le plus performant utilisable, et permet d’espérer une impulsion spécifique proche de 200s.
Il vient alors le rapport de masse Rm = exp(Vinit /(I sp.g)) = 3.23. Il reste donc à construire une
structure de fusée capable de résister à la masse totale Mtot = 3.23 ∗ (Mcharge + M struct + Mmot )
et le moteur correspondant, capable de soulever cette masse totale. La solution de ce problème en
1944 est une structure de 1 tonne équipée d’un moteur de 2 tonnes, soit une masse au décollage de
presque 13 tonnes, comme indiqué dans le Tableau 9.1.
9.2. Loi du vol d’une fusée 189
191
192 Chapitre 10. Annexe 2
F. 10.1 – Réseau de lignes caractéristiques pour un jet sous-détendu obtenu par Thompson (143)
avec 5 lignes caractéristiques (8 pour l’expansion initiale), Mach = 1, NPR = 2
F. 10.2 – Réseau de lignes caractéristiques pour un jet sous-détendu obtenu avec 20 lignes ca-
ractéristiques , Mach = 1, NPR = 2
10.1. Méthode des caractéristiques 193
Cette routine est réutilisée pour étudier le comportement du réseau de caractéristiques pour
différents rapports de pression. Afin de prolonger le calcul, si deux caractéristiques se croisent, celle
du coté aval est stoppée. Par ailleurs, le calcul fait apparaı̂tre un choc droit (ligne verticale) dés
que les rapports de pression et de densité entre l’axe et l’extérieur permettent de vérifier la relation
de Rankine-Hugoniot. Les equations utilisées correspondent à un écoulement bidimensionnel plan.
Un autre jeu d’équations accompagné d’un algorithme plus complexe permet de résoudre les cas
axisymétriques.
Dans le cas d’un jet faiblement sous-détendu (Fig. 10.3,a), la structure en diamant semble pou-
voir se répéter à l’infini. Dans un ecoulement visqueux, la turbulence et la viscosité dissipent cette
structure. Pour un rapport de pression un peu plus élevé (Fig. 10.3,b) des cellules ovoı̈des ont rem-
placé les cellules en diamant. La structure est encore répétitive. Des lignes caractéristiques dispa-
raissent à chaque fin de cellule, diminuant la précision du calcul, c’est pourquoi on ne tient pas
compte de la dilatation des cellules ver l’aval. Un jet fortement sous-détendu possède un réseau de
caractéristiques plus complexe (Fig. 10.3,c). En effet, le faisceau de détente commence à coalescer
depuis l’intérieur du jet dès la reflexion sur la ligne de glissement. Cette coalescence se traduit par un
choc faible à l’intérieur de la cellule de détente. Ce choc faible s’intensifie vers l’aval, avant d’être
terminé par un choc droit, ce qui délimite un volume ovoı̈de fermé que l’on appelle parfois ”bouteille
de Mach”. La même structure est visible pour les jets très fortement sous-détendus (Fig. 10.3,d). Le
choc faible remonte jusqu’à l’injection lorsque l’on augmente le rapport de pression.
La méthode des caractéristiques reproduit explicitement dans un temps très faible les principaux
élements non-visqueux d’un jet sous-détendu. Elle donne par ailleurs des informations utiles sur
l’origine de chaque élément non-visqueux présent dans un jet sous-détendu. Cet outil est aujourd’hui
utilisé pour calculer très précisément l’intérieur de la cellule de détente (112). Ce type d’étude
permet de régler l’alimentation en gaz raréfié d’un accélérateur de particules.
194 Chapitre 10. Annexe 2
a) Mach = 1, NPR = 1, 05
b) Mach = 1, NPR = 2
Pi/Po=10.00
c) Mach = 1, NPR = 10
200 x H
l’Eq. 10.1 :
P stat + ρU2
2
M2γ
!
NPRtot = = 1+ NPR (10.1)
P stat ∞ 2
La longueur de ces bouteilles c.a.d. la position normalisée x M /D dépend de ces trois paramêtres
le la facon suivante : xDM ∝ (Me )∼1 (γ)∼0.5 (NPR)∼0.5 . Trois grandes classes de méthodes permettent
d’obtenir ce résultat. Il y a tout d’abord, de multiples corrélations expérimentales, comme celle de
Lewis & Carlson (97), donnée en Eq. 10.2, valide sur la plage NPR ∈ [1, 100], d’une précision de
10%, avec un écart maximum de 31% par rapport aux expériences de l’époque.
x
= 0.69Me (γNPR)0.5 (10.2)
D
Ensuite on trouve les méthodes de calculs aux caractéristiques comme celle de Hanson & Pal-
mer (112) qui profitent du caractère hyperbolique des équations de Naviers-Stokes en écoulement
supersonique. Sur la plage NPR ∈ [1, 1000], un code de calcul prédit la position du disque (Corre-
lation Eq. 10.3) avec une erreur inférieure à la dispersion des résultats expérimentaux.
x
= 0.8736 (NPR)0.52069 (10.3)
D
Enfin, des simulations plus classiques de volumes finis peuvent être utilisées, comme l’ont fait
Cumber (38). Un code volume finis RANS, avec un schéma du deuxième ordre, un modèle de
turbulence type k − et une correction du taux de dissipation proposée par Sarkar, fonctionne sur
la plage PP∞t ∈ [1, 40]. Malgré la dispersion des résultats expérimentaux proposés, le code semble
sur-estimer jusqu’à 30% les dimensions des bouteilles de Mach. Son comportement est représenté
par la corrélation Eq. 10.4.
x
= 1.2182 (NPR)0.4236 (10.4)
D
L’ensemble des sources menant à ces résultats est listé de maniére non exhaustive dans le ta-
bleau 10.1, puis représenté sur le graphique 10.4. Le graphique supérieur de la Fig. 10.4 montre
196 Chapitre 10. Annexe 2
T. 10.1 – Compilation des différentes sources concernant le rapport x M /D. M.O.C : Méthodes aux
caractéristiques, F.V. Méthodes Volumes Finis
des résultats peu dispersés, donnant une erreur de 20% maximum sur l’ensemble des rapports de
pression. Cette dispersion est d’abord présente dans les résultats expérimentaux, notamment avec
les observations de Ewan & Moodie (57), inférieures de 10% aux autres expériences et corrélations.
La tendance globale des expériences montre une lègère accumulation des résultats vers l’enve-
loppe supérieure du nuage de points. La corrélation de Lewis (97) est en bon accord avec résultats
expérimentaux. On note également un bon accord entre les deux simulations aux caractéristiques,
dites ”M.O.C.”(Method Of Characteristics) : Eastman & Radtke vs. Hanson & Palmer. La méthode
de calcul du disque de Mach n’est pas exactement la même. De plus, 38 années séparent ces deux si-
mulations et cette méthode est très sensible à la puissance de calcul mise en œuvre. Malgré cela, les
résultats de Eastman sont tout à fait corrects, ce qui prouve l’excellent rapport Précision/Coût Calcul
des méthodes M.O.C. Les méthodes volumes finis dites ”F.V.”(Finite Volumes) ont fait leurs preuves
dans les écoulements supersoniques. Ici, les simulations de Cumber & al. s’approchent de l’enve-
loppe supérieure des résultats précédents. De la même manière, des expériences, corrélations, et si-
mulations concernent le diamètre du disque de Mach. L’ensemble des sources menant à ces résultats
est listé de maniére non exhaustive dans le tableau 10.2, puis representé sur le graphique 10.4. Dans
le tableau, on voit une très grande disparité dans les corrélations, avec un exposant sur le paramétre
NPR (Exposant B) variant de 0.6 à 1.7. Le graphique inférieur de la fig Fig. 10.4 montre à nou-
veau une faible dispersion des résultats, donnant une erreur de 10% maximum, exception faite des
résustats de Crist, sur l’ensemble des rapports de pression. Les résultats expérimentaux de Crist
s’éloignent des autres résultats de plus de 30%. On notera cependant que cette série d’expériences
est la plus ancienne de l’ensemble bibliographique présenté ici concernant les tailles de disques
de Mach. Les simulations F.V. de Cumber & al. , ainsi que les simulations M.O.C. de Hanson et
Eastman se placent au milieu des résultats expérimentaux.
10.2. Taille des bouteilles de Mach 197
T. 10.2 – Compilation des différentes sources concernant le rapport d M /D. M.O.C : Méthodes aux
caractéristiques, F.V. Méthodes Volumes Finis
198 Chapitre 10. Annexe 2
10.0
Exp. APL
Exp. Addy
Exp. Crist
Exp. Ewan & Moodie 12.7 mm
Exp. Ewan & Moddie 25.4 mm
8.0 Exp. Love
Exp. NACA
Calc. Adamson
Corr. Lewis
M.O.C. Eastman
M.O.C. Hanson
F.V. Cumber
6.0
xM/D
4.0
2.0
0.0 NPR
0.0 20.0 40.0 60.0 80.0
Ptot/Poo
3.57 10 20 40 60 80 100 120
10.0
Exp. Antsupov
Exp. Crist
Exp. Ewan & Moodie 12.7 mm
Exp. Ewan & Moodie 25.4 mm
Corr. Antsupov
8.0 Corr. Lamkin
M.O.C. Hanson
F.V. Cumber
6.0
dM/D
4.0
2.0
0.0 NPR
0.0 20.0 40.0 60.0 80.0
Ptot/Poo
3.57 10 20 40 60 80 100 120
dv
Niveau 3
Niveau 2
Niveau 1
F. 10.5 – Jet de fluide modèle : Niveau 1, conditions dans le réservoir ; Niveau 2, conditions à
l’orifice ; Niveau 3, conditions après détente
En supposant d’abord que les forces visqueuses sont négligeables dans la zone de détente (du
niveau 2 au niveau 3), puis qu’il n’y a pas d’entraı̂nement du fluide ambient, on peut écrire la
conservation de la masse et du moment sous la forme des Eq. 10.5 et Eq. 10.6 :
ρ3 V3 A3 = ρ2 V2 A2C D (10.5)
ρ2 V22 A2C D
2
A3 = (10.7)
P2 − P3 + ρ2 V22C D
2
200 Chapitre 10. Annexe 2
P2 − P3
V3 = V2C D + (10.8)
ρ2 V 2 C D
Ce jeu d’équations a été établi par Hess & al. (1973). En utilisant les relations isentropiques
(Shapiro, 1953), on peut relier le niveau 1 au niveau 2 :
γ
! γ−1
2
P2 = P1 (10.9)
γ+1
! 1
2 γ−1 W
ρ2 = P1 (10.10)
γ+1 RT 1
s !
RT 1 2γ
V2 = (10.11)
W γ+1
Sachant que le fluide retrouve sa température totale après le disque de Mach, on a la ralation
T 3 ∼ T 1 . Ainsi, la vitesse V3 et le diamètre virtuel dv du jet s’expriment en fonction des grandeurs
du niveau 1 par les Eq. 10.12 et Eq. 10.13.
γ−1 !
−γ
P3 2
1 − P1 γ+1
V3 = V2 C D + (10.12)
γC D
s
1
! γ−1
dv P1 2 V2
= CD (10.13)
d P3 γ + 1 V3
Dans le cas des jets fortements sous-détendus, c.a.d. dans le cas de l’Eq. 10.14, l’equation
Eq. 10.13 se simplifie vers la forme de l’equation Eq. 10.15
! −γ
P3 2 γ−1
= NPR (10.14)
P1 γ+1
v
t ! 1
dv 2 γ−1 1
= NPR C D
1/2
(10.15)
d γ+1 γC D + 1
2
L’important dans cette relation est que le diamètre virtuel dv dépend du rapport de pression
totale NPR selon une loi du type : ddv ∝ NPR1/2 Ce résultat peut être comparé avec les différentes
corrélations donnant les tailles des bouteilles de Mach. En supposant un diamètre de disque de Mach
dM xM
d proportionnel à la longueur de la bouteille de Mach d et au diamètre virtuel du jet équivalent
dv
d , on obtient la relation de l’Eq. 10.16, qui est toutà fait compatible avec l’ensemble des résultats
expérimentaux.
dv d M xM
∝ ∝ ∝ NPR1/2 (10.16)
d d d
10.3. Analyse du champ lointain d’un jet sous-détendu 201
dv
U axe
Uo
1
o x ov
C axe 1
C0
0
0 x ov
F. 10.6 – Jet sonique sous-détendu (Trait plein) et son jet sonique adapté équivalent (Trait inter-
rompu). En champ lointain, les décroissances de vitesse et de concentrations sont hyperboliques.
202 Chapitre 10. Annexe 2
15.0
10.0
Vo / Vaxe
5.0
0.0
0.0 100.0 200.0 300.0 400.0
x/d
F. 10.7 – Données expérimentales de vitesse axiale pour un jet sonique fortement sous-détendu
(NPR=31). L’abcisse est normalisée par le diamètre de l’orifice d. La vitesse est normalisée par la
vitesse du son à la pression extérieure (zone détendue) et à la température totale du jet (température
du réservoir).
La Fig. 10.7 montre le profil de vitesse axiale d’un jet sonique fortement sous-détendu. Comme
c’est inverse de la vitesse qui est tracé, une évolution hyperbolique se traduit par une droite. On
observe donc une évolution hyperbolique pour x/d > 200. Dans la zone proche de l’orifice x ∈
[0, 200] apparaı̂t une forte déviation, qui marque l’éloignement d’un régime hypeerbolique.
L’Eq. 10.17 peut aussi s’ecrire sous la forme affine de l’Éq. 10.19. La Fig. 10.7 permet donc de
calculer à partir de résultats expérimentaux les inconnues xov et k.
U0 x d xov d
= − = αx + β (10.19)
Uaxe (x) d kdv d kdv
Le résultat du calcul des deux inconnues xo v et k à partir de résultats expérimentaux est présenté
dans la Fig. 10.8 pour des NPR allant de 1 à 80. Le facteur de déviation est toujours croissant. Par
contre, l’origine virtuelle xov tend vers une limite aux alentours de 50d, atteinte à 20% près dès
les NPR supérieurs à 20. Les Eq. 10.20 et Eq. 10.21 donnent des corrélations qui permettent de
retrouver ces évolutions.
xov
= 17.1971 + 7.43486 Ln(NPR) (10.21)
d
Les profils de concentrations et de vitesses axiales, une fois corrigés avec ces deux paramètres,
permettent la coalescence des résultats expérimentaux des jets sous-détendus, quelque soit le rapport
de pression NPR. La constante de déviation hyperbolique k ne dépend pas du rapport de pression.
On trouve expérimentalement k = 4.83, ce qui est proche de la valeur de cette constante pour les
jets subsoniques : k = 4.89.
10.3. Analyse du champ lointain d’un jet sous-détendu 203
50.0 50.0
40.0 40.0
30.0 30.0
Kdv/d
xov/d
20.0 20.0
10.0 10.0
0.0 0.0
0.0 20.0 40.0 60.0 80.0 0.0 20.0 40.0 60.0 80.0
NPR NPR
F. 10.8 – Evolution du facteur déviation k et de l’origine virtuelle xov en fonction du rapport de
pression NPR
30.0 20.0
15.0
20.0
Uo / Uaxe
Co / Caxe
10.0
10.0
5.0
0.0 0.0
0.0 50.0 100.0 150.0 0.0 50.0 100.0 150.0
(x−xov) / dv (x−xov) / dv
F. 10.9 – Données expérimentales de vitesse axiale et de concentration axiale pour 26 jets soniques
fortement sous-détendus (NPR ∈ [2.47, 75.92]). L’abcisse est normalisée par le diamètre virtuel du
jet. La vitesse est normalisée par la vitesse virtuelle V3 issue de l’équation Eq. 10.12 avec C D = 1.0
204 Chapitre 10. Annexe 2
Chapitre 11
205
206 Chapitre 11. Annexe 3
France ne sont pas mentionnés ici, car ils sont exclusivement dédiés aux prévisions. On notera
également que le l’IN2P3 se distingue par sa puissance de stockage, bientôt de l’ordre du Poctets,
bien plus que par la la puissance de calcul dégagée par ses fermes de calculateurs. Ces structures
sont présentées dans le Tableau 11.1.
T. 11.1 – Descriptif des structures accueillant les supercalculateurs ouverts au calcul scientifique
en France en Nov. 2005.
On trouve également des centres de calculs de taille intermédiaire, ou mésocentres. Ils mettent
en commun des moyens de calculs à l’échelon régional, dégageant des puissances de l’ordre de 0.1 à
1Tflops. Il y a enfin des moyens de calculs locaux installés pour les besoins specifiques d’une équipe
industrielle, militaire, ou de recherche. Le rapport ministériel 2005 fait remarquer que dans l’idéal
l’essentiel des calculs devraient êtres traités dans les mésocentres, tandis que les calculs de très
grande ampleur seraient installés sur les grands supercalculateurs. C’est pourquoi les rapporteurs
pointent du doigt le fractionnement des calculateurs. En effet, sur l’IBM-SP4 du CINES, 10% de la
production utilise au moins la moitié de la puissance de la machine (128 processeurs sur 256), 31%
utilise le quart (64 processeurs), 54% le huitième (32 processeurs). Inversement, plus de la moitié
de la production pourrait se faire sur une machine de 0.23 Tflops.
La conclusion fondamentale du rapport ministériel 2005 concerne le retard de la France dans
le paysage européen et mondial du calcul à hautes performances. Les statistiques de l’organisation
ARCADE3 montrent que la puissance de calcul de la France, même ramenée au produit intérieur
brut ou à la population, se situe aux toutes dernières places de l’union européenne. A l’échelle
mondiale, la liste du ’top 500’ est éloquente. Ce classement des 500 supercalculateurs les plus
puissants, établi par les universités de Mannhein et du Tenessee, est une référence mondiale. Dans
le classement de novembre 2005, la France n’y fait apparaı̂tre que 8 machines, dont la première de
situe à la 62e place.
T. 11.2 – Extrait du classement ’top 500’ des supercalculateurs en activité en Nov. 2005. La
puissance mentionnée est celle dégagé avec le banc d’essai numérique LINPACK
Pe f f = n p /ntot
p ∗ E ∗ S U(n p ) ∗ Pcrete (11.2)
4
Un code paralléle idéal S U = 1 demande une seconde d’exécution sur 10 processeurs lorsque dix secondes sont
nécessaires sur un seul.
Bibliographie
[1] A-H, K. S. Three-dimensional calculations for underexpanded and over expanded
supersonic jet flows. AIAA paper, 2196 (1989).
[3] A, T., N, J. On the structure of jets from highly underexpanded nozzles
into still air. Journal of Aerospace Science 26 (1959), 16–24.
[5] A, A. V. Properties of underexpanded and overexpanded supersonic gas jets. Soviet
Physics : Technical Physics 19 (1974), 234–238.
[6] A, S., S, M., E, G. On streamwise vortices in high reynolds number
supersonic axisymmetric jets. Phys. of Fluids 5, 1 (1993), 187–202.
[7] A, S. Ueber die reaktionsgeschwindigkeit bei der inversion von rohrzucker durch
saeuren. Z. Phys. Chem, 4 (1889), 226–248.
[8] B, T. A., H, T. A., L C, J.-F. Spray and self-ignition visualisation in a DI
diesel engine. SAE Paper 940681 (1994).
[9] B, R. A., A, G. A., H, H. A. Assumed joint probability density
function approach for supersonic turbulent combustion. J. Propul. Power 10, 4 (1994), 473–
483.
[10] B, R. A., A, G. A., H, H. A., D, J. P. An assumed joint-
beta PDF approach for supersonic turbulent combustion. In AIAA, SAE, ASME, and ASEE,
Joint Propulsion Conference and Exhibit, 28th, Nashville (1992).
[11] B, G., H, P., L, J. The proper orthogonal decomposition in the analysis
of turbulent flows. Ann. Rev. Fluid Mech. 25 (1993), 539–575.
[13] B, A. D., B, D. R., C, D. K., H, G. K. Flame stability in underex-
panded natural gas jets. Combus. Sci. and Tech. 58 (1988), 267.
[14] B, A. D., B, D. R., D, M. G., S, F. The structure and concentra-
tion decay of high pressure jets of natural gas. Combus. Sci. and Tech. 36 (1984), 249–261.
[15] B, A. D., H, D. J., S, F. Velocity decay of high pressure jets. Combus.
Sci. and Tech. 52 (1987), 161–171.
210 Chapitre 11. Annexe 3
[16] B, R., S, W., L, E. Transport phenomena. John Wiley, New York, 1960.
[17] B, C., B, C., J́, D. Noise investigation of a high subsonic, moderate reynolds
number jet using a compressible large eddy simulation. Theoretical and Computational Fluid
Dynamics 16, 4 (2003), 273–297.
[18] B, G. L., R, A. On density effects and large structure in turbulent mixing
layers. J. of Fluid Mech. 64 (1974), 775–816.
[19] B, M. C., K, A. P. An analytical and experimental study of supersonic
combustion of hydrogen in vitiated air stream. AIAA Journal 11, 9 (1973), 1217–1218.
[20] C, J. High resolution frequency-wavenumber spectrum analysis. In IEEE (1969), vol. 57,
pp. 1408–1418.
[21] C, I.-S. Mach Reflection, Mach Disc, and the Associated Nozzle Free Jet Flows. PhD
thesis, University of Illinois at Urbana-Champaign, 1945.
[22] C, D. Simulation numérique d’écoulements réactifs transsoniques. PhD thesis, Uni-
versité de Nice, 1991.
[23] C, M. J., D, C., D, J. J. Janaf thermochemical tables. third edition. i,
al–co. ii, cr–zr. J. Phys. Chem. 14 (1985).
[25] C, R. K., O, A. K. Autoignition in methane-hydrogen mixtures. Comb. and
Flame 58 (1984), 125–139.
[26] C, T. S., L, K. S. Numerical simulations of underexpanded supersonic jet and free
shear layer using WENO schemes. Int. J. Heat Fluid Fl. 26, 6 (2005), 755–770.
[27] C, T. S., W, J. A., P, R. W. Simultaneous temperature and multispecies
measurement in a lifted hydrogen diffusion flame. Comb. and Flame 91 (1992), 323–345.
[28] C, T. S., W, J. A., P, R. W., J, O., N, G. B. Raman mea-
surement of mixing and finite-rate chemistry in a supersonic hydrogen-air difusion flame.
Comb. and Flame 99 (1994), 157–173.
[29] C, S. G., L, M. C., F, G. M. Structure of turbulent sonic underexpanded free
jets. AIAA Journal 27 (1989), 549–559.
[30] C, N. T., M, M. G. Large-scale structure and entrainment in the supersonic
mixing layer. J. of Fluid Mech. 284 (1995), 171–216.
[31] C, O. Simulation aux Grandes Echelles de la Combustion Turbulente Prémélangée dans
les Statoréacteurs. PhD thesis, INP Toulouse, 2000.
BIBLIOGRAPHIE 211
[32] C, O., R, M. Development of high-order taylor-galerkin schemes for uns-
teady calculations. J. of Comp. Phys. 162, 2 (2000), 338–371.
[33] C, M.B., I., S, L. Scramjet fuels autoignition study. Journal of Propulsion
and Power 17, 2 (2001), 315–323.
[34] C, A., C, W. A high-wavenumber viscosity for high-resolution numerical me-
thods. J. of Comp. Phys. , 195 (VNR hyperviscosity 2004), 594–601.
[35] C, A., C, W. Hyperviscosity for shock-turbulence interactions. J. of Comp. Phys.
, 203 (2005), 379–385.
[36] C, R., C, J. A. Chemical aspects of the autoignition of hydrocarbon-air mixtures.
Comb. and Flame 60 (1985), 109–123.
[37] C, S., G, D. R., S, P. M. Study of the highly underexpanded sonic jet.
AIAA Journal 4 (1966), 68–71.
[38] C, P., F, M., F, S., G, J. Predictions of the structure of
turbulent, highly underexpanded jets. Journal of Fluids Engineering 117 (1995), 599–604.
[39] C, E. T., M, S., Eds. Scramjet propulsion. No. ISBN 1563473224 in V-189.
AIAA American Institute of Aeronautics and Astronautics, 2000.
[40] C, T. Scramjet engines : The first forty years. J. Propul. Power 17, 6 (2001), 1138–
1148.
[41] C, A. D., W, J. A. An experimental study of a supersonic coaxial jet. AIAA
Paper, 143 (2001).
[42] C, A. P., B, T. A., P, T. Self-ignition and combustion modeling of
initially nonpremixed turbulent systems. Comb. and Flame 124, 1 (2001), 65–81.
[43] D, S. M., W, D. E. Analysis of turbulent underexpanded jets. i - parabolized
navier-stokes model, SCIPVIS. AIAA Journal 23, 4 (1985), 505–514.
[44] D, M., R, W. C., M, N. N. The structure of the compressible reacting
mixing layer : Insights from linear stability analysis. Phys. of Fluids 10, 4 (1998), 993–1007.
[45] D R, A., I, A., B, C., P, V., G, E. Les of supersonic
combustion of H2/vitiated air. In 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference
and Exhibit ; Fort Lauderdale (2004).
[46] D, C. B., K, J. B., M, J. B., S, C. D. Stability of underexpanded
supersonic jet flames burning H2-CO mixtures. Shock Waves 12 (2002), 241–249.
[47] D, W. S., G, D. Spectral simulation of supersonic reactive flows. SIAM J.
Numer. Anal 35, 6 (1998), 2370–2384.
212 Chapitre 11. Annexe 3
[48] D, J., H, H., Y, Y., D, J. M. Measured lenghts of supersonic hydrogen-
air jet flames compared to subsonic flame lengths and analysis. Comb. and Flame 107 (1996),
176–186.
[50] D, J. P., R, R. C., H, M. Y. A detailed numerical model of a super-
sonic reacting mixing layer. In AIAA, ASME, SAE, and ASEE, Joint Propulsion Conference,
22nd, Huntsville (1986).
[51] D, F., C, P., L, M. Large-eddy simulation of transition to turbulence in
a boundary layer developing spatially over a flat plate layer developing spatially over a flat
plate. J. of Fluid Mech. 326 (1996), 1–36.
[52] E, D., R, P. Location of the normal shock wave in the exhaust plume of a
jet. AIAA Journal 1 (1963), 918–919.
[53] E, T., C, J. Direct numerical simulation of autoignition in non-homogeneous
hydrogen-air mixtures. Comb. and Flame 134, 3 (2003), 169–191.
[54] E, D. R., B, R. A., G, M. R. Numerical study of a scramjet combustor
fueled by an aerodynamic ramp injector in dual-mode combustion. In AIAA, Aerospace
Sciences Meeting and Exhibit, 39th, Reno (2001).
[55] E, J. S., S J., C. J. Influence of chemical kinetics and unmixedness on
burning in supersonic hydrogen flames. AIAA Journal 18, 2 (1980), 188–193.
[56] E, J. S., S J., C. J., B, H. Application of a two-dimensional para-
bolic computer program to prediction of turbulent reacting flows. Tech. Rep. 1169, NASA,
1978.
[57] E, B., M, K. Structure and velocity measurements in underexpanded jets. Com-
bus. Sci. and Tech. 45 (1986), 275–288.
[59] F, A. General Teory of High Speed Aerodynamics. Princeton, NJ, 1954, ch. Supersonic
Flows with Shock Waves.
[62] F, B., K., L. S. Numerical investigation of a transverse jet in a supersonic crossflow
using large eddy simulation. AIAA Paper, 2006-3712 (2006).
BIBLIOGRAPHIE 213
[63] F, J. B., K., L. S., M, P. Numerical simulation of a mach 1.92 turbulent jet and
its sound field. AIAA Journal 38, 11 (2000), 2023–2031.
[64] G, F., M, S. Les of supersonic combustion of hydrocarbon spray in a scramjet.
In 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit ; Fort Lauderdale
(2004).
[65] G, P., S, P., B̈, D. An implicit multigrid method for the simulation
of chemically reacting flows. J. of Comp. Phys. 146, 1 (1998), 322–345.
[68] H, M., K, L., Q, C. The autoignition of hydrocarbon fuels at high
temperatures and pressures-fitting of a mathematical model. Comb. and Flame 30 (1977),
45–60.
[69] H, R., T́, D. Autoignition of turbulent non-premixed flames investigated
using direct numerical simulations. Comb. and Flame 128, 1/2 (2002), 22–37.
[70] H, D., P, D. Hybrid tuned center-difference-weno method for large eddy simula-
tions in the presence of strong shocks. J. Comput. Phys. 194, 2 (2004), 435–450.
[71] H, C. Numerical Computation of Internal and External Flows. John Wiley & Sons,
1989.
[72] H, J., C, C., B, R. Molecular theory of gases and liquids. John
Wiley & Sons, New York, 1964.
[73] H, J., W, A., M, P. Eddies, streams, and convergence zones in turbulent flows.
In Summer Program (1988), vol. SEE N89-24538 18-34, NASA, pp. 193–208.
[74] I, H., C, J., L, C. Ignition of hydrogen-air mixing layer in turbulent flows. In
27th Symp. (Int.) on Combustion (1998), pp. 1047–1056.
[75] J, C. J. An analytical study of the hydrogen-air reaction mechanism with applica-
tion to. scramjet combustion. Tech. Rep. TP-2791, NASA, 1988.
[76] J, A., S, W., T, E. Numerical solution of the euler squations by
finite volume methods using runge-kutta time stepping schemes. In 14th Fluid and Plasma
Dynamic Conference (Palo Alto, 1981), A. p. 81-1259, Ed.
[77] J, I.-S., C, J.-Y. Numerical simulation of supersonic combustion for hypersonic
propulsion. In 5th Asia-Pacific Conference on Combustion ; Adelaide, Australia (2005).
[78] J, G.-S., S, C.-W. Efficient implementation of weighted eno schemes. J. Comput.
Phys. 126, 202-226 (1996).
214 Chapitre 11. Annexe 3
[79] J, W., K, M. Pdf modeling of finite-rate chemistry effects in turbulent nonpre-
mixed jet flames - a general calculation procedure. Comb. Flame 115, 1 (1998), 210–229.
[80] K, T., H, T., M, T., T, S., C, N. Mach 6 testing of a scramjet
engine model. J. Propul. Power 13, 4 (1997), 543–551.
[81] K, R., R, F., M, J. CHEMKIN II, a fortran chemical kinetics package for the
analysis of gas-phase chemical kinetics. Tech. Rep. SAND89-8009, Sandia, 1989.
[82] K, R., D, A., C, B., P, T. Comparison of computational me-
thodologies for ignition of diffusion layers. Combus. Sci. and Tech. 175, 10 (2003), 1783 –
1806.
[83] K, M., N, K., H, T., K, T., M, T. Scramjet inlet flow
computations by hybrid grid method. In AIAA, Aerospace Sciences Meeting and Exhibit,
36th, Reno (1998).
[84] K, M., S, T., N, K. Numerical analysis of scramjet combusting flows
by unstructured hybrid grid method. In AIAA, Aerospace Sciences Meeting and Exhibit, 38th,
Reno (2000).
[85] K, A. The local structure of turbulence in incompressible viscous fluid for very
large reynolds numbers. C.R. Acad. Sci. USSR 30 (1941), 301.
[86] K, S., H, Z., R, R. The development and application of a diesel ignition and
combustion model for multidimensional engine simulations. SAE paper 950278 (1995).
[87] K, S., M, C., R, R. Modeling and experiments of HCCI engine combustion
using detailled chemical kinetics with multidimensional CFD. SAE paper 2001-01-1026
(2001).
[88] K, A., S, R. Computation of supersonic mixing layers. Phys. of Fluids 14,
11 (2002), 3790–3797.
[89] K, A. Numerical simulation of scramjet inlet flow fields. Tech. Rep. TP-2517, NASA,
1986.
[91] L, L. D., L, E. M. Course of theoretical physics : Fluid mechanics. Oxford :
Pergamon Press, 1959.
[92] L̂, L. Simulation des grandes échelles de l’écoulement au-dessus d’une cavité.
PhD thesis, Université Paris 6, 2003.
[93] L̂, L., S, P., L̂, T.-H., C, P. Large-eddy simulation of a compres-
sible flow in a three-dimensional open cavity at high reynolds number. J. Fluid Mech. 516
(2004), 265–301.
BIBLIOGRAPHIE 215
[94] L̂, L., S, P., M, I., L, O. Large-eddy simulation of a compressible
flow past a deep cavity. Phys. Fluids 15, 1 (2003), 193–210.
[95] L, S., L, S., M, P. Interaction of isotropic turbulence with shock waves : effect
of shock strength. J. of Fluid Mech. 340 (1997), 225–247.
[96] L, S. Compressibility effects on turbulence. Ann. Rev. Fluid Mech. 26 (1994), 211–254.
[97] L, C., C, D. Normal shock location in underexpanded gas and gas-particle
jets. AIAA Journal 2, 4 (1964), 776–777.
[98] L, J., Y, G., Z, Y., L, Y., Q, D. Experimental studies on self-ignition of hydro-
gen/air supersonic combustion. Journal of Propulsion and Power 13, 4 (1997), 538–542.
[99] L, E, S., G, C. E., L, L. P., W, M. J. Experimental and theoretical
studies of axisymmetric free jets. NASA, 19980228067 (1959).
[101] M, R. R., H, M. E., P, L. A. Structure of supersonic jet flow and its
radiated sound. AIAA Journal 32, 5 (1994), 897–906.
[102] M, E., B, T., P, T. J. Numerical simulations of autoignition in
turbulent mixing flows. Comb. and Flame 109, 1/2 (1997), 198–223.
[103] M, E., C, A. P., B, T. A., P, T. A model for the effects
of mixing on the autoignition of turbulent flows. Combus. Sci. and Tech. 125, 1/6 (1997),
243–282.
[104] M, Z. Instabilité linéaire d’un jet rond supersonique. Master’s thesis, Ecole Nationale
Supérieure des Arts et Métiers de Paris, France, 2005.
[106] M, H., G, P., B, D. Scalar and joint scalar-velocity-frequency
monte carlo pdf simulation of supersonic combustion. Comb. Flame 132, 1 (2003), 3–24.
[107] M, P., S, C.-W. Real gas computation using an energy relaxation method and
high-order weno schemes. ICASE Report, 98-42 (1998).
[108] M, C. J., Z, W., A, B. R. Supersonic combustion simulations using
reduced chemical kinetics mechanisms and ISAT. In 24th Fluid Dynamics meeting, Orlando
(2003), vol. AIAA Paper 2003-3547.
[109] N, F., D, F. Subgrid-scale stress modelling based on the square of the velocity
gradient. Flow Turbul. Combust. 62, 3 (1999), 183–200.
216 Chapitre 11. Annexe 3
[110] N, J. W., E, J. R. Large-eddy simulation of high-speed, turbulent diffusion
flames with detailed chemistry. In AIAA, Aerospace Sciences Meeting and Exhibit, 35th,
Reno (1997).
[111] O, R., C, G., C, J., H, M. Companion to the History of Modern
Science. Routledge, 1996, ch. Science and War, by D.E.H. Edgarton, pp. 934–945.
[112] P, J., H, R. Application of method of characteristics to underexpanded, freejet
flows with vibrational nonequilibrium. AIAA Journal 36, 2 (1998), 193–200.
[113] P, J. On the origin of rockets. T’oung Pao LXXIII (1987), 2–15.
[114] P, J., S, R. G. Measurement of shock structure and shock–vortex interaction
in underexpanded jets using rayleigh scattering. Phys. of Fluids 11, 12 (1999), 3761–3777.
[115] P, D., B, A. Evolution of large eddies in compressible shear
layers. Phys. of Fluids 9, 3 (1997), 756–765.
[116] P, U., F, K., A, G. Self-ignition of diesel-relevant hydrocarbon-air
mixtures under engine conditions. In 26th Symp. (Int.) on Combustion (1996), p. 781.
[117] P, T., A, C., E, F., V, D. Large eddy simulations of
combustion instabilities. In 1st Int. Symp. On Turbulence and Shear Flow Phenomena (Santa
Barbara, Sept 12-15., 1999), pp. 1–6.
[118] P, T., L, S. Boundary conditions for direct simulations of compressible viscous
flows. J. of Comp. Phys. 101 (1992), 104–129.
[119] P, T., V, D. Theoretical and numerical combustion, second edition. R.T.
Edwards, 2005.
[120] P, S. B. A monte carlo method for the pdf equations of turbulent reactive flow. Combus.
Sci. and Tech. 25 (1981), 159–174.
[122] R, A., D, J., H, H., B, R. A. Combustion efficiencies of supersonic
flames. J. Propul. Power 17, 2 (2001), 301–307.
[123] R, L. Analytic solutions of the rayleigh equation for linear density profiles. In Proc.
London. Math. Soc. (1883).
[124] R, D. Analysis of growth rates in three-dimensional, air-to-air supersonic shear layers
using direct numerical simulation. In AIAA, Aerospace Sciences Meeting and Exhibit, 33rd,
Reno (1995).
[125] R, K. E., V, R. T., R, R. C., H, L. D. NASA’s hyper-X scramjet
engine ground test program. In ISABE - International Symposium on Air Breathing Engines,
14th, Florence, Italy (1999).
BIBLIOGRAPHIE 217
[126] R, P. Approximate riemann solvers, parameter vectors, and difference schemes. J. Comput.
Phys. 43 (1981), 357–372.
[127] R, C. J., E, J. R. Numerical simulation of a three-dimensional flame/shock wave
interaction. AIAA Journal 38, 5 (2000), 745–754.
[128] S, P. Introduction à la simulation des grandes échelles pour les écoulements fluides
incompressibles. Springer, 1998.
[129] S, N. D., R, W. C. Three-dimensional simulations of large eddies in the
compressible mixing layer. J. of Fluid Mech. 224 (1991), 133–158.
[130] S̈, T., P, T. Initial and boundary conditions for large eddy simulation of
combustion instabilities. In Annual Research Briefs. Center for Turbulence Research, NASA
Ames/Stanford Univ., 1999, pp. 73–84.
[131] S, J. M., N, T. D. Aerodynamic aspects of shock containing jet plumes. AIAA
Paper, 90-0217 (1980).
[132] S, A. The dynamics and thermodynamics of compressible fluid flow. New York :
Ronald Press, 1953.
[133] S, H., T, C. K. W. Numerical simulation of the generation of axisymmetric mode
jet screech tones. AIAA Journal 36, 10 (1998), 1801–1807.
[134] S, J., H, E., H, A. A History of Technology. London : Oxford U. Press,
1958.
[135] S, J. General circulation experiments with the primitive equations : 1. the basic
experiment. Monthly Weather Review 91, 3 (1963), 99–164.
[136] S, L., C, M.B., I. Ignition delay characteristics of methane fuels. Progress
in Energy and Combustion Science 20 (1994), 431–460.
[137] S, L., TV, J. Autoignition characteristics of aircraft-type fuels. Comb. and
Flame 46 (1982), 283–300.
[138] S, R. A multi-dimensional upwind discretization method for the Euler equations on
unstructured grids. PhD thesis, TU Delft, 1994.
[139] S, D., P, H. Janaf thermochemical tables, 2nd edition. Tech. Rep. NSRDS-
NBS 37, US National Bureau of Standards, 1971.
[140] S, C. J., L, J., Y, G., L, C. K. Chemical kinetics and self-ignition in a model
supersonic hydrogen–air combustor. AIAA Journal 37, 2 (1999), 208–214.
[141] S, F. Ignition of gaseous hydrocarbon fuels in hypersonic ramjets. In ISABE, Marseille
(1972).
218 Chapitre 11. Annexe 3
[142] T, C. K. W. Supersonic jet noise. Ann. Rev. Fluid Mech. 27 (1995), 17–43.
[144] T, P. A computation fluid dynamic model of a supersonic axi-symmetric jet using a beta
probability function combustion model. Tech. rep., Departement of Mechanical Engineering,
Purdue University, 2002.
[145] V L, B. Towards the ultimate conservative difference scheme v. J. Comput. Phys. 32
(1979), 101–136.
[146] V, R., C, J. H., P, R. W. Modeling ideally expanded supersonic turbulent
jet flows with nonpremixed H2-air combustion. AIAA Journal 30, 2 (1992), 395–402.
[147] L, E., K, M. High resolution simulation of supersonic combustion.
In ASME, SAE, and ASEE, Joint Propulsion Conference and Exhibit, 32nd, Lake Buena Vista
(1996).
[148] VN, J., R, R. A method for the numerical calculation of hydrodyna-
mic shocks. Journal of Applied Physics 21(3) (1950), 232–237.
[149] V, P. Modeling and simulation of transport and combustion phenomena in a su-
personic mixing layer. PhD thesis, Ecole Centrale de Paris, 1992.
[151] W, C., P, H., P, N. Simulation of autoignition delay and location of fuel
sprays under diesel-engine relevant conditions. SAE paper 97590 (1997).
[152] W, J. A., M, J. A pseudo-temporal multi-grid relaxation scheme for solving
the parabolized navier-stokes equations. AIAA Paper, 99-3360 (1999).
[153] W, G., MC, R. Modeling supersonic combustion using a fully implicit
numerical method. AIAA Journal 30, 4 (1992), 1008–1015.
[154] X, J., P, S. B. Pdf calculations of turbulent nonpremixed flames with local extinction.
Comb. Flame 123 (2000), 281–307.
[155] Y, R., D, F., R, H. A comprehensive reaction mechanism for applications
in combustion systems. Combus. Sci. and Tech. 79 (1991), 97–128.
[156] Y, B., L, K., L, M., M, M., B, R., D, R. Visualization of
a supersonic underexpanded jet by planar rayleigh scattering. Phys. Fluids A 1, 9 (1989),
1449.
[157] Y, Y., D, J. M., D, J. F. Blowout stability limits of a hydrogen jet flame in
a supersonic, heated, coflowing air stream. Combus. Sci. and Tech. 97, 1/3 (1994), 137–156.
BIBLIOGRAPHIE 219
[158] Y, A. Statistical theory for compressible turbulent shear flows, with the application
to subgrid modeling. Phys. of Fluids 29, 7 (1986), 2152–2164.
[159] Z, V., K, V. I., L, A. V. Experimental investigation of shock wave
structure and streamwise vortices in supersonic jet. In 4th Symp. on Aerothermodynamics for
Space Vehicles (2001), pp. 603–610.
[160] Z, V., S, A. Spatial structure of flow in the initial section of a super-
sonic underexpanded jet. Preprint no. 23-88, udk 533.6.011, Academy of Sciences USSR,
Siberian section, Institute of Theoretical and Applied Mechanics, 1988.
[161] Z, V. I., K, N. P., P, A. A. Effect of streamline curvature on intensity
of streamwise vortices in the mixing layer of supersonic jets. J. of Applied Mechanics and
Technical Physics(Russia) 45, 3 (2004), 335–343.
[162] Z, Y.-T., S, C.-W. High-order weno schemes for hamilton-jacobi equations on
triangular meshes. SIAM J. Sci. Comput. 24, 3 (2003), 1005–1030.
[163] Z, Z. C., Y, S. T. Shock capturing without riemann solver- a modified space-time
CE/SE method for conservation laws. AIAA paper, 99-0904 (1999).
[164] Z, L. L., B, K. N. C. The application of new combustion and turbulence models
to H2 - air non-premixed supersonic combustion. Comb. and Flame 99 (1994), 440–448.
INSTITUT NATIONAL POLYTECHNIQUE DE TOULOUSE
Doctorat d’université, spécialité dynamique des fluides
16 juin 2006
Antoine Dauptain
Résumé
Les lanceurs spatiaux ont aujourd’hui besoin de moteurs-fusées cryotechniques capables de s’al-
lumer plusieurs fois au cours du même vol. L’allumage étant un mécanisme extrêmement délicat
dans les conditions en vol, il est nécessaire de développer et d’utiliser des outils précis et fiables
pour aider au développement de cette technologie. Cette thèse développe la simulation des grandes
échelles (SGE) pour traiter les écoulements transitoires supersoniques réactifs. Différents aspects
sont abordés : cinétique chimique de l’auto-allumage et diffusion différentielle, traitement numérique
des écoulements supersoniques, combustion. Des comparaisons avec des données expérimentales
sur des configurations académiques permettent de valider les développements réalisés et de com-
prendre en détail le mécanisme d’auto-allumage. Sur la base de ces résultats, des simulations SGE
de configurations industrielles sont maintenant envisageables, afin d’étudier les différents régimes
transitoires d’allumage.
Mots clefs : Simulation des grandes échelles, combustion, auto-allumage, supersonique
Abstract
Today, space launchers require cryotechnic rocket engines able to reignite during flight. The ig-
nition phases in flight conditions are particularly critical and the development of restartable engines
needs accurate and reliable tools. The present thesis develops a Large Eddy Simulation approach
(LES) for the study of unsteady supersonic reactive flows. Several aspects are treated : chemical
kinetics, auto-ignition and differential diffusion, numerical methods suited to supersonic flows and
their discontinuities, combustion. Comparisons with experimental data on academic test cases va-
lidate the models, and give detailed insights into the auto-ignition process. Based on these achie-
vements, LES of industrial configurations may be now envisaged, allowing the study of unsteady
ignition regimes and the optimization of devices.
Key words : Large eddy simulation, combustion, auto-ignition, supersonic
CERFACS
42, Avenue Gaspard Coriolis. 31057 Toulouse Cedex 01. France.