La distribution des points du maillage à l’intérieur est gouvernée par les fonctions de

condensation « Stretching functions » aux frontières.

II.2.1 Fonctions de condensation « Stretching functions » unidimensionnelles

Les fonctions de condensation sont largement utilisées pour la distribution des points
le long des frontières particulières ainsi que des régions spécifiques (à forts gradients) du
domaine qui nécessitent un traitement avec une grande précision (couches limites au
voisinage de la paroi, ondes de choc,…etc.) Pour les écoulements externes autour des corps, il
est nécessaire d’introduire ces fonctions pour résoudre les problèmes où les gradients des
variables physiques sont larges.

Il est désirable d’exprimer les variables dépendantes et indépendantes dans la fonction

de condensation sous une forme normalisée ([Link].1).
Dans ce cas, la variable independante appropriée est :

η −η1
η *= (II.1)
η 2−η1

avec : 0≤η ≤1 η ≤η ≤η
*
ou
1 2

Une fonction de condensation « stretching » proposée par Robert (1971) et modifiée

par Eisemann (1979), est la suivante :

 
 tanhQ1−η  
 *

S = Pη
*
+(1− P)1−    
 tanh (Q ) 
(II.2)
 
 
Où : P et Q sont des paramètres de condensation des lignes coordonnées.

II.2.2. Techniques de deux parois

Cette méthode est généralement illustrée par un canal bidimensionnel courbe

([Link].1). Il est souvent nécessaire que les fonctions de condensation SAD(η*) et SBC(η*)

2
soient définies pour contrôler la distribution des points au-dessus et au-dessous des frontières.
Des équations équivalentes à celles de (II.2) pour générer SAD(η*) et SBC(η*). Pour obtenir la
valeur de S entre les surfaces AD et BC, une simple interpolation linéaire est recommandée,
d’où :

S= SAD +ξ* (SBC-SAD) (II.3)

ξ −ξ 1
Où ξ *= .
ξ 2−ξ 1

η=η2 D

y B

ξ=ξ2

x A
η=η1 ξ=ξ1

Fig.(II.1) Canal bidimensionnel curviligne

De la même façon, la distribution des points le long des frontières AB et CD contrôlée

par les fonctions de condensation rDC(ξ*) et rAB(ξ*). Si rAB et rDC sont considérées comme étant
les coordonnées généralisées mesurées le long de la surface, alors XAB(rAB) et YAB(rAB) se
suivent directement. Il est de même pour XDC(rDC) et YDC(rDC).

La technique des deux parois utilise des moyens d’interpolation pour générer les
points intérieurs entre les deux frontières AB et DC. Une simple interpolation linéaire est
donnée dans ce qui suit :

X (ξ,η )=(1−S) X AB(r AB )+ S X DC (r DC )

Y (ξ,η )=(1− S)Y AB(r AB )+ S Y DC (r DC )
(II.4)

où S est donnée par(II.2).

3
 La difficulté dans l’interpolation précédente réside dans le fait que les points du
maillage, adjacents à la frontière, peuvent devenir tordus si les points sur les frontières
correspondantes (XAB, YAB) et (XDC, YDC) sont en dehors de l’alignement. Par conséquent, il
est préférable de remplacer (II.4) par :

X (ξ,η )= µ1(s )X AB(r AB )+ µ 2(s )X DC (r DC )+T 1µ 3(s ) Y AB (r AB )+T 2µ 4(s ) Y DC (r DC )

d d
dr AB dr DC
Y (ξ,η )= µ1(s )Y AB(r AB )+ µ 2(s )Y DC (r DC )+T 1µ 3(s ) (r AB )+T 2µ 4(s ) (r DC )
d X AB d X DC
dr AB dr DC

µ (s )=2s −3s²+1,
1
3
µ (s )=−2s +3s²
1
3 (II.5)

où:
µ (s )=2s −3s²+1,
1
3
µ (s )=−2s +3s²
1
3

II.2.3. Technique des surfaces intermédiaires (multisurfaces)

L’introduction des surfaces intermédiaires entre les deux surfaces limites AB et CD
peut améliorer d’une manière considérable le contrôle de la distribution des nœuds internes.

La méthode multisurfaces permet d’obtenir un maillage localement orthogonal surtout

