WYZNACZANIE RUCHU CIECZY LEPKIEJ METODĄ SZTUCZNEJ ŚCIŚLIWOŚCI NA SIATKACH NIERÓWNOMIERNYCH

(1)

32, s. 261-266, Gliwice 2006

WYZNACZANIE RUCHU CIECZY LEPKIEJ METODĄ SZTUCZNEJ ŚCIŚLIWOŚCI NA SIATKACH NIERÓWNOMIERNYCH

Z

BIGNIEW

K

OSMA

B

OGDAN

N

OGA

P

RZEMYSŁAW

M

OTYL

Instytut Mechaniki Stosowanej, Politechnika Radomska

Streszczenie. Przedstawione zostały własne algorytmy numeryczne, przeznaczone do wyznaczania dwu- i trójwymiarowego, stacjonarnego ruchu cieczy lepkiej metodą sztucznej ściśliwości w obszarach zagłębień: kwadratowego i sześciennego z jedną poruszającą się ścianką. Zagadnienia te rozwiązywano metodą prostych, sprowadzając zagadnienia początkowo-brzegowe dla układów równań różniczkowych cząstkowych do zagadnień początkowych dla układów równań różniczkowych zwyczajnych. Obliczenia wykonano na nierównomiernych siatkach 30 × 30 oraz 30 × 30 × 30 dla Re ≤ 1000.

1. METODA SZTUCZNEJ ŚCIŚLIWOŚCI

Wyznaczanie ruchu cieczy lepkiej metodą sztucznej ściśliwości [1] polega na rozwiązywaniu zagadnienia początkowo-brzegowego dla układu równań różniczkowych cząstkowych – utworzonego ze zmodyfikowanego równania ciągłości oraz równań Naviera- Stokesa. W ogólnym przypadku ruchu trójwymiarowego przybiera on następującą postać:

Re , 1

2 2 2 2 2 2

G W F E

W D





∂ + ∂

∂

= ∂

∂ +∂

∂

z y z x

y x

t (1)

gdzie:

, ] , , ,

[p u v w ^T

=

W E =[u β, p+u²,uv,uw]^T, ,

] , ,

,

[v uv p+v² vw ^T

= β

F G=[w β,uw,vw,p+w^{2 T}] ,

[

⁰^,¹^,¹^,¹

]

^,

=diag D

a w przypadku ruchu płaskiego jest zmodyfikowanym układem równań (1):

Re , 1

2 2 2 2

F W E

W D





∂ + ∂

∂

= ∂

∂ +∂

∂

y y x

x

t (2)

(2)

gdzie:

, ] , , [p u v ^T

=

W E=[u β, p+u²,uv]^T, ,

] ,

,

[v uv p+v^{2 T}

= β

F ^D⁼^diag

[

⁰^,¹^,¹

]

^.

W układach równań (1) – (2) niewiadomymi funkcjami są: ciśnienie p oraz składowe prędkości [u,v,w], ∇r²

oznacza operator laplasjanu, Re jest liczbą Reynoldsa. Parametr relaksacyjny β wynika z założenia zmiennej gęstości cieczy lepkiej i zastąpienia równania ciągłości dla cieczy nieściśliwej zlinearyzowanym równaniem ciągłości dla gazu.

Zaletą metody sztucznej ściśliwości w rozważanej postaci jest fakt pojawienia się w układach równań opisujących ruch cieczy lepkiej (1) – (2) pochodnej ciśnienia względem pseudo-czasu, a jej zasadniczą wadą – ograniczenie stosowania tylko do wyznaczania ruchu stacjonarnego.

2. SFORMUŁOWANIE ZAGADNIEŃ

Rozważymy zagadnienia wyznaczania ruchu cieczy lepkiej w zagłębieniach: kwadratowym (rys. 1a) i sześciennym (rys. 1b) z jedną poruszającą się ścianką, będące modelami przepływu wody w zagłębieniu w ziemi przy silnie wiejącym wietrze. Zagadnienia te są bardzo często rozwiązywane w celu testowania efektywności i dokładności różnych algorytmów obliczeniowych, przy czym wyniki prezentowane w pracach [2 – 3] są uznawane za rozwiązania dokładne i z nimi są porównywane wyniki innych obliczeń numerycznych.

a) b)

