Résolution numérique d'équa diff

[Yoghourt]

Bonjour,

Le transport de chaleur par conduction dans un bout de paroi hétérogène s'exprime:

sachant:

- la conductivité thermique n'est pas spatialement constante
- Masse volumique et chaleur spécifique sont des constantes

Contrairement à ce qui est discuté ici, le but n'est pas de faire un calcul instationnaire 1D des transports couplés chaleur/humidité, mais de faire un calcul stationnaire 2D. (pour le 3D, on verra ensuite)

L'équation de transport se simplifie donc en:

sachant qu'on travaille en 2D.

Connaissant température intérieure, composition de la paroi, et absence de flux de chaleur sur les bords de la tranche (), quelle est la méthode numérique à utiliser?

Je me souviens avoir fait un TP de physique sur ordi à programmer ça dans le cas dynamique, pour une paroi homogène (échauffement d'une plaque métallique), quand j'étais taupin. Maintenant que je suis vieux et bête, j'enrage de ne plus savoir comment on fait...

But non avoué: calculer et représenter des ponts thermiques
-----

[invite6b1e2c2e]

Salut,

Je ne comprends pas bien le problème. Tu connais tout et tu veux calculer un truc ? Ca me paraît étrange, il y a quelque chose que je n'ai pas du saisir.

Je me demande si tu ne fais pas un problème inverse: Tu connais la température sur un sous domaine, et tu te demandes comment calculer la température dans tout le domaine. Si c'est ça, c'est un problème difficile, et les mots clés que tu cherches sont impedance tomography, inverse problem for elliptic problems, ce genre de choses.

__
rvz, qui répond probablement à coté

[Yoghourt]

Je connais la composition d'un mur hétérogène, la température (homogène) en x=0 et en x=épaisseur totale du mur, l'absence de flux en y=-largeur/2 et y=+largeur/2.
Je cherche typiquement le champ de températures.

[Yoghourt]

En cherchant ce que scilab sait faire en résolution d'EDP parabolique, je suis tombé sur ça:
http://www.freefem.org/ff++/index.htm

Yapluka comprendre un peu et tester (aspirine, mon amie...)

[A voir en vidéo sur Futura]

[Rincevent]

salut,

ton équation n'est pas parabolique mais elliptique si tu t'intéresses à la situation stationnaire. Et y'a pas une méthode numérique, mais plusieurs... Tu peux par exemple la réécrire sous la forme d'une équation de Poisson avec une source dépendant de la fonction inconnue et résoudre par itération.

[Yoghourt]

A force de penser à l'EDP parabolique 1D en parallele de mes cogitations sur l'elliptique 2D, je commence à m'emmeler les pinceaux...

Si je développe selon ton conseil, l'équation de transfert stationnaire de chaleur en milieu inhomogène sans source de chaleur ni convection devient alors:


A resoudre dans le cas présent dans un espace , avec des conditions aux limites de type

• (condition de Dirichlet) ou bien
• , où est un vecteur normal à (condition de Neumann).

En l'occurence, est une section de mur, c'est à dire un rectange ExL (épaisseur en x). Je passe en coordonnées cartésiennes pour ré-exprimer complètement le problème initial:

• Equation de transport
  
  
• Température constante à l'intérieur
  
• Température constante à l'extérieur
  
• pas de flux de chaleur sur les tranches de la section de la paroi
  
• la composition de la paroi est connue, mais n'est a priori pas homogène ni ne répond à une loi simple
  

Jusqu'ici ça va.

Dans le manuel freefem++, je constate qu'une ré-écriture sous forme d'intégrale est effectuée. Mais je patauge complet pour comprendre les exemples et comment transposer au cas qui m'intéresse.

De mes recherches sur le net, l'algo multigrid semble jouir d'une excellente réputation à résoudre les problèmes elliptiques. Dans "numerical recipees in C", il y a même le code de cet algo.
Mais le solveur n'est pas tout. Encore faut-il:

• faire le mesh
• comprendre comment injecter le membre de droite dans le solveur
• traduire le problème dans le formalisme du solveur
• traduire les résultats du solveur dans le mesh et les représenter
• valider tout le bousin

[invitebf65f07b]

bonsoir,

freefem me semble un bon moyen de résoudre ton problème. c'est gratuit et tout à fait apte à traiter ce genre de problème.

D'ailleurs, il y a même un exemple complet sur l'équation de la chaleur en stationnaire qui est traité dans la documentation en ligne (p.24).

De ce que je vois, il est plus simple de traiter un cas où, même si elle n'est pas uniforme dans tout le milieu, la conductivité soit constante par morceaux. Ceci me semble assez compatible avec l'application visée, donc pas de problème de ce côté là.

Enfin, pour ce qui est de la "ré-écriture sous forme d'intégrale" comme tu dis, va voir du côté de ce qu'on appelle la formulation variationnelle ou formulation faible, c'est de cela qu'il s'agit.

[Yoghourt]

'soir,
Voici la modélisation que j'ai faite:

mesh Th=square(30.2,32,[-0.1+(0.1+20+10+0.1)*x, -16+32*y], flags=1);
fespace Vh(Th,P1); Vh u,v,w;
Vh lambda=0.1/0.06*(x<0)
         + 0.04*(x>=0)*(x<10)
         +(80-0.04)*(x>=0)*(x<0.2)*(y>=-3)*(y<3)
	 +(80-0.04)*(x>=0.2)*(x<10)*(y>=-0.4)*(y<0.4)
         + 1.2*(x>=10)*(x<30)
         + (80-1.2)*(x>=10)*(x<15)*(y>=-0.4)*(y<0.4)
         + 0.1/0.11*(x>=30);
Th=adaptmesh(Th,lambda, err=0.005);
plot(Th,wait=1);
solve a(u,v)= int2d(Th)(lambda*(dx(u)*dx(v)+dy(u)*dy(v)))
                +on(4,u=-10)+on(2,u=16);
real[int] isoValues(53);
for(int i=0;i<53;i++)
  isoValues(i) = i/2-10;
plot(Th,u,wait=1,fill=true, viso=isoValues, value=true);
plot(Th,u,wait=1,fill=true, viso=isoValues);

Ci-joint le mesh. Le résultat du champ de températures est dans ce post.

[Yoghourt]

Pour l'instant, je ne sait pas trop comment calculer avec ledit logiciel le flux de chaleur

[emmanuel30]

Bonjour, et merci à toi pour ce travail.

Bon, s'il y a une erreur dans les formules, ce n'est pas moi qui les trouverai, j'ai du louper quelques cours dans ma jeunesse.

J'avoue qu'avec les shémas, j'ai du mal à quantifier les pertes.

Il serait peut être intéressant, d'avoir les mêmes schémas avec hi et he = 0

Ce serait pour visualiser la différence qu’il y a quand on les néglige.

Le but final étant de savoir si on doit tout mettre en œuvre pour limiter les ponts thermiques dus aux vis ou si l’on peu les ignorer.

[Philou67]

Sauf erreur, h_i et h_e représentent les résistances superficielles correspondantes aux fines couches d'air de part et d'autre de la paroi... En quoi le fait de les annuler permettrait de vérifier si les ponts thermiques peuvent être négligés ?

[emmanuel30]

Salut philou, je me suis trompé de post, je voulais répondre sur le fil intitulé Ponts thermiques des vis en acier

Pourrais tu déplacer ma réponse au bon endroit car elle n'a rien à faire ici, si ce n'est, prêter à confusion.

En lisant le début de l'autre post tu devrais avoir une réponse aux questions que je me pose.

Pour résumer si on néglige hi et he le pont thermique sur une porte en ALU est monstrueux et si on les prend en compte ce pont devient raisonnable, en est il de même avec les vis qui traversent tout l'isolant ?

[Philou67]

Bah, on est dans un forum de Math, et je n'ai aucun pouvoir, mais je vais demander à un collègue

[Yoghourt]

Pour l'instant, je ne sait pas trop comment calculer avec ledit logiciel le flux de chaleur

Dis donc dis donc, j'étais un peu intelligent il y a une douzaine d'années. Mes engrenages sont bien moins neufs maintenant formellement, il manque un signe -, p'tit rigolo

Dire que j'avais oublié cette petite prise de main de freefem++ et l'avait "jeté le bébé avec l'eau du bain"...

Voici un petit code de test d'un cas calculable facilement à la main, dont le fameux calcul de flux que je ne voyais pas comment faire à l'époque.

int bordExt=1;
int bordInt=2;
int bordAutre=3;

border bas(t=-0.1,0.1){x=t;y=-0.2; label=bordAutre;}
border droit(t=-0.2,0.2){x=0.1;y=t; label=bordInt;}
border haut(t=0.1,-0.1){x=t;y=0.2;label=bordAutre;}
border gauche(t=0.2,-0.2){x=-0.1;y=t; label=bordExt;}

int nbnoeuds=64;
plot(bas(nbnoeuds) + droit(nbnoeuds) + haut(nbnoeuds) + gauche(nbnoeuds), wait=true);

mesh Th1 = buildmesh(bas(nbnoeuds) + droit(nbnoeuds) + haut(nbnoeuds) + gauche(nbnoeuds));
plot(Th1, wait=true);

func lambda = 0.04*(x<0)
			+ 1.2*(x>=0);
mesh Th2 = adaptmesh(Th1,lambda,err=0.0001);
plot(Th2, wait=true);


fespace Vh(Th2, P1);
Vh T, v;
real Text = -2.0;
real Tint = 18.0;

solve thermiqueStationnaire(T, v, eps=1e-6)
    = int2d(Th2)(
          lambda*(
              dx(T)*dx(v)
            + dy(T)*dy(v)
        )
    )
    +on(bordExt, T=Text)
    +on(bordInt, T=Tint);

plot(T, wait=true, value=true, fill=true);
Vh jx = -lambda*dx(T);
Vh jy = -lambda*dy(T);
plot(jx, wait=true, value=true, fill=true);
plot(jy, wait=true, value=true, fill=true);
// attendu : jx = -20/(0.1/0.04+1.2/0.04) = -7.7419355 W/m² et jy = 0
cout << "jx(-0.05,0)=" << jx(-0.05,0) << endl;
cout << "jx(+0.05,0)=" << jx(+0.05,0) << endl;
cout << "jy(-0.05,0)=" << jy(-0.05,0) << endl;
cout << "jy(+0.05,0)=" << jy(+0.05,0) << endl;
real J = int1d(Th2,bordInt)(-lambda*dx(T));
cout << "jx moyen = " << J/0.4 << endl;

Résultat :
jx(-0.05,0)=-7.74281
jx(+0.05,0)=-7.74527
jx moyen (+0.1, y)= -7.75339
On est à 0,15% d'erreur vis à vis du calcul analytique. Remarquable!

jy(-0.05,0)=0.0051972
jy(+0.05,0)=-0.00261852
On attendait 0, la précision du calcul numérique est là aussi très bonne.

Temps de calcul? 3 secondes! Mon pauvre programme scilab en résolution explicite et relaxation adaptative est totalement enterré.

Deux petits hics pas bloquants :
- je n'ai pas trouvé comment forcer la main de freefem++ pour obtenir une plus grande précision. J'ai joué avec le nombre initial de noeuds, le facteur de précision d'adaptmesh, l'epsilon du solveur, le type de mesh P2 (P1dc crashe le script)... Rien de flagrant, quelques régressions, rien compris!
- pour l'affichage, j'ai essayé de jouer avec hsv et viso pour obliger freefem a utiliser une palette style jetcolormap pour aller du rouge=chaud au bleu=froid en passant par orange, jaune, vert, cyan. Rien compris au biniou, j'ai abandonné.

Bon ben youpi, yapluka!
Y.

P.S. je laisse tomber l'environnement graphique freefem++-cs : ça plante brutalement plus souvent qu'à son tour à l'exécution d'un script quelconque. Beuark!

[Yoghourt]

1, 2, 3, voilà ce qui arriva...
Non pas un bikini rouge et jaune à p'tits pois, mais un bon mal de crâne avec freeFEM++. Y'a-t-il un expert dans la salle?

Je cherche à migrer le calcul stationnaire inhomogène décrit ici de mon "vieux" programme scilab vers un truc plus fiable, plus rapide et tout et tout.
En théorie, ça ne devrait pas poser de problème avec le défrichage ci-dessus.
En pratique, je me prends poteau sur poteau :
- pour que le mesh respecte les interfaces de changement de conductivité thermique comme dans l'exemple heat exchanger de la doc, ça demanderait une décomposition folle à écrire à la main
- adaptmesh guidé par lambda(x,y) s'avère inutilisable. Soit des zones sont maillées trop gros, soit il perd les pédales (message 'n vertices losts, sorry') quand ça devient fin et fait des trous. Et en plus le résultat final n'est pas stable
- pour empiler 3x3 scenarii, quelle galère! Je n'ai pas trouvé comment poser une description de base puis la customiser. La seule solution que j'ai trouvée, c'est en sens inverse : poser les cas particulier en 1er puis descendre dans les couches vers le cas général, un peu comme un front to back sur un schéma power point
- Quand je définis un contour fluv et fluh aux travers desquels je veux calculer le flux de chaleur int1d(Tmesh,fluv)(-lambda*dx(T)) et int1d(Tmesh,fluh)(-lambda*dy(T)) c'est à dire les déperditions de l'espace intérieur via le mur et le plancher, j'obtiens 0. Alors que si je fais un plot du champ vectoriel de flux de chaleur ou si je sors à l'écran le flux en un point donné, ça semble ok. La seule façon que j'ai trouvée pour l'instant, quand ça veut bien, est de calculer ce flux sur un bout de bord extérieur du domaine utilisé dans buildmesh(). Or en bordure extérieure on est justement sur une discontinuité, donc pour un résultat réaliste et fiable du calcul de flux, bof bof
- Si je passe le solveur en P1dc ou en P2, ça fait n'importe quoi
- A part en réduisant la taille des mailles mais pas trop (sinon à nouveau le message "vertices lost sorry"), les contrôle de précision du solveur err= ou errg= ne jouent en rien sur la stabilité du champ de température obtenu. Le flux -probablement mal- calculé ne tend pas vers une valeur fixe, monte lentement au lieu de se stabiliser.
- La doc est surement très bien pour un mathématicien spécialiste de ces trucs, mais pour un développeur utilisateur, difficile de s'y retrouver voire impossible (par exemple carte de couleur hsv=... d'un plot incompréhensible après moult essais pour avoir des couleurs style jetcolormap() et pas leur truc par défaut contre-intuitif)

Description, mesh, calcul, visualisation... Après des heures passées à batailler et une dizaine de variantes pour tenter ceci ou cela, je ne suis pas hyper optimiste. S'il y a des connaisseurs, un coup de main ne serait pas de refus

Chuss,
Y.