Convection magnétohydrodynamique dans un fluide

non-newtonien saturant un milieu poreux

Cyrine Chahtour

To cite this version:

Cyrine Chahtour. Convection magnétohydrodynamique dans un fluide non-newtonien saturant un
milieu poreux. Autre. Université de Picardie Jules Verne; Université de Tunis El Manar, 2018.
Français. �NNT : 2018AMIE0029�. �tel-03692082�

HAL Id: tel-03692082

https://theses.hal.science/tel-03692082v1
Submitted on 9 Jun 2022

HAL is a multi-disciplinary open access L’archive ouverte pluridisciplinaire HAL, est

archive for the deposit and dissemination of sci- destinée au dépôt et à la diffusion de documents
entific research documents, whether they are pub- scientifiques de niveau recherche, publiés ou non,
lished or not. The documents may come from émanant des établissements d’enseignement et de
teaching and research institutions in France or recherche français ou étrangers, des laboratoires
abroad, or from public or private research centers. publics ou privés.
 Thèse de Doctorat en cotutelle
Mention sciences de l'ingénieur
Spécialité Génie Mécanique

présentée à l'Ecole Doctorale en Sciences Technologie et Santé (ED 585)

de l’Université de Picardie Jules Verne

par

Cyrine CHAHTOUR

pour obtenir le grade de Docteur de l’Université de Picardie Jules Verne et de

l’Université de Tunis El-Manar, Faculté des Sciences de Tunis

Convection magnétohydrodynamique dans un fluide non-

newtonien saturant un milieu poreux

Soutenue le 21/11/2018, après avis des rapporteurs, devant le jury d’examen :

M. Rachid BENNACER, Professeur Rapporteur

M. Sadok BEN JABRALLAH, Professeur Rapporteur
M. Mohammed BENLAHSEN, Professeur Examinateur
M. Ezeddine SEDIKI, Professeur Examinateur
M. Hassen BEJI, Professeur Directeur de thèse
M. AmenAllah GUIZANI, Professeur Directeur de thèse
M. Haykel BEN HAMED, Maître de Conférences Co-directeur
2
 Remerciements

Ce travail a été effectué en cotutelle entre l’Université de Picardie Jules Verne en

France au sein du Laboratoire des Technologies Innovantes sous la direction de Monsieur
Hassen BEJI et l’Université de Tunis El-Manar en Tunisie au Laboratoire des Procédés
Thermiques sous la direction de Monsieur AmenAllah GIZANI. Ce travail n'aurait pu être
mené à bien sans l'aide de différents financeurs : CAMPUS France et le Ministère de
l'Enseignement Supérieur et de la Recherche Scientifique de la Tunisie.
J’aimerais remercier particulièrement mes directeurs de thèse :
Monsieur Hassen BEJI Professeur à l’Université de Picardie Jules Verne et directeur du
Laboratoire des Technologies Innovantes, de m’avoir accueillie au laboratoire et permis de
réaliser ce travail dans les meilleures conditions. Je le remercie sincèrement pour sa
disponibilité, ses conseils et précieuses remarques qui ont permis l’aboutissement de ce
travail. J’ai toujours trouvé le personnage de père et de facilitateur des difficultés courantes
que je trouve que ce soit au travail ou ailleurs. Je le remercie pour son aide, sa patience, pour
toutes les connaissances scientifiques qu’il a su me transmettre et surtout pour sa confiance.
Monsieur AmenAllah GIZANI, professeur à l’Université de Tunis El-Manar et directeur
du Laboratoire des Procédés Thermiques de m’avoir accueillie au laboratoire. Je lui suis
également reconnaissante pour le temps conséquent qu’il m’a accordé, ses qualités
pédagogiques et scientifiques, sa franchise et sa sympathie. C’est un homme que j’admire
beaucoup. Il n’a jamais épargné un effort pour m’aider sur tous les plans.
Je tiens également à remercier Monsieur Haykel BEN AHMED Maître de Conférences
à l’Université de Picardie Jules Verne pour son encadrement, son attention de tout instant sur
mes travaux et pour ses conseils avisés et son écoute qui ont été prépondérants pour la bonne
réussite de cette thèse. Son énergie et sa confiance ont été des éléments moteurs pour moi. J’ai
pris un grand plaisir à travailler avec lui. Ces quelques lignes ne peuvent témoigner de toute
ma gratitude et de mon amitié.
Je voudrais remercier les rapporteurs de cette thèse Monsieur Rachid BENNACER
Professeur à l’École Normale Supérieure de Cachan et Monsieur Sadok BEN JABRALLAH
Professeur à l’Université de Tunis, pour l’intérêt qu’ils ont porté à mon travail.
J'associe à ces remerciements Monsieur Mohammed BENLAHSEN Professeur à
l’Université de Picardie Jules Verne, pour avoir accepté d’examiner mon travail.
J’adresse de sincères remerciements à Monsieur Ezeddine SEDIKI Professeur à
l’Université de Tunis El-Manar pour ses remarques constructifs qui m’ont beaucoup aidé pour
améliorer ce manuscrit et pour avoir accepté d’examiner mon travail.
Je désire en outre remercier tous les membres du laboratoire LTI et au personnel du
département de génie civil pour leur sympathie, leur amitié. J’ai eu beaucoup de plaisir à
travailler avec eux. Le cadre de travail était idéal.
Ma plus profonde gratitude va à l’égard de ma chère mère qui a toujours cru en moi, en
ma capacité de surmonter les problèmes et les défis ; son soutien et son encouragement à
3
travers les mots doux et les prières m’ont permis de ne jamais baisser les bras et d’atteindre
mon but. Que Dieu lui accorde une longue vie pleine de santé, de la joie et du bonheur.
Je remercie ma chère sœur Chahrazed et je lui souhaite tout ce qu’il y a de meilleur.
Je souhaite remercier spécialement Faissal pour son soutien et sa patience tout au long
de la thèse.

4
 Résumé

Le travail de la thèse est consacré à une étude analytique et numérique de l'instabilité

thermique dans une cavité poreuse horizontale saturée par un fluide non-newtonien
électriquement conducteur de type loi de puissance. Un champ magnétique externe, uniforme
et constant est appliqué parallèlement à la gravité. Les parois parallèles aux plans (𝑜𝑦’𝑧′) et
(𝑜𝑥′𝑧’) sont maintenues imperméables et adiabatiques. Dans cette étude, nous avons envisagé
deux types de conditions aux frontières thermiques imposées aux parois actives de la cavité
(parallèles au plan (𝑜𝑥′𝑦’)) , les conditions aux frontières thermiques de type Neumann (flux
uniforme de chaleur) et de type Dirichlet (Température constante). De plus, un modèle
rhéologique de type loi en puissance a été utilisé pour modéliser le comportement non-
newtonien du fluide.

Pour l’étude de l’instabilité thermique, nous avons résolu les équations non linéaires
complètes en utilisant deux codes de calcul numériques, l’un est basé sur la méthode des
différences finies et l’autre sur la méthode des volumes finis. Par ailleurs, une solution
analytique, basée sur l’approximation de l’écoulement parallèle, a été développée dans le cas
de cavités minces (A>>1) soumises à un flux de chaleur constant. La comparaison des
valeurs de l’intensité de l’écoulement et le taux de transfert trouvés analytiquement et
numériquement sont en bon accord.

Nous avons utilisé l'analyse de stabilité linéaire pour prédire l'apparition des
mouvements convectifs en déterminant le nombre de Rayleigh critique en fonction de l'indice
de loi en puissance, caractérisant l’aspect non-newtonien du fluide, et du nombre de Hartmann
qui définit l’intensité du champ magnétique.

Les résultats obtenus ont permis de mettre en évidence l’influence du nombre

d’Hartmann ℋ et de l’indice de loi en puissance 𝑛, sur l’intensité de l’écoulement et sur le
transfert de chaleur. Ils ont montré que la présence du champ magnétique modifie les
résultats des travaux antérieurs concernant les fluides newtoniens et non newtoniens de type
loi de puissance. Par ailleurs, pour des champs magnétiques très élevés, il a été montré que la
dissipation d'énergie par effet Joule domine la dissipation d'énergie par contrainte de
cisaillement et confère au fluide un caractère non visqueux.

5
 Abstract

The work of the thesis is devoted to an analytical and numerical study of the thermal
instability in a horizontal shallow porous cavity saturated by a non-Newtonian power-law
liquid. An external, uniform and constant magnetic field is applied parallel to the gravity. The
walls parallel to the planes (𝑜𝑦′𝑧′) and (𝑜𝑥′𝑧′) are kept impermeable and adiabatic. In this
study, we considered two types of boundary conditions at the active walls of the cavity
(parallel to the plane (𝑜𝑥′𝑦 ′)), thermal boundary conditions of Neumann type (constant heat
flow) and Dirichlet type (constant temperature), a rheological model of the power law type
was used to model the non-Newtonian behavior of the fluid.
For the study of thermal instability, we have solved the complete nonlinear equations
using two numerical codes, one based on the finite difference method and the other on the
finite volume method. In addition, an analytical solution, based on the parallel flow
approximation, has been developed in the case of thin cavities (A>>1) subjected to a constant
heat flow. The comparison of the values of the flow intensity and the transfer rate found
analytically and numerically are in good agreement.
We used linear stability analysis to predict the onset of convective motion by
determining the critical Rayleigh number as a function of the power law index, characterizing
the non-Newtonian aspect of the fluid, and Hartmann number that defines the intensity of the
magnetic field.
The results obtained made it possible to highlight the influence of the number of
Hartmann and the index of law in power, on the intensity of the flow and on the transfer of
heat. They showed that the presence of the magnetic field modifies the results of previous
work concerning Newtonian and non-Newtonian fluids of the power law type. Moreover, we
also show that in the limit of very strong magnetic field, the dissipation of energy by Joule
effect dominates the dissipation of energy by shear stress and gives to the liquid an inviscid
character.

6
 Nomenclature
Symbole latin
𝐴 Rapport de forme selon l’axe 𝑥 de la cavité, (= 𝐻 ′ ⁄𝐿′ )

a Rapport de forme selon l’axe 𝑦 de la cavité, (= 𝑑 ′ ⁄𝐻 ′ )

⃗⃗⃗⃗
𝐵0′ Champ magnetique [𝑘𝑔𝐴−1 𝑠 −2 ]

𝐶 Gradient de température dans la direction 𝑥

𝑐 Célérité de la lumière [𝑚. 𝑠 −1 ]

𝑑 Diamètre des pores [𝑚]

𝑑′ Largeur de la cavité [𝑚]

Da Nombre de Darcy
E Module d’élasticité
⃗⃗⃗⃗
𝐸′ Champ électrique [𝑘𝑔. 𝑚. 𝐴−1 𝑠 −3 ]

𝑔 Accélération de la pesanteur [𝑚. 𝑠 −2 ]

𝐻′ Hauteur de la cavité [𝑚]

𝑛−1
𝐾 𝑑′
ℋ Nombre de Hartmann (= 𝛾𝐵′2 ( ) )
𝜂 𝜅

⃗𝐽′ Densité du courant électrique [𝐴. 𝑚−2 ]

𝐾 Perméabilité [𝑚1+𝑛 ]

𝑘𝑚 Conductivité thermique du milieu poreux [𝑊𝑚−1 𝐾 −1 ]

𝑘 Nombre d’onde

𝐿′ Longueur de la cavité [𝑚]

𝑙 Entier positf

ℒ Opérateur différentiel linéaire

𝑛 Indice de la loi en puissance, paramètre réel

𝑁𝑢 Nombre de Nusselt

𝑞′ Flux de chaleur constant par unité de surface [𝑊. 𝑚−2 ]

𝑃 Pression [𝑃𝑎]

7
𝒫 Péclet number (𝑢0′ 𝑑′ ⁄𝜅)
𝑛
𝑞′𝑑′ 𝐾 𝑑 ′
ℜ Nombre de Rayleigh (= 𝜌0 𝑔𝛽𝑇′ ( ) )
𝑘𝑚 𝜂 𝜅

𝑅𝑒 Nombre de Reynolds
s Variation temporelle des perturbations

𝑇′ Temperature [𝐾]

𝑉 Vitesse [𝑚. 𝑠 −1 ]

𝑢 Composante de la vitesse adimensionnelle, direction x

𝑣 Composante de la vitesse adimensionnelle, direction y

𝑤 Composante de la vitesse adimensionnelle, direction z

𝑥, 𝑦, 𝑧 Coordonnées cartésiennes à partir du centre de la cavité

𝑌 Vecteur des composantes de la perturbation

SYMBOLES GRECS

𝛽𝑇′ Coefficient d’expansion thermique[𝐾 −1 ]

𝜖 paramètre (<<1)

𝜀 Porosité du milieu poreux

𝜀̂ Paramètre du modèle de la loi en puissance en milieu poreux, [𝑃𝑎. 𝑠 𝑛 ]

𝜌 Masse volumique du fluide [𝑘𝑔. 𝑚−3 ]

𝜅 Diffusivité thermique [𝑚2 . 𝑠 −1 ]

𝜏 Contrainte de cisaillement [𝑃𝑎]

𝛾 Conductivité électrique [𝐴²𝑠 3 𝑚−3 𝑘𝑔−1 ]

𝛾̇ Taux de déformation [𝑠 −1 ]

𝜂 Facteur de consistance [𝑃𝑎. 𝑠 𝑛 ]

𝜓 Fonction du courant adimensionnelle

𝜗′ Potentiel électrique [𝑉]

𝜃 Distribution de température

𝜇𝑎 Viscosité apparente adimensionnelle

8
𝜇′ Viscosité dynamique d’un fluide newtonien [𝑃𝑎. 𝑠]

𝛼 Angle d’inclinaison

𝜈 Viscosité cinématique du fluide [𝑚2 . 𝑠 −1 ]

𝜔 Fréquence de l’onde

Φ Fonction des composantes spatiales

SOUSCRIPTS
0 Etat de référence
a Réfère à un échelle macroscopique

𝐵 Etat de base

𝐶 Réfère aux conditions critiques

𝑓 Réfère aux fluides

𝑠 Réfère aux solides

𝑚 Réfère à un milieu poreux

p Plastique

SUPERSCRIPT
′ Quantité dimensionnelle

𝑠𝑢𝑝 Supercritique

𝑠𝑢𝑏 Souscritique

̃ Petit écart

INDICES
𝑚𝑎𝑥 Maximum
VER Volume élémentaire représentatif

9
 Sommaire

Résumé ............................................................................................................................ 5

Abstract ........................................................................................................................... 6

Nomenclature .................................................................................................................. 7

Listes des figures........................................................................................................... 14

Listes des tableaux ........................................................................................................ 18

Contexte historique ....................................................................................................... 19

Introduction ................................................................................................................... 22

Chapitre 1 État de l’art ............................................................................................... 25

1.1 Généralités ...................................................................................................... 25

1.1.1 Magnétohydrodynamique MHD ............................................................... 25

1.1.2 Convection naturelle, forcée et mixte ....................................................... 25

1.1.3 Convection de Rayleigh-Bénard (RB) ...................................................... 26

1.1.4 Milieu poreux ............................................................................................ 27

1.1.5 Fluide newtonien ....................................................................................... 31

1.1.6 Fluide non-newtonien ................................................................................ 32

1.1.7 Modèles mathématiques ............................................................................ 39

1.2 Revue bibliographique .................................................................................... 40

1.2.1 Convection MHD dans les fluides newtoniens ......................................... 40

1.2.2 Convection MHD dans les fluides non-newtoniens .................................. 42

1.2.3 Convection naturelle dans les fluides non-newtoniens ............................. 43

1.3 Conclusion ...................................................................................................... 45

Chapitre 2 Formulation mathématique et analyse de la stabilité linéaire................... 46

10
 2.1 Description du problème considéré ................................................................ 46

2.1.1 Hypothèses simplificatrices....................................................................... 47

2.2 Formulation mathématique du problème ........................................................ 48

2.2.1 Forme adimensionnelle des équations....................................................... 52

2.2.2 Taux de transfert thermique ...................................................................... 54

2.3 Analyse de stabilité linéaire ............................................................................ 55

2.3.1 Solution de base (Solution de conduction) ................................................ 56

2.3.2 Équations aux perturbations ...................................................................... 57

2.3.3 Conditions aux frontières des perturbations : ............................................ 58

2.3.4 Décomposition en modes normaux ........................................................... 59

2.3.5 Stabilité marginale stationnaire 𝝎 = 𝟎 ..................................................... 61

2.4 Conclusion ...................................................................................................... 62

Chapitre 3 Convection MHD dans une cavité poreuse soumise à une condition du
flux constant 63

3.1 Introduction .................................................................................................... 63

3.2 Stabilité linéaire .............................................................................................. 64

3.2.1 Convection de Darcy-Bénard (DB) ........................................................... 65

3.2.2 Analyse des résultats linéaires................................................................... 67

3.3 Résolution analytique du problème non-linéaire : Approximation de

l’écoulement parallèle .......................................................................................................... 70

3.3.1 Equations de base simplifiées ................................................................... 72

3.3.2 Nombre de Nusselt .................................................................................... 73

3.3.3 Solutions de 𝝍, 𝑻 et 𝑵𝒖 ............................................................................ 73

3.3.4 Évaluation de la constante C ..................................................................... 77

3.4 Simulations numériques directes du problème non-linéaire........................... 78

3.4.1 Validation du code numérique .................................................................. 78

11
 3.5 Analyse des résultats non linéaires ................................................................. 80

3.5.1 Cas d’une couche poreuse saturée par un fluide newtonien ..................... 80

3.5.2 Cas d’une couche poreuse saturée par un fluide non-newtonien ............. 81

3.6 Conclusion ...................................................................................................... 88

Chapitre 4 Convection MHD dans une cavité poreuse soumise à une condition de
température constante ............................................................................................................... 90

4.1 Introduction .................................................................................................... 90

4.2 Stabilité linéaire .............................................................................................. 91

4.2.1 Convection de Darcy-Bénard (DB) ........................................................... 93

4.2.2 Analyse des résultats non linéaires............................................................ 94

4.3 Simulations numériques directes du problème non-linéaire........................... 99

4.3.1 Validation du code numérique ................................................................ 100

4.3.2 Analyse des résultats non linéaires.......................................................... 102

 Influence de l’indice de la loi en puissance 𝒏 .............................................. 102

 Influence du nombre de Hartmann 𝓗 .......................................................... 106

 Influence du nombre de Rayleigh ................................................................. 111

 Effet de Péclet ............................................................................................... 112

 Le comportement asymptotique ................................................................... 113

4.4 Conclusion .................................................................................................... 115

CONCLUSION GENERALE ET PERSPECTIVES.................................................. 117

Références bibliographiques ....................................................................................... 120

Annexe A .................................................................................................................... 126

A.1 Méthode des différences finies ........................................................................ 126

A.1.1 Algorithmes de résolution ......................................................................... 127

A.1.2 Discrétisation des équations ..................................................................... 128

12
 A.2 Méthode des volumes finis .............................................................................. 131

A.2.1 Résumé de la démarche numérique .......................................................... 132

A.2.2 Discrétisation des équations .......................................................................... 132

A.2.3 Accélération de la convergence .................................................................... 136

Annexe B .................................................................................................................... 137

B.1 Principe variationnel ........................................................................................ 137

B.1.1 Discrétisation de la formulation variationnelle ......................................... 137

B.1.2 Matrices d’assemblage .............................................................................. 139

B.1.3 Détermination du nombre Rayleigh critique ............................................. 141

13
 Listes des figures

Figure 0.1: Orientation des globules rouges sous l’effet d’un champ magnétique
appliqué parallèlement à l'écoulement du sang [9]. ................................................................. 21
Figure ‎0.2: Organigramme du rapport de la thèse……………………………………………23

Figure 1.1 : Cellules hexagonales observées dans des conditions similaires à celles de
Rayleigh Bénard obtenues par Koschmieder et Pallas [16]. .................................................... 26

Figure 1.2 : (a) Matériaux poreux granulaires utilisés dans l'industrie de la

construction, sphères de Liapor de diamètre 0,5 cm, (b) Calcaire broyé de diamètre 1 cm [17],
(c) Sable de Fontainebleau (phase solide en noir et espace des pores en blanc), (d) Mousse
métallique NiCrAl 𝜀 = 0.94, (e) Matériaux fibreux utilisés pour les couches de diffusion
d’une pile à combustible à membrane d’échange de protons et (f) Matériau céramique utilisé
dans des filtres à particules des moteurs diesels [18]. .............................................................. 28

Figure 1.3 : Le volume élémentaire représentatif (VER): la figure illustre la taille

intermédiaire par rapport aux dimensions du domaine d'écoulement et des pores [19]. ......... 29

Figure 1.4 : Cisaillement d’un fluide entre deux plaques parallèles espacées d’une
distance h ; la plaque supérieure se déplace à la vitesse 𝑢. ...................................................... 32

Figure 1.5 : Classification des fluides complexes [26]. ................................................ 33

Figure 1.6 : Rhéogrammes (a) et viscosités (b) en fonction du taux de cisaillement

obtenus pour trois exemples de liquides rhéofluidifiants : solution aqueuse de polyacrylamide,
du sang et d’une suspension de particules de polymère (latex) dans l’eau [27]. ..................... 35

Figure 1.7 : Comportement rhéologique des fluides non-newtoniens à comportement

indépendant du temps [25]. ...................................................................................................... 36

Figure 1.8 : Fluide antithixotrope : (a) au repos, (b) après une forte agitation [28]. .... 37

Figure 1.9 : Résistance du miel à l’écoulement en cisaillement (a) lorsqu'une

contrainte est appliquée (b) le miel présente une déformation qui augmente linéairement avec
le temps (c) [28]. ...................................................................................................................... 38

Figure 1.10 : Classification des fluides non-newtoniens [25]....................................... 38

Figure 1.11 : Une comparaison des résultats expérimentaux, obtenue par Nakagawa
[35], et théoriques, obtenue par Chandrasekhar [32], représentant l’évolution du nombre de
Rayleigh critique en fonction du nombre d’Hartmann. ............................................................ 41

14
 Figure 2.1 : Configuration géométrique dans le cas de condition du flux uniforme et
constant (a) et le cas de condition de température uniforme et constante (b) .......................... 47

Figure 2.2 : Représentation schématique du problème correspondant à l'état conducteur

de base ...................................................................................................................................... 56

Figure 2.3 : Expressions des conditions aux perturbations sur les frontières. .............. 59

Figure 3.1 : Représentation schématique du bilan énergétique dans le volume de

contrôle. .................................................................................................................................... 63

Figure 3.2 : Evolution du nombre de Rayleigh critique ℜ𝑐 en fonction de 𝑛 pour

différentes valeurs de Péclet 𝒫 et pour deux valeurs du nombre d’Hartmann ℋ : (a) 𝒫 ≤ 1
ℋ2 = 0, (b) ≥ 1 ℋ2 = 25, (c) 𝒫 ≤ 1 et ℋ2 = 0 et (d) 𝒫 ≥ 1 et ℋ2 = 25 (d)................... 68

Figure 3.3 : Evolution du nombre de Rayleigh critique ℜ𝑐 en fonction du nombre

d’Hartmann ℋ, pour 𝒫 → 0 et pour 𝑛 > 1 (fluides dilatants), 𝑛 < 1 (fluides
pseudoplastiques) et 𝑛 = 1 (fluide newtonien). ....................................................................... 70

Figure 3.4 : Diagrammes de bifurcation en termes de 𝜓𝑚𝑎𝑥 en fonction du nombre de

Rayleigh ℜ et de l'indice de la loi de puissance 𝑛 pour ℋ = 0. .............................................. 75

Figure 3.5 : Gradient de température horizontal 𝐶ℜ, ℋ pour les fluides newtoniens
(𝑛 = 1) soumis à un champ magnétique vertical. ................................................................... 80

Figure 3.6 : Isothermes prédites par la solution numérique des équations complètes
pour ℜ = 70, 𝑛 = 1 et pour différentes valeurs du nombre de Hartmann ℋ. ........................ 81

Figure 3.7 : Gradient de température horizontal 𝐶 en fonction du nombre de Rayleigh

ℜ et indice de loi de puissance 𝑛 pour ℋ = 0. ........................................................................ 82

Figure 3.8 : Evolution du nombre de Nusselt 𝑁𝑢 et l’intensité de l’écoulement

maximale 𝜓𝑚𝑎𝑥 en fonction du nombre de Rayleigh ℜ et de l'indice de loi de puissance 𝑛
pour ℋ = 0. ............................................................................................................................. 84

Figure 3.9 : ℜ𝑐 pour le début du mouvement convectif en fonction de l'indice de loi de

puissance 𝑛 pour ℋ2 = 2. La solution présentée en traits pointillés est le résultat de la
stabilité linéaire. ....................................................................................................................... 85

Figure 3.10 : 𝜓𝑚𝑎𝑥 et 𝑁𝑢 en fonction de ℜℋ2 et 𝑛 pour: (a) ℋ2 = 2 et (b) ℋ2 =

103. .......................................................................................................................................... 87

Figure 4.1 : Représentation schématique de la configuration géométrique en 3D. ...... 90

15
 Figure 4.2 : Naissance des structures convectives sous forme de rouleaux
longitudinaux 𝑅 ⊥ ou transversaux 𝑅 ∥ selon le type de fluide : (a) pseudoplastique 𝑛 = 0.5
et (b) dilatant 𝑛 = 1.5............................................................................................................... 93

Figure 4.3 : Evolution du nombre de Rayleigh critique ℜ𝑐 en fonction de 𝑛 pour

différentes valeurs de Péclet 𝒫 et pour deux valeurs du nombre d’Hartmann ℋ : (a) ≤ 1
ℋ2 = 0, (b) 𝒫 ≥ 1 ℋ2 = 25, (c) 𝒫 ≤ 1 et ℋ2 = 0 et (d) 𝒫 ≥ 1 et ℋ2 = 25 (d). ............. 96

Figure 4.4 : Evolution du nombre d’onde critique 𝑘𝑐 (a) et de la longueur d’onde

critique 𝜆𝑐 (b) en fonction de 𝒫, pour différentes valeurs de 𝑛 et pour ℋ2 = 2. ................... 97

Figure 4.5 : Evolution du nombre de Rayleigh critique ℜ𝑐 en fonction de ℋ, pour

𝒫 → 0 et pour 𝑛 > 1 (fluide dilatant), 𝑛 = 0.4 et 0.6 (fluides pseudoplastiques) et 𝑛 = 1
(fluide newtonien). ................................................................................................................... 99

Figure 4.6 : La variation du taux de transfert de chaleur Nu en fonction du nombre de

Rayleigh ℜ pour ℋ=0 et A=6. ............................................................................................... 101

Figure 4.7 : La variation du taux de transfert de chaleur Nu en fonction du nombre de

Rayleigh ℜ pour ℋ=0 et A=6. ............................................................................................... 102

Figure 4.8 : Lignes de courant au voisinage du point critique pour ℜ = 150, ℋ = 0,

n = 1 , et n = 1.6. ................................................................................................................. 103

Figure 4.9 : Vecteurs vitesses (gauche) et isothermes (droite) pour ℜ = 200, ℋ = 0

(a) 𝑛 = 1, (b) 𝑛 = 1.3, (c) 𝑛 = 1.6. ....................................................................................... 105

Figure 4.10 : Isothermes pour ℜ = 200, 𝑛 = 1 (gauche) 𝑛 = 1.5 (droite) et pour trois
valeurs de ℋ = 0,2 et 2. ........................................................................................................ 107

Figure 4.11 : Lignes de courant pour ℜ = 200, 𝑛 = 1 (gauche) 𝑛 = 1.5 (droite) et

pour trois valeurs de ℋ = 0,2 et 2. ........................................................................................ 108

Figure 4.12 : Nusselt locale dans le plan (X-Y) pour ℜ = 200, ℋ = 0, 2 et 3, 𝑛 = 1

(gauche) et 𝑛 = 1.5 (droite). .................................................................................................. 109

Figure 4.13 : Evolution du nombre de Nusselt moyen dans le plan (X-Y) en fonction
du nombre de Hartmann ℋ et du nombre de Rayleigh ℜ et pour 𝑛 = 1 et 𝑛 = 1.5. ............ 110

Figure 4.14 : Evolution de la température et de la composante de vitesse U en fonction

de z pour différentes valeurs de ℜ et 𝑛. ................................................................................. 111

Figure 4.15 : Structures convectives en convection naturelle 𝒫 = 0 et en convection

mixte 𝒫 = 10 pour ℜ = 200, ℋ = 1.................................................................................... 113
16
Figure 4.16 : 𝑁𝑢 en fonction de ℜℋ2 et 𝑛 pour: (a) ℋ2 = 2 et (b) ℋ2 = 104. ...... 115
Figure A.1 : Représentation du maillage de notre modèle physique

………………...1156

Figure A.2 : Volume élémentaire de contrôle…………………………………...…...133

17
 Listes des tableaux

Tableau 1.1 : Porosité de quelques matériaux [21]. ...................................................... 30

Tableau 1.2 : Viscosités dynamiques de quelques substances à la température ordinaire

[24]. .......................................................................................................................................... 32

Tableau 3.1 : Les expressions du nombre de Rayleigh critique 𝕽𝒄 et 𝒌𝒄 pour le

problème de DB en fonction de l’indice de loi en puissance 𝒏 et du nombre de Hartmann 𝓗
(Condition de Flux constant). .................................................................................................. 66

Tableau 3.2 : Comparaison des résultats numériques obtenus par la présente étude avec
ceux obtenus par Ben Khelifa et al. [61] en milieu poreux avec A=4 et ℜ=50. ...................... 79

Tableau 3.3 : Comparaison entre les résultats numériques et les résultats analytiques
dans le cas d’un fluide newtonien, pour ℜ =50, ℋ=1 et A=4. ................................................ 79

Tableau 4.1 : Les expressions du nombre de Rayleigh critique ℜ𝑐 et du nombre d’onde

critique 𝑘𝑐 pour le problème de DB (Condition de Dirichlet). ............................................... 93

Tableau 4.2 : Effet du nombre de Rayleigh, du nombre d’Hartmann 𝓗 et de l’indice de

la loi en puissance 𝒏 sur la valeur de Nusselt moyen. ........................................................... 112

Tableau A.1 : Identification des termes sources et des coefficients pour les équations
de transport sous la forme adimensionnelle………………………………………………... 133

Tableau A.2 : Expressions des coefficients 𝑎𝐸 , 𝑎𝑊 , 𝑎𝑁 , 𝑎𝑆 , 𝑎 𝑇 et 𝑎𝐵 ………….…… 135

18
 Contexte historique

L’élaboration de la théorie électromagnétique a été mise en jeu par un grand nombre de

physiciens: Oersted, Ampère, Arago, Faraday, Foucault, Henry, Lenz, Maxwell, Weber,
Helmholtz, Hertz, Lorentz et bien d’autres. Si elle débuta en 1820 avec Oersted, elle ne fut
mise en équations par Maxwell qu’en 1873 et ne trouva d’explication satisfaisante qu’en
1905, dans le cadre de la théorie de la relativité d’Einstein.

Il existe certains fluides qui conduisent l'électricité : le mercure, le sodium liquide, etc.,
lorsqu’ils sont placés dans un champ magnétique, leurs mouvement induisent des champs
électriques ces champs induisent à leur tour des courants électriques, qui interagissent avec le
champ magnétique et produisent des forces, qui affectent le mouvement du fluide. Mais il faut
attendre le XXème siècle pour qu’une théorie complète apparaisse, la théorie de la
magnétohydrodynamique. Le physicien suédois Hannes Alfvén fut le premier à employer le
terme magnétohydrodynamique (MHD) en 1942. Il a reçu le prix Nobel de physique en 1970
pour ses travaux sur ce sujet. La MHD a pour objet l’étude du mouvement d’un fluide
électriquement conducteur (liquide ou gaz ionisé appelé plasma) et soumis à un champ
magnétique. Il existe un couplage entre le mouvement du fluide conducteur et les effets
électromagnétiques.

En réalité, tous les fluides sont soumis généralement dans un champ magnétique
naturel. On peut citer l’exemple de champ magnétique des planètes ou des champs
magnétiques induits. Le couplage de l’écoulement des fluides et en particuliers les fluides non
newtoniens et les champs magnétiques est très peu étudié.

La magnéto-rhéologie a vu le jour au cours de la même période que la MHD, et a été

initialement introduite par Rabinow en 1948 [1]. Les fluides magnéto-rhéologiques (MR) sont
des matériaux sensibles aux champs magnétiques. Ils présentent des changements rapides,
spectaculaires, et réversibles des propriétés lorsqu’ils sont soumis à un champ d’induction.
D’après Rabinow [1], les fluides MR sont constitués de particules microscopiques contenant
du fer en suspension dans une matrice liquide. Lors de l'application d'un champ magnétique,
ces particules acquièrent un moment dipolaire et s’alignent pour former des chaînes.

19
La recherche scientifique dans ce domaine contribue d’une manière significative dans la
maîtrise et le contrôle actif de la viscosité. Les applications sont nombreuses, on peut citer par
exemple les embrayages [1], les visco-coupleurs, le confort de suspensions des véhicules de
luxe [2], le contrôle des vibrations sismiques, les prothèses, la précision de polissage [3], les
fluides de forage [4], etc. Dans quelques années le marché du transport, déjà orienté vers
l’énergie électrique, va certainement s’approprier ces technologies qui touchent les 7 × 109
humains. D’un autre côté, les fluides MR subissent des phénomènes de sédimentations dues à
la grande différence entre la densité des particules contenant les atomes de fer et la solution
newtonienne. La sédimentation dans les liquides magnéto-rhéologiques peut être totalement
inhibée en employant des fluides non-newtoniens [5,6] ou des gels thixotropes tels que des
nanoparticules de silice [7,8].

