Modeling the sawtooth oscillation in Tokamaks with a hybrid diffusion equation

(1)

Modeling the sawtooth oscillation in

Tokamaks with a hybrid diffusion equation

Jelmer de Boer

6079431

(2)

The bachelorthesis and -project, worth 12 EC, were supervised by: Menno Lauret (Eindhoven University of Technology)

(3)

1 Summary

Nuclear fusion involves the merging of lighter atomic nuclei into heavier ones. In the Tokamak reactor the plasma is free to move in the toroidal direction, giving rise to various instabilities, of which one is the ’sawtooth instability’ or the ’sawtooth oscillation’. The sawtooth instability is a phenomenon which leads into the mixing of the hot plasma core with the colder outer regions. This mixing can lead to the disruption of the fusion process, but, at the same time, can possibly be used to regulate the refueling of the core and, therefore, prevent the burning plasma core from being exhausted. Experiments show that the disruption does not so much depend on the amplitude of the insta-bility, but more on the period: longer periods are more likely to disrupt the fusion process. It is therefore desirable to learn how to control this period, so as to optimize between the benefits and drawbacks of this instability. There are various ways to manipulate this period. The first way to manipulate this period is to vary the diffusion constant, a. If doable, this might be a very efficient way to manipulate the period, because of the exponential influence the diffusion constant has on this period. The way this might be achieved is to vary the density and/or the temperature of the plasma, which, in turn, can be achieved by varying the applied magnetic field on the plasma in the reactor. Another way to manipulate the sawtooth period is to vary the value of scriticalor c1, which might also be done by manipulation of the density and

temperature of the plasma. The amplitude of the electric currents that run through the plasma is, based on this research, of no specific influence on the average sawtooth period. Because of the influence this amplitude has on the minimum sawtooth period, this amplitude can be used to enhance the refu-eling of the burning plasma core. The frequencies of the currents that run through the plasma can also be used for manipulation of the sawtooth period. These frequencies can be selected, so as to function as driving frequencies of the sawtooth period. This method might be the simplest for applications, because the variation of the temperature or density of the plasma could have different effects on different parameters. This method might, therefore, be the most applicable method that can be concluded from this research.

(4)

2 Samenvatting

Fossiele brandstoffen (aardgas, steenkool en aardolie) worden steeds schaarser en er wordt vermoed dat ze een negatieve invloed op het klimaat hebben. Hierdoor is er in toenemende mate behoefte aan alternatieve energiebronnen. Een veelbelovende alternatieve energiebron is kernfusie. Kernfusie behelst het samensmelten van twee lichte atoomkernen in zwaardere atoomkernen. Als het eindproduct van een samensmelting lichter is dan een ijzeratoom-kern, zal dit product een aanzienlijke hoeveelheid energie opleveren. Er is veel brandstof aanwezig voor kernfusie, wat kernfusie een duurzaam proces maakt. Ook wordt er relatief erg veel energie geproduceerd bij een fusiepro-ces, veel meer dan, bijvoorbeeld, bij het verbranden van olie. Deze twee redenen samen zorgen ervoor dat kernfusie een veelbelovende alternatieve energiebron is. De belangrijkste reden dat kernfusie vandaag de dag nog niet gebruikt wordt om energie mee op te wekken, is het feit dat het proces bij een extreem hoge temperatuur plaatsvindt. Deze extreem hoge tempe-ratuur zorgt ervoor dat de brandstof niet met de reactorwand in contact mag komen. Dit probleem kan omzeild worden door de brandstof met een heel sterk magneetveld uit de buurt van de reactorwand te houden. De be-langrijkste configuratie waarin dit gebeurt, heeft de vorm van een donut en wordt de Tokamak genoemd. Omdat bij de Tokamak een brandstof met een extreem hoge temperatuur in een heel erg sterk magneetveld beweegt, treden er verschillende interessante effecten op. Een van de processen die plaatsvindt, is de zogenaamde ’zaagtandoscillatie’. Deze periodieke oscillatie is het effect van een abrupte beweging in de brandstof, waardoor de warmere delen van de brandstof vrijwel instantaan met de koudere delen mengen. In de hetere delen van de brandstof vindt het fusieproces plaats en de zaagtan-doscillatie kan daarom het fusieproces verstoren. Voornamelijk de periode van deze oscillatie is van invloed op de verstoring van het fusieproces en het is daarom wenselijk deze periode te kunnen manipuleren. Dit onderzoek heeft als doel de invloed van verschillende parameters op de periode van de zaagtandoscillatie te onderzoeken.

(5)

3 Introduction

Fairly recently, plasma physics has, also quite literally, become a hot topic. This was caused by the fact that fossil fuels (gas, coal and oil) have become less abundant and are believed to have a negative influence on the climate. For these reasons, the demand for alternative energy sources has increased and a rapid development has taken place. One of the most promising, but also most challenging, alternative energy sources is nuclear fusion.

Nuclear fusion involves the merging of lighter atomic nuclei into heavier ones. If the heavy nucleus is lighter than Iron, the most stable nucleus, the sum of the weight of the lighter nuclei will be more than the weight of the heavy nucleus. Due to Einstein’s mass-energy equivalence, E = ∆mc2_{, this}

process will release a considerable amount of energy1.

Under normal, daily circumstances, nuclei will be surrounded by a cloud of electrons, which prevent them from fusing. The fusion fuels must be given such a temperature as to lose those electrons (in the order of millions Kelvin). The fuel is then said to be in the plasma state, the fourth state, next to solid, gas and liquid. After the fuels are brought into a plasma state, the temperature, pressure and density must be high enough to enhance the chance of two nuclei fusing.

