**1. Introduction**

There are number of physical situations where plasmas neutrality breaks down through boundary layers called plasma sheaths, which are either free or in contact with a wall. The plasma sheaths transition problems are at the heart of an industrial revolution whose theme is the design of matter on the molecular scale.The study of the charge separation at a plasma edge requires generally the solution of the kinetic equations of plasmas which, for a collisionless plasma, usually reduce to the well-known Vlasov equation. Some examples for the solution of the Vlasov equation for sheaths transition problems have been presented in Shoucri, 2008a, 2009a. A problem of interest is the problem involving the generation of radial electric fields and poloidal flows to achieve radial force balance at a steep density gradient in the presence of an external magnetic field. This problem is of great importance in the steep density gradients pedestal of the high confinement mode (H-mode) in tokamaks, since it largely affects the edge physics of the H-mode. In the present work, we shall study the problem of the generation of a charge separation and the associated electric field at the edge of a cylindrical plasma column, in the presence of an external magnetic field directed along the cylinder axis. In previous publications on this problem (Shoucri *et al.*, 2003, 2004, 2008b, 2009b), we have considered the case where the electrons were frozen by the magnetic field lines, with a constant density profile which changes rapidly along the gradient over an ion orbit size. Along the gradient the electrons bound by the magnetic field cannot move across this field to exactly compensate the ion charge which results from the finite ions' gyroradius. This effect is especially important for large values of the ratio / *i De* , where *<sup>i</sup>* is the ions' gyroradius and *De* is the Debye length. Accurate calculation of the charge separation is important for the accurate calculation of the self-consistent electric field. This requires also an accurate calculation of the exact ion orbits using a kinetic equation. In the present work, we use an Eulerian Vlasov code to study the

© 2012 Shoucri, licensee InTech. This is an open access chapter 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 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.

evolution of the ion distribution function for the problem of the charge separation at the edge of a cylindrical plasma, in the presence of an external magnetic field directed along the cylinder axis. Eulerian codes have the advantage of a very low noise level and make it possible to measure accurately a very small charge separation (Shoucri *et al.*, 2003, Shoucri *et al.*, 2004 & Shoucri, 2009a,b), and allows accurate results in the low density regions of the phase-space.

Charge Separation and Electric Field at a Cylindrical Plasma Edge 81

As mentioned above, the ions are described by a kinetic Vlasov equation and the electrons are described by a guiding center equation. These equations are solved numerically using a method of characteristics. There have been important advances in the last few decades in the domain of the numerical solution of hyperbolic type partial differential equations using the method of characteristics. The application of the method of characteristics for the numerical solution of the kinetic equations of plasmas and for the guiding center equations has been recently discussed in several publication (see for instance Shoucri, 2008a,c, 2009a). These methods are Eulerian methods which use a computational mesh to discretize the equations on a fixed grid, and have been successfully applied to different important problems in plasma physics involving kinetic equations, such as laser-plasma interaction (Ghizzo *et al.*, 1990, 1992, Strozzi *et al*., 2006, Shoucri, 2008d), the calculation of an electric field at a plasma edge (Shoucri, 2002, 2008b, 2009b, Shoucri *et al.*, 2000, 2003, 2004), the applications of gyro-kinetic codes to study edge physics problems in plasmas (Manfredi *et al.*, 1996, Shoucri, 2001, Shoucri *et al.*, 2005, Pohn & Shoucri, 2008) and to collisional plasmas (Batishchev *et al.*, 1999). These methods present the great advantage of having a low noise level, and allow accurate results in the low density regions of the phase-space (Shoucri, 2008c). In the applications presented (Shoucri, 2008a,c, 2009a), the computation was usually done on a fixed grid, so no dynamical grid adjustment was necessary, and interpolation was restricted to the use of a cubic spline, which compared favourably with other methods (see, for instance, Pohn *et al.*, 2005), so altogether the method was accurate and remained relatively simple. Interpolation in several dimensions using a tensor product of cubic *B*-splines has been also successfully applied (Sonnendrücker *et al.*, 1999, Shoucri, 2008a, 2009a, 2011). The method of characteristics has been also successfully applied in fluids to problems having shock wave solution (Shoucri & Shoucri, 2007). Large Courant-Frederichs-Levy (CFL) computation parameter is possible, and therefore the time-step numerical limitation by large velocities can be removed, if the physics makes it possible. A more complete study on splines can be found in the book of Ahlberg *et al.*, 1967, and an important theoretical study on the method of characteristics can be found, for instance, in the book of Abbott, 1966. More applications to plasma physics problems can be

found in the several references we have cited.

toroidal direction, and

direction is given by:

**2. The relevant equations and the numerical method** 

is the poloidal direction. The constant magnetic field in the *z* 

(1)

 *=* 

*rR rR* . We consider a deuterium plasma,

direction. The radial direction in the cylindrical plasma is normal to a vessel surface which is located at *r* = *rmax*=*R* . The constant magnetic field is in the *z* direction which represents the

> 1 cos *<sup>z</sup> B e <sup>B</sup>*

*Rmaj* is the major radius of the tokamak. We assume that at *r* = *rmax*=*R* we have / *R Rmaj* =0.2.

, ) . The plasma is uniform in the *z*

*/2,* and / *maj* 

*r R* where

We consider the cylindrical coordinate system (*r z* ,

<sup>0</sup>

In this case we can also write / 0.2 / *maj*

where *B0* is the magnetic field along the axis of the cylinder at

It was pointed out in the analysis of H-mode power threshold in a tokamak by Groebner *et al.*, 2001, that the changes in the electron density *ne* and *ne* in the transition to H-mode are small, and changes in *Te* are barely perceptible in the data. Electrons and ions have a density profile which varies rapidly along the gradient over an ion orbit size. In previous publications (Shoucri *et al.*, 2003, 2004, Shoucri, 2002, 2008a,b, 2009a,b) the electrons density and temperature profiles were kept constant. In the present work we will allow the electrons to move, and it will be sufficient for the purpose of our study to describe the motion of the electrons, having a small gyroradius, by a guiding center equation (Shoucri *et al.*, 1997). This allows a more accurate description of the contribution of the electrons, with respect to the approximation previously used in Shoucri *et al.*, 2003, 2004, 2009a, where the profile of the electrons was assumed constant in time. The electrons motion across the magnetic field in the gradient region is limited by the guiding center equation, and the magnetized electrons cannot move sufficiently across the magnetic field to compensate the charge separation which results from the ions motion due to the finite ion orbits. To determine this charge separation at the plasma edge along the gradient, it is important to calculate the ion orbits accurately by using an Eulerian Vlasov code. The larger the ions' gyroradius, the bigger the charge separation and the self-consistent electric field at the edge. (Hence the important role played by even small fractions of impurity ions, because of their large gyroradii, in enhancing the electric field at the plasma edge).

In full toroidal geometry, there are "neoclassical" effects which can play a role in this problem, such as the neoclassical enhancement of the classical ion polarization drift, or the neoclassical damping of poloidal flows (Stix, 1973, Hirshman 1978, Waltz *et. al.*, 1999). We focus for simplicity in the present work on a cylindrical geometry for the problem of the generation of an electric field and poloidal flow at a plasma edge due to the finite ions' gyroradius, when the external magnetic field is applied along the axis of the cylinder. The solution we present is a two-dimensional (2D) solution in a cylindrical geometry, with the external magnetic field assumed uniform along its axis. We compare the radial electric field calculated along the gradient with the macroscopic values calculated from the kinetic code for the gradient of the ion pressure in the radial direction, and we find that this quantity balances the radial electric field fairly well, a result similar to what has been presented in one-dimensional geometry (Shoucri, 2002, Shoucri *et al*. 2003, 2004). The contribution of the Lorentz force term along the gradient is negligible. For the parameters we use in the present work, the solution allows for a small value of *E* to exist, especially on the high field side of the cylinder as it will be discussed below, which causes a small oscillation in the radial direction.

As mentioned above, the ions are described by a kinetic Vlasov equation and the electrons are described by a guiding center equation. These equations are solved numerically using a method of characteristics. There have been important advances in the last few decades in the domain of the numerical solution of hyperbolic type partial differential equations using the method of characteristics. The application of the method of characteristics for the numerical solution of the kinetic equations of plasmas and for the guiding center equations has been recently discussed in several publication (see for instance Shoucri, 2008a,c, 2009a). These methods are Eulerian methods which use a computational mesh to discretize the equations on a fixed grid, and have been successfully applied to different important problems in plasma physics involving kinetic equations, such as laser-plasma interaction (Ghizzo *et al.*, 1990, 1992, Strozzi *et al*., 2006, Shoucri, 2008d), the calculation of an electric field at a plasma edge (Shoucri, 2002, 2008b, 2009b, Shoucri *et al.*, 2000, 2003, 2004), the applications of gyro-kinetic codes to study edge physics problems in plasmas (Manfredi *et al.*, 1996, Shoucri, 2001, Shoucri *et al.*, 2005, Pohn & Shoucri, 2008) and to collisional plasmas (Batishchev *et al.*, 1999). These methods present the great advantage of having a low noise level, and allow accurate results in the low density regions of the phase-space (Shoucri, 2008c). In the applications presented (Shoucri, 2008a,c, 2009a), the computation was usually done on a fixed grid, so no dynamical grid adjustment was necessary, and interpolation was restricted to the use of a cubic spline, which compared favourably with other methods (see, for instance, Pohn *et al.*, 2005), so altogether the method was accurate and remained relatively simple. Interpolation in several dimensions using a tensor product of cubic *B*-splines has been also successfully applied (Sonnendrücker *et al.*, 1999, Shoucri, 2008a, 2009a, 2011). The method of characteristics has been also successfully applied in fluids to problems having shock wave solution (Shoucri & Shoucri, 2007). Large Courant-Frederichs-Levy (CFL) computation parameter is possible, and therefore the time-step numerical limitation by large velocities can be removed, if the physics makes it possible. A more complete study on splines can be found in the book of Ahlberg *et al.*, 1967, and an important theoretical study on the method of characteristics can be found, for instance, in the book of Abbott, 1966. More applications to plasma physics problems can be found in the several references we have cited.

### **2. The relevant equations and the numerical method**

80 Numerical Simulation – From Theory to Industry

enhancing the electric field at the plasma edge).

work, the solution allows for a small value of *E*

direction.

phase-space.

evolution of the ion distribution function for the problem of the charge separation at the edge of a cylindrical plasma, in the presence of an external magnetic field directed along the cylinder axis. Eulerian codes have the advantage of a very low noise level and make it possible to measure accurately a very small charge separation (Shoucri *et al.*, 2003, Shoucri *et al.*, 2004 & Shoucri, 2009a,b), and allows accurate results in the low density regions of the