On peut citer un autre exemple d’écoulement de fluide non-newtonien tel que le sang humain
qui s’infiltre dans les tissus poreux soumis à un champ magnétique. Les globules rouges ont
une certaine élasticité permettant de faciliter l'écoulement du sang dans les veines et les
rétrécissements artériels qui minimisent le risque de thrombose. La viscosité du sang joue un
rôle majeur dans les maladies cardiaques. L’augmentation de la viscosité du sang, entraine
l’endommagement des vaisseaux sanguins, est responsable du risque d’AVC (Accident
vasculaire cérébral). Récemment, des travaux expérimentaux menés par Tao et Huang [9] ont
montré qu’en appliquant pendant une minute un champ magnétique parallèle à la direction de
l’écoulement sanguin et d’intensité de 1,3 tesla, on constate une diminution de la viscosité
sanguine diminue de 20 à 30 %. On remarque sur la figure 0.1, grâce à un champ magnétique
orienté parallèlement à l'écoulement du sang, la taille moyenne des particules en suspension
augmente et les globules s'agglomèrent en petites chaînes parallèles aux parois du vaisseau
(des chapelets de deux, trois, quatre, etc. globules), ce qui diminue la viscosité dans cette
direction. Au bout de deux ou trois heures, la viscosité remonte lentement et retrouve son état
initial.

20
Figure ‎0.1: Orientation des globules rouges sous l’effet d’un champ magnétique appliqué
parallèlement à l'écoulement du sang [9].

21
 Introduction

Au cours des dernières décennies, l’écoulement convectif de fluides non-newtoniens

dans un milieu poreux a suscité un intérêt considérable d’un grand nombre de chercheurs, non
seulement parce que les fluides non-newtoniens sont très présents dans la nature mais aussi
dans de nombreuses applications industrielles. Cependant, le sujet de recherche des fluides
non-newtoniens est assez complexe et nécessite des connaissances dans plusieurs disciplines
(mathématique, physique et chimique). La plupart de ces fluides possèdent des propriétés
rhéologiques très variées. Nous allons nous intéresser à l’étude de transfert de chaleur d’un
milieu poreux saturé par un fluide non-newtonien en particulier le fluide de type loi en
puissance. Ce type de fluide est caractérisé par sa viscosité qui varie avec le taux de
cisaillement. Dans la littérature, d’un point de vue rhéologique, il existe plusieurs lois
destinées à décrire la viscosité des liquides de type loi en puissance en fonction du taux de
cisaillement. Nous choisirons dans notre travail la loi en puissance proposée par Ostwald-de
Waele [10].

Ces dernières années, on s’est intéressé de plus en plus à l’étude de l’application

combinée de la magnétohydrodynamique et de milieux poreux. Étant donné que les domaines
d'applications sont aussi extrêmement vastes et diversifiés. Les équations de Darcy couplées
aux équations de Maxwell et appliquées aux fluides non-newtoniens, présentent l’un des
problèmes les plus riches en dynamiques, et met l’accent sur des interactions très peu étudiées
dans la littérature.

Contenu de la thèse

Dans nos travaux, nous allons nous intéresser à une étude numérique et analytique de
l’apparition de la convection de type Darcy-Bénard d’un fluide non-newtonien saturant une
couche poreuse très mince et soumise à un champ magnétique extérieur, constant et uniforme.
Pour la résolution numérique du problème, nous avons utilisé la méthode des différences
finies et la méthode des volumes finis. Pour l’étude analytique, nous avons considéré
l’approximation de l’écoulement parallèle.

22
 L’objectif consiste à étudier l’influence de l’intensité d’un champ magnétique
caractérisé par le nombre de Hartmann, sur la structure de l’écoulement et le transfert de
chaleur au sein d’un milieu poreux saturé par un fluide non newtonien. Les résultats obtenus
seront comparés à ceux trouvés dans les travaux antérieurs pour les fluides newtoniens.

Le présent travail est structuré en quatre chapitres, une conclusion générale et deux annexes.

Le premier chapitre est consacré d’une part à rappeler les notions fondamentales sur la
magnétohydrodynamique, la convection naturelle, forcée et mixte dans les milieux poreux et à
présenter les modèles rhéologiques de différents types de fluides non-newtoniens. D’autre
part, une étude bibliographique focalisée sur les effets magnétiques sur la convection naturelle
a été présentée. Une synthèse des études portant sur le transfert de chaleur en milieu poreux
pour les fluides newtoniens et non-newtoniens a été faite.

Le deuxième chapitre est consacré à la présentation des deux configurations

géométriques considérées et à l’exposition des équations modélisant le problème traité dans
cette thèse. Les équations adimensionnelles régissantes et les paramètres caractéristiques sont
dérivés. Nous présentons aussi dans ce chapitre l’analyse de la stabilité linéaire. L’objectif
essentiel étant de déterminer les valeurs critiques des paramètres de l'écoulement au-dessus
desquelles la convection apparaît. Nous définissons d’abord la solution de conduction
(solution de base) puis nous traitons les équations adimensionnelles régissant les
perturbations. Enfin, nous définissons la stabilité marginale stationnaire. Le résultat est
présenté sous la forme générale qui servira dans les chapitres suivants à examiner en détail
notre problème dans le cas du problème de Darcy-Bénard (DB).

Le troisième chapitre est consacré à l’étude du cas d’une cavité poreuse soumise à une
condition thermique de type Neumann. Dans une première partie, nous examinons en détail la
stabilité linéaire dans le cas du problème de Darcy-Bénard (DB) pour déterminer le nombre
de Rayleigh critique marquant la naissance des mouvements convectifs. Dans une deuxième
partie, nous développons une étude analytique basée sur l’approche de l’écoulement parallèle
afin de trouver les solutions analytiques de la vitesse et de la température. Afin de valider ces
solutions, nous faisons une comparaison avec une solution numérique, basée sur la méthode
de différences finies, développée dans une troisième partie. Nous montrons aussi dans ce
chapitre, l’effet du champ magnétique, de l’indice de la loi en puissance et du nombre de
Rayleigh sur les écoulements et les transferts de chaleur résultants.
23
 Le quatrième chapitre, il s’agit de la même démarche d’étude que le troisième chapitre
en considérant une cavité poreuse soumise aux conditions thermiques de type Dirichlet.
Dans la première partie, nous déterminons l’expression du nombre de Rayleigh critique qui
permet de prédire l’apparition de la convection au sein du fluide. En se basant sur les données
de la stabilité linéaire, une étude numérique tridimensionnelle, basée sur la méthode des
volumes finis, est présentée dans la deuxième partie. Une expérimentation numérique est
poussée afin de préciser l’influence de tous les paramètres de contrôle sur les écoulements et
les transferts de chaleur.

On termine par une conclusion générale mettant en évidence les résultats pertinents ainsi que
les perspectives de ce travail de recherche.

Figure ‎0.2: Organigramme du rapport de la thèse

24
Chapitre 1 État de l’art

1.1 Généralités
Dans ce chapitre nous présentons brièvement les définitions des différents mots clefs
afin de les familiariser avec le langage qui sera utilisé dans ce travail.

1.1.1 Magnétohydrodynamique MHD

La magnétohydrodynamique (MHD) est consacrée à l’étude des interactions entre le

champ de vitesse et le champ d’induction magnétique. Elle est donc une généralisation de
l'hydrodynamique (définie par les équations de Navier-Stokes) couplée à l'électromagnétisme
(définie par les équations de Maxwell).

On peut définir la MHD comme une discipline scientifique qui décrit le comportement d'un
fluide conducteur du courant électrique (liquide ou gaz ionisé appelé plasma) en présence
d’un champ électromagnétique. Les écarts de vitesse entre espèces chargées constituant le
plasma (électrons et ions) donnent lieu à des courants électriques. Ceux-ci vont donc
engendrer un champ électromagnétique qui aura à son tour une influence sur la dynamique du
plasma. On rencontre de tels plasmas dans les milieux naturels (magnétosphères planétaires,
étoiles, milieu interstellaire, pulsars, environnement de trous noirs) mais également dans
certaines réalisations expérimentales (propulseurs, lasers puissants, explosions de bombes,
machines à fusion).

1.1.2 Convection naturelle, forcée et mixte

La convection naturelle désigne l'ensemble des mouvements internes dans un fluide dus
aux variations locales de la masse volumique du fluide en fonction de la température et/ou de
la concentration d’un soluté. Cette différence de masse volumique implique une différence de
la poussée d'Archimède et donc crée un mouvement.

Contrairement à la convection naturelle, lorsque les mouvements convectifs apparaissent sous

l’action d’une source externe (exemple : une pompe), on parle alors de convection forcée où
l’écoulement persiste même en absence d’un gradient de température. La présence d’une

25
source externe et d’un gradient de température et/ou de concentration, donne lieu à une
convection mixte.

1.1.3 Convection de Rayleigh-Bénard (RB)

Décrite en 1916 par Rayleigh, la convection de Rayleigh Bénard (dite aussi: instabilité
de Rayleigh Bénard) est un phénomène physique qui a attiré l'intention de plusieurs
scientifiques dans le monde entier et qui a fait l’objet de nombreux travaux expérimentaux,
théoriques et notamment numériques.

la convection de Rayleigh Bénard apparait dans une couche mince de fluide maintenue entre
deux plaques horizontales portées à des températures différentes, la température de la plaque
inférieure étant plus élevée que celle de la plaque supérieure. Sous l'effet du chauffage de la
paroi inférieure, la masse volumique des particules fluides situées à proximité de la paroi
inférieure décroit et les particules remontent vers la paroi supérieure sous l'effet de la poussée
d'Archimède. Lorsque la différence de température entre les deux plaques dépasse une valeur
critique, l’énergie ne peut plus être évacuée par simple conduction par conséquence le
système devient instable et génère des rouleaux thermo-convectifs, appelés cellules de RB
(Figure ‎1.1).

Figure ‎1.1 : Cellules hexagonales observées dans des conditions similaires à celles de
Rayleigh Bénard obtenues par Koschmieder et Pallas [16].
L’équivalent du problème de RB en milieu poreux est celui de Darcy-Bénard (DB), également
connu sous le nom de problème de Horton-Rogers-Lapwood, étudié en premier par Horton et
Rodgers [11] et par la suite par Lapwood [12]. Ce problème a fait l'objet d'une grande quantité
d'œuvres durant les années soixante [13] et qui est encore étudié [14,15].
26
1.1.4 Milieu poreux

Les milieux poreux jouent un rôle important dans de nombreux secteurs d’activité
industrielles : géologie (exploitation pétrolière, structure des sols), mécanique et
biomécanique (structures alvéolaires), chimie (catalyse, séparation, milieux désordonnés),
anatomie (structure osseuse), agro-alimentaire (pates, pain), médecine, etc. Ce vaste domaine
d’applications pratiques a suscité l’intérêt de plusieurs chercheurs à étudier la convection dans
des espaces confinés.

 Définition

On appelle milieu poreux un corps composé d’une matrice solide à l’intérieur de

laquelle se trouvent des vides appelés pores de tailles et de géométries différentes, plus ou
moins interconnectés. On peut trouver des pores ne débouchant pas (pores aveugles) ou
occlus. L'interconnexion des pores permet l'écoulement d'un ou plusieurs fluides à travers le
matériau. Dans la situation la plus simple (écoulement monophasique), le vide est saturé par
un seul fluide. Dans l’écoulement diphasique, un liquide et un gaz partagent l'espace vide.

On trouve de nombreux matériaux naturels dans cette catégorie : le sable de plage, la plupart
des roches, le calcaire, le bois et les poumons humains. Certains matériaux artificiels
requièrent d’être poreux soit dans le processus de fabrication soit dans leur finalité pour jouer
un rôle de filtre ou apporter des propriétés macroscopiques particulières (Figure ‎1.2).

D’une manière générale, les milieux poreux sont caractérisés principalement par deux
propriétés macroscopiques liées entre elles :

(1) Porosité : le matériau doit contenir de petits espaces vides, appelés pores, délimités par
une matrice solide ;

(2) Perméabilité : le matériau doit être perméable à un écoulement de fluide (gaz ou liquide).

27
 Figure ‎1.2 : (a) Matériaux poreux granulaires utilisés dans l'industrie de la construction,
sphères de Liapor de diamètre 0,5 cm, (b) Calcaire broyé de diamètre 1 cm [17], (c) Sable de
Fontainebleau (phase solide en noir et espace des pores en blanc), (d) Mousse métallique
NiCrAl 𝜺 = 𝟎. 𝟗𝟒, (e) Matériaux fibreux utilisés pour les couches de diffusion d’une pile à
combustible à membrane d’échange de protons et (f) Matériau céramique utilisé dans des
filtres à particules des moteurs diesels [18].

 Volume Élémentaire Représentatif (VER)

Pour étudier l’écoulement dans un milieu poreux, il est nécessaire d’avoir un milieu
continu pour établir des relations aux dérivées partielles. Dans les milieux poreux, les
propriétés physiques (porosité, perméabilité, etc.) sont discontinues au niveau microscopique.
A cette échelle, les grandeurs locales (la pression et la vitesse) varient très irrégulièrement
d’un point à l’autre du domaine. On est donc amené à effectuer une moyenne spatiale de ces
grandeurs. Cette moyenne s’effectue sur de nombreux pores par l’intermédiaire d’un Volume
Élémentaire Représentatif (VER) du milieu (Figure ‎1.3). A l’intérieur du VER, les propriétés
moyennes des fluides et des matériaux sont supposées uniformes et continues. L'échelle de
longueur du VER de dimension caractéristique L doit être beaucoup plus grande que l'échelle
des pores de dimension 𝑑, mais considérablement plus petite que l'échelle de milieu poreux,
ou échelle macroscopique de dimension 𝐿𝑎 (Figure ‎1.3). Cette dimension 𝐿 doit satisfaire à la
double inégalité suivante : 𝑑 ≪ 𝐿 ≪ 𝐿𝑎

28
 Figure ‎1.3 : Le volume élémentaire représentatif (VER): la figure illustre la taille
intermédiaire par rapport aux dimensions du domaine d'écoulement et des pores [19].

 Porosité
La porosité d’un milieu poreux notée 𝜀, est une valeur numérique définie comme la fraction
du volume moyen du milieu occupé par les espaces vides sur le volume total soit :

volume du vide
𝜀= (1.1)
volume total
Cette grandeur 𝜀 varie donc entre 0 (solide plein) et 1 (volume complètement vide). Puisqu’il
s’agit d’un rapport des mêmes propriétés, la porosité n’a pas d’unité et elle est souvent
exprimée en pourcentage.

On peut distinguer les pores par la taille (Tableau ‎1.1). Rouquerol et al. [20] proposent la
classification suivante pour les porosités :

- Microporosité : relative à des pores dont le diamètre n'excède pas les 2 nanomètres ;
- Mésoporosité : relative à des pores dont le diamètre est compris entre 2 et 50
nanomètres ;
- Macroporosité : relatif à des pores dont le diamètre est supérieur à 50 nanomètres.

29
 Matériaux Porosité
Matériau mousseux 0,98
Fibre de verre 0,88 - 0,93
Fil à tisser 0,68 - 0,76
Grains de silice 0,65
Poudre d’ardoise noire 0,57 - 0,66
Cuir 0,56 - 0,59
Catalyseur 0,45
Granulé de pierres 0,44 - 0,45
Terre 0,43 - 0,54
Sable 0,37 - 0,50
Sphère bien empilée 0,36 - 0,43
Filtre de cigarettes 0,17 - 0,49
Briques 0,12 - 0,34
Poudre de cuivre 0,09 - 0,34
Pierre à chaud, Dolomite 0,04 - 0,10
Houille 0,02 - 0,07

Tableau ‎1.1 : Porosité de quelques matériaux [21].

 La loi de Darcy: Perméabilité

En considérant un écoulement monophasique à très faible vitesse dans un milieu
poreux, et en faisant des recherches sur l'hydrologie de l'alimentation en eau de Dijon,
l’ingénieur Henry Darcy a établi en 1856 une relation entre le débit et la différence de
pression appliquée. Cette relation est ensuite écrite, sous une forme généralisée, comme suit :

𝐾
⃗ =−
𝑉 ⃗⃗⃗⃗⃗
∇𝑃 (1.2)
𝜇
Cette relation est connue sous forme de loi de Darcy, où 𝜇 désigne la viscosité dynamique du
fluide. Le coefficient de proportionnalité 𝐾 est appelé perméabilité c’est une grandeur qui
exprime l'aptitude du milieu poreux à laisser passer le fluide. Elle est indépendante de la
nature du fluide mais elle dépend de la géométrie du milieu.

L’unité de la perméabilité K est le Darcy (1 Darcy est la perméabilité d’un milieu poreux de 1
cm² de section, 1cm de longueur, soumis à une différence de pression de 1 bar et traversé par
un fluide dont la vitesse de filtration est de 1 cm/s) 1 Darcy = 9.87x10-9 m² . La valeur de la
perméabilité est généralement déterminée expérimentalement. Plusieurs modèles empiriques
ont été proposées :

30
 - La relation de Kozeny-Carman [22] pour un milieu poreux consolidé constitué
d’éléments identiques de géométrie simple donne une estimation de la perméabilité K
qui s’écrit sous la forme suivante:
𝑑2 𝜀 3
𝐾= (1.3)
36𝐶0 (1 − 𝜀)²
où 𝑑 désigne une dimension caractéristique des éléments constituant la matrice
poreuse et 𝐶0 est une constante qui dépend des formes des grains (3.6 ≤ 𝐶0 ≤ 5).
- La relation d’Ergun [23], établit une expression semblable à l’équation de Kozeny-
Carman (en prenant 𝐶0 = 5) en considérant un écoulement unidirectionnel d’un fluide
incompressible, au sein d’une colonne poreuse constitué de particules sphériques de
diamètre 𝑑. La colonne est soumise à un gradient de pression. Cette relation s’écrit :
𝑑2 𝜀 3
𝐾= (1.4)
150(1 − 𝜀)²
Quand on modélise la matrice poreuse comme un faisceau de tubes capillaires
parallèles, l’expression de la perméabilité est :
𝑑2
𝐾=𝜀 2 (1.5)
32
Où 𝜀 = 𝑛𝜋𝑑²/4 et n le nombre de capillaires par unité de surface perpendiculaire au
sens d’écoulement et 𝑑 le diamètre des capillaires.

1.1.5 Fluide newtonien

On appelle fluide newtonien un fluide dont la contrainte de cisaillement 𝜏 est

proportionnelle au gradient de vitesse dans l’écoulement. La constante de proportionnalité est
appelée viscosité 𝜇. Dans le cas d’écoulements parallèles, le taux de cisaillement, aussi appelé
vitesse de déformation et noté 𝛾̇ . Ainsi, Newton a montré que lorsqu’on cisaille un fluide
(Figure ‎1.4), il apparait suite à la diffusion de quantité de mouvement :

– une force de résistance du fluide contre cette action de cisaillement ;

– cette force est proportionnelle au taux de cisaillement, ici U/h [1/s].

31
 u

ez

𝐻′
ex

Figure ‎1.4 : Cisaillement d’un fluide entre deux plaques parallèles espacées d’une distance h ;
la plaque supérieure se déplace à la vitesse 𝑢.
Ainsi la contrainte de cisaillement τ est définie comme suit:

𝑈
𝜏=𝜇 = 𝜇𝛾̇ avec 𝜇 = constante (1.6)
𝐻′
La viscosité dynamique 𝜇 (en Pa.s) est indépendante du taux de déformation et dépend
uniquement de la température et de la pression considérée. Cela signifie que le fluide continu
de s’écouler indépendamment des forces extérieures qui agissent sur lui. Par exemple, l’eau
est un fluide newtonien parce que sa viscosité reste constante quelle que soit la vitesse à
laquelle elle est agitée. Les solutions aqueuses, les huiles à faibles viscosités, la plupart des
solvants, l’air et de nombreux gaz sont des exemples de fluides newtoniens. Le Tableau ‎1.2
regroupe les viscosités dynamiques pour quelques produits courants.

Fluide 𝝁 (Pa.s)
air 2.10-5
pétrole 0,65.10-3
eau 10-3
mercure 1,56.10-3
sang 4-25.10-3
huile d’olive 0.1
miel 1-10
sirop d’érable 100
bitume 108

Tableau ‎1.2 : Viscosités dynamiques de quelques substances à la température ordinaire [24].

1.1.6 Fluide non-newtonien

Le fluide newtonien fournit le modèle le plus simple de comportement fluide. La

definition d’un fluide newtonien est assez restrictive : les contraintes de cisaillement sont

32
proportionnelles au gradient de vitesse, ce qui implique que la viscosite est independante du
temps et les contraintes s’annulent immediatement lorsque l’écoulement s’arrête.

En revanche, tous les fluides ne vérifient pas cette règle ou bien la vérifient partiellement. La
viscosité d’un fluide peut dépendre de deux variables : du cisaillement et du temps. Dans ce
cas, il s’agit d’un comportement non-newtonien. La description de ce comportement et son
interprétation, en relation avec la structure microscopique du fluide, constitue la discipline
appelée rhéologie. La Figure ‎1.5 présente la classification de grandes familles de liquides de
comportement non-newtonien.

Figure ‎1.5 : Classification des fluides complexes [26].

Chhabra et Richardson [25] proposent une classification selon trois grandes catégories
(Figure ‎1.10) ;

a. Fluides à comportement indépendant du temps

Ce sont des fluides dont les taux de déformation dépendent uniquement de la contrainte de
cisaillement appliquée. Ils sont généralement appelés fluides purement visqueux. Il existe
trois catégories de fluides :

33
  Fluides pseudoplastiques ou rhéofluidifiants
Les fluides pseudoplastiques sont caractérisés par une viscosité apparente qui diminue au fur
et à mesure que la contrainte de cisaillement à laquelle est soumis le fluide augmente. On
parle alors de comportement rhéofluidifiant (shear-thinning en anglais). Ce comportement est
également observé dans les suspensions de particules solides : plus le cisaillement est rapide,
plus ces particules s’orientent dans le sens de l’écoulement et leurs interactions de frottement
diminuent. Ce phénomène est observé dans les suspensions de vésicules déformables comme
le sang. La Figure ‎1.6 ci-dessous, montre les rhéogrammes (relation contrainte - taux de
cisaillement) et l’évolution de la viscosité en fonction du taux de cisaillement pour trois
fluides rhéofluidifiants. On retrouve ce type de fluide dans les solutions et les suspensions
aqueuses concentrées (certains polymères, purée de fruit, moutarde, Ketchup, etc.).

34
 Figure ‎1.6 : Rhéogrammes (a) et viscosités (b) en fonction du taux de cisaillement obtenus
pour trois exemples de liquides rhéofluidifiants : solution aqueuse de polyacrylamide, du sang
et d’une suspension de particules de polymère (latex) dans l’eau [27].

 Fluides dilatants ou rhéo-épaississants

Ce sont des fluides caractérisés par une viscosité qui augmente avec la contrainte de
cisaillement. Certaines suspensions concentrées (amidon de maïs par exemple) et le sable
mouillé ont un comportement rhéo-épaississant. En effet, le sable mouillé devient très dur
quand on appuie dessus rapidement. La même chose se produit avec de l’eau auquel on a
ajouté une bonne quantité de farine de maïs (Maïzena) : le remuer lentement avec une cuillère
est possible, mais dès qu’on la remue plus vite, la pâte devient très dure. En effet, le fluide
mouille les particules et les colle les unes aux autres, ce qui permet d’amortir la contrainte et
par conséquent d’augmenter la résistance au cisaillement. Le comportement rhéo-épaississant
est rarement observé comparé au comportement rhéofluidifiant.

 Fluides viscoplastiques
Ce sont des fluides caractérisés par une contrainte de cisaillement résiduelle 𝜏0 . Si la
contrainte appliquée au fluide est inférieure à 𝜏0 , aucune déformation ne se produit, aucun
mouvement n’apparaît au sein du fluide. On parle alors de fluide à seuil ou viscoplastique.
Comme il a été indiqué par Chhabra et Richardson [25] (Figure ‎1.7), ces fluides peuvent avoir
un comportement linéaire ou non linéaire. Les fluides viscoplastiques à comportement linéaire

35
sont appelés « fluides plastiques de Bingham ». Les fluides de Bingham sont des quasi-
solides. Exemples : les boues de forage de gisements pétroliers, la mayonnaise, le dentifrice
(qui commence à extruder qu’à partir d’une certaine contrainte appliquée sur le tube), etc.

La représentation la plus simple d’un fluide à seuil est le modèle de Bingham qui donne la
relation suivante entre la contrainte 𝜏 et le taux de cisaillement 𝛾̇ :

