1.2 Energy and particle confinement/transport

1.2.1. Turbulence and transport phenomena

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

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

1. Hamiltonian perturbation methods have been applied in a perturbed guiding-centre Hamiltonian describing particle motion in an axisymmetric magnetic equilibrium under the presence of stochastic static and wave fields. A hierarchy of evolution equations for a partially averaged particle distribution function has been derived.
2. The hierarchy of evolution equations derived in the first part of year 2009 has been considered for the investigation of the interplay between the intrinsic complexity of the particle motion and the external stochasticity due to random static and/or wave fields. The last member of the hierarchy describes the collective particle behaviour in the action space. The three action variables correspond to radial position, magnetic moment and parallel velocity. The equation describing the evolution of the particle distribution function in action space is a Fokker-Planck equation. The respective diffusion operator has a constant positive part describing irreversible dynamics related to external stochasticity due to random static and/or wave fields and a time-dependent part describing the intrinsic complexity of the particle motion under non-axisymmetric deterministic perturbations.

1.2.2 Particle transport model

Background and Objectives: It is well known that the radial cross field plasma transport through the edge (arising mainly by intermittent events in the Scrape Off Layer) which is dominated by turbulence, strongly influences the location and strength of the heat and particle flux to the wall, as well as, the processes of recycling, impurity influx and He ash removal. The modelling of these processes is performed by using complicated numerical codes that are based on the solution of Branginskii’s fluid equations (some models are using the more extended Balescu treatment) for the characterisation of the collisional plasma of Scrape Of Layer (SOL), coupled with the neutral transport (generally is a kinetic description) close to the divertor and the walls. The code SOLPS is going to be used extensively in the area of edge plasma modelling and mainly to investigate the methodology that must be followed so that the Dispersed Phase Structure Dimensionality (DPSD) concept, originally developed for hydrodynamic turbulence, can be incorporated within the aforementioned multiphysics codes used for the simulation of plasma and edge turbulence.

Work performed in year 2009:

1. The existing local C++ solver has been reorganised and the conservation equations for momentum and energy for the charged particles have been incorporated and the implementation validated. The characterisation of the collisional plasma in the Scrape off Layer of a Tokamak device which is given mainly by the solution of Branginskii equations is going to be performed via the SOLPS code. The code uses a fluid description (B2 part) for the plasma area far from the walls and a kinetic one (EIRENE part) in the plasma-surface interface sheath. Access to the aforementioned code has been obtained and so that future development within the edge plasma modelling framework will be incorporated in SOLPS. In parallel with the use of SOLPS we implement a multidimensional cold plasma numerical model that could describe the phenomena in Scrape off Layer of a Tokamak. This includes the solution of the Poisson and continuity equation for the charged particle dynamics within the plasma. The model takes also into account the phenomenon of photoionisation by solving the Helmholtz equations, describing in this way self-consistently the whole phase of formation and evolution of a gas discharge. The neutral gas particle dynamics given by the solution of the compressible Navier-Stokes equation are intended to extend in a three dimensional framework and couple with the gas discharge model. (For more details see Annex 1.)

1.3 MHD stability and plasma control

1.3.1. ECRH for MHD control

Background and Objectives: This activity focuses on the improvement of the theoretical understanding and the interpretation of experiments on MHD stabilisation with the use of localised EC heating and current drive.

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. One of the main aims is to construct a self-consistent numerical model for NTM dynamics and stabilisation with ECRH/ECCD. The primary attempt includes the coupling of different kinds of wave propagation solvers (ray tracing, beam tracing, and FDTD) with the modified Rutherford equation, as well as more sophisticated MHD/kinetic models for the description of the NTM dynamics, in realistic magnetic geometry including the non-axisymmetric perturbations (in the form of magnetic islands) generated by the NTMs.

An additional objective of the present project has been the re-evaluation of the suppression rate of rotating magnetic islands by ECCD, using more advanced theoretical models, since the previously conducted theoretical evaluation of the suppression rate of rotating islands, provided unrealistically high suppression rates, which were not compatible with the experimental findings.

