Súrlódó folyadék stacionárius áramlása kétdimenzióban Koonin: Computational Physics^{∗} könyve alapján

Bartha Ferenc

e-mail: barthaf@physx.u-szeged.hu, http : www.jate.u-szeged.hu/˜ barthaf/

Szegedi Tudományegyetem, Elméleti Fizikai Tanszék (May 9, 2002)

I. A HIDRODINAMIKAI ALAPEGYENLETEK

Az áramló folyadékok jellemzésére a folyadék ρ (r, t) sűrűségét és a V (r, t) sebességmezőt használjuk. Az ezekre vonatkozó kontinuitási egyenlet

\frac{∂}{∂t}ρ + div ρV = 0  (1)

a tömeg megmaradását fejezi ki: a tér minden pontja körül a sűrűség az oda befolyó, illetve az onnan kifolyó anyagmennyiség miatt változik meg. A súrlódó folyadékok áramlását leíró másik egyenlet (Navier-Stokes)

ρ \{\frac{∂V}{∂t} + (V·grad)V\} = f − grad p + η · ΔV + (η + η′) grad div V  (2)

lényegében az impulzus megmaradását jelenti. Az egyenletben a p (r, t, T) nyomás gradiense szerepel. Ha a hőmérséklet változik, akkor a nyomást meghatározó állapotegyenletet is fel kellene írnunk. Most erre nem lesz szükségünk.

II. INKOMPRESSZIBILIS STACIONÁIUS ÁRAMLÁS KÉT DIMENZIÓBAN

Meg kívánjuk határozni, hogy viszkózus folyadék hogyan áramlik az útjába helyezett akadály mellett.

• Satcionárius áramlásra szorítkozunk, így az iménti egyenletekben minden idő szerinti derivált nulla.

• A folyadékot összenyomhatatlannak tekintjük, azaz a sűrűség a térben állandó.

• A z-koordinátában a problémát transzlációra invariánsan fogalmazzuk meg, úgy, hogy elegendő legyen az (x, y) síkmetszetben megoldani az egyenleteket.

• A folyadék hőmérsékletét térben és időben állandónak tekintjük, ezáltal a nyomás pusztán a megadott egyenletekből meghatározható.

Ilyen feltételek mellett a kontinuitási egyenletből

div V(x, y) = \frac{∂V_{x}}{∂x} + \frac{∂V_{y}}{∂y} = 0  (3)

míg a másikból

(V_{x} \frac{∂}{∂x} + V_{y} \frac{∂}{∂y} ) V_{x} = −\frac{1}{ρ} \frac{∂p}{∂x} + ν (\frac{∂^{2}}{∂x^{2}} + \frac{∂^{2}}{∂y^{2}} ) V_{x} és (V_{x} \frac{∂}{∂x} + V_{y} \frac{∂}{∂y} ) V_{y} = −\frac{1}{ρ} \frac{∂p}{∂y} + ν (\frac{∂^{2}}{∂x^{2}} + \frac{∂^{2}}{∂y^{2}} ) V_{y} (4)

ahol bevezettük a ν = η/ρ kinematikai viszkozitást. A három egyenlet a peremfeltételek tisztázása után ebben a fomájában megoldható a keresett V_{x}(x, y), V_{y}(x, y) és p(x, y) függvényekre.

Mégis szokás egyéb ’segédfüggvényeket bevezetni’. Ilyen a Ψ(x, y) áramlási függvény

∗The Benjamin/Cummings Publishing Company, 1986

V_{x} = \frac{∂}{∂y} Ψ(x, y) és V_{y} = −\frac{∂}{∂x}Ψ(x, y)  (5)

melynek létezése biztosított az adott differenciálegyenlet mellett. Szemléletes jelentése van, ugyanis szintvonalai az áramfonalakat adják, másképp mondva a sebesség mindenhol érinti a szintvonalakat: V·grad Ψ(x, y) = 0. A másik függvény a χ(x, y) örvényesség

