Interakce laserového pulsu s plazmatem v souvislosti s inerciální fúzí zapálenou rázovou vlnou Laser-Plasma Interaction in the Shock Ignition Context

České vysoké učení technické v Praze Fakulta jaderná a fyzikálně inženýrská Katedra fyzikální elektroniky

Interakce laserového pulsu s plazmatem v souvislosti s inerciální fúzí zapálenou rázovou vlnou Laser-Plasma Interaction in the Shock Ignition Context BAKALÁŘSKÁ PRÁCE Autor práce: Petr Valenta Vedoucí práce: Ing. Ondřej Klimo, Ph.D. Konzultanti: Prof. Ing. Jiří Limpouch, CSc. Ing. Milan Holec Dr. Stefan Weber Akademický rok: 2013/2014

ČESKÉ VYSOKÉ UČENÍ TECHNICKÉ V PRAZE

FAKULTA JADERNÁ A FYZIKÁLNĚ INŽENÝRSKÁ Katedra fyzikální elektroniky

ZADÁNÍ BAKALÁŘSKÉ PRÁCE

Student:

Petr V a l e n t a

Obor:

Inženýrská informatika

Zaměření:

Informatická fyzika

Školní rok:

2013/2014

Název práce:

Interakce laserového pulsu s plazmatem v souvislosti se inerciální fúzí zapálenou rázovou vlnou Laser-Plasma Interaction in the Shock Ignition Context

Vedoucí práce:

Ing. Ondřej Klimo, PhD.

Konzultant:

Prof. Jiří Limpouch, CSc. Dr. Stefan Weber Ing. Milan Holec

Cíl práce: Cílem této práce je studium interakce laserového záření s plazmatem pro podmínky současných experimentů (zejména těch na laseru PALS) studujících možnosti zapálení inerciální termojaderné fúze silnou rázovou vlnou – tzv. Shock ignition. Bude se jednat zejména o porovnání interakce z hlediska absorpce a absorpčních procesů, rozptylu a spektra rozptýleného záření a vzniku rychlých elektronů pro různé frekvence (především základní a třetí harmonickou laseru PALS) případně i intenzity dopadajícího laserového pulsu. Studium bude probíhat pomocí částicových simulací a jako jejich vstup budou použity výsledky hydrodynamických výpočtů poskytnutých Ing. Holcem.

Pokyny pro vypracování: 1. Prostudujte teorii interakce laserového záření s plazmatem. Zaměřte se na srážkovou absorpci a na parametrické nestability (zejména stimulovaný Ramanův a Brillouinův rozptyl). 2. Prostudujte poskytnutý Particle-in-Cell (PIC) kód LPIC a implementujte okrajovou podmínku pro efektivní absorpci horkých elektronů a otestujte na modelových příkladech. 3. Dále použijte PIC kód k simulaci interakce laserového záření s plazmatem. Profily hustoty, teploty a rychlosti expanze plazmatu získáte z hydrodynamických simulací Ing. Holce. Tyto profily aproximujte a použijte v PIC simulacích. Studujte absorpci a rozptyl záření, spektrum rozptýleného záření a vznik rychlých elektronů v závislosti na parametrech laserového svazku. Výsledky porovnejte a zasaďte do kontextu současného výzkumu zapálení inerciální fúze pomocí silné rázové vlny. Literatura:

1) W. L. Kruer, The Physics of Laser Plasma Interactions, Addison-Wesley Publishing, Redwood City, 1988. 2) S. Atzeni, J. Meyer-ter-Vehn, The Physics of Inertial Fusion: Beam Plasma Interaction, Hydrodynamics, Hot dense Matter, Clarendon Press, Oxford, 2004. 3) S. Eliezer, The Interaction of High-Power Lasers with Plasmas, Series in Plasma Physics , Taylor & Francis, 2002. 4) R. Lichters, R. E. W. Pfund and J. Meyer-ter Vehn, LPIC++ a parallel onedimensional relativistic electromagnetic particle-in-cell-code for simulating laserplasma interactions, Max-Planck Institute fur Quantenoptik, Garching, Tech. Rep. 225, (1997). 5) M. Mašek, Eulerova Vlasovova metoda pro laserové plazma, disertační práce na MFF UK, Praha 2006. Datum zadání:

říjen 2013

Datum odevzdání:

7.červenec 2014

…………………………. Vedoucí katedry

………………………….. Děkan

Prohlášení Prohlašuji, že jsem tuto práci vypracoval samostatně a použil jsem pouze podklady uvedené v přiloženém seznamu. Nemám závažný důvod proti použití tohoto školního díla ve smyslu § 60 Zákona č. 121/2000 Sb., o právu autorském, o právech souvisejících s právem autorským a o změně některých zákonů (autorský zákon). ........................................ V Praze dne .................... Petr Valenta

Název práce: Autor: Druh práce: Obor: Zaměření:

Interakce laserového pulsu s plazmatem v souvislosti s inerciální fúzí zapálenou rázovou vlnou Petr Valenta Bakalářská práce Inženýrská informatika Informatická fyzika

Vedoucí práce:

Ing. Ondřej Klimo, Ph.D. Katedra fyzikální elektroniky, Fakulta jaderná a fyzikálně inženýrská, České vysoké učení technické v Praze

Konzultanti:

Prof. Ing. Jiří Limpouch, CSc. Katedra fyzikální elektroniky, Fakulta jaderná a fyzikálně inženýrská, České vysoké učení technické v Praze Ing. Milan Holec Katedra fyzikální elektroniky, Fakulta jaderná a fyzikálně inženýrská, České vysoké učení technické v Praze Dr. Stefan Weber Fyzikální ústav Akademie věd České republiky, v. v. i.

Abstrakt:

Klíčová slova:

Práce představuje výsledky PIC simulací interakce laserového impulsu s plazmatem v souvislosti s inerciální fúzí zapálenou rázovou vlnou. Byly provedeny simulace svazku laseru PALS na základní vlnové délce 1,315 µm s intenzitou 1 · 1016 W/cm2 a na vlnové délce 0,438 µm odpovídající třetí harmonické základní frekvence s intenzitou 5 · 1015 W/cm2 . Oba impulsy trvaly 100 ps. K simulacím byl využit výpočetní kód LPIC++. Počáteční profily hustoty a teploty byly zvoleny z hydrodynamických výpočtů. Celková reflektivita pro první harmonickou činila 43 %, většina energie byla absorbována hustotními kavitami a parametrickými nestabilitami, v důsledku čehož docházelo k produkci vysokého množství rychlých elektronů. Celková reflektivita pro třetí harmonickou byla 59 %, dominantním mechanismem absorpce byly coulombické srážky mezi elektrony a ionty. Energie rychlých elektronů byla v případě simulace první harmonické 170 keV, pro třetí harmonickou byla 25 keV. V případě třetí harmonické by míra předehřátí palivového terče těmito elektrony pravděpodobně nebyla příliš významná, v případě první harmonické by značná část rychlých elektronů prošla komprimovanou slupkou palivového terče, což je pro zapálení inerciální fúze rázovou vlnou nevhodné. Laser, plazma, inerciální fúze, Particle-in-Cell

Title: Author: Type of work: Branch of study: Specialization:

Laser-Plasma Interaction in the Shock Ignition Context Petr Valenta Bachelor thesis Engineering informatics Computational physics

Supervisor:

Ing. Ondřej Klimo, Ph.D. Department of Physical Electronics, Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University in Prague

Consultants:

Prof. Ing. Jiří Limpouch, CSc. Department of Physical Electronics, Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University in Prague Ing. Milan Holec Department of Physical Electronics, Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University in Prague Dr. Stefan Weber Institute of Physics, Academy of Sciences of the Czech Republic, v. v. i.

Abstract:

Keywords:

The thesis presents the results of PIC simulations of laser-plasma interaction in the context of shock ignition. We have made simulations of laser system PALS at the fundamental wavelength 1.315 µm with intensity 1 · 1016 W/cm2 and at the wavelength 0.438µm corresponding to the third harmonic of the fundamental frequency with intensity 5 · 1015 W/cm2 . Pulse duration was 100 ps for both simulations. We have used a parallel code LPIC++. The initial density and temperature profiles were obtained from the hydrodynamic simulations. The total reflectivity for the first harmonic was 43 %, where most of the energy has been absorbed by density cavities and parametric instabilities, which leads to the production of high quantities of hot electrons. The total reflectivity of the third harmonic was almost 59 %, the dominant mechanism of absorption was Coulombic collisions. Energy of hot electrons for the simulation of first harmonic was 170 keV, for the third harmonic was 25 keV. The energy of these electrons in case of third harmonic is not sufficient to preheat the fuel significantly. In case of first harmonic, high quantities of hot electrons would probably passed into the hot-spot, which may be inappropriate for shock ignition. Laser, plasma, inertial fusion, Particle-in-Cell

Obsah Úvod

11

1 Inerciální fúze 1.1 Princip jaderné syntézy . . 1.2 Lawsonovo kritérium . . . 1.3 Přímé a nepřímé zapálení 1.4 Nové metody . . . . . . . 1.5 Dosavadní výsledky . . . .

. . . . .

13 13 14 15 16 17

. . . . . . . . .

19 19 21 24 26 28 29 31 33 34

. . . . . .

37 38 41 43 45 47 48

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

2 Laser-plazma interakce 2.1 Základní parametry plazmatu . . . . . . . 2.2 Způsoby popisu plazmatu . . . . . . . . . 2.3 Šíření elektromagnetických vln v plazmatu 2.4 Srážková absorpce . . . . . . . . . . . . . . 2.5 Landauův útlum . . . . . . . . . . . . . . 2.6 Ponderomotorická síla . . . . . . . . . . . 2.7 Parametrické nestability . . . . . . . . . . 2.7.1 Stimulovaný Ramanův rozptyl . . . 2.7.2 Stimulovaný Brillouinův rozptyl . . 3 Numerické simulace 3.1 Metoda Particle-In-Cell . . . . 3.1.1 Integrátor polí . . . . . 3.1.2 Váhování částic a polí 3.1.3 Integrátor částic . . . . 3.1.4 Stabilita metody . . . 3.2 Výpočetní kód LPIC++ . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

. . . . .

. . . . . . . . .

. . . . . .

4 Výsledky 51 4.1 Okrajová podmínka . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.2 Simulace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

Závěr

61

Literatura

63

Přílohy

69

A Vstupní soubory 69 A.1 Simulace první harmonické . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 A.2 Simulace třetí harmonické . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 B Obsah CD

71

Úvod Potřeba nových zdrojů energie pro budoucnost nikoliv vzdálenou odstartovala v padesátých letech minulého století intenzivní bádání směřující k ovládnutí jaderného slučování. Snaha experimentátorů, kteří se pokoušeli realizovat potřebné podmínky ke sloučení jaderného paliva, byla v té době převážně orientována na plazma držené ve vhodně konfigurovaném magnetickém poli [16]. Rychlý úspěch, který byl vědeckou komunitou očekáván, se však nedostavil. Tento fakt také přispěl k tomu, že se v šedesátých letech zrodila myšlenka využít k zapálení termojaderné fúze tehdy zcela nový zdroj intenzivního koherentního záření, laser. Okamžitě ale bylo také vzneseno varování, že může docházet k celé řadě jevů, které mohou dosažení termojaderné fúze tímto způsobem vážně zkomplikovat. Posléze se tak skutečně stalo. Jedním z nejvíce nepříznivých jevů z hlediska inerciálně držené fúze je generace rychlých elektronů, které vznikají v plazmatické koroně, v níž dochází k přeměně energie laserového záření na kinetickou energii plazmatu. Tyto elektrony dosáhnou centrální oblasti komprimovaného terče ještě před příchodem rázové vlny a značně jej předehřejí. Stlačení terče na požadované hodnoty termojaderného slučovaní je potom mnohem obtížnější. Proto je interakce laserového impulsu s plazmatem v této souvislosti intenzivně zkoumána, současně se hledají nové metody jak deponovat co největší množství energie do terče při minimální produkci rychlých elektronů. Tato práce je zaměřena právě na interakci laserového záření s plazmatem pro podmínky současných experimentů Badatelského centra PALS v Praze, které studují možnosti zapálení inerciálně držené fúze pomocí jedné z nejmodernějších metod, využitím silné rázové vlny. Jedná se zejména o porovnání interakce z hlediska absorpce, absorpčních procesů, rozptylu a spektra rozptýleného záření a samozřejmě také z hlediska vzniku rychlých elektronů pro základní frekvenci a třetí harmonickou jódového fotodisociačního laserového systému PALS. Práce je členěna následovně. První kapitola přináší stručný náhled na problematiku jaderné fúze, včetně podmínek jejího dosažení pro oba hlavní přístupy k udržení plazmatu. Větší část je samozřejmě věnována inerciálnímu udržení. Tato kapitola mírně zabrousí i do historie výzkumu jaderné fúze, ale nevynechá ani nejaktuálnější výsledky. V druhé kapitole jsou shrnuty nejzákladnější poznatky z oblasti fyziky plazmatu a interakcí laserového impulsu s plazmatem. Jedná se o srážkovou absorpci, Landauův útlum, 11

Úvod ponderomotorickou sílu a zejména pak o dvojici parametrických nestabilit – stimulovaný Ramanův a Brillouinův rozptyl. Třetí kapitola je věnována numerickým simulacím, jakožto základnímu nástroji ke studiu nejen interakcí laserového impulsu s plazmatem. Dále bude důkladně popsána jedna z nejpopulárnějších metod ve fyzice plazmatu – Particle-in-Cell. Na závěr tato kapitola přináší stručný popis výpočetního kódu LPIC++, který byl použit k simulacím v rámci této práce. V poslední kapitole budou představeny výsledky a přínosy této práce. Jedná se o implementaci okrajové podmínky ve výpočetním kódu LPIC++ a otestování její funkčnosti na modelových příkladech. Dále pak o již zmíněné simulace a jejich porovnání z hlediska využitelnosti ve výzkumu inerciálně držené fúze. Ač většina odborných publikací z oboru fyziky plazmatu užívá cgs elektrostatické jednotky, pro celou tuto práci jsme zvolili soustavu jednotek SI (pokud nebude uvedeno jinak).

12

Kapitola 1 Inerciální fúze I přes nejrůznější úsporná opatření se energetická spotřeba lidstva neustále zvyšuje. Zásoby fosilních paliv, ropy a zemního plynu se pomalu ztenčují a obnovitelné zdroje poptávku po energii pokrýt nedokážou. Tento energetický deficit by se v budoucnu mohl stát velkou překážkou v trvale udržitelném rozvoji lidské společnosti [30]. Proto je nezbytné najít alternativní zdroj energie. Nejlépe takový, který dokáže definitivně vyřešit již zmíněné problémy a jeho využití bude navíc ekologicky šetrné. V současné době se zdá, že jediným kandidátem je energie ukrytá v atomovém jádře. Tu lze získávat dvěma způsoby. První spočívá v štěpení těžších jader na lehčí. Štěpná reakce však s sebou nese jistá bezpečnostní rizika. Její produkty jsou navíc silně radioaktivní a kvůli tomu vznikají vážné problémy s jejich dlouhodobým skladováním. Druhou, široké veřejnosti méně známou cestou je proces zcela opačný, tedy slučování neboli fúze. Tato kapitola přináší stručný náhled na problematiku jaderné fúze, speciálně pak přibližuje jeden z hlavních přístupů k jejímu dosažení – inerciálně drženou termojadernou fúzi.

1.1

Princip jaderné syntézy

Jaderná fúze je reakce, při které se sloučí jádra lehčích prvků, a vytvoří tak jádro nové, těžší. Pokud je hmotnost produktu menší než součet hmotností reaktantů, nutně se množství energie odpovídající tomuto hmotnostnímu schodku uvolní. Aby se však atomová jádra mohla sloučit, musí se k sobě vzájemně přiblížit na vzdálenosti srovnatelné s jejich rozměry, teprve pak se uplatní silná interakce zodpovědná za soudržnost nukleonů. Nicméně tomu brání odpudivé coulombické síly mezi souhlasně nabitými protony. Otázkou tedy je, jak tuto potenciálovou bariéru překonat. Můžeme například s výhodou využít kinetickou energii chaotického tepelného pohybu částic. Tento přístup, jenž vede ke sloučení atomových jader, se nazývá termojaderná fúze. 13

Kapitola 1. Inerciální fúze Reakce, která by měla být využitelná k produkci energie na Zemi, musí mít vysoký účinný průřez při relativně nízké požadované energii. Měla by být také přátelská k životnímu prostředí. Zároveň je nutné, aby prokázala ekonomickou konkurenceschopnost, tedy aby byla schopna snadno pokrýt energetické náklady vynaložené k jejímu zapálení. Zapálením je označován moment, kdy vlastní ohřev pomocí fúzních produktů vyrovná energetické ztráty. V tomto směru se jeví jako nejvýhodnější reakce izotopů vodíku deuteria a tritia 2 H + 3 H → n (14,06 MeV) + 4 He (3,52 MeV). (1) Deuterium lze poměrně snadno získat extrakcí z mořské vody, je tudíž okamžitě k dispozici. Naopak tritium je nestabilní a na Zemi se vyskytuje jen v malém množství, proto musí být v rámci palivového cyklu nějakým způsobem vyráběno. Jednou z možností je využít neutrony z fúzních reakcí k ozařování lithia. Přírodní lithium je hojně zastoupeno v zemské kůře [28]. Oba druhy paliva jsou tedy snadno dostupné, rovnoměrně geograficky rozložené a v podstatě nevyčerpatelné. Produktem jejich reakce je neutron a jádro helia, které není radioaktivní. Tyto atributy předurčují jaderné fúzi stát se globálním zdrojem energie.

