## 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 (Energy, momentum and particle transport in the presence of coherent and turbulent spectra)

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 advanced 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 2012: (in co-operation with PSFC-MIT, USA):

(i) Studies on the evolution of distribution function in the action space on the basis a Fokker-Planck equation with time-dependent operator and on the basis of a functional mapping have been performed. The time-dependent operators and functional mappings have been constructed for various cases of discrete and continuous spectra corresponding to the presence of coherent structures and turbulent fluctuations. Their form has been studied for both cases with respect to the expected effect to the evolution of the distribution function for the transient phase as well as for the asymptotic, steady state.

(ii) Comparisons with the standard quasilinear approach as well as particle simulations of wave-particle interactions have been performed for a simplified 1-D case. More specifically, cases of periodic waves having discrete spectra as well as localized solitary waves having continuous spectra have been considered. A key issue has been shown to be the selection of appropriate values of the free parameters of the two schemes, which are a characteristic time interval for the time-dependent Fokker-Plank operator and the time step of the functional mapping. The comparisons will be continued in order to determine selection rules for these parameters depending on the perturbation spectra. The approach considered here has also been applied for the construction of functional mappings for the study of the evolution of the ion velocity distribution function under interaction with EC beat waves (see ANNEX 1). Comparisons with particle simulations have been performed for this case.

### 1.2.2 Extension of the Dispersed Phase Structure Dimensionality (DPSD) concept for conditions relevant to Edge Turbulence

Background and Objectives: Coherent turbulence structures play an important role in the dynamics of turbulence and modify the transport of (passive) scalars and contribute to the preferential concentration of transported particles. In the case of 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 approach based on the dispersed phase structure dimensionality tensor (DPSD), so that it can become applicable to plasma turbulence. In 2010 and 2011, the DPSD concept was tested using plasma turbulence data coming from the ATTEMPT code.

Work performed in 2012:

(i) Activities in the first part of 2012 were devoted to completing the validation of the use of the DPSD concept in edge turbulence and in the scrape-off layer. For this purpose, data from the ESEL code was used in order to compute the one-point turbulence structure tensors. Specifically, turbulence structures of density and electric potential fluctuations were identified using the DPSD concept. The predicted statistical characteristics of the turbulence structures were compared with visualizations of the actual fields coming from ESEL. In the second part of 2012, the software developed was optimized to deal with the larger grids used in ESEL. The code optimization will allow the integration of the DPSD module with the turbulence Consistent Physical Object (CPO).

### 1.2.3 Integration of the Dispersed Phase Structure Dimensionality with Mulitphysics Edge plasma model

Background and Objectives: Building on the outcomes of Task 1.2.2, this work is devoted to the development of a code that will allow the DPSD concept to be used as a diagnostic in existing edge plasma codes. Once the software is tested against representative experimental results, a Kepler actor will be created and made publically available.

Work performed in year 2012:

(i) The first steps were taken to make the DPSD code into a Kepler module capable of performing diagnostics. A Python Kepler module that can interact with the turbulence CPO to compute structure tensors was created in early 2012. In the later part of the year, the Python software was cleaned up and optimized, and initial work to integrate the ITM UAL and CPO concepts has begun. This work will be continued in 2013.

### 1.2.4 Validation of SOL turbulence model

Background and Objectives: It is important that SOL- capable codes are validated against experimental Langmuir probe data. The objective of this task is to develop a Fortran code to import data (from JET, ASDEX, or COMPASS, in cooperation with EDRG, and EFDA Physics actions concerning edge physics).

Work performed in 2012:

(i) The work carried out in 2012 is a continuation of work initiated in 2011. In early 2012, the basic infrastructure for the CPOs was extended to read the SOL parameters from the Langmuir probe data and a Fortran code to import JET data was developed. In the second half of 2012, new subroutines were developed for reading the electron and ion temperatures, and the electron density from the Langmuir CPO. Subroutines were also developed to compute the stroke parameters from the probe position data. The last closed flux surface is approximated from the equilibrium CPO data, and the probe position is recomputed in normalized units relative to the LCFS. These parameters are then used as input to a SOL turbulence code (ESEL) that outputs the results as an HDF5 file. With the completions of this work, we are keen to obtain experimental profile/equilibrium data from any experiment and to obtain standard cases with the imprimatur of ITER IO/F4E.

### 1.3 MHD stability and plasma control

### 1.3.1 ECRH/ECCD for MHD control

