Discussion :

[feynman]

C est ce que je fait

[jeyar]

d'acord, je vais vous passer un code écris par un professeur connu dans le domaine de la simulation numérique, c'est le professeur Toro, c'est un code en 2D (maillage rectangulaire sturcuté) qui résout les équations des écoulements d'eau peu profonde - shallow water - ( les équations de Saint venant ).
pour l'éxecution vous devez besoin d'un fichier de donnée je vais essayer de l'attaché avec cette réponse

bon courage

NB : ce code utilise la méthode des volumes finis

[haydiamet]

Bonjour tout le monde,
moi je traite pareil le même problème. 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 de le code en coordonnées cartésiennes) et le code simule, 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.
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ésous 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
 
!###########################################################################################################################################################################
!Résolution de l'équation de chaleur instationnaire 1-D en coordonnées sphériquess sans terme source avec lambda 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 31 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 cartésiennes (selon x) avec λ constante et avec la méthode de différence finie.
! 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: position qui indique les différentes positions là où on veut calculer la température.
! fichier temperature.data: affiche la température à chaque itération (dt) et dans tous les points 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
double precision :: dx  ! dx est le pas d'espace
double precision :: tps !tps est le temps
double precision :: tpsc !tps est le temps
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! x est le vecteur position, T est le vecteur temperature à l'instant t, Tn est le vecteur température à l'instant t+dt
 
 
 
!######################################################
!Déclaration des paramètres
!######################################################
 
integer :: i,j, it ! i est le 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=293K
Tg=2000.D0 !Tg=300K
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) !ρ=0,0013.10^3kg/m3 ; Cp=1.01KJ/Kg.K ; λ=0.026W/m.K ; D=20.10^-6 m²/s
dt=1.D-7
p=31415.D-4
dx=R/N !dx=1e-3
h=(D*dt)/(dx**2)
 
 
 
!##################################################
!Initialisation des vecteurs temperature (definition des positions) 
!##################################################
do i=0,N
x(i)=i*dx
xc(i)=x(i)
T(i)=Tp
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
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
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))!équation de la chaleur discrétisée sans adimensionnement
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 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.
Merci,