# Le modèle non-linéaire $k - \varepsilon - f_\mu$ et son application à l'étude
De l'écoulement turbulent autour de profil aérodynamique

Ahmed Ezzarfi* et Bousselham Harbouch†

Laboratoire d'Energétique, Mécanique des Fluides et Matériaux, UAE
Faculté des Sciences, B.P. 2121, Avenue Sebta, Mhannech II, Tétouan, Maroc

(reçu le 11 Juillet 2011 – accepté le 26 Décembre 2011)

**Résumé** - Le but de ce travail est d'évaluer l'aptitude du modèle de turbulence non-linéaire $k - \varepsilon - f_\mu$, mis au point par Park et al. (2003) [5], pour reproduire correctement l'anisotropie des contraintes normales du tenseur de Reynolds en proche – paroi. Les résultats obtenus sont comparés à ceux du modèle linéaire $k - \varepsilon - f_\mu$ et du modèle standard $k - \varepsilon$. L'application concerne l'étude de l'écoulement turbulent bidimensionnel autour d'un profil aérodynamique NACA 0012, soumis à un angle d'incidence de 0° et un nombre de Reynolds de 500.000. La méthode des volumes finis est utilisée pour résoudre les équations de conservation de la masse et de la quantité de mouvement. La prédiction de l'anisotropie est fortement liée aux corrections dans la région proche paroi. Le présent modèle est parfaitement réalisable et reproduit bien la relation non linéaire dans la région proche - paroi. A la lumière des résultats obtenus, on montre que le modèle linéaire $k - \varepsilon - f_\mu$ et le modèle standard $k - \varepsilon$ sont incapables à prédire correctement l'anisotropie des contraintes normales au voisinage de la paroi.

**Abstract** - The aim of this study is to evaluate the ability of non-linear model $k - \varepsilon - f_\mu$ recently developed by Park et al. (2003) [5], to reproduce correctly the anisotropy of normal stresses tensor Reynolds in the near-wall. The results are compared with the linear model $k - \varepsilon - f_\mu$ and standard model $k - \varepsilon$. The application concerns the study of two-dimensional turbulent flow around the airfoil NACA 0012, with an incidence of 0° and a Reynolds number of 500.000. In the present study, a numerical simulation, based on the finite - volume procedure is constructed and applied for solving the conservation equations of mass and momentum. The prediction of the anisotropy is strongly related to the corrections in the near-wall region. In light of the results, we show that the non-linear model reproduces correctly the anisotropy of normal stresses in the near-wall, while the linear model $k - \varepsilon - f_\mu$ and standard model $k - \varepsilon - f_\mu$ are unable to predict them in the near the wall region.

**Keywords:** Nonlinear model $k - \varepsilon - f_\mu$ - Airfoils NACA - Anisotropy - Aerodynamic.

## 1. INTRODUCTION

Le développement des modèles de fermeture turbulente dans le cadre d'une approche de (RANS) a fait l'objet de nombreuses études approfondies [6, 8]. L'évaluation de ces modèles présente beaucoup de difficultés dues à la complexité des hypothèses physiques faites pour assurer la fermeture et l'adaptation des conditions de réalisabilité. En général, les divers modèles de turbulence mis en œuvre sont appliqués, sans modifications, à des écoulements de natures différentes.

<footnote>* zarfi5@yahoo.fr Adresse: Hay Salam Groupe ‘O’, Rue 8, N°24, Ksar El Kebir, Maroc</footnote>

<footnote>† kbousselham@yahoo.fr</footnote>
Néanmoins, leur modélisation correcte permet d'assurer une bonne représentation d'un écoulement turbulent homogène ou non. Les quantités turbulentes en proximité de la paroi sont fortement anisotropes car les fluctuations normales de vitesse turbulente à la paroi sont amorties par la présence de celle-ci.

Les modèles linéaires de turbulence se fondent sur l'hypothèse de Boussinesq qui introduit une relation de colinéarité entre les tensions turbulentes et le tenseur des taux de déformation moyens grâce à un concept de viscosité de turbulence. Cette approche, très largement utilisée, peut conduire à des prédictions convenables dans des configurations simples.

Néanmoins cette loi constitutive présente un certain nombre de limitations. La dépendance linéaire des contraintes turbulentes vis-à-vis du tenseur de déformation ne permet pas de capturer efficacement l'anisotropie des tensions normales, y compris dans des configurations de référence telles que l'écoulement de couche limite turbulente (Bradshaw, 1973 [4]).

Cette relation linéaire peut conduire à une prédiction erronée des quantités turbulentes, en particulier dans les régions de brusques modifications du tenseur de déformation. Par exemple une surproduction d'énergie cinétique turbulente non physique peut être observée au niveau de point d'arrêt dans l'écoulement autour d'une aile, cette région correspondant à un lieu de déformation intense (Jin et al. 1994 [16]).

Afin de pallier certaines limitations des modèles linéaires, des lois constitutives non-linéaires ont été développées. Elles consistent à inclure, dans la relation constitutive, des tenseurs d'ordre supérieurs issus de combinaisons des tenseurs des taux de déformation et de rotation moyens (Pope, 1975 [21], par exemple). Plusieurs approches ont été envisagées par analogie aux lois constitutives considérées en rhéologie pour les fluides visco-élastiques ou à effet de mémoire donnant lieu aux Non Linear Eddy Viscosity Models (NLEVM) (Shih et al. 1993 [17], par exemple). Ces modèles sont ainsi plus sujets aux instabilités numériques.

De nombreuses versions de ces modèles ont été développées en tant que modèles d'extensions aux modèles linéaires isotropes à viscosité turbulente, on peut citer les travaux de Apsley et al. 1998 [11], Craft et al. 1996 [25], Gatski et al. 1993 [20], Girimaji, 1996 [22], Wallin et al. 2000 [23]. Ce type d'approche a notamment pour objectif une meilleure capture de l'anisotropie des contraintes normales du tenseur de Reynolds en proche-paroi par rapport aux modèles linéaires, fondés sur l'hypothèse de Boussinesq qui prédit un comportement isotrope de ces contraintes.

L'objectif de ce présent travail, est d'évaluer l'aptitude du modèle de turbulence non linéaire $k - \epsilon - f_\mu$, récemment mis au point par Park et al. (2003) [5], pour détecter l'anisotropie des contraintes normales du tenseur de Reynolds en proche – paroi.

