Discussion :

[Waldoon]

Bonjour,
Je résous numériquement l'équation de la chaleur en 1D:

(1/viscosité).dT/dt = (1/r).d/dr(r.dT/dr) en cylindrique avec hypothèse d'axisymétrie du problème
(1/viscosité).dT/dt = d²T/dr²) en cartésien sur x

Avec la viscosité constante, le problème est donc linéaire en T, en cartésien et en cylindrique.

On considère une température nulle initialement sur tout le domaine et l'objectif est de connaître
le profil de température à l'équilibre (solution stationnaire) en réponse à un échelon de température
en r=0.
Les conditions aux limites sont de type Dirichlet en r=0 (T=20) et r=L (T=0).
J'utilise les différences finies en discrétisant de manière centrée les opérateurs dérivés.
Le schéma numérique est implicite en temps.
Le système d'équation est linéaire et tridiagonal, j'ai donc utilisé l'algo TDMA pour résoudre le système.


J'obtiens le bon profil et unicité de la solution en coordonnées cartésienne quel que soit le type de maillage (uniforme et non-uniforme) et
le nombre de mailles.

Par contre, en coordonnées cylindrique, je n'obtiens pas de solution stationnaire unique, le profil change si le maillage
spatial est différent!!


Le problème étant linéaire en T je pense qu'il y a nul besoin de converger sur T à chaque pas de temps.
J'ai validé cette hypothèse en utilisant un solveur de système non-linéaire qui ne donne pas de solution unique après convergence.
J'ai vérifié moult fois la matrice, les conditions aux limites...

