Simulation de la température dans le sol — Physique des ondes 2 : Phénomènes de propagation linéaires : absorption et dispersion

🖥️ Simulation de la température dans le sol

Notebook Capytale de cet exercice : b601-11678469

Dans le sol, la température vérifie une équation de diffusion de coefficient $𝐷 \approx 1$ ⋅10−6 m2 s−1. Les données à la surface sont acquises régulièrement par des stations météorologiques. Les données pour Quimper sont disponibles à l’adresse suivante : https://nuage03.apps.education.fr/index.php/s/LgXjiwkxJxcrZmz. Le temps est donné en secondes depuis le 1er janvier 1970 à 00:00:00 UTC. La température est donnée en degrés Celsius. Les mesures sont effectuées toutes les heures.

On peut importer les données dans Python grâce aux instructions suivantes

import numpy as np
data = np.load("températures Quimper.npy")
t = data[:,0] # temps en secondes
Tz0 = data[:,1] # température à la surface (en z=0) en degrés Celsius

1/ Tracer la température à la surface en fonction du temps.

Corrigé

import matplotlib.pyplot as plt
plt.figure("Température à la surface")
plt.plot(t, Tz0)
plt.show()

On souhaite simuler la température dans le sol sur une profondeur de $10$ ⁠⁠ $m$ à l’aide de la méthode d’Euler explicite généralisée aux équations aux dérivées partielles.

La température à l’instant $𝑡_{𝑖} = 𝑡_{0} + 𝑖 Δ 𝑡$ à la profondeur $𝑧_{𝑗} = 𝑗 Δ 𝑧$ est notée $𝑇_{𝑖, 𝑗}$ . On choisit $101$ points de profondeur.

2/ Calculer le pas de profondeur $Δ 𝑧$ . La condition de stabilité de la méthode d’Euler $2 𝐷 Δ \frac{𝑡}{Δ} 𝑧^{2} < 1$ est-elle vérifiée ?

Coup de pouce 1

Considérons les valeurs (0,2,4). Combien y a-t-il de valeurs ? Quel pas sépare deux valeurs successives ? Quelle est l’étendue totale entre la première et la dernière valeur ? Quel lien relie le pas, l’étendue et le nombre de valeurs ?

Corrigé

Le pas de profondeur est donné par $Δ 𝑧 = \frac{10}{100} = 0,1 m$ .

2 𝐷 \frac{Δ 𝑡}{Δ 𝑧^{2}} = 2 \times 1

⋅10−6×36000,12=7,2⋅10−1<1

La condition de stabilité est bien vérifiée.

3/ Exprimer $𝑇_{𝑖, 𝑗}$ en fonction de $𝑇_{𝑖 - 1, 𝑗}$ , $𝑇_{𝑖 - 1, 𝑗 - 1}$ , $𝑇_{𝑖 - 1, 𝑗 + 1}$ , $𝐷$ , $Δ 𝑧$ et $Δ 𝑡$ grâce à la méthode d’Euler explicite.

Coup de pouce 1

Appliquer la formule de Taylor à l’ordre 2 à

𝑇_{𝑖, 𝑗 + 1}

𝑇_{𝑖, 𝑗 - 1}

Coup de pouce 2

Combiner les deux relations de Taylor pour exprimer

\frac{𝜕^{2} 𝑇}{𝜕 𝑥^{2}}

Coup de pouce 3

Si les indices ne correspondent pas à ce que l’on cherche, faire un changement d’indice pour trouver la relation demandée.

Corrigé

L’équation de diffusion est de la forme $\frac{𝜕 𝑇}{𝜕 𝑡} = 𝐷 \frac{𝜕^{2} 𝑇}{𝜕 𝑧^{2}}$ .

