ČESKÁ ZEMĚDĚLSKÁ UNIVERZITA V PRAZE FYZIKÁLNĚ-CHEMICKÉ ASPEKTY PROCESŮ V PROSTŘEDÍ

ČESKÁ ZEMĚDĚLSKÁ UNIVERZITA V PRAZE FAKULTA ŽIVOTNÍHO PROSTŘEDÍ

FYZIKÁLNĚ-CHEMICKÉ ASPEKTY PROCESŮ V PROSTŘEDÍ Marek Vach

2007 0

Předmluva........................................................................................................................... 2 1.

Základní teze rovnovážné termodynamiky ................................................................... 3 1.1

Základní pojmy ...................................................................................................... 3

1.2

Ideální plyn ............................................................................................................. 5

1.3

První věta termodynamická .................................................................................... 6

1.4

Druhá věta termodynamická ................................................................................ 10

1.5

Gibbsova energie .................................................................................................. 15

1.6

Termodynamické potenciály ................................................................................ 17

1.6

Třetí věta termodynamická................................................................................... 18

2.

Rovnovážná konstanta chemické reakce a její relace s Gibbsovou energií ................. 19

3.

Vybrané aspekty chování iontů v roztocích ................................................................. 23 3.1.

Aktivta a aktivitní koeficient ................................................................................ 23

3.2

Debye-Hückelova teorie iontové atmosféry ......................................................... 25

4. Oxidačně – redukční procesy, vztah Gibbsovy energie a oxidačně-redukčního potenciálu ............................................................................................................................. 29 5.

Kinetika fyzikálně-chemických procesů v prostředí .................................................... 31 5.1

Rychlost, rozsah a řád procesů ............................................................................. 31

5.2

Procesy nultého a prvního řádu ............................................................................ 32

5.3

Procesy druhého řádu ........................................................................................... 34

5.4

Procesy n-tého řádu .............................................................................................. 35

5.5

Simultánní procesy ............................................................................................... 36

5.6

Příklad systému kinetických procesů ................................................................... 41

6.

Nelineární dynamické systémy – úvodní poznámky.................................................... 44

7.

Příklad řešení kinetiky iontové výměny při průtoku dráhou s aktivním povrchem ..... 48

Dodatek A – Základní numerické metody řešení problémů s počátečními podmínkami, tj. obyčejných diferenciálních rovnic a jejich soustav .............................................................. 54 Dodatek B – Numerické řešení parciálních diferenciálních rovnic – aproximace Galerkinovou metodou ......................................................................................................... 62 Použitá literatura .................................................................................................................. 68

1

Předmluva Předkládaný soubor několika vybraných témat dotýkajících se problematiky fyzikální chemie, fyziky a matematiky není koncipován jako systematický jednooborový výklad. Na to by daný rozsah 67 stran zajisté nestačil. Tento učební text se svojí redukovanou formou pouze pokouší nenásilně naznačit čtenáři některé možnosti náhledu na procesy v prostředí perspektivou fyzikálního obrazu a matematického modelování. Výchozím tématem jsou základní teze rovnovážné termodynamiky. Neopomenutí této klasické disciplíny je motivováno její filozofickou hloubkou, která si zasluhuje velkého respektu. Vždyť např. struktura druhé věty termodynamické může někdy nabídnout zákoutí, jež nemusela být při jejím prvotním studiu čtenářem probádána. První statí, která se přímo obrací ke konstrukci matematického modelu, je Debye-Hückelova teorie iontové atmosféry. Toto téma je samozřejmě do jisté míry vytrženo z rozsáhlé mozaiky tvořící problematiku fyzikální chemie roztoků, nicméně nabízí ucelený obraz modelového řešení daného typu sféricky symetrického problému. Významná pozornost je v textu věnována dynamickým procesům popsaným obyčejnými diferenciálními rovnicemi. Výklad se odvíjí od základních jednoduše analyticky řešitelných kinetických úloh až po některé klasické příklady nelineárních dynamických systémů. Uvedené principy konstrukce dynamických modelů jsou obecně použitelné pro časově závislé děje libovolného typu, byť byly primárně rozpracovány převážně pro chemické procesy. Závěrem je prezentován příklad modelového řešení specifického problému kinetiky iontové výměny s 1D prostorovou dimenzí. Tato stať snad může jistou měrou rozšířit úhel pohledu na modelování procesů daného typu, i když není „učebnicového původu“. Důležitou částí textu jsou přiřazené dodatky, jejichž obsahem je stručný popis vybraných základních metod numerického řešení obyčejných a parciálních diferenciálních rovnic. Tento učební text se vztahuje k předmětu stejného názvu, jenž je vyučován pro studenty 1. ročníku magisterského studia oboru Environmentální modelování. Forma výkladu je volena se snahou o srozumitelnost i pro čtenáře, který nemusí mít vždy do důsledků vystavěn ucelený systém vyšší matematiky. Mnohde se samozřejmě může jevit míra této „empatie“ jako neadekvátní. V textu jsou např. vypisovány i některé úpravy (postupy substituce v diferenciálních rovnicích ap.), které standardní učebnice nechávají na čtenáři. Uváděná vysvětlení a naznačení kontextu mohou vést čtenáře často spíše k intuitivnímu než ke korektnímu a matematickým formalismem podloženému zažití problému. Zasvěcenější čtenář nechť projeví k občasné naivitě textu shovívavost.

Autor

2

1. Základní teze rovnovážné termodynamiky Rovnovážná termodynamika je konzistentní disciplínou, v mnoha aspektech důsledně filosoficky konstruovanou.

Základní pojmy

1.1

Termodynamika nezkoumá mikrostrukturu, rozlišuje obecně vymezenou soustavu – systém determinovaný (jako celek) stavovými parametry. Mezi systémem a jeho okolím může probíhat výměna hmoty a energie. Systém (soustava) -

izolovaný - nevyměňuje hmotu ani energii s okolím

-

uzavřený - nevyměňuje hmotu, vyměňuje energii

-

otevřený - vyměňuje hmotu i energii (např. živý organizmus)

-

otevřený s ustáleným tokem

Formy energie vyměňované mezi soustavou a okolím -

teplo Q - výměna probíhá v důsledku teplotního rozdílu

-

práce W - výměna probíhá působením sil. Tj. práce může být objemová reprezentovaná změno objemu media tvořícího soustavu a mechanická.

Vždy je odčítán počáteční stav soustavy od konečného (obecně platné - v matematice viz. určitý integrál). Pak platí znaménková konvence: + energie dodaná do systému - energie odebraná ze systému W > 0 systém odebírá práci W < 0 systém koná práci Q > 0 systém odebírá teplo

- endotermický děj

Q < 0 systém produkuje teplo

- exotermický děj

Q = 0 adiabatický děj Fáze, složky, skupenské přeměny soustavy Fáze je homogenní oblast soustavy - vykazující pouze spojité změny vlastností v prostoru. Soustava je tedy homogenní nebo heterogenní - složená z více fází. Každá fáze může být složena s více chemických individuí – složek. Jako složky jsou brána pouze chemická individua, jejichž koncentrace lze nezávisle měnit. Tj. nikoli látky vznikající z nezávislých složek chemickou reakcí, resp. látky koncentračně závislé na složkách prostřednictvím chemických rovnováh.

3

Skupenská přeměna je nespojitá změna vlastnosti soustavy nastávající při specifických p,T podmínkách Skupenství: plyn

tekutiny kondenzovaný stav

kapalina tuhá látka plazma Stav a rovnováha soustavy

Stav soustavy je determinován stavovými veličinami, nezávisí na způsobu – cestě jakou se soustava do daného stavu dostala. Např. p, T, V, n Vlastnosti soustavy vyjádřené stavovými veličinami: -

Extenzivní – mají aditivní charakter (jsou součtem částí) m, V, celková E

-

Intenzivní – nezávisí na látkovém množství – p, T, ρ, c, U

Teplo Q a práce W nejsou stavovými veličinami vztahují se k termodynamickému ději probíhajícímu z počátečního stavu do stavu rovnovážného. Termodynamická rovnováha: -

Mechanická (vyrovnání tlaků)

-

Tepelná (vyrovnání teplot)

-

Koncentrační

-

Fázová (mění-li se skupenský stav)

-

Chemická

Termodynamický děj Do rovnovážného stavu může soustava teoreticky přejít dějem vratným (reverzibilním), který je definován jako cesta nekonečně malými kroky nekonečně dlouhou dobu přes stavy nekonečně blízké rovnováze, přičemž směr děje lze kdykoli obrátit - např. expanze plynu při nekonečně pomalém snižování vnějšího tlaku. Nebo dějem nevratným (ireverzibilním), kdy přechod probíhá v konečném čase a v soustavě, resp. v jejím okolí dochází k nevratným změnám - např. expanze stlačeného plynu do vnějšího prostředí s atmosférickým tlakem. Reálné děje v přírodě probíhají nevratně. Typy dějů: - izotermický konstantní teplota - izobarický

konstantní tlak

[T] [p]

- izochorický konstantní objem

[V]

- adiabatický teplo se nevyměňuje s okolím

Q=0

- izoentalpický

konstantní entalpie

[H]

- izoentropický

konstantní entropie

[S] 4

1.2

Ideální plyn

Ideální plyn je jako zcela homogenní medium vymezen stavovými parametry: f(p,T,V,n) = 0,

f(p,T,Vm) = 0,

p = p(T, Vm)

Ve směsích pak: f(p,T,V,ni) = 0,

Vm je molární objem.

f(p,T,Vm,xi) = 0

Na rozdíl od reálného plynu je zcela zanedbána molekulová struktura a mezimolekulární interakce. Stavová rovnice ideálního plynu Pro stavové chování ideálního plynu platí: pV=nRT

resp.

p Vm = R T

(1.1)

R = 8,314 J. K-1.mol-1 je plynová konstanta Takže platí: p1 V1 = p0 V0

[T]

izoterma

(1.2)

V1 / T1 = V0 / T0

[p]

(1.3)

p1 / T1 = p0 / T0

[V]

(1.4)

Stavové chování ideálního plynu Stavové chování ideálního plynu lze charakterizovat parametry roztažnosti, stlačitelnosti a rozpínavosti. Koeficient izobarické roztažnosti:

 nRT  1  V  1 p p     V  T  p V  T  

   1  T  p

(1.5)

Koeficient izotermické stlačitelnosti:

T  

1  V  V  p

  T

 nRT  1 p  V  p  

    1 nRT  1  V p2 p  T

(1.6)

Koeficient izochorické rozpínavosti:

 nRT   p   V      V  T V  T   Samozřejmě platí:

    nR  p V T   V

V 

p T

(1.7)

(1.8)

5

Pro směsi ideálních plynů dále platí Daltonův zákon – celkový tlak je sumou parciálních tlaků složek. pn

RT RT  V V

pi  ni

k

k



 n    n i 1

i

i 1

i

RT  k    pi V  i 1

RT RT  xi n  xi p V V

(1.9) (1.10)

Amagatův zákon – celkový objem směsi je sumou parciálních objemů složek.

V n

RT RT  p p

k  RT  k  ni   Vi n    i p  i 1 i 1 i 1  k

(1.11)

Objemová práce ideálního plynu

W  F .x

(1.12) pF s F  p.s

(1.13)

dW   p.s.dx   p.dV

(1.14)

x2

W   Fdx x1

Objemová práce dodaná soustavě: V2

W    pdV

(1.15)

V1

Dodání objemové práce reprezentuje zmenšení objemu pracovního media soustavy (ideálního plynu). Ve smyslu obecně platného principu, kdy je vždy odečítán počáteční stav od konečného, je toto stlačení změnou zápornou. ΔW = - p ΔV

1.3

(1.16)

První věta termodynamická

První věta termodynamická je univerzálním principem, jehož podstatou je poznání, že jedinými formami výměny energie mezi soustavou a okolím může být teplo a práce. Teplo je mikrofyzikální - neuspořádaná forma výměny energie, práce je makrofyzikální, tj. uspořádaná forma výměny energie. Zvýšení vnitřní energie ΔU uzavřené soustavy odpovídá z okolí dodanému teplu a dodané práci. ΔU = U2 – U1 = W + Q

(1.17)

Totéž platí samozřejmě v diferenciálním tvaru pro nekonečně malou (infinitezimální) změnu vnitřní energie: 6

dU = dQ + dW

(1.18)

Resp. při dodané objemové práci: dU = dQ – pdV

(1.19)

Práce a teplo jsou dle první věty termodynamické ekvivalentní, změna stavu soustavy vyjádřená ΔU nezávisí na tom, v jakém poměru se na ní obě formy výměny energie podílí. Pro cyklický děj platí definičně:

 dU  0

(1.20)

Entalpie Využívání rozdílu stavů soustavy daného vnitřní energií ΔU pro bilanční vyhodnocování průběhu chemických reakcí za standardních izobarických podmínek (s možností libovolné změny objemu) je však z hlediska měřitelnosti objemové práce problematické. Zřejmým východiskem je vyhodnocování chemických reakcí na základě reakčního tepla entalpie H, která je dobře měřitelná (kalorimetricky). Výsledné reakční teplo (uvolněné či spotřebované) za konstantního tlaku odpovídá nejen změně vnitřní energie, ale i vykonané či dodané (spotřebované) objemové práci. dQ = dU + d(pV) = d(U + pV) = dH

(1.21)

za konst. tlaku samozřejmě platí: pdV = d(pV)

Mezi dQ a dH je tedy rovnítko a entalpie je označována jako reakční teplo, nicméně stále je nutno rozlišovat, že entalpie je na rozdíl od tepla definována jako stavová veličina. Pro entalpii tedy platí: dH = dU + pdV

(1.22)

ΔH = ΔU + p ΔV

(1.23)

resp.: Přírůstek entalpie je roven teplu přijatému za konstantního tlaku soustavou, nekonala-li se jiná práce než objemová. Entalpie je pro vyhodnocování procesů (zjm. chemických reakcí) výhodná – je měřitelná. V obecnosti může být s okolím vyměněna i jiná forma práce než práce objemová. Tepelné kapacity Tepelná kapacita je definována na základě diferenciálního množství tepla dQ potřebného k diferenciálnímu zvýšení teploty dT.

C

dQ dT

(1.24)

Tepelná kapacita je definována za konst. objemu - CV a konst. tlaku - Cp.

CV 

dQV  U    dT  T V

(1.25)

7

Cp 

 H    dT  T  p

dQ p

(1.26)

Za konstantního objemu se teplo nespotřebovává na objemovou práci - na zvýšení teploty tedy stačí méně tepla, tj. obvykle platí C p  CV (nerovnost neplatí např. u vody v teplotním rozmezí 1 – 4 oC, kdy s růstem teploty dochází k teplotní kontrakci). Rozdíl mezi tepelnou kapacitou za konstantního tlaku Cp a konst. objemu CV je vyčíslen:

 H   U   U   V   U  C p  CV         p     T  p  T V  T  p  T  p  T V

(1.27)

Pro totální diferenciály stavových funkcí U, V v prostoru nezávisle proměnných platí:

 U   U  dU    dV    dT  V T  T V

(1.28)

 V  V  dV    dT    T  p  p

(1.29)

  dp T

Po dosazení za dV z (1.29) do (1.28) a derivaci podle T je obdrženo:

 U   U   V   U           T  p  V T  T  p  T V

(1.30)

Výsledný vztah pro rozdíl tepelných kapacit má tedy tvar:

  U   V  C p  CV   p       V T  T  p 

(1.31)

Aplikace první věty termodynamické na ideální plyn Pro ideální plyn vedle stavové rovnice (1.1) platí:

 U    0  V T

(1.32)

Neboť dodané (resp. odebrané) teplo při izotermické expanzi (resp. kompresi) je u ideálního plynu vždy ekvivalentní objemové práci, kterou tento plyn vykonal (resp. která mu byla dodána). Totální diferenciál vnitřní energie ideálního plynu je tedy redukován na:

 U  dU    dT  CV dT  T V

(1.33)

Resp. vnitřní energie ideálního plynu závisí pouze na teplotě. Tepelnou kapacitu ideálního plynu při konst. objemu lze tedy vyjádřit: CV 

dU dT

(1.34)

8

A vztah pro rozdíl tepelných kapacit je u ideálního plynu s ohledem na diferencovaný tvar stavové rovnice (1.1) redukován na konstantu:

 V  C p  CV  p   nR  T  p

(1.35)

Pro Změnu U a H při teplotním přechodu z T1 na T2 samozřejmě vychází: T2

U   CV dT

(1.36)

T1

T2

H   C p dT

(1.37)

T1

Při libovolné izotermické stavové změně ideálního plynu zůstává jeho vnitřní energie konstantní – viz (1.32).

 U   U  dU    dT    dV  0  T V  V  p

(1.38)

Pak tedy:

dQ  pdV  nRT

dV V

(1.39)

Vratný přechod ideálního plynu ze stavu 1 do stavu 2, tj. izotermickou změnu objemu nebo tlaku lze vyjádřit: 2

2

1

1

 dQ   nRT Q  nRT ln

dV V

V2 p  RT ln 1 V1 p2

(1.40)

Naopak pro případ vratné adiabatické expanze ideálního plynu dle první věty termodynamické (dQ = 0, dU = -pdV) platí:

 pdV  CV dT

(1.41)

Při uvážení stavové rovnice (1.1):

CV

dT dV  nR 0 T V

(1.42)

Po integraci pro přechod ze stavu 1 do stavu 2: CV ln

T2 V  nR ln 2  0 T1 V1

(1.43)

Po uvážení (1.35):

9

ln

T2 V    1 ln 2  0 , T1 V1

T1  V2    T2  V1 



Cp CV

(1.44)

 1

(1.45)

Pro stavové chování ideálního plynu platí:

T1 pV  1 1 T2 p2V2

(1.46)

Takže pro výslednou adiabatu (v p, V prostoru) ideálního plynu je z (1.45) a (1.46) obdrženo: 



p1V1  p2V2 ,

1.4

pV   konst.

(1.47)

Druhá věta termodynamická

Další aspekty energetické bilance procesů probíhajících v prostředí (včetně chemických reakcí) již souvisejí s druhou větou termodynamickou, která vyjadřuje principiální vlastnost samovolných procesů probíhat ve směru k větší neuspořádanosti. Např. odrazy míčku jsou při dopadech na pevný podklad samovolně tlumeny v důsledku přeměny práce na teplo (míček i podklad se zahřívá). Opačný proces - spontánní zahájení vertikálních odrazů zahřátého míčku (ze zahřátého podkladu) až do jeho navrácení do výchozího bodu je nemožné, i kdyby byl celý proces důsledně adiabaticky izolován. Tato nemožnost, která je zřejmá sama o sobě, obecně souvisí s uspořádaností dějů. Vzniklé teplo je sice dle první věty termodynamické ekvivalentní počáteční potenciální energii míčku, nicméně možnost, že by v rámci chaotického tepelného pohybu částic (molekul) systému nastala zcela uspořádaná konfigurace vedoucí k samovolnému odskoku míčku od podkladu, je krajně nepravděpodobná.

