Modélisation LES de la combustion turbulente

Télécharger le fichier original (Mémoire de fin d’études)

Variabilité de la combustion, transitoires et combustions anormales

Les moteurs à allumage commandé n’effectuent pas des cycles parfaitement stables [2]. Ces dif-férences d’un cycle à l’autre sont le résultat de variations dans le processus de combustion, c’est ce que l’on appelle la variabilité cyclique de la combustion (VCC). Ces variations se manifestent entre autre par des variations de la pression dans le cylindre d’un cycle à l’autre comme sur la figure 1.2. Trois points de fonctionnement sont présentés dont deux présentent des enveloppes de pression bien plus larges que le premier, cela se traduit par un écart-type plus ou moins important sur la PMI. On définit un moteur stable comme présentant un écart-type sur la PMI inférieur à 5%. Ces variations peuvent être à l’origine de complications supplémentaires, en effet les cycles dans la partie haute de l’enveloppe présentent un risque d’auto-inflammation plus important. Sur la figure 1.3 un cycle en cliquetis est mis en avant, il correspond à une auto inflammation du mélange qui génère une résonance dans la chambre. Ces cycles limitent les propriétés d’inflam-mabilité du carburant (indice d’octane) ainsi que le taux de compression du moteur. À l’opposé les cycles du bas de l’enveloppe risquent des combustions incomplètes très pénalisantes en terme d’émissions voire de ratés d’allumage. Lorsque l’auto inflammation du mélange se produit avant l’allumage par bougie et pendant la phase de compression on parle de préallumage, cette com-bustion anormale doit absolument être évitée car elle cause d’importants dégâts au moteur.

La simulation aux grandes échelles dans les moteurs

L’objectif de cette thèse est l’utilisation de la simulation aux grandes échelles pour le calcul de l’écoulement dans la chambre de combustion du moteur. C’est une technique qui a déjà démontré son efficacité à simuler les phénomènes acycliques présents dans les moteurs [5]. Le code AVBP a déjà été utilisé dans de nombreuses études sur la VCC dans les moteurs [6, 7] qui ont démontré l’intérêt de la SGE pour simuler la VCC et ont mis en évidence le manque de données pour la définition des conditions aux limites du système. La base de données SGEmac [3] dédiée à l’étude de l’utilisation de la SGE dans les moteurs à allumage commandé a donné naissance à une série d’études [8, 9] pour lesquelles un grand nombre de données expérimentales permet la validation des calculs que ce soit dans les lignes d’admission et d’échappement ou dans le cylindre. Dans toutes ces études l’intégralité du banc moteur a été simulée, ce que la complexité d’un moteur industriel rend impossible. Dans l’optique de simuler uniquement le cylindre en 3D, une tentative de couplage entre un code 1D (Gasdyn [10, 11]) et AVBP n’a pu être menée à son terme. Finalement des résultats de simulations système ont été utilisés pour définir les conditions aux limites [12] mais cette solution ne permet pas de prendre pleinement en compte les interactions et l’impact des phénomènes acycliques sur le système entier ne peut pas être simulé. Une étude sur les combustions anormales [13] a été menée sur un autre moteur académique dont la configuration est très proche de celle du système simulé dans la deuxième partie de nos travaux, chapitre 5. L’utilisation du code AVBP offre ainsi un choix de modèles adaptés au formalisme SGE comme les modèles de combustion CFM-LES [14] et TFLES [15] ou le modèle d’allumage par bougie ISSIM-LES [16] qui sont tous validés sur des simulations de moteurs.
Les groupes motopropulseurs actuels sont des ensembles de sous-systèmes extrêmement complexes qui comportent beaucoup d’éléments en interaction. Les turbocompresseurs relient l’écoulement à l’échappement à celui de l’admission par un lien mécanique et des parties tour-nantes ; de la même manière le système EGR lie l’échappement à l’admission par une vanne EGR pilotée par le calculateur. La bougie et les injecteurs sont aussi pilotés par le calculateur, permettant un contrôle dynamique du système moteur. Par exemple, retarder le moment de l’allu-mage permet d’éviter le phénomène de cliquetis et un capteur sur le moteur permet d’asservir le système. Des techniques similaires existent pour les injecteurs : un surplus d’essence est injecté pour les fortes charges et permet de réduire la température.
Le grand nombre d’éléments, de capteurs et de leviers pour optimiser tel ou tel comporte-ment du moteur rendent la simulation système indispensable que ce soit pour la simulation ou le contrôle dynamique. Dans cette optique le couplage entre la simulation système et la SGE per-mettrait d’obtenir un outil de calcul précis permettant des configurations proches des systèmes réels.

