Univerzita Karlova v Praze Matematicko-fyzikální fakulta

Univerzita Karlova v Praze Matematicko-fyzikáln´ı fakulta

´ PRACE ´ DIPLOMOVA Vojtˇech Hrubý Studium interakce plazma-pevn´ a l´ atka pomoc´ı hybridn´ıho modelov´ an´ı Katedra fyziky povrch˚ u a plazmatu Vedouc´ı diplomové práce: Prof. RNDr. Rudolf Hrach, DrSc. Studijn´ı program: Fyzika Studijn´ı obor: Fyzika povrch˚ u a ionizovaných prostˇred´ı

2009

Dˇekuji panu prof. RNDr. Rudolfu Hrachovi, DrSc. za obˇetavou pomoc a vynikaj´ıc´ı spolupráci bˇehem mého studia a pˇr´ıpravy diplomové práce. Pan´ı doc. RNDr. Vˇeˇre Hrachové, CSc. dˇekuji za ˇcetné konzultace a za výbornou výuku fyziky plazmatu. Panu Mgr. Lukáˇsi Schmiedtovi dˇekuji za pomoc pˇri mˇeˇren´ı. Dˇekuji své rodinˇe a pˇrátel˚ um za duchovn´ı i materiáln´ı pomoc, bez které bych tuto práci nezvládl. Pan´ı MUDr. Eliˇsce Rouˇckové dˇekuji za návrat mezi zdravé.

Prohlaˇsuji, ˇze jsem svou diplomovou práci napsal samostatnˇe a výhradnˇe s pouˇzit´ım citovaných pramen˚ u. Souhlas´ım se zap˚ ujˇcován´ım práce. V Praze dne 17. dubna 2009

Vojtˇech Hrubý

Obsah ´ 1 Uvod

6

2 Z´ akladn´ı poznatky fyziky plazmatu 2.1 Definice plazmatu . . . . . . . . . . . . . . . 2.2 Klasifikace plazmatu . . . . . . . . . . . . . 2.3 Kinetika plazmatu . . . . . . . . . . . . . . 2.4 Boltzmannova rovnice . . . . . . . . . . . . 2.5 Diagnostika n´ızkoteplotn´ıho plazmatu . . . . 2.6 Interakce plazmatu s povrchy pevných látek

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

7 . 7 . 8 . 9 . 11 . 12 . 14

3 Poˇ c´ıtaˇ cov´ a fyzika

17

4 Hybridn´ı modelov´ an´ı 4.1 Vztahy metod . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ˇ asticové modelován´ı . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 C´ 4.3 Spojité modelován´ı . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.1 Numerické metody ˇreˇsen´ı soustav parciáln´ıch diferenciáln´ıch rovnic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Obvyklé hybridn´ı modely . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Energetický model . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Prostorový model . . . . . . . . . . . . . . . . . . . . . . . 4.4.3 Iteraˇcn´ı model . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.4 Dalˇs´ı hybridn´ı modely . . . . . . . . . . . . . . . . . . . .

20 . 20 . 21 . 26 . . . . . .

5 C´ıle diplomov´ e pr´ ace 6 Srovn´ an´ı vybran´ ych hybridn´ıch 6.1 Energetický model . . . . . . 6.2 Prostorový model . . . . . . . 6.3 Iteraˇcn´ı model . . . . . . . . . 6.4 Kritéria srovnán´ı . . . . . . .

28 32 32 33 35 36 37

model˚ u . . . . . . . . . . . . . . . . . . . . . . . .

3

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

38 38 40 41 42

4

OBSAH 6.5 6.6

Výsledky srovnán´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Diskuze výsledk˚ u srovnán´ı . . . . . . . . . . . . . . . . . . . . . . . 44

7 Hybridn´ı model ve 2D 7.1 Spojitá ˇcást . . . . . ˇ asticová ˇcást . . . . 7.2 C´ 7.3 Sráˇzkové procesy . . 7.4 Zdroj ˇcástic . . . . . 7.5 Výsledky modelu . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

48 48 50 51 52 57

8 Hybridn´ı model ve 3D a jeho aplikace 61 8.1 Sonda koneˇcných rozmˇer˚ u . . . . . . . . . . . . . . . . . . . . . . . 61 8.2 Substrát s nerovným povrchem . . . . . . . . . . . . . . . . . . . . 68 9 Z´ avˇ er

78

Literatura

80

Pˇ r´ılohy

83

N´ azev pr´ ace: Studium interakce plazma-pevná látka pomoc´ı hybridn´ıho modelován´ı Autor: Vojtˇech Hrubý Katedra: Katedra fyziky povrch˚ u a plazmatu Vedouc´ı diplomov´ e pr´ ace: Prof. RNDr. Rudolf Hrach, DrSc. e-mail vedouc´ıho: [email protected] Abstrakt: V této práci je studována interakce plazmatu s vnoˇreným substrátem metodami poˇc´ıtaˇcové fyziky. Prezentovány jsou r˚ uzné zjednoduˇsené jednodimenzionáln´ı hybridn´ı modely válcové Langmuirovy sondy v plazmatu pozitivn´ıho sloupce doutnavého výboje v argonu. Efektivita a pˇresnost tˇechto model˚ u je porovnávána s ˇcásticovým a spojitým modelem. Dále je rozv´ıjen iterativn´ı hybridn´ı model zaloˇzený na kombinaci spojitého modelu ˇreˇseného metodou koneˇcných prvk˚ u a ˇcásticového modelu Monte Carlo. Modelem jsou ˇreˇseny u ´ lohy vyˇzaduj´ıc´ı plnˇe tˇr´ıdimenzionáln´ı pˇr´ıstup. Výsledky výpoˇctu jsou z´ıskány v pˇrijatelném ˇcase, zachována je do znaˇcné m´ıry také informace o jednotlivých ˇcástic´ıch. Na základˇe namˇeˇrených parametr˚ u plazmatu je modelována koneˇcná Langmuirova sonda v plazmatu s driftem. Dále je studována interakce plazmatu a kovového substrátu s nerovnostmi ve tvaru polokoule r˚ uzného polomˇeru. Na základˇe modelu skupiny polokoul´ı pravidelnˇe rozm´ıstˇených na substrátu jsou diskutovány zmˇeny elektrického pole v bl´ızkosti substrátu v závislosti na vzdálenosti nerovnost´ı. Kl´ıˇ cov´ a slova: Interakce plazma-pevná látka, hybridn´ı modelován´ı, n´ızkoteplotn´ı plazma Title: Study of plasma-solid interaction by means of hybrid modelling Author: Vojtˇech Hrubý Department: Department of Surface and Plasma Science Supervisor: Prof. RNDr. Rudolf Hrach, DrSc. Supervisor’s e-mail address: [email protected] Abstract: In this thesis the interaction of plasma with an immersed substrate is studied by means of the computational physics. Various simplified one-dimensional hybrid models of the cylindrical Langmuir probe in the plasma of the positive column of argon glow discharge are presented and their efficiency and accuracy are compared with particle and fluid models. An iterative hybrid model based on combination of FEM fluid model and Monte Carlo particle model is developed and used to solve problems demanding fully three-dimensional approach. The results are obtained in a reasonable time preserving a good deal of particle information. A finite cylindrical Langmuir probe in plasma with drift is modelled according to measured plasma properties. Interaction of plasma and a metal substrate with one semi-spherical protrusion of various radii is investigated. Groups of semi-spherical protrusions with different distances on the substrate are studied and the changes of distribution of electric field near the surface of substrate are discussed. Keywords: Plasma-solid interaction, hybrid modelling, low-temperature plasma

Kapitola 1 ´ Uvod Pojem plazma se v nedávné dobˇe dostal do povˇedom´ı ˇsiroké veˇrejnosti d´ıky velkoploˇsným plazmovým obrazovkám, pˇrestoˇze se s plazmatem setkáváme jiˇz mnoho let napˇr´ıklad v obyˇcejných záˇrivkách. Ve fyzice byl tento pojem zaveden v roce 1928 Irvingem Langmuirem [1], který tak nazval oblast ionizovaného plynu, ve které je vyrovnaný náboj iont˚ u a elektron˚ u, takˇze výsledný prostorový náboj je velmi malý. Fyzika plazmatu proˇsla do dneˇsn´ı doby prudkým vývojem. Pomohla nám lépe porozumˇet pro ˇclovˇeka tak vzdáleným jev˚ um jako polárn´ı záˇre nebo sluneˇcn´ı v´ıtr a pˇrinesla nám ˇradu nových technologi´ı, jejichˇz produkty dennˇe pouˇz´ıváme. Studium plazmatu nen´ı lehký u ´ kol. Plazma nám dokonce dovede aktivnˇe bránit ve zkoumán´ı jeho parametr˚ u. Stejné principy lze na druhou stranu vyuˇz´ıt v plazmatických technologi´ıch. To plat´ı také o interakci plazmatu s pevnou látkou. Od raných dob studia plazmatu se pouˇz´ıvá sondová diagnostika plazmatu. Pˇrestoˇze sonda plazma naruˇsuje, na základˇe teoretických, experimentáln´ıch a modelových pˇredpoklad˚ u lze do jisté m´ıry urˇcit parametry plazmatu nenaruˇseného. Z´ıskané poznatky ovˇsem nacház´ı své vyuˇzit´ı napˇr´ıklad v plazmatickém opracováván´ı povrch˚ u pevných látek. Procesy v plazmatu a pˇri interakci plazmatu s pevnou látkou jsou velmi sloˇzité, a proto je ˇrada u ´ loh analyticky neˇreˇsitelná. Poˇc´ıtaˇcové modelován´ı nám umoˇzn ˇ uje ˇreˇsit mnohem ˇsirˇs´ı okruh u ´ loh, který se s rostouc´ım výkonem poˇc´ıtaˇc˚ u neustále zvˇetˇsuje. Pˇresto mnohé u ´ lohy nelze zat´ım ani s pouˇzit´ım poˇc´ıtaˇc˚ u vyˇreˇsit snadno, a proto hledáme sloˇzitˇejˇs´ı, ale efektivnˇejˇs´ı metody. Hybridn´ı modelován´ı nám pˇrináˇs´ı ˇradu takových metod. Nˇekteré z nich budou v této práci prezentovány a porovnány. Dále ukáˇzeme ˇreˇsen´ı konkrétn´ıch u ´ loh vybranou metodou v plnˇe trojrozmˇerné geometrii.

6

Kapitola 2 Z´ akladn´ı poznatky fyziky plazmatu 2.1

Definice plazmatu

V knize [2] je plazma definováno následuj´ıc´ım zp˚ usobem: Plazma je kvazineutráln´ı plyn nabitých a neutr´ aln´ıch ˇc´ astic, který vykazuje kolektivn´ı chován´ı. Pojmy kvazineutráln´ı a kolektivn´ı chov´ an´ı je potˇreba rozebrat hloubˇeji. V neutráln´ım plynu na sebe ˇcástice p˚ usob´ı pouze pˇri sráˇzkách, protoˇze gravitaˇcn´ı s´ıla je zanedbatelná a ˇzádná elektrická s´ıla mezi nimi nep˚ usob´ı. Sráˇzky tak plnˇe rozhoduj´ı o pohybu ˇcástic. V plazmatu se vyskytuj´ı i nabité ˇcástice, které mohou vytváˇret prostorové náboje, které vedou ke vzniku elektrického pole. Coulombova s´ıla je dalekodosahová a náboje v plazmatu tak mohou na sebe p˚ usobit i na velké vzdálenosti. Pohyb ˇcástic v plazmatu tedy nezávis´ı jen na lokáln´ıch podm´ınkách, ale také na stavu plazmatu ve vˇetˇs´ıch vzdálenostech. Této vlastnosti ˇr´ıkáme kolektivn´ı chován´ı. D´ıky existenci pohyblivých nabitých ˇcástic plazma dokáˇze odst´ınit vloˇzené elektrické potenciály. Toto odst´ınˇen´ı vˇsak nen´ı dokonalé d´ıky tepelným pohyb˚ um nabitých ˇcástic. V okol´ı kaˇzdého elektrického potenciálu se vytvoˇr´ı st´ın´ıc´ı vrstva, jej´ıˇz ˇ ast nabitých ˇcástic totiˇz efektivn´ı polomˇer je dán tepelnou energi´ı nabitých ˇcástic. C´ m˚ uˇze z potenciálové jámy uniknout právˇe d´ıky tepelné energii, a proto malé potenciály mohou pronikat hloubˇeji do plazmatu a zp˚ usobovat elektrická pole i mimo st´ın´ıc´ı vrstvu. Na základˇe podobných u ´ vah lze odvodit takzvanou Debyeovu délku λD =

s

7

ε0 kT ∗ , ne2

(2.1)

´ Í POZNATKY FYZIKY PLAZMATU KAPITOLA 2. ZAKLADN

8

kde ε0 je permitivita vakua, k je Boltzmannova konstanta, n je koncentrace elektron˚ u v nenaruˇseném plazmatu, e je elementárn´ı náboj a 1 1 1 = + , ∗ T T+ T−

(2.2)

kde T+ teplota kladnˇe nabitých ˇcástic a T− je teplota zápornˇe nabitých ˇcástic. Tento vztah plat´ı v pˇr´ıpadˇe, ˇze kladnˇe i zápornˇe nabité ˇcástice maj´ı srovnatelné hmotnosti. V pˇr´ıpadˇe elektropozitivn´ıho plazmatu, které je tvoˇreno elektrony, kladnými ionty a neutrály, je st´ınˇen´ı vykonáváno pˇredevˇs´ım elektrony. Jejich teplota je pak pro Debyeovu délku rozhoduj´ıc´ı a T ∗ je pˇribliˇznˇe rovno Te . Debyeova délka je m´ırou st´ın´ıc´ı vzdálenosti. V této práci je ˇrada výpoˇct˚ u provedena pro pozitivn´ı sloupec doutnavého výboje v argonu pˇri tlaku 133 Pa s parametry pˇribliˇznˇe n = 1,0 · 1015 m−3 a Te = 23 200 K. Tomu odpov´ıdá Debyeova délka λD = 0,33 mm. Systém lze povaˇzovat za kvazineutr´ aln´ı, pokud jeho rozmˇer L výraznˇe pˇrevyˇsuje Debyeovu délku λD . To znamená, ˇze zavedeme-li do systému vnˇejˇs´ı potenciál nebo objev´ı-li se v systému lokáln´ı prostorový náboj, pˇreváˇzná ˇcást systému z˚ ustane d´ıky st´ınˇen´ı bez velkých elektrických pol´ı. Souˇcasnˇe mus´ı být koncentrace nabitých ˇcástic dostateˇcná, aby ve st´ın´ıc´ı vrstvˇe bylo velké mnoˇzstv´ı ˇcástic. Kdyby bylo ve vrstvˇe jen nˇekolik nabitých ˇcástic, Debyeovo st´ınˇen´ı by ztrácelo smysl. Na plazma klademe jeˇstˇe jednu podm´ınku: Frekvence plazmových oscilac´ı ωp =

s

e2 ne , ε0 me

(2.3)

kde e je elementárn´ı náboj, ne je koncentrace elektron˚ u a me je hmotnost elektronu, mus´ı být vˇetˇs´ı neˇz sráˇzková frekvence elektron˚ u s neutrály. V opaˇcném pˇr´ıpadˇe je dominantn´ı kinetika neutráln´ıho plynu a elektromagnetická interakce je potlaˇcená.

2.2

Klasifikace plazmatu

Plazma m˚ uˇze m´ıt mnoho r˚ uzných podob, které sice splˇ nuj´ı kritéria plazmatu, ale pˇresto se liˇs´ı natolik, ˇze vyˇzaduj´ı odliˇsné zp˚ usoby studia. Základn´ı dˇelen´ı plazmatu pˇrevezmeme z publikac´ı [3] a [2]. Podle m´ıry ionizace dˇel´ıme plazma na slabˇe ionizované a silnˇe ionizované . Ve slabˇe ionizovaném plazmatu je koncentrace nabitých ˇcástic o nˇekolik ˇrád˚ u menˇs´ı neˇz koncentrace neutráln´ıch ˇcástic. Nabité ˇcástice se proto sráˇz´ı pˇredevˇs´ım s neutráln´ımi ˇcásticemi. V silnˇe ionizovaném plazmatu je naopak koncentrace neutráln´ıch ˇcástic velmi malá, aˇz nulová. V tomto pˇr´ıpadˇe jsou dominantn´ı interakce mezi nabitými ˇcásticemi. Toto dˇelen´ı má velký fyzikáln´ı význam: Sráˇzky nabitých


9

ˇcástic s neutrály ve slabˇe ionizovaném plazmatu jsou dané krátkodosahovými silami, zat´ımco nabité ˇcástice v silnˇe ionizovaném plazmatu na sebe p˚ usob´ı Coulombovou silou, která je dalekodosahová. D´ıky tomu se vlastnosti slabˇe a silnˇe ionizovaného plazmatu výraznˇe liˇs´ı. Dˇelen´ı podle teploty je zaloˇzené sp´ıˇse na konvenci. Plazma povaˇzujeme za vysokoteplotn´ı, pokud je silnˇe ionizované ve smyslu Sahovy rovnice [2] ni T 3/2 ∼ 2,4 · 1021 exp(−Ui /kB T ), nn ni

(2.4)

která udává pomˇer hustot ionizovaných ni a neutráln´ıch nn atom˚ u plynu pˇri teplotˇe T , Ui je ionizaˇcn´ı energie plynu a kB je Boltzmannova konstanta v jednotkách SI. Sahova rovnice plat´ı pro plyn v tepelné rovnováze, kdy k ionizaci docház´ı sráˇzkami neutráln´ıch ˇcástic. Zvyˇsujeme-li teplotu p˚ uvodnˇe neutráln´ıho plynu, stupeˇ n ionizace z˚ ustává n´ızký, dokud se Ui nestane jen nevelkým násobkem kB T . Potom pomˇer ni /nn zaˇcne prudce nar˚ ustat, protoˇze atomy maj´ı dostateˇcnou tepelnou energii k odtrˇzen´ı elektronu pˇri vzájemné sráˇzce. Dalˇs´ım zvyˇsován´ım teploty dosáhneme aˇz plné ionizace. S vysokoteplotn´ı plazmatem se setkáváme v kosmu, v laboratorn´ıch podm´ınkách pak napˇr´ıklad pˇri termonukleárn´ı f´ uzi. N´ızkoteplotn´ı plazma tvoˇr´ı napˇr´ıklad pozitivn´ı sloupec doutnavého výboje. Podle teploty lze dˇelit plazma jeˇstˇe odliˇsným zp˚ usobem — na izotermické a neizotermické . Plazma totiˇz m˚ uˇze m´ıt ,,nˇekolik teplot najednou“. Kaˇzdá sloˇzka plazmatu (elektrony, r˚ uzné druhy iont˚ u a neutrály) m˚ uˇze m´ıt vlastn´ı energetické rozdˇelen´ı, kterému odpov´ıdá urˇcitá teplota nebo stˇredn´ı energie v pˇr´ıpadˇe, ˇze rozdˇelen´ı nen´ı maxwellovské. Pˇr´ıˇcinou je pˇr´ıliˇs malá frekvence sráˇzek mezi ˇcásticemi r˚ uzného druhu, která nestaˇc´ı k dosaˇzen´ı celkové tepelné rovnováhy. Teploty r˚ uzných sloˇzek se mohou liˇsit i o nˇekolik ˇrád˚ u. S t´ımto jevem se ˇcasto setkáváme v n´ızkoteplotn´ım slabˇe ionizovaném plazmatu. Vysokoteplotn´ı plazma je naopak izotermické. V této práci se budeme zabývat pˇredevˇs´ım n´ızkoteplotn´ım ˇcásteˇcnˇe ionizovaným neizotermickým plazmatem.

2.3

Kinetika plazmatu

Kinetická teorie plazmatu nám pˇrináˇs´ı d˚ uleˇzitý a uˇziteˇcný matematický aparát pro popis jev˚ u v plazmatu. Kl´ıˇcovým pojmem je rozdˇelovac´ı funkce fi (~r, ~v, t), která udává poˇcet ˇcástic i-tého druhu v objemovém elementu konfiguraˇcn´ıho prostoru dΩ = dxdydz a v objemovém elementu rychlostn´ıho prostoru dC = dvx dvy dvz v ˇcase t. Poloha element˚ u je dána vektory ~r a ~v . Koncentrace ˇcástic i-tého druhu


10

v bodˇe ~r je integrálem rozdˇelovac´ı funkce pˇres celý rychlostn´ı prostor C ni (~r, t) =

Z

C

fi (~r, ~v, t) dC.

(2.5)

Pomoc´ı rozdˇelovac´ı funkce definujeme stˇredn´ı hodnotu veliˇcin, které charakterizij´ı vlastnosti ˇcástic plazmatu. Stˇredn´ı hodnota veliˇciny ϕ(~r, ~v, t) je definována 1 hϕi = n

Z

C

ϕf dC.

(2.6)

Velmi d˚ uleˇzitou veliˇcinou pro naˇse dalˇs´ı výpoˇcty je rychlost ~v . Pomoc´ı stˇredn´ı hodnoty lze rychlost rozdˇelit na dvˇe sloˇzky ~v = ~c + ~u,

(2.7)

kde ~c je chaotická rychlost a ~u je driftov´ a rychlost. Pro tyto rychlosti plat´ı vztahy h~c i = 0,

a ~u = h~v i.

(2.8)

S chaotickou a driftovou rychlost´ı u ´ zce souvis´ı stˇredn´ı hodnota kinetické energie ˇcástic, pro kterou plat´ı hEi = =

D

