UNIVERSITA’ DEGLI STUDI DI BARI - uniba...

UNIVERSITA’ DEGLI STUDI DI BARI

FACOLTA’ DI SCIENZE MATEMATICHE FISICHE E NATURALI

CORSO DI LAUREA IN FISICA ___________________________________________________

TESI DI LAUREA

DINAMICA E CINETICA DI UN PLASMA RF

TRA PIATTI PARALLELI IN OSSIGENO

Relatori: Laureando:

Prof. Mario Capitelli Daniele Marinazzo

Prof. Savino Longo

___________________________________________________

ANNO ACCADEMICO 2000 - 01

INDICE INTRODUZIONE ..................................................................................1 CAPITOLO 1 CARATTERISTICHE DI UN PLASMA RF.........................................5

1.1 Il plasma RF ..............................................................................................5 1.2 Lunghezza di Debye e moti collettivi..............................................................6 1.3 Sheath e bulk .............................................................................................9 1.4 Modello circuitale dello sheath......................................................................10 1.5 Diffusione ambipolare ...................................................................................12 1.6 Riscaldamento elettronico .............................................................................14

CAPITOLO 2 CARATTERISTICHE MATEMATICHE E FISICHE DEL MODELLO...........................................................................................17

2.1 Modello numerico per la simulazione: PIC/MCC.........................................17 2.1.1 PIC..........................................................................................................17 2.1.2 MCC .......................................................................................................20

2.2 Modifica del ciclo temporale.........................................................................24 2.3 Accoppiamento del modello con la cinetica..................................................26 2.4 Dinamica delle collisioni ...............................................................................29

CAPITOLO 3 REAZIONI ED ENERGIA...................................................................33

3.1 Specie neutre presenti nel modello................................................................33 3.2 Sezioni d’urto e costanti di velocità ..............................................................34 3.3 Reazioni da impatto elettronico.....................................................................36

3.3.1 Eccitazioni verso stati metastabili ...................................................36 3.3.2 Diseccitazioni ..................................................................................37 3.3.3. Dissociative attachments O2 ............................................................40 3.3.4 Ionizzazioni .....................................................................................41 3.3.5 Eccitazione dissociativa...................................................................42 3.3.6 Perdite di 4.5, 6.0, 8.4 e 9.97 eV .....................................................44 3.3.7 Momentum transfer .........................................................................45 3.3.8 Dissociative attachment O3..............................................................46 3.3.9 Ricombinazioni...............................................................................46 3.3.10 Collisioni tra specie pesanti.............................................................47 3.3.11 Scambio di carica............................................................................48 3.3.12 Momentum transfer ........................................................................49 3.3.13 Emissione secondaria ...........................................................................50

3.4 Verso l’equazione del calore .........................................................................50 3.5 L’equazione del calore...................................................................................53

3.5.1 Il termine sorgente ..................................................................................53 3.5.2 Conduttività ............................................................................................56

CAPITOLO 4 DISCRETIZZAZIONE ED IMPLEMENTAZIONE DELL’EQUAZIONE DEL CALORE..................................................57

4.1 Forma generale dell’equazione da discretizzare............................................57 4.2 Regole di base................................................................................................59 4.3 Nonlinearità ...................................................................................................65 4.4 Condizioni al contorno ..................................................................................67 4.5 Risoluzione delle equazioni algebriche .........................................................71 4.6 Sovrarilassamento e sottorilassamento ..........................................................75 4.7 Forma finale dell’ equazione discretizzata ....................................................77

CAPITOLO 5 RISULTATI..........................................................................................79

5.1 Temperatura...................................................................................................79 5.2 Potenziale.......................................................................................................80 5.3 Densità ...........................................................................................................82 5.4 Energia media ................................................................................................86 5.5 Densità, energia e velocità in un ciclo RF .....................................................93 5.6 Funzione di distribuzione dell’energia degli elettroni (EEDF) .....................95 5.7 Funzione di distribuzione di energia degli ioni (IEDF).................................97 5.8 Potenza depositata e reazioni.......................................................................106 5.9 Ionizzazione.................................................................................................114

CONCLUSIONI .................................................................................115 APPENDICE LIVELLI ENERGETICI DELLA MOLECOLA DI OSSIGENO.....117 BIBLIOGRAFIA ................................................................................119

1

INTRODUZIONE

I plasmi da scarica chimicamente reattivi sono ampiamente usati

per modificare le proprietà superficiali dei materiali . Processi

superficiali basati sul plasma sono indispensabili nel costruire circuiti

integrati su ampia scala usati dall’ industria elettronica . Tali processi

sono anche di vitale importanza per l’industria aerospaziale,

automobilistica, biomedica, per le acciaierie e per lo smaltimento di

rifiuti tossici .

Col plasma si possono fabbricare materiali e strutture superficiali

altrimenti non ottenibili e le proprietà superficiali dei materiali possono

essere modificati in maniera unica. Ad esempio si possono incidere

solchi larghi 0.2 µm e profondi 4 µm in pellicole o substrati di silicio

per isolare un dispositivo o costruire un condensatore per

l’immagazzinamento di carica in un circuito integrato (il diametro di un

capello è tra i 50 e i 100 µm). Tramite il plasma sono anche stati

prodotti materiali quali i film di diamante e il silicio amorfo .

Le scariche in ossigeno sono usate per depositare tramite

sputtering alluminio, tungsteno o pellicole superconduttrici ad alta

temperatura; per crescere film di SiO2 sul silicio; per rimuovere

selettivamente film di silicio; per rimuovere pellicole fotoresistenti o

pellicole di polimeri .

Per la microcostruzione di un circuito integrato, circa un terzo

delle decine di migliaia di passi di fabbricazione sono basati

sull’utilizzo di plasmi .

Inoltre, una caratteristica delle scariche a bassa pressione che

sono trattate in questo lavoro è che il plasma, e così anche il sistema

plasma-substrato , non è in equilibrio termico, come sarà mostrato in

2

seguito. Ciò permette temperature del substrato relativamente basse,

pur mantenendo tassi adeguati di deposizione o di etching anche

laddove si devono trattare materiali termolabili quali i polimeri

organici. I plasmi di non equilibrio sono caratterizzati da una elevata

energia elettronica e da un’alta concentrazione di particelle cariche ed

eccitate.

L’ottimizzazione di tutti i processi che riguardano i plasmi di non

equilibrio richiede la soluzione di due problemi intrecciati: la

determinazione dei parametri interni (concentrazione di particelle

cariche ed eccitate e di prodotti di reazioni chimiche, e loro

distribuzioni di energia) dati i parametri esterni della scarica

(composizione del gas, pressione; frequenza, potenza ed intensità della

corrente si scarica, geometria del reattore), e il calcolo della cinetica

delle reazioni chimiche nel volume del plasma, sugli elettrodi e sulle

pareti.

Alla luce di ciò si evince l’utilità di un lavoro teorico che vada

fianco a fianco con quello sperimentale. Una simulazione infatti

permette il controllo rapido ed immediato, oltre che poco dispendioso,

di tutti i parametri di cui si è parlato in precedenza, ed è anche un

utilissimo mezzo per comprendere il comportamento delle particelle

all’ interno del plasma, consentendo anche un conteggio preciso di ogni

tipo di collisione.

Questo lavoro è rivolto alla simulazione di una scarica a

radiofrequenza in ossigeno puro, che è stata solo raramente considerata

come caso test per modelli fluidi e cinetici di plasmi di non equilibrio,

lasciando relativamente aperto il campo di ricerca.

Il metodo numerico per la simulazione prende il nome di

PIC/MCC (Particle In Cell with Monte Carlo Collisions). Nel modello

qui in esame è stata prestata particolare attenzione alla cinetica degli

3

elettroni e degli ioni. All’interno di questo modello si sono studiate

dettagliatamente tutte le reazioni che hanno avute per protagoniste le

specie chimiche da noi simulate, con particolare attenzione all’energia

assorbita o rilasciata nel plasma dopo ognuna di esse.

Parallelamente sono stati monitorati tutti i processi di semplice

collisione e le conseguenti variazioni di energia cinetica. Queste

energie sono state raccolte ed utilizzate come sorgente in una

equazione allo stato stazionario per la diffusione del calore, dando

origine ad un profilo di temperatura.

Questo tipo di studio non era finora mai stato condotto, e si vede

come il profilo di temperatura modifica l’energia e la densità delle

specie presenti nel nostro modello.

Nel primo capitolo della presente tesi si farà riferimento alle

caratteristiche generali di un plasma di non equilibrio.

Il secondo capitolo illustra il modello utilizzato per la

simulazione e descrive l’accoppiamento tra cinetica e modello

PIC/MCC.

Il terzo capitolo è una esposizione delle reazioni considerate,

delle rispettive sezioni d’urto e del loro apporto energetico. Inoltre

saranno gettate le basi per la costruzione dell’equazione di diffusione

del calore.

Il quarto capitolo descrive in dettaglio l’implementazione della

equazione della diffusione del calore.

Infine il quinto capitolo è dedicato alla discussione dei risultati.

Il codice PIC/MCC usato per la simulazione fa parte dell’ultima

generazione di modelli che sono stati sviluppati a partire dai primi anni

Novanta presso il Dipartimento di Chimica dell’ Università di Bari e il

Centro di Studi per la Chimica dei Plasmi del CNR. Si vedano in

4

proposito i lavori di M.Capitelli, S.Longo [1,2,3] e la tesi di Laurea in

Fisica di D.Iasillo[4].

Il lavoro sulla cinetica è stato svolto in collaborazione con

K.Hassouni e W.Morscheidt dell’ Università di Paris 13.

5

CAPITOLO 1 CARATTERISTICHE DI UN PLASMA RF

1.1 Il plasma RF

Un plasma generato tramite scarica a radiofrequenza è solo uno

dei numerosissimi plasmi esistenti in natura o che possono essere

prodotti in laboratorio, che spaziano su un range enorme di densità e di

temperatura. Esso viene prodotto in un reattore costituito da un cilindro

le cui basi sono gli elettrodi.

Questo tipo di plasma è caratterizzato da una temperatura

elettronica dell’ordine di qualche eV e da una densità di particelle

cariche n* tra 108 e 1013 cm-3. Per questo motivo un plasma di non

equilibrio è detto anche plasma freddo o a bassa pressione e la scarica

che lo produce è detta glow discharge, dal bagliore diffuso dovuto alla

radiazione emessa dalle specie che si diseccitano.

Un sistema siffatto è caratterizzato dalla molteplicità di

canali che variano in maniera diversa in diversi intervalli di energia e si

influenzano l’un l’altro. Questi processi comportano non equilibrio,

non linearità, non stazionarietà ed instabilità.

La densità di energia nel plasma è tale che in esso si

verificano continuamente processi chimici ed è quindi necessario, data

la dipendenza dall’energia di questi processi, tenere sempre controllo la

variazione di energia e tener conto dei contributi dovuti agli stati

eccitati delle particelle e dei prodotti di reazione.

Lo stato stazionario di nonequilibrio è sostenuto dal flusso

attraverso il sistema dell’energia presa dal campo elettrico esterno ed

* Si fa l’ipotesi di quasineutralità per cui n rappresenta due quantità circa uguali, la densità elettronica e quella ionica

6

usato in primo luogo per scaldare gli elettroni. Come risultato della

collisione con particelle pesanti, l’energia elettronica è usata per

eccitare i vari livelli, per la dissociazione delle molecole e la

ionizzazione delle molecole e degli atomi nei processi inelastici, ed è

trasferita ai gradi di libertà traslazionali nelle collisioni elastiche. Il

rilassamento delle energie rotazionali e vibrazionali porta al

trasferimento di questa energia ai gradi di libertà traslazionali delle

particelle pesanti e al riscaldamento del gas.

1.2 Lunghezza di Debye e moti collettivi

Le forze che agiscono su una particella carica possono essere

divise in due tipi: quelle dovute a cariche distanti, che danno un campo

lentamente variabile in prossimità della particella e contribuiscono così

al moto collettivo delle particelle; e quelle dovute a particelle vicine

che la fanno muovere in maniera non correlata rispetto ad esse. I moti

correlati sono governati da un tempo detto periodo di plasma, tp. ��

��

+ - + - - + - + + - + - - + - + + - + - - + - +

++++++

------

��

+ - + - + - - + - + - + + - + - + - - + - + - + + - + - + - - + - + - +

��

+ - + - - + - + + - + - - + - + + - + - - + - +

++++++

------

Fig. 1.1 Oscillazioni di plasma

Consideriamo una «fetta» di plasma, soggetta ad un campo

elettrico che separa leggermente le cariche e crea due strati sottili di

polarità opposta, ognuno di spessore x come indicato in figura 1.1. Un

7

foglio infinito uniformemente carico con densità di carica origina un

campo uniforme σ/2ε0 nelle sue vicinanze. Dato che la densità di carica

superficiale in ogni strato ha modulo nex, dove e è la carica elementare,

e le due distribuzioni si rinforzano a vicenda, il campo totale è nex/ε0.

Ogni elettrone, di massa m, è soggetto ad una forza –ne2x/ε0 e quindi

compie un moto armonico con frequenza angolare, detta frequenza di

plasma,

ωP = ne2

mε0

12

(1.1)

corrispondente ad un periodo di plasma

tP = 2πne2

mε0

−12

(1.2) .

La distanza entro la quale le interazioni non correlate tra

particelle sono considerate importanti è detta lunghezza di Debye, che

indicheremo con D. Se consideriamo una sfera di Debye, una sfera di

raggio D che circonda la particella, l’effetto della distribuzione di

carica oltre una distanza D è sostanzialmente ridotto dallo schermo

costituito dalle particelle circostanti.

Consideriamo infatti una carica di prova q circondata da

elettroni in equilibrio termico. Risolvendo l’equazione di Poisson si

può determinare il potenziale che agisce sulla particella [5] :

8

φ = q4πε0r

exp − rλD

(1.3)

con

nekT

D 20ελ = (1.4)

è la lunghezza di Debye. E’ evidente che per r<<D il campo è

dovuto alla carica di prova, ma per r dell’ordine di D o maggiori il

campo è sostanzialmente ridotto (schermato) dagli elettroni.

Per distanze minori della lunghezza di Debye i campi

dovuti a particelle cariche individuali sono predominanti è il plasma si

comporta come un insieme di particelle non correlate. Su grandi

distanze il comportamento individuale è schermato dall’azione

collettiva di di molte particelle. Infatti, per tali distanze, si può avere

solo un comportamento collettivo, in cui le particelle si comportano in

maniera coerente, sotto forma di un’onda di plasma. La simulazione

tramite particle-in-cell ha senso in quest’ultima regione.

Infatti, il comportamento delle particelle nel regime a corto

raggio (individuale) è caratterizzato da effetti collisionali che si

differenziano da quelli che si incontrano nei gas per l’importanza

dell’interazione simultanea della particella di prova con tanti altri corpi.

Il problema è adesso quello di decidere quanto sono importanti i moti

indotti dalle interazioni locali rispetto ai moti collettivi dovuti alle

particelle più lontane. Un indicatore per questo problema è il rapporto

tra tc e tP dove tc è il tempo medio tra due collisioni in cui una

particella è deviata di un angolo maggiore di π/2 e tP è il periodo di

9

plasma. Si può dimostrare [6] che questo rapporto è dell’ ordine del

numero di particelle nella sfera di Debye, ND.

Per questo motivo, quando ND è >>1 si può portare avanti

la simulazione per un gran numero di periodi di plasma senza che siano

importanti gli effetti dovuti alle collisioni. In queste condizioni il

sistema può essere considerato non collisionale e il modello può essere

costruito sulla base di valori medi del campo in prossimità di ogni

particella.

Ovviamente il nostro è il modello di un plasma in queste

condizioni, detto plasma ideale. 1.3 Sheath e bulk

I plasmi, che sono quasineutri, sono collegati alle pareti del

reattore da strati sottili carichi positivamente, detti sheaths . Per capire

come mai, si noti in primo luogo che la velocità termica degli elettroni

(eTe/m)1/2 è almeno 100 volte maggiore della velocità termica ionica

(eTi/M)1/2 in quanto m/M<<1 e Te<<Ti.

Consideriamo un plasma di larghezza l con ne=ni inizialmente

confinato tra due pareti assorbenti a potenziale zero (fig.1.2). Dato che

la densità di carica netta =e(ni-ne) è zero, il potenziale elettrico e il

campo elettrico Ex sono nulli ovunque. Quindi gli elettroni veloci non

sono confinati e si perderanno rapidamente contro le pareti. Su una

scala di tempi molto piccola alcuni elettroni vengono persi in

prossimità delle pareti, portando alla situazione mostrata in figura 2.1.

Si formano delle guaine (sheaths) sottili (s<<l) vicino a ciascuna parete

in cui ni>>ne. La densità di carica netta positiva all’interno degli