à la surface limite AB du profil d’aile. En principe, il n’existe pas une limitation concernant le
nombre des surfaces intermédiaires à introduire. En pratique, un contrôle efficace des points
internes peut être obtenu avec seulement deux surfaces intermédiaires. Une adjacente à la
surface limite intérieure AB et l’autre adjacente à la surface limite extérieure DC.

II.2.4 Type des profils étudiés

Les profils considérés dans cette étude sont de type NACA quatre chiffres
NACAxxyy.
Le premier chiffre ‘x’ indique la cambrure (la flèche relative) donnée en pourcentage de la
corde. Le deuxième ‘x’ désigne la position en dizaine de la corde de la flèche relative, tandis
que les deux derniers chiffres ’yy’ désignent l’épaisseur relative maximale donnée en
pourcentage de la corde.

En premier lieu, on s’intéresse seulement aux profils symétriques dont la ligne

moyenne se confond avec la corde du profil. Ces profils sont désignés par NACA00yy.
La distribution de l’épaisseur est donnée par la fonction suivante :

4
 y e=e(a1x1/ 2+a 2 x + a3 x 2+a 4 x3+ a5 x4 ) (II.6)
où :
a1= 1.4779155, a2= -0.624424, a3= -1.727016, a4 =1.384087, a5= -0.510563
« e » représente l’épaisseur maximum du profil par rapport à la corde (e=yy).

Pour notre cas test nous avons choisi le profil symétrique NACA0012 ([Link].2) avec
une épaisseur relative de l’ordre de 12%. Bien que sa forme n’est pas idéale dans le domaine
du dessin aérodynamique, il est extrêmement utile comme référence standard parce qu’il a été
amplement testé numériquement et expérimentalement.

0.30

0.26

0.22

0.18

0.14

0.10

0.06

0.02

-0.02
0.00 0.20 0.40 0.60 0.80 1.00
-0.06

-0.10

-0.14

-0.18

-0.22

-0.26

-0.30

[Link].2. NACA0012

II.3 Adimensionnement des variables

Dans la dynamique des fluides, l’adimensionnement des équations gouvernant le
problème physique est d’une grande importance. Il est convenable, avant toute resolution
numérique d’obtenir des équations adimensionnelles et ça pour diminuer les erreurs
d’arrondis dues aux valeurs dimensionnelles de grandeurs physiques. Par conséquent, tous les
paramètres sont adimensionnés en respectant les grandeurs caractéristiques de l’écoulement, à

5
savoir la vitesse de référence U∞ (vitesse à l’infini amont loin du profil), la longueur
caractéristique C (corde de profil), la densité ρ∞ et la viscosité de référence µ∞.
Il faut noter que l’adimensionnement des variables n’affecte pas la forme des équations,
c’est-à-dire on obtient le même système d’équations original (voir rapport N°1).

Les paramètres adimensionnés utilisés dans notre étude sont :

~ ρ ~ U j , k~ = k , ε~ = εC , µ~ = µ p− p∞
, ν~ = ν ℜe∞ , ~
x j ~ tU
ℜ
−1 −1
x j = C , t = C∞ , ρ = , U j =
~ p= ,
ρ∞ U∞ U∞
2
U∞
3
µ∞ e∞
ν∞ ρ ∞U ∞2
où :
C
ℜ e∞
est le nombre de Reynolds à l’infini amont, ℜ e∞
=U ∞ , et p est la pression
ν∞
statique locale.
Pour ne pas alourdir la notation, on élève le telda (∼) par la suite. L’adimensionnement décrite
ci-dessus permet d’obtenir exactement les mêmes équations originales suivantes :

Equation de continuité :

∂U i
=0 (II.7)
∂ xi

Equations de quantité de mouvement :

∂U i ∂  ∂ ∂U  
+U j U i =− 1 ∂P + ∂ ν  U i + j −u iu j  (II.8)
∂t ∂ x j ρ ∂ xi ∂ x j   ∂ x j ∂ x i  

Le modèle k-ε :

∂k +U ∂k = ∂ ν + ν t  ∂k  −u u ∂U i −ε

∂ x j ∂ x j  σ k  ∂ x j  i j ∂ x j
(II.9a)
∂t j

∂ε +U ∂ε = ∂ ν + ν t  ∂ε −c ε u u ∂U i −c ε²