Couplage de codes

Il convient de réaliser un couplage entre deux codes de dynamique des fluides qui utilisent des équations différentes. L’un se base sur les équations d’Euler, l’autre de Navier-Stokes avec des réactions chimiques. L’existence de solutions d’un tel système est démontrée dans le cas Eu-ler/Euler avec des lois d’état différentes [17] mais pas dans notre cas. De plus les schémas nu-mériques qui permettent d’obtenir une approximation de la solution sont différents. Enfin on aura forcément des pertes d’information au niveau de l’interface de couplage lorsqu’on passe d’une dimension à trois dimensions ou qu’on réduit le nombre d’espèces en créant des espèces composées. Le positionnement de l’interface dans un système peut permettre de réduire la perte d’information [18], en effet il parait évident que la solution sera moins perturbée par une interface située dans une zone où l’écoulement est quasi monodimensionnel que dans une zone de recir-culation. Dans la littérature, on retrouve beaucoup de couplages RANS-SGE en 3D, quelques couplages 1D-RANS et très peu de couplages 1D-SGE. Au niveau des méthodes employées, on peut caractériser quatre grandes tendances, l’utilisation de la méthode de Godunov pour mettre à jour les grandeurs à l’interface [11] (utilisation d’un solveur de Riemann), le passage par la méthode des caractéristiques [19], le chevauchement des domaines à coupler qui permet à l’in-terface d’être entièrement mise à jour par les schémas numériques des codes [20–22] ou encore des méthodes itératives [23].

Simulation d’écoulements réactifs Codes AVBP et Flow1D

Ce chapitre présente les modèles utilisés pour la simulation d’écoulements gazeux réactifs et les équations qui en découlent. Les deux codes utilisés au cours de cette thèse y sont présentés : AVBP et Flow1D. Le code AVBP est utilisé pour la simulation aux grandes échelles (LES) tridi-mensionnelle. Le code Flow1D permet de simuler des écoulements monodimensionnels et a été développé pendant cette thèse afin de tester et utiliser les méthodes de couplage avec AVBP. Les équations de Navier-Stokes sont d’abord présentées sous leur forme générale. Les différents mo-dèles et hypothèses concernant la thermodynamique et la diffusion sont introduits pour les codes AVBP et Flow1D. Le principe de la LES est ensuite présenté. Enfin les méthodes numériques associées aux deux codes ainsi que les modèles particuliers à la combustion turbulente et à la modélisation monodimensionnelle de systèmes moteurs sont décrits.

Description numérique du code AVBP

Discrétisation, méthode Cell-Vertex

