Simulation Mathématique de la Circulation de l'Air dans une Enceinte avec Ouvertures et Sources Chaudes Localisées


MARC DUPUIS et J.-CLAUDE METHOT

Département de génie chimique, Université Laval, Ste-Foy, QC, Canada GIK 7P4

et

EDGAR DERNEDDE

Laboratoires Arvida et Centre de génie expérimental, Alcan International Limitée, Jonquière, Québec,
Canada G7S 4K8


La simulation mathématique de la circulation de l'air et du profil de température dans une enceinte rectangulaire avec plusieurs ouvertures et sources chaudes localisées a été réalisée sur ordinateurs. Une méthode numérique aux différences finies, combinée à un modèle simple de turbulence, a été utilisée pour résoudre les équations différentielles de la conservation de la masse, de la quantité de mouvement et de l'énergie. Les calculs ont été effectués à un nombre de Prandtl de 0,7, un nombre d'Archimède d'environ 34 et le nombre de Grashof a atteint des valeurs de l'ordre de 3 x 109 . Une comparaison entre les résultats numériques et ceux observés dans une maquette permet de conclure que le modèle mathématique reproduit assez fidèlement l'écoulement de l'air et le profil de température.

L'effet combiné de la convection naturelle au voisinage des sources chaudes et de la convection forcée engendre une circulation complexe dans l'enceinte. Ces conditions d'opération ont été choisies dans le but d'étudier les conditions de la circulation de l'air et de la température dans les salles de cuves d'électrolyse de l'aluminium.

A mathematical simulation of the air circulation and the temperature inside a rectangular enclosure with multiple openings and with localised rectangular heat sources was made. Finite difference techniques and a simple turbulence model were used to solve the time-averaged partial differential equations for the conservation of mass, momentum and energy. The computations were carried out for a Prandtl number of 0.7, for an Archimedes number of about 34 and for Grashof numbers up to 3 x 109. Comparison of the numerical results to flow visualization and temperature studies in an experimental apparatus indicated good agreement.

The combined effect of natural convection from the heat sources and the forced convection through the enclosure creates a complex flow pattern inside the enclosure. These conditions were simulated in an attempt to study the air circulation pattern and the air temperature inside the building of an aluminum reduction plant.



Introduction et revue de la literature

La qualité de l'air dans une salle de cuves d'électrolyse de l'aluminium dépend des émissions de gaz provenant des cuves et de la circulation de l'air dans la salle. Cette circulation est causée par la ventilation, souvent forcée, de la bâtisse et par la convection naturelle autour des cuves, lesquelles dégagent une grande quantité de chaleur.

Pour améliorer la qualité de l'air dans une salle de cuves, on peut étudier la circulation en usine, sur une maquette ou par simulation mathématique. Les études en usine sont dispendieuses et difficiles. Les études sur maquette prennent beaucoup de temps en laboratoire. Heureusement, le récent progrès dans les méthodes numériques et les ordinateurs ouvre une vole intéressante sur l'utilisation de la simulation mathématique. La simulation mathématique de la convection naturelle et de la convection forcée a déjà fait l'objet d'un grand nombre de publications.

L'équipe Churchill (Chu et Churchill, 1977 et Wilkes et Churchill, 1966) fut probablement la première à adopter la méthode numérique ADI (Alternating Direction Implicit) aux ordinateurs modernes pour simuler la convection naturelle de l'air dans une enceinte rectangulaire et fermée. Bien que révolutionnaire à l'époque, la méthode ADI, combinée à la technique CDS (Central Difference Scheme), devient passablement instable lorsque le nombre de Grashof excède la valeur de 1,5 X 105 environ.

Toutefois, ces premiers travaux dans ce domaine ont été à l'origine d'un grand nombre de travaux similaires, rendus possibles grâce à l'amélioration des nouvelles méthodes numériques et au rendement accru des ordinateurs de plus en plus sophistiqués et puissants. Szekely et Todd (1971) ont amélioré la stabilité de la méthode ADI en introduisant des itérations sur la vorticité aux parois à chaque intervalle de temps. Malgré cet ajout, les calculs ont été instables pour un grand nombre de Grashof supérieur à 6,0 X 105. Pepper et Harris (1977) ont utilisé la méthode SIP (Strongly Implicit Procedure) avec itérations dans une enceinte fermée et ils ont réussi à obtenir une solution stable pour un nombre de Grashof égal à 1,0 x 106.


