Introduction

Les ingénieurs et urbanistes doivent-ils s'intéresser à la CFD en ville ?

Les vagues de chaleur ne sont plus des événements rares. La canicule européenne de 2003 a tué plus de 70 000 personnes ; celle de 2022 a battu des records de température en France, au Royaume-Uni et en Espagne dans la même semaine tandis que l'Europe fait de nouveau face à de multiples canicules records en mai, juin et juillet 2026. Ces épisodes deviennent le défi climatique définissant notre époque, transformant les environnements urbains denses en pièges thermiques où les températures peuvent dépasser celles des zones rurales environnantes de 5 à 10°C. Ce phénomène, l'effet d'îlot de chaleur urbain (ICU), tue des milliers de personnes chaque été et met les réseaux électriques à rude épreuve.

La question est : peut-on le simuler avec suffisamment de précision pour concevoir des villes plus résilientes face aux vagues de chaleur ?

Une étude publiée en 2022 dans Sustainability par Y. Qu, M. Milliez, L. Musson-Genon et B. Carissimo répond par un oui retentissant, en utilisant code_saturne, le solveur CFD open-source développé par EDF R&D. Leurs travaux simulent le cycle thermique diurne complet d'un quartier réel du centre-ville de Toulouse, en France, le valident contre un jeu de données expérimental exceptionnel, puis l'utilisent pour révéler comment la chaleur urbaine redistribue fondamentalement la dispersion des polluants.

Dans cet article, nous vous présentons leur méthodologie, leur modèle et, surtout, leurs campagnes expérimentales et leurs résultats obtenus avec code_saturne.

Le problème central : pourquoi l'hypothèse d'atmosphère neutre est dangereuse

La plupart des études micro-météorologiques sur les écoulements urbains et la dispersion des polluants supposent une atmosphère neutre. Cela signifie que les variations de température sont ignorées et que les forces de flottabilité sont nulles. Cette hypothèse est commode sur le plan numérique, mais physiquement erronée, en particulier lors de journées estivales où les toits exposés au soleil atteignent des températures bien supérieures à l'air ambiant. Qu et al. démontrent que cette hypothèse n'est pas une simple approximation ; c'est une erreur qui :

  • Fausse les températures de surface de plusieurs dizaines de degrés
  • Sous-estime l'énergie cinétique turbulente au niveau des toits
  • Modifie complètement la forme prédite des panaches de polluants

La solution qu'ils proposent est un modèle thermo-radiatif-convectif entièrement couplé, implémenté dans le module atmosphérique de code_saturne, appliqué à un quartier réel de Toulouse lors de la campagne de terrain CAPITOUL (Canopy and Aerosol Particles Interactions in TOulouse Urban Layer).

Présentation de la physique résolue par le modèle

La physique de la chaleur urbaine

Dans les villes, le transfert de chaleur fait principalement intervenir trois phénomènes qui ne peuvent pas être traités indépendamment :

  • Le transfert radiatif : le rayonnement solaire chauffe les surfaces, celles-ci ré-émettent un rayonnement infrarouge thermique, et les bâtiments échangent du rayonnement entre eux et avec le ciel.
  • Le transfert convectif : le vent transporte la chaleur à travers les canyons de rue et autour des bâtiments, avec un coefficient de transfert thermique variant spatialement sur chaque facette de paroi.
  • Le transfert conductif : la chaleur diffuse à travers les murs, les toits et le sous-sol.

La contribution centrale de cette étude est d'encoder ces trois phénomènes dans un cadre de simulation unique et cohérent.

Le bilan d'énergie de surface

Pour chaque surface du domaine urbain (chaque facette de mur, chaque tuile de toit, chaque segment de chaussée...): le modèle résout un Bilan d'Énergie de Surface (BES) :

Qcond+QH=QQ_{cond} + Q_H = Q^*

où :

  • QcondQ_{cond} (W.m2)(W.m^{-2}) est le flux de chaleur conductif reliant la température de surface à la température intérieure du bâtiment ou à la température profonde du sol.
  • QHQ_H (W.m2)(W.m^{-2}) est le flux de chaleur sensible échangé avec l'écoulement local.
  • QQ^* (W.m2)(W.m^{-2}) est le flux radiatif net (ondes courtes solaires + ondes longues émises par les bâtiments).

Le flux radiatif net QQ^* se décompose en ses composantes ondes courtes et ondes longues :

Q=S+LQ^* = S^* + L^*

où :

  • S=(1α)(SD+Sf+Se)S^* = (1 - \alpha)(S_D + S_f + S_e) (W.m2)(W.m^{-2}) est le flux d'ondes courtes net : l'énergie solaire effectivement absorbée par la surface après réflexion, avec α\alpha l'albédo de surface et SDS_D, SfS_f, SeS_e les composantes solaires directe, diffuse atmosphérique et réfléchie par l'environnement respectivement.
  • L=ε(La+LeσTsfc4)L^* = \varepsilon(L_a + L_e - \sigma T_{sfc}^4) (W.m2)(W.m^{-2}) est le flux d'ondes longues net : le bilan entre le rayonnement grandes longueurs d'onde reçu du ciel (LaL_a) et des bâtiments environnants (LeL_e), et le rayonnement émis par la surface elle-même en tant que corps gris (εσTsfc4\varepsilon \sigma T_{sfc}^4).

À noter que le flux de chaleur anthropique (émis par les activités humaines) et le flux de chaleur latente sont négligés dans cette étude. Cette simplification devrait être levée pour un modèle pleinement général d'une ville vivante, mais elle est appropriée pour les conditions estivales sèches de la campagne CAPITOUL.

Schéma du bilan d'énergie de surface en ville.
Schéma du bilan d'énergie de surface en ville.