Druhá věta termodynamická v jedné z mnoha možných formulací konstatuje: Nelze přeměnit teplo samovolným dějem na ekvivalentní práci. Carnotův cyklus a jeho účinnost Důsledky tohoto principu lze rozpracovat v rámci popisu účinnosti tzv. Carnotova tepelného stroje. Jde o čistě abstraktní konstrukci cyklicky pracujícího tepelného stroje reprezentovanou dvěma tepelnými zásobníky s teplotami T2 > T1 a vlastní soustavou vyměňující si s těmito zásobníky teplo a odevzdávající či přijímající práci od okolí. Tepelné toky mezi soustavou a oběma zásobníky nikterak neovlivňují teploty zásobníků, T2 a T1 jsou konstantní. V rámci jednoho cyklu proběhnou 4 reverzibilní děje a míra neuspořádanosti soustavy je na konci cyklu rovna počáteční. Pracovním mediem je ideální plyn (1 mol). Carnotův cyklus je popsán následovně: 1. Izotermická vratná expanze při teplotě T2 – plyn expanduje z objemu V1 na V2, přičemž z tepelného zásobníku o teplotě T2 odebere teplo +Q2 a vykoná (do okolí odevzdá) práci –W1.

10

Obr. 1.1. Carnotův cyklus

2.

Adiabatická vratná expanze – plyn expanduje bez výměny tepla s okolím na objem V3, přičemž vykoná (do okolí odevzdá) práci -W2 a ochladí se na teplotu T1.

3.

Izotermická vratná komprese při teplotě T1 – plyn je vynaložením práce +W3 (tj. přijetím práce od okolí) stačen na objem V4, přičemž chladnějšímu zásobníku odevzdá teplo -Q1.

4.

Adiabatická vratná komprese – plyn se bez výměny tepla s okolím stlačí vynaložením práce +W4 zpět na objem V1, přičemž se ohřeje na výchozí teplotu T2

Stroj pracuje cyklicky, vnitřní energie pracovního plynu se po proběhnutí cyklu vrátí na výchozí hodnotu ΔU = 0 - viz (1.20). Účinnost tohoto stroje, jenž přeměňuje teplo na práci, je samozřejmě dána jako poměr vykonané práce k teplu odebranému z tepelného zdroje (teplu přijatému pracovním mediem).



W Q2

(1.48)

Lze doplnit, že vztah (1.48) je obecnou relací vyjadřující účinnost (necyklického) tepelného stroje v závislosti na velikosti podílu odebraného tepla, který je přeměňován na práci -W odevzdávanou do okolí spolu se zbytkovým – nevyužitým teplem. Při vyhodnocování účinnosti Carnotova cyklicky pracujícího stroje je podstatné, že celková vykonaná práce je (dle první věty termodynamické) ekvivalentní tepelné bilanci pracovního cyklu. Platí:

U  Q  W  Q1  Q2  W  0

(1.49)

 W  Q1  Q2

(1.50)

Vyjádřená celková práce má v (1.50) záporné znaménko, jde tedy o práci soustavou vykonanou, která je uvažována v relaci pro účinnost (1.48). Tato vykonaná práce je s ohledem na znaménko Q1 rozdílem mezi teplem soustavě dodaným v prvním ději, tj. Q2 a teplem odevzdaným v třetím ději Q1 (množství tepla vyměňovaného při izotermické expanzi či kompresi ideálního plynu je samozřejmě závislé na nastavené teplotě  Q1  Q2 ). Tj. práce vykonaná v pracovním cyklu Carnotova tepelného stroje je jednoduše vyjádřená prostřednictvím tepelné bilance bez toho, že by bylo nutno vyhodnocovat práce -W1, -W2, W3, W4 realizované v jednotlivých krocích.

11

Teplo Q2 přijaté plynem v prvním kroku odpovídá vykonané objemové práci:

Q2  W1 

V2

V2

 pdV   RT

2

V1

V1

dV V  RT2 ln 2 V V1

(1.51)

Obdobně – plynem dodané teplo Q1 odpovídá práci: Q1  W3  RT1 ln

V4 V3

(1.52)

Pro vratné adiabatické děje platí pro ideální plyn obecný vztah mezi objemem a teplotou TVγ-1=konst. (γ=Cm,p/Cm,V je Poissonova konst.). Takže pro kroky 2 a 4 můžeme psát: T2V2 1  T1V3 1

T1V4 1  T2V1 1

(1.53)

V2 V3  V1 V4

(1.54)

Tj.:

Sloučením rovnic (1.52) a (1.54) pak pro Q1 vychází: Q1   RT1 ln

V2 V1

(1.55)

Následně tedy lze psát: Q1  Q2  RT2  T1  ln

V2 V1

(1.56)

Účinnost pak můžeme s uvážením vztahů (1.49), (1.51) a (1.56) vyjádřit:



 W Q1  Q2 T2  T1   Q2 Q2 T2

(1.57)

Vychází tedy, že účinnost Carnotova tepelného stroje je determinována pouze teplotami zásobníků. Např. situaci T2 = 2 T1 odpovídá 50% účinnost (100% by odpovídalo buď T2   nebo T1  0 , což je v obou případech nemožné). Dále je nutno dodat, že vztah (1.57) vyjadřuje maximální teoreticky možnou účinnost tepelného stroje pracujícího v teplotním intervalu T1 až T2. Resp. maximální účinnost může být dosažena pouze pro teoretickou situaci, kdy všechny děje proběhnou vratně. Ze vztahu (1.57) lze vyvodit relaci: Q1 Q2  0 T1 T2

(1.58)

Vztah (1.58) má z hlediska studovaných souvislostí zásadní význam. Jak již je zdůrazňováno výše, všechny čtyři kroky popsaného Carnotova cyklu probíhají reverzibilně. Resp. lze říci, že pokud je suma zlomků na levé straně rovnice (1.58) rovna nule, jsou všechny procesy probíhající v rámci daného cyklu vratné a v soustavě po proběhnutí cyklu nenarůstá neuspořádanost. Carnotův cyklus lze vyjádřit i jako diferenciál s infinitezimálními úseky obou izoterm:

12

dQ1 dQ2  0 T1 T2

(1.59)

Zobecnění vztahu (1.59) pro libovolný vratný cyklický proces s libovolným počtem jednotlivých infinitezimálních Carnotových cyklů vede k integrálnímu tvaru:



dQ 0 T

(1.60)

Probíhá-li některý s dílčích dějů v Carnotově cyklu nevratně, je nevratný celý proces a jeho celková účinnost musí být nižší oproti primárně uvažované ideální situaci, kdy jsou všechny čtyři dílčí děje vratné.



T2  T1 T2

(1.61)

Tato nerovnost může být rozepsána obdobně vztahu (1.57), kdy pravá strana reprezentuje teploty T1 , T2 tepelných zásobníků. Tj teploty, které by měla soustava (1 mol ideálního plynu) v každém okamžiku izotermické expanze či komprese při jejich vratném průběhu. Při nevratném průběhu (jednoho nebo více dílčích dějů) není soustava a tepelný zásobník v každém okamžiku v teplotní rovnováze a účinnost odpovídající soustavě samé vyjádřená vyměňovanými teply je menší než účinnost odpovídající vratnému procesu vyjádřená teplotami vnějších zásobníků. Q1  Q2 T2  T1  Q2 T2

(1.62)

Po úpravě: Q1 Q2  0 T1 T2

(1.63)

Tj. faktor korespondující s teplem dodaným soustavě Q2 při teplotě T2 vychází menší než zlomek odpovídající teplu soustavou odevzdanému –Q1/T1. Lze říci, že v případě nevratného průběhu jednoho nebo více dílčích dějů Carnotova cyklu musí soustava odevzdat více neuspořádané formy energie, než by odevzdala, pokud by všechny děje proběhly vratně. Tj. ne všechno přijaté teplo se zužitkuje na práci, zůstává zbytkové. V tomto smyslu tedy v soustavě narůstá neuspořádanost, resp. Carnotův tepelný stroj se logicky musí zastavit. Tuto skutečnost lze jednoduše vykreslit např. na adiabatické expanzi, kdy by plyn při nevratném průběhu vykonal práci menší než by odpovídala ději vratnému a jeho teplota by byla v souvislosti se zbytkovým teplem vyšší než T1. V následujícím kroku – izotermické kompresi by plyn musel do tepelného zásobníku odevzdat více tepla Q1, než kdyby adiabatická expanze proběhla vratně. Stejný princip úbytku uspořádané formy energie by se jistě realizoval i u zbývajících tří (nevratných) dějů (dorovnávání teploty pracovního plynu na teploty T1 a T2 by samozřejmě probíhalo pouze v rámci izotermických dějů). Stejně jako pro vratně cyklické procesy lze zobecněním dojít k integrální podobě principu vyjádřeného v (1.63). Pro nevratný obecný cyklus (kdy v soustavě roste neuspořádanost) tedy platí Clausiova nerovnost:



dQ 0 T

(1.64)

13

Lze shrnout, že abstraktní konstrukce Carnotova tepelného stroje má pro interpretaci principu II. věty termodynamické zásadní význam v tom, že jeho účinnost lze vyjádřit vztahem zahrnujícím pouze pracovní teploty uvažovaných izotermických dějů. Takto vyvozená účinnost je maximální – odpovídající vratným dějům, neboť obě teploty jsou udržovány vně soustavy vysoce kapacitními tepelnými zásobníky a jsou konstantní. Druhou možností pro vyjádření účinnosti Carnotova stroje je relace zahrnující tepla přijímaná a odevzdávaná při izotermických dějích, která se vztahuje k pracovnímu plynu – vyhodnocuje účinnost soustavy samotné a může tedy zahrnovat dočasnou nerovnováhu s teplotami T1 a T2. Při porovnání účinností definovaných těmito dvěma způsoby je postulováno, že účinnost vyvozená z tepel je v případě nevratnosti některého z dějů Carnotova stroje nižší, což vystihuje podstatu druhé věty termodynamické. Pro exaktní (matematickou) formulaci II. VT jsou pak k dispozici právě dvě základní termodynamické veličiny vyhodnocující účinnost Carnotova tepelného stroje – teplo a teplota. Pojem entropie V kontextu skutečnostmi vyjádřenými vztahy (1.59) a (1.60) a úvahami v předchozí stati (a také proto, že poměr dQ/T je totálním diferenciálem) lze zavést novou stavovou funkci nazvanou entropie - S charakterizující míru neuspořádanosti termodynamické soustavy. Pro vratný děj platí: dS 

dQ T

(1.65)

Definice entropie tedy vyjadřuje infinitezimální výměnu tepla mezi soustavou a okolím v poměru k teplotě soustavy, která je konstantou. Tento princip v zásadě koresponduje se situací skákajícího míčku v příkladu z úvodu kapitoly. Deformace míčku a podkladu reprezentují vznik tepla, teplota soustavy se však znatelně nezmění. Pádem míčku naroste entropie. Vratný přechod soustavy ze stavu A do B je dán: B

dQ T A

S  

(1.66)

Tj. pro výše analyzovaný Carnotovův cyklus je změna entropie při vratné izotermické expanzi nebo kompresi ekvivalentní zlomku Q2/T2, resp. -Q1/T1, vratné adiabatické objemové změny jsou ději izentropickými (dQ =0, resp. S  0 ). Změnu entropie soustavy tedy lze jednoduše vyčíslit (na základě vztahu (1.66)) pouze pro vratný děj, pro nevratný proces platí nerovnost: dS 

dQ T

(1.67)

Pro nevratný přechod ze stavu A do B: B

dQ T A

S  

(1.68)

Jde-li o nevratný adiabatický proces, platí jednoduše: dS  0 , S  0

(1.69)

Z těchto nerovností samozřejmě vyplývá zákonitost, že entropie soustavy při průběhu nevratného děje vždy vzrůstá (= II.VT). Nerovnost (1.67), resp. (1.68) logicky plyne 14

z předchozího. Nárůst entropie nevratného adiabatického děje lze vyvodit i z úvah o přechodu tepla mezi dvěma tělesy o různých teplotách, která dohromady tvoří tepelně izolovanou soustavu. Pro změnu entropie této soustavy související s odevzdáním tepla -dQ1 tělesem o teplotě T1 tělesu o teplotě T2, které přijme teplo dQ2 (dQ2= -dQ1) lze psát: dS  dS1  dS 2 

dQ1 dQ2  1 1  T  T1     dQ1  2 dQ1 T1 T2 T1T2  T1 T2 

Z čehož resultuje dS  0 , neboť teplo přechází samovolně (nevratně) pouze z tělesa teplejšího na chladnější. Tj. dQ1 < 0 při T2-T1 < 0, nebo naopak dQ1 > 0 při T2-T1 > 0, jiné možnosti nejsou. Stejný výsledek vyplývá i z Clausiovy nerovnosti (1.63). Clausiova nerovnost platí pro cyklus, v kterém alespoň jeden s dílčích dějů proběhne nevratně (nevratný je pak celý cyklus). Nejjednodušším modelovým příkladem je cyklus, kdy soustava přejde ze stavu A do stavu B nevratně (a bez výměny tepla s okolím), zpět do výchozího stavu A se pak musí samozřejmě vrátit "po jiné cestě" - vratně a s možností výměny tepla s okolím. Rozpis nerovnosti (1.64) pro tento případ je následující: B

A

dQ dQ A T ir  B Trev  0

(1.70)

První integrál je roven nule, neboť při přechodu A → B je soustava izolovaná a dQir = 0, druhý integrál je dle (1.61) roven rozdílu entropií SA - SB . Má-li být tento rozdíl menší než nula, musí být entropie soustavy ve stavu B vyšší než ve stavu A - entropie v průběhu nevratného děje A → B vzrůstá. Tento závěr by samozřejmě resultoval i pokud by bylo pořadí dějů zvoleno obráceně.

1.5

Gibbsova energie

Pro průběh samovolného procesu soustavy (tepelně neizolované) platí: dS 

dQ T

(1.71)

Probíhá-li tento děj za izotermicko-izobarických podmínek lze dQ nahradit změnou entalpie dH. Po úpravě je obdržena nerovnost podmiňující samovolný průběh děje:

dH  TdS  0

(1.72)

Tato nerovnost vede k zavedení termodynamické stavové funkce konzistentně zahrnující obecné termodynamické principy umožňující vyhodnocovat samovolnost průběhu dějů za izotermicko-izobarických podmínek - Gibbsovy energie. ΔG = ΔH - T ΔS dG = dH - T dS

(1.73)

Hodnota ΔG je tedy kritériem (resp. potenciálem) samovolného průběhu termodynamických dějů. Chemická reakce (nebo jakýkoli jiný proces) proběhne, vychází-li celková bilance Gibbsovy energie (konečný minus počáteční stav procesu) záporná. U vícesložkových procesů, jako jsou chemické reakce ap., je celková Gibbsova energie jakožto potenciál průběhu (do rovnováhy) ve směru reaktanty → produkty dána rozdílem sumárních

15

slučovacích ΔGo reaktantů a produktů (ΔGoprodukty - ΔGoreaktanty). Hodnoty ΔGo pro jednotlivé reakční složky (sloučeniny) jsou samozřejmě dostupné v tabulkách. Záporná hodnota Gibbsovy energie odpovídá (vedle změny entropie) teplu uvolněnému v průběhu reakce za izobaricko-izotermických podmínek. Při reálných reakcích zpravidla není teplo odvedeno beze zbytku - ohřejí se i reagující látky, nicméně fakt, že soustava není adiabaticky izolovaná, opravňuje hovořit o izotermických podmínkách. Gibbsova energie je definována tak, že reflektuje finální tepelnou a entropickou bilanci procesu. Tj. včetně tepla spotřebovaného či uvolněného v souvislosti s případnou objemovou prací. Lze dodat, že bilanční relace (1.72) v principu neznamená, že samovolnost libovolného procesu je automaticky podmíněna zvýšením neuspořádanosti (ΔS > 0). Tento vztah říká pouze to, že případný nárůst uspořádanosti soustavy (snížení entropie - minus ΔS resp. záporný faktor TΔS) je menší než hodnota, jež by odpovídala záporné entalpii procesu. Viz. druhá věta termodynamická, která nevylučuje přeměnu tepla na práci, říká pouze, že tato přeměna nemůže být ekvivalentní. Tj. snížení neuspořádaností - snížení entropie není překážkou dílčího samovolného procesu. Výše zmíněný obecný princip zvyšování neuspořádanosti (princip nevratnosti) dějů probíhajících v přírodě samovolně ve své obecnosti neřeší důsledky možného rozdílu vnitřní energie počátečního a konečného stavu soustavy. S tímto principem však bezesporu koresponduje samovolný cyklický proces definovaný relací (1.20), který může zahrnovat i dílčí děje popisovaného typu s poklesem entropie. Pro cyklické procesy tedy princip samovolnosti říká, že soustava se po absolvování cyklu zahrnujícího libovolný počet dílčích procesů, při kterých může docházet k změně vnitřní energie, výměně tepla a práce s okolím a případně k zvýšení či snížení neuspořádanosti, vrátí do výchozího stavu více neuspořádaná. Pro samovolně proběhlý cyklus samozřejmě platí (1.20), přičemž celkové ΔG < 0 a ΔS > 0. Zde se dozajista nabízí bezpočet příkladů figurujících např. v rámci biogeochemických cyklů jednotlivých prvků. Např. uhlík se může nacházet ve stavu, kdy je obsažen v biomase, jež je následně za uvolnění tepla spálena - C přejde ve formě CO2 (jenž má mnohem nižší vnitřní energii než organická forma uhlíku v biomase) do atmosféry. Následně se tento CO2 může při kontaktu s rostoucí zelenou biomasou opět vrátit při pohlcení odpovídajícího množství energie slunečního záření do organické formy - cyklus je uzavřen.

Popsaný proces je stejně jako všechny ostatní přirozené děje, které probíhají v prostředí využívajíc energie slunečního záření, samovolný. Tj. v průběhu cyklu dojde vždy (dle druhé věty termodynamické) k nevratnému zvýšení neuspořádanosti - nárůstu entropie. Je-li princip průběžného růstu entropie vyplývající z druhé věty termodynamické pro procesy na Zemi a ve vesmíru bezezbytku platný, jsou jeho důsledky samozřejmě zásadní. Otázky spojené s tímto tématem však již vyvstávají spíše v rovině filozofické než striktně fyzikálněchemické. Helmholzova energie Helmholzova volná energie ΔA se vztahuje k izochorickému ději, kdy platí ΔU = ΔQ, tj. nekoná se objemová práce. Obdobně jako u Gibbsovy energie je Helmholzova energie potenciálem průběhu izotermicko – izochorického děje: ΔA = ΔU - T ΔS dA = dU - T dS

(1.72)

16

Jako kritérium (potenciál) průběhu chemických reakcí probíhajících nejčastěji za atmosférického tlaku tj. izobaricky je samozřejmě vhodnější Gibbsova energie.

1.6

Termodynamické potenciály

Pro diferenciály veličin reprezentujících termodynamické potenciály – U, H, G, A lze vyvodit následující relace. Vnitřní energie – spojením vztahů vyplývajících z I. a II. druhé věty termodynamické (1.19) a (1.64), resp. T dS = dQ (vratný děj) je obdrženo: dU = -pdV + TdS

přirozené proměnné (V, S)

(1.73)

Entalpie – dle první věty termodynamické je dán úplný diferenciál dH = dU + pdV+Vdp, kombinací se vztahem (1.73) vzniká: dH = Vdp + TdS