Rys. 1. Zagłębienia z jedną poruszającą się ścianką: a) kwadratowe, b) sześcienne

Przyjęto, że boki zagłębień są równoległe do osi lub płaszczyzn układu współrzędnych, przy czym ścianki y=1 poruszają się ze stałą prędkością – równoległą do osi x, a warunki brzegowe na pozostałych granicach obszarów wyrażają nieprzenikalność i brak poślizgu na ściankach. Ciśnienie spełnia na wszystkich ściankach zagłębień warunek brzegowy wynikający

(3)

Re ,

1 2

Ω

Ω ∂

∂

∇

∂∂ = n V

r r r

p (3)

przy ustalonej wartości ciśnienia w środku poruszających się ścianek.

3. OPIS PRZEBIEGU OBLICZEŃ

Obliczenia zostały wykonane na siatkach nierównomiernych, otrzymanych po transformacji obszaru pierwotnego do obszaru kanonicznego we współrzędnych ξ , lub η ξ,η,ζ :

, 1 , ,

0 ≤ ξ η ζ ≤ zdefiniowanych zależnościami:

, , ,

2 3

ζ ζ ζ

η η η

ξ ξ ξ

c b a z

c b a y

c b a x

+ +

=

+ +

=

+ +

=

ze współczynnikami a,b,c zawierającymi parametry p i ₀ p : _N

. ,

1

, 2

0 0

0

p c p a b

p p

a _N

=

−

=

− +

=

Przy rozwiązywaniu zagadnień początkowo-brzegowych dla układów równań (1) – (2) wszystkie pochodne względem zmiennych przestrzennych aproksymowano przy wykorzystaniu klasycznych ilorazów różnicowych drugiego rzędu dokładności [4]. Przy zachowaniu czasu t jako zmiennej niezależnej ciągłej uzyskano zagadnienia początkowe dla układów równań różniczkowych zwyczajnych postaci:

( )

^, d

d F U

t

U = (4)

gdzie ^U ⁼

[

^p^,^u^,^v

]

^T^lub ^U ⁼

[

^p^,^u^,^v^,^w

]

^T są wektorami zmiennych zależnych, a F jest różniczkowym operatorem przestrzennym.

Zagadnienia początkowe dla układów równań różniczkowych zwyczajnych (4) rozwiązywano jednokrokową metodą predyktor-korektor wstecznego różniczkowania [5].

Wartości ciśnienia na ściankach zagłębień określano z warunków brzegowych (3), po uprzedniej aproksymacji pochodnych ciśnienia względem normalnej do brzegów obszarów klasycznymi, jednostronnymi ilorazami różnicowymi O(h²).

(4)

4. WYNIKI OBLICZEŃ NUMERYCZNYCH

Obliczenia przepływów cieczy lepkiej w zagłębieniach z jedną poruszającą się ścianką zostały wykonane na nierównomiernych siatkach 30 × 30 oraz 30 × 30 × 30 dla parametrów:

. 25 . 0 ,

5 .

0 ₀

0 = p_N = p = p_N =

p

Algorytmy obliczeń testowano dla dwóch największych liczb Reynoldsa: Re = 400 i Re = 1000, dla których możliwe było uzyskanie rozwiązań stacjonarnych na siatkach równomiernych. Przyjęto następujące dane sterujące przebiegiem obliczeń: parametr relaksacyjny β =0.8, krok czasowy – ∆t=1⋅10⁻³, dokładność korekcji metody wstecznego różniczkowania – ε₂ =1⋅10⁻⁸, dokładność ustalania się modułów pochodnych układów równań (4) – ε₁ =1⋅10⁻¹⁰. Warunkami początkowymi były rozwiązania z siatek równomiernych. Rozwiązania ustalone uzyskano po wykonaniu około 30 000 ÷ 100 000 kroków czasowych. Na rys. 2 ÷ 5 zamieszczono porównania wyznaczonych składowych prędkości: u na linii x = 0.5, v na linii y = 0.5 – odpowiednio, w płaszczyznach: przepływu oraz z = 0.5 – z wynikami analogicznych obliczeń zamieszczonych w publikacjach [2 – 3].

