Sultan, Romain Jamal Arthur

Justine Mathieu

Jamel Simon

Fatiha Enfin

Au cours de ma thèse, j'ai été amené à côtoyer de nombreuses personnes qui m'ont accueilli, guidé et aidé. Sans elles, cette thèse n'aurait tout simplement pas pu voir le jour. A travers ces quelques lignes de remerciements, j'espère parvenir à leur témoigner toute ma gratitude et ma reconnaissance.

Mes remerciements vont

en premier lieu à mon directeur de thèse Chérif Nouar et à ma codirectrice de thèse Christel Métivier qui m'ont proposé un sujet intéressant, ambitieux et stimulant. Tout au long de ces quatre années de thèse, j'ai pu compter sur leurs conseils avisés, leur rigueur scientifique et leur grande disponibilité. Je les remercie également pour leur prévenance et leur soutien sans faille dans les périodes de doute. J'exprime ma profonde gratitude à Daniel Henry et Marc Médale pour avoir accepté d'être les rapporteurs de ce travail de thèse. Leurs remarques judicieuses et leurs questions pertinentes ont permis d'améliorer sensiblement ce manuscrit. Je remercie chaleureusement Alain Bergeon de m'avoir fait l'honneur de présider mon jury de soutenance de thèse. Merci également à Anne Davaille pour avoir examiné attentivement mon travail. Je remercie Célien, mon collègue mais avant tout mon ami, qui a partagé mon bureau pendant toutes ces années et qui a toujours fait preuve à mon égard d'une extrême gentillesse et d'une grande bienveillance. Je le remercie pour toutes les discussions passionnées autour de la musique classique et pour toutes les sorties au cinéma que nous avons faites. Je remercie Waseem, un autre collègue et ami qui m'a accueilli dans son bureau, pour sa bonne humeur au quotidien et sa joie de vivre communicative. Merci à lui également pour toutes les conversations scientifiques que nous avons eues. Toute ma gratitude va à Benoît, le dernier arrivant du bureau -006J qui a immédiatement trouvé sa place et qui a fait preuve dès son arrivée d'une grande bonté d'âme et d'une profonde sollicitude à l'égard de ses collègues. Je remercie bien évidemment mon ami Alexandre qui m'a toujours impressionné par la diversité de ses connaissances. Les conversations que nous avons eues se sont toujours révélées intéressantes et stimulantes. Bravo et merci également pour son engagement auprès des doctorants. Je remercie aussi Morgan pour son sens de l'humour aiguisé, ses mots d'esprit, sa connaissance encyclopédique des épisodes de ≪ Kaamelott ≫ et son amabilité. Les dis-cussions autour de la méthode de Monte Carlo utilisée dans sa thèse étaient passionnantes. Je remercie Naima pour sa jovialité, son caractère enjoué et sa bonne humeur. C'était un plaisir d'être son voisin de bureau.

Table des matières

TABLE DES MATI ÈRES

Introduction générale

Un fluide peut se mettre en mouvement sous l'effet d'un gradient thermique parallèle à la force de gravité. Ce phénomène physique, appelé convection thermique, est observé dans des contextes très divers et se manifeste à des échelles extrêmement variées. A l'échelle planétaire, la convection thermique est par exemple à l'origine des mouvements internes de l'atmosphère et entre en jeu dans de nombreux phénomènes météorologiques comme les orages, les cyclones ou encore les moussons. Elle est également la source des différents courants marins qui parcourent les océans de la planète et qui régulent les climats terrestres. A l'échelle des temps géologiques, l'accroissement de la température de la Terre avec la profondeur est responsable de la convection du manteau terrestre, l'élément moteur de la tectonique des plaques. Enfin, la convection thermique intervient également dans un grand nombre de procédés industriels et ce jusqu'à des échelles millimétriques comme dans certains dispositifs de refroidissement de composants électroniques par exemple. Du fait du nombre conséquent de situations où elle joue un rôle majeur, la convection thermique a été très étudiée depuis le début du XXème siècle. La configuration la plus utilisée est celle de Rayleigh-Bénard qui fait intervenir un fluide placé dans le champ de gravité ĝ et confiné entre deux plaques horizontales. La plaque inférieure est chauffée alors que la plaque supérieure est refroidie.

Les premières études expérimentales visant à étudier la configuration précédente furent réalisées aux alentours de l'an 1900 par Henri Bénard (1874-1939) au cours de sa thèse ayant pour titre ≪ Les Tourbillons cellulaires dans une nappe liquide ≫ [1]. Dans ses expériences, la surface supérieure était libre et les fluides utilisés étaient la plupart du temps de la cire ou du blanc de baleine qui sont liquides entre 50˚C et 100˚C. Le choix de ces fluides vient du fait qu'ils sont peu volatils dans la gamme de températures précédente : leur utilisation permet d'obtenir des couches de fluide uniformes. Les épaisseurs de ces couches étaient de l'ordre du millimètre. Grâce à différentes méthodes optiques, Henri Bénard a été en mesure d'observer différentes structures périodiques comme des polygones, le plus souvent des hexagones (Figure 1), mais aussi des rouleaux de convection. Il a notamment pu mesurer précisément la longueur d'onde des structures qui émergeaient ainsi que les variations dans le profil de l'interface. Toutefois, il ne disposait pas d'interprétations physiques permettant d'expliquer ses observations. Figure 1 -Photographie issue des expériences de Henri Bénard tiré de [1] En 1916, dans un article [2] intitulé ≪ On the convective currents in a horizontal layer of fluid when the higher temperature is on the under side ≫, Lord Rayleigh explique les résultats de Henri Bénard par l'existence de forces de flottabilité induites par les variations de masse volumique dues au gradient de température. Selon lui, deux phénomènes physiques s'opposent :

-D'un côté, la poussée d'Archimède tend à faire s'élever les couches de fluide basses plus chaudes donc moins denses.

-De l'autre, les forces de frottement dues à la viscosité du fluide et la conduction thermique s'opposent à la mise en mouvement. L'explication physique fournie par Rayleigh a été acceptée pendant plus de 40 ans, jusqu'en 1958, où Pearson a démontré dans [3] le rôle prédominant de la tension superficielle dans les résultats expérimentaux de Bénard. En effet, des inhomogénéités Pour des couches de fluides de faible épaisseur, les forces induites par les gradients de tension superficielle sont prépondérantes et on parle de convection de Bénard-Marangoni.

En revanche, pour des couches de fluide de plus grande épaisseur, ce sont les forces de flottabilité qui dominent et on parle alors de convection de Rayleigh-Bénard. A titre indicatif, pour de l'eau pure à 20 ˚C, l'épaisseur en dessous de laquelle l'instabilité de Bénard-Marangoni l'emporte sur celle de Rayleigh-Bénard est de 2.6 cm d'après [4].

Théoriquement, Rayleigh est parvenu à déterminer, à l'aide d'une analyse linéaire de stabilité, les conditions critiques du système à savoir la valeur du nombre de Rayleigh Ra c au delà de laquelle survient l'instabilité ainsi que le nombre d'onde k c des structures qui émergent. Il considère que les plaques sont parfaitement conductrices de la chaleur, ont une extension horizontale infinie et suppose que les surfaces délimitant le fluide sont libres.

Les calculs sont menés dans le cadre des approximations de Boussinesq :

-la masse volumique ne varie que dans le terme de flottabilité de l'équation de conservation de la quantité de mouvement, -la dissipation visqueuse ainsi que la dépendance de l'énergie interne avec la pression sont négligées dans l'équation de l'énergie, -les coefficients thermodynamiques sont considérés constants et la masse volumique ne varie qu'avec la température (pas de compressibilité).

On donne ici les principales étapes de l'analyse linéaire de stabilité qui a permis à Rayleigh de déterminer les conditions critiques (k c , Ra c ).

La couche de fluide est décrite mathématiquement par les équations générales de conservations de la masse, de la quantité de mouvement et de l'énergie. Dans le cadre utilisé par Rayleigh, ces dernières s'écrivent :

∇. V = 0 (2) ρ0 ∂ V ∂ t + ∇ V. V = -∇ P + ρ( T )ĝ + μ0 ∆ V (3) ρ0 Ĉ0 ∂ T ∂ t + ∇ T . V = λ0 ∆ T (4)
où ρ( T ) = ρ0 (1 -β( T -T0 )).

Dans les équations précédentes, V, P et T désignent respectivement les champs de vitesse, de pression et de température. Ĉ0 et λ0 désignent respectivement la capacité calorifique et la conductivité thermique du fluide. La masse volumique ρ ne varie linéairement avec la température que dans le terme de poussée conformément à l'approximation de Boussinesq.

Dans l'état de base, le fluide est au repos : le champ de vitesse est donc nul soit Vb = 0 et les équations décrivant cet état conductif sont alors :

∆ Tb = 0 (5)

-∇ Pb + ρ0 (1 -β( Tb -T0 ))ĝ = 0 (6)
L'extension latérale de la couche de fluide étant supposée infinie, le champ de température dans l'état conductif ne dépend que de ẑ. L'intégration de l'équation (5)

Pour étudier la stabilité de la couche de fluide, on superpose aux grandeurs de l'état de base des perturbations en température θ, en pression p et en vitesse v soit :

T = Tb + θ P = Pb + p V = vb + v = v = ue x + we z (8)
En introduisant ces différentes expressions dans les équations (2), (3) et ( 4) et en retranchant ( 5), ( 6) et (7), on obtient un système d'équations aux perturbations. Ce dernier peut être écrit sous forme adimensionnelle en utilisant respectivement ∆ T , d, d2 /κ 0 , κ0 / d et μ0 κ0 / d2 comme échelles caractéristiques de température, de distance, de temps, de vitesse et de pression. Le système d'équations aux perturbations s'écrit alors :

∇.v = 0 (9)
1 Pr ∂v ∂t + ∇v.v = -∇p + Ra θe z + ∆v (10) ∂θ ∂t + ∇θ.v = w + ∆θ (11) Deux nombres adimensionnels entrent en jeu dans le système d'équations aux perturbations ( 9)-( 11) :

-le premier est le nombre de Rayleigh Ra = ρ0 ĝ β∆ T d3 μ0 κ0 déjà introduit précédemment.

-le second est le nombre de Prandtl Pr = μ0 ρ0 κ0 qui compare les temps de diffusion thermique et de diffusion visqueuse.

Dans l'analyse linéaire de stabilité, on considère des perturbations infinitésimales de la forme :

(u, v, w, p, θ) = δ(u 1 , v 1 , w 1 , p 1 , θ 1 ) + O(δ 2 ) où δ << 1

On introduit ces perturbations dans le système ( 9)- (11) et on ne conserve que les termes d'ordre δ ce qui conduit à un système d'équations linéarisées :

∂u 1 ∂x + ∂w 1 ∂z = 0 ( 12
)
1 Pr

∂u 1 ∂t = - ∂p 1 ∂x + ∆u 1 (13)
1 Pr

∂v 1 ∂t = - ∂p 1 ∂y + ∆v 1 (14)
1 Pr

∂w 1 ∂t = - ∂p 1 ∂z + Ra θ 1 + ∆w 1 (15)
∂θ 1 ∂t = w 1 + ∆θ 1 (16)
Ce système étant invariant par translation suivant x et y, on peut chercher des solutions sous la forme de modes normaux :

        u 1 v 1 w 1 p 1 θ 1         =         ũ1 (z) ṽ1 (z) w1 (z) p1 (z) θ1 (z)         exp(σt) exp(i(k x x + k y y)) (17)
Dans l'expression précédente, k = k x e x + k y e y est le vecteur d'onde et σ est le taux de croissance des perturbations. Quand Re(σ) > 0, les perturbations s'amplifient avec le temps et l'état conductif est donc instable. A l'inverse, quand Re(σ) < 0, les perturbations s'atténuent au cours du temps et l'état conductif est stable. En introduisant la décomposition sous forme de modes normaux (17) dans le système (12)-( 16) et en éliminant ũ1 , ṽ1 et p1 , on obtient un problème aux valeurs propres pour les variables w1 et θ1 :

Pr(D 2 -k 2 ) 2 w1 -Pr Ra k 2 θ1 = σ(D 2 -k 2 ) w1 (18)
w1 + (D 2k 2 ) θ1 = σ θ1 (19) où

k 2 = k 2 x + k 2 y et D = d dz .
Dans le cas de parois libres, Lord Rayleigh est parvenu à déterminer analytiquement la courbe de stabilité marginale associée à ( 18)- (19). Cette dernière correspond au lieu Ra(k) des points (k,Ra) pour lesquels σ = 0. Elle délimite les domaines stable et instable de l'état conductif : ce dernier est stable en dessous de cette courbe et instable au dessus.

Elle présente un minimum en deçà duquel l'état conductif est nécessairement stable : les coordonnées de ce point (k c , Ra c ) sont appelées conditions critiques. Les valeurs obtenues par Lord Rayleigh pour ces dernières sont :

k c = π √ 2 et Ra c = 27 4 π 4 ≈ 657.5
La convection ne peut ainsi se développer que pour des valeurs du nombre de Rayleigh Ra ≥ Ra c = 657.5 quand on dispose de parois libres.

Dans le cas plus réaliste de parois rigides, les équations ( 18)-( 19) n'admettent pas de solution analytique. Une résolution numérique a donc été menée par Pellew et Southwell dans [5] en 1940 afin de déterminer la courbe de stabilité marginale et les conditions critiques associées. Ces dernières ont été évaluées à :

k c = 3.11 et Ra c = 1707.76
Les courbes de stabilité marginales obtenues dans le cas de parois libres et de parois rigides sont données sur la Figure 3 avec les conditions critiques correspondantes. La configuration avec deux parois libres est intéressante parce qu'elle permet des développements analytiques. En revanche, sa réalisation s'avère difficile d'un point de vue expérimental. Toutefois, une manière de s'en approcher est donnée par Goldstein et Graham dans [6]. L'idée de ces auteurs est de considérer comme bonne approximation d'une surface libre, une interface entre deux fluides non miscibles de viscosités très différentes. Ces auteurs ont donc utilisé une huile de silicone de forte viscosité reposant sur une couche de mercure liquide et surmontée d'une couche d'hélium gazeux, toutes deux de viscosité beaucoup plus faible (respectivement d'un facteur 10 6 et 10 4 ). Du point de vue thermique, la forte conductivité du mercure permet de considérer que la température est imposée à l'interface huile-mercure. A l'interface huile-hélium, le flux de chaleur est quasiment nul du fait de la faible conductivité de l'hélium. Avec leur dispositif, Goldstein et Graham trouvent une valeur du nombre de Rayleigh critique Ra c = 596 ce qui représente une erreur de 10% avec la valeur théorique calculée par Lord Rayleigh.

D'autres conditions aux limites peuvent être envisagées. Par exemple, on peut considérer que les plaques limitant le fluide ne sont pas infiniment conductrices de la chaleur, ce qui est le cas dans la pratique. Des auteurs se sont ainsi intéressés à l'influence de la nature finie de la conductivité thermique des plaques sur les conditions critiques.

L'étude la plus ancienne a été menée par Hurle, Jakeman et Pike dans [7] en 1967.

Ces auteurs ont notamment montré que les conditions critiques Ra c et k c diminuaient d'autant plus que les plaques étaient isolantes.

Plus récemment, des études [8] et [9] se sont intéressées à l'impact d'un possible glissement survenant au niveau des parois délimitant le fluide. Ce phénonème est en effet susceptible de survenir dans des systèmes où le nombre de Knudsen est très faible, comme dans des gaz raréfiés par exemple. [8] et [9] ont montré que plus le phénomène de glissement était important et davantage les conditions critiques étaient réduites.

Dans tous les cas, l'analyse linéaire de stabilité permet de déterminer les conditions critiques (k c ,Ra c ) à partir desquelles se développent des motifs de convection mais ne donne pas d'indication sur la nature de la structure qui émerge au seuil de l'instabilité.

En effet, n'importe quel mode ou combinaison linéaire de modes de nombre d'onde k c est solution du problème linéarisé. Une analyse faiblement non linéaire de stabilité s'avère donc nécessaire afin d'étudier les motifs de convection.

Malkus et Veronis sont les premiers à proposer en 1958 dans [10] une méthode permettant de déterminer des solutions stationnaires aux équations descriptives du problème ( 9)-( 11) sous forme de motifs de convection. Leur approche repose sur un développement asymptotique des inconnues du problème : la vitesse, la perturbation en température et le nombre de Rayleigh sont ainsi écrits en série d'un petit paramètre δ ≪ 1 :

v = δv 1 + δ 2 v 2 + δ 3 v 3 + O(δ 4 ) (20) θ = δθ 1 + δ 2 θ 2 + δ 3 θ 3 + O(δ 4 ) (21)
Ra = Ra c + δRa 1 + δ 2 Ra 2 + O(δ 3 ) (22)
En introduisant ( 20)-( 22) dans ( 9)- (11), on aboutit à un système d'équations différentielles ordinaires pour v p et θ p à chaque ordre δ p pour p ≥ 1.

A l'ordre δ, on retrouve le problème linéaire qui admet une infinité de solutions de nombre d'onde k c . Différents motifs de convection peuvent alors être considérés.

Mathématiquement, les champs de vitesse v 1 et de température θ 1 peuvent être décrits comme différentes séries de Fourier selon la nature de la structure considérée :

-pour un rouleau, v 1 et θ 1 sont caractérisés par un seul mode.

-pour un carré, deux modes de Fourier avec des vecteurs d'ondes orthogonaux sont nécessaires.

-pour un hexagone, trois modes avec des vecteurs d'ondes inclinés de 2π/3 les uns par rapport aux autres doivent être utilisés.

dans l'approximation de Boussinesq et ce, que l'on considère adhérence ou glissement parfait aux parois.

En revanche, quand les gradients thermiques sont trop importants pour que les hypothèses de Boussinesq restent valables, d'autres motifs de convection peuvent être stables. Dans un article [13] paru en 1967, Busse montre comment la convection d'un fluide thermodépendant s'organise au fur et à mesure qu'on s'éloigne des conditions critiques.

A proximité du seuil de l'instabilité, les motifs hexagonaux sont les seules structures stables. Puis, quand la convection gagne en intensité, des rouleaux émergent aux côtés des hexagones donnant ainsi naissance à un domaine de bistabilité. Enfin, quand les écarts relatifs au seuil deviennent trop importants, les hexagones perdent leur stabilité et il ne subsiste que des rouleaux. Dans le cas de fluides très fortement thermodépants, des études menées par Jenkins [14] et Busse [15] ont montré que les carrés étaient les motifs de convection stables au détriment des rouleaux.

L'approche faiblement non linéaire mise en place par Malkus et Veronis et reprise par Busse est désormais très fréquemment utilisée dans de nombreux travaux. Dans les études les plus récentes, cette méthode est employée pour déterminer les équations d'amplitude du système. Plus précisément, ces équations sont obtenues à partir de l'application de l'alternative de Fredholm pour déterminer les inconnues Ra 1 , Ra 2 . Ce formalisme décrit la dynamique du système à proximité des conditions critiques et a l'avantage de permettre d'étudier simplement la stabilité des motifs. A titre d'exemple, l'équation d'amplitude pour un rouleau de convection est donnée à l'ordre cubique par :

dA dt = g 0 A -g 1 A|A| 2 (23)
où A est l'amplitude de la perturbation suivant l'axe du rouleau et g 0 , g 1 des constantes.

Les différents motifs de convection qui émergent au seuil de la convection vont ensuite perdre leur stabilité par le biais d'instabilités dites ≪ secondaires ≫. Pour comprendre l'origine de ces dernières, on peut effectuer un développement en série de Taylor du taux de croissance σ au voisinage du seuil de la convection suivant :

σ(k, Ra) = 1 τ 0 Ra -Ra c Ra c - ξ 2 0 τ 0 (k -k c ) 2 (24)
Dans l'expression précédente, τ 0 désigne le temps de relaxation des perturbations et ξ 0 est la longueur de cohérence. Ces deux grandeurs sont évaluées aux conditions critiques et sont définies par :

1 τ 0 = Ra c ∂σ ∂Ra ξ 2 0 τ 0 = - 1 2 ∂ 2 σ ∂k 2
En fonction des valeurs du paramètre de contrôle Ra, plusieurs situations, représentées sur la Figure 4, sont possibles :

-si Ra ≤ Ra c , tous les nombres d'onde ont un taux de croissance négatif et l'état conductif est stable.

-si Ra = Ra c , il existe un nombre d'onde k c pour lequel le taux de croissance est nul et on a bifurcation de l'état conductif vers un état convectif ordonné.

-si Ra ≥ Ra c , toute une bande de nombres d'onde centrée sur k c ont un taux de croissance positif. Les différentes instabilités secondaires délimitant le domaine de stabilité des rouleaux de convection ont d'abord été identifiées dans la limite des grands nombres de Prandtl par Busse d'abord théoriquement dans [16] puis expérimentalement dans [17]. Le domaine de stabilité des rouleaux est ainsi appelé ≪ Ballon de Busse ≫ : pour de fortes valeurs du nombre de Prandtl, il est délimité (cf Figure 3.4 ) par :

-l'instabilité d'Eckhaus.

-l'instabilité Zigzag.

-l'instabilité en rouleaux croisés.

-l'instabilité bimodale. L'instabilité en rouleaux croisés est une instabilité directement liée à l'intensité de la convection. Un ensemble de rouleaux instables est progressivement remplacé par un autre réseau de rouleaux qui se développe orthogonalement au détriment du premier. La Figure 8 illustre le développement de cette instabilité. Comme pour les instabilités précédentes,

20000• \ OE lOO<.•> \ \ 1 GUOLC ! SOC>> ï \' 151' ... \ V'<tf<Olt ><C? r 20C') ï G ^ G M M セ @ lOC•> 1 .
l O l t :I4 SS 7e Il s'agit des ≪ oscillations de Busse ≫ et de l'instabilité variqueuse oblique, toutes deux mises en évidence par Busse et ses collaborateurs. Dans [19] et [20], il est montré que la vorticité n'est pas nulle à faible nombre de Prandtl, ce qui engendre un ≪ écoulement de dérive ≫. Le couplage entre cet ≪ écoulement de dérive ≫ et les perturbations en vitesse et en température se révèle être à l'origine de ces deux instabilités. Dans la thèse, on s'est intéressé uniquement aux fluides à grands nombres de Prandtl, on ne détaille donc pas davantage la manière dont ces deux instabilités se manifestent ou l'allure de la section du ≪ Ballon de Busse ≫ obtenue à faible nombre de Prandtl. Le lecteur intéressé pourra trouver de plus amples informations dans [18] ou [21].

La stabilité des autres motifs de convection vis à vis d'une modulation spatiale a également été considéré dans un certain nombre de travaux. La stabilité des carrés a ainsi été étudiée par Hoyle dans [22] et par Holmedal dans [23]. En ce qui concerne les structures hexagonales, leur stabilité vis à vis d'une modulation spatiale a fait l'objet d'une série d'études dont les plus significatives sont dues à Sushchik et Tsimring [24],

Hoyle [25], Echebarria et Pérez-Garcia [26] et Kuske et Milewski [27]. Les principaux résultats de ces travaux sont présentés dans les deux premiers chapitres de la thèse, on ne les explicite donc pas ici. La première étude numérique du démarrage de la convection et de l'évolution du coefficient de transfert de chaleur pour un fluide rhéofluidifiant a été effectuée par [32]

N u N u newt = 0.87n 2 p -2.28n p + 2.41 (25)
qui peut être utilisée pour Ra ≤ Ra c ≤ 2 Ra c avec 0.5 ≤ n p < 1. Un accord qualitatif avec les résultats expérimentaux [29] est observé.

La décroissance du nombre de Rayleigh critique avec la diminution de l'indice de structure pour un fluide en loi puissance a été confirmée par des simulations numériques 2D

effectuées par [34] et [35]. Dans cette dernière référence, les auteurs montrent clairement que pour un nombre de Rayleigh donné, l'intensité de convection augmente avec le caractère rhéofluidifiant. C'est ce qui explique l'augmentation du coefficient de transfert de chaleur.

[36] a étudié numériquement, la structure de l'écoulement et du champ thermique en convection de Rayleigh-Bénard pour un fluide rhéofluidifiant avec un nombre de Prandtl très élevé de façon à négliger les termes d'inertie. Le comportement rhéologique du fluide est décrit par un modèle en loi puissance. Les résultats obtenus mettent en évidence des modifications fondamentales dans la structure de l'écoulement lorsque l'indice de structure n p est inférieur à 0.3. Ces modifications se caractérisent en particulier par le développement de zones de fluide stagnant avec une forte viscosité. Pour corréler l'augmentation du coefficient de transfert avec la diminution de l'indice de structure, Parmentier utilise dans [36] un nombre de Rayleigh basé sur une viscosité moyenne pondérée par le deuxième invariant du tenseur des taux de déformation :

μ = S µΓdS S ΓdS , (26)
où S définit le domaine délimité verticalement par les deux parois et horizontalement par une période. Cette modélisation se révèle pertinente pour des fluides faiblement à modérément rhéofluidifiant (n p > 0.5).

Balmforth et Rust [37] ont effectué une analyse faiblement non linéaire de la convec- Bouteraa et al [38] ont étendu les calculs de [37] au cas général de glissement avec frottement. En situation de non-glissement, les auteurs trouvent α c ≈ 2.15 10 -4 . Ils montrent aussi que la variation non linéaire de la viscosité avec le cisaillement favorise une convection sous forme de rouleaux contrairement à ce qui a été indiqué par [39].

Les simulations numériques de convection de Rayleigh-Bénard pour un fluide en loi puissance dans une cavité horizontale allongée ont été reprises par [40]. The study of rolls stability shows that : (i) for low shear-thinning effects, the band of stable wavenumbers is bounded by zigzag instability and cross-roll instability. Furthermore, the marginal cross-roll stability boundary enlarges with increasing shear-thinning properties ;

(ii) for high shear-thinning effects, Eckhaus instability becomes more dangerous than crossroll instability. For square patterns, the wavenumber selection is always restricted by zigzag instability and by 'rectangular Eckhaus' instability. In addition, the width of the stable wavenumber decreases with increasing shear-thinning effects. Numerical simulations of the planform evolution are also presented to illustrate the different instabilities considered in the paper.

Introduction

Studies on patterns formation and their stability in Rayleigh-Bénard convection for Newtonian fluids have been considered in several papers. A review can be found in books of Getling [1] and Koschmieder [2] and more recently in Bodenschatz et al. [3] where the most significant progress in the field is identified. Comparatively to the Newtonian case, only a limited number of studies were devoted to non-Newtonian fluids and still fewer to nonlinear developments. Yet, non-Newtonian fluids intervene in a very broad range of industrial processes such as polymer and foodstuffs processing and in complex physical phenomena such as the convective movements in the earth's mantle. Here, we focus on the shear-thinning behavior, i.e. non linear decrease of the effective viscosity with the shear-rate, which is the most common property of non-Newtonian fluids. In recent articles ( [4], [5] and [6]), the nature and the stability of patterns which emerge in Rayleigh-Bénard convection for shear-thinning fluids have been studied using a weakly non linear analysis.

Boussinesq approximations have been adopted and the slabs have been considered as

1.2. INTRODUCTION 29
perfectly conducting. Using Carreau model to describe the shear-thinning behavior of the fluid, it has been shown in [6] that : (i) rolls are the only stable convective patterns at threshold and (ii) there is a critical value of the shear-thinning degree α defined by Eq.(1.12) above which the bifurcation becomes subcritical.

Most analyzes consider ideal situations where the bounding horizontal surfaces are perfect conductors of heat. However in many laboratory experiments and in engineering and geophysical problems, the slabs have a finite conductivity and they are not better conductor than the fluid itself. In this case, the temperature disturbances do not vanish on the boundaries. The thermal boundary conditions that have to be satisfied are the continuity of temperature and heat flux. According to Cerisier et al. [7], the temperature fluctuation occurring in the liquid close to a nearly insulating slab distorts the temperature distribution. This temperature distortion can lead to an instability of the fluid layer. As a consequence, at threshold, the temperature gradient is small and the fluid organizes in a pattern with a small wavenumber. Furthermore, theoretical and experimental studies show that squares may be the convection patterns at the onset instead of rolls. Experimental evidence of square patterns was reported by Legal Pocheau and Croquette [8] and Legal and Croquette [9]. The competition between roll and square patterns for a Newtonian fluid has been examined in weakly supercritical Rayleigh-Bénard convection by Busse and Riahi [10], Proctor [11] and Jenkins and Proctor [12]. The results are presented in terms of the Prandtl number Pr, and the ratio χ of the thermal conductivities of slabs and fluid. It has been shown for instance that for Pr ≥ 10 and for slabs with thickness equal to the width of the fluid layer, that the convective pattern at threshold is in form of squares when χ < χ c = 1. These studies were extended recently to shear-thinning fluids by Bouteraa and Nouar [13]. It has been found that the critical value χ c , below which squares are stable, decreases with increasing shear-thinning effects.

Recently, experiments were done by Kebiche [14], using carboxymethylcellulose (CMC) solutions as shear-thinning fluid. In the Rayleigh-Bénard setup, the slabs are made of polycarbonate with a ratio of thermal conductivities χ ≈ 0.25. PIV measurements done in one vertical section do not allow the determination of convection pattern type.

These studies are valid only in the immediate vicinity of the threshold with perfectly periodic pattern. However, as the Rayleigh number Ra is increased above the onset Ra c , the growth-rate of the perturbation is positive for any wavenumber within a band √ ǫ with ǫ = (Ra -Ra c )/Ra c , around the critical wavenumber k c . A wavepacket centered on k c can be constructed with the unstable modes. The corresponding convective pattern is modulated spatially on a scale of O (1/ √ ǫ). The stability analysis of these convective patterns with respect to long wavelength perturbations is of great interest. It allows in particular, the determination of the range of stable wavenumbers. Typical instability mechanisms are Eckhaus (E), zigzag (ZZ) and cross-roll (CR) instabilities [15]. Eckhaus instability is a phase instability which acts on the roll phase to change the wavelength, CHAPITRE 1.

compressing or dilating the pattern. Zigzag is also a phase instability which arises from perturbations parallel to the roll axes : it creates undulations along the roll axes when the wavelength is too large. Cross-roll instability is an amplitude instability which consists of a set of rolls growing perpendicularly to the original pattern. For a Newtonian fluid with perfectly heat conductive slabs, the instability mechanisms which tend to limit the stability region of rolls depend on the Prandtl number and on the boundary conditions, rigid or stress-free boundary conditions, as shown by Busse [15] and Bolton and Busse [16]. Square patterns are also subject to long wavelength instabilities. In the case of poorly conducting slabs and for Newtonian fluids, Hoyle [18] has shown that the range of stable wavenumbers is restricted by zigzag instability and by Eckhaus rectangular instability.

According to Hoyle [18], this latter instability has a three-dimensional character since the system responds differently in each of the two horizontal directions. It behaves like, one of two rolls, that constitute square pattern, grows locally at the expense of the other. This is why, Holmedal et al. [19] called this instability "Long wavelength cross-roll instability".

This study was extended by Holmedal [19] to the general case of slabs with different finite conductivities and thicknesses of the slabs.

The objective of this paper is to investigate the influence of shear-thinning effects on the stability of the convective patterns and the width of the stable band of wavenumbers in Rayleigh-Bénard convection with slabs of finite conductivity and arbitrary thickness. The rheological law introduces additional nonlinearity and coupling between flow variables. This additional nonlinearity will induce stronger interactions between the two sets of rolls that constitute square patterns, than in the Newtonian case. Therefore, shear-thinning effects will modify not only the range of stable wavenumber but also the more restrictive instability mechanism. A weakly nonlinear analysis based on the amplitude equations formalism ( [18], [20]) is adopted as a first approach to examine nonlinear effects arising from the rheology.

To our knowledge, the present study is the first one which considers the influence of the rheology on secondary instabilities. The structure of the paper is as follows. In §2 the problem is formulated. Linear stability analysis is briefly considered in §3. The nature of the primary bifurcation and pattern selection are discussed in §4. It is observed that