1.2

Lawsonovo kritérium

První pokusy, jak uskutečnit fúzní reakci v pozemských podmínkách, směřují do poloviny třicátých let minulého století. Jednoduchá myšlenka spočívala v urychlení svazku iontů deuteria a jejich následné srážce s terčem z pevného tritia. Záhy se však zjistilo, že tento přístup kýžené ovoce nepřináší. Většina energie se totiž ztrácí v pružných coulombických srážkách, jejichž účinný průřez je o několik řádů vyšší než fúzní. Terč se sice ohřívá, dokonce dochází k fúzním reakcím, mnohem více energie se ale spotřebuje k urychlení částic [41]. Po odtajnění výzkumu fúze na obou stranách železné opony publikoval J. D. Lawson kritérium pro kladnou energetickou bilanci termojaderné reakce. Z něho jasně vyplývá, že součin hustoty a doby udržení dané reakce musí být větší než hodnota funkce závisející pouze na teplotě. Nutno podotknout, že tato funkce má pro každý druh reaktantů specifický průběh, nicméně minima nabývá vždy pro teploty značně vysoké. To je zřejmé, neboť je potřeba, aby částice získaly dostatečnou kinetickou energii k překonání již zmiňované coulombické bariéry, a vzrostla tak četnost fúzních reakcí. I přesto, že je potřebná teplota systému snížena existencí kvantově mechanického tunelového jevu, je při ní většina látek ve stavu plně nebo alespoň částečně ionizovaného plazmatu. Materiály s pevnou strukturou se však při kontaktu s ním začnou tavit. Stěžejní je tedy otázka, jak toto plazma udržet. Vezmeme-li v úvahu fakt, že plazma je směs volných elektricky nabitých částic, nabízí se využít vhodnou konfiguraci magnetického pole. Výsledné trajektorie pohybu částic jsou 14

1.3. Přímé a nepřímé zapálení pak určeny šroubovicemi podél magnetických siločár. Nicméně velikost intenzity magnetického pole je omezena strukturální pevností látky, lze tedy udržet pouze plazma s velmi nízkou hustotou. Z Lawsonova kritéria potom okamžitě plyne, že fúzní reakce musí probíhat relativně dlouho. Druhou možností, jak udržet horké plazma, je využít setrvačnosti hmoty neboli inerce. Díky ní lze zapálení termojaderné reakce provést dříve, než se objem plazmatu vlastním tlakem rozptýlí do prostoru. Setrvačnost udrží hmotu pohromadě po dobu závislou na rychlosti zvukové vlny, nejedná se tedy o žádný externí druh udržení. Doba udržení reakce je již z principu velmi krátká, a tudíž je potřeba dosáhnout extrémně vysokých hustot paliva [17].

1.3

Přímé a nepřímé zapálení

Uvolnění jaderné energie za pomoci inerciálního udržení bylo demonstrováno krátce po druhé světové válce explozí vodíkové bomby. V takovémto rozsahu jsou však devastační účinky uvolněné energie neoddiskutovatelné. Pro smysluplné průmyslové využití je žádoucí ji získávat kontrolovaně, v maximální míře odpovídající pouze několika miligramům směsi deuteria a tritia. V kryogenním stavu se pak takovéto množství paliva s rezervou vejde do sféry o průměru jednoho milimetru [30]. Do takto malé oblasti však nebylo možné deponovat dostatečné množství energie pro zapálení termojaderné fúze, a proto se dále o inerciálním udržení neuvažovalo. Vše se změnilo až na počátku šedesátých let, kdy T. H. Maiman zkonstruoval první fungující laser, emitující záření ve viditelné části spektra. Koncept zapálení termojaderné fúze pomocí laseru vychází ze zákona zachování hybnosti. Intenzivní laserové svazky jsou fokusovány na kulovou slupku obsahující směs jaderného paliva, tzv. peletu. Její vnější povrch se okamžitě odpaří v plazma a vysokou rychlostí začne expandovat do okolního vakua. Tento proces se nazývá ablace. Laserem generované plazma pak předá impuls neodpařené části slupky k silné implozi. Palivo je tedy komprimováno ablačním tlakem, nikoli tlakem laserového záření, jak by se na první pohled mohlo zdát. Pokud se peleta ozáří rovnoměrně ze všech stran, lze dosáhnout požadované komprese paliva a následně zapálení fúzní reakce. Při použití tohoto uspořádání mluvíme o přímo zapalované fúzi. Přímé zapálení s sebou však nese celou sérii problémů. V první řadě je to velmi vysoká náchylnost terče k hydrodynamickým nestabilitám, tedy deformacím. Ty vznikají jednak vlivem nesymetrického ozařování v důsledku konečného počtu oddělených, přirozeně koherentních svazků, dále pak také nehomogenitou a nedokonalou hladkostí povrchu terče. Dalším nezanedbatelným problémem je rychlé předehřívání termojaderného paliva. To je způsobeno generací populace tzv. rychlých elektronů v místech, kde laserové záření interaguje s ven proudícím řídkým plazmatem neboli koronou. Tyto horké elektrony pronikají do dosud nestlačeného paliva, a narušují tak adiabatičnost komprese [34]. 15

Kapitola 1. Inerciální fúze Později byl pro eliminaci nerovnoměrností v ozáření navržen prostředek, který využívá radiační transport v pouzdru obepínajícím fúzní palivo. Zhruba řečeno, laserové svazky se namíří na vnitřní stěnu dutiny vyrobené z materiálu s vysokým nukleonovým číslem, pro níž se vžil německý název hohlraum. Povrch této dutiny se okamžitě odpaří a přemění na husté plazma, které transformuje laserové záření na krátkovlnné měkké rentgenové záření. To komprimuje palivo a zároveň velmi dobře stírá nerovnoměrnosti v původním rozložení energie. Tento princip je velmi podobný konfiguraci vodíkové bomby a nazývá se nepřímé zapálení [28]. Evidentní nevýhodou nepřímo zapalované fúze je nižší účinnost konverze energie vnějšího zdroje na kinetickou energii implodující slupky. Jen malá část rentgenového záření je využita ke kompresi paliva, zbytek uniká vstupními otvory nebo se ztrácí radiační difúzí do stěn dutiny. To však nejsou jediné energetické ztráty. Plazma v okolí vstupních otvorů hohlraumu je řídké a téměř homogenní, což jsou vlastnosti, které vytváří ideální podmínky pro vznik parametrických nestabilit, zejména pak zpětného Ramanova a Brillouinova rozptylu [27], [34]. Potíže těchto nestabilit budou přiblíženy v závěru druhé kapitoly této práce.

1.4

Nové metody

Za více než padesát let výzkumu a vývoje v oblasti laserové techniky byl zaznamenán nesmírný pokrok. Charakteristiky moderních laserů se rapidně zlepšily, impulsy se i nadále zkracují a jejich intenzita narůstá do extrémních hodnot. Díky tomu otevírá tento zdroj koherentního záření cestu k řadě průlomových aplikací v nejrůznějších oblastech lidského bádání. Jak již bylo uvedeno, je to také jeden z mála dostupných nástrojů ke studiu inerciální fúze. Není tedy náhodou, že se v poslední době stále častěji objevují nové, důmyslnější metody jejího zapálení. Z předchozího vyplývá, že palivový terč je potřeba stlačit a zahřát na zápalnou teplotu. Nicméně čím více se terč zahřívá, tím roste jeho vnitřní tlak, který činí kompresi obtížnější [32]. Bylo by tedy vhodné oba procesy od sebe oddělit. Dále se ukázalo, že není potřeba zahřát na zápalnou teplotu celý objem palivového terče, což by bylo energeticky neúnosně náročné. Krátce před okamžikem stagnace, tedy v poslední fázi komprese, dochází totiž přirozeným způsobem uprostřed terče ke vzniku malé, silně zahřáté oblasti, označované jako hot spot. Pokud se slupka terče dostatečně urychlí, vznikne uvnitř při implozi vysoký tlak. Jeho vlivem se hot spot zahřeje natolik, že dojde k zapálení termojaderné reakce. Fúzní produkty pak nesou dostatečně velké množství energie k ohřevu okolního paliva pružnými coulombickými srážkami [1]. V souladu s těmito poznatky byla v nedávné době navržena metoda rychlého zapálení (Fast ignition), která efektivně využívá vliv rychlých elektronů. Palivový terč je komprimován stejně jako při konvenčním způsobu zapálení nanosekundovými laserovými impulsy, 16

1.5. Dosavadní výsledky přičemž nyní není nutností dosáhnout přesné symetrie a rovnoměrnosti. Nejde nám totiž o zapálení termojaderné reakce z centrálního hot spotu, stačí pouze dosáhnout jisté úrovně komprese. Hot spot se posléze vytvoří uměle na povrchu stlačeného terče pomocí elektronového svazku, který lze urychlit aplikací pikosekundového laserového impulsu [18], [34]. Tento přístup má však také své slabé stránky. Stále není objasněno, jak docílit toho, aby svazek rychlých elektronů prošel skrz plazmatickou koronu, aniž by ztratil významnou část své energie a aby její většinu deponoval v relativně malé oblasti, kde by mělo dojít k zapálení fúzní reakce [1]. Uvažuje se buď o zavedení třetího intenzivního laserového svazku, který by vytvořil kanál v nadkritickém plazmatu, nebo o umístění zlatého kuželu do blízkosti palivového terče, pomocí něhož lze plazmatickou koronu odstínit. Mezi nejmodernější metody ovšem patří koncept neizobarického zapálení inerciální fúze silnou rázovou vlnou (Shock ignition), který spočívá v jiném časovém tvarování laserových impulsů. Téměř adiabatickou kompresi palivového terče lze zajistit nanosekundovými impulsy o nižší hustotě výkonu. V okamžiku stagnace na ně navazuje časově velmi krátký, relativně intenzivní impuls, který vyvolá silnou sféricky konvergentní rázovou vlnu. Ta se během průchodu do středu terče při správném načasování sráží s vracejícími se primárními vlnami, určenými ke kompresi palivové směsi. Jejich vzájemná interakce prudce zesílí rázovou vlnu, a ta následně způsobí zapálení fúzní reakce v centru hot spotu [34], [2]. Je dokázáno, že tato metoda umožňuje zapálit palivo při třikrát nižší požadované energii vnějšího zdroje než při konvenčním izobarickém způsobu, což samozřejmě implikuje také vyšší zisky [37]. Dalším pozitivem je zmírnění Rayleigh-Taylorovy a dalších hydrodynamických nestabilit. Na druhou stranu, aplikací intenzivního laserového impulsu na již stlačený terč se nevyhneme nelineární interakci záření s plazmatickou koronou, která vyvolává převážně třísvazkové parametrické nestability [20], [21]. Ty způsobují odraz a rozptyl laserového záření a některé vedou ke vzniku rychlých elektronů. Tyto jevy lze z hlediska fúze považovat za škodlivé nebo alespoň nebezpečné, protože mohou omezovat efektivní absorpci laserového záření a stlačitelnost terče. Studium těchto nestabilit je proto klíčové ke zvládnutí zapálení inerciálně držené fúze.

1.5

Dosavadní výsledky

Na podzim roku 2013 se údajně podařilo vědcům z americké Lawrence Livermore National Laboratory dosáhnout významného milníku, poprvé získali pomocí inerciálně držené termojaderné fúze více energie, než do ní vložili na vstupu. Tyto zprávy byly potvrzeny a publikovány v únoru následujícího roku. V experimentu bylo využito 192 svazků laserového zařízení National Ignition Facility o celkové energii 1,8 MJ [14]. 17

Kapitola 1. Inerciální fúze Fúzní zisk je ovšem měřen v poměru k energii, které byla absorbována palivovým terčem, přičemž současná elektrická účinnost laserových systémů dosahuje pouze jednotek procent. Je to tedy energetická bilance z hlediska fyzikálního mechanismu na té nejnižší úrovni. Pro skutečnou produkci elektrické energie v hypotetické fúzní elektrárně bude potřeba, aby hodnota tohoto zisku byla ještě zhruba o dva řády vyšší [5]. Ač se tedy o žádnou senzaci nejedná, jsou i tyto výsledky významným průlomem v celém termojaderném výzkumu. Pro další kroky k cíli je ale zásadní kvalitnější informovanost politických představitelů široké laické veřejnosti. Výstavba velkých laserových zařízení je značně nákladná a nebýt jejich významného vojenského využití, pravděpodobně by nebyla zafinancována. Nicméně velké laserové systémy jsou široce využitelné také v řadě společensky významných mezioborových aplikací, mohou sloužit například pro vývoj protonových zdrojů určených k léčbě nádorových onemocnění [3]. I kdyby se tedy nakonec nedokázalo využít energii z inerciální fúze, jsou takto vynaložené prostředky účelné a smysluplné.

18

Kapitola 2 Laser-plazma interakce Dopadne-li intenzivní laserový svazek na terč, prakticky okamžitě dochází vlivem silného elektromagnetického pole ke vzniku vysoce ionizovaného plazmatu při jeho povrchu. Během tohoto procesu se vytváří plazmatická korona expandující směrem od terče rychlostí přibližně rovnou rychlosti zvuku. Při fúzních experimentech je nesmírně důležité deponovat co nejvíce energie laseru do palivového terče. Míra přenosu této energie ovšem velmi závisí na procesech probíhajících při průchodu laserového svazku plazmatickou koronou. A právě této problematice je věnována následující kapitola.

2.1

Základní parametry plazmatu

Plazma je kvazineutrální směs elektricky nabitých a neutrálních částic vykazující kolektivní chování [6]. Tato definice si vyžaduje vysvětlení některých pojmů. Kvazineutralitou rozumíme fakt, že ačkoliv v mikroskopických měřítkách tvoří nabité částice lokální elektrická pole, v makroskopickém měřítku se plazma jeví jako neutrální. To znamená, že celkové množství kladného a záporného náboje je v rovnováze. Matematicky lze tuto skutečnost vyjádřit jako X qs ns ≈ 0, (2) s

kde qs a ns je náboj resp. hustota částic druhu s. Suma probíhá přes každý druh částic vyskytující se v plazmatu. Pojmem kolektivní chování označujeme vzájemné působení nabitých částic zprostředkované makroskopickým elektromagnetickým polem. Coulombická interakce ovlivňuje částice na relativně velké vzdálenosti, plazma tak získává bohatý repertoár možných pohybů. Kolektivní chování je podstatné, nicméně nemusí vždy dominovat. 19

Kapitola 2. Laser-plazma interakce Jedním z nejdůležitějších parametrů, který nám umožňuje přesněji předpovídat chování plazmatu, je stupeň jeho ionizace. Ten lze snadno zjistit pro plyn v termodynamické rovnováze ze Sahovy rovnice, která je nejčastěji udávána v následující podobě, 3 εk+1 − εk nk+1 ne ≈ 2,4 · 1021 T 2 exp − . nk kB T





(3)

Zde nk je hustota k-násobně ionizovaných atomů, ne je hustota elektronů, T je teplota plynu, εk je ionizační energie k-té energetické hladiny atomu a kB Boltzmannova konstanta. Z rovnice jasně vyplývá, že existence plně ionizovaného plazmatu je podmíněna přítomností velmi vysoké teploty (vzhledem k ionizační energii). Charakteristickým rysem chování plazmatu je jeho schopnost odstínit elektrické potenciály. Proto se zavádí další důležitý parametr nazývaný Debyeova stínící délka, s

λD =

ε0 kB Te . e2 ne

(4)

Konstanta ε0 značí permitivitu vakua, e je velikost elementárního náboje a Te teplota elektronů. Často se stává, že mají různé druhy částic v plazmatu rozdílné teploty, i přesto však mohou být ve své vlastní tepelné rovnováze [6]. Debyeova délka vyjadřuje vzdálenost, na které potenciál pole poklesne na 1/e své původní hodnoty (zde symbol e označuje Eulerovo číslo). Je tedy jasné, že podmínka kvazineutrality může být splněna pouze tehdy, je-li charakteristický rozměr plazmatu L  λD . V opačném případě nemůžeme ionizovaný plyn nazvat plazmatem. V souvislosti s Debyeovou stínící délkou je definován velmi důležitý tzv. plazmatický parametr ND . Ten vyjadřuje počet elektronů vyskytujících se ve sféře o poloměru Debyeovy délky, platí tedy 4 (5) ND = πne λ3D . 3 Zdůrazněme, že mechanismus Debyeova stínění je platný pouze tehdy, je-li ND  1. V takovém případě hovoříme o ideálním plazmatu, pro jehož složky se používá stavová rovnice ideálního plynu. Při posunutí elektronů proti homogennímu iontovému pozadí, například pomocí intenzivního laserového impulsu, vzniká mimořádně silné elektrické pole, kterým hmotnější ionty táhnou elektrony zpět do původní polohy pro obnovení kvazineutrality plazmatu. Výsledkem je tlumené harmonicky oscilující elektrické pole nabitých částic kmitajících s tzv. elektronovou plazmovou frekvencí ωpe , s

