Les équations régissant le phénomène étudié sont des équations aux dérivées partielles (edp) non-linéaires, dont la résolution analytique ne peut être possible
télécharger 87.86 Kb.
|
| CHAPITRE III MÉTHODE NUMÉRIQUE III.1.Introduction Les équations régissant le phénomène étudié sont des équations aux dérivées partielles (EDP) non-linéaires, dont la résolution analytique ne peut être possible au moyen des outils d’analyse mathématique contemporain. Mais, une solution numérique peut être possible en transformant ces équations différentielles en systèmes d'équations algébriques linéaires par une méthode de discrétisation avant de résoudre ce système par des méthodes directes ou par itérations. Pour notre présente étude, nous avons choisi la méthode des volumes finis pour discrétiser les équations du modèle mathématique. Pour obtenir une solution numérique, le problème étudié doit être discrétisé en transformant les équations différentielles en un système d'équations algébriques. Il existe plusieurs méthodes de discrétisation telles que: la méthode des volumes finis, des différences finies et des éléments finis, ….etc. Parmi ces méthodes, nous avons choisi la méthode des volumes finis. III.2. Présentation de la Méthode Le principe de la méthode des volumes finis consiste à intégrer les équations de transport sur un ensemble discret de volume finis jointifs couveront le domaine physique. Le résultat de la discrétisation en un point est une équation algébrique liant la valeur d’une variable aux valeurs des variables des points voisins et d’autres considérés constants. Sa simplicité, sa facilité de linéarité les termes sources et son garantit de conserver la quantité de masse, de mouvement et d’énergie dans tout le domaine étudié donne un choix judicieux d’adopter cette méthode.  Fig(III-1) :Schémas représentant un volume de contrôle III.3 Maillage Le domaine physique est discrétisé en domaine de calcul suivant un maillage uniforme ou non uniforme dans les deux directions, horizontales verticales (figureIII-2). Une suite géométrique est utilisée afin de raffiner d'avantage le maillage au niveau où on a des variables dépendantes (u, v, p, θ, ρ). Dans ces régions les frontières du domaine coïncident avec les faces des volumes de contrôle, ce qui facilite l'incorporation des conditions aux limites. Chaque noeud du maillage est repéré par deux indices I et J, donnant sa position suivant les directions x et y. Le nombre total des noeuds est NI et NJ, tel que : NI nombre total suivant x. NJ nombre total suivant y. III.3.1.Stockage des Variables Les variables dépendantes scalaires (Ρ, θ, ρ) sont stockées aux nœuds du maillage, tan disque les variables dépendantes vectorielles (u, v) sont stockées au milieu des segments reliant les nœuds. Si on note par P le centre du volume de contrôle de la variable Φ et Ε, W, Net S les nœuds voisins des volumes de contrôle adjacents comme illustre les figures (III-2) et (III-3). Ces nœuds seront les lieux de stockage des variables scalaires (Ρ,Τ, ρ).Au milieu de chaque segment reliant deux nœuds adjacents on note par е, w, n et s où sont stockées les variables vectorielles.  Fig(III-3) : Localisation des variables III.3.2.Maillage Décalé (Staggered grid) La discrétisation d'une équation de transport de diffusion sur un volume de contrôle par la méthode des volumes finis fait intervenir les valeurs des vitesses aux interfaces des volumes (ue, uw, vn, vw) .Il est donc intéressant de les calculer directement sur les interfaces sans avoir à effectuer d'interpolation. D'autre part, la discrétisation de l'équation de continuité et du gradient de pression avec l'utilisation d'une interpolation linéaire peut induire des erreurs importantes du fait qu'une répartition de pression ou de vitesse en "damier" est vue comme un champ uniforme. Pour contourner ces difficultés on préfère utiliser des grilles décalées "staggered grid". Une grille principale est construite sur laquelle on calcule la pression et la température. Deux grilles décalées vers la droite et vers le haut respectivement sont utilisées pour le calcul des vitesses horizontales et verticales (Fig. (III-4)).  Fig. (III-4) : Maillage décalé III.4 Equation générale de transport L’équation générale de transport d’une variable φ pour un écoulement axisymétrique et incompressible s’écrit comme suit :  1 : Représente le terme transitoire. 2 : Représente le terme de transport par convection. 3 : Représente le terme de transport par diffusion. 4 : Représente le terme de source. Avec : U : Composante radiale de la vitesse. V : Composante axiale de la vitesse. Γ : Coefficient de diffusion. Dans le tableau suivant, nous donnons la définition de φ , Γ et φ S pour les équations qui gouvernent notre problème générale. Tableaux (III.1) : Les variables et les coefficients des équations de transport adimensionnelles.
III.4 Discrétisation des équations mathématiques III.4.1 Intégration l’équation générale de transport Intégrons l'équation (III-4) à travers le volume de contrôle décrit par la figure (III-1).  La division par  , nous donne :Le terme transitoire :  Le terme convectif :   Le terme diffusif :   Le terme source :  L’équation (III.3) s’écrira alors :  + + =  + + En arrivant à ce stade, il faudra exprimer les termes des flux convectifs et diffusifs aux interfaces des volumes de contrôle. Afin de surmonter à ce problème, on fait appel aux schémas de discrétisation (différences centrées, exponentiel, Power-Law, hybride, …). Ces schémas différents par la façon avec laquelle on prend en compte les termes de convection et de diffusion. La présentation de la forme générale de l’équation algébrique discrétisée s’écrit comme suit :   Avec :        Où      Et       , , ,  Sont les coefficients correspondants, respectivement, aux nœuds Est,Ouest, Nord, Sud et centre du volume. b : est un terme de source supposé être constant dans le volume de contrôle.  ont respectivement les termes convectifs et diffusifs aux faces Est, Ouest, Nord, Sud. désignent le rapport du flux convectif au flux diffusif (nombres de Peclet) aux différents faces du volume de contrôleIII.4.2 Différentes schéma de discrétisation Les profils approximatifs décrivant la variation de  entres les nœuds, sont exprimés par la fonction  spécifique pour chaque schéma numérique. S. Patankar [38], à cité les cinq schémas suivants:Tableau 5.2
III.4.3 Discrétisation de l’équation de quantité de mouvement suivant X L’intégration de l’équation de quantité de mouvement suivant X sur un volume de contrôle décalé vers la droite (voir Figure 3.3.a) donne l’équation algébrique :  (III.24)Avec:   (III.25)    Discrétisation du terme source :  Après l’intégration, le terme source devient :  Les termes convectifs :    (III.28)  Les termes diffusifs :    (III.29)  III.4.4 Discrétisation de l’équation de quantité de mouvement suivant Y L’intégration de l’équation de quantité de mouvement suivant Y sur un volume de contrôle décalé vers le haut (voir Figure 3.3.b) donne l’équation algébrique :  Avec: (III.25)Discrétisation du terme source :    Après l’intégration, le terme source devient :    Les termes convectifs :    (III.28)  Les termes diffusifs :    (III.29)   Figure 3.3 : (a) Volume de contrôle décalé vers la droite. (b)Volume de contrôle décalé vers le haut III.4.5 Discrétisation de l’équation de l’énergie L’intégration de l’équation de l’énergie adimensionnelle (2.14) sur un volume de contrôle typique (voir figure 3.1) donne l’équation algébrique :  Avec: (III.25)Discrétisation du terme source :  Les termes convectifs :    (III.28)  Les termes diffusifs :    (III.29)  III.4.6 Discrétisation de l’équation de concentration L’intégration de l’équation de l’énergie adimensionnelle (2.15) sur un volume de contrôle typique (voir figure 3.1) donne l’équation algébrique :  Avec: (III.25)Discrétisation du terme source :  Les termes convectifs : (III.28) Les termes diffusifs : (III.29) III.5 Algorithme SIMPLER On a vu que la discrétisation a remplacé les équations aux dérivées partielles par des systèmes d’équations algébriques. La résolution de ce système présente quelques difficultés, parce que: • Les coefficients apparaissant dans l’équation de discrétisation dépendent des variables donc l’équation n’est pas linéaire. • Les termes de sources des équations de quantités de mouvement, impliquent les gradients de pression or nous ne disposons pas d’équation pour cette variable. Grâce à l’algorithme SIMPLER (Semi Implicit Method for Pressure-Linked Equation Revised) (Patankar[48] (1980), parce que c’est une méthode itérative qui permet justement le calcul des vitesses et la pression. Après convergence de la solution, les champs de vitesse et de pression doivent satisfaire simultanément l’équation de continuité et l’équation de quantité de mouvement. Après intégration sur un volume fini et discrétisation des différents termes on aboutit à la forme algébrique des équations de quantité de mouvement    S ur la base d’un champ de pression estimée  , ces équations donnent un champ de vitesses  et  qui ne satisfont pas l’équation de continuité. On a les relations :  Pour aboutir à des champs corrects de vitesse et de pression  il faut corriger Comme suit : Comment déterminer les corrections de vitesse  et la pression ?Soustrayons (III-14) de (III-13) on obtient :   A ce stade les termes  et  sont omis. A noté que la solution finale des champs de vitesse et de pression ne contiendra pas d’erreur due à cette omission puisque tous les thermes de ces équations tendent vers zéro. On obtient donc, pour les 4 faces e, w, n, s du volume central.     : Surface   : Surface  En fait les vitesses peuvent être corrigées à condition que l’on ait une estimation du champ des corrections des pressions. C’est l’équation de continuité qui va être transformée pour donner une équation des corrections de pression. L’intégration et la discrétisation de l’équation de continuité par rapport à un volume de contrôle central (pour un fluide compressible) donne.  Substituons maintenant à la place des vitesses les expressions reliant les corrections de pression   Ou   La séquence des opérations essentielles constituant l’algorithme SIMPLER [62] est la suivante : 1- Estimer un champ de vitesse. 2-Calculer les coefficients de l’équation de quantité de mouvement et par conséquent calculer les pseudo-vitesses  donnée par les équations :(III-22a)-(III-22d)3-calculer les coefficients de l'équation de pression (III-24), et on la résoudre pour obtenir le champ de pression 4-Considérer le champ de pression comme estimation ., et résoudre l’équation de mouvement pour obtenir  .5-Calculer le terme source bm (Eq. III-20) et par conséquent résoudre l’équation de correction de pression  .6-Corriger le champ de vitesse en utilisant l’équation (III-16), mais ne pas corriger la pression. 7-Résoudre : - l’équation de l’énergie, et obtenir  .- l’équation de concentration, pour obtenir  .- l’équation de potentiel électrique  8- Retourner à l’étape 2, avec le nouveau champ de vitesse jusqu’à convergence. III.6 Méthode itérative de résolution (Algorithme TDMA) La résolution directe du système d’équations algébriques est compliquée, pour cela, on utilise la technique de balayage qui est une méthode de résolution semi-itérative. Elle consiste à déterminer les valeurs de la variable φ sur chaque ligne du domaine d’étude indépendamment des autres lignes, donc le système se transforme d’un système d’équations algébriques multidimensionnelles en un système unidimensionnel, en ajoutant à la source de la dimension choisie des termes des autres dimensions. Le système d’équations obtenu est représenté par une matrice tridiagonale et peut être résolu par l’algorithme de Thomas. L’équation algébrique s‘écrit pour le noeud P du maillage comme suit : Le système d’équations obtenu peut se mettre sous la forme : (III.46) : est une matrice de (NxN) élements. : Vecteur des inconnues  avec  and  Pour déterminer les valeurs de φ sur une colonne « i » on suppose que les valeurs de cette dernière sont connues sur les colonnes « i −1 » et « i + 1 ». L’équation algébrique (III-45) écrire pour chaque noeud de la colonne « i » est alors réduire à une équation qui contient seulement trois inconnues (    ).Pour le nœud (i, j) du maillage, l’équation peut être écrite sous la forme d’une équation unidimensionnelle :  Et en posant :     On obtient l’équation (III.48) sous la forme suivante :  avec : Pour tous les nœuds [j = 2, JL −1]de la colonne i , l’équation (III.48) donne un système de la forme :   Les valeurs de 1 φ et JL φ sont connues (conditions aux limites). La matrice associée au système est tridiagonale. On utilisera l’algorithme TDMA [48] (algorithme de THOMAS) pour la résolution en réarrangeant toutes les équations du système (III.49) sous la forme :  On obtient ce qui suit :  Détermination of  Pour le nœud (i, j −1) , on a :  (III.51)En remplaçant (III.55) dans (III.49) on trouve :  D’où on a :  . Remarquons que pour : j = 1 on a : = 0 j c ; l’équation (III.58) pour j = 1 se réduit à :  Ce qui correspond à la forme de l’équation (III.50). Donc, on a :   Aussi, pour  we have so and L’algorithme de THOMAS se résume comme suit :
1-Calculer  et  de (III.62). 2-Calculer à partir de (III.59) et (III.60) les coefficients j P et J Q avec : j = 2, JL .
3. On pose  4. On utilise la relation (III.54) pour j = JL −1 , JL − 2 , ........1 afin d’obtenir les valeurs  |