D'autres auteurs ont employé la méthode numérique explicite, laquelle se prête bien à l'adaptation de la technique UDS (Upwind Difference Scheme). Hellums et Churchill (1961), Barakat et Clark (1966), MacGregor et Emery (1969) ont réussi à maintenir la stabilité numérique même à un nombre de Grashof de 1,0 x 107. À cette valeur du nombre de Grashof, I'écoulement près des parois chaudes devient turbulent. Torrance et Rockett (1969) ont étudié le cas de l'enceinte fermée avec une source de chaleur locale sur le plancher pour un système d'axes cylindriques. Ces derniers ont utilisé la méthode explicite avec la technique UDS du deuxième ordre et ont obtenu des résultats stables jusqu'à un nombre de Grashof de 4 x 1010.

Gosman et al. (1969) ont aussi adopté la technique UDS pour aborder le cas épineux de la convection forcée en écoulement turbulent, tout en introduisant un nouveau modèle de turbulence à deux équations différentielles. Cette équipe a étudié la circulation forcée dans une enceinte ouverte avec source de chaleur (Nielsen et al., 1979) étendue à tout le plancher.

Dans la même ligne de pensée que celle des travaux ci-haut mentionnés, nous présentons un résumé des travaux de recherche qui ont conduit à la simulation mathématique de la circulation d'un fluide visqueux et turbulent dans une enceinte rectangulaire et ouverte, pour un nombre d'Archimède égal à 34 et pour une géométrie relativement complexe par rapport aux travaux antérieurs.

Théorie des écoulements visqueux

La circulation de l'air dans une salle de cuves s'effectue essentiellement en deux dimensions où la pression statique a peu d'intérêt. Par conséquent, la formulation des équations bidimensionnelles en fonction de la vorticité et de la fonction de courant représente un avantage marqué.

Les équations différentielles adimensionnelles représentant un écoulement visqueux bidimensionnel sont les suivantes (Roache, 1982; White, 1970):

1. Équation de la vorticité:

Où ve et vo représentent les viscosités effective et laminaire respectivement.

2. Équation de l'énergie T:

3. Équation de la fonction de courant:

Les vitesses horizontales et verticales de l'écoulement sont obtenues à partir de la définition de la fonction de courant comme suit:

Ces équations sont applicables aux écoulement laminaires (ve = vo) ainsi qu'aux écoulements turbulents, à la condition d'utiliser les valeurs moyennes des variables vortivité, température et fonction de courant.

L'écoulement de l'air dans les salles de cuves étant de nature turbulente, un terme ve a été introduit dans les équations de la vorticité et de l'énergie. Ce nouveau terme représente la viscosité cinématique effective, c'est-à-dire, la viscosité cinématique laminaire à laquelle vient s'ajouter un terme qui tient compte de l'augmentation de la viscosité du fluide, causée par la turbulence. L'hypothèse de Prandtl pour les écoulements libres et turbulents a été appliquée pour évaluer ce terme qui peut être représenté dans un système à deux dimensions par l'équation suivante (Gosman et al., 1969):

Doù d représente une longueur de mélange quasi constante à déterminer en fonction de la profondeur des jets turbulents. Comme la turbulence est plus grande dans les jets d'air où la vitesse est relativement élevée, on a retenu après de nombreux essais les définitions suivantes pour dv et dT :

Un terme a aussi été ajouté pour tenir compte de la turbulence causée par la convection naturelle.

Les méthodes numériques

Les trois équations différentielles du système peuvent être résolues par la méthode des différences finies et la présente méthode numérique s'inspire principalement des travaux de Wilkes et Churchill (1966). En effet, la méthode ADI a été utilisée pour résoudre les deux équations paraboliques (vorticité et énergie), tandis que la méthode SOR (Successive Overrelaxation) se prêtait mieux à la résolution de l'équation elliptique (fonction de courant).

À partir de la fonction de courant, les vitesses horizontales et verticales ont été déterminées par des différences centrées partout dans le maillage.

Les deux méthodes choisies, ADI et SOR, sont bien connues pour leur stabilité et leur convergence. Mais, la stabilité de la méthode ADI, théoriquement assurée, est limitée en pratique par les conditions aux parois. Cette restriction est similaire à celle de la méthode explicite (Roache, 1982).


Afin d'assurer la stabilité numérique des deux équations paraboliques à des nombres de Grashof élevés, les termes convectifs ont été discrétisés par la technique UDS. Cette technique, décrite entre autres par Torrance et Rockett (1969), consiste à évaluer les termes convectifs en utilisant les informations en amont du courant. De cette façon on obtient les équations suivantes:

