Retour

L'équation de la chaleur

Un peu d'histoire

Une première chose : la dénomination "équation de la chaleur" pour désigner le phénomène et l'équation que nous allons aborder dans cette page n'est pas très appropriée, bien que consacrée. Nous allons voir que l'équation porte sur la variation de température dans le temps et l'espace et non sur une quantité un peu floue nommée "chaleur". Il est bien plus précis de parler d'équation de diffusion thermique, terme qui est d'ailleurs utilisé maintenant dans les classes prépa. De même, on ne parle plus de chaleur, mais d'énergie thermique.

L'équation de diffusion thermique est historiquement liée à Joseph Fourier. Ce dernier naquit en 1768 à Auxerre, où il étudie dans une école militaire (lui aussi !). Fin 1794, il est élève à Normale Sup (Ulm, il n'y avait que celle-ci) et en 1795-1796, il enseigne la physique à Normale Sup et à l'X !! Il débute ses travaux sur la chaleur en 1802 lorsqu'il est nommé préfet de l'Isère par Napoléon : douce vie de fonctionnaire. Il publie sa théorie sur la chaleur (et l'analyse de Fourier) en 1822. Il est élu à l'Académie française en 1826 et décéde à Paris en 1830.

Un peu de physique

Energie interne et transfert thermique d'énergie

Sur la cheminée à coté de mon bureau est posé un cube d'allume-feu, vous savez cet objet dont on se sert pour allumer le feu dans la cheminée ou le barbecue. Cet objet possède une énergie potentielle non nulle dans le référentiel de mon bureau puisqu'il est posé sur la cheminée, en hauteur. Son énergie cinétique Ec est nulle dans ce référentiel. Son énergie mécanique est donc égale à son énergie potentielle Ep. Mais est-ce son énergie totale ?

Si vous allumez ce cube, sa combustion se déclenchera quasi immédiatement et vous pourrez noter une brusque élévation de la température ambiante à sa proximité. Il y a donc dégagement d'énergie thermique. La combustion du cube est une réaction chimique d'oxydation des molécules qui composent le cube. Le cube stockait donc une énergie interne qui est "libérée" au cours de la combustion. Cette énergie interne que l'on note classiquement U est donc une des composantes de l'énergie totale du cube : son énergie potentielle Ep, son énergie cinétique Ec et son énergie interne U.

Cette énergie interne d'un corps peut prendre plusieurs formes. Elle est due par exemple à sa température, si le corps n'est pas au zéro absolu. Comme vous le savez, la température mesure l'agitation des molécules ou des atomes au niveau microscopique. Mais il peut aussi s'agir de l'énergie chimique susceptible d'être libérée lors d'une réaction chimique (notre cas) ou bien d'une énergie nucléaire pour les corps radioactifs.

Le premier principe de la thermodynamique nous apprend qu'une variation de cette énergie interne ne peut se traduire que par un travail W de forces externes ou/et par une libération, un transfert thermique d'énergie, Q, que l'on appelait autrefois la chaleur. Cela s'écrit \( \delta U = W + Q\). Dans l'expression "transfert thermique d'énergie", l'adjectif thermique s'applique bien au transfert et non à l'énergie. En effet, il n'existe pas de forme particulière d'énergie "thermique", comme il existe une énergie cinétique ou potentielle. Il n'existe qu'une énergie interne U, dont les variations peuvent se traduire sous forme de travail ou de transfert thermique.