1 mv 2 2

1 mu2 2

E

+

= 12 m(hu2 i + hc2 i + h~u · ~c i) = 12 m(u2 + hc2 i)

(2.9)

1 mhc2 i. 2

Stˇredn´ı hodnota energie je tedy souˇctem energie driftového a chaotického pohybu. Jiˇz nˇekolikrát jsme pouˇzili pojem teplota, ale nedefinovali jsme jej, protoˇze u ´ zce souvis´ı s rozdˇelovac´ı funkc´ı a se stˇredn´ı hodnotou energie. Je-li plyn v tepelné rovnováze, nejpravdˇepodobnˇejˇs´ı rozdˇelen´ı rychlost´ı ˇcástic je Maxwellovo rozdˇelen´ı. V jednorozmˇerném pˇr´ıpadˇe, napˇr´ıklad pro komponentu vx , je podle [2] pˇr´ısluˇsná rozdˇelovac´ı funkce dána vztahem f (vx ) = n

s

m 1 exp − mvx2 /kB T , 2πkB T 2

(2.10)

kde m je hmotnost ˇcástice, kB je Boltzmannova konstanta a T je teplota. Vypoˇc´ıtáme-li pro tento pˇr´ıpad stˇredn´ı hodnotu kinetické energie ˇcástice, dostaneme 1 hEi = kB T. 2

(2.11)

Pro tˇr´ırozmˇerný pohyb s izotropn´ım rozdˇelen´ım je Maxwellovo rozdˇelen´ı dáno vztahem 3/2 m 1 2 2 f (v) = n 4π v exp − mv /kB T , (2.12) 2πkB T 2


11

který lze odvodit ze vztahu (2.10) integrac´ı pˇres vˇsechny smˇery v tˇr´ırozmˇerném prostoru. Pro tento pˇr´ıpad vypoˇc´ıtáme stˇredn´ı hodnotu kinetické energie ˇcástice 3 hEi = kB T. 2

(2.13)

Obecnˇe plat´ı, ˇze na kaˇzdý stupeˇ n volnosti pˇripadá energie 21 kB T . Ve fyzice plazmatu je obvyklé udávat teplotu ˇcástic v jednotkách energie. S ohledem na závislost na dimenzi se uvád´ı energie odpov´ıdaj´ıc´ı právˇe kB T , napˇr´ıklad pro kB T = 1 eV je T = 11 600 K.

2.4

Boltzmannova rovnice

Rozdˇelovac´ı funkce f (~r, ~v, t) kaˇzdé komponenty plazmatu mus´ı splˇ novat Boltzmannovu rovnici F~ δf ∂f + ~v · ∇f + · ∇c f = , (2.14) ∂t m δt kde ∇c je gradient v rychlostn´ım prostoru, m je hmotnost ˇcástice, F~ je makroskopická s´ıla p˚ usob´ıc´ı na ˇcástice komponenty a ˇclen na pravé stranˇe zahrnuje sráˇzky s ostatn´ımi ˇcásticemi, tedy pˇredevˇs´ım mikroskopické interakce. Pokud je sráˇzkový ˇclen zanedbatelný a s´ıla F~ je výhradnˇe elektromagnetická, rovnice (2.14) se nazývá Vlasovova. ˇ sen´ı Boltzmannovy rovnice je velmi sloˇzité. Velmi v´ Reˇ yznamné jsou takzvané momenty Boltzmannovy rovnice, které z´ıskáme integrac´ı Boltzmannovy rovnice pˇres rychlostn´ı prostor. Tyto momenty jsou pouze d˚ usledky Boltzmannovy rovnice. Nultý moment z´ıskáme integrac´ı samotné Boltzmannovy rovnice pˇres rychlostn´ı prostor. Výpoˇcet je proveden napˇr´ıklad v [2]. Pˇredpokládáme-li pouze elektromagnetickou makroskopickou s´ılu, dostaneme rovnici kontinuity ∂n + ∇ · (n~u) = 0. ∂t

(2.15)

Docház´ı-li v plazmatu ke sráˇzkám, pˇri kterých se mˇen´ı poˇcet ˇcástic urˇcitého druhu, napˇr´ıklad ionizace a rekombinace, na pravé stranˇe rovnice kontinuity se objev´ı takzvané ˇcleny vzniku a zániku ˇc´ astic. Tyto ˇcleny maj´ı obvykle tvar souˇcinu ±k

Y

nj ,

(2.16)

j

kde k je rychlostn´ı koeficient interakce a souˇcin prob´ıhá pˇres vˇsechny druhy ˇcástic, které se sráˇzky u ´ˇcastn´ı.


12

Prvn´ı moment Boltzmannovy rovnice z´ıskáme tak, ˇze ji vynásob´ıme vektorem rychlosti a opˇet vyintegrujeme pˇres celý rychlostn´ı prostor. Výsledek nakonec vynásob´ıme m a dostaneme podle [2] pohybovou rovnici pro tekutinu !

