**Modelling of Global Heat Transfer by Local Cooling in Fusion Plasmas in Fusion Plasmas**

**Modelling of Global Heat Transfer by Local Cooling**

Mikhail Tokar Additional information is available at the end of the chapter

Mikhail Tokar

Additional information is available at the end of the chapter 10.5772/60642

http://dx.doi.org/10.5772/60642

### **1. Introduction**

Global reaction on very localized but intensive disturbances is a problem of interest for diverse physical media and systems. In devices for magnetic nuclear fusion such conditions arise, e.g., if a significant inflow of neutral particles, either of the working hydrogen isotopes or of impurities, enter the plasma. These inflows can occur both spontaneously and be created deliberately. Even particles enter the plasma through spots with dimensions much smaller than those of the plasma boundary their density is so high that the local plasma can be perturbed very substantially. One of the most important channels for perturbations is the local cooling of the plasma through the energy losses on the excitation and ionization of inflowing neutral particles. Although such distortions are triggered very locally, they can spread in the plasma far from the position of their origin. Thus, the heat conduction along magnetic surfaces towards the region, where the plasma is directly cooled down by intruding particles, may decrease the temperature in remote plasma. The drop of the pressure, induced by cooling in the neutral cloud, triggers plasma flow resulting in a reduction of the plasma density far away from the perturbation source. If neutral particles are injected for diagnostic purposes [1] it is desirable to keep distortions in the plasma as small as possible, to get trustful information about the original plasma state. In other cases, e.g., by a massive gas injection (MGI) [2], used to mitigate plasma disruptions in Tokamak devices [3], the impact has to be maximized. Generally, the changes induced in the plasma parameters are determined both by the duration of the injection and by the number of injected particles. Depending on the application type these characteristics change in extremely broad ranges. Thus by laser-induced ablation spectroscopy, LIAS [4], an intense laser pulse irradiates a small part of the device wall during an extremely short time in the nanosecond-range. By measuring the intensity of impurity line emission in the plasma one can assess the amount of particles emitted by the radiation, normally of the order of 1016−17, and judge about the wall composition. This diagnostics is designed as a passive non-invasive one and negligible

the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2015 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

©2012 prezimena autora, kod vise prvi et al., licensee InTech. This is an open access chapter distributed under

modifications in the background plasma are expected by LIAS. In MGI experiments up to 1023−<sup>24</sup> impurity neutral particles are injected through valves if special magnetic coils record precursors of a disruption. This results in a strong cooling of the whole plasma through radiation during several *ms* and the generation of relativistic electrons, being extremely dangerous for the wall, can be avoided. Besides a deliberate injection strong localized sources of neutrals can arise in fusion plasmas spontaneously, often if some critical phenomena, limiting the discharge performance, develop. For example at the so called "density limit" where the confinement of energy in the plasma deteriorates suddenly by approaching to some critical density [5], inflows of neutrals from certain wall spots increase tremendously. Multifaceted asymmetric radiation from the edge (MARFE) is a prominent example of such phenomena [6].

Understanding, modelling and prediction of the behaviour of plasmas affected by intense localized inflows of neutral particles is of importance for designing of injection systems and for controlling of conditions critical for performance. This is, however, a very demanding task, including a self-consistent description of non-linear processes of extremely different scales in time and space, from *ns* and *mm* till *s* and *km*, respectively. A consistent approach has to be based on time-dependent three-dimensional simulations. Presently there are, however, principal limitations to perform such calculations, mostly because various physical mechanisms are non-linearly interrelated and parallelization of time-dependent calculations is impracticable. Estimates [7] show that a realization of such a comprehensive approach to obtain reasonably accurate solutions during a time shorter than a year is still out of possibilities even of the most modern computers. To overcome technical problems outlined and to obtain an approximate but accurate enough model for the phenomena in question a new reduced approach has been elaborated [7–9]. It is based on reduction of three-dimensional fluid transport equations for particle, parallel momentum and energy transfer for all neutral and charged components of the plasma to a set of one-dimensional equations describing the time variation of the radial profiles for several parameters characterizing 3-D solutions although crudely but exhaustively enough.

For impurity spreading from a localized source it takes into account that the cross-sections by magnetic surfaces of the regions occupied by different neutral and charged species represent nested shells with dimensions increasing with the particle charge. This approach, named "shell model" [7, 8], allows assessing the characteristic density of the species and dimensions of their shells along and across magnetic field lines. By comparing with direct numerical solutions [8, 9], the approximations done in the "shell model" have been proven to be sufficiently accurate. By describing the main plasma components, a two zone approximation (TZA) has been elaborated presuming the separation of magnetic surfaces into cooled and hot zones. In the former one the energy is dissipated in direct interactions with injected particles, and the latter one is cooled down by the transfer, predominantly with the parallel heat conduction, to the cooled zone [9]. The radial profiles of the temperatures of plasma components in both zones are governed by one-dimensional heat balance equations interrelated non-linearly through the heat and particle fluxes between zones. A firm determination of these fluxes on the basis of two-dimensional transport equations [10] tremendously slows down numerical calculations. Therefore an assessment on the basis of a minimum entropy production principle [11] has been utilized in the TZA approach [9, 12]. It allows to relate analytically the heat fluxes between cooled and hot zones with the average temperatures in the zones. In addition to extreme savings in the

**Figure 1.** One-dimensional plasma stripe with an intense heat sink llocalized in the small segment of the stripe *Lx* − *δ<sup>c</sup>* ≤ *x* ≤ *Lx*.

CPU time, this approach allows to implement straightforwardly the so called heat flux limit [13], restricting the heat flux to a fraction of the free-streaming one if the temperature decay length becomes small compared to the mean free path length of electrons. The present paper is mostly devoted to a sophisticated test of the models being on the ground of the TZA, by comparing its predictions with the results of direct numerical solutions in oneand two-dimensional configurations. As examples for application the plasma edge cooling by argon MGI into JET and plasma distortions by LIAS in TEXTOR are modelled and the results of calculations are compared with the experimental data.

#### **2. One-dimensional configuration**

#### **2.1. Numerical solution**

2 Heat Transfer Studies and Applications

phenomena [6].

modifications in the background plasma are expected by LIAS. In MGI experiments up to 1023−<sup>24</sup> impurity neutral particles are injected through valves if special magnetic coils record precursors of a disruption. This results in a strong cooling of the whole plasma through radiation during several *ms* and the generation of relativistic electrons, being extremely dangerous for the wall, can be avoided. Besides a deliberate injection strong localized sources of neutrals can arise in fusion plasmas spontaneously, often if some critical phenomena, limiting the discharge performance, develop. For example at the so called "density limit" where the confinement of energy in the plasma deteriorates suddenly by approaching to some critical density [5], inflows of neutrals from certain wall spots increase tremendously. Multifaceted asymmetric radiation from the edge (MARFE) is a prominent example of such

Understanding, modelling and prediction of the behaviour of plasmas affected by intense localized inflows of neutral particles is of importance for designing of injection systems and for controlling of conditions critical for performance. This is, however, a very demanding task, including a self-consistent description of non-linear processes of extremely different scales in time and space, from *ns* and *mm* till *s* and *km*, respectively. A consistent approach has to be based on time-dependent three-dimensional simulations. Presently there are, however, principal limitations to perform such calculations, mostly because various physical mechanisms are non-linearly interrelated and parallelization of time-dependent calculations is impracticable. Estimates [7] show that a realization of such a comprehensive approach to obtain reasonably accurate solutions during a time shorter than a year is still out of possibilities even of the most modern computers. To overcome technical problems outlined and to obtain an approximate but accurate enough model for the phenomena in question a new reduced approach has been elaborated [7–9]. It is based on reduction of three-dimensional fluid transport equations for particle, parallel momentum and energy transfer for all neutral and charged components of the plasma to a set of one-dimensional equations describing the time variation of the radial profiles for several

parameters characterizing 3-D solutions although crudely but exhaustively enough.

For impurity spreading from a localized source it takes into account that the cross-sections by magnetic surfaces of the regions occupied by different neutral and charged species represent nested shells with dimensions increasing with the particle charge. This approach, named "shell model" [7, 8], allows assessing the characteristic density of the species and dimensions of their shells along and across magnetic field lines. By comparing with direct numerical solutions [8, 9], the approximations done in the "shell model" have been proven to be sufficiently accurate. By describing the main plasma components, a two zone approximation (TZA) has been elaborated presuming the separation of magnetic surfaces into cooled and hot zones. In the former one the energy is dissipated in direct interactions with injected particles, and the latter one is cooled down by the transfer, predominantly with the parallel heat conduction, to the cooled zone [9]. The radial profiles of the temperatures of plasma components in both zones are governed by one-dimensional heat balance equations interrelated non-linearly through the heat and particle fluxes between zones. A firm determination of these fluxes on the basis of two-dimensional transport equations [10] tremendously slows down numerical calculations. Therefore an assessment on the basis of a minimum entropy production principle [11] has been utilized in the TZA approach [9, 12]. It allows to relate analytically the heat fluxes between cooled and hot zones with the average temperatures in the zones. In addition to extreme savings in the Our consideration begins with a one-dimensional model for the transport in a plasma stripe along the magnetic field oriented in the direction *x*, as it is shown in figure 1. This situation is typical for the scrape-off layer (SOL) in a limiter tokamak like TEXTOR [14]. We consider the reaction of the whole plasma stripe on a cooling which is suddenly initiated at time *t* = 0 in a small fraction of the stripe through, e.g., energy losses on the excitation of injected impurity. The plasma is characterized by the temperature *T* (*t*, *x*) and the deviation of *T* from its initial homogeneous level *T*<sup>0</sup> is governed by the following heat conduction equation

$$1.5\partial\_{\!\!t}T - \partial\_{\!\!x} \left(\chi \partial\_{\!\!x} T\right) = \nu\_{\!\!h} \left(T\_0 - T\right) - \nu\_{\!\!I} E\_{\!\!I} \tag{1}$$

Here *χ* = *κ*/*n*, with *n* being the plasma density and *κ* the heat conduction of plasma electrons; the latter is assumed here in the Spitzer-Härm approximation [15], i.e. *κ* = *AeT*2.5 with *Ae* = 1020*eV*−2.5*cm*−1*s*−1; the former term on the right hand side (rhs) is the heat source, maintaining *T* (*t* < 0, *x*) = *T*<sup>0</sup> in a steady equilibrium; the latter term is the localized energy loss activated at *t* ≥ 0 in the region *Lx* − *δ<sup>c</sup>* ≤ *x* ≤ *Lx* with *δ<sup>c</sup> Lx*; this is assumed dependent on the electron temperature, *ν<sup>l</sup>* = *ν*<sup>0</sup> exp (−*El*/*T*), to take into account that the excitation rate drops sharply if the electron temperature reduces below the excitation energy *El* of injected particles. The boundary condition correspond zero heat conduction through the borders of the computation domain, *∂xT* (*t*, *x* = 0) = *∂xT* (*t*, *x* = *Lx*) = 0.

Henceforth the temperature profile is characterized by its values averaged over the hot, 0 ≤ *x* ≤ *Lx* − *δc*, and cooled, *Lx* − *δ<sup>c</sup>* ≤ *x* ≤ *Lx*, zones of the plasma stripe:

$$\langle T \rangle\_{\hbar} = \int\_{0}^{L\_{\text{x}} - \delta\_{\text{c}}} T \left( t, \mathbf{x} \right) \frac{d\mathbf{x}}{L\_{\text{x}} - \delta\_{\text{c}}}, \\ \langle T \rangle\_{\text{c}} = \int\_{L\_{\text{x}} - \delta\_{\text{c}}}^{L\_{\text{x}}} T \left( t, \mathbf{x} \right) \frac{d\mathbf{x}}{\delta\_{\text{c}}} \tag{2}$$

**Figure 2.** The time evolution of the temperature values averaged over the hot (a,c) and cooled (b,d) zones, computed by solving equation (1) numerically for *ν<sup>h</sup>* = 104*s*−<sup>1</sup> (a,b) and *ν<sup>h</sup>* = 102*s*−<sup>1</sup> (c,d), with different cooling rates *ν*0.

Figure 2 demonstrates the time variation of *T<sup>h</sup>*,*<sup>c</sup>* computed by solving equation (1) with numerical methods outlined in Ref.[16], for different combinations of the parameters *ν<sup>h</sup>* and *ν*0; other parameters involved into the problem are typical for the SOL plasma in TEXTOR with injected carbon impurity: *Lx* = 10*m*, *δ<sup>c</sup>* = 0.1*m* and *El* = 1.26*eV*. The main conclusion from these results is the existence for sufficiently large *ν*<sup>0</sup> of phase with sudden acceleration in the temperature drop in the cooled zone. This fact is related to a cooling instability predicted in references [17, 18] on the basis of zero-dimensional models and discussed deeper below.

