Bonjour à tous, j'ai besoin de votre aide. Je suis en stage dans le domaine de la mécanique des fluides (je ne suis pas informaticien de formation) et j'ai un problème sur mon code de résolution de l'équation de la chaleur à 2 dimensions avec terme advectif.
 
Le problème est le suisant:
 
J'ai un fluide (de l'eau) qui entre dans un canal à 20°C et qui est en contacte avec une paroi chauffante. Je résoud l'équation de la chaleur par différence finis en coordonnées cartésiennes. Mes conditions aux limites sont:
 
Une rampe de température de la paroi de 20°C à 100°C à en une seconde puis maintient à 100°C le reste du temps (bord gauche du maillage).
Une température d'entrée de 20°C
Une conditions de flux nul pour le bord droit du maillage.
En sortie je fais une extrapolation de la température gràce aux valeurs de températures précédentes
 
Voila mon code, je me suis inspiré d'un code que j'ai trouvé sur le net résolvant un problème de diffision de particule
 

 
Je n'arrive pas à visualiser correctement le champ de température. L 'écran est totalement noir. J'ai fait des prints pour regarder l'allure de mon tableau ui et les résultats sont physiques (le fluide se chauffe bien au contacte de la paroi) pourtant je n'arrive pas à les visualiser. Il faudrait que les zones chaudes se colorent en rouge au cours du calcul. Sauriez vous d'où vient le problème ?
 
Merci d'avance,
 
Fauster
 
edit du 25 Juillet : 8h13  prise en compte des remarques de Tentacle sur la stabilité et le schéma pour la dérivée d'ordre 1 (terme d'advection)
                                    Modification sur le calcul de u[i,j]. u[i,j] remplace u[1:-1, 1:-1] etc..
                                    Bonne prise en compte des conditions aux limites.
 
J'ai réussi à faire apparaître une image du champs de température (en bidouillant un peu)  mais je ne parvient toujours pas à itérer cette image au cours du calcul:
 

 
PS: Le résultat est bien physique: apparition d'une couche limite de température pendant le transitoire. (champs de température normalisé, le fluide rentre à gauche et va vers la droite)
 

Hum, j'ai peut être des problèmes de convergences...j'éclaircis cela

Si tu pouvais surtout encadrer ton code de la balise appropriée (balise "code" ), ce serait bien

Quoiqu'il en soit, lorsque le calcul converge, ça ne marche toujours pas. Je suis pas sur que ce soit ça la balise "code"

Bonjour,
 

Si j'ai bien compris ton problème, tu as effectivement des problèmes de convergence mais ça ne résoudra peut-être pas ton problème d'affichage (quoique, en cas d'instabilité, les valeurs explosent ce qui donne un affichage saturé) :
 

• dtm représente apparemment ton pas de temps maximal assurant la stabilité de ton schéma. Pour la diffusion, j'aurais plutôt mis min(dx2,dy2)*0.5/alphaf mais ton expression en est un minorant donc ok ; Le problème est que tu ne prends pas en compte la condition de stabilité de ton schéma d'advection qui devrait être (après correction du schéma) dt <= dy/abs(vo). Vu les paramètres de ton problème, tu devrais avoir dt <= 0.001 ce qui n'est pas le cas. Tu peux faire un truc du genre :

le coefficient 0.9 étant là pour contrôler le pas de temps par rapport à sa limite maximale (sachant que plus il est bas, plus ton schéma sera artificiellement diffusif).
 

• concernant ton schéma d'advection, il correspond à un schéma d'Euler explicite décentré mais il est décentré du mauvais côté : pour une vitesse positive, il faut décentrer à gauche (en amont par rapport au vecteur vitesse) et ici tu décentres à droite (le schéma est alors inconditionnellement instable). Il faudrait plutôt mettre :

à la place du dernier terme de ton schéma.
 
En espérant que ça fera avancer le smilblick
 
PS : tu affiches l'ensemble du vecteur ui or sur les bords du domaine, il vaudra toujours 0 au moment de l'affichage car il est copié à partir de u qui lui-même n'est calculé que sur les indices [1:-1, 1:-1]. Soit tu n'affiches que cette sous-matrice, soit il faut compléter ces valeurs avec les contraintes de ton problème.
 
PS2 : Je ne m'y connais pas beaucoup en python mais si j'ai bien compris, la méthode idle_add de gobject appelle updatefig tant que celle-ci renvoie True. Tu devrais peut-être faire en sorte de pouvoir contrôler ton résultat numérique toutes les n itérations ?

Merci à toi Tentacle. J'ai également lu un document qui expliquait pourquoi le terme advectif engendrait une instabilité vu son évaluation.  
 
J'ai un peu modifié le calcul de u|i,j] pour que cela soit plus clair, cela reste équivalent à l'ancienne version, bien que le calcul soit moins rapide.
 
Je cherche toujours à résoudre le problème d'affichage.  
 
J'ai tenté de normaliser mon profil de température-->ui/umax = ui/100 mais ecla ne fonctionne pas

Le calcul vectoriel est probablement plus rapide que faire des boucles sous python. D'ailleurs pourquoi l'avoir modifié ?  
Sur le coup, tu mélanges des parties vectorielles (les conditions aux bords) avec le schéma numérique calculé en chaque i,j.
De plus, tu t'es trompé en convertissant le terme advectif :

Autre point : il y a un problème d'indentation dans ton code (ou du moins pour l'affichage sur le forum) ce qui le rend ambigu (surtout pour du python).
Comme on est dans les considérations de codage, tu pourrais ne mettre que la définition des conditions aux bords dans le test m*dt<0.1, le schéma numérique étant indépendant du temps.
 

Tu devrais peut-être mettre plus en avant, dans ton post, la partie affichage (problème avec l'utilisation de manager.canvas.draw de matplotlib et gobject.idle_add) et reléguer au second plan la partie numérique. Il te reste quoi comme problème d'affichage ? C'est toujours noir ou c'est juste que tu ne vois pas l'évolution du calcul ?
Le titre n'est d'ailleurs pas du tout clair.
 
À part ça, l'image que tu obtiens ressemble bien à la solution stationnaire que tu devrais obtenir. En attendant de résoudre tes soucis d'affichage, tu devrais d'ailleurs diminuer la vitesse vo pour obtenir une solution stationnaire moins proche de l'état initial et ainsi mieux visualiser l'évolution de ta solution.

 
C'est corrigé, merci.
 

 
Oui en effet...je modifie ça
 

 
Toujours noir

Le code suivant fonctionne:
 

 
Par contre ça ne fait pas ce que je veux. Je n'ai pas des conditions aux limites qui changent avec le pas de temps. Manifestement il y a un problème avec ça. Je rappel que je veux pour ui[0,:] une rampe de température de 0 à 1 seconde (20°C à 100) puis maintient en température à 100 °C le reste deu temps

 
En regardant la doc de imshow, la palette de couleur est automatiquement adaptée au min/max des données fournies mais uniquement au moment de l'appel de cette méthode. Dans ton code, quand imshow est appelée, ui est constant et vaut 20.0 ce qui donne une palette de couleurs réduite à sa plus petite valeur (noir).
 
Pour éviter ça, tu peux spécifier les bornes de tes données :

 
PS: faut attendre un peu que la température monte pour voir quelquechose.

WOUOU ca marche. Merci encore à toi

Au passage, je te conseille *fortement* ( pour l'avoir testé ) de repasser à la version vectorielle de ton code, tu devrais voir la différence de performance (tu peux aussi n'afficher le résultat que toutes les n itérations)
 
 

Oui c'est franchement bcp bcp plus rapide