import numpy as np
import matplotlib.pyplot as plt

L=2e-2 # plateau de la table de 2 cm d'épaisseur
N=40
dx=L/N

duree=3*3600

lamb=0.18 # bois conductivité thermique en W/m-1/K
D=0.15e-7    # coefficient de diffusivité

h=1       # coefficient conducto-convectif bois/air en W/m2/K

T1=30     # température imposée face supérieure
T0=20     # face inférieure

T=np.zeros(N)   # tableau stockant le profil de température
T2=np.zeros(N)  # et sa copie

# Profil à t=0 : conditions initiales
T[1:N]=T0
T[0]=T1

# convergence de la méthode si D*dt/(dx**2)<0.5
P=int(3*D*duree/dx**2) # on prend du coup un peu plus d'itérations (facteur 3)
dt=duree/P             # calcul du "pas" de temps

print("nombres d'iterations : "+str(P))
x=np.linspace(0,L*1000,N) # on créé x de 0 à L en mm

plt.figure(figsize=(10,8))
# méthode des éléments finis
for k in range(P):
    T2[0]=T[0]        # condition de Dirichlet en 0
    for i in range(1,N-1):
        T2[i]=T[i]+D*dt/(dx**2)*(T[i+1]+T[i-1]-2*T[i]) # équation de diffusion
                                                       # discrète : on
                                                       # travaille sur la copie
                                                       # T2

    T2[N-1]=20  #(h*T0+lamb/dx*T[N-2])/(h+lamb/dx)     # condition limite
                                                       # Dirichlet ou Newmann
    T=T2.copy()     # on recopie dans T

    if k%20==0:
        plt.clf()                 # efface le graphe
        plt.xlabel("en mm")
        plt.ylabel("en °C")
        plt.title("Profil de température à t = "+str(round(k*dt))+" s")
        plt.plot(x,T)
        plt.pause(0.01)           # affiche avec une pause de 0,01 s
plt.show()