Effect of L-H and H-L Transitions on Tokamak-reactor Operation

The effect of L-H and H-L transitions on the tokamak-reactor operation is considered. Both initial modes are considered as quasi-equilibrium states with the same thermal energy for constant total toroidal currents. A method has been developed for quantification the change in neutron yield in a tokamakreactor during these transitions occurring over times much shorter than the plasma energy confinement time. The method is based on the use of duality of solutions of the GradShafranov equation. The arbitrary functions included in this equation were found as a result of approximation of the normalized plasma pressure profiles, presented versus on the radial flow coordinate obtained at the DIII-D facility. To calculate changes in neutron fluxes during L-H and back H-L transitions, we used these plasma pressure distributions for the ITER device parameters presented in Cartesian coordinates. A numerical calculation showed that in the back H-L transition, a large spike on the global neutron production is possible, which was previously discovered experimentally (ALCATOR-C-Mode, 2001). Since such an increase in neutron fluxes during tokamak-reactor ITER operation poses a serious threat to both the personnel and the facility itself, it is necessary to exclude the possibility of such transitions. Thus, it is necessary to develop such a reactor design that would make it possible to obtain a self-sustaining thermonuclear reaction in the L-mode operation.


Introduction
For many years, tokamaks worked in a mode with not very good plasma confinement, which were later called L-modes [1]. The energy confinement time L E τ in this mode was not enough to fulfill the Lawson criterion [2] and, therefore, to create a tokamak reactor.
In the 80s of the last century, the ASDEX setup [3][4] produced an enhanced confinement mode called the H mode, which gave hope for a self-sustaining fusion reaction.
These two modes are obtained experimentally on all tokamaks. The corresponding plasma energy confinement times L E τ and 2 H L E E τ τ ≈ are well described by ITER-like tokamaks scaling obtained from a large number of installations [5]. The energy confinement time is due to complex turbulent processes depending on the level of fluctuations, profiles, and plasma pressure gradients and radial electric field [6].
The transition from the L mode to the H mode usually occurs when the amount of thermal energy introduced into the plasma exceeds a certain threshold value. We assume that this condition is satisfied. In experiments, the back H -L transition is also observed, which can be spontaneous [6] or caused by some external action [7]. The physical nature and dynamics of transitions are not considered in the framework of this work, since the transitions occurring over times much smaller E τ are analyzed.
It should be noted that although it is assumed that it is possible to operate in the H-mode for a long time due to the continuous input of the required amount of power [8], the operation scenarios of the ITER international tokamak reactor nevertheless imply the possibility of a back HL transition, for example, during discharge quenching in the installation [7-10].

The Thermonuclear Regime in a Tokamak
It is known that the intensity of thermonuclear reactions in a plasma depends on the space distribution of plasma parameters [11,12]. The greater the plasma peaked profiles, the more intense the thermonuclear reaction takes place and the more thermonuclear neutrons are produced. The power released per unit volume of the plasma torus due to the thermonuclear reaction is equal to MeV is the energy released as a result of one DT reaction, D n and T n are the deuteron and tritium ion densities, n is the electron plasma density, > is the D-T fusion reactivity, which is proportional 2 T in the ion temperature range of 10 -20 keV. It follows from (1) that The number of neutrons released in this case from a unit volume of plasma per unit time is The total neutron yield will be From formulas (2) -(4) it follows that for the same amount of thermal energy in both modes of operation of the facility ( L H P P = ), with a different radial distribution of plasma parameters, the total neutron yield will be different.
The coefficient of change in the magnitude of the neutron flux in the back H-L transition is Since in L mode the radial distribution is more picked than in H mode, in the L-H transition 1 LH K ≤ , and in the back H-L transition 1 HL K ≥ .

GS Equation Used for the Experimental Data Approximation
Assuming that the plasma is in a quasi-stationary state in the L and H operating modes of the device, to determine the distribution of its parameters, we will use the GS equation. The solutions of the GS equation will be used to approximate the experimental data.
The GS dimensional equation, which allows one to calculate axially symmetric equilibrium magnetic configurations of a toroidal current-carrying plasma ring with an arbitrary cross-sectional shape of the boundary magnetic surface (TOKAMAK type system), has the form: is a toroidal current density, , , r z ϕ are the cylindrical coordinates. For the convenience of numerical solution of (6), a dimensionless flow is usually introduced values on the magnetic axis and on the plasma boundary. The introduction of such a dimensionless flow at given values max Ψ and b Ψ , and also for a given value of the total current, it allows us to find the only one solution of the GS equation [13].
To find all solutions of the GS equation, it is convenient to introduce a new dimensionless flow in the form having the dimension of the magnetic field of the current, 0 R and a are the major and minor tokamak radii.
In addition, we write are normalization coefficients having dimensions of pressure and magnetic field. Using the introduced values, we write the GS equation in a dimensionless form  (7). With this nondimensionalization, the values of max ψ are not fixed, but are found from the solution of equation (7).
Consider a model that is sometimes used to describe L and H modes in a tokamak. The GS equation for a plasma containing minimal energy while the total toroidal current is constant has the form [14] here ψɶ is the flux function.
To simplify the problem, we represent equation (8) in the form  (9) are found by the following formula To approximate the pressure distribution in the L mode, the L solution of equation (9) should be applied.
Pressure distributions in H modes in the presence of external barriers there are two types. If the profile is concave in the central part of the current-carrying plasma ring, then this is the mode of the first type H1, if it is convex, then this is the mode of the second type H2.
To approximate the pressure distribution in the H mode of the first type, we will use the L solution, and for the pressure distribution in H2 we will use the H solution of the following equation is function approximating the Heaviside function ( ) (11) describes the external barrier.
Here Ω determines the height of the barrier, c ψ -its location, and b characterizes the decrease in the flux near the plasma boundary.
The plasma pressure distributions determined by H and L solutions of equation (11) have the form where for H1 mode 0 c = . In this paper, normalization coefficients are determined from the condition of equality of thermal energy for all the three considered pressure profiles