Another part of the activity is devoted to the investigation of the effect of edge density fluctuations in ECRH/ECCD applications. In the experiments for ECCD-based NTM stabilisation in ITER, the beam will eventually cross the plasma edge where turbulence can be rather strong, so that the associated density gradient scatters the wave. One could expect this deflection to be unimportant; however, since the distance to the EC resonance is large, even a small scattering angle may lead to considerable beam broadening and therefore degrade the stabilising effect. Our study includes both asymptotic and full-wave propagation in a plasma with different types of distributions for the density fluctuations.

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

1. A module for the solution of the modified Rutherford equation was developed, including all the terms widely considered to be important, and its coupling to the ray-tracing code CODERAY was initiated. At the same time, in order to assess the physics basis and the numerical methods related to MHD equilibria, instabilities and their dynamics, as well as in order to benchmark our results, an evaluation of the available codes in this field was performed. The codes of interest that were selected are ESC (Grad-Shafranov solver) and NIMROD (MHD dynamics). Based on the know-how acquired, the construction of a closed MHD model for the description of NTM stabilisation dynamics, including the effect of the ECCD on the background equilibrium, has been ini8tiated. (In co-operation with IPP-Garching and FOM.)
   In addition, the suppression rate of magnetic islands for power absorption at a constant helical angle, which (under realistic experimental conditions) is feasible only for the non-rotating islands, was evaluated using the precise configuration of the current density (inside and outside the island), which was obtained numerically by new more sophisticated techniques. The suppression rate of rotating islands (by continuous ECCD), subsequently, was obtained by the average suppression rate for power absorption at all helical angles. It was found that the suppression rate of rotating islands is always negative (stabilising) for all island widths, but strongly peaked for small islands (i.e. very small for large islands) as we observe at the experiment (for more details, see Annex 2). The effect of the driven current on the classical equilibrium, which was not included in the present work, will be investigated in the future. (In co-operation with IPP-Garching.)
2. An estimate of the effect of edge turbulence on the EC wave deposition in ITER has been obtained on the basis of geometric optics. Assuming that the wave scattering by the edge density fluctuations is a random process, a Fokker-Planck model for the distribution of the wave scattering angles has been developed. The results show that the effect can be important and might cause problems in the use of EC beams for NTM control. For example, in the case of ITER the results imply an effect of additional beam broadening even up to 100%. Comparisons of the simple Fokker-Planck model with ray-tracing numerical calculations, based on the statistics of ensembles of rays in a plasma environment with edge turbulence, showed reasonable agreement (see Annex 3). Full-wave computations will be initiated as soon as the code FWTOR (former ECFW) is updated to use a different scheme for the outgoing-wave boundary conditions (for details see 1.5.1(i)). (In co-operation with U. Warwick.)

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.

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

1. We focused on the development of the full non-linear and compressible MHD version of the code SFELES, which will be used for the study of plasma instabilities and turbulence in toroidal geometries of general cross-section. A set of conservation equations for plasma mass, momentum, temperature and pressure must be simultaneously solved using the mixed spectral/P1-P1 finite element formulation. This set of equations has been already coded and checked. The P1/P1 formulation of the finite elements is very convenient for the study of turbulent flows because it results in coarser linear systems and for this reason it has been adopted in the majority of the existing plasma codes. However, it introduces numerical instabilities and, thus, numerical stabilisation (e.g. the PSPG–pressure stabilised Petrov Galerkin) is needed. The appropriate stabilisation terms for the numerical schemes needed for the code SFELES were coded and tested. However, the numerical stability properties of the code were not satisfactory. For this reason a different public domain code named “Openfoam” is used for the study of plasma.

1.3.3 Effects of rotation on stability of multi-phase MHD turbulence

Background and Objectives: The transport of scalars and of a dispersed phase in MHD turbulence is relevant to both liquid metal flows found in the blanket modules of tokomaks and in fusion plasmas. Coherent structures in turbulence can potentially modify the transport of (passive) scalars as well as contribute to the preferential concentration of transported particles. This activity aims at developing the computation capability that would allow studying the effects of mean shear and frame rotation on particle and scalar transport with the aim to develop a modelling phenomenology.

Work performed in 2009:

1. A new parallel 3D MHD shear code with scalar transport and Lagrangian particle tracking was implemented. The code was used to study passive scalar transport in sheared MHD turbulence, initially in stationary frames and later in the presence of system rotation. Simulation results are in excellent agreement with existing literature for the case of scalar transport in the absence of system rotation. In the case of MHD turbulence sheared in a rotating frame, where previous results were not available, it was found that the system rotation has a strong effect on passive scalar transport because it tends to modify the turbulence structure. By carrying out simulations covering a wide range of parameter values, it was shown that for a given magnetic Reynolds number and magnetic interaction number, passive scalar transport depends strongly on the ratio of the mean shear time scale to the magnetic time scale and the ratio of the frame rotation rate to the mean shear rate. The idea of a scalar structure dimensionality tensor was introduced and it was shown to provide an effective statistical descriptor of anisotropy in passive scalar structure in MHD turbulence.

1.5 Physics of plasma heating and current drive

1.5.1. ECRH-ECCD and transport in fusion plasmas

Background and Objectives: In this activity we investigate the Electron Cyclotron (EC) wave propagation, absorption and current drive, as well as the wave-electron resonant interaction, in wave and plasma geometry relevant to ECRH/ECCD applications in plasma heating, mm-diagnostics and MHD stabilisation. The activity has several parts:

The first part focuses on the importance of wave-beam effects and hot plasma effects in the propagation and absorption of EC waves in tokamak plasmas. In the field of wave propagation in plasmas, the mainstream in theory and applications is oriented towards frequency-domain asymptotic methods, assuming linear plasma response, where the solution of the wave-plasma system is more tractable. However, in many cases of interest (like e.g. O-X-B mode conversion, high-power ECRH) this approach breaks down, the solution becomes questionable and therefore a more sophisticated treatment (full-wave and nonlinear) is called for.

In the second part, we analyse the propagation of EC waves, their absorption, the current drive and the transport of particles when the NTM instability is excited inside a tokamak. The general idea for the stabilisation of NTMs is to drive a current in the vicinity of the magnetic island’s O-point. So far, the analysis of the EC propagation and absorption has been done in axisymmetric magnetic topology, hence ignoring the presence of islands. We study the propagation in the presence of magnetic islands, focusing on the effect of the island topology on the efficiency of the EC absorption and ECCD.

In the third part, the main objective is to address the advantages of the utilisation of non-Gaussian beams in ECRH systems. This is performed in terms of modelling the ECRH/ECCD in simplified and realistic tokamak geometry, using the paraxial WKB beam tracing method for following the wave propagation and the linear hot plasma dielectric tensor for the plasma response. The results are exploited in order to upgrade certain features of the existing TORBEAM code regarding the propagation and absorption of arbitrary-shaped EC beams in realistic magnetic configuration.

Finally, in the fourth 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.

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