0 𝑠𝑖 𝜏 < 𝜏0𝐵
𝛾̇ = { (1.7)
𝜏 − 𝜏0𝐵 ⁄𝜇𝑝 𝑠𝑖 𝜏 ≥ 𝜏0𝐵
où 𝜏0𝐵 et la contrainte seuil et 𝜇𝑝 est appelé la viscosité plastique.

𝝉
Fluide viscoplastique

Fluide Bingham plastique

Fluide dilatant

Fluide newtonien
Contrainte de cisaillement [Pa]

Fluide
pseudoplastique

Taux de déformation [s-1] 𝜸̇

Figure ‎1.7 : Comportement rhéologique des fluides non-newtoniens à comportement

indépendant du temps [25].

b. Fluides à comportement dépendant du temps

La viscosité de ces fluides varie en fonction du temps et dépend des traitements antérieurs.

 Fluides thixotropes

36
Un fluide est dit thixotrope si sous contrainte constante, sa viscosité apparente diminue au
cours du temps. Comme la caractéristique du fluide dépend de son historique, les mesures
rhéologiques deviennent complexes dû aux effets d’hystérésis (effets de mémoire). Parmi les
fluides thixotropes, on peut citer le yaourt, les sables mouvants, certaines encres d'impression,
certaines boues de forage pétrolier, certains gels, etc.

 Fluides antithixotropes
L’antithixotropie est le phénomène inverse de la thixotropie : le fluide durcit par agitation, par
exemple : les suspensions aqueuses de gypse (Figure ‎1.8).

Figure ‎1.8 : Fluide antithixotrope : (a) au repos, (b) après une forte agitation [28].

c. Fluides viscoélastiques

Les fluides viscoélastiques sont des fluides qui présentent des caractéristiques à la fois
visqueuses et élastiques lorsqu'ils subissent une déformation. Les liquides visqueux, comme le
miel, résistent à un écoulement en cisaillement et présentent une déformation qui augmente
linéairement avec le temps lorsqu'une contrainte est appliquée (Figure ‎1.9). La pâte de
silicone est un exemple de fluide viscoélastique, elle peut réagir de deux manières différentes.
Une boule de silicone rebondit sur le sol comme une balle élastique (aux temps courts)
pourtant, si on pose cette boule sur une surface horizontale et on attend quelques minutes, on
voit la pâte de silicone s’étaler comme un fluide visqueux.

Le modèle le plus simple de fluide viscoélastique consiste à additionner les contraintes

d’origine élastique et les contraintes d’origine visqueuse :

37
 𝜏 = 𝜏𝑒𝑙𝑎𝑠𝑡. + 𝜏𝑣𝑖𝑠𝑞. = 𝐸𝛾̇ + 𝜇𝛾̇ (1.8)
où E est un module d’élasticité et 𝛾̇ est le taux de déformation.

Figure ‎1.9 : Résistance du miel à l’écoulement en cisaillement (a) lorsqu'une contrainte est
appliquée (b) le miel présente une déformation qui augmente linéairement avec le temps (c)
[28].

Fluide non-newtonien

dépendant du temps indépendant du temps Fluides viscoélastiques

Fluides
Fluides thixotropes
pseudoplastiques

Fluides Fluides dilatants

antithixotropes
Fluides
viscoplastiques

Figure ‎1.10 : Classification des fluides non-newtoniens [25].

38
1.1.7 Modèles mathématiques

Plusieurs modèles mathématiques ont été proposés dans la littérature pour modéliser le
comportement des fluides non-newtoniens. Dans cette thèse, on s’intéresse en particulier aux
fluides à comportement indépendant du temps et plus particulièrement aux fluides dilatants et
pseudoplastiques. Pour ces types de fluides, le modèle empirique de la loi en puissance
proposé par Ostwald-de Waele [10] est le plus utilisé dans la littérature. Cette loi est donnée
par l’expression suivante :

𝜏 = 𝜂(𝛾̇ )𝑛 (1.9)
ou encore

𝜇𝑎 = 𝜂(𝛾̇ )𝑛−1 (1.10)

où 𝜂 et n sont deux paramètres empiriques connus respectivement comme l’indice de
consistance et l’indice de comportement du fluide et 𝜇𝑎 est la viscosité apparente du fluide. Il
est à noter que ces paramètres n’ont aucune signification physique. La variation de l’indice n
permet de définir les classes suivantes :

- pour les fluides pseudoplastiques : 0 < n <1.

- pour les fluides newtoniens : n =1
- pour les fluides dilatants : n >1
En milieux poreux et en adoptant une loi de comportement de type Ostwald, une expression
empirique de la viscosité apparente 𝜇𝑎 du fluide a été proposée par Pascal [29], à savoir :

𝜇𝑎 = 𝜀̂(𝑢′2 + 𝑣′2 + 𝑤′2 )(𝑛−1)/2 (1.11)

Où 𝑢′, 𝑣′ et 𝑤′ sont, respectivement, les composantes de la vitesse selon 𝑥, 𝑦 et 𝑧.
L‘expression du coefficient 𝜀̂ est donnée par la relation suivante :

(𝑛+1) (𝑛−1)
𝜀̂ = 2𝜇⁄(8 2 (𝐾𝜀) 2 (𝑛/(1 + 3𝑛)𝑛 ) (1.12)

Le paramètre 𝜀̂ dépend du comportement rhéologique, de la viscosité dynamique 𝜇, de la

porosité 𝜀 et de la perméabilité du milieu poreux 𝐾.

39
1.2 Revue bibliographique
Dans cette section, nous proposons une revue bibliographique sur les écoulements
convectifs à travers les fluides newtoniens et non-newtoniens. Nous nous intéressons surtout à
la convection magnétohydrodynamique MHD dans les cavités poreuses.

L'étude de la convection dans un fluide électriquement conducteur saturant un milieu poreux

en présence d'un champ magnétique appliqué a eu un intérêt croissant ces dernières années vu
qu’elle revêt une importance particulière dans plusieurs disciplines comme l'astrophysique
(les taches solaires, le vent solaire, les éruptions solaires, etc.) et la mécanique des fluides
géophysiques (étude de la convection dans la croûte terrestre, qui se comporte comme un
milieu poreux).

1.2.1 Convection MHD dans les fluides newtoniens

Motivé par les applications de la MHD en astrophysique et en géophysique,

Chandrasekhar [30-32] a été le premier à travailler sur ce domaine. Il a étudié l’effet d’un
champ magnétique externe et constant sur l’apparition de l’instabilité thermique dans une
couche de fluide, électriquement conductrice et chauffée par le bas. Il a montré que le champ
magnétique a un effet inhibiteur sur l’apparition de l’instabilité et que les comportements
asymptotiques du nombre de Rayleigh critique ℜ𝑐 et le nombre d'onde critique 𝑘𝑐 associé, en
présence d’un champ magnétique très élevé (ℋ 2 → ∞), sont donnés par :

ℜ𝑐 → 𝜋 2 ℋ 2
} ℋ2 → ∞
𝑘𝑐 → (𝜋²ℋ 2 ⁄2)1⁄6

Où ℋ est le nombre de Hartmann définissant l'intensité du champ magnétique.

Parmi les premières expériences sur l'apparition de l’instabilité thermique [33-35] dans des
couches de mercure chauffées par le bas et soumises à un champ magnétique, Nakagawa [35]
a réalisé des expériences en utilisant des couches de mercure de 3, 4, 5 et 6 cm de profondeur
et un champ magnétique uniforme variant de 250 et 8000 gauss. Il a été capable de réaliser
l’expérience pour des valeurs de ℋ 2 dans la gamme de 102 − 106 . Les résultats
expérimentaux ont montré une croissance exponentielle du nombre de Rayleigh citrique en
fonction du nombre de Hartmann conforme à la prédiction théorique (Figure ‎1.11).

40
 Figure ‎1.11 : Une comparaison des résultats expérimentaux, obtenue par Nakagawa [35], et
théoriques, obtenue par Chandrasekhar [32], représentant l’évolution du nombre de Rayleigh
critique en fonction du nombre d’Hartmann.

Etant donné que l'application d'un champ magnétique peut aider à contrôler à la fois la
dynamique du fluide et le transfert de chaleur, le champ magnétique a été utilisé dans
l’industrie d’ingénierie telle que le confinement du plasma, la récupération de l'énergie
géothermique, l'extraction du pétrole, etc. En conséquence, une attention considérable a été
consacrée à l'étude de la magnétohydrodynamique (MHD) et transfert de chaleur dans les
milieux poreux au cours des dernières années.

La stabilité d’un écoulement newtonien en milieu poreux entre deux plaques horizontales
infinies, en présence d’un champ magnétique constant vertical, a été étudiée par Raptis et al.
[36]. Ils ont constaté que la vitesse diminue lorsque le nombre de Grashof augmente alors
qu’elle augmente avec la perméabilité.

Alchaar et al. [37] ont étudié la stabilité de l’écoulement sous l’effet d’un champ magnétique
de direction aléatoire en milieu poreux dans une cavité rectangulaire horizontale chauffée et
refroidie par des flux de chaleur constants, en utilisant le modèle de Brinkman. Il a été trouvé
qu’un champ magnétique vertical est le plus efficace pour retarder la convection.

41
L’étude numérique de la convection naturelle d’un fluide diamagnétique saturant une enceinte
rectangulaire poreuse et soumis à un champ magnétique de valeur élevée a été effectuée par
Zeng et al. [38]. Le modèle de Brinkman-Forchheimer a été utilisé pour modéliser le milieu
poreux. Les résultats montrent que la force du champ magnétique a un effet significatif sur le
champ d'écoulement et le transfert de chaleur.

Quelques études ont également traité le cas de cavités poreuses inclinées chauffées par les
côtés [39,40].

Robillard et al. [41] ont mené une recherche sur l’effet d’un champ électromagnétique
constant et vertical sur la convection naturelle dans une cavité rectangulaire verticale,
chauffée par un flux de chaleur uniforme imposé sur les parois verticales, remplie avec un
milieu poreux et saturée par un mélange binaire. Une solution analytique a été obtenue afin de
prédire l’impact du nombre de Hartmann sur l’écoulement.

Revnic et al. [42] ont étudié numériquement, en utilisant un schéma de différences finies
implicite, les effets d'un champ magnétique incliné et la génération de chaleur sur la
convection libre instable au sein d'une cavité carrée remplie d’un fluide saturant un milieu
poreux. Des températures constantes sont imposées sur les parois horizontales. Il a été
constaté qu'en augmentant le nombre de Hartmann le transfert de chaleur diffusif arrive à
prédominer, même pour les grandes valeurs du nombre de Rayleigh.

Alloui et al. [43] ont examiné analytiquement et numériquement l'effet d'un champ
magnétique vertical non constant sur la convection naturelle d'un fluide saturant une cavité
poreuse horizontale peu profonde. Le champ magnétique non constant est supposé être généré
par des courants électriques circulant parallèlement aux limites horizontales de la couche. Ils
ont démontré que l'effet retardateur de la force de Lorentz diminue la force du mouvement
convectif au sein de la couche et donc le transfert de chaleur par convection.

En outre, quelques résultats expérimentaux, traitant la MHD dans les milieux poreux, ont été
rapportés dans la littérature [44-46].

1.2.2 Convection MHD dans les fluides non-newtoniens

Considérant leurs applications pratiques dans divers procédés de l'industrie chimique et

des matériaux, l’écoulement des fluides non-newtoniens a été étudié par plusieurs chercheurs.

42
Parmi ces applications on peut citer : le refroidissement des réacteurs nucléaires par des
métaux liquides, la transformation des aliments et la propagation de contaminants dans
l'environnement.

La plupart des problèmes de MHD avec un fluide non-newtonien ont porté sur un fluide de
type spécial tel que micropolaire ou viscoélastique [47,48].

Récemment, Kefayati [49] a étudié numériquement le transfert de chaleur et la génération

d'entropie par convection naturelle et laminaire des nanofluides de type loi en puissance, en
présence d'un champ magnétique horizontal externe dans une cavité carrée. Les résultats
indiquent que l'augmentation de l'indice de loi de puissance entraine la diminution du transfert
de chaleur en absence du champ magnétique, en revanche, le transfert de chaleur augmente
avec l'augmentation de l'indice de loi de puissance en présence du champ magnétique.

Il existe un grand nombre de travaux étudiant l’influence du champ magnétique sur la

convection. La littérature dans ce domaine est concentrée autour de l’étude numérique, de
l’analyse des couches limites ou de la convection mixte dans les nanofluides. Le très petit
nombre de travaux concernant la convection dans un canal poreux, considère généralement un
gradient de température horizontal, avec la prise en compte des effets radiatifs et de l’effet
Hall.

1.2.3 Convection naturelle dans les fluides non-newtoniens

Depuis les travaux pionniers de Chen et Chen [50] beaucoup d'autres documents sur la
convection dans un milieu poreux saturé par un fluide en loi de puissance ont été publiés. Ils
sont concernés par une couche limite le long d'une plaque ou par des flux porteurs dans des
enceintes [51]-[53] ou convection forcée dans un canal [54-55].

D’un autre côté, Amari et al. [51] ont été parmi les premiers à étudier la convection naturelle
d’un fluide non-newtonien au sein d’une cavité poreuse horizontale. Des résultats à la fois
analytiques et numériques ont été rapportés pour le cas d’une couche horizontale chauffée par
le bas ou par les côtés par un flux de chaleur constant, sur la base du modèle de type loi en
puissance proposé par Pascal [29].

43
Une étude théorique et numérique de la convection naturelle induite par double diffusion dans
une cavité rectangulaire poreuse saturée par un fluide en loi de puissance et soumise à une
température constante a été menée par Getachew et al. [56].

Dans une étude récente, Nield [57] a montré que l'application de la théorie de la stabilité
linéaire au problème de Horton-Rogers-Lapwood HRL, pour un fluide de type loi en
puissance, conduit à un problème mathématique singulier sauf lorsque l'indice de la loi en
puissance est égal à l'unité, c'est à dire lorsque le fluide est newtonien. Afin de rendre le
problème non singulier, Nield [57] a suggéré que le terme de traînée dans la version non-
𝑛−1
⃗|
newtonienne de la loi de Darcy est proportionnel à |𝑉 ⃗ , où 𝑉
𝑉 ⃗ est la vitesse d'infiltration
𝑛
(m.s-1). Lorsque n>1 (fluide dilatant) et pour les grandes valeurs de |𝑉
⃗ |, le terme d'ordre |𝑉
⃗|

⃗ |. Donc la loi en puissance habituelle

devient dominant par rapport au terme proportionnel à |𝑉
avec l'indice 𝑛 a été retenue afin que l'ajustement empirique entre la nouvelle loi et l'ancienne
loi soit bon. En revanche, lorsque n<1 (fluides pseudoplastiques), l'approche décrite par cet
auteur [57] n’aboutit à aucune conclusion sur le début des instabilités.

Barletta et Nield [58] ont été motivés pour étudier davantage la nature de cette singularité. Ils
ont considéré le problème Prats [59] avec l'extension du problème de HRL en introduisant,
comme un flux de base, un certain nombre de Péclet 𝒫 caractérisant l'écoulement plutôt qu’un
débit nul. Le problème de HRL est récupéré dans la limite de 𝒫 tendant vers zéro. Les
résultats ainsi obtenus indiquent que le nombre de Rayleigh supercritique est infini pour un
fluide pseudoplastique et nul pour un fluide dilatant.

D’autres travaux ont également abordé le problème HRL avec un fluide non-newtonien en
considérant un fluide de type spécial tel qu’un nanofluide ou un fluide viscoélastique [60].

Ben Khelifa et al. [61] ont travaillé numériquement et analytiquement avec l’approximation
de l’écoulement parallèle de la convection naturelle dans une cavité poreuse peu profonde,
remplie d'un fluide binaire non-newtonien, soumis à des flux de chaleur et de soluté constants.
Les résultats obtenus indiquent que le nombre de Rayleigh critique pour l'apparition de la
convection est infini pour un fluide pseudoplastique et nul pour un fluide dilatant.

Les études mentionnées ci-dessus montrent que l’influence du champ magnétique sur la
convection d’un fluide newtonien a été suffisamment étudiée. Cependant, le cas des

44
écoulements MHD d’un fluide non-newtonien de type loi en puissance dans une couche
poreuse, chauffée par le bas et soumise à un champ magnétique uniforme et vertical n'a pas
été étudié jusqu'à présent.

1.3 Conclusion
Dans ce chapitre, nous avons présenté des généralités sur le phénomène physique étudié
et une synthèse bibliographique focalisée sur la convection naturelle des fluides newtoniens et
non-newtoniens en milieu poreux ainsi que les effets magnétiques sur la convection naturelle.
Nous avons montré que le cas des écoulements MHD d’un fluide non-newtonien de type loi
en puissance n'a pas été étudié jusqu'à présent. Dans ce travail, l’étude d’un écoulement
convectif dans un fluide non-newtonien de type loi de puissance saturant un milieu poreux de
Darcy et soumis à un champ magnétique uniforme et constant sera réalisée.

45
Chapitre 2 Formulation mathématique
et analyse de la stabilité linéaire

Dans ce chapitre, la première partie est consacrée à la présentation des deux

configurations géométriques considérées et à la mise en équations. Nous exposons les
équations modélisant le problème traité dans cette thèse, notamment le choix du système
d’adimensionnement. Cette formulation est fondée sur les lois de bilan ; de conservation de
masse ; de quantité de mouvement adaptée à un modèle rhéologique et d'énergie. En présence
d'un champ magnétique, une force apparaît au sein du liquide (force de Lorentz), une équation
supplémentaire est alors ajoutée pour compléter le système. Ensuite nous présentons les
conditions aux limites appropriées. Elles sont de plusieurs types et seront détaillées selon
l’application envisagée. Dans la deuxième partie, nous allons présenter la théorie de la
stabilité linéaire de la convection dans le cas général. L’objectif essentiel étant de déterminer
les valeurs critiques des paramètres de l'écoulement au-dessus de desquelles la convection
apparaît.

2.1 Description du problème considéré

Nous considérons le cas d’une cavité poreuse saturée par un fluide non-newtonien. La
configuration géométrique du problème étudié est représentée sur la Figure ‎2.1. Il s’agit d’une
enceinte tridimensionnelle de forme d'un parallélépipède de longueur L’ selon (Ox’), de
largeur l’ selon (Oy’) et d’hauteur H’ selon (Oz’). Ce système est supposé d’extension
horizontale infinie. Les parois parallèles au plan (𝑜𝑥′𝑦′) sont soumises à deux types de
conditions thermiques: à un flux de chaleur uniforme et constant 𝑞 ′ ou à une température
uniforme constante chaude 𝑇𝑐 = 𝑇0′ + ∆𝑇 ′ et froide 𝑇𝑓 = 𝑇0′ tandis que les autres parois sont
⃗ 0′ est appliqué
supposées adiabatiques. Un champ magnétique uniforme et constant 𝐵
parallèlement à la gravité 𝑔.

46
Figure ‎2.1 : Configuration géométrique dans le cas de condition du flux uniforme et constant
(a) et le cas de condition de température uniforme et constante (b)

2.1.1 Hypothèses simplificatrices

La modélisation du système étudié est basée sur les hypothèses simplificatrices

suivantes :

1) L'écoulement est laminaire.

47
2) Le fluide considéré est non-newtonien de type loi en puissance, homogène et
incompressible.

3) On suppose que la matrice poreuse est isotrope, homogène, immobile et en équilibre

thermodynamique avec le fluide.

4) Les forces d’inertie et les forces visqueuses, connues sous les noms de termes de
Forchheimer et de Brinkmann, sont négligées.

5) On néglige la propagation des ondes électromagnétiques ainsi que l’effet du champ

⃗⃗⃗⃗′ .
électrique 𝐸

6) La masse volumique varie linéairement avec la température. Selon l'approximation de

Boussinesq [62], elle s’écrit :

𝜌 = 𝜌0 [1 − 𝛽𝑇′ (𝑇 ′ − 𝑇0′ )] (2.1)

où 𝑇 ′ est la température du fluide en un point du système, 𝑇0′ est la température de référence,
𝜌0 représente la masse volumique de référence et 𝛽𝑇′ est le coefficient d’expansion thermique
défini par :

1 𝜕𝜌
𝛽𝑇′ = − ( ) (2.2)
𝜌0 𝜕𝑇 ′ 𝑃′
∆𝜌
L’approximation de Boussinesq est valable si la condition ≤ 0.1 est vérifiée.
𝜌

2.2 Formulation mathématique du problème

Le système d'équation qui gouverne l'écoulement laminaire et le transfert thermique
dans un milieu poreux en coordonnés cartésiennes s'écrit comme suit :

a. Equation de continuité

La conservation de la masse pour un fluide incompressible (ρ=constante) en régime

permanent est régit par l'équation de la continuité donnée par:

𝜕𝑢′ 𝜕𝑣 ′ 𝜕𝑤 ′
⃗ . ⃗⃗⃗⃗
∇ 𝑉′ = + + =0 (2.3)
𝜕𝑥 ′ 𝜕𝑦 ′ 𝜕𝑧 ′

⃗ exprimées selon :
où 𝑢′ , 𝑣 ′ et 𝑤 ′ sont les composantes du vecteur vitesse 𝑉

48
 ′
𝜕Ψ ′ ′
𝜕Ψ ′ (2.4)
𝑢 = et 𝑤 = − ′
𝜕𝑧 ′ 𝜕𝑥

où Ψ ′ représente la fonction de courant.

b. Equation de conservation de quantité de mouvement

On suppose maintenant que le fluide est électriquement conducteur soumis à un champ

magnétique ⃗⃗⃗⃗
𝐵0′ uniforme et constant. Le déplacement des charges dans le fluide induit par leur

mouvement un courant électrique de densité ⃗⃗𝐽′ normal à ⃗⃗⃗⃗ ⃗⃗⃗⃗′ . La force de Lorentz (force
𝑉 ′ et 𝐵 0

électromagnétique) ⃗⃗⃗
𝑓′ exercée sur la charge est :

⃗⃗⃗ = ⃗⃗𝐽′ × ⃗⃗⃗⃗

𝑓′ 𝐵0′ (2.5)

On définit la densité de courant électrique par la loi d’Ohm :

⃗⃗𝐽′ = 𝛾 (−∇ 𝑉 ′ × ⃗⃗⃗⃗

⃗ 𝜗 ′ + ⃗⃗⃗⃗ 𝐵0′ ) (2.6)

Le courant électrique est conservatif

⃗ . ⃗⃗𝐽′ = 0
∇ (2.7)
⃗ 𝜗 ′ = ⃗⃗⃗⃗
où 𝛾 est la conductivité électrique, 𝜗 ′ est le potentiel électrique et −∇ 𝐸 ′ est le champ
électrique associé.