It was pointed out in the analysis of H-mode power threshold in a tokamak by Groebner *et al.*, 2001, that the changes in the electron density *ne* and *ne* in the transition to H-mode are small, and changes in *Te* are barely perceptible in the data. Electrons and ions have a density profile which varies rapidly along the gradient over an ion orbit size. In previous publications (Shoucri *et al.*, 2003, 2004, Shoucri, 2002, 2008a,b, 2009a,b) the electrons density and temperature profiles were kept constant. In the present work we will allow the electrons to move, and it will be sufficient for the purpose of our study to describe the motion of the electrons, having a small gyroradius, by a guiding center equation (Shoucri *et al.*, 1997). This allows a more accurate description of the contribution of the electrons, with respect to the approximation previously used in Shoucri *et al.*, 2003, 2004, 2009a, where the profile of the electrons was assumed constant in time. The electrons motion across the magnetic field in the gradient region is limited by the guiding center equation, and the magnetized electrons cannot move sufficiently across the magnetic field to compensate the charge separation which results from the ions motion due to the finite ion orbits. To determine this charge separation at the plasma edge along the gradient, it is important to calculate the ion orbits accurately by using an Eulerian Vlasov code. The larger the ions' gyroradius, the bigger the charge separation and the self-consistent electric field at the edge. (Hence the important role played by even small fractions of impurity ions, because of their large gyroradii, in

In full toroidal geometry, there are "neoclassical" effects which can play a role in this problem, such as the neoclassical enhancement of the classical ion polarization drift, or the neoclassical damping of poloidal flows (Stix, 1973, Hirshman 1978, Waltz *et. al.*, 1999). We focus for simplicity in the present work on a cylindrical geometry for the problem of the generation of an electric field and poloidal flow at a plasma edge due to the finite ions' gyroradius, when the external magnetic field is applied along the axis of the cylinder. The solution we present is a two-dimensional (2D) solution in a cylindrical geometry, with the external magnetic field assumed uniform along its axis. We compare the radial electric field calculated along the gradient with the macroscopic values calculated from the kinetic code for the gradient of the ion pressure in the radial direction, and we find that this quantity balances the radial electric field fairly well, a result similar to what has been presented in one-dimensional geometry (Shoucri, 2002, Shoucri *et al*. 2003, 2004). The contribution of the Lorentz force term along the gradient is negligible. For the parameters we use in the present

the cylinder as it will be discussed below, which causes a small oscillation in the radial

to exist, especially on the high field side of

We consider the cylindrical coordinate system (*r z* , , ) . The plasma is uniform in the *z* direction. The radial direction in the cylindrical plasma is normal to a vessel surface which is located at *r* = *rmax*=*R* . The constant magnetic field is in the *z* direction which represents the toroidal direction, and is the poloidal direction. The constant magnetic field in the *z*  direction is given by:

$$\vec{B} = \frac{B\_0 \vec{e}\_z}{1 + \varepsilon \cos \theta} \tag{1}$$

where *B0* is the magnetic field along the axis of the cylinder at  *= /2,* and / *maj r R* where *Rmaj* is the major radius of the tokamak. We assume that at *r* = *rmax*=*R* we have / *R Rmaj* =0.2. In this case we can also write / 0.2 / *maj rR rR* . We consider a deuterium plasma,

/ 2x1836 *m m i e* . The ions are described by the normalized 2D Vlasov equation for the ions distribution function ( , , , , ) : *i r f r vvt* 

$$\begin{aligned} \frac{\partial \hat{f}\_i}{\partial t} + \upsilon\_r \frac{\partial \hat{f}\_i}{\partial r} + \frac{\upsilon\_\theta}{r} \frac{\partial \hat{f}\_i}{\partial \theta} + \left( E\_r + \alpha\_{ci} \upsilon\_\theta + \frac{\nu\_\theta^2}{r} \right) \frac{\partial \hat{f}\_i}{\partial \upsilon\_r} \\ + \left( E\_\theta - \alpha\_{ci} \upsilon\_r - \frac{\upsilon\_r \upsilon\_\theta}{r} \right) \frac{\partial \hat{f}\_i}{\partial \upsilon\_\theta} = 0 \end{aligned} \tag{2}$$

where from Eq.(1) the ion cyclotron frequency <sup>0</sup> / (1 cos ) *ci ci* . The electrons are described by the normalized 2D guiding-center equation:

$$\frac{\partial n\_e}{\partial t} + \nabla.\vec{V}\_d n\_e = 0\tag{3}$$

Charge Separation and Electric Field at a Cylindrical Plasma Edge 83

*ci pi* . With the initial distribution

spatially constant, the factor 2*Ti* in

(7)

0

/

<sup>2</sup> <sup>1</sup> 1 ; <sup>20</sup>

 

*e De e ci pi*

*i ii*

 

Eq. (7) in the calculation of the gyroradius takes into account the fact that the perpendicular

Equations (2) and (3) are solved by a method of fractional step (Yanenko, 1971), first applied for a Vlasov equation by Cheng & Knorr, 1976, and Gagné & Shoucri, 1977, coupled to a method of characteristics. For the general case where several dimensions are involved, the fractional step technique allows the reduction of the multi-dimensional equation to an equivalent set of reduced equations. To advance Eq. (2) and Eq.(3) for one time-step *t*, the

> *r f f vf <sup>v</sup> t rr*

*t r*

below. We then solve Poisson's equation in Eq.(5) to calculate the new electric field.

*r ci ci r r <sup>f</sup> v f vv f E v E v t r v r v*

 

2

*<sup>n</sup> n E n E tr r*

Again Eq.(10) is solved by a method of characteristics, to be described below.

 

 

*dr d nnn V V*

Equations (8) and (9) are solved by interpolation along their characteristics, to be described

0 *<sup>i</sup> <sup>i</sup> r i*

1 cos 1 1 cos *<sup>e</sup> e e r*

**Step 3.** We repeat **Step1** for *t*/2, and then solve Eq.(5) to calculate the new electric field

   

0 0

 

(11)

*ci ci*

 

(8)

(9)

 

 

(10)

In the present calculations the parameter 0 / 0.1 / 2

.

0 *ii i*

0 *ee e*

For the solution of Eq. (8), we first solve the characteristic equations:

*i T*

*m*

function for the ions in Eq. (6), we have *i ir i TTT*

*r*

*v vv*

splitting of the equations is applied as follows.

**Step 1.** We solve for *t*/2 the equations:

**Step 2.** We solve for *t* the equations:

*n* 1 *<sup>r</sup> <sup>E</sup>* , *<sup>n</sup>* <sup>1</sup> *<sup>E</sup>*

.

**2.1. The solution for Step 1** 

This completes a one time-step cycle *t*.

temperature is 2 22 <sup>2</sup> *<sup>i</sup>*

*T T T T* 

where the drift velocity <sup>2</sup> x / *V EB B <sup>d</sup>* . Equation (3) can be developed in our 2D system to give the following equation:

$$\frac{\partial n\_{\varepsilon}}{\partial t} + V\_{dr} \frac{\partial n\_{\varepsilon}}{\partial r} + V\_{d\theta} \frac{\partial n\_{\varepsilon}}{\partial \theta} = -n\_{\varepsilon}E\_{\theta} \frac{\partial}{\partial r} \frac{1 + \varepsilon \cos \theta}{o o\_{\varepsilon 0}} + n\_{\varepsilon}E\_{r} \frac{1}{r} \frac{\partial}{\partial \theta} \frac{1 + \varepsilon \cos \theta}{o o\_{\varepsilon 0}} \tag{4}$$

where 0 1 cos *dr ci V E* ; 0 1 cos *<sup>r</sup> d ci <sup>E</sup> <sup>V</sup> r* 

This system is coupled to Poisson's equation:

$$\frac{1}{r}\frac{\partial}{\partial r}(r\frac{\partial \phi}{\partial r}) + \frac{1}{r^2}\frac{\partial^2 \phi}{\partial \theta^2} = -(n\_i - n\_e);\ E\_r = -\frac{\partial \phi}{\partial r};\ E\_\theta = -\frac{1}{r}\frac{\partial \phi}{\partial \theta} \tag{5}$$

Time is normalized to the inverse plasma frequency <sup>1</sup> *pi* , velocity is normalized to the acoustic velocity / *C Tm s ei* , and length is normalized to the Debye length / *De s pi C* . *Te* is the electron temperature and *mi* is the ion mass. The potential is normalized to / *<sup>e</sup> T e* , and the electric field is normalized to / ( ) *e De T e* , and the density is normalized to the peak initial central density. The ion cyclotron frequency *ci* as previously defined is normalized to *pi* . The system is solved over a length *Lr* = 175 Debye lengths in front of the vessel surface, with an initial ion distribution function for the deuterons over the domain *r RrR <sup>L</sup>* , given in our normalized units, by:

$$f\_i(r\_r \, \theta \, \upsilon\_{r\,\,\,\,} \, \upsilon\_\theta) = n\_i(r) \frac{e^{-\left(v\_r^2 + v\_\theta^2\right)/2 \, T\_i}}{2\pi T\_i} \; ; \tag{6}$$

( ) *<sup>i</sup> n r* and ( ) *<sup>e</sup> n r* are the initial ion and electron density profiles respectively, with () () *e i nr nr* in the initially neutral system. *R* is the radius of the cylinder as we previously mentioned. We also use the following parameters:

Charge Separation and Electric Field at a Cylindrical Plasma Edge 83

$$\frac{T\_i}{T\_e} = 1 \; ; \quad \frac{\rho\_i}{\lambda\_{De}} = \sqrt{\frac{2T\_i}{T\_e}} \frac{1}{o\_{ci0} / \; o\_{pi}} = 20\tag{7}$$

In the present calculations the parameter 0 / 0.1 / 2 *ci pi* . With the initial distribution function for the ions in Eq. (6), we have *i ir i TTT* spatially constant, the factor 2*Ti* in Eq. (7) in the calculation of the gyroradius takes into account the fact that the perpendicular temperature is 2 22 <sup>2</sup> *<sup>i</sup> r i T v vv m* .

Equations (2) and (3) are solved by a method of fractional step (Yanenko, 1971), first applied for a Vlasov equation by Cheng & Knorr, 1976, and Gagné & Shoucri, 1977, coupled to a method of characteristics. For the general case where several dimensions are involved, the fractional step technique allows the reduction of the multi-dimensional equation to an equivalent set of reduced equations. To advance Eq. (2) and Eq.(3) for one time-step *t*, the splitting of the equations is applied as follows.

**Step 1.** We solve for *t*/2 the equations:

82 Numerical Simulation – From Theory to Industry

distribution function ( , , , , ) : *i r*

give the following equation:

0 1 cos

This system is coupled to Poisson's equation:

*rr r r* 

and the electric field is normalized to / ( ) *e De T e*

*r RrR <sup>L</sup>* , given in our normalized units, by:

mentioned. We also use the following parameters:

initial central density. The ion cyclotron frequency

*ci*

;

where

to 

( )

*dr*

*V E* *f r vvt* 

described by the normalized 2D guiding-center equation:

. 0 *<sup>e</sup>*

/ 2x1836 *m m i e* . The ions are described by the normalized 2D Vlasov equation for the ions

*ii i i r r ci*

*<sup>f</sup> f vf <sup>f</sup> v E t rr r v*

 

*t* 

0 0

*r*

2 2 2 1 1 ( ) ( ); *i e <sup>r</sup> n n*

 

*<sup>E</sup> <sup>V</sup>*

*d*

Time is normalized to the inverse plasma frequency <sup>1</sup>

*nnn V V nE n E t r r r*

1 cos *<sup>r</sup>*

where from Eq.(1) the ion cyclotron frequency <sup>0</sup> / (1 cos )

*ci r*

*d e <sup>n</sup> V n*

where the drift velocity <sup>2</sup> x / *V EB B <sup>d</sup>* . Equation (3) can be developed in our 2D system to

1 cos 1 1 cos *ee e dr d e e r*

  ; *<sup>r</sup> E r* 

> *pi*

