# -*- coding: utf-8 -*- """ @author: Melzani M., www.mmelzani.fr Résolution numérique de l'équation de Laplace. Chaque exemple est mis sous la forme d'une fonction. Par exemple pour exécuter celui de la partie 2.1 de l'article, taper ex_2_1() dans le shell. """ ## -------------------------------------------------------------------------- import numpy as np import matplotlib.pyplot as plt plt.close('all') # ------------------------- Fonctions générales ----------------------------- def graphe_bords(B): """Affiche un graphique des bords du domaine. """ Nx,Ny = B.shape plt.figure() plt.imshow(B,origin='lower',cmap='binary', interpolation='Nearest') plt.colorbar() plt.title('matrice B des bords du domaine') plt.show() def graphe_equipot(B,f,a=0,nb_equipot=25): """Affiche les bords du domaine B en noir, et des surfaces équipotentielles de f. Prendre a=1 si on souhaite également afficher un colorplot des valeurs de f.""" Nx,Ny = B.shape plt.figure() x = np.linspace(0,Nx-1,Nx) # de 0 à Nx-1 inclus en 10 points y = np.linspace(0,Ny-1,Ny) plt.imshow(B,origin='lower',cmap='binary',interpolation='nearest') if (a==1): plt.imshow(f,origin='lower') # pour tracer un colorplot de f plt.colorbar() # idem cont=plt.contour(y,x,f,nb_equipot,colors='k') plt.title('equipotentielles') plt.clabel(cont,fmt='%1.1f') # le 1.1 controle le nombre de chiffres après la virgule sur l'affichage plt.show() def w_opt(Nx,Ny): return 2./(1.+np.pi/(Nx*Ny)*((Nx**2+Ny**2)/2.)**0.5) # -------- Suite : une fonction par exemple du texte --------------- #------------------------------------------------------------------- def ex_2_1(): """Exemple de la partie 2.1 : condensateur dans une boite.""" def initialisation_condensateur(B,V,a,b): """Initialise un condensateur plan, dont les plaques sont au potentiel a et b""" for j in range(Ny//2-Ny//10,Ny//2+Ny//10): B[Nx//2+Nx//20,j] = 1 V[Nx//2+Nx//20,j] = b B[Nx//2-Nx//20,j] = 1 V[Nx//2-Nx//20,j] = a def une_iteration_GS(B,f): f2 = f.copy() ecart = 0. for i in range(1,Nx-1): # on va de 1 à N-2 pour ne pas traiter le cadre (conditions de Dirichlet) for j in range(1,Ny-1): if B[i,j]==0: # on n'est alors pas sur un bord f[i,j] = (1.-w)*f[i,j] + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4. ecart = ecart + (f[i,j] - f2[i,j])**2 return (ecart/(Nx*Ny))**0.5 def iterations_GS(B,f,eps): ecart = eps+1. nb_iterations = 0 while ecart > eps: ecart = une_iteration_GS(B,f) nb_iterations += 1 return nb_iterations Nx, Ny = 100, 100 w = w_opt(Nx,Ny) B = np.zeros((Nx,Ny),bool) V = np.zeros((Nx,Ny)) initialisation_condensateur(B,V,-10.,10.) iterations_GS(B,V,1.e-5) graphe_bords(B) graphe_equipot(B,V,1) #------------------------------------------------------------------- def ex_2_2(): """Exemple de la partie 2.2 : étude de l'effet de pointe.""" def initialisation_effet_pointe(B,V,a,b,h,d): """Initialise un cadre avec en bas V=a, en haut V=b, à droite et à gauche V qui va linéairement de a jusqu'a b, et une tige de hauteur h et largeur d, au potentiel a, qui part du milieu du bord bas, et dont le bout est arrondi. Fonctionne correctement pour le bout si d est *impair*.""" # bord bas : B[0,:] = 1 V[0,:] = a # bord haut : B[Nx-1,:] = 1 V[Nx-1,:] = b # bord gauche : B[:,0] = 1 V[:,0] = np.linspace(a,b,Nx) # bord droit : B[:,Ny-1] = 1 V[:,Ny-1] = np.linspace(a,b,Nx) # Pointe : B[0:h,Ny//2-d//2:Ny//2+(d+1)//2] = 1 V[0:h,Ny//2-d//2:Ny//2+(d+1)//2] = a # Bout arrondi pour la pointe : i0 = h j0 = Ny//2 for i in range(Nx): for j in range(Ny): if (i-i0)**2+(j-j0)**2 <= (d//2)**2: # cercle de centre i0, j0 et de rayon d//2 B[i,j] = 1 V[i,j] = a def une_iteration_GS(B,f): f2 = f.copy() ecart = 0. for i in range(1,Nx-1): # on va de 1 à N-2 pour ne pas traiter le cadre (conditions de Dirichlet) for j in range(1,Ny-1): if B[i,j]==0: # on n'est alors pas sur un bord f[i,j] = (1.-w)*f[i,j] + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4. ecart = ecart + (f[i,j] - f2[i,j])**2 return (ecart/(Nx*Ny))**0.5 def iterations_GS(B,f,eps): ecart = eps+1. nb_iterations = 0 while ecart > eps: ecart = une_iteration_GS(B,f) nb_iterations += 1 return nb_iterations def calcul_E(V): Ex = np.zeros((Nx,Ny)) Ey = np.zeros((Nx,Ny)) norme_E = np.zeros((Nx,Ny)) for i in range(1,Nx-1): for j in range(1,Ny-1): Ex[i,j] = -(V[i+1,j]-V[i-1,j])/(2.*delta) Ey[i,j] = -(V[i,j+1]-V[i,j-1])/(2.*delta) norme_E = (Ex**2+Ey**2)**0.5 return Ex, Ey, norme_E Nx, Ny = 150, 150 w = w_opt(Nx,Ny) delta = 3.e-2 # pas de la grille en m B = np.zeros((Nx,Ny),bool) V = np.zeros((Nx,Ny)) initialisation_effet_pointe(B,V,0.,(Nx-1)*delta*100.,50,13) iterations_GS(B,V,1.e-3) Ex, Ey, norme_E = calcul_E(V) max_E = np.max(norme_E) print("Valeur maximale de E : " + str(max_E) + "V/m") graphe_equipot(B,V,0) #------------------------------------------------------------------- def ex_2_3(): """Exemple de la partie 2.3 : solide avec T fixée sur les bords.""" def initialisation_cadre(B,T,T0,T1): """Exemple de la section 4.1. Fixe des températures sur le bord du domaine : T0 en bas, à gauche et à droite, et T1 en haut.""" B[:,0] = 1 B[:,Ny-1] = 1 T[:,0] = T0 T[:,Ny-1] = T0 B[0,:] = 1 B[Nx-1,:] = 1 T[0,:] = T0 T[Nx-1,:] = T1 def une_iteration_GS(B,f): ecart_max = 0. # actualisation sur tout le domaine, i allant de haut en bas, j de gauche à droite for i in range(Nx): for j in range(Ny): if B[i,j]==0: temp = f[i,j] f[i,j] = (1.-w)*temp + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4. ecart_max = max(ecart_max, abs(temp-f[i,j])) return ecart_max def iterations_GS(B,f,eps): ecart = eps+1. nb_iterations = 0 while ecart > eps: ecart = une_iteration_GS(B,f) nb_iterations += 1 return nb_iterations Nx, Ny = 20, 60 w = w_opt(Nx,Ny) T0 = 60. # degrés, temperature bords droit, gauche, inferieur T1 = 20. # degrés, temperature bord superieur Tinitiale = (T1+T0)/2. # temperature initiale B = np.zeros((Nx,Ny),bool) T = np.zeros((Nx,Ny)) + Tinitiale initialisation_cadre(B,T,T0,T1) iterations_GS(B,T,1.e-5) # on peut mettre eps=1.e-10. graphe_equipot(B,T,0) #------------------------------------------------------------------- def ex_2_4_i(): """Exemple de la partie 2.4(i) : barre calorifugée sur les côtés, T fixé à gauche, T fixé à droite.""" def initialisation_cadre(B,T,T0,T1): """Exemples de la section 4.2""" B[:,0] = 1 B[:,Ny-1] = 1 T[:,0] = T0 T[:,Ny-1] = T1 B[0,:] = 1 B[Nx-1,:] = 1 T[0,:] = T0 T[Nx-1,:] = T0 def une_iteration_GS(B,f): ecart_max = 0. for i in range(Nx): for j in range(Ny): temp = f[i,j] if B[i,j]==0: # on n'est alors pas sur un bord f[i,j] = (1.-w)*temp + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4. if i==0 and j>0 and j0 and j eps: ecart = une_iteration_GS(B,f) nb_iterations += 1 return nb_iterations Nx, Ny = 10, 100 w = w_opt(Nx,Ny) T0 = 100. # degrés, temperature de la face chaude T1 = 20. # degrés, temperature de la face froide Tinitiale = (T0+T1)/2. # temperature initiale B = np.zeros((Nx,Ny),bool) T = np.zeros((Nx,Ny)) + Tinitiale initialisation_cadre(B,T,T0,T1) iterations_GS(B,T,1.e-5) graphe_equipot(B,T,0) plt.figure() x = np.linspace(0,Nx-1,Nx) y = np.linspace(0,Ny-1,Ny) plt.subplot(211) plt.plot(y,T[Nx//2,:],'o', label='T fin de simulation') plt.plot(y,(T1-T0)/(Ny-1)*y+T0, label='T analytique') plt.xlabel('j') plt.ylabel(u'T (degrés)') plt.legend(loc='best') plt.grid() plt.subplot(212) plt.plot(y,T[Nx//2,:]-((T1-T0)/(Ny-1)*y+T0),'o', label='T[Nx//2,:] - Tanalytique') plt.xlabel('j') plt.ylabel(u'Delta T (degrés)') plt.legend(loc='best') plt.grid() plt.show() plt.figure() plt.plot(x,(T[:,Ny//2]-max(T[:,Ny//2])),'o', label='T - Tmax dans la tranche y=Ny//2') plt.xlabel('i') plt.ylabel('T - Tmax dans la tranche y=Ny//2') plt.title('T - Tmax dans la tranche y=Ny//2') plt.legend(loc='best') plt.grid() plt.show() #------------------------------------------------------------------- def ex_2_4_ii(): """Exemple de la partie 2.4(ii) : barre calorifugée sur les côtés, T fixé à gauche, flux fixé à droite.""" def initialisation_cadre(B,T,T0): B[:,0] = 1 B[:,Ny-1] = 1 T[:,0] = T0 T[:,Ny-1] = T0 B[0,:] = 1 B[Nx-1,:] = 1 T[0,:] = T0 T[Nx-1,:] = T0 def une_iteration_GS(B,f): ecart_max = 0. # actualisation sur tout le domaine, i allant de haut en bas : for i in range(Nx): for j in range(Ny): temp = f[i,j] if B[i,j]==0: # on n'est alors pas sur un bord f[i,j] = (1.-w)*temp + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4. if j==Ny-1 and i>0 and i0 and j0 and j eps: ecart = une_iteration_GS(B,f) nb_iterations += 1 return nb_iterations Nx, Ny = 10, 100 w = w_opt(Nx,Ny) T0 = 100. # degrés, temperature de la face chaude Tinitiale = T0 # temperature initiale de l'ailette delta = 1.e-2 # pas de la grille en m lambd = 400. # W.K^-1.m^-1, conductivité thermique jth = 15.*80. # flux thermique sur le bord droit de l'ailette jth_ad = delta*jth/lambd # flux thermique adimensionné B = np.zeros((Nx,Ny),bool) T = np.zeros((Nx,Ny)) + Tinitiale initialisation_cadre(B,T,T0) iterations_GS(B,T,1.e-5) graphe_equipot(B,T,0) plt.figure() x = np.linspace(0,Nx-1,Nx) y = np.linspace(0,Ny-1,Ny) pente = jth_ad # pente théorique plt.subplot(211) plt.plot(y,T[Nx//2,:],'o', label='T fin de simulation') plt.plot(y,-pente*y+T0, label='T analytique') plt.xlabel('j') plt.ylabel(u'T (degrés)') plt.legend(loc='best') plt.grid() plt.subplot(212) plt.plot(y,T[Nx//2,:]+pente*y-T0,'o', label='T[Nx//2,:] - Tanalytique') plt.xlabel('j') plt.ylabel(u'Delta T (degrés)') plt.legend(loc='best') plt.grid() plt.show() plt.figure() plt.plot(x,(T[:,Ny//2]-max(T[:,Ny//2])),'o', label='T - Tmax dans la tranche y=Ny//2') plt.xlabel('i') plt.ylabel('T - Tmax dans la tranche y=Ny//2') plt.title('T - Tmax dans la tranche y=Ny//2') plt.legend(loc='best') plt.grid() plt.show() #------------------------------------------------------------------- def ex_2_4_iii(): """Exemple de la partie 2.4(iii) : barre calorifugée sur les côtés, T fixé à gauche, flux conducto-convectif à droite.""" def initialisation_cadre(B,T,T0): """Initialise le domaine et ses bords (matrice B).""" B[:,0] = 1 B[:,Ny-1] = 1 T[:,0] = T0 T[:,Ny-1] = T0 B[0,:] = 1 B[Nx-1,:] = 1 T[0,:] = T0 T[Nx-1,:] = T0 def une_iteration_GS(B,f): """Réalise une itération via la méthode de Gauss-Seidel. B matrice des bords, f matrice de la fonction. Retourne l'écart (e dans l'article).""" ecart_max = 0. # actualisation sur tout le domaine, i allant de haut en bas, j de gauche à droite : for i in range(Nx): for j in range(Ny): temp = f[i,j] if B[i,j]==0: # on n'est alors pas sur un bord f[i,j] = (1.-w)*temp + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4. if j==Ny-1 and 0 flux fixé par conductoconvectif if i==0 and 0 flux nul if i==Nx-1 and 0 flux nul ecart_max = max(ecart_max, abs(temp-f[i,j])) f[Nx-1,Ny-1] = 0.5*(f[Nx-2,Ny-1]+f[Nx-1,Ny-2]) # traiter les coins (inutile mais f[0,Ny-1] = 0.5*(f[1,Ny-1]+f[0,Ny-2]) # plus joli sur le rendu graphique) return ecart_max def iterations_GS(B,f,eps): """Enchaine les itérations jusqu'à ce que le critère de convergence eps soit atteint. B matrice des bords, f matrice de la fonction. eps = epsilon dans l'article.""" ecart = eps+1. nb_iterations = 0 while ecart > eps: ecart = une_iteration_GS(B,f) nb_iterations += 1 return nb_iterations Nx, Ny = 10, 100 # taille du domaine w = w_opt(Nx,Ny) # calcul du paramètre de sur-relaxation omega optimal Text = 10. # degrés, temperature exterieure T0 = 100. # degrés, temperature de la face chaude Tinitiale = T0 # degrés, temperature initiale de l'ailette delta = 1.e-2 # pas de la grille en m lambd = 400. # W.K^-1.m^-1, conductivité thermique h = 15. # W.m^-2.K^-1, coefficient conducto-convectif de Newton a = delta*h/lambd # paramètre intervenant dans les CL conducto-convectives B = np.zeros((Nx,Ny),bool) # matrice des bords (0 ou 1) T = np.zeros((Nx,Ny)) + Tinitiale # matrice contenant la température initialisation_cadre(B,T,T0) # on initialise le domaine iterations_GS(B,T,1.e-7) # on effectue les itérations graphe_equipot(B,T,0) # tracé des équipotentielles plt.figure() x = np.linspace(0,Nx-1,Nx) y = np.linspace(0,Ny-1,Ny) pente = (T0-Text)/(Ny-1+1./a) # pente théorique pour comparaison plt.subplot(211) plt.plot(y,T[Nx//2,:],'o', label=u'T fin de simulation') plt.plot(y,-pente*y+T0, label=u'T analytique') plt.xlabel('j') plt.ylabel(u'T (degrés)') plt.legend(loc='best') plt.grid() plt.subplot(212) plt.plot(y,T[Nx//2,:]+pente*y-T0,'o', label=u'T[Nx//2,:] - Tanalytique') plt.xlabel('j') plt.ylabel(u'Delta T (degrés)') plt.legend(loc='best') plt.grid() plt.show() plt.figure() plt.plot(x,(T[:,Ny//2]-max(T[:,Ny//2])),'o', label='T - Tmax dans la tranche y=Ny//2') plt.xlabel('i') plt.ylabel('T - Tmax dans la tranche y=Ny//2') plt.title('T - Tmax dans la tranche y=Ny//2') plt.legend(loc='best') plt.grid() plt.show() #------------------------------------------------------------------- def ex_2_5(): """Exemple de la partie 2.5 : ailette de refroidissement.""" def initialisation_cadre(B,T,T0): B[:,0] = 1 B[:,Ny-1] = 1 T[:,0] = T0 T[:,Ny-1] = T0 B[0,:] = 1 B[Nx-1,:] = 1 T[0,:] = T0 T[Nx-1,:] = T0 def une_iteration_GS(B,f): ecart_max = 0. # actualisation sur tout le domaine, i allant de haut en bas for i in range(Nx): for j in range(Ny): temp = f[i,j] if B[i,j]==0: # on n'est alors pas sur un bord f[i,j] = (1.-w)*temp + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4. if j==Ny-1 and i>0 and i0 and j0 and j eps: ecart = une_iteration_GS(B,f) nb_iterations += 1 return nb_iterations Nx, Ny = 32, 100 w = w_opt(Nx,Ny) Text = 10. # degrés, temerature exterieure T0 = 100. # degrés, temperature de la face chaude Tinitiale = T0 # degrés, temperature initiale de l'ailette delta = 1.e-3 # pas de la grille en m lambd = 400. # W.K^-1.m^-1, conductivité thermique h = 15. # W.m^-2.K^-1, coefficient conducto-convectif de Newton a = delta*h/lambd # paramètre intervenant dans les CL conducto-convectives delta_p_ad = ((Nx-1)/(2.*a))**0.5 # facteur exponentiel dans le cas ailette alpha = 1./(a*delta_p_ad) beta = (Ny-1)/delta_p_ad A1 = 1./(1.+(alpha-1.)/(alpha+1.)*np.exp(-2.*beta)) B1 = 1./(1.+(alpha+1.)/(alpha-1.)*np.exp(+2.*beta)) print(" Parametres de la solution analityque 1D : ") print(" alpha = " + str(alpha)) print(" beta = " + str(beta)) print(" delta_p_ad = " + str(delta_p_ad)) print(" A1 = " + str(A1) + "et B1 = " + str(B1)) B = np.zeros((Nx,Ny),bool) T = np.zeros((Nx,Ny)) + Tinitiale initialisation_cadre(B,T,T0) iterations_GS(B,T,1.e-7) graphe_equipot(B,T,0) plt.figure() x = np.linspace(0,Nx-1,Nx) y = np.linspace(0,Ny-1,Ny) plt.subplot(211) plt.plot(y,(T[Nx//2,:]-Text)/(T0-Text),'o',label=u'fin de simulation') plt.plot(y,np.exp(-y/delta_p_ad), label=u'analytique si delta_p=infty') plt.plot(y,A1*np.exp(-y/delta_p_ad)+B1*np.exp(+y/delta_p_ad), label=u'analytique si delta_p fini') plt.xlabel('j') plt.ylabel(u'(T - T0)/(T0-Text)') plt.legend(loc='best') plt.grid() plt.subplot(212) plt.plot(y,(T[Nx//2,:]-Text)/(T0-Text)-A1*np.exp(-y/delta_p_ad)-B1*np.exp(+y/delta_p_ad), label=u'(T[Nx//2,:]-Text)/(T0-Text) - (Tanalytique-Text)/(T0-Text)') plt.xlabel('j') plt.ylabel(u'"Delta T "') plt.legend(loc='best') plt.grid() plt.show() plt.figure() plt.plot(x,(T[:,Ny//2]-max(T[:,Ny//2])),'o', label='T - Tmax dans la tranche y=Ny//2') plt.xlabel('i') plt.ylabel('T - Tmax dans la tranche y=Ny//2') plt.title('T - Tmax dans la tranche y=Ny//2') plt.legend(loc='best') plt.grid() plt.show() #------------------------------------------------------------------- def ex_2_6(): """Exemple de la partie 2.6 : pont thermique.""" def initialisation_bords(B,T): # mur tout à gauche : B[:,0] = 1 T[:,0] = Tdehors # bout du mur tout en haut et tout en bas : B[0,0:e] = 1 B[Nx-1,0:e] = 1 T[0,0:e] = Tdehors T[Nx-1,0:e] = Tdehors # bord droit du mur vertical, les CL seront conductoconvectives : B[0:Nx//2-h//2, e-1] = 1 B[Nx//2+h//2+1:Nx, e-1] = 1 # bord verticalde la dalle tout à droite, les CL seront conductoconvectives : B[Nx//2-h//2:Nx//2+h//2, Ny-1] = 1 # limites haute et basse de la dalle, les CL seront conductoconvectives : B[Nx//2-h//2,e:] = 1 B[Nx//2+h//2,e:] = 1 # pièce intérieure : on ne fera pas de calcul dans cette zone B[:Nx//2-h//2, e:] = 2 # 2 = on ne fait pas de calcul dans cette zone. B[Nx//2+h//2+1:, e:] = 2 # 2 = on ne fait pas de calcul dans cette zone. T[:Nx//2-h//2, e:] = Tmaison # on initialise tout de même une température, T[Nx//2+h//2+1:, e:] = Tmaison # même si elle ne sert pas. def une_iteration_GS(B,f): ecart_max = 0. # actualisation sur tout le domaine, i allant de haut en bas for i in range(0,Nx): for j in range(0,Ny): temp = f[i,j] if B[i,j]==0: # on n'est alors pas sur un bord f[i,j] = (1.-w)*temp + w*(f[i+1,j]+f[i,j+1]+f[i-1,j]+f[i,j-1])/4. elif B[i,j]==1: # on est sur un bord # Les trois 1er if et elif sont la limite entre le mur et le coté extérieur (à 0°C par ex.). # Les commenter si on préfère imposer T fixe (à 0°C par ex.). if j==0 and 0 eps: ecart = une_iteration_GS(B,f) nb_iterations += 1 return nb_iterations Nx,Ny = 300,180 w = w_opt(Nx,Ny) Tmaison = 20. # degrés, temperature de l'intérieur de la maison Tdehors = 0. # degrés, température à l'extérieur de la maison Tinitiale = 0. # degrés, temperature initiale pour l'algorithme delta = 1.e-2 # pas de la grille en m h = 50 # épaisseur de la dalle, en nb de cellules e = 50 # épaisseur du mur, en nb de cellules lambd = 2. # W.K^-1.m^-1, conductivité thermique hNewton = 15. # W.m^-2.K^-1, coefficient conducto-convectif de Newton a = delta*hNewton/lambd # paramètre intervenant dans les CL conducto-convectives B = np.zeros((Nx,Ny))#,bool) T = np.zeros((Nx,Ny)) + Tinitiale initialisation_bords(B,T) iterations_GS(B,T,1.e-4) T[B==2] = np.nan T[B==1] = np.nan B[B==2] = np.nan graphe_equipot(B,T,1,15)