Analysis of the Thermonuclear Instability including Low-Power ICRH Minority Heating in IGNITOR
AAnalysis of the Thermonuclear Instability including Low-PowerICRH Minority Heating in IGNITOR Alessandro CARDINALI ENEA for EUROfusion, Via E. Fermi 45, 00044 Frascati - ItalyEmail: [email protected] SONNINO , Department of Theoretical Physics and MathematicsUniversit´e Libre de Bruxelles (ULB), Campus Plain CP 231Boulevard de Triomphe, 1050 Brussels, Belgium.& Royal Military School (RMS)Av. de la Renaissance 30 1000 Brussels - Belgium.Email: [email protected]
Abstract
The nonlinear thermal balance equation for classical plasma in a toroidal geometry is analyti-cally and numerically investigated including ICRH power. The determination of the equilibriumtemperature and the analysis of the stability of the solution are performed by solving the energybalance equation that includes the transport relations obtained by the classical kinetic theory. Anestimation of the confinement time is also provided. We show that the ICRH heating in the IGN-ITOR experiment, among other applications, is expected to be used to trigger the thermonuclearinstability. Here a scenario is considered where IGNITOR is led to operate in a slightly sub-criticalregime by adding a small fraction of He to the nominal 50%-50% Deuterium-Tritium mixture.The difference between power lost and alpha heating is compensated by additional ICRH heating,which should be able to increase the global plasma temperature via collisions between He minorityand the background D − T ions.PACS numbers: 28.52.-s, 28.52.Av a r X i v : . [ phy s i c s . p l a s m - ph ] J un . INTRODUCTION Tokamak with a strong magnetic field like IGNITOR operates on the low-temperaturebranch of the ignition boundary [1],[2] and [3] making impossible a stationary fusion reactiondue to the thermonuclear instability. As a consequence of the instability, the self-heatingof the plasma by alpha particles induces a significant rise of its temperature accompaniedby an increase in the pressure, which in its turn will reinforce the thermal instability ofthe plasma. There has been a great effort, in the last decades, in investigating the variousmechanism proposed for controlling the fusion thermal instability [4]. In some work it wasproposed that the balance in the growth of thermonuclear power be stabilized by increasingthe energy losses from the plasma by changing the major radius R . Increasing R ( i.e. ,adiabatic expansion) there will be a reduction of the plasma temperature. However thereare serious engineering difficulties with this approach. Thus the possibility of significantlychanging the large radius increases the volume of the chamber, which, obviously, will increasethe volume of the magnet system [5]. In other references it was suggested that α -power couldbe regulated, by injecting pellets of fuel [6]. This method has significant advantages dueto the technological progress of these last years in injecting a fuel tablet up to the centerof the plasma column; this is connected with the fact that the tablet reaches a relativelyhigh velocity [7]. In addition, the required injection rate ( ∼ Hz ) is technically easyto achieve. The only difficulty remains that after the tablet is injected, the decrease inthe cross section of the D − T fusion reaction ( < σv > ∼ T ) will be compensated bythe increase in the density, and the intensity of the thermonuclear reactions will remainunchanged. Control with modulation of the fueling rate and high- Z impurity injection hasalso been demonstrated as an effective means for controlling the fusion thermal instability[8], especially when auxiliary power modulation cannot be used. The effects of a numberof other phenomena on controlling the fusion thermal instability have been examined, andare: i) transport losses due to toroidal magnetic field ripple via the τ E term; ii) impurityinjection; iii) the poloidal divertor; iv) a soft beta limit; v) compressing or decompressingthe plasma; vi) an ergodic magnetic limiter; vii) modulation of divertor pumping; viii) modification of alpha-particle transport; ix) saw tooth oscillations; and x) radial motion.A very exhaustive reference can be found in Ref. [4]. In other works it was proposed thatthe power of the thermonuclear burning be stabilized at a fixed level by regulation of the2ower of additional heating [9], [10] and [11]. In particular in Ref. [11] it was proposedthat the reactor is operating in the sub-critical regime, i.e. , the parameters of the plasmaare chosen so that the power of the thermonuclear reactions is slightly less than the powerlost, for example by adding to the Deuterium-Tritium mixture a small fraction of He (fewpercent); this small fraction of impurity unbalances the ideal ignition condition (50% − D − T ), and the difference is compensated by additional heating. ICRH is, in fact, able toheat directly the minority species (ICRH minority heating) and by collision to transfer thepower to the main species of the plasma: electrons and deuterium-tritium ions, by increasingthe plasma temperature. The ICRH power acts to regulate the thermonuclear power vianegative feedback. In this work this approach is accurately studied by solving the energybalance equation including the additional ICRH heating. Here the problems of ensuringstability of burning and the quality of transient processes for different confinement laws arestudied.The main purpose of our work is to estimate the equilibrium temperature and the energyconfinement time by assuming that the transport is governed by the classical kinetic theory.However, it is well known that there is almost always a strong anomalous diffusion in theouter plasma, for which a number of heuristic experimental scalings exist. So the classicalthermal diffusion should be understood as an ideal shape , since the coefficient varies stronglywith inverse temperature. Hence, we certainly do not claim to have estimated the confine-ment time in real experimental conditions, but only to show the results obtained by usingthe classical kinetic theory, applied to plasmas in toroidal geometry, without the aid of an ad hoc transport model. The next step, will be then to consider all the three collisionaltransport regimes (ı.e., the classical, Pfirsch-Schl¨uter and banana transport regimes) and tocompare the analytical results with the solutions obtained by adopting a turbulent transportmodel like the gyro-Bohm model. This will be subject of future works.The manuscript is organized as follows. In section I the thermal energy balance equationin the framework of the neoclassical theory (see, for example, Ref. [12]) is recalled andall the terms of power gain and lost are specified. The energy-gain and -loss terms areevaluated in Sec. II. Sec. III is devoted to the determination and the discussion of theequilibrium solutions. The stationary thermal profile is determined by making use of theplasma dynamical equations and of the transport relations, which are rigorously obtainedby kinetic theory. The estimation of the confinement time and the stability of the steady3tate thermal solution can be found in Sec. IV and in Sec. V, respectively (a simplifiedcalculation of the unstable modes is reported in the Appendix). The role of the ICRH powermodulation in stabilizing/destabilizing the phenomenon is also herein discussed. Conclusionsare reported in Sec. VI. II. ONE-DIMENSIONAL TOROIDAL PLASMADYNAMICAL EQUATIONS INTHE STANDARD MODEL AND EVALUATION OF THE ENERGY-GAIN AND-LOSS TERMS
Our first objective is to determine the electron temperature profiles of species ι ( ι = e forelectrons and ι = i for ions). This task will be accomplished by considering the balanceequations (mass and energy) for species ι in toroidal geometry and by adopting the validityof the following standard model for the magnetic configuration B = B (cid:16) (cid:15)ρq ( ρ ) e θ + 11 + (cid:15)ρ e φ (cid:17) (1)Here, e φ and e θ are the versors in the toroidal and poloidal directions, respectively, (cid:15) is theinverse of the aspect ratio: (cid:15) = a/R (with a = 47 cm and R = 132 cm denoting the minorand the major radius of IGNITOR, respectively), and B is the toroidal magnetic field atthe magnetic axis (which in IGNITOR increases from B = 8 T to 13 T , during the ramp-upphase). ρ and q ( ρ ) denote the normalized minor radius ( ρ ≡ r/a ) and the safety factor,respectively. The safety factor profile we shall use for IGNITOR-plasma in this work iscompatible with a plasma current of 11 M A at the end of the ramp-up phase and is greaterthen 1 on the plasma magnetic axis. The equations of one-dimensional plasma dynamics, intoroidal geometry, by assuming the validity of the standard model, can be brought into theform (see, for example, Ref. [12]) ∂n e ∂t = − r ∂∂r (cid:16) r < γ er > (cid:17) ∂p∂t + 1 r ∂∂r (cid:104) r (cid:16) < q e > + < q i > + 52 (1 + Z − ) T e < γ er > (cid:17)(cid:105) = c π E B φ Rr ∂∂r (cid:16) r q ( r ) (cid:17) + S gain − loss (2)where p e and p i are the plasma pressure due to the electrons and ions, respectively, and n e , T e and Z are electron density, electron temperature and the ion charge number, respectively.4ere, < · · · > denotes the surface-average operation. < q ι > and < γ er > are the averagedradial heat flux of species ι and the averaged electron flux, respectively. E is the externalelectric field at ρ = 0, and S gain − loss is the source term, i.e. the loss and energy-gain.Equation (2) must be completed with the transport equations, i.e. with the thermodynamicflux force relations, in order to close the plasma dynamical equations. We make now severalassumptions and approximations for reducing Eqs (2) to a much simpler form. First, weassume that, in Eqs. (2), the contributions related to the averaged electron flux < γ er > andthe external electric field E , may be neglected with respect to the other terms. Second,the fuel is assumed to consist of a 50%-50% mixture of Deuterium ( D ) and Tritium ( T ),with a negligible concentration of α -particles and He (2-3%). Third, the temperature ofthe plasma is the same for all species: T e = T D = T T = T . Then Eqs (2), reduce to ∂∂t (cid:16) p e (cid:88) i = D,T p i (cid:17) + 1 r ∂∂r (cid:104) r (cid:0) < q e > + < q i > (cid:1)(cid:105) = S gain − loss (3)From the local electro-neutrality condition we get n e = Z T n T + Z D n D = n (4)and we have taken into account that Z D = Z T = 1. In the calculation, we chose the followingprofile for the electron density: n e ( ρ ) = n ae (1 − ρ ) + n be , with n ae = 8 . × ( cm − ) and n be = 5 × ( cm − ), respectively. The total hydrodynamic pressure term is provided by thestate equation p = p e + p D + p T = 2 nT (5)and Eq. (3) may be rewritten as32 ∂p∂t + 1 r ∂∂r (cid:104) r (cid:0) < q e > + < q D > + < q T > (cid:1)(cid:105) = S gain − loss (6)Before discussing on the structure of the heat flux term (to which we will dedicate Sec. III),we determine the structure of the loss-gain terms on the r.h.s. of Eq. (6). Note that in ourequations we are assuming the physical quantities expressed in cgs units (unless differentlyspecified) so that the pressure is given in terms of [ m ][ l ] − [ t ] − , and the heat flux [ M ][ t ] − ,in this manner Eq. (6) has the dimension of a power density, and it represents the powerdensity balance. The term on the r.h.s. of Eq. (6) (the power gain-loss term) is specified asfollows: S gain − loss = Q α + Q b + Q add (7)5here Q α is the alpha heating power, Q b is the radiation loss (Bremsstrahlung), and Q add isthe additional ICRH heating respectively. The alpha heating power density is given by thefollowing formula Q α = n < σv > D − T E α (8)where E α is the energy at which the alpha particles are created (3 . M eV ), σ is the reactioncross section and it is a measure of the probability of a fusion reaction as a function ofthe relative velocity of the two reactant nuclei, given in barn [ l ] − . If the reactants havea distribution of velocities, e.g. a thermal distribution with thermonuclear fusion, then itis useful to perform an average over the distributions of the product of cross section andvelocity i.e. < σv > D − T in units [ l ] [ t ] − . The reaction rate (fusions per volume per time)is < σv > times the product of the reactant number densities, (1 + δ ij ) − n i n j < σv > D − T , ifa species of nuclei (deuterium) is reacting with another species (tritium), such as the D − T reaction at 50%, then the product n i n j must be replaced by n < σv > D − T , which increasesfrom virtually zero at room temperatures up to meaningful magnitudes at temperatures of10 − keV . At these temperatures, well above typical ionization energies (13 . eV inthe hydrogen case), the fusion reactants exist in a plasma state. In our calculation we haveassumed that the dependence of the cross-section σ on temperature is given by an analytical3-parameter fitting [13], σ ( E ) = π (2 µ/ (cid:126) ) E ( m b /m a + m b ) 1 θ × ( − C )( C + C /E ) + ( C − /θ ) (9)where E is the energy of incident particles in the laboratory system, θ = 12 π (cid:104) exp (cid:104) π (cid:112) µE/ (cid:126) [ (cid:126) / ( µZ i Z j e )] (cid:105) − (cid:105) is the Gamow penetration factor, µ is the reduced mass, (cid:126) is the Planck constant, and thefitting points for the 3-parameter fit formula in Deuterium Tritium reaction are: C = 0 . C = 0 . C = 0 . < σv > in terms of theplasma temperature. The calculation of < σv > D − T can be performed by averaging thecross section over the relative velocities of the reactants keeping the relative ion energydistribution in plasma to be Maxwellian < σv > D − T = 1 m (cid:90) ∞ σ ( E ) exp (cid:16) − EK B T (cid:17) dE where m is the mass particle, and K B the Boltzmann constant.6he term Q b is the radiation loss (Bremsstrahlung) Q b = − . × − n ( cm − ) Z eff T / ( eV ) (10)This formula gives directly the power density loss for bremsstrahlung in cgs units with thetemperature is measured in eV and the number density in cm − . Finally the term due to theICRH power absorption corresponds to Q add = Q ICRH ( r ). To this end it is useful to recallthat the power deposition profile in IGNITOR ignited scenario with ramping magnetic field[from 8 to 13 T ] has been analyzed in detail in Ref. [14] by using the ICRH full wave codeTORIC [15], coupled to the quasi-linear Fokker-Planck routine SSFPQL [16]-[17]. The bestICRH scenario to achieve an efficient absorption rate is the minority heating: in the case ofIGNITOR, when a small fraction of He (2-3%) is added to the D − T mixture, the first passabsorption on the ions near the center of the plasma column is very efficient. The remainingcoupled power is damped on the electrons over a broad radial interval of the plasma column.The fundamental harmonic of He is located in a radial interval − . < r/a < +0 . a denoting the minor radius of the tokamak) when the magnetic field is varied from 9 to13 T and the antenna frequency is f = 115 M Hz . For example, calculations of the powerabsorption level for three different external magnetic field (11, 12 and 13 T ), correspondingto three different times of the discharge evolution, for a plasma formed by D (50%), T(50%) and a small fraction of He ( (cid:39) − × m − , and from 4 to 6 KeV respectively, show that the peak of absorption is locatedat the fundamental harmonic of He and second harmonic of Tritium. In Fig. 3 the RF powerdeposition for a coupled power of 1 . M W is shown vs ρ when the magnetic field is rampingup from 11 to 13 T . It is possible to observe that the deposition is mainly concentratedat the fundamental harmonic of He , a small fraction being given to the electrons viaLandau damping depending on the minority fraction; moreover, when the field increases,the resonance layer moves towards the periphery, but still remaining in the bulk of the plasmaat 13 T . The power absorbed by the He (minority heating) is quasi- linearly redistributedon the collisional time essentially to the Deuterium and Tritium bulk ions, with a fractionto the electrons. The consequence is that the plasma temperature increases accelerating theattainment of ignition. An analytical expression for the power profiles inside the plasmacan be deduced by fitting the numerical results giving a Q add = Q ICRH ( r ) that is essentially7 =11TB=12TB=13T Q ICRH (MW/m ) r ê a5101520 Q ICRH H MW ê m L ! r/a FIG. 1:
RF power deposition (on the minority He (2%)) in M W/m vs r/a when the magneticfield is ramping up from (black line) to (red line) Tesla for the Ignitor plasma parameters inthe ignited scenario. The applied frequency is M Hz . independent on the bulk temperature Q add = β exp (cid:104) α B ( ρ ICRH ) B (cid:105) exp (cid:104) − ( ρ − ρ ICRH ) ∆ (cid:105) (11)The expression B/B may be estimated by adopting the validity of the standard model.We get [see Eq. (1)] B ( ρ ICRH ) /B (cid:39) / (1 + (cid:15)ρ ICRH ). The expression in Eq. (11) fits verywell the numerical curve (obtained by running TORIC+SSQLFP) by setting α = 15 . β = 6 . × − M W/m and ∆ = 0 . B = 13 T ( ρ ICRH = 0 . P ICRH , injected into the plasma P add = P ICRH = (cid:90) dV Q add = β (cid:90) dV exp (cid:104) α B ( ρ ICRH ) B (cid:105) exp (cid:104) − ( ρ − ρ ICRH ) ∆ (cid:105) (cid:39) . M W (12)which coincides with the power input in the numerical code TORIC.
III. EVALUATION OF THE THERMAL LOSS
The energy balance equation (6) should be completed with the transport equations relatingthe averaged thermodynamic flows < q e > and < q i > with the thermodynamic forces − T − e ∂ r T e and − T − i ∂ r T i . The complete transport relations are composed by the sum of8hree terms: the classical, the Pfirsch-Schl¨uter and the banana contributions. However inthis work, at the first and simplest approximation, we shall study the case where the closureequations are provided solely by the classical term, appropriately estimated for a plasma ina toroidal geometry. The general situation, where all the transport contributions are takeninto account, will be subject of a future work. Under this approximation, by kinetic theorywe find that, for a plasma in toroidal geometry, the averaged total heat flow is related tothe temperature gradient by the following equation [12] < q tot > = − (cid:88) ι = i,i (cid:16) < κ er ( T ) > CL + < κ ir ( T ) > CL (cid:17) ∂T∂r (13)where we have taken into account that T − i d r T i = T − e d r T e . The expression for the electronand ion thermal conductivities < κ ιr > ( ι = e, i ) can be brought into the form < κ ιr > CL = 52 n ι T ι m ι < ˜ κ ιr > CL (14)where τ ι and < ˜ κ ιr are the collision time and the dimensionless averaged thermal conductivityof the species ι , respectively. By kinetic theory we know that, in toroidal geometry, theclassical contribution to the transport coefficients coincides exactly with the asymptotic limitof the perpendicular transport coefficients estimated by the classical theory, averaged overa magnetic surface [12]. For asymptotic limit we mean the value of the classical transportcoefficients, estimated for Ω ι τ ι >> ι and τ ι denoting the Larmor frequency andthe collision time of species ι , respectively. In other terms, < κ er > CL = c e < e τ e ) > ; < κ ir > CL = c i < i τ i ) > = (cid:16) m i m e (cid:17) / c i c e < ˜ κ er > CL (15)Here, c ι are the (dimensionless) coefficients of the linear collision matrix of species ι , i.e. , c e = (13 + 42) /
10 (with Z = 1) and c i = 2 √ /
3. Index i in Eq. (15) stands for the effectiveion with mass m i = ( m D + m T ) / Z = 1. Note that, Eq. (15) takes into account the toroidal geometry of the Tokamak, but not the inhomogeneity and curvature of the magneticfield. As known, the latter is matter of the neoclassical theory. Moreover, in Eqs (14)and (15), we have taken into account that the nonlinear corrections to the linear classicaltransport coefficients may be neglected [18], [19] and [20]. At the steady state, we get − r ddr r (cid:16) < κ tot ( T ) > CL ddr T (cid:17) = Q α + Q b + Q add (16)9 IG. 2:
Surface-averaged, total classical thermal conductivity coefficient for IGNITOR plasmas , < κ tot > ( cm − sec − ) vs the minor radius and temperature . where < κ tot ( T ) > CL = < κ e ⊥∞ > + < κ i ⊥∞ > , and Q α , Q b and Q add are given by Eqs. (8),(10) and (11), respectively. Equation (16) is the simplest version of the steady state powerbalance where the terms corresponding to the loss of energy density due to the expansionof the fluid and to the loss of energy density due to diffusive processes are neglected. Thesesimplifications are justified from the fact that a magnetic fusion reactor is (almost) a steadystate system with small and negligible flows. In addition, we are dealing with stronglymagnetized plasmas and the turbulent effects are, therefore, notably reduced. Hence, at thefirst approximation, the time derivative as well as the convection and compression terms mayreasonable be neglected. Fig. (2) shows the total thermal conductivity (sum of the electronand the reduced ion thermal conductivities). In Figs (3) and (4) the total source-densityare reported, against the minor radius and temperature, in absence ( P ICRH (cid:39)
0) and inpresence of additional source ( P ICRH (cid:39) . M W ), respectively.Finally, Eq. (16) results to be a highly non-linear second order ordinary differential equationin the radial variable r , submitted to the boundary conditions for the equilibrium tempera-ture. The equilibrium temperature has been obtained by solving Eq. (16) numerically withthe following conditionsand d r T | r =0 = 0 T | r = a = pedestal = 20 eV (17)The first condition derives from the symmetry T ( r ) = T ( − r ) close to the center of theTokamak (we assume that, at r = 0, the derivative of the temperature exists and does notdiverge), and the choice of edge temperature T = 20 eV is reasonable for several types ofIgnitor L-mode plasmas. Figures (5) and (6) show the temperature profiles, against the10 IG. 3:
Total source profile versusthe minor radius and temperature in ab-sence of additional source ( P ICRH (cid:39)
Total source profile ver-sus the minor radius and tempera-ture in presence of ICRH ( P ICRH (cid:39) . M W ).
10 20 30 40 r (cid:72) cm (cid:76) (cid:72) Kev (cid:76)
FIG. 5:
Equilibrium Temperature profile whenno power RF is provided. This solution has beenobtained by solving numerically the steady stateenergy balance equation, with ( P ICRH (cid:39) submitted to the boundary conditions Eqs (17).
10 20 30 40 r (cid:72) cm (cid:76) (cid:72) Kev (cid:76)
FIG. 6:
Equilibrium Temperature profile whenICRH power is injected in the plasma ( P ICRH (cid:39) . M W ). This solution has been obtained bysolving numerically the steady state energy bal-ance equation, submitted to the boundary con-ditions Eqs (17). minor radius r , without RF power ( P ICRH (cid:39) P ICRH (cid:39) . M W ), respectively. Figures (7) and (8) illustrate the sourceprofiles against the normalized radius ρ for temperatures when no power RF is providedand when the power ICRH is injected into the plasma, respectively. These profiles havebeen obtained by inserting the equilibrium temperatures into the r.h.s. of Eq. (16), when11 .2 0.4 0.6 0.8 1.0 Ρ(cid:45) (cid:72) MW (cid:144) (cid:130) (cid:76) FIG. 7:
Total source profile versus the normal-ized minor radius when ( P ICRH (cid:39) Ρ (cid:72) MW (cid:144) m (cid:76) FIG. 8:
Total source profile versus the normal-ized minor radius when ( P ICRH (cid:39) . M W ). Q add = 0 and Q add (cid:54) = 0, respectively. Figs. (5) (6) (7), and (8), identify the region of theplasma where the profiles are not negligible (here referred to as the core of the plasma ).This region ranges from 0 ≤ r ≤ r , with r (cid:39) , cm ) and r (cid:39) , cm ) in absenceand in presence of ICRH (with P ICRH (cid:39) . M W ) respectively. Note that r is solution ofthe equation T ( r ) = ¯ T ≡ V − (cid:82) T dV = (2 /a ) (cid:82) a xT ( x ) dx . In the core of the Tokamak,the average temperature ¯ T core , defined as ¯ T core = V core (cid:82) core T dV = (2 /r ) (cid:82) r xT ( x ) dx , are¯ T core = 1 . KeV and ¯ T core = 6 . KeV , for P ICRH = 0 and P ICRH (cid:39) . M W , respectively.Hence, in absence of ICRH, the average temperature in the core of the plasma does notreach the desired ignition temperature, which, as known, should be ¯ T core > . KeV . Theadditional ICRH power significantly increases the equilibrium temperature in the core of theTokamak allowing the plasma to reach the ideal ignition temperature.
IV. CALCULATION OF THE STATIONARY ENERGY CONFINEMENT TIME
As known, the energy confinement time is defined as ratio between the total thermal energy W e in the plasma over the total energy rate through the boundary Γ E | boundary i.e. , τ E = W e Γ E | boundary (18)with W e = 3 V − (cid:82) nT dV and Γ E | boundary = V − (cid:82) ∇ · q tot dV . In our opinion, the mostconvenient way to estimate the energy confinement time correctly is to follow the procedureindicated by Freidberg [2]. The calculation begins by recalling the steady state 0 − D plasma12ower balance relation. For this, let us reconsider the stationary power balance relation,Eq. (19) n < σv > D − T E α − c B n T / + Q add = ∇ · q tot (19)where we assumed that the particle density is constant, n = n = const. , and c B , denotesthe Bremsstrahlung constant [see Eq. (10)]. The 0 − D steady state equation is obtainedby assuming that the temperature profile is constant across the plasma cross section withmagnitude equal to its average value ¯ T = V − (cid:82) T dV . We get[23]Γ E | boundary = Γ E | r = r = V − (cid:90) ∇ · q dV = n E α < σv > D − T ( ¯ T ) − c B n ¯ T / + Q add = W e τ E = 3 n ¯ Tτ E (20)where Eq. (18) has been used for the definition of the energy confinement time, and n ≡ n ( r )and Q add ≡ Q add ( r ) , respectively (recall that r = 25 . cm and T ( r ) = ¯ T ). Hence, τ E = 12 n ¯ TE α n < σv > D − T ( ¯ T ) − c B n ¯ T / + 4 Q add (21)Finally, in presence of ICRH, with ( P ICRH (cid:39) . M W ), the estimated energy confinementtime is τ E (cid:39) . sec . It should be stressed that the analysis is underestimating the diffusivelosses and this result should be regarded as just lower limits to the energy loss and hence anupper limit to the energy confinement time. Fig. (9) shows the profiles of pτ E against theminor radius, in presence of ICRH, with ( P ICRH (cid:39) . M W ). In conclusion, the additionalheating is required during the startup transient phase in order to heat the plasma from itslow initial temperature to the desired ignition temperature.
V. ANALYSIS OF THE THERMONUCLEAR INSTABILITY
The study of the stability of the solution, based directly on Eq. (6), is quite complex.In the Appendix, we report a simplified analysis of stability where calculations are greatlysimplified by eliminating the appearance of < γ er > and by assuming that the profiles (excepttemperature) are flat. These approximations will enable to determine the modes unstable.In this section we shall proceed in an even simpler way. We shall analyze the stability of thesolution in the core of the plasma through the time-dependent form of 0 − D power balanceequation and by exploiting the (approximately) uniformity, in the core of the Tokamak,13 .3 atm sec Critical Lawson ValueIGNITOR10 20 30 40 r (cid:72) cm (cid:76) Τ E (cid:72) atm sec (cid:76) FIG. 9:
The Lawson variable, pτ E , against the normalized minor radius for IGNITOR-plasmasubject to ICRH. At the core of the plasma, the minimum values of T and pτ E , required to satisfythe Lawson criterion for ignition, i.e., T ≥ . KeV and pτ E ≥ . atm sec , are attained (see alsoFig. (6). of the density profile. This approach will provide with the desired indications on thermalstability of the solution in the core of the plasma.The space-time plasma dynamical equation, expressing the conservation energy relation forIGNITOR plasma in the standard model, can be brought into the form3 ∂∂t ( nT ) = Q α + Q b + Q add − Q κ (22)where Q κ ≡ r ∂∂r (cid:2) r ( < q e > + < q D > + < q T > ) (cid:3) By averaging Eq. (22) over space, we obtain the corresponding time dependent 0 − D powerbalance equation: 3 n ddt T = S α + S b + S add − S κ (23)where S ξ ≡ V − core (cid:82) core dV Q ξ (with ξ = α, b, add, κ ) and we have assumed that n = n core = n (cid:39) const. The goal is now to examine the time dependance of a small perturbation δT ( t ) of the equilibrium temperature ( i.e. , T ( t ) = ¯ T core + δT ( t ) with δT ( t ) / ¯ T core << n = n = const. and S add is a fixedquantity independent of temperature (i.e., dS add /dT = 0) [2]. At the leading order, a smallperturbation δT ( t ) satisfies the evolution equation : dδTdt = 112 E α n (cid:16) ddT < σv > D − T − ˜ c b T − / − E α τ E (cid:17) T = ¯ T core δT (24)14 Α (cid:61)
15 10 15 20 25 30 35 T core (cid:72) Kev (cid:76) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) Λ FIG. 10:
Generic behavior of the critical eigen-value λ versus the average temperature (in ourcase, ¯ T core ), at the ignition value f α = 1. Thedashed line corresponds to λ profile estimated byneglecting the Bremsstrahlung radiation. T core (cid:61) f Α (cid:45) (cid:180) (cid:45) (cid:45) (cid:180) (cid:45) (cid:45) (cid:180) (cid:45) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) (cid:180) (cid:45) Λ FIG. 11:
Critical eigenvalue λ against f α , at ¯ T core = 6 . keV . The dashed line corre-sponds to the estimation of λ by neglectingthe Bremsstrahlung effect . where ˜ c B ≡ c B /E α . The stability condition can be further simplified by considering that(at linear order) the critical eigenvalue should be estimated at equilibrium . To this end,we recall that at the equilibrium temperature S κ = S α + S add and S add = [(1 − f α ) /f α ] S α .Hence, we find S κ = S α /f α , with f α denoting the fraction f α of the total heating power i.e. , f α = S α / ( S α + S add ) [2] (so, f α = 1 corresponds to ignition and f α = 0 to no α -power [24]).Now, by taking into account that S κ | T = ¯ T core = n ¯ T core /τ E and S α = n E α < σv > D − T / / ( E α τ E ) = < σv > D − T / ( f α ¯ T core ) and Eq. (24) finally simplifies to dδTdt = 112 E α n (cid:16) ddT < σv > D − T − ˜ c b T − / − < σv > D − T T f α (cid:17) T = ¯ T core δT (25)Hence, the solution is stable if λ ( ¯ T core , f α ) < λ ( ¯ T core , f α ) >
0, where λ ( ¯ T core , f α ) ≡ ddT < σv > D − T T = ¯ T core − ˜ c B ¯ T / core − f α ¯ T − core < σv > D − T T = ¯ T core (26)Fig. (10) reports on the generic profile of the critical eigenvalue λ against the average temper-ature (in our case, ¯ T core ) at the ignition value f α = 1. Fig. (11) shows the critical eigenvalue λ against the fraction of the total heating power, f α , estimated at ¯ T core = 6 . KeV . Thedashed lines refer to the critical eigenvalue estimated by neglecting the Bremsstrahlung ra-diation. In line with our expectations, at ¯ T core = 6 . KeV and f α = 1 (ignition), the coreof the plasma is unstable. The Bremsstrahlung effect provides a negligible contribution.15 I. CONCLUSIONS
One of the objectives of this work is to study in detail the equilibrium, and stability proper-ties, of the temperature evolution in burning fusion IGNITOR-plasma, in presence of ICRH.Although in this respect many manuscripts already appeared in literature, our aim is todetermine the equilibrium temperature and to study the stability of the solution, by makinguse of the plasma dynamical equations and transport relations, which are rigorously obtainedby kinetic theory. In addition, our approach gives some new insights concerning the thermalrunaway problem and, in particular, the relation between the 0 − D and 1 − D models. Herea scenario is considered where IGNITOR is led to operate in a slightly sub-critical regime byadding a small fraction of He to the nominal 50%-50% Deuterium-Tritium mixture. At thefirst step, we considered the simplest case where the transport coefficients are determinedby kinetic theory applied to classical plasma in a toroidal geometry. The obtained resultsmay be sketched as follows. i) We determined the temperature equilibrium profile solely by kinetic theory i.e., withoutthe auxilium of ad hoc models for the transport coefficients; ii)
We showed that in the core of the plasma the thermal solution is unstable, and weestimated that the value of the confinement time is τ E (cid:39) . sec. ; iii) The additional heating,ICRH, is required during the start up transient phase in orderto heat the plasma from its low initial temperature to the desired ignition temperature; iv)
We showed that the ICRH heating in the IGNITOR experiment is expected to triggerthe thermonuclear burning by means of the RF coupled power. The use of the ICRHcan be switched on and off along with the plasma parameter evolution and in particularwith the temperature. If we apply ICRH to a plasma, characterized by a subcriticalignition regime, we have shown that it is possible to trigger a thermal instability, byswitching off the ICRH the regime can be recast to a subcritical one. This meansthat in the subcritical regime the difference between power lost and alpha heating iscompensated by additional ICRH heating, which should be able to increase the globalplasma temperature via collisions between He minority and the background D − T ions. 16t is well known that a realistic estimation of energy confinement time should account the tur-bulent contributions and, in particular, the strong anomalous diffusion in the outer plasma.However, this is a very complex task. Our analysis is obviously valid in the core of theplasma, corresponding to the region 0 ≤ r ≤ r with T ( r ) = ¯ T . The clear identificationof the core of the plasma is due to our choice of boundary conditions, i.e. , a pedestal tem-perature with a thermal derivative that is zero at the edge. In the core the presence of theauxiliary ICRH heating is responsible of the triggering of the instability.We mention another aspect concerning the instability problem. Here, we have considered anIGNITOR (type) device, characterized by a large B -field and small dimensions. Of course,these conditions enormously simplified calculations, since a large- B field tends to freeze turbulent effects. In IGNITOR, indeed, we have evaluated that the non-linear (turbulent)contribution to the transport is not dramatic owing to the fact that IGNITOR operateswith a very strong external magnetic field. The magnetic field has a stabilizing effect on theturbulence. Evidence of this fact can be deduced by the calculation done by means of theTFT code [18] where an evaluation of the strength of the non linear contribution has beenestablished for the electron and ion fluxes. The result is that the difference is sufficiently weakand the linear theory can be used safely. In addition, the peculiarity of IGNITOR is that,since this reactor works at sufficiently low temperature (positive slope of temperature curve),the instability can develop as soon as the criteria of ignition are met. Another argument isthat in IGNITOR the collisionality regime is essentially banana for most of the dischargeradius being of the Pfirsch-Schl¨uter type only in a small portion at the edge and at the center.However we decided to dedicate a publication per se to study in a deeper manner all thetransport regimes of IGNITOR by covering the various radial zones. Other Reactor tokamakdesigns, based on low- B field and large dimensions have also been analyzed, but only interms of heuristic experimental scalings. This obviously does not hold for Tokamak Reactorlike ITER or DEMO, where the turbulence can play a crucial role in the determination oftransport coefficient. The main difference with respect to IGNITOR is that DEMO (which ischaracterized by low magnetic field, large dimensions, and very high temperature) is far fromdeveloping a thermal instability. DEMO in fact is characterized by a negative slope of thetemperature curve and for this reason is thermally stable [2]. In addition, for these reactors,a realistic estimation should take into account the strong anomalous diffusion in the outerplasma and, under this conditions, the temperature profiles estimated by using the classical17hermal diffusion would probably result in highly unrealistic shapes. Anyhow, a deeperanalysis in which a comparison between both approaches (high field and low temperature)and low field and large dimension will be performed more extensively in a dedicated work.We would like also clarify another crucial point: the role of accumulation of reaction ashes He , which may eventually quench the thermonuclear process. In reality, in this work weconsidered the emergence and development of the thermal instability just at the end of theflattop. In this scenario, the presence of the alpha particle is still too low to give someevaluable effect on the dynamic of the reaction. Obviously, during the flattop, the presenceof a consistent fraction of He could induce the quench of the thermonuclear reaction belowthe useful threshold. In fact the ashes play the same role of the impurities by unbalancingthe good ratio of the reactant (50% Deuterium and 50% Tritium). Also in this case, in ouridea, the ICRH power turns to be a useful tool in giving a boost at the plasma temperatureto compensate the presence of the impurities that are degrading the reaction rate.Now, we should proceed step-by-step. In the next step we shall consider the general situation,where the transport coefficients are determined by considering all the collisional transportregimes (ı.e., the classical, Pfirsch-Schl¨uter and banana transport regimes), and the nonlinearcontributions are no longer neglected. The results will be compared successively with thesolutions obtained by using a turbulent transport model, like the gyro-Bohm model. All ofthis will be subject of future works. VII. ACKNOWLEDGEMENTS
One of us (GS) is indebted to Gy¨orgy Steinbrecher, of the University of Craiova (Roma-nia), for the fruitful discussions concerning the topic presented in the Appendix. GS is alsovery grateful to Alberto Sonnino, of the Karlsruhe Institute of Technology (KIT) - Germanyand the Universit´e Catholique de Louvain (UCL) - Belgium, for his assistance in performingnumerical calculations.
Appendix A: Determination of the Modes Unstable - Simplified Calculations.
In this Section, we report a quite simplified analysis of the stability of the equilibriumtemperature, showing the methodology allowing the determination of the modes unstable. A18emi-quantitatively accurate approximation is to assume that all the profiles, except temper-ature, and the coefficients are flat. This is not a very good approximation because, actually,the profiles and the coefficients (in particular) are not flat, but the approximation greatlysimplified the analysis. Hence, this Appendix should be understood only as an example ofcalculation with a view to illustrating the procedure.Let us put T ( r, t ) = T ( r ) + δT ( r, t ) where T ( r ) is the equilibrium temperature, solution ofEq. (19), and δT ( r, t ) the temperature perturbation. From Eqs (6), (13) and (19), and takinginto account the expression of S gain − loss , we find the evolution equation for the perturbation δT ( r, t )3 n ∂δT∂t − < κ tot > T = T ( r ) ∂ δT∂r = − n (cid:16) c B T ( r ) − / − E α ∂∂T < σv > T = T ( r ) (A1) − n ∂ T ∂r ∂∂T < κ tot > T = T ( r ) (cid:17) δT where we have taken into account that ∂ T Q add = 0 and we have neglected terms higher thanthe first order in δT . In addition, we have supposed that the contribution r − ( ∂ r T ) ∂ r ( r <κ tot > ) may be neglected with respect to < κ tot > ∂ T∂r . We assume now that the thermalconductivity < κ tot > coefficient and − n (cid:0) c B T ( r ) − / − E α ∂ T < σv > | T = T ( r ) (cid:1) are constantand estimated at the average temperature T = ¯ T , i.e. < κ tot > T = T ( r ) = < κ tot > T = ¯ T = const. − n (cid:16) c B T ( r ) − / − E α ∂∂T < σv > T = T ( r ) (cid:17) = (A2) − n (cid:16) c B ¯ T − / − E α ∂∂T < σv > T = ¯ T (cid:17) = const. In order to have an idea on the validity of first approximation in Eq. (A2), we report inFigs (12) and (13) the total average thermal coefficient < κ tot > versus the normalizedminor radius ρ , in absence and in presence of ICRH, respectively. These profiles have beenobtained by putting the equilibrium temperature-profiles, given in Figs (5) and (6), into < κ tot > CL , respectively. Hence, Eq. (A1) takes the form ∂∂t δT = D ∂ ∂r δT + βδT (A3)with D ≡ < κ tot > | T = ¯ T n and β = − n (cid:16) c B ¯ T − / − E α ∂∂T < σv > T = ¯ T (cid:17) (A4)19 .2 0.4 0.6 0.8 1.0 Ρ (cid:180) (cid:180) (cid:180) (cid:180) (cid:60) K tot (cid:62) FIG. 12:
Total average thermal conduc-tivity (measured cm − sec − ), against thenormalized minor radius, in absence of ad-ditional sources. This profile has been ob-tained by putting the equilibrium tempera-ture given in Fig. (5) into < κ tot > . Ρ (cid:180) (cid:180) (cid:180) (cid:180) (cid:60) K tot (cid:62) FIG. 13:
Total average thermal conduc-tivity (measured cm − sec − ), against thenormalized minor radius, in presence ofICRH, with ( P ICRH (cid:39) . M W ). Thisprofile has been obtained by putting theequilibrium temperature given in Fig. (6)into < κ tot > . The boundary conditions may be determined by imposing that both temperature and itsderivative do not fluctuate at the boundary. So, we have to solve Eq. (A1) subject to ddr δT r =0 = 0 and δT r = a = 0 (A5)By setting δT ( r, t ) = e − ωt f ( r ), we get D d f ( r ) dr + ( β + ω ) f ( r ) = 0 (A6)with d r f | r =0 = f | r = a = 0. The solution of Eq. (A6) can be brought into the form f ( r ) = (cid:80) nk =0 ˆ f k cos( kr ). We find ˆ f = 0 (for k = 0) − Dk + ( β + ω ) = 0 (for k (cid:54) = 0) (A7)The boundary conditions (A5) provide the relation between ω and the modes n . Indeed,cos( ka ) = 0 = ⇒ ka = π nπ ( n = ± , ± , · · · ) (A8)By substituting Eq. (A8) into Eq. (A6), we get ω ( k ) = − β + Dk ⇒ ω ( n ) = − β + D (cid:16) πa (cid:17) (cid:16) n + 12 (cid:17) (A9)20y taking into account Eq. (A4), we find that the modes unstable satisfy the inequality (cid:16) n + 12 (cid:17) a n π < κ tot > | T = ¯ T (cid:16) E α ∂∂T < σv > T = ¯ T − c B ¯ T − / (cid:17) (A10)In particular, the Goldstone mode ( n = 0) is unstable if a n E α ∂∂T < σv > T = ¯ T − a n c B ¯ T − / − π < κ tot > T = ¯ T > ribbon modes have to be considered in order to envision and predict how a condition ofglobal ignition can be reached [22]. [1] B. Coppi, M. Nassi and L. E. Sugiyama, Physica Scripta , , 112 (1992).[2] J. P. Freidberg, Plasma Physics and Fusion Energy, Cambridge University Press , Cambridge,USA, (2007).[3] R. G. Mills,
The Problem of Control of Thermonuclear Reactors , Los Alamos report , LA-4250 ,B1-1-B1-5 (1969).[4] W. M. Stacey,
Fus. Science Techn. , , 29 (2007).[5] S. V. Putvinskii, Sov. J. Plasma Phys. , , 694 (1980).[6] S. G. Bespoludennov, V. I. Pistunovich, A. I. Mel’dianov and S. A. Galkin, Techniques ofthe fusion burn control, Proceedings of the 4th Technical Committee Meeting and Workshopon Fusion Reactor Design and Technology ; Yalta (USSR); 26 May-6 June 1986; InternationalAtomic Energy Agency, Vienna (Austria) (1987),
Fusion reactor design and technology , .IAEA-TC-392.3/38 (1986).[7] S. Migliori, A. Frattolillo, S.K. Combs, L.R. Baylor, G. Roveta, F. Bombarda, R. Foust, D. T.Fehling, J. M. McGill, J.B.O. Caughman and J.C. Thomas, The Compact Multiple Barrel HighSpeed Pellet Injector for the Ignitor Experiment, Proceedings of 21rst IEEE/NPS Symposiumon Fusion Engineering , Sept. 2005, IEEE, doi: 10.1109/FUSION.2005.252921 (2005).[8] J. Mandrekas and W. M. Stacey,
Evaluation of Different Control Methods for the ThermalStability of the International Thermonuclear Experimental Reactor , Fusion Technology , ,57 (1991).
9] L. Hartch, V. Fuchs, and A. Bers,
Nucl. Fus. , , 833 (1980).[10] Ya. I. Kolesnichenko, V.V. Lutsenko and S. N. Reznik, Fusion Technology , , 84 (1994).[11] A. Cardinali and B. Coppi, Bull. Am. Phys. Soc. , , 73 (2009).[12] R. Balescu, Transport Process in plasmas Vol. I & II , Elsevier Science Publication, North-Holland, (1988).[13] Xing Z. Li, Qing M. Wei and Bin Liu,
Nucl. Fusion , , 125003 (2008).[14] A. Cardinali, ENEA Technical Report RT/2009/37/FUS , Frascati, Italy (2009).[15] M. Brambilla,
Plasma Physics and Controlled Nuclear Fusion , , 1 (1999).[16] M. Brambilla, Nucl. Fusion , , 1121 (1994).[17] M. Brambilla and R. Bilato, Nucl. Fusion , , 085004 (2009).[18] G. Sonnino and P. Peeters, Physics of Plasmas , , 062309 (2008).[19] G. Sonnino, Phys. Rev. E , , 051126 (2009).[20] G. Sonnino and A. Sonnino, J. of Thermodynamics & Catalysis , doi: 10.4172/2157-7544.1000129 (2014).[21] B. Coppi, Tri-dimensional Ribbon Burning Modes in Igniting Plasmas , Bull. Am. Phys. Soc. , (cid:63)(cid:63)(cid:63) (2014).[22] B. Coppi, Comments on Plasma Phys. and Cont. Fus. , , 2 (1977).[23] Note that, for magnetically confined plasma in toroidal geometry, in the classical regime, thethermal diffusion strongly varies with the inverse of temperature. In particular, within thedomain 0 ≤ r ≤ a , the vector field is not of class C everywhere (see FIG. 2). As known,under this circumstance, the divergence theorem does not apply. When the vector field is notof class C , we should proceed with the direct estimation of the integral of the divergenceof the field over the volume: P Q = (cid:82) dV ∇ · q tot = (cid:82) dV ( P Q + P b + P ICRH ) (cid:54) = 0. In ourproblem, the correct way to perform calculations when q tot = − < κ > dT /dr is not of class C everywhere, is to proceed as follows: 1) Take into account the balance equation, Eq. (19),and 2) Follow the procedure indicated by Freidberg i.e. , to develop the integral by retainingonly the leading order, according to the equations: f ( T ) ≈ f ( ¯ T ) + f ( ¯ T ) δT + O ( δT ). Hence, V − (cid:82) f ( T ) dV ≈ f ( ¯ T ) + O ( δT ). By following the above-mentioned procedure, we arrive atEq. (21), which should be correct.[24] Note that the definition for ignition f α = S α / ( S α + S add ) is satisfied for S add = 0 at anytemperature and nonzero fusion power S α . It only makes sense when the temperature exceeds ome limit, such as the ideal ignition temperature T ≥ KeV , so S = S b ..