Project-Calvi:Application of Vlasov codes to magnetic fusion

Inria / Raweb 2009
Presentation of the Project Calvi

Section: New Results

Application of Vlasov codes to magnetic fusion

Participants : Pierre Bertrand, Nicolas Besse, Jean-Philippe Braeunig, Nicolas Crouseilles, Daniele Del Sarto, Sever Hirstoaga, Giovanni Manfredi, Michel Mehrenberger, Etienne Gravier, Guillaume Latu, Jean Roche, Amar Mokrani, Simon Labrunie, Hocine Sellama, Eric Sonnendrücker.

The computation of turbulent thermal diffusivities in fusion plasmas is of prime importance since the energy confinement is determined by these transport coefficients. Fusion plasmas are thus prompt to instabilities and required a self-consistent analysis of the plasmas and the electromagnetic fields. Since collisions between particles play a negligible role, the kinetic plasma behavior is described by the well-known Vlasov model. A key issue is the resonant interaction between particles and waves, which has to be accurately described. The self consistent coupling of the electromagnetic fields then exhibits resonant (Landau) interactions between charged particles and electromagnetic waves, a feature that cannot be addressed in the fluid limit.

Development of a 5D gyrokinetic code

The collaboration around the optimization of the GYSELA code used for gyrokinetic simulations of turbulence in tokamaks went on. The upgrade from four to five dimensions of the phase space is now efficient on several thousands of processors. Last year, several developments of the code enabled to achieve this task.

The aim of this work is to improve the accuracy and the physical relevance of the GYSELA 5D code. Mainly, we would like to perform simulations using a mesh based on curvilinear coordinates in such a way mesh lines are aligned with magnetic field lines. This will diminish numerical diffusion, reduce the size of the computations and improve results accuracy. Moreover, a conservative scheme for plasma turbulence simulations will be used. Using a conservative formalism is an asset for physical relevance of simulations and allows a directional splitting of the Vlasov equation, what is very convenient to deal with curvilinear meshes.

We started in this way by studying conservative numerical methods with curvilinear coordinates for the Guiding-Center (GC) model, which is a reduced 2D conservative model for plasma turbulence. We use the Parabolic Spline Method (PSM), which is a semi-lagrangian conservative scheme using a cubic spline reconstruction. It is a similar method to the Backward Semi-Lagrangian (BSL) scheme [88] (exactly the same scheme for linear constant advection), but based on the conservative form of the Vlasov equation, obviously for conservative models. Different curvilinear meshes are tested on linear advection benchmarks and on simulation of the growth of unstable turbulent modes. This work has been presented at Conference "Numerical Models for Controlled Fusion" and described in the Conference Proceeding Braeunig et al [39] .

Afterward, the conservative semi-lagrangian PSM scheme has been integrated in the GYSELA code. We first have faced numerical difficulties about the velocity field characteristics. We noticed that cells geometric evolution through the scheme should conserve their volumes, what is equivalent to have a divergence free velocity field at the discrete level, and not only at the continuous level. It is important to obtain a robust scheme in the sense of getting nearby a maximum principle condition. In the same way of improving robustness, we propose an "entropic flux limiter" for the PSM scheme to cut off the spurious oscillations that may appear when dealing with small scales turbulent structures in the flow, below the mesh size. A different way of MPI parallel computing is used for this scheme (transposition), achieved by Guillaume Latu. Drift kinetic 4D benchmarks have been performed and compared with the former BSL scheme. This work is described in an INRIA report Braeunig et al [47] . Further validation will be done, especially for the full 5D gyrokinetic model in the GYSELA code which is already functional with the PSM scheme.

Physicists have improved the Gysela code in order to add heat flux driven ion turbulence. It involves large displacements during the advection phase of the Semi-Lagrangian method. In order to take care of this new feature, another parallel algorithm is needed to replace the local spline method that requires a too strong CFL-like condition. A 4D transposition parallel operator has been found that provides a solution. Even if the communication time is growing with this new algorithm, it remains less than the advection computing time, even for big runs. The overall parallel relative efficiency of Gysela decreased in this new configuration, nevertheless the asymptotic behaviour (for large number of processors) is still very scalable.

Numerical study of the gyroaverage operator