L'application du modèle s'intéresse à l'écoulement turbulent bidimensionnel autour d'un profil aérodynamique NACA 0012, soumis à une incidence de $\theta$ et un nombre de Reynolds de 500.000. La méthode des volumes finis est utilisée pour résoudre les équations de conservation de la masse et de la quantité de mouvement. La formulation non linéaire du modèle est obtenue en se basant sur le théorème de Cayley – Hamilton (Pope, 1975) [21] pour un écoulement homogène.
La prédiction de l'anisotropie est fortement liée aux corrections dans la région proche paroi. Le présent modèle est parfaitement réalisable et reproduit correctement l'anisotropie en proximité de la paroi.

A la lumière des résultats obtenus, on montre que le modèle linéaire $k - \varepsilon - f_\mu$ et le modèle standard $k - \varepsilon$ restent incapables à prédire l'anisotropie des contraintes normales au voisinage de la paroi, ce qui n'est pas le cas pour le modèle non-linéaire $k - \varepsilon - f_\mu$ qui est bien adapté pour l'étude des écoulements turbulents autour des profils aérodynamiques et que sa performance est généralement très satisfaisante.

## 2. FORMULATION MATHEMATIQUE

### 2.1 Equations de conservation

Pour un écoulement turbulent incompressible, bidimensionnel, les équations de conservation de la masse et de la quantité de mouvement sont écrites dans le système de coordonnées cartésiennes $(x_1, x_2)$ avec des notations tensorielles sous la forme suivante:

$$ \frac{\partial U_i}{\partial x_i} = 0 \qquad (1) $$

$$ U_j \frac{\partial U_i}{\partial x_j} = - \frac{1}{\rho} \frac{\partial P}{\partial x_i} + \frac{\partial}{\partial x_j} \left[ v \frac{\partial U_i}{\partial x_j} - \overline{u_i u_j} \right] \qquad (2) $$

Où $U_j$, représente la vitesse moyenne, $u_j$, la fluctuation de vitesse. $P$, $\rho$ et $v$, sont respectivement la pression moyenne, la masse volumique et la viscosité moléculaire de fluide.

### 2.2 Modèles de turbulence

#### 2.2.1 Modèle non - linéaire $k - \varepsilon - f_\mu$

Selon la formulation de Park et al. [5], le tenseur de Reynolds $(\overline{u_i u_j})$ de l'équation (2) s'écrit sous la forme suivante:

$$ \overline{u_i u_j} = \overline{u_i} \overline{u_j}^{\text{linear}} + \tau_{ij}^{\text{N}} \qquad (3) $$

$$ \text{avec, } \overline{u_i u_j}^{\text{linear}} = -v_t \times \left( \frac{\partial U_i}{\partial x_j} + \frac{\partial U_j}{\partial x_i} \right) + \frac{2}{3} k $$

$\tau_{ij}^{\text{N}}$ est le terme non linéaire du tenseur de Reynolds défini par:

$$ W_{i1}^* \times W_{im}^* \times S_{mj}^* + S_{il}^* \times W_{im}^* \qquad (4) $$

où,
$$S_{ij}^* = \frac{S_{ij}^k}{\epsilon}; W_{ij}^* = \frac{W_{ij}^k}{\epsilon}; S^* = \sqrt{2S_{ij}^* \times S_{ij}^*};$$

$$W^* = \sqrt{2W_{ij}^* \times W_{ij}^*} \text{ et III}_s = S_{im}^* \times W_{mn}^* \times W_{nl}^*$$

$S^*$ et $W^*$ représentent respectivement les normes des tenseurs de déformation et de rotation. Les coefficients $\beta_i$ selon Park et al. [5], sont donnés par:

$$\begin{aligned}
\beta_1 &= \alpha_1 / (1 - 2/3\alpha_2^2 \times S_{mn}^* \times S_{mn}^* + 2\alpha_3^2 \times W_{mn}^* \times W_{mn}^*); \quad \beta_2 = 2\alpha_2 \times \beta_1; \\
\beta_3 &= -\alpha_3 \times \beta_1; \quad \beta_4 = \alpha_2 \times \beta_3 - \alpha_3 \times \beta_2; \quad \beta_5 = -\alpha_3 \times \beta_2
\end{aligned} \quad (5)$$

Avec, $\alpha_1$, $\alpha_2$ et $\alpha_3$ sont réarrangés à partir de l'évaluation de Gatski et al. [20] sous la forme suivante:

$$\alpha_1 = \tilde{\alpha}_1 / g = -0.48; \quad \alpha_2 = \tilde{\alpha}_2 / g = -0.375; \quad \alpha_3 = \tilde{\alpha}_3 / g = -0.8$$

$g$ est une variable positive pour toutes les valeurs $S^*$ et $W^*$, elle dépend essentiellement des conditions de l'écoulement. La viscosité turbulente $v_t$ définit par:

$$v_t = C_\mu \times f_\mu \times \frac{k^2}{\epsilon} \quad (6)$$

$f_\mu$ est une fonction correctrice tenant compte des effets de la paroi. Elle se décompose comme produit de deux fonctions, c'est-à-dire $f_\mu = f_{\mu 1} \times f_{\mu 2}$, où $f_{\mu 1}$ tient compte des effets approximatifs au voisinage de la paroi, et elle est obtenue en faisant la résolution de l'équation elliptique de Helmotz en $f_\omega$ sous forme:

$$f_{\mu 1} = \left( 1 + f_D \times R_t^{-3/4} \right) \times f_\omega^2 \quad (7)$$

$$\frac{\partial^2 f_{\omega}}{\partial x_j \partial x_j} = \frac{(f_{\omega} - 1) \times R_t^{3/2}}{A^2 L^2} \quad (8)$$

La fonction $f_{\mu 2}$ représente les effets de déformation [12], donnée par:

$$f_{\mu 2} = \frac{\frac{15}{3}(1+g)}{(g^2 + C_{\mu} \times g^3 + A_s)} \quad (9)$$