1. The implementation of the full-wave code FWTOR (former ECFW) in 2D was completed and benchmarked for simple cases of wave propagation. The benchmarks showed that, for accuracy reasons, there is a need to use a different scheme for the outgoing-wave boundary conditions, and a related improvement was initiated. The numerical code uses a realistic magnetic equilibrium, based on a Hamiltonian field-line tracing model, and assumes a linear, weakly-relativistic hot-plasma response to EC waves. Regarding the implementation of quasilinear and nonlinear models for the plasma response, a feasibility study for the cases of interest showed that Fokker-Planck codes or simplified nonlinear models (e.g. macro-particle approximation) can be coupled to the FDTD code and run for reasonable CPU times only in parallel implementation. A feasibility study of ECRH/ICRH and mm reflectometry/interferometry scenarios of interest for ITER was performed, and the resulting scenarios will be simulated after the code FWTOR is updated.
2. The calculation of the volume of flux surfaces corresponding to magnetic islands was performed as the first step towards the inclusion of a non-axisymmetric equilibrium in the codes CODERAY and TORBEAM. For the specific purpose, a FORTRAN routine was developed. The input of the routine are the parameters which characterise a magnetic island, and the output is the volume of a flux tube shell around the O-point of the island. The routine was included in the mentioned codes in order to calculate the power density profile when ECRH is deposited in the area around the island. A thorough study of the effect of the magnetic island geometry on the ECRH/ECCD deposition was performed (see Annex 4). In addition, the case of stochastic magnetic fields generated by the overlapping of the magnetic islands of two NTMs has been studied. It was shown that the wave propagation is not affected significantly by the complex structure of the magnetic field, while, regarding the deposition, one has to resolve the following issues (currently under investigation): (a) Which are the exact positions of the islands’ O-points? (b) What is the proper definition for the deposition volume? (c) Is there a mechanism for the simultaneous stabilisation of two NTMs? Finally, the study of the perturbed version of axisymmetric equilibria based on analytical solutions of the Grad-Shafranov equation has been initiated. (In co-operation with IPP-Garching.)
3. The coupling of the FORTRAN routine NGBT with the TORBEAM code, which had been completed in the previous year, was upgraded to include the case of experimental input for the initial EC beam amplitude profile. Then, the influence of non-Gaussian beams on the EC power deposition profiles was investigated, focusing on particular issues of interest for ASDEX Upgrade and ITER. It was found that the influence of spurious higher-order modes on a beam designed to be purely Gaussian can in some cases be significant, but also that improvement in the localisation of the absorption by convenient superposition of Gaussian and higher-order modes is feasible (see Annex 5). (In co-operation with IPP-Garching.)
4. Wave-particle interaction with localised wave packets and their result on ECRH-ECCD and transport have been investigated for waves obliquely propagating with respect to the confining magnetic field. The canonical perturbation theory has been applied in order to calculate approximate invariants of the motion. The latter describe all the essential features of the phase space of the particle motion. Moreover, they provide useful information for the collective particle behaviour as described by their distribution function or averaged quantities. The dependence of the interaction on the form and width of the wave packets and a significant role of finite Larmor radius effects are incorporated in the respective analytical results. (In co-operation with PSFC-MIT, USA.)
5. The application of the diffusion equation with non-singular and time-dependent diffusion tensor, derived in 2008, has been considered for describing ECRH, ECCD and transport. It has been shown that this equation is more appropriate than the standard quasilinear equation for describing cases where the rf fields are of small to moderate amplitude, which is actually the case relevant to ITER and other fusion devices. However, the respective diffusion tensors are not positive definite in the entire action and time domain, resulting in a challenging numerical problem. Preliminary numerical results have been obtained for a characteristic one-dimensional case and significant differences with respect to the standard quasilinear theory have been shown. For more details see Annex 6. (In co-operation with PSFC-MIT, USA.)
6. The canonical perturbation method with the utilisation of Lie series has been applied to the case of particle interactions with multiple waves with different frequencies and wave-vectors and approximate integrals of the single particle motion have been calculated. The results have been used in order to study collective particle dynamics as described by their distribution functions or averaged quantities. (In co-operation with PSFC-MIT, USA.)

1.7 Theory and modelling for ITER

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

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

Work performed in year 2009 (in co-operation with ULB):