Le soleil émet des rayonnements en ondes courtes, absorbés puis ré-émis en ondes longues par les bâtiments : ces échanges sont encodés dans QQ^*. La chaleur sensible QHQ_H est la chaleur échangée avec l'écoulement local, qui l'advecte ensuite localement. Enfin, les bâtiments et le sol ont leur propre température qui interagit avec l'air par des processus conductifs QcondQ_{cond}. Image tirée de [1]. Le flux de chaleur anthropique et le flux de chaleur latente (chaleur libérée ou absorbée lors des changements de phase) sont négligés ici.

Cette équation n'est pas résolue une seule fois pour une surface moyenne. Elle est résolue à chaque maille de surface du maillage, à chaque pas de temps, avec QHQ_H calculé à partir du champ de vent 3D local et QQ^* calculé par le solveur de transfert radiatif 3D. C'est ce qui rend le modèle véritablement tridimensionnel.

Du bilan thermique à la température de surface TsfcT_{sfc}

L'un des aspects les plus élégants de cet article est le schéma thermique hybride : la température du sol et la température de surface des bâtiments sont traitées par des modèles différents, chacun choisi pour correspondre au mieux à la physique sous-jacente.

Température du sol : la méthode force-restore

Pour le sol, les auteurs utilisent la méthode classique force-restore (Deardorff 1978), largement utilisée dans les modèles météorologiques :

Tsfct=2ωμt(L+SQH)ω(TsfcTg/b)\frac{\partial T_{sfc}}{\partial t} = \sqrt{2\omega} \, \mu_t (L^* + S^* - Q_H) - \omega(T_{sfc} - T_{g/b})

où :

  • TsfcT_{sfc} est la température de surface du sol.
  • ω=2π/86400\omega = 2\pi/86400 (s1)(s^{-1}) est la fréquence angulaire de la Terre.
  • LL^* (W.m2)(W.m^{-2}) est le flux radiatif net en ondes longues introduit ci-dessus.
  • SS^* (W.m2)(W.m^{-2}) est le flux radiatif net en ondes courtes introduit ci-dessus.
  • QHQ_H (W.m2)(W.m^{-2}) est le flux de chaleur sensible.
  • μt\mu_t (J.m2.s0.5.K1)(J.m^{-2}.s^{-0.5}.K^{-1}) est l'admittance thermique du sol.
  • Tg/bT_{g/b} (K)(K) est la température profonde du sol (ou intérieure du bâtiment) jouant le rôle de température de rappel.

Le terme de « force » (premier terme à droite) pilote la température de surface en réponse au forçage énergétique net. Le terme de « rappel » (le second) empêche la surface de s'écarter trop de sa température d'équilibre profonde aux échelles de temps diurnes. C'est une manière élégante, physiquement motivée et peu coûteuse en calcul de traiter l'inertie thermique du sol.

Température de surface des bâtiments : le modèle thermique de paroi

Pour les murs et les toits des bâtiments, un modèle de conduction monocouche est utilisé. Le bilan d'énergie à la surface d'une paroi est :

λe(TsfcTint)+hf(TsfcTa)=L+S=ε(La+LeσTsfc4)+(1α)(SD+Sf+Se)\frac{\lambda}{e}(T_{sfc} - T_{int}) + h_f(T_{sfc} - T_a) = L^* + S^* = \varepsilon(L_a + L_e - \sigma T_{sfc}^4) + (1 - \alpha)(S_D + S_f + S_e)

Décodons cette équation terme par terme ; c'est le cœur du modèle thermique :

Membre gauche :

  • λe(TsfcTint)\frac{\lambda}{e}(T_{sfc} - T_{int}) : flux conductif à travers la paroi, où λ\lambda (W.K1.m1)(W.K^{-1}.m^{-1}) est la conductivité thermique moyenne et ee (m)(m) l'épaisseur de la paroi. Ce terme lie la température de surface à la température intérieure du bâtiment TintT_{int}.
  • hf(TsfcTa)h_f(T_{sfc} - T_a) : flux convectif vers l'air extérieur, où hfh_f (W.m2.K1)(W.m^{-2}.K^{-1}) est le coefficient de transfert thermique local et TaT_a la température de l'air extérieur.

Membre droit :

  • L=ε(La+LeσTsfc4)L^* = \varepsilon(L_a + L_e - \sigma T_{sfc}^4) : flux radiatif net en ondes longues. LaL_a est le rayonnement atmosphérique descendant en ondes longues, LeL_e est le rayonnement en ondes longues provenant des surfaces des bâtiments environnants (l'échange mur-à-mur crucial), et εσTsfc4\varepsilon \sigma T_{sfc}^4 est l'émission propre de la surface en tant que corps gris.
  • S=(1α)(SD+Sf+Se)S^* = (1-\alpha)(S_D + S_f + S_e) : rayonnement solaire absorbé, où α\alpha est l'albédo de surface et SDS_D, SfS_f, SeS_e sont les composantes solaire directe, diffuse atmosphérique et réfléchie par l'environnement.

Cette équation est résolue implicitement pour TsfcT_{sfc} à chaque facette de surface et à chaque pas de temps, couplant simultanément le solveur de rayonnement, le solveur CFD et le modèle de conduction pariétale. Il faut à présent expliquer comment les termes de transfert radiatif LaL_a, LeL_e, SDS_D, SfS_f et SeS_e sont calculés.

Transfert radiatif QQ^*

Résolution du problème 3D complet par la méthode des ordonnées discrètes (DOM)

Les flux radiatifs sont calculés à l'aide de la méthode des ordonnées discrètes (DOM), qui résout l'équation du transfert radiatif (ETR) pour des milieux gris, non diffusants et semi-transparents. L'atout clé de l'implémentation de la DOM dans code_saturne est que l'ETR est résolue dans tout le domaine fluide, et non seulement aux faces solides comme ce serait le cas avec les méthodes de facteurs de forme. Cela signifie :

  • Le champ radiatif est calculé de manière volumique, ce qui rend l'approche naturellement extensible aux milieux non transparents (brouillard, aérosols, pollution). On suppose ici que l'atmosphère est entièrement transparente à l'échelle microscopique de ces simulations.
  • La discrétisation spatiale utilise le même maillage que le solveur CFD : aucune interpolation entre grilles différentes.
  • La discrétisation angulaire utilise 32 ou 128 directions discrètes, selon la précision requise.

Les rayonnements en ondes courtes (solaires) et en ondes longues (infrarouge thermique) sont traités séparément, ce qui est physiquement essentiel : les surfaces des bâtiments absorbent le rayonnement solaire avec un albédo spectral α\alpha et émettent le rayonnement thermique avec une émissivité ε\varepsilon, et ces deux grandeurs sont des propriétés matérielles indépendantes.

Rayonnement en ondes courtes : comment le soleil chauffe la ville

Le flux d'ondes courtes net absorbé par une surface est :

S=(1α)(SD+Sf+Se)S^* = (1 - \alpha)(S_D + S_f + S_e)

Un défi de modélisation central est le calcul des ombres : SDS_D est nul sur une surface ombrée par un bâtiment voisin. Le solveur DOM, opérant sur le maillage 3D complet, gère naturellement cela en propagant le faisceau solaire à travers le domaine et en laissant la géométrie le bloquer. C'est l'une des raisons pour lesquelles l'article utilise le même maillage pour les calculs CFD et radiatifs : les frontières d'ombre sont résolues à la même résolution spatiale que le champ d'écoulement (jusqu'à 0,5 m au centre du domaine).