Otrzymane pola prędkości dla Re = 1000 na siatkach N =30 z parametrami p₀ = p_N =0.25 zostały zamieszczone na rys. 6.

Rys. 2. Porównania rozkładów składowych prędkości w kwadratowym zagłębieniu dla Re = 400: u na linii x = 0.5, v na linii y = 0.5; siatka 30 × 30, p₀ = p_N =0.5

Rys. 3. Porównania rozkładów składowych prędkości w kwadratowym zagłębieniu dla Re = 1000: u na linii x = 0.5, v na linii y = 0.5; siatka 30 × 30, p₀ = p_N =0.25

(5)

Rys. 4. Porównania składowych prędkości w sześciennym zagłębieniu w płaszczyźnie z = 0.5:

u na linii x = 0.5, v na linii y = 0.5; Re = 400, siatka 30 × 30 × 30, p₀ = pN =0.5

Rys. 5. Porównania składowych prędkości w sześciennym zagłębieniu w płaszczyźnie z = 0.5:

u na linii x = 0.5, v na linii y = 0.5; Re = 1000, siatka 30 × 30 × 30, p₀ = p_N =0.25

a) b)

Rys. 6. Pola prędkości w zagłębieniach: a) kwadratowym – siatka 30 × 30, b) sześciennym w płaszczyźnie z = 0.5 – siatka 30 × 30 × 30; Re =1000, p₀ = pN =0.25

(6)

5. PODSUMOWANIE

Skuteczność opracowanych algorytmów numerycznych, polegających na sprowadzaniu zagadnień początkowo-brzegowych dla układów równań różniczkowych cząstkowych do zagadnień początkowych dla układów równań różniczkowych zwyczajnych (metoda prostych), została w pełni potwierdzona przy wyznaczaniu zarówno dwu-, jak i trójwymiarowego ruchu cieczy lepkiej. Algorytmy te są bardzo proste koncepcyjnie i pozwalają na szybkie wyznaczanie laminarnych przepływów cieczy lepkiej w różnych obszarach płaskich i przestrzennych.

Stwierdzono bardzo duży wpływ zastosowanego sposobu zagęszczania siatek obliczeniowych przy ściankach zagłębień, gdzie występują duże gradienty składowych prędkości, na dokładność obliczeń. Na rzadkich siatkach nierównomiernych okazało się możliwe uzyskiwanie rozwiązań z bardzo dużą dokładnością, wymagającą przyjęcia dwu- lub trzykrotnie gęstszych siatek równomiernych.

LITERATURA

1. Hirsch Ch.: Numerical computational of internal and external flows. Vol 2: Computational methods for inviscid and viscous flows. New York: John Wiley & Sons,1990.

2. 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.

3. Shu C., Wang L., Chew Y.T.: Numerical computation of three-dimensional incompressible Navier-Stokes equations in primitive variable form by DQ method. Int. J. Numer. Meth.

Fluids, 43, 2003, s. 345-368.

4. Kosma Z.: Metody numeryczne dla zastosowań inżynierskich. Radom: WPR, 2004.

5. Hairer E., Norsett S.P., Wanner G.: Solving ordinary differential equations I. Nonstiff problems. Berlin: Springer-Verlag, 1993.

NUMERICAL SIMULATIONS OF VISCOUS FLUID MOTIONS ON NON-UNIFORM GRIDS USING THE ARTIFICIAL

COMPRESSIBILITY METHOD

Summary. The artificial compressibility method is designed for computation of stationary viscous incompressible flows in the square and cubic cavities. The spatial derivatives and the boundary conditions are discretized by means of the classical second-order finite-difference schemes on non-uniform grids, while preserving the time-variable continuos. The resulted system of ordinary differential equations has been integrated using the one-step backward-differentiation predictor-corrector method. Calculations have been made for Reynolds number values of 400 and 1000 on the 30 × 30 and 30 × 30 × 30 non-uniform grids.