Accueil

Diffusion bidimensionnelle

Simulation numérique de la diffusion de la chaleur

1 Diffusion unidimensionnelle

1 Diffusion unidimensionnelle

2 Diffusion bidimensionnelle

3 Diffusion tridimensionnelle

4 Application aux ponts thermiques

5 Simulations dynamiques

6 Changement d'état

7 Déconvolution numérique

8 Défloutage photo


  Bien qu'elle soit traitée ici sous l'angle de la conduction thermique, la diffusion est un phénomène  général, impliqué en chimie, en métallurgie, en fabrication micro électronique ...
  Elle se développe en milieu tridimensionnel, mais se traite souvent sur un mode unidimensionnel : chauffage d'un mur composite ou non, d'une plaque ...

  La diffusion de la chaleur et celle des espèces chimiques dans les solides ou les gels sont régis par la même équation aux dérivées partielles :

les solutions analytiques et les simulations numériques développées pour la chaleur seront donc applicables à la diffusion des espèces chimiques dans les solides en remplaçant T la température par la concentration.

  Quand D est une diffusivité thermique, elle est donnée par :

  Le flux de chaleur traversant un plan perpendiculaire à l'axe x est donné par :


1.1 Différences finies appliquées à la diffusion unidimensionnelle

  Lorsque les conditions expérimentales conduisent à des isothermes plans et parallèles, les flux thermiques selon y et z sont nuls, le problème relève d'un traitement unidimensionnel et l'équation de la chaleur se réduit à :

(1.1)

  Ce chapitre concerne donc la diffusion de la chaleur dans un milieu semi infini mais aussi dans des parallélépipèdes ou des cylindres isolés sur leurs parois latérales.

  Considérons une barre conductrice, de grande longueur, parfaitement isolée, à température initiale nulle.
  Au temps t = 0 chauffons sa partie centrale pendant un très court instant. Il s'agit d'un cas classique : la température se répartira selon des gaussiennes de surfaces constantes et de sigma croissants.

  Le problème sera traité par la méthode des différences finies, FDTD ( Finite Difference Time Domain ), appliquée sans recours extérieur.

  Pour cela, la barre est subdivisée en cellules de longueurs égales.Au temps zéro la cellules centrale est chauffée en lui apportant, par exemple, 256 unités de chaleur.

0 0 0 0 256 0 0 0 0

  Simulons maintenant la dispersion de cette quantité de chaleur. Admettons pour cela qu'après un temps t une cellule quelconque aura perdu 1/4 de sa chaleur vers la cellule située à sa gauche et un autre 1/4 vers la cellule située à sa droite. Il lui restera donc la moitié de sa chaleur initiale :

1/4 1/2 1/4

  Ainsi, dans la barre, après ce temps  t la nouvelle répartition des quantités de chaleur deviendra :

      64 128 64      

  4 étapes seront nécessaire pour calculer la répartition après un nouveau laps de temps t :

- 1- les 64 unités de la 4ème colonne se partageront en 16 , 32 et 16 unités

      64 128 64      
    1/4 1/2 1/4        
    16 32 16        

- 2- les 128 unités de la 5ème colonne se partageront  en 32 , 64 et 32 unités

      64 128 64      
      1/4 1/2 1/4      
      32 64 32      

- 3- les 64 unités de la 6ème colonne se partageront en 16 , 32 et 16 unités

      64 128 64      
        1/4 1/2 1/4    
        16 32 16    

- 4 - La somme de ces 3 contributions donneront les quantités de chaleur au temps 2 t

    16 32 16        
      32 64 32      
        16 32 16    
    16 64 96 64 16    

  En poursuivant, pour 3t :

  4 24 60 80 60 24 4  

  Puis, pour 4t

1 8 28 56 70 56 28 8 1

Qui sont les coefficients du binôme, ceux de la huitième ligne du triangle de Pascal.

...

  Les 3 pondérations  { 1/4, 1/2, 1/4 }assurent déjà des simulations satisfaisantes, elles seront plus fidèles en utilisant ces trois autres, { 1/6, 4/6, 1/6 }, selon la méthode présentée par Jean LEGRAS dans :
    " TECHNIQUES DE RÉSOLUTION DES ÉQUATIONS AUX DÉRIVES PARTIELLES " , pages 33 et 34 .

 Au temps t, le développement en série de Taylor exprime les températures Tx+h,t et Tx-h,t  en fonction de la température et de ses dérivées en x 

d’où, par addition, une approximation où seules les dérivées d'ordres pairs supérieurs à 2 sont négligées :

(1.2)

  La température d'un élément aux temps ( t + k ) s'exprime aussi en fonction de la température au temps t :

d'où une approximation négligeant les dérivées d'ordres supérieurs à 2 :

(1.3)

  En reportant les approximations (1.2) et (1.3) dans l’équation de la diffusion (1.1) il vient :

