## CHAPTER B - B1 Provision of Support to the Advancement of the ITER Physics Basis

**1.2 Energy and particle confinement/transport**

1.2.1. Turbulence and transport phenomena

Background and Objectives: The presence of turbulent fluctuations along with coherent structures and wave-fields is studied for their effect on transport phenomena. The main aim of this activity is to provide improved and more accurate evolution equations for particle distribution functions in comparison to standard approaches. The latter commonly involve simplifying statistical assumptions for the underlying particle dynamics which, in general, do not hold for realistic cases occurring in ITER. The analytical tools utilised in this activity are powerful Hamiltonian methods allowing for the direct connection between the phase space structure of the particle dynamics and the evolution of their distribution function describing collective particle behaviour.

Work performed in year 2010 (*in co-operation with PSFC-MIT, USA*):

- Analytical results for the variations of particle velocity and position and their dependence on the characteristics of the localised RF coherent structures (frequency, wave-number, group velocity and spatial width) have been obtained for the case of a single particle interaction with a single localised RF structure (see
**Annex 1**). These results have been utilised for the study of collective particle behaviour under interaction with a statistically distributed set of such field structures by iterating the analytical mapping relating particle position and momentum before and after the interaction with the localised structure as well as by studying the evolution of the particle distribution function in terms of a Fokker-Planck type of equation. - The effect of the intrinsic complexity of particle motion on transport has been modelled in terms of an appropriate evolution equation for the particle distribution function. The proper way to introduce the effect of stochastic fields has been identified and the interplay between the effects of the intrinsic complexity and field stochasticity has been studied in terms of a Fokker-Planck equation in the action space with special attention to the passing-trapped particle boundary.

1.2.2. Development of the Dispersed Phase Structure Dimensionality (DPSD) concept for MHD Turbulence with particle transport

Background and Objectives: Coherent turbulence structures play an important role in the dynamics of turbulence, modify the transport of (passive) scalars and contribute to the preferential concentration of transported particles. In the case of edge turbulence in fusion plasmas, it is important to be able to model the effect of coherent structures on density and electric potential fields. The turbulence structure dimensionality tensor is a one-point tensor that carries non-local information providing a statistical description of coherent structures. This activity is aimed at further developing the dispersed phase structure dimensionality tensor so that it can become applicable to edge turbulence.

Work performed in year 2010:

- In the first part of 2010, a new parallel 3D MHD shear code with scalar transport and Lagrangian particle tracking was completed and a series of small-scale simulations were carried out that verified the applicability of the DPSD (Dispersed Phase Structure Dimensionality) concept in MHD turbulence at low magnetic Reynolds numbers. A number of numerical implementation issues were identified that were successfully resolved during the second half of 2010.
- A series of large-scale direct numerical simulations of particle dispersion in MHD turbulence subjected to mean shear at moderately high magnetic Reynolds numbers (Re_m ~ 20) was initiated that will be completed in the first part of 2011. Based on the results obtained so far (see
**Annex 2**), these simulations demonstrate the ability of the DPSD tensor to provide a concise description of the structure of coherent fields under the influence of magnetic field. In addition, in the second half of 2010, we initiated the effort to test the DPSD concept using edge turbulence data. Instantaneous profiles of density and electric potential fluctuations were generated and the solutions were then interpolated on to a simple toroidal grid (*r*,*θ*,*ψ*). Turbulence structures of density and electric potential fluctuations were calculated on this grid using a newly developed code. Our preliminary results show promise that the DPSD method could be used for benchmarking and diagnostics in various turbulence codes that are currently being developed for edge turbulence.

1.2.4. Numerical Tools for a comprehensive ITER modelling

Background and Objectives: Unstructured Triangular mesh generation approaches in a Tokamak vessel configuration are very promising because they can be adapted easily to geometrical discontinuities, such as the ones that exist in the divertor’s region. Furthermore, they can be extended relatively easily for three-dimensional numerical models. A major requirement that has to be fulfilled is that the triangular mesh has to be aligned with the magnetic surfaces in the vessel. This requirement is very restrictive and satisfying it in real computation is challenging. The objective of this task is to formulate an approach for generating unstructured triangular meshes that satisfy the alignment requirement as best as possible.

Work performed in year 2010:

- In the first part of 2010, we focused in the anisotropic mesh generation techniques. An anisotropic mesh is commonly generated as a quasi-uniform mesh in the metric determined by a tensor specifying the shape, size, and orientation of elements. The metric field is defined either analytically or by the Hessian of a solution. In our work, we have used an open source code (Ani2D) that performs the anisotropic mesh generation. By giving the Hessian of the solution, Ani2D constructs an anisotropic mesh that is aligned as much as possible with the solution. After several tests we realised that the alignment of the triangular meshes with magnetic surfaces was poor. The anisotropic mesh generator could not perform well in real vessel geometries. So our plans were modified accordingly. In the later part of 2010, a second approach that performs a classical Delauny triangulation of the mesh was adopted. This method supplies additional information to the mesh generator in order for the mesh to be as aligned as possible. The strategy we have adopted is to extract information, in the form of isolines, from the equilibrium data of a 2D Tokamak simulation and then demand from a classical mesh generator to construct a mesh by respecting the distances between the isolines. This approach gave well-aligned results and, in the context of IMP12, we have initiated the process of translating the code from Python to C, as a first step towards integrating the code in Gateway with an appropriate CPO interface.

1.2.5. Development of multi-physics computational models for plasma simulations

Background and Objectives: This activity is aimed at investigating tools for the better characterisation of turbulence transport in the edge of Tokamak configurations. The main objective is to couple the SOLPS code for the characterisation of the edge plasma with the Dispersed Phase Structure Dimensionality (DPSD) concept, which provides a statistical description of coherent structures in turbulent fields. In particular, the turbulence structure dimensionality tensor is a one-point tensor that carries non-local information providing a statistical description of coherent structures.

Work performed in year 2010:

- In the first part of 2010, the code SOLPS for edge plasma modelling was analysed in order to identify its limits concerning the turbulence transport as well as the way the integration of the DPSD concept in the SOLPS could be achieved. The implementation of the DPSD tensors in SOLPS was attempted in the second half of 2010, but was only partially achieved due to the complexity of the code and the lack of previous experience in working with the code. In the mean time, parallel work carried out under TRA and ITM (see 1.2.2) has progressed well and in the future is expected to supersede this task.
- Due to pending issues regarding the integration of the DPSD concept in the SOLPS code, the assessment of the integration of high temperature plasma modelling methodologies in existing MHD codes was not performed. Future work in this direction is being superseded by work initiated in 2010 under 1.2.2.

1.2.6. Particle transport model

Background and Objectives: The analysis of realistic problems within the tokamak require the complete solution of the charged and neutral gas particle conservation equations, together with the electric and magnetic field calculations such as minimum assumptions are made. In the past, we had coupled the charged particle mass continuity equations, together with the Navier-Stokes and the Poisson equations, whereas the conservation of momentum and energy for the charged particles were not included. In order to calculate the velocities of charged particles in the past in atmospheric plasmas, the transport properties, which are functions of E/N, i.e. the ratio of the electric field per neutral gas density per unit volume, were used. To simulate collisional plasmas efficiently in the tokamak reactor, one option is to calculate the charged particle transport properties of the plasma fluid model by including the conservation equations for momentum and energy for the charged particles as well.

Work performed in year 2010:

- In the first part of 2010, an in-house plasma code was transformed into a modular Finite Element Flux Corrected Transport (FE-FCT) format (using programming language C++). The transformation of the existing plasma code into a modular FE-FCT method was deemed necessary in order to simplify the implementation of the complicated conservation equations that are needed. A configuration file, that dispatches the number of conservation equations to be solved, the different terms and the way the terms are coupled, both at half-time and at full-time prediction steps of the two-step Lax-Wendroff Finite Element technique, was used. In the later part of 2010, the implementation of the FE-FCT module using the two-step Lax-Wendroff technique has been competed and the ability to develop complicated modular conservation models is now fully in place. Testing of the model using square wave tests in all three directions to validate the accuracy of the code is nearing completion.

**1.3 MHD stability and plasma control**

Background and Objectives: We focus on the improvement of the theoretical understanding and the interpretation of experiments on MHD stabilisation by ECRH/ECCD. Research is conducted in three parts: In the first part, the main target is to include the effect of magnetic island topology on EC wave absorption and current drive in the beam-tracing code TORBEAM. The inclusion is performed on the basis of computing the deposited power density in terms of the volumes of the perturbed flux surfaces, whereas usually the computation is done in terms of axisymmetric flux surfaces. Such an upgrade in TORBEAM will aid in the theoretical prediction of the minimum amount of current needed for NTM stabilisation, taking into account the established criteria for j_{ECCD}/j_{bs} and the novel ECCD computation. In the second part, we develop a self-consistent model of NTM dynamics and stabilisation with ECCD. This is achieved by coupling different wave propagation solvers with different models for the NTM dynamical evolution that account for the non-axisymmetric island geometry. Finally, the third part is devoted to the investigation of the deleterious effect of edge turbulence in ECRH/ECCD applications. In the experiments for ECCD-based NTM control in ITER, the beam has to travel through the plasma edge where turbulence can be rather intense, so that the associated density gradient may scatter the wave. One may expect this effect to be minor, however, since the distance to the EC resonance is long, even a small scattering angle could lead to considerable beam broadening and result in degradation of the stabilising effect. Our study includes both ray and full-wave propagation in plasma with different types of spatial distributions for the density fluctuations.

