Návrh modelu proudění plynu Design of gas flow model

TECHNICKÁ UNIVERZITA V LIBERCI

Fakulta mechatroniky a mezioborových inženýrských studií

Studijní program: M 2612 – Elektrotechnika a informatika Studijní obor: 3901T025 – Přírodovědné inženýrství

Návrh modelu proudění plynu Design of gas flow model

Diplomová práce

Autor: Jakub Plešinger

Vedoucí práce: Ing. Jan Šembera, Ph.D.

V Liberci 9. 5. 2007

Prohlášení

Byl(a) jsem seznámen(a) s tím, že na mou diplomovou práci se plně vztahuje zákon č. 121/2000 o právu autorském, zejména § 60 (školní dílo).

Beru na vědomí, že TUL má právo na uzavření licenční smlouvy o užití mé diplomové práce a prohlašuji, že s o u h l a s í m s případným užitím mé diplomové práce (prodej, zapůjčení apod.).

Jsem si vědom(a) toho, že užít své diplomové práce či poskytnout licenci k jejímu využití mohu jen se souhlasem TUL, která má právo ode mne požadovat přiměřený příspěvek na úhradu nákladů, vynaložených univerzitou na vytvoření díla (až do jejich skutečné výše).

Diplomovou práci jsem vypracoval(a) samostatně s použitím uvedené literatury a na základě konzultací s vedoucím diplomové práce a konzultantem.

Datum 27.10.2007

Podpis

Poděkování

Na tomto místě bych rád poděkoval vedoucímu mé diplomové práce panu Ing.

Janu Šemberovi, Ph.D. za odborné vedení, trpělivost, cenné rady a připomínky, poskytnuté informace a zapůjčení odborné literatury.

Děkuji

Anotace

Tato práce se zabývá implementací a numerickým řešením úlohy vazkého izotermického proudění stlačitelné tekutiny pomocí metody konečných diferencí (FDM – Finite Difference Method) a metody konečných objemů (FVM – Finite Volume Method).

Pro metodu konečných diferencí je navrženo explicitní a implicitní schéma a pro metodu konečných objemů je to explicitní Lax-Friedrichsovo schéma.

Pro matematický popis proudění jsou použity Navier–Stokesova rovnice, která popisuje zákon zachování hybnosti vazké tekutiny, rovnice kontinuity, která vyjadřuje zákon zachování hmoty, a stavová rovnice plynu.

Všechna schémata se porovnávají na navržených testovacích úlohách, na dvojdimenzionální obdélníkové oblasti bez překážek a jednou vstupní hranou, jednou výstupní hranou a dvěma nepropustnými hranami. V práci se porovnává vliv velikosti časového kroku, prostorového kroku a dynamické viskozity na stabilitu jednotlivých schémat.

Klíčová slova

metoda konečných diferencí, metoda konečných objemů, izotermické proudění

Annotation

The topic of this work is numerical solution and implementation of isothermal viscous compressible flow using the Finite Difference Method and Finite Volume Method.

There are different schemes proposed for both methods – implicit and explicit scheme for the Finite Difference Method and the Lax-Friedrichs explicit scheme for the Finite Volume Method.

There is a mathematical description of the flow, which is based on the following three equations: Navier–Stokes equation which represents the momentum conservation law, continuity equation which represents the mass conservation law and the gas state equation.

All of those schemes are compared on testing problems. The flow is solved in a two-dimensional rectangular space without barriers. There are one inlet side and one outlet side in the area. Two other sides are impermeable. In this work, the effect of size of time and spatial discretization step and size of dynamic viscosity parameter on results of individual schemes is compared.

Keywords

finite volume method, finite difference method, isothermal flow

Obsah

Použité symboly...7

Úvod ...8

1. Fyzikální popis...9

1.1 Řídící rovnice...9

1.2 Počáteční a okrajové podmínky...9

2. Matematický model ...11

2.1 Operátory v rovnicích ...11

2.2. Rovnice ...11

2.3. Metoda konečných diferencí...12

2.3.1 Síť pro metodu konečných diferencí ...14

2.3.2 Časová diskretizace ...15

2.3.3 Explicitní schéma...16

2.3.4 Implicitní schéma...17

Newtonova metoda ...18

2.4 Metoda konečných objemů...21

2.4.1 Síť pro metodu konečných objemů...24

2.4.2 Lax-Friedrichsovo schéma nahrazení numerického toku H ...24

3. Implementace...30

3.1. Proměnné v programu...30

3.2. Popis funkcí v programu...32

3.3. Zobrazování výsledků...32

4. Testovací úlohy...34

4.1 Nulové proudění ...34

4.2 Ustálené proudění ve směru osy x ...35

4.3 Neustálené proudění ve směru osy x ...35

4.4 Ustálené proudění ve směru osy x se zavedenou chybou...35

4.5 Spočítané testovací úlohy ...36

5. Diskuse výsledků ...48

6. Závěr ...51

Literatura ...52

Příloha - Obsah CD...53

Použité symboly

Ω oblast v prostoru R²

Γ D hranice oblasti, na které je předepsaná Dirichletova okrajová podmínka Γ N hranice oblasti, na které je předepsaná Neumanova okrajová podmínka ρ hustota plynu, kg · m^-3

u v

= ⎜ ⎟⎛ ⎞

v ⎝ ⎠ vektor rychlosti plynu, m · s^-1 p tlak plynu, Pa

e tenzor rychlosti deformace

R plynová konstanta, R = 8, 31441 J · K^-1· mol^-1 T teplota, K

λ druhá viskozita, N · s · m^-2 µ dynamická viskozita, N · s · m^-2

∆x diskretizační parametr sítě ve směru osy x, m

∆ y diskretizační parametr sítě ve směru osy y, m

∆t časový diskretizační parametr, s

Úvod