*r i vv T*

*i*

(4)

 

0

acoustic velocity / *C Tm s ei* , and length is normalized to the Debye length /

*Te* is the electron temperature and *mi* is the ion mass. The potential is normalized to / *<sup>e</sup> T e* ,

*pi* . The system is solved over a length *Lr* = 175 Debye lengths in front of the vessel surface, with an initial ion distribution function for the deuterons over the domain

 2 2 /2 , , , () ; <sup>2</sup>

*<sup>i</sup> n r* and ( ) *<sup>e</sup> n r* are the initial ion and electron density profiles respectively, with () () *e i nr nr* in the initially neutral system. *R* is the radius of the cylinder as we previously

*<sup>e</sup> f r v v nr <sup>T</sup>*

*ir i*

*ci*

 

*vv f E v*

 

  2

 

*r i*

*r v*

*ci ci* 

*ci ci*

 

<sup>1</sup> *<sup>E</sup> r*

0

 

(3)

 

, and the density is normalized to the peak

(6)

(5)

*De s pi C* .

, velocity is normalized to the

*ci* as previously defined is normalized

*r*

(2)

. The electrons are

$$\frac{\partial f\_i}{\partial t} + \upsilon\_r \frac{\partial f\_i}{\partial r} + \frac{\upsilon\_\theta}{r} \frac{\partial f\_i}{\partial \theta} = 0 \tag{8}$$

$$\frac{\partial \mathfrak{n}\_{\varepsilon}}{\partial t} + V\_{dr} \frac{\partial \mathfrak{n}\_{\varepsilon}}{\partial r} + V\_{d\theta} \frac{\partial \mathfrak{n}\_{\varepsilon}}{\partial \theta} = 0 \tag{9}$$

Equations (8) and (9) are solved by interpolation along their characteristics, to be described below. We then solve Poisson's equation in Eq.(5) to calculate the new electric field.

**Step 2.** We solve for *t* the equations:

$$\frac{\partial f\_i}{\partial t} + \left(E\_r + v\_\theta \,\,\phi v\_{ci} + \frac{v\_\theta^2}{r}\right) \frac{\partial f\_i}{\partial v\_r} + \left(E\_\theta - \alpha\_{ci} \,\, v\_r - \frac{v\_r v\_\theta}{r}\right) \frac{\partial f\_i}{\partial v\_\theta} = 0 \tag{10}$$

$$\frac{\partial n\_e}{\partial t} = -n\_e E\_\theta \frac{\partial}{\partial r} \frac{1 + \varepsilon \cos \theta}{o o\_{ci0}} + n\_e E\_r \frac{1}{r} \frac{\partial}{\partial \theta} \frac{1 + \varepsilon \cos \theta}{o o\_{ci0}} \tag{11}$$

Again Eq.(10) is solved by a method of characteristics, to be described below.

**Step 3.** We repeat **Step1** for *t*/2, and then solve Eq.(5) to calculate the new electric field *n* 1 *<sup>r</sup> <sup>E</sup>* , *<sup>n</sup>* <sup>1</sup> *<sup>E</sup>*.

This completes a one time-step cycle *t*.

#### **2.1. The solution for Step 1**

For the solution of Eq. (8), we first solve the characteristic equations:

$$\frac{dr}{dt} = \begin{array}{c} \upsilon\_r \\ \end{array}; \qquad \qquad \qquad \frac{d\theta}{dt} = \frac{\upsilon\_\theta}{r} \tag{12}$$

Charge Separation and Electric Field at a Cylindrical Plasma Edge 85

. The following

is at the grid point *j*

(20)

 

 

(18)

(19)

is the point where the characteristic is originating at *nt* (not necessarily a

2 *j nt* 

( )

> 

(21)

(22)

, and which are solved by iteration. Usually, two

0

0

*r*

*r*

*j r*

*ci*

*r*

*j j rj j j ci*

*r*

1 cos (, )

 

1 cos

*j*

*r*

*r*

*j*

(23)

1

0

*ci*

*j*

1 cos( )

1

(24)

. We get at the first

Equations (17) are solved using an iterative process. We assume that at the time-step

*j* and

*r r j n n n j nj n dr dr r rt r rt t*

*<sup>r</sup> j n n n j nj n d d*

*t r rt t*

 

 

(, ) <sup>4</sup>

4

*r*

*j j*

*r r*

*<sup>t</sup> <sup>R</sup> E r*

*E r <sup>t</sup> <sup>R</sup> r*

 

 

(25)

1/2 / 2 *n n t tt* , *r* is at the grid point *<sup>r</sup>*

*t*

*t* 

 

we can rewrite Eqs.(18,19) as follows:

which are implicit equations for ( , ) *<sup>r</sup>*

<sup>1</sup> (, ) <sup>4</sup> *<sup>r</sup> r dr j j*

<sup>1</sup> (, ) <sup>4</sup> *<sup>r</sup> dj j*

where ( ( ), ( )) *n n rt t* 

grid point). Put:

iteration:

Second iteration:

leapfrog scheme can be written for the solution of Eq.(17):

1/4 1/4 ( ) () () (, )( , ) / 2 2 2

1/4 1/4 ( ) () () (, )( , ) / 2 2 2

*r j n*

> *t V r*

*t V r*

*r r rt*

(, ) <sup>4</sup> *<sup>r</sup> r dr j r j*

(, ) <sup>4</sup> *<sup>r</sup> d j rj*

*t V r*

*t V r*

 

2 11 1 1

*r r*

*r dr j r j j rj*

(, ) (, ) 4 4

*t t <sup>R</sup> V r E r*

iterations are sufficient for convergence. We start with 0 0 ( 0, 0) *<sup>r</sup>*

*V r V*

*V r V*

( ) ; <sup>2</sup>

at a given ( , ) *<sup>r</sup> v v* in velocity space. The solution of Eqs (12) originating at ( , *o o r* ) at time *t,*  and reaching the grid point (*r*, ) at *tt t* / 2 can be written as follows, for a half timestep *t/2*:

$$r\_o = r - \upsilon\_r \Delta t / 2 ; \; \theta\_o = \theta - \frac{\upsilon\_\theta}{\upsilon\_r} \ln \frac{r}{r - \upsilon\_r \Delta t / 2} \tag{13}$$

For *vr t*/2 « 1, the second equation reduces to / 2 *<sup>o</sup> v t r* . Therefore the solution of Eq. (8) can be written, for half a time-step *t*/2, as follows:

$$f\_i^\* \left( r\_r \theta\_r v\_{r'} v\_{\theta'}, t + \Delta t / 2 \right) = f\_i \left( r - v\_r \, \Delta t / 2, \theta - \frac{v\_\theta}{v\_r} \ln \frac{r}{r - v\_r \Delta t / 2}, v\_{r'}, v\_{\theta'}, t \right) \tag{14}$$

The right hand side of Eq. (14) is calculated by interpolation using a tensor product of cubic *B*-splines, in which is periodic (Shoucri *et al.,* 2004, Shoucri, 2008a, 2009a). For each *<sup>r</sup> v* and *v* , we write for the interpolation function *s r*(, ) :

$$\text{cs}(r,\theta) = \sum\_{j=0}^{N\_r} \sum\_{k=0}^{N\_\theta} \eta\_{jk} B\_j(r) B\_k(\theta) \tag{15}$$

taking into account that is periodic. ( ) *<sup>j</sup> B r* and ( ) *<sup>k</sup> B* are cubic splines.The cubic *B*-spline is defined as:

$$B\_{j}(\mathbf{x}) = \frac{1}{6h^{3}} \begin{cases} (\mathbf{x} - \mathbf{x}\_{j})^{3} & \mathbf{x}\_{j} \le \mathbf{x} < \mathbf{x}\_{j+1} \\ h^{3} + 3h^{2}(\mathbf{x} - \mathbf{x}\_{j+1}) + 3h(\mathbf{x} - \mathbf{x}\_{j+1})^{2} \\ & -3(\mathbf{x} - \mathbf{x}\_{j+1})^{3} \\ h^{3} + 3h^{2}(\mathbf{x}\_{j+3} - \mathbf{x}) + 3h(\mathbf{x}\_{j+3} - \mathbf{x})^{2} \\ & -3(\mathbf{x}\_{j+3} - \mathbf{x})^{3} \\ & -3(\mathbf{x}\_{j+3} - \mathbf{x})^{3} & \mathbf{x}\_{j+2} \le \mathbf{x} < \mathbf{x}\_{j+3} \\ (\mathbf{x}\_{j+4} - \mathbf{x})^{3} & \mathbf{x}\_{j+3} \le \mathbf{x} < \mathbf{x}\_{j+4} \end{cases} \tag{16}$$

and () 0 *<sup>j</sup> B x* otherwise. The grid size *h* (either *r* or in our notation), is assumed uniform. For the calculation of the coefficients *jk* of the cubic *B*-spline interpolation function *s r*(, ) for periodic in Eq.(15) see details in Shoucri *et al.,* 2004 and Shoucri, 2008a, 2009a.

We then solve Eq.(9) for a half time-step *t/2* along the characteristics:

$$\frac{dr}{dt} = V\_{dr}(r, \theta) = E\_{\theta} \frac{1 + \varepsilon \cos \theta}{o o\_{ci0}} \ ; \qquad \qquad \frac{d\theta}{dt} = V\_{d\theta}(r, \theta) = -\frac{E\_r}{r} \frac{1 + \varepsilon \cos \theta}{o o\_{ci0}} \tag{17}$$

Equations (17) are solved using an iterative process. We assume that at the time-step 1/2 / 2 *n n t tt* , *r* is at the grid point *<sup>r</sup> j* and is at the grid point *j* . The following leapfrog scheme can be written for the solution of Eq.(17):

$$\frac{r\_{j\_r} - r(t\_n)}{\Delta t / 2} = V\_{dr}(r^{n+1/4}, \theta^{n+1/4}) = V\_{dr}(\frac{r\_{j\_r} + r(t\_n)}{2}, \frac{\theta\_{j\_\theta} + \theta(t\_n)}{2}) \tag{18}$$

$$\frac{\theta\_{j\_\theta} - \theta(t\_n)}{\Delta t / 2} = V\_{d\theta}(r^{n+1/4}, \theta^{n+1/4}) = V\_{d\theta}(\frac{r\_{j\_r} + r(t\_n)}{2}, \frac{\theta\_{j\_\theta} + \theta(t\_n)}{2}) \tag{19}$$

where ( ( ), ( )) *n n rt t* is the point where the characteristic is originating at *nt* (not necessarily a grid point). Put:

$$\Delta\_r = \frac{r\_{j\_r} - r(t\_n)}{2}; \ \Delta\_\theta = \frac{\theta\_{j\_\theta} - \theta(t\_n)}{2} \tag{20}$$

we can rewrite Eqs.(18,19) as follows:

84 Numerical Simulation – From Theory to Industry

and reaching the grid point (*r*,

/ 2; *o r*

at a given ( , ) *<sup>r</sup> v v*

step *t/2*:

; *r*

*v*

*rr t* 

For *vr t*/2 « 1, the second equation reduces to / 2 *<sup>o</sup>*

Eq. (8) can be written, for half a time-step *t*/2, as follows:

, we write for the interpolation function *s r*(, )

( )

 

For the calculation of the coefficients *jk*

3 3 2

3

*B*-splines, in which

taking into account that

*j*

*B x*

defined as:

periodic

*v* *dr d v*

/ 2 *<sup>o</sup> r r*

 

 

> 

*jk j k*

3

3

4 3 4

1 cos 1 cos (, ) ; (, ) *<sup>r</sup>*

 

*x*

*j j j*

*x x x xx*

(12)

(13)

. Therefore the solution of

) at *tt t* / 2 can be written as follows, for a half time-