D'autre part, nous savons que les lois de conservation nous imposent une logique très simple. La variation d'énergie temporelle est égale à : (ce qui rentre - ce qui sort à travers la surface d'interface) + la génération interne volumique (par réaction chimique ou nucléaire par exemple). Dans notre cas, ce dernier terme sera nul. Voici une ligne de conduite simple !

Dans cette page, nous allons nous intéresser au cas où la variation d'énergie interne de notre système (objet quelconque) est entièrement due à un transfert thermique d'énergie, sans travail, alors que notre système ne sera pas isolé puisqu'en contact avec un thermostat, c'est à dire un dispositif qui lui impose une température.

La conduction thermique

Le cours de physique de TS nous apprend qu'il existe trois modes de transfert thermique :

Nous avons tous l'expérience de la conduction thermique qui nous occupera ici : un objet chaud plongé dans l'air ambiant de notre cuisine se refroidit ! Et a priori, nous n'avons jamais vu une casserole froide se rechauffer spontanément en puisant de la chaleur dans l'air ambiant de la cuisine.

Nous pouvons donc émettre l'hypothèse qu'il existe un flux d'énergie qui circule des régions de température chaude vers les régions de température froide, et qui tend à uniformiser la température entre ces régions.

Au niveau microscopique, ce transfert d'énergie d'une région chaude vers une région froide se produit différemment selon le milieu.

La conduction thermique dans les gaz et les liquides

Dans les gaz et les liquides la mobilité des molécules est grande, plus grande dans les gaz que dans les liquides. Ces molécules possèdent une énergie cinétique et s'entrechoquent. Le transfert thermique d'énergie s'effectue dans ces milieux par le transfert d'énergie cinétique qui se produit au moment du choc entre une molécule rapide, plus "chaude" et une molécule plus lente, plus "froide". J'ai mis des guillemets parce que les termes "chaud" et "froid" sont impropres pour une molécule, ce sont des termes qui s'appliquent statistiquement à une population de molécules.

La conduction thermique dans les solides

Dans les solides, les molécules ou les atomes sont figés dans un réseau maillé qui empêchent les grands déplacements. Le seul mouvement possible est un mouvement de vibration autour de leur position d'équilibre dans le réseau. Le transfert thermique d'énergie se traduit par une augmentation plus ou moins grande de l'amplitude de ces vibrations. Lorsque l'énergie apportée est suffisante pour que l'amplitude de la vibration dépasse une certaine valeur dépendante du réseau, alors les molécules se libèrent du réseau et le solide fond et s'évapore...

Le cas des métaux est un peu particulier : eux possèdent des électrons de conduction, qui circulent librement sur le réseau du solide. Ces électrons se comportent comme les molécules d'un gaz (on parle de gaz d'électrons) et dans ce cas, à l'augmentation de l'amplitude des vibrations sur le réseau s'ajoute le transfert d'énergie cinétique des électrons rapides, "chauds" vers les électrons lents, "froids".

La loi de Fourier

Notion de champ de température

Nous avons déjà rencontré la notion de champ en électromagnétisme, notion qui est maintenant au programme de 1ere S. Pour rappel, un champ associe à chaque point désigné par ses coordonnées (x,y,z) de l'espace eucliden de la physique classique une grandeur scalaire ou vectorielle. En l'occurence, le champ de température associe à chaque point de l'espace une grandeur scalaire nommée "température", qui mesure l'agitation thermique en ce point.

Rappelons aussi que la température se mesure en Kelvin (K).

Ce champ peut être localement uniforme, par exemple dans une région où toutes les températures sont identiques. Mais en général, il est variable, présentant des régions plus chaudes et d'autres régions plus froides. Entre une région chaude et une région plus froide, il existe donc une différence de température, plus ou moins intense et brutale. Le gradient de température mesure cette différence. Dans une région où la température est constante, le gradient de température est nul.

La variation de température dans un champ de température peut dépendre à la fois de l'espace mais aussi du temps, ce qui complique un peu les choses... Nous verrons un peu plus bas comment traiter ces variations multiples.

Notion de flux d'énergie

La notion de flux d'énergie est très courante en physique. Nous l'avons déjà rencontré sous le nom d'intensité du courant électrique. Dans ce cas, il s'agit du flux de charges électriques qui traversent une surface unitaire donnée par unité de temps. Nous l'avons noté \( I_q = \frac{dQ}{dt} \).

Par analogie, on définit le flux thermique, plus exactement le flux d'énergie interne par \( I_u = \frac{dU}{dt} \). C'est la quantité d'énergie interne U qui traverse, sans travail, une surface unitaire par unité de temps.

