#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Thu Nov 4 21:00:47 2021 @author: maximebouchereau """ import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation from mpl_toolkits.mplot3d import Axes3D from scipy.linalg import block_diag import sys # Résolution de l'équation de la chaleur en deux dimensions d'espace: dyke & encaissant # Paramètres physiques [à ajuster] L = 20 # Longueur du domaine H = 10 # Hauteur du domaine Hcouche2= 5 # Hauteur de la couche 2 Ldyke = 2 # Longueur du dyke Hdyke = 5 # Hauteur du dyke Tdyke = 1200 # Température du dyke kdyke = 3 # Conductivité thermique du dyke rhodyke = 2900 # Masse volumique du dyke cdyke = 1000 # Capacité thermique massique du dyke Tencaissant = 300 # Température de l'encaissant kencaissant = 1.5 # Conductivité thermique de l'encaissant (cas où les deux couches ont même conductivité) kcouche1 = 1.5 # Conductivité thermique de la couche 1 kcouche2 = 0.3 # Conductivité thermique de la couche 2 rhoencaissant = 2400 # Masse volumique de l'encaissant cencaissant = 1000 # Capacité thermique massique de l'encaissant tau = 1.577E7 # Durée de la simulation (en s - égal à 6 mois) # Paramètres de la simulation numérique [à ajuster] J = 50 # Nombre de subdivisions en espace selon chaque direction N = 200 # Nombre de subdivisions en temps # Paramètres calculés xx = np.linspace(0,L,J+1) # Discrétisation de l'intervalle des x yy = np.linspace(0,H,J+1) # Discrétisation de l'intervalle des y Ddyke = kdyke/(rhodyke*cdyke) # Coefficient de diffusion thermique du dyke Dencaissant = kencaissant/(rhoencaissant*cencaissant) # Coefficient de diffusion thermique de l'encaissant Dcouche1 = kcouche1/(rhoencaissant*cencaissant) # Coefficient de diffusion thermique de la couche 1 Dcouche2 = kcouche2/(rhoencaissant*cencaissant) # Coefficient de diffusion thermique de la couche 2 MatDiff = np.zeros((J+1,J+1)) # Donne le coefficient de diffusion thermique en fonction de l'endroit sur lequel on se trouve for i in range(J+1): for j in range(J+1): if abs(xx[i]-L/2) <= Ldyke/2 and yy[j] <= Hdyke: # Test de vérification si on se trouve sur le dyke MatDiff[i,j] = Ddyke else: # On n'est pas sur le dyke if yy[j] <= Hcouche2: # Test de vérification si on se trouve sur la couche 2 MatDiff[i,j] = Dcouche1 else: # C'est qu'on se trouve sur la couche 1 (hors du dyke) MatDiff[i,j] = Dcouche2 MatDiffVect = MatDiff.T.reshape((J+1)*(J+1)) # Transforme MatDiff en vecteur (changement de dimension) MatDiffNum = np.diag(MatDiffVect) # Transforme le vecteur précédent en matrice diagonale utilisable dans le schéma numérique hx = L/J # Pas de subdivision en espace selon x hy = H/J # Pas de subdivision en espace selon y ht = tau/N # Pas de subdivision en temps # Construction des matrices impliquées dans le schéma numérique # A est la matrice qui va imiter le laplacien en 1D selon x (avec conditions de Neumann) A = np.diag(np.ones(J),-1) + np.diag(np.ones(J),1) - np.diag(2*np.ones(J+1),0) A[0,0] = -1 A[J,J] = -1 # B est la matrice qui va imiter le laplacien en 2D (avec conditions de Neumann) F = (1/hx**2)*A-(2/hy**2)*np.eye(J+1,J+1) Fbis = (1/hx**2)*A-(1/hy**2)*np.eye(J+1,J+1) B = Fbis print(" ") print(" ") print("Construction de la matrice diagonale par blocs...") for j in range(J-1): sys.stdout.write("\r%d "%int(100*(j+1)/(J-1))+"%") sys.stdout.flush() #print("Assemblage de la matrice diagonale par blocs... {}% \r".format(str(int(100*(j+1)/(J-1))).rjust(3)), end="") B = block_diag(B,F) # Construction de la "diagonale" de B par blocs B = block_diag(B,Fbis) # On termine par le dernier bloc qui est différent B = B + np.diag((1/hy**2)*np.ones(J*(J+1)),J+1) + np.diag((1/hy**2)*np.ones(J*(J+1)),-J-1) # M est la matrice qui sera utilisée à chaque itération dans le schéma numérique I = np.eye((J+1)**2,(J+1)**2) # Matrice identité M = np.linalg.inv(I-ht*MatDiffNum@B) # Inversion à cause du schéma implicite # Calcul de la solution approchée de l'équation de la chaleur U = np.zeros(((J+1)**2,N+1)) # Collection de vecteurs contenant les solutions aux temps t_n # Construction de la condition initiale V = np.zeros((J+1,J+1)) print(" ") print(" ") print("Construction de la matrice associée à la condition initiale...") for i in range(J+1): for j in range(J+1): sys.stdout.write("\r%d "%int(100*(i*(J+1)+j+1)/((J+1)**2))+"%") sys.stdout.flush() if abs(xx[i]-L/2) <= Ldyke/2 and yy[j] <= Hdyke: # Test de vérification si on se trouve sur le dyke V[i,j] = Tdyke else: V[i,j] = Tencaissant U[:,0] = V.T.reshape((J+1)**2) # Calcul de la solution à chaque temps print(" ") print(" ") print("Calcul de la solution...") for n in range(N): sys.stdout.write("\r%d "%int(100*(n+1)/N)+"%") sys.stdout.flush() U[:,n+1] = M@U[:,n] # Calcul de la surface ayant dépassé les 600°C durant le processus nb600 = 0 # Compte le nombre de points où la température a dépassé les 600°C dans le domaine discrétisé for i in range((J+1)**2): if max(U[i,:]) >= 600: nb600 = nb600+1 print(" ") print(" ") print("Surface de l'auréole métamorphique:",format(nb600/((J+1)**2)*L*H - Ldyke*Hdyke,'.2f'),"m^2") # On retire Ldyke*Hdyke, car ce domaine n'est pas contitué d'encaissant # Tracé de la solution au bout de six mois plt.figure() plt.title("Température au bout de six mois [°C]") plt.xlabel("x") plt.ylabel("y") plt.imshow(np.flipud(U[:,N].reshape(J+1,J+1)),cmap="jet",aspect="equal",extent=(0,L,0,H)) plt.colorbar(format='%.f') plt.show() # Animation du tracé de la solution au cours du temps (bonus) deltat = 5000/N # Durée entre deux frames (en ms) fig=plt.figure() im = plt.imshow(np.flipud(U[:,0].reshape(J+1,J+1)),cmap="jet",aspect="equal",extent=(0,L,0,H)) plt.title("Evolution de la température [°C]") plt.xlabel("x") plt.ylabel("y") plt.colorbar(format='%.f') def animate(n): im.set_array(np.flipud(U[:,n].reshape(J+1,J+1))) return [im] anim = animation.FuncAnimation(fig, animate, frames=N, blit=True, interval=deltat, repeat=True) #%matplotlib qt plt.show()