IN
DEGREE PROJECT ELECTRICAL ENGINEERING, SECOND CYCLE, 30 CREDITS
,
STOCKHOLM SWEDEN 2020
Implementation of Collisions in
FEMIC for Modelling of the RF
Heating with a Lower Hybrid
Resonance
LÉO BELLEIL
KTH ROYAL INSTITUTE OF TECHNOLOGY
Implementation of Collisions
in FEMIC for Modelling of the
RF Heating with a Lower
Hybrid Resonance
LÉO BELLEIL
Master in Electromagnetics, Fusion and Space Engineering, Date: December 18, 2020
Supervisor: Thomas Jonsson and Björn Ljungberg Examiner: Jan Scheffel
School of Electrical Engineering and Computer Science, Division of Fusion Plasma Physics
Swedish title: Implementering av kollisioner i FEMIC för
iii
Abstract
Fusion energy is a potential sustainable solution to provide clean energy for the human society. The fusion reaction requires that the ionised fuel (plasma) is heated to extreme temperatures. Magnetic confinement is used to keep the plasma away from the wall and the components mounted on it. The heating method that will be discussed in this thesis is the Ion Cyclotron Resonance Heating (ICRH). A wave is launched into the plasma with a frequency equal to the ion cyclotron frequency of one of the ions presents in the plasma. The wave transfers its energy to the plasma at the location where the two frequencies match.
This thesis concerns the modelling of RF wave propagation with a lower hybrid resonance using the FEMIC code (Finite Element Method for Ion Cy-clotron resonance heating). We aim to model the propagation of the slow mag-neto sonic wave (SW) in a specific region of a fusion device, the scrape off-layer (SOL). This region is located between the core plasma and the wall and is characterised by a low plasma density.
However, the modelling of the SW in this region is difficult numerically because the Lower Hybrid Resonance (LHR) is expected here. There is a pole (singularity) in the dispersion relation and the wave vector diverges at this location. The corresponding wavelength would then be infinitely small. The introduction of a new quantity, the collision frequency, can solve this problem by moving the pole to the imaginary plane. This collision frequency can also be increased artificially to increase the wavelength of the SW. It can be useful to reduce the computational power required for the modelling of the SW, since the mesh size can be increased as well. There is however a problem related to this method, an artificial damping is induced by the introduction of the artificial collision frequency.
In this thesis, the possibility of having a lower hybrid resonance in the scrape-off layer has been assessed for specific values of the density and tem-perature. Moving the pole into the complex plane using an artificial collision frequency has been proven to be successful while aiming to reduce the val-ues of the perpendicular wave number. The damping related to this artificial collision frequency limits the artificial decrease of the wave number (if the damping is too strong the wave is absorbed before the resonance).
The results of this thesis are encouraging since the propagation of the SW towards the LHR have been successfully modelled in ITER. With sufficient computational power, it is possible to expand the domain of study to the entire SOL and to resolve the SW at the LHR.
iv
Sammanfattning
Energi från fusionsreaktioner är en kandidat till att förse framtida samhällen med en hållbar och ren energi. Fusionsreaktionen kräver dock att bränslet är joniserat, dvs i form av ett plasma, och att det värms upp till extrema tempe-raturer. Bränslet kan inneslutas med hjälp av ett magnetfält som håller plas-mat borta från en kringliggande vägg, samt komponenter som är monterade på väggen. Detta arbete behandlar en metod för att värma plasmat som kallas joncyklotronuppvärmning (ICRH), vilket bygger på att en våg skickas in i plas-mat med en frekvens som motsvarar cyklotronfrekvensen hos ett av jonslagen i plasmat. När dessa två frekvenser är lika absorberas vågens energi av plasmat. I denna uppsats beskrivs modellering med FEMIC-koden (Finite Element Method for Ion Cyclotron resonance heating) av hur radiovågor propagerar när det finns en så kallad lägre-hybridresonans i plasma. Målet är att nume-riskt modellera hur en långsam magnetosonisk våg (SW) propagerar i den del av en fusionsanläggning som kallas avskapningslagret (SOL). Denna region ligger mellan huvuddelen av plasmat och väggen, och karakteriseras av låg plasmatäthet.
Numerisk modellering av en SW i denna region är dock speciellt utmanan-de om utmanan-det finns en lägre-hybridresonans i regionen. Denna resonans iutmanan-dentifie- identifie-ras med en pol i dispersionsrelationen, varvid vågtalet divergerar. Motsvarande våglängd är oändligt kort. Problemet kan lösas genom att introducera kollisio-ner, vilket flyttar polen ut i det komplexa planet. I numerisk modellering kan kollisionsfrekvensen förstärkas artificiellt för att förlänga våglängden ytterli-gare. Detta kan vara ett effektivt sätt att minska antalet noder i det numeriska systemet och därmed minska de beräkningsmässiga kraven. Detta kan med-föra problem då den förstärkta kollisionsfrekvensen ger upphov till artificiell dämpning av vågen.
I denna uppsats har möjligheten att ha en läghybridresonance i SOL re-gionen utvärderas för ett fall med givna täthets- och temperaturprofiler. Att flytta polen ut i det komplexa planet genom en artificiellt hög kollisionsfre-kvens har visat sig vara en bra metod för att reducera de våglängder som ska upplösas. Men, den dämpning som de artificiella kollisionerna orsakar kan innebära att vågen dämpas innan den når den lägre-hybridresonansen.
Resultaten visar att man kan beskriva numeriskt hur SW propagerar mot den lägre-hybridresonansen i en mindre region av fusionsanläggningen ITER. Med tillräckliga beräkningsresurser kommer det att vara möjligt att modellera SW i en större del av avskrapningslagret, samt att lösa upp SW vid den lägre-hybridresonansen.
Contents
1 Introduction 1
1.1 The world’s energy production . . . 1
1.2 The fusion reaction . . . 2
1.2.1 The tokamak . . . 2
1.2.2 Heating processes . . . 4
1.2.3 The scrape-off-layer . . . 5
2 Wave and plasma physics 7 2.1 The wave equation . . . 7
2.2 The dielectric tensor in a hot plasma . . . 7
2.3 Dispersion relation . . . 10
2.3.1 Lower Hybrid resonance frequency . . . 12
2.4 The LHR in FEMIC . . . 12
3 Collisions 14 3.1 Theory . . . 14
3.1.1 Calculation of the Coulomb logarithm . . . 15
3.2 Implementation in the Dielectric Tensor . . . 15
4 Modelling in FEMIC 17 4.1 Meshing . . . 17
4.2 Tilting of the magnetic field . . . 18
5 Results 20 5.1 JET . . . 20
5.1.1 Collisions . . . 21
5.1.2 The dispersion relation . . . 22
5.2 ITER . . . 25
5.2.1 The dispersion relation . . . 26
5.2.2 Modelling of the slow wave in FEMIC . . . 28
vi CONTENTS
5.2.3 Tilting of the antenna . . . 31
6 Discussion 32
7 Conclusions 34
A Wave equation 37
B Third component of the dielectric tensor 38 C Lower hybrid resonance: estimation of the density 42 D Validity of the collisional model 44
Chapter 1
Introduction
1.1
The world’s energy production
The society is nowadays facing new challenges. One of them is to find a viable and clean way to produce energy. The current global energy consumption cor-responds approximately to 1800 TWh annually [1], and it is going to increase with the access to technology of developing and underdeveloped countries. The main problem is the pollution caused by the energy production. Indeed, 85 percent of the energy produced comes from fossil fuels (Figure 1.1). Burn-ing fossil fuels emits a lot of CO2in the atmosphere, increasing the greenhouse
effect, which leads to global warming. Other alternatives to produce energy need to be found to satisfy the increasing consumption.
Figure 1.1: The world power production per fuel in percentage 2019 [1]
One such alternative is nuclear fusion. It has been developed since shortly
2 CHAPTER 1. INTRODUCTION
after World War 2 and is now approaching commercialisation with projects like ITER (International Thermonuclear Experimental Reactor).
Fusion energy has many advantages over other sources of energy. Indeed, the location of the power plant is not that much restricted by the environment, which means that it can provide a local supply of energy and reduce losses due to transport. In terms of safety, there is no risk of an uncontrolled chain reaction because the fusion process is hard to maintain. Therefore, the reaction will stop if anything unusual happen. Moreover, fusion radioactive waste is much less of a problem than in a fission reactor.
1.2
The fusion reaction
For a fusion reaction to occur two particles of same charge needs to collide with each other’s. However, the Coulomb force repel two particles of same charge. The Coulomb barrier is the minimum energy that two particles need to overpass to collide. If they gain enough energy, they will come close enough to each other so that the nuclear force become stronger than the Coulomb force. To gain energy the particles are heated up to temperature exceeding 100 million ° K, so high that the electrons leave the orbit of ions and move freely. The resulting ensemble of charged particles is quasi-neutral and shows a collective behaviour, it is a plasma. Since all the electrons leave the ion it is called a fully ionised plasma. During a fusion reaction, the reactants fuse together and create other elements called the fusion products. For example, two atoms of hydrogen fuse into an atom of helium. As a result, the hydrogen atoms loose mass while fusing. As explained by Einstein’s famous relation(1.1), the mass is energy
E = ∆mc2, (1.1) where ∆m is the difference of mass between the reactants and the products.
In a fusion reactor, deuterium and tritium are the reactants. The fusion re-action of deuterium and tritium is showed Figure 1.2. It produces total energy output of 17.6 MeV, which can be divided in the energy of an alpha particle (3.5 MeV) and the energy of a neutron (14.1 MeV) [2].
1.2.1
The tokamak
The plasma reaches extreme temperature and it would melt any other materials in direct contact. Therefore, a specific property of the plasma is exploited
CHAPTER 1. INTRODUCTION 3
Figure 1.2: Scheme of the fusion between deuterium and tritium
to keep itself far from the wall. Since the plasma is an ensemble of charge particles, it can be contained by a magnetic field. In this thesis, the simulation will be focus on a device called tokamak, which use magnetic confinement. The shape of a tokamak is commonly referred as a torus. Figure 1.3 is a scheme of tokamak with a similar confinement system as ITER/JET.
Figure 1.3: Scheme of a tokamak
For a better plasma confinement, two magnetic fields are induced: a poloidal magnetic field and a toroidal magnetic field (Figure 1.3). The combination of
4 CHAPTER 1. INTRODUCTION
the two fields result in a helical magnetic field (Figure 1.3), which reduces the particle drift towards the wall.
The magnetic fields keep the plasma separated from the wall. The limit of the plasma is called the separatrix and the gap between the separatrix and the wall is called the scrape-off-layer (SOL). At the separatrix, there is large gra-dients in density and temperature since the core of the plasma is dense and hot whereas the region near the wall is cold and low in density. This transitional region is called a pedestal. The density profile showing the pedestal can be seen Figure 1.4.
Figure 1.4: Electron density profile at the separatrix
1.2.2
Heating processes
As discussed above, the plasma needs to be heated to extreme temperature. ITER will rely on a combination of three heating methods [3].
The first one is the ohmic heating, a current is driven through the plasma and coupled to the plasma resistivity it produces heating following Ohms law PΩ = ηj2, (1.2)
with η the plasma resistivity and j the current density. The main issue with this method is that the higher the temperature of the plasma is the lower the resistivity become. It means that over a certain temperature the ohmic heating is not efficient anymore.
CHAPTER 1. INTRODUCTION 5
The second method is the neutral beam injection (NBI), a beam of neu-tral particle is launched with a high velocity into the plasma. This beam will accelerate and give energy to the ions by colliding with them.
The last method is the radio-frequency heating (RFH), which can be split in three categories: ion cyclotron resonance heating (ICRH), electron cyclotron resonance heating (ECRH) and lower hybrid resonance heating (LHRH). In this master thesis, the focus will be on ICRH. An electromagnetic wave is launched into the plasma by an antenna. The antenna is located outside of the plasma, mounted on the wall. The ions gyrate around the magnetic field lines at a specific frequency, the ion cyclotron frequency. When an electromagnetic wave propagates into the plasma at the same frequency the wave is damped and transfers its energy to the ions. This is represented in figure 1.5.
Figure 1.5: Scheme of ICRH. The black line represent the wave at the fre-quency ω and the dot red line represent the gyration of an ion around one of the magnetic field line at the frequency ωc.
1.2.3
The scrape-off-layer
In Fusion device the frequency of the antenna is fixed at approximately 50 MHz, since this frequency matches the ion cyclotron frequency of deuterium at the centre of the plasma.
However, the density changes by at least an order of magnitude between the plasma centre and the SOL. Hence, the frequency of the antenna could also match another characteristic frequency of the plasma: the lower hybrid frequency (LHF) [4]. The LHF depends on the electron density. ITER will be bigger and the density might be lower in the SOL than it has been for JET. For example the electron density is expected to be lower than 3.7 × 1017m−3
[5]). Consequently, the lower hybrid resonance (LHR) might be located in the SOL and heat the plasma in this region. In this thesis, we will study the
6 CHAPTER 1. INTRODUCTION
conditions required for the LHR to happen in the SOL and resolve the effect of this resonance.
This thesis does not aim to study the lower hybrid resonance heating that take place at the core of the plasma for a frequency of 5 GHz [5] but to verify that the resonance does not stop the propagation of a wave of great amplitude towards the centre of the plasma. If the wave is trapped in between the reso-nance and the wall, the energy is lost but most importantly this energy could also damage the equipment on the wall (ex: the antenna).
In that regard, the FEMIC code [6] (Finite Element Model for Ion Cy-clotron heating) will be used. This code has been developed with the objec-tive to build a 2D axisymmetric high resolution model that can describe the wave physics in the core plasma and the regions outside. The code is using a combination of two software, the calculation of the dielectric tensor is done in MATLAB and the rendering is done in COMSOL.
FEMIC handles the wave equation in the core of the plasma and in the SOL. Nonetheless, the code is not showing yet the effect of the LHR in the SOL. The main problem is that the wave vector goes to infinite values at the resonance. The corresponding resolution would then be infinitely small and this is not realistic to use very small mesh size due to numerical requirements. Not being able to fulfil this condition would result in being enable to see the impact of the LHR while studying the power absorption for example. To solve that problem, one can introduce a small quantity in the imaginary plane. This can be done by introducing another quantity, the collision frequency.
To conclude, the goal with this master thesis is to: • Implement collisions into the FEMIC code.
• Implement the collisions into the dispersion relation, study the disper-sion relation for the slow wave in the SOL and estimate the resolution required to resolve the effect of the LHR.
• Model the propagation of the slow wave in the scrape off layer using FEMIC.
• Use the FEMIC code to study the impact of the lower hybrid resonance on the power coupling between the antenna and the plasma.
Chapter 2
Wave and plasma physics
Two different wave models are used in this thesis, the cold and the hot plasma model. The particles, in the cold plasma model, do not have any kinetic ther-mal motion of their own. Their motion is induced by the action of the elec-tromagnetic field of a wave, if there is none, particles are considered at rest. On the contrary, the hot plasma model considers the effects of temperature. Therefore, the particles are not at rest when there is no electromagnetic field. They are already moving and have a kinetic thermal motion on their own.
In the following sections, if not specified the units are by default SI units.
2.1
The wave equation
Most of the plasma waves can be described by the wave equation. The wave equation 2.1 is derived from Maxwell’s equations (see appendix), and given by
∇ × (∇ × E) −ω
2
c2K · E = iωµ0Jant, (2.1)
where ω is the frequency, c the speed of light, µ0the vacuum permeability, E
the electric field, K the dielectric tensor and Jantthe current provide by the
antenna (external current).
2.2
The dielectric tensor in a hot plasma
The Dielectric tensor describes the response of the plasma to an electric field within it. The hot plasma model can be used to describe the dielectric tensor [6]. The coordinate system is defined so that the z component is along the
8 CHAPTER 2. WAVE AND PLASMA PHYSICS
magnetic field, thus kz = kk. The coordinate axis y is chosen to be
perpendic-ular to the wave vector, such that kx = k⊥, and ky = 0. This is a customary
rotation of the coordinate system that makes the symmetries of the dielectric tensor more apparent.
The dielectric tensor can be expressed as
K = K1 K2 K4 −K2 K1+ K0 −K5 K4 K5 K3 , (2.2) with the tensor elements
CHAPTER 2. WAVE AND PLASMA PHYSICS 9 K0 =2 X j ωpj2 e−λj ωkzvkj ∞ X n=−∞ λj(In− In0) 1 − kzv0j ω Z (ζnj) +kzvkj ω 1 −T⊥j T1j Z0(ζ nj) 2 (2.3) K1 =1 + X j ω2 pje −λj ωkzv`j ∞ X n=−∞ n2I n λj 1 − kzv0j ω Z (ζnj) +k2v`j ω 1 −T⊥j T1j Z0(ζ nj) 2 (2.4) K2 =i X j jωpj2 e−λj ωkzv`j ∞ X n=−∞ n (In− In0) 1 −kzv0j ω Z (ζnj) +kzv`j ω 1 −T⊥j Tkj Z0(ζ nj) 2 (2.5) K3 =1 − X j ω2 pje −λj ωkzv`j ∞ X n=−∞ In ζnj + v0j v`j × 1 + nωcj ω 1 − Tkj T⊥j Z0(ζnj) +2nωcjTkjv0j ωT⊥jv`j " Z (ζnj) + ζnj+ v0j v`j −1#) (2.6) K4 = X j k⊥ωpj2 e−λj kzωωcj ∞ X n=−∞ nIn λj nωcjv0j ωv`j Z (ζnj) + T⊥j Tkj − nωcj ω 1 − T⊥j T1j Z0(ζ nj) 2 (2.7) K5 =i X j k⊥jωpj2 e−λj kzωωcj ∞ X n=−∞ (In− In0) nωcjv0j ωv`j Z (ζnj) + T⊥j Tkj − nωcj ω 1 − T⊥j Tkj Z0(ζ nj) 2 , (2.8) where ζnj = ω0 + nωcj− kzv0j kzv`j , (2.9) Z(ζ) ≡ √1 π Z ∞ −∞ e−ξ2 ξ − ζdξ. (2.10) Here, Tkj and T⊥j are the temperature components along the parallel and
10 CHAPTER 2. WAVE AND PLASMA PHYSICS
species j and vlj the thermal velocity in the longitudinal direction. Inis the
modified Bessel function of the first kind which use λj as an argument,
λj = 1 2k 2 ⊥ρ 2 L, (2.11)
with ρLthe Larmor radius.
The dielectric tensor components are expressed in a slightly different form that in Swanson’s book [7]. The frequency expressed inside the term ζnas ω0
is derived from the momentum equation. On the contrary, ω comes from the derivative in time of some other terms. In the derivation performed here, ω and ω0are identical but in the chapter 2.4 a collision term will be added to ω0.
2.3
Dispersion relation
The dispersion relates the wave vector and the frequency of the wave. It shows the effect of the dispersion on the properties of the wave in a medium. The dispersion relation is obtained from the wave equation (2.1). To simplify the wave equation, the refractive index can be used
n = kc
ω . (2.12)
Assuming plane wave and no external current, the solution of the wave equation (2.1) can be found using
n × (n × E) + K · E = 0, (2.13) using the wave operator
Λ · E = 0, (2.14) There are two solutions to the wave equation, the trivial solution where the electric field is zero and the non-trivial solutions where the determinant of the matrix Λ is zero. det(Λ) = det Kxx− n2z− n2y Kxy + nxny Kxz + nxnz Kyx+ nynx Kyy − nz2− n2x Kyz+ nynz Kzx+ nznx Kzy+ nzny Kzz − n2x− n2y = 0. (2.15) The coordinate system used in the previous section implied that ky = 0,kx =
CHAPTER 2. WAVE AND PLASMA PHYSICS 11 det(Λ) = det K1− n2k K2 K4+ nkn⊥ −K2 K1+ K0 − n2 −K5 K4+ nkn⊥ K5 K3− n2⊥ = 0. (2.16)
In a typical fusion plasma with modest temperature effects, an ordering of the dielectric tensor elements can be found when ω ∼ ωc. The ordering may
be express as K3 >> K1, K2 >> K5, K4, K0.
Then K3can be evaluated in comparison to the refractive index to identify
the slow wave from the fast wave. Firstly, if K3 is large compared to n2⊥, it
means that the determinant of the wave operator can be simplified as
det(Λ) = K1− n2k K2 −K2 K1+ K0− n2 = 0. (2.17) From this equation and the relation 2.12 the perpendicular wave vector can be found as k⊥2 = ω 2 c2 K1+ K0− n 2 k+ K2 2 K1− n2k ! . (2.18) This simplification leads to the Fast magnetosonic Wave (FW) solution. The parallel wave vector in a tokamak geometry can be approximated by
kk =
nφ
R, (2.19)
with nΦ the toroidal mode number and R the major radius of the tokamak.
In the second approximation, K3 is considered of the same order of n2⊥.
After some algebra, the determinant of the wave vector can be written as
det(Λ) = K1− n2k nkn⊥ nkn⊥ K3− n2⊥ = 0. (2.20) It leads to the Slow magnetosonic Wave (SW) solution
k⊥2 = ω 2 c2K3 1 − n2k K1 ! . (2.21) As long as k2
⊥ > 0, the wave is propagating. But if k⊥2 is negative, it means that
the value of k⊥is imaginary. The wave is then an evanescent wave and decays
with a decay length Im(k⊥). The transition between a evanescent wave and
a propagating wave is called a cut-off. It occurs when the index of refraction squared is equal to zero (when k2
12 CHAPTER 2. WAVE AND PLASMA PHYSICS
2.3.1
Lower Hybrid resonance frequency
The lower hybrid frequency (LHF) is characteristic frequency of the plasma and leads to a resonance when matching the frequency of the wave. The res-onance is a singularity, or a pole in the function
k⊥(ω) =
1 ω − ωLH
, (2.22)
where ωLH the LHF and ω the frequency of the wave.
It is expected that,during ICRH simulation, the LHR will be present in the SOL for the SW [4]. Since the temperatures in the SOL are low enough, the particles are assumed to have no kinetic thermal motion of their own. The cold plasma model can be used to derive the LHF.
The LHR happens for an electric field parallel to the magnetic field and a wave vector perpendicular to the magnetic field. In Swanson [7], the LHF is given by ωLH2 = ΩeΩi ω2 pe+ ΩeΩi ω2 pe+ Ω2e , (2.23) where Ωαis the absolute value of the cyclotron frequency of the species α and
ωpethe electron plasma frequency.
2.4
The LHR in FEMIC
As it has been said previously, the FEMIC code handles the wave equation in the core of the plasma and in the SOL. However, the code do not model the lower hybrid resonance correctly, where the wave vector goes to infinity at resonance. The corresponding resolution would then be infinitely small, which is not achievable numerically. Not being able to fulfil this condition would limits the study on the impact of the LHR. The problem can be resolved by moving the pole into the complex plane. It is done by adding an imaginary term to the frequency ω.
For example, let’s take the resonance that happens at ω = ωLH.
1 ω − ωLH
. (2.24)
By replacing ω by ω + iν, we obtain 1 ω + iν − ωLH = 1 ω + iν − ωLH ×ω − iν − ωLH ω − iν − ωLH = ω − iν − ωLH (ω − ωLH) 2 + ν2. (2.25)
CHAPTER 2. WAVE AND PLASMA PHYSICS 13
The denominator in equation 2.25 does not go to zero anymore, so the sin-gularity disappears and the study of the LHR is then possible. By introducing collisions, an imaginary term proportional to the collision frequency appears in the momentum equation . This imaginary term will move the pole in the complex plane. The implementation of the collisions in FEMIC, is explained in the chapter 3.1.
Chapter 3
Collisions
In this chapter, the theory of the collisions will be described, as well as its implementation in the dielectric tensor. Collision is the mechanism chosen to move the pole of the LHR to the imaginary plane.
3.1
Theory
In this model, the collision frequency for one species is considered as the inter-action of one particle with a cloud of particles (from all species). The Coulomb electric field forces, produced by individual background particles, are small. But, the effect of successive long-range elastic Coulomb collisions on to a sin-gle particle, with a background of charged particles, leads to the modification of the single particle’s motion.
In this thesis, the formalism of the book "Collisional Transport in Mag-netized Plasma" by Per Helander and Dieter J. Sigmar [8] has been chosen to express the collision frequency.
The total collision frequency of a single species a is
νa= X b nbe2ae2bln Λab 3π3/22 0m2avT a3 , (3.1)
with nb the density of species b, ea the charge of a, eb the charge of b, ma
the mass of a, vT a the thermal velocity of specie a and ln Λab the Coulomb
logarithm for species a and b. The Coulomb logarithm is described chapter 3.1.1.
CHAPTER 3. COLLISIONS 15
3.1.1
Calculation of the Coulomb logarithm
The NRL formulary 2019 [9] is used to calculate the Coulomb logarithm. The density nα is expressed in CGS units, the temperature Tα in eV and the mass
mαin kg.
For electron-electron collisions the Coulomb logarithm is: ln Λee= 23.5 − ln n1/2e T −5/4 e − 10 −5 + (ln Te− 2)2/16 1/2 . (3.2) For electron-ion and ion-electron collisions the Coulomb logarithm is
ln Λei = λie = 23 − lnn1/2e ZTe−3/2 , Timmei < Te < 10Z2eV; 24 − lnn1/2e Te−1 , Timme i < 10Z 2eV < T e; 16 − lnn1/2i Ti−3/2Z2µ, Te < Timmei, (3.3) with µ = mi
mp the ratio between the ion mass and the proton mass. For ion-ion
collisions the Coulomb logarithm is
ln Λii0 = ln Λi0i = 23 − ln ZZ0(µ + µ0) µTi0 + µ0Ti niZ2 Ti +ni0Z 02 Ti0 !1/2 . (3.4)
3.2
Implementation in the Dielectric Tensor
Here, the implementation of the collision frequency in the dielectric tensor will be shown. For simplification, the cold plasma model is used. In the first place, let’s assumed that there are no collisions. Then, the momentum equation is written as
m∂v
∂t = q (E + v × B) . (3.5) Equation (3.5) can be simplified by assuming plane waves and introducing a perturbation.
E = E1ei(k.r−ωt) (3.6)
v = v1ei(−ωt) (3.7)
B = B0+ B1ei(k.r−ωt). (3.8)
Here, the background magnetic field B0 is assumed to be along z and much
16 CHAPTER 3. COLLISIONS
electric field and velocities. The equation of motion becomes −iωδij − ij3 qB0 m vj = qEi m . (3.9) From equation (3.9), the conductivity tensor and therefore the dielectric tensor can be obtained as described in (2.8).
Now, the equation of motion is written considering the collisions
m∂v
∂t = q (E + v × B) − νmv, (3.10) where ν is the collision frequency. Using the same assumptions and notation as in equation (3.9), results in (−iω + ν)δij − ij3 qB0 m vj = qEi m . (3.11) In this equation we can do a change of variable and consider the collision frequency as a part of an effective wave frequency
ω0 = ω + iν. (3.12) If the ω0 of the equation (3.12) is implemented in the equation (3.5) it gives
the equation of motion with collisions (3.11).
The same substitution can be done during the derivation of the warm plasma dielectric tensor. The collisions frequency only influence one of the term of the different dielectric tensor components (2.8). This term is the argument of the plasma dispersion relation function ζn. It is rewritten in all components as
ζn =
ω + iνj+ nωcj− kzv0
kzv`
Chapter 4
Modelling in FEMIC
The FEMIC code [6] will be used to solve the wave equation. This code has been developed with the objective to build a 2D axisymmetric high resolution model that can describe the wave physics in the core plasma and the SOL. The code is using a combination of two software, the calculation of the dielectric tensor is done in MATLAB and the calculation of the electric field is done in COMSOL using the finite element method.
This section will focus on the modification of the COMSOL configura-tion, with the aim to study the SW and LHR after the implementation of the collisions.
4.1
Meshing
The mesh of the ITER model is designed accordingly to the need of this the-sis. Since the mesh size required to resolve the SW is (according to chapter 2.4) very small, it has been decided to build a customised mesh to reduce the computation time and improve efficiency. Free triangular meshes are used in the plasma core and SOL. Boundary mesh is applied at the separatrix and at the antenna straps. Boundary meshes in this thesis are define with the same size as the free triangular one. Two other domains are created near one of the antenna straps in the SOL and the other around one antenna strap.
In the SOL and the core plasma the mesh size is set at approximately 0.06 m. It is defined by the minimum wavelength required to resolve the FW inside the plasma. Near the antenna, the element size is 0.02 m. A smaller domain around one of the antenna strap is created with the aim to show the LHR and the propagation of the SW. The mesh size of this domain will be define ac-cordingly to the wavelength of the wave at this location (i.e. 5 grid elements
18 CHAPTER 4. MODELLING IN FEMIC
per wavelength). This will be explained in more details chapter5.2. The corresponding mesh is given figure 4.1.
Figure 4.1: Description of the meshing in ITER. The vertical axis correspond to the height z of the tokamak and the horizontal axis corresponds to the Major radius R. On the left side, the meshing in core plasma is shown with a sparse meshing in the SOL. On the right side, a zoom on the domain of study near the antenna is done.
4.2
Tilting of the magnetic field
In reality, the magnetic field and the antenna are not perfectly perpendicular to each other. The magnetic field is slightly tilted and have components in both the toroidal and poloidal directions, being respectively perpendicular and parallel to the antenna. The tilting can be represented as a small angle α. This angle is related to the ratio between the toroidal and poloidal magnetic fields which is approximately 1/10. Then, α = arctanBθ Bφ = arctan 1 10 = 5.71°
It is not possible to do a simple modification of the model to tilt the mag-netic field. However, the geometry of the antenna can be changed. The current density inside the antenna is tilted by an angle α = 5.71°. It leads to a more accurate interaction between the magnetic field and the antenna.
CHAPTER 4. MODELLING IN FEMIC 19
Figure 4.2: Scheme of the tilting between the antenna and the magnetic field. The angle α is shown, being the angle between the poloidal and toroidal arc lengths with r the minor radius and R the major radius
Chapter 5
Results
In this chapter, we first present the implementation of the collisions using the JET model. Then, the influence of the collision on the dispersion relation will be analysed and the required resolution will be determined. This will be replicated for the ITER model. It follows a study of the propagation of the SW toward the resonance using the finite element method.
5.1
JET
The LHR is, in reality, unlikely to appear for JET model. We will built a theoretical and ideal scenario to observe the collision frequency and its impact over the dispersion relation. The purpose of using JET, before investigating for ITER, is to prepare the study of the dispersion relation and slow wave in ITER with an easy analysis of the LHR without any kinetic mode effect interfering. The temperature and density profiles are set in the SOL such that the LHR appears at the major radius 3.9 m. The pedestal in temperature and density profiles have been observed to have different width and to appear at slightly different locations. Such difference can cause rapid variation in the collisions frequency since ν ∼ n
T3/2. We want to avoid those rapid variations because it
results in high spikes in the collision frequency near the pedestal. The study of the impact of the collision frequency on the dispersion relation would then be difficult in this region. So, it has been chosen to use the same pedestals for the temperature and the density. The density and temperature profiles for the electrons are presented figure 5.1a and 5.1b. The temperature is isotropic. The ion species follow the same profiles except that the maximum temperature for ions is usually assumed to be smaller than for the electron (in this case 3200 eV).
CHAPTER 5. RESULTS 21
(a) (b)
Figure 5.1: Electron density profile (a) and temperature profile (b) at the equa-torial plane for an isotropic temperature.
Four different species are introduced in the modelling below: beryllium, deuterium, hydrogen and electron. In this scenario, the collisions are inves-tigated for positively and negatively charged particles but also for varying masses and charge. The beryllium will be tested for a density equal to zero to check the validity of the model (i.e. the collision frequency of a species that has zero density does not impact the total collision frequency).
The antenna is almost at the equatorial plane. The equatorial plane is de-fined by the vertical position of the magnetic axis, at z = 0.34 m in JET and z = 0.509m in ITER. The study of the dispersion relation will be done at this location, since this is where the electric field is the strongest.
5.1.1
Collisions
The collision frequency in the JET tokamak has properties that may appear counter-intuitive; the collision frequency is higher by a factor 100 in the SOL than at the core of the plasma as shown figure 5.2. The density increases towards the centre of the plasma. Therefore, one would expect the curve for the collision frequency to have the same shape. However, the temperature is also increasing and the collision frequency in proportional to the inverse of the temperature at the power of 1.5 (equation (3.1)). As expected, the collision frequency of the ions is lower than the one of the electrons due to the difference of mass.
The maximum collisions frequency is located in the SOL, whereas the frequency in the centre of the plasma remains low. The sum of the maximum
22 CHAPTER 5. RESULTS
Figure 5.2: Collision frequency in the equatorial plane for the electron on the left (νe) and for the ion Deuterium on the right (νD).
Figure 5.3: Coulomb logarithm for Deuterium-electron collisions.
collision frequencies of each species in the SOL is ν ≈ 2.7 × 106 Hz.
The Coulomb logarithm is most of the time assumed constant in the liter-ature [10] and could, in this situation, be approximated by 17. However, the figure 5.3 shows that the Coulomb logarithm in the SOL and the core differ.
5.1.2
The dispersion relation
The resolution in FEMIC is limited and it causes problems when trying to resolve short wavelength near the LHR. Thus, waves with high magnitude of k⊥cannot be modelled.
By introducing the collision frequency, the pole at the resonance is moved into the imaginary plane and a continuous curve is obtained for k⊥. Reducing
the value of |k⊥|would allow a coarser mesh since the required mesh size is
CHAPTER 5. RESULTS 23
ignore the parallel component of the wave vector kk. The aim is to see how the
maximum wave number can be limited by adding collisions, such that the LHR can be resolved. Hence, the dispersion relation is first analysed in a program called DISPREL. This program focuses on the calculation of the dispersion relation only.
For now, it has been chosen to simplify the simulation and look at the dispersion relation for only two species: the hydrogen ions and the electrons. The code is calculating all the solutions of the dispersion relation using the dielectric components as shown section 2.2.
To identify the effect the LHR, it is required to determine the location of the phenomena. The LHR happens for low density and it is possible to estimate the corresponding electron density using equation (C.6) in appendix C. It is found for ne ≈ 4.6 × 1016m−3. According to the density profile (figure 5.1a),
the LHR should appears for a major radius R = 3.9 m
No collisions
Figure 5.4 shows the solutions to the dispersion relation in the SOL (R > 3.85 m). Here, the collisions are neglected. Here, k⊥is diverging between 3.869 m
and 3.870 m. This is the discontinuity expected from the chapter 2.3.
Figure 5.4: Dispersion relation squared plotted against the major radius R when there is no collisions
Natural Collisions
The collision frequency ω0 = ω + iν is implemented in the dielectric tensor
seen in section 2.3. The collision frequency used is the one calculated numer-ically in the section 5.1.1. The dispersion relation obtained is similar to figure
24 CHAPTER 5. RESULTS
5.4. The curve is not going down and the discontinuity remains. The curve might be going down at a higher value of k⊥ which does not appear on the
figure. Since the domain over which we are operating is between -5000 and 5000 m−1, the value of k
⊥would then be greater than 5000 m−1which means
that the resolution required would be
λ = 2π k⊥
= 2π
5000 ≈ 0.0012m. (5.1) The finite element method will be used in the next step and in order to get reasonable results, there should be 5 grid element per wavelength. The corre-sponding mesh size is equal to 0.00024 m. Therefore, there is no need to check for values of k⊥ > 5000m−1since the corresponding mesh-size is too low to
get a reasonable computation time. It is not possible to resolve the LHR with this resolution.
The collision frequency in the SOL for JET is of the order of 106Hz, which
is not enough. The next step was to try an artificial collision frequency higher than this value.
Artificial collisions
Since the "natural" collision frequency of the plasma is too small, an artificial collision frequency can be introduced to increase the minimum wavelength that needs to be resolved in FEMIC. Using a frequency equal to or higher than 6.5 × 107Hz allows to see a relevant effect of the artificial collision frequency over the dispersion relation. An example is given in figure 5.5.
Figure 5.5: Dispersion relation squared plotted against the major radius R for an artificial collision frequency νart = 6.5 × 107Hz
CHAPTER 5. RESULTS 25
The maximum value at the resonance is lower than the value correspond-ing to the resolution in FEMIC. The maximum, k⊥at the resonance is
approx-imately 300 m−1. Using the relation (5.1), the obtained wavelength is:
λ = 2π
300 ≈ 0.02m, (5.2) with 5 grid element per wavelength, the mesh size is 0.004 m. It is a reasonable mesh size, hence the LHR can be modelled in JET.
The dispersion relation gives a backwards propagation of the SW. When there is no collisions, the propagation seems to be forward, but it is due to the imaginary part of k⊥being 0 at this location. DISPREL cannot properly
eval-uate the direction of propagation and set it to forward by default. By assuming cold plasma model and a propagation of the wave from the antenna towards the plasma at the location of the SW, the calculation of the group velocity gives a negative value. So, the wave is indeed supposed to propagate backwards.
5.2
ITER
The ITER tokamak will be bigger and have different parameters values than in JET. The temperature at the core of the plasma can reach few tens of keV and the density will be around 1020m−3 [11]. In the SOL, the density can be
estimated to be of the order of 1017-1018m−3. Moreover, the magnetic field is
on average stronger than in JET. It implies that the LHR electron density needs to be calculated again. For hydrogen only, the formula does not apply since the frequency of the wave remains lower than the LHR for any electron density. So, the LHR cannot be model in ITER for a hydrogen plasma. However, by introducing a deuterium minority, the cyclotron ion frequency changes (mass two times heavier) and the LHR electron density is found to be ne≈ 8.6×1016
m−3. A typical fusion plasma also uses a third species. The chosen third
species is Helium 3 and have density equal to 3% of the total density of the plasma.
The "natural" collisions frequency in the plasma is essentially the same as in JET (8.5 × 105) . But it is strongly affected by the profile in temperature
and density. Since those profile are not well known and that in section 5.1.2 this frequency has been shown to have a really small impact on the dispersion relation, it is better to leave this collision aside for the moment. Therefore, the following plots shows a LHR for an artificial collision frequency.
26 CHAPTER 5. RESULTS
5.2.1
The dispersion relation
In this section, the aim is to show that the LHR also appears for ITER using DISPREL, and to find an appropriate artificial collision to resolve the propa-gation of the SW to the resonance and its behaviour at this location.
The first step has been to plot the dispersion relation for the cold plasma model using DISPREL. Indeed, for low temperature the kinetic effects are not very important and allows to see the LHR clearly. This first step allows us to check that the lower hybrid resonance exists for this model.
In figure 5.6 the dispersion relation for the ITER edge and SOL region is shown.
Figure 5.6: Dispersion relation squared against the major radius R for a cold plasma and for an artificial collision frequency.
Figure 5.6 shows the lower hybrid resonance at R ≈ 8.27 m for an artificial collision frequency of νart = 6×106Hz. As one might think, this artificial
col-lision frequency implies that a relatively high resolution in FEMIC is required since the k⊥≈ 1000m−1. The corresponding wavelength is λ = 10002π = 0.006
m. Since there should be 5 grid element per wavelength, the grid size should be approximately 1 mm. It can be modelled in FEMIC by refining the mesh and consequently increasing the computation time. The main reason to refine the mesh instead of increasing the artificial collision frequency (like in the JET model), is the already high value of k⊥away from the resonance. By
increas-ing the artificial frequency, the wavelength away from the resonance will be affected. Moreover, a damping is induced by the introduction of the artificial collision frequency.
The cold plasma approximation leads to results that match our expecta-tions. However, the dispersion relation does not behave the same for hot plasma.
CHAPTER 5. RESULTS 27
Indeed, what has been assumed to be high kinetic modes appears in hot plasma. It is plotted on the green curve for values of k2
⊥ above 107 m−2 between 8.1
m and 8.2 m figure 5.7. This was not the scope of this thesis to investigate those modes. Those modes should not affect the study and are for that rea-son kept from appearing in the SOL. This is done by maintaining a really low temperature in the SOL (few eV).
Figure 5.7: Dispersion relation squared against the major radius R for a hot plasma and for an artificial collision frequency of 5 × 106Hz
In figure 5.7, the LHR can clearly be seen between 8.27 m and 8.3 m and has the same behaviour as in JET. The mesh size corresponding to the maxi-mum wavelength is again approximately 0.001 m.
According to the results section 5.2.1, the artificial collision is enough to resolve the resonance. However, even if the real part of k⊥remains low enough
when using the high artificial collision, its imaginary part becomes too high and the wave is damped before the wave reaches the resonance. A reasonable damping length would be few centimetres, which corresponds to =(k⊥) ≈ 15
m−1. The artificial collision frequency needs to be reduced to ν = 3 × 105
Hz to produce the desired damping. This is shown figure 5.8, the value of k imaginary remains above 15 m−1 between the antenna and the resonance. It
gives an e-folding length around 6 cm.
It is currently impossible to numerically model the behaviour of the SW at the LHR since it requires a mesh size below 1 mm. We did not have at the lab enough computational power to run such simulation in reasonable amount of time. Nonetheless, the propagation of the SW away from the LHR can be modelled since the corresponding mesh size remains slightly higher than 1 mm for ν = 3 × 105 Hz.
28 CHAPTER 5. RESULTS
Figure 5.8: Dispersion relation for an artificial collision of 3 × 105 Hz (hot
plasma model). On the left side , a zoom on the real part <(k⊥)shows the
wavelength of the propagation of the SW. On the right side, the zoom is done on the imaginary part =(k⊥)which describe the damping of the SW.
5.2.2
Modelling of the slow wave in FEMIC
In this section, we aim to show the propagation of the SW in the SOL towards the resonance and quantify the fraction of energy that belongs to the SW. The grid from section 4.1 is used. The mesh size of the domain of study (with the finer grid) is known from section 5.2 to be approximately 1 mm. The real value of the parallel component of the electric field is given in figure 5.9a.
When there is no collisions, see figure (5.9b), the wave is not damped any-more due to the introduction of a collision frequency. The discretization do not allow to model the propagation of the wave towards the other domain (SOL and plasma). Indeed, the mesh size used outside of the domain of study do not allow to model the propagation of the SW. Similarly at the resonance the wavelength of the SW decreases, which means that we do not maintain the requirement of 5 grid per wavelength anymore. And since the LHR is not resolved, the wave is not damped at this location but reflected. The wave is propagating between the antenna and the resonance but is reflected every time it reaches a domain with a larger mesh size. The wave should propagate every-where in the SOL and reach the resonance, every-where it gets absorbed. Due to the many reflections, the electric fields is much larger and it creates an eigenmode pattern.
The reflection at the LHR seen in figure 5.9b is the reason why it is im-portant to include artificial collisions. In fact, damping the wave before the resonance limits the reflections and allows a qualitative study of the SW. Fig-ure 5.9a shows the propagation of a wave from the antenna towards the plasma.
CHAPTER 5. RESULTS 29
(a) (b)
Figure 5.9: Parallel component of the Electric field for a tilted angle α = 5.71°. For an artificial collision frequency 3 × 105 Hz (a). Without collisions (b)
Figure 5.10: 1D plot of Ek against the major radius at the cut line (red line
30 CHAPTER 5. RESULTS
In order to study of the damping of the SW, a 1D plot of the real part of the parallel electric field is shown in figure 5.10). The electric field is plotted along the cut line shown in red in figure 5.9a. The damping of the wave and the wavelength can be measured on this figure. The wavelength is approximately 1 mm and the damping length is approximately 5 cm. These results are similar to the value obtained earlier using DISPREL in section 5.2.
The Faraday screening should reflect the slow wave back to the antenna and reduce the amount of power that belongs to the SW to reach the plasma. It effectively reduces the coupling to the SW. Behind the white squares in figure 5.9a, the wave decreases. In between the squares, the wave remains strong. It can be assumed that the wave is only weakly affected by the Faraday screen at this location.
The efficiency of the artificial damping is also proven since the wave is almost completely damped before the critical location of 8.29 m. This critical location is where the wave vector start to increase drastically near the LHR, the resolution does not match the wavelength of the wave. This is shown figure 5.11.
There is the possibility to test the damping without the Faraday screen at all since it is part of the geometry of the ITER model in FEMIC. However, it would need more work than what we wanted to put into it.
Figure 5.11: Dispersion relation against the major radius for an artificial col-lision of 3 × 105 Hz. The picture on the left is obtained using DISPREL and
the one on the right is obtained in FEMIC (Two notation are used R and r but both corresponds to the major radius).
Figure 5.11 shows a comparison of the wave vector predicted by DISPREL and an estimate of the wave vector from the FEMIC results. The latter is ob-tained from the assumption Ez = E0exp(R ik(r)dr), thus =[k(r)] can be
estimated as
k(r) = −i d
CHAPTER 5. RESULTS 31
In FEMIC, the wave stops its propagation around 8.29 m. The wave vector at this location reaches 4500 m−1, which is 4.5 times higher than the k
⊥used
to define the number of grid element per wavelength. It means that instead of having 5 grid elements per wavelength there is only 1. When there is less than 5 grid points per wavelength the estimate becomes unreliable. When the number of grid point per wavelength becomes less than 2 the wave cannot propagate. Between 1 and 2 grid point per wavelength we may have a combination of damping and reflection, while the wave field patterns tend to be irregular and oscillate from grid-cell to grid-cell.
5.2.3
Tilting of the antenna
In section 4.2, we described the tilting of the antenna and the reason why it is needed. We want to compare the solution between the non-tilted antenna and the tilted antenna. The Poynting flux is a good criteria to evaluate the sensitivity of the power inside the plasma to the angle. Indeed, the Poynting flux is proportional to the energy of the wave. The value obtained for a tilted angle α = 5.71° gives a flux 10 times higher than for an angle α = 0°. Thus, the tilting of the antenna is critical to include in any future work on the power coupling to the slow wave.
Chapter 6
Discussion
In this thesis, the propagation of the SW and the presence of the LHR has been demonstrated in ITER. However, there are limitations from the results obtained in the previous chapter that need to be discussed.
In the simulations presented above, the SW is reflected before reaching the lower hybrid resonance due to the insufficient resolution. An artificial damping (e-folding length around 6 cm) is necessary to limit the effects of the reflection and to show the propagation of the SW towards the resonance. However, there are consequences to the introduction of artificial damping. But the artificial damping does not only affect the SW but also of the FW. So, in a simulation involving the entire tokamak, having a strong artificial collision would impact the results at the centre of the plasma as well. The coupling between the SW and the plasma could also be affected.
In future simulations, it would be helpful to have more resources to use a smaller grid size. The smaller grid would allow us to resolve the slow wave near the LHR in the SOL. To be able to see the behaviour of the SW at this location, the propagation of the wave needs to be modelled from the antenna toward the resonance and at the LHR. Since the wave vector can reach values up to 5000 m−1, the grid size should be 5 times smaller than the current one.
We are working with a two dimensional geometry, so it would correspond to a number of degree of freedom 25 times higher. When resolving for the propagation of the SW before the resonance, we almost reached the 256 GB of RAM available. To run the simulation with a the finer grid in the same amount of time, it requires 6.4 TB of RAM. Improvements can be done. For example, the resolution could be locally enhance by creating a new sub-domain at the resonance. The propagation of the SW is modelled using the 1 mm resolution in the SOL and, when approaching the resonance, enters a domain of higher
CHAPTER 6. DISCUSSION 33
resolution (0.2 mm mesh size).
When the value of k⊥ increases, the mesh size needs to decrease
accord-ingly. in agreement with the previous results section 5.2.2, the wave vector is expected to increase drastically at R = 8.3 m. The maximum value of k⊥ at
the resonance will depend on the artificial collision. A compromise has to be done between limiting the maximum value of the wave vector and keeping the artificial damping as small as possible. Indeed, the SW is naturally damped when approaching the resonance but if the artificial damping is too strong the wave could be completely absorbed before reaching the LHR. The best options would be to use the natural collision frequency of the plasma, which is sensi-tive to the values of the density and temperature. Unfortunately, those values are currently not well-known.
The model could also be run for bigger domain but the size of the domain of study drastically increases the number of degrees of freedom. For example, including a full antenna straps would increase the degrees of freedom by a factor 4, and 4 straps by a factor 10.
Furthermore, the solver was using quadratic discretization of the electric field for numerical efficiency. This is a fairly good simplification when the antenna is tilted. However, tilting the magnetic field itself is much more accu-rate. Running FEMIC including a poloidal component of the magnetic fields requires a cubic representation of the electric field [12]. The cubic representa-tion of the electric field requires more memory and can therefore not be applied to study the slow wave without significant enhancement of our computational resources (by a factor 2 or 3).
The power coupled to the slow wave can be measured using the Poynting flux. Some preliminary results shows that the tilting of the antenna is important but more work is needed to quantitatively estimate the coupled power. Some test of the effects of the artificial collision on the coupling between the SW and the plasma and the absorption of the LHR could be done. It is expected that the coupling is not affected by the change in the artificial collisions.
Another issue is the presence of kinetic modes at high k⊥∼ 7 × 103m−1.
However, these modes are not modelled in FEMIC since the grid size at the core of the plasma is too small. To model them the grid size should be decrease below 0.001 m. This is currently impossible to expand such fine grid size to the whole plasma with the current computational power. Moreover, the coupling between the "high kinetic mode" wave and the plasma is expected to be relatively small. This is why it has been chosen not to study the kinetic modes in this thesis.
Chapter 7
Conclusions
During this thesis, the possibility of having a Lower Hybrid Resonance (LHR) in the Scrape-Off Layer (SOL) has been assessed in for specific values of tem-perature and density. The introduction of collisions in FEMIC shows a good agreement with the theory and the propagation of the Slow Wave (SW) at the LHR have been proven to be numerically solvable thanks to the introduction of an artificial collision frequency. The modelling of the SW in FEMIC at the location of the resonance has not been successful, since our present com-putational resources does not allow us to resolve the highest values of k⊥.
Nonetheless, the results of this thesis are encouraging since the propagation of the SW towards the LHR have been successfully modelled in ITER. With sufficient computation power, it is possible to expand the domain of study to the entire SOL and to resolve the SW at the LHR. There is a number of out-comes of this thesis. The first one is the improvement of the FEMIC code through the introduction of a new variable the collision frequency. The imple-mentation of collisions allows the study of the LHR. Secondly, it is important to know what will happen in front of the antenna. There are strong field near the antenna since the power emitted is around 20 × 106W [11]. If this power
is reflected, it could cause arcing and damage the antenna. If such damage is avoided, it might still cause problem for the transmission of the energy to the core of the plasma. Thus, extension to this study could be used to protect the antenna and improve heating efficiency.
Bibliography
[1] “BP Statistical Review of World Energy 2019, 68th edition.”
[2] J. Scheffel and P. Brunsell, Fusion Physics- Introduction to the Physics
Behind Fusion Energy, 6th Edition. Fusion Plasma Physics, Alfvén
Laboratory, KTH, 2016.
[3] “Iter external heating systems.” [Online]. Available: https://www.iter. org/mach/Heating
[4] M. Usoltceva, R. Ochoukov, W. Tierens, A. Kostic, K. Crombe, S. Heuraux, and J.-M. Noterdaeme, “Simulation of the ion cyclotron range of frequencies slow wave and the lower hybrid resonance in 3D in RAPLICASOL,” Plasma Physics and Controlled Fusion, vol. 61, no. 11, 2019.
[5] O. Meneghini and S. Shiraiwa, “Full wave simulation of lower hybrid waves in iter plasmas based on the finite element method,” Plasma and
Fusion Research, vol. 5, 01 2010.
[6] P. Vallejos, “Modeling RF waves in hot plasmas using the finite element method and wavelet decomposition: Theory and applications for ion cy-clotron resonance heating in toroidal plasmas,” Ph.D. dissertation, KTH Royal institute of Technology, 2019.
[7] D. G. Swanson, Plasma Waves, 2nd Edition. IoP Publishing, 2003. [8] P. Helander and D. J. Sigmar, Collisional Transport in Magnetized
Plasma. Cambridge University Press, 2005.
[9] A.S.Richardson, NRL Plasma Formulary. Pulsed Power Physics Branch,Plasma Physics Division, Naval Research Laboratory, Washing-ton,DC20375, 2019.
36 BIBLIOGRAPHY
[10] J. Wesson and D. Campbell, Tokamaks. Oxford: Clarendon Press, 2004. [11] A. M. et al, “Performance of the ITER ICRH system as expected from TOPICA and ANTITER II modelling,” Nuclear Fusion, vol. 50, 2010. [12] B. Ljundberg, private communication.
Appendix A
Wave equation
The wave equation is derived from the following vector identity
∇ × ∇ × E = ∇ (∇ · E) − ∇2E. (A.1) Using Maxwell-Faraday’s and Maxwell-Ampere’s equations, we get:
∇ × ∇ × E = ∇ × −∂B ∂t = −µ0 ∂J ∂t − 1 c2 ∂2E ∂t2 . (A.2)
The current J can be expressed as internal and external current. The external corresponds in our case to the current of the antenna.
J = Jint+ Jant. (A.3)
The internal current depends of the Electric field.
Jint= σ · E. (A.4)
In Fourier space and because of the wave Ansatz , the time derivative is sim-plified with a multiplication by −iω.
∇ × ∇ × E −ω
2
c2E − iωµ0σE = iωµ0Jant, (A.5)
K = I + i ω0
σ (A.6)
∇ × ∇ × E −ω
2
c2K · E = iωµ0Jant. (A.7)
Appendix B
Third component of the
dielec-tric tensor
First of all, we assume that Φ = 0, so k⊥ = kxand ky = 0. So, the dielectric
tensor can be written as
K = K1 K2 K4 −K2 K1+ K0 −K5 K4 K5 K3. (B.1)
In this section we are only interested in the component K3. Since the
dielectric tensor is derived from the expression of the current, this element will only depend on Jze~z.
J ≈ K · E. (B.2)
Jze~z = (K4Ex+ K5Ey + K3Ez) ~ez. (B.3)
Now, let’s assume a Maxwellian distribution so that:
f0(v⊥, vz) = F (vz) πv2 t e−v⊥2/vt2, (B.4) with v2
t = 2κT⊥/mthe transverse thermal speed.
In Swanson’s book [7], the following relation is obtained for the current.
J =X j njqj Z ∞ −∞ [hvxf1ji⊥e~x+ hvyf1ji⊥e~y + vzhf1ji⊥e~z]d3v, (B.5)
only the last term is along ~ez and corresponds to Jze~z.
APPENDIX B. THIRD COMPONENT OF THE DIELECTRIC TENSOR 39 hf1i⊥= q m Z ∞ 0 dτ eφ iv 2 t 2 (Axax+ Ayay) − αz , (B.6) where φ = i (ω − kzvz) τ − λ (1 − cos ωcτ ).
The dielectric tensor element K3 is defined by its multiplication with the
z component of the electric field. So,
(vzhf1i⊥)Ez = vzivt2qe −λ 2ωcm ∞ X n=−∞ (αxkx+ αyky) In λ + i(αykx− αxky)(In− I 0 n) −2αzωc In v2 t Ez ω + ωc− kzv0 . (B.7) Since we assume Φ = 0, then k⊥ = kx, ky = 0and αy = 0, we have
(vzhf1i⊥)Ez = vzivt2qe −λ 2ωcm ∞ X n=−∞ −k 2 x ω F0+ 2vz v2 t F In λ + 2F 0 ωc In v2 t ] Ez ω + ωc− kzv0 , (B.8) then σzzEz = Z ∞ −∞ (vzhf1i⊥)Ezd 3v. (B.9)
Using equation B.10 and B.11, the conductivity tensor can be rewritten equa-tion (B.12). F2(ζn) = v` kz ζn 2 + v0 v` Z0(ζn) − v2 0 v2 ` Z (ζn) . (B.10) Z ∞ −∞ dvz vzF0(vz) ω + nωc− kzvz = − 1 kzv` ζn+ v0 v` Z0(ζn) . (B.11) σzz = −iqe−λ 2ωcm ∞ X n=−∞ −v2 t kzvl ζn+ v0 vl Z0(ζn) + 2v2l kzvl ζn 2 + v0 vl Z0(ζn) − v20 v2 l Z(ζn) k2 xIn ωλ − 2 kzvl ζn+ v0 vl ωcInZ0(ζn). (B.12) Let’s consider all the component multiplied by Z0(ζ
n). Using the relation
λ = 12k2 ⊥ρ2⊥ = 12k 2 x v2 t ω2 c, we get:
40 APPENDIX B. THIRD COMPONENT OF THE DIELECTRIC TENSOR (σzz)Z0(ζ n)= −iqe−λ 2ωcm ∞ X n=−∞ − v 2 t kzvl ζn+ v0 vl 2ω2 cIn ωv2 t − 2 kzvl ζn+ v0 vl ωcIn +2v 2 l kzvl ζn 2 + v0 vl 2ω2 cIn ωv2 t Z0(ζn). (B.13) (σzz)Z0(ζ n)= −iqe−λ m ∞ X n=−∞ 1 kzvl In ζn+ v0 vl −ωc ω − 1 + ζn 2 + v0 vl 2v2 lωc ωv2 t Z0(ζn). (B.14) A simplification can be done on the last term:
2 ζn 2 + v0 vl = ζn+ v0 vl + v0 vl . (B.15) The transverse and longitudinal thermal velocity can written in function of the temperature. vl2 v2 t = Tk T⊥ . (B.16) (σzz)Z0(ζ n)= −iqe−λ m ∞ X n=−∞ 1 kzvl In ζn+ v0 vl 1 + ωc ω 1 − Tk T⊥ + v0Tkωc vlT⊥ω Z0(ζn). (B.17) The last term, v0Tkωc
vlT⊥ωZ
0(ζ
n)doesn’t appear in Swanson’s expression. The
reason is that Z0(ζ
n)is expressed as Z0(ζn) = −2[1 + ζnZ(ζn)].
Now we will look at the second part of the conductivity tensor,i.e. all the terms multiplied by Z(ζ) (σzz)Z(ζn) = −iqe−λ 2ωcm ∞ X n=−∞ −v2 0 v2 l 2k2 xωc2 ωk2 xvt2 InZ(ζn) 2v2 l kzvl . (B.18) (σzz)Z(ζn)= iqe−λ m ∞ X n=−∞ 2v2 0ωc kzvlωT⊥ Inz(ζn). (B.19)
Adding the term from (σzz)Z0(
zetan)rewritten as a function of Z(ζn).
(σzz)Z(ζn)= iqe−λ m ∞ X n=−∞ 2v2 0ωc ωT⊥ + 2v0Tkωc vlT⊥ω ζn In kzvl Z(ζn)+ 2 kzvl In v0Tkωc vlT⊥ω . (B.20)
APPENDIX B. THIRD COMPONENT OF THE DIELECTRIC TENSOR 41 (σzz)Z(ζn) = iqe−λ m ∞ X n=−∞ In kzvl 2v2 0ωc ωT⊥ + 2v0Tkωc vlT⊥ω ζn Z(ζn) + 2v0Tkωc vlT⊥ω . (B.21) Using ζn = ω+nωkzcv−kl zv0. (σzz)Z(ζn) = iqe−λ m ∞ X n=−∞ In kzvl 2 vlωT⊥ kzv2lv02ωc+ v0Tk(ω + nωc− kzv0 kzvl Z(ζn)+ 2v0Tkωc vlT⊥ω × kzvl ω + nωc ω + nωc kzvl . (B.22) (σzz)Z(ζn) = iqe−λ m ∞ X n=−∞ In kzvl ω + nωc kzvl 2v0Tkωc vlT⊥ω Z(ζn) + kzvl ω + nωc . (B.23) In equation (B.23) we can replace ω+nωc−kzv0
kzvl by ζn. (σzz)Z(ζn) = iqe−λ m ∞ X n=−∞ In kzvl ζn+ v0 vl " 2v0Tkωc vlT⊥ω Z(ζn) + ζn+ v0 vl −1!# . (B.24) The last step is to combine (σzz)Z(ζn)and (σzz)Z0(ζn).
K = 1 + i 0ω
σ. (B.25)
Using the relation (B.25 ) between the dielectric tensor and the conductivity tensor, we get: K3 =1 − X j ω2 pje−λj ωkzv`j ∞ X n=−∞ In ζnj + v0j v`j × 1 + nωcj ω 1 − Tkj T⊥j Z0(ζnj) +2nωcjTkjv0j ωT⊥jv`j " Z (ζnj) + ζnj+ v0j v`j −1#) . (B.26)
Appendix C
Lower hybrid resonance:
esti-mation of the density
The antenna is set at a certain frequency which has been first assumed to be 50 MHz. It implies that ωLH = 2π × 50 × 106 The LHR frequency can be
find using ωLH2 = ΩeΩi ω2 pe+ ΩeΩi ω2 pe+ Ω2e , (C.1) with Ωe = |q|B(R) me , Ω i = |q|B(R) mi and ω 2 pe = e2n e 0me. The magnetic B is
known in FEMIC through the relation
B(R) = R0B0
R (C.2)
The equation C.1 can be rewritten as ωLH2 × ω2 pe+ Ω2e ΩeΩi = ωpe2 + ΩeΩi . (C.3) ωpe2 1 − ω 2 LH ΩeΩi = Ωe ω2 LH Ωi − Ωi . (C.4) So, ωpe2 = Ωe ω2 LH Ωi − Ωi 1 − ω2LH ΩeΩi . (C.5)
It implies that ne is equal to
APPENDIX C. LOWER HYBRID RESONANCE: ESTIMATION OF THE DENSITY 43 ne= 0meΩe e2 ω2 LH Ωi − Ωi 1 − ω2LH ΩeΩi . (C.6) If a constant magnetic field B=1.4 T is assumed then the density is approxi-mately n ≈ 4.6 × 1016m−3.
Appendix D
Validity of the collisional model
To check the validity of the model, one method is to compare it to an other existing experimental, analytical or numerical solution. Here the model will be compare to the solution plotted in the book Wesson Tokamak page 730 [10]. The data used will be the collision time for electron-ion collision. This is done for a hydrogen plasma. Quasi-neutrality is assumed and only the values for an electron density of 1 × 1021will be used.
The way the collision time is calculated in the code makes it hard to ex-tract the exact collision time for electron-ion collisions. However, the theory can help us to simplify the process. The collision time of the electron comes from the inverse of the sum of collisions frequencies. In our case, the only ion species is the hydrogen, which means that the collision frequency for the electron, νe, is composed by νeeand νei. Since the ion is simplified by a
pro-ton, the charge Ze becomes equal to the electron charge. Therefore the quasi-neutrality leads to ne = ni. Over all the approximation that those two collision
frequencies are equal can be made. It is not exactly true whereas the Coulomb logarithm changes slightly but the approximation holds since the effects on the frequency are negligible compare to an approximation of the Coulomb loga-rithm to 17. Then the collision time of the electron can be divided by two and gives the collision time of the electron-ion collision.
Wesson formula correspond to the ion-electron collision time:
τe = 12π3/2 √ 2 ε2 0m 1/2 e Te3/2 niZ2e4ln Λ (D.1)
This equation can be found using (3.1) and the thermal velocity as (D.2):
vα = r 2Tα mα , (D.2) 44
APPENDIX D. VALIDITY OF THE COLLISIONAL MODEL 45
with b =ion and α =electron.
The Coulomb logarithm in Wesson’s book [10] is assumed to be constant in function of the temperature and approximately equal to 17, which has been seen earlier to be imprecise. An other information important to highlight is that the collision time calculated in Wesson is between an hydrogen ions and electrons. Hence, the simulation was run again for this specific case with a constant density set as 1021m−3 and tested for different temperatures.
Figure D.1: Electron-ion (only) collision time [10]
(a) Collision frequency (b) Coulomb logarithm for electron-ioncollisions
Figure D.2: Difference between FEMIC and Wesson