∂t j
∂ x j ∂ x j  σ ε  ∂ x j  ε1 k i j ∂ x j ε 2 k (II.9b)

Le système d’équations (II.7-II.9) gouvernant l’écoulement incompressible turbulent en

utilisant le modèle k-ε peut être écrit sous la forme conservatrice suivante :

∂ (ρU φ )= ∂  Γ ∂φ + S (II.10)

∂x j j
∂x j ∂x j 

Le tableau (II.1) résume les différentes variables et les termes sources :

6
 Equation φ Γ S
continuité 1 0 0
Quantité de mouvement Ui µ eff ∂p
−
∂ xi

Energie turbulente k k µt Pk-ρε

µ+
σk
Dissipation de l’énergie ε µt ε  −c ρε 
µ+
σε k  cε 1P k ε2 

∂U j ∂U j ∂U i  µ = c µ ρk ² / ε + µ
P =µ
k t ∂
 + 
x i  ∂ xi ∂ x j 
eff

[Link].1 Les efférents paramètres dans l’équation générale de transport

Modèle k-ε

Dans notre cas (écoulement incompressible) et sous sa forme adimensionnelle on a ρ=1. Par
conséquent µ=ν.

II.4 Transformations des équations en coordonnées généralisées (cas

bidimensionnel)

Souvent, avant d’appliquer les algorithmes de résolution numérique, l’équation

générale de transport ( II.10) doit être transformée du domaine physique de coordonnées
cartésiennes (x,y) vers le domaine computationnel (de calcul) de coordonnées (ξ,η).
La raison principale de ce changement de domaine et l’application de la transformation des
variables indépendantes à l’équation générale de transport est de transformer plusieurs
surfaces géométriques complexes dans le problème à des lignes de coordonnées constantes
(maillage rectangulaire uniforme). Ainsi, l’application des conditions aux limites peut être
acquise sans difficulté.

7
 (x, y)↔ (ξ ,η )

Plan computationnel (ξ ,η ) plan physique (x,y)

[Link].3 Transformation de maillage du plan physique au plan computationnel

Pour un maillage stationnaire, en bidimensionnel, cette transformation est exprimée

par :
ξ =ξ (x, y )
η =η(x, y ) (II.11)

Dans ce cas ξ,η représentent des axes orthogonaux dans le plan computationnel.
Avant que les équations soient résolues dans le domaine (ξ,η) doivent être écrites comme
étant fonctions des nouvelles variables en espace. Pour atteindre cet objectif, on considère la
règle suivante :
∂ =ξ ∂ +η ∂
∂x x ∂ξ x ∂η (II.12)
∂ =ξ ∂ +η ∂
∂y y ∂ξ y ∂η

 = ∂ξ ,....etc. .
 ξ x ∂x 
 
L’équation (II.12) permet de transformer l’équation générale de transport au plan
computationnel mais, il faut tout d’abord expliciter les paramètres métriques de la
transformation (ξx, ηy,…etc.) à l’aide des équations différentielles suivantes :

dx= xξ dξ + xη dη
dy = yξ dξ + yη dη (II.13)

Ainsi de même pour :

dξ =ξ x dx+ξ y dy
dη =η x dx+η y dy (II.14)

Sous la forme matricielle (II.13) et (II.14) s’écrivent :

8
 dx  xξ xη   dξ 
 = .  (II.15)
dy   yξ yη  dη 


 dξ  ξ x ξ y  dx 
  = .  (II.16)
dη  η x η y  dy 

Puisque les transformations sont inverses l’un par rapport à l’autre, on a :

−1

ξ x ξ y   xξ xη   yη − yξ 
 =  = J.  (II.17)
η x η y  − xη xξ 
   yξ yη   
 
Où J est le jacobien de la transformation.
Afin d’évaluer les quantités paramétriques, les identités numériques suivantes sont
demandées :

(
J = xξ yη − xη yξ )
−1

ξ x= J yη , η x=−J yξ (II.18)

ξ y =−J xη , η y = J xξ