$$g = \begin{cases} \frac{C_0}{3} + (P_1 + \sqrt{P_2})^{1/3} + \text{Sign}(P_1 - \sqrt{P_2}) \times |P_1 - \sqrt{P_2}|^{1/3} & P_2 \ge 0 \\ \frac{C_0}{3} + 2(P_1^2 + P_2^2)^{1/6} \cos\left(\frac{1}{3}\arccos\left(\frac{P_1}{\sqrt{P_1^2 - P_2^2}}\right)\right) & P_2 < 0 \end{cases} \quad (10)$$
La fonction $F_D$ est définie par $F_D = 10 \exp(-(R_t / 120)^2)$, L représente l'échelle de longueur de turbulence (Durbin *et al.*, 1996 [7]), qui est défini par

$$L^2 = \frac{k^3}{\varepsilon^2} + 70^2 \sqrt{\frac{v^2}{\varepsilon}}.$$

Les constantes du modèle sont données par: $C_\mu = 0.09$, $A = 8.4$ et $C_0 = 2.5$.

$P_1$ et $P_2$ sont données par:

$$P_1 = C_0 \left[ \frac{C_0^2}{27} - \frac{(A_s + \tilde{\alpha}_1 \eta^2 + \frac{1}{2})}{6} \right] \text{ et } P_2 = P_1^2 - \left[ \frac{C_0^2}{9} - \frac{(A_s + \tilde{\alpha}_1 \eta^2)}{3} \right] \quad (11)$$

$$\text{où, } \quad A_s = \tilde{\alpha}_3^2 \times \xi^2 - \frac{\tilde{\alpha}_2^2 \eta^2}{3}$$

Afin de tenir compte des effets approximatifs au voisinage de la paroi, les taux moyens de déformation et rotation sont modifiés comme suit: $\eta = f_\omega S^*$ et $\xi = f_\omega W^*$, où $f_\omega$ est utilisée comme fonction correctrice en proximité de la paroi.

Les équations de transport de l'énergie cinétique turbulente $k$ et son taux de dissipation $\varepsilon$ sont écrites dans le système de coordonnées cartésiennes sous la forme suivante:

$$U_j \frac{\partial k}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( v + \frac{v_t}{\sigma_s} \right) \frac{\partial k}{\partial x_j} \right] + P_k - \varepsilon \quad (12)$$

$$U_j \frac{\partial \varepsilon}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( v + \frac{v_t}{\sigma_s} \right) \frac{\partial \varepsilon}{\partial x_j} \right] + \frac{(C_{\varepsilon_1}^* P_k - C_{\varepsilon_2} f_2 \varepsilon)}{T} + C_{s2} (1 - f_\omega) v v_t \left( \frac{\partial^2 U_i}{\partial x_j \partial x_k} \right)^2 \quad (13)$$

Les constantes du modèle sont données par: $\sigma_k = 1.1$, $\sigma_{\varepsilon_1} = 1.3$, $\sigma_{\varepsilon_2} = 1.9$, $\sigma_{\varepsilon_3} = 0.8$, le terme production de l'énergie cinétique turbulente $P_k$ est défini par: $P_k = -\overline{u_i u_j} \times \frac{\partial U_i}{\partial x_j}$, la fonction $f_2$ du modèle, qui décrit l'effet de dégradation de la turbulence (Coleman *et al.* 1992 [19]), est définie par: $f_2 = 1 - \frac{1}{9} \exp(-0.33 \times R_t^{1/2})$.

L'échelle de temps turbulent, $T$ (Durbin *et al.* [7]) est défini par:

$$T = \sqrt{\left(\frac{k}{\varepsilon}\right)^2 + 36\left(\frac{v}{\varepsilon}\right)}, \text{ le coefficient } C_{\varepsilon_1}^* \text{ est employé comme valeur constante dans de nombreux modèles de type } k - \varepsilon.$$

Afin de mettre en exergue les effets de déformation, ce coefficient est modélisé selon Park *et al.* 1997, [18] de la manière suivante:
$$C_{\varepsilon_1}^* = 1.42 + \frac{C_\mu}{(1 + 5f_{\mu_2} \eta^2)} \quad (14)$$

## Modification des coefficients du modèle

Afin de détecter l'anisotropie des contraintes normales en proximité de la paroi, les coefficients $\beta_i$ de l'équation (5) doivent être modifiés selon la formulation de Park et al. [5], sous forme:

$$\beta_2 = \tilde{\beta}_2 + \tilde{\beta}_{2,wall} \times C_{\omega}; \beta_3 = \tilde{\beta}_3 + \tilde{\beta}_{3,wall} \times C_{\omega}$$

$$\beta_4 = \tilde{\beta}_4; \beta_5 = \tilde{\beta}_5$$

La constante $C_\omega$ du modèle est donnée par:

$$C_{\omega} = 1 \text{ pour } i=j \text{ et } C_{\omega} = 0 \text{ pour } i \neq j$$

Les coefficients $\tilde{\beta}_i$, qui sont liés aux taux de déformations et rotations, sont données par:

$$\tilde{\beta}_1 = C_\mu \times f_{s_1} \times f_\omega^2; \tilde{\beta}_3 = C_\mu \times f_{s_2} \times f_\omega^2$$

$$\tilde{\beta}_{2,wall} = (1 - f_{\omega}) \times 2.5 \frac{S_B}{S_{\omega}}; \tilde{\beta}_{3,wall} = (1 - f_{\omega}) \times \left( 1.8 \times \frac{S_B}{S_{\omega}} - \tilde{\beta}_3 \right)$$

$$\tilde{\beta}_4 = -C_\mu (f_{s_1} \times f_\omega)^2; \tilde{\beta}_5 = -C_\mu (f_{s_2} \times f_\omega)^2$$

$$\text{Où, } f_{s_1} = \frac{(1 + \eta_s f_{\mu_2})}{(1 + 4\eta_s + 4\eta_s^2)}; \quad f_{s_2} = \frac{(1 + \eta_s f_{\mu_2})}{(1 + \eta_s + 8\eta_s^2)}$$

$$\eta_s = C_\mu \text{Max}(\eta, \xi); S_\omega = 1 + [\text{Max}(S^*, W^*)]^2; S_B = \frac{2S_{mn} S_{mn} (1 - \delta_{mn})}{S^{*2}}$$

La performance du modèle non-linéaire sera comparée aux résultats des deux modèles de turbulences, qui sont explicités ci-dessous.

### 2.2.2 Modèle linéaire $k - \varepsilon - f_\mu$

Les constantes et les fonctions de ce modèle sont résumées comme suit:

