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 = d2T/dr2) 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