Finalement l’équation générale de transport (II.10) en bidimensionnel peut être écrite sous la
forme :
φ
1 ∂ (ρφ )+ ∂ (ρU φ )+ ∂ (ρV φ )= ∂  ΓJα ∂φ −ΓJβ ∂φ + ∂  ΓJγ ∂φ −ΓJβ ∂φ − S (II.19)
J ∂t ∂ξ c
∂η c
∂ξ  ∂ξ ∂η  ∂η  ∂η ∂ξ  J
où :
U c =U yη −V xη
(II.20)
V c =V xξ −U yξ

Uc et Vc désignent les vitesses contravariantes selon les axes ξ et η respectivement, et Sφ est le

terme source généralisé dans le nouveau système de coordonnées (ξ,η).
Les paramètres métriques avec le jacobien de la transformation sont exprimés par :

9
 α = xη2 + yη2

γ = xξ2 + yξ2
(II.21)
β = xξ xη + yξ yη

(
J = xξ yη − xη yξ )
−1

Remarque :
Pour le schéma numérique nous avons appliqué la méthode des volumes finis avec un
maillage décalé dont les principaux détails seront évoqués dans le prochain rapport.
Nous présentons dans ce qui suit quelques résultats obtenus concernant la génération de
maillage et les résultats relatifs au modèle k-e.

Conclusion
La partie de génération de maillage permet d’obtenir différentes configurations de la
limite externe telles que : maillage en C et maillage en O. Ainsi, on peut obtenir des maillages
triangulaires utilisés en particulier dans les méthodes des éléments finis. Le raffinement de
maillage à la paroi et au sillage est possible à l’aide des fonctions de condensation.

Les quelques résultas présentés à savoir la distribution de cœfficient de pression par

exemple montrent bien que nous avons un problème au bord de fuite. Alors donc, avant toute
validation avec d’autres résultas il est convenable de régler tout d’abord ce problème.

10
[Link].4. Maillage uniforme en C 11x11 pour l’extrados de NACA0012

[Link].5 Maillage uniforme en C 21x11 pour le profil d’aile NACA0012

11
 [Link].6 Raffinement du maillage (101x51) au sillage et à la paroi
(maillage type en C)

[Link].7 Zoom de maillage raffiné à la paroi et au sillage (bord de fuite)

12
 a)

b)

[Link].8 Maillage quadrilatéral structuré en O

13
 a)231 nœuds et 400 triangles

b) 5151 nœuds et 10000 triangles

[Link].9 Maillage triangulaire structuré en C

14
 a) 231 nœuds et 400 triangles

b)5151 nœuds et 10000 triangles

[Link].10 Maillage triangulaire structuré en O

15
 -1.00

0.00

Cp

1.00

2.00

0.00 0.20 0.40 0.60 0.80 1.00

x/c

[Link].11 Coefficient de pression Cp(x/C)

Angle d’attaque α=0°. (Re∞=105)

[Link].12 Carte isopression p/p∞

Angle d’attaque α=0°. (Re∞=105)

16
 Cp
-1.00 Extrados
Intrados

0.00
Cp

1.00

2.00

0.00 0.20 0.40 0.60 0.80 1.00

x/c

[Link].13 Coefficient de pression Cp(x/C)

Angle d’attaque α=2°. (Re∞=105)

Frame 001  04 Oct 2005  p

1.01226
1.01176
1.01126
1.01076
1.01026
0.999232 1.00975
1 1.00925
0.998731
1.00875
1.00825
1.00775
1.00725
0.5
1.00675
1.00625
0.9

7
22
99

97

1.00575
73

0.997729
0.9
3

1.00124 1.00524
0 1.00474
y

1.00424
0.998731 1.00374
1.00

1.00324
02

1.00274
3

-0.5 1.00224
1.00174
23

1.00124
1.000

0.99
9 733
1.00073
1.00023
-1
0.99973
0.99923
0.99873
0.99823
0 1 2 0.99772
x 0.99722
0.99672

[Link].14 Carte isopression p/p∞

Angle d’attaque α=2°. (Re∞=105)

17
 Cp
-2.00 Extrados
Intrados

-1.00

0.00
Cp

1.00

2.00

3.00

0.00 0.20 0.40 0.60 0.80 1.00

x/c

[Link].15 Coefficient de pression Cp(x/C)

Angle d’attaque α=6°. (Re∞=105)

[Link].16 Carte isopression p/p∞

Angle d’attaque α=6°. (Re∞=105)

18