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