přirozené proměnné (p, S)

(1.74)

Gibbsova energie – dle (1.71) je úplný diferenciál ve tvaru dG = dH – TdS - SdT, spojením se vztahem (1.74) vzniká: dG = Vdp - SdT

přirozené proměnné (p, T)

(1.75)

Helmholtzova energie – dle (1.72) musí být úplný diferenciál dán ve tvaru dA = dU – TdS SdT, spojením se vztahem (1.73) vzniká: dA = -pdV - SdT

přirozené proměnné (V, T)

(1.76)

Maxwelovy vztahy Vztahy pro termodynamické potenciály (1.73 – 76) lze jednoduše psát pomocí parciálních diferenciálních kvocientů:

 U   U  dU    dV    dS  V  S  S V

(1.77)

 H   H   dp   dH    dS  S  p  p  S

(1.78)

 G   G   dp   dG    dT  T  p  p T

(1.79)

 A   A  dA    dV    dT  V T  T V

(1.80)

Takže vychází:

 U   U     -p,   T  V  S  S V

(1.81)

17

 H   H     V ,   T  S  p  p  S

(1.82)

 G     V ,  p T

(1.83)

 G     S  T  p

 A   A      p,    S  V T  T V

(1.84)

Pro parciální diferenciální kvocienty v (1.81) lze dále vyvodit:

   U    U   p     -p,         V  S  S V  S  V  S V

(1.85)

 U   T     U      T,         S V  V  S  V  S V  S

(1.86)

Pořadí derivací je zaměnitelné – platí Eulerův reciproční vztah:

 T   p        V  S  S V

(1.87)

Stejným způsobem jsou obdrženy další výsledné Maxwelovy vztahy:

1.6

 T   V        p  S  S  p

(1.88)

 S   V        T  p  p T

(1.89)

 p   S       T V  V  T

(1.90)

Třetí věta termodynamická

Obsahem třetí věty termodynamické je tvrzení: Žádným postupem, ať jakkoliv idealizovaným, nelze u žádné soustavy dosáhnout snížení její teploty na absolutní nulu konečným počtem operací. Entropie látky je při teplotě absolutní nuly limitní – nabývá konečné nejnižší možné hodnoty S0 , resp. dle Plancka je nulová (což pravděpodobně platí pouze pro čisté látky ve stavu dokonalého krystalu). lim S  S 0 , resp lim S  0 (1.91) T 0

T 0

V limitě T → 0 však zcela jistě platí ΔS0 = 0 pro libovolný vratný děj. Snížení teploty při ochlazovací operaci (kroku) v limitě T → 0 bude nulové bez ohledu na vynaloženou energii. 18

Princip třetí věty termodynamické se vedle klasických experimentů založených na izotermickém stlačování a adiabatické expanzi pracovního plynu prakticky projevuje při snižování teploty adiabatickou demagnetizací soli prvků vzácných zemin s vysokými hodnotami magnetické susceptibility. Tímto experimentálním postupem byly dosud dosaženy teploty, které se mohou jevit velmi blízké absolutní nule (zlomky K, nicméně dosažení 0K je nemožné). Křivka B = Bi na obr. 1.2 odpovídá zmagnetizovanému stavu s orientovanou (tj. uspořádanou) strukturou, B = 0 je stav bez působení magnetického pole. Adiabatická demagnetizace - krok 2→3 apod. je izoentropický děj, tj. zvýšení strukturní neuspořádanosti po vypnutí magnetického pole musí být kompenzováno poklesem teploty, resp. strukturně neuspořádanější stav odpovídá stejné entropii pouze po snížení teploty. Z grafu je zřejmé, že v blízkosti absolutní nuly (0 K) je snížení teploty dosahované demagnetizačním krokem blízké nule.

Obr. 1.2 S – T graf adiabatické demagnetizace

2. Rovnovážná konstanta Gibbsovou energií

chemické

reakce

a

její

relace

s

Smysl pojmu chemická rovnováha souvisí se skutečností, že chemické reakce obecně neprobíhají do úplného vymizení výchozích látek, nýbrž dospějí do rovnovážného stavu, kdy jsou vedle produktů přítomny i reaktanty, byť často v nepatrném množství. Chemickou reakci A+B↔X+Y lze chápat jako dva dílčí procesy - jeden probíhající ve směru zleva doprava rychlostí v1 a druhý jdoucí v opačném směru rychlostí v2. Jde-li o reakci exotermní s hladkým průběhem, bude samozřejmě pro počátek procesu platit v1 >> v2. Rovnováha nastane v okamžiku, kdy se obě rychlosti vyrovnají. Tyto dílčí rychlosti lze vyjádřit na základě Guldbergova a Waageho zákona o působení "aktivní hmoty". Tento zákon reflektuje skutečnost, že okamžitá rychlost chemické přeměny libovolné látky musí být přímo úměrná jejímu relativnímu množství, resp. koncentraci. Pro dvě reagující látky pak bude tato rychlost logicky úměrná součinu jejich aktuálních koncentrací (pravděpodobnost setkání molekul příp. atomů reagujících látek determinující reakční rychlost je samozřejmě úměrná součinu jejich koncentrací). Faktorem kvantifikujícím tuto úměrnost je rychlostní konstanta, která je pro každý proces specifická. Tj. platí:

v1  k1c AcB v2  k 2 c X cY

(2.1)

Pro stav rovnováhy pak: 19

k1c AcB  k2c X cY

(2.2)

Poměr rychlostních konstant k1/k2 je brán jako rovnovážná konstanta dané reakce K.

k1 c c K X Y k2 c AcB

(2.3)

Vztah (2.3) lze zobecnit pro reakci s libovolnou stechiometrií aA + bB ↔ xX + yY a přepsat do korektnějšího tvaru za použití aktivit místo koncentrací.

K

a Xx .aYy a Aa .aBb

(2.4)

Aktivita vystihuje reálný potenciál chování látky v roztoku. U koncentrovanějších roztoků se aktivita od molární koncentrace (odpovídající molárnímu množství látky nasypanému do roztoku) výrazněji odchyluje směrem k nižším hodnotám, neboť každý iont je ve svých chemických projevech "tlumen" působením opačně nabitých iontů okolních, tzn. vliv iontové atmosféry – viz následující kapitola 3. Hodnoty rovnovážných konstant jsou pro standardní chemické reakce prováděné za izobarických podmínek k dispozici v tabulkách. Rovnovážná konstanta chemické reakce samozřejmě umožňuje determinovat rovnovážné koncentrace (resp. aktivity) látek účastnících se reakce pro různé podmínky určené koncentracemi (aktivitami) ostatních složek. Lze jen dodat, že vše platí pouze pro ustálený rovnovážný stav. Vyšší hodnota K reprezentuje hladký průběh reakce zleva doprava, kdy je po dosažení rovnováhy v reakčním systému výrazný přebytek produktů. Tomuto případu samozřejmě přísluší záporný rozdíl Gibbsovy energie rovnovážného (konečného) a výchozího stavu. Tato souvislost naznačuje zřejmou provázanost mezi K a ΔG. Odvození vztahu mezi Gibbsovou energií a reakční konstantou vychází ze vztahu (1.71), z kterého pro izotermické podmínky vyplývá, že změna Gibbsovy energie při změně tlaku je ekvivalentní objemu –viz (1.83)  dG     V  dp T

Při uvážení (1.1) vychází pro izotermické podmínky:

dG  nRT

dp p

(2.5)

Pro přechod plynu ze standardního stavu (tj. stavu za standardního tlaku 101,325 kPa) charakterizovaného veličinami Go a po do libovolného jiného stavu vyjádřeného jako G, p platí:  p  dp  nRT ln  o   nRT ln pr p p  po p

G  G  nRT o



(2.6)

Tj. podíl p/po je vyjádřen jako relativní tlak pr. Rovnici (2.6) lze zobecnit pro směs k složek - ideálních plynů, jenž spolu budou chemicky reagovat.

20

k

Gr  Gr0  RT  i ln pr ,i

(2.7)

i 1

V tomto případě se již nejedná o přechod soustavy do libovolného stavu. Konečný stav charakterizovaný členem ΔGr odpovídá reakční směsi ideálních plynů v rovnováze, tj. po zreagování, kdy chemický proces dospěl do energetického minima. Počáteční stav je vyjádřen bilancí výchozích energetických stavů, tj. standardních slučovacích Gibbsových energií ΔGosl jednotlivých složek reakční směsi. Reaktanty se v průběhu děje spotřebovávají - jejich ΔGosl je tedy započítávána se záporným znaménkem. l

m

i 1

i 1

0 0 Gr0   i Gsluč . prod.,i   i Gsluč .reakt.,i

(2.8)

Člen ΔGor označuje standardní reakční Gibbsovu energii (potenciál průběhu reakce), která je vztažena na na 1 mol základních reakčních obratů, tj. v rovnici (2.7) již nefiguruje symbol n označující počet molů. Pro reakční směs naopak přibyly stechiometrické koeficienty νi jednotlivých složek. Rovnice (2.7) tedy popisuje proces, na jehož počátku je směs k ideálních plynů, jejichž výchozí parciální tlaky poi odpovídají jedničce (relativní parciální tlaky pr,i samozřejmě taktéž, na počátku pi = poi). V případě ΔGor < 0 budou parciální tlaky složek, jež jsou brány jako produkty, narůstat a pr,i reaktantů budou klesat. Děj bude za pozvolného snižování rychlosti chemických přeměn spět do rovnovážného stavu, okamžitá Gibbsova energie procesu ΔG se bude vzdalovat od počáteční hodnoty ΔGor a bude spět k nule (zdola). Pro stav rovnováhy platí ΔGr = 0 (za konst. T a p). S přihlédnutím k bilanční rovnici (2.8) lze psát: m  l  Gr0   RT   i ln pr . prod,i   i ln pr .reakt,i  (2.9) i 1  i 1  Přičemž relativní parciální tlaky reakčních složek pr v rovnici (2.9) odpovídají konečnému rovnovážnému stavu. Chemickou přeměnu každé složky (ideálního plynu) reakční směsi lze tedy brát jako samostatný proces se změnou parciálního tlaku odpovídající ΔGosluč. této složky. Vztah (2.9) lze samozřejmě převést do tvaru:

G   RT ln 0 r

pr .1prod,1  pr .l prod,l pr .1reakt,1  pr .mreakt,m

  RT ln K p

(2.10)

Výsledkem je tedy relace standardní reakční Gibbsovy energie a rovnovážné reakční konstanty Kp (pro chemické reakce v plynném skupenství). Pro reakce v roztoku lze vedle podmínek izotermických uvažovat i neměnnost tlaku – podmínky izobarické. Oba členy vyplývající z formulace dG jako obecného termodynamického potenciálu (1.79) jsou tedy pro izotermické děje v roztoku nulové. V souvislosti s chemickými změnami je změna Gibbsovy energie ekvivalentní změně molového množství jednotlivých látek n v součinu s jejich chemickým potenciálem μ.

dG  SdT  Vdp    i dni    i dni i 1

[T, p]

(2.11)

i 1

Pro změnu molového množství reakční složky platí:

dni   i di

(2.12)

21

 je rozsah reakce definovaný jako poměr látkového množství libovolné vznikající či spotřebovávané reakční složky a k jejímu stechiometrickému koeficientu.

i 

ni

(2.13)

i

Tj.:

  dG    i  i d  i 1   G      i  i    T , p i 1

[T, p] (2.14)

Pro rovnováhu platí:

 G     0    T , p    0  Gr   i  i   i 1  rovn

(2.15)

Pro chemický potenciál látky ai (reakční složky) jako pro ekvivalent Gibbsovy energie - viz (2.11) platí: (2.16)  i   i0  RT ln ai Vztah (2.16) je obdobou (2.6) (aktivita v roztoku koresponduje s relativním tlakem v směsi plynů). Pro aktivity reakčních složek odpovídající rovnovážnému stavu lze psát: k





k

k

k

Gr   i   RT ln ai  i   RT ln  ai  G  RT ln  ai i i 1

0 i

i 1

0 i

i 1

i

0 r

(2.17)

i 1

S uvážením znaménkové konvence vyjádřené v (2.8) je obdržena výsledná relace mezi standardní reakční Gibbsovou energií ΔGor a rovnovážnou konstantou Ka pro chemické reakce v roztoku:

G   RT ln 0 r

l 1 aprod ,1  a prod,l m 1 areakt ,1  areakt, m

  RT ln K a

(2.18)

V kontextu uvedeného lze standardní reakční Gibbsovu energii ΔGor vymezit jako změnu volné energie, která nastane při chemické reakci složek, jejichž počáteční aktivity jsou rovny jedné (resp. jsou takové, aby zlomek v (2.18) byl roven jedné). Pro doplnění lze dále uvést, že standardní slučovací Gibbsovy energie sloučenin - tj. jednotlivých reakčních složek ΔGosl, jež slouží k bilančnímu odvození ΔGor, jsou definovány jako Gibbsovy energie (vztažené na 1 mol základních reakčních obratů) pro reakce, jimiž by dané sloučeniny vznikly přímo z prvků. Výchozím referenčním stavem jsou tedy prvky, jejichž ΔGosl je logicky nulová. Logaritmický vztah (2.18) je samozřejmě platný i pro procesy disociační. Pro disociaci vody H2O ↔ H+ + OHplatí:

22

Gr0   RT ln

aH  aOH  a H 2O

  RT ln K v

(2.19)

Kv označuje disociační konstantu vody. Míra disociace čisté vody je velmi nízká (za standardních podmínek se aktivita - koncentrace H+ a OH- blíží hodnotě 10-7 mol.l-1). Aktivita nedisociované molekuly je tedy prakticky maximální, tj. = 1. Pro disociační konstantu vody je změřeno: K v  aH  aOH   1014 Konstantní hodnota tohoto iontového součinu je samozřejmě platná i pro veškeré vodné roztoky v rovnováze za standardních podmínek. Z uvedených skutečností lze nahlédnout, že vhodným parametrem pro bilanční vyhodnocování chemických a fyzikálně-chemických dějů je vždy funkce charakterizující stav dané soustavy, který je logaritmicky závislý na parametru vymezeném relativním množstvím (molárním) každé ze složek. V této souvislosti a také z praktických důvodů daných naměřenou hodnotou disociační konstanty vody (celá záporná mocnina čísla 10), byl zaveden parametr pH = -log aH+, který charakterizuje stav vodných roztoků z hlediska míry kyselosti resp. zásaditosti. Pro mírně kyselé či mírně alkalické prostředí lze aktivitu H+ nahradit koncentrací - lze psát: pH = -log [H+]. Pro popis stavových vlastností alkalických roztoků bývá často používán obdobný parametr pOH = -log aOH-, příp. pOH = -log [OH-].

3. Vybrané aspekty chování iontů v roztocích 3.1.

Aktivta a aktivitní koeficient

Definici ideálního roztoku, jehož fyzikálně-chemický stav je determinován rovnicemi termodynamických rovnováh (např. ΔG0 = -RT ln K) zahrnujícími přímo koncentrace rozpuštěných látek, se reálně blíží pouze zředěné roztoky neelektrolytů, kde se neuplatňují coulombické interakce. Roztoky elektrolytů, tj. látek ve vodě disociovaných na ionty, jsou roztoky neideálními. Rovnice termodynamických rovnováh zde platí pro aktivity složek (iontů) těchto neideálních roztoků. Aktivita a vyjadřuje reálné (s termodynamickými rovnicemi konzistentní) působení rozpuštěné látky i (iontu) odpovídající její koncentraci korigované zahrnutím „stínících“ coulombických interakcí vyjádřených aktivitním koeficientem γ.

a i  ci  i

(3.1)

Při nekonečném zředění roztoku elektrolytu je γ → 1 a roztok se blíží ideálnímu. Relace aktivitního koeficientu a Gibbsovy energie (tj. stavové veličiny) vyplývá z porovnání chemického potenciálu i-té složky ideálního a neideálního roztoku.

iid  i0  RT ln ci i  i0  RT ln ci  RT ln  i

(3.2)

23

Člen RT ln γi odpovídá změně molární Gibbsovy energie ΔG i-té složky (iontu) roztoku při přechodu z ideálního do neideálního stavu (ΔG rovnovážného procesu je u roztoků ekvivalentní chemickému potenciálu násobenému stechiometrickým koeficientem), γi je individuální iontový aktivitní koeficient. Z měření příp. výpočtu lze s ohledem na elektroneutralitu (vždy jsou samozřejmě přítomny oba ionty) obdržet střední aktivitní koeficient γ± daného elektrolytu. Změna Gibbsovy energie pro přechod 1 molu obecného elektrolytu

K xA y   xK zK  yA zA z ideálního do neideálního stavu - tzv. dodatková Gibbsova energie je vyjádřena:

resp.

G  RT x ln    y ln   

(3.3)

G  RT x  y ln  

(3.4)

Z porovnání (3.3) a (3.4) plyne vztah pro střední aktivitní koeficient elektrolytu:

    x . y 

1  x y 

(3.5)

Střední aktivita i molární koncentrace jsou zavedeny obdobně jako střední aktivitní koeficient: a  x  y ax .ay  c  

(3.6)

c  x  y cx .cy  c  x  y x x . y y

(3.7)

Ve zředěnějších roztocích nezávisí střední aktivitní koeficient silného elektrolytu na druhu iontů, ale pouze na koncentraci ci a nábojovém čísle zi všech iontů v roztoku. V této souvislosti je pro vyjádření vlivu iontů v roztoku zavedena iontová síla Ic:

Ic 

1 ci zi2  2 i

(3.8)

Závislost mezi středním aktivitním koeficientem a iontovou silou je pro zředěné elektrolyty vyjádřena Debye-Hückelovým limitním zákonem:

 log    A zK zA I c

(3.9)

24

3.2

Debye-Hückelova teorie iontové atmosféry

Debye-Hückelova teorie odvozuje analytický vztah pro výpočet potenciálu iontové atmosféry φ’ (dodatkového potenciálu), na jehož základě lze vypočítat dodatkovou Gibbsovou energii roztoku elektrolytu, tj. φ’ je exaktním východiskem pro vyčíslení aktivitního koeficientu.

Obr. 3.1 Středový kladný ion obklopený iontovou atmosférou se záporných a kladných iontů. Všechny ionty mají poloměr a. Popsaný model (Debye Hückelova teorie) počítá dodatkový potenciál v místě středového iontu, resp. ve vzdálenosti a od středu. Primárním východiskem je Poissonova rovnice popisující vztah potenciálu φ k objemové hustotě náboje ρ pro sféricky symetrické uspořádání.

1 d  2 d   r  2 r dr  dr   0 r

(3.10)

Objemová hustota náboje v okolí centrálního iontu samozřejmě není pro řešenou situaci konstantou - je závislá na potenciálu všech iontů, jenž v daném místě svým dosahem figurují. Je tedy nutno nalézt příslušnou závislost ρ = f(φ). Použitelným východiskem je Boltzmannův teorém (3.11) vymezující počet částic i-tého druhu C’i s energií o hodnotu E větší než je energie průměrná, která odpovídá tepelnému pohybu kT (pokud E = 0, C’i = Ci). Ci  Ci e  E kT