sheaths porta ad un profilo del potenziale (x) che è positivo all’interno

10

del plasma e diventa rapidamente zero avvicinandosi alle pareti. Questo

profilo funge da «valle» di potenziale per gli elettroni e da «collina»

per gli ioni, poichè i campi elettrici all’interno degli sheaths vanno dal

plasma alle pareti. Quindi la forza –eEx che agisce sugli elettroni è

diretta verso l’interno del plasma; questo porta ad una riflessione

all’indietro degli elettroni che si muovono verso le pareti. Al contrario,

gli ioni che dal plasma entrano negli sheaths sono accelerati verso le

pareti. Se il potenziale del plasma (rispetto alle pareti) è VP , ci si

aspetta che VP sia pari a qualche eV per confinare la maggior parte

degli elettroni. L’energia degli ioni che bombardano le pareti è quindi

essa stessa pari a qualche eV.

Il potenziale uniforme nella regione centrale della scarica è

detto invece potenziale di plasma, e la regione centrale è chiamata

glow.

Fig. 1.2 Formazione dello sheath 1.4 Modello circuitale dello sheath

La frequenza del voltaggio applicato ha un valore intermedio tra

la frequenza di oscillazione degli elettroni e quella degli ioni. Sotto

11

queste ipotesi si assume che lo sheath sia puramente capacitivo e la

regione centrale sia equipotenziale. Ragionando in termini di circuiti

elettrici, questo equivale a dire che gli sheath hanno resistenza infinita

e il glow abbia resistenza trascurabile. Il plasma RF può essere

rappresentato come in figura:

Fig. 2.2 Circuito equivalente del modello dello sheath capacitivo

Cs1 e Cs2 sono capacità che rappresentano gli sheath, D1 e D2

sono diodi che fanno in modo che il potenziale di plasma Vs sia sempre

positivo rispetto agli elettrodi.

Nell'ipotesi in cui il comportamento dello sheath sia puramente capacitivo, il potenziale di plasma Vp(t) è espresso come: Vp (t) = V p + ∆Vp sen(ωt) (1.5)

dove V p è la media temporale del potenziale di plasma e ∆Vp

l’ampiezza della modulazione dello stesso potenziale. ∆Vp è

determinato dalla partizione della tensione Vrf secondo la:

∆Vp =Cs1

Cs1 + Cs2

Vrf (1.6)

dove Cs1 e Cs2 sono le capacità che rappresentano rispettivamente

lo sheath di fronte all’elettrodo sotto tensione e quello di fronte

12

all’elettrodo a terra. Il potenziale di plasma istantaneo è sempre più

elevato rispetto al potenziale agli elettrodi, ovvero: Vp ,max = V p + ∆Vp ≥ Vrf

Vp ,min = V p − ∆Vp ≥ 0 (1.7)

Una volta ogni periodo il potenziale di plasma uguaglia tanto il

potenziale dell’elettrodo sotto tensione, quanto quello dell’elettrodo a

terra per bilanciare la corrente netta attraverso lo sheath, ovvero le 1.7

diventano uguaglianze.

I valori delle capacità si sheath, anch’essi mediati nel

tempo, sono determinati dall’area degli elettrodi ma sono anche

influenzati dal modulo delle cadute di potenziale nei due sheaths.

1.5 Diffusione ambipolare

Ai campi di sheath e bulk occorre aggiungere un terzo

elemento, il cosiddetto campo ambipolare, che è dovuto allo stesso tipo

di considerazioni, ma in assenza di campo applicato.

Si fa nel plasma l’ipotesi di congruenza, ovvero che il

flusso di elettroni e quello di ioni, dovuti tanto alla diffusione nel

campo elettrico quanto a quella per gradiente di concentrazione, siano

uguali ovunque. Questa ipotesi continua a valere in presenza di

collisioni ionizzanti, che producono un ugual numero di entrambe le

specie. Dato che gli elettroni sono più leggeri, e tenderebbero a lasciare

il plasma più velocemente, deve nascere un campo elettrico che

mantenga l’equilibrio locale del flusso.

13

Questo campo ambipolare ha l’effetto di impedire

l’accelerazione nel campo di sheath agli elettroni di bassa energia.

Questi elettroni, confinati nella regione centrale del plasma, saranno

soggetti alla sola azione del campo bulk.

Aumentando il potenziale applicato all’elettrodo aumenta il

potenziale di accelerazione degli ioni verso le pareti. La caduta di

potenziale in prossimità degli elettrodi dipende dalla superficie degli

stessi: in particolare il potenziale decade in modo inversamente

proporzionale alla superficie elettrodica. Anche la pressione influenza

l’energia del bombardamento ionico poichè un suo aumento determina

un maggior numero di urti nello sheath che rallentano gli ioni.

Inoltre, se la frequenza della scarica è minore della

frequenza con cui gli ioni attraversano lo sheath, questi ultimi sono

sensibili al potenziale istantaneo e il bombardamento ionico risulta

estremamente energetico. Quando questa condizione non è verificata

gli ioni non riescono a seguire le oscillazioni del potenziale e sono

sensibili solo al potenziale medio. Quindi, operando nella regione dei

MHz il moto delle particelle pesanti è controllato dal campo medio.

Considerando insieme gli effetti dovuti alla diffusione

ambipolare e ai campi elettrici di bulk e di sheath si osserverà un

profilo di concentrazione ionica che non cambia molto nel tempo a

causa dell’elevata inerzia degli ioni rispetto a quella degli elettroni. Per

gli elettroni si osserva invece una oscillazione.

14

1.6 Riscaldamento elettronico

Una delle questioni chiave riguardo la struttura e la

dinamica delle scariche RF è il meccanismo di riscaldamento

elettronico e il corrispondente meccanismo di ionizzazione o di

sostentamento della scarica.

Sono stati identificati tre meccanismi [7] di riscaldamento

elettronico nelle scariche glow:

1. gli elettrodi emettono elettroni secondari, e questi elettroni si

riversano negli elevati campi di sheath ed entrano nel plasma con

notevole energia;

2. gli elettroni nel plasma sono riscaldati dal campo elettrico

oscillante di sheath;

3. i campi elettrici all’interno del plasma scaldano gli elettroni che

trasportano la corrente di scarica.

Il secondo meccanismo è una prerogativa esclusiva delle scariche

RF ed è stato diffusamente studiato [8,9,10,11] . Si può distinguere tra

due casi limite di riscaldamento per oscillazione dello sheath:

riscaldamento stocastico in cui gli elettroni non subiscono collisioni

con le particelle neutre mentre interagiscono col campo elettrico

oscillante di sheath, e moto dominato dalle collisioni. Il primo è detto

riscaldamento stocastico di sheath e il secondo riscaldamento ohmico

di sheath.

Il riscaldamento ohmico è detto anche regime «wave-riding»[7].

Infatti in questo caso, una volta che l’elettrone è stato «intrappolato»

nel campo elettrico durante l’espansione dello sheath, «cavalca» lo

sheath in espansione, e percepisce un campo elettrico proporzionale

alla velocità di espansione dello sheath.

15

Oltre alla distinzione tra i diversi tipi di riscaldamento

dovuti all’oscillazione dello sheath, va fatta anche una distinzione tra

scariche sostenute dall’emissione di elettroni secondari e scariche

sostenute dal riscaldamento elettronico nello sheath o nel bulk. Sono

dette scariche alfa quelle in cui il meccanismo dominante di

ionizzazione è il riscaldamento degli elettroni all’interno del plasma e

scariche gamma quelle in cui la ionizzazione avviene prevalentemenet

ad opera degli elettroni secondari.

E’ possibile passare da un regime di scarica ad un altro

variando parametri quali potenza, pressione e frequenza di scarica. Si

va dal regime alfa al regime gamma aumentando potenza e pressione,

mentre le scariche con frequenza minore di 1 MHz sono sostenute

prevalentemente dagli elettroni secondari[12].

17

CAPITOLO 2 CARATTERISTICHE MATEMATICHE E

FISICHE DEL MODELLO

2.1 Modello numerico per la simulazione: PIC/MCC

2.1.1 PIC

La maggior parte dei plasmi di non equilibrio di interesse è

caratterizzata da un gran numero di elettroni nella sfera di Debye: in

queste condizioni l’interazione fra particelle cariche può essere

approssimata da una interazione tra particella e campo di carica

spaziale . Questo è proprio il punto di vista del metodo PIC.

(Particle In Cell). L’ approccio del PIC consiste nel risolvere

l’equazione di Vlasov:

∂∂t

+ v ⋅ ∇ r + e∇ rϕm

⋅ ∇ v

f r ,v ,t( )= 0 (2.1)

e l’equazione di Poisson:

∇2ϕ(r) = −ρε0

(2.2)

assumendo una soluzione a N particelle:

f (r,v,t) = w δ(r − r i )

i=1

N

∑ δ(v − v i)

ρ(r , t) = −ew S(r − r i) + ρ' (r )i =1

N

∑ (2.3)

18

dove ri e vi sono le posizioni delle singole particelle, che

soddisfano le equazioni di Newton, δ è la funzione di Dirac, S è il

fattore di forma delle particelle (una δ con larghezza finita) e ρ’ è la

distribuzione di carica dovuto alle altre particelle.

Le equazione vengono risolte per un gran numero di elettroni

(104-106 ) tenendo conto del campo elettrico locale risultante

dall’interpolazione locale della carica spaziale all’interno di una cella

di un reticolo matematico. E’ praticamente impossibile effettuare la

simulazione con il numero reale di particelle, per cui le particelle

simulate sono in realtà ‘superparticelle’ con un peso matematico w

(pari al rapporto tra particelle simulate e particelle reali), con carica e

massa uguali rispettivamente a –ew ed mw. In questo modo il rapporto

q/m rimane invariato e le superparticelle si muovono come elettroni

sotto l’influenza del campo elettrico e del campo magnetico. Dopo un

passo di calcolo delle equazioni del moto, viene determinata la carica

elettrica in ogni cella del reticolo dal numero di elettroni presenti in

essa, opportunamente pesati. Dalla densità di carica elettrica,

risolvendo l’equazione di Poisson all’interno del reticolo stesso, si

ricavano il potenziale e il campo elettrico.

L’assegnazione della carica a ciascun punto, che determina

la forma di S ) consiste in una interpolazione tra i punti della griglia più

vicini alla particella. Nel caso monodimensionale esistono due ordini di

assegnazione. Nell’assegnazione di ordine zero, o NGP (nearest grid

point) si conta il numero di particelle in ciascuna cella e si assegna

questo numero al punto di griglia rappresentativo della cella.

Una assegnazione che, sebbene più complessa della NGP, è

molto meno rumorosa, è l’assegnazione di prim’ordine, o CIC (cloud in

cell). In questo caso si assegna la carica della particella ai due punti più

vicini alla cella in frazioni dipendenti dalla posizione della particella

19

stessa; ovvero, con i indice del punto di griglia e q carica della

particella:

( )

xxixqq

xxxiqq

xixxi

i

i

∆∆−=

∆−∆+=

∆+≤≤∆

