Univerzita Karlova v Praze Matematicko-fyzikální fakulta DIPLOMOVÁ PRÁCE. Studium interakce plazma-pevná látka při středních tlacích

Univerzita Karlova v Praze Matematicko-fyzikální fakulta

DIPLOMOVÁ PRÁCE Studium interakce plazma-pevná látka při středních tlacích Štěpán Roučka Katedra fyziky povrchů a plazmatu

Vedoucí diplomové práce: Prof. RNDr. Rudolf Hrach, DrSc. Studijní program: Fyzika Studijní obor: Fyzika povrchů a ionizovaných prostředí

2008

Děkuji především svému vedoucímu Prof. RNDr. Rudolfu Hrachovi DrSc. za cenné informace poskytnuté při konzultacích. Také děkuji ostatním pracovníkům KFPP, kteří svými připomínkami přispěli ke zvýšení kvality této práce. Za nezbytnou podporu při studiu děkuji své rodině. V neposlední řadě patří můj dík mému počítači tetě Bertě, který mi tak věrně sloužil.

Tímto prohlašuji, že jsem svou diplomovou práci napsal samostatně a výhradně s použitím citovaných pramenů. Souhlasím se zapůjčováním práce. V Praze dne

Štěpán Roučka

i

Obsah 1 Úvod

1

2 Současné poznatky fyziky plazmatu 2.1 Plazma . . . . . . . . . . . . . . . . . 2.1.1 Nízkoteplotní plazma . . . . . 2.1.2 Argonové a kyslíkové plazma . 2.1.3 Sondová diagnostika plazmatu

. . . .

2 2 4 4 7

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

10 10 11 12 13 15 17 18 19 21 23 23 26 28 29 30 32

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

3 Modelování ve fyzice plazmatu 3.1 Boltzmannova kinetická rovnice . . . . . . . . . . . . . . . 3.1.1 Srážkový člen . . . . . . . . . . . . . . . . . . . . . 3.2 Rozvoj rozdělovací funkce do řady . . . . . . . . . . . . . . 3.3 Tekutinové modely . . . . . . . . . . . . . . . . . . . . . . 3.4 Chemická kinetika . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Řešení rovnic chemické kinetiky . . . . . . . . . . . 3.4.2 Numerické řešení obyčejných diferenciálních rovnic 3.4.3 Rozšíření modelu o difúzi . . . . . . . . . . . . . . . 3.5 Částicové modelování . . . . . . . . . . . . . . . . . . . . . 3.6 Particle-in-cell . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.1 Vzorkování a interpolace . . . . . . . . . . . . . . . 3.6.2 Řešení Poissonovy rovnice . . . . . . . . . . . . . . 3.6.3 Dynamika částic . . . . . . . . . . . . . . . . . . . 3.6.4 Srážky . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.5 Generování pseudonáhodných čísel . . . . . . . . . 3.7 Hybridní modely . . . . . . . . . . . . . . . . . . . . . . .

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

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

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

4 Cíle práce

33

5 Kinetický model kyslíkového plazmatu 5.1 Výpočet elektronové rozdělovací funkce . . . . . . . . . . . . . . . 5.2 Kinetický model . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34 34 35

ii

6 Částicový model pro vyšší tlaky 6.1 Model založený na generování náhodné volné dráhy 6.2 Obecný model srážek . . . . . . . . . . . . . . . . . 6.2.1 Referenční implementace . . . . . . . . . . . 6.2.2 Alternativní implementace I . . . . . . . . . 6.2.3 Alternativní implementace II . . . . . . . . 6.3 Vliv diskrétního času na srážkové procesy . . . . . . 6.3.1 Střední hodnoty veličin v modelu PIC . . . 6.4 Model srážek nezávislý na volbě časového kroku . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

39 39 41 44 46 48 49 51 54

7 Výsledky modelu 7.1 Verifikace modelu MCC . . . . . . . 7.2 Srovnání neselfkonzistentních modelů 7.2.1 Tlak 266 Pa . . . . . . . . . . 7.2.2 Tlak 1330 Pa . . . . . . . . . 7.3 Srovnání selfkonzistentních modelů . 7.3.1 Tlak 266 Pa . . . . . . . . . . 7.3.2 Tlak 665 Pa . . . . . . . . . . 7.4 Sondové charakteristiky . . . . . . . 7.5 Ostatní výsledky . . . . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

57 59 61 61 63 65 65 66 70 72

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

8 Závěr

77

A Popis programu

79

iii

Název práce: Studium interakce plazma-pevná látka při středních tlacích Autor: Štěpán Roučka Katedra (ústav): Katedra fyziky povrchů a plazmatu Vedoucí bakalářské práce: Prof. RNDr. Rudolf Hrach, DrSc. e-mail vedoucího: [email protected] Abstrakt: Předmětem práce je částicové modelování plazmatu při středních tlacích. Hlavní výsledky práce se týkají modelování srážek v algoritmu typu particlein-cell. V práci jsou diskutovány různé metody simulace srážek a je diskutována jejich přesnost a rychlost výpočtu. Dále byla navržena nová modifikace algoritmu particle-in-cell pro použití při středních tlacích. Numerické testy uvedené v závěrečné kapitole ukazují výhody nového algoritmu. Klíčová slova: plazma, modelování, particle-in-cell, srážky, Monte Carlo Title: Study of plasma-solid interaction at medium pressures Author: Štěpán Roučka Department: Department of Surface and Plasma Science Supervisor: Prof. RNDr. Rudolf Hrach, DrSc. Supervisor’s e-mail address: [email protected] Abstract: The thesis is aimed at particle simulation of plasma at medium pressures. Main results concern modelling of collisions in particle-in-cell algorithm. Different methods of modelling of collisions are discussed and their accuracy and efficiency are investigated. Furthermore a new modification of particle-in-cell algorithm for modelling at medium pressures was proposed. Numerical benchmarks in the final chapter show advantages of the new algorithm. Keywords: plasma, modelling, particle-in-cell, collisions, Monte Carlo

iv

Kapitola 1 Úvod Diplomová práce je zaměřena na rozšiřování možností počítačových modelů při středních tlacích. Nízkoteplotní plazma při středních tlacích nachází v současné době řadu uplatnění. Na jedné straně jde o technologické využití například při povrchové úpravě materiálů, leptání, čištění nebo nanášení tenkých vrstev, využívá se v zobrazovacích technologiích, v plynových laserech nebo při spalování odpadu. Na straně druhé je plazma využíváno v základním výzkumu, kde jeho prostřednictvím například studujeme elementární procesy a získáváme cenné poznatky z atomové fyziky, fyziky povrchů a podobně. Bohužel, přesná sondová diagnostika plazmatu při středních tlacích je velmi obtížná, nebot’ teoretický popis sondy v plazmatu při středním tlaku je taktéž obtížný až nemožný. Možnou cestou k pochopení procesů na sondě a interpretaci měření jsou právě počítačové modely. Avšak i pro počítačové modely tvoří oblast středních tlaků významný problém. V takovém plazmatu totiž dochází k ohromnému množství srážek, a proto se obvykle využívá aproximace tekutinového modelu. Vypovídací hodnota i přesnost tekutinových modelů jsou však omezené. Nejpřesnější třída počítačových modelů, což jsou částicové modely, musí zase důsledky každé srážky jednotlivě započítávat a to značně zpomaluje jejich výpočet. Vytvoření kvalitních modelů přitom otevírá cestu ke zpřesnění diagnostiky, vytváření nových metod a také k celkovému rozvoji fyziky plazmatu při středních tlacích. Proto je potřeba hledat vylepšení stávajících algoritmů vedoucí ke zpřesnění a především urychlení výpočtů, což je zaměření této práce.

1

Kapitola 2 Současné poznatky fyziky plazmatu 2.1

Plazma

Pojem plazma pochází pravděpodobně z řeckého slova πλασµα, které označuje tvárnou nebo rosolovitou hmotu. Do moderní terminologie zavedl tento termín v devatenáctém století významný český vědec Jan Evangelista Purkyně. Označil tak tekutou složku krve. Tato práce se však namísto krevní plazmy bude zabývat plazmatem, kterému dal jméno až Irving Langmuir (1928). Podle Langmuira je plazma taková část objemu ionizovaného plynu, kde je vyrovnaný náboj elektronů a iontů. Moderní definice plazmatu se sice různí, avšak pro potřeby této práce si vystačíme například s definicí dle (Chen, 1984, s. 19). „Plazma je kvazineutrální plyn nabitých a neutrálních částic vykazující kolektivní chování“ Smysl této definice nyní podrobně vysvětlíme a také uvedeme některá další kritéria vymezující pojem plazmatu. Kolektivní chování je způsobeno přítomností nabitých částic v plazmatu. Pohybem a změnou koncentrace nabitých částic totiž vzniká elektromagnetické pole, jehož dosah je mnohem větší než v případě například Van der Waalsovy interakce, kterou interagují neutrální částice. Lokální změna v rozložení částic tedy může způsobit odezvu i ve vzdálených oblastech plazmatu. Kvazineutralita znamená, že množství volného kladného náboje je přibližně rovné množství volného záporného náboje. V případě elektropozitivního plazmatu, neboli plazmatu obsahujícího pouze kladné ionty, s jedním druhem iontů tedy platí ne ≈ ni , kde ne je koncentrace elektronů a ni je koncentrace iontů. Tento požadavek je v souladu s původní Langmuirovou definicí. Splnění kvazineutrality je v plazmatu zajištěno mechanismem stínění, který nyní vysvětlíme. Mechanismus stínění se uplatňuje, pokud do plazmatu vložíme elektrický náboj nebo vodič na určitém potenciálu. Do okolí vloženého náboje jsou totiž přita2

3

2.1. PLAZMA

hovány částice opačné polarity, zatímco částice stejné polarity jsou naopak odpuzovány. Tento proces způsobí odstínění elektrického pole vloženého náboje. Pokud tedy dojde v plazmatu z nějakého důvodu k narušení kvazineutrality, vzniklé elektrické pole přitáhne částice opačné polarity a kvazineutralita se opět obnoví. Podle Debyeovy a Hückelovy teorie (Debye – Hückel, 1923) klesá elektrický potenciál náboje v plazmatu exponenciálně, přičemž charakteristickou vzdálenost poklesu nazýváme Debyeova vzdálenost s ε0 kB , (2.1) λD = P 2 i qi ni /Ti

kde suma probíhá přes všechny nabité složky plazmatu, qi označuje náboj, ni označuje koncentraci a Ti označuje teplotu dané složky. Teplota v tomto vztahu může být různá pro různé složky plazmatu. Plazma totiž často není v termodynamické rovnováze, a tak pojem teploty v obvyklém smyslu nelze definovat. Termodynamickou definici teploty proto nahrazujeme vztahem 3 hEi = kB T , 2

(2.2)

kde hEi je střední energie částic. Tento vztah, jenž je pro plyny v termodynamické rovnováze splněn identicky, se zde stává definicí. V případě, že takto definovaná teplota je pro různé složky plazmatu různá, říkáme, že je plazma neizotermické. Vrat’me se nyní k problému stínění. Pokud uvážíme, že vložený náboj se pohybuje, případně je ponechán v plazmatu pouze po krátkou dobu, může se stát, že ionty, které mají vyšší setrvačnost a nižší termální rychlosti oproti elektronům, nestihnou zaujmout rovnovážnou pozici. Potom se stínění účastní převážně elektrony a vztah pro Debyeovu délku nabývá běžně používaného tvaru s ε0 kB Te λD = . (2.3) qe2 ne Kritéria existence plazmatu Jelikož mechanismus stínění funguje na vzdálenostech řádově λD , je zřejmé, že pro zachování kvazineutrality musí být lineární rozměry plazmatu L podstatně větší, neboli L ≫ λD .

(2.4)

Proces stínění dále pro svoji funkci vyžaduje, aby ve stínící vrstvě byl dostatek nabitých částic. Označíme-li ND počet nabitých částic ve sféře o poloměru λD , potom musí platit ND ≫ 1 . (2.5)

Poslední kritérium zajišt’uje, aby vliv elektromagnetické interakce převažoval nad vlivem srážek s neutrálními částicemi. Dle (Chen, 1984) lze toto kritérium formulovat jako ωτ > 1 , (2.6)

2.1. PLAZMA

4

kde ω je frekvence typických oscilací plazmatu a τ je střední doba mezi srážkami s neutrálními částicemi. Vliv zvyšování tlaku Pod pojmem zvyšování tlaku rozumíme v této práci přechod k tlakům řádu 100 Pa až 1000 Pa. Důsledkem zvyšování tlaku je několik jevů. I když stupeň ionizace může klesat, objemová koncentrace nabitých částic s tlakem obvykle roste, protože roste celková koncentrace částic. Dále s rostoucím tlakem dochází ke změně v poměru koncentrací jednotlivých složek a s tím souvisí i změna v zastoupení jednotlivých srážkových procesů. Z hlediska této práce však nejdůležitější změna spočívá ve značném nárůstu počtu srážek. Dříve uvedené změny totiž spočívají pouze ve změně některých parametrů modelu, zatímco rostoucí počet srážek představuje významný problém z hlediska výpočetní náročnosti modelu. Rostoucí počet srážek také může zapříčinit porušení vztahu (2.6), a tak může plazma při středních tlacích do jisté míry vybočovat z kritérií dle (Chen, 1984).

2.1.1

Nízkoteplotní plazma

V této práci se budeme zabývat výhradně nízkoteplotním plazmatem, proto zde tento typ plazmatu krátce představíme. Nízkoteplotní plazma je pojem označující neizotermické plazma vznikající mimo jiné v takzvaném doutnavém výboji. K jeho základním vlastnostem patří nízká teplota iontů. Ta je blízká teplotě neutrálních částic respektive pokojové teplotě a obvykle nepřesahuje stovky kelvinů. Na druhou stranu teplota elektronů, které zodpovídají za ionizaci částic plynu, dosahuje typicky elektronvoltů, neboli desítek tisíc kelvinů. Je tedy zřejmé, že toto plazma je daleko od termodynamické rovnováhy, avšak může se nalézat ve stacionárním stavu, což znamená, že makroskopické parametry plazmatu jsou v čase konstantní. Jak již bylo uvedeno, nízkoteplotní plazma můžeme generovat v doutnavém výboji. To je výboj v plynu obvykle za sníženého tlaku. Mechanismem udržení doutnavého výboje je potenciálová a kinetická emise elektronů z katody a následná lavinová ionizace v 2. katodovém temném prostoru (Martišovitš, 2004). Schéma struktury doutnavého výboje je zpodobněno na obrázku 2.1. Z obrázku je zřejmé, že doutnavý výboj vykazuje poměrně komplikovanou strukturu. Oblastí našeho zájmu je však pouze kladný sloupec, kde je homogenní a relativně slabé elektrické pole (Sugawara – Stansfield, 1998). Plazma zde obsažené odpovídá definici nízkoteplotního plazmatu. Navíc v případě, že výboj není narušen nestabilitami, je plazma v této oblasti homogenní v podélném směru.

2.1.2

Argonové a kyslíkové plazma

Nízkoteplotní plazma má obecně mnoho technologických aplikací a kyslíkové plazma díky své vysoké reaktivitě patří k často používanému typu v oblastech, kde je potřeba chemicky aktivní plazma. Argon, jakožto inertní plyn, netvoří natolik

5

2.1. PLAZMA

Temný 1. katodový prostor: 2. katodový Faradayův

Anodový

+

− Katodové světlo Negativní doutnavé světlo

Kladný sloupec

Anodové světlo

Obrázek 2.1: Schéma doutnavého výboje, převzato z (Roučka, 2006).

reaktivní plazma, avšak jeho aplikace jsou také široké. Jmenujme nyní několik zajímavých technologických využití kyslíkového a argonového plazmatu. Jde především o plazmové leptání, povrchovou úpravu materiálů, depozici tenkých vrstev, použití ve směsích plynových laserů, ale také například destrukce těkavých organických sloučenin ve výboji při atmosférickém tlaku (Grossmannová et al., 2007), čištění vzorků pro SEM (Mueller et al., 2007) nebo moderní metody desinfekce a destrukce nežádoucích tkání v medicíně (Laroussi – Lu, 2005). Jak již bylo uvedeno, argon je inertní plyn, a díky tomu jsou procesy probíhající v argonovém plazmatu relativně jednoduché. Nízkoteplotní argonové plazma je typicky tvořeno převážně neutrálními atomy a v závislosti na podmínkách výboje může být významná část argonových atomů v excitovaných metastabilních stavech, především jako Ar(3 P2 ). Nabité částice jsou zastoupeny prakticky výhradně elektrony a jednonásobnými kladnými ionty Ar+ . Vícenásobná ionizace je v nízkoteplotním plazmatu nepravděpodobná. Elektronová konfigurace atomu v základním stavu je tvaru [Ne]3s2 3p6 a můžeme ji popsat pomocí termu 1 S0 . Dále jsou v nízkoteplotním plazmatu důležité čtyři stavy s elektronovou konfigurací 3p5 4s. Dva z těchto stavů, 3 P0 a 3 P2 , jsou metastabilní s dlouhou dobou života, zatímco zbylé dva, 3 P1 a 1 P1 , jsou rezonanční, tedy s krátkou dobou života. Rozložení energetických hladin argonu je schematicky znázorněno v obrázku 2.2. Nejnižší excitační energii má metastabilní stav 3 P2 (1s5 v Paschenově notaci) ležící 11,55 eV nad základním stavem. První ionizační energie argonu je 15,76 eV. Pro dvojnásobnou a trojnásobnou ionizaci jsou již potřeba energie 27,7 eV a 40,8 eV, pročež jsou tyto procesy značně nepravděpodobné. Metastabilní stavy vznikají excitací elektronem a zanikají obvykle další excitací, deexcitací na stěně nebo mechanismem Penningovy ionizace. Ionty vznikají běžně ionizací atomu v základním nebo metastabilním stavu elektronem anebo Penningovou ionizací. Účinné průřezy některých interakcí elektronů s argonem jsou vyneseny v grafu 2.3. Pochopení procesů v kyslíkovém plazmatu je podstatně složitější než v případě argonu. Kyslík je totiž molekulární plyn, kde dochází k excitaci různých

6

2.1. PLAZMA

Obrázek 2.2: Důležité energetické hladiny argonového atomu, převzato z (Brok, 1998).

vibračních a rotačních hladin při relativně nízkých energiích, dále také elektronicky excitované stavy kyslíkové molekuly mají nižší prahové energie a podstatně složitější strukturu. Navíc je kyslíkové plazma elektronegativní, což znamená, že se v něm vyskytují i záporné ionty. Nejdůležitější složky kyslíkového plazmatu shrnuje následující seznam: • e− elektrony • O2 molekula kyslíku v základním stavu • O atom kyslíku v základním stavu, silně reaktivní, reasociuje na stěnách • O3 ozón − • O− , O− 2 , O3 záporné ionty, vznikají díky elektronegativitě kyslíku

• O+ , O+ 2 kladné ionty

7

• O2 (v) vibračně excitované stavy, účinné průřezy pro excitaci čtyřech nejnižších vibračních stavů s energiemi 0,19, 0,38, 0,57 a 0,75 eV jsou uvedeny v obrázku 2.4 • O2 (a1 ∆g ), O2 (b1 Σ+ g ) singletní elektronicky excitované stavy s energiemi 0,977 eV a 1,627 eV, účinné průřezy jejich excitace jsou uvedeny v obrázku 2.5 • O(1 D) metastabilní excitovaný stav atomárního kyslíku s energií 1,97 eV

7

2.1. PLAZMA

σ [10−20 m2 ]

102

pružný rozptyl excitace 3 P2 ionizace

101 100 10−1 10−2 10−3

10−2

10−1

100

101

102

103

104

E [eV] Obrázek 2.3: Účinné průřezy interakcí elektronů s atomy argonu, převzato z (Phelps).

• O∗2 mnoho elektronicky excitovaných stavů, například skupiny stavů na energiích přibližně 4,5, 6,0, 8,4 a 9,97 eV a stavy excitované elektronicky i vibračně Všechny uvedené složky spolu navzájem interagují velkým množstvím reakcí. V nízkoteplotním plazmatu hrají důležitou roli především interakce s elektrony, nebot’ pouze elektrony mají dostatek energie k tomu, aby vybudily elektronicky excitované stavy nebo způsobily ionizaci. Účinné průřezy pro excitaci vibračních stavů kyslíkové molekuly jsou znázorněny v grafu 2.4. Ostatní účinné průřezy zahrnující rotační a elektronickou excitaci a ionizaci jsou vyneseny v grafu 2.5. Další interakce v kyslíkovém plazmatu budou podrobněji probrány v kapitole 5 v souvislosti s modelem chemické kinetiky kyslíkového plazmatu.

2.1.3

Sondová diagnostika plazmatu

Sondová diagnostika plazmatu je jednou z technicky nejjednodušších a zároveň nejrozšířenějších metod diagnostiky plazmatu. Princip této metody spočívá v zasunutí jedné nebo více vodivých sond do objemu plazmatu a následném měření elektrického proudu tekoucího na sondu, případně na různé části sondy, v závislosti na napětí přivedeném na sondu. Vodivou sondu nazýváme také Langmuirova sonda. Vložené sondy mohou mít různé tvary, což potom ovlivňuje interpretaci výsledků. Asi nejčastěji se využívá válcová sonda, neboli drát. Schéma experimentálního zapojení pro měření metodou dvou sond je zakresleno v obrázku 2.6. Měření metodou jedné sondy se liší tím, že jednu ze sond nahrazujeme jinou elektrodou, obvykle anodou. Ačkoliv je metodika měření poměrně jednoduchá, následná analýza výsledků může být poměrně komplikovaná. Existují různé teorie popisující měření pro různé

8

2.1. PLAZMA

σ [10−20 m2 ]

101

v v v v

100

=1 =2 =3 =4

10−1 10−2 10−3

100

101 E [eV]

Obrázek 2.4: Účinné průřezy vibrační excitace molekul O2 , převzato z (Phelps).

101

pružná rotační O2 (a1 ∆g ) O2 (b1 Σ+ g) O2 (4,5) O2 (6,0) O2 (8,4) O2 (9,97) − O+ 2 +e 130 nm disociace

σ [10−20 m2 ]

100

10−1

10−2

10−3 10−1

100

101

102

103

104

E [eV] Obrázek 2.5: Účinné průřezy interakcí elektronů s molekulami O2 , převzato z (Phelps).

2.1. PLAZMA

9

Obrázek 2.6: Schéma zapojení pro měření metodou dvou sond.

parametry měření. O použité teorii rozhodují poměry Debyeovy délky λD , poloměru sondy rp a střední volné dráhy λ. Vezměme si například sondu o poloměru rp = 1 · 10−4 m v plazmatu o tlaku 665 Pa a koncentraci iontů 1 · 10−15 m−3 , kterou budeme modelovat v některých našich výpočtech. Pro tu platí platí přibližně λD = 4 · 10−4 m a λ = 2,5 · 10−5 m. Protože volná dráha je podstatně menší než Debyeova délka, říkáme, že stínící vrstva je srážková. Teoretická analýza měření se srážkovou stínící vrstvou je komplikovaná a výsledky jsou vždy pouze přibližné, proto zde může modelování sloužit jako vhodný doplněk experimentálních technik. Cílem práce není analýza sondových měření, a tak zde nebudeme rozebírat příslušné teorie, prezentovaný model však může k tomuto účelu sloužit.

