1
STABILITA A DYNAMIKA
PŘÍRODNÍCH SYSTÉMŮ
JOSEF ZEMAN, PROSINEC 2009
2
OBSAH
1 Úvod ......................................................................................................................................................................4
2 Dynamika procesů .................................................................................................................................................5
2.1 Koncepce rezervoárů ...................................................................................................................................6
2.2 Jednorezervoárový systém...........................................................................................................................9
2.2.1 Časový vývoj.............................................................................................................................................9
2.2.2 Stacionární stav......................................................................................................................................11
2.2.3 Doba zdržení ..........................................................................................................................................12
2.2.4 Doba odezvy ..........................................................................................................................................12
2.3 Dvourezervoárový systém..........................................................................................................................13
2.3.1 Analytické řešení....................................................................................................................................14
2.3.2 Doba zdržení ..........................................................................................................................................16
2.3.3 Doba odezvy ..........................................................................................................................................17
2.4 Třírezervoárový systém..............................................................................................................................18
2.5 Vícerezervoárové systémy ....................................................................... Chyba! Záložka není definována.
2.6 Nelineární systémy................................................................................... Chyba! Záložka není definována.
3 Příklady ................................................................................................................................................................21
3.1 Příklad 1......................................................................................................................................................21
3.1.1 Zadání: ...................................................................................................................................................21
3.1.2 Řešení ....................................................................................................................................................21
3.2 Příklad 2......................................................................................................................................................25
3.2.1 Zadání ....................................................................................................................................................25
3.2.2 Řešení ....................................................................................................................................................25
3.3 Příklad 3......................................................................................................................................................29
3.3.1 Zadání ....................................................................................................................................................29
3.3.2 Řešení ....................................................................................................................................................29
3
4 Vztah mezi stacionárním stavem a termodynamickou rovnováhou ...................................................................36
4.1 Arrheniova rovnice.....................................................................................................................................37
4.2 Termodynamické odbočení........................................................................................................................37
4.3 Dynamika procesů......................................................................................................................................39
4.4 Gibbsova aktivační funkce..........................................................................................................................41
4.5 Další interpretace.......................................................................................................................................44
4.6 Vztah mezi rychlostí rozpouštění a srážení ................................................................................................45
5 Numerické řešení diferenciálních rovnic .............................................................................................................47
5.1 Využití numerického řešení diferenciálních rovnic v geochemii................................................................50
5.2 Příklad 1......................................................................................................................................................50
5.2.1 Zadání ....................................................................................................................................................50
5.2.2 Řešení ....................................................................................................................................................50
5.3 Příklad 2......................................................................................................................................................59
5.3.1 Zadání ....................................................................................................................................................59
5.3.2 Řešení ....................................................................................................................................................60
6 Seznamy...............................................................................................................................................................64
6.1 Seznam obrázků .........................................................................................................................................64
6.2 Seznam tabulek..........................................................................................................................................67
4
1 ÚVOD
Země jako celek a její jednotlivé části podléhají po celou dobu jejího geologického vývoje změnám. V zásadě
mohou být její jednotlivé části od makroměřítka až po mikroměřítko považovány za rezervoáry, mezi kterými
dochází k přenosu (toku) hmoty. Takovými rezervoáry jsou například v makroměřítku atmosféra, oceány, zemská
kůra a biota. Samozřejmě je možné měřítko rezervoárů zjemňovat až na úroveň pórů v horninách, které
obsahují horninovou vlhkost. Přestože dochází mezi rezervoáry k výměně hmoty, zůstávají koncentrace či
množství hmoty v rezervoárech dlouhodobě stabilní.
Základní přístup k přírodním systémům začíná jejich termodynamickým popisem. Základní principy (zákony)
termodynamiky říkají, jaké procesy a kterým směrem se v přírodě odehrávají. Podle prvního principu se
v přírodě odehrávají jen takové procesy, při kterých se zachovává veličina, kterou označujeme jako energie.
Druhý princip ukazuje, kterým směrem se procesy ubírají – vždy tím směrem, ve kterém roste celková entropie
systému. Příroda se vždy snaží dosáhnout stavu s maximální možnou celkovou entropií. Ukazatelem na celkovou
změnu entropie systému je Gibbsova funkce. Pokud Gibbsova funkce klesá, pak celková entropie roste a
naopak. Termodynamika má při popisu přírodních systémů jednu výraznou výhodu – nezávisí na cestě, po
které se systém z výchozího do konečného stavu dostal.
Tento fakt je však také hlavní nevýhodou termodynamiky – přestože víme, do kterého stavu systém směřuje,
termodynamika nemůže předpovědět, za jak dlouho a zda vůbec do tohoto stavu dojde. O tom rozhoduje
dynamika procesů, tedy jejich rychlost. Ta může být tak velká, že se přeměny odehrají prakticky okamžitě,
nebo tak nízká, že ani v geologické časové škále k přeměnám, které předvídá termodynamika, nedojde.
Kromě toho často dochází k tomu, že se koncentrace látek v rezervoárech udržují na stabilních hodnotách, ty
však neodpovídají žádné z možných termodynamických rovnováh. Příčinou jsou stacionární stavy, které se
ustavují mezi vstupním a výstupními toky do a z rezervoárů. Právě těmto procesům je věnován následující text.
Pro studium uvedené problematiky je třeba mít základní znalosti termodynamiky, matematických operací a
přírodních systémů. Předkládaný text představuje základní přístup k řešení dynamiky přírodních systémů. Ukazuje,
jakým způsobem je možné ke studiu dynamiky přistupovat, uvádí definice základních parametrů a je doplněn
ukázkovými příklady. Zároveň je konstruován jako „samonosný“ tak, aby při jeho studium nebylo třeba
dalších textů či dalších pomůcek. Nutný matematický aparát, použitý pro numerické řešení příkladů, je uveden
v příloze. Kromě toho je do textu zařazena i problematika, která ukazuje propojení mezi termodynamikou a
kinetikou procesů a tím vytváří jednotný základ pro studium přírodních systémů.
Při studiu přeji hodně napínavých chvil.
Prosinec 2009
Josef Zeman
5
2 DYNAMIKA PROCESŮ
Systém je jakákoliv část Vesmíru, kterou si pozorovatel či badatel zvolil pro své studium. Může být velký nebo
malý, jednoduchý nebo složitý – může se pohybovat od jednotlivých atomů až po celý Vesmír. Může jím být
vzorek horniny, pórová voda, jezero, oceán, celá planeta. Obvykle je studovaný systém součástí většího systému
a tak tvoří jeho podsystém. Pochopení funkce studovaného systému – podsystému většího systému – není
obvykle možné bez jeho zasazení do kontextu jeho vnějších vazeb na větší systémy.
Stav přírodních systémů je spolehlivě určen jeho teplotou, tlakem a složením. Termodynamika je schopna na
základě analýzy stavu systému spolehlivě zjistit, zda je systém ve stabilním stavu nebo ve stavu, který se bude
nebo může změnit. Stabilní stav se za daných podmínek vyznačuje nejvyšší možnou celkovou entropií látek
v něm obsažených. Je to zároveň stav s nejnižší hodnotou Gibbsovy funkce, která je měřítkem celkové entropie.
Takový stav se označuje jako stav termodynamické rovnováhy.
Pokud tomu tak není, nachází se systém v nestabilním stavu a příroda bude mít tendenci tento stav změnit na
stav rovnovážný (stabilní). Hnací silou této změny je, že tímto procesem získá další entropii, která je rovna
rozdílu mezi entropií stávajícího nestabilního stavu a rovnovážného stavu.
To, že je systém v nerovnovážném (nestabilním) stavu ovšem neznamená, že do něj systém automaticky přejde.
Zda k této přeměně dojde, závisí na rychlosti procesů, kterými se přechod odehrává. Rychlost některých
procesů je tak vysoká, že k nim dojde prakticky okamžitě (neutralizace kyseliny zásadou), u jiných je třeba několik
hodin nebo dnů (rozpouštění plynů ve vodě), některé se odehrávají v časové škále, která je srovnatelná
s geologickou historií Země (rozpad některých radioaktivních prvků) a k některým vůbec nedojde.
Při studiu přírodních systémů můžeme rozlišit dvě základní situace:
– systém je ve stabilním stavu a bude v něm setrvávat
– systém je v nestabilním stavu a bude mít tendenci tento stav změnit.
Pokud je v nestabilním stavu, pak je zajímavé nejen to, kam bude směřovat (do stavu s nejvyšší možnou celkovou
entropií), ale také zda a za jak dlouho tohoto stavu dosáhne. Zatímco směr změn je možné určit pomocí
termodynamického aparátu, studiem rychlosti změn a jejich mechanismem se obecně zabývá dynamika a detailněji
u chemických přeměn kinetika.
Zde je třeba poznamenat, že dosažení termodynamické rovnováhy v přírodních systémech je z dlouhodobého
hlediska výjimečným stavem. Termodynamickou rovnováhu můžeme najít v malých systémech, jako jsou například
v plynokapalné uzavřeniny v minerálech, pórové vody hluboko uložených hornin, jednotlivé části půdního
profilu či sedimentu.
Planeta Země a její jednotlivé součásti se v čase mění nebo jsou dlouhodobě udržovány v nerovnovážném
stavu. Příčinou jsou vnější podmínky. Prakticky nikde nejsou jednotlivé části atmosféry právě nasyceny vůči
vodní páře. V některých oblastech je atmosféra vůči vodní páře nenasycená a voda se z vodních ploch i krajiny
odpařuje, v jiných částech je naopak přesycená, kondenzuje a je z atmosféry vylučována v podobě srážek. Příčinou
je nerovnoměrné a stále se měnící prohřívání povrchu Země. Většina minerálů vyvřelých hornin je
v povrchových podmínkách nestabilní (živce, slídy) a postupně se mění na minerály stabilní (jílové minerály).
Tak by bylo možné pokračovat dále prakticky s jakoukoliv součástí planety. Funkce některých přírodních systémů
je právě umožněna udržováním nerovnováhy. Příkladem může být výše zmíněný koloběh vody v přírodě
a nakonec i funkce všech živých organismů. Ty začnou směřovat k dosažení termodynamické rovnováhy až po
jejich smrti.
6
2.1 KONCEPCE REZERVOÁRŮ
Jako velmi užitečná se při studiu přírodních systémů ukazuje koncepce rezervoárů. Obecně můžeme jakýkoliv
studovaný systém rozdělit na rezervoáry, které obsahují látky a ve kterých se odehrávají chemické přeměny.
Mezi těmito rezervoáry pak dochází k toku látek – rezervoáry si látky vyměňují. Koncept rezervoárů není nijak
omezen celkovou velikostí systémů.
a b
c d
Obr. 2.1 Koncepce systémů a rezervoárů. (a) Celý zemský systém může být rozdělen na čtyři velké rezervoáry, mezi kterými se vyměňují
látky – atmosféra, litosféra, hydrosféra a litosféra. (b) Celý systém hydrosféry je možné rozdělit na jednotlivé rezervoáry, které obsahují
vodu a znázornit vzájemné toky mezi nimi. (c) Každý rezervoár je charakterizován hmotným obsahem látky a toky jsou vyjádřeny jako
přenesené množství látky za časovou jednotku. (d) Globální biogeochemický cyklus uhlíku v podobě „krabičkového“ modelu. Krabičky
představují jednotlivé rezervoáry uhlíku (různě vázaný uhlík), šipky znázorňují toky uhlíku mezi rezervoáry. Čísla uvnitř krabiček značí
odhadované množství uhlíku v Gt (1015
g C). Čísla u jednotlivých toků jsou odhadované toky uhlíku mezi rezervoáry v Gt C za rok. Jednotlivé
odhady se velmi liší a dosud neexistuje všeobecně přijímaný model. Odhadované antropogenní ovlivnění globálního cyklu je
znázorněno v levé části diagramu a představuje spalování fosilních paliv a vypalování pralesů. Pro konstrukci modelu globálního uhlíkového
cyklu byla modifikována data Siegenthalera a Sarmienta (1993) a Kwona a Schnoora (1995).
Planetu Zemi můžeme rozdělit na čtyři velké rezervoáry, mezi kterými se vyměňuje hmota. Jsou jimi atmosféra,
litosféra, hydrosféra a biosféra (obr. 2.1a). Každý ze studovaných systémů je možné rozdělit na jednotlivé rezervoáry,
které jsou obvykle znázorněny obdélníky a šipkami znázornit jednotlivé toky mezi nimi. Globální
model pro hydrosféru je znázorněn na obr. 2.1b. Obvykle je hmotný obsah rezervoáru uveden přímo v symbolu
rezervoáru (obdélníku) a toky jsou uvedeny u příslušných šipek. Pro část hydrosféry jsou uvedeny odhadované
7
obsahy a toky na obr. 2.1c. Toto znázornění systému v podobě rezervoárů a toků se označuje jako „box“ (krabičkový)
model systému. Jako příklad box modelu studovaného systému je uveden globální cyklus uhlíku (obr.
2.1d).
Při členění systému na jednotlivé rezervoáry a toky mezi nimi záleží vždy na tom, co studujeme a co nás na
chování systému zajímá. Jako příklad je uvedena konstrukce box modelu pro cirkulaci vody na ostrově (obr.
2.2a). Nejprve je třeba identifikovat jednotlivé rezervoáry a toky. Na modelovém ostrově se po srážkách část
vody přímo z krajiny odpaří (evapotranspirace), část odteče do jezer a nakonec do oceánu. Do atmosféry přispívá
také kromě evapotranspirace výpar z jezer a oceánu. Identifikované toky a rezervoáry jsou uvedeny na
obr. 2.2b a jejich převedení do box modelu na obr. 2.2c. Pro úplnost je třeba box model doplnit číselnými údaji
(obr. 2.2d). Takto připravený model pak slouží k vytvoření numerického modelu studovaného systému a studiu
jeho dynamiky.
a b
c d
Obr. 2.2 Konstrukce box modelu pro studovaný systém. (a) Cirkulace vody na studovaném ostrově je určována několika základními toky. (b)
Tyto toky zajišťují přenos vody mezi čtyřmi rezervoáry. (c) Identifikovaných toků a rezervoárů do box modelu. (d) Konečný box model cirkulace
vody na ostrově s vyznačení obsahu rezervoárů a velikosti toků.
Hranice mezi jednotlivými rezervoáry může být fyzická, ale může být také pouze virtuální. Jako příklad je možné
uvést model karbonátových látek ve vodě. Ten je tvořen čtyřmi formami karbonátových látek (speciemi),
mezi kterými je přenášen uhlík
 
31 2
1 2 3
2
2 2 2 3 3 3CO H O H CO H HCO H CO
rr r
r r r
aq
 
  
   
   (2.1)
Specie uhlíkatých látek ve vodě můžeme považovat za samostatné rezervoáry, mezi kterými dochází k toku
uhlíkatých látek. Místo hmotného obsahu rezervoáru je pak pro rezervoár charakteristická koncentrace příslušné
specie.
Tok látky z rezervoáru je určitým způsobem úměrný obsahu rezervoáru. Čím více bude vody v krajině, tím více
se jí odpaří, čím více bude vody v jezeře, tím více jí odteče a čím vyšší bude koncentrace rozpuštěného CO2(aq),
odtok
z jezer
a podzemní
vody
výpar
z oceánu
evapotranspirace
výpar z jezer
výpar z toků
srážky
srážky
evapotranspirace výpar
výpar
odtokodtok
jezera a toky oceán
atmosféra
krajina
jezera
toky
atmosféra
krajina oceán
srážky
evapotranspirace výpar
výpar
odtokodtok
J = 0,2 t
A = 5 t
K = 6 t O = 30 000 t
rKA = 0,6 t den–1
rKJ = 0,3 t den–1
rJO = 0,3 t den–1
rJA = 0,05 t den–1
rJA = 0,5 t den–1
rKA = 0,5 t den–1
8
tím více se jej přemění na kyselinu uhličitou H2CO3. Nejjednodušším vyjádřením tohoto vztahu je přímá úměra
(lineární vztah) v podobě
Yr kX (2.2)
kde rY je tok látky X z rezervoáru Y, který je úměrný jejímu hmotnému obsahu nebo koncentraci X. Konstanta
úměrnosti k se nazývá rychlostní konstanta. Například pro odtok vody z jezera do oceánu bude platit v případě
přímé úměry
jo jo jr k m (2.3)
kde rjo je množství vody, které odteče z jezera do oceánu za jednotku času, mj je množství vody v jezeře a kjo je
rychlostní konstanta pro odtok vody z jezera. Pro přeměnu rozpuštěného oxidu uhličitého CO2(aq) na nedisociovanou
kyselinu uhličitou H2CO3 bude platit
  22+1 1 H OCO aq