Figure 1: The Tokamak configuration, also used in the ITER-project. The temperature of the plasma makes it impossible to put it into direct contact with a container, but the charged nuclei can be confined in a magnetic field. This fact is used in several reactor configurations, of which the Toka-mak, depicted in figure 1, is presently the most important (Toroidal’naya kameras magnitnymi katushkami, Russian: toroidal chamber with magnetic coils). This type of reactor will be used in the ITER-project (formerly: International Thermonuclear Experimental Reactor; iter, Latin: journey),

1_{For a slightly more detailed introduction on the subject, see the introductions of [1]}

(7)

an international cooperation building the largest fusion test reactor of the world in the Cadarache, France. This type of fusion reactor is therefore likely the first to be used for private energy consumption.

With the increasing demand for alternative energy sources, fusion research has become more important. This is due to the fact that nuclear fusion has the potential of delivering enough energy to sustain (or even increase) today’s standard of energy consumption. There are few to none other alternative energy sources that share this potential with nuclear fusion. What is also important: the fuels used for near-future fusion reactors are overwhelmingly abundant. In ITER, Deuterium and Tritium are used, both heavy isotopes of Hydrogen, which are obtained from, respectively, seawater and Lithium. There is enough Lithium available (in ore deposits, but mainly in seawater) to sustain today’s standard of energy consumption for over a million years. It is not (yet) possible to fuse Deuterium with Deuterium, for that requires a higher temperature in the reactor, but once that is made possible, there will be enough Deuterium available to increase today’s energy consumption by at least a factor 1000 for the coming million years2_{. So it is fair to say,}

that researching the possibility of nuclear fusion is indeed worth the effort.

(8)

4 Introduction into the topic of the thesis

In the Tokamak the plasma is free to move in the toroidal direction. There-fore charge will be moving, giving rise to a magnetic field of its own. This magnetic field leads to various instabilities, of which one is the ’sawtooth instability’ or the ’sawtooth oscillation’. The sawtooth instability3 _{is a}

phe-nomenon which leads into the mixing of the hot plasma core with the colder outer regions. Typical sawtooth behavior consist of an increasing or de-creasing temperature profile (inde-creasing for the core, dede-creasing for the outer regions), followed by a sudden ’crash’, in which the temperature suddenly rises or drops due to the mixing (drops for the core, rises for the outer re-gions). This mixing can lead to the disruption of the fusion process, but, at the same time, can possibly be used to regulate the refueling of the core and, therefore, prevent the burning plasma core from being exhausted. Experi-ments show that the disruption does not so much depend on the amplitude of the instability, but more on the period: longer periods are more likely to disrupt the fusion process. It is therefore desirable to learn how to control this period, so as to optimize between the benefits and drawbacks of this instability.

Figure 2: Schematic visualization of typical sawtooth behavior. The tempera-ture of the hot plasma core gradually rises, but then drops due to the mixing. The reverse process takes place in the colder outer regions

The mechanism behind the sawtooth instability is not yet fully under-stood. Several models have been proposed to simulate the behavior of the

(9)

plasma during and between crashes. This thesis is based on a model by Porcelli4_{, which, in large aspect-ratio Tokamaks, can be further simplified.}

The magnetic field lines in a Tokamak can be separated in a θ- and a φ-component. For large aspect-ratio Tokamaks, the toroidal Bφcan be assumed

constant. The safety factor, q, and the magnetic shear, s, defined as5_:

q(r, t) = rBφ RBθ(r, t) and s(r, t) = r q(r, t) ∂q(r, t) ∂r = 1 − r Bθ(r, t) ∂Bθ(r, t) ∂r ,

are therefore only dependent from the poloidal Bθ, where R is the

Toka-mak’s major radius and r = 0 is the center of the TokaToka-mak’s minor radius. Because of the large aspect ratio, Bθ(r, θ, t) becomes independent of θ (so:

Bθ(r, θ, t) → Bθ(r, t)). Between crashes, the magnetic field diffuses following

Maxwell’s equations: ∂B

∂t = −a

2_{∇ × (∇ × B − S} input),

where Sinput is an extra input signal. This equation simplifies to an

inhomo-geneous cylindrical diffusion equation for Bθ(r, t):

∂Bθ(r, t) ∂t = a 2₍∂ 2_B θ(r, t) ∂r2 + 1 r ∂Bθ(r, t) ∂r ) + f (r, t), with a ∈ R.

The Porcelli model now states that if s(rq=1) ≥ scritical, a crash will occur,

where rq=1 is the surface on which q = 1. In the model, scritical can depend

on several parameters, but is assumed constant in this thesis in order to further reduce the complexity of the model. If s(rq=1) ≤ scritical, Bθ evolves

according to the cylindrical diffusion equation and since s(r, t) depends on Bθ, the following sections will study the diffusion equation, in order to predict

when the critical shear condition is met.

4_{For a more detailed description of this model: see [1].} 5_{Definitions and formulas taken as in [1].}

(10)

5 Diffusion equation in Cartesian coordinates

5.1 Homogeneous diffusion equation

The diffusion equation is a partial differential equation of the form: ∂u

∂t = a

2_∆u, ₍₁₎

with a ∈ R, the diffusion constant.

In Cartesian coordinates this equation reduces to: ∂u(x, t)

∂t = a

2∂2u(x, t)

∂x2 . (2)

This equation can be solved by the method of separation of variables6_{. This}