D’après Chandrasekhar [32] ; étant donné que la vitesse des charges électriques du fluide 𝑉 ′
est très faible par rapport à la célérité de la lumière (𝑉 ′ << c), nous ne serons pas préoccupés
par les effets qui sont liés d'aucune façon à la propagation des ondes électromagnétiques, le
⃗ ) et le courant de déplacement peut être
⃗⃗⃗⃗′ ≈ 0
courant électrique varie très lentement (𝐸
négligé, d’où nous obtenons,

⃗⃗⃗⃗′ × ⃗⃗⃗⃗
⃗⃗𝐽′ = 𝛾 (𝑉 𝐵0′ ) (2.8)

Pour modéliser les écoulements de fluide dans les milieux poreux, on a choisi le modèle de
Darcy, le plus couramment utilisé dans les travaux antérieurs [63]. Ce modèle est valable
quand les conditions suivantes sont satisfaites :

49
 𝐾 𝑉̅ √𝐾
𝜀 ≤ 0.8 ; 𝐷𝑎 = ≤ 10−6 ; 𝑅𝑒 = <1 (2.9)
𝐻 ′2 𝜈

où 𝜀 est la porosité du milieu poreux, 𝐾 est la perméabilité, 𝐷𝑎 est le nombre de Darcy, 𝑅𝑒 est
le nombre de Reynolds, 𝑉̅ est la vitesse moyenne de filtration des particules fluides à travers
la matrice solide et 𝜈 est la viscosité cinématique du fluide.

En adoptant ce modèle et en présence d’un champ magnétique et gravitationnel terrestre,

l’équation de conservation de quantité de mouvement tenant compte de l’hypothèse de
Boussinesq s’écrit :

𝜂 ′ 𝑛−1 ′
⃗⃗⃗⃗ | ⃗⃗⃗⃗
|𝑉 ⃗ 𝑃′ + 𝜌0 𝑔(1 − 𝛽𝑇′ (𝑇 ′ − 𝑇0′ )) + ⃗⃗⃗
𝑉 = −∇ 𝑓′ (2.10)
𝐾

où 𝜂 (en Pa.sn) est le facteur de consistance et 𝑛 est l’indice de la loi en puissance. Pour n=1,
on a un fluide newtonien et le facteur de consistance coïncide avec la viscosité dynamique et
⃗ 𝑃′ désigne le gradient de pression.
∇

d. Equation de conservation de l’énergie

Dans un milieu poreux saturé par un fluide, le transfert de chaleur est assuré à la fois par la
phase fluide et la phase solide. Or ces deux phases ne possèdent ni la même capacité
thermique (respectivement (𝜌𝑐)𝑓 et (𝜌𝑐)𝑠 ), ni la même conductivité thermique
(respectivement 𝑘𝑓 et 𝑘𝑠 ). Donc Combarnous et Bories [64] avaient proposé un modèle de
deux équations d’énergie, établies à partir du premier principe de la thermodynamique,
décrivant l’évolution de la température des deux phases.

En prenant les moyennes des températures sur un VER, nous avons :

 pour la phase solide

𝜕𝑇𝑠′
(1 − 𝜀)(𝜌𝑐)𝑠 = ∇. (k ∗𝑠 ∇𝑇𝑠′ ) (2.11)
𝜕𝑡 ′
 pour la phase fluide
𝜕𝑇𝑓′
𝜀(𝜌𝑐)𝑓 ′ + (𝜌𝑐)𝑓 ⃗⃗⃗⃗
𝑉 ′ . ⃗∇𝑇𝑓′ = ∇. (k𝑓∗ ∇𝑇𝑓′ ) (2.12)
𝜕𝑡

50
où 𝑇𝑠′ et 𝑇𝑓′ représentent respectivement la température de la phase solide et fluide et les
scalaires k ∗𝑠 et k𝑓∗ représentent les coefficients de conductivité thermique équivalente et
dépendent des coefficients de conductivité thermique propre 𝑘𝑓 et 𝑘𝑠 et de la porosité 𝜀.

En supposant que l’équilibre thermique local est atteint (𝑇𝑠′ = 𝑇𝑓′ = 𝑇 ′ ) et après l’addition des
équations (2.11) et (2.12), nous déduisons le modèle de transfert de chaleur le plus
couramment utilisé pour les milieux poreux :

𝜕𝑇 ′
(𝜌𝑐)∗ ′
+ (𝜌𝑐)𝑓 ⃗⃗⃗⃗
𝑉 ′ . ⃗∇𝑇 ′ = ∇(𝑘 ∗ ∇𝑇 ′ ) (2.13)
𝜕𝑡
avec

(𝜌𝑐)∗ = 𝜀(𝜌𝑐)𝑓 + (1 − 𝜀)(𝜌𝑐)𝑠 (2.14)

et 𝑘 ∗ = k ∗𝑠 + k𝑓∗ (2.15)
sont, respectivement, la chaleur spécifique équivalente et la conductivité thermique globale.

En général, la conductivité thermique globale d'un milieu poreux dépend de manière

complexe de la géométrie du milieu. On la prendra constante dans la suite (𝑘 ∗ = 𝑘𝑚 ). On peut
en donner une approximation assez simple. Parmi les modèles les plus usuels [60], on
distingue :

- si la structure et l'orientation du milieu poreux sont telles que la conduction de chaleur

se déroule en série, avec un flux de chaleur passant à la fois par la phase solide et
fluide,

1−𝜀 𝜀
𝑘∗ = + (2.16)
𝑘𝑠 𝑘𝑓
- si la conduction thermique dans les phases solide et fluide se produit en parallèle,

𝑘 ∗ = (1 − 𝜀)𝑘𝑠 + 𝜀𝑘𝑓 (2.17)

En simplifiant par (𝜌𝑐)𝑓 , l’équation (2.13) peut s’écrire de la façon suivante :

𝜕𝑇 ′
𝜎 + ⃗⃗⃗⃗
𝑉 ′ . ⃗∇𝑇 ′ = 𝜅∇2 𝑇 ′ (2.18)
𝜕𝑡 ′
avec 𝜅 = 𝑘𝑚 ⁄(𝜌𝑐)𝑓 est la diffusivité thermique et 𝜎 = (𝜌𝑐)∗ ⁄(𝜌𝑐)𝑓 .

51
 e. Conditions aux frontières

Conditions aux frontières hydrodynamiques : conditions d'imperméabilité :

𝐿′
𝑥′ = ± 𝑢′ = 0
2
′
𝑑
𝑦′ = ± 𝑣′ = 0 (2.19)
2
′
𝐻′
{𝑧 = ± 2 𝑤′ = 0

Conditions aux frontières thermiques de type Neumann :

′
𝐿′ 𝜕𝑇 ′
𝑥 =± =0
2 𝜕𝑥 ′
′
𝑑 𝜕𝑇 ′
𝑦′ = ± =0 (2.20)
2 𝜕𝑦 ′
′
𝐻 𝜕𝑇 ′ 𝑞′
𝑧′ = ± = −
{ 2 𝜕𝑦 ′ 𝑘𝑚
Conditions aux frontières thermiques de type Dirichlet :

′
𝐿′ 𝜕𝑇 ′
𝑥 =± =0
2 𝜕𝑥 ′
′
𝑑 𝜕𝑇 ′
𝑦′ = ± =0
2 𝜕𝑦 ′
(2.21)
𝐻′
𝑧′ = − ′
𝑇 = 𝑇0′ + ∆𝑇 ′
2
′
𝐻
′
{𝑧 = + 2 𝑇 ′ = 𝑇0′

2.2.1 Forme adimensionnelle des équations

Pour faire apparaître les paramètres de contrôle du problème étudié il est nécessaire
d'introduire les variables caractéristiques du problème. La dimension 𝐻 ′ , entre les deux parois
horizontales actives de la cavité, est considérée comme longueur de référence. La diffusivité
thermique 𝜅est utilisée pour adimensionner la vitesse et la fonction du courant. La
température est adimensionnée respectivement par rapport à la différence caractéristique de
température ∆𝑇 ′ dans le cas de condition de température imposée et par rapport à 𝑞 ′ 𝐻′⁄𝑘𝑚
dans le cas de condition du flux imposée.

Les variables adimensionnelles utilisées sont les suivantes:

52
 𝑥′ 𝑦′ 𝑧′
Coordonnées spatiales (𝑥, 𝑦, 𝑧) = ( , , )
𝐻′ 𝐻′ 𝐻′
𝐻′
Vitesse (𝑢, 𝑣, 𝑤) = (𝑢′ ′
, 𝑣 , 𝑤′)
𝜅
′
𝑡𝜅
Temps 𝑡=
𝜎𝐻 ′ 2
(𝑇 ′ − 𝑇0′ )
𝑐ondition de température
∆𝑇 ′
Température 𝑇=
𝑘𝑚 (𝑇 ′ − 𝑇0′ )
condition du flux
{ 𝑞 ′ 𝐻′

Tableau 2.1 : Variables sans dimension

En introduisant les variables sans dimensions dans les équations (2.3), (2.10) et (2.18) et en
prenant le rotationnel de l’équation de mouvement (2.10), nous obtenons le système
d’équation adimensionnel suivant :

⃗∇. ⃗V = 0 (2.22)
𝑛−1 𝜕𝑇 𝜕𝑇 𝜕𝑣 𝜕𝑢 𝜕𝑢 𝜕𝑣 (2.23)
⃗ × (|V
∇ ⃗| ⃗V) = ℜ ( 𝑒𝑥 −
⃗⃗⃗ 𝑒𝑦 ) + ℋ 2 ( ⃗⃗⃗
⃗⃗⃗⃗ 𝑒𝑥 − 𝑒𝑦 + ( − ) ⃗⃗⃗
⃗⃗⃗⃗ 𝑒 )
𝜕𝑦 𝜕𝑥 𝜕𝑧 𝜕𝑧 𝜕𝑦 𝜕𝑥 𝑧
𝜕𝑇
⃗ .∇
+V ⃗ 𝑇 = ∇²𝑇 (2.24)
𝜕𝑡
Les nombres adimensionnels qui apparaissent dans les équations ci-dessus sont le nombre de
Rayleigh thermique ℜ et le nombre de Hartmann ℋ qui caractérisent respectivement les
effets convectifs et magnétiques au sein du fluide.

𝑛 𝑛−1
𝐾 𝐻′ 𝐾 𝐻′
ℜ= 𝜌0 𝑔𝛽𝑇′ ∆𝑇 ′ ( ) 2
ℋ = 𝛾𝐵 ( ) ′2
𝜂 𝜅 𝜂 𝜅 (2.25)

Sous leur forme adimensionnée, les conditions aux limites associées s'écrivent :

𝐴 𝜕𝑇
𝑥=± u=0 Ψ=0 ; =0
2 𝜕𝑥 (2.26)
𝑎 𝜕𝑇
𝑦=± 𝑣=0 Ψ=0 ; =0
2 𝜕𝑦
𝜕𝑇
1 = −1 𝐶𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 𝑑𝑢 𝑓𝑙𝑢𝑥
𝑧=± 𝑤 = 0 Ψ = 0 ; { 𝜕𝑧
2
{ 𝑇 = 1,0 𝐶𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛 𝑑𝑒 𝑡𝑒𝑚𝑝é𝑟𝑎𝑡𝑢𝑟𝑒

53
A est le rapport de forme de la cavité défini par : 𝐴 = 𝐿′ ⁄𝐻 ′ et 𝑎 est le rapport 𝑎 = 𝑑 ′ ⁄𝐻 ′

2.2.2 Taux de transfert thermique

L’étude du transfert de chaleur dans la cavité soumise à des conditions thermiques,

nécessite la détermination des taux de transfert de chaleur sur les parois actives de la cavité à
une position 𝑥 donnée. Pour l’ensemble des problèmes convectifs, les échanges de chaleur en
paroi se mesurent à l’aide du nombre de Nusselt 𝑁𝑢 qui représente le rapport entre le transfert
thermique total et le transfert par conduction. Si la conduction est le principal mode de
transfert, alors le nombre de Nusselt sera de l'ordre de l'unité. En cas de présence de
convection, le transfert thermique s'effectuera principalement par déplacement du fluide et la
valeur du nombre de Nusselt augmentera (𝑁𝑢 > 1).

Le nombre de Nusselt est défini comme suit :

 Cavité est soumise à un flux de chaleur

𝑞′ 1
𝑁𝑢(𝑥) = = (2.27)
𝑘∆𝑇 ′ /𝐻′ ∆𝑇(𝑥)
La valeur moyenne du nombre de Nusselt le long des parois actives est exprimée par :

1 𝐴/2
̅̅̅̅ =
𝑁𝑢 ∫ 𝑁𝑢(𝑥). 𝑑𝑥 (2.28)
𝐴 −𝐴/2

 Cavité est soumise à une température constante

𝜕𝑇
𝑁𝑢 = | (2.29)
𝜕𝑧 𝑧=±1/2
La valeur moyenne du nombre de Nusselt le long des parois actives est exprimée par :

𝐴/2 𝑎/2
𝜕𝑇
̅̅̅̅
𝑁𝑢 = ∫ ∫ | 𝑑𝑥𝑑𝑦 (2.30)
𝑥=−𝐴/2 𝑦=−𝑎/2 𝜕𝑧 𝑧=±1/2

54
2.3 Analyse de stabilité linéaire
Contrairement aux équations linéaires qui admettent des solutions bien déterminées,
l’équation de Navier-Stockes (2.10) est une équation non linéaire qui peut avoir plusieurs
solutions pour des paramètres d'écoulement (ℜ, ℋ,...) de valeurs identiques. La question est
donc de savoir quelle est la solution réelle ? La théorie qui traite ce genre de problème est la
théorie de la stabilité linéaire qui est basée sur l’introduction de notion de perturbations de très
faible amplitude dans les équations. Les équations deviennent linéaires pour les perturbations:
le produit de deux perturbations (ou plus) devient négligeable par rapport aux autres termes.

Le but de la stabilité linéaire est de déterminer la valeur critique des paramètres de

l'écoulement au-dessus de laquelle on est sûr que l'écoulement est instable. La présentation de
cette théorie est l’objet de cette partie.

L’analyse de stabilité linéaire se fait sur quatre étapes principales :

- étape 1: détermination du profil de base Φ(𝑥, 𝑦, 𝑧, 𝑡) = (𝑢, 𝑣, 𝑤, 𝑇).

- étape 2: écriture des équations de perturbations en Φ(𝑥, 𝑦, 𝑧, 𝑡) = (𝑈, 𝑉, 𝑊, 𝜃).

- étape 3: il est convenable de décomposer les solutions convectives spatio-temporelles en un

produit d’une exponentielle temporelle et d’une fonction spatiale. Les solutions Φ qui
vérifient les conditions aux limites sous la forme de modes normaux s’écrivent sous la forme
suivante,

̃ (x, y, 𝑧)𝑒 𝑠𝑡
Φ(x, y, z) = Φ (2.31)
̃ représente la fonction des composantes spatiales selon les directions confinées et 𝑠
Φ
représente la variation temporelle des perturbations et constitue la valeur propre du problème
linéaire.

- étape 4: Détermination du mode critique (𝑘𝑐 , ℜ𝑐 ), au-delà duquel l’écoulement cesse d’être
stable.

55
2.3.1 Solution de base (Solution de conduction)

Figure ‎2.2 : Représentation schématique du problème correspondant à l'état conducteur de

base
Nous supposons un écoulement de base stationnaire, parallèle et uniforme dans la
couche poreuse, à savoir,

𝑢𝐵 𝒫𝑐𝑜𝑠𝛼
𝑢 𝑣
⃗ 𝐵 ( 𝐵 ) = (𝒫𝑐𝑜𝑠𝛼 ) (2.32)
𝑤𝐵 0

Où 𝒫 est le nombre de Péclet obtenu par la division la vitesse d'écoulement de base 𝑢𝐵 par la
vitesse de référence 𝜅⁄𝐻 ′ et 𝛼 désigne l'angle d'inclinaison de la vitesse de base par rapport à
l'axe des 𝑥 (Figure ‎2.2). L'indice "B" signifie « solution de base ».

Physiquement, le nombre de Péclet, en rhéologie et dans l’étude de la diffusion, permet

d’évaluer l’effet couplé de la convection et de la diffusion. Lorsque 𝒫 ≫ 1, la convection
l’emporte sur la diffusion. Les particules sont donc transportées (advectées) par le fluide.
Dans le cas contraire, lorsque 𝒫 ≪ 1, la diffusion l’emporte sur la convection.
Mathématiquement et dans notre cas, l'utilisation du nombre de Péclet 𝒫, dans l’analyse de la
stabilité linéaire, est une technique artificielle pour éviter la difficulté mathématique à
résoudre le système caractéristique linéaire.

Le choix de l'hypothèse d'un angle d'inclinaison arbitraire α par rapport à l'axe des 𝑥 permet
d'étudier la stabilité linéaire de la solution de base en considérant des perturbations linéaires
𝜋
se propageant le long de l'axe des 𝑥 et en variant 𝛼 dans l'intervalle [0, ].
2

56
 𝜋
Ici, 𝒫 est positif. Avec 𝛼 ∈ [0, ], cela signifie que toutes les directions possibles de la vitesse
2

de base par rapport au vecteur d'onde de perturbation peuvent être autorisés [58]. Ainsi, le
profil de température correspondant à l'équation adimensionnelle (2.24) qui satisfait les
conditions aux limites (2.26) est,

𝑇𝐵 = 1 − 𝑧 (2.33)

2.3.2 Équations aux perturbations

Cette approche permet de décomposer la solution recherchée en une partie moyenne,

solution de l’écoulement dit de base (solution conductive), et d’une perturbation dont
l’amplitude est supposée très faible devant celle de l’écoulement de base (solution
convective).

La solution générale du système (2.22-24) sera recherchée sous la forme :

̃;
𝑢 = 𝑢𝐵 + 𝜖𝑈 𝑇 = 𝑇𝐵 + 𝜖𝜃̃ (2.34)
̃ = (𝑈
Telle que 𝜖 ≪ 1 est un petit paramètre et 𝑈 ̃, 𝑉̃, 𝑊
̃ ).

En remplaçant les relations (2.34) dans le système d’équation (2.22)-(2.24), en négligeant les
termes non linéaires, d’ordre supérieur à 1, et en tenant compte des équations (2.32) on
obtient le système d’équations aux perturbations linéarisées suivant :

̃ 𝜕𝑉̃ 𝜕𝑊
𝜕𝑈 ̃
+ + =0 (2.35)
𝜕𝑥 𝜕𝑦 𝜕𝑧
̃
𝜕𝑊 𝜕𝑉̃ ̃
𝜕𝑈 ℜ 𝜕𝜃̃ ℋ² 𝜕𝑉̃
− (𝑛𝑠𝑖𝑛²𝛼 + 𝑐𝑜𝑠²𝛼) − (𝑛 − 1) 𝑠𝑖𝑛𝛼 𝑐𝑜𝑠𝛼 = 𝑛−1 + 𝑛−1 (2.36)
𝜕𝑦 𝜕𝑧 𝜕𝑧 𝒫 𝜕𝑦 𝒫 𝜕𝑧
𝜕𝑈̃ 𝜕𝑊 ̃ 𝜕𝑉̃ ℜ 𝜕𝜃̃ ℋ² 𝜕𝑈̃
(𝑛𝑐𝑜𝑠²𝛼 + 𝑠𝑖𝑛²𝛼) − + (𝑛 − 1) 𝑠𝑖𝑛𝛼 𝑐𝑜𝑠𝛼 = − 𝑛−1 − 𝑛−1 (2.37)
𝜕𝑧 𝜕𝑥 𝜕𝑧 𝒫 𝜕𝑥 𝒫 𝜕𝑧
𝜕𝑉̃ 𝜕𝑈̃ ̃ 𝜕𝑉̃
𝜕𝑈
(𝑛𝑠𝑖𝑛²𝛼 + 𝑐𝑜𝑠²𝛼) − (𝑛𝑐𝑜𝑠²𝛼 + 𝑠𝑖𝑛²𝛼) + (𝑛 − 1) ( − ) 𝑠𝑖𝑛𝛼 𝑐𝑜𝑠𝛼
𝜕𝑥 𝜕𝑦 𝜕𝑥 𝜕𝑦
(2.38)
ℋ² 𝜕𝑈 ̃ 𝜕𝑉̃
= 𝑛−1 ( − )
𝒫 𝜕𝑦 𝜕𝑥

57
 𝜕𝜃̃ 𝜕𝜃̃ 𝜕𝜃̃
+ 𝒫 (𝑐𝑜𝑠𝛼 ̃ = ∇²𝜃̃
+ 𝑠𝑖𝑛𝛼 ) − 𝑊 (2.39)
𝜕𝑡 𝜕𝑥 𝜕𝑦
̃ , 𝜃̃ ) et les variables
̃, 𝑉̃, 𝑊
Il s’agit d’un système de 5 équations dont les inconnues sont (𝑈
(𝑥, 𝑦, 𝑧, 𝑡). L’écoulement de base est supposé connu et est noté (𝑢𝐵 , 𝑣𝐵 , 𝑤𝐵 , 𝑇𝐵 ). Pour
l’instant, aucune approximation n’a été faite concernant le champ moyen et les perturbations ;
le système conserve donc un caractère général, tridimensionnel et instationnaire. Néanmoins,
il est maintenant nécessaire de procéder à quelques hypothèses et simplifications pour pouvoir
mener à bien l’analyse.

D’après le théorème de squire [65], pour arriver à une perturbation tridimensionnelle il faut
passer par une perturbation bidimensionnelle; on peut alors se limiter à une perturbation
bidimensionnelle. Le système (2.35)-(2.39) peut donc être simplifié en 𝑦, en supprimant tous
les termes faisant intervenir toute dérivée par rapport à y. Donc cet écoulement peut s’écrire,

̃=𝑈
𝑈 ̃(𝑥, 𝑧, 𝑡), 𝑉̃ = 𝑉̃(𝑥, 𝑧, 𝑡), 𝑊 ̃ (𝑥, 𝑧, 𝑡), 𝜃̃ = 𝜃̃ (𝑥, 𝑧, 𝑡)
̃ =𝑊 (2.40)
Ainsi, nous obtenons

̃
𝜕2𝑊 ̃
𝜕2𝑊 ℜ 𝜕 2 𝜃̃
+ 𝑆² =
𝜕𝑥 2 𝜕𝑧 2 𝒫 𝑛−1 𝜕𝑥 2 (2.41)
𝜕𝜃̃ 𝜕𝜃̃ 𝜕 2 𝜃̃ 𝜕 2 𝜃̃
+ 𝒫𝑐𝑜𝑠𝛼 − ̃ =
𝑊 +
{ 𝜕𝑡 𝜕𝑥 𝜕𝑥 2 𝜕𝑧 2
Avec
ℋ²(𝑛𝑐𝑜𝑠²𝛼 + 𝑠𝑖𝑛²𝛼) + 𝑛𝒫 𝑛−1
𝑆² = + ℋ²𝒫1−𝑛 (2.42)
ℋ² + (𝑛𝑠𝑖𝑛²𝛼 + 𝑐𝑜𝑠²𝛼)𝒫 𝑛−1
et
(𝑛 − 1)𝑠𝑖𝑛𝛼 𝑐𝑜𝑠𝛼 (2.43)
𝑉̃ = − ̃
𝑈
ℋ²
𝑛𝑠𝑖𝑛²𝛼 + 𝑐𝑜𝑠²𝛼 + 𝑛−1
𝒫

2.3.3 Conditions aux frontières des perturbations :

Les conditions aux frontieres hydrodynamiques et thermiques pour la solution

convective (ou perturbation) diffèrent par rapport à celles des grandeurs de base. Elles sont
déterminées à partir des équations (2.26) et (2.33) et interpretées comme suit :

̃ = 0 : Condition d’impérméabilité
1) 𝑊
58
 ̃ ⁄𝜕𝑧² = 0 : Condition de glissement
2) 𝜕²𝑊
3) 𝜕𝜃̃ ⁄𝜕𝑧 = 0 : Condition du flux constant
4) 𝜃̃ = 0 : Condition de température constante

̃ 𝜕𝜃̃
𝜕²𝑊 = 0 𝐹𝑙𝑢𝑥 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡
𝑧 = ±0.5 ∶ ̃ =
𝑊 =0 ; { 𝜕𝑧 (2.44)
𝜕𝑧²
𝜃̃ = 0 𝑇𝑒𝑚𝑝𝑒𝑟𝑎𝑡𝑢𝑟𝑒 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒

̃
𝜕²𝑊 𝜕𝜃̃
̃ =
𝑊 =0 ; =0 𝑜𝑢 𝜃̃ = 0
𝜕𝑧² 𝜕𝑧

𝑧
𝑦

̃
𝜕²𝑊 𝜕𝜃̃
̃ =
𝑊 =0 ; =0 𝑜𝑢 𝜃̃ = 0
𝜕𝑧² 𝜕𝑧

Figure ‎2.3 : Expressions des conditions aux perturbations sur les frontières.

2.3.4 Décomposition en modes normaux

La troisième étape consiste à décomposer les solutions en modes normaux. L'analyse en

modes normaux est désormais classique pour une étude de stabilité linéaire.

En effet, comme le système est considéré comme infini, homogène et isotrope dans la
direction 𝑥 avec des conditions aux limites qui sont indépendantes de 𝑥 et 𝑡. Dans ces
̃ et 𝜃̃ sous la forme de mode de Fourier, dits
conditions, on peut chercher les solutions de 𝑊
modes normaux, suivant 𝑥 et oscillant dans le temps à la fréquence ω. Ces solutions sont
développées sur une base de fonctions adaptées à la géométrie du système.

Ainsi le système d’équation (2.41) admet des solutions de la forme :

59
 ̃ (𝑥, 𝑧, 𝑡) = 𝑊
𝑊 ̃ (𝑧)𝑒 𝑖(𝑘𝑥+𝜔𝑡)
{ (2.45)
𝜃̃ (𝑥, 𝑧, 𝑡) = 𝜃̃ (𝑧)𝑒 𝑖(𝑘𝑥+𝜔𝑡)

̃ (𝑧) et 𝜃̃ (𝑧) décrivent respectivement les variations inconnues par rapport

où les amplitudes 𝑊
à 𝑧 de la vitesse verticale et de la température et dépendent respectivement des conditions
hydrodynamiques et thermiques aux surfaces horizontales.

(𝜔, 𝑘) ∈ ℂ² désignent respectivement la pulsation complexe et le nombre d’onde complexe.

Ces nombres complexes peuvent se décomposer en : 𝜔 = 𝜔𝑟 + 𝑖𝜔𝑖 et 𝑘 = 𝑘𝑟 + 𝑖𝑘𝑖 avec :

 𝑘𝑟 : le nombre d’onde
 −𝑘𝑖 : le taux de croissance spatiale (à un temps fixé t, lorsque 𝑘𝑖 < 0 l’instabilité
s’amplifie dans l’espace pour 𝑥 > 0 sinon elle s’amortit).
 𝜔𝑟 : la pulsation ou la fréquence de l’onde
 𝜔𝑖 : le taux de croissance temporelle (à une position fixée. Lorsque 𝜔𝑖 > 0 l’instabilité
s’amplifie sinon elle s’amortit).
L’introduction de cette forme modale (2.45) dans le système d’équations précédent (2.41)
̃ (𝑧), 𝜃̃ (𝑧)) :
permet d’obtenir un nouveau jeu d’équations pour les fonctions d’amplitude (𝑊

ℜ
(𝑆²𝐷² − 𝑘²)𝑊̃ (𝑧) + 𝑘² 𝜃̃ (𝑧) = 0
{ 𝒫 𝑛−1 (2.46)
−𝑊̃ (𝑧) + (𝑘 2 − 𝐷2 + 𝑖Ω)𝜃̃ (𝑧) = 0

où D est une notation simplifiée pour 𝜕⁄𝜕𝑧 et

Ω = ω + 𝒫kcosα (2.47)
En écriture matricielle le système aux perturbations (2.46) devient :

ℒ. 𝑌 = 0 (2.48)
Où 𝑌 est le vecteur des deux composantes de la perturbation et ℒ est un opérateur différentiel
linéaire qui dépend des paramètres de contrôle ℜ, 𝑆, 𝒫, k, Ω et n :

60
 2 2 2
𝑘 2 ℜ𝒫1−𝑛 ] ̃
ℒ = [𝑆 𝐷 − 𝑘 ; 𝑌 = [𝑊 ] (2.49)
−1 𝑘 2 − 𝐷2 + 𝑖Ω 𝜃̃
Une solution non triviale n'existe que si le déterminant de ℒ conduit à une équation de
dispersion complexe écrite comme suit :

𝑑𝑒𝑡(𝑘, ω) = (𝑆²𝐷² − 𝑘²)(𝑘 2 − 𝐷2 + 𝑖Ω) + 𝑘²ℜ𝒫1−𝑛 = 0 (2.50)

La manière dont la solution du problème de stabilité doit être déterminée est claire. Pour un 𝑘
donné, nous devons déterminer la valeur caractéristique la plus basse pour ℜ; la valeur
minimale en fonction de 𝑘 de la valeur caractéristique obtenue est le nombre de Rayleigh
auquel l'instabilité se manifestera.

Deux approches de stabilité linéaire sont adoptées. Lorsque la perturbation est supposée être
étendue dans tout le système, une approche temporelle est suffisante. En revanche, la réponse
du système à une perturbation localisée nécessite une analyse spatio-temporelle.

2.3.5 Stabilité marginale stationnaire 𝝎 = 𝟎

Comme l'écoulement n'est pas tout le temps instable, on s'intéresse pour l'instant à la
naissance des premiers modes déstabilisant la solution de conduction. L'évolution temporelle
de la perturbation sera étudier en supposant :

𝑘 ∈ ℝ et 𝜔 ∈ ℂ

La croissance et la décroissance des perturbations dépendent de la partie imaginaire de la

pulsation 𝜔𝑖 ; pour 𝜔𝑖 < 0 la distribution est stable et pour 𝜔𝑖 > 0 la distribution instable.
Lorsque la perturbation n'est ni amplifiée ni atténuée, nous sommes dans la condition de
stabilité marginale qui est atteinte pour 𝜔𝑖 = 0, elle est donc utile pour déterminer les modes
les plus instables.

Dans ce cas la partie réelle 𝜔𝑟 contrôle le reste de la dépendance temporelle :

- 𝜔𝑟 = 0: Une instabilité stationnaire s'installe.

- 𝜔𝑟 ≠ 0: Un comportement oscillatoire apparaît.

61
2.4 Conclusion
Dans ce chapitre, nous avons présenté une modélisation mathématique du problème et
nous avons posé les équations dynamiques ainsi que l’équation de Darcy généralisée aux
fluides de type loi en puissance. Par ailleurs, les paramètres adimensionnées pertinents du
problème (ℜ, ℋ 2 ) ont été déterminés. Ce problème admet une solution dite solution de
conduction. A partir des équations développées, il est possible d’étudier la stabilité de l’état
de conduction en fonction des paramètres adimensionnés du problème. Pour cela nous avons
présenté la théorie de la stabilité linéaire de la convection dans le cas général. Nous avons
perturbé le système par un débit forcé et nous avons réduit le problème en un problème
caractéristique. Par la suite, nous allons étudier la stabilité linéaire de la convection naturelle
et mixte avec les conditions aux frontières thermiques de type Neumann au chapitre suivant
ensuite nous allons reprendre l’étude avec des conditions aux frontières thermiques de type
Dirichlet dans le quatrième chapitre.

62
Chapitre 3 Convection MHD dans une
cavité poreuse soumise à une
condition du flux constant

3.1 Introduction
Ce chapitre porte dans une première partie sur l’étude de la stabilité linéaire et dans une
deuxième partie sur l’application de l’approximation de l’écoulement parallèle au cas d’une
cavité à très grande extension, soumise à un flux de chaleur imposé sur les parois horizontales
tandis que les parois verticales sont adiabatiques. La configuration géométrique étudiée, les
conditions aux frontières et les axes de coordonnées sont montrés par la Figure ‎3.1. Afin de
valider ces solutions, nous avons fait une comparaison avec la solution numérique développée
dans une troisième partie.

L’objectif principal de cette étude est de prédire l’effet du champ magnétique et de

l’indice de loi en puissance n sur l’apparition de la convection, l’intensité des mouvements
convectifs et le transfert de chaleur résultant.

Figure ‎3.1 : Représentation schématique du bilan énergétique dans le volume de contrôle.

63
3.2 Stabilité linéaire
̃ (𝑧) et 𝜃̃ (𝑧) qui
Dans le cas de condition du flux imposée, les fonctions d’amplitude 𝑊
vérifient les conditions aux limites (2.44) sont :

𝑁 𝑁
̃ (𝑧) = ∑ 𝑊𝑙 cos(𝑙𝜋𝑧)
𝑊 𝑒𝑡 𝜃̃ (𝑧) = ∑ 𝜃𝑙 cos(2(𝑙 − 1)𝜋𝑧) (3.1)
𝑙=1 𝑙=1

Le problème aux valeurs propres (2.46) est résolu numériquement en utilisant le

développement (3.1) jusqu’à un ordre N, suffisamment grand pour assurer la convergence.

Pour déterminer le mode qui se déstabilise en premier, on doit calculer le nombre d’onde qui
minimise l’équation de dispersion (2.50). Comme ℜ est une fonction croissante de 𝑙, le
minimum est atteint pour 𝑙 = 1 et on obtient le mode critique (𝑘𝑐 ; ℜ𝑐 ), au-delà duquel
l’écoulement cesse d’être stable et fait place à la convection.
En appliquant le principe variationnel (Annexe B) au premier ordre (𝑙 = 1) et avec
Ω = ω + 𝒫kcosα, nous multiplions chaque variable par une fonction de test qui vérifie les
̃ est 𝑐𝑜𝑠(𝜋𝑧) et celle utilisée pour 𝜃̃ est
conditions aux limites. La fonction test utilisée pour 𝑊
1. Après intégration et en écriture matricielle le système (2.46) donne:
1 2
− (𝑆²𝜋² + 𝑘²) 𝑘²ℜ𝒫1−𝑛 𝑊 0
[ 2 𝜋 1
].[ ] = [ ] (3.2)
2 𝜃1 0
− 𝑘² + 𝑖Ω
𝜋
Le système (3.1) admet une solution non triviale (𝑘1 ≠ 0) si et seulement si 𝑑𝑒𝑡(𝑘, ω) = 0.
Donc l’équation de dispersion complexe (2.50) s’écrit sous la forme suivante :

1 4
𝑑𝑒𝑡(𝑘, ω) = − (𝑆 2 𝜋 2 + 𝑘 2 )(𝑘 2 + 𝑖Ω) + 𝑘²ℜ𝒫1−𝑛 = 0 (3.3)
2 𝜋²
Donc

𝜋² (𝑆 2 𝜋 2 + 𝑘 2 )(𝑘 2 + 𝑖Ω)
ℜ=𝒫 𝑛−1 (3.4)
8 𝑘²
On notera toutefois que le paramètre de contrôle ici est le nombre de Rayleigh
ℜ(ℋ, 𝒫, 𝑘, 𝜔, 𝑛) qui doit être réel (on impose 𝐼𝑚𝑔(ℜ) = 0), l’expression de ℜ s’écrit alors
sous la forme suivante :

64
 𝜋² 2 2
ℜ = 𝒫 𝑛−1 (𝑆 𝜋 + 𝑘 2 ) (3.5)
8

La valeur du nombre d'onde critique 𝑘𝑐 est obtenue en imposant 𝜕ℜ⁄𝜕𝑘 = 0. Ce qui donne
𝑘𝑐 = 0 (3.6)
On obtient ainsi le nombre ℜ𝑐 qui marque la naissance de la convection stationnaire (dite
aussi supercritique) :

𝜋4
ℜ𝑐 = 𝒫 𝑛−1 𝑆² (3.7)
8
On peut conclure qu’à partir des résultats ci-dessus, le nombre d’onde critique pour
l’apparition de la convection correspond à 𝑘𝑐 = 0, c’est-à-dire à une longueur d’onde infinie.
Cette valeur critique montre clairement que la convection se produit pour donner lieu à un
écoulement monocellulaire. En outre, les perturbations linéaires se propageant parallèlement à
l'axe des 𝑥 donc les perturbations les plus instables sont des rouleaux transversaux et 𝛼 = 0.
L’expression de S² se simplifie alors et nous obtenons :

𝑆² = 𝑛 + 𝒫1−𝑛 ℋ 2 𝑝𝑜𝑢𝑟 𝑛 ≠ 1 (3.8)

{
𝑆² = 1 + ℋ 2 𝑝𝑜𝑢𝑟 𝑛 = 1 (3.9)

3.2.1 Convection de Darcy-Bénard (DB)

Dans cette partie, nous examinons en détail notre problème dans le cas simple du
problème de Darcy-Bénard (DB) également connu sous le nom du problème de Horton-
Rogers-Lapwood; il fût traité tout d’abord par Horton et Rogers en (1945) [11] et
indépendamment par Lapwood en (1948) [12]. Ce problème consiste en une couche
horizontale infinie d’un milieu poreux isotrope saturé d’un fluide newtonien et chauffée par le
bas. Ce problème a fait l'objet d'une grande quantité d'œuvres depuis les années soixante [13-
15].

Dans la limite du nombre de Péclet tendant vers zéro (𝒫 → 0), la solution de base est relative
à un fluide au repos, et d’une stabilité marginale (𝜔 = 0), notre problème coïncide avec le
problème de DB. Dans ces conditions et à partir des équations (3.6) et (3.7), nous obtenons les

65
expressions du nombre de Rayleigh critique ℜ𝑐 et du nombre d’onde critique 𝑘𝑐 regroupés
dans le Tableau ‎3.1.

𝓟→ 𝟎 Flux constant
ℜ𝑐 𝑘𝑐
𝒏=𝟏 𝜋 4 ⁄8 (1 + ℋ²) 0
𝒏<𝟏 ∞ 0
𝒏>𝟏 ℋ²𝜋 4 ⁄8 0
Tableau ‎3.1 : Les expressions du nombre de Rayleigh critique 𝕽𝒄 et 𝒌𝒄 pour le problème de
DB en fonction de l’indice de loi en puissance 𝒏 et du nombre de Hartmann 𝓗 (Condition de
Flux constant).
Dans le problème de DB, les résultats présentés dans le Tableau ‎3.1 montrent que,
indépendamment de 𝑛, le nombre d’onde critique est égale à zéro (𝑘𝑐 = 0) ce qui donne lieu
à une longueur d’onde critique 𝜆𝑐 = 2𝜋⁄𝑘𝑐 → ∞ donc l’écoulement dans la cavité est
monocellulaire. En d’autre terme, pour une cavité peu profonde (A >> 1), l’écoulement est
essentiellement parallèle dans la direction horizontale. La solution numérique développée
dans la section suivante s'appuiera sur ces résultats.

D'autre part, nous déduisons que l’apparition de la convection dépend du nombre de

Hartmann ℋ pour le cas des fluides newtoniens (𝑛 = 1) ou fluides dilatants (𝑛 > 1). En effet
𝜋4
pour 𝑛 = 1 et ℋ = 0, on obtient ℜ𝑐 = ≈ 12.17. On retrouve le résultat trouvé par Nield
8

[13] pour la convection naturelle d'un fluide newtonien saturant un milieu poreux. Alors que
𝜋4
pour ℋ ≠ 0, l'apparition de l'instabilité convective s’effectue pour ℜ𝑐 > , elle dépend du
8

nombre de Hartmann ℋ. Et pour 𝑛 > 1 et ℋ = 0, on obtient ℜ𝑐 = 0 le système est

inconditionnellement instable. Alors que pour ℋ ≠ 0, l'apparition de l'instabilité convective
ne s’effectue qu’à partir de ℜ𝑐 = ℋ²𝜋 4 ⁄8.

Toutefois, pour le cas des fluides pseudoplastiques (𝑛 < 1), l’écoulement est
inconditionnellement stable. La valeur du nombre de Rayleigh critique pour l’apparition de la
convection, indépendamment de ℋ, se produit à ℜ𝑐 = ∞. Cependant, les résultats analytiques
et numériques des travaux de Amari et al. [51] montrent que la convection d’amplitude finie
est possible pour les fluides pseudoplastiques. Le critère ℜ𝑐 = ∞, prédit par la théorie de
stabilité linéaire n’est pas réaliste.

66
La simulation numérique ainsi que notre modèle basé sur la théorie non linéaire de
l’écoulement parallèle, vont être examinés dans les sections §‎3.3 et §‎3.4 pour voir s’il donne
un autre critère pour le début de la convection.

3.2.2 Analyse des résultats linéaires

 Influence du nombre de Péclet 𝓟 et de l’indice de loi en puissance 𝒏

Le comportement de ℜ𝑐 en fonction de 𝑛 pour les grandes et les faibles valeurs de 𝒫 est
représenté sur la Figure ‎3.2. L'utilisation du nombre de Péclet 𝒫, dans l’analyse de la stabilité
linéaire, est une technique artificielle pour éviter la difficulté mathématique à résoudre le
système caractéristique linéaire. Toutefois, cela donne un résultat important. Une faible valeur
de Péclet 𝒫 ≪ 1 donne une bonne approximation de la convection naturelle. On remarque sur
la Figure ‎3.2 qu'un débit filtrant très faible (𝒫 = 0.01) est suffisant pour voir un équilibre
instable.

Un résultat similaire a été rapporté récemment par Barletta et Nield [58] en considérant la
même configuration mais avec des conditions aux limites thermiques de type Dirichlet et sans
champ magnétique. La variation de ℜ𝑐 semble très similaire entre les Figure ‎3.2-(a) et (c)
lorsque 𝒫 ≤ 1 et entre les Figure ‎3.2-(b) et (d) lorsque 𝒫 ≥ 1.

Lorsque 𝒫 ≤ 1, les Figure ‎3.2-(a) et (c) montrent clairement que l’évolution de Rayleigh
critique ℜ𝑐 en fonction de 𝑛, pour un 𝒫 donné, est non monotone. Le débit filtrant 𝒫 donne
une simple transformation galiléenne dans le cas de fluide newtonien. En effet, pour 𝑛 = 1,
ℜ𝑐 est égal au seuil de convection naturelle pour toute valeur de Péclet. Cependant, pour les
fluides de loi de puissance, le nombre de Péclet rompt la transformation galiléenne connue en
raison de la dépendance de ℜ𝑐 en 𝒫 𝑛−1 (Eq (3.7)).

Pour 𝒫 = 1, nous obtenons ℜ𝑐 = 𝜋 4 (1 + ℋ²)⁄8 pour chaque 𝑛 ≥ 1. Pour de très grandes

valeurs de 𝒫 et pour 𝑛 < 1, l’apparition des mouvements convectifs se produit à un nombre
de Rayleigh sous critique ℜ𝑐 = 𝜋 4 ℋ 2 ⁄8 ≈ 304, avec un saut soudain vers l'infini quand
𝑛 ≥ 1. Par contre, pour les très faibles valeurs de 𝒫, l’évolution de ℜ𝑐 en fonction de 𝑛 a
tendance à être uniformément égale à zéro pour 𝑛 > 1, avec un saut soudain vers l'infini
lorsque 𝑛 = 1.

67
La raison physique de ce comportement particulier s’explique par l’adoption du modèle de loi
en puissance utilisé, Eq. (1.10). Il prédit une viscosité effective du fluide donnée par 𝜂(𝛾̇ )𝑛−1
qui dépend à la fois de l'indice de loi de puissance n et de la vitesse. En fait pour des faibles
valeurs de Péclet 𝒫 < 1, le taux de déformation 𝜸̇ est très faible à l'échelle des pores ceci rend
le fluide moins visqueux pour 𝑛 > 1. En conséquence, l'apparition de l'instabilité convective
aura lieu faiblement, entraînant ainsi une disparition de ℜ𝑐 . Tandis qu’à un débit élevé 𝒫 ≫
1, le taux de déformation 𝛾̇ est très grand et donc le fluide devient plus visqueux, par
conséquent l'apparition de l'instabilité convective sera complètement inhibée, ceci entraine
une valeur infinie de ℜ𝑐 . Les mécanismes de dissipation varient avec l'intensité du débit
d’écoulement 𝒫, d'où nous observons la rupture de la transformation galiléenne. C'est un
aspect très important des fluides de type loi en puissance.

30
10000
(a) H 2= 0 (b) 1000 100
H 2= 0
P=0.01
20
5

c 0.1 1 P=1
10
0.2
0.4
0.6 1
0.4 5
0.6 0.2 100
1 0.1 1000
0 0.01 10000

10000
328 (c) H 2= 25 (d) 1000
H 2= 25
100
P=0.01
320 5

c
312 0.1 1 P=1
0.2 1
0.6
0.4 0.4
0.6 5
304 0.2 100
1 0.1 1000
0.01 10000

0.0 0.5 1.0 1.5 2.00.0 0.5 1.0 1.5 2.0

n n

Figure ‎3.2 : Evolution du nombre de Rayleigh critique 𝕽𝒄 en fonction de 𝒏 pour différentes

valeurs de Péclet 𝓟 et pour deux valeurs du nombre d’Hartmann 𝓗 : (a) 𝓟 ≤ 𝟏 𝓗𝟐 = 𝟎, (b)
≥ 𝟏 𝓗𝟐 = 𝟐𝟓, (c) 𝓟 ≤ 𝟏 et 𝓗𝟐 = 𝟎 et (d) 𝓟 ≥ 𝟏 et 𝓗𝟐 = 𝟐𝟓 (d).
68
L’effet du nombre de Hartmann ℋ sur le nombre de Rayleigh critique ℜ𝑐 , pour l’apparition
des mouvements convectifs, est également présenté dans la Figure ‎3.2. On observe un grand
décalage sur l’axe de ℜ𝑐 lorsque ℋ² = 25. Ce décalage est dû au champ magnétique qui
retarde le seuil de convection d'un facteur de 12ℋ² tel qu'obtenu dans l'analyse de stabilité
linéaire.

 Le comportement asymptotique

A partir des résultats de la stabilité linéaire, nous avons présenté l’influence du nombre
de Hartmann sur la variation du nombre de Rayleigh critique ℜ𝑐 sur la Figure ‎3.3, pour
différents type de fluides. On peut constater à partir de la figure l’effet inhibiteur du champ
magnétique sur l’apparition de l’instabilité. L’évolution du nombre de Rayleigh critique ℜ𝑐
pour les fluides dilatants, pseudoplastiques et newtoniens présentent les mêmes
caractéristiques générales. Ce comportement particulier est dû à la dépendance asymptotique
du nombre de Rayleigh critique en ℋ 2 , qui découle de l’équation (3.7) ;

𝜋 4 (1 + ℋ 2 )⁄8 pour n = 1 (3.10)

ℜ𝑐 → {
𝜋 4 (𝑛𝒫 𝑛−1 + ℋ 2 )⁄8 pour n ≠ 1 (3.11)
La dépendance de ℜ𝑐 en ℋ², lorsque ℋ → ∞, obéit à la loi de proportionnalité ℜ𝐶 ⟶ 12ℋ 2 .
Dans cette limite, en remplaçant les équations (2.25) dans l’équation (3.7), on peut déduire
que le flux de chaleur critique 𝑞𝐶′ est donné par;

𝑘𝑚 𝛾𝜅𝐵2
𝑞𝐶′ = 12 (3.12)
𝜌0 𝑔𝛽𝑇′ 𝑑′ 2
Nous notons que le flux de chaleur critique 𝑞𝐶′ devient indépendant de la viscosité, de la
perméabilité et de l’indice de loi de puissance. Il dépend plutôt de la conductivité électrique,
de l'intensité du champ magnétique, de la gravité, de la hauteur de la couche et des propriétés
thermiques du liquide. Par conséquence, la dissipation d'énergie par effet Joule doit
prédominer la dissipation d'énergie par contrainte de cisaillement. La présence d'un fort
champ magnétique confère au nouveau comportement liquide, et on ne peut plus le décrire
comme newtonien ou non-newtonien, puisque la viscosité n'a pas d'importance.

69
 105

104

0.25

103
c 0.5

0.75
102
n=1
12 H 2
n>1
101
100 101 102

Figure ‎3.3 : Evolution du nombre de Rayleigh critique 𝕽𝒄 en fonction du nombre

d’Hartmann 𝓗, pour 𝓟 → 𝟎 et pour 𝒏 > 𝟏 (fluides dilatants), 𝒏 < 𝟏 (fluides
pseudoplastiques) et 𝒏 = 𝟏 (fluide newtonien).

D’autre part, La Figure ‎3.3 montre que les courbes dilatantes adhèrent rapidement à la ligne
asymptotique en raison de l'effet d'épaississement par cisaillement, mais les fluides
pseudoplastiques nécessitent un champ magnétique plus puissant pour inhiber l'effet de la
viscosité et de la rhéologie non-newtonienne. Par exemple, dans le cas de 𝑛 = 0.75, un
nombre de Hartmann d'environ 10 est suffisant pour créer un mouvement non visqueux, mais
ℋ = 50 et ℋ = 80 sont respectivement nécessaires pour les cas 𝑛 = 0.5 et 𝑛 = 0.25.
Lorsque ℋ < 100, la dissipation par viscosité et effet Joule agissent ensemble.

3.3 Résolution analytique du problème non-linéaire :

Approximation de l’écoulement parallèle
Dans la partie précédente, nous avons vu que l’écoulement convectif dans le système
considéré (Figure ‎2.1-(a)) est un écoulement monocellulaire et parallèle aux longues parois

70
(𝑘𝑐 = 0, 𝜆𝑐 = ∞). Ainsi dans cette partie, on présente le développement d’une solution
analytique basée sur le concept de l’écoulement parallèle, valide pour le cas limite d’une cavité
de grande extension (A>>1).

Nous avons établi dans le Chapitre 2 le système d’équations qui régit les phénomènes
d’écoulements et de transferts thermiques au sein d’un milieu poreux saturé par un fluide non-
newtonien et en présence d’un champ magnétique. Ce système est composé d’équations
différentielles aux dérivées partielles non linéaires, et la résolution d’un tel système de manière
analytique demeure extrêmement compliquée, sauf pour certaines hypothèses pour lesquelles
elles se simplifient considérablement. Dans le cas de cavités de grande extension (A >>1),
lorsque les parois actives sont exposées à un flux constant de chaleur, il est possible d’obtenir,
en utilisant le concept de l'écoulement parallèle, des solutions analytiques approximatives. Le
modèle résultant donne un aperçu général sur l’écoulement et le champ de température ainsi
qu’une prédiction de certains phénomènes physiques impliqués dans la convection naturelle.

L’approche de l’écoulement parallèle suppose, pour les cavités à grand rapport de forme
( 𝐴 >> 1) ou( 𝐴 << 1), que l'écoulement résultant est parallèle relativement aux longues
parois. Ceci permettra de négliger la composante de vitesse perpendiculaire à ces parois. Cette
approche a été utilisée dans le passé pour étudier la convection naturelle dans plusieurs
configurations en milieux fluides ou poreux confinés soumises à différentes conditions aux
frontières, en utilisant une méthode d’expansion asymptotique, par Cormak et al. [66] et, une
méthode de type intégrale, par Bejan et Tien [67]. D’autres approches similaires ont été
utilisées successivement par [68]-[74]. En utilisant le modèle de Darcy simplifié, pour une
couche poreuse verticale saturé par un fluide de type loi en puissance, le problème a été résolu
analytiquement par Bian et al. [75] dans la limite d'une couche mince. Dans une étude récente,
Ben Khelifa et al. [58] ont mené une investigation sur la convection naturelle dans une cavité
poreuse horizontale peu profonde saturée par un fluide binaire non-newtonien, sujette à des
flux constants de chaleur et de soluté. Des solutions dépendant de l’indice de loi en puissance
n ont été obtenues analytiquement. Un bon accord a été trouvé entre le modèle analytique et la
solution numérique des équations complètes régissant le problème.

Nous citons, de plus, des publications récentes effectuées par différents auteurs tels que Alloui
et al. [76] qui ont également utilisé cette méthode pour l’étude de la convection MHD dans un
fluide micropolaire. Ainsi que celle effectuée par Alloui et al. [43] sur la convection naturelle

71
dans une couche poreuse horizontale chauffée par le bas et soumise à un champ magnétique
non-constant. L’approximation de l’écoulement parallèle a été utilisée par ces auteurs pour
prédire analytiquement le nombre de Rayleigh critique pour le début de la convection.

Dans notre cas, la cavité présente un grand rapport de forme, 𝐿′ >> 𝐻′ , (voir Figure ‎3.1) les
lignes de courant au centre de la cavité deviennent parallèles à l’axe des 𝑥. Ceci permet de
négliger la composante de la vitesse perpendiculaire aux parois horizontales donc la vitesse
devient fonction de la coordonnée 𝑧 seulement, de telle sorte que :

𝑢(𝑥, 𝑧) ≈ 𝑢(𝑧) (3.13)

𝑤(𝑥, 𝑧) ≈ 0 (3.14)
Le profil de la température est alors donné par la somme d’un terme définissant une variation
longitudinale linéaire en 𝑥 et d’un autre terme donnant la distribution transversale en 𝑧 :

𝑇(𝑥, 𝑧) ≈ 𝐶𝑥 + 𝜃(𝑧) (3.15)

Où 𝐶 est une constante qui exprime le gradient de température horizontal, selon la direction 𝑥.
La valeur de C peut être déterminée en intégrant l’équation d’énergie sur un volume de
contrôle arbitraire près des régions d'extrémité de la cavité (Figure ‎3.1), tout en tenant compte
des conditions aux frontières (2.26). Ce qui conduit après manipulation à l’expression
suivante :

0.5
𝐶=∫ 𝑢𝑇𝑑𝑧 (3.16)
−0.5

3.3.1 Equations de base simplifiées

Dans cette partie, nous présentons les solutions en fonction 𝜓(𝑧). En remplaçant les
composantes des vitesses par leurs expressions en 𝜓 (Eq. (2.4)) et en substituant les équations
(3.13)-(3.15) respectivement dans les équations (2.23) et (2.24), nous obtenons le système
d’équations différentielles simplifiées suivant:

𝜕 𝜕𝜓 𝑛−1 𝜕𝜓 (3.17)
[(| | + ℋ 2 ) ] = −ℜ𝐶
𝜕𝑧 𝜕𝑧 𝜕𝑧
𝜕²𝜃 𝜕𝜓
=𝐶 (3.18)
𝜕𝑧² 𝜕𝑧
𝜕𝜃 (3.19)
{ 𝑧 = ±0.5 𝜓=
𝜕𝑧
+1=0

72
3.3.2 Nombre de Nusselt

Lorsque les parois actives de la cavité sont exposées à un flux de chaleur constant, le taux de
transfert de chaleur dans celle-ci s’exprime par le nombre de Nusselt 𝑁𝑢, tel que,

𝑞′ 1 1 (3.20)
𝑁𝑢(𝑥) = = =
𝑘∆𝑇′ ∆𝑇(𝑥) 𝑇(𝑥, − 1⁄2) − 𝑇(𝑥, 1⁄2)
𝐻′

3.3.3 Solutions de 𝝍, 𝑻 et 𝑵𝒖

La résolution analytique du système ci-dessus (3.17)-(3.18) avec les conditions aux frontières
(3.19) donne une solution générale, pour tout 𝑛, lorsque ℋ = 0. Toutefois, il n'y a pas de
solution analytique à l'exception de certains cas particuliers lorsque ℋ > 0.
 ℋ = 0: Fluide de type loi en puissance sans champ magnétique
o pour tout 𝑛 > 0
 ℋ > 0: Quelques cas avec champ magnétique
o 𝑛 = 1 Fluide newtonien
o 𝑛 = 2 fluide dilatant
o 𝑛 = 1/2 fluide pseudoplastique

a. Solutions pour 𝓗 = 𝟎 ; 𝒏 ≠ 𝟏

Dans cette section, nous considérons l'effet unique du flux thermique vertical sur le
mouvement convectif au sein de la couche poreuse saturée par un fluide de type loi en
puissance, les solutions sont données par :

1
(ℜ𝐶)𝑛 1 (3.21)
𝜓(𝑧) = − [|𝑧|𝑝 − 𝑝 ]
𝑝 2
1 𝐶𝑧 |𝑧|𝑝 1
𝜃(𝑧) = −(ℜ𝐶)𝑛 [ − 𝑝] − 𝑧 (3.22)
𝑝 𝑝+1 2
En substituant l’équation (3.22) dans l’équation (3.15), nous obtenons :

1 𝐶𝑧 |𝑧|𝑝 1
𝑇(𝑥, 𝑧) = 𝐶𝑥 − (ℜ𝐶)𝑛 [ − 𝑝] − 𝑧 (3.23)
𝑝 𝑝+1 2
où 𝑝 = (1 + 𝑛)⁄𝑛 et le coefficient 𝐶 est déterminé à partir de l’expression (3.16).

73
En remplaçant l’expression (3.23) dans l’équation (3.20), le nombre de Nusselt s’écrit sous la
forme suivante :

1 (3.24)
𝑁𝑢 = 1
1
1− 𝑝 ℜ𝑛 𝐶 𝑝
(𝑝 + 1)2
Ces solutions sont en accord avec les résultats rapportés par Ben Khelifa et al. [61] et qui a
étudié les cas particuliers 𝑛 = 1, 𝑛 = 1/3 et 𝑛 = 3.

Dans le cas d’un fluide newtonien pur 𝑛 = 1, on retrouve les résultats rapportés par Vasseur
et al. [77] en considérant les effets thermiques verticaux purs,

ℜ𝐶 2 1
𝜓(𝑧) = − [𝑧 − ] (3.25)
2 4
ℜ𝐶² 𝑧² 1
𝑇(𝑥, 𝑧) = 𝐶𝑥 − 𝑧[ − ] − 𝑧 (3.26)
2 3 4
1
𝐶=± √10(ℜ − 12) (3.27)
ℜ
1
𝑁𝑢 = (3.28)
1 10
+
6 ℜ
ℜ ici est le nombre de Rayleigh newtonien donné par

𝜌0 𝑔𝛽𝑇′ ∆𝑇 ′ 𝐾𝐻 ′
ℜ= (3.29)
𝜂𝜅
Dans la discussion qui suit, nous allons montrer que le critère de l’apparition de la convection
ℜ𝑐 peut-être prédit par l’approximation de l’écoulement parallèle tel que considérée dans la
présente étude. La Figure ‎3.4 représente un diagramme de bifurcation en termes de 𝜓𝑚𝑎𝑥 en
fonction du nombre de Rayleigh ℜ et de l'indice de loi de puissance 𝑛 pour ℋ = 0. Ces
résultats sont obtenus à partir de la solution analytique de 𝜓 (Eq. 3.21) pour différentes
valeurs de n. On constate qu’il existe différentes types de bifurcations selon la valeur de
l’indice de loi en puissance n. Les lignes en pointillées et continues dans ces graphiques sont
la prédiction de l'approximation de l’écoulement parallèle. Les symboles représentent les
solutions numériques des équations complètes et qui sont en bon accord avec les solutions
analytiques.

74
Les signes dans l’équation (3.27) explique l’apparition de deux branches à partir du point de
bifurcation. L’écoulement peut être dans le sens horaire ou antihoraire.

Le point de bifurcation correspondant à un fluide newtonien montre que la convection est

possible pour toute valeur de Rayleigh critique au-dessus de 12. Ce résultat est en accord avec
la prédiction de Nield [13], obtenue sur la base de la théorie de la stabilité linéaire.

Pour un fluide dilatant (𝑛 = 2), nous constatons que la convection est possible pour toute
valeur de ℜ ≥ 0. Ainsi, la théorie de l'écoulement parallèle non linéaire confirme la
prédiction de la stabilité linéaire montrant que les fluides dilatants sont thermiquement
instables pour tout flux de chaleur. Alors que pour les fluides pseudoplastiques (𝑛 = 0.5), on
montre le mouvement convectif est possible au-dessus d'un nombre de Rayleigh sous-critique
𝑠𝑢𝑝
ℜ𝑐𝑠𝑢𝑏 donc le critère ℜ𝑐 = ∞, prédit par la théorie de stabilité linéaire n’est pas réaliste. Le
diagramme de bifurcation présente deux branches stables pour lesquelles le seuil se produit
avec une convection à amplitude finie, et deux branches instables esquissées en fines lignes
discontinues.

n=0.5
2 n=1
n=2

max 0
branches instables

-1

-2

0 5 10 15 20 25 30 35 40

Figure ‎3.4 : Diagrammes de bifurcation en termes de 𝝍𝒎𝒂𝒙 en fonction du nombre de

Rayleigh 𝕽 et de l'indice de la loi de puissance 𝒏 pour 𝓗 = 𝟎.

75
 b. Solutions pour 𝓗 > 𝟎 ; 𝒏 = 𝟏

Pour le cas newtonien (𝑛 = 1) et en présence du champ magnétique ℋ > 0, les solutions sont
données par :

1 ℜ (3.30)
𝜓(𝑧) = 𝜓𝑚𝑎𝑥 (1 − 4𝑧 2 ) avec 𝜓𝑚𝑎𝑥 = √10 ( − 12)
8 1 + ℋ2
1
𝑇(𝑥, 𝑧) = 𝐶𝑥 + 𝐶𝜓𝑚𝑎𝑥 (3𝑧 − 4𝑧3 ) − 𝑧 (3.31)
3
1
𝑁𝑢 =
2 (3.32)
1 − 𝐶𝜓𝑚𝑎𝑥
3

En remplaçant les expressions (3.30) et (3.31) dans l’équation (3.16) et en intégrant, nous
pouvons montrer que :

1 + ℋ2 ℜ
𝐶= √10 ( − 12) (3.33)
ℜ 1 + ℋ2

A partir de l’équation (3.33) on peut déduire que la condition d’existance de mouvement de

fluide est obtenue pour : ℜ > 12(1 + ℋ 2 ) = ℜ𝑐 . Au-dessous de cette valeur de Rayleigh
critique, le fluide est au repos et nous avons l’état de conduction. Ce résultat montre que
l’approximation de l’écoulement parallèle est en accord avec le résultat de la stabilité linéaire
(Tableau ‎3.1).

c. Solutions pour 𝓗 > 𝟎 ; 𝒏 = 𝟐 ; 𝒛 < 𝟎

1 3 1 1 (3.34)
𝜓(𝑧) = [𝐴1 − (ℋ 4 − 4ℜ𝐶𝑧)2 ] − ℋ 2 (𝑧 + )
12ℜ𝐶 2 2
4
5
2
ℋ2 𝐶 2
10 (3.35)
𝑇(𝑥, 𝑧) = 𝐶𝑥 + 𝐴3 ((ℋ − 4ℜ𝐶𝑧) − ℋ ) − 𝑧 + ( 𝐴 2 − 1) 𝑧
4
1
𝑁𝑢 = (3.36)
1 − 𝐴7

76
 d. Solutions pour 𝓗𝟐 > 𝟎 ; 𝒏 = 𝟏/𝟐 ; 𝒛 < 𝟎

3
1 2 2
(1 − 4ℋ 2 ℜ𝐶|𝑧|)2 (3.37)
𝜓(𝑧) = − (ℋ ℜ𝐶𝑧 + − 𝑧 − 𝐴4 )
2ℋ 4 6ℋ 2 ℜ𝐶

ℜ𝐶2 3
𝐶 2 2 2
5 (3.38)
𝑇(𝑥, 𝑧) = 𝐶𝑥− 2𝑧 +
( )
4 𝑧 − 𝐴 (1 − 4ℋ ℜ𝐶𝑧) + 𝐴6 − 1 𝑧 + 𝐴5
6ℋ 4ℋ 5
1
𝑁𝑢 = (3.39)
1 − 𝐴8
avec

−3ℜ𝐶ℋ 2 + (ℋ 4 + 2ℜ𝐶)3⁄2
𝐴1 = 2ℜ𝐶 + ℋ 4 𝐴2 =
12ℜ
1 1 (1 + 2ℋ 2 ℜ𝐶)3/2 ℋ 2 ℜ𝐶
𝐴3 = 𝐴4 = + +
120ℜ2 𝐶 2 6ℋ 2 ℜ𝐶 4
3
1 6ℜ𝐶ℋ 2 −2(1 + 2ℋ 2 ℜ𝐶)2 + 3ℋ 4 ℜ2 𝐶 2
𝐴5 = 𝐴6 =
120ℜ2 𝐶ℋ 8 24ℜℋ 6

𝐴7 = 𝐴3 [(ℋ 4 − 2ℜ𝐶)5⁄2 −(ℋ 4 + 2ℜ𝐶)5⁄2 + 5𝑅𝐶((ℋ 4 + 2ℜ𝐶)3⁄2 + (ℋ 4 − 2ℜ𝐶)3⁄2 )]

𝐴8 = 𝐴5 [(1 − 2ℋ 2 ℜ𝐶)5⁄2 + (1 + 2ℋ 2 ℜ𝐶)5⁄2 + 10ℋ 6 ℜ3 𝐶 3 − 2
+ 5ℋ 2 ℜ𝐶((1 − 2ℋ 2 ℜ𝐶)3⁄2 − (1 + 2ℋ 2 ℜ𝐶)3⁄2 )]

Les solutions (3.34)-(3.39) sont trouvées avec la condition 𝑧 = −0.5, nous pouvons trouver
d'autres solutions lorsque 𝑧 > 0, (𝑧 = 0.5).

3.3.4 Évaluation de la constante C

Le coefficient C est défini comme étant le gradient de température horizontale. Il

dépend du nombre de Rayleigh ℜ, de l’indice de comportement 𝑛 et du nombre de Hartmann
ℋ. Pour certaines valeurs remarquables de 𝑛, déterminer la constante 𝐶 pour un (𝑛, ℜ, ℋ 2 )
donné est un problème algébrique sévère. Pour cette raison, le système (3.17-18) doit être
modifié de manière appropriée. Intégrer l'équation (3.18) et prenant en compte l'équation
(3.19) donne:

77
 𝜕 𝜕𝜓 𝑛−1 𝜕𝜓 (3.40)
[(| | + ℋ 2 ) ] = −ℜ𝐶
𝜕𝑧 𝜕𝑧 𝜕𝑧
0.5 0.5
𝐶=∫ 𝜓 𝑑𝑧⁄(1 + ∫ 𝜓2 𝑑𝑧) (3.41)
−0.5 −0.5
𝜕𝜃 (3.42)
𝑧 = ±0.5 𝜓= +1=0
{ 𝜕𝑧
L’équation (3.41) est résolue à l'aide d'un schéma numérique explicite. Pour un ℜ donné, avec
une valeur de conjecture de 𝐶, on résout l'équation (3.40) en utilisant la méthode Müller pour
obtenir 𝜓, puis on intègre (3.41) en utilisant la méthode Simpson pour mettre à jour 𝐶. Les
itérations continuent jusqu'à la convergence du schéma explicite.

3.4 Simulations numériques directes du problème non-linéaire

Dans cette section, nous allons présenter quelques simulations numériques directes des
équations non-linéaires. De manière générale, les approches numériques permettent de
résoudre de manière approchée des problèmes pour lesquels la solution analytique n’existe
pas. En effet, on a vu dans la partie précédente que la résolution analytique des équations aux
dérivées partielles non-linéaires est impossible dans certain cas. Deux approximations
s’imposent : l’approximation spatiale des opérateurs aux dérivées partielles et celle, si
nécessaire, des opérateurs temporels. Cependant, il est possible de résoudre ces équations
numériquement après une discrétisation du domaine on obtient un système d'équations
algébriques linéaires. Ensuite le système est résolu par des méthodes directes ou itératives.

La méthode semblant la mieux appropriée à nos objectifs, parmi la variété des méthodes
existantes, est la méthode des différences finies.

3.4.1 Validation du code numérique

Pour valider le code numérique élaboré, il est indispensable de comparer les résultats
obtenus à ceux existant dans la littérature. Cependant, en absence de travaux antérieurs
publiés traitant le cas Darcy-Bénard avec un fluide non-newtonien de type loi en puissance en
présence d’un champ magnétique appliqué. Nous proposons de le confronter dans des
situations traitant les effets séparés suivants:

- Fluide non-newtonien en milieu de Darcy pour le cas de la convection naturelle de

Darcy-Bénard (DB) avec condition thermique de type Neumann.

78
 - Fluide newtonien en milieu de Darcy pour le cas de la convection naturelle de Darcy-
Bénard (DB) en présence d’un champ magnétique et avec condition thermique de type
Neumann.
Nous présentons dans ce qui suit les différentes confrontations effectuées en termes du
nombre de Nusselt 𝑁𝑢 et de l’intensité d’écoulement maximale 𝜓𝑚𝑎𝑥 .

 Cas de la convection naturelle d’un fluide non-newtonien saturant un milieu

poreux soumis à une condition thermique de type Neumann et en absence d’un
champ magnétique :
Nous avons confronté avec succès nos résultats à ceux obtenus par Ben Khelifa et al [61] dans
le cas d’une cavité rectangulaire poreuse saturée par un fluide non-newtonien et soumise à un
flux de chaleur constant pour A=4 et ℜ=50.

Ben Khelifa et al [61] Présente travail

𝒏 𝝍𝒎𝒂𝒙 𝑵𝒖 𝝍𝒎𝒂𝒙 𝑵𝒖
0.6 4.29 4.78 4.29 4.78
1.0 2.44 2.73 2.46 2.73
1.6 1.41 1.64 1.44 1.64

Tableau ‎3.2 : Comparaison des résultats numériques obtenus par la présente étude avec ceux
obtenus par Ben Khelifa et al. [61] en milieu poreux avec A=4 et ℜ=50.

 Cas de la convection naturelle d’un fluide newtonien saturant un milieu poreux

soumis à une condition thermique de type Neumann et en présence d’un champ
magnétique :
Ce code permet aussi de reproduire avec une excellente concordance les résultats d’une
solution analytique développée dans la section précédente (§‎3.3) pour le cas d’une cavité
rectangulaire poreuse saturée par un fluide newtonien et soumise à un flux de chaleur et
champ magnétique constant.

Solution numérique Solution analytique

𝝍𝒎𝒂𝒙 1.46 1.43
𝑵𝒖 1.76 1.76
Tableau ‎3.3 : Comparaison entre les résultats numériques et les résultats analytiques dans le
cas d’un fluide newtonien, pour ℜ =50, ℋ=1 et A=4.

79
Ces comparaisons montrent que notre code donne des résultats très satisfaisants. En effet,
l’écart relatif maximal est de 2.08% pour la fonction de courant maximal
𝜓𝑚𝑎𝑥 et de 0% pour le nombre de Nusselt 𝑁𝑢.

3.5 Analyse des résultats non linéaires

Nous traitons dans cette partie, l’influence de l’indice de la loi en puissance, du nombre
de Rayleigh, qui caractérise les mouvements thermiques au sein du fluide, et du nombre de
Hartmann, qui caractérise les mouvements magnétiques au sein du fluide, sur le
déclenchement de la convection et faire la comparaison avec la théorie de stabilité linéaire et
la solution analytique non linéaire. Et enfin nous allons clore cette partie en étudiant le
comportement des fluides de type loi en puissance pour des faibles et des grandes valeurs du
nombre de Hartmann.

3.5.1 Cas d’une couche poreuse saturée par un fluide newtonien

L’influence du nombre de Rayleigh ℜ et du nombre de Hartmann ℋ sur le gradient de

température horizontale 𝐶 est représentée sur la Figure ‎3.5.

0.5

0.4

C 0.3

H 2= 0
0.2
H 2= 1
H 2= 2
0.1
0 40 80 120 160 200

Figure ‎3.5 : Gradient de température horizontal 𝐶(ℜ, ℋ) pour les fluides newtoniens (𝑛 = 1)
soumis à un champ magnétique vertical.

80
Comme l’indiquent les équations (3.25) et (3.26), la circulation du fluide dans la cavité poreuse
dépend du gradient longitudinal de la température C. Comme déjà discuté dans l’analyse de
stabilité linéaire (Tableau ‎3.1) le champ magnétique a un effet inhibiteur sur la convection et
l’apparition de la convection est retardé d'un facteur de (1 + ℋ²). La Figure ‎3.5 montre que
le gradient de température horizontal C existe pour les nombres de Rayleigh supérieurs à
12(1 + ℋ²) et augmente avec le nombre de Hartmann. Cela signifie que le champ
magnétique, malgré le retard du mouvement du fluide et la réduction de l'intensité du courant,
contribue à ce que les zones de fluide les plus chaudes soient concentrées dans l'extrémité du
canal, comme le montre la Figure ‎3.6. Ces résultats sont en accord avec ceux obtenus par
Alloui et al. [43].

Figure ‎3.6 : Isothermes prédites par la solution numérique des équations complètes pour
𝕽 = 𝟕𝟎, 𝒏 = 𝟏 et pour différentes valeurs du nombre de Hartmann 𝓗.

3.5.2 Cas d’une couche poreuse saturée par un fluide non-newtonien

 En absence du champ magnétique 𝓗 = 𝟎

L’influence du nombre de Rayleigh ℜ et de l’indice de loi de puissance 𝑛 sur le gradient

de température horizontal 𝐶 pour ℋ = 0 est représentée sur la Figure ‎3.7. Les solutions

81
numériques sont représentées par des symboles et les solutions analytiques sont représentées
par des lignes continues. Nous observons que les solutions analytiques de l’intensité de
courant 𝜓(𝑧) et la température 𝑇(𝑥, 𝑧) dépendent du gradient de température horizontal 𝐶 (éq
(3.34)-(3.39)). Il n'est pas facile de trouver une expression explicite de 𝐶 même pour des
valeurs particulières de 𝑛, et 𝐶 est calculé alors en utilisant numériquement l'algorithme décrit
au §‎3.3.4.

Mathématiquement, l'existence de la fonction 𝐶(ℜ, 𝑛) donne l'apparition d'un mouvement

possible. Physiquement, les valeurs de 𝐶 renseignent sur la distorsion des isothermes.

0.5

0.4

0.3
C
0.2
n = 0.5
n=1
0.1 n=2

0.0
0 20 40 60 80 100 120


Figure ‎3.7 : Gradient de température horizontal 𝐶 en fonction du nombre de Rayleigh ℜ et
indice de loi de puissance 𝑛 pour ℋ = 0.

Il est constaté que pour 𝑛 = 1, la convection existe pour chaque nombre de Rayleigh
supérieur à 12. Pour 𝑛 > 1, la constante C est strictement positive comme présenté sur la
Figure ‎3.7, et par conséquent la convection existe pour chaque nombre de Rayleigh positif.
Pour 𝑛 < 1 et contrairement à l'analyse de stabilité linéaire qui prédit une surstabilité ℜ𝑐 = ∞
pour les fluides pseudoplastiques, les solutions numériques ainsi que notre modèle basé sur la

82
théorie de l’écoulement parallèle non linéaire donnent un autre critère pour le début de la
convection. Nous remarquons à partir de la Figure ‎3.7, que pour 𝑛 = 0.5 nous avons des
valeurs réelles de C à partir de ℜ𝑐 ≃ 9, donc il existe un mouvement convectif sous-critique
plus réaliste que la surstabilité.

D’autre part, l’effet de l’indice 𝑛 sur le gradient de température horizontal apparaît dans la
Figure ‎3.7. Il est observé que les fluides dilatants présentent un gradient de température
horizontal supérieur au cas newtonien pour des nombres de Rayleigh ℜ > 75 et que les
fluides pseudoplastiques présentent un gradient de température horizontal supérieur au cas
newtonien, pour des nombres de Rayleigh ℜ > 100. Dans un sens, cela signifie que les
fluides de type loi de puissance présentent des régions plus chaudes à la fin de la couche
lorsque ℜ est suffisamment élevé.

Les diagrammes de bifurcation correspondant aux fluides dilatants en termes de 𝜓𝑚𝑎𝑥 et Nu
en fonction de ℜ et n sont illustrés pas les Figure ‎3.8-(a) et (b). À titre de référence, le cas
d’un fluide newtonien ( 𝑛 = 1) est aussi représenté dans les graphiques. Ces résultats sont
obtenus grâce à l’approximation de l’écoulement parallèle et la simulation numérique qui
présentent un très bon accord. Ils indiquent que le transfert de chaleur et l’intensité de
l’écoulement dépendent considérablement du nombre de Rayleigh et de l'effet de l’indice 𝑛.
Pour ce type des fluides et comme a été vu précédemment, la convection est possible pour
toute valeur du nombre de Rayleigh supérieure à zéro quel que soit la valeur de n 1.

L'effet de l'indice de puissance 𝑛 sur le transfert de chaleur par convection est très différent.
Pour les fluides dilatants (𝑛 = 2 et 𝑛 = 3), l’intensité de l’écoulement (𝜓𝑚𝑎𝑥 ) et le taux de
transfert de chaleur (Nu), sont renforcés lorsque la valeur de n est de plus en plus élevée. Nos
résultats sont parfaitement en accord avec l’ensemble des travaux existants dans la littérature.

83
 (a)
3.0

2.5

2.0

max 1.5

1.0
n=1
n=2
0.5 n=3

0.0
0 10 20 30 40 50 60


(b)
1.75

1.50

Nu

1.25 n=1
n=2
n=3
1.00
0 10 20 30 40 50 60


Figure ‎3.8 : Evolution du nombre de Nusselt 𝑁𝑢 et l’intensité de l’écoulement maximale
𝜓𝑚𝑎𝑥 en fonction du nombre de Rayleigh ℜ et de l'indice de loi de puissance 𝑛 pour ℋ = 0.

84
  En présence d’un champ magnétique 𝓗 ≠ 𝟎

Nous envisageons maintenant un champ magnétique parallèle à la gravité régnant dans

un fluide de type loi en puissance.

L’effet de l’indice de loi en puissance n sur le nombre de Rayleigh critique ℜ𝑐 pour

l’apparition des mouvements convectifs est présenté dans la Figure ‎3.9. Dans le cas des
fluides dilatants, la présence du champ magnétique augmente la valeur de Rayleigh
𝑠𝑢𝑝
supercritique de zéro à ℜ𝑐 = 12ℋ 2 , ce qui est un résultat important. Nous avons déjà vu
dans la Figure ‎3.3 que 12ℋ 2 est une limite asymptotique pour toute valeur de 𝑛.

sup
c
= sup = 12 H 2
c
24

c
16

sub
c

8
sub
c

0.0 0.5 1.0 1.5 2.0

Figure ‎3.9 : ℜ𝑐 pour le début du mouvement convectif en fonction de l'indice de loi de

puissance 𝑛 pour ℋ 2 = 2. La solution présentée en traits pointillés est le résultat de la
stabilité linéaire.

Dans le cas des fluides pseudoplastiques (n 1), la théorie de stabilité linéaire (Tableau ‎3.1)
indique que le nombre de Rayleigh supercritique marquant la naissance des mouvements
𝑠𝑢𝑝
convectifs est ℜ𝑐 = ∞,présenté en lignes en pointillés, qui est la même limite supercritique
en absence du champ magnétique. Cette valeur montre que le système est
inconditionnellement stable pour des perturbations infinitésimales. Cependant, comme la
montre la Figure ‎3.9, le présent modèle non linéaire prédit que l’apparition de la convection se

85
fait par une bifurcation souscritique présentée en cercles creux, ce sont les résultats de
simulations numériques directes. Les lignes pleines sont des interpolations reliant entre les
cercles. En augmentant n, le nombre de Rayleigh souscritique ℜ𝑐𝑠𝑢𝑏 continue à augmenter
jusqu’à une valeur maximale ℜ𝑐𝑠𝑢𝑏 = 24 à 𝑛 ≃ 1 puis la ligne souscritique présente une
discontinuité à 𝑛 = 1 et ℜ𝑐𝑠𝑢𝑏 fait un saut soudain vers 6 puis continue à décroitre en
augmentant n. Cette discontinuité est due à la discontinuité de la viscosité effective donnée
par 𝜂(𝛾̇ )𝑛−1 . Pour le fluide pseudoplastique, le taux de cisaillement 𝛾̇ sera très faible à
l'échelle du pore et par conséquent les fluides seront extrêmement visqueux, ce qui rend ℜ𝑐
tendant vers l'infini. Cependant, les fluides dilatants seront non visqueux, donc ℜ𝑐 tendra
vers zéro.

 Le comportement asymptotique

Tel que discuté dans l’analyse strictement linéaire c’est à dire la partie de la stabilité
linéaire de la convection de Darcy-Bénard (DB), ℜ𝑐 a un comportement asymptotique pour
les grandes valeurs de ℋ 2 . Dans cette discussion, nous essayons de répondre à la question
suivante: le comportement asymptotique est-il toujours valide si nous considérons les termes
non linéaires?

En considérant l’approximation de l’écoulement parallèle, nous pouvons écrire l'équation

(3.43) déduite de (3.17),

𝜕 𝜕𝜓 𝑛−1 𝜕𝜓 (3.43)
(| | )
𝜕𝑧 𝜕𝑧 𝜕𝑧 𝜕²𝜓 ℜ
+ = − 2𝐶
ℋ2 𝜕𝑧² ℋ
Pour un ℋ 2 élevé et comme on peut le montrer à partir de l’équation (3.43), le système
puisse dépendre uniquement de ℜ⁄ℋ 2 si le terme de loi de puissance disparaît. Nous
effectuons également des investigations numériques afin d'approfondir la question. Les
résultats obtenus sont présentés sur la Figure ‎3.10 en termes de 𝜓𝑚𝑎𝑥 et 𝑁𝑢 en fonction de
ℜ⁄ℋ 2 , pour un champ magnétique faible (ℋ 2 = 2) sur la Figure ‎3.10-(a) et un fort champ
magnétique (ℋ 2 = 103 ) dans la Figure ‎3.10-(b). Les lignes continues représentent la solution
analytique tandis que les symboles sont les résultats numériques des équations complètes
régissant le problème.

86
 2
(a) H = 2
3.0 3.0

2.5
n=0.5 max 2.5
n=1 Nu
2.0 n=2 2.0

max 1.5 1.5 Nu

1.0 1.0

0.5 0.5

0.0 0.0
10 20 30 40 50 60
/H 2

2 3
(b) H = 10
4.5 4.5
4.0 4.0
n=0.5 max
3.5 3.5
n=1 Nu
3.0 n=2 3.0
2.5 2.5
max Nu
2.0 2.0
1.5 1.5
1.0 1.0
0.5 0.5
0.0 0.0
0 20 40 60 80 100
/H 2

Figure ‎3.10 : 𝝍𝒎𝒂𝒙 et 𝑵𝒖 en fonction de 𝕽⁄𝓗𝟐 et 𝒏 pour: (a) 𝓗𝟐 = 𝟐 et (b) 𝓗𝟐 = 𝟏𝟎𝟑 .

87
Lorsque la valeur du nombre de Hartmann ℋ 2 = 2, c’est à dire pour un faible champ
magnétique, la Figure ‎3.10-(a) montre une différence nette entre les fluides newtoniens,
pseudoplastiques et dilatants. Dans un sens, cela signifie que la dissipation par viscosité
prédomine à une telle plage du nombre de Hartmann. Cependant, la Figure ‎3.10-(b) montre
une dépendance parfaite de 𝜓𝑚𝑎𝑥 et 𝑁𝑢 sur ℜ⁄ℋ 2 aux nombres de Hartmann élevés ℋ 2 =
103 . De telle sorte que la résistance des effets visqueux devient faible en présence d'un fort
champ magnétique. En effet, nous voyons que les trois courbes (𝑛 = 0.5, 𝑛 = 1 𝑒𝑡 𝑛 = 2) se
trouvent sur la même ligne, les trois types de fluides s'écoulent exactement de la même
manière et on ne peut plus décrire le fluide comme newtonien ou non-newtonien, puisque la
viscosité n'a plus d'importance, elle a été inhibée par l'énergie magnétique. Les résultats
numériques confirment que les déductions de l'analyse de stabilité linéaire peuvent être
étendues au domaine non linéaire.

3.6 Conclusion
Dans ce chapitre, la convection au sein d’une une couche poreuse horizontale d’un
fluide non-newtonien dans le cas de condition thermique de type Neumann a été étudiée. Les
parois horizontales de la cavité sont soumises à un flux de chaleur uniforme et celles
verticales sont adiabatiques. Un champ magnétique constant est appliqué parallèlement à la
gravité 𝑔. Un modèle de de type loi en puissance a été utilisé pour caractériser le
comportement non-newtonien du fluide. L’effet de l’indice de la loi en puissance, du nombre
de Hartmann, du nombre de Rayleigh et du nombre de Péclet sur les écoulements et les
transferts de chaleur ont été présentés pour une large gamme des paramètres régissant le
problème.

L’objective principale était de déterminer des solutions explicites exactes basées sur
l’approximation de l’écoulement parallèle. Le problème semble être plus difficile que prévu,
et seules quelques solutions particulières ont été trouvées (𝑛 = 0.5, 𝑛 = 1 et 𝑛 = 2). En
conséquence, l'analyse de stabilité linéaire et la résolution numérique par différences finies
ont été utilisées en plus.

Les résultats montrent que le champ magnétique a un effet inhibiteur. En effet, pour les
fluides newtoniens, la présence du champ magnétique retarde l’apparition de la convection
d'un facteur de (1 + ℋ 2 ), ce qui est en accord avec les résultats donnés par Chandrasekhar

88
[32]. Pour les fluides dilatants, la présence du champ magnétique augmente la valeur
𝑠𝑢𝑝
classique du nombre de Rayleigh pour l’apparition de la convection de zéro à ℜ𝐶 = 12ℋ 2 .
Pour les fluides pseudoplastiques, selon la théorie de stabilité linéaire, le système est
inconditionnellement stable pour des perturbations infinitésimales. Cependant, la théorie non
linéaire a prédit que l’apparition de la convection se produit à un certain nombre de Rayleigh
𝑠𝑢𝑝
souscritique ℜ𝐶 , qui dépend des valeurs de 𝑛 et ℋ, ce qui est en accord avec les résultats
des travaux de Ben Khelifa et al. [61] et Nield [13].

En présence d'un fort champ magnétique, un résultat important a été trouvé concernant le
caractère non visqueux des fluides pseudoplastiques et dilatants, car la dissipation d'énergie
par effet Joule domine la dissipation de l'énergie par la contrainte de cisaillement.

89
Chapitre 4 Convection MHD dans une
cavité poreuse soumise à une
condition de température constante

4.1 Introduction
L’objectif de ce chapitre consiste à étudier la stabilité linéaire et la simulation
numérique de la convection au sein d’un fluide non-newtonien saturant une cavité poreuse
horizontale avec des conditions thermiques de type Dirichlet. Nous allons étudier le cas de
condition de température imposée sur les parois horizontales (parallèle au plan (𝑜𝑥’𝑦’)).
Expérimentalement, ce type d’écoulement est difficilement réalisable. Toutefois, il est
théoriquement intéressant puisqu’il permet de trouver une solution analytique. Pour cela, nous
présentons l’analyse de la stabilité linéaire à travers laquelle nous déduisons l’expression du
nombre de Rayleigh critique qui permettra de prédire l’apparition de la convection au sein du
fluide.

Figure ‎4.1 : Représentation schématique de la configuration géométrique en 3D.

90
 La configuration géométrique étudiée ainsi que les conditions aux frontières sont
montrées par la Figure ‎4.1. L’objectif principal de cette étude est de prédire l’influence du
champ magnétique et de l’indice de loi en puissance n sur l’apparition et l’intensité des
mouvements convectifs et le transfert de chaleur résultant.

4.2 Stabilité linéaire

̃ (𝑧) et 𝜃̃ (𝑧) qui
Dans le cas de condition de température imposée, les fonctions d’amplitude 𝑊
vérifient les conditions aux limites (2.44) sont :

𝑁 𝑁
̃ (𝑧) = ∑ 𝑊𝑙 cos(𝑙𝜋𝑧)
𝑊 𝑒𝑡 𝜃̃ (𝑧) = ∑ 𝜃𝑙 cos(𝑙𝜋𝑧) (4.1)
𝑙=1 𝑙=1

Comme discuté dans le chapitre précèdent, pour déterminer le mode critique (𝑘𝑐 , ℜ𝑐 )
marquant la naissance des mouvements convectifs, on doit calculer le nombre d’onde qui
minimise l’équation de dispersion (2.50) pour 𝑙 = 1.
En introduisant les développements (4.1) dans le système (2.46). A l’ordre 𝑙 = 1, le système
s’écrit sous forme matricielle comme suit :

−(𝑆 2 𝜋 2 + 𝑘 2 ) 𝑘²ℜ𝒫1−𝑛 𝑊1 0
[ ] . [ ] = [ ] (4.2)
−1 𝑘² + 𝑖Ω + 𝜋 2 𝜃1 0
Le système (4.2) admet une solution non triviale si et seulement si 𝑑𝑒𝑡(𝑘, ω) = 0. Donc
l’équation de dispersion complexe (2.50) s’écrit sous la forme suivante :

𝑑𝑒𝑡(𝑘, ω) = −(𝑆 2 𝜋 2 + 𝑘 2 )(𝑘 2 + 𝑖Ω + 𝜋 2 ) + 𝑘²ℜ𝒫1−𝑛 = 0 (4.3)

Donc

(𝑆 2 𝜋 2 + 𝑘 2 )(𝑘 2 + 𝑖Ω + 𝜋 2 )
ℜ = 𝒫 n−1 (4.4)
𝑘²
On notera toutefois que le nombre de Rayleigh ℜ(ℋ, 𝒫, 𝑘, 𝜔, 𝑛) est un paramètre de contrôle
qui doit être réel (on impose 𝐼𝑚𝑔(ℜ) = 0). La valeur du nombre d'onde critique 𝑘𝑐 est
obtenue en imposant 𝜕ℜ⁄𝜕𝑘 = 0. Ce qui donne

𝑘𝑐 = 𝜋√𝑆 (4.5)

91
On obtient ainsi le nombre ℜ𝑐 qui marque la naissance de la convection stationnaire (dite
aussi supercritique) :

ℜ𝑐 = 𝒫 𝑛−1 𝜋²(𝑆 + 1)² (4.6)

Nous notons à partir de ces résultats, de l’expression de 𝑆² (Eq. (2.42)) et de la Figure ‎4.2 que
𝑆 est une fonction monotone croissante de 𝛼 et que les perturbations les plus instables sont
des rouleaux transversaux 𝑅∥ (𝛼 = 0), pour les fluides pseudoplastiques (𝑛 < 1). Par contre
pour les fluides dilatants (𝑛 > 1), 𝑆 présente une faible dépendance de 𝛼. La Figure ‎4.2
montre la cohabitation des rouleaux transversaux 𝑅∥ et longitudinaux 𝑅⊥ . Pour les fluides
newtoniens (𝑛 = 1), nous avons 𝑆² = 1 + ℋ 2 donc l'analyse de stabilité est indépendante de
l'angle d'inclinaison 𝛼, toutes les orientations des rouleaux obliques sont équivalents. Nous
obtenons alors,

𝑆² = 𝑛 + 𝒫1−𝑛 ℋ 2 𝑝𝑜𝑢𝑟 (𝑛 < 1) → 𝑅∥

{ 𝑆² = 1 + 𝒫1−𝑛 ℋ 2 𝑝𝑜𝑢𝑟 (𝑛 < 1) → 𝑅∥ + 𝑅⊥ (4.7)
𝑆² = 1 + ℋ 2 𝑝𝑜𝑢𝑟 (𝑛 = 1)
La connaissance de la structure d’écoulement obtenu va faciliter la démarche
d’expérimentation numérique.

(a) n=0.5

1.2 90°
0°
1.1

1.0

0.9
S
0.8

0.7

0.6

0.5
0 1 2 3 4 5 6
H

92
 (b) n=1.5
1200
90°
1000 0°

800

600
S
400

200

0 1 2 3 4 5 6
H

Figure ‎4.2 : Naissance des structures convectives sous forme de rouleaux longitudinaux 𝑹⊥
ou transversaux 𝑹∥ selon le type de fluide : (a) pseudoplastique 𝒏 = 𝟎. 𝟓 et (b) dilatant
𝒏 = 𝟏. 𝟓.

4.2.1 Convection de Darcy-Bénard (DB)

Tel que discuté pour le cas de la condition de type Neumann, les critères (𝑘𝑐 , ℜ𝑐 ) pour
l’apparition des mouvements convectifs de Darcy-Bénard (DB) peuvent être prédits en
prenant 𝒫 → 0 et 𝜔 = 0 (stabilité marginale).

En conséquence des équations (4.5) et (4.6), nous obtenons les expressions du nombre de
Rayleigh critique ℜ𝑐 et du nombre d’onde critique 𝑘𝑐 pour le problème de DB regroupés dans
le Tableau ‎4.1.

Température constante
𝓟→ 𝟎
ℜ𝑐 𝑘𝑐
2
𝒏=𝟏 𝜋² (√1 + ℋ² + 1) 𝜋√1 + ℋ²

𝒏<𝟏 ∞ 𝜋𝑛1⁄4
𝒏>𝟏 ℋ². 𝜋² ∞
Tableau ‎4.1 : Les expressions du nombre de Rayleigh critique ℜ𝑐 et du nombre d’onde
critique 𝑘𝑐 pour le problème de DB (Condition de Dirichlet).
93
Dans le problème de DB avec les conditions thermiques de type Dirichlet, les résultats
présentés dans le Tableau ‎4.1 montrent que pour le cas des fluides pseudoplastiques (𝑛 < 1),
la valeur du nombre de Rayleigh critique pour l’apparition de la convection, indépendamment
de ℋ, se produit à ℜ𝑐 = ∞. Ces résultats sont similaires avec les résultats rapportés dans le
Tableau ‎3.1 dans le cas des conditions thermiques de type Neumann. Comme il a été déjà
mentionné dans le chapitre 3, le critère ℜ𝑐 = ∞, prédit par la théorie de stabilité linéaire n’est
pas réaliste.

D'autre part, nous déduisons que l’apparition de la convection dépend du nombre de

Hartmann ℋ pour le cas des fluides newtoniens (𝑛 = 1) ou fluides dilatants (𝑛 > 1). En
effet, pour 𝑛 = 1 et ℋ = 0 on obtient ℜ𝑐 = 4𝜋², le même résultat a été rapporté récemment
par Nield [57] en considérant la même configuration avec des conditions aux limites
thermiques de type Dirichlet mais en absence de champ magnétique. Alors que pour ℋ ≠ 0,
l'apparition de l'instabilité s’effectue pour ℜ𝑐 > 4𝜋², elle dépend du nombre de Hartmann
ℋ.Pour 𝑛 > 1 et ℋ = 0, on obtient ℜ𝑐 = 0 le système est inconditionnellement instable.
Alors que pour ℋ ≠ 0, l'apparition de l'instabilité ne s’effectue qu’à partir de ℜ𝑐 = ℋ². 𝜋².

4.2.2 Analyse des résultats non linéaires

En examinant les expressions du nombre de Rayleigh critique ℜ𝑐 (Tableau ‎4.1), on voit

facilement que l’apparition des mouvements convectifs dépend de plusieurs paramètres de
contrôle tels que le nombre de Hartmann ℋ, l’indice de la loi en puissance 𝑛 et le nombre de
Péclet 𝒫. Les résultats prédits par l’analyse de la stabilité linéaire seront présentés dans cette
partie en fonction de ces paramètres.

Dans la plus part des résultats, les fluides dilatants sont présentés pour n=1.6, les fluides
pseudoplastiques pour 𝑛 = 0.6 et évidement 𝑛 = 1 pour les fluides newtoniens. On note que
lorsque le nombre de Péclet est très faible 𝒫 → 0, nous obtenons une bonne approximation de
la convection naturelle.

 Effet du nombre du Péclet et du nombre d’Hartmann sur

l’apparition de la convection

94
 Dans cette partie de l'étude, les effets dus à la modification des fluides de type loi en
puissance, à l'augmentation de la valeur du nombre de Hartmann et le nombre de Péclet sur
l’apparition de la convection sont représentés.

En vue d'apprécier le rôle joué par le nombre de Péclet 𝒫, nous avons tracé l’évolution de ℜ𝑐
en fonction de 𝑛 pour des grandes et des faibles valeurs de 𝒫 sur la Figure ‎4.3. Un bref aperçu
de cette figure montre que le comportement qualitatif des caractéristiques linéaires de la
stabilité pour des conditions aux limites thermiques de type Dirichlet est similaire à celui
observé pour des conditions aux limites thermiques de type Neumann largement discuté dans
le chapitre précédent. La variation de comportement de ℜ𝑐 semble très similaire entre les
Figure ‎4.3-(a) et (c) lorsque 𝒫 ≤ 1 et entre les Figure ‎4.3-(b) et (d) lorsque 𝒫 ≥ 1. On
remarque sur la Figure ‎4.3 qu'un débit filtrant très faible (𝒫 = 0.01) est suffisant pour voir un
équilibre instable. Un résultat similaire a été rapporté récemment par Barletta et Nield [58] en
considérant la même configuration avec les mêmes conditions aux limites thermiques de type
Dirichlet mais en absence d’un champ magnétique.

La dépendance du nombre de Rayleigh critique vis-à-vis de l’indice de loi en puissance 𝑛 est

également illustrée sur la Figure ‎4.3. C’est clairement indiqué que l’évolution de Rayleigh
critique ℜ𝑐 en fonction de 𝑛, pour un 𝒫 donné, est monotone. En tenant compte du système
d’équations (4.7), nous pouvons dire que nous avons des rouleaux transversaux ou
longitudinaux lorsque 𝑛 > 1 et des rouleaux transversaux lorsque 𝑛 < 1.

Le débit filtrant 𝒫 avec condition de glissement donne une simple transformation galiléenne
dans le cas fluide d’un newtonien. En effet, pour 𝑛 = 1, ℜ𝑐 est égal au seuil de convection

naturelle pour toute valeur de Péclet. Et pour 𝒫 = 1, nous obtenons ℜ𝑐 = 𝜋² (√1 + ℋ² +

2
1) pour chaque 𝑛 ≥ 1. Cependant, pour les fluides de loi de puissance, le petit nombre de

Péclet rompt la transformation galiléenne en raison de la dépendance de ℜ𝑐 en 𝒫 𝑛−1 (Eq.

(4.6)). Puisque la viscosité dépend à la fois de la vitesse et de l’indice de loi de puissance 𝑛.
Un petit nombre de Péclet rend le liquide plus visqueux pour 𝑛 < 1, et moins visqueux pour
𝑛 > 1. Les mécanismes de dissipations visqueuses varient avec l'intensité du débit filtrant,
d'où nous observons la rupture de la transformation galiléenne. C'est un aspect très important
des fluides de type loi en puissance.

95
L’effet du nombre de Hartmann ℋ sur le nombre de Rayleigh critique ℜ𝑐 pour l’apparition
des mouvements convectifs est également présenté sur la Figure ‎4.3. On observe un grand
décalage sur l’axe de ℜ𝑐 lorsque ℋ² = 25. Ce décalage est dû à l’application d’un champ
magnétique qui retarde le seuil de convection d'un facteur de 𝜋 2 ℋ 2 c’est la limite
asymptotique pour toute valeur de 𝑛 (Figure ‎4.5).

80
(a) 0.01 H 2= 0 (b) 10000 H 2= 0 80
0.1 100
0.2 10
60 60
2

0.4 P=1
c 40 P=1 40

0.6 0.6
P=1 0.4
20 0.2 P=1 20

0.1 2
0.01 10 100
0 10000 0

600
(c)
0.01
H 2= 25 (d) n
10000 H 2= 25
600

0.1 100

500 500
0.2 10

0.4 2
c 400 0.6 400
P=1 P=1 P=1
P=1 0.6
0.4 2
300 0.2 10 300
0.1 100
0.01 10000
200 200
0.0 0.5 1.0 1.5 2.0
0.0 0.5 1.0 1.5 2.0

n n

Figure ‎4.3 : Evolution du nombre de Rayleigh critique ℜ𝑐 en fonction de 𝑛 pour différentes

valeurs de Péclet 𝒫 et pour deux valeurs du nombre d’Hartmann ℋ : (a) ≤ 1 ℋ 2 = 0, (b)
𝒫 ≥ 1 ℋ 2 = 25, (c) 𝒫 ≤ 1 et ℋ 2 = 0 et (d) 𝒫 ≥ 1 et ℋ 2 = 25 (d).

La dépendance du nombre d'onde critique 𝑘𝑐 et de la longueur d’onde critique 𝜆𝑐 vis-à-vis du

nombre de Péclet est illustrée sur la Figure ‎4.4. Les résultats de la stabilité linéaire (Eq. (4.7))
montre que l'évolution du nombre d'onde critique 𝑘𝑐 dépend du nombre de Péclet 𝒫, de

96
l’indice de comportement 𝑛 et du nombre de Hartmann ℋ. A partir de l’expression de 𝑘𝑐 , nous
pouvons déduire l’expression de la longueur d’onde critique 𝜆𝑐 donné par :

2𝜋 2 (4.8)
𝜆𝑐 = =
𝑘𝑐 √𝑆
La Figure ‎4.4 illustre ce comportement pour trois valeurs données de l’indice de la loi en
puissance 𝑛 et pour des valeurs de 𝒫allant de 0.01 à 10. Les mécanismes de dissipation
visqueuse varient avec l'intensité du débit filtrant, d'où nous observons la rupture de la
transformation galiléenne. En effet, pour 𝑛 = 1, 𝜆𝑐 et 𝑘𝑐 sont constants pour toute valeur de
Péclet. Cependant, pour les fluides de type loi en puissance, le nombre de Péclet rompt la
transformation galiléenne connue. L'équation (4.5) montre que le nombre d'onde critique
dépend de 𝒫 𝑛−1 . Un nombre de Péclet plus élevé rend le liquide moins visqueux pour 𝑛 < 1,
et plus visqueux pour 𝑛 > 1. C'est un aspect très important des fluides de puissance.

En outre, il apparaît sur la Figure ‎4.4-(b) que, pour 𝑛 < 1, la valeur critique de 𝜆𝑐 est une
fonction décroissante de 𝒫, cela signifie que le débit filtrant a un effet déstabilisateur. Pour
𝑛 > 1, la valeur critique de 𝜆𝑐 est une fonction croissante de 𝒫, de sorte que le débit filtrant a
un effet stabilisateur. Pour 𝑛 = 1, il n'y a pas d'effet sur les conditions critiques dues au débit
filtrant, en accord avec les résultats rapportés par Prats [59].

(a) (b)
2.00
2
5.0
H =2 H 2= 2

n=0.6 n=1.6
1.75
4.5

kc n=1 c 1.50 n=1

4.0

n=0.6
1.25
3.5 n=1.6

3.0 1.00
0 2 4 6 8 10 0 2 4 6 8 10
P P

Figure ‎4.4 : Evolution du nombre d’onde critique 𝒌𝒄 (a) et de la longueur d’onde critique 𝝀𝒄
(b) en fonction de 𝓟, pour différentes valeurs de 𝒏 et pour 𝓗𝟐 = 𝟐.
97
  Comportement asymptotique
Nous avons représenté sur la Figure ‎4.5, l’évolution du nombre de Rayleigh critique ℜ𝑐
en fonction du nombre de Hartmann ℋ, pour 𝒫 → 0 et pour différentes valeurs de 𝑛. On
observe que les résultats sont qualitativement identiques avec ceux du cas traité dans le
chapitre précèdent.

Cette figure montre que pour des petites valeurs de ℋ, l’évolution du nombre de Rayleigh
critique ℜ𝑐 dépend de ℋ et est fortement sensible à la valeur prise par l’indice de loi en
puissance 𝑛. La dépendance de la valeur critique du nombre de Rayleigh critique ℜ𝑐 vis-à-vis
l’indice de loi en puissance 𝑛 est illustré sur cette figure. Il est à noter que ℜ𝑐 diminue lorsque
𝑛 augmente, signalant la nature déstabilisante des effets visqueux. Lorsque ℋ < 100, la
dissipation par viscosité et effet Joule agissent ensemble. En revanche pour des valeurs de ℋ
très élevées (ℋ > 100), l’évolution du nombre de Rayleigh critique à différentes valeurs de 𝑛
tendent vers une valeur asymptotique qui ne dépend que de ℋ. Dans ce cas, on observe une
prédominance de la dissipation d'énergie par effet Joule. Cette dépendance en ℋ, lorsque
ℋ → ∞, découle de l’équation (4.6) et obéit à la loi de proportionnalité ℜ𝐶 ⟶ 𝜋 2 ℋ 2 . Dans
cette limite, en remplaçant les équations (2.25) dans l’équation (4.6) on peut déduire que la
différence de température critique ∆𝑇𝐶′ pour le début de l'instabilité devient indépendante de la
viscosité, de la perméabilité et de l’indice de loi de puissance. Elle dépend plutôt de la
conductivité électrique, de l'intensité du champ magnétique, de la gravité, de la hauteur de la
couche et des propriétés thermiques du liquide :

4𝜋 2 𝐵′2 𝜅𝛾
∆𝑇𝑐′ = (4.9)
𝜌0 𝑔𝛽𝑇′ 𝐻 ′

On constate que l’effet inhibiteur du champ magnétique sur le début de l’instabilité apparait à
partir de ce résultat. La présence d'un fort champ magnétique confère à un nouveau
comportement du fluide, et on ne peut plus le décrire comme newtonien ou non-newtonien,
puisque la viscosité n'a plus d'importance. Si nous augmentons la valeur de ℋ → ∞, nous
renforçons encore les effets magnétiques en conséquence le seuil de convection diminue et le
gradient de température critique devient beaucoup plus faible. Dans cette limite, la dissipation
d'énergie par effet Joule doit prédominer la dissipation d'énergie par contrainte de
cisaillement.

98
 100000
P 0

10000

c 1000

n=0.4
n=0.6
100
n=1
n=1.6

10
1 10 100
H

Figure ‎4.5 : Evolution du nombre de Rayleigh critique ℜ𝑐 en fonction de ℋ, pour 𝒫 → 0 et

pour 𝑛 > 1 (fluide dilatant), 𝑛 = 0.4 et 0.6 (fluides pseudoplastiques) et 𝑛 = 1 (fluide
newtonien).

4.3 Simulations numériques directes du problème non-linéaire

Nous avons vu au cours de la section précédente que l’étude des équations linéarisées
autour de l’état conductif de repos permettait de déterminer le seuil d’instabilité de cet état,
c’est-à-dire le mode critique (𝑘𝑐 ; ℜ𝑐 ) au-delà duquel l’état conductif est instable et fait place à
la convection.

L’étude de stabilité linéaire consiste à analyser la réponse d’un système par rapport à
des faibles perturbations. Lorsque les perturbations s’amplifient et deviennent plus
importantes, le traitement linéaire des instabilités n’est plus valable. Pour traiter le
comportement de ces instabilités dans le fluide, on considère l’approche d’analyse faiblement
non-linéaire. Ce choix est justifié par le fait que les termes d’ordre supérieur dans les
équations de perturbations augmentent plus rapidement que les termes d’ordre inférieur.

99
D’autre part, l’analyse linéaire s’avère insuffisante pour la détermination de la forme de
l’écoulement, même lorsque l’écart au seuil est infinitésimal. La résolution de ce type de
problème est possible par la voie numérique.

Dans cette section, nous allons présenter quelques simulations numériques directes des
équations non-linéaires pour une couche de dimensions 𝐴 × 𝑎 × 1 = 6 × 2 × 1.

Les simulations numériques ont été effectuées avec un code 3D basé sur la méthode des
volumes finis. Vu le grand nombre de nœuds que présente une simulation tridimensionnelle et
afin d’accélérer le calcul, une technique de maillage multigrille a été utilisé (Annexe A).

4.3.1 Validation du code numérique

En absence d’étude numérique dans la littérature traitant le cas Darcy-Bénard avec un

fluide non-newtonien de type loi en puissance et en présence d’un champ magnétique
appliqué. Le code élaboré a été validé en prenant comme référence l’étude théorique de
Barletta et Nield [58] traitant le cas Darcy-Bénard avec un fluide non-newtonien de type loi
en puissance et dans le cas des conditions aux limites thermiques de type Dirichlet. Nous
présentons dans ce qui suit les différentes confrontations effectuées en termes du nombre de
Nusselt 𝑁𝑢 en fonction du nombre de Rayleigh ℜ.

D’après la fFigure ‎4.6, on remarque que pour n=1, le nombre de Nusselt commence à croitre à
partir de ℜ𝑐 ≈ 40, valeur critique de la naissance de l’instabilité thermique dans le cas de la
convection naturelle d'un fluide Newtonien saturant un milieu poreux et dans le cas des
conditions aux limites de type Dirichlet. Pour 𝑛 = 1.3, le nombre de Nusselt commence à
croitre à partir de ℜ𝑐 = 0. Au-dessous de ces valeurs du nombre de Rayleigh critique, le taux
de transfert de chaleur 𝑁𝑢 est égal à 1 et donc le transfert de chaleur s’effectue uniquement
par conduction. Ces résultats sont en accord avec les résultats rapporté récemment par
Barletta et Nield [58].

100
 c c
2.2
n=1.3
2.0
n=1

1.8

1.6
Nu
1.4

1.2

1.0

0.8
0 20 40 60 80


Figure ‎4.6 : La variation du taux de transfert de chaleur Nu en fonction du nombre de

Rayleigh 𝕽 pour 𝓗=0 et A=6.

En présence d’un champ magnétique, le code a été validé à l’aide des résultats de la stabilité
linaire (Tableau ‎4.1) pour ℋ = √2. D’après la Figure ‎4.7, on remarque que le nombre de
Nusselt croît avec le nombre de Rayleigh. Cette augmentation commence avec la valeur du
nombre de Rayleigh critique ℜ𝑐 ≈ 74 pour n=1 (fluide newtonien) et ℜ𝑐 ≈ 20 pour n=1.3
(fluide dilatant). Au-dessous de ces valeurs de Rayleigh critique, le nombre de Nusselt à
l’intérieur de la couche poreuse est égal à 1 et le fluide est au repos.

101
 c c
1.8
n=1.3
n=1
1.6

1.4

Nu
1.2

1.0

0.8
0 20 40 60 80 100 120


Figure ‎4.7 : La variation du taux de transfert de chaleur Nu en fonction du nombre de

Rayleigh 𝕽 pour 𝓗=0 et A=6.

4.3.2 Analyse des résultats non linéaires

Les paramètre de contrôle tels que le nombre de Rayleigh ℜ, le nombre de Hartmann

ℋ, l’indice de la loi en puissance 𝑛 et le nombre de Péclet 𝒫 ont des effets considérables sur
le transfert de chaleur par convection. Les résultats obtenus par l’étude numérique seront
présentés dans cette partie.

Les résultats présentés dans cette étude couvrent la plage (0<ℜ<300), (0<ℋ< 104 ),
(0< 𝒫<10) et (0< 𝑛< 2). Les résultats numériques sont très proches des résultats de la stabilité
linéaire.

 Influence de l’indice de la loi en puissance 𝒏

 Structure de l'Ecoulement:

Pour étudier l’effet de l’indice de la loi en puissance 𝑛 sur la structure de l’écoulement

dans une cavité peu profonde, de rapport de forme 𝐴 = 6 et 𝑎 = 2, nous avons réalisé des
simulations numériques pour ℜ = 150, ℋ = 0, n = 1 et n = 1.6 respectivement. Les

102
résultats sont présentés en termes de lignes de courant sur la Figure ‎4.8. On remarque que la
convection, produite par des conditions thermiques de type Dirichlet, donne lieu à un
écoulement multicellulaire. En absence d’un champ magnétique (ℋ = 0) et dans le cas d’un
fluide newtonien (𝑛 = 1), les perturbations les plus instables sont des rouleaux transversaux
𝑅∥ (𝛼 = 0). Par contre dans le cas d’un fluide dilatant (𝑛 = 1.6), les perturbations les plus
instables sont des rouleaux longitudinaux 𝑅⊥ . Ces résultats sont en accord avec ceux trouvés
par Barletta et Nield [58].

Figure ‎4.8 : Lignes de courant au voisinage du point critique pour ℜ = 150, ℋ = 0, n = 1 ,

et n = 1.6.

D’autre part, l’effet de l’indice de loi en puissance est également présenté dans la Figure ‎4.9.
Les résultats sont présentés en termes de vecteurs vitesses (gauche) et isothermes (droite) pour
ℜ = 200, ℋ = 0 (a) 𝑛 = 1, (b) 𝑛 = 1.3, (c) 𝑛 = 1.6.
103
En comparant les isothermes produites avec différents indices de loi de puissance 𝑛, il est
clair que tous les résultats ont donné la même règle d'évolution. Les isothermes présentent des
stratifications diagonales dans la cavité et des gradients extrêmes près de la partie inférieure
de la paroi chaude et de la partie supérieure de la paroi froide dans tous les cas. Puisque la
cavité est horizontale, les particules de fluide chauffées descendent, sous l’action de la force
de flottabilité, près de la paroi chaude. En conséquence, les isothermes sont plus proches l'une
de l'autre près de la paroi chaude, indiquant un flux thermique de surface légèrement courbé
vers la paroi froide.

Bien que les isothermes produites avec les différents indices de loi de puissance 𝑛
apparaissent des modèles d'évolution similaires, les résultats obtenus avec le plus faible indice
de puissance semblent plus forts que les simulations produites avec l'indice de puissance plus
grand. Ainsi, le mouvement convectif devient moins important pour des fluides dilatants et
plus important pour des fluides newtoniens (𝑛 = 1). Cette notion est étroitement liée à la
variation de la viscosité effective du fluide donnée par 𝜂(𝛾̇ )𝑛−1 . En effet les fluides dilatants
sont des fluides rhéoépaississants, leur viscosité augmente lorsque le taux de cisaillement
augmente.

104
Figure ‎4.9 : Vecteurs vitesses (gauche) et isothermes (droite) pour 𝕽 = 𝟐𝟎𝟎, 𝓗 = 𝟎 (a)
𝒏 = 𝟏, (b) 𝒏 = 𝟏. 𝟑, (c) 𝒏 = 𝟏. 𝟔.

105
  Influence du nombre de Hartmann 𝓗

L’influence du nombre de Hartmann ℋ sur la structure de l’écoulement est illustré sur les
Figure ‎4.10 et Figure ‎4.11 pour les valeurs de ℋ = 0, √2 et 2. Ces figures montrent
l'importante réduction dans l’intensité globale de l’écoulement après l’application d’un champ
magnétique. Après l'augmentation de la valeur du nombre de Hartmann, une modification
remarquable de la structure des isothermes existe et de nouvelles recirculations apparaissent
entre la paroi froide et la paroi chaude.

Sur la Figure ‎4.10, on remarque que les isothermes sont stratifiées en diagonale pour ℋ = 0
et √2. tandis que pour ℋ = 2, les isothermes sont presque horizontalement stratifiées dans la
région centrale de la cavité. Une règle de changement similaire peut être observée pour les
lignes de courant de la Figure ‎4.10. De plus, comme nous pouvions nous y attendre, les
valeurs du seuil de convection sont beaucoup moins élevées en appliquant un champ
magnétique et dans le cas d’un fluide dilatant (𝑛 = 1.5). En effet, non seulement les effets
magnétiques tendent à stabiliser le système mais aussi les effets visqueux élevés; il faudra
donc fournir une beaucoup plus grande quantité d’énergie pour que les effets déstabilisants
l’emportent et que l’instabilité s’amplifie.

Le comportement dynamique de l'écoulement est illustré sur la Figure ‎4.11. On constate qu’en
augmentant le champ magnétique, les rouleaux tendent à s’installer dans la direction des
parois les plus allongés de la cavité. Par contre, si la valeur du nombre de Hartmann est nulle
ou très faible, aucune direction n’est privilégiée, dans ce cas les rouleaux s’installent en
diagonale et nous constatons un chevauchement des rouleaux longitudinaux et transversaux.
Pour une valeur de ℋ plus élevé ℋ = 2, nous constatons que le fluide choisit la forme des
rouleaux longitudinaux pour se déplacer, les rouleaux s’allongent sur le coté de rapport
d’aspect le plus grand. Ceci est traduit par un nombre minimum de rouleaux ce qui permet
d’économiser l’énergie et stabiliser l’écoulement. Ceci est dû à l'effet de la force
électromagnétique sur les vecteurs vitesses.

106
Figure ‎4.10 : Isothermes pour 𝕽 = 𝟐𝟎𝟎, 𝒏 = 𝟏 (gauche) 𝒏 = 𝟏. 𝟓 (droite) et pour trois
valeurs de 𝓗 = 𝟎, √𝟐 et 𝟐.

107
Figure ‎4.11 : Lignes de courant pour 𝕽 = 𝟐𝟎𝟎, 𝒏 = 𝟏 (gauche) 𝒏 = 𝟏. 𝟓 (droite) et pour
trois valeurs de 𝓗 = 𝟎, √𝟐 et 𝟐.

108
 Figure ‎4.12 : Nusselt locale dans le plan (X-Y) pour 𝕽 = 𝟐𝟎𝟎, 𝓗 = 𝟎, 𝟐 et 3, 𝒏 = 𝟏
(gauche) et 𝒏 = 𝟏. 𝟓 (droite).

D’autre part, les effets combinés du nombre de Hartmann et de l’indice de loi en puissance sur
le taux du transfert de chaleur sont également présentés sur les Figure ‎4.12 et Figure ‎4.13.

L’évolution du nombre de Nusselt local dans le plan (X-Y), pour différentes valeurs de ℋ et
𝑛, est représenté sur la Figure ‎4.12. Cette figure montre l'importante réduction dans le taux de
transfert de chaleur après l’application d’un champ magnétique. On constate qu’en
augmentant le nombre de Hartmann, les mouvements convectifs sont presque horizontalement
stratifiés dans la région centrale de la cavité et tendent à s’installer dans la direction des parois
les plus allongées. Lorsque le nombre de Hartmann augmente, le champ magnétique réduit
considérablement le nombre de Nusselt en supprimant les courants de convection.

109
 5.0
=250
4.5 =200
n=1.5
4.0
n=1
3.5

Numoy 3.0
2.5

2.0

1.5

1.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0
H

Figure ‎4.13 : Evolution du nombre de Nusselt moyen dans le plan (X-Y) en fonction du
nombre de Hartmann ℋ et du nombre de Rayleigh ℜ et pour 𝑛 = 1 et 𝑛 = 1.5.

Pour plus de précision et pour bien mettre en évidence l’influence du nombre de Hartmann,
nous avons présenté l’évolution du nombre de Nusselt moyen 𝑁𝑢 en fonction du nombre de
Hartmann, pour différentes valeurs de Rayleigh ℜ et pour 𝑛 = 1 et 𝑛 = 1.5 sur la Figure ‎4.13.
Les résultats indiquent clairement que la valeur du nombre de Nusselt moyen 𝑁𝑢 diminue
lorsque ℋ augmente. En outre, lorsque le nombre de Rayleigh augmente le nombre moyen de
Nusselt augmente également mais, et comme discuté ci-dessus, le mouvement convectif
devient moins important pour des fluides dilatants (𝑛 = 1.5) et plus important pour des
fluides newtoniens (𝑛 = 1). Cette notion est étroitement liée à l’augmentation de la viscosité
du fluide qui freine l’écoulement.

Le Tableau ‎4.2 confirme les résultats mentionnés ci-dessus. En augmentant la valeur de ℋ et

de 𝑛, le taux de transfert de chaleur diminuent considérablement.

110
  Influence du nombre de Rayleigh

Dans cette partie, nous étudions l'influence du nombre de Rayleigh sur la température et la
vitesse. Pour donner une analyse quantitative, nous résumons les simulations montrées sur la
Figure ‎4.14 en calculant la vitesse U et la température T le long de l’axe des Z.

1.0
1.0
=150 =150
=200 =200
0.8
=250 0.8 =250

0.6 0.6

Z Z
0.4 0.4

0.2 0.2

0.0
0.0
-3 -2 -1 0 1 2 0.0 0.2 0.4 0.6 0.8 1.0
U T
(a) n=1 (b)

1.0
1.0
=150 =150
=200 =200
0.8
=250 0.8 =250

0.6 0.6

Z Z
0.4 0.4

0.2 0.2

0.0
0.0
-0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 0.0 0.2 0.4 0.6 0.8 1.0

U T
(c) n=1.6 (d)

Figure ‎4.14 : Evolution de la température et de la composante de vitesse U en fonction de z

pour différentes valeurs de ℜ et 𝑛.
La Figure ‎4.14 montre que l'amplitude de U augmente notablement lorsque la valeur du
nombre de Rayleigh augmente, pour une valeur donnée de l'indice de loi en puissance,
indiquant l'augmentation du mouvement du fluide, qui a été démontrée qualitativement plus
111
tôt dans la Figure ‎4.9. De plus, l'augmentation de ℜ entraîne une diminution de l'épaisseur de
la couche limite thermique (c'est-à-dire une augmentation du gradient de température à
proximité des parois actives 𝑧 = 0).

Les valeurs du nombre de Nusselt moyen 𝑁𝑢 pour différentes valeurs du nombre de Rayleigh
de ℜ=200, 250 et 300 en présence du champ magnétique avec différentes valeurs de l’indice
de loi de puissance ont été analysées dans le Tableau ‎4.2. Comme il a été mentionné, les
valeurs de Nusselt moyen diminuent avec l'augmentation de l’indice de loi de puissance et du
nombre d’Hartmann ℋ, à toutes les valeurs du nombre de Rayleigh ℜ.

𝕽 𝒏 𝓗 𝑵𝒖 𝒏 𝓗 𝑵𝒖 𝒏 𝓗 𝑵𝒖
200 1 0 3.95 1.3 0 3.28 1.5 0 3.09
√2 2.76 √2 2.57 √2 2.51
2 2.15 2 2.13 2 2.1
250 1 0 4.55 1.3 0 3.64 1.5 0 3.36
√2 3.16 √2 2.91 √2 2.81
2 2.55 2 2.47 2 2.43
300 1 0 4.96 1.3 0 4.09 1.5 0 3.67
√2 3.46 √2 3.21 √2 3.02
2 2.88 2 2.76 2 2.68
Tableau ‎4.2 : Effet du nombre de Rayleigh, du nombre d’Hartmann 𝓗 et de l’indice de la loi
en puissance 𝒏 sur la valeur de Nusselt moyen.

 Effet de Péclet

La Figure ‎4.15 montre l’influence de 𝒫 sur les mouvements convectifs au sein du fluide pour
ℜ = 200, ℋ = 1. Les résultats obtenus illustrent l’existence d’un deuxième paramètre de
contrôle sur la stabilité, qui est le nombre de Péclet 𝒫 = 𝑅𝑒 ∗ 𝑃𝑟. Le nombre de 𝒫 a une
action stabilisante sur l’écoulement et les mouvements convectifs sont de moins en moins
importants.

Nous notons également, comme vu dans la stabilité linéaire, pour 𝑛 > 1, un débit de plus en
plus élevé entraine l’augmentation de la valeur de longueur d’onde critique 𝜆𝑐 , de sorte que le

112
débit filtrant a un effet stabilisateur. Nous observons un tel effet stabilisant parce que les
rouleaux transversaux ont disparu ainsi la zone d’entrée thermique occupe des rouleaux
longitudinaux plus étendus. Ces conclusions confirment les résultats de la stabilité linéaire.

Figure ‎4.15 : Structures convectives en convection naturelle 𝓟 = 𝟎 et en convection mixte

𝓟 = 𝟏𝟎 pour 𝕽 = 𝟐𝟎𝟎, 𝓗 = 𝟏.

 Le comportement asymptotique

Nous effectuons également des investigations numériques afin de traiter le

comportement asymptotique de ℜ pour les grandes valeurs de ℋ. Les résultats obtenus sont
présentés sur la Figure ‎4.16 en termes de 𝑁𝑢 en fonction de ℜ⁄ℋ 2 , pour un champ
magnétique faible (ℋ 2 = 4) sur la Figure ‎4.16-(a) et un champ magnétique fort (ℋ 2 = 104 )
dans la Figure ‎4.16-(b).

113
La Figure ‎4.16-(a) montre clairement que pour une faible valeur de ℋ (ℋ 2 = 4), il existe une
différence nette entre l’évolution des fluides newtoniens et dilatants en fonction de ℜ⁄ℋ 2 .
Cependant, pour une valeur très élevée de ℋ (ℋ 2 = 104 ), la Figure ‎4.16-(b) montre une
dépendance parfaite de 𝑁𝑢 sur ℜ⁄ℋ 2 , indépendamment de 𝑛. En effet, nous voyons que les
deux courbes ( 𝑛 = 1 𝑒𝑡 𝑛 = 1.6) se trouvent sur la même ligne, les deux types de fluides
s'écoulent exactement de la même manière et on ne peut plus décrire le fluide comme
newtonien ou non-newtonien. A cette limite, la viscosité est inhibée par l'énergie magnétique,
cela signifie que, à une telle plage de nombre de Hartmann, la dissipation par effet Joule
prédomine.

On peut déduire à partir de ces résultats que le comportement asymptotique est toujours valide
même si nous considérons les termes non linéaires ce qui confirme les résultats de l'analyse de
stabilité. Un scénario similaire a été également révélé dans le ‎Chapitre 3.

2
(a) H = 4
5 5

n=1.6 Nu
4 n=1 4

3 3

Nu
2 2

1 1

0 0
0 20 40 60 80

/H 2

114
 2 4
(b) H = 10

3.5
n=1.6 Nu
3.0
n=1

2.5

2.0
Nu
1.5

1.0

0.5

0.0
12 14 16 18 20

/H 2

Figure ‎4.16 : 𝑁𝑢 en fonction de ℜ⁄ℋ 2 et 𝑛 pour: (a) ℋ 2 = 2 et (b) ℋ 2 = 104 .

4.4 Conclusion
Dans ce chapitre, la convection au sein d’une une couche poreuse horizontale d’un
fluide non-newtonien dans le cas de condition thermique de type Dirichlet a été étudiée. Des
conditions de températures sont imposées sur les parois horizontales (parallèle au plan
(𝑜𝑥’y’) et les autres parois sont adiabatiques. Un champ magnétique constant est appliqué
parallèlement à la gravité 𝑔. Un modèle de de type loi en puissance a été utilisé pour
caractériser le comportement non-newtonien du fluide. L’effet de l’indice de la loi en
puissance, du nombre de Hartmann, du nombre de Rayleigh et du nombre de Péclet sur les
écoulements et les transferts de chaleur ont été présentés pour une large gamme des
paramètres régissants le problème.

L'analyse de stabilité linéaire a été utilisé pour déterminer le mode critique (𝑘𝑐 , ℜ𝑐 ), au-
delà duquel l’écoulement cesse d’être stable. Les résultats trouvés sont similaires à ceux
observés dans le chapitre précédent dans le cas des conditions aux limites thermiques de type

115
Neumann. En effet, pour les fluides newtoniens, l'apparition de l'instabilité convective
s’effectue pour ℜ𝑐 > 4𝜋², elle dépend du nombre de Hartmann ℋ. La présence du champ
magnétique retarde l’apparition de la convection. Pour les fluides dilatants, la présence du
champ magnétique augmente la valeur classiques du nombre de Rayleigh pour l’apparition de
𝑠𝑢𝑝
la convection de zéro à ℜ𝐶 = 𝜋 2 ℋ 2 . Pour les fluides pseudoplastiques, indépendamment de
ℋ, le système est inconditionnellement stable pour des perturbations infinitésimales.

En présence d'un fort champ magnétique, un résultat important a été trouvé concernant
le caractère non visqueux des fluides de type loi en puissance. Une solution numérique basée
sur la méthode des volumes finis est venue conforter les résultats analytiques.

116
 CONCLUSION GENERALE ET
PERSPECTIVES

Dans cette thèse nous avons étudié l’effet d’un champ magnétique constant sur la
convection naturelle et mixte d’un fluide non-newtonien de type loi en puissance
électriquement conducteur saturant une couche poreuse très mince d’une configuration
géométrique rectangulaire. Deux conditions aux limites thermiques ont été étudiées, à savoir
le cas d’un écoulement thermo-convectif induit par un flux constant (Condition de type
Neumann) et celui induit par une température constante (Condition de type Dirichlet),
imposées sur les parois horizontales de la cavité.

Le phénomène de la convection naturelle a été modélisé en se basant sur

l’approximation de Boussinesq et en utilisant un modèle rhéologique de type loi en puissance
pour simuler le comportement visqueux du fluide non-newtonien.

Pour résoudre notre problème nous avons développé trois approches de calcul :

1) Etude de la stabilité linéaire 3D qui montre que les structures dynamiques sont
tridimensionnelles dans le cas des conditions de température, et bidimensionnelles de
type écoulement parallèle dans le cas des conditions de flux.
2) Une méthode numérique, basée sur la méthode des différences finies et la méthode des
volumes finis, sous différentes conditions aux limites (température et flux) des
régimes : sous critique et supercritique. Les codes de calcul (en Fortran) élaborés pour
résoudre les équations régissant le problème ont été validés avec des résultats
numériques et analytiques disponibles dans la littérature.
3) Une méthode analytique de type écoulement parallèle valable pour les cavités mince
(A>>1) soumises à un flux constant. Quelques solutions particulières ont été trouvées
(n=0.5, n=1 et n=2). Dans la plupart des cas la solution n’est possible que par des
voies numériques.

117
L’analyse de stabilité linéaire, l’approximation de l’écoulement parallèle ainsi que les
simulations numériques utilisées comme outils d’investigation, ont permis d’aboutir aux
conclusions suivantes :

- Pour les fluides pseudoplastiques, la couche est inconditionnellement stable. Le

nombre de Rayleigh critique pour l’apparition de la convection, indépendamment de
ℋ, est égal à l’infini. En revanche, les résultats numériques et l’approximation de
l’écoulement parallèle montrent que ce résultat prédit par la théorie de stabilité linéaire
n’est pas réaliste et que le mouvement convectif est possible au-dessus d'un nombre de
Rayleigh sous-critique ℜ𝑐𝑠𝑢𝑏 .
- La présence du champ magnétique retarde l’apparition de la convection (joue un rôle
stabilisant) et le nombre de Rayleigh critique dépend du ℋ pour le cas d’un fluide
dilatant.
- Dans la limite d'un champ magnétique très puissant, le nombre de Rayleigh critique
ℜ𝑐 a une dépendance asymptotique vis-à-vis du nombre de Hartmann ℋ,
indépendamment de l’indice de loi de puissance 𝑛. Cela implique que, dans cette
limite, le gradient de température critique pour l'apparition de l'instabilité devient
indépendant de la viscosité et dépend plutôt du coefficient de conductivité électrique 𝛾
et de la force du champ magnétique car la dissipation d'énergie par effet Joule
prédomine la dissipation d'énergie par contrainte de cisaillement.
- Dans la limite d'un champ magnétique très puissant, le problème montre que la
viscosité a une dépendance asymptotique d’une nouvelle variable ℜ⁄ℋ 2 . Le
comportement du fluide devient indépendant de la viscosité, et les forces de Lorentz
deviennent l’acteur principal de dissipation de l’énergie.

Perspectives
Les perspectives envisagées sont nombreuses :

- La continuité des travaux de recherche entamés sur le milieu poreux en utilisant le

modèle de Darcy-Brinkman-Forcheimer qui tient compte des effets visqueux, lorsque
la perméabilité est importante, et des effets d’inertie.

118
- Viser à étudier l’effet stabilisateur d’un champ magnétique sur la convection d’autres
modèles de fluides non-newtoniens, comme les fluides viscoélastiques.
- Prendre en compte l’effet d’inclinaison d’une cavité et l’effet Soret sur ces
écoulements.
- Considérer des configurations géométriques autres que la cavité rectangulaire. Par
exemple le cas de plaques, de cylindres, etc. Il serait aussi intéressant de considérer le
cas de la convection induite par un cylindre en rotation.
- Vérification des résultats obtenus par voie expérimentale.

119
 Références bibliographiques

[1] Rabinow J. (1948), The magnetic fluid clutch. AIEE Trans vol.67:1308-1315.
[2] Spencer BF., Dyke SJ., Sain MK., Carlson JD. (1997), Phenomenological model for
magnetorheological dampers. J Eng Mech vol.123:230-238.
[3] Kordonski WI., Golini D. (1999), Fundamentals of magnetorheological fluid utilization in
high precision finishing. J Intell Mater Syst Struct vol.10:683-689.
[4] Zitha PLJ. (2004), Method of drilling with magnetorheological fluid. US patent number
7021406.
[5] Rankin PJ., Horvath AT., Klingenberg DJ. (1999), Magnetorheology in viscoplastic
media. Rheol Acta vol.38:471-477.
[6] Park BO., Park BJ., Hato MJ., Choi HJ. (2011), Soft magnetic carbonyl iron microsphere
dispersed in grease and its rheological characteristics under magnetic field. Colloid Polym Sci
vol.289:381-386.
[7] De Vicente J., Lopez-Lopez MT., Gonzalez-Caballero F., Duran JDG. (2003), Rheological
study of the stabilization of magnetizable colloidal suspensions by addition of silica
nanoparticles. J Rheol vol.47:1093-1109.
[8] Lopez-Lopez MT., Zugaldia A., Gonzalez-Caballero F., Duran JDG. (2006),
Sedimentation and redispersion phenomena in iron-based magnetorheological fluids. J Rheol
vol.50:543-560.
[9] Tao R., Huang K. (2011), Reducing blood viscosity with magnetic fields. Phys. Rev. E
vol.84:1-5.
[10] Ostwald W., Kolloid Z. (1925), 36-99.
[11] Horton CW., Rogers FT. (1945) Convection currents in a porous medium. J. Appl. Phys.
vol.16:367-370.
[12] Lapwood ER. (1948), Convection of a fluid in a porous medium. Proc. Cambridge
Philos. Soc. vol.44:508-521.
[13] Nield D.A. (1968), Onset of thermohaline convection in a porous medium. Water
Ressources Res. vol.3:553-560.
[14] Ben Hamed H., Bennacer R. (2008), Analytical development of disturbed matrix
eigenvalue problem applied to mixed convection stability analysis in Darcy media. C. R.
Mécanique, vol.336:656-663.
[15] Ben Hamed H., Bennacer R., Beji H., Tanglet T. (2011), Linear stability analysis of
Horton Rogers Lapwood problem under Soret effect. Int. J. of Dynamics of Fluids vol.7:55–
74.

120
[16] Koschmieder E.L., Pallas S.G. (1974), Heat transfer through a shallow horizontal
convecting fluid layer. Int. J. heat Mass Transfer vol.17:991-1002.
[17] Bejan A. (1984), Convection Heat Transfer. Wiley, New York.
[18] Bories S., Mojtabi A., Prat M., Quintard M. (2010), Transfert de chaleur dans les milieux
poreux- conduction, convection, rayonnement. Technique de l’ingénieur, Réf. 42214210.
[19] Nield A., Bejan A. (2013), Convection in Porous Media. 4rd ed. Springer, New York.
[20] Rouquerol J. et al. (1994), Recommendations for the characterization of porous solids
(Technical Report). Pure & Appl. Chem. vol.66:1739-1758.
[21] Kaviany M. (1995), Principles of heat transfer in porous media. Springer (cf. p. 5).
[22] Kozeny J. (1927), Flow in porous media.S.B. Akad. Wiss, 2a: 126.

[23] Ergun S. (1927), Fluid flow through packed columns. Chem. Engineering Prog., vol.48:
89-94.
[24] Guyon É., Hulin J.P., Petit L. (2001), Hydrodynamique physique. EDP Sciences, CNRS
Éd.
[25] Chhabra R.P., Richardson J.F. (1999), Non-Newtonian flow in the process industries.
Fundamentals and Engineering Application.
[26] Source web : S. Poncet Initiation à la rhéologie, “http ://l3mgp.l3m.univmrs.
fr/site/SitePersoPoncet/enseignement/cours_rheologie.pdf”
[27] Cours de Mécanique des fluides - José Bico, Marc Fermigier, Pascal Kurowski, ESPCI,
Paris.

[28] Boger D.V., Walters K. (1993), Rheological phenomena in focus. Elsevier, Amsterdam.
[29] Pascal H. (1983), Rheological behaviour effect of non-Newtonian fluids on steady and
unsteady flow through porous media. Int. J. Numerical Anal. Methods Geomech. vol.7:207-
224.
[30] Chandrasekhar S. (1954), On the inhibition of convection by magnetic field. II, Phil.
Mag. Ser. 7, vol. 45:1177-1191.
[31] Chandrasekhar S. (1956), On the stability of simplest solution of the equation of
hydromagnetics. Proc. Nat. Aca. Sci. vol.42: 273-276.
[32] Chandrasekhar S. (1961), Hydrodynamic and hydromagnetic stability. Oxford University
press.
[33] Jirlow K. (1956), Experimental investigation of the inhibition of convection by a
magnetic field. ibid. vol.8: 252-253.
[34] Lehnert B., Little N. C. (1957), Experiments on the effect of inhomogeneity and
obliquity of a magnetic field in inhibiting convection. Tellus vol.9:97-103.

121
[35] Nakagawa Y. (1957), Experiments on the inhibition of thermal convection by a magnetic
field. Proc. Roy. Soc. (London) A, vol.240: 108-113.
[36] Raptis A., Massalas C., Tzivanidis G. (1982), Hydromagnetic free convection through a
porous medium between two parallel plates. Physics letters A, vol.90: 288-289.
[37] Alchaar S., Vasseur P., Bilgen E. (1995), Natural convection heat transfer in a
rectangular enclosure with a transverse magnetic field. J. Heat Transfer, vol.117: 668-673.
[38] Zeng M., Wang Q.W., Ozoe H. (2009), Natural convection of a diamagnetic fluid in an
enclosure filled with a porous medium under magnetic field. Prog Comput Fluid Dyn
vol.9:77-85.
[39] Bian W., Vasseur P., Bilgen E., Meng F. (1996), Effect of an electromagnetic field on
natural convection in an inclined porous layer. Int J Heat Fluid Flow vol.17:36-44.
[40] Ece M.C., Buyuk E. (2006), Natural-convection flow under a magnetic field in an
inclined rectangular enclosure heated and cooled on adjacent walls. Fluid Dyn Res
vol.38:564-590
[41] Robillard L., Bahloul A., Vasseur P. (2006), Hydromagnetic natural convection of a
binary fluid in a vertical porous enclosure. Chem. Eng. Comm., under press.
[42] Revnic C., Grosan T., Pop I., Ingham D.B. (2011), Magnetic effect on the unsteady free
convection flow in a square cavity filled with a porous medium with a constant heat
generation. Int J Heat Mass Transf vol.54:1734-1742.
[43] Alloui Z., Vasseur P. , Costa V. A. F. , Sousa A. C. M. (2013), Effect of a non-constant
magnetic field on natural convection in a horizontal porous layer heated from the bottom. J
Eng Math vol.81:141-155.
[44] Pangrle BJ, Walsh EG, Moore SC, Dibiasio’o D (1992), Magnetic resonance imaging of
laminar flows in porous tube and shell systems. Chem Eng Sci vol.47:517-526.
[45] McWhiter J.D., Crawford M.E., Klein D.E. (1998), Magnetohydrodynamic flows
through porous media II: experimental results. Fusion Sci Technol vol.34:187-197.
[46] Khuzhir P., Bossis G., Bashtovoi V., Volkova O. (2003), Flow of magnetorheological
fluids through porous media. Eur J Mech B vol.22:331-343.
[47] Malvandi A., Ganji D.D. (2014), Magnetic field effect on nanoparticles migration and
heat transfer of water/alumina nanofluid in a channel. J.Mag.Magn.Mater. vol.362: 172-179.
[48] Sheikholeslami M., GorjiBandpy M., Ellahi R., Hassan M., Soleimani S. (2014), Effects
of MHD on Cu–waternanofluid flow andheat transfer by meansof CVFEM. J.Magn.Magn.
Mater. vol.349:188-200.
[49] Kefayati GH.R. (2016), Simulation of heat transfer and entropy generation of MHD
natural convection of non-Newtonian nanofluid in an enclosure. Int J Heat Mass Transf
vol.92:1066-1089.
[50] Chen H. T., Chen C. K. (1988), Free convection flow of non-Newtonian fluids along a
vertical plate embedded in a porous medium. ASME Trans. J. Heat Transfer vol.110:257-260.
122
[51] Amari B., Vasseur P., Bilgen E. (1994), Natural convection of non-Newtonian fluids in a
horizontal porous layer. Heat Mass Transfer vol.29:185-193.
[52] Bian W., Vasseur P., Bilgen E. (1994), Boundary-layer analysis for natural convection in
a vertical boundary layer filled with a non-Newtonian. Int. J. Heat Fluid Flow vol.15: 384-
391.
[53] Pascal J. P., Pascal H. (1997), Free convection in a non-Newtonian fluid saturated porous
medium with lateral mass flux. Int. J. Non-linear Mech. vol.32:471-482.
[54] Chen G., Hadim H.A. (1998), Forced convection of a power-law fluid in a porous
channel Numerical solutions. Heat Mass Transfer vol.34:221-228.
[55] Nield D.A., Kuznetsov A.V. (2005), Thermally developing forced convection in a
channel occupied by a porous medium saturated by a non-Newtonian fluid. Int. J. Heat Mass
Transfer vol.48:1214-1218.
[56] Getachew D., Poulikakos D., Minkowycz (1998), Double diffusion in a porous cavity
saturated with a non-Newtonian fluid. J. Thermophysics and Heat Transfer vol.38:437-446.
[57] Nield D.A. (2011), A further note on the onset of convection in a layer of a porous
medium saturated by a non-Newtonian fluid of power-law type. Transp. Porous Med.
vol.88:187-191.
[58] Barletta A., Nield D.A. (2011), Linear instability of the horizontal through flow in a
plane porous layer saturated by a power-law fluid. Phys. Fluids vol.23, 013102-7.
[59] Prats M. (1966), The effect of horizontal fluid flow on thermally induced convection
currents in porous mediums. J. Geophys. Res. vol.71:4835-4838.
[60] Nield A. and Bejan A. (2006), Convection in Porous Media. 3rd ed. Springer, New York.
[61] Ben Khelifa N., Alloui Z., Beji H., Vasseur P. (2012), Natural convection in a horizontal
porous cavity filled with a non-Newtonian binary fluid of power-law type. J. non-Newtonian
Fluid Mechanics, Vol. 169-170, 15-25.
[62] Gray D. D., Giorgini A. (1976), The Validity of the Boussinesq Approximation for
Liquids and Gases. Int J Heat Mass Transf, vol.19:545-551.
[63] Nield D., Bejan A. (1999), Convection in Porous Media. 2nd ed. Spinger, New York.
[64] Combarnous M., Bories S. (1974), Modelisation de la Convection Naturelle au sein
d’une Couche Poreuse Horizontale à l’aide d’un Coefficient de Transfert Solide-Fluide. Int J
Heat Mass Transf vol.17:505-515.

[65] Squire H. B. (1933), On the stability for three-dimensional disturbances of viscous fluid
flow between parallel walls. Proc. Royal Soc.
[66] Cormack D.E., Leal L.G., Imberger J. (1974), Natural convection in a shallow cavity
with differentially heated end walls. Part.1, Asymptotic theory, J. Fluid Mechanics
vol.65:209-230.

123
[67] Bejan A., Tien, C. L. (1978), Natural Convection in a Horizontal Porous Medium
Subjected to an End-to-End Temperature Difference. Journal of Heat Transfer, vol.100:191-
198.
[68] Trevisan O., Bejan A., (1986), Mass and heat transfer by natural convection in a vertical
slot filled with porous medium. Int. J. Heat Mass Transfer vol.29:403-415.
[69] Sen M., Vasseur P., Robillard L. (1987), Multiple Steady for Unicellular Natural
Convection in a Inclined Porous Layer. Int. J. Heat Mass Transfer, vol. 30: 2097-2113.
[70] Vasseur P., Wang CH., Sen M. (1989), The Brinkman model for natural convection in a
shallow porous cavity with Uniform Heat Flux. Numerical Heat Transfer vol.15: 221-242.
[71] Alavyoon F., Masuda Y., Kimura S. (1994), On natural convection in vertical porous
enclosures due to opposing fluxes of heat and mass prescribed at the vertical walls. Int .J.
Heat Mass Transfer vol.37:195-206.
[72] Alavyoon F. (1993), On natural convection in vertical porous enclosures due to
prescribed fluxes of heat and mass at the vertical boundaries. Int. J. Heat Mass Transfer
vol.36:2479-2498.
[73] Mamou M., Vasseur P., Bilgen E. (1995), Multiple solutions for double diffusive
convection in a vertical porous enclosure. Int. J. Heat Mass Transfer vol.38:1787-1798.
[74] Kalla L., Mamou M., Vasseur P., Robillard L. (1999), Multiple steady states for natural
convection in a shallow porous cavity subject to uniform heat fluxes. Int. Comm. Heat Mass
Transfer vol.26:761-770.
[75] Bian W., Vasseur P., Bilgen E. (1994), Natural convection of non-Newtonian fluids in an
inclined porous layer. Chem. Eng. Commun. vol.129:79-97.
[76] Alloui Z., Vasseur P. (2012), Onset of Rayleigh–Bénard MHD convection in a
micropolar fluid. Int. J. Heat Mass Transfer vol.55:1164-1169.
[77] Vasseur P., Satish M.G., Robillard, L. (1987) Natural convection in a thin, inclined,
porous layer exposed to a constant heat flux. Int. J. Heat Mass Transfer vol.30:537-548.

[78] Zhou X., Yu M. (2012), The structural flow in pipe containing porous medium saturated
with power-law fluid. Journal of Hydrodynamics: 55.

[79] Mamou M., Vasseur P. Bilgen E., Galerkin A. (1998), Finit element study of the onset of
double-diffusive convection in an inclined porous enclosure. Int. J. Heat Mass Transfer,
vol.41:1513-1529.
[80] Boutana N., Bahloul A., Vasseur P., Joly F. (2004), Soret-driven and double-diffusive
natural convection in a vertical porous cavity. Journal of Porous Media, vol.7(1):1-17.
[81] Frankel, S. P. (1950), Convergence rates of iterative treatments of partial differential
equations. Math Tables and Other Aids to Computations vol.4:65-75.
[82] Roache P.J. (1985), Computational fluid dynamics. Hermosa publishers.

124
[83] Patankar N. V. (1980), Numerical Heat Transfer and Fluid Flow. Hemisphere,
Washington.
[84] Bennacer R., Tobbal A., Beji H., Vasseur P. (2001), Double diffusive convection in a
vertical enclosure filled with anisotropic porous media. Inter. J. Ther. Sci. 2001 vol.40:30-41
[85] Hackbusch W. (1985), Multi-Grid Methods and applications. New York:Springer-
Verlag.
[86] Hackbusch W., Trottenberg U. (1991), Multigrid Methods III. Boston: Birkhauser.
[87] Horthmann M., Peric M., Scheuerer G. (1990), Finite volume multigrid prediction of
laminar natural convection: Bench-mark solutions. Int. J. Numer. Methods Fluids vol.11:189.
[88] Wesseling P. (1992), An introduction to multigrid Method. New York: Wiley.
[89] Pellew A., Southwell R. V. (1940) Proc. Roy. Soc. A, 176: 312.

125
 Annexe A
A.1 Méthode des différences finies
La méthode des différences finies a été décrite en détails par plusieurs auteurs [79],
[80]. Pour cette raison seulement les points saillants de ce code seront discutés ci-dessous.

Le choix de la méthode des différences finies est basé sur ses avantages dans sa grande
simplicité d’écriture et son faible coût de calcul, bien qu’elle présente l’inconvénient liée à sa
limitation à des géométries simples.

Dans le cas 2D instationnaire, la méthode aux différences finies consiste à discrétiser le

système physique par le biais d’un maillage uniforme suivant les deux directions de l’espace
(figure A.1), décomposé en 𝑁 × 𝑃 noeuds (𝑥𝑖 , 𝑧𝑗 ) répartis régulièrement avec un pas d'espace
∆𝑥 dans la direction 𝑥 et ∆𝑧 dans la direction 𝑧. Ainsi, on calcule les valeurs de la
température, de la fonction de courant et de la viscosité apparente sur chaque nœud, en
𝑛
résolvant les équations gouvernantes. On notera 𝑢𝑖𝑗 la valeur discrète de la grandeur 𝑢(𝑥, 𝑧, 𝑡)
au nœud (𝑥𝑖 , 𝑧𝑗 ) et au temps 𝑛∆𝑡.

Figure A.1 : Représentation du maillage de notre modèle physique

126
A.1.1 Algorithmes de résolution

L’équation de l'énergie est résolue par la méthode implicite aux directions alternées
(A.D.I) modifiée par Peaceman et Rachford en 1955. En utilisant un schéma centré de
deuxième ordre et la forme conservatrice pour les termes convectifs. La méthode ADI donne
lieu à un système matriciel tridiagonal à résoudre. La solution est obtenue en balayant le
domaine de calcul dans la direction des 𝑥 puis dans celle de 𝑧. La résolution de l’équation de
la fonction de courant est assurée par la méthode de surrelaxation successive SOR, proposée
par Frankel [81].

Pour obtenir la solution de ces équations, après un nombre suffisant d'itérations, un

critère de convergence relatif à la fonction de courant et la température doit être satisfait. Ce
critère est défini par :

𝑚+1 𝑚
∑𝑖 ∑𝑗[𝜙𝑖,𝑗 − 𝜙𝑖,𝑗 ]
𝑚+1 ≤𝜀 (A.1)
∑𝑖 ∑𝑗 𝜙𝑖,𝑗
Où 𝜙 représente 𝜓 et T, m représente le nombre d’itérations et 𝜀 = 10−6 . La faible valeur de
𝜀 ne change rien dans le résultat final.

Les principales étapes de l’algorithme de calcul sont :

1. Génération du maillage

2. Introduction des conditions initiales sur la température T sur la fonction de courant 𝜓

3. Evaluation du champ de vitesse initiale 𝑢 et 𝑤

4. Evaluation du champ de viscosité apparente 𝜇𝑎

5. Début de la boucle dans le temps

6. Résoudre l’équation de l'énergie par la méthode ADI

7. Résoudre l'équation de la fonction de courant par la méthode SOR

8. Calcul des champs de vitesse et de viscosité apparente

9. Retour à l’étape 5 et répéter les opérations 6-8 jusqu’à la convergence.

127
A.1.2 Discrétisation des équations

Dans le domaine discrétisé, les dérivées partielles du premier et deuxième ordre sont
approchées selon un schéma aux différences finies centrées, ils prennent les expressions
suivantes:

𝜕𝑓 𝑓𝑖+1,𝑗 − 𝑓𝑖−1,𝑗 𝜕𝑓 𝑓𝑖,𝑗+1 − 𝑓𝑖,𝑗−1

) = ) =
𝜕𝑥 𝑖,𝑗 2∆𝑥 𝜕𝑧 𝑖,𝑗 2∆𝑧
(A.2)
𝜕2𝑓 𝑓𝑖+1,𝑗 − 2𝑓𝑖,𝑗 + 𝑓𝑖−1,𝑗 𝜕2𝑓 𝑓𝑖,𝑗+1 − 2𝑓𝑖,𝑗 + 𝑓𝑖,𝑗−1
2
) = ) =
𝜕𝑥 𝑖,𝑗 ∆𝑥 2 2
𝜕𝑧 𝑖,𝑗 ∆𝑧 2

La méthode ADI donne lieu à des systèmes matriciels tridiagonaux dans les directions x et z.
Le système dans la direction de x est obtenu en faisant une discrétisation implicite par rapport
à x et explicite par rapport à z, tandis que le système dans la direction de z est obtenu en
faisant une discrétisation implicite par rapport à z et explicite selon x.
 Discrétisation de l’équation de l’énergie
L’équation d’énergie (2.16) se réécrit de la façon suivante :

𝜕𝑇 𝜕(𝑢𝑇) 𝜕(𝑤𝑇) 𝜕 2 𝑇 𝜕 2 𝑇
+ + = 2+ 2 (A.3)
𝜕𝑡 𝜕𝑥 𝜕𝑧 𝜕𝑥 𝜕𝑧

-Schéma implicite en x et explicite en z

C’est un schéma aux différences centrées, au temps (𝑛 + 1/2) pour les dérivées en x et au
temps (n) pour les dérivées en z, d’où :

𝑛+1/2 𝑛
𝜕𝑇 𝑇𝑖,𝑗 − 𝑇𝑖,𝑗
) =
𝜕𝑡 𝑖,𝑗 (∆𝑡⁄2)

𝑛 𝑛+1/2
𝑛 𝑛+1/2 𝑛 𝑛 𝑛 𝑛
𝜕(𝑢𝑇) 𝑢𝑖+1,𝑗 𝑇𝑖+1,𝑗 − 𝑢𝑖−1,𝑗 𝑇𝑖−1,𝑗 𝜕(𝑤𝑇) 𝑤𝑖,𝑗+1 𝑇𝑖,𝑗+1 − 𝑤𝑖,𝑗−1 𝑇𝑖,𝑗−1
) = ) = (A.4)
𝜕𝑥 𝑖,𝑗 2∆𝑥 𝜕𝑧 𝑖,𝑗 2∆𝑧

𝑛+1/2 𝑛+1/2 𝑛+1/2 𝑛 𝑛 𝑛

𝜕2𝑇 𝑇𝑖+1,𝑗 − 2𝑇𝑖,𝑗 + 𝑇𝑖−1,𝑗 𝜕2𝑇 𝑇𝑖,𝑗+1 − 2𝑇𝑖,𝑗 + 𝑇𝑖,𝑗−1
2
) = ) =
𝜕𝑥 𝑖,𝑗 ∆𝑥 2 𝜕𝑧 2 𝑖,𝑗 ∆𝑧 2

