II. Theory, Simulation and Modeling

For the NEXT (Numerical Experiment of Tokamak) project in 2007, a new conservative gyrokinetic full-f Vlasov code has been developed in order to realize a long time simulation. Zonal field generation was investigated in finite beta ion temperature gradient driven turbulence by global Landau-fluid simulations. Also, nonlinear MHD simulations found the Alfven resonance effects on the evolution of magnetic islands driven by an externally applied perturbation. The distribution of eigenvalues of the resistive MHD equations in the complex plane has been re-investigated for smaller resistivity than the previous works, and a numerical matching scheme for linear MHD stability analysis was proposed in a form offering tractable numerical implementation. A conjugate variable method was formulated in order to apply the Hamilton-Lie perturbation theory to a system of ordinary differential equations that does not have the Hamilton dynamic structure.

The development of the integrated modeling is promoting in order to study complex behaviors of JT-60U plasmas and to predict the performance of future burning plasmas. Several kinds of element codes have been developed for the integrated modeling; MHD stability code, integrated SOL/divertor code, fundamental SOL/divertor code of the particle model, heating and current drive analysis code, and integrated transport code TOPICS-IB. Achievements in 2007 with the use of these codes are described below in section 2.

1. Numerical Experiment of Tokamak (NEXT)
1.1 Magnetohydrodynamic (MHD) Theory and Simulation
1.1.1 Alfven Resonance Effects on an Externally Driven Magnetic Island in Rotating Plasma

The Alfven resonance effects on the magnetic island evolution driven by an externally applied perturbation were investigated. In a low viscosity regime, perturbed current sheets form at the Alfven resonance surfaces (Fig.II.1.1-1), which differ from the radial position of the magnetic neutral surface. Therefore, the sheets exist outside the inner non-ideal layer, defined for the non-rotating plasma. According to this perturbed current sheet profile, the total torque, which affects the plasma, extends wider than the radial position of the Alfven resonance. It is found that the radial position of an Alfven resonance moves to the magnetic neutral surface and causes the rapid growth of a driven magnetic island (Fig.II.1.1-2). This is because the background flow is damped by the appearance of the magnetic island, which operates torque on the plasma. These features are inconsistent with former theoretical assumptions, which enable using the asymptotic matching method to estimate the force balance and the critical value of the external perturbation, beyond which the driven magnetic island grows rapidly [1.1-1].

Fig.II.1.1-1      Fig.II.1.1-2

1.1.2 MHD Spectrum of Resistive Modes

The resistive MHD spectrum is investigated in detail by solving the eigenvalue problem of the reduced MHD equations in cylindrical tokamak plasmas, in particular for asymptotically smaller resistivity than the previous works. In the presence of the resistivity, the eigenvalues of linear resistive MHD modes are not necessarily purely real, but have an imaginary part. The eigenvalues are classified as the continuum spectrum of dumping modes on the negative real axis, the discrete spectrum of oscillatory and dumping modes on the complex plane, and the discrete spectrum of unstable modes on the positive real axis. For a wide range of the resistivity parameter (η=10-4 ~ 1.5×10-6), the shape and location of eigenvalue distribution for m/n=1/1 modes which have the same resonant surface at q=1 is almost independent of the resistivity, although the density of the eigenvalues increases. However, it is found that for a further small resistivity regime, the eigenvalue distribution changes sensitively depending on the resistivity [1.1-2].

1.1.3 Numerical Matching Scheme for Linear MHD Stability Analysis

A new matching scheme for linear MHD stability analysis is proposed in a form offering tractable numerical implementation. This scheme divides the plasma region into outer regions and inner layers, as in the conventional matching method. However, the outer regions do not contain any rational surface at their terminal points; an inner layer contains a rational surface as an interior point. The Newcomb equation is therefore regular in the outer regions. The MHD equation employed in the layers is solved as an evolution equation in time, and the full implicit scheme is used to yield an inhomogeneous differential equation for space coordinates. The matching conditions are derived from the condition that the radial component of the solution in the layer is smoothly connected to those in the outer regions at the terminal points. The proposed scheme was applied to the linear ideal MHD equation in a cylindrical configuration, and was proved to be effective from the viewpoint of a numerical scheme [1.1-3].

1.1.4 Conjugate Variable Method in the Hamilton-Lie Perturbation Theory -Applications to Plasma Physics-