Pour l'exprimer plus précisément, on introduit un nouvel objet par anologie avec le courant électrique. Il s'agit du vecteur de courant volumique d'énergie interne sans travail noté classiquement \( \vec{J_u}\) et souvent nommé improprement "vecteur courant thermique", ce qui nous permet d'écrire l'expression du flux de ce vecteur à travers une surface dS orientée vers l'extérieur par le vecteur normal \( \vec{n} \), ce qui nous donne \( I_u = \int_S \vec{J_u} . \vec{n} dS \). Le signe de \( I_u \) dépend du sens du flux à travers la surface. Il est négatif pour le flux entrant et positif pour le flux sortant (pensez au signe du produit scalaire sous l'intégrale...).

Pour information, le courant volumique d'énergie interne sans travail a pour dimension W.m-2, comme le vecteur de Poynting...

L'expression de la loi de Fourier

Après ces préambules utiles, venons en à la loi de Fourier proprement dite. Pour la petite histoire, il découvrit cette loi au tout début du XIX siècle (1805 ou 1807 selon les sources) et la publia en 1822. Selon d'autres sources, ce serait Jean-Baptiste Biot qui l'aurait exhibé mathématiquement en 1804. Mais ce n'est pas très important.

Ce qu'il est important de se souvenir, c'est qu'il s'agit d'une loi tirée de l'expérience et qu'elle possède des limites d'application à ne pas dépasser. En particulier, il ne faut pas que les variations de température, le gradient, soit trop fortes ou trop faibles. On trouve en physique d'autres lois analogues, dont la plus connue est sans doute la loi d'Ohm, qui possède exactement la même forme.

Pour faire simple, ici et dans la suite, on va se placer dans le cadre d'un problème unidimensionnel, c'est à dire que le transfert thermique d'énergie se fait sur une dimension Ox.

La loi de Fourier relie, après constat expérimental, le vecteur de courant volumique d'énergie interne sans travail \( \vec{J_u}\) avec le gradient de température \( \vec{grad} T \). La relation est linéaire et s'écrit \( \vec{J_u} = - \kappa \vec{grad} T \). La linéarité de l'équation n'est due qu'aux approximations que nous avons fixé à son domaine de validité, comme pour la loi d'Ohm.

Le paramètre \( \kappa \) est appelé conductivité thermique, toujours positif, de dimension W.m-1.K-1. Notez la présence du signe -, qui résulte du second principe de la thermodynamique : le flux d'énergie va des régions aux températures les plus hautes vers les régions aux températures les plus basses.

Dans notre hypothèse d'un problème unidimensionnel, en développant le gradient, on obtient l'équation \( \vec{J_u} = - \kappa \frac{\partial T}{\partial x}\vec{u} \), soit en projetant sur Ox \( J_u(x,t) = - \kappa \frac{\partial T(x,t)}{\partial x} \). Nous obtenons une équation à deux variables, la position x et le temps t, avec une dérivée partielle.

Pour le lecteur qui ne serait pas familier avec les équations aux dérivées partielles, je l'invite à aller les découvrire dans cette page.

L'équation de diffusion thermique

Comment obtenir cette équation

Nous allons l'établir en utilisant la loi de Fourier décrite ci-dessus et le principe de conservation d'énergie. Nous resterons dans le cadre du problème unidimensionnel, sachant que l'extension en dimension 2 ou 3 n'est pas très compliquée, mais trop lourde pour notre étude.

Considérons un élément de volume dV orienté selon l'axe Ox de propagation du flux d'énergie, limité par deux surfaces dS, l'une entrante et l'autre sortante, et d'épaisseur dx. Appelons \( \vec{J_uE} \) le vecteur de courant d'énergie volumique entrant et \( \vec{J_uS} \) le vecteur de courant d'énergie volumique sortant. Supposons que ces deux vecteurs soient normaux aux surfaces entrantes et sortantes.

Cet élément de volume dV est immobile et son énergie potentielle n'est pas modifiée par hypothèse. Selon le premier principe de la thermodynamique, la variation d'énergie interne dU n'est donc attribuable qu'à la variation dQ, le transfert thermique d'énergie.

Calculons la variation d'énergie interne dans ce volume dV, de masse volumique \( \rho \), en utilisant la définition de \( \vec{J_u} \) donnée plus haut dans le cas unidimensionnel. Nous obtenons après simplification, en égalant dU et dQ :
\( \frac{\partial(\rho u)}{\partial t}dx = -J_{u,x}(x + dx,t) + J_{u,x}(x,t) \)
En remarquant que \( J_{u,x}(x + dx,t) - J_{u,x}(x,t)= \frac{\partial J_{u,x}}{\partial x}dx \), on obtient finalement :
\( \frac{\partial(\rho u)}{\partial t} = - \frac{\partial J_{u,x}}{\partial x} \).

Dans cette dernière équation, remplaçons dans le terme de droite \( J_{u,x} \) par sa définition donnée par la loi de Fourier, on obtient : \( \frac{\partial(\rho u)}{\partial t} = - \frac{\partial}{\partial x}(- \kappa \frac{\partial T}{\partial x}) \) ou en condensant l'écriture :
\( \frac{\partial(\rho u)}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2} \)

