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, where is the ions’ gyroradius and 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 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
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
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
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
2. The relevant equations and the numerical method
We consider the cylindrical coordinate system (. The plasma is uniform in the
where from Eq.(1) the ion cyclotron frequency. The electrons are described by the normalized 2D guiding-center equation:
where the drift velocity. Equation (3) can be developed in our 2D system to give the following equation:
This system is coupled to Poisson’s equation:
Time is normalized to the inverse plasma frequency, velocity is normalized to the acoustic velocity, and length is normalized to the Debye length. T
In the present calculations the parameter. With the initial distribution function for the ions in Eq. (6), we have spatially constant, the factor 2T
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.
We solve for t/2 the equations:
We solve for t the equations:
Again Eq.(10) is solved by a method of characteristics, to be described below.
We repeat Step1 for t/2, and then solve Eq.(5) to calculate the new electric field,.
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:
at a given in velocity space. The solution of Eqs (12) originating at () at time t, and reaching the grid point () at can be written as follows, for a half time-step t/2:
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 and, we write for the interpolation function:
taking into account that is periodic. andare cubic splines.The cubic B-spline is defined as:
and otherwise. The grid size h (eitheror in our notation), is assumed uniform. For the calculation of the coefficients of the cubic B-spline interpolation function 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:
Equations (17) are solved using an iterative process. We assume that at the time-step, r is at the grid point and is at the grid point. The following leapfrog scheme can be written for the solution of Eq.(17):
where is the point where the characteristic is originating at (not necessarily a grid point). Put:
we can rewrite Eqs.(18,19) as follows:
which are implicit equations for, and which are solved by iteration. Usually, two iterations are sufficient for convergence. We start with. We get at the first iteration:
2.2. The solution for Step 2
We go to Step2 and solve for Eqs.(10,11). We first solve the equation:
The electric field remained small, since, and r ~ R =5000 towards the edge in our calculation.
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, by the expressions:
The density is calculated by integrating Eq.(11) as follows:
where is calculated in Eq.(27). To calculate and, we then repeat Step 1 for the solution of Eqs. (8,9) for using and calculated from Eq.(30) and Eq.(33). The electric field is then updated to calculate and 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:
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 atare collected by a floating potential cylindrical vessel. Therefore we have the relations for the charge collected:
By integrating the relation, we get a relation between the charge collected on the cylindrical wall and the charge appearing in the initially neutral plasma:
where is the width of the plasma slab, which extends from to at the edge of the cylindrical plasma of radius. is the charge appearing in the system. From Eq.(40), we get the following relation:
We assume that the plasma particles are allowed to enter or leave at the boundary at. So the difference between the electric fields at the boundary and at r = R must be such that Eq. (40) must be satisfied at every time-step in every direction. In the present simulation, 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 and solve Poisson equation in Eq.(35) for the potential. Then the resulting electric field at, calculated at from Eq. (5), must satisfy Eq. (41).
The initial density profiles at the neutral plasma edge are given by:
The profiles in Eq. (42) situate the steep gradient to be centered at a distance of 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, and R=5000 for the radius of the cylinder (we expect this radius to be large in a tokamak, so in the domain, the variation of the quantity 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 in our normalized units, with T
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 (normalized to), then the gyroperiod. 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 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 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) :
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 the values of -0.1475 and -0.1399, very close to the values of -0.1455 and -0.1389 calculated for at =0 and /2 as we previously mentioned. And the quantitywas negligible by an order of magnitude (and at =0 and /2 respectively).
In our code, we allow at for the possibility of plasma to flow across the boundary. Note that for our present set of the parameters the electric field 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. In this case the electric field is calculated from Eq.(41). At t=500, the value calculated for the right hand side of Eq.(41) is, while the value calculated by the code is (which is the value we see in the right figure in Fig.(1) for the = curve ). The agreement is very good.
Figure(2) shows the potential at = 0 (full curve), = /2 (broken line) and = (dashed-dotted 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 E
Figure(3) gives at time t = 500 and for = 0 the electric field E
The initial density profiles are given in Eq.(42), so the initial charge is zero. Fig.(5) shows the charge density at the time t = 500 and (full curve), (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). The stable in time steep density profile changes in space rapidly along the gradient over an ion orbit size (), and the relaxation of the steep gradients during the simulation determines the charge density. 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).
We have noted in Fig.(1) the constant value of the 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 which, as we mentioned before, is small for the parameters we are using. Figures (6) presents this electric field at t=485 and t=495. It shows no oscillation in at = 0 (full curve, essentially equal to zero), a negligible oscillation at = /2 (broken curve), and a small oscillation for the = (dashed-dotted curve), on the high field side of the cylinder. The curve at = is interesting, it shows the oscillating curve for having a constant flat value in the gradient region, and a linear variation in the inner region.
Figure (7) gives at t = 500 the temperature T
Figure(10) presents the plot of the component of the drift, which in our normalized units is written, 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 E
We get indeed a negligible value for. (note that).
We note that due to the smallfield (see Fig.6), the drift has a small oscillating component in the radial direction, with a curve similar to what is presented in Fig.(6). This radial oscillation is negligible around and, and is small on the high field side around, as previously discussed for in Fig.(6).
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 (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 to exist, shows the presence of a small oscillation of on the high field side around (see Fig.(6)), to which is associated a small radial oscillation due to the radial component of the drift, equal to 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).
Abbott B. A. 1966An Introduction to the Method of Characteristics; Thames and Hudson, London.
Ahlberg J. H. Nilson E. N. Walsh J. L. 1967The Theory of Splines and their Applications; Academic Press, New York, N.Y.
Batishchev O. Shoucri M. Batishcheva A. Shkarofsky I. 1999Fully Kinetic Simulation of Coupled Plasmas and Neutral Particles in Scrape-off Layer Plasmas of Fusion Devices. J. Plasma Physics, 61 347 364
Cheng C. Z. Knorr G. 1976The Integration of the Vlasov Equation in Configuration Space. J. Comp. Phys., 22 330 351
Crouseilles N. Respaud T. Sonnendrücker E. 2009A Forward Semi-Lagrangian Method for the Numerical Solution of the Vlasov Equation. Comp. Phys. Comm., 180 1730 1745
Gagné R. Shoucri M. 1977A Splitting Scheme for the Solution of a One-Dimensional Vlasov Equations. J. Comp. Phys., 24 445 449
Ghizzo A. Bertrand P. Shoucri M. Johnston T. Feix M. Fijalkow E. 1990A Vlasov Code for the Numerical Simulation of Stimulated Raman Scattering. J. Comp. Phys., 90 431 457
Ghizzo A. Bertrand P. Shoucri M. Johnston T. Fijalkow E. Feix M. Demchenko V. V. 1992Study of Laser-Plasma Beat Wave Current Drive with an Eulerian Vlasov Code. Nucl. Fusion, 32 45 65
Groebner R. J. Thomas D. M. Deranian R. D. 2001Evidence for Edge Gradients as Control Parameters of the Spontaneous High-Mode Transition. Phys. Plasmas, 8 2722 2732
Hirschman S. P. 1978The Ambipolarity Paradox in Toroidal Diffusion, Revisited. Nucl. Fus., 18 917 927
Knorr G. Joyce G. Marcus J. 1980Fourth-Order Poisson Solver for the Simulation of Bounded Plasmas. J. Comput. Phys., 38 227 237
Manfredi G. Shoucri M. Dendy R. O. Ghizzo A. Bertrand P. 1996Vlasov Gyrokinetic Simulations of Ion-Temperature-Gradient Driven Instabilities. Phys. Plasmas, 3 202 217
Pohn E. Shoucri M. Kamelander G. 2005Eulerian Vlasov Codes. Comp. Phys. Comm., 166 81 93
Pohn E. Shoucri M. 2008Collisionless Diffusion of Particles across a Magnetic field at a Plasma Edge in the Presence of Turbulence. Proc. Vlasovia workshop, in Commun. Nonlinear Sci. Numer. Simul., 13 183 188
Shoucri M. Lebas J. Knorr G. Bertrand Ghizzo. A. Manfredi G. Christopher I. 1997Effect of Viscous Dissipation on the Generation of Shear Flow at a Plasma Edge in the Finite Gyroradius Guiding Center Approximation. Physica Scripta, 55 617 628
Shoucri M. Pohn E. Knorr G. Bertrand P. Kamelander G. Manfredi G. Ghizzo A. 2000Charge Separation at a Plasma Edge in the Presence of a Density Gradient. Phys. Plasmas, 7 2517 2525
Shoucri M. 2001Numerical Simulation of Plasma Edge Turbulence due to ExB Flow Velocity Shear. Czech. J. Phys., 51 1139 1151
Shoucri M. 2002Comments on ‘Evidence for Edge Gradients as Control Parameters of the Spontaneous High- Mode Transition’, Phys. Plasmas, 9 731 734
Shoucri M. Gerhauser H. Finken K. H. 2003Numerical Simulation of the Formation of an Electric Field at a Plasma Edge in the Presence of a Density Gradient. Proceedings of EPS 2003 30th Conference on Contr. Fusion and Plasma Phys., ECA 27A 1 170St. Petersbourg, Russia, July 7-11
Shoucri M. Gerhauser H. Finken K. H. 2004Study of the Generation of a Charge Separation and Electric Field at a Plasma Edge using Eulerian Vlasov Codes in Cylindrical Geometry. Computer Physics Communications, 164 138 149
Shoucri M. Pohn E. Kamelander G. 2005Numerical Simulation of the Collisioless Diffusion of Particles across a Magnetic Field at a Plasma Edge. Proceedings of EPS 2005 32nd Conference on Contr. Fusion and Plasma Phys., ECA 32D 5Tarragona, Spain, June 27-July 1
Shoucri M. 2008aNumerical Solution of Hyperbolic Differential Equations, Nova Science Publishers, 978-1-60456-459-4New York
Shoucri M. 2008bNumerical Simulation For the Formation of a Charge Separation and an Electric Field at a Plasma Edge. Proceedings of EPS 2008 35th Conference on Contr. Fusion and Plasma Phys., ECA 32D 4Hersonissos, Crete, Greece, June 9-13
Shoucri M. 2008cEulerian Code for the Numerical Solution of the Vlasov Equation. Proc. Vlasovia workshop, in Commun. Nonlinear Sci. Numer. Simul., 13 174 182
Shoucri M. 2008dNumerical Simulation of Wake-Field Acceleration using an Eulerian Vlasov Code. Commun. Comp. Phys., 4 703 718
Shoucri M. 2009aThe Application of the Method of Characteristics for the Numerical Solution of Hyperbolic Differential Equations, In: Numerical Simulation Research Progress, S.P. Colombo & C.L. Rizzo (Ed.), Nova Science Publishers, 978-1-60456-783-0New York
Shoucri M. 2009bThe Formation of a Charge Separation and an Electric Field at a Steep Plasma Edge. Proceedings of EPS 2009 36th Conference on Contr. Fusion and Plasma Phys., ECA 32D 4Sofia, Bulgaria, June 29- July 3
Shoucri M. 2011In:Eulerian Codes for the Numerical Solution of the Kinetic Equations of Plasmas. M. Shoucri (Ed.), Nova Science Publishers, 978-1-61668-413-6New York
Shoucri R. Shoucri M. 2007Applications of the Method of Characteristics for the Study of Shock Waves in Models of Blood Flow in the Aorta. Cardiovasc. Eng., 7 1 6
Sonnendrücker E. Roche J. Bertrand P. Ghizzo A. 1999The Semi-Lagrangian Method for the Numerical Resolution of the Vlasov Equation. J. Comp. Phys., 149 201 220
Stix T. 1973Decay of Poloidal Rotation in a Tokamak Plasma. Phys. Fluids, 16 1260 1269
Strozzi D. Shoucri M. Bers A. Williams E. A. Langdon A. B. 2006Vlasov Simulations of Trapping and Inhomogeneity in Raman Scattering. J. Plasma Phys., 72 1299 1302
Waltz R. E. Candy J. Rosenbluth M. N. Hinton F. L. 1999Progress on a Full Radius Electromagnetic Gyrokinetic Code. Proceedings of EPS 1999 26th Conference on Contr. Fusion and Plasma Phys., ECA 23J 1233 1236Maastricht, Holland, June 14-18
Yanenko N. N. 1971The Method of Fractional Steps, Springer-Verlag, New York.