Équ. aux dérivées partielles: schéma explicite centré

Transport de la chaleur

La chaleur se propage dans un matériau au repos par diffusion. Le flux de chaleur est proportionel au gradient de température. Si le milieu est homogène, la distribution de température obéit à l'équation

∂u/∂t - k ∂²u/∂x² = 0

La même équation peut décrire d'autres systèmes diffusifs, comme par exemple la dilution de colorant dans de l'eau. La constante k s'appelle la diffusivité thermique du matériau dans notre cas, et constante de diffusion en général.

Application 1

Considérons la cuisson d'un panino qui est initialement à température ambiante (disons 20°C) et qui est placé entre deux fers chauffants à 220°C. Les fers du gauffrier se maitiennent à cette température et imposent donc la température sur les deux faces du panino avec lesquels ils sont en contact. En supposant le panino nettement plus plat que large ou long, en combien de temps l'intérieur dépassera-t-il les 80°C qu'il faut par exemple pour faire fondre le fromage ? On supposera que le panino fait 1cm d'épaisseur, et possède une diffusivité thermique de 0.0004 cm²/s. On négligera les enthalpies de cuisson et d'évaporation de l'eau contenue dans le panino.

Application 2

Considérons un anneau (un tore) de longueur 1 qui à été coupé en 2 parties identiques, selon un plan de coupe qui comprend l'axe de symétrie de l'anneau. Les 2 parties sont chauffées à deux températures différentes (disons 1 et 2), puis recollées ensemble à l'instant t=0. La thermodynamique nous enseigne que l'ensemble va évoluer vers un état où la température sera uniforme. Calculez la distribution de température le long de l'anneau pendant tout le transitoire, numériquement et analytiquement (adimensionnez l'équation tel que k=1).

L'étude pas à pas

Simulation numérique

Pour simuler le système on discrétise l'espace et le temps: on ne décrit la température qu'en des points j*Δx, j=0,1,2,... distants de Δx, à des instants t=n*Δt, n=0,1,2,... sépares de Δt. Si L est l'épaisseur de notre panino, la distribution de température à un instant donné est décrite par un ensemble de L/Δx + 1 points u_jⁿ, j=0,1,2,...,L/Δx.

Écrivez un programme qui calcule la diffusion de la chaleur en utilisant le schéma centré explicite suivant:

(u_jⁿ⁺¹ - u_jⁿ) / Δt - k (u_j+1ⁿ - 2*u_jⁿ + u_j-1ⁿ) / Δx² = 0

La condition de stabilité de ce schéma est

k Δt / Δx² <= 0.5

Vous pouvez vous servir du squelette de programme suivant comme point de départ. Affichez la distribution de température à intervalles réguliers.

#include <stdio.h>
#include <stdlib.h>

#define LONGUEUR 1.0
#define DUREE 0.1
#define NX 100
#define NT 10000
#define K 1.0

/* xmalloc: remplacement de la fonction malloc */
/* l'appel à malloc est encapsulé dans divers vérifications; xmalloc
   préfère mourir (et entraîner le programme dans la mort :-P) plutôt
   qu'avouer son échec */
void * xmalloc(size_t len)
{
  void *buf;
  if (len == 0) return NULL;
  buf = malloc(len);
  if (buf == NULL) {
    fprintf(stderr, "Erreur d'allocation de mémoire, fatale.\n");
    exit(EXIT_FAILURE);
  }
  return buf;
}

/* xfree: pendant de xmalloc pour libérer la mémoire */
void xfree(void *buf)
{
  if (buf != NULL) free(buf);
}

int main (void) {
  double *u;
  u = xmalloc(NX*sizeof(double));

  /* ... à compléter ... */

  xfree(u);
  return EXIT_SUCCESS;
}

Simulez la cuisson du panino avec les paramètres donnés en introduction, en faisant attention de satisfaire la condition de stabilité, et répondez à la question posée.
Changez le pas de temps pour avoir k Δt/Δx² = 0.5 et k Δt/Δx² = 0.55 et relancez la simulation. Que ce passe-t-il ?

Deuxième cas: le tore de Fourier

Modifiez le programme pour tenir compte de la géométrie périodique de l'anneau, l'absence de sources et les bonnes conditions initiales. Simulez dans des conditions stables du schéma utilisé.

Solution exacte et comparaison

Pour comparer, nous allons calculer une solution analytique approchée. On peut tenter d'écrire la fonction u(x,t) sous forme d'une série de Fourier:

u(x,t) = a₀/2 + ∑_p=1,2,...[ ( a_p cos(2π p x) + b_p sin(2π p x) ) u_p(t) ]

Montrez que les fonctions u_p(t) satisfont l'EDO
```
du_p/dt + (p π/L)² u_p = 0
```
qui a la solution
```
u_p(t) = c_p exp[ - (p π/L)² t ]
```
(on peut absorber la valeur de c_p dans les coefficients a_p et b_p, c'est-à-dire prendre c_p=1)
Determinez les coeffecients a_p et b_p à partir de la condition initiale pour la température du tore:
```
a_p = ∫₀^L u(x,0) cos(2π p x) dx
b_p = ∫₀^L u(x,0) sin(2π p x) dx
```
Utilisez les dix premiers termes (ou plus, si vous avez la solution générale) pour calculer la solution "exacte" dans votre simulation en plus de la résolution numérique. Comparez les deux.

Auteur(s) : A. Daerr. Dernière modification : Sat Oct 25 19:59:41 2008. [validate XHTML]