Kapitola 3 Modelování ve fyzice plazmatu Pod pojmem modelování rozumíme jednak hledání rovnic popisujících zkoumaný fyzikální systém, jednak jejich následné řešení. Jak brzy uvidíme, rovnice popisující plazma jsou vesměs velmi komplikované a jejich analytické řešení není až na výjimky možné. Proto se obvykle uchylujeme k numerickému řešení těchto rovnic. V případě, že rovnice řešíme s použitím počítače, hovoříme o počítačovém modelu. Ve výzkumu plazmatu hrají počítačové modely významnou roli, nebot’ umožňují pochopení procesů probíhajících v experimentech i návrh experimentů nových. Nyní uvedeme Boltzmannovu rovnici, která může sloužit jako základ většiny modelů plazmatu.

3.1

Boltzmannova kinetická rovnice

Plazma je systém mnoha částic, nacházející se často daleko od termodynamické rovnováhy, a proto k jeho popisu využíváme metod nerovnovážné statistické fyziky. Stav souboru identických částic je určen rozdělovací funkcí f (x, v, t)dxdv ,

(3.1)

která při vhodném normování udává pravděpodobnost výskytu částice v elementu fázového prostoru (v1 , v1 + dv1 ) × (v2 , v2 + dv2 ) × (v3 , v3 + dv3 ) × (x1 , x1 + dx1 ) × (x2 , x2 + dx2 ) × (x3 , x3 + dx3 ) . (3.2) Pod symbolem dx zde rozumíme součin dx1 dx2 dx3 a obdobně pro dv. V případě, že systém sestává z více složek, je každá popsána samostatnou rozdělovací funkcí fi . Rozdělovací funkce je obvykle normovaná na koncentraci částic, neboli Z f (x, v, t)dv = n(x, t) . (3.3)

Hodnota f (x, v, t)dxdv při tomto normování vyjadřuje střední počet částic v objemu dv dx. Nyní pouze nastíníme odvození Boltzmannovy rovnice, podrobnější 10

3.1. BOLTZMANNOVA KINETICKÁ ROVNICE

11

rozbor této problematiky může čtenář najít například v (Kracík et al., 1974, kap. 3). V nepřítomnosti srážek se částice pohybují spojitě fázovým prostorem, a proto musí hustota pravděpodobnosti splňovat rovnici kontinuity 3 dxi dvi ∂ ∂f X ∂ + f + f = 0. ∂t ∂x dt ∂v dt i i i=1

(3.4)

Za totální derivace v předchozím vztahu můžeme dosadit z Newtonových pohybových rovnic a přepsat do vektorového tvaru ∂f F + ∇x · (vf ) + ∇v · f = 0. (3.5) ∂t m V případě síly nezávislé na rychlosti nebo speciálně v případě Lorentzovy síly můžeme ještě rovnici přepsat do tvaru ∂f F + v · ∇xf + · ∇v f = 0 . ∂t m

(3.6)

Rovnice (3.6) je Vlasovova rovnice popisující vývoj rozdělovací funkce bezesrážkového systému.

3.1.1

Srážkový člen

Srážky částic způsobují skokovou změnu polohy částice ve fázovém prostoru, nebo dokonce vznik a zánik částic. Z toho důvodu přidáváme do rovnice kontinuity respektive Vlasovovy rovnice srážkový člen vyjadřující bilanci částic v objemovém elementu vlivem srážek X δfi F ∂fi . (3.7) + v · ∇xfi + · ∇v fi = ∂t m δt j j Srážkový člen (δfi /δt)j na pravé straně předchozí rovnice symbolizuje vliv srážek částic typu i s částicemi typu j a jeho význam nyní objasníme. Pravděpodobnost, že se částice typu i o rychlosti v i srazí v okamžiku dt se svazkem částic typu j o rychlosti v j a vyletí do prostorového úhlu ω je dána diferenciálním účinným průřezem σij (ω, gij ) podle vztahu σij (ω, gij )gij nj dω dt ,

(3.8)

kde g ij = v j − v i , a nj je koncentrace částic typu j, element prostorového úhlu sin θ dθ dφ jsme označili dω. V případě, že částice typu j podléhají obecnému rozdělení, nahradíme nj hustotou pravděpodobnosti fj dv j a pravděpodobnost vypočteme integrací přes rychlosti Z dω dt dv j σij (ω, gij )gij fj (v j ) . (3.9) vj

3.2. ROZVOJ ROZDĚLOVACÍ FUNKCE DO ŘADY

12

Nás však zajímá celkový počet částic typu i opouštějících element fázového prostoru, a proto předchozí vztah vynásobíme množstvím částic typu i v elementu a přeintegrujeme přes všechny možné úhly rozptylu Z Z − Nxvt dxdv i dt = fi (v i )dxdv i dt dω dv j σij (ω, gij )gij fj (v j ) . (3.10) ω

vj

V případě pružných srážek lze pro počet částic přicházejících do elementu fázového prostoru vlivem srážek ze symetrie zákonů zachování odvodit analogický vztah tvaru Z Z + Nxvt dxdv i dt = dxdv i dt dω dv j σij (ω, gij )gij fi (v ′i )fj (v ′j ) , (3.11) ω

vj

kde v ′i a v ′j jsou rychlosti před srážkou takové, aby výsledná rychlost po rozptylu pod úhlem ω byla rovna v i . Výrazy (3.11) a (3.10) můžeme sjednotit do jednoho integrálu a po vydělení elementárním objemem dv i dx a časovým elementem dt dostaneme konečně srážkový člen pro pružné srážky mezi částicemi typu i a j elastic Z Z δfi = dω dv j (fi (v ′i )fj (v ′j ) − fi (v i )fj (v j ))σij (ω, gij )gij . (3.12) δt j ω vj Tento srážkový člen platí za podmínek molekulárního chaosu, což znamená, že interagující částice musí být nekorelované. Pro popis obecných nepružných srážek potom existují různé modifikace srážkového členu určené speciálně pro daný typ interakce. Veškeré modely plazmatu můžeme nyní interpretovat jako různé formy řešení Boltzmannovy rovnice. Exaktní řešení této rovnice je s výjimkou několika elementárních případů nemožné, a proto hledáme různé přibližné metody řešení. Tyto metody lze rozdělit do několika skupin, které si nyní stručně popíšeme.

3.2

Rozvoj rozdělovací funkce do řady

Postup, který je nejbližší přesnému řešení Boltzmannovy rovnice, avšak má také velmi přísné předpoklady, spočívá v rozvoji rozdělovací funkce do kulových funkcí v rychlostním prostoru. Nejčastěji používaným kompromisem mezi složitostí výpočtu a přesností výsledku je použití prvních dvou členů rozvoje, jde o takzvanou Lorentzovu aproximaci. Navíc obvykle předpokládáme prostorově homogenní problém a rozdělovací funkce tedy nabývá tvaru (Holstein, 1946) f (v) = f0 (v) +

v · f 1 (v) . v

(3.13)

S užitím těchto předpokladů lze Boltzmannovu rovnici převést na jednodušší integrodiferenciální rovnici pro izotropní část rozdělovací funkce a tu je možné numericky řešit (Frost – Phelps, 1962; Morgan – Penetrante, 1990).

3.3. TEKUTINOVÉ MODELY

13

Pro výpočet rozdělovací funkce v této práci užíváme program ELENDIF (Morgan – Penetrante, 1990). Ze znalosti rozdělovací funkce potom můžeme určit některé makroskopické veličiny, jako je driftová rychlost v závislosti na magnetickém poli nebo pohyblivost. Typickou aplikací této metody potom je korekce účinných průřezů na základě srovnání s experimentálními hodnotami makroskopických veličin. Tímto způsobem byla vypočtena také sada účinných průřezů (Phelps), kterou budeme v této práci používat pro výpočet elektronové rozdělovací funkce. Další důležitou aplikací je výpočet rychlostních konstant reakcí, jež jsou potom použity jako vstup do modelu chemické kinetiky. Ten bude popsán v dalším odstavci.

3.3

Tekutinové modely

V případě, že předmětem zájmu není samotná rozdělovací funkce, nýbrž prostorové rozdělení koncentrace a toku složek plazmatu, je výhodné použít takzvané tekutinové modely. Jak již název napovídá, tyto modely popisují plazma jako směs tekutin řízenou rovnicemi pro dynamiku kontinua. Tekutinové modely nejsou předmětem této práce, a proto uvedu pouze velmi stručný popis. Výchozím bodem pro sestavení evolučních rovnic tekutinového modelu je opět Boltzmannova kinetická rovnice. Abychom řešení zjednodušili, násobíme Boltzmannovu rovnici n-tou mocninou hybnosti, a poté středujeme přes prostor hybnosti. Tak dostaneme n-tý moment Boltzmannovy rovnice. Nultý moment potom popisuje zachování hmoty, první moment zachování hybnosti, druhý moment zachování energie atd. Nejjednodušší z této třídy modelů jsou magnetohydrodynamické modely, které nahlížejí na plazma jako na jednosložkovou vodivou tekutinu. Tok plazmatu je určen prvním momentem Boltzmannovy rovnice respektive Navierovými-Stokesovými rovnicemi například ve tvaru ∂u nm + (u · ∇)u = ∇p + nmg + j × B , (3.14) ∂t kde g označuje intenzitu gravitace a Lorentzovu sílu reprezentuje pouze člen od magnetického pole, což je případ nerelativistické magnetohydrodynamiky. Tlak považujeme za skalární veličinu a tím zanedbáváme viskozitu plazmatu. Elektrický proud tekoucí plazmatem je dán zobecněným nerelativistickým Ohmovým zákonem j = σ(E + u × B) . (3.15)

Dále se k popisu přidává rovnice kontinuity respektive nultý moment Boltzmannovy rovnice ∂n + ∇ · (nu) = 0 (3.16) ∂t

14

3.3. TEKUTINOVÉ MODELY

a typicky může do modelu vstoupit také rovnice pro energii a některé vztahy z termodynamiky. Popis magnetického pole vychází z Maxwellových rovnic, a při uvážení nerelativistické limity a vysoké vodivosti dostáváme vztah ∂B ∆B = ∇ × (u × B) + , ∂t µ0 σ

(3.17)

který popisuje časový vývoj magnetického pole. Magnetohydrodynamický popis je vhodný především pro modelování silně ionizovaného, vysoce vodivého plazmatu a jeho častou aplikací je plazma meziplanetárního prostoru, jež tyto vlastnosti splňuje. Obecnější tekutinové modely popisují plazma jako směs tekutin, které odpovídají jednotlivým složkám plazmatu. Stejně jako v případě magnetohydrodynamického modelu je popis založen na rovnici kontinuity a Navier-Stokesových rovnicích. Tentokrát však tyto rovnice píšeme zvlášt’ pro elektrony, různé druhy iontů, případně i pro neutrální částice. Na rozdíl od magnetohydrodynamického modelu vystupuje v rovnici kontinuity člen vzniku a zániku částic δni /δt, což vyjadřuje, že částice dané složky se nezachovávají, nýbrž mohou vznikat a zanikat například při ionizaci a rekombinaci δni ∂ni + ∇ · (ni ui ) = . ∂t δt

(3.18)

Pohybová rovnice plazmatu, neboli první moment Boltzmannovy rovnice, je vyjádřena vztahem ∂ui + (ui · ∇)ui = ni qi (E + ui × B) + nimi g − ∇pi + mi ni νi ui , (3.19) mi ni ∂t kde νi označuje srážkovou frekvenci a člen mi ni νi ui je přibližným vyjádřením ztráty hybnosti vlivem srážek. Dále bychom mohli dostat rovnici pro energetickou bilanci výpočtem druhého momentu Boltzmannovy rovnice. Vlastností tohoto řetězce rovnic však je, že počet neznámých je větší než počet rovnic. Například rovnice kontinuity nám popisuje koncentraci částic, avšak k jejímu řešení potřebujeme znát rozložení rychlosti. To bychom mohli určit z pohybové rovnice, avšak k jejímu řešení bychom potřebovali znát tlak a tak dále. K ukončení tohoto řetězce rovnic je tedy potřeba použít nějaký další vztah a pro tento účel používáme poznatky termodynamiky. Konkrétně tlak určíme ze vztahu p = nkB T .

(3.20)

Hlavní výhodou spojitého přístupu je velmi vysoká rychlost výpočtu ve srovnání s částicovými modely a tím pádem také schopnost modelovat poměrně komplikované systémy. Nevýhody však plynou ze zanedbání podstatné části informace oproti částicovému modelu. Středováním Boltzmannovy rovnice přes rychlosti totiž ztratíme tvar rychlostní rozdělovací funkce. Tvar rozdělovací funkce je

3.4. CHEMICKÁ KINETIKA

15

sice do modelu zanesen zvenčí na základě určitých předpokladů o stavu systému. Často je však sama rozdělovací funkce předmětem studia a v tom případě nám spojitý model nemůže sloužit. Další problém je, že spojité modely nedostatečně popisují trajektorie částic. Jevy které jsou způsobeny pohyby částic po dlouhých volných drahách proto nejsou ve spojitých modelech dobře reprezentovány.

3.4

Chemická kinetika

Při modelování plazmatu (obzvláště toho chemicky aktivního) nebo obecněji při zkoumání složitých chemických procesů se často setkáváme s následujícím problémem: máme směs mnoha složek, které navzájem interagují prostřednictvím ještě řádově většího množství chemických reakcí. Chemické reakce je sice teoreticky možné zahrnout do obecnějšího tekutinového modelu, avšak, pomineme-li náročnost tekutinového modelu o desítkách složek, často nám stačí znát pouze průběh koncentrací za předpokladu prostorově homogenního problému. V těchto případech tedy používáme metodu chemické kinetiky. Metoda chemické kinetiky, stejně jako všechny výše uvedené metody, vychází z Boltzmannovy rovnice, respektive z rovnice kontinuity ∂n = −∇ · (nu) + I − R , (3.21) ∂t kde I označuje vznik dané složky vlivem chemických procesů a naopak člen R označuje její zánik. Členy I a R jsou integrály srážkového členu Boltzmannovy rovnice. Pro začátek budeme předpokládat, že plazma je prostorově homogenní a zanedbáme vliv toku ∇ · (nu). Později budeme tento předpoklad podrobně diskutovat. Z parciální diferenciální rovnice tak získáváme obyčejnou diferenciální rovnici, v případě vícesložkového plazmatu pak soustavu obyčejných diferenciálních rovnic dni = Ii − Ri ; i = 1 . . . K , (3.22) dt kde K označuje počet složek. Jednotlivé členy vzniku a zániku jsou určeny probíhajícími procesy a jejich rychlostními konstantami. V obecném případě chemické kinetiky, kde jedna reakce může například popisovat komplexní proces složený z několika dílčích reakcí, mohou být tyto členy libovolnou funkcí koncentrací reaktantů a některých dalších proměnných, jako například koncentrace katalyzátoru. V našem modelu plazmatu si však vystačíme pouze s elementárními procesy, jejichž rychlost r lze vyjádřit integrací jednoduchého srážkového členu Boltzmannovy rovnice. Uvažme nyní dvoučásticovou interakci, jejíž pravděpodobnost popisuje diferenciální účinný průřez definovaný vztahem (3.8). Četnost probíhajících reakcí bude určena integrací výrazu (3.10) přes v i Z Z Z dtdx dv i dv j dωσij (ω, gij )gij fi (v i )fj (v j ) . (3.23)

16


Tento výraz nyní pro zjednodušení vydělíme dxdt, čímž dostaneme počet reakcí v jednotkovémR objemu za jednotku času. Dále definujeme totální účinný průřez jako σ(gij ) = σ(ω, gij )dω a dostáváme rychlost reakce Z Z r= dv i dv j σij (gij )gij fi (v i )fj (v j ) . (3.24) Z podmínek molekulárního chaosu vyplývá, že stavy interagujících částic jsou nezávislé, a proto součin fi (v i )fj (v j ) po normování koncentracemi udává hustotu pravděpodobnosti výskytu páru částic o rychlostech v i a v j , neboli P (vi , v j ) =

fi (v i )f (v j ) . ni nj

(3.25)

Vztah (3.24) můžeme tedy chápat jako výpočet střední hodnoty přes všechny páry částic a výsledkem je r = ni nj hσ(gij )gij i .

(3.26)

Definujeme nyní rychlostní konstantu reakce k ≡ hσ(gij )gij i .

(3.27)

Ukazuje se, že vztah (3.26) lze zobecnit i na případ tříčásticových a jednočásticových reakcí vztahem Y Y [Xj ] . (3.28) nj ≡ k r=k reaktanty

reaktanty

Zde jsme zavedli označení Xi ;

i = 1...K

(3.29)

pro jednotlivé složky plazmatu, kde symbol [Xj ] vyjadřuje koncentraci složky Xj . Uvedené zobecnění platí například pro tříčásticové reakce, kde nejprve v důsledku jedné srážky vznikne krátce žijící stav a ten potom vlivem další srážky přejde ve finální produkt. Příkladem jednočásticových interakcí mohou být spontánní přechody, zjednodušený popis vícečásticových interakcí nebo formální započtení vlivu difúze. Příspěvek jedné reakce ke členu vzniku nebo zániku lze tedy vypočíst jako součin koncentrací reaktantů a rychlostní konstanty dané reakce. V případě, že do reakce vstupuje více částic některého reaktantu, vstupují koncentrace tohoto reaktantu do výpočtu rychlosti reakce taktéž vícenásobně. Celý postup ještě objasníme na jednoduchém příkladu. Uvažme následující dvě reakce k0 O+ + e− − → O + hν , k

1 O+ + e− + e− − → O + e− .

17


Číslem nad šipkou zde označujeme, jak je obvyklé, rychlostní konstantu dané reakce. V uvedených reakcích figurují tři chemické složky a kvantum elektromagnetického záření, kterým se v této práci zabývat nebudeme. Proto dostaneme tři diferenciální rovnice pro jejich koncentrace d[e− ] = −k0 [O+ ][e− ] − 2k1 [O+ ][e− ][e− ] + k1 [O+ ][e− ][e− ] , dt d[O] = +k0 [O+ ][e− ] + k1 [O+ ][e− ][e− ] , dt d[O+ ] = −k0 [O+ ][e− ] − k1 [O+ ][e− ][e− ] . dt Tyto soustavy rovnic jsou obvykle natolik složité, že je dokážeme řešit pouze numericky, a i numerické řešení je relativně komplikované. Formulace konkrétního fyzikálního problému bude rozvedena v kapitole 5.

3.4.1

Řešení rovnic chemické kinetiky

Nyní stručně nastíníme metodu řešení rovnic chemické kinetiky. Pro popis složitějšího systému o několika složkách, ve kterém může probíhat mnoho chemických reakcí, je vhodné zavést určitý formalismus popisu takové soustavy rovnic pro potřeby jejich následného numerického řešení. Jak již bylo uvedeno, máme plazma skládající se ze složek Xi . Dále necht’ probíhá N chemických reakcí s rychlostními konstantami ks , s = 1 . . . N. Symbolem R(s) označíme množinu indexů reaktantů s-té reakce a symbol P (s) bude značit produkty. Příklad pro vysvětlení: Pokud bude s-tá reakce popsána rovnicí k

s X1 + X3 − → X5 + X6 ,

(3.30)

tak pro množiny R(s) a P (s) platí R(s) = {1, 3} , P (s) = {5, 6} .

(3.31)

V případě, že jsou některé reaktanty vícenásobné, budou v odpovídající násobnosti obsaženy také v množinách R a P . Obecnou chemickou rovnici potom můžeme formálně zapsat ve tvaru X X ks Xi − → Xi . (3.32) i∈R(s)

i∈P (s)

Diferenciální rovnice popisující vývoj koncentrace i-té složky dále nabývá tvaru X Y X Y dni = ks nj − ks nj . dt s,i∈P (s)

j∈R(s)

s,i∈R(s)

j∈R(s)

(3.33)

18


Při řešení diferenciální rovnice implicitními metodami potřebujeme ještě znát jakobián pravé strany rovnice (3.33), kterou označíme fi . X X ∂fi N(j, s) Y N(j, s) Y nk − ks nk . = ks ∂nj nj nj s,i∈P (s) j∈R(s)

k∈R(s)

s,i∈R(s) j∈R(s)

(3.34)

k∈R(s)

Symbol N(j, s) značí násobnost výskytu čísla j v množině reaktantů R(s).

3.4.2

Numerické řešení obyčejných diferenciálních rovnic

Soustava rovnic (3.33) nemá, s výjimkou některých triviálních případů, analytické řešení. Z toho důvodu musíme tuto soustavu řešit numericky. Ukazuje se, že chemické procesy popsané těmito rovnicemi probíhají na řádově odlišných časových škálách, což vyžaduje použití proměnného časového kroku. Navíc vykazuje tato soustava velikou tuhost, což znamená, že běžné numerické algoritmy mají problémy se stabilitou. Nejjednodušší numerickou metodou řešení diferenciálních rovnic je Eulerova metoda. Máme-li soustavu rovnic y ′ = f (y) ,

(3.35)

dostáváme její přibližné řešení ze vztahu y n+1 = y n + hf (y n ) ,

(3.36)

kde h je délka integračního kroku. Rovnice (3.36) je příkladem takzvaného explicitního diferenčního vzorce, což znamená, že hodnota y n+1 je ve vztahu explicitně vyjádřena. Eulerova metoda se sice pro svoji nízkou přesnost příliš často nepoužívá, avšak i pokročilejší a hojně využívaná metoda Rungeho-Kutty je založena explicitním diferenčním schématu. Běžně používané metody založené na explicitním diferencování mají však tu nevýhodu, že délka jejich časového kroku je z důvodu stability omezena nejrychleji probíhajícím dějem (Press et al., 2002, kap. 16-6). V případě soustav s velkou tuhostí je tento požadavek značně omezující. Řešením těchto problémů je přechod k implicitnímu diferenčnímu schématu. Například implicitní verze Eulerova vzorce je dána vztahem y n+1 = y n + hf (y n+1 ) . (3.37) Toto schéma je v případě lineárních diferenciálních rovnic absolutně stabilní, a i pro nelineární diferenciální rovnice dává dobré výsledky. Chceme-li však ze vztahu (3.37) určit novou hodnotu y n+1 , musíme obecně řešit soustavu nelineárních rovnic. Musíme však uvážit, že změna y je malá, a díky tomu můžeme funkci f rozvinout v řadu v bodě y n a linearizovat. Potom již snadno vypočteme, že −1 ∂f y n+1 = y n + h 1 − h f (y n ) , (3.38) ∂y

19


což je takzvaná semi-implicitní Eulerova metoda. Ze vztahu (3.38) je zřejmé, že k výpočtu nové hodnoty y n+1 potřebujeme znát nejenom pravou stranu rovnice, ale také její jakobián. Více viz (Press et al., 2002, kap. 16-6). Na podobném principu jako semi-implicitní Eulerova metoda je založena také Rosenbrockova metoda. V našem modelu využíváme právě tuto metodu s adaptivním určením integračního kroku implementovanou dle (Press et al., 2002, kap. 16-6).

3.4.3

Rozšíření modelu o difúzi