(3.11)

Přičemž Ci je střední počet iontů i-tého druhu v objemové jednotce roztoku. Aplikace (3.11) na řešenou situaci předpokládá, že energie E může odpovídat práci potřebné pro přivedení iontu o náboji Q do místa potenciálu φ, tj. součinu Qφ. Ci  Ci e Qi

kT

(3.12)

Vztah (3.12) tedy vymezuje počet iontů i-tého druhu (vztažený na objemovou jednotku) působících svým potenciálem v okolí centrálního iontu ohraničeném iontovou atmosférou. Objemová hustota náboje je pak rovna celkovému náboji všech působících iontů (náboj je sumován přes všechny druhy iontů).

   CiQi   CiQi eQ  kT i

(3.13)

Po dosazení (3.13) do (3.10) vychází:

25

1 d  2 d  1 r  2 r dr  dr   0 r

C Q e i

Qi kT

i

(3.14)

Řešením této sféricky symetrické diferenciální rovnice lze obdržet hledanou funkci pro potenciál středového iontu a iontové atmosféry φ = f(r). Pro zachování možnosti nalezení analytického řešení, o což v popisovaném modelu jde, je však nutno u rovnice (3.14), která je nelineární, zavést jistá zjednodušení. Exponenciální faktor na pravé straně (3.14) lze, s ohledem na relaci Qiφ << kT platnou pro zředěné roztoky, aproximovat Taylorovou řadou. Q 1  Q   1 i   i   kT 2!  kT  2

e

Qi kT

(3.15)

Při zanedbání všech vyšších členů řady (3.15) mimo prvních dvou vychází pro objemovou hustotu náboje lineární závislost na potenciálu (se zápornou směrnicí).

   Ci Qi 

 kT

C Q i

2 i

(3.16)

První člen je navíc v souvislosti s elektroneutralitou nulový. Vyjádření náboje iontu Qi prostřednictvím nábojového čísla zi × elektronvolt e pak vede k jednoduchému vztahu:

 

e 2  Ci zi2 kT

(3.17)

Dosazením zjednodušené lineární závislosti objemové hustoty náboje na potenciálu (3.17) do (3.10) je obdržena linearizovaná Poissonova-Boltzmannova rovnice:

1 d  2 d  e 2 r  Ci z i2    2 r dr  dr   0 r kT

(3.18)

Parametry, jenž jsou z hlediska hledaného řešení Poissonovy-Boltzmannovy rovnice konstantní, lze shrnout jako parametr b. d  2 d  2 2 r b r  dr  dr 

(3.19)

e2 Ci zi2   0 r kT

(3.20)

b2 

Substitucí u = r φ přejde rovnice (3.19) na tvar:

d 2u  b 2u 2 dr

(3.21)

Neboť:

du dr d ,   r dr dr dr

d 1  du      dr r  dr 

Pravá strana (3.19) je tedy:

26

d  2 d  d  du d2u d  du  r    r 2   r  r   r dr  dr  dr  dr dr dr  dr du d2u du d2u   r 2     r 2 dr dr dr dr Obecná forma (3.21) umožňuje hledat její řešení ve tvaru obecně vyhovujícím dané rovnosti: u  Ae -br  Be br

(3.22)

Pro původní (3.19) je pak k dispozici tvar:



A -br B br e  e r r

(3.23)

Integrační konstanty A a B vyplynou z okrajových podmínek. Vyvození B je triviální: Pro r → ∞ musí platit φ = 0, tj. 0 

A - B  e  e , což platí pouze pro B = 0. Použitelné řešení  

je tedy ve tvaru:



A -br e r

(3.24)

Výpočet A je založen na vyjádření objemové hustoty náboje ρ = f(r). Po dosazení (3.24) jako vypočítané funkce závislosti potenciálu φ na r do (3.17) a při uvážení (3.20) vychází pro objemovou hustotu náboje jako funkci r:

 

Ab 2 0 r br e r

(3.25)

Vztah (3.25) reflektuje objemovou hustotu náboje v uvažovaném iontovém oblaku (v závislosti na r). Celkový náboj iontového oblaku, tj. integrál přes sférické povrchy obklopující centrální iont (viz obr. 3.1) v intervalu vzdáleností od a do nekonečna, pak musí být až na znaménko roven náboji středového iontu: 

 4r

2

 (r )dr   zi e

(3.26)

a

Po dosazení výrazu pro nábojovou hustotu (3.25) 

Ab 2 0  r  4r 2 e br dr  z i e

(3.27)

a

a integraci je z této modelové úvahy obdržena hodnota konstanty A.

A

eba 4 0 r 1  ba zi e

(3.28)

Pro výsledný potenciál tedy z (3.24) a (3.28) vychází:



eba e br 4 0 r 1  ba r zi e

(3.29)

27

Vypočítaný potenciál φ zahrnuje coulombický příspěvek samotného středového iontu zi e a příspěvek iontové atmosféry φ’ . Z relace φ = φc + φ’ vyplývá pro φ’: c  4 0 r r

 

 eba br   e  1 4 0 r r  1  ba  zi e

(3.30)

Potenciál iontové atmosféry v místě středového (stíněného) iontu vychází po dosazení r = a (a je nejmenší možná vzdálenost odpovídající poloměru středového iontu).

 ra 

 zi e b 4 0 r 1  ba

(3.31)

Pro extrémně zředěné roztoky je ba << 1 a výraz (3.31) lze zjednodušit:

 ra 

 zi eb 4 0 r

(3.32)

Vypočítaný dodatkový potenciál φ’ determinovaný iontovou atmosférou koresponduje s dodatkovou Gibbsovou energií roztoku elektrolytu, tj. je v relaci s aktivitním koeficientem. Dodatková ΔG pro nabití iontu „stínícím“ nábojem odpovídajícím vlivu iontové atmosféry je vyjádřena: ze

G    dQ  o

 bQ z 2 e 2b d Q   o 4 0 r 8 0 r

ze

(3.33)

Pro dodatkovou Gibbsovu energii na jeden iont platí:

G  kT ln  i

(3.34)

Vztah (3.34) principielně koresponduje s (3.4), který je však vztažen k molárnímu množství. Sloučením (3.33) a (3.34) vychází:

ln  i  

z 2 e 2b 8 0 r kT

(3.35)

Vztah (3.20) pro parametr b lze převést do formy pro molární koncentrace zastoupených iontů (platí R = NA k, ci = Ci /NA).  N 2e2 b   A   0 r RT

12

  c z   2 i i

B I

(3.36)

Sloučením (3.35) a (3.36) a dalšími úpravami vznikne standardní tvar Debye-Hückelova limitního zákona (3.9) popisující závislost středního aktivitního koeficientu na iontové síle a nábojových číslech kationu a anionu elektrolytu.

 log    A z K z A

Ic

Odvozený parametr A obsahuje pouze obecné konstanty – v rámci uvažovaného zjednodušení použitelného pro extrémně zředěné roztoky není aktivitní koeficient závislý na poloměrech přítomných iontů.

28

4. Oxidačně – redukční procesy, vztah Gibbsovy energie a oxidačně-redukčního potenciálu Standardní reakční Gibbsova (volná) energie ΔGor je termodynamickým potenciálem průběhu chemických dějů. U oxidačně-redukčních reakcí spojených s výměnou elektronů, se potenciál průběhu děje projevuje ve formě elektrického napětí, které je měřitelné při specifickém uspořádání jako elektrochemický článek. Např. oxidačně-redukční reakci Zno + CuSO4 → Cuo + ZnSO4 probíhající zleva doprava, lze uspořádat tak, že její průběh znamenající výměnu elektronů bude limitován (velkým) odporem na voltmetru připojeném na elektrodách z předmětných kovů - viz obr 3.2. Tohoto efektu je dosaženo tím, že probíhající dílčí oxidačně-redukční děje Zn2+/Zno a Cu2+/Cuo jsou odděleny polopropustnou přepážkou, tj. reakce může probíhat pouze prostřednictvím elektronů přenášených vnějším vodičem s voltmetrem. Daný elektrochemický článek je popsán jako článek Danielův.

Obr. 4.1 Danielův článek Elektrické napětí má rozměr energie – pro vztah mezi oxidačně-redukčním potenciálem a Gibbsovou energií jako obecným termodynamickým potenciálem platí: ∆G = -nfE

(4.1)

Budou-li tedy aktivity Cu2+ a Zn2+ v prezentovaném Danielově článku rovny jedné (resp. bude-li jejich poměr roven jedné; aktivity kovů Zno a Cuo tvořících elektrody jsou maximální, tj. rovny jedné), bude napětí naměřené na voltmetru (s nekonečným odporem) odpovídat standardní Gibbsově energii dané reakce ΔGor vydělené Faradayovou konstantou a dvojkou (2 vyměňované elektrony). Pro elektrické napětí, jež můžeme popsaným způsobem pro oxidačně-redukční reakci naměřit, tj. okamžitý redox potenciál této reakce platí stejná principiální východiska jako pro Gibbsovu energii ΔG. Tj. v počátečním stavu, kterému odpovídá jedničková aktivita reakčních složek, je E = -ΔGor/nf a v rovnováze (konečný stav, který se ustálí, jsou-li

29

elektrody propojeny vodičem) je redox potenciál reakce roven nule - článek je vybit. Tj. pro obecnou oxidačně-redukční reakci s m reaktanty a l produkty platí Nernstova rovnice: l 1 RT a prod ,1  a prod,l EE  ln 1 m nf areakt,1  areakt ,m

(4.2)

Gr0 RT  ln K a nf nf

(4.3)

0

Resp.:

E0  

Vyhodnocování průběhu oxidačně-redukčních dějů se provádí bilancováním potenciálů dílčích poločlánků, tj. elektrodových reakcí. V rámci výše uvedené reakce probíhají dva dílčí děje Zn2+ + 2e- ↔ Zn

-0,7628 V

Cu2+ + 2e- ↔ Cu

+0,337 V

Vpravo jsou uvedeny standardní elektrodové potenciály Eo obou elektrodových reakcí. Hodnoty Eo jsou zjištěny měřením jednotlivých poločlánků (elektroda z příslušného kovu ponořená do roztoku jeho kationů o jednotkové aktivitě) oproti standardní vodíkové elektrodě, jejíž standardní elektrodový potenciál je vzat za rovný nule. Elektrodové reakce jsou konvenčně psány ve směru redukce. Zinek má tendenci přecházet do oxidované formy (do roztoku), ΔGor elektrodové reakce pro tento prvek (zapsané ve směru redukce) je kladná - reakce v daném konvenčně určeném směru probíhat nebude. Standardní elektrodový potenciál Eo je pro tuto reakci dle (4.3) záporný. Měď je ušlechtilý kov - nemá tendenci se oxidovat (rozpouštět), tj. předmětná elektrodová reakce bude v daném směru probíhat - ΔGor < 0, Eo > 0. Oxidačně-redukční reakce probíhající v Danielově článku zahrnuje elektrodovou reakci zinku v opačném směru. V bilančním výpočtu standardního elektromotorického napětí tohoto článku bude tedy zahrnut EoZn+/Zn s opačným znaménkem (+0,7628 + 0,337 = 1,0998 V). Uvedená hodnota odpovídá napětí, jež bude naměřeno, bude-li aktivita Zn2+ a Cu2+ rovna jedné.

Pro jednotlivé elektrodové reakce lze obecně psát:

E  E0 

  RT ared1 ,1  aredl ,l ln 1 nf aox,1  aoxm,m

(4.4)

Převážná většina elektrodových reakcí (pro kovy všechny) zahrnuje v reakční směsi pouze dvě složky - redukovanou a oxidovanou formu, přičemž stechiometrické koeficienty obou jsou rovny jedné. Pro tento typ dějů lze Nernstovu rovnici psát ve zjednodušené podobě: E  E0 

a RT ln ox nf ared

(4.5)

Resp. tvoří-li redukovanou formu kov (tj. poločlánek s elektrodou z příslušného kovu), jehož aktivita je maximální - rovna jedné, lze psát:

E  E0 

RT ln aox nf

(4.6)

30

Oxidačně - redukční potenciál prostředí - Eh Oxidačně-redukční procesy v prostředí (přírodní vody) jsou reprezentovány složitým systémem tvořeným větším počtem jednotlivých párových iontů (Fe3+/Fe2+, Mn4+/Mn2+, PO43/PO33- ap.) - poločlánků. Pro takovéto prostředí je zaveden celkový - sumárním oxidačněredukční potenciál označovaný jako Eh. Hodnota Eh by byla vyčíslitelná jako sumární ze všech dílčích elektrodových potenciálů vypočtených pro jednotlivé párové ionty při dosazení aktuálně zjištěných koncentrací těchto iontů do uvedené rovnice. Zjišťování Eh libovolného vodného prostředí se prakticky provádí měřením - elektrodou z lesklé platiny propojenou s vhodnou referenční. Zpravidla bývá naměřena nenulová hodnota, Eh > 0 reprezentuje oxidační vlastnosti a naopak. Přírodní (vodné) prostředí je zvláště u netekoucích vod považováno za ustálené v rovnováze. Pro rovnovážné podmínky je však obecně ΔG = 0, Eh = 0 (systém v rovnováze nemá potenciál – je „vybitý“). Tento princip (ΔGr = 0, Erovn = 0) se vztahuje k obecnému náhledu chemického děje, kdy jsou počáteční aktivity reakčních složek uvažovány jako jedničkové, příp. mající s ohledem na stechiometrii takové hodnoty, aby zlomek na levé straně (4.4) byl roven jedné. Pro výchozí stav teoreticky platí: ΔG = ΔGo, E = Eo, „Eh = Eho“. Ve vodných systémech v přirozeném prostředí toto obecné východisko samozřejmě dodrženo není a nulová hodnota Erovn nemůže být pro stav, kdy zlomek na levé straně (4.4) odpovídá rovnováze, dosažena. Naměřené Eh > 0 reálného systému zpravidla odpovídá koncentračnímu přebytku rozpuštěného kyslíku. Systém je v energetickém minimu, do obecného rovnovážného stavu ΔG = 0, tj. Eh = 0 však nemůže v souvislosti s koncentračním nepoměrem oxidovatelných složek vůči rozpuštěnému O2 dojít.

5. Kinetika fyzikálně-chemických procesů v prostředí 5.1

Rychlost, rozsah a řád procesů

Kinetika procesů v prostředí, zjm. pak chemických reakcí studuje průběh těchto termodynamicky možných dějů v čase, tj. jejich rychlost. Obecně může v soustavě probíhat jeden izolovaný proces nebo více procesů simultánních, jako jsou děje zvratné, boční, následné. Dále jsou rozlišovány procesy (zjm. chemické reakce) homogenní, kdy jsou všechny složky systému v jedné fázi a procesy heterogenní, např. řízené difúzí na rozhraní fází. Předmětné atributy – rychlost, řád, rozsah jsou primárně zavedeny pro chemické reakce – v dalším textu je tedy dodržována terminologie vztahující se k chemickým dějům. Nicméně popisované principy mají obecnější význam a korespondují i s jinými typy dynamických procesů v prostředí (jaderné, biologické, fyzikálně-chemické ap.). Pro obecnou reakci: αA + βB → ηY + ωZ je definován diferenciál rozsahu reakce – viz (2.13): 31

d  

dnA





dnB





dnY





dnZ



(5.1)

α, β, η, ω jsou stechiometrické koeficienty, nA, nB, nX, nZ jsou látková množství reakčních složek. Pro infinitezimální časovou změnu rozsahu reakce, tj. reakční rychlost w platí:

w

d 1 dnA 1 dnB 1 dnY 1 dnZ     dt  dt  dt  dt  dt

(5.2)

Vztažení w k jednotkovému objemu vede k definici rychlosti standardní reakce v roztoku v.

vwV 

d 1 dc A 1 dcB 1 dcY 1 dcZ     V dt  dt  dt  dt  dt

(5.3)

Okamžitá koncentrace reagující látky souvisí s její počáteční koncentrací vztahem:

ci  ci0   i x

(5.4)

kde x   V je rozsah reakce vztažený na objemovou jednotku a νi je stechiometrický koeficient látky i. Stechiometrické koeficienty jsou pro výchozí látky brány jako záporné, pro produkty jako kladné. Tj.:

c A  c A0  x, c B  c B0  x, cY  cY0  x, cZ  cZ0  x (5.5) Pro rychlost reakce (probíhající v konstantním objemu) pak platí:

v

1 dci dx   i dt dt

(5.6)

dx 1 dc A 1 dcB 1 dcY 1 dcZ     dt  dt  dt  dt  dt Rychlost reakcí je závislá na teplotě a koncentraci reagujících látek. Za izotermických podmínek platí:

v  kcAr cBs

(5.7)

Reakční rychlost je úměrná mocninám okamžitých koncentrací reagujících látek. Důležitým pojmem je řád reakce, který je definován jako součet exponentů okamžitých koncentrací reagujících látek zjištěný měřením reakční rychlosti. Může být zdůrazněno, že exponenty determinující reakční rychlost jsou v obecnosti nezávislé na stechiometrii.

5.2

Procesy nultého a prvního řádu

U procesů nultého řádu je rychlost v průběhu děje konstantní nezávisí na koncentraci reagujících látek. Neuplatňuje se zde princip (5.7), resp. mocniny okamžitých koncentrací reagujících látek v (5.7) jsou rovny nule. v

dc A k dt

(5.8)

32

Pro časový průběh koncentrace (obecně relativního množství) spotřebovávané (výchozí) látky A tedy platí:

c A  c A0  kt

(5.9)

Procesy prvního řádu probíhají rychlostí závisející na okamžité koncentraci (relativním množství) spotřebovávané (výchozí) složky. Jde o rychlost odbourávání této složky A, tj. se záporným znaménkem. k A  Y

dc A  kcA dt



(5.10)

Nejjednodušší možná obyčejná diferenciální rovnice (5.10) je snadno integrovatelná po separaci proměnných. dc A    kt cA

(5.11)

ln c A  kt  C

(5.12)



Integrační konstanta je jednoduše vyčíslena v bodě t = 0. V počátku procesu platí:

C  ln c A0

(5.13)

A výsledné řešení je ve tvaru:

ln

c A0  kt cA

(5.14)

Resp.:

c A  c A0 exp  kt

(5.15)

Pro okamžitou koncentraci výchozí složky A vyjádřenou obecně pomocí (5.4) je:

c A0  kt c A0  x



(5.16)

x  c A0 1  exp  kt

(5.17)

ln



Vedle časového průběhu koncentrace výchozí složky A je samozřejmě vyjádřitelná i funkce pro koncentraci produktu Y.

cY  cA0 1  exp  kt  cY0

(5.18)

Dalším parametrem, který bývá pro kinetické procesy často řešen, je poločas rozpadu, tj. čas, v kterém je koncentrace výchozí složky přesně poloviční oproti počáteční hodnotě:

c A  0,5c A0

(5.19)

Po dosazení (5.19) do (5.14) je pro poločas rozpadu obdrženo:

ln 2  t1 2 k

(5.20)

33

5.3

Procesy druhého řádu

Procesy druhého řádu pro dvě reagující složky k A B  produkty