∂~u ~ + ~u × B) ~ − ∇ · P + P~ij , mn + (~u · ∇)~u = qn(E ∂t

(2.17)

kde P je tenzor napˇet´ı a P~ij vyjadˇruje zmˇenu hybnosti vlivem sráˇzek. Vynásoben´ım Boltzmannovy rovnice ˇclenem 12 mv 2 a integrac´ı bychom analogicky dostali rovnici pro tok tepla. Prvn´ı dva aˇz tˇri momenty tvoˇr´ı základn´ı rovnice spojitých model˚ u plazmatu. Jejich ˇreˇsen´ı je ménˇe nároˇcné neˇz ˇreˇsen´ı Boltzmannovy kinetické rovnice. Na druhou stranu tyto výsledky nedávaj´ı u ´ plnou informaci o rozdˇelovac´ı funkci.

2.5

Diagnostika n´ızkoteplotn´ıho plazmatu

Pro mˇeˇren´ı parametr˚ u plazmatu byla vyvinuta ˇrada metod. Základn´ı metody lze hrubˇe rozdˇelit na ˇctyˇri skupiny: • sondové, • vysokofrekvenˇcn´ı, • optické, • korpuskulárn´ı. Sondová diagnostika je velmi bl´ızká tématu této práce, protoˇze je zaloˇzena na interakci plazmatu s kovovou sondou. Metody sondové diagnostiky, které navrhl Irving Langmuir (viz napˇr´ıklad [4]), patˇr´ı mezi nejstarˇs´ı metody diagnostiky plazmatu. Sondy mohou m´ıt r˚ uzné geometrie. Nejˇcastˇeji se setkáváme s rovinnou, válcovou a kulovou sondou, které jsou schematicky znázornˇeny na obrázku 2.1. Struˇcnˇe pop´ıˇseme základy sondové diagnostiky podle [5] a [3]. Vloˇz´ıme-li do plazmatu vodiˇc na daném potenciálu, separac´ı náboje se vytvoˇr´ı v jeho okol´ı st´ın´ıc´ı vrstva (sheath). Zmˇeˇr´ıme-li voltampérovou charakteristiku jedné sondy podle obrázku 2.2 nebo dvou sond podle 2.3, m˚ uˇzeme pˇribliˇznˇe urˇcit nˇekteré parametry sheathu a nenaruˇseného plazmatu. Nejjednoduˇsˇs´ı analýza voltampérové charakteristiky zaloˇzená na takzvané bezesráˇzkové teorii vyˇzaduje splnˇen´ı následuj´ıc´ıch poˇzadavk˚ u: 1. V sheathu se nabité ˇcástice pohybuj´ı bez sráˇzek. 2. Sonda neemituje.


rovinn´ a

v´ alcov´ a

13

kulov´ a

Obrázek 2.1: Nejˇcastˇejˇs´ı sondy pro diagnostiku n´ızkoteplotn´ıho plazmatu. Modrá barva znázorˇ nuje izolaci, napˇr´ıklad sklo aparatury. 3. Za vrstvou je nenaruˇsené plazma. 4. Proudy na sondu jsou nezávislé. 5. Elektrony a ionty maj´ı v plazmatu Maxwellovo rozdˇelen´ı. Prvn´ı poˇzadavek pˇredstavuje zásadn´ı omezen´ı této metody jen na n´ızké tlaky. Pˇr´ıklad namˇeˇrené voltampérové charakteristiky jedné sondy je uveden na obrázku 2.4. Mˇeˇren´ı bylo provedeno se sondou o polomˇeru 22,5 µm a délce 2,5 mm pˇri tlaku 10 Pa v pulzn´ım magnetronu. Z charakteristiky lze podle bezesráˇzkové teorie urˇcit v prvé ˇradˇe teplotu elektron˚ u a koncentraci nabitých ˇcástic. Pˇredpokládáme plazma tvoˇrené elektrony, jednonásobnými kladnými ionty a neutráln´ımi ˇcásticemi. Teplotu elektron˚ u z´ıskáme ze smˇernice závislosti logaritmu elektronového proudu na napˇet´ı ln |I− | = ln |I−p | +

e(U − Up ) , kB Te

(2.18)

kde I−p je nasycený elektronový proud, Up je potenciál plazmatu a Te je teplota elektron˚ u. Lineárn´ı závislost odpov´ıdá maxwellovskému rozdˇelen´ı energie elektron˚ u. Nen´ı-li závislost lineárn´ı, lze skuteˇcnou rozdˇelovac´ı funkci pˇribliˇznˇe urˇcit analýzou této závislosti. Koncentraci iont˚ u urˇc´ıme z iontového proudu I+f , který je extrapolac´ı iontového proudu do bodu U = Uf . Koncentrace iont˚ u je dána vztahem |I+f | ni = eS

s

mi , kB Te

kde e je elementárn´ı náboj, S je plocha sondy a mi je hmotnost iontu.

(2.19)

´ Í POZNATKY FYZIKY PLAZMATU KAPITOLA 2. ZAKLADN Vd

14

R

A

K

A V

Obrázek 2.2: Schéma jednosondové diagnostiky v pozitivn´ım sloupci doutnavého výboje Koncentraci elektron˚ u z´ıskáme na základˇe nasyceného elektronového proudu |I−p | ne = 4 eS

s

πme , 8kB Te

(2.20)

kde me je hmotnost elektronu. Koncentrace elektron˚ u a iont˚ u by si mˇely být rovny z d˚ uvodu kvazineutrality plazmatu. Ve skuteˇcnosti se výsledky liˇs´ı kv˚ uli nepˇresnostem mˇeˇren´ı a nedostateˇcnému splnˇen´ı podm´ınek platnosti teorie. Ve voltampérové charakteristice dvou sond je dominantn´ı vliv iont˚ u. Jsou-li sondy shodné, charakteristika je symetrická, tvoˇrená iontovými vˇetvemi a pˇrechodovou oblast´ı v okol´ı nulového proudu. V tomto pˇr´ıpadˇe m˚ uˇzeme obdobnˇe urˇcit koncentraci iont˚ u a teplotu elektron˚ u.

2.6

Interakce plazmatu s povrchy pevn´ ych l´ atek

Nabité ˇcástice z plazmatu, které projdou st´ın´ıc´ı vrstvou, dopadaj´ı na povrch vnoˇrené pevné látky – sondy, substrátu v plazmochemických technologi´ıch, na stˇenu výbojky, apod. Co se s tˇemi ˇcásticemi dále dˇeje, závis´ı na vlastnostech vnoˇrené látky. Pokud se jedná o kovovou elektrodu zapojenou do vnˇejˇs´ıho obvodu, ve zjednoduˇsené teorii sondové diagnostiky se pˇredpokládá, ˇze nabitá ˇcástice bude z plazmatu odstranˇena a jej´ı náboj pˇrispˇeje k sondovému proudu. Elektron z plazmatu

´ Í POZNATKY FYZIKY PLAZMATU KAPITOLA 2. ZAKLADN Vd

15

R

A

K

A V

Obrázek 2.3: Schéma dvousondové diagnostiky v pozitivn´ım sloupci doutnavého výboje vstoup´ı do elektrody a dále bude procházet vnˇejˇs´ım obvodem. Iont, kladný nebo záporný, bude na povrchu sondy neutralizován, jeho náboj vstoup´ı do vnˇejˇs´ıho obvodu a v plazmatu z˚ ustane z iontu neutráln´ı ˇcástice. Pˇri pˇresnˇejˇs´ı analýze je vˇsak nutné uvaˇzovat i dalˇs´ı fyzikáln´ı procesy, zejména pro elektrony, z nichˇz nejd˚ uleˇzitˇejˇs´ı je elektronová sekundárn´ı emise ze sondy pod dopadem rychlejˇs´ıch elektron˚ u z plazmatu, která m˚ uˇze jak zeslabovat, tak i zesilovat sondový proud. Pokud se jedná o substrát dielektrický nebo i o kovový, jenˇz ale nen´ı pˇripojen na vnˇejˇs´ı obvod, dojde k nab´ıjen´ı povrchu substrátu dopadem nabitých ˇcástic z plazmatu. Pˇritom se projev´ı rozd´ılná pohyblivost jednotlivých komponent plazmatu, takˇze napˇr. v elektropozitivn´ım plazmatu se substrát nab´ıj´ı zápornˇe. D˚ usledkem bude lokáln´ı elektrické pole, jeˇz ovlivn´ı plazma v okol´ı substrátu. Tento efekt je vˇzdy pˇr´ıtomný ve sklenˇených výbojkách, kde vede k vytvoˇren´ı radiáln´ıho elektrického pole a k ambipolárn´ı difuzi nabitých ˇcástic plazmatu. Speciáln´ım pˇr´ıpadem jsou r˚ uzné plazmochemické technologie nebo sondová diagnostika v chemicky aktivn´ım plazmatu. Zde je substrátem p˚ uvodnˇe kovová elektroda nebo sonda zapojená do vnˇejˇs´ıho obvodu, na n´ıˇz je vloˇzeno urˇcité pˇredpˇet´ı maj´ıc´ı za u ´ kol vytáhnout z plazmatu pˇr´ısluˇsnˇe nabité ionty. Pokud ale tyto ˇcástice zaˇcnou na povrchu substrátu chemicky reagovat a vytváˇret dielektrickou vrstvu, povrch pˇrestane odvádˇet náboj do vnˇejˇs´ıho obvodu a zaˇcne se nab´ıjet. Proti ˇcistˇe dielektrickému substrátu nebo stˇenˇe výbojky je ale situace odliˇsná v tom, ˇze vytváˇrená vrstva je obvykle tenká a má nenulovou elektrickou vodivost. Výsledný


16

0,002 Namˇeˇren´ a data

0,001 I [A]

Ip

0

−0,001 −20

Up

Uf

−10

0

10

U [V]

Obrázek 2.4: Pˇr´ıklad namˇeˇrené voltampérové charakteristiky jedné sondy. Data poskytl Mgr. Jan Klusoˇ n. Uf je plovouc´ı potenciál a Up je potenciál plazmatu s odpov´ıdaj´ıc´ım proudem Ip . náboj na povrchu rostouc´ı vrstvy je proto dán kombinac´ı náboje pˇrivádˇeného z plazmatu a náboje odvádˇeného dielektrikem do vnˇejˇs´ıho obvodu. Výsledkem bude napˇet’ový spád v dielektriku a pokles pˇredpˇet´ı substrátu. Tento proces se mus´ı brát v u ´ vahu pˇri plazmochemických technologi´ıch a naopak m˚ uˇze být vyuˇz´ıván pˇri sondové diagnostice chemicky aktivn´ıho plazmatu. Analýza sondového proudu v takovém pˇr´ıpadˇe je vˇsak velmi obt´ıˇzná a bez aplikace poˇc´ıtaˇcového modelován´ı tˇechto proces˚ u prakticky vylouˇcená. Vedle uvedených u ´ vah obecnˇe platných existuj´ı i speciáln´ı pˇr´ıpady, kdy procesy na elektrodách nabývaj´ı zvláˇstn´ı d˚ uleˇzitost. Pˇr´ıkladem m˚ uˇze být sondová diagnostika pomoc´ı ˇzhavé sondy, kdy rozhoduj´ıc´ım procesem je termoemisi elektron˚ u ze sondy.

Kapitola 3 Poˇ c´ıtaˇ cov´ a fyzika Ve dvacátém stolet´ı se k experimentáln´ı a teoretické fyzice pˇridal tˇret´ı výrazný smˇer — poˇc´ıtaˇcová fyzika. Pˇresná definice tohoto pojmu zˇrejmˇe jeˇstˇe nen´ı ustálená. Vˇetˇsina fyzik˚ u dnes ke své práci pouˇz´ıvá poˇc´ıtaˇc. Poˇc´ıtaˇcovou fyzikou vˇsak nazýváme pouze takovou vˇedeckou práci nebo metodiku vˇedecké práce, pˇri které je práce poˇc´ıtaˇce dominantn´ı metodou ˇreˇsen´ı problému. P˚ uvodn´ı vzájemnou spolupráci teoretické a experimentáln´ı fyziky poˇc´ıtaˇcová fyzika rozˇs´ıˇrila o mnoho nových moˇznost´ı: Pro teoretickou fyziku se stala v´ıtaným pomocn´ıkem tam, kde analytické výpoˇcty nejsou moˇzné, experimentáln´ı fyzice dala moˇznost dopˇredu simulovat nároˇcné experimenty, které se tˇeˇzko opakuj´ı. S rostouc´ım výkonem dostupné poˇc´ıtaˇcové techniky a s rozvojem nových metod se spolupráce dále prohlubuje. V krátké historii poˇc´ıtaˇcové fyziky vznikla ˇrada d´ılˇc´ıch smˇer˚ u: • Poˇc´ıtaˇcové modelován´ı V poˇc´ıtaˇci vytváˇr´ıme zjednoduˇsený model studovaného fyzikáln´ıho jevu. Model ˇreˇs´ıme ˇradou numerických metod. Obvykle se jedná o ˇreˇsen´ı obyˇcejných a parciáln´ıch diferenciáln´ıch rovnic, ˇreˇsen´ı soustav diferenciáln´ıch rovnic a statistické výpoˇcty. Na základˇe výsledk˚ u modelu m˚ uˇzeme pˇredpov´ıdat výsledky experiment˚ u, optimalizovat technologické procesy a podobnˇe. • Zpracován´ı obrazu Z obrazu, který je výsledkem pozorován´ı, z´ıskáváme poˇc´ıtaˇcovými metodami fyzikáln´ı, biologické a dalˇs´ı informace. Zpracován´ı obrazu se ˇsiroce uplatˇ nuje v medic´ınˇe, ˇr´ızen´ı experimentu, automatizaci výroby, atd. • Poˇc´ıtaˇcová grafika Výsledky mˇeˇren´ı nebo výpoˇctu obvykle poˇzadujeme znázornit v grafické formˇe. Za t´ımto u ´ˇcelem jiˇz vznikla ˇrada uˇziteˇcných program˚ u, které obvykle plnˇe postaˇcuj´ı naˇsim potˇrebám.

17

ˇ ÍTACOV ˇ ´ FYZIKA KAPITOLA 3. POC A

18

• Integráln´ı transformace ˇ Radu u ´ loh lze efektivnˇeji ˇreˇsit tak, ˇze je pˇrevedeme integráln´ı transformac´ı do obrazu, se kterým provedeme vhodnou operaci, a provedeme zpˇetnou integráln´ı transformaci. Popularita integráln´ıch transformac´ı prudce vzrostla objeven´ım rychlé Fourierovy transformace v roce 1965. V dneˇsn´ı dobˇe se stále v´ıce uplatˇ nuje tzv. waveletová transformace. • Symbolické manipulace Na rozd´ıl od numerických metod symbolickými manipulacemi z´ıskáme výsledek ve tvaru matematického vzorce. Vˇetˇsinu fyzikáln´ıch problém˚ u nicménˇe t´ımto zp˚ usobem ˇreˇsit nelze. • Modern´ı smˇery poˇc´ıtaˇcové fyziky Do této kategorie patˇr´ı napˇr´ıklad evoluˇcn´ı programován´ı a neuronové s´ıtˇe. Smˇery poˇc´ıtaˇcové fyziky se stále vyv´ıjej´ı a výˇse uvedený seznam proto nen´ı koneˇcný au ´ plný. Nejˇcastˇeji pouˇz´ıvaným smˇerem je poˇc´ıtaˇcové modelov´ an´ı, kterým se zabývá také tato práce. Pˇr´ıbuzným oborem je matematické modelov´ an´ı, jehoˇz c´ılem je matematická formulace fyzikáln´ıch jev˚ u a jejich matematická analýza — ˇreˇsitelnost, stabilita, konvergence a podobnˇe. C´ılem poˇc´ıtaˇcové modelován´ı, na rozd´ıl od matematického, je nalezen´ı konkrétn´ıho ˇreˇsen´ı problému pomoc´ı poˇc´ıtaˇce. ˇ s´ıme-li fyzikáln´ı problém poˇc´ıtaˇcovým modelován´ım, procház´ıme podle [6] Reˇ následuj´ıc´ımi kroky: 1. formulace problému, 2. vytvoˇren´ı modelu, 3. ˇreˇsen´ı modelu, 4. analýza výsledk˚ u. Poˇc´ıtaˇcové modelován´ı má ovˇsem také ˇradu u ´ skal´ı. Prakticky kaˇzdý model je zjednoduˇsený, protoˇze nev´ıme vˇse a moˇznosti poˇc´ıtaˇc˚ u jsou omezené. Pˇri zjednoduˇsován´ı nikdy nemáme jistotu, ˇze jsme nezanedbali podstatný rys problému. Výsledky modelu nevypov´ıdaj´ı o studovaném jevu pˇr´ımo, odpov´ıdaj´ı pouze naˇsemu modelu. Z tˇechto d˚ uvod˚ u je velmi d˚ uleˇzité srovnávat výsledky s experimentem a ˇ také sledovat, zda odpov´ıdaj´ı poznatk˚ um teoretické fyziky. Casto tato ,,zpˇetná vazba“ vede k upˇresnˇen´ı modelu, pˇr´ıpadnˇe k iteraˇcn´ımu postupu. Poˇc´ıtaˇcové modelován´ı zahrnuje ˇradu r˚ uzných technik. Nejˇcastˇejˇs´ı metody lze shrnout do následuj´ıc´ıch kategori´ı:

ˇ ÍTACOV ˇ ´ FYZIKA KAPITOLA 3. POC A

19

ˇ asticové metody • C´ Studovaný jev popisujeme na mikroskopické u ´ rovni. Podrobnˇe sledujeme chován´ı jednotlivých ˇcástic. Nejˇcastˇeji se setkáváme s deterministickou metodou molekulárn´ı dynamiky a stochastickou metodou Monte Carlo. Výsledky tˇechto metod jsou velmi detailn´ı, v ˇradˇe pˇr´ıpad˚ u jsou vˇsak takové výpoˇcty pˇr´ıliˇs nároˇcné na strojový ˇcas. • Spojité metody Jev pop´ıˇseme makroskopicky pomoc´ı veliˇcin jako hustota, teplota, rychlost proudˇen´ı a ˇreˇs´ıme parciáln´ı diferenciáln´ı rovnice, které vyjadˇruj´ı zákony zachován´ı — nejˇcastˇeji energie, hybnosti a rovnici kontinuity. Výsledky jsou podstatnˇe ménˇe detailn´ı neˇz výsledky ˇcásticových model˚ u. Rychlost výpoˇctu vˇsak bývá výraznˇe vyˇsˇs´ı. • Hybridn´ı metody Pˇredstavuj´ı kompromisn´ı ˇreˇsen´ı, kdy kombinujeme v´ yˇse uvedené metody modelován´ı k dosaˇzen´ı co nejpˇresnˇejˇs´ıch výsledk˚ u v krátkém ˇcase. V následuj´ıc´ı kapitole rozebereme detailnˇeji ˇcásticové, spojité a hybridn´ı modely pouˇz´ıvané ve fyzice plazmatu.

Kapitola 4 Hybridn´ı modelov´ an´ı Se slovem hybridn´ı se setkáváme v mnoha oblastech lidské ˇcinnosti. Zpravidla jej pouˇz´ıváme tam, kde chceme zd˚ uraznit, ˇze jsme vyuˇzili ty lepˇs´ı vlastnosti z kaˇzdé ˇcásti celku — pˇri kˇr´ıˇzen´ı v zahradnictv´ı, v automobilovém pr˚ umyslu a podobnˇe. Jako snadno pˇredstavitelný pˇr´ıklad nám poslouˇz´ı hybridn´ı trolejbus. Trolejbus i autobus vˇsichni známe a v´ıme, jaké jsou jejich výhody a nevýhody. Trolejbus je tichý a ˇsetrný k mˇestskému ovzduˇs´ı, ale je závislý na trolejovém veden´ı. Autobus je nezávislý, ale v dopravn´ı zácpˇe dokáˇze vydatnˇe pˇrispˇet ke smogu. Zkusme spojit výhody autobusu a trolejbusu. Od takového vozidla budeme ˇcekat alespoˇ n ˇ ˇcásteˇcnou nezávislost i ˇsetrnost k prostˇred´ı mˇesta. Reˇsen´ı m˚ uˇze m´ıt r˚ uzné podoby — napˇr´ıklad hybridn´ı trolejbus. Do trolejbusu zkrátka zabudujeme spalovac´ı motor. Ve mˇestˇe bude trolejbus napájený z trolejového veden´ı a na vesnici pojede na naftu. Zdánlivˇe je to jednoduché, ale k u ´ spˇeˇsnému v´ ysledku vede nelehká cesta. Mus´ıme sestrojit sloˇzitˇejˇs´ı pˇrevodovku, vyrovnat se s vyˇsˇs´ı celkovou hmotnost´ı a podobnˇe. S obdobnými problémy se setkáváme v modelován´ı ve fyzice. Nˇekteré metody jsou velmi pˇresné a ˇcasovˇe nároˇcné, jiné dávaj´ı výsledky nesrovnatelnˇe rychleji, ale s velkou chybou. Hybridn´ı modely by mˇely vést k co nejlepˇs´ım výsledk˚ um v pˇrimˇeˇreném ˇcase. Dosahuj´ı toho kombinac´ı r˚ uzných pˇr´ıstup˚ u k jednomu fyzikáln´ımu problému. Obvykle vˇsak nelze pˇrevz´ıt jen ty lepˇs´ı vlastnosti, a proto mus´ıme hledat kompromisn´ı ˇreˇsen´ı. V následuj´ıc´ım textu budou popsány nejprve základn´ı metody modelován´ı v n´ızkoteplotn´ım plazmatu. Plynule na nˇe naváˇze popis hybridn´ıch model˚ u.

4.1

Vztahy metod

Základn´ı srovnán´ı metod pouˇz´ıvaných v poˇc´ıtaˇcovém modelován´ı n´ızkoteplotn´ıho plazmatu lze naj´ıt napˇr´ıklad v [6]. N´ıˇze diskutované metody jsou struˇcnˇe zobrazeny 20

´ Í KAPITOLA 4. HYBRIDNÍ MODELOVAN

21

na diagramu na obrázku 4.1. Tato práce se zabývá pˇredevˇs´ım ,,velkým hybridn´ım modelem“, který je vymezen modrou barvou.

Molekulárn´ı dynamika ˇ asticov´ C´ y model Monte Carlo Plazma

ˇ sen´ı BKR Reˇ Spojit´ y model

ˇ Obrázek 4.1: Metody modelován´ı v plazmatu. Cervenou barvou je vyznaˇcen ,,malý hybridn´ı model“, modrou barvou ,,velký hybridn´ı model“.

4.2

ˇ asticov´ C´ e modelov´ an´ı

ˇ asticové modely poˇc´ıtaj´ı pohyb a interakce jednotlivých ˇcástic. Podrobnˇe se této C´ problematice vˇenuje napˇr´ıklad [7]. Ve fyzice plazmatu obvykle postaˇcuje klasická Newtonova pohybová rovnice d2~r F~ = m 2 , (4.1) dt kde F~ je s´ıla p˚ usob´ıc´ı na ˇcástici, m je hmotnost ˇcástice a ~r je polohový vektor ˇcástice. Známe-li s´ıly p˚ usob´ıc´ı na kaˇzdou ˇcástici a poˇcáteˇcn´ı polohy a rychlosti vˇsech ˇcástic, m˚ uˇzeme tuto obyˇcejnou diferenciáln´ı rovnici ˇreˇsit. Polohy a rychlosti ˇcástic postupnˇe poˇc´ıtáme v ˇcasových kroc´ıch ∆t. Pro výpoˇcet existuje ˇrada algoritm˚ u: 1. Euler˚ uv algoritmus: Patˇr´ı k metodám prvn´ıho ˇrádu pˇresnosti v ˇcase. Nové polohy, rychlosti a s´ıly ˇcástic z´ıskáme ze schématu ~rik+1 = ~rik + ~vik ∆t ~vik+1 = ~vik + m1i F~ik ∆t F~ik+1 = . . . Index k odpov´ıdá ˇcasu, i je index ˇcástice.

(4.2)


22

2. Verletovy algoritmy: Diferenˇcn´ı schéma v této práci pouˇz´ıvaného Verletova algoritmu druhého ˇrádu pˇresnosti v ˇcase je 1 ~k ~rik+1 = ~rik + ~vik ∆t + 2m Fi ∆t2 i F~ik+1 = . . . 1 ~ k + F~ik+1 ∆t ~vik+1 = ~vik + 2m F i i

(4.3)

Verletovy algoritmy pˇredstavuj´ı dobrý kompromis mezi pˇresnost´ı a nároˇcnost´ı výpoˇctu a jsou vhodné pro systémy mnoha ˇcástic. Nevýhodou je podm´ınka, ˇze s´ıla nesm´ı záviset na rychlosti ˇcástic. To je ovˇsem v této práci splnˇeno. 3. Algoritmy vyˇsˇs´ıch ˇrád˚ u: Nejsou v této práci pouˇz´ıvány, protoˇze jsou výpoˇcetnˇe nároˇcné a jejich pˇresnost by mimo jiné d´ıky sráˇzkovým proces˚ um v plazmatu z˚ ustala nevyuˇzita. V celé práci je pouˇz´ıván výˇse popsaný Verlet˚ uv algoritmus. Pro ˇreˇsen´ı pohybových rovnic dále potˇrebujeme znát s´ıly, které p˚ usob´ı na ˇcástice. V n´ızkoteplotn´ım plazmatu jsou významné dva druhy interakc´ı ˇcástic: Coulombické silové p˚ usoben´ı nabitých ˇcástic a jejich sráˇzky s neutráln´ımi ˇcásticemi. Z pohledu poˇc´ıtaˇcové fyziky tyto interakce ˇreˇs´ıme po ˇradˇe metodou molekulárn´ı dynamiky a metodou Monte Carlo. Spojen´ı tˇechto metod nˇekdy bývá nazýváno ,,malým hybridn´ım modelem“. Výpoˇcet elektrických sil p˚ usob´ıc´ıch na ˇcástice lze ˇreˇsit ˇradou algoritm˚ u. V této práci pouˇz´ıváme standardn´ı algoritmus PIC (Particle-in-Cell ) ve variantˇe CIC (Cloud-in-Cell ), viz napˇr. [7]. V pracovn´ı oblasti sestroj´ıme s´ıt’ vhodné geometrie a v uzlových bodech vypoˇc´ıtáme potenciál elektrického pole ϕ ˇreˇsen´ım Poissonovy rovnice na základˇe hustoty náboje ˇcástic v s´ıti ̺ ̺ (4.4) ∆ϕ = − . ε0 Intenzitu elektrického pole pak vypoˇc´ıtáme vztahem ~ = −∇ϕ. E

(4.5)

Existuje ˇrada moˇznost´ı, jak ˇreˇsit Poissonovu rovnici. Ve vyˇsˇs´ıch dimenz´ıch pˇredstavuje ˇreˇsen´ı Poissonovy rovnice ˇcasovˇe velmi nároˇcný u ´ kol. Podle práce [8] v 1D a 2D tvoˇr´ı výpoˇcet sil, resp. ˇreˇsen´ı Poissonovy rovnice ménˇe neˇz 50% celkové doby výpoˇctu. V 3D výpoˇctu je to pˇribliˇznˇe 99%. V 1D pˇr´ıpadˇe se obvykle pouˇz´ıvá pˇr´ımé ˇreˇsen´ı Thomasovým algoritmem [9]. V 2D zvláˇstˇe dˇr´ıve byla ˇcasto pouˇz´ıvána superrelaxaˇcn´ı metoda (Succesive Overrelaxation, SOR, viz [10]), která vyniká snadnou implementac´ı, ale postupnˇe byla nahrazena efektivnˇejˇs´ımi metodami. Pˇr´ıkladem je metoda LU dekompozice ˇr´ıdké matice s knihovnou UMFPACK [11].


23

V plnˇe trojrozmˇerném modelu se setkává v´ıce nesnáz´ı. Optimalizaci ˇreˇsen´ı Poissonovy rovnice je potˇreba vˇenovat velkou pozornost. V´ ychodiskem mohou být silné metody jako multigridy, konjugované gradienty a metoda koneˇcných prvk˚ u. Souˇcasnˇe má vˇsak 3D model velmi vysoké nároky na poˇcet ˇcástic v modelu. Mámeli krychlovou s´ıt’ o rozmˇeru 200 × 200 × 200 uzl˚ u a poˇzadujeme pr˚ umˇernˇe deset 7 ˇcástic v kaˇzdé buˇ nce, poˇcet ˇcástic v modelu bude 8 · 10 , coˇz je na hranici kapacity pamˇeti bˇeˇzného osobn´ıho poˇc´ıtaˇce. Bohuˇzel tyto poˇzadavky jsou minimáln´ı, má-li ´ ernˇe tˇemto poˇzadavk˚ m´ıt model v˚ ubec smysl. Umˇ um také roste doba výpoˇctu a výpoˇcet i relativnˇe malého modelu m˚ uˇze trvat ˇrádovˇe mˇes´ıce. Odtud mimo jiné plyne snaha vytváˇret hybridn´ı modely, které by mˇely v´ yraznˇe niˇzˇs´ı hardwarové a ˇcasové nároky. Mezi slabiny metody molekulárn´ı dynamiky patˇr´ı takzvaný nefyzikáln´ı ohˇrev. Pˇri numerickém ˇreˇsen´ı pohybových rovnic se d´ıky koneˇcnému ˇcasovému kroku m´ırnˇe zvyˇsuje kinetická energie ˇcástic. Vol´ıme proto dostateˇcnˇe malý ˇcasový krok, aby nefyzikáln´ı ohˇrev bˇehem poˇc´ıtaného ˇcasového intervalu byl zanedbatelný. Jiˇz bylo uvedeno, ˇze sráˇzky nabitých ˇcástic s neutráln´ımi ˇcásticemi modelujeme metodou Monte Carlo. Kaˇzdý druh sráˇzky je charakterizován pˇredevˇs´ım sráˇzkovým pr˚ uˇrezem, který je obvykle funkc´ı vzájemné energie interaguj´ıc´ıch ˇcástic. V grafu 4.2 jsou vyneseny závislosti sráˇzkového pr˚ uˇrezu na energii pro pruˇzné, excitaˇcn´ı a ionizaˇcn´ı sráˇzky elektron˚ u s neutrály argonu pˇri tlaku 133 Pa podle Bogaertsové [12]. Známe-li sráˇzkové pr˚ uˇrezy, mus´ıme naj´ıt vhodný algoritmus, kterým rozhodneme, ve kterém okamˇziku se konkrétn´ı ˇcástice sraz´ı. V oboru n´ızkých a stˇredn´ıch tlak˚ u se nejˇcastˇeji pouˇz´ıvá takzvaná ,,metoda nulové sráˇzky“ popsaná v [13]. Vycház´ı z u ´ vahy, ˇze náhodnou volnou dráhu ˇcástice m˚ uˇzeme metodou Monte Carlo generovat podle vzorce s = −λ log γ, (4.6) kde λ = 1/nσ je stˇredn´ı volná dráha ˇcástice, n je koncentrace terˇc˚ u, σ je sráˇzkový pr˚ uˇrez a γ je náhodné ˇc´ıslo z intervalu (0, 1) s rovnomˇerným rozdˇelen´ım, za pˇredpokladu, ˇze λ je konstantn´ı. Tento pˇredpoklad obvykle splnˇen nen´ı. Nicménˇe zavedeme-li umˇelý konstantn´ı sráˇzkový pr˚ uˇrez σ0 , m˚ uˇzeme generovat odpov´ıdaj´ıc´ı náhodné volné dráhy s s t´ım, ˇze pˇri kaˇzdé takové události metodou výbˇeru rozhodneme, ke kterému typu sráˇzky doˇslo. Moˇzné sráˇzky vˇsak mus´ıme rozˇs´ıˇrit o ,,nulovou sráˇzku“, pˇri které ˇcástice pokraˇcuje v pohybu, jako by k ˇzádné sráˇzce nedoˇslo. Sráˇzkový pr˚ uˇrez ,,nulové sráˇzky“ je pak právˇe takový, aby pro libovolnou energii platilo X σ0 = σi (ǫ), (4.7) i

kde i bˇeˇz´ı pˇres vˇsechny druhy sráˇzek vˇcetnˇe ,,nulové sráˇzky“. Metoda nulové sráˇzky je vhodná tam, kde stˇredn´ı je doba mezi sráˇzkami výraznˇe vˇetˇs´ı neˇz ˇcasový krok pohybu ˇcástice. Pro vyˇsˇs´ı tlaky pouˇz´ıváme odliˇsné


24

20 Pruˇzn´ a sr´ aˇzka Excitace Ionizace

σ [10−20 m2 ]

15

10

5

0

0

5

10

15 Er el [eV]

20

25

30

Obrázek 4.2: Sráˇzkové pr˚ uˇrezy pro interakci elektron˚ u s neutrály argonu metody, které jsou efektivnˇejˇs´ı a nevyˇzaduj´ı neˇzádouc´ı zkracován´ı ˇcasového kroku pohybu ˇcástice. V práci [14] je popsán model sráˇzek nezávislý na volbˇe kroku, který je pro vyˇsˇs´ı tlaky vhodnˇejˇs´ı. Tato metoda dále pˇredpokládá, ˇze terˇce jsou nehybné. V pˇr´ıpadˇe sráˇzek elektron˚ u s neutrály je tento pˇredpoklad pˇribliˇznˇe splnˇený, protoˇze rychlost elektron˚ u je mnohem vˇetˇs´ı neˇz rychlost neutrál˚ u. Sráˇzky iont˚ u s neutrály tento pˇredpoklad nesplˇ nuj´ı, protoˇze rychlosti ˇcástic jsou srovnatelné a terˇce tedy nelze povaˇzovat za nehybné. Moˇzné ˇreˇsen´ı tohoto problému je uvedeno v práci [14]. Rychlost terˇce se pˇri sráˇzce negeneruje na základˇe maxwellovského rozdˇelen´ı, ale z upraveného rozdˇelen´ı, které chybu kompenzuje. Výˇse popsaný ˇcásticový model oznaˇcujeme jako selfkonzistentn´ı, tj. s´ am sobˇe odpov´ıdaj´ıc´ı. Nejen v oblasti hybridn´ıho modelován´ı se vˇsak také setkáváme s neselfkonzistentn´ım ˇcásticovým modelem. V selfkonzistentn´ım modelu se ˇca´stice pohybuj´ı pod vlivem vnˇejˇs´ıch a vzájemných sil, které jsou za bˇehu modelu pr˚ ubˇeˇznˇe poˇc´ıtány. V neselfkozistentn´ım modelu jsou s´ıly pevnˇe dány napˇr´ıklad potenciálem elektrického pole a nabité ˇcástice se pohybuj´ı bez vzájemné interakce. Pˇrednosti neselfkonzistentn´ıho modelu vynikaj´ı pˇri modelován´ı v tˇr´ıdimenzionáln´ım prostoru. Obvyklý selfkozistentn´ı model je zat´ım sp´ıˇse za hranic´ı naˇsich výpoˇcetn´ıch moˇznost´ı. V neselfkonzistentn´ı modelu jsou s´ıly dány, a proto se vyhneme ˇcasovˇe nároˇcnému výpoˇctu sil. Nav´ıc nen´ı nutné uchovávat v pamˇeti infor-


25

maci o vˇsech ˇcástic´ıch najednou, protoˇze spolu v modelu neinteraguj´ı. Mezi dalˇs´ı výhody patˇr´ı snadná a efektivn´ı paralelizace výpoˇctu. Kaˇzdé vlákno poˇc´ıtá pohyb jedné ˇcástice nezávisle a vlákna se ovlivˇ nuj´ı pouze v okamˇziku pr˚ ubˇeˇzného ukládán´ı výsledk˚ u. Neˇz se rozhodneme pouˇz´ıt neselfkonzistentn´ı pˇr´ıstup, mus´ıme zváˇzit, zda je pro naˇsi u ´ lohu vhodný. Neselfkonzistentn´ı model lze pouˇz´ıt tam, kde nejsou podstatné s´ıly mezi jednotlivými ˇcásticemi. To v plazmatu m˚ uˇze být splnˇeno d´ıky st´ınˇen´ı. Pro pohyb d´ılˇc´ıch ˇcástic sledovaný ve vˇetˇs´ım mˇeˇr´ıtku jsou pak rozhoduj´ıc´ı prostorové náboje a vnoˇrené substráty. N´ıˇze uvád´ıme nˇekolik pˇr´ıklad˚ u pouˇzit´ı neselfkonzistentn´ıch model˚ u. Souˇcasnˇe také zodpov´ıme otázku, jak zjistit s´ıly p˚ usob´ıc´ı na ˇcástice. • Kombinace s analytick´ ym v´ ypoˇ ctem

V nˇekterých jednoduchých pˇr´ıpadech dokáˇzeme vypoˇc´ıtat p˚ usob´ıc´ı s´ıly analyticky na základˇe teorie nebo experimentáln´ıch poznatk˚ u. Tématu této práce se bl´ıˇz´ı napˇr´ıklad publikace [15], která se zabývá studiem pohybu elektron˚ u v okol´ı Langmuirovy sondy. Potenciál elektrického pole je jednoduˇse zadán aproximativn´ı funkc´ı v závislosti na vzdálenosti od nekoneˇcné válcové sondy. Pro kaˇzdou hodnotu pˇredpˇet´ı na sondˇe je sledován pohyb 10 000 ˇcástic, coˇz umoˇzn ˇ uje provést výpoˇcet v krátkém ˇcase pro velké mnoˇzstv´ı hodnot pˇredpˇet´ı. Zásadn´ı nevýhodou této metody je omezen´ı na pˇr´ıpady, kde jsme schopni s´ıly p˚ usob´ıc´ı na ˇcástice takto jednoduˇse popsat.

• Kombinace s metodou PIC

Neselfkonzistentn´ı model m˚ uˇzeme pouˇz´ıt k upˇresnˇen´ı výsledk˚ u selfkonzistentn´ıho modelu zaloˇzeného na metodˇe PIC. Poˇcet ˇcástic, se kterým pracuje PIC model, bývá obvykle menˇs´ı, neˇz jaký je potˇreba k urˇcen´ı napˇr´ıklad energetického rozdˇelen´ı ˇcástic dopadaj´ıc´ıch na sondu. Mezi moˇzná ˇreˇsen´ı patˇr´ı výpoˇcet pomoc´ı neselfkonzistentn´ıho modelu, ve kterém poˇc´ıtáme pohyb ˇcástic v potenciálu daném PIC modelem. Poˇcet ˇcástic v této fázi m˚ uˇze být výraznˇe vyˇsˇs´ı neˇz v PIC modelu bez zvýˇsených nárok˚ u na operaˇcn´ı pamˇet’. Postup je výhodný z hlediska poˇc´ıtaˇcového i fyzikáln´ıho, protoˇze v obou modelech pouˇz´ıváme mnoho shodných metod, napˇr´ıklad ˇreˇsen´ı pohybových rovnic, ˇreˇsen´ı sráˇzek s neutrály a okrajové podm´ınky.

• Kombinace se spojit´ ym modelem

Pˇredevˇs´ım v tˇr´ıdimenzionáln´ıch u ´ lohách je pˇredchoz´ı pˇr´ıklad stále pˇr´ıliˇs nároˇcný na výpoˇcetn´ı prostˇredky. Selfkonzistentn´ı PIC model m˚ uˇzeme nahradit spojitým modelem, který je sice výraznˇe ménˇe nároˇcný, ale také ménˇe


26

pˇresný. Neselfkonzistentn´ım modelem m˚ uˇzeme z´ıskat velmi detailn´ı výsledky, které ale odpov´ıdaj´ı nepˇresnému potenciálu elektrického pole. Nicménˇe tyto výsledky m˚ uˇzeme pouˇz´ıt k opravˇe p˚ uvodn´ıho spojitého modelu. Tato myˇslenka tvoˇr´ı základ iteraˇcn´ıho hybridn´ıho modelu, který bude n´ıˇze detailnˇeji popsán. Z tˇechto pˇr´ıklad˚ u je zˇrejmé, ˇze neselfkonzistentn´ı ˇcásticový model je sp´ıˇse pomocným prostˇredkem, který sám nem˚ uˇze plnˇe nahradit selfkonzistentn´ı výpoˇcet.

4.3

Spojit´ e modelov´ an´ı

Spojitý model nahrazuje diskrétn´ı ˇcástice spojitými veliˇcinami, které popisuj´ı jejich chován´ı jako tekutinu. M´ısto souˇradnic konkrétn´ıch ˇcástic, jejich energi´ı a podobnˇe zavád´ıme odpov´ıdaj´ıc´ı koncentraci, hustotu energie a pod. Prostorový a ˇcasový vývoj veliˇcin urˇcuj´ı odpov´ıdaj´ıc´ı parciáln´ı diferenciáln´ı rovnice: momenty Boltzmannovy kinetické rovnice (2.14), Maxwellovy rovnice, materiálové vztahy a podobnˇe. Minimáln´ı sada rovnic pro popis n´ızkoteplotn´ıho dvousloˇzkového plazmatu se skládá z následuj´ıc´ıch rovnic: • Rovnice kontinuity pro elektrony ∂ne + ∇ · J~e = re . ∂t

(4.8)

∂ni + ∇ · J~i = ri . ∂t

(4.9)

~ − De ∇ne . J~e = ne µe E

(4.10)

~ − Di ∇ni . J~i = ni µi E

(4.11)

e (ni − ne ). ε0

(4.12)

• Rovnice kontinuity pro ionty

• Tok elektron˚ u • Tok iont˚ u • Poissonova rovnice • Potenciál elektrického pole

∆ϕ = −

~ = −∇ϕ. E

(4.13)


27

V rovnic´ıch výˇse jsou ne a ni koncentrace elektron˚ u a iont˚ u, µe a µi jsou pohyblivosti (se znaménkem), De a Di jsou koeficienty difuze a re a ri pˇredstavuj´ı ˇcleny vzniku a zániku pˇr´ısluˇsných ˇcástic nepruˇznými sráˇzkami. Tyto parciáln´ı diferenciáln´ı rovnice s odpov´ıdaj´ıc´ımi okrajovými podm´ınkami lze ˇreˇsit r˚ uznými metodami. Pro jednoduˇsˇs´ı geometrická uspoˇrádán´ı lze pouˇz´ıt metodu koneˇcných diferenc´ı. Pro sloˇzitˇejˇs´ı geometrie a pˇredevˇs´ım ve tˇrech dimenz´ıch je vhodnˇejˇs´ı metoda koneˇcných prvk˚ u. Abychom mohli tuto soustavu numericky ˇreˇsit, potˇrebujeme znát zm´ınˇené koeficienty difuze D a pohyblivosti µ. M˚ uˇzeme je z´ıskat pˇr´ımo z experimentu, nebo je m˚ uˇzeme vypoˇc´ıtat z jiných parametr˚ u plazmatu. V [16] jsou pro elektrony tyto koeficienty odvozeny z Boltzmannovy kinetické rovnice vztahy D=

4π 3

Z

∞

0

c4 f0 (c)dc ν

(4.14)

a

q 4π Z ∞ c3 df0 (c) · dc, (4.15) m 3 0 ν dc kde q je naboj elektronu, m je hmotnost elektronu a f0 (c) je nenaruˇsená rozdˇelovac´ı funkce. Koneˇcnˇe ν je sráˇzková frekvence daná obecnˇe vztahem µ=

ν = ng cσ(c),

(4.16)

kde ng je koncentrace neutráln´ıch ˇcástic a σ(c) je u ´ˇcinný pr˚ uˇrez pruˇzné sráˇzky elektronu s neutráln´ımi ˇcásticemi. Je-li sráˇzková frekvence nezávislá na rychlosti a rozdˇelen´ı maxwellovské, snadno dostaneme vztahy D= a

kT mν

(4.17)

q (4.18) mν Tyto vztahy jsou velmi ˇcasto uˇz´ıvány, ale je nutné pamatovat na to, ˇze podm´ınky, za kterých byly odvozeny, jsou velmi pˇr´ısné. Je-li to moˇzné, je vhodnˇejˇs´ı pouˇz´ıt vztahy (4.14) a (4.15). Pro ionty budeme povaˇzovat koeficient difuze a pohyblivost za konstanty. V modelech totiˇz budeme pracovat s kladným pˇredpˇet´ım, d´ıky kterému se ionty budou vyskytovat pˇredevˇs´ım v oblastech slabˇe naruˇseného plazmatu. Kdybychom poˇzadovali pˇresnˇejˇs´ı popis difuze a driftu iont˚ u, lze v literatuˇre naj´ıt vhodné aproximativn´ı závislosti parametr˚ u na elektrickém poli, pro argon viz napˇr´ıklad [17]. Základn´ım nedostatkem spojitých model˚ u je jejich malá vypov´ıdac´ı schopnost. Výsledkem výpoˇctu jsou totiˇz závislosti veliˇcin jako koncentrace, stˇredn´ı energie a µ=


28

podobnˇe na poloze v ˇcase. Na rozd´ıl od ˇcásticového modelu nem˚ uˇzeme zjistit parametry jednotlivých ˇcástic. To je souˇcasnˇe d˚ uvod nepˇresnosti tˇechto výpoˇct˚ u. Obvykle mus´ıme totiˇz a priori pˇredpokládat rozdˇelovac´ı funkce ˇcástic, abychom mohli urˇcit koeficienty reakc´ı a difuze a pohyblivosti ˇcástic. Odtud se jiˇz dostáváme opˇet k myˇslence hybridn´ıho modelu, který nám m˚ uˇze poskytnout chybˇej´ıc´ı informace.

4.3.1

Numerick´ e metody ˇ reˇ sen´ı soustav parci´ aln´ıch diferenci´ aln´ıch rovnic

V následuj´ıc´ım textu struˇcnˇe uvedeme základn´ı metody, které se obvykle pouˇz´ıvaj´ı pro ˇreˇsen´ı parciáln´ıch diferenciáln´ıch rovnic spojitých model˚ u. Výˇse uvedenou soustavu uprav´ıme do vhodného tvaru dosazen´ım z (4.10) do (4.8) a z (4.11) do (4.9). Pro jednoduchost zde budeme povaˇzovat De , µe , Di a µi za konstanty a ˇcleny re a ri rovné nule. Výsledkem jsou následuj´ıc´ı rovnice ∂ne ~ − De ∆ne = 0, + ∇ · (µe ne E) (4.19) ∂t ∂ni ~ − Di ∆ni = 0. + ∇ · (µi ni E) (4.20) ∂t Do soustavy patˇr´ı rovnˇeˇz výˇse uvedená Poissonova rovnice. Rovnice (4.19) a (4.20) se liˇs´ı pouze v koeficientech a nazýváme je transportn´ımi rovnicemi. Naˇs´ım c´ılem ~ r , t), respektive ϕ(~r, t) pro zadané poˇcáteˇcn´ı je nalézt funkce ne (~r, t), ni (~r, t) a E(~ a okrajové podm´ınky. Pro jednoduché geometrické podm´ınky postaˇcuje metoda koneˇcných diferenc´ı. V pracovn´ı oblasti sestroj´ıme mˇr´ıˇz a spojité prostorové souˇradnice nahrad´ıme diskrétn´ımi indexy uzl˚ u mˇr´ıˇze. Obdobnˇe provedeme diskretizaci ˇcasu. Diferenciáln´ı operátory nahrad´ıme diskrétn´ımi výrazy, napˇr´ıklad ∂n nj − nji−1 ∼ i+1 , ∂x 2∆x

∂n nj+1 − nji ∼ i ∂t ∆t

a

∂2n nji−1 − 2nji + nji+1 ∼ , (4.21) ∂x2 ∆x2

kde i je index prostorové souˇradnice x a j je index ˇcasu. Parciáln´ı diferenciáln´ı rovnice tak pˇrevedeme na soustavu lineárn´ıch algebraických rovnic. Pouˇzit´ı tˇechto výraz˚ u vede k tzv. FTCS (forward time, centered space) schématu — prostorové diference urˇcujeme z hodnot v j-tém ˇcasovém kroku prostorovˇe symetricky, ˇcasové diference z krok˚ u j a j + 1. Pˇrestoˇze toto schéma vede ke snadnému algoritmu, jeho pouˇzit´ı je omezené velmi nepˇr´ıznivými podm´ınkami konvergence, podrobnˇeji viz [18]. Z hlediska konvergence jsou pro transportn´ı rovnici vhodnˇejˇs´ı implicitn´ı schémata. V tˇechto schématech se obecnˇe v diskretizaci prostorových derivac´ı objevuj´ı hodnoty v ˇcasovém kroku j + 1. Jako pˇr´ıklad uvád´ıme Crankovo-Nicolsonovo


29

schéma podle [18] ∂n 1 ∼ ∂x 2

nji+1 − nji−1 nj+1 − nj+1 i−1 + i+1 , 2∆x 2∆x !

(4.22)

j+1 1 nji−1 − 2nji + nji+1 nj+1 + nj+1 ∂2n i−1 − 2ni i+1 ∼ + . (4.23) ∂x2 2 ∆x2 ∆x2 Diskretizace ˇcasové derivace je shodná s pˇredchoz´ım pˇr´ıpadem. Nevýhodou implicitn´ıch schémat je výraznˇe vyˇsˇs´ı výpoˇcetn´ı nároˇcnost. V kaˇzdém kroku ˇreˇs´ıme velkou soustavu lineárn´ıch algebraických rovnic, jej´ıˇz matice je ˇr´ıdká. V souˇcasné dobˇe jsou ovˇsem k dispozici kvalitn´ı knihovny pro manipulaci s ˇr´ıdkými maticemi, napˇr´ıklad jiˇz zmiˇ novaná knihovna UMFPACK. ˇ Rada autor˚ u pouˇz´ıvá takzvaný Scharfetter˚ uv-Gummel˚ uv algoritmus, jehoˇz odvozen´ı a pouˇzit´ı v jednorozmˇerném modelu je struˇcnˇe popsáno v [19]. Dvoudimenzionáln´ı válcová varianta tohoto algoritmu je podrobnˇe popsána v [20]. ˇ sen´ı transportn´ı rovnice pˇredstavuje pouze ˇcást výpoˇctu. V praxi ˇreˇs´ıme Reˇ soustavu minimálnˇe tˇr´ı rovnic, která zahrnuje transportn´ı rovnice pro elektrony a ionty a Poissonovu rovnici. V kaˇzdém ˇcasovém kroku ˇreˇsen´ım transportn´ıch rovnic z´ıskáme nové hodnoty koncentrac´ı elektron˚ u a iont˚ u, na základˇe kterých ˇreˇs´ıme Poissonovu rovnici. Z´ıskaný potenciál elektrického pole v dalˇs´ım kroku vstupuje do transportn´ıch rovnic. C´ılem této práce je tvorba plnˇe tˇr´ıdimenzionáln´ıch model˚ u. V takovém pˇr´ıpadˇe je velmi vhodná metoda koneˇcných prvk˚ u (Finite Element Method, FEM ). Na rozd´ıl od metody koneˇcných diferenc´ı, která snadno popisuje problémy s jednoduchou geometri´ı, metoda koneˇcných prvk˚ u dokáˇze pracovat na velmi sloˇzité s´ıti, kterou lze snadno pˇrizp˚ usobit i jinak protich˚ udným poˇzadavk˚ um. Modelujemeli napˇr´ıklad intarakci relativnˇe malé sondy s plazmatem, poˇzadujeme souˇcasnˇe jemnou s´ıt’ v okol´ı sondy a velkou pracovn´ı oblast. Kdybychom chtˇeli pouˇz´ıt rovnomˇernou pravo´ uhlou s´ıt’, na jakou jsme zvykl´ı z metody koneˇcných diferenc´ı, problém by byl s dostupnými výpoˇcetn´ımi prostˇredky témˇeˇr neˇreˇsitelný. V metodˇe konˇcných prvk˚ u vytvoˇr´ıme v okol´ı geometricky sloˇzité sondy velmi jemnou s´ıt’ a ve vzdálenˇejˇs´ıch oblastech ji ponecháme výraznˇe hrubˇs´ı. T´ımto zp˚ usobem lze výpoˇcet podstatnˇe zefektivnit. Pˇr´ıklad dvoudimenzionáln´ı s´ıtˇe pro okol´ı válcové sondy je uveden na obrázku 4.3. V této práci se nebudeme podrobnˇe zabývat teoretickými detaily metody koneˇcných prvk˚ u. Pro ˇreˇsen´ı konkrétn´ıch problém˚ u budeme pouˇz´ıvat komerˇcn´ı software COMSOL Multiphysics, pˇresto povaˇzujeme za vhodné alespoˇ n struˇcnˇe uvést základn´ı myˇslenky této metody. Postupovat budeme podle [21]. Zvol´ıme si obecnou u ´ lohu na oblasti Ω, která se ˇcasto nazývá doména, s hranic´ı ∂Ω. V oblasti Ω poˇzadujeme splnˇen´ı diferenciáln´ı rovnice

!

a(u) = f,

(4.24)


30

Obrázek 4.3: Pˇr´ıklad 2D s´ıtˇe v okol´ı nekoneˇcné válcové sondy pro metodu koneˇcných prvk˚ u kde u je hledaná funkce. Na hranici má ˇreˇsen´ı splˇ novat Dirichletovu hraniˇcn´ı podm´ınku b0 (u) = g0 (4.25) nebo Neumannovu podm´ınku b1 (u) = g1 .

(4.26)

Na rozd´ıl od metody koneˇcných diferenc´ı, ve které hledáme pˇribliˇzné ˇreˇsen´ı jen v uzlových bodech mˇr´ıˇze, výsledkem metody koneˇcných prvk˚ u je aproximativn´ı funkce uˆ, kterou lze vyjádˇrit lineárn´ı kombinac´ı b´ azových funkc´ı φk uˆ =

N X

φk u k .

(4.27)

k=1

Vzhledem k tomu, ˇze uˆ nen´ı pˇresné ˇreˇsen´ı, definujeme podle zadán´ı u ´ lohy rezidua rΩ = a(ˆ u) − f,

r0 = b0 (ˆ u) − g0

a r1 = b1 (ˆ u) − g1 .

(4.28)


31

Aproximativn´ı ˇreˇsen´ı uˆ hledáme takové, aby rezidua byla v urˇcitém smyslu minimáln´ı. Toho dosáhneme splnˇen´ım nulovosti souˇctu integrál˚ u Z

Ω

wl rΩ dΩ +

Z

∂Ω0

wl0r0 d∂Ω

+

Z

∂Ω1

wl1 r1 d∂Ω = 0,

(4.29)

kde wl jsou váhové funkce. Rovnice (4.29) se nazývá slab´ a formulace. Soustava rovnic (4.24), (4.25) a (4.26) se analogicky nazývá siln´ a formulace. Slabá formulace nen´ı ekvivalentn´ı silné formulaci. Kaˇzdé ˇreˇsen´ı silné formulace je souˇcasnˇe ˇreˇsen´ım slabé formulace, nicménˇe slabou formulaci splˇ nuj´ı i dalˇs´ı ˇreˇsen´ı. Typickým rozd´ılem formulac´ı jsou r˚ uzné poˇzadavky na spojitost ˇreˇsen´ı a jeho derivac´ı. Zat´ım jsme nespecifikovali, jakou konkrétn´ı podobu maj´ı bázové funkce φk a váhy wl . Velmi vhodné jsou funkce s kompaktn´ım nosiˇcem, které jsou nenulové jen nad malým poˇctem element˚ u mˇr´ıˇze. Integrály v rovnici (4.29) pak m˚ uˇzeme pˇrepsat jako sumu integrál˚ u pˇres jednotlivé elementy. Bázové funkce vol´ıme jako interpolaˇcn´ı funkce v rámci element˚ u — nejˇcastˇeji lineárn´ı nebo kvadratické. Váhy wl jsou obvykle shodné s bázovými funkcemi φk (tzv. Galerkinova formulace), ale nen´ı to nutné. Zaveden´ım tˇechto funkc´ı pˇrevedeme rovnici (4.29) na soustavu lineárn´ıch rovnic pro koeficienty uk . Konkrétn´ı metodu ˇreˇsen´ı soustavy je potˇreba volit s ohledem na poˇzadovanou pˇresnost, rychlost a dostupné prostˇredky. COMSOL nab´ız´ı mimo jiné metody LU dekompozice (UMFPACK, SPOOLES), pˇr´ımé ˇreˇsiˇce (PARADISO), iterativn´ı metody (GMRS, konjugované gradienty, BiCGStab). Na základˇe koeficient˚ u uk podle vztahu (4.27) z´ıskáme aproximativn´ı ˇreˇsen´ı u ´lohy uˆ. f (ǫ)

Spojitý model

0

ˇ asticový model C´

ǫ0

ǫ

Obrázek 4.4: Schéma rozdˇelen´ı elektron˚ u na ,,rychlé“ (nad ǫ0 ) a ,,pomalé“ (pod ǫ0 )


4.4

32

Obvykl´ e hybridn´ı modely

Moˇznost´ı, jak kombinovat metody modelován´ı plazmatu, je obrovské mnoˇzstv´ı. Základn´ı pˇrehled je uveden napˇr´ıklad v [22]. V této práci bude provedeno srovnán´ı tˇr´ı hybridn´ıch model˚ u, které oznaˇc´ıme názvy iteraˇcn´ı, energetický a prostorový model .

4.4.1

Energetick´ y model

Obrázek 4.5: Zjednoduˇsené schéma energetického modelu Touto metodou se zabývá napˇr´ıklad [19]. Hranice mezi ˇcásticovým a spojitým modelem je dána nˇejakou významnou hodnotou energie elektronu. Napˇr´ıklad v pozitivn´ım sloupci doutnavého výboje v argonu je to nejniˇzˇs´ı energie elektronu potˇrebná k excitaci atomu argonu v základn´ım stavu. Elektrony s vyˇsˇs´ı energi´ı, neˇz je tato prahová energie, podstupuj´ı nepruˇzné excitaˇcn´ı a ionizaˇcn´ı sráˇzky, coˇz naruˇsuje elektronovou rozdˇelovac´ı funkci, která je pro niˇzˇs´ı energie napˇr´ıklad max-


33

wellovská. Elektrony tedy rozdˇel´ıme na ,,pomalé“ a ,,rychlé“ podle jejich energie. Pro pˇresnˇe maxwellovskou rozdˇelovac´ı funkci je situace znázornˇena na obrázku 4.4. Tato metoda je velmi výhodná, pokud v pracovn´ı oblasti jen velmi málo elektron˚ u pˇresahuje prahovou energii danou vlastnostmi neutrál˚ u. Nevýhodou je nutnost pˇredpokládat rozdˇelovac´ı funkci pomalých elektron˚ u, coˇz je problematické napˇr´ıklad v bl´ızkosti sondy na relativnˇe vysokém potenciálu v˚ uˇci teplotˇe elektron˚ u. Elektrony pˇri svém pohybu mohou vlivem elektrického pole a sráˇzek zvyˇsovat i sniˇzovat svou energii. Pokud energie ,,rychlého“ elektronu klesne pod prahovou energii, je pˇresunut do spojité ˇcásti modelu. Model vˇsak mus´ı zahrnovat také opaˇcný proces — ,,pomalé“ elektrony mohou zvýˇsit svou energii nad prahovou energii. V takovém pˇr´ıpadˇe je potˇreba tento elektron zaˇradit do ˇcásticové ˇcásti modelu. Nicménˇe ze spojité ˇcásti modelu z principu neznáme parametry jednotlivých elektron˚ u, a proto vyuˇz´ıváme pravdˇepodobnostn´ıho pˇr´ıstupu. Na základˇe velikosti elektrického pole, které pˇredevˇs´ım zp˚ usobuje urychlován´ı elektron˚ u, urˇc´ıme poˇcet elektron˚ u, které pˇrejdou do ˇcásticového modelu, a jejich konkrétn´ı parametry rozehrajeme metodou Monte Carlo. V pr˚ ubˇehu výpoˇctu, který je uveden na obrázku 4.5, se v kaˇzdém ˇcasovém kroku stˇr´ıdá ˇcásticový a spojitý model. V kaˇzdém kroku je dále potˇreba provést výˇse uvedenou výmˇenu ˇcástic mezi modely podle jejich energie. Jsou-li ˇcasové zmˇeny veliˇcin spojitého modelu dostateˇcnˇe malé, lze spojitou ˇcást poˇc´ıtat s delˇs´ım ˇcasovým krokem s t´ım, ˇze na jeden krok spojitého modelu pˇripadne nˇekolik krok˚ u ˇcásticového modelu podle pomˇeru délek ˇcasových krok˚ u.

4.4.2

Prostorov´ y model E(r )

ˇ asticový model C´

Spojitý model

0

rz

l

r

Obrázek 4.6: Rozdˇelen´ı pracovn´ı oblasti prostorového hybridn´ıho modelu


34

Obrázek 4.7: Zjednoduˇsené schéma prostorového modelu V ˇradˇe ˇreˇsených problém˚ u v n´ızkoteplotn´ım plazmatu pozorujeme, ˇze jen relativnˇe malá ˇcást zkoumané oblasti vyˇzaduje výpoˇcet ˇcásticovým modelem. Nab´ız´ı se tedy myˇslenka rozdˇelit pracovn´ı oblast na dvˇe ˇcásti. Tam, kde poˇzadujeme pˇresnˇejˇs´ı výpoˇcet, pouˇzijeme ˇcásticový model a zbytek oblasti pop´ıˇseme spojitým modelem. Pokud ˇcásticový model pracuje na malé oblasti, m˚ uˇzeme výpoˇcet t´ımto zp˚ usobem výraznˇe urychlit. Jako pˇr´ıklad nám poslouˇz´ı válcová Langmuirova sonda. Na obrázku 4.6 je nakresleno schéma modelu interakce sondy s plazmatem ve válcové geometrii 1d3v (tj. prostorová souˇradnice r a sloˇzky rychlosti vr , vϕ a vz ). Polomˇer válcové pracovn´ı oblasti je l a sonda o polomˇeru rp má osu shodnou s osou souˇradného systému. Na obrázku je ˇcervenˇe vynesena intenzita elektrického pole v okol´ı sondy. Dále od urˇcité vzdálenosti od poˇcátku menˇs´ı neˇz rz je elektrické pole velmi malé a jen velmi málo naruˇsuje plazma. Naopak bl´ıˇze k sondˇe se nacház´ı st´ın´ıc´ı oblast, ve které je plazma silnˇe naruˇsené. V praxi se ukazuje, ˇze z d˚ uvodu stability výpoˇctu je vhodné um´ıstit hranici mezi ˇcásticovým modelem (r < rz ) a spojitým modelem (r > rz ) dostateˇcnˇe daleko od výraznˇe naruˇsené oblasti.


35

Na obrázku 4.7 je schematicky uveden postup výpoˇctu prostorového modelu. V kaˇzdém ˇcasovém kroku provedeme spojitý i ˇcásticový výpoˇcet v odpov´ıdaj´ıc´ıch oblastech a zajist´ıme výmˇenu ˇcástic mezi oblastmi. Stejnˇe jako v energetickém hybridn´ım modelu nemus´ı být délka ˇcasového kroku stejná pro obˇe ˇcásti modelu.

4.4.3

Iteraˇ cn´ı model

Tato metoda je podrobnˇe popsána v [23], dˇr´ıve také napˇr´ıklad v [24]. Základn´ı schéma bˇehu iteraˇcn´ıho modelu je uvedeno na obrázku 4.8. Nejprve provedeme

Obrázek 4.8: Zjednoduˇsené schéma iteraˇcn´ıho modelu (podle [23]) výpoˇcet jednoduchým spojitým modelem, který je popsán v kapitole 4.3. Výsledné prostorové rozdˇelen´ı potenciálu elektrického pole pˇrevezmeme do neselfkonzistentn´ıho ˇcásticového modelu. Nyn´ı necháme elektrony (pˇr´ıpadnˇe ionty) pohybovat se v pevnˇe daném elektrickém poli a na základˇe analýzy jejich pohybu upˇresn´ıme koeficienty v rovnic´ıch spojitého modelu. Pracovn´ı oblast m˚ uˇzeme vhodnˇe rozdˇelit na buˇ nky, ve kterých zjist´ıme rozdˇelovac´ı funkci. Rozdˇelovac´ı funkci totiˇz potˇrebujeme znát pro výpoˇcet koeficientu difuze a pohyblivosti (vzorce (4.14) a (4.15)). Tyto upˇresnˇené koeficienty následnˇe vstupuj´ı do spojitého modelu. Tento proces opakujeme, aˇz dosáhneme ustáleného stavu. Tento postup je vhodný pˇredevˇs´ım pro hledán´ı ustáleného stavu. Pˇr´ıkladem takového problému je interakce substrátu s plazmatem v pozitivn´ım sloupci dout-


36

naváho výboje. Metoda nen´ı vhodná pro ˇreˇsen´ı problém˚ u s výraznými zmˇenami v ˇcase, napˇr´ıklad iontová implantace skokovou zmˇenou napˇet´ı.

4.4.4

Dalˇ s´ı hybridn´ı modely

V literatuˇre se lze setkat s ˇradou variant hybridn´ıch model˚ u, které ˇcasto vznikaj´ı jako kombinace výˇse popsaných model˚ u. Napˇr´ıklad v [25] je popsán model argonového výboje s pˇr´ımˇes´ı vod´ıku. Celkem tˇrináct druh˚ u ˇcástic je modelováno pˇreváˇznˇe energeticky dˇeleným modelem a nalezneme zde nav´ıc mimo jiné spojitý model pro neutráln´ı ˇcástice H a H2 , který zapoˇc´ıtává produkci a zánik neutrál˚ u v reakc´ıch a dif´ uzi.

Kapitola 5 C´ıle diplomov´ e pr´ ace C´ılem této práce je rozvoj hybridn´ıch metod modelován´ı ve fyzice plazmatu a jejich aplikace na konkrétn´ı problémy. Tyto modely nacházej´ı uplatnˇen´ı tam, kde výpoˇcty ˇcistˇe ˇcásticovými modely maj´ı pˇr´ıliˇs velké nároky na hardware a dobu výpoˇctu a pomˇernˇe nenároˇcné spojité modely nepodávaj´ı dostateˇcnˇe detailn´ı a pˇresné výsledky. Problémy interakce plazmatu s pevnou látkou, které vyˇzaduj´ı plnˇe tˇr´ıdimenzionáln´ı pˇr´ıstup, jsou typickým pˇr´ıkladem. Vzhledem k tomu, ˇze existuje ˇrada pˇr´ıstup˚ u k hybridn´ımu modelován´ı, vyzkouˇs´ıme nejprve nˇekolik nejˇcastˇejˇs´ıch variant a porovnáme je se základn´ım ˇcásticovým a spojitým modelem. Budeme diskutovat jejich pˇrednosti a nevýhody na základˇe aplikace na zjednoduˇsený problém nekoneˇcné válcové sondy. Vybranou metodu hybridn´ıho modelován´ı rozpracujeme dále tak, abychom mohli ˇreˇsit plnˇe tˇr´ıdimenzionáln´ı u ´ lohy. Hlavn´ı principy modelu pˇredstav´ıme ve dvoudimenzionáln´ı geometrii vzhledem k vˇetˇs´ı názornosti obrázk˚ u a graf˚ u. Výsledným modelem budou ˇreˇseny dva problémy — koneˇcná válcová sonda v plazmatu s driftem a interakce plazmatu s vodivým substrátem s nerovnostmi.

37

Kapitola 6 Srovn´ an´ı vybran´ ych hybridn´ıch model˚ u Typické hybridn´ı modely, které byly pˇredstaveny v kapitole 4, zrealizujeme a porovnáme jejich výsledky a efektivitu. Do srovnán´ı zahrneme také ˇcásticový a spojitý model. Abychom se mohli soustˇredit pˇredevˇs´ım na metody hybridn´ıho modelován´ı, budeme ˇreˇsit pomˇernˇe jednoduchý problém — nekoneˇcnˇe dlouhou válcovou sondu vnoˇrenou do n´ızkoteplotn´ıho argonového plazmatu. Tato u ´ loha je v prostoru jednorozmˇerná, coˇz výraznˇe usnadˇ nuje poˇc´ıtaˇcové ˇreˇsen´ı spojité ˇcásti model˚ u.

6.1

Energetick´ y model

Pr˚ ubˇeh výpoˇctu energetickým modelem je uveden ve schématu 4.5. Podle potˇreby si stanov´ıme hraniˇcn´ı energii mezi spojitou a ˇcásticovou ˇcást´ı modelu, viz obrázek 4.4. Model m˚ uˇzeme rozdˇelit do ˇctyˇr blok˚ u: 1. ˇreˇsiˇc Poissonovy rovnice (4.12), 2. spojitý ˇreˇsiˇc pro ,,pomalé“ elektrony, 3. ˇcásticový ˇreˇsiˇc pro ,,rychlé“ elektrony, 4. spojitý ˇreˇsiˇc pro ionty. ˇ siˇc Poissonovy rovnice na základˇe prostorové hustoty náboje a okrajových Reˇ podm´ınek vypoˇc´ıtá potenciál elektrického pole v uzlových bodech mˇr´ıˇze. V jednodimenzionáln´ım pˇr´ıpadˇe pouˇz´ıváme Thomas˚ uv algoritmus [9]. Spojitý ˇreˇsiˇc pro ,,pomalé“ elektrony ˇreˇs´ı transportn´ı rovnici ∂ne + ∇ · (−ne µe ∇ϕ − De ∇ne ) = 0, ∂t 38

(6.1)

´ Í VYBRANYCH ´ ˚ KAPITOLA 6. SROVNAN HYBRIDNÍCH MODELU

39

která plyne z rovnic (4.8) a (4.10). Metodou koneˇcných diferenc´ı tuto rovnici pˇrevedeme na soustavu lineárn´ıch algebraických rovnic, které opˇet ˇreˇs´ıme Thomasovým algoritmem. Analogicky postupojeme v pˇr´ıpadˇe spojitého ˇreˇsiˇce pro ionty. ˇ asticový ˇreˇsiˇc je zaloˇzen na metodách Monte Carlo a molekulárn´ı dynamiky. C´ Pohyb ,,rychlých“ elektron˚ u je dán elektrickým polem, jehoˇz potenciál poˇc´ıtá ˇreˇsiˇc ´ cinné pr˚ Poissonovy rovnice, a sráˇzkami s atomy neutráln´ıho pozad´ı. Uˇ uˇrezy tˇechto sráˇzek podle [12] jsou uvedeny na obrázku 4.2. Po startu modelu nastav´ıme poˇcáteˇcn´ı podm´ınky. Spojitým ˇreˇsiˇcem pro ionty a pro ,,pomalé“ elektrony vypoˇc´ıtáme zmˇenu koncentrace iot˚ u a ,,pomalých“ elekˇ asticovým ˇreˇsiˇcem vypoˇc´ıtáme nové polohy ,,rychlých“ elektron˚ tron˚ u za ˇcas ∆t. C´ u za ˇcas ∆t. Následuje výmˇena ˇcástic mezi spojitým a ˇcásticovým ˇreˇsiˇcem elektron˚ u, pro kterou jsou uˇziteˇcné následuj´ıc´ı poznatky. Pro zjednoduˇsen´ı zavedeme bezrozmˇernou energii E . kT

ǫ=

(6.2)

Maxwellovské rozdˇelen´ı pro tuto bezrozmˇernou energii má tvar 2 √ f (ǫ)dǫ = √ ǫ exp(−ǫ)dǫ. π

(6.3)

Energie ǫ0 tvoˇr´ı zm´ınˇenou hranici. Na obrázku 4.4 vid´ıme, ˇze ,,rychlých“ elektron˚ u je obvykle výraznˇe ménˇe neˇz ,,pomalých“. Souˇcasnˇe pro popis ,,pomalých“ elektron˚ u postaˇcuje spojitý model poˇc´ıtaj´ıc´ı s maxwellovským rozdˇelen´ım. Naopak ,,rychlé“ elektrony je nutné modelovat ˇcásticovˇe, protoˇze vstupuj´ı do nepruˇzných sráˇzek, které výraznˇe mˇen´ı jejich energii. Pro realizaci modelu této kategorie potˇrebujeme nˇekolik vzorc˚ u, pomoc´ı kterých generujeme ,,rychlé“ elektrony na poˇcátku i v pr˚ ubˇehu výpoˇctu. Koncentrace ˇcástic, které maj´ı vyˇsˇs´ı energii neˇz ǫ0 = E0 /kT , je Γ(1,5, ǫ0 ) 2 , n(ǫ0 ) = n0 √ Γ(1,5, ǫ0 ) = n0 π Γ(1,5)

(6.4)

kde Γ(a, x) je horn´ı ne´ uplná gamma funkce definovaná Γ(a, x) =

Z

∞

x

ta−1 e−t dt.

(6.5)

V programech byla pouˇzita funkce gamma q z knihovny Boost [26] definovaná gamma q(a, x) =

Γ(a, x) . Γ(a)

(6.6)


40

V metodˇe Monte Carlo potˇrebujeme rozehrávat energie ,,rychlých“ elektron˚ u. Vyuˇzijeme k tomu inverzn´ı ne´ uplnou horn´ı gamma funkci. Snadno generujeme náhodná ˇc´ısla ξ s rovnomˇerným rozdˇelen´ım na intervalu (0, ξ0 ), kde ξ0 =

Γ(1,5, ǫ0 ) . Γ(1,5)

(6.7)

Energie elektronu odpov´ıdaj´ıc´ı danému ξ je dána inverzn´ı ne´ uplnou horn´ı gamma funkc´ı ǫ = Γ−1 (1,5, ξ · Γ(1,5)). (6.8) Tuto funkci poskytuje rovnˇeˇz knihovna Boost. Vygenerované energie pˇriˇrad´ıme jednotlivým elektron˚ um a na poˇcátku vol´ıme jejich smˇer náhodnˇe. V pr˚ ubˇehu simulace vyˇrazujeme z ˇcásticového modelu elektrony, jejichˇz energie klesla pod ǫ0 . Nutné je zohlednit také opaˇcný proces — ,,pomalé“ elektrony z´ıskávaj´ı energii v elektrickém poli a jejich ˇcást pˇrecház´ı do ˇcásticového modelu. Elektrické pole v bl´ızkosti sondy naruˇsuje energetické rozdˇelen´ı elektron˚ u. Pˇr´ımo z modelu vˇsak m˚ uˇzeme zjistit pouze ˇcást rozdˇelovac´ı funkce — pro ǫ > ǫ0 . Pro menˇs´ı energie mus´ıme rozdˇelovac´ı funkci zadat. Nejjednoduˇsˇs´ı moˇznost´ı je uvaˇzovat pro ,,pomalé“ elektrony lokáln´ı termodynamickou rovnováhu (local thermodynamic equilibrium, LTE). V takovém pˇr´ıpadˇe maj´ı ,,pomalé“ elektrony maxwellovské rozdˇelen´ı a jejich teplota vlivem elektrického pole roste. V podm´ınkách typických pro n´ızké aˇz stˇredn´ı tlaky, zvláˇstˇe v oblastech silného elektrického pole (napˇr. v okol´ı sond), elektrony lokáln´ı rovnováhy zcela nedosahuj´ı, protoˇze jejich stˇredn´ı volná dráha je pˇr´ıliˇs dlouhá. Dále je vhodné pˇrej´ıt od LTE k ˇcásteˇcné lokáln´ı termodynamické rovnováze (partial LTE, PLTE , podrobnosti napˇr´ıklad v [27]), která zahrnuje v bilanˇcn´ı rovnici energie v´ıce zdroj˚ u. Nyn´ı známe koncentraci iont˚ u a celkovou koncentraci elektron˚ u po uplynut´ı ˇcasového kroku. Z koncentrac´ı vypoˇc´ıtáme prostorovou hustotu náboje, která vstupuje do ˇreˇsiˇce Poissonovy rovnice. Jej´ım vyˇreˇsen´ım z´ıskáme potenciál elektrického pole a postupujeme do dalˇs´ıho ˇcasového kroku. Opakován´ım tohoto postupu vypoˇc´ıtáme pr˚ ubˇeh parametr˚ u plazmatu v ˇcase podle poˇzadavk˚ u zadán´ı.

6.2

Prostorov´ y model

Základn´ı bloky prostorového modelu jsou podobné jako v pˇr´ıpadˇe energetického modelu: 1. ˇreˇsiˇc Poissonovy rovnice (4.12), 2. spojitý ˇreˇsiˇc pro elektrony,


41

3. ˇcásticový ˇreˇsiˇc pro elektrony, 4. spojitý ˇreˇsiˇc pro ionty. Spojitý a ˇcásticový ˇreˇsiˇc pro elektrony pracuj´ı na r˚ uzných ˇcástech pracovn´ı oblasti — ˇcásticový ˇreˇsiˇc v bl´ızkosti sondy (r < rz na obrázku 4.6) a spojitý ˇreˇsiˇc ve zbytku oblasti (r > rz ). V bodˇe rz mus´ım provést navázán´ı spojité a ˇcásticové ˇcásti. Pˇrechod z ˇcásticové ˇcásti do spojité ˇreˇs´ıme tak, ˇze koncentraci elektron˚ u ˇcásticového modelu v malém intervalu (rz −dr, rz ) pouˇzijeme jako Dirichletovu okrajovou podm´ınku spojité ˇcásti. Délku intervalu dr zvol´ıme shodnou s elementem osy r. Pˇr´ıpadnˇe je moˇzné provést extrapolaci do bodu rz na základˇe pr˚ ubˇehu koncentrace ve v´ıce bodech ˇcásticové ˇcásti modelu bl´ızko rz . Souˇcasnˇe mus´ıme zajistit pˇrechod ˇcástic z oblasti spojitého modelu do ˇcásticového. Ze spojitého modelu m˚ uˇzeme urˇcit pouze koncentraci elektron˚ u, ale neznáme jejich rychlostn´ı rozdˇelen´ı. Prvn´ı pokusy zaloˇzené na pˇredpokladu maxwellovského rozdˇelen´ı byly ne´ uspˇeˇsné, protoˇze rozdˇelen´ı rychlost´ı je naruˇsené elektrickým polem. Rychlosti ˇcástic vstupuj´ıc´ıch do ˇcásticového modelu proto generujeme pomocným neselfkonzistentn´ım modelem, který je virtuálnˇe um´ıstˇen vpravo od bodu rz — tedy vnˇe ˇcásticové ˇcásti hybridn´ıho modelu. V pomocném neselfkonzistentn´ım modelu poˇc´ıtáme pohyb elektron˚ u pod vlivem elektrického pole daného spojitým modelem a sráˇzek s neutrály. Okrajové podm´ınky jsou periodické a kopie ˇcástic souˇcasnˇe vstupuj´ı do ˇcásticové ˇcásti hybridn´ıho modelu. Takto zaruˇc´ıme velmi dobré navázán´ı spojitého a ˇcásticového modelu. Výpoˇcet pohybu elektron˚ u je opˇet proveden metodami molekulárn´ı dynamiky a Monte Carlo. Elektrická s´ıla p˚ usob´ıc´ı na elektrony je z´ıskávána z rozdˇelen´ı potenciálu elektrického pole, který je poˇc´ıtán ˇreˇsiˇcem Poissonovy rovnice, který pracuje stejným zp˚ usobem jako v energetickém modelu. Spojitý ˇreˇsiˇc pro ionty je shodný s energetickým modelem. Také v tomto pˇr´ıpadˇe postupujeme po jednotlivých ˇcasových kroc´ıch poˇzadovaného ˇcasového intervalu.

6.3

Iteraˇ cn´ı model

Iteraˇcn´ı model je se skládá z následuj´ıc´ıch blok˚ u: 1. standardn´ı spojitý model, 2. neselfkonzistentn´ı ˇcásticový model elektron˚ u. Standardn´ı spojitý model byl popsán v kapitole 4.3. Rovnice (4.8) aˇz (4.11) byly po diskretizaci implicitn´ım schématem ˇreˇseny Thomasovým algoritmem. Po-


42

issonova rovnice byla diskretizována podle schématu (4.21) a rovnˇeˇz ˇreˇsena Thomasovým algoritmem. V neselfkonzistentn´ım ˇcásticovém modelu poˇc´ıtáme pohyb elektron˚ u v elektrickém poli, které z´ıskáváme ze spojitého modelu. Elektrony vpouˇst´ıme do pracovn´ı oblasti ,,po jednom“, a proto má výpoˇcet niˇzˇs´ı nároky na operaˇcn´ı pamˇet’ neˇz ˇcásticový selfkonzistentn´ı model. Elektrony se v pracovn´ı oblasti sráˇz´ı s neutráln´ımi atomy argonu. Uvaˇzujeme sráˇzky pruˇzné, excitaˇcn´ı a ionizaˇcn´ı, jejichˇz sráˇzkové pr˚ uˇrezy jsou uvedeny v grafu 4.2. Na základˇe vzorkován´ı rychlost´ı elektron˚ u v buˇ nkách mˇr´ıˇze postupnˇe sestavujeme rozdˇelovac´ı funkci pro kaˇzdou buˇ nku. Z rozdˇelovac´ı funkce na konci neselfkonzistentn´ıho výpoˇctu integracemi (4.14) a (4.15) z´ıskáme nové koeficienty difuze a pohyblivosti pro kaˇzdou buˇ nku samostatnˇe. Nové koeficienty pˇredáváme zpˇet do spojitého modelu a poˇc´ıtáme vývoj koncentrac´ı ˇcástic a potenciálu v dalˇs´ım ˇcasovém intervalu, na jehoˇz konci pˇredáme potenciál elektrického pole do dalˇs´ı iterace neselfkonzistentn´ıho modelu. Takto postupujeme do konce celkového ˇcasového intervalu.

6.4

Krit´ eria srovn´ an´ı

Zde srovnávané hybridn´ı modely se od sebe podstatnˇe liˇs´ı, a proto je obt´ıˇzné naj´ıt kritéria, podle kterých bychom je mohli srovnat. Nejjednoduˇsˇs´ı a asi nejˇzádanˇejˇs´ı je srovnán´ı celkové doby bˇehu programu. Pro toto srovnán´ı byly pouˇzity k ˇreˇsen´ı odpov´ıdaj´ıc´ıch problém˚ u stejné algoritmy. Mnohé ˇcásti program˚ u jsou shodné pro vˇsechny modely. Vˇsechny programy byly zámˇernˇe napsané bez pouˇzit´ı paralelizace, protoˇze moˇznosti paralelizace povaˇzujeme za samostatné kritérium. Velmi d˚ uleˇzité je kritérium pˇresnosti výsledk˚ u. Bohuˇzel výsledky tˇechto model˚ u nelze pˇr´ımo srovnat s experimentem. Vhodnou náhradou je srovnán´ı s klasickým ˇcásticovým selfkonzistentn´ım modelem. Srovnáme pˇredevˇs´ım pr˚ ubˇeh potenciálu elektrického pole v okol´ı válcové sondy. Obt´ıˇzné je hodnotit vhodnost metody pro daný problém. Kaˇzdý model nám dává jiné informace. Prostorový model detailnˇe popisuje bl´ızké okol´ı sondy, ale zbytek pracovn´ı oblasti jen velmi povrchnˇe. Energetick´ y model v celé oblasti sleduje rychlé elektrony, ale informace o pomalých elektronech je velmi zkreslená. Hodnocen´ı modelu pak závis´ı souˇcasnˇe na tom, zda popisuje daný problém dostateˇcnˇe fyzikálnˇe správnˇe a zda jeho výsledky odpov´ıdaj´ı na naˇse otázky. Toto hodnocen´ı proto nem˚ uˇze být obecné. Modely srovnáváme ve velmi jednoduché podobˇe, v praxi vˇsak potˇrebujeme mnohem sloˇzitˇejˇs´ı modely. Velký význam má otázka, jaké budou nároky daného modelu ve vyˇsˇs´ıch dimenz´ıch nebo napˇr´ıklad za vyˇsˇs´ıch tlak˚ u. Mezi odpovˇedi patˇr´ı


43

napˇr´ıklad oˇcekávané pamˇet’ové nároky a moˇznosti paralelizace kódu.

6.5

V´ ysledky srovn´ an´ı

Nejprve uvedeme konkrétn´ı parametry srovnávaných model˚ u. • Spoleˇ cn´ e parametry vˇ sech model˚ u: Modely popisuj´ı interakci válcové sondy s n´ızkoteplotn´ım argonovým plazmatem. V souˇradnici jsou modely jednodimenzionáln´ı, v rychlostech jsou tˇr´ıdimenzionáln´ı (tzv. 1d3v model). Osa r je ekvidistantnˇe diskretizována na 200 bod˚ u s´ıtˇe. Délka pracovn´ı oblasti je l = 1 cm. Polomˇer Langmuirovy válcové sondy je rp = 0,1 mm. Na sondu je pˇrivedeno pˇredpˇet´ı U = 5 V. Okraj modelu je tvoˇren nenaruˇseným plazmatem s koncentrac´ı elektron˚ u a iont˚ u ne = ni = 1 · 1015 m−3 . Energetické rozdˇelen´ı elektron˚ u a iont˚ u v nenaruˇseném plazmatu povaˇzujeme za maxwellovské s teplotami Te = 23200 K a Ti = 300 K. Uvedené výsledky jsou vzorkem v ˇcase t = 1 · 10−4 s. ˇ • Iteraˇ cn´ı model: Casov´ y krok ve spojitého modelu je shodný pro elektrony −10 ˇ i ionty ∆tf = 1 · 10 s. Casov´ y krok ˇcásticového modelu elektron˚ u je ∆te = −12 5 · 10 s. V kaˇzdé iteraci vstupuje do pracovn´ı oblasti 1 · 105 elektron˚ u. ˇ • Energetick´ y model: Casov´ y krok je v ˇcásticové i spojité ˇcásti modelu stejný, ale liˇs´ı se pro elektrony (∆te = 5·10−12 s) a pro ionty (∆ti = 5·10−9 s). Poˇcet elektron˚ u v ˇcásticovém modelu silnˇe závis´ı na aktuáln´ım stavu. V ustáleném stavu je jejich poˇcet pˇribliˇznˇe 5·104, coˇz odpov´ıdá ˇcásticovému modelu s pˇribliˇznˇe 7 · 106 elektron˚ u. Maximáln´ı poˇcet elektron˚ u v pr˚ ubˇehu ustalován´ı je pˇribliˇznˇe 1,5 · 105 . • Prostorov´ y model: Model pracuje s r˚ uznými ˇcasovými kroky pro elektrony (∆te = 5 · 10−12 s) a ionty (∆ti = 5 · 10−9 s). Pro zjednoduˇsen´ı vazby mezi ˇcásticovou a spojitou ˇcást´ı modelu jsou pouˇzity tyto kroky i ve spojité ˇcásti, pˇrestoˇze by postaˇcoval delˇs´ı krok. Hranice mezi ˇcásticovým a spojitým moˇ asticový delem se nacház´ı ve vzdálenosti 0,25 cm od poˇcátku souˇradnic. C´ model je tedy omezen na 25% pracovn´ı oblasti, zbylých 75% zab´ırá spojitý model. Poˇcet ˇcástic v ˇcásticové ˇcásti je promˇenn´ y a pohybuje se okolo 6 · 105 ˇcástic. Koeficienty dif´ uze a pohyblivosti jsou povaˇzovány za konstantn´ı v celé oblasti spojitého modelu. ˇ ˇ asticov´ • C´ y model: Casov´ y krok je r˚ uzný pro elektrony (∆te = 5 · 10−12 s) a ionty (∆ti = 5 · 10−9 s). Celkový poˇcet ˇcástic je v modelu pˇribliˇznˇe 4 · 106 podle aktuáln´ıch fyzikáln´ıch podm´ınek.


44

• Spojit´ y model: Jde o základn´ı spojitý model s konstantn´ımi koeficienty difuze a pohyblivosti De = 520 m2 s−1 , Di = 0,033 m2 s−1 , µe = −170 m2 V−1 s−1 , µi = 0,45 m2 V−1 s−1 . • Spojit´ y model s LTE: Na rozd´ıl od základn´ıho spojitého modelu jsou zde koeficienty difuze a pohyblivosti promˇenné. Zjednoduˇsenˇe pˇredpokládáme, ˇze energie, kterou elektrony z´ıskávaj´ı v elektrickém poli, je d´ıky sráˇzkám lokálnˇe pˇrerozdˇelována na maxwellovské rozdˇelen´ı. Teplota elektron˚ u se mˇen´ı a podle n´ı se mˇen´ı také koeficienty difuze a pohyblivosti. Doba bˇehu jednotlivých zkoumaných hybridn´ıch model˚ u je uvedena na obrázku 6.1. Pro srovnán´ı uvád´ıme téˇz u ´ daj pro základn´ı ˇcistˇe spojitý model a ˇcásticový model. Na obrázku 6.2 jsou grafy pr˚ ubˇehu koncentrace elektron˚ u a iont˚ u. Pr˚ ubˇehy potenciálu jsou vyneseny do graf˚ u 6.3 a 6.4. ˇ asticový C´

440 min

Prostorový

82 min

Iteraˇcn´ı

35 min

Energetický

8 min

Spojitý

1 min 0

100 t [min]

Obrázek 6.1: Výsledky srovnán´ı vybraných model˚ u.

6.6

Diskuze v´ ysledk˚ u srovn´ an´ı

Porovnán´ı doby výpoˇctu splnilo naˇse oˇcekáván´ı: Výpoˇcet spojitým modelem je o nˇekolik ˇrád˚ u rychlejˇs´ı neˇz ˇcásticovým modelem jiˇz v jednodimenzionáln´ım pˇr´ıpadˇe. Pˇri pˇrechodu do v´ıce dimenz´ı se tento rozd´ıl dále prohlubuje. Hybridn´ı modely jednoznaˇcnˇe pˇrináˇs´ı velkou ˇcasovou u ´ sporu oproti ˇcásticovému modelu, nicménˇe jejich vzájemné porovnán´ı je problematické. Kaˇzdý z nich pˇrináˇs´ı proti základn´ımu spojitému modelu jiné vylepˇsen´ı — napˇr´ıklad prostorový model je zamˇeˇrený na


45

oblast sheathu, zat´ımco iteraˇcn´ı model je globáln´ı. Proto je potˇreba pˇri srovnáván´ı hybridn´ıch model˚ u brát v u ´ vahu dalˇs´ı kritéria. Srovnáme-li modely podle pr˚ ubˇeh˚ u koncentrace elektron˚ u a iont˚ u, viz obrázek 6.2, pozorujeme, ˇze nejv´ıce se výsledk˚ um ˇcásticového modelu bl´ıˇz´ı prostorový a iteraˇcn´ı model. Naopak výsledky energetického modelu nejsou uspokojivé. Energetický model totiˇz pˇredpokládá, ˇze rozdˇelovac´ı funkce elektron˚ u je bl´ızká maxwellovské a liˇs´ı se pˇredevˇs´ım ve vyˇsˇs´ıch energi´ıch. Zvolené pˇredpˇet´ı je nicménˇe dostateˇcnˇe velké, aby v okol´ı sondy byla rozdˇelovac´ı funkce podstatnˇe v´ıce naruˇsená. Prostorový model a iteraˇcn´ı model takové pˇredpoklady na rozdˇelovac´ı funkci nemaj´ı. Pr˚ ubˇehy potenciálu elektrického pole v grafech 6.3 a 6.4 rovnˇeˇz ukazuj´ı, ˇze nejlepˇs´ı shody s ˇcásticovým modelem dosáhl iteraˇcn´ı model a prostorový model. ˇ asticová Z hlediska moˇznost´ı paralelizace je nejvýhodnˇejˇs´ı iteraˇcn´ı model. C´ fáze výpoˇctu je totiˇz tvoˇrena neselfkonzistentn´ım ˇcásticovým modelem, jehoˇz paralelizace je velmi jednoduchá a efektivn´ı. Jednotlivá vlákna výpoˇctu jsou témˇeˇr nezávislá a m˚ uˇze jich v jeden okamˇzik pracovat velké mnoˇzstv´ı — v praxi jsme u ´ spˇeˇsnˇe vyzkouˇseli 8 vláken na ˇctyˇrjádrovém procesoru Intel i7. Paralelizace energetického a prostorového modelu by byla sloˇzitˇejˇs´ı a ménˇe efektivn´ı, protoˇze v tˇechto pˇr´ıpadech je ˇcásticový výpoˇcet selfkonzistentn´ı. Pˇri pˇrechodu k 2D a 3D výpoˇct˚ um má iteraˇcn´ı model také významné výhody. Stˇr´ıdán´ı ˇcásticové a spojité ˇcásti totiˇz umoˇzn ˇ uje pouˇz´ıt zcela odliˇsné prostˇredky pro jejich ˇreˇsen´ı. V tˇr´ıdimenzionáln´ım pˇr´ıpadˇe lze spojitou ˇcást ˇreˇsit profesionáln´ımi nástroji — napˇr´ıklad COMSOL Multiphysics. Naopak ˇcásticový model je vhodné naprogramovat v jazyce C++ ˇci ve Fortranu s paralelizac´ı. Pˇredáván´ı dat mezi fázemi iteraˇcn´ıho modelu je ve srovnán´ı s ostatn´ımi hybridn´ımi modely velmi jednoduché a nepˇredstavuje velkou pˇrekáˇzku pˇri pouˇzit´ı komerˇcn´ıho softwaru. Na základˇe výsledk˚ u srovnán´ı a zkuˇsenost´ı z realizace jednotlivých variant hybridn´ıch model˚ u jsme pro dalˇs´ı rozvoj vybrali iteraˇcn´ı model.


1,0

Konc. elektron˚ u Konc. iont˚ u Dif´ uzn´ı ˇreˇsen´ı

0,5

n [1015 m−3 ]

n [1015 m−3 ]

1,0

0,5

Z´ akladn´ı spojitý 0

0

0,005 r [m]

Spojitý LTE 0

0,010

0

0,5

0

0,005 r [m]

Prostorový hybridn´ı 0

0,010

0

0,005 r [m]

0,010

1,0 n [1015 m−3 ]

1,0

0,5

0,5

ˇ asticový C´

Iteraˇcn´ı hybridn´ı 0

0,010

0,5

Energetický hybridn´ı

n [1015 m−3 ]

0,005 r [m]

1,0 n [1015 m−3 ]

n [1015 m−3 ]

1,0

0

46

0

0,005 r [m]

0,010

0

0

0,005 r [m]

0,010

Obrázek 6.2: Výsledky srovnán´ı vybraných model˚ u: Závislost koncentrace elektron˚ u a iont˚ u na vzdálenosti od sondy.


47

5,0 Iteraˇcn´ı hybridn´ı Energetický hybridn´ı Prostorový hybridn´ı Z´ akladn´ı spojitý Spojitý LTE ˇ asticový C´

4,0

U [V]

3,0

2,0

1,0

0

0

0,005 r [m]

0,010

Obrázek 6.3: Výsledky srovnán´ı vybraných model˚ u: Pr˚ ubˇeh potenciálu. 5,0 Iteraˇcn´ı hybridn´ı Energetický hybridn´ı Prostorový hybridn´ı Z´ akladn´ı spojitý Spojitý LTE ˇ asticový C´

4,0

U [V]

3,0

2,0

1,0

0

0

0,001 r [m]

0,002

Obrázek 6.4: Výsledky srovnán´ı vybraných model˚ u: Detail pr˚ ubˇehu potenciálu.

Kapitola 7 Hybridn´ı model ve 2D Na základˇe výsledk˚ u kapitoly 6 a praktických zkuˇsenost´ı pˇri tvorbˇe hybridn´ıch model˚ u jsme se rozhodli pokraˇcovat ve vývoji iteraˇcn´ıho hybridn´ıho modelu. Vzhledem k ˇcasové nároˇcnosti výpoˇct˚ u ve tˇrech dimenz´ıch zaˇrazujeme jeˇstˇe kapitolu vˇenovanou dvoudimenzionáln´ımu modelu. Zde uvedené metody se totiˇz jiˇz jen málo liˇs´ı od plnˇe 3D modelu, ale výpoˇcty jsou ménˇe ˇcasovˇe nároˇcné a výsledky lze snáze zobrazit. V kapitole 8 pak uvedeme jen m´ırné u ´ pravy, které je potˇreba uˇcinit pˇri pˇrechodu k tˇr´ırozmˇerným model˚ um. Popisované metody budeme demonstrovat na modelu sondy vnoˇrené do plazmatu s driftem. Mˇejme nekoneˇcnou válcovou sondu o polomˇeru r = 0,1 mm ve ˇctvercové pracovn´ı oblasti s hranou délky L = 10 mm v pozitivn´ım sloupci doutnavého výboje v argonu pˇri tlaku 133 Pa. Podél osy x zavedeme vnˇejˇs´ı elektrické pole o velikosti E = 160 V m−1 . Koncentrace elektron˚ u a iont˚ u v nenaruˇseném 15 −3 plazmatu je ne = ni = 1,0 · 10 m . Elektrony budou podstupovat pruˇzné, excitaˇcn´ı a ionizaˇcn´ı sráˇzky s neutrály argonu se sráˇzkovými pr˚ uˇrezy podle grafu 4.2. Ionty budeme modelovat výhradnˇe spojitˇe s koeficientem difuze Di = 0,012 m2 s−1 a µi = 0,46 m2 V−1 s−1 podle [23]. Pˇredpˇet´ı na sondˇe bude U0 = 5,0 V.

7.1

Spojit´ aˇ c´ ast

Spojitou ˇcást modelu tvoˇr´ı tˇri parciáln´ı diferenciáln´ı rovnice vycházej´ıc´ı ze vztah˚ u (4.8) aˇz (4.13). ∂ne ~ − De ∇ne ) = re , + ∇ · (ne µe E (7.1) ∂t ∂ni ~ − Di ∇ni ) = ri , + ∇ · (ni µi E (7.2) ∂t e ∆U = − (ni − ne ) (7.3) ε0 48

KAPITOLA 7. HYBRIDNÍ MODEL VE 2D

49

~ = −∇U, kde ne je koncentrace elektron˚ a vztah E u a ni je koncentrace jednonásobných argonových iont˚ u. Tuto soustavu ˇreˇs´ıme metodou koneˇcných prvk˚ u programem COMSOL Multiphysics. Pro rovnice (7.1) a (7.2) zvol´ıme aplikaˇcn´ı mód ,,PDE, General Form (g)“, pro Poissonovu rovnici (7.3) je vhodný mód ,,Electrostatics (es)“. V grafickém prostˇred´ı navrhneme geometrii domény, ve které budeme provádˇet výpoˇcet, a sestroj´ıme s´ıt’ pro metodu koneˇcných prvk˚ u. Vzhledem k malému polomˇeru kˇrivosti válcové sondy je vhodné v jej´ım okol´ı s´ıt’ zjemnit, viz obrázek 4.3. Ve srovnán´ı s modelem popsaným v [23] náˇs model neobsahuje rovnici pro tok energie, protoˇze potˇrebné informace z´ıskává z detailnˇejˇs´ıho ˇcásticového modelu. Na hranic´ıch domény mus´ıme definovat hraniˇcn´ı podm´ınky. Na vnˇejˇs´ı hranici poˇzadujeme parametry nenaruˇseného plazmatu: • U = 0,8 V na hranˇe x = −5 mm, • U = −0,8 V na hranˇe x = 5 mm, ~ = 0 na ostatn´ıch hranách, • ~n · E • ne = 1,0 · 1015 m−3 , • ni = 1,0 · 1015 m−3 . Vnitˇrn´ı hranice je tvoˇrena sondou s parametry • U = 5 V, ~ · ~n − 1 vth ne , • −~n · ~je = −ne µe E 4 • ni = 0 m−3 . Pro koncentraci elektron˚ u pouˇz´ıváme Neumannovu podm´ınku podle [23], ve které ~n je normála k hranici smˇerem ven z domény, je je tok elektron˚ u podle (4.10), µe je pohyblivost elektron˚ u se znaménkem a koneˇcnˇe vth =

s

8kB Te πme

(7.4)

je tepelná rychlost elektron˚ u v bl´ızkosti sondy. Koncentraci iont˚ u v bl´ızkosti elektrody pˇredpokládáme nulovou vzhledem k jejich n´ızké energii, která je nedostateˇcná k pr˚ uniku iont˚ u k sondˇe.


7.2

50

ˇ asticov´ C´ aˇ c´ ast

ˇ asticovou ˇcást problému ˇreˇs´ıme samostatným programem v jazyce C++. Ze spoC´ jitého modelu z´ıskáváme elektrické pole, které COMSOL pro zjednoduˇsen´ı výpoˇctu interpoluje ze s´ıtˇe FEM na pravo´ uhlou rovnostrannou mˇr´ıˇz. Kaˇzdá sloˇzka elektrického pole je uloˇzena v samostatném textovém souboru. Program ˇcásticového modelu naˇcte do operaˇcn´ı pamˇeti elektrické pole a provede výpoˇcet pohybu elektron˚ u v pracovn´ı oblasti. Elektrony do oblasti vstupuj´ı ze zdroje, jehoˇz podrobnˇejˇs´ı popis uvedeme dále, a jejich pohyb je ukonˇcen bud’ dopadem na sondu, nebo opuˇstˇen´ım pracovn´ı oblasti. Pohyb elektron˚ u je urˇcen lokáln´ım elektrickým polem a sráˇzkami s ˇcásticemi neutráln´ıho pozad´ı. Pr˚ ubˇeˇznˇe sledujeme rychlost elektron˚ u, abychom na konci výpoˇctu mohli pˇredat do spojité ˇcásti modelu prostorové rozdˇelen´ı koeficient˚ u difuze a pohyblivosti, rychlostn´ı koeficienty reakc´ı a podobnˇe. K tomu obecnˇe potˇrebujeme znát rozdˇelovac´ı funkci elektron˚ u v kaˇzdé buˇ nce pravo´ uhlé mˇr´ıˇze, coˇz ovˇsem výraznˇe zvyˇsuje nároky na operaˇcn´ı pamˇet’. Proto zpracováváme data pr˚ ubˇeˇznˇe. V kaˇzdé buˇ nce poˇc´ıtáme ˇctyˇri parametry: • Koeficient difuze (4.14): De = • Sráˇzková frekvence • Koeficient excitace • Koeficient ionizace

4π 3

Z

Z

∞

νe =

∞

0

c3 f0 (c)dc. ng σela (c)

ng cσela (c)f0 (c)dc.

0

Z

∞

kexc =

Z

∞

kion =

0

0

(7.5)

(7.6)

cσexc (c)f0 (c)dc.

(7.7)

cσion (c)f0 (c)dc.

(7.8)

Funkce σela , σexc a σion jsou sráˇzkové pr˚ uˇrezy elastické, excitaˇcn´ı a ionizaˇcn´ı sráˇzky, viz graf 4.2. Sráˇzkovou frekvenci nepotˇrebujeme pˇr´ımo, naopak potˇrebujeme pohyblivost elektron˚ u µe . Pˇr´ımý výpoˇcet pohyblivosti podle (4.15) je pro náˇs model nevhodný, protoˇze obsahuje derivaci rozdˇelovac´ı funkce. Proto poˇc´ıtáme nejprve sráˇzkovou frekvenci, ze které následnˇe urˇc´ıme koeficient pohyblivosti podle (4.18). Tento postup pochopitelnˇe nen´ı pˇresný, nicménˇe velmi usnadˇ nuje výpoˇcet. Nav´ıc chyba zp˚ usobená numerickou derivac´ı rozdˇelovac´ı funkce v pˇr´ımé metodˇe podle (4.15) by byla znaˇcná.


51

Vˇsechny uvedené parametry z´ıskáváme integrac´ı pˇres rozdˇelovac´ı funkci. Napˇr´ıklad pro difuzi pˇrevedeme integrál na souˇcet 1 4π X c3i De = · , N 3 i ng σela (ci )

(7.9)

kde i bˇeˇz´ı pˇres vˇsechny vzorky dané buˇ nky a N je poˇcet vzork˚ u. Tento souˇcet m˚ uˇzeme provádˇet pr˚ ubˇeˇznˇe. Pro kaˇzdou buˇ nku potˇrebujeme ˇctyˇri ˇcásteˇcné souˇcty a poˇcet vzork˚ u v buˇ nce. Kdybychom v kaˇzdé buˇ nce ukládali rozdˇelovac´ı funkci, pamˇet’ové poˇzadavky by byly pˇribliˇznˇe desetinásobné. Výhodou tohoto neselfkonzistentn´ıho ˇcásticového modelu je snadná paralelizace. Bˇeh jednotlivých ˇcástic pracovn´ı oblast´ı je totiˇz témˇeˇr nezávislý. Pouze v okamˇziku vzorkován´ı dat do bunˇek mˇr´ıˇze a pˇri z´ıskáván´ı nové ˇcástice ze zdroje se vlákna výpoˇctu ovlivˇ nuj´ı. Výpoˇcty byly provádˇeny na ˇctyˇrjádrovém procesoru Intel i7. Vyuˇzili jsme souˇcasnˇe hyperthreading, který umoˇznil pracovat s osmi vlákny. T´ımto pˇr´ıstupem jsme výraznˇe urychlili výpoˇcet oproti p˚ uvodn´ı variantˇe popsané v [23], která byla zaloˇzená na komerˇcn´ım softwaru MATLAB. Vzhledem k tomu, ˇze v ˇcásticovém modelu pouˇz´ıváme metodu Monte Carlo, bylo nutné zvolit vhodný generátor náhodných ˇc´ısel. Ve vˇsech výpoˇctech jsme pouˇzili generátor Mersenne Twister MT19937 popsaný v publikaci [28], ve které je generátor doporuˇcen právˇe pro výpoˇcty metodou Monte Carlo. Generátor jsme u ´ spˇeˇsnˇe otestovali programem Diehard [29].

7.3

Sr´ aˇ zkov´ e procesy

Pro tuto práci jsou nejvýznamnˇejˇs´ı sráˇzky elektron˚ u s atomy neutráln´ıho argonového pozad´ı. Uvaˇzujeme sráˇzky pruˇzné, excitaˇcn´ı a ionizaˇcn´ı. Pˇr´ıˇsluˇsné sráˇzkové pr˚ uˇrezy jsou uvedeny na obrázku 4.2 podle [12]. Programové ˇreˇsen´ı sráˇzkových proces˚ u je zaloˇzné na metodˇe nulové sráˇzky popsané v kapitole 4.2. Pˇri pruˇzné sráˇzce docház´ı k pˇrerozdˇelen´ı energie elektronu a neutráln´ıho atomu. Má-li elektron pˇred sráˇzkou kinetickou energii E0 , jeho energie po sráˇzce bude E = E0

!

me 1−2 γ , mg

(7.10)

kde me je hmotnost elektronu, mg je hmotnost neutrálu a γ je náhodné ˇc´ıslo z intervalu (0, 1). Aby mohlo doj´ıt k excitaˇcn´ı nebo ionizaˇcn´ı sráˇzce, mus´ı m´ıt elektron energii vyˇsˇs´ı, neˇz je prahová energie pˇr´ısluˇsné interakce. Pˇri sráˇzce o tuto ˇcást energie pˇricház´ı. V pˇr´ıpadˇe ionizace nav´ıc ˇcást své energie pˇredává sekundárn´ımu elektronu. Prahová energie uvaˇzované excitace je Ee = 11,55 eV, pro ionizaci pak


52

´ Ei = 15,76 eV. Uhlov´ e rozdˇelen´ı elektronu po sráˇzce povaˇzujeme pro jednoduchost za izotropn´ı. Problematice sráˇzkových proces˚ u se vˇenujeme podrobnˇeji v publikaci, která tvoˇr´ı pˇr´ılohu této diplomové práce. Publikace je v tisku v ˇcasopise Vacuum. Studovali jsme v obecnˇejˇs´ı rovinˇe vliv tˇeˇzkých záporných iont˚ u a r˚ uzného sloˇzen´ı plynu neutráln´ıho pozad´ı na utváˇren´ı st´ın´ıc´ı vrstvy v okol´ı válcové sondy. Publikované výsledky byly z´ıskány selfkonzistentn´ım ˇcásticovým modelem, plánujeme vˇsak pˇrenesen´ı pouˇzitých technik do hybridn´ıho modelu.

7.4

Zdroj ˇ c´ astic

Bˇeˇznou souˇcást´ı ˇcásticových model˚ u je takzvaný zdroj ˇcástic. Obvykle totiˇz pracovn´ı oblast modelu zab´ırá jen malou ˇcást objemu plazmatu. Plazma mimo pracovn´ı oblast v takovém pˇr´ıpadˇe nahrazujeme zdrojem ˇcástic, ze kterého vcház´ı ˇ astice odcházej´ıc´ı z pracovn´ı oblasti jsou z modelu ˇcástice do pracovn´ı oblasti. C´ odstraˇ novány. Na zdroj ˇcástic jsou kladeny r˚ uzné poˇzadavky, které vedou k r˚ uzným zp˚ usob˚ um realizace zdroje ˇcástic v modelu: • Zdroj popsan´ y vzorcem Rychlostn´ı rozdˇelen´ı je popsáno vzorcem (napˇr´ıklad tok Maxwellova rozdˇelen´ı). Tato varianta je nejménˇe flexibiln´ı, protoˇze vhodné analytické vzorce existuj´ı jen pro velmi omezené mnoˇzstv´ı pˇr´ıpad˚ u. • Zdroj se zadanou rozdˇ elovac´ı funkc´ı Rychlostn´ı rozdˇelen´ı je popsáno tabulkou. Je vˇsak nutné generovat na základˇe tohoto rozdˇelen´ı tok ˇcástic. Metoda je vhodná, pokud máme dané rychlostn´ı rozdˇelen´ı (napˇr´ıklad z experimentu). • Zdroj s pomocn´ ym modelem Rychlostn´ı rozdˇelen´ı je výsledkem pomocného ˇcásticového modelu. Pˇred samotným výpoˇctem do pomocného modelu nasad´ıme ˇcástice s libovolným rychlostn´ım rozdˇelen´ım a necháme je interagovat s okol´ım, dokud nenastane rovnováha. Teprve po tomto ustálen´ı ˇcástice necháme vstupovat do modelu. Metoda je vhodná, pokud máme dané sráˇzkové procesy, vnˇejˇs´ı pole, a pod. V pr˚ ubˇehu vývoje model˚ u jsme vyzkouˇseli po ˇradˇe vˇsechny tˇri metody a ukázalo se, ˇze zdroj s pomocným modelem je pro náˇs problém nejv´ yhodnˇejˇs´ı. Jedinˇe v tomto pˇr´ıpadˇe se totiˇz podaˇrilo výraznˇe potlaˇcit vliv existence hranice modelu na výsledky výpoˇctu. To je velmi dobˇre patrné na pr˚ ubˇehu koeficientu excitace kexc argonových neutrál˚ u elektronem. V tabulce 7.4 pˇrevzaté z práce [23] jsou uvedeny

KAPITOLA 7. HYBRIDNÍ MODEL VE 2D l [m] 0,000 – 0,001 0,001 – 0,002 0,002 – 0,003 0,003 – 0,004 0,004 – 0,005 0,005 – 0,006 0,006 – 0,007 0,007 – 0,008 0,008 – 0,009 0,009 – 0,010 0,010 – 0,011 0,011 – 0,012 0,012 – 0,013 0,013 – 0,014 0,014 – 0,015 0,015 – 0,016 0,016 – 0,017 0,017 – 0,018 0,018 – 0,019 0,019 – 0,020

53

kexc [m3 s−1 ] 3,45990 · 10−16 1,28772 · 10−17 2,55010 · 10−18 4,06769 · 10−19 2,82774 · 10−19 3,94909 · 10−19 2,88464 · 10−19 2,98691 · 10−19 5,72011 · 10−19 7,68077 · 10−19 7,87768 · 10−19 8,80815 · 10−19 1,17279 · 10−18 1,22711 · 10−18 1,37210 · 10−18 1,51285 · 10−18 2,13103 · 10−18 2,81576 · 10−18 3,93644 · 10−18 1,16404 · 10−17

Tabulka 7.1: Pˇr´ıklad závislosti koeficientu excitace na vzdálenosti od válcové sondy. Data pˇrevzata z práce [23]. výsledky ˇcásticové ˇcásti hybridn´ıho modelu v ustáleném stavu. Oˇcekávali bychom, ˇze v silnˇe naruˇseném plazmatu (tj. pro malé l) se bude koeficient excitace výraznˇe mˇenit a v slabˇe naruˇseném plazmatu (pro velké l) se bude mˇenit velmi pomalu. Naˇse oˇcekáván´ı je velmi dobˇre splnˇeno v oblasti od sondy do poloviny pracovn´ı oblasti. Dále od sondy, kde je elektrické pole prakticky zanedbatelné, ovˇsem docház´ı k výraznému nár˚ ustu témˇeˇr o dva ˇrády. K tomuto jevu docház´ı, pokud zdroj nen´ı dobˇre pˇrizp˚ usobený modelu. Elektrony generované zdrojem pˇri pr˚ uchodu pracovn´ı oblast´ı podstupuj´ı nepruˇzné sráˇzky a mˇen´ı se s rostouc´ı vzdálenost´ı jejich rozdˇelovac´ı funkce. Výsledky pak silnˇe závis´ı na volbˇe velikosti pracovn´ı oblasti, coˇz je neˇzádouc´ı. Chceme-li vliv tohoto jevu podstatnˇe zmenˇsit, je nutné zajistit, aby zdroj generoval ˇcástice s takovým rozdˇelen´ım, které se nebude v pracovn´ı oblasti samovolnˇe mˇenit. Právˇe to umoˇzn ˇ uje snadno výˇse uvedená tˇret´ı varianta zdroje. Zdroj je tvoˇren malým ˇcásticovým modelem objemu plazmatu. V jeho malém objemu se nacház´ı ˇcástice, které podstupuj´ı stejné interakce jako ˇcástice v pracovn´ı oblasti


54

modelu. Pro kompenzaci nepruˇzných sráˇzek je nutné dodávat z vnˇejˇsku energii. Toho pˇrirozenˇe doc´ıl´ıme zaveden´ım konstatn´ıho podélného elektrického pole, které je typické napˇr´ıklad pro pozitivn´ı sloupec doutnavého výboje. Takto ˇreˇsený zdroj ˇcástic jsme testovali v jednoduchém 2D modelu s výˇse uvedenými parametry pro vnˇejˇs´ı elektrické pole 160 V m−1 a 800 V m−1 . Správné pˇrizp˚ usoben´ı zdroje je doloˇzeno v grafech 7.1 a 7.2, na kterých je vynesena závislost koeficientu excitace na souˇradnic´ıch v ustáleném stavu. Na ˇrezu 7.2 je zˇrejmé, ˇze rozhran´ı zdroje a modelu nemá vliv na velikost koeficientu excitace v rámci ˇsumu, který zámˇernˇe nebyl filtrován. Naopak sonda s osou v poˇcátku souˇradnic má na pr˚ ubˇeh koeficientu excitace ve svém okol´ı pˇredpokládaný vliv ve shodˇe s [23].


55

Obrázek 7.1: Koeficient excitace v 2D pracovn´ı oblasti modelu pro E = 160 V·m−1 (nahoˇre) a E = 800 V · m−1 (dole)


56

10,0 Koeficient excitace, E = 160 V · m−1 Koeficient excitace, E = 800 V · m−1

kexc [10−16 m3 s−1 ]

8,0

6,0

4,0

2,0

0

0

0,001

0,002 0,003 x [m]

0,004

0,005

kexc [10−18 m3 s−1 ]

10,0 Koeficient excitace, E = 160 V · m−1

8,0 6,0 4,0 2,0 0

0

0,001

0,002 0,003 x [m]

0,004

0,005

Obrázek 7.2: Koeficient excitace, ˇrez pro x > 0 a y = 0 z grafu 7.1. V doln´ım grafu je vynesen detail.


7.5

57

V´ ysledky modelu

V následuj´ıc´ıch grafech shrneme podstatné výsledky modelu. Potenciál elektrického pole v celé pracovn´ı oblasti je uveden na obrázku 7.3. Na obrázku 7.4 uvád´ıme potenciál v ˇrezech podle osy x (podél vnˇejˇs´ıho elektrického pole) a y (kolmo na ˇ podle osy y je symetrický, a proto uvád´ıme jen pravou vˇetev. elektrické pole). Rez

Obrázek 7.3: Potenciál elektrického pole v pracovn´ı oblasti Pr˚ ubˇehy koncentrace elektron˚ u a iont˚ u jsou uvedeny opˇet v ˇrezech na obrázku 7.5. Výsledky jsou zˇretelnˇe v souladu s odpov´ıdaj´ıc´ımi pr˚ ubˇehy na obrázku 6.2. Na dalˇs´ıch obrázc´ıch uvád´ıme výsledky ˇcásticové ˇcásti modelu — koeficient excitace na obrázku 7.6, koeficient difuze elektron˚ u na obrázku 7.7 a pohyblivost elektron˚ u na obrázku 7.8. Z graf˚ u je zˇrejmé, ˇze v okol´ı sondy se tyto parametry výraznˇe mˇen´ı. To vysvˇetluje chybné výsledky základn´ıho spojitého modelu, ve kterém je povaˇzujeme za konstanty.

58

5,0

5,0

4,0

4,0

3,0

3,0

2,0

2,0

1,0

1,0

0 −1,0 −0,005

U [V]

U [V]


0 0 x [m]

0,005

−1,0 0,005

0 y [m]

1,0

1,0

0,5

0,5

0 −0,005

0 0,005 x [m] Koncentrace elektron˚ u

0

0 0,005 y [m] Konc. iont˚ u

Obrázek 7.5: Koncentrace elektron˚ u a iont˚ u, ˇrez podle os x a y

n [1015 m−3 ]

n [1015 m−3 ]

Obrázek 7.4: Potenciál elektrického pole, ˇrez podle os x a y.

59

4,0

4,0

3,0

3,0

2,0

2,0

1,0

1,0

0 −0,005

0 x [m]

0,005

0 0,005

0 y [m]

Obrázek 7.6: Koeficient excitace, ˇrez podle os x a y

Obrázek 7.7: Koeficient difuze elektron˚ u v pracovn´ı oblasti

kexc [10−16 m3 s−1 ]

kexc [10−16 m3 s−1 ]



Obrázek 7.8: Pohyblivost elektron˚ u v pracovn´ı oblasti

60

Kapitola 8 Hybridn´ı model ve 3D a jeho aplikace Plnˇe tˇr´ırozmˇerný hybridn´ı model, jehoˇz výsledky budou tvoˇrit hlavn´ı náplˇ n této kapitoly, je zaloˇzen na dvourozmˇerném modelu, který jsme popsali v pˇredchoz´ı kapitole. Tento pˇrechod nepˇrináˇs´ı, ponˇekud pˇrekvapivˇe, ˇzádné výrazné komplikace. Pochopitelnˇe vzrostou nároky na operaˇcn´ı pamˇet’ a na ˇcas potˇrebný k výpoˇctu, a proto mus´ıme vˇenovat podstatnˇe vˇetˇs´ı pozornost konstrukci s´ıtˇe pro metodu koneˇcných prvk˚ u. Standardn´ı s´ıt’ generovaná programem COMSOL Multiphysics je obvykle naprosto nedostateˇcná. Modelován´ı interakce plazmatu s pevnou látkou totiˇz vyˇzaduje velmi jemnou s´ıt’ v oblasti st´ın´ıc´ı vrstvy. Takto jemnou s´ıt’ vˇsak nelze pouˇz´ıt v celé oblasti, protoˇze nároky na operaˇcn´ı pamˇet’ by byly pˇr´ıliˇs velké. COMSOL Multiphysics umoˇzn ˇ uje interaktivn´ı tvorbu s´ıt’ˇe. Následuj´ıc´ım jednoduchým postupem jsme z´ıskali ve vˇetˇsinˇe pˇr´ıpad˚ u dostateˇcnˇe vhodnou s´ıt’: Na povrˇs´ıch substrátu, sondy a geometricky sloˇzitˇejˇs´ıch objekt˚ u jsme vytvoˇrili jemnou ’ ’ s´ıt (oznaˇcenou Finer ). S´ıt na zbývaj´ıc´ı hranici oblasti a v objemu jsme sestrojili jako normáln´ı (oznaˇcenou Normal ), v nˇekterých pˇr´ıpadech jsme nav´ıc sn´ıˇzili hodnotu parametru Element growth rate. Vzhledem k rozsahu s´ıtˇe jsme pro ˇreˇsen´ı soustavy lineárn´ıch rovnic pouˇzili iterativn´ı ˇreˇsiˇc FGMRES (flexible generalized minimum residual ).

8.1

Sonda koneˇ cn´ ych rozmˇ er˚ u

V této podkapitole ukáˇzeme výpoˇcet st´ın´ıc´ı vrstvy v okol´ı sondy koneˇcných rozmˇer˚ u pro reálné parametry plazmatu z´ıskané ze sondového mˇeˇren´ı. Mˇeˇren´ı jsme provedli v laboratoˇri doc. Hrachové na Katedˇre fyziky povrch˚ u a plazmatu MFF UK v Praze. Parametry mˇeˇric´ı aparatury a výboje byly následuj´ıc´ı:

61

KAPITOLA 8. HYBRIDNÍ MODEL VE 3D A JEHO APLIKACE Výboj: Výbojový proud: Materiál výbojové trubice: Diagnostika: Vzdálenost sond: Polomˇer sondy: Aktivn´ı délka sondy:

62

Doutnavý v argonu, p = 150 Pa Id = 10 − 40 mA Pyrex Dvousondová d = 15 mm rp = 0,05 mm lp = 5,0 mm

Metodou dvou sond jsem urˇcili základn´ı parametry plazmatu v pozitivn´ım sloupci doutnavého výboje pro ˇctyˇri výbojové proudy. Namˇeˇrené sondové charakteristiky jsou uvedeny na obrázku 8.1, zjiˇstˇené parametry v plazmatu uvád´ıme v tabulce 8.1. 10 0 Id Id Id Id

I [10−6 A]

−10

= 10 = 20 = 30 = 40

mA mA mA mA

−20 −30 −40 −50 −60 −20

0

20

40

60 U [V]

80

100

120

140

Obrázek 8.1: Namˇeˇrené sondové charakteristiky pro r˚ uzné výbojové proudy Id Budeme modelovat interakci plazmatu s pouˇzitou sondou pˇri výbojovém proudu 20 mA. Na obrázku 8.2 je znázornˇena geometrie u ´ lohy. Pracovn´ı oblast je tvoˇrena krychl´ı o stranˇe 20 mm. Válcová sonda výˇse uvedených rozmˇer˚ u je uchycena ve válcovém dielektriku o polomˇeru 0,5 mm. Na sondu je pˇriveden potenciál 5 V. Pˇredpokládáme, ˇze sonda neemituje. Na rovinˇe x = −10 mm je definován potenciál U+ = 1,96 V, v rovinˇe x = 10 mm je U− = −1,96 V. T´ım je dáno

KAPITOLA 8. HYBRIDNÍ MODEL VE 3D A JEHO APLIKACE Výbojový proud Id [mA] Teplota elektron˚ u Te [K] Konc. elektron˚ u ne [m−3 ] Elektrické pole E [V m−1 ]

63

10 20 30 40 15 200 18 400 18 540 19 810 16 16 16 4,9 · 10 7,1 · 10 8,7 · 10 9,4 · 1016 220 196 127 152

Tabulka 8.1: Experimentálnˇe urˇcené parametry plazmatu

Obrázek 8.2: Geometrie koneˇcné válcové sondy vnˇejˇs´ı pole o velikosti E = 196 V m−1 . Na povrchu dielektrika pˇredpokládáme pro jednoduchost nulový náboj. Koncentrace elektron˚ u a iont˚ u v nenaruˇseném plazmatu je ne = ni = 7,1 · 1016 m−3 . Na sondˇe je koncentrace iont˚ u nulová. Pro koncentraci elektron˚ u na povrchu sondy pouˇzijeme Neumannovu podm´ınku shodnˇe s kapitolou 7. Výpoˇcet prob´ıhá stejným zp˚ usobem jako v modelu popisovaném v kapitole 7, pouze je rozˇs´ıˇren do plnˇe tˇr´ıdimenzionáln´ı podoby. Na obrázc´ıch 8.3, 8.4 a 8.5 uvád´ıme pr˚ ubˇehy koncentrace elektron˚ u a iont˚ u a potenciálu elektrického pole v pracovn´ı oblasti. Na obrázku 8.6 je uvedena rozdˇelovac´ı funkce elektron˚ u dopadaj´ıc´ıch na sondu zvláˇst’ pro hrot (maximálnˇe 1 mm od konce sondy) a pro válcovou ˇcást (zbytek povrchu sondy).

KAPITOLA 8. HYBRIDNÍ MODEL VE 3D A JEHO APLIKACE

Obrázek 8.3: Koncentrace elektron˚ u v okol´ı válcové sondy (ˇrezy a detail)

64


Obrázek 8.4: Koncentrace iont˚ u v okol´ı válcové sondy (ˇrezy a detail)

65


66

Obrázek 8.5: Potenciál elektrického pole v okol´ı válcové sondy (ˇrezy a detail)


0,3

f (E) [a. u.]

V´ alcov´ a ˇc´ ast Hrot 0,2

0,1

0

0

2

4

6

8 10 E [eV]

12

14

16

18

Obrázek 8.6: Rozdˇelovac´ı funkce elektron˚ u dopadaj´ıc´ıch na sondu

67


8.2

68

Substr´ at s nerovn´ ym povrchem

Plnˇe tˇr´ırozmˇerný model interakce plazmatu s pevnou látkou je vhodný pro ˇreˇsen´ı ˇrady problém˚ u se sloˇzitou povrchovou strukturou substrátu. V literatuˇre se setkáváme nejˇcastˇeji s dvourozmˇernými modely, které jsou vˇsak omezeny na zjednoduˇsené geometrické uspoˇrádán´ı — napˇr´ıklad [30] a [31]. Náˇs model takové omezen´ı nemá. Nejprve jsme provedli ˇradu výpoˇct˚ u v jednoduˇsˇs´ı geometrii. Substrát je tvoˇren nekoneˇcnou vodivou deskou v rovinˇe z = 0 s jednou nerovnost´ı ve tvaru polokoule. Stˇred podstavy polokoule je um´ıstˇen v poˇcátku souˇradnic. Pracovn´ı oblast modelu je tvoˇrena krychl´ı o stranˇe L = 10 mm, jej´ıˇz spodn´ı stˇena leˇz´ı v rovinˇe desky. Na substrát je pˇriveden potenciál 5 V. Koncentraci elektron˚ u a iont˚ u v nenaruˇseném argonovém plazmatu pˇri tlaku 133 Pa uvaˇzujeme ne = ni = 1,0 · 1015 m−3 . Sráˇzky elektron˚ u s neutrály argonu modelujeme podle popisu v kapitole 7.3. Výpoˇcet jsme provedli pro tˇri r˚ uzné polomˇery polokoule: 2,5 mm, 0,5 mm, 0,1 mm a pro substrát bez nerovnosti. Výsledky jsou shrnuty na následuj´ıc´ıch grafech. Na obrázku 8.7 uvád´ıme pr˚ ubˇeh potenciálu elektrického pole v ˇrezu rovinou y = 0. Velikost elektrického pole pro vˇsechny konfigurace shrnuj´ı obrázky 8.8 a 8.9, na kterém uvád´ıme pole podél osy z. Z graf˚ u je patrné, ˇze nerovnosti výraznˇe mˇen´ı elektrického pole v tˇesné bl´ızkosti povrchu. Na vrcholu nerovnost´ı docház´ı k podstatnému zes´ılen´ı pole.


69

Obrázek 8.7: Potenciál elektrického pole nad rovinným substrátem bez nerovnosti a s nerovnost´ı ve tvaru polokoule o polomˇeru r = 0,5 mm


10

z [mm]

z [mm]

10

5

0 −5

0 x [mm]

0 x [mm]

5

0 x [mm]

5

10

z [mm]

z [mm]

5

0 −5

5

10

5

0 −5

70

0 x [mm] 0

2

5 4

5

0 −5

6 8 10 E [103 V m−1 ]

12

14

Obrázek 8.8: Elektrické pole nad rovinným substrátem s nerovnost´ı ve tvaru polokoule o polomˇeru 2,5 mm, 0,5 mm, 0,1 mm a bez nerovnosti


71

12 Rovinn´ a sonda Polokoule r = 0,1 mm Polokoule r = 0,5 mm Polokoule r = 2,5 mm

E [103 V m−1 ]

10 8 6 4 2 0

0

0,005 z [m]

0,010

Obrázek 8.9: Rovinný substrát s nerovnost´ı ve tvaru polokoule o r˚ uzném polomˇeru r.


72

V dalˇs´ıch výpoˇctech jsme studovali substrát s v´ıce nerovnostmi v pravidelném uspoˇrádán´ı. Vzdálené nerovnosti se vzájemnˇe neovlivˇ nuj´ı a lze je studovat oddˇelenˇe. Pˇribl´ıˇz´ıme-li nerovnosti tak, ˇze jejich vzdálenost je srovnatelná s jejich velikost´ı, vykazuj´ı sloˇzitˇejˇs´ı chován´ı. Na následuj´ıc´ıch obrázc´ıch prezentujeme výsledky dvou konfigurac´ı. V obou pˇr´ıpadech jsme na substrát um´ıstili ˇsestnáct polokoul´ı o polomˇeru r = 0,5 mm. Uspoˇrádán´ı je ˇctvercové. Konfigurace se liˇs´ı vzdálenost´ı stˇred˚ u podstav polokoul´ı. V prvn´ım pˇr´ıpadˇe, který oznaˇc´ıme jako tˇesné uspoˇrádán´ı, je tato vzdálenost r1 = 1,0 mm, v druhém pˇr´ıpadˇe, který nazveme volné uspoˇrádán´ı, pak je r2 = 1,25 mm. Na obrázku 8.10 uvád´ıme jako pˇr´ıklad vizualizaci s´ıtˇe na substrátu s volným uspoˇrádán´ım polokoul´ı. Velikost elektrického pole nad substrátem s tˇesným uspoˇrádán´ım shrnuj´ı obrázky 8.11 a 8.12, pro volné uspoˇrádán´ı uvád´ıme odpov´ıdaj´ıc´ı výsledky na obrázc´ıch 8.13 a 8.14. V pˇr´ıpadˇe tˇesného uspoˇrádán´ı docház´ı k lokalizaci silného elektrického pole nad vnitˇrn´ımi polokoulemi a v mezerách elektrické pole výraznˇe klesá aˇz k nule. Nad okrajovými polokoulemi se pr˚ ubˇeh elektrického pole bl´ıˇz´ı k výˇse uvedenému pˇr´ıkladu se samostatnou polokoul´ı. Výsledky konfigurace s volným uspoˇrádán´ım polokoul´ı ukazuj´ı, ˇze jiˇz pˇri relativnˇe malém oddálen´ı polokoul´ı docház´ı k výraznému pˇribl´ıˇzen´ı k pˇr´ıpadu samostatné polokoule. Rozd´ıly mezi vnitˇrn´ı a okrajovou polokoul´ı se sniˇzuj´ı, coˇz je patrné ze srovnán´ı kˇrivek A (nad okrajovou polokoul´ı) a B (nad vnitˇrn´ı polokoul´ı) na obrázku 8.14 s tˇesným uspoˇrádán´ım na obrázku 8.12. Výpoˇcet jedné konfigurace trval obvykle nˇekolik hodin. Pˇresnou dobu lze tˇeˇzko stanovit, protoˇze závisela na mnoha faktorech. Jedna iterace spojitého ˇcásti trˇ asticový movala nejvýˇse nˇekolik minut, coˇz je v celkovém ˇcase zanedbatelné. C´ del dokázal d´ıky paralelizaci na osm vláken na ˇctyˇrjádrovém procesoru Intel i7 s hyperthreadingem zpracovat 106 ˇcástic za 15 minut. V rámci jedné iterace jsme poˇzadovali 3 · 106 aˇz 10 · 106 ˇcástic. Po tˇret´ı iteraci jiˇz nedocházelo k podstatn´ ym zmˇenám výsledk˚ u. Dalˇs´ı iterace slouˇzily pˇredevˇs´ım k z´ıskán´ı detailnˇejˇs´ıch výsledk˚ u. Do celkové doby výpoˇctu je potˇreba zahrnout také pˇr´ıpravu geometrie, konstrukci vhodné s´ıtˇe metody koneˇcných prvk˚ u a pˇr´ıpadné u ´ pravy v kódu ˇcásticového modelu. Tyto ˇcinnosti v nˇekterých pˇr´ıpadech trvaly delˇs´ı dobu neˇz samotný výpoˇcet. Z´ıskané výsledky budou prezentovány na mezinárodn´ı konferenci CIP 09 v Marseille a poslány do tisku do ˇcasopisu European Physical Journal - Applied Physic. Abstrakt je pˇr´ılohou diplomové práce.


73

Obrázek 8.10: Vizualizace s´ıtˇe metody koneˇcných prvk˚ u na nerovném substrátu.


74

Obrázek 8.11: Velikost elektrického pole nad substrátem s nerovnostmi tˇesného uspoˇrádán´ı (nahoˇre: pole nad nerovnostmi, dole: pole nad mezerou)


75

y [mm]

5

0

−5 −5

0 x [mm]

5

E [V m−1 ]

10000 A B C D 5000

0

0

0,005 z [m]

0,010

Obrázek 8.12: Pr˚ ubˇeh velikosti elektrického pole ve smˇeru osy z nad substrátem s nerovnostmi tˇesného uspoˇrádán´ı (nahoˇre: geometrie substrátu s vyznaˇcenou polohou vzork˚ u)


76

Obrázek 8.13: Velikost elektrického pole nad substrátem s nerovnostmi volného uspoˇrádán´ı (nahoˇre: pole nad nerovnostmi, dole: pole nad mezerou)


77

y [mm]

5

0

−5 −5

0 x [mm]

5

E [V m−1 ]

10000 A B C D 5000

0

0

0,005 z [m]

0,010

Obrázek 8.14: Pr˚ ubˇeh velikosti elektrického pole ve smˇeru osy z nad substrátem s nerovnostmi volného uspoˇrádán´ı (nahoˇre: geometrie substrátu s vyznaˇcenou polohou vzork˚ u)

Kapitola 9 Z´ avˇ er V této práci jsme se vˇenovali hybridn´ım metodám modelován´ı interakce plazmatu a pevné látky, které nám umoˇzn ˇ uj´ı studovat geometricky sloˇzité konfigurace vyˇzaduj´ıc´ı plnˇe tˇr´ıdimenzionáln´ı pˇr´ıstup pˇri pomˇernˇe krátké dobˇe výpoˇctu. V kapitole 6 jsme provedli srovnán´ı r˚ uzných metod hybridn´ıho modelován´ı se základn´ımi metodami. Vˇsemi metodami byla ˇreˇsena stejná u ´ loha — nekoneˇcná válcová sonda v n´ızkoteplotn´ım argonovém plazmatu — a byly porovnány výsledky modelu a doba výpoˇctu. Zjiˇstˇené doby výpoˇct˚ u hybridn´ımi modely byly podle naˇseho oˇcekáván´ı výraznˇe kratˇs´ı neˇz výpoˇcet ˇcásticovým modelem a delˇs´ı neˇz spojitým modelem. Provedli jsme porovnán´ı pr˚ ubˇeh˚ u potenciálu a koncentrac´ı elektron˚ u a iont˚ u v pracovn´ı oblasti a doˇsli jsme k závˇeru, ˇze se výsledky hybridn´ıch model˚ u výraznˇe bl´ıˇz´ı k výsledk˚ um ˇcásticového modelu, který povaˇzujeme za referenˇcn´ı. Na základˇe výsledk˚ u srovnán´ı v kapitole 6 a praktick´ ych zkuˇsenost´ı pˇri vývoji model˚ u jsme se rozhodli dále rozv´ıjet iteraˇcn´ı hybridn´ı model. Navázali jsme na práci [23] a vytvoˇrili jsme plnˇe tˇr´ıdimenzionáln´ı iteraˇcn´ı hybridn´ı model. Náˇs model ve srovnán´ı s [23] pˇrinesl ˇradu inovac´ı. Spojitá ˇcást je ˇreˇsena metodou koneˇcných prvk˚ u v programu COMSOL Multiphysics. Koeficienty difuze a pohyblivosti vystupuj´ıc´ı ve spojité ˇcásti modelu jsou z´ıskávány z ˇcásticového neselfkonzistentn´ıho modelu ve velmi husté mˇr´ıˇzi, d´ıky které dostáváme velmi detailn´ı výsledky pˇredevˇs´ım v bl´ızkosti substrátu. Z´ıskán´ı tˇechto výsledk˚ u je vˇsak podm´ınˇeno vysokým poˇctem ˇcástic v ˇcásticovém modelu. Toho jsme dosáhli kompletn´ım pˇrepracován´ım ˇcásticového modelu, který je novˇe tvoˇren samostatným programem v jazyce C++ s velmi efektivn´ı paralelizac´ı, která umoˇzn ˇ uje plné vyuˇzit´ı dostupných poˇc´ıtaˇcových prostˇredk˚ u. V kapitole 7 jsme vyvinutý model prezentovali na jednoduchém dvoudimenzionáln´ım pˇr´ıkladu. Podrobnˇeji jsme se také vˇenovali problematice zdroje ˇcástic v ˇcásticovém modelu. Ukázali jsme, ˇze pˇrizp˚ usoben´ı zdroje podm´ınkám v pracovn´ı oblasti modelu vede k výraznému zlepˇsen´ı výsledk˚ u na okraji oblasti. 78

´ ER ˇ KAPITOLA 9. ZAV

79

Aplikace plnˇe tˇr´ıdimenzionáln´ıho hybridn´ıho modelu jsme uvedli v kapitole 8. Na základˇe namˇeˇrených parametr˚ u plazmatu jsme vytvoˇrili model koneˇcné válcové sondy s nevodivým uchycen´ım k aparatuˇre v n´ızkoteplotn´ım plazmatu s driftem. Výsledky výpoˇctu jsou prezentovány na ˇradˇe graf˚ u. Dále jsme model pouˇzili k výpoˇctu st´ın´ıc´ı vrstvy v okol´ı rovinného substrátu s nerovnostmi ve tvaru polokoule. Studovali jsme vliv polomˇeru kˇrivosti polokoul´ı a prostorového uspoˇrádán´ı pˇredevˇs´ım na velikost elektrického pole v b´ızkosti substrátu. Vyvinutý model budeme dále rozˇsiˇrovat. Oˇcekáváme, ˇze model bude vhodný pro oblast vyˇsˇs´ıch tlak˚ u. Rozvoj t´ımto smˇerem je limitovaný pˇredevˇs´ım ˇcásticovou ˇcást´ı modelu. Dále budeme roˇsiˇrovat model do oblasti multikomponentn´ıho plazmatu. Pˇr´ılohou diplomové práce je publikace z této oblasti. Výsledky diplomové práce budeme prezentovat na mezinárodn´ıch konferenc´ıch a budeme usilovat o jejich publikaci v odborných ˇcasopisech.

Literatura [1] Langmuir, I.: Oscillation in ionized gases, Proc. Nat. Acad. Sci. U.S., 14 628, 1928. ´ [2] Chen, F. F.: Uvod do fyziky plazmatu, Academia, Praha, 1984. [3] Martiˇsovitˇs, V.: Základy fyziky plazmy, Fakulta matematiky, fyziky a informatiky UK Bratislava, Bratislava, 2004. [4] Langmuir, I., Mott-Smith, H. M.: The Theory of Collectors in Gaseous Discharges, Phys. Rev., 28 727, 1926. [5] Swift, J. D., Schwar, M. J. R.: Electrical Probes for Plasma Diagnostics, Iliffe Books Ltd., London, 1970. ´ ı nad Labem, 2003. [6] Hrach, R.: Poˇc´ıtaˇcová fyzika, PF UJEP, Ust´ [7] Hockney, R. W., Eastwood, J. W.: Computer Simulation Using Particles, Taylor and Francis, New York, 1988. ˇ [8] Simek, J.: Dizertaˇcn´ı práce: Rozvoj metod poˇc´ıtaˇcové fyziky pro fyziku plazmatu a fyziku tenkých vrstev, MFF UK v Praze, Praha, 2006. [9] Birdsall, C. K., Langdon, A. B.: Plasma physics via computer simulation, Adam Hilger, Bristol, 1991. [10] Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P.: Numerical Recipes, The Art of Scientific Computing, Third Edition, Cambridge University Press, New York, 2007. [11] Davis, T. A.: Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method, ACM Transactions on Mathematical Software, 30 196–199, 2004. [12] Bogaerts, A., Gijbels, R.: Monte Carlo model for the argon ions and fast argon atoms in a radio-frequency discharge, IEEE Trans. Plasma Science, 27 1406, 1999. 80

LITERATURA

81

[13] Skullerud, H. R.: The stochastic computer simulation of ion motion in a gas subjected to a constant electric field, J. Phys. D: Appl. Phys., 1 1567, 1968. ˇ Diplomová práce: Studium interakce plazma-pevn´ [14] Rouˇcka, S.: a l´ atka pˇri stˇredn´ıch tlac´ıch, MFF UK v Praze, Praha, 2008. [15] Trunec, D. et al.: Monte Carlo Simulations of the Electron Currents Collected by Electrostatic Probes, Contrib. Plasma Phys., 44 577–581, 2004. [16] Huxley, L. G. H., Crompton, R. W.: The Diffusion and Drift of Electrons in Gases, v ruském jazyce, Mir, Moskva, 1977. [17] Robertson, S., Sternovsky, Z.: Monte Carlo model of ion mobility and diffusion for low and high electric fields, Phys. Rev. E , 67 046405, 2003. [18] Fletcher, C. A. J.: Computational Techniques for Fluid Dynamics, Volume 1, Second Edition, Springer-Verlag Berlin Heidelberg, 1991. [19] Jel´ınek, P.: Dizertaˇcn´ı práce: Pokroˇcilé techniky poˇc´ıtaˇcového modelov´ an´ı ve fyzice plazmatu, MFF UK v Praze, Praha, 2007. [20] Passchier, J. D. P., Goedheer, W. J.: A two-dimensional fluid model for an argon rf discharge, J. Appl. Phys., 74 3744, 1993. [21] Wendt, J. F.: Computational Fluid Dynamics, Springer-Verlag Berlin Heidelberg, 1992. [22] Kim, H. C. et al.: Particle and fluid simulations of low-temperature plasma discharges: benchmarks and kinetic effects, J. Phys. D: Appl. Phys., 38 R283– R301, 2005. [23] Bartoˇs, P.: Dizertaˇcn´ı práce: Hybridn´ı modelov´ an´ı ve fyzice plazmatu, MFF UK v Praze, Praha, 2007. [24] Bogaerts, A., Gijbels, R., Goedheer, W.: Hybrid Modeling of a Capacitively Coupled Radio Frequency Glow Discharge in Argon: Combined Monte Carlo and Fluid Model, Jpn. J. Appl. Phys., 38 4404–4415, 1999. [25] Bogaerts, A., Gijbels, R.: Effect of small amounts of hydrogen added to argon glow discharges: Hybrid Monte Carlo-fluid model, Phys. Rev. E., 65 056402, 2002. [26] Boost C++ Libraries [online], posledn´ı revize 23. 4. 2008 [cit. 7. 11. 2008]. Dostupné z: .

LITERATURA

82

[27] Kroesen, G. M. W. et al.: The Energy Balance of a Plasma in Partial Local Thermodynamic Equilibrium, IEEE Trans. Plasma Science, 18 985, 1990. [28] Matsumoto, M., Nishimura, T.: Mersenne Twister: A 623-Dimensionally Equidistributed Uniform Pseudo-Random Number Generator, ACM Trans. on Modeling and Computer Simulation, 8 3–30, 1998. [29] The Marsaglia Random Number CDROM including the Diehard Battery of Tests of Randomness [online], posledn´ı revize v r. 1995 [cit. 27. 3. 2009]. Dostupné z: . [30] Demokan, O., Filiz, Y.: Ion-matrix sheaths related to planar targets with semicylindrical grooves, J. Appl. Phys., 93 83, 2003. [31] Miyagawa, Y. et al.: Plasma analysis for the plasma immersion ion implantation processing by a PIC-MCC simulation, Comp. Phys. Comm., 177 84–87, 2007.

Pˇ r´ılohy K diplomové práci pˇrikládáme výsledky naˇs´ı publikaˇcn´ı ˇcinnosti: ˇ • Hrubý, V., Hrach, R., Simek, J.: Computational study of plasma composition influence on sheath formation Proc. JVC-12/EVC-10 Balatonalm´ adi, Mad’arsko, 22. aˇz 26. 9. 2008 Vacuum (v tisku) • Hrach, R., Hrubý, V.: Computational analysis of plasma interaction with uneven substrates Proc. CIP 2009 Marseille, Francie, 22. aˇz 26. 6. 2009

83

Univerzita Karlova v Praze Matematicko-fyzikální fakulta

Recommend Documents