WYZNACZANIE RUCHU CIECZY LEPKIEJ METODĄ SZTUCZNEJ ŚCIŚLIWOŚCI NA SIATKACH NAKŁADAJĄCYCH SIĘ

MODELOWANIE INŻYNIERSKIE ISNN 1896-771X 32, s. 267-272, Gliwice 2006

WYZNACZANIE RUCHU CIECZY LEPKIEJ METODĄ SZTUCZNEJ ŚCIŚLIWOŚCI NA SIATKACH NAKŁADAJĄCYCH SIĘ

Z

K

B

N

P

M

Instytut Mechaniki Stosowanej, Politechnika Radomska

Streszczenie. Przedstawione zostały własne algorytmy numeryczne, przeznaczone do wyznaczania dwuwymiarowego, stacjonarnego ruchu cieczy lepkiej metodą sztucznej ściśliwości w obszarach kanałów prostoliniowych, podwójnie zagiętych.

Kanały wlotowe i wylotowe mają stałe długości, różne są długości kanałów pośrednich. 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 trzech równomiernych, nakładających się siatkach 30 L × 30 w kanałach o długościach L = 2 ÷ 9 dla Re ≤ 200.

1. METODA SZTUCZNEJ ŚCIŚLIWOŚCI

Wyznaczanie płaskiego 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

1 =0



∂ +∂

∂ + ∂

∂

∂

y x u

p v

β

τ (1) oraz równań Naviera-Stokesa, zapisanych w postaci diwergentnej:

Re . 1 ) (

) (

Re , 1 ) ( ) (

2 2

2 2

v v v

v

v

∇

∂+ = +∂

∂ +∂

∂

∂

∇

∂ = +∂

∂+ +∂

∂

∂

r r

y p x

u t

y u u x

p u t u

(2)

Układ równań (1) – (2) zawiera trzy funkcje niewiadome: ciśnienie p(x,y,τ) oraz dwie składowe prędkości u(x,y,t) i v(x,y,t), ∇r²

oznacza operator laplasjanu, Re jest liczbą Reynoldsa. Występujący w równaniu ciągłości (1) 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ładzie równań opisujących ruch cieczy lepkiej (1) – (2) pochodnej ciśnienia względem pseudo-czasu τ (przyjęto τ =t), a jej zasadniczą wadą – ograniczenie stosowania tylko do wyznaczania ruchu stacjonarnego.

2. SFORMUŁOWANIE ZAGADNIENIA

Rozważymy zagadnienie wyznaczania ruchu cieczy lepkiej w obszarach kanałów prostoliniowych, podwójnie zagiętych z dwoma wklęsłymi narożami [2 – 5]. Przyjęto (rys. 1), że wszystkie kanały mają jednakową wysokość – H =1, a ich długości są następujące: kanał górny – L₁=2, kanał dolny – L₃ =9, kanał pośredni – L₂ =(2÷4)H.

y

x 0 1

−1 1

−2

−3

2 3 4 5 6 7 8 9 10

H H H

L₁

L2

L₃

Rys. 1. Prostoliniowy kanał podwójnie zagięty

Dla określenia warunków brzegowych dla układu równań (1) – (2) w zakresie przepływu laminarnego przyjmiemy paraboliczne rozkłady składowej prędkości u na wlocie do kanałów i wylocie z kanałów oraz ponadto założymy, że wydatek przepływającego strumienia cieczy jest równy Q=1. Stąd mamy

. 0 ),

