Discussion :

[haydiamet]

Bonjour tout le monde,
je traite le problème suivant:résolution de l'équation de chaleur en coordonnées sphériques (1D) en régime transitoire. j'ai résolu l'équation de chaleur en coordonnées cartésiennes (1D) avec la méthode de différences finies sans problème. Pour les coordonnées sphériques, j'ai codé le problème avec fortran (le même principe du code en coordonnées cartésiennes) et le code simule, mais le problème c'est que j'arrive pas à conserver le flux de chaleur à différentes positions en régime permanent (ce qui est faux physiquement)!!!
J'ai essayé de rester sur la méthode des différences finies en changeant plusieurs schéma de discrétisation mais toujours sans succès (équation adimensionnée et non adimensionnée à plusieurs adimensionnements).
Mon équation est la suivante: (1/D) dT/dt=d²T/dr²+(2/r)*dT/dr . les conditions aux limites sont simples: Températures imposée au centre de la sphère =Tg (constante=2000) et température imposée à la paroi (Tp=constante =450)
Je résoud l'équation en 1D suivant le rayon de la sphère R.
Voilà le code:

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
!###########################################################################################################################################################################
!Résolution de l'équation de chaleur instationnaire 1-D en coordonnées sphériquess sans terme source avec la conductivité thermique  constante: schéma explicite et conditions aux limites : températures imposées. prise en compte de la condition de courant de Friedriches-Lewy pour la convergence (dt<R²dx/D<0.045s)!! équation non adimensionnée!!!!!!!!!!!!!
!exemple avec Tg(centre)=300K et Tp(paroi)=293K.Nombre de positions x(i)=30 et tfinal=dt*itmax !!!!équation non adimensionnée
!exemple: application au refroidissement de la lune (Guillaume Marmin, 2009-2010) 
!###########################################################################################################################################################################
!on définit l'axe de propagation de la température de longueur R (3cm) et donc les 1001 positions x(i) et on fixe les températures aux limites qui restent constantes au cours du temps ainsi que les températures initiales au centre x(0). la résolution consiste à faire des itérations dans le temps pour pouvoir calculer les températures à toutes les positions x à chaque pas de temps dt. On résoud l'équation de chaleur (∂T/∂t)-(λ/ρCv)ΔT=0 en coordonnées sphériques (selon x) avec λ constante et avec la méthode de différences finies.
! on fixe les températures T et Tn égales, à l'initialisation (=Tp) puis au premier dt (itération) on impose des températures différentes aux limites x(0): Tg et x(R):Tp qui reste toujours constante.
!pour la convergence de la méthode: on fixe un critère de convergence h qui doit être <=0.5(voir méthodes numériques).
!La sortie du programme: fichier: données.data qui indique les  températures en fonction de x et du temps chaque 10000 itérations (sinon si c'est toutes les itérations on aura un fichier volumineux).
! des fichiers qui  affichent la valeur des flux conductifs à chaque itération (dt) à différentes positions x(i).
!###########################################################################################################################################################################
program calculeqchainst
 
implicit none
! Déclarations des variables et des paramètres
integer, parameter :: N = 1000 !N est le nombre de segments (N+1: nombre de noeuds selon x)
integer, parameter :: itmax=10000000 !itmax est le nombre maximal d'itérations dans le temps
double precision :: dx ! dx est le pas d'espace
double precision :: tps !tps est le temps (en secondes)
double precision :: tpsc !tpsc est le temps pour le traçage des courbes avec "gnuplot"
double precision :: dt !dt est le pas de temps
double precision :: rho !masse volumique de l'azote
double precision :: v !capacité calorififique en J/kg.K=741J/Kg.K
double precision :: dens !densité du gaz (n dans la formule de Carminati)
double precision :: lambda !conductivité thermique du gaz
double precision :: m !masse moléculaire du gaz
double precision :: R !R est le rayon de la chambre de combustion
double precision :: p
double precision :: Tp !Tp est la température de la paroi constante T(r=R,t): condition aux limites
double precision :: Tg !Tg est la température T(r=0,t): condition aux limites!!!! à vérifier car la température au centre de la chambre n'est pas toujours constante
double precision :: D !D=λ/ρC D est la diffusivité thermique qui est constante (car λ est constante)
double precision :: h !h=D*dt/dx**2 critère de convergence dans la mèthode de différence finie
double precision x(0:N), xc(0:N), T(0:N), Tn(0:N),Tc(0:N),Q2(1:itmax),Q500(1:itmax),Q200(1:itmax)
double precision Q800(1:itmax),Qcd(1:itmax), Qcdtot(1:itmax),Q999(1:itmax) ! x est le vecteur position, T est le vecteur temperature à l'instant t, 
!Tn est le vecteur température à l'instant t+dt! xc est le vecteur position pour le traçage des courbes sous gnuplot, Tc est le vecteur température à l'instant t+dt, (Q2,Q200,Q500, Q800, Qcdtot) sont les vecteurs flux aux noeuds, respectivement, (2,200, 500, 800, 1000=paroi), Qcd: densité de flux à la paroi (noeud 1000)
 
 
 
!######################################################
!Déclaration des paramètres
!######################################################
 
integer :: i,j, it ! i,j :compteur position. it est le compteur d'itération dans le temps
 
R=3.D-2 !rayon de la chambre=3cm
Tp=450.D0 !Tp=450K
Tg=2000.D0 !Tg=2000K
dens=25.D24 !n=densité=25.D24! densité de l'azote à pression athmosphérique pour un gaz à l'équilibre à 300K
m=46.D-27! masse d'une molécule d'azote
v=741.D0 !Cv: capacité calorifique à volume constant
rho=dens*m !rho=1.15Kg/m3
lambda=26.D-3 !=0.026W/m.K conductivité thermique de l'azote (égale à celle de l'air)
D=lambda/(rho*v) !diffusivité thermique constante
dt=1.D-7
p=31415.D-4
dx=R/N !dx=3e-5
h=(D*dt)/(dx**2)
 
 
 
