Chauffer une maison restée froide : Euler explicite et temps de mise en régime

La stratégie. Le régime final se calcule par les résistances ; le régime transitoire, non : on discrétise l'équation de la chaleur et l'on traite les faces par la loi de Newton. On compare ensuite le temps de mise en régime à l'échelle $\frac{e^2}{D}$, et l'énergie stockée au bilan.

1) Le régime final. En série : $r=\frac{1}{h_i}+\frac{e}{\lambda }+\frac{1}{h_e}=0,130+0,114+0,040=0,284{\mathrm{m}}^2\mathrm{K}{\mathrm{W}}^{-1}$, d'où $$j_{\infty }=\frac{20-5}{0,284}=52,8\mathrm{W}{\mathrm{m}}^{-2},T_{si}=20-\frac{52,8}{7,7}=13,1{}^{\circ }\mathrm{C},T_{se}=5+\frac{52,8}{25}=7,1{}^{\circ }\mathrm{C}.$$

2) Le schéma. Les développements de Taylor $T(x\pm \Delta x)=T\pm \Delta x\frac{\partial T}{\partial x}+\frac{\Delta x^2}{2}\frac{{\partial }^{2}T}{\partial x^2}\pm \dots$ donnent, en les ajoutant, $\frac{{\partial }^{2}T}{\partial x^2}\approx \frac{T_{i+1}-2T_i+T_{i-1}}{\Delta x^2}$ ; Euler explicite approche $\frac{\partial T}{\partial t}$ par $\frac{T_i^{n+1}-T_i^n}{\Delta t}$. En reportant dans $\frac{\partial T}{\partial t}=D\frac{{\partial }^{2}T}{\partial x^2}$, on obtient le schéma annoncé : la nouvelle valeur se calcule explicitement à partir des anciennes.

3) Les bords. Le nœud $0$ reçoit $h_i(T_{ai}-T_0)$ de l'air et cède au nœud $1$, par Fourier, $\lambda \frac{T_0-T_1}{\Delta x}$ ; sans capacité, ces deux flux sont égaux : $$T_0=\frac{h_i{T}_{ai}+\frac{\lambda }{\Delta x}{T}_1}{h_i+\frac{\lambda }{\Delta x}},T_{N-1}=\frac{h_e{T}_{ae}+\frac{\lambda }{\Delta x}{T}_{N-2}}{h_e+\frac{\lambda }{\Delta x}}.$$ On calcule d'abord les nœuds intérieurs, puis les bords. Le prix de la simplification : les deux demi-mailles de bord ne stockent rien.

4) La stabilité. $T_i^{n+1}=r{T}_{i-1}^n+(1-2r)T_i^n+r{T}_{i+1}^n$ : les coefficients ont pour somme $1$ et sont positifs si $r⩽\frac{1}{2}$. La nouvelle valeur est alors comprise entre les anciennes : aucun extremum n'apparaît, comme dans la diffusion réelle. Si $r>\frac{1}{2}$, le coefficient central est négatif ; une perturbation $\epsilon (-1{)}^i$ devient $\epsilon (-1{)}^i(1-4r)$ à chaque pas, amplifiée si $|1-4r|>1$, c'est-à-dire $r>\frac{1}{2}$ : facteur $1,2$ par pas pour $r=0,55$. Les erreurs d'arrondi suffisent à la déclencher. Ici $D=7,6\times 10^{-7}{\mathrm{m}}^2{\mathrm{s}}^{-1}$ et $\Delta t=\frac{r\Delta x^2}{D}=13,1\mathrm{s}$ pour $r=0,4$.

5) Le script et son résultat.

import numpy as np
lam, rho, c, e = 1.75, 2300.0, 1000.0, 0.20  # beton ; e en m
hi, he, Tai, Tae, T0 = 7.7, 25.0, 20.0, 5.0, 5.0
D, N = lam/(rho*c), 41                        # 41 noeuds
dx = e/(N - 1)
g = lam/dx                        # conductance d'une maille

def simule(r, duree):
    dt = r*dx**2/D
    T = np.full(N, T0)                        # mur a 5 C
    t, E, hist = 0.0, 0.0, []
    while t < duree:
        T[1:-1] = T[1:-1] + r*(T[2:] - 2*T[1:-1] + T[:-2])
        T[0] = (hi*Tai + g*T[1])/(hi + g)     # Newton = Fourier
        T[-1] = (he*Tae + g*T[-2])/(he + g)
        jin, jout = hi*(Tai - T[0]), he*(T[-1] - Tae)
        t += dt
        E += (jin - jout)*dt                  # energie stockee
        hist.append((t, jin, jout))
    return T, np.array(hist), E, dt

