Thermal analysis of the FILD head in the ASDEX Upgrade Tokamak

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 Fin de Grado Grado en Ingeniería en Tecnologías Industriales - Intensificacíon energética Thermal analysis of the FILD head in the ASDEX Upgrade Tokamak Autor: Javier Hidalgo Salaverri Tutor: Manuel Toscano Jiménez Tutor: Juan Manuel Ayllón Guerola Dpto. Física Aplicada III Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Trabajo Fin de Grado Grado en Ingeniería en Tecnologías Industriales - Intensificacíon energética Thermal analysis of the FILD head in the ASDEX Upgrade Tokamak Autor: Javier Hidalgo Salaverri Tutor: Manuel Toscano Jiménez Profesor Titular de Universidad Tutor: Juan Manuel Ayllón Guerola Profesor Invitado Dpto. Física Aplicada III Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2018 Trabajo Fin de Grado: Thermal analysis of the FILD head in the ASDEX Upgrade Tokamak Autor: Javier Hidalgo Salaverri Tutor: Manuel Toscano Jiménez Tutor: Juan Manuel Ayllón Guerola 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: Agradecimientos Am i hermano y a mis padres, por escucharme hablar de ingeniería durante cuatro años aunque no me entiendan (y los que quedan) y por su apoyo incondicional. A María, que ha hecho más llevadera la carrera de fondo del grado. A mis tutores: Manuel Toscano, por creer en mi y apoyarme durante estos cuatro años y a Juan Manuel Ayllón por recibirme con los brazos abiertos en el mundo de la fusión. Saber es noche I VIII List of Figures 3.8 End results for FILD 4 with 2 cm insertion. Using a mesh with 8958 faces versus the 747 faces of the auxiliary mesh. Yellow represents wet triangles and blue dry ones 26 3.9 Schematic of a limiter shadowing a FILD. A similar situation to this figure happens to FILD 4 27 4.1 a) Magnetic flux expansion near the divertor b) Geometric flux expansion [22] 31 4.2 Basic overview of the simulation software. In bold are scripts are in red, end results 33 4.3 Inner operation of field_ line_ tracing.m script. Used to get the trajectories born in each triangle and discern whether is dry or wet. Scripts in bold and results in red 33 5.1 a) Mesh used for ANSYS, more detailed to get better end results. b) Mesh used for MatLab, rougher to lower computational time 35 5.2 Correction to MatLab data done using a rigid transformation throught ANSYS’s External data component system 36 5.3 ANSYS mesh triangulated with a) relevance = 25; b) relevance = 100 36 5.4 Validation results via reverse mapping for a) relevance = 25; b) relevance = 100 37 5.5 For FILD 4 with 3 cm of insertion: a) wet (yellow) or dry (blue) mesh; b) heat flux distribution for totally parallel convective flux and no shadowing; c) heat flux distribution accounting for shadowing and intersection angle of the colliding particles 38 5.6 Cyclic operation of FILD 4 in intervals of one second reading at 3 cm insertion and one second cooling at parking position. FILD heat map for instants 3 s, 8 s and 17 s are also represented. Figures with the FILD heat map at times A, B and C are shown in figure 5.7 39 5.7 FILD 4 temperature profile during a 20 second cycle for times: a) 3 seconds with Tmax = 237.48ºC; b) 8 seconds with Tmax = 253.12ºC; c) 17 seconds with Tmax = 424.00ºC 39 5.8 FILD 4 Tmax evolution during a 20 seconds exposition and temperature profile at the end of the simulation 40 List of Tables 1.1 ASDEX Upgrade technical data [11] 10 3.1 Comparison between using or not the OcTree algorithm in a medium refined mesh at three different positions for all FILDs. ASDEX Upgrade shot #34570 was used 23 3.2 Comparison between using or not the frontier algorithm in addition to OcTree in a highly refined mesh at three different positions for all FILDs. ASDEX Upgrade shot 34570 was used. 26 IX 1 Introduction The aim of this work is to develop a flexible tool to study the thermal behavior of fast ions loss detectors under a wide range of circustances. Before entering inside the study of FILD itself an introduction to nuclear fusion, plasma physics, the tokamak reactor and fast ions is needed. 1.1 Nuclear Fusion Nuclear fusion in a general sense is the collision of two light elements to birth a heavier one liberating a great amount of energy. This requires that the binding energy per nucleon (the energy the nucleus has in order to maintain its inner stability) of the resulting nucleus would be bigger than the lighter ones. This is easily seen in figure 1.1 where the maximum is around the iron (with A = 56), for lighter elements this energy would be generated by fusing two light elements meanwhile for heavier ones, energy would be produced by the fission of a heavy element. Figure 1.1 Binding energy per nucleon vs. mass number for the dominant form of each chemical element. [1]. Different reactions can be used in order to achieve fusion. As the figure 1.1 suggests the best solution would be to choose the lightest element in order to get the maximum energy 1 2Chapter 1. Introduction from the reaction. Therefore the focus will be in hydrogen and helium, specifically in the reactions: D-D, D-T and D-He 3 (being D, deuterium and T, tritium) for being those the optimal reactions. In figure 1.2 it is shown the cross section of each of those reactions in function of the kinetic energy of D, in other words, of its temperature. As the figure shows the D-T reaction has the biggest cross section and the lower requirement in terms of temperature. That is the reason behind the election of this reaction as the first to be studied in fusion reactors and, therefore, the one that would be taken into account for now on. Figure 1.2 Experimentally measured cross sections for the D–T, D–He3, and D–D fusion reactions as a function of deuterium kinetic energy (KD). [2]. The D-T reaction occurs as: D+T→α(3.5MeV)+n(14.1MeV)(1.1) This reaction produces around 3.52 MeV/nucleon which translates to the macroscopic world to 338x10 6 MJ/kg. For comparison, the energy released in a U 235 fission reaction accounts for 0.88 MeV/nucleon and 84x10 6 MJ/kg. But for the fusion to happen it requires around 100 keV (around 100 million K) this means the hydrogen gas will be ionized and therefore it would have become a plasma. For fusion to take place without any external help (this is, a self sustainable reaction) the Lawson criterion must be accomplished. Lawson criterion (also known as the triple product) gives the conditions were a plasma can exist as a function of its temperature (T), plasma pressure (p) and the confinement time ( τe ) in a simplified way this inequation can be written as [1]: p∗τe≥0.11 ∗T2 <σ∗v>(1.2) With <σ*v> the average value of the cross-section times the relative velocity. 1.2 Plasma 3 This criterion opens different ways of approaching fusion plasmas. Depending on how the parameters are chosen the main options are: •Gravitational Fusion: the only fusion source in nature. The problem is addressed with high pressures. Pressures that can only be obtained by the tremendous gravitational forces that are found in the nucleus of a star. Therefore this way of achieving fusion is not in the reach of technology. •Inertial Confinement Fusion: D-T fuel in the form of a 2 mm sphere is highly compressed and then hit by numerous lasers in a brief period reaching temperature of keV. Though the confinement time in this case is very small (in the order of µ s), high pressure and specially high energies allow the fusion to happen. [3] •Magnetic Confinement Fusion: D-T fuel is in the form of plasma. It is confined for few seconds thanks to a series of magnetic fields. Temperatures of around 100 kev are achieved thanks to Ohmnic heating and neutral beams between others. To close the triple product the pressure is supplied by the pressure of the own plasma and the one applied by the electromagnetic field. This is the kind of fusion that takes place in the ASDEX Upgrade Tokamak and therefore the one that will be the subject of this work. 1.2 Plasma Plasma is commonly known for being the forth state of matter. It occurs when a gas is heated to high enough temperature that the atom overcomes its ionization threshold taking place the separation between the ionized nuclei and their electrons. This results to a charged cloud of ions and electrons. For an hydrogen atom to ionize an energy of 13.6 eV [ 4 ] must be supplied and as it was stated in the previous section temperatures of 100 keV are needed for the fusion reaction D-T to succeed, so it is assure the mix will be ionized. The fusion therefore occurs when kinetic energy overcomes Coulomb’s barrier and strong nuclear forces come into play fusing both nuclei. Plasma is generally characterized by: •Ionized gas: nuclei and electrons are no longer linked due to the high kinetic energies overcoming the Coulomb attractive forces. •Quasi-neutrality: in a microscopic scale plasma is charged due to the separation between charges caused by the ionization but, if seen macroscopically, it appears to be neutral. This phenomenon occurs because each ion is surrounded by a cloud of electrons with the same absolute charge which gives birth to neutrally charged spheres. These spheres, called Debye’s spheres, when translated into the macroscopic world result into an apparently non-charged plasma. •Collective behavior: plasma being an ionized gas is governed by Maxwell’s laws, so the electric potential should promote the rapprochement between ions and electrons 4Chapter 1. Introduction due to Coulomb’s force. However this is not the case in plasma. As it has already been said before, plasma temperatures are quite high and therefore its kinetic energy will be as well. Comparing kinetic and potential energy means that the potential is neglected against the huge value of kinetic energy, so the effect of Coulomb’s force is not enough to join ions and electrons. The only moment when the potential could overcome the high kinetic energy would be when an electron randomly got close enough to the nucleus making the potential even bigger than kinetic (after all potential energy grows as 1/r) but this situation won’t happen often due to the low densities in plasma. Given all of this, the electric density can be represented as a smooth, continuous distribution of electrons and ions and therefore long range interactions along the plasma will occur. [1] •High electrical conductivity: low collisions and high kinetic energy means low resistivity and therefore high conductivity. As the plasma temperature is increased conductivity increases, this is a problem when talking about heating up the plasma as one of the means to do it is by Ohmic heating using the plasma as the secondary of a transformer. The electromagnetic properties of fusion plasma all together with its high temperatures makes working with it an important challenge both for engineering and for physics. 1.3 Particle motion in plasma Particles in plasma are affected by an electric and a magnetic field. Suppose first the particle is only affected by an uniform magnetic field, therefore it will have a circular motion in a perpendicular plane to those magnetic field lines as shown in figure 1.3 centered around the called guiding center. This movement is characterized by its ciclotronic frequency ( ωc ) and its Larmor radius (rL) or gyroradius: ωc=|q|∗B m(1.3) rL=m∗v⊥ |q|∗B=v⊥ ωc (1.4) Particle motion is also commonly represented using the pitch angle, defined as the ratio between the parallel component of the velocity and the total velocity. Λ=v|| vt =cos(α)(1.5) An electric field perpendicular to the magnetic one is added, this means that the particle accelerates so Larmor radius is modified. Due to the presence of both fields a drift happens 1.4 Magnetic confinement 5 Figure 1.3 Circular motion of a particle in a field line [5]. denoted as ExBdrift, shown for ions and electrons in figure 1.4. ~vd=~ Ex~ B B2(1.6) Figure 1.4 ExB drift in ions (left) and electrons (right) [6]. Therefore the real motion will be the sum of the two of them: an helicoidal movement with a guiding center traveling at v d speed, ρ gyroradius and ωc ciclotronic frequency [ 5 ]. This description of particle motion in plasma is also known as full orbit. 1.4 Magnetic confinement Currently magnetic confinement is the most advanced way of fusion. Even though no Tokamak has been able to produce positive net energy yet, in Cadarache (France) ITER project is being built as the culmination of a world wide project lead by the European Union 6Chapter 1. Introduction with the goal of producing 30 times more energy than the fed to the system. This landmark would demonstrate the technological and physical feasibility of fusion energy as a way of producing commercial energy. Magnetic confinement takes advantage of the ionized nature of plasma to confine it using a series of powerful magnets. Using this concept two main designs have been raised: Tokamak and Stellarator, both of them toroidal reactors. But before entering in the differences between both designs, why toroidal? Figure 1.5 Comparison between linear and toroidal confinement in a Tokamak.[7] . Figure 1.5 shows the difference between linear and toroidal confinement. At first sight a linear fusion device could seem like the best option as the toroidal geometry brings to the table harder challenges for engineers and physicists, indeed this was historically the first design in magnetic confinement. But the problem with this design is that it has some intrinsic losses that are not acceptable. As the device is finite, at the end of it particles will collide with the walls, this would be disastrous not only for the plasma stability but for the walls themselves that would be severely damaged in a few shots. Also in a linear configuration, between two consecutive magnets some of the field lines won’t be closed and therefore particles that would follow them won’t be confined neither. These reasons make the toroidal geometry the best option, as it has no ends where the particles could collide (at least if the confinement were ideal) and inside it magnetic field lines are perfectly closed. But toroidal geometry has its own flaws too. With this configuration the distance between two magnets is not constant anymore. In the inner radius magnets will be closer than in the exterior part and therefore as B α 1/r, the magnetic field would experience a gradient within that would make the device unsuitable for magnetic confinement. In how this problem is approached is where the difference between Tokamak and Stellarator resides: •Tokamak: Russian design developed in 1956, it is shown in the right part of figure 1.5. The magnets configuration by themselves will produce the magnetic gradient issue but this is fixed by a central solenoid that induces a current in the plasma canceling the magnetic drift. The saturation current of the solenoid has become one of the limiting factors to get higher confinement times in this kind of reactors. Thanks 1.5 Tokamak and ASDEX Upgrade 7 Figure 1.6 Stellarator conceptual design [7]. to its relatively easy geometry, Tokamaks presents axisymmetry which makes its analysis quite easier. •Stellarator: Created by the USA in 1951, represented in figure 1.6. In this case the magnetic drift is solved by the disposition of the own magnets, for this to be accomplish numerical simulations must be run in order to obtain the very complex shapes that each magnet must have. Even though this solves the magnetic drift problem in a clean way (from a purely physics point of view), some other issues are created along the way. The most obvious one is a manufacturing problem that involves serious engineering challenges; the second one passes on the effectivity of the own reaction, the complex geometry of the magnets create little free volume for the plasma and therefore for a given reactor size, Stellarator will have less plasma than a Tokamak. Therefore net positive energetic reactions are harder to obtain with this design. Another problem comes from the changing form of the plasma, which means that there is no axisymmetry and therefore each section will have to be studied separately. The benefits and problems enumerated for both designs have tip the scale in favor of the Tokamak design that currently is the most advanced solution for nuclear fusion. And therefore will be the design that will be used for ITER as the first try to get to an almost commercial state. 1.5 Tokamak and ASDEX Upgrade In this section Tokamak characteristics will be defined but with a special focus on the peculiarities of ASDEX Upgrade Tokamak. As it was explained in the previous section Tokamaks achieve magnetic stability by a series of magnets and a central electrical conductor. Those magnets are positioned in two ways: poloidal magnets and toroidal magnets. In figure 1.7 both directions are defined and also the poloidal magnets are shown, normally toroidal magnets are not depicted in schematic figures as it tend to make figures harder to understand. Now that toridal and poloidal directions are defined, a description of the general parts of 14 Chapter 2. The Fast Ion Loss Detector (FILD) • Good time resolution to follow the evolution and frequency of losses due to MHD instabilities. •Particle flux that impacts the vessel. • Detector must be flexible and resistant to stand a hostile environment with high neutron flux, gamma rays and high temperatures. FILD is a scintillator based detector. This means particles, fast ions in this case, collide with a sensible sheet, called scintillator, that registers the impact. Concretely in FILD, fast ions go through a pinhole (figure 2.2) where a collimator is set to prevent the entry of thermal ions. Inside the device, the collision takes place with the scintillator who registers the location of the collision, this data will be necessary to reconstruct the trajectory of the fleeing ion. The scintillator is formed by a sheet of stainless steel that serves as the base for a covering of phosphorus. Figure 2.2 Fast ion entering through FILD pinhole [14]. Once fast ions collide with the scintillator, light is emited. This light is collected by an optic system which divides it in two equal beams. One is then received by a high spatial resolution camera (CCD camera) and the other by a high temporal resolution one (APD camera). This leads to a good definition of the escaping ions [ 15 ]. This layout is represented by figure 2.3. Inside ASDEX Upgrade a total of five FILDs are installed, following the distribution of figure 2.4. FILD 1 and 2 are aligned are at the same poloidal position but in different toroidal locations. FILD 3 and 5 are close to the upper and lower divertor respectively. FILD 4 is able to perform fast movements thanks to a Magnetically Driven Reciprocation System (MDRS). FILD 5 is situated between two magnetic coils which facilitates the study of the effect of those coils on the fast ions [ 16 ]. FILD measuring positions are a compromise solution between getting good readings and not receiving high heat loads. 15 Figure 2.3 Basic layout of the optical collecting system of FILD in MAST-U Tokamak. Light emitted by the scintillator is divided and then received by a CCD camera and a APD camera. [15]. Figure 2.4 Poloidal distribution of FILDs in ASDEX Upgrade [16]. In section 1.4 axisymmetry was presented as one of the pros of the Tokamak design. Meanwhile this is true as a first approximation, if the calculation is more accurate the nouniformities in the magnetic field produced by discrete magnets (also discussed in section 1.4) disrupt the assumption. Therefore having FILD 1 and 2 in the same toroidal line allows studying its differences which is a measurement of the strength of the axisymmetryc hypothesis. And with a poloidal plane covered thanks to the rest of FILD a good enough idea of the 3D behavior of fast ions along the Tokamak can be obtained. FILD head is manufactured in graphite and coated with tungsten, a heat resistant material 16 Chapter 2. The Fast Ion Loss Detector (FILD) but not enough for plasma edge temperatures and the ion flux. This is why the head is not always exposed to plasma. Normally the device is in a parking position behind the first wall where it still receives thermal loads but in a bearable levels. When it is time to get some readings FILD head is introduced in the vacuum vessel where it starts collecting fast ions and with them high loads are deposited over the head surface. This situation can’t be sustained too long without structural damages. Therefore a good thermomechanical simulation of FILD during its working time must be achieved in order to get FILD working optimally and without damages, this is the motivation of this study. 3 Modeling the particle collisions with FILD In this section the motion of fast ions and electrons will be simulated together with their collisions with FILD head. To get whether an intersection happens or not, Möller- Trumbore ray-triangle algorithm will be used which will require a previous refinement of the given FILD mesh to get well-defined results. This section corresponds with the core of this work, creating a series of tools that permits this simulations to be flexible and efficient. MatLab software will be used for this chapter. 3.1 Mesh refinement A 3D model of ASDEX Upgrade vacuum vessel and of FILD head has been provided by CNA (Centro Nacional de Aceleradores) of University of Seville. Figure 3.1 FILD head mesh: a) initially b) after triangulation with triangle surface minor than 10 cm2. Head model was given as shown in figure 3.1 a). Its geometry is quite simple, a big 17 18 Chapter 3. Modeling the particle collisions with FILD cylinder with some flat surfaces, so the mesh can be defined by triangles of a considerable size. This rough geometry would give inaccurate results on which zones of the FILD head receive collisions and therefore heat flux. To solve this, the initial triangles must be subdivided into smaller ones until all triangles get small enough to achieve the desired precision. For this matter a MatLab program was created. A maximum area is given to the program as an input, then all triangles are checked if their surface is bigger than the reference. If this happens those triangles are subdivided into four equal-area smaller ones, this continues in a loop until no triangle is larger than the reference, at the end with a reference of 10 mm 2 the result is shown in figure 3.1 b). To make this subdivision code from MathWorks repository was used, concretely subdivide_ tri.m [17]. 3.2 Field Line Tracer Not the whole exposed surface of FILD head receives impact of particles. Part of those particles collide first with other elements of the reactor. This phenomenon is called shadowing [ 18 ]. Taking into account that FILD surface is triangulated gives birth to two different types of triangles: plasma wet triangles, which receive collisions of particles and dry triangles, which receive no collisions. This creates zones with a higher heat flux than the bordering ones with the consequent thermal stress. In section 1.3 full orbit motion was described. Applying this in MatLab suppose a high load in terms of computational resources. For this tool, particle motion will be simulated as if they were following the magnetic field lines on the plasma. This reduces code complexity and computational time and doesn’t harm the final goal, that is developing a tool to obtain the heat flux map over FILD in numerous situations, not a precise simulation of plasma particles. Therefore a magnetic field line tracer must be developed in order to register which triangles are wet or dry. The most straight forward approach to this problem is initiate n particles in random positions in plasma and after some simulation time see how many have collide or not. This obviously is not a really good idea, as the size of n should be big enough to assure that every possible trajectory is followed. So the contrary approach results as an easier path. Particles will start in each triangle (defined by its center of masses) and will go ”back in time” until it collides with the reactor walls or other detector or a given time passes uneventfully. If time passes and no collision is registered then that trajectory would be feasible and therefore that triangle would be wet. Otherwise, if a collision takes place, that trajectory could not take place so the triangle would be dry. For this collisions, particles going through the positive and negative direction of the magnetic field will be registered. Particle trajectories are also subjected to the magnet and plasma configuration in the moment of the simulation. This magnetic field data has been obtained from ASDEX Upgrade inner database of plasma shots, where for a given shot and time a magnetic field 3.3 Triangle-Ray Intersection 19 3D mesh is provided. Having the magnetic field mesh the particle motion simulation can start. The simulation starts in a node of FILD head, in given x y z position at time t 0 . As the particle will follow magnetic field lines, the objective is to define the magnetic field in its position and use that vector as its velocity to calculate its position in the next instant. Integrating this way results in a complete trajectory in the integration interval. First, particle parameters in the given position are calculated: α=tan−1(y/x)(3.1) Radio =px2+y2(3.2) With this parameters magnetic field vector in particle position can be obtained by interpolating the 3D mesh, for this, Radio is used. The result is in cylindrical coordinates and to convert it to Cartesian, a rotation matrix is used:   Bx By Bz  =  Bradial Btoroidal Bz  ∗  cos(α)sin(α)0 −sin(α)cos(α)0 0 0 1  (3.3) With the magnetic field vector defined in the particle position, integration can start and the whole trajectory recorded. Those trajectories are expected to be closed ones as those particles must be confined within the reactor. Once the collision check is run it will be seen which ones lose this confinement due to this collisions with the 3D structures of the reactor. 3.3 Triangle-Ray Intersection Once all trajectories for a given FILD at a given measuring position are recorded, it is time to check if any collision occurs before the integration time is due. This is done by checking every segment of a particle trajectory with all the triangles that form the vacuum vessel and the five FILDs. A highly vectorized MatLab code of Möller and Trumbore [ 19 ] intersection algorithm is used. This MatLab version was obtained from the MathWorks repository, the function is called ”TriangleRayIntersection.m” [20]. The algorithm uses an unitary triangle perpendicular to the x=0 plane and expressed in barycentric coordinates. In barycentric coordinates, points inside a triangle are written as T = f(u,v,w) with each of those coordinates being equal to 1 in a vertex a 0 in the other two. This creates a contour condition inside the triangle (eq. 3.4) that converts the triangle from 20 Chapter 3. Modeling the particle collisions with FILD 3D to 2D. An equilateral triangle in barycentrial coordinates is represented in figure 3.2. u+v+w=1(3.4) Figure 3.2 Equilateral triangle expressed in barycentrical coordiantes [21]. In this coordinate system a point within the triangle is defined by: T(u,v) = (1−u−v)∗V0+u∗V1+v∗V2 u≥0v≥0eq. 3.4 →u+v≤1(3.5) Being V0,1,2 the vertices of the triangle. A ray is then define as equation 3.6, passing by the origin and with normalized direction D. R(t) = 0+t∗D(3.6) The intersection point (if there were one) will be the equality between equations 3.5 and 3.6: 0+t∗D= (1−u−v)∗V0+u∗V1+v∗V2→−D V1−V0V2−V0∗  t u v  =0−V0→ →M∗  t u v  =0−V0 (3.7) So the coordinates from the intersection point (u,v) as well as the normalized distance 3.4 Optimizing the code 21 from the origin (t) are both the results of linear equations with the light computational load this means. For any triangle in a 3D environment it will have to be first converted. This is done by translating it to the origin and then by multiplying it by the inverse of the transformation matrix M, from equation 3.7. This process is summarized in figure 3.3. Figure 3.3 Steps to follow from a 3D triangle to a Möller-Trumbore-ready triangle [19]. For further information about the algorithm or its general application refer to the original paper [19]. This algorithm is applied to every segment of each particle trajectory. During the checks if one segment registers an intersection with any triangle of the environment, the function for that trajectory stops and the origin of that trajectory is recorded as a wet one and so on the triangle it represents. With this algorithm together with the field line tracer the triangles that form FILD head are already defined as wet or dry. For the wet ones its position and angle of incidence is recorded as they will be necessary in the upcoming sections. Being this angle of incidence, the angle between the field line and the normal to the surface. 3.4 Optimizing the code In some scenarios where FILD head is highly exposed, simulation time can takes even days for a well refined mesh. Getting down this time, saving computational resources, is essential in order to make the code more user friendly and useful. In this section the optimization methods used are presented. 3.4.1 OcTree Vessel render (figure 3.4 alone accounts for 162000 triangles and still it would have to be added up the amount of faces in each of the five FILDs. This means that each segment of 22 Chapter 3. Modeling the particle collisions with FILD all trajectories would have to check for collision over 162000 times, this is computationally unacceptable. OcTree algorithm has been used to cut down this computational load to more acceptable levels. This algorithm has already been used in other studies with similar profiles as [18]. Figure 3.4 Vessel mesh used for MatLab simulations. Peaks represent ports for detector or other devices. OcTree intends to decrease total computational time by reducing the amount of triangles each segment has to be compared with. The algorithm divides the to-be-compared-with mesh into two submeshes (called bins) with the same amount of points, this process is repeated with each new smaller mesh until a certain criteria is found (by default, the amount of points inside a bin is equal or smaller than a tenth of the total amount). Now the intersection would only be checked with points inside the same bin as the segment, reducing drastically the computational time. The problem is that the triangle-ray intersection does not check with points but with triangles. So as the triangle is not completely defined it won’t be checked with the segment in question. This idea is illustrated in figure 3.5. To solve this, an index relating each vertex with the face it belongs has been created. This index has the same length as the amount of points in the mesh, in each row the faces where this vertex is participating are stored. This means that for the points contained in target bin, this index returns all faces related to said bin. These faces will be then checked with the ray-triangle intersector. Also the situation were the initial point of the segment and the final one are not in the same bin must be covered. In this case it may happen that segment passes by a third bin. This could be solved by giving the field line a very small step but this will result in an increasing of the computational time. A faster solution is to not use the OcTree algorithm in this particular case were initial and final points of the segment are in different bins. 3.4 Optimizing the code 23 Figure 3.5 The green colored zone of the bin won’t be checked if only full triangles inside the bin are taken into account. The efficiency of the algorithm depends on the configuration of the simulation. Generally a configuration that would increase the computational time, will also increase the effect of OcTree. For instance, if FILD is parked all the trajectories will collide almost at the beginning of the simulation, so the algorithm won’t be called as many times as if the head was totally exposed and, due to this, most of its triangles were wet. So there is no way to give an accurate number of the impact of this algorithm in the overall computational time. However, a demonstration case has been prepared to illustrate the benefits of using the aforementioned algorithm. For this example a medium refined mesh was used with a total of 747 faces. Simulations were run for each FILD in three positions: at parking position and at 4, 8 cm of insertion. Results are shown in table 3.1. Plasma shot used was ASDEX Upgrade #34570. Table 3.1 Comparison between using or not the OcTree algorithm in a medium refined mesh at three different positions for all FILDs. ASDEX Upgrade shot #34570 was used. FILD 1 FILD 2 FILD 3 Insertion [m] 0 0.04 0.08 0 0.04 0.08 0 0.04 0.08 OcTree [s] 623.33 1690.4 2240.2 855.61 745.93 633.85 1477.7 1070.9 1350.3 no OcTree [s] 739.48 2339.6 3009.9 811.53 820.3 731.28 2052.3 1217.1 1739.4 Savings [%] 15.707 27.748 25.572 -5.432 9.066 13.323 27.998 12.012 22.370 FILD 4 FILD 5 Insertion [m] 0 0.04 0.08 0 0.04 0.08 OcTree [s] 916.17 1368.4 1589 1593.8 1727.7 1983.4 no OcTree [s] 1162.3 1850.1 2249.3 1949.2 1859.1 2194.1 Savings [%] 21.176 26.036 29.356 18.233 7.068 9.603 Even though as it was stated before no real percentage of saved time can be defined, a similar value can be obtained with the mean of the results as those simulations cover a wide range of situations. Therefore applying the OcTree algorithm results in a mean saved 30 Chapter 4. Thermal Model for plasma facing components Stangeby [ 26 ] relates the perpendicular and parallel influx of particles as a way of checking the availability of this model for certain situations: φ⊥ φ|| =qLf/L(4.3) Being Lf, the length front face and L = π*RFILD φ⊥ φ|| =p0.1/(π∗1.6)→φ|| φ⊥(4.4) Now that the feasibility of cosine model for this study has been proven, equation 4.1 can be dissected: • q 0 : upstream peak heat flux, is the heat in the separatrix. Sieglin [ 22 ] and Mitteau [23] gives it a value of around 100 MW/m2. •δ: distance of the target from the separatrix. δ=s−s0. •λq: e-folding length of the heat flux. Experimentally the power that escapes convectively the reactor to impact with the first wall is known, via power balance [ 22 ]. Knowing this power and the total surface, q 0 is obtained by imposing that this power is equal to the received by all wet surfaces. The cosine model is a experimental model based in the knowledge that the convective parallel heat is the predominant and that is evolves as an exponential. The λq parameter is not well known yet and therefore a range of parameters has to be studied in order to find one that adjusts to reality. Different studies [ 18 ], [ 23 ],[ 27 ] situates the range of values of λqfrom 1 mm to 1 cm. Equation 4.1 is defined in the outside mid plane (omp) as the results varies poloidally. This change can be input directly by calculating the heat flux as a function of ρpoloidal (this concept will be defined in next section) or by adding an new parameter f x . The first case will be the used for this study. Even though heat flux will be calculated as a function of ρpoloidal , f x parameter is explained below as the physical idea of the difference between omp and other planes is better explained using this parameter. As it was previously stated the convective model follows the magnetic field lines, therefore if the density of field lines decrease the rate of variation of the heat flux will be slower. This is what happens when moving polloidally to the divertors, magnetic field lines will open gradually as represented in figure 4.1 a). This is accounted by the 4.2 Model application 31 magnetic flux expansion (fx, mag). fx,mag =Rmp ∗Bmp pol Rx∗Bx pol (4.5) Figure 4.1 a) Magnetic flux expansion near the divertor b) Geometric flux expansion [ 22 ]. Also the same effect can take place with the incidence angle between field lines and surfaces. This is illustrated in figure 4.1 b). This variation is represented with f x, geo . This geometric effect can also be account separately by multiplying q 0 by the cosine of the angle between the normal surface and the magnetic field line. This idea was already presented in equation 4.2. fx,geo =1 sin(αpol)(4.6) Knowing that fx=fx,mag ∗fx,geo , the flux expansion is integrated in the cosine model is represented in equation: q|| =q0∗e−δ λq∗fx(4.7) 4.2 Model application This model will be applied using the MatLab software. As it was stated in the previous section, the magnetic flux expansion won’t be explicitly represented. The magnetic effect 32 Chapter 4. Thermal Model for plasma facing components will be taken into account by transforming equation 4.1 from coordinate δ to ρpoloidal . Meanwhile the geometric part will be calculated with the pitch angle between field line and wet surface. As this model requires the impact between magnetic field lines and FILD faces, only wet triangles will receive it. First the poloidal normalized radius ( ρpol ) must be defined. This normalized radius has its zero value at the magnetic axis of the plasma and is equal to the unit at the separatrix. As the separatrix real position varies in the z-axis, ρpol = 1 will change its Cartesian position. This variation is equivalent to the shown by the magnetic flux expansion. For this task kk library from ASDEX Upgrade has been used. kkMap function transforms a point from Cartesian coordinates to poloidal ones for a given plasma shot. With this function the position of the separatrix (s 0 ) can be obtained by imposing ρpol = 1. Therefore equation 4.1 is already defined in the omp. To get the rest of points the normalized position of all FILDs must be obtained. Again the kkMap function is used, this time with all wet points of the FILD heads. Power deposition over FILD will come from the interpolation between the heat flux value at ρpol, omp at ρpol,FILD. Finally the geometric flux expansion is applied by multiplying the previous result by the cosine of the pitch angle that was already calculated in section 3.3. By this point the heat flux [W/m 2 ] in each wet triangle is already known, so by multiplying each one by its own surface the total power received by the FILD head at that position is known. 4.3 Numerical tool description Now that the whole simulation is finished and results ready to use in ANSYS has been obtained, an overview of the inner operation of the software is given as a user-friendly fast way of getting into its usage. First inputs are defined, the only output is the heat flux through each of the triangles that form FILD: •FILD mesh. •Refinement for the auxiliary and target FILD meshes. •Vacuum vessel mesh. •Magnetic field mesh. •FILD position and target FILD for the simulation. Figure 4.2 represents the basics of the operation. First, the FILD mesh is refined into 4.3 Numerical tool description 33 two different steps, an auxiliary mesh to find the frontier between dry and wet triangles and a well refined one that will be the real target of the simulation. Those two meshes together with the vessel mesh are carried into field_ line_ tracing.m where a traject structure is obtained. Within this structure information of each trajectory is kept as well as if each triangle is dry or wet. traject is run in wet_ triangles.m where the intersection angle ( α ) is calculated for each wet triangle. traject and α are finally used in q_ map.m to get the resulting heat flux in each triangle. So until q_ map.m, the program would be related to chapter 3 meanwhile the rest would linked with this very chapter. Figure 4.2 Basic overview of the simulation software. In bold are scripts are in red, end results. As field_ line_ tracing.m represents the major part of the software, a more detailed version of its inner operations is shown in figure 4.3. Figure 4.3 Inner operation of field_ line_ tracing.m script. Used to get the trajectories born in each triangle and discern whether is dry or wet. Scripts in bold and results in red. 34 Chapter 4. Thermal Model for plasma facing components First the OcTree algorithm is applied with the vessel mesh, giving birth to the OT class which contains all the information related with the OcTree bins. This class in conjunction with the auxiliary FILD mesh are used in FLIntersection.m, script that checks for collisions, to obtain the limits of the frontier. Again the same script is used but this time to analyze the target mesh between the limits of the frontier. The result is the traject structure. 5 Results In this chapter one specific FILD position is studied in ANSYS software in order to obtain the temperature profile in different scenarios. The aim of this chapter is not to discern whether the FILD will be damaged or not in that situation but to show the potential uses of the results from the numerical tool developed in this work. The studied case is FILD 4 with an insertion of 3 cm. From MatLab a TXT file is returned, composed by the node ID, XYZ coordinates and heat flux in W/m 2 . This data is fed to ANSYS Transient Thermal analysis system through the External Data component system. FILD mesh used for ANSYS is not the same as the used for MatLab, which was a rougher one, a comparison between both meshes can be found in figure 5.1. Figure 5.1 a) Mesh used for ANSYS, more detailed to get better end results. b) Mesh used for MatLab, rougher to lower computational time. ANSYS mesh is more detailed in order to get more accurate results, on the other hand MatLab mesh looks for efficiency and therefore has little details. The highlights in their differences are the pinhole, that is treated as a plane in MatLab and the end part were ANSYS mesh is shorter and more detailed, also in general corners are not as smooth in MatLab as in ANSYS. To transfer the data from MatLab to ANSYS faithfully, surplus 35 36 Chapter 5. Results points have been cut out as well as the points in the detailed bottom zone; this is easily justifiable as those points have never gone wet during all done simulations. Pinhole points have not been touched as their influence in the final results are meaningless. Also data from MatLab is not in the same position as in ANSYS this has been corrected via the Rigid Transformation field of External Data. Numbers used in the correction are shown in figure 5.2 Figure 5.2 Correction to MatLab data done using a rigid transformation throught ANSYS’s External data component system. To check if there is a good correlation between MatLab and ANSYS mesh reverse mapping utility from ANSYS has been used. From the given mesh ANSYS develops its own triangulation where the software will work by interpolating the heats of the MatLab mesh were located in the barycenter of the triangles composing that mesh. Refinement of ANSYS mesh can be controlled with the relevance parameter, in figure 5.3 FILD mesh is shown with relevance of 25 and 100. The similarities between both ANSYS meshes and MatLab one are then checked with reverse mapping, which checks the differences between source and mapped data [28], results are shown in figure 5.4. Figure 5.3 ANSYS mesh triangulated with a) relevance = 25; b) relevance = 100. Results, as expected, are better with a higher relevance being the most conflictive zone around the rounded corner of the front of FILD. In future versions, if needed, a more refined mesh of that conflictive surface can be prepared for MatLab to assure better results 37 Figure 5.4 Validation results via reverse mapping for a) relevance = 25; b) relevance = 100. there. Nonetheless that area is not critical for this scenario as the highest thermal flux is located in a well mapped surface as will be shown later. For both situations a simulation with a 20 seconds long exposition has been run, getting the relevance-25 mesh a maximum temperature of 730.74ºC and relevance-100, 732.87ºC; a 0.29% difference that shows the convergence of the mesh. Even though the difference is quite small and relevance-25 mesh could be used, relevance-100 will be the chosen one due to it having slightly better results with no significant time increase. In figure 5.5 it is shown separately the wet/dry mesh (5.5 a)); the heat map without taking into account the intersection angle α , this is the whole heat flux being totally parallel (5.5 b)) and, finally, the addition of both situations to give birth to the heat distribution for FILD 4 at 3 cm insertion (5.5 c)). Two scenarios have been prepared for FILD 4. A 20 seconds cyclic simulation and an insertion during another 20 seconds. For both simulations an initial temperature of 100 º C has been chosen as it is common for FILD to stay hot due to a prior operation of the detector and the lack of active cooling. Also a contour condition has been set as 100 º C for the back of FILD in order to dissipate the incoming heat flux. First scenario represents the normal operation of a FILD. During a plasma shot FILD goes out and into parking position several times to control its temperature. This simulation has 1 second periods of no thermal load which mean the parking position and then another second of full heat load which is the inserted period. It could be discussed whether the duration of the simulation is unrealistic for current plasma shots but in the future ITER this duration may be a reality, also the aim of this section is not to give a realistic simulation of a FILD operation (i.e. there is no transition between parking, 0 cm, and 3 cm insertion) but to present the possibilities this work opens in matter of thermal simulations of FILD. Results for the first scenario are shown graphically in figure 5.6 where the maximum temperature at the FILD head is represented as a function of time. As it was expected the cyclic and discrete nature of the simulation is translated into sawtoothed graph. The 38 Chapter 5. Results Figure 5.5 For FILD 4 with 3 cm of insertion: a) wet (yellow) or dry (blue) mesh; b) heat flux distribution for totally parallel convective flux and no shadowing; c) heat flux distribution accounting for shadowing and intersection angle of the colliding particles. three intermediate states represented, help showing the irregular heating of FILD that was already advanced in figure 5.5. 39 Figure 5.6 Cyclic operation of FILD 4 in intervals of one second reading at 3 cm insertion and one second cooling at parking position. FILD heat map for instants 3 s, 8 s and 17 s are also represented. Figures with the FILD heat map at times A, B and C are shown in figure 5.7. Figure 5.7 FILD 4 temperature profile during a 20 second cycle for times: a) 3 seconds with Tmax = 237.48 º C; b) 8 seconds with Tmax = 253.12 º C; c) 17 seconds with Tmax = 424.00ºC. Second scenario (represented in figure 5.8) is a simpler one, just a 20 second long