1.3. PROBLEM FORMULATION
shear-thinning effects favor formation of rolls. The stability of the convective patterns with respect to inhomogeneous spatial perturbations is investigated in §5, using the amplitude equations formalism. Influence of shear-thinning effects is highlighted. In §6, time evolution of the convective pattern is illustrated from the numerical simulation of amplitude equations. The paper ends with a conclusion where the relevant results are summarized.

Problem formulation

We consider a shear-thinning fluid layer of infinite horizontal extent which is heated from below and cooled from above. We assume that the rigid slabs which enclose the fluid have arbitrary conductivities and thicknesses. The thermal conductivity and diffusivity are noted K and κ for the fluid and Kp and κp for the slabs. We define χ as the ratio of Kp and K and we assume as in [11] and [21] that χ = Kp K = κp κ . This assumption is reasonable for several couples (fluid, slab) where the ratio of the thermal capacities of the fluid and the slabs r = (ρ Ĉp ) f luid (ρ Ĉp ) slabs remains of order 1 : for instance, r(water, copper) = 1.22 and r(glycerin, glass) = 1.67 .

Dimensional quantities are denoted with the symbol hat (ˆ). In the following, we note d the depth of the fluid layer, ∆ T = T1 -T2 , the temperature difference between the outer surfaces of the upper and lower slabs. Because of the thermal expansion, the temperature difference between the two plates induces a vertical density stratification. Heavy cold fluid is above a warm light fluid. For small ∆ T , the fluid remains motionless and the heat is transferred by conduction with a linear temperature profile across the fluid layer. In the fluid, -d/2 < ẑ < d/2, the hydrostatic solution and the temperature profile are :

d P dẑ = -ρ 0 ĝ 1 -β T -T0 and Tcond = T0 - ∆ T 1 + 2Λ/χ ẑd , (1.1)
where ĝ is the acceleration due to gravity and Λ the dimensionless thickness of slabs.

The z-axis is directed upwards with the origin located at the middle of the fluid layer.

The reference temperature T0 is the temperature at the middle of the fluid layer, ρ0

the fluid density at T0 and μ0 is the zero-shear rate viscosity at T0 . The temperature difference between the top and the bottom of the fluid layer is ∆ Tf = ∆ T [1 + 2Λ/χ]. The temperature profiles in top and bottom plates are

Tcond = T0 + ∆ T 2Λ + χ 1 2 - 1 2 χ - ẑd ; d 2 ≤ ẑ ≤ 1 2 + Λ d (1.2)
and

Tcond = T0 + ∆ T 2Λ + χ 1 2 χ - 1 2 - ẑd ; -Λ + 1 2 d ≤ ẑ ≤ - d 2 . (1.3) CHAPITRE 1.
When the top and bottom plates are poor thermal conductors, a large part of ∆ T occurs across the plates, and remains only a small part ∆ Tf of ∆ T , acting as a driving force for the convection. When ∆ Tf exceeds a critical value, the convection sets in and convective patterns emerge. The stability of the hydrostatic solution is considered by introducing temperature and pressure perturbations as well as a fluid motion. The fluid is incompressible and Boussinesq approximations are adopted. We use d, μ0 , d2 κ , κ d , ρ0 κ2 d2 and ∆ T Ra as characteristic scales of length, viscosity, time, velocity, pressure and temperature. Using these scales, the perturbation equations read

∇.v = 0, (1.4)
1 Pr

∂v ∂t + (v.∇)v = -∇p + θe z + ∇.τ , (1.5
)

∂θ ∂t + (v.∇) θ = Ra w + ∆θ, (1.6
)

∂θ p ∂t = χ∆θ p , (1.7)
where v = ue x + ve y + we z , p are respectively the velocity and pressure perturbations, θ and θ p are temperature perturbations in the fluid and the slabs respectively and τ the deviatoric of the stress tensor. The Prandtl Pr and Rayleigh Ra numbers are defined by : Pr = μ0 ρ0 κ Ra = ρ0 ĝ β∆ T d3 κμ 0 .

We consider a purely viscous shear-thinning fluid,

τ = µ(Γ) γ, (1.8)
where Γ is the second invariant of the strain-rate tensor :

Γ = 1 2 γij γij ; γij = ∂v i ∂x j + ∂v j ∂x i (1.9)
where v i are the components of the velocity and x i are the spatial coordinates.

Concerning the nonlinear rheological law µ(Γ), we have used a Carreau model [22].

In dimensional form, it is given by :

μ -μ∞ μ0 -μ∞ = (1 + λ2 Γ) nc-1 2 , (1.10)
where μ∞ is the infinite-shear viscosity, μ0 the zero-shear viscosity, λ a characteristic time of the fluid, n c the shear-thinning index. For several polymer solutions, μ∞ << μ0 [23].

Hence neglecting μ∞ with respect to μ0 , we have in dimensionless form :

µ = (1 + λ 2 Γ) nc-1 2 , (1.11)

PROBLEM FORMULATION

33

A Taylor series expansion of µ around the base state (where the fluid is at rest) allows to define the degree of shear-thinning of the fluid as :

α = dµ dΓ Γ=0 = 1 -n c 2 λ 2 .
(1.12)

No-slip and non penetration boundary conditions as well as continuity of temperature and heat flux at the interface slabs-fluid read :

v z = ± 1 2 = 0, (1.13
)

θ z = ± 1 2 = θ p z = ± 1 2 , (1.14)
∂θ ∂z z = ± 1 2 = χ ∂θ p ∂z z = ± 1 2 . (1.15)
Temperatures of the outer surfaces of the upper and lower slabs are fixed, thus :

θ p z = 1 2 + Λ = θ p z = - 1 2 -Λ = 0 (1.16)
In the momentum equations, the pressure field can be eliminated using the curl of Eq.

(2.3). We then take the curl of Eq. (2.3) one more time. Using the continuity equation, and projecting onto e z , we get the following evolution equations for the vertical vorticity ζ and the vertical velocity w :

1 Pr

∂ζ ∂t + e z • ∇ × [(v • ∇) v] = ∆ζ + e z • ∇ × [∇ • (µ -1) γ] , (1.17)
1 Pr

∂∇ 2 w ∂t -e z • [∇ × ∇ × [(v.∇) v]] = ∆ 2 w + ∇ 2 H θ - (1.18) [∇ × ∇ × [∇ • (µ -1) γ]] • e z , ∂θ ∂t + (v • ∇) θ = Ra w + ∇ 2 θ , (1.19
)

∂θ p ∂t = χ∇ 2 θ p , (1.20)
where

ζ = ∂v ∂x - ∂u ∂y and ∇ 2 H = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 .
From the continuity equation and the vertical vorticity definition, one can deduce the horizontal velocity components (u , v) :

∇ 2 H u = - ∂ 2 w ∂x∂z - ∂ζ ∂y ; ∇ 2 H v = - ∂ 2 w ∂y∂z + ∂ζ ∂x . (1.21) CHAPITRE 1.
The boundary conditions for w are :

w = Dw = 0 at z = ±1/2, (1.22)
For the temperature, the boundary conditions are :

θ p = 0 at z = ± (1/2 + Λ) , (1.23)
θ = θ p at z = ±1/2, (1.24
)

Dθ = χDθ p at z = ±1/2. (1.25)
Five dimensionless parameters appear in the governing equations : the Rayleigh number Ra, the Prandtl number Pr, the thermal conductivities ratio χ, the dimensionless thickness of the slab Λ, and the shear-thinning degree α. In the present study, Pr = 10. Actually, our results do not vary significantly with Pr when Pr ≥ 10.

1.4

Linear stability analysis

   =    F 11 (z) G 11 (z) G p11 (z)    exp (i k • r + s t) , (1.30)
1.4. LINEAR STABILITY ANALYSIS

35

where r = (x, y) is the vector position in the horizontal plane, and k is the wave-vector.

Substituting (1.30) into (1.27)-(1.29) leads to the differential equations

s Pr -1 D 2 -k 2 F 11 = -k 2 G 11 + D 2 -k 2 2 F 11 , (1.31)
s G 11 = Ra F 11 + (D 2 -k 2 )G 11 , (1.32)
s G p11 = χ(D 2 -k 2 ) G p11 , (1.33)
with k = |k|. The boundary conditions are :

F 11 = DF 11 = 0 at z = 0, 1 (1.34)
G p11 = 0 at z = -Λ, 1 + Λ, (1.35)
G 11 = G p11 at z = 0, 1 (1.36)
DG 11 = χDG p11 at z = 0, 1 .

(1.37)

The eigenvalue problem (1.31)-(1.33) with the boundary conditions (1.34)-(1.37) is solved using a Chebyshev collocation method. The marginal stability curve Ra(k) is determined by the condition s = 0. The minimum of the marginal stability curve gives the critical Rayleigh number Ra c and critical wavenumber k c . We recover the results of [13] for Λ = 1

and we extend them to other thicknesses Λ on Fig. 1.1. We observe that k c and Ra c decrease with decreasing the ratio χ of conductivities. Actually, k c and Ra c vary from respectively, 3.11 and 1708 to 0 and 720 ( [24], [11]). An explanation of this evolution can be found in [25] in the limit of perfectly insulating slabs. From a physical point of view, such configuration means that the temperature field is fixed in the solid (or evolves on a very long time scale compared to that of the fluid). As a consequence, the temperature gradient, and therefore the energy flux, is fixed in the solid. Hence, the temperature fluctuations at the interface do not propagate inside the solid and primary bifurcation needs less energy to occur which explains the decrease of Ra c with decreasing χ.

Remark

Linear stability analysis gives the critical Rayleigh number Ra c for instability onset and determines the modulus k c of the critical wave-vector k of the unstable modes. The direction of k is arbitrary. This orientation degeneracy is related to the isotropy of the horizontal plane [26]. There is also a pattern degeneracy that results from the linear theory itself ; indeed, any superposition of normal modes

[w(r, z), θ(r, z), θ p (r, z)] = ℓ c ℓ exp (ik ℓ • r) [F 11 (z), G 11 (z), G 11p (z)] (1.38)
with |k ℓ | = k c and where the c ℓ 's are constant coefficients, is also a solution of the linear problem with a zero growth rate at criticality.

(2)

(1) (4)

Characteristic time of the instability

Near the onset of convection, the growth rate s of the perturbation may be approximated using Taylor expansion

s = ǫ τ 0 + O(ǫ 2 ) with ǫ = Ra -Ra c Ra c , (1.39)
where τ 0 is the characteristic time for the instability to grow. It is given by

τ 0 = Ra c ∂s ∂Ra kc,Rac
. It is obtained from the curve, temporal amplification rate s versus Rayleigh number Ra, at the critical conditions. Its evolution is represented in Fig. 1.2. The increase of τ 0 with decreasing χ is related to the increase of the thermal diffusion time in the slab as kp diminishes. For χ > 10, we recover the value corresponding to a perfect heat conductor, τ 0 = 0.053. Note that τ 0 does not depend on the rheological parameters.

Pattern selection at the onset of convection

The selection of the convective pattern is determined by the non linearities of the problem, i.e. nonlinear inertial and nonlinear viscous terms. A weakly nonlinear analysis based on amplitude expansion method similar to that considered in [27], [28], [29], [30], [6] is used as a first approach to investigate nonlinear effects on the competition between convective patterns. Actually, the patterns that emerge near the onset of convection are either rolls or squares. Further calculations show that hexagons are unstable ( [6], [12]).

CHAPITRE 1.

Principles of the amplitude expansion method : case of square

For a square pattern, the fundamental solution in the linear regime is A e ik 1 •r + B e ik 2 •r ψ 11 , with k 2 orthogonal to k 1 (the two vectors k 1 and k 2 have the same modulus), A and B are the complex amplitudes of the perturbation along the two wavevectors and ψ 11 (z) stands for F 11 (z), G 11 (z) or G p11 (z). The interaction of the fundamental solution with itself, through the quadratic nonlinear inertial terms produces the first harmonic A 2 e 2ik 1 •r + B 2 e 2ik 2 •r ψ 22 and a coupling between modes k 1 and k 2 ,

ABe i(k 1 +k 2 )•r ψ AB .
The interaction of the fundamental with its complex conjugate leads to an other coupling between modes k 1 and k 2 , AB * e i(k 1 -k 2 )•r ψ AB * , where (.) * denotes the complex conjugate and a correction of the base state,

| A 2 | + | B 2 | ψ 02
. The feedback at the cubic order on the fundamental solution through nonlinear inertial and viscous terms

is | A 2 | + | B 2 | A e ik 1 •r + B e ik 2 •r ψ 13 .
From this cascade of nonlinear interactions, the nonlinear solution can be written as

ψ (r, z, t) = A(t) e ik 1 •r + B(t) e ik 2 •r ψ 11 (z) + c.c. + (1.40) A 2 (t) e 2ik 1 •r + B 2 (t) e 2ik 2 •r ψ 22 (z) + A(t)B(t)e i(k 1 +k 2 )•r ψ AB (z) + c.c. + | A 2 (t) | + | B 2 (t) | ψ 02 (z) + A(t)B * (t)e i(k 1 -k 2 )•r ψ AB * (z) + c.c. + | A 2 (t) | + | B 2 (t) | A(t) e ik
dA dt = ǫ τ 0 A -(g 1 |A| 2 + β|B| 2 )A, (1.41)
dB dt = ǫ τ 0 B -(g 1 |B| 2 + β|A| 2 )B, (1.42)
where g 1 and β are respectively self-saturation and crossed-saturation coefficients. The form of the amplitude equations (1.41) and (1.42) is completely determined by the rules of invariance via symmetry by rotation of an angle π/2 and by translation [31], [32].

Substituting (1.40)-(1.42) into (1.18)-(1.20) yields after some algebra to a set of differential equations for each mode that are solved sequentially. To avoid secular terms at the cubic order, compatibility conditions have to be enforced using the Fredholm alternative. The latter states that the resonating forcing terms have to be orthogonal to the kernel of the adjoint of the linear operator. This allows the determination of landau saturation coefficients g 1 and β.

1.5. PATTERN SELECTION AT THE ONSET OF CONVECTION

Nature of the primary bifurcation

As shown in [6], the self-saturation g 1 and crossed saturation β coefficients can be written as the sum of Newtonian (N superscript) and non-Newtonian contributions (nN superscript) :

g 1 = g N 1 -αg nN 1 , (1.43)
β = β N -αβ nN . (1.44)
It is therefore possible to define a critical value α c of the shear-thinning degree above which the bifurcation becomes subcritical.

In the case of rolls, β = 0 and

α c = g N 1 g nN 1 .
(1.45)

In the case of squares, β = 0 and

α c = g N 1 + β N g nN 1 + β nN .
(1.46)

Variations of α c with χ for different Λ are depicted in Fig. 1.3 (a) and (b) for rolls and squares respectively. For large χ, the asymptotic limit of α is α c = 2.15×10 -4 in agreement with [6] and [4]. With decreasing Kp (decreasing χ), the intensity of convection decreases, therefore, it is not surprising that stronger shear-thinning effects are needed to obtain a subcritical bifurcation.

Convective patterns at threshold

A linear stability analysis of stationary roll and square solutions of Eqs. (41) and (42) allows to show that squares are stable when β < g 1 , i.e. when the coupling between the two orthogonal modes that describe the square pattern is weak enough. By contrast, when β > g 1 , the coupling is too strong, the squares lose their stability and rolls are the stable nonlinear state. A similar phenomenological description can be found in [26]. Figure 1.4 shows that with increasing α, β g 1 increases and thus the interaction between the two orthogonal modes k 1 and k 2 becomes stronger. A possible interpretation may be related to the reduction of viscosity with increasing shear-thinning effects, which leads to an increase of the convection intensity. Nonlinearities and coupling between modes become stronger which favor roll patterns.

Using shear-thinning decomposition of g 1 and β (Eqs. (1.43) and (1.44)), it is found that rolls are stable when α > α S-R , with :

α S-R = β N -g N 1 β nN -g nN 1 .
(1.47) For the limit of Newtonian fluids i.e. α = 0, we recover the results of the literature [12] : for Pr ≥ 10 and Λ = 1, rolls are stable patterns provided that χ > χ c = 1 . We notice that the domain of stability of squares shrinks when Λ decreases which is understandable.

Indeed, the thinner the slabs are, the weaker the thermal resistance is. Then the problem is closer to the case of perfectly conducting slabs where rolls are the preferred patterns.

For Λ > 1, the dependence of α S-R with respect to Λ is weak.

Finally, we observe that α S-R increases as χ decreases. Poorly conducting slabs favor square patterns as shown in [10], [11] and [12], so stronger shear-thinning effects are necessary so that rolls become the preferred planform. This last result is in agreement with [13].

Flow structure, viscosity and temperature fields

In this section, features of the flow, temperature distribution and shear-thinning effects on the viscosity field in a roll and a square solutions are studied for highly and poorly conducting slabs.

Case of highly conducting walls : χ = 100, Λ = 1

The flow structure and the viscosity field for a roll solution are illustrated by Fig. 1.6.

The interior of the roll is practically isoviscous with µ ≈ 1. The viscosity is minimal at the walls where the shear-rate γxz is maximal. It is also weakly reduced at the four corners square cell are shown in Fig. 1.9. It is worthy to note that the temperature perturbation does not vanish at the walls and the vertical thermal gradient is weak.

Comparison between roll and square solutions for fixed χ and λ

For fixed values χ and λ, velocity and viscosity fields are determined for roll and square solutions. It is observed that the maximum of shear rate and therefore the minimum of viscosity occurs for the stable pattern.

Heat transfer

The heat transfer through the horizontal fluid layer is described by the Nusselt number, Nu, the ratio of the total heat to the purely conductive heat flux, i.e. when the fluid is at rest.

Nu = 1 - ∂ θ ∂z z=-1/2 = 1 -| A | 2 + | B | 2 (DG 02 ) z=-1/2 . (1.48)
where |A| 2 + |B| 2 G 02 is the modification at the second order of the conductive temperature profile due to the interaction of the fundamental mode with its complex conjugate.

The overbar denotes the average over one wavelength. Using the stationary solutions of the amplitude equations, one obtains :

Nu r = 1 - ǫ τ 0 g 1 (DG 02 ) z=-1/2 for rolls, (1.49) Nu s = 1 - 2ǫ τ 0 (g 1 + β) (DG 02 ) z=-1/2 for squares. (1.50)
The Figure 1.10 shows the variation of the Nusselt number as a function of the shearthinning degree at ǫ = 0.1, χ = 0.01 and four different values of Λ. The Nusselt number increases with increasing shear-thinning effects in agreement with [33], [34], [6]. This is a consequence of the increase of the rolls amplitude. As expected, the Nusselt number decreases significantly with increasing the slab thickness. The difference between Nu rolls and Nu squares is small. Nevertheless, Nu is larger for the stable convective pattern in agreement with the maximum heat transfer principle : "the only stable solution is the one of maximum heat transport [35], [36]". This is also illustrated by Fig. 1.11 where Nu is represented as a function of ǫ for given α = 4 × 10 -4 and χ = 10 -2 . At Λ = 1 and Λ = 0.1 squares are stable and rolls are unstable, whereas at Λ = 0.01, rolls are stable. Other principles can be considered to predict the stable pattern such the maximum entropy production or the maximum viscous dissipation [37], [38]. Indeed, it can be shown that for a steady solution, Ra (Nu -1) = Ω τ ij γij dΩ [39], where Ω is a domain delimited by the top and bottom walls and one wavelength in the x-and y-directions. For shear-thinning fluids one can consider the principle of maximum viscosity reduction, as indicated in the previous section.

0 2 4 6 8 10 -4 10 -3 10 -2 10 -1 (1) (2) (3) (4)

Secondary instabilities

Departing from the critical conditions, a band of wavenumbers of width O( √ ǫ), centered on k = k c will now have a positive growth-rate. The wavepacket centered on the most unstable wavenumber can be considered as a monochromatic wave, with complex amplitude modulated in space and time.

In a square lattice, and up to third order in the perturbations, the spatio-temporal evolution of the amplitudes is described by a set of two coupled Ginzburg-Landau equations derived by Newell & Whitehead [20] and Segel [40] :

∂A ∂t = ǫ τ 0 A + ξ 2 0 τ 0 ∂ ∂x - i 2k c ∂ 2 ∂y 2 2 A -(g 1 |A| 2 + β|B| 2 )A, (1.51)
∂B ∂t = ǫ τ 0 B + ξ 2 0 τ 0 ∂ ∂y - i 2k c ∂ 2 ∂x 2 2 B -(g 1 |B| 2 + β|A| 2 )B, (1.52)
where the coherence length ξ 0 is defined by

ξ 2 0 = 1 2Ra c ∂ 2 Ra ∂k 2 Rac,kc
. It does not depend on rheological properties and can be calculated from the curve of the growth-rate σ versus (kk c ). For χ ≤ 10, the coherence length ξ 0 varies between 0.375 and 0.415 . For χ > 10, we obtain ξ 0 = 0.386 which is in agreement with the literature [41].

Stability of a roll modulated solution

In the following, we consider first the case where rolls emerge at primary bifurcation (α > α S-R ). We look for a stationary solution of the system (1.51)-(1.52), of the form :

A 0 = R 0 exp(iqx) where q = k -k c , (1.53)
B 0 = 0 . (1.54)
Substituting the above expressions into (1.51) leads to :

R 0 = ǫ -ξ 2 0 q 2 g 1 τ 0 . (1.55)
Thereafter, we examine the stability of the stationary solution (1.53), (1.54) with respect to infinitesimal perturbations, in terms of amplitudes r A and r B and phases Φ A and Φ B .

1.6. SECONDARY INSTABILITIES 47
The perturbed solution can be written as :

A(x, y, t) = (R 0 + r A (x, y, t)) exp( i (Φ A (x, y, t) + qx)) (1.56) B(x, y, t) = r B (x, y, t) exp( i Φ B (x, y, t)) (1.57)
Substituting expressions (1.56) and (1.57) into (1.51) and (1.52), we obtain after linearization and separating the real and imaginary parts of the equations :

∂r A ∂t = -2g 1 R 2 0 r A + ξ 2 0 τ 0 ∂ 2 r A ∂x 2 -2qR 0 ∂Φ A ∂x + q k c ∂ 2 r A ∂y 2 + R 0 k c ∂ 3 Φ A ∂x∂y 2 - 1 4k 2 c ∂ 4 r A ∂y 4 , (1.58) ∂Φ A ∂t = ξ 2 0 τ 0 ∂ 2 Φ A ∂x 2 + 2q R 0 ∂r A ∂x + q k c ∂ 2 Φ A ∂y 2 - 1 k c R 0 ∂ 3 r A ∂x∂y 2 - 1 4k 2 c ∂ 4 Φ A ∂y 4 ,
(1.59)

∂r B ∂t = (g 1 -β) R 2 0 r B + q 2 ξ 2 0 τ 0 r B + ξ 2 0 τ 0 ∂ 2 r B ∂y 2 - 1 4kc 2 ∂ 4 r B ∂x 4 .
(1.60)

Equation associated with ∂Φ B /∂t does not contain linear terms so Φ B does not intervene at the first order. Using normal mode decomposition i.e. Ψ(x, y, t) = Ψ exp(σt

+ i(Q 1 x + Q 2 y))
where Ψ stands for r A , Φ A and r B , an eigenvalue problem is derived :

LX = σX (1.61) where X = r A , Φ A , r B T
is the eigenvector, σ the eigenvalue and L a 3 × 3 square matrix arising from equations (1.58), (1.59) and (1.60). Note that the eigenvalue problem corresponding to (1.60) can be solved independently from the whole system.

We consider the long wavelength limit where Q 1 → 0 and Q 2 → 0. In that case, the relevant eigenvalues of the former system are given by :

σ 1 = -2g 1 R 2 0 + O (Q 1 ) , (1.62)
σ 2 = - ξ 2 0 τ 0 Q 2 1 1 - 2ξ 2 0 τ 0 q 2 g 1 R 2 0 - ξ 2 0 τ 0 qQ 2 2 k c + O Q 2 1 Q 2 2 , (1.63)
σ 3 = ξ 2 0 τ 0 q 2 + (g 1 -β)R 2 0 + O(Q 4 1 , Q 2 2 ). (1.64)
The eigenvector ( r A , Φ A ) associated with the first eigenvalue

σ 1 = -2g 1 R 2 0 < 0 is (O(1/Q 1 ), 1)
). Therefore, σ 1 describes the relatively rapid relaxation of the amplitude perturbation r A to its equilibrium value. The eigenvector ( r A , Φ A ) associated with σ 2 is (O(Q 1 ), 1). The second eigenvalue describes the evolution of the phase perturbation Φ A .

The third root σ 3 describes the evolution of rolls growing perpendicularly to the original ones.

The eigenvalue σ 2 can also be derived using the phase approximation. This ap-CHAPITRE 1.

proach described in [32] and [42] relies on the fact that the amplitude r A relaxes quickly with time, it can be considered to be adiabatically slaved to the phase Φ A . This comes down to writing ∂r A /∂t = 0. Furthermore, in the long wavelength limit, spatial derivatives are very small compared to the variables themselves. Therefore, the amplitude r A is approximately given by its adiabatic value

r A = - qξ 2 0 g 1 R 0 τ 0 ∂Φ A ∂x . (1.65)
This expression is substituted in (1.59) to determine the evolution of φ A . A phase-diffusion equation is then derived :

∂Φ A ∂t = D ∂ 2 Φ A ∂x 2 + D ⊥ ∂ 2 Φ A ∂y 2 .
(1.66)

The longitudinal D and transverse D ⊥ phase diffusion coefficients are given by

D = ξ 2 0 τ 0 1 - 2q 2 g 1 R 2 0 ξ 2 0 τ 0 and D ⊥ = ξ 2 0 τ 0 q k c . (1.67)
Equation (1.66) shows that a perturbation of the wavenumber leads to a readjustment of the system through a phase diffusion process. The eigenvalue which stems from (1.66) is the same as σ 2 .

Eckhaus instability

For perturbations that vary only in the x-direction (Q 2 = 0), the eigenvalue σ 2 (1.63) reduces to

σ 2 = -Q 2 1 D + O Q 4 1 .
(1.68)

For positive longitudinal phase-diffusion coefficient D , σ 2 is negative, the perturbation is

damped and the roll solution (1.53) is stable. Using (1.55), the stability is satisfied if

ǫ > ǫ E = 3q 2 ξ 2 0 , (1.69)
where the subscript E means "Eckhaus". Note that, this instability does not depend on the rheological parameters.

Zigzag instability

For perturbations that vary only in the y-direction (Q 1 = 0), the eigenvalue σ 2 (1.63) reduces to

σ 2 = -Q 2 2 D ⊥ + O Q 4 2 (1.70) 1.6. SECONDARY INSTABILITIES 49
For negative D ⊥ , i.e. when q is negative (k < k c ), rolls at wavelength greater than the critical one, the eigenvalue σ 2 is positive, the perturbation is amplified and the roll solution (1.53) is unstable. In this case, the rolls will saturate into bends that decrease the wavelength.

Cross-Roll instability

The eigenvalue σ 3 (1.64) corresponds to the cross-roll (CR) instability. The system is CR stable if σ 3 < 0. Using (1.55), the system is CR stable if

ǫ > ǫ CR = β β -g 1 q 2 ξ 2 0 , (1.71)
where the subscript CR means " cross-roll ". When ǫ < ǫ CR , the stationary roll solution (1.53) becomes unstable : new rolls expand perpendicularly. It can be shown straightforwardly that the cross-roll is a more restrictive instability than the Eckhaus instability, when

β < 3 2 g 1 . (1.72)
-For a Newtonian fluid at Pr = 10 and in the case of perfectly conducting walls, β/g 1 = 1.242. Decreasing the conductivity of the wall, χ, will decrease (β/g 1 ) and will give a narrow band of stable rolls. The width of this band vanishes as the singularity

β/g 1 → 1 + is approached.
-For a shear-thinning fluid, β/g 1 increases with increasing shear-thinning effects. Combining (1.72) with (1.43) and (1.44), we can define a shear-thinning degree α CR-E , below

which cross-roll is the more restrictive instability :

α CR-E = 3g N 1 -2β N 3g nN 1 -2β nN .
(1.73)

For α > α CR-E , the Eckhaus instability takes over as the most restrictive of the two. The