r k a a (2.4)
kdy je tok tedy rychlost přeměny rozpuštěného oxidu uhličitého na kyselinu uhličitou úměrný aktivitě rozpuštěného
CO2 a aktivitě vody. Konstantou úměrnosti je rychlostní konstanta k+1. Systémy, ve kterých jsou toky
přímo úměrné množství látky nebo její koncentraci, se označují jako lineární.
Ne všechny toky musí být přímo úměrné množství nebo koncentraci látky. Mezi tokem a koncentrací či množstvím
látky může panovat i nelineární vztah například v podobě
2
Yr kX (2.5)
kde je tok látky X z rezervoáru Y úměrný druhé mocnině množství látky X. Takové systémy se pak označují jako
nelineární. Skutečný zájemný vztah mezi tokem a množstvím či koncentrací látky je třeba určit experimentálně
nebo na základě pozorování reálného přírodního systému.
Z jednotlivých toků, které vstupují a vystupují do a z rezervoáru je pak možné sestavit celkovou bilanci toku
látky pro rezervoár
výsledný i j
i j
r r r   (2.6)
kde ri jsou vstupní toky a rj jsou výstupní toky. Pokud bude výsledný tok kladný, pak se látka v rezervoáru hromadí,
v opačném případě látka v rezervoáru ubývá. Výsledný tok je zároveň roven změně obsahu rezervoáru za
jednotku času. Například pro rezervoár C pak platí
C
i j
i j
dm
r r
dt
   (2.7)
kde dmC/dt je změna množství látky v rezervoáru C za jednotku času. Obdobný vztah bude platit i pro změnu
koncentrace či aktivity.
9
2.2 JEDNOREZERVOÁROVÝ SYSTÉM
Nejjednodušším systémem je jednorezervoárový systém s konstantním vstupem látek a výstupem, který je
úměrný hmotnému obsahu látky rezervoáru.
V rezervoáru je určité množství látky X. Do rezervoáru vstupuje látka X konstantní rychlostí – za jednotku času
vstoupí vždy stejné množství. Výstup z rezervoáru je závislý na množství látky X v rezervoáru a je tomuto
množství přímo úměrný. Znamená to, že čím více bude látky X v rezervoáru, tím více jí vystoupí a naopak. Uvedený
systém s vyznačením jednotlivých veličin je schematicky uveden na obr. 2.3.
Obr. 2.3 Jednorezervoárový systém s konstantním vstupem a výstupem úměrným množství látky X v rezervoáru.
2.2.1 ČASOVÝ VÝVOJ
Protože se vstupující množství látky s časem nemění, můžeme je označit jako konstantu A v jednotkách např.
mol s
–1
(do rezervoáru vstoupí A molů látky X za sekundu). Počáteční množství látky X v rezervoáru označíme
jako X0 v jednotkách mol. Tok látky X do rezervoáru rvstup je pak roven
vstupr A (2.8)
Tok látky z rezervoáru rvýstup je přímo úměrný množství látky v rezervoáru, což je možné vyjádřit vztahem
výstupr kX (2.9)
kde k je rychlostní konstanta. Změna obsahu látky X v rezervoáru v závislosti na čase bude dána rozdílem toku
látky X do rezervoáru a toku látky X z rezervoáru
výsledný vstup výstupr r r  (2.10)
Kladná hodnota výsledného toku rvýsledný bude odpovídat hromadění látky X v rezervoáru, záporná jejímu ubývání.
V diferenciální podobě je možné změnu množství látky X v rezervoáru vyjádřit rovnicí
vstup výstup
dX
r r
dt
  (2.11)
Dosazením za jednotlivé toky z rovnic (2.8) a (2.9) do rovnice (2.11) dostáváme
dX
A kX
dt
  (2.12)
X
rX =
rvstup – rvýstup
rvýstuprvstup
10
Časový vývoj koncentrace látky X v rezervoáru získáme integrací uvedené diferenciální rovnice v příslušných
mezích. Nejprve provedeme separaci proměnných
dX
dt
A kX


(2.13)
a následně můžeme integrovat
0 0
X t
X
dX
dt
A kX

  (2.14)
kde X0 je množství látky X v počátečním čase t = 0. Integrací v příslušných mezích dostáváme
   0
0
1
ln
X t
X
A kX t
k
     (2.15)
     0ln ln 0A kX A kX k t      (2.16)
Dalšími úpravami pak postupně obdržíme množství látky X v rezervoáru jako funkci času
0
ln
A kX
kt
A kX

 

(2.17)
0
ktA kX
e
A kX



(2.18)
 0
kt
A kX A kX e
   (2.19)
 0
kt
kX A kX e A
    (2.20)
0
ktA A
X X e
k k
 
   
 
(2.21)
Ostatní veličiny, které ve vztahu (2.21) vystupují, zůstávají konstantní.
11
Obr. 2.4 Vývoj množství látky v jednorezervoárovém systému s konstantním vstupem a výstupem přímo úměrným množství látky
v rezervoáru. Vstupní parametry: rychlost vstupu látky X do rezervoáru A = 2 g hod–1
, rychlostní konstanta výstupu látky z rezervoáru k
= 0,5 hod–1
, počáteční množství látky X v rezervoáru X0 = 18 g. Doba odezvy 1 je definována jako pokles odchylky od stacionárního stavu
na 10 % původní hodnoty, doba odezvy 2 je definována jako pokles odchylky na 10 % hodnoty stacionárního stavu.
2.2.2 STACIONÁRNÍ STAV
Z jednoduchého posouzení podmínek v systému je zřejmé, že množství látky X v rezervoáru bude směřovat do
stavu, kdy se vstupní a výstupní toky vzájemně vyrovnají a množství látky X v rezervoáru se dále nebude měnit.
Pro tento stav pak musí platit
vstup výstupr r (2.22)
S 0výsledný
dX
r A kX
dt
    (2.23)
Odtud pak snadno získáme pro množství látky X v rezervoáru za stacionárního stavu vztah
S
A
X
k
 (2.24)
Ke stejnému výsledku bychom měli dojít i z rovnice (2.21) pro dostatečně dlouhý čas. Pokud se bude čas t blížit
k nekonečnu t (nebo bude alespoň dostatečně dlouhý), pak se bude hodnota členu
kt
e
blížit k nule
1 0kt kt
e e
  a hodnota X se bude blížit k A/k. Tím je hodnota stacionárního stavu potvrzena. Rovnici
(2.21) je možné také přepsat do tvaru
 S 0 S
kt
X X X X e
   (2.25)
kde X0 – XS je odchylka množství látky X v rezervoáru na počátku tedy v čase t = 0.
stacionární
stav
doba
odezvy 1
doba
odezvy 20
2
4
6
8
10
12
14
16
18
20
0 5 10 15 20
X(g)
t (hod.)
12
2.2.3 DOBA ZDRŽENÍ
Doba zdržení tzdržení (residence time) je definována jako čas, za který se vymění celý obsah rezervoáru vzhledem
ke sledované látce. Je zřejmé, že má význam ji uvádět pro ustálený tedy stacionární stav, kdy jsou vstupy a
výstupy vyrovnané. Pokud by tomu tak nebylo, pak se bude v průběhu výměny látky v rezervoáru zároveň
měnit i její množství a tím i doba nutná pro její výměnu. Pro dobu zdržení za stacionárního stavu pak bude
platit
S S
zdržení
vstup výstup
X X
t
r r
  (2.26)
Pro jednorezervoárový systém s konstantním vstupem a výstupem přímo úměrným hmotnosti látky
v rezervoáru získáme dobu zdržení látky tak, že do rovnice (2.26) dosadíme za vstupní nebo výstupní tok a
stacionární stav příslušné hodnoty
S 1
zdržení
A
X kt
A A k
   (2.27)
nebo
S
S
1
zdržení
X
t
kX k
  (2.28)
Doba zdržení je nepřímo úměrná rychlostní konstantě. Znamená to, že čím bude rychlostní konstanta větší, tím
vyšší bude celková rychlost procesu a tím kratší bude doba zdržení.
2.2.4 DOBA ODEZVY
Dobu odezvy todezvy (response time) můžeme definovat jako čas, za který se vychýlení od stacionárního stavu
„dostatečně“ zmenší. Záleží na nás, co budeme považovat za dostatečné zmenšení odchylky a jaká bude volba
srovnávací hladiny. Buď může být měřítkem (i) zmenšení samotné odchylky nebo může být srovnávací hladinou
(ii) poměr velikosti odchylky vůči stacionárnímu stavu:
(i) V prvním případě budeme hodnotit velikost druhého členu v rovnici (2.21) nebo (2.25). Velikost
odchylky od stacionárního stavu je dána výrazem
 0 0 S
kt ktA
X e X X e
k
  
   
 
(2.29)
Její velikost s časem klesá úměrně výrazu
kt
e
. Pokud má například dosáhnout hodnoty deseti
procent původního vychýlení, pak musí platit
0,1odezvykt
e

 (2.30)
ln0,1odezvykt  (2.31)
ln0,1 2,303
odezvyt
k k
   (2.32)
13
Po tomto čase bude původní odchylka od stacionárního stavu utlumena na jednu desetinu původní
hodnoty. V tomto případě nebude doba odezvy závislá na vlastní hodnotě stacionárního
stavu. Pokud budeme za dostatečné utlumení původní výchylky považovat jinou hodnotu, stačí ji
doplnit do vztahu (2.30).
(ii) Pokud bude měřítkem velikost odchylky ve srovnání se stacionárním stavem, ke kterému systém
směřuje, pak je její relativní velikost dána podílem odchylky od stacionárního stavu a hodnoty
stacionárního stavu. Pro její relativní velikost musí platit
 0 S 0
S S
1
kt
kt
rel
X X e Xvelikost odchylky
X e
stacionární stav X X

  
     
 
(2.33)
Pro relativní odchylku 10 % od stacionárního stavu pak musí platit
0
S
0,1 1 odezvyktX
e
X
 
  
 
(2.34)
0 S
S
0,1 odezvyktX X
e
X
 
  
 
(2.35)
S
0 S
0,1 odezvyktX
e
X X



(2.36)
S
0 S
0,1
ln odezvy
X
kt
X X
 

(2.37)
0 S
S
1
ln
0,1
odezvy
X X
t
k X

 (2.38)
Obě doby odezvy jsou, stejně jako doba zdržení, nepřímo úměrné rychlostní konstantě. Čím větší bude hodnota
rychlostní konstanty, tím rychleji bude utlumeno vychýlení od stacionárního stavu. Systém má tendenci se
do stacionárního stavu vracet.
2.3 DVOUREZERVOÁROVÝ SYSTÉM
Složitějším systémem je dvourezervoárový systém, ve kterém dochází k obousměrnému toku látky mezi rezervoáry.
Tok látky je v obou směrech přímo úměrný množství látky v rezervoárech (obr. 2.5). Pro tok látek (přenesené
množství látky za jednotku času) mezi rezervoáry pak platí
AB AB Ar k m (2.39)
BA BA Br k m (2.40)
kde rAB je tok látky z rezervoáru A do rezervoáru B a kAB je příslušná rychlostní konstanta, rBA je tok látky
z rezervoáru B do rezervoáru A a kBA je opět rychlostní konstanta. Výsledný tok mezi rezervoáry je z hlediska
rezervoáru A roven
A AB BA AB A BA Br r r k m k m      (2.41)
14
a z hlediska rezervoáru B
B BA AB BA B AB Ar r r k m k m      (2.42)
Pro změnu množství látky v rezervoáru A platí
A
AB A BA B
dm
k m k m
dt
   (2.43)
a pro změnu nožství látky v rezervoáru B
B
AB A BA B
dm
k m k m
dt
  (2.44)
Protože se jedná o uzavřený systém, zůstává celkové množství látky v systému konstantní a po celou dobu je
její celkové množství v systému rovno součtu počátečních množství látky v rezervoárech
A0 B0 A Bm m m m m    (2.45)
Obr. 2.5 Dvourezervoárový systém se vzájemným přenosem látky mezi rezervoáry. Toky látky mezi rezervoáry jsou přímo úměrné jejich
množství.
2.3.1 ANALYTICKÉ ŘEŠENÍ
Ve kterémkoliv okamžiku musí pro množství látky v rezervoáru B podle rovnice (2.45) platit vztah
B Am m m  (2.46)
Dosazením za mB do rovnice (2.43) dostáváme
 A
AB A BA A
dm
k m k m m
dt
    (2.47)
Další úpravou
A
mA
rAB = –kAB × mA
rBA = kBA × mB
B mB
rAB = kAB × mA
rBA = –kBA × mB
15
 A
AB BA A BA
dm
k k m k m
dt
    (2.48)
 
A
AB BA A BA
dm
dt
k k m k m

  
(2.49)
 
A
A0
A
0
AB BA A BA
m t
m
dm
dt
k k m k m

    (2.50)
 
 
AB BA A BA
AB BA AB BA A0 BA
1
ln
k k m k m
t
k k k k m k m
  
 
   
(2.51)
 
 
 AB BA A BA
AB BA
AB BA A0 BA
ln
k k m k m
k k t
k k m k m
  
  
  
(2.52)
 
 
 AB BAAB BA A BA
AB BA A0 BA
k k tk k m k m
e
k k m k m
   

  
(2.53)
     AB BA
AB BA A BA AB BA A0 BA
k k t
k k m k m k k m k m e
 
       (2.54)
 AB BABA BA
A A0
AB BA AB BA
k k tk k
m m m m e
k k k k
  
   
  
(2.55)
Obdobně dostaneme pro množství látky v rezervoáru B výraz
 AB BAAB AB
B B0
AB BA AB BA
k k tk k
m m m m e
k k k k
  
   
  
(2.56)
Pro stacionární stav látky v rezervoáru A pak musí platit
 A
AB BA AS BA0
dm
k k m k m
dt
     (2.57)
BA
AS
AB BA
k
m m
k k


(2.58)
a obdobně pro množství látky v rezervoáru B
AB
BS
AB BA
k
m m
k k


(2.59)
Tento výsledek je v souladu s rovnicemi (2.55) a (2.56), kdy se pro nekonečný čas množství látky mA a mB
asymptoticky blíží k těmto hodnotám. Celkové množství látky v obou rezervoárech se nezávisle na počátečních
hodnotách přerozdělí v opačném poměru jejich rychlostních konstant.
16
Obr. 2.6 Vývoj množství látky ve dvourezerovárovém systému, kde jsou výstupy z rezervoárů přímo úměrné množství látky
v rezervoárech. Vstupní podmínky: počáteční množství látky v rezervoáru A tj. mA0 = 3 mol l–1
, počáteční množství látky v rezervoáru B
tj. mB0 = 7 mol l–1
, rychlostní konstanta přestupu látky z rezervoáru A do rezervoáru B tj. kAB = 0,1 hod–1
, rychlostní konstanta přestupu
látky z rezervoáru B do rezervoáru A tj. kBA = 0,4 hod–1
. Doba odezvy 1 je definována jako pokles odchylek od stacionárního stavu na 10
% původních hodnot a je pro oba rezervoáry stejná. Doba odezvy 2 je definována jako pokles odchylek na 10 % hodnot stacionárních
stavů a pro jednotlivé rezervoáry se liší.
2.3.2 DOBA ZDRŽENÍ
Doba zdržení v rezervoárech je rovna času, za který jsou schopny stacionární toky látku v rezervoáru vyměnit
AS AS
A
A B
zdržení
m m
t
r r
  (2.60)
AS AS
A
AB AS BA BS AB
1
zdržení
m m
t
k m k m k
   (2.61)
a podobně pro rezervoár B
BS BS
B
A B
zdržení
m m
t
r r
  (2.62)
BS BS
B
AB AS BA BS BA
1
zdržení
m m
t
k m k m k
   (2.63)
Doba zdržení látky v rezervoárech je nepřímo úměrná rovnovážným konstantám procesů, kterými látka
z rezervoárů „odtéká“. Čím bude jejich hodnota větší, tím rychleji látka odtéká a tím rychleji bude vyměněna.
cAS
doba
odezvy 1
doba
odezvy 2A
cBS
doba
odezvy 2B
0
1
2
3
4
5
6
7
8
9
0 5 10 15
cA,cB(moll–1)
t (hod.)
17
2.3.3 DOBA ODEZVY
Pro dobu odezvy můžeme opět použít obdobná kriteria jako u jednorezervoárového systému. Vychýlené
množství látky v rezervoáru od stacionárního stavu je dáno druhým členem v rovnicích (2.55) resp. (2.56), které
mají po dosazení za členy, které odpovídají stacionárním množstvím látek, tvar
   AB BA
A AS A0 AS
k k t
m m m m e
 
   (2.64)
   AB BA
B BS B0 BS
k k t
m m m m e
 
   (2.65)
Odchylky jsou rovny hodnotám
   AB BA
A0 AS
k k t
m m e
 
 (2.66)
   AB BA
B0 BS
k k t
m m e
 
 (2.67)
kde člen v závorce představuje vychýlení v čase t = 0 (rozdíl původní koncentrace a koncentrace ve stacionárním
stavu) a exponenciální člen ukazuje „vymírání“ této odchylky s časem. Pokud má výchylka dosáhnout například
10 % své původní hodnoty, pak musí platit
 AB BA
0,1odezvyk k t
e
 
 (2.68)
AB BA
1
ln0,1odezvyt
k k
 

(2.69)
a doba odezvy bude pro oba rezervoáry stejná.
Při kriteriu 10% odchylky z hodnoty stacionárního stavu musí platit pro rezervoár A
   AB BAA0 AS
AS
0,1
k k tm m
e
m
 
 (2.70)
AS
A
AB BA A0 AS
0,11
lnodezvy
m
t
k k m m

 
 