En substituant les expressions (A.4) dans l'équation (A.3), on obtient le système algébrique :

128
 1 1 1
𝑛+ 𝑛+2 𝑛+
𝐴𝑖 𝑇𝑖−1,𝑗2 − 𝐵𝑖 𝑇𝑖,𝑗 + 𝐶𝑖 𝑇𝑖+1,𝑗2 = 𝐷𝑖 (A.5)

Avec

𝑛 𝑛
𝑢𝑖−1,𝑗 1 2 2 1 𝑢𝑖+1,𝑗
𝐴𝑖 = + 2 𝐵𝑖 = − − 2 𝐶𝑖 = 2 −
2∆𝑥 ∆𝑥 ∆𝑡 ∆𝑥 ∆𝑥 2∆𝑥
(A.6)
𝑛 𝑛
1 𝑤𝑖,𝑗−1 𝑛
2 2 𝑛
1 𝑛
𝑤𝑖,𝑗+1
𝐷𝑖 = ( − ) 𝑇𝑖,𝑗−1 + (− + ) 𝑇𝑖,𝑗 + (− + ) 𝑇𝑖,𝑗+1
∆𝑧 2 2∆𝑧 ∆𝑡 ∆𝑧 2 ∆𝑧 2 2∆𝑧

-Schéma implicite en z et explicite en x