\begin{matrix} 𝑇_{𝑖, 𝑗 + 1} = 𝑇 (𝑡_{𝑖}, (𝑗 + 1) Δ 𝑧) \approx 𝑇_{𝑖, 𝑗} + Δ 𝑧 \frac{𝜕 𝑇}{𝜕 𝑧} + \frac{Δ 𝑧^{2}}{2} \frac{𝜕^{2} 𝑇}{𝜕 𝑧^{2}} \\ 𝑇_{𝑖, 𝑗 - 1} = 𝑇 (𝑡_{𝑖}, (𝑗 - 1) Δ 𝑧) \approx 𝑇_{𝑖, 𝑗} - Δ 𝑧 \frac{𝜕 𝑇}{𝜕 𝑧} + \frac{Δ 𝑧^{2}}{2} \frac{𝜕^{2} 𝑇}{𝜕 𝑧^{2}} \end{matrix}

D’où

𝑇_{𝑖, 𝑗 - 1} + 𝑇_{𝑖, 𝑗 + 1} \approx 2 𝑇_{𝑖, 𝑗} + Δ 𝑧^{2} \frac{𝜕^{2} 𝑇}{𝜕 𝑧^{2}}

Et donc

\frac{𝜕^{2} 𝑇}{𝜕 𝑧^{2}} \approx \frac{𝑇_{𝑖, 𝑗 - 1} + 𝑇_{𝑖, 𝑗 + 1} - 2 𝑇_{𝑖, 𝑗}}{Δ 𝑧^{2}}

Pour la dérivée temporelle, on a

\begin{matrix} 𝑇_{𝑖 + 1, 𝑗} = 𝑇 (𝑡_{0} + (𝑖 + 1) Δ 𝑡, 𝑗 Δ 𝑧) \approx 𝑇_{𝑖, 𝑗} + Δ 𝑡 \frac{𝜕 𝑇}{𝜕 𝑡} \\ \approx 𝑇_{𝑖, 𝑗} + Δ 𝑡 𝐷 \frac{𝜕^{2} 𝑇}{𝜕 𝑧^{2}} \\ \approx 𝑇_{𝑖, 𝑗} + \frac{Δ 𝑡}{Δ 𝑧^{2}} 𝐷 (𝑇_{𝑖, 𝑗 - 1} + 𝑇_{𝑖, 𝑗 + 1} - 2 𝑇_{𝑖, 𝑗}) \end{matrix}

L’énoncé demande $𝑇_{𝑖, 𝑗}$ , on fait donc le changement d’indice $𝑖 \to 𝑖 - 1$ :

𝑇_{𝑖, 𝑗} \approx 𝑇_{𝑖 - 1, 𝑗} + \frac{Δ 𝑡}{Δ 𝑧^{2}} 𝐷 (𝑇_{𝑖 - 1, 𝑗 - 1} + 𝑇_{𝑖 - 1, 𝑗 + 1} - 2 𝑇_{𝑖 - 1, 𝑗})

4/ Compléter le code suivant implémentant la méthode d’Euler explicite pour simuler la température dans le sol.

Nx = ... # nombre de points de profondeur
L = 10 # longueur de la colonne de sol simulée (en m)
x = np.linspace(..., ..., ...) # liste des profondeurs des points simulés (en m)

Nt = len(t) # nombre d'instants simulés
T = np.zeros((Nt,Nx)) + Tz0[0] # initialisation de la matrice de température
dt = 3600 # pas de temps (en s)
dx = L/(Nx-1) # pas de profondeur (en m)
D = 1e-6 # coefficient de diffusion (en m^2/s)
    
for i in range(1, Nt):
    T[i,0] = ... # condition à la surface
    for j in range(1, Nx-1):
        T[i,j] = ... # schéma d'Euler explicite
    T[i,-1] = T[i,-2] # le dernier point n'a pas de voisin de droite, on ne peut pas appliquer
                      # le schéma d'Euler explicite, on reprend la température du point
                      # précédent (condition de Neumann)

Coup de pouce 1

La température à la surface est contenue dans la variable Tz0.

Corrigé

Nx = 101
L = 10
x = np.linspace(0, L, Nx)

Nt = len(t)
T = np.zeros((Nt,Nx)) + Tz0[0]
dt = 3600
dx = L/(Nx-1)
D = 1e-6
    