À Toulouse, les auteurs ont en outre classifié les surfaces des bâtiments en quatre catégories d'albédo selon la couleur de la peinture murale, observée sur Google Maps : rose, gris, badigeon et blanc. Cela reflète une réalité physique : un mur peint en blanc exposé au soleil absorbe dramatiquement moins d'énergie solaire qu'un mur de brique sombre, et démontre qu'un détail de modélisation apparemment mineur (la couleur des murs) a un effet mesurable sur les températures de surface simulées.

outgoing city energy
outgoing city energy

Visualisation par le modèle de SS^\uparrow, à comparer et calibrer avec les données expérimentales.

Le flux solaire ascendant quittant une surface, mesurable par un pyranomètre sur le mât, est :

S=α(SD+Sf+Se)S^\uparrow = \alpha (S_D + S_f + S_e)

Cette grandeur est directement comparée aux mesures de la campagne CAPITOUL, avec un très bon accord (surestimation de 20 W.m2\sim 20\ W.m^{-2} à midi).

Rayonnement en ondes longues : les bâtiments qui rayonnent les uns vers les autres

En parallèle, chaque surface émet un rayonnement thermique (ondes longues, infrarouge) proportionnel à la puissance quatrième de sa température. Le flux d'ondes longues net à une surface est :

L=ε(La+LeσTsfc4)L^* = \varepsilon \left( L_a + L_e - \sigma T_{sfc}^4 \right)

LaL_a est le rayonnement descendant en ondes longues provenant du ciel (traité comme un corps noir à la température effective du ciel), LeL_e est le rayonnement en ondes longues reçu des surfaces des bâtiments environnants, et εσTsfc4\varepsilon \sigma T_{sfc}^4 est l'émission propre de la surface en tant que corps gris. Le terme d'échange mur-à-mur LeL_e est ce qui rend l'approche DOM si puissante : elle prend automatiquement en compte le fait qu'un mur orienté au sud chauffé par le soleil à 60°C irradie en permanence le mur opposé orienté au nord de l'autre côté de la rue, le maintenant plus chaud qu'il ne le serait isolément.

Le flux d'ondes longues ascendant quittant une surface est :

L=εσTsfc4L^\uparrow = \varepsilon \sigma T_{sfc}^4

Cela le rend extrêmement sensible à la température de surface : une erreur de 3°C sur TsfcT_{sfc} à 320 K se propage en une erreur d'environ 20 W.m220\ W.m^{-2} sur LL^\uparrow. La simulation montre une sous-estimation d'environ 20 W.m220\ W.m^{-2} la nuit et une surestimation allant jusqu'à 100 W.m2100\ W.m^{-2} à midi, toutes deux directement imputables aux biais correspondants de température de surface.

Température de brillance : ce que mesurent réellement les caméras infrarouges

Une innovation clé dans la stratégie de validation est l'utilisation de la température de brillance TbrT_{br} plutôt que de la température de surface réelle. Les caméras infrarouges aéroportées mesurent la radiance (et non directement la température), laquelle dépend à la fois de la vraie température de surface et de l'émissivité. Pour une surface de corps gris dans un canyon urbain, en tenant compte des réflexions infinies entre les parois, la température de brillance est :

Tbr=(εTsfc4+(1ε)Lσ)1/4T_{br} = \left( \varepsilon T_{sfc}^4 + \frac{(1-\varepsilon) L^\downarrow}{\sigma} \right)^{1/4}

ε\varepsilon est l'émissivité en ondes longues, σ=5,67×108\sigma = 5{,}67 \times 10^{-8} (W.m2.K4)(W.m^{-2}.K^{-4}) est la constante de Stefan-Boltzmann, et LL^\downarrow (W.m2)(W.m^{-2}) est le rayonnement incident en ondes longues. La température de brillance est toujours inférieure à la vraie température de surface lorsque ε<1\varepsilon < 1. L'article utilise une approximation communément admise, valable lorsque le terme de réflexion (1ε)L/σ(1-\varepsilon)L^\downarrow/\sigma est petit devant εTsfc4\varepsilon T_{sfc}^4, c'est-à-dire lorsque TsfcT_{sfc} est élevée ou (1ε)(1-\varepsilon) faible :

TbrTsfcε1/4T_{br} \approx T_{sfc} \cdot \varepsilon^{1/4}
infrared city data vs simulation
infrared city data vs simulation

L'article montre que par un après-midi d'été, cet écart atteint environ 3°C sur les toits : non négligeable lors d'une validation contre des images infrarouges aéroportées. Comparer des températures de surface brutes à des mesures IRT sans cette correction surestimerait systématiquement les températures.

