MODELOWANIE STRUKTUR KAWITACJI PRZY OPŁYWIE PROFILU

(1)

MODELOWANIE INŻYNIERSKIE 2017 nr 63, ISSN 1896-771X

MODELOWANIE STRUKTUR KAWITACJI PRZY OPŁYWIE PROFILU

Dorota Homa

^1a

, Włodzimierz Wróblewski

^1b

1Instytut Maszyn i Urządzeń Energetycznych, Politechnika Śląska

adorota.homa@polsl.pl, ^bwlodzimierz.wroblewski@polsl.pl

Streszczenie

W artykule przedstawiono wyniki symulacji numerycznych zjawiska kawitacji występującego przy opływie ło- patki profilu ClarkY wodą. Przepływ odbywał się w komorze o przekroju prostokątnym. Kąt natarcia łopaty wy- nosił 8°. Symulacje przeprowadzono przy użyciu kodu OpenFOAM. Dla modelu przepływu dwufazowego prowadzono obliczenia niestacjonarne przy założeniu modelu kawitacji Schnerr & Sauer. Przebadano i porównano dwa modele: model 2D oraz model 3D. Wyniki obliczeń porównano pod względem kształtu i dynamiki zmian chmur kawitacyjnych oraz procentowego udziału pary wodnej w domenie w funkcji czasu. Porównano uśrednione w jednym okresie wartości ciśnienia i udziałów objętościowych pary w punktach monitorujących. Obliczenia dla modelu 3D pozwoliły zaobserwować zmienność struktur kawitacyjnych w kierunku prostopadłym do kierunku przepływu wody.

Słowa kluczowe: kawitacja, opływ łopatki, przepływ dwufazowy, OpenFOAM

CAVITATION STRUCTURES MODELING IN CASE OF FLOW OVER THE FOIL

Summary

The results of numerical simulation of water cavitating flow are shown in the paper, in case of flow over the ClarkY foil. The calculations were performed for flow through a rectangular channel. The angle of attack of the blade was 8°. The simulation was conducted in OpenFOAM software. For the two-phase model the transient calculations were carried out with Schnerr & Sauer cavitation model. The two models: 2D and 3D, were investigated and compared. The results of the calculation were compared regarding the shape and dynamics of cavitation structures and percentage of vapor in the domain in time. The time-averaged pressure and vapor volume fraction values in monitor points located over the foil were compared. The 3D calculations provided information about the fluctuation of cavitation structures in the direction perpendicular to the water flow.

Keywords: cavitation, flow over a foil, two phase flow, OpenFOAM

1. WSTĘP

Symulacje numeryczne dotyczące zjawiska kawitacji występującego podczas przepływu wody pozwalają poszerzać wiedzę na temat rozkładu, kształtu i dynamiki zmian pojawiających się struktur pary wodnej. Dla użytkowników maszyn wirnikowych – pomp, turbin wodnych – istotnym przypadkiem przepływu z kawitacją jest opływ łopatki profilu, ze względu na geometryczne podobieństwo do przepływu w wirnikach maszyn energetycznych. Występowanie w przepływie z kawitacją zapadania się pęcherzy pary prowadzi do erozji kawitacyjnej wirników i łopat zarówno pomp jak i turbin

wodnych [3-4, 8]. W związku z tym eksperymenty oraz symulacje numeryczne dotyczące opływu łopatki w warunkach sprzyjających rozwojowi kawitacji stały się tematem wielu opracowań naukowych. Istotnym zagad- nieniem jest identyfikacja miejsca powstawania struktur parowych oraz obszarów, gdzie struktury odrywają się od powierzchni łopatki i zapadają. Analiza modelu 2D może dostarczyć takich informacji [1, 6]. W [7] rozważa- no przepływ z kawitacją dla profilu NACA0015 dla dwuwymiarowego modelu numerycznego. Zaproponowa- no model, który uwzględniał interakcje między kawita-

(2)

cyjnymi strukturami wirowymi oraz pęcherzami pary zbitymi w klastrze. Wyniki obliczeń zostały przedsta- wione jedynie na płaszczyźnie równoległej do kierunku przepływu cieczy (wzdłuż profilu). Saito [14] przeprowa- dził obliczenia dla profilu NACA0015 z modelem 3D.

Przeanalizowano dwa warianty: z występowaniem po bokach profilu ścian oraz z ich brakiem. Wyniki przed- stawione w [14] wskazują, że ściany ograniczające kanał przepływowy mają duży wpływ na kształt powstającej chmury kawitacyjnej. Dzięki uwzględnieniu ścian bocznych udało się zaobserwować strukturę kawitacyjną o kształcie litery „U”, co nie jest możliwe w przypadku zastosowania modelu dwuwymiarowego.

Celem opisanych w artykule badań numerycznych było porównanie wyników symulacji przepływu z kawitacją dla modelu 2D profilu ClarkY oraz modelu 3D, w któ- rym po jednej stronie profilu przyjęto założenie gładkiej ściany, a po drugiej – symetrii. Obydwa modele zostały przebadane w takich samych warunkach przepływowych oraz z założeniem tego samego modelu kawitacji – Schnerr & Sauer. Dzięki temu można było zaobserwować wpływ rodzaju warunku brzegowego zadanego na ścianie na powstające struktury kawitacyjne oraz zmienność parametrów przepływu w kierunku prostopadłym do kierunku przepływu.

Przeprowadzono symulacje dla warunków pozwalających zaobserwować kawitację warstwową. Wielkością charak- teryzującą przepływ kawitacyjny jest liczba kawitacji σ, definiowana jako:

σ = _. (1)

p – ciśnienie statyczne (absolutne), Pa