Le code AVBP a été développé de manière à utiliser des maillages non-structurés (la compo-sition des éléments et les positions de ses nœuds ne peuvent pas être déduites d’une règle de construction) et hybrides (composés de différents types d’éléments). Ce type de maillage permet une grande flexibilité pour la discrétisation de géométries complexes et de choisir le type d’élé-ment le plus adapté au mouvement de maillage. La thèse de Enaux [32] montre que les tétraèdres sont les éléments les plus adaptés à la simulation d’écoulement moteur avec le code AVBP. On commence par discrétiser le domaine fluide en éléments C appelés cellules ou mailles.
Th étant la triangulation du domaine . Les méthodes numériques implantées dans AVBP sont des méthodes volumes finis de type Cell-Vertex. Cette orientation a démontré une meilleure ro-bustesse aux déformations de maillage que les méthodes Cell-Centered ou Vertex-Centered dans les travaux de Rudgyard [33] et Morton et al. [34]. Les méthodes centrées (Centered) consistent à stocker les variables aux centres des volumes de contrôle : en Cell-Centered le volume de contrôle est cellule primale et les variables sont stockées aux centres des mailles ; en Vertex-Centered le volume de contrôle est la cellule duale et les variables sont stockées aux nœuds du maillage. Dans le cas Cell-Vertex les variables sont stockées aux nœuds du maillage alors que le volume de contrôle est la maille primale. Cette technique nécessite donc de répartir l’information des nœuds vers les centres de maille et réciproquement.
Le schéma d’un maillage 2D est présenté sur la figure 2.1. On utilise l’indice j pour désigner un nœud du maillage, C une maille primale et Cj la maille duale associée au nœud j. On appelle Dj l’ensemble des cellules primales qui contiennent le nœud j.
En intégrant l’équation 2.4 sur un volume de contrôle C de normale sortante ~n, en négligeant les termes sources et en y appliquant le théorème de Green-Gauss on peut écrire : @ ZC W dv + Z@C F:~nds = 0 (2.47)
La première étape de la méthode Cell-Vertex consiste à calculer le résidu aux cellules qui corres-pond à la divergence du flux de l’équation 2.4 : 1 Z@C Fh:~nds (2.48) rC = jCj où jCj est le volume de la maille C et Fh l’approximation numérique du flux F . Comme les variables sont stockées aux nœuds, pour calculer ce résidu il est nécessaire de définir l’approxi-mation numérique du flux Fh et les normales aux nœuds. La normale nj;C au nœud j associée à la cellule C est une combinaison linéaire des normales aux faces adjacentes à j dans C.

Schéma Lax-Wendroff

Le schéma Lax-Wendroff [35, 36] combine une discrétisation spatiale et temporelle pour donner un schéma volume fini centré d’ordre 2 en temps et en espace. Comme tous les schémas centrés d’ordre 2, il produit des oscillations parasites dans les zones de fort gradient. Les méthodes de stabilisation utilisées sont décrites au chapitre 2.2.3. Il s’écrit dans la formulation Cell-Vertex Wjn+1 = Wjn t X Dj Dj;CLW rCn jCj (2.53)

Stabilisation : Viscosité artificielle

Les schémas numériques pour la convection centrés d’ordre 2 en espace génèrent des oscilla-tions dans les zones présentant de forts gradients. Ces oscillations n’ont aucun sens physique et posent des problèmes de stabilité aux schémas numériques, il est donc nécessaire de les éliminer. Dans AVBP ce sont des viscosités artificielles appliquées dans les zones de forts gradients qui permettent de stabiliser les schémas. Une itération du schéma numérique s’écrit : Wjn+1 = Wjn t rj + dj(2) + dj(4) (2.55)
L’équation 2.58 traduit le fait que cet opérateur est appliqué partout où la viscosité artificielle d’ordre 2 ne l’est pas et est ajustable par l’utilisateur par le biais de la variable (4).
Les senseurs C ont été introduits par Jameson [38]. Ils sont calculés pour chaque cellule et leurs valeurs sont comprises entre 0 et 1. Leur objectif est de n’être non nul que dans les zones où des problèmes de stabilité numérique apparaissent. Deux senseurs sont disponibles dans AVBP, celui de Jameson [38] développé pour des simulations stationnaires et celui de Colin [15] qui est une adaptation moins dissipative du précédent pour les calculs de LES transitoires.
Dans AVBP on appelle modèle de viscosité artificielle la combinaison des variables à partir desquelles les valeurs des senseurs sont calculées et les équations dans lesquelles sont appliqués les opérateurs de viscosité artificielle. Le modèle associé au senseur de Colin ayant pour but d’effectuer des calculs LES, utilise l’énergie totale et les fractions massiques pour calculer le senseur. Les viscosités artificielles sont appliquées dans les équations de conservation de la masse et de l’énergie mais seule la viscosité artificielle d’ordre 2 est appliquée à l’équation de quantité de mouvement. Ainsi aucune dissipation de fond n’est introduite dans cette équation.