hf=ρCpuκfhσtln ⁣(d+z0z0T)fmh_f = \frac{\rho C_p u_* \kappa f_h}{\sigma_t \ln\!\left(\dfrac{d + z_0}{z_{0T}}\right) \sqrt{f_m}}

où :

  • CpC_p est la chaleur spécifique de l'air (J.kg1.K1)(J.kg^{-1}.K^{-1}).
  • uu^* est la vitesse de frottement locale (m.s1)(m.s^{-1}), calculée par le modèle kkε\varepsilon.
  • κ=0,41\kappa = 0{,}41 est la constante de von Kármán.
  • σt\sigma_t est le nombre de Prandtl turbulent.
  • dd (m)(m) est la distance du centre de la maille à la paroi.
  • z0z_0 et z0Tz_{0T} sont les longueurs de rugosité dynamique et thermique (m)(m).
  • fmf_m et fhf_h sont les fonctions de stabilité de Louis (1979), qui valent 1 en conditions neutres.

Ceci est implémenté via une loi logarithmique modifiée qui tient compte de la stratification atmosphérique à travers les fonctions de stabilité : une approche explicite qui évite la méthode itérative de similarité de Monin-Obukhov, la rendant numériquement tractable en 3D. La conséquence physique de ce choix est profonde : un toit exposé au soleil dans une zone venteuse perdra sa chaleur beaucoup plus efficacement qu'un mur abrité dans une zone de recirculation. Seul un solveur d'écoulement 3D peut capturer cette variabilité.

Le modèle de turbulence : kkε\varepsilon dans la canopée urbaine

Les auteurs utilisent une approche RANS avec une fermeture turbulente kkε\varepsilon. Le choix est délibéré et honnête : le modèle kkε\varepsilon ne peut pas résoudre les grands tourbillons dépendants de la géométrie que capturerait une LES, et il surestime l'énergie dissipée. Cependant, il fournit une précision acceptable à un coût de calcul raisonnable pour ce type de simulation diurne à l'échelle du quartier. La viscosité turbulente est :

μt=ρCμk2ε\mu_t = \rho C_\mu \frac{k^2}{\varepsilon}

avec Cμ=0,09C_\mu = 0{,}09. La caractéristique critique pour les simulations thermiques urbaines est le terme de production par flottabilité dans l'équation de l'énergie cinétique turbulente :

Gb=μtρσtβgiρxiG_b = -\frac{\mu_t}{\rho \sigma_t} \beta g_i \frac{\partial \rho}{\partial x_i}

Ce terme est ce qui couple le champ de température au champ de turbulence. Lorsque les surfaces sont chauffées, Gb>0G_b > 0 et la turbulence est amplifiée par la flottabilité: en particulier au niveau des toits où l'air chaud monte depuis les surfaces chauffées. Ce n'est pas un terme correctif : dans des conditions de vent faible et d'ensoleillement intense, il peut être le mécanisme de production dominant de turbulence dans la canopée.

Flux de chaleur sensible : le moteur turbulent de l'atmosphère urbaine

Le flux de chaleur sensible au niveau du mât est estimé à partir de la solution kkε\varepsilon comme :

QH=ρCpwθ=ρCpKtθzQ_H = \rho C_p \overline{w'\theta'} = \rho C_p K_t \frac{\partial \theta}{\partial z}

où la diffusivité turbulente Kt=νt/PrtK_t = \nu_t / Pr_t est directement disponible depuis le modèle de turbulence. L'accord modèle-observation est bon en journée (différences <40 W.m2< 40\ W.m^{-2}) mais présente un biais nocturne de 45-70 W.m2W.m^{-2}. L'explication : la nuit, avec un vent calme et une faible différence de température surface-air, le flux sensible observé est quasi nul. Cependant, la température de toit initialisée dans le modèle crée un gradient de température non nul, forçant un flux de chaleur non nul (ici spurieusement négatif).

Vitesse de frottement : relier la turbulence à la quantité de mouvement

La vitesse de frottement à la surface est définie par :

u=(τwρ)1/2u^* = \left(\frac{|\tau_w|}{\rho}\right)^{1/2}

et au niveau du mât est estimée à partir des flux de quantité de mouvement cinématiques dans les deux directions horizontales :

u=(uw2+vw2)1/4u^* = \left(u'w'^2 + v'w'^2\right)^{1/4}

où :

uw=νtuz,vw=νtvzu'w' = \nu_t \frac{\partial u}{\partial z}, \quad v'w' = \nu_t \frac{\partial v}{\partial z}

Le modèle capture correctement l'ordre de grandeur de uu^* dans la plupart des conditions, avec des difficultés apparaissant aux extrêmes (valeurs très faibles ou très élevées). C'est une limitation connue des fermetures RANS dans les géométries urbaines complexes.

Le domaine de calcul : le centre-ville de Toulouse, pas une géométrie jouet

L'étude simule un quartier urbain réel de 891×963×200891 \times 963 \times 200 m centré sur l'intersection des rues Alsace-Lorraine et Pomme dans le centre de Toulouse, un centre-ville européen dense du XIXe siècle avec :

  • Des bâtiments de 4 à 5 étages (environ 20 m de haut)
  • Des murs principalement en brique et en pierre
  • Des toits en tuile de terre cuite
  • Une végétation très clairsemée

La base de données urbaine 3D a été fournie par la ville de Toulouse au format AutoCAD et prétraitée dans ICEM CFD. Ce prétraitement était non trivial : la base de données brute contenait des murs intérieurs, des hauteurs variables, pas de surface de sol, et un excès de détails fins. Les étapes de préparation ont inclus :

  • La suppression des surfaces intérieures.
  • La simplification des détails architecturaux.
  • La création du plan de sol et la projection des bâtiments sur celui-ci.
  • L'application d'une stratégie de détail géométrique progressif : détail complet au centre Alsace-Lorraine/Pomme, blocs urbains simplifiés dans les zones environnantes, et une valeur de rugosité élevée aux limites du domaine.