ps – ciśnienie nasycenia dla temperatury cieczy, Pa ρ∞ – gęstość cieczy w niezakłóconym przepływie, kg/m³ u∞ – prędkość niezakłóconego przepływu, m/s

Kawitacja warstwowa jest obserwowana, gdy liczba kawitacji spada do wartości ok. 1.4 [16]. Ten typ kawitacji charakteryzuje się częstotliwością powstawania i zapadania się struktur kawitacyjnych równą 200 Hz (okres zmian wynosi 5 ms) [16]. Powstające pęcherze pary przy krawędzi natarcia układają się w kształt litery U. Przepływ w obszarze krawędzi spływu profilu jest silnie niestacjonarny. Struktury kawitacyjne przekształ- cają się w wiry, które odrywają się okresowo.

2. MODEL MATEMATYCZNY KAWITACJI

Model kawitacji zastosowany w badaniach należy do grupy modeli typu one-fluid. Charakteryzują się one tym, że płyn jest traktowany jako mieszanina fazy dyspersyjnej (gazowej) i ciągłej. W modelach tego typu rozwiązywane są równania zachowania dla mieszaniny

(zachowania masy, pędu) oraz równanie zachowana masy dla fazy gazowej. W użytych modelach założono brak poślizgu między fazami, co oznacza, że prędkości fazy ciekłej i gazowej są równe. Udział pary wodnej w płynie jest wyznaczany z zależności [12]:

+ ∇ αρ u = R − R (2) gdzie:

α – udział objętościowy pary wodnej, ρv – gęstość pary wodnej, kg/m³ u – prędkość mieszaniny, m/s t – czas, s

Re, Rc – człony źródłowe równania zachowania masy, kg/m³s

Człony źródłowe opisują strumień wymiany masy mię- dzy fazami odniesiony do jednostki objętości fazy ciekłej.

W modelu Schnerr & Sauer określane są na podstawie równania Rayleigha-Plesseta (RP), które opisuje dyna- mikę wzrostu i zaniku pęcherzy pary wodnej [2]. Używa się uproszczonej formy równania RP z pominięciem napięcia powierzchniowego oraz członów drugiego rzędu.

Równanie to w uproszczonej formie ma postać:

= !^"_#^$ _%^$ (3)

rB – promień pęcherza pary wodnej, m pB – ciśnienie w pęcherzu pary wodnej, Pa ρl – gęstość fazy ciekłej, kg/m³

Człon źródłowy Re można wyznaczyć z [15]:

R = ^% (4)

Natomiast zależność między udziałem objętościowym pary i promieniem pęcherza opisuje równanie [10]:

α =_*+&^& ^'⁽⁾' ⁽

()⁽ (5)

nB – liczba pęcherzy na jednostkę płynu, 1/m³

Na podstawie zależności (5) określany jest promień pęcherza:

,-= ._{* /}^/ ₀₁^#₂^*₃4⁵⁽ (6) Ostatecznie człony źródłowe równania zachowania masy dla pary przybierają postać [15]:

67= ^{8 9}α 1 − α _;^#₃!^"_# ₈ =>≥ = (7) 6@= ^{8 9}α 1 − α _;^#₃!^"_# ₈ =>< = (8)

(3)

MODELOWANIE STRUKTUR KAWITACJI PRZY

3. MODEL SYMULACJI 2D i 3D

Tabela 1. zawiera wybrane parametry symulacji 2D i 3D. Obliczenia prowadzono przy użyciu oprogramowania OpenFOAM. Wybrany solver to interPhaseChangeFoam – uwzględniający zmianę fazy z ciekłej na gazową z zadanym modelem kawitacji [5]. Obliczenia były nie cjonarne, z krokiem czasowym wynoszącym Wartość kroku czasowego została dobrana na podstawie wyników wcześniejszych symulacji, uwzględniających warunek maksymalnej wartości liczby Couranta poniżej 1. Parametry modelu zostały dobrane na podstawie badań literaturowych [6, 9-11, 13].

Tab. 1. Podstawowe parametry modelu Długość cięciwy c 70 mm

Kąt natarcia 8°

Przepływ ciepła Przepływ izotermiczny Temperatura płynu 20°C

Model turbulencji k-ω SST Warunek brzegowy

ściana górna/dolna

Ściana gładka

Warunek brzegowy wylot 72 000 Pa Warunek brzegowy wlot Prędkość 10 m/s Intensywność turbulencji

na wlocie

5%

Liczba Reynoldsa 700 000

nB 1.6 x 10¹³

rB 10^-6m

Krok czasowy Δt 10^-5 s

Istotną różnicą między modelami numerycznymi była geometria domeny obliczeniowej oraz warunek brzegowy zadany dla ścian bocznych kanału przepływowego.

obydwu przypadkach generowanie siatki rozpoczęto od utworzenia płaskiej siatki blokowej. Następnie siatkę wyciągnięto w kierunku prostopadłym na szerokość 1 mm (jeden element) dla obliczeń 2D oraz na szerokość 35 mm dla obliczeń 3D. Przy ścianie siatkę zagęszczono w celu poprawnego odwzorowania zjawisk przyście nych. Rys. 1. przedstawia siatkę powierzchniową z zaznaczonym obszarem wokół profilu.

Rys. 1. Geometria przyjęta do obliczeń

MODELOWANIE STRUKTUR KAWITACJI PRZY OPŁYWIE PROFILU

SYMULACJI 2D i 3D

parametry symulacji 2D i Obliczenia prowadzono przy użyciu oprogramowania OpenFOAM. Wybrany solver to interPhaseChangeFoam uwzględniający zmianę fazy z ciekłej na gazową z . Obliczenia były niestacjonarne, z krokiem czasowym wynoszącym 10^-5 s.