Work performed in year 2010 (*in co-operation with the institutes indicated*):

- The numerical routine for the computation of the volume between adjacent magnetic flux surfaces was benchmarked successfully for the simple case of circular magnetic equilibrium using the ray-tracing code CODERAY. The specific routine has been constructed with the scope of improving the accuracy in the estimation of the EC power deposition and current drive when magnetic islands are present. Following the benchmark, the routine was incorporated in the TORBEAM code and a detailed investigation of the effect of the non-axisymmetric magnetic island geometry on the localisation of the ECRH/ECCD deposition has been initiated. First results show that the EC current density is computed 2-4 times larger than in the axisymmetric case, which shows that the magnetic island geometry, through the formation of small-volume flux surfaces, enhances the localisation of the ECCD around the island’s O-point. (
*In co-operation with IPP-Garching and FOM*) - The construction of a semi-closed numerical model for the description of NTM stabilisation dynamics with ECCD has been completed. The model is based on the coupling of an existing ray-tracing code for the description of EC propagation (CODERAY) with a Fortran routine solving the modified Rutherford equation for the quantification of the NTM dynamical evolution. The ECCD current density along propagation is computed by CODERAY, then it is provided as one of the input parameters to the modified Rutherford equation solver. Using this tool, we have traced the effect of the perturbed magnetic geometry to the NTM stabilisation, arising from the modification of the ECCD current density in comparison to the computation based on the unperturbed flux surface volumes (for more details see
**Annex 3**). - The modification of the numerical scheme for the boundary conditions of the FDTD code FWTOR has been successfully completed (this issue is analysed in 1.5.1(i)), and the full-wave computation of EC propagation in the presence of edge density fluctuations has been initiated. In the corresponding simulations, the plasma response is described in terms of the linear hot-plasma dielectric tensor in axisymmetric magnetic equilibrium, whereas the distribution of the edge density blobs along the wave propagation path is random, i.e. their geometrical centre, characteristic length and density amplitude are uniformly random numbers. (
*In co-operation with Uni. Bayreuth*)

1.3.2. Development of a two-fluid MHD solver for a toroidal geometry with a general cross-section

Background and Objectives: This activity focuses on using computational MHD in order to understand the physics of plasma instabilities, which are one of the major obstacles in achieving a feasible fusion reaction. The new computational fluid dynamics MHD model will be based on non-spectral methods for the simulation of resistive plasma instabilities and micro-turbulence in tokamaks. Three-dimensional validation of the model for the sawtooth instability in a small tokamak is being attempted. The model will be expanded to include various features of a tokamak, for example the divertor, and advanced numerical techniques, such as large-eddy simulation models, for the microturbulence.

Work performed in year 2010 (*in co-operation with CEA and ULB*):

- Various modules of the new MHD code based on non-spectral methods for the simulation of plasma stability and turbulence in tokamaks were completed. These concern mainly the modification of the open source library called “openFOAM” (a general numerical package for the solution of partial differential equations) in order to be adapted to the ITER geometry and plasma flow conditions. In particular, the basic compressible resistive non-linear MHD equations were implemented into “openFOAM” and the new code was successfully tested in simple MHD flows against existing analytical solutions, such as Orzang-Tang vortex flow. Various numerical issues have been addressed, such as assurance of divergence-free magnetic fields, construction of a 3D mesh for the tokamak torus, global mesh refinement tools, and possible adaptive local refinement tools to handle large gradients. The benchmark solution of the 3D sawtooth cycle in the CDX-U tokamak with the present MHD model and the comparison with existing non-linear results from the M3D and NIMROD codes has been obtained (see
**Annex 4**). We have also focused on the possibility of data and object input/output and the connection of the present MHD model with the KEPLER interface within the IMP12 activities of ITM.