Nyní do modelu zpětně zavedeme vliv difúze, který jsme v předchozí diskusi zanedbali. Budeme uvažovat rovnici kontinuity respektive rovnici difúze bez členů vzniku a zániku. Toto zanedbání je zcela oprávněné v případě, že zdrojové členy jsou lineární funkcí koncentrace. V opačném případě bude výsledek pouze přibližný. Uvážíme-li navíc, že difúzní koeficient nebude závislý na poloze docházíme k difúzní rovnici ve tvaru ∂n = D∆n . (3.39) ∂t Tuto rovnici je možno řešit separací proměnných. Jak známo, obecné řešení prostorové části této rovnice lze popsat pomocí množiny vlastních funkcí, takzvaných difúzních modů, které splňují vztah nj ∆nj = − 2 , (3.40) Λj kde Λj je vlastní číslo, takzvaná difúzní délka, odpovídající modu nj . Snadno nahlédneme, že časová závislost takového řešení je dána exponenciální závislostí 2

n(x, t) = e−Dt/Λj nj (x) .

(3.41)

Jelikož funkce nj tvoří ortogonální systém, lze řešení s libovolnou počáteční podmínku rozepsat do řady z členů typu (3.41). Nyní musíme určit okrajovou podmínku odpovídající konkrétnímu fyzikálnímu problému a na ní vyřešit rovnici (3.39). Často používanou okrajovou podmínkou je nulová koncentrace částic na stěnách aparatury. Tato podmínka nemůže být nikdy přesně splněna, nebot’ nulová koncentrace by implikovala nulový tok na stěny, je však dobrou aproximací krátké volné dráhy částic a vysokého koeficientu absorpce stěn. Přesnější okrajovou podmínku získáme rozborem toku částic na stěny aparatury. Definujeme nyní usměrněný tok částic na stěnu Γw . Tato veličina popisuje částice, které se na stěně zachytávají nebo rekombinují. Dále definujeme tok částic na stěnu Γ+ a od stěny Γ− , pro které platí Γw = Γ+ − Γ− .

(3.42)

Usměrněný tok částic na stěnu můžeme určit z Fickova zákona Γw = −D∇n|w · k ,

(3.43)

20


kde k je normálový vektor ke stěně a |w označuje vyhodnocení veličiny na stěně. V případě, že rychlostní rozdělovací funkce částic bude blízká izotropní rozdělovací funkci, bude celkový tok částic na stěnu dán vztahem (Hays et al., 1974) 1 1 Γ+ = nw v¯ + Γw , (3.44) 4 2 kde první člen vyjadřuje chaotický pohyb částic a druhý člen vyjadřuje korekci na usměrněný pohyb částic. Nyní definujeme pravděpodobnost zachycení částice na stěně γ = Γw /Γ+ (3.45) a s užitím této definice eliminujeme Γ+ v rovnici (3.44). Po jednoduché úpravě dostáváme γ nw v¯ Γw = . (3.46) 2−γ 2 Dosazením z Fickova zákona (3.43) konečně dostáváme okrajovou podmínku difúzní rovnice γ v¯ 1 ∇n|w · k = − (3.47) nw 2 − γ 2D Uvedené výsledky nyní aplikujeme na řešení konkrétního fyzikálního problému, válcové výbojové trubice. Uvažme válec o výšce h0 a poloměru podstavy r0 . Potom řešení Poissonovy rovnice s nulovou Dirichletovou okrajovou podmínkou nabývá ve válcových souřadnicích tvaru 2,405r πz J0 , (3.48) n0 = n0 (0, 0) cos h0 r0

kde J0 je Besselova funkce nultého řádu. Uvádíme zde pouze nejnižší difúzní mod. Počátek soustavy souřadné je umístěn uprostřed válce. Difúzní délka tohoto řešení je určena vztahem 2 2 1 π 2,405 = + (3.49) Λ20 h0 r0 Okrajovou podmínku (3.47) lze nyní splnit dle (Chantry, 1987) přeškálováním výše uvedeného řešení. Tím se samozřejmě změní taktéž hodnota difúzní délky Λ a jak také ukazuje (Chantry, 1987), lze novou hodnotu Λ s chybou maximálně 11 % aproximovat vztahem Λ2 = Λ20 + l0 λ . (3.50) Zde l0 vyjadřuje v jistém smyslu rozměr aparatury a pro kulový, válcový a kvádrový tvar aparatury je definován vztahem V , (3.51) A kde V označuje objem aparatury a A je povrch aparatury. Veličina λ vychází z okrajové podmínky (3.47) a platí pro ni l0 =

λ=

2 − γ 2D . γ v¯

(3.52)

3.5. ČÁSTICOVÉ MODELOVÁNÍ

21

Dosazením tohoto řešení do difúzní rovnice a vyhodnocením Laplaceova operátoru dostaneme rovnici ∂n D = − 2n . (3.53) ∂t Λ Člen na pravé straně rovnice vyjadřuje úbytek částic vlivem difúze, avšak v kontextu chemické kinetiky je formálně shodný se členem úbytku od rovnice prvního řádu s rychlostní konstantou D/Λ2. Docházíme tak k závěru, že vliv difúze můžeme zahrnout do modelu chemické kinetiky zavedením fiktivní reakce prvního řádu s rychlostní konstantou k danou rovnicí 1 Λ2 V 2−γ = 0 +2 (3.54) k D A v¯γ

3.5

Částicové modelování

Částicové modely, jak již název napovídá, přistupují k plazmatu jako k souboru částic. Jde tedy o modely popisující plazma pomocí mikroskopických veličin, ze kterých potom můžeme zjistit prakticky libovolnou informaci o studovaném plazmatu. Dříve jsme uvedli, že všechny modely plazmatu lze interpretovat jako různé metody řešení Boltzmannovy rovnice, a ani částicové modely nejsou výjimkou. Chceme-li například vypočítat vývoj rozdělovací funkce f0 (x, v) v čase, můžeme daný výpočet provést metodou Monte Carlo: Nejprve vytvoříme náhodný soubor systémů částic. Ze znalosti f0 určíme statistické váhy, neboli pravděpodobnosti vzniku těchto systémů. Dále částicově vypočteme vývoj oněch systémů a střední hodnota přes mnoho vzorků potom bude konvergovat k řešeni Boltzmannovy rovnice. Při náhodném generování částicových systémů obvykle užíváme tzv. „importance sampling“, což v ideálním případě znamená, že systémy generujeme s pravděpodobností rovnou pravděpodobnosti vzniku danou rozdělovací funkcí. Díky tomu se věnujeme pouze statisticky významným variantám a můžeme tak ušetřit mnoho výpočtů. Více o této metodě viz (Landau – Binder, 2003). Existuje ještě alternativní interpretace částicových modelů jako řešení Boltzmannovy rovnice. Představme si, že bychom chtěli numericky řešit Boltzmannovu rovnici. Postup je takový, že musíme provést diskretizaci fázového prostoru. Při výpočtu časového vývoje se potom hustota pravděpodobnosti přesouvá ve fázovém prostoru v souladu s pohybovými rovnicemi a tento tok pravděpodobnosti odpovídá pohybu částic v simulaci. Nyní objasníme principy částicových modelů. Obvyklý přístup je takový, že v závislosti na čase počítáme trajektorie částic, jež jsou určeny Newtonovými pohybovými zákony. Působící síla, která vystupuje v pohybové rovnici, je dána součtem vnější vtištěné síly a vnitřní síly vytvářené samotnými částicemi. Právě výpočet onoho silového působení bývá co do náročnosti kritickou částí modelu. Uvažujeme-li totiž párové interakce mezi částicemi, je náročnost výpočtu silového působení v celém systému úměrná kvadrátu počtu částic, neboli asymptotická složitost algoritmu je třídy O(N 2 ). Při běžném počtu částic řádově 106 potom

3.5. ČÁSTICOVÉ MODELOVÁNÍ

22

dochází k neúnosnému zpomalení simulace. Při modelování pevných látek se tento problém řeší započtením síly pouze od částic v nějakém malém okolí, tento trik však v plazmatu kvůli nelokální povaze elektromagnetické interakce nemůžeme použít. Používá se však řada metod pro urychlení tohoto výpočtu, podle kterých můžeme modely rozdělit do několika skupin: Stromové algoritmy Mocnou metodou výpočtu silového působení jsou takzvané stromové algoritmy. Princip jejich funkce spočívá v tom, že interakci na krátké vzdálenosti počítají přesněji (obvykle jako obyčejnou párovou interakci), zatímco interakci se vzdálenějšími oblastmi počítají pouze přibližně, vystředovanou přes větší oblasti plazmatu. Tyto algoritmy přinášejí výrazné urychlení výpočtu síly, nebot’ jejich asymptotická složitost je třídy O(N log N) pro algoritmus Barnese-Huta (Barnes – Hut, 1986) nebo dokonce O(N) pro rychlou multipólovou metodu (Greengard – Rokhlin, 1987). Implementací a porovnáním stromových algoritmů se zabývá například (Entlicher, 2004). Nevýhodou těchto modelů je náročnější implementace, odměnou nám však je kromě vysoké rychlosti i relativní přesnost ve srovnání s particle-in-cell modely, které uvedeme dále. Ewaldova sumace Metoda Ewaldovy sumace je založena na rozkladu interakčního potenciálu částic na krátkodosahovou a dalekodosahovou složku. Síla vyplývající z krátkodosahové interakce je vypočtena klasicky jako párová interakce. Naproti tomu dalekodosahovou interakci vypočítáváme ve Fourierově obrazu systému. Nejprve určíme nábojovou hustotu z rozložení částic, provedeme Fourierovu transformaci, vypočteme interakční potenciál a transformujeme zpět. Tvar dalekodosahového interakčního potenciálu volíme hladký, abychom mohli zanedbat vysoké prostorové frekvence při výpočtu. Různé varianty tohoto algoritmu vykazují složitost O(N log N) až O(N 3/2 ). Particle-in-cell Ve fyzice plazmatu je zřejmě nejrozšířenější metoda particlein-cell (PIC). Ta je do jisté míry podobná Ewaldově metodě, nebot’ také převádí interakci částic na interakci s mříží. Liší se však v tom, že nerozlišuje mezi krátkodosahovou a dalekodosahovou interakcí. Veškeré silové působení mezi částicemi tedy určujeme tak, že vypočítáme nábojovou hustotu, vypočteme potenciál (podle potřeby můžeme použít při výpočtu Fourierovu transformaci) a z potenciálu potom určíme síly. Složitost tohoto algoritmu je opět O(N log N). V této práci modelujeme právě s využitím metody PIC a jejímu podrobnému popisu je proto věnována kapitola 3.6. Neselfkonzistentní algoritmy Zvláštní třídu částicových algoritmů tvoří neselfkonzistentní algoritmy. Ty určují působící sílu nezávisle na aktuálním rozložení částic. Síla je určena potenciálem, který je známý z jiné simulace nebo jde o nějaký teoretický odhad. Tyto simulace neumožňují výpočet potenciálu, avšak jsou vhodné například k výpočtu energetických a úhlových rozdělení. Obvykle

23

3.6. PARTICLE-IN-CELL

je využíváme v kombinaci s jiným, at’ už spojitým nebo částicovým modelem. Výhodou této metody je, že díky absenci interakce částic stačí provádět simulaci jednočásticově, což výrazně redukuje pamět’ové i výpočetní nároky.

3.6

Particle-in-cell

Princip modelu particle-in-cell byl nastíněn již v předchozím odstavci, nyní probereme detaily implementace. Běh modelu můžeme znázornit schématem 3.1, které zpodobňuje úkony provedené při výpočtu jednoho časového kroku ∆t. Nejprve z poloh částic určíme rozložení nábojové hustoty ρ, dále numericky vyřešíme Poissonovu rovnici pro známé ρ, z výsledného potenciálu φ určíme numerickou derivací elektrické pole E a se znalostí elektrického pole potom můžeme řešit pohybové rovnice. Spojitý prostor

Vzorkování

Mříž

-

Dynamika částic Newtonovy rovnice Monte Carlo t → t + ∆t

-

Výpočet ρij

-

1

Interpolace

Výpočet gradientu E ij = −∇Uij

4

Výpočet F (x)

Řešení Poissonovy rce ∆Uij = −ρij /ε0

3

2

Obrázek 3.1: Schéma jednoho časového kroku v modelu particle-in-cell

Uvedené kroky dobře ilustrují princip algoritmu, což je interakce pole známého na diskrétní prostorové mříži s částicemi, které se pohybují ve spojitém prostoru. Výměna informace mezi polem a částicemi probíhá v krocích 1 a 3. V kroku 1 jde o vzorkování spojité funkce v diskrétních bodech, zatímco krok 3 zahrnuje interpolaci funkce definované na mříži do poloh částic. V následujícím textu probereme detaily modelu particle-in-cell. Symboly hx a hy označují vzdálenosti sousedních mřížových bodů v souřadnici x a y. Dále definujeme mřížový vektor xij ≡ (ihx , jhy ). Veličiny definované na mříži označujeme obecně symbolem Xijk = X(ihx , jhy , k∆t). Pro jednoduchost v textu vynecháváme nepodstatné indexy.

3.6.1

Vzorkování a interpolace

Při výpočtu nábojové hustoty na mříži postupujeme tak, že každému mřížovému bodu xij přiřadíme oblast na okolí tohoto bodu, neboli mřížovou buňku. Hustotu ρij v bodě xij potom vypočteme jako podíl náboje obsaženého v této buňce a objemu buňky Vij . Pokud bychom uvažovali pouze bodové částice, dojde k tomu,

24


že při spojitém pohybu částice se nábojová hustota bude skokově měnit při přechodu mezi jednotlivými buňkami a vznikne tak nefyzikální šum. Tento problém eliminujeme nahrazením bodové částice nábojovým oblakem. Tvar nábojového oblaku udává funkce S(x). Nábojová hustota pocházející od částice o náboji q v bodě x′ je potom dána vztahem ρ(x) = qS(x − x′ ) .

(3.55)

V případě bodové částice platí jednoduše S(x) = δ(x). Nyní požadujme minimalizaci chyby výpočtu potenciálu na velkých vzdálenostech, dále aby náboj na mříži a tím pádem i síla byly hladkou funkcí polohy a nakonec aby náš model zachovával hybnost, což jsou rozumné požadavky vyplývající z chování reálného systému. Lze ukázat (Hockney – Eastwood, 1988), že z těchto požadavků vyplývá celá hierarchie možných tvarů nábojového oblaku a odpovídajících schémat pro přiřazování náboje na mříž. Nejnižším prvkem této hierarchie je opět bodová částice. V tomto schématu se celý náboj přiřazuje nejbližšímu mřížovému bodu, a proto se nazývá NGP (Nearest Grid Point). Chyba v určení potenciálu touto metodou je řádu ∆p , což je vzdálenost částice od nejbližšího mřížového bodu. Síla vypočtená touto metodou není spojitou funkcí polohy částice. Dalším členem je schéma CIC (Cloud in Cell). Zde má částice tvar obdélníkové funkce ⊓(x) definované v jedné dimenzi jako  |x| > 1/2  0 1/2 |x| = 1/2 ⊓(x) = (3.56)  1 |x| < 1/2 . Rozměry částice jsou shodné s charakteristickou funkcí mřížové buňky, tedy x S(x) = ⊓ . (3.57) hx

Zobecnění do dvou dimenzí je součin jednorozměrných funkcí x y S(x) = ⊓ ⊓ . hx hy

(3.58)

Chyba v určení potenciálu je řádu ∆2p . Síla je již spojitou funkcí polohy částice, avšak její derivace není. Funkce S zasahuje v jednorozměrném případě do dvou nejbližších buněk, a proto se náboj rozpočítává na dva nejbližší mřížové body. V hierarchii dále následuje schéma TSC (Triangular Shaped Cloud). Jak již název napovídá, tvar nábojového oblaku je určen trojúhelníkovou funkcí definovanou jako 0 |x| > 1 ∧(x) = . (3.59) 1 − |x| |x| ≤ 1

25


Máme tedy oblak tvaru

x S(x) = ∧ hx

.

(3.60)

Chyba v určení potenciálu je řádu ∆3p . Síla je nyní spojitou funkcí polohy částice včetně své první derivace. Funkce S zasahuje do tří buněk, a proto se náboj rozpočítává na tři nejbližší mřížové body. Jak lze vypozorovat z prvních třech členů této posloupnosti, tvar nábojového oblaku lze vypočíst jako vícenásobnou konvoluci obdélníkových funkcí. Platí totiž vztahy SPIC = H1 ⊓ ∗SNGP , STSC = H1 ⊓ ∗SPIC . S každým dalším členem se o jednotku zvyšuje řád přesnosti dalekodosahových sil a hladkost silového působení se zvětšuje o jednu derivaci v závislosti na poloze. Na druhou stranu se však také zvětšuje „rozmazání“ nábojového oblaku a tím pádem roste počet buněk, do kterých se náboj rozpočítává a také se snižuje přesnost výpočtu pro malé vzdálenosti. Chceme-li omezit rozmazání nábojové hustoty u metod vyššího řádu, je nutné aplikovat na nábojovou hustotu před výpočtem potenciálu „zaostřovací operátor“. Volba vhodného schématu je kompromisem mezi výpočetní náročností a přesností výpočtu. Ve fyzice plazmatu se často dává přednost rychlosti výpočtu před přesností, nebot’ výpočetní náročnost modelů je vysoká, a navíc nedokonalá znalost vstupních parametrů často způsobuje podstatně větší chybu než použité numerické metody. Proto je obvyklou volbou metoda CIC a nejinak tomu bude i v této práci. Na celou problematiku sumace náboje existuje ještě alternativní pohled. Můžeme totiž zadefinovat váhovou funkci W (x) jako příspěvek částice s jednotkovým nábojem k hodnotě náboje v mřížovém bodě o vzdálenosti x. Funkci W (x) určíme jako integrál váhové funkce posunuté do bodu x přes objem mřížové buňky. V případě, že váhová funkce je sudá a mříž je kartézská, můžeme ji symbolicky zapsat v jednorozměrném případě pomocí konvoluce W = ⊓∗S.

(3.61)

Zobecnění vztahu pro více rozměrů je zřejmé. Nábojovou hustotu ρij v bodě xij potom vypočteme jako N 1 X ρij = W (xij − xk )qk , Vij k=1

(3.62)

kde index k označuje částici, N je počet částic a xk je poloha k-té částice. Funkce na pravé straně je definovaná na celé pracovní oblasti a uvedený výraz tedy můžeme chápat jako vzorkování funkce v mřížových bodech xij . Na druhou stranu, při interpolaci hodnot elektrostatického potenciálu předpokládáme, že interpolovaná hodnota je lineární kombinací hodnot v mřížových bodech. Potom můžeme také pro interpolaci zavést váhovou funkci W ′ (x), která

26


bude udávat interpolovanou hodnotu výrazem XX E(x) = W ′ (x − xij )Eij . i

(3.63)

j

Je vhodné zvolit stejné váhové funkce pro interpolaci i vzorkování, tedy W = W ′ . Tato volba totiž, jak ukazuje (Hockney – Eastwood, 1988, s. 141), je spolu s použitím symetrických diferenčních schémat postačující podmínkou pro to, aby v modelu platil zákon zachování hybnosti.

3.6.2

Řešení Poissonovy rovnice

Poissonovu rovnici řešíme metodou sítí (Rektorys et al., 2000). Náš model využívá kartézský souřadnicový systém, a také použitá sít’ má čtvercové buňky. To značně zjednodušuje veškeré výpočty související s mříží. Pro řešení Poissonovy rovnice ρ ∆U = − (3.64) ε0 je nutné vyjádřit operátor druhé derivace v diferenční formě. V našem modelu používáme diferenční schémata přesnosti O(h2 ) fi+1/2 − fi−1/2 dfi ≈ , dx h

d2 fi fi−1 − 2fi + fi+1 ≈ . 2 dx h2

(3.65)

Dosazením tohoto schématu převedeme Poissonovu rovnici na diferenční rovnici a její aplikací na mřížové body dostaneme soustavu nehomogenních lineárních rovnic pro hodnoty potenciálu. V našem případě dvourozměrného modelu na kartézské mříži s konstantním h půjde o rovnice ρij 4Uij − Ui−1,j − Ui+1,j − Ui,j−1 − Ui,j+1 = . 2 h ε0

(3.66)

Tyto rovnice platí pro všechny body, které neleží na okraji pracovní oblasti a v jejichž okolí je kartézská mříž. Předmětem našeho modelu bude objem plazmatu s vnořenou vodivou sondou. Na rozhraní pracovní oblasti s neporušeným plazmatem i na rozhraní se sondou uvažujeme Dirichletovu okrajovou podmínku pro potenciál. Splnění této podmínky na okrajích plazmatu je triviální, zvolíme-li okraje tak, aby procházely mřížovými body. Na rozhraní se sondou musíme za tímto účelem zavést nepravidelnou sít’, aby mřížové body opět ležely na rozhraní. Diferenční schéma pro derivaci na nepravidelné síti je dáno vztahem d2 fi fi+1 2 1 1 fi−1 , (3.67) ≈ − fi ( + ) + dx2 l+r r r l l kde l = xi − xi−1 a r = xi+1 − xi jsou vzdálenosti levého a pravého mřížového bodu. Z tohoto schématu vyplývá Poissonova rovnice na nepravidelné pravoúhlé

27


mříži ve tvaru 2 2 rl + tb Ui+1,j Ui−1,j Ui,j+1 Ui−1 ρij + −2 + + Uij = , (3.68) l+r r l t+b t l rltb ε0 kde t = yi,j+1 − yij a b = yij − yi,j−1 jsou vzdálenosti horního a dolního přilehlého mřížového bodu, jak ilustruje obrázek 3.2. j +1 t l

j

r b

j −1 i−1

i

i+1

Obrázek 3.2: Ilustrace k diferencování na nepravidelném rozhraní

Rovnice (3.66) a (3.68) v kombinaci s Dirichletovou okrajovou podmínkou Uokraj = U 0 tvoří soustavu rovnic určující výsledný elektrostatický potenciál. Tuto soustavu můžeme formulovat také maticově ve tvaru Aφ = ρ ,

(3.69)

kde matice A má ve dvou dimenzích nenulové prvky pouze na hlavní diagonále a na čtyřech vedlejších diagonálách. Zde dočasně značíme potenciál φ, aby nedošlo k záměně s maticí U vystupující v LU dekompozici. Tato soustava nabývá v případě vícerozměrné mříže obrovských rozměrů. Strany mříže mívají totiž typicky stovky bodů. Pro dvourozměrnou mříž potom celkový počet mřížových bodů, a tedy i rovnic nabývá řádově 104 . Matice této soustavy tedy má 108 až 109 koeficientů. V závislosti na implementaci může řešení Poissonovy rovnice tvořit významnou část výpočtu. Pro svoji jednoduchost se často k řešení používá metoda superrelaxace, neboli SOR (Successive Over Relaxation) (Ralston, 1973; Press et al., 2002). Tato metoda však zdaleka není optimální z hlediska rychlosti výpočtu. Navíc jde o metodu iterační a její přesnost je tedy dána počtem iterací. Metodu SOR využívala první verze našeho modelu. V současnosti je jednou z nejsilnějších metod řešení metoda založená na LU dekompozici matice soustavy. LU dekompozice, respektive trojúhelníkový rozklad, je typ maticové dekompozice, při které čtvercovou matici A rozložíme na součin A = LU , (3.70)

28


kde matice L je dolní trojúhelníková matice a U je horní trojúhelníková matice. Jednu z těchto matic můžeme zvolit tak, aby měla jedničky na hlavní diagonále, a potom je tento rozklad jednoznačný. Rovnici (3.69) užitím tohoto rozkladu převedeme na dvě maticové rovnice Uφ = x Lx = ρ .

