Discussion :

[FR119492]

Salut!

On a alors la relation de récurrence

Tu as été un pas trop loin: il fallait discrétiser selon r, mais pas selon t. Ainsi tu aurais obtenu un système différentiel ordinaire avec t comme variable indépendante et Theta1(t), Theta2(t), etc. comme variables d'état. En fait, en écrivant ta relation de récurrence, tu as, sans t'en rendre compte, appliqué la méthode d'Euler du 1er ordre, méthode dont on dit: "Admirez mais n'appliquez pas!"

j'ai Simulink mais je ne me souviens plus comment faire une dérivée

Simulink est un logiciel bien adapté pour ton problème, mais n'utilise jamais le bloc "dérivation"!

Jean-Marc Blanc

[Franzis]

Bonsoir;

Ok j'ai compris le principe: j'ai discrétisé suivant r seulement j'obtiens donc une équation dTheta/dt en fonction de Theta(i), où i est la position suivant r du type V'=f(t,V) que je résoud avec ode45 de Matlab par exemple, où V est le vecteur colonne de la température pour chaque position différente (i=1 à r=dR, i=2 à r=2dr... jusqu'à i=M à r=R) et la fonction 'f' c'est celle que j'écris avec l'équation:



Le soucis que je rencontre c'est dans l'écriture de la fonction si je la fais entre i=1 et i=M, il y a un problème aux bornes car Matlab doit trouver Theta(0) qui n'existe pas... c'est la valeur de la température en r=0 mais comment l'écrire?

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
function dX=simu_chaleur(t,X)
 
a=lambd*M*M/(Rho*Cp*R*R);
 
dX=[];             %Pré-remplissage
 
for i=2:1:M
     dX(i)=a*(-2*X(i)+(1-1/(2*i))*X(i-1)+(1+1/(2*i))*X(i+1));
end

Erreurs renvoyées:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
 
??? Input argument "X" is undefined.
 
Error in ==> simu_chaleur at 19
     dX(i)=a*(-2*X(i)+(1-1/(2*i))*X(i-1)+(1+1/(2*i))*X(i+1));
 
??? Attempted to access X(2); index out of bounds because numel(X)=1.

Et puis il dit qu'il ne connait pas X alors que c'est ce qu'il doit trouver... Soit je m'y prend mal soit je n'ai pas bien compris. Dans tous les cas merci encore de m'éclairer.

[FR119492]

Salut!

si je la fais entre i=1 et i=M

En effet, tes rayons partent de 0.0 (réel) et tes indices de 1 (entier). Tu as donc r_i=(i-1)*DR/Rmax .
Jean-Marc Blanc

[Franzis]

Donc pour récapituler:

j'ai une ligne des espaces de i=1 (r=0), i=2 (r=Dr)... jusqu'à i=M+1 (r=R) et un vecteur temps de t=0 à temps_final.

J'ai aussi

• un vecteur colonne X à M+1 lignes qui est le vecteur des variables d'état: les températures Théta pour chaque position: Théta1=Théta(r=0)...
• un vecteur colonne à M+1 lignes qui est le vecteur définissant les conditions initiales.

et une fonction qui définie dX, le solveur de Matlab ode45 donnera alors une matrice où les lignes correspondent aux positions dans l'espace et les colonnes le temps. Jusque là OK sur le principe?

Quand j'écris ma fonction, pour i=1 Théta1(r=0) est connu, puis j'écris dX(pour i=2) en fonction de Theta1 (en autre) et ainsi de suite jusqu'à dX(i=M), avec en plus la condition limite en r=R (i=M+1) de gradient nul. Ce qui donne:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
 
dX=zeros(M+1,1);    % Pré-remplissage
dX(M+1)=0;          % condition limite en r=R gradiant nul
 
for i=2:1:M
    X(1,i)=E(i,1)'; % condition limite en r=0
    dX(i)=a*(-2*X(i)+(1-1/(2*i*(i-1)))*X(i-1)+(1+1/(2*(i-1)))*X(i+1));
end

Je sais pas comment écrire pour i=1. Est-ce que

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

X(1,i)=E(i,1)';

est correct? Est-ce qu'il faut que ça soit dans la boucle? Car avant que Matlab fasse le calcul X est un vecteur colonne puis après une matrice...

En tout cas merci pour ton temps et ta patience, je sens que je ne suis plus très loin de la solution.

[kwatz]

Bonjour,

Je suis dans un cas proche du tien Franzis, sauf que je suis en coordonnées cartésiennes (et sans génération de chaleur interne pour le moment).

J'ai discrétisé l'équation de la chaleur comme tu l'as posté le 3 mars, et je cherche la solution à l'instant t+dt grâce à la solution à l'instant t.

Cependant, dans mon cas, le pas de temps dt dépend de mon maillage et des caractéristiques matériaux, et je pense que dans le tiens aussi. L'expression de mon pas de temps est :

dt=dx²*rho*Cp/(2*lambda) avec dx qui représente la longueure d'une maille.
Ce pas de temps représente physiquement le temps de diffusion de la chaleur sur une distance dx.

En espérant t'avoir aidé!

[Franzis]

Bonjour,

Quelqu'un aurait une idée pour m'aider à écrire les conditions limites?

Merci d'avance

[Ehouarn]

Bonjour,

Franzis: Pour écrire la condition en r=R sur le flux nul, tu peux écrire que le polynôme d'interpolation de degré deux (puisque tu es en différences finies d'ordre 2), donc basé sur les points i=M-1,M et M+1 à une dérivée qui s'annule. Un exemple de cette méthode est donné ici.