method assumes that the solution can be written in the form: u(x, t) = X(x)T (t). Equation 2 then reduces to:

1 a2_{T (t)} d dtT (t) = 1 X(x) d2 dx2X(x).

Since the left side of the equation is independent from x and the right side of the equation is independent from t, both sides must be equal to a constant, k. For T (t) the equation becomes:

1 a2_{T (t)}

d

dtT (t) = k. Which has the solution:

T (t) = Cea2kt, with C ∈ R. For X(x) the equation becomes7_:

1 X(x)

d2

dx2X(x) = k.

Which has the solution:

X(x) = Ae

√

kx_{+ Be}−√kx_,

6_{The diffusion equation and the methods for solving this equation are given in [2].} 7_{It may be noted that the assumption that u(x, t) can be written in the form X(x)T (t),}

is justified by the fact that the equation for X(x) is a Sturm-Liouville problem, of which the solution is complete and orthogonal, so that every function of x can be expressed as a linear combination of X(x). In the case of the diffusion equation, the coefficients of X(x) are a function of t, sinceR_−∞∞ u(x, t)X(x)dx is independent from x. For more information on Sturm-Liouville problems: see [2].

(11)

with A and B ∈ R. The final solution then, is of the form: u(x, t) = (Ae

√

kx_{+ Be}−√kx_)Cea2_kt .

The solution is bounded if k ≤ 0. If k = 0, u(x, t) is a constant. For k < 0, name k = −λ2_{, with λ ∈ R, with which the solution becomes:}

u(x, t) = e−λ2a2t(Eeiλx+ F e−iλx),

with E, F respectively AC and BC. The specific solution is determined by the boundary and initial conditions of the problem. Important is to note that the function X(x) is periodic, so there exist infinitely many λ that satisfy a given boundary condition. The unique solution is a linear combination of all these solutions, so:

u(x, t) =X

n

e−λ2na2t_(E

neiλnx+ Fne−iλnx).

5.2 Inhomogeneous diffusion equation

The inhomogeneous diffusion equation is a partial differential equation of the form:

∂u ∂t = a

2

∆u + f, (3)

with a ∈ R. In Cartesian coordinates this reduces to: ∂u(x, t)

∂t = a

2∂2u(x, t)

∂x2 + f (x, t). (4)

This equation can be solved by means of the homogeneous diffusion equation with the same boundary conditions. As noted, the Xn(x) form a complete

set, so f (x, t) can be written as P

nfn(t)Xn(x), with fn(t) =

R f (x,t)Xn(x)dx

R X2

n(x)dx . Since the Xn(x) are complete, also u(x, t) can be written as

P

nTn(t)Xn(x),

from which follows, together with _X(x)1 d2_dxX(x)2 = −λ2 ⇔

d2_X(x)

dx2 = −λ2X(x) and equation 4, that:

dTn(t)

dt + λ

2_a2_T

n(t) − fn(t) = 0,

which has the solution for Tn(t):

Tn(t) = cne−(λna) 2_t + Z t 0 e−(λna)2(t−τ )_f n(τ )dτ ,

where the cn can be derived from the initial conditions. The final solution