1.3.3 Explore new techniques for the suppression of tearing modes by localised Electron Cyclotron Current Drive

Background and Objectives: The suppression of neoclassical tearing modes (NTMs) by localised electron cyclotron current drive (ECCD) is a well-known technique, which by now has been successfully applied in every high-beta tokamak. The previously conducted theoretical evaluation of the suppression rate of magnetic islands, however, was focused on the so called “helical effect” of the driven current, while its classical effect (due to gradient of the driven current) was overlooked. The objective of the present theoretical study, which actually initiated in the framework of the 2009 program, has been to evaluate this effect.

Work performed in year 2010 (*in co-operation with IPP-Garching*):

- The improved stability of the tearing modes, which was observed in JT60-U by early application of ECCD, is attributed to the classical effect of the driven current which was evaluated using the precise configuration of the current density inside and outside the magnetic island. It was found that the classical effect of ECCD (due to the gradient of the driven current) is stronger by one or two orders of magnitude than the “helical” effect (due to the helical deformation of the driven current) when EC absorption is at the X and O points of a locked island, and negligible compared to the helical effect when EC absorption is between the X and O points (for more details, see
**Annex 5**). It was also found that the total effect of the driven current (i.e. the classical and the helical effect, which is finally what matters at the experiment) on locked islands with the EC absorption at the O point is stronger by one order of magnitude than the total effect on rotating islands. The achievement of improved stability of tearing modes by early application of ECCD in ASDEX Upgrade was postponed for the 2011 campaign because of the priority of the theoretical study.

**1.4 Power and particle exhaust, plasma-wall interaction**

1.4.1. Modelling and large scale MD simulations of chemical erosion mechanisms

Background and Objectives: In order to assess the erosion rate of beryllium walls in future nuclear fusion reactors one needs sputtering yields, which can be conveniently obtained with numerical molecular dynamics (MD) simulations. In 2010, our objective was to perform a massive series of MD simulations of beryllium sputtering due to low-energy beryllium and helium bombardment.

Work performed in year 2010:

In early 2010, we have developed a library of scripts to perform classical MD simulations with the LAMMPS code (Large-scale Atomic/Molecular Massively Parallel Simulator). This code was installed and optimised on a modern computer cluster with 216 cores and 512 Gb random access memory. This allows us to simulate atomic systems composed of hundreds of thousands of atoms and for time ranges nearing nanoseconds. Then, we performed a series of MD experiments to study Be self-sputtering and Be sputtering by He bombardment with the classic pairwise potential (Ueda, Ohsaka, Kuwajima, J.Nucl.Mat, 258-263, (1998), 713-718, 283-287, (2000), 1100-1104) and a modern three-body Tersoff potential (Bjoerkas et al, J. Phys.: Cond. Matt, 21 (2009) 445002), see **Annex 6**. In each individual experiment, we calculated the sputtering yield by analysing anywhere from 512 to 2048 MD trajectories. By carrying out a series of different experiments (each at fixed incident energy and polar angle) we were able to vary the incident atom energy in the range 20-300 eV and the incident polar angle in the range 0-90 degrees. The total number of performed MD experiments was a few million and were completed in late 2010. (For comparison, the previous MD results by Ueda et al. contained altogether a few dozens of thousands of MD trajectories). The key findings from these numerical experiments are the following:

*Be-Be system.*To simulate Be self-sputtering both the newer Tersoff and earlier pairwise potential were applied. The pairwise potential gives the higher assessment for the sputtering yield in contrast to the Tersoff potential. Because the former potential applied for MD simulations is approximately an order of magnitude faster when compared to the newer Tersoff potential, the pairwise potential can be safely used for a fast assessment of the upper limit of Be self-sputtering.*He-Be system.*It was found that previous results by Ueda et al are in satisfactory agreement with our results. The largest difference occurs for small and moderate polar incident angle, where the earlier yields are underestimated compared to ours. We believe this is due to poor statistics (~ 10^{2}MD samples compared to ~10^{3 }in our simulation) taken by Ueda et al. for averaging the sputtering yield.

**1.5 Physics of plasma heating and current drive**

1.5.1. ECRH-ECCD and transport in fusion plasmas

Background and Objectives: In this activity we investigate the Electron Cyclotron (EC) wave propagation, absorption and current drive as well as the resonant interaction of the EC wave with the charged particles of the plasma, in wave and plasma geometry relevant to heating and current drive applications, mm-wave diagnostics and MHD stabilisation. Our research is conducted in several parts:

The first part focuses on the importance of wave effects (like e.g. beam localisation, diffraction, interference and mode conversion) and plasma effects (e.g. hot plasma dispersion, nonlinear absorption and non-local response) in the propagation and absorption of EC waves in tokamak plasmas. Such effects could possibly play a leading role in many cases of interest, like e.g. in incomplete cyclotron-resonant absorption, high-power ECRH/ECCD and EM wave reflectometry/interferometry, where the conventional formalism involving asymptotic wave methods (ray/beam tracing) and local linear dielectric plasma response breaks down and a more sophisticated treatment is called for.

In the second part, the main objective is to ascertain the advantages of the utilisation of non-Gaussian wave beams in ECRH systems. This is formulated in terms of modelling the EC wave propagation in simplified tokamak geometry, using the paraxial WKB beam tracing method for following the wave propagation and the local linear hot-plasma dielectric tensor for the electron plasma response. The results are exploited in order to upgrade certain features of the existing numerical code TORBEAM, regarding the propagation, absorption and current drive of EC beams with arbitrary-shaped electric field profile in tokamak magnetic configuration.

In the third part, nonlinear dynamics of wave-particle interactions in fusion plasmas are considered from the point of view of applications to ECRH and its role in the stabilisation of the NTM, ECCD and plasma diagnostics through interactions with rf pulses. The Hamiltonian formalism along with a set of accompanying tools such as the phase space analysis, the symplectic integration, the Canonical Perturbation Method and Lie transforms are used in order to extend the theory beyond the quasilinear approximation for realistic fusion plasma configurations and general form of the wave spectra.

Finally, in the fourth part the deflection of RF waves off density fluctuation at the plasma edge (blobs) especially in the EC and LH range of frequencies is considered on the basis of idealized blob characteristics (spheres or ellipsoids elongated along the magnetic field lines) and for mild blob density increment (with respect of the ambient plasma density).

Work performed in year 2010 (*in co-operation with the Institutes indicated*):