Cílem diplomové práce je navrhnout model vazkého izotermického proudění stlačitelné tekutiny pomocí metody konečných diferencí (Finite Difference Method - FDM) a pomocí metody konečných objemů (Finite Volume Method - FVM) a na vhodně formulovaných úlohách jej otestovat a porovnat.

Na Ústavu nových technologií a aplikované informatiky FM TUL je řešen problém s časovou diskretizací při výpočtu proudění tekutin. Jednou z možností řešení časového vývoje proudění popsaného nelineárními rovnicemi se zabývá tato práce. Obsahuje celkem tři různé návrhy řešení soustavy. Plně implicitní schéma a plně explicitní schéma pro metodu konečných diferencí a explicitní Lax-Friedrichsovo schéma pro metodu konečných objemů.

Výpočet proudění je prováděn na obdélníkové oblasti orientované ve směru os souřadného systému. Oblast představuje jednoduchý kanál s jednou vstupní a jednou výstupní stěnou a dvěma nepropustnými stěnami.

Práce obsahuje fyzikální popis proudění pomocí Navier-Stokesových rovnic, rovnice kontinuity a stavové rovnice plynu, přechod k matematickému popisu pomocí metody konečných diferencí a metody konečných objemů, interpretaci okrajových podmínek a porovnání všech uvedených metod. Dále je zde velice stručně popsána implementace v Matlabu. V další části zprávy jsou pak navrženy základní testovací úlohy a uvedeny dílčí výsledky některých z nich. Zdrojový text programu a podprogram pro zobrazování výsledků v Matlabu a výsledky spočítaných testovacích úloh jsou na přiloženém CD-ROMu.

1. Fyzikální popis

Model proudění vzduchu je navržen pro dvourozměrnou oblast. Proudění předpokládáme stlačitelné, izotermické a vazké. Pro popis proudění plynu jsme použili rovnici kontinuity, která vyjadřuje zákon zachování hmoty, Navier–Stokesovu rovnici, která popisuje zákon zachování hybnosti vazké tekutiny, a stavovou rovnici plynu, která dává do souvislosti tlak a hustotu. Rovnice jsou v následujících tvarech:

1.1 Řídící rovnice

Rovnice kontinuity

( )

div 0

t

ρ ρ

∂ + =

∂ v . (1.1)

Navier–Stokesova rovnice

( )

^{1 p 1}

(

^div

) (

²

)

⁰

t λ µ

ρ ρ

∂ + ⋅ + − − =

∂

v v grad v grad grad v div e , (1.2)

kde u

v

= ⎜ ⎟⎛ ⎞

v ⎝ ⎠ je vektor rychlosti proudění vzduchu, ρ je hustota vzduchu, p je tlak vzduchu, λ je druhá viskozita, µ je dynamická viskozita, e je tenzor rychlosti deformace. Obě viskozity nejsou navzájem nezávislé a lze je dát do vztahu 3λ+2µ= 0 Stavová rovnice plynu – definuje vztah mezi hustotou, tlakem a teplotou

p 0

ρ⁻RT ^{= , (1.3)}

kde R je plynová konstanta a T je teplota, kterou budeme považovat za konstantu.

1.2 Počáteční a okrajové podmínky

Počáteční podmínkou je rychlostní pole a hustota vzduchu ve sledované oblasti v čase t = 0. Konkrétní hodnoty zadané v počátečních podmínkách se liší podle typu počítané úlohy. Obecně je můžeme zapsat ve tvaru

( )

0, = 0

( )

, ρ

( )

0, =ρ0

( )

, p

( )

0, = p0

( )

=ρ0

( )

RT

v x v x x x x x x . (1.4)

Pomocí okrajových podmínek definujeme potřebné veličiny na hranici oblasti.

Budeme předepisovat dva druhy okrajových podmínek.

- Dirichletova okrajová podmínka

D D

D D

D D D

na na

p p RT na

ρ ρ ρ

= Γ

= Γ

= = Γ

v v

(1.5)

- Neumanova okrajová podmínka

0

0

N N

N

N

v na na

p RT na

ρ

ρ

⋅ = Γ

⋅ = Γ

⋅ = ⋅ = Γ

v n grad n

grad n grad

(1.6)

kde n je vnější normálový vektor oblasti.

V modelu jsou použity následující výše uvedené okrajové podmínky. Dirichletovy okrajové podmínky (1.5) na částech hranice, kde plyn vstupuje do oblasti Ω a kde plyn z oblasti vystupuje. Protože přes výstupní hranici plyn z oblasti vystupuje, musí zde být předepsána rychlost vN s opačným znaménkem než na vstupní hranici. Neumanovy okrajové podmínky (1.6) jsou předepsány na nepropustných částech hranice a tedy vN = 0.

Obr. 1.1. Oblast, ve které je úloha proudění řešena

2. Matematický model

Pro snadnější programování modelu proudění plynu rovnice rozepíšeme po složkách a všechny operátory v rovnicích rozepíšeme do jednotlivých derivací.

2.1 Operátory v rovnicích

Za zmínku stojí čtyři operátory. V Navier-Stokesově rovnici (1.2) se nacházejí následující operátory, z nichž první se nazývá konvektivní derivace

( )

u u

u v

u x y

u v

v v v

x y

u v

x y

∂ ∂

⎛ + ⎞

⎜ ∂ ∂ ⎟

⎛ ∂ ∂ ⎞⎛ ⎞ ⎜ ⎟

⋅ =⎜⎝ ∂ + ∂ ⎟⎜ ⎟⎠⎝ ⎠=⎜⎜⎝ ∂∂ + ∂∂ ⎟⎟⎠

v grad v , (2.1)

( )

2 2

2

2 2

2

div

u v

x x y

u v

x y y

λ λ

⎛∂ ∂ ⎞

⎜∂ +∂ ∂ ⎟

⎜ ⎟

= ⎜⎜ ∂ +∂ ⎟⎟

∂ ∂ ∂

⎝ ⎠

grad v , (2.2)

( )

2 2 2

2 2

2 2 2

2 2

2 2

2