then, is of the form: u(x, t) = X n (cne−(λna) 2_t + Z t 0 e−(λna)2(t−τ )_f

(12)

6 Diffusion equation in cylindrical coordinates

6.1 Homogeneous diffusion equation in cylindrical

co-ordinates with θ-independence

In cylindrical coordinates equation 1 becomes: ∂u(r, θ, t) ∂t = a 2 (∂ 2_{u(r, θ, t)} ∂r2 + 1 r ∂u(r, θ, t) ∂r + 1 r2 ∂2_{u(r, θ, t)} ∂θ2 ). (6)

This equation reduces to: ∂u(r, t) ∂t = a 2 (∂ 2_{u(r, t)} ∂r2 + 1 r ∂u(r, t) ∂r ), (7)

when u(r, θ, t), is not a function of θ, so u(r, θ, t) → u(r, t). This simplification is justified by the fact that Bθ is also θ-independent. This equation can be

solved with the same method of separation of variables, assuming that u(r, t) can be written as T (t)R(r). Equation 7 then reduces to:

1 a2_{T (t)} d dtT (t) = 1 R(r)( d2 dr2R(r) + 1 r d drR(r)).

Since the left side is independent from r and the right side is independent from t, both must equal a constant, −k2_{. For T (t), the equation becomes:}

1 a2_{T (t)}

d

dtT (t) = −k

2_.

T (t) = Ce−(ak)2t, with C ∈ R. For R(r) the equation becomes8:

d2 dr2R(r) + 1 r d drR(r) + k 2_{R(r) = 0.}

R(kr) = AJ0(kr) + BY0(kr),

with J0 and Y0 being Bessel-functions of zeroth order and A and B ∈ R. The

final solution then, is of the form:

u(r, t) = (EJ0(kr) + F Y0(kr))(e−(ak)

2_t ),

(13)

with E and F respectively AC and BC. The solution is bounded if k > 0. The specific solution is determined by the boundary and initial conditions of the problem. Important is to note that the function R(r) is periodic, so there exist infinitely many k that satisfy a given boundary condition. The unique solution is a linear combination of all these solutions, so:

u(r, t) =X

n

(EnJ0(knr) + FnY0(knr))(e−(akn)

2_t

). (8)

6.2 Inhomogeneous diffusion equation in the sawtooth

model

Between crashes, so for s(rq=1) ≤ scritical, the diffusion equation in the

Toka-mak is of the form9_:

Diffusion equation : ∂u(r, t)

∂t = a 2 (∂ 2_{u(r, t)} ∂r2 + 1 r ∂u(r, t) ∂r ) + f (r, t),(9) Boundary conditions : u(L, t) = H_∂u(0,t)

∂r = 0

, (10)

Initial conditions : u(r, 0) = Bint(r), (11)

with a and H ∈ R, and Bint an undefined function of r. This equation

can be solved by, first, substitution, and then, by using the solution of the homogeneous diffusion equation with the same boundary conditions. At first the substitution u(r, t) = v(r, t) + H is made, with which equations 9, 10 and 11, can be rewritten as:

Diffusion equation : ∂v(r, t) ∂t = a 2 (∂ 2_{v(r, t)} ∂r2 + 1 r ∂v(r, t) ∂r ) + f (r, t),(12) Boundary conditions : v(L, t) = 0_∂v(0,t) ∂r = 0 , (13)

Initial conditions : v(r, 0) = Bint(r) − H. (14)

Since the Rn(r) form a complete set, f (r, t) and v(r, t) can be written as

a linear combination of them. It follows that f (r, t) can be written in the form P nfn(t)Rn(r), with fn(t) = RL 0 f (r,t)Rn(r)rdr RL 0 R2n(r)rdr . Similarly, v(r, t) = P

nTn(t)Rn(r), which, together with d

2 dr2R(r) + 1 r d drR(r) + k 2_{R(r) = 0 ⇔} d2 dr2R(r) + 1 r d drR(r) = −k

2_{R(r) and equation 12, leads to the conclusion:}

Tn(t) = cne−(λna) 2_t + Z t 0 e−(λna)2(t−τ )_f n(τ )dτ ,

(14)

the cn being derived from condition 14: cn = RL 0 (Bint− H)Rn(r)rdr RL 0 R2n(r)rdr .

From the conditions 13 and equation 8, it follows that Rn(r) has the form

EnJ0(knr), with kn = α_Ln, where αn is the nth zero of the J0-function. This

leads to the final result for u(r, t):

u(r, t) = X n cne−(kna) 2_t + Z t 0 e−(kna)2(t−τ )_f n(τ )dτ J0( αn L r) + H, (15) fn(t) = RL 0 f (r, t)J0( αn Lr)rdr RL 0 J 2 0( αn Lr)rdr , (16) cn = RL 0 (Bint(r) − H)J0( αn Lr)rdr RL 0 J 2 0(αLnr)rdr . (17)

(15)

7 Research

7.1 Method

The aim of the research was to investigate which influence some specific parameters would have on the period of the sawtooth oscillation. Because of the (numerical) complexity of the sawtooth model, the simpler, Cartesian model has been studied in this research. As can be seen in equation 5 with no extra input, every next term decreases faster than the previous one, since 0 < λ1 < λ2 < . . ., according to Sturm-Liouville’s theorem. Therefore, first the

system was studied in which only the first term was taken into account, then the more complex system was considered, in which the calculations involved the first ten terms. Throughout the research, the input, f (x, t), was supposed of the form εsin(2πϕt), because of the sinoidal currents that run through the plasma and are used as extra input10_{. The critical shear condition of this}

system was taken: if _∂x∂B(xB=1, t) ≤ scritical, a crash occurs. After the crash,

a reset takes place, in which the function is replaced by another solution of the diffusion equation with a new initial condition. Throughout this research, the new initial condition is set the same as the first initial condition. The reset does not take place for the extra input, and the input will, therefore, have no discontinuous behavior after the crash. The boundary conditions of the system were taken B(0, t) = B(1, t) = 0, from which follows that λn = nπ. In all cases, the first twenty sawtooth periods were calculated

using Mathematica. The calculated periods are also represented in graphs, made in Origin.

Figure 3: Visualization of the critical shear condition in this system. On the intersection between the line y = 1 and the function, the derivative of the function should not become less than a certain value, scritical.

(16)

For these conditions, the first situation, in which only the first term is taken into account, can partly be done analytically. Equation 5 reduces to:

B(x, t) = (c1e−(πa)

2_t +

Z t

0

e−(πa)2(t−τ )f1(τ )dτ )(E1eiπx+ F1e−iπx),

which, in turn, further reduces to: B(x, t) = (c1e−(πa) 2_t + Z t 0 e−(πa)2(t−τ )f1(τ )dτ )sin(πx).

Since the integral f1(t) =

R εsin(2πϕt)sin(πx)dx

R sin2_(πx)(x)dx can be done analytically ([3]), this answer is then reduced to:

B(x, t) = G(t)sin(πx), with: G(t) = c1e−(πa) 2_t +2ε π (

(πa)2_{sin(2πϕt) − (2πϕ)cos(2πϕt) + (2πϕ)e}−(πa)2_t

(πa)4_{+ (2πϕ)}2 ).

First be noted that:

B(x, t) = G(t)sin(πx) = 1 ⇒ x = arcsin(

1 G(t))

π ,

the critical shear condition becomes, that if: ∂

∂xB(xB=1, t) = πG(t)cos(arcsin( 1

G(t))) ≤ scritical,

a crash occurs11_{. The influence of the diffusion constant, a, the critical shear}

constant, scritical, the initial amplitude of the function, c1, the amplitude of

the input, ε, and the frequency of the input, ϕ, on the crash period, τ , was investigated.

Several things can be done analytically for the second situation, in which the first ten terms are taken into account. Equation 5 simplifies to:

u(x, t) =X n (cne−(λna) 2_t + Z t 0 e−(λna)2(t−τ )_f n(τ )dτ )sin(nπx), with: cn= R1 0 Bintsin(nπx)dx R1 0 sin2(nπx)dx ,

11_{This inequality, ≤, was considered an equality, =, and numerically solved for t in}

(17)

and: fn(t) = R1 0 εsin(nπϕ)sin(nπx)dx R1 0 sin 2_(nπx)dx ,

where fn(t) can be evaluated analytically. Though, this situation was

ap-proached mostly numerically, since the calculations were, numerically speak-ing, simpler in their unevaluated form. With the first ten terms taken into account, it is also not possible to obtain an analytical expression for xB=1

and there is, therefore, no other way to express the critical shear condition12. The influence of the diffusion constant, a, the critical shear constant, scritical,

the initial shape of the function, Bint, the amplitude of the input, ε, and the

frequency of the input, ϕ, on the crash period, τ , was investigated.

12_{Also in this case, the ’≤’ in the critical shear condition was replaced by an ’=’ and}

solved numerically for x and t, although the t was the main interest. For the code: see appendix.

(18)

7.2 Results

7.2.1 Dependence of the sawtooth oscillation period on the diffu-sion constant, a

The first conclusion that can be drawn, since two cases were considered, is that only the first term is important for the calculations. The influence of the other terms is, in most cases, a factor 100 or 1000 smaller than the influence of the first term on τ , as can be seen after comparing the tables 1 and 2. Table 1 contains more data points than table 2, due to numerical problems.

a n scritical c1 ε ϕ τaverage τmax τmin

1.00 1 1 5 1 1 0.158329 0.16367 0.153069 1.05 ... ... ... ... ... 0.143586 0.147973 0.139217 1.10 0.131159 0.134566 0.127048 1.20 0.109983 0.112646 0.107251 1.50 0.0705501 0.0715006 0.0691417 2.00 0.0396054 0.0399325 0.0391659 10.00 0.00158187 0.00158193 0.00158181 11.00 0.00130731 0.00130735 0.00130728 12.00 0.0010985 0.00109852 0.00109848

Table 1: The average, maximum and minimum of the first twenty periods for different values of the diffusion constant, a.

a n scritical Bint ε ϕ τaverage τmax τmin

1.00 10 1 5sin(πx) 1 1 0.158444 0.169002 0.148587 1.05 ... ... ... ... ... − − − 1.10 − − − 1.20 − − − 1.50 0.0707772 0.726335 0.0680973 2.00 0.03966 0.0402944 0.0388214 3.00 0.0176819 0.0177239 0.0175834 10.00 0.00158193 0.00158204 0.00158181 11.00 0.00130735 0.0013074 0.00130728

Table 2: The average, maximum and minimum of the first twenty periods for different values of the diffusion constant, a, and n = 10.

The graphs 4 and 5 show exponential behavior. This can well be explained by the fact that in equation 5 the a is indeed situated in the power of the

(19)

Figure 4: The average, maximum and minimum of the first twenty periods for different values of the diffusion constant, a.

Figure 5: The average, maximum and minimum of the first twenty periods for different values of the diffusion constant, a, and n = 10.

exponent. Although the a is also placed in the exponent of the integral, this should be of no influence, since the average of the input is zero. The asymptote of the function should be the a−axis, as suggested in both graphs. The interpretation of these graphs would be, that a larger coefficient, a, causes a shorter oscillation period, which would make the oscillation period

(20)

a function of the density of the plasma.

7.2.2 Dependence of the sawtooth oscillation period on different shapes of the initial funstion, c1 and Bint

a n scritical c1 ε ϕ τaverage τmax τmin

1 1 1 5.00 1 1 0.158329 0.16367 0.153069 .. . ... ... 5.05 ... ... 0.159368 0.164682 0.154143 5.10 0.160399 0.165684 0.155185 5.20 0.162434 0.167655 0.157097 5.50 0.168292 0.17333 0.162717 6.00 0.177172 0.18208 0.171543 7.00 0.192395 0.197646 0.187157 50.00 0.391521 0.396855 0.386369 100.00 0.461969 0.467044 0.456594

Table 3: The average, maximum and minimum of the first twenty periods for different values of the initial amplitude of the function, c1.

Figure 6: The average, maximum and minimum of the first twenty periods for different values of the initial amplitude of the function, c1.

Graph 6 shows exponential behavior, which can also be explained by equation 5, since c1 is the coefficient of the exponent. Again, the extra

(21)

of the input is zero. As can be seen, a larger initial amplitude causes a longer oscillation period, which would make the sawtooth oscillation period a function of magnetich field strength and, therefore, the temperature of the plasma.

a n scritical Bint ε ϕ τaverage τmax τmin

1 10 1 5sin(πx) 1 1 0.158444 0.169002 0.148587 .. . ... ... 5x ... ... − − − 5 0.183419 0.193085 0.173069 5sin(7πx) 0.003330 0.003336 0.003326 5(sin(πx) + sin(7πx)) 0.158444 0.169002 0.148587 5P n=1,3,7sin(nπx) 0.158443 0.169001 0.148585 5sin(20πx) − − − 5P10 n=1sin(nπx) 0.158524 0.169051 0.148696 5P10 n evensin(nπx) 0.040908 0.042542 0.038609 5P10 n unevensin(nπx) 0.158443 0.169001 0.148585 5sin(πx) + e90_sin(7πx) _0.186522 _0.186596 _0.186444

Table 4: The average, maximum and minimum of the first twenty periods for different shapes of the initial function, Bint, and n = 10.

The results for different initial conditions in the case n = 10, are shown in table 4. The test case, 5sin(20πx), linear independent from the ten terms, indeed gave no answer, which enhances the reliability of the calculations. The line 5x could not be calculated due to numerical problems. The other initial conditions listed in this table show that, in most calculations, it is completely justified to only take the first term into account. The exceptions to this rule are the cases where there is no first Fourier component in the initial condition, or the other component(s) have a coefficient of different order of magnitude. Both cases are listed in table 4.

7.2.3 Dependence of the sawtooth oscillation period on the criti-cal shear constant, scritical

The first conclusion that can be drawn, since two cases were considered, is that only the first term is important for the calculations. The influence of the other terms is, in most cases, a factor 100 or 1000 smaller than the influence of the first term on τ , as can be seen after comparing the tables 5 and 6.

The behavior of graphs 7 and 8 is almost linear, but not quite. Expected was exponential behavior, because of the dominant exponent in equation 5, which is not excluded, nor confirmed by means of these graphs. As shown in

(22)

a n scritical c1 ε ϕ τaverage τmax τmin 1 1 1.00 5 1 1 0.158329 0.16367 0.153069 .. . ... 1.05 ... ... ... 0.157839 0.163166 0.152607 1.10 0.15733 0.162642 0.152142 1.20 0.156262 0.161538 0.151207 1.50 0.152698 0.15782 0.147811 2.00 0.145904 0.150594 0.141409 3.00 0.13072 0.134294 0.126339 5.00 0.0991623 0.102105 0.0963364 10.00 0.041328 0.0426847 0.0393537

Table 5: The average, maximum and minimum of the first twenty periods for different values of the critical shear constant, scritical.

Figure 7: The average, maximum and minimum of the first twenty periods for different values of the critical shear constant, scritical.

the graphs, a larger critical shear constant causes shorter oscillation periods. Something not shown in the graphs, is the fact that, contrary to a, c1, ε and

(23)

a n scritical Bint ε ϕ τaverage τmax τmin 1 10 1.00 5sin(πx) 1 1 0.158444 0.169002 0.148587 .. . ... 1.05 ... ... ... 0.157939 0.168459 0.148147 1.10 0.157416 0.167895 0.147716 1.20 0.15632 0.166708 0.146888 1.50 0.152701 0.162725 0.143626 2.00 0.145953 0.155052 0.137581 3.00 0.131125 0.138102 0.122948 5.00 0.0991617 0.104993 0.0937847 10.00 0.0416044 0.04460096 0.0375427

Table 6: The average, maximum and minimum of the first twenty periods for different values of the critical shear constant, scritical, and n = 10.

Figure 8: The average, maximum and minimum of the first twenty periods for different values of the critical shear constant, scritical, and n = 10.

7.2.4 Dependence of the sawtooth oscillation period on the am-plitude of the input, ε

The first conclusion that can be drawn, since two cases were considered, is that only the first term is important for the calculations. The influence of the other terms is, in most cases, a factor 100 or 1000 smaller than the influence of the first term on τ , as can be seen after comparing the tables 7 and 8. Only for relatively large ε, the influence of the other terms on the longest and shortest sawtooth period, can be a factor 10 smaller.

(24)

a n scritical c1 ε ϕ τaverage τmax τmin 1 1 1 5 0.00 1 0.158181 0.158181 0.158181 .. . ... ... ... 0.50 ... 0.158256 0.160883 0.155591 1.00 0.158329 0.16367 0.153069 1.05 0.158336 0.163953 0.152821 1.10 0.158343 0.164237 0.152573 1.20 0.158357 0.164809 0.152079 1.50 0.158398 0.166542 0.150611 2.00 0.158465 0.169497 0.148214 3.00 0.158589 0.175647 0.143592 7.00 0.159002 0.203461 0.127012

Table 7: The average, maximum and minimum of the first twenty periods for different values of the amplitude of the input, ε.

Figure 9: The average, maximum and minimum of the first twenty periods for different values of the amplitude of the input, ε.

The graphs 9 and 10 show linear behavior in the average and minimum period and almost linear behavior in the maximum period. It was well ex-pected that the average period should not have been influenced by a bigger amplitude of the extra input, because of the argument stated above. What is unexpected, is the fact that, although the average is unaltered by bigger amplitudes, the maximum period increases faster than the minimum period decreases. What is also unexpected, is, an effect more explicitly shown in graph 10, that the behavior of the maximum period is almost, but not

(25)

en-a n scritical Bint ε ϕ τaverage τmax τmin 1 10 1 5sin(πx) 0.00 1 0.158181 0.158181 0.158181 .. . ... ... ... 0.50 ... 0.158319 0.16343 0.153267 1.00 0.158444 0.169002 0.148587 1.05 0.158456 0.169576 0.148131 1.10 0.158468 0.170154 0.147677 1.20 0.158491 0.171319 0.146775 1.50 0.158558 0.174887 0.14412 2.00 0.158659 0.181058 0.139848 3.00 0.158834 0.193996 0.131821 7.00 0.161143 0.277033 0.104593

Table 8: The average, maximum and minimum of the first twenty periods for different values of the amplitude of the input, ε, and n = 10.

Figure 10: The average, maximum and minimum of the first twenty periods for different values of the amplitude of the input, ε, and n = 10.

tirely, linear and the behavior of the minimum period is perfectly linear. It seems that in this system, there is a substantial difference between the effect of ε on the maximum and minimum period, but the cause of this difference is not clear from equation 5. A solution could be that for this frequency and these intervals, the input is slightly more often a positive function, but further research is required. What is also shown, in graphs 11, 12, 13 and 14, is how the sawtooth oscillation correlates with the oscillation of the input. The correlation is stronger for higher values of ε, which is to be expected

(26)

Figure 11: A graph showing period of the sawtooth oscillation and the os-cillation of the input, for a = 1, n = 1, scritical = 1, c1 = 1, ε = 1 and

ϕ = 1.

Figure 12: A graph showing period of the sawtooth oscillation and the os-cillation of the input, for a = 1, n = 1, scritical = 1, c1 = 1, ε = 7 and

ϕ = 1.

(27)

Figure 13: A graph showing period of the sawtooth oscillation and the oscil-lation of the input, for a = 1, n = 1, scritical = 1, Bint = 5sin(πx), ε = 1

and ϕ = 1.

Figure 14: A graph showing period of the sawtooth oscillation and the oscil-lation of the input, for a = 1, n = 1, scritical = 1, Bint = 5sin(πx), ε = 7

and ϕ = 1.

7.2.5 Dependence of the sawtooth oscillation period on the fre-quency of the input, ϕ

The results for ϕ are depicted in graph 15 and 16. The first conclusion that can be drawn, since two cases were considered, is that only the first term is important for the calculations. The influence of the other terms is, in most cases, a factor 100 or 1000 smaller than the influence of the first term on τ , as can be seen after comparing the tables 9 and 10. Only for very small ϕ, the influence of the other terms on the longest and shortest sawtooth period, can be a factor 10 smaller. Further, these graphs show steeply increasing behavior for small frequencies, something that can be explained by the fact

(28)

a n scritical c1 ε ϕ τaverage τmax τmin 1 1 1 5 1 0 0.158181 0.158181 0.158181 .. . ... ... ... ... 0.5 0.159309 0.164215 0.152415 1.00 0.158329 0.16367 0.153069 1.50 0.158389 0.162776 0.153811 2.00 0.158411 0.162363 0.154409 3.00 0.158365 0.161431 0.155367 5.00 0.158199 0.160063 0.156333 6.00 0.158354 0.15977 0.156673 6.50 0.158419 0.159649 0.156967 7.00 0.158099 0.15954 0.156827 8.00 0.158214 0.159387 0.156993

Table 9: The average, maximum and minimum of the first twenty periods for different values of the frequency of the input, ϕ.

Figure 15: The average, maximum and minimum of the first twenty periods for different values of the frequency of the input, ϕ.

that the extra input starts as a positive function, and will remain so for a long time for low frequencies. This extra positive input adds to the amplitude and therefore, for these low frequencies, the period is elongated. The graphs also show exponentially decreasing and increasing behavior for the maximum

(29)

a n scritical Bint ε ϕ τaverage τmax τmin 1 10 1 5sin(πx) 1 0.00 0.158181 0.158181 0.158181 .. . ... ... ... ... 0.50 0.160282 0.170139 0.147378 1.00 0.158444 0.169002 0.148587 1.50 0.158581 0.16715 0.149985 2.00 0.158639 0.166346 0.15109 3.00 0.158564 0.164476 0.152953 5.00 0.158298 0.161695 0.154784 6.00 0.158901 0.161102 0.155777 6.50 0.157721 0.160845 0.155761 7.00 0.157949 0.160628 0.155773 8.00 0.15821 0.160305 0.156079 Table 10: The average, maximum and minimum of the first twenty periods for different values of the frequency of the input, ϕ, and n = 10.

Figure 16: The average, maximum and minimum of the first twenty periods for different values of the frequency of the input, ϕ, and n = 10.

and minimum period, something that is not easily shown from equation 5, but is intuitively clear. The asymptote for the maximum and minimum period should be the horizontal line with the value of the period of the system with no extra input, as suggested in both graphs. Another feature,

(30)

nicely depicted in both graphs, is the concept of driving frequencies13. The natural period of this system, as can be read in tables 9 and 10, is 0.158181, which corresponds to a natural frequency of 0.63 . . .. The function should be slowed down by an input with a frequency slightly less than the natural frequency and accelerated by an input that has a frequency slightly higher than the natural frequency. As stated, this feature can be seen near the point ϕ = 0.63 . . ..

(31)

8 Discussion

The results of this research are reasonably satisfying: there are few to none unexpected results; most of the graphs behave sensibly; very few research needs to be done to complete this research on this system. The main thing that will be under discussion, is the usefulness of this research, since the calculations were done in different coordinates and with different conditions than the original, cylindrical problem.

First be noted that serious attempts have been made to tackle the original problem, but these attempts failed, probably due to the numerical complexity of the final solution (equation 15) and the critical shear condition. Second: be noted that the usefulness of the original model might not be guaranteed as well, because all the simplifications, due to the large aspect-ratio of the Tokamak, are not analytical. These simplifications may, therefore, very well be inapplicable to bigger Tokamak reactors, like the one used in the ITER-project, but may not even be applicable to the smaller system considered in this research. The whole model presents us, at best, the qualitative features of the sawtooth oscillation.

If the original model is assumed to give a good, qualitative description, the simpler, Cartesian model could, probably, give this as well. Although this research was done in different coordinates and with different conditions than the cylindrical problem, the global behavior of the functions, and therefore the graphs, is expected to be the same, due to the dominant exponent in both equation 15 and 5. Therefore, the graphs of a, scritical, c1, ε and ϕ are

expected to look the same for the same type of critical shear condition. As can be seen, the critical shear condition of the original and the studied system are alike, but not the same. First of all, the function q(r, t) is proportional to _B(r,t)r , and therefore the condition in the studied problem, B(x, t) = 1, is different from the condition q(r, t) = 1. Secondly, the function s(r, t) is not only proportional to _∂r∂B(r, t), but also proportional to _B(r,t)r , which makes the condition in the studied problem, _∂x∂ B(x, t) ≤ scritical, different from the

original condition. Further research needs to be done, with a different critical shear condition, to improve the applicability of the simpler, Cartesian model. A few things need to be noted considering the research. In experiments, the electric currents are also used to increase the temperature of the plasma14_,

and, therefore, the influence of ε and ϕ on the sawtooth period is more complex than assumed in this research. Since it was not completely clear from this research how the sawtooth oscillation period depends on scritical,

further research needs to be done on the behavior of scritical.

(32)

9 Conclusion

If we assume the original model to give a good, qualitative description and, further, a fair similarity between the original and the studied problem, a few different methods can be concluded from this research to manipulate the sawtooth period. The first way to manipulate this period is to vary the diffusion constant, a. If doable, this might be a very efficient way to manipulate the period, because of the exponential behavior, shown in graphs 4 and 5. The way this might be achieved is to vary the density and/or the temperature of the plasma, which, in turn, can be achieved by varying the applied magnetic field.

Another way to manipulate the sawtooth period is to vary the value of scritical or c1, which might also be done by manipulation of the density and

temperature of the plasma.

The amplitude of the currents that run through the plasma is, based on this research, of no specific influence on the average sawtooth period. Because of the influence this amplitude has on the minimum sawtooth period, this amplitude can be used to enhance the refueling of the burning plasma core. Experiments might lead to different results because of the influence of the electric currents on the temperature of the plasma.

The frequencies of the currents that run through the plasma can also be used for manipulation of the sawtooth period. These frequencies can be selected, so as to function as driving frequencies of the sawtooth period. This method might be the simplest for applications, because the variation of the temperature or density of the plasma could have different effects on different parameters. This method might, therefore, be the most applicable method that can be concluded from this research. Again, experiments may lead to slightly different results because of the influence of the electric currents on the temperature of the plasma.

(33)

10 Acknowledgements

The author would like to thank Menno Lauret. Thank you very much for all the patience and the help you supported me with. Also Peter Schall is sincerely thanked for the trust given in the project and suggestions made.

(34)

11 Bibliography

[1] G. Witvoet (2011). Feedback control and injection locking of the sawtooth oscillation in fusion plasmas. Ph. D. thesis, Eindhoven University of Technology.

[2] E. Kreyszig (2006). Advanced engineering mathematics. John Wiley & Sons, ninth edition.

[3] R. A. Adams & C. Essex (2010). Calculus: a complete course. Pearson Addison Wesley, seventh edition.

[4] D. J. C. MacKay (2009). Sustainable energy - without the hot air. UIT Cambridge Ltd., first edition.

[5] J.Cullerne a.m.o. (2009). Dictionary of physics. Oxford University Press, sixth edition.

[6] D. De Lazzari (2011). Stabilization of magnetic islands in tokamaks by localized heating and current drive: a numerical approach. Ph. D. thesis, Eindhoven University of Technology.

(35)

12 Appendix

12.1 Mathematica code

Bx[t_, \[Epsilon]_, \[Phi]_, A_, c_, L_, teind_] := Re[(Pi/L)*(A*\[ExponentialE]^(-t*(c*Pi/L)^2) +

(2*\[Epsilon]/Pi)*

FullSimplify[Integrate[

\[ExponentialE]^(-(t + teind - \[Tau])*(c*Pi/L)^2)* Sin[2*Pi*\[Phi]*\[Tau]], {\[Tau], 0, t + teind}]] )* Cos[ArcSin[ Re[1/(A*\[ExponentialE]^(-t*(c*Pi/L)^2) + (2*\[Epsilon]/Pi)* FullSimplify[ Integrate[

\[ExponentialE]^(-(t + teind - \[Tau])* (c* Pi/L)^2)*

Sin[2*Pi*\[Phi]*\[Tau]],

{\[Tau], 0, t + teind}]] )]]]]

Bxn[x_, t_, n_, c_, L_, Bbegin_, fQ\[Tau]_, teind_] := Sum[( (N[Integrate[Bbegin*Sin[Q*p*Pi/L], {Q, 0, L}]]/ N[Integrate[Sin[Q*p*Pi/L]^2, {Q, 0, L}]])*\[ExponentialE]^(-t*(c*p*Pi/L)^2) + Integrate[ (N[Integrate[fQ\[Tau]*Sin[Q*p*Pi/L], {Q, 0, L}]]/ N[Integrate[Sin[Q*p*Pi/L]^2, {Q, 0, L}]])* \[ExponentialE]^

(-(t + teind - \[Tau])*(c*p* Pi/L)^2), {\[Tau], 0, t + teind}]

)*Sin[x*p*Pi/L], {p, 1, n}]

Sxn[x_, t_, n_, c_, L_, Bbegin_, fQ\[Tau]_, teind_] := ReplaceAll[

(36)

Evaluate[D[Bxn[y, t, n, c, L, Bbegin, fQ\[Tau], teind], y]], y -> x]

tmax[x_] := Max[Table[t[i] - t[i - 1], {i, 1, x}]] tmin[x_] := Min[Table[t[i] - t[i - 1], {i, 1, x}]] tave[x_] := (1/x)*t[x]

Modeling the sawtooth oscillation in Tokamaks with a hybrid diffusion equation