(3.71) (3.72)

Díky tomu, že matice těchto soustav jsou trojúhelníkové, lze je snadno řešit zpětnou substitucí. Výpočet matic L a U je poměrně náročný, nebot’ pro obecné matice řádu n vyžaduje 2/3n3 operací. Ovšem v případě, že potřebujeme řešit rovnici (3.69) mnohokrát pro stejnou matici, vyplatí se již tuto dekompozici na začátku výpočtu provést. Matice soustavy A vyskytující se při řešení Poissonovy rovnice v modelu PIC se vyznačuje jednak svou velikostí, ale také tím, že drtivá většina jejích prvků je nulových. Jde o takzvanou řídkou matici. Inverze řídké matice sice již obecně není řídká matice, což znemožňuje její výpočet pro velké soustavy rovnic, avšak LU dekompozice obvykle řídkost matice zachovává, a pro dvourozměrné modely tak není problém ji vypočítat na běžném hardwaru. V našem modelu využíváme pro řešení Poissonovy rovnice knihovnu UMFPACK (Davis, 2004a). Tato knihovna je optimalizovaná právě pro práci s řídkými maticemi. Pro zvýšení přesnosti a rychlosti výpočtu je matice A před výpočtem škálována diagonální maticí R a je provedena permutace jejích řádků a sloupců maticemi P a Q. Matice P RAQ je potom rozložena na součin LU . Volbu transformačních matic taktéž obstarává knihovna UMFPACK a její algoritmy jsou detailně popsány v (Davis – Duff, 1999; Davis, 2004b). Pro maticové operace dokáže UMFPACK využít optimalizované knihovny BLAS, konkrétně jsme zvolili pro akademické účely volně dostupnou implementaci GotoBLAS (Goto). Jak ukazuje (Pekárek et al., 2007), je tato metoda řádově 10× rychlejší než superrelaxace. V typickém dvourozměrném PIC modelu se tak výpočetní náročnost řešení Poissonovy rovnice stává zanedbatelnou a další optimalizace v této oblasti tím pádem nemají význam. Jiná situace ovšem nastává v případě třírozměrných modelů, kde řešení Poissonovy rovnice stále představuje významný problém.

3.6.3

Dynamika částic

V každém časovém kroku musíme vypočítat změnu polohy a rychlosti částice během času ∆t. Tyto veličiny jsou v našem modelu určeny Newtonovými rovnicemi d2 x F = . 2 dt m

(3.73)

29


Nejjednodušší metodou numerického řešení této rovnice je Eulerův algoritmus, který popisuje přechod o jeden časový krok následujícími rovnicemi F k = −q∇U xk+1 = xk + v k ∆t 1 v k+1 = v k + F k ∆t . m

(3.74)

Přesnost Eulerova algoritmu je řádu O(∆t) a pro svoji nízkou přesnost není často používán. Z algoritmů s přesností druhého řádu je velmi rozšířený Verletův algoritmus a leap-frog algoritmus. Rychlostní Verletův algoritmus je definován vztahy Fk 2 ∆t m F k+1 + F k ∆t . = vk + 2m

xk+1 = xk + v k ∆t + v k+1

(3.75)

Nevýhodou tohoto algoritmu je, že vyžaduje znalost síly ve dvou časech. V částicových modelech plazmatu je často používán algoritmus leap-frog ve tvaru F k = −q∇U

1 k F ∆t m = xk + v k+1/2 ∆t .

v k+1/2 = v k−1/2 + xk+1

(3.76)

Tento algoritmus je jednoduchý a vyžaduje dokonce méně operací než Eulerův algoritmus. Jeho přesnost řádu O(∆t2 ) je přitom pro model PIC dostačující. Algoritmus sice vyžaduje počáteční hodnotu rychlosti v čase −∆t/2, avšak náš model generuje počáteční rychlosti nezávisle na poloze a čase, a proto nepůsobí tento požadavek vzhledem k počátečním podmínkám žádné obtíže. Problém ovšem může nastat při simulaci plazmatu se srážkami, nebot’ v okamžiku srážky potřebujeme znát rychlost. Při srážce se navíc rychlost změní a museli bychom vypočítat v k+1/2 nějakou alternativní metodou, abychom v dalším kroku mohli navázat leap-frog algoritmem. V případě plazmatu s malým množstvím srážek respektive s krátkým časovým krokem tento problém zanedbáváme. V sekci 6.4 ovšem odvodíme model s nekonstantním časovým krokem pro případ, že srážky mají dominantní vliv na dynamiku částic. V onom modelu použijeme Eulerovu a modifikovanou Verletovu metodu integrace, nebot’ pohyb částic je především řízen srážkami a přesnost integrace není tak důležitá.

3.6.4

Srážky

Pod pojmem srážky rozumíme především různé kvantové procesy, jako jsou excitace různých elektronických, vibračních a rotačních stavů, dále také pružné srážky


30

nabitých částic s neutrálními částicemi a přenos náboje mezi ionty a neutrály. Coulombická interakce je sice v modelu již zahrnuta, avšak vlivem interpolace na mříž jsou krátkodosahové účinky této interakce vyhlazeny. Proto musíme, v případě že je coulombický rozptyl do velkých úhlů významný, zahrnout tuto interakci do srážek. V modelu typu particle-in-cell simulujeme srážky metodou Monte Carlo (MCC neboli Monte Carlo Collisions). Důvodem je jednak náhodná povaha kvantových procesů, a potom také fakt že neutrální atomy a molekuly v modelu nejsou částicemi, nebot’ hlediska PIC modelu nemají neutrální částice žádný vliv na potenciál. Princip metody MCC spočívá v simulaci srážek pro náhodně vybrané částice v náhodných okamžicích. Také druh probíhající interakce, rychlost interagujícího partnera a výsledek srážky jsou náhodné. Rozdělovací funkce uvedených náhodných veličin však musíme pečlivě volit s ohledem na modelovaný systém. Vzhledem k tomu, že metody simulace srážek nejsou v literatuře kvalitně popsány, provedeme v této práci podrobný rozbor metodiky s odvozením exaktních metod pro simulaci srážek v sekci 6.2.

3.6.5

Generování pseudonáhodných čísel

V souvislosti s metodou Monte Carlo jsme zmínili, že srážky modelujeme pro náhodně vybrané částice v náhodných okamžicích. Realizovat takový náhodný výběr na počítači, který můžeme v ideálním případě považovat za deterministický stroj, představuje komplikovaný problém. Jednou z možností je získávat náhodná data z fyzikálního procesu náhodné povahy jako je třeba jaderný rozpad. Levnější a rychlejší alternativou jsou softwarové generátory pseudonáhodných čísel. Ty generují určitou složitou, ovšem jednoznačně určenou, posloupnost celých čísel s velkou periodou. Snahou je, aby tato posloupnost byla pomocí různých statistických testů nerozeznatelná od náhodné posloupnosti s rovnoměrným rozdělením. Velmi rozšířená je například sada testů Diehard (Marsaglia, 1995). Autor těchto testů George Marsaglia uvádí také řadu generátorů náhodných čísel (Marsaglia, 1999), patří mezi ně například rychlý generátor s posuvnými registry uint32 jz, jsr = seed; uint32 shr3() { return jz=jsr, jsr^=(jsr<<13), jsr^=(jsr>>17), jsr^=(jsr<<5), jz+jsr; }; kde jsme použili syntaxe jazyka C a typ uint32 označuje 32 bitové celé číslo bez znaménka. Další variantou uvedenou v (Marsaglia, 1999) je také generátor mwc() typu multiply-with-cary definovaný kódem uint32 w=seed1, z=seed2; #define znew (z=36969*(z&65535)+(z>>16))


31

#define wnew (w=18000*(w&65535)+(w>>16)) inline uint32 mwc(){ return (znew<<16)+wnew;}; Jde vlastně o kombinaci dvou generátorů náhodných čísel. Tento generátor úspěšně projde všemi testy Diehard (Marsaglia, 1999, 1995) a má periodu přibližně 260 . Generátory shr3() a mwc() používáme v našem modelu pro výpočet pseudonáhodných čísel. Pseudonáhodná čísla budeme v dalším textu označovat jen jako náhodná. Pokud chceme generovat náhodná čísla z nějakého složitějšího pravděpodobnostního rozdělení, získáváme je transformací vstupního rovnoměrného rozdělení. Nejjednodušší transformací je vydělení náhodného celého čísla jeho maximální dosažitelnou hodnotou. Tak získáme náhodné rozdělení na intervalu h0, 1i značené jako U(0, 1). Z rovnoměrného rozdělení na intervalu h0, 1i lze získat další rozdělení metodou inverzní distribuční funkce. Princip metody spočívá v tom, že analyticky či numericky určíme inverzní funkci ke kumulativní distribuční funkci požadovaného rozdělení (CDF neboli Cumulative Distribution Function). Náhodný výběr z U(0, 1) potom dosadíme do inverzní CDF a výsledná funkční hodnota je již náhodnou proměnnou podléhající požadovanému rozdělení. Metodou inverzní distribuční funkce snadno odvodíme vztah pro generování náhodné veličiny s exponenciálním rozdělením, které určuje střední dobu mezi následujícími nezávislými náhodnými událostmi v Poissonovském procesu. Náhodnou veličinu ξ z exponenciálního rozdělení f (ξ) = −1/λ exp(−ξ/λ) určíme jako ξ = −λ ln γ , γ ∈ U(0, 1) . (3.77)

Pomocí složitější Boxovy-Mullerovy transformace (Box – Muller, 1958) můžeme také generovat dvě čísla ξ1,2 z normálního rozdělení s jednotkovým rozptylem p ξ1 = cos(2πγ1) −2 log γ2 , γ1 , γ2 ∈ U(0, 1) (3.78) p ξ2 = sin(2πγ1 ) −2 log γ2 . Rovnice (3.78) však vyžadují výpočet čtyř matematických funkcí a jejich vyhodnocení je tedy poměrně pomalé. Proto je výhodnější použití Boxovy-Mullerovy transformace v Marsagliově polárním tvaru float x1, x2, w, y1, y2; do { x1 = 2.0 * uni() - 1.0; x2 = 2.0 * uni() - 1.0; w = x1 * x1 + x2 * x2; } while ( w >= 1.0 ); w = sqrt( (-2.0 * ln( w ) ) / w ); y1 = x1 * w; y2 = x2 * w;

3.7. HYBRIDNÍ MODELY

32

Tato transformace vyžaduje pouze dvě vyhodnocení matematických funkcí pro generování dvou náhodných čísel. Existují však ještě podstatně efektivnější metody generování náhodných čísel z normálního rozdělení. Jednou z takových metod je algoritmus Ziggurat (Marsaglia – Tsang, 2000), jenž je vhodný pro generování libovolné náhodné proměnné s monotónní rozdělovací funkcí. Článek (Marsaglia – Tsang, 2000) uvádí implementaci tohoto algoritmu pro generování náhodných čísel z exponenciálního a normálního rozdělení, které používáme v této práci. Generátor typu Ziggurat vyžaduje v 99% většině případů pouze generování jednoho náhodného celého čísla a načtení dvou čísel z paměti, pročež je velmi rychlý.

3.7

Hybridní modely

Jak již bylo uvedeno, techniky modelování nelze striktně rozdělit na částicové a spojité. Průnik těchto oblastí tvoří takzvané hybridní modely. Cílem těchto modelů je využití spojitého přístupu k urychlení výpočtu a částicového modelu pro zpřesnění popisu kinetiky částic. Jako hybridní modely jsou někdy také označovány modely typu particle-incell využívají Monte Carlo simulaci srážek, tento typ hybridního modelu zde však nepopisujeme, nebot’ byl probrán v předchozí sekci jde prakticky o standard v particle-in-cell simulacích. Hybridní modely existují v mnoha variantách a často jsou vytvářeny speciální jednoúčelové varianty hybridního modelu pouze pro konkrétní aplikace. Asi nejčastější přístup je takový, že některé složky plazmatu, například určité ionty nebo neutrální částice, popisujeme spojitě, zatímco jiné složky, například elektrony, popisujeme částicově. Například práce (Cartwright et al., 1998) uvádí model kombinující PIC model iontů kombinovaný popisem elektronů pomocí Boltzmannova vztahu n(x) = n0 exp(eU(x)/kB T ), který řeší problém rozdílných časových škál iontů a elektronů a umožňuje simulaci na delších časových škálách. Alternativně můžeme například nízkoenergetické elektrony popisovat spojitě a vysokoenergetické elektrony, které díky své dlouhé střední volné dráze zodpovídají za určité nelokální jevy, popisovat částicově (Kolev – Bogaerts, 2004). Další variantou je simulace určité oblasti prostoru spojitě a částicová simulace v některých kritických oblastech. Nakonec uveďme ještě variantu, která kombinuje spojitý a částicový model v iteračním hledání stacionárního stavu, jež je popsána v (Sommerer – Kushner, 1992) a zabývá se jí také (Bartoš, 2007).

Kapitola 4 Cíle práce Cílem této práce je především rozvoj metod modelování plazmatu a optimalizace modelů pro výpočty při středních tlacích. Protože při středních tlacích hrají velmi významnou roli srážky částic, zaměřujeme se jednak na zpřesnění výpočtu srážek a jednak na celkové urychlení výpočtu. V oblasti zpřesňování výpočtu navazujeme na bakalářskou práci (Roučka, 2006), kde byly rozebrány problémy jedné z metod simulace srážek a zde si klademe za cíl nalezení vhodnější metody. Pro zrychlení výpočtu se pokusíme modifikovat algoritmus particle-in-cell v kombinaci s Monte Carlo výpočtem srážek a přinést tak významné zkrácení výpočetního času. Dále v práci rozvíjíme model chemické kinetiky kyslíkového plazmatu vzniklý v rámci seminární práce na KFPP. Model bude doplněn o výpočet elektronové rozdělovací funkce a aktuálnější hodnoty rychlostních konstant.

33

Kapitola 5 Kinetický model kyslíkového plazmatu Matematická stránka kinetického modelu je popsána v sekci 3.4. Zásadní práce při tvorbě modelu chemické kinetiky však spočívá v sestavení vhodné sady rychlostních konstant. Mnoho důležitých rychlostních konstant lze vypočítat přímo z účinných průřezů dané interakce dle vztahu (3.27). Pro výpočet střední hodnoty dle (3.27) je nutná znalost elektronové rozdělovací funkce a tu nyní určíme rozvojem rozdělovací funkce popsaným v sekci 3.2.

5.1

Výpočet elektronové rozdělovací funkce

Náš model vychází ze sady účinných průřezů sestavené v (Phelps). Vzhledem k tomu, že měření účinných průřezů jsou často neúplná, případně jejich chyba dosahuje desítek procent, je při tvorbě modelu nutné pečlivé vybírat, a někdy dokonce „odhadovat“ účinné průřezy simulovaných interakcí. Obvyklý postup je takový, že je sestaven model na základě dostupných experimentálních dat a následně jsou hodnoty účinných průřezů korigovány tak, aby model dosahoval shody s pozorovanými makroskopickými veličinami jako například driftová rychlost a pohyblivost elektronů. Tento postup můžeme interpretovat jako upřesňování naměřených hodnot. Jednou z takto vyladěných sad účinných průřezů je právě (Phelps), kterou použiji při výpočtu elektronové rozdělovací funkce. Uvažované reakce jsou zapsány v tabulce 5.1. Účinné průřezy pro vibrační excitace jsou graficky znázorněny v obrázku 2.4, ostatní reakce jsou znázorněny v obrázku 2.5. S uvážením výše uvedených interakcí byly poté pomocí programu ELENDIF (Morgan – Penetrante, 1990) vypočteny rozdělovací funkce elektronů při daném vnějším elektrickém poli. Ukázka výsledků pro různé intenzity elektrického pole je zakreslena v grafu 5.1. Rozdělovací funkce jsou v grafu vynásobeny funkcí E −1/2 , jak je v některých publikacích zvykem. Touto transformací se maxwellovské rozdělení redukuje na lineární závislost v logaritmickém měřítku a je tak možné snáze 34

35

5.2. KINETICKÝ MODEL

přenos hybnosti vibrační excitace excitace singletních stavů vyšší elektronické hladiny ionizace disociativní excitace

O2 + e− → O2 + e− O2 + e− → O2 (v = 1 . . . 4) + e− O2 + e− → O2 (a1 ∆g ) + e− − O2 + e− → O2 (b1 Σ+ g)+e − O2 + e → O2 (4,5, 6,0, 8,4, 9,97) + e− − − O2 + e− → O+ 2 +e +e O2 + e− → O + O + e− + hν(130 nm)

Tabulka 5.1: Interakce uvažované při výpočtu elektronové rozdělovací funkce v kyslíku.

porozumět výsledným grafům. Grafy byly vypočteny při obdobných podmínkách jako v práci (Capitelli et al., 1979) a výsledky jsou také velmi podobné, což může sloužit jako verifikace našich výpočtů.

5.2

Kinetický model

Vypočtená energetická rozdělení elektronů nám již dávají dostatek informace, abychom mohli vypočítat rychlostní konstanty probíhajících reakcí podle vztahu (3.27). Použitá sada chemických reakcí vychází z několika prací zabývajících se chemií kyslíkového plazmatu. Jako základní vodítko nám sloužila práce (Odehnal, 1996). Nejprve jsme sestavili model s rychlostními konstantami přesně podle této práce. Původní sadu více než 100 reakcí jsme dále redukovali na přibližně 40 nejdůležitějších reakcí tak, aby toto zanedbání nemělo vliv na výsledné koncentrace v ustáleném stavu. Jako míru důležitosti i-té reakce pro bilanci j-tého prvku jsme definovali veličinu Qij vztahem Qij =

Iij Rij + , Ij Rj

(5.1)

kde Rij je člen zániku prvku j způsobený reakcí i ve vztahu (3.22). Obdobně pro člen vzniku I. Dále jsme se pokusili získat nejaktuálnější hodnoty rychlostních konstant a zpřesnit tak výsledky modelu. Bohužel, znalost rychlostních konstant je velmi omezená. Často se setkáváme s tím, že rychlostní konstanty udávané v různých pramenech se několikanásobně liší. Jako ukázkový případ uvedeme například reakci O + O2 + O3 → 2O3 . (5.2)

Zde (Odehnal, 1996) uvádí k = 1 · 10−29 cm3 · s−1 , zatímco (Kossyi et al., 1992) 1,25 udává hodnotu k = 6,9 · 300 · 10−34 cm3 · s−1 , která při pokojové teplotě T nabývá hodnoty k ≈ 6,9 · 10−34 cm3 · s−1 . Dalším příkladem budiž reakce − O + O− 3 → 2O2 + e .

(5.3)


36

Obrázek 5.1: Elektronové rozdělovací funkce v O2 v závislosti na E/n.

Pro její rychlostní konstantu platí dle (Kossyi et al., 1992) k = 3 · 10−10 cm3 · s−1 , dále (Ionin et al., 2007) uvádí k = 1 · 10−13 cm3 · s−1 a konečně (Odehnal, 1996) udává k = 1 · 10−11 cm3 · s−1 . Udávané rychlostní konstanty se tedy napříč prameny rozcházejí o několik řádů a podobná situace, i když ne vždy tak dramatická, nastává i u mnoha dalších reakcí. Nepříjemnou vlastností kinetického modelu je jeho velká citlivost vůči změně byt’ i jen jediné rychlostní konstanty. Pravděpodobně vlivem nelinearity diferenciálních rovnic, které tvoří základ tohoto modelu, vykazuje tento systém do jisté míry chaotické chování, což znamená, že malá změna vstupních parametrů (rychlostních konstant) způsobí velkou změnu výsledků (koncentrace v ustáleném stavu). Při doplňování aktuálních rychlostních konstant do modelu jsme tak často narazili na problém, že po změně jedné rychlostní konstanty začal model dávat fyzikálně nesmyslné výsledky. Ukázalo se tedy, že sestavení smysluplného kinetického modelu kyslíkového plazmatu není možné bez důkladného porovnávání výsledků s experimentem. Rychlostní konstanty je potřeba ručně nastavit tak, aby se výsledky modelu shodovaly s experimentálními daty a teprve potom můžeme model použít ke zkoumání probíhajících fyzikálních procesů. V rámci této práce však podobná interakce s experimentem nebyla k dispozici, a především tvorba kinetického modelu nebyla její hlavní náplní. Proto jsme

37


O2 O O2 (∆) O2 (Σ) O3 O(D) O+ 2 e− O− O− 2 O− 3 O+

n [cm−3 ]

1015

1010

105

100 10−15

10−10

10−5

100

105

t [s] Obrázek 5.2: Časový vývoj koncentrací jednotlivých složek kyslíkového plazmatu při elektrickém poli E/n = 80 Td, teplotě plynu T = 300 K a tlaku p = 133 Pa. Složky jsou v legendě seřazeny podle koncentrace v rovnovážném stavu.

se omezili pouze na hrubou kontrolu výsledků, kterou jsme zajistili, aby poměry koncentrací v ustáleném stavu řádově odpovídaly publikovaným výsledkům například ve (Vahedi – Surendra, 1995; Odehnal, 1996). Výsledkem naší práce tedy je databáze účinných průřezů a funkčních závislostí a tabulek rychlostních konstant. Na základě vstupních parametrů probíhá automatický proces, který spočívá nejprve ve výpočtu elektronových rozdělovacích funkcí, na jejich základě jsou potom vyhodnoceny rychlostní konstanty a ty slouží jako vstupní data kinetického modelu, který nám nakonec vrátí časovou závislost koncentrace jednotlivých složek, případně koncentrace v ustáleném stavu. V tabulce 5.2 jsou uvedeny rychlostní konstanty našeho modelu vyhodnocené pro elektrické pole E/n = 80 Td a teplotu plynu T = 300 K. Na základě těchto konstant byl jako ukázka schopností kinetického modelu vypočten časový průběh koncentrací jednotlivých složek kyslíkového plazmatu. Výsledky jsou graficky znázorněny v grafu 5.2.

38


reakce O2 + e− → e− + O2 (a1 ∆g ) O3 + O2 (a1 ∆g ) → 2O2 + O(D) O + O2 (b1 Σ+ g ) → O + O2 − 1 e + O2 (a ∆g ) → O2 + e− O2 + O(D) → O + O2 (b1 Σ+ g) − − O2 + e → O + O O2 + e− → e− + O2 (b1 Σ+ g) − O2 + e− → O+ + 2e 2 O− + O2 (a1 ∆g ) → O + O− 2 O + e− → O+ + 2e− O2 + O+ → O + O+ 2 O2 + 2O2 (a1 ∆g ) → 2O3 O + O− → O2 + e− O3 + O− → O + O− 3 − O + O− 2 → O2 + O + O2 + e− → 2O 2O + O2 → 2O2 − O + O− 2 → O3 + e − O+ 2 + e → O + O(D) − 1 O + O2 (a ∆g ) → O3 + e− O2 + e− → 2O + e− 2O2 + O− → O2 + O− 3 3O → O + O2 O + O2 + O3 → 2O3 O2 + O(D) → O + O2 − O + O− 3 → O2 + O2 − O + O− 3 → 2O2 + e − + O2 + e → O + O + 2e− 2O + O2 → O + O3 O− + O+ 2 → O + O2 O3 + O2 (b1 Σ+ g ) → O + 2O2 − 1 O2 + O2 (a ∆g ) → 2O2 + e− O + O2 (a1 ∆g ) → O + O2 O2 + O3 + e− → O2 + O− 3 + → O + O O− + O 2 3 2 3 O3 + e− → O + O2 + e− O− + e− → O + 2e− 1 O2 (a ∆g ) + e− → O + O−