u v u

x x y y

v u v

y x y x

µ µ

⎛ ∂ + ∂ +∂ ⎞

⎜ ∂ ∂ ∂ ∂ ⎟

⎜ ⎟

= ⎜⎜ ∂ + ∂ +∂ ⎟⎟

∂ ∂ ∂ ∂

⎝ ⎠

div e (2.3)

a v rovnici kontinuity (1.1) je to operátor divergence

div( ) div( ) ( ) u v

u v

x y x y

ρ ρ

ρv = ⋅ρ v + ⋅v grad ρ =ρ^⎛⎜⎝^∂∂ +^∂∂ ^⎞⎟⎠+ ^∂∂ + ^∂∂ . (2.4)

2.2. Rovnice

Když rovnice popisující proudění (1.1) a (1.2) rozepíšeme po složkách a dosadíme stavovou rovnici (1.3) a operátory rozepsané v předcházející části, dostaneme výsledné rovnice ve tvaru:

2 2 2 2 2

2 2 2

2 2 2 2 2

2 2 2

2

0 2

u u u v u v u

u u v

x y RT x x x y x x y y

t

v u v v v u v v u v

t x y y x y y y x y x

ρ

λ µ

ρ ρ

⎛ ⎞ ⎛ ⎞

∂ ∂ ∂ ∂ ∂ ∂ ∂

⎛ ⎞ ⎛∂ ⎞

⎛∂ ⎞ ⎜ + ⎟ ⎜ ⎟ ⎜ + ⎟ ⎜ + + ⎟

⎜∂ ⎟+⎜ ∂ ∂ ⎟+ ⎜∂ ⎟− ⎜∂ ∂ ∂ ⎟− ⎜ ∂ ∂ ∂ ∂ ⎟=

⎜∂ ⎟ ⎜ ∂ ∂ ⎟ ⎜∂ ⎟ ⎜ ∂ ∂ ⎟ ⎜ ∂ ∂ ∂ ⎟

⎜ ⎟ + + + +

⎜∂ ⎟ ⎜ ∂ ∂ ⎟ ⎜∂ ⎟ ⎜ ⎟ ⎜ ⎟

⎝ ⎠ ⎝ ⎠ ⎝ ⎠ ⎝∂ ∂ ∂ ⎠ ⎝ ∂ ∂ ∂ ∂ ⎠

(2.5)

u v 0

u v

t x y x y

ρ ρ^{⎛ ⎞} ρ ρ

∂∂ + ⎜⎝∂∂ +∂∂ ⎟⎠+ ∂∂ + ∂∂ = (2.6)

2.3. Metoda konečných diferencí

Základní myšlenku metody konečných diferencí, zvané též metoda sítí, ukážeme na následujícím příkladu. Uvažujme jednorozměrnou oblast ve směru osy x a prostorovou diskretizaci s n+1 diskretizačními body v intervalu <a,b>. V tomto intervalu hledáme řešení dané okrajové úlohy. Pro jednotlivé body xk platí předpis

,x_k = +a kh k=0,1,...,n, (2.7)

kde

h b a n

= − (2.8)

je diskretizační krok a n je přirozené číslo. V bodech této sítě se splní daná diferenciální rovnice přibližně, a to tak, že derivace, které se v nich vyskytují, se nahradí diferenčními podíly (tj. lineárními kombinacemi funkčních hodnot v okolních bodech). Vztahy pro nahrazení derivace diferenciálním vzorcem vycházejí z Taylorova rozvoje.

Pokud funkci y(x+ h) nahradíme Taylorovým rozvojem do 1. řádu, dostaneme

( ) ( )

^'

( ) ( )

²

y x h+ = y x + ⋅h y x +O h (2.9)

a pro derivaci platí vztah

( ) ( ) ( ) ( )

' y x h y x

y x O h

h

+ −

= + , (2.10)

kde O(h) je diskretizační chyba omezená násobkem první mocniny diskretizačního kroku.

Zlomek na pravé straně (2.10) nazýváme diskretizační podíl.

Uvažujme výše uvedený příklad jednorozměrné oblasti s n+1 diskretizačními body xk a diskretizačním krokem h. Funkční hodnoty v těchto bodech značíme y(xk) a aproximaci funkčních hodnot označíme yk. První derivaci funkce y v bodě xk pak můžeme definovat mnoha způsoby. Vycházíme-li ze vztahu

( ) (

¹

) ( ) ( )

'

k

k k

k

x x

y x y x

y x dy O h

dx h

+

=

⎛ ⎞ −

=⎜⎝ ⎟⎠ = + , (2.11)

pak dostaneme diferenční podíl, který nazýváme dopředná diference

( )

¹

' _k y^k y^k

y x h

+ −

≈ , (2.12)

vycházíme-li z jiného vztahu pro první derivaci

( ) ( ) ( )

¹

( )

'

k

k k

k

x x

y x y x

y x dy O h

dx h

−

=

⎛ ⎞ −

=⎜⎝ ⎟⎠ = + , (2.13)

dostáváme vztah pro zpětnou diferenci

( )

¹

' _k y^k y^k

y x h

− −

≈ . (2.14)

Oba vztahy jsou aproximacemi prvního řádu a jsou nesymetrické vzhledem k bodu xk. Pokud chceme použít symetrického vzorce, užijeme vztah pro centrální diference, která je již aproximací druhého řádu

( ) (

¹

) ( )

¹

( )

^{2 1 1}

' _k 2 2

k k k k

k

x x

y x y x y y

y x dy O h

dx h h

+ − + −

=

− −

⎛ ⎞

=⎜ ⎟ = + ≈

⎝ ⎠ . (2.15)

Pokud chceme aproximovat derivace vyšších řádů, použijeme opět Taylorova rozvoje a odvodíme diferenční vzorec pro příslušnou derivaci. Příkladem může být vztah pro druhou derivaci. Nejprve vyjádříme Taylorovy rozvoje pro y(x + h) a y(x - h) do 3.

