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

(i) Studies on transport of particles interacting with a statistically distributed set of localised wave structures have been performed. Averaged momentum and position variations have been calculated for ensembles of particles for the case of localised structures with parameters within a statistical set. The averaging is over the particle position along the magnetic field and the gyro-angle as well as over the respective statistical parameters of the wave structures. For simplicity of the calculations, Gaussian spatial profiles of the wave structures have been considered with statistically varying spatial width along and across the magnetic field, amplitude, phase and group velocity and temporal width. The interaction has been further considered as a random-walk problem and a Fokker-Planck equation has been derived. The actual duration of each interaction of a particle with a wave-packet, necessary for the calculation of the respective diffusion coefficient, has been calculated. The latter has been shown to depend on both the initial particle velocity as well as the statistical properties of the group velocity, spatial and temporal width of the localised structures.

(ii) Studies on the effect of turbulent wave-spectra on barely circulating and barely trapped particles in a tokamak have been performed. The role of turbulent wave-spectra acting as perturbations on the separatrix between circulating and trapped particles has been investigated. The deformation of the separatrix and the respective formation of heteroclinic tangles by the stable and unstable manifolds of the saddle points have been studied in terms of the Melnikov’s method. A canonical mapping for the particle dynamics in this area has been constructed for the description of the chaotic particle trapping and de- trapping. This mapping can be further used for the calculation of particle and momentum transport across passing-trapped boundary.

(iii) We have investigated the interaction of ions in a uniform collisionless magnetised plasma with extended wavepackets. The wavepackets comprise of two beating X-mode electromagnetic waves propagating perpendicularly to the magnetic field. Using Lie perturbation theory, we have calculated a second order mapping describing the evolution of the ion distribution function (see **Annex 1**). The results of the application of this mapping show very good agreement with numerical simulations of ion orbits. The purpose of this work is to investigate the values of the parameters of the wavepacket that optimise the energy transfer from the wave to an ensemble of ions.

### 1.2.2.Extension of the Dispersed Phase Structure Dimensionality (DPSD) concept for conditions relevant to edge turbulence (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 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 (DPSD), so that it can become integrated in the edge code ATTEMPT. In 2010, the DPSD concept was tested using plasma turbulence data by generating instantaneous profiles of density and electric potential fluctuations and interpolating them on to a simple toroidal grid (r,theta,psi). Turbulence structures of density and electric potential fluctuations were calculated on this grid. The objective in 2011 was to extend the initial validation work from 2010 to higher resolution grids, and verify the results.

Work performed in 2011:

(iv) During 2011, proof-of-concept simulations were completed for the case of a dispersed phase in MHD turbulence. These simulations verified the suitability of the DPSD concept for use as a diagnostic for the case of particle dispersion in MHD turbulence.

(v) Following these encouraging results, the DPSD concept was applied to model structures in edge turbulence using the software called ATTEMPT. Turbulence structures of density and electric potential fluctuations were identified using the DPSD concept. During the first half of 2011, simulations were done using a grid size of 16x16x32, and during the second half of 2011, simulations were done using grid sizes of 32x32x64 and 64x64x128. The main results from the three simulations are similar, thus showing that the DPSD computations are grid independent. The code needs to be sped up further because the computation of the DPSD tensor during 64x64x128 simulations takes a long time. Thus, work in late 2011 and early 2012 is focused on optimising the DPSD implementation.

### 1.2.4.Modular particle transport model (numerical tools for a comprehensive ITER modeling)

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 computations 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 2011:

(i) TRIMESH, a triangular mesh generator that constructs unstructured aligned meshes in a tokamak vessel geometry has been developed and made available on the ITM gateway. TRIMESH is a combination of a Python module and the open source mesh generator Triangle. Input parameters of the code are the vessel geometry and a snapshot of equilibrium data in the vessel. These input parameters are supplied as external data files. The main requirement that TRIMESH fulfills is to produce mesh elements aligned with the equipotential surfaces in a tokamak vessel. Because this requirement is very restrictive, special algorithmic approaches have been developed in order to construct aligned meshes. The algorithm starts with a file containing equilibrium data. Based on these data, we try to find an optimal arrangement of contour lines that is afterwards introduced in the Triangle mesh generator as an input parameter to produce the mesh of the vessel. The choice of the equipotential lines is very crucial because these lines will guide the Triangle mesh generator to fulfil the alignment. The Python part of the TRIMESH constructs the required contour levels in stages. In the first stage, we plot one contour level by splitting the contour range in two equal parts; in the second stage, we plot three contour levels by splitting the contour range in four equal parts, and so on. In this manner, at each stage, we produce contour lines that are placed between the lines from the previous stage. Using this method, we can arrange the contours in a hierarchy and we get better coverage over the entire domain, and subsequently a more aligned mesh. An additional parameter that TRIMESH takes into account is the relation between the maximum triangle area, and the spacing between the nodes. This allows the reduction of the number of nodes that Triangle introduces and leads to a better-optimised mesh.

### 1.2.5.Integration of the dispersed phase structure dimensionality with mulitphysics edge plasma model (development of multi-physics computational models for plasma simulations)

Background and Objectives: 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 dispersed phase structure dimensionality tensor (DPSD) as a structure diagnostic tool, so that it can become applicable to edge turbulence, in conjunction with 1.2.2. While in 1.2.2 the DPSD was only computed aposteriori for validation purposes, in 1.2.5 the DPSD module was integrated in the code so that the DPSD diagnostic is computed in real time during the calculations.

Work performed in 2011:

(i) During the later part of 2011, the structure tensor computation code has been partially implemented in ATTEMPT in order to allow DPSD calculations in real time during computations. Further work on this will continue in 2012. Preliminary tests completed in 2011 indicated that the code, as currently implemented, takes too long to compute on the 64x64x128 grid, while the results were not very different from the 32x32x64 grid. A strategy currently being explored in order to improve the speed of the DPSD computation is to use only alternate grid points for the structure tensor computation on higher resolution grids.

## 1.3 MHD stability and plasma control

### 1.3.1.ECRH for MHD control

Background and Objectives: This activity aims at the improvement of theory, modelling and the interpretation of experiments on MHD stabilisation with the use of localised EC wave heating and current drive. Research is conducted in two parts:

In the first part, the main target is to model the effect of magnetic island geometry on the EC power deposition and current drive and the consequence to the efficiency of NTM stabilisation. The inclusion of this effect in wave codes is performed on the basis of computing the deposited power density in terms of the flux volumes of the perturbed flux surfaces, and it is currently performed in our ray-tracing code CODERAY and the beam-tracing code TORBEAM. Such an upgrade to the TORBEAM code will aid in the theoretical prediction of the minimum amount of current needed for NTM stabilisation, in terms of the proposed criteria for jECCD/jbs and the novel ECCD computation. Furthermore, this opens the road for a fully self-consistent numerical model of the ECCD-driven NTM dynamics in terms of coupling the wave propagation solver with the equations for the description of the NTM evolution, in equilibrium geometry that includes the non-axisymmetric magnetic islands generated by the NTMs.

The second part is devoted to the investigation of the deleterious effects of edge density fluctuations in ECRH/ECCD applications. In the experiments for ECCD-based NTM control in ITER, the wave 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 to be less important, however, since the distance to the resonance layer is long enough, even a small scattering angle could lead to considerable beam broadening and result in degradation of the ECCD stabilising effect. Our study consists of ray and full-wave modelling of the propagation assuming different types for the spatial distribution of the density fluctuations.

Work performed in year 2011 (in cooperation with the Institutes indicated):