Wartość kroku czasowego została dobrana na podstawie , uwzględniających warunek maksymalnej wartości liczby Couranta poniżej zostały dobrane na podstawie

Przepływ izotermiczny

ciana gładka

Prędkość 10 m/s

numerycznymi była oraz warunek brzegowy zadany dla ścian bocznych kanału przepływowego. W obydwu przypadkach generowanie siatki rozpoczęto od tworzenia płaskiej siatki blokowej. Następnie siatkę padłym na szerokość dla obliczeń 2D oraz na szerokość siatkę zagęszczono w celu poprawnego odwzorowania zjawisk przyścien-

siatkę powierzchniową rofilu.

Na rys. 2. przedstawiony został

siatką do obliczeń 3D. Wymiary kanału wynoszą: x=700 mm, y=189 mm, z=70 mm. Domena obl

stanowi połowę kanału (szerokość

brzegowe zadane na powierzchniach ograniczających:

wlot płaszczyzna YZ x=0, wylot

x=700 mm; ściana płaszczyzny XZ y=0 oraz y=189 mm, ściana płaszczyzna XY z=0;

zna XY z=35 mm.

Rys. 2. Geometria do obliczeń 3D

4. WYNIKI

4.1. SYMULACJA 2D

W celu określenia okresu wzrostu i zapadania się struktur kawitacyjnych podczas obliczeń monitorowano całkowitą objętość pary wodnej w

czasu. Na rys. 3. pokazano przebieg zawartości objętości pary wodnej w

zmian to 3 ms. Do dalszej analizy wybrano okres zmian ograniczony na rys. 3. pionowymi liniami.

Rys. 3. Przebieg procentowego udziału pary wodnej w przepł wie w funkcji czasu, obliczenia 2D

Na rys. 4. przedstawiono przebieg zjawiska wybranego okresu (T – okres zmian, t początku okresu). Można zauważyć

chmury kawitacyjnej w pobliżu krawędzi natarcia, jej wzrost do ok. 0.4 okresu oraz stopniowe zapadanie się.

Podczas gdy jedna struktura ulega zapadnięciu w pobliżu krawędzi spływu, kolejna zaczyna rosnąć przy krawędzi natarcia. Na rys. 5. przedstawiono rozkład

OPŁYWIE PROFILU

kontur kanału wraz z Wymiary kanału wynoszą: x=700

=70 mm. Domena obliczeniowa szerokość 35 mm). Warunki brzegowe zadane na powierzchniach ograniczających:

wylot płaszczyzna YZ płaszczyzny XZ y=0 oraz y=189

=0; symetria płaszczy-

W celu określenia okresu wzrostu i zapadania się struktur kawitacyjnych podczas obliczeń monitorowano całkowitą objętość pary wodnej w przepływie w funkcji s. 3. pokazano przebieg procentowej pary wodnej w domenie. Okres Do dalszej analizy wybrano okres zmian

ys. 3. pionowymi liniami.

Przebieg procentowego udziału pary wodnej w przepły-

rzebieg zjawiska w trakcie okres zmian, t – czas od . Można zauważyć pojawienie się w pobliżu krawędzi natarcia, jej oraz stopniowe zapadanie się.

jedna struktura ulega zapadnięciu w pobliżu krawędzi spływu, kolejna zaczyna rosnąć przy ys. 5. przedstawiono rozkład

(4)

objętościowego udziału pary w punktach monitorujących rozłożonych wzdłuż profilu. Objętościowy udział pary wodnej został uśredniony w czasie jednego okresu zmian.

Na osi odciętych umieszczono bezwymiarową odległość x/c, czyli stosunek współrzędnej odciętej do

cięciwy profilu. Najwyższa wartość udziału objętościowego pary wodnej zostaje osiągnięta dla x/c wynoszącego 0.2. Następnie udział pary wodnej maleje, jedynie przy krawędzi spływu, na końcu profilu można zauważyć jego wzrost.

Rys. 4. Rozkład udziału fazy ciekłej w czasie wybranego okresu kawitacji warstwowej, symulacja 2D

objętościowego udziału pary w punktach monitorujących rozłożonych wzdłuż profilu. Objętościowy udział pary został uśredniony w czasie jednego okresu zmian.

Na osi odciętych umieszczono bezwymiarową odległość nej odciętej do długości Najwyższa wartość udziału objętościowego pary wodnej zostaje osiągnięta dla x/c szącego 0.2. Następnie udział pary wodnej maleje, jedynie przy krawędzi spływu, na końcu profilu można

Rozkład udziału fazy ciekłej w czasie wybranego okresu

Rys. 5. Rozkład objętościowego udziału pary wodnej 2D

4.2. SYMULACJA 3D

Obliczenia powtórzono dla siatki uwzględniającej ścian boczną kanału przepływowego.

w przypadku symulacji 2D, monitorowano objętość pary wodnej powstającej w

rys. 6. przedstawiono przebieg procentowego udziału objętości pary wodnej dla obliczeń

przypadku 2D, po odcinku początkowy ustalił się na poziomie 3 ms.

Rys. 6. Przebieg procentowego udziału pary wodnej w wie w funkcji czasu, obliczenia 3D

Przebieg zjawiska w trakcie wybranego okresu zmian (zaznaczonego na rys. 6. pionowymi liniami) przedstawiony na rys. 7. Płaszczyzna cięcia domeny obliczeniowej w tym przypadku znajdowała się na płaszczyźnie z zadanym warunkiem

Od poczatku okresu następuje systematyczny wzrost obłoku kawitacyjnego, jednak dynamika zmian jest odmienna niż w przypadku obliczeń 2D.