La même discrétisation s'applique dans l'autre direction. Cependant, on doit signaler que l'emploi de cette technique introduit un terme de diffusion numérique qui pourrait fausser la solution (Patankar, 1980; Raithby, 1976; Roache, 1982).

Grâce au modèle de turbulence, la valeur de la viscosité est augmentée dans plusieurs zones locales. Cette augmentation réduit l'effet de la diffusion numérique.

Les propriétés de transfert de l'air, lesquelles n'étaient pas constantes à cause de la turbulence, étaient évaluées par les moyennes arithmétiques (Patankar, 1980).

L'enceinte et les conditions aux limites

L'enceinte comprend quatre ouvertures, trois sur le plancher pour les entrées d'air (convection forcée) et un canal au plafond, constituent la seule sortie. Deux rectangles solides, représentant les cuves d'électrolyse et placés séparément au plancher, constituent les sources de chaleur (convection naturelle) en plus de faire partiellement obstacle à l'écoulement de l'air (figure 1). L'enceinte est symétrique par rapport à son axe vertical central.


Figure 1 - Géométrie de l'enceinte et évaluation de la vorticité autour d'un coin.

a) LA VORTICITÉ

Dans les problèmes réels d'écoulement de fluide, les conditions aux limites pour la vorticité ne sont pas données explicitement.

Selon la précision désirée, on peut obtenir plusieurs équations pour la vorticité à la surface (Roache, 1982). Dans le cas présent, l'équation suivante a été utilisée.

L'évaluation de la vorticité aux limites sur un coin est particulièrement difficile à cause de la discontinuité de la surface. La salle de cuves contient plusieurs coins du même genre: les coins des cuves, les coins aux entrées et à la sortie de l'air. Plusieurs solutions à ce problème particulier ont été proposées (Roache, 1982). À l'usage, la formulation qui consiste à utiliser la valeur moyenne de la vorticité de part et d'autre du coin (Figure 1) a été retenue:

Dans les entrées, loin des surfaces, on a supposé que la vorticité était nulle. À la sortie du canal, loin des surfaces, on a posé que le gradient de vorticité perpendiculaire à la sortie était nul. La vorticité à la sortie est un des résultats de la simulation.

Pour obtenir une solution numérique stable, il y a une précaution à prendre au sujet de ces équations. Comme le calcul des valeurs de la vorticité aux parois est retardé d'un intervalle de temps par rapport aux valeurs intérieures, il faut minimiser l'effet de ce décalage dans le temps par l'itération des équations de la fonction de courant et de la vorticité à chaque intervalle de temps (Szekely et Todd, 1971).

b) LA FONCTION DE COURANT

Ici, trois genres de conditions aux limites prévalent pour l'enceinte étudiée. Aux parois, la valeur de la fonction de courant demeure constante. Dans les ouvertures (des entrées et sorties), la valeur varie linéairement en fonction du débit:

où u est la vitesse d'entrée ou de sortie.

À la sortie du canal, le gradient de fonction de courant perpendiculaire à la sortie a été fixé comme nul.

Par la suite, les valeurs de la fonction de courant peuvent être définies en tout point aux parois. À cette fin, il suffit de choisir un point de référence où la fonction de courant = 0. Pour le présent travail, ce point de référence a été fixé à l'axe central vertical. II s'ensuit qu'il y a une ligne de courant de valeur nulle sur cet axe puisqu'aucun débit d'air ne la traverse, dans le cas où le vent n'a aucun effet sur le système (débits d'air symétriques dans les entrées).

c) LA TEMPÉRATURE

Les conditions aux limites pour la température sont fixées en fonction de la position sur la paroi. Les murs, le plancher et le plafond de l'enceinte sont considérés comme parfaitement isolés, ce qui se traduit par les équations suivantes:

Les entrées d'air ont une température égale à la température de référence adimensionnelle:


La température adimensionnelle des surfaces chaudes a été fixée comme suit, indépendamment du temps:

Par contre, pour limiter le taux de transfert de chaleur, la chaleur requise était générée dans les éléments différentiels près des surfaces chaudes et ces surfaces ont été considérées comme parfaitement isolées. Cette technique a permis d'obvier à la difficulté de simuler le transfert de chaleur sur les surfaces chaudes par convection lorsque le maillage est grossier.

La température de l'air à la sortie de l'enceinte constitue un des résultats de la simulation; on a fait l'hypothèse que le gradient de température de l'air à la sortie de l'enceinte était nul en direction perpendiculaire.

La maquette