1. A FORTRAN 90/95 module for modelling MHD flows was optimised (numerical stability) during 2009 (after an extensive validation during 2007 and 2008 for fully-developed turbulent HD and strong MHD flows).
2. The code was used to model flow in a fringing magnetic field with a maximum Hartmann number of 7000 in a pipe with conducting walls. The results of the simulation were compared against asymptotic approximation solutions and experimental results. The results match the asymptotic approximation results very well before the drop in magnetic field, which suggests that simplifications in the asymptotic approximation method detrimentally affect the downstream solution. The comparison with the analogous case in an insulating pipe (performed in 2008) showed that the over-speed areas are reduced by the conductivity of the wall, henceforth delaying the onset of turbulence downstream. Further details are given in Annex 7.
3. An analysis of one-dimensional spatially varying magnetic fields, commonly used for numerical simulations, was done. A numerical iterative solution was suggested to convert one-dimensional fits (of the real experimental magnetic field) into two-dimensional discretely consistent fields.
4. The generated magnetic field compared qualitatively well with other analytical methodologies to construct consistent magnetic fields and, therefore, the validation was carried out successfully. The effect of the inconsistency of the magnetic field revealed itself as an important actor in the historically under-predicted peak of the transverse pressure gradient for liquid-metal flow leaving a strong magnetic field. Further details are given in Annex 7.
5. An assessment of existing CFD codes for MHD flows is being carried out in connection with the subtask (vi) below and the next subtask 1.7.2 (i). The recently developed quasi-static MHD version of the code SFELES was used together with other spectral accurate codes to predict the 3D MHD flow in ducts subjected to high magnetic fields. The flow is not expected to be laminar because of the high magnetic field, but rather turbulent. In this stage, the codes are being validated against the Hunt and Hartmann laminar solutions under high magnetic fields and various wall electric conductivities. The work has been focused on the comparison of the laminar sidewall jets that are formed at various Reynolds and Hartmann numbers. The pressure driven flow in rectangular ducts under the effect of a uniform transverse magnetic field is a case where strong Joule dissipation will probably suppress motion. Thus, most studies have been focused on stationary inertialess approximations related to laminar flow. However, this assumption is not always true and experiments have shown that turbulent fluctuations are still present under very strong magnetic field conditions. In contrast to the usual channel flow theory, turbulence in this high Hartmann number flow is localised in the sidewall layers, parallel to the magnetic field. Prior to the transition to turbulence, wakes are formed according to the semi-analytical solution proposed by Hunt. The transition occurs for a certain combination of high Hartmann and Reynolds numbers according to existing experimental data. The laminar results have been compared to the semi-analytic solution of Hunt and the analytical solution of Hartmann for Hartmann numbers up to 10000. The heavy numerical simulations have been conducted at the Julich supercomputer centre.
6. Attempts to extend the range of Hartmann numbers are being made (see also Annex 8). A new implementation of fully-developed inflow boundary conditions for three-dimensional ducts has been implemented in order to more accurately transfer turbulent information to new computational domains, saving henceforth computational resources. The discretisation of the Lorentz force is the most challenging aspect in the Direct Numerical Simulation (DNS) of liquid-metal flow at high Hartmann numbers. The discretisation has been improved in order to enhance numerical stability during simulations of very high Ha MHD flows in three-dimensional duct flows. This development made possible the simulation of MHD duct flows at Ha of the order or 7000 to 10,000. (Also in co-operation with KIT.)
7. Preliminary tests carried out in early 2009 demonstrated the need for improved grid generation techniques for the simulation of very high Ha flows. This led to the implementation of a structured near-wall boundary layer with grid refinement close to the walls (able to reproduce sharp velocity gradients present in Hartmann and side layers), in combination with a triangle-based unstructured core away from the walls. This approach, in combination with an improved discretisation of the Lorentz force, resulted in improved numerical performance. The improved structured/unstructured grid design and code enhancements have been successfully tested for 3D duct flows for Ha numbers up to 10000 under electrically insulating and poorly-conducting walls. 

1.7.2 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 underlying physics. Due to many particular problems related to the magnetohydrodynamic flows, for example, MHD turbulence, heat transfer, transition and stability, a series of individual accurate solvers based on computational fluid dynamics techniques are being developed.

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

1. The 3D flow of conducting fluids in ducts subjected to high magnetic fields is being investigated using accurate Direct Numerical Simulation (DNS) and large-eddy simulation (LES) methods. For high magnetic fields, the flow is expected to be rather turbulent than laminar. The fully developed turbulence is confined in the sidewall layers where the fluid is not subjected to high resistance by Lorentz forces. This phenomenon is being studied together or separately from the conductivity of the walls. The recently developed quasi-static MHD version of the code SFELES is used together with other accurate codes. The work was focused on the study of the effect of the Reynolds and Hartmann numbers to the onset of side-layer MHD instability and the transition to turbulence using DNS. The simulations are being repeated by increasing the Hartmann number until the transition to turbulence occurs. The same approach is being followed for other Reynolds numbers and a complete mapping of the transition is being formed. Full visualisation of the flow at the moment of transition is obtained and the results are compared, where possible, to the experiments (see also Annex 8). (In co-operation with ULB, KIT.)