Vous aurez noté ici, si vous êtes attentifs, que j'ai considéré que \( \kappa\) était constant puisque je l'ai sorti de la dérivée sans autre forme de procés ! C'est un peu osé, et vrai seulement si le milieu est isotrope (le matériaux est homogène) et si l'on ne chauffe pas trop fort ou trop vite, parce que sinon, il devient dépendant de la température.

Reste maintenant à traiter le terme de gauche. Pour ce faire, je vais faire appel à l'expression de la capacité calorique qui relie les variations de l'énergie avec les variations de température. On peut donc écrire que, si \( \rho u \) désigne l'énergie interne volumique :
\( \frac{\partial(\rho u)}{\partial t} = \rho c_v \frac{\partial T}{\partial t} \)
avec \( c_v\) la capacité thermique massique à volume constant.

En reportant cette expression dans le terme de gauche de notre équation, j'obtiens :
\( \rho c_v \frac{\partial T}{\partial t} = \kappa \frac{\partial^2 T}{\partial x^2} \)
soit en regroupant les termes constants à droite de l'équation :
\( \frac{\partial T}{\partial t} = \frac{\kappa}{\rho c_v} \frac{\partial^2 T}{\partial x^2} \)

Pour simplifier l'écriture, je vais appeler D le rapport \( \frac{\kappa}{\rho c_v} \). Ce paramètre est la diffusivité thermique du matériau constituant notre élément de volume. La dimension de D est m2.s-1. Finalement, j'obtiens l'équation :
\( \frac{\partial T}{\partial t} = D \frac{\partial^2 T}{\partial x^2} \)
qui constitue l'équation de diffusion thermique !

Notions de conditions initiales et de conditions aux limites

Vous vous souvenez sans doute que lorsque nous cherchons une solution particulère à une EDO, nous fixons une ou plusieurs conditions initiales qui fixent la valeur de la ou des constantes de la solution générale de l'EDO. Il y a autant de constantes et donc de conditions initiales que l'ordre de l'EDO. Par exemple, l'équation de Newton, qui est une EDO d'ordre 2 requiert 2 conditions initiales : la position initiale et la vitesse initiale, initiale voulant dire dans ce cas à t=0.

Nous retrouvons le même principe lorsqu'il s'agit d'équation aux dérivées partielles (EDP), auquel il faut ajouter la nécessité de définir les conditions aux limites du domaine de définition. Dans le cas de l'équation de diffusion thermique, la condition initiale est la température du corps en début d'expérience et les conditions aux limites sont les températures qui règnent aux extrémités de la barre, voire à sa périphérie si la barre n'est pas isolée sur sa longueur. Etablir les conditions initiales et les conditions aux limites d'un problème modélisé par une EDP est délicat, car cela a des conséquences mathématiques et physiques non négligeables. D'ailleurs, les mathématiciens ont classé les problèmes de physique en plusieurs classes caractérisées par ces conditions (voir la page au sujet de cette classification).