(1.4)

  Un choix judicieux des pas en x et en t annulera la contribution des 2 dérivées.

  Dériver l’équation de la diffusion (1.1) par rapport au temps ou deux fois par rapport à x donne :

et

d’où

(1.4) devient alors:

(1.5)

  En prenant :

(1.6)

(1.5) se réduit à :

(1.7)

 Ces 3 pondérations optimales s'imposeront seulement pour les modèles simulés par un petit nombre de cellules ou, en multidimensionnel, ceux comportant des détails représentés par quelques cellules seulement.

 En revanche, un maillage fin autorise l'emploi de pondérations différentes sous réserve de respecter :

        H >= 2

(1.8)

  Connaissant les températures dans les différentes cellules, le flux thermique entre deux cellules adjacentes est donné par :

(1.9)

  Le flux thermique s'annule quand le gradient s'annule, par exemple au contact d'un isolant idéal. Avec les différences finies le flux en x s'annulera quand Tx,t =  Tx+h,t


1.2 Technique opératoire

  Pour une simulation portant sur N cellules on utilisera une variable dimensionnée T(2,N).
  Les températures initiales seront introduites dans T(1,I) .
  Pour calculer la répartition des températures au temps k, les transferts thermiques entre cellules seront simulés en appliquant (1.7) à chacune des cellules, les résultats étant provisoirement stockés dans T(2.I) .A la fin de ce premier cycle toutes ces valeurs sont recopiées dans T(1.I) avant d'attaquer un second cycle qui, par le même processus, conduira à la répartition au temps 2k. ...

  Périodiquement, un tracé du profil des températures atteintes permettra d'évaluer l'avancement du travail.

  En simulation unidimensionnelle, indépendamment de la prise en compte des conditions aux limites, l'essentiel des opérations se limite à 3 boucles :

     For J = 1 To Duree
          For I = 1 To Longueur - 1
               T(2, I) = (T(1, I - 1) + 4 * T(1, I) + T(1, I + 1))
          Next I
           For I = 1 To Longueur
                T(1, I) = T(2, I) / 6
          Next I
     Next J

1.3 Conditions aux limites

1.3.1 Température imposée

  Aux points où la température est fixe, nulle ou variable, le programme inscrira continuellement la valeur imposée.

1.3.2 Extrémité isolée

  Par définition, sur une extrémité parfaitement isolée le flux thermique est nul. Il en serait de même si cette extrémité était en contact thermique parfait avec un objet symétrique soumis au même régime thermique : pour simuler la présence de cet objet virtuel il suffira de constituer une cellule supplémentaire, extérieure à l'objet, à laquelle on attribuera systématiquement la température de l'avant dernière cellule (son symétrique par rapport à la limite).

  Une technique équivalente consiste à appliquer, à chaque cycle :

 

ou

selon que la paroi isolante se trouve à gauche ou à droite.

1.3.3 Température fonction du temps

  Ce sera souvent le cas dans les applications concrètes : les températures de surface seront réajustées après chaque cycle de simulation.
  Le comportement d'un mur soumis à des variations sinusoïdales de températures, indépendantes sur chacune de ses deux faces ses deux faces, est traité en quatrième page. Passée la période transitoire, la fiabilité du procédé se confirme.


1.4 Comparaison avec une solution analytique

1.4.1 Solutions analytiques de la diffusion dans un milieu semi infini

  Deux cas classiques :

- Température de surface constante 

  Condition initiale : Tx,t = 0 pour x>0 et t<0
  Condition en surface : T0,t= 1 pour t>0

  L'équation différentielle (1,1) admet alors la solution

où, en posant :   

(1.10)

  La fonction erreur, erf(x), et la fonction erreur complémentaire, erfc(x) = 1 - erf(x), sont calculées par développement en série :

Afficher une table de la fonction erreur complémentaire

- Quantité de chaleur constante

  Condition initiale : Tx,t = 0 pour x>0 et t<0 et injection d'une unique quantité de chaleur en x = 0 au temps t = 0

  La température en fonction de la profondeur et du temps évolue comme :

( 1.11 )