In this work, we are concerned with numerical approximation of the gyroaverage operators arising in plasma physics to take into account the effects of the finite Larmor radius corrections. Several methods are proposed in the space configuration and compared to the reference spectral method. In particular, a new method is proposed; the basic point is the expansion of the function on a basis (such as polynomial basis). Computing the gyroaverage of a function then reduces to compute the gyroaverage of its basis, which are known analytically. Hence, in the same way as finite element formulation, the method can be formulated into a matrix form. The approach presents other several advantages. On the one side, the number of quadrature points is automatically determined as the intersection between the mesh and the Larmor circle, which provides the necessary adaptivity to reach a good accuracy even for large Larmor circles. Contrary to previous works, the basis function are not evaluated at the quadrature points but integrated on a circle arc, which turns out to be more accurate. the gyroaverage operator also includes an integration with respect to the Larmor radius which belongs to [0, + $\infty$ ] . Hence a good accuracy when large Larmor radius are considered is required. But, it turns out that it is not sufficient since the numerical integration with respect to the Larmor radius is also important in order to reach a good accuracy for the gyroaverage operator. For example, traditional quadrature rules (such as trapezoidal or Laguerre ones) are not sufficiently efficient. In this context, we develop a new approach which appears to be very efficient for arbitrary wave numbers. The coupling with the guiding-center model leads to very good results compared to the analytical results in the linear regime and to the spectral method for the gyroaverage operator. We then investigate the influence of the different approximations considering the coupling with some guiding-center models available in the literature.

Plasma-wall interactions – application to ELM modes on JET

The vast majority of plasmas produced in the laboratory are in contact with a material surface. In fusion devices, the surface can be either the material vessel that contains the plasma, or some ad-hoc device (limiter or divertor) specifically designed to optimize the interaction with the charged particles. These surfaces are eroded by ion and neutral bombardment and may thus see their lifetime considerably reduced. This is one of the key issues for the ITER tokamak, and even more so, for future commercial power plants.

The heat load on the divertor plates can be considerably strained during violent events known as ELMs (edge-localized modes). Although the precise origin of ELMs is still under debate, in practice, these events result in the ejection of energetic charged particles beyond the magnetic separatrix. These particles travel along the magnetic field lines and end up hitting the divertor plates.

In order to model such a scenario, we have developed a 1D Vlasov-Poisson code with open boundaries, for both ions and electrons. The spatial co-ordinate represents the distance along the magnetic field line, whereas the absorbing boundaries are supposed to mimic the divertor plates. Appropriate diagnostics are set up to compute the particles and heat fluxes on the plates.