- In this period, the implementation of a CPML layer in the FDTD code FWTOR for the outgoing-wave boundary conditions was completed (see
**Annex 7**). The replacement of the boundary condition scheme previously used in the code was shown to be necessary in order to improve the numerical accuracy in the case of 2D/3D hot-plasma propagation, in terms of reducing the artificial electromagnetic field generated by the reflection of the propagating fields from the problem space boundary. Furthermore, in order to introduce flexibility in treating different kinds of wave propagation problems, a convenient Cartesian coordinate system was formulated as a local coordinate system for the computational space. All the above completed the upgrade of the FWTOR code, initiated during the previous period, which was benchmarked successfully for simple cases of wave propagation. Following that, we initiated the simulation of the first two problems of interest for ECRH/ECCD in ITER: (a) The effect of edge density fluctuations on the stabilisation of NTMs by ECCD [see 1.3.1(iii)], (b) The effect of hot plasma dispersion on the EC wave propagation. - The numerical computation of the beam electric field amplitude along propagation, in terms of the amplitude transport equation, was incorporated in the routine NGBT and its coupling to TORBEAM was finalised. Furthermore, the corresponding benchmarking of this updated version of the NGBT routine in slab and tokamak geometry was successfully completed. Following these actions, the numerical investigation of the possible benefits on the power deposition and current drive profiles when utilising non-Gaussian beams in the design of the ITER EC wave launchers and mm-wave diagnostics was initiated. First results show that the deleterious effect of higher-order “error” modes in a Gaussian beam might in some cases become important, but also that a more localised power deposition profile may be achieved by a convenient superposition of Gaussian and higher-order beam modes. (
*In co-operation with IPP-Garching*). - The development of the non-singular, time dependent diffusion tensor to a Fokker-Planck solver has been initiated. A one-dimensional code has been implemented and tested in several cases. Furthermore, the development has been continued and the code has been extended to two dimensions (see
**Annex 8**). Several cases have been tested and alternative implementations are still under investigation. - The kinetic theory has been extended in order to include the effect of RF waves on phase space transport of alpha particles. An evolution equation for the alpha particle distribution function has been obtained with the utilisation of the Hamiltonian techniques. In addition, ensemble averaged variations of particle momenta and guiding centre positions have been obtained for the case of alpha particle interactions with localised waves (see
**Annex 9**). - A Fokker-Planck solver based on a non-singular, time dependent diffusion tensor has been implemented for the one-dimensional and the two-dimensional case [see task 1.5.1 (iii)]. The aim is to apply this code along with the full wave FDTD code FWTOR [completed in the first part of year 2010, see task 1.5.1 (i)], among others, in studying the off-axis current drive and rotation. Preliminary studies on issues regarding the coupling of the two codes, such as interfacing in terms of the selection of common parameters, have been initiated.
- Analytical results for the interaction of a particle with a localised wave-packet in the case of a single transit have been utilised in order to construct a mapping which can be iterated for the study of multiple interactions. The analytical mapping has been used for the study of the case of multiple transits of particles (electrons, protons, deuterons and alphas) through localised wave packets and the effect on heating, current drive, and transport (for more details see
**Annex 1**). Moreover, averaged variations of particle momenta and guiding centre positions have been obtained analytically for ensembles of particles under a single interaction with the wave-packet (see**Annex 9**). - An extensive literature survey on the refraction of RF waves generated in a homogeneous gyrotropic plasma by a single gyrotropic sphere of higher density than the ambient plasma density has been completed. Extension to the case where mild density gradients of the ambient plasma and the confining magnetic field is in progress.
- The refraction of RF waves (in the range of EC and LH range of frequencies) on a single blob was extended to the case of statistically distributed in space and size blobs on the basis of a Fokker-Planck approach, assuming mildly higher densities in the spherical blobs than the ambient density. On this basis, the ray expected values of the ray deflection and parallel wave number changes were estimated. The problem was also extended to the case of ellipsoidal blobs elongated along the magnetic field lines (see
**Annex 10**. A further extension with the incorporation of ambient plasma density and magnetic field gradients is now in progress. Additionally, the investigation was also extended to the direction of incorporating the diffraction which has been left out so far through the aforementioned geometrical optics approach. Along these lines, the formulation of the Mie scattering by gyrotropic spherical blobs (a fundamental problem of electromagnetics not fully studied in the past) is near completion. - The interaction of a single ion in a uniform magnetic field with a wave-packet whose envelope is in resonance with the unperturbed motion of the ion has been studied both analytically and numerically. The conditions of preferable energisation of the particle are determined in the case of two beating waves (for more details see
**Annex 11**). The aim of this study is to investigate the possibility of utilising two EC waves (which can efficiently access the core) with a frequency difference in the range of local IC or LH resonance and estimate the momentum transfer and the efficiency of the scheme as far as the current drive capabilities are concerned.

**1.7 Theory and modelling for ITER**

1.7.1 Development of computational fluid dynamics solvers for liquid-metal flows relevant

to blanket modules (DEMO incl.)

Background and Objectives: The reactor blanket and the circulating liquid metal contained within are exposed to the strong magnetic fields used to confine the plasma. When the liquid metal enters or exits the area of influence of the magnetic field, the forces acting upon the liquid metal change, which in turn can modify the characteristics of the flow. The goal of this activity is to accurately predict the liquid-metal flow in poorly conducting pipes when leaving strong magnetic fields and to analyse the effects of the wall-conductivity in such phenomena. The modelling is done using a parallel, unstructured, Navier-Stokes solver developed in UCY.

Work performed in year 2010 (*in co-operation with the institutes indicated*):

- Further developments on the FORTRAN 90/95 software code have been implemented to combine different strategies concerning the numerical discretisation of the Lorentz force in order to optimise accuracy and numerical stability for high Ha number simulations. In 2010, Hartmann numbers flows up to 10000 have been successfully addressed.
- Furthermore, full Direct Numerical Simulations of the turbulence modification of liquid-metal flow entering low, intermediate and strong magnetic fields have been carried out for insulating pipes. Several inflow numerical strategies have been implemented to overcome development lengths issues. A very detailed description of the turbulence decay and of additional flow instabilities observed is provided in
**Annex 12**. Additional code modifications have been carried out to model the effect of conducting boundaries on the turbulence modification of liquid-metal flow entering magnetic fields. Preliminary simulations at low Hartmann number for conducting pipes have been performed. Further details are given in**Annex 12**. - Finally, a feasibility study concerning a wall-function treatment for the resolution of the well-known Hartmann and parallel layers (also known as Shercliff, Roberts layers) of MHD flows has been carried out. The study concludes that defining implicit boundary conditions for velocity and electric potential, thus imposing the mass flow and the electrical currents in near wall-cells, is feasible and greatly advantageous from a computational point of view. The implementation of the wall-functions will clearly enhance the capability and numerical stability of the in-house parallel-unstructured Navier-Stokes solver developed in UCY to address very high Hartmann numbers (Ha>10000) relevant for fusion applications.
- A numerical code was developed in order to simulate the 3D time dependent MHD flow in rectangular ducts with thin finite conductive walls in the presence of strong external magnetic fields. At first approximation, the numerical treatment of the Hartmann and side walls of the duct was based on various wall-layer approximations for the velocity and the electric potential. In a second approach, the finite conductive thin walls were also treated by using an immersed boundary technique. In this approach, the boundary conditions (BC) for the wall-normal electric current (Newmann BC) and the electric potential (Dirichlet BC) are asymmetrically satisfied at the interface. Alternatively, a mixed type boundary condition is used to treat numerically the interface between the electrically conducting fluid and the thin conducting fluid. Several indicative numerical simulations were also performed in order to verify the aforementioned numerical approaches against existing numerical, theoretical or experimental results for the laminar hydrodynamic and MHD flows in square or rectangular ducts with a) four non-conducting walls, b) two non-conducting and two perfectly conducting walls, c) two non-conducting and two thin finite conducting walls, and d) four thin finite conducting walls. Finally, a parametric numerical study covering representative values of Reynolds and Hartmann numbers, wall conductivity and aspect ratio of the rectangular duct has been initiated in the end of year 2010. (
*In co-operation with ULB*.) - There was no activity in this subtask due to a tied work schedule focusing to other subtasks. It will be reinitiated next year. (
*In co-operation with ULB.*) - There was no activity in this subtask because the principal investigator (Dr. M. Vlachomitrou) left the project.