h

i

cm3 s−1 cm6 s−1 6,5 ·10−10

k

1 ·10−11 4 ·10−14 2,9 ·10−09 3,2 ·10−11 3,7 ·10−11 1,5 ·10−10 7 ·10−13 3,3 ·10−10 1,1 ·10−15 2 ·10−11 1 ·10−31 5 ·10−10 5,3 ·10−10 3,3 ·10−10 2,2 ·10−07 7,2 ·10−33 3 ·10−10 2 ·10−07 3 ·10−10 5,5 ·10−11 1,1 ·10−30 4,8 ·10−33 6,9 ·10−34 8 ·10−12 3,2 ·10−10 1 ·10−13 2,3 ·10−15 6,9 ·10−34 2 ·10−07 1,5 ·10−11 2 ·10−10 7 ·10−16 1,8 ·10−31 2 ·10−07 8,8 ·10−10 4 ·10−08 1,2 ·10−10

zdroj

pozn.

(Phelps) (Odehnal, 1996) (Ionin et al., 2007) (Odehnal, 1996) (Kossyi et al., 1992) (Eliasson, 1986) (Phelps) (Odehnal, 1996) (Stoffels et al., 1995) (Odehnal, 1996) (Kossyi et al., 1992) (Odehnal, 1996) (Kossyi et al., 1992) (Ionin et al., 2007) (Stoffels et al., 1995) (Odehnal, 1996) (Odehnal, 1996) (Ionin et al., 2007) (Odehnal, 1996) (Ionin et al., 2007) (Odehnal, 1996) (Vasiljeva et al., 2004) (Ionin et al., 2007) (Kossyi et al., 1992) (Vasiljeva et al., 2004) (Vasiljeva et al., 2004) (Ionin et al., 2007) (Eliasson, 1986) (Kossyi et al., 1992) (Kossyi et al., 1992) (Ionin et al., 2007) (Ionin et al., 2007) (Kossyi et al., 1992) (Kossyi et al., 1992) (Kossyi et al., 1992) (Ionin et al., 2007) (Ionin et al., 2007) (Stoffels et al., 1995)

EEDF

f (T ) EEDF EEDF

f (T )

f (T ) f (T ) f (T )

EEDF f (T ) f (T )

f (Te ) f (T )

Tabulka 5.2: Rychlostní konstanty interakcí v kyslíkovém plazmatu. V poznámce je uvedeno, jestli je rychlostní konstanta skutečně konstanta, nebo jestli byla vyhodnocena z funkční závislosti nebo účinného průřezu.

Kapitola 6 Částicový model pro vyšší tlaky Vzhledem k tomu, že se zabýváme plazmatem při středních a vyšších tlacích, musíme při našich výpočtech klást velký důraz na přesné zohlednění srážkových interakcí v plazmatu. Nyní provedeme teoretické zdůvodnění použitého modelu srážek typu Monte Carlo Kromě toho, že srážky mají velký vliv na výsledek simulace, tvoří jejich výpočet většinu doby výpočtu při středních tlacích. Proto také na konci této kapitoly vypočteme, jaká chyba vzniká při nadměrném zvyšování časového kroku simulace v zájmu urychlení výpočtu a navrhneme modifikaci metody particle-in-cell, která umožňuje podstatné prodloužení časového kroku.

6.1

Model založený na generování náhodné volné dráhy

Při modelování srážek částice se obvykle vychází z představy, že sledovaná částice se pohybuje mezi náhodně rozloženými nehybnými rozptylovými centry – částicemi. Pravděpodobnost, že se částice srazí na elementární dráze dx, neboli počet srážek na dráze dx, je p = nσ dx , (6.1) kde n je koncentrace rozptylových center a σ je účinný průřez srážky. Vyšleme-li do takového prostředí svazek částic o počáteční intenzitě I0 , bude pro intenzitu v závislosti na poloze platit I(x + dx) = I(x)(1 − nσ dx) .

(6.2)

Z tohoto vztahu vyplývá pro intenzitu diferenciální rovnice I ′ = −nσI, jejímž řešením je I = I0 exp(−nσx) . (6.3) Výše uvedený vztah můžeme po normování interpretovat jako rozdělení pravděpodobnosti očekávané volné dráhy částice. Zavedeme-li označení pro střední hodnotu 39

6.1. MODEL ZALOŽENÝ NA GENEROVÁNÍ NÁHODNÉ VOLNÉ DRÁHY

40

tohoto rozdělení

1 , σn můžeme rozdělení náhodné volné dráhy psát jako 1 l f (l)dl = exp − dl . λ λ λ≡

(6.4)

(6.5)

Pro pozdější referenci ještě definujeme střední dobu mezi srážkami τ≡

λ 1 = . v σnv

(6.6)

Účinný průřez interakce je obecně závislý na relativní rychlosti částic, a proto je také očekávaná volná dráha závislá na rychlosti. To přináší problémy, protože nedokážeme dopředu určit očekávanou volnou dráhu částice, jejíž rychlost se mění. Z toho důvodu přidáváme do modelu fiktivní reakci, jejíž účinný průřez je takový, aby totální účinný průřez byl konstantní. Následky této fiktivní interakce jsou však nulové, a proto se tato metoda se nazývá metoda nulové srážky. Náhodnou volnou dráhu potom můžeme generovat nezávisle na rychlosti. Pohyb částice tedy počítáme tak, že náhodně vygenerujeme volnou dráhu částice z rozdělení (6.5), dále počítáme trajektorii o délce l z Newtonových rovnic, poté započteme důsledky srážky a vše se znovu opakuje. V okamžiku srážky musíme náhodně generovat druhou interagující částici. V případě, že probíhá více různých druhů interakcí, počítáme střední volnou dráhu z totálního účinného průřezu, jenž je součtem všech účinných průřezů. O probíhající interakci potom rozhodujeme v okamžiku srážky, přičemž pravděpodobnosti jednotlivých srážek jsou úměrné jejich účinným průřezům. Tento postup je však ve fyzice plazmatu platný jen částečně. V úvodu jsme totiž zanedbali pohyb okolních částic a jak ukazuje (Šimek, 2006), může toto zanedbání vést k velké chybě obzvláště v případě pohybu iontů. Obrázek 6.1 ukazuje rychlostní rozdělení iontů v plazmatu vypočítané uvedeným algoritmem. Model uvažoval interakci přenosu náboje a pružnou srážku s neutrální částicí. Interagující partneři byli generováni z Maxwellova rozdělení. Nebyly uvažovány interakce s elektrony ani s vnějším elektrickým polem. V tomto modelu by se ionty měly uvést do rovnováhy s neutrálními částicemi, jejichž rozdělení je maxwellovské. Výsledky modelu však ukazují výrazné odchylky od maxwellovského rozdělení, dochází k nefyzikálnímu ochlazování. Důsledky tohoto ochlazování byly vysvětleny v (Roučka, 2006). Problém je jednak to, že očekávaná volná dráha se s klesající rychlostí snižuje a pro nehybnou částici je dokonce nulová a jednak, že interagující částice nejsou výběrem z Maxwellova rozdělení. Generováním částic z nemaxwellovského rozdělení bylo v (Roučka, 2006) toto ochlazování částečně potlačeno. V následující sekci odvodíme algoritmus, který tento problém zcela odstraňuje.

41

6.2. OBECNÝ MODEL SRÁŽEK

f (v) [10−3 · s · m−1 ]

2,5

Maxwellovo rodělení rozdělení iontů

2 1,5 1 0,5 0 -1000

-500

0 v [m ·

500

1000

s−1 ]

Obrázek 6.1: Rychlostní rozdělení iontů v plazmatu vypočtené nevhodným algoritmem, který zanedbává pohyb neutrálních částic. Výsledek je značně odlišný od očekávaného Maxwellova rozdělení.

6.2

Obecný model srážek

Budeme nyní zkoumat pohyb částice – iontu, v plynu neutrálních částic. Model lze samozřejmě použít na pohyb libovolné částice v libovolném plynu, například při výpočtu elektron-elektronových srážek. Sledovaná částice interaguje s plynem interakcí o účinném průřezu σ(v). Abychom mohli modelovat srážky metodou Monte Carlo, potřebujeme znát, jak často dochází ke srážkám a jaké je rozdělení interagujících částic. Uvažujme tedy iont v plynu. Náš iont budiž umístěn v počátku soustavy souřadné. Určíme nyní, kolik částic s rychlostí v z intervalu (vx , vx + dvx ) × (vy , vy + dvy ) × (vz , vz + dvz ) narazí v časovém intervalu (t, t + dt) na sledovanou částici o účinném průřezu σ(v). Budou to částice ze vzdálenosti (−tv, −(t + dt)v) od iontu z válce o objemu σ(v)|v|dt. Jejich množství tedy je Nvt dv dt = n(−tv)f (v)|v|σ(v)dv dt ,

(6.7)

kde n(x) je koncentrace a f (v)dv je rychlostní rozdělovací funkce částic. Předpokládáme homogenní rozdělení částic, a proto nahradíme n(x) konstantou n. Stejný výsledek dostaneme také, představíme-li si plyn s určitým rychlostním rozdělením jako superpozici svazků částic. Počet srážek iontu se svazkem částic o koncentraci f (v)n odvodíme ze vztahu (6.1) a získáme opět vztah (6.7) Pokud uvážíme případ částice pohybující se rychlostí v 0 napříč plynem, jehož rychlostní rozdělení v laboratorní soustavě je f0 (v), můžeme určit f (v) jako

42


f (v)dv = f0 (v+v0 )dv. Do vztahu (6.7) nyní dosadíme rychlostní rozdělení plynu a vztahem v lab = v + v 0 transformujeme tento vztah do laboratorní soustavy Nvt dv lab dt = nf0 (v lab )|vlab − v 0 |σ(v lab − v 0 )dv lab dt .

(6.8)

Tento vztah nám říká, jak často se budou částice o určité rychlosti srážet s naší sledovanou částicí, což je již veškerá informace potřebná pro modelování srážkových procesů. Celkovou srážkovou frekvenci Nt v závislosti na rychlosti částice můžeme určit integrací vztahu (6.8) přes celý rychlostní prostor a rychlost interagující částice potom bude výběrem z rozdělení Nv dv =

1 f0 (v)|v − v 0 |σ(v − v 0 )dv . Nt

(6.9)

Zde jsme provedli přeznačení v ≡ v lab , nebot’ další výpočty budou probíhat v laboratorní soustavě. Uvedený postup je sice správný, avšak pro použití v částicovém modelu nevhodný, nebot’ vyžaduje neustálé vyčíslování srážkové frekvence Nt a generování rychlosti částice výběrem z vícerozměrného parametrického rozdělení jež je programátorsky komplikované a výpočetně náročné. Abychom však potvrdili správnost tohoto modelu, aplikujeme jej na pohyb částice v plynu s maxwellovským rozdělením a ukážeme jednoduchý algoritmus, který již správně určuje energetické rozdělení částic. Maxwellovské rozdělení rychlostí je určeno vztahem r 2 1 2kB T v f0 (v)dv = 3 3/2 exp − 2 dv ; vm = . (6.10) vm π vm m Předpokládejme nyní, že náš problém je izotropní. Díky tomu můžeme vystředovat rychlostní rozdělení (6.8) přes úhlové souřadnice, aniž bychom, alespoň v případě modelu s interakcí přenosu náboje, ovlivnili výsledné rychlostní rozdělení. Výraz (6.8) nabývá po převedení do polárních souřadnic tvaru r v0 v2 (6.11) Nvθφt dv dθ dφdt = nf0 (v) 1 + 02 − 2 cos θ · v 3 σ sin θ dv dθ dφdt , v v kde bez újmy na obecnosti předpokládáme, že v 0 = (0, 0, v0 ) a účinný průřez jsme nahradili konstantou za užití metody nulové srážky. Po přeintegrování přes θ a φ dostáváme rozdělení v celkové rychlosti  4π v0 v  3  dv dt . nf0 (v)σv 3 +  v < v0 : 3 v0 v 2 Nvt dv dt = (6.12) 4π v0  3  nf0 (v)σv 3 + 2 dv dt .  v > v0 : 3 v

Pro rychlost částice v 0 = 0 nabývá toto rozdělení tvaru rozdělení maxwellovského toku, což bylo možné očekávat, jelikož na nehybnou částici musí dopadat právě

43


Maxwellovo rozdělení

Nt [s−1 ]

3e+07 2e+07 1e+07

srážková frekvence

λ [m]

0 6e-05 4e-05

očekávaná volná dráha 1/nσ vτ

2e-05 0

0

500

1000 v [m

1500

2000

· s−1 ]

Obrázek 6.2: Rychlostní závislost očekávané volné dráhy částice v plynu s maxwellovským rozdělením. Horní panel ukazuje rozdělení plynu, prostřední srážkovou frekvenci a dolní očekávanou volnou dráhu. Parametry přibližně odpovídají argonu při pokojové teplotě a tlaku 133 Pa.

maxwellovský tok. Na druhou stranu pro velké rychlosti v 0 se toto rozdělení redukuje na obyčejné maxwellovské rozdělení, což je také předvídatelný výsledek, nebot’ při velké rychlosti iontu přestává být termální rychlost částic plynu důležitá a problém se redukuje na dříve diskutovaný problém částice a nehybných rozptylových center. Integrováním tohoto rozdělení přes rychlost a přes jednotkový časový interval obdržíme počet částic dopadajících za jednotku času, tedy srážkovou frekvenci Nt respektive střední dobu mezi srážkami τ 2 1 vm v0 v02 1 vm 1 1+ . (6.13) exp − 2 + erf Nt = = nσv0 √ τ vm vm 2 v02 π v0 Ze srážkové frekvence můžeme určit očekávanou volnou dráhu vztahem −1 2 1 vm 1 vm v0 1 v02 v0 √ 1+ . = exp − 2 + erf λ(v0 ) = Nt nσ vm vm 2 v02 π v0

(6.14)

44


Závislost očekávané volné dráhy na rychlosti je vynesena v grafu 6.2 zároveň se srážkovou frekvencí a rychlostním rozdělením částic plynu. Z grafu je patrné, že při velkých rychlostech sledované částice je volná dráha konstantní a lze tedy simulovat srážky pomocí náhodné volné dráhy. To je případ srážek elektronů a iontů v nízkoteplotním plazmatu. Při velmi nízkých rychlostech je naopak konstantní doba života částice.

6.2.1

Referenční implementace

Nyní uvedeme referenční algoritmus implementující model srážek prezentovaný v předchozí sekci. Pravděpodobnost, že se částice během časového kroku ∆t alespoň jednou srazí určíme ze vztahu (6.3) jako ∆x v∆t ∆t pi = 1 − exp − = 1 − exp − = 1 − exp − . (6.15) λ(v) λ(v) τ (v) Zde již uvažujeme rychlostně závislou očekávanou volnou dráhu λ(v) dle vztahu (6.14) respektive očekávanou dobu života τ (v) dle (6.13). K simulaci srážky je dále potřeba náhodně vygenerovat druhou interagující částici. Rozdělení rychlostí interagujících částic je určeno vztahem (6.9). Jak jsme již uvedli, náš model je izotropní, a proto se výsledné rychlostní rozdělení nezmění, použijeme-li rozdělení přeintegrované přes úhlové souřadnice. Vystředované rozdělení rychlostí získáme tudíž normováním rozdělení (6.12) f (v)dv =

Nvt dv . Nt

(6.16)

Po dosazení ze vztahů (6.12) a (6.13) dostáváme −1 2 1 1 vm 1 vm v0 v02 √ f (v)dv = 1+ · exp − 2 + erf v0 vm vm 2 v02 π v0  v v  0 2  pro v < v0 .  · 3 + 4 v3 v v v0 · √ 3 exp − 2 dv v02  3 π vm vm  · 3 + pro v > v0 .  v2

(6.17)

Rychlosti částic z tohoto rozdělení jsou generovány von Neumannovou zamítací metodou. Postup spočívá v tom, že generujeme rovnoměrně rozložené body na ploše pod křivkou M(x), přičemž křivka M(x) je majorantou rozdělovací funkce f (x). V případě, že vygenerovaný bod leží pod křivkou f (x), je x-ová souřadnice bodu výběrem z rozdělení f (v). V opačném případě vygenerovaná čísla zamítneme a postup opakujeme. Jako majorantu rozdělení f (v) jsme zvolili distribuční funkci Rayleighova rozdělení ve tvaru 2 2b v , (6.18) M(v) = v exp − 2 a a2

45


f (v) [10−3 · s · m−1 ]

2.5

majoranta v0 = 1 v0 = 200 v0 = 300 v0 = 500 v0 = 1000 v0 = 2000

2 1.5 1 0.5 0

0

200

400

600

800

v [m ·

1000

1200

1400

s−1 ]

Obrázek 6.3: Rozdělení rychlosti interagujících neutrálních částic f (v)dv pro různé rychlosti iontu v0 a jejich majoranta (6.20).

kde parametr a ovlivňuje pouze „natažení“ funkce bez vlivu na plochu pod křivkou, zatímco parametr b přímo udává plochu pod křivkou, neboli Z ∞ M(v)dv = b . (6.19) 0

Důvodů pro tuto volbu je několik. Asymptotické chování této funkce pro velké hodnoty v je shodné s rozdělením f (v), zatímco pro malé v je pokles dokonce pomalejší. Navíc dokážeme jednoduše generovat body rovnoměrně rozložené pod touto křivkou. Hodnoty parametrů a a b jsme nalezli vizuálním porovnáním grafů rozdělení (6.17) pro různé hodnoty v0 , abychom dostali majorantu tvaru 2 · 1,52 v2 . (6.20) M(v) = v exp − 2 (vm · 1,55) (vm · 1,55)2 Ilustrace této křivky s několika tvary rozdělení (6.17) je zakreslena v grafu 6.3. Velikost parametru b = 1,52 nám říká, že pro vygenerování K náhodných čísel zamítací metodou musíme vygenerovat přibližně bK náhodných čísel. Body ležící rovnoměrně rozloženy pod majorantou (6.20) jsou generovány za užití metody inverzní distribuční funkce ze vztahů p ξ2 = M(ξ1 )ζ2 , (6.21) ξ1 = − ln(ζ1 ) · vm · 1,55 ;

kde ζ1 , ζ2 ∈ U(0, 1) jsou výběrem z rovnoměrného rozdělení na intervalu (0, 1). Generování zamítací metodou je ilustrováno v grafu 6.4.

46


f (v) [10−3 · s · m−1 ]

2.5

majoranta f (v)

2 1.5

zamítnuté body přijaté body

1 0.5 0

0

200

400

600

800

v [m ·

1000

1200

1400

s−1 ]

Obrázek 6.4: Ilustrace generování náhodných čísel z rozdělení f (v) von Neumannovou zamítací metodou.

Algoritmus jsme opět testovali výpočtem energetického rozdělení argonových iontů interagujících s neutrálními částicemi interakcí přenosu náboje. Vypočtené rozdělení je vykresleno v grafu 6.5. Z grafu je zřejmá dobrá shoda s teoretickým Maxwellovým rozdělením na rozdíl od výsledků metody s konstantní volnou dráhou zobrazených v grafu 6.1. Nevýhodou této metody je nutnost výpočtu pravděpodobnosti srážky v každém kroku, což lze eliminovat zavedením horního odhadu pro pravděpodobnost srážky Pmax a následným výpočtem pravděpodobnosti pouze u náhodně vybraných Pmax N částic. Další komplikací je, v případě že nás například zajímají úhlová rozdělení částic, nutnost generovat náhodné částice z poměrně složitého anizotropního rozdělení (6.8), což dále zpomaluje výpočet.

6.2.2

Alternativní implementace I

V literatuře lze najít různé modely, které se zaměřují na srážkové procesy mezi ionty a neutrálními částicemi, liší se však implementací algoritmu a také efektivitou výpočtu. Jeden z modelů iont-neutrálních srážek popisuje (Nanbu – Kitatani, 1995). Pro popis srážek se používá pravděpodobnost srážky v časovém intervalu ∆t daná vztahem P = ngσ(g)∆t , (6.22) kde g = |v 0 − v| je relativní rychlost, v 0 je rychlost iontu a v je náhodně vygenerovaná rychlost druhé neutrální částice z Maxwellova rozdělení. Pravděpodobnost, že se iont srazí s částicí o rychlosti v z intervalu (v, v + dv) je součin

47


f (v) [10−3 · s · m−1 ]

2

Maxwellovo rodělení rozdělení iontů

1,5 1 0,5 0 -1000

-500

0 v [m ·

500

1000

s−1 ]

Obrázek 6.5: Rychlostní rozdělení iontů v plazmatu vypočtené referenčním algoritmem zohledňujícím pohyb neutrálních částic.

pravděpodobnosti vygenerování takové částice a pravděpodobnosti srážky, tedy Pvt dv = P f0 (v)dv = nf0 (v)|v 0 − v|σ(|v0 − v|)dv∆t ,

(6.23)

což dává vztah podobný vztahu (6.8). Rychlostní rozdělení určíme normováním této pravděpodobnosti Pv = Kf0 (v)|v0 − v|σ(|v0 − v|)

(6.24)

a vidíme že rozdělení rychlosti interagujících částic je ve shodě se vztahem (6.9). Vypočteme nyní ještě střední hodnotu pravděpodobnosti srážky Z hP i = f0 (v)P dv . (6.25) Tento integrál je zcela analogický integrálu rozdělení (6.8), který je pro Maxwellovo rozdělení vyhodnocen ve vztahu (6.13). Proto můžeme psát hP i = Nt ∆t =

∆t . τ

(6.26)

Tato hodnota pravděpodobnosti se liší od veličiny v referenční implementaci (6.15). Jak ukážeme v kapitole 6.3, důsledkem použití takovéto pravděpodobnosti srážky může teoreticky být řádově větší chyba, ovšem v praxi nejsou rozdíly tak dramatické, jak ukazuje kapitola 7.2. Výraznou nevýhodou této metody je nutnost generovat částici z Maxwellova rozdělení a počítat relativní rychlost v každém časovém kroku, i když i zde by pomohl horní odhad pravděpodobnosti srážky.

48


6.2.3

Alternativní implementace II

Nyní odvodíme přesnou a relativně efektivní metodu simulace srážkových procesů. Vyjdeme ze základního vztahu (6.8), který zde pro připomenutí zopakujeme Nvt dv dt = nf0 (v)|v − v 0 |σ(v − v 0 )dv dt.

(6.8)

Tento výraz zjednodušíme tak, že nahradíme člen |v − v 0 |σ(v − v 0 ) konstantou. Onu konstantu určíme jako maximální hodnotou M = max{|v|σ(v)} , v

(6.27)

kde maximum počítáme přes všechny relativní rychlosti, které mohou v simulaci významnou měrou nastat. To odpovídá zavedení nového srážkového procesu s účinným průřezem M σ0 (v) = − σ(v) . (6.28) |v| Tento proces se nazývá nulovou srážkou a jeho zavedením zajistíme, že totální účinný průřez bude nepřímo úměrný relativní rychlosti částic. V případě ionneutrálních interakcí může taková situace nastat i přirozeně, nebot’ účinný průřez ion-neutrálních interakcí je dle Langevinovy teorie právě nepřímo úměrný rychlosti. Dostaneme tedy nové rozdělení Nvt dv dt = nf0 (v)M dv dt .

