CFD Urbaine

Modélisation de la turbulence

Approches globales

En théorie, il est possible de résoudre tout problème de mécanique des fluides grâce aux équations de Navier-Stokes . Cependant cette approche (appelé DNS – Direct Numerical Simulation ), est aujourd’hui encore beaucoup trop coûteuse en ressources numériques, et demanderait d’avoir une discrétisation spatiale extrêmement fine pour pouvoir capter toutes les échelles de turbulences et ce sur l’intégralité du domaine étudié.

La première solution pour réduire cette demande en calcul est d’avoir recours à la méthode appelée LES – Large Eddy Simulation . L’idée est de ne simuler (représenter) que les échelles les plus larges de turbulence (les plus “gros tourbillons”), et prendre en compte ce qui se passe aux échelles inférieures à l’aide de modèles (voir figure ci dessous). Cependant de nombreuses discussions sont encore en cours sur le type de modélisation optimal [Sagaut, 2006]¹.

Figure D’après [Lu, 2008]²

La méthode la plus largement utilisé pour la résolution des écoulements turbulents est aujourd’hui la méthode RANS – Reynolds Averaged Navier Stokes [Frank et al., 2004]³. Cette méthode moyenne temporellement les équations d’écoulements sur toutes les échelles de turbulences. De façon analogue aux méthodes LES, il convient ensuite de modéliser les effets de la turbulence sur le flux moyen : c’est le rôle des modèles de turbulence discutés dans la suite.

Les simulations se font en régime établi (Steady State), du fait que les données météo sont moyennées au minimum sur des périodes allant de 15 à 30 minutes.

Modélisation de la turbulence

Au sein des approches de types RANS, les modèles de turbulence ont pour but de définir les variables turbulentes du flux en fonction des variables globales du flux moyen. Deux grandes familles d’approches existent : une première est fondée sur la viscosité des tourbillons, et dérive les contraintes turbulentes à partir de la vitesse moyenne ; et la seconde considère des équations différentielles supplémentaires.

Historiquement, l’approche la plus utilisé reste l’approche linéaire $k - \epsilon$ :

• $k$ : énergie cinétique turbulente,
• $\epsilon$ : dissipation de l’énergie cinétique turbulente.

Un problème connu de ce modèle est la surproduction d’énergie cinétique turbulente dans les zones stagnantes. De nombreux auteurs ont tenté d’améliorer ce modèle, ce qui a permis de façon générale d’améliorer la prédiction des coefficients de pression ($C_p$), mais diminue la qualité des prédictions de vitesses, en particulier à l’arrière des obstacles.

Ce modèle, et les différentes variantes qui ont suivi peuvent se retrouver ici :

• $k - \epsilon$ :
  • Problèmes de “surplus d’énergie” dans les zones stagnantes
• Modifications du $k - \epsilon$ :
  • Bons coefficients de pressions sur les façades, mais mauvaises prédictions des vitesses en aval des obstacles
• RNG $k-\epsilon$ (ReNormalization Group) ou Realizable $k - \epsilon$ :
  • Atténuation des problèmes de stagnation
  • Meilleurs résultats en aval des obstacles
  • Semble fonctionner le mieux pour la simulation d’environnements complexes (nombreux bâtiments).
• SST $k - \omega$ (Shear Stress Transport)
  • A priori la meilleure variante du $k - \epsilon$.

Globalement, tous les modèles linéaires souffrent de l’hypothèse d’isotropie des contraintes normales de Reynolds. Il faudrait utiliser en théorie des modèles non-linéaires, mais encore trop peu d’applications existent à ce jour.

En résumé, il vaut mieux éviter le modèle $k - \epsilon$ tel quel et utiliser au moins Realizable $k - \epsilon$ , ou mieux encore le SST $k - \omega$. Dans le cas de modélisation de multiples bâtiments, le modèle RNG $k - \epsilon$ semble aussi fournir de bons résultats.

La couche limite atmosphérique

Origines

Lorsqu’un fluide s’écoule le long d’une paroi immobile, les vitesses sur la paroi sont nulles, alors que loin de la paroi elles sont égales à la vitesse de l’écoulement non perturbé, comme illustré sur la figure 2. Cette zone de gradient de vitesse se nomme Couche Limite.

Selon une normale à la surface, la vitesse va donc varier entre $0$ [m.s^-1] et un maximum. La loi de variation dépend de la viscosité du fluide qui induit des frottements entre les couches voisines.