!##################################################
!Initialisation des vecteurs temperature (definition des positions) 
!##################################################
do i=0,N
x(i)=i*dx !les positions 
xc(i)=x(i)
T(i)=Tp ! tout le rayon est maintenue à une température constante égale à celle cde la paroi=Tp=450
Tc(i)=T(i) ! à t0 la température dans tout le domaine est égale à Tp
Tn(i)=Tp ! Tn n'a pas d'importance à t0 car il n'y a pas d'itérations!! !Tn est le vecteur température à toutes les positions x(i) au temps suivant, c'est la solution de notre schéma de discrétisation
end do
 
 
!#######################################################
!ouverture des fichiers pour posttraitement
!######################################################
open (unit = 17,file = 'Donneescourbe.data')! fichier qui contient les valeurs de xc, Tc et le temps pour les itérations désirées ds le temps 
open (unit = 19,file = 'fluxconductif.data')! densité de flux conductif à la paroi à chaque pas de temps
write(19,*) 'Timestep=',dt,' flux conductif à la paroi au cours du temps'!Qcd pour chaque itération de temps
open (unit = 20,file = 'densitédefluxconductifcourbe.data', status='unknown')! densité de flux conductif à la paroi en fonction du temps pour gnuplot
open (unit = 21,file = 'fluxconductifcourbe.data', status='unknown')! flux conductif à la paroi en fonction du temps pour gnuplot
open (unit = 23,file = 'flux200courbe.data', status='unknown')! flux conductif àN=200 en fonction du temps pour gnuplot à N=200
open (unit = 24,file = 'flux800courbe.data', status='unknown')! flux conductif en fonction du temps pour gnuplot à N=800
open (unit = 25,file = 'flux999courbe.data', status='unknown')! flux conductif en fonction du temps pour gnuplot à N=999
open (unit = 26,file = 'flux2courbe.data', status='unknown')! flux conductif en fonction du temps pour gnuplot à N=2
open (unit = 27,file = 'flux500courbe.data', status='unknown')! flux conductif en fonction du temps pour gnuplot à N=500
 
!######################################################
 
!######################################################
!calcul avancé en temps
!######################################################
tps=0.
do it=1,itmax !début boucle principale (iteration dans le temps)
tps=tps+dt!le temps augmente
T(0)=Tg ! on force les conditions aux limites! le changement de température est effectué au centre de la sphère (càd à r=0)
T(N)= Tp !la température de la paroi reste constante
Tn(0)=Tg
Tn(N)=Tp
 
do i=1,N-1
!if (h.le.0.5) then ! critère de convergence du schéma explicite: h<=0.5
Tn(i)= ((h+(2*h/i))*T(i+1))+((1-2*h-(2*h/i))*T(i))+(h*T(i-1))![/B][/B]équation de la chaleur discrétisée sans adimensionnement
!schéma décentré pour dT/dr avec dT/dr=T(i+1,j)-T(i,j)/dx
end do
 
!calcul de flux conductif à travers la paroi à chaque pas de temps
 