kawitacyjne są bardziej rozległe się przy górnej powierzchni profilu.

zanikanie chmury, jednak proces ten zachodzi wolniej niż w przypadku obliczeń 2D (dla t/T = 0.6 dla obliczeń 2D struktury są już jedynie szczątkowe, dla obliczeń 3D – wciąż rozległe). Na rys. 8. przedstawiono

objętościowy udziału pary w punktach monitorujących.

Rozkład objętościowego udziału pary wodnej dla symulacji 3D w zakresie długości cięciwy profilu od 0.1 do 0.4 znacznie odbiega od przebiegu tego parametru dla obliczeń 2D. W przypadku obliczeń z

objętościowego udziału pary wodnej, symulacja

Obliczenia powtórzono dla siatki uwzględniającej ścianę kanału przepływowego. Podobnie jak przypadku symulacji 2D, monitorowano całkowitą objętość pary wodnej powstającej w przepływie. Na rzedstawiono przebieg procentowego udziału pary wodnej dla obliczeń 3D. Podobnie jak w początkowym, okres zmian

Przebieg procentowego udziału pary wodnej w przepły-

trakcie wybranego okresu zmian ys. 6. pionowymi liniami) został Płaszczyzna cięcia domeny obliczeniowej w tym przypadku znajdowała się na iem brzegowym symetrii.

