Transportní procesy. Učební text. 23. září Milan Hokr. Fakulta mechatroniky a mezioborových inženýrských studií Technická univerzita v Liberci

Transportní procesy Učební text

Milan Hokr

23. září 2005 Fakulta mechatroniky a mezioborových inženýrských studií Technická univerzita v Liberci

OBSAH Předmluva . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

1. Úvod: fyzikální a matematický charakter konvekčně-difuzních procesů .

5

2. Transportní procesy jako zákony zachování ve fyzice . . . . . . . . 2.1 Vyjádření základních transportních procesů . . . . . . . . . . . 2.2 Konvekční a difuzní toky – odvození . . . . . . . . . . . . . . . 2.3 Převod transportních rovnic do bezrozměrné formy . . . . . . . 2.4 Okrajové podmínky pro úlohy transportu . . . . . . . . . . . . . 2.5 Cvičení: Řešení jednoduchých úloh integrálními transformacemi 2.5.1 Označení a důležité vztahy . . . . . . . . . . . . . . . . . 2.5.2 Difuzní rovnice (vedení tepla) v 1D nekonečné oblasti . . 2.5.3 Difuzní rovnice v polonekonečné oblasti . . . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

7 8 9 10 11 12 12 13 16

3. Numerické metody pro konvekčně-difuzní rovnici . . . . . . . . . . . 3.1 Analytická řešení za speciálních podmínek . . . . . . . . . . . . 3.1.1 Rovnice konvekce . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Rovnice konvekce-difuze . . . . . . . . . . . . . . . . . . . 3.2 Metoda konečných diferencí . . . . . . . . . . . . . . . . . . . . 3.2.1 Základní charakteristika . . . . . . . . . . . . . . . . . . . 3.2.2 Diferenční vzorce a schémata . . . . . . . . . . . . . . . . 3.2.3 Numerické vlastnosti – konzistence, stabilita, konvergence . 3.2.4 Von Neumannova metoda vyšetřování stability . . . . . . 3.3 Vlastnosti schémat MKD pro konvekčně-difuzní rovnici . . . . . 3.3.1 Vlastnosti jednoduchých explicitních schémat . . . . . . . 3.3.2 Implicitní schémata . . . . . . . . . . . . . . . . . . . . . . 3.3.3 Obecné parametrické schéma . . . . . . . . . . . . . . . . . 3.3.4 Parametrizace pomocí „uměléÿ difuze . . . . . . . . . . . . 3.4 Formulace metody konečných prvků pro konvekčně-difuzní rovnici 3.4.1 Variační formulace a metoda konečných prvků . . . . . . . 3.4.2 Diskrétní formulace v 1D s lineárními bázovými funkcemi . 3.4.3 Petrov-Galerikova metoda – upwinding . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . .

18 19 19 19 20 20 21 23 26 27 27 32 33 34 36 36 37 38

Obsah

3.5 Metoda rozkladu operátoru . . . . . . . . . . . . . . . . . 3.6 Cvičení: Výpočet podmínek stability Fourierovou metodou 3.6.1 Centrální schéma pro konvekční rovnici . . . . . . 3.6.2 Upwind schéma pro konvekční rovnici . . . . . . . 3.6.3 Centrální schéma pro difuzní rovnici . . . . . . . . 3.6.4 Upwind schéma pro konvekčně-difuzní rovnici . . 3.6.5 Centrální schéma pro konvekčně-difuzní rovnici . . 3.7 Cvičení: Demonstrační numerické výpočty . . . . . . . . 3.7.1 Rovnice konvekce . . . . . . . . . . . . . . . . . . . 3.7.2 Konvekčně-difuzní rovnice . . . . . . . . . . . . . 3.7.3 Použití interaktivního programu JBONE . . . . . .

3

. . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

39 42 42 42 44 44 45 46 47 48 49

. . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

53 53 54 54 55 56 57 58 59 61 62 66 66

5. Nelineární a vícerozměrné transportní úlohy . . . . . . . . . . . 5.1 Eulerovy a Navier-Stokesovy rovnice v konzervativním tvaru 5.2 Převod systému rovnic na „kanonickýÿ tvar . . . . . . . . . 5.3 Nalezení řešení rovnice pomocí charakteristik . . . . . . . . . 5.3.1 Rázová vlna . . . . . . . . . . . . . . . . . . . . . . . 5.4 Model dopravního proudu . . . . . . . . . . . . . . . . . . . 5.5 Příklady úpravy systému rovnic na kanonický tvar . . . . . . 5.5.1 Vlnová rovnice . . . . . . . . . . . . . . . . . . . . . 5.5.2 Rovnice mělké vody . . . . . . . . . . . . . . . . . . . 5.6 Konvekčně-reakčně-difuzní rovnice . . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

73 74 75 77 77 79 84 84 86 87

Doporučená literatura . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89

4. Transportní procesy v porézním prostředí . . . . . . . . . . . 4.1 Definice spojitého popisu porézního prostředí . . . . . . . 4.2 Proudění kapaliny v porézním prostředí . . . . . . . . . . 4.2.1 Darcyho experiment a základní veličiny . . . . . . 4.2.2 Darcyho zákon ve 3D . . . . . . . . . . . . . . . . . 4.2.3 Rovnice bilance hmoty . . . . . . . . . . . . . . . . 4.2.4 Okrajové a počáteční podmínky . . . . . . . . . . 4.3 Transport rozpuštěných látek . . . . . . . . . . . . . . . . 4.3.1 Advekce a disperze . . . . . . . . . . . . . . . . . . 4.3.2 Počáteční a okrajové podmínky . . . . . . . . . . 4.3.3 Sorpce . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.4 Vliv chemických a dalších reakcí . . . . . . . . . . 4.4 Motivace: aplikace modelů na úlohy podzemního proudění

PŘEDMLUVA Tento dokument slouží jako studijní materiál pro předmět „Transportní procesyÿ na Fakultě mechatroniky a mezioborových inženýrských studií Technické univerzity v Liberci. Učební text byl průbežně doplňován během prvních let výuky předmětu a zatím existuje jen v elektronické verzi, která je k dispozici ke stažení na stránkách Katedry modelování procesů. V současné době text pokrývá zhruba 90% probírané látky. Autor uvítá jakékoli připomínky a náměty k doplnění textu.

Kapitola 1

ÚVOD: FYZIKÁLNÍ A MATEMATICKÝ CHARAKTER KONVEKČNĚ-DIFUZNÍCH PROCESŮ Spojujícím tématem jednotlivých kapitol a nosným tématem předmětu jsou tzv. konvekčnědifuzní procesy. Typickými příklady jsou transport rozpuštěné látky nebo tepla v pohybující se kapalině a proudění viskózní tekutiny. S určitou mírou abstrakce lze tyto procesy vyjádřit jednotně tzv. konvekčně-difuzní rovnicí ∂u − ∇ · (D∇u) + ∇ · (vu) = 0 ∂t

(1.1)

což je parciální diferenciální rovnice 2.řádu (u(x, t) je neznámá, D a v jsou parametry). Náplní předmětu (a obsahem tohoto dokumentu) jsou odvození a popis rovnice pro různé fyzikální případy, její zobecnění na složitější jevy a matematické metody řešení (analytické a numerické). Mezi parciálními diferenciálními rovnicemi (PDR) jde o důležitou speciální třidu úloh – vyjadřují konkrétní úlohy reálné praxe a jejich numerické řešení vyžaduje speciální techniky. Příčina obtíží je v přítomnosti dvou procesů odlišného fyzikálního a matematického charakteru – konvekce a difuze, první popsaný hyperbolickou PDR 1. řádu, druhý popsaný parabolickou PDR 2. řádu. Rozdílný charakter obou procesů a vlastnosti řešení řídící rovnice si můžeme snadno demonstrovat na jednoduché úloze transportu rozpuštěné látky v proudící vodě v 1D. Mechanismus tohoto procesu si přirozeně představíme tak, že látka je unášena spolu s proudící vodou (tzv. konvekce, funkce popisující rozložení látky v prostoru se jen „posouváÿ v čase) a zároveň dochází k přenosu látky z míst s vyšší koncentrací do míst s nižší koncentrací (tzv. difuze, funkce popisující rozložení látky v prostoru se „rozmazáváÿ, snižuje se strmost a maxima a minima se přibližují). Každému z procesů odpovídá jeden ze členů rovnice (1.1), samostná konvekce je popsána ∂u − ∇ · (vu) = 0 (hyperbolická 1.řádu), samotná difuze ∂u − ∇ · (D∇u) = 0 ∂t ∂t (parabolická 2.řádu). Poznamenáme i rozdíl z pohledu termodynamického: konvekce je proces vratný (rovnice má smysl i po obrácení směru toku času), difuze je proces

Kapitola 1. Úvod: fyzikální a matematický charakter konvekčně-difuzních procesů

6

nevratný (rovnice po obrácení směru času, ev. znaménka koeficientu D není „rozumnáÿ – nelze formulovat korektně definovanou úlohu). Numerické řešení se stává obtížné v případě, kdy konvekce je dominantní nad difuzí, což se vyjadřuje bezrozměrným Pécletovým číslem Pe = vL , kde L je charakteristický D rozměr úlohy (v metrech). Velké Pe znamená významnější konvekci, malé Pe významnější difuzi. Typickými nežádoucími jevy při numerickém výpočtu jsou nefyzikální oscilace a tzv. numerické difuze (ve výsledku nuemrického řešení se projevuje větší difuze než jaká odpovídá koeficientu v rovnici). Další komplikované situace jsou spojeny s nelineárními členy v rovnici – tj. pokud např. materiálové koeficienty závisí na neznámé veličině. V některých případech je chování řešení kvalitativně odlišné a to je nutno vzít v úvahu při formulaci problému i při numerickém výpočtu.

Kapitola 2

TRANSPORTNÍ PROCESY JAKO ZÁKONY ZACHOVÁNÍ VE FYZICE Transportním procesem rozumíme děj, při kterém dochází k přesunu a transformaci nějaké veličiny v prostoru a čase. Běžnými případy jsou přenos1 tepla (vedení tepla, pohyb kapaliny s nehomogenní teplotou) – transportovanou veličinou je teplota a šíření rozpuštěné látky (difuze, konvekce) – transportovanou veličinou je koncentrace látky v kapalině. V závislosti na řídícím fyzikálním mechanismu se popis procesu bude lišit, ale vždy je založen na vyjádření bilance transportované veličiny v určitém místě prostoru. Rovnice bilance má pro každý typ procesu obdobný tvar, ale obecně může být různá pohybová rovnice (tj. závislost toku veličiny na jejím rozložení v prostoru) – ta už je daná konkrétním typem procesu (konvekce, difuze, . . . ). Jednotně můžeme vyjádřit bilanci veličiny X takto2 : akumulace X ve V = přítok X do V − odtok X z V + produkce X ve V dt

(2.1)

přičemž rozměr rovnice je X/s. Pro hustotu toku veličiny X zavedeme označení ΦX (rozměr X/m2 /s). Níže vyjádříme bilanci v přesné matematické formě (integrální a diferenciální), transportovanou veličinu X vyjádříme ve formě objemové hustoty XV (jednotka X/m3 ). pøítok

odtok

X objem V

1

Termín přenos by byl plně český ekvivalent termínu transport, ve starší literatuře se používá označení „přenosové procesy/jevyÿ. 2 Označení X bude použito pro veličinu samotnou i její jednotku, význam je však vždy ze souvislosti zřejmý. Důvodem je přizpůsobení značení použitého dle [2] ostatnímu textu.

Kapitola 2. Transportní procesy jako zákony zachování ve fyzice

2.1

8

Vyjádření základních transportních procesů

Ukážeme, jak jednoduché procesy transportu odpovídají základním zákonům zachování ve fyzice (ZZ hmoty, energie a hybnosti) a lze je popsat strukturálně identickými matematickými vztahy (pohybové rovnice a koeficienty). Všechny vztahy jsou „difuzníhoÿ charakteru – tok veličiny je úměrný gradientu veličiny ΦX = konst ·

∂XV . ∂x

(2.2)

přičemž rozměr konstanty je m2 /s. Vztahy po řadě popisují transport rozpuštěné látky (vyjádřený koncentrací: hmotnost/objem), transport tepla (vyjádřeno hustotou energie: energie/objem, přičemž jako stavovou veličinu lze energii nahradit teplotou) a transport hybnosti, tj. vazkost proudící kapaliny. Poslední případ je poněkud méně názorný z toho důvodu, že bilancované/transportovaná veličina je vektorová a je nutno uvažovat speciálně definovanou úlohu pro to, aby transport hybnosti byl řízen pouze uvedeným vztahem a ne ještě konvekcí. ∂c ∂x ∂(%cp T ) ZZ energie / vedení tepla ΦE = −a ∂x ∂(%v) ZZ hybnosti / viskozita Φp = −ν ∂x ZZ hmoty / difuze látky

Φm = −D

[kg/m2 /s = m2 s · kg/m3 /m] [J/m2 /s] [kgm−2 s−1 /m3 /s = Pa]

Význam značek je následující: c T v % cp D a ν

koncentrace teplota rychlost hustota tepelná kapacita koeficient molekulární difuze koef. teplotní vodivosti kinematická viskozita

Vidíme analogii jednak stavových veličin („něcoÿ na jednotkový objem – koncentrace c, hustota energie %cp T a hustota hybnosti %v). Dále je analogie mezi koeficienty: všechny mají rozměr m2 /s. Druhou a třetí rovnici je obvyklé formulovat též s jinou „stavovouÿ veličinou a jiným koeficientem: ∂v ∂T Φp = −η (2.3) ΦE = −λ ∂x ∂x


9

kde λ = %cp a [W/m/K] se nazývá tepelná vodivost (vs. teplotní vodivost a; anglicky heat conductivity λ vs. thermal diffusivity a) a η = %ν je dynamická viskozita. Transport hybnosti difuzního charakteru (přenos ve směru gradientu) lze názorně vyjádřit pro proudění kapaliny s polem rychlosti se stejným směrem a různou velikostí (např. laminární proudění v trubce uprostřed a na kraji). Transport hybnosti pak probíhá ve směru příčném (kolmém na směr rychlosti), přičemž tok hybnosti je vlastně smykové napětí (rozměr rovnice je Pa!). Uvedený vztah platí pro tzv. Newtonovskou kapalinu (viz předmět mechanika tekutin a příklad ve cvičení - kvadratické rozložení teploty resp. rychlosti napříč trubky).

2.2

Konvekční a difuzní toky – odvození

Z diskrétního vyjádření zachování veličiny X rovnice (2.1) odvodíme diferenciální rovnici [bude doplněn podrobný výpočet] ∂XV + ∇ · ΦX − rX = 0 (2.4) ∂t V této podobě platí rovnice pro libovolný mechanismus transportu (nejen lineární konvekce a difuze), na jejím základě budeme vyjadřovat nelineární rovnice v kapitole 5. Dva základní typické mechanismy transportu, podstatně odlišného charakteru, jsou konvekce a difuze. Pro konvekci3 (přenos veličiny spolu s nosným médiem – např. rozpuštěná látka s proudící vodou) platí Φkonv = vXV X

(2.5)

a pro difuzi (rozptyl – přenos ve směru gradientu hustoty) Φdif X = ∇XV

(2.6)

Po dosazení do rovnice bilance (celkový tok jako součet difuzního a konvekčního toku) dostaneme konvekčně-difuzní rovnici: ∂XV + ∇ · (XV v) − ∇ · (D∇XV ) − rX = 0 (2.7) ∂t Rovnice je parabolického typu (obsahuje první derivaci podle času a druhou derivaci podle prostoru). Přítomnost první derivace („konvektivníÿ člen ∇ · (XV v)) však hraje v rovnici podstatnou roli: pokud je difuzní člen s druhou derivací ∇ · (D∇XV ) malý, rovnice se vlastnostmi blíží k hyperbolickému typu (pouze první derivace). To je důležité např. při formulaci okrajových podmínek a při numerickém řešení rovnice. Poměr významnosti jednotlivých členů rovnice, a tedy celkový charakter fyzikálního děje, lze vyjádřit pomocí tzv. podobnostních čísel (bezrozměrných konstant). 3

Též advekce, termín konvekce se někdy užívá v užším smyslu (např. stoupání teplého vzduchu v atmosféře).


Cµ plyn 2

10

difúze konvekce

plyn 1 u x = konst.

x

Obrázek 2.1:

2.3

Převod transportních rovnic do bezrozměrné formy

jednotlivé příklady (cvičení) – ukázky definice podobnostních čísel Pecletovo, Reynoldsovo – poměr konvekce/difuze v přenosu tepla: Fourierovo (pronikání tepla do média), Nusseltovo (vedení tepla vs. přestup tepla)

Poměr konvekce/difuze – Pécletovo číslo Rovnice (2.7) pro koncentraci v bezrozměrné formě: ∂C v 1 + ∇ · (C ) − ∇ · ( ∇C) = 0 ∂T |v| Pe

(2.8)

kde

x |v|t c vL T = C= Pe = . L L C0 D Zavedli jsme charakteristickou délku L, charakteristickou koncentraci C0 a jedno bezrozměrné podobnostní číslo: Pécletovo Pe. X=

Příklad: polonekonečná trubice ústící do nádrže (interval (0, +∞)) – obrázek 2.1. Trubicí proudí čistá voda (směrem „z nekonečnaÿ do nádrže) rychlostí v a z nádrže s udržovanou konstantní koncentrací látky c∞ se látka šíří difuzí (koeficient D) proti směru proudu vody a konvekcí je transportována zpět. Dojde k vytvoření ustáleného stavu. Rozložení látky v ustáleném stavu je popsáno rovnicí (nulová časová derivace) 0 = −D Řešení rovnice je

∂c − vc ∂x

vx c = e− D = e−Pe c∞

(2.9)


11

Bezrozměrné Pécletovo číslo dává do souvislosti intenzitu konvekce a difuze a měřítko pozorování (souřadnice x – charakteristické délka). Názorně: rychlost poklesu koncentrace je dána poměrem konvekce a difuze, ale zároveň délkovým rozměrem resp. měřítkem pozorování (podle toho jestli uvažujeme nekonečnou oblast nebo konkrétní těleso) – ukázka na obrázku.

Vliv vazkosti kapaliny – Reynoldsovo číslo Reynoldsovo číslo je obdobou Pécletova čísla (poměr mezi „konvekčnímÿ a „difuznímÿ členem) pro Navier-Stokesovy rovnice (tj. konvekčně-difuzní transport hybnosti). Viz předmět „Mechanika tekutinÿ. vd Re = ν kde d je charakteristický rozměr, např. průměr trubky s proudící vodou nebo velikost obtékaného tělesa. v2 Při zahrnutí vlivu gravitace se uplatňuje Froudovo číslo Fr = dg

2.4

Okrajové podmínky pro úlohy transportu

Fyzikální význam okrajových podmínek (zadána hodnota nebo tok), příklady okr.p. v úlohách pro různé procesy – vedení tepla, přenos látky.

Dirichletova podmínka Zadaná hodnota transportované veličiny: • teplo: dokonalé chlazení/zahřívání – pevná teplota na hranici • zadaná koncentrace a teplota ve vodě přitékající do řešené oblasti (při dominantní konvekcí) – „vstupÿ do systému

Neumannova podmínka Zadána derivace veličiny na hranici, tj. tok: Typický příklad homogenní podmínka (nulový tok) – izolovaná hranice.

Newtonova/Cauchyova podmínka Typicky jde o vyjádření interakce: tok je funkcí hodnoty veličiny: • tok přes polopropustnou vrstvu na hranici je úměrný rozdílu veličiny uvnitř (neznámá fce) a vně (zadáno) – pro teplo i rozpuštěnou látku


2.5

12

Cvičení: Řešení jednoduchých úloh integrálními transformacemi

K řešení okrajových a počátečních úloh pro parciální diferenciální rovnice (PDR) lze použít integrální transformace, např. Laplaceovu nebo Fourierovu. Tyto transformace zobrazují funkce (určité třídy) na jiné funkce, přičemž operace derivování se transformují na algebraické operace. Tím v případě obyčejných diferenciálních rovnic převedeme úlohu na řešení algebraické rovnice a v případě parciální rovnice dostaneme po transformaci obyčejnou diferenciální rovnici nebo PDR o jedno nižší dimenze. Ve výkladu předpokládáme, že čtenář má základní znalosti o použití zmíněných transformací (na FM v předmětu Aplikovaná matematika). Transformaci je třeba volit v závislosti na intervalu proměnné v dané úloze tak, aby odpovídal definici transformace. Typicky se Fourierova transformace použije pro prostorové proměnné (interval (−∞, +∞)) a Laplaceova transformace pro časovou proměnnou (interval (0, +∞)).

2.5.1

Označení a důležité vztahy

Derivace: ut =

∂u ∂t

uxx =

∂2u ∂x2

etc.

(2.10)

Funkce normálního rozdělení pravděpodobnosti Používají se následující označení funkce erf(x) a erfc(x) (chybová funkce a komplementární chybová funkce), až na aditivní a multiplikativní konstanty jde o distribuční funkci normálního rozdělení ( erf(x)+1 = Fnorm (x)). 2 2 erfc(x) = √ π

Z∞ x

−u2

e

2 du = 1 − √ π

Zx 0

2

e−u du = 1 − erf(x)

(2.11)

Laplaceova transformace L[f (t)] =

Z∞

f (t)e−st dt = F (s)

0

Operace: linearita derivace v t jiná derivace

L[c1 f1 (t) + c2 f2 (t)] = c1 F1 (s) + c2 F2 (s) L[f 0 (t)] ≡ L[ft (t)] = s · F (s) + f (0+) L[fx (x, t)] = (F (x, s))x

(2.12)


13

Fourierova transformace 1 F [f (x)] = √ 2π Operace: linearita derivace v x jiná derivace

Z∞

f (x)e−iωx dx = F (ω)

(2.13)

−∞

F [c1 f1 (x) + c2 f2 (x)] = c1 F1 (ω) + c2 F2 (ω) F [f 0 (x)] ≡ F [fx (x)] = iω · F (ω) F [ft (x, t)] = (F (ω, t))t