χ(x, y) = \frac{∂V_{x}}{∂y} − \frac{∂V_{y}}{∂x} = − (rot V)_{z}  (6)

Az ezekre érvényes egyenletek

ΔΨ = χ ; Δχ = ν (\frac{∂χ}{∂x} \frac{∂Ψ}{∂y} − \frac{∂Ψ}{∂x} \frac{∂χ}{∂y}) ; Δp = 2ρ ((\frac{∂^{2}Ψ}{∂x^{2}}) (\frac{∂^{2}Ψ}{∂y^{2}}) − (\frac{∂^{2}Ψ}{∂x∂y} )^{2})  (7)

Örömünkre szolgál, hogy az egyenletek részben szétcsatolódtak. Meg kell előbb oldanunk az első két csatolt egyeneletet, majd az így kapott Ψ ismeretében a (ha egyáltalán érdekel bennünket a nyomás) a harmadikat.

III. A KONKRÉT PEREMÉRTÉKPROBLÉMA

Balról folyadék áramlik szabadon a tartományba V (x = −∞, y) = [V_{0}, 0]. A középen elhelyezett merev akadályt megkerülve a jobboldalon V (x = ∞, y) = [V_{0}, 0] sebességeloszlással távozik.

\begin{aligned}(−∞) & → \\ & → \\ & → & ↗ \\ & → \\ & → \\ V_{0} & → & → \\ & → \\ & → \\ & → & ↘ \\ & → \\ & → \end{aligned} − y = 0 − \begin{matrix}· & · & · \\ · & · & · \\ · & · & · \\ · & · & · \\ · & · & · \\ · & · & · \end{matrix}  \begin{matrix}→ (∞) \\ → \\ → \\ → \\ → \\ → V_{0} \\ → \\ → \\ → \\ → \\ → \end{matrix}  (8)

Az áramlási kép szimmetrikus lesz a felező (y = 0) egyenesre: V_{x}(x, y) = V_{x}(x, −y) és V_{y}(x, y) = −V_{y}(x, −y) következésképpen

V_{y}(x, 0) = 0 és \frac{∂}{∂y} V_{x}(x, y)|_{y=0} = 0  (9)

Az előbbiből

V_{y}(x, 0) = − \frac{∂}{∂x}Ψ(x, y)|_{y=0} = 0 ⟹ Ψ(x, 0) = konstans (speciálisan legyen Ψ(x, 0) = 0)  (10)

A konstanst megválaszthatjuk így, hiszen az áramlási függvénynek csak a deriváltjai érdekesek, bevezetése egy állandó erejéig határozatlanná teszi. Az akadályra merőlegesen folyadék nem folyhat, tehát az áramló folyadék sebessége párhuzamos az akadály oldalaival. Az akadályt így egyetlen áramvonal veszi körbe, ami egyúttal szintvonala az áramlási függvénynek. Ebből következik, hogy Ψ = 0 az akadály mentén. A szimmetria másik következménye, hogy

χ(x, 0) = \frac{∂V_{x}}{∂y} − \frac{∂V_{y}}{∂x} |_{y=0} = 0  (11)

a tengelyen.

A χ-re vonatkozó további határfeltételek vizsgálatára vegyük figyelembe, hogy Ψ állandó az akadály falai mentén. Így

