OPTYMALIZACJA ALGORYTMÓW WYZNACZANIA RUCHU CIECZY LEPKIEJ METODĄ DEKOMPOZYCJI POLA PRĘDKOŚCI

36, s. 181-186, Gliwice 2008

OPTYMALIZACJA ALGORYTMÓW WYZNACZANIA RUCHU CIECZY LEPKIEJ METODĄ DEKOMPOZYCJI POLA PRĘDKOŚCI

Z

BIGNIEW

K

OSMA

, P

RZEMYSŁAW

M

OTYL Instytut Mechaniki Stosowanej, Politechnika Radomska e-mail: zbigniew.kosma@pr.radom.pl, p.motyl@pr.radom.pl

Streszczenie. Celem pracy było poszukiwanie i optymalizacja efektywnych algo- rytmów obliczeniowych wyznaczania ruchu cieczy lepkiej w obszarach płaskich i przestrzennych - konkurencyjnych do kodów komercyjnych. Wyznaczano ruch cieczy lepkiej metodą dekompozycji pola prędkości, opisywaną równaniami w zmiennych fizycznych: składowe prędkości, ciśnienie. Nowe algorytmy nume- ryczne zaadaptowano do rozwiązywania zagadnień modelowych, ze względu na moŜliwość porównywania wyników własnych obliczeń numerycznych z wynikami prezentowanymi w publikacjach i rezultatami badań eksperymentalnych.

1. WSTĘP

Zgodnie z podstawowym postulatem metody dekompozycji pola prędkości, tzn. rozszcze- pienia pól prędkości i ciśnienia, obliczenia są wykonywane w dwu etapach. W pierwszym etapie obliczeń w kaŜdym kroku czasowym wyznaczane jest pomocnicze pole prędkości, w drugim - rozwiązywane jest zagadnienie Neumanna dla ciśnienia obliczeniowego i dokonywana jest korekcja składowych prędkości. Przy rozwiązywaniu zagadnienia począt- kowo-brzegowego dla układu równań róŜniczkowych cząstkowych w pierwszym etapie obli- czeń stosowano metodę prostych, polegającą na zachowaniu czasu jako zmiennej niezaleŜnej ciągłej i sprowadzaniu go do zagadnienia początkowego dla układu równań róŜniczkowych zwyczajnych dla nieznanych wartości składowych prędkości pomocniczej w kaŜdym węźle wewnętrznym wygenerowanej, równomiernej siatki. Wszystkie pochodne względem zmien- nych przestrzennych aproksymowano klasycznymi ilorazami róŜnicowymi drugiego rzędu do- kładności, zagadnienie początkowe dla układu równań róŜniczkowych zwyczajnych całkowa- no metodą Heuna drugiego rzędu. Zagadnienie Neumanna dla ciśnienia obliczeniowego roz- wiązywano metodą róŜnic skończonych z wykorzystaniem metody Gaussa-Seidela. Stwier- dzono przy tym istotny fakt, Ŝe ciśnienie obliczeniowe nie musi być wyznaczane do pełnej zbieŜności - ponadto dodatkowe skrócenie czasu obliczeń uzyskano przy ich realizacji na dwóch siatkach. Wykonano szereg obliczeń testowych dla zagadnień ruchu cieczy lepkiej w zagłębieniach z jedną poruszającą się ścianką: kwadratowym i sześciennym oraz w płaskim kanale z uskokiem jednej ścianki. Opracowane algorytmy obliczeń okazały się bardzo efek- tywne, uzyskano kilkukrotne skrócenie czasu obliczeń w porównaniu z czasami wyznaczania rozwiązań tych zagadnień za pomocą pakietu Fluent.

2. METODA DEKOMPOZYCJI POLA PRĘDKOŚCI

Po wprowadzeniu zmiennych bezwymiarowych i pominięciu pola sił masowych jednost- kowych układ równań opisujący niestacjonarny ruch cieczy lepkiej w postaci zachowawczej ma postać [1]:











∂

∂

∂ + ∂

∂

− ∂

∂ = + ∂

∂

∂

∂ =

∂