Zobrazení konvoluce:   Z∞ 1 F √ f (x − ξ)g(ξ)dξ  = F [f (x) ∗ g(x)] = F (ω) · G(ω) 2π

(2.14)

−∞

2.5.2

Difuzní rovnice (vedení tepla) v 1D nekonečné oblasti

Řešíme rovnici ut = Duxx

(2.15)

pro neznámou funkci u(x, t) v oblasti − ∞ < x < +∞

t>0

(2.16)

s počáteční podmínkou u(x, 0) = δ(x) ,

(2.17)

kde δ(x) je Diracova funkce. Úlohu transformujeme Fourierovou transformací v proměnné x, obraz funkce označíme F [u(x, t)] = U(ω, t) (proměnná t má zde význam parametru). Z diferenciální rovnice (2.15) dostaneme dU(ω, t) = −D · ω 2 U(ω, t) (2.18) dt a z počáteční podmínky (2.17) dostaneme 1 U(ω, 0) = √ 2π

(2.19)

tj. počáteční podmínku pro (nyní již obyčejnou) diferenciální rovnici (2.18). Převedli jsme tedy PDR s proměnnými x a t na ODE s proměnnou t (v tomto okamžiku uvažujeme t jako proměnnou a ω jako parametr).


14

Řešení (2.18) s počáteční podmínkou (2.19) určíme např. metodou separace proměnných: 1 2 U(ω, t) = √ · e−Dω ·t (2.20) 2π Nyní toto řešení převedeme zpětnou transformací do původní proměnné x (opět se mění význam – t je parametr a ω je proměnná). S využitím převodního vztahu 2 −1 − ω 2 (F ) 1 √ 4a → e a 2

2 x2

e−a

(zde a =

√1 ) 2 Dt

dostaneme

x2 1 1 u(x, t) = √ √ e− 4Dt . 2 π Dt

(2.21)

Zobecněním této úlohy bude: • obecná počáteční podmínka • zahrnutí konvekce (člen vux v rovnici) Obecná počáteční podmínka Uvažujeme nyní tutéž rovnici (2.15) pro x ∈ R, t ∈ R+ , ale s obecnou počáteční podmínkou, danou funkcí φ(x), tj. u(x, 0) = φ(x),

x∈R

(2.22)

Postupujeme analogicky zobrazením diferenciální rovnice i počáteční podmínky Fourierovou transformací v proměnné x a dostáváme dU(ω, t) = −D · ω 2 U(ω, t) dt

(2.23)

a z počáteční podmínky dostaneme U(ω, 0) = Φ(ω)

(2.24)

kde Φ = F (φ), přičemž uvedenou operaci lze použít pouze pro takovou počáteční podmínku, pro níž existuje Fourierův obraz (tj. např. už ne pro Heavisideovu skokovou funkci). Řešením této úlohy je funkce (z pohledu řešení ODE jde o tutéž situaci jako u (2.19) - počáteční podmínka je konstanta, pouze závislá na parametru) U(ω, t) = Φ(ω) · e−Dω

2 ·t

(2.25)


15

Zpětná transformace však již je složitější, neboť nyní je ω proměnná a jde tedy o součin dvou funkcí, který se zpětnou transformací převádí na operaci konvoluce, tj. 2

u(x, t) = F −1 [Φ(ω) · e−Dω ·t ] = x2 1 1 = φ(x) ∗ √ √ e− 4Dt = 2 Dt +∞ Z ξ2 1 = √ φ(x − ξ)e− 4Dt dξ 2 πDt

(2.26)

−∞