(6.29)

Rychlostní rozdělení interagujících atomů je nyní stejné jako rozdělení rychlostí všech atomů neutrálního pozadí nezávisle na rychlosti iontu, což značně zjednodušuje generování interagujících částic. Výpočet celkové srážkové frekvence integrací rozdělení (6.29) je nyní triviální, nebot’ rychlostní rozdělení neutrálního plynu f0 je normované: Z Nt dt = Nvt dv = nM dt . (6.30) Opět tedy dostáváme hodnotu nezávislou na rychlosti částice. Pravděpodobnost, že dojde k alespoň jedné srážce v časovém úseku ∆t, je potom dána vztahem ∆t . (6.31) P = 1 − exp(−Nt ∆t) = 1 − exp − τ

Tento model uvádí také (Vahedi – Surendra, 1995), avšak odvozuje jej za nesprávných předpokladů, jelikož zanedbává rychlost narážejících částic. V případě, že na základě vztahu (6.31) rozhodneme, že dochází k interakci, generujeme neutrální atom z Maxwellova rozdělení. Na základě relativní rychlosti částic v rel vypočteme účinné průřezy {σi }ni=1 jednotlivých interakcí respektive odpovídající srážkové frekvence νi = nvrel σi (vrel ). Potom s pravděpodobností úměrnou velikostem srážkových frekvencí vybereme jednu ze srážkových interakcí. Je pozoruhodné, že tato metoda v sobě ukrývá generování rozdělení (6.8).

49

6.3. VLIV DISKRÉTNÍHO ČASU NA SRÁŽKOVÉ PROCESY

model volná dráha (6.1) referenční (6.2.1) implementace II (6.2.3)

T 9s 33 s 22 s

N 5 · 106 9 · 106 16 · 106

N∅ 2000 3000 7 · 106

Nreal 5 · 106 9 · 106 9 · 106

T /N 1.6 µs 3.6 µs 1.3 µs

T /Nreal 1.6 µs 3.6 µs 2.3 µs

Tabulka 6.1: Srovnání efektivity různých algoritmů výpočtu srážek. Tabulka ukazuje výpočetní čas T simulace pohybu jedné částice po dobu 0.5 s pro různé algoritmy. Symboly N , N∅ a Nreal vyjadřují celkový počet srážek v simulaci, počet nulových srážek a počet nenulových srážek.

Výhodou této metody je kromě konstantního rozdělení interagujících částic také konstantní pravděpodobnost srážky pro všechny částice nezávisle na jejich rychlosti. To umožňuje efektivní výběr interagujících částic. Tato metoda je vlastně modifikací metody nulové srážky zavedené v kapitole 6.1. Dříve uvedená metoda zajišt’uje pomocí nulové srážky, aby celkový účinný průřez a tím i očekávaná volná dráha byly konstantní. Toho však nelze dosáhnout, pokud rychlosti částic plynu nejsou zanedbatelné. Na druhou stranu tato modifikovaná metoda zajišt’uje pomocí nulové srážky, aby srážková frekvence byla konstantní. Tento postup již funguje nezávisle na rychlostním rozdělení částic plynu. Pomocí jednorozměrného neselfkonzistentního programu, který byl použit k výpočtu grafů 6.1 a 6.5 jsme provedli ilustrační srovnání rychlosti výpočtu různými algoritmy. Srovnali jsme klasický výpočet pomocí generování náhodné volné dráhy 6.1, dále referenční implementaci 6.2.1 a nakonec algoritmus 6.2.3. Simulovali jsme dynamiku jednoho kladného iontu v plynu o tlaku 133 Pa po dobu 0,5 s. Potenciál byl předem určen pomocí jednorozměrného selfkonzistentního modelu s rovinnou symetrií popsaného v (Roučka, 2006). Výpočty byly provedeny na procesoru Intelr CoreTM 2 T5500 @ 1.66GHz bez paralelizace. Tabulka ukazuje, že „nový“ algoritmus je přibližně dvakrát pomalejší, než algoritmus generující náhodnou volnou dráhu. Tento rozdíl však není způsoben pomalostí nového algoritmu, nýbrž hrubou nepřesností původního algoritmu, který zanedbává srážky pomalých částic. Uvedené výsledky jsou pouze ilustrační, nebot’ jsou vypočteny za použití dnes již zastaralé sady účinných průřezů, které jsou velmi blízké konstantě a tím dále zvýhodňují algoritmus s náhodnou volnou dráhou.

6.3

Vliv diskrétního času na srážkové procesy

Jak bylo uvedeno v úvodní kapitole, modely particle-in-cell pracují s diskrétním časovým krokem. Délka časového kroku je dána několika faktory. Především musí být časový krok kratší, než je typická časová škála zkoumaných jevů, což může být například perioda oscilací v plazmatu a pod. Obecněji lze říci, že časový krok volíme tak, aby chyba řešení pohybových rovnic byla pro nás malá. Další


50

podmínkou je, že časový krok musí být menší než střední doba mezi srážkami. V praxi obvykle požadujeme, aby časový krok byl mnohem menší než uvedené časové škály. Uvedená kritéria jsou však pouze orientační a vhodný časový krok obvykle najdeme tak, že provedeme simulaci s různými délkami časového kroku, a potom vybereme nejdelší časový krok, při kterém ještě nedochází ke změně výsledků. Problém simulací při vyšších tlacích však je, že doba mezi srážkami je velmi malá, což nás nutí k nadměrnému zkracování časového kroku simulace. Pro urychlení výpočtu je výhodné volit časový krok blízký srážkové frekvenci částic. Jaké jsou však důsledky takové volby pro energetické a rychlostní rozdělení částic? To si nyní ukážeme na zjednodušeném modelu srážkových procesů. V následující sekci rozumíme pod pojmem PIC model s diskrétním časem, vlivy elektromagnetické interakce částic do naší diskuse nezahrnujeme. Pravděpodobnost, že v časovém kroku ∆t dojde k alespoň jedné srážce je dána obvyklým vztahem ∆t P1 = 1 − exp − , (6.32) τ

kde τ označuje střední dobu mezi srážkami. V článku (Vahedi – Surendra, 1995) bylo odvozeno množství zanedbaných srážek následující úvahou: Předpokládáme, že srážkové procesy jsou na sobě nezávislé, proto pravděpodobnost, že proběhne n srážek v čase ∆t, je určena součinem pravděpodobností jednotlivých srážek P1n . Množství vynechaných srážek v jednom kroku je dáno vztahem r′ =

∞ X n=2

P1n =

P12 . 1 − P1

(6.33)

Zkusíme-li však tento vztah vyhodnotit například pro ∆t = τ , dostaneme r = 1,09, což je očividně špatný výsledek, nebot’ střední počet srážek je 1 a nemůžeme zanedbat více srážek, než je jejich celkové množství. Chyba spočívá v tom, že s událostí srážky nemůžeme při výpočtu zacházet jako s hodnotou náhodné proměnné a tak je nesmyslné počítat pravděpodobnost dvou srážek podle vztahu pro podmíněnou pravděpodobnost dvou nezávislých proměnných. Ve skutečnosti podléhá množství srážek v jednom časovém kroku Poissonovu rozdělení a pravděpodobnost n srážek v intervalu ∆t je dána vztahem n ∆t 1 ∆t (6.34) exp − n! τ τ Množství zanedbaných srážek bychom tedy mohli získat sečtením řady definované na základě této pravděpodobnosti srážky. Existuje však podstatně jednodušší cesta ke správnému výsledku: Když střední doba mezi srážkami je rovna τ , přihodí se za jednotku času v průměru 1/τ srážek a za čas ∆t tedy dochází v průměru


51

k ∆t/τ srážkám. Rozdíl této hodnoty oproti pravděpodobnosti srážky nám dává počet zanedbaných srážek ∆t ∆t ∆t r= . (6.35) − P1 = − 1 + exp − τ τ τ

Podstatné je však spíše relativní množství vynechaných srážek −1 ∆t ∆t 1 − exp − −1 r/P1 = τ τ

(6.36)

Pokud bychom chtěli zachovat celkový počet srážek, mohli bychom použít výraz pro pravděpodobnost ∆t P2 = , (6.37) τ jenž je použit například v metodě (Nanbu – Kitatani, 1995) prezentované v sekci 6.2.2. Tento vztah dle teorie v následující sekci nedává dobré výsledky při výpočtu energie, avšak výsledky numerických simulací jsou nejednoznačné.

6.3.1

Střední hodnoty veličin v modelu PIC

Uvažujme nyní plyn, ve kterém je srážková frekvence sledovaných částic nezávislá na rychlosti. Dochází pouze k jedinému druhu srážky, jenž se projevuje zastavením částice. Na částice přitom působí konstantní homogenní síla. Tento model je sice pouze ilustrační idealizací, avšak podobné vlastnosti vykazují ionty při srážkách s neutrálním plynem. Nyní vypočteme driftovou rychlost částice v takovém plynu. Dobu od poslední srážky označíme T . Potom rychlost v této době vypočteme jako v(T ) = aT ,

(6.38)

kde a označuje zrychlení vlivem vnější síly. Celkovou střední rychlost dále určíme středováním přes rozdělení doby od poslední srážky. To je určeno vztahem T 1 dT . (6.39) f (T )dT = exp − τ τ

Pro střední hodnotu dostáváme vztah Z ∞ hvi = v(T )f (T )dT ,

(6.40)

0

který vypočteme integrací per partes po dosazení z rovnic (6.38) a (6.39) hvi = aτ .

(6.41)

v 2 (T ) = a2 T 2 ,

(6.42)

2 v = 2a2 τ 2 .

(6.43)

Obdobně budeme postupovat i pro kvadrát rychlosti, neboli energii. Výpočet je analogický, a proto uvedu pouze výsledky

52


Rychlost Nyní se podíváme, jaké budou střední hodnoty těchto veličin v modelu particle-in-cell. Zásadní problém je, jak vůbec definovat střední hodnotu veličiny. Vezmeme-li si totiž velký statistický soubor ve stacionárním stavu, zjistíme, že střední hodnoty jeho parametrů jsou závislé na „čase“. Při každém posunu o jeden časový krok ∆t totiž získávají energii od vnějšího pole, načež dojde ke skokové ztrátě energie systému vlivem srážek. Vypočteme tedy střední hodnoty parametrů částic v systému středované přes časový úsek ∆t. Dále také střední hodnoty parametrů částic v okamžiku před výpočtem srážek metodou MCC a v okamžiku po výpočtu srážek MCC. Nejprve definujeme pravděpodobnost Q = 1 − P , že částice se za čas ∆t ani jednou nesrazí. Nyní vypočteme střední rychlost hvii částice, která před i časovými kroky utrpěla srážku a od té doby nerušeně letí Z (i+1)∆t Z (i+1)∆t 1 1 1 hvii = v(t)dt = atdt = a∆t(2i + 1) . (6.44) ∆t i∆t ∆t i∆t 2 Dále vypočteme střední hodnotu hvii přes soubor částic. Pravděpodobnost, že se částice v posledním kroku srazila, je rovna P . Pravděpodobnost, že se srazila částice v dřívějším kroku a poté se šířila bez srážky, je rovna P Q atd. Střední hodnota tedy bude dána řadou hviPIC = P hvi0 + QP hvi1 + Q2 P hvi2 · · · ,

(6.45)

kterou můžeme převést do tvaru hviPIC = P

∞ X i=0

Qi hvii .

(6.46)

Dosazením ze vztahu (6.44) dostaneme hviPIC = P a∆t

∞ X i=0

∞

∞

X 1 1X i Qi (i + ) = P a∆t Q Q Qi−1 i + 2 2 i=0 i=1

!

.

(6.47)

První sumu v tomto výrazu převedeme na geometrickou řadu a vysčítáme ∞ X i=1

Qi−1 i =

∞ ∞ X d i d X i 1 d Q Q = = , Q = 2 dQ dQ dQ 1 − Q (1 − Q) i=1 i=1

(6.48)

zatímco druhá řada je přímo geometrickou řadou a její součet je roven 1/P . Po dosazení dostáváme střední rychlost částic hviPIC =

a∆t (Q + P/2) . P

(6.49)

Závislost relativní odchylky oproti driftové rychlosti v modelovému systému (6.41), definované jako | hvi−hviPIC |/ hvi, je zakreslena v grafu 6.6 pro obě volby pravděpodobnosti srážky.

53


relativní chyba hvi

100

P1 před MCC P1 po MCC P1 střední P2 střední

10−1 10−2 10−3 10−4 10−5

5

10

15

20

25

30

35

40

45

50

τ /∆t Obrázek 6.6: Relativní chyba výpočtu hvi v závislosti na délce časového kroku.

Chceme-li nyní určit střední rychlost částic v okamžiku těsně před výpočtem srážek MCC, stačí v předchozím výpočtu nahradit střední rychlost hvii okamžitou rychlostí před výpočtem srážek v((i+1)∆t) = a(i+1)∆t. Postup výpočtu je zcela analogický a výsledkem je a∆t . (6.50) hvipřed PIC = P Zde stojí za povšimnutí, že při použití pravděpodobnosti srážky P2 je střední rychlost přesně shodná s modelovým systémem. Rychlost po srážce lze dát do souvislosti s rychlostí před srážkou vztahem před hvipo PIC = Q hviPIC =

a∆t Q. P

(6.51)

Relativní odchylka rychlosti od modelového systému je opět zakreslena v obrázku 6.6. Z grafu je patrné, že při použití pravděpodobnosti P2 je střední hodnota rychlosti zatížena řádově větší chybou. Stojí také za povšimnutí chyba rychlosti v okamžiku před srážkou, která je relativně velká. Při časovém kroku 5× kratším než doba mezi srážkami dosahuje chyba přibližně 10 %, což může znatelně ovlivnit následky srážkového procesu. Energie Dále zjistíme, jaké budou hodnoty energie respektive kvadrátu rychlosti v modelu PIC. Stejně jako v případě rychlosti částic nejprve vypočteme střední hodnotu přes jeden časový krok

2 1 v i= ∆t

Z

(i+1)∆t

i∆t

1 v (t)dt = ∆t 2

Z

(i+1)∆t

i∆t

1 a2 t2 dt = a2 ∆t2 (i2 + i + ) . 3

(6.52)

6.4. MODEL SRÁŽEK NEZÁVISLÝ NA VOLBĚ ČASOVÉHO KROKU

54

Střední hodnotu hv 2 iPIC nyní vypočteme z analogie vztahu (6.46) jako # "∞ ∞ ∞ ∞ X X X X

2

1 v PIC = P Qi + . (6.53) Qi v 2 i = P a2 ∆t2 Qi i2 + Qi i + 3 i=0 i=0 i=0 i=0 Řady v předchozím vztahu opět vyhodnotíme převodem na geometrickou řadu a dostáváme výsledek

2 Q 1 2 2 Q(1 + Q) v PIC = a ∆t . (6.54) + + P2 P 3 Relativní odchylka modelu PIC oproti fyzikálnímu modelu (6.43) je zakreslena v grafu 6.7. Obdobně jako v případě rychlosti ještě určíme střední hodnotu kvadrátu rychlosti před výpočtem srážek a po výpočtu srážek

2 před 1+Q v PIC = a2 ∆t2 , P2

(6.55)

2 po Q(1 + Q) v PIC = a2 ∆t2 . (6.56) P2 Výsledky jsou taktéž zakresleny v grafu 6.7 a jsou podobné jako pro rychlost. Odchylka energie v okamžiku před srážkou je ještě větší než u rychlosti. Požadujemeli odchylku energie částice vstupující do srážky menší než 10 %, musíme volit časový krok kratší než přibližně desetina střední doby mezi srážkami. Dokonce i když zvolíme časový krok 50× kratší, bude odchylka stále přibližně 2 %. To je však již zanedbatelná chyba v porovnání s přesností vstupních dat.

6.4

Model srážek nezávislý na volbě časového kroku

Výsledky předchozí kapitoly ukazují, že model srážek Monte Carlo vyžaduje, aby délka časového kroku byla několikanásobně kratší než doba mezi srážkami. Délka časového kroku v modelu při vyšším tlaku tak bude shora omezena právě požadavky Monte Carlo simulace. Model při tlaku řádově 1000 Pa tak vyžaduje délku časového kroku kolem 10−13 s. Při tak krátkém časovém kroku jsou výpočty extrémně časově náročné, jak uvidíme v kapitole 7. Toto omezení jsme eliminovali modifikováním metody particle-in-cell. Vytvořili jsme metodu, která počítá srážky nezávisle na časovém kroku a čas je tedy z hlediska srážek spojitý. Klasický model PIC integruje běžnými numerickými metodami Newtonovy pohybové rovnice a potom v diskrétních okamžicích testuje, jestli došlo ke srážce. Naše modifikace zachovává obvyklé schéma modelu

55


100

relativní chyba v 2

10−1 P1 před MCC P1 po MCC P1 průměr P2 průměr

10−2 10−3 10−4 5

10

15

20

25

30

35

40

45

50

τ /∆t

Obrázek 6.7: Relativní chyba výpočtu v 2 pro různé pravděpodobnosti srážky v závislosti na délce časového kroku.

PIC znázorněné v kapitole 3.6, ovšem nahrazuje výpočet pohybových rovnic zobecněným výpočtem dynamiky částic, jenž zahrnuje také vliv srážek. Výpočet srážek na konci časového kroku potom již není potřeba. Ona zobecněná dynamika částic je vlastně řešení transportního problému metodou Monte Carlo. Postupujeme tak, že každé částici přiřadíme nový parametr, čas do srážky. Na počátku výpočtu přiřadíme každé částici náhodně vygenerovaný čas do srážky výběrem z exponenciálního rozdělení se střední hodnotou rovnou střední době mezi srážkami definované v kap. 6.2. V případě, že čas do srážky určité částice je větší než délka časového kroku, počítáme trajektorii částice klasicky přes čas ∆t a odečteme ∆t od času do srážky. V opačném případě počítáme trajektorii z Newtonových rovnic do času do srážky. Vypočteme důsledky srážky, opět vygenerujeme nový čas do srážky a tento cyklus opakujeme dokud nedosáhneme času ∆t. Zjednodušeně můžeme říci, že srážky probíhají v náhodných časech daných srážkovou frekvencí zcela nezávisle na diskrétním časovém kroku metody particle-in-cell. Tento princip je také znázorněn v obrázku 6.8. Obrázek ukazuje, jak srážky v klasickém algoritmu probíhají pouze v násobcích ∆t, zatímco v novém algoritmu mohou probíhat kdykoliv. Tato modifikace algoritmu má také důsledky pro metodu integrace pohybových rovnic. Časový krok integrace je totiž určen náhodným výběrem z exponenciálního rozdělení a navíc v okamžiku srážky potřebujeme znát rychlost i polohu zároveň, což se vylučuje s použitím algoritmu leap-frog. Proto musíme hledat alternativní metodu integrace pohybových rovnic. Nejjednodušší variantou je modelování pohybu mezi srážkami jednodušším Eulerovým algoritmem s přesností prvního řádu. Vzhledem k tomu, že dominantní vliv na dynamiku částic mají


klasický

∆t

∆t

∆t

∆t

∆t

τ3

τ4

56

∆t

nový τ1

τ2

Obrázek 6.8: Ilustrace rozdílu mezi novým a klasickým algoritmem simulace srážek. Horní šipka znázorňuje časovou osu historie částice v klasickém algoritmu, zatímco dolní časová osa představuje nový algoritmus. Čtverečky označují srážky.

právě srážkové procesy, není přesnost integrace pohybových rovnic tak důležitá a Eulerův algoritmus by mohl pro potřeby našeho modelu postačovat. Jako další variantu jsme zkoušeli modifikaci Verletovy metody pro potřeby našeho algoritmu. Protože pohyb částice simulujeme až v několika náhodných krocích během jednoho časového kroku, neznáme změnu působící síly během jednoho kroku a rychlostní Verletova metoda není v běžné formě použitelná. Tento problém jsme obešli předpokladem konstantní síly v průběhu jednoho časového kroku a Verletův algoritmus jsme modifikovali do tvaru xk+1 = xk + v k ∆t + v k+1 = v k +

Fk 2 ∆t m

(6.57)

Fk ∆t . m

Jaké jsou výsledky této metody integrace pohybových rovnic ukážeme v kapitole 7. Srovnání numerických výsledků tohoto modelu s klasickým modelem typu particle-in-cell pro různé konfigurace experimentu a především pro různé délky časového kroku také uvedeme v kapitole 7. Pro zjednodušení zápisu budeme v dalším textu a především v grafech označovat tento algoritmus jako advance2, zatímco klasický algoritmus typu particlein-cell a model srážek s pravděpodobností srážky P1 označíme jako advance1. Do některých výpočtů ještě zahrneme model srážek s pravděpodobností srážky P2 označený jako advance0.

Kapitola 7 Výsledky modelu Na základě teoretických výsledků sekcí 6.2.3 a 6.4 byl implementován model typu particle-in-cell, jehož základní charakteristiky byly popsány v kapitole 3.6. Některé detaily implementace upřesníme v následujících odstavcích. Jak již bylo uvedeno, jde o dvourozměrný model s kartézskou mříží. Okrajové podmínky odpovídají simulaci neporušeného plazmatu. Implementaci okrajových podmínek znázorňuje obrázek 7.1. Obdélníková pracovní oblast je obložena zdrojovými oblastmi, což jsou oblasti, ve kterých probíhá neselfkonzistentní simulace objemu plazmatu s periodickými okrajovými podmínkami (PBC) za daného vnějšího elektrického pole. Vždy, když některá částice zdroje prolétne rozhraním s pracovní oblastí, dochází k její duplikaci, protože jednak pokračuje v letu do pracovní oblasti a jednak se v souladu s periodickými okrajovými podmínkami přesouvá na opačnou stranu zdrojové oblasti. Vzhledem k tomu, že jsme experimentovali s délkou časového kroku, byl kladen důraz na správnou funkci zdroje nezávisle na délce časového kroku. Použitý zdroj tak pracuje bezchybně, i když zdrojové částice během jednoho časového kroku prolétnou mnohokrát napříč zdrojovou oblastí. Výhodou použitých okrajových podmínek je, že na rozdíl od analyticky definovaného generátoru toku částic zde nezadáváme do modelu žádný předpoklad o tvaru rozdělovací funkce elektronů ani o jejich teplotě. Rozdělení přilétajících elektronů je v okrajové oblasti vypočteno konzistentně s rozdělením v pracovní oblasti. Model srážkových procesů odpovídá teoretickým závěrům sekce 6.2.3. Program umožňuje simulaci multikomponentního plazmatu s mnoha různými interakcemi. Účinné průřezy σi jsou do programu zadávány buď v podobě textových tabulek, nebo jako funkční předpis, ten je však nutné zadat přímo do zdrojového kódu. V inicializační fázi programu jsou numericky nalezeny hodnoty max{σv} a na základě koncentrací interagujících složek je vypočtena střední doba mezi srážkami pro každou složku. Vzhledem k tomu, že modelujeme v oblasti vysokých tlaků, jsme se při provádění optimalizací zaměřili právě na výpočet srážek, který zabírá podstatnou část výpočetního času. 57

58

Obrázek 7.1: Schéma dvourozměrného modelu a jeho okrajových podmínek. optimalizace žádné Ziggurat Ziggurat + tabelování σ

doba výpočtu [s] 155,1 72,9 28,0