Le problème que nous allons aborder est un problème de Dirichlet (mathématicien allemand 1805-1859) : la température initiale de la barre est donnée et la température de ses deux extrémités est identique et constante pendant toute l'expérience. Le problème de Dirichlet consiste à trouver une fonction u(x,t) sur un domaine spatial [0,L] pour t > 0 telle que:
\( u_t(x,t) = Du_{xx}(x,t)\) pour \(0 \le x \le L \) et \(t > 0\)
\( u(x,0) = u_0(x) \) pour \(0 \le x \le L \), qui est la condition initiale
\( u(0,t) = T_0 \) et \(u(L,t) = T_L \) pour \(t > 0\), les conditions aux limites

Une remarque : j'ai utilisé la notation anglo-saxonne pour les dérivées partielles : \( u_t(x,t) \) signifie \( \frac{\partial u(x,t)}{\partial t}\) et \( u_{xx}(x,t) \) signifie \( \frac{\partial^2 u(x,t)}{\partial x^2}\). C'est plus rapide !

Quelques remarques sur l'équation de diffusion thermique

Invariance par rapport au temps

Il est facile de constater que notre équation n'est invariante par renversement du temps, car sa dérivée par rapport au temps est du premier ordre. Cette remarque est fondamentale car elle prouve que le phénomène de diffusion thermique est un phénomène irréversible par essence.

Sa linéarité en T

Notre équation est linéaire en T et donc on peut appliquer le théorème de superposition aux solutions. Mais attention dans ce cas aux considérations sur les conditions aux limites, qui doivent prendre en compte la superposition.

Unicité de la solution

Cette équation possède un ensemble de solutions analytiques bien connu. En effet, les solutions dépendent des conditions aux limites que l'on impose au modèle.

Résolution numérique de l'équation de diffusion thermique

Le modèle physique

Nous allons considérer une barre métallique homogène de faible diamètre, de telle sorte que nous puissions négliger ses dimensions spatiales autre que sa longeur. Autrement dit, je m'arrange pour avoir un modèle approximativement 1D.

Je vais considérer que cette barre est isolée sur toute sa longeur et que seules ses deux faces sont en contact avec l'extérieur. Cela pour pouvoir considérer que les flux d'énergie thermiques ne concernent que les deux parois extrémales de la barre. Bien sur, il n'existe aucune production d'énergie interne, par réaction chimique ou nucléaire. Notre dispositif ressemble à cela :

Schéma barre

Je chauffe de manière uniforme cette barre dans un four à une température de 500 K pour fixer les idées. Puis je la sors du four et je la laisse refroidir à l'air ambiant, disons à 293 K toujours pour fixer les idées.

Comme vous l'avez sans doute deviné, je pose les conditions initiales et les conditions aux limites de mon problème ! Voyons maintenant comment le traiter.

Adimensionnement de l'équation

Nous avons déjà abordé le sujet : en simulation numérique, il est de rigueur d'adimensionner les équations pour faciliter leur analyse. Cela revient à supprimer les paramètres physiques par un changement de variables astucieux.

Dans le cas de l'équation de la chaleur, je vais procéder à deux changements de variables. Je vais poser x = X/L, avec X la variable d'abscisse et L la longuer de la barre. J'obtiens donc une variable en abscisse qui n'a plus de dimension. Et comme je vais fixer arbitrairement L = 1, j'aurais donc la même valeur pour x et X.

C'est un peu plus compliqué pour le temps, parce que je veux virer la constante D qui me gêne, et dont la dimension est [L][L]/[T]. Je vais donc définir une variable \( \tau \) telle que \( t = \frac{t'}{\tau}\). Il faut donc que \( \tau \) est la dimension d'un temps. Nous sommes dans le domaine de la diffusion et il existe une constante bien connue qui est le temps caractéristique de diffusion : ça tombe ! Une très rapide analyse dimensionnelle me conduit à la valeur \( \tau = \frac{L^2}{D}\).