C’est un schéma aux différences centrées, au temps (𝑛 + 1) pour les dérivées en 𝑧 et au temps
(𝑛 + 1/2) pour les dérivées en 𝑥, d’où :

1
𝑛+1 𝑛+2
𝜕𝑇 𝑇𝑖,𝑗 − 𝑇𝑖,𝑗
) =
𝜕𝑡 𝑖,𝑗 (∆𝑡⁄2)

1 1 1 1
𝑛+ 𝑛+ 𝑛+ 𝑛+ 𝑛+1/2 𝑛+1/2 𝑛+1/2
𝜕(𝑢𝑇) 𝑢𝑖+1,𝑗 𝑇𝑖+1,𝑗2
2
− 𝑢𝑖−1,𝑗 𝑇𝑖−1,𝑗2
2
𝜕2𝑇 𝑇𝑖+1,𝑗 − 2𝑇𝑖,𝑗 + 𝑇𝑖−1,𝑗
) = ) =
𝜕𝑥 𝑖,𝑗 2∆𝑥 𝜕𝑥 2 𝑖,𝑗 ∆𝑥 2 (A.7)

𝑛+1/2
𝑛+1 𝑛+1 𝑛+1/2 𝑛 𝑛 𝑛
𝜕(𝑤𝑇) 𝑤𝑖,𝑗+1 𝑇𝑖,𝑗+1 − 𝑤𝑖,𝑗−1 𝑇𝑖,𝑗−1 𝜕2𝑇 𝑇𝑖,𝑗+1 − 2𝑇𝑖,𝑗 + 𝑇𝑖,𝑗−1
) = ) =
𝜕𝑧 𝑖,𝑗 2∆𝑧 𝜕𝑧 2 𝑖,𝑗 ∆𝑧 2