Gestion du maillage mobile : ALE et ITC

Le domaine fluide dans un moteur à combustion interne présente de fortes déformations : le vo-lume de gaz dans le cylindre varie d’un facteur 10 (taux de compression) au cours du cycle et les soupapes passent d’un position fermée à une levée d’environ 8 mm. Les méthodes numé-riques utilisées doivent être adaptées et assez robustes pour traiter de telles déformations. En mécanique des milieux continus il y a deux représentations mathématiques de l’écoulement d’un fluide : la description Lagrangienne où l’on s’intéresse aux trajectoires de particules individuelles et la description Eulérienne où l’on s’intéresse aux grandeurs du fluide en un point fixe de l’es-pace (référentiel du laboratoire). C’est le formalisme généralement utilisé en CFD ; toutes les équations du présent chapitre utilisent cette description.
Par analogie on parle de méthode de mouvement de maillage Lagrangienne quand les nœuds suivent le mouvement des points matériels et Eulérien quand les mailles sont fixes. La mé-thode ALE (Arbitrary Langragian Eulerian) apporte une correction au schéma numérique afin de prendre en compte le mouvement des nœuds du maillage. Elle permet au sein d’un même maillage d’avoir des mailles fixes (Eulériennes), des mailles se déplaçant à la vitesse du mi-lieu (Lagrangiennes) et des mailles avec un déplacement arbitraire. Initialement proposée par Noh [39] et Hirt et al. [40] elle a pour but de combiner les avantages de ces deux méthodes en minimisant leurs inconvénients. On obtient ainsi un champ de déplacement lisse pondéré par les valeurs aux parois présenté sur la figure 2.2(a). Plus le maillage est déformé plus sa qualité est faible, lorsque la qualité du maillage est trop faible (environ 30% du déplacement du piston ou des soupapes) on interpole la solution vers un nouveau maillage que l’on peut recommencer à déformer. Notons que les hexa-èdres présentent ici un gros avantage sur les tétraèdres car ils permettent des taux de déformation plus importants avant d’avoir recours à l’interpolation. Ainsi un cycle moteur complet se décom-pose en plusieurs phases de calcul représentées sur la figure 2.2(b). Un cycle moteur utilise 41 maillages au chapitre 4 et 66 au chapitre 5.

Conditions aux limites

Les conditions aux limites permettent de fermer un problème mathématique en spécifiant les valeurs que doivent prendre la solution ou ses dérivées aux frontières du domaine de calcul. La qualité de la solution dépend donc fortement du choix des conditions aux limites. Pour la simu-lation des équations de Navier Stokes, des conditions aux limites assurant d’avoir un problème mathématiquement bien posé ne donnent pas assez d’informations pour pouvoir les appliquer dans un code de calcul. En effet d’autres variables que celles déterminées par la condition aux limites sont nécessaires au schéma numérique. On les appelle conditions aux limites numériques. De plus la simulation d’écoulements compressibles en LES nécessite un contrôle précis des pro-priétés acoustiques des conditions aux limites D’un point de vue numérique, où l’espace et le temps sont discrétisés, l’avancée temporelle d’un nœud se fait en utilisant les valeurs aux nœuds adjacents.
Wjn+1 = Wjn − rj (Wk∈Dj )Δt (2.66)
L’équation 2.66 exprime l’avancée en temps des variables conservatives Wj au nœud j par le schéma numérique (cf. 2.2.2), Dj étant l’ensemble des nœuds voisins de j. Les nœuds du bord ne peuvent alors pas être mis à jour par le schéma numérique car les valeurs des nœuds hors du domaine manquent pour calculer le résidu rj2@ . Dans la suite nous appellerons fantômes les éléments de maillage hors du domaine de calcul.
Dans AVBP le schéma numérique calcule les résidus du bord dans un premier temps sans les valeurs des nœuds fantômes. On obtient alors des résidus dits prédits rjp2@ . Il faut ensuite corriger ces valeurs rjc2@ afin d’obtenir la condition aux limites désirée. Le code peut ensuite faire la mise à jour temporelle des variables aux nœuds du bord avec 2.66.
Il reste à déterminer la manière dont les résidus vont être corrigés. La technique classique consiste à imposer les grandeurs cibles de la condition aux limites et à laisser les valeurs prédites pour toutes les autres variables. Dans cette thèse deux méthodes sont utilisées pour effecteur cette opération : les conditions aux limites caractéristiques [42] dont le principe est abordé dans la section 2.3.3 et les lois de paroi [43–45].