(

Stability of a square modulated solution

In the case where the convection starts with perfect square patterns, a stationary solution is given by :

A 0 (x) = R 0 exp(iqx) where q = k -k c , (1.74) B 0 (y) = R 0 exp(iqy) . (1.75)
Replacing these expressions in (1.51) and (1.52) leads to :

R 0 = ǫ -ξ 2 0 q 2 (g 1 + β)τ 0 . (1.76)
As in the case of rolls, we carry out a linear stability analysis of the stationary square solution (1.74)-(1.76). A perturbation of the form

A(x, y, t) = (R 0 + r A (x, y, t)) exp( i (Φ A (x, y, t) + qx)) , (1.77) B(x, y, t) = (R 0 + r B (x, y, t)) exp( i (Φ B (x, y, t) + qy)) , (1.78) 52 CHAPITRE 1.
is introduced. Substituting, the amplitudes A and B by their expressions (1.77), (1.78) into (1.51) and (1.52) leads after linearization :

∂rA ∂t = -2g1R 2 0 rA -2βR 2 0 rB + ξ 2 0 τ0 ∂ 2 rA ∂x 2 -2qR0 ∂ΦA ∂x + q kc ∂ 2 rA ∂y 2 + R0 kc ∂ 3 ΦA ∂x∂y 2 - 1 4k 2 c ∂ 4 rA ∂y 4 , (1.79)
∂ΦA ∂t = ξ 2 0 τ0 ∂ 2 ΦA ∂x 2 + 2q R0 ∂rA ∂x - 1 kcR0 ∂ 3 rA ∂x∂y 2 + q kc ∂ 2 ΦA ∂y 2 - 1 4k 2 c ∂ 4 ΦA ∂y 4 , (1.80) ∂rB ∂t = -2g1R 2 0 rB -2βR 2 0 rA + ξ 2 0 τ0 ∂ 2 rB ∂y 2 -2qR0 ∂ΦB ∂y + q kc ∂ 2 rB ∂x 2 + R0 kc ∂ 3 ΦB ∂x 2 ∂y - 1 4k 2 c ∂ 4 rB ∂x 4 , (1.81)
∂ΦB ∂t = ξ 2 0 τ0 ∂ 2 ΦB ∂y 2 + 2q R0 ∂rB ∂y - 1 kcR0 ∂ 3 rB ∂x 2 ∂y + q kc ∂ 2 ΦB ∂x 2 - 1 4k 2 c ∂ 4 ΦB ∂x 4 . (1.82)
Using a normal modes decomposition, i.e. Ψ(x, y, t) = Ψ exp (σt + i (Q

1 x + Q 2 y))
, where Ψ stands for r A , Φ A , r B , Φ B , the following eigenvalue problem is derived :

MX = σX . (1.83) In Equation (1.83), X = r A , Φ A , r B , Φ B T
is the eigenvector, σ the eigenvalue and M the 4 × 4 square matrix arising from equations (1.79)-(1.82). The eigenvalues and eigenvectors can be determined numerically using Matlab. Examples of results are shown in Appendix 1.A.

Actually, we are particularly interested by the long wavelength limit approach, i.e. Q 1 → 0 and Q 2 → 0. In this approach, the eigenvalues σ 1 and σ 2 associated with the amplitudes r A and r B respectively, are given by

σ 1 = -2R 2 0 (g 1 + β) + O Q 2 1 , Q 2 2 ; σ 2 = -2R 2 0 (g 1 -β) + O Q 2 1 , Q 2 2 (1.84)
Since, in this case β and g 1 are positive and β < g 1 , the amplitude modes r A and r B decrease quickly with time and can be considered adiabatically slaved to the phase modes Φ A and Φ B . They can be approximated by their adiabatic values :

r A = q R 0 (g 2 1 -β 2 ) ξ 2 0 τ 0 β ∂Φ B ∂y -g 1 ∂Φ A ∂x , (1.85)
r B = q R 0 (g 2 1 -β 2 ) ξ 2 0 τ 0 β ∂Φ A ∂x -g 1 ∂Φ B ∂y . (1.86)
Substituting these expressions in (1.80) and in (1.82) leads to the following diffusion equations of phases Φ A and Φ B :

∂Φ A ∂t = D ∂ 2 Φ A ∂x 2 + D ⊥ ∂ 2 Φ A ∂y 2 + D xy ∂ 2 Φ B ∂x∂y , (1.87) ∂Φ B ∂t = D ∂ 2 Φ B ∂y 2 + D ⊥ ∂ 2 Φ B ∂x 2 + D xy ∂ 2 Φ A ∂y∂x , (1.88)
where the coefficients phase-diffusion have the following expressions

D = ξ 2 0 τ 0 1 - 2q 2 g 1 R 2 0 (g 2 1 -β 2 ) ξ 2 0 τ 0 ; D ⊥ = ξ 2 0 τ 0 q k c ; D xy = 2q 2 β R 2 0 (g 2 1 -β 2 ) ξ 2 0 τ 0 2 . (1.89) 1.6. SECONDARY INSTABILITIES 53
Using normal mode decomposition, an eigenvalue problem is derived. The eigenvalues are

σ 3 = -Q 2 1 D + Q 2 2 D ⊥ -Q 1 Q 2 D xy ; σ 4 = -Q 2 1 D + Q 2 2 D ⊥ + Q 1 Q 2 D xy (1.90)
Phase instabilities : case where

Q 1 = Q 2
Considering the case where Q 1 = Q 2 = Q and a long wavelength limit, i.e. Q → 0, the eigenvalues (1.90) reduce to

σ 3 = -Q 2 D + D ⊥ + D xy = ξ 2 0 τ 0 Q 2 2 ξ 2 0 τ 0 q 2 R 2 0 (g 1 + β) - k c + q k c , (1.91
)

σ 4 = -Q 2 D + D ⊥ -D xy = ξ 2 0 τ 0 Q 2 2 ξ 2 0 τ 0 q 2 R 2 0 (g 1 -β) - k c + q k c . (1.92)
Square Eckhaus instability .

The eigenvalue σ 3 (1.91) states that the system is stable provided that :

ǫ > ǫ SE = 3k c + q k c + q ξ 2 0 q 2 , (1.93)
where ǫ SE is the boundary of the square Eckhaus instability. When q << k c , the universal expression ǫ = 3ξ 2 0 q 2 is recovered. In the phase approximation, the eigenvector corresponding to σ 3 is (1, 1). The wavenumbers in the x-and y-directions evolve in the same way. Note that like for two-dimensional rolls, ǫ SE does not depend on the rheological parameters.

Rectangular Eckhaus instability .

The eigenvalue σ 4 (1.92) states that the system is stable provided that :

ǫ > ǫ RI = ξ 2 0 q 2 1 + 2 g 1 + β g 1 -β k c k c + q . (1.94)
When q << k c , we recover the expression given by Holmedal [19] and Hoyle [18], i.e.

ǫ RI = (3g 1 + β)/(g 1 -β
). An eigenvector associated with σ 4 (1.92) is (1, -1), therefore the wavenumbers in the x-and y-directions don't have the same time evolution. Hoyle [18] denoted ǫ RI (1.94) as rectangular Eckhaus instability. An other point of view was given by Holmedal [19]. According to this author, since the eigenvector is (-1, 1) one of the two rolls will grow at the expense of the other at a particular horizontal location. The decreasing mode will be the growing one at another location. This instability is denoted by Holmedal [19] as "long wavelength cross-roll instability". Concerning the influence of the rheological parameters, it can be shown straightforwardly that ǫ RI boundary shrinks with increasing shear-thinning effects.

CHAPITRE 1.

Phase instabilities : case where either Q 1 or Q 2 is zero

In the case where either Q 1 or Q 2 is zero, the phase equations (1.87) and (1.88) reduce to

∂Φ A ∂t = D ∂ 2 Φ A ∂x 2 and ∂Φ B ∂t = D ⊥ ∂ 2 Φ B ∂x 2 .
(1.95)

The eigenvalues σ 3 and σ 4 of phase-diffusion equations (1.95) are then

σ 3 = -Q 2 1 D and σ 4 = -Q 2 1 D ⊥ . (1.96) (1.97)
Zigzag instability .

The eigenvalue σ 4 (1.96) which is independent of the rheological parameters causes an instability of the squares if q < 0, i.e. k < k c . This is similar to the condition for zigzag instability of the rolls.

2D Eckhaus instability .

The eigenvalue σ 3 (1.96) leads to another phase instability boundary given by

ǫ 2DE = ξ 2 0 q 2 3g 1 -β g 1 -β , (1.98)
which can be considered as a 2D Eckhaus instability [18].

Finally, it can be shown straightforwardly that ǫ RI > ǫ 2DE > ǫ SE when β < g 1 .

Therefore, the stability boundaries of the square pattern are given by the condition k > k c (zigzag instability) and (1.94), i.e. Eckhaus rectangular instability.

(

A ℓp e -i(kxx ℓ +kyyp) , k x , k y = 2π L -N 2 , ..., N 2 -1 (1.99)
Leaving the time stepping in Fourier space, gives the following system of ODEs

d dt Âkxky = ǫ τ 0 - ξ 2 0 τ 0 k 2 x + k x k 2 y k c + k 4 y 4k 2 c
Âkxky -N 1,kxky (A, B) , (1.100)

d dt Bkxky = ǫ τ 0 - ξ 2 0 τ 0 k 2 y + k y k 2 x k c + k 4 x 4k 2 c Bkxky -N 2,kxky (A, B) , (1.101)
with Fourier transformed initial conditions. The nonlinear terms N 1,kxky and N 2,kxky are evaluated in physical space and then transformed to Fourier space :

N 1,kxky (A, B) = -g 1 F |A(x, y, t)| 2 A(x, y, t) -βF |B(x, y, t)| 2 A(x, y, t) , (1.102)
and

N 2,kxky (A, B) = -g 1 F |B(x, y, t)| 2 B(x, y, t) -βF |A(x, y, t)| 2 B(x, y, t) , (1.103)
where F designates the 2D discrete Fourier transform. For the temporal discretization, the time domain [0, t max ] is discretized with equal time step of width ∆t as t m = m∆t, m = 0, 1, 2.... Exponential Time Differencing method of second order (ETD2) proposed by Cox and Matthews [43] is used. Additional details can be found in [44]. The pseudo-spectral method is implemented in Matlab. Finally, to check the convergence, several simulations are carried out with increasing numbers of grid points and refining the time step. The CHAPITRE 1.

stability properties of ETD2 are given in Appendix 1.B.

Numerical results

Instability of a roll solution This is illustrated by Figs. 1.17 and 1.18 where the planform function f (x, y, t) is represented. At point P 3 , the initial state is a uniform roll solution with q = kk c = 0.4, ǫ = 0.05 and α = 1.6 10 -4 . It is in CR-stable and in Eckhaus unstable region, where the pattern wavelength is too short. A small random perturbation is added in the x-and y-directions.

Figure 1.17 shows the time evolution of the convective pattern. The system eliminates two pairs of rolls in order to augment its wavelength. After the local elimination of the wavelength, the system readjusts through a process of phase diffusion. The system reaches a wavenumber inside the Eckhaus stable region. Actually, the final wavenumber is close to k c . Note that, unlike the one dimensional situation, where the defect exists only for an instant while a pair of rolls is created or eliminated, in two-dimensional situation, the defects persist for some time, as shown by Figs. 1.17(b) and (c). At point P 4 , the initial state is a uniform roll with q = 0.5. In this case, the system is CR and Eckhaus unstable. Figure 1.18 shows the time evolution of the structure. In the first stage, a competition between CR and Eckhaus instability mechanisms is observed, before a phase diffusion process.

Instability of a square solution The Nusselt number is roughly the same for rolls and squares and the stable structure has the highest Nusselt number, in agreement with the maximum heat transfer principle.

Subsequently, the stability of modulated rolls and squares with respect to inhomogeneous spatial perturbation is analyzed. The influence of shear-thinning effects is clearly highlighted. In the case of modulated rolls, contrary to the Newtonian case where cross-roll instability is always dominant, except at low Prandtl number, it is shown that for a non-Newtonian shear-thinning fluid, this instability prevails only when α is less than a critical value denoted α CR-E . Under this condition, the domain of stable rolls is bounded by zigzag instability for k < k c and by cross-roll instability for k > k c . Furthermore, the marginal cross-roll curve enlarges with increasing shear-thinning effects. For sufficiently strong shear-thinning effects, here α > α CR-E , Eckhaus instability which is independent of the rheology becomes dominant. The stable rolls are bounded by zigzag and Eckhaus instabilities.

In the case of modulated squares, i.e. α < α S-R , it is observed that the rectangular instability is dominant and the width of the stable wavenumbes band decreases as α increases.

The time evolution of the convective pattern initially in the unstable part of the stability-diagram is obtained from the numerical computation of the amplitude equations.

The instability mechanisms are illustrated, and for all the cases considered, it is observed at the final state the structure reaches a wavenumber very close to the critical value.

CHAPITRE 1.

In this study, the variation of the fluid properties and particularly the viscosity with temperature is not taken into account. Generally, non-Newtonian fluids are highly viscous and thermodependent. The thermodependency of the fluid properties leads to hexagonal patterns at the onset [START_REF] Busse | The stability of finite amplitude cellular convection and its relation to an extremum principle[END_REF], [START_REF] Hoard | Experiments on the cellular structure in Bénard convection[END_REF], [START_REF] Darbouli | Natural convection in shear-thinning fluids : Experimental investigations by MRI[END_REF]. Analysis of the stability of this convective pattern is the next step of our work, dealing with the influence of the rheology on the pattern selection.

An other direction of the present study is the determination of the Busse Balloon for highly shear-thinning fluids. This is particularly interesting, since the lower part of this balloon is delimited by zigzag and Eckhaus boundaries. increases with increasing shear-thinning effects.

The case where (q + k c , ǫ), is zigzag unstable (q < 0) and Eckhaus stable is illustrated by Fig. 1.23. The eigenvalues σ 1 and σ 2 are negative and their variations on Q 1 and Q 2 are similar to those in Fig. 1.22 and are therefore not represented. The eigenvalue σ 3 is negative and may be associated with the Eckhaus stability, whereas σ 4 is positive and is associated with the zigzag instability which is maximum along the axis.

1.B Stability of the ETD2 scheme : Extension to two coupled Ginzburg-Landau equations Using the second order exponential time differencing scheme for the time discretization leads to :

u n+1 = u n e c∆t + F n (1 + c∆t)e c∆t -1 -2c∆t c 2 ∆t + F n-1 1 + c∆t -e c∆t c 2 ∆t (1.105)
To evaluate the stability domain of this scheme, we have adopted the same approach as in [43] and in [START_REF] Ashi | Numerical methods for stiff systems[END_REF]. We suppose that there is a fixed point u 0 , so that cu 0 + F (u 0 ) = 0.

Linearizing about this fixed point leads to

du(t) dt = cu(t) + λu(t), (1.106)
where u is now the perturbation to u 0 and λ = F ′ (u 0 ). The fixed point, u 0 , is stable provided that Re(c + λ) < 0.

It can be shown straightforwardly that the fixed points of the ETD2 scheme are the same as those of the ODE (1.104) : consequently, the linear stability analysis of the ETD2 scheme can be performed by replacing F n by λu n in (1.105). A recurrence relation involving u n+1 , u n and u n-1 is obtained :

u n+1 = u n e c∆t + λu n (1 + c∆t)e c∆t -1 -2c∆t c 2 ∆t + λu n-1 1 + c∆t -e c∆t c 2 ∆t (1.107)
Defining r = u n+1 /u n , x = λ∆t, y = c∆t, the following quadratic equation for the factor r by which the solution is multiplied after each step is derived :

y 2 r 2 -r(y 2 e y + x[(1 + y)e y -2y -1]) + (e y -1 -y)x = 0 (1.108)
r is the factor by which the solution is multiplied after each step so ETD2 scheme is stable provided that |r| < 1. In general, both c and λ are complex and consequently, the stability domain of the ETD2 scheme is four dimensional. To simplify, we choose to determine the stability region in the real plane (Re(x), Re(λ)). The boundaries of this domain are obtained for the values r=1 and r=-1 and correspond to the curves :

y = -x and x = -y 2 (1 + e y ) (y + 2)e y -3y -2
.

(1.109)

1.B.2 Extension to two coupled Ginzburg-Landau equations

In our case, we have the following set of ODE for each couple (k x , k y ) :

d Âkx,ky (t) dt = ǫ τ 0 - ξ 2 0 τ 0 k 2 x + kxk 2 y kc + k 4 y 4k 2 c Âkx,ky (t) + N 1,kxky (1.110) d Bkx,ky (t) dt = ǫ τ 0 - ξ 2 0 τ 0 k 2 y + kyk 2 x kc + k 4 x 4k 2 c
Bkx,ky (t) + N 2,kxky (1.111)

1.B. STABILITY OF THE ETD2 SCHEME

67

where Âkx,ky (t) = F(A(x, y, t)) and Bkx,ky (t) = F(B(x, y, t)). F(.) designates the 2D discrete Fourier transform. We add small perturbations a and b to the initial stationary solutions A 0 and B 0 . Replacing A = A 0 + a and B = B 0 + b in the former equations and after linearization, we get :

dâ kx,ky (t) dt = ǫ τ 0 - ξ 2 0 τ 0 k 2 x + kxk 2 y kc + k 4 y 4k 2 c -2g 1 |A 0 | 2 -β |B 0 | 2 âkx,ky (t) (1.112) -g 1 F A 2 0 ā kx,ky -βF A 0 B0 b kx,ky -βF A 0 B 0 b kx,ky d bkx,ky (t) dt = ǫ τ 0 - ξ 2 0 τ 0 k 2 y + kyk 2 x kc + k 4 x 4k 2 c -2g 1 |B 0 | 2 -β |A 0 | 2 bkx,ky (t) (1.113) -g 1 F B 2 0 b kx,ky -βF B 0 Ā0 a kx,ky -βF (B 0 A 0 ā) kx,ky
Since we consider infinitesimal perturbations, we assume that we can simplify the former equations according to :

dâ kx,ky (t) dt = ǫ τ 0 - ξ 2 0 τ 0 k 2 x + kxk 2 y kc + k 4 y 4k 2 c -2g 1 |A 0 | 2 -β |B 0 | 2 âkx,ky (t) (1.114) d bkx,ky (t) dt = ǫ τ 0 - ξ 2 0 τ 0 k 2 y + kyk 2 x kc + k 4 x 4k 2 c -2g 1 |B 0 | 2 -β |A 0 | 2 bkx,ky (t) (1.115)
We deduce that the values of c and λ in this approximation are given by :

c a (k x , k y ) = ǫ τ 0 - ξ 2 0 τ 0 k 2 x + kxk 2 y kc + k 4 y 4k 2 c ; λ a (k x , k y ) = -2g 1 |A 0 | 2 -β |B 0 | 2 (1.116) c b (k x , k y ) = ǫ τ 0 - ξ 2 0 τ 0 k 2 y + kyk 2 x kc + k 4 x 4k 2 c ; λ b (k x , k y ) = -2g 1 |B 0 | 2 -β |A 0 | 2 .(1.117)
We have checked that for our time step ∆t and for each couple (k x , k y ), the points (x a , y a ) =

(c a (k x , k y )∆t, λ a (k x , k y )∆t) and (x b , y b ) = (c b (k x , k y )∆t, λ b (k x , k y )∆t)

Introduction

Rayleigh-Bénard convection in a fluid layer which viscosity depends on temperature is an important problem, because of its interest in many industrial and natural systems such the convection in the Earth mantle and in magma chambers. The viscosity, µ, is usually much lower near the heated lower plate than near the cold upper plate. This spatially varying viscosity modifies the onset of convection. Furthermore, the variation of µ with the temperature, T , causes an additional nonlinear coupling between the velocity and the temperature fields and breaks the up-down reflection symmetry with respect to the mid-plane of the fluid layer. These features will affect the onset of convection and the selection of the pattern convection.

For Newtonian fluids, the effect of a temperature dependent viscosity on the onset of convection was first studied by Palm [1] in the case of free-free boundary conditions. [1] assumed that the kinematic viscosity ν varies as ν = ν 1 + ∆ν cos (b (T -T 1 )) where ∆ν is the difference in the viscosity between the top and the bottom boundaries, b a constant and T 1 is the temperature at the bottom of the fluid layer. In this analysis, it is required that ∆ν/ν 1 << 1. It is found that the critical Rayleigh number Ra c defined with the average viscosity ν 0 as well as the critical wavenumber k c decreases with increasing the viscosity variation ∆ν. They differ by O (∆ν/ν 0 ) 2 from that obtained with constant viscosity. The decrease of Ra c and k c with increasing ∆ ν was confirmed by Stengel et al [2] in Free-Free and Rigid-Rigid boundary conditions when a cosine law is used for µ(T ). Busse and Frick [3] assumed, for numerical convenience, a linear dependence of the viscosity on temperature. The onset of convection is determined in the case of rigid boundary conditions.

The variation of Ra c and k c as a function of the viscosity ratio r = µ max /µ min is quite similar to that obtained by [1] using cosine law for µ(T ). As pointed out by [3], for cosine and linear functions µ(T ), the viscosity at the midplane equals to the average viscosity of the static layer, this is why Ra c decreases with increasing r. However, if an exponential viscosity variation is used, the average viscosity exceeds the value used in the definition of Ra c . In this case, the critical Rayleigh number Ra c increases, reaches a maximum of Ra c ≈ 2200 at a viscosity ratio r ≈ 3000 and then decreases [2]. This result was confirmed by White [4]. It can be explained by a simple physical argument based on the idea that convection begins first in the sublayer with maximum Rayleigh-number. Actually, for a large viscosity contrast, the convection is confined to the sublayer near the hot boundary, and a stagnant zone develops near the cold (top) boundary. [2,5,6]. Whereas, for cosine and linear laws µ(T ), the convection occurs throughout the entire fluid layer. The onset of two-dimensional convection with strongly temperature dependent viscosity has been also considered by Bottaro et al. [7], assuming Arrhenius law. In this case, the viscosity ratio depends on the temperature difference across the fluid layer and on the temperature level, while for exponential law, the viscosity ratio depends only on the temperature difference. [7] found that depending on the reference temperature, the dependence of the critical Rayleigh-number Ra c on the viscosity ratio across the layer, may have one of the two behaviors described previously. Either, Ra c decreases with increasing the viscosity ratio as predicted by [1] and [3], or Ra c increases initially with increasing the viscosity ratio, reaches a maximum and then decreases as predicted by [2]. Actually, there are two controlling factors that play opposing roles. The reduced thickness of the active layer on one hand requires a larger Rayleigh number for the onset of convection. On the other hand, the fluid layer near the heated wall is less stable because of the decrease of the viscosity.

The effect of weakly temperature-dependent viscosity on the planform near the critical conditions, has been studied by [1], [8], [9]. They found that near the onset, the convection occurs in the form of hexagons, when r is small. The occurrence of hexagonal pattern

CHAPITRE 2.
can be explained by the triadic wavevector interactions enabled by the quadratic term in the amplitude equation. If the Rayleigh number is increased slightly, hexagons become unstable to rolls solution. A formula for the range of Rayleigh numbers where the hexagons are stable was proposed by [8]. It is shown that the extent of the interval where the hexagons are stable is proportional to the inverse of the fluid depth to the sixth power.

Experimental planform studies performed by Hoard et al [10] and Stengel [2] confirm qualitatively the results of the weakly nonlinear theory at low values of r. The quantitative disagreement may be due partly to the influence of the lateral walls and also to the fact that the Rayleigh numbers considered are probably outside the domain of validity of the theory.

Besides rolls and hexagons, a new planform in the form of squares was observed when the viscosity contrast between upper and lower boundaries exceed a value of order ten, [2] and [4]. The planform selection problem between rolls and squares was analyzed by [3] with the assumption that the viscosity varies linearly with temperature. They found that near the critical conditions, rolls are preferred for low values of r, but squares are preferred for large values of r. The change from rolls to squares occurs at r ≈ 2. Jenkins [11] used a weakly nonlinear method to investigate the stability of squares. In the case of a linear variation of the viscosity with temperature, he found that the transition from rolls to squares occurs at r ≈ 3.2. The disagreement with [3] was not clarified in the literature. For exponential fluids, Jenkins [11] found that the transition occurs at r ≈ 3.

Above onset, the wavenumber k of the pattern is selected from a continuous band of allowed wavenumbers, that is bounded by secondary instabilities. In [12] and [13], the authors identified the mechanisms of secondary instability that restrict the region of stable straight convection rolls in Rayleigh-Bénard convection. The nature of secondary instabilities in more complex patterns such as squares or hexagons is not as well studied as rolls. In the case of hexagonal pattern, it is shown that the secondary instability is induced by long wavelength modulation of the phase of the pattern. In Bénard-Marangoni problem, estimates of the size and shape of stable band of wavenumbers have been made by Echebarria and Pérez-Garcia [14] and Young and Riecke [15] using amplitude equations.

To our knowledge, the influence of viscosity temperature dependence of a non-Newtonian fluid on the pattern selection and the phase instabilities that delimit the stable band of wavenumbers in Rayleigh-Bénard convection problem has not been considered in the literature. Shear-thinning fluids are the most kind of non-Newtonian fluids and are characterized by a decreasing viscosity with increasing shear rate. Furthermore, their viscosity may vary strongly with temperature. The decrease of the viscosity with temperature modifies the flow structure. The shear rate is higher near the heated wall. This leads to a further decrease of the viscosity near the heated wall because of the shear-thinning behavior of the fluid. Hence, the degree of asymmetry with respect to the midplane becomes more significant.

The objective of the present paper is to investigate the influence of shear-thinning effects on the pattern selection and their stability, taking into account the dependence of the viscosity on temperature. This paper is organized as follows. We start with the governing equations in Sec. 2.3. The linear stability analysis is presented in Sec. 2.4. In Sec. 2.5, the weakly nonlinear analysis using a multiple scale method is presented, the amplitude equations for hexagons are derived and the different coefficients are determined as a function of shear-thinning effects and the viscosities ratio. First homogeneous hexagons and rolls are considered and their relative stability is discussed in Sec. 2.6. In the Section 2.7, the limit value of the viscosity ratio above which rolls become unstable to square is determined as a function of shear-thinning effects. The stability of hexagons with respect to long wavelength perturbations is addressed in Sec. 2.8 in which the phase equations are derived. Finally, a brief summary of the results is given in Sec.2.9.

Basic equations

Hereafter, quantities with hats are dimensional quantities. We consider a layer of shearthinning fluid of depth d confined between two impermeable horizontal plates, infinite in extent, which are perfect heat conductors. The bottom and top plates are kept at constant temperatures, respectively T0 +δ T /2 and T0 -δ T /2, with δ T > 0. The fluid has density ρ, thermal diffusivity κ, thermal expansion coefficient β and viscosity μ0 at zero shear-rate.

In the absence of convection, the stationary heat conducting state is described by

û = 0 , d P dẑ = -ρ 0 ĝ 1 -β Tcond -T0 , Tcond -T0 = δ T 2 1 - 2ẑ d , (2.1)
where û is the fluid velocity, P the pressure and T0 the mean of the boundary temperatures. The z-axis is directed upwards, with its origin located at the bottom plate. The stability of the hydrostatic solution is considered by introducing temperature and pressure perturbation as well as a fluid motion. Using the units d2 /κ, d, κ/ d and δ T for time, length, velocity and temperature, the dimensionless perturbation equations are :

∇ • u = 0 , (2.2) 1 P r ∂u ∂t + (u • ∇) u = -∇p + Ra θ e z + ∇ • τ , (2.3)
∂θ ∂t + u • ∇θ = u • e z + ∇ 2 θ . (2.4)
Here, e z denotes the unit vector in the vertical direction, p(x, t) and θ(x, t) represent the pressure and temperature deviations from their values in the conductive state. The Boussinesq approximations are taken into account, i.e., the variation of the density is neglected except in the buoyancy term. Denote (x, y, z) the components of the position vector x, and (u, v, w) the components of the velocity vector u. The Rayleigh number Ra and the Prandtl number P r are Ra = ρ0 ĝ β δ T d3 κ μ0 ; P r = μ0 ρ0 κ .

(2.5)

At the second order Taylor expansion of 1 + λ 2 Γ nc-1

2 , a relevant rheological parameter, i.e., the 'degree of shear-thinning' appears :

α = dµ dΓ Γ=0 = 1 -n c 2 λ 2 .
(2.11)

Boundary conditions

For the velocity field, no-slip boundary conditions (NSBC) are considered. Isothermal conditions are used at slabs. For the temperature deviation θ, the thermal conductivity of the boundaries is assumed much larger than that of the fluid, so that slabs' temperature remains 'fixed'. The boundary conditions are :

θ = u = v = w = 0 on z = 0, 1.
(2.12)

Reduction : elimination of the pressure

Applying twice the curl to momentum equations (2.3) and using the continuity equation, we get the following evolution equations for the velocity components w, u and v :

1 P r ∂ ∂t ∆w = 1 P r ∂ 2 ∂y∂z N (v) + ∂ 2 ∂x∂z N (u) -∆ H N (w) + Ra∆ H θ + µ b ∆ 2 w + 2 dµ b dz ∆ ∂w ∂z + d 2 µ b dz 2 ∂ 2 w ∂z 2 -∆ H w + (2.13) ∆ H N V z - ∂ 2 ∂x∂z N V x - ∂ 2 ∂y∂z N V y , 1 P r ∂ ∂t ∆ H u + ∂ 2 w ∂x∂z = 1 P r ∂ 2 ∂x∂y N (v) - ∂ 2 ∂y 2 N (u) + µ b ∆ ∆ H u + ∂ 2 w ∂x∂z + dµ b dz ∂ ∂z ∆ H u + ∂ 2 w ∂x∂z - ∂ ∂y N V z (2.14) 1 P r ∂ ∂t ∆ H v + ∂ 2 w ∂y∂z = 1 P r ∂ 2 ∂x∂y N (u) - ∂ 2 ∂x 2 N (v) + µ b ∆ ∆ H v + ∂ 2 w ∂y∂z + dµ b dz ∂ ∂z ∆ H v + ∂ 2 w ∂y∂z + ∂ ∂x N V z , (2.15)
where the 'horizontal Laplacian' is

∆ H = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 .
The nonlinear inertial terms N (•) and nonlinear viscous terms N V x are defined by

N (•) = (u • ∇) (•) ; N V x = [∇ • ((µ -µ b ) γ)] • e x , (2.16
)

CHAPITRE 2.
similarly for N V y and N V z . The boundary conditions are

θ = w = ∂w ∂z = u = v = 0 at z = 0, 1 (2.17)
The condition ∂w ∂z = 0 comes from the incompressibility condition.

In a matrix notation, the above system can be written formally as

M ∂Ψ ∂t = LΨ + N I + N V , (2.18)
where Ψ = (w, u, v, θ) T , the operators M , L, N I and N V represent the weight matrix, the linear operator, the nonlinear inertial operator and the nonlinear viscous operator respectively. The nonlinear operators can also be decomposed as

N I = [N I w , N I u , N I v , N I θ ] t and N V = [N V w , N V u , N V w , 0] t . (2.19)
2.4 Linear stability analysis

Onset of convection : direct eigenvalue problem

In the linear theory, u and θ are assumed infinitesimal. By neglecting the nonlinear terms in (2.13)-(2.15) and (2.4), one obtains the linear problem :

1 P r ∂ ∂t ∆w = µ b ∆ 2 w + 2 dµ b dz ∆ ∂w ∂z + d 2 µ b dz 2 ∂ 2 w ∂z 2 -∆ H w + Ra∆ H θ , (2.20) 1 P r ∂ ∂t ∆ H u + ∂ 2 w ∂x∂z = µ b ∆ ∆ H u + ∂ 2 w ∂x∂z + dµ b dz ∂ ∂z ∆ H u + ∂ 2 w ∂x∂z (2.21) 1 P r ∂ ∂t ∆ H v + ∂ 2 w ∂y∂z = µ b ∆ ∆ H v + ∂ 2 w ∂y∂z + dµ b dz ∂ ∂z ∆ H v + ∂ 2 w ∂y∂z (2.22) ∂θ ∂t = w + ∆θ. (2.23)
At this order, no non-Newtonian effects enter the problem and the thermodependency appears through the viscosity profile of the base state µ b (z). The disturbance quantities w, u, v, θ are assumed periodic and of the form :

(w, u, v, θ) = (F 11 (z), U 11 (z), V 11 (z), G 11 (z)) f (x, y) exp (st) (2.24)
with f (x, y) = exp (ik x x + ik y y), k = (k x , k y , 0) the horizontal wavenumber and s = s r +is i 2.4. LINEAR STABILITY ANALYSIS 81 a complex number. This leads to the differential equations

s P r -1 D 2 -k 2 F 11 = µ b D 2 -k 2 2 F 11 + 2Dµ b D 2 -k 2 DF 11 + (2.25) D 2 µ b D 2 + k 2 F 11 -k 2 RaG 11 , s P r -1 -k 2 U 11 + ik x DF 11 = µ b D 2 -k 2 -k 2 U 11 + ik x DF 11 + (2.26) Dµ b -k 2 DU 11 + ik x D 2 F 11 s P r -1 -k 2 V 11 + ik y DF 11 = µ b D 2 -k 2 -k 2 V 11 + ik y DF 11 + (2.27) Dµ b -k 2 DV 11 + ik y D 2 F 11 s G 11 = F 11 + (D 2 -k 2 )G 11 , (2.28)
with D the derivative with respect to z and k the norm of the vector k. The boundary conditions are with F 11 = 0. For this class of eigenmodes, it can be shown that the eigenvector. It can be written

F 11 = DF 11 = U 11 = V 11 = G 11 = 0 at z = 0, 1 . ( 2
s i = 0 and s r = - k 2 µ b |U 11 | 2 + |V 11 | 2 + µ b |DU 11 | 2 + |DV 11 | 2 |U 11 | 2 + |V 11 | 2 < 0. ( 2
s M • X 11 = L • X 11 . (2.31)
It is easy to show that the principle of exchange of stability still holds, i.e. s i = 0, when the viscosity profile is not uniform. Since any multiple of the eigenvector X 11 is also a solution of (2.31) X 11 has to be normalized. We have adopted the same normalization as in [18] :

G 11 (z = 1/2) = 1. ( 2

.32)

A spectral Chebyshev method is used to determine the critical Rayleigh number and the critical wave number [18]. The marginal stability curve Ra(k) is obtained by the condition s(Ra, k) = 0. Using 20 Chebyshev polynomials, the first eigenvalue, i.e. that for which the real part is the largest, is calculated with an accuracy of 10 increases rapidly. In this regime, the convection is governed by a sublayer that is more unstable than the full layer [2].

The results of [3] can be reproduced by a linear approximation of the exponential law (2.8). The viscosity contrast between the top and the lower boundaries will be reinforced by the shear-thinning effects. The temperature perturbation becomes more confined near the heated wall. As expected, when c = 0, the eigenfunctions, F 11 (z) and G 11 (z), are symmetric with the respect to the midplane of the fluid layer.

Adjoint eigenvalue problem : Adjoint mode

In the analysis developed in the next section, it will be necessary to eliminate secular terms in non homogeneous differential equations, i.e. the solvability condition must be applied. It is therefore necessary to determine the linear adjoint of the direct problem at the critical conditions. For vector fields f and g, one defines an inner product between two vector functions f (z) and g(z) by

f , g = 1 0 f * • gdz , (2.33)
where f * is the complex conjugate of f . To the direct eigenvalue problem (2.31) corresponds the following adjoint problem

s M + • X ad = L+ • X ad with X ad = (F ad , G ad ) , (2.34)
where the adjoint operators M + and L + are defined by

X ad , M • X 11 = M + • X ad , X 11 , X ad , L • X 11 = L+ • X ad , X 11 , (2.35)
where X 11 fulfills the 'linear' boundary conditions (2.29). By integrating by part we get the linear adjoint problem and the corresponding boundary conditions

s P r -1 D 2 -k 2 F ad = µ b D 2 -k 2 2 F ad + 2Dµ b D 2 -k 2 )DF ad + D 2 µ b D 2 + k 2 F ad + G ad , (2.36
)

s G ad = -k 2 Ra F ad + D 2 -k 2 G ad , (2.37)
with

F ad = 0 , DF ad = 0 , G ad = 0 at z = 0, (2.38a
)

F ad = 0 , DF ad = 0 , G ad = 0 at z = 1. (2.38b) CHAPITRE 2.
The solution of these equations is obtained using the same method as for the direct problem. Similarly, the normalization adopted for the adjoint mode is

G ad (z = 1/2) = 1.
(2.39)

At Ra = Ra c , the so-called adjoint critical mode does not depend on the Prandtl number.

Amplitude equations in a hexagonal lattice

The i.e. there is also a pattern degeneracy. In this section, we consider the case where the wavevectors lie on a hexagonal lattice

(w, θ) = 3 p=1 A p (F 11 , G 11 ) exp (ik p • r) + c.c. + h.o.t. , (2.40)
where "c.c." denotes the complex conjugate of the prior expression and "h.

Multiple scales method

As the Rayleigh number is increased above the onset Ra c , the growth-rate of the perturbation is positive for any wavenumber within a band √ ǫ around the critical wavenumber, where ǫ = (Ra -Ra c ) /Ra c is the distance from the onset. Indeed, Taylor expansion of the dispersion relation near its maximum shows that s ∝ ǫ, kk c ∝ √ ǫ. For ǫ > 0, emergent patterns are described by an infinite sum of unstable modes (in a continuous band) of the form exp

ǫt τ 0 exp (ik c x) exp i √ ǫx ξ 0 .
Here, τ 0 is the characteristic time for the instability to grow and ξ 0 is the coherence length. For small ǫ, we can separate the dynamics into fast eigenmodes and slow modulation of the form exp ǫ t τ 0 exp i √ ǫx ξ 0 .

A similar reasoning can be done for the y-direction.

Let us denote δ = √ ǫ. The multiple-scales approach is used to obtain the amplitude equation, which describes the slow temporal and spatial variation of the field variables. The

∂ ∂t -→ ∂ ∂t + δ 2 ∂ ∂T , ∂ ∂x -→ ∂ ∂x + δ ∂ ∂X , ∂ ∂y -→ ∂ ∂y + δ ∂ ∂Y , ∂ ∂z -→ ∂ ∂z . (2.42)
The fast spatial variables vary on the order of a typical wavelength. The slow variables describe, the temporal and the spatial modulations of these fast variables. Furthermore, as the marginal mode is stationary, then

∂ ∂t -→ δ 2 ∂ ∂T . (2.43)
The solution of the nonlinear problem in the neighborhood of the critical conditions, corresponding to the onset of convection is developed with respect to the parameter δ by

Ψ = δΨ (1) + δ 2 Ψ (2) + δ 3 Ψ (3) + O(δ 4 ) , (2.44
)

Ra = Ra c + δRa (1) + δ 2 Ra (2) + O(δ 3 ) . (2.45)
The Taylor expansion is also applied to the operators

M = M (0) + δM (1) + O δ 2 , (2.46)
L = L (0) + δL (1) + δ 2 L (2) + O δ 3 , (2.47)
N I = δ 2 N I (2) + δ 3 N I (3) + O δ 3 , (2.48
)

N V = δ 2 N V (2) + δ 3 N V (3) + O δ 3 .
(2.49)

The explicit expressions of M , L, N I and their sub-scales are given in appendices. The expressions of N V and its sub-scales are too lengthy, and thus are not shown. ). After ordering according to the power of δ, a sequence of systems of equations is obtained. In the following, the first three orders are determined.

Derivation of the Ginzburg-Landau equation

Solution at order δ

At the first order of δ, the linearized problem is obtained

L (0) Ψ (1) = 0 .
(2.50)

The system (2.50) corresponds to the linear problem discussed in Sec.2.4. However, now Ψ (1) is also a function of the slow variables X, Y and T . These variables do not appear in the linear stability analysis section. For hexagon patterns, the first order solution

Ψ (1) = CHAPITRE 2.
w (1) , u (1) , v (1) , θ (1) t is

w (1) = F 11 (z) [A 1 exp (ik 1 • r) + A 2 exp (ik 2 • r) + A 3 exp (ik 3 • r)] + c.c. , (2.51) θ (1) = G 11 (z) [A 1 exp (ik 1 • r) + A 2 exp (ik 2 • r) + A 3 exp (ik 3 • r)] + c.c. , (2.52) u (1)
H = DF 11 k 2 ∇ Hx [A 1 exp (ik 1 • r) + A 2 exp (ik 2 • r) + A 3 exp (ik 3 • r)] + c.c. , (2.53)
where ∇ Hx denotes the horizontal gradient for the fast variables, u H = (u, v) the horizontal velocity components, and A p the amplitude of the perturbation :

A p = A p (X, Y, T ) p = 1, 2, 3 . (2.54)
Solution at order δ 2

At the next order δ 2 , we have The calculation leads to a result of the form

L (0) Ψ (2) = -L (1) Ψ (1) -N I (2) -N V (2) . ( 2
A * 2 A * 3 1 P r 1 0 F ad 2DF 11 D 2 F 11 + F 11 D 3 F 11 -3k 2 c F 11 DF 11 dz - A * 2 A * 3 1 0 F ad [N V w ] (2) A * 2 A * 3 dz - 2 i (k 1 • ∇ HX ) A 1 2 1 0 µ b D 2 F 11 -k 2 c F 11 F ad dz + 2 1 0 dµ b dz dF 11 dz F ad dz - 2 i (k 1 • ∇ HX ) A 1 - 1 0 d 2 µ b dz 2 F 11 F ad dz + Ra c 1 0 G 11 F ad dz - k 2 c Ra 1 A 1 1 0 G 11 F ad dz = 0 .
Ra 1 A 1 + b A * 2 A * 3 = 0 . (2.57)
Again, two other similar relations are obtained by circular permutation. These expressions allow to determine the solution at the second order.

Solution at order δ 3

At this order, we obtain the equation for the evaluation of Ψ (3) :

L (0) Ψ (3) = M (0) ∂Ψ (1) ∂T -L (1) Ψ (2) -L (2) Ψ (1) -N I (3) -N V (3) (2.58)
The solvability condition at this order, yields an equation for Ra 2 . To obtain the amplitude equations at cubic order, we use equation (2.45), combined with ǫ = (Ra -Ra c )/Ra c , the departure from the linear threshold. We have

ǫA 1 = δ Ra c Ra (1) + δ 2 Ra c Ra (2) A 1 . (2.59)
We substitute in (2.59), Ra (1) and Ra (2) by their expressions derived from the solvability conditions at orders δ 2 and δ 3 . Finally, returning to the fast variable

δA j (X, Y, T ) = A ′ j (x, y, t), ∂ ∂X = 1 δ ∂ ∂x , .
.., we obtain

∂A 1 ∂t = ǫα l A 1 + α d (k 1 • ∇ Hx ) 2 A 1 + γA * 2 A * 3 - (2.60) g 1 |A 1 | 2 A 1 -g 2 |A 2 | 2 + |A 3 | 2 A 1 + iβ 1 [A * 2 (k 3 • ∇ Hx ) A * 3 + A * 3 (k 2 • ∇ Hx ) A * 2 ] + iβ 2 [A * 2 (k 2 • ∇ Hx ) A * 3 + A * 3 (k 3 • ∇ Hx ) A * 2 ] .
Companion equations for A 2 and A 3 are obtained by subindex permutation.

Following [14], it is useful to express the derivatives in equation (2.61) in terms of unitary vectors of the corresponding mode. One obtains

∂A 1 ∂t = ǫ τ 0 A 1 + ξ 2 0 τ 0 (n 1 • ∇ Hx ) 2 A 1 + γA * 2 A * 3 - (2.61) g 1 |A 1 | 2 A 1 -g 2 |A 2 | 2 + |A 3 | 2 A 1 + iα 1 [A * 2 (n 3 • ∇ Hx ) A * 3 + A * 3 (n 2 • ∇ Hx ) A * 2 ] + iα 2 [A * 2 (τ 3 • ∇ Hx ) A * 3 + A * 3 (τ 2 • ∇ Hx ) A * 2 ] , CHAPITRE 2.
where n i is the unitary vector in the direction of k i and τ i orthogonal to n i . As indicated by [19,20], there is no Lyapunov functional for Eqs (2.61) and it is possible that the steady state cannot be reached. In contrast, when α 1 and α 2 vanish, a Lyapunov functional can be written down in a straightforward way.

The characteristic time for the instability to grow τ 0 and the coherence length ξ 0 are shown in figure 2.3 as a function of r and for different values of P r. The coefficient γ arises from non Oberbeck-Boussinesq effects. It increases with increasing the viscosity ratio, since γ ∝ c = ln(r). The coefficient γ, also increases with the Prandtl number as it is shown in figure 2.4. From P r ≈ 50, there is no significant effect of P r. The coefficient g 1 refers to the self-saturation coefficient and g 2 the cross-saturation coefficient. It can be shown straightforwardly that g 1 and g 2 can be written as the sum of two contributions. The first one (g N 1 , g N 2 ) similar to that obtained for a Newtonian fluid arises from the nonlinear inertial terms and the thermodependency of the viscosity. The second contribution (g V 1 , g V 2 ) arises from the nonlinear variation of the viscosity with the shear-rate :

g 1 = g N 1 + g V 1 with g V 1 = -αg N N 1 , (2.62)
and similarly for g 2 . Variations of

g N 1 , g N N 1 , g N 2 and g N N 2
as a function of the viscosity ratio for different values of P r are displayed in Fig 2 .6. The coefficients g 1 and g 2 increase significantly with Prandtl number up to P r = 50. Whereas, the dependency of g 1 and g 2 on r is quiet modest. The coefficients α 1 and α 2 are real. The term with α 1 accounts for distortions in the directions of rolls and therefore corresponds physically to wavenumber dilatation. The coefficient α 1 is positive and takes values of the same order as γ. Note also that α 1 vanishes when r = 1, and increases with increasing r. The terms with α 2 account for distortions in the hexagonal form. The coefficient α 2 , is negative and smaller (in absolute value) than α 1 . For the set of coefficients discussed above, the following correlations can be used for 1¡r¡3 and 1 < P r < ∞ : τ 0 = 0.0508 + 0.026P r -1 ; ξ 0 = 0.385 -3.6510 -4 r -1.70810 -4 r 2 , (2.63)

γ = 9.80 -2.9P r -1 + 0.22P r -2 log(r) , g N 1 = 254.3 1.0037 -0.4722P r -1 + 0.1392P r -2 1.0067 -0.0037r -0.002r 2 g N N 1
= 11.86 10 5 1.0038 -0.4808P r -1 + 0.1422P r -2 0.9733 + 0.0128r + 0.0094 r 2 , g N 2 = 375.9 1.0029 -0.3545P r -1 + 0.1020P r -2 1.0091 -0.0047 r -0.0029 r 2 , g N N 2 = 1.343 10 6 1 -0.476P r -1 + 0.146P r -2 0.9756 + 0.0117r + 0.0086r 2 α 1 = 12.5 1 -0.504P r -1 + 0.148P r -2 log(r) , α 2 = -2.1 1.003 -0.238P r -1 + 0.069P r -2 log(r) .

Amplitude instabilities

We study the stability of homogeneous and stationary solutions of equations (2.61) by including slightly off-critical wavenumber in the amplitudes.

(i) Roll solution with a wavenumber slightly off-critical k = k c + q. It is given by

A 1 = R 0 exp (iqx)
, A 2 = A 3 = 0, and any circular permutation with R 0 = ǫξ 2 0 q 2 τ 0 g 1 .

A linear stability analysis of this solution with respect to uniform perturbations

A 1 = (R 0 + r 1 ) exp (iqx 1 ), A 2 = r 2 exp (iqx 2 ) and A 3 = r 3 exp (iqx 3 )
shows that the roll solution is stable when g 2 > g 1 and

ǫ > ǫ r (q) = τ 0 g 1 (γ + 2α 1 q) 2 (g 2 -g 1 ) 2 + ξ 2 0 q 2 .
(2.64) (ii) Hexagon solution is given by three sets of rolls of equal amplitude, A p = H 0 exp (i q x p ) with p = 1, 2, 3 and

H 0 = (γ + 2α 1 q) + (γ + 2α 1 q) 2 + 4 (g 1 + 2g 2 ) (ǫ -ξ 2 0 q 2 )/τ 0 2 (g 1 + 2g 2 ) , (2.65)
called up-hexagons, that correspond to up flow in the center, and

H 0 = -(γ + 2α 1 q) + (γ + 2α 1 q) 2 + 4(ǫ -ξ 2 0 q 2 ) (g 1 + 2g 2 ) /τ 0 2 (g 1 + 2g 2 ) , (2.66)
called down hexagons, that correspond to down flow motion in the center.

Solutions called up-hexagons exist for

ǫ > ǫ a = - (γ + 2α 1 q) 2 τ 0 4 (g 1 + 2g 2 ) + ξ 2 0 q 2 , (2.67)
and are linearly stable for ǫ a < ǫ < ǫ h , with

ǫ h = τ 0 (γ + 2α 1 q) 2 (2g 1 + g 2 ) (g 1 -g 2 ) 2 + ξ 2 0 q 2 (2.68)
Solutions called down hexagons exist for ǫ > ξ 2 0 q 2 and are linearly unstable.

(iii) The "mixed states" given by

A 1 = γ -2α 1 q g2 -g 1 , A 2 = A 3 = (ǫ -ξ 2 0 q 2 )/τ 0 -g 1 A 2 1 g 1 + g 2 (2.69)
and any circulation permutation exist for ǫ > (γ + 2α 1 q) 2 g 1 (g 1g 2 ) 2 + ξ 2 0 τ 0 and are linearly unstable. An example of amplitude stability curves in (ǫ, q) space and the associated bifurcation diagram for q = 0 is given in figures 2.7 and 2.8. Hexagons bifurcate transcritically from the conductive state where they are unstable. Both hexagons and the conductive state are stable in the range ǫ a ≤ ǫ ≤ 0 and both hexagons and rolls are stable in the range ǫ r ≤ ǫ ≤ ǫ h . In this range, rolls and hexagons are linked via a branch of mixed modes which are always unstable.

Variations of ǫ a , ǫ r , ǫ h and (ǫ hǫ r ) as a function of the viscosities ratio, r for different values of the shear-thinning degree are depicted in Fig. 2.9. Overall, the thermodependency of the viscosity favors convection in form of hexagons and their stability whereas shear-thinning effects favor convection in form of rolls and their stability. For instance, in Fig. 2.9(b), the domain of stability of hexagons increases with increasing r and decreases with increasing shear-thinning effects. In the same way, the domain of bistability rolls and hexagons shrinks with increasing α, and increases with increasing the viscosities ratio, i.e. the thermodependency effect.

Competition between rolls and squares

From a limit value of the viscosity ratio r ℓ , rolls become unstable with respect to squares as it was shown theoretically by [3] and [11] and experimentally by [4]. In this section, we investigate the influence of shear-thinning effects on this limit value r ℓ . Here, we consider only competition between perfect rolls and squares without spatial modulation.

In a square lattice, the solution at order δ is The derivation of amplitude equations without spatial terms, for the two modes A 1 and A 2 forming an angle of 90 o follows the same procedure as in §2.5. They are given by

w (1) = F 11 (z) [A 1 exp (ik 1 • r) + A 2 exp (ik 2 • r)] + c.c. , (2.70)
θ (1) = G 11 (z) [A 1 exp (ik 1 • r) + A 2 exp (ik 2 • r)] + c.c. , (
dA 1 dt = sA 1 -g 1 |A 1 | 2 + g 2s |A 2 | 2 A 1 + ... (2.72) dA 2 dt = sA 2 -g 1 |A 2 | 2 + g 2s |A 1 | 2 A 2 + ... (2.73)
Note that Eqs (2.72), (2.73) can be obtained using symmetries introduced by the square lattice [21,22]. A linear stability analysis of stationary rolls and squares solutions of Eqs (2.72), (2.73) allows to show that squares are stable when g 2s < g 1 i.e when the crosscoupling between two orthogonal modes that describe the square pattern is weak enough.

The numerical results are displayed in figure where stability domains of squares and rolls are represented in the plane (α, r) for P r = 10. The curve is a boundary between squares and rolls : below the boundary, squares are stable and above the boundary, rolls are the stable convective pattern. It can be observed that r ℓ increases with increasing shearthinning effects.

Phase instabilities

The stability of the stationary up-hexagons solution (2.65) is investigated by adding small amplitude and phase perturbations of the form of the perturbation respectively. Substitution of (2.74) into (2.61) and linearizing with respect to r p and φ p leads to the following set of equations

A p = (H 0 + r p ) exp [i (qx p + φ p )] , p = 1, 2, 3 , (2
∂r p ∂t = -[2g 1 H 0 + (γ + 2α 1 q)] H 0 r p + ξ 2 0 τ 0 ∂ 2 r p ∂x 2 p + (2.75) [(γ + 2α 1 q) -2g 2 H 0 ] H 0 (r j + r k ) -2qH 0 ξ 2 0 τ 0 ∂φ p ∂x p + α 1 + α 2 √ 3 H 2 0 ∂φ j ∂x j + ∂φ k ∂x k + 2α 2 √ 3 H 2 0 ∂φ j ∂x k + ∂φ k ∂x j ∂φ p ∂t = -(γ + 2α 1 q) H 0 (φ p + φ j + φ k ) + ξ 2 0 τ 0 ∂ 2 φ p ∂x 2 p + 2q H 0 ξ 2 0 τ 0 ∂r p ∂x p + (2.76) α 1 + α 2 √ 3 ∂r j ∂x j + ∂r k ∂x k + 2α 2 √ 3 ∂r j ∂x k + ∂r k ∂x j .
In the following, we consider long-wavelength perturbations, in the region of parameters where hexagons are stable with respect to spatially homogeneous disturbances. In this case, it can be shown that to the leading order, i.e, in the limit of zero wavenumber, the amplitudes r p and the total phase Φ = φ 1 + φ 2 + φ 3 are fastly decaying variables, whereas the remaining two phase modes are neutral (the growth rate vanishes). Actually, a spatially uniform modification of the phase corresponds to a shift of a periodic pattern, and thus is neutral because of translational invariance in space [23]. When the spatial modulation is considered in the long wavelength limit, the two phase modes are no longer neutral, but evolve on long space and time scales [24]. The dynamics are therefore dominated by these two phase modes. Instead of using φ 2 and φ 3 , it is convenient to consider φ x = -(φ 2 + φ 3 )

and φ y = 1 √ 3 (φ 2 + φ 3 ), which are related to the two translational symmetries in the

x and y directions respectively [25,14]. Using the phase approximation, the amplitudes perturbations and the global phase Φ are slaved adiabatically to the slow phase modes.

Adiabatic elimination, reduces then (2.75), (2.76) to

∂φ ∂t = D ⊥ ∇ 2 φ + D -D ⊥ ∇ (∇ • φ) , (2.77)
where the phase vector φ = (φ x , φ y ), D ⊥ and D are the transverse and longitudinal phase diffusion coefficients, given by

D ⊥ = 1 4 ξ 2 0 τ 0 - q 2 2u ξ 2 0 τ 0 2 + H 2 0 8u α 1 - √ 3 α 2 2
(2.78)

D = 3 4 ξ 2 0 τ 0 - q 2 (4u + v) 2uv ξ 2 0 τ 0 2 + H 2 0 8u α 1 - √ 3 α 2 2 - (2.79) H 2 0 α 1 v α 1 + √ 3 α 2 + H 0 q v ξ 2 0 τ 0 3α 1 + √ 3 α 2 , 2.9. CONCLUSION 97 with u = H 2 0 (g 1 -g 2 ) + (γ + 2α 1 q) H 0 , (2.80) v = 2H 2 0 (g 1 + 2g 2 ) -(γ + 2α 1 q) H 0 .
(2.81)

The phase equation (2.77) allows to split the phase vector φ into a longitudinal φ ℓ and a transverse φ t modes, φ = φ ℓ + φ t , that satisfy ∇ × φ ℓ = 0 and ∇ • φ t = 0 respectively. This leads to the uncoupled phase diffusion equations [25,14,26,27]

∂φ ℓ ∂t = D ∇ 2 φ ℓ , ∂φ t ∂t = D ⊥ ∇ 2 φ t .
(2.82)

The normal modes φ t and φ ℓ correspond to Eckhaus rectangular and rhomboidal perturbations, respectively. The hexagons are stable to phase modes in the domain defined by D > 0 and D ⊥ > 0. In the figure 2.11, we show the phase stability diagrams for a Newtonian fluid and two different values of the viscosities ratio r. Curves ( 1) and ( 2) correspond to D ⊥ = 0 and D = 0, respectively. Curve ( 4) is the upper stability amplitude where a bifurcation to rolls occurs. The minimum of the curve ( 4) is located in the region q < 0. Below curve (3), no hexagons exist. Hexagons are stable in the shaded region. For viscosities ratio 1 ≤ r ≤ 2, the region of stability to amplitude and phase modes is closed.

Whereas for larger values of r, the stability domain is open. Note also that the domain of stability is decentered towards the right. It is delimited mainly by the stability amplitude curves and the phase instability boundary. Nevertheless, the numerical results show that close to the threshold, the longitudinal mode is the relevant destabilizing mode. The region where the longitudinal mode destabilizes the pattern increases with increasing r. The influence of shear-thinning effects is clearly illustrated by figures 2.12 and 2.13. The stability domain of hexagons reduces with increasing shear-thinning effects. The region where the longitudinal mode is the relevant destabilizing mode remains small, even if the numerical results show that it increases with increasing shear-thinning effects. For comparison, we have represented in figure 2.13(b), the stability region of hexagons when α 1 = α 2 = 0.

This region is symmetrical with respect to the vertical axis.

Conclusion

We have investigated the influence of shear-thinning effects on Rayleigh-Bénard convection for a Carreau fluid, taking into account the variation of the viscosity with temperature.

The dependence of the viscosity on temperature was assumed of exponential type. The critical Rayleigh number Ra c for the onset of convection was determined as a function of the viscosity ratio, r, across the fluid layer. The results are in very good agreement with those given in the literature for a Newtonian fluid. The shear-thinning behavior does not play a role in the linear stability analysis of the conduction state.

In order to investigate nonlinear effects on bifurcation characteristics of convection onset,

(1) (2)

(3) (6) (7) (9)
(a) (b)

Figure 2.13 -Influence of shear-thinning effects and viscosities ratio on the hexagons stability diagram at P r = 100. Curves (1), ( 4) and ( 7) correspond to a Newtonian fluid and three different values of r : 2.5, 2.0 1.5 respectively. Curves (2), ( 5) and ( 6) correspond to a shear-thinning fluid with α = 5 × 10 -5 . Curves (3), ( 6) and ( 9) correspond to a shearthinning fluid with α = 10 -4 . (a) With quadratic spatial terms (b) Without quadratic spatial terms a weakly nonlinear analysis is performed using a multiple scale method. The procedure to derive amplitude equations is illustrated for a hexagonal lattice. Generalized Ginzburg-Landau equations are obtained including spatial nonvariational terms which describe distortion of hexagons. The stability analysis of the spatially homogeneous rolls and hexagons is first performed. It is shown that the threshold value of the relative Rayleigh-number ǫ h above which hexagons become unstable to rolls increases with the viscosity ratio and decreases with increasing shear-thinning effects. Similar evolution is observed for the threshold value ǫ r above which the rolls are stable. This behavior is qualitatively in agreement with [18]. The stability of hexagons with respect to long-wave perturbations is then addressed. Phase equation is derived and the stable band of wavenumbers is determined.

For moderate viscosity ratios, the domain of stable hexagons is open and decentered. It is clearly shown that shear-thinning effects reduce significantly the stable domain of hexagons. As for a Newtonian fluid, the stable hexagons domain is delimited mainly by the transverse phase instability. The longitudinal phase instability is relevant only for a very small range of ǫ.

CHAPITRE 2.

2.A Operators and matrix coefficients 2.A.1 The operator M

M =       M 11 0 0 0 M 21 M 22 0 0 M 31 0 M 33 0 0 0 0 1       (2.83)
with

M 11 = P r -1 ∆ , M 21 = ∂ 2 ∂x∂z , M 22 = ∇ 2 Hx , (2.84)
M 31 = ∂ 2 ∂y∂z , M 33 = ∇ 2 Hx .
The sub-scale M 0

The coefficients of M 0 in eq. (2.46) are

M (0) 11 = P r -1 ∆ Hx + ∂ 2 ∂z 2 , (2.85) M (0) 21 = ∂ 2 ∂x∂z , M (0) 22 = ∆ Hx , M (0) 31
= ∂ 2 ∂y∂z , M (0) 33 = ∆ Hx , M (0) 44 = 1 ,
The sub-scale M 1

The coefficients of M 1 in eq. (2.46) are

M (1) 11 = 2 P r -1 ∇ Hx • ∇ HX , (2.86) M (1) 21 = ∂ 2 ∂X∂z , M (1)
22 = 2 ∇ Hx • ∇ HX , M (1) 31
= ∂ 2 ∂Y ∂z , M (1)
33 = 2 ∇ Hx • ∇ HX , M (1
L =       L 11 0 0 L 14 L 21 L 22 0 0 L 31 0 L 33 0 1 0 0 ∆       (2.88) with L 11 = µ b ∆ 2 + 2 dµ b dz ∆ ∂ ∂z + d 2 µ b dz 2 ∂ 2 ∂z 2 -∆ H , L 14 = Ra ∆ H , (2.89)
L 21 = µ b ∆ ∂ 2 ∂x∂z + dµ b dz ∂ ∂z ∂ 2 ∂x∂z , L 22 = µ b ∆∆ H + dµ b dz ∂ ∂z ∆ H , L 31 = µ b ∆ ∂ 2 ∂y∂z + dµ b dz ∂ ∂z ∂ 2 ∂y∂z , L 33 = µ b ∆∆ H + dµ b dz ∂ ∂z ∆ H .

2.A.3 sub-scale L (0)

The components of L (0) in eq. (2.47) are

L (0) 11 = µ b ∇ 2 Hx + ∂ 2 ∂z 2 2 + 2 dµ b dz ∇ 2 Hx + ∂ 2 ∂z 2 ∂ ∂z + d 2 µ b dz 2 ∂ 2 ∂z 2 -∇ 2 Hx , L (0) 14 = Ra c ∇ 2 Hx , (2.90) L (0) 21 = µ b ∇ 2 Hx + ∂ 2 ∂z 2 ∂ 2 ∂x∂z + dµ b dz ∂ ∂z ∂ 2 ∂x∂z , L (0) 22 = µ b ∇ 2 Hx + ∂ 2 ∂z 2 ∇ 2 Hx + dµ b dz ∂ ∂z ∇ 2 Hx L (0) 31 = µ b ∇ 2 Hx + ∂ 2 ∂z 2 ∂ 2 ∂y∂z + dµ b dz ∂ ∂z ∂ 2 ∂y∂z , L (0) 33 = L (0) 22 , L (0) 41 = 1 , L (0) 44 = ∇ 2 Hx + ∂ 2 ∂z 2 .
CHAPITRE 2.

2.A.4 sub-scale L (1)

The components of L (1) in eq. (2.47) are

L (1) 11 = 4 µ b ∇ 2 Hx + ∂ 2 ∂z 2 (∇ Hx • ∇ HX ) + 4 dµ b dz (∇ Hx • ∇ HX ) ∂ ∂z -2 d 2 µ b dz 2 (∇ Hx • ∇ HX ) L (1) 14 = Ra (1) ∇ 2 Hx + 2Ra c ∇ Hx • ∇ HX , (2.91) L (1) 21 = µ b ∇ 2 Hx + ∂ 2 ∂z 2 ∂ 2 ∂X∂z + 2µ b (∇Hx • ∇ HX ) ∂ 2 ∂x∂z + dµ b dz ∂ 2 ∂z 2 ∂ ∂X L (1) 22 = 2µ b ∂ 2 ∂z 2 + 2 ∂ 2 ∂x 2 + ∂ 2 ∂y 2 (∇ Hx • ∇ HX ) L (1) 31 = µ b ∇ 2 Hx ∂ 2 ∂Y ∂z + 2 ∇ Hx • ∇ HX + ∂ 2 ∂z 2 ∂ 2 ∂Y ∂z + dµ b dz ∂ 2 ∂z 2 ∂ ∂Y L (1) 41 = 0 , L (1)
44 = 2∇ Hx • ∇ HX .

2.A.5 sub-scale L (2)

The components of L (2) in eq. ( 2.47) are (2.93)

L (2) 14 = Ra c ∇ 2 HX + 2Ra (1) ∇ Hx • ∇ HX + Ra (2) ∇ 2 Hx , (2.92)
L (2) 21 = µ b 2 (∇ Hx • ∇ HX ) ∂ 2 ∂X∂z + ∇ 2 HX ∂ 2 ∂x∂z L (2) 22 = µ b ∂ 2 ∂z 2 + ∇ 2 Hx ∇ 2 HX + 4µ b (∇ Hx • ∇ HX ) 2 + µ b ∇ 2 HX ∇ 2 Hx , L (2
Here w

(2)

1 means the first component of the second order solution. Similarly, for the horizontal velocity, we have

u (2) 1 = v (2) 1 = 0 .
(2.94)

2.B. SECOND-ORDER SOLUTION (HEXAGONS)

103

There is no velocity for the zero mode. The correction of the conductive temperature profile can be written as θ

(2)

1 = T 1 (z) |A 1 | 2 + |A 2 | 2 + |A 3 | 2 ,

2.B.2 Solution proportional to

A 2 p exp(2ik p • r)
The second component of the second order solution, proportional to A 2 p E 2 p , where E p = exp(ikp • r) represents the first harmonic of the fundamental. We have

w (2) 2 , θ (2) 2 = (W 2 (z), T 2 (z))) A 2 1 E 2 1 + A 2 2 E 2 2 + A 2 3 E 2 3 + c.c. , (2.97)
with

µ b D 2 -4k 2 c 2 + 2 dµ b dz D 3 -4k 2 c D + d 2 µ b dz 2 D 2 + 4k 2 c W 2 -4k 2 c Ra c T 2 = 2 P r F 11 D 3 F 11 -DF 11 D 2 F 11 -[N V w ] (2) A 2 p E 2 p , ( 2
.98)

W 2 + D 2 -4k 2 c T 2 = F 11 DG 11 -G 11 DF 11 .
(2.99)

The boundary conditions on W 2 and T 2 are identical to those on F 11 and G 11 , (2.29).

Concerning the horizontal velocity, we have

∇ 2 Hx u (2) H2 + ∇ Hx ∂w (2) 2 ∂z = 0 . (2.100)
We obtain, u

H2 = DW 2 4k 2 c ∇ Hx A 2 1 E 2 1 + A 2 2 E 2 2 + A 2 3 E 2 3 + c.c. . (2.101) (2)

2.B.3 Solution proportional to

A p A * q E p E * q
The third component of the second order solution, proportional to

A p A * q E p E * q , reads w (2)
3 , θ

(2) 3 = (W 3 (z), T 3 (z)) (A 1 A * 2 E 1 E * 2 + A 1 A * 3 E 1 E * 3 + A 2 A * 3 E 2 E * 3 ) + c.c. , (2.102) CHAPITRE 2.
with

µ b D 2 -3k 2 c 2 + 2 dµ b dz D 3 -3k 2 c D + d 2 µ b dz 2 D 2 + 3k 2 c W 3 -3k 2 c Ra c T 3 = 3 P r F 11 D 3 F 11 -k 2 c F 11 DF 11 -[N V w ] (2) EpE q *
(2.103)

W 3 + D 2 -3k 2 c T 3 = 2F 11 DG 11 -G 11 DF 11 .
(2.104)

Boundary conditions on (W 3 , T 3 ) are the same as the ones on (F 11 , G 11 ).

The horizontal velocity components satisfy

∇ 2 Hx u (2) H3 + ∇ Hx ∂w (2) 3 ∂z = 0 . (2.105)
We obtain u

(2)

H3 = DW 3 3k 2 c ∇ Hx (A 1 A * 2 E 1 E * 2 + A 1 A * 3 E 1 E * 3 + A 2 A * 3 E 2 E * 3 ) + c.c. (2.106)

2.B.4 Solution proportional to exp (ik p • r)

The fourth component of the second order solution is proportional to exp (ik p • r)

(resonant term). The solution is achieved using the solvability condition.

w (2) 4 , θ (2) 4
= (W 41 , T 41 ) For the horizontal velocity, we have

E 1 + (W 42 , T 42 ) E 2 + (W 43 , T 43 ) E 3 + c.c. , (2.107) with µ b ∂ 2 ∂z 2 -k 2 c 2 + 2 dµ b dz ∂ 2 ∂z 2 -k 2 c ∂ ∂z + d 2 µ b dz 2 ∂ 2 ∂z 2 + k 2 c W 41 -Ra c k 2 c T 41 = -4iµ b D 2 F 11 -k 2 c F 11 -4i dµ b dz DF 11 + 2i d 2 µ b dz 2 F 11 (k 1 • ∇ HX ) A 1 + k 2 c Ra 1 G 11 A 1 -2iRa c G 11 (k 1 • ∇ HX ) A 1 + (2.108) 1 P r F 11 D 3 F 11 + 2DF 11 D 2 F 11 -3k 2 c F 11 DF 11 A * 2 A * 3 -N V (2) w E 1 A * 2 A * 3 , W 41 + D 2 -k 2 c G 11 + T 41 = -2iG 11 (k 1 • ∇ HX ) A 1 + (2F 11 DG 11 + G 11 DF 11 ) A * 2 A * 3 . (2.
u (2) H4 = U H41 E 1 + U H42 E 2 + U H43 E 3 + c.c , (2.110) with µ b ∂ 2 ∂z 2 -k 2 c -k 2 c U H41 + i (k 1 • r) ∂W 41 ∂z + (2.111) dµ b dz -k 2 c ∂ ∂z U H41 + i (k 1 • r) ∂ 2 W 41 ∂z 2 = 2 dµ b dz D 2 F 11 ∇ HX A 1 + µ b D 3 F 11 -k 2 c DF 11 + dµ b dz D 2 F 11 2 k 1 • r k 2 c • (k 1 • ∇ HX A 1 ) -∇ HX A 1

Equations gouvernant le problème

On considère un fluide rhéofluidifiant non thermodépendant, placé entre deux plaques séparées d'une distance d, qui est chauffé par le bas à une température T1 et refroidi par le haut à une température T2 (Figure 3.1).

On note μ0 sa viscosité dynamique à cisaillement nul et κ0 sa diffusivité thermique. Sa masse volumique ρ varie linéairement avec la température dans le terme de poussée, conformément à l'hypothèse de Boussinesq :

ρ( T ) = ρ0 (1 -β( T -T0 )) avec T0 = T1 + T2 2 (3.1)
Le champ de température et la structure de l'écoulement sont décrits de manière classique par les équations de conservations de l'énergie et de la quantité de mouvement. Le champ de pression est éliminé en appliquant deux fois le rotationnel à l'équation de la quantité de

EQUATIONS GOUVERNANT LE PROBL ÈME

Méthode Pseudo-Spectrale

Décomposition de Fourier-Chebyshev

Les structures périodiques que l'on cherche à étudier sont des rouleaux de convection.

Ces motifs émergeant avec un nombre d'onde k c , il est alors pertinent de décomposer les champs de température et de vitesse en série de Fourier comme suit :

w(x, z) θ(x, z) = n=N n=-N w n (z) θ n (z) exp(ink c x) (3.13)
Comme w et θ sont réels, alors w -n = w * n et θ -n = θ * n (où * désigne le complexe conjugué). Chaque mode de Fourier est projeté sur la base des polynômes de Chebyshev de première espèce :

w n (z)

θ n (z) = M m=0 a nm b nm T m (Z) où Z = 2z (3.14)
avec T m (Z) le polynôme de Chebyshev de degré m.

En combinant les Eqs.(3.13) et (3.14), les approximations spectrales des champs de température et de vitesse sont données par :

w S (x, z) θ S (x, z) = N n=-N M m=0 a nm b nm T m (Z) exp(ink c x) (3.15) La condition w -n = w * n et θ -n = θ * n se traduit par : a -nm = a * nm et b -nm = b * nm (3.16)
Ainsi, on peut considérer uniquement les modes tels que n ≥ 0.

On peut également noter que l'approximation spectrale de u s'obtient directement à partir de celle pour w en considérant l'équation de continuité (3.2).

Ces différentes approximations spectrales sont évaluées aux points de collocation

M ÉTHODE PSEUDO-SPECTRALE

111 de Gauss-Lobatto définis par :

Z j = cos jπ M , j = 0...M (3.17)
Ce choix permet en effet d'éviter l'apparition du phénomène de Runge qui peut survenir avec des points équirépartis (cf [2] ou [3]).

Pour un mode E n = exp(ink c x), on note a n et b n les vecteurs des coefficients spectraux associés aux champs de vitesse et de température aux différents points de Gauss-Lobatto (3.17) :

a n b n = [a n0 , a n1 , ...a nM ] [b n0 , b n1 , ...b nM ] , ∀ 0 ≤ n ≤ N. (3.18)
Les valeurs prises par les approximations spectrales et leurs dérivées peuvent être reliées aux coefficients spectraux en introduisant les matrices de dérivation D0, D1, D2, D3 et D4 par rapport à z, définis par :

(Dα) i,j = 2 α T (α) i (Z j ) , α = 0, 1, 2, 3, 4 (3.19) où T (α) i
désigne la α ième dérivée du polynôme T i .

On adopte la notation [Φ S ] E n pour désigner la composante de Φ S suivant le mode E n pour tous les points de Gauss-Lobatto. Toutes les expressions du type [Φ S ] E n sont donc des vecteurs de taille M+1.

Les approximations spectrales et leurs dérivées pour le mode E n sont alors données par : -termes non linéaires visqueux NV, Ainsi :

∂ α ∂ β ∂z α ∂x β w S θ S E n = (ink c ) β Dα × a n Dα × b n . ( 3
F = F w F θ = L w + NI w + NV L θ + NI θ (3.23)
Dans les paragraphes suivants, on explicite la manière dont ces différents termes sont projetés dans l'espace spectral de Fourier-Chebyshev.

Projection des termes linéaires

Les termes linéaires L w et L θ sont évalués directement dans l'espace spectral de Fourier-Chebyshev. D'après l'Eq. (3.20), les expressions des termes linéaires projetées sur le mode E n = exp(ink c x) sont données par :

[L w ] E n = (D4 -2n 2 k 2 c D2 + n 4 k 4 c D0) × a n -n 2 k 2 c RaD0 × b n (3.24) [L θ ] E n = -D0 × a n -(D2 -n 2 k 2 c D0) × b n (3.25)

Projection des termes non linéaire d'inertie

Les termes non linéaires d'inertie présentent une non linéarité quadratique, il est possible d'obtenir directement leur décomposition dans l'espace spectral de Fourier-Chebyshev à partir des équations (3.5) et (3.7).

Ainsi la projection du terme NI θ sur un mode E n s'obtient par :

[NI θ ] E n = N l=n-N [u] E l ∂θ ∂x E n-l + [w] E l ∂θ ∂z E n-l (3.26)
En utilisant (3.20), on obtient :

[w] E l = D0 × a l ∂θ ∂z E n-l = D1 × b n-l [u] E l =    i l k c D1 × a l si l = 0 0 si l = 0 ∂θ ∂x E n-l = i(n -l)k c D0 × b n-l 3.3. M ÉTHODE PSEUDO-SPECTRALE 113
On peut procéder de la même manière pour le terme NI w . En explicitant tous les calculs, on obtient :

[NI w ] E n = N l=n-N 2 ∂u ∂x E l ∂ 2 u ∂x∂z E n-l + ∂u ∂z E l ∂ 2 u ∂x 2 E n-l + [u] E l ∂ 3 u ∂z∂x 2 E n-l + ∂ 2 w ∂x∂z E l ∂u ∂z E n-l + ∂w ∂x E l ∂ 2 u ∂z 2 E n-l + [w] E l ∂ 3 u ∂x∂z 2 E n-l + ∂w ∂z E l ∂ 2 u ∂x∂z E n-l - ∂ 2 u ∂x 2 E l ∂w ∂x E n-l -2 ∂u ∂x E l ∂ 2 w ∂x 2 E n-l -[u] E l ∂ 3 w ∂x 3 E n-l - ∂ 2 w ∂x 2 E l ∂w ∂z E n-l -2 ∂w ∂x E l ∂ 2 w ∂x∂z E n-l -[w] E l ∂ 3 w ∂x 2 ∂z E n-l
Chaque terme entre crochets est ensuite exprimé en fonction des coefficients spectraux et des matrices de dérivation de manière analogue au terme [NI θ ] E n . On obtient ainsi :

∂u ∂x E l = -D1 × a l si l = 0 0 sinon ∂ 2 u ∂x∂z E n-l = -D2 × a n-l si n = l 0 sinon ∂u ∂z E l =    i lk c D2 × a l si l = 0 0 sinon ∂ 2 u ∂x 2 E n-l = -i(n -l)k c D1 × a n-l ∂ 3 u ∂z∂x 2 E n-l = -i(n -l)k c D2 × a n-l ∂ 2 w ∂x∂z E l = ilk c D1 × a l ∂u ∂z E n-l =    i (n -l)k c D2 × a n-l si n = l 0 sinon ∂w ∂x E l = ilk c D0 × a l ∂ 3 u ∂x∂z 2 E n-l = -D3 × a n-l si n = l 0 sinon ∂w ∂z E l = D1 × a l ∂ 2 w ∂x 2 E n-l = -(n -l) 2 k 2 c D0 × a n-l ∂ 3 w ∂x 3 E n-l = -i(n -l) 3 k 3 c D0 × a n-l ∂ 3 w ∂x 2 ∂z E n-l = -(n -l) 2 k 2 c D1 × a n-l 114 CHAPITRE 3.

Projection des termes non linéaire visqueux

A cause de la non linéarité complexe (non polynomiale) du modèle de

: Ω = (x, z) ∈ [0, L c ] × -1 2 , 1 2 .
Les termes non linéaires visqueux NV sont évalués à l'aide des approximations spectrales de u S et w S (3.15) et des matrices de dérivation Dα (3.19) aux points de grille : projection dans l'espace spectral de Fourier-Chebyshev (cf [4] et [5] ). Cette règle stipule qu'on doit avoir :

(x i , z j ) = iL c N d , Z j 2 où Z j = cos jπ M d -1 (3.
N d ≥ 3 2 (2N + 1) et M d ≥ 2M (3.
F (x, z) = N n=0 M m=0 f nm T m (Z)E n (3.29)
Les fonctions E n où n = 0...N sont orthogonales vis à vis du produit scalaire défini par :

<f (x), g(x)> E = 2π kc 0 f (x)g(x) * dx (3.30)
En effet, on a :

<E n , E p > E =    0 si n = p 2π k c sinon (3.31)
De la même manière, les polynômes de Chebyshev sont orthogonaux vis à vis de l'application :

<F (Z), G(Z)> T = 1 -1 F (Z)G(Z) √ 1 -Z 2 dZ (3.32)
On montre que (cf [3]) :

<T m , T q > T = 0 si m = q c q π 2 si m = q où c q = 2 si q = 0 1 sinon (3.33)

On réalise le produit scalaire, défini par (3.30), de F (x, z) et d'un mode E p soit :

F p = <F (x, z), E p > E = N n=0 M m=0 f nm T m (z)<E n , E p > E (3.34)
D'après (3.31), on obtient :

F p = <F (x, z), E p > E = M m=0 f pm T m (z) × 2π k c (3.35)
En appliquant maintenant < . , . > T défini par (3.32) à l'expression précédente et à un polynôme T q (Z) quelconque, on obtient :

<F p , T q (Z)> T = 2π k c M m=0 f pm <T m (Z), T q (Z)> T (3.36)
soit avec (3.31) :

<F p , T q (Z)> T = 2π k c π 2 c q f pq (3.37) CHAPITRE 3.
Ainsi, le coefficient f pq s'obtient par :

f pq = k c 2π 2 π 1 c q <F p , T q (Z)> T , ∀0 ≤ p ≤ N et 0 ≤ q ≤ M (3.38)
Or par définition,

F p = <F (x, z), E p > E = 2π kc 0 F (x, z) exp(-ipk c x)dx (3.39)
donc finalement en substituant l'expression (3.39) dans (3.38), on obtient pour f qp , (∀ 0 ≤ p ≤ N et 0 ≤ q ≤ M ) :

f pq = k c 2π 2 π 1 c q 2π kc 0 1 -1 F (x, z)T q (Z) exp(-ipk c x) √ 1 -Z 2 dxdZ (3.40) avec Z = 2z et c q = 2 si q = 0 1 sinon
Dans la pratique, on ne dispose que de l'évaluation de F (x, z) aux points de la grille définis par (3.27). Pour évaluer (3.40), on est donc amené à utiliser une formule de quadrature. Dans notre cas, où on interpole les grandeurs par des polynômes de Chebyshev évalués aux points de Gauss-Lobatto, la version discrète de (3.40) est donnée par la formule suivante (cf [6] et [5]) :

f pq = 1 γ q 1 N d N d -1 i=0 M d -1 j=0 F (x i , z j )T q (Z j ) exp(-ipk c x i )W j (3.41)
où W j sont les poids de la formule de quadrature définis par :

     W 0 = W M d -1 = π 2(M d -1) W j = π M d -1 , 1 ≤ j ≤ M d -2 (3.42)
Le coefficient γ q est quant à lui défini par :

     γ q = π 2 c q si q < N où c q = 2 si q = 0 1 sinon γ N = π (3.
u ∂θ ∂x + w ∂θ ∂z = w + ∂ 2 θ ∂x 2 + ∂ 2 θ ∂z 2 ⇔ u ∂θ ∂x - ∂ 2 θ ∂x 2 + w ∂θ ∂z - ∂ 2 θ ∂z 2 = w ⇔ ∂(uθ) ∂x -θ ∂u ∂x + ∂(wθ) ∂z -θ ∂w ∂z - ∂ 2 θ ∂x 2 - ∂ 2 θ ∂z 2 = w ⇔ ∂(uθ) ∂x - ∂ 2 θ ∂x 2 + ∂(wθ) ∂z - ∂ 2 θ ∂z 2 -θ ∂u ∂x + ∂w ∂z =0 = w ⇔ ∂ ∂x uθ - ∂θ ∂x + ∂ ∂z wθ - ∂θ ∂z = w 118 CHAPITRE 3.
On peut intégrer l'équation précédente sur une période :

Lc ∂ ∂x uθ - ∂θ ∂x dx + Lc ∂ ∂z wθ - ∂θ ∂z dx = Lc wdx
Compte tenu de la décomposition des variables en série de Fourier (3.13), la première intégrale du membre de gauche s'annule par périodicité. La vitesse verticale moyenne étant nulle dans le fluide (w 0 = 0), le membre de droite s'annule également. On obtient :

Lc ∂ ∂z wθ - ∂θ ∂z dx = 0, ∀ - 1 2 < z < 1 2
En intervertissant les opérateurs de dérivation et d'intégration, on obtient finalement :

Lc wθ - ∂θ ∂z dx = cste, ∀ - 1 2 < z < 1 2
En évaluant cette expression en z = -1 2 où w = θ = 0, on a :

Lc wθ - ∂θ ∂z dx = - Lc ∂θ ∂z z=0 dx, ∀ - 1 2 < z < 1 2
On en déduit donc que :

Nu = 1 + 1 L c Lc wθ - ∂θ ∂z dx, ∀ - 1 2 < z < 1 2
En tenant compte du fait que w n = w -n , θ n = θ -n dans (3.13) et de la périodicité de l'écoulement, on a finalement :

Nu = 1 - ∂θ 0 ∂z + w 0 θ 0 + 2 N n=1 w n θ n , ∀ - 1 2 < z < 1 2
On est ainsi en mesure de calculer le nombre de Nusselt en utilisant les données en différents points de l'écoulement. Cela permet d'avoir une garantie supplémentaire sur la validité des calculs numériques.

Procédure de continuation

On a vu que la méthode pseudo-spectrale décrite précédemment conduisait à l'obtention d'un système d'équations non linéaires de la forme :

F(X, Ra) = 0 avec F : R NT × R → R NT et NT = 2(N + 1)(M + 1) (3.44)
où les inconnues du problème X k , k = 1...NT étaient les coefficients spectraux associés aux champs de vitesse et de température.

Étudier les caractéristiques de l'écoulement quand on s'éloigne du seuil de la convection revient donc à résoudre le système (3.44) pour différentes valeurs du nombre de Rayleigh.

Généralement, un problème tel que (3.44) admet des branches continues de solutions, fonctions du paramètre de contrôle à savoir ici le nombre de Rayleigh. Les procédures de continuation désignent l'ensemble des techniques permettant de suivre ces branches. Ce sont des processus itératifs permettant de déterminer de proche en proche les solutions du problème au fur et à mesure que le paramètre de contrôle varie : on obtient des couples (X i , Ra i ) solutions de (3.44).

De manière générale, de telles procédures se font en 4 étapes (Figure 3.2) :

1. Une étape de prédiction (X i+1 , Ra i+1 ) de la solution recherchée à partir de la solution précédente (X i , Ra i ) sur la branche.

2. Une étape de paramétrisation de la branche de solution X(Ra) telle que F(X(Ra), Ra) = 0.

3. Une étape de correction permettant de déterminer la solution exacte (X i+1 , Ra i+1 ) à partir de la prédiction (X i+1 , Ra i+1 ).

4. Une étape de contrôle de la taille du pas utilisé lors de la phase de prédiction suivante.

La mise en place des différentes étapes est détaillée dans les paragraphes suivants. On Il définit alors les opérateurs matriciels

(X i ,Ra i ) (X i 1 ,Ra i 1 ) (X i
∂ X F et ∂ λ F respectivement de dimensions (NT, CHAPITRE 3.
NT) et (NT, 1) tels que :

[∂ X F] pq = ∂F p ∂X q [∂ λ F] p = ∂F p ∂λ ∀1 ≤ p, q ≤ NT (3.45)
L'opérateur ∂ X F est la matrice jacobienne associé à F(X, λ).

(X 0 , λ 0 ) est un point de bifurcation à condition que :

1. F(X 0 , λ 0 ) = 0 ; 2. rang (∂ X F(X 0 , λ 0 )) = NT -1 ;
3. ∂ λ F(X 0 , λ 0 ) appartienne au sous espace vectoriel engendré par les colonnes de

∂ X F(X 0 , λ 0 ) ;
4. deux branches de solutions se croisent avec des tangentes distinctes.

Un point de retournement (X 0 , λ 0 ) vérifie, quant à lui, les propriétés suivantes :

1. F(X 0 , λ 0 ) = 0 ; 2. rang (∂ X F(X 0 , λ 0 )) = NT -1 ;
3. ∂ λ F(X 0 , λ 0 ) n'appartient pas au sous espace vectoriel engendré par les colonnes de

∂ X F(X 0 , λ 0 ) ; 4. il existe une paramétrisation X(σ) et λ(σ) telle que X(σ 0 ) = X 0 , λ(σ 0 ) = λ 0 et d 2 λ(σ 0 )
dσ 2 = 0. On peut retenir de ces deux définitions que la matrice jacobienne ∂ X F n'est inversible ni aux points de retournement ni aux points de bifurcation puisque rang (∂ X F(X 0 , λ

0 )) = NT -1. B R 1 R 2 X λ

Prédiction

La phase de prédiction consiste à définir une approximation (X i+1 , Ra i+1 ) de la solution à partir du point précédent (X i , Ra i ) sur la branche. Le premier point (X 0 , Ra 0 ) provient d'une analyse faiblement non linéaire de stabilité semblable à celle utilisée dans les chapitres précédents.

Il existe différents prédicteurs. Le plus simple d'entre eux, dénommé prédicteur trivial, consiste à choisir comme approximation en Ra :i+1 la solution obtenue précédemment sur la branche soit :

(X i+1 , Ra i+1 ) = (X i , Ra i+1 ) (3.46)
Dans la pratique, ce prédicteur conduit à une approximation très grossière de la solution ce qui s'avère très pénalisant pendant la phase de correction où de nombreuses itérations sont alors nécessaires pour obtenir la solution exacte.

Pour obtenir une meilleure approximation, on peut commencer par calculer la différentielle de (3.44) ce qui conduit à la relation :

dF = 0 = (∂ X F)dX + (∂ Ra F)dRa (3.47)
où X est le vecteur des coefficients spectraux et où ∂ X F et ∂ Ra F sont les opérateurs matriciels définis par (3.45).

A partir de (3.47), on a :

dX dRa = -(∂ X F) -1 (∂ Ra F) (3.48)
On peut alors obtenir une approximation de la solution en Ra i+1 en intégrant cette dernière relation, soit :

X i+1 = X i + (Ra i+1 -Ra i ) dX dRa i (3.49)
Cette estimation de la solution se révèle bien meilleure que celle obtenue dans le cas du prédicteur trivial : elle permet des incréments plus grands du paramètre de contrôle Ra au cours de la procédure de continuation. Toutefois elle présente l'inconvénient important d'échouer aux points de retournement où (∂ X F) -1 n'est pas définie.

Le prédicteur de la tangente présenté dans [8] permet de s'affranchir de ce problème. On introduit un vecteur tangent T défini par :

T k = dX k pour 1 ≤ k ≤ N T T NT+1 = dRa
On peut alors réécrire (3.47) sous la forme : Une fois le vecteur tangent T calculé, l'approximation (X i+1 , Ra i+1 ) de la solution est donnée par :

(∂ X F|∂ Ra F)T = 0 (3.50) CHAPITRE 3.
(X i+1 , Ra i+1 ) = (X i , Ra i ) + σ i T (3.53)
La grandeur σ i qui intervient dans la formule précédente (3.53) est le pas de la méthode de prédiction. Le choix d'une valeur adéquate de ce paramètre s'avère délicat. Considérer un faible pas va aboutir à avoir une bonne estimation de la solution et peu d'itérations lors de la phase de correction ; en revanche, on risque de peu se déplacer le long de la branche de solution. A l'inverse si on choisit un pas trop grand, l'estimation de la solution va être mauvaise et la phase de correction risque d'échouer ou d'être excessivement longue.

Différentes stratégies mises en oeuvre dans le choix de ce paramètre sont présentées ultérieurement.

Dans le code numérique mis en place, c'est le prédicteur de la tangente qui est utilisé. Il est à noter que d'autres prédicteurs basés sur une approximation polynomiale locale de la courbe de solutions peuvent également être employés. Le plus commun d'entre eux, appelé prédicteur de la sécante, consiste à choisir comme approximation :

(X i+1 , Ra i+1 ) = (X i , Ra i ) + σ i (X i -X i-1 , Ra i -Ra i-1 ) (3.54)
D'après [8], l'erreur sur la solution fournie par ce prédicteur est du même ordre de grandeur que celle obtenue dans le cas de la tangente. Elle requiert cependant les données en deux points précédents : c'est pour cette raison que la tangente lui a été préférée. Une branche de solutions peut être décrite mathématiquement par l'intermédiaire d'une grandeur appelée paramètre : on parle alors de paramétrisation. Une fois la paramétrisation effectuée, chaque point de la courbe de solutions est parfaitement identifié par une valeur du paramètre.

Pour une courbe donnée, différentes paramétrisations peuvent être possibles. Le paramètre le plus naturel pour décrire une branche de solutions est le paramètre de contrôle du problème, dans notre cas le nombre de Rayleigh Ra. Ce choix permet de conserver la dimension physique du problème mais s'avère inapproprié aux points de retournement :

en conséquence, on choisit généralement d'autres paramétrisations.

De manière générale, utiliser une grandeur γ comme paramètre revient à considérer que les solutions de F(X, Ra) = 0 sont des fonctions continues de γ soit : Plusieurs équations de paramétrisation sont alors possibles. Une des plus fréquentes repose sur une approximation par différences finies des dérivées de (3.58) soit :

X = X(γ) Ra = Ra(γ) ( 3
p 1 (X, Ra, s) = NT k=1 (X k -X k (s i )) dX ds s i + (Ra -Ra(s i )) dRa ds s i -(s -s i ) 2 = 0 (3.59)
Cette paramétrisation fait intervenir les dérivées de X et de Ra par rapport à l'abscisse curviligne s évaluées au point précédent. Une façon pratique de les évaluer est présentée dans [4]. Toutefois, il semblerait qu'il existe des configurations dans lesquelles cette évaluation se révèle imprécise et où la procédure de continuation échoue (cf [10]).

Une autre équation de paramétrisation fréquemment utilisée est donnée par : Enfin, on peut signaler qu'il existe des alternatives aux paramétrisations par longueur d'arc et pseudo longueur d'arc. On peut notamment mentionner les paramétrisations locales qui permettent également de franchir les points de retournement. Dans ces approches, le paramètre utilisé, une composante de X ou Ra, est déterminé localement à chaque pas de la méthode de continuation. On pourra se tourner vers [8] pour de plus amples informations.

p 2 (X, Ra, s) = NT k=1 (X k -X k (s i )) 2 + (Ra -Ra(s i )) 2 -(s -s i ) 2 = 0 (3.

Correction

La phase de correction a pour but de déterminer une solution exacte (X i+1 , Ra i+1 ) du système étendu F (X, Ra, s) = 0 à partir de l'estimation (X i+1 , Ra i+1 ) obtenue lors de la phase de prédiction.

La correction à faire à partir de l'estimation est déterminée en appliquant la méthode de Newton au système étendu. Cette dernière appartient à la classe des méthodes de points fixes et permet de déterminer les zéros d'une fonction de manière itérative. Ainsi, à partir d'un point initial (X 0 i , Ra 0 i ) = (X i+1 , Ra i+1 ), suffisamment proche de la solution du problème, on construit une suite de points (X 1 i , Ra 1 i ), (X 2 i , Ra 2 i ), ..., (X j i , Ra j i ) ..., qui

PROC ÉDURE DE CONTINUATION

125 converge vers une solution de F (X, Ra, s) = 0 soit :

lim j→∞ X j i Ra j i = X i+1 Ra i+1 où F (X i+1 , Ra i+1 , s i+1 ) = 0 (3.61)
Les itérations de Newton successives s'obtiennent par la relation de récurrence suivante :

X j+1 i Ra j+1 i = X j i Ra j i -J -1 (X j i ,Ra j i ) F p 2 (X j i ,Ra j i ) (3.62) où J = ∂ X F ∂ Ra F ∂ X p 2 ∂ Ra p 2
est la matrice jacobienne du système étendu (3.57). Son inverse J -1 est bien défini aux points de retournement contrairement à la jacobienne (∂ X F) -1 du système de base (3.44).

L'opérateur matriciel ∂ X F de dimensions (NT,NT) intervenant dans la jacobienne J est calculé par différences finies centrées :

[∂ X F] pq = ∂F p ∂X q = lim δ→0 F p (X + δe q ) -F p (X -δe q ) δ (3.63) où [e q ] k = δ qk pour k = 1...NT.
Les autres opérateurs ∂ Ra F, ∂ X p 2 respectivement de dimensions (NT, 1) et (1, NT) ainsi que le scalaire ∂ Ra p 2 sont quant à eux calculés analytiquement.

La convergence de la méthode de Newton est quadratique (cf [11]) :

||X i+1 -X j+1 i || 2 ≤ C||X i+1 -X j i || 2 2 (3.64) où C est une constante et où ||X|| 2 = NT k=1 X 2 k 1 2
.

Dans la pratique, on a considéré que la fonction F s'annulait lorsque sa norme devenait inférieure à une valeur ζ choisie très faible. Les itérations de Newton sont ainsi menées jusqu'à ce que :

||F (X j+1 i , Ra j+1 i , s i+1 )|| 2 ≤ ζ (3.65)
Dans les simulations numériques, on a fixé : ζ = 10 -6 .

Plusieurs difficultés peuvent survenir dans l'application de la méthode de Newton.

La première, inhérente à la méthode, réside dans le fait que la convergence est locale : si le point initial (X 0 i , Ra 0 i ) = (X i+1 , Ra i+1 ) est trop éloigné de la solution du problème, les itérés risquent de diverger. La seconde difficulté provient de l'impossibilité de calculer de manière analytique le terme ∂ X F : on est obligé de procéder par différences finies.

Numériquement, définir une valeur pour la grandeur δ considérée comme infiniment petite dans (3.63) peut se révéler ardu : si la valeur est trop faible, les variations de la fonction CHAPITRE 3.

F risquent d'être trop faibles pour être perceptibles alors qu'à l'inverse une valeur trop forte mène à une estimation grossière de ∂ X F. Par tâtonnement, la valeur retenue a été fixée à δ = 10 -7 . Par ailleurs, déterminer ∂ X F numériquement nécessite des temps de calcul importants : la raison en incombe aux multiples projections dans l'espace spectral nécessaires pour évaluer F aux différents points.

Contrôle du pas

Le comportement des branches de solution peut varier fortement en fonction des valeurs du paramètre de contrôle. Dans certaines gammes de valeurs, les branches vont avoir une allure monotone alors que dans d'autres, leurs évolutions seront erratiques. Dans le premier cas, utiliser de grands incréments du paramètre de contrôle permet à la fois de bien capter la forme de la branche et de réduire les temps de calcul. En revanche dans le second cas, de plus petits incréments sont nécessaires afin de capter l'allure des branches de solutions.

Il paraît alors nécessaire d'adapter les incréments du paramètre de contrôle en fonction du comportement des branches de solutions au cours de la procédure de continuation : c'est l'objectif de l'étape de contrôle du pas.

On utilise généralement des arguments empiriques pour esquisser l'allure des branches. On peut par exemple considérer le nombre d'itérations effectuées lors de la phase de correction comme un bon indicateur de leurs variations : un faible nombre d'itérations indique une évolution monotone tandis qu'un fort nombre d'itérations correspond à une évolution brusque. Dans la pratique, on a estimé que la situation optimale permettant de bien saisir l'allure de la branche tout en minimisant les temps de calcul correspondait à faire une dizaine d'itérations par étape de correction. Ainsi, quand le nombre d'itérations effectuées N it est trop en dessous de cette valeur N opt = 10, on choisit d'augmenter les incréments d'abscisses curvilignes ∆s = s i+1s i et le pas σ i utilisé lors de l'étape de prédiction. De manière similaire, on diminue ces grandeurs quand on se trouve trop au dessus de N opt .

Plus précisément, dans le code mis en place, on a choisi :

-si N it ≤ 5, ∆s i+1 = 1.1∆s i et σ i+1 = 1.1σ i -si N it ≥ 20, ∆s i+1 = 1 2 ∆s i et σ i+1 = 1 2 σ i

Résultats pour des fluides Newtoniens

Dans cette section, on présente les résultats obtenus en appliquant la procédure de continuation au cas Newtonien. On donne les évolutions des différents modes de Fourier, des champs de vitesse et de température et des transferts thermiques au fur et à mesure qu'on s'éloigne des conditions critiques. Les résultats sont comparés avec ceux issus de différentes études qui sont présentées en premier lieu.

Convergence de la méthode de continuation en Newtonien

La procédure de continuation a été menée en utilisant différents nombres N de modes de Fourier afin de s'assurer de sa convergence dans le cas Newtonien. Les résultats numériques montrent que le nombre de modes nécessaire à la convergence ne dépend pas du nombre de Prandtl lorsque Pr ≥ 1. On présente donc uniquement le cas Pr = 10 dans ce paragraphe.

Le nombre de Nusselt Nu a été utilisé comme grandeur physique représentative de la convergence de la méthode. La

\ | セ @ j 1 Z Z Z セ @ V'ISTAE!t tNSnsLE ., \ \ 20C() j" 1 IOC•O 0 • ' • $ • T • ,_ ǫ 0 2 4
E(w n ) = 1 2 -1 2 w 2 n (z)dz E(θ n ) = 1 2 -1 2 θ 2 n (z)dz (3.66)
On définit aussi les valeurs relatives E r (w n ) et E r (θ n ) comme suit : des résultats de W. Plows. Enfin, on peut signaler que nos résultats sont en bon accord avec les résultats expérimentaux obtenus par Rossby [14].

E r (w n ) = E(w n ) N n=0 E(w n ) E r (θ n ) = E(θ n ) N n=0 E(θ n ) (3.
Des comparaisons ont également été réalisées avec les travaux de W. Plows [7] Récemment, des travaux numériques ont été menés afin de suivre la branche stable de solution sous-critique mise en évidence dans [20] et [1] lorsque les effets rhéofluidifiants sont suffisamment importants. Une première étude a été réalisée par O. Benouared, M. Mahmoud et N.A. Messaoudene en 2014 [21]. Ces auteurs ont considéré un fluide rhéofluidifiant confiné dans une cavité rectangulaire. La rhéologie a été décrite en utilisant le modèle de Carreau-Yasuda (cf [17]) : Les travaux de Benouared et al. [21] et de Jenny et al. [22] se sont concentrés sur la bifurcation sous-critique qui survenait pour de forts effets rhéofluidifiants : en conséquence, la gamme de nombres de Rayleigh explorée est relativement restreinte (0 ≤ Ra ≤ 10000 dans [21] et 0 ≤ Ra ≤ 5000 dans [22]). Notre démarche vise davantage à étudier l'influence d'un effet rhéofluidifiant sur le devenir des rouleaux de convection quand on s'éloigne des conditions critiques : on cherche donc à mener la procédure de continuation dans une plage de nombres de Rayleigh plus étendue que dans [21] et [22].

α = 1 -n c 2 λ 2 . ( 3
µ -µ ∞ 1 -µ ∞ = (1 + (λ √ Γ) a ) nc-1 a , ( 3

Convergence de la méthode de continuation pour un fluide rhéofluidifiant

La principale difficulté dans l'application de la procédure de continuation à un fluide rhéofluidifiant provient de la projection des termes non linéaires visqueux dans l'espace de Fourier Chebyshev suivant la procédure décrite dans 3.

Conclusion et perspectives

Dans ce chapitre, une procédure de continuation a été mise en place afin de déterminer des branches de solutions stationnaires en rouleaux apparaissant dans la convection de Rayleigh-Bénard. Afin de tester sa validité, la méthode a d'abord été utilisée en considérant un fluide Newtonien. Les résultats obtenus se sont avérés en excellent accord avec plusieurs articles [12], [14], [7], [13].

Le code a ensuite été utilisé dans le cas d'un fluide rhéofluidifiant. Là encore, on est parvenu à retrouver les principaux résultats de la littérature [21], [22] :

-en diminuant le paramètre n c ou en augmentant le paramètre λ du modèle de Dans la deuxième configuration, on considère un fluide rhéofluidifiant dont la viscosité varie avec la température. On a étudié la compétition rouleaux/carrés et rouleaux/hexagones pour différentes valeurs du degré de rhéofluidification, α, et du rapport de viscosité, r. Dans le cas Newtonien, le motif stable au démarrage de la convection passe de rouleaux à carrés à r = r c ≈ 3 en accord avec [1]. Pour un fluide rhéofluidifiant, 154 le motif stable de convection passe de rouleaux à carrés pour des valeurs de r qui augmentent avec le degré de rhéofluidification. Cette variation est en accord avec les conclusions de [2]. Pour des valeurs modérées de r (disons r < 3), l'étude de la compétition rouleaux/hexagones montre que l'étendue du domaine de bistabilité rouleaux/hexagones diminue de manière significative avec le degré de rhéofluidification at augmente avec r.

En s'écartant des conditions critiques, une gamme continue de nombres d'onde d'étendue √ ǫ participe à la convection conduisant à une modulation spatiale de la solution stationnaire. La stabilité des hexagones présentant une modulation spatiale est étudiée dans le cadre des équations de Ginsburg-Landau qui incluent les dérivées spatiales quadratiques de l'amplitude. Après avoir évalué les différents coefficients qui interviennent dans ces équations, une analyse de stabilité par rapport à des perturbations de grande longueur d'onde est effectuée. On montre que pour des valeurs modérées du rapport de viscosité, le domaine du plan (ǫ, q) où q = kk c où les hexagones sont stables est ouvert et décentré vers les valeurs positives de q. Il est délimité essentiellement par une instabilité de phase transverse. Le caractère rhéofludifiant, encore une fois, réduit ce domaine de stabilité. La résolution numérique des équations de Ginzburg-Landau montre en particulier la formation de défauts probablement du type heptagone-pentagone liés à l'instabilité de phase transverse [3].

Après avoir considéré les motifs de convection à proximité du seuil, nous nous sommes inétéressés à leur devenir loin des conditions critiques. Pour cela, une approche fortement non linéaire reposant sur une procédure de continuation a été mise en place afin de déterminer de proche en proche les caractéristiques de l'écoulement. On a considéré lors de cette approche un fluide purement rhéofluidifiant confiné entre deux plaques parfaitement conductrices : dans cette configuration, seuls les rouleaux sont stables. Le comportement rhéologique est décrit par le modèle de Carreau. La méthode mise en oeuvre a d'abord été validée par comparaison avec différents travaux réalisés en Newtonien tels que [4], [5], [6] et [7]. On a pu mettre en évidence l'influence de la rhéologie : (i) quand on augmente l'effet rhéofluidifiant (que ce soit en diminuant le paramètre n c ou en augmentant le paramètre λ du modèle de Carreau), on réduit les frottements visqueux ce qui a pour conséquence de faciliter la convection thermique et les échanges de chaleur dus à cette dernière, (ii) à partir d'une certaine valeur du degré de rhéofluidification alpha, la bifurcation de l'état conductif à l'état convectif devient sous-critique. Ces résultats sont en adéquation avec ceux obtenus dans la littérature [8], [9].

Plusieurs perspectives s'ouvrent suite aux travaux réalisés dans cette thèse : -la masse volumique ne varie avec la température que dans le terme de flottabilité de l'équation de la conservation de quantité de mouvement.

-la dissipation visqueuse ainsi que la dépendance de l'énergie interne avec la pression sont négligées dans l'équation de l'énergie.

-les coefficients thermodynamiques sont considérés constants et la masse volumique ne varie qu'avec la température (pas de compressibilité). Pour étudier les motifs de convection loin des conditions critiques, on a recours à une procédure de continuation permettant de déterminer de proche en proche les caractéristiques de l'écoulement comme les champs de vitesse ou de température ainsi que le nombre de Nusselt.

Mots clés : Convection, Rayleigh-Bénard, Rhéologie, Bifurcation, Instabilités secondaires.

  Introduction générale Bibliographie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Instabilities of convection patterns for a shear-thinning fluid between slabs of finite conductivities 1.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Linear stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Critical conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Characteristic time of the instability . . . . . . . . . . . . . . . . . . 1.5 Pattern selection at the onset of convection . . . . . . . . . . . . . . . . . . 1.5.1 Principles of the amplitude expansion method : case of square . . . . 1.5.2 Nature of the primary bifurcation . . . . . . . . . . . . . . . . . . . 1.5.3 Convective patterns at threshold . . . . . . . . . . . . . . . . . . . . 1.5.4 Flow structure, viscosity and temperature fields . . . . . . . . . . . . 1.5.5 Comparison between roll and square solutions for fixed χ and λ . . . 1.5.6 Heat transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6 Secondary instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.1 Stability of a roll modulated solution . . . . . . . . . . . . . . . . . . 1.6.2 Stability of a square modulated solution . . . . . . . . . . . . . . . . 1.7 Numerical solutions of amplitude equations . . . . . . . . . . . . . . . . . . 1.7.1 Numerical simulation . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.A Eigenvalues versus Q 1 and Q 2 . . . . . . . . . . . . . . . . . . . . . . . . . . 1.B Stability of the ETD2 scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 1.B.1 Case of a single ODE . . . . . . . . . . . . . . . . . . . . . . . . . . 1.B.2 Extension to two coupled Ginzburg-Landau equations . . . . . . . . Bibliographie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . TABLE DES MATI ÈRES 3.5.1 Prédiction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Paramétrisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.3 Correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.4 Contrôle du pas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Résultats pour des fluides Newtoniens . . . . . . . . . . . . . . . . . . . . . 3.6.1 Les principaux travaux existants . . . . . . . . . . . . . . . . . . . . 3.6.2 Convergence de la méthode de continuation en Newtonien . . . . . . 3.6.3 Évolution des différents modes de Fourier avec l'écart au seuil de la convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.4 Évolution des champs de vitesse et de température avec l'écart au seuil de la convection . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.5 Influence du nombre de Prandtl . . . . . . . . . . . . . . . . . . . . . 3.7 Résultats pour des fluides rhéofluidifiants . . . . . . . . . . . . . . . . . . . 3.7.1 Les principaux travaux existants en non-Newtoniens . . . . . . . . . 3.7.2 Convergence de la méthode de continuation pour un fluide rhéofluidifiant . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.3 Évolution des différents modes de Fourier avec l'écart au seuil de la convection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.4 Évolution des champs de vitesse, de température et de viscosité avec l'écart au seuil de la convection . . . . . . . . . . . . . . . . . . . . . 3.7.5 Influence de l'indice de rhéofluidification n c sur les échanges de chaleur145 3.7.6 Influence du paramètre λ sur les échanges de chaleur et la nature de la bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Conclusion et perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliographie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion Bibliographie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  La compétition entre ces deux mécanismes antagonistes est quantifiée par une grandeur adimensionnelle appelée nombre de Rayleigh et défini par : expression précédente, ρ0 , μ0 et κ0 désignent respectivement la masse volumique, la viscosité dynamique et la diffusivité thermique du fluide. L'indice 0 signifie que ces grandeurs sont évaluées à une température de référence, en général celle qui règne au centre de la couche de fluide. La grandeur β désigne le coefficient de dilatation thermique du fluide et est à l'origine de la poussée d'Archimède. d est l'épaisseur de la couche de fluide. ∆ T = T1 -T2 est l'écart de température entre les plaques avec T1 et T2 qui désignent respectivement les températures chaude et froide. A partir d'une différence de températures assez grande donc pour une valeur suffisamment élevée du nombre de Rayleigh notée Ra c , la poussée d'Archimède l'emporte sur les forces de friction et le fluide initialement au repos se met en mouvement : on passe d'un régime purement conductif à un régime convectif (Figure 3.1). C'est la convection de Rayleigh-Bénard.

Figure 2 -

 2 Figure 2 -Illustration de la convection de Rayleigh-Bénard.

Figure 3 -

 3 Figure 3 -Courbes de stabilité marginale obtenues dans le cas Libre-Libre et Rigide-Rigide

Figure 4 -

 4 Figure 4 -Taux de croissance des perturbations à proximité du seuil.

Figure 6 -Figure 7 -

 67 Figure6-Illustration de l'instabilité de compression/ dilatation d'Eckhaus tirée de[18]

4 .

 4 Les différents travaux précédents, aussi bien ceux s'intéressant aux conditions critiques que ceux cherchant à déterminer la nature du motifde convection et sa stabilité, font intervenir des fluides Newtoniens. Toutefois, dans de nombreuses situations, les fluides rencontrés ont un comportement rhéologique plus complexe. C'est le cas dans de nombreux procédés industriels ou en géophysique où la majorité des fluides rencontrés sont non-Newtoniens. Le caractère rhéofluidifiant, décroissance non linéaire de la viscosité avec le cisaillement est le caractère le plus commun. Dans ce qui suit, on présente une revue bibliographique des travaux théoriques numériques et expérimentaux consacrés à l'influence du comportement rhéofluidifiant sur la convection de Rayleigh-Bénard et le démarrage de celle-ci. Dans les études théoriques et numériques, le caractère rhéofluidifiant peut être isolé. Par contre dans les études expérimentales, les fluides peuvent être à un certain degré viscoélastiques ou thixotropes. Les résultats expérimentaux doivent être interprétés avec précaution. A notre connaissance, la première étude expérimentale de convection de Rayleigh-Bénard dans une couche de fluide rhéofluidifiant confinée entre deux parois a été effectuée par [28]. Les fluides utilisés sont des solutions aqueuses de Methocel à 1% en masse, de Carbopol 934 à 0.5, 0.75 et 0.1 % en masse. Le comportement rhéologique de ces fluides est décrit par le modèle en loi puissance, τ = KΓ (np-1)/2 γ, avec un indice de rhéofluidification, n p , compris entre 0.4 et 1. Ici, τ est le déviateur du tenseur des contraintes, K, la consistance, γ le tenseur des taux de déformation, et Γ son deuxième invariant. Le dispositif expérimental est cylindrique avec un rapport de diamètre sur épaisseur du fluide égal à Les résultats sont présentés sous forme de corrélations reliant le nombre de Nussselt aux nombres de Rayleigh et de Prandtl dans la gamme 10 5 ≤ Ra ≤ 10 6 . Les nombres de Rayleigh et de Prandtl sont définis en utilisant une viscosité nominale déduite de la mise sous forme adimensionnelle des équations. La première étude théorique portant sur la détermination du nombre de Rayleigh critique pour un fluide rhéofluidifiant, à partir duquel la convection démarre a été menée par [29]. Le comportement rhéologique du fluide est décrit par un modèle en loi puissance. Pour contourner la singularité de ce modèle à cisaillement nul, qui ne permet pas d'effectuer une analyse linéaire de stabilité classique, les auteurs étendent sans justification, l'approche énergétique développée par [30] pour un fluide Newtonien à un fluide en loi puissance. De plus, les auteurs utilisent, les mêmes fonctions propres que celles obtenues en Newtonien ! Moyennant ces approximations, les auteurs trouvent que la valeur critique du nombre de Rayleigh généralisé décroit lorsque l'indice de rhéofluidification diminue. Cette tendance a été vérifiée expérimentalement par les auteurs en utilisant des solutions aqueuses de CarboxyMéthylCellulose (CMC) à différentes concentrations. Par ailleurs, Liang et Acrivos [31] ont étudié expérimentalement la convection de Rayleigh-Bénard pour des solutions aqueuses de polyacrylamide (Separan AP30) à 0.1 et 0.5% en masse. Ces fluides ont l'avantage de présenter clairement un plateau Newtonien pour les très faibles cisaillements. Les auteurs constatent que la valeur critique du nombre de Rayleigh défini avec la viscosité à cisaillement nul est la même que celle pour un fluide Newtonien. Ce résultat n'est pas étonnant puisque l'état de base correspond à une situation où le fluide est au repos i.e à cisaillement nul. Dans [31], les auteurs montrent que le coefficient de transfert de chaleur augmente avec le caractère rhéofluidifiant.

  Les auteurs montrent clairement que la bifurcation est sous critique et que l'amplitude de la perturbation pour le démarrage de la convection augmente avec le caractère rhéofluidifiant. Cette étude a été reprise dans [41] pour un fluide de Carreau dans une cavité ayant un rapport d'aspect de 1 ou de 10. La constante de temps adimensionnelle dans le modèle rhéologique est fixée à λ = 0.4. Dans le même ordre d'idées, Jenny et al [42] se sont focalisés sur les fluides fortement rhéofluidifiants, où la constante de temps adimensionnelle du fluide λ >> 1. Les auteurs définissent, une viscosité caractéristique de telle façon que le nombre de Rayleigh nécessaire pour le démarrage de la convection dans le cas d'une bifurcation sous-critique soit égale à la valeur critique en Newtonien. Une corrélation est proposée pour la viscosité effective : µ e ≈ n 2.2411 λ n-1 . Les études expérimentales de thermoconvection pour des fluides rhéofluidifiants sont rares. Ceci est probablement lié aux difficultés rencontrées dans le choix de fluide modèle où d'autres caractères non-Newtoniens tels que la viscoélasticité ne soient pas dominants. En outre, les fluides non-Newtoniens sont généralement très visqueux. Un écart de température assez important doit donc être appliqué entre les deux parois pour que la convection démarre. Les hypothèses de Boussinesq risquent de ne plus être valables, au démarrage même de la convection. Récemment des études expérimentales [43] [44] ont été effectuées au laboratoire. Les fluides utilisés sont des solutions aqueuses de Xanthane-Gum (polymère semi-rigide) à différentes concentrations. La rhéologie de ces fluides est décrite par le modèle de Carreau. Le dispositif expérimental est cylindrique avec un diamètre de 120 mm et une hauteur (épaisseur de la couche fluide) de 15 mm et 20 mm. Ce faible rapport d'aspect a été imposé par des raisons techniques liées à l'utilisation d'un instrument d'Imagerie par Résonance Magnétique pour mesurer le champ de vitesse à l'intérieur de la cellule. Les motifs de convection sont mis en évidence par ombroscopie. La figure 9 montre un exemple de résultat obtenu pour une solution de Xanthan à 0.1% en masse. Bien que la taille de la cellule soit du même ordre que la période spatiale des structures convectives, elle montre clairement la complexité dans l'organisation des motifs.

Figure 9 -

 9 Figure 9 -Transition hexagones-rouleaux pour une solution de Xanthan à 0.1% avec un rapport d'aspect de 3. (1) ǫ = 0 ; (2) ǫ = 1 ; (3) ǫ = 1.40 ; (4) ǫ =1.44 ; (5) ǫ = 1.45 ; (6) ǫ = 1.46 ; (7) ǫ = 1.48 ; (8) ǫ = 1.50 ; (9) ǫ = 1.52

Rayleigh-

  Bénard convection in a horizontal layer of a non-Newtonian fluid between slabs of arbitrary thickness and finite thermal conductivity is considered. The first part of the paper deals with the primary bifurcation and the relative stability of convective patterns at threshold. Weakly nonlinear analysis combined with Stuart-Landau equation is used. The competition between squares and rolls, as a function of the shear-thinning degree of the fluid, the slabs' thickness and the ratio of the thermal conductivity of the slabs to that of the fluid, is investigated. Computations of heat transfer coefficients are in agreement with the maximum heat transfer principle. The second part of the paper concerns the stability of the convective patterns towards spatial perturbations and the determination of the band width of the stable wavenumber in the neighborhood of the critical Rayleigh number. The approach used is based on the Ginzburg-Landau equations.

Figure 1 . 1 -Figure 1 . 2 -

 1112 Figure 1.1 -Evolution of the critical wave number k c (a) and the critical Rayleigh number Ra c (b) as a function of ratio of the thermal conductivities χ for different values of the slab thickness Λ : (1) Λ = 0.01, (2) Λ = 0.1, (3) Λ = 1, (4) Λ = 10.

Figure 1 . 3 -Figure 1 . 4 -

 1314 Figure 1.3 -Critical degree of shear-thinning, above which the primary bifurcation becomes subcritical, versus χ for rolls (a) and squares (b). (1) Λ = 0.01, (2) Λ = 0.1, (3) Λ = 1, (4) Λ = 10.

Figure 1 . 5 -Figure 1 . 6 -Figure 1 . 7 -Figure 1 . 8 -

 15161718 Figure 1.5 -Stability domains of rolls and squares as a function of χ and α for different values of Λ : (1) Λ = 0.01, (2) Λ = 0.1, (3) Λ = 1, (4) Λ = 10.

Figure 1 . 9 -

 19 Figure 1.9 -Squares, χ = 0.01, Λ = 1. Temperature distribution and velocity vectors in (a) a lateral section which delimits the square cell, y = π/k c and (b) in a vertical diagonal section. Case of Carreau fluid with α = 10 -4 at ǫ = 0.1.

Figure 1 . 10 -

 110 Figure 1.10 -The Nusselt number as a function of the shear-thinning degree α at ǫ = 0.1, χ = 0.01, Pr = 10 and different values of Λ : (1) Λ = 0.01, (2) Λ = 0.1, (3) Λ = 1 and (4)Λ = 10. (thick lines) Rolls, (thin lines) Squares.

Figure 1 . 11 -

 111 Figure 1.11 -Nusselt number vs ǫ at α = 4 × 10 -4 , χ = 10 -2 and four different values of Λ. (thick lines) Rolls, (thin lines) Squares.

Figure 1 .Figure 1 . 12 -

 1112 Figure1.12 shows the variation of α CR-E as a function of the conductivity of the walls and for different values of the thickness Λ. We have also represented the boundaries α S-R and α c (limit of subcritical bifurcation). As expected, α CR-E increasing with decreasing χ. In the Figure1.13, we have represented in the plane (k, Ra) the curves which delimit the stability domain of rolls with respect to (i) cross-roll instability for different α and (ii) Eckhaus instability (which is independent of α). With increasing shear-thinning effects (β/g 1 increases), the CR stability boundary enlarges and becomes less restrictive than Eckhaus instability for α > α CR-E .

Figure 1 . 13 -

 113 Figure 1.13 -Influence of shear thinning effects on stability boundaries for convection rolls as a function of the Rayleigh number Ra and the wavenumber k for Λ = 1, Pr = 10 and χ = 100. (1) Marginal stability curve, (2) Eckhaus boundary, (3) CR boundary for a Newtonian fluid, (4) CR boundary for a shear-thinning fluid with α = 0.2 × α c = 4.3 10 -5 , (5) CR boundary for α = 0.7 × α c = 1.5 10 -4 , (6) CR boundary for α = 0.9 × α c = 1.93 10 -4 . The vertical dotted line is the zigzag (ZZ) boundary. The conductive state labeled Cond., is stable below curve (1).

  Figure 1.14 depicts the domain of stability of the square stationary solution for different values of the shear thinning degree α. The range of stable wavenumber for square patterns decreases with increasing shear-thinning effects, in contrast with the case of roll patterns. 1.7 Numerical solutions of amplitude equations 1.7.1 Numerical simulation The secondary instabilities described in the previous section are studied here by solving numerically the Ginzburg-Landau equations. For the numerical integration of Eqs. (1.51) and (1.52), we employed a Fourier pseudo-spectral method on a square mesh with periodic boundary conditions. The square domain [-L/2, L/2] × [-L/2, L/2] is discretized into N × N uniformly spaced grid points M ℓp = (x ℓ , y p ) with x ℓ = -L/2 + ℓ∆x, (similarly fot y p ), ∆x = ∆y = L/N and N even. Given A M ℓp = A ℓp , ℓ, p = 1, 2, ...N (similarly for B ℓp ), 1.7. NUMERICAL SOLUTIONS OF AMPLITUDE EQUATIONS 55

Figure 1 . 14 -

 114 Figure 1.14 -Influence of shear thinning effects on stability boundaries of squares as a function of the Rayleigh number Ra and the wavenumber k for Λ = 1, Pr = 10 and χ = 10 -3 . (1) Marginal stability curve, (2) Square Eckhaus boundary, (3) RI boundary for a Newtonian fluid, (4) RI boundary for Carreau fluid with α = 0.5 × α S-R = 4.25 × 10 -4 , (5) = α = 0.75 × α S-R = 6.37 × 10 -4 and (6) α = 0.5 × α S-R = 7.65 × 10 -4 . (ZZ) is the zigzag boundary.

Figure 1 . 15 -Figure 1 . 16 -

 115116 Figure 1.15 -Rolls : Points P 1 , P 2 , P 3 , P 4 , where numerical simulations were performed. (a) Low or moderate shear thinning effects, (b) strong shear thinning effects. (M) Marginal stability curve, (E) Eckhaus boundary, (CR) cross-roll boundary, (Z) zigzag boundary.

Figure 1 .Figure 1 . 17 -Figure 1 . 18 -Figure 1 . 19 -Figure 1 . 20 -Figure 1 . 21 -

 1117118119120121 Figure 1.19 shows two representative points, denoted by the symbole (+), where the integration of amplitude equations (1.51) and (1.52) is performed. At point P 5 , we have the following parameters : q = 0.2, ǫ = 0.05, α = 5 × 10 -5 , Λ = 0.01 and χ = 0.01. With these parameters, square solution is RI unstable and SE stable. In addition |D | < D xy . The Figure 1.20 shows the evolution of the convective pattern with time. We have represented f (x, y, t) = A(x, y, t)e ikcx + B(x, y, t)e ikcy + c.c.at four different chosen times. Initially, the convective pattern is a perfect square. A small random perturbation is added in the x-and y-directions. After elimination and adjustment of wavelengths, the structure reaches a stable state with a wavenumber very close to the critical value.At point P 6 , the parameters are q = 0.4, ǫ = 0.05, α = 5 × 10 -5 , Λ = 0.01 and χ = 0.01. In this case, the square solution is RI and SE unstable, furthermore |D | > D xy .The time evolution of the structure is shown in Fig.1.21. It is not surprising that the dynamics is faster that in the previous case, since P 6 is farther for IR stability curve than P 5 . The process of wavelength elimination is also quite different. This could be related to the fact that at P 5 ,|D | < D xy , whereas at P 6 , |D | > D xy .

Figure 1 . 22 -Figure 1 . 23 -

 122123 Figure 1.22 -Variations on Q 1 and Q 2 of the four eigenvalues in the situation where the square pattern is RI unstable only. Here ǫ = 0.1, q = 0.3, α = 0.5 × α c , Λ = 1, χ = 10 -3 . (a) σ 1 , (b) σ 2 , (c) σ 3 , (d) σ 4

Figure 1 . 24 -

 124 Figure 1.24 -The stability domain of ETD2 scheme and the repartition of the points (x a , y a ) and (x b , y b ) along the vertical line.

  In a recent article(Bouteraa et al. J. Fluid Mech., vol. 767, 2015, pp. 696-734) a weakly nonlinear analysis of the Rayleigh-Bénard convection for shear-thinning fluids between two horizontal plates heated from below has been performed. A Landau-Stuart equation is used. Competition between homogeneous rolls, squares and hexagons is investigated. The authors found that in the supercritical regime, only rolls are stable near onset. In this paper, the variation of the viscosity with temperature is taken into account. The variation of the viscosity in the convection layer breaks the up-down symmetry causing a subcritical convection onset to hexagonal pattern. This convection pattern is studied in the framework of amplitude equations. Ginzburg-Landau equations including nonvariational quadratic spatial terms are established using a multiple-scale expansion method. The stability of hexagonal patterns towards homogeneous perturbations (amplitude instabilities) and long wavelength perturbations (phase instabilities) is analyzed for different values of the shearthinning degree α (Eq. 2.11) and the ratio r of the viscosities at the top and bottom walls. The rheological behavior of the fluid is described by the Carreau model and the relationship between the viscosity and the temperature is of exponential type. It is shown that the range of the reduced Rayleigh number ǫ = Ra -Ra c Ra c where both hexagons (with up flow at the center) and rolls are stable decreases with increasing shear-thinning effects, whereas the viscosity contrast has the opposit effect. The stability diagram of hexagons determined from the phase equations is closed for low values of r and becomes open and asymmetric for moderate values of r. Shear-thinning effects reduce the range of stable wave numbers. For the parameters considered here, the stable hexagons domain remains delimited mainly by the transverse phase instability.

. 30 )

 30 These eigenmodes are thus always damped. The first class of eigenmodes corresponds to the eigenvalue problem (2.25) and (2.28) where s is the eigenvalue and X 11 = (F 11 , G 11 )

- 4 .

 4 The minimum of the marginal stability curves gives the critical Rayleigh number Ra c and critical wave number k c . In the case of exponential fluids, figure 2.1 displays the variation of the criticalCHAPITRE 2.

Figure 2 . 1 -Figure 2 . 2 -

 2122 Figure 2.1 -Exponential fluid. Critical Rayleigh number (a) and critical wavenumber (b) as function of the viscosity ratio. (1) NSBC, (2) SFBC.

Figure 2 .

 2 Figure 2.2 displays the profiles of the vertical velocity eigenfunction and the temperature perturbation at the first order for different values of the thermodependency coefficient.

  critical Rayleigh number for the onset of convection determined from the linear stability analysis depends only on the norm k c of the wavevector. Because of the isotropy of the extended horizontal plane, the direction of the wavevector is arbitrary. In addition any linear combination of modes A p exp (ik p • r) (F 11 (z), G 11 (z)) where r = (x, y), k p = (k px , k py ), |k p | = k c and A p 's are constant coefficients is a solution of the linear problem,

  o.t." means "higher order terms". The hexagon patterns are made of three pairs of wavevectors at 2π/3 angles apart : k 1 = k c e x , k 2 = k c -e x /2 + √ 3/2 e y and k 3 = k c -e x /2 -√ 3/2 e y . The objective is to determine the spatio-temporal evolution of the amplitude A p , above threshold, due to different nonlinearities of the problem.

  slow scales X = δx , Y = δy and T = δ 2 t with δ << 1 (2.41) 2.5. AMPLITUDE EQUATIONS IN A HEXAGONAL LATTICE 85 are treated as independent of the fast scales x, y and t. The derivatives with respect to the new variables are

( 1 ) 2 A * 3 1 0G

 121 .55) The forcing terms in the right-hand side of equation (2.55) are computed by introducing the first order solution (2.51)-(2.53). This forcing is composed of four components : Terms proportional to |A p | 2 (p = 1, 2, 3), with the wavenumber modulus |k| = 0, due to the interaction of the eigenmode with its complex conjugate. (2) Terms proportional to A 2 p exp(2ik p • r), |k| = 2k c , due to the interaction of the eigenmode with itself, (3) terms proportional to A p A * q exp (i (k pk q ) • r), |k| = √ 3k c and (4) resonant forcing with wavevector k ℓ (ℓ = 1, 2, 3 and |k ℓ | = k c ). Four separate sets of non homogeneous differential equations are then derived for each component. They are given in Appendix 2.B. For the fourth component, the right-hand side of the non-homogeneous differential contains secular terms. A solvability condition, known as the Fredholm alternative should then be applied for a solution to exist, i.e. the left hand side of equation (2.55) has to be orthogonal to the null-space of the adjoint operator. We obtain A * ad (2F 11 DG 11 + G 11 DF 11 ) dz + (2.56)

2. 5 .

 5 AMPLITUDE EQUATIONS IN A HEXAGONAL LATTICE 87 Two other similar relations are obtained by circular permutation of the indices. The integrals in (2.56) are evaluated numerically by means of the Clenshaw and Curtis method.

Figure 2 . 3 -Figure 2 . 4 -

 2324 Figure 2.3 -Variation of the characteristic time τ 0 (a) and the coherence length ξ 0 (b) as a function of the viscosities ratio, for different values of Prandtl number. (1) P r = 100 ; (2) P r = 50 ; (3) P r = 10 ; (4) P r = 5 ; (5) P r = 2 ; (6) P r = 1.

Figure 2 . 5 -Figure 2 . 6 -

 2526 Figure 2.5 -(a) "Newtonian" and (b) non-Newtonian contribution to the first Landau coefficient and to the cross-saturation coefficient (c) and (d) respectively as a function of r for different values of P r. (1) P r = 100 ; (2) P r = 50 ; (3) P r = 10 ; (4) P r = 5 ; (5) P r = 2 ; (6) P r = 1.

Figure 2 . 7 -

 27 Figure 2.7 -Amplitude stability curves in (ǫ, q) plane at r = 2.5 and two different values of α : (a) α = 0, (b) α = 10 -4

Figure 2 . 8 -

 28 Figure 2.8 -Bifurcation diagram for hexagons in the case where g 2 > g 1 , r = 2.5 and α = 1.5 10 -4 . The amplitude |A 1 | is plotted against the distance to the threshold ǫ, for the roll-solution branch (labeled R), for the mixed mode branch (labeled M ) and for the two hexagon-solution branches, up-and down-hexagons. Solid lines indicate stable solutions and dashed lines represent unstable solutions.

Figure 2 . 9 -Figure 2 . 10 -

 29210 Figure 2.9 -Variations of ǫ a , ǫ r , ǫ h and (ǫ hǫ r ) versus r for three values of the shearthinning degree α. The Prandtl number is fixed, P r = 100. (1) Newtonian fluid, α = 0 ; (2) Carreau fluid with α = 1.07 10 -4 ; (3) Carreau fluid with α = 1.43 10 -4 . The dashed line is the correlation proposed par [9] for a Newtonian fluid.

Figure 2 . 11 -Figure 2 . 12 -

 211212 Figure 2.11 -Hexagon stability diagram for a Newtonian fluid at P r = 100 and two different values of the viscosities ratio : (a) r = 1.5 and (b) r = 2.5. Curve (1) : D ⊥ = 0, curve (2) D =0 , curve (3) bifurcation from the conductive state to convection with hexagons, curve(4) bifurcation from hexagons to rolls.

) 31 = 2 HX 2 .

 3122 µ b 2 (∇ Hx • ∇ HX ) B Second-order solution (hexagons) 2.B.1 Solution proportional to |A p | 2 (zero mode) The first component of the second order solution, proportional to |A p | 2 , provides a correction of the basic state. Considering the w-equation, it is shown that the factor of |A 1 | 2 , |A 2 | 2 and |A 3 | 2 in the non linear inertial N I

  109) Similar equations are obtained for the amplitudes A 2 and A 3 . The boundary conditions on W 41 and T 41 are the same as the ones on (F 11 , G 11 ).

Bibliographie 3 . 1 Introduction

 31 Les approches faiblement non linéaires mises en oeuvre dans les chapitres précédents permettent d'étudier la compétition et la stabilité des motifs de convection à proximité du seuil de l'instabilité. En revanche, elles se révèlent inappropriées pour étudier les structures dès lors qu'on s'éloigne trop des conditions critiques. En effet, elles conduisent à considérer des développements asymptotiques à des ordres de plus en plus élevés et donc des équations de plus en plus complexes au fur et à mesure que la convection gagne en intensité et que les effets non linéaires deviennent de plus en plus importants. La difficulté calculatoire qui en résulte s'avère rédhibitoire et on est contraint d'employer une approche alternative appelée procédure de continuation pour étudier l'évolution des structures. Cette méthode, fréquemment utilisée dans les problèmes de bifurcation, consiste à suivre de proche en proche la branche de solutions stationnaires lorsque le paramètre de contrôle du problème varie. Dans ce chapitre, les principes de cette technique sont tout d'abord exposés avant de l'appliquer au cas d'un fluide rhéofluidifiant non thermodépendant où les rouleaux sont les seuls motifs de convection stables[1]. On s'intéresse notamment à l'influence de la rhéofluidification d'une part sur les champs de vitesse, de température et de viscosité et d'autre part sur les échanges de chaleur dans le fluide. Le choix des paramètres numériques utilisés ainsi que les performances en terme de temps de calcul sont également discutés.

Figure 3 . 1 - 8 ) 2 ( 3 . 9 ) 3 . 12 )

 318239312 Figure 3.1 -Configuration de Rayleigh-Bénard.

. 20 )

 20 En introduisant les expressions de w S et de θ S (Eq.(3.15)) dans les équations (3.3) et (3.4), on aboutit à un système d'équations algébriques non linéaires pour chaque mode de Fourier et en chaque point de Gauss-Lobatto. Ce système d'équations peut être écrit formellement sous la forme : F(X, Ra) = 0 avec F : R NT × R → R NT où NT = 2(N + 1)(M + 1) (3.21) où X est le vecteur des coefficients spectraux de dimension NT : X = [a 0 , a 1 , ...a N , b 0 , b 1 , ...b N ] (3.22) La fonction F peut se décomposer en : -termes linéaires L w , L θ , -termes non linéaires d'inertie NI w , NI θ , CHAPITRE 3.

  27) avec i = 0, ..., N d -1 et j = 0, ..., M d -1 où N d et M d désignent les nombres de points dans les directions tangentielle et normale aux parois. On a donc une discrétisation uniforme dans la direction horizontale et une discrétisation non-uniforme utilisant les points de Gauss-Lobatto dans la direction verticale. Les nombres de points N d et M d doivent au minimum respecter la règle d'Orszag permettant d'éviter les erreurs de repliement spectral (erreurs d'aliasing) lors de la

43 )dx

 43 La démarche précédente consistant à évaluer des termes dans l'espace physique avant de les projeter dans l'espace spectral, peut être appliquée aussi bien aux termes d'inertie qu'aux termes visqueux. Cependant, elle se révèle très coûteuse en temps de calcul et requiert parfois un grand nombre de points pour être précise. En l'utilisant pour projeter les termes d'inertie à la place de la procédure donnée dans 3.3.3, on perd en précision et 3.4. NOMBRE DE NUSSELT 117 on augmente fortement la durée des calculs. C'est pour ces raisons qu'on la réserve au traitement des termes visqueux où il n'existe pas d'alternative.3.4 Nombre de Nusselt Le nombre de Nusselt Nu est la grandeur adimensionnelle qui caractérise l'influence de la convection sur les échanges de chaleur dans la couche de fluide. C'est le rapport entre le flux de chaleur total et le flux de chaleur dans l'état conductif de base. Le champ de température étant périodique suivant la direction horizontale de période Lc = 2π k c d , on peut calculer le nombre de Nusselt en considérant uniquement une longueur d'onde. A la paroi inférieure, le nombre de Nusselt est défini par : température T correspond à la superposition de la perturbation en température au champ conductif de l'état de base soit T = T b + θ où T b (z) = -z + T 0 . Le nombre de Nusselt s'écrit alors : Il est possible d'évaluer l'intégrale intervenant dans l'expression du nombre de Nusselt en utilisant les valeurs des variables en différents points de l'écoulement i.e. pour différentes valeurs de z comprises entreen rendre compte, on peut réécrire l'équation de l'énergie (3.4) en tenant compte du fait que le fluide est incompressible :

Figure 3 . 2 -

 32 Figure 3.2 -Les différentes étapes des méthodes de continuation.

Figure 3 . 3 -

 33 Figure 3.3 -Points de bifurcation (B) et de retournement (R 1 et R 2 )

  où (∂ X F|∂ Ra F) est l'opérateur matriciel de dimensions (NT, NT+1) obtenu par concaténation des deux opérateurs ∂ X F et ∂ Ra F. Le système (3.50) faisant intervenir NT équations pour NT + 1 inconnues, on doit imposer une condition de normalisation pour que T soit déterminé de manière unique. De manière très générale, une telle condition s'écrit c tr T = 1 où c est un vecteur ≪ approprié ≫ et où tr désigne l'opération de transposition. Un choix possible est de contraindre une composante du vecteur tangent à valoir 1 : on prend alors c = e k où e k désigne le vecteur unitaire de taille NT dont toutes les composantes sont nulles à l'exception de la k ième qui vaut 1 soit : (e tr k )T = T k = 1 (3.51) En associant (3.50) et (3.51), on détermine le vecteur tangent par résolution du système linéaire : ∂ X F|∂ Ra F e tr k T = e NT+1 (3.52) A condition que rang(∂ X F|∂ Ra F) = NT, le vecteur tangent est bien défini le long de la branche. En particulier, il est possible de le déterminer aux points de retournement ce que ne permettent ni (3.46) ni (3.49).

  .55) Dans cette approche, Ra devient donc une inconnue du problème fonction de γ. Le système F(X, Ra) = 0 admet alors NT équations pour NT + 1 inconnues (X, Ra). Il est alors nécessaire d'ajouter une équation supplémentaire liant X, Ra et γ. Une telle équation est appelée équation de paramétrisation et a la forme générale suivante : p(X, Ra, γ) = 0 (3.56) En l'associant avec (3.44), on obtient un système étendu qui admet bien autant d'équations que d'inconnues : F (X, Ra, γ) = F(X, Ra) p(X, Ra, γ) = 0 (3.57) Pour résoudre le système étendu (3.57), il reste à définir une paramétrisation. Une manière très courante de procéder est de choisir la longueur d'arc s comme paramètre. Cette dernière vérifie la relation : l'inconvénient de faire intervenir les dérivées de X et de Ra. Mathématiquement, utiliser (3.58) comme équation de paramétrisation conduit à faire de F (X, Ra, s) = 0 un système de NT + 1 équations à 2(NT + 1) inconnues puisqu'on doit considérer à la fois les inconnues et leurs dérivées. Pour cette raison, on ne travaille pas directement avec (3.58) mais avec une approximation de cette dernière. s correspond alors à une approximation de la longueur d'arc et on parle de continuation par pseudo longueur d'arc. Cette technique introduite par Keller en 1977 [9] est très fréquemment utilisée dans les problèmes de continuation . CHAPITRE 3.

  60) Elle présente l'avantage de ne pas faire intervenir de dérivées : c'est pour cette raison qu'elle a été choisie dans le code numérique mis en place. La plupart des autres équations de paramétrisation pour la longueur d'arc sont basées sur (3.59) et (3.60) : elles introduisent un facteur d'atténuation 0 ≤ ξ ≤ 1 qui permet de donner plus de poids à X ou à Ra dans l'équation de paramétrisation. Par souci de simplicité, ces équations de paramétrisation n'ont pas été retenues pour le code numérique.

3. 6 . 3 . 6 . 1 3 .Figure 3 . 4 -

 6361334 Figure 3.4 -Ballon de Busse pour Pr → ∞

Figure 3 .

 3 5 montre ainsi ses variations avec l'écart relatif au seuil ǫ = Ra -Ra c Ra c pour différentes valeurs de N (des valeurs numériques précises sont répertoriées dans le tableau 3.1). On constate que pour 0 ≤ ǫ ≤ 2, les nombres de Nusselt obtenus sont identiques quelque soit le nombre de modes de Fourier utilisés : 4 modes sont donc suffisants pour décrire les champs de vitesse et de température dans cette gamme. Pour des écarts au seuil plus importants, les nombres de Nusselt ne convergent vers une valeur fixe qu'au delà d'un certain nombre de modes : par exemple, pour ǫ = 4 et ǫ = 8, il faut respectivement 5 et 6 modes d'après le tableau 3.1. On constate que ce nombre minimal de modes de Fourier à considérer pour avoir convergence augmente quand l'écart relatif au seuil s'accroît. Pour converger, la procédure de continuation requiert donc d'autant plus de modes de Fourier que l'écart relatif au seuil est important : cela s'explique par le fait que les harmoniques issus des interactions non linéaires entre modes de Fourier deviennent de plus en plus importants au fur et à mesure que la convection gagne en intensité.

3. 6 . 3

 63 Évolution des différents modes de Fourier avec l'écart au seuil de la convection On s'intéresse ici aux évolutions des différents modes de Fourier avec ǫ, l'écart relatif au seuil de la convection. Comme précédemment, on se limite au cas Pr = 10 et on considère uniquement l'intervalle 0 ≤ ǫ ≤ 6 où 7 modes de Fourier (N=6) sont suffisants pour garantir la convergence de la procédure de continuation comme développé à la section

9 Figure 3 . 5 -

 935 Figure 3.5 -Convergence de la procédure de continuation pour Pr = 10 pour un fluide Newtonien

Figure 3 . 6 -Figure 3 . 7 -

 3637 Figure 3.6 -Modes de Fourier associés à w pour un fluide Newtonien

et 3 . 9 .

 39 Ces dernières montrent notamment une saturation du mode fondamental θ 1 : l'amplitude et l'énergie de ce dernier augmentent jusqu'à ǫ = 1 avant de décroître progressivement pour des valeurs plus importantes de l'écart relatif au seuil. Les amplitudes et énergies des autres modes ne font, quant à elles, qu'augmenter avec ǫ. Avant de commencer à saturer, on observe que c'est le mode θ 1 qui présente les plus grandes amplitude et énergie. En revanche, une fois la saturation amorcée, c'est le mode θ 0 qui devient dominant. Son amplitude et son énergie continuent ensuite d'augmenter fortement avec ǫ. En conséquence, la prédominance de θ 0 sur les autres modes ne fait alors que se renforcer quand on s'éloigne du seuil de la convection. Seuls θ 0 et θ 1 contribuent significativement à l'énergie du champ de température : les participations des autres modes θ l , l = 2...6 demeurent négligeables.

Figure 3 . 8 - 1 ≤ 6 Figure 3 . 9 -

 381639 Figure 3.8 -Modes de Fourier associés à θ pour un fluide Newtonien

Figure 3 . 10 -

 310 Figure 3.10 -Champs de vitesse w pour un fluide Newtonien

Figure 3 . 11 -

 311 Figure 3.11 -Champs de vitesse u pour un fluide Newtonien

Figure 3 . 12 - 1 .

 3121 Figure 3.12 -Champs de température θ pour un fluide Newtonien

= 7 ×

 7 .71) Pour des nombres de Prandtl Pr ≥ 1 et des parois rigides, il est montré que la bifurcation est : -surcritique si α ≤ α c = 2.15 × 10 -4 , -sous-critique si α > α c = 2.15 × 10 -4 . Dans le cas de parois libres, M. Bouteraa et al. trouvent que la transition d'une bifurcation surcritique à sous-critique a lieu pour α c = 1.025 × 10 -4 ce qui est en accord avec la valeur trouvée par N.J. Balmforth et A.C. Rust dans [19]. En revanche, elle ne correspond pas à celle obtenue par B. Albaalbaki et R.E. Khayat dans [20] qui est de α AK c 10 -4 . Les trois études s'accordent donc sur l'existence d'une bifurcation sous-critique pour un caractère rhéofluidifiant suffisamment marqué mais s'opposent quant à la valeur de α c nécessaire. En ce qui concerne les motifs de convection, M.Bouteraa et al. montrent que seuls les rouleaux sont stables à proximité du seuil de la convection et ce quelque soit la rhéologie du fluide, infirmant ainsi les résultats de B. Albaalbaki et R.E. Khayat.

3 . 4 .

 34 Comme mentionné précédemment, cette projection nécessite d'abord l'évaluation des termes visqueux sur une grille de l'espace physique constituée de M d points de Gauss-Lobatto dans la direction verticale et de N d points équirépartis dans la direction horizontale. En fonction de la finesse du maillage, la projection sera plus ou moins précise : une série de simulations a donc été menée afin de déterminer les nombres de points M d et N d nécessaires à une décomposition rigoureuse dans l'espace de Fourier Chebyshev. Ces simulations ont été menées en considérant un fluide rhéofluidifiant pour lequel n c = 0.5 et λ = 0.02 et en utilisant 7 modes de Fourier (N = 6). Pour raffiner le maillage, on a choisi d'augmenter uniquement N d ; le nombre de points de Gauss-Lobatto dans la direction verticale a été fixé à M d = 42.

3. 7 .Figure 3 . 14 - 3 . 7 . 3

 7314373 Figure 3.14 -Convergence nombre de points N d pour n c = 0.5, λ = 0.02

Figure 3 . 15 - 5 Figure 3 . 16 - 6 Figure 3 . 17 - 3 .Figure 3 . 18 - 6 Figure 3 . 19 -

 3155316631733186319 Figure 3.15 -Convergence nombre de modes pour n c = 0.5, λ = 0.02

Figure 3 . 20 -Figure 3 . 21 -Figure 3 . 22 -

 320321322 Figure 3.20 -Champs de vitesse w pour n c = 0.5, λ = 0.02

Figure 3 . 23 -

 323 Figure 3.23 -Champs de viscosité µ pour n c = 0.5, λ = 0.02

Figure 3 .

 3 Figure 3.24 montre ainsi que µ P et µ M décroissent d'autant plus rapidement avec l'écart relatif au seuil que n c est faible. Pour de faibles valeurs de n c , les frottements visqueux sont donc moindres et en conséquence, la convection et les échanges de chaleur induits par cette dernière sont supérieurs.

Figure 3 . 26 - 3 . 3 . 7 . 6 5 λ 3 Figure 3 . 27 -

 326337653327 Figure 3.26 -Tentatives d'obtention d'une courbe maîtresse en utilisant les viscosités µ P et µ M .

-

  A court terme, on souhaite complexifier le modèle faiblement non linéaire mis en place pour décrire la convection d'un fluide thermodépendant à proximité du seuil en considérant un développement au second ordre de la masse volumique avec la température. En effet, les fluides rhéofluidifiants utilisés expérimentalement sont très visqueux et les écarts de température nécessaires pour déclencher la convection sont importants dans la pratique (de l'ordre d'une dizaine de degrés). En conséquence, on peut s'attendre à ce que les hypothèses de Boussinesq qui Résumé de la thèse en français On désigne par convection thermique le phénomène physique par lequel un fluide se met en mouvement sous l'effet d'un gradient thermique. Pour l'étudier, on considère comme configuration de référence un fluide placé entre deux plaques horizontales chauffé par le bas et refroidi par le haut dans le champ de gravité. Pour un gradient thermique suffisamment important, le fluide initialement au repos se met en mouvement. Les premiers à étudier expérimentalement et théoriquement cette transition d'un état conductif à un état convectif dans cette configuration furent Lord Rayleigh et Henri Bénard ; aussi on parle de convection de Rayleigh-Bénard pour qualifier ce phénomène. Henri Bénard a notamment montré que la convection du fluide se faisait de manière ordonnée, on constate ainsi l'émergence de différents motifs de convection tels des rouleaux, des carrés ou des hexagones. L'origine de la convection fut quant à elle donnée par Lord Rayleigh. Ce dernier montre que, dans la configuration utilisée, deux mécanismes antagonistes entrent en jeu : -D'un côté, les forces de poussée d'Archimède induites par les variations de la masse volumique avec la température tendent à faire s'élever les couches de fluides basses plus chaudes donc moins denses. -De l'autre, les forces de frottement visqueux et la conduction thermique s'opposent à la mise en mouvement du fluide. Le rapport de force entre ces deux phénomènes peut être quantifié par une grandeur adimensionnelle appelée nombre de Rayleigh Ra. Ainsi pour une valeur suffisamment grande de ce paramètre, les forces de poussée d'Archimède l'emportent sur les forces de friction et la diffusion thermique ce qui se traduit par le déclenchement de la convection : on parle alors de bifurcation primaire. En considérant des surfaces libres comme conditions limites, Lord Rayleigh a déterminé théoriquement la valeur Ra c à partir de laquelle le fluide se met en mouvement ainsi que le nombre d'onde k c des structures qui se développent. En s'éloignant des conditions critiques, les différents motifs de convection vont à leur tour subir des instabilités qualifiées de ≪ secondaires ≫ (par opposition à la bifurcation primaire) qui vont engendrer une modulation spatiale des structures. Pour les rouleaux de convection et dans la limite des grands nombres de Prandtl, la littérature répertorie à proximité des conditions critiques (k c , Ra c ) trois mécanismes d'instabilités distincts : l'instabilité d'Eckhaus, l'instabilité Zigzag et l'instabilité rouleaux croisés. L'instabilité d'Eckhaus est une instabilité de grande longueur d'onde qui se manifeste par une dilatation ou une compression des rouleaux de convection. L'instabilité Zigzag est elle aussi une 158 instabilité de grande longueur d'onde mais elle opère dans la direction perpendiculaire à l'axe des rouleaux donnant ainsi lieu à l'apparition d'ondulations. L'instabilité rouleaux croisés, quant à elle, se traduit par le développement progressif d'un réseau de rouleaux perpendiculairement à un premier ensemble de motifs. Comparativement aux rouleaux, les instabilités secondaires des autres structures de convection (carrés, hexagones) ne sont que peu étudiées. Que ce soit au niveau de la bifurcation primaire ou des instabilités secondaires, les principaux travaux de la littérature font intervenir des fluides newtoniens. Or dans de nombreux domaines tels que la géologie ou les industrie alimentaire et pétrochimique, les fluides rencontrés ont un comportement rhéologique plus complexe. Le cas le plus fréquemment rencontré est celui où les fluides ont une viscosité qui diminue avec le taux de cisaillement ce qui correspond au comportement rhéofluidifiant. L'objectif de la thèse est ainsi de déterminer l'influence d'une telle rhéologie sur les motifs de convection et leur stabilité à la fois à proximité du seuil de la convection et loin des conditions critiques. La thèse compte trois chapitres : les deux premiers sont consacrés aux instabilités secondaires des motifs de convection à proximité des conditions critiques tandis que le dernier traite du devenir des structures loin du seuil de la convection. Dans le premier chapitre de la thèse, on étudie le cas où les plaques confinant le fluide ont une conductivité thermique finie. Les équations descriptives du comportement du fluide à savoir les équations de conservation de quantité de mouvement, de continuité et de la chaleur sont écrites dans le cadre des approximations de Boussinesq. Celles ci stipulent que :

  Les plaques sont considérées rigides et le comportement rhéofluidifiant est décrit par le modèle de Carreau. Une étude linéaire de stabilité est d'abord menée pour déterminer les conditions critiques (k c , Ra c ) à partir desquelles la convection s'amorce. Pour ce faire, on considère des perturbations infinitésimales et on linéarise les équations descriptives du système. On cherche ensuite une solution sous forme de modes normaux ce qui permet d'aboutir à un problème aux valeurs propres. Ce dernier est résolu numériquement en utilisant une méthode spectrale reposant sur les polynômes de Chebyshev et les points de collocation de Gauss-Lobatto. En cherchant les couples (k, Ra) qui annule le taux de croissance des perturbations, on détermine la courbe de stabilité marginale Ra(k) dont le minimum correspond aux conditions critiques. L'évolution de ces conditions critiques en fonction de la conductivité des plaques a été déterminée pour différentes épaisseurs de ces dernières. Signalons que la rhéologie n'intervient pas à l'ordre linéaire. Dans la limite BIBLIOGRAPHIE 159 de plaques parfaitement conductrices, on retrouve les valeurs k c = 3.11, Ra c = 1707.76 bien connues dans la littérature. Quand la conductivité des plaques diminue ou que leur épaisseur augmente, les conditions critiques diminuent jusqu'à valoir k c = 0 et Ra c = 720 dans le cas de plaques parfaitement isolantes. Physiquement, de telles évolutions s'expliquent par le fait que dans le cas de parois isolantes, les perturbations en température ne sont pas éliminées aux parois internes ce qui se traduit par un effet d'accumulation des perturbations dans la couche de fluide. De ce fait, la bifurcation primaire a besoin de moins d'énergie pour se produire ce qui explique la diminution de Ra c avec la réduction de la conductivité des plaques. Pour étudier les motifs de convection au seuil de l'instabilité, on utilise une approche multi-échelle. Cette dernière consiste à décomposer les différentes variables du problème en série d'un petit paramètre δ. La substitution de ces décompositions dans les équations descriptives du système conduit à un ensemble d'équations différentielles ordinaires pour chaque puissance de δ que l'on vient résoudre. L'application d'une condition de compatibilité (alternative de Fredholm) permet de déterminer les équations de Ginzburg-Landau décrivant la dynamique du système à proximité des conditions critiques. En menant une analyse linéaire de stabilité sur ces équations, on est en mesure de déterminer la nature du motif émergeant à la bifurcation primaire ainsi que les différentes instabilités secondaires subies par ce dernier. En ce qui concerne les structures émergeant au seuil de la convection, il est montré que le caractère isolant favorise l'émergence de motifs en carrés tandis que le comportement rhéofluidifiant avantage les rouleaux de convection. Ces motifs subissent différentes instabilités secondaires. Pour les rouleaux, trois instabilités viennent limiter la gamme de nombres d'onde stables : il s'agit des instabilités d'Eckhaus, Zigzag et rouleaux-croisés. A proximité du seuil de la convection, il est montré que la rhéologie a une influence uniquement sur l'instabilité rouleaux croisés : en effet, la gamme de nombres d'onde stables est d'autant plus étendue que le caractère rhéofluidifiant est prononcé. Contrairement au cas Newtonien, il est montré que l'instabilité d'Eckhaus peut être plus restreignante que rouleaux-croisés à condition que l'effet rhéofluidifiant soit suffisamment marqué. Pour les carrés, on met en évidence deux instabilités secondaires : l'instabilité d'Eckhaus et l'instabilité rectangulaire. Le comportement rhéofluidifiant n'a d'influence que sur l'instabilité rectangulaire et limite d'autant plus la gamme de nombre d'ondes que la rhéofluidification est forte. Par ailleurs, on constate que l'instabilité rectangulaire est toujours dominante devant Eckhaus. Ce premier chapitre se termine par une simulation numérique directe des équations de Ginzburg-Landau qui permet d'une part de confirmer les résultats précédents et d'autre part d'observer l'évolution temporelle des diverses instabilités secondaires.Dans le deuxième chapitre de la thèse, on s'intéresse aux motifs de convection qui se développent dans la convection de Rayleigh-Bénard lorsqu'on considère un fluide thermodépendant et rhéofluidifiant. La thermodépendance est prise en compte en considérant une variation exponentielle de la viscosité du fluide avec la température tandis que le comportement rhéofluidifiant est décrit par le modèle de Carreau comme dans le chapitre précédent. L'approximation de Boussinesq est utilisée pour la masse volumique. On adopte la même démarche que précédemment : on commence par écrire les équations de conservation de la quantité de mouvement, de continuité et de l'énergie puis on procède à une analyse linéaire de stabilité. On détermine ainsi l'évolution des conditions critiques avec le rapport r des viscosités aux plaques supérieure et inférieure. L'analyse linéaire de stabilité révèle que (i) pour de faibles rapports de viscosité 1 ≤ r ≤ 1.5, k c et Ra c sont presque constants, (ii) pour des valeurs de r modérées 1.5 ≤ r ≤ 8, Ra c augmente avec r tandis que k c varie peu, (iii) pour des rapports de viscosités importants, la convection est gouvernée par une sous-couche plus instable que la couche complète ce qui conduit à une diminution de Ra c avec r et une augmentation de k c avec r. L'étude des motifsde convection et des instabilités secondaires se fait grâce à une approche multi-échelles semblable à celle du chapitre précédent. A proximité du seuil de la convection, on montre que deux types de structures entrent en compétition : ce sont les rouleaux de convection et les motifs hexagonaux. Pour la bifurcation primaire, on met en évidence différents domaines d'existence des structures en fonction de l'écart relatif au seuil ǫ = Ra -Ra c Ra c . Ainsi, on montre que pour de faibles valeurs ǫ soit ǫ a ≤ ǫ ≤ ǫ r , les hexagones sont les seules structures stables. Signalons que ǫ a est négatif ce qui indique que la bifurcation primaire est sous critique. Pour ǫ r ≤ ǫ ≤ ǫ h , on a un domaine de coexistence des deux motifs. Enfin quand ǫ ≥ ǫ h , les hexagones perdent leur stabilité et il ne subsiste que des rouleaux. L'évolution de la taille de ces différents domaines a été étudiée en fonction des caractères plus ou moins marqués de la thermodépendance et de la rhéofluidification. En faisant croître le rapport des viscosités r, on constate que ǫ r et ǫ h augmentent. Les rouleaux de convection apparaissent donc à partir d'une valeur supérieure de l'écart relatif au seuil tandis que les hexagones perdurent jusqu'à des valeurs plus élevées de ǫ : la thermodépendance favorise donc les motifs hexagonaux au détriment des rouleaux. A l'inverse, en augmentant les effets rhéofluidifiants, on constate que ǫ r et ǫ h diminuent : la rhéofluidification favorise donc les rouleaux de convection aux dépens des hexagones. L'impact de la rhéologie sur les instabilités secondaires et la gamme de nombres d'ondes stables est étudié en utilisant, comme dans le chapitre précédent, le formalisme des équations de Ginzburg-Landau. On se place dans la limite des perturbations de grandes longueurs d'ondes et on procède à une élimination adiabatique. Une telle approche consiste à mettre en évidence les variables lentes qui dirigent la dynamique du système et les variables rapides qui s'adaptent automatiquement aux variations des variables lentes (elles sont≪ éliminées adiabatiquement ≫). On montre que les phases des perturbations sont des variables lentes tandis que les amplitudes et la phase globale sont des variables rapides. La dynamique du système se réduit alors à deux équations de diffusion de phase qui permettent de déterminer le domaine de stabilité des hexagones, délimité par les instabilités secondaires, dans le plan (q, Ra) où q désigne l'écart au nombre critique. Les calculs montrent que ce domaine de stabilité est d'autant plus vaste que le caractère thermodépendant est important et qu'il est d'autant plus réduit que la rhéofluidification BIBLIOGRAPHIE 161 est forte : la thermodépendance et la rhéofluidifiaction jouent donc des rôles antagonistes. Le dernier chapitre est consacré à l'étude des motifs loin du seuil de la convection. Un fluide purement rhéofluidifiant confiné entre deux plaques rigides est considéré. Dans cette configuration, les motifs se développant sont des rouleaux de convection. L'objectif de ce chapitre est d'étudier le devenir de ces structures quand on s'éloigne des conditions critiques. Pour cela, on a recours à une méthode de continuation. Cette dernière est un processus itératif qui permet de déterminer de proche en proche les solutions stationnaires des équations descriptives du problème au fur et à mesure que le paramètre de contrôle, ici le nombre de Rayleigh Ra, varie. Les différentes étapes de cette procédure (prédiction, paramétrisation, correction et contrôle du pas) sont présentées. Les solutions du problème sont cherchées en séries de Fourier et les différents modes associés sont décomposés sur une base de polynômes de Chebyshev. Les termes linéaires et non linéaires d'inertie des équations descriptives du système sont calculés directement dans l'espace spectral de Fourier-Chebyshev tandis que les termes non linéaires visqueux sont d'abord évalués dans l'espace physique puis projetés dans l'espace spectral. Le code numérique mis en place permet de déterminer les champs de vitesse et de température ainsi que les échanges de chaleur via le calcul du nombre de Nusselt Nu. Le cas d'un fluide Newtonien est d'abord considéré. Une étude de convergence est menée pour déterminer le nombre adéquat de modes de Fourier : on constate qu'il est logiquement d'autant plus élevé qu'on s'éloignedes conditions critiques. On montre que le champ de vitesse croît d'abord rapidement à proximité du seuil puis que ses augmentations sont plus modérées. Au fur et à mesure que la convection gagne en intensité, on constate que les élévations de température se concentrent dans la partie supérieure de la couche de fluide tandis que les réductions ont lieu à proximité de la plaque chaude. L'évolution du nombre de Nusselt avec l'écart relatif au seuil ǫ est déterminée pour différentes valeurs du nombre de Prandtl 1 ≤ Pr ≤ ∞. Pour de faibles valeurs de l'écart relatif au seuil 0 ≤ ǫ ≤ 2, on observe que toutes les courbes sont confondues ce qui signifie que le nombre de Prandtl n'a pas d'impact sur les échanges de chaleur à proximité des conditions critiques. Loin de ces dernières, les courbes Nu = f(ǫ) tendent vers une courbe unique quand Pr ≥ 10 ce qui permet de mettre en évidence une limite des grands nombres de Prandtl. Les nombres de Nusselt obtenus se révèlent en très bon accord avec ceux issus de divers travaux de la littérature ce qui permet de valider le code en Newtonien. On s'est ensuite intéressé au cas de fluides rhéofluidifiants. Comme pour le cas Newtonien, plusieurs simulations sont d'abord menées en faisant varier la finesse du maillage de l'espace physique ainsi que le nombre de modes de Fourier pour déterminer les conditions de convergence. On observe que plus les écarts relatifs au seuil sont importants et plus la finesse du maillage doit être grande et le nombre de modes de Fourier important pour s'assurer de la convergence de la méthode. Pour une valeur donnée de l'écart relatif au seuil, on observe que les amplitudes des champs de vitesse et de température sont plus importants pour un fluide rhéofluidifiant. En effet, les forces de frottement sont moindres pour un tel fluide ce qui engendre une convection plus forte. Le champ de viscosité présente une organisation spécifique. On distingue ainsi des zones en forme d'≪ anneaux ≫ où les cisaillements sont très faibles et qui sont donc le siège de fortes viscosités. Ces dernières sont entourées de zones fortement cisaillées où règnent de faibles viscosités. On observe que les échanges de chaleur sont d'autant plus importants que le comportement rhéofluidifiant est prononcé ce qui semble cohérent puisqu'une telle rhéologie favorise la convection et donc les échanges de chaleur induits par cette dernière. Conformément à ce qu'avaient prédites différentes approches faiblement non linéaires, la bifurcation de l'état conductif à l'état convectif devient sous critique si les effets rhéofluidifiants sont suffisamment forts. Signalons que dans son état actuel, le code ne permet pas de suivre la branche sous critique au delà de son point de retournement. En conclusion, la thèse aura permis dans un premier temps de mettre en évidence l'influence d'un effet rhéofluidifiant sur les instabilités secondaires des motifs de convection à proximité des conditions critiques. Dans un second temps, elle aura permis de poser les jalons nécessaires à l'étude des structures convectives loin du seuil de la convection à travers la détermination des branches de solutions stationnaires des équations descriptives du comportement d'un fluide purement rhéofluidifiant. Plusieurs perspectives s'offrent à la suite de ce travail. En premier lieu, on pourra améliorer le modèle mathématique développé dans l'approche faiblement non linéaire en considérant des variations à l'ordre cubique de la masse volumique avec la température. Cela devrait permettre de mieux prendre en considération les gradients de températures importants qui peuvent exister dans les dispositifs expérimentaux faisant intervenir des fluides rhéofluidifiants. Dans le cas de fluides très fortement thermodépendants, différents articles de la littérature rapportent l'existence de structures en carrés : un autre enjeu serait alors d'étudier la compétition entre ces derniers et les hexagones qui existent déjà lorsque les effets thermodépendants sont moindres. En second lieu, on cherchera à développer l'approche fortement non linéaire. On pourra commencer par étudier la stabilité des branches de solutions stationnaires mises en évidence et observer l'influence du caractère rhéofluidifiant sur le domaine de stabilité des rouleaux de convection. On pourra poursuivre en mettant en place une procédure permettant de franchir le point de retournement des branches sous critiques de solutions ce qui nous permettra d'avoir accès à l'intégralité des branches. Enfin, on pourra étendre l'approche fortement non linéaire au cas de fluides rhéofluidifiants et thermodépendants et étudier la stabilité d'autres motifs de convection comme les hexagones. On pourra ainsi apporter un nouvel éclairage au phénomène d'≪ hexagones ré-entrants ≫ à savoir l'existence d'hexagones loin des conditions critiques ce que ne prédit pas une approche faiblement non linéaire. Instabilités secondaires dans la convection de Rayleigh-Bénard pour des fluides rhéofluidifiants. Dans la configuration de Rayleigh-Bénard, on considère une fine couche de fluide placée entre deux parois horizontales, chauffée par le bas et refroidie par le haut. Cette couche peut être le siège d'une instabilité si le gradient thermique est suffisamment important : on passe alors de l'état conductif à l'état convectif et on parle de bifurcation primaire pour qualifier cette première transition. Cette mise en mouvement du fluide se fait de manière ordonnée : on constate ainsi l'émergence de différents motifs de convection comme des rouleaux, des carrés ou encore des hexagones. Ces structures vont à leur tour subir des instabilités qualifiées de ≪ secondaires ≫ qui vont limiter la gamme de nombres d'onde stables. On étudie ici théoriquement ces instabilités d'une part à proximité du seuil de la convection grâce à une approche faiblement non linéaire, d'autre part loin des conditions critiques grâce à une approche fortement non linéaire. Le fluide est rhéofluidifiant, ce qui correspond au comportement rhéologique le plus fréquemment rencontré, et est décrit par le modèle de Carreau. A proximité du seuil, on considère deux situations : la première correspond au cas où les plaques ont une conductivité finie, la seconde à celui d'un fluide thermodépendant. Dans chaque cas, l'influence du caractère rhéofluidifiant sur la nature du motif émergeant à la bifurcation primaire et sur les instabilités secondaires est mise en évidence.

  . Leur démarche sera d'ailleurs reprise par la suite par plusieurs auteurs. Les équations aux perturbations en situation bidimensionnelle sont écrites en termes de fonction de courant et vorticité. Elles sont résolues par une méthode aux différences finies dans un domaine II est le deuxième invariant du déviateur du tenseur des contraintes. Globalement, l'approche utilisée consistait à déterminer une solution stationnaire à un nombre de Rayleigh élevé, ensuite faire décroître le nombre de Rayleigh jusqu'à atteindre le régime conductif. Les auteurs montrent que la valeur critique du nombre de Rayleigh à partir de laquelle la convection démarre diminue lorsque l'indice de structure décroit. La notion de bifurcation sous critique n'est pas évoquée. Il est aussi montré que le coefficient de transfert de chaleur augmente avec le caractère rhéofluidifiant.

	carré. Deux modèles rhéologiques ont été considérés : modèle en loi puissance et modèle d'Ellis 1 η = 1 η 0 1 + ne-1 √ τ II τ 12 , où τ

Une corrélation reliant le nombre de Nusselt à l'indice de structure n p est proposée

[33]

  sur le devenir des motifs de convection loin des conditions critiques. Les deux premiers chapitres de la thèse répondent au premier point tandis que le troisième chapitre est consacré au second.Dans le premier chapitre de la thèse, on considère des plaques de conductivité finie où les motifs de convection sont soit des rouleaux soit des carrés. Une analyse multi-échelles est menée afin de mettre en évidence l'influence de la rhéologie sur les différentes instabilités secondaires limitant la gamme de nombres d'onde stables des différentes structures.

	Chapitre 1
	Dans le deuxième chapitre de la thèse, on s'intéresse à un fluide thermodépendant et rhéofluidifiant pour lequel la convection s'organise en rouleaux ou en hexagones. A nouveau, une approche faiblement non linéaire est utilisée pour étudier l'impact du Instabilities of convection patterns
	comportement rhéofluidifiant sur la stabilité des motifs vis à vis d'une modulation spatiale à proximité du seuil. Dans le troisième et dernier chapitre, une approche fortement non for a shear-thinning fluid between
	slabs of finite conductivities
	[5] A. Pellew and Southwell R.V. On maintained convective motion in a fluid heated
	from below. Proc. R. Soc. Lond. 176, pages 312-343, 1940.
	[6] RJ Goldstein and DJ Graham. Stability of a horizontal fluid layer with zero shear
	boundaries. The Physics of Fluids, 12(6) :1133-1137, 1969.

linéaire reposant sur une procédure de continuation est mise en oeuvre. On considère un fluide purement rhéofluidifiant placé entre deux plaques parfaitement conductrices.

Dans cette configuration, les rouleaux de convection sont les motifs rencontrés. Il s'agit de déterminer de proche en proche les grandeurs caractéristiques de l'écoulement comme les champs de vitesse, de température ou encore le nombre de Nusselt au fur et à mesure que la convection gagne en intensité.

Le premier chapitre a fait l'objet d'un article intitulé ≪ Instabilities of convection patterns in a shear-thinning fluid between plates of finite conductivity ≫, paru dans Physical Review E en Novembre 2017 : il est ainsi écrit en anglais. Le deuxième chapitre qui constitue un projet d'article l'est également. Bibliographie [1] H. Bénard. Les tourbillons cellulaires dans une nappe liquide. Rev. Gén. Sci. pures et appl. 11, pages 1261-1271 & 1309-1328, 1900.

[2] Lord Rayleigh. On the convective currents in a horizontal layer of fluid when the higher temperature is on the under side. Phil. Mag. 32, pages 529-546, 1916. [3] JRA Pearson. On convection cells induced by surface tension. Journal of fluid mechanics, 4(5) :489-500, 1958. [4] M Rieutord. Fluid Dynamics : An Introduction. Springer, 2015. [7] DTJ Hurle, Eric Jakeman, and Edward Roy Pike. On the solution of the bénard problem with boundaries of finite conductivity. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, 296(1447) :469-475, 1967.

Sommaire 1.1 Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 1.3 Problem formulation . . . . . . . . . . . . . . . . . . . . . . . . . 31 1.4 Linear stability analysis . . . . . . . . . . . . . . . . . . . . . . . 34 1.4.1 Critical conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 34 1.4.2 Characteristic time of the instability . . . . . . . . . . . . . . . . 37 1.5 Pattern selection at the onset of convection . . . . . . . . . . . 37 1.5.1 Principles of the amplitude expansion method : case of square . . 38 1.5.2 Nature of the primary bifurcation . . . . . . . . . . . . . . . . . 39 1.5.3 Convective patterns at threshold . . . . . . . . . . . . . . . . . . 39 1.5.4 Flow structure, viscosity and temperature fields . . . . . . . . . . 41 1.5.5 Comparison between roll and square solutions for fixed χ and λ . 43 1.5.6 Heat transfer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

1.6 Secondary instabilities . . . . . . . . . . . . . . . . . . . . . . . . 46 1.6.1 Stability of a roll modulated solution . . . . . . . . . . . . . . . . 46 1.6.2 Stability of a square modulated solution . . . . . . . . . . . . . . 51 1.7 Numerical solutions of amplitude equations . . . . . . . . . . . 54 1.7.1 Numerical simulation . . . . . . . . . . . . . . . . . . . . . . . . 54 1.7.2 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 1.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 1.A Eigenvalues versus Q 1 and Q 2 . . . . . . . . . . . . . . . . . . . . 64 1.B Stability of the ETD2 scheme . . . . . . . . . . . . . . . . . . . . 64 1.B.1 Case of a single ODE . . . . . . . . . . . . . . . . . . . . . . . . 64 1.B.2 Extension to two coupled Ginzburg-Landau equations . . . . . . 66 28 CHAPITRE 1. Bibliographie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 1.1 Abstract

  1 •r + B(t) e ik 2 •r ψ 13 (z) + c.c. + ...

	In Eq. (1.40), ψ (r, z, t) stands for the vertical velocity perturbation, w, or the temperature
	perturbation θ or θ p . For the vertical velocity perturbation, ψ ij is denoted F ij , and for the
	temperature perturbation, ψ ij is denoted G ij .
	In the square lattice, time evolution of the amplitude perturbations is governed by Stuart-
	Landau amplitude equations,

  Solution proportional to |A p | 2 (zero mode) . . . . . . . . . . . . 102 2.B.2 Solution proportional to A 2 p exp(2ik p • r) . . . . . . . . . . . . . . 103 2.B.3 Solution proportional to A p A * q E p E * q . . . . . . . . . . . . . . . . 103 2.B.4 Solution proportional to exp (ik p • r) . . . . . . . . . . . . . . . 104 Bibliographie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

	74	CHAPITRE 2.
	2.1 Abstract
	Bibliographie

2.B.1

  where T 1 (z) stasfiesD 2 T 1 = 2 [G 11 (DF 11 ) + F 11 (DG 11 )) ,

		(2.95)
	with
	T 1 = 0 at z = 0 and z = 1.	(2.96)
	As for the linear problem, equation (2.95) with the boundary conditions (2.96) is solved
	numerically using a spectral Chebyshev collocation method.

  Carreau, on ne peut pas calculer directement les termes non linéaires visqueux dans l'espace spectral de

	Fourier-Chebyshev comme on l'a fait pour les termes linéaires et les termes non linéaires
	d'inertie. Pour obtenir la projection de ces termes, on est amené à utiliser une méthode
	pseudospectrale : elle consiste à évaluer d'abord les termes visqueux dans l'espace phy-
	sique avant de les projeter dans l'espace spectral de Fourier-Chebyshev. La procédure est
	détaillée dans les paragraphes suivants.
	Evaluation des termes visqueux dans l'espace physique
	L'espace physique Ω considéré est une grille bidimensionnelle construite sur une longueur d'onde L c = 2π k c suivant x et délimitée par les deux parois en z = -1 2 et z = 1 2

  Dans cette section, on explique comment une fonction F (x, z) quelconque évaluée aux points (3.27) est projetée dans l'espace de Fourier Chebyshev. Concrètement, cela revient, connaissant F (x, z), à déterminer les coefficients f nm tels que :

	3.3. M ÉTHODE PSEUDO-SPECTRALE	115
	Projection dans l'espace spectral de Fourier Chebyshev
		28)
	Toutefois dans la pratique, il est souvent nécessaire de considérer des nombres de
	points bien supérieurs aux valeurs minimales fournies par la règle d'Orszag. C'est
	particulièrement le cas pour le nombre de points N d . Cela s'explique par le fait que les
	termes non linéaires visqueux NV varient fortement sur de faibles distances. Si on ne
	considère pas suffisamment de points, on ne parvient pas à capter ces variations ce qui
	conduit à une estimation trop grossière et finalement à une mauvaise projection dans
	l'espace spectral de Fourier-Chebyshev.

Table 3 .

 3 'on se réfère au tableau 3.1, 6 modes de Fourier sont suffisants pour garantir la convergence de la méthode. Toutefois, en Newtonien, tous les calculs peuvent être menés en restant dans l'espace de Fourier-Chebyshev ce qui rend la procédure de continuation rapide. Pour cette raison, on a choisi d'utiliser 12 modes de Fourier afin d'avoir les résultats les plus précis possibles. Les évolutions des nombres de Nusselt avec ǫ, issus des différentes simulations, sont présentées sur la Figure 3.13. Cette dernière montre que les courbes Nu = f (ǫ), les échanges de chaleur soient indépendants du nombre de Prandtl au delà d'une valeur suffisamment grande de ce dernier : c'est effectivement ce que montre la Figure 3.13 où les courbes Nu = f (ǫ) sont identiques pour Pr ≥ 10. 2 -Comparaison des nombres de Nusselt obtenus par différents auteurs dans la limite des grands nombres de Prandtl

	"	(•/' 1).:1 (• . " Ofi .:::t • : 4	(1.4 0.2 • .. ッ N セ @ ....	" C2 セ I N R @

セ I @ .. l

Table 3 .

 3 3 -Comparaison des nombres de Nusselt obtenus par différents auteurs pour Pr = 1.

  et de A. Chorin[12] pour Pr = 1. Les valeurs des nombres de Nusselt obtenues par ces auteurs sont comparées à celles provenant de la procédure de continuation dans le tableau 3.3 pour différents ǫ. Là encore, les écarts entre nos résultats et ceux de ces deux auteurs sont minimes. Toutefois, on peut signaler que nos résultats se rapprochent davantage de ceux de W. Plows qui a repris la démarche de A. Chorin en considérant des maillages plus fins ce qui est censé conduire à des simulations plus précises. le premier est le modèle en loi puissance, aussi appelé modèle de Ostwald-de Waele, et le second est le modèle d'Ellis (cf[17]). Comme le soulignent ces deux auteurs, le modèle en loi puissance a l'avantage d'être simple mais conduit à une viscosité infinie pour un taux de cisaillement nul. C'est cette raison qui les a amenés à considérer également le modèle d'Ellis qui ne présente pas cet inconvénient. Quelque soit le modèle adopté, H. Ozoe et S.T. Churchill ont montré que les échanges de chaleur sont d'autant plus importants que le comportement rhéofluidifiant est marqué. viscosité de Parmentier ≫ , notée µ P , définie comme la moyenne de la viscosité pondérée par Γ dans la couche de fluide Ω : Les calculs de M. Bouteraa et al. effectués à la fois dans le cas de parois rigides et libres indiquent également que la transition d'une bifurcation surcritique à sous-critique survient quand le caractère rhéofluidifiant est suffisamment marqué. Plus précisément, cette transition ne dépend que d'un seul paramètre rhéologique α appelé ≪ degré de rhéofluidification ≫ et défini par :

	3.7. R ÉSULTATS POUR DES FLUIDES RH ÉOFLUIDIFIANTS	CHAPITRE 3. 137
	soit Ra sc < Ra c : la bifurcation de l'état conductif à l'état convectif devient donc
	sous-critique.
	En 2011, B. Albaalbaki et R.E. Khayat se sont plus spécifiquement intéressés à la stabilité
	des motifs. Dans [20], ils montrent que :
	-les rouleaux de convection ne sont stables à proximité du seuil que pour une faible
	rhéofluidification,
	-les hexagones et les carrés sont stables quand le caractère rhéofluidifiant est très
	Ce résultat a été confirmé en 1978 par E.M. Parmentier dans [16]. Cet auteur s'est notam-prononcé.
	ment interrogé sur la possibilité d'établir une corrélation entre les échanges de chaleur et Une étude menée en 2015 par M. Bouteraa et al. [1] vient confirmer, nuancer et infir-
	les paramètres rhéologiques m (consistance) et n (indice de structure) d'un fluide décrit mer certains des résultats établis par [20]. Elle vient notamment conforter le fait que
	par une loi puissance : les conditions critiques restent effectivement les mêmes qu'en Newtonien conformément µ = mΓ n-1 2 , (3.68) à [19] et [20].
	où Γ désigne le deuxième invariant du tenseur des taux de déformations.
	En se basant sur des considérations énergétiques, il a introduit une grandeur appelée
	µ P = Ω µΓdxdz Ω Γdxdz	(3.69)
	et il a construit un ≪ nombre de Rayleigh de Parmentier Ra p ≫ :
	Ra p =	Ra µ P	.	(3.70)
	E.M. Parmentier a ensuite montré que Nu = f (Ra p ) était une courbe maîtresse permet-
	tant de décrire les évolutions des transferts de chaleur indépendamment des paramètres
	m et n de la loi puissance (3.68).
	D'autres études, plus récentes, se sont quant à elles intéressées à l'influence d'un effet
	rhéofluidifiant sur les motifs de convection et leurs stabilités à proximité du seuil de la
	convection. Dans ces travaux, le comportement rhéofluidifiant est décrit par le modèle de
	Carreau (3.9) introduit dans [18] qui permet d'éviter les singularités à taux de cisaille-
	3.7 Résultats pour des fluides rhéofluidifiants ment nul et qui a l'avantage de reposer sur des considérations physiques et théoriques.
	Les différentes structures convectives sont étudiées grâce à des approches faiblement non
	3.7.1 Les principaux travaux existants en non-Newtoniens linéaire similaires à celles menées dans les chapitres précédents. Le travail le plus ancien
	a été réalisé par N.J. Balmforth et A.C. Rust en 2008 dans [19] en considérant des parois A notre connaissance, les premières études numériques concernant la convection de Rayleigh-Bénard pour des fluides non-Newtoniens remontent aux années 1970 avec les libres. Ces auteurs sont parvenus à établir un certain nombre de résultats :
	travaux de H. Ozoe, S.T. Churchill [15] et de E.M Parmentier [16]. Ces travaux, reposant -les conditions critiques (k c , Ra c ) à partir desquelles se développe la convection sont
	sur la méthode des différences finies, se sont focalisés sur l'influence d'un comportement indépendantes du comportement rhéofluidifiant décrit par le modèle de Carreau,
	rhéofluidifiant sur les évolutions des transferts de chaleur quand la convection s'intensifie. -pour un effet rhéofluidifiant suffisamment fort, la convection débute à une valeur du
	H. Ozoe et S.T. Churchill ont considéré deux modèles rhéologiques dans leur approche : nombre de Rayleigh Ra sc inférieure à celle prédite par l'analyse linéaire de stabilité

≪

  .72) CHAPITRE 3. qui est une généralisation du modèle de Carreau (3.9). Les caractéristiques de l'écoulement ont été calculées à l'aide de la méthode des différences finies en considérant différents types de conditions limites (température imposée, flux imposé, parois rigides ou libres) et différents confinements. Quelle que soit la configuration envisagée, O. Benouared, M. des éléments finis. Ces auteurs ont travaillé dans la limite λ → ∞ où le modèle de Carreau(3.9) approche la loi puissance (3.68) avec m = λ n-1 . Ils ont confirmé que, dans cette limite, la courbe Nu = f (Ra p ) obtenue par E.M. Parmentier pour une loi puissance était bien une courbe maîtresse permettant de caractériser les transferts de chaleur. Le principal apport de leur travail réside dans l'établissement de corrélations permettant de déterminer le nombre de Rayleigh sous-critique Ra sc ainsi que l'évolution du nombre de Nusselt Nu en fonction des paramètres rhéologiques λ et n c du modèle de Carreau.

	Mahmoud et N.A. Messaoudene montrent que le nombre de Rayleigh sous-critique Ra sc	à
	partir duquel démarre la convection diminue lorsque les effets rhéofluidifiants s'intensifient
	(quand n c diminue ou λ augmente dans la loi (3.72)).
	Une seconde étude [22] réalisée par M. Jenny, E. Plaut et A. Briard en 2015 vient confirmer
	ce résultat. Les équations (3.2)-(3.12) y sont simulées à l'aide du logiciel Freefem++, basé
	sur la méthode

Table 3 .

 3 4 -Comparaison des nombres de Nusselt obtenus pour Pr = 10 avec différents nombres N de modes de Fourier pour n c = 0.5 et λ = 0.02. A ǫ donné, les valeurs en bleu correspondent à des écarts relatifs avec la valeur obtenue pour N ≥ 11 inférieurs à 1%. Celles en rouge correspondent à des écarts supérieurs à 1%.

	11

  les mêmes valeurs de n c et λ, en l'occurence n c = 0.5 et λ = 0.02, notre courbe Nu = f (ǫ) est identique à celle obtenue par Jenny et al. dans[22] dans leur gamme de ǫ. Cette influence de n c sur les échanges de chaleur se justifie physiquement. Quand n c diminue, les effets rhéofluidifiants deviennent plus importants et la viscosité globale est réduite : la

	0.6 t	,	'	• .... t O.E	•' J.: pot .u c 3 •'a::::	•) 1).5"	.. ,•, p:>UT;-l :i	'	•) tJ.5
			2	• . . . . c	2	•)		2	•)
					'		'

  'objectif de cette thèse était d'étudier l'influence d'un effet rhéofluidifiant, qui correspond à une diminution non linéaire de la viscosité quand le cisaillement augmente, sur les motifs de convection se développant dans la convection de Rayleigh Bénard. Pour ce faire, on a commencé par considérer la convection à proximité des conditions critiques. Deux configurations ont été étudiées : la première fait intervenir des plaques de conductivité finie et la seconde fait intervenir un fluide thermodépendant. Pour de faibles valeurs de ce paramètre, c'est l'instabilité rouleaux croisés qui est dominante comme en Newtonien. En revanche pour des valeurs élevées de α, c'est l'instabilité d'Eckhaus qui l'emporte. A noter qu'on constate dans les deux cas la présence de l'instabilité Zigzag qui ne dépend pas des paramètres rhéologiques à proximité du seuil et rend instables les nombres d'ondes k tels que k ≤ k c . Pour les carrés, c'est l'instabilité rec-

	Carreau, les effets rhéofluidifiants sont plus importants ce qui réduit les frottements
	visqueux et renforce en conséquence la convection et les échanges de chaleur induits
	par cette dernière.
	-quand le caractère rhéofluidifiant est suffisamment prononcé, la bifurcation de l'état
	conductif à l'état convectif devient sous critique.
	été considéré.
	Enfin, on pourrait étendre la procédure de continuation à d'autres motifs de convection
	comme les hexagones qui se développent quand on considère des fluides thermodépendants.
	Pour de tels fluides, la théorie faiblement non linéaire classique prévoit l'existence d'hexa-
	gones uniquement à proximité du seuil de la convection, ces derniers laissant place à des
	rouleaux quand la convection gagne en intensité. Toutefois, des travaux expérimentaux tangulaire qui est dominante et qui restreint d'autant plus la gamme de nombres d'ondes
	[23], [24] rapportent la réapparition d'hexagones pour des écarts relatifs au seuil suffi-stables que le caractère rhéofluidifiant est prononcé.
	samment importants : on parle d'≪ hexagones ré-entrants ≫. La seule étude théorique
	justifiant ce phénomène a été menée par Madruga et al. dans [25] : ces auteurs ont montré

