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,