OPTYMALIZACJA ALGORYTMÓW WYZNACZANIA RUCHU CIECZY LEPKIEJ METODĄ SZTUCZNEJ ŚCIŚLIWOŚCI

(1)

36, s. 187-192, Gliwice 2008

OPTYMALIZACJA ALGORYTMÓW WYZNACZANIA

RUCHU CIECZY LEPKIEJ METODĄ SZTUCZNEJ ŚCIŚLIWOŚCI

Z

BIGNIEW

K

OSMA

, B

OGDAN

N

OGA Instytut Mechaniki Stosowanej, Politechnika Radomska e-mail: zbigniew.kosma@pr.radom.pl, b.noga@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ą sztucznej ściśliwości, opisywaną równaniami w zmiennych fizycznych: składowe prędkości, ciśnienie. Nowe algorytmy numeryczne zaadap- towano do rozwiązywania zagadnień modelowych, ze względu na moŜliwość po- równywania wyników własnych obliczeń numerycznych z wynikami prezentowanymi w publikacjach i rezultatami badań eksperymentalnych.

1. WSTĘP

Zasadniczą ideą metody sztucznej ściśliwości jest przyłączenie pochodnej ciśnienia wzglę- dem czasu do równania ciągłości, co umoŜliwia sprzęŜenie ciśnienia z prędkością. Równa- niami wyjściowymi do wyznaczania przepływów cieczy lepkiej metodą sztucznej ściśliwości jest układ równań róŜniczkowych cząstkowych utworzony ze zmodyfikowanego równania ciągłości i równań Naviera-Stokesa. Przy rozwiązywaniu sformułowanego zagadnienia po- czątkowo-brzegowego dla ciśnienia i składowych prędkości wszystkie pochodne względem zmiennych przestrzennych aproksymowano przy wykorzystaniu klasycznych ilorazów róŜni- cowych drugiego rzędu dokładności na równomiernych siatkach obliczeniowych. Zachowując czas jako zmienną niezaleŜną ciągłą, uzyskano zagadnienie początkowe dla układu równań róŜniczkowych zwyczajnych dla nieznanych wartości tych funkcji w kaŜdym węźle we- wnętrznym wygenerowanej siatki obliczeniowej. Zagadnienie to rozwiązywano metodą Ga- lerkina-Rungego-Kutty trzeciego rzędu. 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 efektywne, uzyskano kilkukrotne skrócenie czasu obliczeń w porównaniu z czasami wyznaczania rozwiązań tych zagadnień za pomocą pakietu Fluent.

2. METODA SZTUCZNEJ ŚCIŚLIWOŚ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]:

(2)











∂

∂ + ∂

∂

− ∂

∂ = + ∂

∂

∂ =

∂

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.

W metodzie sztucznej ściśliwości sprzęŜenie pola ciśnienia z polem prędkości następuje w wyniku zastąpienia równania ciągłości dla cieczy nieściśliwej (1a) zlinearyzowanym rów- naniem ciągłości dla gazu [2]

1 0

~ =

∂ + ∂

∂

i i

x V t

p

β (2)

Optymalne wartości parametru relaksacyjnego β są zwykle wyznaczane niezaleŜnie dla kaŜdego problemu po wykonaniu szeregu eksperymentów numerycznych. W wielu przypad- kach zadawalające wyniki daje przyjęcie stałej wartości tego parametru w całym obszarze przepływu.

3. ZAGADNIENIA OBLICZENIOWE

Opracowane algorytmy numeryczne przystosowano do symulacji numerycznej ruchu cieczy 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 ka- nału warunki brzegowe dla prędkości wyraŜają nieprzenikalność i brak poślizgu, na jego wlo- cie przyjęto znikanie gradientu ciśnienia, na wylocie: p=1.

(3)

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) - (1b) 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 (3) w których U =

[

p,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.

Otrzymane w ten sposób algorytmy obliczeniowe są algorytmami uniwersalnymi, prostymi koncepcyjnie i łatwymi do realizacji numerycznej, pozwalającymi na ominięcie konieczności linearyzacji członów nieliniowych. Ich wadą jest sztywność układów równań róŜniczkowych zwyczajnych oraz moŜliwość występowania niestabilności numerycznych.

W celu przetestowania efektywności i dokładności zastosowanych algorytmów numerycznych podzielono je na dwie grupy. Pierwszą grupę tworzą algorytmy, w których pochodne względem zmiennych przestrzennych aproksymowano róŜnicami skończonymi drugiego rzę- du dokładności (MRS), w algorytmach grupy drugiej pochodne względem zmiennych przestrzennych aproksymowano za pomocą kompaktowych schematów róŜnicowych szóstego rzędu dokładności (KSR) - tabela 1.

Zagadnienia początkowe dla układów równań róŜniczkowych zwyczajnych (3) całkowano:

1) jednokrokową metodą predyktor-korektor Adamsa-Moultona (A-M),

2) jedno- i dwukrokowymi metodami predyktor-korektor opartymi na wykorzystaniu wzo- rów metody wstecznego róŜniczkowania (WR 1, WR 2),