+1

)1(1

(2.4)

Le stesse formule devono essere usate per interpolare il campo

elettrico che agisce su ciascuna particella per soddisfare la

conservazione del momento per valori molto piccoli del passo

temporale.

Le equazioni di Newton nel PIC vengono risolte utilizzando un

semplice schema per posizione e velocità, il cosiddetto metodo

Leapfrog[1,13]:

tvrr

trEmqvv

∆+=

∆+=

''

)(' (2.5)

La differenza rispetto al metodo esplicito di Eulero è

semplicemente la sostituzione di v’ con v nella seconda equazione.

Bisogna comunque notare che r è calcolato al tempo t e v è calcolato al

tempo t+∆t.

Una dettagliata descrizione del PIC si trova in Birdsall and

Langdon (1985,1991) [13,14], Vahedi,[23]Cohen[16]Longo et al.[2].

20

2.1.2 MCC

Sulla base di quanto detto finora riguardo al PIC è evidente come

esso sia in grado di riprodurre solo le interazioni tra le particelle

cariche. Ma una descrizione accurata di un plasma collisionale, quale è

quello RF, non può trascurare gli urti tra particelle cariche e specie

neutre. A tal fine il Particle in Cell è stato accoppiato con il metodo

Monte Carlo secondo lo schema a blocchi in figura 2.1.

Ad ogni passo di calcolo il programma risolve le equazioni del

campo e muove tutte le particelle.

Il Monte Carlo viene inserito nel ciclo calcolando la

probabilità che una particella subisca una collisione; in caso

affermativo la particella dopo la collisione rientrerà nel ciclo principale

con una nuova velocità.

��

��

��

Integration of equations of motion, moving particles

Fi v i’ x i

��

��

Monte Carlo Collisionsv i’ v i

��

Weighting

(x ,v)i ρj

Integration of field equations of grid

Ej

Weighting

FiEj

∆∆∆∆t

ρj

Figura 2.1 Schema a blocchi dell’accoppiamento PIC/MCC

Inserendo le collisioni il numero totale di particelle varia.

Poichè l’errore relativo su quantità macroscopiche, dovuto a

fluttuazioni statistiche, diminuisce con il numero di particelle*, è

* L’errore decresce comunque piuttosto lentamente col numero n di particelle (tipicamente come n-1/2)

21

necessario mantenere elevato quest’ultimo nella simulazione. Tuttavia,

il tempo di calcolo pone un limite al numero di particelle simulate. Nel

codice è quindi inserito un controllo che elimina metà particelle a caso

e moltiplica per due il peso statistico qualora le particelle diventino più

del doppio di quelle iniziali, mentre le raddoppia dimezzando il peso

quando sono ridotte a meno della metà rispetto a quante erano

all’inizio.

Il metodo Monte Carlo è basato sul calcolo dei tempi di

volo per la particella testata tra due successive collisioni con le altre

particelle. Questi tempi sono calcolati generando numeri casuali presi

da una distribuzione appropriata modellata secondo la fisica che regola

i processi. Si suppone che i processi siano binari ed istantanei.

I numeri casuali possono essere calcolati utilizzando il

seguente teorema della teoria delle probabilità: dato un insieme di

numeri casuali ri distribuiti uniformemente tra zero ed uno, un insieme

di numeri casuali yi distribuiti secondo la funzione f(y) nell’ intervallo

(a,b) è dato da:

f (x)dxa

yi

∫ = ri f (x)dxa

b

∫ (2.6)

Dopo aver calcolato l’integrale, l’equazione 2.3 diventa una

equazione in yi.

Dato che l’insieme di numeri ri può essere fornito per mezzo di

una routine generatrice di numeri casuali dal calcolatore, si può

simulare ogni processo fisico casuale caratterizzato da distribuzioni

note di probabilità.

Il tempo che trascorrerà fino alla successiva collisione va

calcolato conoscendo la frequenza di collisione in funzione della

velocità elettronica.Il modo migliore per risolvere questo problema è di

22

introdurre una specie nulla in modo che in una collisione tra elettrone e

specie nulla (collisione nulla) la velocità elettronica rimane

invariata[17].

Figura 2.2 Probabilità di collisione e collisione nulla

La sezione d’urto per la collisione nulla è scelta al fine di

ottenere una frequenza totale di collisione νtot costante, data da :

ν tot = maxx ∈V , 0≤ε ≤ ε lim

2εme

12

Np(x)σ p(ε)p

∑

(2.7)

dove la somma è fatta su ogni processo di collisione p, σp(ε) è la

sezione d’urto, Np è la densità del partner collisionale, V è l’estensione

spaziale del sistema e εlim è il limite della scala di energia cinetica usata

per calcolare le sezioni d’urto. Si è usata la seguente approssimazione:

dato che gli elettroni si muovono molto più velocemente delle

23

molecole, possiamo semplificare il calcolo delle frequenze di collisione

assumendo che le particelle obiettivo siano ferme.

Con una frequenza di collisione costante, i tempi di

collisione sono distribuiti secondo [18].

)exp()( ttf tottot νν −= (2.8)

Il tempo prima della successiva collisione è ottenuto usando

l’equazione 2.3 nella forma

ην

ln1

totct −= (2.9)

dove η è un elemento di un insieme di numeri casuali

uniformemente distribuiti tra zero e uno. Il tipo di collisione (incluso

quello nullo) che si verifica dopo il tempo di volo libero viene scelto in

base alle corrispondenti frequenze di collisione : sia

ν k(ε ) =2εme

12

σ k(ε )nk(r ) (2.10)

il contributo dell’elettrone k-esimo alla frequenza totale di

collisione. Dopo aver generato il numero casuale r, il processo di

collisione scelto sarà l’n-esimo con n tale che

νk

ν totk=1

n−1

∑ < r <νk

ν totk=1

n

∑ (2.11)

24

Figura 2.3 Selezione del processo collisionale 2.2 Modifica del ciclo temporale

Consideriamo le modifiche al normale ciclo tempo

temporale che sono richieste al fine di incorporare in Monte Carlo nel

PIC.

Questa modifica è fondamentale per poter permettere che

le collisioni avvengano anche più volte all’ interno del passo di calcolo

PIC [19].Quest’ultimo è il tempo che trascorre tra due successivi

aggiornamenti del campo.

25

Figura 2.4 Modifica del ciclo temporale

Dopo aver calcolato per la prima volta il tempo di volo

libero, le particelle vengono fatte muovere utilizzando il metodo

Leapfrog in cui, però, il passo di calcolo ∆t è sostituito dal tempo di

volo libero tc. A questo punto, se tc è minore del passo di calcolo ∆t la

particella subisce una collisione e tc viene sottratto a ∆t. Si determina il

tipo di collisione e si calcola un nuovo tempo di volo libero che viene

confrontato con la parte restante del passo di calcolo.

Se, ancora una volta, tc risulterà più piccolo, la particella subirà

una nuova collisione e si ripeterà la procedura già vista. Diversamente,

le particelle vengono fatte muovere per la parte restante del passo di

calcolo ∆t e l’eccesso di tc viene immagazzinato come tempo di volo

libero iniziale nel successivo passo di calcolo.

26

2.3 Accoppiamento del modello con la cinetica

A questo punto è utile una precisazione: durante i processi

collisionali, le ∼100000 particelle della simulazione sono considerate

come un campione preso dalla reale distribuzione di tutte le particelle

nel plasma, mentre durante il calcolo del campo le particelle assumono

lo status di superparticelle, in modo da apportare il giusto contributo di

carica.

Un’altra importante caratteristica delle particelle neutre (e

del gas in generale) nel modello PIC/MCC è che si assume che esse

siano distribuite uniformemente nello spazio, e quindi non sono seguite

come particelle. Vedremo invece adesso come, tramite la cinetica

reattiva e diffusiva, le neutre saranno considerate come particelle e la

loro cinetica sarà accoppiata alla dinamica delle particelle cariche.

Questo accoppiamento è insito nella concezione di plasma:

l’equazione 2.4 mostra che la frequenza di collisione dipende dalla

composizione chimica del gas, ovvero la cinetica chimica influenza la

dinamica. Allo stesso tempo nel prossimo capitolo si vedrà come la

velocità di reazione dipenda fortemente dalla funzione di distribuzione

dell’energia elettronica*.

L’accoppiamento che si deve realizzare deve essere

autoconsistente, ovvero la dinamica e la cinetica chimica devono essere

risolte tenendo conto della loro reciproca connessione.

Il diagramma in figura 2.5 illustra la linea guida dell’intero

modello

Durante i calcoli, le densità delle differenti specie devono

essere aggiornate risolvendo le equazioni appropriate. In particolare, le * La forma non Maxwelliana della eedf infatti non permette che il coefficiente di velocità di una reazione venga espresso in termini di parametri macroscopici quali ad esempio la temperatura elettronica.

27

densità delle particelle cariche sono ottenute direttamente dalla

simulazione PIC/MCC mentre quelle dei neutri vengono aggiornate

risolvendo le equazioni reattive e diffusive tenendo conto del fatto che i

parametri che rientrano in queste equazioni non sono costanti ma

dipendono dalla funzione distribuzione locale dell’energia degli

elettroni, che pertanto deve essere calcolata per prima.

Figura 2.5 Schema a blocchi del modello

Durante la simulazione PIC/MCC, la eedf viene campionata su di

una griglia spaziale impiegata anche per la cinetica diversa da quella

utilizzata nella risoluzione della equazione di Poisson. Le due griglie

hanno una densità di punti diversa: la griglia per la cinetica ha 50 punti

mentre quella per l’equazione di Poisson ne ha 400. Questo permette di

ridurre le fluttuazioni statistiche nel calcolo dei coefficienti di velocità

delle reazioni chimiche considerate.

Composizione del gas

Densità particelle cariche

CONDIZIONI AL

CONTORNO

DINAMICA

DELLE PARTICELLE

CARICHE

EQUAZIONE DI

POISSON

Carica spaziale

Campo elettrico

EQUAZIONI DI

REAZIONE E

DIFFUSIONE

eedf

28

Come vedremo in seguito, anche la temperatura viene

calcolata sulla griglia per la cinetica, come è da attendersi essendo in

questo lavoro interessati alla temperatura del gas.

Il metodo utilizzato consiste nel risolvere la cinetica, fino

al raggiungimento dello stato stazionario, a tempi diversi tk durante la

simulazione PIC/MCC. Questa tecnica adiabatica sebbene opposta alla

più naturale (essendo i tempi caratteristici di plasma molto più piccoli

rispetto a quelli della cinetica) si rivela molto efficace poichè i tempi di

calcolo necessari per la risoluzione delle equazioni della cinetica e

della diffusione sono trascurabili in confronto a quelli richiesti per la

descrizione delle particelle cariche nel plasma.

Ad ogni tk i valori dei coefficienti di velocità e le densità

elettroniche sono quelle che si ottengono dalla simulazione PIC/MCC,

mediati nel periodo tk-tk-1.

L’utilizzo di questo metodo di accoppiamento rende quindi

impossibile una descrizione temporale sensata prima che venga

raggiunta la stabilità.

29

2.4 Dinamica delle collisioni

Per quanto riguarda il trattamento delle collisioni binarie, se

consideriamo due particelle che collidono, con masse m1 ed m2 e

velocità v1 e v2 , il trattamento è semplificato scomponendo le velocità

in velocità del centro di massa vc e velocità relativa vr, cioè

rc

rc

vmm

mvv

vmm

mvv

21

12

21

11

+−=

++=

(2.12)

Il processo di collisione può soltanto cambiare la velocità

reòativa, secondo il guadagno di energia interna ∆ε e gli angoli di

scattering θ e ϕ. In questo modo può essere introdotto ogni processo

collisionale che lascia invariato il numero totale di particelle: la

generalizzazione ad altri processi quali la ionizzazione o l’attachment

elettronico non presenta particolari difficoltà.

Per affrontare il problema del trasferimento di energia

definiamo due sezioni d’urto: la sezione d’urto totale elastica

σ te = 2π σ ϑ ,ε( )0

π

∫ sinϑdϑ (2.13)

e la sezione d’urto di momentum transfer

σm = 2π σ ϑ,ε( )0

π

∫ 1 − cosθ( )sinϑdϑ (2.14)

30

dove (1 - cosθ) è la frazione di momento iniziale mv persa dalla

particella incidente e σ(θ,ε) è la sezione d’urto differenziale.

Nella collisione elastica di un proiettile di massa m1 e

velocità v1 con un bersaglio stazionario di massa m2 , la conservazione

del momento lungo la direzione di v1 e perpendicolarmente ad essa e la

conservazione dell’energia possono essere scritte nel sistema

laboratorio come

222

211

211

222111

22211111

'21'

21

21

sin'sin'0cos'cos'

vmvmvm

vmvmvmvmvm

+=

−=+=

θθθϑ

(2.15)

dove gli apici denotano i valori dopo la collisione. Risolvendo il

sistema eliminando v1’ e θ1 si ottiene

22

221

21211

222 cos

)(4

21'

21 θ

mmmmvmvm

+= (2.16)

Guardando l’equazione 2.13 si vede che la frazione di energia

persa dal proiettile nel sistema laboratorio è

=Lζ 22

221

21 cos)(

4θ

mmmm

+ (2.17)

Trasformando l’angolo per portarci nel sistema del CM si ottiene

)cos1()(

42

21

21 Θ−+

=mmmm

Lζ (2.18)

31

A questo punto si media sulla sezione d’urto differenziale per

ottenere la perdita media:

te

mL mm

mmσσζ 2

21

21

)(2

+=

Θ (2.19)

Gli elettroni quindi trasferiscono poca energia nelle collisioni

elastiche con particelle pesanti, facendo in modo che Te>>Ti in una

tipica scarica. Invece, per m1≈m2 si ottiene ζ L Θ =12 , che comporta un

forte scambio di energia tra particelle pesanti e quindi una temperatura

comune.

33

CAPITOLO 3 REAZIONI ED ENERGIA

3.1 Specie neutre presenti nel modello

Le specie neutre presenti nel nostro modello sono sette:

O2(X3ΣΣΣΣ-

g) stato fondamentale

O2(a1∆∆∆∆g) primo stato metastabile[O2*]

O2(b1ΣΣΣΣ+g) secondo stato metastabile[O2

**]

O3

O(3P) stato fondamentale

O(1D) primo stato metastabile[O*]

O(1S) secondo stato metastabile[O**]

All’ inizio della simulazione queste sette specie sono state

immesse con le seguenti densità e relative pressioni parziali* : specie densità(m-3) pressione parziale O2 3.42×1021 1.062×10-1 O2* 1.66×1020 5.135×10-3 O2** 1.34×1019 4.160×10-4 O3 9.50×1017 2.950×10-5 O 2.44×1019 7.580×10-4 O* 5.12×1017 1.590×10-5 O** 3.53×1015 1.100×10-7

Inoltre il numero iniziale di elettroni è stato posto uguale a

5x1015

Una volta partita la simulazione, gli elettroni vengono accelerati

dalla tensione alternata applicata e, collidendo con le specie pesanti,

* Si è usata per le pressioni parziali la formula p=ρRT, con R costante dei gas e T temperatura del gas

34

vanno a eccitarle o diseccitarle ma anche a ionizzarle, producendo altri

protagonisti della simulazione: gli ioni O2+,O-,O2

-,O+.

A loro volta le componenti pesanti, ionizzate o meno,

interagiranno tra di loro in reazioni di ricombinazione, attachment e

detachment che vedremo in dettaglio in seguito.

3.2 Sezioni d’urto e costanti di velocità

Tutte le reazioni che vedono un elettrone collidere con l’ossigeno

neutro atomico o molecolare sono caratterizzate da una sezione d’urto

dipendente dall’energia.

In questo modo si ha una descrizione a livello microscopico del

sistema, per cui una reazione ha una maggiore o minore proabilità di

verificarsi (corrispondente ad un intervallo più o meno ampio entro cui

deve cadere il numero casuale generato) a seconda dell’ energia dell’

elettrone incidente*.

Le reazioni tra gli ioni e le reazioni da impatto elettronico

con gli ioni e con l’ozono sono invece descritte dalla costante di rate K,

definita in modo che il rate della reazione (il numero di reazioni

elementari del tipo e-+A→prodotti che avvengono per unità di tempo e

di volume) sia dato da:

AenKnv = (3.1)

dove nA è la densità molare di A e ne è la densità numerica degli

elettroni.

* vedremo che questa energia è presa da un profilo non maxwelliano della eedf (electron energy distribution function)

35

La costante di rate ha pertanto le dimensioni di cm3 mol-1 sec-1 e

ha il significato fisico di una integrazione sulla distribuzione di energia

elettronica e sulla sezione d’urto del processo:

K =2

me

12

εf (ε)σ(ε )dεεth

∞

∫ (3.2)

dove f è la eedf, me è la massa dell’elettrone e σ è la sezione

d’urto del processo elementare.

Il secondo membro dell’ equazione 3.2 dipende dalla eedf e a

causa della forma decisamente non Maxwelliana di quest’ultima nei

plasmi freddi non può essere scritta come una funzione dei parametri

macroscopici quali la temperatura elettronica.[20]. I dati da noi

utilizzati provengono infatti da dati sperimentali.

Da quanto specificato nel secondo capitolo, le reazioni per

cui non è disponibile il set di sezioni d’urto non possono essere trattate

con il metodo delle collisioni nulle.

Si procede allora calcolando la probabilità che il

protagonista della reazione scompaia nel passo di calcolo ∆t attraverso

la seguente relazione [14]:

( )tP coll∆−−= νexp1 (3.3)

dove νcoll è la frequenza di collisione. Per ogni particella νcoll sarà

la somma della frequenza di collisione per tutti i processi di cui essa è

protagonista. I processi che noi consideriamo sono tutti processi del

secondo ordine, per cui la definizione di frequenza è:

36

νcoll = Ki S[ ]i∑ (3.4)

dove la somma è estesa a tutti i processi di cui la particella è protagonista, con Ki costante di velocità del processo i-esimo e [ ]S

densità numerica della specie compagno.

Si estrae quindi un numero casuale compreso tra zero ed

uno e lo si confronta con P. Se è minore, si verifica un processo di

perdita della particella protagonista. 3.3 Reazioni da impatto elettronico

Sono riportate qui di seguito le reazioni presenti nella

simulazione, insieme alle rispettive sezioni d’urto o costanti di rate e

alla fonte sperimentale. Inoltre in corrispondenza di ogni reazione si

trova il guadagno o la perdita di energia nel plasma che essa comporta. 3.3.1 Eccitazioni verso stati metastabili • e- + O2

* → O2** + e- [21] perdita di 0.64 eV

• e- + O2 → O2

* + e- [21] perdita di 0.98 eV • e- + O2 → O2

** + e- [21] perdita di 1.627 eV • e- + O → O*+ e- [22] perdita di 1.97 eV • e- + O → O** + e- [22] perdita di 4.18 eV

Le eccitazioni generano diversi stati metastabili dell’ossigeno

molecolare (a1∆g, b1Σg, c1Σ-u, A3Σ+

u) e atomico e sono importanti

meccanismi di perdita di energia elettronica, oltre a permettere la

ionizzazione dagli stati metastabili.

37

3.3.2 Diseccitazioni • e- + O2

* → O2 + e- [21] guadagno di 0.98 eV • e- + O2

** → O2* + e- [21] guadagno di 0.64 eV

• e- + O2

** → O2 + e- [21] guadagno di 1.627 eV • e- + O* → O + e- [22] guadagno di 1.97 eV • e- + O** → O + e- [22] guadagno di 4.18 eV

Le diseccitazioni sono delle collisioni superelastiche, in cui un

elettrone urta una specie pesante che si trova in uno stato eccitato e

guadagna l’energia di soglia del processo, mentre la specie pesante si

diseccita.

Le sezioni d’urto per questi processi vengono calcolate a

partire da quelle dei corrispondenti processi anelastici attraverso la

seguente relazione:

( ) ( ) ( )εεσεεεεσ ∆+∆+= angg *sup (3.5)

dove g e g* sono, rispettivamente, le degenerazioni dello stato

fondamentale ed eccitato della specie pesante e ∆ε è l’energia di soglia

del processo.

38

Figura 3.1 Diseccitazione(1) ed Eccitazione(2) tra lo stato fondamentale ed il

primo metastabile dell’O2

Figura 3.2 Diseccitazione(1) ed Eccitazione(2) tra lo stato fondamentale ed il secondo metastabile dell’O2

10-3

10-2

10-1

100

101

σ(cm -16)

Energia(eV)

1

2

1 e- + O2

*(Σ) -> O

2 + e-

2 O2 + e- -> e- + O

2

*(Σ)

10-3

10-2

10-1

100

6 10-18 10

-110

03 10

05 10

07 10

0

σ(cm -16)

Energia(eV)

1 e- + O2

*(∆ ) -> O

2 + e-

2 O2 + e- -> e- + O

2

*(∆ )

1

2

39

Figura 3.3 Diseccitazione(1) ed Eccitazione(2) tra i primi due stati metastabili dell’O2

Figura 3.4 Diseccitazione(1) ed Eccitazione(2) tra lo stato fondamentale ed il primo metastabile dell’ossigeno atomico

10-3

10-2

10-1

100

10-3

10-2

10-1

100

101

σ(cm -16)

Energia(eV)

1 2

1 O2

*(Σ)+e- ->e- + O

2

*(∆)

2 e- + O2

*(∆) -> O

2

*(Σ)+e-

10 -2

10 -1

10 0

10 1

10 -3 10 -2 10 -1 10 0 10 1

σ(cm -16)

Energia(eV)

1

2

1 e-+ O*(1D) -> O(3P) + e-

2 O(3P) + e- -> e-+ O*(1D)

40

Figura 3.5 Diseccitazione(1) ed Eccitazione(2) tra lo stato fondamentale ed il secondo stato metastabile dell’ossigeno atomico

3.3.3. Dissociative attachments O2 • e- + O2 → O + O- [21] perdita di 4.2 eV • e- + O2

* → O + O- [21] perdita di 3.22 eV • e- + O2

** → O + O- [21] perdita di 2.573 eV

Viene riportata la sezione d’urto solo per il primo dei tre

processi. Le sezioni d’urto per gli altri due si ottengono traslando le

curve verso valori decrescenti dell’energia di 0.98 e 1.627 eV, pari alle

distanze in eV tra il ground state e i primi due livelli eccitati, rendendo

ovviamente più agevole il processo.

10-2

10-1

100

101

10-3

10-2

10-1

100

101

σ(cm -16)

Energia(eV)

1

2

1 e-+ O

*(1D) -> O(

3P) + e

-

2 O(3P) + e

- -> e

-+ O

*(1D)

41

La reazione di attachment dissociativo è il principale

meccanismo di creazione di ioni negativi e di perdita di energia

elettronica nel bulk. L’elettrone incidente perde l’energia di soglia di

4.2 eV ed è assorbito dalla molecola di ossigeno per formare uno ione

O2- che si dissocia in O e O-. Si assume che la transizione elettronica

sia rapida su scala nucleare [24] e ciò che rimane dell’energia incidente

dell’elettrone è divisa tra i frammenti. Dato che la sezione d’urto per

questo processo ha il picco intorno ai 6.5 eV i frammenti avranno

energie dell’ordine di 1-2 eV.

3.3.4 Ionizzazioni • e- + O2 → O2

+ + 2e- [21] perdita di 12.06 eV • e- + O2

* → O2+ + 2e- [21] perdita di 11.08 eV

• e- + O2** → O2

+ + 2e- [21] perdita di 10.43 eV • e- + O → O+ + 2e- [22] perdita di 13.61 eV • e- + O* → O+ + 2e- [25] perdita di 11.65 eV • e- + O** → O+ + 2e- [25] perdita di 9.43 eV

La ionizzazione è l’unica reazione che produce ioni positivi nel

sistema. Si crea una coppia elettrone-ione e il bilancio energetico si

scrive:

ionNinciprodscat εεεεεε −+=++ (3.6)

dove εinc, εscat e εprod sono le energie rispettivamente dell’elettrone

incidente prima e dopo la collisione e dell’elettrone rilasciato dallo

ione; εi e εN somo le energie dello ione creato e della neutra bersaglio e

42

εion è l’energia di soglia della ionizzazione. A causa dell’elevato valore

del rapporto tra le masse dello ione e dell’elettrone si può supporre che

l’elettrone incidente “strappi” un elettrone dalla particella neutra e che

questa diventi uno ione, continuando la sua traiettoria indisturbata.

Possiamo riscrivere l’equazione 3.6 come

Ni

ionincprodscat

εεεεεε

=

−=+ (3.7)

L’energia dopo la collisione viene distribuita tra i due

elettroni generando un numero casuale R ∈[0,1]. L’energia

dell’elettrone prodotto sarà:

2ioninc

prod R εεε −= (3.8)

e la 3.7 dà l’energia dell’elettrone scatterato.

3.3.5 Eccitazione dissociativa •••• e- + O2 → e- + O + O (3P) [21] perdita di 14.7 eV

Nelle dissociazioni l’elettrone perde l’energia di soglia di

6, 8.4, 9.97 o 14.7 eV ed eccita elettronicamente la molecola di O2 in

uno stato che si dissocia in due O. Queste reazioni sono i principali

meccanismi responsabili di produzione di ossigeno atomico nel bulk.

Vengono riportate le sezioni d’urto solo per i processi che

hanno come specie di partenza lo stato fondamentale dell’O2. Le

sezioni d’urto per i processi che hanno come specie bersaglio i primi

43

due stati metastabili si ottengono traslando le curve verso valori

decrescenti dell’energia di 0.98 e 1.627 eV, pari alle distanze in eV tra

il ground state e i primi due livelli eccitati, rendendo ovviamente più

agevole il processo.

Figura 3.6 Sezioni d’urto per attachment dissociativo, ionizzazione e eccitazione dissociativa dell’O2

Figura 3.7 Sezioni d’urto per le ionizzazioni dei primi tre stati dell’ossigeno

atomico

10 -3

10 -2

10 -1

10 0

10 1

10 -2 10 -1 10 0 10 1 10 2

σ(cm -16)

Energia(eV)

1

2

3

1 Attachment dissociativo (4.2eV)2 Ionizzazione (12.06 eV)3 Eccitazione dissociativa(14.7 eV)

1

10 -1

10 0

10 0 10 1 10 2

σ(cm -16)

Energia(eV)

1

2 3

1 e- + O*(1S) -> 2e- + O+

2 e- + O*(1D) -> 2e- + O+

3 e- + O(3P) -> 2e- + O+

44

3.3.6 Perdite di 4.5, 6.0, 8.4 e 9.97 eV

L’identificazione delle suscritte perdite di energia con processi

specifici non è definita con precisione e non c’è accordo tra i vari

sperimentatori[26].I candidati più probabili per queste perdite di

energia sono[15]: • e- + O2 → e- + O2(c1Σ-

u, A3Σ+u) [21] perdita di 4.5 eV

• e- + O2 → e- + O + O [21] perdita di 6.0 eV • e- + O2 → e- + O + O* [21] perdita di 8.4 eV • e- + O2 → e- + O* + O* [21] perdita di 9.97 eV

Processi analoghi, con sezioni d’urto e perdite di energia

opportunamente shiftate, si hanno per impatto di elettroni con i primi

due stati metastabili dell’O2.

Figura 3.8 Sezioni d’urto per i processi di energy loss

10 -3

10 -2

10 -1

10 0

10 1

10 -2 10 -1 10 0 10 1 10 2

σ(cm -16)

Energia(eV)

3

4

1

2

1 4.5 eV Loss

2 6.0 eV Loss

3 8.4 eV Loss

4 10 eV Loss

45

3.3.7 Momentum transfer • e- + O2 → e- + O2

In questo caso l’equazione 2.16 ci dice che la perdita di energia

media è dell’ordine di 10-4 a causa del rapporto tra le masse, per cui in

una scarica tipica si ha Te>>Ti. Sebbene questa perdita di energia sia

bassa, questa collisione elastica è l’unico meccanismo di perdita di

energia per gli elettroni con energia minore delle soglie di ionizzazione

ed eccitazione, ed è anche abbastanza frequente a causa della sua

elevata sezione d’urto.

Figura 3.9 Sezione d’urto per il momentum transfer e-/O2

10-2

10-1

100

101

102

10-2

10-1

100

101

102

103

104

σ(cm -16)

Energia(eV)

e- + O

2 -> e

- + O

2

46

3.3.8 Dissociative attachment O3 • e- + O3 → O-

2 +O perdita di 0.61 eV K = 6.02 × 108 cm3/mol sec [27] 3.3.9 Ricombinazioni • 2e- + O2

+ → O2 + e- guadagno di 12.14 eV K = 2.61 × 109 cm3/mol sec [28] • 2e- + O+ → O + e- guadagno di 11.05 eV K = 2.61 × 109 cm3/mol sec [28]

Tra le possibili eccitazioni da impatto elettronico abbiamo

tralasciato le eccitazioni rotazionali e vibrazionali. Per quanto riguarda

le prime, hanno una sezione d’urto minore di alcuni ordini di grandezza

rispetto agli altri processi e sono relativamente poco importanti in un

quadro globale. Le eccitazioni vibrazionali invece sono molto

importanti ai fini della funzione di distribuzione dell’energia

elettronica. Infatti questi processi hanno una sezione d’urto diversa da

zero solo per piccoli intervalli di energia, e questo comporta una forma

della eedf decisamente non maxwelliana. Ai fini della raccolta di

energia, però, si vede che il rilassamento vibrazionale è molto rapido

nell’ossigeno[29] e che tutta l’energia che andrebbe ad eccitare

vibrazionalmente la molecola viene rapidamente rilasciata

47

prevalentemente mediante la reazione O2(v) + O → O2(v-1) + O (dove

v indica il numero quantico vibrazionale). Questa rapidità di

diseccitazione si riflette peraltro sulla forma della eedf, che comunque

anche nel nostro modello si rivela non maxwelliana in accordo con le

previsioni e con i dati sperimentali.

3.3.10 Collisioni tra specie pesanti • O- + O → O2 + e- guadagno di 3.66 eV K = 3.00×108 [30] • O- + O2

** → O + O2 + e- perdita di 1.37 eV K = 4.16×108 [30] • O- + O2

** → O3 + e- perdita di 0.32 eV K = 1.81×108 [28] • O- + O2

* → O2- + O perdita di 0.13 eV

K = 6.02×107 [30] • O-

2 + O+2 → 2 O2 guadagno di 11.7 eV

K = 1.2×1011 [30]

48

• O-

2 + O+ → O2 + O guadagno di 10.61 eV K = 1.2×1011 [30] • O- + O+

2 → O2 + O guadagno di 10.68 eV K = 1.15×1011 [28] • O-

2 + O → O- + O2 guadagno di 1.02 eV K = 1.99×108 [28] • O-

2 + O → O3 + e- guadagno di 0.61 eV K = 1.99×108 [28] • O-

2 + O+2 → O2 + 2 O guadagno di 6.58 eV

K = 2.90×1010 [28] • O2 + O2

- → 2 O2 + e- guadagno di 0.46 eV K = 1.20×108 [28] • O+ + O- → 2 O guadagno di 9.59 eV K = 1.20×108 [28] 3.3.11 Scambio di carica • O+

2 + O2 → O2 + O+2 [21]

• O-2 + O2 → O2 + O-

2 [21]

In un processo di scambio di carica si assume che un

elettrone salti dal neutro sullo ione, per diventare uno ione con velocità

uguale a zero nel sistema di riferimento del neutro. Ritornando nel

sistema di riferimento del laboratorio, il nuovo ione possiede la velocità

49

del neutro incidente, e il nuovo neutro acquista la velocità dello ione

incidente.

Nel nostro caso, lo ione prodotto è orientato a caso con

velocità corrispondente all’energia traslazionale termica 3/2kT:

ionemkTv 3= (3.9)

dove k è la costante di Boltzmann e T è la temperatura

traslazionale del gas. Si vedrà quindi come gli ioni prodotto acquistino

maggiore velocità al centro della scarica, dove la temperatura è

maggiore. 3.3.12 Momentum transfer • O- + O2 → O- + O2 • O+ + O2 → O+ + O2

Figura 3.10 Sezioni d’urto per scambio di carica e momentum transfer tra

specie pesanti

10 -3

10 -2

10 -1

10 0

10 1

10 2

10 3

10 4

10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3

σ(cm -16)

Energia(eV)

momentum transfer

O-/O2

scambio di carica

O2+/O2

50

3.3.13 Emissione secondaria

Un altro meccanismo per cui nuovi elettroni possono

entrare nel plasma è quello di emissione secondaria in seguito al

bombardamento del catodo da parte degli ioni positivi. In input viene

dato il coefficiente di emissione secondaria γ, pari a 0.30. Ancora una

volta si estrae un numero casuale: se è minore di γ si verifica una

emissione.

3.4 Verso l’equazione del calore

La figura 3.11 illustra i canali di dissipazione dell’energia in una

scarica RF in ossigeno così come sono considerati nel presente

modello. Il punto di partenza è la densità totale di potenza elettrica

dissipata per unità di volume nel bulk del plasma, mentre il punto di

arrivo è l’energia dissipata sulle pareti e per irradiazione.

Figura 3.11 Canali di dissipazione dell’energia

51

Il primo passo della dissipazione di potenza è

l’accelerazione degli elettroni per mezzo del campo elettrico di scarica.

Il secondo passo sono le collisioni inelastiche degli

elettroni con le molecole di O2, che hanno come conseguenza la

dissociazione o l’eccitazione.

Durante la dissociazione da impatto elettronico, una

frazione dell’energia è rilasciata in energia in eccesso dei frammenti

atomici. Questi atomi di ossigeno “caldi” rilasciano rapidamente la loro

energia cinetica in collisioni elastiche o inelastiche, andando

nuovamente ad eccitare le molecole.

Gli ioni sono poi accelerati verso le pareti. Durante questo

trasporto diffusivo gli ioni possono dar luogo a reazioni esotermiche

con altri ioni o molecole, per cui una frazione dell’entalpia dello stato

fondamentale degli ioni sarà rilasciata nel gas [31].

I canali di perdita di questa energia sono i seguenti:

- dissociazione da impatto elettronico

- reazioni chimiche

- diffusione

- energia trasferita ai fotoni

In particolare è stato stimato che il contributo alla dissipazione di

potenza degli elettroni termici, dell’accelerazione degli elettroni

secondari, e dagli ioni negli sheaths sia pari rispettivamente al 15%, al

5% e all’80% rispettivamente [32].

Si trascura il contributo dell’energia trasferita ai fotoni:la

densità di potenza dell’emissione del plasma tra 300 e 900 nm,

misurata con un radiometro, è pari ad una parte su 105 di tutta la

potenza dissipata [33].

52

Nel presente modello la temperatura delle pareti è fissata ed è

una condizione al contorno per l’equazione di diffusione del calore nel

corpo del plasma, che è uno degli obiettivi di questo lavoro.

Ciononostante si è fatta una stima di quanta sia l’energia persa nelle

collisioni ioni-parete.

Ricordiamo inoltre che ci troviamo nel regime per cui il

meccanismo dominante di dissipazione di potenza elettrica è il

fenomeno di “wave riding” dovuto all’alternarsi di espansione e

contrazione dello sheath.

Questa energia persa si configura come un flusso di calore

q dal plasma agli elettrodi. Esperimenti a diverse pressioni hanno

mostrato che q è costante rispetto al tempo ad alla temperatura delle

pareti [33]. Dal fatto che q è costante deriva che il contributo al flusso

di calore delle molecole con elevata energia cinetica che colpiscono le

pareti è piccolo. Gli elettroni che arrivano agli elettrodi hanno peraltro

energie bassissime. Si può concludere quindi che il flusso di calore dal

plasma agli elettrodi è causato quasi esclusivamente dal

bombardamento di ioni positivi.

Dato che l’energia degli ioni in prossimità delle pareti è

dell’ordine di qualche centinaio di eV, possono essere tralasciate

correzioni riguardanti la perdita di energia degli elettroni secondari, che

hanno energie inferiori a pochi eV.

Da quanto detto, si può concludere che il calcolo della

temperatura nel plasma è un metodo valido per calcolare il flusso di

energia dovuto agli ioni, una quantità che non è facilmente ottenibile

con altri metodi.

53

3.5 L’equazione del calore

3.5.1 Il termine sorgente

L’obiettivo è adesso quello di determinare come la

dissipazione di energia elettrica influenzi la temperatura del plasma;

dobbiamo cioè trovare una equazione da inserire nel programma di

simulazione. L’equazione che fa al caso nostro è ovviamente

l’equazione di diffusione del calore. Essa si presenta in due forme, una

dipendente dal tempo, l’altra indipendente.

L’equazione di cui ci serviremo è quella allo stato

stazionario (vedi discussione alla fine del paragrafo 2.3). Inoltre, per gli

stessi motivi ivi esposti, essa interverrà non certo in ogni ciclo

temporale (dell’ordine di 10-11secondi), ma ogni 50000, parallelamente

all’inserimento della cinetica, con l’avanzamento delle particelle neutre

e il calcolo delle funzioni di distribuzione dell’energia.

La forma di questa equazione è

ddx

kdTdx

+ S = 0 (3.10)

dove k è la conduttività termica e S è il termine sorgente.

La sorgente ha le dimensioni di un flusso di energia, ovvero

Joule/m3sec. Quindi l’energia netta proveniente dalle collisioni va

divisa per il tempo trascorso dall’ultima volta che essa è stata raccolta

(nel nostro caso 50000 volte il time step) e per l’intervallo spaziale

monodimensionale in cui questa energia è stata raccolta. Ricordiamo

che essendo il nostro un modello unidimensionale l’energia raccolta

nelle collisioni è già una energia per metro quadro, essendo le cariche

stesse rappresentate come se fossero dei fogli di carica.

54

Per quanto riguarda le collisioni degli elettroni con O2 e O,

e i processi di scambio di carica e momentum transfer tra ioni, ovvero

quei processi che sono stati descritti tramite le loro sezioni d’urto,

l’energia di attivazione o quella rilasciata è stata rispettivamente

sottratta o aggiunta all’energia cinetica degli elettroni dopo la

collisione. La variazione positiva o negativa dell’energia nel gas è stata

quindi determinata confrontando l’energia cinetica delle particelle in

ingresso ed in uscita dalla subroutine che tratta queste collisioni.

Il bilancio di energia per gli altri processi è stato calcolato

aggiungendo o sottraendo energia ogniqualvolta che si verificava ogni

singolo processo.

Il termine sorgente non è però costituito esclusivamente

dai contributi provenienti dalle collisioni; deve infatti contenere anche

l’altro canale di dissipazione dell’energia: l’irraggiamento.

Studi sperimentali [33] hanno infatti mostrato come esso

giochi un ruolo cruciale nella dispersione del calore.

Il termine usato per tener conto dell’irraggiamento è la formula

dell’irraggiamento da corpo nero, per cui il calore irraggiato da una

superficie di area unitaria è dato da

( )40

4 TT −σ (3.11)

dove T0 è la temperatura del reattore e σ è la costante di Stefan-

Boltzmann, pari a 5.67 × 10-8 W m-2 K-4.

Nel considerare gli ultimi due contributi al termine

sorgente, ovvero i processi caratterizzati dalla costante di velocità e

l’irraggiamento, si è dovuto tener conto dell’effettivo volume

tridimensionale della scarica. Le proprietà del plasma che andiamo ad

55

investigare con un modello unidimensionale sono grandezze intensive

(densità, energia e velocità di ogni singola particella), quindi

l’assunzione ovviamente irrealistica che la scarica misuri un metro

quadro in due piani perpendicolari alla direzione lungo la quale avviene

la simulazione non inficia il valore dei risultati ottenuti.

Quando però si va ad inserire il termine dovuto

all’irraggiamento, già l’analisi dimensionale mostra che, per ottenere

un flusso di energia, espresso in W/m3, l’equazione 3.11 va moltiplicata

per l’inverso di una distanza, ovvero per un rapporto tra la superficie

che irradia e il volume in essa contenuto.

Allo stesso modo, nel considerare l’energia rilasciata o

assorbita da una reazione in una regione del nostro dominio

unidimensionale, dobbiamo tener conto delle effettive dimensioni

trasverse in m2 della scarica.

Per questi dati, come del resto per la temperatura da fissare

alle pareti, si è fatto riferimento ai reattori in [32] e [33], che sono del

tipo e della misura più frequentemente utilizzati per la produzione del

tipo di plasmi qui considerato e le cui misure nella direzione del

diametro degli elettrodi sono state impiegate in modelli bidimensionali

[34].

I reattori di cui sopra hanno elettrodi circolari aventi

diametro di 30 cm spaziati di 2 cm (ed è proprio lungo la direzione che

collega i due elettrodi che è sviluppata la nostra simulazione).

56

3.5.2 Conduttività

La conduttività termica dipende dalla temperatura secondo la

relazione [35] :

k = k0

TT0

nJ

m ⋅ s ⋅ °K (3.12)

con k0 = 2.4476 × 10-2 J/m⋅s ⋅°K e n = 0.87.

Dato che il grado di ionizzazione e di dissociazione in questo

plasma è basso, si può assumere che le variazioni relativamente piccole

della composizione del gas mentre esso si trova nello stato di plasma

abbiano effetti trascurabili sulle proprietà di conduzione di calore del

gas.

57

CAPITOLO 4 DISCRETIZZAZIONE ED

IMPLEMENTAZIONE DELL’EQUAZIONE DEL CALORE

4.1 Forma generale dell’equazione da discretizzare

La conduzione di calore stazionaria unidimensionale è governata

dalla equazione

ddx

kdTdx

+ S = 0 (4.1)

dove k è la conduttività termica, T è la temperatura e S è il tasso

di generazione di calore per unità di volume .

Per discretizzare l’equazione si impiega la griglia di figura 4.1.

Figura 4.1 Agglomerato di punti griglia per il problema unidimensionale

Focalizziamo l’attenzione sul punto P, che ha come vicini i punti

E e W, rispettivamente nella direzione delle x crescenti e decrescenti.

Le linee tratteggiate indicano le facce del volume di controllo ,

denotati con e e w.

58

Per il problema unidimensionale che stiamo considerando si

assume spessore unitario nelle direzioni y e z. Quindi il volume di

controllo ha volume ∆x×1×1.

Integrando l’equazione 4.1 sul volume di controllo otteniamo

kdTdx

e− k

dTdx

w+ Sdx

w

e

∫ = 0 (4.2)

Allo stato attuale si assegna ad ogni punto della griglia un valore

di temperatura che diventa il valore assunto in tutto il volume di

controllo.

Questa assegnazione restituisce un profilo a scala come in figura

4.2(a), per il quale la pendenza dT/dx non è definita ai lati del volume

di controllo . Un profilo che non soffre di questa limitazione è quello

lineare a tratti della figura 4.2(b) dove si è usata una interpolazione

lineare tra i punti di griglia .

Tuttavia una griglia abbastanza fitta da permettere stabilità in

problemi ben più complessi quale è la nostra può garantire la quasi

linearità del profilo di temperatura e possiamo quindi calcolare le

derivate nell’ equazione 4.2 ottenendo

(4.3)

dove S è il valore medio di S sul volume di controllo .

E’ utile riscrivere l’ equazione discretizzata nella forma : aPTP = aETE + aWTW +b (4.4)

0)(

)()(

)( =∆+−−− xSx

TTkx

TTk

w

WPw

e

PEe

δδ

59

dove

,)(

,)(

w

ww

e

eE

xka

xka

δ

δ

=

=

(4.5)

xSbaaa WEP

∆=

+= ,

Figura 4.1 Due semplici tipologie di profilo. (a) a scala; (b) lineare a tratti 4.2 Regole di base

A. Consistenza ai lati del volume di controllo

Il flusso attraverso una faccia comune a due volumi di controllo

deve essere rappresentato dalla stessa espressione nelle equazioni

discretizzate per i due volumi .

60

B. Coefficienti positivi

Il valore di una variabile dipendente in corrispondenza di un

punto di griglia è influenzato dai valori ai punti confinanti tramite il

processo di diffusione .

Ovviamente un aumento del valore ad un punto , in assenza di

altri cambiamenti , deve portare ad un aumento del valore al punto

confinante , ovvero i coefficienti a devono avere lo stesso segno , che

assumiamo positivo .

C. Somma dei coefficienti confinanti Se l’ equazione 4.1 è soddisfatta da T essa deve essere anche

essere soddisfatta da T+c , con c costante arbitraria . Questa proprietà

dell’ equazione differenziale deve essere riflessa nell’ equazione

discretizzata , ovvero l’equazione 4.1 deve rimanere valida quando TP ,

TE e TW sono aumentati di una costante , da cui segue

aP=aE+aW (4.6) La nostra equazione ovviamente soddisfa questa regola , la quale implica che TP è una media pesata dei valori confinanti .

Passiamo ora a vedere in dettaglio altri aspetti del problema della conduzione del calore in una dimensione .

D. Spaziatura della griglia

Per i punti di griglia mostrati in figura 4.1 non è necessario che le

distanze (δx)e e (δx)w siano uguali , anzi a volte l’uso di una griglia

nonuniforme è indicato per risparmiare tempo di macchina . In generale

61

si otterrà una soluzione accurata solo quando la griglia è

sufficientemente fine , ovvero tanto piu fine quanto più è rapida la

variazione di T con x .

E. Conduttività all’ interfaccia

Nell’ equazione 4.3 ke e kw rappresentano i valori della

conduttività sulle facce e e w sulle facce del volume di controllo .

Quando k è una funzione di x spesso conosciamo i valori di k in

corrispndenza dei punti di griglia W , P , E e così via .

Serve quindi un metodo per calcolare la conduttività all’

interfaccia , diciamo ke , in termini dei valori sui punti della griglia .

Ovviamente questa discussione è irrilevante quando la

conduttività è uniforme.

Disuniformità nella conduttività possono sorgere dalla

disomogeneità del materiale, o, come nel nostro caso, dalla

distribuzione di temperatura.

La procedura più immediata per ottenere la conduttività ke all’

interfaccia è assumere una variazione lineare di k tra i punti P ed E :

Eepee kfkfk )1( −+= (4.7)

dove il fattore di interpolazione fe è un rapporto definito in

termini delle distanze mostrate in figura 4.3 :

e

ee x

xf)()(

δδ += (4.8)

62

Se l’ interfaccia e fosse esattamente a metà tra i punti di griglia ,

fe sarebbe 0.5 e ke sarebbe la media aritmetica di kP e kE .

Questo semplice approccio porta ad implicazioni scorrette e non

può affrontare correttamente gli improvvisi cambi di conduttività che si

verificano nei materiali composti .

Figura 4.3 Distanze associate all’interfaccia e

Nel ricercare un’ alternativa ricordiamo che non è il valore locale

della conduttività all’ interfaccia e che ci interessa direttamente . Il

nostro obiettivo principale è di ottenere una buona rappresentazione del

flusso di calore qe all’ interfaccia tramite

e

EPee x

TTkq)(

)(δ

−= (4.9)

che è stata usata per ottenere l’ equazione discretizzata .

L’ espressione cercata per ke è quella che porta ad un qe

“corretto” .

Immaginiamo che il volume di controllo che circonda il punto P

sia riempito con un materiale di conduttività uniforme kP e quello

intorno ad E con un materiale di conduttività kE . Per la “fetta” di

63

materiale compresa tra i punti P ed E , una analisi stazionaria

unidimensionale (senza sorgenti ) porta a

EePe

EPe kxkx

TTq/)(/)( +− +

−=

δδ (4.10)

La combinazione delle equazioni 4.8 e 4.10 dà

ke =1 − fe

kP

+fe

kE

−1

(4.11)

Quando l’ interfaccia è posta a metà tra P ed E si ottiene fe = 0.5 ,

ovvero

EP

EPe kk

kkk+

= 2 (4.12)

L’ equazione 4.12 mostra che ke è la media armonica di kP e kE

piuttosto che quella aritmetica che darebbe la 4.7 quando fe = 0.5 .

L’ uso dell eq. 4.11 nelle definizioni dei coefficienti 4.5 porta

alla seguente espressione per aE :

aE = (δx)e −

kP

+ (δx)e +

kE

−1

(4.13)

Una espressione analoga sarà scritta per aW . Ovviamente aE

rappresenta la conduttanza del materiale tra i punti P ed E .

64

La validità di questa formulazione è confermata dai sue seguenti

casi limite :

1. Sia kE→0 . Quindi dall’equazione 4.11 ,

ke→ 0. (4.14)

Questo implica che il flusso di calore attraverso una faccia di un

isolante diventa zero, come dovrebbe, laddove invece la formulazione

basata sulla media aritmetica avrebbe dato un flusso non nullo .

Sia kP>>kE . Quindi :

ke →kE

fe 4.15)