*r r t*

*r r*

*v rvt*

(15)

are cubic splines.The cubic *B*-spline is

in our notation), is assumed uniform.

 

 

1

 

*v t r* 

is periodic (Shoucri *et al.,* 2004, Shoucri, 2008a, 2009a). For each *<sup>r</sup> v* and

 

1 12

)

of the cubic *B*-spline interpolation function *s r*(, )

*j jj*

3 2 3

*j j j*

*x x x xx*

0 0

   

*ci ci*

(17)

*x x x xx*

*r* 

 ) at time *t,* 

(14)

(16)

for

*dt dt r*

ln

*<sup>v</sup> <sup>r</sup> <sup>f</sup> r vvt t f rv t vvt*

The right hand side of Eq. (14) is calculated by interpolation using a tensor product of cubic

:

*<sup>j</sup> B r* and ( ) *<sup>k</sup> B*

*j j j*

*x x xxx*

3

in Eq.(15) see details in Shoucri *et al.,* 2004 and Shoucri, 2008a, 2009a.

0 0 (, ) () ( ) *N Nr*

*j k s r B rB* 

1 1

*j j*

<sup>2</sup>

is periodic. ( )

3 2 2

*h h x x hx x*

3 ( ) 3( )

*j j*

*dr d*

*dr <sup>d</sup> <sup>E</sup> Vr E V r dt dt r*

We then solve Eq.(9) for a half time-step *t/2* along the characteristics:

3( )

3