Background and Objectives: This activity aims to the improvement of the modeling of MHD stabilization with the use of localized ECRH/ECCD. The main target is to model the effect of magnetic island geometry on the EC power deposition and the consequence to the efficiency of NTM stabilization. The inclusion of this effect in wave codes is performed by computing the absorbed power density in terms of the volumes of the perturbed flux surfaces (currently available in the codes CODERAY and TORBEAM). Such an upgrade to the codes aids in estimating the minimum amount of current needed for NTM stabilization, based on the known criteria and the novel ECCD computation, and also opens the road for a fully self-consistent numerical model of the ECCD-driven NTM dynamics in terms of coupling the wave code with the solver of the NTM growth and rotation, in equilibrium geometry that includes the magnetic islands. Another issue studied is the deleterious effects of edge density fluctuations in ECRH/ECCD applications. In the experiments for ECCD-based NTM control in ITER, the beam has to propagate through the plasma edge where turbulence can be rather intense, so that the associated density gradient may induce wave scattering and diffraction. One may expect this effect is less important, but, as the distance to the resonance layer is long enough, even a small scattering angle could lead to considerable beam broadening and a degradation of the stabilizing effect. Our study consists of ray and full-wave modeling of the propagation assuming different types for the spatial distribution of the density fluctuations.

Work performed in year 2012: (in cooperation with IPP-Garching (i), Uni. Bayreuth, MIT (ii)):

(i) The implementation of the effects due to the magnetic island geometry in the beam tracing code TORBEAM, as well as the self-consistent loop of the ray tracing code CODERAY with the modified Rutherford equation, has been finalized. The modified version of TORBEAM includes the effect of the magnetic island only in the structure and magnitude of the flux surface volumes in the island region, in contrast to CODERAY which includes also the non-axisymmetric form of the magnetic field and the flattening of the plasma density and temperature profiles. The modified TORBEAM, running in ray-tracing mode, has been benchmarked successfully against CODERAY. Regarding the coupling of CODERAY with the modified Rutherford equation, the first results show that the stabilization occurs faster if the computation includes self-consistently the geometric effect mentioned above (see ANNEX 2). In this context, the coupling of the modified Rutherford equation with TORBEAM (actually an incorporation of TORBEAM in the MRE), together with a simple model for the island rotation, will be pursued in 2013.

(ii) The full-wave modelling of EC wave beam propagation in ITER in the presence of edge turbulence, using the FDTD full-wave code FWTOR, is in progress. The computation of the wave field and its spectrum in wave-vector space before and after crossing the blob has been completed, and the comparison to ray tracing results in a consistent manner is under current investigation.

### 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 implementation is based on the full non- linear, resistive and compressible MHD equations for the plasma flow that are being integrated in the numerical library OpenFOAM. The major objective, after the successful validation of the model, is the extension to a multi- fluid plasma flow equations model that will be able to consider possible impurity transport.

Work performed in year 2012 (in co-operation with CEA, ULB, ITM2, ENEA, MedC):

(i) The incompressible and fully non-linear MHD model that has been developed in OpenFOAM was integrated to a compressible one in order to account better for the fusion plasma. These developments concern mainly the adaptation and modulation for the ITER plasma geometry and flow conditions. In particular, the basic compressible resistive non-linear magnetohydrodynamic equations were coded into the object oriented library and we are now in the process to test the new module for the solution of simple compressible MHD flows that have analytical solutions. Numerical issues like the divergence-free magnetic fields especially in the compressible case where magnetic shocks are developed and propagated into the tokamak are being treated. Moreover, global mesh refinement of the three-dimensional numerical grid, and also, possible adaptive local refinement techniques based on strong local current gradients are being considered. A well-known 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 is on the way.

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

Background and Objectives: The suppression of (neoclassical) tearing modes by localised electron cyclotron current drive (ECCD) is a well-known technique, which by now has been successfully applied in every high-beta tokamak. It is still, however, a very challenging issue, given the remaining difficulties of the detection and the localisation of the magnetic islands. Magnetic islands could be localised, and then continuously targeted at the O-point by the EC beam, if they are locked to Resonant Magnetic Perturbations (RMPs) generated by internal (in-vessel) coils. The major milestones of the 2012 work-programme had been to achieve the locking of the m/n=2/1 (N)TM, by the n=1 RMP at a convenient toroidal position for continuous-wave ECCD at the stationary “O” point, and to prevent the mode formation, and the subsequent disruption, by pre-emptive continuous-wave ECCD at the expected position of the stationary “O” point.

Work performed in year 2012 (in cooperation with IPP-Garching):

(i) The disruptive m/n=2/1 magnetic island was successfully locked to the n=1 RMP in L-mode discharges (with low toroidal plasma rotation) at the toroidal fields of 2.5T and 1.9T, which are resonant to the EC frequencies of 140GHz and 105GHz. The locking position (toroidal angle) was controlled by the proper phasing of the current in the internal coils, and was also evaluated by a new locked-mode detection diagnostic. One of the several locking positions, which were feasible, had been also convenient for continuous-wave ECCD at the stationary “O” point.

(ii) With the application of pre-emptive continuous-wave ECCD at the expected position of the stationary “O” point, the disruptive m/n=2/1 magnetic island had grown to a much smaller amplitude, than in the absence of ECCD (see ANNEX 3) and the disruption was completely avoided. On the contrary with the application of pre-emptive continuous-wave ECCD at the expected position of the stationary “X” point, the disruptive m/n=2/1 magnetic island had grown to a much larger amplitude, and disrupted the plasma much earlier than in the absence of ECCD.

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