Nos travaux se distinguent de ceux de Benouared et al.

[21]

et de Jenny et al.

[22]

parce qu'ils font intervenir des valeurs de λ plus faibles et plus couramment rencontrées et que la gamme d'écarts relatifs au seuil explorée est beaucoup plus large.

Plusieurs perspectives sont possibles à la suite de ce travail. La première serait de déterminer les branches instables de la bifurcation sous-critique observée dans le cas de forts effets rhéofluidifiants. La méthode de continuation permettrait alors d'obtenir l'intégralité des branches. La deuxième serait d'étudier la stabilité des rouleaux de convection qui se développent vis à vis de perturbations tridimensionnelles. C'est d'ailleurs l'objectif principal de cette étude. En Newtonien, ce travail a déjà été réalisé dans

[13]

: un diagramme appelé ≪ ballon de Busse ≫, en hommage à son créateur, représentant le domaine de stabilité des rouleaux dans le plan (k, Ra) a ainsi été établi. Nous souhaitons examiner l'influence du caractère rhéofluidifiant sur la gamme des nombres d'ondes stables et les mécanismes d'instabilité qui délimitent les frontières du domaine de stabilité. A notre connaissance, l'impact de la rhéologie sur la stabilité des rouleaux loin du seuil n'a jamais qu'une analyse faiblement non linéaire permettait de mettre en évidence des ≪ hexagones ré-entrants ≫ à condition de considérer un fluide dont les propriétés thermophysiques varient explicitement avec l'écart relatif au seuil. Cette approche apporte une justification