Afin de vérifier expérimentalement les résultats du modèle mathématique, des essais ont été réalisés à l'échelle réduite à l'aide d'une maquette. Les figures 2 et 7 donnent un aperçu du montage expérimental dont les principales dimensions sont les suivantes: longueur 0,914 m, hauteur 0,508 m, profondeur 0,229 m. Les deux sources de chaleur au plancher sont des éléments chauffants dans des blocs d'aluminium de 0,051 m de hauteur, 0,089 m de longueur et de pleine profondeur. La température de chacune de ces sources de chaleur est de 100°C et la température de l'air aux entrées est égale à 26°C. La longueur des deux entrées latérales est égale à 0,026 m et elles sont situées à une distance de 0,102 m des murs. L'entrée centrale a une longueur égale à 0,051 m.


Figure 2 - Simulation physique de la circulation de l'air: cuves chaudes.

La vitesse de l'air dans les trois entrées au plancher est évaluée à 0,25 m/s. L'air chaud est aspiré à la sortie centrale au plafond par un ventilateur. Le débit a été évalué à 0,54 m3/min. La différence de température entre l'entrée et la sortie d'air est évaluée à 8°C. L'énergie fournie par les éléments chauffants et transmise à l'air circulant dans la maquette est égale à 85 W.

Les conditions d'opération dans la maquette sont telles que le nombre d'Archimède est égal à 34.

Afin de mettre en évidence les effets séparés de la convection naturelle causée par les sources chaudes et de la convection forcée occasionnée par le ventilateur, deux séries d'essais ont été réalisées: avec, puis sans chaleur fournie aux cuves. Pour visualiser l'écoulement de l'air dans la maquette, de la fumée de cigarettes a été injectée par de petite trous à la surface des cuves chaudes. Des lignes de courant ont été tracées approximativement à la main après étude détaillée de plusieurs films et photographies de la maquette en opération. À 70 endroits pré-sélectionnés, la vitesse de l 'air a été évaluée directement à l'intérieur de la maquette avec un anémomètre (TSI, modèle 1650). La température de l'air dans la maquette a été évaluée à 332 endroits différents avec un petit thermocouple.

L'écoulement de l'air froid dans les entrées et de l'air chaud au?dessus des cuves chaudes est turbulent comme le démontre la circulation de la fumée. Le nombre de Reynolds de ces jets d'air est approximativement égal à 700. Cette valeur, supérieure à 30, confirme que 1'écoulement de 1'air dans ces jets est en écoulement turbulent (Schlichting, 1968). Près du plafond, des murs et du plancher, la vitesse de l'air est faible et le nombre de Reynolds est approximativement égal à 1600. Cette valeur, inférieure à 100,000, indique que l'écoulement de l'air près de ces parois reste laminaire (Schlichting, 1968).

Les résultats

Les résultats de la simulation mathématique ont été obtenus avec un maillage grossier (dx = 0,0214, dy = 0,0179) de 1653 noeuds (29 lignes et 57 colonnes). Ce maillage a été choisi comme compromis en tenant compte du coût sur ordinateurs et de la précision numérique. La valeur de dt = 0,5 x 106 utilisée était très petite. L'expérience a démontré que l'emploi d'une valeur de dt plus grande était plus coûteux, puisqu'un grand nombre d'itérations sur la vorticité aux limites étaient nécessaires avec les méthodes ADI et SOR. Même avec un coefficient de surrelaxation optimal F = 1,75, la méthode SOR a nécessité 75% du temps de calcul à chaque pas de temps. À chaque pas de temps, une vingtaine d'itérations sont nécessaires pour résoudre l'équation de la fonction de courant par la méthode SOR. Pour les dix premiers pas de temps, deux ou trois itérations sont nécessaires pour évaluer correctement les conditions de vorticité aux limites. Par la suite, étant donné la petitesse du pas de temps, une seule itération est requise pour atteindre le même but.

La convergence vers l'équilibre à partir de zéro nécessitait 450 pas dans le temps, représentant près de 45 minutes de temps machine sur un ordinateur IBM 4341.

Les résultats des simulations mathématiques sont comparés à ceux de la maquette.

a) AVEC CHALEUR FOURNIE AUX CUVES

La circulation de l'air dans une salle de cuves d'électrolyse a été étudiée dans le cas où les cuves étaient chaudes, à l'aide du modèle mathématique et de la maquette. Les lignes de courant expérimentales, c'est-à-dire à partir des essais à l'échelle réduite, ont été obtenues de deux façons différentes: en observant la circulation de la fumée de cigarettes à l'intérieur de la maquette, comme le montre la figure 2, et à partir des vitesses évaluées à 1'intérieur de la maquette, comme le montre la figure 3. Ces deux méthodes donnent approximativement les mêmes résultats.


