Equation de la chaleur en instationnaire

Le but de ce TP est de calculer l'évolution de la température dans une barre au cours du temps. L'équation non-stationnaire de la chaleur en 1D sera modélisée à l'aide d'une équation aux dérivées partielles avec les conditions aux limites suivantes:

où f désigne la source de chaleur à t=0 et est le coefficient de diffusion.

Question 1

Tout comme au TP précédent on construit un maillage du domaine Ω par N+1 segments de taille

Les N points intérieurs sont:

Calculer les points de discrétisation xi et stocker les dans un vecteur de taille N. Calculer également la condition initiale suivante :

Dans un module, écrire une routine WriteText qui prend en argument deux vecteurs x, y et un nom de fichier. La routine doit écrire dans le fichier les vecteurs x et y sur deux colonnes :

x_1 y_1
x_2 y_2
...
x_N y_N

Appeler ensuite cette sous-routine dans le programme principal :

call WriteText(x, F, "Un0000.dat")

Une fois le code exécuté, on visualisera la condition initiale sous ipython.

ipython --pylab
V = loadtxt("Un0000.dat")
plot(V[:, 0], V[:, 1])

Question 2

On choisira comme pas de temps :

où cfl < 1 est le coefficient de Courant, Friedrichs, Lewy. On utilise le schéma aux différences finies centré en espace et un schéma d'Euler explicite en temps :

Ce schéma reste vrai pour les deux extrêmités car du fait de la condition de Dirichlet, on note par convention :

Ecrire une sous-routine qui implémente le schéma d'Euler explicite. Tester ce schéma avec

Afficher le résultat final avec ipython pour différents temps finaux (variable Tmax), pour Tmax=0.0005, 0.001, 0.005, 0.01 par exemple.

Question 3

Dans le cas de la condition initiale gaussienne sur un intervalle infini, la solution analytique vaut :

Ecrire une sous-routine qui prend les valeurs de x et le temps t en paramètres et génère u. Calculer alors l'erreur au temps final Tmax = 0.001 entre la solution numérique et la solution analytique. Observer que cette erreur est divisée par 4 lorsqu'on multiplie par 2 le nombre d'intervalles (la cfl restant identique par ailleurs). Quel est l'ordre du schéma ?

Question 4

Afin de voir l'évolution de la solution, on pourra imprimer la solution tous les Ndisplay itérations en écrivant dans le programme principal

Ndisplay = 10
// boucle principale en temps
do nt=1, nb_iter

  // on avance le schema
  call AdvanceEuler(Un, Un_next, ...)
  Un = Un_next

  // on ecrit l'instantane
  if (mod(nt, Ndisplay) .eq. 0) then
     call WriteText(x, Un, "Un" // EntierToString(nt/Ndisplay) // ".dat")
  end if
end do

La fonction EntierToString convertit un entier en chaine de 4 caractères (en mettant des zéros au début) :

// convertit l'entier n en chaine de 4 caracteres
// par exemple EntierToString(23) renvoie '0023'
function EntierToString(n) result(chaine)
  
  implicit none
    
  integer :: n, n0
  character(len=4) :: chaine
    
  n0 = n 
    
  chaine(4:4) = achar(48+mod(n0,10))
  n0 = n0/10
    
  chaine(3:3) = achar(48+mod(n0,10))
  n0 = n0/10

  chaine(2:2) = achar(48+mod(n0,10))
  n0 = n0/10

  chaine(1:1) = achar(48+mod(n0,10))
    
end function EntierToString

Vous pouvez insérer cette fonction dans le module de votre choix. Télécharger ensuite le fichier film.py dans un terminal (en vous placant dans le répertoire de travail):

wget http://www.math.u-bordeaux.fr/~durufle/teaching/tp_chaleur_time/film.py

Dans ipython, vous taperez alors :

from film import *
film1D("Un", ".dat", 1, 41, -0.1, 1.1)

Vous devriez voir défiler les frames Un0001.dat jusqu'à Un0040.dat avec une échelle de [-0.1, 1.1] en y. La première ligne réalise un import du fichier film.py, et ne doit étre exécutée que si vous sortez de la console et que vous la relanciez.

Visualiser alors les résultats pour les différentes conditions initiales suivantes afin de retrouver les différents phénomènes vus en cours.

Afficher le résultat (l'animation) sous ipython pour Tmax=0.01 par exemple. Choisir ensuite cfl = 1.1 au lieu de 0.9 et N = 100. Visualiser le comportement de la solution numérique pour les différentes conditions initiales.

Question 5

On considère maintenant le schéma d'Euler implicite :

qui s'écrit sous la forme matricielle :

oû A est la matrice tridiagonale implémentée au TP précédent. En remarquant que pour μ > 0, est une matrice symétrique définie positive, vous pouvez réutiliser vos sous-routines du TP précédent pour résoudre le système matriciel. Implémenter le schéma d'Euler implicite dans une routine en utilisant ces sous-routines. On fera attention de ne factoriser la matrice qu'une seule fois, cette étape étant plus coûteuse que les descentes et remontées.

Reprenez les questions de la deuxième partie et comparez les comportements en termes d'erreur et de stabilité. Considérez plusieurs valeurs de cfl (1, 10, 50 et 100 par exemple), avec N=1000.

Question 6

Implémenter et comparer au schéma de Crank-Nicolson qui s'écrit matriciellement

Ce schéma est-il plus précis que le schéma d'Euler implicite ?