Této kategorii kinetických procesů odpovídají jedničkové exponenty ve vztahu (5.7) a rychlost reakce závisí na okamžité koncentraci obou výchozích složek. S ohledem na dvě výchozí reagující složky je rychlost procesu vyjádřená obecně ve smyslu (5.6), tj.:





dx  k c A0  x c B0  x dt



(5.21)

Po separaci proměnných:

dx  kdt c  x cB0  x





0 A



(5.22)

Integrace levé strany (5.22) je podmíněna rozkladem na parciální zlomky.

c

1 0 A



x c x 0 B











M N M cB0  x  N c A0  x   c A0  x cB0  x c A0  x cB0  x









(5.23)

Z (5.23) plyne rovnost:



 

1  M cB0  x  N cA0  x



(5.24)

Vyčíslením (5.24) pro počáteční stav x = 0 vychází:

Mc B0  Nc A0  1

(5.25)

Pro libovolné x > 0 tedy musí s ohledem na (5.24) (roznásobené) a (5.25) platit: M N 0

(5.26)

Kombinací (5.25) a (5.26) je obdrženo: 1  M c  cB0

N





0 A

(5.27)

Tvar pravé strany (5.21) je po dosazení (5.27) do (5.23) již vhodný pro integraci



1 dx 1 dx  0  kdt 0 0 0 0 c  cB c A  x c A  cB cB  x 





0 A

 





(5.28)

Která (včetně pravé strany (5.28)) vychází po triviální substituci v tvaru:

 

 

 



1 1 ln c A0  x  0 ln cB0  x  kt  C 0 c  cB c A  cB0



0 A

(5.29)

Integrační konstanta je stejně jako u procesů prvního řádu jednoduše vyčíslena pomocí počátečního stavu t = 0. C

1 1 ln c A0  0 ln cB0 0 0 c  cB c A  cB



0 A







(5.30)

Nyní se lze od obecného rozsahu reakce (vztaženého k jednotce objemu) x opět vrátit ke koncentracím reagujících složek - c A0  x  c A , cB0  x  cB . Tj.:

34



 





1 1 ln c A  ln c A0  0 ln cB  ln cB0  kt 0 c  cB c A  cB0



0 A

1 c 1 c ln A0  0 ln B0  kt 0 0 c  cB c A c A  cB cB



0 A







(5.31)

Výsledný vztah zahrnující časový průběh koncentrací reagujících složek při procesu druhého řádu je tedy ve tvaru:

c B0 c A 1 ln  kt c A0  c B0 c A0 c B





(5.32)

Okamžité koncentrace reagujících složek samozřejmě nejsou nezávislé - viz (5.5). Procesy druhého řádu s jednou výchozí látkou Procesy druhého řádu se samozřejmě mohou týkat i odbourávání jedné složky, kdy je tento děj obecně podmíněn interakcí dvou částic této složky A. k 2A   produkty

Pro rychlost odbourávání složky A tedy platí: 

dc A  kcA2 dt

(5.33)

Po integraci a vyčíslení integrační konstanty (v bodě t = 0) je obdrženo: 1 1  kt  0 cA cA

c A0  c A  kt c A0 c A

(5.34)

Poločas rozpadu výchozí složky A vychází:

1  t1 2 c A0 k

5.4

Procesy n-tého řádu

Obecný kinetický vztah pro procesy libovolného řádu n = r + s + ..., který je uvažován jako proměnná může být platný pouze pro situaci, kdy výchozí koncentrace reagujících složek přesně odpovídají stechiometrii reakce. Pouze v tomto případě zůstávají poměry okamžitých koncentrací reagujících složek v čase konstantní. Např. pro reakci:

k αA  βB  γC   produkty

lze psát







dx c A0  x cB0  x cC0  x k dt   



(5.35)

35

Tj. obecný vztah pro proces n-tého řádu



dx  k c0  x dt



n

(5.36)

může platit pouze pro popisovanou podmínku pro počáteční koncentrace reagujících složek, tj.: c 0  c A0   cB0   cC0 

(5.37)

(5.36) je po triviální separaci proměnných jednoduše integrovatelný

 c

dx 0

x

 kdt  

(5.38)

n

Integrační konstanta je samozřejmě vyčíslena prostřednictvím počátečního stavu (t = 0).





 

1 n  1 1 c0  x  kt  c0 1 n 1 n

1 n 

(5.39)

Po úpravě je obdržen obecný vztah pro libovolný řád procesu.

1  1 1   n 1  0 n 1   kt 1  n  c c 

 

(5.40)

Experimentální určování řádu procesu Řád procesu by samozřejmě mohl být jednoduše vyčíslen z (5.40) při známém k a c0 po dosazení změřené koncentrace reagující složky v libovolném čase t v průběhu procesu (při známém k a c0). Resp. při změření koncentrace pro dvě různé počáteční hodnoty by byla determinována i rychlostní konstanta k. Při experimentálním určení řádu procesu tímto způsobem by počáteční koncentrace reagujících složek každopádně musely odpovídat podmínce (5.37) pro platnost (5.40). V této souvislosti se pro určování řádu procesu používá praktičtější způsob založený na vyčíslení koncentračních diferencí v počátku procesu t → 0 pro dvě různé c0. Výchozím principem je tedy (5.36), která je pro počátek děje použitelná i při eventuelním nedodržení (5.37).

 

n

dx dt  x t  k c 0

(5.41)

Jsou-li koncentrační úbytky libovolné složky pro počátek procesu Δx odečteny pro její dvě různé c0 je (při stejném Δt), lze obě kinetické rovnice (5.41) vydělit a obdržet výsledný vztah pro výpočet řádu procesu.

x 2 x 1  c20 5.5

c10



n

(5.42)

Simultánní procesy

V rámci simultánních procesů lze v obecnosti hovořit o zejména dějích (reakcích) bočných, zvratných a následných.

36

Bočné procesy

k1

a)

A k2

b)

c)

Y Z

k1 A  B  Y k2 A  C  Z

A B  Y A C  Z

Bočné procesy mohou probíhat jako rozvětvené – ad a), konkurenční – ad b), příp. nezávislé – ad c). Pro případ procesů rozvětvených platí:









dc dx  k1 c A0  x  k 2 c A0  x   A  c A k1  k 2  dt dt

(5.43)

Pro koncentraci výchozí látky A tedy vychází:

 ln c A  k1  k2 t  ln c A0 ln

c A0  k1  k 2 t cA

(5.44)

Rychlost růstu koncentrací produktů je v rámci uvažovaného prvního řádu obou procesů úměrná cA.

dcY  k1c A dt dc Z  k2c A dt

(5.45)

Při uvážení (5.44):

dcY  k1c A0 e -k1  k2 t dt

(5.46)

dcZ  k 2c A0 e -k1  k2 t dt

Vyčíslené vztahy (5.46) již nejsou diferenciálními rovnicemi, ale pouze rovnostmi nezávislých derivací. Po integraci a vyčíslení integrační konstanty je obdrženo:

cY 



k1c A0 1 - e -k1  k2 t k1  k 2



k 2 c A0 cZ  1 - e -k1  k 2 t k1  k 2

 

(5.47)

37

Přičemž integrační konstanta ve tvaru

C

k1c A0 k1  k 2

(5.48)

je vyčíslena pro počáteční podmínky cY, cZ = 0, při t = 0.

Obr. 5.1 Možný trend koncentrací složek rozvětveného procesu při k2 > k1 Struktura dalších možných typů bočných procesů - ad b), c) víceméně kopíruje skutečnosti již uvedené v rámci této stati 5. Zvratné procesy k1

A

B

k2

Pro oba směry procesu je uvažován první řád. Pro rovnovážný stav pak není předpokládán výrazný přebytek produktu B či výchozí složky A. Rychlostní rovnice je logicky ve tvaru:







dx  k1 c A0  x  k 2 c B0  x dt



(5.49)

Definice zvratného procesu prakticky koresponduje s principem definice rovnovážného stavu, který lze využít pro vyčíslení jedné z rychlostních konstant pomocí konstanty rovnovážné, která je často k dispozici. Pro rovnováhu platí:









k1 c A0  x rov  k2 cB0  x rov

(5.50)

S ohledem na (2.3):

c c

0 B 0 A

x x

 

rov

rov



k1 K k2

(5.51)

Dosazení (2.3), resp. (5.51) do (5.49) umožňuje obdržení tvaru vhodného pro jednoduchou integraci.

 Kc 0  c B0  dx  k 2 Kc A0  Kx  c B0  x  k 2 K  1 A  x  dt  K  1 





(5.52)

38

dx  k 2 K  1dt  Kc  cB0    x   K  1 

(5.53)

0 A

Po integraci a vyčíslení integrační konstanty (pro t = 0, x = 0):

 Kc 0  c 0   Kc 0  c 0   ln  A B  x   k 2 K  1t  ln  A B   K  1   K  1   K  1 x   k K  1t  ln 1  0 0  2  Kc A  cB 

(5.54)

Výsledný vztah pro x = f(t) je ve tvaru: x

Kc A0  c B0 1  e k2  K 1t  K  1





(5.55)

Obr. 5.2 Možný trend koncentrací složek zvratného procesu.

Následné procesy k1 k2 A  B  Y

Odbourávání výchozí složky A za vzniku složky B je uvažováno jako proces prvního řádu, tj.:



dc A  k1c A , c A  c A0 e -k1t dt

Pro časovou změnu cB tedy platí: 

dc B  k 2 c B  k1c A  k 2 c B  k1c A0 e -k1t dt

(5.56)

Vznik konečného produktu Y je taktéž brán jako jednoduchý proces prvního řádu, musí tedy platit:

dcY  k 2 cB dt

(5.57)

Rovnice (5.56) je dobře analyticky řešitelná, byť zde nelze provést jednoduchou separaci proměnných jako v předchozích případech. K řešení vede použití substituce cB ≡ u∙v.

39

dc B   uv   u v  uv  dt

Z porovnání s (5.56) plyne:

u v  k 2 c B , uv  k1c A0 e -k1t u v  k 2 uv , u   k 2 u ,

u  C1e  k2t

v 

k1c A0 e -k1t k1c A0 e -k1t k1c A0 k2-k1 t   e u C1 C1e -k 2t

v

k1c A0 k2-k1 t k1c A0 e d t  e k2-k1 t  C 2  C1 C1 k 2 -k1 

Výsledná funkce je tedy ve tvaru:

 c0 k  uv  c B  e k2t  A 1 e k2 k1 t  C   k 2  k1 

(5.58)

Integrační konstanta C = C1C2 je samozřejmě obdržena pro t = 0.

C

c A0 k1 k 2  k1

Pro časový průběh koncentrace složky B tedy platí:

cB 



c A0 k1 k1t e  e  k 2t k2  k1



(5.59)

Vztah (5.59) je pak dosazen do (5.57) a pro nárůst koncentrace výsledného produktu Y lze psát:



dcY c0 k  k2 A 1 e k1t  e k2t dt k2  k1



(5.60)

Rychlostní rovnice (5.60) je stejného typu jako (5.46) a je snadno integrovatelná.

cY  k 2

c A0 k1  1 k2t 1 k1t   e  e   C k 2  k1  k 2 k1 

C  c A0 při t = 0 tj:   k2 k1 cY  c A0 1  e k1t  e k2t  k 2  k1  k 2  k1 

(5.61)

40

Obr. 5.3 Možný trend koncentrací složek následného procesu.

5.6

Příklad systému kinetických procesů

Chemické, fyzikálně-chemické, biologické ap. procesy, které probíhají v přírodním prostředí nebo jsou experimentálně zkoumány v laboratorních podmínkách mají většinou složitější strukturu než prezentované případy, resp. mohou být popsatelné jako ucelený systém dílčích kinetických dějů. Jediným nezávislým parametrem však stále zůstává čas, tj. tyto systémy jsou označovány jako dynamické. Příkladů možných dynamických systémů popisujících různé procesy, resp. jejich provázané struktury je samozřejmě bezpočet. Pro ilustraci lze zmínit systém kinetických procesů přibližně simulující možné děje probíhající při vyluhování sloučenin arsenu z povrchu částic elektrárenského popílku působením zředěné kyseliny dusičné – viz prezentované modelové schéma. Procesy znázorněné jednosměrnými šipkami jsou považovány za nevratné s jednoznačným průběhem v naznačeném směru. Procesy sorpční jsou modelovány jako vratné. Řídícím parametrem systému je pH resp. koncentrace iontů H3O+ které vstupují do sorpčních procesů a do procesů rozrušování (rozpouštění) povrchu popílku. Hodnota pH je v modelu určována počáteční koncentrací HNO3 a koncentrací iontů kovů, jenž přešly do roztoku v důsledku rozrušování povrchu popílku. Při výpočtu výsledného pH je samozřejmě brána do úvahy také hydrolýza iontů uvažovaných kovů.

41

Prezentovanému schématu odpovídá soustava šesti diferenciálních a dvou bilančních rovnic. Je-li časový nárůst či úbytek vyjádřený levou stranou dle schématu závislý na dvou parametrech, je obecně uvažováno s druhým řádem vyjádřeným součinem těchto parametrů na pravé straně rovnic.

Použité symboly: inert počet nehydratovaných (z hlediska sorpce inertních) sorpčních míst  celkový počet sorpčních míst na daném povrchu n(Ass) množství As volně rozptýleného ve formě As2O3 na povrchu popílku, po smočení As rychle přechází do roztoku fyzikálním rozpouštěním n0(Ass) celkové množství As v systému, tj. počáteční množství As na povrchu popílku sorp počet hydratovaných volných sorpčních míst (AsIII), (AsV) množství sorbovaného AsIII, AsV, tj. počet sorpčních míst blokovaných As III, AsV c(AsIII), c(AsV) koncentrace AsIII, AsV v roztoku c(Me) koncentrace sumárního kladného náboje iontů kovů v roztoku, rychlost nárůstu tohoto parametru není limitována množstvím uvažovaných kovů na povrchu popílku (modelováno jako neměnné) rychlostní konstanty: kh hydratace povrchu popílku kfr fyzikální rozpouštění As2O3 po smočení popílku kir, k2ir přechod sorbovaného AsIII, AsV do roztoku v důsledku rozpouštění (rozrušování) povrchu ks, k2s sorpce AsIII, AsV na volná soprční místa povrchu popílku kd, k2d přechod sorbovaného AsIII, AsV do roztoku v důsledku desorpce kox oxidace AsIII na AsV, probíhá pouze v roztoku kmr přechod uvažovaných iontů kovů do roztoku v důsledku rozrušování povrchu

42

Dynamický systém tvořený prezentovanou soustavou obyčejných diferenciálních rovnic je pouze výchozí formulací popisující daný problém. Výsledné časové průběhy jednotlivých parametrů v závislosti na zvolených počátečních podmínkách (jejich počátečních stavech) však jsou (na rozdíl od obecných základních případů uvedených v předchozím) obdrženy numerickým výpočtem.

Obr. 5.4 Modelovaný časový průběh koncentrace As v roztoku v porovnání s naměřenými hodnotami. Těchto výsledků je dosaženo při nastavení rychlostních konstant rozrušování povrchu kir a sorpce ks v poměru 1 : 650.

Tyto výsledky jsou tedy obdrženy v číselné podobě s únosnou (resp. malou) odchylkou od přesného řešení, které by teoreticky mohlo být obdrženo analyticky. Pro daný případ, kdy jsou na pravých stranách rovnic obsaženy součiny počítaných závislých proměnných, tj. nelinearita, je analytické řešení (které je bezpochyby vždy elegantnějším) problematické, pokud vůbec existuje. Míra aproximace, tj. velikost odchylky numerického řešení od přesného obecně souvisí s použitou numerickou metodou. Smyslem prezentovaného příkladu není podrobná analýza fyzikálních a matematických aspektů daného problému. Např. popis sorpce jako heterogenního procesu pouze prostřednictvím rychlostní konstanty je až neúnosně zjednodušující. Popsaný problém však ilustruje elementární provázanost předpokládaného chemického (či fyzikálně-chemického) schématu probíhajících dějů a příslušejícího dynamického systému tvořeného soustavou obyčejných diferenciálních rovnic. Tento přístup k popisu časově závislých procesů je označován jako black box model. Časté jsou jeho aplikace v oborech zkoumajících procesy v přírodním prostředí, kdy jsou k dispozici jen dílčí měření či pozorování. V těchto případech může formulovaný dynamický systém přinést vedle informací o možné struktuře probíhajících procesů i chybějící neměřitelná data.

43

6. Nelineární dynamické systémy – úvodní poznámky Oscilační Belousov-Žabotinského chemická reakce V rámci výzkumu chemických či biologických dynamických procesů bylo u některých specifických případů pozorováno chaotické chování vyznačující se časovou oscilací některých parametrů. Významným přelomem v oblasti chemické kinetiky byl objev oxidačně-redukční reakční soustavy vykazující průběžnou koncentrační oscilaci některých reagujících složek. Tzv. oscilační reakce nazvaná po svých objevitelích jako reakce Belousov-Žabotinského je velmi složitá. V jedné z modifikací této reakce lze odhadovat následující hlavní děje: CH2(COOH)2 + 6 Ce4+ + 2 H2O → 2 CO2 + HCOOH + 6 Ce3+ + 6 H+ 10 Ce3+ + 2 HBrO3 + 10 H+ → 10 Ce4+ + Br2 + 6 H2O CH2(COOH)2 + Br2 → CHBr(COOH)2 + HBr Zjištěný oscilační průběh obecně nekorespondoval s předpokladem, že koncentrace reagujících složek chemického systému by měly spět do rovnovážného ustáleného stavu v rámci jednoznačně definovatelných trendů, tak jak je diskutováno v předchozí kapitole. Obecné termodynamické principy samozřejmě platí i zde. Rovnovážný stav spojený s ukončením oscilace Ce3+↔Ce4+ nastává okamžiku, kdy jsou převážně spotřebovány výchozí látky. Oscilační chování tohoto chemického systému odpovídá stavu vzdálenému od rovnováhy. Lze říci, že časový průběh koncentrací některých složek systému je na cestě k rovnováze chaotický. Model Brusselator V souvislosti se snahou o objasnění experimentálně zjištěného oscilačního (či chaotického) charakteru chování popsaného chemického systému byl navržen zjednodušený kinetický model, resp. dynamický systém: I.

A  X

II.

k1 2X  Y  3X

III.

k2 B  X  YC

IV.

X  D

Jsou-li složky A, B brány jako konstanty, čemuž by reálně odpovídal jejich velký přebytek (tj. kinetický model se týká stavu vzdáleného od rovnováhy), lze uvedené schéma vyjádřit dvěma diferenciálními rovnicemi.

  1  k  1X  k X 2 Y X 2 1

(6.1)

  k X  k X2Y Y 2 1

(6.2)

Pro rychlost procesu I. A IV. je tedy zjednodušeně počítáno s jedničkovou konstantou a stejně tak konstantní koncentrace složek A a B jsou brány jako jedničkové. Časové derivace na levých stranách jsou vyjádřeny zkráceně - symbolikou, která se nejčastěji užívá. Symboly označující složky systému X, Y, A, B, C, D mají v rovnicích význam přímo jejich koncentrací, příp. obecněji relativních množství. Tento model je dle místa svého vzniku označován jako Brusselator. Obě rovnice obsahují nelineární člen a analytické řešení této soustavy se nenabízí. Numerické řešení není 44