1.7.2 Implementation of an Edge plasma model by solving the Braginskii’s fluid equations for collisional plasmas

Background and Objectives: Edge plasma modelling is very crucial for a better understanding of the plasma wall interaction as well as of the continuous mode of a Tokamak. The physics of the edge region are described by the Braginskii’s equations that are applied at the interface between the hot plasma core and the walls. In that region the plasma is assumed ‘cold’ and the form of these equations is similar to the cold plasma differential equations with the addition of a magnetic field component. Our aim was to apply an existing very efficient cold plasma finite element code in the edge region. This code was optimised for hyperbolic equations with efficiency on discontinuity capturing.

Work performed in 2010:

- From the beginning of 2010 we have been working on the adaptation of an existing 2D cold plasma code to a Tokamak configuration. The modifications that had to be done were of high complexity and especially regarding the coordinates system transformations as well as at the insertion of the viscosity tensor terms that take into account the magnetic field component. In order to obtain a clearer picture on edge plasma modelling we decided to get involved with a mature code that is already functional on edge plasma area (SOLPS - David Coster). For that we contacted the developers of the code and we participated in the code camp of Garching for a better understanding of the SOLPS capabilities. Our involvement started gradually and a first step was to perform few benchmark simulations. In the code camp of Garching, after several contacts with IMP3 members, we also decided to use SOLPS for verification and validation of edge plasma modelling, a task that is in priority in the ITM framework.

At the same time, and during the second part of 2010, a 3D version of the cold plasma code has been under development. The development of fully three-dimentional codes that simulate the whole scrape off layer of a Tokamak vessel could overcome the restriction of using fully aligned meshes in edge codes and could allow a better and more realistic representation of the geometric details of the device. The most restrictive parameter in developing three-dimensional codes in a Tokamak device is the enormous computational load that they introduce. However, such difficulties could be bypassed by adapting highly efficient parallel computing techniques in the supercomputing machines of EFDA infrastructures. By having successfully implemented the 3D version of the cold plasma code we succeeded to suppress the complexities of the coordinates transformations that the 2D models need and to work in the inclusion of the magnetic field component and viscosity terms as well.

1.7.3. Development of computational fluid dynamics solvers viscous for MHD flows (DEMO incl.)

Background and Objectives: The study of conductive and viscous fluid flow and transport phenomena in magnetic fields is of fundamental interest for the understanding of the underlying physics. Due to many particular problems related to the magnetohydrodynamic flows, for example, MHD turbulence, heat transfer, transition and stability, a series of individual accurate solvers based on computational fluid dynamics techniques are being developed.

Work performed in year 2010 (*in co-operation with the Institutes indicated*):