The conjugate variable method, which is an essential ingredient in the path-integral formalism of classical statistical dynamics, was used in order to apply the Hamilton-Lie perturbation theory to a system of ordinary differential equations that does not have the Hamilton dynamic structure. The method endows the system with the Hamilton dynamic structure by doubling the unknown variables; hence the canonical Hamilton-Lie perturbation theory becomes applicable to the system. The method was applied to two classical problems known in plasma physics to demonstrate the effectiveness and to study the property of the method: one is a non-linear oscillator that can explode; the other is the guiding center motion of a charged particle in a magnetic field [1.1-4].

1.1-1 Ishii, Y., et al., "Formation and long term evolution of the externally driven magnetic island in rotating plasmas," to be published in Plasma and Fusion Research.
1.1-2 Matsumoto, T., et al., Proc. 6th General Scientific Assembly of the Asia Plasma and Fusion Association (India, 2007).
1.1-3 Kagei, Y., et al., "Numerical matching scheme for linear magnetohydrodynamic stability analysis," to be published in Plasma and Fusion Research.
1.1-4 Tokuda, S., et al., "Conjugate variable method in the Hamilton-Lie perturbation theory -applications to plasma physics-," submitted to Plasma and Fusion Research.

1.2 Plasma Turbulence Simulation
1.2.1 Development of Non-Dissipative Conservative Finite Difference Scheme for Gyrokinetic Vlasov Simulation

A new conservative gyrokinetic full-f Vlasov code is developed using a finite difference operator which conserves both the L1 and L2 norms. The growth of numerical oscillations is suppressed by conserving the L2 norm, and the code is numerically stable and robust in a long time simulation. In the slab ITG turbulence simulation, the energy conservation and the entropy balance relation are confirmed, and solutions are benchmarked against a conventional δf particle-in-cell (PIC) code. The results show that the exact particle number conservation and the good energy conservation in the conservative Vlasov simulation are advantageous for a long time micro-turbulence simulation (Fig. II.1.2-1). In the comparison, physical and numerical effects of the v// nonlinearity are clarified for the Vlasov and PIC simulations [1.2-1, 2].

1.2.2 Zonal Field Generation and Its Effects on Zonal Flow and Turbulent Transport

Global Landau-fluid simulations of ion temperature gradient (ITG) driven turbulence have been performed for finite beta tokamak plasmas, where the beta value is a ratio of plasma pressure to magnetic pressure. The ITG turbulence can drive zonal magnetic fields as well as zonal flows in finite beta cases. The Reynolds stress drives zonal flows and the geodesic transfer effect acts as a sink for zonal flows usually. It is found that the Reynolds stress and the geodesic transfer effect change their roles at low order rational surfaces where zonal magnetic fields are generated the most strongly. This is not observed in electrostatic (zero beta) simulations. Effects of the zonal magnetic fields on the zonal flows and the turbulent transport are limited in a small region around a low order rational surface at least in a low beta regime where the ITG mode is dominant. The zonal magnetic fields, however, may affect the zonal flows and the turbulent transport in a high beta regime because amplitude of the zonal magnetic fields increases with the beta value [1.2-3].

1.2-1 Idomura, Y., et al., J. Comput. Phys. 226, 244 (2007).
1.2-2 Idomura, Y., et al., Commun. Nonlinear Sci. Numer. Simul. 13, 227 (2007).
1.2-3 Miyato, N., et al., Proc. 34th EPS Plasma Phys. Conf. (Poland, 2007), p4.043.

2. Integrated Modeling
2.1 MHD Stability - Effect of Equilibrium Properties on the Structure of the Edge MHD Modes in Tokamaks -

The ideal MHD stability code MARG2D has been extended to estimate a growth rate of the MHD mode under the incompressible assumption by introducing the plasma inertia [2.1-1]. With this extension, MARG2D realizes not only to identify the stability boundary of ideal MHD modes, but also to investigate physical properties of unstable MHD modes in detail.

By using this extended MARG2D, effects of the pressure profile and the current density profile inside the top of pedestal and that of the plasma shape on the expansion of the structure of the unstable edge MHD mode are investigated numerically [2.1-2]. The radial structure of the edge MHD mode is expanded by spreading the envelope of the edge ballooning mode due to increasing the pressure gradient inside the top of pedestal. Moreover, the increase of the current density induces the decrease in the toroidal mode number of the most unstable mode, and this decrease also expands the structure of the unstable mode.

