Neutral Beam Injector design for the SMall Aspect Ratio Tokamak (SMART) of the University of Seville

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo de Fin de Máster Máster en Ingeniería Aeronáutica Neutral Beam Injector design for the SMall Aspect Ratio Tokamak (SMART) of the University of Seville Author: Lina Velarde Gallardo Supervisors: Manuel Toscano Jiménez Juan Francisco Rivero Rodríguez Dpto. Física Aplicada III Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Trabajo de Fin de Máster Máster en Ingeniería Aeronáutica Neutral Beam Injector design for the SMall Aspect Ratio Tokamak (SMART) of the University of Seville Author: Lina Velarde Gallardo Supervisors: Manuel Toscano Jiménez Profesor Titular Juan Francisco Rivero Rodríguez Culham Fusion Research Fellow Dpto. Física Aplicada III Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 Trabajo de Fin de Máster: Neutral Beam Injector design for the SMall Aspect Ratio Tokamak (SMART) of the University of Seville Autora: Lina Velarde Gallardo Tutores: Manuel Toscano Jiménez y Juan Francisco Rivero Rodríguez El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: A mi hermano, por mantener siempre el espíritu bien alto y ponerme una sonrisa en la cara cada día. Resumen 1. INTRODUCCIÓN Y FUNDAMENTO TEÓRICO 1.1. Importancia de la fusión como fuente de energía El aumento incesante de población, el creciente acceso a electricidad en los países en vías de desarrollo y el cambio a industrias más dependientes en la energía está llevando a un aumento drástico de la energía consumida cada año en el mundo, como se observa en la Figura 1. Puede comprobarse cómo la energía consumida procede en más del 80% de combustibles fósiles, que contribuyen al 75% de la emisión de gases de efecto invernadero. Por tanto, se hace urgente la inversión en otras fuentes de energía, que sean sostenibles y abundantes en todas las partes del mundo. Figure 1 Consumo de energía por fuente en el mundo. En este punto entra en juego la fusión nuclear, siendo una fuente de energía: • Sostenible. Sus combustibles son el deuterio y tritio (obtenido del litio), que están disponibles en todas partes del mundo. El deuterio se puede extraer del océano, y el litio es abundante en la corteza terrestre. •Sin emisión de CO2. Por tanto, no contribuye a la emisión de gases de efecto invernadero. III XResumen Se ha llevado a cabo asimismo una comparativa entre el porcentaje de partículas que chocan con la pared para cada una de las versiones, así como del shine-through, como se muestra en la Tabla 2. Table 2 Comparación numérica de los resultados obtenidos con ambas versiones de BBNBI. Pared FILD Shine-through BBNBI5 0.0148% 75.5% 0.15% BBNBI4 0.0027% 100.0% 0.16% De este modo, se ha concluido que el NBI modelado mediante BBNBI5 es muy similar al obtenido con la versión anterior y presentado en la citada tesis. Figure 7 Histograma de la coordenada r para las distintas versiones de BBNBI y las dos configuraciones del NBI de MAST-U. 3.2. SMART Como se ha mencionado anteriormente, los inputs para la pared, el campo magnético, y los perfiles de densidad y temperatura del plasma fueron generados para esta sección. En cuanto a los marcadores a simular, se han llevado a cabo dos estudios distintos. En primer lugar, se simularon cantidades pequeñas de partículas para realizar un estudio de sensibilidad, comprobando que los inputs se habían generado adecuadamente. De esta forma, es posible almacenar la órbita completa en ASCOT, y analizar los resultados obtenidos. Una vez se comprobó que la dirección seguida por las partículas y por las derivas era la esperada, se procedió a generar los marcadores partiendo de las partículas inyectadas a través de los distintos sistemas NBI testeados, con cantidades más altas (50000 partículas, concretamente). La Figura 9 muestra la disposición del sistema NBI de SMART, mostrando el parámetro de optimización elegido, el radio de tangencia (R t ). Dicho radio ha sido variado hasta encontrar la configuración que minimizaba las pérdidas y maximizaba la corriente inducida. El resultado óptimo se muestra en la Figura 10, donde se pueden observar las proyecciones poloidal y toroidal del Resumen XI Figure 8 Histograma de la coordenada z para las distintas versiones de BBNBI y las dos configuraciones del NBI de MAST-U. haz inyectado mediante el NBI. Para la configuración elegida, se han realizado distintos análisis, estudiando el perfil de deposición entre otras distribuciones, como la de energía frente a ángulo de inclinación. En la Figura 11 se muestran tres órbitas representativas para cada condición de parada: aquellas que finalizaron la simulación por alcanzar la mínima energía al colisionar con otras partículas, y las que chocaron con la pared. Figure 9 Vista toroidal de la configuración del NBI de SMART. Los círculos negros representan las paredes y las rojas la proyección toroidal de la separatrix. La línea amarilla simboliza el haz inyectado. 4. CONCLUSIONES Y TRABAJO FUTURO Como se ha comentado anteriormente, la fusión se posiciona hoy en día como la posible solución a los problemas energéticos en la Tierra. El sistema de confinamiento magnético es hoy en día la mejor configuración, con reactores como ITER a pocos años de su primer plasma, o JET llevando a cabo su segunda campaña de deuterio-tritio en la actualidad. Simultáneamente, otras configuraciones están siendo estudiadas en reactores en todas partes del mundo, entre las que destacan los tokamaks esféricos, ámbito dentro del cual se engloba este trabajo. Se ha trabajado con distintas herramientas, como MATLAB para realizar simulaciones numéricas XII Resumen (a) Vista poloidal. (b) Vista toroidal. Figure 10 Distribución de nacimiento de las partículas generadas mediante el NBI. de las órbitas seguidas por un protón sometido al campo magnético de distintos artefactos. Se ha proseguido con los códigos de Monte-Carlo BBNBI y ASCOT, que han permitido en primer lugar la implementación del NBI de MAST-U en la última versión de ASCOT, ASCOT5, así como la familiarización con los parámetros de ambos códigos. Esto ha permitido generar de cero los inputs necesarios para SMART, el cual no había sido simulado antes en ASCOT5. El diseño del NBI de SMART se ha centrado principalmente en minimizar las pérdidas, provenientes tanto de shine-through como de prompt losses, aunque se han tenido en cuenta otros parámetros como la corriente inducida o la distribución de slowing-down. En la Tabla 3 se muestran los porcentajes totales de pérdidas para las distintas configuraciones estudiadas, pudiéndose concluir que se ha obtenido el diseño óptimo buscado. Resumen XIII (a) Parada por finalización de la simulación. (b) Parada por colisión con la pared. Figure 11 Órbitas representativas en el plano poloidal para dos condiciones de parada distintas. Table 3 Porcentaje total de pérdidas para las distintas configuraciones estudiadas. ID configuración Radio de tangencia Shine through Prompt losses Pérdidas totales 1 0.515m 20.89% 0.01% 20.90% 2 0.461m 16.10% 0.02% 16.12% 3 0.408m 13.49%0.03%13.52% 4 0.355m 12.55% 1.11% 13.66% 5 0.303m 11.98% 8.46% 20.44% 6 0.252m 12.79% 19.81% 32.60% Los siguientes pasos en el diseño del NBI será probar configuraciones que incluyan variar la altura a la que se inyecta, alejándose del plano ecuatorial, y asimismo incluir nuevos parámetros a estudiar, como el par de torsión inducido por el NBI. Finalmente, se realizará el mismo estudio con los parámetros correspondientes a la tercera fase de operación de SMART, que cuenta con perfiles de densidad considerablemente mayores a los usados en este trabajo. Así, se ha comprobado que es posible reducir el shine-through en un 80-90 Abstract Several alternatives are being developed to create a device that allows the energy production through the fusion process. One of the most promising studied configurations is the tokamak, a device that, using a set of magnetic field coils in different configurations –toroidally and poloidally–, keeps the charged particles forming the plasma confined, at high enough temperatures to achieve fusion. These temperatures are attained through different methods, among which is the Neutral Beam Injection (NBI). This system injects neutral particles in the plasma, heating it up via coulomb collisions, and generating fast ions. The correct confinement of these ions is of vital importance, as they do not only maintain the plasma temperature, but also induce current. Furthermore, if a high number of particles is lost, they can pose a threat to the integrity of the plasma facing components of the wall. The Plasma Science and Fusion Technology (PSFT) group at the University of Sevilla is currently building the SMall Aspect Ratio Tokamak, SMART. The goal of this master thesis is the design of the NBI system for this device. For this purpose, three main numerical tools have been used. First, MATLAB has been used to simulate the orbits of protons subject to the magnetic field generated through different devices, such as an infinite wire or a combination of a coil and a toroid, in order to achieve helical orbits similar to those of tokamaks. The other two tools are BBNBI and ASCOT, both developed by the Aalto University, in Finland. The first allows the calculation of the ionization profiles of the particles injected via the NBI. The latter uses those particles as input in simulations that follow their trajectory inside the tokamak, allowing a detailed study of the fast-ion losses. In this way, these tools have been used for the optimization of the NBI design for SMART. For the design, the data for the magnetic equilibrium and the temperature and density profiles have been provided by the PSFT group, corresponding to the second operational phase of SMART. The obtained design has reduced the total losses to 13.52%, for an injection energy of 15keV, a very satisfactory result. In the future, the use of the data corresponding to the third phase of operation is expected to lead to a reduction of the losses by 80-90%. XV Contents Resumen III Abstract XV List of Figures XVIII List of Tables XXI 1 Introduction 1 1.1 Importance of fusion as an energy source 1 1.2 Present of nuclear fusion: Tokamaks 2 1.3 Alternative designs: Spherical Tokamaks 5 1.3.1 MAST Upgrade 6 1.3.2 SMART 6 1.4 Motivation and scope of the project 8 2 Theoretical Background 11 2.1 Physics of tokamaks 11 2.2 Particles motion 13 2.2.1 Drift caused by an Inhomogeneous Magnetic Field 14 2.2.2 Particle orbits 15 2.3 Fast-ions in tokamaks 16 2.3.1 Neutral Beam Injection 17 Optimisation parameters 17 Beam neutralisation 19 2.3.2 Fast Ion Loss Detector 19 3 Numerical tools 21 3.1 Preliminary simulations in MATLAB 21 3.1.1 Infinite wire 21 Numerical resolution 23 3.1.2 Coil 23 Effect of the initial velocity in the OZ axis 27 Effect of the initial velocity in the OY axis 30 Effect of the initial velocity in the OX axis 33 Conclusions 34 3.1.3 Toroid and Coil 35 XVII XVIII Contents First scenario: α=1,β=037 Second scenario: α=0,β=138 3.2 BBNBI 40 3.3 ASCOT 41 3.3.1 Inputs and outputs 42 4 ASCOT simulations for the design of the NBI system: MAST-U and SMART 45 4.1 The NBI injector script: from ASCOT4 to ASCOT5 BBNBI benchmark in MAST-U 45 4.1.1 MAST-U NBI parameters 45 4.1.2 Verification of the BBNBI5 MAST-U model 45 4.2 SMART: Following orbits with ASCOT 51 4.2.1 Inputs 51 4.2.2 Sensitivity study 52 4.3 SMART: NBI design 58 4.3.1 Resulting preliminary design 59 5 Summary and future work 63 5.1 Summary 63 5.2 Future work 64 References 65 Agradecimientos 69 List of Figures 1 Consumo de energía por fuente en el mundo III 2 Principales elementos de un tokamak V 3 Primer plasma del reactor MAST-U VI 4 Esquema de SMART y la configuración de sus imanes VII 5 Esquema 3D del toro y la espira (en rojo) VIII 6 Vista 3D de las órbitas en el interior del toroide IX 7 Histograma de la coordenada r para las distintas versiones de BBNBI y las dos configuraciones del NBI de MAST-U X 8 Histograma de la coordenada z para las distintas versiones de BBNBI y las dos configuraciones del NBI de MAST-U XI 9 Vista toroidal de la configuración del NBI de SMART. Los círculos negros representan las paredes y las rojas la proyección toroidal de la separatrix. La línea amarilla simboliza el haz inyectado XI 10 Distribución de nacimiento de las partículas generadas mediante el NBI XII 11 Órbitas representativas en el plano poloidal para dos condiciones de parada distintas XIII 1.1 Energy consumption by source in the world 1 1.2 Evolution of the cross-section with the temperature 3 1.3 The interior of the JET tokamak 4 1.4 Upending of the first vacuum vessel sector 5 1.5 Spherical tokamaks in comparison to the conventional doughnut-shaped tokamaks 6 1.6 Different divertor configurations that the MAST Upgrade device can use 7 1.7 First plasma achieved by MAST-U in the Culham Centre for Fusion Energy 7 1.8 Sketch of SMART and the coil configuration 8 2.1 Different sketches of the tokamak configuration and coordinate system 12 2.2 Poloidal cross section of a tokamak showing flux surfaces 13 2.3 Surface S used to define the poloidal flux ψin Equation 2.4 13 2.4 Poloidal view of the orbit topology in a tokamak 16 2.5 Magnetic field lines of a simple axisymmetric magnetic mirror 16 2.6 Sketch of an NBI system 18 2.7 Evolution of the neutralisation efficiency with energy for systems based on positive and negative ions 19 2.8 CAD design of the FILD1 in the ASDEX Upgrade tokamak 20 XIX 2Chapter 1. Introduction This is where nuclear fusion comes in [4]: •Sustainability . Fusion relies on deuterium and tritium (obtained from lithium) as fuels, which are available all around the globe. Deuterium can be obtained from the ocean, and lithium is an abundant metal on Earth’s crust, making them an almost infinite resource. •No CO2emission . This means no contribution to increasing the greenhouse gases that go into the atmosphere. •No risk of meltdown . Accidents like Fukushima cannot happen in a fusion device. High amounts of energy are needed to overcome the nuclei repulsion in order to achieve fusion. Hence, the reaction would stop if the energy flux stops. Furthermore, as the quantity of fuel inserted in the device is only enough for a few seconds, there is no risk of a chain reaction. If any disturbance occurs, the plasma will rapidly cool down, stopping the reaction. Finally, operations are carried in ultra-high vacuum, so any problem with the device would break the vacuum and the plasma would stop its course. •No long-lived radioactive waste . No high activity, long-lived nuclear waste is produced. Tritium has a half-life of around 12 years. All this aspects have made fusion the perfect option to replace fossil fuels as the world’s primary energy source. However, as of today, scientists have not been able to achieve energy production through nuclear fusion, remaining a complex challenge. 1.2 Present of nuclear fusion: Tokamaks As it will be further discussed in Chapter 2, thermonuclear fusion is a process where energy is produced due to the combination of two light nuclei (reactants) into a heavier element, releasing a large amount of energy following Einstein’s equation, E=mc2 , as the products of the reaction have a total mass slightly lower than the reactants. Fusion is widely known as the energy of the stars, as it is in fact the source that powers them. In the Sun, hydrogen atoms fuse together to form helium at very high temperatures, at which the particles get ionised, and their nuclei together with electrons form a gas of charged particles called plasma . These temperatures are necessary for the nuclei to overcome the electrostatic repulsive forces and get close enough to collide and hence fuse, leading to a release of energy. On Earth, elements need to be kept at temperatures of the order of 150·106◦C and confined for an amount of time large enough that will allow the nuclei to fuse. The reaction produces high-energy neutrons, which are absorbed by blankets that surround the core of the device and get heated up. A coolant element (water, helium or Li-Pb eutectic) is flowed through the blanket to cool it down, collecting the energy which will be used to generate electricity by conventional methods [5]. As it was previously mentioned, one of the main advantages of fusion energy is that it relies on elements that are widely available: deuterium (D) and tritium (T), two isotopes of hydrogen. But there is more than one fusion reaction that has the potential for energy production. Some of them are: 2 1D+3 1T→4He (3.5MeV ) + n(14.1MeV)(1.1) 2 1D+2 1D→3 1T(1.0MeV) + p(3.0MeV)(1.2) 2 1D+2 1D→3He (0.82MeV ) + n(2.45MeV)(1.3) 1.2 Present of nuclear fusion: Tokamaks 3 2 1D+3He →4He (3.6MeV) + p(14.7MeV )(1.4) Each of these reactions can be characterised by its cross-section , σ , a measure of the reaction probability. In Figure 1.2, the cross-sections for 1.1-1.4 are shown. The DT reaction has the larger cross-section, peaking at the lowest temperature. This is the reason why most of nowadays research has focused in this process, even if deuterium and hydrogen are still the most common elements used in current working devices. However, despite the D-T reaction needing the lowest temperature, the elements will still need to be heated up to hundreds of millions of Kelvin, a temperature at which matter turns into its fourth state, the plasma, as it was previously introduced. Figure 1.2 Evolution of the cross-section with the temperature for the different reactions mentioned above. Both DD reactions have very similar probability of occurring, the orange line shows the total cross section for both processes [6]. It was also mentioned that tritium can be obtained from lithium. This will be possible through a process called breeding, that can take place in the reactor itself: when the neutron released from the reaction (Eq. 1.1) hits the walls, which are covered with a layer with embed lithium called the blanket, a fission reaction is induced, producing a tritium atom (Eq. 1.5). However, it will perhaps still be necessary to introduce tritium from outside. n6 3+Li →4 2He +3 1T(1.5) The problem that fusion presents is the need to develop a device that can maintain the plasma confined long enough and at high enough temperatures so that ignition can be achieved. This means that the process becomes self-sustained, as enough fusion reactions occur. Once this milestone is reached, the process will be able to yield more energy than it receives, hence becoming a viable source of energy. According to the European Fusion Roadmap, achieving a positive energy output is planned by 2050, with nuclear fusion potentially contributing to the energy mix in the long term [7]. As of today, devices using two types of approaches are being studied: magnetic and inertial confinement. Magnetic confinement relies on a compound of coils to induce a current to the plasma and keep the particles confined, while inertial confinement consists of compressing a small pellet using strong lasers or beams. However, in this project, the focus will be laid on the first type. In 4Chapter 1. Introduction Magnetic Confinement Fusion (MCF), magnetic fields are used to confine the particles, holding them away from the reactor walls and preventing the plasma from damaging the wall materials. The most effective configuration is a toroidal magnetic field, combined with a poloidal one, resulting in force lines following helical paths. These forces are able to confine and control the plasma. The most spread toroidal confinement systems are tokamaks , but there are other types of devices being studied, such as stellarators. Tokamaks were invented by two Russian physicists, Igor Tamm and Andrei Sakharov [8], and the name stands for TOroidalnaya KAmera ee MAgnetnaya Katushka [5], which in Russian means toroidal chamber with magnetic coils [9]. The process that takes place inside these devices will be further discussed in Chapter 2. To maintain the high temperatures needed, tokamaks use different methods, such as Neutral Beam Injection (NBI). The NBI system design for a tokamak that is being built in the University of Sevilla will be the scope of this thesis. In the 1970s, European countries came together and began the design of the Joint European Torus (JET), which currently is the largest functioning tokamak, and can be found in Culham, Oxford, UK. In 1983 the first plasmas were achieved in this device [10]. Figure 1.3 The interior of the JET tokamak. Also in the 1980s, the General Secretary of the URSS Gorbachev proposed to US president Reagan the idea of a collaborative international project to develop fusion energy for peaceful purposes at the Geneva Superpower Summit: the International Thermonuclear Experimental Reactor (ITER, “the way” in latin) [10]. Meanwhile, during the 1990s, the first experiments using tritium as fuel were carried out in JET, being the first reactor to ever run a 50-50% mix of tritium and deuterium. In 1997, it also set a record for the highest Q ever achieved, at 0.67 [10], where Q is the ratio of thermal fusion power over input power, Q=Pf us/Pin , – i.e., to achieve fusion energy the Q value must be greater than 1 –, when Q=1 , it is called breakeven. The aim of ITER will be to achieve Q≥10 [11] and to test the technologies to produce tritium from lithium. At present, the second ever DT experimental campaign (the first one was in 1997) is taking place in JET, where it is expected to break the record of fusion energy in tokamaks. Furthermore, ITER has reached the assembly phase and is scheduled to deliver its First Plasma for December 2025. In January 2021, project execution to First Plasma was at 72.4% [12]. In Figure 1.4, the upending of the first vacuum vessel sector performed on March 26, 2021 can be observed. After ITER, the next step in the European Roadmap is DEMO, the demonstration reactor that should produce electricity from the fusion process for the first time, while operating with a closed tritium fuel-cycle (it will produce its own fuel). Meanwhile, ways to improve its design using new technology and physics are being studied. The 1.3 Alternative designs: Spherical Tokamaks 5 Figure 1.4 Upending of the first vacuum vessel sector [13]. International Fusion Materials Irradiation Facility (IFMIF) is a projected materials test facility in which candidate materials for the use in an energy producing fusion reactor (such as DEMO) can be fully qualified [14]. IFMIF is composed of two phases [15]: •EVEDA (Engineering Validation and Engineering Design Activity). It will validate the main technical elements of IFMIF through prototypes. •DONES (DEMO-Oriented Neutron Source). A reduced version of IFMIF, coordinated by EUROfusion. The Spanish government approved on March 16, 2021 the fundings for the Spanish candidature to locate the IFMIF-DONES project in Granada, Andalucía [16]. Apart from the mentioned devices, there are around 250 tokamaks throughout the world [17], each of them having made possible different kinds of progress in the nuclear fusion field, and having provided building blocks for the international project that ITER represents. 1.3 Alternative designs: Spherical Tokamaks This project will focus on the design of the Neutral Beam Injection (NBI) system for the tokamak SMART (SMall Aspect Ratio Tokamak) currently being built in Sevilla, Andalucía, by the Plasma Science and Fusion Technology (PSFT) Group [18], [19]. Its specific characteristics will be further detailed in Section 1.3.2. As its name suggests, SMART will be a small aspect ratio tokamak (ratio between the major and minor radius, see Eq. 1.6). These devices are usually referred to as spherical tokamaks (STs), and are characterised for having a fairly small major radius (R) and rather large minor radius (a), as it can be seen in Figure 1.5. A=R a(1.6) where A is the aspect ratio, which is usually below 2 for spherical tokamaks. Spherical Tokamaks also tend to have large plasma elongation, κ , i.e. the ratio between the plasma height and its minor radius. Their small aspect ratio makes the region of the plasma in the area near to the central column to be close enough to the plasma at the other side of the column to be affected by the poloidal magnetic field it generates. This results in a larger κ than that in conventional tokamaks. In fusion, the fusion power increases with κ5 . Hence, one of the main 6Chapter 1. Introduction Figure 1.5 Spherical tokamaks in comparison to the conventional doughnut-shaped tokamaks [23]. advantages of STs is that, with smaller magnetic fields, they can produce the same fusion power as a conventional tokamak, thus leading to smaller construction costs. These characteristics have various advantages [20]: • Being smaller on average allows smaller magnets to be used, and locates them on average closer to the plasma, which has several benefits: it can achieve higher pressures at a given magnetic field, leading to higher output powers (fusion power scales with the square of the pressure). • Spherical tokamaks may also eliminate the need for superconducting magnets, as they allow to stabilize plasmas with magnetic fields that are much smaller than those of non-spherical tokamaks (up to a factor of ten) [20]. •These advanced configurations can also improve stability properties (see [21] and [22]). However, spherical tokamaks present also some disadvantages, as copper coils have to be replaced regularly and dissipate a considerable amount of power, leading to a decrease in efficiency [20]. This innovative design allows a more compact configuration, and can be key in the future stages of fusion power, as they are hoped to be one day used as fusion reactors. 1.3.1 MAST Upgrade The Mega Amp Spherical Tokamak Upgrade (MAST-U) from the Culham Centre for Fusion Energy (CCFE, Oxford) will be used as part of this project to get started with ASCOT simulations in a spherical tokamak. This is a device based on the original MAST, which ran from 2000 to 2013. Its upgrade will explore plasma exhaust with its new divertor configurations (Figure 1.6), such as the Super-X divertor or the Snowflake configuration, in order to remove excess heat from the machine, helping to preserve the materials. It creates a magnetic configuration that spreads the heat loads at the divertor area of the machine, and could be a viable concept for fusion power stations [24]. Present divertors cannot endure the high power loads created in large tokamaks such as ITER, and would require replacement every few years, which would have an important economical impact. If successful, Super-X could be a prototype for commercial fusion reactors [25]. MAST-U achieved its first plasma on 29 October 2020 (Figure 1.7). 1.3.2 SMART SMART is the SMall Aspect Ratio Tokamak that is being built in the University of Sevilla. Its aim is to contribute to the technological and scientific development of fusion devices by testing innovative divertor configurations, plasma facing components and diagnostics, while providing training for 1.3 Alternative designs: Spherical Tokamaks 7 Figure 1.6 Three different divertor configurations that the MAST Upgrade device can use [26]. Figure 1.7 First plasma achieved by MAST-U in the Culham Centre for Fusion Energy. future engineers and physicists. SMART will also study plasma confinement in both positive and negative triangularities (i.e. plasma shaped as a triangle rather than as an ellipse), as the effects that negative triangularity has on plasma properties are still unknown for STs. The equilibria for the baseline configuration and both triangularities have been presented by S. Doyle et al. [19]. The development of SMART will comprise three phases, as complexity increases. The designed NBI will be introduced in the second phase, as shown in Table 1.1, where the main parameters are shown. The magnetic field and plasma current will be achieved by 12 toroidal coils, 8 poloidal coils (4 Poloidal Field coils and 4 Divertor coils) and a central solenoid, as presented in Figure 1.8(a), where the solenoid is displayed in red around the central column of the tokamak, the plasma in pink rotating around it, the poloidal coils in yellow surrounding the plasma and the toroidal coils in orange in the exterior of the device. 8Chapter 1. Introduction Table 1.1 Main parameters for the three phases of the SMART tokamak. The geometrical parameters are maintained constant through all the phases. Phase 1 Phase 2 Phase 3 Major radius, R0[m] 0.4 Minor radius, a[m] 0.25 Aspect ratio, A< 2 Plasma elongation, κ< 2.1 Toroidal magnetic field [T] 0.1 0.3 1 Plasma current, IP[kA] 30 100 400-500 Pulse length, τ[ms] 20 100 500 Microwave heating power, PECRH [kW] 6 (2.45 GHz) 6 (7.5 GHz) 200 NBI heating power, PNBI [kW] - 600 600 0.0 0.5 1.0 Radius R[m] 1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1.00 Height Z[m] 0 100 200 300 RGrid [cells] 0 50 100 150 200 250 ZGrid [cells] 1.6 m 0.8 m Div1 Div2 PF2 PF1 (a) (b) Sol Figure 1.8 Sketch of SMART and the coil configuration. Figure (b) shows a poloidal section with the 2D mesh employed [19], displaying the location of the poloidal coils in yellow and the solenoid, in red around the central column of the tokamak, with the plasma in pink rotating around it. The toroidal coils are shown in orange in the exterior of the device. 1.4 Motivation and scope of the project The goal of this thesis is to obtain a preliminary design of the NBI system that will be built for the tokamak SMART of the University of Sevilla. The NBI plays a crucial role in heating the plasma up to the needed temperatures, and it is therefore one of the key elements of any designed tokamak. In Section 2.3, the parameters that enable the optimisation of such a design, which will be studied in this thesis in order to find the best configuration possible, will be presented. To perform the simulations that will allow this analysis, different numerical tools will be used: BBNBI5 [27], ASCOT5 [28] and MATLAB. Hence, the contents of this project will be structured as follows: In order to become familiar with the equations and physics behind tokamaks, introduced in 1.4 Motivation and scope of the project 9 Chapter 2 , these equations were programmed in MATLAB, where simulations of proton orbits subjected to the magnetic field of different devices were performed, as it will be shown in Chapter 3 . This allowed a better understanding of how the Monte-Carlo code ASCOT works, which will be briefly introduced in Chapter 3, along with the ionisation Monte-Carlo code BBNBI. Then, Chapter 4 will include the different results obtained as part of this thesis. It will first cover the simulations ran in ASCOT using the MAST-U geometry and magnetic field configuration, with the intention to benchmark the NBI modelled with BBNBI5 for this thesis against the model obtained via BBNBI4 in [29]. It will also include the first orbits ever simulated in ASCOT5 for the SMART tokamak, and finally the optimisation of the SMART NBI configuration, where some additional results for the chosen design will be presented. Finally, in Chapter 5 , the outcome of this thesis will be discussed, and the future work that will be carried out in the NBI design will be described. 2 Theoretical Background In this Chapter, the physics behind tokamaks will be introduced, describing the main theoretical concepts that will be used later in this thesis, such as the orbit topology or the main parameters used to characterise those orbits. Finally, the concept of fast ions will be introduced, focusing on the description of the NBI system as the central element of this thesis. 2.1 Physics of tokamaks As it was introduced in the previous chapter, tokamaks are one of the devices used to confine the plasma. The coordinate system in such a device is shown in Figure 2.1a, where θ is the poloidal angle and ϕ the toroidal one. In the tokamak, the magnetic configuration results from the superposition of a toroidal and a poloidal magnetic field, as it can be seen in Figure 2.1b. The toroidal field, Bϕ , is created by a series of toroidal coils evenly spaced around the reactor (in blue), resulting in a magnetic field that depends on the distance to the toroidal axis with 1/R . A transformer coil in the center of the device (in green) induces a current in the plasma in the toroidal direction (green arrow) that generates a poloidal field, Bθ (also in green). This leads to helical magnetic field lines (in black), which are needed to confine the plasma, overcoming the orbit drifts caused by the magnetic gradient, as will be seen later. The outer poloidal field coils, a system of horizontal coils outside the toroidal magnet structure (in gray) contribute to the poloidal field and help shaping and positioning the plasma. The plasma confinement is hence a result of the Lorentz force, as particles follow the magnetic field lines in helical trajectories. In tokamaks, the current flowing through the plasma also serves to heat up the plasma via the Joule effect. This is usually referred to as ohmic heating (OH). However, this method is limited, as it is not effective to heat the plasma up to fusion relevant temperatures, because plasma resistivity decreases with temperature. Thus, additional heating systems are needed to achieve the temperatures necessary for fusion. The most important one is Neutral Beam Injection (NBI), whose design in SMART will be the scope of this project. There are also other systems such as Electron/Ion Cyclotron Resonance Heating (E/ICRH). In addition, the DT reaction produces He atoms, also called α particles, that use Coulomb collisions to transfer their energy as well. Both NBI and ICRH generate fast-ion populations in the plasma (ions that have an energy above the bulk plasma), which transfer part of their energy to the bulk plasma through Coulomb collisions, heating both ions and electrons. The confinement of these fast ions is essential in the performance of a tokamak. Additionally, if not handled correctly, they can pose a threat to the wall materials. 11 18 Chapter 2. Theoretical Background Figure 2.6 Sketch of an NBI system [38]. •Shine through . Percentage of beam power that is directly lost, because it does not get ionised by the plasma, and collides onto the first wall, without contributing to plasma heating. This parameter must be calculated (and minimised), as it is a measure of the total power deposited in the plasma. •Prompt losses . Percentage of ionised neutrals whose first poloidal transit intersects the wall, undesiredly heating it up. This parameter must also be minimised as this power might damage the wall and affect diagnostic components located near it. •Current drive . As providing current drive is one of the objectives of NBI systems, this parameters will be maximised. •Slowing down distribution . This is represented in a histogram that integrates the state of each particle in time, showing how the energy distribution is scattered as the fast ions exchange their energy with the ions and electrons in the plasma, heating them. The goal here will be to locate the pitch of the more energetic ions as close to 1 as possible, in order to have more passing orbits. Another important parameter to take into account is the Alfvén velocity, i.e. the velocity of the Alfvén waves in the tokamak. The ions moving at that velocity will be resonant and exchange high amounts of energy. In order to allow the study of the Alfvén Eigenmodes, the fast ions generated through the NBI should be injected at energies high enough so they will reach this velocity when slowing down. It can be calculated as: vA=Bt √µ0ni (2.22) where Bt is the toroidal magnetic field, µ0 the magnetic permeability and ni the ion density. As its value only depends on the magnetic equilibrium and the plasma density, it will remain constant for the different NBI configurations tested. Hence, this parameter will help set a minimum for the main injection energy, as it is directly related to the velocity. In a first approach, minimising both the shine through and the prompt losses can be considered the main objective of the optimisation. 2.3 Fast-ions in tokamaks 19 Beam neutralisation The beam neutralisation before injection is a fundamental aspect of the NBI system, which can be based both on positive and negative ions (e.g. for a hydrogen beam): •Ionisation by ions. Hb+H+ p→H+ b+H+ p+e− •Ionisation by electrons. Hb+e−→H+ b+2e− •Charge exchange. Hb+H+ p→H+ b+Hp where Hb would be a neutral and the elements with the “+” superscript are ionised particles, whereas the subscripts band prefer to the injected beam and the plasma, respectively. The evolution of the neutralisation efficiency with energy is shown in Figure 2.7. At low beam energies, of the order of 10 keV, both types of ionisation result in very similar efficiencies. However, positive ions are commonly used since they are technologically easier to achieve. Ionisation by electrons becomes more important at high energies, of the order of 500-1000 keV [37], being able to maintain an efficiency of ∼ 0.6 even for the highest energies. Hence, for large, dense plasmas, where high energies are needed, negative ions beams must be used. There is an energy at which the heating rate of the NBI is the same for both ions and electrons, the critical energy [39]. Figure 2.7 Evolution of the neutralisation efficiency with energy for systems based on positive and negative ions [37]. Generally, beams with negatively charged particles during the acceleration process are known as Negative NBI (NNBI) systems, whereas those with positive ions will be Positive NBI (PNBI). As it can be observed in Figure 2.7, the latter process produces also molecular species such as H+ 2 or D+ 3 . This causes the final beam injected to have energetic ions at three different energies: E0,E0/2 and E0/3 , where E0 would be the main injection energy . For the MAST-U NBI, this energy is of 75keV. For the NBI designed in this thesis, however, the energy will never be higher than 40keV, so only PNBI will be considered. 2.3.2 Fast Ion Loss Detector The Fast Ion Loss Detector (FILD) is a scintillator-based detector [40] that is used to measure fast-ion losses, giving information of their velocity space. Its principle of operation is shown in Figure 2.8, where the head of the detector is displayed, with trajectories from lost particles penetrating through the collimator (which only allows certain trajectories to get through) and eventually reaching the 20 Chapter 2. Theoretical Background Figure 2.8 CAD design of the FILD1 in the ASDEX Upgrade tokamak. The blue lines represent the fast-ions trajectories, while the pink one has been blocked by the collimator. The scintillator is shown in green [41]. scintillator plate, which emits light when irradiated by charged particles. Inside the detector, it is possible to see how particles hit the scintillator in different positions, giving a complete information of the velocity space of the lost particles. A FILD detector for MAST-U was recently developed (see J.F. Rivero-Rodríguez’s work, [42]), and the simulations of the detector’s measurements there presented will be used here as a benchmark in order to ensure the quality of the SMART NBI design. 3 Numerical tools Firstly in this Chapter, the equations introduced in the previous one will be implemented in MATLAB to achieve a more comprehensive understanding of the equations behind tokamaks. For this purpose, orbit simulations will be performed, starting with a proton around an infinite wire and finishing with the orbits in a toroid. Then, the main numerical tools used in this thesis will be introduced: BBNBI, that allows the calculation of the ionization profiles of the NBI particles, and ASCOT, that uses those particles as inputs and simulates them until they reach any of the established limits (a maximum simulation time, a minimum energy or the collision with any element of the wall). 3.1 Preliminary simulations in MATLAB As it was already mentioned, all the simulations in this section will use a proton. The equations for each of the cases will be first presented, followed by the simulations and analysis performed. 3.1.1 Infinite wire An infinite wire (OZ direction) with a current I = 5 MA produces a magnetic field B , following the Ampère’s law. Cylindrical coordinates, ρ,φ,z, will be used. The magnetic field produced by a wire can be expressed as: B(ρ) =   0 µ0·I 2π·ρ 0 (3.1) where µ0=4π·10−7NA2. An analysis of the orbit followed by a proton subjected to the resulting magnetic field will be now performed. In order to run the simulation, it is necessary to find the equations of motion in cylindrical coordinates. It is possible to apply Lorent’z Law to calculate the force applied to it. F=q˙ x×B=m¨ x(3.2) First, the position will be derived, in order to find ˙ xand ¨ x. 21 22 Chapter 3. Numerical tools Provided xis x=ρer+zez(3.3) Then, ˙ x=˙ ρer+ρ˙ φeφ+˙zez(3.4) ¨ x= ( ¨ ρ−ρ˙ φ2)er+(2˙ ρ˙ φ+ρ¨ φ)eφ+¨zez(3.5) Inserting both expressions in 3.2, the differential equation to be used can be obtained. q  −˙zB0/ρ 0 ˙ ρ/ρB0  =m  ¨ ρ−ρ˙ φ2 2˙ ρ˙ φ+ρ¨ φ ¨z (3.6) Finally, the differential equations obtained are:          dρ dt dφ dt dz dt d˙ ρ dt d˙ φ dt d˙z dt          =         ˙ ρ ˙ φ ˙z ρ˙ φ2−q m˙zB0/ρ 2˙ ρ˙ φ/ρ q m˙ ρ/ρB0         (3.7) Besides, in order to calculate the energy of the proton in each moment of time, magnitude that will be used to measure the error committed, the next equation will be used: Ec=1 2m(˙ ρ2+ρ2˙ φ2+˙z2)(3.8) Once the equations have been obtained, the initial conditions have been set as follows: v0=c/300 ρ0=1m c : lightspeed ρ=ρ0z=0φ=0 vρ=v0vz=0vφ=v0/10 The simulations will be ran using different integration methods, in order to compare the accuracy of each of them: •Euler’s method. •2nd order Runge-Kutta’s method. •4th order Runge-Kutta’s method. The simulations have been run for 100 cyclotronic periods ( 100τc ), calculated using the initial magnetic field (B0=1T): τc=2πm qB0 (3.9) 3.1 Preliminary simulations in MATLAB 23 Figure 3.1 Comparison between the different integration methods used, for N=107steps. This value of B0 results in a τc=6.549 ·10−8s . Consequently, the movement of the proton will be simulated for 6.549 ·10−6seconds. Numerical resolution The result obtained for the different methods, with a number of steps of 107 ( ∆t=2.0847 ·10−13 s) is shown below in Figure 3.1. All three methods result in the same orbit, as this number of time steps is high enough for all of them to have a good resolution. If a significantly smaller number was used, such as 104 , it would be possible to start observing some changes and inaccuracies, as reflected in Figure 3.2. For this time step ∆t=2.0847 ·10−10 s, both Runge-Kutta methods still work fine and deliver the orbit, maintaining the kinetic energy constant in a cyclotronic period, with errors under one over a thousand. However, if the number of time steps was downgraded even more, to 103 , resulting in a time step ∆t=2.0847 ·10−9 s, Euler’s method completely fails to solve the equations, and the 2nd order Runge-Kutta starts showing some inaccuracies, as it is shown in Figure 3.3. Recall that the cyclotronic period is τc=6.549 ·10−8seconds. Hence, it can be concluded that the 4th order Runge-Kutta’s method is the one that offers the best performance at an equal number of time steps, as expected. 3.1.2 Coil A magnetic field B is produced by a 1-turn coil (OXY plane) that has an intensity I=1 MA and a radius a=1m. Cylindrical coordinates, ρ,φ,z, will be used once again. The magnetic field produced by a coil will be obtained through the Biot-Savart equation, and can be expressed as: B=µ0I 4πIˆut׈ur r2dl (3.10) 24 Chapter 3. Numerical tools Figure 3.2 Comparison between the different integration methods used, for N=104steps. Figure 3.3 Comparison between the two Runge-Kutta methods for N=103steps. where dl is a differential element of current and ˆut and ˆur unit vectors that point in the direction of the current and towards the point where the value of the magnetic field will be evaluated, P, respectively. Figure 3.4 shows a sketch including all these components. As it can be seen in this figure, r is the distance between the differential element of current dl =a·dϕ and the point of interest P. If dl is located at (acosϕ,asinϕ,0)and P at (0,y,z), the distance rcan be expressed as: r=pa2+y2+z2−2aysinϕ(3.11) 3.1 Preliminary simulations in MATLAB 25 Figure 3.4 Sketch of the elements involved in the Biot-Savart equation for a coil [43]. And hence ˆur: ˆur=−acosϕ·ˆ i+(y−asinϕ)·ˆ j+z·ˆ k r(3.12) While ˆutwill be: ˆut=−sinϕ·ˆ i+cosϕ·ˆ j(3.13) So the cross product can be performed, obtaining: ˆut׈ur=  zcosϕ zsinϕ a−ysinϕ (3.14) This leads to the three components of the magnetic field in cartesian coordinates, for any point in the YZ plane. Bx=µ0I 4πaz Z2π 0 cosϕ (a2+y2+z2−2aysinϕ)(3/2)dϕ=0(3.15) By=µ0I 4πaz Z2π 0 sinϕ (a2+y2+z2−2aysinϕ)(3/2)dϕ(3.16) Bz=µ0I 4πaZ2π 0 a−ysinϕ (a2+y2+z2−2aysinϕ)(3/2)dϕ(3.17) where dl has been substituted by a·dϕ. The first integral is of course zero, as for each differential element of current dl in one side of the YZ plane there is another one symmetric on the other side, mutually cancelling their effect. However, due to the symmetry of this problem, there is no need to restrict the field point to the YZ plane [44]. The coordinate in ˆ j can be written as ˆ ρ , as the value of By can be written as a generic radial component, Bρ , as it does not depend of the ϕ angle and Bϕ=0 . Hence, the components of 26 Chapter 3. Numerical tools the Bfield can be written as: Bϕ=0(3.18) Bρ=2µ0I 4πaz Zπ/2 −π/2 sinϕ (a2+y2+z2−2aysinϕ)(3/2)dϕ(3.19) Bz=2µ0I 4πaZπ/2 −π/2 a−ysinϕ (a2+y2+z2−2aysinϕ)(3/2)dϕ(3.20) Expressions obtained by changing the limits of integration and multiplying the integrals by two, profiting from their symmetry. In order to solve the equations analytically, a change of variable is performed, θ=π/2−ϕ: Bρ=2µ0I 4πaz Zπ 0 cosθ (a2+ρ2+z2−2aρsinθ)(3/2)dθ(3.21) Bz=2µ0I 4πaZπ 0 a−ρcosθ (a2+ρ2+z2−2aρsinθ)(3/2)dθ(3.22) Equations 3.21 and 3.22 have now the form of two elliptical integrals [44]: Zπ 0 dφ (b±cosφ)(3/2)=m 2−2m√2mE (3.23) Zπ 0 ±cosφdφ (b±cosφ)(3/2)=√2mK −2−m 2−2m√2mE (3.24) where K and E are the elliptic integrals of the first and second kind, respectively: K=Zπ/2 0 dφ (1−msin2φ)(1/2)(3.25) E=Zπ/2 0 (1−msin2φ)(1/2)dφ(3.26) So equations 3.21 and 3.22 can be rewritten as: Bρ=µ0I 2π az (2aρ)(3/2)(−√2mK +2−m 2−2m√2mE)(3.27) Bz=µ0I 2π a (2aρ)(3/2)(am 2−2m√2mE +ρ√2mK −ρ2−m 2−2m√2mE)(3.28) where m=2 1+band b= (z2+a2+ρ2)/(2aρ). This way, the expressions for both components of the magnetic field are obtained in terms of ρ and z, allowing the calculation of their value at any point in space. However, as the expressions shown above would diverge for ρ=0 , the analytical solution at the OZ axis can be obtained by using the equation: Bz(ρ=0) = µ0I a2 2(z2+a2)(3/2)(3.29) 3.1 Preliminary simulations in MATLAB 27 as there is only field in the OZ direction along the OZ axis. An analysis of the equations used to describe the motion of the proton around the wire will be performed, implementing the necessary elliptical functions in MATLAB. Once the equations have been obtained, the orbit followed by the proton subjected to the effect of the resulting magnetic field will be analysed using cartesian coordinates, x , y , z . The orbit followed by the proton around the wire will be obtained by applying the Lorent’z Law, as it was previously done in Section 3.1.1, but in this case the magnetic field has components in all three space coordinates. The obtained differential equation will be integrated with the Runge-Kutta Method of the 4th order, with N=104 time steps, always checking that the kinetic energy is conserved in time. This time, the simulation will be run for 50 cyclotronic periods (50τc). In these simulations, also the gyro-centre will be displayed. In order to obtained the path followed by the gyrocentre, its position in time, G, will be obtained as: PG =RLuF(3.30) where P is the position of the proton and RL is the Larmour radius, also called the cyclotronic radius, which can be obtained through the second Newton’s law: F=|qv×B|=mv2 ⊥ RL (3.31) On the other hand, uFis the unitary vector that points in the direction of the force: uF=F F=qv×B |qv×B|(3.32) Hence, being both the velocity and the magnetic field known, it is possible to easily obtain the Larmour radius and find the position of the gyrocentre. The position of the orbit in each time step will be referred to as the point P , as previously mentioned. So, in order to plot the gyrocentre, the calculated vector PG will be added to the position of the orbit, OP. Effect of the initial velocity in the OZ axis The initial conditions can now be set in cartesian coordinates: v0=c/1000 where c=lightspeed r(0) = [0,0,−a](3.33) v(0) = v0·[0,1,1](3.34) In this first scenario, the proton will have initial velocity different to zero in both the OY and OZ directions. The results can be observed in Figure 3.5, where the orange line represents the gyrocentre, and the yellow circle is the coil, located in the OXY plane. As the coil has a 1m radius, it is difficult to properly see the orbit described by the proton. In Figure 3.6a, a zoomed representation is shown. Finally, in Figure 3.6b the projection on the XY axis can be observed. 34 Chapter 3. Numerical tools results in a pitch angle of 90o , i.e. all the velocity is perpendicular. Note that the guiding-centre is not presented in this figure for a clearer visual result. The limit value of α can be found at 0.880, for which the orbit is still trapped. The trajectory described by the proton in this case is shown in Figure 3.16c. As it can ben observed, the result differs from the previous one in the length of the banana. Being this the limit situation, the orbit described by the proton is nearly completed. The 2D projection in Figure 3.16d clearly demonstrates that it’s almost a closed orbit. On the other hand, when setting α=0.881 , the orbit changes, and the full orbit is now described, as shown in Figure 3.16e. The 2D projection in Figure 3.13f shows a clearer view of the completed orbit. Conclusions An interesting result that can be extracted from the previous analysis is the pitch angle at which the topology of the orbit changes. Table 3.1 shows a comparison of the different pitch angles for the studied cases. It is possible to see how the extreme cases for the two scenarios (varying the OX or the OY velocity) all show a very similar pitch, which indicates that the orbit topology is only a function of the pitch angle: under ∼50.6o, the orbit will always be passing. Table 3.1 Pitch angle for the different configurations studied. βPitch angle [o]αPitch angle [o] α=1 0 15.59 β=1 0 90 1.140 50.5670 0.880 50.4517 1.141 50.5904 0.881 50.4819 3.1 Preliminary simulations in MATLAB 35 (a) 3D plot of the orbit described around the coil, α=0. (b) 2D projection of the orbit, α=0. (c) 3D plot of the orbit described around the coil, α=0.880. (d) 2D projection of the orbit, α=0.880. (e) 3D plot of the orbit described around the coil, α=0.881. (f) 2D projection of the orbit, α=0.881. Figure 3.16 Orbit of the proton and guiding-centre plotted along the coil for the different initial velocities in the OX axis. 3.1.3 Toroid and Coil A toroid will now be added to the coil used in the previous section. A toroid is composed by a curved solenoid with its ends united, forming a circular solenoid. The direction of the magnetic field produced in such a device is shown in Figure 3.17. The red dots symbolize each coil forming 36 Chapter 3. Numerical tools the solenoid. Figure 3.17 Position of the coils forming the toroid (section in red) and direction of the magnetic field. The value of the magnetic field produced by this toroid can be easily obtained applying Ampère’s law: B=µ0Ni 2πr(3.40) where N is the number of coils and i the intensity in each of them. As seen in Figure 3.17, this magnetic field is oriented in the direction of the angle φ of the cylindrical coordinates. Until now, there was no field in this direction, which is the feature introduced by the toroid. Note that the field is only present inside the torus, by Ampère’s law. By adding these two elements, a combination of toroidal and poloidal magnetic field, resulting in a helical magnetic field, is achieved, similar to that present in tokamaks. The coil will maintain the radius a=1 m and intensity I=1 MA, while the toroid will have an internal radius rset to 0.75m. These parameters are shown in Figure 3.18. (a) Geometrical parameters of the torus. (b) 3D view of the torus and the coil (in red). Figure 3.18 Sketches of the torus. The number of coils around the torus and their intensity will be variable, as the goal here will 3.1 Preliminary simulations in MATLAB 37 be to study the changes in the orbit as the intensity of the toroid increases. For this purpose, a parameter KIwill be used: KI=Ni I(3.41) This parameter will be varied from 0 (only the coil) to 1 (the total intensity of the toroid is equivalent to that of the coil), in order to study the changes in the orbit. Cartesian coordinates will once again be used, x , y , z . In this case, the initial velocity in the OZ axis is set to zero, and the components in x and y will be progressively modified to check the effect on the type of orbit obtained. The simulations will be run for 20 µs , close to the time interval that will be used in the ASCOT simulations. Moreover, using a final simulation time in terms of the cyclotronic period was not possible here as varying the Ni parameter of the torus would vary that property. The number of time steps used is n=7·104 , keeping the error, calculated via the kinetic energy, of the order of 1e−7 . As for the initial conditions , they will be the same ones as in the previous analysis performed, in terms of two parameters, αand β: v0=c/1000 where c=lightspeed r(0) = [a,0,a/5](3.42) v(0) = v0·[α,β,0](3.43) First scenario: α=1,β=0 The first simulations will be run with an initial velocity only in the direction of the OX axis. Figure 3.19 shows the orbit described by the proton in this new magnetic field. This first case corresponds to a passing orbit, as clearly presented in Figure 3.19a, where ρ is again calculated as the square root of the x and y coordinates squared. The passing orbit is expected, as it was the case in the previous analysis, since the initial velocity is mostly parallel (small pitch angle). Both these figures present the results for the orbit produced by the coil alone (KI=0) and with a toroid with KI=0.1 . While the trajectory in the ρ−z plane does not vary much, as it can be seen in Figure 3.20, where a zoom was performed to clearer see the orbit around the initial point, the added intensity of the torus results in a smaller toroidal period, as shown in Figure 3.19b. This is the expected result, as the toroid increases the magnetic field in the toroidal component, and so the orbit travels a longer path in that direction. The comparison between both extreme cases will now be performed. Figure 3.21 shows the orbit described by the proton for both extreme situations considered, KI=0 and KI=1 , and also for an intermediate case, KI=0.5 . Figure 3.21a, shows how the orbit topology remains the same, a passing orbit, as the initial velocity stays unchanged. Once again, the trajectory in the ρ−z plane does not vary much, as it can be seen in Figure 3.21c, where once again a zoom was performed to clearer see the orbit around the initial point. For these three situations, the changes in the X−Y plane projection are even larger than before, as the intensity of the torus has been increased. The orbit described by the extreme case KI=1 covers a quarter of the equatorial plane. This is better shown in Figure 3.21d, where the 3D view of the torus is presented along the three orbits, showing the effect of the toroid intensity. As it was previously introduced, this figure clearly shows how 38 Chapter 3. Numerical tools (a) 2D projection of the orbit in the ρ−zplane. (b) 2D projection of the orbit in the X−Yplane. Figure 3.19 Orbit of the proton for the cases KI=0,0.1, initial velocity in the OX direction. Figure 3.20 Zoom around the initial position of the 2D projection of the orbit in the ρ−z plane, initial velocity in the O direction. the superposition of the coil and the toroid results in helical trajectories, very similar to those in tokamaks. However, it is important to note that, in order to achieve realistic comparisons with tokamaks, configurations with KI>1 – i.e., with a toroid intensity larger than that of the coil, should be analysed. Second scenario: α=0,β=1 The second set of simulations will be run with an initial velocity only in the direction of the OY axis. The results shown in this section will be parallel to those of the previous one. First, the orbits for the KI=0,0.1 case will be presented in Figure 3.22. With the new initial conditions the orbit topology has changed, as shown in Figure 3.22a. Again, this can be explained by looking at the initial velocity, with a pitch angle of 90o , which means that all the velocity is perpendicular at the beginning. In this scenario, the trajectory in the ρ−z plane is still very similar for both configurations shown, but again this is not the case for the X−Y plane. However, the change is not as pronounced as in the previous scenario. The comparison between both extreme cases is also performed. Figure 3.23 shows the orbit described by the proton for both extreme situations, KI=0 and KI=1 , and also for an intermediate 3.1 Preliminary simulations in MATLAB 39 (a) 2D projection of the orbit in the ρ−zplane. (b) 2D projection of the orbit in the X−Y plane. (c) Zoomed 2D projection of the orbit in the ρ−z plane. (d) 3D view of the orbit inside the torus. Figure 3.21 Orbit of the proton for the cases KI=0,0.5,1, initial velocity in the OX direction. (a) 2D projection of the orbit in the ρ−zplane. (b) 2D projection of the orbit in the X−Yplane. Figure 3.22 Orbit of the proton for the cases KI=0,0.1, initial velocity in the OY direction. case, KI=0.5 . Figure 3.23a, shows how the orbit topology is again trapped. Despite the trajectory in the ρ−z plane not varying much, in this case the orbits for the largest KI go beyond the initial point, as it can be seen in Figure 3.23c, where a zoom was performed to clearer see the orbit. For these three situations, the changes in the X−Y plane projection are more noticeable than before. 40 Chapter 3. Numerical tools This time, the orbits not only change in the length of the path, but they also start leaning towards the initial point, as it can be seen in Figure 3.23d, where the 3D view of the torus is presented along the three orbits. (a) 2D projection of the orbit in the ρ−zplane. (b) 2D projection of the orbit in the X−Y plane. (c) Zoomed 2D projection of the orbit in the ρ−z plane. (d) 3D view of the orbit inside the torus. Figure 3.23 Orbit of the proton for the cases KI=0,0.5,1, initial velocity in the OY direction. 3.2 BBNBI As it was previously introduced, BBNBI is a beamlet-based neutral beam ionisation model coupled to ASCOT, designed to satisfy its needs while being a completely independent tool that can be used with other codes. BBNBI follows neutrals from the injector until they are ionised due to the collision with the particles in the plasma, generating a set of fast ions. It is therefore a Monte-Carlo code that models the initial test particle ensemble, i.e. the injection and deposition of neutral beams. This tool was developed by the Aalto University in 2014, and further information can be found in [27]. Coupled to ASCOT, it ensures that the magnetic and plasma inputs used are the same both for the calculation of the neutral beam ionisation and for the fast-ion orbit tracing. This aspect has gained importance with the use of three-dimensional magnetic fields. Moreover, BBNBI can take into account ionisation even outside the last closed surface, which is useful when ASCOT is intended to study fast particle wall loses. In BBNBI, the beam is modelled as an ensemble of sub-beams, that receive the name of beamlets. It follows the particles from the last accelerator grid (see Section 2.3.1), and allows to realistically 3.3 ASCOT 41 Figure 3.24 Front (a) and side (b) view of a NBI grid. The angle formed by each semi-grid and the vertical is the tilt. Source: O. Asunta et al. [27]. take into account the geometry of each beamlet individually by establishing their location and direction. To do so, parameters such as the orientation of the grid, the vertical and horizontal focal lengths, the tilt angle or the divergence of the beam are defined when modelling the NBI system. Figure 3.24 shows the grid as seen by BBNBI, where each of the dots represents a beamlet. The vertical focal length (the vertical distance where the beamlets are focused) and the divergence parameters are also shown in this figure, where it is also possible to appreciate the tilt angle. All the advantages of this approach where beamlets can be defined individually are detailed in [27]. Furthermore, BBNBI allows to define adjustable parameters for each beamlet, such as the particle species that are being injected, the total power of each NBI injector, the main energy of the beam (what was referred to as E0 in Section 2.3.1), or the fraction of particles that will have the different energy components: E0 , E0/2 and E0/3 . These are all parameters that will be included in the preliminary design for the SMART NBI system. To generate an NBI particle, a neutral is assigned a velocity in the direction of the beamlet, and simulated until it enters the vacuum chamber or hits an obstacle, such as the edges of the aperture through which the beam is injected. 3.3 ASCOT ASCOT stands for Accelerated Simulation of Charged particles Orbits in Tori and it is a particlefollowing Monte Carlo code. It has been developed since the early 1990s at Aalto University (formerly HUT, Helsinki University of Technology) and VTT (Technical Research Centre of Finland) [32]. This code follows the orbit of a particle in a tokamak (or stellarator) device, in a domain that ranges from the core plasma to the first wall, and models the plasma-particle interaction, adding physics such as Coulomb collisions by means of Monte Carlo techniques. It allows the 3D study of many complex phenomena by using diverse MC operators [36]. The entire ASCOT code was rewritten in 2013-2014, and the whole philosophy behind it was reviewed (see [45]), resulting in what is known as ASCOT4. In 2019, the code was updated to the version that was used in this thesis, ASCOT5. This was also the process for BBNBI. The upgrades were developed to take full advantage of present and future many-core CPUs [28], and they were benchmarked against their previous versions. For more detailed information about the evolution of the ASCOT code, see [45]. 42 Chapter 3. Numerical tools ASCOT discretizes the set of equations described by the Lorentz Law, making use of Equation (3.44) to obtain an ODE, with x being the position of the particle and v its velocity, and advances them in time. To solve the full gyro motion, the modified leap-frog method [46] is used. In the absence of an electric field, this method conserves the magnitude of the particle velocity, hence conserving the energy. ˙ x=v(3.44) This tool has been used in several studies such as the calculation of fast ion power loads on the ITER and ASDEX Upgrade first walls, or the analysis of the impurity deposition pattern for JET and ASDEX Upgrade [45]. In this thesis, the work with ASCOT will focus on its capability to simulate ions generated by neutral beam injection, as it includes the NBI geometries of JET, ITER, ASDEX Upgrade and MAST-U among others [45]. Moreover, the code provides the distribution of ions that collide with the first wall, including information about their energy and direction. This data is key when designing an NBI system, as the heat loads induced by the beam can severely damage the first wall. 3.3.1 Inputs and outputs To run a simulation, several inputs are needed, as described in the ASCOT guide [47]. Functions are provided in the python library A5PY to generate ASCOT5 inputs from the raw input data. The main ASCOT5 inputs are as follows: •Options . Simulations parameters including simulation mode (gyro-orbit, guiding centre or hybrid) and time-step, end condition, physics involved (collisions or MHD perturbations on/off, etc.), and specifications about the outputs that will be collected, such as the orbits of all the particles or distribution histograms. •Magnetic field . For a simulation to be ran, the magnetic field in the tokamak domain up to the first wall is a mandatory requirement. Both the toroidal field, Bφ , and the poloidal flux, ψ , are required inputs. The poloidal field components of the equilibrium magnetic field ( Br and Bz ) are calculated from the latter. Perturbations of the magnetic field can also be included in the three components. The magnetic field input can be provided in different formats: –Spline interpolated axisymmetric. It requires both ψand Bφon a uniform (R,z) grid. – Spline interpolated non-axisymmetric. It requires ψ on a uniform ( R,z ) grid, and Bφ , Brand Bzfrom coils and external perturbations on a (R,φ,z) grid. –Analytic tokamak. –Stellarator. –Cartesian. •Electric field. It can be spline-interpolated or cartesian. •Plasma profiles. The temperature and density profiles of the plasma species are a required input when running BBNBI, and for simulations where Coulomb collisions have been activated, as they model the interaction between the generated particles and the plasma. It can be lineary interpolated or spline interpolated, and is always 1D (radial, temperature and density functions of ρpol). •First-wall model. As described in the ASCOT wiki: “CAD Model to check collisions between markers and the wall.” This is a fundamental input in order to obtain realistic results 3.3 ASCOT 43 of particle losses. It can be axisymmetric (2D) or 3D via triangular elements. •Markers . “Markers to be simulated.” The generated markers can be either particles, guiding centres or magnetic field line tracers. The first one will be used in this thesis, and requires the particle ( R,φ,z ) initial position and ( vR,vφ,vz ) initial velocities, along with the mass and charge information and a weight factor, which indicates how many real particles each test particle represents. Other parameters such as the density of the neutral species, time-dependent electromagnetic backgrounds or Fourier decomposed MHD perturbations can be included as inputs in ASCOT simulations. They are, however, out of the scope of this thesis. Furthermore, each of these parameters can be automatically generated as a dummy input, for cases where some information is missing or is not relevant. With regards to the outputs , ASCOT provides the user with information about the initial and final state of each marker. It is also possible to store the full orbit of the particle and different 5D and 6D distributions as histograms. 50 Chapter 4. ASCOT simulations for the design of the NBI system: MAST-U and SMART Figure 4.6 Histogram of the z coordinate for the different versions of BBNBI and the two configurations of the MAST-U NBI. (a) Initial state. (b) Final state. Figure 4.7 Energy histograms at the initial and final state of the particles obtained with BBNBI5. Additionally, the energy and pitch measured by the FILD detector are presented in Figure 4.9. The particles hitting FILD with half and third the main injection energy can be explained as the BBNBI5 simulation registers more particles hitting FILD than the BBNBI4 simulation, as previously mentioned. Nonetheless, both distributions show very similar intensities in the main injection energy range. Finally, the shine through obtained with the two codes and shown in Table 4.2 can be compared. Both results are very similar, and they even show a small reduction in the percentage from one version to the other that could be related to the differences shown in the histograms. From all the results presented in this section, it can be concluded that the modelled NBI in BBNBI5 is very similar to the model obtained with the previous version, and hence the design of the SMART NBI system can be carried out with this tool following a similar procedure. 4.2 SMART: Following orbits with ASCOT 51 (a) Final state of all the particles that hit the wall. (b) Initial state of the particles that hit the wall but not in the FILD surroundings. Figure 4.8 Location of the BBNBI5 particles whose simulation ended because they hit the wall. The colorbars represent the number of particles in each cell. (a) BBNBI5. (b) BBNBI4. Figure 4.9 Energy vs Pitch distribution for the particles that hit FILD. 4.2 SMART: Following orbits with ASCOT The next step towards the SMART NBI design is to generate the needed inputs. In contrast to the previous section, where the MAST-U inputs were given, SMART has never been simulated in ASCOT5 before, so the first stage will be to develop those inputs, as detailed in Section 4.2.1, and then verify that they work properly. This is presented in Section 4.2.2 by testing different orbits and performing a sensitivity study, in order to analyse the change in the orbits when the initial parameters such as the pitch angle are modified. 4.2.1 Inputs As previously mentioned, in order to optimise the SMART NBI configuration with ASCOT, the inputs for this device needed to be generated, as it is the first time that SMART is simulated in ASCOT5. 52 Chapter 4. ASCOT simulations for the design of the NBI system: MAST-U and SMART Figure 4.10 Sketch of the SMART wall. As introduced in Section 3.3.1, the main inputs needed are the wall and the magnetic field. For the wall , the configuration and measures presented in Figure 4.10 are used. As for the magnetic field , the ψ map for SMART has been calculated with FIESTA (see [19]). FIESTA is an axysimmetric Grad-Shafranov force balance solver, that calculates the magnetic equilibria in tokamaks, and was the tool used in [19] to calculate the baseline and positive and negative triangularity configurations’ discharge equilibria for SMART Phase 2. The baseline is a reference configuration, and is the one used for the simulations presented in this master thesis. This data included information on the poloidal flux ψ and the toroidal field, Bφ , as it was detailed in Section 3.3.1, along with the perturbations of the magnetic field. In order to use the data obtained from FIESTA, the convention sign had to be changed from COCOS7 (tokamak COordinates COnventionS, see [48]) as used in FIESTA to COCOS3 used in ASCOT. The flux surfaces obtained for this baseline configuration are shown in Figure 4.11, with results in concordance to those presented in [19]. For the plasma input, the density and temperature profiles for both electrons and ions were obtained via ASTRA [49], a solver for transport equations in toroidal geometry. Note that these profiles were calculated without external heating. The data used for the simulations of SMART is shown in Figure 4.12. 4.2.2 Sensitivity study As previously mentioned, the generated inputs need to be tested before proceeding to the NBI design. To pursue this objective, several markers were generated by injecting them at the separatrix ( r=0.6m,z=0m ) and modifying the initial pitch. The energy used in these simulations was set to 15keV, the same energy that will be used in the NBI design, as a scan in energy was performed to 4.2 SMART: Following orbits with ASCOT 53 Figure 4.11 Contour distribution of the poloidal magnetic flux function for SMART. The separatrix is shown in red and the wall and the magnetic axis in black. (a) Density profile. (b) Temperature profile. Figure 4.12 Density and temperature profiles for electrons (e) and ions (i). analyse how the particles injected with different energies were confined. As the energy is increased, so is the velocity, and hence the gyroradius increases (see Equation (2.10)). This sets a limit to the maximum injection energy that can be used, as the size of the vessel is finite. From the analysis performed, 15keV is therefore expected to be the optimum one. 54 Chapter 4. ASCOT simulations for the design of the NBI system: MAST-U and SMART As ASCOT needs the three components of the velocity as inputs, in the first place, the energy can be used to calculate the modulus of the velocity: E=1 2mv =15keV (4.1) where mis the mass of the particle, which will be Hydrogen atoms. Then, the definition of the pitch angle is used to obtain the components of the parallel and perpendicular velocities. As introduced in Section 2.2, the definition of the pitch angle is related to the parallel velocity. In this master thesis, the negative sign will on occasion be included in the definition for convenience. Λ=cosλ=±vk v(4.2) Following this definition, and knowing that vk is parallel to the magnetic field, it can be calculated as: vk=v·cosλBr Ber+Bφ Beφ+Bz Bez(4.3) The perpendicular velocity, however, could point in any direction perpendicular to the magnetic field, so for this first exercise the direction ez×B was chosen, i.e. the gyrophase is defined arbitrarily. This way, the expression for v⊥is: v⊥=v·sinλ  Bφ qB2 φ+B2 r er−Br qB2 φ+B2 r eφ (4.4) As the total velocity vector will be the sum of the parallel and perpendicular velocities, v= vk+v⊥ , the desired ( vR,vφ,vz ) initial velocities are the sum of the components obtained in the previous expressions.   vr vφ vz  =     vcosλBr B+vsinλBφ qB2 φ+B2 r vcosλBφ B−vsinλBr qB2 φ+B2 r vcosλBz B      (4.5) The pitch angle is known, and so are the components of the magnetic field, than can be evaluated using the ASCOT5 library A5PY. In this case, the particle will be initially simulated for 10 µ s, with a fixed time step of 1·10−8s and the simulation mode set to gyro-orbit. To obtain passing orbits, a large pitch should be used, so the parallel velocity is high enough to overcome the magnetic mirror. The first two orbits, obtained for Λ= [0.9,0.5] are presented in Figure 4.13, for both the poloidal and the toroidal projections. It can be appreciated how the gyroradius increases its radius as the pitch decreases. It is therefore possible to perform a comparison between the analytical calculation and the result obtained with the simulation. As previously introduced: rL=v⊥m |Ze|B(4.6) 4.2 SMART: Following orbits with ASCOT 55 where all the parameters are the same for both orbits except for the perpendicular component of the velocity, that will be calculated for the initial state via the given pitch. Hence: rL0.5 rL0.9 =v⊥0.5 v⊥0.9 (4.7) With equations 4.1 and 4.2, and knowing that v=qv2 k+v2 ⊥ , the perpendicular velocity can be obtained, resulting in the analytical solution: rL0.5 rL0.9 =1.99 (4.8) By measuring the gyroradius in Figure 4.13a, it is clear how the radius corresponding to the green orbit, λ=0.5 , doubles that of the purple orbit, matching the result obtained analytically. It is also possible to analyse the drift that the orbits undergo in the high field side, where they are shifted to the right as a result of the drift pointing upwards. Trapped orbits are obtained when setting the initial pitch to lower values. The results obtained for Λ= [0.12,0.10,0.08,0.02] are presented in Figure 4.14. In this case, the toroidal projections may seem unexpected, so a further analysis will be performed with the help of Figures 4.15a and 4.15b. It is interesting to see that, for the same simulation time, the orbits with the higher pitch (higher parallel velocity) travel a longer path in the toroidal projection. Figure 4.15a shows the toroidal projection of the the orbit obtained with Λ=0.12 , in purple, along with a blue trace, that represents the result obtained running the same simulation but using the guiding-centre simulation mode in ASCOT, that considers the magnetic moment derivative to be zero. As introduced in Section 2.2, this variable can be considered constant when the evolution is adiabatic. However, when the magnetic field gradient in a gyro-radius is high enough, this condition is not fulfilled, and thus the guiding-centre formulation is inexact. This is the case for SMART (and also for MAST-U), as it is a spherical tokamak, where the magnetic gradient is proportional to 1 r2 , with a high gradient for small gyro-radius. This explains why in Figure 4.15a the blue trace does not match the expected guiding-centre of the purple orbit, confirming that the guiding-centre approximation is not valid in spherical tokamaks. Figure 4.15b shows the orbit obtained with Λ=0.12 , but ran for a simulation time of 50 µs , in order to be able to analyse how the orbit evolves with time. This way, the pattern shown in Figure 4.14b is repeated along the toroidal projection. From Figures 4.14a and 4.15b, it can be analysed that for these trapped orbits, the inner part of the orbit in the poloidal projection (Figure 4.14a) has a much shorter leg in the toroidal projection (Figure 4.15b). That behaviour is presented in this last figure as several turns grouped together in a short path, as the particle encounters the magnetic mirror and bounces back and forth in a brief amount of time, and continues on the outer side of the poloidal projection up to the next mirror. By checking the q profile of the plasma in Figure 4.16, it is possible to understand this behaviour, as the q parameter increases considerably towards the separatrix. This explains why the inner leg describes a shorter path in the toroidal projection. 56 Chapter 4. ASCOT simulations for the design of the NBI system: MAST-U and SMART (a) Projection in the poloidal plane. (b) Projection in the toroidal plane. Figure 4.13 Sensitivity study of the orbit trajectory to its decreasing initial pitch. Passing orbits obtained for E = 15keV and λ =[0.9,0.5], simulated for 10 µs . Both the wall (in black) and the separatrix (in red) are shown for reference. (a) Projection in the poloidal plane. (b) Projection in the toroidal plane. Figure 4.14 Study of the sensitivity of the orbit trajectory to its decreasing initial pitch. Trapped orbits obtained for E = 15keV and λ=[0.12,0.10,0.08,0.02], simulated for 10 µs. Both the wall (in black) and the separatrix (in red) are shown for reference. 4.2 SMART: Following orbits with ASCOT 57 (a) The orbit presented in Figure 4.14b is shown here along with the result obtained performing the simulation with the guiding-centre approximation. (b) By continuing the simulation up to 50 µs, it is possible to see how the pattern shown in Figure 4.14b is repeated along the full toroidal path. Figure 4.15 The orbit obtained for λ =0.12 can be used to perform a further analysis of the results. Figure 4.16 Evolution of the safety factor qwith the normalised ρpol coordinate. 58 Chapter 4. ASCOT simulations for the design of the NBI system: MAST-U and SMART 4.3 SMART: NBI design In this section, the different configurations studied to optimize the NBI system will be introduced. Its main parameters are presented in Table 4.3. For a detailed description of the elements shown here, see Section 4.1.1. As opposed to MAST-U, SMART will only work with hydrogen in its initial phases. Table 4.3 SMART NBI main parameters. Hydrogen Energy [keV] 15 E fraction [0.62, 0.27, 0.11] Power [MW] 0.3 A 1 Z 1 For the NBI design, several configurations were tested generating sets of 50000 particles. To generate the different configurations, a key parameter was modified: the tangency radius of the beam centre line, as shown in Figure 4.17. The aim of this analysis was to optimise the parameters that were presented in Section 2.3.1: the percentages of prompt-losses and shine-through need to be minimised, while the current drive should be as high as possible. However, as ASCOT5 does not allow a direct calculation yet, an estimation will be used (see [50]): Itot ∼vk v·(1−fst ) , where fst is the shine through fraction. Finally, the slowing down distribution will also be studied. To generate these distributions, the simulations will be run with Coulomb collisions active, allowing the energy exchange. Therefore, the energy histograms of the initial and final state will not be identical as they were in the previous section. The concept of the Alfvén velocity was introduced in Section 2.3.1. For SMART, at ρpol =0.5 , vA=1.42·106 m/s. As the injection energy used is of E0=15 keV, the obtained velocity is 1.70·106 m/s. Therefore, the main energy set is high enough to allow the study of Alfvén Eigenmodes. Figure 4.17 Sketch of the toroidal view of the SMART NBI configuration. The black circles represent the vessel walls and the red ones the toroidal projection of the separatrix. The yellow line represents the injected beam. The results for the different configurations studied by varying this parameter are shown in Table 4.4. As the separatrix projection in the equatorial plane is located at 0.205m and 0.647m, Rt has been modified from 0.252m to 0.515m. The toroidal view of the birth distribution of the particles 4.3 SMART: NBI design 59 generated with these two extreme configurations are shown in Figure 4.18. As the tangency radius is decreased, the beam travels a longer path through the plasma (up to a certain level), which explains the decrease in shine through shown in Table 4.4. (a) Configuration 1, Rt=0.515m. (b) Configuration 6, Rt=0.252m. Figure 4.18 NBI birth distribution for the two extreme configurations. By looking at the results in this table, it is already possible to start optimising the design in this preliminary study. The shine through variation is high in the two first configurations, but stabilizes in the other four, changing only 1% between them. Thus, outlines 1 and 2 could be discarded based on this parameter. Nonetheless, the percentage is still slightly higher than desired. In order to minimise the shine-through, doubling the plasma density has shown to reduce this value by 80-90%. The prompt losses for the third set-up onwards suffer an important increase from one configuration to the next. Hence, it could be established that the last two configurations should be discarded based on this parameter. Finally, the current drive increases continuously as the tangency radius is decreased. However, the difference between the values obtained for configurations 3 and 4 is not substantial. The slowing down distributions for these two outlines are presented in Figure 4.19. These images show how the pitch is maximised in configuration 3 for the higher energies. Therefore, as this is also the configuration with a significant lower percentage of prompt losses, it will be chosen as the optimum one. Table 4.4 Results for the different configurations studied for the SMART NBI system. Configuration ID Tangency radius Shine through Prompt losses Current drive 1 0.515m 20.89% 0.01% 0.49 2 0.461m 16.10% 0.02% 0.54 3 0.408m 13.49% 0.03% 0.59 4 0.355m 12.55% 1.11% 0.62 5 0.303m 11.98% 8.46% 0.66 6 0.252m 12.79% 19.81% 0.69 4.3.1 Resulting preliminary design The main geometrical parameters of the optimum NBI obtained are listed in Table 4.5, and different histograms will be presented in the following. Figure 4.20 shows the birth distribution for the chosen configuration, both in the toroidal and 66 References [15] Sobre la fusión, Accessed: 2021-04-06. [Online]. Available: https:// www.ifmif.org/ ?page_ id=28. [16] El gobierno autoriza la financiación de la candidatura española al ifmif-dones en granada, Accessed: 2021-04-06, 2021. [Online]. Available: https:// www.ciencia.gob.es/ portal/ site/ MICINN/ menuitem.edc7f2029a2be27d7010721001432ea0/ ?vgnextoid=7c051f0235b38710VgnVCM1000001d04140aRCRD& vgnextchannel=cf091f4368aef110VgnVCM1000001034e20aRCRD. [17] Elsewhere: 250 tokamaks throughout the world, Accessed: 2021-04-06. [Online]. Available: https:// www.irsn.fr/ EN/ publications/ thematic-safety/ ITER/ Pages/ 3-250-tokamaks- throughout-the-world.aspx. [18] A. Mancini, J. Ayllon-Guerola, S. Doyle, et al., “Mechanical and electromagnetic design of the vacuum vessel of the smart tokamak,” Fusion Engineering and Design, vol. 171, p. 112 542, 2021, issn: 0920-3796. doi:https:// doi.org/ 10.1016/ j.fusengdes.2021.112542. [Online]. Available: https:// www.sciencedirect.com/ science/ article/ pii/ S0920379621003185. [19] S. J. Doyle, D. Lopez-Aires, A. Mancini, et al., “Magnetic equilibrium design for the smart tokamak,” 2021, Paper Submitted to Elsevier. [20] Nstx upgraded to the world’s most powerful spherical tokamak, Accessed: 2021-04-07. [Online]. Available: https:// fusenet.eu/ node/ 1124. [21] G. R. et al, “Microinstability properties of small-aspect-ratio tokamaks,” vol. 3, no. 1667, 1998. doi:10.1063/ 1.871686. [Online]. Available: https:// aip.scitation.org/ doi/ abs/ 10. 1063/ 1.871686. [22] C. Ribeiro and F. Team, “The concept of the compact ultra-low aspect ratio tokamak culart,” 46thEPS Conference on Plasma Physics, 2019. [Online]. Available: http:// ocs.ciemat.es/ EPS2019PAP/ pdf/ P4.1065.pdf . [23] Spherical tokamaks in comparison to the conventional doughnut-shaped tokamaks, Accessed: 2021-04-07. [Online]. Available: https:// ccfe.ukaea.uk/ . [24] Mast upgrade turns on the switch, Accessed: 2021-04-14. [Online]. Available: https:// www. iter.org/ newsline/ -/3513. [25] Mast upgrade is the uk’s national fusion experiment, Accessed: 2021-04-14. [Online]. Available: https:// ccfe.ukaea.uk/ research/ mast-upgrade/ . [26] Mast upgrade research plan 2019, Accessed: 2021-07-06. [Online]. Available: https:// ccfe. ukaea.uk/ wp-content/ uploads/ 2019/ 12/ MAST-U_RP_2019_v1.pdf . [27] O. Asunta, J. Govenius, R. Budny, et al., “Modelling neutral beams in fusion devices: Beamletbased model for fast particle simulations,” Computer Physics Communications, vol. 188, pp. 33–46, 2015, issn: 0010-4655. doi:https:// doi.org/ 10.1016/ j.cpc.2014.10.024. [Online]. Available: https:// www.sciencedirect.com/ science/ article/ pii/ S0010465514003701. [28] J. Varje, K. Särkimäki, J. Kontula, et al.,High-performance orbit-following code ascot5 for monte carlo simulations in fusion plasmas, 2019. arXiv: 1908.02482 [physics.plasm-ph] . [29] J. Rivero-Rodríguez, “Development and installation of a scintillator based detector for fastionlosses in the mast-u tokamak,” [30] X. Hao, L. Ming-Hai, C. Ming, et al., “Radial magnetic field in magnetic confinement device.,” vol. 24(9), [Online]. Available: http:// cpb.iphy.ac.cn/ article/ 2015/ cpb_24_9_095202.html. [31] M. Kruskal and R. Kulsrud, “Equilibrium of a magnetically confined plasma in a toroid,” [Online]. Available: http : / / www - naweb . iaea . org / napc / physics / 2ndgenconf / data / Proceedings\%201958/ papers\%20Vol31/ Paper27_Vol31.pdf . References 67 [32] O. Asunta, “Fast ions in fusion plasmas – towards numerical tokamaks,” Ph.D. dissertation, Aalto University, Department of Applied Physics, 2015. [33] F. Dini, R. Baghdadi, R. Amrollahi, and S. Khorasani, “An overview of plasma confinement in toroidal systems,” in Sep. 2010, pp. 159–279, isbn: 978-1-61668-352-8. [34] H. Järleblad, L. Stagner, M. Salewski, et al., “Fast-ion orbit sensitivity of neutron emission spectroscopy diagnostics,” Review of Scientific Instruments, vol. 92, no. 4, p. 043 526, 2021. doi:10.1063/ 5.0040696. [Online]. Available: https:// doi.org/ 10.1063/ 5.0040696. [35] R. Kleiber and R. Bilato, IPP Summer University for Plasma Physics. Max-Planck-Institut Für Plasmaphysik, 2011. [36] A. Salmi, “Fast ions and momentum transport in jet tokamak plasmas,” Ph.D. dissertation, VTT Technical Research Centre of Finland, 2012. [37] Plasma heating - ohmic and by neutral beams, Accessed: 2021-04-27. [Online]. Available: https:// www.youtube.com/ watch?v=rfYCrewjCLU. [38] Neutral beam injection and its importance for tokamaks, Accessed: 2021-04-27. [Online]. Available: https:// www.youtube.com/ watch?v=Hc1jcpsp8Zs. [39] J. Wesson and D. J. Campbell, Tokamaks. Oxford University Press, 2004, vol. 118. [40] S. Zweben, “Pitch angle resolved measurements of escaping charged fusion products in tftr,” Nuclear Fusion, vol. 29, Jan. 1989. doi:10.1088/ 0029-5515/ 29/ 5/ 012. [41] J. G. - Q. et al, “Velocity-space sensitivity and tomography of scintillatorbased fast-ion loss detectors,” Plasma Physics and Controlled Fusion, vol. 60, no. 10, 2018, issn: 0741-3335. doi:10. 1088 / 1361- 6587/ aad76e. [Online]. Available: http:// stacks.iop .org / 0741- 3335/ 60/ i=10/ a=105005?key=crossref.b42633c5760bd48362347bf45a33e689. [42] J. F. Rivero-Rodriguez, M. Garcia-Munoz, R. Martin, et al., “A rotary and reciprocating scintillator based fast-ion loss detector for the mast-u tokamak,” Review of Scientific Instruments, vol. 89, no. 10, p. 10I112, 2018. doi:10.1063/ 1.5039311. [43] Ángel Franco García, Campo magnético producido por una corriente circular en un punto fuera de su eje, 2016. [Online]. Available: http:// www.sc.ehu.es/ sbweb/ fisica3/ magnetico/ espira/ espira_1.html#:~:text=secci\%C3\%B3n\%20no\%20nula-,Campo\%20magn\ %C3 \ %A9tico \ %20producido \ %20por \ %20una \ %20corriente \ %20circular, punto \ %20fuera\%20de\%20su\%20eje&text=Donde\%20dl\%20es\%20un\%20elemento,se\ %20calcula\%20el\%20campo\%20magn\%C3\%A9tico.\%EF\%BF\%83\%EF\%BE\ %A7. [44] R. H. Good, “Elliptic integrals, the forgotten functions,” vol. 22, 2001. [45] E. Hirvijoki, O. Asunta, T. Koskela, et al., “Ascot: Solving the kinetic equation of minority particle species in tokamak plasmas,” Computer Physics Communications, vol. 185, no. 4, pp. 1310–1321, 2014, issn: 0010-4655. doi:https:// doi.org/ 10.1016/ j.cpc.2014.01.014. [Online]. Available: https:// www.sciencedirect.com/ science/ article/ pii/ S0010465514000277. [46] J. P. Boris, “Relativistic plasma simulation-optimization of a hybrid code,” Proceeding of Fourth Conference on Numerical Simulations of Plasmas, 1970. [47] J. Varje, Guide to ascot5. [Online]. Available: https:// wiki.aalto.fi/ display/ ASCOTCo/ Guide+to+ASCOT5. [48] Sauter, O., Medvedev, and S. Yu., “Tokamak coordinate conventions: Cocos,” Computer Physics Communications, vol. 184, pp. 293–302, 2013, issn: 0010-4655. doi:http:// dx.doi. org/ 10.1016/ j.cpc.2012.09.010. [Online]. Available: https:// www.sciencedirect.com/ science/ article/ pii/ S0010465512002962?via\%3Dihub. 68 References [49] P. Y. G.V. Pereverzev, ASTRA - Automated System for TRansport Analysis. Max-Planck- Institut Für Plasmaphysik, 2002, http:// w3.pppl.gov/ ~hammett/ work/ 2009/ Astra_ocr.pdf . [50] A. Salmi, M. Gryaznevich, P. Buxton, M. Nightingale, and T. Tala, “Neutral beam optimisation for the spherical tokamak st40,” Fusion Engineering and Design, vol. 117, pp. 14–19, 2017. doi:http:// dx.doi.org/ 10.1016/ j.fusengdes.2017.02.001. Agradecimientos Antes que nada, agradecer al grupo de Plasma Science and Fusion Technology (PSFT) de la Universidad de Sevilla la oportunidad de haber empezado en un campo tan apasionante y necesario como es la fusión. Gracias a Juan Francisco Rivero por la paciencia, las explicaciones con y sin Señor de los Anillos y por haber compartido conmigo tu experiencia aún cuando no nos habíamos conocido, sin ese primer email no estaría aquí ahora. Gracias a Manolo Toscano por esa primera conversación hace casi un año donde me cambiaron los planes vitales por completo en cuestión de minutos, y por tu cariño y guía desde entonces. Gracias a Eleonora Viezzer y a Manolo García Muñoz por la confianza depositada en mi trabajo y el apoyo brindado estos meses. Agradecer a mis compañeros la increíble acogida. Habéis hecho que me sienta como en casa, a veces literalmente. Gracias a Mica, Kiera, los Javis (y Lau), Jose, Pilar, Antonio, Pablo... y una larga lista. Habría sido imposible sobrevivir estos meses sin vuestra ayuda y guía a todos los niveles, dentro y fuera de la oficina. Sin las conversaciones sin sentido en la Alameda, los collages y los deportes extremos. A mis niñas, por la paciencia infinita todos estos años, manteniéndome siempre en el camino de la luz del sol, recordándome que la vida no es solo ingeniería. Gracias también a Marta, Antonio e Inma, este trabajo es un poquito vuestro. Por último, gracias, gracias y mil gracias a mis padres por haber hecho todo lo que estaba en vuestras manos para que yo pudiera llegar donde quisiera, apoyando todas mis decisiones. Y por supuesto, gracias a mi hermano Ángel por ser la persona en la que siempre puedo contar, la que siempre me va a sacar una sonrisa y la que siempre me va a decir la verdad. Y a mi otro hermano en la práctica, Lucas. Saber que pase lo que pase vas a estar ahí es la mayor de las suertes. 69