ωpe = kde me je hmotnost elektronu. 20

e2 ne , ε0 me

(6)

2.2. Způsoby popisu plazmatu Podobně bychom mohli zavést i iontovou plazmovou frekvenci ωpi , ale vzhledem k tomu, že ionty jsou několikanásobně hmotnější než elektrony, nestíhají na vysokofrekvenční oscilující pole reagovat. Bereme je tedy většinou jako nehybné neutralizující pozadí. Pokud je ovšem frekvence vnějšího zdroje nebo módů indukovaných v plazmatu blízká této frekvenci, nemůžeme si tento zjednodušující předpoklad dovolit. Příkladem může být stimulovaný Brillouinův rozptyl, jehož odvození lze najít v samotném závěru této kapitoly. Vraťme se na okamžik k definici plazmatu. K oddálení nábojů dochází jen na určitou krátkou dobu. Kvazineutralitu má tedy smysl uvažovat jen v případě, že je charakteristický čas děje τ  1/ωpe . Má-li dále dominovat kolektivní působení nad binarním, musí být elektronová plazmová frekvence ωpe větší než frekvence coulombických srážek elektronů s ionty. Srážky mezi nabitými částicemi v plazmatu charakterizuje srážková frekvence νc , která je definována jako převrácená hodnota střední doby, za níž částice ztratí původní orientaci rychlosti. Relativně přesný výpočet srážkové frekvence elektronů s ionty lze získat ze vztahu [26] λD Ze4 ne ln Λ, Λ= . (7) νei = 2 2 3 4πε0 me v b0 Veličina Z značí stupeň ionizace plazmatu, v je vzájemná rychlost srážejících se částic a ln Λ je tzv. Coulombův logaritmus. Jde o poměr Debyeovy a Landauovy délky. Landauova délka b0 je srážkový parametr odpovídající rozptylu na úhel 90◦ . Typická hodnota Coulombova logaritmu v ideálním plazmatu je 5 − 20. V přítomnosti homogenního magnetického pole se nabité částice v plazmatu pohybují po kružnicích nebo šroubovicích s poloměrem rL a frekvencí ωc , rL =

v⊥ , ωc

ωc =

~ |q|kBk . m

(8)

~ je vektor magnetické indukce a v⊥ je kladná konstanta označující rychlost v roZde B ~ Symbolem k.k se rozumí euklidovská norma vektoru. Hodnotu rL pak vině kolmé na B. nazýváme Larmorův poloměr a ωc označujeme jako cyklotronovou frekvenci.

2.2

Způsoby popisu plazmatu

Existují tři základní přístupy k popisu plazmatu. Každý z nich má své výhody a omezení vyplývající ze zjednodušených předpokladů vhodných pouze pro určité jevy a časová měřítka [33]. Nejvytříbenější přístup k popisu plazmatu podává kinetická teorie. Ta bere v úvahu mikroskopický pohyb všech částic vyskytujících se v plazmatu. Matematicky toho lze dosáhnout úplnou časovou derivací jejich hustoty ve fázovém prostoru. S využitím faktu, 21

Kapitola 2. Laser-plazma interakce že se hustota podél trajektorií částic ve fázovém prostoru nemění, dostaneme tzv. Klimontovičovu rovnici [29]. Její tvar zde neuvádíme, protože se v praxi příliš nepoužívá. Obvykle nás totiž nezajímají přesné trajektorie jednotlivých částic a spokojíme se pouze s jistými středními hodnotami, které podávají informaci o jejich průměrném chování jakožto celku. Statistické rozdělení jednotlivých druhů částic přináší do kinetické teorie distribuční funkce fs (~x, ~v , t), kde veličiny ~x, ~v a t označují po řadě polohu, rychlost a čas. Středováním Klimontovičovy rovnice tak dostaneme základní rovnici, kterou musí splňovat distribuční funkce fs , rovnici Boltzmannovu, F~ ∂fs ∂fs + ~v · ∇fs + · = ∂t ms ∂~v

∂fs ∂t

!

.

(9)

c

Zde F~ je síla působící na částice a člen na pravé straně vyjadřuje změnu fs za jednotku času v důsledku srážek. V případě horkého plazmatu můžeme tento člen zanedbat, a je-li navíc síla F~ výhradně elektromagnetická, získáme tak nejméně přesnou, avšak nejčastěji používanou rovnici kinetické teorie plazmatu,  ∂fs qs  ~ ~ · ∂fs = 0. + ~v · ∇fs + E + ~v × B ∂t ms ∂~v

(10)

~ (~x, t) a B ~ (~x, t) je makroskopické elektrické resp. magnetické pole. Rovnice (10) se Zde E nazývá Vlasovova. V případě, kdy je nutné brát v úvahu binární interakce, je možné do této rovnice zahrnout některý ze srážkových členů. Rovnice kinetické teorie jsou propojeny se sadou Maxwellových rovnic pro elektromagnetické pole ~ ~ = − ∂B , ∇×E (11) ∂t ~ ~ = µ0 J~ + 1 ∂ E , ∇×B (12) c2 ∂t ~ = ρ, ∇·E (13) ε0 ~ =0 ∇·B (14) (kde µ0 je permeabilita vakua a c rychlost světla ve vakuu) hustotou náboje ρ (~x, t) a proudovou hustotou J~ (~x, t), ρ=

X s

qs

Z

fs d~v ,

J~ =

X

qs

Z

fs~v d~v .

(15)

s

Společně s Maxwellovými rovnicemi tvoří rovnice kinetické teorie zcela přesný, ovšem velmi komplikovaný systém, který většinou nemá analytické řešení. 22

2.2. Způsoby popisu plazmatu Druhým způsobem je hydrodynamický popis plazmatu. Jedná se o model užívající mechaniku tekutin. Spočívá tedy ve studiu pohybu elementů tekutiny nebo více prolínajících se tekutin, které ovšem obsahují elektrické náboje. Hlavní výhodou tohoto přístupu je jeho jednoduchost. V tekutinové teorii jsou závisle proměnné veličiny funkcemi pouze tří prostorových souřadnic a času. Rozdělení rychlostí považujeme za maxwellovské v každém bodě prostoru. Rovnice pro hydrodynamický popis plazmatu lze získat středováním Vlasovovy rovnice (10) přes tři první momenty rychlosti, ∂ns + ∇ · (ns~us ) = 0, ∂t "

ms ns

(16)

#

  ∂~us ~ + ~us × B ~ − ∇ps , + (~us · ∇) ~us = qs ns E ∂t

(17)

∂ 1 3 1 5 ~ s = qs ns~us · E ~ (18) ns ms u2s + ps + ∇ · ns ms~us u2s + ps~us + Q ∂t 2 2 2 2 První moment Vlasovovy rovnice dává rovnici kontinuity, která vyjadřuje zákon zachování počtu částic každého druhu, přičemž ~us (~x, t) označuje rychlost tekutiny. Druhá rovnice vyjadřuje zákon zachování hybnosti, kde ps (~x, t) je skalární tlak, který se v obecném případě nahrazuje tenzorem napětí. Třetí rovnice vyjadřuje zákon zachování energie, kde ~ s = −χ(T )∇T je tepelný tok a χ(T ) koeficient tepelné vodivosti. Je-li tento tok zaneQ dbatelný, uzavírá se systém fluidních rovnic termodynamickou stavovou rovnicí 







ps = C (ms ns )κs ,

(19)

kde κs je poměr specifických tepel Cp /CV a C libovolná konstanta. V opačném případě je nutné použít stavovou rovnici, v níž vystupuje teplota jako proměnná. Fluidní rovnice jsou provázány s Maxwellovými rovnicemi hustotou náboje ρ a proudovou hustotou J~ následovně, ρ=

X s

qs ns ,

J~ =

X

qs ns~us .

(20)

s

Dohromady podávají tyto rovnice kompletní, ale dosti hrubý popis plazmatu. Pomocí hydrodynamického modelu lze však vysvětlit většinu jevů pozorovaných v reálných experimentech. Pokud je ovšem narušena lokální termodynamická rovnováha, tento přístup není vhodný. Selhávají-li oba předchozí přístupy, je nutné použít částicový popis plazmatu. Jak už název napovídá, stav a vývoj systému je určen sledováním trajektorií jednotlivých částic. Jejich pohyb musí probíhat v souladu s Newtonovými pohybovými rovnicemi s 23

Kapitola 2. Laser-plazma interakce Lorentzovou silou,  d~v qs  ~ d~x ~ . = ~v , = E + ~v × B (21) dt dt ms Z toho plyne, že využití této metody pro větší systémy závisí na přístupu k vysokovýkonným počítačovým klastrům.

2.3

Šíření elektromagnetických vln v plazmatu