Questo risultato ha due implicazioni; la prima più immediata, la

seconda più sottile. L’ equazione 4.15 indica che la conduttività

all’interfaccia ke non dipende affatto da kP . Ciò è prevedibile dato che

il materiale ad alta conduttività attorno al punto P offre una resistenza

trascurabile rispetto a quella offerta dal materiale attorno ad E la

formulazione basata sulla media aritmetica avrebbe trattenuto l’ effetto

di kP su kE ) .

L’altra implicazione è che ke non è uguale a kE , bensì il loro

rapporto è 1/fe .

Il nostro scopo è di ottenere un valore corretto di qe tramite l’ eq.

4.9 . L’uso dell’eq. 4.15 dà

65

+

−=

e

EPEe x

TTkq)(

)(δ (4.16)

Quando kP>>kE , la temperatura TP prevarrà fino all’ interfaccia

e , e la caduta di temperatura TP – TE avrà luogo lungo (δ x)e+ . Perciò il

valore corretto del flusso sarà quello dato dall’ eq. 4.16 . In altre parole

il fattore fe nell’eq. 4.15 compensa l’ uso della distanza nominale (δ x)e

nell’ eq. 4.9 .

La considerazione di questi due casi limite mostra che la

formulazione può far fronte a repentini cambi di conduttività senza

richiedere una griglia eccessivamente fine in prossimità della

discontinuità .

La formula 4.11 per la conduttività all’ interfaccia è basata sulla

situazione unidimensionale stazionaria senza termine sorgente in cui la

conduttività varia secondo un profilo a gradini da un volume di

controllo ad un altro .

Ma anche in situazioni come la nostra con sorgente non nulla e

con variazione continua di conduttività questa formula approssima

molto meglio il reale comportamento [36]. 4.3 Nonlinearità

L’ equazione di discretizzazione 4.4 è un’equazione algebrica

lineare, risolvibile con i metodi appositi per tali equazioni . Tuttavia nel

problema della conduzione del calore si incontrano spesso situazioni in

cui si perde la linearità.

Questo accade ad esempio quando k dipende da T , come nel

presente caso, o quando la sorgente S è una funzione nonlineare di T .

66

In questi casi i coefficienti nella equazione di discretizzazione

dipenderanno essi stessi da T . In queste situazioni la procedura da

adottare è quella dell’ iterazione , che consiste dei seguenti passi :

1. Si parte da una stima o da un valore a caso per i valori di T in

corrispondenza di tutti i punti di griglia .

2. Da questi valori di T si calcolano dei valori di prova dei coefficienti

nell’ equazione di discretizzazione .

3. Si risolve l’ insieme , nominalmente lineare , di equazioni

algebriche per ottenere nuovi valori di T .

4. Con questi nuovi valori di T come stime migliori si ritorna al passo

2 e si ripete il processo fino a che successive iterazioni cessano di

ptodurre variazioni significative nei valori di T.

Lo stato finale senza variazioni è chiamato convergenza delle

iterazioni. *

La soluzione convergente è proprio la soluzione corretta delle

equazioni nonlineari , anche se ci si è arrivati risolvendo equazioni

lineari .

E’ in ogni caso possibile che successive iterazioni non

convergano mai ad una soluzione . I valori di T possono allontanarsi

sempre più o oscillare con ampiezza crescente . Questo processo è

detto divergenza , e un buon metodo numerico deve minimizzare le

possibilità che esso si verifichi. Di certo il rispetto delle tre regole

derivate in precedenza è un buon passo verso la convergenza .

* Alle volte il termine convergenza è usato per indicare il processo in seguito al quale successive regolazioni della griglia portano la soluzione numerica più vicina alla soluzione esatta . Questo aspetto sarà chiamato in questo lavoro ‘accuratezza’ della soluzione numerica , riservando la parola convergenza allo stato finale delle iterazioni qui descritte.

67

4.4 Condizioni al contorno

Assumiamo che , nel nostro problema unidimensionale sia stata

scelta la serie di punti di griglia mostrata in figura 4.4. C’è un punto di

griglia su ciascuno dei due confini . Gli altri punti sono detti punti

interni , attorno a ciascuno dei quali è mostrato un volume di controllo .

Per ognuno di questi volumi di controllo può essere scritta una

equazione discretizzata come la 4.4. Se consideriamo l’ equazione 4.4

come una equazione per TP, abbiamo le equazioni necessarie per tutte

le temperature incognite in corrispondenza dei punti di griglia interni .

Due di queste però contengono le temperature dei punti di griglia sui

contorni . Attraverso il trattamento di queste temperature al contorno si

introduconi le condizioni al contorno nello schema della soluzione

numerica.

Poniamo l’ attenzione sul punto B sul confine sinistro , che è

adiacente al primo punto interno I , come mostrato in figura.

Figura 4.4 Volumi di controllo per i punti interni e per il contorno

In genere si incontrano tre tipi di condizioni al contorno nei

problemi di conduzione del calore . Essi sono :

68

1. Data temperatura sul contorno

2. Dato flusso di calore in corrispondenza del contorno

3. Flusso di calore in corrispondenza del contorno specificato da un

coefficiente di trasmissione del calore e dalla temperatura del fluido

circostante

Se è data la temperatura del contorno ( ovvero si conosce il

valore di TB) , non ci sono particolari difficoltà e non sono richieste

ulteriori equazioni .

Quando non è specificata la temperatura al contorno bisogna

costruire una equazione supplementare per TB integrando l’equazione

differenziale sul ‘mezzo’ volume di controllo adiacente al contorno

(fig.4.5) .

Figura 4.5 Mezzo volume di controllo adiacente al contorno

Si ottiene

0=∆+− xSqq iB (4.17)

dove il flusso di calore q sta per –k dT/dx .

Il flusso all’ interfaccia qi può essere scritto come si è fatto

nell’equazione 4.9 . Il risultato è :

69

0)(

)(=∆+

−− xS

xTTk

qi

IBiB δ (4.18)

Una ulteriore implementazione di questa equazione dipende da

quanto è dato circa il flusso di calore al contorno qB . Se il suo valore è

dato, l’ equazione per TB diventa

bTaTa IIBB += (4.19)

dove

IB

B

i

iI

aaqxSb

xka

=+∆=

=

,

,)(δ

(4.20)

Se il flusso di calore qB è specificato in termini di un coefficiente

di trasmissione di calore h e una temperatura del fluido circostante Tf

tale che

),( BfB TThq −= (4.21)

l’ equazione per TB diventa

,bTaTa IIBB += (4.22)

70

dove

haahTxSb

xka

IB

f

i

iI

+=

+∆=

=

,

,)(δ

(4.23)

71

4.5 Risoluzione delle equazioni algebriche

I metodi per risolvere le equazioni algebriche si dividono in

diretti o iterativi .

Il più semplice dei metodi diretti è il metodo TDMA

(TriDiagonal Matrix Algorithm), cosiddetto perchè scrivendo la

matrice dei coefficienti delle equazioni , tutti i valori diversi da zero si

allineano su tre diagonali della matrice *.

Questo metodo usa una relazione di ricorrenza che permette di

risolvere le equazioni dall’ultima alla prima .

I metodi diretti sono accettabili nella risoluzione di equazioni

lineari , ma quando le equazioni vanno risolte ripetutamente con

coefficienti aggiornati essi portano ad un elevato spreco di tempo e di

memoria della macchina . Inoltre diventano eccessivamente complicati

e dispendiosi passando a più dimensioni .

L’alternativa è data dai metodi iterativi , che richiedono molta

meno memoria e sono particolarmente indicati nel trattare la

nonlinearità . Infatti in un problema non lineare non è necessario nè

saggio portare la soluzione delle equazioni algebriche alla convergenza

finale per un insieme di valori fissi per i coefficienti . Una volta scelto

un primo insieme di questi valori , bastano poche iterazioni

dell’algoritmo risolvente prima che venga effettuato l’aggiornamento

dei coefficienti . In generale ci vuole un certo equilibrio tra lo sforzo

richiesto per calcolare i coefficienti e quello richiesto per risolvere le

equazioni . Non è certo sensato spendere troppo nel risolvere equazioni

basate su coefficienti che sono solo un tentativo .

La base di tutti i metodi iterativi è il metodo di Gauss-Seidel in

cui i valori della variabile sono calcolati visitando in un dato ordine * Ciò avviene perchè ci si riferisce ai valori nel punto in interesse e nei suoi due punti adiacenti

72

tutti i punti della griglia [37]. Nella memoria del calcolatore rimane

solo un insieme di valori di T per volta . All’inizio questi rappresentano

la stima iniziale o i valori che vengono dalla precedente iterazione .

Man mano che si arriva su ogni punto della griglia , il corrispondente

valore di T nella memoria della macchina è variato in questo modo :

- se l’equazione discretizzata è scritta

bTaTaTa WWEEPP ++= (4.24)

il valore di T al punto P sarà dato da :

P

WWEEP a

bTaTaT ++=**

(4.25)

dove T*E,B sono i valori in corrispondenza dei punti limitrofi

immagazzinati nella memoria . Per punti vicini che sono già stati

visitati nell’iterazione corrente , T* è il valore appena calcolato , per

punti vicini che devono ancora essere visitati , T* è il valore che viene

dalla iterazione precedente . Quando tutti i punti sono stati visitati in

questo modo , è stata completata una iterazione del metodo di Gauss-

Seidel .

Per illustrare il metodo , consideriamo due semplici esempi .

Equazioni :

.12.04.0

12

21

+=+=

TTTT

(4.26)

73

Soluzione :

Iterazione n. 0 1 2 3 4 5 ... ∞

T1 0 0.2 0.68 0.872 0.949 0.980 ... 1.0

T2 0 1.2 1.68 1.872 1.949 1.980 ... 2.0

Si può vedere che , partendo con un valore arbitrario ci si è potuti

avvicinare all’esatta soluzione delle equazioni . Un aspetto interessante

dei metodi iterativi è che la precisione dei calcoli non deve essere

necessariamente elevata nei passaggi intermedi . I calcoli approssimati ,

e anche gli errori , tendono a essere eliminati , dato che i valori

intermedi sono usati semplicemente come stime per la successiva

iterazione . Ricaviamo ulteriori informazioni dal seguente esempio .

Equazioni :

5.05.2,1

12

21

−=−=TT

TT (4.27)

Soluzione :

Iterazione n. 0 1 2 3 4

T1 0 -1 -4 -11.5 -30.25

T2 0 -3 -10.5 -29.25 -76.13

Non sembra che ci siano molte speranze . Qui il processo di

iterazione è andato incontro a divergenza .

74

Dopotutto le equazioni 4.27 altro non sono che versioni

riarrangiate delle equazioni 4.26, per le quali si era ottenuta la

convergenza , che quindi non è sempre garantita dal metodo di Gauss-

Seidel .

Scarborough ha formulato un criterio che, se soddisfatto, assicura

la convergenza .

Criterio di Scarborough . Una condizione sufficiente per la

convergenza del metodo di Gauss-Seidel è

av∑

aP≤1 (4.28)

per tutte le equazioni tranne almeno una per cui l’uguaglianza

deve valere in senso stretto , con av che rappresenta il coefficiente di T

per i punti vicini , in questo caso E e W .

Commenti :

1. Il criterio è una condizione sufficiente , non necessaria . Ovvero può

essere talvolta violato ottenendo comunque la convergenza.

2. Alcune delle regole base, che sono state motivate da considerazioni

fisiche, soddisfano le richieste del criterio di Scarborough . Ad

esempio se uno dei coefficienti fosse negativo, aP , che nel nostro

caso è uguale a Σav , avrebbe un modulo minore di Σ|av|, portando

ad una violazione del criterio.

3. Quando, come nel nostro caso, aP = Σav, e tutti i coefficienti sono

positivi, otteniamo, per tutte le equazioni, Σ|av|/|aP| = 1 . In che caso,

75

quindi, il rapporto è minore dell’unità? La risposta è nelle

condizioni al contorno. Affinchè il problema abbia una determinata

soluzione, la temperatura deve essere specificata per almeno un

punto sul contorno. L’equazione discretizzata in cui questo punto

appare come uno dei punti vicini implica Σ|av|/|aP|<1. Questo perchè

Σ|av| va calcolato come somma dei coefficienti dei vicini incogniti;

d’altro canto aP è la somma di tutti i coefficienti dei punti confinanti

incluso quello del punto sul contorno. 4.6 Sovrarilassamento e sottorilassamento

Nella soluzione iterativa delle equazioni algebriche o nello

schema iterativo usato per trattare la nonlinearità è spesso conveniente ,

anche alla luce di quanto si è detto sull’equilibrio tra sforzi per

calcolare i coefficienti e sforzi per calcolare i valori di T , accelerare o

rallentare i cambiamenti dei valori della variabile dipendente da una

iterazione all’altra . Questo processo è detto sovrarilassamento o

sottorilassamento a seconda che le variazioni siano accelerate o

rallentate . Il sovrarilassamento è spesso usato in coppia con il metodo

di Gauss-Seidel prendendo il nome di SOR (Successive Over-

Relaxation) . Il sottorilassamento è uno strumento molto utile per

problemi non lineari [38] .

Vediamo come viene modificata la nostra equazione . Si parte da

P

WWEEP a

bTaTaT ++= (4.29)

76

Sommando e sottraendo il valore TP* (che come si è detto

rappresenta il valore di TP proveniente dalla precedente iterazione) dal

secondo membro si ottiene :

TP = TP* + aETE + aWTW + b

aP

− TP*

(4.30)

dove il contenuto della parentesi è proprio la variazione di TP

prodotta dall’iterazione corrente . Questa variazione può essere

modificata introducendo un fattore di rilassamento α , ovvero

TP = TP* +α aETE + aWTW + b

aP− TP

*

(4.31a)

cioè

*)1( PP

WWEEPP TabTaTaTa

αα

α−+++= (4.31b)

Notiamo in primo luogo che quando le iterazioni convergono,

ovvero TP diventa uguale a TP*, la 4.31a implica che i valori di T

soddisfano l’equazione originale 4.29 .

Questa è una proprietà che ogni schema di rilassamento deve

possedere; la soluzione finale deve soddisfare l’equazione di

discretizzazione comunque si sia arrivati alla convergenza .

77

Quando il fattore di rilassamento α è compreso tra 0 e 1 il suo

effetto è il sottorilassamento , ovvero i valori di TP sono più vicini a TP*

. Quando α è maggiore di 1 si ha sovrarilassamento . Non vi sono

regole generali per la scelta del migliore valore di α . In questo caso si

è scelto il sottorilassamento trattandosi di un caso di nonlinearità , e si è

posto α uguale a 0.8 . 4.7 Forma finale dell’ equazione discretizzata

Partiamo dall’equazione 4.4 . Combinando la 4.5 con la 4.12 si

ottiene :

aP = 2∆x

k(TP)k(TE )k(TP ) + k(TE )

+ k(TP )k(TW )k(TP ) + k(TW )

aW =2

∆xk(TP)k(TW )

k(TP ) + k(TW)

aE =2

∆xk(TP )k(TE )

k(TP ) + k(TE )

(4.32)

e la forma finale dell’equazione discretizzata è quindi :

78

2∆x

k(TP)k(TE )k(TP) + k(TE )

+ k(TP)k(TW )k(TP ) + k(TW )

TP

α=

= 2∆x

k(TP )k(TW)k(TP) + k(TW )

TW +

+ 2∆x

k(TP )k(TE)k(TP ) + k(TE )

TE + SP∆x +

+ (1− α)α

2∆x

k(TP )k(TE)k(TP ) + k(TE )

+ k(TP )k(TW)k(TP ) + k(TW )

TP

*

(4.33)

79

CAPITOLO 5 RISULTATI

La simulazione è stata condotta utilizzando due versioni

differenti del programma: una in cui la temperatura del gas era fissata

pari a 300 ºK e una in cui è stata inserita l’equazione del calore

discretizzata, che conduce ad un profilo non costante di temperatura.

Parte di questo capitolo sarà dedicata alle differenze riscontrate

confrontando i risultati ottenuti, laddove queste ultime sono più

significative.

Inoltre saranno discussi aspetti relativi alle scariche a

radiofrequenza in ossigeno evidenziati da entrambi i programmi.

Sono stati simulati circa 1900 cicli RF, sufficienti a raggiungere

la stabilità, condizione quest’ultima assicurata dal fatto che il numero

di particelle simulate non varia che in minima parte all’aumentare dei

passi di calcolo.

Si passano adesso in rassegna i diversi parametri di plasma.

5.1 Temperatura

La temperatura che si ottiene quando l’energia depositata o

assorbita a seguito delle varie reazioni e collisioni viene inserita come

sorgente nell’equazione allo stato stazionario per la diffusione del

calore mostra il profilo riportato in figura 5.1. L’analisi di Fourier

assicura che la forma del profilo di temperatura è la stessa del profilo

del termine sorgente. Si evince che la maggior parte dell’energia viene

depositata al centro del glow, e questo viene spiegato da una analisi

comparata dei profili di densità e di temperatura delle varie specie che

mostra che per quanto gli ioni abbiano maggiore energia in

80

corrispondenza degli sheath, la loro densità, così come quella degli

elettroni, è molto maggiore al centro. A tale proposito si veda in

seguito il paragrafo dedicato all’energia depositata dalle diverse

reazioni.

Con le pareti del reattore fisse a 300 ºK , i risultati ottenuti

sono in perfetto accordo con quelli misurati sperimentalmente [33].

Figura 5.1 Temperatura del gas 5.2 Potenziale

L’andamento del potenziale medio è mostrato in figura

5.2. Esso è stato ottenuto mediando tutti i valori del potenziale calcolati

dal PIC ad ogni ciclo. Il profilo è simmetrico in accordo con la

geometria del reattore qui simulato.

2,5 102

3,0 102

3,5 102

4,0 102

4,5 102

5,0 102

0 100

5 10-3

1 10-2

2 10-2

2 10-2

Temperatura del gas (Kelvin)

posizione (m)

81

In corrispondenza degli elettrodi il potenziale medio è

chiaramente nullo. A causa della maggiore velocità degli elettroni, che

diffondono più rapidamente verso le pareti, il plasma al centro del glow

si ritrova con una prevalenza di carica positiva e quindi con un

potenziale medio più elevato rispetto agli elettrodi. Peraltro,

considerando la velocità con cui le cariche positive e negative arrivano

agli elettrodi, si ha un flusso di carica netto pari a zero.

Il potenziale di plasma Vp si assesta ad un valore pari a

circa 0.4 VRF, in accordo con la teoria[39].

Figura 5.2 Potenziale medio in funzione dalla posizione

0

20

40

60

80

100

0 100

1 10-2

2 10-2

Potenziale

potenziale medio (V)

posizione (m)

82

5.3 Densità

In figura 5.3 è riportata la densità degli elettroni.

L’avvallamento al centro è dovuto all’inerzia degli ioni negativi al

centro del plasma, per cui gli elettroni “vedono” una maggiore carica

negativa al centro e non vi si posizionano per assicurare la neutralità

[40].

Il numero di elettroni si riduce drasticamente in

corrispondenza degli elettrodi a causa della diffusione ambipolare e

della formazione dello sheath.

Per quanto riguarda lo studio a due diverse temperature, si

vede che il numero di elettroni al centro della scarica è circa del 20%

maggiore nel caso di temperatura variabile, quindi laddove il profilo di

temperatura raggiunge il massimo, mentre la densità si mantiene

costante nelle regioni di diffusione. Ci sono due processi principali che

conseguono all’aumento di temperatura, e che producono effetti

contrastanti: da un lato ad una temperatura maggiore corrisponde una

maggiore diffusione, e questo porterebbe ad una diminuzione del

numero di elettroni nel centro del plasma; d’altra parte con la

temperatura cresce anche il numero di collisioni, che portano ad un

intrappolamento degli elettroni nella regione ambipolare. Nel caso in

analisi prevale quest’ultimo effetto.

83

Figura 5.3 Densità elettronica a due temperature

Le figure da 5.4 a 5.7 mostrano la densità degli ioni.

Anche in questo caso la densità è maggiore al centro, in quanto gli ioni

vengono prodotti in primo luogo per impatto elettronico. L’effetto della

temperatura è pressochè ininfluente sugli ioni negativi, mentre porta ad

una diversa densità gli ioni positivi, in particolare gli O+ nei pressi

degli sheath. Questi fenomeni trovano spiegazione nelle seguenti

considerazioni:

- gli ioni O2- sono prodotti per dissociazione dell’O3 o per scambio di

carica tra O- e O2: l’ozono è una specie relativamente pesante e

quindi poco mobile; gli ioni negativi, dopo essere stati prodotti,

vengono sospinti verso il centro del plasma dal campo ambipolare .

Da questo segue che gli ioni negativi vengono poco influenzati dalle

variazioni di temperatura.

- gli ioni positivi sono invece maggiormente accelerati verso gli

elettrodi, in seguito agli effetti cooperativi del campo ambipolare e

1013

1014

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

densità elettroni

T variabileT costante

elettroni (m

-3 )

posizione (m)

84

della maggiore velocità termica. Gli ioni O+, più leggeri degli O2+ si

perdono con maggiore velocità.

Figura 5.4 Densità degli ioni O2

+

1014

1015

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

densità ioni O2

+

T costante

T variabile

ioni O

2

+ (m-

3 )

posizione (m)

1011

1012

1013

1014

1015

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

densità ioni O-

T costanteT variabile

ioni O

- (m-

3 )

posizione (m)

85

Figura 5.5 Densità degli ioni O+

Figura 5.6 Densità degli ioni O2

-

108

109

1010

1011

1012

1013

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

densità ioni O+


ioni O

+ (m-

3 )

posizione (m)

1011

1012

1013

1014

1015

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

densità ioni O2

-

T costanteT variabile

ioni O

2

- (m-

3 )

posizione (m)

86

Figura 5.7 Densità degli ioni O+ 5.4 Energia media

Le figure da 5.8 a 5.10 riportano l’energia media delle tre

specie più rappresentate nella presente simulazione: elettroni, ioni O2+ e

ioni O-.

Osservando l'energia cinetica media delle particelle,

notiamo ancora come la fenomenologia sperimentale viene riprodotta

da principi primi: gli elettroni hanno energia media dell'ordine dell'eV

nel centro della scarica, dove gli ioni sono invece ad energie termiche.

Questo accade perché gli elettroni, molto leggeri, traggono energia dal

forte campo elettrico nella regione dello sheath e la trasportano per

diffusine nel centro della scarica. Si ha quindi un regime di

riscaldamento non locale, a differenza di quanto avviene nelle scariche

108

109

1010

1011

1012

1013

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

densità ioni O+


ioni O

+ (m-

3 )

posizione (m)

87

ad elevata pressione (in cui la diffusione è poco importante) dove

l'energia media degli elettroni è determinata dal campo elettrico locale.

Gli ioni, molto più pesanti, si muovono più lentamente e sono

caratterizzati da più bassi coefficienti di diffusione per cui il trasporto

di energia nel bulk è molto meno importante. Inoltre gli ioni

trasferiscono molta energia negli urti elastici con le altre particelle

pesanti, a differenza di quanto accade per gli elettroni. Questo spiega

perché nel bulk gli ioni abbiano energie termiche.

Nelle regioni periferiche si ha una moderata diminuzione

dell'energia elettronica ed un vistosissimo aumento dell'energia degli

ioni positivi, che vengono accelerati verso gli elettrodi dal campo di

sheath ed ivi raggiungono la massima energia.

Figura 5.8 Energia media degli elettroni

Per quanto detto in precedenza, gli elettroni acquistano maggiore

energia cinetica lungo tutta la scarica a causa della aumentata

temperatura.

0,0 100

5,0 100

1,0 101

1,5 101

2,0 101

2,5 101

3,0 101

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

energia degli elettroni


energia (eV)

posizione (m)

88

Figura 5.9 Energia media degli ioni O2+

Figura 5.10 Energia media degli ioni O-

Gli ioni O2+ sono poco influenzati dall’aumento di temperatura,

al contrario di quanto avviene per gli ioni O- nei pressi degli sheath.

Questo aumento di energia cinetica può indicare un abbassamento della

10-2

10-1

100

101

102

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

energia ioni O2

+


energia (eV)

posizione (m)

10-1

100

101

102

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

energia degli ioni O-


energia (eV)

posizione (m)

89

buca di diffusione ambipolare, al quale bordo peraltro gli ioni arrivano

con velocità maggiore in quanto trovano un gas maggiormente rarefatto

in seguito alla temperatura più elevata.

5.5 Velocità di deriva

I grafici della velocità di deriva, da figura 5.11 a figura 5.15,

mostrano chiaramente che la scarica a piatti paralleli in esame opera in

regime di diffusione ambipolare, in cui il movimento degli elettroni è

condizionato da quello degli ioni che, essendo più pesanti, diffondono

più lentamente. Infatti ioni ed elettroni hanno ovunque velocità di

deriva con lo stesso segno, contrariamente a quanto avverrebbe se

avessimo invece un trasporto convettivo dominato dall'energia locale.

Figura 5.11 Velocità di deriva degli elettroni

Guardando il grafico della velocità degli elettroni, si nota che

questi ultimi non sono accelerati dall’aumento di temperatura, in

-6,0 105

-4,0 105

-2,0 105

0,0 100

2,0 105

4,0 105

6,0 105

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

velocità degli elettroni


velocità (m/s)

posizione (m)

90

accordo con quanto si era detto a proposito della densità, ovvero che la

diffusione è fortemente limitata dalle collisioni.

Figura 5.12 Velocità di deriva degli ioni O2+ lungo tutta la scarica

Figura 5.13 Velocità di deriva degli ioni O2+ in prossimità del centro della

scarica

-4,0 104

-3,0 104

-2,0 104

-1,0 104

0,0 100

1,0 104

2,0 104

3,0 104

4,0 104

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

velocità ioni O2

+


velocità (m/s)

posizione (m)

-4,0 102

-3,0 102

-2,0 102

-1,0 102

0,0 100

1,0 102

2,0 102

3,0 102

4,0 102

6,0 10-3

8,0 10-3

1,0 10-2

1,2 10-2

1,4 10-2

velocità ioni O2

+


velocità (m/s)

posizione (m)

91

Dalla figura 5.12 si evince che gli ioni O2+, pur entrando

nello sheath con una energia maggiore quando la temperatura è più

elevata, non raggiungono velocità maggiori rispetto al caso di

temperatura costante, ma sono addirittura rallentati dalle collisioni.

La figura 5.13 mostra l’andamento della velocità

all’interno del bulk. È interessante notare come nel caso di temperatura

variabile siano evidenziati dei cambi di segno della derivata. Questo

fenomeno evidenzia l’esistenza di un doppio strato di carica dovuto al

richiamo di cariche negative conseguente all’espansione e

restringimento dello sheath. Gli ioni negativi pesanti vengono

richiamati più lentamente degli elettroni, andando così a formare un

secondo sheath.

Per quanto riguarda gli ioni O-, essi vengono accelerati

molto di più all’aumentare della temperatura, come mostra la figura

5.14. In questo caso, oltre alla minore massa rispetto agli ioni O2+, le

collisioni e le ricombinazioni sono molto meno presenti.

Figura 5.14 Velocità di deriva degli ioni O-

-2,0 106

-1,0 106

0,0 100

1,0 106

2,0 106

3,0 106

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

velocità ioni O-


velocità (m/s)

posizione (m)

92

Figura 5.15 Velocità di deriva degli ioni O- a temperatura costante

Per maggiore chiarezza la figura 5.15 mostra l’andamento della

velocità degli ioni O- a temperatura costante, che non era ben

evidenziato nella figura precedente a causa della diversa scala.

-3,0 104

-2,0 104

-1,0 104

0,0 100

1,0 104

2,0 104

3,0 104

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

velocità ioni O-

T costante

velocità (m/s)

posizione (m)

93

5.5 Densità, energia e velocità in un ciclo RF

I parametri di plasma sopra riportati si riferiscono a medie

su valori campionati per un periodo di tempo corrispondente a diversi

cicli di radiofrequenza. È interessante vedere come questi parametri

variano all’interno di un singolo ciclo.

0.00 0.20 0.40 0.60 0.80 1.000.00

0.50

1.00

1.50

2.00posizione (cm)

cicli rf

1e+11 2.16e+14 (log) elettroni per m3

Figura 5.16 Variazione della densità degli elettroni in un ciclo RF

Come è evidente dalla figura, gli elettroni rispondono

istantaneamente al campo applicato, e penetrano completamente in

ognuno dei due sheath una volta per ogni ciclo.

Gli ioni invece mantengono una densità costante. La figura

5.17 conferma questa previsione per gli ioni O2+ (analoga situazione si

ha per gli altri ioni).

0.00 0.20 0.40 0.60 0.80 1.000.00

0.50

1.00

1.50

2.00posizione (cm)

cicli rf

8.88e+13 1.33e+15 (log) ioni o2+ per m3

Figura 5.17 Densità degli ioni O2+ in un ciclo RF

94

Analogo andamento ci si aspetta dall’energia e dalla

velocità, ovvero le informazioni interessanti provengono soltanto dallo

studio di questi parametri negli elettroni.

0.00 0.20 0.40 0.60 0.80 1.000.00

0.50

1.00

1.50

2.00posizione (cm)

cicli rf

1 15 (log) energia elettroni (eV)

posizione (cm)

cicli rf

Figura 5.18 Energia degli elettroni in un ciclo RF

0.00 0.20 0.40 0.60 0.80 1.000.00

0.50

1.00

1.50

2.00posizione (cm)

cicli rf

-5e+05 5e+05 velocità degli elettroni

Figura 5.19 Velocità degli elettroni in un ciclo RF

95

5.6 Funzione di distribuzione dell’energia degli elettroni (EEDF)

Come detto in precedenza, ci si aspetta in questo tipo di

plasma una eedf non maxwelliana. Questo perché durante la

contrazione dello sheath e il conseguente avanzamento del fronte

elettronico, vengono prodotti degli elettroni “caldi” ad alta velocità. La

eedf dovuta a questi elettroni si combina con quella degli elettroni

“freddi”, portando ad un cambiamneto di pendenza nella funzione di

distribuzione.

La figura 5.20 riporta la eedf in diverse posizioni.

Figura 5.20 eedf in diverse posizioni

Si evince dalla figura che la eedf, per quanto non maxwelliana,

non mostra una elevata differenza tra le due temperature, e questa

differenza tende a diminuire allontanandosi dallo sheath, dove gli

10 -3

10 -2

10 -1

0,0 100 1,0 101 2,0 101

eedf

x= 0.12 cmx= 0.36 cmx= 0.6cmx= 0.84 cm

eedf (eV-3/2 )

energia (eV)

96

elettroni sono più energetici e le due temperature si equilibrano

maggiormente.

Si nota una competizione fra fenomeni di collisione

inelastica, che tendono a fare abbassare la cosiddetta "coda" della

distribuzione ad alta energia, ed il processo di scambio di energia tra

elettroni e campo di sheath, che produce un innalzamento della coda.

Tra i processi di collisione, quelli inelastici sono

naturalmente più efficaci poiché sottraggono all'elettrone gran parte

della sua energia. La diversità delle sezioni d'urto dei processi e la

disomogeneità del campo elettrico nella scarica RF generano la

ricchezza di comportamenti osservabile nelle figure.

Spostandosi dagli sheath verso il bulk si nota un

innalzamento della coda della eedf, ovvero il numero di collisioni

inelastiche che subiscono gli elettroni è minore mentre tende a

prevalere lo scambio di energia.

0.00 0.20 0.40 0.60 0.80 1.000

10

20

30

40energia(eV)

cicli rf

1e-05 0.3 (log) eedf (eV -3/2)

Figura 5.21 EEDF in un ciclo RF

Nel corso di un ciclo RF si nota una lieve modulazione

della eedf per alti valori di energia, corrispondente agli elettroni “caldi”

prodotti nelle ondulazioni dello sheath.

97

5.7 Funzione di distribuzione di energia degli ioni (IEDF)

Il modello del trasporto degli ioni positivi attraverso gli

sheats è basato su tre assunzioni fondamentali, che sono verificate in

questa simulazione: - la densità ionica ni nello sheath è costante nel tempo;

- gli elettroni hanno mobilità infinita: reagiscono istantaneamente al

campo elettrico oscillante;

- il processo dominante tra ione e neutra nello sheath è lo scambio di

carica.

Le prime due assunzioni sono giustificate dal fatto che la

frequenza rf usata in questa simulazione è alta rispetto alla frequenza di plasma degli ioni ωp+ ma è bassa rispetto alla frequenza di plasma degli

elettroni ωp-. Accade quindi che la densità degli ioni è determinata da

potenziali mediati nel tempo mentre gli elettroni seguono all’istante

l’alternare del potenziale (fig. 5.16).

Per quanto riguarda lo scambio di carica, da un’analisi

della sezione d’urto e dell’energia degli ioni si vede che anche a

pressioni relativamente basse (100 mTorr) e temperatura inferiore ai

500 °K il libero cammino medio è dell’ordine del millimetro, ovvero il

processo è molto probabile.

In seguito alla modulazione del potenziale nello sheath gli

elettroni dal plasma penetrano periodicamente nello sheath e compensano il flusso di ioni positivi sul catodo. Sia xe(t) la posizione

del fronte di elettroni. Per l’infinita mobilità degli elettroni, il campo elettrico E(x) è zero per x≤xe e il potenziale V(x) è uguale al potenziale

di plasma. Per x>xe la densità elettronica ne(x) va a zero e la

distribuzione del potenziale è determinata dall’equazione di Poisson,

98

)(0

2

2

xnedx

Vdiε

−= . (5.1)

I processi di scambio di carica portano alla creazione di ioni

termici nello sheath lasciando invariato il flusso ionico totale. Se uno

ione viene creato da un processo di scambio di carica nella posizione S0, raggiungerà il catodo con una probabilità exp[-α(1-S0)] dove α

=d/λ è il rapporto tra lo spessore dello sheath ed il libero cammino medio. L’energia di impatto Eimp(S0,,ϕ0) di tale ione è una funzione

della posizione di partenza S0 e dell’angolo di fase rf ϕ0 al momento

della creazione. Se invece uno ione entra nello sheath dal plasma e

raggiunge il catodo senza collisioni, l’energia di impatto sarà funzione della velocità iniziale v0 e nuovamente dell’angolo di fase.

La sezione d’urto dello scambio di carica può ritenersi pressochè

indipendente dall’energia, per cui la velocità di creazione di ioni

termici è costante nello sheath. Ne segue una espressione per la

distribuzione dell’energia del tipo :

I(E)dE ≈ α exp −α 1− S0( )[ ]

Eimp∈ E, E+ dE[ ]∫ dS0dϕ0 +

+ exp(−α ) f (v0 )dv0 dϕ0Eimp∈ E, E+ dE[ ]

∫ (5.2)

dove f(v0) è la distribuzione di velocità degli ioni normalizzata ad

1 in corrispondenza del confine tra plasma e sheath. Il primo termine

descrive il contributo degli ioni creati nello sheath, il secondo

rappresenta gli ioni che entrano nello sheath dal plasma e raggiungono

l’elettrodo senza collisioni. I due tipi di ioni sono descritti come ioni

99

secondari e primari rispettivamente. Il contributo degli ioni primari

porta ad un profilo “a sella” della IEDF, dove i due picchi

rappresentano il massimo e il minimo dell’energia raggiunta dagli ioni

che seguono in misura minima le oscillazioni del potenziale applicato.

Come detto sopra, l’energia d’impatto degli ioni secondari

dipende dalla posizione e dall’angolo di fase di partenza. Ci si aspetta

un picco nella distribuzione di energia degli ioni quando si annullano le derivate prime dEimp/dS0 e dEimp/dϕ0. La dipendenza di Eimp da ϕ0 è

mostrata in fig.5.22 [40]. Per S0 fissata, l’energia di impatto è una

funzione periodica di ϕ0. Oltre ai due estremi Emin ed Emax si trova una

regione centrata attorno a ϕ0=0 in cui dEimp/dϕ0 =0. Questa regione ∆ϕ0

corrisponde a ioni che sono stati creati nella regione oltre il fronte

elettronico in cui il campo elettrico totale è uguale a zero. Per questi

ioni quindi c’è un ritardo tra la creazione e l’accelerazione, che porta l’energia d’impatto ad esere indipendente dall’angolo di fase ϕ0 quando

quest’ultimo è circa uguale a zero. Ci sono quindi tre diverse energie

che soddisfano la condizione necessaria per avere un picco nella IEDF dEimp/dϕ0 =0. Queste tre energie sono indicate con Emax , Emin e Eϕ=0.

Figura 5.22 Dipendenza dell’energia ionica dall’angolo di fase iniziale

100

La fig. 5.23 mostra la dipendenza di Emax , Emin ed Eϕ=0

dalla posizione di partenza S0. Mentre le prime due decrescono in maniera continua all’aumentare di S0, Eϕ=0 mostra una serie di estremi

per i quali sono verificate entrambe le condizioni dEimp/dS0=0 e

dEimp/dϕ0=0. A questi estremi si possono far corrispondere picchi nella

IEDF, e ricondurre in ultima analisi l’origine dei picchi agli ioni creati per ϕ0=0, ovvero quando il campo elettrico nello sheath è trascurabile.

Ovviamente il numero di picchi è strettamente correlato al

numero di cicli rf impiegati dagli ioni per attraversare lo sheath, quindi

è possibile ottenere una informazione diretta sul tempo di transito

osservando la IEDF.

Figura 5.23 Dipendenza dell’energia ionica dalla posizione di partenza e dal

tempo di transito

101

Le figure riportano le funzioni di distribuzione dell’energia degli ioni O2

+ risultanti dalla simulazione, in diverse posizioni e per due temperature.

Figura 5.24 Picchi nella IEDF degli O2+ in prossimità delle pareti

Figura 5.25 Picchi nella IEDF degli O2+ allontanandosi dallo sheath

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

8,0 100

1,6 101

2,4 101

3,2 101

iedf ioni O2

+

x= 0.4 cmx= 0.48 cmx= 0.56 cm

iedf (eV

-3/2 )

energia (eV)

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

8,0 100 1,6 101 2,4 101 3,2 101

iedf ioni O2

+

x= 0.12 cm

x= 0.2 cm

x=0.28 cm

iedf (eV

-3/2)

energia (eV)

102

Figura 5.26 Progressiva perdita di struttura della IEDF degli O2+ verso il bulk

Dalla figura 5.24 è evidente che i picchi hanno maggiore

intensità per energie minori, in accordo con tutte le distribuzioni di

energia in questi plasmi.

I picchi sono inoltre più alti e più ravvicinati man mano

che ci si allontana dalla parete: questo fenomeno è dovuto alla diversa

velocità degli ioni lungo lo sheath. La regolarità della variazione di

spaziatura tra picchi successivi per diverse posizioni è indice di una

correlazione di fase tra velocità di ingresso nello sheath e angolo di

fase RF.

Le figure 5.35 e 5.26 mostrano la perdita della struttura

della IEDF quando ci si allontana dallo sheath: a centro plasma gli ioni

sono praticamente tutti termici, come si evince anche dai grafici

sull’energia.

Per quanto riguarda lo studio a due diverse temperature, si

trova una diminuzione dell’ampiezza dei picchi ed un aumento della

10-2

10-1

100

0 100

4 100

8 100

1 101

iedf ioni O2

+

x = 0.64 cmx = 0.76 cmx = 0.88 cm

iedf (eV

-3/2 )

energia (eV)

103

loro spaziatura nel caso di temperatura maggiore, similmente a quanto

avviene per una temperatura fissata avvicinandosi alle pareti. Questo

conferma che la spiegazione di questo fenomeno risieda ancora una

volta nella energia acquistata dagli ioni a causa della maggiore

temperatura.

La IEDF degli ioni O+ mostra dei picchi molto più stretti e

frequenti rispetto a quelli degli O2+. La maggiore frequenza è indice

della maggiore velocità del processo di scambio di carica con l’atomo

O. I picchi sono inoltre molto ben definiti grazie alla maggiore mobilità

degli ioni leggeri.

Figura 5.27 Picchi nella IEDF degli ioni O+ nello sheath

Nella figura 5.27, in corrispondenza dei valori più alti di

energia si nota la “sella” dovuta alla modulazione dello sheath, a cui gli

ioni O+ sono più sensibili data la loro minore massa.

Allontanandosi dallo sheath si nota una “spalla” per i

valori più alti di energia, che regredisce sempre più verso valori minori

man mano che si perde la struttura a picchi (figure 5.28 e 5.29). Questa

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

2,5 10-2

3,0 10-2

1 101

2 101

3 101

4 101

iedf ioni O+

x = 0.24 cmx = 0.32 cm

iedf (eV

-3/2 )

energia (eV)

104

caratteristica mostra chiaramente come la probabilità di dissociazione

degli ioni creati per scambio di carica nello sheath sia molto minore

della probabilità di dissociazione degli ioni nel plasma.

Figura 5.28 IEDF degli ioni O+ e progressivo rientro verso valori più bassi di energia

Figura 5.29 Perdita di struttura e assestamento verso valori di energia termica per gli ioni O+ a centro plasma

0,0 100

2,0 10-2

4,0 10-2

6,0 10-2

8,0 10-2

1,0 10-1

0 100

1 101

2 101

3 101

iedf degli ioni O+

x = 0.44 cmx = 0.56 cm

iedf (eV

-3/2 )

energia (eV)

0,0 100

1,0 10-1

2,0 10-1

3,0 10-1

4,0 10-1

5,0 10-1

6,0 10-1

7,0 10-1

4 100

1 101

2 101

3 101

iedf degli ioni O+

x = 0.68 cmx = 0.8 cm

iedf (eV

-3/2)

energia (eV)

105

Per quanto detto all’inizio del paragrafo, la IEDF degli ioni

negativi non è soggetta ad una struttura a picchi. Le figure 5.30 e 5.31

confermano questa ipotesi.

Figura 5.30 IEDF degli ioni O2-

Figura 5.31 IEDF degli ioni O-

10-3

10-2

10-1

100

101

102

2,0 100

6,0 100

1,0 101

1,4 101

iedf degli ioni O-

x = 0.12 cmx = 0.4 cmx = 0.68 cmx = 0.96 cm

iedf (eV

-3/2 )

energia (eV)

10 -3

10 -2

10 -1

10 0

10 1

10 2

4 100 1 101 2 101

iedf ioni O2

-

x = 0.16 cmx = 0.4 cmx = 0.6 cmx = 0.88 cm

iedf (eV

-3/2)

energia (eV)

106

Gli ioni O- si stabiliscono verso valori termici di energia verso il

centro del plasma, mentre gli ioni O2- mantengono una IEDF pressochè

costante nelle varie posizioni essendo maggiormente confinati nel

centro dal campo ambipolare.

5.8 Potenza depositata e reazioni

Al termine sorgente nell’equazione del calore contribuisce

l’energia depositata o assorbita in ogni reazione che si verifica nella

simulazione. La raccolta è stata effettuata distinguendo tra l’energia

depositata nelle collisioni che hanno come protagonisti gli elettroni,

quella depositata nelle reazioni di scambio di carica tra ioni, e quella

depositata in ogni tipo di reazione tra specie pesanti. Nell’equazione

del calore il termine S compare come flusso di energia, ovvero come

potenza per metro cubo. La potenza depositata in funzione della

posizione nei diversi canali di raccolta è riportata nelle figure seguenti.

Figura 5.32 Potenza depositata nelle reazioni da impatto elettronico

102

103

104

0 100

1 10-2

2 10-2

potenza depositata (W/m

3)

posizione (m)

107

Figura 5.33 Potenza depositata dagli ioni per scambio di carica

Figura 5.34 Potenza depositata dalle reazioni tra specie pesanti

102

103

104

105

0 100

1 10-2

2 10-2


3)