$$C_{\mu} = 0.09; C_{\varepsilon_1}^* = 1.42 + \frac{C_{\mu}}{(1 + 5f_{\mu_2} \eta^2)}; \sigma_k = 1.1; \sigma_{\varepsilon} = 1.3; \sigma_{\varepsilon_2} = 1.9;$$

$$C_{\varepsilon_3} = 0,8; f_{\mu_1} = (1 + f_D \times R_t^{-3/4}) \times f_\omega^2; f_{\mu_2} = \frac{\frac{15}{3}(1 + g)}{(g^2 + C_\mu \times g^3 + A_s)}$$

$$f_2 = 1 - \frac{1}{9} \exp(-0.33 \times R_t^{1/2}); T = \sqrt{\left(\frac{k}{\varepsilon}\right)^2 + 36\left(\frac{v}{\varepsilon}\right)} \quad (15)$$
### 2.2.3 Modèle standard $k - \varepsilon$

Les constantes et les fonctions de ce modèle sont résumées au **Tableau** suivant:

<table><thead><tr><th>C<sub>μ</sub></th><th>C<sub>ε2</sub></th><th>C<sup>*</sup><sub>ε1</sub></th><th>C<sub>ε3</sub></th><th>σ<sub>k</sub></th><th>σ<sub>ε</sub></th><th>f<sub>μ</sub></th><th>f<sub>2</sub></th><th>T</th></tr></thead><tbody><tr><td>0.09</td><td>1.92</td><td>1.44</td><td>0.0</td><td>1.0</td><td>1.3</td><td>1</td><td>1</td><td>k/ε</td></tr></tbody></table>

## 3. ETUDE NUMERIQUE

Dans le système de coordonnées cartésiennes $(x_1, x_2)$, les équations de transport (1), (2), (10) et (11) qui gouvernent l'écoulement s'écrivent sous la forme générale suivante:

$$ \frac{\partial}{\partial x_j} \left( \rho \times U_j \times \phi - \Gamma^{\phi} \right) = S_{\phi} \qquad (16) $$

Où $\phi = [1, U_1, U_2, k, \varepsilon]$ et $\Gamma^{\phi} = \left[ 0, v_e, v_e, v + \frac{v_t}{\sigma_k}, v + \frac{v_t}{\sigma_{\varepsilon}} \right]$ sont respectivement variables de transport et coefficients de diffusion, $S_{\phi} = [S_1, S_{U_1}, S_{U_2}, S_k, S_{\varepsilon}]$ sont des termes sources donnés par:

$$ S_1 = 0 ; \quad v_e = v + v_t $$

$$ S_{U_1} = -\frac{\partial P}{\partial x_1} + \frac{\partial}{\partial x_1} \left( v_e \frac{\partial U_1}{\partial x_1} \right) + \frac{\partial}{\partial x_2} \left( v_e \frac{\partial U_2}{\partial x_1} \right) - \frac{2}{3} \frac{\partial k}{\partial x_1} - \frac{\partial \tau_{11}^N}{\partial x_1} - \frac{\partial \tau_{12}^N}{\partial x_2} $$

$$ S_{U_2} = -\frac{\partial P}{\partial x_2} + \frac{\partial}{\partial x_2} \left( v_e \frac{\partial U_2}{\partial x_2} \right) + \frac{\partial}{\partial x_1} \left( v_e \frac{\partial U_1}{\partial x_2} \right) - \frac{2}{3} \frac{\partial k}{\partial x_2} - \frac{\partial \tau_{21}^N}{\partial x_1} - \frac{\partial \tau_{22}^N}{\partial x_2} $$

$$ S_{\varepsilon} = \frac{(C_{\varepsilon1}^* \times P_k - C_{\varepsilon2} \times f_2 \times \varepsilon)}{T} + C_{\varepsilon3} (1 - f_{\omega}) v v_t \left( \frac{\partial^2 U_i}{\partial x_j \partial x_k} \right)^2 $$

$$ S_k = P_k - \varepsilon $$

En utilisant les règles de transformation (Thompson et al. 1985) [24], l'équation de conservation générale (13) sera écrite, après transformation, dans le système de coordonnées curvilignes $(\xi_1, \xi_2)$ comme suit:

$$ \frac{\partial}{\partial \xi_i} \left( G_i \times \phi - \Gamma^{\phi} \times g^{ij} \times \frac{\partial \phi}{\partial \xi_j} \right) = S_{\phi} \times J \qquad (17) $$

Les paramètres suivants: $G_i$, la vitesse contravariante, $g^{ij}$, le tenseur métrique contravariant et $J$, le Jacobien de transformation, sont définis par les expressions ci-après:
$$G_i = U_j \times \frac{\partial \xi_i}{\partial x_j} \times J ; \quad g^{ij} = \frac{\partial \xi_i}{\partial x_k} \times \frac{\partial \xi_j}{\partial x_k} ; \quad J = \frac{\partial \xi_1}{\partial x_1} \times \frac{\partial \xi_2}{\partial x_2} - \frac{\partial \xi_2}{\partial x_1} \times \frac{\partial \xi_1}{\partial x_2} \quad (18)$$

## 3.1 Discrétisation spatiale

La méthode numérique utilisée, pour intégrer les équations de la forme générale (15) en variable $\phi$, est celle de volumes finis [10]. Ces équations sont intégrées sur un volume de contrôle, sachant que le maillage adopté pour un système de coordonnées curvilignes généralisées, est de type collocatif (non-staggered grid) dont lequel toutes les variables $\phi$ sont calculées aux nœuds principaux des mailles non décalées. Le couplage pression-vitesse a été traité en utilisant la méthode de substitution de la pression 'PSM' [15].

Le principe de cette méthode consiste à substituer l'équation discrétisée de quantité de mouvement dans l'équation discrétisée de continuité en faisant établir l'équation implicite de la pression sans faire appel à aucune méthode de correction.

Les termes d'advection sont évalués à l'aide de schéma quadratique 'QUIK' de Leonard (1979) [9], alors que les termes de diffusion sont approximés par un schéma centré aux différences finies de second ordre. Le système constitué par les équations finales discrétisées écrites en chaque nœud du domaine, ainsi que les conditions aux limites, est résolu itérativement grâce à l'application de l'algorithme de Thomas 'TDMA'.