astępuje systematyczny wzrost obłoku kawitacyjnego, jednak dynamika zmian jest nieco przypadku obliczeń 2D. Chmury oraz dłużej utrzymują profilu. Występuje stopniowe zanikanie chmury, jednak proces ten zachodzi wolniej niż w przypadku obliczeń 2D (dla t/T = 0.6 dla obliczeń 2D struktury są już jedynie szczątkowe, dla obliczeń 3D ys. 8. przedstawiono rozkład objętościowy udziału pary w punktach monitorujących.

Rozkład objętościowego udziału pary wodnej dla symulacji 3D w zakresie długości cięciwy profilu od 0.1 do 0.4 znacznie odbiega od przebiegu tego parametru dla obliczeń 2D. W przypadku obliczeń z uwzglęnieniem

(5)

ściany bocznej wartość tego parametru załamuje się dla współrzędniej x/c wynoszącej 0.1. Oznacza to wystąpienie krótszej (zajmującej mniej długości profilu) struktury przy krawędzi natarcia profilu łopatki.

W obydwu przypadkach obliczeń 2D i zaobserwowano wzrost udziału pary wodnej w przepływie w okolicach krawędzi spływu łopatki, czyli dla x/c od 0.8 do 1.

Rys. 7. Rozkład udziału fazy ciekłej w czasie wybranego okresu kawitacji warstwowej, symulacja 3D

ściany bocznej wartość tego parametru załamuje się dla współrzędniej x/c wynoszącej 0.1. Oznacza to wystąpienie krótszej (zajmującej mniej długości profilu) struktury przy krawędzi natarcia profilu łopatki.

obydwu przypadkach obliczeń 2D i 3D zaobserwowano wzrost udziału pary wodnej w przepływie w okolicach krawędzi spływu łopatki, czyli

Rozkład udziału fazy ciekłej w czasie wybranego okresu

Rys. 8. Rozkład objętościowego udziału pary wodnej wzdłuż profilu (uśredniony w jednym okresie)

4.3. SYMULACJA 3D –

W PŁASZCZYZNACH XY

Na rys. 9. przedstawiono przebieg jednego okresu kawitacji dla trzech różnych płaszczyzn XY. Pozwoliło to na przeprowadzenie analizy otrzymywanych wyników wzdłuż kierunku prostopadłego do kierunku przepływu cieczy (kierunek z).

W związku z tym przecięto domenę obliczeniową trzema równoległymi płaszczyznami XY. Pierwsza (A) znajduje się przy płaszczyźnie, na której zadano warunek ściana (z = 5 mm). Druga (B) została ulokowana w połowie domeny obliczeniowej (z = 17.5 mm). Trzecia ( znajduje się na płaszczyźnie na której zadano warunek symetrii (z=35 mm). Na rys. 9. przedstawiono przebieg jednego okresu kawitacji dla trzech różnych przekrojów wzdłużnych.

Rozkład objętościowego udziału pary wodnej wzdłuż profilu (uśredniony w jednym okresie)

– ROZKŁADY W PŁASZCZYZNACH XY

9. przedstawiono przebieg jednego okresu kawitacji dla trzech różnych płaszczyzn XY. Pozwoliło ie analizy otrzymywanych wyników wzdłuż kierunku prostopadłego do kierunku przepływu

W związku z tym przecięto domenę obliczeniową trzema równoległymi płaszczyznami XY. Pierwsza (A) znajduje się przy płaszczyźnie, na której zadano warunek ściana (z = 5 mm). Druga (B) została ulokowana w połowie domeny obliczeniowej (z = 17.5 mm). Trzecia (C) znajduje się na płaszczyźnie na której zadano warunek 9. przedstawiono przebieg jednego okresu kawitacji dla trzech różnych przekrojów

(6)

Rys. 9. Rozkład udziału fazy ciekłej w czasie wybranego okresu kawitacji warstwowej Analiza obrazów struktur parowych przedstawionych na

rys. 9. dostarcza informacji na temat wpływu ściany bocznej na ich wielkość i kształt. W okolicy ściany struktury te są bardziej jednolite. W czasie t/T = 0 kształt i zakres chmury parowej na każdym z przekrojów jest niemal identyczny. Dla t/T = 0.1 można zauważyć, że dla przekroju najbardziej oddalonego od ściany obszaru pary wodnej staje się silnie

Tendencja ta pogłębia się, dla t/T = 0.3 kształt struktury parowej dla przekroju C jest znacznie bardziej poszarpany niż dla przekrojów A i B. Od t/T = 0.4 można także zauważyć różnicę między przekro

B. Kształt struktury przy ścianie jest dalej jednolity, podczas gdy dla przekroju B można zauważyć podział chmury parowej w okolicach krawędzi spływu. W tym czasie w przekroju C obłok parowy dzieli się na kilka

iekłej w czasie wybranego okresu kawitacji warstwowej, symulacja 3D Analiza obrazów struktur parowych przedstawionych na

dostarcza informacji na temat wpływu ściany na ich wielkość i kształt. W okolicy ściany W czasie t/T = 0 kształt i zakres chmury parowej na każdym z przekrojów jest niemal identyczny. Dla t/T = 0.1 można zauważyć, że dla przekroju najbardziej oddalonego od ściany zarys silnie nieregularny.

Tendencja ta pogłębia się, dla t/T = 0.3 kształt struktury parowej dla przekroju C jest znacznie bardziej poszarpany niż dla przekrojów A i B. Od t/T = 0.4 można także zauważyć różnicę między przekrojami A i B. Kształt struktury przy ścianie jest dalej jednolity, podczas gdy dla przekroju B można zauważyć podział chmury parowej w okolicach krawędzi spływu. W tym czasie w przekroju C obłok parowy dzieli się na kilka

mniejszych o różnej wielkości. W poł

t/T = 0.5 w przekroju A pojawia się podział chmury kawitacyjnej przy krawędzi spływu, natomiast w przekrojach B i C następuje stopniowe zanikanie struktur kawitacyjnych. Dla t/T = 0.8 w przekroju A chmura kawitacyjna wciąż rozciąga się na

profilu, podczas gdy w przekrojach B i C pozostają już jedynie jej niewielkie części. Dla t/T = 1 obrazy struktur kawitacyjnych dla wszystkich trzech przekrojów są podobne.

W celu lepszego zobrazowania zjawiska zorientowano widok izometrycznie, dzięki czemu możliwa jest obserwacja całej chmury kawitacyjnej na górnej ścianie profilu łopatki. Na rys. 10. przedstawiono widok na całą łopatkę z zaznaczonymi powierzchniami o stałej wartości udziału objętościowego pary. Dla kolumny

W połowie okresu, dla 0.5 w przekroju A pojawia się podział chmury zy krawędzi spływu, natomiast przekrojach B i C następuje stopniowe zanikanie struktur kawitacyjnych. Dla t/T = 0.8 w przekroju A chmura kawitacyjna wciąż rozciąga się na całą długość profilu, podczas gdy w przekrojach B i C pozostają już jedynie jej niewielkie części. Dla t/T = 1 obrazy struktur kawitacyjnych dla wszystkich trzech przekrojów

W celu lepszego zobrazowania zjawiska zorientowano , dzięki czemu możliwa jest obserwacja całej chmury kawitacyjnej na górnej ścianie przedstawiono widok na całą łopatkę z zaznaczonymi powierzchniami o stałej wartości udziału objętościowego pary. Dla kolumny D udział

(7)

objętościowy pary wynosi 0.5, dla kolumny

0.9. Przepływ odbywa się w prawo, usytuowanie warunków brzegowych jest takie samo jak na

Rys. 10. Rozkład objętościowego udziału pary wodnej (uśre niony w jednym okresie) (D α=0.5, E α=0.9)

Zarówno dla udziału objętościowego pary wynoszącego 0.5, jak i 0.9 można zauważyć, że im dalej od ściany bocznej tym izopowierzchnia staje się bardziej nieregularna i poszarpana. Od połowy okresu (t/T = 0.5) tendencja ta nie jest już tak wyraźna ze względu na to, że przy ścianie bocznej, za krawędzią wlotu łopatki następuje oderwanie chmury. Na rys. 11.

uśrednione w czasie wartości udziału objętościowego pary wodnej dla przekroju A, B i C. Można zauwa duże rozbieżności między wartością tego param

poszczególnych przekrojów. Jedynie w połowie długości profilu (x/c = 0.45) wartości udziału pary były niemal identyczne. Najwyższa wartość objętościowego udziału pary wodnej, równa 0.87, została zanotowana dla przekroju A dla współrzędnej x/c wynos

wy pary wynosi 0.5, dla kolumny E natomiast odbywa się w prawo, usytuowanie warunków brzegowych jest takie samo jak na rys.2.

Rozkład objętościowego udziału pary wodnej (uśred-

=0.9)

udziału objętościowego pary wynoszącego jak i 0.9 można zauważyć, że im dalej od ściany bocznej tym izopowierzchnia staje się bardziej

poszarpana. Od połowy okresu (t/T = 0.5) tendencja ta nie jest już tak wyraźna ze względu na

, za krawędzią wlotu łopatki 11. przedstawiono uśrednione w czasie wartości udziału objętościowego

Można zauważyć duże rozbieżności między wartością tego parametru dla

Jedynie w połowie długości u (x/c = 0.45) wartości udziału pary były niemal

Najwyższa wartość objętościowego udziału tała zanotowana dla przekroju A dla współrzędnej x/c wynoszącej 0.1.

Rys.11. Rozkład objętościowego udziału pary wodnej (uśre niony w jednym okresie)

4.4 PRZEBIEG CIŚNIENIA W CZASIE W PUNKTACH

MONITORUJĄCYCH

Analiza kształtu i dynamiki zmian struktur kawitacyjnych dostarcza wielu informacji o charakterze przepływu kawitacyjnego. W celu pogłębienia tej analizy sporządzono wykresy zmian wartości ciśnienia

w punktach monitorujących ulokowanych wokół prof dla przekrojów A, B i C oznaczonych na

Wybrano cztery współrzędne x/c. Wykres 1 sporządzono dla x/c wynoszącego 0.01, wykres 2 dla x/c = 0.2, wykres 3 dla x/c = 0.57 oraz wykres 4 dla x/c

Na wykresach przedstawiony jest przebieg ciśnien ciągu jednego okresu zmian przedstawionego na za pomocą pionowych linii. W okolicy krawędzi natarcia (rys. 11., wykres 1) nie zanotowano dużych ciśnienia. Podczas obliczeń ciśnienie we wszystkich trzech punktach znajdujących się

tej samej wartości współrzędnej x/c nie przekroczyło wartości 40 kPa. Krzywe przebiegu ciśnienia dla przekrojów A, B i C są do siebie bardzo zbliżone.

miarę oddalania się od krawędzi natarcia przebieg ciśnienia staje się nieregularny,

gwałtowne zmiany (rys.11., wykres 2).

ciśnień dla różnych przekrojów wzdłużnych zaczynają od siebie odbiegać, skoki ciśnienia wystepują w innym momencie cyklu. Najwięcej gwałtownych zmian ciśnienia można zaobserwować dla przekroju C.

ciśnienie wzrasta do 100 kPa. Dla współrzędnej x/c = 0.57 (rys.11., wykres 3) odnotowano po jednym gwałtownym wzroście ciśnienia na każdym

wzdłużnym. Dla przekroju C ciśnienie wzrosło aż do 300 kPa. W okolicy krawędzi spływu dochodzi do największych zmian wartości ciśnienia (

wykres 4).

Rozkład objętościowego udziału pary wodnej (uśred-

PRZEBIEG CIŚNIENIA W CZASIE

MONITORUJĄCYCH

Analiza kształtu i dynamiki zmian struktur kawitacyjnych dostarcza wielu informacji o charakterze przepływu kawitacyjnego. W celu pogłębienia tej analizy wartości ciśnienia w czasie w punktach monitorujących ulokowanych wokół profilu dla przekrojów A, B i C oznaczonych na rys. 9.

Wybrano cztery współrzędne x/c. Wykres 1 sporządzono 0.01, wykres 2 dla x/c = 0.2,

= 0.57 oraz wykres 4 dla x/c = 0.89.

Na wykresach przedstawiony jest przebieg ciśnienia w przedstawionego na rys. 6.

W okolicy krawędzi natarcia ykres 1) nie zanotowano dużych zmian ciśnienie we wszystkich trzech punktach znajdujących się na jednej krawędzi o tej samej wartości współrzędnej x/c nie przekroczyło 0 kPa. Krzywe przebiegu ciśnienia dla są do siebie bardzo zbliżone. W miarę oddalania się od krawędzi natarcia przebieg nieregularny, można zauważyć jego ykres 2). Krzywe przebiegu ciśnień dla różnych przekrojów wzdłużnych zaczynają od ciśnienia wystepują w innym Najwięcej gwałtownych zmian ciśnienia a przekroju C. Maksymalnie Dla współrzędnej x/c = odnotowano po jednym łtownym wzroście ciśnienia na każdym przekroju wzdłużnym. Dla przekroju C ciśnienie wzrosło aż do 300 W okolicy krawędzi spływu dochodzi do największych zmian wartości ciśnienia (rys. 11,

(8)

Rys. 11. Przebieg ciśnienia w punktach monitorujących dla różnych przekrojów wzdłużnych; 1) x/c=0.01; 2) x/c=0.2; 3) x/c=0.57;