problémem, výsledkem je oscilační průběh parametrů X a Y vskutku korespondující s experimentálními poznatky.

Obr. 6.1 Výsledný oscilační průběh proměnných X, Y obdržený numerickým řešením modelu Brusselator (6.1), (6.2). Podstatou tohoto složitého dynamického chování je nelinearita řešené soustavy obyčejných diferenciálních rovnic. Modely tohoto typu jsou tedy označovány jako nelineární dynamické systémy a jejich řešení vykazuje oscilační (resp. chaotický) charakter. Neobsahují-li soustavy rovnic žádný parametr s dopředu determinovaným časovým průběhem, jsou tyto systémy označovány jako autonomní. Model Oregonator Model Oregonator, jehož název taktéž souvisí s místem vzniku (University of Oregon) byl vyvinut rovněž jako teoretická obdoba oscilační chemické reakce, byť je zřejmé, že jeho význam je stejně jako u modelu Brusselator obecnější. Modelové kinetické schéma je tvořeno pěti dílčími ději, v nichž figurují tři proměnné, další výchozí či produkované složky jsou obdobně jako v předchozím případě považovány za konstantní. I.

k1 A  Y  X

II.

k2 X  Y  P

III.

k3 B  X  2X  Z

IV.

k4 2X  Q

V.

k5 Z  fY

Odpovídající soustava diferenciálních rovnic má tedy podobu:   k AY  k XY  k BX  k X 2 X 1 2 3 4

(6.3)

  k AY  k XY  k Z Y 1 2 5

(6.4)

  k BX  k Z Z 3 5

(6.5)

Symboly složek systému v rovnicích mají opět význam přímo jejich koncentrací, příp. relativních množství. Nelinearita rovnic (6.3) a (6.4) vede k typickému oscilačnímu charakteru řešení. Obdobně jako u předchozího případu nelineárního dynamického systému vykazuje řešení modelu Oregonator protichůdnou oscilaci parametru X vzhledem k Y, proměnná Z kopíruje X.

45

Obr. 6.2 Výsledný oscilační průběh proměnných X, Y, Z obdržený numerickým řešením modelu Oregonator (6.3), (6.4), (6.5).

Lorentzův model konvekce v atmosféře Tento model, který je považován za průkopnickou práci v oblasti nelineárních dynamických systémů, popisuje konvekci v zemské atmosféře při slabě nadkritické hodnotě Rayleighova čísla, tj. problém na rozdíl od předchozích případů ryze fyzikální. Ve stručnosti lze uvést, že tento nelineární dynamický systém tvořený třemi obyčejnými diferenciálními rovnicemi, tj. zahrnující tři proměnné X, Y, Z byl vytvořen aplikací Galerkinovy metody (viz dodatek B tohoto textu) na výchozí nelineární rovnice pro tepelnou konvekci v tekutině (odvozené s obecných N-S rovnic) v prostoru tří vlnových módů charakterizujících idealizovaný obraz konvektivních (tj. vertikálním teplotním gradientem podmíněných) dějů v tekutině. Proměnná X odpovídá rychlosti konvektivního proudu, Y je ekvivalentní rozdílu teplot mezi vzestupným a sestupným proudem v zezdola zahřívané tekutině a Z reflektuje odchylku vertikálního teplotního profilu od profilu lineárního.   X  Y X

(6.6)

  XZ  rX  Y Y

(6.7)

  XY  bZ Z

(6.8)

Pro úplnost lze doplnit, že σ je Prandtlovo číslo, r je poměr Rayleighova čísla k jeho kritické hodnotě Ra/Rakrit (rozpracování konkrétních fyzikálních, resp. meteorologických aspektů atmosférické konvekce je nad rámec tohoto textu). Při této kritické hodnotě se tekutina s jistým vertikálním teplotním profilem stává nestabilní a vznikají konvektivní pohyby. Výsledky numerického řešení Lorentzova modelu jsou standardně zobrazovány jako tzv. fázový portrét, tj. (fázový) tok řešení v prostoru tří parametrů X, Y, Z. Časová osa tedy není standardně zobrazena jako u předchozích případů, je zde vlastně souběžná s trajektorií fázového toku řešení. Fázovým portrétem Lorenzova modelu jako nelineárního dynamického systému je tzv. chaotický (bývá označován i jako podivný) atraktor – viz obr. 6.3. Zobrazení chaotického oscilačního průběhu jednotlivých parametrů v časové ose – jako u předchozích případů je samozřejmě taktéž jednoduše představitelné. Lorentzův model byl ve své době (1962) vytvořen v rámci rozboru problému numerické předpovědi počasí. Výsledky jednoznačně ukazují základní vlastnost chování nelineárních dynamických systémů (popisujících nejen konvektivní procesy v atmosféře či tekutinách obecně) a to extrémní citlivost toku řešení na počátečních podmínkách.

46

Obr. 6.3 Fázový portrét numerického řešení Lorentzova modelu (6.6), (6.7), (6.8) v prostoru parametrů X, Y, Z – chaotický (či podivný) atraktor

Jak je patrno, tok řešení Lorentzova modelu je atrahován do určité podoby, resp. oscilace jednotlivých parametrů vykazují jistý řád a pravidelnost. Na druhé straně již nepatrná změna počátečních podmínek vede již po nevelkém počtu oscilací (či cyklů) ke zcela rozdílné poloze bodu v trajektorii toku řešení ve fázovém prostoru X, Y, Z. Tato principiální vlastnost nelineárních dynamických systémů označovaná jako deterministický chaos je podstatou např. praktické nemožnosti dlouhodobě předpovídat počasí. Pro úplnost lze dodat, že možností pro formulaci nelineárních dynamických systémů je samozřejmě bezpočet, zde prezentované tři případy jsou snad nejvíce v povědomí vedle modelu Lotka-Volterra popisujícího biologický (či ekologický) systém dravec – kořist. Tento nejjednodušší nelineární dynamický systém vykazující oscilační průběh parametrů byl formulován již v roce 1925 (lze jej najít např. v Begon, Harper: Ecology).

47

7. Příklad řešení kinetiky iontové výměny při průtoku dráhou s aktivním povrchem Fyzikální úlohy řešené v prostoru kartézských souřadnic (tj. R2, R3) a v čase jsou v obecnosti popisovány parciálními diferenciálními rovnicemi. Problém kinetiky iontové výměny při průtoku dráhou s aktivním povrchem tento přístup s ohledem na jednu prostorovou dimenzi a uvažovanou konstantní rychlost průtoku nevyžaduje a je prezentován jako příklad specifických úloh tohoto typu. dx

0

v

X

Obr. 7.1 Uvažovaná představa průtoku profilových elementů drahou x Formulace úlohy je prezentována na obr. 7.1. Roztok s obsahem vyměňovaných iontů, tj. vodíkových kationtů H+ a kationtů kovů Me+, protéká rychlostí v v 1D dráze o délce X s aktivním povrchem (naznačeném tečkami) s navázanými H+, Me+ označovanými v této souvislosti jako HS a MeS. Mezi protékajícím roztokem a statickým aktivním povrchem je tedy uvažována iontová výměna: k1

MeS  H 

Me  H S

k2

Fixace iontů na aktivním povrchu, která reálně souvisí se sorpčními faktory, je v rámci uvažované formulace problému vyjádřena zjednodušeně jako obecný proces druhého řádu závislý na koncentraci těchto iontů v roztoku a na relativním počtu (sorpčních) míst pro výměnu na aktivním povrchu, tj. míst obsazených vyměňovaným iontem - c S , c S . Tj.lze Me

H

psát:

dc

Me 

dt

 k1c  c H

Me S

 k2c

Me 

c

(7.1)

HS

Modelový výpočet výsledného průběhu c

Me 

 f ( x, t ) a z něj odvoditelných funkcí pro další

složky naznačené iontové výměny lze v tomto případě založit na představě profilového elementu šířky dx protékajícího v dráze x a interagujícího s aktivním povrchem – viz obr. 7.1. S ohledem na konstantní rychlost průtoku v může být časový krok „odměřen“ odpovídajícím posunem profilového elementu v dráze dx = vdt. Změna koncentrace kovu v roztoku vyplňujícím profilový element pak může být vyjádřena v závislosti na posunu v dráze:

dc

Me 

dx





1 k1c  c S  k 2 c  c S H Me Me H v



(7.2)

Pro řešení této výchozí rovnice je samozřejmě nutno eliminovat neznámé na pravé straně. Koncentraci vodíkových kationtů v profilovém elementu lze s ohledem na uvažovanou

48

iontovou výměnu (probíhající při průchodu tohoto elementu dráhou s aktivním povrchem) vyjádřit jednoduchou bilancí:

c

H



( x, t )  cH (0, t )   . cMe  (0, t )  c

Me 

( x, t )



(7.3)

Faktor  reflektuje možný relativní úbytek H+ v roztoku související s případným vytěsňováním dalších kovů z aktivního povrchu do roztoku paralelně s kovem Me+, jehož koncentrace je počítána. Míra tohoto paralelního vymývání kovů je v rámci použitého zjednodušení shodná s chováním Me+. Z hlediska statického aktivního povrchu reprezentuje daný posun profilového elementu změnu stavu tohoto povrchu (determinovaného c S a c S ) v prostorovém kroku dx = vdt o Me

H

časový krok dt. V celé dráze toku tedy platí: dc

Me S

dt

 v

dc

Me 

(7.4)

dx

Sumární změna „koncentrace“ iontu kovu v aktivním povrchu c

Me S

v celé dráze x  <0,X> v

časovém kroku tedy odpovídá nárůstu či poklesu (v důsledku zpětné fixace) koncentrace c  v profilovém elementu, který protekl touto dráhou. Me

X

d c

Me S

dx c

0

dt

Me 

(0, t )  c

Me 

( X , t)

(7.5)

Následující profilový element reprezentuje pro aktivní povrch další časový krok ap. – viz obr. 7.1. Pak tedy hodnota c S fixovaná v bodě x v dráze toku odpovídá celkovému času t = τ, Me

který je „odměřován“ integrálem proteklých (a ionty vyměňujících) profilových elementů. 

c

Me

0 S ( x, )  c Me S  v 

dc

Me 

( x, t )

dx

0

dt

(7.6)

0 0 Pro počáteční cMe S a c H S jsou uvažovány konstantní hodnoty v celé dráze x. Obdobně může

být vyjádřena c

HS

( x, ) : 

c

H

0 S ( x, )  c H S  v  0

dc

Me 

( x, t )

dx

dt

(7.7)

Dosazením (7.3), (7.6) a (7.7) vzniká výchozí diferenciální rovnice, jejíž řešení, tj. průběh c  (prostřednictvím profilového elementu) v dráze x odpovídající jednomu časovému Me

kroku, je hledáno:

dc

Me 

dx

   dc         0 Me      dt   k1  c H  (0, t )    c Me  (0, t )  c Me    c Me S  v      dx   0 1         dc v    Me   0  k 2 c Me   c H S  v  dx dt   0    

(7.8)

Po úpravě:

49

dc

Me 

( x,  1) dx



  dc Me  ( x, t )   0 1  k1  c H  (0,  1)  c Me  (0,  1)  c Me  v dt   S    v d x 0  

 dc   ( x, t ) 1 0 0 c   k 2  k1  Me dt   .k1c Me S  k 2 cH S Me d x v  0



  



(7.9)

Integrály přes časový interval <0,τ> (odpovídající (7.6) a (7.7)) se vztahují k aktuálnímu stavu aktivního povrchu (v bodě x), jenž bude interagovat s řešeným profilovým elementem, tj. elementem, reprezentovaným derivací c  na levé straně (7.9). Tento aktuální stav povrchu Me

vznikl interakcí s profilovými elementy „proteklými“ před elementem řešeným. Tj. časové integrály načítají změny c S v bodě x reprezentované derivacemi c  do stavu v čase τ, Me

zatímco derivace c

Me 

Me

( x,  1) na levé straně (7.9) již vyjadřuje stav v dalším časovém kroku.

Integrály přes časový interval mohou tedy být v (7.9) v obecnosti brány jako nezávislé funkce x a tato rovnice tedy odpovídá obyčejné diferenciální rovnici typu: y  g ( x)  yf ( x)

(7.10)

Obecné řešení rovnice (7.10) je obdrženo jednoduchým postupem za použití substituce y = uv (viz postup v kap. 5.5 - Následné procesy):

y ( x)  e  F ( x )  g ( x)e F ( x ) dx ,

F ( x)   f ( x)dx

(7.11)

Analytické řešení integrálů v (7.11) pro konkrétní funkce f(x), g(x) samozřejmě vede k následnému vyčíslení příslušných integrálních konstant. V rámci obecného tvaru lze tyto integrály přepsat na určité, přičemž interval integrace (na ose x) <0, xl > bude odpovídat bodu v argumentu počítané y(xl). Pro integrální křivku jdoucí bodem [x0, y0] pak vychází: y ( xl )  e

 F ( xl )

xl xl    y 0  g ( x)e F ( x ) dx  , F ( xl )  f ( x)dx 0 0    

(7.12)

Na základě (7.12) lze tedy hledané řešení (7.9) pro průběh c

Me 

v tekoucím profilovém

elementu v dráze x psát ve tvaru:     c Me  (0,  1)   1  c  ( xl ,  1)  e  z ( xl , )  k1 c H  (0,  1)   .c Me  (0,  1)  Me v   xl    dc ( x, t )  z ( x , ) Me   c0  v 0 dx dt e dx   0  Me S    





(7.13)

Přičemž:  dc    ( x, t ) 1 0 0 z ( xl , )    k 2  k1  Me dt   .k1c Me S  k 2 c H S dx  dx v   0  0



xl



  .k 2  k1  c Me  ( xl , t )  c Me  (0, t )dt  0







(7.14)

1 0 0  .k1c Me S  k 2 c H S xl v

50

Primitivní funkci k integrandu v (7.13) obsahujícímu derivaci neznámé funkce v součinu s ez nelze na rozdíl od integrálu v (7.14) (kde figuruje pouze tato derivace) vyjádřit, tj. lze provést pouze úpravu:  dc ( x, t )  z ( x , )  0 Me  c S  v 0  Me 0 dx dt e dx   

xl

(7.15)

xl  dc ( x, t )  1 1  z ( xl , ) 0 z ( x , ) Me    c Me e   v e dt dx S       z ( x ,  ) z ( 0 ,  ) d x l   0 0

Ve tvaru (7.15) se odráží (7.14), resp. nulová hodnota z v bodě počátku dráhy x = 0. Řešení problému (7.13) v dráze x tedy nabývá tvaru:

c

Me 

( xl ,  1)  e  z ( xl , ) c Me  (0,  1) 

1 0  1 1      e  z ( xl , )  c Me S  z (0, )    z ( xl , ) v  k1 c H  (0,  1)   .c Me  (0,  1)   xl  dc ( x, t ) z ( x , ) Me  e  z ( xl , ) 0 0 dx dt e dx  





(7.16)

Kde samozřejmě: 

z ( x, )   k 2  k1 

dc

Me 

( x, t )

dx

0

dt 



1 0 0  .k1cMe S  k 2 cH S v



(7.17)

Rovnice (7.16) (vč. (7.14) a (7.17)) se samozřejmě vztahuje pouze k časovému kroku τ + 1, tj. vyjadřuje průběh c  ( x,  1) v profilovém elementu v dráze x, která je ve stavu τ po Me

„protečení“ předchozího elementu reprezentovaného c

Me 

( x, ) , resp. derivací c

Tato situace tedy odpovídá rekurenci. Výchozí průběh c

Me 

Me 

( x, ) .

( x,1) (+ derivace) pro první

profilový element interagující s povrchem dráhy v počátečním homogenním stavu je snadno vyjádřitelný: 1  K1 x  K2 K  c  ( x,1)   e v  cMe  (0,1)  2  Me K1 K1  

dc

Me 

( x,1)

dx

Kde:

1



1  K1x  e v K 2  K1cMe  (0,1) v

0 0 K1   .k1cMe S  k 2 cH S ,

(7.18)



(7.19)



0 K 2  cMe S k1 c H  (0,1)   .c Me  (0,1)



Relativně „použitelný“ analytický tvar s vyčíslenou primitivní funkcí integrandu v (7.15) lze obdržet ještě pro následující časový krok:

51

c

Me 

( x,2)  e  z ( x ,1) c Me  (0,2) 

1 0  1  e  z ( x ,1)     c   Me S    z ( x,1) z (0,1)  v  (7.20)  1   k1 c H  (0,2)   .c Me  (0,2)   K1 x   e v e  z ( x ,1)  K  K c  (0,1)    1 Me 1 1   2    z ( x,1)  K 1 z (0,1)  K 1   v v   









Kde: 1  K1 x   K2  1 K  z ( x,1)   .k 2  k1   e v  cMe  (0,1)  2   cMe  (0,1)  K1 x K1    K1  v 1





1  K1x 1 z ( x,1)   k 2  k1  e v K 2  K1cMe  (0,1)  K1 v v

Je zřejmé, že analytický tvar pro n-tý časový krok je teoreticky vyjádřitelný v podobě binárního stromu. Je-li zavedena zobecněná struktura problému v časové ose, kdy průběh c  v profilovém elementu v n-tém časovém kroku může být jednoduše označován jako cn a Me

odpovídající funkce (tj. řešení) (7.13), resp. (7.16) jako fn, lze ve smyslu dané rekurence psát: n 1

c n  f n ( ci )

(7.21)

i 1

ci figurují na pravé straně (7.13) a (7.16) samozřejmě také jako suma jejich derivací v dráze x nicméně toto je v rámci uvažovaného zobecnění detailem neprojevujícím se v jeho struktuře. Např. řešení pro c3 je jednoduše dáno:

c3  f 3 ( f 2 (c1 )  c1 ) Dále pak:

c4  f 4 ( f 3 ( f 2 (c1 )  c1 )  f 2 (c1 )  c1 ) c5  f 5 ( f 4 ( f 3 ( f 2 (c1 )  c1 )  f 2 (c1 )  c1 )  f 3 ( f 2 (c1 )  c1 )  f 2 (c1 )  c1 )

(7.22)

apod. Tj. struktura řešení v časové ose odpovídá binárnímu stromu typu (7.22). Hledání explicitního analytického tvaru řešení pro n-tý časový krok zjevně není praktické. Celkový počet výskytů funkcí odpovídajících jednotlivým časovým krokům v daném binárním stromě pro vyjádření cn je vyčíslen jako (samozřejmě c1 = f(c0)):

m fi  2 ni 1 , i = 1, n-1

(7.23)

Při praktických aplikacích tohoto konceptu je rekurentní načítání v časové ose řešeno v číselné podobě na počítači (jako rekurentní posloupnost s počáteční hodnotou c1) . Nejde zde o numerickou aproximaci - krok dx je s ohledem na existenci analytického řešení (7.16) v dráze x infinitezimálním, tj. dt dle (7.4) taktéž. Numerické přiblížení generující nepřesnost je nutno použít pouze pro výpočet integrálu (7.15), jehož explicitně vyjádřený integrand má podobu binárního stromu (7.22). Lze dodat, že rozpracovaný koncept poskytující semianalytické řešení kinetiky iontové výměny při průtoku 1D dráhou s aktivním povrchem v závislosti na vtokových koncentracích