Kwatz: Ton "pas de temps" n'est qu'une estimation de l'ordre de grandeur de l'échelle de temps de diffusion (n'est-ce pas?), ce qui ne signifie pas qu'il soit suffisamment petit pour obtenir une "bonne" solution numérique (on n'échappe généralement pas à une exploration "expérimentale" pour évaluer si le pas de temps utilisé était acceptable).

Ehouarn

[kwatz]

Pour les conditions aux limites, je ferais un bilan de chaleur tout simple.

Tu écris que la somme des flux entrant sur un élément est égale à la chaleur qu'il absorbe. Comme ca, si tu as des parois adiabatique, ou un coef de convection quelquepart, tu devrais t'en sortir facilement.



En ce qui concerne le temps de diffusion, c'est plus qu'une approximation. Dans le cas d'une résolution explicite, ce qui est le cas ici je crois, soit le pas de temps est suffisament petit, soit ce n'est pas le cas et les résultats divergent.

Dans le cas de mon étude, j'ai par contre remarqué que les résultats sont "maillage dépendant", c'est à dire que les mailles doivent être suffisament petites pour que la solution soit valable. Mais comme le pas de temps est lié à la taille des mailles... Ce que je fais pour vérifier mon modèle, c'est que les résultats avec ceux d'un logiciel de résolution par élément finis dans des conditions simples.

C'est vrai que pour obtenir une bonne simulation numérique, on n'échappe pas à une étude expérimentale, mais je pense que c'est plus pour ajuster les conditions aux frontières, non?

[Franzis]

Mes conditions aux limites je les connais déjà:

• en r=0, la température est connue quelque soit le temps
• en r=R le gradient est nul (condition de Neumann).

C'est sur la deuxième CL que porte ma question: comment l'écrire dans mon programme Matlab?

[Ehouarn]

Bonjour

Kwatz:

soit le pas de temps est suffisament petit, soit ce n'est pas le cas et les résultats divergent.

Non, ce n'est pas aussi simple. Avec un pas de temps trop grossier les résultats divergent (c.-à-d. tendent vers l'infini), sinon les résultats ne divergent pas (c.-à-d. restent bornés) et accessoirement plus le pas de temps sera petit, plus on tendra vers la solution du problème continu. On peut très bien ne pas diverger et être loin d'une approximation "acceptable" de la solution du problème continu.
C'est d'ailleurs exactement la même chose pour le pas en espace (la taille des mailles); si on en prend trop peu et qu'ils sont donc trop grands, on ne risque pas de représenter correctement des phénomènes de petite échelle (s'il y en a), ce qui ne veux pas forcément dire que le calcul divergera; on aura seulement une solution trop grossière (et conséquemment assez "fausse" du point de vue physique).