(2.71)
a pro rezervoár B obdobně
BS
B
AB BA B0 BS
0,11
lnodezvy
m
t
k k m m

 
 
(2.72)
Absolutní hodnota ve jmenovateli logaritmu je užita proto, že nás zajímá absolutní odchylka od stacionárního
stavu a logaritmická funkce je definována jen pro kladná čísla. Doba odezvy je opět nepřímo úměrná součtu
rychlostních konstant, bude se však pro jednotlivé rezervoáry lišit v závislosti na tom, jak dalece byly koncentrace
v jednotlivých rezervoárech vychýleny od stacionárního stavu.
18
2.4 TŘÍREZERVOÁROVÝ SYSTÉM
Uvedený předchozí průtočný systém (viz kap. 2.2) a systém s obousměrnými toky (kap. 2.3), které obsahují
konstantní a v čase proměnlivé toky, můžeme téměř libovolně kombinovat tak, abychom co nejlépe vystihli
podstatné vlastnosti studovaného systému. Trojrezervoárový průtočný systém s konstantním vstupem, obousměrnými
toky mezi rezervoáry a jedním výstupem je uveden na obr. 2.7. Vstup do systému tří následných
rezervoárů je konstantní, ostatní toky jsou přímo úměrné obsahu látky v příslušném rezervoáru.
Obr. 2.7 Diagram třírezervoárového průtočného systému se zpětnými toky mezirezervoáry.
Pro toky látek mezi rezervoáry platí
VAr A (2.73)
AB AB Ar k m (2.74)
BA BA Br k m (2.75)
BC BC Br k m (2.76)
CB CB Cr k m (2.77)
VC VC Cr k m (2.78)
a pro bilance látky v jednotlivých rezervoárech
A
A VA BA AB BA B AB A
dm
r r r r A k m k m
dt
       (2.79)
B
B AB BA CB BC AB A BA B CB C BC B
dm
r r r r r k m k m k m k m
dt
         (2.80)
C
C BC CB VC BC B CB C VC C
dm
r r r r k m k m k m
dt
       (2.81)
Po ustálení stacionárního stavu v systému se množství látky v rezervoárech nemění
A
BA BS AB AS 0
dm
A k m k m
dt
    (2.82)
A
rA = rVA + rBA
– rAB rBA
B
rB = rAB + rCB
– rBA – rBC
C
rC = rBC – rCB
– rVC
rVCrVA
rCB
rBCrAB
19
B
AB AS BA BS CB CS BC BS 0
dm
k m k m k m k m
dt
     (2.83)
C
BC BS CB CS VC CS 0
dm
k m k m k m
dt
    (2.84)
Řešením tohoto systému tří lineárních rovnic pro stacionární množství látky v rezervoárech dostáváme
 CB VC BC VC
AS
AB VC
BA
BC
Ak k k Ak k
m
k k k
 
 (2.85)
 CB VC
BS
BC VC
A k k
m
k k

 (2.86)
CS
VC
A
m
k
 (2.87)
Doby zdržení v rezervoárech jsou pak
AS
A
AB AS AB
1
zdržení
m
t
k m k
  (2.88)
 
BS
B
BA BC BS BA BC
1
zdržení
m
t
k k m k k
 
 
(2.89)
 
CS
C
CB VC CS CB VC
1
zdržení
m
t
k k m k k
 
 
(2.90)
Integrace diferenciálních rovnic (2.79) až (2.81) je natolik složitá, že je výhodnější využít některou z metod
numerické integrace nebo využít specializovaného programového vybavení. Z toho také vyplývá, že dobu odezvy
jednotlivých rezervoárů a celého systému není možné určit obecně, ale je nutné vycházet z konkrétní numerické
simulace.
Vývoj daného systému pro zvolené parametry je uveden na obr. . Množství látky v rezervoárech nemusí do
stacionárního stavu směřovat přímo, ale může projít složitějším vývojem. Jak je patrné u numerického řešení
vývoje systému pro zvolené parametry, jsou výchozí množství látky v rezervoárech B a C vyšší, než je konečný
stacionární stav, Při směřování do stacionárního stavu však množství látek v rezervoárech B a C nejprve poklesne
pod hodnotu stacionárního stavu a teprve poté se k němu začne asymptoticky blížit.
20
3 ZÁVĚR
Obdobně je možné simulovat vývoj i u vícererervoárových systémů. Zatímco vyjádření jednotlivých toků a
bilance toku hmoty v rezervoárech je relativně jednoduchá, analytické řešení (integrace diferenciálních rovnic)
je obtížnou a v mnoha případech nemožnou úlohou. Numerické řešení systému diferenciálních rovnic pak
poskytuje velmi výhodnou alternativu, kterou je možné realizovat i v běžném tabulkovém kalkulátoru.
Důležitou vlastností dynamických systémů, a takových je v přírodě většina, je jejich výrazná tendence směřovat
ke stacionárnímu stavu a v tomto stavu setrvávat. Tyto systémy mají vždy jeden stacionární stav, ke kterému
směřují. Ten je určen vzájemným poměrem rychlostních konstant. Pokud je celkové systém uzavřený, to znamená,
že si rezervoáry vyměňují látky pouze mezi sebou a z vnějšího prostředí žádná látka do celého systému
rezervoárů nevstupuje nebo do vnějšího prostředí neodchází, pak systém nakonec spěje do termodynamické
rovnováhy.
21
4 PŘÍKLADY
4.1 PŘÍKLAD 1
Dynamika jednorezervoárového systému s konstantním vstupem a výstupem přímo úměrným množství (koncentraci)
látky v rezervoáru.
4.1.1 ZADÁNÍ:
Proveďte simulaci časového vývoje množství látky X v rezervoáru, který má konstantní vstup A, výstup je přímo
úměrný množství látky v rezervoáru a konstanta úměrnosti je rovna k. Vstupní parametry jsou uvedeny v tab.
4.1.
Tab. 4.1 Vstupní parametry pro model jednorezervoárového systému s konstantním vstupem a výstupem přímo úměrným množství
látky v rezervoáru. Význam jednotlivých parametrů: A – konstantní vstupní tok do rezervoáru, k – rychlostní konstanta výstupního toku
z rezervoáru, X0 – počáteční koncentrace látky X v rezervoáru.
parametr hodnota jednotky
A = 2 g hod
–1
k = 0,5 hod
–1
X0 = 18 g
Vypočítejte dobu zdržení látky X v rezervoáru a dobu odezvy. Pro výpočet doby odezvy použijte jako kritérium
utlumení „vzruchu“ pokles na 10 % původní výchylky a na 10 % hodnoty stacionárního stavu.
4.1.2 ŘEŠENÍ
Model jednorezervoárového systému s vyznačením vstupních parametrů je uveden na obr. 4.1.
Obr. 4.1 Jednorezervoárový systém s konstantním vstupem a výstupem úměrným koncentraci (množství) látky X v rezervoáru.
Změna množství látky v rezervoáru je popsána diferenciální rovnicí
vstup
A =
2 g
rezervoár
X0 = 18 g
výstup
k×X
k =
0,5 hod.–1
22
dX
A kX
dt
  (3.1)
Její integrací získáme pro časový vývoj množství látky X v rezervoáru rovnici
 S 0 S
kt
X X X X e
   (3.2)
kde X je množství látky v rezervoáru v čase t, XS je stacionární množství látky v rezervoáru a X0 je počáteční
množství látky v rezervoáru tedy v čase t = 0. Časový vývoj v systému je uveden na obr. 4.2.
a b
Obr. 4.2 Časový vývoj množství látky X v rezervoáru. (a) Modrá linie ukazuje časový vývoj vypočítaný podle rovnice (2.25), červené
kroužky numerickou simulaci vývoje Eulerovou metodou s časovým krokem 0,5. (b) Vyznačení stacionárního stavu a časů odezvy definovaných
jako 10 % původního vychýlení (1) a jako 10 % odchylky z hodnoty stacionárního stavu (2).
Pro numerickou simulaci časového vývoje množství látky X v rezervoáru použijeme Eulerovu metodu, která
vychází z rovnice (3.1). Místo nekonečně malých změn množství látky v rezervoáru dX za nekonečně malý přírůstek
času dt přejdeme ke konečným změnám X za určitý čas t. Přitom předpokládáme, že pro daný časový
úsek t se množství látky X mění s časem lineárně, což je pro dostatečně malý časový úsek i u „rozumně“ se
chovající nelineární rovnice splněno. Úpravou na tvar
 X A kX t    (3.3)
získáváme změnu množství látky v rezervoáru. Výpočet změn v jednotlivých krocích pak poskytuje pro dané
podmínky a t = 0,5 hodnoty uvedené v tab. 4.2.
Tab. 4.2 Numerická simulace vývoje množství látky X v rezervoáru pro výchozí podmínky uvedené v tab. 4.1. Zvolen byl časový krok 0,5
hod.
krok čas X X
hod. g g
0 0,0 18 -3,50
0
2
4
6
8
10
12
14
16
18
20
0 5 10 15 20
X(g)
t (hod.)
stacionární
stav
časodezvy
1
časodezvy
20
2
4
6
8
10
12
14
16
18
20
0 5 10 15 20
X(g)
t (hod.)
23
1 0,5 14,50 -2,63
2 1,0 11,88 -1,97
… … … …
29 14,5 4,00 0,00
30 15,0 4,00 0,00
Hodnoty vývoje množství látky v rezervoáru jsou pro časový úsek 0–15 hod. vyneseny na obr. 4.2. Je patrné, že
je celkově dosaženo dobré shody analytického a numerického řešení, přesto převyšuje rozdíl skutečné a Eulerovou
metodou numericky simulované hodnoty v počátečních fázích vývoje více než 0,5 g.
Stacionárního stavu množství látky X je dosaženo po vyrovnání vstupního a výstupního toku
SA kX (3.4)
1
S 1
2 g hod
4 g
0,5 hod
A
X
k


   (3.5)
Doba zdržení je čas, za který toky ve stacionárním stavu vymění obsah rezervoáru
S
1
1 1
2 hod
0,5 hod
zdržení
A
X kt
A A k 
     (3.6)
Doba odezvy může být definována jako čas, za který odchylka od stacionárního stavu poklesne na určitou hodnotu.
Pokud to má být 10 % původní hodnoty, pak pro uvedený systém platí
0,1odezvykt
e

 (3.7)
ln0,1 ln0,1
4,61 hod.
0,5
odezvyt
k
     (3.8)
Pokud bude doba odezvy určena dosažením hodnoty, která se bude lišit o 10 % stacionárního stavu, pak pro
uvedený systém platí
 0 S
S
0,1
odezvykt
X X e
X


 (3.9)
0 S
S
1 1 18 4
ln ln 7,11 hod.
0,1 0,5 0,1 4
odezvy
X X
t
k X
 
  

(3.10)
Množství látky X v rezervoáru se asymptoticky blíží k hodnotě stacionárního stavu 4, doba zdržení látky X
v rezervoáru je rovna 2 g. Za čas 4,61 hod. poklesne původní výchylka od stacionárního stavu na 10 % původní
hodnoty a za 7,11 hod. bude odchylka tvořit 10 % hodnoty stacionárního stavu. Všechny tři hodnoty jsou vyznačeny
na obr. 4.2b.
4.1.2.1 POZNÁMKA
24
Pokud budeme modelovat vývoj koncentrace látky X v rezervoáru, pak musíme brát v úvahu ještě objem rezervoáru.
Toky do a z rezervoáru, které jsou definovány jako látková množství za jednotku času, je pak nutné
přepočítávat na příslušný objem, abychom dostali koncentrace látky X v rezervoáru.
25
4.2 PŘÍKLAD 2
Dvourezervoárový systém s obousměrným tokem látky přímo úměrným koncentraci látky v rezervoárech
4.2.1 ZADÁNÍ
Mezi dvěma rezervoáry je přenášena látka oběma směry a tok látky je přímo úměrný koncentraci látky
v rezervoárech (obr. 4.3). Pro jednoduchost budeme předpokládat, že oba rezervoáry mají jednotkový objem 1
l (pokud by byl objem odlišný, museli bychom změny množství látky v rezervoárech, způsobené toky, přepočítávat
na odpovídající koncentrace). Počáteční hodnoty jsou uvedeny v tab. 4.3.
Tab. 4.3 Vstupní parametry pro model dvourezervoárového systému s obousměrnými toky přímo úměrnými koncentracím látky v
rezervoárech.
parametr hodnota jednotky
cA0 3 mol l
–1
kAB 0,1 hod
–1
cB0 7 mol l
–1
kBA 0,4 hod
–1
4.2.2 ŘEŠENÍ
Uvedený dvourezervoárový systém s vyznačením rychlosti jednotlivých toků je uveden na obr. 4.3.
Obr. 4.3 Dvourezervoárový systém se vzájemným přenosem látky mezi rezervoáry. Toky látky mezi rezervoáry jsou přímo úměrné jejich
koncentracím.
Časový vývoj koncentrace složky v rezervoáru A je popsán diferenciální rovnicí
 A
AB A BA A
dc
k c k c c
dt
    (3.11)
A
cA
rAB = –kAB × cA
rBA = kBA × cB
B
cB
rAB = kAB × cA
rBA = –kBA × cB
26
a časový vývoj koncentrace složky v rezervoáru B diferenciální rovnicí
 B
BA B AB B
dc
k c k c c
dt
    (3.12)
kde c je celková koncentrace složky v systému a je rovna
A0 B0 A Bc c c c c    (3.13)
Stacionární stav, ke kterému bude systém směřovat, je charakterizován vyrovnanými toky, tedy
 A
AB BA AS BA0
dc
k k c k c
dt
     (3.14)
 B
AB BA BS AB0
dc
k k c k c
dt
     (3.15)
Odtud pak
1BA
AS
AB BA
0,4
10 8 mol l
0,1 0,4
k
c c
k k

   
 
(3.16)
1AB
BS
AB BA
0,1
10 2 mol l
0,1 0,4
k
c c
k k

   
 
(3.17)
Integrací rovnic (3.11) a (3.12) získáme pro časový vývoj koncentrace látky v rezervoárech A a B vztahy
   AB BA
A AS A0 AS
k k t
c c c c e
 
   (3.18)
   AB BA
B BS B0 BS
k k t
c c c c e
 
   (3.19)
Odchylka od stacionárních koncentrací cAS a cBS „odumírá“ v obou rezervoárech úměrně velikosti členu
 AB BAk k t
e 
, který obsahuje obě rychlostní konstanty a s časem exponenciálně klesá. Vývoj koncentrací je pro
dané výchozí podmínky uveden na obr. 4.4a.
Numerickou simulaci vývoje provedeme Eulerovou metodou s využitím upravených vztahů (3.11) a (3.12)
 A AB A BA Ac k c k c c t        (3.20)
 B BA B AB Bc k c k c c t        (3.21)