Evaluation of Neutron Flow Change on Back H-L Transition
The neutron yield coefficient for the back H-L transition in this work will be estimated by formula (5), using the distribution of squared pressure profiles before and after the transition.
In [15], a thorough comparison of experimental data for plasma with various parameters measured using the D-IIID device with the results of calculations using various codes was performed. These results allow us to simulate the plasma parameters expected at the ITER tokamak ( p β =0.65, A =3.1, elongation is 1.7, triangularity is δ =0.33). The eigenvalues of equation (11) for H-mode depend not only on the tokamak parameters, but also on the quantities characterizing the barrier.
The least squares method was used to approximate the experimental dependences of the plasma pressures. With this approximation, variation was carried out for all four parameters ( , , c λ ψ Ω and b ). Figure 2 shows the dependence of the normalized plasma pressure distribution versus coordinate * ρ for L mode. Points are the data from [15]. The approximation is performed according to the formula (10) with 0.4 L λ = . Figure 2. The normalized pressure in the L mode as a function of ρ*. Points are the data from [15], the curve is result of approximation. Figure 3 shows the dependence for the normalized pressure distribution in the H1 mode carried out using the L solution found by formula (11) Figure 4 shows the dependence of the normalized pressure for the H2 mode. The data were taken from [10].
For this approximation, we used the H solution of the equation (11) where Hev is the Heaviside function. Figure 3. The normalized pressure in H1 mode as a function of ρ*. Points are the data from [15], the curve is result of approximation.
To calculate the gain of the neutron flux, it is necessary to have a distribution of the square of the plasma pressure, described by a two-dimensional function in real coordinates (x, y). Figures 5-7 shows these two-dimensional plasma pressure distributions plotted in the coordinates (x, y) ( Figure  5 -L mode, Figure 6 -H1 mode, Figure 7 -H2 mode).
It can be seen from the figures that the space distribution of plasma pressure in tokamak has a complex shape. The maximum value of the plasma pressure in the center of the plasma ring in L mode exceeds this value in H1 mode by about 6 times, and in H2 mode by 20 times.
In this case, the thermal energy in all three modes is the same and the full toroidal current is preserved.
The increase in neutron yield in the reverse H-L transition was estimated using formula (5), which includes the squared pressure profiles in real space before and after the transition.
The neutron flux increase coefficient was calculated for distributions simulating transitions in ITER operating in the paramagnetic region.   If the transport stiffnes or advance scenario is implemented, then H1-L or H2-L occurs and the neutron flux increase coefficient is 3 HL K = . In 2001, the so-called Enhanced Neutron Mode ITB was discovered on the ALCATOR C-Mod tokamak [16]. In this mode, during the reverse H-L transition, a pulsed increase in the neutron flux occurs. The duration of this pulse is several cycles of sawtooth oscillations or approximately equal to one plasma energy confinement time. With such a dynamic transition, the neutron flux increases in the range from 1.3 to 8 times.
Since the duration of a neutron pulse is of the order of time of E τ formula (5) for estimating the coefficient HL K , strictly speaking, in this case it is not applicable, nevertheless, there is a qualitative agreement between the experimental data and the results of theoretical calculations. Since the tokamak design and biological protection at the ITER device are not designed for such thermal and neutron loads, even a single back H-L transition can lead to the destruction of the installation itself and pose a serious danger to the environment.

Conclusions.
Thus, it is necessary to exclude the possibility of back H-L transitions. It should be noted that in the literature it is discussed the thermonuclear reactor designes which can operate in L mode [17,18].

Appendix
FlexPDE program. In this work, to solve the Grad-Shafranov equations, the FlexPDE program was used. This program is available at www.pdesolution.com/v6student.html and other similar sites.
This program is designed to solve nonlinear partial differential equations using the finite element method. Such equations may have one, several, or no solutions.
The program allows us to find all solutions of such equations.
As a test of the program, we consider the next equation [14]   It should be noted that the value ZNI for GS equations (9) and (11) it is necessary to select for each set of parameters included in these equations ( , A, , , , c b c λ η ψ , Ω ). In the problems solved in this paper, for determining distributions in the H mode, always 1 ZNI < . For L mode, depending on the parameters of the problem ZNI varies from 1.2 to about 3.