Je n'arrive pas à trouver mon erreur, si le problème est bien posé ou non.
Si quelqu'un a déjà résolu ce problème ou voit mon erreur, je suis tout ouï
Voici le code f90 du programme en question, j'ai essayé de le commenter correctement:

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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
PROGRAM calor 
  implicit none
  integer :: i,j,k
  integer, parameter :: Nr=1000 !Nombre de maille
  integer :: niter !Nombre d'itération en temps
  integer :: flg1,flg2,flg3,flg4,flg5,flg6,flg7,flg8,flg9,flg10
  real(kind=8), parameter :: L=1.0,PI=acos(-1.),visco=0.1
  real(kind=8), dimension (0:Nr) :: Lk,Dk,Uk,Sigma !Elément de matrice
  real(kind=8), dimension (0:Nr) :: SM,Tp
  real(kind=8), dimension (0:Nr) :: x
  real(kind=8), dimension (0:Nr) :: dx,dx12u,dx12b,x12u,x12b
  real(kind=8) :: t,dt,Time,test
 
  !Période simulée (s)
  Time = 10.0
  !Pas de temps (s)
  dt = 1E-3
  !Nb d'itérations nécessaire
  niter=Time/dt
  print*, 'dt=',dt
  print*, 'niter=',niter
 
  !Maillage
  do j=0,Nr
     x(j)=j*L/Nr !régulier
     !x(j) = L*(j**2)/(Nr**2) !irrégulier
     !x(j) = L*(j**4)/(Nr**4) !irrégulier
  enddo
 
  !Intervalles et coordonnées demi-mailles
  do j=1,Nr-1
     x12u(j)=0.5*(x(j)+x(j+1))
     x12b(j)=0.5*(x(j-1)+x(j))       
     dx12u(j)=x(j+1)-x(j)
     dx12b(j)=x(j)-x(j-1)
     !dx(j)=0.5*(dx12u(j)+dx12b(j))
     !dx(j)=x(j+1)-x(j) 
     dx(j)=x12u(j)-x12b(j)
  enddo
 
 
  !Initialisation  
  Tp(:) = 0.0
  SM(:) = Tp(:)
 
  !flags pour sauvegarde des profils de température à différents instants
  flg1=0;flg2=0;flg3=0;flg4=0;flg5=0;flg6=0;flg7=0;flg8=0;flg9=0;flg10=0
 
  t = 0.0
  !résolution du système sur j=1,Nr-1 (hors CL (points 0 et Nr))
  do while(t .LE. Time)
     t=t+dt
     !Coef de la matrice tridiagonale 
     Lk(:) = 0.0
     Dk(:) = 0.0
     Uk(:) = 0.0
 
     do j=1,Nr-1
        !Discrétisation de l'équation
 
        !Cylindrique 1D sur r (axie-symétrie)
        !Sigma(j) = visco*dt/(x(j)*dx(j))  
        !Lk(j) = -Sigma(j)*(x12b(j)/(dx12b(j)))
        !Dk(j) = 1+Sigma(j)*(x12u(j)/(dx12u(j)) + x12b(j)/(dx12b(j)))
        !Uk(j) = -Sigma(j)*(x12u(j)/(dx12u(j)))        
 
        !Cartésien en 1D sur x
        Sigma(j) = visco*dt/dx(j)
        Lk(j) = -Sigma(j)/dx12b(j)
        Dk(j) = 1+Sigma(j)*(1/dx12u(j) + 1/dx12b(j))
        Uk(j) = -Sigma(j)/dx12u(j)
 
        !Condition pour que l'algo TDMA soit pertinent        
        test = Dk(j)-(abs(Lk(j))+abs(Uk(j)))
        if (test .LT. 0) then
           print*, "Matrice non digonale-dominante!"
           print*, "j",j,"Lk=",Lk(j),"Dk",Dk(j),"Uk",Uk(j),test
        endif
     enddo
 
 
     !Conditions aux limites (ici Dirichlet)
     Tp(Nr) = 0.0
     Tp(0) = 20.0
     SM = Tp
     !Si CL de Dirichlet non-nulles, modif du second membre du système d'équation nécéssaire: 
     SM(1) = Tp(1) - Lk(1)*Tp(0)
     SM(Nr-1) = Tp(Nr-1) - Uk(Nr-1)*Tp(Nr)
 
 
     !Résolution du système
     call TDMA(Nr,Lk,Dk,Uk,SM,Tp)
 
     !Sauvegarde des profils
     if(int(t*100).EQ.1 .AND. flg1.NE.1) then
        flg1=1
        print*, "0.1s",t
        open(12,file="01s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
     if(int(t*100).EQ.2 .AND. flg2.NE.1) then
        flg2=1
        print*, "02s",t
        open(12,file="02s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
     if(int(t*100).EQ.3 .AND. flg3.NE.1) then
        flg3=1
        print*, "0.3s",t
        open(12,file="03s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
     if(int(t*100).EQ.4 .AND. flg4.NE.1) then
        flg4=1
        print*, "0.4s",t
        open(12,file="04s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
     if(int(t*100).EQ.5 .AND. flg5.NE.1) then
        flg5=1
        print*, "05s",t
        open(12,file="05s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
     if(int(t*100).EQ.6 .AND. flg6.NE.1) then
        flg6=1
        print*, "06s",t
        open(12,file="06s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
     if(int(t*100).EQ.1 .AND. flg7.NE.1) then
        flg7=1
        print*, "07s",t
        open(12,file="07s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
 
     if(int(t*100).EQ.8 .AND. flg8.NE.1) then
        flg8=1
        print*, "08s",t
        open(12,file="08s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
     if(int(t*100).EQ.500 .AND. flg9.NE.1) then
        flg9=1
        print*, "09s",t
        open(12,file="09s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
     if(int(t*100).EQ.1000 .AND. flg10.NE.1) then
        flg10=1
        print*, "1s",t
        open(12,file="1s.txt",status="unknown")
        do j=0,Nr
           write(12,*) x(j),Tp(j)
        enddo
        close(12)
     endif
 
  enddo
 
END PROGRAM calor
 
 
SUBROUTINE TDMA(N,A,B,C,D,X)
  integer, intent(in) :: N 
  real(kind=8), dimension(0:N) :: A,B,C,D
  real(kind=8), dimension(0:N) :: X
  real(kind=8) :: xmult
  integer :: i
 
  do i = 2,N-1
     xmult = A(i)/B(i-1)
     B(i) = B(i) - xmult*C(i-1)
     D(i) = D(i) - xmult*D(i-1)
  end do
  X(N-1) = D(N-1)/B(N-1)
  do i = N-2,1,-1
     X(i) = (D(i) - C(i)*X(i+1))/B(i)
  end do
END SUBROUTINE TDMA

[Ehouarn]

Bonjour,

Peux-tu préciser ce que tu entends par:

Le problème étant linéaire en T je pense qu'il y a nul besoin de converger sur T à chaque pas de temps.

Le schémas temporel est implicite soit; cela signifie simplement qu'il ne divergera pas, pas qu'il est précis.
Par contre on est bien d'accord que la solution stationnaire (asymptotique en temps) doit être atteinte.
Et si elle ne l'est pas, c'est forcément qu'il y a un problème dans l'algo. Coté code à proprement parler, je ne vois rien de choquant.
Coté méthode par contre je ne comprends pas bien comment tu construis ton système, notamment en ce qui concerne l'évaluation des dérivées spatiales; il faut qu'elles soient (au moins) d'ordre 2 et centrés; est-ce bien ce que tu tentes de faire? Je ne vois pas apparaître de terme en dérivée première qui apparaît pourtant dans l'équation du cas cylindrique.

[pvincent]

Je n'ai pas regardé le détail du code mais il me semble que la différence fondamentale entre coordonnées cartésiennes et cylindriques, c'est que si r tend vers zero, 1/r tend vers l'infini, donc il y a lieu d'examiner ce qui se passe dans l'équation (1) par un passage à la limite (qui peut éventuellement dépendre du maillage) avant d'imposer une condition en ce point.

Edit: la théorie des distributions de L. Schwartz est un moyen de le faire, mais ce n'est pas le seul.

Une petite remarque sur le code: en écrivant le nom du fichier à ouvrir dans une chaîne de caractères, on peut sérieusement limiter le nombre de lignes de code utilisées pour enregistrer les valeurs intermédiaires.

[Waldoon]

Je n'ai pas regardé le détail du code mais il me semble que la différence fondamentale entre coordonnées cartésiennes et cylindriques, c'est que si r tend vers zero, 1/r tend vers l'infini, donc il y a lieu d'examiner ce qui se passe dans l'équation (1) par un passage à la limite (qui peut éventuellement dépendre du maillage) avant d'imposer une condition en ce point.

Edit: la théorie des distributions de L. Schwartz est un moyen de le faire, mais ce n'est pas le seul.

En effet en imposant une limite non-nulle et assez supérieure à 0 en r1 (ex 1E-2), j'obtiens bien la solution analytique prévue et ce quel que soit le maillage et le pas d'espace. Lorsque r(1)=1E-6, le terme 1/r génère des erreurs sur le gradient de température en r1. Je vais chercher des références qui traitent de ces type de singularité en méthode numérique.

Une petite remarque sur le code: en écrivant le nom du fichier à ouvrir dans une chaîne de caractères, on peut sérieusement limiter le nombre de lignes de code utilisées pour enregistrer les valeurs intermédiaires.

Oui une boucle sur des strings c'est plus élégant.



Merci Ehouran et pvincent pour vos posts!

[Waldoon]

Ehouarn, non-linéaire en T, je voulais dire qu'il n'y a pas de terme dépendant de T (ex: viscosité(T) ou T²) devant chaque Tj dans l'équation discrétisée.
Les dérivées sont bien d'ordre deux et centrées dans mon schéma et je ne vois pas quel doit être le terme en dérivée première en cylindrique?

[Ehouarn]

Bonjour,

Es-tu bien certain de ta discrétisation? Je ne vois pas pourquoi dans les commentaires du code il est question de "demi-mailles" (je n'ai pas essayé de reprendre et suivre tous tes calculs non plus).
Pour ce qui est de la dérivée d'ordre 1, je parle du développement de l'équation qui le fait apparaître:
(1/viscosité).dT/dt = (1/r).d/dr(r.dT/dr) = (1/r).(dT/dr)+d2T/dr2
Il faut donc discrétiser la dérivée première et la dérivée seconde pour construire le système tridiagonal à résoudre.

Pour ce qui est de la singularité en r=0, ce n'est pas en problème avec des conditions aux limites Dirichlet, vu qu'on ne résout pas en ce point.
Et pour r petit, l'utilisation de la double précision, comme tu le fais devrait éviter les soucis, sauf circonstances exceptionnelles.