Q2(it)=(lambda*(Tn(1)-Tn(2))/dx)*4*p*(x(2)**2) !flux au noeud 2
Q200(it)=(lambda*(Tn(199)-Tn(200))/dx)*4*p*(x(200)**2) !flux de chaleur au noeud 200
Q500(it)=(lambda*(Tn(499)-Tn(500))/dx)*4*p*(x(500)**2) !flux au noeud 500
Q800(it)=(lambda*(Tn(799)-Tn(800))/dx)*4*p*(x(800)**2) !flux au noeud 800
Q999(it)=(lambda*(Tn(998)-Tn(999))/dx)*4*p*(x(999)**2) !flux au noeud 999
Qcd(it)=lambda *(Tn(N-1)-Tp)/dx!densité de flux à la paroi à chaque pas de temps
Qcdtot(it)=Qcd(it)*4*p*(R**2) !flux total à la paroi rayon extérieur de la chambre à chaque pas de temps
 
do i=0,N
T(i)=Tn(i) !changer la température à chaque pas de temps
end do
 
!########################################################"
!sauvegarde des résultats à chaque pas de temps: à chaque pas de temps on affiche le vecteur position et le vecteur température
!######################################################"
 
!write(*,*)'iter=',it,';temps=',tps
!write(*,*)x
!write(*,*)T
 
if (MOD(it,10000)==0) then
do j=0,N
xc(j)=x(j)
Tc(j)=T(j)
tpsc=tps 
write(17,'(3(f0.24, 1x))')xc(j),Tc(j),tpsc !cette structure permet d'afficher xr, Tr et le temps et de les regrouper par 3, ainsi de les afficher en tant que
!réels de 24 chiffres après la virgule et sans espace depuis le début de la ligne "f0.24"; ces données sont séparés par un seul espace "1x" 
 
enddo
write(17,'(/)')!permet de créer une ligne blanche
endif
 
write(19,*)'iter=',it,'temps=',tps,dx,h,Qcd(it),Qcdtot(it)
 
end do ! fin boucle principale 
 
tps=0.
do it=1,itmax,10000
tps=tps+(10000*dt)!le temps augmente
write(20,'(2(f0.24, 1x))')tps,Qcd(it)
write(21,'(2(f0.24, 1x))')tps,Qcdtot(it)
write(23,'(2(f0.24, 1x))')tps,Q200(it)
write(24,'(2(f0.24, 1x))')tps,Q800(it)
write(25,'(2(f0.24, 1x))')tps,Q999(it)
write(26,'(2(f0.24, 1x))')tps,Q2(it)
write(27,'(2(f0.24, 1x))')tps,Q500(it)
enddo
close(17) 
close(19)
close(20)
close(21)
close(23)
close(24)
close(25)
close(26)
close(27)	
 
 
 
 
end program calculeqchainst !fin programme principal

Quelqu'un a une idée de solution pour remédier à ce problème (car si ça marche en coordonnées cartésiennes, pourquoi il bloque en coordonnées sphériques!!).
Merci d'avance,

[François L.]

Bonjour,