Franzis: Je n'utilise pas Matlab (ceci dit je doute que ta boucle fonctionne comme tu le voudrais; il faut normalement poser le système linéaire à résoudre puis le résoudre avec les outils appropriés); je peux seulement préciser la méthode.
Typiquement il y a une séparation entre la discrétisation à l'intérieur du domaine et sur les frontières: à l'intérieur du domaine on discrétise l'équation différentielle (ce qui donne des relations entre valeurs entre voisins comme tu l'a donné plus haut), mais pour déterminer les valeurs sur les frontières on utilise des équations issues de la discrétisation des conditions aux frontières, comme cela est détaillé dans le lien que j'ai donné plus haut.
En clair, pour la condition de Neumann, avec des différence finies d'ordre 2, sur les trois "derniers" points i=M-1,M et M+1 (i=M+1 est la frontière) espacés d'un pas (constant) h, la dérivée y' de la solution au point M+1 est telle que:
y'_{M+1} = (1/2h).(1.y_{M-1} -4.y_{M} +3.y_{M+1}).
où les y_{...} représentent les valeurs aux points i=M-1,M et M+1.
Ta condition de flux nul implique y'_{M+1} = 0 et tu a donc l'équation ci-dessus qui te donne la relation entre les valeur aux trois points, qu'il n'y a plus qu'à résoudre...

[FR119492]

Attention!
Lorsqu'on veut intégrer une équation de Fourier en discrétisant par rapport aux variables d'espaces par différences finies, puis en intégrant pas à pas le système différentiel obtenu, il faut prendre garde au fait que ce système peut être très raide (stiff): avec les méthodes explicites usuelles (Runge-Kutta ou autres), si le pas de temps est trop grand, le système diverge et si le pas est trop petit, le temps de calcul devient prohibitif. C'est pourquoi il faut avoir recours à des méthodes implicites.
Jean-Marc Blanc

[kwatz]

Bonjour,
je ne crois pas qu'une méthode implicite fonctionne ici car il y a une source de chaleur interne.

Et Ehouarn, je persiste mais d'après mon expérience (certe courte! ), le pas de temps est soit bon, soit mauvais (pour l'équation de la chaleur, je ne m'avancerais pas pour le reste ^^). Je résout un problème assez proche, sauf que je suis en cartésien, et la taille de mes mailles influence clairement la qualité des résultats, mais pour le pas de temps, c'est tout ou rien.

Et Franzis, pour ta condition de flux nul à la paroi, je pense que tu devrais ecrire le bilan de flux, et ensuite discrétiser l'équation que tu auras. Je n'ai plus en tête l'équation de fourier en coordonnées cylindrique, mais dans le cas générale, tu as :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

m.Cp.(dT/dt) = Flux_nord + Flux_sud + Flux_est + Flux_ouest

Suivant la paroi où tu es, un ou deux de ces quatre flux seront nuls. Tu mets ensuite l'expression du flux à partir de la formule de fourier, et tu discretise tout.

[Ehouarn]

Bonjour,