La couche limite se définit comme la région de l’écoulement où les effets visqueux sont au moins aussi importants que les effets inertiels. De façon générale on définit $\Delta (x)$ l’épaisseur de la couche limite telle que :

$$ u\left( \Delta (x) \right) = 0,99 \times u_{\text{inf}} $$

où $u_{\text{inf}}$ est la vitesse à l’infini (écoulement non perturbé), comme illustré dans la Figure 3 qui suit.

Définitions

La couche limite atmosphérique (CLA) représente la partie de l’atmosphère directement sensible aux effets de la surface terrestre (continentale ou océanique). Au-delà de cette zone, on considère que l’influence du sol est négligeable. On peut schématiquement diviser cette couche en 3 tranches, comme sur la Figure 3 (inspirée de [Grazzini, 1999]⁴).

La première sous-couche rugueuse au voisinage immédiat du sol est une zone de mélange des sillages des obstacles rencontrés par le vent. Les champs de vitesses y sont très perturbés, et les forces de frottement y sont prépondérantes. Son épaisseur peut varier de quelques millimètres (surface de la mer) jusqu’à une dizaine de mètres (zones urbaines).

Immédiatement au-dessus se situe la Couche Limite de Surface (CLS), aussi appelée couche de mélange ou couche de Prandtl. Elle se définit comme la région où la température diminue rapidement avec l’altitude en journée, et plus généralement comme la région où les flux de quantité de mouvement et de chaleur (sensible et latente) sont conservatifs et égaux à ceux du sol. Au sein de cette zone, la turbulence est homogène, et la direction du vent ne varie pas (ou peu) avec l’altitude. La vitesse, elle, varie proportionnellement au logarithme de l’altitude. Elle représente environ 10% de la CLA, et s’étend sur des hauteurs allant de la dizaine à la centaine de mètres.

Enfin, on trouve la couche d’Ekman (aussi appelé couche limite dynamique), où la force de Coriolis (induite par la rotation de la terre) devient comparable aux forces de frottement au sol. Cela entraine des déviations de la direction du vent avec l’altitude, pouvant atteindre 30 à 40° dans la direction anticyclonique (vers la droite dans l’hémisphère Nord, et vers la gauche dans l’hémisphère Sud).

Caractérisation

La représentation du profil moyen du vent est établie à partir de résultats issus de la mécanique des fluides. En pratique, deux types de profils sont utilisés, l’un suit une logarithmique et l’autre une loi puissance, dépendant de l’altitude z .

Profil logarithmique

Bien qu’il n’existe pas à ce jour de description exacte de cette CLA, de nombreuses approximations existent. La plus répandue, préconisée par différents auteurs (e.g. [Frank et al, 2004]³ ou [Blocken et al., 2012]⁵) est la description proposé par [Richards et al., 1993]⁶. Celle-ci fait l’hypothèse d’une vitesse verticale nulle, d’une pression constante dans la direction verticale et la direction d’écoulement, ainsi que de contraintes de cisaillement constantes dans la CLA. Ainsi le profil vertical de vitesse peut s’obtenir par la loi logarithmique suivante :

$$ U(z)=\frac{u_*}{K} \ln \left( \frac{z+z_0}{z_0}\right) $$

où $U(z)$ est la vitesse à une altitude $z$, $z_0$ est la rugosité de surface, $K$ la constante de von Kármán (usuellement la valeur est autour de 0,4), et enfin \(u_{*}\) la vitesse de frottement. Cette dernière s’obtient par la relation :

$$ u_* = \frac{KU_h}{\ln \left( \frac{h+z_0}{z_0} \right)} $$

Elle est calculée à partir de la vitesse de référence \(U_h\) mesurée à une altitude \(h\). C’est typiquement la valeur mesurée dans les stations météorologiques à 10 m de hauteur.

Profil loi puissance

Dans certain cas, par soucis de simplicité, une approximation plus simple peut être utilisée et permet une évaluation plus rapide des vitesses d’écoulements en fonction de la hauteur. Il s’agit de la “loi puissance”, qui s’exprime comme suit :

$$ \frac{u(z)}{u_0} = \left( \frac{z}{H_0} \right)^m $$

Où $u_0$ est la vitesse d’écoulement d’air mesurée à une hauteur \(H_0\), \(z\) est l’altitude et enfin \(m\) un paramètre qui caractérise la distorsion du profil. Sa valeur est usuellement comprise entre 0.1 et 0.2.