### 1.4.1 MD simulation of the He bubbles formed in Be and W under He bombardment

Background and Objectives: Beryllium is a candidate plasma-facing material, thus, it is important to understand its behaviour in H-He plasma. For this purpose, molecular dynamics (MD) simulation is a promising numerical technique. Simple numerical experiments completed in previous years (2009, 2010), examined only the case when incident particles bombarded a virgin Be crystal. In 2011, this work was further extended by looking at the effect of successive collisions. During the work carried in 2011, we had observed at the end of long simulation runs, when the MD Be cell saturated beyond some threshold by He -atoms, the formation of He-bubbles. These He-Bubbles eventually led to the sputtering of big Be-fragments that contained few dozen Be-atoms, in contrast to single He-atoms sputtered from the MD cell at low He-concentrations. The time to reach the saturation of the Be-cell was found to be dependent on the incident He-fluxes. The main goal in 2012 was to clarify whether these previously obtained molecular dynamics (MD) are universal or dependent on the particulars of the MD simulations. In MD simulations, the frequency of issuing He -atoms towards the Be-surface is artificially high; this is dictated by practical considerations, otherwise, because of limited computational power, the simulation is not feasible. The implication of this fact is that the simulated incident fluxes are in 5 to 6 orders of magnitude higher than in reality. However, if it can be shown that during the time gap between adjacent He-Be collisions, the MD Be-cell reaches equilibrium, then, one can expect the difference between the theoretical and experimental incident fluxes to be unimportant. Under these conditions, one may rely on the MD results to draw useful conclusions. The answer to this question is important, because if the He-bubble formation is proven to be artificial, the peeling effect is unlikely to happen in reality. Otherwise, special care is needed using Be as a plasma-facing material. Thus, during 2012 our MD simulations were directed towards clarifying whether the observed He-bubble formation is an artifact of the exceedingly high He-fluxes in our simulation or not.

Work performed in year 2012:

(i) Systematic MD experiments, requiring very long simulations times, were performed for large MD-cells and various and boundary conditions. The main achievement of the study is that we have checked the destruction of the Be-surface, successively bombarded by He-atoms, at various conditions, such as different BC size of the MD-cell, and different He-Be interaction potentials. We have confirmed that the destruction of the cell starts at almost constant ratio between the number of He-atoms accommodated by the Be-cell and the total number of Be-atoms. Also, the time to reach destruction scales linearly with the size of the MD-cell (see ANNEX 4).

(ii) Through a series of MHD experiments, we have also confirmed the general atomic scenario for Be swift erosion that has been previously suggested as a possible mechanism for Be erosion.

(iii) For the large Be-cells we have not reached destruction barrier with the multibody Tersoff potential due to computational limitations: one integration step with the Tersoff potential is few orders of magnitude more expensive than that of ordinary pairwise potential. So, at the moment it’s impossible to state whether the destruction is also possible with the Tersoff potential. The effort to clarify the effect of using the Tersoff potential will be continued in 2013.

### 1.5 Physics of plasma heating and current drive

### 1.5.1 ECRH/ECCD and transport in tokamak plasmas

Background and Objectives: In this activity we investigate the electromagnetic wave propagation, absorption and current drive in the EC, IC and LH regime, as well as the resonant interaction of the wave with the plasma electrons/ions in geometry relevant to heating, current drive and diagnostics applications. Our research is conducted in several parts: The first part focuses on the role of wave effects (beam localization, diffraction, interference and mode conversion) and plasma effects (hot-plasma dispersion, nonlinear absorption and non-local response) to the propagation and absorption of electromagnetic waves. Such effects could possibly be important in many cases of interest, like incomplete cyclotron-resonant absorption, high-power ECRH/ECCD and mm-wave reflectometry/interferometry, where the conventional formalism involving asymptotic wave methods (ray/beam tracing) and local linear dielectric plasma response may break down. A second part is to clarify the advantages of the utilisation of non-Gaussian beams in ECRH systems, formulated in terms of the paraxial WKB beam tracing method for following the wave propagation and the linear hot-plasma dielectric 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. A third part has been established on the assessment of the capabilities of EC H&CD in DEMO on the basis of power balance calculations with the transport code ASTRA and fully-relativistic computations of the wave deposition and current drive with TORBEAM. For defined sets of DEMO parameters, the plasma profiles required for the wave computation and the part of the current which should be driven by external means are computed by ASTRA. Then, one synthesizes the required current profile for the nullification of the loop voltage profile, followed by the determination of the minimum-power setup that leads to non-inductive equilibrium. Another very important issue is the coupling of EC waves to thermonuclear plasma which is widely considered to be very efficient: Even a relatively small (of the order of a few %) fraction of uncoupled radiation can lead to machine damage and reduced efficiency. Therefore there is a need to estimate the potential for machine damage and efficiency loss due to density perturbations near the plasma edge and due to mode impurity and non-ideal polarisation of launched EC waves. The density fluctuations and blobs present in the edge region of magnetic fusion devices can scatter radio frequency waves through refraction, reflection, diffraction, and coupling to other plasma waves. This, in turn, affects the spectrum of the RF waves and the electromagnetic power that reaches the core of the plasma. There are two basic approaches to studying the effects of density blobs and edge fluctuations on the coupling of electron cyclotron (EC) radio frequency waves to the core of tokamak plasmas. The first is the geometric optics approach in which the effect of fluctuations is to change the refractive properties of the EC beam or rays. There are two consequences of refractive scattering – diffusion in real space leading to a spatial deflection of the rays and diffusion in wave vector space leading to the broadening of the launched spectrum. The geometric optics approach is limited to small density fluctuations of 10% or less. The second approach to studying the effect of blobs on EC fields is using the full wave approach. This approach extends the range of validity well beyond that of geometric optics; however, it is theoretically and computationally much more challenging.

