#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Mon Oct 24 09:06:32 2022 @author: adomps """ #!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Created on Mon Jul 18 10:48:03 2022 @author: adomps """ ###################################################################### # ## Résolution de l'éq de la chaleur 1d par méthode explicite ## conformément aux capacités numériques du programme officiel de 2022 # ####################################################################### # # # Le terme source, les conditions de bords et tous les paramètres # du problème sont traités comme des variables globales. # Les conditions de bord sont exprimés par diverses fonctions. # # Une fonction d'animation est importée depuis un fichier séparé. # # ############################################################################## import numpy as np import matplotlib.pyplot as plt import time from matplotlib import animation from math import erf #from fonction_animation import studio_poinca v_erf = np.vectorize(erf) def pdt_naif(T_anc) : ''' Fait évoluer l'intérieur du domaine. Les conditions de bord sont traitées à part. ''' T_nouv = np.zeros(N) for i in range(1, N-1) : T_nouv[i] = T_anc[i] + prefac * (T_anc[i+1] + T_anc[i-1] - 2 * T_anc[i]) \ + q * p_source[i] return T_nouv def pdt_vecto(T_anc) : T_nouv = np.zeros(N) T_nouv[1:-1] = T_anc[1:-1] + prefac * (T_anc[0:-2] + T_anc[2:] - 2 * T_anc[1:-1]) \ + q * p_source[1:-1] return T_nouv def construit_matrice(N, prefac) : M = np.eye(N) # avec des 1 dans les coins, les bords sont fixés for i in range(1, N-1) : M[i, i-1] = prefac M[i, i+1] = prefac M[i, i] -= 2 * prefac return M def pdt_matrice(T_anc) : T_nouv = np.dot(Mat, T_anc) T_nouv[1:-1] += q * p_source[1:-1] return T_nouv def bords_fixes(T_anc, T_nouv) : T_nouv[0] = T_anc[0] T_nouv[-1] = T_anc[-1] def impose_j_gauche(T_anc, T_nouv) : T_nouv[0] = T_anc[1] + j0 * delta_x / Lambda T_nouv[-1] = T_anc[-1] def evolution(T_init, N_pas, fonction_pdt, fonction_bords) : champ_T = np.zeros((N_pas+1, N)) champ_T[0] = T_init for i in range(N_pas) : champ_T[i+1] = fonction_pdt(champ_T[i]) fonction_bords(champ_T[i], champ_T[i+1] ) return champ_T def representation(x, champ_T, modulo, num_fig) : ''' Représente le champ de température en fonction de x, à diverses dates multiplies de modulo. ''' plt.figure(num_fig) N_pas = len(champ_T) for j in range(0, N_pas, modulo) : plt.plot(x, champ_T[j]) return # Champ de température initial : plusieurs variantes au choix def T_unif(x, T_interieur, T_gauche, T_droite) : ''' Température unfirme à l'intérieur, valeurs à choisir aux deux bords. ''' N = len(x) T = T_interieur * np.ones(N) T[0] = T_gauche T[-1] = T_droite return T def T_mode(x, t, n, A ): ''' Solution en mode propre. Pour l'utiliser comme condition initiale, prendre t=0. On peut bien sûr faire la somme de plusieurs modes.''' return A * np.sin(n * np.pi * x / L) * np.exp(- D * t * n** 2 * np.pi**2 / L**2) def source_uniforme(x, p_0) : ''' Définit une source uniforme de puissance volumique p_0. ''' N = len(x) return p_0 * np.ones(N) def source_portion(x, p_0, x_g, x_d): ''' Défini une source valant p_0 entre les abscisse x_g et x_d, et 0 ailleurs. ''' N = len(x) p = np.zeros(N) masque = (x_g <= x) * (x <= x_d) # multiplication de tableaux booléens, équivaut à AND p[masque] = p_0 return p def solution_1(x, t, T_paroi, T_avant) : ''' Solution exacte dans un milieu semi-infin initialement à T_avant dont la paroi x=0 est porté à T_paroi. ''' return T_paroi + (T_avant - T_paroi) * v_erf(x / (2 * np.sqrt(D * t))) ######################################################################### #### Paramètres numériques du problème #### traités comme des variables globales ######################################################################### L = 1. N = 51 Lambda = 400 # conductivité du cuivre D = 1.2e-4 # diffusivité thermique du cuivre x = np.linspace(0, L, N) delta_x = x[1] - x[0] delta_t = 0.1 prefac = 0.4 # en dimension 1, la condition de stabilité est prefac < 1/2 #prefac = D * delta_t / delta_x**2 print('prefac', prefac) delta_t = prefac * delta_x**2 / D q = D * delta_t / Lambda N_pas = 100000 Mat = construit_matrice(N, prefac) #p_source = source_portion(x, 100, 0.3, 0.8) ##################################################################### ### Applications ################################################################ ## Température imposée à une extrémité. p_source = source_uniforme(x, 0) T_init = T_unif(x, 0, 10, 0) t1 = time.time() champ_T = evolution(T_init, N_pas, pdt_naif, bords_fixes) print('durée', time.time() - t1) representation(x, champ_T, 100, 1) studio_poinca(x, champ_T, delta_t) # Température imposée à une extrémité et comparaison # avec la solution en erf. T_init = T_unif(x, 0, 10, 0) p_source = source_uniforme(x, 0) t1 = time.time() champ_T = evolution(T_init, N_pas, pdt_vecto, bords_fixes) print('durée', time.time() - t1) representation(x, champ_T, 100, 2) i_choisi = 200 # choisir un des indices de temps affichés t_choisi = i_choisi * delta_t # en secondes, prendre une date T_instant = solution_1(x, t_choisi, 10, 0) plt.plot(x, T_instant, linestyle = 'dashed', linewidth = 2, color = 'red') studio_poinca(x, champ_T, delta_t, intervalle = 100) # Température identique fixée aux deux bords, et condition # initiale en combinaison de modes propres. t1 = time.time() T_init = 10 + T_mode(x, 0, 1, 3) + T_mode(x, 0, 3, 3) champ_T = evolution(T_init, N_pas, pdt_matrice, bords_fixes) print('durée', time.time() - t1) representation(x, champ_T, 100, 3) i_choisi = 100 t_choisi = delta_t * i_choisi # en secondes T_instant = 10 + T_mode(x, t_choisi, 1, 3) + T_mode(x, t_choisi, 3, 3) plt.plot(x, T_instant, linestyle = 'dashed', linewidth = 2, color = 'red') studio_poinca(x, champ_T, delta_t) ################################################################### ## Exemple avec flux constant imposé à l'extrémité gauche S = 4.6e-4 # cf TP sur la barre chauffante Phi_0 = 20 j0 = Phi_0 / S T_init = T_unif(x, 0, 0, 0) champ_T = evolution(T_init, N_pas, pdt_vecto, impose_j_gauche) representation(x, champ_T, 100, 4) studio_poinca(x, champ_T, delta_t, intervalle = 1) ##################################################################### #### Exemple avec terme source ## avec source uniforme et bords fixes, on attent Delta T_max = p0 * L**3 / (8 * Lambda) T_init = T_unif(x, 0, 0, 0) p0 = 1000 print('Delta T max : ', p0 * L**2 / (8 * Lambda)) p_source = source_uniforme(x, p0) champ_T = evolution(T_init, N_pas, pdt_vecto, bords_fixes) representation(x, champ_T, 100, 5) studio_poinca(x, champ_T, delta_t)