(i) In the current period, a generalisation of the analytic calculation of the volumes between adjacent magnetic flux surfaces was performed, on the basis of introducing an update to the existing magnetic field model, which represents the island in a more realistic way and incorporates inhomogeneities of the background magnetic field such as the 1/R dependence of the toroidal component, as well as of a new definition for the flux-surface- label function under the assumption of arbitrary magnetic equilibrium (see **Annex 2**). The new flux-surface labeling, which is based on the isostathmic contours of the perturbed poloidal flux function, improves the precision in the determination of the flux surface volumes and leads to increased accuracy in the calculation of the power deposition and current drive. The novel analytic results have been included and tested successfully in the volume computation routine of the ray-tracing code CODERAY, whereas the incorporation in the pWKB beam tracing code TORBEAM is under development (in cooperation with IPP-Garching, FOM).

(ii) The full-wave computation of the propagation of electron-cyclotron waves in the presence of edge plasma electron density fluctuations, using the FDTD code FWTOR, is under evolution. In the simulations, the electron plasma response is formulated in terms of the linear cold-plasma dielectric tensor assuming axisymmetric magnetic equilibrium, whereas the distribution of the density blobs along the propagation path is set-up with a uniformly random position and magnitude generator. The consistent interpretation of the results requires the computation of the wave spectrum in wave-vector space prior and after crossing the blob region, which is currently ongoing (in cooperation with Uni. Bayreuth, PSFC-MIT).

(iii) The model has been extended to the case of smooth density inhomogeneities of scales greater than the involved wavelengths (EC). The case has been formulated analytically and it is a valid approximation for ITER conditions. The numerical implementation of the analytical results will be implemented in the framework of the Work Plan 2012-2013. Two relevant Master’s Theses have been completed (see the details of the formulation in **Annex 3** and details on the Master’s Theses in **Annex 4**) (in cooperation with PSFC-MIT).

### 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 2011 (in co-operation with CEA, ULB, ITM2, U. Torino):

(i) A major part of the implementation of the plasma model using the tensor library OpenFOAM, as well as, the study of the particular characteristics of the setup of the sawtooth instability has been completed. The anisotropic thermal conductivity model for the plasma heat transport is considered. Preliminary numerical simulations revealed several computational issues that needed to be resolved before extending the OpenFOAM based numerical model with the full non-linear compressible MHD equations for plasma and validating it with the sawtooth instability. For this reason, extensive tests of our MHD implementation were conducted for simpler flow cases, such as the Orszag-Tang vortex revealing a good performance. For example, the incompressible version of our numerical model illustrates the emergence of the small scale structures through magnetic reconnection and current sheet formation in agreement with other published studies (see **Annex 5**). This also reveals that, for Mach numbers less than unity, the compressible turbulent plasma flows can be represented in terms of incompressible MHD equations.

### 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. It is still, however, a very challenging issue, given the remaining difficulties of the detection and the localisation of the magnetic islands. In a recent study (and proposal/application to the European Commission for preferential support) by IPP-Garching, magnetic islands could be more efficiently targeted by the EC beam, if they are locked to Resonant Magnetic Perturbations (RMPs) generated by internal (in-vessel) coils. The major milestones of the 2011 work-programme had been the observation of improved stability of tearing modes by early application of ECCD, and the suppression of locked tearing modes (to Resonant Magnetic Perturbations generated by the first set of the in-vessel coils) by CW-ECCD at the “O”-point, in ASDEX Upgrade.

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

(i) Given the significance and the highest priority of the second milestone, no experimental time was allocated to the observation of improved stability of tearing modes by early application of ECCD. This milestone has received minor attention, in favour of the more significant second milestone.

(ii) A dedicated experimental program was conducted in ASDEX Upgrade, to lock tearing modes to Resonant Magnetic Perturbations (RMPs) generated by the first set of the internal coils. The experiments were focused on the m/n=2/1 island, which is closer to the boundary and, hence, the intensity of the RMP is larger. Although deceleration of the rotation of the mode had been observed, the desirable mode-locking was not achieved. This is a very challenging investigation which has never been conducted before. It is planned to continue the experimental program with the second set of the internal coils, which will be available in ASDEX Upgrade for the experimental campaign of 2012.

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