Du fait du couplage et de la non-linéarité des équations à résoudre, la stabilité du processus itératif ne peut être assurée que par l'utilisation de coefficients de sous-relaxation et la linéarisation des termes sources.

Le maillage adopté dans le cadre de cette étude est de type H qui est illustré dans la figure 1. Il est constitué par les lignes équipotentielles et les lignes de courant de fonctions respectivement $\xi_1(x_i)$ et $\xi_2(x_i)$ correspondantes à un écoulement potentiel de fluide parfait autour de même profil.

<figure><img src="image_1.png" /><figcaption>Fig. 1: Maillage autour de profil NACA 0012</figcaption></figure>

L'étude de cet écoulement à pour but de faciliter la résolution numérique des équations gouvernant l'écoulement réel en faisant calculer les fonctions $\xi_1(x_i)$ et $\xi_2(x_i)$ en tout point du domaine de calcul. Ces deux fonctions, qui sont utilisées
ultérieurement comme coordonnées d'espace, seront calculées par la méthode de singularités de types (sources - tourbillon) [13, 14].

L'avantage d'adopter ce type de maillage, d'une part la simplicité de calculer ces deux fonctions dans tout le domaine, d'autre part la condition d'orthogonalité de ces deux fonctions est vérifiée en tout point du domaine étudié.

Dans ce travail, le nombre total de mailles est de l'ordre de 250421 dont 168 nœuds disposés sur la paroi du profil. La distribution de nœuds est non uniforme, ils sont serrés dans les zones à fort gradient (i.e. bord d'attaque et bord de fuite), cette distribution est assuré par l'application de la technique d'interpolation par spline cubique. La taille de la première maille en proche paroi est donnée en échelles d'unité de paroi comme: $y^+ = 2$ et $x^+ = 25$. La dimension de mailles considérées, est généralement, suffisante pour détecter les grandeurs physiques de l'écoulement dans la couche limite.

Le maillage réalisé dans le cadre de cette étude est caractérisé par le fait qu'il est lié au profil, et l'incidence de l'écoulement étudié est introduite par l'intermédiaire des conditions aux limites.

## 3.2 Les conditions aux limites

Vu la nature elliptique des équations, les conditions aux limites sont nécessaires sur toutes les frontières du domaine de calcul. Les conditions appliquées sont données par:

-En amont: $U = U_\infty \cos \alpha$; $V = U_\infty \sin \alpha$; $P = P_\infty$; $f_\omega = 1$;

$$k_\infty = 1m^2/S^2; \quad \varepsilon_\infty = (C_\mu k_\infty^{3/2})/(2.51)$$

-Sur la paroi: $U = V = 0$; $\partial P / \partial n = 0$; $k = 0$; $\partial \varepsilon / \partial n = 0$; $f_\omega = 0$

avec, n: normale à la paroi.

-En aval: Le gradient de toutes les variables transportables $\phi$ est nul.

-Sur les deux frontières externes nord et sud, les flux normaux à ces surfaces sont nuls.

# 4. RESULTATS ET DISCUSSIONS

## 4.1 Performance de modèle non linéaire $k - \varepsilon - f_\mu$

Le comportement des coefficients $\tilde{\beta}_2$ et $\tilde{\beta}_3$, en fonction du taux de déformation adimensionné $S^*$, pour l'écoulement cisaillé autour du profil étudié, est illustré dans la figure 2. La comparaison des résultats est en bonne accord avec les travaux de Park et al. [5] pour l'étude de l'écoulement à l'intérieur d'un canal.

## 4.2 Grandeurs moyennes de l'écoulement

Le coefficient de pression surfacique moyen $C_p = \frac{(P - P_\infty)}{(0.5 \rho U_\infty^2)}$ autour du profil est un paramètre primordial, car il détermine l'amplitude du gradient de pression autour du profil, et est donc relié directement au comportement de la couche limite et en particulier à l'emplacement de la zone de transition. La figure 3 montre l'allure ce coefficient de pression et le compare à des valeurs mesurées par Lee et al. [2].
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Comportement des coefficients $\tilde{\beta}_2$ et $\tilde{\beta}_3$ pour l'écoulement cisaillé simple</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3: Distribution de coefficient de pression $C_p$</figcaption></figure>

On constate un excellent accord entre les valeurs expérimentales et les résultats de la simulation obtenus par les trois modèles non-linéaire, linéaire et le modèle standard $k - \varepsilon$. D'après l'allure de la courbe, on remarque qu'à partir de la position $x / c = 0.15$ vers le bord de fuite, la couche limite est soumise à un gradient de pression défavorable.

Le coefficient moyen de trainée vaut $C_d = 0.00921$ dans cette simulation, ce qui est légèrement plus petit que la valeur de 0.012174 simulée à l'aide du code de calcul CFL3D [27] de type RANS.

Les figures 4, 5 et 6 montrent respectivement $U(x,y)/U_\infty$, $V(x,y)/U_\infty$ et $C_p(x,y)$ en deux sections verticales $x/c = 0.0$ et $x/c = 0.5$. Les résultats obtenus par le modèle non-linéaire $k - \varepsilon - f_\mu$, sont en bon accord avec ceux calculés à l'aide du code de calcul CFL3D [27] de type RANS.
<figure><img src="image_5.png" /><figcaption>Fig. 4: Vitesse moyenne longitudinale, selon la verticale en $x / c = 0.0$ et $x / c = 0.5$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5: Vitesse moyenne transversale, selon la verticale en $x / c = 0.0$ et $x / c = 0.5$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6: Coefficient de pression $C_p$ selon la verticale en $x / c = 0.0$ et $x / c = 0.5$</figcaption></figure>
Pour le champ de pression, on remarque que le pic de surpression calculé par le présent modèle au bord d'attaque (x/c = 0.0) s'approche légèrement à la valeur maximale (Cp = 1.0) calculée par CFL3D, alors que le champ de vitesse au bord d'attaque est nul, c'est le point de stagnation. La symétrie de champ de vitesse longitudinale et celui de pression par rapport à la corde du profil se voit clairement, alors que le champ de vitesse transversale est antisymétrique.

<figure><img src="image_8.png" /><figcaption>Fig. 7: Profils de vitesse u/U∞ selon les quatre stations en extrados de NACA 0012</figcaption></figure>