je ne sais pas si ça va changer beaucoup de choses, mais il ne faut pas utiliser un schéma décentré (qu'il soit amont ou aval) pour ce problème: la conduction a lieu dans les deux directions (ça c'est l'explication physique) et la matrice résultante est à diagonale dominante si on prend un schéma centré.

Les schémas décentrés sont à utiliser si on a des phénomènes convectifs (c'est assez rare dans un solide ), lorsque la convection domine la conduction (critère sur le Péclet de maille).

Et c'est une différence notable par rapport au cartésien, ce qui me fait dire que ça peut venir de là.

Bon, bien sûr, ça n'a rien a voir avec le fortran

[haydiamet]

Bonjour François,
En ce qui concerne le schéma centré, il marche pas dans mon cas (tu peux essayer de le vérifier à la main, juste en calculant la T(1,1): càd: Température au neoud 1 à la première itération): On s'aperçoit facilement qu'avec les CL et CI que j'ai introduit, la valeur de la température reste constante partout (càd égale à T(0,x)=Tp). le code a pu le faire.
Pour le traitement de la convection, dans mon cas je résoud l'éqaution de la chaleur à petit échelle dans un gaz et non pas dans un solide comme tu l'as compris, ce qui permet de considérer un transfert conductif dans un gaz.
Je pense que le problème vient des coordonnées sphériques ainsi que les CL. Vu qu'en coordonnées cartésiennes a bien donné des résultats et qu'en coordonnées sphériques le r=0 représente un point de singularité, même si je résoud pas mon système à r=0 dont je définis une température constante. La non conservation de flux conductif en régime permanent vient peut être de la présence de r² dans la formule de flux (Flux=-lambda*dT/dr*(4*pi*r²)), ainsi pour des r trop petits, on doit avoir un dT trop grand pour assurer un flux constant et vu que la température à r=0 est fixée, ce qui empêche le code à simuler bien le phénomène et donc à donner un grand écart de température aux premiers pas d'espace!

[Marlan]

Salut,

es-tu sûr de ton opérateur de dérivation ? Pourquoi ne pas avoir écrit le problème sous forme matrcielle ?

Tu devrais normalement avoir quelque chose comme :
Discrétisation temporelle :
T^{n+1}=T^n+dt*(d²T/dr²+(2/r)*dT/dr), où dt est le pas de temps
Discrétisation spatiale :
d²T_i/dr²=(T_{i+1}-2T_i+T_{i-1})/dr², avec dr constant le pas en espace,
dT_i/dr=(T_{i+1}-T_{i-1})/(2dr) pour la dérivée 1ère d'ordre 2
Conduisant à l'équation globale :
T_i^{n+1}=T_i^n+dt*[(T_{i+1}-2T_i+T_{i-1})/dr²+(2/r_i)*(T_{i+1}-T_{i-1})/(2dr)] où r_i=i*dr, pour i=1,N-1, où N est ton nombre de points : le problème étant symétrique, on a symétrie par rapport au centre de la sphère --> ce qui se passe pour r=]0;+infini[ se passe également pour =]0;-infini[.

Ensuite s'ajoute la condition aux limites (température constante au centre) soit :
Pour i=0 : T_0=constante=T_g
Pour i=N : T_N=constante=T_p

Est-ce bien ce que tu as codé ? Car la formule : Tn(i)= ((h+(2*h/i))*T(i+1))+((1-2*h-(2*h/i))*T(i))+(h*T(i-1)) ne me semble pas correspondre si ?

Bonne journée,

Marlan

[xflr6]

Bonjour,

Je suis aussi d'avis que les schémas centrés sont nécessaires dans ce type de cas.

Pour le problème susmentionné, je n'aurais pas osé :
- prendre des schémas différents pour les dérivées spatiales 1ère et 2nde alors que ces dernières sont les éléments du même opérateur Laplacien ;
- faire apparaître à l'instant initial (edit: initial ou pas d'ailleurs) une discontinuité au point singulier.

A mon avis, le schéma numérique n'est pas en cause (enfin si, dans le cas présent, prendre des schémas de dérivation différents dans l'expression d'un même opérateur uniquement parce que ça ne marche pas dans le cas le plus fail-safe, ce n'est pas une bonne idée du tout quand même), c'est plutôt le problème qui me semble mal posé.

Un thermostat sur deux sphéres concentriques mais avec la sphère intérieure à r>0 marcherait correctement.
Après, et ça, je ne sais pas si c'est une bonne idée, j'aurais en plus aussi fait commencer le réseau de pas dr à dr/2 pour éviter la singularité et s'en rapprocher.

Cordialement,
xflr6

[haydiamet]

Bonjour,
@Marlan: j'ai bien posé mon problème comme tu le décrivais, et j'ai essayé de prendre différents schémas. Le schéma que tu as développé je l'ai déjà essayé (schéma centré) et il ne marche pas (c'est ce que j'ai déjàexpliqué auparavant, même à la main, ça donne toujours une température constante et tu peux le vérifier à la main juste pour la première iteration)!

[haydiamet]

@ xflr6
Bonjour,

Je suis aussi du même avis que le fait de prendre des schémas d'ordre différents pour un même opérateur n'est pas une bonne alternative (mais c'était un moyen pour essayer de résoudre le problème).

- Pour ce que tu as signalé: "faire apparaître à l'instant initial (edit: initial ou pas d'ailleurs) une discontinuité au point singulier": Je sais qu'il vaut mieux éviter d'imposer des CL aux points singuliers (r=0 dans ce cas), mais c'est le phénomème physique qui impose ces conditions : je traite la réaction de combustion dans des chambres sphériques. C'est pour cette raison d'ailleurs que j'essaye de changer le type des CL et de prendre des conditions de type neumann (au centre). Ce problème alors se résoud par formulation matricielle et non pas par la méthode itérative. J'essaye aussi de commencer par un dr/2 (comme tu l'as mentionné) et voir ce que ça donne.
Bref, du boulot



Cordialement,
haydiamet

[ELHARAKI]

meme situation ,je prépare mon projet de fin d'études sur les trempes (refroidissement d'une bulle sphérique) alors que je dois faire la résolution numérique de l'équation de chaleur 1D en coordonnées sphérique et en régime transitoire (instationnaire) avec des conditions aux limites de convection
mon problème c'est que j'arrive pas à créer le code Fortan

SVP quelqu'un pour m'aider
Merci