# CAPACITE 2 : EQUATION DE LA DIFFUSION THERMIQUE
#################################################################################################################################
import numpy as np
import matplotlib.pyplot as plt
#################################################################################################################
# Exemple 4 : Résolution de l'équation de la diffusion thermique avec température imposée en x = 0 (condition de  Dirichlet) et flux surfacique nul en x = L (condition de Neumann) du à une paroi adiabatique

# Ecriture de la fonction de résolution de l'équation de la diffusion thermique

def Equationdiffusion(M,N,Tini,T0,jQL):    

# Grandeurs d'entrée :
    # M : nombre d'intervalles de la plage de temps [0,tmax] (donc M+1 valeurs du temps t)
    # N : nombre d'intervalles de la plage d'abscisse [0,L]  (donc N+1 valeurs de l'abscisses x)
    # Tini : température à t = 0 pour tout x  (condition initiale)
    # T0 : température en x = 0  (condition aux limites en x = 0) 
    # jQL : flux surfacique en x = L  (condition aux limites en x = L) 
           
  
    T = np.zeros((M+1,N+1))  # Initialisation du tableau des températures T(ti,xk)
    T[0,:] =                 # Première ligne (indice i = 0) : température initiale des points xk          A COMPLETER
    
    for i in range (0,M):    # Itération en temps ti = idt  (M valeurs de 0 à M-1)
    
        T[i+1,0] =                                           # Condition aux limites en x = 0              A COMPLETER
        T[i+1,N] =                                           # Condition aux limites en x = L              A COMPLETER
        
        for k in range(1,N):   # Pour chaque abscisse de x1 à xN-1 (N-1 valeurs de k de 1 à N-1)
            # Discrétisation équation de la diffusion
                                                                                                         # A COMPLETER     
   
    return T
   
# Paramètres physiques

lambd = 1.65                    # Conductivité thermique du mur        (W.m-1.K-1)
µ = 2150                        # Masse volumique du mur               (kg.m-3)  
c = 1000                        # capacité thermique massique du mur   (J.K-1.kg-1) 
D = lambd/(µ*c)                 # Diffusivité thermique du mur         (m2.s-1)
L = 0.40                        # Epaisseur du mur en m                (m)
tmax = 172800                   # Durée de l'étude                     (s)

# Paramètres numériques

N = 100                         # Nombre d'intervalles en abscisses 
dx = L/N                        # Pas en abscisse 
x = np.linspace(0,L,N+1)        # Plage d'abscisse

M = 17280                       # Nombre d'intervalles de temps 
dt = tmax/M                     # Pas de temps 

# Conditions initiales et aux limites

T0 =                            # Condition aux limites en x = 0      (°C)                                 A COMPLETER
jQL =                           # Conditions aux limites en x = L     (W.m-2)                              A COMPLETER  
Tini= 10                        # Température initiale                (°C)

# Constante intervenant dans l'équation de diffusion discrétisée

C = D*dt/dx**2                 

if C < 0.5:        # Test sur la condition de convergence de l'algorithme

    T = Equationdiffusion(M,N,Tini,T0,jQL)       # Résolution de l'équation de la diffusion 
                                                          
    plt.xlabel('x')
    plt.ylabel('T', rotation=0)
    plt.title('Equation de la diffusion thermique 1D - Exemple 1')

    # Tableau des itérations choisies pour le tracé de la température à différents instants
    it = np.array([0,int(M/100),int(M/10),int(2*M/10),int(4*M/10),int(6*M/10),int(8*M/10),M]) 
    t = it*dt                                 # Instants associés aux itérations choisies

    for j in range(0,len(it)):
        plt.plot(x,T[it[j],:],label='t={} s'.format(t[j]))

    plt.grid()
    plt.legend()
    plt.show()

else:
    print("Les pas sont mal choisis car la constante C est supérieure à 0,5")