Conclusion et perspectives

LDans le cas de plaques de conductivité finie, on montre que les motifs de convection émergeant sont soit des rouleaux soit des carrés en fonction de l'épaisseur des plaques, de leur conductivité et du degré de rhéofluidification α. On montre qu'une faible conductivité des plaques favorise l'apparition de structures carrées. En effet, plus la conductivité est faible et moins la convection est forte ce qui conduit à un couplage entre les modes décrivant les carrés suffisamment faible pour que ces derniers soient stables. De manière semblable, une grande épaisseur de plaques favorise les carrés au détriment des rouleaux. A l'inverse, un effet rhéofluidifiant prononcé conduit à une convection plus marquée et donc à des couplages non linéaires supérieurs entre les différents modes caractérisant les carrés. Pour un effet rhéofluidifiant suffisamment marqué, ces couplages deviennent assez forts pour déstabiliser les carrés : le caractère rhéofluidifiant favorise donc les rouleaux de convection. L'influence du caractère rhéofluidifiant sur les instabilités secondaires a également été étudié. Pour les rouleaux, il a été démontré que l'instabilité secondaire qui restreignait le plus la gamme de nombres d'ondes stables dépendait du degré de rhéofluidification α.

Aux ordres δ p où p ≥ 2, les systèmes obtenus admettent des seconds membres comportant, entre autres, les grandeurs Ra p-1 qui sont des inconnues du problème et des termes dits résonants. Ces derniers conduisent lors de la résolution du système d'équations à des termes ≪ séculaires ≫ à savoir des termes croissant linéairement avec le temps. Ces termes séculaires ne pouvant exister physiquement, les coefficients Ra p-1 doivent donc être déterminés de sorte à annuler les différents termes résonants. Une condition de solvabilité doit donc être appliquée afin de résoudre les systèmes d'équations différentielles ordinaires obtenus aux différents ordres δ p , p ≥ 2. La plus couramment utilisée est l'alternative de Fredholm : elle stipule que ≪ toute équation de la forme Lu(t) = f (t) n'a de solution que si la fonction f (t) du membre de droite est orthogonale au noyau de l'opérateur adjoint de L ≫( cf[11] ).En appliquant la démarche précédente, Malkus et Veronis sont parvenus à déterminer de proche en proche les inconnues v p , θ p pour des motifs de convection en rouleaux, en carrés et en hexagones dans le cas de parois libres et dans le cadre des approximations de Boussinesq. Une fois ces différentes solutions stationnaires déterminées, ces auteurs se sont interrogés sur le critère permettant de prédire l'émergence de tel ou tel motif. Leurs calculs suggèrent que la structure qui se développe est celle donnant lieu au plus grand gradient thermique quadratique moyen. Ainsi, selon ces auteurs, les carrés qui maximisent cette quantité, seraient les motifs les plus probables.En 1965, Schlüter, Lortz et Busse reprennent la démarche adoptée dans[10] et l'étendent au cas de parois rigides. Leurs calculs viennent infirmer les résultats de Malkus et Veronis concernant la stabilité des différents motifs de convection. En effet, Schlüter et al montrent dans[12] que les rouleaux de convection sont les seules structures stables