kde cA a cB jsou změny koncentrací látky v rezervoárech za čas t. Výsledky numerické simulace jsou uvedeny
na obr. 4.4a a část číselných hodnot v tab. 4.4.
a b
27
Obr. 4.4 Vývoj koncentrací ve dvourezervoárovém systému s obousměrným tokem přímo úměrným koncentrace látky v rezervoáru.
Vývoj koncentrace v rezervoáru A je vyznačen modrou linií, vývoj koncentrace látky v rezervoáru B je znázorněn zelenou linií. (a) Numerická
simulace vývoje s časovým korkem 0,5 hod. je vyznačena kroužky. (b) Doby odezvy v rezervoárech – doba odezvy 1 je čas, za který
poklesne odchylka od stacionárního stavu na 10 % původní hodnoty (je stejná pro oba rezervoáry, doba odezvy 2 je čas, za který se
koncentrace v rezervoárech liší o 10 % hodnoty stacionárního stavu (je pro rezervoáry A a B odlišná).
Tab. 4.4 Výsledky numerické simulace vývoje koncentrací ve dourezervoárovém systému s obousměrným tokem.
krok t cA num dcA num cB num dcB num
hod. mol l
–1
mol l
–1
mol l
–1
mol l
–1
0 0,0 3 1,25 7 -1,25
1 0,5 4,25 0,94 5,75 -0,94
2 1,0 5,19 0,70 4,81 -0,70
… … … … … …
28 14,0 8,00 0,00 2,00 0,00
29 14,5 8,00 0,00 2,00 0,00
30 15,0 8,00 0,00 2,00 0,00
Doba zdržení látky v rezervoárech je dána za stacionárního stavu výrazy
AS
A
AB AS AB
1 1
10 hod.
0,1
zdržení
c
t
k c k
    (3.22)
BS
B
BA BS BA
1 1
2,5 hod.
0,4
zdržení
c
t
k c k
    (3.23)
Doba odezvy, definovaná jako pokles odchylky od stacionárního stavu na 10 % původní hodnoty je pro oba
rezervoáry stejná
0
1
2
3
4
5
6
7
8
9
0 5 10 15
cA,cB(moll–1)
t (hod.)
cAS
doba
odezvy 1
doba
odezvy 2A
cBS
doba
odezvy 2B
0
1
2
3
4
5
6
7
8
9
0 5 10 15
cA,cB(moll–1)
t (hod.)
28
 
AB BA
1 1
ln0,1 2,303 4,61
0,1 0,4
odezvyt
k k
      
 
(3.24)
Doba odezvy, definovaná jako odchylka 10 % z hodnoty stacionárních stavů koncentrací látky v rezervoárech,
je rovna
AS
A
AB BA A0 AS
0,11 1 0,1 8
ln ln 3,67 hod.
0,1 0,4 3 8
odezvy
c
t
k k c c

    
   
(3.25)
BS
B
AB BA B0 BS
0,11 1 0,1 2
ln ln 6,44 hod.
0,1 0,4 7 2
odezvy
c
t
k k c c

    
   
(3.26)
29
4.3 PŘÍKLAD 3
Systém se třemi za sebou jdoucími rezervoáry s obousměrnými toky mezi nimi
4.3.1 ZADÁNÍ
Nefelín se rozpouští ve vodě v kyselém prostředí při pH = 3,65 a teplotě 60 °C za uvolnění iontů hliníku podle
rovnice
NaAlSiO4 + 4 H
+
 2 H2O + Na
+
+ Al
3+
+ SiO2(aq) (3.27)
Rozpuštěný hliník se pak může srážet v podobě gibbsitu podle rovnice
Al
3+
+ 3 H2O  Al(OH)3 + 3 H
+
(3.28)
Simulujte časový vývoj koncentrace hliníku v 1 l vody:
(a) do které bylo přidáno 284,1 g nefelínu
(b) do které bylo přidáno 78 g gibbsitu
(c) do které bylo přidáno 284,1 g nefelínu a 78 g gibbsitu.
Povrch nefelínu je 0,1 m
2
g
–1
a povrch gibbsitu je 0,3 m
2
g
–1
. Rychlostní konstanty rozpouštění a srážení nefelínu
a gibbsitu jsou následující:
rychl. konstanta hodnota jednotky
knefelín rozp 2,332 mg m
–2
hod
–1
knefelín sráž 0,000234 m
–2
l
–1
hod
–1
kgibbsit rozp 0,0700 mg m
–2
hod
–1
kgibbsit sráž 0,0757 m
–2
l
–1
hod
–1
Pro jednoduchost předpokládejte, že jsou koncentrace rozpuštěného sodíku a SiO2(aq) udržovány na konstantní
hodnotě.
POZNÁMKA
Příklad vychází z upravené analogie rozpouštění nefelínu, který uvádí Lasaga (1981).
Lasaga A. C. (1981): Rate laws of chemical reactions. In: A.C. Lasaga and J.R. Kirkpatrick, Editors, Kinetics of
Geochemical Processes Vol. 8, Mineralogical Society of America (1981), pp. 1–68: str. 10, fig. 1 (Tole, pers.
comm.).
4.3.2 ŘEŠENÍ
Rozpouštění a srážení nefelínu a gibbsitu můžeme znázornit pomocí třírezervoárového systému uvedeném na
obr. 4.5.
30
Obr. 4.5 Třírezervoárový model rozpouštění nefelínu a srážení gibbsitu. Vyznačeny jsou toky, které určují koncentraci rozpuštěného
hliníku ve vodě.
Nejprve rozebereme případy jednotlivých minerálů.
4.3.2.1 ROZPOUŠTĚNÍ NEFELÍNU
Rozpouštění nefelínu probíhá podle rovnice
NaAlSiO4 + 4 H
+
 2 H2O + Na
+
+ Al
3+
+ SiO2(aq) (3.27)
Pokud je udržováno konstantní pH, pak závisí rozpouštění nefelínu a množství uvolněného hliníku výhradně na
povrchu reagujícího minerálu. Tok hliníku do roztoku rozpouštěním nefelínu je
nefelín rozp nefelín nefelín rozpr A k (3.29)
Tok hliníku do nefelínu jeho zpětným srážením
2 H2O + Na
+
+ Al
3+
+ SiO2(aq)  NaAlSiO4 + 4 H
+
(3.30)
je při konstantních koncentracích sodíku a rozpuštěného SiO2(aq) v roztoku závislý pouze na koncentraci hliníku
v roztoku a povrchu minerálu (ten poskytuje místa, kam se může hliník z roztoku připojovat)
3
Alnefelín sráž nefelín nefelín srážr A k c  (3.31)
Výsledný tok hliníku do roztoku je pak dán jejich rozdílem
3+ 3
Al roztok(n) Alnefelín rozp nefelín sráž nefelín nefelín rozp nefelín nefelín srážr r r A k A k c     (3.32)
Tento vztah zároveň udává změnu koncentrace hliníku v závislosti na čase při rozpouštění nefelínu
3
3
Al
Alnefelín nefelín rozp nefelín nefelín sráž
dc
A k A k c
dt

  (3.33)
Stacionární stav je pak dán vztahem
nefelín
rnefelín rozp =
= –knefelín rozp×Anefelín
rnefelín sráž =
=Anefelín ×
×knefelín sráž×cAl3+
roztok
cAl3+
rnefelín rozp =
=knefelín rozp×Anefelín
rnefelín sráž =
= –Anefelínknefelín sráž×cAl3+
rgibbsit rozp =
= kgibbsit rozp×Agibbsit
rgibbsit sráž =
= –Agibbsitkgibbsit sráž×cAl
gibbsit
rgibbsit rozp =
= –kgibbsit rozp×Agibbsit
rgibbsit sráž =
= Agibbsit×
×kgibbsit sráž×cAl3+
31
3
3
SAl
Al
0 nefelín
nefelín nefelín rozp nefelín nefelín sráž
dc
A k A k c
dt

   (3.34)
3
S 1
4Al
2,332
9 972 mg l
2,339 10
nefelín nefelín rozp
nefelín sráž
k
c
k



  

(3.35)
Pro reakci (3.27) je rovnováha určena rovnovážnou konstantou
 
 
   
23 3
2 2 2 3
4 4
4
Na Al SiO Na Al SiO H O
Al
NaAlSiO H H
sp podmíněná
aq aq
K K
   

 
                          
      
(3.36)
která se označuje jako součin rozpustnosti (solubility product). Výrazy v hranatých závorkách jsou aktivity příslušných
složek, v prvním přiblížení je možné použít koncentrace, aktivita (či koncentrace) nefelínu jako čisté
pevné fáze je jednotková. Pokud jsou aktivity vody, sodíku, rozpuštěného SiO2(aq) a pH (tedy aktivita protonů)
udržovány na konstantní hodnotě, pak je můžeme sloučit do konstanty Kpodmíněná a nasycení roztoku vůči nefelínu
je určeno výhradně aktivitou rozpuštěného hliníku. Termodynamické modelování rozpustnosti nefelínu
v programu Geochemist´s Workbench poskytne stejnou hodnotu.
Roztok v rovnováze s nefelínem bude mít koncentraci 9 972 mg l
–1
hliníku. Časový vývoj koncentrace hliníku
v roztoku při rozpouštění nefelínu je možné získat integrací rovnice (3.33) nebo numerickou simulací. Postupné
kroky integrace
3
3
Al
Alnefelín rozp nefelín sráž nefelín nefelín rozp nefelín nefelín sráž
dc
r r A k A k c
dt

    (3.37)
3
3
Al
Alnefelín nefelín rozp nefelín nefelín sráž
dc
dt
A k A k c




(3.38)
   
3
3
0
3Al
Al
0Al
1
ln
t
nefelín nefelín rozp nefelín nefelín sráž
c
nefelín nefelín sráž
A k A k c t
A k



    
(3.39)
 
3Al
3 0
3Al
Al
ln
c
nefelín nefelín rozp nefelín nefelín sráž nefelín nefelín sráž
c
A k A k c A k t



     (3.40)
3
3
Al
0
Al
ln
nefelín nefelín rozp nefelín nefelín sráž
nefelín nefelín sráž
nefelín nefelín rozp nefelín nefelín sráž
A k A k c
A k t
A k A k c



 

(3.41)
 3 3
0
Al Al
nefelín nefelín srážA k t
nefelín sráž nefelín rozp nefelín rozp nefelín srážk c k k k c e 

    (3.42)
3 3
0
Al Al
nefelín nefelín srážA k tnefelín rozp nefelín rozp
nefelín sráž nefelín sráž
k k
c c e
k k
 
 
    
 
(3.43)
3 3
0
Al Al
nefelín nefelín srážA k tnefelín rozp nefelín rozp
nefelín sráž nefelín sráž
k k
c c e
k k
 
 
    
 
(3.44)
32
Podíl rychlostních konstant v rovnici (3.44) je roven podle rovnice (3.35) stacionární koncentraci. Dosazením
3
S
Al
nefelínnefelín rozp
nefelín sráž
k
c
k
 (3.45)
dostáváme
 3 3 3 3
S 0 S
Al Al Al Al
nefelín nefelín srážA k t
c c c c e   

   (3.46)
Numerická simulace jednoduchou Eulerovou metodou poskytuje pro výchozí parametry stejný výsledek (obr.
4.6a).
a b
Obr. 4.6 Simulace vývoje koncentrace rozpuštěného hliníku ve vodě podle analytického řešení (linie) a Eulerovou metodou (body). (a)
Rozpouštění nefelínu: cS-nef – stacionární koncentrace hliníku, která je zároveň rovnovážnou koncentrací. (b) Rozpouštění gibbsitu: cS-nef
– stacionární koncentrace hliníku, která je zároveň rovnovážnou koncentrací.
4.3.2.2 ROZPOUŠTĚNÍ GIBBSITU
Rozpouštění gibbsitu probíhá podle rovnice
Al(OH)3 + 3 H
+
 Al
3+
+ 3 H2O (3.47)
Zároveň dochází k jeho zpětnému srážení podle rovnice
Al
3+
+ 3 H2O  Al(OH)3 + 3 H
+
(3.28)
Celkový tok hliníku mezi gibbsitem a roztokem je dán rozdílem obou procesů.
gibbsit rozp gibbsit gibbsit rozpr A k (3.48)
cS-nef
0,1
1
10
100
1000
10000
0 100 200 300 400 500
cAl3+(mg/l)
čas (dny)
cS-nef nefelín Numericky
cS-gib
0,01
0,1
1
0 0,5 1 1,5 2
cAl3+(mg/l)
čas (dny)
cS-gib analyticky numericky
33
3
Algibbsit sráž gibbsit gibbsit srážr A k c  (3.49)
3+ 3
Al roztok(g) Algibbsit rozp gibbsit sráž gibbsit gibbsit rozp gibbsit gibbsit srážr r r A k A k c     (3.50)
Změna koncentrace rozpuštěného hliníku je pak určena diferenciální rovnicí
3
3
Al
Algibbsit gibbsit rozp gibbsit gibbsit sráž
dc
A k A k c
dt

  (3.51)
Pro stacionární stav pak z rovnice (3.51) odvodíme
3
3
SAl
Al
0 gibbsit
gibbsit gibbsit rozp gibbsit gibbsit sráž
dc
A k A k c
dt

   (3.52)
3
S 1
Al
0,0700
0,925 mg l
0,0757
gibbsit gibbsit gibbsit rozp gibbsit rozp
gibbsit gibbsit sráž gibbsit sráž
A k k
c
A k k


    (3.53)
Pro reakci rovnováhu rozpouštění gibbsitu (3.47) platí součin rozpustnosti
 
 
 
3 33 3
2 2 3
3 3
3
Al H O Al H O
Al
Al OH H H
sp podmíněnáK K
 

 
            
         
(3.54)
Pokud budou aktivita vody a pH udržovány konstantní, pak je lze zahrnout do jedné konstanty Kpodmíněná a stav
nasycení roztoku vůči gibbsitu bude určován výhradně aktivitou rozpuštěného hliníku. Termodynamické modelování
rozpustnosti gibbsitu v programu Geochemist´s Workbench poskytne stejnou hodnotu. Časový vývoj
koncentrace hliníku v roztoku při rozpouštění gibbsituje možné získat integrací rovnice (3.52) nebo numerickou
simulací. Postupné kroky integrace jsou analogické, jako u nefelínu, analytické vyjádření má pak tvar
 3 3 3 3
S 0 S
Al Al Al Al
gibbsit gibbsit srážA k t
c c c c e   

   (3.55)
Výsledky analytického řešení a numerické simulace rozpouštění gibbsitu jsou uvedeny na obr. 4.6b.
4.3.2.3 ROZPOUŠTĚNÍ NEFELÍNU ZA VZNIKU GIBBSITU
Zatím bylo vyřešeno rozpouštění nefelínu až do vzniku roztoku nasyceného vůči nefelínu a rozpouštění gibbistu
až do vzniku roztoku nasyceného vůči gibbsitu. Koncentrace rozpuštěného hliníku se navzájem diametrálně liší
– zatímco nasycení vůči nefelínu bude za daných podmínek dosaženo až při koncentraci rozpuštěného hliníku
9 972 mg l
–1
, vůči gibbsitu bude roztok nasycen už při necelém 1 mg l
–1
rozpuštěného hliníku. Intuitivně budeme
očekávat, že se nefelin bude měnit na gibbsit a koncentrace rozpuštěného hliníku bude udržována na hodnotě,
která je charakteristická pro rozpustnost gibbsitu za daných podmínek, tedy kolem 1 mg l
–1
.
To by platilo (i) v případě dosažení termodynamické rovnováhy, kdy se nakonec všechen nefelín přemění na
gibbsit, nebo (ii) v případě, že bude rychlost srážení rozpuštěného hliníku řádově vyšší, než je rychlost jeho
uvolňování do roztoku z nefelínu. Pokud budou obě rychlosti srovnatelné, pak dosáhne koncentrace hliníku
stacionárního stavu, který bude ležet někde mezi oběma hodnotami.
Tento předpoklad prověříme vyhodnocením toků hliníku do a z roztoku a jeho časovou simulací. Hliník je do
roztoku uvolňován dvěma procesy – rozpouštěním nefelínu a gibbsitu podle rovnic (3.27) a (3.47)
34
NaAlSiO4 + 4 H
+
 2 H2O + Na
+
+ Al
3+
+ SiO2(aq) (3.27)
Al(OH)3 + 3 H
+
 Al
3+
+ 3 H2O (3.47)
Z roztoku je hliník vázán dvěma procesy – srážením nefelínu a gibbsitu podle rovnic (3.30) a (3.28)
2 H2O + Na
+
+ Al
3+
+ SiO2(aq)  NaAlSiO4 + 4 H
+
(3.30)
Al
3+
+ 3 H2O  Al(OH)3 + 3 H
+
(3.28)
Změna koncentrace rozpuštěného hliníku v důsledku těchto čtyř procesů je pak dána diferenciální rovnicí
3
3 3
Al
Al Alnefelín nefelín rozp nefelín nefelín sráž gibbsit gibbsit rozp gibbsit gibbsit sráž
dc
A k A k c A k A k c
dt

     (3.56)
 3
3
Al
Alnefelín nefelín rozp gibbsit gibbsit rozp nefelín nefelín sráž gibbsit gibbsit sráž
dc
A k A k A k A k c
dt

    (3.57)
Integrací rovnice (3.57) dostáváme pro časový vývoj koncentrace hliníku v roztoku vztah
3
3
Al
0
Al
nefelín nefelín rozp gibbsit gibbsit rozp
nefelín nefelín sráž gibbsit gibbsit sráž
nefelín nefelín rozp gibbsit gibbsit rozp
nefelín nefelín sráž gibbsit gibbsit sráž
A k A k
c
A k A k
A k A k
c
A k A k



 

 
  
 nefelín nefelín sráž gibbsit gibbsit srážA k A k t
e
 


(3.58)
Z rovnice (3.57) snadno odvodíme, že pro stacionární stav platí
3
S
Al
nef gib nefelín nefelín rozp gibbsit gibbsit rozp
nefelín nefelín sráž gibbsit gibbsit sráž
A k A k
c
A k A k





(3.59)
Dosazením do rovnice (3.58) dostáváme pro vývoj koncentrace rozpuštěného hliníku
   
3 3 3 3
S 0 S
Al Al Al Al
nefelín nefelín sráž gibbsit gibbsit srážA k A k t
c c c c e   
 
   (3.60)
Pro zadané podmínky bude stacionární koncentrace hliníku rovna 38,2 mg l
–1
, což je výrazně více, než je nasycení
vůči gibbsitu a zároveň výrazně méně, než je stav nasycení vůči nefelínu za daných podmínek. Výsledky
analytického řešení a numerické simulace jsou na obr. 4.7a.
a b
35
Obr. 4.7 Vývoj koncentrace rozpuštěného hliníku při rozpouštění nefelínu a srážení gibbsitu. (a) Pro zadané podmínky. (b) Pokud se
změní poměr rychlostních konstant, pak se stacionární stav posouvá směrem ke stavu nasycení vůči minerálu, který má větší rychlostní
konstantu rozpouštění. Podobný vliv má aktuální povrch minerálů. Pro řešený příklad byl poměr rychlostních konstant knefelín rozp/kgibbsit
rozp = 33,3. Pro srovnání jsou uvedeny stacionární stavy při dalších poměrech rychlostních konstant.
Pokud se změní poměr rychlostních konstant nebo poměr dalších parametrů, které ovlivňují rychlost (množství
minerálů v systému a s tím související povrch; velikost zrn, která opět určuje celkový reagující povrch), pak se
bude také posouvat hodnota stacionárního stavu (obr. 4.7b). Stacionární stav se bude vždy posouvat směrem
k minerálu, který bude mít větší celkovou rychlost uvolňování hliníku do roztoku, ať už k ní dojde v důsledku
změny rychlostní konstanty, množství minerálu, agregátního stavu atd.
Uvedený příklad a model je značně zjednodušený a nezachycuje proces zvětrávání nefelínu v přírodních podmínkách.
V přírodních systémech bude vznikat celá řada dalších meziproduktů, zejména však nebudou fixovány
koncentrace složek, které se reakce účastní a v důsledku reakcí rozpouštění a srážení se bude měnit také pH.
Model však ukazuje jednu důležitou vlastnost přírodních systémů. Ty jsou dynamické, jednotlivé rezervoáry
jsou charakterizovány mnoha vstupy a mnoha výstupy. Zemská kůra je v neustálém pohybu a po celou dobu
vývoje se mění. To jenom my někdy změny nevnímáme, protože jejich časová škála řádově přesahuje délku
lidského života.
Termodynamické modelování přírodních systémů v posledním desetiletí velmi pokročilo a jsme schopni studovat
a modelovat systémy, které se v základních rysech blíží reálným systémům. Přesto je relativně zřídka dosaženo
shody modelu a reálných hodnot, které zjišťujeme v přírodních systémech (například složení povrchových
a podzemních vod v kontaktu s horninami, složení vod v kontaktu s atmosférou atd.). Jejich složení je dáno
ustavením stacionárních stavů a nikoliv termodynamickou rovnováhou, proto často nacházíme vody nenasycené
či přesycené vůči minerálům a horninám, které jsou s nimi v kontaktu. Přesto je termodynamické modelování
velice cenné a užitečné. Ukazuje, kterým směrem se budou změny v přírodních systémech ubírat a
v jakých hranicích se budou pohybovat. A to není málo.
cstac. stav nef-gib
0,1
1,0
10,0
100,0
0 0,5 1 1,5 2
cAl3+(mg/l)
čas (dny)
nef-gib cstac. stav nef-gib numericky
cS-nef
cS-gib
0,1
1
10
100
1000
10000
0 0,5 1 1,5 2
cAl3+(mg/l)
čas (dny)
33,3 1 500 5000
36
5 VZTAH MEZI STACIONÁRNÍM STAVEM A TERMODYNAMICKOU ROVNOVÁHOU
Dvourezervoárový model s toky přímo úměrnými koncentracím (přesněji aktivitám) reagujících látek
v rezervoárech je zároveň modelem jednoduché chemické reakce
A + 2 B C
k
k


(4.1)
s tím, že se vše odehrává v jednom rezervoáru a látky A a B se podle stechiometrie reakce (4.1) mění na látku C
a naopak. Pro toky dopředné a zpětné reakce pak platí
2
A Br k a a  (4.2)
Cr k a  (4.3)
a výsledný tok látek je roven jejich rozdílu
2
A B Cvýslednýr r r k a a k a       (4.4)
Pro stacionární stav, kdy jsou dopředný a zpětný tok vyrovnány, pak platí
2
AS BS CS0výslednýr k a a k a    (4.5)
a
CS
2
AS BS
a k
a a k


 (4.6)
kde index S platí pro aktivity reagujících látek po dosažení stacionárního stavu tedy pro stav, kdy jsou dopředný
a zpětný tok vyrovnány a z makroskopického hlediska v systému nedochází ke změnám (koncentrace zůstávají
v čase konstantní).
Dosažený stav termodynamické rovnováhy neznamená, že chemická reakce ustane, ale že se rychlosti obou
reakci
A + 2 B C
k
 (4.7)
C A + 2 B
k
 (4.8)
vyrovnají. Termodynamická rovnováha je tedy zvláštním případem stacionárního stavu v jednorezervoárovém
systému. Dopředný proces vyprodukuje právě takové množství látek, které je zpětným procesem opět spotře-
bováno.
37
5.1 ARRHENIOVA ROVNICE
Podle empirických měření je velikost rychlostní konstanty pro celou řadu chemických reakcí závislá na teplotě
podle vztahu
aE
RT
k Ae

 (4.9)
kde A je tzv. preexponenciální (frekvenční) faktor a Ea je aktivační energie. Tuto rovnici navrhl Arrhenius již
v roce 1889. Interpretace je následující: Aby došlo k chemické reakci, musí reagující látky překonat určitou
energetickou bariéru o velikosti Ea (obr. 5.1a). Čím bude teplota vyšší, tím vyšší mají reagující molekuly energii
a tím snáze tuto energetickou barieru překonají. Rychlostní konstanta přitom podle vztahu (4.9) roste exponenciálně
s teplotou (pro reakce s aktivační energií kolem 50 kJ/mol velmi zhruba dvakrát při každém zvýšení
teploty o 10 °C).
a b
Obr. 5.1 (a) Látka A musí při chemické reakci přeměny na látku B překonat energetickou bariéru o velikosti Ea A, která se označuje jako
aktivační energie. Rychlostní konstanta této reakce je podle Arrheniovy rovnice úměrná velikosti aktivační energie. Při opačné reakci
musí překonat aktivační energii Ea B, která je větší než Ea A, a proto bude rychlostní konstanta této reakce menší. (b) Vznik aktivovaného
komplexu je spojen s růstem Gibbsovy funkce *
G°+ a *
G°– (poklesem celkové entropie). Rychlostní konstanta přeměny látek A a B na C
bude větší a proto bude při přeměně preferován vznik látky C. Stejný výsledek poskytne i čistě termodynamické posouzení rovnováhy
mezi látkami A a B na jedné straně a látkou C na druhé straně, tedy čistě termodynamický popis pomocí Gibbsovy reakční funkce.
5.2 TERMODYNAMICKÉ ODBOČENÍ
Před další úpravou a odvozením nových rovnic je vhodné připomenout některé termodynamické vztahy. Gibbsova
reakční funkce pro reakci (4.1) je rovna
r p vG G G   (4.10)
reakční
profil
Ea A
A
B
Ea B
energie
reakční koordináta
reakční
profil
*G°+
aktivovaný
komplex
A + 2B
C
*G°–
G°r
Gibbsovafunkce
reakční koordináta
38
kde Gp je hodnota Gibbsovy funkce produktů (látek na pravé straně) a Gv je hodnota Gibbsovy funkce výchozích
látek (látek na levé straně), tedy
C A B2rG G G G    (4.11)
Čárka nad písmenem G značí, že se jedná o hodnotu Gibbsovy funkce pro jeden mol látky. Ta je za určité teploty
a tlaku závislá na složení systému tj. na aktivitě látky
lni i iG G RT a
  (4.12)
kde hodnota s indexem ° označuje standardní hodnotu tj. hodnotu Gibbsovy funkce při jednotkové aktivitě
látky. Dosazením do rovnice (4.11) obdržíme
   C C A A B Bln ln 2 lnrG G RT a G RT a G RT a  
       (4.13)
2
C A B C A2 ln ln lnr BG G G G RT a RT a RT a  
       (4.14)
C
2
A
lnr r
B
a
G G RT
a a

    (4.15)
Argument logaritmu se označuje jako reakční kvocient Q (alternativně také aktivitní kvocient)
C
2
A B
a
Q
a a
 (4.16)
Dosazením do rovnice (4.15) získáváme
lnr rG G RT Q
    (4.17)
Minima Gibbsovy funkce systému je dosaženo v případě, že je pro danou reakci Gibbsova reakční funkce rovna
nule
C
2
A B
0 ln r
r
r r
a
G RT
a a

   (4.18)
kde index r u aktivit látek vyjadřuje, že se jedná o aktivity po dosažení rovnováhy. Pak můžeme rovnici (4.18)
přepsat do tvaru
0 lnrG RT K
   (4.19)
kde
C
2
A B
r
r r
a
K
a a
 (4.20)
Minimum Gibbsovy funkce systému zároveň představuje rovnováhu, systém nebude mít tendenci na svém
stavu nic měnit a bude v něm setrvávat. Reakční kvocient Q je v tomto případě roven rovnovážné konstantě K.
Pro rovnovážnou konstantu pak z rovnice (4.18) platí
39
rG
RT
K e



 (4.21)
Pro reakční kvocient, když není v reakci dosaženo rovnováhy, lze odvodit z rovnice (4.17) výraz
r rG G
RT
Q e

 
 (4.22)
5.3 DYNAMIKA PROCESŮ
Porovnání výrazů pro rovnovážnou konstantu rovnice (4.1)
C
2
A B
r
r r
a
K
a a
 (4.20)
kde aXr označuje rovnovážnou aktivitu látky X a stacionárního stavu podle rovnice (4.6)
CS
2
AS BS
a k
a a k


 (4.6)
kde aXS označuje aktivitu látky X po dosažení stacionárního stavu, ukazuje, že rovnovážné aktivity látek jsou
totožné s aktivitami látek po dosažení stacionárního stavu. To vede k propojení rovnovážné konstanty reakce K
s rychlostními konstantami dopředné a zpětné reakce k+ a k–
k
K
k


 (4.23)
Podle rovnice (4.4) je výsledný tok látek roven
2
A B Cvýslednýr r r k a a k a       (4.4)
Podle rovnice (4.23) můžeme rychlostní konstantu k– nahradit výrazem
k
k
K

  (4.24)
a z rovnice (4.16)
C
2
A B
a
Q
a a
 (4.16)
získáme pro aktivitu aC
2
C A Ba a a Q (4.25)
Dosazením rovnic (4.24) a (4.25) do rovnice (4.4) dostáváme pro výsledný tok látky v reakci (4.1)
2 2 2
A B A A B 1výsledný B
k Q
r k a a a a Q k a a
K K

 
 
    
 
(4.26)
40
Pokud je součin aktivit látek, které se účastní reakce, menší, než odpovídá rovnovážné konstantě, bude podíl
Q/K menší než 1 a výraz v závorce bude mít kladnou hodnotu. Tok látek zleva doprava bude větší než tok látek
zprava doleva a celkově bude reakce probíhat zleva doprava. Pokud bude naopak podíl Q/K větší než 1, bude
tok zleva doprava menší než tok zprava doleva a celkově bude reakce probíhat zprava doleva. Při rovnosti
Q = K jsou aktivity látek v reakci rovny rovnovážným aktivitám, tok látek zleva doprava je stejný, jako tok látek
zprava doleva a celkový tok je nulový.
Termodynamické veličiny jako jsou Gibbsova reakční funkce rG a standardní Gibbsova reakční funkce rG

jsou v rovnici pro výslednou rychlost procesu obsaženy nepřímo ve výrazech pro reakční kvocient Q, viz rovnice
(4.22) a pro rovnovážnou konstantu, viz rovnice (4.21).
Přestože ve výrazu rovnice (4.26) vystupuje pouze rychlostní konstanta pro dopřednou reakci a aktivity látek
na levé straně, z odvození je patrné, že jsou v něm prostřednictvím reakčního kvocientu Q a rovnovážné konstanty
K obsaženy i rychlostní konstanta zpětné reakce a aktivity látek na pravé straně, viz rovnice (4.4) a (4.26)
.
Pro vztah mezi výslednou rychlostí reakcí a Gibbsovou reakční funkcí je možné odvodit i přímý vztah. Pokud do
rovnice (4.26) dosadíme za reakční kvocient Q z rovnice (4.22) a za rovnovážnou konstantu z rovnice (4.21),
dostáváme často používaný vztah
2 2 2
A B A B A B1 1 1
r r
r r r r
r
G G
G G G GRT
RT RT RT
výsledný G
RT
e
r k a a k a a e e k a a e
e

 

 
   
  

 
    
               
 
(4.27)
Záporná hodnota rG poskytne hodnotu exponenciální funkce menší než 1, výraz v závorce bude kladný a
celkový tok látky bude směřovat zleva doprava. Pro kladnou hodnotu rG bude hodnota exponenciální funkce
větší než 1, výraz v závorce bude záporný a výsledný tok látek bude směřovat zprava doleva. Pro nulovou
hodnotu rG , která odpovídá rovnováze, bude hodnota exponenciální funkce právě rovna 1 a výraz v závorce
bude nulový. Tím bude nulový i výsledný tok látek – tok látek zleva doprava je vyrovnán s tokem látek zleva
doprava.
41
5.4 GIBBSOVA AKTIVAČNÍ FUNKCE
Teorie přechodového stavu (transition state theory – TST) nebo také teorie aktivovaného komplexu (activated
complex theory) předpokládá, že při reakci látek podle rovnice (4.1) musí látky nejprve vytvořit aktivovaný
komplex (obr. 5.1b). Pro rychlostní konstantu vzniku aktivovaného komplexu a jeho přechod na látku C je pak
možné použít vztah
*
G
RT
k e



  (4.28)
kde 
*
G°+ je Gibbsova aktivační funkce pro reakci
A + 2 B C
k
 (4.7)
Obdobně je možné pro rychlostní konstantu reakce přeměny látky C na látky A a B
C A + 2 B
k
 (4.8)
psát
*
G
RT
k e



  (4.29)
Čím větší bude hodnota Gibbsovy aktivační funkce v daném směru, tím menší je rychlostní konstanta a tím
pomalejší bude reakce.
Výsledný tok látek pro reakci (4.1)
A + 2 B C
k
k


(4.1)
je roven
2
A B Cvýslednýr k a a k a   (4.4)
Po úpravě dosazením z rovnice (4.25) dostáváme
2 2
A B A Bvýslednýr k a a k a a Q   (4.30)
Pro přechod látek přes aktivovaný komplex pak ve směru zleva doprava platí
* * *
* * * 2
A B A A B BAK AK AK
2 ln ln 2 lnG G G G G RT a G RT a G RT a  
          (4.31)
** * AK
2
A B
ln
a
G G RT
a a

     (4.32)
Pro přechod látek v opačném směru podobně obdržíme
* * *
* * *
C C CAK AK AK
ln lnG G G G RT a G RT a 
       (4.33)
42
** * AK
C
ln
a
G G RT
a

     (4.34)
Pro celkovou změnu zleva doprava pak
* ** * * *AK AK
2
A B C
ln lnr
a a
G G G G RT G RT
a a a
 
              (4.35)
*
*
* * * * CAK
2
A B AK
lnr
a a
G G G G G RT
a a a
 
             (4.36)
* * * *C
2
A B
ln lnr
a
G G G RT G G RT Q
a a
   
              (4.37)
Odtud pak pro reakční kvocient platí
* *
rG G G
RT
Q e
 
   
 (4.38)
Dosazením do rovnice (4.30) z rovnic (4.28), (4.29) a (4.38) obdržíme
* * * *
2 2
A B A B
rG G G G G
RT RT RT
výslednýr e a a e a a e
   
       
 
  (4.39)
a dalšími úpravami dostáváme rovnici
* * * *
2 2
A B A B
rG G G G G
RT RT RT
výslednýr e a a e a a e
  
       
 
  (4.40)
*
2 2
A B A B1 1
r rG G G
RT RT RT
výslednýr e a a e k a a e

  


   
      
   
(4.41)
která je totožná s výrazem, který byl odvozen z čistě termodynamického popisu, viz rovnice (4.27). Za rovnováhy
pak platí
2
A B C 0r r rk a a k a   (4.42)
a b
43
Obr. 5.2 (a) Vznik aktivovaného komplexu z látek A a B je spojen se změnou standardní Gibbsovy aktivační funkce *
G°+ a vznik aktivovaného
komplexu z látky C je spojen se změnou standardní Gibbsovy aktivační funkce *
G°–. (b) Za rovnováhy se změní hodnota Gibbsovy
funkce látek na levé a pravé straně v důsledku změn aktivit reagujících látek tak, až jsou hodnoty Gibbsovy funkce stejné. Hodnota
Gibbsovy reakční funkce – rozdíl hodnot Gibbsovy funkce látek na pravé a levé straně – je pak nulová.
Dosazením za rychlostní konstanty k+ a k– z rovnic (4.28) a(4.29) dostáváme
* *
2
A B C 0
G G
RT RT
r r re a a e a
 
  
 
  (4.43)
* *
2
A B C
G G
RT RT
r r re a a e a
 
  
 
 (4.44)
*
*
C
2
A B
G
RT
r
G
r r RT
a e
a a
e








 (4.45)
* *
G G
RT RT
K e e
 
  

 (4.46)
* *
G G
RT
K e
 
  
 (4.47)
* * *
rG G G
RT RT
K e e
  
   
 
  (4.48)
Tento vztah je totožný s výrazem pro rovnovážnou konstantu, který je možné odvodit z čistě termodynamických
dat, tedy z hodnot Gibbsových funkcí zúčastněných látek. Relativní rozdíl v rychlostech obou procesů je
tedy určen rozdílem standardních hodnot Gibbsovy funkce reagujících látek.
Vznik aktivovaného komplexu je vždy spojen s růstem hodnoty Gibbsovy funkce a tedy poklesem celkové entropie
systému. Reagující látky musí mít dostatečnou zásobu energie tak, aby mohly její určitou část „utratit“
reakční
profil
*G°+
*G°AK*
A + 2B
C
*G°–
G°r
Gibbsovafunkce
reakční koordináta
G°A+2G°B
G°C
RT ln aC
RT ln aA aB
2
reakční
profil
*G+
*G°AK*
A + 2B C
*G–
Gibbsovafunkce
reakční koordináta
G°A+2G°B
G°C
44
v podobě rozptýleného tepla q (sníží se jejich entalpie) a tedy vyprodukovat entropii o velikosti q/T, která bude
tento pokles celkové entropie při vzniku aktivovaného komplexu kompenzovat.
Tak je propojena termodynamika s dynamikou do jednoho celku a vytváří jednolitý obraz principů chemických
přeměn. Změna celkové entropie systému při přeměnách (měřítkem je hodnota Gibbsovy funkce) určuje nejen
to, kterým směrem budou chemické přeměny probíhat – ve směru růstu celkové entropie, tedy poklesu Gibbsovy
funkce –, ale také jak rychle – čím nižší bude hodnota Gibbsovy aktivační funkce, tím bude reakce rychlejší.
Zároveň se také ukazuje, proč některé přeměny neprobíhají vůbec, přestože je standardní hodnota Gibbsovy
funkce pro danou reakci záporná a vzrostla by celková entropie. Pokud nemají reagující látky za daných podmínek
dostatek energie na to, aby mohly uvolnit dostatek tepla na kompenzaci poklesu celkové entropie vznikem
aktivovaného komplexu, pak nebude vůbec aktivovaný komplex vznikat a reakce neprobíhá.
5.5 DALŠÍ INTERPRETACE
Definiční rovnice pro Gibbsovu funkci
G H TS  (4.49)
platí i pro rozdíly termodynamických veličin a tedy i pro vznik aktivovaného komplexu
* * *
G H T S  
       (4.50)
Dosazením do rovnice (4.28) pro rychlostní konstantu dopředné reakce
*
G
RT
k e



  (4.28)
dostáváme
* * * *
H T S H S
RT RT R
k e e e
   
       
 
   (4.51)
kde 
*
H°+ je entalpie aktivace a je 
*
S°+ entropie aktivace.
Entalpie aktivace je teplo, které se zabaví nebo uvolní při vzniku aktivovaného komplexu. Pokud se teplo uvolňuje,
pak je entalpie aktivovaného komplexu nižší, než je entalpie výchozích látek a změna entalpie je záporná.
Čím větší bude teplo, které se uvolní, tím větší bude hodnota první exponenciální funkce, která obsahuje standardní
změnu entalpie reakce a tím vyšší bude rychlostní konstanta reakce. Reakce bude v daném směru probíhat
tím rychleji, čím více tepla se při vzniku aktivovaného komplexu uvolní (zabavení tepla ji bude naopak
zpomalovat).
Entropie aktivace je změna v uspořádání látek, které spolu reagují a vytvářejí aktivovaný komplex. Obvykle
dochází vznikem aktivovaného komplexu k poklesu vlastní entropie reagujících látek – většinou vzniká uspořádanější
přechodný komplex, který se až následně rozpadá nebo mění na produkty. Pokles entropie pak způsobuje,
že hodnota druhého členu klesá a čím je aktivovaný komplex uspořádanější, tím pomalejší bude reakce.
Hnací silou reakce je tedy obvykle teplo, které se při vzniku aktivovaného komplexu uvolňuje a je spojeno
s růstem entropie okolí o velikosti
*
ok
H
S
T


   (4.52)
45
Čím bude větší, tím rychleji bude reakce probíhat. Pokles entropie vlastních látek při vzniku aktivovaného
komplexu naopak rychlost reakcí snižuje. Čím bude aktivovaný komplex uspořádanější (složitější), tím pomaleji
bude reakce probíhat. Tento vliv se označuje jako sterický efekt. Nepříznivá kombinace změny entalpie a entropie
při vzniku aktivovaného komplexu může vést k tomu, že rychlostní konstanta bude mít tak nízkou hodnoty,
že reakce prakticky neprobíhá.
5.6 VZTAH MEZI RYCHLOSTÍ ROZPOUŠTĚNÍ A SRÁŽENÍ
Obecně závisí rychlost rozpouštění a srážení minerálů na velikosti jejich povrchu, rychlostních konstantách a
koncentraci složek v roztoku, které tvoří minerál, případně dalších podmínkách, jako je přítomnost látek, které
určitý proces katalyzují nebo brzdí atd.
Pro jednoduchou rovnici rozpouštění minerálu
C D C D
rozp
sráž
k
y x
x y
k
x y 
 (4.53)
můžeme pro rychlost rozpouštění v nejjednodušším případě psát
rozp rozpr Ak (4.54)
C D
x y
sráž srážr Ak a a (4.55)
kde A je povrch minerálu, k jsou příslušné rychlostní konstanty a a jsou aktivity složek v roztoku. Celkový tok
složek C a D do a nebo z roztoku je dán rozdílem obou rychlostí
C D
x y
výsledný rozp srážr Ak Ak a a  (4.56)
Pokud bude první člen v rovnici (4.56) větší než druhý, bude se minerál rozpouštět. V opačném případě bude
docházet k jeho srážení. Po dosažení vyrovnaných toků rozpouštění a srážení (tedy za rovnováhy) jsou oba toky
vyrovnány a platí
CS DS
x y
rozp srážAk Ak a a (4.57)
kde index S označuje stacionární aktivity, které jsou totožné s rovnovážnými aktivitami. Pro vzájemný vztah
obou rychlostních konstant pak z rovnice (4.57) vyplývá
CS DS
rozp
sráž x y
k
k
a a
 (4.58)
Dosazením do rovnice (4.56) dostáváme
C D
C D
CS DS CS DS
1
x y
rozp x y
výsledný rozp rozpx y x y
k a a
r Ak A a a Ak
a a a a
 
    
 
(4.59)
1výsledný rozp
Q
r Ak
K
 
  
 
(4.60)
46
kde Q je reakční kvocient (součin koncentrací či aktivit složek, které vytvářejí minerál, umocněné na příslušné
stechiometrické koeficienty) a K je rovnovážná konstanta.
Z rovnice (4.60) vyplývá, že rychlost srážení (či rozpouštění) bude tím větší, čím větší bude rychlostní konstanta,
čím větší bude povrch minerálu a čím větší bude přesycení, charakterizované poměrem Q/K. Pokud je Q/K
větší než jedna, pak je aktivita složek v roztoku větší, než odpovídá rovnovážné konstantě a minerál se bude
srážet celkovou rychlostí rvýsledný (výraz v závorce bude záporný). Při hodnotě Q/K menší než jedna bude roztok
vůči minerálu naopak nenasycený a minerál se bude rozpouštět – výsledný tok rvýsledný bude kladný (výraz
v závorce bude kladný). Pokud se bude Q/K rovnat jedné, je koeficient Q roven rovnovážné konstantě K a výsledný
tok bude nulový – minerálu nebude ani přibývat, ani ubývat (výraz v závorce bude roven nule). To
ovšem neznamená, že procesy rozpouštění a srážení ustaly. Probíhají dále, jen se jejich rychlosti vyrovnaly.
Při modelování dynamiky geochemických reakcí a procesů se obvykle místo rovnice (4.56), kterou jsme použili
pro řešení dynamiky rozpouštění nefelínu a srážení gibbsitu, užívá rovnice (4.60), protože umožňuje souběžně
vyhodnocovat změny koncentrací všech složek a jejich vliv na celkovou rychlost a není omezena na jednu složku
a fixované podmínky.
47
6 NUMERICKÉ ŘEŠENÍ DIFERENCIÁLNÍCH ROVNIC
Diferenciální rovnice se v geochemii užívají většinou k vyjádření časového vývoje koncentrace nebo obsahu
látky v systému. Numerické řešení diferenciálních rovnic, známé také jako numerická integrace nebo jako výpočet
integrálu se užívá v případech, kdy není možné diferenciální rovnici řešit analyticky (integrovat ji) nebo je
její integrace obtížná a složitá. Takových rovnic je v geochemii většina.
Pro vysvětlení podstaty numerického řešení diferenciálních rovnic užijeme modelového příkladu vývoje nelineární
funkce
 y f x (5.1)
tj. průběhu závisle proměnné y na nezávisle proměnné x v podobě polynomu
3 2
y ax bx cx d    (5.2)
a její geometrické interpretace. Pokud v bodě [X,Y] sestrojíme k této funkci tečnu, pak bude její směrnice rovna
y/x (obr. 6.1).
Obr. 6.1 Průběh funkce y = f (x) a tečna k této funkci v bodě [X, Y]. Obr. 6.2 Numerická simulace funkce y = f (x) Eulerovou metodou
pro různé velikosti kroků x. Použitý krok je uveden číslem u
křivky.
Pro zvolený úsek x se průběh funkce y velmi výrazně liší od průběhu tečny. Pokud však budeme úsek x
zmenšovat, bude se průběh tečny blížit průběhu funkce v daném úseku. V ideálním případě nekonečně malého
úseku x, který označujeme dx, budou oba průběhy totožné a směrnice tečny je ve kterémkoliv bodě křivky
rovna derivaci funkce y podle x
 
 
f xdy
f x
dx dx
  (5.3)
která má pro zvolený polynom hodnotu
[X,Y]
x
y
0
10
20
30
40
50
60
70
0 2 4 6 8 10 12
y
x
1
0,5
0,1
0
10
20
30
40
50
60
70
0 2 4 6 8 10 12
y
x
48
  2
3 2
dy
f x ax bx c
dx
    (5.4)
To, že se se zmenšováním úseku x blíží průběh nelineární funkce její směrnici v daném bodě, je možné využít
pro numerickou simulaci průběhu nelineární funkce. Pokud bude úsek x dostatečně malý, bude odchylka
zanedbatelná a průběh funkce můžeme v daném úseku nahradit přímkou, jejíž sklon bude roven směrnici tečny.
Přitom je směrnice tečny dána rovnicí (5.4).
Pro výpočet přírůstku hodnoty funkce y z nějakého bodu y0 při zvýšení hodnoty x z x0 na hodnotu x0 + x použijeme
upravenou rovnici (5.4)
2
3 2
y
ax bx c
x

  

(5.5)
 2
3 2y ax bx c x     (5.6)
Celý průběh simulace průběhu funkce (5.2) je pak možné popsat následujícím předpisem (algoritmem):
krok x y y

0 x0 y0 y0 = f´(x0) x
1 x1 = x0 + x y1 = y0 + y0 y1 = f´(x1) x
2 x2 = x1 + x y2 = y1 + y1 y2 = f´(x2) x
… … … …
n–1 xn–1 = xn–2 + x yn–1 = y n–2 + y n–2 yn–1 = f´(xn–1) x
n xn = xn–1 + x yn = y n–1 + y n–1
kde je x krok a x0 a y0 jsou výchozí hodnoty. Tato metoda numerického řešení diferenciálních rovnic se označuje
jako Eulerova metoda. Výsledky simulace modelové funkce pro různé hodnoty kroku jsou uvedeny na obr.
6.2. Dlouhý krok vede k tomu, že se průběh simulované funkce značně liší od jejího skutečného průběhu. Je to
způsobeno tím, že v některých úsecích se mění směrnice funkce rychleji, než aby ji bylo možné bez ztráty přesnosti
nahradit přímkou. Zkracování kroku vede k přibližování simulovaného a skutečného průběhu funkce a je
patrné, že zvolení dostatečně malého kroku povede k dosažení požadované přesnosti. Na druhou stranu vede
zkracování kroku k rychlému růstu nutných výpočtů. Proto je vhodné volit krok co nejdelší a přitom takový, aby
nedocházelo ke ztrátě přesnosti.
Významného zpřesnění souhlasu simulace funkce a jejího skutečného průběhu je možné dosáhnout přesnějším
odhadem směrnice funkce v daném intervalu x. Při Eulerově metodě je použito k výpočtu směrnice v daném
intervalu a tedy změny funkce y pro zvolené x počáteční hodnoty intervalu
49
 1n n ny y f x x
   (5.7)
tedy pro zvolenou funkci
 2
1 3 2n n n ny y ax bx c x      (5.8)
Obr. 6.3 Porovnání numerické simulace průběhu funkce Eulerovou
metodou a metodou středového bodu. Modrá linie je skutečný
průběh funkce. Krok simulace x byl v obou případech roven 1.
Obr. 6.4 Porovnání jednoho kroku numerické simulace Eulerovou
metodou, metodou středu intervalu a Rungeovou-Kuttovou
metodou.
Pokud použijeme k odhadu přírůstku funkce střední hodnotu x v daném intervalu, bude směrnice a tedy i změna
funkce lépe odpovídat změně funkce v daném intervalu
1
2
n n n
x
y y f x x
     
 
(5.9)
a pro zvolenou funkci
2
1 3 2
2 2
n n n n
x x
y y a x b x c x
     
          
     
(5.10)
Jinak zůstává předpis pro simulaci funkce (algoritmus) stejný. Jak je patrné z obr. 6.3, dojde zvolením středu
příslušného intervalu pro výpočet změny y k významnému zpřesnění simulace. Tato metoda bývá také označena
jako Rungeova-Kuttova metoda druhého řádu (RK2). Často se využívají i další metody, jako je RungeovaKuttova
metoda čtvrtého řádu (RK2) a další modifikace uvedeného postupu (viz dále).
Eulerova
metoda
metoda
středového
bodu
0
10
20
30
40
50
60
0 2 4 6 8 10
y
x
Eulerova
metoda
[xn, yn]
metoda
středu
intervalu
RungeovaKuttova
m.
(RK4)
tečna – střed
intervalu
tečna – konec
intervalu
tečna – začátek
intervalu
20
22
24
26
28
30
32
34
36
38
40
4,5 5 5,5 6 6,5 7 7,5 8 8,5 9
y
x
50
6.1 VYUŽITÍ NUMERICKÉHO ŘEŠENÍ DIFERENCIÁLNÍCH ROVNIC V GEOCHEMII
Využití numerického řešení diferenciálních rovnic v geochemii je možné nejlépe demonstrovat na následujících
příkladech. Zatímco jsme schopni velmi dobře určit časovou změnu sledované veličiny, například koncentrace
nebo množství látky v podobě x/t, kde x je sledovaná veličina, a její závislost na podmínkách, mnohem obtížněji
vyjadřujeme závislost sledované veličiny x na čase v podobě závislosti x = f(t), pokud to vůbec dokážeme.
6.2 PŘÍKLAD 1
6.2.1 ZADÁNÍ
Do rezervoáru vtéká konstantní rychlostí látka X, která opět z rezervoáru odtéká. Rezervoár má jednotkový
objem 1 l, výstup látky X je závislý na koncentraci látky v rezervoáru, nemusí to však být lineární závislost. Pro
vývoj koncentrace látky v rezervoáru byly naměřeny hodnoty uvedené v tab. 6.1.
Tab. 6.1 Časový vývoj koncentrace látky X v rezervoáru.
poř. č. t Xnaměřené poř. č. t Xnaměřené
dny mol l
–1
dny mol l
–1
1 0,0 0,000 16 7,5 9,05
2 0,5 0,997 17 8,0 9,22
3 1,0 1,97 18 8,5 9,35
4 1,5 2,91 19 9,0 9,47
5 2,0 3,80 20 9,5 9,56
6 2,5 4,62 21 10,0 9,64
7 3,0 5,37 22 10,5 9,70
8 3,5 6,04 23 11,0 9,76
9 4,0 6,64 24 11,5 9,80
10 4,5 7,16 25 12,0 9,84
11 5,0 7,62 26 12,5 9,87
12 5,5 8,00 27 13,0 9,89
13 6,0 8,34 28 13,5 9,91
14 6,5 8,62 29 14,0 9,93
15 7,0 8,85 30 14,5 9,94
31 15,0 9,95
Určete velikost vstupního toku a rychlostní konstantu výstupního toku.
6.2.2 ŘEŠENÍ
Uvedený systém je možné zobrazit pomocí schématu uvedeného na obr. 6.5.
51
Obr. 6.5 Jednorezervoárový systém s konstantním vstupem a výstupem úměrným koncentraci látky X v rezervoáru.
Podle zadání můžeme předpokládat konstantní vstup do rezervoáru
vstupr A (5.11)
Neznáme však jeho hodnotu. Výstupní tok z rezervoáru je úměrný koncentraci látky X v rezervoáru, nevíme
však, jaká je podoba úměry a neznáme rychlostní konstantu, proto musíme zapsat výstupní tok v podobě rov-
nice
n
výstupr kX (5.12)
Přitom bude pravděpodobné, že exponent u koncentrace látky X bude roven malému celému číslu, případně
jednoduchému zlomku. V úvahu připadá především n = 1 nebo 2.
Obr. 6.6 Vývoj koncentrace látky X v rezervoáru.
Vývoj koncentrace látky X v rezervoáru v závislosti na čase je zobrazen na obr. 6.6. Pro změnu koncentrace
látky X v rezervoáru pak můžeme psát
vstup
A
rezervoár
X
výstup
k×Xn
0
2
4
6
8
10
12
0 5 10 15 20
X(moll–1)
t (dny)
52
n
vstup výstup
dX
r r A kX
dt
    (5.13)
Hodnota dX/dt je směrnicí ke křivce závislosti koncentrace látky X v rezervoáru na čase. Číselné hodnoty v tab.
6.1 umožňují pro jednotlivé kroky vypočítat směrnice v určitém časovém kroku, kdy v prvním přiblížení předpokládáme,
že je časový interval tak krátký, že je možné pro daný úsek nahradit křivku přímkou. Pak bude
platit
1
1
n n
n n
X XX
t t t




 
(5.14)
Tedy například pro pořadové měření 5 to bude
1 15 4
5 4
3,80 2,91 0,89
1,78 mol l den
2,0 1,5 0,5
X XX
t t t
  
   
  
(5.15)
Pro zjištění vstupního toku, exponentu a rychlostní konstanty v rovnici (5.13) vyneseme nejprve závislost
X
A kX
t

 

(5.16)
tedy hodnoty X/t proti X. Pokud by byl výstupní tok látky X z rezervoáru přímo úměrný koncentraci látky X
v rezervoáru, pak bychom měli dostat přímku, jejíž směrnice bude rovna záporné hodnotě rychlostní konstanty
–k a úsek na ose y bude roven hodnotě vstupního toku A. Hodnoty X/t pro jednotlivá měření jsou doplněny
v tab. 6.2 a závislost je vynesena na obr. 6.7a.
Tab. 6.2 Vývoj koncentrací látky X v rezervoáru a vypočítané hodnoty směrnice v jednotlivých měřených krocích včetně hodnot X2
.
poř. č. t Xnaměřené X/t Xnaměřené
2
poř. č. t Xnaměřené X/t Xnaměřené
2
dny mol l
–1
mol l
–1
den
–1
mol
2
l
–2
dny mol l
–1 mol l
–1
den
–1
mol
2
l
–2
1 0,0 0,000 16 7,5 9,05 0,40 81,90
2 0,5 0,997 1,99 0,994 17 8,0 9,22 0,34 85,01
3 1,0 1,97 1,95 3,88 18 8,5 9,35 0,26 87,42
4 1,5 2,91 1,88 8,47 19 9,0 9,47 0,24 89,68
5 2,0 3,80 1,78 14,44 20 9,5 9,56 0,18 91,39
6 2,5 4,62 1,64 21,34 21 10,0 9,64 0,16 92,93
7 3,0 5,37 1,50 28,84 22 10,5 9,70 0,12 94,09
8 3,5 6,04 1,34 36,48 23 11,0 9,76 0,12 95,26
9 4,0 6,64 1,20 44,09 24 11,5 9,80 0,08 96,04
10 4,5 7,16 1,04 51,27 25 12,0 9,84 0,08 96,83
11 5,0 7,62 0,92 58,06 26 12,5 9,87 0,06 97,42
12 5,5 8,00 0,76 64,00 27 13,0 9,89 0,04 97,81
13 6,0 8,34 0,68 69,56 28 13,5 9,91 0,04 98,21
14 6,5 8,62 0,56 74,30 29 14,0 9,93 0,04 98,60
15 7,0 8,85 0,46 78,32 30 14,5 9,94 0,02 98,80
31 15,0 9,95 0,02 99,00
a b
53
Obr. 6.7 Vynesení směrnice závislosti koncentrace látky X na čase proti (a) hodnotě X a proti (b) hodnotě X2
. Zároveň jsou uvedeny
koeficienty rovnice proložení.
Z obr. 6.7a je patrné, že zvolený typ závislosti neodpovídá procesům, které určují koncentraci látky X v rezervoáru.
Proto zkusíme další typ závislosti, který připadá v úvahu
2X
A kX
t

 

(5.17)
Tentokrát budeme vynášet hodnoty X/t proti hodnotám X
2
(obr. 6.7b). Proložení závislosti ve zvolených
souřadnicích poskytuje přímku s koeficienty k = 0,0203 den
–1
a A = 2,06 mol l
–1
den
–1
.
Rovnice (5.17) má pro řešený systém tvar
2
2,06 0,0203
X
X
t

 

(5.18)
Je zřejmé, že výstup látky X z rezervoáru je úměrný druhé mocnině koncentrace látky v rezervoáru. S takto
zjištěnými hodnotami vstupního toku a rychlostní konstanty výstupu je pak možné pro daný systém numericky
simulovat nejrůznější situace. Pro účely simulace je možné vzhledem k přesnosti analytických údajů použít jako
parametry hodnoty k = 0,02 den
–1
a A = 2,00 mol l
–1
den
–1
.
Numerickou simulaci časového vývoje systému provedeme pomocí Eulerovy a Rungeovy-Kuttovy metody čtvrtého
řádu (RK4). Výchozí parametry jsou následující:
parametr hodnota jednotky
X0 = 0 mol l
–1
A = 2 mol l
–1
den
–1
k = 0,02 den
–1
t = 3 den
X/t = –0,232 X + 2,56
R² = 0,9488
0,00
0,50
1,00
1,50
2,00
2,50
0 2 4 6 8 10 12
X/t
X
DX/Dt proložení
X/t = –0,0203X2 + 2,06
R² = 0,9987
0,00
0,50
1,00
1,50
2,00
2,50
0 20 40 60 80 100
X/t
X2
DX/Dt proložení
54
Pro Eulerovu metodu pak bude platit následující předpis:
 2
1 1n n nX X A kX t     (5.19)
Například pro X1 platí
 2
1 0 0X X A kX t    (5.20)
a dosazením
 2
1 00 2 0,02 0 3 6X       (5.21)
Další hodnoty numerické simulace jsou uvedeny v tab. 6.3.
Pro Rungeovu-Kuttovu metodu čtvrtého řádu (RK4) platí předpis:
 1 1 2 3 42 2
6
n n
t
X X k k k k

     (5.22)
kde
 
 
1
2 1
3 2
4 3
,
1 1
,
2 2
1 1
,
2 2
,
n n
n n
n n
n n
k f t X
k f t t X tk
k f t t X tk
k f t t X tk

      
 
      
 
    
(5.23)
a
  2
,
X
f t X A kX
t

   

(5.24)
Veličina tn označuje čas, který byl dosažen postupným přirůstáním času o intervaly t, a podle rovnice (5.24)
můžeme vypočítat pro daný čas směrnici změny koncentrace X.
Přírůstek hodnoty X v rovnici (5.22) je vypočítán jako součin váženého průměru směrnice v daném intervalu a
velikosti intervalu. Význam jednotlivých koeficientů, definovaných vztahy (5.23) je následující:
k1 – je směrnice funkce na počátku intervalu, tedy v bodě tn, tato směrnice je totožná se směrnicí Eulerovy
metody
k2 – je směrnice funkce uprostřed intervalu, tedy v bodě tn + ½ t, přitom je použito směrnice k1 pro výpočet
hodnoty X v tomto bodě
k3 – je směrnice funkce uprostřed intervalu, tedy opět v bodě tn + ½ t, nyní je však použito směrnice k2 pro
výpočet hodnoty X
k4 – je směrnice funkce na konci intervalu, tedy v bodě tn + t, pro výpočet hodnoty X je použito směrnice k3.
55
Při průměrování těchto čtyř směrnic je významně větší váha přisouzena směrnicím uprostřed intervalu
 1 2 3 4
1
2 2
6
vážená směrnice k k k k    (5.25)
 1 2 3 42 2
6
t
vážený přírůstek k k k k

    (5.26)
Poznámka
Metoda středu intervalu bývá také označována jako Rungeova-Kuttova metoda druhého řádu (RK2). Platí pro ni
předpis
1 2n nX X tk    (5.27)
kde
 1
2 1
,
1 1
,
2 2
n n
n n
k t f t X
k f t t X k
  
     
 
(5.28)
Pro porovnání Eulerovy a Rungeovy-Kuttovy metody čtvrtého řádu (RK4), byl záměně zvolen dlouhý krok 3 dny.
Předpis pro Rungeovu-Kuttovu metodu vypadá složitě, ale dá se snadno implementovat v běžném tabulkovém
kalkulátoru. Pro kontrolu je první krok RK4 uveden dosazením do vzorců
     
 
 
2 2
1 1 1 1
2
2
2 1 1 1 1
2
2
3 1 1 2 1
4 1
, 2 0,02 0 2,00
1 1 1
, 3 2,00 2 0,02 0 3 1,82
2 2 2
1 1 1
, 3 1,82 2 0,02 0 2,73 1,85
2 2 2
k f t X A kX
k f t t X tk A k X
k f t t X tk A k X
k f t
      
                              
                              
      
2 2
1 3 1, 3 1,82 2 0,02 0 5,55 1,38t X tk A k X               
   
(5.29)
Změna koncentrace látky X je pak v kroku 1 (čase t = 0 dnů) rovna
 1 2,00 2 1,82 2 1,85 1,38 5,36
6
t
X

        (5.30)
a hodnota X bude v čase t2 = t1 + t = 0 + 3 = 3 rovna
2 1 1 0 5,36 5,36X X X      (5.31)
Hodnoty koeficientů metody RK4, přírůstky X a vývoj koncentrací v čase jsou uvedeny v tab. 6.3 a na obr. 6.8.
56
Tab. 6.3 Numerická simulace vývoje koncentrace látky X v rezervoáru Eulerovou a Rungeovou-Kuttovou metodou čtvrtého řádu.
Časová krok byl zvolen 3 dny.
metoda Euler met. RK4
poř. č. t Xnum. simulace X Xnum. simulace X k1 k2 k3 k4
dny mol l
–1
mol l
–1
mol l
–1
mol l
–1
1 0 0,00 6,00 0 5,36 2,00 1,82 1,85 1,38
2 3 6,00 3,84 5,36 2,94 1,42 0,88 1,11 0,49
3 6 9,84 0,19 8,30 1,13 0,62 0,29 0,47 0,11
4 9 10,03 -0,04 9,44 0,38 0,22 0,09 0,17 0,03
5 12 9,99 0,01 9,82 0,12 0,07 0,03 0,05 0,01
6 15 10,00 9,94
Obr. 6.8 Vývoj koncentrace látky X v rezervoáru. Modrou linií jsou uvedeny naměřené hodnoty, hnědé body vyznačují numerickou
simulaci (numerické řešení) Eulerovou metodou, zelené kroužky ukazují jednotlivé kroky metody RK4. Záměrně byl zvolen dlouhý krok,
aby byl vidět rozdíl v obou metodách.
Poznámka
Výhody numerické simulace řešeného problému jsou patrné také při jeho porovnání s analytickým řešení uvedené
úlohy. Z naměřených hodnot jsme zjistili, že vývoj koncentrace látky X v rezervoáru je dán vztahem
2X
A kX
t

 

(5.17)
V diferenciální podobě jej můžeme přepsat do tvaru
2dX
A kX
dt
  (5.32)
0
2
4
6
8
10
12
0 5 10 15 20
X(moll–1)
t (dny)
Xnaměřené Euler met. RK4
57
Tuto diferenciální rovnici ještě dokážeme integrovat
2
dX
dt
A kX


(5.33)
0
2 0
X t
X
dX
dt
A kX

  (5.34)
0
2 0
X t
X
dX
dt
A kX

  (5.35)
 
0
0
1
ln
2
X
t
X
A k X
t
kA A k X
 
  
  
(5.36)
0
0
1
ln ln
2
A k XA k X
t
kA A k X A k X
 
   
   
(5.37)
0
0
ln 2
A k X
A k X
kAt
A k X
A k X


 


(5.38)
20
0
kAtA k XA k X
e
A k X A k X


 
(5.39)
Označením členu na pravé straně rovnice (5.39)
20
0
kAtA k X
e B
A k X



(5.40)
Dosazením a další úpravou
 A k X A k X B   (5.41)
A k X AB k XB   (5.42)
A AB k XB k X   (5.43)
   1 1A B k B X   (5.44)
dostáváme pro vývoj koncentrace látky X v rezervoáru výraz
58
 
 
1
1
A B
X
k B



(5.45)
Pro dostatečně dlouhý čas t se hodnota členu B blíží k nule B0 a hodnota X se blíží k
 
 
1 0
1 0
A A
X
kk

 

(5.46)
Stacionární stav koncentrace látky X v rezervoáru je dán podmínkou
2
S 0
dX
A kX
dt
   (5.47)
a odtud pak
S
A
X
k
 (5.48)
Měli jsme štěstí, že byla diferenciální podoba vývoje koncentrace látky X v rezervoáru integrovatelná (to je
možné jen v nejjednodušších případech geochemických systémů), přesto je analytické vyjádření rychlosti procesu,
který určuje koncentraci látky X v rezervoáru značně složité.
Numerické řešení dynamiky geochemických procesů je přímočaré. Známé vyjádření rychlosti procesů v diferenciální
podobě vede po vložení příslušných vztahů do tabulkového kalkulátoru k přímé simulaci vývoje. Kromě
toho existuje celá řada specializovaných programových kódů, které jsou schopny po vložení diferenciální
podoby rychlostních rovnic numericky řešit i vysoce nelineární systémy (tzv. stiff rovnice).
59
6.3 PŘÍKLAD 2
6.3.1 ZADÁNÍ
Koncentrace rozpuštěných látek v důlní vodě je v principu na jedné straně určena rychlostí jejich uvolňování
z horninového prostředí (pórová voda, sekundární produkty) a na druhé straně jejich vyvazováním do dalších
produktů, odváděním důlních vod z dolů a jejich čištěním. Principiální schéma je uvedeno na obr. 6.9.
Obr. 6.9 Schéma procesů, které určují koncentraci složek v důlních vodách. Složka je do rezervoáru důlních vod DV uvolňována ze zdrojového
rezervoáru Z určitou rychlostí. Z rezervoáru důlních vod je naopak odčerpávána vznikem dalších minerálů či odváděním důlních
vod na povrch a jejich čištěním.
Vývoj koncentrace složky v důlních vodách může být často popsán diferenciálními rovnicemi
1
dx
k x
dt
  (5.49)
1 2
dy
k x k y
dt
  (5.50)
kde x je koncentrace sledované složky ve zdroji, y je koncentrace této složky v důlní vodě, k1 je rychlostní konstanta
uvolňování složky do důlní vody a k2 je rychlostní konstanta jejího odvádění z důlní vody. Simulujte vývoj
koncentrace složky v důlní vodě, parametry modelu jsou uvedeny v tab. 6.4.
1Zdroj
Z
pórové
vody,
sekundární
produkty –
Fe(OH)3(s),
…
2Důlnívoda
DV
rozpuštěné
složky –
Fe2+, Mn2+,
UO2
2+, SO4
2– ,
…
3Produkty
P
minerály
UO2(s),
FeS2, …
čištění
…
60
Tab. 6.4 Vstupní hodnoty pro simulaci
vývoje koncentrace složky v důlní vodě.
parametr hodnota jednotky
x0 = 190 mg l
–1
y0 = 0 mg l
–1
k1 = 1 den
–1
k2 = 0,15 den
–1
Význam vstupních parametrů je následující: x0 je počáteční koncentrace složky ve zdroji, y0 je počáteční koncentrace
složky v důlní vodě. Pro simulaci použijte Eulerovu metodu s vhodným krokem a porovnejte dosažené
výsledky se simulací metodou RK2.
6.3.2 ŘEŠENÍ
Pro simulaci časového vývoje koncentrace složky v důlní vodě Eulerovou metodou použijeme následující před-
pis:
1 1 1n n nx x k x t    (5.51)
 1 1 1 2 1n n n ny y k x k y t      (5.52)
Pro počáteční krok jsou při časovém kroku t = 0,5 dne dosazeny a vypočítány hodnoty
0,5 0 1 0 190 1 190 0,5 95x x k x t        (5.53)
   0,5 0 1 0 2 0 0 1 190 0,15 0 0,5 0,5y y k x k y t           (5.54)
kdy pro jednoduchost index u hodnot x a y vyjadřuje dosažený čas. Další hodnoty jsou uvedeny pro časový
interval do šesti dnů v tab. 6.5.
Přestože nás zajímá pouze koncentrace složky v důlní vodě y, musíme při numerické simulaci počítat i okamžitou
koncentraci složky ve zdroji x, protože ji potřebujeme pro výpočet koncentrace složky v důlní vodě. Porovnání
výsledků Eulerovy metody pro časové kroky t = 1, 0,5 a 0,1 dne je uvedeno na obr. 6.10a.
61
Tab. 6.5 Vývoj koncentrací složek ve zdroji a v důlní vodě v jednotlivých
krocích Eulerovy metody pro časový krok t = 0,5 dne.
t x y x y
0,0 190,0 0 -95,0 95,0
0,5 95,0 95,0 -47,5 40,4
1,0 47,5 135,4 -23,8 13,6
1,5 23,8 149,0 -11,9 0,7
2,0 11,9 149,7 -5,9 -5,3
2,5 5,9 144,4 -3,0 -7,9
3,0 3,0 136,5 -1,5 -8,8
3,5 1,5 127,8 -0,7 -8,8
4,0 0,7 118,9 -0,4 -8,5
4,5 0,4 110,4 -0,2 -8,1
5,0 0,2 102,3 -0,1 -7,6
5,5 0,1 94,7 0,0 -7,1
6,0 0,0 87,7 0,0 -6,6
a b
Obr. 6.10 Numerická simulace vývojek koncentrace složky v důlní vodě. (a) Porovnání výsledků pro časové kroky t = 1, 0,5 a 0,1 dne.
Skutečný vývoj pro zadané parametry je vynesen modrou linií. (b) Vývoj koncentrace složky v důlní vodě simulovaný Eulerovou metodou
a metodou RK2 s časovým krokem t = 0,5 dne.
I v tomto případě je systém dvou diferenciálních rovnic tak jednoduchý, že je integrovatelný a umožňuje porovnat
výsledky numerické simulace s analytickým řešením. Dostatečně malý krok umožňuje dosáhnout poža-
1
0,5
0,1
0
20
40
60
80
100
120
140
160
180
200
0 1 2 3 4 5 6
koncentrace(mgl–1)
t (dny)
Eulerova
metoda
RK2
0
20
40
60
80
100
120
140
160
0 1 2 3 4 5 6
koncentrace(mgl–1)
t (dny)
62
dovaného výsledku. Při numerické simulaci, kdy neznáme analytické řešení, zkracujeme krok tak dlouho, až je
mezi dvěma po sobě jdoucími simulacemi rozdíl menší, než je požadovaná přesnost simulace.
Pro principiální porovnání Eulerovy metody a metod Rungeho-Kutty použijeme metodu RK2, abychom ilustrovali
princip a zároveň mohli posoudit výhody RK metod a zvolíme relativně velký krok t = 0,5 dne. Pro metodu
RK2 platí obecný předpis
 1
2 1
,
1 1
,
2 2
n n
n n
k t f t X
k f t t X k
  
     
 
(5.28)
1 2n nX X tk    (5.27)
který bude při přepisu pro řešený systém vypadat následovně
 
 
1 1
1 1 2
2 1 1
2 1 1 2 1
1
2
1 1
2 2
x n
y n n
x n x
y n x n y
k t k x
k t k x k y
k k x k t
k k x k k y k t
   
   
 
    
 
    
        
    
(5.28)
1 2
1 2
n n x
n n y
x x k
y y k


 
 
(5.27)
Indexy x a y jsou u koeficientů metody RK2 použity pro odlišení proměnných, kterých se týkají a pro odlišení od
rychlostních konstant. Pro první krok simulace jsou do předchozích vzorců dosazeny konkrétní hodnoty, další
kroky jsou uvedeny v tab. 6.6.
 
 
1
1
2
2
0,5 1 190 95,0
0,5 1 190 0,15 0 95,0
1
1 190 95 0,5 71,2
2
1 1
1 190 95 0,15 0 95 0,5 67,7
2 2
x
y
x
y
k
k
k
k
    
     
 
        
 
    
          
    
(5.28)
0,5 0 2
0,5 0 2
190 71,2 118,8
0 67,7
x
y
x x k
y y k
    
   
(5.27)
Tab. 6.6 Simulace vývoje koncentrace složky v důlní vodě metodou RK2 s časovým krokem t = 0,5
dne.
t x y k1x k1y k2x k2y
den mg l
–1
mg l
–1
mg l
–1
mg l
–1
mg l
–1
mg l
–1
63
0,0 190 0,0 -95,0 95,0 -71,3 67,7
0,5 118,8 67,7 -59,4 54,3 -44,5 37,4
1,0 74,2 105,1 -37,1 29,2 -27,8 18,9
1,5 46,4 124,0 -23,2 13,9 -17,4 7,6
2,0 29,0 131,5 -14,5 4,6 -10,9 0,8
2,5 18,1 132,4 -9,1 -0,9 -6,8 -3,1
3,0 11,3 129,3 -5,7 -4,0 -4,2 -5,3
3,5 7,1 124,0 -3,5 -5,8 -2,7 -6,4
4,0 4,4 117,5 -2,2 -6,6 -1,7 -6,9
4,5 2,8 110,6 -1,4 -6,9 -1,0 -7,0
5,0 1,7 103,6 -0,9 -6,9 -0,6 -6,9
5,5 1,1 96,8 -0,5 -6,7 -0,4 -6,6
6,0 0,7 90,2 -0,3 -6,4 -0,3 -6,3
Porovnání simulace vývoje koncentrace složky v důlní vodě Eulerovou metodou a metodou RK2 ukazuje, že
Rungeova-Kuttova metoda druhého řádu dává i při relativně dlouhém časovém kroku výrazně lepší výsledky
(obr. 6.10b).
64
7 SEZNAMY
7.1 SEZNAM OBRÁZKŮ
Obr. 2.1 Koncepce systémů a rezervoárů. (a) Celý zemský systém může být rozdělen na čtyři velké rezervoáry,
mezi kterými se vyměňují látky – atmosféra, litosféra, hydrosféra a litosféra. (b) Celý systém hydrosféry je
možné rozdělit na jednotlivé rezervoáry, které obsahují vodu a znázornit vzájemné toky mezi nimi. (c) Každý
rezervoár je charakterizován hmotným obsahem látky a toky jsou vyjádřeny jako přenesené množství látky za
časovou jednotku. (d) Globální biogeochemický cyklus uhlíku v podobě „krabičkového“ modelu. Krabičky
představují jednotlivé rezervoáry uhlíku (různě vázaný uhlík), šipky znázorňují toky uhlíku mezi rezervoáry.
Čísla uvnitř krabiček značí odhadované množství uhlíku v Gt (10
15
g C). Čísla u jednotlivých toků jsou
odhadované toky uhlíku mezi rezervoáry v Gt C za rok. Jednotlivé odhady se velmi liší a dosud neexistuje
všeobecně přijímaný model. Odhadované antropogenní ovlivnění globálního cyklu je znázorněno v levé části
diagramu a představuje spalování fosilních paliv a vypalování pralesů. Pro konstrukci modelu globálního
uhlíkového cyklu byla modifikována data Siegenthalera a Sarmienta (1993) a Kwona a Schnoora (1995). .................6
Obr. 2.2 Konstrukce box modelu pro studovaný systém. (a) Cirkulace vody na studovaném ostrově je určována
několika základními toky. (b) Tyto toky zajišťují přenos vody mezi čtyřmi rezervoáry. (c) Identifikovaných toků a
rezervoárů do box modelu. (d) Konečný box model cirkulace vody na ostrově s vyznačení obsahu rezervoárů a
velikosti toků. ................................................................................................................................................................7
Obr. 2.3 Jednorezervoárový systém s konstantním vstupem a výstupem úměrným množství látky X
v rezervoáru...................................................................................................................................................................9
Obr. 2.4 Vývoj množství látky v jednorezervoárovém systému s konstantním vstupem a výstupem přímo
úměrným množství látky v rezervoáru. Vstupní parametry: rychlost vstupu látky X do rezervoáru A = 2 g hod
–1
,
rychlostní konstanta výstupu látky z rezervoáru k = 0,5 hod
–1
, počáteční množství látky X v rezervoáru X0 = 18 g.
Doba odezvy 1 je definována jako pokles odchylky od stacionárního stavu na 10 % původní hodnoty, doba
odezvy 2 je definována jako pokles odchylky na 10 % hodnoty stacionárního stavu.................................................11
Obr. 2.5 Dvourezervoárový systém se vzájemným přenosem látky mezi rezervoáry. Toky látky mezi rezervoáry
jsou přímo úměrné jejich množství. ............................................................................................................................14
Obr. 2.6 Vývoj množství látky ve dvourezerovárovém systému, kde jsou výstupy z rezervoárů přímo úměrné
množství látky v rezervoárech. Vstupní podmínky: počáteční množství látky v rezervoáru A tj. mA0 = 3 mol l
–1
,
počáteční množství látky v rezervoáru B tj. mB0 = 7 mol l
–1
, rychlostní konstanta přestupu látky z rezervoáru A
do rezervoáru B tj. kAB = 0,1 hod
–1
, rychlostní konstanta přestupu látky z rezervoáru B do rezervoáru A tj.
kBA = 0,4 hod
–1
. Doba odezvy 1 je definována jako pokles odchylek od stacionárního stavu na 10 % původních
hodnot a je pro oba rezervoáry stejná. Doba odezvy 2 je definována jako pokles odchylek na 10 % hodnot
stacionárních stavů a pro jednotlivé rezervoáry se liší................................................................................................16
Obr. 2.7 Diagram třírezervoárového průtočného systému se zpětnými toky mezirezervoáry. ..................................18
Obr. 3.1 Jednorezervoárový systém s konstantním vstupem a výstupem úměrným koncentraci (množství) látky
X v rezervoáru..............................................................................................................................................................21
Obr. 3.2 Časový vývoj množství látky X v rezervoáru. (a) Modrá linie ukazuje časový vývoj vypočítaný podle
rovnice (2.25), červené kroužky numerickou simulaci vývoje Eulerovou metodou s časovým krokem 0,5. (b)
65
Vyznačení stacionárního stavu a časů odezvy definovaných jako 10 % původního vychýlení (1) a jako 10 %
odchylky z hodnoty stacionárního stavu (2)................................................................................................................22
Obr. 3.3 Dvourezervoárový systém se vzájemným přenosem látky mezi rezervoáry. Toky látky mezi rezervoáry
jsou přímo úměrné jejich koncentracím......................................................................................................................25
Obr. 3.4 Vývoj koncentrací ve dvourezervoárovém systému s obousměrným tokem přímo úměrným
koncentrace látky v rezervoáru. Vývoj koncentrace v rezervoáru A je vyznačen modrou linií, vývoj koncentrace
látky v rezervoáru B je znázorněn zelenou linií. (a) Numerická simulace vývoje s časovým korkem 0,5 hod. je
vyznačena kroužky. (b) Doby odezvy v rezervoárech – doba odezvy 1 je čas, za který poklesne odchylka od
stacionárního stavu na 10 % původní hodnoty (je stejná pro oba rezervoáry, doba odezvy 2 je čas, za který se
koncentrace v rezervoárech liší o 10 % hodnoty stacionárního stavu (je pro rezervoáry A a B odlišná). ...................27
Obr. 3.5 Třírezervoárový model rozpouštění nefelínu a srážení gibbsitu. Vyznačeny jsou toky, které určují
koncentraci rozpuštěného hliníku ve vodě..................................................................................................................30
Obr. 3.6 Simulace vývoje koncentrace rozpuštěného hliníku ve vodě podle analytického řešení (linie) a
Eulerovou metodou (body). (a) Rozpouštění nefelínu: cS-nef – stacionární koncentrace hliníku, která je zároveň
rovnovážnou koncentrací. (b) Rozpouštění gibbsitu: cS-nef – stacionární koncentrace hliníku, která je zároveň
rovnovážnou koncentrací. ...........................................................................................................................................32
Obr. 3.7 Vývoj koncentrace rozpuštěného hliníku při rozpouštění nefelínu a srážení gibbsitu. (a) Pro zadané
podmínky. (b) Pokud se změní poměr rychlostních konstant, pak se stacionární stav posouvá směrem ke stavu
nasycení vůči minerálu, který má větší rychlostní konstantu rozpouštění. Podobný vliv má aktuální povrch
minerálů. Pro řešený příklad byl poměr rychlostních konstant knefelín rozp/kgibbsit rozp = 33,3. Pro srovnání jsou
uvedeny stacionární stavy při dalších poměrech rychlostních konstant.....................................................................35
Obr. 4.1 (a) Látka A musí při chemické reakci přeměny na látku B překonat energetickou bariéru o velikosti Ea A,
která se označuje jako aktivační energie. Rychlostní konstanta této reakce je podle Arrheniovy rovnice úměrná
velikosti aktivační energie. Při opačné reakci musí překonat aktivační energii Ea B, která je větší než Ea A, a proto
bude rychlostní konstanta této reakce menší. (b) Vznik aktivovaného komplexu je spojen s růstem Gibbsovy
funkce 
*
G°+ a 
*
G°– (poklesem celkové entropie). Rychlostní konstanta přeměny látek A a B na C bude větší a
proto bude při přeměně preferován vznik látky C. Stejný výsledek poskytne i čistě termodynamické posouzení
rovnováhy mezi látkami A a B na jedné straně a látkou C na druhé straně, tedy čistě termodynamický popis
pomocí Gibbsovy reakční funkce.................................................................................................................................37
Obr. 4.2 (a) Vznik aktivovaného komplexu z látek A a B je spojen se změnou standardní Gibbsovy aktivační
funkce 
*
G°+ a vznik aktivovaného komplexu z látky je je spojen se změnou standardní Gibbsovy aktivační
funkce 
*
G°–. (b) Za rovnováhy se změní hodnota Gibbsovy funkce látek na levé a pravé straně v důsledku
změn aktivit reagujících látek tak, až jsou stejné. Hodnota Gibbsovy reakční funkce – rozdíl hodnot Gibbsovy
funkce látek na pravé a levé straně – je pak nulová....................................................................................................43
Obr. 5.1 Průběh funkce y = f (x) a tečna k této funkci v bodě [X, Y]. ...........................................................................47
Obr. 5.2 Numerická simulace funkce y = f (x) Eulerovou metodou pro různé velikosti kroků x. Použitý krok je
uveden číslem u křivky. ...............................................................................................................................................47
Obr. 5.3 Porovnání numerické simulace průběhu funkce Eulerovou metodou a metodou středového bodu.
Modrá linie je skutečný průběh funkce. Krok simulace x byl v obou případech roven 1..........................................49
Obr. 5.4 Porovnání jednoho kroku numerické simulace Eulerovou metodou, metodou středu intervalu a
Rungeovou-Kuttovou metodou...................................................................................................................................49
66
Obr. 5.5 Jednorezervoárový systém s konstantním vstupem a výstupem úměrným koncentraci látky X v
rezervoáru. ..................................................................................................................................................................51
Obr. 5.6 Vývoj koncentrace látky X v rezervoáru.........................................................................................................51
Obr. 5.7 Vynesení směrnice závislosti koncentrace látky X na čase proti (a) hodnotě X a proti (b) hodnotě X
2
.
Zároveň jsou uvedeny koeficienty rovnice proložení. .................................................................................................53
Obr. 5.8 Vývoj koncentrace látky X v rezervoáru. Modrou linií jsou uvedeny naměřené hodnoty, hnědé body
vyznačují numerickou simulaci (numerické řešení) Eulerovou metodou, zelené kroužky ukazují jednotlivé kroky
metody RK4. Záměrně byl zvolen dlouhý krok, aby byl vidět rozdíl v obou metodách...............................................56
Obr. 5.9 Schéma procesů, které určují koncentraci složek v důlních vodách. Složka je do rezervoáru důlních vod
DV uvolňována ze zdrojového rezervoáru Z určitou rychlostí. Z rezervoáru důlních vod je naopak odčerpávána
vznikem dalších minerálů či odváděním důlních vod na povrch a jejich čištěním.......................................................59
Obr. 5.10 Numerická simulace vývojek koncentrace složky v důlní vodě. (a) Porovnání výsledků pro časové
kroky t = 1, 0,5 a 0,1 dne. Skutečný vývoj pro zadané parametry je vynesen modrou linií. (b) Vývoj
koncentrace složky v důlní vodě simulovaný Eulerovou metodou a metodou RK2 s časovým krokem t = 0,5
dne...............................................................................................................................................................................61
67
7.2 SEZNAM TABULEK
Tab. 3.1 Vstupní parametry pro model jednorezervoárového systému s konstantním vstupem a výstupem
přímo úměrným množství látky v rezervoáru. Význam jednotlivých parametrů: A – konstantní vstupní tok do
rezervoáru, k – rychlostní konstanta výstupního toku z rezervoáru, X0 – počáteční koncentrace látky X
v rezervoáru.................................................................................................................................................................21
Tab. 3.2 Numerická simulace vývoje množství látky X v rezervoáru pro výchozí podmínky uvedené v tab. 3.1.
Zvolen byl časový krok 0,5 hod....................................................................................................................................22
Tab. 3.3 Vstupní parametry pro model dvourezervoárového systému s obousměrnými toky přímo úměrnými
koncentracím látky v rezervoárech. ............................................................................................................................25
Tab. 3.4 Výsledky numerické simulace vývoje koncentrací ve dourezervoárovém systému s obousměrným
tokem...........................................................................................................................................................................27
Tab. 5.1 Časový vývoj koncentrace látky X v rezervoáru.............................................................................................50
Tab. 5.2 Vývoj koncentrací látky X v rezervoáru a vypočítané hodnoty směrnice v jednotlivých měřených
krocích včetně hodnot X
2
.............................................................................................................................................52
Tab. 5.3 Numerická simulace vývoje koncentrace látky X v rezervoáru Eulerovou a Rungeovou-Kuttovou
metodou čtvrtého řádu. Časová krok byl zvolen 3 dny. ..............................................................................................56
Tab. 5.4 Vstupní hodnoty pro simulaci vývoje koncentrace složky v důlní vodě.........................................................60
Tab. 5.5 Vývoj koncentrací složek ve zdroji a v důlní vodě v jednotlivých krocích Eulerovy metody pro časový
krok t = 0,5 dne. ........................................................................................................................................................61
Tab. 5.6 Simulace vývoje koncentrace složky v důlní vodě metodou RK2 s časovým krokem t = 0,5 dne................62