Lorsque les cuves sont chauffées, les conditions de circulation de l'air dans la maquette sont les suivantes:

- le jet central se divise en deux parties qui se dirigent vers les cuves;
- les jets latéraux sont aspirés vers les cuves;
- au-dessus des deux cuves, un fort courant d'air chaud monte vers le plafond à 0,13 m/s, se sépare en deux parties, l'une prenant la direction de la sortie tandis que l'autre redescend au plancher près des murs à une vitesse de 0,03 m/s;
- les courants d'air chaud génèrent quatre grandes boucles de recirculation par l'entraînement d'air le long des cuves.


Figure 3 - Simulation physique des lignes de courant: cuves chaudes.


Figure 4 - Simulation mathématique des lignes de courant: cuves chaudes.

La figure 4 montre les résultats de la simulation mathématique pour les mêmes conditions d'opération. Une comparaison entre les figures 3 et 4 permet de conclure que le modèle mathématique reproduit bien l'écoulement de l'air dans l'enceinte. La vitesse maximale de l'air au-dessus des cuves a été calculée à l'aide du modèle mathématique à 0,136 m/s en accord avec les résultats de la maquette.

Le profil de température dans la maquette présentée à la figure 5 a les caractéristiques suivantes:

- la température de l'air augmente rapidement près du plancher et très lentement ailleurs en fonction de la hauteur;
- les courants d'air chaud au-dessus des cuves ont une température de 1°C plus élevée que l'air adjacent à ces courants;
- même près des cuves, ces courants d'air chaud ont une température de 8°C seulement plus élevée que l'air adjacent à ces courants.

La figure 6 présente les isothermes de 1'air dans 1'enceinte calculés par le modèle mathématique. Ces isothermes donnent les caractéristiques similaires à ceux de la maquette.


Figure 5 - Simulation physique des lignes de courant: cuves chaudes.


Figure 6 - Simulation mathématique des lignes de courant: cuves chaudes.

b) SANS CHALEUR FOURNIE AUX CUVES

Le cas des cuves froides représente peu d'intérêt dans cette étude, mais il démontre la flexibilité du modèle mathématique. La comparaison est donc plutôt qualitative et les lignes de courant ont été tracées selon les circulations de la fumée émise par les cuves.

Comme le montre la figure 7, lorsqu'aucune chaleur n'est fournie aux cuves, les conditions de circulation d'air dans la maquette sont les suivantes:


- le jet d'air central se dirige vers la sortie en suivant une trajectoire rectiligne, parallèle à l'axe de symétrie;
- les deux jets d'air latéraux s'inclinent graduellement vers les murs de la maquette, longent ceux-ci et le plafond avant de s'engager vers la sortie;
- les jets d'air ont un effet d'entraînement très important; sur le plancher, ils aspirent le traceur et engendrent six boucles de recirculation dont les quatre principales se situent au-dessus des cuves.

La figure 8 présente les lignes de courant calculées par le modèle mathématique pour les mêmes conditions d'opération. Dans l'ensemble, le modèle mathématique reproduit bien la complexité de l'écoulement identifié sur la maquette, à deux exceptions près.

Les jets d'air latéraux se dirigent plus rapidement vers la sortie sans s'incliner vers les murs du système et les boucles de recirculation près des murs sont plus grandes que celles de la maquette.

Les boucles de recirculation près du jet central sont plus petites.


Figure 7 - Simulation physique de la circulation de l'air: cuves froides.


Figure 8 - Simulation mathématique des lignes de courant: cuves froides.

Conclusion

La méthode numérique aux différences finies, basée sur les travaux de Wilkes et Churchill (1966), ainsi que sur ceux de Szekely et Todd (1971), Torrance et Rockett (1969), et Gosman et al. (1969) a permis d'obtenir une solution stable pour un problème relativement complexe.

De par ses conditions aux limites complexes et de par ses écarts de température considérables, le cas traité ici est relativement plus compliqué que ceux présentés dans la littérature jusqu'à ce jour.

Dans le cas des cuves chaudes, lorsque les paramètres du modèle de turbulence sont bien ajustés, on arrive à reproduire avec une précision inespérée les profils d'écoulement et de température observés expérimentalement.