The mode structure is subject to expanding in strongly shaped plasmas. This is because the pressure gradient inside the top of pedestal can approach to the ballooning mode stability boundary and the current density increases enough to reduce the toroidal mode number of the most unstable mode. These increases of the pressure gradient and the current density destabilize the edge MHD mode and expand the mode structure.

2.1-1 Aiba, N., et al., J. Plasma Fusion Res. 2, 010 (2007).
2.1-2 Aiba, N., et al., to be published in J. Phys. Conf. Series (2008).

2.2 SOL-Divertor

Divertor of tokamak reactors has four major functions, heat removal, helium ash exhaust, impurity retention, and density control. Such divertor performance strongly depends on the various physics, i.e. plasma transport, kinetic effects, atomic processes, and plasma-wall interactions. In order to understand complicated divertor physics and to predict divertor performance, JAEA has developed a series of divertor codes, onion-skin modeling, SOLDOR, NEUT2D, IMPMC, PARASOL, 5-point divertor code coupled with a core transport code (TOPICS-IB), and Core-SOL-Divertor model (CSD).

The benchmark test of SOLDOR/NEUT2D code with B2/EIRENE code was attempted. The simulation study of JT-60SA divertor was carried out with B2/EIRENE and the difference between single-null and double-null configurations was confirmed [2.2-1].

2.2.1 Development of Integrated SOL/Divertor Code and Simulation Study

To predict the heat and particle controllability in the divertor of tokamak reactors and to optimize the divertor design, comprehensive simulations by integrated modeling allowing for various physical processes are indispensable. SOL/divertor codes have been developed in Japan Atomic Energy Agency for the interpretation and the prediction on behavior of SOL/divertor plasmas, neutrals and impurities [2.2-3]. The code system consists of the 2D fluid code SOLDOR, the neutral Monte-Carlo (MC) code NEUT2D, and the impurity MC code IMPMC. Their integration code "SONIC" is almost completed and examined to simulate self-consistently the SOL/divertor plasmas in JT-60U. In order to establish the physics modeling used in fluid simulations, the particle simulation code PARASOL has also been developed.

Simulation studies using those codes are progressed with the analysis of JT-60U experiments and the divertor designing of JT-60SA. The X-point MARFE in the JT-60U experiment is simulated. It is found that the deep penetration of chemically sputtered carbon at the dome for the detached phase causes the large radiation peaking near the X-point as shown in Fig. II.2.2-1. The pumping capability of JT-60SA is evaluated through the simulation. A guideline to enhance the pumping efficiency is obtained in terms of the exhaust slot width and the strike point distance. Transient behavior of SOL/divertor plasmas after an ELM crash is characterized by the PARASOL simulation; the fast-time-scale heat transport is affected by collisions while the slow-time-scale behavior is affected by the recycling.

2.2.2 Extension of IMPMC Code toward Time Evolution Simulation

As a self-consistent modelling of divertor plasma and impurity transport, the SONIC code package has been developed. The key feature of this integrated code is to incorporate the impurity MC code, IMPMC. The MC approach is suitable for modelling of interactions between impurities and walls, including kinetic effects, and the complicated dissociation process of hydro-carbons. The MC modelling, however, has the disadvantage for long computational time, large MC noise, and assumption of steady state. The first and second difficulties were solved by developing a new diffusion model and optimizing with a Message Passing Interface (MPI) on the massive parallel computer. The third subject is solved by extension of IMPMC code toward time evolution simulation. In time-dependent simulation with the MC code, a serious problem to increase number of test particles. The particle reduction method consisting of sorting the weights, pairing and Russian roulette has been proposed [2.2-4]. Sorting of the weights is indispensable to suppress the MC noise.

The divertor configuration in JT-60SA has been optimized from a viewpoint of the neutral recycling with SOLDOR/NEUT2D. In the near future, it will be further optimized from a viewpoint of the impurity control with the SONIC code package coupled with the above extended IMPMC.

2.2-1 Suzuki, Y., et al., Contrib. Plasma Phys. 48, 169 (2008).
2.2-2 Hiwatari, R., et al., Contrib. Plasma Phys. 48, 174 (2008).
2.2-3 Kawashima, H., et al., Plasma Phys. Control. Fusion 49, S77 (2007).
2.2-4 Shimizu, K., et al., Contrib. Plasma Phys. 48, 270 (2008).