posizione (m)

-2,0 102

0,0 100

2,0 102

4,0 102

6,0 102

8,0 102

0 100

1 10-2

2 10-2


3)

posizione (m)

108

La forma della curva che descrive la raccolta energetica

nelle reazioni da impatto elettronico richiama la distribuzione di

energia e di posizione degli elettroni. Allo stesso modo la figura 5.33

conferma che lo scambio di carica avviene nello sheath, laddove gli

ioni riescono a raggiungere energia sufficiente a fare avvenire la

reazione; l’energia depositata al centro del plasma è infatti nulla.

Che lo scambio di carica sia praticamente l’unica reazione

(a parte la collisione elastica tra elettroni e tra elettroni e specie pesanti)

lo si vede dalla figura 5.34, che rappresenta la somma della potenza

raccolta da tutte le reazioni tra specie pesanti. Tutte queste reazioni,

siano esse esotermiche o endotermiche, avvengono nella zona centrale

della scarica.

Sono riportati in seguito i grafici relativi alle singole

reazioni.

Figura 5.35 Potenza depositata

La figura 5.35 mostra come l’O2 venga rigenerato al centro

del plasma a causa delle ricombinazioni. Le prime due reazioni

0,0 100

5,0 103

1,0 104

1,5 104

0 100

1 100