En substituant les expressions (A.7) dans l'équation (A.3), on obtient le système algébrique :

𝑛+1 𝑛+1 𝑛+1

𝐴𝑗′ 𝑇𝑖,𝑗−1 − 𝐵𝑗′ 𝑇𝑖,𝑗 + 𝐶𝑗′ 𝑇𝑖,𝑗+1 = 𝐷𝑗′ (A.8)

Avec

129
 𝑛+1/2 𝑛+1/2
𝑤𝑖,𝑗−1 1 2 2 1 𝑤𝑖,𝑗+1
𝐴𝑗′ = + 2 𝐵𝑗′ =− − 2 𝐶𝑗′ = 2−
2∆𝑧 ∆𝑧 ∆𝑡 ∆𝑧 ∆𝑧 ∆𝑧
(A.9)
𝑛+1/2 𝑛+1/2
1 𝑢𝑖−1,𝑗 𝑛+1/2 2 2 𝑛+1/2 1 𝑛+1/2
𝑢𝑖+1,𝑗
𝐷𝑗′ = ( 2
− ) 𝑇𝑖−1,𝑗 + (− + 2 ) 𝑇𝑖,𝑗 + (− 2 + ) 𝑇𝑖+1,𝑗
∆𝑥 2∆𝑥 ∆𝑡 ∆𝑥 ∆𝑥 2∆𝑥

 Discrétisation de l’équation de mouvement