1.4.2 Simulation 

  Soit une diffusion unidimensionnelle (isothermes plans) simulée par un alignement d'une trantaine de cellules.
  Tous les contenus sont initialement nuls,.Seule la cellule centrale reçoit une quantité de chaleur à t = 0.

  En adoptant H = 4 , soit les trois pondérations [ 1/4 , 1/2 , 1/4 ] , les écarts entre valeurs théoriques et valeurs simulées seraient déjà acceptables. Après 32 cycles de simulations , les différences secondes s'annulent en x = -4 et x = +4 soit un écart type de 6

  Mais la simulation sera optimisée en appliquant (1.7) c'est à dire en adoptant H = 6 et en portant le nombre de cycles à 48 :

  Les écarts deviennent 100 fois moindres, avec des erreurs relatives de quelques 10-5 dans les zones les plus significatives, garantissant ainsi une fiabilité suffisante pour la plupart des applications.

  Ceci suppose que la simulation s'effectue dans un matériaux unique. Lorsque la diffusion concernera divers matériaux en contact thermique, il deviendra indispensable de leur attribuer des valeurs de H , supérieures, appropriée aux diverses diffusivités.

  Cette méthode de simulation restera applicable à la plupart des structures comportant plusieurs matériaux, bien que les désaccords soient alors centuplés. Dans tous les cas la précision sera dépendante de la largeur des pas, limitée par la durées des calculs.


1.5 Application aux composites

  L'équation de la chaleur s'appliquera dans chaque milieu avec ses propres paramètres

  Les flux thermiques étant égaux de part et d'autre de l'interface, le rapport des gradients y sera proportionnel au rapport des diffusivités.

  En l'absence de résistance thermique entre les deux constituants les 2 températures d'interface resteront constamment identiques.

  Ayant à traiter divers matériaux se distinguant par leurs coefficients de diffusion mais aussi par leurs chaleurs spécifiques et leurs masses volumiques, l'équation différentielle unidimensionnelle :

devient :

dont les paramètres sont liés par la relation :

D : diffusivité thermique ( m2 / s )
K : conductivité thermique ( W / m deg)
  : masse volumique ( kg / m3 )
C : chaleur spécifique (J /kg deg )

  A l'interface, le flux thermique sortant et le flux entrant, exprimés en W/m², sont identiques et liés au gradient thermique. 


1.5.1 Analogie électrique

  Le développement du numérique à supplanté les autres techniques de simulation dont les analogies électriques qui gardent cependant l'avantage de faciliter la description des processus.
  La classique équation des télégraphistes exprime l'évolution de la tension sur une ligne de transmission :

où L, C, R et G sont respectivement l'inductance, la capacité, la résistance et la perditance par unité de longueur.

  Si l'inductance est nulle et l'isolation parfaite, L et G sont nuls, l'équation se réduit à :

analogue à l'équation unidimensionnelle de la chaleur

  R et C représentent respectivement un résistance par unité de longueur et une capacitance par unité de longueur, leur produit a même dimension que l'inverse d'une diffusivité, c'est à dire : unité de temps par unité de surface.
  Plus prosaïquement, la capacité pour un matériaux d'emmagasiner de la chaleur est assimilable à la faculté d'un condensateur d'accumuler des charges électriques. Parallèlement, la diffusivité de 
ce matériau, qui exprime son aptitude à transmettre la chaleur, est équivalente à l'admittance d'un conducteur électrique, l'admittance étant l'inverse de la résistance.
  Un barre conductrice de la chaleur est donc assimilable à une ligne de transmission continue, schématisée sous forme d'un circuit équivalent, constitué d'éléments discrets représentatifs d'une unité de longueur de ligne.

  Chaque cellule est constituée de 3 éléments, selon deux types de cellules quadripolaires envisageables :

formant une chaîne de quadripôles soit en Pi soit en T.

  Dans la simulation numérique de la diffusion unidimensionnelle en milieu homogène, l'évolution de la température d'une cellule est établi en fonction des tensions antérieures de cette cellule et des deux cellules adjacentes.
  En présence de deux ou plusieurs milieux différents, chacun sera représenté par plusieurs cellules mais une cellule limitrophe aura une voisine différente aussi bien à sa gauche qu'à sa droite. Pour parer à toutes les éventualités il convient donc de prévoir des structures comportant des quadripôles de caractéristiques différentes.

 

Cellules en Pi

Cellules en T

  Pendant un temps k, ces deux courants chargeront, ou déchargeront, les deux condensateurs en parallèle à gauche ou le condensateur Cx à droite. Les  tensions résultantes au temps Ex,t+k seront respectivement :

1.5.2 Programmation

  Appliquée à la simulation de la diffusion de la chaleur, les tensions deviennent des températures et 1/RC correspond à D, la diffusivité.

  Dans une application, chaque constituant sera décrit par une suite de cellules. Avant exécution de la simulation  les caractéristiques de chacune des cellules seront mémorisées . Quelles valeurs attribuer aux cellules d'interface ?

  Le modèle en Pi, interface entre deux cellules, ne pose aucune difficulté de programmation. En revanche, dans le modèle en T, l'existence de cellule mixtes  nécessiterait un traitement spécifique pour les cellules limitrophes et leurs deux voisines.
  Bien qu'elle présente l'avantage de fournir directement les températures centrales des cellules, l'option en Pi lui sera préférée.