4) x/c=0.89

Dla wszystkich trzech przekrojów krzywe ciśnienia mają nieregularny przebieg. Zmiany wartości ciśnienia dochodzą do ponad 350 kPa pod koniec okresu zmian.

Wysokie, gwałtowne zmiany wartości ciśnienia są jedną z cech charakterystycznych przepływu, w którym występuje zjawisko kawitacji. Wynikają one z nagłych zmian parametrów przy przejściu z fazy ciekłej na gazową i odwrotnie. Występowanie przy krawędzi spływu profilu przepływu silnie turbulentnego, o dużej dynamice zmian parametrów zostało także odnotowane podczas eksperymentów [16]. W związku z

stwierdzić, że przyjęty model kawitacji dobrze odwzorowuje główne, charakterystyczne

przepływu kawitacyjnego.

5. PODSUMOWANIE

W ramach prac przedstawionych w

przeprowadzono symulację przepływu niestacjonarnego kawitacją zachodzącą na górnej, podciśnieniowej powierzchni profilu łopatki. W symulacji użyto modelu kawitacji Schnerr & Sauer. Obliczenia przeprowadzono w środowisku OpenFOAM, przy użyciu sol interPhaseChangeFoam.

W pierwszym etapie prowadzonych prac przeprowadzono symulację 2D przepływu.

że po obu stronach profilu łopatki założono warunek symetria. Symulacja 2D pozwoliła

cykliczny wzrost i zapadanie się struktur k

Przebieg ciśnienia w punktach monitorujących dla różnych przekrojów wzdłużnych; 1) x/c=0.01; 2) x/c=0.2; 3) x/c=0.57;

Dla wszystkich trzech przekrojów krzywe ciśnienia mają wartości ciśnienia dochodzą do ponad 350 kPa pod koniec okresu zmian.

Wysokie, gwałtowne zmiany wartości ciśnienia są jedną tycznych przepływu, w którym Wynikają one z nagłych zmian parametrów przy przejściu z fazy ciekłej na Występowanie przy krawędzi spływu profilu przepływu silnie turbulentnego, o dużej ian parametrów zostało także odnotowane W związku z tym można stwierdzić, że przyjęty model kawitacji dobrze , charakterystyczne cechy