52

vyměňovaných iontů je v obecnosti aplikovatelný na řadu reálných dějů. V laboratorních podmínkách jde samozřejmě o modelování interakcí v kolonách různého typu, v přírodním prostředí je daný model velmi dobře použitelný např. pro simulace acidifikačních či jiných epizod ve vodních tocích.

53

Dodatek A – Základní numerické metody řešení problémů s počátečními podmínkami, tj. obyčejných diferenciálních rovnic a jejich soustav Jde o systémy se soustředěnými parametry, které jsou popsány soustavami obyčejných diferenciálních rovnic pro stavové proměnné v závislosti na čase. Eulerova jednokroková metoda explicitní Numerické řešení obecné úlohy s počáteční podmínkou vyjádřené obyčejnou diferenciální rovnicí y (t )  f (t , y(t ))

(A.1)

y(t0 )  y0 je založeno na aproximaci derivace hledané funkce na levé straně (A.1) její diferencí ve zvoleném časovém kroku h.

y(t  h)  y (t ) h

(A.2)

y(t  h)  y(t ) h f (t , y(t ))

(A.3)

y(t )  Tj.:

Pro výpočet nové hodnoty y(t+h) se tedy nabízí známá hodnota y(t), resp. f(t,y(t)) z předchozího kroku, což je principem Eulerovy explicitní metody. y n1  y n h f (t n , y n )

(A.4)

Numerická aproximace se tedy blíží přesnému řešení při h→0. Lze dodat, že (A.4) dává rovnítko mezi diferenci (yn+1–yn)/h a funkci pravé strany v bodě tn, yn. Charakter chyby související s touto skutečností znázorňuje obr. A.1. Hodnota funkce pravé strany v bodě tn+1/2, yn+1/2, která by zřejmě nejlépe odpovídala dané diferenci, samozřejmě není k dispozici.

Obr. A.1 Chyba vznikající při numerickém řešení Eulerovou jednokrokovou explicitní metodou. Diference v prvním kroku je pro daný tvar funkce přesného řešení nadhodnocena, což se projeví i v dalších krocích – posloupnost numerického řešení vůbec „netuší“, že přesné řešení je níže. Po čtvrtém kroku je pak situace opačná.

Eulerova jednokroková metoda implicitní (zpětná) V rámci Eulerovy jednokrokové metody lze vycházet i z funkce pravé strany v bodě tn+1, yn+1.

54

y(t  h)  y(t ) h f (t  h, y(t  h))

(A.5)

y n1  y n h f (t n1 , y n1 )

(A.6)

Neznámá yn+1 je nyní obsažena i na pravé straně (A.6). Hodnota funkce pravé strany obdržená po vyjádření neznámé yn+1 na stranu levou je samozřejmě k dispozici pro libovolný bod na časové ose (tj. i pro tn+1).

Obr. A.2 Chyba vznikající při numerickém řešení Eulerovou jednokrokovou implicitní metodou. Situace je logicky právě opačná oproti schématu explicitnímu. Stejná diferenční aproximace (Eulerova explicitní, implicitní) může být samozřejmě užívána i v případě soustav obyčejných diferenciálních rovnic. V aktuálně řešené rovnici soustavy lze pro funkci pravé strany, která zpravidla zahrnuje ostatní proměnné, volit u některých z těchto proměnných mezi stavem z předchozího kroku n a stavem n+1, který je již k dispozici, jde-li parametry počítané z rovnic předcházejících této aktuálně řešené. Runge – Kutta metody Jde opět o jednokrokový přístup, který však poskytuje významně přesnější numerické řešení oproti elementární metodě Eulerově. Tento koncept, jehož podstatou je propracovaný odhad hodnoty funkce pravé strany v bodě korektně odpovídajícím aproximující diferenci, navrhli němečtí matematikové Carle David Tolmé Runge, Martin Wilhelm Kutta. Výchozí numerické schéma řešené obyčejné diferenciální rovnice y   f (t , y) je ve tvaru:

y n1  y n  h(t n , y n , h)

(A.7)

Přičemž Φ(tn,yn,h) je hledaný odhad funkce pravé strany. Lokální chybu vztahující se k (A.7) jako k obecné jednokrokové metodě lze vyjádřit: L( y(t ); h)  y(t  h)  y(t )  h(t , y(t ), h)

(A.8)

Relativní přírůstek přesného řešení může být v diskrétní podobě dán diferencí levé strany při nenulovém kroku h.

(t , y, h) 

y (t  h)  y (t ) , h≠0 h

(A.9)

V obecnosti (h = 0) samozřejmě odpovídá funkci pravé strany: (t , y, h)  f (t , y) ,

h=0

(A.10)

Pro lokální chybu tedy platí:

L( y(t ); h)  h(t , y(t ), h)  (t , y(t ), h)

(A.11) 55

S ohledem na (A.10) odpovídá krokový přírůstek přesného řešení hΔ(t,y,h) Taylorovu rozvoji funkce pravé strany f(t,y(t)) přes jednotlivé mocniny kroku h.

1 1 h(t , y, h)  y(t  h)  y(t )  hf (t , y)  h 2 f   ...  h p f  p 1 (t , y)  (h p 1 ) (A.12) 2 p! (h p 1 ) je reziduální chyba v řádu příslušné mocniny h. Pro relativní přírůstek přesného řešení pak platí: y(t  h)  y(t ) 1 1 (A.13) (t , y, h)   f (t , y)  hf   ...  h p 1 f  p 1 (t , y )  (h p ) h 2 p! Nejlepším odhadem funkce pravé strany je tedy rozvoj: 1 1 (A.14) (t , y, h)  f (t , y )  hf   ...  h p 1 f  p 1 (t , y) 2 p! Vyčíslení příslušných derivací f(t,y) v (A.14) (vyjádřených indexem v závorce) není v rámci jednokrokové metody k dispozici. Optimální hodnota funkce pravé strany pro vyčíslení relativního přírůstku přesného řešení je samozřejmě f(t,y) v bodě topt, y(topt) ležícím v intervalu tn ≤ topt ≤ tn+1. Φ tedy může být hledána jako lineární kombinace funkčních hodnot f z tohoto intervalu. Jsou-li tyto funkční hodnoty vyjádřeny jako členy ki, lze psát: s

(t , y, h)   wi k i

(A.15)

i 1

V rámci výchozího numerického schématu jednokrokové metody je samozřejmě k dispozici pouze hodnota yn (z aktuálně vyčísleného kroku n), tj. funkce pravé strany f je vyčíslena v bodě tn, yn. Podstatou metod Runge Kutta je rekurentní konstrukce aproximujících členů funkčních hodnot pravé strany ki. Primární výchozí hodnotou je známá f(tn, yn) (dále psána jako f(t,y)). Tj.:

k1  f (t , y)

(A.16)

Další aproximující členy jsou konstruovány dle principu: i 1

k i  f (t   i h, y  h  ij k j ) ,

i = 2,...,s

(A.17)

j 1

Toto schéma přímo plyne s definice diferenciální rovnice, neboť obecný krokový posun v ose y je dán jako y   hy , tj. y   hf (t , y) , f (t , y)  k1 atd. Resp. v ki je vnořeno i – 1 (v rámci sumy v argumentu také i – 2, ...,1) cyklů typu prediktor – korektor. Lze dodat, že jednou ze zavedených metod numerického řešení obyčejných diferenciálních rovnic je právě postup prediktor – korektor, kdy je výchozí hodnota yn+1, která je vypočtena standardně Eulerovou explicitní metodou, následně dosazena do funkce pravé strany f(t,y) a je vyčíslena nová korigovaná hodnota yn+1. Tento cyklický postup může být opakován např. dokud není dosaženo zvoleného konvergenčního kritéria. Aproximující funkční hodnoty pravé strany ki v navrženém tvaru (A.17) jsou v rámci počítaného kroku s ohledem na rekurentní charakter (k1 = f(t,y) je k dispozici) vyčíslitelné za předpokladu nalezení hodnot faktorů αi, βij. Dále je nutno vyvodit hodnoty váhových faktorů wi, aby mohla být vyčíslena hledaná Φ(t,y,h). Je zřejmé, že kritériem správnosti hodnot těchto faktorů je soulad Φ vyjádřené dle (A.15) s Taylorovým rozvojem Φ (A.14) jako jejím nejlepším odhadem. Resp. wi a jejich součiny s αi, βij by měly být pro členy reprezentující

56

příslušné mocniny hp-1 v (A.14) ekvivalentní faktorům 1/p!. V jednotlivých ki jsou h, αi, βij součástí argumentu a do požadované formy koeficientů mohou být převedeny (Taylorovým) rozvojem funkcí ki pro daný posun t + αih, y + βhki-1. Pro k2 tedy lze s ohledem na tvar totální první a druhé derivace řešené rovnice y   f t (t , y)  f y (t , y) y   f t (t , y)  f y (t , y) f (t , y)

(A.18)

 y    f t (t , y )  f y (t , y ) y   f tt (t , y )  f ty (t , y ) y   f yt (t , y ) y   f yy (t , y ) y  2 

(A.19)

f tt (t , y )  2 f ty (t , y ) f (t , y )  f yy (t , y ) f (t , y ) 2 psát (derivace y’ vzniklé z argumentu f(t,y) jsou konstantou):

k 2  f (t   2 h, y   21hk1 )  f (t   2 h, y   21hf (t , y )) 





f (t , y )  h  2 f t (t , y )   21 f y (t , y ) f (t , y ) 

(A.20)





h2 2  2 f tt (t , y )  2 2  21 f ty (t , y ) f (t , y )   212 f yy (t , y ) f (t , y ) 2  (h 3 ) 2 Příslušné derivace jsou vyjádřeny dolními indexy. Nejjednodušší aplikace Runge - Kutta přístupu je vyčíslení wi, αi, βij pro druhý řád, tj. lineární kombinaci:

(t , y, h)  w1k1  w2 k 2

(A.21)

V tomto případě je dostačujícím rozvoj k2 v řádu h1, tj:

(t , y, h)  w1 f (t , y)  w2 f (t   2 h, y   21hf (t , y)) 

w1  w2  f (t , y)  w2 h 2 f t (t , y)   21 f y (t , y) f (t , y)

(A.22)

Porovnáním s odpovídajícími členy mocnin h (tj. h0, h1) Taylorova rozvoje (A.14) vychází:

w1  w2  1 ,

w2 

1 , 2

w2  

1 2

(A.23)

Soustava tří rovnic (A.23) obsahující čtyři neznámé má samozřejmě nekonečně mnoho řešení. Pro každé z těchto řešení existuje libovolné p ≠ 0, resp. řešení je dáno nastavením tohoto parametru. Lze psát:

w1  1  p ,

w2  p ,



1 , 2p



1 2p

(A.24)

Parametr p je nejčastěji volen p1 = 1 a p2 = ½, z čehož plyne:

1  f (t n  12 h, y n  12 hf (t n , y n ))

(A.25)

2 

(A.26)

1 2

f (t n , y n )  12 f (t n  h, y n  hf (t n , y n ))

Resp.:

yn1  yn  hf (tn  12 h, yn  12 hf (tn , yn ))

(A.27)

y n1  y n  12 h f (t n , y n )  f (t n  h, y n  hf (t n , y n ))

(A.28)

57

Vyčíslení koeficientů wi, αi, βij pro vyšší řády Runge – Kutta metod je již obsáhlejší. Pro třetí řád lze psát: k 3  f (t   3 h, y  h 31k1   32 k 2 )  f (t , y )  h  3 f t   31k1   32 k 2  f y 





(A.29)





h2 2  3 f tt  2 3  31k1   32 k 2  f ty   31k1   32 k 2 2 f yy  (h 3 ) 2

Derivované členy ft, fy, ftt, fty, fyy jsou v (A.29) brány jako nezávislé na konkrétním bodě v intervalu , resp. . Tento bod, ke kterému se derivace f vztahuje, figuruje v součinu s fy (tj. v tvaru derivace složené funkce f), kde y’ = (β31k1 + β32k2), resp. (y’)2 = (β31k1 + β32k2) 2 v součinu s fyy ap. Dále lze psát:

k 3  f (t , y ) 





h  3 f t   31   32  f y f (t , y )   32 hf y  2 f t   21 f y f (t , y )  



(A.30)



h2 2  3 f tt  2 3  31   32  f ty f (t , y )   31   32 2 f yy f (t , y ) 2  (h 3 ) 2 Tj. v členu rozvoje k3 v h1 je uvažován rozvoj k2 taktéž v řádu h1. V dalším členu rozvoje k3 (v h2) je k2 „rozvinuto“ pouze do h0, tj. zde k2 = k1 = f(t,y) – pak tedy platí:

 31k1   32k 2 2 f yy   31   32 2 f yy f (t, y) 2 Tato postupná redukce rozvoje k2 figurujícího v k3 je logická v souvislosti s tím, že pro Runge-Kutta metody třetího řádu jsou dostačující členy rozvoje (A.14) do h2. Pro k3 tedy po zpřehledňující úpravě (A.30) platí:

k 3  f (t , y ) 





h  3 f t   31   32  f y f (t , y ) 

(A.31)

  32 2 f t f y   32  21 f y2 f (t , y )     ( h 3 ) h 1 2 2 1 2  2  3 f tt   3  31   32  f ty f (t , y )  2  31   32  f yy f (t , y )  2

Aproximující funkci pravé strany Φ pro třetí řád Runge-Kutta metod lze vyjádřit na základě (A.20) a (A.31):

 (t , y, h)  w1 k1  w2 k 2  w3 k 3  w1 f (t , y ) 





 f (t , y )  h  2 f t   21 f y f (t , y )   w2  2 1 2  1 2 2  h 2  2 f tt   2  21 f ty f (t , y )  2  21 f yy f (t , y )   f (t , y )  h  3 f t   31   32  f y f (t , y )     2   w3  2   32 2 f t f y   32  21 f y f (t , y )   h  1 2 2 1 2    2  3 f tt   3  31   32  f ty f (t , y )  2  31   32  f yy f (t , y )  









(A.32)

58

Přičemž užitečným je následující rozpis:  (t , y, h)  w1  w2  w3  f (t , y ) 









hw2  2 f t   21 f y f (t , y )  w3  3 f t   31   32  f y f (t , y ) 





2 w2 12  22 f tt   2  21 f ty f (t , y )  12  21  f yy f (t , y ) 2      h 2    32 2 f t f y   32  21 f y2 f (t , y )   w3  1 2 2 1 2       f      f f ( t , y )     f f ( t , y )    3 31 32 ty 31 32 yy 2   2 3 tt

 w1  w2  w3  f (t , y ) 

w2 2  w3 3  f t   h  w2  21  w3  31   32  f y f (t , y )  12 w2 22  w3 32  f tt      w2 2  21  w3 3  31   32  f ty f (t , y )    h2 1  2 2 2  2 w2  21  w3  31   32  f yy f (t , y )     2 w3  32 2 f t f y  w3  32  21 f y f (t , y ) 





(A.33)

Porovnáním (A.33) s (A.14) tedy vychází:

w1  w2  w3  1

(A.34)

1 2

(A.35)

w2 2  w3 3 

w2  21  w3  31   32  

w2 22  w3 32 

1 2

(A.36)

1 3

(A.37)

w2 2  21  w3 3  31   32   2 w2  21  w3  31   32   2

1 6

1 3

(A.38) (A.39)

w3 2  32 

1 6

(A.40)

w3  21 32 

1 6

(A.41)

Z (A.40) a (A.41)) plyne rovnost:

 2   21

(A.42)

Tj. pro (A.35) a (A.36) vyplývá identita a platí:

 3   31   32

(A.43)

59

Identické jsou tedy i (A.37) a (A.39) a dále (A.38) je s (A.39) či s (A.37) v rozporu a nelze ji pro výpočet neznámých parametrů uplatnit. Tato rozpornost (A.38) není bez souvislosti – tato rovnice (a samozřejmě ta také (A.63), (A.39)) ve výsledném tvaru nefiguruje, je-li s relacemi (A.42), (A,43) počítáno již při odvození k2, k3 a Φ. Toto korektnější vyjádření Φ má podobu:  (t , y, h)  w1  w2  w3  f (t , y ) 









hw2  2 f t   21 f y f (t , y )  w3  3 f t   31   32  f y f (t , y ) 





2 w2 12  22 f tt   2  21 f ty f (t , y )  12  21 f yy f (t , y ) 2   h 2   12  32 f tt   3  31   32  f ty f (t , y )  12  31   32 2 f yy f (t , y ) 2 w3  2    32 2 f t f y   32  21 f y f (t , y )



      

(A.44)



 12 w2 22  w3 32 G    w1  w2  w3  f (t , y )  hw2 2  w3 3 F   h    w3  32 2 F  2

Přičemž: F  f t  f y f (t , y) ,

G  f tt  2 f ty f (t , y)  f yy f (t , y) 2

K dispozici tedy jsou čtyři rovnice (A.34), (A.35), (A.37), (A.40) při šesti neznámých (w1-3, α1, α2, β32) - řešení je podobně jako u Runge- Kutta metod druhého řádu nejednoznačné. Lze dojít např. k:

w1 

1 1 1 2 1 2 , w2  , w3  ,  2   21  ,  3  ,  31   ,  32  1 3 4 2 4 3 3

(A.45)

w1 

1 2 1 1 , w2  , w3  ,  2   21  ,  3  1 ,  31  1 ,  32  2 6 3 6 2

(A.46)

Nejčastěji však bývá aplikována sekvence:

1 3 1 2 2 , w2  0 , w3  ,  2   21  ,  3  ,  31  0 ,  32  4 4 3 3 3

(A.47)

y n1  y n  14 hf (t n , y n )  34 hf (t n  23 h, y n  23 hf (t n  13 h, y n  13 f (t n , y n )))

(A.48)

w1  Tj.:

Odvození parametrů Runge-Kutta metod čtvrtého řádu, který je pro numerické řešení diferenciálních rovnic nejefektivnější (vyšší řády již neznamenají významné zpřesnění), si může provést čtenář sám. Výsledné schéma je následující:

y n1  y n  16 hk1  2k 2  2k 3  k 4  k1  f (t n , y n )

k 2  f (t n  12 h, y n  12 hk1 )

(A.49)

k 3  f (t n  12 h, y n  12 hk 2 ) k 4  f (t n  h, y n  hk 3 )

60

Hodnoty odvozovaných parametrů v (A.49) reprezentují samozřejmě pouze jedno z nekonečného počtu možných řešení (i pro čtvrtý řád je řešení nejednoznačné), byť často užívané. V softwarových aplikacích pro numerické řešení diferenciálních rovnic (resp. jejich soustav) je často k dispozici také volba následujícího setu parametrů:

y n1  y n  18 hk1  3k 2  3k 3  k 4 

k1  f (t n , yn ) k 2  f (t n  13 h, y n  13 hk1 )

(A.47)

k 3  f (t n  23 h, y n  13 hk1  hk 2 ) k 4  f (t n  h, y n  hk1  hk 2  hk3 )

61

Dodatek B – Numerické řešení parciálních diferenciálních rovnic – aproximace Galerkinovou metodou Obecná východiska Obecná parciální diferenciální rovnice může být vyjádřena předpisem: Au~  f  0