2 100

O2- + O2+ -> 2O2O- + O2 -> O2 + OO- + O -> O2 + e-02- + O+ -> O2 + O


3)

posizione (m)

109

depositano molta più potenza a causa dell’elevata concentrazione delle

specie protagoniste. La reazione O2-+O+ → O2 +O pur essendo

fortemente esotermica, deposita poca potenza a causa della bassa

concentrazione di reagenti.

Figura 5.36

La scarsità delle specie protagoniste ed il basso guadagno

(0.61 eV) fa si che la reazione di figura 5.36 depositi poca potenza.

0,0 100

5,0 101

1,0 102

1,5 102

2,0 102

2,5 102

3,0 102

0 100

5 10-1

1 100

2 100

2 100

O2- + O -> O3 + e-

potenza depositata (J/m

-3 )

posizione (m)

110

Figura 5.37

Figura 5.38

Le reazioni di figura 5.38 depositano una quantità di

potenza trascurabile.

0,0 100

2,0 103

4,0 103

6,0 103

8,0 103

1,0 104

0 100

5 10-1

1 100

2 100

2 100

O2- + O2+ -> O2 + 2O


3)

posizione (m)

0,0 100

1,0 100

2,0 100

3,0 100

4,0 100

5,0 100

6,0 100

0 100