1. a) Direct numerical simulations (DNS) were performed to study the natural MHD convection in an annulus with the end-walls considered thermally adiabatic and the cylindrical surfaces isothermal and at different temperatures. Homogeneous volumetric fluid heating was also considered with an external homogeneous magnetic field in the direction normal to the cylinder axis. (For more details see Annex 9 and Annex 10.)
   b) Direct numerical simulations were performed of the turbulent-transient natural convection cooling of a hot low-Prandtl number fluid contained in a vertical cylinder in the presence of a vertical and constant magnetic field. The analysis focused on the effect of the magnetic field, expected to reduce the turbulent fluctuations, suppress the fluid motion and, consequently, decrease the cooling rates. It was observed that, in the initial stage of boundary layer development on the cylindrical wall, there is no influence of the magnetic field. In the next stage of thermal intrusion of cold fluid, a significant reduction of heat transfer is observed, while at the final stage of thermal stratification, the magnetic field decelerates or accelerates the cooling of the low-Prandtl number fluid for low or high Rayleigh numbers, respectively. This dependence of the magnetic field effect on the Rayleigh number is rather unusual and it was found to be related to the cold vortices emanating from the vertical boundary layer. The effect of a vertical magnetic field on natural convection cooling of a fluid in a cylindrical container was studied for the transient and laminar regime. In general, the flow is dominated by two major stages, i.e. the development of a boundary layer as the fluid descents at the cold wall and the subsequent stage of the cooling-down. The fluid motion activates forces that are connected with the magnetic field and as a consequence interfere with the heat transfer and the duration of cooling-down. Numerical simulations were performed for a range of Rayleigh and Hartmann numbers where laminar flow occurs and a range of Pr<1, similar to known electrically conductive fluids. The time for the boundary layer to reach steady-state, the Nusselt number on the cylindrical wall and the time for full cool-down were determined. The dominant heat transfer mechanisms, conduction and convection, in relation to the magnetic field were studied. (See also Annex 11.) (In co-operation with ULB and KIT.)
2. Three-dimensional stability analysis was carried out of the experiments on Rayleigh-Bernard convection in the presence of a magnetic field by Burr & Mueller (JFM 2002). A parametric study was performed for large Ha values, accounting for the effect of wall conductivity. In this context, progressive thinning of Hartman layers and appearance of quasi-2D structure were observed. For intermediate Ha values, the 3D stability analysis gives better predictions of the experimental observations than the asymptotic solutions of Burr & Muller.
3. The establishment of connectivity between different eigenmodes and the verification of critical conditions was pursued. The numerical code that performs the stability analysis via continuation on specific eigenvectors was extended in order to follow specific modes for the construction of neutral stability diagrams. The base flow and stability problems are solved sequentially. Tests on parallelisation of finite elements using GMRES on windows and LINUX systems were conducted. (See also Annex 12.) An extension of the above code in order to solve the base and stability problems has been implemented with the base flow and stability problems solved sequentially. Work on alternative finite element formulations, with respect to basis functions, was not performed because the primary researched is not in the programme anymore. (In co-operation with ULB.)

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

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

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

1. The work regarding the immersed boundary method was focused on the incorporation of a hydrodynamic version of the method to existing numerical tools. The numerical model developed was based on the one described by Fadlun et al. [J. Comput. Physics, 161, 2000]. Appropriate forcing schemes for current densities have been developed for arbitrarily shaped immersed surfaces. A velocity and pressure reconstruction algorithm has been implemented to calculate the forces exerted on immersed bodies. Finally, the immersed boundary method was successfully tested against existing experimental and numerical results based on either body-fitted or other immersed boundary techniques in several flows. (In co-operation with ULB.)
2. The verification and validation of the developed methodology has been completed for homogeneous and inhomogeneous magnetic fields against body-fitted DNS results and experimental studies. Using the proposed velocity and pressure reconstruction schemes the two-dimensional and three-dimensional MHD flows around a circular cylinder have been studied in detail (see Annex 13). These test cases were considered as a basic benchmark problem in order to verify the method for various intensities and directions of the magnetic field. They were chosen because they involve interesting phenomena such as flow unsteadiness, vertical motion and time varying separation, which are also expected to significantly modify heat transfer rates due to the obstruction.
3. The onset of unsteady flow regime, its dependence on critical Hartmann numbers and the variation of the force coefficients has been examined and compared against previously reported studies for a wide range of dynamic parameters. Detailed three-dimensional simulations for cylindrical obstructions inside rectangular ducts and channels have been initiated.