2.3 ELM Transport - Effect of Radial Transport Loss on the Asymmetry of ELM Heat Flux -

Large heat load on the divertor plate intermittently produced by ELMs is one of crucial issues for the tokamak fusion reactor research and development. The imbalance in the ELM heat loads on in-out divertor plates is also the problem. It has been reported that the ELM heat deposition to the outer plate is larger than that to the inner plate in JT-60U, while in JET and ASDEX Upgrade the inner-plate heat deposition becomes larger than the outer-plate heat deposition. To understand the physics background of such complex behaviors of the ELM heat flux, the dynamics of SOL-divertor plasmas after an ELM crash is studied with a one-dimensional particle simulation code, PARASOL [2.3-1].

The ELM crash occurs off-centrally in the SOL region, and the ELM heat flux to the near divertor plate and that to the far divertor plate become asymmetric. The peak heat flux to the near plate is larger as compared to the far plate. The asymmetry in the peak heat flux increases with the connection-length ratio. The imbalance in the heat deposition, however, is small. The radial transport loss of ELM flux creates the asymmetry in the heat deposition, but the imbalance is still not large even for the large radial transport loss rate. The electron heat flux to the far plate brought by the SOL current is one of the causes of a small imbalance in the heat deposition. Another cause is the asymmetric SOL flow and its convective heat flux, whose stagnation point stays for a long period near the ELM crash location.

2.3-1 Takizuka T., et al., Contrib. Plasma Phys. 48, 207 (2008).

2.4 Heating and Current Drive - Electron Cyclotron Current Drive in Magnetic Islands of Neo-classical Tearing Mode -

Electron cyclotron current drive (ECCD) is the effective method for stabilization of neo-classical tearing modes (NTM). The driven current is evaluated conventionally by the bounce averaged Fokker-Planck equation (BAFP), where the magnetic field is assumed to be axi-symmetric. When the magnetic islands are formed by NTMs, however, this assumption is incomplete and the validity of ECCD analysis based on the BAFP equation becomes questionable.

The ECCD in the magnetic island is studied numerically by a particle simulation [2.4-1]. Drift motion of electrons with Coulomb collisions and velocity diffusion due to the EC waves is tracked by using a Monte-Carlo technique. An EC resonance region is located around the O-point and localized along the toroidal direction. It is found that the driven current is well confined in the helical flux tube including the EC resonance region and the current channel looks like a helical "Snake". Figure II.2.4-1 shows the dependence of EC current drive efficiency I/P on the EC power P. Plasma parameters are the followings; the magnetic island of m/n = 2/1 is located around 0.5m in the minor-radius direction, where ne = 3×1019m-3, Te = 10keV and Zeff = 1. The major radius is 3.5m and the EC wave frequency is 125GHz (2nd harmonic resonance). The current drive efficiency with the magnetic island is larger than that without magnetic island and is enhanced by the increase of the EC power. The driven current profile in the magnetic island becomes steep around the O-point with the increase of P. These results are caused by the good confinement of EC resonant electrons inside the island like "Snake" and the nonlinear kinetic effect due to the high EC power density.

2.4-1 Hamamatsu K., et al., Plasma Phys. Control. Fusion 49, 1955 (2007).

2.5 Integrated Simulation

Integrated simulation models have been developed on the basis of the research in JT-60U experiments and first-principle simulations in order to clarify complex features of reactor-relevant plasmas. The integrated model of edge-pedestal, SOL and divertor clarified that the steep pressure gradient inside the H-mode pedestal top enhances the ELM energy loss [2.5-1]. A new one-dimensional core transport code, which can describe the radial electric field and plasma rotations, has been developed [2.5-2]. Success in these analysis and development leads to the further effective study of complex plasmas and methods to control the integrated performance.

2.5.1 Integrated ELM Simulation with Edge MHD Stability and Transport of SOL-Divertor Plasmas

The energy loss due to ELMs has been investigated by using an integrated simulation code TOPICS-IB based on a 1.5 dimensional core transport code with a stability code for the peeling-ballooning modes and a transport model for SOL and divertor plasmas. In the previous study, the TOPICS-IB successfully simulated transient behaviors of an H-mode plasma and clarified the mechanism of the collisionality dependence of the ELM energy loss. The ELM energy loss, however, was less than 10% of the pedestal energy.