#### **2.2. Two-zone approximation**

As it was claimed above the main aim of the present study is to develop reduced models to describe global impacts of local cooling in *three-dimensional* configurations. This would give a possibility to model phenomena in question without direct 3-D non-stationary calculations, being presently out of possibilities of even the most modern computers. An important

4 Heat Transfer Studies and Applications

0

100 6 7 T (eV) <sup>h</sup>

0

20

40

60

80

0 10 20 30 40 50

0.3 x109

0 10 20 30 40 t ( s) m 50

0.6 x109

1.8 x109

**2.2. Two-zone approximation**

n0=0.06 x109 -1 s

t ( s) <sup>m</sup> <sup>0</sup>

(c)

20

0

**Figure 2.** The time evolution of the temperature values averaged over the hot (a,c) and cooled (b,d) zones, computed by solving equation (1) numerically for *ν<sup>h</sup>* = 104*s*−<sup>1</sup> (a,b) and *ν<sup>h</sup>* = 102*s*−<sup>1</sup> (c,d), with different cooling rates *ν*0.

Figure 2 demonstrates the time variation of *T<sup>h</sup>*,*<sup>c</sup>* computed by solving equation (1) with numerical methods outlined in Ref.[16], for different combinations of the parameters *ν<sup>h</sup>* and *ν*0; other parameters involved into the problem are typical for the SOL plasma in TEXTOR with injected carbon impurity: *Lx* = 10*m*, *δ<sup>c</sup>* = 0.1*m* and *El* = 1.26*eV*. The main conclusion from these results is the existence for sufficiently large *ν*<sup>0</sup> of phase with sudden acceleration in the temperature drop in the cooled zone. This fact is related to a cooling instability predicted in references [17, 18] on the basis of zero-dimensional models and discussed deeper below.

As it was claimed above the main aim of the present study is to develop reduced models to describe global impacts of local cooling in *three-dimensional* configurations. This would give a possibility to model phenomena in question without direct 3-D non-stationary calculations, being presently out of possibilities of even the most modern computers. An important

20

1.8 x109

40

60

80

100 6 7 T (eV) <sup>c</sup> 1.8 x109

40

60

0.3 x109

(a) (b)

0 10 20 30 40 t ( s) m 50

0 10 20 30 40 t ( s) m 50

0.6 x109

0.3 x109

0.6 x109

80

n0=0.06 x109 -1 s

n0=0.06 x109 -1 s

(d)

100 6 7 T (eV) <sup>c</sup>

0.3 x109

0.6 x109

1.8 x109

n0=0.06 x109 -1 s

20

40

60

80

100

6 7 T (eV) <sup>h</sup>

**Figure 3.** The time evolution of the temperature values averaged over the hot (solid curves, crosses and circles) and cooled (dashed curves, boxes and triangles) zones, computed by solving equation (1) numerically (curves) and by using TZA equations (symbols) for *ν<sup>h</sup>* = 104*s*−<sup>1</sup> (a) and *ν<sup>h</sup>* = 102*s*−<sup>1</sup> (b), with different cooling rates *ν*0.

element of this complex of models is a two-zone approximation (TZA) of the transport on magnetic surfaces. The latter is dominated by transfer along magnetic field lines and therefore it is instructive to start with elaborating of such an approximation in the one-dimensional case. By integrating equation (1) over the hot and cooled zones we get the following ODE, governing the time evolution of *T<sup>h</sup>* and *T<sup>c</sup>*:

$$\left(1.5\partial\_{\hbar}\left\_{\hbar} + \frac{q\_{\hbar c}}{n\left(L\_{\chi} - \delta\_{c}\right)} = \nu\_{\hbar}\left(T\_{0} - \left\_{\hbar}\right),\tag{3}$$

$$1.5\partial\_{\!t} \left< T \right>\_{\mathcal{C}} - \frac{q\_{\!\!\!
th}}{n\delta\_{\!\!\!
\!
}} = \nu\_{\!\!\!
} \left( T\_{0} - \left< T \right>\_{\mathcal{C}} \right) - \left< \nu\_{\!\!\!
} \right>\_{\mathcal{C}} E\_{\!\!\!
} \tag{4}$$

Here *qhc* is the density of the heat flux conducted from the hot zone to the cooled one. To assess this approximately we apply approaches developed in reference [19] and proceed from equation (1) in the form:

$$
\partial\_{\chi} \left( \chi \partial\_{\chi} T \right) = Q \tag{5}
$$

with the term *Q* combining all other contributions in equation (1) besides the heat conduction one. The main assumption of the two-zone approximation is: in both zones *Q* can be approximated by some *a priory* unknown functions of time, *Qh* (*t*) and *Qc* (*t*), respectively, being **independent** of *x*. Under this assumption equation (5) can be analytically integrated for *<sup>χ</sup>* <sup>∼</sup> *<sup>T</sup>β*. By exploiting the continuity of *<sup>T</sup>* and *<sup>∂</sup>xT* at the zone interface, we get:

$$T\left(t,\ 0 \le \mathbf{x} \le L\_{\mathbf{x}} - \delta\_{\mathbf{c}}\right) = T\_h\left(t\right)\phi\_h\left(\xi \equiv \frac{\mathbf{x}}{L\_{\mathbf{x}} - \delta\_{\mathbf{c}}}\right) \tag{6}$$

$$T\left(t, L\_{\mathbf{x}} - \delta\_{\mathbf{c}} \le \mathbf{x} \le L\_{\mathbf{x}}\right) = T\_{\mathbf{c}}\left(t\right)\phi\_{\mathbf{c}}\left(\xi \equiv \frac{L\_{\mathbf{x}} - \mathbf{x}}{\delta\_{\mathbf{c}}}\right) \tag{7}$$

where *Th* ≡ *T* (*t*, *x* = 0), *Tc* ≡ *T* (*t*, *x* = *Lx*) and the functions *φ<sup>h</sup>* and *φ<sup>c</sup>* of 0 ≤ *ξ* ≤ 1 are given by the relations:

$$\phi\_h = \left[1 - \left(1 - \theta^{\theta+1}\right) \left(1 - \frac{\delta\_c}{L\_\chi}\right) \xi^2\right]^{1/(\theta+1)},\\\phi\_c = \left[1 + \left(1/\theta^{\theta+1} - 1\right) \frac{\delta\_c}{L\_\chi} \xi^2\right]^{1/(\theta+1)}\tag{8}$$

with *θ* = *Tc*/*Th*. The temperatures at the stripe edges, *Th* and *Tc*, have to be related to *T<sup>h</sup>* and *T<sup>c</sup>* according to the definitions (2). With the profiles obtained above we get

$$T\_{h\mathcal{L}} = \langle T \rangle\_{h\mathcal{L}} \,/\,\,\langle \Phi \rangle\_{h\mathcal{L}}$$

and *θ* is determined from the transcendent equation

$$\theta = \frac{\langle T \rangle\_c \langle \phi \rangle\_h \left( \theta \right)}{\langle T \rangle\_h \left\langle \phi \right\rangle\_c \left( \theta \right)}$$

Finally, for the Stirzer-Härm formula for the heat conduction, with *β* = 2.5:

$$q\_{\rm hc} = \frac{4A\_{\rm c}}{7} \frac{T\_{\rm h}^{3.5} - T\_{\rm c}^{3.5}}{L\_{\rm x}} \tag{9}$$

In figure 3 we compare the time variation of *T<sup>h</sup>* and *T<sup>c</sup>* calculated from equations (3) and (4) with those found by direct numerical integration of equation (1). The agreement between two approaches is nearly perfect, both for low and for high level of the cooling. Thus the relation (9) for the heat outflow from the hot zone to the cooled one provides a firm basis to describe local and global heat transfer induced by a strongly localized cooling. By concluding this subsection we discuss shortly the phenomenon of cooling instability. It develops since the heat conduction influx into the cooled zone becomes saturated with reducing *Tc*, see equation (9), and cannot compete with the local cooling there if *ν*<sup>0</sup> is large enough. The existence of a threshold for cooling instability is demonstrated in figure 4 by the time variation of *T<sup>c</sup>* computed in the TZA for different magnitudes of *ν*0: there is a characteristic change in the behaviour of the cooled zone temperature occurring by a slight change of *<sup>ν</sup>*<sup>0</sup> in the range 7.2 <sup>−</sup> 7.8 <sup>×</sup> 107*s*−1.

#### **2.3. Assessment of other approximations**

Consider now the role of physical effects neglected in equation (1): heat transfer by electron-ion coulomb collisions, particle convection caused by the pressure drop in the cooled plasma zone, the induced change in the plasma density, etc. For this purpose we use now the following system of equations, describing the evolution of the electron and ion temperature *Te* and *Ti*, respectively:

6 Heat Transfer Studies and Applications

given by the relations:

<sup>1</sup> <sup>−</sup> *<sup>θ</sup>β*+<sup>1</sup>

and *θ* is determined from the transcendent equation

change of *<sup>ν</sup>*<sup>0</sup> in the range 7.2 <sup>−</sup> 7.8 <sup>×</sup> 107*s*−1.

**2.3. Assessment of other approximations**

*Te* and *Ti*, respectively:

<sup>1</sup> <sup>−</sup> *<sup>δ</sup><sup>c</sup> Lx ξ*2

*φ<sup>h</sup>* = 1 −  *T* (*t*, *Lx* − *δ<sup>c</sup>* ≤ *x* ≤ *Lx*) = *Tc* (*t*) *φ<sup>c</sup>*

where *Th* ≡ *T* (*t*, *x* = 0), *Tc* ≡ *T* (*t*, *x* = *Lx*) and the functions *φ<sup>h</sup>* and *φ<sup>c</sup>* of 0 ≤ *ξ* ≤ 1 are

with *θ* = *Tc*/*Th*. The temperatures at the stripe edges, *Th* and *Tc*, have to be related to *T<sup>h</sup>*

*Th*,*<sup>c</sup>* = *Th*,*<sup>c</sup>* / *φh*,*<sup>c</sup>*

*<sup>θ</sup>* <sup>=</sup> *<sup>T</sup><sup>c</sup> <sup>φ</sup><sup>h</sup>* (*θ*) *T<sup>h</sup> φ<sup>c</sup>* (*θ*)

> *T*3.5 *<sup>h</sup>* <sup>−</sup> *<sup>T</sup>*3.5 *c Lx*

In figure 3 we compare the time variation of *T<sup>h</sup>* and *T<sup>c</sup>* calculated from equations (3) and (4) with those found by direct numerical integration of equation (1). The agreement between two approaches is nearly perfect, both for low and for high level of the cooling. Thus the relation (9) for the heat outflow from the hot zone to the cooled one provides a firm basis to describe local and global heat transfer induced by a strongly localized cooling. By concluding this subsection we discuss shortly the phenomenon of cooling instability. It develops since the heat conduction influx into the cooled zone becomes saturated with reducing *Tc*, see equation (9), and cannot compete with the local cooling there if *ν*<sup>0</sup> is large enough. The existence of a threshold for cooling instability is demonstrated in figure 4 by the time variation of *T<sup>c</sup>* computed in the TZA for different magnitudes of *ν*0: there is a characteristic change in the behaviour of the cooled zone temperature occurring by a slight

Consider now the role of physical effects neglected in equation (1): heat transfer by electron-ion coulomb collisions, particle convection caused by the pressure drop in the cooled plasma zone, the induced change in the plasma density, etc. For this purpose we use now the following system of equations, describing the evolution of the electron and ion temperature

, *φ<sup>c</sup>* =

 1 + 

1/(*β*+1)

and *T<sup>c</sup>* according to the definitions (2). With the profiles obtained above we get

Finally, for the Stirzer-Härm formula for the heat conduction, with *β* = 2.5:

*qhc* <sup>=</sup> <sup>4</sup>*Ae* 7

*<sup>ξ</sup>* <sup>≡</sup> *Lx* <sup>−</sup> *<sup>x</sup> δc*

1/*θβ*+<sup>1</sup> <sup>−</sup> <sup>1</sup>

 *δc Lx ξ*2 1/(*β*+1)

(7)

(8)

(9)

**Figure 4.** The time evolution of the temperature value averaged over the cooled zone and computed with slightly different cooling rates *ν*<sup>0</sup> in the vicinity of the cooling instability threshold.

$$1.5\partial\_l P\_\varepsilon + \partial\_\mathbf{x} \left(2.5\Gamma T\_\varepsilon - \kappa\_\varepsilon \partial\_\mathbf{x} T\_\varepsilon\right) = m\nu\_h \left(T\_0 - T\_\varepsilon\right) - Q\_{\rm ei} - \nu\_l n E\_l \tag{10}$$

$$(1.5\partial\_t P\_i + \partial\_x \left(2.5\Gamma T\_i - \kappa\_i \partial\_x T\_i\right) = m\nu\_h \left(T\_0 - T\_i\right) + Q\_{\text{ei}} \tag{11}$$

where *Pe*,*<sup>i</sup>* = *nTe*,*<sup>i</sup>* are the pressures of the plasma components, *Qei* = <sup>3</sup>*me mi n <sup>τ</sup>ei* (*Te* − *Ti*) is collisional heat exchange between electrons and ions, with *me* and *mi* being the electron and ion masses, correspondingly, and *<sup>τ</sup>ei* <sup>∼</sup> *<sup>T</sup>*1.5 *<sup>e</sup>* /*n* the time between electron-ion coulomb collisions. The plasma density *n* and the particle flux density Γ are governed by particle continuity and momentum conservation equations:

$$
\partial\_t n + \partial\_\mathbf{x} \Gamma = \nu\_\mathbf{n} \left( n\_0 - n \right) \tag{12}
$$

$$
\partial\_l \Gamma + \partial\_\mathbf{x} \left( \frac{\Gamma^2}{n} + n \frac{T\_\ell + T\_i}{m\_i} \right) = -\nu\_\Gamma \Gamma \tag{13}
$$

with the terms in the rhs of equations (12) and (13) ensuring the stability of the stationary state with *n* = *n*<sup>0</sup> and Γ = 0 before the initiation of the local cooling. The boundary conditions to equations (12) and (13) correspond to the impenetrability of the stripe borders for particles, Γ (*t*, *x* = 0) = Γ (*t*, *x* = *Lx*) = 0.

The TZA for the system of equations (10)-(13) follows from the averaging of them over the hot and cooled zones:

$$1.5\partial\_{\boldsymbol{l}}\left<\boldsymbol{P}\_{\boldsymbol{\varepsilon}}\right>\_{\boldsymbol{h}} + \frac{q\_{\boldsymbol{h}c,\boldsymbol{\varepsilon}}}{L\_{\boldsymbol{x}} - \delta\_{\boldsymbol{c}}} = \nu\_{\boldsymbol{h}}\left(\left<\boldsymbol{n}\right>\_{\boldsymbol{h}}T\_{0} - \left<\boldsymbol{P}\_{\boldsymbol{\varepsilon}}\right>\_{\boldsymbol{h}}\right) - \frac{3m\_{\boldsymbol{\varepsilon}}}{m\_{\boldsymbol{l}}} \left<\frac{\boldsymbol{n}}{\tau\_{\boldsymbol{c}\boldsymbol{l}}}\left(T\_{\boldsymbol{\varepsilon}} - T\_{\boldsymbol{l}}\right)\right>\_{\boldsymbol{h}}\tag{14}$$

**Figure 5.** The time evolution of the plasma parameters averaged over the hot (solid curves and stars) and cooled (dashed curves and boxes) zones computed numerically (curves) and with TZA equations (symbols) for *ν<sup>h</sup>* = 102*s*−<sup>1</sup> and *<sup>ν</sup>*<sup>0</sup> <sup>=</sup> 1.8 <sup>×</sup> <sup>10</sup>9*s*−1: electron temperature (a), ion temperature (b) and plasma density (c).

<sup>208</sup> Heat Transfer Studies and Applications Modelling of Global Heat Transfer by Local Cooling in Fusion Plasmas 9 10.5772/60642 Modelling of Global Heat Transfer by Local Cooling in Fusion Plasmas http://dx.doi.org/10.5772/60642 209

$$1.5\partial\_{\mathcal{l}}\left\_{\mathcal{L}} - \frac{q\_{\rm lc.}\epsilon}{\delta\_{\mathcal{L}}} = \nu\_{\rm ll}\left(\left\_{\mathcal{L}}T\_{0} - \left\_{\mathcal{L}}\right) - \frac{3m\_{\mathcal{E}}}{m\_{\mathcal{I}}} \left<\frac{n}{\tau\_{\rm el}}\left(T\_{\mathcal{E}} - T\_{\mathcal{i}}\right)\right>\_{\mathcal{L}} - \left<\nu\_{\rm l}n\right>\_{\mathcal{L}}E\_{\mathcal{l}}\tag{15}$$

8 Heat Transfer Studies and Applications

20

20

0

0.8

and *<sup>ν</sup>*<sup>0</sup> <sup>=</sup> 1.8 <sup>×</sup> <sup>10</sup>9*s*−1: electron temperature (a), ion temperature (b) and plasma density (c).

1.0

1.2

1.4

1.6

1.8

2.0

40

60

80

100 6 7 T (eV) i h,c

6 7 n (10 m ) h,c

0 10 20 t ( s) m

**Figure 5.** The time evolution of the plasma parameters averaged over the hot (solid curves and stars) and cooled (dashed curves and boxes) zones computed numerically (curves) and with TZA equations (symbols) for *ν<sup>h</sup>* = 102*s*−<sup>1</sup>

<sup>19</sup> -3 (c)

40

60

80

100 6 7 T (eV) e h,c

(a)

(b)

$$1.5\partial\_{\boldsymbol{l}}\left<\boldsymbol{P}\_{\boldsymbol{l}}\right>\_{\boldsymbol{h}} + \frac{q\_{\boldsymbol{h}c\boldsymbol{j}}}{L\_{\boldsymbol{x}} - \delta\_{\boldsymbol{c}}} = \nu\_{\boldsymbol{h}}\left(\left<\boldsymbol{n}\right>\_{\boldsymbol{h}}T\_{0} - \left<\boldsymbol{P}\_{\boldsymbol{l}}\right>\_{\boldsymbol{h}}\right) + \frac{3\boldsymbol{m}\_{\boldsymbol{c}}}{\boldsymbol{m}\_{\boldsymbol{l}}} \left<\frac{\boldsymbol{n}}{\tau\_{\boldsymbol{c}\boldsymbol{l}}}\left(T\_{\boldsymbol{c}} - T\_{\boldsymbol{l}}\right)\right>\_{\boldsymbol{h}}\tag{16}$$

$$1.5\partial\_{\mathfrak{t}}\left\_{\mathcal{L}} - \frac{q\_{\hbar c\bar{\mathcal{I}}}}{\delta\_{\mathcal{L}}} = \nu\_{\hbar}\left(\left<\eta\right>\_{\mathcal{L}}T\_{0} - \left\_{\mathcal{L}}\right) + \frac{3m\_{\mathcal{E}}}{m\_{\mathrm{i}}} \left<\frac{n}{\pi\_{\mathfrak{e}\bar{\mathcal{I}}}}\left(T\_{\mathfrak{e}} - T\_{\mathrm{i}}\right)\right>\_{\mathcal{L}}\tag{17}$$

$$
\partial\_{\hbar} \left< \mathfrak{n} \right>\_{\hbar} + \frac{\Gamma\_{\hbar c}}{L\_{\mathfrak{X}} - \delta\_{\mathfrak{c}}} = \nu\_{\hbar} \left( \mathfrak{n}\_{0} - \left< \mathfrak{n} \right>\_{\hbar} \right) \tag{18}
$$

$$
\Delta \partial\_t \left< \eta \right>\_{\mathcal{L}} - \frac{\Gamma\_{\hbar \mathcal{c}}}{\delta\_{\mathcal{c}}} = \nu\_{\hbar} \left( \eta\_0 - \left< \eta \right>\_{\mathcal{L}} \right) \tag{19}
$$

$$\left< \partial\_{\mathbf{l}} \left< \Gamma \right> + \frac{P\_{\rm ef} + P\_{\rm ic} - P\_{\rm ell} - P\_{\rm ill}}{m\_{\rm i} L\_{\rm x}} = -\nu\_{\Gamma} \left< \Gamma \right> \tag{20}$$

where <sup>Γ</sup>*hc* is the density of the particle outflow from hot to cooled zone and <sup>Γ</sup> <sup>=</sup> *Lx* <sup>0</sup> <sup>Γ</sup> (*t*, *<sup>x</sup>*) *dx Lx* . To relate the parameters at the zone edges, *fh*,*c*,*hc*, with those averaged over zones, *f<sup>h</sup>*,*c*, we apply an approximate approach similar to that outlined in the previous section. We start by relating Γ and Γ*hc* = Γ (*t*, *x* = *Lx* − *δc*). For this purpose equation (12) is represented as *∂x*Γ = *S* where the source/sink term *S* is assumed as **independent** of *x* in each of the zones. With zero boundary conditions at *x* = 0, *Lx* this results in:

$$\begin{aligned} \Gamma\left(\mathbf{0} \le \mathbf{x} \le L\_{\mathbf{x}} - \delta\_{\mathbf{c}}\right) &= \Gamma\_{\hbar \mathbf{c}} \mathbf{x} / \left(L\_{\mathbf{x}} - \delta\_{\mathbf{c}}\right) \\ \Gamma\left(L\_{\mathbf{x}} - \delta\_{\mathbf{c}} \le \mathbf{x} \le L\_{\mathbf{x}}\right) &= \Gamma\_{\hbar \mathbf{c}} \left(L\_{\mathbf{x}} - \mathbf{x}\right) / \delta\_{\mathbf{c}} \end{aligned}$$

and Γ*hc* = 2 Γ. Equation (13) is represented as *∂<sup>x</sup>* [*n* (*Te* + *Ti*)] = *F*. The term on the rhs is mostly determined by the flux density Γ and the *x*-dependence of *F* in the zones is the same as in the relations above. Straightforward calculations result in the following density profile:

$$\begin{aligned} \text{In } (0 \le \mathbf{x} \le L\_{\mathbf{x}} - \delta\_{\mathbf{c}}) &= n\_h \frac{T\_{\mathbf{c}h} + T\_{\bar{l}h}}{T\_{\mathbf{c}} + T\_{\bar{l}}} \left(1 - \varrho\_h\right) + n\_c \frac{T\_{\mathbf{c}c} + T\_{\bar{l}c}}{T\_{\mathbf{c}} + T\_{\bar{l}}} \varrho\_h \\\\ \text{In } (L\_{\mathbf{x}} - \delta\_{\mathbf{c}} \le \mathbf{x} \le L\_{\mathbf{x}}) &= n\_h \frac{T\_{\mathbf{c}h} + T\_{\bar{l}h}}{T\_{\mathbf{c}} + T\_{\bar{l}}} \varrho\_{\mathbf{c}} + n\_c \frac{T\_{\mathbf{c}c} + T\_{\bar{l}c}}{T\_{\mathbf{c}} + T\_{\bar{l}}} \left(1 - \varrho\_c\right) \end{aligned}$$

where *ϕ<sup>h</sup> <sup>ξ</sup>* <sup>=</sup> *<sup>x</sup> Lx*−*δ<sup>c</sup>* = <sup>1</sup> <sup>−</sup> *<sup>δ</sup><sup>c</sup> Lx ξ*<sup>2</sup> and *ϕ<sup>c</sup> ξ* = *Lx*−*<sup>x</sup> δc* = *<sup>δ</sup><sup>c</sup> Lx <sup>ξ</sup>*2. With relations (6) and (7) assumed for the temperature profiles of electrons and ions we get:

$$\begin{aligned} n\_{h} &= \frac{\langle n \rangle\_{h} \left(\mu\_{c}^{0} - \frac{\delta\_{c}}{L\_{x}}\mu\_{c}^{2}\right) - \langle n \rangle\_{c} \frac{T\_{c\leftarrow} + T\_{h}}{T\_{ch} + T\_{h}} \left(1 - \frac{\delta\_{c}}{L\_{x}}\right) \mu\_{h}^{2}}{\left(\mu\_{h}^{0} - \mu\_{h}^{2}\right)\mu\_{c}^{0} + \frac{\delta\_{c}}{L\_{x}} \left(\mu\_{h}^{2}\mu\_{c}^{0} - \mu\_{h}^{0}\mu\_{c}^{2}\right)} \\\\ n\_{c} &= \frac{\langle n \rangle\_{c} \left[\mu\_{h}^{0} - \left(1 - \frac{\delta\_{c}}{L\_{x}}\right)\mu\_{h}^{2}\right] - \langle n \rangle\_{h} \frac{T\_{ch} + T\_{h}}{T\_{c\leftarrow} + T\_{c}} \frac{\delta\_{c}}{L\_{x}}\mu\_{c}^{2}}{\left(\mu\_{h}^{0} - \mu\_{h}^{2}\right)\mu\_{c}^{0} + \frac{\delta\_{c}}{L\_{x}} \left(\mu\_{h}^{2}\mu\_{c}^{0} - \mu\_{h}^{0}\mu\_{c}^{2}\right)} \end{aligned}$$

where

$$\mu\_h^0 = \int\_0^1 \frac{T\_{\varepsilon h} + T\_{ih}}{T\_{\varepsilon h} \phi\_{\varepsilon h} + T\_{ih} \phi\_{i\hbar}} d\xi\_{\prime}^\varepsilon \ \mu\_c^0 = \int\_0^1 \frac{T\_{\varepsilon c} + T\_{ic}}{T\_{\varepsilon c} \phi\_{\varepsilon c} + T\_{ic} \phi\_{i\varepsilon}} d\xi\_{\prime}^\varepsilon$$

$$
\mu\_{\rm li}^2 = \int\_0^1 \frac{T\_{\rm eh} + T\_{\rm ih}}{T\_{\rm eh} \phi\_{\rm eh} + T\_{\rm ih} \phi\_{\rm ih}} \xi^2 d\xi \, \,\mu\_{\rm c}^2 = \int\_0^1 \frac{T\_{\rm cc} + T\_{\rm ic}}{T\_{\rm cc} \phi\_{\rm cc} + T\_{\rm ic} \phi\_{\rm ic}} \xi^2 d\xi \, \,\mu\_{\rm c}^2
$$

Finally, for the pressures of the plasma components we obtain:

$$
\langle P\_{\mathfrak{e},i} \rangle\_{\hbar} = (P\_{\mathfrak{e}\hbar} + P\_{\mathfrak{i}\hbar}) \left( \eta\_{\mathfrak{e},i\hbar}^{0} - \eta\_{\mathfrak{e},i\hbar}^{2} \right) + \left( P\_{\mathfrak{e}\mathfrak{e}} + P\_{\mathfrak{i}\mathfrak{e}} \right) \eta\_{\mathfrak{e},i\hbar'}^{2}
$$

$$\left\langle P\_{\mathfrak{e},\dot{\mathfrak{i}}} \right\rangle\_{\mathfrak{c}} = \left(P\_{\mathfrak{c}\mathfrak{c}} + P\_{\dot{\mathfrak{c}}\mathfrak{c}}\right) \left(\eta^{0}\_{\mathfrak{e},\dot{\mathfrak{c}}\mathfrak{c}} - \eta^{2}\_{\mathfrak{e},\dot{\mathfrak{c}}\mathfrak{c}}\right) + \left(P\_{\mathfrak{e}\mathfrak{h}} + P\_{\dot{\mathfrak{h}}\mathfrak{h}}\right) \eta^{2}\_{\mathfrak{e},\dot{\mathfrak{c}}\mathfrak{c}}$$

with

$$\eta^{0}\_{\varepsilon,\mathrm{i}\hbar} = \int\_{0}^{1} \frac{d\mathfrak{F}}{1 + \frac{T\_{i,\mathrm{c}\hbar\Phi\mathrm{i},\mathrm{c}\hbar}}{T\_{\varepsilon,\mathrm{i}\hbar}\Phi\_{\varepsilon,\mathrm{i}\hbar}}} \prime \ \eta^{0}\_{\varepsilon,\mathrm{i}\varepsilon} = \int\_{0}^{1} \frac{d\mathfrak{F}}{1 + \frac{T\_{i,\mathrm{c}\hbar}\Phi\_{\mathrm{i},\mathrm{c}\hbar}}{T\_{\varepsilon,\mathrm{i}\hbar}\Phi\_{\mathrm{c},\mathrm{i}\mathrm{c}\hbar}}} \prime$$

and

$$\eta\_{\varepsilon,ih}^2 = \left(1 - \frac{\delta\_c}{L\_{\chi}}\right) \int\_0^1 \frac{\mathfrak{F}^2 d\mathfrak{F}}{1 + \frac{T\_{i,ch}\Phi\_{i,ch}}{T\_{\varepsilon,ih}\Phi\_{\varepsilon,ih}}} \prime \eta\_{\varepsilon,ic}^2 = \frac{\delta\_c}{L\_{\chi}} \int\_0^1 \frac{\mathfrak{F}^2 d\mathfrak{F}}{1 + \frac{T\_{i,c}\Phi\_{i,c}}{T\_{\varepsilon,ic}\Phi\_{\varepsilon,ih}}}$$

The ion heat conduction has the same temperature dependence as the electron one, *κ<sup>i</sup>* = *AiT*2.5 *<sup>i</sup>* , but *Ai* is by the factor <sup>√</sup>*me*/*mi* <sup>≈</sup> 0.016 for deuterium plasmas smaller than *Ae*. Therefore the ion temperature, being reduced in the cooled zone by collisions with electrons cannot drop in the hot zone so fast as the electron one. Through *e* − *i* collisions in the hot area this may slow down the drop of *Te<sup>h</sup>*. However due to smallness of the mass ratio *me*/*mi* coulomb collisions are very ineffective for the heat transfer between the plasma components. This confirmed by the results of calculations presented in figure 5 and by comparing them with those in figure 3. Again, the agreement between the time evolutions of *Te<sup>h</sup>* and *Te<sup>c</sup>* computed by direct numerical integration of equations (10)-(13) and that

**Figure 6.** The tokamak geometry: *x* is the direction of the magnetic field whose lines of force create a closed magnetic surface of the minor radius *r*. All physical quantities are periodic in the toroidal and poloidal directions, *z* and *y*, correspondingly.

found from TZA equations (14)-(20) is nearly perfect. The results for *Ti<sup>h</sup>*, *Ti<sup>c</sup>*,*n<sup>h</sup>* and *n<sup>c</sup>* are less satisfactory, but also in this case the difference between the exact and TZA solutions is tolerable and does not exceed 20 %.

#### **3. Two-dimensional configuration**

10 Heat Transfer Studies and Applications

where

with

and

*AiT*2.5

*nh* <sup>=</sup> *<sup>n</sup><sup>h</sup>*

*nc* <sup>=</sup> *<sup>n</sup><sup>c</sup>*

0

0

*µ*0 *<sup>h</sup>* = 1

*µ*2 *<sup>h</sup>* = 1

> *Pe*,*<sup>i</sup>*

> > *Pe*,*<sup>i</sup>*

*η*2 *<sup>e</sup>*,*ih* = <sup>1</sup> <sup>−</sup> *<sup>δ</sup><sup>c</sup> Lx*

*η*0 *<sup>e</sup>*,*ih* = 1

 *µ*0 *<sup>c</sup>* <sup>−</sup> *<sup>δ</sup><sup>c</sup> Lx <sup>µ</sup>*<sup>2</sup> *c* 

 *µ*0 *<sup>h</sup>* <sup>−</sup> *<sup>µ</sup>*<sup>2</sup> *h µ*0 *<sup>c</sup>* + *<sup>δ</sup><sup>c</sup> Lx µ*2 *hµ*0 *<sup>c</sup>* <sup>−</sup> *<sup>µ</sup>*<sup>0</sup> *hµ*2 *c* 

 *µ*0 *h* − <sup>1</sup> <sup>−</sup> *<sup>δ</sup><sup>c</sup> Lx µ*2 *h* − *n<sup>h</sup>*

 *µ*0 *<sup>h</sup>* <sup>−</sup> *<sup>µ</sup>*<sup>2</sup> *h µ*0 *<sup>c</sup>* + *<sup>δ</sup><sup>c</sup> Lx µ*2 *hµ*0 *<sup>c</sup>* <sup>−</sup> *<sup>µ</sup>*<sup>0</sup> *hµ*2 *c* 

*Teh* + *Tih Tehφeh* + *Tihφih*

*Teh* + *Tih Tehφeh* + *Tihφih*

Finally, for the pressures of the plasma components we obtain:

*<sup>h</sup>* = (*Peh* + *Pih*)

*<sup>c</sup>* = (*Pec* + *Pic*)

0

 1

0

<sup>−</sup> *<sup>n</sup><sup>c</sup> Tec*+*Tic Teh*+*Tih*

*dξ*, *µ*<sup>0</sup> *<sup>c</sup>* = 1

*ξ*2*dξ*, *µ*<sup>2</sup>

 *η*0 *<sup>e</sup>*,*ih* <sup>−</sup> *<sup>η</sup>*<sup>2</sup> *e*,*ih* 

 *η*0 *<sup>e</sup>*,*ic* <sup>−</sup> *<sup>η</sup>*<sup>2</sup> *e*,*ic* 

*dξ* 1 + *Ti*,*ehφi*,*eh Te*,*ihφe*,*ih* , *η*<sup>0</sup> *<sup>e</sup>*,*ic* = 1

*ξ*2*dξ* 1 + *Ti*,*ehφi*,*eh Te*,*ihφe*,*ih*

The ion heat conduction has the same temperature dependence as the electron one, *κ<sup>i</sup>* =

*<sup>i</sup>* , but *Ai* is by the factor <sup>√</sup>*me*/*mi* <sup>≈</sup> 0.016 for deuterium plasmas smaller than *Ae*. Therefore the ion temperature, being reduced in the cooled zone by collisions with electrons cannot drop in the hot zone so fast as the electron one. Through *e* − *i* collisions in the hot area this may slow down the drop of *Te<sup>h</sup>*. However due to smallness of the mass ratio *me*/*mi* coulomb collisions are very ineffective for the heat transfer between the plasma components. This confirmed by the results of calculations presented in figure 5 and by comparing them with those in figure 3. Again, the agreement between the time evolutions of *Te<sup>h</sup>* and *Te<sup>c</sup>* computed by direct numerical integration of equations (10)-(13) and that

0

*<sup>e</sup>*,*ic* <sup>=</sup> *<sup>δ</sup><sup>c</sup> Lx* 1

, *η*<sup>2</sup>

0

0

*<sup>c</sup>* = 1  <sup>1</sup> <sup>−</sup> *<sup>δ</sup><sup>c</sup> Lx µ*2 *h*

*Teh*+*Tih Tec*+*Tic*

*Tec* + *Tic Tecφec* + *Ticφic*

*Tec* + *Tic Tecφec* + *Ticφic*

+ (*Pec* + *Pic*) *η*<sup>2</sup>

+ (*Peh* + *Pih*) *<sup>η</sup>*<sup>2</sup>

*dξ* 1 + *Ti*,*ecφi*,*ec Te*,*icφe*,*ic*

0

*δc Lx <sup>µ</sup>*<sup>2</sup> *c*

*dξ*,

*ξ*2*dξ*.

*<sup>e</sup>*,*ih*,

*e*,*ic*

,

*ξ*2*dξ* 1 + *Ti*,*ecφi*,*ec Te*,*icφe*,*ic* In this section we consider transport along closed magnetic surfaces in fusion devices to a small cooled area on the surface and generalize the TZA on this case.

#### **3.1. 2-D heat conduction equation for heat transport on a closed magnetic surface**

In magnetic fusion devices, e.g., tokamaks, electric currents, flowing both in external coils and inside the plasma, produce magnetic field components both in the toroidal and poloidal directions and lines of force generate closed nested magnetic surfaces, see figure 6. The charged plasma constituents, i.e. electrons, deuterium-tritium and impurity ions, are magnetized. This allows to confine them inside the plasma volume during times exceeding by many orders of magnitude those spent by particles, moving with their thermal velocities, to wrap magnetic surfaces very tightly. Due to this fact the temperatures of the plasma components are normally very homogeneous on the surfaces. However, if there is a very localized position where plasma is suddenly cooled down, the temperatures may become strongly inhomogeneous and it is of interest to know how promptly the whole surface will get cold in response to the local distortion. Even under the dominance of heat transfer along field lines the heat conduction to the locally cooled area cannot be straightforwardly described by a one-dimensional equation (1). This is because the majority of magnetic surfaces is non-resonant so that a field line does not close on itself and makes infinite number of turns around the surface before two points on the line come infinitesimally close each other. Thus, the coordinate *x* along the magnetic field is not periodic and no boundary conditions can be imposed for equation (1). Nonetheless, one can introduce two periodic coordinates on the torus: *z*, aligned along the toroidal direction, and *y* - along the poloidal one. The periods for this coordinates are *Lz* = 2*πR* and *Ly* = 2*πr*, respectively, with *r* and *R* being the minor and major radius of the toroidal surface. In the coordinates *z* and *y* equation (1) has the following form:

**Figure 7.** The time evolution after the start of MGI of the electron temperature averaged over the hot zone on the JET separatrix, computed with Spitzer-Härm formula for heat conduction for *ν<sup>h</sup>* = 102*s*−<sup>1</sup> the solid curve has been obtained by solving equation (21) numerically, the dashed curve and symbols - by using TZA equations (22)-(24) without and with the correction of Sc in equation (24), respectively.

$$\begin{aligned} 1.5\partial\_{\mathcal{V}}T-\cos^{2}\psi\partial\_{\mathcal{Z}}\left(\chi\partial\_{\mathcal{Z}}T\right)-\sin^{2}\psi\partial\_{\mathcal{Y}}\left(\chi\partial\_{\mathcal{Y}}T\right)-\\ -\sin\psi\cos\psi\left[\partial\_{\mathcal{Y}}\left(\chi\partial\_{\mathcal{Z}}T\right)+\partial\_{\mathcal{Z}}\left(\chi\partial\_{\mathcal{Y}}T\right)\right]=\nu\_{\mathcal{h}}\left(T\_{0}-T\right)-\nu\_{\mathcal{I}}E\_{\mathcal{I}}\end{aligned} \tag{21}$$

where *ψ* 1 is the inclination angle of the magnetic field to the toroidal direction *z*. The last cooling term is localized in the area |*Lz*/2 − *z*| , *Ly*/2 − *y* ≤ *δc*/2. By discretizing the derivatives with respect of one of coordinates, e.g., *y*, one reduces equation (21) to a set of one-dimensional PDE, similar to equation (1), which are solved by the same numerical approach. To construct a TZA for equation (21) we start again from equation (1). The only difference to the one-dimensional situation is that in the present case one does not know *a priory* the length *Lx* of the field line connecting the cooled zone with the hot rest of the surface. In reference [12] this is estimated by employing a principle of minimum entropy production [11]. By neglecting viscous and friction forces, the volume density of the entropy production rate is given as *χ* (*∂<sup>l</sup>* ln *T*) <sup>2</sup> [20]. For fixed *Th*,*<sup>c</sup>* the total entropy production rate in the parallel stripe is Θ ≈ *δcqhc* with the heat flux density from hot to the cooled zone *qhc* given by equation (9). Thus Θ is decreasing with increasing *Lx* and approaches its minimum when *Lx* reaches its maximum value corresponding to a parallel stripe covering tightly the whole magnetic surface. This corresponds to *Lx* = *LzLy*/ (2*δc*) = *Ss*/ 2 <sup>√</sup>*Sc* , where *Ss* = 4*π*2*rR* is the total area of the magnetic surface and *Sc* = *δ*<sup>2</sup> *<sup>c</sup>* the area of the cooled zone. For the time evolution of the temperature values averaged over the hot and cooled zones on the magnetic surface we get:

$$\left(1.5\partial\_{\mathfrak{h}}\left\_{\mathfrak{h}} + \frac{Q\_{\mathfrak{h}c}}{n\left(\mathcal{S}\_{\mathfrak{s}} - \mathcal{S}\_{\mathfrak{c}}\right)} = \nu\_{\mathfrak{h}}\left(T\_0 - \left\_{\mathfrak{h}}\right) \tag{22}$$

**Figure 8.** The calculated two-dimensional profile of the electron temperature at *t* = 0.3*µs*

$$1.5\partial\_{l}\left\_{c} - \frac{Q\_{\rm hc}}{nS\_{\rm c}} = \nu\_{\rm h}\left(T\_{0} - \left\_{c}\right) - \left<\nu\_{l}\right>\_{c}E\_{l} \tag{23}$$

with the total heat outflow to the cooled zone:

12 Heat Transfer Studies and Applications

6 7 T (eV) e h

0

<sup>−</sup> sin *<sup>ψ</sup>* cos *<sup>ψ</sup>*

last cooling term is localized in the area |*Lz*/2 − *z*| ,

magnetic surface. This corresponds to *Lx* = *LzLy*/ (2*δc*) = *Ss*/

1.5*∂<sup>t</sup> T<sup>h</sup>* <sup>+</sup> *Qhc*

the total area of the magnetic surface and *Sc* = *δ*<sup>2</sup>

with the correction of Sc in equation (24), respectively.

production rate is given as *χ* (*∂<sup>l</sup>* ln *T*)

surface we get:

50

100

150

200

250

20 40 t ( s) m

 *χ∂yT* 

 *Ly*/2 − *y*

<sup>=</sup> *<sup>ν</sup><sup>h</sup>* (*T*<sup>0</sup> <sup>−</sup> *<sup>T</sup>*) <sup>−</sup> *<sup>ν</sup>lEl*

<sup>2</sup> [20]. For fixed *Th*,*<sup>c</sup>* the total entropy production rate in

 2 <sup>√</sup>*Sc* 

*<sup>n</sup>* (*Ss* <sup>−</sup> *Sc*) <sup>=</sup> *<sup>ν</sup><sup>h</sup>* (*T*<sup>0</sup> <sup>−</sup> *Th*) (22)

*<sup>c</sup>* the area of the cooled zone. For the time

− (21)

≤ *δc*/2. By discretizing the

, where *Ss* = 4*π*2*rR* is

0 60 80 100

**Figure 7.** The time evolution after the start of MGI of the electron temperature averaged over the hot zone on the JET separatrix, computed with Spitzer-Härm formula for heat conduction for *ν<sup>h</sup>* = 102*s*−<sup>1</sup> the solid curve has been obtained by solving equation (21) numerically, the dashed curve and symbols - by using TZA equations (22)-(24) without and

> *χ∂yT*

where *ψ* 1 is the inclination angle of the magnetic field to the toroidal direction *z*. The

derivatives with respect of one of coordinates, e.g., *y*, one reduces equation (21) to a set of one-dimensional PDE, similar to equation (1), which are solved by the same numerical approach. To construct a TZA for equation (21) we start again from equation (1). The only difference to the one-dimensional situation is that in the present case one does not know *a priory* the length *Lx* of the field line connecting the cooled zone with the hot rest of the surface. In reference [12] this is estimated by employing a principle of minimum entropy production [11]. By neglecting viscous and friction forces, the volume density of the entropy

the parallel stripe is Θ ≈ *δcqhc* with the heat flux density from hot to the cooled zone *qhc* given by equation (9). Thus Θ is decreasing with increasing *Lx* and approaches its minimum when *Lx* reaches its maximum value corresponding to a parallel stripe covering tightly the whole

evolution of the temperature values averaged over the hot and cooled zones on the magnetic

1.5*∂tT* <sup>−</sup> cos<sup>2</sup> *ψ∂<sup>z</sup>* (*χ∂zT*) <sup>−</sup> sin<sup>2</sup> *ψ∂<sup>y</sup>*

*∂<sup>y</sup>* (*χ∂zT*) + *∂<sup>z</sup>*

$$Q\_{h\varepsilon} = \frac{16A\_{\varepsilon}}{7} \frac{S\_{\varepsilon}}{S\_{\varepsilon}} \left( T\_h^{3.5} - T\_c^{3.5} \right) \tag{24}$$

As an example of applications to 2-D configurations we consider experiments with massive gas injection into the tokamak JET to mitigate plasma disruptions [21, 22]. Calculations have been performed for the conditions on the magnetic separatrix between the confined plasma region and the SOL, with *R* = 3*m* and *r* = 1*m*, the initial values of the electron temperature and density as in the shot 77808 with MGI into the H-mode plasma [22]: *n*<sup>0</sup> = 1019*m*−<sup>3</sup> and *T*<sup>0</sup> = 250*eV*; the cooled area is a square with a side of 1*m*, i.e. *Sc* ≈ 0.01*Ss*, where the density of argon atoms is increased instantaneously up to 2 <sup>×</sup> 1024*m*−3. Curves in figure 7 represent the computed time evolution of the electron temperature values averaged over the hot zone found by solving 2-D equation (21) and 0-D TZA-equations (22)-(24). There is a noticeable difference between the results of both approaches, much larger than in the one-dimensional case discussed in the previous section. The cause of this discrepancy can be understood by looking on the calculated two-dimensional profile of the electron temperature displayed in figure 8. One can see that there is a very sharp gradient of the temperature at the border of the cooled zone, corresponding to *LT* ≡ |*T*/*∂xT*| *δc*. For the heat conduction *<sup>κ</sup>* <sup>∼</sup> *<sup>T</sup><sup>β</sup>* the TZA approach predicts *LT* <sup>≥</sup> *<sup>δ</sup><sup>c</sup>* (*<sup>β</sup>* <sup>+</sup> <sup>1</sup>) /2 at this position, i.e. *LT* <sup>≥</sup> <sup>7</sup>*δc*/4 for the Spitzer-Härm approximation used above. This means that the area of the cooled zone is actually reduced in the TZA compared with that following from 2-D calculations. One can improve TZA predictions by increasing *Sc* in equation (24). Indeed, the TZA results obtained for the cooled zone enlarged by a factor of 4 are presented in figure 7 with stars and agree with 2-D calculations very good. As a rule of thumb one can estimate the enhancement factor for *Sc* in equation (24) as 1 + (*β* + 1) <sup>2</sup> /4.

#### **3.2. Heat flux limit**

Calculations in the previous section predict that after the initiation of the localized cooling the whole magnetic surface is cooled down during a time shorter than 10*µs*. This is in a striking disagreement with the experimental data showing that such a global cooling requires at least half a millisecond. A possible cause for this discrepancy may be rooted in the fact that the Spitzer-Härm formula *κSH* has been used for the electron heat conduction. This approximation is deduced in the limit of a very small mean free path length between coulomb collisions for thermal particles, *λth*, compared to the characteristic dimension for the temperature change *LT*; it results in the heat flux density *qe* = *qSH* ≡ −*κSH∂xT*. Since *<sup>κ</sup>SH* <sup>∼</sup> *<sup>T</sup>*5/2 the heat conduction increases very fast with rising temperature. However, *<sup>λ</sup>th* <sup>∼</sup> *<sup>T</sup>*<sup>2</sup> for coulomb collisions of charged particles. Thus, the SH-approximation can be violated in the case of strong local cooling where the temperature in the cooled zone drops much faster than that in the hot one and a region with strong temperature gradient arises. Actually, the SH-approximation is violated already if the free path of particles carrying heat, *λhcp*, exceeds *LT*. The energy of such particles is by an order of magnitude larger than *T* and *λhcp* ≈ 100*λth*. In a collision-less limit, *λhcp LT*, heat is transported because more hotter particles escape from the region in question compared to colder particles entering it. This would lead to a free-streaming heat flux with the density *qe* <sup>=</sup> *qFS nT*3/2/ <sup>√</sup>*<sup>m</sup>* [23, 24]. However, interpretation of laser fusion experiments [13] and experiments with heating of magnetic islands in TEXTOR [25, 26] led to the conclusion that for electrons this heat flux level is strongly reduced by the so called heat flux limit (HFL) factor, 0.02 ≤ *ξHFL* ≤ 0.1. Two physical mechanisms can explain such a strong reduction of *qe*. First, non-local effects in collision-less plasma reduce the perturbation in the distribution function caused by the temperature gradient and resulting in heat conduction [27]. Second, light electrons are braked by the ambipolar electric field, arising because heavy ions lag behind electrons and a charge separation is induced [24]. A smooth transition between collisional and collision-less limits can be described by the formula [13], *qe* = *ξHFLqFSqSH*/ (*ξHFLqFS* + *qSH*), and resulting in

$$\kappa\_{\varepsilon} \approx \kappa\_{HFL} = \frac{\kappa\_{SH}}{1 + \lambda\_{\text{lcp}}/L\_T} \tag{25}$$

The dependence of *κ<sup>e</sup>* in equation (25) on the temperature gradient makes the procedure of numerical integration of equation (1) very unstable. This is extremely exaggerated by the fact that *qe* becomes actually independent of the temperature gradient for *λhcp LT*, i.e., the heat flow is convective and the heat conduction equation (1) becomes actually an equation of the first order. Therefore, often the relation (25) is approximated by a simple decrease of *κSH*. Here for rough estimates we assume *LT* ≈ *Lx*. For the parameters on the JET separatrix we get *Lx* <sup>≈</sup> <sup>55</sup>*<sup>m</sup> <sup>λ</sup>hcp* <sup>≈</sup> <sup>6250</sup>*m*, and equation (25) provides *<sup>χ</sup>* <sup>=</sup> *<sup>κ</sup>*/*<sup>n</sup>* <sup>≈</sup> 109*m*−1*s*1. Thus, the characteristic time for cooling of the magnetic surface *tcool* <sup>≈</sup> *<sup>L</sup>*<sup>2</sup> *<sup>x</sup>*/*χ* ≈ 0.3*ms*, in good agreement with the experiment [21, 22], especially by taking into account that the electron temperature and *χ* are reducing as the cooling progresses. Another argument for the relevance of the equation (25) is the fact that for *<sup>λ</sup>hcp LT* it provides *<sup>χ</sup>* <sup>∼</sup> <sup>√</sup>*<sup>T</sup>* and, thus, *tcool* depends weakly on the experimental conditions. This agrees with the experimental observation [21, 22] showing that *tcool* is practically the same for the magnetic surfaces with the initial temperatures of 250*eV* and 400*eV*. For *κ<sup>e</sup>* ≈ *κSH* this would result in *tcool* differing by a factor of 3.

**Figure 9.** The time evolution after the start of MGI of the electron temperature averaged over the hot zone on the JET separatrix, computed with the heat flux limit for *ν<sup>h</sup>* = 102*s*−1; the solid curve has been obtained by solving equation (21) numerically, the dashed curve and symbols - by using TZA equations (22)-(24) without and with the correction of Sc in equation (24), respectively.

To generalize the TZA by including the HFL, we estimate *λhcp* at *T* = (*Th* + *Tc*) /2 and assume 1/*LT* <sup>=</sup> (*Th* <sup>−</sup> *Tc*) /*Lx*/*<sup>T</sup>*. This results in *<sup>λ</sup>hcp*/*LT* <sup>=</sup> 0.05 *T*2 *<sup>h</sup>* <sup>−</sup> *<sup>T</sup>*<sup>2</sup> *c* / (*nLx*), where the temperatures are measured in electron-volts, plasma density - in 1019*m*−<sup>3</sup> and *Lx* in *<sup>m</sup>*. Since now *<sup>κ</sup>* <sup>∼</sup> <sup>√</sup>*T*, one would expect that a significantly smaller, compared to the temperature dependence of *<sup>κ</sup>SH* <sup>∼</sup> *<sup>T</sup>*5/2, increase of the cooled zone area in the TZA approach is necessary to reproduce well the 2-D calculations. Figure 9 shows the time evolution of the *T<sup>h</sup>* calculated for the conditions at the JET separatrix, by taking into account the approximate assessment above for the heat flux limit. One can see that both 2-D and corrected TZA approaches reproduce well the found experimentally [21, 22] characteristic time of 0.5*ms* for the temperature drop far from the cooled area.

#### **4. Three-dimensional configuration**

#### **4.1. Basic equations**

14 Heat Transfer Studies and Applications

Calculations in the previous section predict that after the initiation of the localized cooling the whole magnetic surface is cooled down during a time shorter than 10*µs*. This is in a striking disagreement with the experimental data showing that such a global cooling requires at least half a millisecond. A possible cause for this discrepancy may be rooted in the fact that the Spitzer-Härm formula *κSH* has been used for the electron heat conduction. This approximation is deduced in the limit of a very small mean free path length between coulomb collisions for thermal particles, *λth*, compared to the characteristic dimension for the temperature change *LT*; it results in the heat flux density *qe* = *qSH* ≡ −*κSH∂xT*. Since *<sup>κ</sup>SH* <sup>∼</sup> *<sup>T</sup>*5/2 the heat conduction increases very fast with rising temperature. However, *<sup>λ</sup>th* <sup>∼</sup> *<sup>T</sup>*<sup>2</sup> for coulomb collisions of charged particles. Thus, the SH-approximation can be violated in the case of strong local cooling where the temperature in the cooled zone drops much faster than that in the hot one and a region with strong temperature gradient arises. Actually, the SH-approximation is violated already if the free path of particles carrying heat, *λhcp*, exceeds *LT*. The energy of such particles is by an order of magnitude larger than *T* and *λhcp* ≈ 100*λth*. In a collision-less limit, *λhcp LT*, heat is transported because more hotter particles escape from the region in question compared to colder particles entering it. This

would lead to a free-streaming heat flux with the density *qe* <sup>=</sup> *qFS nT*3/2/

However, interpretation of laser fusion experiments [13] and experiments with heating of magnetic islands in TEXTOR [25, 26] led to the conclusion that for electrons this heat flux level is strongly reduced by the so called heat flux limit (HFL) factor, 0.02 ≤ *ξHFL* ≤ 0.1. Two physical mechanisms can explain such a strong reduction of *qe*. First, non-local effects in collision-less plasma reduce the perturbation in the distribution function caused by the temperature gradient and resulting in heat conduction [27]. Second, light electrons are braked by the ambipolar electric field, arising because heavy ions lag behind electrons and a charge separation is induced [24]. A smooth transition between collisional and collision-less limits can be described by the formula [13], *qe* = *ξHFLqFSqSH*/ (*ξHFLqFS* + *qSH*),

*<sup>κ</sup><sup>e</sup>* <sup>≈</sup> *<sup>κ</sup>HFL* <sup>=</sup> *<sup>κ</sup>SH*

Thus, the characteristic time for cooling of the magnetic surface *tcool* <sup>≈</sup> *<sup>L</sup>*<sup>2</sup>

The dependence of *κ<sup>e</sup>* in equation (25) on the temperature gradient makes the procedure of numerical integration of equation (1) very unstable. This is extremely exaggerated by the fact that *qe* becomes actually independent of the temperature gradient for *λhcp LT*, i.e., the heat flow is convective and the heat conduction equation (1) becomes actually an equation of the first order. Therefore, often the relation (25) is approximated by a simple decrease of *κSH*. Here for rough estimates we assume *LT* ≈ *Lx*. For the parameters on the JET separatrix we get *Lx* <sup>≈</sup> <sup>55</sup>*<sup>m</sup> <sup>λ</sup>hcp* <sup>≈</sup> <sup>6250</sup>*m*, and equation (25) provides *<sup>χ</sup>* <sup>=</sup> *<sup>κ</sup>*/*<sup>n</sup>* <sup>≈</sup> 109*m*−1*s*1.

good agreement with the experiment [21, 22], especially by taking into account that the electron temperature and *χ* are reducing as the cooling progresses. Another argument for the relevance of the equation (25) is the fact that for *<sup>λ</sup>hcp LT* it provides *<sup>χ</sup>* <sup>∼</sup> <sup>√</sup>*<sup>T</sup>* and, thus, *tcool* depends weakly on the experimental conditions. This agrees with the experimental observation [21, 22] showing that *tcool* is practically the same for the magnetic surfaces with the initial temperatures of 250*eV* and 400*eV*. For *κ<sup>e</sup>* ≈ *κSH* this would result in *tcool* differing

1 + *λhcp*/*LT*

<sup>√</sup>*<sup>m</sup>* [23, 24].

(25)

*<sup>x</sup>*/*χ* ≈ 0.3*ms*, in

**3.2. Heat flux limit**

and resulting in

by a factor of 3.

The TZA elaborated in previous sections for 1-D and 2-D situations provides a basis for modelling with reasonable accuracy and CPU time of the global response on a local cooling in a three-dimensional case, e.g., by injection of impurities into a tokamak plasma. The 3-D heat transport equation for electrons has the form:

$$\begin{aligned} \left[\partial\_t \left(1.5\boldsymbol{n}\_{\boldsymbol{\varepsilon}}\boldsymbol{T}\_{\boldsymbol{\varepsilon}}\right) + \partial\_r \left(r\boldsymbol{q}\_{\boldsymbol{r}}\right)/r - \cos^2\psi \partial\_z \left(\kappa\_{||\boldsymbol{\beta}\_{\boldsymbol{z}}\boldsymbol{T}\_{\boldsymbol{\varepsilon}}}\right) - \sin^2\psi \partial\_{\boldsymbol{\beta}}\left(\kappa\_{||\boldsymbol{\beta}\_{\boldsymbol{\beta}}\boldsymbol{T}\_{\boldsymbol{\varepsilon}}}\right) \\ - \sin\psi \cos\psi \left[\partial\_{\boldsymbol{\beta}}\left(\kappa\_{||\boldsymbol{\beta}\_{\boldsymbol{z}}\boldsymbol{T}\_{\boldsymbol{\varepsilon}}}\right) + \partial\_z \left(\kappa\_{||\boldsymbol{\beta}\_{\boldsymbol{z}}\boldsymbol{T}\_{\boldsymbol{\varepsilon}}}\right)\right] = Q\_h - Q\_{\rm loss} \end{aligned} \tag{26}$$

Here we take into account that in fusion devices the electron density *ne* varies with the third coordinate, the minor radius *r* of magnetic surfaces, and the parallel heat conduction of electrons, *<sup>κ</sup>*||, is used in this equation instead of the heat conductivity *<sup>χ</sup>* <sup>=</sup> *<sup>κ</sup>*/*<sup>n</sup>* used in equation (1). The heat flux density across the magnetic surfaces includes both conduction and convection contributions:

$$q\_r = -\kappa\_r \partial\_r T\_\varepsilon + 1.5 \Gamma\_r T\_\varepsilon$$

where *κ<sup>r</sup>* and Γ*<sup>r</sup>* are the heat conduction and density of the electron flux in the radial direction; *Qh* is the density of heating power due to different mechanisms, such as Ohmic dissipation of the plasma current, energetic neutral beams, radio-frequency electromagnetic waves, etc.; henceforth the *r*-dependences of *κr*, Γ*<sup>r</sup>* and *Qh* are prescribed and assumed constant in time. The localized cooling is described by the term *Qloss*. Before its initiation at *t* = 0 the quasi-stationary state is maintained by the balance of heating and radial heat outflow toward the plasma edge. In the SOL region, *rs* ≤ *r* ≤ *rw*, where *rs* is the radius of the last closed magnetic surface (LCMS) and *rw* that of the device wall, the heat is additionally lost along the magnetic field lines to limiter or divertor target plates. The boundary conditions to equation (26) are imposed at the plasma axis, *r* = 0, where *∂rT* = 0, and at the wall, *r* = *rw*, where the so called *e*-folding length *δ<sup>T</sup>* is fixed and *∂rT* = −*T*/*δT*.

By deducing the TZA-equations for *Th*,*<sup>c</sup>* (*t*,*r*) we have to take into account a possible variation of the cooled and hot zone areas, *Sc* and *Sh* = *Ss* − *Sc*, respectively, with *t* and *r*. The averaging of equation (26) over the zones results in:

$$\begin{aligned} \partial\_{\hbar} \left( 1.5 \mathcal{S}\_{\hbar} n\_{\varepsilon} \left< T\_{\varepsilon} \right>\_{\hbar} \right) &+ \mathcal{S}\_{\hbar} \partial\_{r} \left[ r \left( -\kappa\_{r} \partial\_{r} \left< T\_{\varepsilon} \right>\_{\hbar} + 1.5 \Gamma\_{r} \left< T\_{\varepsilon} \right>\_{\hbar} \right) \right] / r = \\ &= \mathcal{S}\_{\hbar} \left( Q\_{\hbar} - Q\_{\rm SOL} \right) - Q\_{\hbar \varepsilon} \end{aligned} \tag{27}$$

$$\begin{aligned} \partial\_t \left( 1.5 S\_{\mathbb{C}} \left< n\_{\varepsilon} T\_{\ell} \right>\_{\mathcal{L}} \right) + S\_{\mathbb{C}} \partial\_r \left[ r \left( -\kappa\_r \partial\_r \left< T\_{\ell} \right>\_{\mathcal{L}} + 1.5 \Gamma\_r \left< T\_{\ell} \right>\_{\mathcal{L}} \right) \right] / r = \\ = S\_{\mathbb{C}} \left( Q\_h - Q\_{\text{SOL}} \right) Q\_{hc} - W\_{\text{loss}} \end{aligned} \tag{28}$$

where *QSOL* is non-zero for *rs* ≤ *r* and gives the heat loss to the target plates, see reference [28]; the heat outflow to the cooled zone is computed by taking into account the heat flux limit and the correction due to smoothing in the TZA of the temperature gradient at the border between zones:

$$Q\_{hc} = 4A\_c \frac{T\_h^{3.5} - T\_c^{3.5}}{\frac{S\_s}{S\_c} + 0.1 \frac{T\_h^2 - T\_c^2}{n \delta\_l}} \tag{29}$$

the plasma density in the hot zone is assumed unchanged, but if the local cooling is caused by the injection of neutral particles, *ne* may be modified in the cooled one due to the ionization of these particles; *Wloss* is the integral of the loss term *Qloss* over the cooled zone.

**Figure 10.** Spreading in the plasma of impurity (carbon) particles released from the wall by LIAS.

16 Heat Transfer Studies and Applications

and convection contributions:

border between zones:

*r*. The averaging of equation (26) over the zones results in:

= *Sh* (*Qh* − *QSOL*) − *Qhc*

= *Sc* (*Qh* − *QSOL*) *Qhc* − *Wloss*

Here we take into account that in fusion devices the electron density *ne* varies with the third coordinate, the minor radius *r* of magnetic surfaces, and the parallel heat conduction of electrons, *<sup>κ</sup>*||, is used in this equation instead of the heat conductivity *<sup>χ</sup>* <sup>=</sup> *<sup>κ</sup>*/*<sup>n</sup>* used in equation (1). The heat flux density across the magnetic surfaces includes both conduction

*qr* = −*κr∂rTe* + 1.5Γ*rTe* where *κ<sup>r</sup>* and Γ*<sup>r</sup>* are the heat conduction and density of the electron flux in the radial direction; *Qh* is the density of heating power due to different mechanisms, such as Ohmic dissipation of the plasma current, energetic neutral beams, radio-frequency electromagnetic waves, etc.; henceforth the *r*-dependences of *κr*, Γ*<sup>r</sup>* and *Qh* are prescribed and assumed constant in time. The localized cooling is described by the term *Qloss*. Before its initiation at *t* = 0 the quasi-stationary state is maintained by the balance of heating and radial heat outflow toward the plasma edge. In the SOL region, *rs* ≤ *r* ≤ *rw*, where *rs* is the radius of the last closed magnetic surface (LCMS) and *rw* that of the device wall, the heat is additionally lost along the magnetic field lines to limiter or divertor target plates. The boundary conditions to equation (26) are imposed at the plasma axis, *r* = 0, where *∂rT* = 0, and at the wall, *r* = *rw*, where the so called *e*-folding length *δ<sup>T</sup>* is fixed and *∂rT* = −*T*/*δT*. By deducing the TZA-equations for *Th*,*<sup>c</sup>* (*t*,*r*) we have to take into account a possible variation of the cooled and hot zone areas, *Sc* and *Sh* = *Ss* − *Sc*, respectively, with *t* and

*∂<sup>t</sup>* (1.5*Shne Teh*) + *Sh∂<sup>r</sup>* [*r* (−*κr∂<sup>r</sup> Te<sup>h</sup>* + 1.5Γ*<sup>r</sup> Teh*)] /*r* = (27)

*∂<sup>t</sup>* (1.5*Sc neTec*) + *Sc∂<sup>r</sup>* [*r* (−*κr∂<sup>r</sup> Te<sup>c</sup>* + 1.5Γ*<sup>r</sup> Tec*)] /*r* = (28)

*<sup>h</sup>* −*T*<sup>2</sup> *c nδ<sup>c</sup>*

(29)

where *QSOL* is non-zero for *rs* ≤ *r* and gives the heat loss to the target plates, see reference [28]; the heat outflow to the cooled zone is computed by taking into account the heat flux limit and the correction due to smoothing in the TZA of the temperature gradient at the

> *T*3.5 *<sup>h</sup>* <sup>−</sup> *<sup>T</sup>*3.5 *c*

*Ss Sc* <sup>+</sup> 0.1 *<sup>T</sup>*<sup>2</sup>

the plasma density in the hot zone is assumed unchanged, but if the local cooling is caused by the injection of neutral particles, *ne* may be modified in the cooled one due to the ionization

*Qhc* = 4*Ae*

of these particles; *Wloss* is the integral of the loss term *Qloss* over the cooled zone.

#### **4.2. Example of application: Plasma distortion by laser-induced ablation spectroscopy**

Laser-based diagnostics have been proposed to measure and monitor *in situ* the composition of layers on plasma facing components (PFC) [4]. By laser-induced ablation spectroscopy (LIAS) a short, in nanosecond-range, intense laser pulse is directed at the PFC during a plasma discharge. Particles of the PFC material, released into the plasma, are excited and ionized by electrons. By measuring the intensity of impurity line radiation, one can assess the amount of particles emitted and judge about the wall composition. To interpret LIAS measurements local plasma parameters, especially the electron density and temperature, have to be known. Normally the parameters, measured somewhere at the plasma edge before the LIAS application, are used. However, processes with particles released from the wall may lead to significant disturbances in the plasma: the energy loss on the particle excitation leads to reduction of the electron temperature, the delivery of new electrons by the ionization to the increase of their density. An assessment of such perturbations is of importance for quantitative interpretation of measurements [29].

Consider the spreading of carbon particles released by the laser pulse from the wall, see figure 10. The penetration depth of these species is usually significantly smaller than the radius *rw* of the wall element in question. Therefore, it is convenient to use an orthogonal coordinate system with the axes *r*, directed from the plasma to the injection spot, *x* and *y*, being tangential to the magnetic surface and oriented parallel and perpendicular to field lines, respectively, see figure 10. The velocity distribution of neutral particles is identical in the directions *x* and *y* and the distribution function can be characterized by its dependence on the velocity components *Vr* and *V<sup>ρ</sup>* in a cylindrical reference system with axes *r* and *<sup>ρ</sup>* = *<sup>x</sup>*<sup>2</sup> + *<sup>y</sup>*2. At the wall the fraction of such particles in the total density *<sup>n</sup>*<sup>0</sup> of neutrals is *f Vr*, *V<sup>ρ</sup> dVrdVρ*, where the distribution function *f Vr*, *V<sup>ρ</sup>* is assumed as a Maxwellian one at the temperature *T*<sup>0</sup> with respect to *V<sup>ρ</sup>* and a one-side Maxwellian, shifted with the drift velocity *Vm*, for *Vr* [30]:

$$f\left(V\_{r\prime}, V\_{\rho}\right) = \frac{4V\_{\rho}}{\sqrt{\pi}\left[1 + \operatorname{erf}\left(V\_{m}/V\_{th}\right)\right]V\_{th}^{3}}\exp\left[-\frac{\left(V\_{r} + V\_{m}\right)^{2} + V\_{\rho}^{2}}{V\_{th}^{2}}\right] \tag{30}$$

Here *Vth* <sup>=</sup> <sup>√</sup>2*T*0/*m*, with *<sup>m</sup>* being the particle mass, and the form-factor arises due to the normalization <sup>0</sup> <sup>−</sup><sup>∞</sup> *dVr* <sup>∞</sup> <sup>0</sup> *f Vr*, *V<sup>ρ</sup> dV<sup>ρ</sup>* = 1. The density *dn*<sup>0</sup> of neutrals with the velocity components in the ranges *dVr* and *dV<sup>ρ</sup>* is determined from the continuity equation:

**Figure 11.** The radial profiles of the electron density (solid curve), temperature (dashed curve) and heating power density (dash-dotted curve) used by the modelling of the spreading process of carbon particles released by LIAS.

$$
\partial\_l dn\_0 + V\_r \partial\_r dn\_0 + (V\_\rho / \rho) \times \partial\_\rho (\rho d n\_0) = -k\_{ion}^0 n\_\epsilon dn\_0 \tag{31}
$$

where *k*<sup>0</sup> *ion* is the ionization rate coefficient. An approximate solution of equation (31) is searched for in the form:

$$dn\_0(t, r, \rho) = d\eta\_0(t, r) \times \exp\left[-\rho^2 / \lambda^2(t, r)\right] \tag{32}$$

This form of the solution mimics the fact that on each magnetic surface the density of injected neutral particles is localized in some vicinity of the ejection position and decays to zero far from this due to ionization.

Introduce new dependent variables:

$$dN\_0(t,r) = \int\_0^\infty dn\_0 2\pi\rho d\rho = \pi\lambda^2 d\eta\_0\tag{33}$$

$$d\Lambda\_0(t,r) = \int\_0^\infty dn\_0 2\pi\rho^2 d\rho = \frac{\sqrt{\pi}}{2}\lambda dN\_0$$

Equations for *dN*0(*t*,*r*) and *d*Λ0(*t*,*r*) follow from the integration of equation (31) with respect to *ρ* with the weights 2*πρ* and 2*πρ*2, respectively, and with *dn*<sup>0</sup> substituted in the form (32):

$$
\partial\_t dN\_0 + V\_r \partial\_r dN\_0 = -k\_{ion}^0 n\_\epsilon dN\_0 \tag{34}
$$

$$
\partial\_t d\Lambda\_0 + V\_r \partial\_r d\Lambda\_0 = V\_\rho dN\_0 - k\_{ion}^0 n\_\epsilon d\Lambda\_0
$$

18 Heat Transfer Studies and Applications

where *k*<sup>0</sup>

searched for in the form:

from this due to ionization.

Introduce new dependent variables:

0

0 0.1 0.2 0.3 0.4 0.5

**Figure 11.** The radial profiles of the electron density (solid curve), temperature (dashed curve) and heating power density (dash-dotted curve) used by the modelling of the spreading process of carbon particles released by LIAS.

This form of the solution mimics the fact that on each magnetic surface the density of injected neutral particles is localized in some vicinity of the ejection position and decays to zero far

*dn*02*πρ*2*dρ* =

Equations for *dN*0(*t*,*r*) and *d*Λ0(*t*,*r*) follow from the integration of equation (31) with respect to *ρ* with the weights 2*πρ* and 2*πρ*2, respectively, and with *dn*<sup>0</sup> substituted in the form (32):

*ion* is the ionization rate coefficient. An approximate solution of equation (31) is

<sup>−</sup>*ρ*2/*λ*<sup>2</sup> (*t*,*r*)

<sup>√</sup>*<sup>π</sup>* <sup>2</sup> *<sup>λ</sup>dN*<sup>0</sup>

*ionned*Λ<sup>0</sup>

*dn*02*πρdρ* = *πλ*2*dη*<sup>0</sup> (33)

*ionnedN*<sup>0</sup> (34)

*<sup>∂</sup>tdn*<sup>0</sup> <sup>+</sup> *Vr∂rdn*<sup>0</sup> + (*Vρ*/*ρ*) <sup>×</sup> *∂ρ*(*ρdn*0) = <sup>−</sup>*k*<sup>0</sup>

*dn*0(*t*,*r*, *ρ*) = *dη*0(*t*,*r*) × exp

∞

0

∞

0

*<sup>∂</sup>tdN*<sup>0</sup> <sup>+</sup> *Vr∂rdN*<sup>0</sup> <sup>=</sup> <sup>−</sup>*k*<sup>0</sup>

*<sup>∂</sup>td*Λ<sup>0</sup> <sup>+</sup> *Vr∂rd*Λ<sup>0</sup> <sup>=</sup> *<sup>V</sup>ρdN*<sup>0</sup> <sup>−</sup> *<sup>k</sup>*<sup>0</sup>

*dN*0(*t*,*r*) =

*d*Λ0(*t*,*r*) =

r(m)

*ionnedn*<sup>0</sup> (31)

(32)

n(4x10 m ), T (keV), Q 19 -3 <sup>2</sup> -3

e h(10 kW m

0.2

0.4

0.6

0.8

1

**Figure 12.** Time evolution of the cumulative energy radiated by *C*<sup>+</sup> ions calculated with (solid curve) and without (dashed curve) taking into account the plasma distortion by LIAS.

The boundary conditions at the wall take into account that no particles are released after the irradiation pulse, i.e. *dN*<sup>0</sup> (*t* > 0,*rw*) = *d*Λ<sup>0</sup> (*t* > 0,*rw*) = 0, since the release occurs during a time by orders of magnitude shorter than a typical time step *τ* in calculations. The latter has to be, however, chosen so that the distance *s* = *τ* |*Vr*| covered by particles in the radial direction during one time step is significantly smaller than their free path before collisions with electrons. In such a case the actual initial radial profile of *dn*<sup>0</sup> is of no importance. Here we assume that at *t* = 0 all particles fill homogeneously the first sell, *rw* − *h* ≤ *r* ≤ *rw*, but practically the same results are obtained for other assumptions, e.g. for the initial density profile decaying exponentially with *s* assumed as the *e*-folding length; along the wall the source region is localized in the irradiated spot with the radius *ρ*0. For the variables *dN*<sup>0</sup> and *d*Λ<sup>0</sup> this results in the following initial conditions:

$$dN\_0 \left( 0, r\_{\overline{w}} - h \le r \right) = \frac{N\_{\text{tot}}}{h} f \left( V\_{r\prime} V\_{\rho} \right) dV\_r dV\_{\rho}; \ d\Lambda\_0 \left( 0, r \right) = \frac{\sqrt{\pi}}{2} \rho\_0 d N\_0 \left( 0, r \right)$$

where *Ntot* is the total amount of injected particles. With known *dN*<sup>0</sup> (*t*,*r*) and *d*Λ<sup>0</sup> (*t*,*r*) one can obtain the original parameters characterizing the local density of neutrals:

$$d\eta\_0\left(t,r\right) = \frac{\left(dN\_0\right)^3}{\left(d\Lambda\_0\right)^2};\ \lambda\left(t,r\right) = \frac{2}{\sqrt{\pi}}\frac{d\Lambda\_0}{dN\_0} \tag{35}$$

The *ρ*-profile of the total density of neutral particles is approximated analogously to the relation (32)

$$m\_0(t, r, \rho) = \eta\_0(t, r) \times \exp\left[-\rho^2 / \lambda\_0^2(t, r)\right] \tag{36}$$

**Figure 13.** Calculated time evolution of the electron temperature (solid curve) and density (dashed curve) induced by LIAS in the cooled zone at the last closed magnetic surface.

where the maximum density *<sup>η</sup>*<sup>0</sup> (*t*,*r*) = *<sup>d</sup>η*<sup>0</sup> *Vr*, *V<sup>ρ</sup>* and the radius of the localization region is defined by conserving the total number of particles on the magnetic surface, *<sup>N</sup>*<sup>0</sup> (*t*,*r*) = *dN*<sup>0</sup> *Vr*, *V<sup>ρ</sup>* :

$$
\lambda\_0 \left( t, r \right) = \sqrt{\frac{N\_0 \left( t, r \right)}{\pi \eta\_0 \left( t, r \right)}} \tag{37}
$$

The area of the cooled zone is assumed henceforth to be equal to:

$$S\_{\mathcal{C}}\left(t,r\right) = \pi\lambda\_0^2\left(t,r\right) \tag{38}$$

The characteristic time for LIAS is of 10*µs*. According to figures 5b,c it is too short for a noticeable reaction of the main ions and their density *ni* is considered as unperturbed. The electron density *ne* may be, however, significantly distorted by the generation of new electrons through the ionization of impurity neutrals; due to the plasma quasi-neutrality one has *ne* = *ni* + *n*1. Here *n*<sup>1</sup> is the density of singly charged impurity ions assessed by solving fluid transport equations:

$$
\partial\_t n\_1 + \partial\_r \Gamma\_{1r} + \partial\_\mathbf{x} \Gamma\_{1\mathbf{x}} = k\_{ion}^0 n\_\mathcal{e} n\_0 - k\_{ion}^1 n\_\mathcal{e} n\_1 \tag{39}
$$

$$\begin{aligned} \partial\_t \Gamma\_{1, \mathbf{x}} + \partial\_r (\Gamma\_{1r} \Gamma\_{1\mathbf{x}} / n\_1) + \partial\_\mathbf{x} (\Gamma\_{1\mathbf{x}} \Gamma\_{1\mathbf{x}} / n\_1 + n\_1 \Gamma\_1 / m) &= \\ \partial\_t = -k\_{\mathrm{ion}}^1 n\_\varepsilon \Gamma\_{1, \mathbf{x}} + en\_1 E\_\mathbf{x} / m\_I \end{aligned} \tag{40}$$

where Γ1*<sup>x</sup>* and Γ1*<sup>r</sup>* are the components of the impurity ion flux density along the magnetic field and in the radial direction, respectively, *mI* is the ion mass and the parallel electric field *Ex* is determined from the electron force balance

20 Heat Transfer Studies and Applications

1022

6 7 n (m ) e c -3

1020

1019

1018

LIAS in the cooled zone at the last closed magnetic surface. where the maximum density *<sup>η</sup>*<sup>0</sup> (*t*,*r*) = *<sup>d</sup>η*<sup>0</sup>

> *Vr*, *V<sup>ρ</sup>* :

*<sup>N</sup>*<sup>0</sup> (*t*,*r*) = *dN*<sup>0</sup>

fluid transport equations:

<sup>=</sup> <sup>−</sup>*k*<sup>1</sup>

0 10 20 t ( s) m

**Figure 13.** Calculated time evolution of the electron temperature (solid curve) and density (dashed curve) induced by

region is defined by conserving the total number of particles on the magnetic surface,

*λ*<sup>0</sup> (*t*,*r*) =

*∂tn*<sup>1</sup> + *∂r*Γ1*<sup>r</sup>* + *∂x*Γ1*<sup>x</sup>* = *k*<sup>0</sup>

*ionne*Γ1,*<sup>x</sup>* + *en*1*Ex*/*mI*

*Sc* (*t*,*r*) = *πλ*<sup>2</sup>

The characteristic time for LIAS is of 10*µs*. According to figures 5b,c it is too short for a noticeable reaction of the main ions and their density *ni* is considered as unperturbed. The electron density *ne* may be, however, significantly distorted by the generation of new electrons through the ionization of impurity neutrals; due to the plasma quasi-neutrality one has *ne* = *ni* + *n*1. Here *n*<sup>1</sup> is the density of singly charged impurity ions assessed by solving

where Γ1*<sup>x</sup>* and Γ1*<sup>r</sup>* are the components of the impurity ion flux density along the magnetic field and in the radial direction, respectively, *mI* is the ion mass and the parallel electric field

The area of the cooled zone is assumed henceforth to be equal to:

 *Vr*, *V<sup>ρ</sup>*

*N*<sup>0</sup> (*t*,*r*)

*ionnen*<sup>0</sup> <sup>−</sup> *<sup>k</sup>*<sup>1</sup>

*∂t*Γ1,*<sup>x</sup>* + *∂r*(Γ1*r*Γ1*x*/*n*1) + *∂x*(Γ1*x*Γ1*x*/*n*<sup>1</sup> + *n*1*T*1/*m*) = (40)

<sup>1021</sup> <sup>40</sup>

50

6 7 T (eV) e c

30

20

10

0

and the radius of the localization

*πη*<sup>0</sup> (*t*,*r*) (37)

<sup>0</sup> (*t*,*r*) (38)

*ionnen*<sup>1</sup> (39)

$$e n\_{\ell} E\_{\mathcal{X}} = -\partial\_{\mathcal{X}} \left( n\_{\ell} T\_{\ell} \right)$$

In this study, by solving equations (39) and (40), we apply the 'shell" approximation, see references [7, 8], to find the profiles of *n*<sup>1</sup> and Γ1,*<sup>x</sup>* along the magnetic field. Due to the rapid decay of the plasma distortion induced by the outburst of impurity neutrals the temperatures of the main and impurity ions do not change noticeably, and only the change of the electron temperature is modelled by solving the TZA equations (27)-(29). The energy losses due to inelastic collisions with impurity species is assessed as

$$\mathcal{W}\_{\rm loss} = \left< n\_{\varepsilon} \right>\_{\mathcal{L}} \sum\_{j=0,1} N\_{\circ} \left( L\_{\circ}^{j} + k\_{\rm ion}^{j} E\_{\rm ion}^{j} \right),$$

where *N*1(*t*,*r*) is the total numbers of impurity singly charged species on the magnetic surface per unit length in the direction *r*, provided by the "shell" model for impurity spreading; *L<sup>j</sup> <sup>c</sup>* and *<sup>E</sup><sup>j</sup> ion* are the cooling rate and the ionization energy of the species, respectively [31]. Calculations have been performed for the conditions of LIAS in Ohmic TEXTOR discharges, see reference [4] and the radial profiles of the plasma parameters before the LIAS application are presented in figure 11. The laser radiation has been concentrated on a wall spot with fine grain graphite bulk material of 0.15*cm*<sup>2</sup> area and *<sup>ρ</sup>*<sup>0</sup> <sup>≈</sup> 0.22*cm*; typically *Ntot* <sup>∼</sup> 1017 carbon atoms were released per pulse. Time of flight measurements are well interpreted by the velocity distribution (30) with *<sup>T</sup>*<sup>0</sup> <sup>≈</sup> 1.5*eV* and *Vm* <sup>≈</sup> 8.7*km s*−1. The results of LIAS are normally quantified by measuring the total radiation emitted by *C*<sup>+</sup> species with a particular wave length [4]. To calculate this one has to apply a firm collision-radiation model that is out of scope of the present study. In figure 12 we display the cumulative energy radiated till time *t* by singly charged carbon ions, *Wrad*, calculated without and with the plasma reaction taken into account. One can see that in the latter case the rise of *Wrad* is significantly delayed but the overall emitted energy is much higher. The delay is explained by the plasma cooling through losses on radiation and ionization of impurity neutrals and significant decrease of the ion excitation rate; the larger cumulative ion radiation - by the increased electron density and growth of the ion cooling rate *L*<sup>1</sup> *<sup>c</sup>* as the temperature in the cold zone recovers due to the heat transfer from the hot zone. The time variation of *ne<sup>c</sup>* and *Te<sup>c</sup>* at *r* = *rs* = 0.465 is shown in figire 13. The results presented show that by interpreting the LIAS measurements one has to take into account distortions in the plasma parameters induced by the diagnostic itself.

#### **5. Conclusion**

Very localized injection of neutral particles of the working gas and impurities routinely happens in fusion devices, both in deliberate and accident ways. For example, by using laser-based diagnostics the particles resided in thin surface layer on the wall are released within a short intense bunch. In the plasma they are excited by electrons, emit light and, by measuring this, one can assess the total number of particles ejected and, thus, judge about the wall composition. Puffing of the working gas for the plasma density control is also done through inlets of much smaller dimensions than that of the wall. The time variation and three-dimensionality of the problem in question requires normally the application of extremely time demanding modelling tools. In the present paper we elaborate a two zone approximation, an approach which allows to reduce the problem to solving of one-dimensional equations, describing the time evolution of the radial profiles of the electron temperature values averaged over the cooled and hot zone, where energy is dissipated in direct interactions with injected particles and transferred to the cooled zone with parallel heat conduction, respectively. The elaborated TZA approach is tested by comparing its predictions with numerical solution of heat conduction equation in oneand two-dimensional configurations. In the latter case both numerical solutions and TZA demonstrate the importance of the heat flux limit by interpreting the edge plasma cooling caused by a massive gas injection in the JET tokamak. As an example of applications for realistic three-dimensional configurations the penetration of carbon atoms released into a TEXTOR deuterium plasma by a short laser pulse is modelled. It is demonstrated that for a firm interpretation of measurements with the laser-based LIAS diagnostics one has to take into account the modifications induced in the plasma by impurity particles released.

#### **Author details**

#### Mikhail Tokar

Institute for Energy and Climate Research - Plasma Physics, Research Centre Jülich, Germany

#### **References**


[8] Tokar MZ, Koltunov M. "Shell" approach to modelling of impurity spreading from localized sources in plasma. International Journal of Modelling, Simulation, and Scientific Computing 2014; 5(1) 1441005.

22 Heat Transfer Studies and Applications

**Author details**

Mikhail Tokar

**References**

2013; 53(9) 093002.

Fusion 2002; 44(8) R27-R80.

Nuclear Fusion 2013; 53(9) 093014.

control is also done through inlets of much smaller dimensions than that of the wall. The time variation and three-dimensionality of the problem in question requires normally the application of extremely time demanding modelling tools. In the present paper we elaborate a two zone approximation, an approach which allows to reduce the problem to solving of one-dimensional equations, describing the time evolution of the radial profiles of the electron temperature values averaged over the cooled and hot zone, where energy is dissipated in direct interactions with injected particles and transferred to the cooled zone with parallel heat conduction, respectively. The elaborated TZA approach is tested by comparing its predictions with numerical solution of heat conduction equation in oneand two-dimensional configurations. In the latter case both numerical solutions and TZA demonstrate the importance of the heat flux limit by interpreting the edge plasma cooling caused by a massive gas injection in the JET tokamak. As an example of applications for realistic three-dimensional configurations the penetration of carbon atoms released into a TEXTOR deuterium plasma by a short laser pulse is modelled. It is demonstrated that for a firm interpretation of measurements with the laser-based LIAS diagnostics one has to take

into account the modifications induced in the plasma by impurity particles released.

Institute for Energy and Climate Research - Plasma Physics, Research Centre Jülich, Germany

[1] McCormick K, Fiedler S, Kocsis G, Schweinzer J, Zoletnik S. Edge density measurements with a fast Li beam probe in tokamak and stellarator experiments.

[2] Hollmann EM, Jernigan TC, Groth M, et al. Measurements of impurity and heat dynamics during noble gas jet-initiated fast plasma shutdown for disruption mitigation

[4] Philipps V, Malaquias A, Hakola A, et al. Development of laser-based techniques for in situ characterization of the first wall in ITER and future fusion devices. Nuclear Fusion

[5] Greenwald M. Density limits in toroidal plasmas. Plasma Physics and Controlled

[6] de Vries PC, Rapp J, Schüller FC, Tokar MZ. Influence of Recycling on the Density

[7] Tokar MZ, Koltunov M. 'Shell' model for impurity spreading from a localized source.

Fusion Engineering and Design 1997; 34-35 125-134.

[3] Wesson J. Tokamaks. Third edition. Oxford: Claredon Press; 2004.

Limit in TEXTOR-94. Physical Review Letters 1998; 80(16) 3519-22.

in DIII-D. Nuclear Fusion 2005; 45(9) 1046-55.