Work performed in year 2012 (in cooperation with MIT (i) and IPP-Garching (ii, iii)):

(i) The 2D FP model has been further developed with the inclusion of a temporal delay in the associated time-ependent RF diffusion operator as the original kinetic description and basis of the formulation of the latter predicts. Other numerical techniques, apart from the Chebyshev pseudospectral methods, such as finite difference methods has also been employed. For the 3D implementation an efficient tool for handling the transformation from velocity and configuration space description to the canonical action-angle description was needed. Therefore, late in 2012 we have started investigating along this direction which has been proved quite non-trivial. The parallel version of the FP codes has benn, therefore postponed for 2013.

(ii) Since detailed runs for ITER/DEMO relevant parameters are expected for the full 3D implementation and the construction of the aforementioned transformation tool has been prioritized, these runs are expected to be initiated within 2013.

(iii) The problem of scattering of a cylinder the eigenfunction expansion method has been formulated (see ANNEX 5, ANNEX 6).

(iv) Various ways of field visualisation has been explored for the case of scattering off a spherical blob. The estimation of the spectrum of the radial component of the Poynting vector was considered as the most informative especially as far as the forward and backward scattering is concerned. Results are summarized in ANNEX 7 where, for ITER relevant densities, size and magnetic field, the Poynting flux is presented (logarithmic scale; normalized to the incident flux of an O-mode at about 90 degrees) on the z-x plane (z is the direction of the magnetic field). It is evident that the forward scattering becomes more "focused" as the size increases along with an intensified activity at the "poles". Activity on the blob surface is intense as expected. This increased EM activity (especially along the magnetic field lines for relatively larger sizes) may trigger parametric effects. Extension to an RF-beam (multiple k's) has also been initiated.

(v) Numerical implementation for the cylindrical blob has been postponed for early 2013 due to other commitments within WP2012.

(vi) The evolution of ion distribution functions in the presence of two beating X-mode electromagnetic wave is analyzed. When the amplitude of the wave is small, the nonlinear interaction with the envelope of the modulated wave can result in a coherent energization of the ions. The purpose of this work was to investigate the values of the parameters of the modulated wave that optimize the energy transfer from the wave to an ensemble of ions as well as to determine the amplitude threshold for chaotic motion and the extent of the chaotic orbits on the phase space.

(vii) Through numerical simulations we found out that the energy exchange is sensitive to a small detuning between the beat frequency and the ion cyclotron frequency, which is increases with increasing carrier frequency. The effect of the group velocity of the waves also found to be significant. In most cases the energization is faster and more efficient, when the envelope velocity is equal to the phase velocity of the waves. A comparison of the energization process for cold IC beat waves has been carried out and is presented in the ANNEX 1. For larger amplitudes, chaotic energy transport takes place. The wave polarization seems to play a crucial role on both the wave amplitude threshold for chaotic transport as well as on the maximum energy acquired through this process.

(viii) FP kinetic modeling of alpha particle dynamics has been initiated before the end of 2012. Tests-comparisons with what have been reported so far in the literature have been postponed for 2013.

(ix) In this period, the development and improvement of the full-wave code FWTOR, based on the FDTD method, was continued. We progressed on the parallelisation of the code and initiated the development of the linear dielectric response tensor for ions in interaction with IC and LH waves, in the scope of the simulation of RF heating, current drive and diagnostics in AUG and ITER. Taking into account the properties of the computers available for simulations within the fusion community (HPC-FF, IFERC-CSC), it was decided that a mixed MPI/OpenMP logic should be followed. A few details remain to complete the parallelization, and then the benchmark against the asymptotic codes of the IMP5 project will be put forward, and then initiate the extension of the code to the IC/LH range of frequencies.

(x) During the current period, a series of amendments took place in the pWKB beam tracing theory for the evolution of non-Gaussian beams, in the spirit of clarifying the physics of the computation of the phase-shift and the amplitude of the higher-order modes. The computation of the phase-shift incorporating the novel theory findings in the numerical codes and benchmark the results, was put behind in order to devote more time on the activity (iii), which is a PPPT PS task.

(xi) For the current period, the work was to synthesize, with a given ECCD system, the required current profile to achieve the nullification of the loop voltage over the whole radial profile, and then derive the required wave power to drop the loop voltage to zero from its value corresponding to the initial equilibrium. The kinetic plasma profiles, and hence the profiles of the bootstrap current and the electrical conductivity, had to be fixed in this approach, so any effects on the equilibrium from the external heating power were not simulated. The results indicate that the usual localized structure of the ECCD profile (as met e.g. in experiments for MHD stabilization) is not preferable when looking for long-pulse or steady-state operation. On the contrary, the current density profiles have to be much broader, even spanning along the total of the poloidal region, and also appear a peak on the plasma centre. Such profiles could be constructed by injection of broader, weakly-focused beams or combination of more than one wave beams in order to improve the control in the formation of the targeted profile (e.g. one beam for the region close to the plasma centre and another one for the edge). A first set of computations in this direction has been completed (see ANNEX 8).

### 1.6 Energetic particle physics

### 1.6.1 Modelling the dynamics of alpha particles from thermonuclear reactions using fluid approximations

Background and Objectives: The long term objective of this task is to develop a 2-D model in R-Z (cylindrical) geometry of burning plasmas, where the initial distributions of the plasma parameters and the magnetic topology will be taken from an equilibrium. The model equations are the fluid equations of the charged species constituting a burning plasma, derived by taking moments of Boltzmann’s equation, supplemented with Maxwell’s equations. This approximation is usually called Two-Fluid approximation, because of negative and positive charged fluids. Here electrons, ions (D-T), and alphas are treated as three different fluids interacting with electric and magnetic fields. In order to explore and solve the numerical problems that will be encountered with the development of the 2-D code (axisymmetric geometry) for burning plasma, first we develop a 1-D code which computes the model equations of a burning plasma in 1-D cylindrical geometry.

Work performed in year 2012:

(i) A one dimensional cylindrical geometry code which treats electrons, (D-T) ions and alphas as three separated fluids has been developed. The system of equations, which the code solves, are the conservation equations for each fluid (continuity, momentum and energy) supplemented with Maxwell’s equations. The alphas are generated by the thermonuclear reaction of D-T. A brief description of the code and preliminary results are given in ANNEX 9. In order to benchmark the code, shock tube tests have been performed with the code, operating without the alphas production (two-fluid approximation). The results from these shock tube tests have been compared with shock tube tests performed with one dimensional MHD (single fluid) code, and are in good agreement. Momentum transfer and temperature equipartition between electrons and ions have been implemented in the code, collisions with alphas is in progress.

### 1.7 Theory and modelling for ITER and DEMO

### 1.7.1 Application of computational tools for liquid metal Blanket modules

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 these potentially strong magnetic fields, 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 turbulent liquid-metal flow when entering low, strong and very strong magnetic fields. Furthermore, the placement of flow inserts in the blanket ducts can modify the stability and heat transfer characteristics of the flow and for design purposes it is desired to quantify these effects in the presence of strong magnetic fields. The analysis in this work is entirely done using a parallel, unstructured, Navier-Stokes solver developed at University of Cyprus. Earlier work (2010, 2011) has focused on the flow stability during entrance or exit from a strong magnetic field. Starting in 2012, focus has been shifted to the stability of duct MHD flow in the presence of flow inserts.

Work performed in year 2012:

(i) Work in early 2012, was devoted to completing preliminary three-dimensional Direct Numerical Simulations (direct solution of the governing partial differential equations, no modelling involved) of MHD duct flow in the presence of a cylindrical insert. The aim of these simulations was to establish appropriate domain size and grid resolution requirements for the desired range of flow parameters, i.e. Reynolds (Re) and Hartmann (Ha) numbers. These simulations were also used to construct a rough 2D map of flow regimes as a function of Re and Ha that will help guide subsequent high-fidelity simulations.

(ii) In the second half of 2012, a detailed parametric investigation of flow regimes and stability regions over a wide range of Ha numbers for MHD flow around confined cylinders in rectangular ducts was initiated. In particular, a series of simulations was devoted to the detailed analysis of two unexpected and counter-intuitive effects that were observed during the preliminary simulations in early 2012. One of these effects is a non-monotonic dependence of the critical Reynolds number (Re) for the onset of vortex shedding, with respect to the Hartmann number (Ha). The second effect is an increase in the flow unsteadiness with increasing intensity of the magnetic field, associated with the onset of newly identified flow regime. A detailed study of the force coefficients on the flow insert as a function of Ha and Re has led to an unexpected result: the spanwise distribution of the force coefficients along the cylinder is found to become more three-dimensional with increasing Ha. This work will be continued and completed in 2013.

### 1.7.2 ITER pertinent equilibrium and stability of fusion plasmas

Background and Objectives: This activity aims at constructing equilibria and relaxed states of fusion plasmas with flow or/and finite electrical conductivity and viscosity and investigating their linear and nonlinear stability in connection with the ITER project.

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

(i) The parallel-flow version of the HELENA code was benchmarked by means of the Solovèv solution. In parallel the development of the arbitrary-direction-flow version of the code was initiated. The input/output parts were modified, though not in full compliance with the requirements of the ITM-TF due to the absence of the mass density as flux quantity in the CPOs. The solver was also partially modified in order to include the additional term associated with the electric field (see also [Kuiroukidis and Throumoulopoulos, Proc. 39th EPS Conf., 2012]). (In co-operation with IPP.)

(ii) We have obtained ITER-relevant MHD steady states of translational symmetric and axisymmetric plasmas with incompressible sheared flow and examined their stability. The study includes equilibria either with monotonically increasing safety factor profiles pertinent to the L-H transition or with reversed magnetic shear. Specifically, linear and nonlinear solutions of a generalized Grad-Shafranov equation have been constructed analytically, quasi-analytically and numerically. For non parallel flows it turns out that the electric field makes the safety factor flatter and increases the magnitude and shear of the toroidal velocity in qualitative agreement with experimental evidence on the formation of Internal Transport Barriers in tokamaks, thus indicating a potential stabilizing effect of the electric field. For parallel flows the linear stability has been examined by applying a sufficient condition. In this case one equilibrium corresponding to the H-state is potentially stable in the sense that the stability condition is satisfied in an appreciable part of the plasma while another solution corresponding to the L-state does not satisfy the condition. The results confirm the conjecture that the sheared flow in conjunction with the equilibrium nonlinearity play a stabilizing role. Details are given in ANNEX 10, ANNEX 11 and [Throumoulopoulos and Tasso, Phys. Plasmas 19, 2012], [Tasso and Throumoulopoulos, J. Plasma Physics 78, 2012]. (In co-operation with IPP and MEdC.)

(iii) An equilibrium with ITER boundary obtained by the HELENA code was used as initial state for modeling the dynamics of alpha particles from thermonuclear reactions using fluid approximations. Also, during the evolution of the system, collision frequencies for Coulomb interactions between electron-alpha particles and deuteron-alpha particles are involved with different electron and ion densities and temperatures. Calculation of these collision frequencies is under way (In co-operation with FORTH).

### 1.7.3 Development of computational fluid dynamics solvers for viscous MHD flows

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 underline 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 2012 (in co-operation with ULB, KIT, UCLA-USA, IIT Kanpur-India):

(i) A numerical model was developed for the study of the transition of the electromagnetically driven flow of an electrically conducting fluid in a toroidal square duct. The flow transition was mapped for a wide range of Hartmann and Reynolds numbers (see ANNEX 12). The results show two types of Taylor-Couette instabilities occurring in this flow with value Ha = 18 separating each instability. For values Ha < 18, the instability occurs when the double Taylor-Couette vortex pair becomes unstable, while for Ha > 18 the instability is connected to numerous Taylor-Couette vortices that are generated near the sidewalls, especially the concave one. The increase of the mean azimuthal velocity of the fluid (i.e. Reynolds number) destabilizes differently each regime as the fluid inertia becomes dominant. At transition, as the Reynolds number is increased, only localized disturbances near the concave sidewall layers were found in the present small duct. For the higher Hartmann numbers of the current simulations, the critical Re for transition was found at Re=32000 for Ha=200, and Re=112000 for Ha=500. A part of the results has been already published in [13] and two other publications will be submitted soon.

(ii) New features of the turbulent Kolmogorov flow was found for elongated periodic domains where the translational invariance in the streamwise direction appears to be broken and the turbulent statistics found to depend on the aspect ratio of the computational box. In the light of the new findings, the Kolmogorov flow was studied here in the frame of statistical ensemble large-eddy simulations. In particular, statistics of an ensemble of up to 96 independent simultaneous simulations of the flow was collected. The only communication permeated between the simulations was made during the determination of the dynamic Smagorinsky model constant which is evaluated and averaged only between the ensemble. In contrast to the usual practice of the large eddy simulations of the homogeneous turbulence, no local averaging in the spatial directions of each simulation was allowed. The simulation in each instance was initiated either from turbulent hydrodynamic flows and large scale magnetic fields where applied leading to the usual transition to quasi-2D MHD turbulence. Results found to be sensitive to the initial conditions in the sense that the dominant statistical flow structure maybe found as in direct numerical simulations or maybe erased by the averaging between the ensemble. In the present stage we are collecting and post-processing all available data and in the near future we will submit them for publication.

(iii) The effect of an axial magnetic field on the transient laminar natural convection cooling of a low-Prandtl number hot fluid contained in a cold vertical cylinder was studied. Direct numerical simulations were carried out for different values of Prandtl, Rayleigh, and Hartmann numbers at which the flow remained always laminar. It was found that in the presence of the magnetic field, the fluid motion is damped resulting to a deceleration of fluid cooling. The magnetic field has little influence on the initial stage of boundary layer development on the cylindrical wall. The vertical thermal boundary layer becomed fully developed very quickly at times that are independent of Hartmann number for low Prandtl and Rayleigh numbers, while longer times are generally required with increasing Hartmann number for large Prandtl and Rayleigh numbers. For large values of Pr and Ra, the cold fluid intrudes from the vertical to the bottom wall and, finally, to the bulk fluid through the central vessel area, opposing the development of the thermal boundary layer on the sidewall. Due to the damping action of the magnetic field, the cold fluid intrusion is significantly diminished and the thermal stratification is delayed with increasing Ha. Consequently, conduction dominates over convection and heat transfer is reduced, which is more pronounced for high Prandtl and Rayleigh numbers. A larger percentage of the cylinder remains at higher temperature and, thus, the cooling process lasts longer relative to the pure hydrodynamic case. Finally, the presence of the magnetic field does not play a great role during the ending of the fluid stratification process. All the results of the milestone were already published in Numerical Heat Transfer journal. Simulations of an MHD free convection 3D liquid metal flow in annuli were successfully performed using an available in-house code. Different aspect ratios were investigated to determine the A.R. influence in heat transfer and MHD effects. More specifically, a series of direct numerical simulations for the Rayleigh numbers Ra=104 and 105 and Hartmann numbers Ha=0 to 100 were performed with respect to four different aspect ratios (ratio of the outer cylinder length to height) 0.5, 1, 2 and 3. Depending on the value of Ra and Ha numbers and the choice of the aspect ratio, laminar, transitional or turbulent flows are developing (see ANNEX 13). In parallel, the target of the assessment of the open source platform open effect, in potential gradients on jet deflection was ascertained, Lorentz force or Foam in terms of identifying its capabilities in simulating MHD flows was successfully accomplished and the development of an MHD solver was also achieved and implemented to verify the validity of analytical solutions derived for a two-dimensional MHD natural convection flow of an electrically conductive fluid in an internally heated horizontal shallow cavity in the presence of an external vertical magnetic field. The analytical solutions were derived using the method of the matched asymptotic expansions. Both the analytical and numerical results demonstrate that the fluid is decelerated by the external magnetic field leading to the dominance of conduction over convection and, therefore, reducing the heat transfer (see ANNEX 14).

(iv) A linear second-order ODE containing a small parameter multiplying the second derivative, thus generating a boundary layer near one of the boundaries, was solved via the standard finite element methodology using Lagrangian basis functions. The exact error was obtained upon comparing against the exact solution. Next, an a-posteriori error estimate was obtained based on the symmetric part of the operator and performing an element-wise calculation using the residual method with flux equilibration and a higher order finite element description. The mesh was subsequently adapted in regions of the domain where the error estimate acquired values beyond a certain predetermined threshold with respect to the numerical solution. Increased convergence rates were compared and the error estimates were compared with the actual numerical error (see ANNEX 15). The analysis of two dimensional problems was not performed and was left for a separate study in the near future.

### 1.7.4 Development of CFD solvers for liquid-metal flows relevant to blanket modules (DEMO incl.)

Background and Objectives: Computational investigation of novel concepts pertaining to blanket design (DEMO incl.). Several multi-physics flow solvers for liquid-metal blanket modules will be developed in order to investigate a number of related technical issues such as, T-transport and permeation, wall roughness and corrosion, impurities deposition, and cleaning processes.

Work performed in year 2012 (in co-operation with ULB, KIT, WP12-PPPT):

(i) In the context of the WP12-PEX-04 EFDA Task Agreement and the ‘Novel PFC material solutions - Liquid Metals’ research, as the latter is defined in the sections of ‘Proposals for improving technological and scientific knowledge (experiment and modelling)’, our main motivation and as a long term target is the development of a free surface solver capable to model such flows, i.e. MHD liquid metal wall (film liquid metal) flows as alternative to metallic plasma facing components. In this framework and for the current year our main target was to investigate the advantages and disadvantages of some of the available numerical methodologies in simulating free surface flows by performing numerical studies on simplified flow arrangements using an in- house developed code and the open source platform open Foam. Indeed, an assessment of the available free surface solvers were performed exploring different numerical strategies by performing comparisons of the available free surface solvers with experimental data (see ANNEX 16). Based on the preliminary one-dimensional analysis that was carried out, the effect of external electric r r j × B conjunction with the importance of electric stresses in accelerating the jet and reducing its radius. As a result it is anticipated that the jet will be susceptible to an interfacial instability, equivalent to the Rayleigh instability but now due to the interfacial electric stresses, that will generate drops. Stability analysis is required in order to identify the parameter range for which electric stresses will overcome capillarity and pinch–off droplets. The motion of a single emerging droplet was subsequently determined by a similar force balance as for the jet where, despite of their smaller size and owing to spherical symmetry the deflection of the trajectory was found to be of the same order of magnitude as for a jet. The effect of shape asymmetry on the jet/drop trajectory is planned as a future study, aiming at minimizing the deflection of the trajectory (see ANNEX 17).

(ii) The corrosion mechanisms and T-transport in liquid-metal flows for fusion application were reviewed. Several corrosion models were developed and included in an MHD CFD solver in order to analyze the influence of the parameters that affect the corrosion behavior, and the transport and deposition of the corrosion products in simplified concepts of proposed liquid-metal blankets. The MHD solver was further developed by accounting also for T transport and for the electric coupling of a multi-material domain. A parametric numerical study to assess the performance of the corrosion models through comparisons with experiments in 2D MHD flows in rectangular ducts was conducted. A good agreement of the present numerical predictions with experimental and other numerical results was generally observed. The higher corrosion rate in the presence of a magnetic field was be explained by the steep velocity gradient in the Hartmann layer, while an increase of the mean velocity leaded also to an increase of the mass losses due to the corrosion by the liquid-metal flow (see ANNEX 18). An analysis of liquid-metal corrosion based on 3D numerical simulations has also been initiated during 2012. Due to the extremely large demands in computational grids and CPU time, the 3D simulations will be continued in 2013.

(iii) The Rayleigh-Benard case was examined in the presence of a magnetic field with increasing intensity. Linear stability analysis that was developed (see ANNEX 19) captures GrCr~Ha2 with great accuracy in comparison with the quasi-2D analysis findings as well as the experimental recordings pertaining to the onset of steady convection. For large side wall conductivity cS, it was seen that GrCr depends mainly on the Hartmann wall conductivity, cH, as Ha increases; cS, cH correspond to the electric conductivity ratios between the wall material comprising the side (stainless steel) and Hartmann (copper) walls, respectively. In particular, as cH is increased then GrCr is also increased and this finding is important concerning the stabilization of the liquid metal (NaK) flow and ultimately the optimization of heat transfer in rectangular cavities. It is also important to note that the experimental finding of time-dependent convection cannot be attributed to a thermal instability, since the present stability analysis does not support the onset of time-dependent modes.

(iv) Nonlinear analysis is carried out of free convection in a rectangular cavity that is heated from below inside a strong magnetic field, seeking a description of the nonlinear evolution of unstable modes beyond criticality. Coupling of the finite element methodology with a spectral approach is employed for the discretization of the above flow arrangement. Furthermore, a semi-implicit time integration scheme is employed, with the nonlinear convective terms treated explicitly using the second order Adams-Bashforth method and the linear terms treated in an implicit manner using the second order accurate Crank-Nicolson scheme. Thus, decoupling of the different spectral modes is possible while favoring parallel treatment of the solution process. The most time consuming part of the algorithm is that of the evaluation of the different Fourier modes and involves matrix inversion, different parallelization strategies are employed, i.e. in the physical and in the Fourier space. In order to optimize in terms of storage and CPU time requirements the ScaLAPACK (Scalable Linear Algebra PACKage) library and FFTE (Fast Fourier Transform East) package are employed. The MPI (Message Passing Interface) and BLACS (Basic Linear Algebra Communication Subprograms) protocols are adopted in order to coordinate communication between different computational nodes. Domain decomposition is also performed by every processor corresponding to different regions of the cavity cross section and consequently different parts of the banded system matrix. It is important to note that the transformations between physical and Fourier space are equivalent to a parallel matrix transposition and are achieved using the ALLTOALL command of the MPI library (see ANNEX 19).

(v) Due to the significant amount of time required for the other tasks, the work on benchmark calculations in the nonlinear regime was not initiated in 2012.

### 1.7.5 Design of integrated pumping system

Background and Objectives: Develop an integrated algorithm for simulating vacuum pumping systems of any complexity in the whole range of the Knudsen number and implement this s/w tool in the simulation of the pumping systems of ITER as well as of other DT fusion reactors.

Work performed in year 2012 (in co-operation with KIT, CEA, WP12-IPH):

(i) The recently developed in-house algorithm for solving gas distribution systems consisting of pipes of any length and cross section under any vacuum conditions has been successfully implemented to simulate a sample network representing the main ITER divertor pumping system. However, all provided simulation results are based on the 2008 geometrical characteristics of the ITER lower port region and divertor configuration (see ANNEX 20). The latest design and drawings of the aforementioned modules of the ITER reactor have been delivered to UoThly in October 2012. The translation of the flow path inside the cassette and the divertor regions into a network of channels has been initialized and is currently under way. It is expected to be completed early next year and then, once the geometry of the network is set, the network algorithm will be used to run the most recent design of the ITER divertor neutral distribution system and to perform some comparisons with ITERVAC.

Last Updated (Tuesday, 25 February 2014 10:39)