Re , 1 , 0

i j j i

j i j i

j j

x V x x

V p x V t V

x V

(1)

w którym V₁,V₂,V₃ są składowymi prędkości w kierunkach osi kartezjańskiego układu współ- rzędnych x₁,x₂,x₃, p jest ciśnieniem, Re - liczbą Reynoldsa.

Ogólna idea zastosowanego algorytmu metody dekompozycji pola prędkości - będącego modyfikacją schematu rozszczepienia, zaproponowanego pierwotnie w pracy [2] przez Husera i Biringena - polega na wykonywaniu obliczeń na kaŜdej nowej warstwie czasowej w dwóch etapach. W pierwszym etapie obliczeń w przedziale czasu od tⁿ do t~

rozwiązywane jest za- gadnienie początkowo-brzegowe dla pomocniczego pola prędkości ~,

Vi określonego uprosz- czonymi równaniami Naviera-Stokesa dla znanego gradientu ciśnienia

~, Re

~ 1

~

~ ~

i j j i

n j

i j

i V

x x x

V p x V t V

∂

∂

∂ + ∂

∂

−∂

∂ = + ∂

∂

∂ (2)

przy załoŜeniu, Ŝe na granicach obszaru ∂ΩΩΩΩ i w chwili czasowej t =tⁿ pole prędkości po- mocniczej jest identyczne z polem prędkości fizycznej:

n.

i n i i

i V V V

V = _∂ =

∂

, ~

~

Ω Ω Ω Ω Ω

ΩΩ

Ω (3)

W drugim etapie obliczeń dla kaŜdego kroku czasowego ∆t=tⁿ⁺¹−tⁿ w przedziale czasu od

~t

do tⁿ⁺¹ korygowane są wartości składowe prędkości pomocniczej z zaleŜności

,

~

~ 2

~ ¹

1 



 





∂

−∂

∂

− ∂

=

+ +

i n

i n i

n

i x

p x

p V t

V ∆

otrzymanych po scałkowaniu równań sprzęgających pola prędkości fizycznej i pomocniczej z gradientem ciśnienia obliczeniowego, po uprzednim rozwiązaniu zagadnienia Neumanna dla ciśnienia obliczeniowego p~ na warstwie czasowej tⁿ⁺¹

j n j

j j n

j

j x

V p t

x p x

x

x ∂

+ ∂

∂

∂

∂

= ∂

∂

∂

∂

∂ +

2 ~

~

~ ¹

∆ ⁽⁴⁾

z jednorodnymi warunkami brzegowymi dla gradientu ciśnienia obliczeniowego.

183

3. ZAGADNIENIA OBLICZENIOWE

Opracowane algorytmy numeryczne przystosowano do symulacji numerycznej ruchu cie- czy lepkiej w zagłębieniach z jedną poruszającą się ścianką: kwadratowym (rys. 1a) i sze- ściennym (rys. 1b) oraz w prostoliniowym kanale z uskokiem jednej ścianki (rys. 1c). Zagad- nienia te są często rozwiązywane w celu testowania efektywności i dokładności róŜnych algo- rytmów obliczeniowych [3 - 6].

a) b) c)

Rys. 1. Zagadnienia obliczeniowe: a) kwadratowe zagłębienie, b) sześcienne zagłębienie, c) kanał z uskokiem jednej ścianki

Ścianki górne zagłębień poruszają się ze stałą prędkością, równoległą do osi x.

W przekrojach: wlotowym i wylotowym kanału przyjęto paraboliczne rozkłady prędkości, wynikające z równości wydatku

(

Q=0.5

)

przepływającego strumienia cieczy lepkiej. Na ściankach kanału warunki brzegowe dla prędkości wyraŜają nieprzenikalność i brak poślizgu, na jego wlocie przyjęto znikanie gradientu ciśnienia, na wylocie: p=1.

4. ALGORYTMY OBLICZENIOWE

Wszystkie symulacje numeryczne w obszarach zagłębień z jedną poruszającą się ścianką oraz w kanale z uskokiem ścianki wykonano na równomiernych siatkach obliczeniowych o kwadratowych oczkach.