The effect of the pressure gradient inside the pedestal top, P'inped, on the ELM energy loss is examined [2.5-1]. Figure II.2.5-1(a) shows profiles of the total pressure, P, at the ELM onset. The transport is reduced to the ion neoclassical level in the pedestal region for case A, and the transport is additionally reduced inside the pedestal top for other cases B-D. The pedestal top is located at ρ = 0.925 for all cases and P'inped becomes steeper in order of A, B, C, D. Even for the case A, P'inped is as the same as that observed in JT-60U. Profiles of the ELM enhanced diffusivity, χELM, in the cases A-D are shown in Fig. II.2.5-1(b). The steep pressure gradient broadens eigenfunctions of unstable modes and the region of the ELM enhanced transport. Figure II.2.5-1(c) shows the ELM energy loss, ΔWELM, normalized by the pedestal energy, Wped, as a function of P'inped/P'ped where P'ped is the pedestal pressure gradient. In the case A, the ELM energy loss is less than 10% of the pedestal energy and is comparable with those in JT-60U. The steep pressure gradient inside the pedestal top enhances the ELM energy loss. The density collapse, which is not considered here, enhances the values of ΔWELM/Wped by about 50% under the assumption of the similar collapse to the temperature one. The ELM energy loss in the simulation becomes larger than 15% of the pedestal energy, as is shown in the database of multi-machine experiments.

2.5.2 Dynamic Transport Simulation Code Including Plasma Rotation and Radial Electric Field

A new one-dimensional multi-fluid transport code named TASK/TX has been developed [2.5-2]. This code is able to describe dynamic behavior of tokamak plasmas, especially the time-evolution of the radial electric field and the plasma rotations. A set of flux-surface averaged equations is solved self-consistently in the cylindrical coordinates; Maxwell's equations, continuity equations, equations of motion, heat transport equations, momentum transfer equations for fast particles, and two-group neutral diffusion equations.

The finite element method with a piecewise linear interpolation function is employed. The Streamline Upwind Petrov-Galerkin method is also incorporated for numerically robust calculation. Despite solving the very nonlinear equations, the code shows a good convergence performance.

Modification of a density profile during neutral beam injection (NBI) is presented. We found the density peaking for the counter-NBI and the density broadening for the co-NBI. The balance between neoclassical and turbulent effects defines the status of the density profile. We have confirmed that the TASK/TX well reproduces the profiles observed in JFT-2M.

In the presence of ion orbit losses, the code predicts the generation of the intrinsic (spontaneous) rotation in the counter direction with the inward radial electric field. The non-ambipolar loss breaks quasi-neutrality and the plasma instantaneously generates the inward radial current near the periphery to compensate the ion loss current, inducing the more negative radial electric field and the torque toward the counter direction. Other conventional transport codes assuming the quasi-neutrality cannot follow these processes. It is the very special characteristic of the TASK/TX code that there is no need to impose an explicit quasi-neutrality condition.

2.5-1 Hayashi, N., et al., Contrib. Plasma Phys. 48, 196 (2008).
2.5-2 Honda, M., et al., J. Comput.Phys. 227, 2808 (2008).

3. Atomic and Molecular Data

We have been producing, collecting and compiling cross-section data for atomic and molecular collisions and spectral data relevant to fusion research.

The electron capture and the electron loss cross-sections of singly ionized tungsten by collision with CH4 and C2H6 have been measured at collision energies of 27.2 and 54.4 eV/amu. The state selective charge transfer cross-section data of B5+ and C6+ by collision with H* ( n = 2 ) in the collision energy range between 62 eV/amu and 6.2 keV/amu have been calculated with a molecular-bases close-coupling method. The cross-section data for 42 processes of collisions of He, He*, He-, He+, He2+ and 3He2+ with H, H-, H2, He and He+ have been complied. The recommended cross-section data are expressed with analytic functions to facilitate practical use of the data. The compiled data are in preparation for the Web at the URL http://www-jt60.naka.jaea.go.jp/english/JEAMDL/. The charge transfer data published in 2007 have been collected, and the database for the chemical sputtering yield data of graphite materials with hydrogen isotope collisions have been established.

Page Top