Uvod do fyziky plazmatu interaktivně Část II: Praktická cvičení Adam Obrusník, Lenka Zajíčková Copyright © 2013 Masaryk University Licensed under the Creative Commons Attribution-NonCommercial 3.0 Unported License (the "License"). You may not use this file except in compliance with the License. You may obtain a copy of the License at http: //creativecommons .org/licenses/by-nc/3.0. Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License. First printing, September 2014 1 Úvod........................................................ 5 1.1 Čím tento text je 5 1.2 Čím tento text není 5 1.3 Napojení na přednášku 6 1.3.1 2. a 3. týden: Příprava........................................... 6 1.3.2 4. až 6. týden: Částice v polích.................................... 6 1.3.3 7. a 8. týden: Zpracování dat..................................... 6 1.3.4 9. a 10. týden: Rozdělovači funkce................................. 6 1.3.5 11. až 13. týden: závěrečná úloha ................................. 6 2 Základy Matlabu ............................................ 7 2.1 Instalace Matlabu na MUNI 7 2.1.1 Licence ...................................................... 7 2.1.2 Získání instalačních souborů ...................................... 7 2.1.3 Instalace ..................................................... 8 2.1.4 Spuštění Matlabu............................................... 9 2.2 Základy Matlabu 9 2.2.1 Proměnné a indexování ........................................ 10 2.2.2 M-soubory ................................................... 11 2.2.3 Kam dál..................................................... 12 3 Pohyb částic v el. a mag. polích ............................ 15 3.1 Částice v konstantních elektrických a magnetických polích 15 3.1.1 Související rovnice ............................................. 15 3.1.2 Řešení soustavy ODR v Matlabu .................................. 16 3.1.3 Implementace................................................ 17 3.1.4 Cvičení...................................................... 19 3.2 Van Allenovy radiační pásy 21 3.2.1 Související rovnice ............................................. 21 3.2.2 Implementace................................................ 22 4 Zpracování dat............................................. 29 4.1 Úvod 29 4.2 Smyčky a podmínky v Matlabu 29 4.2.1 Smyčka for................................................... 30 4.2.2 Podmínky if... elif... else........................................ 32 5 Interakce částic v plazmatu................................. 35 5.1 Maxwell-Boltzmannovo rozdělení 35 5.1.1 Střední rychlost................................................ 35 5.1.2 Tvar Maxwell-Boltzmannova rozdělení.............................. 38 5.2 Časový vývoj rozdělovači funkce 38 5.3 Rychlostní koeficienty 42 6 Rovnováha částic v plazmatu............................... 45 6.1 Formulace problému 45 6.1.1 Srážkový člen................................................. 45 6.1.2 Reakční schéma .............................................. 47 6.2 Implementace 48 6.3 Parametrická studie 52 6.4 Závěr 53 Bibliography ............................................... 55 Books 55 Articles 55 Online resources 55 1.1 Čím tento text je Tento text je doplňujícím materiálem ke cvičením z předmětu ¥5170: Úvod do fyziky plazmatu. Zatímco přednáška obsahuje mnoho teoretických odvození a úvah, cílem tohoto studijního materiálu je zejména zprostředkovat intuitivní zkušenost s fyzikou plazmatu pro studenty, kteří se s plazmatem setkávají poprvé. V průběhu textu budete jako studenti postaveni před řadu úloh, z nichž většina musí být vyřešena s pomocí počítače. Mimo jiné se naučíte jak řešit soustavy obyčejných diferenciálních rovnic, jak hromadně zpracovat data nebo jak numericky integrovat. Většinu schopností a dovedností, které se v tomto kurzu naučíte upotřebíte nejen ve fyzice plazmatu, ale také v jiných oblastech fyziky. Praktická cvičení jsou navržena tak aby studentům pomohly vizualizovat a intuitivně pochopit těžce představitelné jevy probíhající v plazmatu. V textu je zohledněno, že většina studentů nemá rozsáhlé předchozí zkušenosti s programováním. Z tohoto důvodu je veškerý zdrojový kód bohatě komentován a dokumentován a studenti by měli být schopni pochopit základy programování v Matlabu z praktických příkladů v tomto textu. 1.2 Čím tento text není V první řadě musí být řečeno, že tento text si neklade za cíl být vyčerpávající učebnicí programování v Matlabu. Existuje nespočet kvalitních knih a manuálů týkajících se Matlabu nebo jeho balíčků. Ty nejrelevantnější a nejsnáze dostupnejšou vyjmenovány v oddílu 2.2.3. Očekávejte, že textu může místy být příliš stručný, krkolomně napsaný nebo může mít nečekané logické členění. Důvodem je, že tento text není navržen tak, aby byl použitelný sám o sobě a nejlépe bude fungovat v kombinaci se cvičeními k předmětu F5170: Úvod do fyziky plazmatu. Tento text také není vhodným zdrojem, který by jeho uživatel měl citovat v jakýchkoliv vědeckých publikacích nebo diplomových pracích. Vždy je nejlepší odkázat se přímo na knihy nebo články, ze kterých tento text vychází. Ačkoliv se může zdát, že hledat původní reference jen abyste ověřili jeden vzorec nebo obrázek, je zbytečně složité, jedná se o schopnost, která vám usnadní vaši budoucí vědeckou práci. Dále je nutno zmínit, že tento text není zaměřen na matematické principy numerických metod. Řešiče diferenciálních rovnic, které jsou v tomto textu použity, jsou provazovány za černé 6 Chapter 1. Úvod skříňky a princip jejich fungování není dále analyzován. Autoři se také vyhnuli transformaci rovnic do bezrozměrného tvaru, protože používání smysluplných fyzikálních jednotek vám umožní kromě kvalitativního náhledu také získání kvantitativní představy o důležitých veličinách ve fyzice plazmatu. 1.3 Napojení na přednášku Přednáška je rozložena do 13 týdnů podzimního semestru na Masarykově Univerzitě. Tento text odráží obsah přednášky a zhruba jednou za dva týdny byste měli pokročit o jednu kapitolu. V prvním týdnu semestru cvičení většinou neprobíhá. 1.3.1 2. a 3. týden: Příprava Vzhledem k tomu, že přednáška vyžaduje oproti cvičení určitý náskok, jsou celé první dva týdny vyhrazeny tomu, abyste si nainstalovali potřebný software a naučili se s ním zacházet. Jak uživatelé systému Widnows tak uživatelé Linuxových systémů si nainstalují univerzitní licenci Matlabu (vizte kapitolu 2). Studenti by si také měli nastudovat druhou část této kapitoly, aby získali základní představu o logice a syntaxi Matlabu. 1.3.2 4. až 6. týden: Částice v polích V týdnech 4 až 6 byste měli pochopit všechny Matlabové programy představené v kapitole 3 a v šestém týdny byste měli dokončit všechna relevantní cvičení. Cvičné programy se týkají pohybu částic v elektrických a magnetických polích a měly by vám pomoci pochopit povahu různých driftů v plazmatu. 1.3.3 7. a 8. týden: Zpracování dat Na příkladech z fyziky plazmatu se naučíte se základy hromadného zpracování dat. Konkrétně se jedná o interpolaci a integraci diskrétních dat a také hromadný převod dat do jiných jednotek. Ačkoliv se budou úlohy týkat dat souvisejících s fyzikou plazmatu, budou se vám získané znalosti jistě hodit i v jiných oblastech fyziky. 1.3.4 9. a 10. týden: Rozdělovači funkce Tato krátká kapitola by Vám měla pomoci pochopit rozdělovači funkce. Za tímto účelem vám bude poskytnuto několik programů, které jsou schopny vizualizovat časový vývoj rozdělovačích funkcí. Kromě toho budete analyzovat tvar Maxwell-Boltzmannova rozdělení za různých podmínek. 1.3.5 11. až 13. týden: závěrečná úloha V závěrečné úloze využijete schopností, které jste získali v předchozích kapitolách. Budete muset vytvořit případně dotvořit program, který řeší rovnováhu částic v laboratorním argonovém plazmatu. Budete pozorovat jak se s měnícím tlakem či teplotou elektronů nebo těžkých částic mění složení plazmatu. Instalace Matlabu na MUNI Licence Získání instalačních souborů Instalace Spuštění Matlabu Základy Matlabu Proměnné a indexování M-soubory Kam dál (l. Základy Matlabu Matlab (původně Matrix Laboratory) je komerční software a skřip to vací jazyk vyvíjený společností MaťhWorks. Na rozdíl od jiných tradičních programovacích jazyků umožňuje Matlab intuitivní práci s maticemi a vektory. Uživatelé si také mohou definovat své vlastní funkce s takřka libovolným vstupem a výstupem. 2.1 Instalace Matlabu na MUNI 2.1.1 Licence K roku 2014 vlastní Masarykova univerzita licenci nejnovějšrho Matlabu se současným přístupem pro nejméně 250 uživatelů. Aby instalace fungovala a Matlab šel spustit, musí být uživatel připojen k síti VPN Masarykovy univerzity nebo se nacházet v prostorách univerzity. Licence je omezena na výzkumné nekomerční účely. 2.1.2 Získání instalačních souborů Instalační soubory Matlab 2014b jsou k dispozici na intranetu Masarykovy univerzity. Instalační soubory mají zhruba 7.5 GB, ale v rámci MU (například z knihovny) by mělo stahování být extrémně rychlé. Od verze 2013 je Matlab k dispozici pouze pro 64bitové operační systémy. Pokud máte starší počítač vybavený 32bitovým procesorem, kontaktujte cvičícího. Návod 2.1.1 — Stažení souborů. Pro stažení souborů postupujte následovně: 1. V prohlížeči otevřete inet.muni.cz 2. Přihlaste se pomocí U CO a primárního hesla 3. Jděte do Software —> Nabídka softwaru 4. V seznamu najděte Matlab a klikněte Získat 5. Po odsouhlasení licenční smlouvy si do počítače uložte • souborlicense.dat, • instalační obraz, R201xx_Windows. iso pro Windows nebo R201xx_UNIX. iso pro Linuxové systémy. • Autorizační kód - uložte si jej například do textového souboru, instalátor jej bude vyžadovat. 8 Chapter 2. Základy Matlabu 2.1.3 Instalace Po stažení ISO obrazuje musíte připojit jako virtuální mechaniku a spustit instalaci. Alternativně můžete vypálit obraz ISO na DVD a nainstalovat Matlab z něj. Návod 2.1.2 — Připojení ISO ve Windows. Pravděpodobně nejčastěji používaný software na připojování ISO souborů (tedy vytvoření virtuální mechaniky, která bude obsahovat data uložená v ISO souboru) je pravděpodobně Daemon Tools Lite. Přestože je nejlepším softwarem pro připojování virtuálních mechanik, nutí Daemon tools Lite uživatele do instalace nechtěného software. Pokud se chcete vyhnout instalaci nechtěného panelu v prohlížeči, čtěte pozorně instrukce během instalace a ve správnou chvíli zvolte, že nechcete instalovat software třetích sttan. Po instalaci Daemon Tools připojte DVD s Matlabem. 1. Klikněte pravým tlačítkem na ikonu Daemon Tools v oznamovací oblasti. 2. Zvolte Virtual DVD —> Device 0: —> Mount Image. 3. Najděte na svém disku ISO obraz a potvrďte. 4. Nyní uvidíte DVD s Matlabem v Počítač Návod 2.1.3 — Připojení ISO v Linuxu. Většina moderních linuxových distribucí obsahuje zabudované funkce pro připojování ISO souborů. Většinou stačí kliknout pravým tlačítkem na ISO soubor a kliknout na "Mount" či "Připojit". Pokud tuto funkci ve vaší distribuci nemůžete najít, můžete použít program AcetonelSO, který by se měl nacházet v repozitářích. Pokud dáváte přednost příkazové řádce, nemusíte instalovat nic a stačí do terminálu napsat následující: mount -o exec R201xx_UNIX.iso /mnt/disk i kde /mnt/disk je adresář, kam chcete ISO soubor připojit. Instalátor Matlabu spustíte ve Windows přes soubor setup. exe a v Linuxu pomocí binárního souboru Jinstall. Oba tyto soubory jsou uloženy v kořenové složce instalačního DVD. Pokud se Vás instalátor zeptá, zda si přejete Matlab aktivovat, klikněte na Yes. Na Linuxu se doporučuje instalovat Matlab jako superuživatel. V opačném případě se nevytvoří odkaz v /usr/local/bin a nebude možné Matlab spouštět pomocí příkazu matlab z příkazové řádky. Samotná instalace Matlabu je přímočará a uživatelský přívětivá. Nejprve se Vás zeptá na Autorizační kód následně na umístění, kam se má Matlab instalovat a nakonec Vás požádá o licenční soubor. Těsně předtím, než se na disk začnou soubory kopírovat budete dotázáni na instalaci toolboxů. Toolboxy jsou dodatečné knihovny funkcí, které rozšiřují funkčnost Matlabu. Pokud potřebujete ušetřit místo na disku, nemusíte instalovat toolboxy z následujících kategorií: • Test and measurement, • Computational finance, • Computational Biology, • Code generation. 2.2 Základy Matlabu 9 2.1.4 Spuštění Matlabu Uživatelé Windows mohou Matlab spustit přes zástupce na ploše, uživatelé Linuxu pomocí příkazu matlab (pokud jste software instalovali jako superuživatel) nebo z jeho instalačního adresáře. /path/to/Matlab/R20xx/bin/matlab Po spuštění Matlabu se již ukáže hlavní okno, které je už identické na všech systémech. MATLAB R2014a is^ _ m _ b ©I ~p r-^ ■ {]. J |j dg, New Variable Analyze Code Open Variable ^ jff Run and Time Data Workspace Clear Workspace ^ ^ clear Commando ^ Library VARIABLE CODE SIMULINK New New Open y>j Compare Import Save Script ~ ▼ iSSj = SJmulink Layout Q Set Path jjjj Parallel -ENVIRONMENT ^ i±3 lp Ö / ► data » Skola » Teaching ► Pla5ma_FRMU * AdamllPPprog * Progl_partide_motion 3 Command Window _I Nam e L. ^ cline.rn J EBdrift.eps _J EBdrift.bd f*\ odefun.rn Q Output.eps solve.m (Script) Workspace q St tf ti 5] timespan im-1 lvalue 1000 1000x6 double [0;D;D;1D0;0;D] I OOOxI double 0.0025 0 IxlOOO double (i) New to MATLftB? Watch this Vi-:l-=■:■■ see E.-:amples, or read Getting 'Started, » solve » solve Figúre 2.1: Hlavní okno Matlabu Hlavní okno má tři důležitá podokna: • Current folder zobrazuje obsah složky, ve které právě pracujete. To jsou zejména všechny matlabové skripty, které můžete spustit. • Workspace obsahuje všechny proměnné (numerické hodnoty nebo matice), které jsou definovány. Dvojitým kliknutím na proměnnou ze zobrazí její hodnota a můžete ji editovat. • Command window je nejzásadnější částí. Zde můžete definovat nové proměnné, volat funkce nebo spouštět Matlabové skripty, o kterých se toho více dozvíte v dalším oddílu. 2.2 Základy Matlabu Tento oddíl stručně vysvětluje základy syntaxe a logiky Matlabu. Vysvětlení nejsou velmi vyčerpávající a rigorózní, protože se předpokládá, že si poznatky osvojíte na praktických příkladech níže. Na konci tohoto oddílu byste měli vědět zejména odpovědi na následující otázky: • Co je proměnná v Matlabu a jak ji definovat. • Jak přistupovat k prvkům vektorů a matic, jak vybírat submatice. • Jaký je rozdíl mezi soubory skriptů a soubory funkcí. • Co je v Matlabu funkce a jak ji volat. Pro některé z Vás budou informace v tomto oddíle triviální, ale tento text myslí i ty, kteří se zatím neměli možnost setkat s programováním. 10 Chapter 2. Základy Matlabu Proměnné a indexování Práce s proměnnými v Matlabu je stejně jednoduchá jako na papíře, jen je důležité si uvědomit, že Matlab provádí numerické operace, nikoliv analytické. Co to znamená v praxi? Zkuste si definovat proměnnou tím, že do Command window napíšete: v = 10 a předpokládejme, že tato proměnná vyjadřuje rychlost nějakého objektu. Nyní definujeme hmotnost našeho imaginárního objektu v kilogramech m = 3 a jeho kinetickou energii tedy spočítáme jako Ek = 0.5*m*v~2 Pokud nyní předefinujete rychlost na jinou hodnotu a podíváte se do podokna Workspace, uvidíte, že ačkoliv jsme změnili rychlost, hodnota proměnné Ek odpovídající kinetické energii se nezměnila. To je proto, že tato proměnná obsahuje pouze numerickou hodnotu a neobsahuje žádný odkaz na proměnné, pomocí kterých byla spočtena. Pokud chcete získat aktuální hodnotu Ek, musíte spustit odpovídající příkaz. Pracovat s maticovými a vektorovými proměnnými je v Matlabu také jednoduché. Předpokládejme, že chceme rotovat vektor v 3D kartézském prostoru v xy rovině. Nejprve tedy definujeme vektor, který budeme rotovat: r = [10, 20, 30] a následně matici otočení pomocí vestavěných funkcí sin () ad cos () M = [cos(pi/3) ,-sin(pi/3) ,0; sin(pi/3) , cos(pi/3) , 0; 0,0,1] Jak vidíte jsou matice definovány řádek po řádku. Prvky ve sloupcích jsou odděleny čárkami (,) a řádky jsou odděleny středníky (;). Vektorem se potom v podstatě rozumí matice 1 x n nebo n x 1 Nicméně, pokud se pokusíte Mar vynásobit a napíšete rrot = M*r dostanete chybu, inner dimensions do not agree. To je způsobeno tím, že r je řádkový vektor a matici 3x3 nelze vynásobit vektorem 1x3. Pro násobení tedy musíte napsat příkaz rrot = M*r' kde apostrof značí transpozici vektoru r. Jak jste si všimli, zobrazí se při definování proměnné její hodnota v okně příkazů {Command window). Abyste tomu v budoucnu předešli, ukončujte každý příkaz středníkem (;)■ Standardní operace, které můžete provádět na vektorech a maticích jsou sčítání (+), odčítání (-), násobení (*), dělení (/) a umocnění {"). Kromě toho existují v Matlabu operace po prvcích. Například násobení po prvcích (. *) dvou vektorů rozměru 1x3 funguje následujícím způsobem (1,4,5) .*(2,4,6) = (1*2,4*4,5*6) = (2,16,30). (2.1) Kromě násobení po prvcích existuje také dělení po prvcích (. /) umocňování po prvcích (."). Dále je důležité vědět, jak přistupovat k jednotlivým prvkům matic a vektorů nebo jak vytvořit sub-matice a sub-vektory. V případě vektorů můžete provádět následující operace 2.2 Základy Matlabu 11 a = [1,2,3,4,5,6,7] ; b = a(4) % ziska ctvrty prvek "a" = 4 b = a(l:3) % subvektor od prvku 1 do 3 = [1,2,3] b = a(4:end) °/0 subvektor od prvku 4 do konce = [4,5,6,7] V případě dvourozměrných matic vybíráte prvky podobně, potřebujete nicméně dva indexy kdy první udává řádek a druhý udává sloupec a = [1,2,3;4,5,6;7,8,9] ; b = a(l,2) °/0 prvek na souřadnicích 1x2 = 2 b = a(:,2) °/0 druhy sloupec = [2,5,8] b = a(3,:) °/. treti radek = [7,8,9] b = a(l:2,l:2) °/0 ziska submatici = [1,2;4,5] M-soubory Definování proměnných a práce s nimi v okně příkazů je sice užitečné, ale pokud by Matlab uměl jen to, byl by jen trochu lepší kalkulačkou. Typicky se Matlab (a jiné podobné skriptovací jazyky) používají pro automatizaci repetitivních nebo náročných úkolů, obsahujících určitou posloupnost operací. V Matlabu se posloupnosti příkazů ukládají do dvou typů souborů, do souborů skriptů a souborů funkcí. Trochu matoucí může být fakt, že oba typy souborů mají příponu .m. Soubory skriptů Soubor skriptu obsahuje posloupnost Matlabových příkazů, jeden příkaz na každém řádku. Tyto příkazy se postupně spouštějí, jeden za druhým. Ukážeme si použití skriptu na hypotetické situaci, kdy po delší dobu měříte určitou veličinu a chcete každý den generovat vždy aktuální graf se stejnou úpravou. V tomto případě je pohodlné vytvořit si skript Příklad skriptu v Matlabu day = [1 , 2 , 3 , 4, 5 , 7 , 8 , 10 , 12 , 15] ; °/„ days count = [450, 363, 291, 220, 180, 150, 120, 100, 90, 80]; % počet bakterii mc = max (count) ; °/0 najde maximum vektoru count norm_count = count/mc; °/0 provede normalizaci °/0 vykresli day na osy x a count na osu y plot(day, norm_count , 'rx - '); xlabel ( ' Days ') ; °/» nastavi popisek osy x ylabel ( ' Normalized cell count [CFU]'); °/» nastavi popisek osy y title ('Plasma sterilisation experiment') °/0 nastavi nadpis obrázku print -dpng -r72 ' experiment . png' ; °/0 ulozi obrázek do png Na prvních dvou řádcích jsou uložena měřená data. Proměnná day obsahuje jednotlivé dny měření zatímco proměnná count obsahuje měřené hodnoty v jednotlivé dny. Na řádcích 3 a 4 jsou data normalizována a od řádku 7 dále jsou vykreslena. V naší modelové situaci by Vám takový skript jistě ušetřil čas. 12 Chapter 2. Základy Matlabu Zdrojové kódy výše obsahují komentáře. Komentáře můžete a měli byste používat k anotování veškerého vašeho kódu, abyste si byli schopni vybavit, co který příkaz dělá, i po delší době. Znak procenta v Matlabu tedy dává najevo, že cokoliv za ním (až do konce řádku) není součástí příkazu a Matlab to bude ignorovat. Pokud chcete spustit skript, změníte Current folder na adresář, ve kterém je skript uložen a napíšete název skriptu (bez přípony .m) do okna příkazů. Cvičení 2.1 Napište skript, který vykreslí následující měření koncentrace ozonu v logaritmické škále. Dejte pozor na to, že funkce log O počítá přirozený logaritmus a pro spočítání logaritmu se základem 10 je třeba použít funkci loglOO polO£ ;a [mm] 0.0 0.1 0.2 0.3 0.5 0.6 konc 03 [m-3] 1.4 -1020 1.0 -1020 5.2-1019 1.6 • 1019 4.5-1018 9.1 • 1017 Soubory funkcí M-soubory obsahující funkce jsou poněkud odlišné od skriptů. V první řadě musí být zmíněno, že funkce už není jen prostá posloupnost příkazů ale je to objekt, který má vždy N vstupů a M výstupů, podobně jako funkce matematická. Pro příklad, funkce, která dokáže spočítat Larmorův poloměr a cyklotronovou frekvenci pro elektron by vypadala následovně: Matlab function example: gyro.m function [r, omega] = gyro(B, vp) i q = -1.602e-19; °/0 elementární náboj [C] 2 m = 9.109e-31; °/. hmotnost [kg] 3 4 r = m*vp/(abs(q)*B); °/0 Larmorův poloměr 5 omega = q*B/m; °/0 Cyklotronová frekvence 6 end °/0 konec funkce 7 Chování této funkce je dáno na prvním řádku. Proměnné r a omega jsou výstupní proměnné zatímco B a vp jsou vstupní proměnné. Název funkce na prvním řádku (gyro) se musí shodovat s názvem souboru. Tuto funkci můžete následně volat buď z okna příkazů nebo z jiného matlabového skriptu umístěného ve stejném adresáři. Pro vyzkoušení vyzkoušejte napsat [rO, wO] = gyro(le-4, le3) i do Vašeho okna příkazů. Ve Workspace by se měly objevit dvě nové proměnné, rO a wO. Cvičení 2.2 Napište funkci, která bude mít na vstupu koncentraci elektronů v m~3 a teplotu elektronů v eV. Na výstupu byla měla být Debyeova délka a plazmová frekvence. Potřebné vzorce můžete najít například v [Bit04]. ■ 2.2.3 Kam dál V tomto oddíle byly vysvětleny jen ty nejdůležitější základy programování v Matlabu. Pro detailnější interaktivní přehled syntaxe doporučují autoři nahlédnout na www.mathworks.com/help/matlab/getting-started-with-matlab.html. Stejně jako u ostatních počítačových dovedností platí, že Google je Váš kamarád. Pokud budete hledat odpověď na nějaký problém, všímejte si zejména odkazů na stackexchange.com. Vzhledem k rozšířenosti Matlabu je skoro jisté, že nejste první, kdo narazil na určitý problém a odpověď na Vás čeká na jednom z fór. 2.2 Základy Matlabu 13 Pokud budete chtít používat Matlab i nad rámec tohoto kurzu, například pro vlastní výzkumnou činnost, je vhodné seznámit se i s portálem MatlabCentral www.mathworks.com/matlabcentral. Tento portál obsahuje detailní dokumentaci Matlabových funkcí a většina z nich je doplněna příklady použití. Dále je na MatlabCentral dostupné diskuzní fórum plné zkušených uživatelů a repozitář užitečných funkcí (např. pokročilé litování, jednoduší vykreslování dat, atp.). Částice v konstantních elektrických a magnetických polích Související rovnice Řešení soustavy ODR v Matlabu Implementace Cvičení Van Allenovy radiační pásy Související rovnice Implementace 3. Pohyb částic v el. a mag. polích 3.1 Částice v konstantních elektrických a magnetických polích První úloha, kterou si implementujeme, se týká pohybu nabité částice v elektrickém a magnetickém poli. Tato úloha má dobře známé analytické řešení, ale my ji zde budeme řešit numericky právě proto, abychom demonstrovali, jak se liší analytické řešení obyčejných diferenciálních rovnic (ODR) na papíře od toho numerického. První program, se kterým budeme pracovat bude bohatě komentován a vysvětlen a je tedy důležité, abyste pochopili, co dělá každý řádek kódu. Pohyb částic ve složitých a často neuniformních či časově závislých elektrických a magnetických polích hraje důležitou roli ve fyzice plazmatu. Aplikace, ve kterých je třeba dobře znát pohyb částic jsou například: • Hmotnostní spektrometry, ve kterých jsou částice filtrovány podle poměru hmotnost/náboj. • Fúzní reaktory využívající magnetické pole pro udržení plazmatu [aa ] • Elektronové mikroskopy, kde je elektronový svazek směrován elektrickým polem [FEI10] • Hallův motor používaný pro vesmírný pohon (vizte obrázek 3.1) • a mnoho dalších... 3.1.1 Související rovnice Program vyvinutý v tomto oddíle řeší pohybové rovnice ve trojrozměrných kartézských souřadnicích pro částici s určitou hmotností a nábojem na kterou působí Lorentzova síla. FL = í(E + vxB) = 9(E + řxB) (3.1) kde q je náboj částice, E je elektrické pole, B je magnetické pole, v je rychlost částice a r je její poloha. Pohybová rovnice tedy nabývá tvaru f = -(E + řxB) (3.2) m kde m je hmotnost částice. Pohybová rovnice v tomto tvaruje užitečná při mnohých teoretických úvahách a při hledání analytického či přibližného analytického řešení této úlohy. Pokud ale chceme tuto rovnici řešit numericky, budeme se na ni muset podívat z trochu jiného úhlu. V první řadě musíme zdůraznit, že Matlab (stejně jako většina dalších numerických knihoven) neobsahuje funkce pro přímé řešení ODR druhého řádu, obsahuje ale velmi pokročilé a stabilní 16 Chapter 3. Pohyb částic v el. a mag. polích (a) Schematický pohled, zdroj: David Staack (b) Hallův motor v provozu Figuře 3.1: Hallovy motory využívají elektrické a magnetické pole k urychlení iontů ven z plazmatu. V roce 2014 se již Hallovy motory běžně využívají pro geostacionární družice i extraterestriální mise [GK08, page 15] funkce pro numerické řešení ODR prvního řádu i jejich soustav. Použijeme-li standardní značení, r = (x,y,z), (3.3) ř = {vx,vy,vz), (3.4) B = (Bx,By,Bz), (3.5) Ě = (Ex,Ey,Ez) (3.6) je poměrně jasné, že rovnice (3.2) je ekvivalentní následujícím šesti rovnicím. x = vx, (3.7) ý = vy, (3.8) ž = vz, (3.9) vx = -(Ex + vyBz-vzBy), (3.10) m vy = ^-(Ey + vzBx-vxBz), (3.11) m vz = — (Ez + vxBy - vyBx). (3.12) Jak vidíte, jednoduše jsme transformovali tři ODR druhého řádu popsané rovnicí (3.2) do šesti ODR prvního řádu. Transformování N diferenciálních rovnic druhého řádu do 2N diferenciálních rovnic prvního řáduje ve fyzice poměrně časté. Většina z Vás se s tím jistě setkala v kurzech teo- retické mechaniky, kde se mimo jiné provádí přechod mezi Lagrangeovským a Hamiltonovským formalismem. Z transformace provedené výše je zřejmé, že formalismem vhodnějším pro numerickou analýzu je ten Hamiltonovský. 3.1.2 Řešení soustavy ODR v Matlabu matlabová funkce, kterou budeme používat pro řešení soustavy ODR je pojmenovaná ode45. Tato funkce umí řešit diferenciální rovnice prvního řádu zadané ve tvaru q = /(',q)- (3.13) 3.1 Částice v konstantních elektrických a magnetických polích 17 kde t je nezávislá proměnná (v našem případě čas) a q = q(t) je vektor záviských proměnných. Nyní by už mělo být zřejmé, proč jsme přeformulovali pohybovou rovnici do tvaru popsaného rovnicemi (3.7) až (3.12). Pokud tyto rovnice porovnáte s rovnicí (3.13) vidíte, že v našem případě je q = (x,y,z,vx,vy,vz). V Matlabu potom funkci ode45 voláte následujícím způsobem (3.14) [t,q] = solver(@f,tspan,qO) Prvním argumentem (vstupní proměnnou) je odkaz na funkci f(t,q), která dává hodnoty pravých stran rovnic v naší soustavě. Druhým argumentem je interval nezávislé proměnné (u nás času), na kterém budeme systém ODR řešit a třetím je vektor počátečních podmínek pro q. Pokud zavoláte funkci ode45 tak, jak je uvedeno výše, dostanete jako výstup vektor t = tspan a matici q v následujícím tvaru ( 1 \ ( x y z Vx Vy v A 0.0 0 0 0 100.0 0 0 0.0001 -0.0000 0 99.91 -3.95 0 0.0002 0.0005 -0.0000 0 99.65 -7.91 0 0.0003 q — 0.0007 -0.0000 0 99.21 -11.84 0 0.0004 0.0010 -0.0001 0 98.60 -15.76 0 0.0005 0.0012 -0.0001 0 97.82 -19.64 0 V - > v ... •••/ (3.15) Sloupcový vektor t obsahuje hodnoty naší nezávislé proměnné (času) zatímco sloupce matice q obsahují rychlosti a pozice (tedy závislé proměnné) v odpovídajících časech. 3.1.3 Implementace jak už možná tušíte, budeme pro řešení pohybové rovnice v Matlabu potřebovat dva soubory. První z nich, odef un .m obsahuje funkci, která dává pravé strany našeho systému obyčejných diferenciálních rovnic a druhý, solve. m obsahuje počáteční podmínky pro řešič ode45, časový interval, na kterém se bude soustava ODR řešit a příkazy pro vykreslování. Oba soubory musí být umístěny v tomtéž adresáři a celý program se spouští přes soubor solve .m. Podívejme se nejprve na soubor odef un .m. Program 1: odefun.m function dqdt = odefun(t,q), i °/0 Toto je vstupní funkce pro 2 °/0 resic diferenciálních rovnic, ode45 . 3 °/0 Vstup: nezávislá proměnna t(cas), 6-rozmerny vektor 4 q °/o q = [x ,y , z , v_x , v_y , v_z] 5 °/0 vystup: casova derivace q, dqdt = dq/dt 6 7 qe = -1.602e-19;% Elementárni naboj [C] 8 m = 9.109e-31; °/0 Hmotnost častice [kg] 9 10 °/0 Nyni si definujeme proměnné popisujici polohu a n rychlost 18 Chapter 3. Pohyb částic v el. a mag. polích % Pomoci komponent vektoru q. 12 % Toto neni nezbytný krok, ale pouzivat napriklad 13 °/0 "x" namisto q(l), atd . . cini kod nize 14 % vice citelným a intuitivnim 15 x = q( 1) ; 16 y = q(2) ; 17 z = q(3) ; 18 vx = q(4) ; 19 vy = q(5) ; 20 vz = q (6) ; 21 22 Ex = 0 ; 23 Ey = 2e-6; 24 Ez = 0 ; 25 Bx = 0 ; 26 By = 0 ; 27 Bz = le-7; 28 29 % Nyni spočteme komponenty 6rozmerneho vektoru dqdt 30 % vector dq/dt . 31 dxdt = vx; % dx/dt = v_x 32 dydt = vy; °/0 dy/dt = v_y 33 dzdt = vz; % dz/dt = v_z 34 dvxdt = qe/m*(Ex + vy*Bz - vz*Bx) ; °/0 dv_x/dt 35 dvydt = qe/m*(Ey + vz*Bx - vx*Bz); °/0 dv_y/dt 36 dvzdt = qe/m*(Ez + vx*By - vy*Bx) ; °/0 dv_z/dt 37 % Nakonec sestavime samotný vektor dqdt 38 % z jeho komponent. Tim, ze vektoru dqdt 39 % priradime hodnotu jsme nastavili vystup funkce 40 odefun dqdt = [dxdt; dydt; dzdt; dvxdt; dvydt; dvzdt]; 41 42 end % konec funkce 43 I když je tento soubor obsahující funkci odef un k poměrně bohatě komentován, zaměřme se ještě jednou na jeho nej důležitější části. Jak jste se již naučili v první kapitole, musí každý soubor obsahující funkci začít definicí této funkce, tedy specifikací vstupních a výstupních proměnných. V našem případě máme na vstupu nezávislou proměnnou t a vektor závislých proměnných q a funkce vrací derivaci tohoto vektoru q (v ASCII zápisu dqdt) Na řádcích 8 a 9 jsou definovány elementární náboj a na řádcích 16 až 22 jsou přejmenovány elementy vektoru t tak, aby se nám s nimi dále lépe pracovalo. Toto přejmenování samozřejmě není nutné a mohli bychom dále v kódu místo x psát q(l) atd. na řádcích 35 až 37, kde vyjadřujeme právě derivaci jednotlivých komponent q. Na řádcích 23 až 28 potom definujeme komponenty elektrického a magnetického pole. Definovat konstanty qO a m v těle funkce není z programátorského hlediska optimální, pro lepší čitelnost kódu si to ale můžeme dovolit. Problém spočívá v tom, že během řešení soustavy ODR bude funkce odef un() zavolána mnohokrát a opakovaná definice též proměnné zabírá zbytečně systémové prostředky a navíc nemáte konstanty definovány všechny na jednom místě. Dále v tomto textu se naučíme, jak běh programu optimalizovat pomocí globálních proměnných. 3.1 Částice v konstantních elektrických a magnetických polích 19 Nyní, když jsme definovali funkci dávající pravou stranu systému obyčejných diferenciálních rovnic (3.13), můžeme konečně přistoupit k jejímu řešení a vykreslení výsledků. Potřebné příklady uložíme do souboru solve .m. Program 1: solve.m % Těmto Matlabov skript e pohybovou rovnici pro i % nabitou štici v p tomnosit elektrick ho 2 % a magnetiek ho pole. 3 4 % Nejprve mus me definovat po t e n podm nky 5 °/0 Pro p ipomenut , q = [x , y , z , v_x , v_y , v_z] 6 qO = [0;0 ; 0;le2 ; 0;0]; i 8 % Nyn definujeme asov interval, na kter m 9 % budeme eit na i soustavu ODR. 10 ti = 0; n tf = 2.5e-3; 12 N = 1000; 13 timespan = linspace (ti , tf , N) ; 14 °/0 funkce linspace(ti, tf , N) vytvo vektor 15 °/0 s line rn rostouc i hodnotami od "ti" 16 */. po "tf" s "N" kroky/prvky. 17 18 '/oNynspustme ei samotn. 19 % Z pis Oodefun k matlabu , e prvn m 20 % argumentem funkce ode45 je j in funkce. 21 [t, q] = ode45(@odefun , timespan, qO) ; 22 °/0 V sledkem bude lxN vektor t a 6xN matice q 23 24 % N sleduj c p kazy zobraz spo tanou trajektorii 25 mesh([q(: ,1) q(:,l)], [q(:,2) q(:,2)], [q(:,3) q(:,3)], [t ( :) 26 t(:)],'EdgeColor', 'interp', 'FaceColor', 'none'); 1 vykre slen view(2) °/0 nastav pohled shora 21 set (gca , ' FontSize ' , 16 , ' f ontWeight ' , ' bold ') °/» nastav 28 velikost p srna xlabel('x [m] ') °/0 popisek osy x 29 ylabel('y [m] ') °/0 popisek osy y 30 zlabel('z [m] ') °/0 popisek osy z 31 title(['Particle trajectory, ti=', num2str(ti), ' s, tf=', 32 num2str(tf), ' s' ]) °/0 nadpis grafu 33 print -depsc '0utput.eps' °/0 uloženi obrázku do souboru 34 3.1.4 Cvičení Cvičení 3.1 Zkopírujte si zdrojové kódy z předchozího cvičení do počítače nebo použijte m-soubory v adresáři Progl_Particle_motion. • Spusťte program a ověřte, že jste dostali křivku podobnou té na obrázku 3.2. 20 Chapter 3. Pohyb částic v el. a mag. polích • Jaký drift zde pozorujeme? • Jaký je směr driftové rychlosti elektronu a pozitronu? Particle trajectory, ti=0 s, tf=0.0025 s ^0.004 0.002 0 ......Ĺ... \ 1 \ í \ 1 ' \ 1 1 1 \ ] / \ \ \ / \ / \ r \ \ \ / / / \ J y :' \ < 0 0.01 0.02 0.03 0.04 0.05 Figure 3.2: Trajektorie elektronu v E = (0,1(T6,0) V/m, B = (0,0,1(T7)T aq(0) = (0,0,0,100m/s,0,0,0) Cvičení 3.2 Změňte náboj qe a hmotnost m tak, aby odpovídaly protonu to those of a proton. • Kolikrát musíte zvětšit přibližně časovou škálu, aby jste mohli u protonu pozorovat trajektorii charakteristickou pro tento drift? • Porovnejte amplitudy periodického pohybu elektronu a protonů a velikost jejich driftových rychlostí (data vyčtěte z grafu nebo z numerického řešení daného vektorem t a maticí q). Cvičení 3.3 Upravte program ze cvičení 3.1 tak, aby byla pohybová rovnice řešena bez přítomnosti elektrického pole s magnetickým polem kvadraticky rostoucím ve směru osy y E = (0,0,0) (3.16) B = (0,0,5J (3.17) \ 2 y Bz=B0^yJ (3.18) qo = (0,0,0,V2v0,V2v0,0) (3.19) Zvolte konstantu y q a časový interval (u,tf) tak, aby částice za tento čas urazila zhruba vzdálenost vo- • Pozorujete nějaký drift? Pokud ano, pro jaké hodnoty Bq, yo a vo jej pozorujete? • Jaký je směr driftové rychlosti pro elektron a pro pozitron? Pokročilé cvičení 3.1 Upravte program tak, aby se elektrické pole měnilo harmonicky v 3.2 Van Allenovy radiační pásy 21 čase, například Ex = E0-cos (cot), (3.20) Ey = 0, (3.21) Ez = 0, (3.22) a nastavte počáteční rychlost na v* = 0, (3.23) Vy=VQ, (3.24) vz = 0. (3.25) Nyní pro několik frekvencí co (např. 106, 107, 108 a 109 Hz) pozorujte pohyb protonu a elektronu. • Jak částice reagují na pole? • Porovnejte, jak moc jsou elektron a proton ovlivněny elektrickým polem pro různé frekvence. 3.2 Van Allenovy radiační pásy Objevení van Allenových radiačních pásů se často považuje za první velký objev v historii vesmírných letů. Jejich existence byla předpovězena Jamesem van Allenem a potvrzena družicí Explorer I v roce 1958. Pásy jsou přímým důsledkem magnetického pole Země, které zachycuje částice s vysokou energií a zabrání jim tak aby dopadly na povrch. Žádné stínění samozřejmě není ideální, takže vždy existuje určitý slabý tok částic směrem k povrchu. Tyto částice ionizují a excitují molekuly v nižších vrstvách atmosféry, což pozorujeme jako polární záři (zobrazena v úvodu této kapitoly). Existují dva stabilní van Allenovy pásy. První z nich je umístěn ve výšce cca 1000 to 6 000 km a druhý cca 13000 to 60000 km [Bakl2]. Třetí van Allenův pás byl objeven velmi nedávno a zdá se, že s postupem času mizí a objevuje se [Scil3; Shp+13]. Vzhledem k tomu, že nyní již máte povědomí o všemožných driftech, je formulace, že částice jsou zachyceny v magnetickém poli příliš vágní. V tomto oddíle upravíme Program 1, který jsme použili dříve, a vyřešíme pomocí něj pohyb nabité částice v magnetickém poli Země. V průběhu se naučíme, co jsou globální proměnné a jak je použít v Matlabu. Také se seznámíte s několika pokročilými příkazy pro vykreslování. Tento oddíl je inspirován přednáškami S. Markidise [Marl3].. 3.2.1 Související rovnice Řešené rovnice i metoda řešení jsou identické jako v předchozím oddíle 3.1. Hlavním rozdílem zde bude, že elektrické pole považujeme všude za nulové a magnetické pole bude silně záviset na poloze a budeme jej aproximovat magnetickým polem dipólu s konstantním magnetickým dipólovým momentem m. Považovat geomagnetické pole za čistě dipólové je poměrně hrubá aproximace, které s rozumnou přesností platí pouze do vzdálenosti několika Re (poloměrů Země). Dále od země je geomagnetické pole silně nesymetrické kvůli jeho interakci se solárním větrem. Většina experimentálního i teoretického výzkumu geomagnetického pole, je podporována NASA v rámci programu Living with a star (lws.gsfc.nasa.gov). Pokročilé modely geomagnetického pole a vesmírného počasí potom naleznete například na ccmc.gsfc.nasa.gov. 22 Chapter 3. Pohyb částic v el. a mag. polích Ve vektorovém zápisu je magnetické pole dipólu dáno výrazem »«=£(^-") kde m je magnetický moment, ;Uo permeabilita vakua a r = {x,y,z), (3.27 r=|r|. (3.28 Cvičení 3.4 Budeme pracovat v kartézské soustavě souřadnic se středem ve středu Země a osa z splývá se zemskou osou. V tomto systému je magnetický moment planety • Vyjádřete komponenty magnetického pole B = (Bx,By,Bz) v kartézských souřadnicích. • Jaká je přibližně hodnota M pokud na rovníku naměříme mag neticképole 3.12- 1(T5 T? ■ 3.2.2 Implementace Stejně jako v předchozím případě musíme definovat funkci odef un(), která nám bude dávat pravou stranu naší soustavy ODR. Tato funkce bude skoro stejná jako v předchozím případě. Program 2: odefun.m function dqdt = odefun(t,q) , i °/0 Toto je vstupní funkce pro 2 °/0 resic diferenciálních rovnic, ode45 . 3 °/0 Vstup: nezávislá proměnna t(cas), 6-rozmerny vektor 4 q °/o q = [x ,y , z , v_x , v_y , v_z] 5 °/0 vystup: casova derivace q, dqdt = dq/dt 6 7 °/0 V soubory solve.m, definujeme několik globalnich 8 proměnných. Zde řekneme °/0 Matlabu , chceme ve funkci použit proměnné "m" a "qe 9 " definovane v jiném souboru globál m; globál qe ; 10 n °/0 Přeznačeni poloh a rychlosti . 12 x = q( 1) ; 13 y = q(2) ; 14 z = q(3) ; 15 vx = q(4) ; 16 vy = q(5) ; 17 vz = q(6) ; 18 19 °/0 Elektrické pole bude nulové 20 Ex = 0 ; 21 Ey = 0 ; 22 Ez = 0 ; 23 3.2 Van Allenovy radiační pásy 23 24 °/0 Magnetické pole je nyni pocitano pomoci jine 25 externi funkce, bfield(r), která je funkci pouze polohy a je uložena v souboru bfield.m rO = [q(D, q(2) , q(3)]; °/. vektor polohy 26 B = bfield(rO); °/0 funkce vraci vektor mag. pole 27 % ... a my opet pro názornost komponenty preznacime. 28 Bx = B (1) ; 29 By = B (2) ; 30 Bz = B (3) ; 31 32 % Odsud dale je vse stejné jako v Programu 1 33 dxdt = vx; % dx/dt = v_x 34 dydt = vy; °/0 dy/dt = v_y 35 dzdt = vz; % dz/dt = v_z 36 dvxdt = qe/m*(Ex + vy*Bz - vz*By) ; °/0 dv_x/dt 37 dvydt = qe/m*(Ey + vz*Bx - vx*Bz); °/0 dv_y/dt 38 dvzdt = qe/m*(Ez + vx*By - vy*Bx) ; °/0 dv_z/dt 39 dqdt = [dxdt; dydt; dzdt; dvxdt; dvydt; dvzdt]; 40 41 end % konec funkce 42 První rozdíl oproti předchozímu programu vidíme na řádku 10. Místo toho, abychom proměnným m a qe přiřadili nějakou určitou hodnotu, prohlásili jsme je pouze za globální proměnné. Pokud vám myšlenka globálních a lokálních proměnných není známá, podívejte se na následující poznámku. Globální proměnná je proměnní, které je dostupná ve všech funkcích a programech, ve kterých byla prohlášena za globální. Běžné proměnné jsou naopak lokální a nejsou automaticky sdíleny mezi funkcemi a programy. Představte si situaci, kdy chcete použít proměnnou m jak v solve .m tak v odef un .m • Pokud použijete lokální proměnnou, musíte jí nastavit číselnou hodnotu v obou těchto souborech tím, že do obou napíšete m = 1.627e-27;. • Pokud ale použijete globální proměnnou, nastavíte numerickou hodno tu m = 1.627e-27; v souboru solve.m, který se spouští jako první, a deklarujete, že proměnná je globální tím, že napíšete globál m; jak do odef un .m tak do solve .m. Tím pádem je číselná hodnota uvedena jen na jednom místě. Výhoda druhého přístupu je zjevná. Pokud budete kdykoliv v budoucnu potřebovat hmotnost změnit na jinou hodnotu, stačí je přepsat na jednom místě. Druhý rozdíl v souboru odef un. m naleznete kolem řádku 26. Jak vidíte, magnetické pole zde není zadáno předpisem, ale pomocí další externí funkce uložené v souboru bfield.m, která rovněž používá některé z globálních proměnných. Níže vidíte neúplný zdrojový kód této funkce, do kterého musíte doplnit správné výrazy pro Bx, By, Bz (cvičení 3.4). Program 2: bfield.m function B = bfield(rO), i °/0 Tato funkce spocita geomagnetické pole B na pozici 2 rO. Pracuje s kartézskými souřadnicemi. % Funkce pouziva nasledujici globálni proměnné (vizte 3 solve.m): 24 Chapter 3. Pohyb částic v el. a mag. polích global Re; global q; global m; global muO; global M; 4 5 % opet přeznačeni 6 x = rO(l) ; 7 y = rO(2) ; 8 z = rO(3) ; 9 r = sqrt (x~2 + y~2 + z~2) ; io Ii % Doplňte nasledujici radky na zaklade odpovidajiciho 12 cvičeni . Bx = XXX; 13 By = YYY; 14 Bz = ZZZ; 15 16 B = [Bx; By; Bz] ; 17 end 18 Na závěr potřebujeme soubor se skriptem, obsahujícím všechny nezbytné proměnné, počáteční podmínky a vykreslovací příkazy. Program 2: solve.m °/0 Tento skript resi pohybovou rovnici 1 % pro nabitou častici v magnetickém poli Zeme 2 3 % Abychom program zrychlili, deklarujeme všechny nase 4 °/0 konstanty jako globálni proměnné 5 globál Re; globál m; globál qe; globál muO; globál M; 6 Re = 6378137; °/0 Polomer zeme v metrech 7 m = 1.627e-27; °/0 Hmotnost častice v kg 8 qe = 1.602e-19; °/» Naboj častice v C 9 esc =1; X Skalovaci faktor pro vykreslováni 10 c = 299792458; °/. Rychlost svetla in m/s 11 mu0 = 4*pi*le-7; °/0 Permeabilita vakua, V*s/(A*m) 12 M = 7.94e22; °/0 Magneticky moment Zeme, H/m 13 14 % P o t e n podm nky pro polohu ... 15 % ... definujeme nejprve ve sférických souradnicich 16 rO = 3*Re; 17 phiO = 0; 18 thetaO = pi/2; 19 °/0 ... a transformujeme do kartézských 20 xO = rO*sin(thetaO)* cos(pni 0) 21 yO = r0*sin(theta0)*sin(phi0) 22 zO = r0*cos (thetaO) 23 24 % Pocatecni podminky pro rychlost 25 % Definovane pomoci kin. energie a 2 uhlu 26 Ek_eV = 5e7; °/» energie [eV] 27 Ek = Ek_eV*l . 602e-19; % energie [J] 28 3.2 Van Allenovy radiační pásy 25 v.rO = c*(l+m*c~2/Ek)~-0.5 % velikost rychlosti ( 29 relativisticky) [m/s] v_phiO = 0; % azimut [rad] 30 v_theta0 = pi/4; °/0 polárni uhel [rad] 31 32 % Transformace souřadnic pro rychlost 33 vxO = v_r 0* sin (v_thetaO) * cos (v_phi 0) 34 vyO = v_r 0* sin (v_thetaO) * sin (v_phi 0) 35 vzO = v_r0*cos (v_theta0) 36 37 °/0 Vektor pocatecnich podminek 38 qO = [xO ; yO ; zO ; vxO ; vyO ; vzO] ; 39 40 °/0 Definováni časového intervalu 41 ti = 0; % pocatecni cas [s] 42 tf = 20; % konečný cas [s] 43 N = 10000; °/. počet kroku 44 timespan = linspace (ti , tf , N) ; 45 46 % Reseni systému ODR 47 % Nezapomeňte, ze odefun(t.q) se nyni lisi od Prgramu 1 48 [t, q] = ode45(Oodefun , timespan, qO) ; 49 50 51 °/0 Po vyřešeni rovnice je treba data vykreslit 52 % Zde nemusite nutne rozumět, co dela který 53 % z prikazu nize 54 55 h = figuře ; °/0 Vytvorime novy graf 56 57 % Nasledujici prikazy vykresli magnetické siločáry. [Upraveno 58 z MagLForce skriptu od A. Abokhodaira] n = 4; 59 d2r = pi/180; r2d=l/d2r; 60 tht=d2r*(0:5:360) ' ; 61 phi=d2r*(0:ceil(180/n):180); 62 hh=phi*r2d; 63 A = r0; 64 r=A*sin(tht) .~2; 65 rho=r . * sin (tht) ; 66 x=rho*cos(phi); y = rho * sin (phi ) ; 67 [nR , nC] =size (x) ; 68 u=ones (1 ,nC) ; z=r . * cos (tht) *u ; 69 plot3(x/Re,y/Re,z/Re,'r','LineWidth ' ,1 .0) ; 70 71 hold on; °/0 Tento prikaz rika Matlabu , ze cokoliv budeme dale 72 vykreslovat ma byt součásti tehoz grafu!!! 73 % Vykresleni elipsoidu predstavujiciho Zemi. 74 26 Chapter 3. Pohyb částic v el. a mag. polích [u,v,w]=sphere(30); surf(esc*u, esc*v, esc*0.9*w); colormap('default'); camlight right ; lighting phong; axis equal °/0 Nastavi stejné meritko os . °/0 Vykresleni trajektorie, přidáni popisku os plot3(q(: ,1)/Re , q(:,2)/Re, q(: ,3)/Re,'LineWidth ' , 1 . 0) set (gca , ' FontSize ' , 18) °/0 velikost pisma xlabel('x [R_e]') °/» popisek osy x ylabel('y [R_e]') °/» popisek osy y zlabel('z [R_e]') °/» popisek osy z title(['Proton trajectory, E = 50 MeV, t_i=', num2str(ti) , s, t_f=', num2 str (tf) , ' s ' ] ) °/0 nadpis grafu print -dpng -r200 ' ProtonMagField . png' °/0 ulozi obrázek do souboru ve formátu png s rozlisenim 200 dpi. Obsah souboru solve .m je v principu také podobný tomu předchozímu. Na řádcích 6 až 13 definujeme globální proměnné. Konstanty, které budeme chtít používat jako globální proměnné je nejlepší definovat a přidělit jim číselné hodnoty na začátku prvního souboru, který bude spuštěn. Na řádcích 15 až 39 definujeme počáteční podmínky pro rychlost a polohu. Pro lepší představitelnost definujeme počáteční podmínky nejprve ve sférických souřadnicích a následně je transformujeme do kartézských. Příkaz pro řešení systému ODR je na řádku 49 a od řádku 57 dále jsou vypsány příkazy pro vykreslování. Všimněte si zejména příkazu hold on, který Matlabu říká, že všechny další objekty se budou vykreslovat do téhož grafu. Také si povšimněte, že před vykreslením dat všechny proměnné škálujeme poloměrem Země Re. Pokud jste chybějící výrazy pro mag. pole v souboru bf ield .m uvedli správně, program by měl pro předdefinované počáteční podmínky vykreslit graf podobný tomu na obrázku 3.3. Proton trajectory, E = 50 MeV, ti=0 s, tf=20 s Figuře 3.3: Pohyb protonu s vysokou energií v geomagnetickém poli 3.2 Van Allenovy radiační pásy 27 Cvičení 3.5 Popište pohyb vysoce energetického protonu v geomagnetickém poli, jaké tři komponenty pohybu můžeme pozorovat? Podívejte se na pohyb částice v rovinách xy případně y z. ■ Cvičení 3.6 Přestože je model, který jsme vyvinuli, poměrně jednoduchý, může nám pomoci v pochopení základní struktury van Allenových radiačních pásů. Zkuste změnit počáteční polohu protonu. • Jaká je maximální počáteční vzdálenost od Země, pro kterou má proton ještě stabilní dráhu? • Jaká je minimální počáteční vzdálenost od Zeme, pro kterou proton ještě nedosáhne povrchu? Pokročilé cvičení 3.2 Nahraďte proton v předchozím cvičení elektronem a zkuste najít počáteční podmínky (zejména vzdálenost od Země), pro které bude elektron mít stabilní dráhu. Očekáváte, že to bude dále nebo blíže povrchu? Jak se mění drift elektronu oproti protonu? ■ 4.1 Úvod Znalosti Matlabu, které získáte v této kapitole se Vám budou hodit nejen ve fyzice plazmatu. Kapitola se zabývá zpracováním dat včetně interpolace a integrace diskrétních dat. Kromě toho se naučíte jak používat smyčky a podmínky, jejichž znalost je velmi užitečná při hromadném zpracování velkých souborů dat. 4.2 Smyčky a podmínky v Matlabu Podívejme se nejprve na dva velmi důležité prvky programování (nejen) v Matlabu, smyčky a podmínky. Pokud již máte předchozí zkušenosti s programováním, bude pro Vás tento oddíl poměrně triviální. Pokud jste se se smyčkami a podmínkami v programování zatím nesetkali, bude se Vám jejich znalost jistě hodit ve vaší vědecké kariéře. Smyčky jsou programátorské struktury, které Vám umožní mnohokrát opakovat určitou sekvenci příkazů. Ve fyzice se může jednat například o • vykreslování mnoha datových souborů v témž formátu, • hromadný převod velkého objemu dat, • řešení diferenciálních rovnic pro několik různých počátečních/okrajových podmínek. • ... V Matlabu, stejně jako v mnoha ostatních programovacích jazycích, existují dva druhy smyček, smyčky typu for a smyčky typu while. Smyčku typu for použijete, pokud chcete danou sekvenci příkazů spustit Af-krát (například "zintegruj data v prvních 50 souborech"), zatímco smyčka typu while se bude spouštět neustále dokola, dokud se nesplní určitá podmínka (např. "integruj data, dokud nenajdeš integrál větší než 10"). Ve fyzice je nicméně použití smyčky typu while poměrně omezené a navíc je v případě nepozornosti možné vytvořit nekonečnou smyčku while. Z toho důvodu se zde zaměříme pouze na smyčku typu/or, která by Vám ve většině fyzikálních aplikací měla stačit. Podmínky se používají, má-li se určitá sekvence příkazů spustit jen tehdy, když je určitý výrok pravdivý a obzvláště se Vám budou hodit v kombinaci se smyčkami. Příkladem použití může být například hromadné vykreslení dat, která se vykreslí v logaritmické škále pokud je rozdíl mezi největší a nejmenší hodnotou alespoň dva řády. 30 Chapter 4. Zpracování dat 4.2.1 Smyčka tor Jak už bylo zmíněno, smyčky typu for Vám umožní JV-krát provést určitou operaci. Pokud má smyčka proběhnout desetkrát, použijete ji následujícím způsobem for j=l:10, i °/0 vaše príkazy 2 end 3 Proměnná j je proměnnou cyklu, která postupně nabývá hodnot od 1 do 10 a můžete ji používat uvnitř cyklu. Ukažme si nyní použití cyklů na konkrétním programu, konkrétně vykreslí funkci y = xn pro n G (1,6) a x G (0,1). figuře ; °/0 vytvorí grad i colors = ['r', 'g', 'b', 'm', 'k', 'y', 'C, 'r']; °/„ toto je 2 vektor textových proměnných for j=l:8, % will run for 8 values of j 3 x = linspace (0,1 ,50) ; °/0 vektor délky 50 s komponentami od 0 4 do 1 color = colors (j); °/0 vybere j-ty prvek vektoru colors 5 plot(x, x.~j, color); 6 hold on; °/0 chceme vse v jednom grafu i end 8 Cvičení 4.1 Upravte program výše tak, že vykreslí y=—^ proxG (0,27T)a«G (1,5) (4.1) Nezapomeňte, že x je vektor a musíte tedy použít prvkové operace. ■ Ve fyzice plazmatu i v dalších odvětvích fyziky je často třeba pracovat s daty z literatury a často se stane, že data, která jsou k dispozici jsou v jiných jednotkách, než potřebujete. V dalším cvičení proto vytvoříte program, který bude hromadně konvertovat účinné průřezy zadané v jednotce cm2 do jednotky m2. Abyste tento program mohli napsat, budete potřebovat navíc funkce loadO adlmwriteQ. Funkce load() funguje velmi intuitivně, načte data ze souboru do Matlabové proměnné. Vstupní data mohou mít několik formátů (hodnoty oddělené tabulátorem, středníkem, csv, atd...) a do proměnné je načtete příkazem matrix = load('datafile.dat'); Naopak pokud chcete uložit data do ASCII souboru, použijete funkci dlmwrite (), která funguje podobně jednoduše. dlmwrite('new_datafile.dat', matrix, '\t'); Jak vidíte, bere si funkce tři parametry, název výstupního souboru, proměnnou, kterou chcete uložit a symbol oddělující sloupce dat. Běžně používané oddělovače sloupců jsou ',' (čárka), ';' (středník) a ' \ť (tabulátor). Poslední věc, kterou musíte znát jsou základy práce s řetězci v Matlabu. Řetězec je termín používaný v programování de facto pro textové proměnné a vzhledem k tomu, že řetězec v Matlabu je vlastně jenom vektor znaků s ním můžete pracovat podobně s tím rozdílem, že části řetězce musí být vždy uzavřeny v jednoduchých uvozovkách ('). 4.2 Smyčky a podmínky v Matlabu 31 Ukažme si spojování řetězců na příkladu. Následující skript otevře a vykreslí postupně obsah souborů cskl. dat až csk3. dat (tyto soubory můžete najít v adresáři Prog3_cross_section_convert). colors = ['r', 'g', 'b', 'm', 'k', 'y' , 'C, >r>]; i figure; 2 for j=l:3, 3 color = colors(j); 4 matrix = load(['csk', num2str(j), '. dat '] ) ; 5 xdata = matrix(:,1); °/0 prvni sloupec dat 6 ydata = matrix(:,2); °/0 druhy sloupec dat 7 plot(xdata, ydata, color); 8 hold on; 9 end 10 xlim ( [0 , le6] ) ; °/0 nastavi rozsah osy x 11 Data, která vám tento program vykreslí nejsou náhodné křivky, ale jsou to účinné průřezy srážek elektronů s argonovými atomy. Na ose x je elektronová teplota v jednotce Kelvin a na ose y účinný průřez v m2. Konkrétně soubor cskl. dat obsahuje data pro elastickou srážku e + Ar —> e + Ar, soubor csk2 .dat data pro excitaci argonu e + Ar-^e + Ar*, a soubor csk3. dat data pro ionizaci e + Ar->2e + Ar+. (4.2) (4-3) (4-4) Cvičení 4.2 Program uvedený výše vykreslí graf, který není vzhledem k velmi rozdílným magnitudám účinných průřezů příliš intuitivní. Upravte program tak, aby • byla elektronová teplota vyjádřena v jednotce eV, • a osa y byla logaritmická. Výsledný graf by měl vypadat podobně jako ten následující 10 20 30 40 Electron energy [eV] kde červená křivka odpovídá elastické srážce, zelená křivka excitaci a modá křivka ionizaci. • Vyčtěte z obrázku, jaký je ionizační a excitační práh argonu. Cvičení 4.3 Upravte program tak, že převede teplotu elektronů z Kelvin na elektronvolty a uloží konvertovaná data do souboru csevN. dat (kde N je odpovídající číslo souboru). 32 Chapter 4. Zpracování dat 4.2.2 Podmínky if... elif... else V předchozím oddíle byla popsána smyčka typu for, která umožňuje opakovaně provádět nějaký úkon. Dalším důležitým komponentem pokročilejších programů jsou podmínky. Použitím podmínky můžete dosáhnout toho, že pokud je určitý výrok pravdivý, provede se jiný příkaz než když je výrok nepravdivý. Předpokládejme například, že vyvíjíte určitý složitý program, ve kterém často potřebujete hledat inverzní matice. Jak jistě víte, inverzní matice A-1 k matici A je taková, že A^A = I (4.5) kde / je jednotková matice. Cvičení 4.4 Funkce, která za Vás v Matlabu spočítá inverzní matici se jmenuje inv() a jejím jediným argumentem je matice, kterou chcete invertovat. • Definujte 3x3 singulární matici A a pokuste seji invertovat funkcí inv(). Jak vypadá výsledná matice? ■ Jak jste viděli v předchozím cvičení a jak i jistě víte, snaha spočítat inverzní matici k singulární matici nevede k ničemu příliš užitečnému a pokud k něčemu takovému dojde ve vašem složitém hypotetickém programu, nebudete tušit, kde se stala chyba. Mnohem praktičtější by bylo definovat vlastní funkci, která inverzi provede jen když vstupní matice není singulární. Pokud matice singulární je, měla by Vaše funkce vypsat chybu a ukončit program. Funkce, kterou jsme právě popsali může vypadat následovně. function inverted = myinv(matrix), i rows = length(matrix(:,1)); °/0 spocita délku prvniho 2 sloupce = počet radku cols = length (matrix (1 ,:)) ; °/» spocita delku prvniho 3 radku = počet sloupců 4 if(rows ~= cols), % pokud se rows NEROVNÁ cols 5 disp('The supplied matrix is not a square 6 matrix') return °/0 ukonči skript 7 elseif (det (matrix) == 0); °/0 pokud je determinant nula 8 disp('The matrix is singular.') 9 return °/0 ukonči skript io else, % pokud jsou obe podm nky vyse nepravda n inverted = inv(matrix); °/0 provede inverzi a 12 vrati invertovanou matici end 13 end 14 Syntaxe podmínky if. . .elseif. . .else by měla být zjevná z kódu uvedeného výše. Vidíte, že výraz, jehož pravdivostní hodnota se ověřuje, je uveden v závorkách za if nebo elseif. Matlab postupuje procedurálně a podmínky testuje od první k poslední, tedy 1. Nejdříve zkontroluje tvar matice. Pokud není čtvercová, zobrazí se chyba. 2. Pokud je matice čtvercová, otestuje se její determinant. Pokud je nulový, zobrazí se chyba. 3. Pouze pokud je matice čtvercová a není singulární se provede inverze. 4.2 Smyčky a podmínky v Matlabu 33 Cvičení 4.5 Ověřte, že funkce myint () se chová správně, vyzkoušejte ji na několika maticích, které nejsou čtvercové a na maticích, které jsou singulární. Přidejte další elseif, které vypíše varování pokud je hodnost vstupní matice vyšší než 10. Vyzkoušejte. Pokročilé cvičení 4.1 Podívejte se do adresáře Prog4_cross_section_convert_if, který obsahuje tři soubory s účinnými průřezy Dva z těchto souborů obsahují teplotu elektronů v e V a jeden z nich v Kelvin. • Upravte program ze cvičení 4.2 přidáním smysluplné if. . .else podmínky, která rozhodne, zda by soubor měl být konvertován nebo ne. • Podobný jednoduchý program by byl velice užitečný například při vývoji kolizně -radiačního modelu plazmatu, kde je třeba hromadně předzpracovat tisíce podobných souborů. Jaká jsou omezení tohoto programu? Maxwell-Boltzmannovo rozdělení Střední rychlost Tvar Maxwell-Boltzmannova rozdělení Časový vývoj rozdělovači funkce Rychlostní koeficienty 5. Interakce částic v plazmatu Tato kapitola ukazuje rozdíly mezi rovnovážnými a nerovnovážnými rozdělovacími funkcemi a měla by Vám zprostředkovat to, jak závisí makroskopické proměnné na tvaru rozdělovači funkce. Na příkladech z fyziky plazmatu se zde naučíte základy numerického integrování a interpolace v Matlabu. 5.1 Maxwell-Boltzmannovo rozdělení 5.1.1 Střední rychlost Jak již víte z přednášky, rozdělení rychlosti se v rovnovážném stavu řídí Maxwell-Boltzmannovým rozdělením. Rozdělovači funkce velikosti F (v) rychlosti pro částice se třemi stupni volnosti potom je 9 / m \3/2 /— mv2\ F<,) = w(a*ř) exp(w> <51> Rozdělovači funkce, zejména ty elektronové, Vám toho o plazmatu mohou poměrně hodně prozradit, pokud víte, jak je interpretovat. Především je důležité pochopit, jaký vliv mají rozdělovači funkce na makroskopické proměnné, které vnímáme intuitivně. Vzhledem k tomu, že mnoho makroskopických veličin může být velmi často vyjádřena jako integrály rozdělovači funkce, měli byste se v tomto oddíle naučit, jak numericky integrovat v Matlabu. Některé z integrálů, které v tomto oddíle budeme počítat numericky mají také analytické řešení. Nicméně ve většině laboratorních výbojů platí, že rozdělovači funkce určitých typů částic (především elektronů) nelze popsat analytickým výrazem. Proto zde upřednostníme numerickou integraci před analytickou. Schopnost numericky integrovat data samozřejmě upotřebíte i v jiných odvětvích fyziky. Existuje několik způsobů, jak numericky integrovat data. Nejjednodušším z nich je lichoběžníkové pravidlo, které si pravděpodobně pamatujete z matematické analýzy. Při použití lichoběžníkového pravidla je funkce aproximována lichoběžníky a celková plocha těchto lichoběžníků je potom odhadem určitého integrálu mezi body a ab. Lichoběžníkové pravidlo je schematicky znázorněno na obrázku 5.1(a). Pokročilejší metodou integrace je takzvané Simpsonovo pravidlo, které funkci f(x), kterou chceme integrovat, aproximuje několika po 36 Chapter 5. Interakce částic v plazmatu částech spojitými polynomy Pi{x) (viz obrázek 5.1(b)). Proto dosahuje Simpsonova metoda vyšší přesnosti pro rychle se měnící funkce. (a) Lichoběžníkové pravidlo 9 9-l 9 2 (b) Simpsonovo pravidlo Figuře 5.1: Ilustrace lichoběžníkového pravidla pro numerickou integrace (vlevo) a Simpsonova pravidla (vpravo). V případě Simpsonova pravidla je horní strana lichoběžníků dána polynomem, který je analyticky integrovatelný. Z obrázku 5.1(a) je také zřejmé, že při snižování šířky lichoběžníků poroste přesnost této metody. Pro naše účely bude lichoběžníkové pravidlo mít dostatečnou přesnost. Algoritmus, který provádí integraci lichoběžníkovým pravidlem je v Matlabu implementován pod jménem trapzO Vstupní argumenty funkce jsou dva vektory, první z nich obsahuje nezávislou proměnnou x a druhý funkci, kterou chceme integrovat, f(x). x = linspace(0, pi, 100); fx = sin(x . ~2) ./log(x) ; I = trapz (x , fx) ; Jak vidíte, není u integrálu třeba udávat integrační meze. Je to proto, že integrační meze jsou dány minimální a maximální hodnotou vstupního vektoru x. Jak již chápete, tři řádky kódu výše spočítají následující určitý integrál sin r ■dx - sin r ■ dx (5.2) V navazujícím cvičení se naučíte používat funkci trapz () a otestujete její přesnost na příkladu, který má i analytické řešení. Cvičení 5.1 Program 5 uvedený níže je šablona pro vykreslení Maxwell-Boltzmannova rozdělení a pro spočtení středních rychlostí odpovídajících určité rozdělovači funkci. Program dokončete tím, že • doplníte výraz pro F (v) na řádku 7 (rozdělovači funkce velikosti rychlosti ve 3D, koncentraci částic zvolte rozumně), • spočítáte střední rychlost na řádku 12 a střední kvadratickou rychlost na řádku 13 (obojí pomocí lichoběžníkového pravidla), • doplníte analyticky spočítané výrazy pro střední rychlost, střední kvadratickou rychlost a nejpravděpodobnější rychlost na řádcích 16 až 18. Až dokončíte program, odpovězte na následující otázky • Odpovídají si numericky spočítané a analyticky spočítané hodnoty? 5.1 Maxwell-Boltzmannovo rozdělení 37 • Která ze tří rychlostí v grafu je nejnižší a která nejvyšší? Mění se jejich pořadí pokud měníte teplotu? • Zkuste změnit počet kroků ve funkci linspace (), například na 10, 50, 100, 500. Jak se při tom mění shoda numericky spočítaných hodnot s analytickými? Program 5: Maxwell-Boltzmann distribution and speeds mO = 1.6605e-27; °/0 a.m.u. v kilogramech i m = 40*m0; °/» hmotnost častice v kg 2 kB = 1.380e-23; % Boltzmannova konst. in m~2 kg s~-2 K~-l 3 T = 1000; X teplota v Kelvin 4 5 v = linspace (0 ,4000 ,500) ; °/. x data 6 Fv = . . .; X rozdělovači funkce i plot(v , Fv) ; 8 hold on; 9 10 % numerická integrace n v_mean = . . . ; '/, stredni rychlost 12 v_sq_mean = . . . ; °/» stredni kvadratická rychlost 13 14 % analytické výrazy 15 v_mean_an = . . . ; '/, stredni rychlost 16 v_sq_mean_an = . . . ; '/, stredni kvadratická rychlost 11 v_mp_an = ...; '/, ne j pravdepodobne j si rychlost 18 19 lineht = max (Fv) * 1 . 1 ; °/0 vyska car v grafu 20 plot ([v_mean , v.mean] , [0, lineht], 'r'); °/» stredni 21 rychlost (červena) plot ([ v_sq_mean , v_sq_mean] , [0, lineht], 'g'); °/» stredni 22 kvadratická rychlost (zelena) plot ([ v_mp_an , v_mp_an] , [0, lineht], 'm'); % 23 nej pravdepodobnej si rychlost (modra) hold off 24 Pokročilé cvičení 5.1 V předchozím programu jsme numericky spočítali střední rychlost a střední kvadratickou rychlost, ale nejpravděpodobnější rychlost byla spočtena analyticky. Abyste mohli spočítat nejpravděpodobnější rychlost analyticky, musíte najít maximum rozdělovači funkce (zadaná vektory v a Fv). Zeptejte se Internetu, jak to udělat. Jak víte z přednášek, rozdělovači funkce musí splňovat normovači podmínku / F{v)áv = n. (5.3) Jo Proto, pokud integrujete rozdělovači funkci od určité rychlosti vo do 00, dostanete počet částic, které mají velikost rychlosti větší nebo rovnu vo- Pokud integrujete F (v) od vo do v\, dostanete počet částic s velikostí rychlosti v intervalu (vo, vi). Tohoto se týká cvičení 5.2. 38 Chapter 5. Interakce částic v plazmatu Cvičení 5.2 Hmotnost dusíkových molekul je 28 a.m.u. a jejich koncentrace za standardních podmínek je cca 1.7 • 10 m . Kolik dusíkových molekul v místnosti, kde se právě nacházíte je rychlejších než 50, 500, 1000, 2500, 5000 a 10000 m/s? ■ 5.1.2 Tvar Maxwell-Boltzmannova rozdělení Pravděpodobně jste si všimli, že rovnovážné Maxwell-Boltzmannovo rozdělení je dáno jen dvěma makroskopickými parametry, konkrétně hmotností částice m a teplotou T. Následující cvičení nevyžaduje žádné programování, pokud ale náhodou nevíte odpověď, můžete použít program výše. Cvičení 5.3 Následující obrázek ukazuje tři rovnovážná rozdělení velikosti rychlosti • Pokud je hmotnost částic m stejná ve všech třech případech, která z rozdělovačích funkcí odpovídá nejnižší teplotě a která nejvyšší? • Pokud je konstantní teplota T, které rozdělení odpovídá částicím z nejnižší hmotností a které s nejvyšší? .x 10 500 1000 1500 2000 2500 Velocity magnitude, v [m/s] 3000 5.2 Časový vývoj rozdělovači funkce V předchozím oddíle jsme rozebrali klíčové vlastnosti Maxwell-Boltzmannova rozdělení. Nicméně tohoto rozdělení se dosahuje obecně pouze v termodynamické rovnováze a mimo ni může být rozdělovači funkce v podstatě libovolná. Navíc jsme zatím uvažovali pouze rozdělení velikosti rychlosti, čímž jsme předpokládali homogenní a izotropní rozdělovači funkci rychlosti. Zde budeme stále předpokládat, že prostor je homogenní (nezávisí na souřadnici r), ale rozdělovači funkce již nebude izotropní (tedy může záviset na vektoru rychlosti v, nejen na jeho velikosti). Program, pomocí kterého budeme analyzovat časový vývoj rozdělovači funkce není příliš komplikovaný. Je založen na nejjednodušší nerovnovážné formulaci Boltzmannovy kinetické rovnice, ve které je prostor homogenní, síla působící na částice je nulová F = 0, (5.4) a srážkový člen je vyjádřen v Krookově tvaru f-fo d_f\ ^Aoii (f-fo) (5.5) 5.2 Časový vývoj rozdělovači funkce 39 kde f q je rovnovážná rozdělovači funkce rychlosti,? je relaxační doba a vm je srážková frekvence. Celá Boltzmannova kinetická rovnice se tímto redukuje na ^|^ = -vm-(/(v,0-/o(v)) (5.6) a má následující analytické řešení f(y,t) =/0(v) + [/(v,0)-/0(v)]e-v™ř. (5.7) Pokud přijmeme tuto nejjednodušší aproximaci, není ani nutné řešit Boltzmannovu kinetickou rovnici a stačí jen vykreslit řešení. A právě toto dělá program, jehož zdrojový kód je uvedený níže. Kód je opět komentován, zaměřte se ale zejména na následující řádky. • Řádek 5 nastavuje časový interval, po který je řešení vykreslováno. • Řádek 6 nastavuje rozpětí rychlostí, podobně jako předchozí program. • Řádky 11 až 13 nastavují počáteční rozdělovači funkci rychlosti • Řádek 32 Udává škálu grafů. Pokud jej zakomentujete, budou mít grafy proměnlivou škálu a pokud bude odkomentován, bude škála stále stejná. Program 6: Vykreslení časového vývoje rozdělovači funkce. mO = 1.6605e-27; °/0 a.m.u. v kilogramech m = 40*m0; °/» hmotnost častice e v kg kB = 1.380e-23; °/0 Boltzmannova konstanta v m~2 kg s~-2 K~-l Vm = 5e8; °/0 srážková frekvence v Hz time = linspace(0, le-8, 100); °/0 casovy interval v = linspace (-2000 ,2000 ,200) ; °/. interval rychlosti % pocatecni rozdělovači funkce , f_0 sigma = 10; °/0 sirka vO = 500; % nej pravdepodobnej si rychlost fx_init = l/(sqrt(2*pi)*sigma)*exp(-(v-vO) ."2/(2*sigma~2)) ; fy_init = v*0; fz_init = v*0; °/0 Spočteni rovnovážne teploty M-B rozděleni (plyne ze Z.Z.E) v_sq = sqrt(trapz(v, (fx_init + fy_init + fz_init) . *v . ~2)) ; Teq = l/3*m*v_sq~2/kB; % rovnovážna rozdělovači funkce fx_eq = sqrt(m/(2*pi*kB*Teq)) * exp(-m*v . ~2/(2*kB*Teq)) fy_eq = sqrt(m/(2*pi*kB*Teq)) * exp(-m*v.~2/(2*kB*Teq)) fz_eq = sqrt(m/(2*pi*kB*Teq)) * exp(-m*v.~2/(2*kB*Teq)) fx = fx_init; fy = fy_init; fz = fz_init; hFig = figuře ; °/0 vytvorime graf , nastavime rozmery set(hFig, 'Position', [100 300 1000 300]); for j = 1 : length (time) , °/0 smyčka for dela kroky v case absmax = max([fx, fy, fz])*l.l; °/0 proměnlivá skala osy y 40 Chapter 5. Interakce částic v plazmatu % absmax = max([max(fx_init),max(fy_init),max(fz_init)]); °/0 pro konstantní skalu osy y odkomentujte subplot(1,3,1); % prvni podobrazek - v_x plot(v, fx, 'b', 'LineWidth', 2); ylabel('f(v_x) ') ; xlabel( ' v_x'); ylim([0, absmax]); whitebg('white') subplot (1 , 3 , 2) ; °/0 druhy podobrazek - v_y plot(v, fy, 'r', 'LineWidth', 2); ylabel('f(v_y) ') ; xlabel( ' v_y'); ylim( [0 , absmax]) subplot (1 , 3 , 3) ; °/0 treti podobrazek - v_z plot(v, fz, 'm', 'LineWidth', 2); ylabel('f(v_z) ') ; xlabel( ' v_z'); ylim([0, absmax]); % nastaveni zahlavi obrázku ax = axes('Units','Normal ', 'Position' , [.075 .075 .85 .85],' Visible','off ') ; set(get(ax,'Title '),'Visible' , 'on') title(['t = ', num2str(time(j)) , ' s']); °/0 spočteni nových hodnot fx, fy, fz fx = fx_eq + (fx_init - fx_eq)*exp(-time(j)*Vm); fy = fy_eq + (fy_init - fy_eq)*exp(-time(j)*Vm); fz = fz_eq + (fy_init - fz_eq)*exp(-time(j)*Vm); M(j) = getf rame (gcf) ; °/0 pripoji snimek k proměnné M end movie2avi(M, 'out.avi', 'fps', 7); °/0 ulozi animaci Cvičení 5.4 Spusťte program výše, prohlédněte si výstupní animaci a odpovězte na otázky. • Jak můžeme fyzikálně interpretovat počáteční podmínku pro rozdělovači funkci? • Jaký typ částic by mohla rozdělovači funkce popisovat? • Srážkovou frekvenci jsme nastavili na smysluplnou hodnotu 5 • 108 Hz. Jaký je řádově čas potřebný k dosažení rovnováhy? • Zkuste v programu 6 zvýšit nebo snížit srážkovou frekvenci. Jak se mění doba potřebná k dosažení rovnováhy? 5.2 Časový vývoj rozdělovači funkce 41 Obrázky níže ukazují možný výstup programu 6. Je třeba zdůraznit, že obrázky mají proměnlivé měřítko na ose y. 0.03 0.025 „ 0.02 á- 0.015 0.01 0.005 01— -2000 x 10" 6 t = 0 ns 0.03 0.025 0.02 0.015 0.01 0.005 2000 -2000 x 10" 0 v y 3 ns 6 ě 4 0.03 0.025 0.02 0.015 0.01 0.005 2000 01— -2000 x 10" 6 š" 4 2000 -2000 2000 -2000 2000 -2000 2000 x 10" 1.5 —x 1 0.5 -2000 1.5 x 10" 0.5 t = 6.5 ns x10~3 1.5 —>. 1 0.5 2000 -2000 1.5 t = 10 ns x10~3 0.5 x 10" 1.5 —M 1 > 0.5 2000 -2000 1.5 x 10" 0.5 2000 -2000 2000 -2000 2000 -2000 2000 42 Chapter 5. Interakce částic v plazmatu Rychlostní koeficienty Na závěr si vyzkoušíme, jak různé rozdělovači funkce rychlosti ovlivňují makroskopické proměnné, v našem případě rychlostní koeficient reakce. Rychlostní koeficient reakce spočítáme z jejího celkového účinného průřezu ar jako Ve vztahu výše již předpokládáme, že účinný průřez závisí pouze na rychlosti dopadající částice a nikoliv na úhlu srážky. Při počítání integrálu (5.8) narazíme najeden technický problém, konkrétně na interpolaci účinných srážkových průřezů <7r(v), protože účinné průřezy většiny reakcí nelze vyjádřit analyticky a jsou k dispozici pouze jako tabelovaná data. Vzhledem k tomu, že účinný průřez cr(v) není většinou tabelován pro mnoho hodnot v, je potřeba tabelovaná data interpolovat. V Matlabu lze interpolovat ID data pomocí funkce interpl (). Řekněme, že máte funkci zadanou vektory xdata (nezávislá proměnná) and ydata (funkční hodnoty) a chcete získat funkční hodnoty pro jemnější vektor nezávislé proměnné xinterp. Syntaxe potom bude následující xdata = [1 , 5 , 10 , 15, 20]; ydata = [1, 35, 110, 235, 410]; xinterp = linspace(8 ,17 , 0.5) ; yinterp = interp1(xdata, ydata, xinterp, 'pchip', 0); Jak vidíte, funkce interpl () má 5 vstupních parametrů. První dva z nich udávají data, která budete interpolovat, tedy vektor nezávislých proměnných (xdata) a stejně dlouhý vektor funkčních hodnot (ydata). Vektor (xinterp) je potom jiný vektor nezávislých proměnných (typicky s jemnějším krokem), pro které chcete funkci interpolovat. Čtvrtý vstupní parametr je interpolační metoda a pátý (nula) udává, že se neprovádí extrapolace a všechny body, které jsou mimo interval daný vektorem xdata budou považovány za nulové. Při použití funkce interpl () můžete zvolit z několika interpolačních metod. Nejpop-ulárnějšími volbami jsou 'nearest', u které jsou data interpolovaná na základě hodnoty nejbližšflio souseda, ' linear', kdy je použita lineární interpolace a ' pchip', u které se používá polynomiální interpolace. Zkusme si nyní spočítat rychlostní koeficient ionizace argonového atomu dopadem elektronu. Výpočet provedeme pro dvě různé rozdělovači funkce, pro Maxwell-Stefanovu rozdělovači funkci a pro skoro monoenergetický svazek částic, které tím pádem mají skoro stejnou rychlost. Program by neměl obsahovat nic, co jsme se do této chvíle nenaučili. Jeho výstupem jsou dvě čísla, rychlostní koeficient pro Maxwell-Boltzmannovo rozdělení a rychlostní koeficient pro zmíněný monoenergetický svazek. m = 9.109e-31; °/0 hmotnost elektronu kB = 1.380e-23; % Boltzmannova konst. v m~2 kg s~-2 K~-l TeV = 15; °/0 teplota elektronu v eV (5.8) e + Ar->2e + Ar+. (5.9) Program 7: Ionizace dopadem elektronu 5.3 Rychlostní koeficienty 43 T = TeV*11604; °/, teplota elektronu v K 4 5 % načteni účinného prurezu 6 sigdata = load('sigmaion.dat'); 7 xsig = sigdata(:,1) ; 8 ysig = sigdata(:,2); 9 10 v = linspace (0 , 2e7 , le6) ; °/0 vektor nezávislých proměnných - n velikosti rychlosti sigma = interpl (xsig , ysig, v, 'pchip', 0); °/0 zjištěni 12 funkcnich hodnot sigma(v) 13 °/0 Maxwell - Boltzmannovo rozděleni 14 Fv = 4*pi*v."2 . *(m/(2*pi*kB*T))~1•5.*exp(-m*v."2/(2*kB*T))í % 15 rozdělovači funkce kr_MB = trapz(v, v . *Fv . * sigma) °/0 rychlostni koeficient 16 vmean = trapz(v, Fv.*v); 17 18 °/0 Rozdělovači funkce pro skoro monoenergeticky svazek 19 v0 = vmean ; 20 width = vmean/100; 21 Fv2 = l/(sqrt(2*pi)* width) * exp ( - (v-vO) . ~ 2/(2* width "2)) ; °/0 22 rozdělovači funkce kr_beam = trapz(v, v .* Fv2 .* sigma) °/0 rychlostni koeficient 23 Cvičení 5.5 Spusťte program a porovnejte výsledné rychlostní koeficienty. Všimnete si, že nejsou shodné přestože použité rozdělovači funkce mají stejnou střední rychlost. Vysvětlete, proč tomu tak je. ■ Cvičení 5.6 Upravte program tak, že v průběhu vykreslí <7r(v) a odpovězte na následující otázky • Jaký je ionizační práh v eV? • Na jaké hodnotě dosahuje účinný průřez maxima?. Cvičení 5.7 Spusťte program pro několik hodnot elektronové teploty (definována v e V na řádku 3). • Jak se změní rychlostní koeficienty, když zvyšujete teplotu? • Existuje taková teplota, pro kderou rychlostní koeficient monoenergetického svazku překročí rychlostní koeficient Maxwell-Boltzmannova rozdělení? Poskytněte vysvětlení, proč to je/není možné. Pokročilé cvičení 5.2 Přepište Program 7 tak, aby místo velikosti rychlosti pracoval s energií elektronu. Používání energie je ve fyzice plazmatu mnohem častější. "pN Existuje několik veřejně přístupných databází, ze kterých můžete získat účinné průřezy. Ne- 44 Chapter 5. Interakce částic v plazmatu jvíce vyčerpávající z nich je pravděpodobně LxCat, dostupná na lxcat.net. Další poměrně rozsáhlou databází je databáze ALADDIN, dostupná na www-amdis.iaea.org/ALADDIN/. Formulace problému Srážkový člen Reakční schéma Implementace Parametrická studie Závěr Books Articles Online resources 6. Rovnováha částic v plazmatu Jak víte, typické plazma obsahuje mnoho typů částic. Kromě neutrálních atomů a molekul v základním stavu se mohou vyskytovat excitované částice (popsáno v Úvod do fyziky plazmatu interaktivně: část I), kladné nebo záporné ionty a samozřejmě i elektrony, díky kterým je plazma udrženo. Kromě toho mohou mít různé typu částic různé teploty. Pro různé aplikace plazmatu je třeba znát teploty a koncentrace různých typů částic v plazmatu. Například pro bioaplikace je důležité vytvořit plazma, které bude mít vysokou elektronovou teplotu zatímco teplota těžkých částic (iontů, neutrálů, excitovaných neutrálů) bude nízká, aby nedošlo k tepelnému poškození tkáně. 6.1 Formulace problému V této poslední kapitole tohoto textu implementujeme program, který počítá koncentrace různých typů částic v argonovém výboji za atmosferického tlaku při zadané elektronové teplotě Te a teplotě těžkých částic Tg. Model bude opět implementován v Matlabu a bude nula-rozměrný. Toto prakticky znamená, že koncentrace každého typu částic bude pouze funkcí času, nikoliv polohy 6.1.1 na =na(t), což fyzikálně odpovídá homogennímu plazmatu. (6.1) Srážkový člen V aproximaci homogenního plazmatu, kterou jsme učinili se rovnice kontinuity pro typ částic a redukuje na následující obyčejnou diferenciální rovnici. dt (6.2) Rovnice vypadá formálně poměrně jednoduše, musíme ale vzít v potaz, že srážkový člen Sa je složitou funkcí teplot (které jsou naštěstí v našem modelu konstantní), a koncentrací různých typů částic. Obecně vyjádříme srážkový člen pro typ částic a jako sumu příspěvků od všech reakcí, kde částice typu a figurují jako produkty či reaktanty. Sa=^Sa,r- (6-3) 46 Chapter 6. Rovnováha částic v plazmatu Příspěvek Sa,r udává, kolik částic typu a vznikne nebo zanikne díky reakci r za jednotku času. Konkrétní tvar Sa,r samozřejmě závisí na řádu reakce a v našem modelu použijeme reakce druhého řádu (dvojčasticové reakce) a reakce třetího řádu (trojčástickové reakce). Jak už jejich jména naznačují, v prvním případě se pro uskutečnění reakce musí srazit dvě částice, v druhém případě se musí srazit částice tři. Pamatujte, že řád reakce vám řekne jenom to, kolik reaktantů do reakce vstupuje, neříká nic o tom, kolik má reakce produktů. Ve skutečnosti trojčásticové reakce neprobíhají, ale ve fyzice plazmatu je často používáme jako aproximaci mnoha dějů, zejména za vyšších tlaků. Trojčásticová reakce ve skutečnosti probíhá jako dvě bezprostředně navazující dvojčásticové reakce A + B—>C + D, (6.4) C + E—>F. (6.5) Pokud je druhá z těchto reakcí výrazně rychlejší a pravděpodobnější než ta první, můžeme tento proces popsat zjednodušeně, troj časticovou reakcí A + B + E—>D + F. (6.6) Příkladem dvojčásticové reakce je ionizace argonového atomu dopadem elektronu (později budeme tuto reakci označovat R4). e + Ar* —> 2e + Ar+ (R4) (6.7) V této reakci se srazí elektron s excitovaným atomem argonu a dodá mu energii potřebnou pro jeho ionizaci. Tato reakce se potom promítne do srážkových členů všech jejích reaktantů a produktů, Se,r4 = +h(Te,TĚ) -ne-nAr*, (6.8) ^Ar+,r4 = +h (Te, Tg) ■ ne • «Ar*, (6.9) >W,r4 = -h(Te,Tg)-ne-nAr*. (6.10) Vidíte, že příspěvek k srážkovému členu obsahuje kromě koncentrací částic i rychlostní koeficient reakce, který může záviset na teplotě elektronů i teplotě těžkých částic. Jiné teploty v tuto chvíli uvažovat nebudeme, protože za našich podmínek je relativně rozumné pracovat s dvojteplotním modelem, tedy modelem, ve kterém mají elektrony svou vlastní teplotu, ale všechny ostatní typu částic (neutrály, ionty, excitované neutrály) jsou v termodynamické rovnováze a mají teplotu Tg. Můžete si také všimnout, že jednotlivé příspěvky Sa,r výše se liší znaménkem, jehož volba znaménka je velmi intuitivní. Vzhledem k tomu, že reakce (R4) vytvoří jeden nový elektron a jeden nový iont, jsou znaménka u odpovídajících příspěvků ke srážkovým členům kladná. Na druhou stranu reakce spotřebuje jeden excitovaný argonový atom a odpovídající příspěvek ke srážkovému členu 5at*,r4 je tedy záporný. Příkladem trojčásticové srážky je například vznik molekulárního argonového iontu. Existence této částice vás možná překvapí vzhledem k tomu, že argon je inertní plyn, ale ve vysokotlakém plazmatu hraje tato částice poměrně důležitou roli. Reakce, které vede na vznik tohoto iontu je Ar+ + 2Ar —> Ar+ + Ar. (R8) (6.11) Analogicky k dvoj časticovým reakcím přispívá tato trojčásticová reakce ke srážkovým členům následovně ArJ,r8 +h(Tg) ■ nAľ ■ nAľ ■ nAr+, (6.12) ^Ar+,r8 = +h{Tg) ■ «Ar • «Ar • «Ar+, (6.13) ^Ar,r8 = ~h(Tg) ■ nAr ■ nAr ■ nAľ+. (6.14) 6.1 Formulace problému 47 Tentokrát může rychlostní koeficient záviset jenom na teplotě těžkých částic TĚ, protože do reakce nevstupují elektrony. Navíc se musí dát pozor na to, že rychlostní koeficient bude mít jinou jednotku (viz cvičení níže). Cvičení 6.1 V naší formulaci má srážkový člen vždy jednotku 1 / (m3 • s). Jaká je odpovídající fyzikální interpretace? Jaká je jednotka rychlostního koeficientu kr pro dvojčásticové a trojčásticové reakce? ■ 6.1.2 Reakční schéma Reakční schéma pro náš OD model bylo převzato z článku od Baevy a kol., který se zabývá simulacemi mikrovlnného pochodňového výboje v argonu [Bae+12]. Typy částic, které model obsahuje jsou uvedeny v tabulce 6.1. Table 6.1: Seznam typů částic zahrnutých do modelu argonového plazmatu Ar atom argonu v zákl. stavu (koncentraci získáme ze stavové rovnice) Ar* excitovaný atom argonu (sdružuje resonanční a metastabilní 4s stavy) Ar+ atomový argonový iont ArJ molekulární argonový iont e elektron Mělo by být zdůrazněno, že rovnici kontinuity (6.2) budeme řešit jenom pro 4 typy částic, Ar*, Ar+, ArJ and e. Koncentrace argonového atomu v základním stavu totiž může být spočtena ze stavové rovnice ideálního plynu, která pro jednoatomové plyny platí velmi dobře. Pokud by argon nebyl vůbec ionizován, byla by koncentrace argonových atomů jednoduše P «Ar = i-^r (6-15) &b Tg kde p je tlak, &b je Boltzmannova konstanta a Tg je teplota těžkých částic. Nicméně, rovnice (6.15) neplatí, pokud je část argonu ionizována nebo excitována. Cvičení 6.2 Jaká je skutečná koncentrace argonových atomů v základním stavu pokud koncentrace excitovaných atomů je «Ar*> koncentrace atomových iontů je nAľ+ a koncentrace molekulárních iontů je nAľ+ ? ■ Typy částic uvedené v tabulce 6.1 mohou vstupovat do celkem 11 reakcí. Úplnou sadu těchto reakcí potom naleznete v tabulce 6.2. Toto reakční schéma předpokládá, že rozdělovači funkce elektronů i těžkých částic jsou Maxwellovské. Díky tomuto předpokladu a pokročilému fitování můžeme vyjádřit rychlostní koeficienty reakcí jako analytické funkce Te a Tg. Cvičení 6.3 Reakce (Rl), (R3) a (R4) v tabulce 6.2 obsahují exponenciální část exp(—C/Te), která je nulová při nízké elektronové teplotě a blíží se jedné při vyšších elektronových teplotách. • Jaká je jednotka a fyzikální význam konstant 11.65, 15.76a4.11? • Podle rychlostních koeficientů zkuste odhadnout, který z ionizačních kanálů bude pravděpodobně nej důležitější. 48 Chapter 6. Rovnováha částic v plazmatu Table 6.2: Reakční schéma pro vysokotlaké argonové plazma (převzato z [Bae+12]). Elektronová teplota je vždy vyjádřena v eV a teplota plynu v K. reakce rychlostní koeficient jednotka (Rl) e + Ar—>e + Ar* 4.9- 1(T15 • ^/Te-exp(-U.65/Te) m3/s (R2) e + Ar*—>e + Ar 4.8 -KT16-^ m3/s (R3) e + Ar—>2e + Ar+ 1.27 • 1(T14 • exp (-15.76/7;) m3/s (R4) e + Ar*—>2e + Ar+ 1.37 • 1(T13 • exp (-4.11/Te) m3/s (R5) e + e + Ar+—>e + Ar 8.75 • 10~39 • T~4-5 m6/s (R6) e + Ar+^Ar + Ar* 1.04 • 10'12 • (re/0.026)-°-67 • t^Texp^ís/r.) mVs (R7) e + Ar+—>e + Ar + Ar+ 1.11 • 10 12 • exp (-2.94 - 3 ^„X) mVs (R8) Ar++2Ar—>Ar + Ar+ 2.25 • 10~43(Tg/300)-04 m6/s (R9) Ar + Ar+—>Ar++2Ar 0.522 • ÍO"15^1 exp(-15131/7g) m6/s (RIO) Ar*+Ar*—>e + Ar+ + Ar 6.2-10~16 m3/s (Rll) Ar*+Ar—>Ar + Ar 3.0-10~21 m3/s Cvičení 6.4 Vyjádřete srážkové členy pro následující typy částic: e, Ar*, Ar+, ArJ pomocí rychlostních koeficientů k\ alku a koncentrací jednotlivých typů částic, nAľ,ne,riAľ*,nAľ+,nAľ+. m Jak vidíte, i pro tak jednoduchý plyn jakým je argon je počet reakcí poměrně vysoký. V plazmatu vzroste počet přípustných reakcí velmi dramaticky, pokud pracujeme v molekulárním plynu nebo ve směsi plynů. U molekulárních plynů musíme vzít v potaz rotační a vibrační excitace, disociaci a reasociaci molekuly. Jen pro představu o složitosti, pro minimální popis plazmatu ve směsi O2/N2 byste potřebovali nejméně několik desítek reakcí. Pro popis plazmatu v argonu a vlhkém vzduchu roste počet reakcí na cca 4000 [GB 13]. Implementace Nyní již známe všechno pro to, abychom mohli vyřešit systém obyčejných diferenciálních rovnic v následujícím tvaru d_ di ( "e \ V "Ar+ / SAľ+ \ sa4 ) (6.16) Formálně podobnou soustavu ODR jsme řešili v kapitole 3. Pokud si tedy nejste jistí ohledně čehokoliv, co se týká kódu a zejména řešiče ode45(), zkuste se znovu podívat na cvičení ve zmíněné kapitole. Vzhledem k tomu, že už se nejedná o Váš první program v Matlabu, poskytne vám tento studijní text pouze šablonu pro celý program a chybějící části budete muset doplnit. Spustitelný Matlabový skript, který obsahuje počáteční podmínky a volání řešiče vypadá následovně: Program 8: solve.m % nastaveni globalnich proměnných globál p; globál kB; globál Tg; 6.2 Implementace 49 global Te ; p=le5; '/, tlak v Pa kB = 1.38e-23; '/, Boltzmannova konst . v m~2 kg s~-2 K~-l Tg = 400; °/. Teplota plynu v Kelvin Te = 2; °/0 Teplota elektronu v eV tsteps = 1000; °/0 počet casov ch kroku tspan = logspace (-11 , -6, tsteps); °/0 pro casove kroky nepoužijeme lineárne rostouci hodnoty, ale logaritmicky rostouci °/0 Pocatecni koncentrace v m~-3 nArs_init = lel2; °/, Ar* nArp_init = lel2; °/, Ar + nAr2p_init = lel2; °/, Ar_2 + ne_init = nArp_init + nAr2p_init ; °/0 poc . podminka pro elektrony dana kvazineutralitou °/0 Vektor pocatecnich koncentraci initial = [nArs_init , nArp_init , nAr2p_init, ne_init]; % reseni soustavy 0DR [t, sol] = ode45( ' odefun ' , tspan, initial); % Vykresleni dat close all figure ; hold on; nAr = p./(kB*Tg)-sol(:,2) plot (loglO(t) , loglO(nAr) plot(loglO(t) , logl0(sol( plot(loglO(t) , logl0(sol( plot (loglO(t) , logl0(sol( plot(loglO(t) , logl0(sol( ylim([12, 26]) legend('Ar', 'Ar~*', 'Ar~ + ', 'northwest') 5*sol ( C) ; ,D) ,2)) ,3)) ,4)) ,3)-sol(:, ) ) ; ); ); ); Ar 2~ + electrons 'Location xlabel('Time, log_{10} [s]'); ylabel('Number density, log_{10} [m~{-3}]'); title(['p=', num2str(p), ' Pa, T_g=' , num2str(Tg), ' K, T_e=' , num2str(Te), ' eV ']) ; set(gca,'fontsize ' , 16); Ve zdrojovém kódu výše by vás nic nemělo překvapit. Největší rozdíl oproti našim předchozím kódům řešícím soustavu ODR je na řádku 12. Konkrétně zde nepoužíváme lineárně rostoucí vektor pro nezávislou proměnnou (čas), ale používáme vektor, který roste logaritmicky. Je to proto, že většina veličin v plazmatu se v čase mění velmi rychle. Stejně jako v kapitole 3, volá program solve .m řešič ode45 (), abychom získali časový vývoj soustavy ODR. Funkce, která vrací pravé strany soustavy diferenciálních rovnic (6.16) se opět nazývá odefunQ a je uložena v matlabovém souboru s odpovídajícím jménem. 50 Chapter 6. Rovnováha částic v plazmatu Program 8: odefun.m function dqdt = odefun(t, q), % Použijeme globálni proměnné definovane v solve.m globál p; globál kB; globál Tg; globál Te; °/0 vektor q obsahuje koncentrace Ar*, Ar+, Ar_2+ a elektronu nArs = q(1); nArp = q(2); nAr2p = q(3); ne = q(4) ; °/0 koncentrace argonu je spocitana pomoci stavové rovnice a koncentraci zbyvajicich typu castic nAr = p ./(kB*Tg)-nArs-0.5*nAr2p-nArp; % Rychlostni koeficienty jsou uloženy v oddelených souborech s funkcemi kl = f _kl(Te , Tg) k2 = f _k2(Te , Tg) k3 = f _k3(Te, Tg) k4 = f _k4(Te, Tg) k5 = f _k5(Te, Tg) k6 = f _k6(Te, Tg) k7 = f _k7(Te, Tg) k8 = f _k8(Te, Tg) k9 = f _k9(Te, Tg) klO = f_klO(Te, T kil = f_kll(Te , T % Nyni musime definovat srážkové cleny SArs = +kl*ne*nAr . . . -k2*ne *nArs . . . -k4*ne *nArs . . . + k6*ne*nAr2p . . . -2*klO*nArs*nArs . . . +kll*nArs*nAr; SArp = . . . ; SAr2p = . . . ; Se = ... ; % a na zaver je vrátime jakožto derivace vstupniho vektoru dqdt = [SArs; SArp; SAr2p; Se] ; end Ve funkci odef un() je třeba doplnit srážkové členy pro jednotlivé typy částic. Srážkový člen pro Ar* (proměnná SArs) už byl předvyplněn, ale srážkové členy pro Ar+ (proměnná SArp), ArJ (proměnná SAr2p) a elektrony (proměnná Se) musíte doplnit na základě cvičení 6.4. Navíc skript používá rychlostní koeficienty kl až kil definované pomocí externích funkcí f _kl () až f _kll (), které si berou na vstupu teplotu těžkých částic Tg v Kelvin a teplotu elektronů Te v elektronvoltech. Samozřejmě, ne všechny rychlostní koeficienty budou funkcemi 6.2 Implementace 51 jak Tg tak Te ale pro větší konzistentnost jsou funkcím f _kl () až f _kl 1 () předány obě teploty. Musíte tedy vytvořit 11 jednoduchých funkcí f _kl () až f _kl 1 () na základě reakcí v tabulce 6.2. Cvičení 6.5 Dokončete program výše, konkrétně 1. Do souboru odef un. m doplňte výrazy pro srážkové členy. 2. Vytvořte 11 souborů s funkcemi, f _kl .m až f _kll .m. Každý z těchto souborů bude počítat rychlostní koeficient odpovídající reakce. Vezměte v potaz, že funkce musí být schopny práce s vektorovými argumenty a používejte s Tg a Te pouze prvkové operace (. *, . /, atd...). 3. Spusťte kód s předdefinovanými hodnotami tlaku p, Tg and Te. Ověřte, zda výstup vypadá podobně jako následující obrázek. p=100000 Pa, T =400 K, T =2eV ľ g e 26 r Time, log10[s] Prosím vezměte v potaz, že jste první studenti pracující s tímto textem a není vyloučeno, že autoři neudělali chybu při formulaci srážkových členů. Pokud váš výstup vypadá lehce odlišně, a jste přesvětčeni o správnosti vašich srážkových členů, neváhejte protokol odevzdat. Nyní, když program funguje, můžete zkusit změnit různé parametry, abyste získali kvalitativní představu o tom, jak argonové plazma funguje. Cvičení 6.6 Zkuste zvýšit či snížit elektronovou teplotu a odpovědět na následující otázky. • Jak a proč se mění rovnovážná koncentrace elektronů? • Jak se mění doba potřebná pro zápal plazmatu? • Jaký je dominantní iont ve fázi zápalu a který iont je dominantní po stabilizaci plazmatu? Mění se relativní zastoupení iontů s elektronovou teplotou? Cvičení 6.7 Spusťte program pro vámi zvolenou hodnotu teploty elektronů a pro tlaky 103, 104, 105 a 106 Pa. Odpovězte na tyto otázky • Jak a proč se mění rovnovážná koncentrace elektronů? • Jak se mění doba potřebná pro zápal plazmatu? 52 Chapter 6. Rovnováha částic v plazmatu Parametrická studie Složitější programy podobné tomu, který jsme vyvinuli v tomto oddíle se často ve fyzice plazmatu používají pro zkoumání vlastností plazmatu za různých podmínek, ačkoliv počet reakcí je většinou výrazně vyšší (několik tisíc). Mnohdy je užitečné vědět, jak se bude složení plazmatu měnit pro různé teploty elektronů, teploty těžkých částic a tlaky. Pokročilé cvičení 6.1 Modifikujte program v předchozím oddíle tak, že bude řešit soustavu diferenciálních rovnic pro více hodnot Te or p a po ustálení vykreslí hustotu elektronů. Použijte při tom smyčku typu for. To náročně na tomto cvičení je fakt, že doba zápalu se s teplotou elektronů a tlakem mění. Proto musíte v každém kroku smyčky for upravit délku časového intervalu podle aktuálních hodnot Te a p. V opačném případě bude řešení trvat velice dlouho. Pokud se řeší podobné systémy diferenciálních rovnic pro mnoho typů částic a pro mnoho ^ podmínek, musí se použít speciální numerické postupy, aby řešení trvali přijatelnou dobu. Jedním z jednoduchých triků je transformace rovnice kontinuity do logaritmického tvaru. Nezávislou proměnnou potom není koncentrace částic, ale její logaritmus lj = lnnj. (6.17) Tento trik může konvergenci velmi urychlit, protože zatímco hodnoty n j se mění napříč mnoha řády, hodnoty lj se většinou mění jen v rámci jednoho řádu. Pokročilé cvičení výše může být zbytečně časově náročné zejména pro studenty, kteří neplánují Matlab dále používat. Nicméně, i v případě, že jste toto cvičení nedokončili, dokončete alespoň to následující. Cvičení 6.8 Jak již bylo zmíněno, složení plazmatu se často dramaticky mění s elektronovou teplotou. Obrázek níže ukazuje rovnovážnou koncentraci různých typů částic jako funkci elektronové teploty při konstantním tlaku p = 105 Pa a teplotě Tg = 400 K. Promyslete si následující otázky a odpovězte na ně • Pokud byste chtěli plazma použít jako zdroj světla, jakou elektronovou teplotu byste zvolili a proč? • Poskytněte jednoduché vysvětlení pro pokles koncentrace ArJ s rostoucí Te? p=100000 Pa, T =400 K g 60 0.5 1 1.5 2 2.5 3 Electron temperature, [eV] 6.4 Závěr 53 Závěr Tím se dostáváme na konec tohoto studijního textu. Gratulujeme, pokud jste došli až sem a podařilo se Vám dokončit všechna cvičení. Srdečně doufáme, že Vám cvičení a příklady v tomto textu alespoň trochu pomohly prohloubit pochopení komplexních procesů, které se v plazmatu odehrávají. Také doufáme, že schopnosti programovat v Matlabu, které jste se zde naučili, upotřebíte tom odvětví fyziky, kterým se rozhodnete zabývat. Fields gauged by Poritv MORE PORE ~ > Sociology is Ps^HOLosr is bioí-ocy is w«* '$ ^st TUST APPLIED JtóT APPLIED JUST APPLIED APFUED PKY9CS. PSYCHOLOGY BIOLOGY. CHEMISTRY IT'S MICE TO oh, hd; i didn't 5EEYD0 GľťSňU "IHĚ WW OVEJ? TH£T!E. SOdOLCClSTS PSYCHOLOGISTS BIOLOGISTS CHPIISTS PHYSIÍJSTS matheiwiciakjs Figure 6.1: Staženo z xkcd.org Knihy [Bit04] [GK08] J. A. Bittencourt. Fundamentals of Plasma Physics. 3rd edition. New York: Springer, 2004 (cited on page 12). Dan M. Goebel and Ira Katz. Fundamentals of Electric Propulsion: Ion and Hall Thrusters. 1st edition. Pasadena: California Institute of Technology, 2008 (cited on page 16). Články [Bae+12] [Bakl2] [GB13] [Shp+13] M. Baeva et al. "Modeling of microwave-induced plasma in argon at atmospheric pressure". In: Physical Review E 85.5 (May 2012), page 056404. ISSN: 1539-3755. DOl:10.1103/PhysRevE.85.056404. URL: http://link.aps.org/doi/10. 1103/PhysRevE. 85.056404 (cited on pages 47, 48). Daniel N Baker. "New Twists in Earth's radiation belts". In: American Scientist 102 (2012) (cited on page 21). W Van Gaens and A Bogaerts. "Kinetic modelling for an atmospheric pressure argon plasma jet in humid air". In: Journal of Physics D: Applied Physics 46.27 (2013), page 275201. URL: http: //stacks. iop. org/0022- 3727/46/i=27/a=275201 (cited on page 48). Yuri Y. Shprits et al. "Unusual stable trapping of the ultrarelativistic electrons in the Van Allen radiation belts". In: Nature Physics 9.11 (2013), pages 699-703. ISSN: 1745-2473. DOl: 10 .1038/nphys2760. URL: http://www.nature.com/ doif inder/10.1038/nphys2760 (cited on page 21). Online zdroje [FEI10] FEI. Introduction to Electron Microscopy. July 2010. URL: http: //www. f ei . com/ documents/introduction-to-microscopy-document/ (cited on page 15). 56 Chapter 6. Rovnováha částic v plazmatu [Marl3] Stefano Markidis. Lectures on Computational Plasma Physics. Sept. 2013. URL: https : / /www. pdc . kth . se/education/computational-plasma-physics/ computational-plasma-physics (cited on page 21). [Scil3] Scitechdaily. Scientists Explain the Formation of the Unusual Third Van Allen Radiation Ring. Sept. 2013. URL: http : / /scitechdaily . com/scientists-explain - format ion - unusual - third - van - alien - radiat ion - ring/ (cited on page 21).