Me voilà donc avec une équation de diffusion thermique adimensionnée, égale à \( \frac{\partial T}{\partial t} = \frac{\partial^2 T}{\partial x^2} \), équation que je vais maintenant discrétiser.

Le choix des schémas numériques

Un schéma numérique est un algorihme basé sur des considérations et des résultats d'analyse numérique, qui nous permettra de chercher les solutions numériques à une équation différentielle ordinaire ou aux dérivées partielles.

Dans la suite, nous utiliserons les schémas numériques basés sur la méthode des différences finies. Il s'agit d'une méthode très usitée en physique numérique pour construire les solutions aux EDP. Ces méthodes sont décrites dans cette page, que je vous invite à consulter.

La discrétisation de l'expérience

L'évolution de la température de la barre évolue dans l'espace et dans le temps de façon continue. Notre objectif est de calculer la température en chaque point de la barre à un instant donné de l'expérience. Quand j'écris "en chaque point de la barre", c'est bien entendu une vue très théorique. Je ne peux pas faire le calcul en tous points de la barre: je suis obligé d'en choisir quelques uns. Je transforme l'espace continu de la barre en un ensemble représentatif de points de celui-ci. Je "discrètise" l'espace. Il en est de même pour le temps: je ne peux pas faire le calcul pour tous les instants ! Aussi, je vais également "discrétiser" le temps. En bref donc, je découpe l'espace, en l'occurence la longueur de la barre en un certain nombre d'intervalles de largeur dx et je découpe le temps en intervalles de largeur dt.

Pour fixer les idées, j'imagine une grille dont chaque noeud représente un élement de l'espace-temps de mon expérience, la voici :

Principe Schéma FTCS

La ligne du bas représente l'état de la barre à l'instant t0 de l'expérience. En chaque noeud de cette ligne la température de la barre est définie. Comme cette grille ressemble furieusement à une matrice, je vais employer la notation matricielle pour désigner la température à ce noeud T(x,t), qui se lit "température au noeud d'abscisse x et d'ordonnée t". On appelle cette grille un maillage.

A l'instant t=0, la température de la barre est égale à la température initiale, que je note TBarreInit. Partout, sauf aux extrémités, parce que nous avons convenu que la barre étant à l'air ambiant et que le contact thermique était parfait. Donc, pendant toute la durée de l'expérience, les deux extrémités de la barre sont à une température constante, que je note TExt. Sur la grille cela est schématisé par les noeuds noirs sur ses bords gauche et droit.

Le but de la manoeuvre est de calculer la température T pour tous les noeuds T(x,t) de la grille, pour x compris dans l'intervalle [1,L-1] et pour chaque instant t > 0. Pour ce faire, il existe plusieurs méthodes. Nous allons en voir quelques unes en commençant pour le schéma FTCS.

Utiliser le schéma FTCS

le schéma FTCS

Le schéma FTCS (Forward Time Centred Space) est décrit dans la page sur la méthode des différences finies, mais je vais en rappeler le principe ici. C'est un schéma explicite d'ordre 1 en temps et d'ordre 2 en espace.

Il s'agit de discrétiser les deux dérivées partielles en utilisant un développement de Taylor. Le membre de gauche peut donc s'écrire en approximant :
\( \frac{\partial T(x,t)}{\partial t} \approx \frac{T(x_x, t_t+dt) - T(x_x,t_t)}{dt} = \frac{T(x,t+1) - T(x,t)}{dt}\)

En suivant la même méthode, je discrétise le membre de droite, ce qui me donne :
\( \frac{\partial^2 T(x,t)}{\partial x^2} \approx \frac{T(x+1,t) + T(x-1,t) - 2T(x,t)}{dx^2}\)