Bien que les comparaisons entre l'écoulement d'air dans la maquette et celui prédit par le modèle mathématique soient qualitatives seulement pour les cuves sans chaleur, on peut constater que le modèle mathématique reproduit assez fidèlement le changement du profil d'écoulement de l'air dans l'enceinte lorsque la convection forcée est la seule cause du mouvement de l'air.

La diffusion numérique est approximativement égale à la valeur de ve/vo. Mais quelques essais avec des valeurs de la viscosité turbulente augmentées ou diminuées par un facteur 4 ont très peu changé les résultats de la simulation. En effet, dans l'ensemble de la salle sauf très près des cuves, la convection domine largement la diffusion. On estime donc que la diffusion numérique ne peut fausser de manière significative les résultats.

Suite à ces explorations au sujet du modèle mathématique, on peut confirmer deux remarques déjà soulevées par d'autres auteurs. La méthode ADI a compliqué de beaucoup la programmation; à cause de la présence des cuves et du canal de sortie, plusieurs points se trouvent même en dehors du domaine de calcul. La méthode SOR s'est avérée coûteuse, car elle utilise approximativement 75% du temps de calcul.

Le choix du type de formulation pour évaluer la vorticité aux parois, particulièrement dans les coins, s'est avéré un facteur important dans la simulation mathématique. Plusieurs techniques UDS et des modèles de turbulence plus sophistiqués sont disponibles dans la littérature. Les travaux en cours sur ces trois sujets feront l'objet d'autres publications prochainement.

Remerciements

Les auteurs tiennent à remercier la compagnie Alcan International Limitée pour son appui financier et technique au projet de recherche ainsi que le Conseil de Recherches en Sciences Naturelles et en Génie pour les mêmes raisons et pour une bourse d'études accordée à monsieur Dupuis.

Nomenclature


Symboles grecs

Indices

Exposants

Bibliographie

  • Barakat, H. Z. et J. A. Clark, "On the Solution of the Diffusion Equations by Numerical Methods", J. Heat Trans. 88, 421 - 427 (1966).
  • Chu, H. N. S. et S. W. Churchill, "The Development and Testing a Numerical Method for Computation of Laminar Natural Convection in Enclosures", Computers Chem. Eng. 1, 103 - 108 (1977).
  • Gosman, A. D., W. M. Pun, A. K. Runchal, D. B. Spalding et M. Wolfshtein, "Heat and Mass Transfer in Recirculating Flows", Academic Press, London (1969).
  • Hellums, J. D. et S. W. Churchill, "Computation of Natural Convection by Finite Difference Methods", Int. Developments in Heat Transfer, ASME, New York, 5, 985-994 (1961).
  • MacGregor, R. K. et A. F. Emery, "Free Convection Through Vertical Plane Layers-Moderate and High Prandtl Number Fluids", J. Heat Trans. 91, 391-403 (1969).
  • Nielsen, P. V., A. Restivo et J. H. Whitelaw, "Buoyancy Affected Flows in Ventilated Rooms", Numerical Heat Transfer, 2, 115 - 127 (1979).
  • Patankar, S. V., "Numerical Heat Transfer and Fluid Flow", Hemisphere Publishing Corporation, New York (1980).
  • Pepper, D. W. et S. D. Harris, "Numerical Simulation of Natural Convection in Closed Containers by Fully Implicit Method", J. Fluids Eng. 99, 649-656 (1977).
  • Raithby, G. D., "A Critical Evaluation of Upstream Differencing Applied to Problems Involving Fluid Flow", Computer Methods Appl. Mechanics Eng. 9, 75-103 (1976).
  • Roache, P. J., "Computational Fluid Dynamics", Hermosa Publishers, Albuquerque (1982).
  • Schlichting, H., "Boundary Layer Theory", McGraw-Hill, New York (1968).
  • Szekely, J. et M. R. Todd, "Natural Convection in a Rectangular Cavity Transient Behaviour and Two-Phase Systems in Laminar Flow", Int. J. Heat Mass Trans. 14, 467?482 (1971).
  • Torrance, K. E. et J. A. Rockett, "Numerical Study of Natural Convection in an Enclosure with Localized Heating from below -Creeping Flow to the Onset of Laminar Instability", J. Fluid Mech. 36, 33 - 54 (1969).
  • White, F. M., "Viscous Fluid Flow", McGraw-Hill, New York (1970).
  • Wilkes, J. O. et S. W. Churchill, "The Finite Difference Computation of Natural Convection in a Rectangular Enclosure", AlChE J. 12, 161 - 166 (1966).