Optimization of the plasma start-up of the Seville spherical tokamak (SMART)

University of Seville Optimization of the plasma start-up of the Seville spherical tokamak (SMART) Supervisors: Dr. Carlos Soria del Hoyo & Dr. Eleonora Viezzer Jes´us Poley Sanju´an p[email protected] Abstract The SMall Aspect Ratio Tokamak (SMART) is currently being designed by the Plasma Science and Fusion Technology Group (PSFT) of the University of Seville. Therefore, several simulations and calculations regarding the different aspects of the reactor are needed. SMART is an experimental fusion reactor based on magnetic confinement, which has an aspect ratio (A) of >1.6, that classifies it as a spherical tokamak. Moreover, it has a central solenoid, which induces the toroidal electric field to the plasma, two divertor coil sets and two poloidal field coil sets. During this bachelor thesis simulations of the plasma startup in SMART has been carried out. A study optimizing the current of the solenoid and poloidal field coils has been performed. Furthermore, an analysis of the suitable wall width for the outer, top and bottom walls studying the main parameters involved in the startup process such as the connection length (L) is conducted. Finally, the necessary electric field to have a reliable breakdown is studied. Resumen El grupo ”Plasma Science and Fusion Technology Group” (PSFT) de la Universidad de Sevilla est´a actualmente dise˜nando el reactor de fusi´on ”SMall Aspect Ratio Tokamak” (SMART). Por lo tanto, se necesitan gran cantidad de simulaciones y c´alculos sobre los diferentes aspectos del reactor. La relaci´on de aspecto (A) en SMART es de >1,6, lo que lo clasifica como un tokamak esf´erico. Adem´as, tiene un solenoide central, que induce el campo el´ectrico toroidal al plasma, dos pares de bobinas divertoras y dos pares de anillos de campo poloidal. A lo largo de este trabajo de fin de grado se llevan a cabo distintas simulaciones del arranque del plasma en SMART. Se hace a su vez un estudio para optimizar la corriente de las bobinas del solenoide y del campo poloidal. Asimismo, se realiza un an´alisis de sobre cu´al es el grosor adecuado para las paredes exteriores, superiores e inferiores, estudiando los principales par´ametros que intervienen en el proceso de arranque, como la longitud de la conexi´on (L). Por ´ultimo, se estudiar´a el campo el´ectrico necesario para tener un ”breakdown” que sea fiable. Contents 1 Introduction 7 1.1 Tokamaks ..................................... 8 1.1.1 SphericalTokamaks............................ 10 1.2 Stellarators .................................... 11 1.3 TheSMARTproject ............................... 12 2 Theoretical Background 15 2.1 Motion of charged particles in a homogeneous magnetic field . . . . . . . . . 15 2.2 Magneticfluxfunction .............................. 16 2.3 Hoop force and tire tube force . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.4 Physics of the startup phase . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4.1 Breakdown and avalanche . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4.2 Burn-through phase . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 2.4.3 Ramp-upphase.............................. 21 3 Numerical Background 23 3.1 Introduction to COMSOL Multiphysics . . . . . . . . . . . . . . . . . . . . . 23 3.1.1 PIDControl................................ 25 3.1.2 Field line integration . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 3.1.3 Comparison with FIESTA . . . . . . . . . . . . . . . . . . . . . . . . 27 4 Simulations and Results 29 4.1 Estimation of the toroidal electric field value . . . . . . . . . . . . . . . . . . 29 4.2 PF1 and Solenoid current optimization . . . . . . . . . . . . . . . . . . . . . 31 4.3 Optimization of the walls width . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.3.1 Outer wall with optimization . . . . . . . . . . . . . . . . . . . . . . 36 4.3.2 Optimization of width of top and bottom walls . . . . . . . . . . . . 40 4.3.3 Conclusion on the optimization of the width of the wall . . . . . . . . 42 4.4 Electric field optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.4.1 Conclusion on the electric field optimization . . . . . . . . . . . . . . 48 5 Summary and Conclusions 51 A Appendix Magnetic field produced by the axial current of a finite helical solenoid 52 Bibliography 54 Acknowledgments 57 Chapter 1 Introduction Our current lifestyle necessitates large amounts of energy, which will soon turn into a serious problem as the energy sources we are relying on will soon disappear. The traditional ways of producing energy e.g. burning coal and oil, have been proven to be extremely pollutant. Therefore, the seek for alternative ways of generating energy is one of the main aims of nowadays research. In the last decades, substantial effort has been placed into developing renewable energy sources. However, this energy source may not be enough to satisfy the energetic necessity of the actual society. Hence, an alternative way of producing clean energy comes from nuclear fusion. Nuclear fusion is based on the opposite idea of nuclear fission. In the case of nuclear fission technology the energy is produced by the splitting of heavy atoms like uranium in controlled chain reactions. The main disadvantages of fission energy are the residuals of the process, which are radioactive and long-lasting, and the fuel shortage. These two disadvantages do not show in the case of nuclear fusion. On the other hand, in nuclear fusion two light atoms e.g. hydrogen, are fused together to form a heavier one (helium). The process of fusion is depicted in Fig.(1.1). When the temperature of a material is raised, its state changes from solid to liquid and then to gas. As the temperature is elevated further, the gas is ionized and become a high temperature gaseous state, also called plasma, in which the charge numbers of ions and electrons are almost the same and charge neutrality is satisfied on a macroscopic scale. Moreover, when the ions and electrons move collectively, the Coulomb force, which is a long range force and decays with the square of the distance, mediates the interaction between them. The resultant current flow is a result of the Lorentz force and the motion of the charged particles. The hot plasmas used in fusion are referred to this high temperature ionized gaseous state with macroscopic charge neutrality and collective interaction between the particles [2]. As shown in figure 1.1, deuterium (one proton plus one neutron) combines with tritium (one proton, two neutrons) to form helium and a free energetic neutron [1]. During this process the total mass reduces and as a consequence there is an energy release as kinetic energy. The reaction is described by the following equation. 1.1. TOKAMAKS Figure 1.1: Schematic of the fusion process. Protons (p) are shown in red and neutrons (n) in blue [1]. 1D2+1T3→2He4(3.5 MeV) + 0n1(14.1 MeV) [3],(1.1) where the energies shown are the kinetic energies of the reaction products. The total energy released is 17.6 MeV. The reaction is induced by collisions of the particles. Hence, the cross section for the reaction is fundamental. Due to quantum mechanical tunneling the deuterium-tritium (D- T) fusion occurs at lower energies than that required to overcome the Coulomb barrier [3]. Therefore, as it can be seen in Fig.(1.3) the cross section of the D-T fusion has a higher cross section than D-D and D-He3. Nevertheless, the predominant fusion process in the stars is proton-proton. The gas comprised of hydrogen isotopes is called a plasma. There are three ways to confine the plasma. Gravitational confinement which is based on large gravitational forces. This process is the one used by stars. Inertial confinement compresses the hydrogen gas through a controlled implosion, with inertia then holding the gas together long enough for fusion reactions to occur [1]. The third concept is magnetic confinement which uses magnetic fields to exert a force on the moving particles of the plasma which is completely ionized. This type of confinement is the basis for tokamaks and stellarators, which are the main types of nuclear fusion reactors. Moreover, the energies necessaries to overpass the Coulomb barrier (Fig.(1.2)) lead to the need of having high temperatures. As it can be appreciated in Fig.(1.3), the maximum cross section occurs at 100 MeV. Therefore, the plasma will form at high temperatures. 1.1 Tokamaks The magnetic confinement is a result of the particles following the magnetic field lines as shown in Fig.(1.4). The motion of the charged particles in a homogeneous magnetic field and zero electric field is called Larmor motion, which results in a spiral motion around the magnetic field lines. The radius of the motion is derived from the balance of the centrifugal 8 Chapter 1 University of Seville Figure 1.2: Schematic of the potential energy as function of the nuclear separation [3]. Figure 1.3: Cross sections for the reactions D-T, D-D and D-He3[3]. force and Lorentz force, and is the so called Larmor radius [2]. However, there are drifts of the particles out of the path. Nevertheless, as the path has a helical shape the drifts of the outer and inner part are compensated. Fig.(1.5) shows how the combination of the toroidal and poloidal fields results in a helix path that prevents the escape of the particles from the plasma. Furthermore, the toroidal magnetic field is usually of the order of a few Tesla. In order to produce such a field1in a big volume2the coils are wound with many turns of conductors carrying several kilo amperes of current [1]. To avoid the Joule heating loss, which is produced the intensity that goes through those conductor coils, superconductor materials can be used. The magnetic fields applied produce an external magnetic pressure that helps to balance the inner pressure of the gas resulting from the kinetic motion of the particles. Nevertheless, plasma instabilities can emerge destroying the magnetic field. Tokamak is a toroidal plasma confinement device based on magnetic field confinement. Charged particles gyrate around magnetic field lines due to Lorentz force as it is going to be explained below. The cylindrical geometry results in particle losses. Hence, to reduce these losses the helicity shape is needed in a toroidal geometry. In order to produce such a geometry a combination of a poloidal and a toroidal magnetic field is used. Moreover, to balance the plasma pressure the toroidal magnetic field (Bt) and poloidal magnetic field (Bp) are needed. The toroidal component of the magnetic field is typically an order of magnitude bigger than the poloidal one [3]. This poloidal field is mainly produced by the current that goes though the plasma in the toroidal direction while the toroidal filed is produced by external coils. The plasma pressure is proportional to the product of the particle density and temperature. The pressure that can be confined by a certain magnetic field depends on reactor design. In Fig.(1.6) one of the most important tokamak is shown. The Joint European Torus (JET) in Culham, Oxford UK. one large tokamak, is depicted. Moreover, in Fig.(1.6) the 1As a comparison the Earth magnetic field has an order of magnitude of 10−6Tesla. 2In the case of JET the volume is about 200 m3[4]. 2.2. MAGNETIC FLUX FUNCTION Figure 2.1: Schematics of the Larmor motion of a charged particle in a magnetic field [2]. The solution of these equations is a spiral movement (Fig.2.1) around the magnetic field line with the angular velocity ~ω. This movement is also called Larmor motion [2]. The angular frequency ωis called cyclotron (or angular) frequency. Denoting the radius of the orbit by ρω, then the centrifugal force is mv2 ⊥/ρωand the Lorentz force is qv⊥~ B. From the balance of both forces the following relation is obtained: ρω=mv⊥ |q|~ B.(2.3) This is the so called Larmor radius. Furthermore, the center of the Larmor motion is called guiding center1[2]. 2.2 Magnetic flux function A magnetic field line satisfies the equations: dx Bx =dy By =dz Bz =dl B,(2.4) where lis the length along the magnetic field line : (dl)2= (dx)2+ (dy)2+ (dz)2. Furthermore, as it can be seen in Fig.(2.2) the magnetic flux surface ψ(r) = const. is defined as the one on which all magnetic field lines that lie on that surface satisfy following condition [2]: (5ψ(r)) ·~ B= 0.(2.5) The vector 5ψ(r) is normal to the magnetic flux surface and must also be orthogonal to ~ B. Now, rewriting the magnetic field ~ Bin cylindrical coordinates (r, θ, φ) and considering an axi-symmetric configuration (∂/∂θ = 0) the following relation holds: 1If a uniform electric field ~ Eis applied perpendicular to the uniform magnetic field ~ Ba drift motion of its guiding center appears. In [2] a full description can be found. 16 Chapter 2 University of Seville ψ(r, z) = rAθ(r, z),(2.6) where Aθ(r, z) is the magnetic vector potential. The function ψ(r, z) = rAθ(r, z) is called magnetic flux function. Figure 2.2: Schematics of the magnetic surface ψ=const., the normal 5ψand magnetic field line [2]. 2.3 Hoop force and tire tube force The hoop force is present force in the balance of a tokamak. As shown in Fig.(2.3), when the plasma is in a torus, a net force pointing outward the axisymmetric axis exists. This force results from the plasma and the magnetic pressures, i.e. as they push in the toroidal direction this net force appears pointing outward [14]. The hoop force arises due to the toroidal current flowing in the plasma and it is similar to the force produced by a current-carrying loop of wire. Conservation of flux implies that ψ1=ψ2. In Fig.(2.4)(a) is shown the schematics of the hoop force. ψ1denotes the flux of the inner part, it passes through a smaller area than ψ2(the one of the outer part). One expects the average magnetic field on the inner plasma surface S1to be greater than the one on the outer surface S2:B1> B2. Therefore, Bpushes the plasma outward. The net outward force on the plasma along the radial direction is proportional to Fout ∝C(B2 1S1−B2 2S2)/2µ0[10]. In order to balance this outward hoop force, a vertical magnetic field is usually introduced in the tokamak. The toroidal current will interact with this vertical field generating an inward force. Furthermore, there is another force pointing outward that depends on the pressure times the inner and outer surface. This is the so called tire tube force [10]. It arises as follows. The area S1on the inside of the surface is smaller than the area S2on the outside: S1< S2, on a constant pressure contour (which are the magnetic flux surfaces), as it can be seen in Fig.(2.4). If the pressure on this contour is denoted by pthen there is a net outward force along the radial direction proportional to Fout ∝C(pS1−pS2) [10]. Tanking into account that S1< S2then Fout >0, which is a net outward force. 2.4. PHYSICS OF THE STARTUP PHASE Figure 2.3: Schematic of the hoop force [14]. Figure 2.4: Schematics of the outward forces in a system with purely poloidal fields: (a) the hoop force and (b) the tire tube force. [10]. 2.4 Physics of the startup phase In the case of an inductive startup in a tokamak an electric field is applied in order to accelerate the free electrons. Once they have been accelerated they collide, and ionize the neutrals generating more free electrons. These electrons repeat the process leading to an exponential increase in the number of free electrons. The process can be described by the Townsend breakdown/ avalanche[8]. The resulting plasma is conductive with a current Ip. The poloidal magnetic field created by the plasma starts to dominate over the one created by the stray fields. In order to obtain a successful inductive plasma discharge of the radiation barrier coming from the power loss, the Ohmic heating power has to exceed the 18 Chapter 2 University of Seville power loss by radiation and ionization [8]. The loop voltage applied in one turn is given by Vloop =VcoilM/L, where Mis the mutual inductance between the coil and the plasma and Lis the inductance of the coil. Moreover, the electric field at the radius Rcan be described as follows [9]: E=Vloop 2πR.(2.7) The three processes that constitute the start up, i.e. the breakdown, burn-through and ramp-up are going to be described briefly. 2.4.1 Breakdown and avalanche As mentioned above the breakdown2is dominated by the collisions of the free electrons with the neutrals and can be described by using a Townsend-like model [8]. In this type of model the ions are assumed to be stationary3. The temporal evolution of the electron density (ne) can be expressed as follows: dne dt =ne(νion −νloss),(2.8) where νion and νloss are the ionization and loss rate respectively. Hence, the electron density is given by: ne=ne0exp[(νion −νloss)t],(2.9) where ne0is the initial electron density. Equation (2.9) is only valid for a small ionization rate i.e. when the electron neutrals interaction predominates over the Coulomb interaction. Breakdown occurs when the ionization rate is higher than the loss rate. Furthermore, the ionization rate can be expressed using the First Townsend coefficient (α) and the parallel velocity (uk) [8], νion =ukα, (2.10) where α=Apnexp −Bpn Eφ.(2.11) Eφand pnare the toroidal electric field and the neutral gas pressure respectively. The coefficients Aand Bof eq.(2.11) are obtained experimentally. Moreover, during the initial breakdown electrons are lost due to the stray magnetic field created by the stray currents in 2This bachelor thesis will be focused on this part of the startup. 3This approximation is based on the large mass difference between the mass of neutral particles and electrons. 2.4. PHYSICS OF THE STARTUP PHASE the coils and the eddy currents of the vacuum vessel (VV). The loss rate of electrons along the magnetic field lines can be written as follows[8]: νloss =uk Lk ,(2.12) where Lkis the connection length, which is the distance the electron has to travel along the magnetic field lines to impinge the VV wall and are lost for avalanche [8]. The connection length varies in the VV, being larger when the poloidal magnetic field is smaller. In these regions of minimum magnetic field the breakdown is expected to be initiated. Hence, during startup, the poloidal field null regions are created at the desired breakdown location, which means the region that has been chosen to produce the breakdown, at the breakdown time4. The loss rate of eq.(2.12) can be expressed using the average or effective connection length, Leff =rBφ h|Bp|i,(2.13) where h|Bp|i is the volume averaged poloidal field near the null region, defined as a cylinder of radius r[8]. This field can is given by: h|Bp|i ≈ 2 3 r |∇Bp|.(2.14) Therefore, eq.(2.13) can be rewritten as follows: Leff ≈3 2 Bφ |∇Bp|.(2.15) Finally, setting the right hand side of eq.(2.9) to zero, the condition for the onset of the avalanche is obtained [8]: Apnexp −Bpn Eφ=2 3 |∇Bp| Bφ (2.16) Therefore, from eq.(2.16) it can be concluded that the successful breakdown depends on the neutral gas pressure (pn), the toroidal electric field (Eφ) and the poloidal field gradient |∇Bp|. The pressure and the electric field have to be selected accordingly, i.e. if the pressure is too low the electrons are lost before they can ionize a neutral, and if the pressure is too high their mean free path is too short and collide before getting accelerated. 4In fact, in this bachelor thesis the solenoid and poloidal field coils current will be optimized to produce the null region. 20 Chapter 2 University of Seville 2.4.2 Burn-through phase The burn-through phase begins at the end of the avalanche. During this phase low-Z impurities5, usually originating from the VV walls [9], can limit the temperature and the current ramp-up leading to the failure of the discharge. This radiated power (Prad) is expressed as, Prad =neΣnzf(Z, Te),(2.17) where nzis the impurity density, f(Z, Te) is the cooling rate from impurities and the sum is over the impurity species [9]. Moreover, the power available is limited by the power supplies to E2/η for low Teand by j2ηat high Te, where ηis the plasma resistivity, Ethe electric field and jthe current density. Hence, there must be some power left to increase the plasma current (Ip) or the discharge will fail cooling and collapsing [9]. As shown in Fig.(2.5) the peak cooling occurs at 8 and 20 eV for C and O, respectively. Moreover, the cooling rate for beryllium is ten times smaller than for the case of C and O. The use of metal walls can reduce the low-Z impurities. Auxiliary heating can be used to increase the power available to compensate for losses from low-Z radiation [9]. Figure 2.5: Steady state cooling rate due to impurity radiation from Be, C, and O as a function of Te[9]. 2.4.3 Ramp-up phase The plasma current ramp-up phase starts after the burn-through phase of the main gas but is independent of the impurity radiation. Hence, the burn-through and the ramp-up phase can overlap [8]. In order to obtain the successful ramp-up, the stability of the plasma position and the energy and particle confinement of the plasma have to be considered. In order to adjust the plasma position, the poloidal field (PF) coils are used. Furthermore, as the plasma current increases the hoop force increases too leading to the need of a feedback control of the plasma position. The edge safety factor (qedge) has to be kept above 2 [8]. 5Z denotes the atomic number of the material used. 2.4. PHYSICS OF THE STARTUP PHASE A disruption can appear when the current profile is too broad (the internal plasma inductance, liis too low) or when the current profile is too peaked, i.e. there are too strong gradients (litoo high). This is also related with the electron temperature (Te), as Terises the current penetrates less the plasma resulting in instabilities (low li) [9]. The plasma conductivity increases with the temperature as T3 2([3]) leading to the need of a higher diffusion time through the plasma needed. Therefore, controlling the plasma inductance6a disruption can be avoided [8]. 6This can be done choosing accordingly the current ramp rate [8]. 22 Chapter 3 Numerical Background In this chapter the numerical background used to carry out the simulations and results presented during this thesis are going to be reviewed briefly. Moreover, a comparison to FIESTA, a reference equilibrium code developed in Culham (UK) is presented. 3.1 Introduction to COMSOL Multiphysics The software used to carry out the simulations presented in this bachelor thesis is COMSOL Multiphysics, a finite elements software [13]. The process to create the model is described in the flowchart shown in Fig.(3.1). The main part of the interface i.e. the model builder, is shown in Fig.(3.2). It mainly consists of the parameters section, where all the parameters used in the study will be set, the component, which represents the physical object, the study section where type of study can be set and its parameters (time range, sweeps, etc) and the results section where data sets and different plots can be created. Moreover, inside the component its physical properties can be selected. First, the definitions section is introduced, where the variables and functions are defined. Second, the geometry is created by importing it from an external file or by creating it with the tools provided in the software. Third, the materials which constitute the objects are selected from the library or one can import them. Then, the physics will be set accordingly. In this section the coils are added with their respective current and size. Moreover, the differential equations as well as the events for the field line integration1can be added in this section. Finally, the mesh is created defining its size and distributions along the model. 1The field line integration will be explained later in more detail. 3.1. INTRODUCTION TO COMSOL MULTIPHYSICS Figure 3.1: Flowchart of COMSOL Multiphysics. Figure 3.2: Model builder interface in COMSOL Multiphysics. 24 Chapter 3 University of Seville In Fig.(3.3) the geometry used in the model to carry out the simulation is shown. The squares located outside of the VV are the poloidal field (PF) coils one and two and the divertor coils. The solenoid constitutes of 210 coil turns. The VV is divided in three main regions to separate the part where the plasma will be placed. Inside the central region there is a small square of 14 cm2which is the null region i.e. the region where a minimum of poloidal field is produced. Two horizontal lines are drawn to join this area with the walls in order to create the mesh. On the right, the revolved plot of the VV and the coils including the solenoid can be seen. The walls are made of non-magnetic steel 316 and the coils are made of copper. In the simulations presented here, no gas is inserted in the VV as the plasma is not considered. (a) Geometry used for the simulations. (b) Revolved geometry. Figure 3.3: Geometry used to carry out the simulations. 3.1.1 PID Control In order to carry out the optimization of the solenoid and poloidal field coils current a feedback control feature provided by COMSOL is used. The process to apply this tool is shown in Fig.(3.4). It is a Proportional Integral Derivative (PID) control based on solving the equation shown in the following equation [15]. u(t) = ubias +kp[cset −c(t)] + kiZt 0 [cset −c(τ)]dτ −kd d dtc(t),(3.1) 4.2. PF1 AND SOLENOID CURRENT OPTIMIZATION Both the PF2 coils current and the Div 2 coils current are constant for all the simulations. They are plotted in Fig.(4.4). (a) PF2 coils current. (b) Div2 coils current. Figure 4.4: PF2 and Div 2 coils current. Setting Etequal to 2.4 (V/m) as it mentioned before, and Bpequal to zero, (the poloidal magnetic field has to be below 10 G in order to achieve sufficiently large values for the connection length during the Ohmic startup), the plots for the currents optimized are obtained in Fig.(4.5) and Fig.(4.6). Figure 4.5: PF1 coils current. Figure 4.6: Solenoid current. As mentioned above, the Div1 coils current is the same as the current of the solenoid. The PF1 coils current plot (Fig.(4.5)) shows that it is almost constant during the first 10 ms. For the solenoid the current decreases continuously reaching -5.5 (kA) at 10 ms. Fig.(4.7)(a) shows the Bpcontour at 0 ms. As shown, the poloidal magnetic field values in the null region are below 10 G. Furthermore, for an ohmic startup with a loop voltage (see section 2)(Vloop) of 3 V (which is approximately the one obtained here as shown in Fig.(4.7)) the Bpvalues are similar to the ones shown in Fig.(18) in [17], i.e. the poloidal magnetic field values in the central region are in the range between 2 and 12 Gauss. In addition, from Fig.(4.7)(b) it can be concluded that the null is maintained during the first 10 ms as Bp 32 Chapter 4 University of Seville (a) (b) Figure 4.7: Bpcontour [Gauss] in the central region for different time points: (a) at 0 ms, (b) at 10 ms. (a) (b) Figure 4.8: EtBt Bpcontour [V/m] at the central region: (a) at 0 ms, (b) at 10 ms. is almost 10 G for the null region. This is difficult to achieve, in fact, from VEST results shown in [19] the Bpvalue increases significantly for 10 ms. As shown in Fig.(4.8) the EtBt Bpcontour has been plotted for different time points. The noise comes from numerical noise. In order to achieve a reliable Ohmic breakdown the following condition must be satisfied, which was derived for the Ohmic breakdown in JET [17]: Et Bt Bp ⩾103V/m (4.1) In Fig.(4.8) it is observed that at 0 ms the condition is achieved for almost the entire 4.2. PF1 AND SOLENOID CURRENT OPTIMIZATION central region of the vessel. Furthermore, this condition is satisfied for the null region at 10ms but not for the central region. Nevertheless, the condition has to be achieved mainly at 0 ms, that is when ignition starts. For the connection length study, the field line integration, (see the chapter 3) starting at 5 different points, (which is calculated and plotted for every time step in order to see its evolution), and the empirical formula [20] are compared. The empirical formula used is defined as follows: L∼0.25 ·aeffBt Bp ,(4.2) where aeff is the distance from the center of the null region to the closest wall [20]. In this case aeff is equal to 0.09 m. (a) (b) Figure 4.9: Images taken for the connection length using: (a) empirical formula, (b) field line integration at every time step. The evolution of L as calculated by the field line integration is shown in figure 4.10 (b). Here, the time corresponds to the simulation time. As shown in Fig.(4.9)(a) in the case of the empirical formula some regions in which values of more than 100 m are obtained. For the case of the null region L is approximately 30 m. Comparing this value to the average of the values given by the field line integration for t≤2 ms it can be seen that they are almost equal. Moreover, for t≥2 ms the connection length grows significantly as is shown in Fig.(4.9)(b). Therefore, from Fig.(4.9)(b) it can be concluded that the estimation derived for L (see section 4.1, L = 57 m) is only valid for t≥2 ms. Finally, the magnetic flux (psi) contour is plotted in Fig.(4.10) at 0 ms and at 10 ms. Figure (4.10) also shows the value of the poloidal magnetic field have been also plotted the value of the poloidal magnetic field for three different points of the null region1. From 1The points are (0.19,0), (0.19,0.07) and (0.19,-0.07) m. 34 Chapter 4 University of Seville (a) (b) Figure 4.10: Psi contour [Wb] plot for different time points: (a) at 0 ms, (b) at 10 ms. these values it can again be concluded that the null is maintained at 10 ms. Moreover, as the Bpvalues are lower at 10 ms the connection length grows with time (Fig.(4.10)(b)). In these simulations the plasma has not been considered. In addition, as the electric field was set to be 2.4 V/m the loop voltage is near 3 V. As shown in Fig.(4.11) it is even bigger than 3 V. This is because in the the region Etis larger than 2.4 V/m. Moreover, its value depends on the radial value taken (as it is given by Vloop = 2πrE). In this case ris chosen to be 0.19 m. This value of ris called Rin referring to the radial value of the last closed flux surface. It is derived from the geometrical radius (Rgeo) minus the minor radius (a) i.e. Rin =Rgeo −a.Rgeo is set to be equal to 0.42 m, and ais derived from the aspect ratio, i.e. a=Rgeo/A where A(the aspect ratio) has also been selected. Figure 4.11: Loop voltage [V] plot. Moreover, as Etis constant Vloop is also constant as shown in Fig.(4.11). 4.3. OPTIMIZATION OF THE WALLS WIDTH 4.3 Optimization of the walls width In order to improve the walls width in SMART, in this section the previous results obtained for the geometry shown in Fig.(4.3) are compared with the ones derived from the new widths, which are from 8 mm to 6 mm for the outer wall and from 15 mm to 12 and 8 mm for the top and bottom walls. 4.3.1 Outer wall with optimization First, the outer wall width is changed from 8 mm to 6 mm. The width is chosen smaller in order to reduce the eddy currents. The contours are plotted for t = 0 ms. Figure 4.12: PF1 coils current optimized for a width of 8 mm (blue) and 6 mm (green) for the outer wall. Figure 4.13: Solenoid current optimized for 8 mm width outer wall and for 6 mm width outer wall. 36 Chapter 4 University of Seville The PF1 coils current reach a lower value for the case of 6 mm width (Fig.(4.12)). As shown in Fig.(4.13) the solenoid currents are almost the same for both cases. This is because they have the same electric field value set for the optimization. To analyze the breakdown conditions the plots of the Bpcontour at 0 ms and the EtBt Bp contour are shown in Fig.(4.14) and Fig.(4.15). (a) (b) Figure 4.14: Bpcontour [Gauss] in the central region at 0 ms: (a) for 8 mm, (b) for 6 mm. (a) (b) Figure 4.15: EtBt Bpcontour [V/m] in the central region at 0 ms: (a) for 8 mm, (b) for 6 mm. As shown in Fig.(4.14)(b) for the case of 6 mm for the outer wall width the values of the poloidal magnetic field are lower in the null region. In fact, for the central region Bpis 1 G lower2. The empirical condition given in expression (4.1) is satisfied in the null region for both wall widths as shown in Fig.(4.15). Nevertheless, in the null region (near the inner wall) 2In the case of the 8 mm wall the value is 6.5 G and for even a bigger region in the 6 mm case Bpis 5.7 G. 4.3. OPTIMIZATION OF THE WALLS WIDTH higher values are achieved for the 6 mm width wall (Fig.(4.15)(b)). Therefore, in the case of 6 mm the breakdown will be more reliable [17]. The connection length is compared in figures (4.16) and (4.17). The results provided by the empirical formula are shown in Fig.(4.16), while the connection length resulting from the field line integration is shown in Fig.(4.17). As shown the values are higher near the inner wall (inside the null region) for the 6 mm case (Fig.(4.16)(b)). In this case the average is closer to the estimation done for the connection length. Hence, the electric field value chosen is more suitable. (a) (b) Figure 4.16: Connection length given by the empirical formula: (a) for 8 mm, (b) for 6 mm. (a) (b) Figure 4.17: Connection length given field by line integration: (a) for 8 mm, (b) for 6 mm. In Fig.(4.17) the connection length based on field line integration starting at 5 different points is plotted. Until 2 ms the connection length is smaller than 30 m which does not suited the estimation done. After this, for 2 ms in the case of the 6 mm the connection length has larger values on average. Based on both the empirical formula and the field line integration, the 6 mm wall will have a longer connection length. 38 Chapter 4 University of Seville Finally, the the flux surfaces with Psi contours are plotted at 0 ms. The Bpvalue is shown for three points inside the null region3. (a) (b) Figure 4.18: Psi contour [Wb] plot at 0 ms: (a) for 8 mm, (b) for 6 mm. (a) (b) Figure 4.19: Psi contour [Wb] plot at 10 ms: (a) for 8 mm, (b) for 6 mm. The values of the poloidal magnetic field confirms that for the case of 6 mm width wall (Fig.(4.18)(b)) Bpis lower. This explains the increase of the connection length shown above. Furthermore, Fig.(4.19) shows that the null is maintained in time. From Fig.(4.18) and Fig.(4.19) it can be concluded that at 0 ms Bpis more than 1 G lower for the 6 mm case and, at 10 ms it is more than 2 G lower4. Therefore, this comparison shows that the 6 mm case for the outer wall width is more appropriate in order to have a successful breakdown. 3As before the points chosen are (0.19,0), (0.19,0.07) and (0.19,-0.07) m. 4Since no plasma is considered the magnetic flux lines almost do not penetrate the vessel. 4.3. OPTIMIZATION OF THE WALLS WIDTH 4.3.2 Optimization of width of top and bottom walls In this subsection, the top and bottom walls width is varied in order to study which width is more suitable. Initially, the width is set to 15 mm. This value is compared with a width of 12 mm and 8 mm for both the top and the bottom walls. For all the simulations the outer wall width is set to 6 mm. The optimized PF1 coils current and the solenoid current are plotted in Fig.(4.20) and Fig.(4.21). Figure 4.20: PF1 coils current optimized for 15 mm, 12 mm and 8 mm. Figure 4.21: Solenoid current optimized for 15 mm, 12 mm and 8 mm. Both the PF1 coils current and the solenoid current are similar for the three wall widths. The jump of the green line in Fig.(4.20) may appear because the optimization starts to collapse at 9 ms. As the top and bottom walls width decrease the PF1 coils current increases. From the connection length plot based on field line integration starting at 5 different points5(Fig.(4.22)) it can be seen that the first two cases (15 mm and 12 mm, Fig.(4.22)(a) 5In this case the empirical formula is not used as no clear difference could be observed. 40 Chapter 4 University of Seville and Fig.(4.22)(b)) show a similar behavior. Moreover, until 2 ms the second case has a bigger connection length on average. However, after the first 2 ms the 15 mm case has the largest values compared to 12 and 8 mm6. (a) (b) (c) Figure 4.22: Connection length plot by field line integration: (a) for 15 mm, (b) for 12 mm, (c) 8 mm. They all share the same legend plotted on (a). Comparing the 15 mm wall width with the 8 mm one, it can be concluded that the connection length decreases significantly. In this case on average at 10 ms the connection length is more than an 80% bigger for the 15 mm case as it can be seen in Fig.(4.22)(c). Furthermore, even at 2 ms the 15 mm case presents a higher value of L on average compared to 8 mm case. Hence, regarding the connection length based on field line integration the width of the top and bottom walls which has the largest values for the 15 mm case. The flux surfaces with Psi contours are plotted at 0 ms and at 10 ms for the three cases in order to see in detail the magnetic flux. Moreover, the poloidal magnetic field value is plotted for three points inside the null region, which will clarify more the differences in the poloidal magnetic field. (a) (b) (c) Figure 4.23: Psi contour [Wb] plot at 0 ms: (a) for 15 mm, (b) for 12 mm, (c) 8 mm. 6On average at 10 ms the 15 mm wall has a 20% bigger L than the 12 mm wall. 4.4. ELECTRIC FIELD OPTIMIZATION 4.4.1 Conclusion on the electric field optimization In order to find the most suitable electric field value for a reliable breakdown in SMART the optimization of the electric field was carried out. In Table 4.3 the results obtained during the study are presented. The poloidal field coils current decreases from -270 A to -500 A when increasing the electric field. The same tendency shows the solenoid coils current, which decreases from -1000 A to -2200 A approximately. Moreover, as shown in Table 4.3 the poloidal magnetic field increases more than 2 G each step (from 2.4 to 5 V/m). In the third case (Et= 5 V/m), the poloidal field reaches 12 G, which is still sufficiently small to achieve a reliable Ohmic breakdown [17]. The connection length decreases on average in the null region (as the poloidal field increases) when increasing the electric field. Electric field value 2.4 V/m 3.5 V/m 5 V/m PF1 Current [A] -270 -350 -500 Solenoid Current [A] -1000 -1800 -2200 Bp[G] at Rin 6.11 8.37 12.2 EtBt/BpCondition Yes Yes Yes Lemp [m] at Null region 50 - 100 10 - 120 4 - 50 Table 4.3: Results for the top and bottom walls optimization at 0 ms. For Et= 2.4 V/m a loop voltage of 3.4 V is achieved, while for 3.5 V/m and 5 V/m the Vloop is 5 V and 7 V respectively as shown in Fig.(4.35). A comparison between these values and the ones of the spherical tokamak GlobusM [21] is done in Table 4.4. Et(V/m) Vloop (V) SMART 3.5 - 5 5 - 7 GlobusM 1.8 - 3.1 4.5 - 8 (a) (b) (c) Figure 4.35: Loop voltage [V] evolution: (a) for 2.4 V/m, (b) for 3.5 V/m, (c) 5 V/m. The differences between SMART and GlobusM arise because in this case the loop voltage has been done for Rin (r = 0.19 m), which is the inner radius, while in the case of GlobusM it has been done with RGeo (r = 0.36 m)[22], which is the geometrical radius. 48 Chapter 4 University of Seville To conclude, based on the field line integration of the connection length, despite that for the three cases two times the Paschen condition is satisfied from 2 ms on, the only one that satisfies the condition from 0 ms up to 2 ms is Et= 5 V/m. Moreover, the increase of the poloidal field is below 10 G when increasing the electric field from 2.4 to 5 V/m, but the loop voltage increase is more than twice the first case. Therefore, in order to have a successful breakdown until 2 ms it would be necessary to set the toroidal electric field equal to 5 V/m. Chapter 5 Summary and Conclusions The main goal of this bachelor thesis was to determine most suitable currents for the solenoid and poloidal field 1 coils, vacuum vessel geometry and toroidal electric field values to achieve a reliable breakdown in SMART, a spherical tokamak currently being designed for the University of Seville. To this end, three main studies were carried out. First, the optimization of the solenoid and PF1 coils current was done using the PID controller provided by COMSOL Multiphysics. The parameters of the PID control were selected by testing which combination would provide the best results. Moreover, in order to do the optimization of the solenoid current an initial estimation of the value of the electric field was performed. The values of the connection length were calculated using field line integration in the central region of the vacuum vessel. The area where the highest values of the connection length were obtained was selected as the null region. This way, considering the Paschen curve for the average value of the connection length in the null region and the condition given in [16] the value of the electric field needed was selected (2.4 V/m). The region where the highest values of the connection length obtained (57 m on average) is in agreement with the region obtained using Matlab for the field line integration. The solenoid and the PF1 coils current optimization leads to a null region of below 10 Gauss which is maintained until t=10 ms. In fact, the poloidal field is lower at 10 ms than at 0 ms. This is remarkable as the poloidal field must be below 10 G in order to achieve the breakdown with an Ohmic startup [17]. The analysis of the connection length shows that the estimation of the connection length is only valid for t≥2 ms. Nevertheless, as the electric field value is higher than 2.4 V/m in the inner part of the null region, the breakdown will take place in that part moving outward due to the hoop and tire tube forces. The loop voltage obtained is higher than for VEST [19] and in the range of the values of GlobusM [21]. It is constant during the whole optimization process as the electric field. The optimization of the walls was carried out comparing the new values with the results obtained for the initial geometry (8 mm outer wall, and 15 mm top and bottom walls). First, the outer wall was optimized, comparing the 8 mm width (initial value) with the 6 mm wall width. The currents obtained for the solenoids and PF1 coils were very similar to the initial ones (∼ −280 A for the PF1 coils and ∼ −1000 A for the solenoid at 0 ms). The poloidal field contours and the psi contours show that the poloidal magnetic field decreases more than one Gauss in the central region of the vacuum vessel. This shows that the empirical condition for a reliable breakdown in an Ohmic startup is satisfied at 0 ms. Furthermore, the studies on the connection lengths show that there is an increase in the null region reducing the electrons lost in the breakdown [8]. Therefore, the breakdown conditions are more reliable for the case of the 6 mm width of the outer wall. In order to improve the top and bottom wall widths, taking the new outer wall width, a study analyzing the optimal width of those walls was carried out comparing the 15 mm width with 12 mm and 8 mm. From the connection length plots it can be concluded that the connection length decreases on average when the wall width decreases. Moreover, from the poloidal field values of the psi contours it can be seen that the poloidal field increases 0.5 Gauss in every step which leads to smaller losses during the breakdown favoring the avalanche. Therefore, the most suitable results were obtained for the 15 mm top and bottom walls width (the initial value). Finally, the electric field was varied from 2.4 V/m to 3.5 and 5 V/m in order to observe the differences and find the most suitable value for a reliable breakdown in SMART. The solenoid current decreases significantly as the electric field value increases (from -1000 to -2200 A at 0 ms). Furthermore, as the electric field increases the poloidal field increases from 6.11 G to 12.2 G at Rin (r = 0.19 m ). The connection length decreases on average as the toroidal electric field increases as a result of the increase of the poloidal magnetic field. Therefore, it can be concluded that as the connection length and the poloidal field values do not change drastically, i.e. the connection length is of the same order of magnitude and the difference for Bpis below 10 Gauss, the most suitable case will be the 5 V/m case. This is also supported based on that a higher loop voltage is achieved, favoring the breakdown. In this case the loop voltage will almost reach the limit of GlobusM [21] and will be more than two times the value of VEST [19]. To sum up, the solenoid and PF1 coils current were successfully optimized using a PID controller improving the breakdown conditions. Moreover, the outer wall width was changed from 8 mm to 6 mm as the results obtained indicated a more reliable breakdown for the 6 mm case. Finally, the electric field was studied obtaining the best breakdown conditions for the Et= 5 V/m. 52 A Appendix Magnetic field produced by the axial current of a finite helical solenoid The solenoid axial current is not taken into account in the simulation shown in chapter 4. In order to conclude if this assumption is valid the calculation of the value of the magnetic field created by the axial current of a finite solenoid with the helical distribution shown in Fig.(5.1) is presented here. Figure 5.1: Schematic of a typical helical solenoid of length L, radius a, helical angle α, and N wires [23]. As shown in Fig.(5.1) the total current that goes along the wire is NI0, where I0is the current per wire. It is assumed that the N wires will be equally distributed along the solenoid. The solenoid will be approximated by a current sheet shown in Fig.(5.2). After some algebra [23], the following equation is obtained for the magnetic field created by the axial current. Hθ=JZ 4ξk π(ar)1/2K(k) + (r −a)ξ |(r −a)ξ|λo(φ, k)ξ+ ξ− ,(5.1) where Jz(see Fig.(5.2)) is the current density in the vertical axis, ξis the variable of integration, k2= 4ar/ [ξ2+ (a+r)2], K(k) is the elliptic integral of the first kind, ris the radial coordinate of the point where the field is calculated, ais the radius of the solenoid and λ0is the Heuman lambda function given by the following equation [24]: λ0(φ, k) = F(φ, 1−k) K(1 −k) +2 πK(k)Z(φ, 1−k),(5.2) where Zis the Jacobi zeta function, F(φ, 1−k) and K(k) are the incomplete and complete elliptic integrals of the first kind. The equation was solved using Matlab. The mean value of the magnetic field for points near the central region of the vacuum vessel was 8.1094 ·10−7T. Therefore, as it is 3 orders of magnitude below the poloidal field values obtained from the simulations, it is reasonable not to take it into account. Moreover, the vacuum vessel will make this value even smaller. Figure 5.2: Schematic of the mathematical idealization of the N-wire helical solenoid [23]. 54 Bibliography [1] A. Pironti and M. Walker. Fusion, tokamaks, and plasma control: an introduction and tutorial. IEEE Control Systems Magazine, vol. 25, no. 5, pp. 30-43, Oct. 2005, doi: 10.1109/MCS.2005.1512794. [2] Miyamoto, Kenro. Plasma Physics for Controlled Fusion. Springer Series on Atomic, Optical, and Plasma Physics. [3] Wesson, John. (2004). Tokamaks. Oxford, Clarendon press. [4] T. Coyne1, N. Balshaw,S. Davies JET-EFDA. Experiences and Challenges in Pumping the 200m3 vacuum chamber at JET. Culham Science Centre, Abingdon, OX14 3DB. UK. [5] G. Ambrosino and R. Albanese. Magnetic control of plasma current, position, and shape in Tokamaks: a survey or modeling and control approaches. IEEE Control Systems Magazine, vol. 25, no. 5, pp. 76-92, Oct. 2005, doi: 10.1109/MCS.2005.1512797. [6] Y-K.M. Peng and D.J. Strickler. Features of spherical torus plasmas. 1986 Nucl. Fusion 26 769. [7] Costley A. E. (2019). Towards a compact spherical tokamak fusion pilot plant. Phil. Trans. R. Soc. A.37720170439. [8] SINHA, Joyeeta. (2017). Plasma breakdown and current formation in single core and doublet configurations on TCV. Thesis Nº7825 ´ Ecole Polytechnique F´ed´erale de Lausanne. [9] Mueller, D. (2013). The physics of tokamak start-up. Physics of Plasmas, 20(5), 058101- 05810110. doi:101063/14804416. [10] P. Freidberg, Jeffrey. (2014). Ideal MHD. United States of America, Cambridge University Press, New York. [11] Yuhong Xu. (2016). A general comparison between tokamak and stellarator plasmas. Matter and Radiation at Extremes 1, 192. [12] S. Doyle et al., 31st Symposium on Fusion Technology (SOFT2020). [13] (2020). COMSOL INC. Understand, Predict, and Optimize Physics-Based Designs and Processes with COMSOL Multiphysics. At: https://www.comsol.com/comsolmultiphysics BIBLIOGRAPHY [14] Zheng, Linjin. (2015). Advanced Tokamak Stability Theory. Morgan & Claypool Publishers. [15] Ringh, Magnus. (2020). How to Simulate Control Systems Using the PID Controller Add-In. At: https://www.comsol.de/blogs/how-to-simulate-control-systems-using-the- pid-controller-add-in/ [16] Y. Gribov. (2007). Plasma operation and control et al., Nucl. Fusion 47 S385. [17] B. Lloyd. (1991). Low voltage Ohmic and electron cyclotron heating assisted startup in DIII-D et al., Nucl. Fusion 31 2031. [18] A. Mancini et al., 31st Symposium on Fusion Technology (SOFT2020). [19] YoungHwa An, Jeongwon Lee, HyunYeong Lee, JongGab Jo, Bong-Ki Jung, Kyoung- Jae Chung, Young-Gi Kim, Jungmin Jo, Jeong-hun Yang, Yong-Su Na, T.S. Hahm, Y.S. Hwang. (2015). Plasma start-up design and first plasma experiment in VEST et al., Fusion Engineering and Design 96–97 274–280. [20] ITER Physics Expert Group on Disruptions, Plasma Control, and MHD. (1999). Chapter 8: Plasma operation and control et al., Nucl. Fusion 39 2577. [21] V.K. Gusev. (2001). Plasma formation and first OH experiments in the Globus-M tokamak et al., Nucl. Fusion 41 919. [22] Golant, V.E. & Gusev, Vasily & Minaev, Vladimir & Novokhatsky, A.N. & Podushnikova, K.A. & Sakharov, N.V. & Shakhovetz, K.G. & Uzlov, V.S. & Belyakov, V.A. & Divavin, V.A. & Kasatkin, A.I. & Kavin, A.A. & Korotkov, V.A. & Krylov, V.A. & Kostzov, Yu.A. & Malkov, A.A. & Soikin, V.I. & Utin, Yu.L. & Dvorkin, N.Ya & Peng, Y.-K.M.. (1995). GLOBUS-M. From concept to engineering. 1464 - 1467 vol.2. 10.1109/FUSION.1995.534501. [23] A. R. Sass and James C. Stoll. (1993). MAGNETIC FIELD OF A FINITE HELICAL SOLENOID. NASA TECHNICAL NOTE D-1993. [24] Wolfram Research. (2020). At: https://mathworld.wolfram.com/HeumanLambdaFunction.html 56 Chapter 5 University of Seville Acknowledgments First of all, I would like to express my sincere gratitude to my supervisors Dr. Carlos Soria Hoyo and Dr. Eleonora Viezzer for giving me the opportunity to do this bachelor thesis and the necessary help and feedback to develop it. A special thanks to Daniel L´opez Aires for all the support, feedback and knowledge that he has given to me. Huge thanks to the person who has been enduring all my eccentricities and bachelor thesis talks for hours, without whom this thesis would not be possible, my friend and love, Laura. Last but not least, I would like to express my deepest gratitude to my parents and sister who have always supported me and motivated me to follow my path whatever it may be. Without them nothing would be possible and I would be nothing. Seville, August 29, 2020