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 |
Partager