Kwatz: Je pense que ton expérience se borne à la recherche de solutions stationnaires. Et dans ces cas-là ce qui compte effectivement est d'aboutir à la solution à partir d'une condition initiale (typiquement complètement arbitraire). La bonne nouvelle, dans ces cas de figure, est qu'effectivement, si on converge, on converge vers la bonne solution; mais la mauvaise nouvelle est que le transitoire (l'évolution temporelle depuis la condition initiale vers l'état stationnaire) sera d'autant plus approximatif que le pas de temps est grossier.
Pour ce qui est de ta proposition pour la discrétisation, je ne suis pas sûr qu'utiliser une approche "volumes finis" sur la dernière maille, alors que le reste est formulé en différences finies soit très pratique.

Jean-Marc: Les méthodes implicites sont certes plus stables, mais la mise en œuvre plus complexe (et donc couteuses) vu les systèmes qu'il faut alors résoudre.
Sachant que

si le pas est trop petit, le temps de calcul devient prohibitif

est vrai quelque soit la méthode, c'est plutôt le cout global pour obtenir une solution satisfaisante donnée qu'il faut considérer. Et Kwatz a raison, si il y a un terme source qui a le mauvais gout d'être non-linéaire, une résolution implicite peut devenir cauchemardesque...

[kwatz]

Pour ce qui est de ta proposition pour la discrétisation, je ne suis pas sûr qu'utiliser une approche "volumes finis" sur la dernière maille, alors que le reste est formulé en différences finies soit très pratique.

Que tu discretises l'equation de la chaleur ou que tu raisonnes en bilan d'energie, tu obtiens la même expression, ce qui assez logique d'ailleurs. L'avantage du "volume finis" c'est que les relations sont facilement trouvables avec un peu de sens physique et jolie dessin

Kwatz: Je pense que ton expérience se borne à la recherche de solutions stationnaires. Et dans ces cas-là ce qui compte effectivement est d'aboutir à la solution à partir d'une condition initiale (typiquement complètement arbitraire).

Dans mon cas, je regarde l'évolution de température d'un morceau de plastique, mais je n'arrive jamais à un état stationnaire car dans mes conditions aux limites j'impose une température variante au cours du temps. Cependant, si que tu affirme est vrai, ce dont je doute, cela remet en cause mon approche du problème. D'après toi, comment trouve-t-on le "bon" pas de temps, c'est à dire celui en dessous duquel les résultats restent contants?

[Ehouarn]

Que tu discretises l'equation de la chaleur ou que tu raisonnes en bilan d'energie, tu obtiens la même expression, ce qui assez logique d'ailleurs

Pas exactement; suivant les méthodes utilisées les expressions discrètes diffèrent. Par contre elles convergent toutes asymptotiquement (lorque pas de temps et d'espace tendent vers zéro) vers la solution du problème continu. On est bien d'accord,
après ce n'est essentiellement qu'une question pratique d'implémentation d'une façon de faire ou d'une autre.

comment trouve-t-on le "bon" pas de temps, c'est à dire celui en dessous duquel les résultats restent contants?

Il n'y a pas de miracle, ce "bon" pas de temps ne peut être déterminé a priori (à moins de disposer de la solution analytique du problème, mais dans ce cas on n'en serait pas a faire une résolution numérique...). Il faut procéder par tâtonnement en réalisant des calculs avec des pas différents et évaluer la convergence. En toute rigueur le "taux de convergence" asymptotique vers la solution est dicté par l'ordre du schémas temporel (voir par ex. ici pour une petite description). En comparant l'écart entre des solutions numériques obtenues pour différents pas de temps et celle obtenue pour le plus petit des pas de temps employés, on peut déjà se faire une idée d'où en est la convergence (mais ce ne sera qu'une estimation et on ne peut hélas faire mieux).
Une approche complémentaire, un peu moins rigoureuse mais moins couteuse, est de faire ce genre d'étude sur un cas test aussi proche possible que celui que tu souhaites résoudre, mais pour lequel tu dispose d'une solution analytique. Tu peux dans ce cas très rigoureusement évaluer la convergence et le pas de temps "acceptable" pour ce cas test et ainsi en déduire qu'il faudra au moins ça pour gérer correctement le problème complet.