The main numerical bottleneck is represented by the very large ratio between the system size ( $Im23 ${\#8764 10}$$ m) and the Debye length ( $Im24 ${\#8764 10^{-4}}$$ m). A possible approach is to use the recently developed “asymptotic preserving” schemes [72] , which allows to work with very small Debye length without incurring numerical instabilities. We are currently implementing this type of scheme in the Vlasov-Poisson code. This project is developed in collaboration with W. Fundamenski and S. Devaux (JET).

Gyro-water-bag approach

In the Gyro-Water-Bag (GWB) approach, a discrete distribution function taking the form of a multi-step like function is used in place of the continuous distribution function along the velocity direction. According to Liouville's phase-space conservation property the distribution function remains constant in time between the bag contours. The time evolution of the system is completely described by the knowledge of the contours. We get a set of hydrodynamic equations, where the system behaves as N fluids coupled together by the electromagnetic fields (in our case the quasi-neutrality). As a matter of fact a small bag number (not more than 10) has been shown to be sufficient to correctly describe the Ion Temperature Gradient (ITG) instability observed in fusion plasmas. Thus the water-bag offers an exact description of the plasma dynamics even with a small bag number, allowing more analytical studies and bringing the link between the hydrodynamic description and the full Vlasov one. Encouraging results have been obtained using three nonlinear waterbag codes: a 3D nonlinear self-consistent gyro-water-bag code which is based on a semi-Lagrangian method (this code was subject to comparison with the GYSELA code dealing with curvilinear meshes); a 3D quasilinear self-consistent gyro-water-bag code and a 1D nonlinear self-consistent electromagnetic-water-bag code where the Discontinuous Galerkin methods are used.

Trapped-ion driven turbulence: gyrokinetic modeling

Drift wave instabilities and their subsequent turbulence are believed to be an important component of the experimentally observed anomalous transport of heat in tokamaks. It was known for several decades, that the transport is “anomalous” in the sense that it does not follow the standard magnitude scaling of collisional transport. Expressed in terms of diffusivities the transport of ions and electrons is comparable, in the range of one or fewer times 1m²s^-1 in the core region, which is indeed one order of magnitude larger than the collisional ion transport and three higher than collisional electron transport.

Most previous studies of ion-temperature-gradient-driven drift wave turbulence have utilized a simple one-field, nonlinear fluid model, known as the Hasegawa-Mima (HM) equation. This equation is an important paradigm for the description of drift-wave turbulence in magnetically confined plasmas and has a structure similar to that of the Charney equation, describing geostrophic motions in planetary atmospheres. This model is indeed connected to two-dimensional (2D) turbulence. However it is well known that, for 2D turbulence, the feature of turbulence are different from the standard 3D turbulence.

At this step, different remarks must be pointed out:

(i) At high temperature and taking into account the toroidal effects, collisionless kinetic effects allow an effective short time decoupling of particles and field lines (“defreezing effects”) thereby introducing non fluid effects, as Landau resonances, strong finite Larmor radius effects and trapped particle resonances.

Trapped particles, in general have unfavorable magnetic drift, leading to the trapped particle instability, a localized low-frequency interchange. The most prominent type of long wavelength toroidal microinstability is Trapped Ion driven Mode (TIM) in the presence of a significant ion temperature gradient. The initial work by Kadomtsev and Pogutse 1971 indicates that the TIM mode has a real frequency below the diamagnetic drift frequency $\omega$ _* and below the bounce frequency $\omega$ _b . These instabilities are easier to excite in the core region of tokamak fusion plasmas where the plasma can be considered to be effectively collisionless. Thus these instabilities are characterized by frequencies of the order of the trapped particle precession frequency and radial scales of the order of several banana widths. These properties make TIM modes have a global nature i.e. large radial correlation lengths that can vary on the equilibrium scale length rather than on the ion gyro-radius scale length usually associated with kinetic micro-instabilities. It has been predicted that because of their large scales and in spite of their low frequencies, these instabilities could lead to a huge transport corresponding to a Bohm-like confinement observed in many experiments. Such Trapped-ion driven modes are thus a simple prototype of kinetic instability since they are driven through the resonant interaction between a wave and trapped ions via their precession motion and implies low toroidal numbers (n<100 ), allowing fast numerical gyrokinetic simulations (see Depret el al 2000).

Although TIM have been studied during several decades, there exists a real interest today. Recently trapped particle modes seem to play a key role in nonambipolar transport in ITER-relevant regimes, where strong enhancement of transport is thus predicted by bounce-harmonic resonance.

(ii) Another feature concerns the fact that turbulent transport can be greatly reduced in the presence of a magnetic shear. It is well known that the stabilization effect of the precession drift reversal for trapped ion driven instabilities is efficient in nonlinear regime. Magnetic shear describes the dependence of field line direction on the coordinate across closed flux surfaces. In a toroidal field with circular flux surfaces with minor and major radii r and R , the shear is defined as $Im25 ${s=\mfenced o=( c=) r/q\#8706 q/\#8706 r}$$ . The quantity $Im26 ${q\mfenced o=( c=) r}$$ is the winding number of magnetic field lines on any given flux surface, i.e. the ratio of the displacement in toroidal angle to displacement in poloidal angle between two points on a field line. Thus, in a toroidal plasma, the magnetic field is not uniform. Therefore, its magnitude varies both along and transverse to the field lines. The longitudinal variation traps particles (the projection of the guiding-center orbit onto a cross section at fixed toroidal angle traces out a “banana” shape). While the transverse variation, referred to as magnetic shear, has a strong stabilization effect on collective instabilities. This effect exists both in fluid and kinetic instabilities. For trapped particle drive instabilities (as TIM instabilities), it takes place through the reversal of the precession drift (see Kadomtsev and Pogutse 1971).

(iii) Although 5D gyrokinetic simulations are now possible with realistic geometry, such simulations require huge CPU time and memory resources. The most mature kinetic approach is the particle-in-cell (PIC) method. The method involves to integrate gyrokinetic Vlasov equation in time by advancing marker particles along a set of characteristics within the phase space. It basically involves a lagrangian formulation in which the dynamics of an ensemble of gyro-averaged “particles” are tracked.

Several previous points were addressed in a gyrokinetic model of ITG instabilities for TIM model: forward cascade processes, condensation process of turbulence toward low toroidal modes, turbulence stabilization for TIM modes due to reversed magnetic shear. Trapped-ion driven modes (TIMs) are studied by solving a Vlasov equation in action-angle variables using Hamiltonian formalism (see Depret et al 2000). This Vlasov equation is then averaged over the fast scales: cyclotron and bounce motion for trapped ions. The distribution is then calculated as a function of the poloidal flux $\psi$ (normalized to 2 $\pi$ and which indeed plays the role of a radial coordinate) and the precession angle of trapped ions $Im27 ${\#966 _3=\#981 -q\#952 }$$ , parametrized by the particle energy E and the trapping parameter $\kappa$ (depending of the adiabatic invariant $\mu$ ).

The reduced gyrokinetic Vlasov equation is given by (see Garbet et al 1990, 1992, 1994 and 2001):

$Im28 ${\mfrac {\#8706 f}{\#8706 t}+\#969 _d\mfenced o=( c=) \#954 ;s\mfrac ET_0\mfrac {\#8706 f}{\#8706 \# ... 66 _3}-\mfrac {\#8706 \mover J^_0\mover U\#732 }{\#8706 \#966 _3}\mfrac {\#8706 f}{\#8706 \#968 }=0}$$

(3)

where the trapped ion distribution function is $Im29 ${f=f_{\#954 ,E}\mfenced o=( c=) \#966 _3,\#968 ,t}$$ . The self-consistency constraint is

$Im30 ${C\mfenced o=( c=) U-\mfenced o=〈 c=〉 U_\#968 -C\#945 \#8711 ^2U=\mover n¯_i-n_{pi}}$$

(4)

where we have introduced $Im31 ${U=e_i\mover U\#732 /T_0}$$ , n_pi the passing ion density and the gyro-averaged pseudo-density of trapped ions is given by:

$Im32 ${\mover n¯_i\mfenced o=( c=) \#966 _3,\#968 ,t=\mfrac 2\sqrt \#960 \#8747 _{0}^\#8734 d{(E/T_0)}~\sq ... ax}d\#954 ~\#954 K\mfenced o=( c=) \#954 \mover J^_0f_{\#954 ,E}\mfenced o=( c=) \#966 _3,\#968 ,t.}$$

(5)

In the previous expression we have $Im33 ${C=\mfenced o=( c=) 1+T_i/T_e/f_p}$$ where $Im34 ${f_p=2\sqrt {2\#949 }/\#960 }$$ is the fraction of trapped ions. Consequently we have reduced the six-dimensional Vlasov equation into a set of 2D Vlasov equations (indeed $Im35 ${N_\#954 ×N_E}$$ equations) acting on a two-dimensional phase space $Im36 $\mfenced o=( c=) \#966 _3,\#968 $$ (precession angle $Im37 ${\#966 _3=\#981 -q_0\#952 }$$ , poloidal flux normalized to 2 $\pi$ ), parametrized by the kinetic energy E and the trapping parameter $\kappa$ ( $Im38 ${\#954 \#10877 \#954 _{max}}$$ with a chosen value of $\kappa$ _max = 0.975 for numerical simulations). $Im39 $\mover J^_0$$ is the gyroaverage operator - multiplication by $Im40 ${J_0\mfenced o=( c=) k_\#8869 \#961 _cJ_0\mfenced o=( c=) k_\#968 \mover \#948 ^_2}$$ - in Fourier space, $\rho$ _c being the gyro-radius and $Im41 $\mover \#948 ^_2$$ the banana width. The present implementation of this operator is based on a Padé approximation of the Bessel function J₀ . Eqs (3 ), (4 ) and (5 ) constitute the basis of our Vlasov model for describing TIM instabilities.

Here the density refers to the density of banana centers which do not correspond to the real density of particles. For instance Hahm and Tang 1996 use a Taylor expansion to give an estimation of the potential value averaged along a field line in the form $Im42 ${\mfenced o=〈 c=〉 \mover U\#732 _{banana~orbits}=\mover {U+\mfrac 12\#948 _{b}^2\#8711 ^2\mover U\#732 }\#732 }$$ , which allows to take into account the difference between the real trapped ion particle density and the density of banana centers.

The trapped ion driven mode has been studied by solving a Vlasov equation averaged over the cyclotron and bounce motion of trapped particles in the presence of a magnetic shear. The magnetic shear parameter appears to play a major role in the stabilization process of the instability and allows us to study the occurring of turbulence inside the system. The major mechanism of turbulent transport is large-scale (electrostatic) convection of trapped particles. The dominant driving process is the temperature gradient over some critical value, and one of the stabilizing elements is reversed magnetic shear. Turbulent plasma is obtained when the magnetic shear reaches a limit value which is positive. The introduction of the magnetic shear parameter allows us to study the different regimes of turbulent transport from plasma state described by a small number of toroidal modes till the occurring of Kolmogorov-like cascade which implies the excitation of a large band of toroidal modes (fully developed turbulence). This is achieved by using a semi-Lagrangian Vlasov code which allows a very accurate description of the trapped particle distribution function in phase space.

Full wave modeling of lower hybrid current drive in tokamaks

This work is performed in collaboration with Yves Peysson (DRFC, CEA Cadarache). The goal of this work is to develop a full wave method to describe the dynamics of lower hybrid current drive problem in tokamaks. The propagation and the absorption of the lower hybrid electromagnetic wave is a powerful method to generate current drive by wave-particle resonance. However since the wavelength of the electromagnetic field at the LH frequency is very small, a conventional full wave analysis represents an extremely difficult numerical task, far beyond computer capabilities. The problem is addressed in our work, by full wave calculations based on a mathematical finite element technique, allowing for the first time to treat a large volume of plasma. This method is based on a mixed augmented variational (weak) formulation taking account of the divergence constraint and essential boundary conditions, which provides an original and efficient scheme to describe in a global manner both propagation and absorption of electromagnetic waves in plasmas, without the WKB approximation, as done in the conventional ray tracing techniques. This offers the possibility to investigate specific cases for which the usual description fails and may help therefore a general understanding of the LH physics. A first application to a realistic small tokamak configuration is considered in [42] .

Logo Inria