Do rozwiązywania zagadnienia (2) - (3) w pierwszym etapie obliczeń zastosowano metodę prostych, polegającą na jego sprowadzeniu do układów równań róŜniczkowych zwyczajnych przy zachowaniu czasu jako zmiennej niezaleŜnej ciągłej. Uzyskano w ten sposób układy równań róŜniczkowych zwyczajnych postaci

( )

^,

d

d F U

U =t (5)

w których U

[ ]

V~i ^T

= jest wektorem zmiennych niezaleŜnych, obliczanym w kaŜdym węźle wewnętrznym siatek, F - operatorem róŜniczkowania względem zmiennych przestrzennych.

Po wykonaniu szeregu testów numerycznych najefektywniejszymi okazały się:

• aproksymacja pochodnych względem zmiennych przestrzennych klasycznymi, trzypunk- towymi ilorazami róŜnicowymi drugiego rzędu dokładności,

• całkowanie zagadnień początkowych dla układów (5) metodą Heuna drugiego rzędu.

Najefektywniejsze okazało się równieŜ rozwiązywanie zagadnienia Neumanna dla równa- nia Poissona (4) metodą róŜnic skończonych z aproksymacją pochodnych względem zmien-

nych przestrzennych klasycznymi róŜnicowymi drugiego rzędu dokładności. Otrzymane ukła- dy algebraicznych równań liniowych rozwiązywano metodą Gaussa-Seidela. Przy wykonywa- niu eksperymentów numerycznych stwierdzono, Ŝe wyznaczanie ciśnienia obliczeniowego w drugim etapie obliczeń nie musi być realizowane do pełnej zbieŜności na kolejnych war- stwach czasowych (w całym procesie obliczeń), co istotnie skraca czas obliczeń przy zacho- waniu Ŝądanej dokładności. Dodatkowe skrócenie czasu obliczeń uzyskano po zastosowaniu do rozwiązywania układów algebraicznych równań liniowych metody dwóch siatek.

5. WYNIKI SYMULACJI NUMERYCZNYCH

Symulacje numeryczne zostały wykonane dla danych wyszczególnionych w tabeli 1 - okre- ślających rozmiary siatek oraz wartości liczb Reynoldsa i kroków czasowych. Przyjęto jedna- kowe dokładności ustalania się modułów pochodnych układu równań (2) oraz rozwiązywania zagadnienia Neumanna dla równania (4) metodą róŜnic skończonych - równą 1⋅10⁻¹².

Tabela 1. Dane do obliczeń numerycznych

Parametr obliczeń Kwadratowe zagłębienie Sześcienne zagłębienie Kanał z uskokiem ścianki Rozmiar siatki 50 × 50 ÷ 200 × 200 50 × 50 ÷ 100 × 100 1000 × 50 Liczba Reynoldsa 100 ÷ 10 000 100 ÷ 1000 100 ÷ 800 Krok czasowy 1⋅10⁻² ÷ 1⋅10⁻³ 1⋅10⁻³ 1⋅10⁻² ÷ 1⋅10⁻³

Obliczenia wykonywano albo na jednej siatce z ograniczeniem liczby iteracji (30 iteracji) albo teŜ na dwóch siatkach: wykonywano najpierw 5 iteracji na siatce gęstej, 20 iteracji na siatce rzadkiej i - na koniec - 5 iteracji na siatce gęstej (łącznie kaŜdorazowo 30 iteracji).

Stwierdzono ponadto, Ŝe przyjęta liczba siatek nie ma wpływu na dokładność obliczeń, ma natomiast istotny wpływ na zmniejszenie czasu obliczeń (ok. 20 %). Zastosowanie algorytmu obliczeń na dwóch siatkach przyspiesza zbieŜność procesu wyznaczania ciśnienia obliczenio- wego oraz zmniejsza oscylacje błędu obliczeń składowych prędkości (rys. 2).

a) b)

Rys. 2. Kwadratowe zagłębienie - błędy obliczeń dla Re 100 na siatce 50 × 50:

a) składowych prędkości, b) ciśnienia obliczeniowego