CHAPITRE 2.

The reference viscosity, μ0 , is the zero-shear rate viscosity evaluated at T0 , i.e. the mean of the boundary temperatures.

Rheological model and parameters

The fluid is assumed to be purely viscous and shear-thinning. The viscous stress-tensor writes τ = µ (Γ) γ (2.6) with γ = ∇u + (∇u) T the rate-of-strain tensor and Γ = 1 2 γij γij its second invariant. We assume a Carreau-law fluid where the viscosity depends exponentially on temperature,

with µ 0 = μ0 / μ0 and µ ∞ = μ∞ / μ0 the viscosities at low and high shear rate, b the thermodependency coefficient which measures the sensitivity of viscosity to variation in temperature, n c < 1 the shear-thinning index, λ the characteristic time of the fluid. The characteristic shear rate for the onset of shear-thinning is determined by 1/ λ. The infinite shear viscosity, μ∞ , is generally significantly smaller (10 3 to 10 4 times smaller) than μ0 , [16,17]. The ratio μ∞ /μ 0 will be thus neglected in the following. The exponential model used for the viscosity thermodependency is referred in the literature as Frank-Kamenetski model and can be derived from the Arrhenius law by expanding the arguments of the exponential (in the Arrhenius law) in a Taylor series about the reference temperature T0 [7].