Nyní se zaměříme na to, jakým způsobem se šíří elektromagnetické vlny (tedy i laserový impuls) v plazmatu. Omezíme se na případ šíření elektromagnetické vlny v homogenním ~ = ~0 ). plazmatu (ne = konst.) bez přítomnosti srážek a vnějšího magnetického pole (B Budeme uvažovat, že je tato vlna vysokofrekvenční, a tedy můžeme ionty považovat za nehybné neutralizující pozadí. Hlavním cílem bude nalézt disperzní vztah této vlny a provést jeho analýzu. Využijeme k tomu hydrodynamický popis plazmatu. Aplikací operátoru rotace na rovnici (11) získáme   ~ ~ = −∇ × ∂ B . (22) ∇× ∇×E ∂t Zkombinujeme-li tuto s rovnicí (12) a využijeme-li identitu [13] 

~ ∇× ∇×X







~ − ∆X, ~ =∇ ∇·X

(23)

dostaneme vlnovou rovnici pro elektrické pole, ~ − ε0 µ 0 ∆E

~ ∂ 2E ∂ J~ = µ . 0 ∂t2 ∂t

(24)

Nyní provedeme linearizaci pohybové rovnice (17) pro elektrony. Tím se rozumí, že závisle proměnné veličiny rozdělíme na část rovnovážnou (označenou indexem 0) a poruchovou (označenou indexem 1), přičemž budeme uvažovat malou amplitudu oscilací, a tudíž můžeme členy obsahující vyšší řády amplitudových faktorů zanedbat, ne = ne0 + ne1 ,

~ue = ~ue0 + ~ue1 ,

~ =E ~0 + E ~1 E

(25)

Předpokládali jsme, že než dojde k vychýlení elektronů, je plazma homogenní a neutrální. Platí tedy ~0 ∂ne0 ∂~ue0 ∂E ~ 0 = 0, = = =0 (26) ∇ne0 = ~ue0 = E ∂t ∂t ∂t 24

2.3. Šíření elektromagnetických vln v plazmatu a pohybová rovnice (17) pro elektrony bude mít v přiblížení prvního řádu následující tvar, ∂~ue1 e ~ =− E 1. ∂t me

(27)

Oscilující veličiny budeme brát ve tvaru rovinných harmonických (monochromatických) vln, protože nás zajímají pouze vlastní módy (obecné řešení bychom dostali superpozicí těchto vln), ~

~ue1 ∼ ei(k·~x−ωt) ,

~

ne1 ∼ ei(k·~x−ωt) ,

~ 1 ∼ ei(~k·~x−ωt) , E

(28)

kde ~k je vlnový vektor (udává směr šíření vlny) a ω je úhlová frekvence vlny. Dosazením řešení pohybové rovnice (27) do druhého ze vztahů (20) pak snadno vyjádříme vysokofrekvenční proud, e2 ne0 ~ J~ = −ene0~ue1 = i E1 , (29) me ω kde i je imaginární jednotka. Po provedení Fourierovy transformace vlnové rovnice (24) a využitím rovnosti ε0 µ0 = 1/c dostaneme 2 ωpe ω2 ~ 2~ ~ 1. − k E1 + 2 E1 = 2 E (30) c c Odsud už snadno získáme hledanou disperzní relaci pro příčné elektromagnetické vlny ~ 1 = 0) v plazmatu, (~k · E 2 ω 2 = ωpe + c2 k 2 . (31) 2

Vztah (31) nyní podrobíme důkladné analýze. Bude-li frekvence elektromagnetické vlny ω větší než elektronová plazmová frekvence ωpe , vlnový vektor ~k bude reálný. Elektromagnetická vlna se tedy bude šířit plazmatem. Naopak bude-li ω < ωpe , pak ~k je ryze imaginární a vlna bude exponenciálně tlumena. Charakteristická vzdálenost tlumení se nazývá skinová hloubka a je definována jako δ=

1 . k~kk

(32)

Dosáhne-li frekvence elektromagnetické vlny mezní frekvence, tedy ω = ωpe , pak dochází k úplnému odrazu vlny. Tento mezní případ nastává na tzv. kritické hustotě plazmatu ε0 m e ω 2 nc = . e2

(33) 25

Kapitola 2. Laser-plazma interakce Na závěr definujeme fázovou rychlost ~vf a grupovou rychlost ~vg šíření vlny, ω ~k ~vf = , k~kk k~kk

!

~vg =

∂ω ∂ω ∂ω , , . ∂kx ∂ky ∂kz

(34)

Fázová rychlost vlny vf vyjadřuje rychlost přesunu vlnoploch a v plazmatu je vyšší než rychlost světla ve vakuu c. Grupová rychlost vln vg je mechanickou rychlostí částic k nim přidružených. Vyjadřuje tedy rychlost přenosu energie, a tudíž nemůže být vyšší než c.

2.4

Srážková absorpce

Jedním z nejdůležitějších mechanismů absorpce laserové energie v plazmatu jsou srážkové procesy dané coulombickou interakcí volných částic, nazývané také jako inverzní brzdné záření. Je tedy patrné, že jde o jev, kdy v důsledku srážek nabitých částic mezi sebou dojde k pohlcení fotonu. Pokud se elektron srazí s iontem, změní se jeho uspořádaná oscilační rychlost na tepelnou. Elektromagnetická vlna pak musí elektronu chybějící energii oscilace v poli dodat. Tímto způsobem dojde k absorpci energie vlny a její přeměně na kinetickou energii plazmatu [22]. Rigorózní odvození srážkové absorpce energie elektromagnetické vlny pomocí kinetické teorie lze nalézt v [36]. Zde se omezíme na jednoduchý hydrodynamický model homogenního plazmatu bez vnějšího magnetického pole (podobný postup lze nalézt v [7]). Ionty budeme opět považovat za nehybné. Linearizovaná pohybová rovnice pro elektrony se zahrnutím srážek je ~1 eE ∂~ue1 − νei~ue1 , (35) =− ∂t me Bereme v úvahu pouze srážky elektronů s ionty. Srážky elektronů s elektrony zanedbáme, neboť k absorpci nepřispívají. Řešení budeme hledat opět ve tvaru rovinných harmonických vln. Z pohybové rovnice (35) dostaneme výraz pro proudovou hustotu, J~ = −ene0~ue1 = i

2 ωpe ~ 1. ε0 E ω + iνei

(36)

Řešíme-li nyní vlnovou rovnici (24) pro příčné elektromagnetické vlny, obdržíme disperzní vztah 2 ω ωpe ω2 2 k = 2 − 2 . (37) c c (ω + iνei ) 26

2.4. Srážková absorpce Za předpokladu, že platí νei  ωpe (typické pro plazmatickou koronu) rozvineme disperzní vztah do Taylorovy řady. Ponecháme-li členy do prvního řádu, dostáváme 2 2 ωpe i νei ωpe ω2 2 k ≈ 2 1− 2 + . c ω ω3

!

(38)

2 2 /ω),  νei ωpe Stejným způsobem vztah odmocníme (za předpokladu, že navíc ω 2 − ωpe

ω2 ω k≈± 1 − pe c ω2

!1 " 2

νei 1+i 2ω 



2 ωpe ω2

!

#

1 . 2 /ω 2 1 − ωpe

(39)

Vidíme, že vlivem srážek bude elektromagnetická vlna tlumena. Míru útlumu intenzity vlny κib pak definujeme jako dvojnásobek imaginární části k, νei κib = 2 Im(k) = c 



2 ωpe ω2

!

1 . 2 /ω 2 1 − ωpe

(40)

Je užitečné tento vztah převést do jednotek kritické hustoty (33), νei (nc ) ne κib = c nc 

2 

ne 1− nc

− 1 2

(41)

,

kde νei (nc ) je srážková frekvence vyhodnocená na kritické hustotě. Ze vztahu (41) je okamžitě vidět, že významná část absorpce energie elektromagnetické vlny mechanismem inverzního brzdného záření nastává v blízkosti kritické hustoty nc plazmatu.

Úbytek intenzity I laserového impulsu procházejícího plazmatem ve směru osy x probíhá v souladu se vztahem dI = −κib I. (42) dx Koeficient absorpce αabs laserové energie v plazmatu o charakteristické délce L dostaneme jako   ZL Iout αabs = 1 − = 1 − exp − κib dx , (43) Iin 0

kde Iout /Iin je poměr intenzit vystupujícího a vstupujícího laserového impulsu. Pro nehomogenní plazma je řešení (42) poněkud komplikovanější. Na závěr tedy uvádíme koeficienty absorpce pro lineární a exponenciální profil hustoty plazmatu [7], αabs

32 νei (nc )L = 1 − exp − 15 c

!



pro

ne = nc

x 1− , L 

(44) 27

Kapitola 2. Laser-plazma interakce

αabs

2.5

8 νei (nc )L = 1 − exp − 3 c

!



pro

ne = nc exp −

x . L 

(45)

Landauův útlum

Jak již bylo uvedeno, elektromagnetická vlna šířící se plazmatem rozkmitává elektrony na elektronové plazmové frekvenci ωpe . Vezmeme-li v úvahu tepelné děje, zjistíme, že se tyto oscilace mohou šířit. Vznikají tak podélné elektronové plazmové vlny (bývají také často nazývány po svém objeviteli jako Langmuirovy vlny). Elektrony s rychlostí blízkou fázové rychlosti vf této vlny postupují společně s vlnou, a tak si s ní mohou účinně vyměňovat energii. Pokud je jejich rychlost nepatrně nižší než fázová, získávají energii na úkor vlny. Naopak, pokud je jejich rychlost nepatrně vyšší, energii ztrácí a předávají ji vlně. Vezmeme-li v úvahu maxwellovské rozdělení, bude v plazmatu vždy statisticky více elektronů s nižší rychlostí, díky tomu převládá proces útlumu vlny (projevuje se i bez přítomnosti srážek), jehož přímým důsledkem je produkce rychlých elektronů. Tento jev teoreticky popsal L. D. Landau v roce 1946, později byl ověřen také experimentálně. Landauův útlum není možné odvodit z hydrodynamického modelu, musíme vycházet z kinetické teorie. Cílem bude získat disperzní vztah elektronových plazmových vln. Jak již bylo uvedeno, Landauův útlum se projevuje i bez přítomnosti srážek, vyjdeme tedy z Vlasovovy rovnice (12). Nechť dále mají ionty v prostoru pevné homogenní rozložení a interakce elektronů s plazmovou vlnou zprostředkovává pouze pole elektrické. Ze sady Maxwellových rovnic nám tedy postačí pouze rovnice pro divergenci elektrického pole (13). Tyto rovnice linearizujeme, přičemž oscilující veličiny budeme brát jako harmonické rovinné vlny šířící se ve směru osy x. Snadno získáme disperzní relaci pro elektronové plazmové vlny, +∞ Z dfˆ 1 2 2 dvx , (46) k = ωpe dvx (vx − vf ) −∞

kde fˆ(vx ) symbolizuje jednorozměrnou rozdělovací funkci, již zintegrovanou přes obě zbývající složky rychlosti. Pro další pokračování je nutné vypočíst integrál komplexní proměnné (46) se singularitou v bodě vx = vf . Landau jako první ukázal, že správná integrační cesta musí vést pod pólem. Z komplexní analýzy dostaneme disperzní vztah ve tvaru 2 2 ω2 2 2 ωpe ω dfˆ 2 ω 2 ≈ ωpe + 3 pe v k + iπ ω 2 te k 2 dvx vx =vf



28

kde

2 vte =

kB Te . me

(47)

2.6. Ponderomotorická síla 2 Veličina vte označuje kvadrát tepelné rychlosti elektronů. Druhý člen v rovnosti (47) je tedy způsoben tepelnými procesy. Vystačíme si s takovou přesností, kdy tento člen budeme moci zanedbat a disperzní vztah odmocníme pomocí Taylorova rozvoje,





2 dfˆ π ωpe . ω ≈ ωpe 1 + i 2 k 2 dvx vx =vf



(48)

Je-li fˆ(vx ) jednorozměrné maxwellovské rozdělení, bude s

dfˆ vx2 2 vx . =− exp − 3 2 dvx π vte 2vte !

(49)

Dekrement Landauova útlumu γL potom získáme jako imaginární část ω (kde můžeme velmi dobře aproximovat vf ≈ ωpe /k), !

1 π ωpe γL = Im (ω) ≈ − . 3 exp − 2 (kλD ) 2 (kλD )2 r

(50)

Imaginární část ω je záporná, není tedy pochyb o tom, že se jedná o útlum. Je zřejmé, že je-li kλD  1, pak je tento útlum nepatrný a nehraje téměř žádnou roli. Pokud je ovšem vlnová délka elektronové plazmové vlny srovnatelná s Debyeovou délkou, Landauův útlum se stává velmi silným a dochází k deformaci rozdělovací funkce. Koncept bezesrážkového tlumení, jenž Landau zavedl do moderní fyziky, je v dnešní době stále živé téma. Landauův útlum hraje rozhodující roli v podstatě ve všech odvětvích fyziky plazmatu. A nejen tam, podrobným zkoumáním tohoto jevu bylo ukázáno, že ho lze využít k analýze problémů v hydrodynamice, astrofyzice nebo fyzice vysokých hustot energií [6]. Vysvětlení projevů Landauova útlumu není jednoduché, jeho výklad v různých specifických podmínkách činí stále velké obtíže. Nové a nové studie tohoto jevu se pravidelně objevují i dnes, téměř sedmdesát let po jeho objevu. Přetrvávající zájem o tento v mnoha ohledech paradoxní jev slibuje nové objevy, a to jak v teorii, tak v aplikacích [35].

2.6

Ponderomotorická síla

Laserové svazky působí na prostředí silou, která se interpretuje jako tlak záření. V plazmatu má tato síla v důsledku oscilačního pohybu nabitých částic ještě další složku, která může dosahovat velmi vysokých hodnot. Ta vytlačuje nabité částice z oblastí intenzivního pole do míst s nižší intenzitou. Sílu, která na ně tímto mechanismem působí, nazýváme ponderomotorickou silou. 29

Kapitola 2. Laser-plazma interakce Tuto nelineární sílu nyní odvodíme studiem pohybu elektronu v elektromagnetickém poli v plazmatu (podobný postup lze najít v [7]). Budeme uvažovat pouze nerelativistický případ. Vezměme elektrické pole vlny ve tvaru ~ (~x, t) = E ~ 1 (~x) cos ωt. E

(51)

Integrací Faradayova zákona (11) snadno obdržíme magnetické pole, ~ 1 (~x) sin ωt. ~ (~x, t) = − 1 ∇ × E B ω

(52)

Pohybová rovnice (21) má v přiblížení prvního řádu tvar me

d~v1 ~ 1 (~x0 ) cos ωt, = −eE dt

d~x1 = ~v1 , dt

(53)

~ (pro nerelativistické elektrony je malý) a E ~ (~x, t) vzít kde jsme mohli zanedbat člen ~v × B v počáteční poloze ~x0 . Řešením rovnic (53) dostaneme ~v1 = −

eme ~ E1 (x~0 ) sin ωt, ω

~x1 =

eme ~ E1 (x~0 ) cos ωt. ω2

(54)

Nyní vyjádříme stejnou pohybovou rovnici se členy druhého řádu, me

h i d~v2 ~ 1 (x~0 ) cos ωt + ~v1 × B ~ 1 (x~0 ) sin ωt , = −e (~x1 · ∇) E dt

(55)

~ (~x, t) jsme získali jeho rozvojem do Taylorovy řady v blízkosti počáteční kde druhý řád E polohy ~x0 . Po dosazení (52) a (54) do (55) dostaneme me

   i e2 h ~ d~v2 2 2 ~ ~ ~ =− E (~ x ) · ∇ E (~ x ) cos ωt + E (~ x ) × ∇ × E (~ x ) sin ωt . (56) 1 0 1 0 1 0 1 0 dt me ω 2

Vystředováním (56) přes jednu periodu elektromagnetické vlny (poznamenejme, že hcos2 i = hsin2 i = 1/2) a využitím identity [13] 

   ~ ·∇ X ~ +X ~ × ∇×X ~ = 1 ∇X 2 X 2

(57)

získáme efektivní nelineární sílu působící na jeden elektron, *

me

d~v2 dt

+

=−

e2 ∇E12 . 2 4me ω

(58)

Abychom získali sílu působící na jednotku objemu, musíme výraz (58) přenásobit elek30

2.7. Parametrické nestability tronovou hustotou ne . Dostaneme tzv. ponderomotorickou sílu, kterou můžeme vyjádřit pomocí ωpe (přičemž E12 = 2 hE 2 i), 2 ωpe ne e2 ε0 hE 2 i 2 = − . F~p = − ∇E ∇ 1 4me ω 2 ω2 2

(59)

Ponderomotorická síla působí stejným mechanismem i na ionty, ale jelikož je jejich hmotnost mnohem vyšší, není její okamžitý projev tak patrný. Přímým důsledkem této nízkofrekvenční síly je samofokusace laserových svazků v plazmatu. Laserový impuls vyvolá ponderomotorickou sílu v radiálním směru, která vytlačuje plazma ven ze svazku, a to pak působí jako konvexní čočka [6]. Ponderomotorická síla se také významně podílí na vzniku parametrických nestabilit při vlnových interakcích.

2.7

Parametrické nestability

K parametrickým nestabilitám dochází při nelineární interakci vln v plazmatu. Obvykle se jedná o stimulovaný rozpad původní tzv. čerpající elektromagnetické vlny na elektronové nebo iontové vlny a jinou elektromagnetickou vlnu. Útlum dceřiných vln potom způsobuje ohřev plazmatu. V některých případech se tvoří dvouteplotní rychlostní rozdělení elektronů, v důsledku čehož se narušuje lokální termodynamická rovnováha plazmatu. K parametrickým nestabilitám dochází v případě, kdy amplituda čerpající vlny převyšuje tzv. prahovou hodnotu. V opačném případě jejich vznik potlačuje srážkový nebo Landauův útlum. Rezonanční podmínky jsou pro interakci tří vln vyjádřeny dvěma následujícími vztahy, ω0 = ω1 + ω2 ,

~k0 = ~k1 + ~k2 ,

(60)

kde indexem 0 jsou označeny parametry náležející čerpající vlně, a indexy 1, 2 parametry dceřinných vln. S využitím disperzních vztahů jednotlivých druhů vln v plazmatu můžeme rozlišit několik situací, které mohou nastat v závislosti na frekvenci ω0 čerpající vlny: 1. Je-li ω0 ≈ ωpe , může dojít k rozpadu čerpající vlny na elektronovou plazmovou vlnu a iontovou akustickou vlnu. Tento jev nastává v okolí kritické hustoty nc plazmatu. 2. Je-li ω0 > ωpe , může dojít k rozpadu čerpající vlny na iontovou akustickou vlnu a jinou elektromagnetickou vlnu šířící se stejným nebo opačným směrem. Tento jev nastává na celé oblasti podkritického plazmatu a nazývá se stimulovaný Brillouinův rozptyl. 3. Je-li ω0 = 2ωpe , může dojít k rozpadu čerpající vlny na dvě elektronové plazmové vlny. Tento jev nastává na nc /4 a nazývá se dvouplazmonový rozpad. 31

Kapitola 2. Laser-plazma interakce 4. Je-li ω0 ≥ 2ωpe , může dojít k rozpadu čerpající vlny na elektronovou plazmovou vlnu a jinou elektromagnetickou vlnu, šířící se stejným nebo opačným směrem. Tento jev nastává v řídkém plazmatu o koncentraci menší nebo rovné nc /4 a nazývá se stimulovaný Ramanův rozptyl. V této práci se zaměříme na případ stimulovaného Ramanaova a Brillouinova rozptylu (detailní popis lze najít v [22]). Nejprve odvodíme rovnici vyjadřující rozptyl čerpající vlny. Budeme uvažovat elektromagnetickou vlnu šířící se homogenním plazmatem. Díky tomu, ~ (~x, t) vztahem že je magnetické pole solenoidální, můžeme zavést vektorový potenciál A ~ = ∇ × A. ~ B

(61)

Dosazením (61) do rovnice elektromagnetické indukce (11) dostaneme 



~ ~ + ∂ A  = 0. ∇ × E ∂t

(62)

~ (~x, t) potenciální není, nicméně pole (62) už ano. Lze tedy zavést skalární potenciál Pole E ϕ (~x, t) následovně ~ ~ = −∇ϕ − ∂ A . E (63) ∂t Podotkněme, že takto zavedené potenciály nejsou určeny jednoznačně a lze na ně klást dodatečné tzv. kalibrační podmínky [38]. V dalším se omezíme na coulombickou kalibraci ~ = 0). (∇ · A Dosazením (61) a (63) do Maxwellových rovnic redukujeme soustavu na dvě rovnice, ρ ∆ϕ = − , ε0

!

1 ∂2 ~ = µ0 J~ − 1 ∂ ∇ϕ. −∆ A 2 2 c ∂t c2 ∂t

(64)

Proudovou hustotu J~ (~x, t) můžeme rozdělit na příčnou J~t (~x, t) a podélnou J~l (~x, t) část. Obě složky pak získáme z pohybové rovnice (17) a rovnice kontinuity (16), 2

e ne ~ J~t = − A, me

∂ J~l = ε0 ∇ϕ. ∂t

(65)

Dosazením (65) do druhé z rovnic (64) a provedením několika triviálních úprav obdržíme vlnovou rovnici pro elektromagnetické vlny v plazmatu vyjádřenou pomocí vektorového ~ (~x, t), potenciálu A ! 2 1 ∂2 ~ = −µ0 e ne A. ~ − ∆ A (66) c2 ∂t2 me Dále budeme uvažovat malé poruchy od rovnovážného stavu a rovnici (66) linearizujeme 32

2.7. Parametrické nestability metodou perturbací podobně jako v podkapitole 2.3. V přiblížení prvního řádu dostaneme rovnici, která vyjadřuje rozptyl čerpající vlny vlivem malé fluktuace hustoty ne1 , !

∂2 2 ~ 1 = −ω 2 ne1 A ~ 0. − c2 ∆ + ωpe A pe 2 ∂t ne0

2.7.1

(67)

Stimulovaný Ramanův rozptyl

Nyní odvodíme disperzní vztah pro stimulovaný Ramanův rozptyl. Nejdříve však musíme nalézt rovnici pro fluktuaci hustoty vyvolanou elektronovou plazmovou vlnou. Využijeme hydrodynamický popis plazmatu, přičemž ionty budeme považovat za nehybné. Rozdělíme rychlost elektronové tekutiny ~ue na příčnou ~ut a podélnou ~ul část, ~ue = ~ut + ~ul =

e ~ A0 + ~ue1 . me

(68)

Užitím definic skalárního potenciálu (63), vektorového potenciálu (61) a rovnosti (68) lze linearizovat rovnici pro pohyb tekutiny (17),  2 e 1 e ~ ∇pe ∂~ue1 ~1 . = ∇ϕ1 − ∇ ~ue1 + A0 + A − ∂t me 2 me (ne0 + ne1 ) me

(69)

Poznamenejme, že druhý člen na pravé straně (69) představuje ponderomotorickou sílu. Využitím adiabatické stavové rovnice (pe /n3e = konst.) dojdeme pomocí několika triviálních algebraických úprav až k finálnímu tvaru, ∂~ue1 e e = ∇ϕ1 − ∂t me me 

2

2   ~0 · A ~ 1 − 3vte ∇ne1 . ∇ A ne0

(70)

Podobně dostaneme i linearizovanou rovnici kontinuity (16), ∂ne1 + ne0 ∇ · ~ue1 = 0 ∂t

(71)

Konečně, kombinací časové derivace rovnice (71) a divergence (70) dostaneme rovnici pro fluktuaci hustoty ne1 , !

∂2 e2 ne0  ~ ~  2 2 − 3v ∆ + ω n = − ∆ A0 · A1 . e1 te pe ∂t2 m2e

(72)

Za pomoci Fourierovy transformace rovnic (67) a (72) obdržíme výsledný disperzní 33

Kapitola 2. Laser-plazma interakce vztah pro stimulovaný Ramanův rozptyl, 



2 2 ωpe k 2 vosc 1 1 2 2 2 2   +   , ω − ωpe − 3k vte = 4 D ω − ω0 , ~k − ~k0 D ω + ω0 , ~k + ~k0

(73)

~ 0 /me označuje kde ω a k jsou parametry náležející elektronové plazmové vlně, ~vosc eA  = rychlost elektronů v poli dopadající elektromagnetické vlny a D ω, ~k = ω 2 − k 2 c2 − 2 ωpe . Důkladnou analýzou disperzního vztahu (73) a využitím rezonanční podmínky (60) získáme koeficient růstu nestability γ, − 1 ωpe k vosc  q 2 2 2 2 2 2 2 . γ= ω0 ωpe + 3k vte − ωpe − 3k vte 4

(74)

V případě zpětného stimulovaného Ramanova rozptylu (kdy je koeficient růstu γ maximální) musí být vlnové číslo 2ωpe ω0 1− k = k0 + c ω0 

2.7.2

1

2

.

(75)

Stimulovaný Brillouinův rozptyl

Nyní odvodíme disperzní vztah pro stimulovaný Brillouinův rozptyl. Nejdříve však musíme nalézt rovnici pro nízkofrekvenční fluktuaci hustoty vyvolanou iontovou akustickou vlnou. V tomto případě tedy musíme zahrnout i pohyb iontů. Pohybovou rovnici pro elektronovou tekutinu linearizujeme stejným způsobem jako v případě Ramanova rozptylu. Nyní ale použijeme izotermickou stavovou rovnici (pe = ne Te ), a vzhledem k nízkofrekvenční fluktuaci zanedbáme člen ∂~ue1 /∂t. Dostáváme, e e ∇ϕ1 = me me 

2

2   ~0 · A ~ 1 + vte ∇ne1 . ∇ A ne0

(76)

Linearizované rovnice pro iontovou tekutinu mají následující tvar, ∂ni1 + ni0 ∇ · ~ui1 = 0, ∂t

(77)

∂~ui1 Ze = − ∇ϕ1 . ∂t mi

(78)

Zkombinujeme-li časovou derivaci rovnice (77) a divergenci (78) a za ∇ϕ1 dosadíme z (76) (přičemž Zni0 = ne0 , Zni1 ≈ ne1 ), obdržíme rovnici pro nízkofrekvenční fluktuaci hustoty 34

2.7. Parametrické nestability ne1 ,

!

∂2 Ze2 ne0  ~ ~  2 ∆ A0 · A1 , − c ∆ n = e1 s ∂t2 me mi

(79)

kde cs = (ZTe /mi )1/2 je iontová akustická rychlost. Provedením Fourierovy transformace rovnic (67) a (79) obdržíme výsledný disperzní vztah pro stimulovaný Brillouinův rozptyl, 



2 2 2 ωpi k vosc 1 1 2 2 2   +   , ω − k cs = 4 D ω − ω0 , ~k − ~k0 D ω + ω0 , ~k + ~k0

(80)

q

kde ω a k jsou parametry náležející iontové akustické vlně a ωpi = ωpe Zme /mi je iontová plazmová frekvence. Koeficient růstu nestability γ získáme opět využitím rezonanční podmínky (60). V případě zpětného stimulovaného Brillouinova rozptylu je roven 1 vosc k0 ωpi γ= √ √ . 2 2 ω 0 k 0 cs

(81)

35

Kapitola 3 Numerické simulace Kolektivní chování částic a polí při interakci laserového záření s plazmatem představuje komplexní a silně nelineární problém. Vyšetřování těchto systémů nelze provádět jen prostřednictvím teoretické a experimentální práce. Veškeré analytické modely v důsledku vysokého množství mimořádně složitých simultánních interakcí s mnoha stupni volnosti selhávají nebo se uchylují ke zjednodušením či idealizaci. Experimentálně se zase nedaří zaznamenat řadu významných detailů. Proto je zapotřebí využít jiné nástroje k prohloubení znalostí v této oblasti vědeckého bádání [15]. Numerické simulace jsou v současné době nedílnou součástí základního i aplikovaného výzkumu a mají zásadní vliv na povahu a chápání fyzikálních problémů. Ruku v ruce se zvyšující se výpočetní silou moderních počítačů se dostávají na výsluní a právem se řadí mezi tradiční základní pilíře fyziky, teorii a experiment [10]. S vývojem pokročilých numerických metod vědci odkrývají mnohá tajemství v nejrůznějších oborech. Numerická schémata Newtonových rovnic mohou být implementována pro studium molekulární dynamiky, hydrodynamické simulace umožňují predikovat počasí, techniky používané pro řešení difúzních rovnic zase napomáhají kontrolovat znečištění životního prostředí [31]. Jeden počítačový program může modifikací počátečních podmínek řešit celou škálu fyzikálních problémů. Numerická simulace je často rychlejší a výrazně levnější než skutečný laboratorní experiment. Mnohdy je to také jediná možnost [23]. Simulace umožňují testovat teoreticky navržené modely a poskytují velmi podrobné informace, které lze jen stěží naměřit experimentálně. Numerické modelování se tak stává základním nástrojem při jejich vývoji a napomáhá k porozumění experimentálním výsledkům. V současné době je jasné, že detailní pochopení fyzikálních mechanismů interakce laserového záření s plazmatem může být dosaženo pouze prostřednictvím kombinace teorie, experimentu a simulace. Vývoj vhodných paralelních algoritmů, které vedou ke stabilnímu a přesnému řešení, však patří do nejnáročnějších oblastí výzkumu moderní vědy. Tato kapitola je z velké části zaměřena na detailní popis metody Particle-In-Cell. Jedná se o jeden z nejpopulárnějších numerických algoritmů používaných k simulaci dějů 37

Kapitola 3. Numerické simulace nejen v plazmatu. Bude přiblíženo matematické odvození této metody, popis jednotlivých kroků simulačního cyklu a také kritéria stability metody. Na závěr přináší tato kapitola stručné informace o výpočetním kódu LPIC++, který byl využit k simulacím v rámci této práce.

3.1

Metoda Particle-In-Cell

Označení Particle-In-Cell (PIC) reprezentuje techniku používanou k řešení určité třídy parciálních diferenciálních rovnic. Metoda byla navržena v polovině padesátých let minulého století a velké popularity si takřka okamžitě získala právě ve fyzice plazmatu [12]. Patří mezi částicové simulace, stav a vývoj systému je tedy určen sledováním trajektorií jednotlivých částic. Těch ovšem může být v reálných experimentech skutečně velmi mnoho. Naštěstí se ukázalo, že dostatečně přesný model chování plazmatu lze získat pomocí časového vývoje polohy relativně malého počtu částic konečné velikosti. Pro výrazné urychlení výpočtu se proto zavádí tzv. makročástice, které reprezentují shluk velkého množství skutečných částic v malé oblasti fázového prostoru. Otázkou ovšem je, zdali se takto zavedené makročástice budou chovat stejně, jako by se chovaly jednotlivé částice. Musíme tedy ověřit, zda splňují stejné pohybové rovnice. I přesto, že tento přístup silně redukuje počet simulačních částic, nelze brát v úvahu přímo všechny binární interakce. To si můžeme dovolit při studiu molekulární dynamiky, nikoli však plazmatu. Každý snadno zjistí, že počet operací takového algoritmu by rostl s druhou mocninou počtu makročástic, což je pro větší systémy daleko za horizontem výpočetní udržitelnosti. Při interakci intenzivního laserového impulsu s plazmatem ovšem bývá počet částic v Debyeově sféře dostatečně vysoký, a tedy dominuje kolektivní chování plazmatu. Srážky lze tedy většinou účinně zanedbat. Jednotlivé makročástice se v bezesrážkovém modelu plazmatu mohou prostupovat, proto se na krátké vzdálenosti používá redukovaný potenciál klesající k nule [24]. Pokud by nás ovšem zajímal účinek coulombických srážek, lze je implementovat odděleně jinými metodami. Poloha a rychlost nabitých částic jednoznačně určuje stav systému v libovolném časovém okamžiku. Abychom zjistili, jak se bude dále vyvíjet, je nejvhodnější využít kinetickou teorii. Hydrodynamický popis plazmatu, založený na nerelativistické dynamice a lokální termodynamické rovnováze, je v těchto případech většinou odsouzen k neúspěchu. Pravděpodobnost nalezení částice druhu s v elementu fázového prostoru je dána distribuční funkcí fs (~x, ~v , t), která se v případě bezesrážkového plazmatu řídí Vlasovovou rovnicí (10). Předpokládejme, že celkovou hustotu pravděpodobnosti těchto částic obdržíme superpozicí distribučních funkcí pro makročástice fs (~x, ~v , t) =

X p

38

fp (~x, ~v , t) .

(82)

3.1. Metoda Particle-In-Cell Vezměme nyní pro názornost v úvahu pouze jednu složku polohy a rychlosti. Tvar distribuční funkce jedné makročástice získáme tenzorovým součinem fp (x, v, t) = Np Sx (x − xp (t)) Sv (v − vp (t)) ,

(83)

kde Np je počet fyzikálních částic, které jsou reprezentovány makročásticí v daném objemu fázového prostoru, Sx a Sv jsou její tvarové funkce polohy a rychlosti a xp (t), vp (t) je poloha resp. rychlost makročástice v čase t. Tvarové funkce nelze volit libovolně. Jejich vlastnosti podléhají několika speciálním kritériím. Nechť Sξ je tvarová funkce libovolné fázové proměnné ξ. Potom platí: 1. Nosič tvarové funkce je kompaktní, ∃R > 0, supp Sξ ⊂ (−R, R) .

(84)

2. Integrál tvarové funkce je roven jedné, +∞ Z

Sξ (ξ) dξ = 1.

(85)

Sξ (ξ) = Sξ (−ξ) .

(86)

−∞

3. Tvarová funkce je symetrická,

Přestože tyto podmínky ponechávají stále relativně širokou škálu možností pro volbu tvarových funkcí, tradičně se v PIC simulacích volí jenom z několika jednoduchých funkčních tříd. Tvarová funkce rychlosti se standardně volí jako Diracova δ-funkce Sv (v − vp (t)) = δ (v − vp (t)) .

(87)

Tato volba vyjadřuje skutečnost, že všechny částice obsažené v makročástici mají stejnou rychlost, a tedy během časového vývoje systému zůstávají poblíž sebe. Jako tvarové funkce prostorových souřadnic se dříve používaly také δ-funkce, v současné době se však téměř vždy volí křivky b-spline. B-spline nultého řádu je definován jako ( 1 pro |ξ| < 1/2 b0 (ξ) = . (88) 0 jinak Vyšší řády získáme iterativně pomocí konvoluce, bm (ξ) = (b0 ∗ bm−1 ) (ξ) =

+∞ Z

b0 (ξ − η) bm−1 (η) dη.

(89)

−∞

39

Kapitola 3. Numerické simulace Tvarovou funkci prostorových souřadnic pak definujeme jako !

Sx (x − xp (t)) = bm

x − xp (t) , ∆p

(90)

kde ∆p je velikost nosiče makročástice. Volbou vyšších řádů pro tvarové funkce lze dosáhnout snížení šumu a redukce nefyzikálních jevů v simulacích, samozřejmě na úkor zvýšení výpočetního času. Ve většině současných PIC kódů se nejčastěji používají b-spliny prvního řádu, tedy lineární funkce. Zbývá dodat, že rozšíření distribuční funkce pro makročástici do více dimenzí provedeme triviálně součinem tvarových funkcí pro každou souřadnici,  

Sξ ξ~ = Sξ (ξ1 ) Sξ (ξ2 ) Sξ (ξ3 ) .

(91)

Nyní můžeme odvodit rovnice pro časový vývoj parametrů ~xp , ~vp a Np určujících distribuční funkce makročástic. Díky linearitě Vlasovovy rovnice v distribuční funkci mají rovnice pro hustotu pravděpodobnosti jednotlivých makročástic stejný tvar. Volbou libovolných tvarových funkcí ovšem nedosáhneme zcela přesného řešení. Obvykle stačí, aby bylo splněno pouze několik prvních momentů [19]. Rovnice pro nultý moment a první momenty v prostorové a rychlostní části fázového prostoru mají následující tvar, ZZ

ZZ ZZ 1 ~ ∂ ∂fp F· d~x d~v + ~v · ∇fp d~x d~v + ∂t ms ∂~v

!

fp d~x d~v = 0,

ZZ

ZZ ZZ ∂fp ~x ~ ∂ d~x d~v + ~x (~v · ∇) fp d~x d~v + F· ~x ∂t ms ∂~v

!

ZZ

ZZ ZZ ∂fp ~v ~ ∂ ~v d~x d~v + ~v (~v · ∇) fp d~x d~v + F· ∂t ms ∂~v

!

(92)

fp d~x d~v = 0,

(93)

fp d~x d~v = 0.

(94)

Zde F~ (~x, ~v , t) nahrazuje výraz pro Lorentzovu sílu. S využitím vlastností tvarových funkcí dostáváme úpravou těchto rovnic následující vztahy, dNp = 0, dt

d~xp = ~vp , dt

d~vp F~p = . dt ms

(95)

První výraz vyjadřuje zákon zachování počtu částic reprezentovaných každou jednotlivou makročásticí, zbylé dvě rovnice jsou shodné s Newtonovými pohybovými rovnicemi. Důvodem toho je fakt, že poměr náboje ku hmotnosti je při takto zavedené transformaci invariantní. Zásadní rozdíl je ovšem v tom, že na makročástice působí střední Lorentzova síla   ~ p (t) + ~vp (t) × B ~ p (t) , F~p (t) = qs E (96) 40

3.1. Metoda Particle-In-Cell kde ~ p (t) = E

Z

~ (~x, t) Sx (~x − ~xp ) d~x, E

~ p (t) = B

Z

~ (~x, t) Sx (~x − ~xp ) d~x. B

(97)

Srdcem PIC metody je cyklus, který je znázorněn na obr. 1. Jeho jednotlivé kroky jsou detailně popsány v několika následujících sekcích. Volba časového a prostorového kroku, na níž závisí stabilita metody, je diskutována v závěru této podkapitoly. Integrátor částic F~p → (~x, ~v ) p



Váhování polí  ~ B ~ E, → F~p

Váhování částic   (~x, ~v ) → ρ, J~ p

i,j,k

Integrátor polí    ~ B ~ ρ, J~ → E,

i,j,k



i,j,k

i,j,k

Obr. 1: Základní cyklus metody Particle-In-Cell

3.1.1

Integrátor polí

Připomeňme si sadu Maxwellových rovnic popisující nestacionární elektromagnetické pole buzené pohybujícími se částicemi ve vakuu (11), (12), (13), (14). Metoda PIC ovšem za účelem výrazné redukce požadavků na výpočetní výkon řeší tyto rovnice jen v určitých bodech simulační oblasti. Je tedy nutné provést diskretizaci prostoru souřadnic ~x → xi,j,k , kde (i, j, k) ∈ Z3 . Pro názornost zavedeme pravoúhlou výpočetní síť s konstantní vzdáleností mezi jednotlivými uzly. V případě trojrozměrné kartézské soustavy souřadnic tedy xi,j,k = (i∆x, j∆y, k∆z). Makročástice se skrz takto definovanou výpočetní síť mohou volně pohybovat. Zajímá nás časový vývoj polí, tudíž využijeme některý z modelů pracujících v časové oblasti. Jednou z nejpopulárnějších výpočetních metod pro elektromagnetické problémy je metoda konečných diferencí v časové doméně. K diskretizaci Maxwellových rovnic v čase použijeme schéma Leap-Frog. Je explicitní, velmi rychlé a navíc časově reverzibilní. Schéma počítá se složkami elektrického a magne~ je definován tického pole v různých časových krocích. Vektor intenzity elektrického pole E n n+1/2 ~ v časech t = n∆t, vektor magnetické indukce B v časech t = (n + 1/2)∆t, n ∈ N. 41

Kapitola 3. Numerické simulace Díky tomu nabývá toto schéma druhého řádu přesnosti. Maxwellovy rovnice (11), (12), (13), (14) přecházejí na následující tvar, ~ n+1/2 − B ~ n−1/2 B ~ n, = −∇− × E ∆t

~ n+1 − E ~n 1E ~ n+1/2 − µ0 J~ n+1/2 , = ∇+ × B c2 ∆t

ρn n ~ ∇ ·E = , ε0 +

~ n+1/2 = 0. ∇− · B

(98) (99)

Zde diskrétní operátory (∇+ ) a (∇− ) působí na skalární pole fi,j,k způsobem !

fi+1,j,k − fi,j,k fi,j+1,k − fi,j,k fi,j,k+1 − fi,j,k , , , ∆x ∆y ∆z

+

∇ fi,j,k = −

∇ fi,j,k =

fi,j,k − fi−1,j,k fi,j,k − fi,j−1,k fi,j,k − fi,j,k−1 , , ∆x ∆y ∆z

(100)

!

(101)

a mají následující vlastnosti, ∇− · ∇− × = ∇+ · ∇+ × = 0,

∇− · ∇+ = ∇+ · ∇− = ∆.

(102)

Symbol ∆ označuje diskrétní Laplaceův operátor v centrálních diferencích, fi−1,j,k + 2fi,j,k + fi+1,j,k fi,j−1,k + 2fi,j,k + fi,j+1,k fi,j,k−1 + 2fi,j,k + fi,j,k+1 + + . ∆x2 ∆y 2 ∆z 2 (103) Před samotným řešením Maxwellových rovnic je nutné si uvědomit, že tato soustava není nezávislá. V trojrozměrném případě se jedná o osm diferenciálních rovnic prvního řádu, ovšem pouze se šesti neznámými vektorovými komponentami. Aplikací operátoru (∇− ·) resp. (∇+ ·) na první resp. druhou z rovnic (98) dostáváme s využitím záměny pořadí derivování ~ n+1/2 − ∇− · B ~ n−1/2 ∇− · B = 0, (104) ∆t ρ n+1 − ρ n + ∇+ · J~ n+1/2 = 0. (105) ∆t To znamená, že stačí vyřešit rovnice (98), přičemž na divergenční rovnice (99) lze nahlížet jako na počáteční podmínky v čase t = 0. V tomto případě ale musí být splněna diskrétní rovnice kontinuity (105). Metoda konečných prvků v časové doméně definuje jednotlivé komponenty vektoru intenzity elektrického pole ve středech stěn buňky a složky vektoru magnetické indukce ve středech hran buňky, ∆fi,j,k =





2 3 ~ = E1 E i,j+1/2,k+1/2 , Ei+1/2,j,k+1/2 , Ei+1/2,j+1/2,k ,

42

(106)

3.1. Metoda Particle-In-Cell 



2 3 ~ = B1 B i+1/2,j,k , Bi,j+1/2,k , Bi,j,k+1/2 .

(107)

Složky proudové hustoty jsou definovány ve stejných bodech jako složky vektoru intenzity elektrického pole a hustota náboje je definována uprostřed buňky, 



1 2 3 J~ = Ji,j+1/2,k+1/2 , Ji+1/2,j,k+1/2 , Ji+1/2,j+1/2,k ,

(108)

ρ = ρi+1/2,j+1/2,k+1/2 .

(109)

Toto schéma, známé jako Yee síť [42], se ukázalo být velmi robustní. Ilustraci jedné kartézské Yee výpočetní buňky používané v PIC kódech lze najít na obr. 2.

(i + 1, j + 1, k + 1)

(i, j + 1, k + 1)

(i + 1, j, k + 1)

(i, j, k + 1)

Ex Bz

Ey (i, j + 1, k)

Ez

By (i, j, k)

(i + 1, j + 1, k)

(i + 1, j, k)

Bx

Obr. 2: Kartézská Yee výpočetní buňka používaná v PIC kódech

3.1.2

Váhování částic a polí

Nyní se zaměříme na interpolaci zdrojových členů Maxwellových rovnic na výpočetní síť a následně zpětnou interpolaci hodnot elektrického a magnetického pole na místa, kde jsou lokalizovány makročástice. Hustota náboje v libovolném bodě výpočetní sítě je definována jako průměrná hodnota v jeho okolí o velikosti jedné buňky. V jedno-dimenzionálním případě tedy 1 ρi (t) = ∆x

xi +∆x/2 Z

ρ (x, t) dx.

(110)

xi −∆x/2

43

Kapitola 3. Numerické simulace Využijeme předpisu pro hustotu náboje pomocí tvarových funkcí prostorových souřadnic jednotlivých makročástic, ρ (x, t) =

X

qp Sx (x − xp (t)) ,

qp = qs Np .

(111)

p

Charakteristickou funkci množiny, přes kterou integrujeme, můžeme napsat jako b-spline nultého řádu s nosičem o velikosti buňky ∆x. Dostáváme +∞   x − xi 1 X Z b0 qp Sx (x − xp (t)) dx. ρi (t) = ∆x p ∆x

(112)

−∞

Nyní můžeme snadno definovat interpolační funkci Wi (xp (t)) =

+∞ Z



b0

−∞

x − xi Sx (x − xp (t)) dx, ∆x 

(113)

kterou lze opět rozšířit do více dimenzí součinem interpolačních funkcí pro každou souřadnici,   (114) Wi,j,k ξ~ = Wi,j,k (ξ1 ) Wi,j,k (ξ2 ) Wi,j,k (ξ3 ) . Hustotu náboje ρi,j,k (t) v libovolném bodě výpočetní sítě xi,j,k tedy můžeme vyjádřit bez nutnosti numerické integrace, ρi,j,k (t) =

1 X qp Wi,j,k (~xp (t)) , ∆V p

(115)

kde ∆V = ∆x∆y∆z je objem buňky. Obdobně získáme i proudovou hustotu J~i,j,k (t), 1 X J~i,j,k (t) = qp~vp (t) Wi,j,k (~xp (t)) . ∆V p

(116)

Nutno podotknout, že tento triviální způsob interpolace proudové hustoty nesplňuje přesně diskrétní rovnici kontinuity (105). Proto se v PIC kódech používají pokročilejší metody, jejich popis lze najít např. v [40], [8] nebo [39]. Po vyřešení Maxwellových rovnic na výpočetní síti je nutné elektrické a magnetické pole zespojitit. Pro libovolnou složku např. elektrického pole (opět pro názornost v jedné dimenzi) toho docílíme následovně, E (x, t) =

X i

44



Ei (t) b0

x − xi . ∆x 

(117)

3.1. Metoda Particle-In-Cell Stejný vztah bude samozřejmě platit i pro pole magnetické. Předpokládá se tedy, že jsou jednotlivé komponenty polí v okolí o velikosti jedné buňky konstantní a rovné průměrné hodnotě. Dále využijeme rovnost (97) a dostaneme střední hodnoty polí působících na makročástice. Libovolnou složku vektoru intenzity elektrického pole Ep a vektoru magnetické indukce Bp (nyní již v obecném případě) získáme snadno, Ep (t) =

X

Ei,j,k Wi,j,k (~xp (t)) ,

Bp (t) =

i,j,k

X

Bi,j,k Wi,j,k (~xp (t)) .

(118)

i,j,k

Zbývá dodat, že je důrazně doporučováno využívat stejné interpolační funkce pro interpolaci elektromagnetických polí a zdrojových členů Maxwellových rovnic. Jsou pro to dobré důvody. Používání stejných interpolačních funkcí eliminuje silové působení makročástic na sebe samé (tedy zamezuje tomu, aby se jednotlivé makročástice samy urychlovaly) a zajišťuje zachování hybnosti [9].

3.1.3

Integrátor částic

Zjistili jsme, že časový vývoj volných parametrů ~xp a ~vp distribuční funkce pro makročástice je určen Newtonovými pohybovými rovnicemi (95) se střední Lorentzovou silou (96). Při interakci intenzivního laserového svazku s plazmatem se ovšem makročástice obvykle pohybují s rychlostmi blízkými rychlosti světla ve vakuu c, provedeme tedy relativistické zobecnění rychlostí pomocí faktoru γ, ~up = γ ~vp ,

γ=

v u u t

~up 1+ c

!2

(119)

.

Dostáváme následující tvar rovnic, ~up d~xp = , dt γ

!

d~up qs ~ ~up ~p . = Ep + ×B dt ms γ

(120)

K diskretizaci pohybových rovnic (120) v čase použijeme opět schéma Leap-Frog. Nyní budou definovány složky polohy makročástice ~xp v časech t n a složky zobecněné rychlosti ~up v časech t n+1/2 . Po nahrazení spojitých veličin diskrétními hodnotami tedy obdržíme následující tvar rovnic, ~xpn+1 − ~xpn ~u n+1/2 = pn+1/2 , (121) ∆t γ ~upn+1/2 − ~upn−1/2 qs ~ n ~upn+1/2 + ~upn−1/2 ~n . = Ep + ×B p n ∆t ms 2γ !

(122)

Přestože tyto rovnice vypadají velmi jednoduše, jejich řešení představuje časově nejnáročnější část výpočtu. Rovnice je totiž nutné řešit pro každou makročástici zvlášť. 45

Kapitola 3. Numerické simulace Jakákoliv optimalizace výpočtu tedy výrazně redukuje simulační čas. V PIC algoritmech je v současné době de facto standardem elegantní Borisova metoda [4], která separuje účinek elektrického a magnetického pole. J. P. Boris využil toho, že magnetické pole způsobuje pouze rotaci vektoru rychlosti, nikoliv změnu jeho velikosti. Definujeme-li ~upn−1/2

=

~up−

~ n ∆t qs E p − , ms 2

~upn+1/2

=

~up+

~ n ∆t qs E p + ms 2

(123)

a následně tyto vztahy dosadíme do (122), efektivně se zbavíme elektrického pole,  ~up+ − ~up− qs  + − ~ pn . ×B + ~ u = ~ u p ∆t 2γ n ms p

(124)

Z rovnice (124) nyní okamžitě vidíme, že během jednoho časového kroku ∆t dojde v rovině ~ n k otočení ~u − k ~u + . Důležité je si uvědomit, že velikost úhlu θ, který kolmé na pole B p p p oba vektory ~up− a ~up+ svírají, závisí na cyklotronové frekvenci ωc způsobem θ = ωc ∆t.

(125)

Nejprve tedy přidáme k původnímu vektoru ~upn−1/2 polovinu účinku elektrického pole ~ n . Ze ~ n , provedeme rotaci o úhel θ a nakonec přidáme zbývající druhou polovinu účinku E E p p základní geometrie Boris odvodil následující posloupnost vztahů, vedoucí až k hledanému vektoru ~upn+1/2 . Z první z rovnic (123) vyjádříme vektor ~up− , ~up− = ~upn−1/2 +

~ pn ∆t qs E ms 2

(126) 



a využijeme ho ke konstrukci pomocného vektoru ~up0 , který bude kolmý na ~up+ − ~up− a ~ pn zároveň, B ~ n ∆t qB ~t = s p ~up0 = ~up− + ~up− × ~t, . (127) 2γ n ms ~ pn a jeho velikost je Vektor ~t, který zde vystupuje, musí být logicky rovnoběžný s B ~ pn ) rovnotan (θ/2) ≈ θ/2 pro malé úhly. Dále využijeme toho, že je nyní vektor (~up0 × B  běžný s ~up+ − ~up− a vyjádříme ~up+ , ~up+ = ~up− + ~up0 × ~s, 46

~s =

2~t . 1 + t2

(128)

3.1. Metoda Particle-In-Cell ~ n a je normován tak, aby byla splněna podmínka k~u + k = Zde vektor ~s je rovnoběžný s B p p − − + k~up k. Přechod od ~up k ~up lze názorněji zapsat po složkách pomocí matice, 

~up+



1 − s2 t2 − s3 t3 s2 t1 + s3 s3 t1 − s2  1 − s1 t1 − s3 t3 s3 t2 + s1  =  s1 t2 − s3 up− . ~ s1 t3 + s2 s2 t 3 − s1 1 − s1 t1 − s2 t2

(129)

Vektor ~up+ nakonec dosadíme do druhé z rovnic (123), ~upn+1/2

=

~up+

~ n ∆t qs E p . + ms 2

(130)

Nyní už nám nic nebrání k tomu, abychom získali složky vektoru polohy makročástic v následujícím časovém kroku, ~xpn+1 = ~xpn +

3.1.4

~upn+1/2 ∆t, γ n+1/2

γ n+1/2 =

v   u n+1/2 2 u ~up u  . t1 + 

c

(131)

Stabilita metody

Stabilita standardní PIC metody přímo závisí na velikosti časového kroku a buněk výpočetní sítě. Velikost buněk výpočetní sítě musí být dostatečně malá, abychom byli schopni rozlišit Debyeovu stínící délku. V opačném případě dochází k tzv. aliasingu neboli vyhlazování různých Fourierových módů. To vede k nestabilitám a většinou také k numerickému ohřevu [19]. Velikost časového kroku musí být taková, aby bylo umožněno šíření elektromagnetických i elektrostatických vln bez ohledu na jejich relevanci na měřítku, které nás zajímá. To znamená, že imaginární část jejich frekvence musí být nulová. To je splněno pouze, pokud platí 1 1 1 1 > + + . (132) 2 2 2 2 c ∆t ∆x ∆y ∆z 2 Nerovnice (132) vyjadřuje tzv. Courantovu (nebo Courant–Fridrichs–Levyho, CFL) podmínku ve třech dimenzích [15]. Tato podmínka limituje rozsah pohybu makročástic během jednoho časového kroku. Zajišťuje, aby jakýkoliv objekt pohybující se rychlostí světla ve vakuu c nemohl urazit vzdálenost přesahující velikost jedné buňky. Respektuje tedy omezení uložená během řešení pohybových rovnic. Na druhou stranu, volba příliš malého časového kroku vede také ke změnám disperzního vztahu, a tedy k nefyzikálním výsledkům. Proto se doporučuje volit jeho velikost tak velkou, jak jen to umožňuje Courantova podmínka. 47

Kapitola 3. Numerické simulace Stručně řečeno, standardní PIC metoda musí být schopna postihnout procesy odehrávající se na nejmenších možných časových a prostorových škálách.

3.2

Výpočetní kód LPIC++

LPIC++ [25] je výpočetní Particle-in-Cell kód zaměřený na simulaci interakce laserového záření s plazmatem. Byl vyvinut v Max–Planck–Institut für Quantenoptik v německém Garchingu. Program je založen na relativistickém elektromagnetickém algoritmu, což znamená, že řeší celou sadu Maxwellových rovnic. Je jedno-dimenzionální v prostorové souřadnici, ale obsahuje všechny tři složky rychlosti, tudíž lze modelovat i kruhově polarizované laserové svazky. Dále byl obohacen o simulaci pružných coulombických srážek metodou Monte Carlo. Kód je napsán v programovacím jazyce C++ a je objektově orientován, aby byl snadno rozšiřitelný. Paralelizace a komunikace mezi procesory je zajištěna pomocí knihovny MPI a její implementace spočívá v automatickém rozdělení simulační oblasti na řadu sousedících domén tak, aby byl v každé přibližně stejný počet částic. Díky tomu klesá simulační čas téměř lineárně s hardwarovou výpočetní silou. Kód LPIC++ lze využít ke studiu absorpce laserového záření v plazmatu, produkce rychlých elektronů, šíření laserového svazku v podkritickém plazmatu, generace magnetického pole nebo také inerciálně držené fúze. Program používá bezrozměrné jednotky. Lze se tak vyhnout práci s fyzikálními konstantami, což vede k lepšímu výkonu simulačního kódu. Výsledky jsou navíc obecnější tím, že nejsou vázané na konkrétní geometrii. Prostorová a časová souřadnice jsou vztaženy k vlnové délce λ resp. periodě τ laserového impulsu, x t t 0 := . (133) x 0 := , λ τ Amplituda elektrického pole a intenzita laserového impulsu jsou udávány v následujících jednotkách, eE0 a0 = , Iλ2 = a20 · 1, 37 · 1018 Wµm2 /cm2 , (134) me ωc kde ω je frekvence laserového impulsu. Rychlosti jsou normovány rychlostí světla ve vakuu c, ~v ~v 0 := . (135) c Hmotnost a náboj částic jsou udávány v jednotkách hmotnosti elektronu me resp. elementárního náboje e, q m m 0 := , z := . (136) me e 48

3.2. Výpočetní kód LPIC++ Nábojové a proudové hustoty jsou vztaženy ke kritické hustotě nc , ρ 0 :=

ρ , nc e

~j 0 :=

~j . enc c

(137)

Elektrické a magnetické pole dostáváme v následujících tvarech, ~ ~ 0 := eE , E me ωc

~ ~ 0 := eB . B me ω

(138)

Jednotky vektorového a skalárního potenciálu mají tuto podobu, ~ ~ 0 = eA , A me c

ϕ0 =

eϕ . me c2

(139)

49

Kapitola 4 Výsledky V této kapitole budou představeny výsledky dosažené v rámci této práce. Jedná se o implementaci okrajové podmínky pro efektivnější absorpci rychlých elektronů ve výpočetním kódu LPIC++ a testování její funkčnosti na modelových příkladech. Dále pak o PIC simulace interakce laserového impulsu s plazmatem pro podmínky současných experimentů Badatelského centra PALS představujícího jódový fotodisociační laserový systém pracující v blízké infračervené oblasti na základní vlnové délce. Tyto výsledky budou porovnány a zasazeny do kontextu současného výzkumu inerciální fúze zapálené silnou rázovou vlnou.

4.1

Okrajová podmínka

Jakékoliv ukončení simulační oblasti je nefyzikální. Tuto skutečnost je nutné mít neustále na mysli. Speciálně pak při delších simulacích je potřeba jí věnovat zvýšenou pozornost. Elektrony urychlené při interakci laserového impulsu s plazmatickou koronou prostupují skrz terč a zcela určitě dosáhnou okraje simulační oblasti. V případě, že dovolíme, aby se elektron odrazil zpět, vrátí se do interakční oblasti a ovlivňuje zde distribuční funkci, což může ve výsledku vést k významné změně fyzikálních jevů. Pokud zvolíme variantu absorpce elektronu, dochází na okraji simulační oblasti k hromadění náboje, a tedy generaci vysokého elektrického pole, které opět ovlivňuje výsledek simulací (obr. 3-a). V ideálním případě by z důvodu zachování kvazineutrality měly být rychlé elektrony nahrazeny zpětným proudem tepelných elektronů s přibližně počáteční teplotou plazmatu. V kódu LPIC++ byla implementována okrajová podmínka metodou Monte Carlo. Před startem simulace je zvolena funkce (hustota pravděpodobnosti), jejíž hodnota závisí na čísle výpočetní buňky. V každém kroku cyklu se při pohybu částice generují náhodná čísla od 0 do 1. Pokud je toto číslo menší než hodnota funkce v dané buňce, dojde k tomu, že se rychlost částice nahradí tepelnou rychlostí. Tento proces má nahradit termalizaci 51

Kapitola 4. Výsledky (a)

(b) 12

f(x) = 1/2 − 1/2cos(πx/3000) 0.8

3

0.75

2

0.5 1

0.25

0.6 0.4 0.2

0

0 0

1

f(x)

ni/nc (−)

1

ni Ex

E (V/m)

1.25

x 10 4

25

50

75

x (µm)

100

0 0

125

500

1000

1500

2000

pocet bunek (−)

2500

3000

Obr. 3: (a) Profil iontové hustoty (modře) a podélné elektrické pole (červeně) v modelové simulaci interakce intenzivního laserového impulsu s plazmatem bez okrajové podmínky. Svazek dopadá zleva, nc je kritická hustota. (b) Graf funkce pro okrajovou podmínku nastavenou na posledních 3000 buněk, která byla využita v simulacích. Testováno pro případ, kdy na jednu vlnovou délku laseru připadá 120 výpočetních buněk. (a)

(b)

1.4

0 1000 2000

namerenéhhodnoty

1

ni/nch(−)

ni/nc (−)

1.3

analytickáhaproximace

1.2

1.2 1.1

0.8 0.6 0.4

1 0.2 0.9 0

500

1000

pocet bunek (−)

1500

2000

500

1000

1500

2000

pocethbunekh(−)

2500

3000

Obr. 4: (a) Porucha v profilu iontové hustoty způsobená vlivem okrajové podmínky. Modře je označena varianta bez využití okrajové podmínky, červeně resp. zeleně je zobrazen profil hustoty v případě okrajové podmínky nastavené na posledních 1000 resp. 2000 výpočetních buněk. (b) Exponenciální závislost absolutní hodnoty poruchy v hustotním profilu na délce okrajové podmínky.

částice v hustém plazmatu dále v terči, kterou však z důvodu výpočetní náročnosti nejsme schopni simulovat. Vzhledem k tomu, že tímto způsobem implementace okrajové podmínky měníme vlastnosti prostředí na relativně malé vzdálenosti, dochází zde ke vzniku poruchy hustoty šířící se zpět do míst s nižší hustotou. Příklad takovéhoto narušení hustotního profilu terče lze 52

4.2. Simulace (a)

(b) do terce zpet

pocet elektronu (−)

3

10

2

10

1

10

0

10

4

do terce zpet

10

pocet elektronu (−)

4

10

3

10

2

10

1

10

0

200

400

600

800

E (keV)

1000

1200

1400

10

0

200

400

600

800

E (keV)

1 000 1 200 1 400

Obr. 5: (a) Případ, kdy okrajová podmínka nefunguje. Do terče se vrací rychlé elektrony a způsobují nefyzikální jevy v simulaci. (b) V tomto případě okrajová podmínka funguje a do terče se vrací pouze elektrony s výrazně nižší energií.

vidět na obr. 4-a. Velikost poruchy lze minimalizovat rozložením okrajové podmínky na větší počet výpočetních buněk (obr. 4-b). Tím docílíme toho, že se vlastnosti prostředí budou měnit pozvolněji. V průběhu práce byly otestovány různé lineární, exponenciální a goniometrické funkce. V simulacích byla nakonec použita modifikovaná goniometrická funkce (obr. 3-b), jež poskytovala nejlepší výsledky. Okrajová podmínka byla testována pro kolmo dopadající p-polarizovaný laserový svazek, který byl charakterizován parametrem Iλ2 = 4, 95·1016 Wµm2 /cm2 , přičemž na jednu vlnovou délku připadalo 120 výpočetních buněk. Počáteční teplota pro elektrony byla 530 eV, pro ionty 60 eV. Coulombické srážky nebyly v simulaci uvažovány. Na závěr ještě uvádíme ukázku distribučních funkcí elektronů pohybujících se směrem do terče a zpět pro případ, kdy okrajová podmínka nefunguje (obr. 5-a) a naopak, kdy lze její chování považovat za bezchybné (obr. 5-b).

4.2

Simulace

V rámci této práce byly provedeny jedno-dimenzionální simulace pomocí výpočetního kódu LPIC++ pro základní vlnovou délku λ = 1,315 µm laserového systému PALS a vlnovou délku λ = 0,438 µm odpovídající třetí harmonické základní frekvence. Jako vstup těchto simulací byly použity výsledky dvou-dimenzionálních jednotekutinových hydrodynamických výpočtů poskytnutých Ing. M. Holcem a Ing. M. Kuchaříkem, Ph.D. Počáteční podmínky pro PIC simulace byly zvoleny z hydrodynamických simulací v čase 50 ps před maximem impulsu uprostřed laserového svazku. Profily elektronové hustoty plazmatu a teploty pro první i třetí harmonickou v tomto čase jsou znázorněny na obr. 6. Profil hustoty plazmatu pro první harmonickou byl aproximován exponenciální 53

Kapitola 4. Výsledky (a)

(b)

0.8

0.6

2

0.4

1

0 0

50

100

150

200

x (µm)

250

300

350

400

1 T ne aproximace

0.8

3

0.6

2

0.4

0.2

1

0.2

0

0 0

T (keV)

3

4

ne/nc (−)

T (keV)

4

5

1 ne aproximace T

50

100

150

200

250

ne/nc (−)

5

0

x (µm)

Obr. 6: (a) Profil elektronové hustoty a teploty z hydrodynamických výpočtů pro první harmonickou. (b) Profil elektronové hustoty a teploty z hydrodynamických výpočtů pro třetí harmonickou.

funkcí s charakteristickou délkou 90 µm. Teplotu jsme uvažovali konstantní, rovnou 4,5 keV pro elektrony i ionty zároveň. Pro třetí harmonickou byla hustota elektronů aproximována exponenciálou s charakteristickou délkou 50 µm a teplotu jsme brali rovněž jako konstantu, ovšem s hodnotou 3 keV. Energie impulsu byla v hydrodynamických simulacích 400 J pro první resp. 200 J pro třetí harmonickou. Intenzita laserového impulsu pro první harmonickou byla 1·1016 W/cm2 , což odpovídá maximální intenzitě v ohnisku soustředěného paprsku laseru PALS. Intenzitu pro třetí harmonickou jsme zvolili poloviční v souladu s 55% účinností konverze hlavního svazku na třetí harmonickou, tedy 5 · 1015 W/cm2 . Ačkoli je délka impulsu laserového systému PALS 350 ps, zvolili jsme délku impulsu v simulaci 100 ps. Naší snahou bylo totiž modelovat interakci jen pro vysoké intenzity kolem maxima laserového impulsu. Simulační čas byl tedy nastaven na 100 ps pro obě simulace, velikost simulační oblasti byla zvolena podle hustotních profilů na 800 µm v případě první harmonické, 350 µm v případě třetí harmonické. Velikost buněk byla konstantní, rovná 19 nm pro první harmonickou a 5 nm pro třetí harmonickou. Obě hodnoty jsou srovnatelné s minimální Debyeovou délkou. Obě simulace obsahují pouze jeden druh iontů, jejichž stupeň ionizace Z = 3,5 a dále zahrnují coulombické srážky s konstantním Coulombovým logaritmem Λ = 8. Další důležité parametry simulací jsou uvedeny v příloze A, kompletní vstupní soubory pro výpočetní kód LPIC++ lze nalézt na přiloženém CD. Prahová hodnota nestabilit stimulovaných laserovým zářením se zvyšuje při zkracování vlnové délky. Lze proto očekávat, že parametrické nestability budou hrát mnohem významnější roli v simulaci první harmonické. Zda-li tomu tak skutečně je zjistíme z obr. 7, kde jsou zobrazeny spektrální hustoty rozptýleného záření v závislosti na čase pro první a třetí harmonickou. Oba grafy si vyžadují důkladný rozbor. Záření na frekvenci 54

4.2. Simulace (b)

100

−3

100

80

−4

80

60

−5

40

−6

40

20

−7

20

0 0

0.25

0.5

0.75

ω/ω0 (−)

1

1.25

1.5

−8

t (ps)

t (ps)

(a)

−4

−5

60 −6

0 0

−7

0.25

0.5

0.75

ω/ω0 (−)

1

1.25

−8

Obr. 7: (a) Časová závislost spektra rozptýleného záření v logaritmickém měřítku pro první harmonickou. (b) Časová závislost spektra rozptýleného záření v logaritmickém měřítku pro třetí harmonickou. Na barevné stupnici jsou vyneseny exponenty normované spektrální hustoty.

vstupujícího laserového svazku ω0 odpovídá odrazu na kritické hustotě a stimulovanému Brillouinovu rozptylu. Závislost frekvence Brillouinova rozptylu na hustotě se řídí následujícím vztahem [20], s

ω = ω0 − 2k0 (cs + u),

ω0 ne k0 = 1− , c nc

(140)

kde u je rychlost expanze plazmatu. Odsud plyne, že je frekvenční posuv stimulovaného Brillouinova rozptylu způsoben frekvencí iontové akustické vlny a Dopplerovým jevem. Z obr. 7 je dále zjevný výskyt záření na polovině frekvence vstupujícího laserového impulsu ω0 , což je stimulovaný Ramanův rozptyl. Dále vidíme záření na čtvrtině ω0 , které lze interpretovat jako sekundární Ramanův rozptyl na již odražené elektromagnetické vlně. Stimulovaný Ramanův rozptyl je v tomto případě nestabilita absolutní. To je způsobeno vysokou počáteční teplotou, která vede k silnému Landauovu útlumu plazmových vln. Absolutní Ramanův rozptyl je prostorově lokalizován, což je patrné z obr. 7-b. V dané oblasti pak dochází k významným změnám parametrů plazmatu. Z tohoto obrázku je také nápadná závislost míry rozptýleného záření na okamžité intenzitě laserového impulsu. Na obr. 7-a je s očekáváním zaznamenán silný frekvenční šum, který lze vysvětlit předchozími úvahami a vysokou intenzitou vstupujícího laserového svazku na základní frekvenci. Dekrement Landauova útlumu není dostatečně vysoký, aby potlačil vznik parametrických nestabilit a dochází zde i k rozsáhlým sekundárním rozptylům. Obr. 8 zachycuje vývoj iontové hustoty plazmatu v čase. V důsledku vysoké počáteční teploty plazma postupně expanduje, ale jeho charakteristická délka se v průběhu simulace významně nemění. Je důležité si uvědomit, že v čase kolem 50 ps nabývá intenzita vstupujícího laserového impulsu svého maxima. Po tomto okamžiku dochází ke snížení 55

Kapitola 4. Výsledky (a)

(b) 140

100

70

100

120

60 80

100

60

80 60

40

t (ps)

t (ps)

80

50

60

40 30

40

40 20 0 0

20 20

20 100

200

300

400

500

0 0

0

600

10 50

100

x (µm)

150

200

250

0

x (µm)

Obr. 8: (a) Časová závislost iontové hustoty pro první harmonickou. (b) Časová závislost iontové hustoty pro třetí harmonickou. Barevná stupnice udává procenta kritické hustoty. (a)

(b)

1.2

1.2 ni 1

0.8

0.8

ne/nc (−)

ne/nc (−)

ni 1

0.6

0.6

0.4

0.4

0.2

0.2

0 0

100

200

300

x (µm)

400

500

600

0 0

50

100

150

200

250

300

x (µm)

Obr. 9: (a) Hustotním profil v čase 50 ps pro první harmonickou. (b) Hustotní profil v čase 50 ps pro třetí harmonickou.

ponderomotorické síly, a tedy je expanze plazmatu výraznější. Na čtvrtině kritické hustoty pak vznikají vlivem Ramanova rozptylu v obou případech, první i třetí harmonické, poruchy v hustotním profilu neboli kavity. V případě simulace první harmonické vznikají tyto kavity též na kritické hustotě (obr. 8-a). Kritická hustota plazmatu je v případě první harmonické přibližně 6, 46 · 1020 částic na cm3 , v případě třetí harmonické 5, 82 · 1021 částic na cm3 . Ještě lépe jsou vidět poruchy v profilu elektronové hustoty na obr. 9. Ten zachycuje hustotní profil v čase 50 ps, kdy dochází k největšímu růstu kavit. Zdůrazněme, že tyto poruchy nemusí plně odpovídat realitě, protože kavitace je trojrozměrný jev [20]. V kavitách se zachycuje určitá část elektromagnetického pole, což lze nahlédnout z 56

4.2. Simulace (a)

(b) −3

0.1

100

0.08

60

0.06

40

0.04

20

0.02

0 0

100

200

300

400

500

4

80

t (ps)

t (ps)

80

x 10 5

100

60

3

40

2

20

1

0 0

0

50

x (µm)

100

150

200

0

x (µm)

Obr. 10: (a) Hustota elektromagnetického pole v logaritmickém měřítku v závislosti na čase a prostoru pro první harmonickou. (b) Hustota elektromagnetického pole v logaritmickém měřítku v závislosti na čase a prostoru pro třetí harmonickou. Na barevné stupnici jsou vyneseny bezrozměrné jednotky používané kódem LPIC++. (a)

(b) 1

incident total 1ω SRS 1ω

0.8

normovaná intenzita (−)

normovaná intenzita (−)

1

0.6 0.4 0.2 0 0

20

40

60

t (ps)

80

100

incident total 3ω SRS 3ω

0.8 0.6 0.4 0.2 0 0

20

40

60

80

100

t (ps)

Obr. 11: (a) Závislost intenzity vstupujícího laserového impulsu na čase (zeleně) a jeho celková reflektivita (modře) resp. reflektivita způsobená stimulovaným Ramanovým rozptylem (červeně) pro první harmonickou. (b) Závislost intenzity vstupujícího laserového impulsu na čase (zeleně) a jeho celková reflektivita (modře) resp. reflektivita způsobená stimulovaným Ramanovým rozptylem (červeně) pro třetí harmonickou.

obr. 10. Takto zachycené záření nemůže v jedno-dimenzionálních simulacích uniknout a je postupně konvertováno na kinetickou energii elektronů. Tímto způsobem tedy také dochází k absorpci energie vstupujícího laserového záření. Okamžitě vidíme, že větší množství elektromagnetického pole je v kavitách zachyceno v případě simulace první harmonické (obr. 10-a). Na obr. 11 je zobrazena reflektivita laserového impulsu. Pro první harmonickou činí 57

Kapitola 4. Výsledky (a)

(b)

4

4

x 10

4 t = 0 ps T = 4.5 keV t = 50 ps T = 5.1 keV

1.5

t = 0 ps T = 3 keV t = 50 ps T = 4.2 keV

3.5

pocet elektronu (−)

pocet elektronu (−)

2

x 10

1

0.5

3 2.5 2 1.5 1 0.5

0 0

5

10

15

20

25

E (keV)

30

0 0

35

5

10

15

20

E (keV)

25

Obr. 12: (a) Energetické spektrum všech elektronů na začátku simulace (modře) a v čase 50 ps (červeně) pro první harmonickou. (b) Energetické spektrum všech elektronů na začátku simulace (modře) a v čase 50 ps (červeně) pro první harmonickou. Byla provedena aproximace hodnot pomocí Maxwell-Boltzmannovy distribuční funkce s teplotou uvedenou v legendě. (a)

(b) 4

10

pocet elektronu (−)

4

pocet elektronu (−)

10

0 ps 50 ps

3

10

2

10

0 ps 50 ps

3

10

2

10

1

10

1

10 −3

−2

−1

0

px / (mec) (−)

1

2

3

−0.6

−0.4

−0.2

0

0.2

px / (mec) (−)

0.4

0.6

0.8

Obr. 13: Rozdělovací funkce hybnosti ve směru osy x všech elektronů na počátku simulace (modře) a v čase 50 ps (červeně) pro první harmonickou. (b) Rozdělovací funkce hybnosti ve směru osy x všech elektronů na počátku simulace (modře) a v čase 50 ps (červeně) pro třetí harmonickou.

celková reflektivita 43 %, reflektivita způsobená stimulovaným Ramanovým rozptylem kolem 10 %. Celková reflektivita pro třetí harmonickou je výrazně vyšší, téměř 59 %, reflektivita způsobená stimulovaným Ramanovým rozptylem je v tomto případě 26 %. Na obr. 11-b je dále markantní závislost množství odraženého záření na okamžité intenzitě laserového impulsu, což je v souladu s obr. 7-b. Podíváme-li se na obr. 12-a, zjistíme, že ačkoli bylo v simulaci pro první harmonickou absorbováno zhruba 57 % energie laserového svazku, srážková absorpce nehraje 58

4.2. Simulace (a)

(b)

5

4

10 do terce zpet 170 keV 5 keV

4

10

pocet elektronu (−)

pocet elektronu (−)

10

3

10

2

10

1

10

200

400

600

E (keV)

800

1000

do terce zpet 25 keV 3.2 keV

3

10

2

10

1

10

20

40

60

E (keV)

80

100

120

Obr. 14: (a) Energie elektronů jdoucích směrem do terče (modře) a zpět (červeně) v místě začátku okrajové podmínky pro první harmonickou. (a) Energie elektronů jdoucích směrem do terče (modře) a zpět (červeně) v místě začátku okrajové podmínky pro třetí harmonickou.

příliš velkou roli. Tímto mechanismem se zvýšila teplota elektronů pouze o 0.6 keV. To lze vysvětlit vysokou intenzitou laserového impulsu a vyšší teplotou. Zároveň tato skutečnost ukazuje na přítomnost nelineárních bezesrážkových disipačních procesů. Všechny tyto procesy jsou ovšem silně selektivní, jinak řečeno, energie je přenášena do relativně malé skupiny rychlých elektronů. Naopak, v případě třetí harmonické (obr. 12-b) srážková absorpce dominuje. Její účinnost je dvojnásobná, absorbovalo se kolem 41 % energie laserového impulsu a teplota plazmatu se zvýšila o 1,2 keV. Vznik rychlých elektronů v případě první harmonické je markantní z rozdělovací funkce na obr. 13-a. V případě simulace třetí harmonické vzniká mnohem menší množství rychlých elektronů (obr. 13-b). Na obr. 14 jsou vykresleny energie elektronů procházejících přes myšlenou hranici v místě, kde začíná okrajová podmínka. Do terče se vrací elektrony s teplotou téměř shodnou s tou, která byla zvolena na počátku v případě obou simulací. Okrajová podmínka tedy funguje podle našich požadavků. Podle simulací [20] se zdálo, že teplota rychlých elektronů příliš nezávisí na intenzitě vstupujícího laserového impulsu a pohybuje se v řádu několika desítek keV. Nedávné studie navíc ukazují, že tyto rychlé elektrony mohou dokonce samy vyvolat silnou rázovou vlnu [11]. Míra předehřátí stlačovaného paliva způsobená těmito elektrony v pozdější fázi komprese při inerciální fúzi by nemusela být významná, protože největší podíl na transportu energie do terče mají elektrony s teplotou vyšší než 100 keV, kterých je relativně málo. Vzhledem k vysokým ztrátám energie při konverzi do třetí harmonické by tato skutečnost mohla kompletně změnit úhel pohledu na parametrické nestability v kontextu inerciální fúze zapálené silnou rázovou vlnou. 59

Kapitola 4. Výsledky Tyto příznivé výsledky se však pro 1. harmonickou (obr. 14-a) laseru PALS nepotvrdily. Podíváme-li se na teplotu rychlých elektronů směřujících do terče, zjistíme, že je v tomto případě příliš vysoká, lze ji odhadnout na 170 keV. V tomto případě by značná část rychlých elektronů pravděpodobně prošla stlačenou slupkou palivového terče, což je pro zapálení inerciální fúze rázovou vlnou nevhodné. Tato vlnová délka je již příliš dlouhá při intenzitě 1 · 1016 W/cm2 , to vede k zapojení procesů, které generují více energetické elektrony. Vysokou teplotu pak lze vysvětlit několika způsoby. Teplota je většinou dána fázovou rychlostí plazmových vln a ta závisí na vlnovém čísle. V případě dopředného Ramanova rozptylu je vlnové číslo nižší, a tedy je fázová rychlost plazmových vln relativně vysoká. Dopředný Ramanův rozptyl má sice menší rychlost růstu, ale to při této intenzitě nemusí hrát velkou roli. Druhou možností je elektromagnetické pole zachycené v kavitách, ~ 0 je 0,58) ukazují na skutečnost, že se jedná jehož vysoké hodnoty (maximální hodnota A ~ o relativistickou interakci, při níž má nemalý význam urychlování silou ~v × B. Teplota rychlých elektronů, směřujících do terče v případě třetí harmonické je zhruba 25 keV, což je v souladu s výsledkem simulací [20]. Tyto a předchozí výsledky jsou rovněž důležité pro interpretaci experimentů probíhajících na laseru PALS.

60

Závěr Předkládaná práce se zabývá studiem fyzikálních procesů při interakci intenzivních laserových impulsů s plazmatem pro podmínky současných experimentů studujících možnosti zapálení inerciální fúze rázovou vlnou. Během této práce byla úspěšně implementovaná okrajová podmínka pro efektivnější absorpci rychlých elektronů ve výpočetním kódu LPIC++. Její důkladné testovaní na modelových simulacích interakce intenzivního laserového impulsu s plazmatem prokazuje její funkčnost. Dále byly provedeny jedno-dimenzionální simulace svazku laserového systému PALS na základní vlnové délce 1,315 µm s intenzitou 1 · 1016 W/cm2 a na vlnové délce 0,438 µm odpovídající třetí harmonické základní frekvence s intenzitou 5 · 1015 W/cm2 . Oba impulsy trvaly 100 ps. K simulacím byl využit výpočetní kód LPIC++ založený na metodě Particle-in-Cell. Počáteční profily hustoty a teploty byly zvoleny z hydrodynamických výpočtů. Celková reflektivita pro první harmonickou činila 43 %, přičemž většina energie byla absorbována hustotními kavitami a parametrickými nestabilitami. Tyto procesy formují dvouteplotní rychlostní rozdělení a způsobují vznik rychlých elektronů. Celková reflektivita pro třetí harmonickou byla výrazně vyšší, téměř 59 %, ovšem dominantním mechanismem absorpce byly coulombické srážky mezi elektrony a ionty. Vznik hustotních poruch a parametrických nestabilit byl potlačen, produkce rychlých elektronů byla výrazně nižší. Teplota rychlých elektronů v případě simulace první harmonické byla odhadnuta na 170 keV, pro třetí harmonickou byla zhruba 25 keV. Je nutné dále zkoumat, jaký vliv by tyto elektrony měly při inerciální fúzi, nicméně v případě 1. harmonické by značná část těchto elektronů pravděpodobně prošla stlačenou slupkou palivového terče, což je pro zapálení inerciální fúze rázovou vlnou nevhodné.

61

Poděkování/Acknowledgements Děkuji Ing. Ondřeji Klimovi, Ph.D. za trpělivé vedení mé práce, nesčetné množství konzultací k tématu a předání mnoha praktických dovedností, které jsou pro mě velkým přínosem a příslibem do budoucna. Dále děkuji Ing. Milanu Holcovi a Ing. Milanu Kuchaříkovi Ph.D za poskytnutí hydrodynamických výpočtů. V neposlední řadě děkuji své rodině za neutuchající podporu v průběhu celého dosavadního studia, které si velmi vážím a cením. Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the programme ”Projects of Large Infrastructure for Research, Development, and Innovations”(LM2010005), is greatly appreciated. Access to the CERIT-SC computing and storage facilities provided under the programme Center CERIT Scientific Cloud, part of the Operational Program Research and Development for Innovations, reg. no. CZ. 1.05/3.2.00/08.0144, is greatly appreciated. Petr Valenta 62

Literatura [1] Atzeni, S. – Meyer-ter-Vehn, J.: The Physics of Inertial Fusion: Beam Plasma Interaction, Hydrodynamics, Hot Dense Matter. Oxford science publications, Clarendon Press, 2004, ISBN 0191524050. [2] Betti, R. a kol.: Shock ignition: A new approach to high gain inertial confinement fusion on the National Ignition Facility. Phys. Rev. Lett., roč. 103, č. 4, 2009, ISSN 1079-7114. [3] Bin, J. a kol.: A laser-driven nanosecond proton source for radiobiological studies. Appl. Phys. Lett., roč. 101, č. 24, 2012, ISSN 0003-6951. [4] Birdsall, C. K. – Langdon, A. B.: Plasma Physics via Computer Simulation. Series in Plasma Physics, CRC Press, 2004, ISBN 0750310251. [5] Bishop, B.: NIF experiments show initial gain in fusion fuel. 2014. URL [6] Chen, F. F.: Úvod do fyziky plazmatu. Academia, 1984, ISBN (neuvedeno). [7] Eliezer, S.: The Interaction of High-Power Lasers with Plasmas. Series in Plasma Physics, CRC Press, 2010, ISBN 9781420033380. [8] Esirkepov, T. Z.: Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor. Computer Physics Communications, roč. 135, č. 2, 2001, ISSN 0010-4655. [9] Fehske, H. – Schneider, R. – Weiße, A.: Computational Many-Particle Physics. Springer, 2008, ISBN 9783540746850. [10] Gould, H. – Tobochnik, J. – Christian, W.: An introduction to computer simulation methods: applications to physical systems. Addison–Wesley series in physics, Addison–Wesley, třetí vyd., 2007, ISBN 1201165032. 63

[11] Guskov, S. a kol.: Ablation Pressure Driven by an Energetic Electron Beam in a Dense Plasma. Phys. Rev. Lett., roč. 109, č. 25, 2012, ISSN 1079-7114. [12] Hockney, R. W. – Eastwood, J. W.: Computer Simulation Using Particles. CRC Press, 2010, ISBN 9781439822050. [13] Huba, J. D.: NRL Plasma Formulary. NRL publication, Naval Research Laboratory, 1994, ISBN 9781234775360. [14] Hurricane, O. A. a kol.: Fuel gain exceeding unity in an inertially confined fusion implosion. Nature, roč. 506, č. 7488, 2014, ISSN 1476-4687. [15] Jaroszynski, D. A. – Bingham, R. – Cairns, R. A.: Laser-Plasma Interactions. Scottish Graduate Series, CRC Press, 2009, ISBN 1584887788. [16] Kálal, M.: Generace rychlých elektronů při interakci intenzivního laserového záření s plazmatem. Pokroky matematiky, fyziky a astronomie, roč. 25, č. 4, 1980, ISSN 0032-2423. [17] Klimo, O.: Fyzika inerciální fúze. 2. lekce – Termojaderná fúze a její udržení. URL [18] Klimo, O.: Fyzika inerciální fúze. 12. lekce – Fast Ignition. URL [19] Klimo, O.: Metody počítačové fyziky. 8. lekce – Řešení Vlasovovy rovnice – Numerical methods for the Vlasov equation. URL [20] Klimo, O. – Tikhonchuk, V. T.: Laser–plasma interaction studies in the context of shock ignition: the regime dominated by parametric instabilities. Plasma Physics and Controlled Fusion, roč. 55, č. 9, 2013, ISSN 1361-6587. [21] Klimo, O. a kol.: Particle-in-cell simulations of laser–plasma interaction for the shock ignition scenario. Plasma Physics and Controlled Fusion, roč. 52, č. 5, 2010, ISSN 1361-6587. [22] Kruer, W. L.: The Physics of Laser-Plasma Interactions. Frontiers in physics, Westview Press, 2003, ISBN 9780813340838. [23] Kuchařík, M.: Metody počítačové fyziky. Užití počítačů ve fyzice. URL 64

[24] Lapenta, G.: Particle In Cell Methods With Application to Simulations in Space Weather. Lecture notes. URL [25] Lichters, R. – Meyer-ter-Vehn, J. – Pfund, R.: LPIC++: A Parallel Onedimensional Relativistic Electromagnetic Particle-In-Cell Code for Simulating Laser– Plasma–Interaction. Max-Planck-Institut für Quantenoptik, MPQ, 1997. [26] Limpouch, J.: Základy fyziky plazmatu. 1. lekce – Úvod. URL [27] Mašek, M.: Eulerova Vlasovova metoda pro laserové plazma. Dizertační práce, Univerzita Karlova v Praze, 2006. [28] McCracken, G. – Stott, P.: Fusion: The Energy of the Universe. Complementary science series, Elsevier Academic Press, druhé vyd., 2005, ISBN 0123846579. [29] Nicholson, D. R.: Introduction to Plasma Theory. Plasma Physics Series, Wiley, 1983, ISBN 9780471090458. [30] Řípa, M. a kol.: Řízená termojaderná fúze pro každého – 4U. Materiály pro nové tisíciletí, Reklamní agentura Daniel s.r.o. Most, čtvrté vyd., 2013, ISBN 978-80-2604785-8. [31] Pang, T.: An introduction to computational physics. Cambridge University Press, druhé vyd., 2006, ISBN 0521825695. [32] Pfalzner, S.: An Introduction to Inertial Confinement Fusion. Series in Plasma Physics, CRC Press, 2006, ISBN 1420011847. [33] Pšikal, J.: Ion Acceleration in Small-size Targets by Ultra-intense Short Laser Pulses (Simulation and Theory). Dizertační práce, Czech Technical University in Prague, Université Bordeaux 1, 2009. [34] Rohlena, K.: Inerciální fúze. Čs. čas. fyz, roč. 60, č. 6, 2010, ISSN 1804-8536. [35] Ryutov, D. D.: Landau damping: half a century with the great discovery. Plasma Physics and Controlled Fusion, roč. 41, č. 3, 1999, ISSN 1361-6587. [36] Simon, A. – Thompson, W. B.: Advances in Plasma Physics. An Interscience series, Wiley, 1976, ISBN 9780471791928. [37] Theobald, W. a kol.: Initial experiments on the shock-ignition inertial confinement fusion concept. Physics of Plasmas, roč. 15, č. 5, 2008, ISSN 1089-7674. 65

[38] Štoll, I.: Elektřina a magnetismus. Vydavatelství ČVUT, druhé vyd., 2003, ISBN 80-01-02693-0. [39] Umeda, T. a kol.: A new charge conservation method in electromagnetic particlein-cell simulations. Computer Physics Communications, roč. 156, č. 1, 2003, ISSN 0010-4655. [40] Villasenor, J. – Buneman, O.: Rigorous charge conservation for local electromagnetic field solvers. Computer Physics Communications, roč. 69, č. 2-3, 1992, ISSN 0010-4655. [41] Weinzettl, V.: Čistá energie tokamaků. Vesmír, roč. 77, č. 4, 1998, ISSN 12144029. [42] Yee, K.: Numerical solution of inital boundary value problems involving maxwell’s equations in isotropic media. IEEE Transactions on Antennas and Propagation, roč. 14, č. 3, 1966, ISSN 0018-926X.

66

Přílohy

67

Příloha A Vstupní soubory V příloze A jsou vypsány nejdůležitější parametry simulací pro první a třetí harmonickou ze vstupních souborů pro výpočetní kód LPIC++. Kompletní soubory i s detailním popisem parametrů lze nalézt na přiloženém CD.

A.1 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Simulace první harmonické

&pulse_front -------------------------Q = 1 amplitude = 0.1123 amplitude2 = 0.0 amplitude3 = 0.0 phase2 = 0 phase3 = 0 angle = 0 polarization = 2 shape = 3 delay = 0 raise = 11410 duration = 22820 pulse_save = 0 pulse_save_step = 0.02

19 20 21

&propagate -------------------------prop_start = 0 prop_stop = 25000

26 27 28 29 30 31 32 33 34

&box

-------------------------cells_per_wl = 68 cells = 42780 cells_left = 10000 cells_plasma = 32780 cells_ramp = 22780 n_ion_over_nc = 0.6 scale_length = 68.5 boundary_cond = 1 emin = 1000 box_save = 1

35 36 37 38 39

&electrons -------------------------fix = 0 ppc = 140 vtherm = 0.0938

41 42 43 44 45

22 23

25

40

17 18

24

46

&ions -------------------------fix = 0 z = 3.5 m = 8568

69

47 48

ppc = 140 vtherm = 0.001 &collisions

A.2 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

19 20 21

&pulse_front -------------------------Q = 1 amplitude = 0.0265 amplitude2 = 0.0 amplitude3 = 0.0 phase2 = 0 phase3 = 0 angle = 0 polarization = 2 shape = 3 delay = 0 raise = 34230 duration = 68460 pulse_save = 0 pulse_save_step = 0.02

24 25 26 27

28 29 30 31 32 33 34

&propagate -------------------------prop_start = 0 prop_stop = 75000

cells_plasma = 61000 cells_ramp = 46000 n_ion_over_nc = 0.372 scale_length = 114.2 boundary_cond = 1 emin = 1000 box_save = 1

35 36 37 38 39 40

&electrons -------------------------fix = 0 ppc = 120 vtherm = 0.0766

41 42 43 44 45 46 47 48

&ions -------------------------fix = 0 z = 3.5 m = 8568 ppc = 120 vtherm = 0.0008272

49

22 23

54

-------------------------Q = 1 lambda = 1.3e-6 coulomb_log = 8.0

Simulace třetí harmonické

17 18

52 53

49 50

51

&box -------------------------cells_per_wl = 94 cells = 76000 cells_left = 15000

70

50 51 52 53 54

&collisions -------------------------Q = 1 lambda = 0.438e-6 coulomb_log = 8.0

Příloha B Obsah CD adresář/soubor

popis

bp.pdf

tato práce

lpic plain

Přeložený kód LPIC++ určený k běhu na jednom výpočetním jádře

lpic parallel

Přeložený kód LPIC++ určený k běhu na více výpočetních jádrech

input 1w.lpi

vstupní soubor kódu LPIC++ pro simulaci první harmonické

input 3w.lpi

vstupní soubor kódu LPIC++ pro simulaci třetí harmonické

src

adresář se zdrojovým kódem LPIC++ obsahující finální úpravy

bash

skripty pro unixový shell Bash sloužící ke spouštění úloh

python

skripty napsané v jazyce Python sloužící pro práci se soubory

matlab

skripty napsané v jazyce MATLAB sloužící převážně k vizualizaci dat

71