Tím jsme odvodili řešení v obecném tvaru, po kontrolním dosazení Diracovy funkce skutečně dostáváme totéž řešení jako v předchozí úloze (integrál se redukuje na bodovou hodnotu). Zajímavé je si všimnout, že integrál (2.26) existuje pro značně širší třídu funkcí φ(x) než pro které existuje Fourierův obraz (tj. pro které bylo toto řešení odvozeno). Vzniká tedy otázka, je-li po dosazení takové funkce φ (která nemá obraz), výraz (2.26) skutečně řešením původní úlohy. Přesnou odpověď dává funkcionální analýza, v jednoduchých případech se lze o vlastnostech výsledné funkce přesvědčit dosazením do rovnice (derivováním). Odvozené řešení je v souladu s principem superpozice – integrál konvoluce vyjadřuje součet příspěvků od více bodových zdrojů (Dirac, ev. spojitě rozložených) v počáteční podmínce (vzhledem k řešení úlohy s Diracovou počáteční podmínkou. Příkladem řešení pro funkci φ, která nemá Fourierův obraz je již zmíněná Heavisideova funkce. Např. pro počáteční podmínku 1 x<0 u(x, 0) = −H(x) = (2.27) 0 x>0 dostáváme řešení 1 u(x, t) = √ 2 πDt

Z+∞ Z+∞ ξ2 1 1 2 − 4Dt dξ = e−ξ dξ = erfc(x) 1·e 2 2 x

(2.28)

x

Rovnice s konvekčním členem K řešení je možno buď použít transformací souřadnic (náhrada x−vt, převedení na čistě difuzní rovnici) nebo rovněž přímo Fourierovou transformací. Ukážeme stručně rozdíl oproti difuzní rovnici, pro Diracovu počáteční podmínku (pro obecnou je to totéž). Zobrazená rovnice má tvar dU(ω, t) = (−ivω − D · ω 2 )U(ω, t) dt

(2.29)


16

a řešení

1 2 U(ω, t) = √ · e−Dω ·t · e−ivωt 2π −iaω Dle pravidla F [f (x − a)] = e · F (ω) dostáváme (x−vt)2 1 u(x, t) = √ e− 4Dt , 2 πDt

(2.30)

(2.31)

tj. vidíme rozdíl v posunu „po charakteristiceÿ (zde není tento termín přesný, jde o parabolickou rovnici a tvar funkce se mění).

2.5.3

Difuzní rovnice v polonekonečné oblasti

Tvar oblasti proměnných u PDR má podstatný vliv na řešení. V případě použití integrálních transformací znamená přizpůsobit volbu transformace a proměnné, eventuálně s přihlédnutím ke konkrétní počáteční a okrajové podmínce (transformovatelnost). V případě difuzní rovnice v polonekonečném intervalu v „prostoruÿ použijeme Laplaceovu transformaci v proměnné t, transformovaná úloha pak bude ODE 2. řádu s proměnnou x (na rozdíl od předchozích úloh s Fourierovou transformací, kde výsledkem byla ODE 1. řádu s proměnnou t). Roli potřebných podmínek (v tomto případě okrajových) bude hrát obraz okrajové podmínky v budě x = 0 a chování funkce v nekonečnu (tedy dvě podmínky). Nyní podrobně. Řešíme tedy rovnici 2.15 pro x ∈ R+ a t ∈ R+ s podmínkami u(0, t) = 1

u(x, 0) = 0

(2.32)

tj. např. na počátku studená tyč ohřívaná od konce a jinak izolovaná. Laplaceovou transformací L[u(x, t)] = U(x, s) dostáváme z diferenciální rovnice s · U(x, s) − u(x, 0+) = D

dU(x, s) dx2

(2.33)

a z okrajové podmínky U(0, s) = 1/s

(2.34)

Díky nulové počáteční podmínce (nyní zahrnuta do transformované rovnice) jde o homogenní rovnici 2.řádu v proměnné x, s konstantní počáteční podmínkou (s je nyní parametr). Obecné řešení tedy je r r s s U(x, s) = k1 exp( x) + k2 exp(− x) (2.35) D D a konstanty k1 , k2 určíme z podmínek U(0, s) = k1 + k2 = 1/s

U(∞, s) = k2 ∞ = [konečné číslo]

(2.36)


tj. k1 = 0 a k2 = 1/s. Provedeme nyní zpětnou transformaci (v proměnné s) x √ x −1 1 −( √D ) s L e = erfc( √ ) s 2 Dt

17

(2.37)

Kapitola 3

NUMERICKÉ METODY PRO KONVEKČNĚ-DIFUZNÍ ROVNICI Zabývat se budeme řešením konvekčně-difuzní rovnice, popisující různé typy fyzikálních transportních procesů jak bylo zmíněno v předchozí kapitole – tj. např. šíření rozpuštěné látky a tepla v proudící kapalině. Nejprve jsou uvedena analytická řešení pro vybrané úlohy se speciálními okrajovými a počátečními podmínkami, jejichž znalost je dobrým vodítkem při vyhodnocování výsledků numerických výpočtů. Hlavním obsahem textu pak je popis základní metody pro řešení parciálních diferenciálních rovnic, metody konečných diferencí (MKD), a prozkoumání jejích vlastností v závislosti na použité numerické variantě a hodnotě koeficientů rovnice. Ve většině případů je uvažována jednorozměrná úloha a rovnice s konstantními koeficienty, kdy je možno popisované jevy názorně demonstrovat. Pro úlohy vyšší dimenze je pak možné příslušné úvahy provést obdobně. Vyjdeme z obecné rovnice pro konvekčně-difuzní transport látky odvozené v předchozí kapitole (rovnice (2.7)), kterou zapíšeme pro jednorozměrnou oblast (úsečku nebo přímku) a jako neznámou funkci uvažujeme koncentraci látky c(x, t). Budeme uvažovat ustálené proudění s rychlostí v = konst, konstantní difuzní koeficient D, žádné chemické reakce ani jiné zdroje/propady. ∂c ∂c ∂2c = −v +D 2 ∂t ∂x ∂x V matematicky zaměřených textech se zavádí jiné značení, např. a posuvná rychlost, b difuzní koeficient a u neznámá funkce. Píšeme tedy ∂u ∂u ∂2u +a = b 2, ∂t ∂x ∂x

(3.1)

kde klademe b ≥ 0 (z fyzikální podstaty difúze-difuze) a a ≥ 0, což odpovídá proudění v kladném smyslu osy x. Jak už je ze vzniku a popisu koeficientů a, b patrné, člen ∂u ∂x 2 charakterizuje konvekci a člen ∂∂xu2 difúzi.

Kapitola 3. Numerické metody pro konvekčně-difuzní rovnici

19

Rovnice (3.1) je obecně parabolického typu, je však nutné vzít speciálně v úvahu případ b = 0, kdy člen s druhou derivací (difuze) není přítomen a rovnice je tedy typu hyperbolického. Z toho důvodu je obvykle nutné tento případ zkoumat zvlášť: zadávají se okrajové podmínky v jiné struktuře a numerické metody mohou mít odlišné vlastnosti.

3.1 3.1.1

Analytická řešení za speciálních podmínek Rovnice konvekce

Uvažujeme rovnici (3.1) pro speciální případ b = 0, tj, s nulovým difuzním členem ∂u ∂u +a =0 ∂t ∂x

(3.2)

Úloha na nekonečné oblasti Ω = (−∞ , ∞) s počáteční podmínkou u(x, 0) = u0(x) x ∈ Ω

(3.3)

má jednoduché řešení resp.

u(x, t) = u0 (x − at) x ∈ Ω , t > 0 u = konst podél přímek x − at = konst

(3.4) (3.5)

což znamená, že určitý signál se posouvá z bodu x = 0 rychlostí a v kladném smyslu osy x. Řešení lze snadno zobecnit i na úlohy na omezeném intervalu se zadanou okrajovou podmínkou na jednom z konců, v závislosti na znaménku rychlosti a (podmínka se zadává v bodě „vtokuÿ).

3.1.2

Rovnice konvekce-difuze

Uvedeme opět řešení úlohy na nekonečné oblasti Ω = (−∞ , ∞), tj. bez zadání okrajových podmínek. Budeme uvažovat dvě základní okrajové podmínky, pro které lze zapsat řešení v přehledném tvaru. Počáteční podmínka ve tvaru Diracovy funkce Je-li uvažována fyzikálně idealizovaná situace, kdy na počátku je veškerá hmota soustředěna v jednom bodě a ve zbytku oblasti je koncentrace nulová, lze tuto situaci zapsat M δ(x) (3.6) u(x, 0) = S


20

kde M je celková hmota, vhodně normalizovaná, např. příčným průřezem S, a δ(x) je Diracova δ-funkce. Řešení této úlohy pro všechny body x ∈ Ω a libovolný čas t > 0 je h (x − at)2 i M √ exp − u(x, t) = 4bt 2S πbt

(3.7)

což je funkce hustoty pravděpodobnosti gaussovského rozdělení. Počáteční podmínka ve tvaru skokové funkce Pro počáteční podmínku ve tvaru u(x, 0) = u0 pro x < 0 u(x, 0) = 0 pro x > 0 je řešení ve tvaru distribuční funkce gaussovského rozdělení h x − at i u0 u(x, t) = erfc √ x ∈ Ω, t > 0 2 2 bt kde

2 erfc(x) = √ π

Z∞

−u2

e

2 du = 1 − √ π

x

Zx

2

e−u du = 1 − erf(x)

(3.8) (3.9)

(3.10)

(3.11)

0

je tzv. komplementární chybová funkce.

3.2 3.2.1

Metoda konečných diferencí Základní charakteristika

Metoda konečných diferencí (MKD), někdy též metoda sítí, spočívá v nahrazení derivací diferencemi , které aproximují derivace pomocí hodnot hledané funkce v několika blízkých bodech. V oblasti Ω × (0, T ), v níž hledáme řešení, si volíme síť konečného počtu bodů a výsledkem metody pak jsou přibližné hodnoty hledané funkce v těchto bodech. Na rozdíl od jiných metod (metoda konečných objemů, metoda konečných prvků) je v případě MKD kladen požadavek tzv. strukturovanosti sítě, což znamená, že body sítě v Rn leží na průsečících křivek sdružených do n skupin přičemž každý bod musí ležet na právě jedné křivce z každé skupiny. Nejčastěji se pak používá pravoúhlá síť, kdy za jednotlivé skupiny volíme systémy rovnoběžek s příslušnou souřadnou osou. Strukturovanost sítě umožňuje nahradit skutečné souřadnice n celočíselnými indexy, udávajícími pořadí průsečíku na příslušné přímce (křivce).


21

Forward difference

y Backward difference

u i+1 ui

y = u(x)

u i-1

Central difference

Dx i-1

Dx i

x i+1

Obrázek 3.1: Grafické znázornění aproximace derivace pomocí diferenčních vzorců – dopředný, centrální, zpětný. V dalším textu popíšeme varianty použití MKD pro konvekčně-difuzní rovnici (3.1) uvedenou v úvodu. Prostorová oblast je uvažována 1D, diferenční vzorce a způsob sestavení schématu jsou pro vyšší dimenze analogické. Řešení rovnice hledáme v dvourozměrném prostoru (proměnné x, t) a použijeme pravoúhlou síť se stejnými kroky a to ∆t ve směru osy t a ∆x ve směru osy x. Potom můžeme označit uni = u(i∆x, n∆t) (3.12)

3.2.2

Diferenční vzorce a schémata

Jak již bylo řečeno, diferenční vzorce jsou vztahy pro aproximaci derivace pomocí funkčních hodnot v bodech sítě (obrázek 3.1). Lze je odvodit více způsoby, v tomto textu pouze stručně naznačíme postup pomocí Taylorova rozvoje funkce v sousedních bodech sítě. Chceme-li např. aproximovat derivaci ∂u ≡ ux , uvažujeme body i + 1 a i − 1 : ∂x ∆x2 (uxx )i + O(∆x3 ) (3.13) 2 ∆x2 ui−1 = ui − ∆x(ux )i + (uxx )i + O(∆x3 ) (3.14) 2 kde pro přehlednost vynecháváme časový index n, který by byl všude stejný. Vyjádřímeli z (3.13) derivaci (ux )i , dostaneme ui+1 = ui + ∆x(ux )i +

(ux )i (ux )i

ui+1 − ui ∆x O(∆x3 ) = − (uxx )i + ∆x 2 ∆x ui+1 − ui + O(∆x) = ∆x

tj. (3.15)


22

Tento vztah se nazývá dopředná diference. Provedeme-li totéž pro (3.14), dostaneme tzv. zpětnou diferenci ui − ui−1 (ux )i = + O(∆x) (3.16) ∆x Společně se vztahy (3.15) a (3.16) nazývají jednostranné diferenční vzorce. Přesnější aproximaci lze pak získat pomocí centrálního diferenčního vzorce, k němuž dospějeme odečtením (3.14) od (3.13): (ux )i =

ui+1 − ui−1 + O(∆x2 ) 2∆x

(3.17)

Mocnina v posledním členu (O()) se nazývá řád aproximace daného vzorce, uvedené jednostranné diference jsou tedy aproximací prvního řadu a centrální diference aproximací druhého řadu. Jestliže vztahy (3.14) a (3.13) sečteme, dostáváme centrální diferenční vzorec pro 2 druhou derivaci ∂∂xu2 ≡ uxx

ui+1 − 2ui + ui−1 + O(∆x2 ) (3.18) 2 ∆x Z Taylorova rozvoje hodnot ui+2 a ui−2 pak získáme obdobnými úpravami i jednostranné diference ui+2 − 2ui+1 + ui (uxx )i = + O(∆x) (3.19) ∆x2 ui − 2ui−1 + ui−2 (uxx )i = + O(∆x), (3.20) ∆x2 které mají opět o řád nižší přesnost a pro námi uvažovanou rovnici nemají význam. Stejným postupem lze dostat vzorce pro časovou derivaci, např. dopředná diference je un+1 − uni (ut )ni = i + O(∆t), (3.21) ∆t a vede na explicitní schémata (viz níže), kterými se budeme obvykle zabývat. (uxx )i =

Numerické schéma explicitní a implicitní Dosadíme-li nyní do diferenciální rovnice za derivace příslušné diferenční výrazy, dostáváme vztah pro výpočet hodnot v uzlových bodech, označovaný výpočetní nebo numerické schéma. Použijeme-li např. pro rovnici konvekce centrální diferenci v x a dopřednou v t, dostáváme schéma

un+1 i

un − uni−1 un+1 − uni i + a i+1 = 0 tj. ∆t 2∆x σ i = 1, . . . , N − 1 = uni − (uni+1 − uni−1 ) 2

(3.22)


23

∆t kde σ = a ∆x . Jak je ihned vidět, schéma umožňuje na základě znalosti hodnot v čase n pouhým dosazením získat hodnoty v čase n + 1, je tedy explicitní v čase. Jiná situace nastává v případě, kdy se ve schématu vyskytují hodnoty pro novou časovou hladinu současně ve více uzlech (např. při použití zpětné diference pro časovou derivaci v předešlém příkladu). Vztahu pro jednotlivé uzly jsou pak provázané a jde o řešení soustavy lineárních rovnic. Takováto metoda se označuje jako implicitní.

Okrajové podmínky Hodnoty v krajní uzlech i = 0 a i = N je třeba vyjádřit pomocí okrajových podmínek původní diferenciální rovnice. Dirichletova podmínka nám dává přímo požadovanou hodnotu, v případě Neumannovy podmínky je třeba pro i = 0 použít dopřednou diferenci a pro i = N zpětnou diferenci un+1 − un+1 ∂u un+1 − un+1 ∂u N N −1 1 0 = + O(∆x) = + O(∆x) ∂x x=0 ∆x ∂x x=N ∆x

Např. při použití explicitního schématu dostáváme pro homogenní Neumannovu podmínku přímo un+1 = un+1 un+1 = un+1 (3.23) 0 1 N N −1

na základě již známých hodnot un+1 , un+1 1 N −1 .

3.2.3

Numerické vlastnosti – konzistence, stabilita, konvergence

Jak se lze přesvědčit na jednoduchých příkladech, výpočetní schéma získané dosazením diferenčních vzorců do diferenciální rovnice nemusí automaticky dávat uspokojivé výsledky. Numerické vlastnosti schématu a přesnost řešení se obvykle vyjadřují pomocí pojmů konzistence, stabilita a konvergence. Konzistence Pojem konzistence charakterizuje vztah mezi diferenčním schématem a řešenou diferenciální rovnicí. Lze jej vyjádřit podmínkou, že při zjemňování sítě (∆x, ∆t → 0) se hodnota výrazu vyskytujícího se v diferenčním schématu blíží hodnotě odpovídajícího výrazu v diferenciální rovnici. Tuto skutečnost je možno ověřit dosazením Taylorova rozvoje v bodě (i, n) za ostatní členy ve schématu a to bez ohledu na to, jakým způsobem bylo schéma odvozeno. Např. pro výše uvedené schéma (3.22) vyjádříme un+1 = uni + ∆x(ut )ni + i

∆t3 ∆t2 (utt )ni + (uttt )ni + O(∆t4 ) 2 6

(3.24)


uni+1 uni−1

∆x2 ∆x3 n = + + (uxx )i + (uxxx )ni + O(∆x4 ) 2 6 3 2 ∆x ∆x (uxx )ni − (uxxx )ni + O(∆x4 ) = uni − ∆x(ux )ni + 2 6 uni

∆x(ux )ni

24

(3.25) (3.26)

a po dosazení dostaneme uni+1 − uni−1 un+1 − uni i +a − (ui + aux )ni = ∆t ∆x ∆x2 ∆t (utt )ni + a(uxxx )ni + O(∆t2 , ∆x4 ) = 2 6

(3.27)

kde na levé straně první dva členy vyjadřují diferenční schéma a třetí pak řešenou diferenciální rovnici. Výraz na pravé straně pak pro ∆x, ∆t → 0 jde k nule a schéma je tedy s danou rovnicí konzistentní. Řád konvergence pravé strany udává řád přesnosti použitého schématu, v našem případě máme schéma 1. řádu v čase a 2. řádu v prostoru. Stabilita Vlastnost stability se vztahuje čistě k diferenčnímu schématu, bez zřetele na původní diferenciální rovnici. Vyjadřuje skutečnost, že numerická chyba (např. vlivem zaokrouhlovacích chyb v počítači), vzniklá v určitém místě výpočtu, nesmí v dalších krocích neomezeně růst. Označíme-li u¯ni hodnotu přesného řešení v bodě (i∆x, n∆t), můžeme chybu zapsat εni = uni − u¯ni (3.28) a podmínka stability pak dostává tvar

lim |εni | ≤ K

n→∞

pro pevné ∆t,

(3.29)

kde K je konstanta nezávislá na n. Obecnější definici stability je možno podat na základě chování samotného řešení s postupujícím časem. Podmínku lze stručně vyjádřit tak, že žádná složka počítaného řešení nesmí neomezeně růst. Přesněji můžeme tento princip vyjádřit s použitím maticového zápisu schématu. Označíme-li U n vektor řešení v čase n∆t  n  u1  ..   .  Un ≡  n   ui  .. . lze libovolné explicitní schéma zapsat ve formě

U n+1 = C · U n

(3.30)


25

kde C je čtvercová matice, jejíž prvky jsou obecně funkcí ∆t a ∆x. Například pro centrální schéma (3.22) má C tvar   .. .   σ 1 − σ2   2   σ σ 1 −2   2  σ σ  1 −  2 2  .. .

Je-li tedy dáno řešení U 0 v čase t = 0 (počáteční podmínka), pak pro řešení v čase n∆t platí U n = Cn · U 0 (3.31) Tím dostáváme podmínku stability ve tvaru omezenosti posloupnosti norem kCkn < K

∀n

(3.32)

kde K je jistá konstanta. Konvergence Výše uvedené pojmy konzistence a stability udávaly vlastnosti schémat zvlášť vzhledem k zmenšování kroků ∆x, ∆t a zvlášť pro rostoucí n při již dané diskretizaci. O konvergenci numerického řešení ke skutečnému pak hovoříme v případě hodnoty v pevně daném bodě (xi , tn ) ≡ (i∆x, n∆t) a je tedy třeba při zjemňování sítě (∆x, ∆t → 0) zároveň zvyšovat i a n tak, aby i∆x a n∆t zůstávaly konstantní. Označíme-li chybu vypočteného řešení v uvažovaném bodě εni = uni − u¯(i∆x, n∆t) pak pro konvergenci dostáváme podmínku lim

∆x,∆t→0

|εni | = 0 při konstantních hodnotách xi = i∆x a tn = n∆t

(3.33)

odpovídající požadavku, že při zjemňování sítě konverguje vypočtené řešení ke skutečnému. Jak lze očekávat, konzistence, stabilita a konvergence spolu souvisí a přesné vyjádření tohoto vztahu dává Laxova věta: Pro korektně formulovanou úlohu a konzistentní diskretizační schéma je jeho stabilita nutnou a postačující podmínkou pro konvergenci Pojem korektně formulované úlohy se týká původní diferenciální rovnice a jejích okrajových podmínek, tj. spojitá závislost řešení na koeficientech. Uvedená věta umožňuje na základě vyšetření stability schématu získat úplnou informaci o případné konvergenci ke správnému řešení a tedy o použitelnosti schématu pro výpočet.


3.2.4

26

Von Neumannova metoda vyšetřování stability

Vyšetřování stability je obecně dosti složitá úloha, většina existujících metod se omezuje na lineární problémy, případně i s konstantními koeficienty. Dobře jsou vyřešeny případy, kdy je možno vyloučit nebo zanedbat vliv hranice. To nastává buď u nekonečné oblasti nebo pro periodické okrajové podmínky u oblasti konečné. Druhým případem se zabývá tzv. Von-Neumannova (též Fourierova) metoda, která vychází z předpokladu, že řešený interval délky L na ose x lze periodicky opakovat a pracuje s rozkladem všech veličin (vlastního řešení a chyby) do konečné Fourierovy řady na intervalu délky 2L. Nejprve vyjádříme, jakým způsobem se chyby při výpočtu šíří. Označíme u¯ni přesné řešení diferenčních rovnic a uni vypočtené řešení, které se může vlivem chyb v počátečních datech (mezi ně je možno zahrnout i chyby vzniklé v předchozích krocích výpočtu) lišit. Platí tedy uni = u¯ni + εni (3.34) kde εni značí chybu v příslušném bodě sítě. Při určitém kroku výpočtu pak rovnici schématu musí vyhovovat přesné řešení i řešení zatížené chybou. Odečtením těchto rovnic pak dostaneme vztah pro chyby v odpovídajících bodech sítě, který má stejný tvar jako původní vztah pro vlastní řešení. Například při použití schématu (3.22) pro rovnici konvekce platí pro chybu εn+1 − εni = −σ(εni+1 − εni−1 ) i

(3.35)

Nyní provedeme rozklad chyby do Fourierovy řady, který je možno interpretovat jako rozklad vektoru chyby v n-tém časovém kroku E n = [εn−N , . . . , εn0 , . . . , εnN ]T L do 2N + 1 prvkové báze. N = ∆x označuje počet diskretizačních kroků na intervalu h0, Li, místo něhož nyní uvažujeme dvojnásobný interval hL, Li. Jako bazické vektory použijeme harmonické funkce s jednotkovou amplitudou na hL, Li. Při uvážení diskretizace mohou jejich frekvence, vyjádřené vlnovým číslem kj = 2π/λj , nabývat hodnot násobků nejmenší frekvence, která odpovídá vlnové délce λmax = 2L, až do maximální N násobné, která odpovídá nejmenší vlnové délce na příslušné síti λmin = 2∆x. Konstantní funkce pak odpovídá vlnovému číslu k0 = 0 a celkově tedy máme N + 1 funkcí odpovídajících hodnotám

kj = j kmin = j

2π π =j 2L N∆x

j = 0, 1, . . . , N

(3.36)

Každou nenulovou frekvenci je pak možno vyjádřit dvojicí exponenciálních funkcí eIkj ·i∆x a e−Ikj ·i∆x , kde I je komplexní jednotka, čímž dostáváme potřebných 2N + 1 bazických


27

vektorů a chybu je možno psát εni =

N X

Ejn eIkj ·i∆x =

j=−N

N X

jπ

Ejn eIi N

(3.37)

j=−N

kde Ejn je amplituda j-té složky v n-tém časovém kroku. Označíme-li ještě fázový úhel φ ≡ kj · ∆x =

jπ N

lze j-tou složku chyby vyjádřit v jednoduchém tvaru Ejn eIiφ . Hodnota φ v blízkosti 0 pak odpovídá nízkým frekvencím a v blízkosti π pak vysokým frekvencím. Analýzu stability pak provedeme dosazením složek do vztahu pro časový vývoj chyby (3.35), přičemž hledáme takové podmínky, aby pro žádné j (tj. ani pro žádné φ) nebyla amplituda Ejn s rostoucím n zesilována.

3.3 3.3.1

Vlastnosti schémat MKD pro konvekčně-difuzní rovnici Vlastnosti jednoduchých explicitních schémat

Uvádíme přehled podmínek stability, především v souvislosti s fyzikální podstatou uvažovaných dějů. Odvození některých vztahů je provedeno v příloze 3.6 a je vhodným námětem k procvičení. Upwind schéma pro rovnici konvekce „Upwindÿ způsob aproximace je fyzikálně přirozený pro rovnici konvekce a obecně pro hyperbolické diferenciální rovnice. Fyzikální proces konvekce (ve smyslu transportu obecné veličiny) je vždy spojen s pohybem v určitém směru a to s danou orientací. Upwind numerické schéma respektuje fyzikální orientaci procesu – prostorová diskretizace je volena nesymetricky tak, že se pro výpočet bere v úvahu hodnota proti směru pohybu. Konstrukce schématu je tak obecně závislá na řešené úloze, konkrétně koeficientu u prostorové derivace, vyjadřujícího rychlost posunu. V 1D případě jde o znaménko koeficientu a v rovnici (3.2). Popíšeme situaci pro a > 0, pro opačné znaménko a použitím symetricky otočeného výrazu vzhledem k indexu i (tj. např. dopřednou diferenci níže) bude analogická. Dosazením zpětné diference za ux dostaneme un+1 = uni − a i

∆t n (ui − uni−1 ) ∆x

(3.38)


∆x - e n+1 e e e n i−1 i a∆t

28

dx =a dt

e

i+1

Obrázek 3.2: Vyjádření CFL podmínky pro upwind schéma. kde, stejně jako výše, označíme

∆t (3.39) ∆x což je tzv. Courantovo číslo (někdy též značené Cr ). Podmínka stability pro upwind schéma (3.38) je ∆x 0<σ≤1 tj. ∆t ≤ (3.40) a a nazývá se CFL podmínka (Courant, Friedrichs, Lewy). Pro danou rychlost a prostorovou diskretizaci jde o omezující podmínku na časový krok. σ=a

Fyzikální význam CFL podmínky V daném časovém kroku je ve schématu (3.86) hodnota un+1 vypočtena z hodnot uni a i uni−1 , bere se tedy v úvahu hodnota v tomtéž bodě na ose x a v sousedním bodě proti směru proudění (tj. princip upwind zmíněný v úvodu). Porovnáme nyní postup numerického výpočtu s tím, jakým způsobem se chová řešení diferenciální rovnice. Hodnota řešení rovnice konvekce (3.2) v bodě (x, t + ∆t) je podle charakteristiky x − at = konst rovna hodnotě v bodě (x − a∆t, t). Při výpočtu pomocí upwind schématu je un+1 počítáno jako vážená kombinace uni−1 a uni . Numerický i Výpočet tedy respektuje fyzikální závislost pouze v případě, že bod xi − a∆t leží mezi body xi a xi−1 . Toto nastane právě v případě splnění CFL podmínky: 0 < a∆t < ∆x.


29

Obdobnou úvahu je možno provést úpravou schématu do tvaru un+1 = (1 − σ) uni + σ uni−1 i

(3.41)

kde vidíme, že hodnota uni je konvexní kombinací výchozích hodnot (v případě splnění CFL podmínky). V opačném případě je jeden z koeficientů záporný, což zjevně odporuje fyzikálnímu charakteru úlohy. Numerická difuze Typickým jevem u upwind schématu a obecně u schémat 1.řádu v prostoru je vznik tzv. numerické difuze. Název je odvozen od skutečnosti, že vypočtené numerické řešení se chová jako řešení rovnice, kde by byl navíc přítomen i difuzní člen (případně difuzní člen s větším koeficientem, pokud by šlo o obecnou konvekčně-difuzní rovnici). Příčinou je aproximace 1.řádu: jestliže chyba aproximace je řádu O(uxx ), má tato chyba právě charakter difuzního členu (druhé derivace). Pomocní Taylorova rozvoje lze velikost tohoto členu odhadnout 3 ∂2u ∂ u ui − ui−1 1 aux = a + a∆x 2 (ξxi + (1 − ξ)xi−1 ) + O (3.42) ∆x 2 ∂x ∂x3 kde parametr ξ ∈ h0, 1i udává polohu hodnoty druhé derivace v intervalu mezi dvěma uzly diskretizace. Hodnotu numerické difuze je možno udat přesněji, a to v závislosti na Courantově čísle 1 Dnum = a∆x(1 − σ) (3.43) 2 a pohybuje se tedy od nulové hodnoty pro časový krok přesně na mezi stability do hodnoty 21 a∆x pro časový krok blížící se k nule. Vztah pro obecné časově i prostorově vážené schema (3.55) lze najít v [11]. Centrální schéma pro rovnici difuze Z rovnice ut = b uxx dostáváme použitím centrální diference (3.18) schéma un+1 = uni + β(uni+1 − 2uni + uni−1 ) i kde β=b

(3.44)

∆t (∆x)2

Podmínka stability je β≤

1 2

(3.45)


30

a její smysl lze intuitivně vidět obdobně jako u CFL podmínky: ve schématu musí být koeficient (1 − 2β) u členu uni (původní hodnoty v uvažovaném bodě) nezáporný. Vyjádřená podmínka pro časový krok je ∆t ≤

(∆x)2 2b

(3.46)

což může znamenat značné omezení při zjemňování sítě, vzhledem ke kvadratické závislosti na ∆x. Upwind schéma pro konvekčně-difuzní rovnici Označení „upwindÿ se týká pouze konvekčního členu, u difuze nemá význam (jde o fyzikálně symetrický děj). Schéma má po úpravě a zavedení výše uvedených bezrozměrných čísel tvar un+1 = uni − σ(uni − uni−1 ) + β(uni+1 − 2uni + uni−1) (3.47) i Podmínka stability je resp. pro časový krok

σ + 2β ≤ 1 1 ∆t = 2b + (∆x)2

a ∆x

(∆x)2 = 2b + a∆x

(3.48) (3.49)

Je vidět, že pro a → 0 nebo b → 0 přechází podmínka ve vztahy získané pro příslušné speciální případy – rovnice konvekce a rovnice difuze. Centrální schéma pro konvekčně-difuzní rovnici Tento případ je zajímavý v tom, že konvekce je diskretizována způsobem, který u rovnice bez difuzního členu vede na nestabilní schéma. Přítomnost difuze tedy schéma jistým způsobem „stabilizujeÿ, lze však očekávat obtíže v případě, kdy koeficient difuze bude blízký nule (a tedy fyzikální proces bude blízký prosté konvekci). Schéma má tvar σ un+1 − uni + (uni+1 − uni−1) = β(uni+1 − 2uni + uni−1 ) (3.50) i 2 Podmínka stability je σ 2 ≤ 2β ≤ 1

(3.51)

zapíšeme-li souhrnně obě podmínky, z nichž vždy jedna je silnější než druhá, v závislosti <

na platnosti 2β ≥ σ (viz odvození v příloze). Pro časový krok má podmínka tvar (∆x)2 2b , (3.52) ∆t ≤ min 2b a2


31

σ 6

σ=1

C

C

C

C

C

CFL

-

Pe = 2

Pe

Obrázek 3.3: Grafické znázornění podmínek stability pomocí Pecletova a Courantova čísla a srovnání s CFL podmínkou. Viz též obr.3.8, vyjádření pomocí diskretizačních kroků ∆x a ∆t. Zavedeme tzv. Pecletovo číslo sítě (grid Peclet number) vztahem Pe =

a∆x σ = β b

(3.53)

tedy jako poměr míry konvekce a difuze,1 zároveň v souvislosti s jemností sítě. Pomocí Pecletova čísla můžeme přehledně vyjádřit dělící hranici mezi situacemi, kdy je silnější jedna nebo druhá podmínka stability. Pro hodnoty Pe < 2 je výraznější vliv difuze a uplatňuje se druhá nerovnost v (3.51), což je podmínka, která by musela být splněna při prosté difuzi. Pro hodnoty Pe > 2, tj. silnější konvekci, se uplatňuje první nerovnost v (3.51), kterou můžeme chápat např. jako omezení na Courantovo číslo σ. Z uvedené nerovností je zřejmá v úvodu zmíněná problematická situace pro úlohy s dominantní konvekcí, tj. s velkým Pecletovým číslem – podmínka stability může být velmi restriktivní. Navíc v tomto případě stabilita ještě nezaručuje uspokojivé numerické výsledky. Ukazuje se, že právě pro Pe > 2 (viz výše zmíněnou dělící podmínku) vykazuje řešení oscilace kolem míst s velkým gradientem hodnot řešení u. Tento jev je demonstrován na příkladě v kapitole 3.7.


     

·

· · −β (2β − 1) −β −β (2β − 1) −β −β (2β − 1) ·  · ·  −σ 1 + σ   −σ 1 + σ   −σ 1 + σ · ·

32

 −β · · 

    

    

Obrázek 3.4: Tvar matice implicitního schématu při použití centrální diference 2.řádu pro difuzi a při použití upwind diference pro konvekci.

3.3.2

Implicitní schémata

Zmíníme dva charakteristické případy: upwind schéma pro konvekci a centrální schéma pro difuzi. Obdobně jako u explicitního schématu jsou pak kombinovány při diskretizaci konvekčně-difuzní rovnice. V obou případech jde o schéma nepodmíněně stabilní, tedy připouštějící použití libovolně velkého časového kroku. Úspora výpočetního času díky eventuálnímu provedení nižšího počtu časových kroků v daném intervalu je však vyvážena nutností řešení soustavy lineárních rovnic v každém kroku. Významný rozdíl co se týče charakteru matice je mezi konvekcí a difuzí. Výpočet samotné difuze vede na řešení soustavy se symetrickou maticí, zatímco konvekce vede na nesymetrickou matici. To je zřejmé již z přítomnosti hodnot v sousedních uzlech sítě – centrální schéma je v tomto smyslu symetrické a upwind schéma (jednostranná diference) nesymetrické. V případě rozsáhlých úloh a velké řídké matice existují pro případ symetrických matic relativně spolehlivé a rychlé algoritmy (např. Krylovovské metody), zatímco řešení nesymetrické matice může způsobovat problémy. Symetrie či nesymetrie algebraické úlohy je sice bezprostředně dána použitým schématem, musíme si však uvědomit, že souvisí již s původním fyzikálním procesem – difuze je symetrického charakteru a konvekce nesymetrického (je dána orientací proudění). Jak jsme již naznačili, zmíněné vlastnosti schémat a vzniklé jevy jsou obecné a obdobným způsobem se projevují i u jiných numerických metod, tedy konečných objemů 1

Připomeňte si jiné definice Pecletova čísla v podobném smyslu (poměr konvekce a difuze) za jiných okolností, např. Pe= vd D pro hydrodynamickou disperzi v porézním prostředí. Uvědomme si rovněž, že vzhledem k požadované bezrozměrnosti takovéhoto parametru, bude vždy k vyjádření vztahu mezi konvekcí a difuzí nezbytné zavedení jisté charakteristické délky – pro numerickou metodu je to parametr sítě.


33

a konečných prvků.

3.3.3

Obecné parametrické schéma

Použití různých typů schémat, jak byly uvedeny výše, lze obecně kombinovat. Motivací jsou mimo jiné rozporné vlastnosti zmíněných základních přístupů, např. numerická difuze u upwind schématu (a obecně u schémat 1.řádu v prostoru) a oscilace u centrálního schématu. Obecné schéma zkonstruujeme tak, že s určitým váhovým koeficientem α nakombinujeme dopřednou a zpětnou diferenci v prostoru a s váhovým koeficientem ω hodnoty v časovém kroku n a n+1. Podrobně je tato úvaha provedena v knize [11]. Za prostorovou derivaci v konvekčním členu tedy dosadíme (ux )i ≈ α

ui − ui−1 ui+1 − ui + (1 − α) ∆x ∆x

(3.54)

Použitím α = 21 jde o centrální schéma, pro α = 0 nebo α = 1 pak v závislosti na znaménku a dostaneme upwind schéma. V případě časové derivace nepůjde o vážení dopředné a zpětné diference, ale o vážení hodnot v diskretizovaném konvekčním a difuzním členu v krocích n a n + 1, které se vyskytují v použité dopředné časové diferenci (pro zpětnou diferenci by byla situace obdobná, šlo by jen o posun o 1 v časovém indexu). Schematicky můžeme zapsat un+1 − uni i = ω(buxx − aux )n+1 + (1 − ω)(buxx − aux )ni i ∆t

(3.55)

Dosazením prostorově vážené diference (3.54) za konvekční člen a centrální diference za difuzní člen můžeme obecné schéma s použitím bezrozměrných čísel Courantova σ a Pecletova Pe (viz (3.39) a (3.53)) zapsat následovně nσ n un+1 − u = (1 − ω) (un − 2uni + uni−1 ) i i Pe i+1 o n n n n − σ [(1 − α)ui + αui+1 − (1 − α)ui−1 − αui ] nσ ω (un+1 − 2un+1 + un+1 i i−1 ) Pe i+1 o n+1 n+1 n+1 n+1 − σ [(1 − α)ui + αui+1 − (1 − α)ui−1 − αui ] (3.56)

Ihned je zřejmé, že pro ω = 0 jde o explicitní schéma a pro ω = 1 o schéma implicitní. Použitím „symetrickéhoÿ váhového parametru ω = 21 dostáváme schéma, označované jako Crank-Nicholsonovo. Toto schéma je druhého řádu přesnosti v čase a nepodmíněně stabilní (tj. pro libovolné kombinace koeficientů a diskretizace). I když jsme jej (vzhledem k výše uvedené formulaci pro hodnoty ω) terminologicky odlišili od implicitního


34

schématu, jde pochopitelně opět o schéma implicitní v tom smyslu, že vyžaduje řešení soustavy lineárních rovnic v každém časovém kroku. Pro úplnost uvedeme vyjádření numerické difuze pro obecné schéma v závislosti na parametrech: 1 1 Dnum = a∆x −α +σ ω− (3.57) 2 2 (všimněte si, že numerická difuze je nulová např. pro Crank-Nicholsonovo schéma nebo pro explicitní upwind při σ = 1 (mez stability)).

3.3.4

Parametrizace pomocí „uměléÿ difuze

Ukážeme některé další zajímavé souvislosti mezi schématy líšícími se různým poměrem mezi upwind a centrální aproximací konvekce. Pro zjednodušení provedeme celé odvození jen pro rovnici konvekce, ale vzhledem k interpretaci výsledků by mělo být zřejmé, že úvahy jsou platné i pro konvekčně-difuzní rovnici. Ze vztahu (3.57) je vidět, že numerická difuze poměrně jednoduchým způsobem souvisí s parametrem α (váhový parametr v aproximaci konvekce). Tuto souvislost vyjádříme přesněji. Uvažujeme obecnou tříbodovou aproximaci konvekčního členu a explicitní schéma v čase un+1 = a−1 uni−1 + a0 uni + a1 uni+1 (3.58) i kde ai jsou koeficienty. Ty pochopitelně budou záviset na koeficientech rovnice a diskretizaci, ale na základě předchozího textu můžeme očekávat, že v jejich volbě bude jistá volnost (předtím jsme v aproximaci konvekce měli volný parametr α). Aby bylo tedy schéma konzistentní s rovnicí, musí platit: 1.

a−1 + a0 + a1 = 1 odpovídá „konzervativnostiÿ schématu (nultý řád přesnosti) – z konstatního rozložení neznámé funkce musíme v dalším kroku dostat tutéž neznámou hodnotu

2.

a−1 − a1 = σ odpovídá aproximaci první derivace (1. řád přesnosti)

Pokud tedy uvažujeme schéma 1. řádu přesnosti, máme předepsány dvě rovnice pro tři neznámé koeficienty a tedy jeden „stupeň volnostiÿ ve volbě parametrů. To je názorně vidět při vyjádření řešení soustavy rovnic jako součtu obecného řešení homogenní soustavy rovnic (s nulovou pravou stranou) a partikulárního řešení původní soustavy. Jednoduchým výpočtem najdeme homogenní řešení γ · (1, −2, 1) a partikulární např. σ σ , 1, − 2 , tj. 2 σ σ , 1, − + γ · (1, −2, 1) (a−1 , a0 , a1 ) = 2 2


35

Parametr γ je zmíněný volný parametr. Vidíme, že volba γ = 0 odpovídá centrální aproximaci a volba γ = σ2 upwind aproximaci. Podstatné na tomto vyjádření je to, že homogenní řešení vlastně odpovídá diferenčnímu vzorci pro druhou derivaci, tedy diskretizaci nějakého difuzního členu. Parametr γ se nazývá umělá difuze – určuje jak moc difuze se objeví ve výpočtu vlivem aproximační chyby v konvekčním členu, tedy výše komentované numerické difuze. Hodnotu gamma však nelze přímo zaměňovat s hodnotou numerické difuze v podobě, jaké byla uváděna v předchozím textu. Význam γ totiž je závislý na zvoleném partikulárním řešením soustavy a numerická difuze pro γ = 0 nemusí být nulová. Je však možno říci, že všechna schémata 1. řádu [rp konvekci se liší právě a jen množstvím umělé/numerické difuze. Hodnotu γ, pro kterou bude numerická difuze nulová (nezávisle na vztahu ∆x, ∆t a a – srovnejte případ upwind a σ = 1) určíme z podmínky konzistance pro schéma 2. řádu a ta je a−1 + a1 = σ 2 2

Řešení soustavy nyní již třech rovnic odpovídá volbě parametru γ = σ2 (tato hodnota umělé difuze odpovídá nulové numerické difuzi). Příslušné schéma nese jméno LaxWendroffovo σ σ un+1 = uni−1 (1 + σ) + uni (1 − σ 2 ) − uni+1 (1 − σ) i 2 2 Schéma je podmíněně stabilní (podobně jako explicitní upwind) a produkuje oscilující řešení (jako každé schéma s některým koeficientem a−1 , a1 záporným. Nyní nám ještě zbývá dát do souvislosti parametrizaci pomocí α a pomocí γ. rozklad matice na symetrickou a antisymetrickou část 1 1 A= A + AT + A − AT 2 2 (α, 1 − 2α, α − 1)

symetrická:

1 1 1 [(α, 1 − 2α, α − 1) + (α − 1, 1 − 2α, α)] = α − , 1 − 2α, α − = 2 2 2 1 ·(1, −2, 1) = α− 2 | {z } −γ, σ antisymetrická: 1 1 1 ui+1 − ui−1 [(α, 1 − 2α, α − 1) − (α − 1, 1 − 2α, α)] = , 0, − ∼ 2 2 2 2∆x σ2 1−σ 1 −α = ⇒ α= γ=σ 2 2 2


3.4

36

Formulace metody konečných prvků pro konvekčně-difuzní rovnici

Předpokádáme znalost základních myšlenek a základů matematické teorie metody konečných prvků (MKP), pro stacionární i nestacionární úlohy. Ukážeme na příkladu formulaci MKP pro konvekčně-difuzní rovnici, tedy především diskretizaci konvekčního členu, který je rozdílem oproti modelovým eliptickým a parabolickým úlohám používaným k výkladu MKP. Dále ukážeme principy tzv. Petrov-Galerkinovy metody, která je zobecněním Galerikovy metody a umožní zavést do konečných prvků upwind techniku, která je jak jsme v předchozích sekcích poznali užitečná pro řešení konvekčně-difuzních úloh.

3.4.1

Variační formulace a metoda konečných prvků

Pro úspornější a přehlednější vyjádření uvažujeme rovnici bez zdrojů a s nulovými okrajovými podmínkami. ∂u − ∇ · (D∇u) + ∇ · (vu) = 0 ∂t

(3.59)

Slabou formulaci této rovnice dostaneme standardní operací násobením testovací funkcí φ a integrací přes oblast řešení úlohy Ω a difuzní člen dále upravíme pomocí Greenovy věty, tj. musí platit pro všechna φ z prostoru H1 (Ω): ∂u ,φ + (D∇u, ∇φ)Ω + (∇(vu), φ)Ω = 0 (3.60) ∂t Ω Princip Galerkikovy metody aproximace této slabé fomulace je, že nepožadujeme splnění pro všechna φ z nekonečnědimenzionálního Sobolevova prostoru H1 , ale pouze pro funkce z vhodně zvoleného konečnědimenzionálního prostoru a zároveň i neznámé funkce hledáme v „stejnémÿ konečnědimenzionálním prostoru. Metodu konečných prvků pak můžeme chápat jako konkrétní postup, jak zvolit bázové funkce pro Galerkinovu metodu – tj. funkce s malým supportem s definovaným stupněm hladkosti. Označíme N dimenzi aproximačního prostoru a φi , i = 1, . . . , N bázové funkce. Potom neznámou funkci a časovou derivaci můžeme vyjádřit u=

N X i=1

N

Ui φi

∂u X ∂Ui = φi ∂t ∂t i=1

kde Ui , i = 1, . . . , N jsou hledané koeficienty – neznámé.

(3.61)

Kapitola 3. Numerické metody pro konvekčně-difuzní rovnici fi

f i-1

x i-1

37

f i+1 x i+2

x i+1

xi

Obrázek 3.5: dfi/dx 1

1/Dxi

fi

x i-1

Dxi

x i-1

xi

x i+1

xi

x i+1

-1/Dxi+1

Obrázek 3.6: Diskrétní (konečněrozměrná) formulace původní úlohy dostaneme dosazením lineárních kombinací do rovnice (3.60) N X ∂Ui i=1

∂t

(φi , φj )Ω +

N X i=1

Ui (D∇φi , ∇φj )Ω +

N X i=1

Ui (∇ · (vφ), φj )Ω

j = 1...N

(3.62)

tj. dostaneme soustavu N rovnic pro N neznámých koeficientů Ui (t). V této podobě jde o soustavu obyčejných diferenciálních rovnic, kterou můžeme řešit libovolnou metodou, v Uin −Uin−1 i nejjednodušším případě dosazením zpětných diferencí za časové derivace ∂U ≈ ∂t ∆t dostaneme soustavu N algebraických rovnic pro hodnoty koeficientů v n-tém časovém kroku (M + ∆tK)U n = MU n−1 (3.63) kde U n = (U1n . . . UNn )T a Mij = (φi , φj )Ω

3.4.2

Kij = (D∇φi, ∇φj )Ω + (∇ · (vφ), φ)Ω

(3.64)

Diskrétní formulace v 1D s lineárními bázovými funkcemi

Na jednoduchém příkladu ukážeme souvislost metody konečných prvků s jednoduchými lineárními bazickými funkcemi a metody konečných diferencí s centrální aproximací konvekčního členu. Uvažujeme 1D úlohu a oblast řešení jako úsečku (interval) (x0 , xN ) rozdělenou na N intervalů ∆xi = xi − xi−1 , i = 1 . . . N. Definujeme bázové funkce jako po částech lineární funkce, které v jednom uzlu diskretizace nabývají hodnoty 1 a ve všech zbylých uzlech hodnoty 0 (obrázek 3.5). Koeficienty lineární kombinace (3.61) pak přirozeně mají význam hodnot aproximované (hledané) funkce v bodech diskretizace xi .


38

Vypočteme nyní složky matic M (matice hmotnosti) a K (matice tuhosti). Na to i potřenujeme vyjádřit součiny funkcí φi a jejich derivací dφ připadající různým dvojicím dx uzlů diskretizace. Průběhy funkcí pro zvolený bod i jsou zakresleny na obrázku 3.6. Podstaným faktem na kterém je MKP založena je, že díky malému supportu funkcí φi bude většina skalárních součinů (prvků matic M a K) nulová. Integrály přes celou oblast Ω (v našem případě interval (x0 , xN )) se pak redukují na integrály přes jeden (v případě prvků matic těsně vedle diagonály – dva sousední body diskretizace) nebo dva (v případě prvků matic na diagonále) konečné prvky – intervaly diskretizace. Vypočteme prvky matic pro případ ekvidistantních intervalů ∆xi = 1 a po elementech (intervalech) konstantní koeficienty D a v (můžeme vytknout před integrál). Integrály výsledných konstantních, lineárních nebo kvadratických funkcí určíme zpaměti, triviálním výpočtem nebo geometrickou analogií. Platí 2 1 1 (φi, φi ) = (φi+1 , φi ) = časový/hmotnostní(3.65) člen (φi−1 , φi) = 6 3 6 (∇φi−1 , ∇φi) = −1 (∇φi , ∇φi ) = 2 (∇φi+1 , ∇φi ) = −1 difuzní člen (3.66) 1 1 (∇φi−1 , φi ) = − (∇φi , φi) = 0 (∇φi+1 , φi ) = konvekční člen (3.67) 2 2 Z hodnot prvků matice pro konvekční a difuzní člen ihned vidíme ekvivalenci s metodou konečných diferencí (centrálním schématem); i-tá rovnice pro neznámé Uin je n−1 n−1 n n n n Ui−1 + 4Uin−1 + Ui+1 Ui+1 − Ui−1 Ui−1 + 4Uin + Ui+1 n n n − D(Ui−1 − 2Ui + Ui+1 ) + v = 6∆t 2 6∆t (3.68) Rozdíl je v časovém členu (matici hmotnosti), kde v případě MKD je v i-té rovnici pouze hodnota v i-tém uzlu (nikoli ve dvou sousedních jako v MKP) – matice by byla diagonální. Při řešení nestacionárních úloh pomocí MKP však lze použít tzv. techniku „mass lumpingÿ (soustředění hmoty), kdy matici hmotnosti nahradíme diagonální maticí, jejíž prvky jsou součty řádků původní matice (součet „hmotyÿ). Tím je levá strana rovnice plně ekvivalentní pro MKP i MKD. Podstatný důsledek pro řešení konvekčně-difuzních úloh však je, že použití standardní MKP (obecnějších bázových funkcí a obecnější geometrie) povede na obdobné numerické obtíže jako použití centrálního diferenčního schématu pro konvekční člen, tj. vznik oscilací při velkém Pecletově čísle. Zavedení upwind (s určitou váhou) aproximace konvekčního členu, která tento jev odstraňuje, i v případě MKP, umožňuje tzv. Petrov-Galerkinova metoda.

3.4.3

Petrov-Galerikova metoda – upwinding

Nejprve dáme do souvislosti termíny Galerkinova metoda a Metoda konečných prvků. Oba se vztahují k definici konečněrozměrné aproximace slabé formulace. Tato aproximace je založena na dvou bodech


39

• neznámou u hledáme v konečněrozměrném prostoru (aproximačním) – to je princip MKP • rovnice platí pro testovací funkce z konečněrozměrného prostoru (testovacího) – to je princip Galerkinovy metody Petrov-Galerikova metoda pak znamená změnu v druhém bodě – jako testovací funkce použijeme jiné funkce než pro aproximaci neznámé, ty označíme ψi (jejich počet však musí být stejný, abychom dostali soustavu rovnic se stejným počtem rovnic a neznámých). Konečněrozměrná Petrov-Galerkinova aproximace tedy je N X ∂Ui i=1

∂t

(φi , ψj )Ω +

N X i=1

Ui (D∇φi, ∇ψj )Ω +

N X

Ui (∇·(vφ), ψj )Ω = 0 j = 1 . . . N (3.69)

i=1

Konstrukci aproximačních a testovacích funkcí ukážeme opět na příkladu v 1D. Ty budou zobecněním po částech lineárních funkcí φj použitých v sekci xx jako příklad bazických funkcí pro 1D MKP. Nyní v Petrov-Galerkinově MKP uvažujeme jako aproximační bazické funkce tytéž φj a jako testovací funkce ψj = φj + ασj (x) kde funkce σj jsou definovány jako  j −x)  3(x−xj−1 )(x x ∈ ej ∆x2j σj = −3(x−xj )(xj+1 −x)  x ∈ ej+1 ∆x2

(3.70)

(3.71)

j+1

a vyjadřují upwind efekt. Koeficient α je váhový parametr, který má význam umělé difuze přidané do schématu, což lze ekvivalentně vyjádřit (viz sekce xxx) jako váha odpovídající parametrizaci diferenčního schématu mezi centrální a plně upwind aproximací. Dosazením bázových funkcí do konečněrozměrné aproximace slabé formulace dostaneme rovnici obdobnou jako (3.68), kde konvekční člen bude v vα (Uj+1 − Uj−1 ) + (Uj+1 − 2Uj + Uj−1 ) = 2 2 tj. přibude umělá difuze s koeficientem

3.5

vα , 2

(3.72)

a pro α = 1 dostáváme upwind schéma.

Metoda rozkladu operátoru

Jde o speciální techniku časové diskretizace, kterou nelze zaředit do nějakého hierarchického třídění numerických metod, naopak – umožňuje např. kombinovat různé numerické


40

metody pro různé fyzikální jevy zahrnuté v jedné úloze. Nejprve ukážeme obecné algebraické vyjádření a potom konkrétní aplikaci na konvekčně-difuzní rovnici. To je typický příklad aplikace techniky rozkladu operátoru – „operátor konvekce-difuzeÿ se rozloží na „operátor konvekceÿ a „operátor difuzeÿ a výpočet časového vývoje probíhá ve střídání kroků konvekce a difuze. Použití různých numerických metod pro konvekci a pro difuzi je přirozený postup, vzhledem k jejich odlišnému fyzikálnímu charakteru, čímž se zmírní výskyt nežádoucích numerických jevů – oscilací a numerické difuze (prezentováno v předchozích kapitolách). Uvažujeme soustavu obyčejných diferenciálních rovnic (takovýto systém typicky reprezentuje časový vývoj nějakého diskrétního systému, včetně speciálního případu, kdy tento diskrétní systém vznikl diskretizací v prostoru spojitého systému např. MKD nebo MKP) Ut = AU, kde U je vektor funkcí ui(t) délky I, A je matice I × I a index t označuje časovou derivaci po složkách. Diskretizací v čase (např. explicitní Euler) dostaneme   un1 Un+1 − Un   = A · Un + O(∆t2 ) Un =  ...  ∆t unI tj.

Un+1 = (I + ∆tA)Un + O(∆t2 )

explicitní

Jde o metodu 1. řádu přesnosti. Provedeme nyní zmíněný rozklad operátoru – v tomto případě to znamená vyjádření matice (operátoru) A jako součet (rozklad) A = A1 + A2 . Schéma pak upravíme následovně: Un+1 = (I + ∆tA1 + ∆tA2 )Un + 0(∆t2 ) Un+1 = ([I + ∆tA1 ] · [I + ∆tA2 ] − ∆t2 A1 A2 )Un + 0(∆t2 ) | {z } | {z } skládání operátorů stejný řád chyby

(3.73)

Součin [I + ∆tA1 ] · [I + ∆tA2 ] („skládáníÿ operátorů) odpovídá dvěma po sobě jdoucím časovým krokům Eulerovy metody, ale s různými maticemi – tedy jeden krok „řízenýÿ rovnicí Ut = A1 U a druhý krok „řízenýÿ rovnicí Ut = A2 U. Reálně však jde o jeden a tentýž časový úsek a aplikaci dvou částí „dynamických vlastnostíÿ systému každé zvlášť. Podstatné je, že schémata s rozkladem a bez rozkladu jsou v jistém smyslu ekvivalentní – mají stejný řád chyby – nemusí však mít stejné numerické vlastnosti (stabilita). Vztah je obdobný jako např. mezi implicitní a explicitní Eulerovou metodou: (I − ∆tA)Un = Un−1 + 0(∆t2 ) (I − ∆tA)−1 = I + ∆tA + ∆t2 A2 + · · · Un = (I + ∆tA)Un−1 + ∆t2 A2 + · · · + 0(∆t2 ) | {z } 0(∆t2 )

implicitní (3.74) pokud ∆tkAk < 1(3.75) explicitní (3.76)


41

Popsanýn rozkladem matice a skládáním operací lze vyjádřit mnoho dalších speciálních metod – metody typu „prediktor–korektorÿ, „metoda střídavých směrůÿ (matice A vzniklá MKD diskretizací se rozloží na prvky odpovídající diskretizaci ve směru x a prvky odpovídající diskretizaci ve směru y. Viz podrobně [7, 9, 10]. Ukážeme nyní konkrétnějčí případ rovnice Ut = AU a rozkladu diskretizaci 1D konvekčně-difuzní rovnice, tj.  · · ·  −β − σ (2β + σ) −β   −β − σ (2β + σ) −β A=  −β − σ (2β + σ) −β · · ·

matice A a to      

(3.77)

kde jsme uvažovali centrální schéma pro difuzi a upwind schéma pro konvekci, β = ∆x ∆t b (∆x) , a, b jsou koeficienty rovnice (viz sekce 3.3). Rozklad matice na 2 a ∆t ≤ a difuzní a konvekční část je přirozeně   · · ·   −β 2β −β     −β 2β −β A1 =  (3.78)    −β 2β −β · · · 

  A2 =   

·

· −σ

σ −σ

σ −σ σ · ·

     

(3.79)

Vyjádření dvou kroků výpočtu (jeden „reálnýÿ časový krok) je (jiný tvar rovnice (3.73)) U˜n+1 = (I + ∆tA1 )Un + 0(∆t2 ) ˜n+1 + 0(∆t2 ) Un+1 = (I + ∆tA2 )U

(3.80) (3.81)

Tento výpočetní postup je pak možno chápat jako diskrétní vyjádření dvou „samostatnýchÿ fyzikálních procesů – první krok odpovídá řešení rovnice difuze ut = buxx v časovém úseku ∆t a druhý řešení rovnice konvekce ut + aux = 0 v časovém úseku ∆t. To je obecnější pohled na metodu rozkladu operátoru – nějaký složitější modelovaný proces vyjádříme jako rozklad jednodušších procesů a na každý aplikujeme samostatnou výpočetní metodu.


3.6

42

Cvičení: Výpočet podmínek stability Fourierovou metodou

Pro jednotlivá schémata (pro obecnou rovnici (3.1) a speciální případy a = 0, tj. difuze a b = 0, tj. konvekce) zmíněná v kapitole 3.2 vypočteme Von Neumannovou metodou podmínky stability. Všechny výpočty jsou vhodným námětem k samostatnému procvičení.

3.6.1

Centrální schéma pro konvekční rovnici

Po dosazení libovolné složky chyby Ejn eIiφ , při vynechání indexu j, dostáváme a E n+1 − E n Iiφ e + (E n eI(i+1)φ − E n eI(i−1)φ ) = 0 ∆t 2∆x po úpravě

σ n Iφ E (e − e−Iφ ) = 0 (3.82) 2 Podmínku stability (3.29) musí splňovat amplituda každé složky chyby (každé j, resp. φ). Definujeme-li E n+1 (3.83) G= En faktor zesílení, dostáváme podmínku ve tvaru (E n+1 − E n ) +

|G| ≤ 1 pro ∀φ (resp. ∀j)

(3.84)

Dělíme-li (3.82) E n , dostáváme vztah pro G G=1−

σ Iφ (e − eIφ ) = 1 − Iσ sin φ 2

platí tedy |G| =

q

1 + σ 2 sin2 φ ≥ 1

(3.85)

což znamená, že podmínka stability není splněna - schéma je nepodmíněně nestabilní.

3.6.2

Upwind schéma pro konvekční rovnici

Schéma dostaneme použitím zpětné diference za ux un+1 − uni + σ(uni − uni−1 ) = 0 i přičemž význam σ je stejný jako v předchozím. Po dosazení chyby pak (E n+1 − E n )eIiφ + σE n (eIiφ − eI(i−1)φ ) = 0

(3.86)


43

η 1

σ 1−σ 0

G

φ

1 ξ

Obrázek 3.7: Grafické znázornění podmínky stability pro G a pro zesílení tedy platí G = 1 − σ + σe−Iφ = 1 − σ + σ cos φ − Iσ sin φ Podmínku stability |G| ≤ 1 můžeme snadno ověřit geometricky v Gaussově rovině. Zapíšeme-li reálnou a imaginární složku G ξ ≡ ReG = (1 − σ) + σ cos φ η ≡ ImG = −σ sin φ

(3.87)

dostáváme parametrické rovnice kružnice s φ jako parametrem. Pro různá φ tedy faktor G leží na kružnici se středem na reálné ose v (1 − σ) a poloměrem σ. Aby bylo schéma stabilní, musí pro libovolné φ být |G| ≤ 1, uvažovaná kružnice tedy musí ležet uvnitř jednotkového kruhu. Jak je vidět z obrázku, toto nastává pro σ ≤ 1. Z definice je σ < 1 a schéma je tedy stabilní, platí-li pro σ (Courantovo číslo) 0<σ≤1

(3.88)

což je tzv. CFL podmínka (viz hlavní text, sekce 3.3.1). Zvolíme-li si tedy prostorový diskretizační krok ∆x, dostáváme pro časový krok podmínku ∆t ≤

∆x a

(3.89)


3.6.3

44

Centrální schéma pro difuzní rovnici

Z rovnice ut = b uxx dostáváme použitím (3.18) schéma un+1 − uni = β(uni+1 − 2uni + uni−1 ) i

(3.90)

∆t kde β = b (∆x) 2 . Opět dosadíme rozklad chyby

E n+1 − E n = βE n e−Iφ − 2βE n + βE n eIφ a je tedy G = 1 − 2β + 2β cos φ = 1 − 2β(1 − cos φ)

(3.91)

Výraz v závorce nabývá hodnot v intervalu < 0, 2 >. Aby byla splněna podmínka |G| ≤ 1, musí i celý výraz 2β(1 − cos φ), čímž dostáváme podmínku stability ve tvaru

1 2 Je-li nyní dán krok ∆x, musí pro časový krok platit β≤

∆t ≤

3.6.4

(∆x)2 2b

(3.92)

(3.93)

Upwind schéma pro konvekčně-difuzní rovnici

Schéma má tvar (un − uni−1) (un − 2uni + uni−1 ) un+1 − uni i +a i = b i+1 ∆t ∆x (∆x)2

(3.94)

Obdobným postupem jako doposud dostaneme (E n+1 − E n ) + σ(E n − E n e−Iφ ) = βE n eIφ − 2βE n + βE n e−Iφ tj.

G = 1 − σ + σ(cos φ − I sin φ) + 2β cos φ − 2β ξ = Re G = 1 − (σ + 2β) + (σ + 2β) cos φ η = Im G = −σ sin φ

(3.95)

σ + 2β ≤ 1

(3.96)

což jsou parametrické rovnice elipsy se středem v 1 − (σ + 2β) na reálné ose a poloosami (σ + 2β) ve směru reálné osy a σ ve směru imaginární osy. Elipsa se tedy dotýká jednotkové kružnice a je protáhlá ve směru reálné osy, neboť podle předpokladu o koeficientech rovnice β ≥ 0 a σ ≥ 0. Celá tedy bude ležet v jednotkovém kruhu, jestliže


45

což je hledaná podmínka stability. Opět ji vyjádříme jako podmínku pro časový krok ∆t =

1 2b + (∆x)2

a ∆x

=

(∆x)2 2b + a∆x

(3.97)

Je vidět, že pro a → 0 nebo b → 0 přechází podmínka ve vztahy získané pro příslušné speciální případy rovnice (3.1).

3.6.5

Centrální schéma pro konvekčně-difuzní rovnici

Použitím centrální diference (3.17) dostaneme σ un+1 − uni + (uni+1 − uni−1) = β(uni+1 − 2uni + uni−1 ) i 2 Provedeme-li obvyklé úpravy, pak σ G − 1 + (eIφ − e−Iφ ) = β(2 cos φ − 2) 2 a ξ = Re G = 1 − 2β + 2β cos φ η = Im G = σ sin φ

(3.98)

(3.99)

Hodnoty G tedy opět leží na elipse se středem na reálné ose a dotýkající se jednotkové kružnice. Nyní však není pevně určeno, která z poloos je větší, a je třeba vyšetřovat oba možné případy samostatně: 1. Je-li 2β ≥ σ, je poloosa ve směru reálné osy delší (popř. jsou obě stejné) a podobně jako v předcházejících případech dostáváme podmínku 1 2

(3.100)

(∆x)2 2b

(3.101)

β≤ pro časový krok pak ∆t ≤

2. Je-li 2β < σ, má elipsa delší poloosu ve směru imaginární osy. Vyjádření skutečnosti, že leží celá uvnitř jednotkového kruhu, pomocí nerovností je nyní poněkud složitější – geometrickou úvahou lze dospět k podmínce, že poloměr křivosti elipsy v bodě dotyku musí být menší než poloměr jednotkového kruhu, tj. 1. Tím dostáváme σ2 ≤1 tj. σ 2 ≤ 2β (3.102) 2β a pro časový krok 2b (3.103) ∆t ≤ 2 a

Kapitola 3. Numerické metody pro konvekčně-difuzní rovnici ∆t

(3)

46

(1)

(2)

∆x

Obrázek 3.8: Grafické znázornění podmínek stability pro ∆t <

Vyjádříme nyní i dělící podmínku 2β ≥ σ pomocí diskretizačních kroků

2b a Snadno se lze přesvědčit, že nastává-li případ (1), je podmínka (3.103) z případu (2) slabší než (3.101) a naopak nastává-li (2), je podmínka v (1) slabší. Souhrnně lze tedy podmínku stability pro ∆t zapsat vztahem (∆x)2 2b ∆t ≤ min , (3.104) 2b a2 >

∆x ≤

Situace je nejlépe vidět z grafu, na obrázku 3.8 křivky (1), (2) znázorňují podmínky v prvním resp. druhém případě a (3) příslušná dělící hranice. Vyšrafována je množina přípustných hodnot ∆x, ∆t.

3.7

Cvičení: Demonstrační numerické výpočty

Vlastnosti vyšetřovaných schémat předvedeme na modelových úlohách s konkrétními hodnotami parametrů. Za oblast řešení zvolíme např. interval (0, 1) a rovnici doplníme příslušnými počátečními a okrajovými podmínkami: u(x, 0) = u0 (x) u(0, t) = u1 (t) ∂u = 0 ∂x x=1

(3.105) (3.106) (3.107)


47

1 pocatecni podminka N=20 k=1 N=20 k=0,5 N=200 k=1 N=200 k=0,5 0.8

0.6

0.4

0.2

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Obrázek 3.9: Řešení rovnice konvekce - upwind schéma Konkrétní hodnoty jsou zvoleny následovně: počáteční podmínka je skoková funkce 1 pro x < 0, 3 u(x, 0) = (3.108) 0 pro x > 0, 3 a okrajová podmínka v x = 0 je u(0, t) = 1. Parametry volíme tak, aby bylo vidět pokud možno co nejlépe chování výsledku. Pro všechny níže popsané příklady použijeme hodnoty a = 1 T = 0, 3

3.7.1

rychlost proudění čas výsledku

(3.109) (3.110)

Rovnice konvekce

V tomto případě přichází v úvahu pouze upwind schéma, na něm pak lze pozorovat závislost na hustotě diskretizační sítě a to konkrétně na počtu bodů N na ose x a poměru k zvoleného časového kroku k maximálnímu, který je dán podmínkou stability (3.89). Pro obě veličiny použijeme vždy dvě různé hodnoty: N = 20 resp. 200 k = 1 resp. 0, 5

(3.111) (3.112)


48

1 pocatecni podminka upwind N=20 k=0,5 upwind N=200 k=0,5 centralni N=200 k=0,5 0.8

0.6

0.4

0.2

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Obrázek 3.10: Srovnání centrálního a upwind schématu Výsledky jsou zobrazeny na obr. 3.9. Pro k = 1 se zachovává skok tak, jak byl v počáteční podmínce – schéma tedy počítá přesně bez ohledu na hodnotu N. Pro k = 0, 5 (i pro libovolné k < 1) dochází k zaoblování skoku, vlivem numerické difuze. Tento efekt je výraznější pro nižší počet bodů sítě, při zjemňování se výsledek blíží žádané skokové funkci, tak jak to vyplývá z teorie (konvergence numerického řešení ke skutečnému).

3.7.2

Konvekčně-difuzní rovnice

Pro testování použijeme stejné kombinace hodnot N a k, jako pro rovnici konvekce, koeficient difúze volíme b = 0, 005 (3.113) Nyní se však kromě pozorování závislosti na N, k nabízí i vzájemné srovnání centrálního a upwind schématu. Z obr. 3.10 vidíme, že i za přítomnosti skutečné difúze zaobluje upwind schéma původní skok více při nižším počtu bodů sítě. U centrálního schématu je pak zaoblení poněkud nižší než u srovnatelného upwind schématu, našimi prostředky však nelze rozlišit, nakolik je způsobeno fyzikální difúzí (koeficientem b) a nakolik samotným schématem, tj. numerickou difuzí. U centrálního schématu pak nastává zajímavý jev v souvislosti se vzájemným vzta-


49

1.4 pocatecni podminka N=200 k=1 N=20 k=1 N=20 k=0,5

1.2

1

0.8

0.6

0.4

0.2

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Obrázek 3.11: Efekty u centrálního schématu hem parametrů a, b, ∆x. Je-li totiž |a|∆x > 2, což v našem případě nastává pro N = 20, b pak v blízkosti skoku vznikají kmity, které se při nižším ∆t poněkud zmenšují (viz závislost chování na Pecletově čísle, sekce 3.3.1). Srovnání pro různé časové a prostorové kroky je na obr. 3.11.

3.7.3

Použití interaktivního programu JBONE

Pro názorné otestování vlastností různých numerických metod a schémat je k dispozici interaktivní program JBONE, který vytvořil prof. André Jaun na Royal Institute of Technology, Stockholm (kurz parciálních diferenciálních rovnic http://pde.fusion.kth.se/). Jeho rozšíření pro předmět Transportní procesy (tento učební text) je možno spustit z adresy http://www.fm.vslib.cz/~kmo/czech/vyuka/trp/jb/run-moje.html. S trochou šikovnosti přímo ze stránky nebo na žádost u autora (pošlu emailem) lze program stáhnout a poté spustit bez připojení na internet. Ovládání je stručně popsáno na zmíněné intenetové adrese a ve většině případů by měl být význam voleb a parametrů zřejmý. První volba je typ úlohy – volíme Advection (= konvekčně-difuzní rovnice) nebo Burgers (viz kapitola 5). Druhá volba je tvar počáteční podmínky (lze dále měnit parametry IC1AMP, IC1POS a IC1WID). Třetí volba je numerické metoda – použijeme Finite differences nebo Finite elements. Čtvrtá volba


50

je varianta numerického schématu (pro různé numerické metody ve volbě 3 je nabídka 4 různá) – schémata zmíněná v tomto textu jsou označena TRP. V pravé části jsou hodnoty parametrů, jejichž význam je zřejmý. Numerická difuze a stabilita u upwind schématu konvekce Pro výpočet použijeme následující sadu parametrů: parametr hodnoty

∆t a 0, 08 1 0, 8 1 1, 2 1

b ∆x obrázek 0 0, 5 3.12 0 0, 5 3.13 0 0, 5 3.14

V prvním případě je splněna podmínka stability a je zároveň je velká numerická difuze (vztah Dnum = 12 a∆x(1 − σ)) – je viditelný rozdíl mezi přesným (modrá křivka) a vypočteným řešením (černá křivka). V druhém případě je rovněž splněna podmínka stability (v tomto případě jsme blíže k mezní hodnotě která je ∆t = 1) a zároveň numerická difuze je menší než v předchozím případě. V třetím případě podmínka stability splněna není – z náhodných šumů vznikají oscilace (viz obrázek), které se jak je možno si vyzkoušet v dalším průběhu výpočtu zesilují do nekonečna (dojde k přetečení a chybě). Schémata pro konvekčně-difuzní rovnici Porovnáme explicitní (upwind aproximace konvekce), implicitní a Crank-Nicholsonovo schéma (obě centrální aproximace konvekce). Tabulka parametrů je: schéma explicit upwind implicit centrální Cranck-Nicholson

∆t a b ∆x obrázek 0, 55 1 0, 5 0, 5 3.15 0, 55 1 0, 5 0, 5 3.17 0, 55 1 0, 5 0, 5 3.16

V prvním případě je schéma pro daný časový krok nestabilní (vidíme vznikající oscilace), zatímco zbylá schémata pro tytéž hodnoty stabilní jsou (nepodmíněně stabilní). U implicitního schámatu je vidět jistá numerická difuze, Crank-Nicholsonovo schéma je druhého řádu přesnosti v prosotru i čase a numerická difuze u něho není.


Obrázek 3.12: Numerická difuze u upwind schématu, malý časový krok.

Obrázek 3.13: Numerická difuze u upwind schématu, velký časový krok.

Obrázek 3.14: Upwind schéma, velký časový krok, nestabilní.

51


Obrázek 3.15: Upwind schéma, nestabilní.

Obrázek 3.16: Implicitní centrální schéma.

Obrázek 3.17: Crank-Nicholsonovo schéma.

52

Kapitola 4

TRANSPORTNÍ PROCESY V PORÉZNÍM PROSTŘEDÍ Úvod V této kapitole uvádíme základní popis transportních procesů v porézním prostředí – vysvětleny a odvozeny jsou rovnice popisující proudění roztoku, transport rozpuštěných látek advekcí, hydrodynamickou disperzí a s vlivem sorpce. Výklad je podán pro obecné úlohy a jako příklady aplikace jsou uváděny úlohy týkající se podzemních vod.

4.1

Definice spojitého popisu porézního prostředí

Porézním prostředím (materiálem) nazýváme strukturu složenou ze zrn nebo vláken pevné látky (matrice) mezi nimiž je volný prostor (póry), který může být zaplněn vzduchem nebo kapalinou. V tomto textu se zabývame pouze situací, při níž je celý pórový prostor zaplněn kapalinou – roztokem (saturované porézní prostředí, jednofázové proudění). Přesný popis dějů v porézním prostředí by byl z důvodu mikroskopické struktury velmi složitý, proto se pracuje s aproximací porézního prostředí jako spojitého (kontinua). Všechny veličiny jsou pak uvažovány vždy jako průměr přes určitý objem. Pracujeme s pojmem reprezentativního elementárního objemu (REV), což je objem dostatečně velký na to, aby pokryl mikroskopickou nehomogenitu (tj. rozměr řádově větší než charakteriskický rozměr zrn pevné látky) a zároveň dostatečně malý vzhledem k rozměru zkoumané oblasti. Definice REV nemusí být vždy jednoznačná, např. v případě víceúrovňové mikroskopické struktury či různého typu makroskopických heterogenit. Obecně je každá veličina (koncentrace, tlak, rychlost, . . . ) v určitém bodě prostoru definována vztahem typu Z 1 αmic dV , (4.1) α= V V

kde V je objem REV se středem v uvažovaném bodě a αmic je hodnota příslušné veli-

Kapitola 4. Transportní procesy v porézním prostředí

54

činy na mikroskopické úrovni. Podle typu veličiny se do objemu zahrnuje buď jen volný prostor mezi pevnou matricí (póry), nebo jen pevná matrice případně celý objem prostředí. V případě rychlosti tak můžeme dostat dvě různé veličiny, jak je uvedeno níže: průměrnou rychlost v pórech a Darcyovskou rychlost (hustotu toku). Stejnou úvahou je pak definována porozita – poměr objemu pórů vzhledem k celému objemu materiálu objem pórů v REV n= (4.2) objem REV pro REV obklopující daný bod prostoru.

4.2

Proudění kapaliny v porézním prostředí

Proudění kapaliny (nebude-li uvedeno jinak, uvažujeme vodní roztok) v plně saturované prostředí je popsáno dvěma diferenciálními rovnicemi - pohybovou rovnicí a rovnicí bilance hmoty (rovnice kontinuity). Pohybovou rovnicí, která určuje závislost pohybu kapaliny na silovém půusobení, je v případě proudění v porézním prostředí Darcyho zákon, přesněji jeho zobecnění pro vektorové veličiny ve 3D.

4.2.1

Darcyho experiment a základní veličiny

V původní podobě byl Darcyho zákon formulován na základě 1D experimentu s prouděním vody podélným porézním filtrem v trubici. Experiment stručně popíšeme, jelikož nám umožní názorně definovat veličiny používané pro popis proudění v porézním prostředí – piezometrickou výšku a tzv. Darcyovskou rychlost (přesněji hustotu toku). Máme-li např. trubku naplněnou porézním materiálem umístěnou v obecně šikmé poloze, spojující dvě nádrže s různou výškou hladiny vody. Pro průtok vody trubkou Q (objem za čas) pak platí S · (φ2 − φ1 ) Q=K· , (4.3) L kde K je koeficient propustnosti (závislý na vlastnostech porézního materiálu i kapaliny), S příčný průřez trubky, L délka trubky a člen v závorce rozdíl výšky hladiny v obou nádržích. Symbol φ jako veličina přiřazená k místům na koncích trubky je pak výše zmiňovaná piezometrická výška, která hraje roli potenciálu. Názorně jde o výšku volné hladiny nad daným místem (tj. do jaké výšky by vystoupila hladina ve vrtu provedeném v daném místě) a je definována vztahem p φ=z+ , (4.4) %g kde z je svislá souřadnice, p tlak, % hustota a g tíhové zrychlení. Ze vztahu je zřejmé, že piezometrická veličina v sobě zahrnuje potenciál tíhového pole (z-ová poloha) a potenciál tlaku (hydrostatický tlak vody nad daným místem).


55

p Samotný člen %g se nazývá tlaková výška (značíme h) a jde v podstatě o tlak vyjádřený v délkových jednotkách (tlak jako výška vodního sloupce), což je praktické zvláště v hydrogeologických aplikacích. Dále přejdeme k lokálnímu popisu proudění vody, jt. vyjádření pohybu v daném místě. Definujeme plošnou hustotu toku vody jako podíl množství protékající vody a (pro výše uvedený příklad). Uvažujeme-li velikosti plochy kolmé na směr proudění q = Q S φ2 −φ1 i limitu v podélném směru, ∇φ = L , dostaneme Darcyho zákon ve tvaru

q = K∇φ .

(4.5)

Veličina q má rozměr rychlosti, proto se také běžně nazývá Darcyovská rychlost, nejde však o makroskopický pohyb nějakého bodu touto rychlostí. Odvodíme nyní, jakou rychlostí se bude pohybovat zvolená částice vody, případně z hlediska níže popsaného transportu konvekcí, jakou rychlostí se bude pohybovat rozpuštěná látka. Uvažujeme úsek trubky z Darcyho experimentu délky L, který má celkový objem V = SL a obsahuje vodu o objemu Vw = n · SL (n je porozita). Při dané hustotě toku (resp. Darcyovské rychlosti) q, proteče tento objem vody kolmým průřezem za dobu t=

Vw nL = Sq q

(4.6)

a je tedy vidět, že se po délce L posunul rychlostí v=

q , n

(4.7)

která je jak je vidět vyšší než Darcyovská rychlost q a nazývá průměrná rychlost vody v pórech, protože jde zároveň o veličinu, kterou bychom vypočetli standardním průměrováním přes REV (jako vektor se zahrnutím obecných směrů pórů): Z 1 v= w v (mic) dV w . (4.8) VREV w VREV

4.2.2

Darcyho zákon ve 3D

V obecném případě jsou piezometrická výška a vektor Darcyovské rychlosti funkce prostorových souřadnic a času a platí mezi nimi vztah q = −K∇φ ,

(4.9)

který budeme pod pojmem Darcyho zákon používat v dalším textu. Připomínáme vý∂φ ∂φ , ). Protože má porézní prostředí obecně znam symbolu ∇ jako gradientu, ∇φ = ( ∂φ ∂x ∂y ∂z


56

anizotropní charakter, je koeficient hydraulické vodivosti tenzor druhého řádu se složkami   Kxx Kxy Kxz (4.10) K =  Kyx Kyy Kyz  , Kzx Kzy Kzz

které jsou v případě nehomogenního prostředí navíc funkcí prostorových souřadnic. Často se vyskytuje situace, kdy hlavní směry anizotropie přibližně odpovídají kartézským souřadnicím (např. vrstevnaté půdní sedimenty) a Darcyho zákon lze použít v podobě soustavy qx = −Kx ∂φ ∂x qy = −Ky ∂φ (4.11) ∂y , ∂φ qz = −Kz ∂z kde Kx , Ky , Kz jsou vodivosti ve směrech os. Koeficient hydraulické vodivosti závisí na vlastnostech porézního prostředí (především mikroskopické geometrické struktuře) i vlastnostech kapaliny. Lze zavést takzvanou propustnost (permealilitu) k (s rozměrem metr), která již závisí jen na porézním materiálu a platí k. % . g K= , (4.12) µ kde % je hustota a µ je dynamická viskozita kapaliny [kg/m/s]. Poznámka: Platnost lineárního Darcyho zákona je omezena zdola i shora. Při velmi malých gradientech tlaku se projeví molekulární síly, které proudění zpomalují nebo úplně znemožní (tj. k pohybu podle Darcyho zákona dojde až po překročení určité prahové hodnoty gradientu tlaku). Při velkých rychlostech ztrácí proudění laminární charakter a lineární závislost nepopisuje děj s dostatečnou přesností. Ve většině případů je však Darcyho zákon plně vyhovující.

4.2.3

Rovnice bilance hmoty

V libovolně zvoleném objemu musí platit, že změna hmotnosti vody bude odpovídat hmotnosti vody prošlé přes hranici a změně v rámci propadů a zdrojů, tj. Z Z Z ∂ ρn dV = − ρq · dS + P ρ dV , (4.13) ∂t V

∂V

V

kde P je hustota zdrojů (+) či propadů (-) vyjádřená jako objem kapaliny vtlačený do jednotkového objemu porézního materiálu za jednotkový čas [m3 /m3 /s] (specifická


57

vydatnost zdroje). Standartní úpravou přes Gaussovu větu, převedení členů na jednu stranu, dostáváme nulovost integrandu, po úpravě ∂(%n) + div(%q) = P % . ∂t

(4.14)

Člen s časovou derivací vyjadřuje změnu hustoty a porozity, které mohou nastat např. vlivem stlačitelnosti porézního prostředí (často existují izolované póry, které zůstanou zaplněny vzduchem), závislosti hustoty kapaliny na rozpuštěných látkách a vlivem teplotní roztažnosti. V případě proudění v nesaturovaném prostředí by v tomto členu přibyla i saturace – podíl vody v celkovém objemu pórů. V tomto textu pro účely popisu transportu rozpuštěných látek budeme nadále uvažovat zjednodušený model s konstantní hustotou a porozitou. Rovnice (4.14) se pak zjednoduší na P div v = . (4.15) n Vidíme, že ani v jedné z popisujících rovnic (4.9) a (4.15) se nevyskytuje čas, jde tedy o úlohu ustáleného (stacionárního) proudění. Za podmínek nestlačitelnosti, konstantní hustoty a saturovaného prostředí totiž neustálená úloha ani nemá smysl – při změně okrajových podmínek (např. tlaku v určitém místě) by se změna musela projevit okamžitě v celé řešené oblasti (prostředí má nekonečnou tuhost a vzruch by se šířil nekonečnou rychlostí). V praktických aplikacích se často úloha formuluje jako série ustálených stavů (při změně vnějších podmínek vznikne nový ustálený stav relativně rychle). Rovnice (4.9) a (4.15) můžeme ponechat ve tvaru systému dvou rovnic 1.řádu pro dvě neznámé funkce q, φ nebo dosazením upravit na jednu rovnici druhého řádu ∇ · (K∇φ) = P

(4.16)

pro neznámou funkci φ(x), x = (x, y, z).

4.2.4

Okrajové a počáteční podmínky

Kromě řídících diferenciálních rovnic, tvaru a struktury řešené oblasti je pro formulaci úlohy nezbytné zadání okrajových a počátečních podmínek. Pro ustálené prouděné (stacionární) mají smysl pouze okrajové podmínky a ty je nutno zadat na celé hranici Γ. Budeme užívat standartní značení, kdy Ω ⊂ R3 označuje oblast, ve které úlohu řešíme (definiční obor neznámých funkcí v, φ) a Γ = ∂Ω její hranici. V úlohách proudění se uplatňují okrajové podmínky všech tří typů:


58

Dirichletova okrajová podmínka (stabilní okrajová podmínka, podmínka prvního typu) udává hydraulickou výšku na části hranice Γ1 ⊂ Γ φ(x) = φD (x)

∀x ∈ Γ1 .

(4.17)

V případě podzemního proudění se hodnota zadává na svislé části hranice, piezometrické výšky jsou známy z monitorovacích vrtů provedených podél hranice. Neumannova okrajová podmínka (nestabilní okrajová podmínka, podmínka druhého typu) udává specifický průtok na části hranice Γ2 ⊂ Γ q · ν = qN (x) ,

(4.18)

kde ν je jednotkový vektor normály k hranici v daném místě x a qN hodnota specifického průtoku v příslušném bodě hranice. Podmínka se uplatňuje nejčastěji v případě nepropustné hranice, kdy je průtok roven nule a platí tedy q · ν = 0, (4.19)

což je homogenní Neumannova podmínka. Newtonova (též Cauchyho) okrajová podmínka (smíšená okrajová podmínka) se používá v případě, že je část hranice Γ3 ⊂ Γ tvořena polopropustnou vrstvou. Velikost specifického průtoku je úměrná rozdílu hydraulických výšek na obou stranách polopropustné překážky q · ν = (φ − φ0 )/C , (4.20) kde φ

je

φ0 C = B/K B K

hydraulická výška uvnitř oblasti [m] (u hranice) hydraulická výška vně oblasti [m] odpor překážky vůči proudění [1/s] její šířka [m] její hydraulická vodivost [m/s]

Je-li proudění neustálené, je třeba zadat též počáteční podmínky. Ty popisují proudění ve sledované oblasti na počátku řešení, tedy v čase t = 0. Jde o zadání hydraulické výšky ve všech bodech sledované oblasti. Okrajové podmínky mohou být v takovém případě též funkcemi času. Jsou-li okrajové podmínky konstantní, je předmětem řešení přechodový děj, tzn. přechod od počátečního stavu ke stavu ustálenému.

4.3

Transport rozpuštěných látek

Transport rozpuštěných látek v porézním prostředí probíhá mnoha mechanismy, z nichž některé probereme jednotlivě v následujících podkapitolách. Mezi základní procesy patří


59

advekce (často se užívá též termín konvekce), disperze, sorpce, chemické reakce a radioaktivní rozpad. Další jevy, které ovlivňují transport látky, se pak objevují též vlivem speciální struktury pevné matrice, typicky v případě dvojí porozity (viz níže). Množství látky v roztoku se popisuje pomocí koncentrace c, vyjádřené jako hmotnost rozpuštěné látky připadající na jednotkový objem roztoku (vody). V některých případech je výhodné vyjádřit množství látky vzhledem k hmotnosti vody (hmotnostní koncentrace), případně vzhledem k objemu celého porézního média. Zmíněné veličiny si lze snadno odvodit, přesná vyjádření uvedeme v sekci sorpce, kde budou nutné k sestavení rovnic.

4.3.1

Advekce a disperze

Advekce1 znamená přenos látky vlivem pohybu celého roztoku. Množství přenesené látky můžeme snadno vyjádřit z darcyovské rychlosti proudění: hmotnost látky prošlá jednotkovou plochou za jednotkový čas je q adv = cq. c Difuzně-disperzní děje jsou takové, kdy dochází k pohybu rozpuštěné látky vlivem gradientu koncentrace, z míst vyšší koncentrace do míst nižsí koncentrace, tj. např. dochází k rozmývání „hranÿ mezi roztokem a sousedící čistou vodou. V porézním prostředí se uplatňují procesy v této podobě: • molekulární difuze (rovněž ovlivněná mikroskopickou strukturou prostředí, ale založená na stejném principu) • mechanická disperze (přímo způsobená nehomogenitou rychlosti v pórech – v některých místech se roztok pohybuje rychleji a v některých pomaleji než je průměrná rychlost v) Mechanická disperze „vznikáÿ až jako důsledek přechodu od mikroskopického popisu k makroskopickému – na mikroskopické úrovni by šlo pouze o advekci a molekulární difuzi, ale v mnohem složitějším rychlostním poli (případně o „ jinouÿ disperzi vlivem turbulencí v proudění). Intenzita molekulární difuze nezávisí na proudění a je dána hodnotou difúzního koeficientu Dm . Ta je pro různé látky odlišná a obecně může záviset na dalších vlivech, např. teplotě. V čisté vodě je pro většinu běžných látek Dm v řádu 10−9 m2 /s. Samotná difuze ve volné kapalině je popsána Fickovým zákonem, který má tvar X ∂2c ∂c = Dm ∆c = Dm ∂t ∂x2i i

(4.21)

se skalárním koeficientem Dm (tj. děj probíhá izotropně). V porézním prostředí je difuze ovlivněna geometrií pórů, vyjádřenou tenzorem tortuozity T , který má pro izotropní 1

Užívá se též pojem konvekce, ve významu užívaném v tomto textu je lze považovat za ekvivalentní.


60

prostředí diagonální tvar se stejnými hodnotami. Zavádí se tenzor molekulární difuze D m = Dm T . Mechanická disperze je popsána stejnou rovnicí jako molekulární difúze, ale s jinak vyjádřeným koeficientem. Nazývá se tenzor mechanické disperze Df a je závislý na rychlosti proudění podle vztahu [Df ]ij = αT .|v|.δij + (αL − αT )

vi vj , |v|

(4.22)

kde koeficienty αL a αT se nazývají podélná a příčná disperzivita (longitudinal/ transversal). Situaci si lze nejsnáze představit tak, že disperze probíhá ve směru proudění s koeficientem DL = αL v a ve směru kolmém na směr proudění s koeficientem DT = αT v (jako jednodimenzionální děj), přičemž αT je zhruba o řád menší než αL . Celkový disperzní děj, složený z molekulární difuze a mechanické disperze se nazývá hydrodynamické disperze a je charakterizována tenzorem hydrodynamické disperze Dh , definovaným D = Dm + Df [m2 /s] . (4.23) Vzájemný poměr vlivu molekulární difuze a mechanické disperze na celkovou disperzi může být v různých situacích velmi odlišný. Záleží na vlastnostech porézního média, řešené úloze (okrajové podmínky → rychlostní pole) i měřítku pozorování. Charakter disperze udává Pecletovo číslo vd , (4.24) Pe = Dm kde v d Dm

je

střední pórová rychlost roztoku střední velikost zrna koeficient molekulární difuze

[m s−1 ] [m] [m2 s−1 ]

Malé Pecletovo číslo (Pe<0.01) odpovídá převládajícímu vlivu molekulární difuze, velké Pecletovo číslo (Pe>104 ) odpovídá převládající mechanické disperzi. Běžně se molekulární difuze zanedbává pro Pe>20. Nyní odvodíme řídící rovnici popisující advekčně-disperzní transport. Celkový hmotnostní tok vlivem advekce a hydrodynamické disperze je q c = n(cv + D h ∇c) .

(4.25)

Pro bilanci látky v roztoku platí (vyjádřenou vzhledem k celkovému objemu prostředí) platí Z Z Z Z ∂ + ∗ − nc dV = − q c · dS + (P c + P c) dV + r dV v libovolném V, (4.26) ∂t V

∂V

V

V


61

kde P + a P − vyjadřuje kladnou a zápornou část funkce hustoty zdrojů roztoku (rovnice (4.14)) a vyjadřují skutečnost, že je vtláčen roztok se zadanou koncentrací c∗ , zatímco čerpán je roztok s koncentrací c odpovídající hledané hodnotě funkce c(x, t) v daném místě. Člen r vyjadřuje množství látky vznikající nebo vstupující do roztoku vlivem dalších dějů – sorpce, chemických reakcí a radioaktivního rozpadu (hmotnost na jednotkový objem porézního média za jednotkový čas). Přirozeně může mít kladné i záporné znaménko v závislosti na „směruÿ děje. Standartní úpravou intergrálů dostáváme diferenciální rovnici 1 r ∂c = −∇ · (cv) + ∇ · (D h ∇c) + (P + c∗ + P − c) + , ∂t n n

(4.27)

která se nazývá advekčně-disperzní (případně advekčně-difuzní nebo konvekčně-difuzní) rovnice. Často se uvádí s advekčním a disperzním členem na levé straně, potom pro nulovou pravou stranu označujeme transport za konzervativní (zachovává se bilance látky) a pro nenulovou pravou stranu transport nekonzervativní (tj. se zdroji a propady různého charakteru). Rovnici lze rovněž upravit derivováním součinu v advekčním členu a dosazením za ∇ · v z rovnice bilance proudění vody (4.15), kdy dojde k vyrušení členů odpovídajících čerpání, na tvar ∂c P+ ∗ r + v · ∇c − ∇ · (D h ∇c) = (c − c) + . ∂t n n

4.3.2

(4.28)

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

Transportní úloha je vždy nestacionární (hledaná koncentrace je funkcí času). Proto je třeba zadat počáteční podmínky, což je pro naši úlohu rozložení koncentrace látky na počátku řešeného časového intervalu c(x, 0) = c0 (x) ∀x ∈ Ω .

(4.29)

V případě nekonzervativního transportu (vlivu chemických reakcí, nerovnovážné sorpce) je třeba obdobným způsobem zadat rovněž koncentrace látek reagujících se sledovanou látkou. Dirichletova okrajová podmínka Dirichletova okrajová podmínka (podmínka 1. druhu) udává koncentraci na dané části hranice Γ1 c(x, t) = cD (x, t) ∀x ∈ Γ1 , (4.30)


62

kde cD je zadaná funkce. Tento typ podmínky se uplatňuje např. když oblast v daném místě sousedí dobře míseným roztokem se známou koncentrací a lze předpokládat, že bude docházet k rychlému vyrovnání koncentrací na obou stranách hranice. V případě procesů s dominantním vlivem advekce se tato podmínka předepisuje na části hranice, kde je směr proudění roztoku směrem dovnitř a zadává se tedy koncentrace přímo vstupující do sledované oblasti. Neumannova okrajová podmínka Narozdíl od úlohy proudění nelze v případě advekčně-disperzního transportu přímo ztotožnit situaci předepsaného toku a Neumannovy podmínky. Tok rozpuštěné látky je dán vztahem (4.25), tj. q c = n(cv + Dh ∇c) , (4.31)

kde se vyskytuje jak hledaná koncentrace, tak její derivace a obecně by tedy šlo o podmínky 3.druhu. Neumannova podmínka se tedy uplatňuje pouze v případě, kdy je zadán pouze disperzní tok (D h ∇c) · ν = qdisp ∀x ∈ Γ2 . (4.32) Podmínka v této formě se uplatňuje nejčastěji v případě, kdy je tok zadán nulový (homogenní Neumannova podmínka) a to buď na nepropustné hranici nebo u úloh s dominantní advekcí na hranici s prouděním roztoku směrem ven, kdy se použije předpoklad, že koncentrace těsně za hranicí je přibližně stejná. Newtonova (Cauchyho) okrajová podmínka Jde o obecně formulovanou podmínku předepsaného toku (cv + D h ∇c) · ν = qdisp

∀x ∈ Γ3 ,

(4.33)

z níž lze odvodit podmínky pro různé specifické situace. V případě nepropustné hranice, při nulovém zadaném toku, dostáváme podmínku stejného tvaru jako při použití Neumannovy podmínky nulového disperzního toku, neboť rychlost proudění (a tedy i advekce) přes nepropustnou hranici je automaticky nulová.

4.3.3

Sorpce

Proces sorpce je příkladem interakce mezi roztokem a pevnou matricí (např. horninou). Jde o zachycování (ulpívání) migrující látky na povrchu pevné fáze (adsorpce) a její zpětné uvolňování do roztoku (desorpce). Obdobný charakter mají též rozpouštěcí a srážecí reakce.


63

Nejprve zavedeme veličiny, které budeme používat pro určení množství zachycené látky. Hmotnost zachycené látky lze normovat buď na hmotnost pevné fáze nebo na objem, a to buď jen samotné pevné fáze nebo celého porézního média. V literatuře se používají různá vyjádření i různá značení, my budeme uvádět veličiny v té formě, aby se v rovnicích vyskytoval pokud možno co nejmenší počet koeficientů. Označíme s hmotnost sorbované látky na jednotkový objem pevné fáze s˜ hmotnost látky na jednotkový objem porézního prostředí, s˜ = (1 − n)s s s˜ s¯ hmotnost látky na jednotkovou hmotnost pevné fáze, s¯ = = %s %b kde %s označuje hustotu pevné fáze (solid) a %b = (1 − n)%s měrnou hmotnost suchého porézního prostředí (bulk density). Jiná označení pro koncentraci sorbované látky jsou např. c¯ či cs . V dalším textu bude v některých případech výhodné uvádět i koncentraci v roztoku vzhledem k celkovému objemu prostředí, v souladu s výše uvedeným značením zavedeme c˜ = nc (celkové množství látky v daném objemu je tedy potom c˜ + s˜) Základní charakteristikou popisující sorpci je vztah mezi koncentrací látky v roztoku a koncentrací látky sorbované. Ten má různou podobu, v závislosti na charakteru procesu a míře zjednodušení v modelu. Z hlediska rychlosti procesu rozlišujeme rovnovážnou nebo nerovnovážnou sorpci. Rovnovážný případ znamená, že např. při změně koncentrace látky v roztoku (přitečením roztoku z jiného místa) dojde k ustavení rovnováhy mezi rozpuštěnou a sorbovanou látkou velmi rychle a lze tedy přímo definovat určitou funkci udávající závislost sorbované koncentrace na koncentraci v roztoku, která neobsahuje čas (sorbovaná koncentrace nezávisí na stavu v dřívějším čase, ale jen na okamžité hodnotě c). Nerovnovážný případ znamená, že změny koncentrace vlivem transportu jsou srovnatelně rychlé jako přesun látky mezi roztokem a pevnou fází: koncentrace na pevné fázi a v roztoku pak nejsou v rovnováze a nelze tedy přímo určit sorbovanou koncentraci v daném okamžiku (závisí též na předchozím stavu). V tomto případě je sorpce popsána velikostí toku látky mezi oběma fázemi, který je funkcí obou koncentrací c i s. Popis rovnovážné sorpce – izotermy Vztah mezi sorbovanou koncentrací s a koncentrací v roztoku c v rovnovážném případě se nazývá sorpční izoterma (podle toho, že se závislost popisuje za předpokladu konstatní teploty). Závislost je obecně složitá a nejčastěji se aproximuje těmito vztahy (obvykle vyjádřované pro s¯): • lineární izoterma s = kD c = %s KD c (tj. s¯ = KD c) • Freundlichova izoterma s¯ = KF ca • Langmuirova izoterma s¯ =

KL s¯max c 1 + KL c


64

Je zřejmé, že lineární izoterma je speciálním případem Freundlichovy. Koeficient KD se nazývá rozdělovací (distribuční) koeficient. Lagmuirova izoterma respektuje to, že může dojít k jistému „nasyceníÿ pevné fáze a koncentrace nepřekročí jistou hodnotu s¯max . Odvodíme nyní rovnici transportu s vlivem sorpce (pro jednoduchost uvažujeme nyní lineární sorpci). Úvahu je možno provést dvěma způsoby: buď vezmeme již odvozenou advekčně-disperzní rovnici (4.28) a vyjádříme produkční člen r pomocí sorpční izotermy nebo vyjádříme bilanci ve formě rovnice (4.26) společně pro roztok i pevnou fázi (na základě této úvahy pak někdy označujeme transport s vlivem sorpce rovněž jako konzervativní). Odvození přes zdrojový člen r: V prvním případě je třeba vyjádřit intenzitu přechodu látky mezi fázemi, což odpovídá časové změně koncentrace. Úbytek sorbované látky rozpočtený na celý objem je ∂(s(1 − n)) ∂c = −kD (1 − n) ∂t ∂t a na jednotkový objem roztoku tj. již jako zdrojový člen v ADE r˜s = −

(4.34)

1 − n ∂c r˜s = −kD . (4.35) n n ∂t Po dosazení vidíme, že ve zdrojovém členu se vyskytuje „tatážÿ časová derivace jako na levé straně rovnice, kde vyjadřuje bilanci v roztoku. Sloučením obou členů dostaneme 1 − n ∂c 1 1 + kD = −∇ · (cv) + ∇ · (Dh ∇c) + (P + c∗ + P − c) . (4.36) n ∂t n Součinitel u časové derivace se nazývá retardační faktor rs =

1−n (4.37) n a vyjadřuje zpoždění advekčně-disperzního děje vzhledem ke stavu, který by probíhal bez vlivu sorpce. Z rovnice je zřejmé, že výsledky budou stejné, jako kdyby rychlost proudění byla Rv (v souladu s odpovídajícím způsobem sníženou intenzitou zdrojů P/R) a disperzní koeficient byl DRh . R = 1 + kD

Odvození přes bilanci v celém objemu porézního prostředí: Druhý způsob je v podstatě zopakování odvození advekčně-disperzní rovnice (4.28) s tím rozdílem, že ve výchozím vztahu (4.26) člen udávající změnu množství látky v daném místě (kontrolním objemu V ) vyjádříme pro látku přítomnou v roztoku i na povrchu pevné fáze, tj. Z Z Z ∂ (˜ c + s˜) dV = − q c · dS + (P + c∗ + P − c) dV , (4.38) ∂t V

∂V

V


65

kde jsme již neuvažovali další produkční člen r (tj. další vlivy kromě sorpce). Dosazením za c˜ = nc a s˜ = (1 − n)s = (1 − n)kD c dostaneme rovnici ve stejném tvaru ∂c 1 = −∇ · (cv) + ∇ · (D h ∇c) + (P + c∗ + P − c) , ∂t n kde R je výše definovaný retardační faktor. R

(4.39)

Nelineární sorpce V případě nelineární sorpce je možné zobecnit první způsob odvození, vztah (4.34) – vyjádření produkčního členu rs v advekčně-disperzní rovnici. Pro obecnou závislost (izotermu) s = f (c) pak platí r˜s = −

∂s ∂c ∂(s(1 − n)) = −(1 − n) · · , ∂t ∂c ∂t

(4.40)

kde, narozdíl od vztahu (4.35), místo konstanty kD lineární sorpční izotermy obecně vystupuje derivace izotermy v daném bodě (pro příslušnou koncentraci c). Stejným postupem pak lze definovat retardační faktor 1 − n ∂s · , (4.41) n ∂c který však již není konstatní, ale závislý na koncentraci c a transport látky je pak popsán nelineární parciální diferenciální rovnicí. R=1+

Nerovnovážná sorpce Proces nerovnovážné sorpce lze určit závislostí intenzity přechodu látky mezi roztokem a pevnou fází, přesněji hmotnosti látky za jednotku času v jednotkovém objemu porézního prostředí, na koncentraci v roztoku c a sorbované koncentraci s. Závislost může být obecně velmi složitá, v nejjednodušší formě ji lze vyjádřit jako lineární na rozdílu koncentrace a odpovídajícího rovnovážného stavu: r˜s = (1 − n) α (s − kD c) .

(4.42)

Vývoj rozložení látky je popsán dvěma funkcemi – koncentrací v roztoku c(x, t) a sorbovanou koncentrací s(x, t). Dosazením výše uvedeného vztahu pro intenzitu výměny do advekčně difuzní rovnice pro c a vyjádřením bilance sorbované látky dostáváme systém rovnic ∂c 1−n = −∇ · (cv) + ∇ · (D h ∇c) + α (s − kD c) , ∂t n ∂s = −α (s − kD c) . ∂t

(4.43) (4.44)


4.3.4

66

Vliv chemických a dalších reakcí

Obecný souhrnný popis všech možných vlivů na roztok v porézním prostředí je nad rámec tohoto textu. Výběr sledovaných dějů velmi závisí na konkrétní aplikaci. Obvykle se zmiňuje vliv radioaktivního rozpadu, biodegradace a chemických reakcí těchto typů: • oxidace / redukce • acidobázické procesy • srážení • rozpouštění • tvorba komplexů • substituce / hydrolýza Chemické reakce lze z hlediska rovnováhy rozdělit na rovnovážné a nerovnovážné, podobně jako výše podrobněji popisovaný proces sorpce. Matematický popis procesů je obecně složitý, z důvodu kombinace mnoha vlivů (teplota, pH, přítomnost dalších látek). Jako základní aproximace se obvykle používají dva reakční modely - reakce nultého a prvního řádu ∂c = k0 , ∂t ∂c r= = k1 c . ∂t

r=

(4.45) (4.46)

Typickým případem reakce prvního řádu je radioktivní rozpad, kdy r = −λc

λ=

ln 2 , T1/2

(4.47)

kde T1/2 je poločas rozpadu.

4.4

Motivace: aplikace modelů na úlohy podzemního proudění

Fyzikální a chemické procesy, probíhající v horninovém prostředí, mají velký význam při řešení ekologických problémů i v technické praxi. Ekologická problematika zahrnuje čištění kontaminovaných vod včetně vytváření systémů, zabraňujících šíření kontaminace a průniku kontaminantů do životního prostředí. Z technických oborů mají podzemní procesy značný význam při těžbě nerostů, při využívání zemského tepla a při stavbě vodohospodářských děl.


67

Fyzikální procesy zahrnují hlavně pohyb podzemní vody, s nímž souvisí i transport rozpuštěných látek. Ten je kromě konvektivního přenosu zprostředkován též disperzí a difúzí. Migrace látek je ovlivněna interakcí mezi roztoky a horninami. Zde jsou nejdůležitější sorpční a desorpční procesy, patří sem však i chemické reakce, zejména oxidačněredukční děje, rozpouštění a srážení minerálů. Kromě toho probíhají chemické reakce i v roztocích v důsledku jejich míchání či ředění při postupu prostředím. Jiným typem fyzikálního procesu je šíření tepla v horninovém masivu. Teplo se přenáší jednak konvekcí (pohybem pórové vody), jednak vedením horninou. Může vznikat při rozpadu radionuklidů. Zpravidla zanedbatelný význam má uvolňování a spotřeba tepla při průběhu chemických reakcí. Abychom mohli procesy v horninovém prostředí řídit, musíme správně chápat jejich podstatu a umět prognózovat další vývoj. To není jednoduché, uvážíme-li, že o procesech, probíhajících často v rozsáhlých areálech, máme k dispozici jen omezený okruh informací, které mají často „bodovýÿ charakter v čase i v prostoru. Nezbytným prostředkem pro pochopení a kvantifikaci podzemních procesů je matematické modelování. Simulační modely popisují proudění podzemních vod, šíření tepla i migraci látek za stanovených podmínek. Pokud některé podmínky neznáme s dostatečnou přesností, můžeme alespoň posoudit míru jejich vlivu na sledovaný proces. Z hlediska řízení procesů je důležitá možnost provádění numerických experimentů. Zde nadefinujeme řídící zásah, převedeme jej do vstupních souborů modelu a zjistíme odezvu v podzemí. Porovnáním výsledků lze ocenit technologický efekt různých zásahů a získat vstupy do optimalizačních výpočtů. Zásahy do podzemí provádíme ze zemského povrchu. Kromě případů jednorázových opatření k odchýlení směru proudění podzemních vod (vytváření nepropustných překážek) navazují na vývoj v podzemí povrchové technologické procesy. Jedná se např. o separaci užitkových složek a zpětné vtláčení roztoků po doplnění potřebných chemikálií, o čištění čerpaných roztoků a vypouštění vyčištěných vod do recipientu, o získávání tepla, které roztoky přenášejí apod. Účinnost povrchových technologií je vázána na objem a vlastnosti zpracovávaných roztoků. To znamená, že je nutno projektovat vrtné práce a technické zařízení pro čerpání a transport roztoků tak, aby byly technologické kapacity efektivně využívány. Přitom je nutno počítat se změnami vydatnosti vrtů, koncentrací rozpuštěných a nerozpuštěných látek, teploty a dalších parametrů, které se mohou projevit v průběhu čerpání. K povrchovým činnostem se váží investiční a provozní náklady. Často jsme postaveni před problém porovnání ekonomické výhodnosti různých zásahů se srovnatelným technologickým efektem nebo obecně před ocenění ekonomicko-technologického efektu odlišných postupů. Pro tyto účely je vhodné sestavit technologické a ekonomické modely jako nadstavbu nad modely podzemních procesů. S jejich pomocí lze hodnotit přímo ekonomický dopad zamýšlených řídících zásahů.


68

Vlastnosti horninového prostředí Horninový masiv představuje nehomogenní a často anizotropní prostředí. Skládá se z pevné fáze, označované jako horninová matrice, a z volného prostoru. Horninová matrice sestává z minerálních zrn. Volný prostor představují póry, pukliny a různé dutiny. Volný prostor může být vyplněn vodou. V nezvětralých krystalických (tj. vyvřelých a metamorfovaných) horninách jsou minerální zrna většinou těsně uspořádána. Volný prostor zde reprezentují hlavně pukliny, vzniklé při chladnutí magmatu nebo při tektonických procesech. V sedimentárních horninách, s výjimkou chemicky (srážením) vzniklých sedimentů (např. vápence), je volný prostor tvořen hlavně póry. Póry vytvářejí systém vzájemně propojených i slepých kanálků různého průměru, délky a zakřivení. Proudění v porézních horninách označujeme jako průlinové, v rozpukaných jako puklinové. Přesný popis pohybu kapaliny v takovém systému by byl velmi komplikovaný. Naštěstí se bez něj obejdeme. Pro matematický popis pohybu podzemní vody chápeme horninový masiv jako kontinuum s určitými vlastnostmi, které zjistíme průměrováním daných veličin v dostatečně velké objemu. Reprezentativní elementární objem (REV) musí být dostatečně velký, aby se v něm již neprojevoval vliv lokální mikrostruktury, a přitom dostatečně malý, aby se do něj nepromítaly nehomogenity, vyplývající z geologické stavby oblasti. Podrobnější výklad lze nalézt v (Císlerová, Vogel, 1998). REV se může značně (až o několik řádů) lišit pro pórovité horniny a pro horniny se systémem puklin. Rozpukané horniny můžeme aproximovat kontinuálním prostředím zpravidla pouze v případě řešení úloh regionálního charakteru, kde rozměry oblasti mnohonásobně převyšují rozměry běžně se vyskytujících puklin. I v takovém případě však musíme významné pukliny interpretovat jako nehomogenity. Důležitými vlastnostmi hornin jsou pórovitost a hydraulická vodivost. Pórovitost vyjadřuje poměr mezi objemem volného prostoru a jednotkovým objemem horniny. V některých případech můžeme do pórovitosti zahrnout i prostor puklin. Slepé (neprůtočné, neaktivní) pórové kanály a pukliny se nepodílejí na proudění vody, mají však význam pro migraci rozpuštěných látek, které difundují mezi aktivními a neaktivními póry. Proto se v případě potřeby uvažuje dvojí pórovitost - aktivní a neaktivní. Orientační hodnoty pórovitosti pro některé zeminy a horniny jsou uvedeny v tabulce 1 (Císlerová, Vogel, 1998). Hydraulická vodivost (filtrační koeficient) K má rozměr rychlosti [m/s]. V izotropním prostředí odpovídá hustotě toku při jednotkovém hydraulickém spádu. Hydraulická vodivost závisí jak na charakteru horniny, tak na vlastnostech pohybující se kapaliny. Přesněji bude pojem hydraulické vodivosti vysvětlen níže, v souvislosti s Darcyho zákonem popisujícím rychlost proudění za různých podmínek.


Materiál štěrk hrubý štěrk jemný písek hrubý písek jemný prach jíl

Pórovitost (%) 24 - 36 25 - 38 31 - 46 26 - 53 34 - 61 34 - 60

pískovec prachovec vápenec, dolomit krasový vápenec

5 - 30 21 - 41 0 - 20 5 - 50

rozpukané krystalinikum nerozpukané krystalinikum zvětralá žula zvětralé gabro

0 - 10 0-5 34 - 57 42 - 45

69

Tab. 4.1: Orientační hodnoty pórovitostí Orientační hodnoty hydraulické vodivosti pro některé typy hornin a zemin jsou uvedeny v tabulce 2. Z porovnání údajů z tabulek 1 a 2 je zřejmé, že vyšší pórovitost neznamená vyšší propustnost. Hydraulická vodivost je totiž závislá také na vlastnostech povrchu pevné fáze. V jemnozrnných materiálech s velkou pórovitostí a s velkým měrným povrchem pevné fáze se významně uplatňují molekulární síly, které brání pohybu vody. Proto je jejich propustnost výrazně nižší než u materiálů hrubozrnných. Hydrogeologické struktury - kolektory a zvodně Hydrogeologickou strukturou rozumíme geologické prostředí, v němž dochází k akumulaci a pohybu podzemní vody. Přirozený pohyb vody směřuje z oblasti infiltrace k oblasti vývěru. Podmínkou je existence souvislého propojení propustných partií horninového prostředí. Hydrogeologická struktura se skládá z jednoho nebo několika kolektorů. Kolektor je propustný útvar (vrstva, souvrství), v němž se může hromadit a relativně snadno pohybovat podzemní voda. Těleso podzemní vody v kolektoru označujeme jako zvodeň. Na existenci a velikosti akumulační oblasti závisí doba zdržení elementárních objemů vody při jejich průchodu podzemím. Je nutno mít na zřeteli, že ve složitěji tvarovaných a strukturovaných oblastech se doby zdržení pro jednotlivé částice mohou lišit o několik


Druh horniny/zeminy štěrk hrubozrnný písek jemný písek jílovité zeminy jíl pískovec prachovitý pískovec prachovec rozpukané krystalické horniny kompaktní krystalické horniny

70

Hydraulická vodivost K [m/s] 10−3 - 10−2 10−4 - 10−3 10−5 - 10−4 10−9 - 10−6 <10−9 10−6 - 10−4 10−8 - 10−6 10−11 - 10−8 10−10 - 10−7 <10−13

Tab. 4.2: Orientační hodnoty hydraulické vodivosti řádů. Jestliže se v oblasti nachází více kolektorů, jsou odděleny hydrogeologickými izolátory, tvořenými nepropustnými horninami. Propustnost je nutno chápat relativně. Absolutně nepropustná hornina neexistuje. Skutečnost, že kolektor je významně propustnější než izolátor, nevypovídá nic o absolutních hodnotách jejich hydraulické vodivosti. Hodnota hydraulické vodivosti při jinak srovnatelných hydraulických podmínkách ovlivňuje rychlost pohybu vody. Proto při řešení praktických problémů často hodnotíme propustnost či nepropustnost z hlediska časové dimenze řešené úlohy. Izolátor brání volnému přetoku z jednoho kolektoru do druhého. V důsledku toho je v jednotlivých zvodních rozdílný tlak, jehož projevem je piezometrická úroveň ve zvodních s napjatou hladinou nebo poloha hladiny ve zvodních s hladinou volnou. Vzhledem k tomu, že izolační schopnost není absolutní, dochází při takto vzniklém hydraulickém spádu k mezikolektorovému přetoku. Ten je z lokálního hlediska většinou zanedbatelný, v rozsáhlých areálech však může docházet k přetékání významného množství vody, zejména je-li izolátor narušen tektonickou nebo vulkanickou činností. V takovém případě mluvíme o poloizolátoru. Mezikolektorový přetok má podle směru přestupu na proudění obdobný vliv jako infiltrace a vývěry. Podobně působí umělá dotace nebo odběry vody z kolektoru. Jak bylo již zmíněno, kolektory lze dělit podle charakteru hladiny. Kolektor s napjatou hladinou je zdola i shora omezen izolátory. Již z toho je patrné, že se většinou jedná o hlubší kolektory. Kontura zvodně se přizpůsobuje ohraničujícím plochám a všechny volné prostory jsou zaplněny vodou. Proudění v takovém kolektoru označujeme jako nasycené (saturované). Při provrtání nepropustného stropu vystoupí hladina na jistou úroveň, označovanou jako piezometrická (výtlačná) výška. Jestliže voda vystoupí nad úroveň terénu, jedná se o artézskou zvodeň. Zjištění piezometrické


71

výšky v dostatečném počtu bodů umožní konstrukci piezometrické plochy, charakterizující tlakové poměry ve zvodni. V dalším textu se omezíme na tento případ, který vede na lineární diferenciální rovnici a je popsán jednodušším způsobem než obecnější úloha s volnou hladinou. Kolektor s volnou hladinou je situován zpravidla bezprostředně pod zemským povrchem. V aridních oblastech nebo v lokalitách silně ovlivněných těžební nebo jinou průmyslovou činností se však může jednat i o hlubší kolektor. Zvodeň takového kolektoru je shora ohraničena volnou hladinou. Nad ní se nachází nenasycená (nesaturovaná) zóna, kde je volný prostor zčásti naplněn plynnou fází - nejčastěji vzduchem a vodními parami. Bezprostředně nad hladinou podzemní vody se nachází kapilární zóna. Její výška je závislá na vlastnostech horninového prostředí. Čím jsou póry jemnější, tím je větší výška kapilárního vzlínání. Orientační hodnoty jsou uvedeny v tabulce č. 3 (Štamberg 1996). V zóně kapilární vody klesá nasycení se vzdáleností od volné hladiny. Spodní část této zóny je prakticky nasycená. Podpovrchový kolektor je přímo dotován vodou z dešťových a sněhových srážek. Jeho hydrologický režim souvisí s povrchovými vodními toky a nádržemi. Těsně pod povrchem se nachází zóna půdní vody. V ní dochází k vertikálnímu pohybu vody směrem dolů při infiltraci a směrem vzhůru při odpařování (evaporaci) a životních procesech rostlin (transpiraci). Není-li hladina podzemní vody blízko pod povrchem, pak se mezi zónami půdní a kapilární vody nachází přechodné pásmo, obsahující vlhkost vázanou hygroskopickými a kapilárními silami. Tímto pásmem prochází občas infiltrující gravitační srážková voda. Skutečné rozložení vody v podpovrchových partiích je značně komplikované, pro naše úlohy však má zanedbatelný význam. Složitost řešení úloh v kolektorech s volnou hladinou spočívá v tom, že poloha volné hladiny a jednotlivých zón není předem známa, může se s časem měnit a jde tedy o další neznámou. Navíc vzniklé rovnice jsou nelineární a jejich numerické řešení vede k mnoha obtížím.


Druh zeminy písky jemné písky hlinité písky sprašové hlíny hlíny jílovité zeminy jíly

Kapilární výška, m 0,03 - 0,10 0,10 - 0,50 0,50 - 2 2-5 5 - 15 15 - 50 až přes 50

Tab. 4.3: Kapilární výška pro několik druhů zemin

72

Kapitola 5

NELINEÁRNÍ A VÍCEROZMĚRNÉ TRANSPORTNÍ ÚLOHY Dosud jsme se zabývali případem konvekčně-difuzní rovnice, kde neznámá (transportovaná) veličina byla skalární a rovnice byla lineární (koeficienty rovnice byly nezávislé na neznámé, přičemž nehomogenita byla přípustná). V této kapitole ukážeme příklady obecnějších úloh, z hlediska obou zmíněných vlastností: transport vícerozměrné veličiny (v první kapitole zmíněný zákon zachování hybnosti – Eulerovy a Navier-Stokesovy rovnice) a úlohy s nelineárním konvekčním členem (ukážeme na 1D Eulerově rovnici a jednu nefyzikální aplikaci: pohyb aut na silnici).

Příklad nelineární rovnice bude doplněno Úvodní kapitola: formulace transportu pro obecnou závislost toku na rozložení veličiny: ∂XV + ∇ · ΦX − rX = 0 (5.1) ∂t V matematické formě v 1D (neznámá fce u(x, t), tok veličiny u označíme f (u):) ut + fx = r

(5.2)

Pro konvekci je f = vu, pro difuzi f = ux = ∂u , nyní uvažujeme obecnou závislost f (u) ∂x (algebraickou, ne diferenciální). Tvar (5.2) se nazývá konzervativní. Lze jej upravit na nekonzervativní tvar provedením derivace složené funkce f (u): ux + v˜(u)ux = r

(5.3)

kde v˜(u) = ∂f∂u(u) . Jak uvidíme, v˜ hraje podobnou roli jako v lineární rovnici konvekce rychlost, s tím rozdílem, že v nelineárním případě je „rychlost vlnyÿ funkcí transportované veličiny, což vede na kvalitativně odlišná řešení a speciální jevy (např. vznik rázových vln).

Kapitola 5. Nelineární a vícerozměrné transportní úlohy

5.1

74

Eulerovy a Navier-Stokesovy rovnice v konzervativním tvaru

Eulerovy a Navier-Stokesovy rovnice jsou zároveň příkladem nelineárních hyperbolických a parabolických rovnic (tj. konvekčních a konvenčně-difuzních) a zároveň příkladem systému rovnic (tedy zobecnění dosud uvažovaných lineárních úloh popsaných jednou rovnicí se skalární proměnnou). Navier-Stokesovy rovnice jsou vyjádřením zákona zachování hybnosti (pohybovou rovnicí) pro viskózní tekutinu (obecně stlačitelnou), tj. kapalinu nebo plyn. Při interpretaci jako transportní proces jde tedy o konvekčně-difuzní transport hybnosti. Pro úplný popis pohybu stlačitelné tekutiny je třeba uvažovat též rovnice bilance (transportu) hmoty a energie (tepla) a stavovou rovnicí (vztah hustota – tlak). Eulerovy rovnice jsou pak speciální případ Navier-Stokesových rovnic pro neviskózní tekutinu. Obvyklý zápis Navier-Stokesových rovnic používaný v mechanice tekutin je 1 ∂v + (v · ∇)v = f − ∇p + ν∆v ∂t % kde

(5.4)

 v − rychlost  % − hustota neznámé  p − tlak

 ν − kinematická viskozita  f − objemová síla vztažená na jednotku hmotnosti tekutiny parametry  (tj. např. tíhové zrychlení)

V tomto případě jde o tři nezávislé rovnice uvažujeme-li proudění v 3D prostoru. Ukážeme nyní formulaci rovnic v tzv. konzervativním tvaru, která má tvar transportní rovnice (5.1), resp. (5.2). Jde v tomto případě o soustavu transportních rovnic (transport hmoty, energie a jednotlivých složek hybnosti), tj. 5 rovnic pro 3D případ. Tu lze vyjádřit ve tvaru „vektorovéÿ transportní rovnice, kde neznámá veličina je sada skalárních veličin (algebraický vektor bez fyzikálního významu vektoru) a počet složek toku je součin počtu transportovaných skalárních veličin a dimenze prostoru (lze si představit jako matici). V jednorozměrném případě (jedna složka rychlosti a hybnosti) je soustava NavierStokesovy pohybové rovnice, rovnice kontinuity a rovnice energie zapsaná v konzervativním tvaru následující: ∂Q ∂E + =0 (5.5) ∂t ∂x kde   % Q =  %v  (5.6) e


75

je vektor neznámých transportovaných veličin (hmota, hybnost, energie), difúzní člen konvekční člen z }| z }| { {  0 %v 4 ∂v  µ E =  %v 2 + p  −  3 ∂x ∂v v(e + p) − κ ∂T − 43 µv ∂x ∂x

(5.7)

je vektor toků zmíněných veličin a hustota energie (vztažená na jednotkový objem) jako součet vnitřní energie a kinetické energie e=

p 1 + %v 2 γ−1 2

(5.8)

Neznámé veličiny jsou tedy hustota %, rychlost (skalár) v a hustota energie e (tedy tři), další neznámá tj. tlak je určena rovnicí (5.8) pomocí hustoty energie a rychlosti, podobně teplota T je určena vnitřní energií. Parametry jsou µ = %ν dynamická viskozita a κ tepelná vodivost.

5.2

Převod systému rovnic na „kanonickýÿ tvar

Ukážeme obecný postup převedení systému hyperbolických rovnic prvního řádu (tj. rovnice transportu vektorové veličiny) na řešení jednotlivých nezávislých rovnic (tj. transport skalární veličiny), tj. převedení obecnější úlohy na úlohu, kterou umíme řešit (nebo alespoň názorně fyzikálně interpretovat). Nejprve pro jednoduchost uvažujeme lineární případ. Systém hyperbolických rovnic v nekonzervativním tvaru má tvar   u1 (x, t) !  u2 (x, t)  a11 . . .   Ut + AUx = 0 U(x, t) =  (5.9)  A= .. ..   . . aN N uN (x, t)

kde ui (x, t) jsou neznámé funkce (N skalárních funkcí) a aij jsou konstanty (nejsou funkcí ui) a derivace jsou definovány po složkách. Po složkách jde o soustavu rovnic N

∂u1 X ∂ui + a1i = 0 ∂t ∂x i=1 .. . ∂uN + ∂t

N X i=1

aN i

∂ui = 0 ∂x

(5.10)


76

Transformace soustavy (vázaných) rovnic na soustavu nezávislých skalárních rovnic je založena na algebraických operacích s maticí A. Pro hyperbolický systém lze matici A diagonalizovat (tato vlastnost se obvykle používá k definici hyperboličnosti systému): X −1 AX = Λ

tj. A = XΛX −1

kde Λ obsahuje na diagonále vlastní čísla  λ1  .. Λ= .

(5.11)

 λN

 

která jsou pro hyperbolický systém reálná (to je diagonalizovatelnost) a matice X obsahuje ve sloupcích vlastní vektory. Následující algebraické úpravy původního systému ∂U ∂U · X−1 zleva +A =0 ∂t ∂x

∂(X−1 U) ∂U + X−1 (X · Λ · X−1 ) =0 ∂t ∂x ∂(X−1 U) ∂(X−1 U) +Λ =0 ∂t ∂x umožní zavedením nového (transformovaného) vektoru neznámých W = X −1 U převést původní rovnici na tvar Wt + ΛWx

tj.

∂wi ∂wi + λi = 0 i = 1...N ∂t ∂x

(5.12)

tedy systém N nezávislých skalárních rovnic konvekce. Vlastní čísla matice A mají význam rychlosti jakou jsou transportovány jednotlivé nové veličiny wi (které jsou lineární kombinací původních veličin ui). Výsledné řešení je tedy superpozice „vlnÿ postupujících obecně různou rychlostí různým směrem. Pro nelineární systém lze transformací rovněž provést a potom i vlastní čísla (rychlosti) jsou funkcí neznámé veličiny U resp. W podobně jako ve skalárním případě (rovnice (5.3)). Např. pro 1D Eulerovy rovnice q jsou u, u − c a u + c, kde u je rychlosti proudění (jedna z neznámých veličin) a c = γp je rychlost zvuku. Vidíme tedy jak je %

podstatné rozlišení podzvukového a nadzvukového proudění (rychlosti mají buď stejné nebo různá znaménka) - to má význam např. pro korektní zadání okrajových podmínek (vtoková vs. výtoková hranice) a pro numerické řešení (konstrukce metody typu upwind).


5.3

77

Nalezení řešení rovnice pomocí charakteristik

Pojem charakteristiky jsme zmínili v kapitole o numerických metodách při fyzikální interpretaci CFL podmínky stability upwind explicitního diferenčního schématu pro rovnici konvekce. V případě úlohy v jednorozměrném prostoru (tj. rovnice s neznámou funkcí u(x, t)), jde o křivky v rovině (x, t), na nichž je hodnota řešení u konstantní. Pro lineární hyperbolickou rovnici prvního řádu ut + vux = 0 (tj. konvektivní transport rychlostí v) jsou charakteristiky přímky o rovnici x − vt = 0, což přesně koresponduje s tím, že řešení rovnice s počáteční podmínkou u0 (x) je u(x, t) = u0 (x − vt). Nyní ukážeme význam charakteristik u nelineární skalární hyperbolické rovnice ut + v˜(u)ux = 0. Charakteristiku lze definovat zcela analogicky jako přímku se směrnicí 1/˜ v(u), přičemž je zjevné, že její směr závisí na hodnotě řešení a tedy obecně v každém bodě jiný. Uvažujeme-li počáteční podmínku u(x, 0), pak z každého bodu x0 vychází charakteristika o rovnici x = x0 + v˜[u(x0 , 0)] t (5.13) na níž je hodnota řešení u(x, t) rovna odpovídající počáteční podmínce u(x0 , 0). Máme tedy postup jak najít přesné řešení nelineární rovnice, za podmínky, že lze jednoznačně určit charakteristiku, která prochází bodem počáteční podmínky (x0 , 0) a bodem hledaného řešení (x, t). Tato podmínka není splněna např. pokud se charakteristiky protínají.

5.3.1

Rázová vlna

Vysvětlíme na příkladu fyzikální smysl případu, kdy se protínají charakteristiky, tj. vznik rázové vlny. Jde o důsledek nelinearity úlohy: řešením úlohy se spojitými parametry je vznikne v určité časové hladině nespojitá funkce. Je nutno poznamenat, že v tom případě pak nejde o řešení v klasickém smyslu (nespojitou funkci nelze „dosaditÿ do diferenciální rovnice), je nutno definovat v jistém smyslu obecnější pojem řešení. Tímto problémem se zde zabývat nebudeme, omezíme se na fyzikální znázornění vzniku nespojitosti (rázové vlny). Uvažujeme rovnici ut + uux = 0

resp. ut + fx = 0 , f (u) =

u2 2

(5.14)

což je tzv. Burgersova rovnice, která je zjednodušením Eulerových rovnic proudění tekutiny (pro 1D), neznámá u(x, t) má fyzikální význam rychlosti a budeme hledat její řešení pro x ∈ R a t > 0 s počáteční podmínkou danou následovně (obrázek 5.1):  x<0  1 1 − x x ∈ (0, 1) u0 (x) = (5.15)  0 x>1


78

u

0

x

1

t

rázová vlna t=1

(1,1)

ch1: t = x-x0 ch3: x=x0 x0

0

x0

1

x

ch2: t =( x-x0)/(1-x0)

Obrázek 5.1: Počáteční podmínka (5.15) pro úlohu s Burgersovou rovnicí. Charakteristiky a vznik rázové vlny pro tuto úlohu. Rychlost posunu vlny, tj. převrácená hodnota směrnice charakteristik, je dána koeficientem v rovnici v˜(u) = u. Podle vztahů pro počáteční podmínku tedy můžeme rozlišit tři skupiny charakteristik: (1) x − 1 · t = k (2) x − (1 − x)t = k (3) x − 0 · t = k

které jsou vykresleny na obrázku 5.1. Pro časy 0 < t < 1 tedy můžeme snadno určit řešení rovnice: hodnoty u = 0 se nepohybují (pro x > 1) zůstává řešení nulové, hodnoty u = 1 se pohybují rychlostí 1 doprava a přechod se postupně stává strmější (obrázek 5.2), až v čase t = 1 dojde k vzniku nespojitosti – rázové vlny (na hranici mezi u = 1 a u = 0 v bodě x = 1, v čase t = 1). To je bod v rovině (x, t), kde se protínají charakteristiky vycházející z intervalu < 0, 1 > počáteční podmínky. Na zjištění „co se bude dít dálÿ nám dosavadní výklad nestačí: za průsečíkem by původní charakteristiky ztratily fyzikální smysl. Výsledek uvedeme bez odvození: Z místa průsečíku bude pokračovat jedna charakteristika, která reprezentuje pohyb rázové vlny. Rychlost pohybu (tj. sklon charakteristiky) je dána vztahem rychlost =

f (uR ) − f (uL) uR − uL

(5.16)

kde uR , resp. uL je hodnota řešení zprava resp. zleva rázové vlny a f (u) je funkce toku

Kapitola 5. Nelineární a vícerozměrné transportní úlohy u

u

t=1/3

0

x

1

u

1

t=2/3

0

x

x

1

u

t=1

0

79

t>1

0

1

x

Obrázek 5.2: Časový vývoj rozložení funkce u před a po vzniku rázové vlny. v řešené rovnici. V případě naší úlohy je f (u) = rychlost =

u2 , 2

uR = 0, uL = 1, tj.

f (0) − f (1) 1 = 0−1 2

(5.17)

Příslušná charakteristika je rovněž zakreslena na obrázku 5.1, průběh řešení v prostoru a v čase na obrázku 5.2.

5.4

Model dopravního proudu

Uvažujeme úlohu pohybu aut na silnici, což lze popsat analogickou technikou jako např. proudění tekutiny. Hlavní rozdíl pak bude ve vyjádření závislosti toku na hustotě, což bude specifikum tohoto procesu ve srovnání s jinými úlohami toku. Pro popis děje budeme používat tyto veličiny: počet aut u – hustota provozu jednot. délka silnice počet aut f – tok daným místem silnice (příčným řezem) čas v – rychlost jízdy aut, f = v · u

Rovnici pohybu aut (tj. např. závislost rozložení aut na silnici v prostoru a čase) dostaneme uplatněním dvou výchozích principů: 1. bilance (zákon zachování) a 2. závislost rychlosti aut (resp. toku) na hustotě provozu. První z nich je popsán univerzální rovnicí (viz úvod kapitoly) ut + fx = 0 Závislost rychlosti na hustotě není dána jednoznačně, je třeba sestavit empirický vztah, vyjadřující ve zvolené míře zjednodušení různé vlivy: technické možnosti jízdy kolony aut po silnici, psychologie chování řidiče apod. Obecně je třeba takovýto vztah (závislost) určit měřením. My použijeme tuto základní úvahu: rychlost pohybu aut je maximální při nulové hustotě (samotné jedno auto není omezováno) a nulová při maximální

Kapitola 5. Nelineární a vícerozměrné transportní úlohy f

v

v

umax

umax

u

u 0

80

0

0

umax

u

Obrázek 5.3: Empirické závislosti rychlosti pohybu aut v, toku aut f a rychlosti „posunu vlnyÿ v˜ na hustotě provozu. hustotě (dojde k zácpě a auta nemohou jet). Předpokládáme, že maximální rychlost je stejná pro všechny auta. Potom u v(u) = A(k − u) = vmax 1 − umax Tok aut je pak

f (u) = u · v(u) = u · vmax 1 −

u

umax Ze vztahu pro tok můžeme jednoduše určit rychlost pohybu vlny (pozor, je třeba odlišit od rychlosti jízdy aut; toto je rychlost jakou se posouvá místo s určitou hustotou provozu) ∂f 2u rychlost vlny: v˜ = = vmax 1 − ∂u umax Všechny tři závislosti (odvozené z našeho vlastního empirického vztahu v(u)) jsou na obrázku 5.3. Příklady: Ukážeme řešení rovnice dopravního proudu pro různé počáteční podmínky. Pro jednoduchost uvažujeme tyto hodnoty parametrů: vmax = 1

umax = 1

Rychlost vlny určující sklon charakteristik pak je v˜ = 1 − 2u

Varianta 1 Uvažujeme počáteční podmínku uvedenou v úloze pro Burgersovu rovnici – rovnice (5.15). Průběh počáteční podmínky, výsledné charakteristiky a schéma pohybu vlny jsou nakresleny v obr. 5.4. Charakteristiky pro x < 0 mají směrnici v˜−1 = (1 − 2 · 1)−1 , pro x > 1 mají směrnici v˜−1 = (1 − 2 · 0)−1 a v intervalu (0, 1) se převrácená hodnota směrnice lineárně mění. Pro x < 21 tedy dochází k pohybu vlny v opačném směru než je rychlost jízdy aut. V tomto případě tedy existuje spojité řešení pro všechny časy t > 0, charakteristiky se neprotínají. Řešení odpovídá přirozené představě, že dochází k postupnému rozjíždění aut a zřeďování přechodové vlny (výjezd aut z parkoviště, rozjíždění kolony).


81

u

0

t

x

1

x - t = konst.

x + t = konst.

charakteristiky x

0 1/2 1

smìrnice= -1 u -v max

smìrnice= +1

v max 0

1

x

Obrázek 5.4: Úloha dopravního proudu, varianta počáteční podmínky č.1 – počáteční podmínka, charakteristiky a vývoj rozložení hustoty v čase.


82

u 1 = u max

x

t

rázová vlna

t=1/2

0 u

-1

x

1

1

-1/2

u

1/2

1

3/2

x

x

Obrázek 5.5: Úloha dopravního proudu, varianta počáteční podmínky č.2 – počáteční podmínka, charakteristiky s vyznačením vzniku rázové vlny, tvar rozložení hustoty v okamžiku vzniku rázové vlny a charakteristika dalšího děje.


83 t

u u max

1/3 -1

1

x

-1

x x-1/3t = 0

Obrázek 5.6: Úloha dopravního proudu, varianta počáteční podmínky č.3 – počáteční podmínka, charakteristiky s vyznačením vzniku rázové vlny, tvar rozložení hustoty v okamžiku vzniku rázové vlny a charakteristika dalšího děje.

Varianta 2 Tato varianta je uvedena na obrázku 5.5, tj. graficky znázorněná počáteční podmínka a vypočtené charakteristiky. V tomto případě dojde ke vzniku rázové vlny (nespojitosti) v čase t = 1/2. Charakteristiky jsou nakresleny na obrázku 5.5. Rychlost pohybu rázové vlny je vR =

f (uR ) − f (uL ) 0−0 = =0 uR − uL 1−0

což koresponduje s představou: hranice stojích aut a žádných aut se nepohybuje. Tento stav však je pouze hypotetický v prvním okamžiku vzniku rázové vlny, neboť v dalších časových okamžicích dochází k protnutí jiných charakteristik a snížení maxima hustoty aut pod hodnotu umax = 1. Popis děje je pak mnohem složitější. Nespojitost (rozhraní nulové hustoty s nenulovou, ale ne maximální hustotou) se pak začne pohybovat, náznak průběhu je na obr. 5.5 dole.

Varianta 3 Tento případ je obměnou varianty 2 pro jinou hodnotu hustoty ( 31 místo 1) v bodě x = 0. Tato hodnota byla zvolena cíleně tak, aby byla menší než 12 a tedy všechny charakteristiky měly kladnou směrnici. Charakteristiky z bodů -1 a 0 jsou (obrázek 5.6) x − t = −1 a jejich průsečík je

1 x− t=0 3

3 1 2 t=1⇒t= , x= 3 2 2


84

I v tomto případě tedy dojde ke vzniku rázové vlny a její rychlost je vR =

f (0) − f (1/3) 2 = 0 − 1/3 3

Děj po vzniku rázové vlny bude obdobný jako u varianty 2 s tím rozdílem, že v okamžiku vzniku rázové vlny nebude její rychlost nulová a nepůjde o stav s nepohybujícími se auty (na obrázku 5.5 dole v podstatě některá z pozdějších fází). Opět připomínáme rozdíl mezi rychlostí aut a rychlostí vlny: rázová vlna nemusí nijak souviset se „zácpouÿ na silnici nebo „náhlýmÿ rozjezdem/zastavením aut, jde o nespojitost v hustotě provozu, v konkrétním případě variant 2 a 3 tedy konec souvisle jedoucí kolony aut za kterým následuje prázdná silnice.

5.5 5.5.1

Příklady úpravy systému rovnic na kanonický tvar Vlnová rovnice

Vlnová rovnice

2 ∂2u 2∂ u = a ∂t2 ∂x2 je základní příklad hyperbolické parciální rovnice druhého řádu a její význam by měl být znám ze základního kurzu fyziky. V této podobě popisuje např. příčnou výchylku u pružné struny v jednotlivých bodech v závislosti na čase. Rovnici lze vyjádřit též jako systém hyperbolických rovnic 1. řádu (derivováním a sečtením těchto dvou rovnic dostaneme původní vlnovou rovnici)

∂v ∂u =a ∂t ∂x ∂v ∂u =a ∂t ∂x Tento systém rovnic lze zapsat ve smyslu sekce 5.2 ve vektorovém tvaru ∂v ∂ v 0 −a ∂x + · ∂u =0 −a 0 ∂t u ∂x

tj. Ut + AUx = 0 kde

U=

v u

v − rychlost u − výchylka

A=

0 −a −a 0

Nyní můžeme aplikovat techniku popsanou v sekci 5.2, tj. transformujeme rovnici na soustavu dvou nezávislých rovnic. Transformace spočívá v nalezení vlastních čísel matice A a vyjádření nových proměnných w1 , w2 pomocí vlastních vektorů A a původních proměnných u a v.


85

Řešením rovnice det (A − λI) = 0 tj.

−λ −a −a −λ

= 0 tj. λ2 − a2 = 0

dostáváme dvě vlastní čísla (reálná) λ1,2 = ±a. K nim příslušné vlastní vektory určíme řešením rovnic (A − λ1,2 I)x = 0 postupně pro obě vlastní čísla −a −a x11 0 1 +a : = ... −a −a x12 0 −1 a −a w21 1 −a : =0 −a a w22 1 Pro sestavení transformační matice X použijeme jednotkové vlastní vektory, tj. 1 −1 1 X=√ 1 1 2

Platnost transformace A = X −1 ΛX můžeme snadno ověřit ručním pronásobením matic 1 1 λ1 −1 1 a 0 1 1 −1 √ = A=X X=√ λ2 1 1 0 −a −1 1 2 2 1 0 −2a −1 1 a a 1 = = 0 1 1 a −a 2 2 −2a Nové proměnné v transformované soustavě pak jsou u v w1 u =X tj. w1 = − √ + √ w2 v 2 2

u v w2 = √ + √ 2 2

Výsledná soustava dvou nezávislých rovnic (z pohledu proměnných w1 a w2 ) je ∂(−u + v) ∂(−u + v) +a = 0 ∂t ∂x ∂(u + v) ∂(u + v) −a = 0 ∂t ∂x Kontrolním sečtením a odečtením se dostaneme k původní soustavě rovnic pro neznámé u a v: ⊕: 2vt − 2aux = 0 : − 2ut + 2avx = 0

původní systém vt = aux ut = avx

Interpretace separované soustavy rovnic jako transportní proces tedy je, že rovnice popisuje transport veličin w1,2 rychlostmi ±a, což je v souladu s fyzikálním procesem vlnění: šíření vlny rychlostí a v kterémkoli směru.


86

h(x,y,t)

z(x,y)

Obrázek 5.7: Schéma popisu proudění mělké vody – neznámá tloušťka vodní vrstvy h a vstupní funkce výšek terénu z.

5.5.2

Rovnice mělké vody

Takzvaná „rovnice mělké vodyÿ (shallow water equation) je rovnicí pohybu neviskózní kapaliny po nevodorovném zemském povrchu. Je vyjádřením Eulerových rovnic pro tenkou vrstvu kapaliny, konkrétně snížením dimenze z 3D na 2D pomocí integrace ve svislém směru. Jako vnější síla se pak uplatňuje gravitace. Rovnice tedy v praxi popisuje proudění vody v řece a různých umělých vodohospodářských zařízeních. Numerické modely založené na jejím řešení se uplatní např. pro optimalizaci manipulace s přehradami nebo pro předpovědi průběhu a následků povodní. Jde o rovnice ve 2D (průmět zemského povrchu) a neznámými funkcemi jsou tloušťka vodní vrstvy h(x, y, t) (obr. 5.7) a vodorovná složka rychlosti proudění vody (předpokládá se ve všech místech podél svislice stejná) v(x, y, t). Rovnice v konzervativním tvaru je ∂G ∂Q ∂F + + =S (5.18) ∂t ∂x ∂y kde 

       hu hv h 0 2  S =  gh(s0x − sf x )  huv Q =  hu  F =  hu2 + gh2  G =  gh2 2 hv gh(s0y − sf y ) huv uv + 2

přičemž vstupní parametry rovnice jsou

• sf x , sf y – tření o zemský povrch vyjádřené jako ekvivalentní snížení sklonu terénu s0x • s0 = = ∇z(x, y) – vektor sklonu zemského povrchu odvozený od funkce s0y výšek terénu (tvar povrchu) z(x, y) (obr. 5.7) Separaci rovnic si ukážeme na zjednodušeném příkladu pro 1D oblast (kanál) s vodorovným povrchem. Soustava rovnic je ∂h ∂u ∂h +u +h =0 ∂t ∂x ∂x


87

∂u ∂h ∂u +u +g =0 ∂t ∂x ∂x pro neznámé skalární funkce h(x, t) (tloušťka vodní vrstvy) a u(x, t) (rychlost). V maticové formě pak h u h U= A= u g u Vlastní čísla matice určíme velmi snadno řešením rovnice dle definice (u − λ)2 − gh = 0 p λ1,2 = u ± gh

Připomínáme, že rovnice jsou nelineární a tedy i příslušná vlastní čísla jsou závislá na neznámých h, u.

5.6

Konvekčně-reakčně-difuzní rovnice

V této sekci navážeme na výklad procesů v porézním prostředí, konkrétně transportu látky s nelineární sorpcí (tj. reakční člen speciálního tvaru). Problematika je však mnohem obecnější, pojem reakčně-difuzních rovnic se uplatňuje v mnoha jiných oborech a vede na složité matematické problémy, jde o popis kvalitativně odlišných jevů (deterministický chaos). Rovnice transportu v upraveném tvaru je R

∂c + v · ∇c − ∇ · (D h ∇c) = 0 , ∂t

R=1+

1 − n ∂s · n ∂c

(5.19)

kde R je retardační faktor a s = f (c) obecná funkce popisující závislost sorbované koncentrace na koncentraci v roztoku – „příčinaÿ nelinearity rovnice. Retardační faktor R > 1 vyjadřuje zpomalení transportu látky vlivem sorpce oproti procesu bez vlivu sorpce R = 1. Dělením R lze rovnici upravit na tvar ∂c ˜ h (c) ∇c) = 0 + v˜(c) · ∇c − ∇ · (D ∂t kde

(5.20)

v ˜ h (c) = Dh D R(c) R(c) který až na difuzní člen odpovídá svojí strukturou rovnici (5.3), jejíž řešení pro konkrétní příklady nelineární funkce v˜(c) jsme si ukázali v sekcích 5.3.1 a 5.4. Obdobně lze tedy určit chování rozpuštěné látky vlivem konvekce a nelineární sorpce – tj. různé koncentrace látky se budou posouvat různou rychlostí v závislosti na tvaru funkce s(c). Přiblížíme stručně charakter řešení pro konkrétní funkce nelineární sorpce – Freundlichovu a Langmuirovu izotermu (viz kapitola Porézní prostředí). v˜ (c) =


88

Freundlichova izoterma Freundlichova izoterma je definována závislostí s = kF caF , kde kF a aF jsou konstanty. Pro aF = 1 jde o lineární izotermu a hodnota R je pak nezávislá na c a rovnice není nelineární – retardace je stejná pro všechny hodnoty koncentrace a řešení je tedy analogické situaci bez sorpce. Pro aF > 1 je derivace podle c rostoucí a tedy i R(c) je rostoucí funkcí c a v˜(c) je klesající funkcí c. To znamená, že transport látky ve vyšší koncentraci bude pomalejší než transport látky v nižší koncentraci. Situaci si lze představit na 1D úloze transportu s počáteční podmínkou ve tvaru gaussovského rozdělení exp(−x2 ): přechod mezi malou a velkou koncentrací bude pozvolnější na straně po směru rychlosti a strmější na straně proti směru rychlosti – pohyb vrcholu „kopečkuÿ je oproti úpatí zbrzděn. Pro reálnější představu je třeba vzít v úvahu i difuzi, tedy strmost přechodu se bude snižovat, sice na obou stranách symetricky, ale různě pro různé hodnoty koncentrace. Pro aF < 1 je derivace podle c klesající a tedy i R(c) je klesající funkcí c a v˜(c) je rostoucí funkcí c. To znamená, že transport látky ve vyšší koncentraci bude rychlejší než transport látky v nižší koncentraci. Pro výše uvedenou gaussovskou počáteční podmínku bude přechod mezi malou a velkou koncentrací pozvolnější na straně proti směru rychlosti a strmější na straně po směru rychlosti. Platí stejná úvaha pro difuzi. Podstatný rozdíl však je, že pro c = 0 je derivace nekonečná a v˜(c = 0) = 0, což znamená, že příslušný konvekční člen z rovnice zmizí (rychlost transportu se sníží na nulu - s látkou se pokud není difuze neděje nic), rovnice je tzv. degenerovaná. Tato vlastnost přináší významné obtíže při numerickém řešení rovnice i při teoretickém studiu jejích vlastností. Při malé difuzi vede řešení na vznik nespojitosti koncentrace (velmi strmý náběh koncentrace na přední straně „vlnyÿ).

Langmuirova izoterma Langmuirova izoterma je definována závislostí s=

kL smax c 1 + kL c

kde kL a smax jsou konstanty, přičemž smax má zřejmý fyzikální význam sorbovaného množství při nasycení (c → ∞). V blízkosti c = 0 se závislost blíží lineární, derivace a tedy i retardační faktor jsou pro tyto hodnoty největší, ale narozdíl od druhého případu Freundlichovy izotermy jde o konečné hodnoty a nedochází tedy k degeneraci rovnice. Naopak pro velké koncentrace je chování systému ovlivněno sorpcí jen velmi málo (malá retardace, v˜(c) ≈ v).

LITERATURA [1] Bear, Verruijt: Modeling Underground Flow and Pollution, Reidel Pbl., 1990. [2] W.J. Beek, K.M.K. Muttzall, J.W. van Heuven Transport Phenomena, Wiley&sons, 1999. [3] M. Brdička, L. Samek, B. Sopko: Mechanika kontinua, Academia, Praha, 2000. [4] Císlerová, Vogel: Transportní procesy, skripta ČVUT (FSv), 1998. [5] S.J. Farlow: Partial differential equations for scientists and engineers, John Wiley & Sons, 1982, ruský překlad Moskva Mir 1985. [6] Hirsch,C.: Numerical Computation of Internal and External Flows, Volume 1 Fundamentals of Numerical Discretization, John Wiley & Sons Ltd., 1991. [7] Pulliam T.H., Lomax, H.: Fundamentals of Computational Fluid Dynamics, Springer, 2001. [8] Štamberk: Modelování migračních procesů v životním prostředí, skripta ČVUT (FJFI), 1998. [9] Vitásek, E.: Numerické metody (Technický průvodce 67), SNTL, 1987. [10] Vitásek, E.: Základy teorie numerických metod pro řešení diferenciálních rovnic, Academia, 1994. [11] Zheng C. and Bennett G.H.: Applied contaminant transport modeling, Van Nostrand Reinhold, New York, 1995.

Transportní procesy. Učební text. 23. září Milan Hokr. Fakulta mechatroniky a mezioborových inženýrských studií Technická univerzita v Liberci

Recommend Documents