D'où mon EDP discrétisée :
\( \frac{T(x,t+1) - T(x,t)}{dt} = \frac{T(x+1,t) + T(x-1,t) - 2T(x,t)}{dx^2} \) ou encore :
\( T(x,t+1) = T(x,t) + \frac{dt}{dx^2}(T(x+1,t) + T(x-1,t) - 2*T(x,t)) \)

Si nous revenons à notre grille, vous comprendrez vite le nom du schéma : je calcule la température en un noeud et à un instant donné en utilisant la température à l'instant précédent pour les noeuds situés en dessous et de part et d'autre du noeud courant.

Ce schéma est très simple, mais présente un inconvénient majeur: il est instable ! Pour certaines valeurs de dx et dt, il explose littéralement en donnant des valeurs astronomiques. Il reste stable sous la condition que \( dt \lt (dx*dx)/(2.0)\), condition qu'il faudra tester dans notre script. Pour information, cette condition s'appelle la condition de Courant-Friedrich-Levy, ou CFL.

L'initialisation des variables

Le script ECPython1.py commence classiquement par l'import des librairies nécessaires au calcul et au tracé des résultats. Rien de spécial à signaler. Puis vient l'initialisation des constantes et des variables. Le critère de stabilité est testé. S'il n'est pas respecté, le script s'arrête.

Nx est le nombre de noeuds spatiaux et Nt le nombre de noeuds temporels. Les différents buffer de travail et la grille de calcul sont définis par les instructions:

X = linspace(0.0,1,Nx)

duree = linspace(0.0, tfin,Nt)

T = zeros((Nt,Nx))

Les conditions initiales sont définies par:

for i in range(1,Nx-1):

    T[0][i] = TBarreInit

Les conditions aux limites par :

for t in range(0,Nt):

    T[t][0] = T[t][-1] = TExt

L'implémentation du schéma FTCS

Elle est très simple : il suffit de coder l'équation ! Voilà ce que cela donne :

c = dt/(dx*dx)

TT = zeros(Nx)

for t in range(0,Nt-1): # boucle de discrétisation temporelle

    TT[0] = TT[-1] = TExt

    for x in range(1,Nx-1): # boucle de discrétisation spatiale

        TT[x] = T[t][x] + c*(T[t][x-1] - 2*T[t][x] + T[t][x+1])

    T[t+1] = TT.copy() # recopie de T[t] dans T[t+1]

Le reste des instructions concerne l'affichage de la surface obtenue avec son codage en fausses couleurs, rien que de très classique et qui est largement réutilisé dans tous les codes de calcul.

Voilà le résultat obtenu avec le script ECPython1.py qui implémente le schéma FTCS :

Schéma FTCS

Vous observez le refroidissement de la barre : la barre reste chaude en son milieu et la température chute en se rapprochant des bords. Il y a diffusion thermique par les extrémités.

Utiliser la décomposition en EDO

Le principe

Cette technique est basée sur le même schéma numérique et utilise une remarque intéressante. En considérant chaque ligne de notre grille de calcul, qui correspont à un instant donné, nous pouvons écrire :
\( \frac{dT}{dt} = \frac{T(x+1,t) + T(x-1,t) - 2T(x,t)}{dx^2}\)
Oui, vous avez bien lu, j'utilise ici une dérivée ordinaire, car dans le membre de droite, à un instant donné, je peux considérer que t est une constante.

Cela revient donc à résoudre deux EDO, l'une en x et l'autre en t, dans un système de deux EDO. Et ça, nous savons faire facilement avec odeint(), une fonction bien connue. Je l'ai implémenté dans le script ECPython2.py

L'implémentation

Je ne reviens pas sur les initialisations de variables, buffers et constantes, qui sont identiques à ECPython1.py.

Les deux points saillants de ce script sont dans la description du système différentiel, qui calcule l'évolution de la température sur le segment à t=constante :

def Xdiscret(T,t):

   # initialisation

   dTdt = zeros(X.shape)

   # calcul sur le segment

   dTdt[1:-1] = (T[:-2] - 2*T[1:-1] + T[2:]) / dx**2

   return dTdt