for i in range(1, Nt):
    T[i,0] = Tz0[i]
    for j in range(1, Nx-1):
        T[i,j] = D * dt/dx**2 * T[i-1,j+1] + D * dt/dx**2 * T[i-1,j-1] + (1 - 2*D * dt/dx**2) *  T[i-1,j]
    T[i,-1] = T[i,-2]

5/ Tracer la température en fonction du temps pour des profondeurs de $0$ ⁠⁠ $m$ , $10$ ⁠⁠ $cm$ , $1$ ⁠⁠ $m$ et $10$ ⁠⁠ $m$ .

Coup de pouce 1

Quel indice correspond à la profondeur de

10

⁠

cm

? De

1

⁠

m

? De

10

⁠

m

Coup de pouce 2

On peut utiliser la notation T[:,j] pour accéder à toutes les lignes d’une colonne j et T[i,:] pour accéder à toutes les colonnes d’une ligne i.

Corrigé

for j in [0, 1, 10, 100]:
    plt.plot(t, T[:,j], label=f"z={x[j]:.2f} m")
plt.legend()
plt.show()

6/ Écrire une suite d’instructions permettant de calculer à quelle profondeur minimale il faut enterrer une conduite d’eau pour qu’elle ne subisse pas le gel.

Coup de pouce 1

On peut commencer par écrire des instructions pour détecter s’il gèle à une profondeur

𝑧_{𝑗}

Coup de pouce 2

On parcourt une colonne de T. Si tous les éléments sont positifs, alors il ne gèle jamais à cette profondeur. Sinon, il gèle au moins une fois à cette profondeur.

Corrigé

for j in range(Nx):
    gèle = False
    for i in range(Nt):
        if T[i,j] < 0:
            gèle = True # il gèle au moins une fois à la profondeur x[j]
            break # inutile de continuer, on passe à la profondeur suivante
    if not gèle: # on a trouvé une profondeur où il ne gèle jamais
        break # on s'arrête, pas la peine d'étudier les profondeurs suivantes
print(x[j])

7/ Compléter le code ci-dessous pour animer l’évolution du profil de température dans le sol. La fonction line.set_data prend les mêmes arguments que la fonction plot : les abscisses et les ordonnées des points à tracer.

import matplotlib.animation as animation
from datetime import datetime, timezone

def formateDate(t): # met en forme le temps pour afficher la date et l'heure
    return datetime.fromtimestamp(t, tz=timezone.utc).strftime("%d/%m/%Y %H:%M")

fig, ax = plt.subplots()
plt.clf()
ax = fig.add_subplot(111)
line, = ax.plot(..., ...) # tracé du profil de température à l'instant initial t0

ax.set_ylim(..., ...) # on choisit les limites de l'axe des ordonnées pour que la courbe soit bien visible
ax.set_xlabel("Profondeur (m)")
ax.set_ylabel("Température (°C)")
title = ax.set_title("")

def animate(i):
    line.set_data(..., ...)  # tracé du profil de température à l'instant ti
    title.set_text(formateDate(t[i]))
    return line, title

ani = animation.FuncAnimation(fig, animate, frames=range(0,Nt,10), interval=10, blit=False, repeat=False)
plt.show()

Corrigé

import matplotlib.animation as animation
from datetime import datetime, timezone

def formateDate(t): # met en forme le temps pour afficher la date et l'heure
    return datetime.fromtimestamp(t, tz=timezone.utc).strftime("%d/%m/%Y %H:%M")

fig, ax = plt.subplots()
plt.clf()
ax = fig.add_subplot(111)
line, = ax.plot(x, T[0,:])

ax.set_ylim(np.min(T), np.max(T))
ax.set_xlabel("Profondeur (m)")
ax.set_ylabel("Température (°C)")
title = ax.set_title("")

def animate(i):
    line.set_data(x, T[i,:])
    title.set_text(formateDate(t[i]))
    return line, title

ani = animation.FuncAnimation(fig, animate, frames=range(0,Nt,10), interval=10, blit=False, repeat=False)
plt.show()