Le maillage discrétisé représentant la ville de Toulouse
Le maillage discrétisé représentant la ville de Toulouse

Le maillage résultant est une grille tétraédrique non structurée d'environ 4,5 millions de mailles, avec une résolution allant de 0,5 m près du centre à 10 m aux bords du domaine. De plus, dix thermomètres infrarouges fixes (IRT) ont été fixés à des balcons et des bras à travers trois canyons de rue (Alsace-Lorraine, Pomme et Rémusat), enregistrant les températures de brillance de surface des murs, chaussées et toits à des positions fixes. Enfin, trois vols IRT aéroportés ont été effectués le 15 juillet 2004 à bord d'un Piper Aztec PA23 volant à 460 m d'altitude, produisant des images thermiques du quartier à une résolution de 1,5 à 3 m :

  • Vol 430 : 0749–0816 UTC (matin)
  • Vol 431 : 1115–1150 UTC (fin de matinée)
  • Vol 432 : 1348–1423 UTC (début d'après-midi)

Ce jeu de données est exceptionnellement riche : il combine des informations spatiales (imagerie aéroportée) et temporelles (séries chronologiques du mât, capteurs IRT fixes), permettant de valider à la fois la distribution spatiale des températures et leur évolution diurne.

Comparaison IRT : le couplage complet l'emporte

Le test le plus direct du modèle est de comparer le champ de température de brillance simulé aux images IRT aéroportées. Trois configurations de modélisation sont comparées au vol 431 (1138 UTC) :

Configuration du modèleErreur température toitErreur mursErreur rue
Rayonnement seul (sans convection)> +25°C (hors échelle)GrandeGrande
Rayonnement + hfh_f constant =12 W.m2.K1= 12\ W.m^{-2}.K^{-1}Surrefroidissement ~15°C4–5°C4–7°C
Couplage radiatif-convectif complet< 10°C partout< 5°C< 3°C en moyenne

Le résultat de la première configuration (rayonnement sans refroidissement convectif) est catastrophiquement faux : sans QHQ_H, toute l'énergie radiative absorbée s'accumule en surface sans voie d'échappement, portant les températures de surface à plus de 25°C au-dessus des observations. Cela seul justifie l'approche couplée.

La seconde configuration (hfh_f constant =12 W.m2.K1= 12\ W.m^{-2}.K^{-1}) est une amélioration significative, mais introduit un surrefroidissement systématique au niveau des toits d'environ 15°C. La raison est simple : la valeur constante choisie surestime le coefficient de transfert thermique réel spatialement moyenné à ce moment de la journée, retirant trop d'énergie de la surface. Sur les murs (ombrés ou ensoleillés), l'erreur est de 4–5°C, et sur les rues de 4–7°C.

Seul le couplage radiatif-convectif complet, avec un hfh_f spatialement variable dérivé à chaque facette de surface à partir de la solution d'écoulement 3D locale, réduit les erreurs à moins de 10°C sur les toits, moins de 5°C sur les murs, et moins de 3°C en moyenne pour les rues. Les erreurs résiduelles sur les rues s'expliquent en partie par la résolution limitée du maillage dans les passages les plus étroits (un minimum de trois mailles par largeur de rue), qui approxime les frontières d'ombre plutôt que de les résoudre exactement.

Un résultat clé, illustré par la comparaison entre les champs de hfh_f constant et variable, est que le coefficient de transfert thermique varie de plus d'un facteur trois à travers le domaine à un instant donné, selon l'exposition locale au vent. Une valeur constante unique ne peut pas représenter cette variabilité, quelle que soit la soin avec lequel elle est choisie.

Évolution diurne de la température de brillance TbrT_{br}

Les capteurs IRT fixes permettent de valider l'évolution temporelle de la température de brillance à des positions spécifiques tout au long de la journée du 15 juillet 2004. Dans l'ensemble, le modèle reproduit correctement le cycle diurne à toutes les positions de capteurs, mais trois écarts instructifs méritent d'être examinés.

Murs du canyon Alsace-Lorraine (faces est et ouest) : le modèle montre un surrefroidissement d'environ 5°C pendant la soirée (1800–2400 UTC) sur les deux faces. Sur la face ouest, le modèle prédit une TbrT_{br} plus élevée qu'observée de 0600 à 1200 UTC. Cette surestimation matinale a également été trouvée par Hénon et al. (2012) avec le modèle SOLENE sur le même jeu de données, et est attribuée à une possible sous-estimation par le capteur IRT lui-même.

Chaussée Alsace-Lorraine : le biais sur la température de brillance de la chaussée peut s'expliquer par l'approximation des frontières d'ombre dans le maillage. La valeur extraite d'une seule maille sélectionnée peut différer significativement de ses voisines lorsque le bord d'ombre tombe à proximité, rendant la comparaison point-à-point sensible au positionnement des ombres à l'échelle sous-maille.

Mur nord-est de Pomme : le désaccord principal se produit en après-midi, où le biais atteint 8°C. D'après des photographies de ce mur (Figure 11 de l'article), la façade est couverte de nombreuses fenêtres équipées de stores blancs. Le capteur IRT peut avoir détecté une position où les stores étaient fermés et donc surchauffés, ou inversement où des fenêtres étaient ouvertes et la ventilation intérieure abaissait la température de surface extérieure. Ni l'un ni l'autre de ces effets n'est représenté dans le modèle, qui traite le mur comme une surface opaque uniforme. Cela rappelle que les hétérogénéités sous-façade dans les bâtiments réels introduisent un bruit irréductible dans la validation point-à-point.