- MHD free convection of conductive fluids was studied in a closed vertical annular cavity in which the upper and bottom surfaces are adiabatic while the cylindrical walls are kept isothermal at different temperatures. Because of the applied horizontal external magnetic field, tangential velocities are developed in the fluid, resulting in a toroidal motion. The laminar and turbulent flow regimes were assessed by performing 3D direct numerical simulations. The results show that in the absence of the magnetic field, turbulent flow develops in most cases, while as the magnetic field increases the flow becomes laminar. The volumetric heating of the fluid resulted in higher temperatures in the upper-central part of the annular cavity, and in the creation of two convection currents as the hot fluid ascends in the central part and descends close to both colder walls. The Hartmann and Roberts layers developing near the walls normal and parallel to the magnetic field, respectively, were found to be responsible for the loss of flow axisymmetry. These results are reported in more detail in
**Annex 13**. (*In co-operation with ULB*.) - A heat transfer numerical module was developed and included successfully in both the Cartesian and Cylindrical versions of a computational fluid dynamics MHD code for incompressible unsteady flow. This was done by solving an extra partial differential equation for the thermal energy of the electrically conducting fluid in the form of a passive scalar quantity (i.e., temperature). Representative numerical simulations of fully developed MHD turbulent channel flows with heat transfer in channels and pipes were performed for the verification of the aforementioned computational numerical module. A reasonable agreement with other (experimental and numerical) results has been found. Direct numerical simulations and large eddy simulations of turbulent channel flows with heat transfer in the presence of uniform external magnetic fields have been conducted for various Prandtl and Hartmann numbers. In addition, a conditional sampling procedure was used to extract the dominant coherent structures in the near-wall region, focusing on the investigation of the magnetic field effect on the organised structures and the associated turbulent heat transfer. Several subgrid scale turbulence models have also been developed and tested by performing additional large eddy simulations (LES) in several representative flows. Results based on the study of 3D effects and particles have also been obtained (see
**Annex 14**). (*In co-operation with KIT.*) - An extensive parametric study was conducted with varying cavity aspect ratio, and Hartmann number for the base flow consisting of two interacting rolls, employing parametric continuation techniques. The parametric study on the effect of volumetric heating and aspect ratio on the dominant eigenmode was then carried out indicating the transition from the Goertler instability to a hyperbolic flow type instability as the heat production increases. The travelling wave mode was identified as the dominant mode in that case. An investigation of different approaches for parallelisation, i.e. frontal techniques vs. ScaLapack, was performed in order to optimise system matrix inversion. ScaLapck is opted for as an optimal choice that is especially suited for application on moderate-high Hartmann flows in square ducts (see
**Annex 15**). - A comparison of 3D linear stability predictions for the onset of 3D convection against experimental findings and predictions of quasi 2D linear studies was performed for Rayleigh-Benard flow in the presence of a magnetic field. It is shown that 3D stability provides more accurate predictions for moderately high Hartmann, Ha~200. The flow arrangement that is relevant to flow in breeder modules of HCLL blankets was identified and compared against available stability studies and the results of in house codes. Agreement in the stability thresholds as calculated with the latter methodology against available experimental studies is ascertained for Ha~200 and is seen to be superior over available predictions in the literature. Appearance of quasi 2D structures pertains to even higher Ha numbers (see
**Annex 16**).

1.7.4. Development of an immersed boundary solver for MHD flow for blanket modules (DEMO incl.)

Background and Objectives: This activity aims at developing an Immersed Boundary (IB) methodology for MHD flows of liquid metals in complex geometries. Such a technique can extend the predictive capability of existing flow solvers and the range of computable wall-bounded MHD flows at a reasonable cost.

Work performed in year 2010:

- The verification of the immersed boundary methodology has been completed for the two dimensional hydrodynamic flow around an unbounded circular cylinder against body fitted results. The verification was based on three-dimensional simulations for the hydrodynamic flow around a circular cylinder at low magnetic Reynolds number in a channel with insulating walls against three-dimensional body fitted DNS results. Detailed three-dimensional simulations were performed for a series of hydrodynamic and magnetohydrodynamic cases over a wide range of Reynolds and Hartmann numbers in order to validate the numerical methodology against a body-fitted DNS code. Flows in channels and ducts were tested to verify the developed immersed boundary solver.
- The MHD flow around a circular cylinder in channels and ducts with insulating walls was examined to study the effects of end walls, aspect and blockage ratios. Three-dimensional simulations around a circular cylinder in a channel and in a duct at Reynolds numbers in the range Re=100-400 and Hartmann numbers Ha=0.1-100 were performed with both a boundary fitted DNS code and the developed immersed boundary method code. An excellent agreement was reached between the two different codes. The effect of Re/Ha ratio on the flow steadiness is the subject of future studies.

Last Updated (Sunday, 01 April 2012 18:32)