Equation de la chaleur

[axuiop]

Bonjour,

Pour un projet je dois résoudre numériquement l'équation de la chaleur en coordonnées cylindriques avec un flux imposé comme conditions au limites. Je vous joins l'énoncé. Le premier exercice correspond à une résolution 1D du problème, j'ai résolu cet exercice par la méthode des différences finies. Je me suis inspiré de cet exercice pour résoudre l'equation en cylindrique mais j'ai du mal a écrire numériquement les conditions aux limites. Le programme que j'ai écrit fonctionne mais la courbe de températures obtenue n'est pas réaliste.
voici l'énoncé Projet.pdf
Voici le script que j'ai écrit :

import numpy as np
import scipy as sp
from math import exp
from math import log
from math import atan
from math import sin
from math import cos

import matplotlib as mpl
import matplotlib.pyplot as plt
a=0.000020
l=150
h1=10
T2=2000
T1=300
h2=500
r1=0.35
r2=0.25
R=0.60
Nz=24
L=0.120
tf=500
dt= 0.1
b=atan(0.35/0.55)
c=atan(0.35/0.65)
dz=L/Nz
dr1=b*dz
dr2=c*dz
Nr1=int(r1/dr1)
Nr2=int(r1/dr2)


Nt=int(tf/dt) # nombre de points dans le temps
     # pas de l'espace
Z=np.linspace(0,L,Nz)
u=np.linspace(0,tf,Nt)
R1=np.linspace(0.25,R,Nr1)

    #on va résoudre pour 0<z<0.55 avec une discrétisation adaptée aux frontières (dr1) puis pour 0,55<z<1.2 (dr2) 
    # on fait en sorte pour que un incrément de dz corresponde à un incrément de dr sur les frontières
    #on subidivise donc la tuyére en deux domaines
    
    #Condtions initiales:
T_1=np.zeros((Nr1,Nz,Nt))

for i in range(Nr1):
    for j in range(Nz):
        T_1[i,j,0]=T1
        L=T1
        

        
    

    #calcul de la température dans le premier domaine
w=0
for k in range(Nt-1):
    d=0
    p=0
    w=w+1
    if w==1000:
        print (w)
    for j in range(11,0,-1):
        d=d+1
        for i in range(d,Nr1-1):
            if d<=(Nr1-1):
                T_1[i,j,k+1]=T_1[i,j,k]+dt*a*((1/R1[i])*(T_1[i+1,j,k]-T_1[i,j,k])/(dr1)+(1/(dr1**2))*(T_1[i+1,j,k]-2*T_1[i,j,k]+T_1[i-1,j,k])+(1/(dz**2))*(T_1[i,j+1,k]-2*T_1[i,j,k]+T_1[i,j,k-1]))
            T_1[Nr1-1,j,k+1]=(h1*T1+l*(T_1[Nr1-2,j,k]/dr1))/((l/dr1)+h1)
        if d==1:
            T_1[d-1,j,k+1]=(h2*T2+l*T_1[d,j,k+1]*(-cos(a))/dr1+l*T2*sin(a)/dz)/(l*(-cos(a))/dr1+l*sin(a)/dz+l*h2)
        else:
            T_1[d-1,j,k+1]=(h2*T2+l*T_1[d,j,k+1]*cos(a)/dr1+l*T2*sin(a)/dz)/(l*cos(a)/dr1+l*sin(a)/dz+l*h2) #CAL sur le bord 
    for t in range(12,23):
        p=p+1
        for i in range(p+1,Nr1-1):
            if p<=(Nr1-1):
                T_1[i,t,k+1]=T_1[i,t,k]+dt*a*((1/R1[i])*(T_1[i+1,t,k]-T_1[i,t,k])/(dr1)+(1/(dr1**2))*(T_1[i+1,t,k]-2*T_1[i,t,k]+T_1[i-1,t,k])+(1/(dz**2))*(T_1[i,t+1,k]-2*T_1[i,t,k]+T_1[i,t,k-1]))

                
            #conditions aux limites
            T_1[Nr1-1,t,k+1]=(h1*T1+l*(T_1[Nr1-2,t,k]/dr1))/((l/dr1)+h1)
        if t==12:
            T_1[p,t,k+1]=(((cos(b)*T_1[p+1,t,k+1])/dr1)-((sin(b)*T_1[p,11,k+1])/dz))/((cos(b)/dr1)-(sin(b)/dz)) #raccordement avec T1
        
        T_1[p,t,k+1]=(((cos(b)*T_1[p+1,t,k+1])/dr1)-((sin(b)*T_1[p,t-1,k+1])/dz))/((cos(b)/dr1)-(sin(b)/dz)) #CAL sur le bord 
            

A=T_1[18,17,:]

plt.plot(u,A)

plt.show()

Pouvez vous me dire ce qui ne va pas dans mon script?
Merci
-----