Toit : le modèle montre un excellent accord avant le lever du soleil et de la mi-après-midi jusqu'à minuit. Cependant, il présente une avance systématique dans le réchauffement matinal : le toit simulé commence à chauffer plus tôt qu'observé. Le capteur était monté sur un toit qui, en réalité, est entouré de petits obstacles (cheminées, gaines de ventilation, acrotères) visibles sur des photographies aériennes mais absents du modèle 3D simplifié. Ces structures ombragent le capteur lors des premières heures d'exposition au soleil. Une fois que le soleil s'est suffisamment levé pour les dépasser, le biais disparaît (cohérent avec l'hypothèse que l'erreur est purement géométrique et non un problème de physique du modèle).

Flux radiatifs sortants : LL^\uparrow et SS^\uparrow

Le mât météorologique mesure à la fois le flux d'ondes longues sortant LL^\uparrow et le flux solaire sortant SS^\uparrow, fournissant une vérification directe du modèle radiatif.

Flux d'ondes longues LL^\uparrow : le flux infrarouge sortant est directement lié à la température de surface via L=εσTsfc4L^\uparrow = \varepsilon \sigma T_{sfc}^4. En raison de la dépendance en puissance quatre, les erreurs sur TsfcT_{sfc} sont amplifiées dans LL^\uparrow : une erreur de 3°C à 320 K correspond à environ 20 W.m220\ W.m^{-2}. Le modèle montre :

  • Une sous-estimation d'environ 20 W.m220\ W.m^{-2} la nuit, directement imputable au biais de température de surface nocturne discuté ci-dessus.
  • Une surestimation allant jusqu'à 100 W.m2100\ W.m^{-2} à midi, imputable à l'avance dans le réchauffement matinal des toits.

La différence entre LL^\uparrow modélisé calculé à la surface du toit et au niveau du mât est généralement inférieure à 20 W.m220\ W.m^{-2}, confirmant que la variation verticale du flux infrarouge sortant entre la surface et 47,5 m est faible dans une atmosphère transparente.

Au-delà de la comparaison scalaire au mât, le modèle permet une visualisation 3D complète de la distribution de LL^\uparrow dans tout le domaine (Figure 14 de l'article). Cela révèle une variabilité spatiale significative : les surfaces horizontales (toits, chaussées) sont nettement plus chaudes que les surfaces verticales (murs) en journée, et émettent donc beaucoup plus de rayonnement infrarouge. Cela ne peut pas être capturé par tout modèle qui représente la surface urbaine comme une couche unique moyennée.

Flux solaire SS^\uparrow : pour le toit gris à la position du mât, un albédo α=0,15\alpha = 0{,}15 a été attribué, cohérent avec la littérature. L'accord entre modèle et observation est très bon tout au long de la journée, avec une surestimation d'environ 20 W.m220\ W.m^{-2} seulement à midi. Cela valide à la fois le schéma de rayonnement solaire et la stratégie de classification des albédos. La visualisation 3D de SS^\uparrow (Figure 16 de l'article) rend également immédiatement visibles les surfaces à albédo plus élevé (plus brillantes sur l'image) et les frontières d'ombre nettes projetées par les arêtes des bâtiments (deux sorties directes du solveur DOM sur le maillage 3D complet).

Flux de chaleur sensible QHQ_H et vitesse de frottement uu^*

Flux de chaleur sensible : le flux de chaleur sensible à la surface du toit est comparé à la mesure par covariance des turbulences sur le mât. L'accord est bon en journée, avec des différences généralement inférieures à 40 W.m240\ W.m^{-2}. Le principal désaccord se produit la nuit : de minuit à 0700 UTC, le QHQ_H observé est quasi nul (l'air juste au-dessus du toit est à une température très proche de celle du toit lui-même). Cependant, le modèle initialise la température du toit à une valeur mesurée inférieure à la température de l'air à cette heure, produisant un flux de chaleur sensible spurieusement négatif (c'est-à-dire de la chaleur s'écoulant de l'air vers la surface). Le biais nocturne de 45–70 W.m2W.m^{-2} est donc une conséquence de la procédure d'initialisation plutôt qu'un échec de la physique du modèle. De 1300 à 1500 UTC, un désaccord secondaire apparaît que les auteurs attribuent à l'influence d'une géométrie de toit non résolue sur l'échange convectif local.

Vitesse de frottement : la vitesse de frottement uu^* fournit une vérification du transfert de quantité de mouvement, indépendante du modèle thermique. Le modèle capture le bon ordre de grandeur tout au long de la journée, avec un bon accord lorsque u<0,2 m.s1u^* < 0{,}2\ m.s^{-1}. À des valeurs plus élevées, le modèle tend à surestimer uu^* au niveau du mât pendant la nuit (1800–2400 UTC). C'est une limitation connue de la fermeture RANS kkε\varepsilon dans les géométries complexes : elle surestime le transfert de quantité de mouvement turbulent dans les écoulements décollés et recirculants, et les estimations de uu^* par anémométrie sonique sont elles-mêmes incertaines lorsque les valeurs sont faibles. La différence entre uu^* calculé à la surface du toit et au niveau du mât est également partiellement un artefact numérique de la méthode d'estimation des gradients.

L'application critique : comment les effets thermiques redistribuent la dispersion des polluants

Ayant établi la confiance dans le modèle thermique grâce aux comparaisons CAPITOUL, les auteurs utilisent le même cadre pour étudier une question de santé publique directe : que se passe-t-il pour un polluant rejeté au niveau du sol dans un centre-ville par une chaude journée d'été, par rapport à l'hypothèse standard d'atmosphère neutre ?