185

Na rys. 3 ÷ 6 przedstawiono wyniki najwaŜniejszych obliczeń numerycznych dla rozwią- zywanych zagadnień oraz ich porównania z wynikami prezentowanymi w publikacjach.

a) b) c)

Rys. 3. Kwadratowe zagłębienie - wyniki obliczeń dla Re = 10 000 na siatce 200 × 200:

a) linie prądu; rozkłady składowych prędkości: b) u - na linii x = 0.5, c) v - na linii y = 0.5

a) b) c)

Rys. 4. Sześcienne zagłębienie - wyniki obliczeń dla Re = 1000 na siatce 100 × 100 × 100:

a) linie prądu w płaszczyźnie z = 0.5; rozkłady składowych prędkości: b) u - na linii x = 0.5, c) v - na linii y = 0.5

Rys. 5. Kanał z uskokiem jednej ścianki - linie prądu dla Re = 800, siatka 1000 × 50

a) b)

Rys. 6. Kanał z uskokiem jednej ścianki: Re = 800, siatka 1000 × 50. Rozkłady składowej prędkości u dla na liniach: a) x = 7, b) x = 15

6. PODSUMOWANIE

We wszystkich rozwaŜanych zagadnieniach uzyskano poprawne wyniki symulacji nume- rycznych w szerokim zakresie liczb Reynoldsa. Stwierdzono dobrą zgodność wartości skła- dowych prędkości na osiach symetrii zagłębień oraz w wybranych przekrojach kanału z uskokiem ścianki z wynikami analogicznych obliczeń - prezentowanymi w publikacjach.

Skuteczność zastosowanych algorytmów: metody prostych w pierwszym etapie obliczeń oraz wyznaczania ciśnienia obliczeniowego w drugim etapie obliczeń została więc w pełni po- twierdzona. Algorytmy te są nieskomplikowane, efektywniejsze od algorytmów wykorzysty- wanych w pakietach komercyjnych. MoŜna je łatwo zmodyfikować do wyznaczania ruchu cieczy lepkiej w obszarach ograniczonych nieregularnymi liniami brzegowymi, stosując struk- tury danych typu lista, graf.

LITERATURA

1. Kosma Z.: Podstawy mechaniki płynów. Radom: WPR, 2007.

2. Huser A.D., Biringen S.: Calculation of two-dimensional shear-driven cavity flows at high Reynolds numbers. Int. J. Numer. Meth. Fluids, 14, 1992, s. 1087-1109.

3. Ghia U., Ghia K.N., Shin C.T.: High-Re solutions for incompressible flow using the Na- vier-Stokes equations and a multigrid method. J. Comp. Phys., 48, 1982, s. 387-411.

4. Erturk E., Corke T.C., Gökçöl C.: Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers. Int. J. Numer. Meth. Fluids, 48, 2005, s. 747-774.

5. Shu C., L. Wang L., Chew Y.T.: Numerical computation of three-dimensional incom- pressible Navier-Stokes equations in primitive variable form by DQ method. Int. J. Nu- mer. Meth. Fluids, 43, 2003, s. 345-368.

6. Gartling D.K.: A test problem for outflow boundary conditions-flow over a backward- facing step. Int. J. Numer. Meth. Fluids, 11, 1990, s. 953-967.

OPTIMIZED ALGORITHMS FOR THE CALCULATIONS OF VISCOUS INCOMPRESSIBLE FLOWS

USING A VELOCITY CORRECTION METHOD

Summary. New version of an algorithm for unsteady motion of a viscous incom- pressible fluid is proposed. The algorithm can be considered as a modification of the projection scheme originally developed by Huser and Biringen. In the first step of calculation the method of lines is adopted, while for solving the Neumann prob- lem for the computational pressure the finite difference method is applied. The numerical experiments showed that the computational pressure had not to be de- termined accurately, moreover two grids method has been found to work very well in a reasonable accelerating the rate of convergence. Test calculations for laminar flows in the square and cubic cavities with one moving wall and the backward- facing step have been performed. The proposed algorithm proved to be very effec- tive for the demanded time of calculations.