3( ) 1 ( ) <sup>6</sup> 3 ( ) 3(

*h h h x x hx*

3

( )

and () 0 *<sup>j</sup> B x* otherwise. The grid size *h* (either *r* or

 

\* , , , , /2 / 2, ln , , , / 2 *i r i r <sup>r</sup>*

in velocity space. The solution of Eqs (12) originating at ( , *o o*

$$
\Delta\_r = \frac{\Delta t}{4} V\_{dr} (r\_{\dot{j}\_r} - \Delta\_{r'}, \theta\_{\dot{j}\_\theta} - \Delta\_\theta) \tag{21}
$$

$$
\Delta\_{\theta} = \frac{\Delta t}{4} V\_{d\theta} (r\_{\dot{j}\_r} - \Delta\_{r'} \theta\_{\dot{j}\_\theta} - \Delta\_\theta) \tag{22}
$$

which are implicit equations for ( , ) *<sup>r</sup>* , and which are solved by iteration. Usually, two iterations are sufficient for convergence. We start with 0 0 ( 0, 0) *<sup>r</sup>* . We get at the first iteration:

$$\Delta\_r^1 = \frac{\Delta t}{4} V\_{dr}(r\_{j\_r}, \theta\_{j\_\theta}) = \frac{\Delta t}{4} E\_\theta(r\_{j\_r}, \theta\_{j\_\theta}) \frac{1 + \frac{r\_{j\_r}}{R} \cos \theta\_{j\_\theta}}{a \theta\_{\text{cl}0}} \tag{23}$$

$$
\Delta\_{\theta}^{1} = \frac{\Delta t}{4} V\_{d\theta} \langle r\_{\dot{j}\_r}, \theta\_{\dot{j}\_\theta} \rangle = -\frac{\Delta t}{4} \frac{E\_r(r\_{\dot{j}\_r}, \theta\_{\dot{j}\_\theta})}{r\_{\dot{j}\_r}} \frac{1 + \frac{r\_{\dot{j}\_r}}{R} \cos \theta\_{\dot{j}\_\theta}}{a o\_{d0}} \tag{24}
$$

Second iteration:

$$\Delta\_{r}^{2} = \frac{\Delta t}{4} V\_{dr} (r\_{\dot{j}\_{r}} - \Delta\_{r}^{1}, \theta\_{\dot{j}\_{\theta}} - \Delta\_{\theta}^{1}) = \frac{\Delta t}{4} E\_{\theta} (r\_{\dot{j}\_{r}} - \Delta\_{r}^{1}, \theta\_{\dot{j}\_{\theta}} - \Delta\_{\theta}^{1}) \frac{1 + \frac{r\_{\dot{j}\_{r}} - \Delta\_{r}^{1}}{R} \cos(\theta\_{\dot{j}\_{\theta}} - \Delta\_{\theta}^{1})}{o \nu\_{\dot{c}0}} \tag{25}$$

$$\Delta\_{\theta}^{2} = \frac{\Delta t}{4} V\_{d\theta} (r\_{\dot{j}\_r} - \Delta\_{r'}^{1} \theta\_{\dot{j}\_\theta} - \Delta\_{\theta}^{1}) = -\frac{\Delta t}{4} \frac{E\_r (r\_{\dot{j}\_r} - \Delta\_{r'}^{1} \theta\_{\dot{j}\_\theta} - \Delta\_{\theta}^{1})}{r\_{\dot{j}\_r} - \Delta\_{r}^{1}} \frac{1 + \frac{r\_{\dot{j}\_r} - \Delta\_{r}^{1}}{R} \cos(\theta\_{\dot{j}\_\theta} - \Delta\_{\theta}^{1})}{\alpha \theta\_{\dot{c}0}} \tag{26}$$

The quantities 1 1 (, ) *<sup>r</sup> rj r j E r* and 1 1 (, ) *<sup>r</sup> j rj E r* are calculated by a cubic *B*-spline interpolation, similar to what is described in Eq.(15). Finally to advance Eq.(9) by / 2 *t* , we calculate:

$$n\_e^\* \{ r\_{\dot{j}\_r}, \theta\_{\dot{j}\_\theta}, t + \Delta t \,/\2 \} = n\_e \{ r\_{\dot{j}\_r} - 2\Delta\_r^2, \theta\_{\dot{j}\_\theta} - 2\Delta\_\theta^2, t \} \tag{27}$$

Charge Separation and Electric Field at a Cylindrical Plasma Edge 87

 

*<sup>e</sup> <sup>n</sup>* , we then repeat **Step 1** for the

by solving Poisson

:

*r*

0 0

 

(33)

(, ) *<sup>n</sup> E r* 

*<sup>e</sup> n* calculated from Eq.(30) and Eq.(33). The

(34)

 

0

*t*

*ci ci*

*r r*

*r* 

 

*ci r*

 

is calculated by integrating Eq.(11) as follows:

1 cos 1 cos <sup>1</sup> ( , , ) exp( ) *e e <sup>r</sup>*

*r r*

*i*

(, ) *<sup>n</sup> <sup>r</sup> E r* and <sup>1</sup>

*f* and *<sup>n</sup>* <sup>1</sup>

*R R nr t t n E t E t*

*f* and \*\*

Equation (5) is solved by first Fourier transforming the equation in the periodic direction

 (, ) () *im m m r re*

 

and the resulting equation is then discretized in the radial direction over a uniform grid following the method of Knorr *et al*., 1980 (also rediscussed more recently in Crouseilles *et al.*

> 2 , 1 , ,1 ,1 , ,1 ( ) ( 10 ) <sup>12</sup> *m mj m mj m mj mj mj mj*

 

, ,, 22 2 ( ) 10 ( ) ( ) 1 ;2 ;1 2 2 <sup>12</sup> <sup>12</sup> <sup>12</sup> *m j m j m j*

*rm r m r rm r A BC*

 

Equation (35) is solved using a tridiagonal algorithm. To get the boundary conditions, we assume in the present calculation that the deuterons and electrons currents hitting the cylindrical wall surface at *r R* are collected by a floating potential cylindrical vessel.

(, ) ( ) or ( , )( )

*E r <sup>J</sup> <sup>J</sup> E r J J dt <sup>t</sup>*

*j j jj j*

*r r rr r* 

> 

 

(35)

2 2 2 2 2 2

(36)

*ri re r ri re rR rR rR rR rR*

(37)

where

where \*

equation.

2009):

where

where

and , ( ) *mj m j* 

 

*r*

1

*<sup>t</sup> a E*

The density \*\*(, , ) *<sup>e</sup> nr t t*

2 *r ci*

\*\* \*

solution of Eqs. (8,9) for *t* / 2 using \*\*

**2.3. The solution of Poisson's equation** 

 *r* ; *i e* 

*r R*

electric field is then updated to calculate <sup>1</sup>

2

; <sup>1</sup>

2

*i*

 

*<sup>r</sup> A BC*

*n n* ; ( , , , ) *ii r r n fr dd*

 

Therefore we have the relations for the charge collected:

*<sup>t</sup> b E*

*r* 

*<sup>e</sup> <sup>n</sup>* is calculated in Eq.(27). To calculate *<sup>n</sup>* <sup>1</sup>

 

where the right hand side of Eq.(27) is calculated again using a cubic *B*-spline interpolation as indicated in Eq.(15).

#### **2.2. The solution for Step 2**

We go to **Step2** and solve for *t* Eqs.(10,11). We first solve the equation:

$$\frac{\partial f\_i}{\partial t} + \left(E\_r + \alpha\_{ci}\nu\_\theta + \frac{\nu\_\theta^2}{r}\right)\frac{\partial f\_i}{\partial \upsilon\_r} + \left(E\_\theta - \alpha\_{ci}\upsilon\_r - \frac{\upsilon\_r\upsilon\_\theta}{r}\right)\frac{\partial f\_i}{\partial \upsilon\_\theta} = 0\tag{28}$$

The electric field *E* remained small, since <sup>1</sup> *<sup>E</sup> r* , and *r* ~ *R* =5000 towards the edge in our calculation.

$$\frac{dv\_r}{dt} = E\_r + v\_\theta \alpha\_{ci} + \frac{v\_\theta^2}{r}; \qquad \frac{dv\_\theta}{dt} = E\_\theta - \alpha\_{ci} v\_r - \frac{v\_r v\_\theta}{r} \tag{29}$$

Following the same iterative steps as described for Eq.(27), to an order <sup>2</sup> *O t* the solution of Eq. (29) yields the following solution to Eq.(28):

$$f\_i^{\star\star} \left( r, \theta, \upsilon\_r, \upsilon\_{\theta'}, t + \Delta t \right) = f\_i^{\star} \left( r, \theta, \upsilon\_r - 2a, \upsilon\_{\theta} - 2b, t \right) \tag{30}$$

Again the 2D interpolation in Eq. (30) is done using a tensor product of cubic *B*-spline, as explained in Eq.(15), and *a* and *b* are calculated with similar iterative steps as in Eqs.(23-26) and are given, to an order <sup>2</sup> *O t* , by the expressions:

$$a = \frac{\Delta t}{2} \left\{ E\_r + \alpha\_{ci} \left( v\_\vartheta - b^1 \right) + \frac{(\nu\_\vartheta - b^1)^2}{r} \right\} \tag{31}$$

$$b = \frac{\Delta t}{2} \left\{ E\_{\theta} - \alpha \iota\_{ci} \left( \upsilon\_r - a^1 \right) - \frac{(\upsilon\_{\theta} - b^1)(\upsilon\_r - a^1)}{r} \right\} \tag{32}$$

$$\text{where } a^1 = \frac{\Delta t}{2} \left\{ E\_r + \alpha\_{ci} \nu\_\theta + \frac{\nu\_\theta^2}{r} \right\}; \ b^1 = \frac{\Delta t}{2} \left\{ E\_\theta - \alpha\_{ci} \nu\_r - \frac{\nu\_\theta \nu\_r}{r} \right\}$$

The density \*\*(, , ) *<sup>e</sup> nr t t* is calculated by integrating Eq.(11) as follows:

$$n\_{\epsilon}^{\prime \ast} \left( r\_{\nu} \theta , t + \Delta t \right) = n\_{\epsilon}^{\ast} \exp(-E\_{\theta} \frac{\partial}{\partial r} \frac{1 + \frac{r}{R} \cos \theta}{o \epsilon\_{\mathrm{c0}0}} \Delta t + E\_{r} \frac{1}{r} \frac{\partial}{\partial \theta} \frac{1 + \frac{r}{R} \cos \theta}{o \epsilon\_{\mathrm{c0}0}} \Delta t \} \tag{33}$$

where \* *<sup>e</sup> <sup>n</sup>* is calculated in Eq.(27). To calculate *<sup>n</sup>* <sup>1</sup> *i f* and *<sup>n</sup>* <sup>1</sup> *<sup>e</sup> <sup>n</sup>* , we then repeat **Step 1** for the solution of Eqs. (8,9) for *t* / 2 using \*\* *i f* and \*\* *<sup>e</sup> n* calculated from Eq.(30) and Eq.(33). The electric field is then updated to calculate <sup>1</sup> (, ) *<sup>n</sup> <sup>r</sup> E r* and <sup>1</sup> (, ) *<sup>n</sup> E r* by solving Poisson equation.

#### **2.3. The solution of Poisson's equation**

Equation (5) is solved by first Fourier transforming the equation in the periodic direction :

$$\phi(r,\theta) = \sum\_{m} \phi\_m(r)e^{im\theta} \tag{34}$$

and the resulting equation is then discretized in the radial direction over a uniform grid following the method of Knorr *et al*., 1980 (also rediscussed more recently in Crouseilles *et al.* 2009):

$$-A\_m \phi\_{m,j+1} + B\_m \phi\_{m,j} - C\_m \phi\_{m,j-1} = \frac{\{\Lambda r\}^2}{12} (\rho\_{m,j+1} + 10 \,\rho\_{m,j} + \rho\_{m,j-1}) \tag{35}$$

where

86 Numerical Simulation – From Theory to Industry

 <sup>2</sup> 1 1 (, ) <sup>4</sup> *<sup>r</sup> d j rj t V r*

The quantities 1 1 (, ) *<sup>r</sup> rj r j E r*

as indicated in Eq.(15).

The electric field *E*

our calculation.

**2.2. The solution for Step 2** 

of Eq. (29) yields the following solution to Eq.(28):

and are given, to an order <sup>2</sup> *O t* , by the expressions:

2

calculate:

 

 

 \* 2 2 ( , , / 2) ( 2 , 2 , ) *r r ej j ej r j nr t t nr <sup>t</sup>* 

We go to **Step2** and solve for *t* Eqs.(10,11). We first solve the equation:

remained small, since <sup>1</sup> *<sup>E</sup>*

2 *r ci*

  2

2

4

 *r*

 and 1 1 (, ) *<sup>r</sup> j rj E r*

interpolation, similar to what is described in Eq.(15). Finally to advance Eq.(9) by / 2 *t* , we

where the right hand side of Eq.(27) is calculated again using a cubic *B*-spline interpolation

0 *<sup>i</sup> <sup>i</sup> r i r ci ci r r <sup>f</sup> <sup>f</sup> vv f <sup>E</sup> E v t r v r v*

 

; *r r r ci ci r dv v dv v v E v E v dt r dt r*

Following the same iterative steps as described for Eq.(27), to an order <sup>2</sup> *O t* the solution

 , ,, , , , 2, 2, *i r ir f r v v t t f r v av bt*

Again the 2D interpolation in Eq. (30) is done using a tensor product of cubic *B*-spline, as explained in Eq.(15), and *a* and *b* are calculated with similar iterative steps as in Eqs.(23-26)

*<sup>t</sup> <sup>b</sup> a E vb*

*ci r <sup>t</sup> b a b E va*

 

 

*r*

   

> 

1 2

1 1

*r*

<sup>1</sup> ( )

<sup>1</sup> ( )( )

*r* 

*r*

 

(30)

*r*

 

*r*

1 1 1 1

1 cos( ) (, )

(27)

 

, and *r* ~ *R* =5000 towards the edge in

 

*j rj r j ci j r*

*r*

*r*

*j r*

0

are calculated by a cubic *B*-spline

(26)

(28)

(29)

(31)

(32)

1

 

 

*E r <sup>t</sup> <sup>R</sup>*

$$\mathbf{A}\_{m,j} = \mathbf{1} + \frac{\Delta r}{2r\_j} - \frac{m^2 \langle \Delta r \rangle^2}{12r\_j^2};\\\mathbf{B}\_{m,j} = \mathbf{2} + \frac{10m^2 \langle \Delta r \rangle^2}{12r\_j^2};\\\mathbf{C}\_{m,j} = \mathbf{1} - \frac{\Delta r}{2r\_j} - \frac{m^2 \langle \Delta r \rangle^2}{12r\_j^2} \tag{36}$$

and , ( ) *mj m j r* ; *i e n n* ; ( , , , ) *ii r r n fr dd* 

Equation (35) is solved using a tridiagonal algorithm. To get the boundary conditions, we assume in the present calculation that the deuterons and electrons currents hitting the cylindrical wall surface at *r R* are collected by a floating potential cylindrical vessel. Therefore we have the relations for the charge collected:

$$\left. \frac{\partial E\_r(r, \theta)}{\partial t} \right|\_{r=R} = -(J\_{ri}|\_{r=R} - J\_{re}|\_{r=R}) \qquad \text{or} \qquad E\_r(r, \theta)|\_{r=R} = -\left. \frac{t}{\imath} (J\_{ri}|\_{r=R} - J\_{re}|\_{r=R}) dt \right. \tag{37}$$

where

$$J\_{ri}(r, \theta) = \left[\upsilon\_r f\_i\left(r, \theta\_r \upsilon\_r, \upsilon\_\theta, t\right) d\upsilon\_\theta d\upsilon\_r \quad ; \qquad J\_{re} = \upsilon\_e V\_{dr} \tag{38}$$

and *Vdr* is defined in Eq.(17). *r r rR R E E* from Eq.(37) is used to obtain for the potential the following boundary condition at *r R* :

$$E\_r \Big|\_{R} = -\frac{\partial \phi}{\partial r} \Big|\_{R} \tag{39}$$

Charge Separation and Electric Field at a Cylindrical Plasma Edge 89

 *=* 

, and we note that the curves for 0<*R-r*<75 have reached a

*L*

the values of -0.1475 and -0.1399,

*<sup>L</sup> <sup>r</sup> R r <sup>E</sup>* was negligible by an order of magnitude

 *= .*

*/2* in Fig.(1). Indeed the code calculate for

 *=0* and

(43)

*ci* (normalized to

shows a small steady state oscillation

*pi* ),

*/2* (broken curve) and

 *= 0* is higher (in

 *= 0, /2*

calculated remained

*/2* as we

. This means that the code has executed more than 5

 *= 0* (full curve),

The code was executed for a sufficiently long time to reach a steady state. We noted around

 (dashed-dotted curve), at time *t* = 495 (left figure) and *t*=500 (right figure). At the edge and along the gradient the electric field *<sup>r</sup> E* is directed towards the center to the interior of the plasma. The electric field at the floating vessel wall at *r = R* was calculated using Eq. (37)

circle. We see from Fig. (1) that the electric field at

 *=* 

are respectively -0.1455, -0.1389 and -0.063 at *t*=500. For the solution of Poisson

*<sup>E</sup> dr R* 

 *,*the charge appearing in the system is not exactly equal to the charge collected

steady state stable equilibrium, with the charge collected at *r=R* remaining constant, while

equation, we set the value of the electric field *<sup>r</sup> r R <sup>E</sup>* at the cylindrical vessel wall exactly equal to the charge collected on the wall according to Eq.(37). If the charge collected on the

*r r R R r*

*L*

In our code, we allow at *<sup>L</sup> r Rr* for the possibility of plasma to flow across the boundary.

at the cylindrical vessel wall, but shows a small oscillation due the fact that there is a small plasma circulation at *<sup>L</sup> r Rr* . In this case the electric field *<sup>L</sup> <sup>r</sup> rRr <sup>E</sup>* is calculated from Eq.(41). At *t*=500, the value calculated for the right hand side of Eq.(41) is <sup>2</sup> 1.144x10 , while the value calculated by the code is <sup>2</sup> 1.155x10 (which is the value we see in the right figure

*/2* (see Fig.(6) below), and small around

curve ). The agreement is very good.

*R*

*R r*

*/2* respectively).

*<sup>E</sup> <sup>E</sup> dr R E* 

 *=0* and

around zero. The charge collected on the floating wall of the cylindrical vessel at

wall is equal to the charge appearing in the system, then according to Eq.(41) :

( )) / ; 0

the charge appearing in the system ( )) /

*R r*

very close to the values of -0.1455 and -0.1389 calculated for *<sup>r</sup> <sup>R</sup> <sup>E</sup>* at

Note that for our present set of the parameters the electric field *E*

*R*

*L*

*t*=500 that the code has indeed reach a steady state. Since 0 0.1 / 2

 *=*

**3. Results** 

 *=* 

and 

then the gyroperiod 0 2 / 88 *ci* 

at 128 points over the 2

absolute value) than the one at

gyroperiods. Figure(1) shows the electric field at

towards the center for 75<*R-r*<175, the curve at

which is what we see for the curves at

previously mentioned. And the quantity

 *=0* and

( <sup>3</sup> 2.x10 and <sup>3</sup> 1.x10 at

 *=0* to 

 *=* 

negligible from

 *=* 

in Fig.(1) for the

Around

By integrating the relation .*E* , we get a relation between the charge collected on the cylindrical wall and the charge appearing in the initially neutral plasma:

$$\left. \mathcal{R} \boldsymbol{E}\_{r} \right|\_{\mathbb{R}} - \left( \boldsymbol{R} - \boldsymbol{r}\_{\mathbb{L}} \right) \boldsymbol{E}\_{r} \Big|\_{\mathbb{R} - \boldsymbol{r}\_{\mathbb{L}}} + \int\_{\mathbb{R} - \boldsymbol{r}\_{\mathbb{L}}}^{\mathbb{R}} \frac{\partial \boldsymbol{E}\_{\theta}}{\partial \theta} d\boldsymbol{r} = \boldsymbol{\sigma}; \text{ where } \boldsymbol{\sigma} = \int\_{\mathbb{R} - \boldsymbol{r}\_{\mathbb{L}}}^{\mathbb{R}} \rho \boldsymbol{r} d\boldsymbol{r} \tag{40}$$

where *Lr* is the width of the plasma slab, which extends from *R rL* to *R* at the edge of the cylindrical plasma of radius *R* . is the charge appearing in the system. From Eq.(40), we get the following relation:

$$\left.E\_{r}\right|\_{R-r\_{L}} = -\frac{\left.\partial\phi\right|}{\left.\partial r\right|\_{R-r\_{L}}}\Big|\_{R-r\_{L}} = \left(\left.RE\_{r}\right|\_{R} - \left(\sigma-\int\_{R-r\_{L}}^{R} \frac{\partial E\_{\theta}}{\partial\theta} dr\right)\right) / \left(R-r\_{L}\right) \tag{41}$$

We assume that the plasma particles are allowed to enter or leave at the boundary at . *<sup>L</sup> r Rr* So the difference between the electric fields at the boundary *<sup>L</sup> r Rr* and at *r* = *R* must be such that Eq. (40) must be satisfied at every time-step in every direction . In the present simulation, *<sup>r</sup> r R <sup>E</sup>* is calculated from Eq. (37), which defines the derivative of the potential at the right boundary to be used in the solution of Eq. (35). We fix the potential to be zero at the boundary at *<sup>L</sup> r Rr* and solve Poisson equation in Eq.(35) for the potential. Then the resulting electric field at *<sup>L</sup> <sup>r</sup> rRr <sup>E</sup>* , calculated at *<sup>L</sup> r Rr* from Eq. (5), must satisfy Eq. (41).

The initial density profiles at the neutral plasma edge are given by:

$$n\_i(r) = n\_c(r) = 0.5(1 + \tanh\left(\left(R - r - 2r\_L / 5\right) / 4\right))\tag{42}$$

The profiles in Eq. (42) situate the steep gradient to be centered at a distance of 2 / 5 *Lr* from the wall of the vessel, which put the plasma relatively close to the floating wall of the vessel. The wall of the vessel will collect the charge coming from the plasma, especially due to the large ions gyroradius. The system is solved for an edge thickness 175 *Lr* , and *R*=5000 for the radius of the cylinder (we expect this radius to be large in a tokamak, so in the domain ( ,) *<sup>L</sup> r RrR* , the variation of the quantity *r R*/ remains very close to 1.

We use *N*=250 grid points in space in the radial direction, and 128 grid points in the azimuthal direction. 80 grid points are used in each velocity direction. The velocity extrema for the ions velocities are 4 / *e i T T* in our normalized units, with *Te/Ti* =1 in the present simulations.

### **3. Results**

88 Numerical Simulation – From Theory to Industry

following boundary condition at *r R* :

By integrating the relation .*E*

cylindrical plasma of radius *R* .

get the following relation:

resulting electric field at

simulations.

*<sup>r</sup> <sup>R</sup>*

( , ) , , , , ; *ri ri r <sup>r</sup> re e dr J r v f r v v t d d J nV*

( ) ; where

*r* 

*r Lr R R r*

cylindrical wall and the charge appearing in the initially neutral plasma:

*L*

( ( )) / ( ) *<sup>L</sup>*

such that Eq. (40) must be satisfied at every time-step in every direction

The initial density profiles at the neutral plasma edge are given by:

( ,) *<sup>L</sup> r RrR* , the variation of the quantity *r R*/ remains very close to 1.

 

*E*

 

*R*

*r* 

(38)

, we get a relation between the charge collected on the

 

is the charge appearing in the system. From Eq.(40), we

(40)

. In the present

*L L*

*R R*

*R r R r*

*R*

(41)

*<sup>L</sup> <sup>r</sup> rRr <sup>E</sup>* , calculated at *<sup>L</sup> r Rr* from Eq. (5), must satisfy Eq. (41).

(39)

 

and *Vdr* is defined in Eq.(17). *r r rR R E E* from Eq.(37) is used to obtain for the potential the

*<sup>E</sup> RE R r E dr rdr* 

where *Lr* is the width of the plasma slab, which extends from *R rL* to *R* at the edge of the

*L L*

We assume that the plasma particles are allowed to enter or leave at the boundary at . *<sup>L</sup> r Rr* So the difference between the electric fields at the boundary *<sup>L</sup> r Rr* and at *r* = *R* must be

simulation, *<sup>r</sup> r R <sup>E</sup>* is calculated from Eq. (37), which defines the derivative of the potential at the right boundary to be used in the solution of Eq. (35). We fix the potential to be zero at the boundary at *<sup>L</sup> r Rr* and solve Poisson equation in Eq.(35) for the potential. Then the

( ) ( ) 0.5(1 tanh (( 2 / 5) / 4)) *i e <sup>L</sup> nr nr Rr r* (42)

The profiles in Eq. (42) situate the steep gradient to be centered at a distance of 2 / 5 *Lr* from the wall of the vessel, which put the plasma relatively close to the floating wall of the vessel. The wall of the vessel will collect the charge coming from the plasma, especially due to the large ions gyroradius. The system is solved for an edge thickness 175 *Lr* , and *R*=5000 for the radius of the cylinder (we expect this radius to be large in a tokamak, so in the domain

We use *N*=250 grid points in space in the radial direction, and 128 grid points in the azimuthal direction. 80 grid points are used in each velocity direction. The velocity extrema for the ions velocities are 4 / *e i T T* in our normalized units, with *Te/Ti* =1 in the present

*rr L R r R R r R r <sup>E</sup> <sup>E</sup> RE dr R r*

 

The code was executed for a sufficiently long time to reach a steady state. We noted around *t*=500 that the code has indeed reach a steady state. Since 0 0.1 / 2 *ci* (normalized to *pi* ), then the gyroperiod 0 2 / 88 *ci* . This means that the code has executed more than 5 gyroperiods. Figure(1) shows the electric field at  *= 0* (full curve),  *= /2* (broken curve) and  *=*  (dashed-dotted curve), at time *t* = 495 (left figure) and *t*=500 (right figure). At the edge and along the gradient the electric field *<sup>r</sup> E* is directed towards the center to the interior of the plasma. The electric field at the floating vessel wall at *r = R* was calculated using Eq. (37) at 128 points over the 2 circle. We see from Fig. (1) that the electric field at  *= 0* is higher (in absolute value) than the one at  *=*, and we note that the curves for 0<*R-r*<75 have reached a steady state stable equilibrium, with the charge collected at *r=R* remaining constant, while towards the center for 75<*R-r*<175, the curve at  *=*  shows a small steady state oscillation around zero. The charge collected on the floating wall of the cylindrical vessel at  *= 0, /2* and are respectively -0.1455, -0.1389 and -0.063 at *t*=500. For the solution of Poisson equation, we set the value of the electric field *<sup>r</sup> r R <sup>E</sup>* at the cylindrical vessel wall exactly equal to the charge collected on the wall according to Eq.(37). If the charge collected on the wall is equal to the charge appearing in the system, then according to Eq.(41) :

$$\left.E\_r\right|\_{\mathbb{R}} \approx (\sigma - \int\_{\mathbb{R}-\eta\_{\mathbb{L}}}^{\mathbb{R}} \frac{\partial E\_\theta}{\partial \theta} dr\,\rangle) / \,\mathbb{R} \quad ; \qquad \left.E\_r\right|\_{\mathbb{R}-\eta\_{\mathbb{L}}} \approx 0 \tag{43}$$

which is what we see for the curves at  *=0* and */2* in Fig.(1). Indeed the code calculate for the charge appearing in the system ( )) / *L R R r <sup>E</sup> dr R* the values of -0.1475 and -0.1399, very close to the values of -0.1455 and -0.1389 calculated for *<sup>r</sup> <sup>R</sup> <sup>E</sup>* at  *=0* and */2* as we previously mentioned. And the quantity *<sup>L</sup> <sup>r</sup> R r <sup>E</sup>* was negligible by an order of magnitude ( <sup>3</sup> 2.x10 and <sup>3</sup> 1.x10 at  *=0* and */2* respectively).

In our code, we allow at *<sup>L</sup> r Rr* for the possibility of plasma to flow across the boundary. Note that for our present set of the parameters the electric field *E* calculated remained negligible from  *=0* to */2* (see Fig.(6) below), and small around  *= .*

Around  *= ,*the charge appearing in the system is not exactly equal to the charge collected at the cylindrical vessel wall, but shows a small oscillation due the fact that there is a small plasma circulation at *<sup>L</sup> r Rr* . In this case the electric field *<sup>L</sup> <sup>r</sup> rRr <sup>E</sup>* is calculated from Eq.(41). At *t*=500, the value calculated for the right hand side of Eq.(41) is <sup>2</sup> 1.144x10 , while the value calculated by the code is <sup>2</sup> 1.155x10 (which is the value we see in the right figure in Fig.(1) for the  *=* curve ). The agreement is very good.

Charge Separation and Electric Field at a Cylindrical Plasma Edge 91

**Figure 2.** Potential at

**Figure 3.** Plot at

, , *ir i nE P* at 0

(0.1 / 2) /(1. 0.2( / )cos ) *v rR* 

*t*=495 (left gigure) and *t*=500 (right figure).

0 (full curve),

 

0 of the electric field *Er* (full curve), the Lorentz force

 

The / *i i P n* term (broken curve) shows a very good agreement along the gradient with the solid curve for *<sup>r</sup> E* , and the Lorentz force appears negligible along the gradient. In a region of about two gyroradii from the wall (around 40 Debye lengths from the wall), we have small irregular oscillations in space (and time), the accuracy of / *i i P n* being degraded by the division with a very small value of the density *<sup>i</sup> n* appearing close to the vessel surface. To avoid this problem, we plot in Fig. (4) the quantities (0.1 / 2) /(1. 0.2( / )cos ) *<sup>i</sup> n v rR*

curve). The curve -*ni*/2 is plotted for reference ( dashed- 3 dotted curve). At time *t* =500.

( note that *i i J nv*

/ 2 (broken curve) and

  (dashed-dotted curve) at

). We see that there is a very nice agreement for

,

(dashed-dotted. curve), and the pressure force / *i i P n* (broken

**Figure 1.** Electric field *Er* at 0 (full curve), / 2 (broken curve), and (dashed-dotted curve), at *t*=495 (left figure) and *t*=500 (right figure).

Figure(2) shows the potential at  *= 0* (full curve),  *= /2* (broken line) and  *=*  (dasheddotted line), at time *t* = 495 (left figue) and *t*=500 (right figure). The curves at  *= 0* and  *= /2*  show negligible oscillation, while the curve at  *=*  shows a small oscillation. Since as indicated in Fig.(1) the *Er* profiles are essentially constant, especially for 0<*R-r*<75 along the gradient, then the oscillation of the curves at  *=*  (dashed-dotted line) is taking place in such a way that the slope of the potential *<sup>r</sup> E r* remains constant for 0<*R-r*<75. Only for 75<*R-r*<175 does the small oscillation of the curve at  *=* in Fig.(2) (dashed-dotted line)

translates into a small oscillation of the electric field *Er.*.

Figure(3) gives at time *t* = 500 and for  *= 0* the electric field *Er* (full curve). Also at  *= 0* the dashed-dotted curve in Fig.(2) gives the Lorentz force, which in our normalized units is given by 0.1 / /(1. 0.2 cos ) 2 *ci pi r v v R* , and the broken curve gives the pressure force / *i i P n* , 0.5 *i i ir i P nT T*, with the following definition:

$$T\_{\rm ir,\theta}(r,\theta) = \frac{1}{n\_i} \int d\upsilon\_r \, d\upsilon\_\theta \, \left(\upsilon\_{r,\theta} - \varsigma \upsilon\_{r,\theta} > \right)^2 f\_i\left(r,\theta,\upsilon\_{r'},\upsilon\_\theta\right) \tag{44}$$

$$<\boldsymbol{v}\_{r,\theta} > = \frac{1}{n\_i} \Big[d\boldsymbol{v}\_r \, d\boldsymbol{v}\_\theta \, \boldsymbol{v}\_{r,\theta} \, \boldsymbol{f}\_i \Big(\boldsymbol{r}\_r \boldsymbol{\theta}, \boldsymbol{v}\_{r'} \, \boldsymbol{v}\_\theta\Big) \; ; \qquad \qquad \boldsymbol{n}\_i(\boldsymbol{r}, \boldsymbol{\theta}) = \Big[d\boldsymbol{v}\_r \, d\boldsymbol{v}\_\theta \, \boldsymbol{f}\_i \Big(\boldsymbol{r}\_r \boldsymbol{\theta}, \boldsymbol{v}\_r \, \boldsymbol{v}\_\theta\Big) \; \tag{45}$$

**Figure 1.** Electric field *Er* at

Figure(2) shows the potential at

show negligible oscillation, while the curve at

gradient, then the oscillation of the curves at

Figure(3) gives at time *t* = 500 and for

*i*

*n*

pressure force / *i i P n* , 0.5 *i i ir i P nT T*

such a way that the slope of the potential *<sup>r</sup> E*

75<*R-r*<175 does the small oscillation of the curve at

translates into a small oscillation of the electric field *Er.*.

given by 0.1 / /(1. 0.2 cos ) 2 *ci pi*

 

*v v*

curve), at *t*=495 (left figure) and *t*=500 (right figure).

0 (full curve),

dotted line), at time *t* = 495 (left figue) and *t*=500 (right figure). The curves at

 

 <sup>2</sup> , , ,

*i*

*n*

 

 *= 0* (full curve),

indicated in Fig.(1) the *Er* profiles are essentially constant, especially for 0<*R-r*<75 along the

*r* 

dashed-dotted curve in Fig.(2) gives the Lorentz force, which in our normalized units is

<sup>1</sup> (, ) ( ) ,, , *ir rrr i r*

 , , <sup>1</sup> , , , ; ( , ) ,, , *<sup>r</sup> r ri r i ri r*

 

*v dv dv v f r v v n r dv dv f r v v*

*T r dv dv v v f r v v*

*r*

*R*

, with the following definition:

 

 

 *=* 

/ 2 (broken curve), and

 *=* 

 *= 0* the electric field *Er* (full curve). Also at

 (44)

(45)

 *=* 

 *=*   

*/2* (broken line) and

(dashed-dotted

 *=* 

 *= 0* and

in Fig.(2) (dashed-dotted line)

shows a small oscillation. Since as

(dashed-dotted line) is taking place in

remains constant for 0<*R-r*<75. Only for

, and the broken curve gives the

(dashed-

 *= /2* 

 *= 0* the

 

**Figure 2.** Potential at 0 (full curve), / 2 (broken curve) and (dashed-dotted curve) at *t*=495 (left gigure) and *t*=500 (right figure).

**Figure 3.** Plot at 0 of the electric field *Er* (full curve), the Lorentz force (0.1 / 2) /(1. 0.2( / )cos ) *v rR* (dashed-dotted. curve), and the pressure force / *i i P n* (broken curve). The curve -*ni*/2 is plotted for reference ( dashed- 3 dotted curve). At time *t* =500.

The / *i i P n* term (broken curve) shows a very good agreement along the gradient with the solid curve for *<sup>r</sup> E* , and the Lorentz force appears negligible along the gradient. In a region of about two gyroradii from the wall (around 40 Debye lengths from the wall), we have small irregular oscillations in space (and time), the accuracy of / *i i P n* being degraded by the division with a very small value of the density *<sup>i</sup> n* appearing close to the vessel surface. To avoid this problem, we plot in Fig. (4) the quantities (0.1 / 2) /(1. 0.2( / )cos ) *<sup>i</sup> n v rR* , , , *ir i nE P* at 0 ( note that *i i J nv* ). We see that there is a very nice agreement for

the relation *ir i nE P* along the gradient (the density / 10 *<sup>i</sup> n* is also plotted in Fig.(4) to locate the different profiles with respect to the gradient). The Lorentz force term remains negligible in the gradient region, and remains very small in the bulk at 75<*R-r*<175. The / *i i P n* term is zero in the bulk since *ni* and *Ti* are essentially flat in that region. We have Charge Separation and Electric Field at a Cylindrical Plasma Edge 93

which, as we mentioned

 *=* (dashed-

is interesting, it shows

 *= 0* (full curve, essentially equal to zero),

  (dash-dot curve),

(dashed-dotted

at

 

**Figure 5.** Charge density at

curve), at *t*=500.

This results in a

a negligible oscillation at

the oscillating curve for *E*

**Figure 6.** Electric field *E*

 at 

at *t*=485 (left figure) and *t*=495 (right figure).

0 (full curve),

 

variation in the inner region.

*t*=485 and *t*=495. It shows no oscillation *E*

 = 

dotted curve), on the high field side of the cylinder. The curve at

 *=*  0 (full curve),

 

variation with an electric field <sup>1</sup> *<sup>E</sup>*

before, is small for the parameters we are using. Figures (6) presents this electric field *E*

 in at 

We have noted in Fig.(1) the constant value of the *<sup>r</sup> E* profiles, especially for 0<*R-r*<75 along the gradient, translates the constant slope of the potential. However, the potential profile at

, although keeping essentially the same slope for 0<*R-r*<75, shows an oscillation in time.

/ 2 (broken curve) and

/2 (broken curve), and a small oscillation for the

*r*

having a constant flat value in the gradient region, and a linear

/ 2 (broken curve), and

 *=* 

seen in this case that the total charge ( )) / *L R R r <sup>E</sup> dr R* appearing in the system and

calculated by the code by integrating the charge as in Eq. (40) and by calculating *E* from Poisson equation, is essentially equal to the electric field *<sup>r</sup> r R <sup>E</sup>* . We note that curves similar to Fig.(4) can be calculated at different angles , showing *i r n E* essentially balanced by the *<sup>i</sup> P* .

**Figure 4.** Plot at 0 of *i r n E* (solid curve), (0.1 / 2) /(1. 0.2( / )cos ) *<sup>i</sup> n v rR* (dash-dot curve), and *<sup>i</sup> P* (broken-curve). The curve –*ni*/10 is plotted for reference (dashed-3 dotted curve). At time *t*=500.

The initial density profiles are given in Eq.(42), so the initial charge is zero. Fig.(5) shows the charge density *i e n n* at the time *t* = 500 and 0 (full curve), / 2 (broken curve) and (dashed-dotted curve). It is interesting to note that the code was able to maintain the steep density gradients (see Fig.(4) for a plot of *<sup>i</sup> n* ). The stable in time steep density profile changes in space rapidly along the gradient over an ion orbit size ( / 20 *i De* ), and the relaxation of the steep gradients during the simulation determines the charge density *i e n n* . The charge density is important along the gradient at the plasma edge. The electrons, described by the guiding center equation given in Eq.(3) cannot compensate along the gradient the charge separation caused by the finite ion gyroradius, which results in the charge density we see in Fig.(5).

92 Numerical Simulation – From Theory to Industry

seen in this case that the total charge ( )) /

to Fig.(4) can be calculated at different angles

*<sup>i</sup> P* .

**Figure 4.** Plot at

time *t*=500.

  charge density we see in Fig.(5).

charge density *i e n n* at the time *t* = 500 and 0

the relation *ir i nE P* along the gradient (the density / 10 *<sup>i</sup> n* is also plotted in Fig.(4) to locate the different profiles with respect to the gradient). The Lorentz force term remains negligible in the gradient region, and remains very small in the bulk at 75<*R-r*<175. The / *i i P n* term is zero in the bulk since *ni* and *Ti* are essentially flat in that region. We have

*L*

*<sup>E</sup> dr R* 

appearing in the system and

, showing *i r n E* essentially balanced by the

  / 2 (broken curve) and

), and the

(dash-dot

from

*R*

*R r*

Poisson equation, is essentially equal to the electric field *<sup>r</sup> r R <sup>E</sup>* . We note that curves similar

0 of *i r n E* (solid curve), (0.1 / 2) /(1. 0.2( / )cos ) *<sup>i</sup> n v rR*

curve), and *<sup>i</sup> P* (broken-curve). The curve –*ni*/10 is plotted for reference (dashed-3 dotted curve). At

The initial density profiles are given in Eq.(42), so the initial charge is zero. Fig.(5) shows the

relaxation of the steep gradients during the simulation determines the charge density *i e n n* . The charge density is important along the gradient at the plasma edge. The electrons, described by the guiding center equation given in Eq.(3) cannot compensate along the gradient the charge separation caused by the finite ion gyroradius, which results in the

changes in space rapidly along the gradient over an ion orbit size ( / 20 *i De*

 (dashed-dotted curve). It is interesting to note that the code was able to maintain the steep density gradients (see Fig.(4) for a plot of *<sup>i</sup> n* ). The stable in time steep density profile

(full curve),

 

calculated by the code by integrating the charge as in Eq. (40) and by calculating *E*

**Figure 5.** Charge density at 0 (full curve), / 2 (broken curve) and (dashed-dotted curve), at *t*=500.

We have noted in Fig.(1) the constant value of the *<sup>r</sup> E* profiles, especially for 0<*R-r*<75 along the gradient, translates the constant slope of the potential. However, the potential profile at  *=*  , although keeping essentially the same slope for 0<*R-r*<75, shows an oscillation in time. This results in a variation with an electric field <sup>1</sup> *<sup>E</sup> r* which, as we mentioned before, is small for the parameters we are using. Figures (6) presents this electric field *E* at *t*=485 and *t*=495. It shows no oscillation *E* in at  *= 0* (full curve, essentially equal to zero), a negligible oscillation at = /2 (broken curve), and a small oscillation for the  *=*  (dasheddotted curve), on the high field side of the cylinder. The curve at  *=*  is interesting, it shows the oscillating curve for *E* having a constant flat value in the gradient region, and a linear variation in the inner region.

**Figure 6.** Electric field *E* at 0 (full curve), / 2 (broken curve), and (dash-dot curve), at *t*=485 (left figure) and *t*=495 (right figure).

Figure (7) gives at *t* = 500 the temperature *Tir* as defined in Eq. (43), for  *= 0* (full curve), = /2 (broken curve) and  *=*  (dashed-dotted curve). The dashed-3 dotted curve is for *niTir*, which is essentially the same for all angles . Figure (8) presents similar results for *Ti* as defined in Eq.(43). Note in Fig.(7) and Fig.(8) the division by the very small density at the edge gives irregular oscillations. However the dashed-3 dotted curve in Fig.(7) and Fig.(8) is very smooth, and is for *niTir* and *niTi* respectively, which removes the problem of the division by the very small density at the edge. We present in Fig.(9) the total pressure 0.5 *i i ir i P nT T* for  *= 0* (full curve), = /2 (broken curve) and  *=*  (dashed-dotted curve). We see that these curves are essentially identical, for the parameters we are using in the present simulation there is no variation in for the total pressure term.

Charge Separation and Electric Field at a Cylindrical Plasma Edge 95

/ 2 (broken curve) and

 = 

, where the diamagnetic drift

component of the <sup>2</sup> *E BB* / drift, which in our normalized

at the edge, and above the acoustic speed at

(46)

 ).

field (see Fig.6), the <sup>2</sup> *E BB* / drift has a small oscillating

in the radial direction, with a curve similar to what is

 *= 0* (full curve),

 

 *=* 

,

/2 (broken curve)

shows a

 *= 0*

units is written 1 0.2( / ) cos 2 / 0.1 *<sup>r</sup> E rR*

result we get if we calculate the poloidal current:

We get indeed a negligible value for *<sup>i</sup> J*

We note that due to the small *E*

component (1 0.2( / )cos ) *E rR* 

at

 (, ) , , , , *i ir r J r v f r v v t dv dv* 

discussed in Fig.(1) for *Er* , while towards the center for 75<*R-r*<175 the curves at

 *=* 

0 (full curve),

, at

small steady state oscillation around zero, as previously mentioned for *Er* in the right figure in Fig.(1). Note in Fig. (10) that this azimuthal (*i.e.* poloidal) <sup>2</sup> *E BB* / drift is flat close to the

 

(dashed-dotted curve). The curves for 0<*R-r*<75 have reached a stable equilibrium as

/2 at the edge (velocities are normalized to the acoustic speed *Cs*). We plot in Fig. (11)

direction) <sup>2</sup> / *i D n E BB v*

is <sup>2</sup> / *D ii i v B n T n eB* . We see that the total current is essentially zero, *i.e.* the profile is

adjusting itself so that the <sup>2</sup> *E BB* / drift and the diamagnetic drift are essentially equal and

opposite (in our units, the total poloidal current is 1 0.2( / ) cos 2 / 0.1 *ir i nE p r R*

 

essentially equal to zero , showing a small oscillation around zero in Fig.(11), at the left boundary). So the <sup>2</sup> *E BB* / drift is balanced fairly well by the diamagnetic drift. This is the

> 

 . (note that *i i J nv* 

 

**Figure 9.** Pressure 0.5 *i i ir i P nT T*

Figure(10) presents the plot of the

edge, and is below the acoustic speed at

the total poloidal current (in the

(dashed-dotted curve), at *t*= 500.

and  *=* 

and = 

**Figure 7.** Temperature *Tir* at 0 (full curve), / 2 (broken curve) and (dashed-dotted curve). The dashed -3 dotted curve is for *niTir*, which is essentially the same for all . Time *t*=500.

**Figure 8.** Temperature *<sup>i</sup> T* at 0 (full curve), / 2 (broken curve) and (dash-dot curve). The dashed-3 dotted curve is for *i i n T* , which is essentially the same for all . Time *t*=500.

very smooth, and is for *niTir* and *niTi*

for

/2 (broken curve) and

0.5 *i i ir i P nT T*

**Figure 7.** Temperature *Tir* at

**Figure 8.** Temperature *<sup>i</sup> T*

The dashed-3 dotted curve is for *i i n T*

 at 0 

 = 

Figure (7) gives at *t* = 500 the temperature *Tir* as defined in Eq. (43), for

 *= 0* (full curve),

 *=* 

*niTir*, which is essentially the same for all angles

0 (full curve),

curve). The dashed -3 dotted curve is for *niTir*, which is essentially the same for all

(full curve),

 

, which is essentially the same for all

 

/ 2 (broken curve) and

/ 2 (broken curve) and

the present simulation there is no variation in

. Figure (8) presents similar results for *Ti*

 *=* 

 

. Time *t*=500.

(dashed-dotted

(dash-dot curve).

. Time *t*=500.

 

respectively, which removes the problem of the

(dashed-dotted curve). The dashed-3 dotted curve is for

/2 (broken curve) and

for the total pressure term.

as defined in Eq.(43). Note in Fig.(7) and Fig.(8) the division by the very small density at the edge gives irregular oscillations. However the dashed-3 dotted curve in Fig.(7) and Fig.(8) is

division by the very small density at the edge. We present in Fig.(9) the total pressure

curve). We see that these curves are essentially identical, for the parameters we are using in

 =   *= 0* (full curve),

(dashed-dotted

**Figure 9.** Pressure 0.5 *i i ir i P nT T* at 0 (full curve), / 2 (broken curve) and (dashed-dotted curve), at *t*= 500.

Figure(10) presents the plot of the component of the <sup>2</sup> *E BB* / drift, which in our normalized units is written 1 0.2( / ) cos 2 / 0.1 *<sup>r</sup> E rR* , at  *= 0* (full curve), = /2 (broken curve) and  *=*  (dashed-dotted curve). The curves for 0<*R-r*<75 have reached a stable equilibrium as discussed in Fig.(1) for *Er* , while towards the center for 75<*R-r*<175 the curves at  *=*  shows a small steady state oscillation around zero, as previously mentioned for *Er* in the right figure in Fig.(1). Note in Fig. (10) that this azimuthal (*i.e.* poloidal) <sup>2</sup> *E BB* / drift is flat close to the edge, and is below the acoustic speed at  *=*  at the edge, and above the acoustic speed at  *= 0* and = /2 at the edge (velocities are normalized to the acoustic speed *Cs*). We plot in Fig. (11) the total poloidal current (in the direction) <sup>2</sup> / *i D n E BB v* , where the diamagnetic drift is <sup>2</sup> / *D ii i v B n T n eB* . We see that the total current is essentially zero, *i.e.* the profile is adjusting itself so that the <sup>2</sup> *E BB* / drift and the diamagnetic drift are essentially equal and opposite (in our units, the total poloidal current is 1 0.2( / ) cos 2 / 0.1 *ir i nE p r R* , essentially equal to zero , showing a small oscillation around zero in Fig.(11), at the left boundary). So the <sup>2</sup> *E BB* / drift is balanced fairly well by the diamagnetic drift. This is the result we get if we calculate the poloidal current:

$$J\_{\theta i}(r,\theta) = \int \upsilon\_{\theta} f\_i \left( r, \theta, \upsilon\_r, \upsilon\_{\theta'}, t \right) d\upsilon\_{\theta} d\upsilon\_r \tag{46}$$

We get indeed a negligible value for *<sup>i</sup> J* . (note that *i i J nv* ).

We note that due to the small *E* field (see Fig.6), the <sup>2</sup> *E BB* / drift has a small oscillating component (1 0.2( / )cos ) *E rR* in the radial direction, with a curve similar to what is

presented in Fig.(6). This radial oscillation is negligible around 0 and / 2 , and is small on the high field side around , as previously discussed for *E*in Fig.(6).

**Figure 10.** Plot of 1 0.2( / ) cos 2 / 0.1 *<sup>r</sup> E rR* for 0 (full curve), / 2 (broken curve) and (dashed-dotted curve). At time *t*=500.

**Figure 11.** Plot of the total poloidal current 1 0.2 ( / )cos 2 / 0.1 *ir i nE P r R* for 0 (full curve), / 2 (broken curve) and (dashed-dotted curve). At time *t*=500.

Charge Separation and Electric Field at a Cylindrical Plasma Edge 97

on the

in

bound to the magnetic field cannot compensate along the gradient the charge separation due to the finite ions' gyroradius. In the cylindrical geometry considered, a fully kinetic equation has been used to describe the ions, and a guiding center equation has been used to describe the electrons bound to the magnetic field. These equations have been solved using the method of characteristics (Shoucri, 2008a,b,c,d, 2009a). The numerical method used for the solution in cylindrical geometry, based on an integration of the equations along the characteristics coupled to a two-dimensional interpolation (Shoucri *et al*., 2004) applied successively in configuration space and in velocity space, is producing accurate results.

The problem of the formation of a charge separation is of great importance in the study of the H-mode physics in tokamaks. We have considered the case where the gradient in the density profiles is located in front of a floating cylindrical vessel. So the charge appearing in the system is essentially equal to the charge collected on the walls of the floating vessel. The solution shows in the radial direction that the electric field along the gradient is balanced by the radial gradient of the pressure, and the total poloidal current is essentially zero. The present results where electrons are described by a guiding center equation are close to what has been previously reported when assuming the electrons stationary and frozen by the magnetic field, and the code shows that a solution with a steep gradient, maintaining a charge separation, where the electron and ion densities vary rapidly over an ion gyroradius, is possible. Also the calculation with the present set of parameters which allow a small value

to exist, shows the presence of a small oscillation of *E*

to the radial component of the <sup>2</sup> *E BB* / drift , equal to *E rR* 1 0.2 ( / )cos 2 / 0.1

The present code is a step closer to a code which will include "neoclassical" effects due to toroidal geometry, which can play a role in this problem, such as the neoclassical enhancement of the classical ion polarization drift, or the neoclassical damping of poloidal

The author is grateful to Dr. Rejean Girard for his support, and to the Centre de calcul scientifique de l'IREQ (CASIR) for computer time for the simulations presented in this work.

Abbott, B.A. (1966). *An Introduction to the Method of Characteristics*; Thames and Hudson,

(see Fig.(6)), to which is associated a small radial oscillation due

of the poloidal field *E*

high field side around

our normalized units.

**Author details** 

**Acknowledgement** 

Magdi Shoucri

**5. References** 

London.

 

flows (Stix, 1973, Hirshman 1978, Waltz *et. al.*, 1999).

*Institut de Recherche Hydro-Québec (IREQ), Varennes, Québec, Canada* 

#### **4. Conclusion**

We have presented in this work the self-consistent kinetic solution for the problem of the generation of a charge separation and an electric field at a plasma edge, under the combined effect of a large ratio of the ions' gyroradius to the Debye length / *i De* (equal to 20 in the simulation we have presented) and a steep density gradient, when the electrons which are bound to the magnetic field cannot compensate along the gradient the charge separation due to the finite ions' gyroradius. In the cylindrical geometry considered, a fully kinetic equation has been used to describe the ions, and a guiding center equation has been used to describe the electrons bound to the magnetic field. These equations have been solved using the method of characteristics (Shoucri, 2008a,b,c,d, 2009a). The numerical method used for the solution in cylindrical geometry, based on an integration of the equations along the characteristics coupled to a two-dimensional interpolation (Shoucri *et al*., 2004) applied successively in configuration space and in velocity space, is producing accurate results.

The problem of the formation of a charge separation is of great importance in the study of the H-mode physics in tokamaks. We have considered the case where the gradient in the density profiles is located in front of a floating cylindrical vessel. So the charge appearing in the system is essentially equal to the charge collected on the walls of the floating vessel. The solution shows in the radial direction that the electric field along the gradient is balanced by the radial gradient of the pressure, and the total poloidal current is essentially zero. The present results where electrons are described by a guiding center equation are close to what has been previously reported when assuming the electrons stationary and frozen by the magnetic field, and the code shows that a solution with a steep gradient, maintaining a charge separation, where the electron and ion densities vary rapidly over an ion gyroradius, is possible. Also the calculation with the present set of parameters which allow a small value of the poloidal field *E* to exist, shows the presence of a small oscillation of *E* on the high field side around (see Fig.(6)), to which is associated a small radial oscillation due to the radial component of the <sup>2</sup> *E BB* / drift , equal to *E rR* 1 0.2 ( / )cos 2 / 0.1 in our normalized units.

The present code is a step closer to a code which will include "neoclassical" effects due to toroidal geometry, which can play a role in this problem, such as the neoclassical enhancement of the classical ion polarization drift, or the neoclassical damping of poloidal flows (Stix, 1973, Hirshman 1978, Waltz *et. al.*, 1999).