Lois de paroi

L’une des difficultés des simulations LES concerne le traitement de l’écoulement proche pa-rois [46]. En effet, dans la couche limite turbulente, il existe des échelles turbulentes susceptibles de produire de l’énergie cinétique. Or dans ces régions, nous faisons le choix de ne pas raffi-ner le maillage pour ne pas pénaliser le temps de calcul. Ainsi le maillage est trop gros pour résoudre ces échelles et les modèles de sous maille utilisés ne prennent en compte que la dissi-pation d’énergie cinétique. Pour éviter de surestimer la dissipation aux parois, il est nécessaire de prendre en compte le mécanisme de production. Dans la formulation Cell-Vertex d’AVBP, les nœuds du bord coïncident avec la paroi. Pour l’utilisation de loi de paroi on fait l’hypothèse que ces nœuds sont en fait à une distance de la paroi suffisamment petite pour être négligeable devant la distance inter nœuds. Ainsi la vitesse au nœud n’est pas égale à celle de la paroi mais calculée par la loi de paroi. Dans la couche limite on fait les hypothèses suivantes :
les variables ne dépendent que de la direction normale à la paroi les propriétés thermodynamiques du fluide sont constantes
le gradient de pression est négligeable
le frottement et le flux de chaleur sont constants
La loi de paroi s’écrit (u+ = 1 ln(y+) + C pour y+ > 11

Modélisation LES de la combustion turbulente

Dans les moteurs à allumage commandé, la combustion est déclenchée par l’arc électrique formé aux électrode d’une bougie. La flamme se propage alors dans les gaz frais jusqu’à s’éteindre contre les parois du cylindre. On parle de combustion prémélangée dans les moteurs à injection indirecte et de combustion partiellement prémélangée ou de flamme de diffusion dans les mo-teurs à injection directe. Quel que soit le type de flamme rencontrée dans un moteur, la turbulence est primordiale car elle permet d’accélérer la phase de combustion et ainsi de consommer tout le carburant à temps. En effet si les soupapes d’échappement s’ouvrent avant la fin de la combus-tion, une partie du carburant est directement perdue et le cycle produit moins d’énergie. Dans le cadre de simulations LES de moteurs, une modélisation précise du phénomène de combustion est nécessaire à l’obtention de résultats prédictifs. Dans cette thèse seuls des cas de combustion pré-mélangée seront étudiés. Les combustions prémélangées turbulentes peuvent présenter plusieurs régimes en fonction du rapport des propriétés de la flamme aux propriétés de la turbulence. En moteur, l’épaisseur de flamme δl est d’environ 0,1 mm, ce qui est petit par rapport à l’échelle inté-grale de la turbulence lt (1 à 10 mm) [14]. La vitesse de flamme laminaire sl est également faible par rapport aux fluctuations de vitesse u , Vermorel et al. [7] estime 1 < u /sl < 10. Ainsi, en se basant sur le diagramme 2.3, l’hypothèse que la combustion se produit en régime de flammelette dans les moteurs est généralement admise.

Table des matières

Introduction générale
1 Introduction
1.1 Moteurs à allumage commandé
1.2 Variabilité de la combustion, transitoires et combustions anormales
1.3 La simulation aux grandes échelles dans les moteurs
1.4 Couplage de codes
1.5 Objectifs
2 Simulation d’écoulements réactifs : codes AVBP et Flow1D
2.1 Équations de la mécanique des fluides
2.1.1 Variables thermodynamiques
2.1.2 Équation d’état
2.1.3 Diffusion moléculaire multi-espèce
2.1.4 Diffusion de la chaleur
2.1.5 Coefficients de transport
2.1.6 Simulation aux grandes échelles (LES)
2.2 Description numérique du code AVBP
2.2.1 Discrétisation, méthode Cell-Vertex
2.2.2 Schéma Lax-Wendroff
2.2.3 Stabilisation : Viscosité artificielle
2.2.4 Gestion du maillage mobile : ALE et ITC
2.2.5 Conditions aux limites
2.2.6 Modélisation LES de la combustion turbulente
2.2.6.1 Le modèle ECFM-LES
2.2.6.2 Le modèle d’allumage ISSIM-LES
2.3 Modélisation 1D : développement du code Flow1D
2.3.1 Équations d’Euler quasi-1D
2.3.2 Schémas numériques
2.3.2.1 Lax-Wendroff two step (LW2S)
2.3.2.2 Lax-Wendroff cell vertex (LWCV)
2.3.3 Conditions aux limites
2.4 Validations du code Flow1D
2.4.1 Tube à choc
2.4.2 Tuyère de Laval
2.4.3 Écoulements pulsés
2.5 Modèles pour la simulation de système moteur 1D
2.5.1 Modèles de jonction
2.5.2 Modélisation des variations de section soudaines
2.5.3 Modèle de papillon d’admission
2.5.4 Arrête flamme
2.6 Conclusion
3 Développement d’une méthode de couplage 1D-3D LES
3.1 Introduction
3.2 La problématique du couplage
3.2.1 Continuité à l’interface
3.2.1.1 Chevauchement des domaines
3.2.1.2 Méthodes itératives
3.2.2 Méthodes des caractéristiques
3.2.3 Méthode de Godunov
3.2.4 Méthodes alternatives
3.3 Choix et implémentation de méthodes de couplage 1D-1D
3.4 Validations des couplages 1D-1D
3.4.1 Tuyère de Laval
3.4.2 Tube à choc
3.4.3 Propagation d’ondes acoustiques
3.4.4 Comparaison des méthodes
3.4.5 Sensibilité aux discontinuités dans les propriétés des gaz
3.4.6 Convergence en maillage
3.5 Implémentation du couplage 1D-3D dans AVBP
3.5.1 Réduction
3.5.2 Extrapolation
3.6 Architecture informatique
3.7 Validation du couplage 1D-3D
3.7.1 Impulsion acoustique
3.7.2 Onde acoustique
3.8 Conclusion
4 Étude du moteur SGEmac : comparaison entre calculs et expériences
4.1 Configuration expérimentale
4.2 Modélisation
4.3 Résultats
4.3.1 Moteur entrainé
4.3.1.1 Comparaison avec la PIV
4.3.2 Point de fonctionnement stable
4.3.3 Points de fonctionnement instables
4.3.3.1 Comparaison avec la LIF
4.3.4 Comparaison avec les calculs entièrement 3D
4.3.5 Temps de calcul
4.4 Conclusion
5 Application à l’étude de transitoires
5.1 Base de données expérimentales ASTRIDE
5.2 Régime stabilisé
5.2.1 Modélisation
5.2.2 Pressions dans l’admission, l’échappement et le cylindre
5.2.3 Aérodynamique interne
5.3 Transitoire de régime
5.4 Transitoire de charge
5.5 Conclusion
6 Conclusion
A Prediction of cyclic combustion variability in internal combustion engines via coupled 1D-3D LES method
B Development and validation of a two-way coupled 1D-3D LES method