MODELOWANIE INŻYNIERSKIE ISNN 1896-771X 32, s. 205-210, Gliwice 2006
MODELOWANIE NUMERYCZNE OPŁYWU PROFILU
KRZYSZTOF JAN JESIONEK
TOMASZ KOZŁOWSKI
Zakład Maszyn Przepływowych, Politechnika Wrocławska
Streszczenie. Praca przedstawia wyniki obliczeń numerycznych opływu profilu o znanych charakterystykach – geometrycznej i aerodynamicznej. Obliczenia realizowano według następujących modeli turbulencji: standardowego k–ε, rozszerzonej wersji modelu standardowego, RNG k–ε, Reynolds stress oraz modelu z niską liczbą Reynoldsa k–ω.
Dla powyższych modeli turbulencji wyznaczono aerodynamiczne charakterystyki obliczeniowe, które porównano z wynikami badań eksperymentalnych.
Przedstawiono analizę przydatności poszczególnych modeli w rozwiązywaniu zagadnień dotyczących opływu płata.
1. WSTĘP
Przeprowadzone obliczenia numeryczne dla wybranych modeli turbulencji – przy wykorzystaniu komercyjnego programu FLUENT – dotyczą opływu profilu NACA 4412, rys.
1. Wyznaczone charakterystyki aerodynamiczne badanego płata porównano z wynikami badań eksperymentalnych. Przeanalizowano przydatność poszczególnych modeli w rozwiązywaniu zagadnień dotyczących opływu płata.
Rys. 1. Geometria profilu NACA 4412, [7]
Biorąc pod uwagę dużą liczbę modeli turbulencji oferowaną przez współczesne kody numeryczne oraz pełną swobodę przy ich wyborze, należy podkreślić, że odpowiedni wybór modelu turbulencji dla danej klasy przepływów może mieć kluczowe znaczenie dla dokładności i zbieżności obliczeń.
W pracy przeanalizowano również wpływ dyskretyzacji członu konwekcyjnego na dokładność obliczeń numerycznych. Przebadano kilka znanych schematów (wstecz, Lax–
Wendroff, Fromm, Beam–Warming, Godunow). Na podstawie obliczeń jednowymiarowej
konwekcji, oraz przepływu przez kanał zbieżno–rozbieżny określono dokładność wybranych schematów dyskretyzacji.
2. OPŁYW PROFILU NACA 4412
2.1. Wybrane modele turbulencji
Obliczenia zrealizowano na podstawie następujących modeli turbulencji – Spalarart, Allmaras
– k–ε – RNG k–ε – k–ω
– naprężeń Reynoldsa.
Poniżej zestawiono dwa z nich (zastosowano konwencję sumacyjną po j). Często obecnie stosowany model k–ε:
σ ρε µ µ ρ
ρ + −
∂
∂
+
∂
= ∂
∂ + ∂
∂
∂ G
x k x
x u k t k
j k t j
j
j (1)
C k kG x C
x u x
t j
t j
j j
2 2 1
ρε ε
ε σε µ µ ρ ε
ρ ε + −
∂
∂
+
∂
= ∂
∂ + ∂
∂
∂ , (2)
j i i j j i
t x
u x u x G u
∂
∂
∂ +∂
∂
=µ ∂ , (3)
ρ ε µ µ k2
t =C , (4)
gdzie: k – energia kinetyczna turbulencji ε – strata energii turbulencji µt – lepkość turbulentna
G – produkcja energii turbulencji
C1,C2,Cµ,σk,σ – współczynniki empiryczne [1], [2], [4], [6]
oraz jednorównaniowy model Spalararta–Allmarasa, w którym wartości współczynników C1, C2, C3, znaleźć można w pracach [1], [2], [4].
+
+
∂ + ∂
∂
∂
+
∂
= ∂
∂ + ∂
∂
∂
d f u C uG x C
C u x
u u x
x u u t u
j j
k j
j
j ρ ρ ρ
σ µ ρ ρ
ρ 2 3
2
1 . (5)
2.2. Wyniki obliczeń
Obliczenia numeryczne przeprowadzono dla profilu NACA 4412 dla danych: U = 88 m/s, ν
= 1.46.10–6 m2/s, ρ = 1,25 kg/m3 i cięciwy profilu b = 0,5 m. Zastosowano niestrukturalną siatkę trójkątną.
Na rys. 3 zestawiono obliczeniową i doświadczalną [7] charakterystykę aerodynamiczną omawianego płata, dla wybranych modeli turbulencji. Z wykresu wynika, że najlepszą zgodność z eksperymentem uzyskuje się dla jednorównaniowego modelu – Spalararta i Allmarasa, (podobnie jak w [1] ) oraz pięciorównaniowego – naprężeń Reynoldsa.
Największe rozbieżności między wynikami obliczeń i eksperymentem obserwuje się dla wy- sokich wartości – zarówno dodatnich jak i ujemnych – kąta natarcia. Należy dodać, że żaden z rozważanych modeli turbulencji nie pozwolił na pełne odzwierciedlenie warunków przepły- wowych panujących przy ujemnych kątach natarcia. Przyczyną może być wykorzystany schemat dyskretyzacji członu konwekcyjnego – WSTECZ – pierwszego rzędu. W przepły- wach, w których występują duże odchylenia wektora prędkości u od głównego kierunku ruchu (a więc przy wirach pojawiających się dla dużych kątów natarcia, Rys. 2) wpływ dyfuzji numerycznej wprowadzanej sztucznie wraz z schematem WSTECZ, może znacznie oddziały- wać na wyniki obliczeń, [6]. Należy dodać, że wpływ ten rośnie wraz ze wzrostem prędkości u, a więc również przy wzroście prędkości wiru. Na tej podstawie można przypuszczać, dlac- zego obliczeniowe ekstrema współczynnika siły nośnej są przesunięte o kilka stopni względem eksperymentalnych.
a) b) c)
Rys. 2. Opływ płata NACA 4412, Re= 3.106, kąt natarcia odpowiednio:
a), α = 00 a), α = 240 a), α = –160
-1,2 -0,8 -0,4 0,0 0,4 0,8 1,2 1,6 2,0
-32 -24 -16 -8 0 8 16 24 32
Rys. 3. Krzywa współczynnika siły nośnej CL
0 0,1 0,2 0,3 0,4 0,5
L,m
EKSPERYMENT WSTECZ QUICK -2
-1 0 1 2
p
Rys. 4. Bezwymiarowy rozkład ciśnienia na profilu NACA 4412, α = 80, Re= 3.106.
Rozwiązaniem omawianego problemu może być wykorzystanie innego schematu, np. QUICK, [2], [6], w którym sztucznie wprowadzana dyfuzja numeryczna jest znacznie mniejsza. Schemat ten – dla pewnych klas przypływów – może się jednak charakteryzować pewnym „przestrzelaniem” wyników. Na rys. 4 zestawiono wyniki obliczeń ciśnienia na profilu NACA 4412 dla schematów dyskretyzacji członu konwekcyjnego: WSTECZ i QUICK.
3. WPŁYW DYSKRETYZACJI CZŁONU KONWEKCYJNEGO NA DOKŁADNOŚĆ OBLICZEŃ
3.1. Jednowymiarowa konwekcja
Wpływ schematu dyskretyzacji na dokładność obliczeń przeanalizowano dla przypadku jednowymiarowej konwekcji
=0
∂ + ∂
∂
∂
x U u t
u . (6)
Dla równania (6) wykorzystano schemat postaci
( ) ( )
−
− +
−
−
= − −
+ n
i n i n
i n i n
n u dt
dxdt U u
u dx Udt u
u 1
2 2 1
1
2 2
1 φ φ . (7)
Tabela 1. Schematy dyskretyzacji członu konwekcyjnego, [8]
WSTECZ FROMM LAX–WENDROFF BEAM–WARMING
=0
n
φi
dx u uin in
n
i 2
1
1 −
+ −
=
φ dx
u uin in
n i
= +1−
φ dx
u uin in
n i
−1
= − φ Dokładny schemat GODUNOVa, wraz z limiterami znaleźć można w [8].
Na rys. 5 przedstawiono porównanie dokładności omawianych schematów konwekcji. Dla aproksymacji WSTECZ wyraźnie widoczne jest rozmywanie rozwiązania, związane z wprowadzaną sztucznie dyfuzją numeryczną. Schematy LAX–WENDROFFa oraz BEAM–
WARMINGa są mniej dyfuzyjne, ale charakteryzują się pewnym przestrzelaniem rozwiązań.
Podobnie schemat FROMMa, który jest liniową kombinacją dwóch powyższych. Najmniejszy błąd uzyskuje się dla aproksymacji GODUNOVa [8], dlatego jest on często stosowany w obliczeniach numerycznych równań konwekcyjnych. Charakteryzuje się małą dyfuzją numeryczną i mniejszym przestrzelaniem rozwiązań niż często stosowany schemat QUICK.
Rys.5. Porównanie dokładności schematów jednowymiarowej konwekcji, z prędkością U = 1 m/s, po czasie 1s. Liczba Couranta C = 0.5
3.2. Kanał zbieżno–rozbieżny
Obliczenia numeryczne przeprowadzono na równaniach zapisanych dla krzywoliniowego układu współrzędnych
=0
∂
∂ J u x
j
j
, (8)
∂∂
( )
+ ∂∂ +{ }
+ ∂∂ = ∂∂ ∂∂ l j
i il i
ji k j i
j
i j j
x u J g x
x J g p u u k i J j
u u Ju x
t u
ν 1 , (9) gdzie:
uj, ui uk – kontrawariantne składowe prędkości
J – jakobian przekształcenia układu współrzędnych
{ }
jik – christoflangij, gil – tensory metryczne wzajemne
na podstawie metody korekcji ciśnienia [1], [4], [5], [6], dla danych geometrycznych:
wysokość kanału H = 1 m, minimalna wysokość przewężenia H0 = 0,5 m, długość kanału L = 1 m.
Obliczenia przeprowadzono w programie MATLAB. Wybrano niską liczbę Reynoldsa, z uwagi na stabilność schematu WSTECZ, zapisanego w podstawowej formie.
Na rys. 6 zestawiono wyniki obliczeń dla kanału zbieżno–rozbieżnego. Ponieważ celem jest jedynie porównanie schematów dyskretyzacji członu konwekcyjnego, jako wzór przyjęto rozwiązanie uzyskane w oparciu o program FLUENT. Z wykresu wynika, że większe wartości błędu bezwzględnego w całym obszarze obliczeń otrzymuje się dla schematu GODUNOVa.
Jednak w strefie występowania wiru (oznaczonej kółkiem), trend linii rozkładu ciśnienia na ściance jest bliższy „wzorcowemu” niż dla schematu WSTECZ. Sugerują to również otrzymane rozkłady prędkości.
Rys.6. Rozkład ciśnień na ściankach kanału zbieżno–rozbieżnego, Re = 70
Rys.7. Kontury prędkości przy przepływie przez kanał zbieżno–rozbieżny, Re = 70
Większa wartość błędu bezwzględnego otrzymana dla schematu GODUNOVa wynika z zaimplementowania schematu w formie stosowanej dla siatki równomiernej, a nie siatki przestawnej (straggered grid), na której prowadzono obliczenia.
4. PODSUMOWANIE
W pracy przedstawiono obliczenia numeryczne opływu profilu NACA 4412, dla wybranych modeli turbulencji. Wyznaczono charakterystyki aerodynamiczne badanego płata, które porównano z wynikami badań eksperymentalnych. Podobnie jak w [1], najlepszą zgodność z eksperymentem zapewnia model Spalararta i Allmarasa, który bywa często stosowany w obliczeniach aerodynamicznych profili. Analiza położenia ekstremów obliczeniowego współczynnika siły nośnej względem wartości eksperymentalnych, pozwala przypuszczać, że na błąd obliczeń znaczny wpływ może mieć wprowadzana wraz ze schematem konwekcji WSTECZ, dyfuzja numeryczna. Niekorzystne zjawisko sztucznej lepkości, intensyfikowane jest wraz ze zwiększeniem odchylenia wektora prędkości u od głównego kierunku ruchu. Dzięki zastosowaniu schematów np.: QUICK lub GODUNOVa można znacznie zmniejszyć lepkość numeryczną, poprawiając tym samym dokładność obliczeń numerycznych.
LITERATURA
1. Chung T.: Computational fluid dynamics. Cambridge University Press, Cambridge 2003.
2. FLUENT 6.1 User's Guide.
3. Fortuna Z., Macukow B., Wąsowski J.: Metody numeryczne. WNT, Warszawa 1998.
4. Gryboś R.: Podstawy mechaniki płynów. Wydawnictwo Naukowe PWN, Warszawa 1998.
5. Hirt C. W., Nicholas B. D., Romero N.C.: Sola – a numerical solution algorithm for tran- sient fluid flows. Report Los Alamos Scientific Laboratory of the University of California, USA 1975.
6. Kazimierski Z.: Numeryczne wyznaczanie trójwymiarowych przepływów turbulentnych.
PAN, Wrocław 1992.
7. Pinkerton R. M.: Calculated and measured pressure distributions over the midspan section of the NACA 4412 Airfoil. Report 563, 1936
8. Warburton T.: Numerical methods for partial differential equations. Department of CAAM at Rice University – www.caam.rice.edu/~caam452.
NUMERICAL CALCULATION OF FLOW AROUND AERODYNAMICS PROFILE
The results of numerical calculations of the flow around aerodynamics profile – for well–known characteristics (geometrical and aerodynamic) – are presented.
These calculations were realized in following turbulence models: the standard k–ε, modified version of standard model, RNG k–ε, Reynolds stress and low Reynolds number k–ω. For above mentioned turbulence models the theoretical aerodynamics characteristics were found, which were compared with results of experimental in- vestigations. The analysis of particular models usefulness in solving aerodynamics profile flow problems is presented.