Le comportement de la couche limite le long de l'extrados du profil a été illustré dans les figures 7a-b-c-d-. Les profils des vitesses moyennes normalisées par la vitesse infinie amont, sont tracés le long de l'extrados suivant la normale à la paroi.

Pour la position x/c = 0.52, on constate que la vitesse varie linéairement en fonction de la distance normale à la paroi, par conséquent le comportement de la couche
limite est laminaire. Une légère différence apparaît entre le profil de vitesse calculé par le modèle non-linéaire et celui du modèle LES d'Olivier et al. [1, 3].

Cette différence est due au fait que pour cette position, la couche limite laminaire commence à transitionner vers l'état turbulent sachant que la transition est très sensible pour les phénomènes de l'écoulement qui sont fortement liés aux conditions expérimentales, en particulier le niveau de l'intensité de la turbulence en amont [26].

De même dans notre calcul, l'écoulement est supposé complètement turbulent. Pour les trois autres positions, un bon accord entre les résultats du modèle non-linéaire et ceux du modèle LES [1, 3], ceci explique bien que pour ces positions la couche limite est turbulente.

Les profils de vitesse moyenne en variables de paroi, définis par: $U^+ = U / U_\tau$, $y^+ = y \times U_\tau / v$ sont illustrés en figure 8, où $U_\tau$, est la vitesse de frottement, y, la distance à la paroi et v, la viscosité cinématique.

En $x/c = 0.52$, le profil de vitesse calculé par le modèle non linéaire, est en bon accord avec la loi linéaire $u^+ = y^+$ indiquant que la vitesse augmente linéairement avec la distance par rapport à la paroi. C'est une caractéristique typique d'une couche limite laminaire, où le comportement du fluide est dominé par les effets de la viscosité.

Dans la zone de transition (après la position $x/c = 0.52$), on remarque que tous les profils de vitesse moyenne augmentent dans un premier temps linéairement jusqu'à $y^+ = 10$, après ils décrochent et augmentent selon une loi quasi- logarithmique.

Cependant, la pente de la croissance dans la zone de transition est supérieure à la pente de la loi logarithmique typique du comportement d'une couche limite turbulente
$$u^+ = 2.44 \times \log(y^+) + 5.2$$
formulée par Prandtl et Von Karman, qui caractérise le comportement de la couche limite turbulente.

Enfin, proche du bord de fuite en $x/c = 0.85$, on constate que le profil de vitesse calculé par le modèle non linéaire suit la loi logarithmique avec la bonne pente. Cela permet ainsi de situer avec une précision raisonnable, les bornes de la transition qui sont sensiblement au voisinage de $\approx [0.52-0.85]$, comme il a été observé expérimentalement par les travaux de Lee et al. [2].

La figure 9 représente la variation de vitesse moyenne en fonction des paramètres unités de paroi pour la position $x/c = 0.63$ en extrados de profil. Le profil de vitesse moyenne obtenu par le modèle non-linéaire est en bon accord avec les données expérimentales.

Cependant pour $y^+ \ge 30$, on constate que les profils de vitesse obtenus respectivement par le modèle linéaire $k - \varepsilon - f_\mu$ et le modèle standard $k - \varepsilon$, ne sont pas conformes aux résultats expérimentaux.

Une meilleure indication de la zone de transition peut être obtenue à partir de la variation de coefficient de frottement local $C_f = \tau_\omega / (0.5\rho \times U_\infty^2)$ illustrée sur la figure 10.
<figure><img src="image_9.png" /><figcaption>Fig. 8: Profils de vitesse en variable de paroi en x / c = 0.65 de l'extrados de profil</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Profils de vitesse en variable de paroi en x / c = 0.65 de l'extrados de profil</figcaption></figure>

En fait, la variation de ce coefficient dans la zone laminaire suit une loi différente que celle dans la zone turbulente. On constate que le frottement décroît presque à partir du bord d'attaque jusqu'à un minimum en x / c = 0.51, puis commence à augmenter jusqu'à x / c = 0.85 avant de décroître vers le bord de fuite.

Cela indique les limites approchées de la zone de transition ≅ [0.52 – 0.85] qui a été déjà mentionnée. La valeur minimale de frottement calculée est C<sub>f</sub> ≈ 6.5 × 10<sup>-4</sup> , elle est légèrement inférieure à la valeur C<sub>f</sub> ≈ 7.5 × 10<sup>-4</sup> mesurée par Lee *et al*. [2].

<figure><img src="image_11.png" /><figcaption>Fig. 10: Coefficient de frottement C<sub>f</sub> le Long de la corde de profil</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11: Profil de déficit de vitesse moyenne à partir du bord de fuite</figcaption></figure>
La figure 11 représente le profil de déficit de la vitesse moyenne adimensionnée,
$$ \frac{(U_\infty - U)}{U_\infty} $$
, le long de l'axe de sillage à partir du bord de fuite. Ce déficit, qui représente la trace de traînée de frottement, diminue au fur est à mesure en s'éloignant du bord de fuite.

## 4.3 Grandeurs turbulentes de l'écoulement

On peut également examiner des grandeurs liées à la turbulence, qui se développe dans la couche limite après le début de la transition. Les profils de fluctuations tracés selon la direction normale à la paroi peuvent également nous renseigner sur l'état de la couche limite. La figure 12 représente les fluctuations de vitesse longitudinale $u'_{rms} = u'_{rms} / U_\tau$ normalisée par la vitesse de frottement $U_\tau$, répartis à $x/c = 0.65$ en surface extrados.

La comparaison des résultats montre que les fluctuations de vitesse calculées par le modèle non linéaire sont en bonne concordance avec les données expérimentales de Lee *et al.* [2], ce qui n'est pas le cas pour le modèle linéaire et le modèle standard $k - \varepsilon$ dans la zone proche paroi ($y^+ < 100$). La valeur maximale de l'intensité turbulente $u'_{rms}$ obtenue par le modèle non-linéaire est sous-estimée d'environ 0.7 % de la valeur expérimentale mesurée par Lee *et al.* [2]. En s'éloignant de la paroi, ($y^+ > 150$), les résultats obtenus par les trois modèles de turbulence sont en bonne concordance avec les données expérimentales de Lee *et al.* [2].