doba výpočtu [%] 100 47 18

Tabulka 7.1: Ukázka vlivu optimalizací Monte Carlo simulace na celkovou dobu simulace. Časy udávají dobu simulace pohybu argonových iontů v plynu o tlaku 1330 Pa s časovým krokem 5 · 10−8 s pomocí modelu srážek nezávislého na časovém kroku. Za daných podmínek docházelo přibližně ke dvaceti srážkám na jeden časový krok, což zvýrazňuje časovou náročnost MCC simulace.

Přínos těchto optimalizací je velmi významný, jak ukazuje tabulka 7.1. Původní program používal generátor náhodných čísel s rovnoměrným rozdělením rand() ze standardní knihovny jazyka C. Náhodná doba života byla generována výběrem z exponenciálního rozdělením metodou inverzní distribuční funkce a náhodné rychlosti částic jsme generovali výběrem z normálního rozdělení s užitím vztahu (3.78). Druhý řádek tabulky ukazuje více než dvojnásobné urychlení výpočtu po implementování generátoru normálního a exponenciálního rozdělení typu Ziggurat (Marsaglia – Tsang, 2000). Pro generování rovnoměrně rozdělené náhodné proměnné byl použit generátor s posuvnými registry shr3() (Marsaglia – Tsang, 2000) popsaný v sekci 3.6.5. Použití generátoru shr3() v kombinaci s metodou Ziggurat bylo podrobeno kritice v práci (Nadler, 2006). Odhalená odchylka od normálního rozdělení dosahuje však v nejhorším případě řádu pouze 10−4 , což je ve srovnání s celkovou přesností našeho modelu a především jeho vstupních dat zanedbatelné. Přesto jsme pro jistotu provedli srovnávací test s mírně pomalej-

7.1. VERIFIKACE MODELU MCC

59

ším generátorem mwc(), který ukázal, že výsledky modelu jsou pro oba generátory v rámci statistické chyby shodné. Dalšího urychlení jsme dosáhli optimalizací výpočtu účinných průřezů, který je nutné provést v každém okamžiku srážky. V původním modelu byly použity aproximace účinných průřezů dle (Phelps). Výpočet takových účinných průřezů vyžadoval několikanásobně vyhodnocení mocninné funkce, což je poměrně náročná výpočetní operace. Proto jsme do inicializační fáze programu přidali automatické tabelování funkcemi zadaných účinných průřezů. Jak ukazuje třetí řádek tabulky 7.1, bylo touto úpravou dosaženo dalšího více než dvojnásobného urychlení. Časové testy v této tabulce byly měřeny na notebooku s procesorem Intelr CoreTM 2 T5500 @ 1.66GHz s pamětí 3 GB DDR2 667 MHz bez paralelizace, stejně jako všechny další měření v této kapitole. Pro potřeby měření bylo využito funkce getrusage() ze systémové knihovny <sys/resource.h> jazyka C. Nadcházející sekce uvedou některé numerické výsledky našeho modelu a především srovnání našeho nového algoritmu s klasickým modelem particle-in-cell. Veškeré výpočty jsou, pokud není uvedeno jinak, provedeny při intenzitě vnějšího elektrického pole 500 V · m−1 . Ve všech výpočtech bylo využitu triku pro urychlení výpočtu, který spočívá v tom, že časový krok iontů volíme delší (konkrétně 1000×) než časový krok iontů. Čas iontů tedy plyne v simulaci tisíckrát rychleji. V případě neselfkonzistentního modelu je tento postup zcela legitimní, nebot’ částice spolu neinteragují. V případě selfkonzistentního modelu nelze touto metodou zkoumat dynamické jevy, avšak v případě stacionárního stavu je tato metoda použitelná. Pokud tedy v následujícím textu nebude časový krok iontů explicitně uveden, znamená to, že je 1000× delší než krok elektronů.

7.1

Verifikace modelu MCC

V kapitole 6.2 v grafu 6.5 jsme již ověřili výsledky MCC modelu pro ionty argonu. Nyní ověříme také výsledky při výpočtu energetického rozdělení elektronů. Výpočty se týkají objemu plazmatu bez přítomnosti sondy a byly určeny neselfkonzistentním modelem. Jako modelový systém jsme použili jednak čisté argonové plazma a jednak čisté kyslíkové plazma. V případě argonového plazmatu jsme brali v úvahu pružnou srážku, excitaci a ionizaci. V kyslíkovém plazmatu jsme uvažovali všechny interakce uvedené v tabulce 5.1. Jako referenční výsledek nám posloužily výsledky programu ELENDIF (Morgan – Penetrante, 1990). V grafech 7.2 a 7.3 jsou zakresleny výsledky metody Monte Carlo a programu ELENDIF, který počítá rozdělovací funkce metodou rozvoje Boltzmannovy rovnice do sférických funkcí. Oba grafy vykazují výbornou shodu obou metod.

60

7.1. VERIFIKACE MODELU MCC

100

ELENDIF MC

10−1 f (E) [eV−1 ]

10−2 10−3 10−4 10−5 10−6 10−7 10−8

0

5

10

15

20

E [eV] Obrázek 7.2: Srovnání výpočtu EEDF metodou Monte Carlo s výsledkem programu ELENDIF při poli E/n = 15 Td v čistém argonu.

100

ELENDIF MC

f (E) [eV−1 ]

10−1 10−2 10−3 10−4 10−5 10−6 10−7

0

2

4

6

8

10

E [eV] Obrázek 7.3: Srovnání výpočtu EEDF metodou Monte Carlo s výsledkem programu ELENDIF při poli E/n = 15 Td v čistém kyslíku.

7.2. SROVNÁNÍ NESELFKONZISTENTNÍCH MODELŮ

7.2

61

Srovnání neselfkonzistentních modelů

Nyní ukážeme, jaké jsou výsledky modelu srážek se spojitým časem ve srovnání s modelem o diskrétním časovém kroku. Demonstrovat budeme na výpočtu střední energie iontů a elektronů. Modelový systém byl tvořen objemem nízkoteplotního argonového plazmatu bez přítomnosti sondy. Simulace byla provedena neselfkonzistentním výpočtem v konstantním elektrickém poli, abychom izolovali vliv simulace srážek od ostatních efektů spojených s diskretizací času, jako je například chyba numerického řešení pohybových rovnic a proměnlivost selfkonzistentního elektrického pole. Chybu metody určujeme vždy vůči určitému výsledku, který je dle dostupných informací nejsprávnější. Tento výsledek je ve většině srovnání odhadnut z grafů jako hodnota, ke které konvergují všechny metody s klesajícím časovým krokem.

7.2.1

Tlak 266 Pa

Graf 7.4 (a) ukazuje výborné výsledky nového algoritmu advance2, nebot’ vypočtená hodnota střední energie elektronů je v rámci přesnosti nezávislá na volbě časového kroku. Na druhou stranu chyba klasických algoritmů dosahuje při časovém kroku ∆t rovném střední době mezi srážkami τ již několika procent. Pokud chceme udržet chybu menší než přibližně 1 %, je vhodné volit časový krok 10−11 s. Výsledky získané pro energii iontů jsou zobrazeny v grafu 7.4 (b). Nový algoritmus opět vykazuje nezávislost vypočtené střední energie na volbě časového kroku. Chyba klasických algoritmů pro ∆t = τ dosahuje hodnoty přibližně 1 %. Pokud tedy chceme omezit chybu výpočtu přibližně na 1 %, zvolíme časový krok 10−8 s. Ukazuje se tedy, že použití 1000× delšího časového kroku pro ionty v neselfkonzistentním modelu je oprávněné. Pokud bychom ovšem chtěli zjišt’ovat například frekvenci výskytu určité interakce, musíme vzít v úvahu vztah 6.36, udávající počet zanedbaných srážek. V případě iontů, kde časový krok 10−8 s přibližně odpovídá jejich střední době mezi srážkami se tak dopouštíme zanedbání téměř 40 % srážek. Je tedy zajímavé, že i přes takovou odchylku v množství srážek je chyba v určení energie poměrně malá. Požadujeme-li chybu v určení počtu srážek taktéž přibližně 1 %, musíme dle vztahu 6.36 volit časový krok pro ionty přibližně 2 · 10−9 s. Případně bychom mohli použít pravděpodobnost srážky P2 , která dovoluje prodloužení kroku až na τ bez zanedbání srážek. Graf 7.4 (c) ukazuje, jak dlouho trvá výpočet dynamiky souboru 200 000 elektronů a 200 000 iontů v časovém rozmezí jedné sekundy času elektronů a 1000 sekund času iontů. Graf ukazuje, že klasický algoritmus je pro dlouhé časové kroky podstatně rychlejší než nový algoritmus, což je však způsobeno podhodnocením počtu srážek. Pro krátké časové kroky jsou rychlosti algoritmů téměř stejné, což dokazuje, že režijní náklady nového algoritmu jsou malé. Konkrétně při časovém kroku 5 · 10−13 s je nový algoritmus pouze o 1 % pomalejší než ten klasický. Jelikož nový algoritmus advance2 dává výsledky nezávislé na časovém kroku,

62


5,55 5,5

advance0 advance1 advance2

τ

E [eV]

5,45 5,4 5,35 5,3 5,25

10−12

10−11

10−10

∆telektron [s] 40

E [meV]

39,8

advance0 advance1 advance2

(a)

τ

39,6 39,4 39,2 39

10−9

10−8

10−7

doba výpočtu [1010 s]

∆tion [s] 5 4,5 4 3,5 3 2,5 2 1,5 1 0,5 0

(b) advance1 advance2

10−12

10−11 ∆telektron [s]

10−10 (c)

Obrázek 7.4: Výsledky neselfkonzistentního modelu při tlaku 266 Pa v závislosti na délce časového kroku. Grafy (a) a (b) ukazují střední energii elektronů a iontů. Graf (c) znázorňuje dobu výpočtu jedné sekundy elektronového času.


63

můžeme zvolit takový krok tak, aby výpočet byl co nejrychlejší. Při časovém kroku 5 · 10−11 s tak dostáváme dobu výpočtu přibližně 1,15 · 1010 s. Časový krok klasického modelu musíme volit s ohledem na chyby zanesené do výpočtu srážek. Pokud bychom se spokojili s 1% chybou ve střední hodnotě energie, můžeme zvolit časový krok 10−11 s pro elektrony a 10−8 s pro ionty. Výpočet s tímto časovým krokem trvá 0,89 · 1010 s. V tomto případě je tedy klasický algoritmus rychlejší, ovšem za cenu podhodnocení počtu srážek. V předchozích odstavcích jsme však odvodili, že když požadujeme 1% chybu v celkovém počtu srážek, musíme volit časový krok 2 · 10−9 s pro ionty a odpovídající 2 · 10−12 s pro elektrony. Doba výpočtu klasickým algoritmem potom dosahuje 2,16 · 1010 s. Nový algoritmus je tedy v tomto případě téměř dvakrát rychlejší a jeho použití se nepochybně vyplatí. Do grafů byl zahrnut také klasický algoritmus s netradiční pravděpodobností srážky P2 . Výsledky takového algoritmu jsou v případě elektronů překvapivě lepší než při použití pravděpodobnosti P1 . Tento výsledek může být námětem pro další práci v podobě podrobného porovnání těchto dvou variant klasického algoritmu.

7.2.2

Tlak 1330 Pa

Nyní se podíváme na analogické výsledky, ovšem při vyšším tlaku, konkrétně 1330 Pa. Graf 7.5 (a) ukazuje, že při časovém kroku 10−11 s je již chyba výpočtu neúnosná a dosahuje více než 20 %. Z grafu lze odečíst, že pro dosažení chyby řádově 1 % je potřeba počítat s časovým krokem 5 · 10−13 s až 10−12 s. Bohužel je z grafu patrné, že křivky se ani při nejkratším zvoleném kroku nescházejí. To si vysvětluji příliš pomalou konvergencí k ustálenému stavu. Soubor elektronů zřejmě nebyl ani po uplynutí 3 µs ve zcela stacionárním stavu. Výrazně odlišné výsledky vykazuje graf 7.5 (b), který znázorňuje střední energii iontů. Chyba v energii je zanedbatelná již při časovém kroku 2 · 10−9 s, který odpovídá střední době mezi srážkami. To je zřejmě způsobeno povahou interakcí iontů s neutrálními atomy. Pružná srážka a především resonanční přenos náboje totiž zprostředkovávají velmi efektivně výměnu energie mezi ionty a atomy, a proto se střední energie iontů udržuje stabilně blízko střední energie atomů. Střední energie tedy pravděpodobně není dobrým ukazatelem přesnosti výpočtu. Vyvstává zde ovšem otázka, jak posoudit přesnost výpočtu. Zde je například odchylka od správné hodnoty energie velmi malá i pro relativně velký časový krok 2 · 10−9 s. Přitom víme, že při tomto časovém kroku dochází k zanedbání více než 30 % srážek. Odpovědí pravděpodobně bude, že časový krok musíme volit s ohledem na požadované výsledky modelu. Zvolíme-li však dlouhý časový krok, musíme mít stále na paměti že zanedbání v oblasti srážek se může nakonec projevit neočekávaným způsobem, a tento postup tedy nelze doporučit. V grafu 7.5 (c) je vynesena závislost doby výpočtu na časovém kroku obdobně jako v grafu 7.4 (c). Výpočet novým algoritmem trvá na použitém hardwaru 5,5 · 1010 s při ∆t = 5 · 10−11 s. Doba výpočtu klasickým algoritmem při časovém kroku elektronů 10−12 s je 6,4·1010 s. Délka časového kroku diktovaná požadavkem

64


5,5

E [eV]

5

advance1 advance2

τ

4,5 4 3,5 3 10−13

10−12

10−11

10−10

10−8

10−7

∆telektron [s] 40 39,8 39,6

advance1 advance2

τ

E [meV]

39,4 39,2 39 38,8 38,6 10−10

10−9 ∆tion [s]


25

advance1 advance2

20 15 10 5 0 10−13

10−12

10−11

10−10

∆telektron [s] Obrázek 7.5: Výsledky neselfkonzistentního modelu při tlaku 1330 Pa v závislosti na délce časového kroku. Grafy (a) a (b) ukazují střední energii elektronů a iontů. Graf (c) znázorňuje dobu výpočtu jedné sekundy elektronového času.

7.3. SROVNÁNÍ SELFKONZISTENTNÍCH MODELŮ

65

na maximálně 1% zanedbání srážek dokonce klesá až k 3 · 10−13 s. Odpovídající doba výpočtu potom dosahuje 14·1010 s, což je již výrazně delší čas oproti novému algoritmu.

7.3

Srovnání selfkonzistentních modelů

Výsledky nového algoritmu v případě neselfkonzistentního algoritmu jsou sice pozitivní, avšak větší význam má implementace nového algoritmu do selfkonzistentního particle-in-cell modelu. Výsledkům nového algoritmu bude věnována tato sekce. Základem našich testů byl výše popsaný dvourozměrný model typu particlein-cell. Plazma je argonové s koncentrací iontů n = 1015 m−3 . Vnější elektrické pole E = 500 V·m−1 . Simulace zahrnovala 2·105 iontů a 2·105 elektronů, rozměry pracovní oblasti byly 1 cm × 1 cm a rozměry mříže byly 200 × 200.

7.3.1

Tlak 266 Pa

První test spočívá v simulaci nekonečné válcové sondy o průměru 0,1 mm, na kterou je přivedeno napětí −5 V oproti potenciálu plazmatu. Tlak plazmatu byl roven 266 Pa. Protože se v okolí sondy vyskytuje vysoce nehomogenní elektrické pole, je tato konfigurace vhodná také k testování integrátoru pohybových rovnic. V případě klasického PIC modelu využíváme vždy leap-frog metodu (3.77). Nový algoritmus není kompatibilní s leap-frog metodou, a proto jsme pro integraci zkoušeli klasický Eulerův algoritmus (3.75) a modifikovaný Verletův algoritmus (6.58). Graf 7.6 (a) ukazuje závislost střední energie iontů dopadajících na sondu, což je důležitý parametr při popisu interakce plazmatu s pevnou látkou. Nový algoritmus s Verletovým integrátorem pohybových rovnic dává patrně nejlepší výsledky, ačkoliv se jeho chyba příliš neliší od varianty s Eulerovým algoritmem. Ovšem výsledky klasického algoritmu jsou zde podstatně horší. Zdá se, že při kroku ∆t = 1 · 10−12 s je odchylka stále ještě kolem 5 %. Takové chyby dosahuje nový algoritmus až při časovém kroku 1 · 10−11 s. Můžeme tedy volit až 10× delší časový krok. Rozdíly v grafu 7.6 (b), který zobrazuje proud iontů na sondu, jsou ještě daleko markantnější. Zde jsou výsledky nového algoritmu s Verletovým integrátorem prakticky nezávislé na volbě časového kroku, odchylka při extrémním časovém kroku 5 · 10−11 s činí pouhé 1 %. Výsledky s Eulerovým integrátorem jsou podstatně horší, při časovém kroku ∆t = 5 · 10−12 s dosahuje odchylka téměř 10 % a pro omezení chyby pod 5 % by byl potřeba časový krok 2 · 10−12 s. Nejhorší výsledky však vykazuje klasický algoritmus, kde již při časovém kroku 2 · 10−12 s je chyba přibližně 9 %. Je tedy nutné volit krok 1 · 10−12 s, pokud požadujeme chybu maximálně 5 %. V případě měření toku iontů na sondu tedy nový algoritmus za daných podmínek prakticky neklade podmínky na časový krok na rozdíl


66

od klasického algoritmu. Výsledky výpočtu proudu elektronů na sondu v grafu 7.6 (c) jsou poněkud odlišné, nebot’ zde nejhorší chování vykazuje nový algoritmus s Eulerovým integrátorem. Tento algoritmus nadále nebudeme diskutovat, nebot’ jde na základě všech testů o horší variantu nového algoritmu. V dalším textu tedy pod novým algoritmem rozumíme jeho variantu s modifikovaným Verletovým integrátorem. Výsledky nového algoritmu při výpočtu elektronového proudu jsou zde prakticky shodné s klasickým algoritmem. Pro dosažení 5% chyby musíme volit časový krok 5 · 10−12 s až 1 · 10−11 s. Závěrem všech srovnání tedy je, že pokud chceme, aby chyba sledovaných parametrů nepřesahovala 5 % musíme u klasického PIC algoritmu volit časový krok přibližně 1 · 10−12 s, na druhou stranu nový algoritmus dovoluje použití časového kroku 5 · 10−12 s nebo dokonce až 1 · 10−11 s. Jak ukazuje graf 7.7, je doba výpočtu nový algoritmem při uvedených volbách časového kroku 4× až 5× kratší.

7.3.2

Tlak 665 Pa

Další obdobnou sadu simulací jsme provedli při poněkud vyšším tlaku 665 Pa. Sonda měla tentokrát průměr 0,2 mm a je na ni přivedeno napětí −3 V. Výsledky výpočtů jsou zpracovány zcela analogickým způsobem v grafu 7.9. Graf 7.9 (a) ukazuje střední energii iontů dopadajících na sondu. Nový algoritmus vykazuje až do časového kroku 2 · 10−12 s chybu pod 1 % a při časovém kroku 5 · 10−12 s chybu 5 %. Klasický algoritmus dosahuje chyby 5 % již při časovém kroku 2 · 10−12 s. Výsledky jednotlivých algoritmů se v tomto grafu příliš neliší. Graf 7.9 (b) ukazuje tok iontů na sondu a zde jsou rozdíly opět nejvýraznější. Odchylky nového algoritmu zde v celém rozsahu časových kroků nepřesahují 1 %. Odchylka klasického algoritmu přitom přesahuje 5 % již při časovém kroku 5 · 10−13 s a pro ∆t = 1 · 10−12 s dosahuje již 13 %. Tok elektronů zobrazený v grafu 7.9 nevykazuje signifikantní závislost na volbě časového kroku pro ∆t ≤ 2 · 10−11 s. Pokud bychom tedy opět chtěli udržet chybu všech sledovaných parametrů pod 5 %, můžeme pro klasický algoritmus zvolit časový krok maximálně 5·10−13 s z důvodu velké chyby při výpočtu iontového proudu. Nový algoritmus přitom umožňuje časový krok o délce až 5 · 10−12 s. To podle grafu 7.8 znamená až sedminásobné urychlení výpočtu.

67


0,7

E [eV]

0,6

τion

0,5 0,4

advance1 advance2 Euler advance2 Verlet

0,3 10−13

10−12

10−11

10−10

∆telektron [s] 0,8

I [µA]

0,7 0,6

(a)

advance1 advance2 Euler advance2 Verlet τion

0,5 0,4 0,3 0,2 10−13

10−12

10−11

10−10

∆telektron [s]

(b)

-20

I [µA]

-25 -30

τe

-35 -40

advance1 advance2 Euler advance2 Verlet

-45 10−13

10−12


10−10 (c)

Obrázek 7.6: Výsledky selfkonzistentního algoritmu při tlaku 266 Pa, předpětí na sondě -5 V a vnějším elektrickém poli 500 V · m−1 pro různé délky časového kroku. Graf (a) ukazuje střední energii iontů dopadajících na sondu. Grafy (b) a (c) ukazují proud iontů a elektronů na sondu.

68



60

advance1 advance2 Verlet

50 40 30 20 10 0 10−13

10−12

10−11

10−10

∆telektron [s] Obrázek 7.7: Časová náročnost selfkonzistentního algoritmu při tlaku 266 Pa, předpětí na sondě -5 V a vnějším elektrickém poli 500 V · m−1 pro různé délky časového kroku. Graf ukazuje dobu výpočtu jedné sekundy elektronového času při simulaci se 400 000 částicemi a mříží 200 × 200.


60


50 40 30 20 10 0 10−13

10−12

10−11

10−10

∆telektron [s] Obrázek 7.8: Časová náročnost selfkonzistentního algoritmu při tlaku 665 Pa, předpětí na sondě -3 V a vnějším elektrickém poli 500 V · m−1 pro různé délky časového kroku. Graf ukazuje dobu výpočtu jedné sekundy elektronového času při simulaci se 400 000 částicemi a mříží 200 × 200.

69


0,22 0,2

advance1 advance2 Verlet τion

E [eV]

0,18 0,16 0,14 0,12 0,1 10−13

10−12

10−11

10−10

∆telektron [s] 0,2 0,18

(a)


I [µA]

0,16 τion

0,14 0,12 0,1 0,08 10−13

10−12

10−11

10−10

I [µA]

∆telektron [s] -40 -42 -44 -46 -48 -50 -52 -54 -56 advance1 -58 advance2 Verlet -60 10−13 10−12

(b)

τe


10−10 (c)

Obrázek 7.9: Výsledky selfkonzistentního algoritmu při tlaku 665 Pa, předpětí na sondě -3 V a vnějším elektrickém poli 500 V · m−1 pro různé délky časového kroku. Graf (a) ukazuje střední energii iontů dopadajících na sondu. Grafy (b) a (c) ukazují proud iontů a elektronů na sondu.

7.4. SONDOVÉ CHARAKTERISTIKY

7.4

70

Sondové charakteristiky

