MASARYKOVA UNIVERZITA PŘÍRODOVĚDECKÁ FAKULTA Ústav fyziky kondenzovaných látek ______________________________________
Modelování a konstrukce uhlíkové elektrody z nanorozměrového materiálu Diplomová práce
Rostislav Turek
Vedoucí diplomové práce: Mgr. Dušan Hemzal, Ph.D.
1
Brno 2011
Bibliografický záznam
Autor:
Bc. Rostislav Turek
Název práce:
Modelování a konstrukce uhlíkové elektrody z nanorozměrového materiálu
Studijní program:
PřF N-FY Fyzika, magisterský studijní program
Studijní obor:
PřF BF Biofyzika
Vedoucí práce:
Mgr. Dušan Hemzal, Ph.D.
Rok obhajoby:
2011
Klíčová slova:
Uhlíkové elektrody, vodivost, metoda
konečných prvků
2
Bibliographic entry
Author:
Bc. Rostislav Turek
Title of thesis:
Modelling and construction of nano-carbon
electrode Degree programme: PřF N-FY Physics, Master's degree programme Field of study:
PřF BF Biophysics
Supervisor:
Mgr. Dušan Hemzal, Ph.D.
Year of defence:
2011
Keywords:
Carbon electrodes, conductivity, finite element
method
3
Rád bych poděkoval panu Mgr. Dušanu Hemzalovi, Ph.D. za vedení, trpělivost a vstřícný postoj při realizaci této práce. Chtěl bych poděkovat i své rodině za podporu.
Prohlašuji, že jsem svou diplomovou práci napsal samostatně a výhradně s použitím citovaných pramenů. Brno
....................... Rostislav Turek
4
Abstrakt : Uhlíkové elektrody patří v biofyzice k dlouhodobě používaným. Vlastnosti uhlíkových elektrod se značně liší v závislosti na mnoha faktorech, jak strukturních tak povrchových. Metoda konečných prvků představuje jeden z modelů schopných popsat elektrické vlastnosti oblasti, kterou není možné celou proměřit. Takovým modelem může býž rozložení potenciálu ve 2D oblasti. Část dat lze zajistit pomocí měření lineární čtyřsondy.
Abstract: Carbon electrodes are using in Biophysics long time. Carbon electrodes properties deviates due to many factors like structure and surface properties. Finite element method (FEM) can help to build a model which can describe electric properties in cases in which we can measure it directly. The model like mentioned can show distribution of electric potential in 2D. As a additional method how to get data for FEM, it can be used measuring by fourprobe method.
5
Obsah 1. Úvod.................................................................................................................................................8 2. Teoretická část..................................................................................................................................9 2.1 Vlastnosti elektrod [1]................................................................................................................9 2.1.1 Alotropy ............................................................................................................................9 2.1.1.1 Grafit a související sp2 hybridy..................................................................................9 2.1.1.2 Uhlíkové vlákna.........................................................................................................9 2.1.1.3 Glassy Carbon - skelný uhlík ..................................................................................10 2.1.1.4 Diamant....................................................................................................................10 2.1.1.5 Carbon Nanotubes....................................................................................................11 2.1.2 Elektronové vlastnosti uhlíkových elektrodových materiálů...........................................11 2.1.3 Spektroskopické vlastnosti...............................................................................................11 2.1.4 Elektrochemické vlastnosti..............................................................................................13 2.1.4.1 Povrchová stuktura uhlíkových elektrodových materiálů........................................13 Terminace ......................................................................................................................13 2.1.5 Příprava povrchu uhlíkových elektrod.............................................................................14 2.1.5.1 Leštění, broušení, čištění..........................................................................................14 2.1.5.2 Aktivace uhlíkových elektrod...................................................................................15 Povrchové čištění.............................................................................................................15 Aktivace laserem.............................................................................................................15 2.2 Maxwellovy rovnice pro makroskopické elektromagnetické pole [2] ...................................16 2.3 Metoda konečných prvků - Teorie Modelu [3].......................................................................18 2.3.1 Metoda konečných prvků (MKP)....................................................................................18 2.3.1.1 Generace sítě prvků..................................................................................................19 2.3.1.2 Aproximace potenciálu na jednotlivých prvcích z uzlových hodnot.......................20 Postup aproximace ve 2D prostoru .................................................................................20 Aproximace potencialu ve 3D úloze................................................................................21 2.3.1.3 Sestavení soustavy rovnic........................................................................................22 2.3.1.4 Vyřešení soustavy.....................................................................................................24 2.3.1.5 Výpočet dalších veličin – postprocesor....................................................................24 2.4 Měření Lineární čtyřsondou [4,5,6] ........................................................................................25 3. Praktická část..................................................................................................................................27 3.1 Ověření Ohmova zákona.........................................................................................................27 3.1.2 Neznámý vzorek uhlíkového materiálu...........................................................................27 3.1.2.1 Sonda v poloze paralelní ke straně AB.....................................................................28 3.1.2.2 Sonda v poloze paralelní ke straně BC.....................................................................29 3.1.2.3 Sonda v poloze paralelní k úhlopříčce BD...............................................................30 3.1.2.4 Sonda v poloze paralelní k úhlopříčce AC...............................................................31 3.1.3 HOPG ..............................................................................................................................33 3.2 Závislost velikosti měrné vodivosti na směru lineární sondy nad vzorkem. .....................35 3.3 Model stacionárního proudového pole [3]...............................................................................36 3.3.1 Laplaceova rovnice..........................................................................................................36 3.3.2 Upřesnění výpočtu MKP..................................................................................................36 3.3.3 Modely elektrod...............................................................................................................38 4. Závěr...............................................................................................................................................42 6
Příloha A.............................................................................................................................................44 Zadání generování sítě do programu Gmsh (příklad)....................................................................44 1) Pro model dvou bodových elektrod......................................................................................44 Skript programu matlab na výpočet matice K (příklad)................................................................45 1) MaticeK pro model 2 bodových elektrod.............................................................................45 2) Řešení soustavy ..................................................................................................................46
7
1. Úvod
Uhlíkové elektrody patří v biofyzice k dlouhodobě používaným. Mezi nejčastější důvody využití uhlíkových elektrod patří jejich nízká cena, široké potenciálové rozmezí, relativně inertní elektrochemické vlastnosti a elektrokatalytická aktivita pro mnohé redoxní reakce. V případě oxidačně-redukčních procesů organických a biologických molekul ve vodných i nevodných prostředích hrají uhlíkové elektrody na poli elektrochemických vlastností významnější roli než ušlechtilé kovy. Jejich využití je však známo i v oblastech výroby a zpracování kovů, uchovávání energie v bateriích a superkapacitorech a v mnohých dalších. V teoretické části jsou popsány vlastnosti uhlíkových materiálů využívaných jako elektrody. Dále je v práci popsána numerická metoda konečných prvků, která umožňuje aproximovat vlastnosti (v tomto případě elektrické) v částech vymezené oblasti, ve kterých je není možné získat přímým měřením. Na jejím základě je vytvořen 2D model rozložení potenciálu včetně pro dvě různé 2D geometrie. V rámci této práce jsou studovány dvě formy grafitu. Jeden blíže neurčený vzorek a vysoce orientovaný pyrolitický grafit (highly oriented pyrolitic graphite - HOPG). Podrobněji jsou zkoumány jejich elektrické vlastnosti a to metodou měření pomocí lineární čtyřsondy, která umožňuje zapojení elektrického obvodu i na malém vzorku a pomáhá tak zjistit měrnou vodivost (měrný odpor) těchto materiálů.
8
2. Teoretická část
2.1 Vlastnosti elektrod [1] Boronem dopovaný diamant, fullereny, a další nové uhlíkové materiály disponují odlišnými vlastnostmi v porovnání s grafitovým uhíkem běžně používaným z kraje 90. let minulého století (klasické uhlíkové materiály: grafit, skelný uhlík,...). Jejich vlastnosti umožňují také nové využití v oblastech snímání, elektrokatalýzy nebo třeba elektroniky. Rozdílnost uhlíkových elektrodových materiálů spočívá v jejich strukturním polymorfismu, chemické stabilitě, povrchové aktivitě a také v silné uhlíkové vazbě často přítomné mezi uhlíkovým materiálem a povrchovým modifikátorem.
2.1.1 Alotropy 2.1.1.1 Grafit a související sp2 hybridy
Nejjednodušším grafitovým materiálem je dvourozměrný grafen, grafenová deska. Nejuspořádanějším trojrozměrným grafitovým materiálem je tzv. HOPG (z angl. highly ordered pyrolytic graphite) a přírozený krystalický grafit. HOPG je vyroben při vysokoteplotní rozložením plynných uhlovodíků jako acetylen za vysokového tlaku. Grafitové materiály se často charakterizují pomocí rozměrů krystalu označovaných L a, pro velikost v rovině krystalu a L c, kolmo k rovivě grafenu. Hodnoty La a Lc pro HOPG a (natural single crystal) přesahují 1 (mikrometr), zatímco polykrystalický grafit, podobný hrotu tužky, bývá v rozsahu od 10 do 100 nm. Uspořádaná hexagonální rovina atomů obsahující osu „ a “ je známa jako bazální rovina (basal plane). Nepravidelný povrch paralelní s osou „ c “ je zase označován jako hrana (edge plane).
2.1.1.2 Uhlíkové vlákna
Jsou vyrobena z malých uhlovodíků, polymerů nebo chemickokatalytickou depozicí par. Průměr vláken se zpravidla pohybuje v rozsahu 5 – 50 μm. Grafitová osa „a“ je orientována ve směru délky vlákna a to díky výrobnímu procesu, který se snaží dosáhnout co největší síly v tahu. Uhlíkové vlákna se dělí do tří hlavních typů : 1. "Radial fibers" – grafenové plochy se rozšiřují ze středu vlákna 9
2. "Onion fibers" – soustředěné válce grafenových ploch 3. "Random fibers" – náhodně orientované grafitové plochy Příprava a morfologie "Onion fibers" je podobná jako u nanouhlíkových trubiček, kde však tyto nanotrubičky mají o hodně menší průměr. Přesto je rozdělení mezi velmi malé „onion fibers“ a MWCNT vytvořeno spíše uměle, ačkoliv z praktického hlediska jsou uhlíkové vlákna více předmětem manuální manipulace. 2.1.1.3 Glassy Carbon - skelný uhlík
Je vyráběn tepelným zpracováním různých polymerů, nejčastěji polyakrylonitrilů. Při teplotě 1000 – 3000 °C se struktura zbavuje heteroatomů. Uhlíkové vazby při těchto teplotách se neporuší a vytvoří grafitové plochy omezených rozměrů, kde L a a Lc jsou v rozsahu 30 – 70 Å. Prostor mezi planárními rovinami je větší než u HOPG, kolem 3,6Å. Hustota GC je zhruba 60 % úrovně HOPG a musí tak obsahovat hodně prázdných míst. 2.1.1.4 Diamant
Celková sp3 hybridizace a tetrahedrální vazby diamantu mají za následek tvrdost a nízkou vodivost, která činí jednoduchý i polykrystalický diamant nezajímavým pro elektrodový materiál. Nicméně, úmyslné zavedení příměsi ,v podobě bóru nebo dusíku, během depozice může zvýšit vodivost na dostatečnou úroveň. Příkladem je všeobecně známý mikrokrystalický bórem dopovaný diamant (Boron-doped diamond – BDD), který se vyrábí metodou depozice vrstev (vapor deposition) z H2 plasmy obsahující methan a zdroj bóru, nejčastěji B2H6. Bór vykazuje elektronový deficit vzhledem k uhlíku a je proto označován jako dopant typu „p“. Jeho výskyt ve výsledné struktuře je dosti vysoký, často v rozsahu 1018- 1020 atomů na cm3, nebo se uvádí jeho poměr k zastoupení uhlíku, tedy 10-5 – 10-3 . Vedle mikrokrystalického existuje také nanokrystalický diamant, který je vyráběn z CH 4/Ar plasmy. Tyto náhodně orientované krystaly mají rozměry několika desítek nanometrů. Přítomnost plynného dusíku nebo zdroje bóru během depozice vede k dopování typu „n“ nebo „p“. Poměr povrchu a objemu je zde výrazně vyšší než pro mikrokrystalický diamant. Povrch obsahuje význámné π-vazby a sp2 hybridizaci, což vede k vyšší elektrické vodivosti než u nedopovaného diamantu a to dokonce i v případě, kdy samotný nanokrystalický diamant není dopován.
10
2.1.1.5 Carbon Nanotubes
Rozdělujeme na jednostěnné (single-walled nanotubes, SWNT) a vícestěnné (multi-walled nanotubes, MWNT) .Nanotrubičky jsou vyráběny buď pomocí metody obloukového výboje (arcdischarge method), kterou se také produkují fullereny jako C 60 nebo C70, nebo chemickou depozicí par na kovových částečcích (chemical vapor deposition on metal particles – CVD), nejčastěji na Fe nebo Ni.
2.1.2 Elektronové vlastnosti uhlíkových elektrodových materiálů Důležitý význam z hlediska elektrodových materiálů má tzv. hustota stavů. (Density of electronic states – DOS). Na rozdíl od hustoty stavů pro kovy, se tvar i velikost rozdělení pro uhlíkové materiály značně mění se strukturou uhlíku. Orbitaly σ a π grafitového uhlíku tvoří obsazené valenční pásy, kdežto antivazebné orbitaly zahrnují vodivostní pás. Pro nekonečnou plochu graphenu a ideální krystalický grafit existuje překryv valenčních a vodivostních pásů, což vede k malé hustotě stavů v blízkosti Fermiho energie. Porucha v mřížce grafitu vede k mnoha defektním stavům s energiemi mezi vodivostním a valenčním pásem. HOPG je považován za polokov, zatímco neuspořádaný grafit vykazuje elektronické chování jako kov s nízkou hustotou stavů. Krystalický diamant je pro změnu polovodič se širokým pásem zakázaných energií větším než 6 eV v jedno krystalu. Dodatečných energiových stavů v tomto pásmu je dosaženo dotováním pomocí Bóru nebo dusíku. Nanokrystalický diamant má dostatečně velký počet defektů a je proto mhohem více vodivý nežli samotný krystal diamantu i bez dotování. Uhlíkové nanotrubičky (Carbon nanotubes - CN) vykazují rozmanitost v rozdělení hustot stavů v závislosti na průměru trubiček. Polokovové trubičky nemají v zakázaném pásmu žádné elektronické stavy, zatímco kovové trubičky mají nenulovou hustotu stavů pod i nad Fermiho energií. Za celkem správné přirovnání lze považovat nepravidelný sp 2 uhlíkový materiál co by polykrystalický grafit nebo skelný uhlík a sp 3 uhlík jako silně bórem dotovaný diamant, který vede elektrony kvalitativně podobně jako kovy.
2.1.3 Spektroskopické vlastnosti Uhlíkové materiály jako grafit a GC absorbují v širokém energetickém spektru, od nízké UV až po 11
radiové frekvence. Uhlík, co by předmět technologického a vědeckého zájmu, byl zkoumán mnoha formami analytické spektroskopie, jako Ramanova a infračervená spektroskopie, (X-ray photoelecton spectroscopy XPS) rentegenová fotoelektronová spektroskopie, UV – vis, NMR, etc. Ramanova spektroskopie se ukázala být mimořádně přínosná vzhledem ke všem uhlíkovým materiálům používaným v elektrochemii. Dokonalý grafit má dva ramanovsky aktivní módy, a to při 47 cm -1 (E2g(1)) a 1582 cm-1 (E2g(2)). Grafit je velmi polarizovatelný a mód při 1582 cm-1 je hlavním rysem grafitových materiálů. Ačkoliv, je mód sám o sobě docela silný, hloubka průniku laseru a Ramanovského rozptylu fotonů je v grafitu krátká (~20 nm), proto je Ramanovský signál slabší než už diamantu. Příkladem může být bazální rovina HOPG. Ostrá E2g fononová vibrace při 1582 cm-1 je jediný prominentní rys v této oblasti. Jestliže je grafit narušen, objeví se nový pás v blízkosti 1360 cm -1. Navíc, poruchy způsobují rozšiřování a mírný posun pásu 1582 cm-1 směrem k vyšším frekvencím, stejně jako změny znaků vyšších řádů . Pásy odpovídající ~1360 a ~1582 cm-1 nedokonalých grafitových materiálů obecně odpovídají „D“ (z angl. disorder) a „G“ (graphite) pásům. „D“ pás byl dříve chybně připsán vibracím podél hran (edges) grafitových krystalů, protože menší rozměr La poskytoval větší poměr intenzity D/G. D pás byl také pozorován v HOPG substitučně dopovaných atomy bóru. Proto je správnější rozhodnout že D pás pochází z prolomení symetrie v grafitové mřížce a výsledek ovlivňuje výběrové ramanovské pravidla. D pás také vykazuje neobvyklý posun frekvence s vlnovou délkou laseru. V souvislosti s elektrochemií, je D pás užitečný k pozorování změn ve struktuře elektrodového povrchu s rozdílnou přípravou jako je třeba laserová aktivace nebo elektrochemická oxidace. Spolu se separací grafitových vrstev, se G pás rozděluje do dvou, 1582 a 1620 cm -1 , s intenzitami závisejícími na stupni interkalace iontů mezi grafitové vrstvy. Při rozpadu grafenových ploch se objevuje D pás. Pro méně uspořádané grafitové materiály než GC, se D a G pás rozšiřují a postupně se začínají spojovat. Přestože tvar a relativní velikost D a G pásu závisí na mikrostruktuře uhlíku, jejich šířka a překrytí způsobují analýzu velmi složitou. Uhlíkové vlákna, (carbon black), polykrystalický grafit vykazují D a G pásy zobrazující stupeň neuspořádanosti stejně jako teplotní historie. Mikrokrystalický diamant ukazuje prominentní phononový pás při 1332 cm-1 , který se nachází v samotném krystalu diamantu. Šířka pásu a změny v ostatní rysech ramanova spektra diamantu jsou
12
dobrými indikátory čirosti a čistoty, což je užitečné pro diagnostiku přípravy diamantových filmů. Diamantový pás 1332 cm-1je široce používán k vyhodnocení relativního množství sp3 a sp2 hybridizovaného uhlíku v diamantových vzorcích. Ramanovský typický vzorek pro sp2 uhlík je přibližně 50 násobek sp3, takže intenzita 1360 cm -1 vůči 1332 cm-1 je citlivý indikátor sp2 nečistot v syntetickém diamantu. Nanokrystalický diamant ukazuje velmi rozšířené spektrum, ve kterém je 1332 cm -1 pás část široké řady pásů rozpínající se od 1100 – 1600 cm -1 . Dopování borem v BDD vede k pozorovatelným změnám v symetrii 1332 cm-1 pásu, který tak může být použit k odvození stupně dopování. Uhlíkové nanotrubičky mají zajímavé Ramanovy spektra díky rezonančnímu efektu a elektronfononové interakci. G pás v oblasti 1600 cm-1 má podobný původ jako v případě HOPG, ale nápadně se mění v umístění a tvaru pro různé nanotrubičky. Nízká frekvence módů blízko 300 cm -1 je běžně připisována tzv. ring breathing modes (RBM) a vychází z vibrací trubičky. Protože nemají tyhle vibrace v grafitu své obdoby, využívá se této vlastnosti k rozlišení nanotrubiček od neuspořádaného grafitového materiálu, který často nanotrubičky doprovází co by nečistoty.
2.1.4 Elektrochemické vlastnosti Výběr materiálu elektrod a přípravy jeho povrchu je spjatý s konkrétním elektrochemickým parametrem, jako například redoxní potenciál. 2.1.4.1 Povrchová stuktura uhlíkových elektrodových materiálů
Uhlíkové materiály mají složitější povrchové složení než je tomu u kovů, a to nejen kvůli základní široce se měnící mikrostruktuře, ale také díky možnosti tvořit množství různých povrchových vazeb a funkčních skupin. Terminace
Když dojde k narušení mikrostruktury povrchu, reakce probíhající v okolní kapalině nebo plynu vedou k terminaci. Pro HOPG je jeho bazální rovina , paralelní k ose „a“, atomicky uspořádaná a nereaguje se vzduchem, jen při vysokých teplotách. Tato rovina je však hladká jen v rámci několika málo desítek nanometrů, a defektům je těžké se vyhnout v oblastech větších než pár mikrometrů. To v případě elektrochemie nestačí, a bazální rovina tak bude obsahovat defekty v podobě schodků. Hustota těchto defektů závisí na historii a přípravě materiálu. Opakem bazální roviny je rovina hran, pro HOPG je drsná a rozstřepená. 13
Když se během rozštěpení nebo leštění tvoří nezaplněné valence grafenových rovin z hran, reagují pak s kyslíkem a vodou a vytváří funkční skupiny obsahující kyslík. Nejběžnější procedurou přípravy uhlíkových elektrod je leštění, obvykle s brusivem jako diamand nebo oxid hlinitý. Bohužel právě obrušování/leštění přináší velké množství povrchových reakcí se vzduchem nebo vodou. V případě GC je možností jak se vyhnout porušení povrchu zapouzdření do epoxidu a ulomenou část takovéto tyčinky vystavit jen elektrolytickému roztoku. Odlomený GC povrch vykazuje neobvykle vysokou electrochemickou reaktivitu a relativně nízké překrytí oxidy. Terminace pomocí vodíku vede k docela stabilnímu povrchu. Vodíkové plasma uzavírá povrch GC C-H vazbami a navíc odstraňuje většinu povrchových oxidů. Ačkoliv je „H-terminace“ efektivní pro stabilizaci povrchu GC jako i jiných forem sp2 uhlíkových materiálů, vyžaduje drahou a neskladnou vakuovou aparaturu a není tak běžně používána. Z BDD se v prostředí H 2 plasmy tvoří automaticky a „H-terminovaný“ diamantový povrch, který s kyslíkem reaguje pomalu. Uhlíkové nanotrubičky jsou tvořeny různými druhy zakončení v závislosti na jejich přípravě. Stěny válců nanotrubiček podobné bazální rovině HOPG nejsou relativně reaktivní, zato jejich konce jsou předmětem oxidací a dalších reakcí.
2.1.5 Příprava povrchu uhlíkových elektrod Povrchová příprava je samozřejmostí pro heterogenní proces k jakému v elektrochmii dochází. Tři alotropy uhlíku, široká variace povrchových struktur a funkčních skupin stejně jako náchylnost uhlíku k adsorpci náhodných nečistot činí přípravu povrchu velmi komplikovanou, ale v mnoha případech může vést k dosažení selektivity nebo modifikaci elektrodové kinetiky. 2.1.5.1 Leštění, broušení, čištění
Leštění se k úpravě povrchu uhlíkových materiálů používá už dlouhou dobu a zatím zůstává nejběžnější metodou úpravy uhlíkových elektrod v elektrochemii. Snad neběžněji používaným uhlíkovým elektrodovým materiálem je skelný uhlík. Leštění pomocí brusného papíru s karbidem křemíku vede k relativně drsnému a reaktivnímu povrchu. V případě mikroelektrod z uhlíkových 14
vláken je hrubozrnné leštění adekvátní a účelné, potenciální znečištění je zde minimalizováno. Osvědčeným dobrým postupem je také provést ultrazvukové čištění po leštění, ovšem s apelací na čistotu dané kapaliny. Leštění obecně vede ke zbytkům jak z leštícího/brusného materiálu tak z uhlíkových částic, což je těžko odrstranitelné a může nepříjemně působit na elektrochemickou reaktivitu.
2.1.5.2 Aktivace uhlíkových elektrod
Termín „aktivace“ se používá k popisu postupu při němž se zvyšuje reaktivita uhlíkových elektrod, často za účelem detekce specifického analytu.
Povrchové čištění
Mnohé z aktivačních postupů jsou ve skutečnosti čistící kroky použité k odstranění povrchových vrstev, které jsou ovlivněny zacházením nebo znečištěním použitým
materiálem k čištění.
Metodou, která dosahuje dokonalosti čištěného povrhu je ulomení části elektrody, v případě běžně používaných GC elektrod však není praktická k běžnému používání. Přijatelně úspěšnou alternativou pro uhlíkové vlákna lze dosáhnout pomocí oříznutí bezprostředně před měřením. Aktivace laserem
V roce 1984 byla poprvé vydána zpráva o překvapivém účinku ozáření ulíkové elektrody impulsem laseru přímo ve zkoumaném elektrolytu. Postup obvykle zahrnuje 7-20 ns laserový impuls z Nd:YAG (1064 nm) nebo dusíkový (337) nm laseru s hustotou výkonu 10 – 100 MW/cm 2 . Prudká tepelná expanze a lokální ohřev způsobují desorpci adsorbátu, a při vyšší hustotách výkonu také narušení uhlíkové mikrostruktury. Pro GC a uhlíkové vlákna hodnoty do 25 MW/cm 2 nepůsobí pozorovatelné morfologické změny v povrchový strukturách, přesto způsobují výrazné zrychlení elektron-transferové kinetiky. Výhody použití laserové aktivace v elektrochemických a analytických aplikacích: –
použití in situ vedením laseru skrze celu, či jiné prostupné prostředí, odstraňuje potřebu vyjmutí elektrody nebo rozbrání cely
–
rychlost a opakovatelnost provedení vede k obnovitelnému povrchu uhlíkové elektrody což je vhodné pro pulsní voltametrii
–
in situ deaktivace po znečištění
–
aktivace mikroskopických nebo nerovinných elektrod nepřístupných pro jiné čištění jako jsou uhlíkové vlákna a mikrodisky
15
2.2 Maxwellovy rovnice pro makroskopické elektromagnetické pole [2] div D=ρ f div B=0 rot E=−
∂B ∂t
rot H =J f + D– ρf
(1) ∂D ∂t
elektrická indukce
[
C N ; ] 2 m Vm
objemová hustota volných nábojů (nezahrnuty vázané náboje)
[
C ] 3 m
[T ;
Wb Vs ; 2] 2 m m
B–
magnetická indukce
E–
intenzita elektrického pole
[
V N ; ] m C
H–
intenzita magnetického pole
[
A ] m
Jf –
hustota volných proudů (nezahrnut vázaný proud)
[
A ] m2
Určení vztahu mezi elektrickou induckí D a intenzitou elektrckého pole E , a také mezi intenzitou magnetického H a magnetickou indukcí B určují tzv. stavové rovnice (materiálové vztahy, které zahrnují vázané náboje a proudy. D( r ,t )=ϵ 0 E (r ,t )+P (r , t) (2)
1 H (r , t)= μ0 B( r , t)−M (r , t) P–
polarizace
M–
magnetizace
ε0 –
permitivita vakua
16
[
F ] m
μ0 –
[
permeabilitia vakua
H ] m
Vektory polarizace a magnetizace jsou obecně dány vztahy 3
P (r , t)=ϵ0∫ d r ' dt ' χ ̂elec (r , r ' ,t , t ' ; E) E( r ' ,t ' ) 1 ̂ (r , r ' , t , t ' ; B) B(r ' , t ' ) M (r , t)= μ 0 ∫ d 3 r ' dt ' χ magn
,
(3)
kde funkce permitivity a permeability jsou nahrazeny integrálem přes obecnější elektrickou a magnetickou susceptibilitu. rovnice kontinuity : J tot =J +
Jtot -
17
∂D ∂t
celková hustota proudu včetně posuvných proudů
(4)
2.3 Metoda konečných prvků - Teorie Modelu [3] Řešení maxwellových rovnic v diferenciálním tvaru je složité a náročné na použitý matematický aparát. Zpravidla se v praxi setkáváme se známým rozložením zdrojů, tedy s potenciály a s proudy zdrojů. Existuje řada metod řešení těchto rovnic, které lze obecně rozdělit na analytické a numerické. Mezi obecné analytické metody patří separace promněnných nebo převedení diferenciálních rovnic na intergální a jejich následné řešení. Mezi numerické metody patří metoda konečných diferencí, metoda konečných prvků a metoda hraničních prvků. Společným rysem numerických metod je náhrada přesného řešení diferenciální nebo integrální rovnice přibližným řešením. Takové řešení nalezne např. v elektrostatické úloze potenciál nebo hustotu náboje ve vybraných bodech prostorové sítě. Pole mezi body pak na získáme vhodnou aproximací. Přesnost řešení záleží na hustotě sítě. Výpočet se výrazně zjednoduší, pokud lze některou ze souřadnic nebo závislost pole na čase zanedbat. Příkladem takového zjednodušení je závislost pole jen na jedné prostorové souřadnici a nezávislost na čase, soustava se poté rozpadne na obyčejné diferenciální rovnice, které lze snadno integrovat. Podle počtu proměnných se rozlišují úlohy jedno-, dvou-, a trojrozměrné, pomocí známých zkratek 1D, 2D, 3D („D“ z angl. „dimense“ - rozměr). Analytické řešení konkrétních úloh lze nalézt jen pro některé 1D a 2D úlohy. Složitější problémy se řeší výhradně numericky. Všechny metody řešení jsou však závislé na zjednodušení obecné soustavy.
2.3.1 Metoda konečných prvků (MKP) Metoda konečných prvků je účinná metoda k řešení všech okrajových úloh inženýrské praxe, popsaných diferenciálními rovnicemi. Principem metody je vavedení uzlů a uzlových potenciálů v oblasti kde se počítá pole. Uzly zde mohou být rozloženy nerovnoměrně a tak se přizpůsobovat tvarům hraničních ploch. V místech, kde se očekává náhlá změna pole, se zavádí vetší hustota sítě. K jednotlivým uzlům se sestaví soustava rovnic pro neznámé uzlové potenciály. Koeficienty matice soustavy se počítají jako integrály přes elementární plošky nebo objemy, v jejichž vrcholech jsou právě uzly. Tyto elementární útvary nazýváme konečné prvky. Pro 2D systémy se používají konečné
18
prvky tvaru trojúhelníku a čtyřúhelníku. Složitější prvky se v rovině nepoužívají. Pro 3D systémy se používají tvary čtyřstěnu, pětistěnu a šestistěnu. Jednotlivé prvky mohou mít uzly i ve středu hran. Postup při použití MKP je následující : 1. Generace sítě prvků s uzly 2. Aproximace potenciálu na jednotlivých prvcích z uzlových hodnot 3. Dosazení zvolené aproximace do diferenciální rovnice nebo jejího ekvivalentu a sestavení soustavy rovnic pro neznámé uzlové hodnoty 4. Vyřešení soustavy 5. Zpracování dodatečných požadavků – výpočet dalších veličin a zobrazení výsledků
2.3.1.1 Generace sítě prvků
Na úsečkách (větvích) uzavírajících zkoumanou oblast se zavádí uzlové body (počáteční, koncový), počet uzlů na větvi a parametr, který předepisuje nelinearitu dělení. Počet uzlů na jednotlivých větvích může být různý. Nelinearita dělení se volí tak, aby v okolí hrany vnitřní elektrody, kde se očekává největší intenzita, bylo dělení nejhustší. Volba okrajových větví je přizpůsobena možnosti pozdějšího přiřazení okrajových podmínek. Pro některé uzly nebo větve je tak určen potenciál jehož hodnota není ovlivněna vnitřním děním systému. Takovéto podmínce se říká „Dirichletova“ a lze si ji představit jako kontakt se měřící elektrodou. Jiným typem podmínky je tzv. „Neumannova“, kterou charakterizuje normálová derivace proměnné, v některých bodech. Jde o určení toku neznámé veličiny hranicí a častým případem je nulová Neumannova podmínka u n≡∇ u . ⃗n=0 , tedy izolace dané oblasti. Vystupují-li podmínky obě, mluvíme o smíšené úloze. Toto třídění je však používanější v matematické fyzice, v inženýrské praxi jsou skutečné okrajové podmínky obecně složitější. Mohou zahrnovat kombinace potenciálu, normálové derivace potenciálu, ale i jejich integrály. Generátory 2D sítí jsou poměrně jednoduché, oproti tomu generace 3D úlohy bývá velmi náročná. Modelování sítě v této práci byl využit program Gmsh. [7] Tento volně šiřitelný generátor sítí konečných prvků se skládá z modulů pro sestavení geometrie modelu, tvorbu sítě, výpočetním zařízením (solver) a následným zpracováním (post-processing). K potřebám této práce postačí první dva. Ovládání programu je umožněno jak pomocí graficky uživatelského prostředí tak pomocí příkazů v textovém režimu, oba kroky se generují současně, od jednoho k druhému přecházet dle potřeby. 19
2.3.1.2 Aproximace potenciálu na jednotlivých prvcích z uzlových hodnot
Metoda konečných prvků je postavena na jednoduchém avšak ne zcela obvyklém principu aproximace hledané funkce. Využívá se co nejnižšho stupně aproximačního polynomu, který po dosazení do diferenciální rovnice představuje ještě netriviální řešení. Pomocí integrace per partes lze snížit řád rovnice o jeden a toho se využívá i u parciální diferenciální rovnice. Proto postačí rovnice 2.řádu s lineární aproximací.MKP tedy nevyužívá polynomy vyšších řádů na dlouhém intervalu, ale naopak na mnoha malých intervalech lineární nebo nejvýše kvadratickou aproximaci. Postup lze demonstrovat na úloze v rovině.
Postup aproximace ve 2D prostoru
Obvykle se využívá lineárního trojúhelníku, který lze charakterizovat uzly u1( x 1, y 1 ) , u2( x 2, y 2) , u3( x 3, y 3) . Na trohúhelník aproximujeme potenciál lineární funkcí ϕ= Ax+By+C
(5)
(e) Zavedeme funkce určující tvar aproximace, tzv. tvarové funkce N (e) 1 ( x , y) , N 2 ( x , y) (e) N 3 ( x , y) , které jsou rovny 1 v příslušném uzlu a nule v sousedních dvou uzlech. Nulové jsou i
mimo prvek (e). Pro kubické a vyšší stupně tvarových funkcí má řešení tendenci oscilovat, proto je nepoužíváme. Potenciál na prvku je poté dán uzlovými potenciály (e)
(e)
(e)
(e)
(6)
ϕ =ϕ1 N 1 ( x , y)+ϕ2 N 2 ( x , y)+ϕ3 N 3 ( x , y) K odvození tvarové rovnice
(e) N 1 ( x , y) využijeme rovnice roviny
(e)
(7)
N 1 =ax+bx +c která prochází body (x 1, y 1,1) ,( x 2, y 2,0 ) ,( x 3, y 3,0) . Řešení pak nalézáme v soustavě ax 1+by 1+c=1 x1 y1 1 a 1 ax 2 +by 2+c=0 nebo x 2 y 2 1 b = 0 ax 3 +by 3+c=0 x3 y3 1 c 0 Determinant soustavy je roven dvojnásobku plochy trojúhelníkového prvku
(8) S Δ tedy Δ=2SΔ .
Kladnou hodnotu determinantu zavedeme číslováním uzlů v kladném smyslu. Koeficienty 20
vypočítané ze soustavy dosadíme do rovnice roviny a dostaneme tak tvarovou funkci
N (e) 1
trojúhelníkového prvku s vrcholy 1,2,3 N (e) 1 =
1 [( y 2− y 3) x+(x 3− x 2) y+x 2 y 3−x 3 y 2 ] 2SΔ
Cyklickou záměnou indexů 1 – 2 – 3 – 1 na prvku dostaneme zbylé tvarové funkce
(9) (e) N 2 ( x , y) ,
(e) N 3 ( x , y) .
Ze všech tvarových funkcí uzlu j lze setrojit aproximační funkci tohoto uzlu jako součet všech tvarových funkcí uzlu j, tedy N j =∑ N j
(e)
Pj
(10)
kde Pj je počet prvků se společným uzlem j. Sestavíme-li aproximační funkci pro každý uzel sítě, dostaneme globální aproximaci potenciálu v celé oblasti NU
ϕa =∑ ϕ j N j ( x , y)
(11)
j=1
kde NU je počet uzlů.
Aproximace potencialu ve 3D úloze
Nejjednodušším prvkem je lineární čtyřstěn, který je určen čtyřmi vrcholy a lze na něm zavést lineární tvarovou funkci se čtyřmi konstantami (e)
N j =ax+by +cz +d
(12)
Stejným postupem jako u výše zmíněného lineárního trojúhelníku se odvodí konstanty a, b, c,d, které jsou řešním soustavy rovnic . Determinant soustavy Δ je roven šestinásobku objemu čtyřstěnu. Hledaná funkce má pak tvar (e)
N1 =
1 (a x+b 1 y+c1 z+d 1) Δ 1
(13)
Potenciál na prvku se aproximuje z tvarových funkcí prvku e a z uzlových hodnot potenciálu (e) (e) (e) ϕ(e)=ϕ1 N (e) 1 ( x , y , z )+ϕ2 N 2 ( x , y , z )+ϕ3 N 3 ( x , y , z )+ϕ 4 N 4 ( x , y , z )
(14)
Aproximační funkce uzlu j je rovna součtu všech tvarových funkcí uzlu j N j ( x , y , z )=∑ N j (x , y , z ) (e)
Potenciál v celé oblasti se vyjádří pomocí takto sestavených aproximačních funkcí a uzlových hodnot ze vztahu 21
(15)
ϕa =∑ ϕ j N j ( x , y , z )
(16)
2.3.1.3 Sestavení soustavy rovnic
Kvůli sestavení rovnic pro neznámé uzlové potenciály ϕ j si zavedeme oblast Ω , která může být 1D (úsečka), 2D (rovina) nebo 3D (objem). Oblast je uzavřena hranicí Γ , kterou tvoří uzavřená plocha u 3D. Hranice se v jednodušších případech dá rozdělit na dvě části Γ=Γ e+Γ n
(17)
kde Γ e představuje elektrody se zadaným potenciálem. Hranice se snažíme volit tak, aby šlo například předpokládat, že části Γ n s jednotkovou vnější normálou u n tvoří siločáry E, a tedy E má jen tečnou složku ke hranici, pak platí E⋅u n=0 na Γ=Γ n tj.
E n =−grad ϕ⋅u n=−
∂ϕ =0 ∂n
(18)
MKP je vhodná k výpočtu polí v uzavřených oblastech. Výpočet otevřených úloh vyžaduje doplňující matematický aparát, který zajistí expanzi pole do neomezené oblasti. Oblast rozdělíme na NE prvků s NU uzly. Potenciál aproximujeme podle (16) NU
ϕ≈ϕa=∑ ϕ j N j ( x , y , z )
(19)
j=1
Jako příklad sestavení rovnic využijeme nyní poissonovu rovnici v následujícím tvaru div D=div ϵ E=div ϵ(−grad ϕ)=ρ
(20)
aproximujeme indukci NU
NU
j=1
j=1
D ≈ −ϵ grad ∑ ϕ j N j ( x , y ,( z ))=−∑ ϕ j ϵgrad N j (x , y ,( z )) Derivace se z rovnice odstraní aplikováním operátoru derivace grad na známou funkci
(21) Nj .
Protože aproximační funkce jsou součtem funkcí tvarových, derivují se právě ony. Za předpokladu konstantní permitivity na prvku a její nezávislosti na teplotě lze dosadit aproximaci z (21) do (20), tato však pro libovolné ϕ j není splněna přesně, ale vykazuje zbytek, tzv reziduum res( x , y , z , ϕ j) , které je funkcí souřadnic a uzlových potenciálů NU
div D−ρ=−div ∑ ϕ j ϵ grad N j ( x , y ,(z ))−ρ=res( x , y , z , ϕ1, ... ϕ NU ) j =1
22
(22)
Úkolem je nalézt takové hodnoty potenciálů, které dají minimální zbytek. Různé způsoby minimalizace zbytku dají soustavy NU rovnic pro hledané uzlové potenciály. Nejrozšířenější metodou v MKP je Galerkinova metoda, ve které se zbytek minimalizuje v okolí i-tého uzlu jeho aproximační funkcí
∫Ω res( x , y , z , ϕ1, ... ϕNU ) N i ( x , y , z ) d Ω=0
i=1, ... , NU
(23)
kde NU je počet uzlů jejichž potenciál není znám. Dosazením (22) do (23) dostaneme
∫ N i div D d Ω−∫Ω N i ρ d Ω=0
(24)
Ω
Úpravami mezi jejichž hlavní část patří snížení řádu derivace pomocí integrace per partes, kdy se operátor derivace přenáší z div D na Ni využitím vztahu N i div D=div(N i D)−grad N i . D
(25)
dosadíme za D aproximaci a získáme rovnici soustavy pro i-tý uzel NU
∑ ϕ j ∫Γ ϵ grad N i⋅grad N j ( x , y , (z ))d Ω=∫ N i ρ d Ω j=1
(26)
Ω
Není třeba psát rovnice pro uzly na Γ e , ale jen pro vnitřní uzly a uzly na Γ n kde potenciál neznáme. Nechť je počet uzlů s neznámým potenciálem NUI a celkový počet NU. Volbou i=1, ... , NU dostaneme přeurčenou soustavu NU rovnic K ϕ=F
(27)
K je čtvercová matice koeficientů NU x NU. Pro koeficienty k ij matice K vyplývá porovnáním (27) a (26) k ij =∫ ϵ grad N i⋅N j d Ω Ω
prvky
(28)
f i vektoru F na pravé straně jsou f i=∫ N i ρ d Ω Ω
(29)
Vektor uzlových potenciálů obsahuje nejprve neznámé potenciály vnitřních uzlů a poté potenciály uzlů na elektrodách, tedy [ϕ1, ... , ϕNUI , ϕ NUI +1 ,... , ϕNU ]T =[ϕx , ϕe ]T
(30)
Soustavu rozdělíme na submatice ve tvaru K xx K ex
23
K xe K ee
ϕx = F x ϕe F e
(31)
Protože ϕe známe , řešíme redukovanou soustavu, která již není singulární K xx ϕ x=F x −K xe ϕ e
(32)
zde první vektor vpravo je od zdrojů hustoty ρ a druhý od známých potenciálů eletrod ϕe . Aproximační funkce
N i , N j jsou nulové mimo prvky, které obsahují současně i-tý a j-tý uzel,
stejně tak i jejich derivace. Proto je součin v rovnici (28) nenulový jen tehdy, patří-li uzly i, j témuž prvku. Většina koeficientů je pak nulová, matice K je řídká a symetrická podle diagonály. Protože N j je vyjádřena na jednotlivých prvcích tvarovými funkcemi
N j , je třeba při numerickém
vyčíslení počítat integrály (28) a (29) na každém prvku samostatně. S použitím rov. (27) získáme příslušný koeficient jako součet příspěvků jednotlivých prvků k ij =∑ k ij
f i= ∑ f i
(e)
(e)
p
p
(33 a,b)
příspěvky jednotlivých prvků pak jsou (e)
(e)
(e)
(e)
(e)
k ij =k ji =∫ ϵ grad N i ⋅grad N j d Ω (e)
Ω (e)
(e)
f i =∫ ϵ ρ d Ω Ω
(34) (35)
(e)
2.3.1.4 Vyřešení soustavy
Soustava rovnic je řídká a dobře podmíněná, což znamená že i při velkém počtu rovnic je řešení stabilní. Přednost před klasickými eliminačními metodami se dávaá metodám iteračním. Eliminační metody se totiž používají pro soustavy nepřevyšující desítku tisíc rovnic. Výhodou iteračních metod je, že uchovávají jenom pole nenulových koeficientů.
2.3.1.5 Výpočet dalších veličin – postprocesor
Z uzlových potenciálů vypočteme potenciál na jednotlivých prvcích jako je (14) . Z něho lze pak vypočíst intenzitu na každém prvku , tedy E =−grad ϕ =−∑ ϕ j grad N j (e)
24
(e)
(e)
(36)
2.4 Měření Lineární čtyřsondou [4,5,6]
Čtyřsondu tvoří čtyři kovové hroty umístěné na přímce ve vzájemné vzdálenosti s.Hroty jsou přitlačeny na povrch vzorku s rovinnou plochou. Proud přivádíme do vnější dvojice hrotů a na vnitřní dvojici měříme napětí. Ve vzorku vzniká elektrické pole stacionárního proudu s kulovými ekvipotenciálami. Za předpokladu, že vzdálenost hranice vzorku od hrotů je mnohem větší než vzdálenost hrotů s, lze považovat objem vzorku považovat za polonekonečný.
Ilustrace 1: Uspořádání elektrod na lineární čtyřsondě přiložené k obecnému vzorku Po přiložení bodového kontaktu, kterým vtéká stacionární proud I, na rovinu poloprostoru s konstantní měrnou vodivostí, vzniká v této rovině elektrické pole stacionárního proudu. Vektorové pole hustoty proudu J(r) a také elektrické intenzity E(r) ve vzdálenosti r od kontaktu jsou radiální a ekvipotenciály pole jsou polokoule o poloměru r se středem v kontaktu. Vytéka-li jiným bodovým kontaktem z poloprostoru proud -I vznikne obdobné elektrické pole, samozřejmě s opačně orientovanou hustotou a elektrickou intenzitou vzhledem k polohovému vektoru r . Vtéká-li jedním kontaktem proud a druhým vytéká, tedy působí-li oba proudy současně, pak vzniklé pole stacionárního proudu je superpozicí obou polí výše zmíněných.
Závislost hustoty proudu na elektrické intenzitě plyne z diferenciálního tvaru Ohmova zákona
25
J =σ E kde
(37)
σ představuje měrnou vodivost. V případě, kdy má hustota proudu ve všech bodech plochy
stejnou velikost i směr, lze redukovat integrál proudu I i=∫S J d S na i
I =J⋅S
(38)
Proud vytéká z polokulovitého povrchu elektrody o obsahu plochy S=2 π r 1
2
(39)
Intenzita elektrického pole v místě kde proud vtéká do bude i
E (r 1)=
I 1 σ 2 π r 12
(40)
potenciál elektrického pole ve vzdálenosti r 1 od prvního kontaktu I 1 I ϕi (r 1)=−∫ E i ( r 1 )dr 1=−∫ dr 1= (41) 2 σ 2π r1 σ2πr1 potenciál pole ve vzdálenosti r 4 od čtvrtého kontaktu (kterým proud vytéká) je pak analogicky I ϕ f (r 4 )= (42) σ 2 π r4 Oba proudové kontakty vytvářejí na povrchu potenciál daný jejich superpozicemi I 1 1 i f ϕ=ϕ +ϕ = ( − ) σ 2 π r1 r4 Rozdíl potenciálů v bodě druhého a třetího kontaktu určuje napětí měřené sondou I 1 1 ϕ 2= ( − ) σ 2 π s 1 s 2+s 3 I 1 1 ϕ 3= ( − ) σ 2 π s 1+s 2 s 3 I 1 1 1 1 U =ϕ2 −ϕ3= ( + − − ) σ 2 π s1 s3 s 2+s3 s1+s 2 Měrnou vodivost pak lze vyjádřit následovně
(43)
(44) (45) (46)
I 1 1 1 1 ( + − − ) (47) U2 π s1 s3 s 2+s3 s1+s 2 V případě, kdy pro jednotlivé vzdálenosti mezi elektrodami platí s1 =s2 =s3=s , zredukuje se vztah pro měrnou vodivost na I σ= (48) U2 π s σ=
26
3. Praktická část 3.1 Ověření Ohmova zákona První část měření byla provedena pro ověření platnosti Ohmova zákona, tj. zjištění zda v oblasti proudů, kde se provádí měření, nezávisí odpor (měřný odpor) na velikosti proudu. K dispozici byl vzorek uhlíkového materiálu blíže nespecifikovaný a vzorek vysoce uspořádaného uhlíkového materiálu HOPG. Oba vzorky měly tvar kvádru s podstavou čtverce. Měření probíhalo na rovinné ploše vzorku, kterou bylo potřeba před samotným měřením uchytit tak, aby elektrody lineární čtyřsondy dosedly rovnoměrně na tuto rovinu a nedocházelo tak k poškození vzorku přítlakem sondy nebo nedokonalostem v kontaktech jednotlivých hrotů elektrody. Aretace vzorku vůči sondě byla provedena nastavením podstavce, který umožňoval naklonění roviny ve dvou na sebe kolmých osách. Nedokonalost kontaktu jednotlivých elektrod s povrchem vzorku oznamoval přístroj generující stejnosměrný proud na svém kontrolním panelu.
3.1.2 Neznámý vzorek uhlíkového materiálu Lineární čtyřsonda byla přikládána na rovinu vzorku se čtvercovou podstavou o hraně a = 5,0 mm. Pomocí nastavitelného podstavce byla rovina uvedena do vodorovné polohy a kolmo k čtyřsondě. Měřila se Voltampérová charakteristika vzorku ve čtyřech orientacích lineární sondy vůči vzorku. Při pojmenování čtverce podstavy ABCD byla provedena jednotlivá měření ve směru paralelním se stranami AB, BC a s úhlopříčkami BD a AC. Sonda byla vycentrována do středu vzorku. Hroty jednotlivých elektrod sondy byly od sebe vzdáleny s = 1,066 mm. Hodnoty měrné vodivosti byly vypočítány dle vzorce .. , měrný odpor je převrácená hodnota měrné vodivosti.
27
3.1.2.1 Sonda v poloze paralelní ke straně AB
I [mA] U(p1) [mV] U(p2) [mV] 5 0,136 0,136 10 0,271 0,272 15 0,405 0,406 20 0,539 0,539 25 0,674 0,674 30 0,809 0,810 35 0,943 0,943 40 1,077 1,078 45 1,213 1,212 50 1,345 1,346 55 1,481 1,481 60 1,615 1,615 65 1,749 1,749 70 1,883 1,883 75 2,017 2,016 80 2,151 2,150 85 2,283 2,283 90 2,419 2,418 95 2,552 2,551 100 2,685 2,684
σ [S/cm] 54,918 55,019 55,257 55,427 55,407 55,359 55,442 55,453 55,439 55,510 55,474 55,496 55,515 55,530 55,558 55,569 55,616 55,588 55,617 55,644
Tabulka 1: Hodnoty V/A charakteristiky a měrné vodivosti pro neznámý vzorek ve směru paralelním s AB. Hodnoty napětí p1 a p2 odpovídají opačným orientacím proudu. U [mV] 3,0 y = 0,0268x + 0,0041
2,5 2,0 1,5 1,0
A-B (p1) A-B (p2)
0,5
Lineární (A-B (p2))
0,0 0
10
20
30
40
50
60
70
80
90 100 I [mA]
Ilustrace 2: Graf V/A charakteristiky pro neznámý vzorek v orientaci sondy ve směru AB. Spojnice trendu byla přidána pro hodnoty napětí p2.
28
3.1.2.2 Sonda v poloze paralelní ke straně BC I [mA] U(p1) [mV] U(p2) [mV] 5 0,315 0,318 10 0,627 0,631 15 0,942 0,946 20 1,254 1,257 25 1,566 1,559 30 1,879 1,882 35 2,194 2,198 40 2,504 2,507 45 2,817 2,820 50 3,131 3,134 55 3,444 3,447 60 3,756 3,759 65 4,066 4,068 70 4,381 4,179 75 4,694 4,696 80 5,002 5,004 85 5,319 5,320 90 5,627 5,628 95 5,941 5,941 100 6,249 6,249
σ [S/cm] 23,599 23,749 23,736 23,796 23,900 23,830 23,808 23,848 23,849 23,843 23,845 23,853 23,874 24,444 23,862 23,886 23,869 23,890 23,886 23,904
Tabulka 2: Hodnoty V/A charakteristiky a měrné vodivosti pro neznámý vzorek ve směru paralelním s BC. Hodnoty napětí p1 a p2 odpovídají opačným orientacím proudu. U [mV] 7,0 6,0 5,0 4,0 3,0 2,0 1,0 0,0 0
y = 0,0623x + 0,0085
B-C (p1) B-C (p2) Lineární (B-C (p2))
10
20
30
40
50
60
70
80
90
100
I [mA]
Ilustrace 3: Graf V/A charakteristiky pro neznámý vzorek v orientaci sondy ve směru BC. Spojnice trendu byla přidána pro hodnoty napětí p2.
29
3.1.2.3 Sonda v poloze paralelní k úhlopříčce BD
I [mA] U(p1) [mV] U(p2) [mV] 5 0,227 0,234 10 0,458 0,463 15 0,688 0,693 20 0,915 0,919 25 1,144 1,150 30 1,376 1,380 35 1,606 1,610 40 1,833 1,838 45 2,064 2,069 50 2,294 2,298 55 2,523 2,526 60 2,751 2,756 65 2,979 2,984 70 3,210 3,214 75 3,441 3,456 80 3,681 3,685 85 3,908 3,912 90 4,139 4,143 95 4,367 4,371 100 4,597 4,601
σ [S/cm] 32,410 32,439 32,450 32,580 32,558 32,520 32,514 32,553 32,528 32,530 32,544 32,550 32,566 32,554 32,487 32,447 32,473 32,465 32,481 32,480
Tabulka 3: Hodnoty V/A charakteristiky a měrné vodivosti pro neznámý vzorek ve směru paralelním s úhlopříčkou BD. Hodnoty napětí p1 a p2 odpovídají opačným orientacím proudu. U [mV] 5,0 y = 0,046x + 0,0002
4,0 3,0 2,0
B-D (p1)
1,0
Lineární (B-D (p2))
B-D (p2)
0,0 0
10
20
30
40
50
60
70
80
90 100 I [mA]
Ilustrace 4: Graf V/A charakteristiky pro neznámý vzorek v orientaci sondy ve směru úhlopříčky BD. Spojnice trendu byla přidána pro hodnoty napětí p2.
30
3.1.2.4 Sonda v poloze paralelní k úhlopříčce AC
I [mA] U(p1) [mV] U(p2) [mV] 5 0,276 0,282 10 0,553 0,559 15 0,829 0,833 20 1,102 1,107 25 1,377 1,382 30 1,654 1,655 35 1,920 1,926 40 2,196 2,202 45 2,472 2,480 50 2,749 2,754 55 3,025 3,029 60 3,300 3,303 65 3,574 3,579 70 3,847 3,854 75 4,122 4,128 80 4,395 4,402 85 4,671 4,678 90 4,944 4,951 95 5,219 5,226 100 5,490 5,496
σ [S/cm] 26,773 26,867 26,963 27,049 27,071 27,086 27,188 27,172 27,149 27,145 27,141 27,147 27,148 27,156 27,159 27,169 27,162 27,173 27,172 27,194
Tabulka 4: Hodnoty V/A charakteristiky a měrné vodivosti pro neznámý vzorek ve směru paralelním s úhlopříčkou AC. Hodnoty napětí p1 a p2 odpovídají opačným orientacím proudu. U [mV] 6,0 5,0
y = 0,0549x + 0,0081
4,0 3,0 2,0
A-C (p1) A-C (p2)
1,0 0,0
Lineární (A-C (p2))
0
10
20
30
40
50
60
70
80
90 100 I [mA]
Ilustrace 5: Graf V/A charakteristiky pro neznámý vzorek v orientaci sondy ve směru úhlopříčky AC. Spojnice trendu byla přidána pro hodnoty napětí p2.
31
Vzhledem k tomu, že se hodnoty měrné vodivosti pro obě polarity zdroje významně nelišily, byla spočítána jedna podle rov (47) jako aritmetický průměr: směr AB
σ=(55,44±0,04)
S cm
směr BC
σ=( 23,86±0,04)
S cm
směr BD
σ=( 32,50±0,01)
S cm
směr AC
σ=( 27,10±0,03)
S cm
Měrný odpor určíme jako převrácenou hodnotu měrné vodivosti: směr AB
ρ=(0,01804±0,00001)Ω. cm
směr BC
ρ=( 0,04191±0,00006)Ω . cm
směr BD
ρ=( 0,03076±0,00001)Ω . cm
směr AC
ρ=(0,03690±0,00003)Ω. cm
Odchylky byly určeny statistickým zpracováním dat při hladině pravděpodobnosti 0,6827.
32
3.1.3 HOPG Stejným způsobem jako v případě neznámého vzorku byl proměřen také HOPG v rozsahu 100 mA. Rozměry tohoto kvádru s podstavou čtverce: 10mm x 9,9mm x 2,5mm. I [mA] 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100
U(p1) AB U(p2) AB U(p1) BC U(p2) BC [mV] [mV] [mV] [mV] 0,013 0,018 0,015 0,014 0,029 0,033 0,030 0,029 0,044 0,049 0,044 0,043 0,059 0,064 0,059 0,059 0,075 0,080 0,073 0,073 0,091 0,095 0,087 0,088 0,107 0,111 0,102 0,103 0,122 0,126 0,117 0,117 0,138 0,141 0,131 0,132 0,154 0,157 0,146 0,146 0,169 0,172 0,161 0,161 0,185 0,187 0,175 0,175 0,200 0,203 0,190 0,190 0,216 0,218 0,205 0,204 0,231 0,234 0,219 0,219 0,247 0,249 0,234 0,233 0,259 0,262 0,248 0,248 0,274 0,277 0,263 0,262 0,290 0,293 0,278 0,277 0,305 0,308 0,292 0,292
Měrná vodivost dle vztahu (47) : směr AB
σ=( 482±3)
S cm
směr BC
σ=(511±1)
S cm
Měrný odpor určíme jako převrácenou hodnotu měrné vodivosti:
33
směr AB
ρ=( 0,002071±0,000003)Ω . cm
směr BC
ρ=( 0,001956±0,000002)Ω. cm
σ AB [S/cm] 481,861 481,861 481,861 485,778 481,861 481,861 479,650 481,861 481,861 480,311 481,861 481,861 481,861 481,861 481,861 481,861 487,410 487,982 486,820 487,363
σ BC [S/cm] 515,092 506,362 515,092 506,362 511,564 512,149 510,067 510,690 511,175 511,564 510,293 512,149 511,026 511,314 511,564 511,783 511,977 512,149 511,380 511,564
I [mA] 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100
U(p1) AC [mV] 0,015 0,031 0,047 0,064 0,080 0,096 0,112 0,129 0,145 0,161 0,178 0,194 0,210 0,226 0,243 0,259 0,275 0,291 0,308 0,324
U(p2) AC [mV] 0,018 0,034 0,050 0,066 0,083 0,099 0,115 0,131 0,148 0,165 0,180 0,196 0,212 0,229 0,245 0,261 0,277 0,293 0,309 0,326
Měrná vodivost dle vztahu (47) :
U(p1) BD [mV] 0,013 0,027 0,041 0,055 0,069 0,084 0,098 0,113 0,126 0,141 0,155 0,169 0,183 0,198 0,212 0,226 0,239 0,254 0,268 0,282
směr AC
σ=(459±1)
S cm
směr BD
σ=(526±1)
S cm
U(p2) BD [mV] 0,016 0,030 0,045 0,059 0,073 0,087 0,100 0,114 0,128 0,142 0,156 0,170 0,184 0,198 0,212 0,225 0,240 0,253 0,267 0,282
σ AC [S/cm] 452,657 459,621 461,990 459,621 458,211 459,621 460,633 459,621 458,837 457,650 458,979 459,621 460,165 459,621 459,150 459,621 460,037 460,408 459,993 459,621
σ BD [S/cm] 515,092 524,129 521,082 524,129 525,975 524,129 528,100 526,438 529,288 527,833 528,342 528,767 529,128 528,100 528,456 529,940 530,147 530,332 530,497 529,705
Měrný odpor určíme jako převrácenou hodnotu měrné vodivosti: směr AC
ρ=( 0,002177±0,000002)Ω. cm
směr BD
ρ=( 0,001898±0,000003)Ω . cm
Výsledky pro HOPG byly určeny statistickým zpracováním dat při hladině pravděpodobnosti 0,6827.
34
3.2 Závislost velikosti měrné vodivosti na směru lineární sondy nad vzorkem. Druhá část měření zjišťovala závislot velikosti měrného odporu na orientaci posazení lineární sondy vůči povrchu vzorků. Neznámý vzorek i HOPG byli nastaveni na aretačním stolečku tak aby je bylo možné otáčet kolem pomyslené svislé osy procházející středem čtvercového tvaru, kde dosedal také střed lineární čtyřsondy. Aretační stoleček byl posazen do jiného uchycení, které umožňovalo otočení kolem svislé osy o 360° a zároveň daný úhel odečíst. Měření bylo provedeno pro každý vzorek vícekrát, neboť v jeho průběhu přitláčením sondy k povrchu docházelo k vadám materiálu a měření tak bylo ovlivněno. Měřilo se napětí v obou polaritách zdroje pro dvě hodnoty stejnosměrného proudu I=40 mA a I=70 mA v průběhu celých 360° a to po 5°. Jde tedy o značné množství dat pro které zde tabulku neuvádím, ale jejichž hodnota je patrna z následujícího grafu. Hodnoty měrného odporu byly počítáne pomocí jako převrácená hodnota rov (47).
ρ [Ω.cm] 0,045000 0,040000 0,035000 0,030000 0,025000 40 mA
0,020000
70 mA
0,015000
HOPG 40 mA
0,010000
HOPG 70 mA
0,005000 0,000000 0
30
60
90
120 150 180 210 240 270 300 330 360 uhel [°]
Z obrázku je patrno, že narozdíl od HOPG , neznámý vzorek vykazuje v určitých směrech větší měrnou rezistivitu a tudíž ani měrná vodivost nebude pro šíření proudu ve všech směrech konstantní. Materiál vykazuje určitou anizotropii.
35
3.3 Model stacionárního proudového pole [3] 3.3.1 Laplaceova rovnice Pole stacionární proudu je polem pohybujících se nábojů konstantní rychlostí. Nadále budeme uvažovat, že tento pohyb není ovlivněn magnetickým polem. Proud procházející vodičem libovolného průřezu v něm vytváří pole vektoru hustoty proudu. Proud I procházející plochou S je definován jako množství náboje prošlé touto plochou za jednotku času. I=
dQ dt
(49)
Pro vodič který vede proud platí stejně jako v elektrostatice, že výsledná hustota kladných a záporných částic uvnitř vodiče je nulová a že volný náboj smí sídlit jen na povrchu. Elektrické pole pak tyto náboje uvádí do pohybu a vzniká tak proud hustoty J. Elektrony se nemohou uvnitř vodiče shlukovat, protože by tak vytvořilz dodatečné pole. Pro libovolnou uzavřenou plochu uvnitř vodiče pak platí
∮S J . d S=0 div J =0
což odpovídá podmínce
(50) (51)
Hustota proudu závisí u některých materiálů lineárně na intenzitě elektrického pole ⃗ ⃗ J =γ E Pole ustáleného proudu je nevírové
(52)
∂B =0 , proto platí ∂dt
∮l E . d l=0
(53)
Stejně jako v elektrostatickém poli můžeme zavést skalární potenciál vztahem E=−grad ϕ , který splňuje rovnici
div γ grad ϕ =0
(54) (55)
Pomocí metody konečných prvků budeme tedy hledat řešení Laplaceovy rovnice
Δ ϕ =0
(56)
3.3.2 Upřesnění výpočtu MKP Rovnice
36
K ϕ=F uvedená v teoretické části bude mít vektor F roven nule a řešení soustavy bude
záviset na hodnotách potenciálu, který budeme považovat za předem stanovený, tedy na Dirichletově podmínce. Redukovaná soustava pak bude vypadat následovně K xx ϕ x=−K xe ϕ e ,
(57)
kde ϕe jsou známé díky okrajové podmínce. Pro výpočet matice K dle vzorce uvedeného v teorii budeme potřebovat znát konkrétní tvar aproximačních funkcí Ni . V matici K budou vystupovat diagonální členy Kii , které se složí z jednotlivých částí obsahujících vrchol i a další nenulové členy budou Kij , kde i a j náleží stejnému elementu. K ii =∑N
i
∫Ω grad N i grad N i dS
(58)
1 2 ∫e grad N grad N dS =( 2 S ) ∫e (( y l − y r )2+(x r −x l )2 ) dxdy e e i
e i
(59)
proměnné y a x s indexy jsou souřadnice patřící dvěma zbývajícím uzlům v trojúhelníkovém elementu kromě uzlu i pro nějž je tatot diagonalní matice počítána. 2.
1 I e=( ) (( y l − y r )2+(x r −x l )2 ) I (0,0) , kde 2 Se
(60)
I s koeficientem (0,0) je momentovým integrálem pro prvek e , tento momentový integrál ale také odpovídá ploše elementu e, kde vystupují souradnice vrcholů trojúhelníku 1 1 S e = (x i y j+ x j y k + x k yi )− ( xi y k +x j y i+x k y j ) 2 2
(61)
Obdobně pro nediagonální člen Kij platí K ij =∑ jϵ J ∫Ω grad N i grad N j dS , kde
(62)
J je hranicí podpory Nj. Příspěvek nediagonálního členu pak tvoří dva integrály přes přilehlé boční trojúhelníky
∑ j ϵ J ∫L +∫R j
()
j
()
grad N i grad N j
(63)
Nediagonální členy příspěvku od j a k uzlu v trojúhelníku pak Kij a Kik vypadají následovně 2.
I ij =(
1 ) (( y k − y j )( y i− y k )+( x j −x k )(x k −x i )) I (0,0) 2 Se
(64)
2.
1 I ij =( ) (( y k − y j )( y i− y k )+( x j −x k )(x k −x i )) I (0,0) 2 Se Koeficienty matice K vypočítá program matlab dle následujícího schématu : Pro každý prvek e ve vygenerované síti
37
(65)
pro každý vnitřní uzel i prvku e 1) Kii + = Ie 2) Kij + = Iij
pro vnejší uzel je Kij násoben potenciálem ϕe
3) Kik + = Iik
pro vnejší uzel je Kij násoben potenciálem ϕe
další vnitřní uzel prvku e Další prvek ... Skript programu matlab stejne jako popis zadaní pro vygenerování sítě konečných prvků je uveřejněn v příloze.
3.3.3 Modely elektrod Vygenerovaná síť, která nám poskytne jednotlivé uzlové body, tedy śouřadnice vrcholů trojúhelníků, byla připravena pro dva případy. V prvním případě se jednalo oveření aproximace rozložení potenciálu pro Laplaceovu rovnici mezi dvěma paralelními elektrodami. Síť byla vygenerována se 482 prvky při 271 uzlech.
Ilustrace 6: Síť prvků mezi paralelními eletrodami, ty leží mezi body 1-2 a 3-4 Za elektrody byly vybrány naproti sobě léžící úsečky mezi body 1-2 a 3-4. Na jednu elektrodu bylo nastaveno napětí 5 V a na druhou -5V. Graficke zobrazení modelu na této jednoduché síti se nachází na následujícím obrázku.
38
Ilustrace 7: Rozložení aproximovaného potenciálu je zde znázorněno stupni šedi mezi +5V a - 5V. Hodnoty aproximovaného potenciálu mezi elektrodami odpovídají předpokládanému rozložení dle Laplaceovy rovnice. Dalším otestováním modelu bylo vygenerování sítě simulující bližší podmínky měření, tedy přiložení dvou prozatím bodových elektrod k ploše materiálu. Počet prvků v dané síti je 1673 při 894 uzlových bodech. Vzorek byl pomyslně rozříznut v místě přiložení elektrod. Na obrázku jsou místa přiložení označena jako body číslo 5 a 8. Opět byla nastavena referenční hodnota napětí + 5V a – 5V. Protože se předpokládá větší spád potencálu na nejkratší cestě mezi body 5 a 8, byla zda nastavena jemnější síť , která směrem od těchto bodů řídne. Okrajové podmínky jsou opět udané jen hodnotou napětí na elektrodách.
39
Ilustrace 8: Síť konečných prvku a uzlových bodů pro případ přiložení dvou bodových eletkrod. Výsledný model elektrody byl zpracován jak ve 2D podobě, kde rozložení potenciálu tvoří škála stupně šedi spolu s ekvipotenciálními hladinami, tak ve 3D zobrazení, kde velikost napětí na jednotlivých bodech plochy x-y byla dána z-složkou.
Ilustrace 9: Rozložení potenciálu mezi 2 bodovými elektrodami. Stupnice šedi opovídá hodnotám ve V.
40
Ilustrace 10: Rozložení potenciálu mezi dvěma elektrodami. Z-souřadnice nahrazuje velikost potenciálu. Hodnoty ve stupních šedi jsou ve V. Z posledního obrázku je patrný velký potenciálový gradient mezi oběma body reprezentujícími přiložené elekterody k materiálu, ale také málé změny potenciálu na opačné straně této modelové plochy.
41
4. Závěr
První část práce zahrnovala výpis vlastností uhlíkových elektrod jak z heldiska strukturních tak elektrochemických. V další části je výpis maxwellových rovnic pro daný systém. V části zabývající se teorií modelu byl naznačen postup jakým metoda konečných prvků aproximuje v oblasti systému dle zadaných okrajových podmínek neznámé hodnoty. V praktické části jsme ověřili jak pro neznámý vzorek uhlíkového materiálu, tak pro HOPG, že měrná vodivost se nemění s velikostí procházejícího proudu. Dále jsme zjištovali zda velikost měrné vodivosti závisí na směru průchodu proudu. Měrná vodivost HOPG v různých směrech ukazuje spíše na izotropní vlastnosti tohoto materiálu. Naopak u zkoumaného vzorku byl zřetelně vidět rozdíl ve vodivosti v různých směrech jak dokládá graf této závislosti. Vypočítaná velikost měrné vodivosti pro HOPG vyplývající z našeho měření se pohybuje mezi 482 - 526 S/cm pro dané směry měření. Měrná vodivost neznámého vzorku vykazovala největší vodivost : σ=(55,44±0,04)
S v jednom ze směrů. A naopak nejmenší vodivost ve směru kolmém k výše cm
uvedenému , a to σ=( 23,86±0,04)
S . cm
V poslední praktické části byl vygenerován model založený na metodě konečných prvků. Tento model byl nastaven pro dva 2D systémy. Jeden pro paralelní elektrody v podobě přímek od sebe vzdálených, kde se potvrdilo rozložení potenciálu dle laplaceovy rovnice na se kterou model počítá. Druhý model simuluje přiložení dvou bodových elektrod. Pro toto nastavení byly zhotoveny dva obrázky, které srozumitelně popisují rozložení potenciálu v nadefinovaném prostředí.
42
Seznam použité literatury [1]
McCreery, R. L. : Advanced Carbon Electrode Materials for Molecular
Electrochemistry.
Chemical Reviews, ročník 108, č. 7, 2008: s. 2646 - 2687 [2]
Maxwell's equations in Wikipedia : the free encyclopedia [online]. Wikimedia Foundation, last modified on 7 May 2011 - [cit. 2011-05-08]. Anglická verze.
Dostupná z WWW:
. [3]
DĚDEK, Libor; DĚDKOVÁ, Jarmila. Elektromagnetismus. 2.vyd. Vysoké učení technické v Brně: Nakladatelství VUTIUM, 2000.232 s. ISBN 80-214-1548-7
[4]
SEDLÁK, Bedřich; ŠTOLL, Ivan. Elektřina a magnetismus. 2.opravené a
rozšířené
vyd. Praha : Academia, 2002. 632 s. ISBN 80-200-1004-1 [5]
Plaček, F.: Měrný odpor polovodičů. 2004
[6]
Valdes, L. B.: Proceedings of the I.R.E. 42 (1954), 420-427
[7]
C. Geuzaine and J.-F. Remacle. Gmsh: a three-dimensional finite element
mesh
generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering, Volume 79, Issue 11, pages
43
1309-1331, 2009
Příloha A Zadání generování sítě do programu Gmsh (příklad) 1) Pro model dvou bodových elektrod Point(1) = {0, 0, 0, 0.0003}; Point(2) = {0.005, 0, 0, 0.0003}; Point(3) = {0.005, 0.004, 0, 0.0003}; Point(4) = {0, 0.004, 0, 0.0003}; Point(5) = {0.001, 0.004, 0, 0.00005}; Point(6) = {0.002, 0.004, 0, 0.00005}; Point(7) = {0.003, 0.004, 0, 0.00005}; Point(8) = {0.004, 0.004, 0, 0.00005}; Line(1) = {5, 4}; Line(2) = {4, 1}; Line(3) = {1, 2}; Line(4) = {2, 3}; Line(5) = {3, 8}; Line(6) = {8, 7}; Line(7) = {7, 6}; Line(8) = {6, 5}; Line Loop(9) = {8, 1, 2, 3, 4, 5, 6, 7}; Plane Surface(10) = {9}; Physical Surface(11) = {10}; Physical Line(12) = {8, 7, 6}; Physical Line(13) = {1, 2, 3, 4, 5};
44
Skript programu matlab na výpočet matice K (příklad) 1) MaticeK pro model 2 bodových elektrod % MaticeK % Pro kazdy element T: % pro kazdy VNITRNI uzel i elementu T: %
1) Kii += prislusny intergral It
%
2)
%
a) je-li uzej j vnitrni: Kij+ = prislusny intergral Ik
%
je-li j vnejsi : 0 = (Fi)j*Ik nebo 0
%
b) je-li uzej k vnitrni: Kik+ = prislusny intergral Ij
%
je-li j vnejsi : 0 = (Fi)jkIj nebo 0
% dalsi vnitrni uzel elementu T % dalsi element. % format long; load vzorek-prvky.msh; P = [vzorek_prvky(:,1), vzorek_prvky(:,6),vzorek_prvky(:,7),vzorek_prvky(:,8)];
% seznam prvku a ktere uzly obsahuji
load vzorek-uzly.msh; M = [vzorek_uzly(:,2),vzorek_uzly(:,3)];
% seznam uzlu a ktere souradnice maji
% momentove integraly I = []; for p=1:1673;
% pro vsechny prvky
u1 = P(p,2); u2 = P(p,3); u3 = P(p,4); xa = M(u1,1); xb = M(u2,1); xc = M(u3,1); ya = M(u1,2); yb = M(u2,2); yc = M(u3,2); i00 = (1/2)*(xa*yb + xb*yc + xc*ya - xa*yc - xb*ya - xc*yb); i10 = (1/3)*(xa + xb + xc)*i00; i01 = (1/3)*(ya + yb + yc)*i00; i11 = (1/12)*(2*(xa*ya + xb*yb + xc*yc) + xa*(yb + yc) + xb*(ya + yc) + xc*(ya + yb))*i00; i20 = (1/6)*(xa^2 + xb^2 + xc^2 + xa*xb + xb*xc + xa*xc)*i00; i02 = (1/6)*(ya^2 + yb^2 + yc^2 + ya*yb + yb*yc + ya*yc)*i00; I = [I; [i00,i10,i01,i11,i20,i02]]; end; save moment-integraly.dat I -ascii; load moment-integraly.dat; I = moment_integraly;
% I = [i00, i10, i01, i11, i20, i02] % pro dany prvek "p" tedy plati
45
% i00 = I(p,1), i10 = I(p,2), % i01 = I(p,3), i11 = I(p,4), % i20 = I(p,5), i02 = I(p,6) K(894,894) = 0;
% ZAKLAD pro vyslednou obecnou matici Kij
for p = 1:1673;
% pro vsechny prvky
U = []; for c=[2 3 4];
% pozice sloupce pro jednotlive uzly tvorici prvek p
U=[U,P(p,c)];
% vytvori matici U obsahujici uzly patrici prvku p
end; for d = U;
% Algoritmus pro Kii
di = [];
% odebere ze seznamu uzlu pro prvek p prave ten uzel pro který se pocita Kii
for n = U; if n ~= d; di = [di, n]; end; end; xl = M(di(1,1),1); yl = M(di(1,1),2); xr = M(di(1,2),1); yr = M(di(1,2),2); St = I(p,1); It = (1/(4*St))*((yl - yr)^2 + (xr - xl)^2);
% Integral pro Kii daneho uzlu d
K(d,d) = K(d,d) + It;
% Pricteni vysledku It do obecne Matice Kij na patricne misto (tedy Kdd)
xi = M(d,1);
% Algoritmus pro Kij
yi = M(d,2); xj = M(di(1,1),1); yj = M(di(1,1),2); xk = M(di(1,2),1); yk = M(di(1,2),2); Sv = I(p,1); Iij = (1/(4*Sv))*((yk - yj)*(yi - yk) + (xj - xk)*(xk - xi)); K(d,di(1,1)) = K(d,di(1,1)) + Iij; % Algoritmus pro Kik Iik = (1/(4*Sv))*((yj - yk)*(yi - yj) + (xk - xj)*(xj - xi)); K(d,di(1,2)) = K(d,di(1,2)) + Iik; end; end; save MaticeK.dat K -ascii
2) Řešení soustavy
K xx ϕ x=−K xe ϕ e
% vyreseni soustavy K*Phi = 0 load MaticeK.dat; A = MaticeK;
% Matice vhodnejsi nazev pro upravy % A zustane jako zaloha pro porovnani zda
46
% byly zmeny provedeny spravne B = A; B ([1:894],[1:894]) = B ([1:4 6 7 9:894 5 8],[1:4 6 7 9:894 5 8]); % v B byly presunuty sloupce a radky na % konec na konec % tyto body odpovidajici mistum pro nez zname % potencial Kxx = B;
% Matice Kxx - pro nitrni uzly jejichz potencial
Kxx ([893 894],:) = [];
% neni znam. Matice Kxx je symetricka 892x892.
Kxx (:,[893 894]) = []; Kxe = B;
% Matice Kxe - 892x2
Kxe (:,[1:892]) = []; Kxe ([893 894],:) = [];
% Kxx * Phi = b % k reseni neznamych potencialu je treba znat vektor b % Kxe * E = b % Vektor E je dan potencialem 5. a 8. uzlu ktere je Phi5 = 5 V, Phi 8 = -5 E (1,1) = 5; E (2,1) = - 5; b = - Kxe*E;
% vypocitane hodnoty byly prevedeny na druhou strany -> znamenko -
x = Kxx\b; % ted je potreba doplnit spravne v seznamu potencialu Phi 5 a 8 uzlu. Phi (1:4,:) = x(1:4,:); Phi (5,:) = E(1,:); Phi (6:7,:) = x(5:6,:); Phi (8,:) = E(2,:); Phi (9:894,:) = x(7:892,:);
47