ΔΨ = χ ⟹ a falaknál \{ \begin{matrix}\frac{∂^{2}Ψ}{∂y^{2}} = χ & ha a fal || x \\ \frac{∂^{2}Ψ}{∂x^{2}} = χ & ha a fal || y \end{matrix}  (12)

Súrlódó folyadék esetén indokolt lenne azt előírni, hogy a folyadék tapad a falhoz, azaz a tangenciális áramlás sebessége nulla. Ezt azonban nem tehetjük meg az áramlási függvénnyel, ugyanis ekkor az elliptikus peremértékprobléma túlhatározottá válna. A tapadást tehát nem írjuk elő közvetlenül Ψ-re, de közvetve majd felhasználjuk, amikor az előbbi egyenleteket diszkretizáljuk a χ-re vonatkozó határfeltételhez. A balról be és jobbra kifolyó örvénymentes áramlás miatt ezeken a határokon

V (x = ±∞, y) = [V_{0}, 0] ⟹ \frac{∂Ψ}{∂x} = 0 , \frac{∂Ψ}{∂y} = V_{0} és χ = 0 ha |x| → ∞  (13)

Végül a felső félteret zárt peremmel látjuk el, ha rendelkezünk arról, hogy az akadálytól y irányban messze milyen az áramlás. Kézenfekvő azt mondani, hogy elegendően nagy y távolságra az akadály hatása elenyésző, így a szabad áramlásra jellemző

\frac{∂Ψ}{∂x} = 0 , \frac{∂Ψ}{∂y} = V_{0} és χ = 0 ha |y| → ∞  (14)

asszimptotikus alakok választhatók.

IV. DISZKRETIZÁLT EGYENLETEK

A számítást x_{i} = (i − 1) ∗ h, i = 1, 2, ..., Nx és y_{j} = (j − 1) ∗ h, j = 1, 2, ..., Ny gridpontokon végezzük el. Célszerű a hosszúság és a sebesség egységeit úgy választani, hogy h = 1 és V_{0} = 1 legyen. Az áramlási függvényt hV_{0}, az örvényességet V_{0}/h a nyomást pedig ρV_{0}^{2} egységekben fogjuk mérni.

A. Áramlási függvény

A diszkretizált differenciálegyenlet

Ψ_{i+1,j} + Ψ_{i−1,j} + Ψ_{i,j+1} + Ψ_{i,j−1} − 4Ψ_{i,j} = χ_{i,j}  (15)

A határfeltételek által közvetlenül nem érintett belső pontokban

Ψ_{i,j} = \frac{1}{4} [Ψ_{i+1,j} + Ψ_{i−1,j} + Ψ_{i,j+1} + Ψ_{i,j−1} − χ_{i,j}]  (16)

• Baloldalt (x = 0)

\frac{∂Ψ}{∂x}|_{x=0} = 0 → Ψ_{2,j} − Ψ_{1,j} = 0  (17)

miatt

Ψ_{3,j} + Ψ_{1,j} + Ψ_{2,j+1} + Ψ_{2,j−1} − 4Ψ_{2,j} = χ_{2,j} → Ψ_{3,j} + Ψ_{2,j+1} + Ψ_{2,j−1} − 3Ψ_{2,j} = χ_{2,j}  (18)

azaz

Ψ_{2,j} = \frac{1}{3} [Ψ_{3,j} + Ψ_{2,j+1} + Ψ_{2,j−1} − χ_{2,j}] és Ψ_{1,j} = Ψ_{2,j}  (19)

• Jobboldalt (x = X)

\frac{∂Ψ}{∂x}|_{x=X} = 0 → Ψ_{X,j} − Ψ_{X−1,j} = 0  (20)
Ψ_{X−1,j} = \frac{1}{3} [Ψ_{X−2,j} + Ψ_{X−1,j+1} + Ψ_{X−1,j−1} − χ_{X−1,j}] és Ψ_{X,j} = Ψ_{X−1,j}  (21)

• Felül (y = Y)

\frac{∂Ψ}{∂y}|_{y=Y} = V_{0} → Ψ_{i,Y} − Ψ_{i,Y −1} = hV_{0} = 1  (22)
Ψ_{i,Y −1} = \frac{1}{3} [Ψ_{i+1,Y −1} + Ψ_{i−1,Y −1} + Ψ_{i,Y −2} + 1 − χ_{i,Y −1}] és Ψ_{i,Y} = Ψ_{i,Y −1} + 1  (23)

• Felső sarkokban

Ψ_{2,Y −1} = \frac{1}{2} [Ψ_{3,Y −1} + Ψ_{2,Y −2} + 1 − χ_{2,Y −1}]  (24)
Ψ_{1,Y −1} = Ψ_{2,Y −1} Ψ_{2,Y} = Ψ_{2,Y −1} + 1 Ψ_{1,Y} = Ψ_{2,Y}  (25)

és

Ψ_{X−1,Y −1} = \frac{1}{2} [Ψ_{X−2,Y −1} + Ψ_{X−1,Y −2} + 1 − χ_{X−1,Y −1}]  (26)
Ψ_{X−1,Y} = Ψ_{X−1,Y −1} + 1 Ψ_{X,Y −1} = Ψ_{X−1,Y −1} Ψ_{X,Y} = Ψ_{X−1,Y}  (27)

B. Örvényesség

A diszkretizációval kapjuk, hogy

[ \begin{matrix}(χ_{i+1,j} − 2χ_{i,j} + χ_{i−1,j}) \\ ↪ + (χ_{i,j+1} − 2χ_{i,j} + χ_{i,j−1}) \end{matrix} ] = \frac{R}{4} [ \begin{matrix}(Ψ_{i,j+1} − Ψ_{i,j−1}) & (χ_{i+1,j} − χ_{i−1,j}) \\ ↪ − (Ψ_{i+1,j} − Ψ_{i−1,j}) (χ_{i,j+1} − χ_{i,j−1}) \end{matrix} ] (28)

ahol R = V_{0}h/ν a rács Reynolds száma.

• A határfeltétel a hátsó falnál legyen most az, hogy

\frac{∂χ}{∂x}|_{x=X} = 0 → χ_{X,j} − χ_{X−1,j} = 0  (29)
[ \begin{matrix}(χ^{↓}_{i,j} −2χ_{i,j} + χ_{i−1,j}) \\ ↪ + (χ_{i,j+1} − 2χ_{i,j} + χ_{i,j−1}) \end{matrix} ] = \frac{R}{4} [ \begin{matrix}(Ψ_{i,j+1} − Ψ_{i,j−1}) & (χ^{↓}_{i,j} −χ_{i−1,j}) \\ ↪  − (Ψ_{i+1,j} − Ψ_{i−1,j}) (χ_{i,j+1} − χ_{i,j−1}) \end{matrix} ] i = X − 1 (30)

ahonnan

[3 + \frac{R}{4} (Ψ_{i,j+1} − Ψ_{i,j−1})] χ_{i,j} = χ_{i−1,j} + χ_{i,j+1} + χ_{i,j−1}  (31)
+\frac{R}{4} (Ψ_{i+1,j} − Ψ_{i−1,j}) (χ_{i,j+1} − χ_{i,j−1}) + \frac{R}{4} (Ψ_{i,j+1} − Ψ_{i,j−1}) χ_{i−1,j} (32)

• Az akadály felső peremén

\frac{∂^{2}Ψ}{∂y^{2}} = χ miatt Ψ_{i,j+1} − 2Ψ_{i,j} + Ψ_{i,j−1} = χ_{i,j}  (33)

lenne, de ebben Ψ_{i−1,j}, mint az akadály belsejében levő pont értelmetlen. Ha most a tapadást figyelembe vesszük, azaz használjuk, hogy

\frac{∂Ψ}{∂y} ≈ 0 ⟹ Ψ_{i,j+1} ≈ Ψ_{i,j−1}  (34)

akkor kapjuk, hogy

χ_{i,j} = 2 (Ψ_{i,j+1} − Ψ_{i,j})  (35)

• Hasonlóan kapjuk az akadály bal és jobb oldalán, hogy

χ_{i,j} = 2 (Ψ_{i−1,j} − Ψ_{i,j}) illetve χ_{i,j} = 2 (Ψ_{i+1,j} − Ψ_{i,j})  (36)