W ramach prac przedstawionych w artykule niestacjonarnego z , podciśnieniowej W symulacji użyto modelu Obliczenia przeprowadzono , przy użyciu solvera prowadzonych prac przeprowadzono symulację 2D przepływu. Oznacza to, że po obu stronach profilu łopatki założono warunek symetria. Symulacja 2D pozwoliła zaobserwować ktur kawitacyjnych

oraz zindentyfikować obszary na profilu

w których występuje najwyższy udział pary wodnej Następnie przeprowadzono obliczenia 3D.

obliczeniowa została poszerzona, po jednej stronie profilu założono warunek ściany. W wyniku obliczeń uzyskano obrazy trójwymiarowe chmur pary wodnej.

Podobnie jak w przypadku obliczeń 2D, zaobserowano cykliczny wzrost i zapadanie się struktur parow o częstotliwości zmian takiej samej jak dla obliczeń dwuwymiarowych. Zaobserwowano znaczne różnice w dynamice zmian chmur kawitacyjnych

prostopadłym do kierunku przepływu cieczy.

ściany obłoki kawitacyjne były wyraźnie ba stabilne oraz jednorodne. Najbardziej nieregularny kształt zauważono w przekroju znajdującym się na płaszczyźnie o zadanym warunku brzegowym symetria, czyli najdalej oddalonym od ścia

obserwowane przy ścianie najprawdopodobniej wynikają z mniejszych turbulencji w tym obszarze.

można także zaobserwować bardziej stabilne chmury kawitacyjne w przypadku obliczeń 3D niż w

obliczeń 2D na ścianie o warunku brzegowym symetria.

Obserwację realizowaną na płaskich przekrojach poszerzono o widok trójwymiarowy powierzchni o stałej wartości udziału objętościowego pary.

także przebieg zmian ciśnienia podczas wybranego cyklu wzrostu i zapadania się struktury parowej

monitorujących rozłożonych wzdłuż profilu w różnych przekrojach wzdłużnych

Przebieg ciśnienia w punktach monitorujących dla różnych przekrojów wzdłużnych; 1) x/c=0.01; 2) x/c=0.2; 3) x/c=0.57;

na profilu badanej łopatki, najwyższy udział pary wodnej.

Następnie przeprowadzono obliczenia 3D. Domena obliczeniowa została poszerzona, po jednej stronie profilu założono warunek ściany. W wyniku obliczeń uzyskano obrazy trójwymiarowe chmur pary wodnej.

Podobnie jak w przypadku obliczeń 2D, zaobserowano adanie się struktur parowych, częstotliwości zmian takiej samej jak dla obliczeń Zaobserwowano znaczne różnice chmur kawitacyjnych w kierunku prostopadłym do kierunku przepływu cieczy. W pobliżu ściany obłoki kawitacyjne były wyraźnie bardziej Najbardziej nieregularny ju znajdującym się na płaszczyźnie o zadanym warunku brzegowym symetria, any. Stabilne struktury obserwowane przy ścianie najprawdopodobniej wynikają z mniejszych turbulencji w tym obszarze. Z tego powodu rwować bardziej stabilne chmury awitacyjne w przypadku obliczeń 3D niż w przypadku 2D na ścianie o warunku brzegowym symetria.

Obserwację realizowaną na płaskich przekrojach poszerzono o widok trójwymiarowy powierzchni o stałej wartości udziału objętościowego pary. Przeanalizowano ciśnienia podczas wybranego cyklu wzrostu i zapadania się struktury parowej w punktach monitorujących rozłożonych wzdłuż profilu w trzech różnych przekrojach wzdłużnych. Zauważono, że im

(9)

dalej od krawędzi spływu łopatki, tym krzywe przebie wartości ciśnienia w czasie dla poszczególnych przekrojów wzdłużnych zaczynają od siebie odbie Zanotowano gwałtowne skoki ciśnienia, dochodzące nawet do 350 kPa. Najbardziej nieregularny przebieg wartości ciśnienia w czasie odnotowano przy krawędzi spływu łopatki. Jest to obszar najbardziej narażony na skutki implozji pęcherzyków pary.

Pracę zrealizowano w ramach projektu finansowanego przez Narodowe Centrum Nauki oraz w ramach badań statusowych BK 287/RIE5/2016/502.

Literatura

1. Arndt R.: Some remarks on hydrofoil cavitation. “ 2. Brennen C. E.: Cavitation and bubble

3. Brennen C. E.: Hydrodynamics of Pumps.

4. Escaler X., Egusquiza E., Farhat M., Avellon F., Coussirat M.

“Mechanical Systems and Signal Processing”

5. Gosset A., Casas V.D. and Pena F.L.:

around a 2D hydrofoil. In: Fifth OpenFOAM Workshop, 6. Huang B., Wang G.: Partially averaged Navier

“Journal of Hydrodynamics” 2011, Vol.23(1) 7. Kubota A., Kato H., Yamaguchi H.

tion on a hydrofoil section, “Journal of Fluid Mechanics”

8. Kumar P., Saini R. P.: Study of cavitation in hydro Reviews” 2010, Vol. 14(1), p. 374-383

9. Li D. Q., Grekula M., Lindell P.: A modified SST k cavitation on 2D and 3D hydrofoils. In:

CAV2009, Ann Arbor, Michigan, USA

10. Li D. Q., Grekula, M., Lindell P.: Towards numerical prediction of unsteady sheet cavitation on hydro 9th International Conference on Hydrodynamics

supplement: p. 741-746.

11. Li Z., Pourquie M., Van Terwisga T.:

In: 9th International Conference on Hydrodynamics Vol.22(5), supplement: p. 770-777.

12. Puffary B.: Numerical modelling of cavitation RTO-EN-AVT-143, 2006, paper 3, p. 3

13. Roohi E., Zahiri A. P. and Pasandideh

hydrofoil using VOF method and LES turbulence. “ 6488.

14. Saito Y., Takami R., Nakamori I.,

a hydrofoil section, “Computational Mechanics” 2007, Vol. 40, p. 85

15. Yuan W., Schnerr G. H., Sauer J.: Modeling and computation of unsteady cavitating flows in injection nozzles.

“Mecanique and Industries Journal” 2001, Vol. 2, p. 383 16. Wang G., Senocak I., Shyy W., Ikohago T. and Cao S.

gress in Aerospace Sciences” 2001, Vol.

Artykuł dostępny na podstawie licencji Creative Commons Uznanie autorstwa 3.0 Polska.

http://creativecommons.org/licenses/by/3.0/pl

tym krzywe przebiegu dla poszczególnych rojów wzdłużnych zaczynają od siebie odbiegać.

Zanotowano gwałtowne skoki ciśnienia, dochodzące nawet do 350 kPa. Najbardziej nieregularny przebieg wartości ciśnienia w czasie odnotowano przy krawędzi Jest to obszar najbardziej narażony na

Na podstawie wyników przeprowadzonych analiz można stwierdzić, że symulacja 3D dostarczyła cennych informacji na temat rozkładu parametrów w kierunku prostopadłym do kierunku przepływu, co z kolei pozwala na lepszą identyfikację obszarów najbardziej narażonych na gwałtowne zmiany parametrów przepływu.

Uproszczenie do przypadku dwuwymiarowego pozwala znacznie zmniejszyć czas obliczeń, nie oddaje jednak w pełni złożoności zjawiska.

Pracę zrealizowano w ramach projektu finansowanego przez Narodowe Centrum Nauki UMO- BK 287/RIE5/2016/502.

remarks on hydrofoil cavitation. “Journal of Hydrodynamics” 2012, Vol(24)3 ubble dynamics. Oxford University Press 2014.

.: Hydrodynamics of Pumps. Oxford University Press 2011.

Egusquiza E., Farhat M., Avellon F., Coussirat M.: Detection of cavitation in hydraulic turbines

“Mechanical Systems and Signal Processing” 2006, Vol. 20, p. 983-1007

V.D. and Pena F.L.: Evaluation of the cavitating foam solver for low Mach OpenFOAM Workshop, 2010, Gothenburg, Sweden

Partially averaged Navier-Stokes methods for time-dependent turbulent cavitating flows.

Vol.23(1), p. 26-33.

Yamaguchi H.: A new modelling of cavitating flows: a numerical study of unsteady cavit , “Journal of Fluid Mechanics” 1992, Vol. 240, p. 59-96

Study of cavitation in hydro – turbines – A review. “Renewable and Sustainable Energy 383

A modified SST k- turbulence model to predict the steady and

itation on 2D and 3D hydrofoils. In: Proceedings of the 7th International Symposium on Cavitation Arbor, Michigan, USA, 2009, Vol. 1, p. 39-51.

Towards numerical prediction of unsteady sheet cavitation on hydro 9th International Conference on Hydrodynamics. Shanghai, China, 2010, “Journal of Hydrodynamics

e M., Van Terwisga T.: A numerical study of steady and unsteady cavitat

9th International Conference on Hydrodynamics. Shanghai, China, 2010, “Journal of Hydrodynamics

merical modelling of cavitation, design and analysis of high speed pumps 143, 2006, paper 3, p. 3-1 – 3.54..

and Pasandideh-Fard M.: Numerical simulation of cavitation around a

LES turbulence. “Applied Mathematical Modelling” 2013, Vol. 37 (9), p. 6469

Ikohagi T.: Numerical analysis of unsteady behavior of cloud cavitation , “Computational Mechanics” 2007, Vol. 40, p. 85-96

: Modeling and computation of unsteady cavitating flows in injection nozzles.

“Mecanique and Industries Journal” 2001, Vol. 2, p. 383-394.

Senocak I., Shyy W., Ikohago T. and Cao S.: Dynamics of attached turbulent cavitating flows Vol. 37, p. 551- 581.

http://creativecommons.org/licenses/by/3.0/pl

przeprowadzonych analiz można stwierdzić, że symulacja 3D dostarczyła cennych informacji na temat rozkładu parametrów w kierunku prostopadłym do kierunku przepływu, co z kolei pozwala na lepszą identyfikację obszarów najbardziej narażonych gwałtowne zmiany parametrów przepływu.

Uproszczenie do przypadku dwuwymiarowego pozwala obliczeń, nie oddaje jednak

-2016/21/B/ST8/01164

Vol(24)3, p. 305-314.

Detection of cavitation in hydraulic turbines,

ow Mach number flow

dent turbulent cavitating flows.

of cavitating flows: a numerical study of unsteady cavita-

“Renewable and Sustainable Energy

teady and unsteady sheet International Symposium on Cavitation

Towards numerical prediction of unsteady sheet cavitation on hydrofoils. In:

Journal of Hydrodynamics” Vol.22(5),

A numerical study of steady and unsteady cavitation on a 2D hydrofoil.

Journal of Hydrodynamics”

umps. Educational Notes

avitation around a two – dimensional Applied Mathematical Modelling” 2013, Vol. 37 (9), p. 6469-

Numerical analysis of unsteady behavior of cloud cavitation on

: Modeling and computation of unsteady cavitating flows in injection nozzles.

Dynamics of attached turbulent cavitating flows. “Pro-