Deux simulations sont réalisées avec des conditions de vent identiques (Uref=3 m.s1U_{ref} = 3\ m.s^{-1} vers l'est à 47,5 m, profil de température neutre à 22°C), ne différant que dans le traitement des effets thermiques :

  • Cas neutre : pas d'effets thermiques, pas de chauffage de surface, pas de flottabilité.
  • Cas thermo-radiatif : couplage 3D radiatif-convectif-conductif complet actif, forçage solaire à partir de 1000 UTC pendant 3 heures.

Un traceur passif est rejeté au niveau du sol au centre du domaine pendant 3 heures.

Effet sur l'énergie cinétique turbulente : dans le cas neutre, la production de TKE est purement mécanique, confinée à l'immédiate proximité des arêtes des bâtiments où l'écoulement se décolle. Avec le couplage thermique actif, les surfaces chauffées injectent continuellement de l'énergie de flottabilité dans l'écoulement. La TKE augmente dans tout le domaine, pas seulement près des bâtiments, et est particulièrement amplifiée au niveau des toits où le terme de production par flottabilité GbG_b est le plus grand. L'écoulement n'atteint jamais un état stationnaire dans le cas thermo-radiatif : les températures de surface continuellement évolutives modifient continuellement le forçage de flottabilité, rendant l'écoulement intrinsèquement instationnaire même sous des conditions aux limites de vent constantes.

Effet sur la concentration au niveau du sol : dans le cas neutre, le panache se déplace vers l'est avec le vent dominant. À l'intérieur des cours et des canyons de rue étroits, les vitesses du vent sont drastiquement réduites, créant des zones de stagnation où le polluant s'accumule et est mal dilué. De fortes concentrations au sol persistent bien en aval de la source.

ground pollution in Toulouse; neutral atmosphere vs full radiative scheme
ground pollution in Toulouse; neutral atmosphere vs full radiative scheme

Avec le transfert thermique actif, le chauffage de surface entraîne un mouvement ascendant buoyant persistant. Le panache s'élève depuis le point d'émission, les concentrations près du sol sont plus faibles en aval, et les concentrations dans les cours sont réduites. Cependant, les concentrations au-dessus du niveau des toits sont plus élevées : la poussée thermique qui protège les piétons au niveau de la rue transfère la pollution vers le haut dans la couche où elle peut être advectée sur de plus longues distances.

Ce ne sont pas de subtiles différences quantitatives. Elles représentent des régimes de dispersion qualitativement différents : le modèle neutre prédit un panache rampant au sol piégé dans les canyons de rue, tandis que le modèle thermiquement couplé prédit un panache ascendant s'échappant de la canopée. Si vous êtes un gestionnaire de qualité de l'air urbain évaluant le risque d'un rejet accidentel au niveau du sol par un chaud après-midi d'été, l'hypothèse neutre vous donne la mauvaise géométrie de la zone de danger : pas seulement de mauvaises concentrations, mais une mauvaise distribution spatiale, une mauvaise population affectée et une mauvaise stratégie d'atténuation.

De la simulation à la conception : ce que cela permet aux urbanistes

La validation CAPITOUL donne confiance dans le cadre thermo-radiatif-convectif de code_saturne comme laboratoire numérique pour les interventions de conception urbaine. Les auteurs eux-mêmes énumèrent les applications :

🏗️ Surfaces à albédo élevé (« toits blancs »)

Dans le modèle thermique de paroi, changer l'albédo α\alpha modifie directement le rayonnement en ondes courtes absorbé :

(1α)(SD+Sf+Se)(1-\alpha)(S_D + S_f + S_e)

Augmenter α\alpha de 0,15 (toit gris, tel que mesuré à Toulouse) à 0,7 (revêtement de toit blanc) réduit l'énergie solaire absorbée par la surface de 65 %. Cela a un effet en cascade : température de surface plus basse → moins de QHQ_H dans l'air de la rue (refroidissement diurne) et moins de QcondQ_{cond} stocké dans la paroi (réduction du relâchement nocturne de chaleur) → température de l'air plus basse le jour comme la nuit. Le modèle peut quantifier cette cascade pour un quartier spécifique avec sa géométrie et son orientation propres.

🌿 Toits verts et espaces verts

Le cadre peut être étendu pour inclure le terme de flux de chaleur latente (actuellement négligé) afin de représenter la végétation. L'évapotranspiration retire de l'énergie de la surface sous forme de chaleur latente plutôt que sensible, refroidissant directement l'air local. Le modèle fournit les champs de vent et de turbulence spatiaux nécessaires pour calculer correctement où ce refroidissement a le plus d'impact — ce qui est fonction des schémas d'écoulement locaux qui varient énormément à l'intérieur d'un seul pâté de maisons.

🌬️ Ventilation des canyons de rue

En exécutant la simulation d'écoulement 3D complète, le modèle identifie quelles rues sont bien ventilées (fort uu^* local, fort hfh_f, évacuation efficace de la chaleur) et lesquelles sont des pièges thermiques (faible vitesse du vent, faible hfh_f, accumulation de chaleur). Cette information peut directement éclairer les décisions sur les ratios hauteur/largeur des bâtiments, les largeurs de rue et l'orientation des bâtiments dans les nouveaux aménagements.

☁️ Points chauds de qualité de l'air

Comme le montrent les résultats de dispersion des polluants, les effets thermiques modifient fondamentalement où se concentre la pollution. Le modèle couplé peut identifier des emplacements spécifiques dans un centre-ville où les émissions au niveau du sol sont le plus susceptibles de s'accumuler dans des conditions typiques de vague de chaleur estivale, et évaluer si une chaussée fraîche ou une rangée d'arbres changerait significativement la qualité de l'air local.

Pourquoi code_saturne ? Une référence de la CFD scientifique open-source

code_saturne n'est pas une boîte noire commerciale. Développé et maintenu par EDF R&D, il est :

Open-source (licence GNU GPL) : chaque équation décrite dans cet article est inspectable, modifiable et entièrement auditable

Volumes finis sur maillages non structurés : gère les géométries extrêmement complexes des centres-villes réels sans simplification

Physiquement complet : le module atmosphérique inclut la turbulence RANS (kkε\varepsilon) et LES, le rayonnement DOM (ondes courtes et longues séparément), le modèle force-restore pour le sol, le modèle thermique de paroi, et la dispersion de scalaire passif, le tout entièrement couplé

Massivement parallèle : passe à l'échelle de milliers de cœurs CPU, rendant les simulations diurnes à l'échelle du quartier tractables

Validé pour les environnements urbains réels : l'étude CAPITOUL est la validation la plus complète à ce jour, mais s'appuie sur une validation préalable contre la campagne de terrain MUST et des cas idéalisés

Gratuit : ce qui signifie que les villes de toutes tailles, les pays en développement et les groupes de recherche académique peuvent accéder aux mêmes outils qu'EDF et le CEREA

Premiers pas : votre première simulation thermique urbaine

Voici une feuille de route pratique pour reproduire ce type d'étude pour votre propre ville :

  • Géométrie : obtenez la base de données urbaine 3D auprès de votre municipalité (de nombreuses villes européennes ont cela au format CityGML ou similaire). Importez-la dans ICEM CFD ou Salome-Meca (open-source). Effectuez le nettoyage de la géométrie : supprimez les murs intérieurs, créez la surface du sol, projetez les bâtiments et simplifiez le détail progressivement du bord du domaine vers le centre d'intérêt.
  • Maillage : construisez un maillage tétraédrique non structuré. Visez une résolution de ~0,5 m dans votre zone d'intérêt, s'épaississant jusqu'à ~10 m aux limites du domaine. Un nombre total de 3 à 5 millions de mailles est un point de départ pratique pour un quartier de ~1 km².
  • Propriétés des matériaux : classifiez vos surfaces de bâtiments par albédo et émissivité en vous basant sur une inspection visuelle ou les archives municipales. Distinguez au minimum les matériaux de toiture (tuile, gravier, bitume) et les couleurs des murs. Définissez la conductivité thermique de la paroi λ\lambda et l'épaisseur ee à partir de la littérature pour votre type de bâtiment : ceux-ci contrôlent la résistance thermique conductrice R=e/λR = e/\lambda et donc l'inertie thermique du bâtiment.
  • Conditions aux limites atmosphériques : utilisez les profils de vitesse du vent, de direction et de température mesurés à la station météorologique la plus proche, mis à jour toutes les 2 heures. Définissez la stabilité atmosphérique à partir des profils de température potentielle.
  • Activer le couplage complet : dans le module atmosphérique de code_saturne, activez le solveur de transfert radiatif DOM (32 directions pour un premier essai, 128 pour la validation), le modèle force-restore pour le sol, et le modèle thermique de paroi. Définissez la position solaire à partir de vos coordonnées géographiques et de la date de simulation.
  • Valider : comparez avec toutes les mesures de température de surface disponibles, idéalement incluant une station météorologique sur un toit. Si possible, utilisez l'imagerie thermique infrarouge par satellite ou drone pour une validation spatiale. N'oubliez pas de comparer les températures de brillance, pas les températures de surface brutes.
  • Exploration de conception : une fois validé, changez une variable à la fois (albédo du toit, matériau de surface, canopée arborée, hauteur de bâtiment) et quantifiez l'effet sur la température de surface, le flux de chaleur sensible et la dispersion des polluants. C'est votre laboratoire urbain numérique.

Conclusion : la ville comme équation résoluble

Les travaux de Qu, Milliez, Musson-Genon et Carissimo sont bien plus qu'une étude académique. C'est une preuve de concept que l'environnement thermique urbain est quantitativement prévisible en trois dimensions, à l'échelle des bâtiments, sur un cycle diurne complet, dans un vrai centre-ville européen.

La physique est encodée dans le bilan d'énergie de surface. Les échanges radiatifs sont encodés dans la méthode des ordonnées discrètes. La température de surface émerge des modèles couplés force-restore et thermique de paroi. La conduction pariétale relie la température de surface extérieure à la masse thermique du bâtiment, encodant le lent relâchement nocturne de chaleur qui entretient l'îlot de chaleur urbain bien après le coucher du soleil. Le refroidissement convectif est calculé localement à chaque facette de surface à partir de la solution d'écoulement 3D. Et code_saturne relie tout cela dans un cadre de simulation résoluble, reproductible et open-source.

Le résultat peut-être le plus préoccupant de l'article est celui sur la dispersion des polluants : supposer une atmosphère neutre — comme le fait encore aujourd'hui la grande majorité des études d'écoulement urbain — vous donne la mauvaise forme de panache, les mauvaises concentrations au sol et la mauvaise identification des zones à risque. Par une chaude journée d'été dans un centre-ville, la flottabilité n'est pas une perturbation. C'est un effet de premier ordre.

Nous avons les outils. Nous avons la physique. Nous avons le logiciel open-source. Nous avons même les données de validation d'une campagne de terrain d'un an à Toulouse. Ce dont nous avons besoin maintenant, c'est de plus d'ingénieurs, d'urbanistes et de chercheurs prêts à lancer les simulations et à apporter des preuves quantitatives, fondées sur la physique, dans les salles où se prennent les décisions de conception urbaine.

Car lorsque la prochaine vague de chaleur tuera des gens dans les centres-villes — et elle le fera — la question à laquelle nous devrons répondre ne sera pas « avions-nous les outils pour comprendre cela ? » (nous les avions clairement) mais « les avons-nous utilisés ? »

Références

📌 Article complet : Qu Y., Milliez M., Musson-Genon L., Carissimo B. — « Modelling Radiative and Convective Thermal Exchanges over a European City Center and Their Effects on Atmospheric Dispersion » — Sustainability, 2022, 14(12), 7295. DOI : 10.3390/su14127295

📌 Jeu de données expérimental CAPITOUL : Masson V. et al. — « The Canopy and Aerosol Particles Interactions in TOulouse Urban Layer (CAPITOUL) Experiment » — Meteor. Atmos. Phys., 2008, 102, 135–157.

[1] : Butera, Federico. (2019). Energy and Resource Efficient Urban Neighbourhood Design Principles for Tropical Countries — A Practitioner's Guidebook.

[2] : Pigeon, G., Moscicki, M. A., Voogt, J., & Masson, V. (2008). Simulation of fall and winter surface energy balance over a dense urban area using the TEB scheme. Archiv für Meteorologie, Geophysik und Bioklimatologie, Serie B, 102, 159–171.