(B.1)

Kde u~ je přesné řešení, tj. funkce např. v prostoru R3, f je libovolná funkce prostorových souřadnic (např. v R3) a případně času (u parabolických diferenciálních rovnic), v řadě případů (hyperbolické a některé eliptické rovnice) je f ≡ 0. A je diferenciální operátor, tj. např.:

   ex  ey  ez   x y z

(B.2)

Nabla operátor ve fyzikálním R3 prostoru, ex, ey, ez jsou jednotkové vektory v tomto prostoru.

2 2 2    2   2 2 2 x y z

(B.3)

Laplaceův operátor v R3. Při popisu časově závislých dějů jsou aplikovány odpovídající diferenciální operátory:

2 2 2 1 2    □ x 2 y 2 z 2 c 2 t 2

(B.4)

D’Alembertův operátor figurující v rovnici vlny šířící se rychlostí světla c.

     ux  uy  uz t x y z

(B.5)

Nelineární operátor, který je součástí rovnic popisujících vektorové pole proudění a další děje v tekutinách. Disketizace parciálních diferenciálních rovnic (dále PDR) je v rámci konceptu tzv. vážených reziduí, kam patří Galerkinova metoda, založená na tom, že přesné řešení je aproximováno lineární kombinací konečného počtu zvolených funkcí. Tj. spektrum těchto tzv. bázových funkcí by mělo být konstruováno tak, aby se jejich specifická lineární kombinace mohla blížit přesnému řešení. n

u~  u   uii

(B.6)

i 1

Tyto funkce φi tvořící bázi jsou tedy lineárně nezávislé, tj. žádná z nich nesmí být vyjádřitelná jako lineární kombinace ostatních. Splnění této podmínky je zaručeno ortogonalitou bázových funkcí φi. Pro skalární součin vyjádřený integrálem přes oblast Ω v R3, ve které jsou zavedeny bázové funkce φi tedy platí:

 ,      d   i

j

i

j

ij

(B.7)



62

δij je Kroneckerovo delta, které nabývá hodnoty1 pro i = j a 0 pro i ≠ j. Funkce báze jsou tedy normované ((φi , φi ) = 1) a vyhovují požadavku: n

 i 1

i

1

(B.8)

Podmínka (B.8) je splněna pro libovolnou trojici (pokud je Ω uvažována v R3) hodnot lokálních proměnných v oblasti Ω dosazenou do bázových funkcí. Lze upřesnit, že bázové funkce jsou s ohledem na (B.7) a (B.8) ortonormální. Tato ortonormální báze je v této souvislosti definována v Hilbertově prostoru H, kde tvoří konečný podprostor bázových funkcí dimenze n, resp. aproximativní řešení PDR u je hledáno v tomto podprostoru bázových funkcí v H. Řešení (hledané dle (B.6)) se v obecnosti bude blížit řešení přesnému u~ v závislosti na početnosti zvolené báze, tj. dimenzi podprostoru bázových funkcí v H. Přibližné řešení u nevyhovuje PDR (B.1) přesně, tj. lze psát: Au  f  0 Au  f  r

(B.9)

Kde r je reziduum vyjadřující „míru“ nepřesnosti u. Řešení u ve tvaru (B.6) se tedy bude nejvíce přibližovat řešení přesnému u~ při nalezení souboru hodnot koeficientů ui reprezentujících ve smyslu (B.9) nejmenší možné reziduum (tj. koeficienty ui by měly v obecnosti dávat větší váhu bázovým funkcím přibližujícím se přesnému řešení u~ a naopak). Reziduum r je samozřejmě funkcí příslušných souřadnic na oblasti Ω a při jeho vyhodnocení, resp. minimalizaci lze vycházet pouze z principů uplatnitelných pro prostor funkcí H. Konkrétně lze reziduum „vážit“ prostřednictvím skalárního součinu (na oblasti Ω) se zvolenou tzv. váhovou funkcí w, kdy je tento skalární součin dvou funkcí (v H) položen rovný nule.

w, r   w, Au  f   0 Tj.:

 w Au  f  d  0 

(B.10)

Tato podmínka tedy znamená ortogonalitu rezidua a váhové funkce w - viz (B.7), nebo lze zjednodušeně řící, že reziduum a váhová funkce se na oblasti Ω „nepřekrývají“ (pro libovolnou (3D) souřadnici na oblasti Ω je vždy jedna z funkcí nulová a celý integrál jejich součinu je tedy nulový). Vedle této skutečnosti je nutno neopomenout, že rovnost (B.10) by byla splněna také v případě, že funkce reprezentující reziduum by byla na celé oblasti Ω nulová a u by tedy byla přesným řešením (w je v každém případě nenulová). V obecnosti pak platí rovnost:

 w Au  f  d   w Au~  f  d  0 



(B.11)

Relace (B.10) tedy vyhovuje v integrálním smyslu dané PDR (B.1), resp. funkce u vyjádřená z (B.10) je tzv. slabým (nebo zobecněným) řešením (B.1) na oblasti Ω. Jedna váhová funkce w samozřejmě může ve smyslu (B.10) pouze „vytěsnit“ reziduum mimo její nenulové hodnoty na oblasti Ω. Pro nalezení řešení u blízkého řešení přesnému u~ se zajisté nabízí možnost vážit reziduum vícekrát použitím více zvolených funkcí wi, resp. odpovídá-li počet váhových funkcí dimenzi n podprostoru bázových funkcí φi, je k dispozici uzavřený systém n integrálních rovnic (v případě časově nezávislého operátoru A) pro 63

výpočet n neznámých koeficientů ui a tedy pro nalezení přibližného řešení u. J-tý řádek dané matice je tedy ve tvaru:

  n  w u i i    j  A  i 1 

  n  f d    u i  w j A i d    w j fd  0  i 1    

(B.12)

Váhové funkce wi musí samozřejmě být vhodně zvoleny. Primárně se pro tuto volbu nabízí logická podmínka, aby výsledné reziduum bylo „vážením“ s n funkcemi wi „vytěsněno“ mimo podprostor bázových funkcí φi. Pak je zajištěno, že slabé řešení u ve tvaru (B.6) bude v rámci dimenze n tohoto podprostoru nejlepší možné – nejvíce se blížící řešení přesnému, resp. že daný soubor bázových funkcí už lepší variantu řešení nabídnout nemůže. Vyřešení soustavy (B.12) (tj. výpočet souboru ui ,i = 1, .., n) vede k reziduu, které musí být ortogonální se všemi wi. Je zřejmé, že jsou-li jako váhové funkce použity přímo funkce báze wi = φi, (i = 1, ..., n), je vyřešením (B.12) dosaženo žádoucího „vytěsnění“ rezidua mimo podprostor bázových funkcí. Tato varianta metody vážených reziduí (wi = φi) je dle svého autora nazývána Galerkinovou metodou. Pro korektnost je ještě třeba dodat, že pojem reziduum (který je zde komentován v jednotném čísle) fakticky reprezentuje funkci tvořenou ortogonálními složkami – rezidui (plurál ve smyslu názvu metody). Jednotlivé skalární součiny s váhovými (= bázovými φj) funkcemi pak nulují vždy příslušnou (při i=j) složku rezidua v podprostoru báze. Slabé řešení u ve tvaru (B.6) je obecně použitelnou aproximací na úrovni numerického řešení diskretizované PDR v rámci kroku korespondujícího s oblastí Ω (ztotožňované s konečným prvkem). Zjm. při řešení PDR s nelineárními operátory – např. (B.5) samozřejmě vyvstává možnost, že přesné řešení u~ bude mít v rámci měřítka Ω složitý fluktuační průběh, který nebude slabým řešením u (vyhovujícím v integrálním smyslu v měřítku Ω) v žádném případě vystižen. Resp. výsledné reziduum (funkce v H mimo podprostor bázových funkcí) bude na oblasti Ω nabývat nezanedbatelných funkčních hodnot. Stacionární úlohy Výsledná soustava algebraických rovnic vzniklá aplikací Galerkinovy metody na časově nezávislou PDR je ve tvaru:

u1  1 A1d  u n  1 A n d   1 fd        u1   n A1d  u n   n A n d    n fd 



(B.14)



Funkce Aφi vzniklé působením diferenciálního operátoru A na bázové funkce nejsou s váhovými (tj. bázovými) funkcemi φi ortogonální a jednotlivé integrály vycházejí v obecnosti nenulové (max. hodnota = 1). Stejně tak f, tj. integrály jejích součinů s φi tvoří nenulový vektor pravé strany (pokud nejde o typ PDR neobsahující f). Je třeba dodat, že výsledná soustava lineárních rovnic (B.14) umožňující po vyčíslení jednotlivých integrálů výpočet hodnot ui např. Gaussovou eliminací odpovídá časově nezávislému lineárnímu operátoru A (tj. jde o řešení stacionárního problému). Pro operátor nelineární (časově nezávislý) by v rámci jednotlivých členů výsledné soustavy figurovaly např. součiny uiuj – šlo by o soustavu nelineárních rovnic řešitelnou vhodnou iterační procedurou.

64

Časově závislé úlohy U časově závislých PDR vede aplikace Galerkinovy metody k soustavě rovnic v podobě:

u1 t u n t

   d  u   A d 

1 1

1



1

1

 n

n





   d  u   

 u n  1 A n d

1



n

  1 fd 





(B.15)

A1d  u n   n A n d    n fd 



Proměnnými v čase jsou tedy hledané koeficienty ui. Člen reprezentující časovou derivaci se samozřejmě v každém řádku objevuje pouze jednou (s ohledem na (B.7)). Integrály pro skalární součiny (φi, φi) v tomto členu jsou rovny jedné a nemusely by být vůbec vypisovány. (B.15) tedy reprezentuje soustavu obyčejných diferenciálních rovnic (dále ODR). Pro parabolické PDR jde o ODR standardního typu u (t )  f (t , u(t )) , pro hyperbolické PDR vznikají soustavy ODR typu u (t )  f (u(t )) . Aplikace Galerkinovy metody tedy převádí časově závislou PDR na soustavu obyčejných diferenciálních rovnic poskytující slabé řešení pro oblast Ω, tj.: n

u (t ,  , ,  )   u i (t ) i ( , ,  )

(B.16)

i 1

kde  , ,  jsou lokální prostorové souřadnice v Ω. Opět je nutno dodat, že nelineární operátor povede k soustavě nelineárních ODR. Soustava ve tvaru (B.15) také nebude vyhovující pro PDR obsahující druhou časovou derivaci, případně první a druhou zároveň. Rozbor těchto eventualit je již nad rámec tohoto textu. Možnosti numerického řešení soustav typu (B.15) jsou již nastíněny v dodatku A. Na tomto místě lze velmi stručně zmínit možné další alternativy. Jednou z nich je analytické řešení použitelné pro soustavu lineárních ODR s nulovým vektorem pravé strany. Tato soustava, která je výsledkem diskretizace lineární PDR s první časovou derivací a f ≡ 0 Galerkinovou metodou má následující tvar:

u1  u1  1 A1d  u n  1 A n d  0   t     u n  u1   n A1d  u n   n A n d  0   t

(B.17)

Řešení reprezentované n funkcemi časového průběhu jednotlivých ui lze předpokládat ve tvaru součinů časově zprůměrovaných u i s obecně vyhovující funkcí času: ui (t )  ui e t ,

(i = 1, ..., n)

(B.18)

Po dosazení (B.18) do (B.17) (a vykrácení) je obdržena homogenní soustava n lineárních algebraických rovnic pro neznámé složky vlastního vektoru u :



u1     1 A1d   u1   n A1d 



 u n  1 A n d     u n      n A n d







0  0

(B.19)

Rozpis (B.19) odpovídá maticovému zápisu:

A  λEu  0

(B.20)

65

Kde prvek Aij matice A je roven   i A j d a E je jednotková matice. Tvar (B.20) je v 

podobě problému vlastních čísel λ matice A, neboť platí: Au   u

(B.21)

Pro nenulové řešení homogenní soustavy (B.19) u musí být vlastní číslo λ voleno tak, aby:

det A  λE  0

(B.22)

Vyjádření determinantu (B.22) vede k charakteristické rovnici matice A pro neznámou λ. Kořeny charakteristické rovnice je obecně n hodnot λ, tj. n vlastních čísel λi matice A. Dílčím výsledkem je tedy n možných časových průběhů vektoru u , který je k dispozici v n podobách jako řešení (B.19) pro jednotlivá vlastní čísla λi, tedy jako n vlastních vektorů matice A. Jsouli obdržené vektory u i normalizovány na jednotkovou velikost, tj. pokud platí:

ui , ui   1

(B.23)

Lze vyjádřit obecné řešení soustavy (B.17): n

u(t )  ui e it  u e t

(B.24)

i 1

Slabé řešení lineární PDR generující soustavu (B.17) pro oblast Ω je tedy ve tvaru: n

u (t ,  , ,  )  e t  u i i ( , ,  )

(B.25)

i 1

Výsledné tvary (B.24) a (B.25) reprezentují v rámci obecného řešení (v čase) identickou evoluci pro časově závislé parametry ui. Tato situace dokumentuje linearitu, resp. výsledek tohoto typu by pro nelineární problém nepřipadal do úvahy. Další možnou alternativou pro řešení časově závislých PDR je přiřazení časové osy k prostorovým dimenzím oblasti Ω, tj. zavedení bázových funkcí zahrnujících závislost na čase. Časový průběh řešení u (samozřejmě pouze pro lokální časový interval oblasti Ω) je tedy počítán přímo v rámci diskretizace PDR Galerkinovou metodou a nikoli následně, jak je tomu u výše rozpracovaných přístupů. Podrobnější popis tohoto konceptu již přesahuje rámec textu.

Rozšíření na konečné prvky Při řešení praktických problémů v R3, zejména pak modelování transportních procesů zpravidla vyvstává potřeba rozšíření numerického řešení PDR na více oblastí. Tento požadavek je reflektován zavedením tzv. konečných prvků. Metoda konečných prvků konstruuje jednotlivé lokální oblasti Ω jako prvky, které vyplňují celou globální oblast řešení - jsou zavedeny prostřednictvím sítě. Galerkinovská diskretizace PDR (nebo soustavy PDR) v celé globální oblasti, jejíž tvar je dán fyzikálním zadáním, tedy generuje velkou řídkou soustavu s nenulovými členy v blízkosti diagonály. Základní způsoby tvorby sítě vedou ve 2D fyzikálním prostoru k trojúhelníkovým nebo čtyřúhelníkovým prvkům. Ve 3D jim odpovídají prvky čtyřboké a šestiboké. Problematika konstrukce konečných prvků a jim příslušejících bázových funkcí je značně široká, zájemce lze odkázat např. na práci R. Černého (1997). V rámci zaměření toho textu lze uvést pouze vybraný (prakticky nejjednodušší) příklad. 66

Čtyřúhelníkový prvek ve 2D

η Lη2 1

4

Lξ1

3

Lξ2 1

ξ

2

Obr. B.1 Čtyřúhelníkový prvek z vyznačeným průběhem lineárních bázových funkcí L1 – L4.

-1 Lη1 1

-1

Prvek je zaveden v lokálních prostorových souřadnicích ξ, η, které nemusí s ohledem na možný neortogonální tvar prvku korespondovat s globálními kartézskými souřadnicemi. Při diskretizaci PDR, které jsou formulovány v kartézských souřadnicích, je nutno provádět příslušné transformační operace mezi lokálními a globálními souřadnicemi. V 2D lokálních souřadnicích lze zavést lineární izoparametrickou bázi:

L1 ( ) 

1 1    , L2 ( )  1 1   , L1 ( )  1 1    , L2 ( )  1 1    2 2 2 2

(B.26)

Funkce (B.26) jsou tedy navrženy jako nejjednodušší polynomy prvního stupně s funkční hodnotou <0,1> v oblasti prvku - viz obr. B.1. Bázové funkce, které jsou „nabízeny“ pro možný průběh řešení PDR, musí být při 2D aplikaci samozřejmě funkcemi dvou prostorových proměnných. Z (B.26) je tedy odvozena bilineární -1 báze:

1  L1 ( ) L1 ( ) 

1 1   1    4

(B.27)

 2  L2 ( ) L1 ( ) 

1 1   1    4

(B.28)

 3  L2 ( ) L2 ( ) 

1 1   1    4

(B.29)

 4  L1 ( ) L2 ( ) 

1 1   1    4

(B.30)

Tj. bázové funkce (B.27 – 30) mají v příslušejících rohových bodech hodnotu jedna. Z obr. B.1 je patrno, že i takto jednoduchá báze čtyř funkcí „nabízí“ základní možnosti variability prostorového průběhu řešení u v oblasti prvku Ω. Pro praktické výpočty jsou k dispozici složitější typy bází. Např. báze bikvadratická, která odvozená ze tří základních typů polynomů druhého stupně, což pro čtvercový prvek ve 2D obecně vede k devíti bázovým funkcím, přičemž tento počet může být účelově redukován (osm funkcí – tzv. serendipovský prvek). Bázové funkce mohou být taktéž odvozovány z Lagrangeových či Hermitových polynomů, pak jde o Lagrangeovské a Hermitovské prvky.

67

Použitá literatura

MOORE, W. J.: Fyzikální chemie. 2. vyd. Praha: SNTL, 1981. 976 s. ISBN 04-604-81. ECKERT, E., HORÁK, J., JIRÁČEK, F.: Aplikovaná reakční kinetika. Praha: SNTL, 1986. 340 s. ISBN 05-088-86. FISCHER, O., a kol.: Fyzikální chemie. 1. vyd. Praha: SPN, 1984. 336 s. ISBN 14-517-84. VITÁSEK, E.: Základy teorie numerických metod pro řešení diferenciálních rovnic. 1. vyd. Praha: Academia, 1994. 409 s. ISBN 80-200-0281-2. ČERNÝ, R.: Řešení transportních jevů na počítači. Praha: ČVUT, 1997. 157 s. ISBN 80-0101580-7. KLÍČ, A., KUBÍČEK, M.: Matematika III. Diferenciální rovnice. Kvalitativní teorie a aplikace. Numerické metody. Praha: VŠCHT, 1992. 242 s. ISBN 80-7080-162-X. HORÁK, J., KRLÍN, L.: Deterministický chaos a matematické modely turbulence. 1. vyd. Praha: Academia, 1996. 445 s. ISBN 80-200-0416-5. HOLODNIOK, M., KLÍČ, A., KUBÍČEK, M.: Metody analýzy nelineárních dynamických modelů. Praha: Academia , 1986. 424 s. ISBN 21-010-86. NAVRÁTIL, T., VACH, M., NORTON, S.A., SKŘIVAN, P., HRUŠKA, J., MAGGINI, L., 2003: Chemical response of a small stream in a forested catchment (central Czech Republic) to a short-term in-stream acidification. Hydrology and Earth System Sciences, 3, (7), p. 411423. VACH, M., SKŘIVAN, P., 1996: Study of the sorption and dissolution processes in the system of fly ash - aqeous solution. Kinetic model of the leaching of arsenic., Conf. „Global Changes and Essential Elements Cycling in the Environment“, SCOPE, Geological Inst. of AS CR, Prague 17-18.9.1996, Conference proceedings, p. 16-17.

68