Na závěr uvedeme některé fyzikální výsledky našeho modelu válcové Langmuirovy sondy. Pro porovnání modelu s experimentem je jednoznačně nejvýhodnější modelovat sondové charakteristiky, nebot’ ty lze získat přímo z měření, a proto zde uvedeme sondové charakteristiky pro různé konfigurace systému. Veškeré výpočty jsme provedli při tlaku 266 Pa, vnějším poli 500 V · m−1 a průměru sondy 0,2 mm. Simulace počítá s nekonečnou válcovou sondou, výsledky jsou ovšem normovány na délku sondy 1 cm, což je řádově typická délka Langmuirovy sondy. Sondové charakteristiky jsme počítali pomocí sekvence simulací. Nejprve byl pouze pro urychlení dalších výpočtů určen ustálený stav plazmatu v přítomnosti sondy na určitém víceméně libovolném potenciálu. Ten sloužil jako počáteční stav pro další simulace. Při výpočtu proudu za určitého napětí jsme nejprve nechali systém po určitou dobu ustalovat a poté jsme statistickým středováním určili hodnotu sondového proudu a jiných parametrů plazmatu. Nejprve byla vypočtena sondová charakteristika v čistém argonovém plazmatu. Výsledek je zakreslen v grafu 7.10. Simulovaná charakteristika vykazuje typické znaky, jako je velmi slabý iontový proud pro záporná napětí na sondě, dále plovoucí potenciál mezi −10 V a −15 V a nakonec exponenciální nárůst elektronového proudu až do potenciálu plazmatu, kde se nachází inflexní bod. Dále v grafech 7.11 a 7.12 ukazujeme obdobnou charakteristiku v kyslíkovém plazmatu. Kyslíkové plazma má obvykle komplikované složení, ve kterém figurují různé druhy kladných i záporných iontů. Správné hodnoty koncentrací pro simulaci by bylo možné do jisté míry zjistit experimentálně nebo pomocí kinetického modelu. Jak jsme již uvedli, kinetický model prezentovaný v kapitole 5 není dost spolehlivý na kvantitativní výpočet koncentrací. Byl ale nastaven tak, aby alespoň pořadí prvků podle koncentrací odpovídalo výsledkům z literatury. Výsledky modelu při poli 7 Td udávají, že nabité částice s nejvyšší koncentrací jsou − ionty O+ 2 a elektrony, následovány zápornými ionty O . Přesné koncentrace nám však nejsou známy, a tak jsme pouze pro ilustraci zvolili dva extrémní případy. Nejprve jsme uvažovali kyslíkové plazma, ve kterém se kromě neutrálních částic vyskytují pouze kladné ionty O+ 2 a elektrony. Model srážek iontů s neutrálními částicemi byl inspirován modelem (Vahedi – Surendra, 1995) a uvažuje pouze interakci přenosu náboje. V kyslíkovém plazmatu dochází k většímu množství interakcí, kterými elektrony ztrácejí energii, a proto je zde i nižší plovoucí potenciál, jehož hodnota nabývá přibližně −5 V. Výsledná charakteristika je zakreslena v grafu 7.11. Nakonec je v grafu 7.12 uvedena sondová charakteristika elektronegativního − a žádné elektrony. kyslíkového plazmatu, které obsahuje pouze ionty O+ 2 a O Pro záporné ionty jsme uvažovali interakci pružné srážky s molekulami o účinném průřezu dle (Vahedi – Surendra, 1995). Je patrné, že elektronová větev sondové charakteristiky zde zanikla a byla nahrazena tokem záporných iontů, který je řádově slabší. Asymetrie pozorované charakteristiky je způsobena jednak nižší hmotností záporných iontů a jednak jinou povahou jejich interakce s plynem.

71

7.4. SONDOVÉ CHARAKTERISTIKY

-0,3

I [mA]

-0,2 -0,1 0 -30

-25

-20

-15

-10

-5

0

5

10

15

U [V]

I [mA]

Obrázek 7.10: Voltampérová charakteristika v argonu při tlaku 266 Pa.

-0,12 -0,1 -0,08 -0,06 -0,04 -0,02 0 0,02 -10

-5

0

5

10

15

U [V]

I [µA]

Obrázek 7.11: Voltampérová charakteristika v kyslíkovém plazmatu obsahujícím ionty O+ 2 a elektrony při tlaku 266 Pa.

-0,07 -0,06 -0,05 -0,04 -0,03 -0,02 -0,01 0 0,01 0,02

-3

-2

-1

0

1

2

3

U [V] Obrázek 7.12: Voltampérová charakteristika v elektronegativním kyslíkovém − plazmatu obsahujícím ionty O+ 2 a O a žádné elektrony při tlaku 266 Pa.

7.5. OSTATNÍ VÝSLEDKY

7.5

72

Ostatní výsledky

Výhodou částicového modelu je, že kromě makroskopických parametrů plazmatu nám umožňuje zjištění prakticky libovolných mikroskopických veličin. Nyní si ukážeme některé výsledky modelu, jež by bylo velmi obtížné zkoumat experimentálními technikami. Následující řadu grafů jsme získali simulací ustáleného stavu argonového plazmatu při tlaku 266 Pa. Na vloženou sondu o poloměru 0,2 mm bylo přivedeno napětí −3 V. Vnější pole bylo opět rovno 500 V · m−1 . Výpočet byl proveden při časovém kroku 5 · 10−12 s pro elektrony a 5 · 10−9 s pro ionty. Prvním z výsledků jsou energetická rozdělení částic v simulaci. Graf 7.13 ukazuje energetické rozdělení elektronů dopadajících na sondu ve srovnání s rozdělením v objemu. Z grafu je patrné, že vlivem brzdného pole sondy v kombinaci se srážkami ve stínící vrstvě je energie při dopadu na sondu menší než u elektronů v objemu plazmatu. Podotýkáme, že tento jev není samozřejmý, nebot’ přinejmenším pro elektrony s maxwellovským rozdělením v bezesrážkovém plazmatu je jejich jejich energetické rozdělení na povrchu sondy nezávislé na velikosti brzdného potenciálu, mění se pouze jejich množství v souladu s MaxwellovýmBoltzmannovým rozdělením (Roučka, 2006, s. 38). Obdobné výsledky pro ionty argonu uvádíme v grafu 7.14. Zde je na rozdíl od elektronů patrný silný vliv urychlujícího pole sondy. Ionty dopadající na sondu tak dosahují několikanásobně vyšších energií než ionty v objemu plazmatu. Graf 7.15 ukazuje úhlové rozdělení elektronů dopadajících na povrch sondy. Vypočtené úhlové rozdělení se dobře shoduje s kosinovým rozdělením, které je v grafu také zakresleno. Tento výsledek nám napovídá, že úhlové rozdělení elektronů je v blízkosti sondy přibližně izotropní. Úhlové rozdělení iontů při dopadu na sondu je vyneseno v grafu 7.16. Zde je patrná významná odchylka od kosinového rozdělení směrem ke kolmému dopadu. To je nepochybně způsobeno přitažlivým působením sondy. Další veličinou, jejíž měření je velmi obtížné až nemožné je elektrostatický potenciál v okolí sondy. Z výpočtu metodou particle-in-cell jej přitom můžeme určit velmi jednoduše a výsledek ukazuje graf 7.17. Nakonec uvedeme ještě dva grafy veličin, které současně s potenciálem tvoří skupinu veličin definovaných na mříži v modelu particle-in-cell. Jde o koncentraci elektronů v grafu 7.18 a koncentraci iontů v grafu 7.19. Koncentrace elektronů vykazuje pokles na jistém okolí sondy, což je způsobeno jejím odpudivým potenciálem a jde o projev mechanismu stínění. Na druhou stranu koncentrace iontů klesá jen pomalu až k samému povrchu sondy. Mírný pokles koncentrace obou složek na velké vzdálenosti od sondy je způsoben „odsáváním“ těchto částic sondou.

73


100

f (E) [eV−1 ]

10−1 10−2 10−3 10−4

na sondě v objemu 0

2

4

6

8

10

12

14

E [eV] Obrázek 7.13: Energetické rozdělení elektronů dopadajících na sondu a elektronů v objemu plazmatu.

f (E) [eV−1 ]

na sondě v objemu 101

100

0

0,1

0,2

0,3

0,4

0,5

E [eV] Obrázek 7.14: Energetické rozdělení iontů dopadajících na sondu a iontů v objemu plazmatu.

74


0

0,1

0,2

0,3

0,4

0,5

rozdělení elektronů kosinové rozdělení Obrázek 7.15: Úhlové rozdělení elektronů dopadajících na sondu ve srovnání s kosinovým rozdělením.

0

0,2

0,4

0,6

0,8

rozdělení iontů kosinové rozdělení Obrázek 7.16: Úhlové rozdělení iontů dopadajících na sondu ve srovnání s kosinovým rozdělením.

75


1

U [V]

-3

0

-2

-1 -2

-1

-3

0

7

1 7

5 y [mm] x [mm]

5

3 3

Obrázek 7.17: Elektrický potenciál v okolí válcové sondy na potenciálu −3 V při vnějším poli 500 V · m−1 .

7 8 · 1014 6 · 1014 5

4 · 1014

4

2 · 1014

3

0 · 100

3

4

5

6

7

x [mm]

Obrázek 7.18: Koncentrace elektronů v okolí válcové sondy.

n [m−3 ]

y [mm]

6

76


7 8 · 1014 6 · 1014 5

4 · 1014

4

2 · 1014

3

0 · 100

3

4

5

6

7

x [mm]

Obrázek 7.19: Koncentrace iontů v okolí válcové sondy.

n [m−3 ]

y [mm]

6

Kapitola 8 Závěr V oblasti fyziky plazmatu při středních tlacích a z našeho pohledu především v oblasti modelování plazmatu při středních tlacích existuje mnoho zatím nevyřešených problémů. Práce se zabývá částicovými modely a jejich hlavní nectností jsou extrémní výpočetní nároky. Na tento problém se, mimo jiné, předložená diplomová práce zaměřuje a pomocí nové modifikace algoritmu particle-in-cell ukazuje možnost několikanásobného zrychlení výpočtů. Kapitola 5 prezentuje model chemické kinetiky kyslíkového plazmatu. Sestavení věrohodného modelu bez důsledné interakce s experimentem se však ukázalo jako obtížné, a proto jsou výsledky vytvořeného modelu platné jen kvalitativně. V sekci 6.2 byla provedena analýza algoritmů pro simulaci srážek metodou Monte Carlo. Byl potvrzen již dříve známý fakt, že pro simulace srážek nelze obecně použít algoritmus založený na náhodné volné dráze. V sekci 6.2.1 byl navržen referenční algoritmus na základě teoretických závěrů sekce 6.2 a bylo ověřeno, že jeho výsledky jsou ve shodě s teorií. U dvou dalších algoritmů vyskytujících se v literatuře jsme matematicky ověřili, že jsou s tím naším ekvivalentní a pro další použití v našem modelu jsme vybrali ten nejefektivnější. Dále jsme v sekci 6.3 analyticky rozebrali různé chyby vzniklé diskretizací času při výpočtu srážek a ukázali jsme, že je vhodné volit časový krok simulace několikrát kratší než střední doba mezi srážkami. Tento omezující požadavek na délku časového kroku jsme eliminovali novou metodou simulace srážek v kapitole 6.4. Nová metoda je kompatibilní s algoritmem particle-in-cell a spočívá v simulaci srážek nezávisle na časovém kroku simulace. Nakonec jsme v kapitole 7 provedli numerické testy našich algoritmů. Verifikovali jsme náš Monte Carlo model pomocí porovnání s alternativní metodou v sekci 7.1. Především jsme však porovnali výsledky nového algoritmu výpočtu srážek oproti klasickému algoritmu se srážkami v diskrétním čase. V případě neselfkonzistentního výpočtu v homogenním konstantním poli jsou výsledky nového algoritmu zcela nezávislé na volbě časového kroku a získané urychlení výpočtu je maximálně dvojnásobné. Stěžejní část výsledků však spočívá v selfkonzistentní 77

78 particle-in-cell simulaci válcové sondy v sekci 7.3. Výpočty ukazují, že ve většině případů je chyba nového algoritmu podstatně menší a při požadavku na stejnou velikost chyby je výpočet novým algoritmem čtyřikrát až sedmkrát rychlejší, což je příznivý výsledek. Nová metoda nám tedy umožňuje modelovat v plazmatu o několikanásobně vyšším tlaku při zachování výpočetních nároků. Tato metoda neprovádí žádná zanedbání v oblasti výpočtu srážek, což je na jednu stranu její výhoda, nebot’ můžeme modelu více důvěřovat, avšak na druhou stranu se to může jevit i jako nevýhoda při přechodu do ještě vyšších tlaků. Množství srážek při vyšších tlacích totiž roste a přímo úměrně tomu se i zpomaluje výpočet částicovou simulací. Pokud chceme rozšířit možnosti částicových modelů do kvalitativně nové oblasti tlaků kolem 10 kPa nebo dokonce až k atmosférickému tlaku, musíme výpočet dynamiky částic dále modifikovat. Možností je například využití Langevinovy dynamiky nebo jiných efektivních metod pro započtení vícenásobných srážek jako například (Iza – Lee, 2006). Na závěr jsme provedli několik simulací s fyzikálními výsledky. Jedním z výsledků jsou sondové charakteristiky v argonovém a kyslíkovém plazmatu. Pomocí sondových charakteristik by byla možná přímá konfrontace modelu s experimentem. Aparatura pro měření při středních tlacích v kyslíkovém a argonovém plazmatu na KFPP sice bohužel nebyla zprovozněna před dokončením práce, avšak měření sondových charakteristik je v plánu během nejbližších dnů. Hlavní výsledky diplomové práce budou prezentovány v září 2008 na mezinárodní konferenci JVC (Joint Vacuum Conference, Balatonalmádi, 22.–26. 9. 2008), kde již pravděpodobně budou doplněny o vazbu na experiment.

Příloha A Popis programu Zdrojové kódy programu jsou volně dostupné na adrese . Dvourozměrný model plazmatu je na webu publikován pod názvem plasma2d. Pro kompilaci je nutné mít kompilátor jazyka C++ (testováno s g++ 4.2.3, 4.1.3), knihovny UMFPACK (Davis, 2004a) a GotoBLAS (Goto) pro rychlé řešení Poissonovy rovnice a dále knihovnu getpot a knihovnu boost-filesystem . Volitelně je možné také nainstalovat program gnuplot pro průběžné vykreslování výsledků při výpočtu. Veškerá konfigurace programu se provádí pomocí konfiguračního souboru. Ke zdrojovému kódu je přiložen ukázkový konfigurační soubor, který slouží zároveň jako dokumentace. Program se spouští z příkazové řádky a jako parametry lze specifikovat konfigurační soubor a výstupní adresář ve tvaru > plasma2d config=config_0.txt output_dir=argon/dt/advance1/2e-13s Program je psán v jazyce C++ s využitím objektově orientovaného programování. Třída Species implementuje veškeré metody potřebné pro modelování určité složky plazmatu. Pro přidání nové složky do plazmatu tak stačí vytvořit novou třídu jako potomka třídy Species. Tato třída si pouze musí uchovávat informaci o základních fyzikálních parametrech modelované složky a musí implementovat metodu scatter(), která se volá v případě srážky. Po zařazení této třídy do globální tabulky složek je nová složka plnohodnotnou součástí modelu.

79

Literatura BARNES, J. – HUT, P. A hierarchical O(N log N) force-calculation algorithm. Nature. 1986, 324, 6096, s. 446–449. BARTOŠ, P. Hybridní modelování ve fyzice plazmatu. Dizertační práce, MFF UK, 2007. BOX, G. – MULLER, M. A note on the generation of random normal deviates. Annals of Mathematical Statistics. 1958, 29, 2, s. 610–611. BROK, W. J. M. The effect of dust on the argon metastable density in an RFdischarge. Trainee report VDF/NG 98-07, Eindhoven University of Technology, Eindhoven, The Netherlands, 1998. CAPITELLI, M. – DILONARDO, M. – GORSE, C. Self consistent electron energy distribution functions in non-equilibrium oxygen. In XIVth International Conference on Phenomena in Ionized Gases, Grenoble, France, 1979. CARTWRIGHT, K. L. – VERBONCOEUR, J. P. – BIRDSALL, C. K. Nonlinear Hybrid Boltzmann-PIC Acceleration Scheme. 1998. CHANTRY, P. J. A simple formula for diffusion calculations involving wall reflection and low density. Journal of Applied Physics. 1987, 62, 4, s. 1141–1148. CHEN, F. F. Úvod do fyziky plazmatu. Praha : Academia, 1984. DAVIS, T. A. Algorithm 832: UMFPACK, an unsymmetric-pattern multifrontal method. ACM Transactions on Mathematical Software. 2004a, 30, 2, s. 196– 199. DAVIS, T. A. A column pre-ordering strategy for the unsymmetric-pattern multifrontal method. ACM Transactions on Mathematical Software. 2004b, 30, 2, s. 165–195. DAVIS, T. A. – DUFF, I. S. A combined unifrontal/multifrontal method for unsymmetric sparse matrices. ACM Transactions on Mathematical Software. 1999, 25, 1, s. 1–20.

80

LITERATURA

81

DEBYE, P. – HÜCKEL, E. Zur Theorie der Elektrolyte. I. Gefrierpunktserniedrigung und verwandte Erscheinungen. Physikalische Zeitschrift. 1923, 24, 9, s. 185–206. ELIASSON, B. Basic data for modelling of electrical discharge in gasses: oxygen. Baden : Brown Boveri Forschungszentrum, 1986. ENTLICHER, M. Nové algoritmy částicového modelování metodou molekulární dynamiky. Dizertační práce, MFF UK, 2004. FROST, L. S. – PHELPS, A. V. Rotational excitation and momentum transfer cross sections for electrons in H2 and N2 from transport coefficients. Physical Review. 1962, 127, s. 1621–1633. GOTO, K. GotoBLAS. Dostupné z: resources/software/#blas>.

GREENGARD, L. – ROKHLIN, V. A fast algorithm for particle simulations. Journal of Computational Physics. 1987, 73, 2, s. 325–348. GROSSMANNOVÁ, H. – CIGÁNEK, M. – KRČMA, F. High-molecular products analysis of VOC destruction in atmospheric pressure discharge. Journal of Physics: Conference Series. 2007, 63, 1, s. 012011. HAYS, G. N. – TRACY, C. J. – OSKAM, H. J. Surface catalytic efficiency of a sputtered molybdenum layer on quartz and Pyrex for the recombination of nitrogen atoms. The Journal of Chemical Physics. 1974, 60, 5, s. 2027–2034. HOCKNEY, R. W. – EASTWOOD, J. W. Computer Simulation Using Particles. Bristol : Adam Hilger, 1988. HOLSTEIN, T. Energy distribution of electrons in high frequency gas discharges. Physical Review. 1946, 70, s. 367–384. IONIN, A. A. et al. Physics and engineering of singlet delta oxygen production in low-temperature plasma. Journal of Physics D: Applied Physics. 2007, 5, s. R25–R61. IZA, F. – LEE, J. K. Moment conserving method for modelling multiple collisions in particle simulations. Journal of Physics D: Applied Physics. 2006, 39, 9, s. 1853–1865. KOLEV, I. – BOGAERTS, A. Numerical Models of the Planar Magnetron Glow Discharges. Contributions to Plasma Physics. 2004, 44, 7–8, s. 582–588. KOSSYI, I. A. et al. Kinetic scheme of the non-enquilibrium discharge in nitrogenoxygen mixtures. Plasma Sources Science and Technology. 1992, 1, s. 207–220.

LITERATURA

82

KRACÍK, J. – ŠESTÁK, B. – AUBRECHT, L. Základy klasické a kvantové fyziky plazmatu. Praha : Academia, 1974. LANDAU, D. P. – BINDER, K. A Guide to Monte Carlo Simulations in Statistical Physics. Cambridge : Cambridge university press, second edition, 2003. LANGMUIR, I. Oscillations in ionized gases. Proceedings of the National academy of Sciences of the United States of America. 1928, 14, s. 627–637. LAROUSSI, M. – LU, X. Room-temperature atmospheric pressure plasma plume for biomedical applications. Applied Physics Letters. 2005, 87, 11, s. 113902. MARSAGLIA, G. Diehard Battery of Tests of Randomness, 1995. Dostupné z: . MARSAGLIA, G. Random numbers for C: End, at last?, 1999. Dostupné z: . MARSAGLIA, G. – TSANG, W. W. The Ziggurat Method for Generating Random Variables. Journal of Statistical Software. 10 2000, 5, 8, s. 1–7. MARTIŠOVITŠ, V. Základy fyziky plazmy. Bratislava : Fakulta matematiky, fyziky a informatiky Univerzita Komenského, 2004. MORGAN, W. L. – PENETRANTE, B. M. ELENDIF: A time-dependent Boltzmann solver for partially ionized plasmas. Computer Physics Communications. 1990, 58, s. 127–152. MUELLER, T. O. – COWAN, J. – SWANSON, E. The Use of Oxygen in SEM Plasma Cleaning Equipment. Microscopy and Microanalysis. 2007, 13, s. 210– 211. NADLER, B. Design Flaws in the Implementation of the Ziggurat and Monty Python methods (and some remarks on Matlab randn), 2006. Dostupné z: . NANBU, K. – KITATANI, Y. An ion-neutral species collision model for particle simulation of glow discharge. Journal of Physics D: Applied Physics. 1995, 28, s. 324–330. ODEHNAL, P. Studium interakce chemicky aktivního plazmatu s povrchy pevných látek. Diplomová práce, 1996. PEKÁREK, Z. – LAHUTA, M. – HRACH, R. Improving performance of multidimensional Particle-In-Cell codes for modelling of medium pressure plasma. Journal of Physics: Conference Series. 2007, 63, 1, s. 012009.

LITERATURA

83

PHELPS, A. V. COMPILATION OF ELECTRON CROSS SECTIONS. Dostupné z: . PRESS, W. H. et al. Numerical Recipes in C. Cambridge : Cambridge university press, second edition, 2002. RALSTON, A. Základy numerické matematiky. Praha : Academia, 1973. REKTORYS, K. et al. Přehled užité matematiky II. Praha : Prometheus, 2000. ROUČKA, Š. Počítačové modelování interakce nízkoteplotního plazmatu s pevnými látkami, 2006. bakalářská práce, MFF UK. SOMMERER, T. J. – KUSHNER, M. J. Numerical investigation of the kinetics and chemistry of rf glow discharge plasmas sustained in He, N2, 02, e/N2/02, He/CF4/02, and SiH4/NH3 using a Monte Carlo-fluid hybrid model. Journal of Applied Physics. 1992, 71, 4, s. 1654–1673. STOFFELS, E. et al. Negative ions in a radio-frequency oxygen plasma. Physical Review E. 1995, 51, 3. SUGAWARA, M. – STANSFIELD, B. L. Plasma Etching: Fundamentals and Applications. : Oxford University Press, 1998. ŠIMEK, J. Rozvoj metod počítačové fyziky pro fyziku plazmatu a fyziku tenkých vrstev. Dizertační práce, MFF UK, 2006. VAHEDI, V. – SURENDRA, M. A Monte Carlo collision model for the particlein-cell method: applications to argon and oxygen discharges. Computer Physics Communications. 1995, 87, s. 179–198. VASILJEVA, A. N. et al. On the possibility of O2 (a1 ∆g ) production by a nonself-sustained discharge for oxygen–iodine laser pumping. Journal of Physics D: Applied Physics. 2004, 37, s. 2455–2468.

Univerzita Karlova v Praze Matematicko-fyzikální fakulta DIPLOMOVÁ PRÁCE. Studium interakce plazma-pevná látka při středních tlacích

Recommend Documents