5 10-1

1 100

2 100

2 100

O+ + O- -> 2O2e- + O+ -> O + e-


3)

posizione (m)

111

Figura 5.39

Le reazioni di figura 5.39, pur non essendo fortemente

endotermiche, assorbono molta potenza a causa dell’elevata

concentrazione degli ioni O- e degli stati metastabili dell’ossigeno

molecolare.

-1,2 103

-1,0 103

-8,0 102

-6,0 102

-4,0 102

-2,0 102

0,0 100

0,0 100

5,0 10-1

1,0 100

1,5 100

2,0 100

O- + O2* -> O2- + O

O- + O2** -> O + O2 + e-

potenza depositata (J / m

3)

posizione (m)

112

Figura 5.40

Si vede che le reazioni endotermiche assorbono molta

meno potenza di quanta ne depositino quelle esotermiche. Il resto della

potenza viene rilasciata per irraggiamento.

Al solito, in un ciclo RF variano l’energia, la densità e la

velocità dei soli elettroni. Pertanto varierà anche la potenza depositata

in seguito a collisioni da impatto elettronico.

-3,0 10-1

-2,5 10-1

-2,0 10-1

-1,5 10-1

-1,0 10-1

-5,0 10-2

0,0 100

0,0 100

5,0 10-1

1,0 100