### 1.4.1.MD modelling of mixed layer erosion processes

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, examined only the case when incident particles bombarded a virgin Be crystal. In 2011, our objective was to perform cumulative MD experiments, where the effect of previous collisions is taken into account.

Work performed in year 2011:

(i) The MD code (LAMMPS) has been completely adjusted and tested to perform cumulative sputtering. The code is flexible and allows us to vary many physical parameters, such as relaxations of borders, incident energy and temperature, orientation of the crystal, etc. We have carried out lengthy time series of MD experiments for He atoms successively bombarding Be surface. Several MD samples were simulated, starting from few thousands up to few dozen thousands of Be atoms.

(ii) The Ueda pairwise He-Be and Be-Be interaction potentials (see **Annex 6**) have been applied. As to Tersoff triple potential, it turns out that the simulation time is incredibly large for the MD cell of necessary size (few dozen thousands of Be atoms), while for the small MD size, the side boundaries are too close to each other and results are artificial. So, simulations with the Tersoff potential have not been performed.

(iii) Boundary conditions were varied in order to exclude artificial numerical boundary effects. In the course of the MD experiments, we recorded and then analysed instantaneous atomic configurations in the bombarded Be crystal. Altogether a few dozen Be samples have been systematically studied for broad energy and temperature ranges, as well as for various combinations of boundary conditions (see accompanying **Annex 6**). The main achievement of the study is that we have numerically demonstrated the mechanism of destruction for a Be surface that is successively bombarded by He atoms. The details of the scenario of swift Be destruction are presented in the **Annex 6**. It is important to stress that this scenario is general and it is independent of the size of MD cell and the imposed boundary conditions.

(iv) As a result of simulations carried out in the later part of 2011, it became clear that subsurface defects play a much less important role than boundary conditions. The crucial factors are the incident dose and the spanwise size of the MD cell. The swift erosion (or macroscopic destruction) of the Be crystal takes place when the size of the gas bubble accumulated inside the substrate reaches the boundaries of the MD cell. In reality, however, instead of the borders there might be inner crystal imperfections that can drastically increase Be erosion rate.

## 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 charged particles of the plasma (electrons/ions), in wave launching and plasma geometry relevant to heating, current drive and mm-wave diagnostics applications. Our research is conducted in several parts:

The first part focuses on the role of wave effects (beam localisation, 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 in tokamak plasmas. 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 breaks down and a more sophisticated treatment is called for.

In the second part, the main objective is to clarify the advantages of the utilisation of non-Gaussian beams in EC wave systems. This is formulated in terms of modelling the wave evolution 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 idealised 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 2011 (in cooperation with the Institutes indicated):

(i) In this period, the development and improvement of the full-wave code FWTOR, based on the FDTD method, was continued. After a series of corrections regarding the incident wave source and the boundary conditions, the code was benchmarked in the cases of 2D plane-wave and Gaussian beam propagation in vacuum and cold plasma, using different implementations of the boundary condition scheme (Mur ORBC and CPML ABC). The results expected for the 1D simulation were recovered successfully in all cases, whereas in 2D the improper nature of the plane-wave source as an incident field was pointed out. During the benchmarks, it was seen in practice that the code, even for simple plasma geometries, has very high computing-resource demands for parameters relevant to modern tokamaks. In this respect, an effort for the parallelisation of the code has been set forward, starting from a feasibility study on the different techniques of parallelisation of the FDTD algorithm. 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 parallelisation strategy should be followed and the relevant effort was initiated. (In cooperation with PSFC-MIT.)

(ii) 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. In the last period, the numerical computation of the beam amplitude was incorporated in NGBT and its coupling to TORBEAM was finalised, however the validity of the computation was not rigorously proven. The computation of the amplitude profile was amended in terms of novel analytical results, indicating that the amplitude transport equation is not invariant under a transformation H΄=f⋅H of the ray Hamiltonian. The proper form of the amplitude transport equation was derived, its solution algorithm was incorporated again in the routine NGBT and the numerical benchmarks for all the test cases were completed successfully. The study moved to amending the existing theory for the computation of the phase-shift with the necessary updates, incorporating the new theoretical findings in the numerical codes and benchmarking the results, which is currently in progress. (In cooperation with IPP-Garching.)

(iii) Test runs on the FP equation with the time dependent diffusion tensor were performed in order to access the conditions for numerical stability. The effort will continue also through 2012. The FP model has also been extended with the inclusion of pitch-angle collisional scattering. The extension involves two actions and thus, it constitutes a 2D extended model. Extension to a 3D model (involving 3 actions) and parallelisation is under way and it is expected to be boosted through an envisioned (discussions are under way) collaboration with CEA (see **Annex 7** and [Y. Kominis, K. Hizanidis, A. Papadopoulos, A. K. Ram, “Quasilinear kinetic formulation of wave-particle interactions in RF heating and current drive”, 53rd Annual Meeting of the APS Division of Plasma Physics, November 14-18, 2011, Salt Lake City, Utah, USA]) (in cooperation with PSFC-MIT, CEA).

(iv) Studies of the effect of EC on the highly energetic alpha particles (mono-energetic) and the respective time dependent RF diffusion tensor of our kinetic approach continue. The main differences in the approach is the action of RF on the time involving electron distribution (in resonance with EC) and additionally on the Maxwellian distribution of the fusing ions as well as with the generated high energy mono-energetic alphas. On the other hand, the difficulty encountered is the collisional scattering of the involving (as far as their distribution is concerned) electrons on the ions and the synergetic effects that this RF induced electron distribution evolution may have in the overall process (in cooperation with PSFC-MIT).

(v) Test runs on the FP equation with a simplified diffusion tensor suitable for LH range of frequencies were performed and the associated driven current was numerically evaluated. In order to estimate deviations from the quasilinear result, the FP model has also been extended with the inclusion of pitch-angle collisional scattering. The extension involves two actions and thus, it constitutes a 2D extended model. Preliminary results indicate higher driven current prediction in the low collisionality regime. Extension to a 3D model (involving 3 actions) and parallelisation is under way and it is expected to be boosted through an envisioned (discussions are under way) collaboration with CEA (for details see **Annex 7**) (in cooperation with PSFC-MIT, CEA).

(vi) The applicability of previously obtained analytical mappings for the change of particle positions and momenta during single interactions with localised wave packets has been considered. It has been concluded that the utilisation of these analytical mappings significantly reduces the computational time for calculations of averaged particle position and momentum variation for the case of multiple transits of an ensemble of particles through a localised wave-packet, as it is the case on a tokamak device. A nonrelativistic mapping has been constructed and preliminary studies for the construction of a relativistic mapping have been performed. The latter will extend the applicability of our approach to wider range of particle (electrons, protons, deuterons and alphas) velocities. The possibility of extending our studies from the case of constant wave-packets to the case of evolving wave-packets with varying amplitudes, spatial width, and phase/group velocities has been identified (in cooperation with PSFC-MIT).

(vii) The numerical code (in Matlab), which evaluates the fields (electric and magnetic) as well as the Poynting vector inside the blob in the vicinity and few blob radii away from the blob, has been completed and it is now ready to be benchmarked against a similar one in FORTRAN which is independently under development at PSFC/MIT. Preliminary results show considerable differences from the simplified geometric optics approach. The latter is incapable of capturing all the interesting new features the diffraction brings into the picture (see **Annex 8**) (in cooperation with PSFC-MIT).

(viii) The work is continuing and results have been obtained [on the basis of the numerical code developed in task (vi)] for choices of blob sizes, inhomogeneity magnitudes and angles of propagation with respect to the magnetic field. The work concerns a single spectral component (k vector) as yet. Extension from the single plain wave component to spatially confined EC beams is under way in order to prepare the grounds for introducing the statistics as far as the blob sizes and positions are concerned. The work will continue in 2012 according to the Work Plan of 2012-2013 (in cooperation with PSFC-MIT).

## 1.7 Theory and modelling for ITER

### 1.7.1.Application of computational tools for 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. The analysis is entirely done using a parallel, unstructured, Navier-Stokes solver developed at University of Cyprus.

Work performed in year 2011:

(i) In early 2011, code modifications were implemented that improved the stability of the numerical code for computation at high Ha numbers.

(ii) After completing the required code modifications for improved stability, the first part of 2011 was devoted to Direct Numerical Simulations of pipe turbulence entering an increasing magnetic field. Both poorly-conducting and insulating wall cases were considered, while a parametric study at different Ha numbers (50 < Ha < 7000) was also completed.

(iii) The implementation of wall-functions for simulations at Ha > 10000 was initiated in the second half of 2011 with the testing of functions in simple geometries. However, numerical tests in more complicated three-dimensional configurations showed the performance of the wall functions not to be satisfactory. Based on the work completed, it has been concluded that the implementation of improved wall functions would require a much more intensive effort, which has not been foreseen in the work plan for 2012. For this reason, this activity will not be continued in 2012.

### 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 2011 (in co-operation with the Institutes indicated): ULB, KIT

(i) The various properties of the statistical subgrid-scale LES model have been assessed. The statistically stationary Kolmogorov flow in which the flow arises when an incompressible fluid is subjected to an artificial, sinusoidal body force was chosen as a test case, in order to investigate the feasibility of the new LES model. Numerical experiments have been conducted for cases without and with strong external magnetic fields. For small computational domains, the translational invariance in the streamwise direction can be broken, resulting in a situation in which the turbulence statistics depend on the computational box aspect ratio. It was found that, this situation can influence the performance of the subgrid-scale turbulence model. Finally, another important property also found was the possibility of statistically erasing of the dominant turbulent features due to the ensemble-averaging process (see **Annex 9**).

(ii) The 2D mixed convection in a horizontal channel with sinusoidal temperature on the lower wall has been studied using the full MHD equations and the low magnetic Reynolds (Rm) approximation. The flow regimes have been investigated regarding the effect of Reynolds, Hartmann, Rayleigh, and Prandtl numbers. It was found that the velocity and temperature are damped by the magnetic field. The wall-normal magnetic field modifies more efficiently the flow and thermal features with respect to the streamwise magnetic field, while the spanwise magnetic field has no impact on the present 2D flow configuration. For high values of Rayleigh, the wall-normal magnetic field affects both recirculations that appear near the lower and upper walls of the horizontal channel, while the streamwise magnetic field seems to influence mostly the recirculation near the upper wall. In these cases, the penetration of the hot plume is stronger covering a larger part of the channel. For low Prandtl fluids, buoyancy becomes stronger and the effect of the magnetic field is more pronounced with respect to the high Prandtl fluids. Finally, it was found that, as the buoyancy becomes stronger, the numerical solution based on the low R_{m} approximation differs from that of the full MHD equations for R_{m}=0.001 or 0.01, which are typical expected values in several duct configurations and devices operating at ITER conditions (see **Annex 10**). Also, the effect of wall electric conductivity on 3D MHD Rayleigh-Benard convection in a duct has been assessed using the thin wall approximation. The thin wall boundary condition is employed at the four walls of the duct following the finite element methodology and performing integration by parts in order to reduce the order of integration. Thus the boundary conditions on the electric potential assume the form of natural BC’s in the finite element formulation. In this fashion previous experimental observations pertaining to the onset of Rayleigh-Bénard convection in the presence of a strong magnetic field, reported by Burr and Müller (JFM 2002) for highly conducting side walls and moderately conducting Hartmann walls, were recovered (see **Annex 11**). The three dimensional stability of two counter rotating vortices in the presence of a magnetic field was examined and it was seen that it is due to the elliptic structure of the streamlines. Over a range of dimensionless heat production rates and Hartmann numbers the dominant eigenmode is a travelling wave mode. The importance of the local Ekman number and flow eccentricity are identified as key factors in determining the stability threshold (see **Annex 11**).

(iii) The effect of curvature on the MHD stability and transition has been studied in a closed toroidal electro-magnetically driven duct. Electrically conducting walls were considered by implementing proper conductive boundary conditions models. In order to verify the properties of different flow regimes, various DNS simulations have been performed to assess the effects of the Hartmann and Reynolds numbers. The parametric study included values of the Reynolds number in the range of 500 < Re < 10^{6} and particularly high Hartmann numbers 3 < Ha < 5000. Two kinds of Taylor instabilities were found for the present electromagnetically driven toroidal duct flow with a critical value of Ha = 18 characterising each instability. For Ha < 18, the instability occurs as the double Taylor vortex pair is getting unstable, while for Ha > 18, the instability is connected with the appearance of numerous Taylor vortices generated near the sidewalls, especially the concave one. As the fluid inertia becomes dominant wiht increasing Reynolds number, each flow regime is destabilised differently. At transition and as Re → ∞ only localised disturbances near the concave sidewall layers were found for the present small toroidal duct (see **Annex 12**).

(iv) Rayleigh-Benard convection in the presence of a magnetic field, which is dominated by a standing wave followed by a travelling wave instability has been studied. The numerical method was adapted for the particular flow arrangement involving two counter rotating vortices. It emerges in the context of the two dimensional analysis of the flow. A parametric study of the range of validity of quasi 2D analysis was performed, in the context of Rayleigh Benard flows at high Ha. As it turns out Gr_{Cr} scales like Ha^{2} corresponding to a balance between buoyancy and magnetic forces. The secondary time dependent instability that was registered in the experiments by Burr and Muller (JFM 2002) is probably due to the elliptic instability of two counter rotating vortices that was above studied (see **Annex 11**).

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

Background and Objectives: Computational investigation of novel concepts pertaining to blanket design (DEMO included). 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 2011 (in co-operation with ULB, KIT):

(i) An immerse boundary method has been developed to account for the effects of electric conductivity and roughness of walls based on appropriate electric boundary conditions, as an alternative approach to the widely used thin wall models. This technique has been successfully tested in several pure hydrodynamic flows and it has been extended properly to simulate MHD flows (see **Annex 13**). A parametric study has been initiated to investigate the effect of wall roughness in MHD flows relevant to blanket modules at relatively high Hartmann values, in order to clarify the limits of feasibility of several thin wall models that have also been developed. Due to limited computational resources, the direct numerical simulations will also be carried out during 2012. Also, the development of a parallelised code for 2D MHD stability in a square cavity using ScaLapack and ParPack with distributed Finite Elements has been completed. Domain decomposition was performed along the horizontal direction of the cross-section, in order to alleviate the storage requirements of the solution. Literature review was performed and the methodology for solving 1D problems containing strong layers was identified, via the finite element method coupled with a posteriori error estimates. The approach adopted is that proposed by Oden and Ainsworth (Comp. Meth. Appl. Mech. and Eng 1994) based on element residual methods. Several 1D and 2D test problems were constructed to be carried out during 2012.

(ii) Testing of the parallel code for high Hartmann numbers in comparison with experiments for Rayleigh Benard flow in ducts with conducting walls has been performed (see **Annex 11**). Extensive verifications of the developed numerical models have been conducted in several 2D and 3D test flow cases during 2011. Because of the extremely large demands in computational grids and CPU time, the 3D direct numerical simulations will be continued in 2012.

Last Updated (Friday, 11 January 2013 12:05)