#!/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) : ... return T_nouv def pdt_vecto(T_anc) : ... ... #recopier le document papier ... return T_nouv def construit_matrice(N, prefac) : M = np.eye(N) # matrice unité ... ... return M def pdt_matrice(T_anc) : ... ... return T_nouv def bords_fixes(T_anc, T_nouv) : ... def impose_j_gauche(T_anc, T_nouv) : ... def evolution(T_init, N_pas, fonction_pdt, fonction_bords) : champ_T = np.zeros((N_pas+1, N)) champ_T[0] = ... for i in range(N_pas) : ... ... 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 directes ################################################################ ## 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 = .... print('durée', time.time() - t1) representation(...) # Comparaison avec la solution exacte pur un milieu semi-infini 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(...) # 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 = .... print('durée', time.time() - t1) representation(...) 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(....) ################################################################### ## 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(...) representation(...) studio_poinca(...) ##################################################################### #### 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 = .... champ_T = evolution(....) representation(...) studio_poinca(...) ########################################################################## #### Réacteur chimique avec possible explosion, exo cf Landau # Selon la valeur du paramètre gamma, un régime stationnaire existe ou # non (cas d'explosion). Valeur critique gamma = 0.88 # Quand le régime stationnaire existe, on peut prédire T_max au centre. ## Se reporter à l'exo et au fichier << resolution_explosion_thermique >>. ########################################################################## def source_reaction_chi(T, p0, alpha, T0) : return p0 * np.exp(alpha * (T - T0)) def evolution_reacteur(T_init) : global p_source ... ... return champ_temp p0 = 40000 alpha = 0.01 ell = L/2 # ATTENTION : demi-longueur du réacteur Lambda = 400 T0 = 10 gamma = p0 * alpha * ell**2 / Lambda # valeur critique gamma = 0.88 T_init = T0 + T_mode(x, 0, 1, 0.6) + T_mode(x, 0, 3, -0.2) N_pas = 30000 champ_T = evolution_reacteur(T_init) from resolution_explosion_thermique_Landau import resolution_theta1 print('gamma', gamma) if gamma > 0.88 : print("Il n'existe pas de régime stationnaire.") else : theta1_attendu = resolution_theta1(gamma) # calculé avec le code << resolution_explosiont_thermique >> T_max_attendu = T0 + theta1_attendu / alpha print('Régime staionnaire avecT_max attendu ', T_max_attendu) representation(x, champ_T, 100, 5) studio_poinca(x, champ_T, delta_t) # studio_poinca(x, champ_T[0:500], delta_t) # Si gamma > 0.88, on prend seulement le début, après on a des NaN ########################################################################## ### Chauffage alternatif comme dans le TP sur la barre ########################################################################## def evolution_onde_flux(T_init) : global j0 ... ... ... return champ_T periode = 300 S = 4.6e-4 # section de la barre chauffante du TP Phi_0 = 20 # puissance en watts. jmax = Phi_0 / S L = .3 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 delta_t = prefac * delta_x**2 / D q = D * delta_t / Lambda print('prefac', prefac) N_pas = 20000 T_init = T_unif(x, 20, 20, 20) champ_T_onde = evolution_onde_flux(T_init) representation(x, champ_T_onde, 123, 6) t = [i * delta_t for i in range(N_pas + 1)] studio_poinca(x, champ_T_onde, delta_t, intervalle = 1) plt.figure(7) plt.clf() plt.plot(t, champ_T_onde[:, 1], color = 'black') plt.plot(t, champ_T_onde[:, 5], color = 'red') plt.plot(t, champ_T_onde[:, 9], color = 'blue') plt.xlim(0, 1500) plt.ylabel('T (°C)') plt.xlabel('t (s)')