The dimensionless effective viscosity is then

where, µ b (z) = exp (c (z -1/2)) is the viscosity profile at quiescent state and c = b δ T is a measure of the viscosity contrast. The Newtonian behavior, μ = μ0 , is obtained by setting n c = 1 or λ = 0. The viscosity ratio across the fluid layer,

depends on b and δ T , but not on the temperature level. For a small amplitude disturbance, the viscosity can be expanded about the hydrostatic solution,

-la viscosité moyenne µ M aussi utilisée dans [22] et définie par :

où Ω désigne le domaine occupé par la couche de fluide.

La Figure 3.24 montre les évolutions du nombre de Nusselt Nu, quantifiant les échanges de chaleur, avec l'écart relatif au seuil ǫ pour différentes valeurs de n c .

Abstract

In the Rayleigh-Bénard configuration, we consider a thin layer of fluid confined between two horizontal slabs which is heated from below and cooled from above. This layer undergoes an instability if the thermal gradient is strong enough : a transition from the conductive state to the convective state and called ≪ primary bifurcation ≫occurs. Moreover, it happens in an ordered way : we can notice the emergence of various convection patterns such as rolls, squares or hexagons. In their turn, these patterns undergo ≪ secondary instabilities ≫ that limit the range of stable wavenumbers. These instabilities are studied theoretically firstly near the threshold thanks to a weakly nonlinear approach, and secondly far from critical conditions thanks to a strongly nonlinear approach. We consider a shear thinning fluid, the most common rheological behavior, which is described by the Carreau model. Near the threshold, two situations are considered : the first corresponds to finite conductivity plates, the second corresponds to a thermodependent fluid. In each case, the influence of the shear thinning effect on the nature of the pattern emerging at the primary bifurcation and on secondary instabilities is highlighted.

To study the convection patterns far from the critical conditions, a continuation procedure is used to determine, step by step, the characteristics of the flow, such as the velocity or temperature fields and the Nusselt number.

Keywords : Convection, Rayleigh-Bénard, Rheology, Bifurcation, Secondary instabilities