<figure><img src="image_13.png" /><figcaption>Fig. 12: Fluctuations de vitesse longitudinale adimensionnée $u'_{rms} / U_\tau$ en $x/c = 0.65$</figcaption></figure>

La figure 13 représente les fluctuations de vitesse longitudinale $u'_{rms} / U_\infty$ dans la couche limite supérieure de profil. Elles sont tracées en pourcentage de la vitesse infinie $U_\infty$. Ces fluctuations sont calculées à une hauteur de $y^+ = 16$ au dessus de la paroi.

Les résultats obtenus par le modèle non-linéaire sont en bonnes concordances avec les calculs de Marsden *et al.* [1, 3] qui utilisent l'approche 'LES'. On observe que les
fluctuations sont quasi-nulles en amont du point $x/c = 0.65$, et croissent à partir de ce point. Il s'agit donc du début de la zone de transition.

La valeur maximale de l'intensité turbulente $u'_{rms}/U_\infty$ obtenue par le modèle non-linéaire est sous-estimée d'environ 6% de la valeur calculée par Marsden *et al.* [1, 3].

La différence est due à l'utilisation du maillage grossier composé de 250421 mailles, alors que Marsden a, en fait, utilisé un maillage plus fin composé de 12.3-millions mailles. En conséquence, le présent modèle permet de donner des résultats satisfaisants même avec un maillage n'est pas très fin.

<figure><img src="image_14.png" /><figcaption>Fig. 13: Fluctuations de vitesse $u'_{rms}/U_\infty$ le long de l'extrados en $y^+ = 16$</figcaption></figure>

La figure 14 montre le comportement des contraintes de Reynolds normales $\overline{u'u'}^+ = \overline{u'u'} / U_\tau^2$ et $\overline{v'v'}^+ = \overline{v'v'} / U_\tau^2$ et contrainte de cisaillement $-\overline{u'v'}^+ = -\overline{u'v'} / U_\tau^2$, normalisées par la vitesse de frottement, calculées à $x/c = 0.65$ en extrados de profil.

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 14: Profils des contraintes normales et de cisaillement à x/c=0.65 en extrados de NACA 0012</figcaption></figure>