Le champ de température obtenu sur l’ensemble des nœuds du domaine est utilisé pour
résoudre l’équation de mouvement. Connaissant le champ de température à l'instant (𝑛 +
1)∆𝑡, l'étape suivante consiste à déterminer le champ de fonction de courant à partir de
l'équation (2.23). Les dérivées partielles du premier et du deuxième ordre sont approximées
selon le schéma centré classique. Ainsi au nœud (i, j) elles s’écrivent :
𝑛+1 𝑛+1 𝑛+1 𝑛+1
𝜕𝑇 𝑇𝑖+1,𝑗 − 𝑇𝑖−1,𝑗 𝜕𝑇 𝑇𝑖,𝑗+1 − 𝑇𝑖,𝑗−1
) = ) =
𝜕𝑥 𝑖,𝑗 2∆𝑥 𝜕𝑧 𝑖,𝑗 2∆𝑧

𝜕𝜇𝑎 𝜇𝑎 𝑛+1
𝑖+1,𝑗
− 𝜇𝑎 𝑛+1
𝑖−1,𝑗 𝜕𝜇𝑎 𝜇𝑎 𝑛+1
𝑖,𝑗+1
− 𝜇𝑎 𝑛+1
𝑖,𝑗−1
) = ) = (A.10)
𝜕𝑥 𝑖,𝑗 2∆𝑥 𝜕𝑧 𝑖,𝑗 2∆𝑧

𝑛+1 𝑛+1 𝑛+1 𝑛+1 𝑛+1 𝑛+1

𝜕2𝜓 𝜓𝑖+1,𝑗 − 2𝜓𝑖,𝑗 + 𝜓𝑖−1,𝑗 𝜕2𝜓 𝜓𝑖,𝑗+1 − 2𝜓𝑖,𝑗 + 𝜓𝑖,𝑗−1
) = ) =
𝜕𝑥 2 𝑖,𝑗 ∆𝑥 2 𝜕𝑧 2 𝑖,𝑗 ∆𝑧 2

En substituant les approximations approchées (A.10) dans l'équation (2.23) et en utilisant la

méthode itérative de sur-relaxation SOR, on obtient l’équation algébrique de la fonction de
courant au nœud (i, j):

𝑘 𝑘+1 𝑘 𝑘+1
𝜇𝑎 𝑖,𝑗 (𝑎1 (𝜓𝑖,𝑗+1 + 𝜓𝑖,𝑗−1 ) + 𝑎2 (𝜓𝑖+1,𝑗 + 𝜓𝑖−1,𝑗 ))
𝑘+1 𝑘 𝜔
𝜓𝑖,𝑗 = (1 − 𝜔)𝜓𝑖,𝑗 + +𝑎3 𝑢𝑖,𝑗 (𝜇𝑎 𝑖,𝑗+1 − 𝜇𝑎 𝑖,𝑗−1 ) + 𝑎4 𝑣𝑖,𝑗 (𝜇𝑎 𝑖+1,𝑗 (A.11)
𝜇𝑎 𝑖,𝑗
[ −𝜇𝑎 𝑖−1,𝑗 ) + 𝑎5 (𝑇𝑖+1,𝑗 − 𝑇𝑖−1,𝑗 ) ]

avec

130
 𝑏∆𝑥² ∆𝑧² ∆𝑥²∆𝑧
𝑎1 = 𝑎2 = 𝑎3 =
2(𝑏∆𝑥² + ∆𝑧²) 2(𝑏∆𝑥² + ∆𝑧²) 4(𝑏∆𝑥² + ∆𝑧²)
(A.12)
∆𝑥∆𝑧² ℜ∆𝑥∆𝑧² ℋ²
𝑎4 = 𝑎5 = 𝑏 =1+
4(𝑏∆𝑥² + ∆𝑧²) 4(𝑏∆𝑥² + ∆𝑧²) 𝜇𝑎

où 𝜔 est le coefficient de sur-relaxation, il est introduit pour accélérer la convergence et donc

pour diminuer le nombre d'itérations. Sa valeur optimale, 𝜔𝑜𝑝𝑡 dépend de la géométrie, du
maillage et des conditions aux limites. Pour un domaine rectangulaire avec un maillage
uniforme, 𝜔𝑜𝑝𝑡 étant (Roache [82]):

(1 − √1 − 𝛿)
𝜔𝑜𝑝𝑡 = 2 (A.13)
𝛿

où
2
𝜋 ∆𝑥 2 𝜋
2𝑐𝑜𝑠 ( ) + 2 𝑐𝑜𝑠 ( )
𝑛𝑥 ∆𝑧 𝑛𝑧
𝛿=( ) (A.14)
∆𝑥 2
(1 + 2 )
∆𝑧

Les expressions des approximations approchées des composantes de la vitesse 𝑢𝑖,𝑗 , 𝑤𝑖,𝑗 et le
champs de la viscosité apparente 𝜇𝑎 𝑖,𝑗 sont respectivement:

𝑛+1 𝑛+1
𝜓𝑖,𝑗+1 − 𝜓𝑖,𝑗−1
𝑢𝑖,𝑗 = (A.15)
2∆𝑧

𝑛+1 𝑛+1
𝜓𝑖+1,𝑗 − 𝜓𝑖−1,𝑗
𝑤𝑖,𝑗 = (A.16)
2∆𝑥

𝑛−1
2
𝜇𝑎 𝑖,𝑗 = (𝑢𝑖,𝑗 2
+ 𝑤𝑖,𝑗 ) 2 (A.17)

On note que l’indice 𝑛 dans l’équation (A.17) est l’indice de la loi en puissance et pas un
indice de discrétisation.

A.2 Méthode des volumes finis

Nous utilisons la méthode des volumes finis pour résoudre numériquement les équations
gouvernantes. Contrairement à la méthode des différences finies qui met en jeu des
approximations des dérivées, la méthode des volumes finis exploite des approximations

131
d'intégrales. L’avantage est que les propriétés de conservation sont, par principe, beaucoup
mieux assurées avec une méthode de volumes finis. Les équations différentielles (2.22)-(2.24)
régissant la convection naturelle d’un fluide de type loi en puissance dans un milieu poreux
sont transformées en un système d'équations algébriques en utilisant l'approche du volume de
contrôle.

A.2.1 Résumé de la démarche numérique

La technique numérique utilisée dans ce travail est caractérisée par ces points principaux:

 Les équations sont discrétisées selon la méthode des volumes finis, se basant sur le
schéma d’interpolation QUICK du troisième ordre.
 L’algorithme de résolution est basé sur le coupleur des équations de pression SIMPLE
[83-84].
 La technique du maillage utilisée est la technique de multi-grille pour obtenir des
résultats plus précis. Cette technique offre l’avantage d’un gain considérable en temps
de calcul.
 Les corrections différentielles sont effectuées pour les valeurs aux faces des volumes
de contrôle.

A.2.2 Discrétisation des équations

Chacune des équations de conservation peut s’écrire sous la forme générale suivante :

𝜕 𝜕 𝜕 𝜕
(𝜏𝜙 𝜙) + (𝑢𝐶𝜙 𝜙) + (𝑣𝐶𝜙 𝜙) + (𝑤𝐶𝜙 𝜙)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧
(A.18)
𝜕 𝜕𝜙 𝜕 𝜕𝜙 𝜕 𝜕𝜙
= (Γ𝜙 ) + (Γ𝜙 ) + (Γ𝜙 ) + 𝑆𝜙
𝜕𝑥 𝜕𝑥 𝜕𝑦 𝜕𝑦 𝜕𝑧 𝜕𝑧

Où 𝜙 désigne la variable à calculer, 𝐶𝜙 est le coefficient d’advection, Γ𝜙 est le coefficient de

diffusion et 𝑆𝜙 est le terme source.

L’écriture des équations sous la forme (A.18) permet une grande flexibilité de
programmation, car elle conduit à la résolution d’une forme unique et généralisée d’équations.

132
 Figure A.2 : Volume élémentaire de contrôle.

Le tableau A.1 définit les quantités correspondant aux différentes équations de conservation.
Ces équations non linéaires et fortement couplées seront résolues numériquement en utilisant
l’algorithme SIMPLE pour le couplage pression-vitesse. La méthode des volumes finis est
utilisée pour exprimer les intégrales sur les volumes finis [83]. Cette méthode s’appuie sur
une discrétisation du domaine de calcul en différents nœuds, chacun d’entre eux étant entouré
d’un volume élémentaire (Figure A.2) sur lequel on recherche la valeur moyenne de la
variable. Les dérivées partielles étant linéarisées en utilisant le théorème de TAYLOR-
YOUNG.

Equation : 𝜙 𝜏𝜙 𝐶𝜙 Γ𝜙 𝑆𝜙

Continuité 𝑉 0 1 0 0
𝜕𝑃
Suivant x 𝑢 0 1 0 − − 𝑢(|𝑉|𝑛−1 + ℋ²)
𝜕𝑥
Quantité de
mouvement

𝜕𝑃
Suivant y 𝑣 0 1 0 − − 𝑣(|𝑉|𝑛−1 + ℋ²)
𝜕𝑦
𝜕𝑃
Suivant z 𝑤 0 1 0 − − 𝑤|𝑉|𝑛−1 + ℜ𝑇
𝜕𝑧
Energie 𝑇 1 1 1 0

Tableau A.1 : Identification des termes sources et des coefficients pour les équations de
transport sous la forme adimensionnelle.

133
Nous pouvons réduire l’équation (A.18) sous la forme :

𝜕𝜙 𝜕𝐽𝑥 𝜕𝐽𝑦 𝜕𝐽𝑧

+ + + = 𝑆𝜙 (A.19)
𝜕𝑡 𝜕𝑥 𝜕𝑦 𝜕𝑧

avec

𝜕𝜙 𝜕𝜙 𝜕𝜙
𝐽𝑥 = 𝑢𝜙 − Γ𝜙 ; 𝐽𝑦 = 𝑣𝜙 − Γ𝜙 ; 𝐽𝑧 = 𝑤𝜙 − Γ𝜙 ; 𝐽 = ∑ 𝐽𝑖 ; 𝑖 = 𝑥, 𝑦, 𝑧 (A.20)
𝜕𝑥 𝜕𝑦 𝜕𝑧
𝑖

En intégrant l’équation (A.19) à travers le volume de contrôle de la Figure A.2, on aboutit à la

somme discrète suivante :

𝜕𝜙
∫ 𝑑Ω + ∫ ∇𝐽𝑑Ω = ∫ S𝑑Ω (A.21)
Ω 𝜕𝑡 Ω Ω

𝜙𝑝𝑛+1 + 𝜙𝑝𝑛
∆𝑉 + (𝐽𝑒 − 𝐽𝑤 )∆𝑧∆𝑦 + (𝐽𝑛 − 𝐽𝑠 )∆𝑧∆𝑥 + (𝐽𝑡 − 𝐽𝑏 )∆𝑥∆𝑦 = 𝑆𝜙 ∆𝑉 (A.22)
∆𝜏

Les quantités 𝐽𝑖 avec 𝑖 = 𝑒, 𝑤, 𝑛, 𝑠, 𝑡, 𝑏 représentent les flux totaux sur les faces des volumes
de contrôles et ∆𝑉 = ∆𝑥∆𝑦∆𝑧 est le volume de contrôle total. Lorsqu’on intègre l’équation de
conservation de la masse sur un volume de contrôle on aboutit à une égalité de la forme :

𝐹𝑒 − 𝐹𝑤 + 𝐹𝑛 − 𝐹𝑠 + 𝐹𝑡 − 𝐹𝑏 = 0 (A.23)

Où 𝐹𝑖 représentent les débits massiques à travers les faces de la maille de contrôle.

𝐹𝑒 𝐹𝑤 𝐹𝑛 𝐹𝑠 𝐹𝑡 𝐹𝑏
𝑢𝑒 ∆𝑦∆𝑧 𝑢𝑤 ∆𝑦∆𝑧 𝑣𝑛 ∆𝑥∆𝑧 𝑣𝑠 ∆𝑥∆𝑧 𝑤𝑡 ∆𝑥∆𝑦 𝑤𝑏 ∆𝑥∆𝑦

On multiplie l’équation (A.23) par 𝜙𝑝 et on soustrait le produit de l’équation (A.22) :

𝜙𝑝𝑛+1 + 𝜙𝑝𝑛
∆𝑉 + (𝐽𝑒 ∆𝑦∆𝑧 − 𝐹𝑒 𝜙𝑝 ) + (𝐽𝑤 ∆𝑦∆𝑧 − 𝐹𝑤 𝜙𝑝 ) + (𝐽𝑛 ∆𝑥∆𝑧 − 𝐹𝑛 𝜙𝑝 )
∆𝜏
(A.24)
+ (𝐽𝑠 ∆𝑥∆𝑧 − 𝐹𝑠 𝜙𝑝 ) + (𝐽𝑡 ∆𝑦∆𝑥 − 𝐹𝑡 𝜙𝑝 ) + (𝐽𝑏 ∆𝑦∆𝑥 − 𝐹𝑏 𝜙𝑝 )
= 𝑆𝜙 ∆𝑉

134
Patankar [83] propose d’écrire les quantités 𝐽𝑖 ∆𝑆 − 𝐹𝑖 𝜙𝑝 sous la forme :

𝐽𝑖 ∆𝑆 − 𝐹𝑖 𝜙𝑝𝑖 = 𝑎𝐼 (𝜙𝑝 − 𝜙𝐼 ) (A.25)

avec 𝑖 = 𝑒, 𝑤, 𝑛, 𝑠, 𝑡, 𝑏 et 𝐼 = 𝐸, 𝑊, 𝑁, 𝑆, 𝑇, 𝐵.

En substituant la quantité (A.25) dans l’équation (A.24), nous obtenons une relation
algébrique entre les quantités de la variable au centre de la maille et aux nœuds adjacents aux
faces du volume de contrôle :

𝑎𝑝 𝜙𝑝 = 𝑎𝐸 𝜙𝐸 + 𝑎𝑊 𝜙𝑊 + 𝑎𝑁 𝜙𝑁 + 𝑎𝑆 𝜙𝑆 + 𝑎 𝑇 𝜙 𝑇 + 𝑎𝐵 𝜙𝐵 + 𝐵 (A.26)

avec

𝐵 = 𝑆𝜙 Δ𝑉 + 𝑎𝑝𝑛−1 𝜙𝑝𝑛−1 où 𝑎𝑝𝑛−1 = Δ𝑉 ⁄Δ𝑡

Les coefficients d’approximation 𝑎𝐸 , 𝑎𝑊 , 𝑎𝑁 , 𝑎𝑆 , 𝑎 𝑇 et 𝑎𝐵 sont définis en fonction de

l’interpolation de la valeur 𝜙 entre deux points voisins du maillage, 𝐵 est un groupement qui
contient le terme source et le terme 𝑎𝑝 𝜙𝑝 est la quantité calculée au pas précédent.

(𝐽𝑒 ∆𝑦∆𝑧 − 𝐹𝑒 𝜙𝑝 ) = 𝑎𝐸 (𝜙𝑝 − 𝜙𝐸 ) 𝑎𝐸 = 𝐷𝑒 . 𝐴(|𝒫𝑒 |) + ‖−𝐹𝑒 , 0‖

(𝐽𝑤 ∆𝑦∆𝑧 − 𝐹𝑤 𝜙𝑝 ) = 𝑎𝑊 (𝜙𝑝 − 𝜙𝑊 ) 𝑎𝑊 = 𝐷𝑤 . 𝐴(|𝒫𝑤 |) + ‖𝐹𝑤 , 0‖

(𝐽𝑛 ∆𝑥∆𝑧 − 𝐹𝑛 𝜙𝑝 ) = 𝑎𝑁 (𝜙𝑝 − 𝜙𝑁 ) 𝑎𝑁 = 𝐷𝑛 . 𝐴(|𝒫𝑛 |) + ‖−𝐹𝑛 , 0‖

(𝐽𝑠 ∆𝑥∆𝑧 − 𝐹𝑠 𝜙𝑝 ) = 𝑎𝑆 (𝜙𝑝 − 𝜙𝑆 ) 𝑎𝑆 = 𝐷𝑠 . 𝐴(|𝒫𝑠 |) + ‖𝐹𝑠 , 0‖

(𝐽𝑡 ∆𝑦∆𝑥 − 𝐹𝑡 𝜙𝑝 ) = 𝑎 𝑇 (𝜙𝑝 − 𝜙 𝑇 ) 𝑎 𝑇 = 𝐷𝑡 . 𝐴(|𝒫𝑡 |) + ‖−𝐹𝑡 , 0‖

(𝐽𝑏 ∆𝑦∆𝑥 − 𝐹𝑏 𝜙𝑝 ) = 𝑎𝐵 (𝜙𝑝 − 𝜙𝐵 ) 𝑎𝐵 = 𝐷𝑏 . 𝐴(|𝒫𝑏 |) + ‖𝐹𝑏 , 0‖

Tableau A.2 : Expressions des coefficients 𝑎𝐸 , 𝑎𝑊 , 𝑎𝑁 , 𝑎𝑆 , 𝑎 𝑇 et 𝑎𝐵

𝐴(|𝒫𝑖 |) est la fonction d’interpolation et 𝒫𝑖 est le nombre de Péclet de maille. Les paramètres
𝐷𝑖 sont les coefficients de diffusion, appelés également conductances, et sont donnés par les
expressions suivantes :

135
 𝐷𝑒 𝐷𝑤 𝐷𝑛 𝐷𝑠 𝐷𝑡 𝐷𝑏
Γ𝑒 Δ𝑥Δ𝑧 Γ𝑤 Δ𝑥Δ𝑧 Γ𝑛 Δ𝑦Δ𝑧 Γ𝑠 Δ𝑦Δ𝑧 Γ𝑡 Δ𝑥Δy Γ𝑏 Δ𝑥Δy
𝛿𝑥)𝑒 𝛿𝑥)𝑤 𝛿𝑥)𝑛 𝛿𝑥)𝑠 𝛿𝑥)𝑡 𝛿𝑥)𝑏

Le nombre de Péclet de maille est défini comme suit :

𝐹𝑖 𝑐𝑜𝑛𝑣𝑒𝑐𝑡𝑖𝑜𝑛
𝒫𝑖 = = , 𝑖 = 𝑒, 𝑤, 𝑛, 𝑠, 𝑡, 𝑏 (A.27)
𝐷𝑖 𝑐𝑜𝑛𝑑𝑢𝑐𝑡𝑖𝑜𝑛

𝒫𝑒 𝒫𝑤 𝒫𝑛 𝒫𝑠 𝒫𝑡 𝒫𝑏
u𝑒 𝛿𝑥)𝑒 u𝑤 𝛿𝑥)𝑤 u𝑛 𝛿𝑥)𝑛 u𝑠 𝛿𝑥)𝑠 u𝑡 𝛿𝑥)𝑡 u𝑏 𝛿𝑥)𝑏
Γ𝑒 Γ𝑤 Γ𝑛 Γ𝑠 Γ𝑡 Γ𝑏
Tableau A.3 : Nombre de Péclet aux différents nœuds.

A.2.3 Accélération de la convergence

Les essais préliminaires, effectués pour tester la sensibilité des résultats obtenus au choix du
maillage, ont montré qu’un maillage de 82 × 42 × 42 (= 144648 ≅ 1.5 × 105 points) est
suffisant pour modéliser correctement l’écoulement du fluide et le transfert de chaleur pour le
cas d’une cavité tridimensionnelle. Ainsi, un tel maillage engendre une amplification
d’opérations logiques dans les machines de simulation et donc un temps de calcul très grand
par rapport à celui des calculs bidimensionnels. Pour cette raison, nous utilisons la technique
multi-grille qui rend le temps de calcul comparables ou même inférieurs à ceux des codes 2D.

Cette technique, due à Hackbusch [85,86], Horthmann [87] et Wesseling [88], consiste à
accélérer la convergence vers la solution en passant par interpolation d’une solution calculée
par un maillage grossier vers une autre calculée par un maillage plus fin.

136
 Annexe B
B.1 Principe variationnel
L’application du principe variationnel pour le problème de stabilité thermique a été
prouvé par Pellew & Southwell [89] dans le cas des fluides délimités par deux plans parallèles
horizontaux infinis. Chandrasekhar [32] a appliqué le principe variationnel pour la même
géométrie lorsque le fluide est un conducteur électrique et en présence d'un champ
magnétique externe, uniforme, orienté arbitrairement et appliqué dans la direction verticale.

Nous montrerons maintenant que la solution de la valeur caractéristique du problème formulé

dans le système (2.46) peut s’exprimer en adoptant le principe variationnel. En effet, en
appliquant le principe variationnel, on peut écrire le système (2.46) sous forme matricielle et
le résoudre numériquement en utilisant les développements (3.1) et (4.1) jusqu’à un ordre N,
suffisamment grand pour assurer la convergence vers la valeur exacte du nombre de Rayleigh
critique. Un programme Fortran a été développé dans ce sens.

B.1.1 Discrétisation de la formulation variationnelle

La discrétisation des équations du système (2.46) consiste à faire le changement des variables
selon un maillage M*N, avec M et N sont respectivement le nombre des lignes et le nombre
des colonnes de la matrice.

Pour simplifier la résolution du problème, nous considérons 𝑧 et 𝑥 compris 0 et 1.

Les solutions (3.1) et (4.1) seront écrites de la manière discrète suivante :

𝐹𝑀 (𝑧) = ∑ 𝑊𝑀 sin(𝑀𝜋𝑧) (B.1)

𝑀

∑ 𝜃𝑀 cos((𝑀 − 1)𝜋𝑧) 𝐹𝑙𝑢𝑥 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡

𝑀
𝐺𝑀 (𝑧) = (B.2)
∑ 𝜃𝑀 sin(𝑀𝜋𝑧) 𝑇𝑒𝑚𝑝é𝑟𝑎𝑡𝑢𝑟𝑒 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒
{ 𝑀

137
On définit les formes linaires par :

𝐿11 (𝑀) = −(𝑆 2 . 𝑀2 . 𝜋 2 + 𝑘 2 )

𝐿12 (𝑀) = 𝑘 2 𝒫1−𝑛
𝐿21 (𝑀) = −1 (B.3)
−((𝑀 − 1)2 . 𝜋 2 + 𝑘 2 ) 𝐹𝑙𝑢𝑥 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡
𝐿22 (𝑀) = { 2 2 2)
{ −(𝑀 . 𝜋 + 𝑘 𝑇𝑒𝑚𝑝é𝑟𝑎𝑡𝑢𝑟𝑒 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒

En substituant ces expressions ainsi que les solutions (B.1) et (B.2) dans le système (2.46),
nous obtenons,

𝐿 (𝑀)𝐹𝑀 (𝑧) + 𝑅𝐿12 (𝑀)𝐺𝑀 (𝑧) = 0

{ 11 (B.4)
𝐿21 (𝑀)𝐹𝑀 (𝑧) + 𝐿22 (𝑀)𝐺𝑀 (𝑧) = 0

La construction d’un système matricielle se base essentiellement sur la méthode de Galerkin,

on multiplie les équations par une fonction test qui vérifie les mêmes conditions aux limites
que les solutions (B.1) et (B.2), après on intègre le long de z. Les fonctions tests seront écrites
des manières discrètes suivantes:

𝐹𝑁 (𝑧) = 𝑠𝑖𝑛(𝑁𝜋𝑧)
{ 𝑐𝑜𝑠((𝑁 − 1)𝜋𝑧) 𝐹𝑙𝑢𝑥 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡 (B.5)
𝐺𝑁 (𝑧) = {
𝑠𝑖𝑛(𝑁𝜋𝑧) 𝑇𝑒𝑚𝑝é𝑟𝑎𝑡𝑢𝑟𝑒 𝑐𝑜𝑛𝑠𝑡𝑎𝑛𝑡𝑒

On multiplie les équations du système (B.4) par 𝐹𝑁 (𝑧) quand la variable est 𝑤 et par 𝐺𝑁 (𝑧)
quand le variable est 𝜃

Le système (B.4) devient,

1 1

∫ 𝐿11 (𝑀)𝐹𝑀 (𝑧)𝐹𝑁 (𝑧) 𝑑𝑧 = −𝑅 ∫ 𝐿12 (𝑀)𝐺𝑀 (𝑧)𝐹𝑁 (𝑧)𝑑𝑧

0 0
1 1 (B.6)
∫ 𝐿21 (𝑀)𝐹𝑀 (𝑧)𝐺𝑁 (𝑧) 𝑑𝑧 = − ∫ 𝐿22 (𝑀)𝐺𝑀 (𝑧)𝐺𝑁 (𝑧)𝑑𝑧
{0 0

En alternant 𝑀 et 𝑁 dans l’équation (B.6), on peut déduire que :

138
 1 1

∫ 𝐹𝑀 (𝑧)𝐹𝑁 (𝑧) 𝑑𝑧 = ∫ 𝐺𝑀 (𝑧)𝐺𝑁 (𝑧)𝑑𝑧 = 0 ; 𝑀≠𝑁 (B.7)

0 0

Les fonctions 𝐹𝐽 (𝑧) et 𝐺𝐽 (𝑧) sont orthogonales.

B.1.2 Matrices d’assemblage

La matrice globale qui regroupe toutes les matrices élémentaires est construite de la manière
suivante,

[𝑊1 ]
[𝜃1 ]

(B.8)

[𝑊2 ] [𝜃2 ]

En isolant les blocs de la matrice globale, les matrices élémentaires apparaissent et elles sont
définies par,

[𝑊1 ] = 𝐿11 (𝑀) ∫ 𝐹𝑀 (𝑧)𝐹𝑁 (𝑧) 𝑑𝑧

0
1

[𝜃1 ] = 𝐿21 (𝑀) ∫ 𝐹𝑀 (𝑧)𝐺𝑁 (𝑧)𝑑𝑧

0
1 (B.9)
[𝑊2 ] = 𝐿12 (𝑀) ∫ 𝐺𝑀 (𝑧)𝐹𝑁 (𝑧)𝑑𝑧
0
1

[𝜃2 ] = 𝐿22 (𝑀) ∫ 𝐺𝑀 (𝑧)𝐺𝑁 (𝑧)𝑑𝑧

{ 0

139
𝐹𝐽 et 𝐺𝐽 doivent être orthogonales.

Exemple d’une matrice globale à l’ordre 4 pour le cas de condition du flux constant imposée :

𝑘² 𝜋² 2𝑘² 2𝑘²
− − 0 0 0 ℜ 0 ℜ 0 𝑊1 0
2 2 𝜋 3𝜋
2
𝑘 4𝑘² 4𝑘²
0 − − 2𝜋² 0 0 0 ℜ 0 ℜ 𝑊2 0
2 3𝜋 5𝜋
𝑘² 9𝜋² 2𝑘² 6𝑘²
0 0 − − 0 ℜ 0 ℜ 0 𝑊3 0
2 2 3𝜋 5𝜋
𝑘2 8𝑘² 8𝑘² 𝑊4 0
0 0 0 − − 8𝜋² 0 ℜ 0 ℜ
2 15𝜋 7𝜋
. = (B.10)
2 2
− 0 − 0 𝑘2 0 0 0 𝜃1 0
𝜋 3𝜋
4 8 𝑘² 𝜋²
0 − 0 − 0 + 0 0 𝜃2 0
3𝜋 15𝜋 2 2
2 6 𝑘² 𝜃3
0 − 0 0 0 + 2𝜋² 0 0
3𝜋 5𝜋 2
4 8 𝑘² 9𝜋² [ 𝜃4 ] [0]
[ 0 0 − 0 0 0 +
5𝜋 7𝜋 2 2 ]

Exemple d’une matrice globale à l’ordre 4 pour le cas de condition de température constante
imposée :

𝑘² 𝜋² 𝑘²
− − 0 0 0 ℜ 0 0 0 𝑊1 0
2 2 2
2
𝑘 𝑘²
0 − − 2𝜋² 0 0 0 ℜ 0 0 𝑊2 0
2 2
𝑘² 9𝜋² 𝑘²
0 0 − − 0 0 0 ℜ 0 𝑊3 0
2 2 2
𝑘2 𝑘² 𝑊4 0
0 0 0 − − 8𝜋² 0 0 0 ℜ
2 2
. = (B.11)
1 𝑘² 𝜋²
−
2
0 0 0
2
+
2
0 0 0 𝜃1 0
1 𝑘²
0 − 0 0 0 + 2𝜋² 0 0 𝜃2 0
2 2
1 𝑘² 9𝜋² 𝜃3 0
0 0 − 0 0 0 + 0
2 2 2
1 𝑘² [ 𝜃4 ] [0]
[ 0 0 0 − 0 0 0 + 8𝜋²]
2 2

140
B.1.3 Détermination du nombre Rayleigh critique

Pour déterminer le nombre de Rayleigh critique, une suite d’opérations est nécessaire pour
ramener notre système matriciel vers un problème de recherche de valeurs propres.
L’instabilité se produit à partir de la valeur minimale de ces nombres obtenus.

Partant du système matriciel (B.8) on peut écrire :

[𝑊1 ](𝑊) = −ℜ[𝜃1 ](𝜃) (B.12)

[𝑊2 ](𝑊) = −[𝜃2 ](𝜃) (B.13)
A partir de (B.13) on peut écrire :

[𝜃2 ]−1 [𝑊2 ](𝑊) = −𝐼𝑑 (𝜃) (B.14)

Où 𝐼𝑑 est la matrice identité.

Remplaçant l’équation (B.14) dans (B.12)

[𝑊1 ](𝑊) = ℜ[𝜃1 ][𝜃2 ]−1 [𝑊2 ](𝑊) (B.15)

1
⇒ [ [𝑊1 ]−1 [𝜃1 ][𝜃2 ]−1 [𝑊2 ] − 𝐼𝑑 ] (𝑊) = (0) (B.16)
ℜ

Posons [𝐸] = [𝑊1 ]−1 [𝜃1 ][𝜃2 ]−1 [𝑊2 ]

L’écriture symbolique du résultat sera :

1 (B.17)
[[𝐸] − 𝐼 ] (𝑊) = (0)
ℜ 𝑑
Le problème est transformé par réduction d’un endomorphisme [𝐸].

Physiquement, l’équation (B.17) admet une solution non triviale si et seulement si le

déterminant |[𝐸] − 𝜆𝐼𝑑 | est nul. Ceci nous conduit à chercher toutes les valeurs propres 𝜆𝑖
vérifiant l’équation (B.17).

1
Où ℜ𝑐 = est la plus petite valeur qui génère la perte de stabilité de l’écoulement et 𝜆𝑚𝑎𝑥
𝜆𝑚𝑎𝑥

est le rayon spectral.

141