# CAPACITE 2 : EQUATION DE LA DIFFUSION THERMIQUE
#################################################################################################################################
import numpy as np
import matplotlib.pyplot as plt
#################################################################################################################
# Exemple 3 : Résolution de l'équation de la diffusion thermique avec flux surfacique imposé en x = O (condition 
# de Neumann) et loi de Newton en x = L (condition mixte)

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

def Equationdiffusion(M,N,Tini,jQ0):    

# 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)
    # jQ0 : flux surfacique en x = 0  (condition aux limites en x = 0) 
  
    T = np.zeros((M+1,N+1))  # Initialisation du tableau des températures T(ti,xk)
    T[0,:] = Tini            # Première ligne (indice i = 0) : température initiale des points xk    
    
    for i in range (0,M):    # Itération en temps ti = idt  (M valeurs de 0 à M-1)
    

        T[i+1,0] = T[i,0] + 2*C*(T[i,1]-T[i,0]+dx*jQ0/lambd)              # Condition aux limites en x = 0 
        T[i+1,N] = T[i,N] +2*C*(T[i,N-1]-T[i,N]-dx*h*(T[i,N]-TF)/lambd)   # Condition aux limites en x = L
        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
            T[i+1,k] = T[i,k] + C*(T[i,k+1] + T[i,k-1] - 2*T[i,k])     
   
    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

jQ0 = 50                        # Conditions aux limites en x = 0     (W.m-2)
TF = 10                         # Température du fluide extérieur     (°C)
h = 50                          # Coefficient d'échange de surface    (W.m-2.K-1)
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,jQ0)       # 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")
##################################################################################################################