Pour la contrainte longitudinale de la figure 15a-, les résultats obtenus par le modèle non-linéaire sont en bonne concordance avec les mesures expérimentales de Lee *et al.* [2]. La valeur maximale du tenseur de Reynolds longitudinal $\overline{u'u^+}$, obtenue par le modèle non-linéaire est sous-estimée d'environ 1.4 % de la valeur mesurée par Lee *et al.* [2].

En conséquence le modèle non-linéaire permet de détecter correctement l'anisotropie au voisinage de la paroi, tandis que le modèle linéaire et le modèle $k - \varepsilon$ standard sont incapables à prédire l'anisotropie de l'écoulement en proche paroi.

<figure><img src="image_17.png" /><figcaption>Fig. 15: Profils de l'énergie cinétique turbulente à x/c=0.65 en extrados de NACA 0012</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16: Variation de taux de dissipation à x/c=0.65 en extrados de NACA 0012</figcaption></figure>

On constate que les profils des contraintes normales et de cisaillement illustrés dans la figure 14 présentent des maximums presque en même distance de la paroi $y^+ = 25$.
On remarque que les maximums calculés par le modèle non linéaire $k - \varepsilon - f_\mu$ sont supérieurs à ceux calculés par le modèle linéaire et le modèle standard $k - \varepsilon$, et ce pour chaque composante du tenseur de Reynolds. En s'éloignant de la paroi ($y^+ \ge 125$), on constate que les trois modèles de turbulence donnent les mêmes résultats pour chaque composante du tenseur de Reynolds.

La figure 15 représente la variation d'énergie cinétique turbulente normalisée par la vitesse de frottement ($k^+ = k / U_\tau^2$) à ($x/c = 0.65$) en extrados de NACA 0012. On constate que la valeur maximale d'énergie cinétique prévue par le modèle non linéaire, en une distance $y^+ \approx 25$ de la paroi, est supérieure à celle calculée par le modèle linéaire et le modèle standard $k - \varepsilon$.

En conséquence, les fluctuations normales de vitesse en proche paroi, prévues par le modèle non-linéaire sont importantes par rapport aux modèles linéaires de turbulence.

En s'éloignant de la paroi ($y^+ \ge 125$), on remarque que l'énergie cinétique prévue par les trois modèles de turbulence est presque constante. Par conséquent, on peut conclure qu'à partir de cette distance, l'effet de la paroi est pratiquement négligeable.

La figure 16 représente variation du taux de dissipation de l'énergie cinétique turbulente $\varepsilon^+$, normalisée par la vitesse à l'infini amont $U_\infty$, ($\varepsilon^+ = \varepsilon \times C / U_\infty^2$), en fonction de paramètre $y^+$.

On remarque qu'au voisinage de la paroi ($y^+ \le 100$), le taux de dissipation prévu par le modèle non-linéaire est supérieur à celui calculé par les deux autres modèles de turbulence. En s'éloignant de la paroi ($y^+ > 100$), on constate que le taux de dissipation reste constant en s'approchant de la valeur nulle et ce pour les trois modèles de turbulence.

## 5. CONCLUSION

Le modèle de turbulence non linéaire $k - \varepsilon - f_\mu$ développé par Park et al. (2003) [5], a été appliqué pour étudier l'écoulement turbulent bidimensionnel autour de profil aérodynamique NACA 0012 soumis à une incidence de 0° et un nombre de Reynolds de 500.000.

Les contraintes de réalisabilités pour ce modèle sont bien vérifiées. Selon les résultats obtenus, on remarque que la prédiction de l'anisotropie des contraintes normales de tenseur de Reynolds au voisinage de la paroi, est fortement liée aux effets de correction des coefficients de modèle dans la région proche de la paroi.

La localisation de la zone de transition laminaire-turbulente se produit presque à moitié de la corde du profil, de même les grandeurs moyennes et turbulentes de l'écoulement prévues par ce modèle sont en bonne concordance avec les résultats expérimentaux disponibles.

La validité du modèle non linéaire, a été testée avec succès, en conséquence ce présent modèle est parfaitement réalisable et reproduit correctement l'anisotropie des
contraintes normales du tenseur de Reynolds au voisinage de la paroi, ce qui n'est pas le cas pour le modèle linéaire $k - \varepsilon - f_\mu$ et le modèle standard $k - \varepsilon$.

Des travaux seront réalisés, dans le futur, en appliquant ce modèle pour étudier les caractéristiques de l'écoulement turbulent autour des profils aérodynamiques aux incidences plus élevées.

## REFERENCES

[1] O. Marsden, C. Bogey, and C. Bailly, ‘Direct Noise Computation Around a 3-D NACA0012 Airfoil’, 2nd AIAA/CEAS Aeroacoustics Conference (27th AIAA Aeroacoustics conference), Cambridge, Massachusetts, 8-10 May 2006.

[2] H.K. Lee and S.H Kang, ‘Flow Characteristics of Transitional Boundary Layer on the NACA0012 Airfoil in Wakes’, Proceedings of FEDSM’98, Washington, DC, 21-25 June, 1998.

[3] O. Marsden, C. Bogey and C. Bailly, ‘Direct Noise Computation of the Turbulent Flow around a Zero-Incidence Airfoil’, AIAA Journal, Vol. 46, N°4, pp. 874–883, 2008.

[4] P. Bradshaw, ‘Effect of Streamline Curvature on Turbulent Flow’, Technical Report, Agardograph, N°169, 1973.

[5] T.S. Park, H.J. Sung and K. Suzuki, ‘Development of Nonlinear Near-Wall Turbulence Model for Turbulent Flow and Heat Transfer’, International Journal of Heat and Fluid Flow, Vol. 24, N°1, pp. 29–40, 2003.

[6] W.C Reynolds, ‘Computation of Turbulent Flow’, Annual Reviews in Fluid Mechanics, Vol. 8, pp. 183–208, 1976.

[7] P.A. Durbin and D. Laurence, ‘Nonlocal Effects in Single Point Closure’, Third Advances in Turbulence Research Conference, Korea University, Korea, pp. 109–120, 1996.

[8] C.G Speziale, ‘Analytical Methods for the Development of Reynolds-Stress Closures in Turbulence’, Annual Reviews in Fluid Mechanics, Vol. 23, pp. 107–157, 1991.

[9] B.P. Leonard, ‘A Stable and Accurate Convective Modelling Procedure Based on Quadratic Upstream Interpolation’, Computer Methods in Applied Mechanics and Engineering, Vol. 19, pp. 19–59, 1979.

[10] S.V. Patankar, ‘A Calculation Procedure for Two-Dimensional Elliptic Situations’, Numerical Heat Transfer, Vol. 4, pp. 409–425, 1981.

[11] D.D. Apsley and M.A. Leschziner, ‘A New Low-Reynolds-Number Nonlinear Two-Equation Turbulence Model for Complex Flows’, International Journal of Heat and Fluid Flow, Vol. 19, N°3, pp. 209–222, 1998.

[12] T.S. Park and H.J. Sung, ‘A Nonlinear Low-Reynolds-Number $k-\varepsilon$ Model for Turbulent Separated and Reattaching Flows – (I) Flow Field Computation’, International Journal of Heat and Mass Transfer, Vol. 38, N°14, pp. 2657–2666, 1995.

[13] J. Bousquet, ‘Méthode des Singularités’, Cepadues–Editions, Septembre 1990.

[14] J.L. Hess, ‘Review of Integral-Equation Techniques for Solving Potential-Flow Problems with Emphasis on the Surface Method’, Computer Methods in Applied Mechanics and Engineering, Vol. 5, N°2, pp. 145–196, 1975.

[15] G.V. Hobson and B. Lakshminarayana, ‘Prediction of Cascade Performance using an Incompressible Navier-Stokes Technique’, Transactions of ASME, Journal of Turbomachinery, Vol. 113, pp. 561–572, 1991.
[16] G. Jin and M. Braza, 'A Two-Equation Turbulence Model for Unsteady Separated Flows Around Airfoils', J. AIAA, Vol. 32, N°11, pp. 2316-2320, 1994.

[17] T.H. Shih, J. Zhu and J.L. Lumley, 'A Realizable Reynolds Stress Algebraic Equation Model', Rapport Technique, 105993, NASA Technical Memorandum, pp. 47, 68, 69, 73, 1993.

[18] T.S Park and H.J. Sung, 'A New Low-Reynolds-Number $k-\epsilon$-$f_\mu$ Model for Predictions Involving Multiple Surfaces', Fluid Dynamics Research, Vol. 20, N°1-6, pp. 97-113, 1997.

[19] G.N. Coleman and N.N. Mansour, 'Simulation and Modelling of Homogenous Compressible Turbulence Under Isotropic Mean Compression', Turbulent Shear Flows, Vol. 8, pp. 269-282, 1992.

[20] T.B. Gatski and C.G Speziale, 'On Explicit Algebraic Stress Models for Complex Turbulent flows', Journal of Fluid Mechanics, Vol. 254, pp. 59-78, 1993.

[21] S.B Pope, 'A More General Effective-Viscosity Hypothesis', Journal of Fluid Mechanics, Vol. 72, pp. 331-340, 1975.

[22] S.S Girimaji, 'Fully Explicit and Self-Consistent Algebraic Reynolds Stress Model', Theoretical and Computational Fluid Dynamics, Vol. 8, N°6, pp. 387-402, 1996.

[23] S. Wallin and A.V. Johansson, 'An Explicit Algebraic Reynolds Stress Model for Incompressible and Compressible Turbulent Flow', Journal of Fluid Mechanics, Vol. 403, pp. 89-132, 2000.

[24] J.F Thompson, Z.U.A Warsi and C.W Mastin, 'Numerical Grid Generation', Foundations and Applications, North Holland, New York, 1985.

[25] T.J. Craft, B.E. Launder and K. Suga, 'Development and Application of a Cubic Eddy - Viscosity Model of Turbulence', International Journal of Heat and Fluid Flow, Vol. 17, N°2, pp. 108-115, 1996.

[26] R.L. Evans, 'Freestream Turbulence Effects on Turbulent Boundary Layers in Adverse Pressure Gradients', AIAA Journal, Vol. 23, N°11, pp. 1814-1816, 1985.

[27] D.P. Lockard, L.S. Luo, S.D. Milder and B.A. Singer, 'Evaluation of Power FLOW for Aerodynamic Applications', Journal of Statistical Physics, Vol. 107, N°1/2, 2002.