řádu:

( ) ( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( ) ( ) ( )

2 3

4

2 3

4

' '' '''

2 6

' '' ''' .

2 6

h h

y x h y x h y x y x y x O h

h h

y x h y x h y x y x y x O h

+ = + ⋅ + ⋅ + ⋅ +

− = − ⋅ + ⋅ − ⋅ +

(2.16)

Tyto vztahy sečteme

( ) ( )

²

( )

^{2 ''}

( ) ( )

⁴

y x h+ +y x h− = y x +h y x⋅ +O h (2.17)

a po drobných úpravách dostaneme vztah pro druhou derivaci

( )

²2

(

¹

) ( ) ( )

2 ¹

( )

^{2 1} 2 ¹

2 2

''

k

k k k k k k

k

x x

y x y x y x y y y

y x d y O h

dx h h

+ − + −

=

− +

⎛ ⎞ − +

=⎜ ⎟ = + ≈

⎝ ⎠ . (2.18)

Výše uvedené diferenční vzorce jsou nejjednodušší příklady nahrazení derivací.

Obecně do vzorců nemusíme používat jen hodnoty v nejbližších bodech, ale pro speciální případy se používají i hodnoty ve vzdálenějších bodech. Pro jinak zvolené diferenční vzorce má i schéma jiné vlastnosti.

Pro parciální diferenciální rovnice používáme vícerozměrnou síť, v našem případě dvourozměrnou. Opět použijeme ekvidistantní rozdělení prostoru. Diskretizační

kroky mohou být pro jednotlivé směry různé a označíme je ∆x, ∆y. Body sítě pak mají souřadnice (xk,yl)

0

0

, 0,1,..., , 0,1,..., .

k

s

x x k x k n

y y l y l m

= + ∆ =

= + ∆ = (2.19)

Sousední body sítě jsou body, které se od sebe liší vždy v jednom indexu právě o 1.

Regulární bod (vnitřní bod) je takový bod, jehož všechny sousední body leží uvnitř oblasti. Parciální derivace nahrazujeme v regulárních bodech pomocí obdobných diferenčních vzorců jako totální derivace v jednorozměrném případě. Jako příklad uvádím již odvozené diferenční vztahy pro centrální diference

1, 1,

, 1 , 1

2 2 ,

k l

k l

k l k l

x x y y

k l k l

x x y y

w w

w

x x

w w

w

y y

+ −

=

=

+ −

=

=

∂ −

⎛ ⎞ ≈

⎜∂ ⎟ ∆

⎝ ⎠

⎛∂ ⎞ −

⎜∂ ⎟ ≈ ∆

⎝ ⎠

(2.20)

a diferenční vztah pro smíšenou derivaci 2.řádu

1, 1 1, 1 1, 1 1, 1

4

k l

k l k l k l k l

x x y y

w w w w

w

x y x y

+ + + − − + − −

=

=

− − +

⎛ ∂ ⎞

⎜∂ ∂ ⎟ ≈ ∆ ∆

⎝ ⎠ , (2.21)

která se ve schématu nachází.

Okrajová podmínka na nepropustné části hranice je popsána nulovou normálovou složkou rychlosti a tečnou složkou stejnou jako v nejbližším regulárním bodě. Na vtokové a výtokové straně se do hraničních bodů zapisují vstupní hodnoty normálové rychlosti a hustoty, tečná rychlost je na vtokové straně nulová. Podrobnější vysvětlení použití metody konečných diferencí pro parciální diferenciální rovnice viz [3].

2.3.1 Síť pro metodu konečných diferencí

Pro obě následující schémata navržená pro metodu konečných diferencí je použita stejná síť, která pokrývá celou počítanou oblast. Parametry sítě jsou ∆x a ∆y. Okrajové podmínky jsou rozlišeny podle stěn, na kterých má být předepsána.

Obr. 2.1. Příklad sítě použité pro metodu konečných diferencí

2.3.2 Časová diskretizace

V soustavě rovnic se vyskytují pouze první časové derivace proměnných. V Navier-Stokesově rovnici je to derivace rychlosti a v rovnici kontinuity je to derivace hustoty. Pro jednodušší popis časové diskretizace můžeme soustavu přepsat ve tvaru

2 2 2

2 2

, , , , , ,t

t x y x y x y

⎛ ⎞

∂ ∂ ∂ ∂ ∂ ∂

= ⎜ ⎟

∂ ⎝ ∂ ∂ ∂ ∂ ∂ ∂ ⎠

w w w w w w

f w , (2.22)

kde

u v ρ

⎛ ⎞⎜ ⎟

= ⎜ ⎟

⎜ ⎟⎝ ⎠

w , (2.23)

kde u,v jsou složky vektoru rychlosti a ρ hustota.

Pro diferenciální funkci , , , ²₂ , ² ₂ , ² ,t

x y x y x y

⎛ ∂ ∂ ∂ ∂ ∂ ⎞

⎜ ∂ ∂ ∂ ∂ ∂ ∂ ⎟

⎝ ⎠

w w w w w

f w sestavíme její diskrétní formu

nahrazením derivací w diferenčními podíly a vztahy vyplývajícími z dosazení okrajových podmínek

( ) ( )

( )

1 D

D D

N D

f

f

⎛ ⎞

⎜ ⎟

= ⎜ ⎟

⎜ ⎟

⎝ ⎠

w f w

w

, (2.24)

kde

1 D

N

w

w

⎛ ⎞

⎜ ⎟

= ⎜ ⎟

⎜ ⎟

⎝ ⎠

w (2.25)

je vektor aproximací funkčních hodnot v jednotlivých bodech sítě. N je počet bodů v síti.

Diskrétní forma také obsahuje okrajové podmínky.

2.3.3 Explicitní schéma

Nahradíme-li časovou derivaci v rovnicích (1.1), (1.2) pomocí následujícího diferenčního vzorce

( ) ( )

t

t t t

t t

+ ∆ −

⎛∂ ⎞ ≈

⎜ ∂ ⎟ ∆

⎝ ⎠

w w

w (2.26)

a pokud všechny proměnné, které se vyskytují v diskretizované rovnici mimo časovou derivaci, dosadíme ze starého kroku tzn. z času t, dostaneme explicitní schéma

( )

1

n n n

D ⁺ = D D ⋅ ∆ −t D

w f w w . (2.27)

Indexy n+1 a n označují nový a starý časový krok.

Použitím tohoto schématu dostaneme systém N rovnic s N neznámými, ale v každé rovnici je zastoupena pouze jedna neznámá. Zbavili jsme se tím nelineární soustavy rovnic a dostali jsme přímo vztahy pro jednotlivé neznámé. Obecnou nevýhodou explicitních schémat je jejich možná nestabilita. Velice záleží na zvolené diskretizaci – na parametrech

∆x, ∆y a ∆t. Explicitní schéma je pro soustavu lineárních rovnic stabilní při splnění určitých podmínek, které určují vztah mezi časovým a prostorovým krokem – takové schéma se nazývá podmíněně stabilní viz [2].

Rovnice pro explicitní metodu

Navier-Stokesova rovnice pro rychlost ve směru osy x:

1, 1, , 1 , 1 1, 1,

1

, , , ,

,

1, , 1, 1, 1 1, 1 1, 1 1, 1

2 ,

2 2 2

2

4

n n n n n n

k l k l k l k l k l k l

n n n n

k l k l k l k l n

k l

n n n n n n n

k l k l k l k l k l k l k l

n k l

u u u u RT

u u t u v

x y x

u u u v v v v

t x x

ρ ρ

ρ λ

ρ

+ − + − + −

+

+ − + + − + + − − −

⎡ ⎛ − ⎞ ⎛ − ⎞ ⎛ − ⎞⎤

= − ∆ ⎢⎢⎣ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⎜⎜⎝ ∆ ⎟⎟⎠⎥⎥⎦+

− + − − +

+∆ +

∆ ∆ ∆

1, , 1, , 1 , , 1 1, 1 , 1 1, ,

2 2

2 2

2 4

n n n n n n n n n n

k l k l k l k l k l k l k l k l k l k l

y

u u u u u u v v v v

tµ ⁺ x ^{− +} y ^{− + + +}x y ⁺

⎛ ⎞

⎜ ⎟+

⎜ ⎟

⎝ ⎠

⎛ − + − + − − + ⎞

+∆ ⎜⎜⎝ ∆ + ∆ + ∆ ∆ ⎟⎟⎠

(2.28)

Navier-Stokesova rovnice pro rychlost ve směru osy y:

1, 1, , 1 , 1 , 1 , 1

1

, , , ,

,

1, 1 1, 1 1, 1 1, 1 , 1 , , 1

,

2 2 2

2 4

n n n n n n

k l k l k l k l k l k l

n n n n

k l k l k l k l n

k l

n n n n n n n

k l k l k l k l k l k l k l

n k l

v v v v RT

v v t u v

x y y

u u u u v v v

t x y y

ρ ρ

ρ λ

ρ

+ − + − + −

+

+ + − + + − − − + −

⎡ ⎛ − ⎞ ⎛ − ⎞ ⎛ − ⎞⎤

= − ∆ ⎢⎢⎣ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⎜⎜⎝ ∆ ⎟⎟⎠⎥⎥⎦+

− − + − +

+∆ +

∆ ∆ ∆ ²

, 1 , , 1 1, , 1, 1, 1 1, 1 1, 1 1, 1

2 2

2 2

2 4

n n n n n n n n n n

k l k l k l k l k l k l k l k l k l k l

v v v v v v u u u u

tµ ⁺ y ^{− +} x ^{− + + − +} x y^{+ − − −}

⎛ ⎞

⎜ ⎟+

⎜ ⎟

⎝ ⎠

⎛ − + − + − − + ⎞

+∆ ⎜⎜⎝ ∆ + ∆ + ∆ ∆ ⎟⎟⎠

(2.29)

Rovnice kontinuity

1, 1, , 1 , 1

1

, , ,

1, 1, , 1 , 1

, ,

2 2

2 2

n n n n

k l k l k l k l

n n n

k l k l k l

n n n n

k l k l k l k l

n n

k l k l

u u v v

t

RT x y

t t

u v

RT x RT y

ρ ρ ρ

ρ ρ ρ ρ

+ − + −

+

+ − + −

⎛ − − ⎞

= − ∆ ⋅ ⎜⎜⎝ ∆ + ∆ ⎟⎟⎠+

⎛ − ⎞ ⎛ − ⎞

∆ ∆

+ ⋅ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⋅ ⎜⎜⎝ ∆ ⎟⎟⎠

(2.30)

2.3.4 Implicitní schéma

Nahradíme-li časovou derivaci v rovnicích (1.1), (1.2) pomocí následujícího diferenčního vzorce

( ) ( )

t t

t t t

t _+∆ t

+ ∆ −

⎛∂ ⎞ ≈

⎜ ∂ ⎟ ∆

⎝ ⎠

w w

w (2.31)

a pokud všechny proměnné, které se vyskytují v diskretizované rovnici mimo časovou derivaci, dosadíme z nového časového kroku tzn. z času t+∆t, dostaneme implicitní schéma

( )

1 1

n n n

D ⁺ − D D ⁺ ⋅ ∆ =t D

w f w w . (2.32)

Indexy n+1 a n označují nový a starý časový krok.

Tím dostaneme mnohem složitější schéma, jehož výhodou je jeho stabilita.

Schéma je stabilnější než explicitní, může být dokonce nepodmíněně stabilní, což platí při řešení soustavy lineárních rovnic. Při řešení soustavy nelineárních rovnic tomu ovšem tak být nemusí. V případě úlohy proudění popsaného rovnicemi (1.1),(1.2),(1.3), která je nelineární, vede tato metoda na řešení N nelineárních rovnic pro N neznámých.

Rovnice pro implicitní metodu

Navier-Stokesova rovnice pro rychlost ve směru osy x:

1 1 1 1 1 1

1, 1, , 1 , 1 1, 1,

1 1 1

, , , 1

,

1 1 1 1 1

1, , 1, 1, 1 1, 1

1 2

,

2 2 2

2

n n n n n n

k l k l k l k l k l k l

n n n

k l k l k l n

k l

n n n n n

k l k l k l k l k l

n k l

u u u u RT

u t u v

x y x

u u u v v

t x

ρ ρ

ρ λ

ρ

+ + + + + +

+ − + − + −

+ + +

+

+ + + + +

+ − + + − +

+

⎡ ⎛ − ⎞ ⎛ − ⎞ ⎛ − ⎞⎤

+ ∆ ⎢⎢⎣ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⎜⎜⎝ ∆ ⎟⎟⎠⎥⎥⎦−

− + −

−∆ +

∆

1 1

1, 1 1, 1

1 1 1 1 1 1 1 1 1 1

1, , 1, , 1 , , 1 1, 1 1, 1 1, 1 1, 1

2 2 ,

4

2 2

2 4

n n

k l k l

n n n n n n n n n n

k l k l k l k l k l k l k l k l k l k l n

k l

v v

x y

u u u u u u v v v v

t u

x y x y

µ

+ +

+ − − −

+ + + + + + + + + +

+ − + − + + − + + − − −

⎛ − + ⎞

⎜ ⎟−

⎜ ∆ ∆ ⎟

⎝ ⎠

⎛ − + − + − − + ⎞

−∆ ⎜⎜⎝ ∆ + ∆ + ∆ ∆ ⎟⎟⎠=

(2.33)

Navier-Stokesova rovnice pro rychlost ve směru osy y:

1 1 1 1 1 1

1, 1, , 1 , 1 , 1 , 1

1 1 1

, , , 1

,

1 1 1 1

1, 1 1, 1 1, 1 1, 1

1 ,

2 2 2

4

n n n n n n

k l k l k l k l k l k l

n n n

k l k l k l n

k l

n n n n

k l k l k l k l k

n k l

v v v v RT

v t u v

x y y

u u u u v

t x y

ρ ρ

ρ λ

ρ

+ + + + + +

+ − + − + −

+ + +

+

+ + + +

+ + − + + − − −

+

⎡ ⎛ − ⎞ ⎛ − ⎞ ⎛ − ⎞⎤

+ ∆ ⎢⎢⎣ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⎜⎜⎝ ∆ ⎟⎟⎠+ ⎜⎜⎝ ∆ ⎟⎟⎠⎥⎥⎦−

− − +

−∆ +

∆ ∆

1 1 1

, 1 , , 1

2

1 1 1 1 1 1 1 1 1 1

, 1 , , 1 1, , 1, 1, 1 1, 1 1, 1 1, 1

2 2 ,

2

2 2

2 4

n n n

l k l k l

n n n n n n n n n n

k l k l k l k l k l k l k l k l k l k l n

k l

v v

y

v v v v v v u u u u

t v

y x x y

µ

+ + +

+ −

+ + + + + + + + + +

+ − + − + + − + + − − −

⎛ − + ⎞

⎜ ⎟−

⎜ ∆ ⎟

⎝ ⎠

⎛ − + − + − − + ⎞

−∆ ⎜⎜⎝ ∆ + ∆ + ∆ ∆ ⎟⎟⎠=

(2.34)

Rovnice kontinuity

1 1 1 1

1, 1, , 1 , 1

1 1

, ,

1 1 1 1

1, 1, , 1 , 1

1 1 1

, , ,

2 2

2 2

n n n n

k l k l k l k l

n n

k l k l

n n n n

k l k l k l k l

n n n

k l k l k l

u u v v

t

RT x y

t t

u v

RT x RT y

ρ ρ

ρ ρ ρ ρ

ρ

+ + + +

+ − + −

+ +

+ + + +

+ − + −

+ + +

⎛ − − ⎞

+ ∆ ⋅ ⎜⎜⎝ ∆ + ∆ ⎟⎟⎠−

⎛ − ⎞ ⎛ − ⎞

∆ ∆

− ⋅ ⎜⎜⎝ ∆ ⎟⎟⎠− ⋅ ⎜⎜⎝ ∆ ⎟⎟⎠=

(2.35)

K řešení této nelineární úlohy lze snadno využít Newtonovu metodu pro řešení nelineárních rovnic.

Newtonova metoda

Nelineární soustavu rovnic (2.33), (2.34), (2.35) přepíšeme do tvaru, abychom na pravé straně dostali nulu. Schématicky lze zapsat

(

u_{k l}ⁿ,⁺¹,v_{k l}ⁿ,⁺¹,ρ_{k l}ⁿ,⁺¹,u_{k l}ⁿ,,v_{k l}ⁿ,,ρ_{k l}ⁿ,

)

=0

F , (2.36)

kde u_{k l}ⁿ_,⁺¹,v_{k l}ⁿ_,⁺¹,ρ_{k l}ⁿ_,⁺¹ jsou hledané hodnoty rychlosti a hustoty ve všech bodech sítě v novém časovém kroku a u_{k l}ⁿ_,,v_{k l}ⁿ_,,ρ_{k l}ⁿ_, jsou známe hodnoty z posledního časového kroku. Pro jednoduchost zápisu přepíšeme rovnici (2.36) do tvaru

( )

⁼⁰

F x (2.37)

Za předpokladu, že existují všechny potřebné derivace, můžeme napsat Taylorův rozvoj funkce F v bodě x^k

( )

⁼

( ) ( )(

^{k + ' k − k}

)

⁺¹₂ ^''

^{( )}

^ξ

(

^{− k}

)

²

F x F x F x x x F x x , (2.38)

kde F x je Jacobiho matice funkce F(x) v bodě x^'

( )

^k k

( )

( ) ( )

( ) ( )

1 1

1

1

'

k k

N k

k k

N N

N

F F

x x

F F

x x

⎛ ∂ ∂ ⎞

⎜ ⎟

∂ ∂

⎜ ⎟

⎜ ⎟

= ⎜ ⎟

⎜∂ ∂ ⎟

⎜ ⎟

⎜ ∂ ∂ ⎟

⎝ ⎠

x x

F x

x x

,

kde N je rozměr úlohy.

V Jacobiho matici se vyskytují derivace jednotlivých rovnic podle všech proměnných. Protože derivace v původních rovnicích jsou nahrazeny diferencemi, které používají pouze hodnoty v sousedních bodech sítě, zjistíme, že většina derivací v matici bude nulová. V každé rovnici je teoreticky maximálně dvacet sedm proměnných. Jedná se o hodnoty rychlosti ve směrech x a y a hustoty v daném uzlu a ve všech jeho nejbližších sousedních bodech. Prakticky ovšem žádná rovnice tolik proměnných nemá. Matice je řídká a volbou uspořádání matice lze dosáhnout různého rozložení nenulových prvků v matici, což by mohlo být důležité pro jednodušší řešení úlohy. Následující tvar matice vznikl při řazení rovnic v pořadí Navier-Stokesova rovnice ve směru osy x (2.33), ve směru osy y (2.34) a rovnice kontinuity (2.35). Řazení proměnných je pak ve tvaru

1,1, 1,1, 1,1, , _m,1, _m,1, _m,1, , _{m n}, , _{m n}, , _{m n},

u v ρ … u v ρ … u v ρ .

Obr. 2.2. Příklad struktury nenulových prvků Jacobiho matice pro síť se 36 uzly (6x6 uzlů).

Soustavu rovnic (2.37) pak můžeme nahradit soustavou lineárních rovnic

( ) ( )(

^{k + ' k − k}

)

⁼⁰

F x F x x x . (2.39)

Řešení této soustavy pak označíme x^k+1 a pokud rozdíl

(

^{x x− k}

)

respektive

(

^xk+1−xk

)

^po

dosazení předpokládaného řešení nazveme h^k, dostaneme soustavu

( ) ( )

' ^{k k} = − ^k

F x h F x (2.40)

Tuto soustavu vyřešíme například pomocí inverse Jacobiho matice a dostaneme vektor hka pomocí něho pak novou hodnotu x^k+1

Tento postup řešení se opakuje až do té doby, dokud není splněna podmínka

( )

^{k ≤ε}

F x , (2.42)

kde ε je přesnost, se kterou chceme soustavu vyřešit. Ne vždy lze dosáhnout vyřešení soustavy s požadovanou přesností a proto se dále také předepisuje maximální počet iterací, aby nedošlo k nekonečnému cyklu při řešení úlohy.

2.4 Metoda konečných objemů

Pro další možné řešení úlohy proudění plynu jsme se rozhodli použít metodu konečných objemů viz [4]. Výhodou metody konečných objemů je, že pro výpočet nepotřebuje pravidelnou síť jako metoda konečných diferencí. Oblast, ve které provádíme výpočet, musíme rozdělit diskretizační sítí na konečný počet podoblastí. Abychom mohli použít metodu konečných objemů, musí zvolená diskretizace splňovat určité podmínky [4]. Ke každému objemu je definován bod, který je použit jako reprezentativní bod tohoto objemu. Metoda počítá průměrné funkční hodnoty v jednotlivých objemech sítě a tyto průměrné hodnoty jsou pak vyjádřeny jako funkční hodnoty v reprezentativních bodech.

V metodě konečných objemů převádíme diferenciální rovnici do integrálně- diferenciální formy v jednotlivých konečných objemech a objemové integrály, které obsahují divergenční výrazy, se převádějí na integrály přes povrch pomocí divergenční věty. Tyto integrály vyjadřují toky přes hranice jednotlivých konečných objemů.

Na odvození schématu pro výpočet proudění metodou konečných objemů nejprve původní soustavu rovnic (1.1), (1.2), (1.3) přepíšeme do tvaru, ve kterém jsou jednotlivé časové a prostorové derivace oddělené od sebe. Upravené rovnice viz [2] pak mají tvar

( ) ( ) ( ) ( )

1 2 1 2

1 2 1 2

t x x x x 0

∂ ∂ ∂ ∂

∂ + + − − =

∂ ∂ ∂ ∂ ∂

f w f w g w g w

w , (2.43)

kde

( ) ( )

2 3

1 2

2 2 3

2

2 1 1 2

1 1 2

3 2

2 3 3

1

1 1

, , ,

w w

w u v

w w

w u w kw u k uv

w w

w v uv v k

w w w

w w kw

ρ ρ ρ

ρ ρ ρ ρ

ρ ρ ρ ρ

⎛ ⎞

⎛ ⎞

⎜ ⎟

⎜ ⎟

⎜ ⎟

⎜ ⎟

⎛ ⎞ ⎛ ⎞ ⎜ ⎟ ⎛ ⎞ ⎜ ⋅ ⎟ ⎛ ⎞

⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟

=⎜⎜⎝ ⎟ ⎜⎟ ⎜⎠ ⎝= ⎟⎟⎠ =⎜⎜⎜⎜⎝ +⋅ ⎟⎟⎟⎟⎠=⎜⎜⎝ + ⎟⎟⎠ =⎜⎜⎜⎜⎜⎝ + ⎟⎟⎟⎟⎟⎠=⎜⎜⎝ + ⎟⎟⎠

w f w f w

( ) ( )

1 2

0 0

2 2 ,

3

2 2

3

u v u v

x y y x

u v u v

y x x y

µ µ

µ µ

⎛ ⎞ ⎛ ⎞

⎜ ⎟ ⎜ ⎟

⎜ ⎟ ⎜ ⎟

⎜ ⎛ ∂ ∂ ⎞⎟ ⎜ ⎛∂ ∂ ⎞ ⎟

⎜ ⎟ ⎜ ⎟

=⎜ ⎜⎝ ∂ −∂ ⎟⎠⎟ =⎜ ⎜⎝∂ +∂ ⎟⎠ ⎟

⎜ ⎟ ⎜ ⎟

⎛∂ ∂ ⎞ ⎛ ∂ ∂ ⎞

⎜ ⎜ + ⎟ ⎟ ⎜ ⎜− + ⎟⎟

⎜ ⎝∂ ∂ ⎠ ⎟ ⎜ ⎝ ∂ ∂ ⎠⎟

⎝ ⎠ ⎝ ⎠

g w g w ,

kde u,v jsou složky vektoru rychlosti, ρje hustota a k =RT.

Dalším krokem pro získání schématu pro metodu konečných objemů je integrace rovnice (2.43) přes celou oblast a v časovém intervalu

(

t t_k, _k+1

)

( ) ( ) ( )

1 1 1

1 1

, 0

k k k

k i k i k i

t t N t N

s s

s s s s

t D t D t D

x t dxdt dxdt dxdt

t x x

+ + +

= =

∂ ∂ ∂

+ − =

∂

∑

∂

∑

∂

∫ ∫

^w

∫ ∫

^{f w}

∫ ∫

^{g w . (2.44)}

Pokud předpokládáme spojitost v celé oblasti a v časovém intervalu

(

t t_k, _k+1

)

, můžeme zaměnit pořadí integrací u prvního členu. Pokud na druhý a třetí člen rovnice použijeme Greenovu větu, tak po dalších úpravách pak dostaneme

( ) ( ) ( )

1 1 1

1 1

, 0

k k k

i k i k i

k

t t N t N

s s s s

s s

D t t t D t D

x t dx n dSdt n dSdt

+ + +

= =

∂ ∂

=

+

∑

−

∑

=

∫

^w

∫ ∫

^{f w}

∫ ∫

^{g w , (2.45)}

kde D je objem, přes který integrujeme, _i ∂ je hranice objemu a N je rozměr úlohy, který D_i je v našem případě N = 2.

Rovnici (2.45) lze dále rozepsat do tvaru

( ) ( )

( ) ( ) ( )

( ) ( )

1

,

1

,

1

( ) 1

( ) 1

, ,

0,

k

i k i j

k

k i j

t N

k k s ij s

j S i s

D t

t N

s ij s

j S i s

t

x t x t dx dS dt

dS dt

+

+

+

∈ Γ =

∈ Γ =

⎛ ⎞

⎜ ⎟

− + −

⎜ ⎟

⎝ ⎠

⎛ ⎞

⎜ ⎟

− =

⎜ ⎟

⎝ ⎠

∫ ∫ ∑ ∑ ∫

∫ ∑ ∑ ∫

w w f w n

g w n

(2.46)

kde ^{S i}

( )

je množina sousedních objemů, Γ je hrana mezi sousedními objemy i a j, a _ij n _ij je normálový vektor hrany mezi sousedními objemy.

Pokud aproximujeme integrál

(

^,

)

^/

i

k i

D

x t dx D

∫

^{w pomocí w ki}

( )

1 ,

i

k

i k

i D

x t dx

≈ D

∫

w w , (2.47)

toky fs pomocí numerických toků H

( ) ( ) ( )

1

, ,

N l l

s ij s i j ij

s=

∑

^{f w n} ≈H w w n , (2.48)

kde l může být k+1 nebo k,

a toky gs pomocí numerických toků G

( ) ( ) ( )

1

, ,

N l l

s ij s i j ij

s=

∑

^{g w n} ≈G w w n , (2.49)

dostaneme již schéma pro metodu konečných objemů

( ) ^{( )} ( )

( ) ^{( )} ( )

1 1 1

1 1

, , 1 , ,

, , 1 , , ,

k k k k k k

i i i j ij i j ij ij

i j S

k k k k

i j ij i j ij ij

i j S

t D t

D

ϑ ϑ

ϑ ϑ

+ + +

∈

+ +

∈

∆ ⎡ ⎤

= − ⎣ + − ⎦ Γ +

∆ ⎡ ⎤

+ ⎣ + − ⎦ Γ

∑

∑

w w H w w n H w w n

G w w n G w w n

(2.50)

kde ϑje parametr, který zle schéma měnit od plně explicitního k plně implicitnímu.

Volbou parametru ϑ =0získáme explicitní schéma

( ) ( )

1 , , , ,

k k k k k k

i i i j ij ij i j ij ij

j S j S

i i

t t

D D

+

∈ ∈

∆ ∆

= −

∑

Γ +

∑

Γ

w w H w w n G w w n . (2.51)

2.4.1 Síť pro metodu konečných objemů

Pro metodu konečných objemů jsme zvolili pravoúhlou síť podle obrázku 2.3. Je patrné, že všechny objemy neležící na hranici oblasti mají rozměry ∆x a ∆y. Objemy ležící na jedné hranici mají rozměry ∆x/2 a ∆y pro vstupní a výstupní část hranice nebo ∆x a

∆y/2 pro nepropustné stěny. Objemy ležící v rozích oblasti pak mají rozměry ∆x/2 a ∆y/2.

Okrajové podmínky jsou opět rozlišeny podle toho, na kterých hranách oblasti daný objem leží. Jak je z obrázku 2.3 vidět, s touto sítí předepisujeme osm tvarů okrajových podmínek v různých částech oblasti, které jsou číselně rozlišeny (1, 2, 3, 4, 13, 14, 23, 24).

Obr. 2.3. Příklad sítě použité pro metodu konečných objemů

2.4.2 Lax-Friedrichsovo schéma nahrazení numerického toku H

Pro nahrazení numerického toku H v rovnici (2.51) je použito Lax-Friedrichsova schématu ve tvaru viz [5]