T, H, E, dt = simule(0.4, 4*86400)
js = (Tai - Tae)/(1/hi + e/lam + 1/he)
print(f"dt = {dt:.1f} s ; j stationnaire = {js:.1f} W/m2")
for h in (1, 6, 15):
    k = np.searchsorted(H[:, 0], 3600*h)
    print(f"t = {h:2d} h : j_in = {H[k, 1]:5.1f}, "
          f"j_out = {H[k, 2]:5.1f} W/m2")
k = np.argmax(H[:, 2] >= 0.9*js)
print(f"j_out = 0,9 j a {H[k, 0]/3600:.1f} h ; "
      f"e^2/D = {e**2/D/3600:.1f} h")
print(f"energie stockee : {E/1e6:.2f} MJ/m2")
T, H, E, dt = simule(0.55, 3600)
print(f"r = 0,55 : max |T| apres 1 h = {abs(T).max():.0e} C")

Résultat obtenu en l'exécutant :

dt = 13.1 s ; j stationnaire = 52.8 W/m2
t =  1 h : j_in =  89.9, j_out =   0.4 W/m2
t =  6 h : j_in =  67.0, j_out =  27.2 W/m2
t = 15 h : j_in =  55.8, j_out =  47.4 W/m2
j_out = 0,9 j a 15.1 h ; e^2/D = 14.6 h
energie stockee : 2.30 MJ/m2
r = 0,55 : max |T| apres 1 h = 6e+11 C

La ligne T[1:-1] = T[1:-1] + r*(...) calcule tout le second membre avec les anciennes valeurs avant de l'affecter : c'est bien le schéma explicite.

Lecture : $j_{\text{out}}$ atteint $90\%$ de $j_{\infty }$ en $15,1\mathrm{h}$, et $\frac{e^2}{D}=14,6\mathrm{h}$ : la loi d'échelle donne le bon ordre de grandeur (le préfacteur proche de $1$ est ici une coïncidence de la définition « $90\%$ »). Énergie stockée exacte : le profil final est affine de $13,1$ à $7,1{}^{\circ }\mathrm{C}$, de moyenne $10,1{}^{\circ }\mathrm{C}$ ; le mur a gagné $\rho ce(10,1-5)=2,36\mathrm{MJ}{\mathrm{m}}^{-2}$. Le script en compte $2,30$, soit $2,5\%$ de moins : exactement la part des deux demi-mailles de bord sans capacité, $2\times \frac{\Delta x/2}{e}=\frac{1}{40}$. Avec $r=0,55$, le calcul explose (des températures de l'ordre de $10^{11}{}^{\circ }\mathrm{C}$ en une heure de temps simulé) : la condition $r⩽\frac{1}{2}$ est impérative.

6) Ce qu'on ressent, ce que ça coûte. Pendant la première heure, il faut fournir $90\mathrm{W}{\mathrm{m}}^{-2}$ ($1,7$ fois le régime final) et la surface intérieure n'est qu'à $T_{si}=20-\frac{89,9}{7,7}=8,3{}^{\circ }\mathrm{C}$ : l'air est à $20{}^{\circ }\mathrm{C}$, mais les murs froids rayonnent, et l'on a froid. Pour $100{\mathrm{m}}^2$ de murs, la mise en température stocke $236\mathrm{MJ}$, soit $66\mathrm{kWh}$, en plus des pertes : c'est le prix de l'inertie.

Ce que l'exercice installe. Le schéma d'Euler explicite : dérivée seconde centrée, $T_i^{n+1}=T_i^n+r(T_{i+1}^n-2T_i^n+T_{i-1}^n)$, stable si $r=\frac{D\Delta t}{\Delta x^2}⩽\frac{1}{2}$ (coefficients positifs), ce qu'on démontre et qu'on constate en exécutant. La loi de Newton comme condition aux limites : égalité du flux reçu de l'air et du flux de Fourier. Le temps de mise en régime suit $\frac{e^2}{D}$ ; un bilan d'énergie contrôle le calcul et mesure ses approximations.