3) metodami Rungego-Kutty: Eulera pierwszego stopnia (R-K 1), Heuna drugiego stopnia (R-K 2), Heuna trzeciego stopnia (R-K 3), Galerkina trzeciego stopnia (R-K-G), zmodyfiko- waną piątego stopnia (R-K 5).

Tabela 1. Algorytmy obliczeniowe

Grupy algoryt-

mów Aproksymacja pochodnych Całkowanie zagadnienia początkowego

MRS Klasyczne ilorazy róŜnicowe drugiego

rzędu dokładności A-M, WR 1, WR 2, R-K 1,

R-K 2, R-K 3, R-K-G, R-K 5 KSR Kompaktowe schematy róŜnicowe szóste-

go rzędu dokładności

Po wykonaniu testowych symulacji numerycznych okazało się, Ŝe szybsze są algorytmy na- leŜące do pierwszej grupy - tabela 2. Dokładniejsze wyniki otrzymywano natomiast w przy- padku zastosowania aproksymacji pochodnych względem zmiennych przestrzennych kompak- towymi ilorazami róŜnicowymi szóstego rzędu dokładności. W obydwu grupach algorytmów

(4)

najefektywniejsze okazało się całkowanie zagadnień początkowych (3) metodą Rungego- Kutty-Galerkina trzeciego stopnia.

Tabela 2. Efektywność algorytmów

Wybór wariantu obliczeń Szybkość algorytmów Dokładność obliczeń

Aproksymacja pochodnych przestrzennych

Klasyczne ilorazy róŜnicowe drugiego rzędu dokładności (MRS)

Kompaktowe schematy róŜnicowe szóstego rzędu dokładności (KSR)

Metoda całkowania zagadnienia początkowego

Metoda Rungego-Kutty-Galerkina trzeciego stopnia (R-K-G)

5. WYNIKI SYMULACJI NUMERYCZNYCH

Symulacje numeryczne zostały wykonane dla danych wyszczególnionych w tabeli 3, okre- ślających rozmiary siatek oraz wartości liczb Reynoldsa i kroków czasowych. Dokładność ustalania się modułów pochodnych układu równań (2) przyjęto równą 1⋅10⁻¹⁴.

Tabela 3. Dane do obliczeń numerycznych

Parametr obliczeń Kwadratowe zagłębienie Sześcienne zagłębienie Kanał z uskokiem ścianki

Rozmiar siatki 150 × 150 100 × 100 × 100 1000 × 50

Liczba Reynoldsa 10 ÷ 7500 10 ÷ 1000 100 ÷ 1000

Krok czasowy 1⋅10⁻² ÷ 1⋅10⁻³ 1⋅10⁻³ 1⋅10⁻³

W metodzie sztucznej ściśliwości bardzo waŜnym parametrem jest współczynnik β, jego optymalna wartość została określona na podstawie szeregu eksperymentów numerycznych.

Okazało się, Ŝe wartość parametru β nie ma istotnego wpływu na dokładności obliczeń, ma natomiast bardzo duŜy wpływ na szybkość uzyskiwania rozwiązań; metoda Rungego-Kutty- Galerkina jest najefektywniejsza dla parametru β = 1.

Na rys. 2 ÷ 6 zostały przedstawione wyniki najwaŜniejszych obliczeń numerycznych dla rozwiązywanych zagadnień oraz ich porównania z wynikami prezentowanymi w publikacjach.

a) b) c)

Rys. 2. Kwadratowe zagłębienie - wyniki obliczeń dla Re = 5000 na siatce 150 × 150:

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

(5)

a) b) c)

Rys. 3. Kwadratowe zagłębienie - wyniki obliczeń dla Re = 7500 na siatce 150 × 150:

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)

6. PODSUMOWANIE

We wszystkich rozwaŜanych zagadnieniach uzyskano poprawne wyniki symulacji numerycznych 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 została więc w pełni potwierdzona.

Algorytmy te są nieskomplikowane, efektywniejsze od algorytmów wykorzystywanych w pakietach komercyjnych. MoŜna je łatwo zmodyfikować do wyznaczania ruchu cieczy lepkiej w obszarach ograniczonych nieregularnymi liniami brzegowymi oraz dla przepływów turbulentnych.

LITERATURA

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

2. Kosma Z.: Symulacja numeryczna ruchu cieczy lepkiej metodą sztucznej ściśliwości.

Monografie. Radom: WPR, 2007.

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.” 1982, 48, p. 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” 2005, 48, p.

747-774.

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

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

OPTIMIZED ALGORITHMS FOR THE CALCULATIONS OF VISCOUS INCOMPRESSIBLE FLOWS

USING THE ARTIFICIAL COMPRESSIBILITY METHOD

Summary. The introduced pseudo-time derivative of pressure directly couples the pressure with the velocity and changes the mathematical character of the con- tinuity equation from elliptic to hyperbolic. A standard method of lines approach is applied in this contribution. The system of partial differential equation is discre- tized in space by central, second-order finite-difference schemes on uniform grids with the same mesh sizes in each direction, the time-variable preserved conti- nuous. The initial-boundary value problem for this system of equations is then re- duced to an initial value problem for a system of ordinary differential equations, and the unknown values of pressure and velocity components in each inner knot of uniform meshes are computed using the Galerkin-Runge-Kutta method of third order. 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 algorithms proved to be very effective for the demanded time of calculations.