Note : Il n’y a pas de justification mathématique sur la validité de cette loi, cependant elle semble être en adéquation avec les profils expérimentaux mesurés, surtout après avoir ajusté le paramètre \(m\).

Turbulence de la CLA

Lorsque l’on utilise des modèles de turbulences à 2 équations (les plus utilisés), les valeurs et profils de \(k\) et \(\epsilon\) s’obtiennent aussi avec l’hypothèse de CLA. Les formules pour les déterminer se trouvent dans [Ricahards et al., 1993]⁶ et sont le résultat d’équations de conservations.

Le traitement des parois

Figure tirée de wikipedia : Law of The Wall

La couche limite en abord de parois est divisée en trois zones caractérisées par la distance à la paroi adimensionnelle \(y^+\) définie telle que \(y^+=\frac{U^+ y}{\nu}\) avec \(U^+\) la vitesse de cisaillement, \(\nu\), la viscosité cinématique, et \(y\) la distance à la paroi. Ces trois zones, illustrées en Figure 5, sont :

• la sous-couche visqueuse: dans cette zone où \(0 \lt y^+ \lt 5\), les effets visqueux dominent.
• la zone tampon où \(y^+\) est tel que \(5 \lt y^+ \lt 30\).
• la couche logarithmique avec \(y^+\) compris entre 30 et 300. Dans cette zone, les effets turbulents dominent et la vitesse suit la loi logarithmique définie par:

$$ U^+=\frac{1}{\kappa}\log y^+ + C $$

Il y a deux approches couramment utilisées en CFD pour modéliser la turbulence en paroi. D’un coté, il y a les modèles de turbulences à “faible nombre de Reynolds” qui calculent l’écoulement dans toute la couche limite. L’inconvénient de cette méthode est qu’elle requiert un maillage très fin avec des premières mailles localisées dans la sous-couche visqueuse. D’autre part, il y a les modèles de turbulence dit à “haut nombre de Reynolds” qui modélisent l’écoulement dans la couche limite via la loi logarithmique. La précision des résultats en abord de paroi est réduite mais ces méthodes diminuent considérablement les temps de calcul.

Les recommandations pour l’aéraulique urbaine

Zone de modélisation

Si l’on s’intéresse à un bâtiment de hauteur $H$, il convient de modéliser les faces latérales et la face supérieure à une distance minimum de $5H$. Idéalement sur la face aval, il convient de modéliser une distance \(15H\) [Frank et al., 2004]³ et [Tominaga et al., 2008]⁷. La distance accrue en aval de l’objet d’étude a pour but d’avoir “assez” d’espace pour obtenir un écoulement stable.

Figure d’après [Frank, 2007]⁸

D’une façon générale, le “blocage” (ratio de surfaces entre la projection de l’objet d’étude et la face amont) doit être inférieur à 10%, et idéalement autour de 3 à 5%.

Discrétisations

Recommandations générales

• D’une façon générale, il est préférable d’utiliser des éléments prismatiques (plutôt que tétraédriques).
• Sauf cas particuliers, afin de limiter les problèmes de diffusion numérique, le maillage doit être construit dans l’axe du vent.
• Pour tous les bâtiments modélisés, il convient d’avoir au moins 10 nœuds par dimensions.
• Pour les “rues” (espace séparant 2 bâtiments voisins), il convient d’avoir un minimum de 6 à 8 nœuds ( < 2m) afin de pouvoir représenter les jets d’air, ainsi que les possibles décollements de parois.
• Dans le cas où l’environnement urbain est aussi modélisé, les « autres » bâtiment doivent avoir au minimum 3 mailles dans toutes les directions
• La résolution verticale doit être fine ($\delta \le 10 m$) jusqu’à une hauteur de 3H.
• La résolution verticale doit être fine près du sol, spécialement si l’on s’intéresse au confort des piétons.

Recommandations spécifiques au bâtiment d’étude

• Sur le toit, et près des arrêtes (coins) du bâtiment, le maillage doit être “assez fin” pour être à même de représenter les décrochements de parois et les zones de recirculations (\(\delta < \frac{L_{\text{max}}}{10}\)) .
• Globalement on veut des maillages de plus en plus fins à l’approche du bâtiment, cependant le facteur d’expansion doit rester inférieur à 1,2 – 1,3 près des surfaces (ratio entre les tailles de 2 cellules adjacentes).