Géodynamique S. Labrosse 27 septembre 2010 2 Introduction Depuis Sadi Carnot qui, dans la première page de son livre (Carnot, 1824) souligne l’importance de la chaleur comme origine de la dynamique de la Terre, on a coutume de dire que la Terre est une machine thermique. Si l’analogie entre le fonctionnement de la Terre et celui des machines à vapeur, dont la compréhension était le principal objectif du travail de Carnot et a été en grande partie à l’origine du développement de la thermodynamique, n’est pas si évident que cela à y regarder de plus près, il est indéniable que l’évolution thermique de la Terre est à l’origine de toute sa dynamique interne. Ou pour être plus précis car une telle dynamique a également pu fonctionner durant des périodes de réchauffement si jamais elles ont existé, c’est le transfert de chaleur qui est a l’origine de la dynamique interne de la Terre. Le lien entre le transfert de chaleur et la dynamique de la Terre se crée parce qu’un des modes de transfert est la convection, c’est à dire le transport par mouvement de matière. Dans la Terre, ce mode se manifeste par la tectonique des plaques qui est à l’origine de la formation des chaı̂nes de montagnes, du volcanisme et des seismes. Cette convection est aussi à l’origine du champ magnétique de la Terre, par le biais de la dynamo à l’œuvre dans le noyau terrestre. Cependant, la convection ne supprime pas la conduction, qui reste le mode de transfert de chaleur principal lorsque le mouvement vertical de la matière est empêché, c’est à dire lorsque l’on s’approches des limites horizontales, soit rigides (dans une expérience de laboratoire) soit libre et imposées par des différences de densité chimiques trop importantes, comme c’est le cas à la frontière noyau-manteau (FNM) et à la surface de la Terre. Le transfert conductif se fait par transmission de l’agitation moléculaire de proche en proche et il s’agit donc d’un processus microscopique, que l’on ne considérera ici que dans ses effets macroscopiques quantifiés par la loi de Fourier et la conductivité thermique. Enfin, un troisième mode de transfert de chaleur existe, le transfert radiatif, et il peut s’avérer important dans la Terre. La chaleur dans ce mode de transfert est transportée par rayonnement électromagnétique, c’est à dire par de la lumière. C’est un mode de transfert qui intéresse particulièrement la dynamique de l’atmosphère mais beaucoup moins la Terre interne, plutôt opaque. Cependant, à haute température, le rayonnement de la matière peut devenir important comme toute personne appréciant le chauffage par un feu de bois le sait. L’intérieur de la Terre étant chaud et les cristaux qui la composent comme l’olivine non-totalement opaque, la possibilité d’un tranport de chaleur radiatif dans la Terre n’est pas totalement exclue. Un tel transport peut alors être pris en compte par une conductivité effective, un point qui sera discuté en fin de cours. Pour des raisons pédagogiques, nous commencerons par traiter quelques problèmes simples de diffusion de la chaleur. Les généralités sur le transfert de chaleur, basées sur les principes de 4 la thermodynamique, seront ensuite présentés. Des méthodes plus générales de résolution des problèmes de diffusion pourront alors être exposées. Nous aborderons alors le problème de la convection de Rayleigh-Bénard et pourront ainsi comparer ses prédictions aux observations terrestre. Les modèles d’évolution thermique de la Terre pourront alors être exposés et discutés. Chapitre 1 Loi de Fourier et premières applications 1.1 Introduction Dans un premier temps, pour introduire les approches utilisées, le cas de la conduction de la chaleur sera discuté, la loi de Fourier introduite. Des applications à quelques problèmes stationnaires permettront de fixer les idées. 1.2 Loi de Fourier La loi de Fourier est une relation empirique entre la différence de température et le flux de chaleur. Il ne s’agit ici que de flux diffusif dans lequel la chaleur est transportée en transmettant de proche en proche l’agitation moléculaire. Selon la théorie cinétique des gaz, la température a pour origine microscopique l’énergie cinétique moyenne des molécules, avec la relation 3kT /2 = m < v 2 > /2, k étant la constante de Boltzmann. Les chocs entres molécules permettent alors de tranmettre l’énergie cinétique et donc la chaleur. Dans le cas des solides, ce sont les vibrations des molécules autour de leur position d’équilibre qui se transmettent via les liaisons inter-atomiques. Des ondes de vibrations se propagent et on parle de phonons. Nous n’irons pas plus loin dans la discussion de la base microscopique de la diffusion thermique et nous nous contenterons d’introduire sa quantification macroscopique via le coefficient de conductivité thermique, k. Considérons un mur, d’épaisseur d, initialement maintenu à une température uniforme T0 , qui sépare une pièce de l’extérieur, les deux étant à la même température T0 . Au temps t = 0, on chauffe la pièce et on la maintient à une température T1 . Comment est la distribution de température dans le mur ? Clairement, celle-ci ne dépend que de la distance aux bords du mur, notée par exemple x, et pas de la position le long du mur. Au cours du temps, la température dans le mur change, celle sur ses faces étant maintenue constante (fig. 1.1). Après un temps très long (infini en fait), la température dans le mur ne change plus et varie linéairement en fonction 6 Loi de Fourier et premières applications 1 t=0,002 t=0.01 t=infini Température 0.8 Fig. 1.1 – Distribution de température dans un mur, initialement à T0 = 0 et dont une face est brusquement portée à T1 = 1. Au cours du temps, la température passe de la courbe rouge à la verte et atteint finalement la courbe bleue après un temps infini. 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 Distance de x. Si l’on arrête de chauffer la pièce, la température va baisser et l’expérience montre que pour maintenir la température dans la pièce avec le profil de température linéaire atteint après un temps long, il faut fournir un flux de chaleur égal à T1 − T0 Q =k , A d (1.1) A étant la surface du mur. Ceci peut se mesurer facilement en regardant sa facture d’électricité ou de gaz. Le paramètre k est appelé conductivité thermique. La relation linéaire entre le flux de chaleur et la différence de température est valide tant que cette différence n’est pas trop grande. Le flux de chaleur Q a pour unité le Watt, W. Une telle expérience peut être faı̂te pour tout type de matériaux, qu’il soit solide, liquide ou gazeux, tant qu’aucun mouvement convectif ne se déclenche et que la radiation est négligeable. Cependant, la valeur de la conductivité thermique dépend du matériau considéré et également de paramètres physiques tels que la température, la pression. Le temps pour atteindre la distribution linéaire de température dans le mur dépend de sa conductivité thermique, de sa capacité calorifique et de son épaisseur. Plus le mur est fin, plus vite cette distribution limite sera établie. À la limite d’un mur infiniment fin, le profil linéaire de température est établit instantanément et on voit que le terme de droite de l’équation (1.1) fait apparaı̂tre la dérivée de la température en fonction de la direction x. Avant d’écrire explicitement cette relation différentielle, il est nécessaire d’introduire le caractère vectoriel du flux de chaleur. Dans l’exemple représenté sur la figure 1.1, la température élevée est du coté des x faibles et le flux de chaleur va vers l’extérieur, c’est à dire vers les x élevés. On voit donc que le vecteur flux de chaleur a une composante sur l’axe des x (uniquement celle-là dans ce cas), qx , qui est positive, c’est à dire que la chaleur va de l’intérieur vers l’extérieur, alors que la dérivée de la température par rapport à x est négative. On a donc la relation suivante qx = −k dT , dx (1.2) qx étant une densité de flux de chaleur, c’est à dire un flux de chaleur par unité de surface, avec pour unité le W m−2 . L’équation (1.2) est la forme unidimensionnelle de la loi de Fourier, valide uniquement dans le cas où la température ne dépend que de x. Dans le cas général, le flux de chaleur a trois 1.3 Diffusion de chaleur en état stationnaire 7 composantes et s’écrit pour un matériau isotrope q = −k∇T. (1.3) Un matériau est dı̂t isotrope si son comportement, en particulier sa conductivité thermique, ne dépend pas de la direction. Les matériaux cristallins ne sont généralement pas isotropes et le manteau terrestre présente par exemple une anisotropie des vitesses de propagations sismiques. il est probable qu’une anisotropie de conductivité soit également présente. Dans ce cas là, la loi de Fourier doit faire intervenir un tenseur d’ordre 2 des conductivités au lieu du paramètre scalaire introduit précédemment. La discussion de ce type de matériau dépasse l’objet de ce cours. Le fait que la conductivité thermique est un paramètre positif est apparu naturellement dans la discussion précédente et le signe négatif de la loi de Fourier (1.2) est une conséquence du fait que la chaleur va du coté chaud vers le coté froid. On verra plus loin que l’on peut relier cette observation au second principe de la thermodynamique. 1.3 Diffusion de chaleur en état stationnaire Le cas d’un mur solide a déjà été abordé (§ 1.2) et nous allons maintenant nous intéresser à un cas plus pertinent pour les sciences de la Terre : la distribution de température dans une sphère chauffée par radioactivité. La résolution du problème de la conduction en état stationnaire pour un problème ayant une symétrie élevée est généralement simple. Il s’agit d’écrire un bilan de chaleur en état stationnaire qui s’écrit sous la forme suivante : chaleur sortante = chaleur entrante + chaleur produite. Ce bilan peut se faire pour toute partie du système considéré, tant que celui-ci est en état stationnaire. On considère donc une sphère pleine dont la surface est maintenue à une température Ts avec un taux de production de chaleur uniforme et constant dans le temps h en Wm−3 . La symétrie sphérique du problème est telle que la température ne dépend que de la distance au centre et le flux de chaleur est dans la direction de l’extérieur. On écrit le bilan de chaleur d’une coquille sphérique comprise entre les sphères de rayons r et r + dr : le flux de chaleur qui sort en r + dr est égal à la somme du flux qui entre en r et de la chaleur produite entre r et r + dr : q(r + dr)4π(r + dr)2 = q(r)4πr2 + h En prenant la limite dr → 0, on trouve d 2
r q = hr2 dr soit  4π  (r + dr)3 − r3 . 3 1 d 2
r q = h. r2 dr (1.4) (1.5) Le terme de gauche de la seconde équation n’est autre que la divergence en coordonnées sphériques d’une fonction q qui ne dépend que du rayon. Cette équation a une solution de la forme q(r) = h C1 r+ 2 3 r (1.6) 8 Loi de Fourier et premières applications et la constante d’intégration est déterminée par la condition au centre : le flux de chaleur doit rester borné et doit même être nul pour satisfaire la symétrie du problème. Donc C1 = 0. Avant d’aller plus loin et de déterminer la température, on peut calculer le flux de chaleur total qui sort à la surface de la sphère, r = R : Qs = 4πR2 q(r) = 4π 3 R h, 3 (1.7) ce qui n’est pas une surprise. En effet, en état stationnaire, le flux à la surface doit être égal au total de la chaleur produite dans le volume. Pour calculer la température, il suffit maintenant d’exprimer le lien entre le flux et la température, à savoir la loi de Fourier (1.3) qui prend, dans le système de coordonnées sphériques du problème, la forme suivante : q(r) = −k dT . dr (1.8) La température a donc la forme T (r) = − h 2 r + C2 6k (1.9) et la constante d’intégration est obtenue en utilisant la condition limite à la surface, T = Ts , qui donne T (r) = h 2 (R − r2 ) + Ts . 6k (1.10) La température maximale est bien sûr obtenue au centre de la sphère et est égale à Tc = Ts + h 2 R . 6k (1.11) Connaissant le flux de chaleur à la surface de la Terre, 44 TW, si l’on suppose que ce flux est maintenu par une production de chaleur uniformément répartie et que l’on est en état stationnaire, on peut calculer la production de chaleur par unité de volume (1.7) : h= 3Qs ≃ 4 10−8 W m−3 , 4πR3 (1.12) et la différence de température entre le centre de la Terre et la surface est Tc − Ts = Qs = 6.8 104 K 8πkR (1.13) où l’on a supposé une conductivité thermique égale à 4 Wm−1 K−1 . Une telle température est tout à fait irréaliste et ne permettrait pas la présence d’une graine de fer solide au centre de la Terre. Il est facile de pointer les nombreuses raisons pour lesquelles ce calcul ne s’applique pas à la Terre et les plus importantes sont l’hypothèse de transfert de chaleur conductif et celle de stationnarité. Exercise 1 Dans le cadre du modèle de température dans la Terre proposé ci-dessus, supposons maintenant qu’en dessous d’une couche supérieure, de 100 km d’épaisseur, la conductivité thermique est N fois plus forte. En dehors de cette discontinuité, on suppose chaque coquille homogène. Calculer la température au centre de la Terre dans ce modèle et le discuter. 1.4 Équation de diffusion non–stationnaire 9 Exercise 2 Considérons une autre variante du modèle de Terre stratifiée : On suppose que toute la radioactivité est maintenant concentrée dans la partie superficielle d’épaisseur δ << R. Pour le même flux de chaleur en surface, comparer la température au centre avec celle obtenue dans le cas d’une Terre homogène. 1.4 Équation de diffusion non–stationnaire Ayant obtenu une équation de diffusion en état stationnaire et en géométrie sphérique, revenons au cas cartésien qui constituait notre premier exemple (fig. 1.1) et établissons l’équation qui permet de calculer l’évolution de la température avec le temps de l’état initial à l’état stationnaire. On suppose que la température ne dépend que de la coordonnée x. Considérons une tranche de fluide d’épaisseur faible δx et de surface A quelconque et faisons le bilan d’énergie. En tout point de notre domaine d’étude, on a un flux de chaleur q(x) = −kdT /dx. Lorsque ce flux est positif, cela signifie que la chaleur va dans la direction des x croissants, le contraire (q(x) < 0) étant également possible. En x, notre tranche de matière reçoit de la chaleur à un taux Aq(x) (ou la donne si ce nombre est négatif) alors qu’en x + δx elle perd de la chaleur à un taux de Aq(x + δx). Un autre gain que l’on peut considérer (et c’est important pour la Terre interne) est la production de chaleur par radioactivité, à un taux h par unité de volume, soit hAδx pour tout le domaine. Si les pertes et les gains sont égaux, la température de la tranche reste constante. Dans le cas contraire, elle augmente ou diminue suivant le signe de la différence. Le lien entre l’apport ou la perte nette de chaleur et l’évolution de la température est donné par un paramètre physique appelé capacité calorifique, C dont l’unité est [C] = J kg−1 K−1 . Pour obtenir une énergie par unité de temps, c’est à dire un flux en Watts, il faut multiplier C par la masse du système considéré et pas le taux de changement de température. On obtient ainsi ∂T ρδxA C hδxA , = Aq(x) − Aq(x + δx) + | {z } | {z } ∂t | {z } | {z } masse radioactivité perte gain (1.14) ρ étant la masse volumique. On n’a pas considéré de perte ni de gain dans les directions y et z car on a supposé que la température ne dépendait que de x et donc le flux de chaleur ne peut se produire que dans cette direction. On voit que cette équation peut être divisée par le volume Aδx pour obtenir ρC q(x + δx) − q(x) ∂T =− + h, ∂t δx (1.15) et en faisant tendre δx vers 0, on fait apparaı̂tre la dérivée du flux dans la membre de droite : ρC ∂q ∂T =− + h. ∂t ∂x (1.16) Remarquons que cette équation est tout à fait générale et ne fait pas intervenir la Loi de Fourier. Il s’agit juste d’un bilan d’énergie. Cependant, elle n’est pas suffisante pour résoudre le problème 10 Loi de Fourier et premières applications en termes de température et il est alors nécessaire de faire intervenir la loi de Fourier (1.2) pour obtenir une équation aux dérivées partielles pour la température :   ∂T ∂ ∂T = k + h. (1.17) ρC ∂t ∂x ∂x Un cas particulier que nous rencontrerons souvent est celui où la conductivité thermique est uniforme. On peut alors la sortir de la dérivée et diviser l’équation par ρC pour faire apparaı̂tre la diffusivité thermique, κ = k/ρC. L’équation est alors   ∂T ∂ ∂T h =κ k + . (1.18) ∂t ∂x ∂x ρC L’unité de la diffusivité est [κ] = m2 s−1 et on rencontrera d’autres grandeurs physiques ayant la même unité et elles pourront toutes être appelées diffusivité (magnétique, de quantité de mouvement, etc). Enfin, remarquons qu’en l’absence de chauffage radioactif (h = 0), une solution stationnaire au problème de diffusion existe et que dans le cas d’une conductivité constante il correspond une variation linéaire de la température (dérivée seconde nulle). Il s’agit de la situation pour t → ∞ représentée en bleu sur la figure 1.1. Dans le cas général où la température dépend des trois directions de l’espace, on procède de manière similaire mais il faut considérer un petit élément de volume (δx, δy, δz) autour du point (x, y, z). On fait le bilan de la même manière mais il faut considérer les trois composantes du flux de chaleur q = (qx , qy , qz ) : ρCδxδyδz ∂T = − δyδz(qx (x + δx, y, z) − qx (x, y, z)) ∂t − δxδz(qy (x, y + δy, z) − qy (x, y, z)) − δxδy(qz (x, y, z + δz) − qz (x, y, z)) + hδxδyδz. (1.19) Une fois encore on divise cette équation par le volume δxδyδz et on fait tendre δx, δy et δz vers 0 pour obtenir ρC ∂T ∂qx ∂qy ∂qz =− − − + h = −∇·q + h. ∂t ∂x ∂y ∂z (1.20) Comme dans le cas uni–dimensionnel, cette équation exprime la conservation de l’énergie et est de ce fait totalement générale. Par contre, elle n’est pas suffisante pour résoudre le problème est une relation phénoménologique entre le flux de chaleur est nécessaire, la loi de Fourier (1.3). On obtient ainsi l’équation ρC ∂T = ∇· (k∇T ) + h. ∂t (1.21) Dans le cas d’une conductivité thermique uniforme, on écrit l’équation de diffusion sous la forme h ∂T = κ∇2 T + . ∂t ρC (1.22) 1.5 Solutions aux exercices 1.5 11 Solutions aux exercices Solution 1 L’équation (1.5) reste valable puisqu’elle ne fait pas intervenir la conductivité thermique. Les conditions aux limites étant inchangées, la solution pour le flux est également la même, c’est à dire q(r) = hr/3. La loi de Fourier (1.8) s’applique dans chaque coquille et la solution s’écrit (i pour interne, s pour surface) : h 2 r + Ci , 6kN h Ts (r) = − r2 + Cs . 6k Ti (r) = − (1.23) Les deux constantes d’intégration doivent être déterminées en utilisant les conditions limites : la température à la surface peut être prise égale à 0, sans perte de généralité, et doit être continue à l’interface entre les deux coquilles (r = Ri = R − δ). Ceci conduit à " #  2  2 1 Ri hR 1− 1− Ci = 6k R N (1.24) 2 hR . Cs = 6k La température est finalement  h
 hR2 N − Ri (N − 1) − 6kN R   T =  hR2 1 − r22 6k R r2 R2 i si r ≤ Ri , si r ≥ Ri . (1.25) La température au centre est égale à Ci et dans la limite où l’épaisseur de la couche supérieure est faible, Ri = R − δ avec δ << R, et N >> 1, elle s’écrit T (0) = hR2 2δ , 6k R (1.26) ce qui permet de voir que cette température est d’autant plus faible que la zone à faible conductivité est fine. Pour que la température au centre soit 5500 K, il suffit qu’elle soit de 250 km, ce qui est raisonnable. Le fait de supposer une conductivité thermique plus forte en profondeur qu’au centre de la Terre permet de modéliser la convection thermique comme moyen de transport dans la Terre. Perry (Perry, 1895a,b,c) avait considéré cette situation comme alternative au modèle de Kelvin (Thomson, 1862, 1899) et montré que l’âge de la Terre pouvait ainsi être grandement augmenté. Solution 2 Le flux de chaleur en surface Qs est équilibré par une production répartie sur un volume plus petit, ayant donc une concentration plus forte, h : Qs =  4π  3 R − (R − δ)3 h, 3 soit au premier ordre en δ/R h= Qs 4πR2 δ (1.27) 12 Loi de Fourier et premières applications L’équation pour le flux de chaleur est la même que celle écrite précédemment (1.5), et dans la partie centrale elle donne d 2 C1 (r q) = 0 ⇒ q = 2 = 0 ⇒ T = Tc dr r (1.28) car le flux doit être nul au centre. Indépendamment de la valeur de la conductivité thermique, la température doit être uniforme dans la partie centrale est donc égale à la température au centre. Dans la partie superficielle, l’équation est la même que dans le cas homogène et la solution a la forme donnée par l’équation (1.10). Par contre, la condition limite n’est pas d’un flux nul en r = 0 qui est en dehors du domaine mais q(R − δ) = 0 qui permet de conserver l’énergie en R − δ. On obtient donc pour la partie superficielle   h dT (R − δ)3 q(r) = = −k . r− (1.29) 2 3 r dr L’intégration pour la température avec, comme condition limite, une valeur nulle à la surface, donne    1 1 h R2 − r 2 3 (1.30) + (R − δ) − T = 3k 2 R r La température doit être continue en R − δ et on a donc     h 3 2 δ3 δ h R2 (R − δ)2 2 Tc = − − (R − δ) = δ − . 3k 2 2 R 3k 2 R (1.31) En ne conservant que le premier terme en δ/R, on trouve Tc = hδ 2 Qs δ . = 2k 8πkR2 (1.32) En gardant le même flux et la même conductivité thermique, on voit que ce résultat égal à δ/R fois le résultat obtenu pour une Terre homogène (eq.(1.13)). En prenant cette valeur suffisamment faible, on obtient une température centrale raisonnable. Chapitre 2 Méthodes de résolution pour la conduction 2.1 Introduction Dans ce chapitre, quelques méthodes pour la résolution du problème de la conduction de la chaleur sont présentées. L’objectif n’est pas la rigueur mathématique mais plutôt de montrer comment ces différentes approches permettent de donner des réponses dont la qualité peut dépendre des conditions du problème. Une démonstration par l’exemple est préférée à la généralisation. Un exposé systématique des différentes approches et des solutions est présenté dans le livre classique de Carslaw & Jaeger (1959). 2.2 Solutions de similitude : l’âge de la Terre selon Fourier et Kelvin Considérons le problème classique du refroidissement de la Terre selon Lord Kelvin (Thomson, 1862; Burchfield, 1975) : la Terre est supposée être formée chaude, à une température T0 , et est soudainement plongée dans un milieu à température nulle, température à laquelle est ainsi maintenue la surface. Quelle est alors la température dans la Terre en fonction du temps et de la profondeur ? Comment évolue le flux de chaleur à la surface en fonction du temps ? On s’intéresse d’abord aux temps courts, seul problème effectivement étudié par Kelvin. La température étant initialement uniforme, la température ne change au départ que près de la surface, de même que dans le cas du mur étudié précédemment (§ 1.2 et fig. 1.1). Tant que l’épaisseur affectée par le refroidissement est faible devant le rayon de la Terre, on peut négliger la sphéricité de celle-ci et considérer un demi-espace infini, la température ne variant qu’avec la profondeur z et le temps. Pour justifier cette hypothèse dans le cadre du modèle de Kelvin, rappelons en le principe. La Terre est initialement à la température T0 homogène et à t = 0 une température nulle est imposée à la surface. Le gradient de température en surface 14 Méthodes de résolution pour la conduction initialement infini diminue avec le temps et on atteint la valeur actuelle G pour un temps égal à l’âge de la Terre (Ajouter un schéma). L’intersection entre la tangente au géotherme à la surface et la température initiale est à la profondeur δ = T0 /G ∼ 60km avec les valeurs utilisées par Kelvin. Ceci nous donne un ordre de grandeur de l’épaisseur de Terre affectée par la baisse de température, épaisseur très faible par rapport au rayon de la Terre, ce qui justifie l’approximation de Terre plate. Le problème est donc de résoudre l’équation ∀t > 0 et z > 0, ∂2T ∂T =κ 2 ∂t ∂z (2.1) avec les conditions limites T (z, t < 0) = T0 , T (z = 0, t) = 0. (2.2) La première chose à remarquer est que l’équation (2.1) est linéaire en température et que l’on peut se débarrasser de l’échelle de température en divisant l’équation par T0 . La solution peut donc s’écrire T = T0 f (z, t), la fonction f étant sans dimension et satisfaisant ∂f ∂2f = κ 2. ∂t ∂z (2.3) On peut aller plus loin dans l’adimensionnement en remarquant que le problème n’a pas de dimension spatiale intrinsèque et en conséquence pas de dimension temporelle intrinsèque. Cela signifie que la solution doit prendre une forme universelle, fonction d’une variable unique sans dimension, appelée variable de similitude. Une autre façon de formuler cet argument est que la fonction f étant sans dimension elle ne peut dépendre que de variables sans dimension. Le problème dépend de deux variables indépendantes, z et t et ne contient qu’un paramètre physique, κ. Avec ces trois quantités, on ne peut former qu’une variable sans dimension, ξ = √ z/ κt, et f (z, t) = f (ξ). L’équation que doit satisfaire f est alors  2 ∂ξ ∂2ξ ∂ξ ′ ′′ ′ f (ξ) = κf (ξ) 2 + κf (ξ) , (2.4) ∂t ∂z ∂z La définition de ξ donne −ξ ∂ξ = , ∂t 2t ∂ξ 1 =√ , ∂z κt ∂2ξ =0 ∂z 2 (2.5) et la fonction f doit satisfaire l’équation différentielle suivante √ ξ ξ 2 f ′′ = − f ′ ⇒ f ′ (ξ) = f ′ (0)e−ξ /4 ⇒ f = f ′ (0) π erf + f (0). 2 2 où l’on a introduit la fonction erreur définie par Z x 2 2 erf x = √ e−u du. π 0 (2.6) (2.7) Les constantes d’intégrations f ′ (0) et f (0) sont déterminées en utilisant les conditions aux limites. En z = 0, T = 0 et cela implique f (0) = 0. Pour tout z > 0, la valeur quand t tend vers 0 est T0 . Cette valeur correspond à la limite quand ξ tend vers l’infini. Or lim erf x = 1, x→∞ (2.8) 2.2 Solutions de similitude : l’âge de la Terre selon Fourier et Kelvin 15 √ ce qui implique que f ′ (0) = T0 / π. La solution s’écrit finalement z T = T0 erf √ . 2 κt (2.9) Le flux de chaleur à la surface est obtenu en calculant la dérivée par rapport à z : q(z = 0) = −k ∂T kT0 , (z = 0) = − √ ∂z πκt (2.10) qui est bien sûr négatif puisque l’axe des z a été orienté vers le bas. Remarquons que ce résultat, √ hormis pour sa constante π aurait pu être obtenu par simple analyse dimensionnelle. En effet, √ la seule échelle de distance qui existe dans le problème est l’échelle diffusive κt et le gradient √ à la surface ne peut que s’écrire comme proportionnel à kT0 / κt. Si maintenant on s’intéresse au flux √ à une profondeur d, la situation est différente puisque deux échelles de distance existent, d et κt. On peut alors écrire que le flux de chaleur s’écrit comme   d kT0 (2.11) G √ q(z = −d) = − d κt et la forme exacte de la fonction G ne peut s’obtenir que par une résolution totale du problème, comme nous l’avons fait précédemment. O .1b .2b .3b .4b .5b .6b .7b .8b .9b A b Y 1° Fahr. per 51 feet 790° Fahr. 400 000 feet .5 P N P' ax1.0 Fig. 2.1 – Figure de Kelvin(Thomson, 1862) pour la température et le flux dans la Terre. ax1.5 R ax2.0 Q ax2.5 X 7000° Fahr. L’expression (2.10), comparée à l’observation du flux géothermique à la surface, a permis à Kelvin de calculer l’âge de la Terre comme étant 100 Ma (Thomson, 1862). (Donner des détails des paramètres utilisés par Kelvin) On sait que ce résultat est incorrect et la raison principale est la non prise en compte de la convection comme mode de transfert de chaleur principal dans l’intérieur de la Terre (Richter, 1986; Harrison, 1987). Un modèle conductif qui modélise la convection à l’aide d’une conductivité plus importante en profondeur permet d’ailleurs d’obtenir un âge pour la Terre beaucoup plus raisonnable (Perry, 1895b,c,a; Heaviside, 16 Méthodes de résolution pour la conduction 1899). La même observation a pu être faite lors de l’exercice 1. Même si le modèle de Kelvin pour le refroidissement de la Terre n’est pas correct, on verra que la solution donnée par l’équation (2.9) reste une très bonne approximation pour la température dans la lithosphère océanique. La technique appliquée ici utilise en fait un théorème important, le théorème Π, qui permet de réduire le nombre de paramètres indépendants aux paramètres sans dimension du problème. Nous n’avons pas ici exposé en détail les méthodes de l’analyse dimensionnelle, qui est un outil très puissant comme on a pu le voir dans cet exemple, et un exposé de cette approche peut être trouvé dans le livre de Barenblatt (1996). 2.3 L’âge de la Terre selon Perry Perry (1895b,c,a), ancien étudiant de Kelvin, propose une modification de son modèle et montre que l’on peut ainsi obtenir un âge pour la Terre beaucoup plus important. Le principe est similaire à celui proposé dans l’exercice ?? du chapitre précédent et repose sur l’idée d’une diffusion plus efficace en profondeur que celle utilisée dans le modèle de Kelvin. L’obtention d’une solution à la diffusion dans un espace infini composite est relativement simple en utilisant les transformées de Laplace mais suivons plutôt ici la méthode de Perry, très subtile. Considérons un semi–espace infini, homogène, ayant une conductivité thermique k1 , une diffusivité thermique κ1 , initialement à la température T1 est dont la surface supérieure z = 0 est maintenue à température nulle. La solution à ce problème a été obtenue juste au dessus et selon l’équation (2.9) en remplaçant les paramètres par ceux du présent problème. Pour un temps court ou une profondeur faible z ≤ z1 , cette solution peut être linéarisée en T0 z T =√ , πκ1 t pour z ≤ z1 . (2.12) Considérons maintenant un autre semi-espace infini ayant une conductivité thermique k2 , etc, au même temps que le précédent. On peut choisir z2 tel que la température et le flux de chaleur à la profondeur z1 du milieu 1 égalent la température et le flux de chaleur du milieux 2 à la profondeur z2 . Pour cela, en utilisant la forme linéarisé de la solution (eq. (2.12)), il suffit de choisir z2 T2 z1 T1 √ =√ κ1 κ2 (2.13) k1 T1 k2 T2 √ = √ κ1 κ2 la première de ces équations permettant d’obtenir l’égalité de température et la seconde l’égalité de flux, pour tout temps, tant que l’approximation linéaire est correcte. Une condition nécessaire qui découle des deux précédente est que √ κ1 T2 k1 z1 = =√ = N. (2.14) k2 z2 κ2 T1 Considérons maintenant un milieu composite dans lequel la partie supérieure [0, z2 ] du second milieu au sommet de la partie inférieure [z1 , +∞[ du premier. On laisse refroidir ce système 2.4 Séries de Fourier 17 et on peut se rendre compte que la fonction fabriquée en superposant les parties des deux solutions discutées précédemment est la bonne solution. En effet, chaque fonction est solution de la diffusion dans son domaine. Par ailleurs, la continuité de la température et du flux en z = z2 est assurée par construction. Et les conditions aux limites sont satisfaites. Le gradient de température à la surface, c’est à dire celui qui est accessible à la mesure, est G= √ T1 T2 = N√ πκ2 t πκ1 t (2.15) On voit ainsi que si la conductivité thermique est N plus grande en profondeur qu’à la surface, le gradient à la surface est N fois plus grand que celui qui prévaut juste sous la discontinuité, ceci afin de conserver l’énergie. Appliquons maintenant cette idée au modèle de Kelvin. La mesure du gradient géothermique et des propriétés physique des roches en surface étant inchangées, Kelvin obtient pour l’âge de la Terre T02 . tK = πκG2 (2.16) Considérons maintenant un autre semi-espace infini, avec un gradient n fois plus faible à la surface et une diffusivité thermique m fois plus forte. Le temps correspondant est donc t= n2 T02 , m πκG2 (2.17) c’est à dire n2 /m fois la valeur de Kelvin. Mais bien sûr le gradient à la surface n’est pas correct. Pour corriger cela, on remplace une couche fine à la surface par une couche ayant la diffusivité de Kelvin et le gradient de Kelvin. Il suffit de bien choisir les rapports de conductivité et d’épaisseur pour la continuité du flux soit satisfaite. En prenant n ≃ m, on peut obtenir l’âge désiré. Perry propose un facteur 56, soit un âge de 5.6 milliard d’années ! Pour justifier une diffusivité plus grande en profondeur, Perry évoque la pression plus grande ou la convection d’un fluide interstitiel. 2.4 Séries de Fourier Une autre méthode classique pour obtenir des solutions de l’équation de la diffusion est d’utiliser un développement en série de Fourier. Rien d’étonnant à cela puisque Fourier a laissé son nom dans ces deux domaines. 18 Méthodes de résolution pour la conduction 2.5 Transformée de Laplace : l’âge de la Terre selon Heaviside 2.6 Solutions approchées par la méthode intégrale Lorsqu’une solution ne peut être trouvé par aucune des méthodes exposées ci-dessus, on a souvent recours à une résolution numérique. Ces méthodes ne seront pas décrites ici mais on peut s’en inspirer pour déterminer des solutions approchées qui permettent souvent d’obtenir le comportement correct de la solution exacte. De manière générale, une solution numérique est obtenue en écrivant la solution comme un développement sur un ensemble de fonctions sur lesquelles on peut travailler facilement. Par exemple, on peut développer la fonction recherchée sur un ensemble fini de modes de Fourier et on sait que lorsque l’on augmente le nombre de modes considérés, la solution partielle va converger vers la solution complète. On peut aussi choisir de décomposer l’espace en domaines qui vont être les support de polynômes sur lesquels on développe la solution. Cela donne lieu à la méthode des éléments finis. La solution étant développée sur l’ensemble de fonction choisit, on écrit une version approchée des équations de conservation pour ce développement et cela conduit à résoudre un système d’équations pour les coefficients du développement. La qualité du développement et la rapidité de convergence vers la solution exacte dépend du choix des fonctions pour le développement et leur qualité intrinsèque pour représenter le phénomène considéré. On a vu plus haut par exemple que le choix d’une solution sous forme d’un développement de Fourier était peu approprié aux temps courts lorsque le problème présente une discontinuité entre la condition initiale et les conditions limites. Par contre, un développement en fonctions erreur est indiqué dans ce cas là. Nous allons ici proposer une approche de simplification extrême qui permet d’obtenir des solutions approchées ayant le bon comportement. La réussite de cette approche dépend de deux ingrédients : une bonne intuition concernant la forme de la solution exacte et l’écriture des équation de bilan à l’échelle globale. Pour être plus spécifique, considérons un exemple. Le modèle de refroidissement de la Terre de Kelvin peut être traité de manière approché par une telle méthode. On considère un demiespace infini, initialement à la température T0 , dont la température est maintenue à T = 0 pour tout t > 0. Un front froid se propage depuis la surface vers la profondeur et on peut modéliser le géotherme par une partie linéaire proche de la surface et une température constante égale à T0 au delà (fig. ??). L’épaisseur de la partie linéaire est notée δ et on s’attache à trouver son évolution temporelle. Entre la situation initiale et la situation au temps t, la température a baissé de ∆T (z), z étant la profondeur, ce qui correspond à l’extraction d’une énergie par unité de surface ∆Es = Z 0 ∞ 1 ρCp ∆T (z)dz = ρCp T0 δ, 2 (2.18) où l’on a supposé que la densité ρ et la capacité calorifique Cp étaient constantes. L’intégrale est simplement l’aire du triangle entre la température initiale et la température au temps t. L’extraction de cette énergie se produit par la surface supérieure, puisque le flux vers le bas est 2.6 Solutions approchées par la méthode intégrale 19 nul. À chaque instant τ , le flux à la surface est q(τ ) = k T0 , δ(τ ) et la conservation de l’énergie stipule que Z t d∆Es q(τ )dτ ⇒ ∆Es = = q(t). dt 0 Ceci donne donc une équation différentielle pour δ : √ δ δ̇ = 2κt ⇒ δ = 2 κt. (2.19) (2.20) (2.21) Le gradient de température à la surface décroı̂t donc avec le temps comme T0 G= √ , 2 κt (2.22) √ qui est assez proche du résultat exact (eq. 2.10) (aussi proche que π est proche de 2). Une telle approche ne permet pas de faire des prédictions quantitatives très précises mais permet d’obtenir un comportement. Cela peut aider à interpréter des résultats plus complets obtenus par exemple par modélisation numérique. Plus loin, nous utiliserons une approche similaire pour construire un modèle de convection thermique d’amplitude finie. Exercise 3 Cet exercice a pour objectif d’explorer le modèle de Perry pour le refroidissement de la Terre et d’en trouver une solution approchée par la méthode intégrale. Perry a apporté une modification au modèle de Kelvin en considérant que la conductivité thermique k utilisée par ce dernier est celles des roches de surface et que la matière dans la Terre profonde (en dessous d’une peau fine d’épaisseur δ) a une conductivité thermique égale à N k. Pour simplifier l’obtention d’une solution, nous allons supposer que la température varie avec la profondeur comme une fonction linéaire par morceaux. Plus précisément, on suppose que la solution a la forme suivante (z étant la profondeur) :  z  si z ≤ δ  δ T1 (t) z−δ T (z) = T1 (t) + d(t)−δ (T0 − T1 (t)) si δ ≤ z ≤ d(t) (2.23)   T0 si z ≥ d(t) T0 étant la température initiale de la Terre (uniforme), et d(t) et T1 (t) devant être déterminés. 1. Représenter schématiquement le profil de température dans la Terre selon ce modèle. 2. Quelles sont les conditions qui doivent être satisfaites à la limite z = δ ? En déduire une relation entre les variables du problème. 3. Écrire le bilan de chaleur de la Terre intégré entre le temps 0 et le temps t. En déduire une relation supplémentaire. 4. Pour résoudre le problème, on fait une hypothèse supplémentaire : on suppose que la peau peu conductrice de chaleur est infiniment fine, soit δ << d/N . En déduire des relations simplifiées entre T1 et d et en déterminer les expressions. 5. Quel est l’âge de la Terre selon ce modèle ? Comparer à celui de Kelvin et discuter de la pertinence de ce modèle. 20 Méthodes de résolution pour la conduction 2.7 Solutions aux exercices Solution 3 1. La forme du profil de température est montrée sur la figure 2.2. La température est continue mais son gradient est discontinu, de façon à conserver l’énergie, comme discuté dans la question suivante. Fig. 2.2 – Profil de température simplifié dans le modèle de Perry. La brisure de pente à z = δ est due à la discontinuité de conductivité thermique. 2. À l’interface z = δ, comme à tout autre interface, on doit satisfaire deux conditions thermiques : continuité de la température et conservation de l’énergie. La première condition est déjà prise en compte dans l’expression de la solution approchée. La seconde s’écrit k T1 T0 − T1 = Nk . δ d−δ (2.24) Cette équation peut être transformée pour donner T1 comme T1 = T0 Nδ . d + (N − 1)δ (2.25) 3. La température étant uniforme pour z > d, aucun flux de chaleur n’est à prendre en compte de ce coté là. Le bilan de chaleur équilibre donc simplement le flux à la surface le refroidissement séculaire du volume ayant une capacité calorifique Cp . Ceci s’écrit de la manière suivante Z t Z d T1 ρCp k dt = ρCp (To − T )dz = [T0 (δ + d) − T1 d] . (2.26) δ 2 0 0 On peut alors dériver cette équation par rapport au temps pour obtenir une équation différentielle qui s’écrit : i h ˙ (2.27) 2κT1 = δ (T0 − T1 ) d − dṪ1 4. On suppose maintenant que δ << d/N < d et l’équation (2.25) se simplifie en T1 = Nδ T0 d ⇒ T1 d = T0 N δ. (2.28) 2.7 Solutions aux exercices 21 On voit donc que le produit T1 d est constant et l’équation différentielle (2.27) qui fait intervenir la dérivée de ce produit se simplifie en 2κT1 = δT0 d.˙ (2.29) En éliminant T1 entre les deux équations on obtient simplement dd˙ = 2N κ. Cette équation admet comme solution, en supposant que d(0) = 0, √ (2.30) d = 2 N κt. et donc Nδ T1 = √ T0 . 2 N κt Le gradient géothermique de surface est alors r 1 N T0 . G= 2 κt (2.31) (2.32) 5. Dans ce modèle (comme dans celui de Kelvin), l’âge de la Terre est obtenu déterminant le temps nécessaire pour le gradient géothermique de surface T1 /δ pour décroı̂tre jusqu’à la valeur observée, G0 , soit t=N T02 . 4κG20 (2.33) Cette valeur est N fois plus grande que celle de Kelvin et permet éventuellement d’obtenir un âge raisonnable si la conductivité en profondeur est suffisamment importante. Physiquement, cette conductivité plus élevée permet de maintenir plus longtemps le gradient de surface en apportant de l’énergie depuis des profondeurs plus grandes dans la Terre √ (d’un facteur N ). Les raisons invoquées par Perry pour justifier une telle hypothèse étaient variées, comme bien sûr des états de la matière inconnus à haute pression mais aussi de la convection. Il pensait plutôt à la convection de fluides intersticiels qu’a de la convection à l’état solide du manteau, mais on comprend maintenant la pertinence de ce modèle. Il est à noter que, dans les années 1950, c’est à dire avant l’avènement de la tectonique des plaques, plusieurs équipes de chercheurs ont essayé d’estimer la part radiative du transfert de chaleur dans le manteau, justement pour trouver un mécanisme qui permette de réconcilier l’âge de la Terre avec son évolution thermique. 22 Méthodes de résolution pour la conduction Chapitre 3 Déformation, contraintes et rhéologie 3.1 Introduction Dans le chapitre précédent, on a vu comment le flux de chaleur était relié, de façon empirique, au gradient de température par la loi de Fourier. Le même type d’approche s’applique à la relation entre contraintes et déformations, la relation entre les deux étant appelée rhéologie. Comme dans le cas de la Loi de Fourier, la rhéologie est une relation phénoménologique et pas une relation fondamentale de la physique comme le sont les équations de conservation qui feront l’objet du chapitre suivant. Le problème est cependant plus complexe que dans le cas de la loi de Fourier qui s’applique remarquablement bien dans la plupart des situations rencontrées. La rhéologie, quant à elle, est non seulement plus complexe du fait d’un ordre tensoriel plus élevé, mais également parce que de nombreux fluides rencontrés dans la nature, notamment en géophysique de la Terre solide, s’écartent considérablement de la relation la plus simple que l’on puisse faire, à savoir la rhéologie linéaire appelée rhéologie newtonienne. 3.2 Approche unidimensionnelle Pour introduire la discussion, commençons par un cas unidimensionnel simple, analogue du cas de la diffusion dans un mur discuté au chapitre précédent. Nous allons ainsi introduire une première rhéologie simple et écrire une équation de conservation simple. Cela permettra d’identifier les étapes de l’analyse d’un problème de dynamique des fluides. V d z x Fig. 3.1 – Écoulement de Couette stationnaire. La vitesse de la plaque inférieure (z = 0) est nulle alors que celle de la plaque supérieure (z = d) est V x̂. En état stationnaire la vitesse du fluide est u = x̂V z/d. 24 Déformation, contraintes et rhéologie Considérons une couche de fluide horizontale (fig. 3.1), comprise entre deux plaques en z = 0 et en z = d, infiniment étendue dans la direction horizontale. On suppose que la plaque du bas est immobile et que celle du haut se déplace dans la direction x̂ à la vitesse V stationnaire. Cette écoulement s’appelle écoulement de Couette, du nom du physicien Maurice Couette (1858– 1943). En supposant l’absence de glissement entre le fluide et les deux plaques, ce qui est une bonne hypothèse pour un fluide visqueux, la vitesse du fluide doit être égale à celle des plaques sur les bords. On introduit ainsi les conditions limites que doit satisfaire le problème : u = 0, u = V x̂, z = 0, z = d. (3.1) (3.2) Considérons d’abord le système en état stationnaire, ce qui implique, comme on le verra plus bas, que la vitesse dans le fluide varie linéairement avec z, soit z ux = V. d (3.3) Pour maintenir ce mouvement il faut fournir un travail pour déformer le fluide qui résiste. La caractérisation de cette résistance est justement l’étude de la rhéologie du fluide. On doit donc appliquer une force sur la plaque supérieure pour la maintenir en mouvement. Cette force est d’autant plus importante que la surface est grande et est proportionnelle à celle-ci. La quantité pertinente est donc la force par unité de surface, quantité appelée contrainte. Cette contrainte, en état stationnaire, est équilibrée par la résistance du fluide à la déformation, qui est quantifiée par le gradient de vitesse. En effet, une vitesse uniforme correspond à un déplacement en bloc du fluide, sans déformation, alors que dans le cas étudié ici, il existe un gradient dans la vitesse horizontale qui correspond à une déformation du fluide à une certaine vitesse. L’hypothèse la plus simple que l’on puisse est que pour une vitesse de déformation donnée V /d, la contrainte F/A à appliquer (F étant la force et A l’aire considérée) lui est reliée linéairement : V F =η , A d (3.4) où l’on a introduit un coefficient, η appelé viscosité dynamique. La contrainte est une force par unité de surface, comme la pression qui est un cas particulier de contrainte (nous y reviendrons), et a donc pour unité le Pascal, Pa = N m−2 . La vitesse de déformation, V /d, a pour unité la s−1 et la viscosité dynamique a donc pour unité le Pa s. V d z x Fig. 3.2 – Mise en route de l’écoulement de Couette. La vitesse de la plaque inférieure (z = 0) est nulle alors que celle de la plaque supérieure (z = d) est V x̂. La partie supérieure se met en mouvement d’abord la quantité se transmet de proche jusqu’à atteindre la plaque inférieure. L’équation (3.4) n’est valide qu’une fois l’état stationnaire atteint. Que se passe-t-il avant d’y arriver ? Dès que l’on met en mouvement la plaque supérieure, la condition de non–glissement (eq.(3.2)) doit être satisfaite et le fluide près du bord prend la vitesse de la plaque, le reste du fluide sous-jacent étant encore au repos. Cela signifie que le fluide est très rapidement 3.2 Approche unidimensionnelle 25 déformé et, la viscosité agit pour transférer la quantité de mouvement vers le bas, de proche en proche (fig.3.2). Clairement, l’équation (3.4) n’est plus la bonne pour déterminer la contrainte nécessaire pour maintenir ce mouvement. Par contre, si l’on considère la tranche de fluide entre d − δ et d, d’épaisseur δ suffisamment faible pour que le profil de vitesse puisse y être considéré comme linéaire, une équation de la même forme peut être écrite. En d’autre terme, la contrainte est reliée linéairement à ∂z ux . Avant d’expliciter cette relation, in convient de stabiliser les notations pour la contrainte. On voit que celle-ci est appliquée sur une surface horizontale, c’est à dire perpendiculaire à la direction ẑ et doit être appliquée dans la direction x̂. On notera donc cette contrainte σxz . La relation entre cette contrainte et la vitesse de déformation du fluide s’écrit maintenant σxz = η ∂ux . ∂z (3.5) Cette loi rhéologique est appelée rhéologie Newtonienne. C’est la plus simple que l’on puisse faire et on verra plus loin comment la généraliser. Écrivons maintenant le bilan de la quantité de mouvement pour une tranche de fluide comprise entre z et z + δz. La loi de Newton de conservation de la quantité de mouvement nous dit que la variation de la quantité de mouvement contenue dans cette tranche est égale à la somme des forces appliquées. Dans le cas d’une région fluide il convient également prendre en compte le transport de la quantité de mouvement par l’écoulement lui-même, mais dans le cas présent, l’écoulement étant horizontal, il n’induit ni gain ni perte de quantité de mouvement pour la tranche de fluide considérée, puisqu’elle est horizontale et d’extension infinie. La loi de Newton est vectorielle mais seule sa composante selon x̂ est intéressante dans le cas présent. La force de gravité n’intervient donc pas et seules les contraintes visqueuses sur les bords horizontaux mettent le fluide en mouvement. La tranche considérée est infinie horizontalement mais on peut l’écrire pour une unité de surface. s’écrit par unité de surface horizontale : ∂ (ρux δz) = σxz (z + δz) − σxz (z). ∂t (3.6) Pour comprendre la différence de signe entre les deux termes du membre de droite, il faut se rappeler que les surfaces sont orientée, de l’intérieur de la zone considérée vers l’extérieur. Donc, en z + δz, la surface est orientée dans le sens de ẑ alors qu’en z, elle est orientée dans le sens de −ẑ. Évidemment, on divise cette équation par δz que l’on fait tendre vers 0 pour obtenir la forme locale de l’équation de conservation :   ∂ux ∂σxz ∂ ∂ η . (3.7) (ρux ) = = ∂t ∂z ∂z ∂z En supposant maintenant que ρ et η sont constants, on peut écrire cette équation sous la forme η ∂ 2 ux ∂ux = . ∂t ρ ∂z 2 (3.8) On voit alors que la vitesse ux est solution d’une équation de diffusion avec pour diffusivité ν = η/ρ, qui est appelé viscosité cinématique. Cela montre que les résultats obtenus pour la température peuvent être transposés au cas de la quantité de mouvement. 26 Déformation, contraintes et rhéologie Ce premier exemple nous a permis de mettre en évidence plusieurs points qu’il s’agit maintenant d’explorer complètement. D’abord, la déformation du fluide est mesurée par le gradient du champ de vitesse et que l’on peut chercher une relation linéaire entre celui-ci et le tenseur des contraintes. Exercise 4 Considérons un écoulement de lave dans un conduit volcanique vertical que nous supposerons parfaitement cylindrique (rayon R). Comment s’écrit la loi de conservation de la quantité de mouvement dans la direction ẑ ? (Pour la résoudre, il faudra écrire une loi rhéologique plus générale que celle déjà écrite). 3.3 Vitesse de déformation Supposons maintenant que l’on puisse se déplacer dans le champ de vitesse pour le mesurer en différents points, sans le perturber. Pour être pratique, cela est réalisable dans une expérience de laboratoire en utilisant une sonde doppler. Donc, supposons que l’on déplace la position du point où l’on mesure la vitesse du fluide, selon une trajectoire r(t). Le champ de vitesse est stationnaire mais non–uniforme. Donc, lorsque la sonde se déplace, elle mesure des vitesses variables. Plus précisément, la composante selon x̂ de la vitesse, ux , varie selon l’équation : dr ∂ux dx ∂ux dy ∂ux dz dux = ·∇ux = + + . dt dt ∂x dt ∂y dt ∂z dt (3.9) Une équation similaire peut être écrite pour les autres composantes et la variation du vecteur complète peut finalement s’écrire     dt x ∂x ux ∂y ux ∂z ux dr du  ∂x uy ∂y uy ∂z uy  ·  dt y  = G· . = (3.10) dt dt dt z ∂x uz ∂y uz ∂z uz G est un tenseur d’ordre (ou rang) 2 et est représenté en tout point de l’espace par une matrice 3 × 3 (ou 2 × 2 en 2D) dont les valeurs varient avec la position. Ce tenseur s’appelle tenseur gradient de vitesse. ou tenseur des vitesses de déformation 1 . Dans le cas d’une vitesse uniforme et donc sans déformation, on voit que ce tenseur est nul. Remarquons que l’on peut le décomposer en une partie symétrique et une partie anti-symétrique : ∂ui 1 Gij = = ∂xj 2  ∂ui ∂uj + ∂xj ∂xi  1 + 2  ∂uj ∂ui − ∂xj ∂xi  . (3.11) On note eij = (∂j ui + ∂i uj )/2 la partie symétrique du gradient de vitesses et ωij = (∂j ui − ∂i uj )/2 sa partie anti–symétrique. Pour visualiser la signification de chaque partie, nous allons considérer une situation bi–dimensionnelle. Considérons pour cela deux points suffisamment 1 Notons que dans la littérature anglo-saxonne, ce tenseur est appelé strain rate tensor que l’on traduit souvent tenseur des taux de déformation en hydrodynamique (voir par exemple ?). Cette dénomination est ambiguë et nous préférerons celles données dans le texte principal. 3.3 Vitesse de déformation 27 proches, A et B, et leurs positions A′ et B ′ après un temps δt. La déformation peut être quantifiée par la différence vectorielle suivante : A′ B′ − AB = A′ A + AB + BB′ − AB = BB′ − AA′ = [u(B) − u(A)] δt. (3.12) Si l’on considère des points initialement proches, c’est à dire AB = (δx, δy)T , le vecteur entre crochets dans le membre de droite de l’équation (3.12) peut être exprimé comme développement au premier ordre en δx et δy : A′ B′ − AB ≃ [∂x uδx + ∂y uδy] δt. En introduisant les deux composantes de la vitesse, u = ux x̂ + uy ŷ, on obtient A′ B′ − AB ≃ [(∂x ux δx + ∂y ux δy) x̂ + (∂y uy δy + ∂x uy δx) ŷ] δt = G · ABδt (3.13) et on retrouve les composantes du tenseur gradient de vitesse. Cas particulier no 1 : ∂y ux = ∂x uy = 0 (pas de terme hors diagonale). Si l’on considère un rectangle dont les cotés sont alignés avec les axes du repère (fig. 3.3), on peut facilement calculer l’évolution de ces cotés. Ainsi, l’application de l’équation (3.13) avec AB = (δx, 0)T donne A′ B′ − AB ≃ ∂x ux δxδtx̂, c’est à dire qu’un vecteur aligné avec l’axe x̂ ne subit pas de rotation. On peut vérifier aisément qu’il en est de même pour le vecteur AC et que l’on peut donc généraliser : un écoulement caractérisé par un gradient de vitesse sans terme hors de la diagonale ne produit aucune rotation. C D δt δt C' D' y x δt A' A δt B' B Fig. 3.3 – Déformation dans le cas où le tenseur gradient de vitesse ne contient des termes non-nuls que sur la diagonale. Les traits discontinus représentent l’évolution pendant le temps δt Considérons maintenant la variation de longueur des cotés du rectangle : ∂x ux δxδt δ(AB) = = ∂x ux δt, AB δx δ(AC) = ∂y uy δt, AC 28 Déformation, contraintes et rhéologie et la variation relative de la surface s’écrit δ(AB) δ(AC) δS = + = (∂x ux + ∂y uy )δt = ∇ · uδt. S AB AC Ce résultat peut facilement être étendu au cas tridimensionnel et on voit ainsi que la trace du tenseur G, qui n’est autre que la divergence du vecteur vitesse, mesure la vitesse d’expansion ou de contraction relative dans un écoulement de fluide. ∂vx δyδt ∂y C' δβ C B' A' y δα ∂vy δxδt ∂x Fig. 3.4 – Déformation dans le cas où le tenseur gradient de vitesse contient des termes non–nuls uniquement hors de la diagonale. x A B Cas particulier no 2 : ∂x ux = ∂y uy = 0 (pas de terme diagonaux). En partant du même rectangle initial que précédemment, on voit maintenant que A′ B′ − AB ≃ ∂x uy δxδtŷ, ce qui signifie que l’on a ajouté au vecteur selon x̂ une composante selon ŷ, c’est à dire qu’il a subit une rotation, d’un angle δα. En supposant que l’angle δα est petit, soit tan δα ≃ δα, on estime cet angle à δα = ∂x uy δt. De même, en considérant le coté AC = (0, δy)T , on trouve qu’il tourne d’un angle δβ = −∂y ux δt et on peut calculer la variation de l’angle γ, initialement égal à π/2, entre AB et AC : δα − δβ δγ =− = −(∂x uy + ∂y ux ) = −2exy . δt δt On fait ainsi apparaı̂tre les termes non–diagonaux de la partie symétrique du gradient de vitesse, qui peuvent être ainsi associés à la vitesse de cisaillement du fluide. En l’absence de termes diagonaux, comme on le suppose ici, il n’y a pas de changement de volume (puisque celui-ci est quantifié par la divergence de la vitesse) et on a affaire à un cisaillement pur. Analyse des termes anti–symétriques. On suppose ici que la composante symétrique est nulle pour isoler la signification de la partie anti–symétrique. En particulier, les termes diagonaux sont nuls et leur somme, la divergence de la vitesse, l’est également. On a déjà vu que les deux cotés du rectangle de départ subissaient une rotation, avec les angles δα = ∂x uy δt et δβ = −∂y ux δt. Or, puisque exy = 0 = (∂y ux + ∂x uy )/2, on voit que δα = δβ (fig. 3.5). Autrement dit, l’angle de rotation est indépendant de la direction initiale : on a affaire à une 3.3 Vitesse de déformation 29 C' C D δα D' Fig. 3.5 – Évolution dans le cas où le tenseur gradient de vitesse est antisymétrique. y A' x δα B' A B rotation pure ou rotation solide. La vitesse angulaire de cette rotation est δα 1 = ∂x uy = (∂x uy − ∂y ux ) = ωyx . δt 2 Remarquons que ωyx = −(∇ × u)z /2 et bien sûr ωxy = (∇ × u)z /2. Le vecteur ~ ω ≡ ∇ × u est appelé vorticité de l’écoulement et est relié à la rotation du fluide sur lui même. La vitesse de rotation instantanée est caractérisée par le vecteur tourbillon Ω = ∇ × u/2. L’utilisation du même symbole pour noter la vorticité et la partie anti–symétrique du tenseur gradient de vitesse ne fait que refléter leur signification très proche, au facteur 1/2 et à l’ordre tensoriel près. On relie ces deux objets par 1 ωk = − εijk ωij , 2 (3.14) où le tenseur d’ordre 3, εijk , est appelé symbole de Levi-Civita a été introduit. Il vaut zéro si deux de ses indices sont égaux, vaux 1 s’ils sont une permutation directe de (1, 2, 3) et −1 dans le cas d’une permutation indirecte, comme par exemple (1, 3, 2). Résumé et conclusion. Le gradient de vitesse peut être décomposé en trois parties ayant chacune une interprétation claire : Gij = eij + ωij = tij + dij + ωij , avec 1 tij = ∇ · uδij 3 et dij = eij − tij . (3.15) tij est le tenseur qui quantifie les changements isotropes2 de volume, dij le reste de la partie symétrique, qui quantifie le cisaillement sans changement de volume et ωij qui quantifie la rotation pure. Dans le cas où tout changement de volume est interdit (fluide incompressible), la trace Gii = tii = ∇ · u = 0. 2 C’est à dire indépendant de la direction 30 Déformation, contraintes et rhéologie 3.4 Contraintes dans un fluide On a définit plus haut (§ 3.2) une contrainte avec deux indice, un correspondant à la direction dans laquelle elle s’exerce et l’autre correspondant à la direction normale à la surface sur laquelle elle s’exerce. Il s’agit ici de généraliser cette définition. Ainsi, en considérant un élément de surface dans le plan (ŷ, ẑ), c’est à dire normal à la direction x̂, trois contraintes peuvent s’exercer, dans les trois direction de l’espace : σzx dans la direction ẑ, σyx dans la direction ŷ et σxx dans la direction x̂. De même pour une surface orientée dans la direction ŷ ou ẑ, ce qui définit neuf coefficients différents au tenseur des contraintes. Nous allons ici étudier les propriétés de ce tenseur. Dans un premier temps, les composantes du tenseur ont été définies par rapport à des éléments de surface particulier, orientés perpendiculairement aux vecteurs de la base x̂, ŷ, ẑ. Comment en déduire la contrainte exercée sur un élément de surface quelconque, δS = δSn̂ ? Sans perte de généralité, la force totale δf exercée sur la surface, proportionnelle à son aire δS, a trois composante et s’écrit δf = δS (σxn x̂ + σyn ŷ + σzn ẑ) . (3.16) Il s’agit maintenant de déterminer les trois composantes σin en fonctions du tenseur σij . y y δy n θ δz δx n δy θ x z δx x z Fig. 3.6 – Tétraèdre formé de l’élément de surface δS et de ses projections sur les plans de la base (gauche) et sa projection dans le plan (x̂, ŷ) (droite). Nous allons écrire le bilan des forces appliquées sur le tétraèdre formé par l’élément de surface δS et ses projections sur les trois plans de la base (fig.3.6). Dans la direction x̂ s’exercent σxn δSx̂ et trois autres forces, sur les trois autres faces du tétraèdre. Prenons, pour commencer, le cas de la face perpendiculaire à ŷ. La contrainte concernée est, par définition, σxy , mais la surface étant orientée vers l’extérieur du tétraèdre, est orientée dans la direction −ŷ. La contribution de cette contrainte apparaı̂tra donc en négatif. Cette face est la projection de l’élément de surface δS sur le plan (x̂, ẑ) et son aire est donc simplement δS · ŷ = δSny (fig.3.6). Donc la force résultante sur la face perpendiculaire à ŷ est −σxy ny δS. De même, on obtient −σxx nx δS et −σxz nz δS sur les faces perpendiculaires à x̂ et ẑ, respectivement. 3.4 Contraintes dans un fluide 31 En plus des forces de surface considérées ci-dessus, le système peut être soumis à une force de volume f δV , comme la gravité, et le bilan de la quantité de mouvement sur notre petit volume s’écrit donc, dans la direction x̂ d (ρux ) δV = fx δV + (σxn − σxx nx − σxy ny − σxz nz ) δS, dt (3.17) et des équations similaires dans les direction ŷ et ẑ. Comme souvent, nous allons faire tendre les longueurs vers 0 et voir ce qu’il se passe. Pour éviter que les directions ne changent, nous faisons tendre la longueur des arrêtes du tétraèdre vers 0 à la même vitesse ou, autrement dit, l’élément de volume tend vers une taille nulle de manière homothétique. On peut diviser l’équation (3.17) par δS pour obtenir   d δV (ρux ) − fx = σxn − σxx nx − σxy ny − σxz nz . (3.18) dt δS Clairement, le volume tend vers 0 plus rapidement que la surface et on en déduit que σxn = σxx nx + σxy ny + σxz nz . Le même travail simplement que  σxx  σn = σyx σzx (3.19) peut être effectué pour les forces dans les direction ŷ et ẑ et on en déduit   nx σxy σxz   ny  σyy σyz nz σzy σzz (3.20) où l’on a utilisé la notation classique d’un produit matriciel. Il est important de réaliser cependant que σ n’est pas une matrice mais un tenseur, le tenseur des contraintes. Montrons y σyz δy σzy σzy σyz x Fig. 3.7 – Volume élémentaire sur lequel faire le bilan des moments pour démontrer la symétrie du tenseur des contraintes. δz δx x z maintenant que le tenseur des contraintes est symétrique. Pour ce faire, nous allons écrire le bilan des moments cinétiques sur une volume élémentaire, comme représenté sur la figure 3.7, centré au point de coordonnées (x, y, z). Considérons le moment des forces appliquées sur le volume par rapport à l’axe x̂ passant en son centre, le cas des autres axes étant traité de la même façon. Parmi les forces de surface appliquées, toutes ne contribuent pas au moment total, comme toutes les forces appliquées dans la direction x̂. Par ailleurs, les forces tangentes appliquées sur les faces perpendiculaires à x̂ ont une résultante qui passe par l’axe de rotation 32 Déformation, contraintes et rhéologie considéré et ne contribuent donc pas au moment. Il ne reste à considérer que σyz et σzy appliqués sur les faces perpendiculaires à ẑ et ŷ, respectivement. Considérons d’abord la face perpendiculaire à ẑ située en z + δz. Le vecteur surface est δS = δxδyẑ et la force à considérer est donc δf = σyz (x, y, z + δz/2)δxδyŷ. Si σyz > 0 cette force a tendance à faire tourner le cube dans le sens indirect par rapport à l’axe x̂. La contribution au moment total sera donc opposée au signe de σyz . La longueur du bras de levier est δz/2 et la contribution est donc −σyz (x, y, z + δz/2)δxδyδz/2. Considérons maintenant la face opposée, en z − δz/2. Le vecteur normal est maintenant −δS = δxδyẑ et la force est δf = −σyz (x, y, z − δz/2)δxδyŷ. On voit donc que si la contrainte est de même signe, la force est de sens opposée à celle considérée précédemment. Par contre, comme elle est appliquée sur une face qui se trouve de l’autre coté de l’axe, la contribution au moment est de même signe et vaut −σyz (x, y, z − δz/2)δxδyδz/2. Considérons à présent les deux forces selon ẑ appliquées sur les faces perpendiculaires à ŷ. La même analyse que précédemment permet de voir que sur la face en y + δy/2, la contribution au moment total est du même signe que σzy et vaut σzy (x, y + δy/2, z)δxδyδz/2. De même en y − δy/2 on obtient σzy (x, y − δy/2, z)δxδyδz/2 si bien que la résultante du moment autour de l’axe x̂ s’écrit  σzy (x, y − δy/2, z) + σzy (x, y + δy/2, z) Γx = δV 2  (3.21) σyz (x, y, z − δz/2) + σyz (x, y, z + δz/2) , − 2 avec δV = δxδyδz. Par la suite nous allons faire tendre le volume vers 0 et cette expression va donc tendre vers Γx = δV (σzy − σyz ) . (3.22) Écrivons maintenant le bilan des moments pour notre volume élémentaire, dont le moment d’inertie autour de l’axe x̂ est noté δIx : δIx dΩx = Γx dt (3.23) avec Ωx la vitesse angulaire autour de l’axe x̂ passant par le centre de notre volume. Aucune force de volume n’est considérée ici puisque son point d’application est sur l’axe et n’exerce donc aucun moment. Le moment d’inertie Ix d’un objet par rapport à l’axe x̂ est définit comme Ix = Z V kx̂ × OMk2 ρdV, (3.24) O étant un point sur l’axe et M le point courant dans le calcul de l’intégrale. Dans le cas de notre volume élémentaire, il n’est pas nécessaire de faire un calcul explicit pour pouvoir affirmer que le moment d’inertie qui nous intéresse s’écrit δIx = a(δy 2 + δz 2 )ρδV , a étant un nombre sans dimension caractérisant la forme de notre volume. Le bilan du moment angulaire (3.23) s’écrit donc après simplification par le volume a(δy 2 + δz 2 )ρ dΩx = σzy − σyz . dt (3.25) 3.5 Rhéologie newtonienne 33 En faisant tendre vers 0 la taille du volume considéré, on voit que le terme de gauche tend vers 0, ce qui implique σzy = σyz . La même démonstration peut être faite pour la rotation autour des autres axes et cela démontre que le tenseur des contraintes est symétrique. 3.5 Rhéologie newtonienne On a déjà introduit la rhéologie newtonienne dans le cas uni–dimensionnel (eq.3.5) et il s’agit maintenant d’établir le cas général d’une rhéologie linéaire reliant le tenseur des contraintes au tenseur des gradients de vitesse. On a vu (§ 3.3) que ce dernier pouvait être décomposé en une partie anti–symétrique, ωij et une partie symétrique, eij . Le tenseur ω est associé à une rotation d’ensemble et pas à une déformation. Il ne doit donc pas intervenir dans la rhéologie. Une autre façon de le voir concerne la symétrie des tenseurs : le tenseur des contraintes est symétrique et on ne peut donc le relier linéairement à un tenseur anti–symétrique. On cherche donc une relation linéaire entre σij et eij . Le tenseur des contraintes que l’on a définit contient les contributions de toutes les forces de surface. Or, on connaissait avant une force de surface, la pression, qui semblait beaucoup plus simple mais qui doit contribuer au tenseur total que l’on vient de définir. Ce qui rend la pression plus simple, c’est qu’elle s’applique de la même manière dans toutes les direction. En effet, pour tout élément de surface δS, la force de pression est simplement −pδS. Autrement dit, la force s’applique perpendiculairement à la surface, indépendamment de la direction de celle-ci. Il s’agit donc d’une contrainte isotrope et un scalaire est suffisant pour la caractériser. Cependant, on peut l’écrire à l’aide d’un tenseur d’ordre 2, −pδ, δ étant le tenseur identité (symbole de Kronecker). Comme on a dit que le tenseur σ contenait toutes les contraintes, il semble utile d’en extraire la partie qui est due à la pression. Celle-ci est définie comme la partie isotrope du tenseur, le reste étant appelé contrainte visqueuse(il faudrait ici ajouter une note sur les contraintes déviatoriques), τ : σ = −pδ + τ . (3.26) La relation linéaire la plus générale entre deux tenseurs d’ordre deux fait apparaı̂tre un tenseur d’ordre 4 τij = ηijkl ekl (3.27) qui est la rhéologie Newtonienne correcte pour un milieu isotrope. C’est sans doute le cas du manteau terrestre qui est un solide cristallin mais nous ne poursuivrons pas plus avant l’étude de cette rhéologie générale. Dans le cas isotrope, pertinent pour la plupart des liquides et généralement utilisé en dynamique des fluides, on doit faire apparaı̂tre deux coefficients :   ∂uk 2 ∂uk ∂ui ∂uj +ζ + − . (3.28) τij = 2ηdij + 3ζtij = η ∂xj ∂xi 3 ∂xk ∂xk η est la viscosité dynamique que l’on a déjà rencontré et ζ est appelé viscosité seconde, ou en anglais bulk viscosity. Elle exprime une résistance visqueuse au changement de volume. Dans le 34 Déformation, contraintes et rhéologie cas où le fluide est incompressible, ∇·u = 0, ce qui implique que l’équation (3.28) se réduit à   ∂ui ∂uj . (3.29) + τij = η ∂xj ∂xi Nous somme maintenant équipés pour obtenir les équations générales de conservation. 3.6 Solutions aux exercices Solution 4 On considère la couronne cylindrique entre z et z + δz et entre r et r + δr. Sa quantité de mouvement totale selon ẑ est ρ2πrδrδzuz et sa vitesse de variation est égale à la somme des forces qui lui sont appliquées. Dans la direction ẑ, la gravité s’applique vers le bas et s’écrit −ρ2πrδrδzg. Les contraintes s’appliquent sur toutes les surfaces, deux horizontales et deux verticales. Le bilan complet s’écrit donc d (ρ2πrδrδzuz ) = − ρ2πrδrδzg dt − σzz (r, z)2πrδr + σzz (r, z + δz)2πrδr − σzr (r, z)2πrδz + σzr (r + δr, z)2π(r + δr)δz. (3.30) Avant d’aller plus loin, plusieurs choses sont à noter ici. D’abord, dans le membre de gauche, il est important de bien considérer une dérivée totale, et pas une dérivée partielle, car la quantité de mouvement est transporté par convection, c’est à dire par la vitesse. Ensuite, il est important de bien faire attention à la surface sur laquelle s’applique chaque contrainte. σzz s’applique à une surface perpendiculaire à ẑ, et les deux surface considérées ont pour aire 2πrδr. σzr s’applique à une surface perpendiculaire à r̂ mais cette fois les deux surface considérées ont une aire différente, 2πrδz et 2π(r + δr)δz. Enfin, il est important de bien comprendre d’où viennent les signes. σzz s’applique dans la direction ẑ (celle qui nous intéresse) mais sur une surface orientée dans la direction de ẑ. C’est le cas de la surface supérieure à notre cylindre (z + δz) mais la surface inférieure est orientée dans la direction opposée et la contrainte appliquée est donc opposée à σzz . Le même raisonnement s’applique pour les deux surfaces verticales. Ceci étant fait, on divise l’équation (3.30) par 2πδrδz et on fait tendre ces deux incréments vers 0. On obtient alors d ∂σzz ∂ (ρruz ) = −ρgr + r + (rσzr ) . dt ∂z ∂r (3.31) Il est important de noter le r dans la dérivée du dernier terme qui doit nous rappeler la divergence en symétrie cylindrique. Pour aller plus loin, il faut exprimer les deux contrainte en fonction des dérivées de la vitesse. Pour l’instant, on ne l’a fait que pour un cas particulier et il est temps de généraliser ce travail. Chapitre 4 Équations de bilan 4.1 Introduction Après avoir traité plusieurs problèmes de transfert de chaleur particuliers, nous sommes près à généraliser l’approche employée. On se place pour cela dans le cadre de la thermodynamique hors équilibre (de Groot & Mazur, 1983) et on écrit les bilans qui doivent être satisfaits. Le bilan qui nous intéresse plus particulièrement et qui a été utilisé précédemment (§ 1.3) est celui de la chaleur mais il ne peut être obtenu proprement sans établir les bilans pour la masse et la quantité de mouvement. De façon générale, nous allons nous intéresser à des fonctions définies en tout point du fluide qui seront des densités massiques de variables extensives, telles que la relation entre la quantité globale F correspondant à un volume V et la densité massique correspondante f s’écrive simplement : Z ρf dV, (4.1) F = V ρ étant la densité de masse. Nous allons écrire des équations de bilan pour plusieurs de ces quantité, en particulier, la masse (M, ρ), la quantité de mouvement, l’énergie (E, e), l’entropie (S, s). L’obtention de ces équations de bilan est tout à fait standard et se trouve dans divers livres de mécanique des fluides (e.g. Landau & Lifchitz, 1994). Nous nous plaçons ici dans le cadre de la thermodynamique hors équilibre (de Groot & Mazur, 1983; Kondepudi & Prigogine, 1998) : On considère ici un système hors équilibre mais dans lequel il existe un équilibre local, si bien que les relation de la thermodynamique classique à l’équilibre sont valide à l’échelle de la parcelle de fluide. On appelle parcelle de fluide un élément de volume de taille suffisamment petite pour qu’elle soit considérée comme élémentaire du point de vue du système global et des équations de bilan (et donc avec une valeur bien défnie de chaque grandeur intensive et donc à l’équilibre) mais suffisamment grande pour contenir un très grand nombre de particules afin que les équations de la physique statistique classique soient applicables. Il s’agit donc d’une échelle intermédiaire entre l’échelle microscopique et la taille du domaine étudié. 36 Équations de bilan 4.2 Conservation de la masse Considérons un volume V constant de fluide, que l’on appelle volume de contrôle. La masse qu’il contient ne peut varier que du fait d’un mouvement de fluide traversant la surface qui l’entoure, notée S (c’est à dire qu’il ne peut y avoir ni création ni disparition de masse). En chaque point de la surface, un vecteur normal n1 peut être défini, dirigé, par convention, vers l’extérieur. Si le champ de vitesse du fluide est noté u ≡ (u, v, w), la masse dm traversant la surface dS ≡ dSn pendant le temps dt est dm = −ρdtu · dS. Le signe de cet expression provient de la convention choisie pour l’orientation du vecteur de surface. En effet, une vitesse de fluide dans la même direction que ce vecteur (c’est à dire vers l’extérieur) correspond à une diminution de la masse. Par ailleurs, · signifie ici le produit scalaire, c’est à dire la réduction d’un indice tensoriel. Pour obtenir le taux de variation total de la masse contenue dans le volume V , il faut intégrer cet élément de masse sur toute la surface : d dt Z V ρdV = Z V ∂ρ dV = − ∂t Z S ρu · dS (4.2) Le terme de droite de cette équation peut être transformé en intégrale de volume en utilisant le théorème de Gauss, pour obtenir  Z  ∂ρ + ∇ · (ρu) = 0 (4.3) ∂t V Cette équation intégrale étant valable pour tout volume de contrôle V , son intégrant doit être nul, soit ∂ρ + ∇ · (ρu) = 0. ∂t (4.4) Cette équation qui exprime la conservation de la masse de manière locale est aussi appelée équation de continuité. La divergence qui apparaı̂t dans le seconde terme de l’équation de continuité peut être développée pour obtenir : Dρ + ρ∇ · u = 0 Dt (4.5) ∂ D ≡ +u·∇ Dt ∂t (4.6) avec la dérivée matérielle. Cette dérivée prend en compte la variation locale sans mouvement (∂/∂t) ainsi que la variation due au mouvement de matière qui emporte ses propriété. Pour cette raison, elle est parfois appelée dérivée totale ou encore dérivée particulaire (c’est à dire en suivant la particule de fluide). 1 Les vecteurs sont notés en caractère gras dans ce document. 4.3 Forme générale d’une équation de bilan 37 Dans le cas d’un fluide incompressible, la dérivée matérielle de la densité est nulle, c’est à dire qu’une parcelle de fluide emportée ne change pas de masse. Cela implique ∇ · u = 0 pour un fluide incompressible. (4.7) Cette condition ne sera généralement pas utilisée pour l’intérieur de la Terre, la compressibilité étant à prendre en compte dans son refroidissement. L’équation de continuité sera couramment utilisée pour simplifier les équations bilan des différentes quantités qui nous intéressent. En effet, si l’on considère une quantité F de densité massique f , on peut écrire ∂ Df (ρf ) + ∇ · (ρf u) = ρ . ∂t Dt (4.8) Pour un volume matériel V immobile, c’est à dire un volume qui est tel qu’aucun mouvement ne traverse sa surface (u · n = 0 en tout point de sa surface), on peut donc écrire Z Z Z ∂ d Df ρf dV = (ρf ) dV = ρ dV pour un volume matériel constant, (4.9) dt V Dt V ∂t V puisque l’intégrale du second terme de l’équation (4.8) peut être transformé en intégrale de surface, qui est alors nulle dans le cas d’un volume matériel. Cette équation de bilan pour un volume matériel est parfois appelé théorème de transport de Reynolds et nous sera très utile pour les versions intégrées sur une volume en convection (par exemple le noyau ou le manteau) des équations de bilan. 4.3 Forme générale d’une équation de bilan Si l’on considère la quantité F ayant pour densité massique f on peut, de manière tout à fait générale, considérer que la variation de son contenu dans le volume V provient de deux contributions : la création (ou destruction) de cette quantité avec un taux volumique χf et l’échange avec l’environnement avec un flux I′f . Ainsi, l’équation de bilan s’écrit Z Z Z Z ∂ d ′ ρf dV ≡ (ρf ) dV = − If · dS + χf dV. (4.10) dt V V ∂t S V On peut alors transformer l’intégrale de surface en intégrale de volume, Z Z Z ∂ ′ χf dV, ∇ · If dV + (ρf ) dV = − V V ∂t V et celle–ci étant valable pour tout volume de contrôle, l’intégrant doit être nul, ce qui nous donne la forme locale de la l’équation de bilan pour f : ∂ (ρf ) = −∇ · I′f + χf . ∂t (4.11) On dit qu’une quantité est conservée si le terme de production χf = 0. Ceci est le cas pour la masse mais, comme on le verra plus loin, pas pour l’entropie, par exemple. Il n’y a donc pas 38 Équations de bilan de conservation de l’entropie et l’équation correspondante est une équation de bilan plutôt que de conservation. Néanmoins, on utilise souvent le terme d’équation de conservation de façon générique, que la quantité considérée soit conservée ou non. Une partie du vecteur de flux provient du transport par mouvement au travers de la surface, I′f = ρf u + If , (4.12) où l’on a noté Ia le reste du flux, c’est à dire ce qui n’est pas dû au mouvement de matière. Introduisant cette expression dans l’équation de bilan locale (4.11), et après ré–arrangement et utilisation de l’équation (4.8), on obtient : Df = −∇ · If + χf . (4.13) Dt Nous allons donc chercher les expressions des flux et des taux de création pour la quantité de mouvement, l’énergie, et l’entropie. Ces équations pourront alors être intégrées sur le volume d’un système convectif comme le noyau ou le manteau pour étudier le rendement de la convection et le refroidissement de la Terre. ρ 4.4 Bilan de la quantité de mouvement La densité massique de quantité de mouvement est simplement la vitesse u, c’est à dire que R la quantité de mouvement contenue dans le volume V est simplement V ρudV . Cette quantité varie pour deux raisons : le transport par l’écoulement et les forces appliquées sur le fluide contenu dans V (principe fondamental de la dynamique). Ces forces sont de deux types : les forces de volume (par exemple le poids) et les forces de surface, dont le taux par unité de surface est appelée contrainte. Le tenseur des contraintes2 total est noté σ et est tel que la force appliquée sur un élément de surface dS est σ · dS. L’équation de bilan de la quantité de mouvement pour le volume de contrôle V est donc simplement : Z Z Z Z ∂ ρgdV. (4.14) (ρu) dV = − ρu (u · dS) + σ · dS + V S V ∂t S Bien sûr, le terme de flux par mouvement (première intégrale dans le membre de droite) peut être rassemblé avec le membre de gauche pour faire apparaı̂tre une dérivée totale, comme expliqué dans la section précédente3 (§ 4.3). Par ailleurs, la seconde intégrale de surface dans le membre de droite est aussi transformée en intégrale de volume en utilisant le théorème de Gauss. Comme il s’agit ici d’un terme tensoriel, explicitons cette procédure en utilisant la notation indicielle : Z Z Z Z σ · dS = σij dSj = ∂j σij dV = ∇ · σdV, (4.15) S S V V ce qui définit la divergence du tenseur des contraintes. Enfin, l’équation obtenue étant valide pour tout volume de contrôle, on obtient une équation concernant les intégrants : ρ Du = ∇ · σ + ρg. Dt (4.16) 2 Les tenseurs d’ordre 2 sont notés comme les tenseurs d’ordre 1 (les vecteurs) par un une police de caractère grasse. Le contexte est suffisant pour faire la différence. 3 La quantité transportée étant ici un vecteur, il faut appliquer la procédure à chaque composante du vecteur. 4.5 Bilan de l’énergie 39 On utilise maintenant la décomposition du tenseur des contraintes en deux parties, la pression et la contrainte visqueuse, comme discuté plus haut (§ 3.5), σ = −pδ + τ . (4.17) Après remplacement de ce tenseur dans l’équation locale de bilan de la quantité de mouvement (4.16) on obtient : ρ Du = −∇p + ∇ · τ + ρg. Dt (4.18) Cette équation est souvent appelée équation de Navier-Stokes. 4.5 4.5.1 Bilan de l’énergie Premier principe de la thermodynamique L’énergie totale E (densité massique ǫ) contenue dans un notre volume de contrôle V est la somme de deux contributions, l’énergie cinétique liée au mouvement du fluide qu’il contient et l’énergie interne :  Z  Z u2 ρ e+ ρǫdV = E≡ dV, (4.19) 2 V V e étant la densité massique d’énergie interne. Le premier principe de la thermodynamique nous dit que la variation de cette quantité provient de l’apport de chaleur et du travail des forces appliquées. Par ailleurs, comme on l’a vu pour les équations de bilan obtenues précédemment, l’apport de chaleur contient une partie provenant des mouvements de matière (flux convectif) que l’on peut tout de suite séparer du reste : Z Z ∂ (4.20) (ρǫ) dV = − ρǫu · dS + Q + W, V ∂t S W étant la puissance des forces appliquées et Q la chaleur apporté au volume autrement que par mouvement. Ce dernier contient deux contributions : un flux de surface non–convectif (conductif et éventuellement radiatif) et une production interne, lié à la radioactivité naturelle : Z Z ρhdV, (4.21) Q = − q · dS + S V q étant la densité de flux de chaleur, généralement appelée plus simplement “flux de chaleur” et h le taux massique de chauffage interne par désintégration d’éléments radioactifs. L’intégrale de surface peut bien sûr être transformée en intégrale de volume, si bien que Z (ρh − ∇ · q) dV. (4.22) Q= V 40 Équations de bilan Les forces dont nous devons maintenant calculer le travail ont déjà été énumérées au moment de considérer la bilan de la quantité de mouvement (§ 4.4) et se séparent en deux groupes : les forces de volume et les forces de surface. Pour obtenir leur puissance associée, il suffit de les multiplier par la vitesse du fluide à l’endroit où elles s’appliquent et de les intégrer. En se basant sur la forme intégrale de l’équation de bilan de la quantité de mouvement (4.14) on obtient facilement W = = Z ZS V u · (σ·dS) + Z V ρu · gdV, (4.23) [∇ · (σ · u) + ρu · g] dV. Pour démontrer le théorème de Stokes dans le cas utilisé ci-dessus, le plus simple est d’utiliser celui-ci sous la forme d’indices : Z Z ∂i qi dV (4.24) qi dSi = S V qui, dans le cas qui nous intéresse s’écrit Z Z Z Z u · (σ·dS) = ui σij dSj = ∂j (σij ui )dV = ∂j (σji ui )dV, S S V (4.25) V où l’on a utilisé la symétrie du tenseur des contraintes (§ 3.4). Remplaçant Q et W donnés par les équations (4.22) et (4.23) dans l’équation de bilan (4.20), on obtient, après avoir regroupé tous les termes de variations de l’énergie dans le terme de gauche comme expliqué en § 4.3 : Z Z Dǫ ρ dV = [ρh − ∇ · q + ∇ · (σ · u) + ρu · g] dV. (4.26) V Dt V Cette équation de bilan étant valable pour tout volume de contrôle, on obtient une équation de bilan locale de bilan de l’énergie :   u2 D ρ e+ = ρh − ∇ · q + ∇ · (σ · u) + ρu · g. (4.27) Dt 2 À nouveau, on identifie la pression et la partie visqueuse dans le tenseur des contraintes pour obtenir la forme finale de l’équation de bilan de l’énergie totale :   u2 D e+ = ρh − ∇ · q − ∇ · (pu) + ∇ · (τ · u) + ρu · g. (4.28) ρ Dt 2 4.5.2 Énergie interne Une deuxième équation de bilan de l’énergie peut maintenant être obtenue dans laquelle l’énergie cinétique n’intervient pas. Pour cela, on utilise l’équation de bilan de la quantité 4.6 Bilan d’entropie de mouvement (4.18) que l’on multiplie par la vitesse,   D u2 = u · (∇ · τ − ∇p + ρg) . ρ Dt 2 41 (4.29) Cette équation peut alors être soustraite à l’équation de bilan de l’énergie totale (4.28) pour obtenir : ρ De = ρh − ∇ · q − p∇ · u + τ : ∇u, Dt (4.30) avec ∇u ≡  ∂ui ∂xj  (4.31) le tenseur des gradients de vitesse. Ce tenseur est d’ordre 2 et le symbole : signifie que l’on contracte les deux indices, c’est à dire τ : ∇u ≡ τij ∂uj . ∂xi (4.32) Là encore, la démonstration de l’équation (4.30) est obtenue le plus simplement en utilisant une notation en indices pour la partie concernant le tenseur des contraintes visqueuses, en se rappelant de sa symétrie : ∂i (τij uj ) − uj ∂i τij = τij ∂i uj . (4.33) On voit que l’équation (4.30) est une équation pour la bilan de l’énergie interne uniquement puisqu’elle ne fait plus apparaı̂tre l’énergie cinétique. Cette équation est indépendante de l’équation de bilan pour l’énergie totale puisqu’elle fait intervenir l’équation de bilan de la quantité de mouvement. Lorsque l’on compare cette équation pour l’énergie interne à l’équation générale de bilan (4.13) on peut identifier Ie = q, χe = ρh − p∇ · u + τ : ∇u. (4.34) (4.35) On voit donc que l’énergie interne n’est pas conservée (même en l’absence de radioactivité) et contient, comme source, le travail de forces de pression et de viscosité. Son seul flux est le flux de chaleur. 4.6 Bilan d’entropie En thermodynamique classique, on sait que l’énergie interne varie avec l’entropie et le volume comme dE = T dS − pdV. (4.36) 42 Équations de bilan Le volume spécifique n’étant pas très pertinent lorsqu’on ne travaille pas avec un gaz, nous utilisons plutôt la masse volumique, qui est son inverse. Par ailleurs, toute variation ici est liée au temps et à l’espace et l’équivalent en thermodynamique hors équilibre de la relation (4.36) est Ds p Dρ De =T + 2 . Dt Dt ρ Dt (4.37) On utilise cette relation pour remplacer la variation d’énergie interne dans l’équation (4.30), ce qui donne : ρT Ds p Dρ = ρh − ∇ · q + τ : ∇u − p∇ · u − . Dt ρ Dt L’équation de conservation de la masse (4.5) nous montre que les deux derniers termes dans le membre de droite se compensent et on obtient : ρT Ds = ρh − ∇ · q + τ : ∇u. Dt (4.38) Cette équation peut être modifiée pour la mettre sous la forme générale d’une équation de conservation (4.11) et isoler le flux et la production d’entropie : ρ 1 q 1 ρh Ds = −∇ · − 2 q · ∇T + τ : ∇u + , Dt T T T T (4.39) ce qui donne Is = q T χs = − (4.40) 1 ρh 1 q · ∇T + τ : ∇u + . 2 T T T (4.41) Le second principe de la thermodynamique nous dit qu’il existe une fonction d’état, appelée entropie, qui est telle qu’elle ne peut qu’augmenter pour un système isolé. Dans le cas d’un système non isolé, qui est le cas généralement considéré dans la thermodynamique hors équilibre et en particulier ici où nous nous intéressons aux échanges de chaleur dans la Terre, on peut écrire pour tout système que la variation de son entropie est due à deux contributions ds = δte s + δti s, dt (4.42) où l’indice e renvoie aux variations dues aux échanges avec l’environnement, alors que l’indice i renvoie aux variations d’origine interne. Le second principe de la thermodynamique dit que les variations d’origine interne sont forcément positives : δti s ≥ 0. (4.43) On a vu au paragraphe 4.3 que le terme d’échange dans une équation de conservation locale s’écrit comme −∇ · I′ , et que le reste est lié à la production interne. Le second principe peut donc s’écrire χs ≥ 0. (4.44) 4.6 Bilan d’entropie 43 Dans le cas considéré ici, il y a trois contributions à la production interne d’entropie. On peut trouver des situations où il seule une de ces contributions existe. Par exemple, en l’absence de production de chaleur par radioactivité et de mouvement, il ne reste que la contribution provenant de la conduction. Le second principe de la thermodynamique reste valide, ce qui signifie : − 1 q · ∇T ≥ 0. T2 (4.45) Le flux de chaleur conductif est lié au gradient de température et on considère souvent (en particulier ici) que cette relation est linéaire, q = −k∇T, (4.46) qui n’est autre que la loi de Fourier, k étant la conductivité thermique. La production d’entropie due à la conduction est alors χk = k (∇T /T )2 ≥ 0, ce qui montre que la conductivité thermique doit être positive. De même, Il peut exister des systèmes fluides en mouvement ayant une température uniforme, en l’absence de toute production de chaleur. Le second principe s’écrit alors : 1 τ : ∇u ≥ 0. T (4.47) La température étant bien sûr (comme dans toute thermodynamique) strictement positive, cela implique que le produit τ : ∇u soit positif. Comme discuté précédemment, on peut décomposer le tenseur gradient de vitesse en trois parties (eq.3.15), Gij ≡ ∂i uj = ωij + dij + tij , avec ωij la partie anti–symétrique, dij la partie symétrique à trace nulle et tij la partie contenant la trace et qui correspond à un changement de volume. De même que dans le cas de la conduction, on peut relier le tenseur des contraintes à celui des taux de déformation, par une rhéologie, comme discuté précédemment (§ 3.5). En supposant une rhéologie newtonienne isotrope, on a τij = 2η(dij ) + 3ζtij (eq.3.28). Il s’agit alors de calculer la contraction des trois tenseur ω, d et t en utilisant leurs propriétés. Vérifions d’abord que le produit d’un tenseur symétrique et d’un tenseur anti–symétrique est nul : τij ωij = τji ωij = τji (−ωji ) = −τij ωij = 0. (4.48) La première égalité vient de la symétrie de τ , la seconde de l’anti–symétrie de ω . Cela signifie que la partie rotation solide du tenseur gradient de vitesse ne contribue pas à la dissipation visqueuse, ce qui est logique puisqu’elle ne correspond pas à une déformation. Voyons maintenant le produit entre d et t : dij tij = ∇·u ∇·u dij δij = dii = 0, 3 3 puisque par définition la trace du tenseur d est nulle. On trouve ainsi  2 η X ∂ui ζ ∂ui τ : ∇u = 2ηdij dij + 3ζtij tij = + (∇·u)2 . + 2 i,j ∂xj ∂xj 3 (4.49) (4.50) Pour un fluide donné, on peut toujours trouver un écoulement tel l’un ou l’autre des termes de droites de l’équation (4.50) soit nul. Le second principe de la thermodynamique implique donc que η > 0 et ζ > 0. 44 Équations de bilan 4.7 Évolution de la température L’équation de bilan pour l’énergie interne est généralement plus utile que celle pour l’énergie totale car ne contient pas l’énergie cinétique. Pour rendre cet avantage plus net, il est intéressant de se ramener à l’évolution de la température, quantité plus pertinente que l’énergie interne lorsque l’on s’intéresse au transfert de chaleur. Ceci nécessite l’introduction d’une capacité calorifique. L’énergie interne peut être exprimée en fonction de la température et le volume massique V sous la forme différentielle suivante :     ∂e ∂e de = dV + dT, (4.51) ∂V T ∂T V qui est une relation classique de la thermodynamique à l’équilibre. Rappelons ici que, même si l’on s’intéresse à des problèmes hors équilibre, les relations locales sont obtenues en écrivant des bilans sur un volume élémentaire suffisamment petit pour pouvoir être considéré à l’équilibre interne, tout en étant en déséquilibre avec son environnement. Cela signifie que les relations de la thermodynamique à l’équilibre sont valable dans leur forme locale différentielle. L’équation (4.51) donne une relation entre les dérivées, totales ou partielles :     DV DT ∂e ∂e De = + . (4.52) Dt ∂V T Dt ∂T V Dt Par définition, la capacité calorifique à volume constant par unité de masse CV est   ∂e CV = , ∂T V (4.53) et une relation classique peut être obtenue pour l’autre dérivée partielle (à développer), si bien que     DV DT ∂p De = −p + T ρ + ρCV . (4.54) ρ Dt ∂T V Dt Dt Pour les phases condensées (liquides et solides), préfère en général utiliser la masse volumique ρ que la volume massique V = 1/ρ et l’on a ρ DV 1 Dρ =− = ∇ · u. Dt ρ Dt On obtient alors l’équation suivante pour l’évolution de la température   DT ∂p ρCV = ρh − ∇ · q − T ∇ · u + τ : ∇u. Dt ∂T V (4.55) (4.56) On préferre souvent travailler avec p et T plutôt que V et T comme variables d’état lorsque l’on étudie des systèmes condensés, liquides ou solides. Dans ce cas là, l’énergie interne s’écrit comme (à développer ou mettre en appendice) DT αT Dp ∇·u De = Cp − −p , Dt Dt ρ Dt ρ (4.57) 4.7 Évolution de la température 45 et l’équation d’évolution de la température est alors ρCp DT Dp = ρh − ∇ · q + αT + τ : ∇u. Dt Dt (4.58) L’équation que l’on vient d’obtenir est générale mais d’autre équations peuvent être obtenues pour des cas particuliers. Celui qui nous intéresse le plus est celui d’un solide indéformable et en particulier dont la densité ne dépend pas de la température. Dans ce cas là, Cp = CV et u = 0. Si de plus on suppose que la loi de Fourier est valide, et que la conductivité thermique est constante, on obtient ρCp ∂T = k∇2 T + ρh. ∂t (4.59) Il s’agit de l’équation de la conduction de la chaleur. L’utilisation de Cp est préferrable à celle de CV pour les matière condensées car c’est un paramètre plus facilement accessible à la mesure dans ce cas là. On définit κ = k/ρCp la diffusivité thermique et on obtient l’équation de diffusion suivante h ∂T = κ∇2 T + . ∂t Cp (4.60) Cette équation a le mérite de faire ressortir immédiatement l’échelle de temps d’évolution de la température. En effet, l’analyse dimentionnelle des deux premiers termes de cette équation nous montre que l’échelle de temps d’évolution de la température à une distance d d’une source de chaleur est d2 /κ. 46 Équations de bilan Chapitre 5 Géophysique des fonds océaniques 5.1 Observations générales Profondeur (m) 7000 3500 1500 500 0 -500 -1500 -2500 -3500 -4500 -5500 -6500 Fig. 5.1 – Topographie des continents et des océans. On voit bien les dorsales océaniques qui forment une chaı̂ne de montagnes sous-marine. En s’en éloignant perpendiculairement, la profondeur augmente. Carte et données disponible à http ://topex.ucsd.edu/marine topo. Les observations des fonds océaniques ont permis de démontrer l’ouverture océanique et la tectonique des plaques. C’est dans ce contexte que de nombreuses observations des fonds océaniques ont été depuis interprétées. En particulier, la topographie des fonds océaniques (fig. 5.1) montre que les dorsales représentent des montagnes sous-marines et que lorsque l’on s’en éloigne la profondeur augmente, de manière quasiment monotone. Cette observation semble simple à expliquer et cela fait partie des objectifs de ce chapitre. 48 Géophysique des fonds océaniques 40 60 80 100 mW 120 220 320 m-2 Fig. 5.2 – Modèle en harmoniques sphériques du flux de chaleur terrestre jusqu’au degré 18. Tiré de Pollack et al (1993). Une autre observation concerne le flux de chaleur terrestre qui est mesuré dans les fonds océaniques aussi bien que sur les continents depuis des décennies. La figure (5.2) représente un modèles en harmoniques sphériques au degrés 18 de ce flux de chaleur et montre que les dorsales sont associées à un flux élevé et que le flux diminue à mesure que l’on s’en éloigne. Là aussi, l’observation est suffisamment simple pour que l’on puisse espérer l’expliquer avec un modèle géophysique. Le fait que les deux observations a priori indépendantes du flux de chaleur et de la topographie sont en fait corrélées nous indique que le même modèle doit être capable d’expliquer les deux observations. Nous allons donc d’abord développer un modèle basé sur l’évolution thermique des fonds océaniques après leur formation aux dorsales. Ce modèle thermique nous donnera directement accès au flux de chaleur et on pourra alors comparer ses prédictions aux données existantes. La variation de densité avec la température, moteur de la convection thermique et donc de la tectonique des plaques, nous fournira un mécanisme physique pour l’approfondissement des fonds océanique lorsque l’on s’éloigne des dorsales. La comparaison avec les observations nous permettra là encore de tester le modèle et de contraindre ses paramètres libres. 5.2 Un modèle de refroidissement de la lithosphère 5.2 49 Un modèle de refroidissement de la lithosphère La lithosphère, par définition, est la zone à la surface de la Terre qui se comporte de manière rigide (de lithos, pierre en grec). La raison pour cette rigidité, par opposition à la faible rigidité de l’asthenosphère, provient de sa faible température. En effet, la viscosité des roches du manteau, à l’instar de celle de beaucoup de matériaux comme le miel, dépend fortement de la température, et celle–ci est la plus faible proche de la surface. C’est ce qui confère cette rigidité à la lithosphère. Pour comprendre les fonds océaniques il est donc nécessaire de s’intéresser à la température. Plus généralement, on considère maintenant que la lithosphère océanique est la couche limite supérieure de la convection dans le manteau, pour laquelle les variations de température sont le moteur. u x=u·a δ(a) T=1300ºC Fig. 5.3 – Modèle schématique de lithosphère océanique. À partir de la dorsale où la matière arrive chaude vers la surface (ce qui induit de la fusion partielle par décompression et produit la croûte océanique), un refroidissement se produit par contact avec l’eau de mer. La zone froide est très visqueuse et forme la lithosphère. La diffusion thermique permet au front froid de se propager vers le bas, ce qui mène à l’épaississement de la lithosphère. Pour calculer la température dans la lithosphère océanique, le principe est expliqué sur la figure 5.3. La croûte océanique est formée par fusion due à la décompression de la matière mantellique chaude qui arrive proche de la surface au niveau de la dorsale. Au contact avec l’eau de mer, la matière se refroidit en même temps qu’elle se déplace horizontalement. Le front froid (matérialisé sur le schéma par l’isotherme 1300◦ C, on verra plus loin pourquoi), inexistant proche de la dorsale où le volcanisme est actif, se propage vers le bas lorsque l’on s’éloigne de la dorsale. Cet enfoncement se produit par diffusion thermique et dépend du temps écoulé depuis la formation de la lithosphère à la dorsale. Ce problème est similaire à celui du refroidissement de la Terre tel que formulé par Fourier et Kelvin et sa solution est bien connue (§ 2.2). Pour faire le lien entre les deux problèmes, écrivons l’équation de diffusion dans sa version la plus simple DT = κ∇2 T Dt (5.1) où, par rapport à l’équation générale (4.58), on a négligé la radioactivité, l’effet de la variation de pression et la dissipation visqueuse et on a utilisé la loi de Fourier en supposant une conductivité thermique constante. La concentration en éléments radioactifs (U, Th, K) de la lithosphère 50 Géophysique des fonds océaniques océanique est suffisamment faible, et la lithosphère suffisamment fine pour que la radioactivité puisse être négligée. La faible épaisseur de la lithosphère, en particulier par rapport à son extension horizontale, permet une autre simplification : la dérivée horizontale de la température ou tout autre quantité est faible devant la dérivée verticale. De ce fait, on ne conserve dans le laplacien du terme de droite que le terme en dérivée verticale. Pour ce qui est du terme de gauche, DT /Dt est la variation de la température en suivant la lithosphère (dérivée entraı̂née) et peut donc être remplacée par la dérivée en fonction de l’âge de la lithosphère, ce qui donne à une dimension : ∂2T ∂T =κ 2, ∂a ∂z (5.2) avec a l’âge de la lithosphère (c’est à dire le temps écoulé depuis sa création au niveau de la dorsale), z la profondeur. En analysant les dimensions de√l’équation (5.2), on voit que la profondeur δ du front froid ne peut qu’être proportionnelle à κa. En prenant l’âge maximum de la lithosphère océanique, 180 Ma, et une diffusivité thermique κ = 10−6 m2 /s, on trouve une épaisseur maximale de lithosphère océanique de l’ordre de δmax ≃ 75 km. Cette épaisseur est faible par rapport à l’épaisseur du manteau et on considère donc un milieu semi–infini. Pour résoudre le problème, il faut ajouter une condition initiale et des conditions limites. Au niveau de la dorsale (a = 0), la température est considérée comme uniforme, égale à TM . Cette température est aussi la température pour z → ∞. À la surface, le contact avec l’océan maintient la température à une valeur constante de l’ordre de 2◦ C (que l’on prendra égale à 0, pour simplifier). La solution à ce problème est connue depuis Fourier (et a été obtenue plus haut § 2.2) :   z √ T = TM erf . (5.3) 2 κa Le flux de chaleur en surface est obtenu en prenant le gradient vertical à la surface, que l’on multiplie par la conductivité thermique, k, (Loi de Fourier) : −q = k kTM ∂T (0, a) = √ . ∂z πκa (5.4) Nous nous intéressons ici à −q, qui est positif et est le flux dans la direction des z décroissants, c’est à dire vers l’extérieur de la Terre. On a ainsi obtenu une expression du flux de chaleur océanique en fonction de l’âge de la lithosphère. 5.3 Âge des fonds océaniques Pour comparer la prédiction de l’équation (5.4) avec les observations, il faut mesurer le flux de chaleur et l’âge de la lithosphère aux mêmes points. La détermination de l’âge d’une roche utilise souvent des méthodes radiométriques, mais cela n’est guère pratique pour les fonds océaniques. On utilise en fait le magnétisme et la compréhension que l’on a de l’expansion océanique. On se rappelle en effet que c’est la découverte des bandes de polarité magnétique successivement normales et inverses enregistrées dans la lithosphère océanique qui a permis de démontrer sans équivoque les inversions du champ magnétique et l’ouverture océanique. La succession irrégulière des inversions du champ magnétique permet en outre de dater les fonds 5.3 Âge des fonds océaniques 51 Fig. 5.4 – Le code barre des fonds océaniques. À gauche, les âges des transitions entre périodes de polarité normale (c’est à dire identique au temps présent, en noir) et les périodes inverses (en blanc) donnés par Cande & Kent (1995). La succession de ces périodes se retrouve dans l’aimantation du plancher océanique (droite, d’après Vine, 1966) et on peut ainsi dater les fonds océaniques en faisant correspondre les limites de polarité observées (haut droite) avec celles données dans l’échelle des âges (gauche). 52 Géophysique des fonds océaniques océaniques. Les datations radiochronologiques des inversions du champ magnétiques ont permis de construire une échelle des temps, véritable “code–barre” qui par comparaison avec le champ magnétique crustal du plancher océanique (fig. 5.4), permet de dater ces fonds. La figure 5.5 montre le résultat sous forme de carte d’un tel travail. 0.0 20.1 40.1 55.9 83.5 126.7 154.3 Ma 180.0 Fig. 5.5 – Carte des âges (en millions d’années, Ma) des fonds océaniques obtenus par le magnétisme. Données tirées de Müller et al (2008) et disponible à http ://www.earthbyte.org/Resources/agegrid2008.html. 5.4 Mesures du flux de chaleur océanique et comparaisons avec le modèle Flux de chaleur (mW m-2) (Stein & Stein, 1992) Fig. 5.6 – Compilation des données de flux de chaleur océanique, rangées en fonction de l’âge du plancher océanique. Les points représentent la moyenne pour chaque intervalle de deux millions d’années et les lignes représentent l’écart type de la distribution. Tiré de Stein & Stein (1992) Âge (Ma) Maintenant que l’on connaı̂t l’âge des fonds océaniques, il nous faut déterminer le flux de 5.4 Mesures du flux de chaleur océanique et comparaisons avec le modèle 53 chaleur. Celui–ci est mesuré lors de campagnes océanographiques dédiées qui utilisent une sonde équipée de thermocouples et de résistances chauffantes que l’on enfonce dans les sédiments peu consolidés du plancher océanique. Cette sonde permet de mesurer au même point la conductivité thermique et le gradient vertical de température. (Ajouter un paragraphe sur la technique de mesure.) Le flux de chaleur est alors obtenu en utilisant la loi de Fourier (première partie de l’équation (5.4)). On peut alors collecter toutes les données obtenues et les classer par ordre d’âge du plancher océanique. La figure 5.6 représente le résultat d’une telle compilation et semble indiquer un accord approximatif avec la théorie mais la dispersion des données est trop importante pour trancher. Le problème est de deux ordres ici : les données anciennes étaient obtenues en ne mesurant que le gradient de température, la conductivité thermique étant estimée, et les mesures peuvent être faussées par la présence de circulation hydrothermale. Les chercheurs qui travaillent sur ce sujet on développé une méthode pour déterminer la conductivité thermique du milieu directement au point de mesure en mesurant sa réponse au chauffage. Ils ont ainsi montré que cette conductivité thermique était fortement variable dans les sédiments océaniques, du fait notamment des variation de la porosité. Il est donc crucial de ne retenir que les données de flux de chaleur pour lesquelles la conductivité thermique a été mesurée en même temps que la température. Pour ce qui concerne le second problème, celui de la circulation hydrothermale, il est nécessaire de vérifier la cohérence spatiale des données et la conservation de l’énergie pour chaque point de mesure. Théorie 300 200 100 0 (Lister etal, 1990) 0 40 80 120 Âge (Ma) 160 200 Flux de chaleur (mW m-2) Flux de chaleur (mW m-2) 400 60 50 40 30 100 120 140 160 180 200 Âge (Ma) Fig. 5.7 – Compilation des données de flux chaleur océanique en fonction de l’âge de la croûte. Chaque boite représente la variabilité des données dans la tranche d’âge considérée. La ligne pointillée est un ajustement selon l’équation (5.4). La figure de droite est un zoom pour les âges élevés et montre un écart systématique avec le modèle. Tirée de Lister et al (1990). Lister et al (1990) ont trié les données qui, dans la base globale, satisfaisaient à ces critères de qualité ont obtenu la figure 5.7. L’accord entre la théorie donnée par l’équation (5.4) et les observations est très bon pour les âges inférieurs à 80 Ma, âge à partir du quel le flux de chaleur semble cesser sa décroissance. Plusieurs paramètres libres dans l’équation nécessitent un ajustement, et en particulier la température, et le résultat est la ligne pointillée de la figure. Cela permet de contraindre la température du manteau sous la lithosphère qui doit être d’environ 1300◦ C. Une vérification directe de la théorie peut également être faite lors d’une campagne de mesures. La figure 5.8 montre les résultats d’une telle campagne, effectuée à proximité de la dorsale Juan de Fuca. On voit que les mesures suivent très bien la théorie, représentée avec son incertitude par les deux lignes pointillées, là où la couverture sédimentaire est suffi- Géophysique des fonds océaniques Flux de chaleur (mW m-2) 54 600 Théorie 400 200 Mesures 0 Profondeur (m) 2000 Océan 2400 Sédiments 2800 Socle 3200 (Davis et al 1999) 0 20 40 60 80 100 Distance à la dorsale (km) Fig. 5.8 – Flux de chaleur (haut) et profil topographique et épaisseur sédimentaire (bas) sur le flanc est de la dorsale Juan de Fuca, en fonction de la distance à la dorsale. Tiré de Davis et al (1999). sante pour que les mesures soient de bonne qualité. Par contre, on voit que lorsque l’épaisseur sédimentaire devient trop faible, le flux mesuré décroche et devient très faible. On met ici en évidence l’importance du flux convectif lié à la circulation hydrothermale. En effet, lorsque la couche sédimentaire est peu épaisse, elle ne se compacte que faiblement et sa perméabilité reste importante. Proche de la dorsale, là ou le flux de chaleur est important, la convection dans le milieu sédimentaire poreux peut démarrer. Dans ce cas là, le transport de chaleur se fait majoritairement par convection et plus par conduction. La mesure du flux conductif ne peut alors fournir qu’une borne inférieure du flux total. On peut aussi comparer le modèle de refroidissement de la lithosphère océanique aux observations sismologiques, ou mieux aux modèles tomographiques qui en sont tirés. La figure 5.9 montre les variations de la vitesse VS des ondes S dans le manteau sous l’océan Pacifique en fonction de la profondeur et de l’âge de la plaque. On retrouve très clairement la structure en température calculée précédemment jusqu’à la profondeur de 200 km environ. Au delà, VS augment avec la profondeur du fait de l’augmentation de pression. Enfin, on peut comparer cette théorie et ces observations à des résultats numériques de modèles de convection. La figure 5.10 montre de tels résultats (Grigné et al, 2005; Labrosse & Jaupart, 2007). Le calcul est ici réalisé en deux dimensions, c’est à dire qu’il représente une tranche du système complet. La partie basse de la figure représente la température (en unités arbitraire entre 0, en bleu, et 1, en rouge) et montre bien, en bleu, une lithosphère froide qui s’épaissit en s’éloignant de la dorsale, et qui finit par subducter. La courbe donnant le flux de chaleur en surface (qtop ) montre bien une diminution lorsque l’on s’éloigne de la dorsale. La courbe donnant l’inverse du carré de ce flux de chaleur (1/q 2 ) montre une évolution linéaire avec la distance à la dorsale. Ceci indique que la vitesse de surface ne varie pas tellement et que la théorie développée pour les fonds océaniques permet également d’interpréter les résultats de 5.4 Mesures du flux de chaleur océanique et comparaisons avec le modèle 55 (Priestley & McKenzie, 2006) Fig. 5.9 – Variations des vitesses sismiques d’ondes S dans le manteau de l’océan Pacifique en fonction de l’âge. Figure tirée de Priestley & MKenzie (2006) 1/q2 Subduction Dorsale Subduction Dorsale Subduction Dorsale Subduction Dorsale 0.1 qtop 0.0 Flux de chaleur à la surface 50 utop 0 Vitesse de surface 0 Température 1 2 4 6 (Grigné et al, 2005 ; Labrosse & Jaupart, 2007) 8 10 12 14 16 18 20 22 24 26 0.2 28 0.4 30 0.6 32 0.8 Fig. 5.10 – Résultat d’un modèle de convection dans le manteau (Grigné et al, 2005). La partie basse représente la température dans le domaine de calcul (en unités arbitraire entre 0, en bleu, et 1, en rouge) et les courbes supérieures montrent la vitesse de surface (utop ), le flux de chaleur à la surface (qtop ) et l’inverse de son carré (1/q 2 ). 56 Géophysique des fonds océaniques modèles numériques de convection. 5.5 Modèle d’approfondissement du plancher et comparaison avec les observations Une autre observation importante concerne la profondeur du plancher océanique, qui augmente avec la distance à la dorsale océanique (figure 5.1). Comme on l’a vu pour le flux de chaleur, la distance à la dorsale ne semble pas être la variable la plus pertinente et il semble plus intéressant de considérer la variation de la profondeur en fonction de l’âge du plancher océanique, ou mieux encore, la racine de cet âge. La figure 5.11 montre justement une compilation de telles données, (Carlson & Johnson, 1994) Root(Age, Ma) Fig. 5.11 – Profondeur du plancher oceánique en fonction de la racine de son âge. Tiré de Carlson & Johnson (1994). océan par océan et en moyenne globale. On voit qu’effectivement la profondeur augmente linéairement avec la racine carrée de l’âge, jusqu’à un âge d’environ 80 millions d’années où la profondeur n’augmente plus. On peut d’ailleurs remarquer que cela correspond également à l’arrêt de la variation du flux de chaleur avec l’âge sur la figure 5.7. Cet écart entre théorie et observation est généralement interprété comme étant dû au démarrage d’une convection à petite échelle sous la lithosphère mais cela de sera pas discuté ici. Pour la plus grande partie du plancher océanique, l’approfondissement proportionnel à la racine de l’âge est observé et nécessite une explication. La raison pour laquelle le plancher océanique voit sa profondeur augmenter avec l’âge est liée à son évolution thermique, telle que décrite plus haut. Lorsque la température baisse, la matière généralement se contracte, où encore sa masse volumique (ou sa densité) augmente. Il s’agit du 5.5 Profondeur des fonds océaniques 57 mécanisme moteur pour la convection thermique : lorsque la matière refroidit à la surface de la Terre, elle devient plus dense que celle du manteau et elle finit par plonger. Ce phénomène de contraction thermique (ou de dilatation dans le cas de l’opération inverse) est quantifié par un paramètre thermodynamique appelé coefficient de dilatation thermique, α. En utilisant la température du manteau, TM , comme référence, on écrit la variation de masse volumique avec la température comme ρ = ρM [1 − α (T − TM )] . (5.5) En se refroidissant, la plaque lithosphérique voit sa densité augmenter, ce qui induit une déformation du manteau sous-jacent. Le plancher s’enfonce de façon à ce que la masse d’une colonne verticale soit identique en tout point. C’est ce qu’on appelle l’isostasie. On note d(a) la différence de profondeur de l’océan entre la dorsale (âge nul) et celle à l’âge a (fig. 5.12), et ρw la masse volumique de l’eau. u x=u a ρw d ρ(Τ) ρ M zc Fig. 5.12 – Modèle schématique de lithosphère océanique prenant en compte, de façon exagérée, l’enfoncement du plancher océanique. Le niveau horizontal rouge, situé sous la lithosphère, est le niveau de compensation. L’équilibre isostatique postule que la pression à une profondeur dı̂te de compensation zc est indépendante de l’âge. En supposant que cette pression est hydrostatique, elle est égale au poids de la matière au dessus. En écrivant cette égalité pour les deux points figurés en rouge sur la figure 5.12, on obtient ρM (d + zc ) = ρw d + Z 0 zc ρdz ⇒ d = ρM 1 − ρw Z 0 zc (ρ − ρM )dz. (5.6) L’écart entre la densité en profondeur (ρM ) et celle (ρ) qui dépend de la température dans la lithosphère est obtenue en utilisant l’équation (5.5), dans laquelle la température donnée par 58 Géophysique des fonds océaniques l’équation (5.3) est utilisée : Z zc ρM α (TM − T ) dz, d= ρ M − ρw 0   Z zc  z ρM αTM √ dz, 1 − erf ρ M − ρw 0 2 λa √ √ Z 2ρM αTM λa zc /2 λa [1 − erf (u)] du, ρM − ρw 0 (5.7) où l’on a fait un changement de variable pour obtenir la dernière égalité. Comme la variation de température, et donc de masse volumique, n’est importante que dans la lithosphère, c’est à dire près de la surface, toute valeur de zc supérieure à l’épaisseur de la lithosphère donne un résultat similaire et on peut calculer l’intégrale en prenant zc → ∞. L’intégrale faisant intervenir la fonction erreur est connue : Z ∞ 1 (1 − erf x)dx = √ , (5.8) π 0 et on obtient 2ρM αTM d= ρM − ρw r λa . π (5.9) On explique ainsi l’augmentation de la profondeur comme la racine carrée de l’âge du fond océanique. 5.6 Carte de flux calculé et flux océanique total 48 60 80 100 120 140 220 320 mWm−2 5000 Fig. 5.13 – Flux de chaleur océanique calculé à partir de la carte des âges du plancher selon l’équation (5.10). 5.6 Carte de flux calculé et flux océanique total 59 Ayant obtenu un bon accord entre théorie et observations du flux de chaleur, on peut utiliser la carte des âges des fonds océaniques (figure 5.5) pour calculer une carte de flux de chaleur, avec la formule 490 q= √ , a (5.10) avec q le flux en mW m−2 et a l’âge en millions d’années (Ma). Le coefficient, CQ = 490 ± 20 provient de l’ajustement fait sur la figure 5.7. On obtient alors la figure 5.13 où l’on voit que le flux de chaleur sortant au niveau des dorsales peut être très important, comme cela est directement mesuré lorsque la couverture sédimentaire est suffisante (figure 5.8). dA/dt (km2/yr) 3 Fig. 5.14 – Densité de probabilité de trouver une lithosphère océanique d’âge donné en fonction de l’âge. Le trait plein est la valeur mesurée avec la barre d’erreur en grisé, et le trait discontinu est le meilleur ajustement linéaire. Tiré de Cogné & Humler (2004). 2 1 0 200 150 100 50 0 Age (Ma) À partir de cette carte, on peut calculer le flux de chaleur total extrait de la Terre au travers des fonds océaniques, dont la surface totale est Aoc . Pour cela on range les éléments de surface en fonction de leur âge, de 0 à la valeur maximale observée amax , et on calcule la somme : Z amax Z dA q(a) (a)da qdA = (5.11) Qoc = da 0 Aoc où l’on a introduit la densité de probabilité dA/da(a) de trouver une lithosphère d’âge a. Cette distribution est obtenue en mesurant la surface occupée par chaque tranche d’âge et il a été observé (Parsons, 1982; Cogné & Humler, 2004) que cette distribution avait une forme grossièrement triangulaire (fig. 5.14), c’est à dire   a dA , = C0 1 − da amax (5.12) avec C0 le taux actuel de création de lithosphère océanique et amax l’âge maximum observé. Ces deux paramètres sont en fait reliés par la surface océanique totale : Z amax dA C0 amax Aoc = (a)da = . (5.13) da 2 0 60 Géophysique des fonds océaniques En injectant ces deux dernières relations dans l’équation (5.11), et en utilisant l’expression du flux de chaleur en fonction de l’âge (eq. 5.4) on obtient la perte de chaleur totale au travers des fonds océaniques : 8 4kTM √ kTM Qoc = √ C0 amax = Aoc √ . 3 πκamax 3 πκ (5.14) Pour obtenir une valeur numérique, il faut préciser les valeurs des différents paramètres qui interviennent dans cette équation. Les paramètres physique peuvent être déterminés en laboratoire si l’on dispose d’échantillons représentatifs mais il peut toujours subsister des doutes à ce sujet. Le moyen le plus simple et le plus juste est d’utiliser les données existantes dont on a √ montré qu’elles étaient bien expliquées par le modèle. Cela revient à considérer que kTM / πκ = 490Wm−2 Ma1/2 et on obtient Qoc = 30TW, c’est à dire 30 1012 W. Bibliographie Barenblatt, G. I. (1996). Scaling, Self-Similarity, and Intermediate Asymptotics. Cambridge University Press, Cambridge. 2.2 Bénard, H. (1900a). Les tourbillons celullaires dans une nappe liquide. deuxième partie : procédés mécaniques et optiques d’examen. lois numériques des phénomènes. Revue générale des Sciences, 12, 1309–1328. 7.1.1, 7.1.2 Bénard, H. (1900b). Les tourbillons celullaires dans une nappe liquide. première partie : description générale des phénomènes. Revue générale des Sciences, 12, 1261–1271. 7.1.1, 7.1.2 Bénard, H. (1901). Les tourbillons celullaires dans une nappe liquide transportant la chaleur par convection en régime permanent. Ann. Chim. Phys., 23, 62–144. 7.1.1, 7.1.2 Bénard, H. & Avsec, D. (1938). Travaux récents sur les tourbillons cellulaires et les tourbillons en bandes. applications à l’astrophysique et à la météorologie. J. Phys. Radium, 9, 486–500. URL http://dx.doi.org/10.1051/jphysrad:01938009011048600 7.1 Block, M. J. (1956). Nature, 178, 650–651. 7.1.1, 7.1.2 Burchfield, J. D. (1975). Kelvin and the age of the Earth. The University of Chicago Press. 2.2 Cande, S. C. & Kent, D. V. (1995). Revised calibration of the geomagnetic polarity timescale for the late cretaceous and cenozoic. JGR, 100, 6093–6095. 5.4 Carlson, R. L. & Johnson, H. P. (1994). On modeling the thermal evolution of the oceanic upper-mantle : An assessment of the cooling plate model. J. Geophys. Res., 99, 3201–3214. 5.11 Carnot, S. (1824). Réflexions sur la puissance motrice du feu et sur les machines propres à développer cette puissance. Bachelier, Paris. Réédité par les éditions Jacques Gabay, 1990, et disponible sur http ://www.gallica.fr. (document) Carslaw, H. & Jaeger, J. (1959). Conduction of heat in solids. Oxford university press. 2.1 Cogné, J.-P. & Humler, E. (2004). Temporal variation of oceanic spreading and crustal production rates during the last 180 My. Earth Planet. Sci. Lett., 227, 427–439. 5.14, 5.6 Davis, E. E., Chapman, D. S., Wand, K., Villinger, H., Fisher, A. T., Robinson, S. W., Grigel, J., Pribnow, D., Stein, J. & Becker, K. (1999). Regional heat flow variations across the sedimented juan de fuca ridge eastern flank : constraints on lithospheric cooling and lateral hydrothermal heat transport. J. Geophys. Res., 104, 17,675–17,688. 5.8 62 BIBLIOGRAPHIE Grigné, C., Labrosse, S. & Tackley, P. J. (2005). Convective heat transfer as a function of wavelength : Implications for the cooling of the Earth. J. Geophys. Res., 110, B03409. Doi :10.1029/2004JB003376. 5.4, 5.10, 7.2.3 de Groot, S. R. & Mazur, P. (1983). Non-equilibrium thermodynamics. Dover. 4.1, 4.1 Harrison, T. M. (1987). Comment on “Kelvin and the age of the earth”. Jour. Geology, 95, 725–727. 2.2 Heaviside, O. (1899). Electromagnetic theory, vol. 2, Chapitre 5 : Mathematics and the age of the Earth. SPON, London. 2.2 Howard, L. N. (1964). Convection at high Rayleigh number. In Proceedings of the Eleventh International Congress of Applied Mechanics (H. Gortler, ed.), 1109–1115. Springer-Verlag, New York. 7.2.2 Jochum, K. P., Hoffmann, A. W., Ito, E., Seufert, H. M. & White, W. M. (1983). K, U and Th in mid-ocean ridge basalt glasses and heat production. K/U and K/Rb in the mantle. Nature, 306, 431–436. 8.1 Kondepudi, D. & Prigogine, I. (1998). Modern thermodynamics. From heat engines to dissipative structures. Wiley, Chichester. 4.1 Krishnamurti, R. (1968). Finite amplitude convection with changing mean temperature. part 1. theory. J. Fluid Mech., 33, 445–455. 7.1.1 Krishnamurti, R. (1973). Some further studies on the transition to turbulent convection. J. Fluid Mech., 60, 285–303. 7.6 Labrosse, S. & Jaupart, C. (2007). Thermal evolution of the earth : Secular changes and fluctuations of plate characteristics. Earth Planet. Sci. Lett., 260, 465–481. 5.4 Landau, L. & Lifchitz, E. (1994). Physique théorique - Mécanique des fluides. Mir-Ellipses, 3e edn. 4.1 Lister, C. R. B., Sclater, J. G., Davis, E. E., Villinger, H. & Nagihara, S. (1990). Heat flow maintained in ocean basins of great age : investigations in the north-equatorial west pacific. GJI, 102, 603–628. URL http://dx.doi.org/10.1111/j.1365-246X.1990.tb04586.x 5.7, 5.4 McKenzie, D. & Weiss, N. O. (1975). Speculations on the thermal and tectonic history of the Earth. Geophys. J. R. Astr. Soc., 42, 131–174. 8.2 Müller, R. D., Sdrolias, M., Gaina, C. & Roest, W. R. (2008). Age, spreading rates, and spreading asymmetry of the world’s ocean crust. G3, 9. 5.5 Parsons, B. (1982). Causes and consequences of the relation between area and age of the ocean floor. J. Geophys. Res., 87, 289–302. 5.6 Perry, J. (1895a). The age of the Earth. Nature, 51, 582–585. 1.5, 2.2, 2.3 Perry, J. (1895b). On the age of the Earth. Nature, 51, 224–227. 1.5, 2.2, 2.3 Perry, J. (1895c). On the age of the Earth. Nature, 51, 341–342. 1.5, 2.2, 2.3 BIBLIOGRAPHIE 63 Pollack, H. N., Hurter, S. J. & Johnson, J. R. (1993). Heat flow from the earth’s interior : Analysis of the global data set. Rev. Geophys., 31. URL http://dx.doi.org/10.1029/93RG01249 5.2 Priestley, K. & MKenzie, D. (2006). The thermal structure of the lithosphere from shear wave velocities. Earth and Planetary Science Letters, 244, 285 – 301. URL http://www.sciencedirect.com/science/article/B6V61-4JFHFB5-1/2/ 127961ae059cfe9be839759f8df382de 5.9 Rayleigh, L. (1916). On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. Phil. Mag., 32, 529–546. 7.1, 7.1.2 Richter, F. M. (1986). Kelvin and the age of the earth. Jour. Geology, 94, 395–401. 2.2 Schubert, G. (1979). Subsolidus convection in the mantles of terrestrial planets. Ann. Rev. Earth Planet. Sci., 7, 289–342. 8.2 Schubert, G. & Young, R. E. (1976). Cooling the Earth by whole mantle subsolidus convection : A constraint on the viscosity of the lower mantle. Tectonophysics, 35, 201–214. 8.2 Sharpe, H. N. & Peltier, W. R. (1978). Parameterized mantle convection and the Earth’s thermal history. Geophys. Res. Lett., 5, 737–740. 8.2 Sotin, C. & Labrosse, S. (1999). Three-dimensional thermal convection of an isoviscous, infinite– Prandtl–number fluid heated from within and from below : applications to heat transfer in planetary mantles. Phys. Earth Planet. Inter., 112, 171–190. 7.4 Stein, C. A. & Stein, S. (1992). A model for the global variation in oceanic depth and heat flow with lithospheric age. Nature, 359, 123–129. URL http://dx.doi.org/10.1038/359123a0 5.6 Stevenson, D., Spohn, T. & Schubert, G. (1983). Magnetism and thermal evolution of the terrestrial planets. Icarus, 54, 466–489. 8.2 Thomson, W. L. K. (1862). On the secular cooling of the earth. Trans. Roy. Soc. Edinburgh, 23, 295–311. 1.5, 2.2, 2.1, 2.2 Thomson, W. L. K. (1899). The age of the earth as an abode fitted for life. Jour. Victoria Instit., 31, 11–35. 1.5 Turcotte, D. L. & Schubert, G. (2001). Geodynamics. Cambridge University Press, Cambridge, UK, 2e edn. 7.2.3 Index des sujets asthenosphère, 49 gradient de vitesse de déformation, 27 bilan de l’énergie, 39–41 de l’entropie, 41–43 de la quantité de mouvement, 38–39, 41 forme générale, 37–38 bilans, 35–45 Boltzmann, constante de, 5 Boussinesq, approximation de, 72, 75 hydrostatique, 72 capacité calorifique, 9, 44 à volume constant, 44 chaleur interne production de, 39 conditions limites, 24 conductivité thermique, 5, 6, 43 conservation de la masse, 36–37 contrainte, 24, 38 contrainte visqueuse, 33, 39 contraintes tenseur des, 38 convection dans le manteau, 49 couche limite, 49, 77, 79 Couette, écoulement de, 24, 80 manteau supérieur, 88 Marangoni, effet, 69 masse volumique, 42 matériel, volume, 37 incompatible, éléments, 88 incompressibilité, 37 isotrope, 7, 33 isotropes, 29 Levi-Civita, Symbole de, 29 lithosphère, 49 Navier-Stokes équation de, 39 Nusselt, nombre de, 76, 83 parcelle de fluide, 35 phonons, 5 Prandtl, nombre de, 73, 76 dérivée matérielle, 36 dérivée particulaire, voir dérivée matérielle dérivée totale, voir dérivée matérielle diffusivité thermique, 10, 45 divergence, 7 Rayleigh, nombre de, 68, 70 de couche limite, 79 Rayleigh-Bénard, convection de, 69 Reynolds théorème de transport de, 37 rhéologie, 23, 43 rhéologie Newtonienne, 25 rotation solide, 29 énergie cinétique, 39 interne, 39 similitude, variable de, 14 stabilité, 68 Stokes, vitesse de, 69 flux de chaleur, 39 Fourier, loi de, 5–7, 10, 43, 50 fusion partielle, 88 tenseur, 26 tenseur des contraintes, 30, 31 thermodynamique classique, 35, 41 hors équilibre, 35, 42 gradient de vitesse, 26, 29, 41 INDEX DES SUJETS premier principe, 39 second principe, 42–43 tourbillon, 29 travail des forces, 39 Urey, nombre d’, 87 viscosité cinématique, 25 viscosité dynamique, 24, 33 viscosité seconde, 33 vitesse de déformation, 26 volume de contrôle, 36 volume spécifique, 42 vorticité, 29 Watt, 6 65