1,5 100

2,0 100

2e- + O2+ ->O2 + e-

e- + O3 -> O2- + O

potenza depositata (J / m

3)

posizione (m)

113

Figura 5.41 Variazione della potenza depositata da reazioni da impatto

elettronico in un ciclo RF

La figura 5.41 rivela chiaramente l’avanzamento e

l’indietreggiamento del fronte elettronico dentro e fuori dallo sheath.

La lunga coda e lo spostamento del massimo verso il centro della

scarica per 1/4 e 3/4 del ciclo RF è dovuto agli elettroni “caldi” prodotti

quando il fronte elettronico penetra completamente nello sheath.

0,0 100

5,0 103

1,0 104

1,5 104

2,0 104

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

.25

.5

.751media


3)

posizione (cm)

114

5.9 Ionizzazione

Le ionizzazioni, anch’esse dovute ad impatto elettronico,

non si verificano ad un tasso uniforme, infatti un un ciclo RF variano la

velocità e la densità degli elettroni, mentre rimane pressochè invariata

la densità delle molecole o atomi che vengono ionizzati. La figura 5.41

riporta la variazione del tasso di ionizzazione dell’O2 in un ciclo RF.

0.00 0.20 0.40 0.60 0.80 1.000.00

0.50

1.00

1.50

2.00posizione (cm)

cicli rf

9e+17 8e+20 tasso di ionizzazione

Figura 5.42 Variazione del tasso di ionizzazione dell’O2 in un ciclo RF

Figura 5.43 Tasso di ionizzazione dell’O2 in quattro frazioni di ciclo RF in funzione della posizione

1021

1022

1023

0,0 100

5,0 10-3

1,0 10-2

1,5 10-2

2,0 10-2

0.250.50.751media

tasso di ionizzazione (m

-3 s

-1 )

posizione (m)

115

CONCLUSIONI

Il codice numerico messo a punto nell’ambito di questa

tesi implementa un nuovo modello auto-coerente per una scarica RF in

ossigeno in un reattore a piatti paralleli. Il numero delle specie presenti

e l’attenzione rivolta alle reazioni in gioco lo rendono completo anche

dal punto di vista della teoria cinetica e della cinetica chimica.

È stata evidenziata una ricca dinamica caratterizzata dal

non equilibrio locale, non solo tra le temperature elettroniche e ioniche,

ma anche all’interno delle stesse distribuzioni traslazionali. Esse sono

infatti lontane dalla distribuzione di Boltzmann, come dimostrano

soprattutto i picchi dinamici trovati nella distribuzione di energia degli

ioni positivi nella regione di sheath.

È emersa inoltre l’impossibilità si trattare separatamente il

trasporto reattivo delle specie neutre e quello delle cariche elettriche,

rendendo necessario l’accoppiamento tra il problema di Vlasov-

Boltzmann-Poisson multispecie e l’equazione di reazione e diffusione

per il trasporto delle particelle neutre.

L’importanza di una buona simulazione in cui i parametri

di plasma rispondano correttamente alle variazioni delle condizioni

applicabili dall’esterno rende questo modello un’ottima base

‘sperimentale’ in senso lato per studi di carattere fondamentale.

Per quanto riguarda gli sviluppi futuri è prevista

l’estensione a due dimensioni spaziali della simulazione per il caso di

simmetria cilindrica, con tre dimensioni nello spazio delle velocità.

116

A conclusione di questo lavoro voglio ringraziare tutti i docenti,

ricercatori, dottorandi e laureandi presso il Centro di Studio per la

Chimica dei Plasmi, per la disponibilità e la simpatia che mi hanno

sempre dimostrato. Ringrazio poi quanti, volendomi bene, mi sono stati

vicini in particolare nel corso di questo impegno.

117

APPENDICE LIVELLI ENERGETICI DELLA MOLECOLA DI

OSSIGENO

Energia potenziale in funzione della distanza internucleare

119

BIBLIOGRAFIA [1] Longo S , Capitelli M 1994 Plasma Chem Plasma Proc 14 1 [2] Longo S , Capitelli M and Hassouni K 1997 J.Physique IV C4 7

271 [3] Longo S, Capitelli M, Hassouni K 1998 J. Thermophysics and

Heat Transfer 12 473 [4] Iasillo D Tesi di Laurea in Fisica , Università degli Studi di Bari

1994 [5] Hockney RW 1966 “Computer Simulation of anomalous plasma

diffusion and numerical solution of Poisson’s equation” Phys.Fluids 9 1826

[6] Woolfson J, Pert C 1999 An Introduction to Computer Simulation, Oxford University Press

[7] Boeuf JP , Belenguer P in Non Equilibrium Processes in Partially Ionized Gases, M Capitelli and JM Bardsley, Eds. New York: Plenum 1990

[8] Lieberman MA 1988 IEEE Trans. Plasma Sci., 16 638 [9] Surendra M , Graves DB , Jellum JM 1990 Appl. Phys. Lett. 56 11

1022 [10] Godyak VA 1972 Sov. Phys. Tech. Phys. 16 7 1073 [11] Kushner MJ 1986 IEEE Trans. Plasma Sci PS 14 188 [12] Belenguer P Tesi di Dottorato, Tolosa 1990 [13] Birdsall CK, Langdon AB 1985 Plasma Physics via Computer

Simulation McGraw-Hill NY [14] Birdsall CK 1991 “Particle in Cell Charged-Particle Simulations,

Plus Monte Carlo Collisions With Neutral Atoms, PIC-MCC” IEEE Trans. Plasma Sci 19 2 65

[15] Vahedi V , Surendra M 1995 “A Monte Carlo Collision Model for the Particle in Cell Method: applications to argon and oxygen discharges” Comput. Phys. Commun. 87 179

[16] Cohen BI , Langdon AB , Hewett DW and Procassini RJ 1989 J.Computat. Phys. 81 151

[17] Skullerud HR 1968 J. Phys. D: Appl. Phys. 1 1567 [18] Spanier J and Gelbard EM 1969 Monte Carlo Principles and

Neutron Transport Problems (Reading, MA: Addison Wesley) [19] Eastwood SW, Hockney RW 1981 Computer Simulation using

Particles, McGraw-Hill, NY [20] Longo S 2000 “Monte Carlo models of electron and ion transport

in non-equilibrium plasmas” Plasma Sources Sci. Technol. 9 468 [21] Phelps AV 1985 Tabulations of cross sections and calculated

transport and reaction coefficients for electron collisions with O2, JILA Information Center Report, University of Colorado, Boulder

120

[22] Rapp D, Briglia D 1965 J. Chem. Phys. 43 1480 [23] Vahedi V, Birdsall CK, Lieberman MA, DiPeso G, Rognlien TD

1993 “Capacitive RF discharges modelled by particle-in-cell Monte Carlo simulation” Plasma Sources Sci. Technol. 2 261, 273

[24] Smirnov BV 1982 Negative Ions (McGraw-Hill, New York) [25] Itikawa Y, Hara S, Kato T, Nakazaki S, Pindzola MS, Crandall

DH 1985 “Electron impact cross sections and rate coefficients for excitation of C and O Ions”, Atomic Data and Nuclear Data Tables 33 149

[26] Lieberman MA, Lichtenberg AJ 1994 Principles of plasma discharges and materials processing, John Wiley & Sons, NY

[27] Lee C, Graves DB, Lieberman MA, Hess DW 1994 J. Elettrochem. Soc., 141 1546

[28] Eliasson B, Kogelschatz U 1986 Basic data for modelling of electrical discharges in gases: oxygen, Report KLR 86-11C, Brown Boveri Konzernforschung, CH-5405 Baden

[29] Capitelli M, Molinari E 1980 “Kinetics of dissociation processes in plasmas in the low and intermediate pressure range”, Topics in Current Chemistry 60

[30] Kossyi IA, Kostinsky AY, Matveyev AA, Silakov VP 1992 Plasma Sources Sci. Technol. 1 207

[31] Perrin J 1993 “Modelling of the power dissipation and rovibrational heating and cooling in SiH4-H2 RF glow discharges” J. Phys. D: Appl. Phys 26 1662

[32] Baggerman JAG, Visser RJ, Collart EJH 1994 “Power dissipation measurements in a low-pressure N2 radio-frequency discharge” J. Appl. Phys. 76 2

[33]Visser RJ 1989 “Determination of the power and current densities in argon and oxygen plasmas by in situ temperature measurements” J. Vac. Sci. Technol. A 7 2 189

[34] Bogaerts A, Gijbels R 2000 “Modelling of radio-frequency and direct current glow discharges in argon” J.Tech.Phys. 41 1 183

[35] Tsederberg NV 1965 Thermal Conductivity of Gases and Liquids (MIT Press Cambridge, Mass)

[36] Runchal AK 1972 Int. J. Num. Methods Eng. 4 541 [37] Patankar SV 1980 Numerical Heat Transfer and Fluid Flow,

McGraw-Hill, NY [38] Press WH, Flannery BP, Teukolsky SA, Vetterling WT 1986

Numerical Recipes-The Art of Computer Simulation , Cambridge [39] Kohler K, Coburn JW, Horne DE, Kay E 1985 “Plasma potentials

of 13.56 MHz rf argon glow discharges in a planar system”, J. Appl. Phys. 57 59

121

[40] Surendra M, Graves DB 1991 “Particle simulations of Radio-Frequency Glow Discharges”, IEEE Transactions on Plasma Science, 19 144

[41] Wild C, Koidl P 1991 “Ion and electron dynamics in the sheath of radio-frequency glow discharges”, J. Appl. Phys. 69 2909

UNIVERSITA’ DEGLI STUDI DI BARI - uniba...

Documents

Transcript of UNIVERSITA’ DEGLI STUDI DI BARI - uniba...