1 (

6 ₀

0 = − ₌ =

= x

x y y

u v

Warunki brzegowe dla składowych prędkości na ściankach kanałów wyrażają nieprzenikalność i brak poślizgu.

Ciśnienie spełnia na ściankach kanałów warunek brzegowy wynikający z równań Naviera- Stokesa

Re ,

1 2

Ω

Ω ∂

∂

∇

∂∂ = n V

r r r

p (3)

w przekrojach wlotowych założono ∂p ∂x=0, w przekrojach wylotowych – p=1.

3. OBLICZENIA NUMERYCZNE

Obszary kanałów podwójnie zagiętych (rys. 1) podzielono na trzy nakładające się prostokątne kanały (rys. 2) o długościach – L = 2 ÷ 9, w których wygenerowano trzy równomierne, nakładające się siatki 30 L × 30.

Przy rozwiązywaniu zagadnień początkowo-brzegowych dla układów równań (1) – (2) w każdym kanale (prostokącie) wszystkie pochodne względem zmiennych przestrzennych aproksymowano przy wykorzystaniu klasycznych ilorazów różnicowych drugiego rzędu dokładności [6]. 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 F U

t

U = (4)

gdzie ^{U =}

[

]

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

Wartości ciśnienia na ściankach kanałów i na wlocie określano z warunków brzegowych (3) oraz ze znikania gradientu ciśnienia, po uprzedniej aproksymacji pochodnych ciśnienia względem normalnej do brzegów obszarów klasycznymi, jednostronnymi ilorazami różnicowymi O(h²).

I

II

III

Rys. 2. Podział kanału podwójnie zagiętego na trzy nakładające się kanały prostokątne

Zagadnienia początkowe dla układów równań (4) rozwiązywano oddzielnie w każdym kanale prostokątnym z przekazywaniem wartości ciśnienia i składowych prędkości na nakładających się granicach w dwóch sąsiadujących ze sobą kanałach. Obliczenia wykonano dla Re ≤ 200. Pierwsza iteracja z warunkami początkowymi określającymi stan spoczynku została pokazana na rys. 3, do ustalenia się pól prędkości i ciśnienia niezbędne było wykonanie około 100 iteracji globalnych i 50 000 ÷ 900 000 iteracji lokalnych na każdej siatce. 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 –

, 10

1 ⁸

2

⋅ −

=

ε dokładność ustalania się modułów pochodnych układów równań (4) – .

10 1 ¹⁰

1

⋅ −

=

ε Rezultaty obliczeń numerycznych dla L₂ =2 i L₂ =4 zostały przedstawione na rys. 4 – 5.

I

II

III

Rys. 3. Warunki początkowe w nakładających się obszarach kanału podwójnie zagiętego

Rys. 4. Obrazy funkcji prądu w kanale podwójnie zagiętym dla L₂ =2

Rys. 5. Obrazy funkcji prądu w kanale podwójnie zagiętym dla L₂ =4

4. PODSUMOWANIE

Stwierdzono dużą skuteczność i szybkość działania opracowanych algorytmów numerycznych wyznaczania ruchu cieczy lepkiej w obszarach o geometrycznie skomplikowanych kształtach, polegających na zastosowaniu siatek nakładających się oraz sprowadzaniu zagadnień początkowo-brzegowych dla układów różniczkowych cząstkowych do zagadnień początkowych dla układów równań różniczkowych zwyczajnych (metoda prostych). Algorytmy te są proste koncepcyjnie i mogą być łatwo przystosowane do wyznaczania laminarnego ruchu cieczy lepkiej w różnych obszarach płaskich.

Uzyskano zadowalającą zgodność położenia stref recyrkulacyjnych w rozważanych obszarach prostoliniowych kanałów podwójnie zagiętych, w porównaniu z wynikami obliczeń numerycznych innych autorów [3 – 5]. Założenie stacjonarności przepływów okazało się istotną trudnością, uniemożliwiającą wykonanie obliczeń dla wyższych liczb Reynoldsa.

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. Tulapurkara E.G., Lakshmana Gowda B.H., Balachandran N.: Laminar flow through slots, J. Fluid Mech., 190, 1988, s. 179-200.

3. Hwang Y.-H.: Arbitrary domain velocity analyses for the incompressible Navier-Stokes equations, J. Comp. Phys., 110, 1994, s. 134-149.

4. Bathe K.-J., Zhang H.: A flow-condition-based interpolation finite element procedure for incompressible fluid flows, Comp. Struct., 80, 2002, s. 1267-1277.

5. Ramšak M., Škerget L.: A subdomain boundary element method for high-Reynolds lami- nar flow using stream function-vorticity formulation, Int. J. Numer. Meth. Fluids, 46, 2004, s. 815-847.

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

7. 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 OVERLAPPING GRIDS USING THE ARTIFICIAL

COMPRESSIBILITY METHOD

Summary. The artificial compressibility method is designed for computation of stationary viscous incompressible flows in double bent channels. The spatial deriva- tives and the boundary conditions are discretized by means of the classical second- order finite-difference schemes on three overlapping uniform grids, while preserv- ing the time-variable continuos. The resulted system of ordinary differential equa- tions has been integrated using the two-step backward-differentiation predictor- corrector method. Calculations have been made for Re ≤ 200 on three uniform 30 L × 30 grids for the channels lengths L = 2 ÷ 9.