1.5.3 Application

  Soit un mur constitué de 2 matériaux, maçonnerie et isolant en contact thermique parfait. La température initiale est nulle et pour t>0 ses deux faces seront maintenues à températures constantes :
 -  par un air extérieur parfaitement brassé à -25°C à gauche
 -  par un chauffage thermostaté à +25°C à droite.

  A l'équilibre les températures seront constantes, la dérivée par rapport au temps sera donc nulle, en unidimensionnel l'équation différentielle se réduit donc à :

  Pour que la dérivée seconde soit nulle la dérivée est constante et les solutions sont nécessairement deux droites. Les flux thermiques étant devenus égaux dans les deux matériaux le rapport des pentes est imposé par le rapport des diffusivités, tout en respectant les températures aux limites.

  La simulations proprement dite travail avec deux types de nombres  :
 - T(I)  où seront stockées les températures successives, dont les températures initiales mentionnées dans le tableau suivant 
 - L(I) contient les conductivités thermiques
 - C(I) contient les produits  masse volumique x chaleur spécifique.

  Conventionnellement, C(I) est nul pour les cellules d'air.

  D'où, pour les paramètres implicites du programme ci-dessous :

I -1 0 1 2 3 4 5 6 7 ... 13 14 15 16 17 18 ... 23 24 25 26 27 28 28 30 31
T(I) -25 -25 -25 -25 -25 -25 -25 0 0 ... 0 0 0 0 0 0 ... 0 0 25 25 25 25 25 25 25
L(I)               2 2 ... 2 2 2 0.5 0.5 0.5 ... 0.5 0.5              
C(I)   0 0 0 0 0 0 1 1   1 1 1 1 1 1   1 1 0 0 0 0 0 0 0

  4 éventualités correspondent aux 4 instructions " If ".
  Les nouvelles températures sont d'abord enregistrées dans le relais " R(I) ", puis recopiées dans T(I) en fin de cycle.
  A l'intérieur du mur, toutes les cellules sont soumises au même traitement qu'elles soient ou non en contact avec une cellule de nature différente. Ceci est réalisé en adoptant le calcul établi pour le schéma électrique équivalent en forme de Pi .

  Ce programme de simulation affiche les températures entre cellules après 10 000 cycles de simulation , ainsi que les écarts aux températures théoriques.
  Le profil initial de température figure en bleu, les pointillés bleus correspondent à des étapes intermédiaires et le profil final, pour lequel l'équilibre thermique est pratiquement atteint, est porté en rouge.

Télécharger l'exécutable

  Pour observer les profils de température en fonction des paramètres ou avant d'entreprendre un programme plus élaboré dans un autre langage.
 

Télécharger Mur_Isol_exe

 

  Ce fichier est compressé mais auto- extractible. Après appel, l'acceptation du chargement ouvrira cette fenêtre.

  En cliquant  "Décompresser"  vous l'installerez dans votre dossier "Program Files" (ou dans celui que vous choisirez).

  Le dossier Mur_Iso_Exe.exe contiendra :
- l'exécutable
- un raccourci (orientant vers C:\Program Files )
- 2 fichiers indispensables si Visual Basic 5 n'est pas installé sur votre disque dur.

  Si ce télé-chargement vous posait quelques problèmes, consultez d'abord le mode opératoire détaillé.

1.5 Chauffage à flux constant

  Dans un autre cas classique, la surface libre du milieu semi infini est chauffée par un flux constant de chaleur, une résistance alimentée à puissance constante par exemple.

1.5.1 Solution analytique de la diffusion à flux constant

  Sous l'effet d'un flux de chaleur F par unité de surface et unité de temps, la première tranche d'épaisseur dx verra sa température croître d'une quantité dT/Q, Q étant la chaleur spécifique du matériau.

  Problème classique présenté dans les livres et les sites didactiques : la température Tx,t à la profondeur x sera au temps t : 

  La température superficielle croît donc avec la racine du temps sans autre limite que la fusion de l'ensemble. En profondeur, elles correspondent aux points rouges portés sur le dernier profil de la prochaine figure.

1.5.2 Simulation numérique la diffusion à flux constant

  Ici le problème admet une solution analytique, son traitement par simulation numérique permet de vérifier sa fiabilité avant de l'appliquer aux cas concrets.

  L'ensemble des cellules est remis à zéro puis, la température de la cellule superficielle T(1,0) est rehaussée avant chaque cycle de simulation. 

  Dans la seconde ligne de la boucle, la contribution T(1, 1) de la deuxième cellule est doublée afin de tenir compte des apports de la cellule virtuelle d'abscisse J = -1 .

  Après programmation de la solution analytique et de la présentation :

  Les valeurs calculés ( points rouges ) s'insèrent correctement dans la courbe établie par simulation numérique, ce qui  valide le procédé pour les applications sans solution analytique.


Extraire la page pour l'enregistrer ou l'imprimer

//