et dans son utilisation , que vous connaissez maintenant par coeur :

Tevol = odeint(Xdiscret, init,t)

La solution se trouve dans le buffer Tevol, que j'affiche avec les mêmes instructions que précédement.

Voilà le résultat obtenu avec le script ECPython2.py qui implémente la décomposition en deux EDO de première ordre :

odeint

Vous noterez l'économie de moyens et le résultat, qui rend cette méthode très intéressante pour le calcul d'un schéma FTCS. On peut aussi l'utiliser pour d'autres schémas numériques.

Utiliser le schéma de Crank-Nicolson

Le schéma de Crank-Nicolson

Je vous renvois à cette page pour la description du schéma de Crank-Nicolson (CN). C'est un schéma basé sur la transformation matricielle du schéma numérique. Ce dernier point présente des avantages de calcul importants lorsqu'on emploie des langages spécialisés dans le calcul matriciel : MatLab, SciLab, ou Python qui possède une excellente libraire matricielle. Le schéma CN est un schéma d'ordre 2 en temps et en espace.

L'implémentation du schéma CN

Elle est un peu plus compliquée que les deux codes précédents. L'initialisation des données est identique à l'exception des conditions aux limites, qui prennent une forme matricielle :

T[:,0] = TExt*ones(Nt+1,float) # bord gauche

T[:,Nx] = TExt*ones(Nt+1, float) # bord droit

Attention aussi à la valeur de c, le coefficient de discrétisation, qui vaut ici :

c = dt/(2.0*(dx*dx))

Pour construire les deux matrices diagonales, j'ai écris une petite fonction qui fait ça, sans appel aux fonctions de scipy :

def MatConst(n,c1,c2):

    M = zeros((n+1,n+1),float)

    M[0,0] = M[n,n] = 1

    for i in range(1,n):

    M[i,i]= c2

    M[i,i-1] = M[i,i+1]= -c1

    return M

La construire les deux matrices tridiagonales caractéristiques du schéma NC se réduit donc à :

A = MatConst(Nx,c,1 + 2*c)

B = MatConst(Nx,-c, 1 - 2*c)

Et enfin, le codage du schéma CN en lui-même, qui est des plus simples grâce aux fonctions matricielles de Python :

for t in range(0,Nt):

    T[t+1,:] = solve(A,dot(B,T[t,:]))

La fonction solve() qui permet de résoudre l'équation matricielle A*X = B et la fonction dot() qui calcule le produit matriciel de deux tableaux.

Les intructions d'affichage restent identiques.

Voilà le résultat obtenu avec le script ECPython3.py qui implémente le schéma de Crank-Nicolson :

Crank-Nicolson

Comme vous le remarquez, nous obtenons à peu près les mêmes résultats qu'avec le schéma FTCS. Cependant, je vous conseille vivement d'utiliser le schéma CN pour vos calculs sur les EDP, car sa stabilité est remarquable ! En effet, j'ai utilisé dans cet exemple un dt = 10^-4, au lieu de 10^-5, pour un dx de 10^-2 dans le schéma FTCS, et mon schéma n'a pas explosé !

Les scripts Python

Les scripts Python étudiés dans cette page sont disponibles dans le package EDTPython.zip :

Pour conclure

A travers l'exemple simple de la résolution numérique de l'équation de diffusion thermique, nous avons abordé les EDP, chapitre important de la physique numérique. On en trouve partout ! Aussi, j'aborderai dans les pages à venir d'autres EDP célébres : Poisson, Laplace, etc...

A l'occasion, vous pouvez creuser tous les aspects mathématiques des schémas numériques, en particulier leur stabilité et leurs conditions aux limites... Amusez-vous bien !

Contenu et design par Dominique Lefebvre - www.tangenteX.com février 2016 -- Pour me joindre par mail ou PhysiqueX ou

Licence Creative Commons

Cette œuvre est mise à disposition selon les termes de la Licence Creative Commons Attribution - Pas d’Utilisation Commerciale - Pas de Modification 3.0 France.