Detailed derivation of axisymmetric double adiabatic MHD equilibria with general plasma flow

Detailed derivation of axisymmetric double adiabatic MHD equilibria with general plasma flow

Roberto A. Clemente

Instituto de Física Gleb Wataghin,

Universidade Estadual de Campinas

13083-970, Campinas, São Paulo, Brazil

Ricardo L. Viana

Departamento de Física,

Universidade Federal do Paraná

81531-990, Curitiba, Paraná, Brazil

Received 11 December, 1998

We derive a system of equations to describe stationary axisymmetric MHD equilibria characterized by toroidal and poloidal flows as well as plasma anisotropy due to strong magnetic field and eventual auxiliary heating methods. The system consists of a nonlinear partial differential equation for the poloidal magnetic flux function and an algebraic Bernoulli-type equation, which relates plasma density with several surface functions. We analyse the ellipticity of the equation and the plasma density bifurcation as the Alfvén-Mach number is changed.

I Introduction

Plasma rotation in tokamaks has been observed during non-ohmic current drive and neutral beam injection, when a net momentum is imparted to the plasma [1], [2]. Plasma rotations with Mach numbers of order unity have been detected in tokamaks at almost all operating regimes, both in toroidal and poloidal directions [3, 4, 5]. Plasma flow has also been observed in reversed-field pinches (RFP), showing a high velocity shear between core and edge plasma rotation [6]. In field-reversed configurations (FRC), strong plasma azimuthal rotation is responsible for an n = 2 instability that may destroy plasma confinement [7].

From the point of view of magnetohydrodynamic (MHD) equilibrium theory, in the absence of plasma resistivity (ideal theory), magnetic field lines lie on magnetic flux surfaces that rotate rigidly with the plasma, according to Alfvén's theorem [8]. Equilibria are characterized by surface functions, i.e., physical quantities that are constant on nested closed surfaces, with the topology of tori. Examples of surface functions include the volume enclosed by a toroidal magnetic surface, the magnetic flux of the poloidal magnetic field component, and the field line rotational transform.

We will assume that the plasma is in a MHD stationary equilibrium state, which means that all partial time derivatives vanish, but allow for stationary velocity profiles, including plasma flows in toroidal and poloidal directions. The momentum balance equation in this case leads to an equilibrium partial differential equation for the poloidal flux, which involves other surface functions. Earlier studies to derive an equilibrium equation with plasma flows were conducted by Zehrfeld and Green [9] and Morozov and Solov'ev [10].

Maschke and Perrin [11] in 1980 were able to solve a MHD equilibrium equation for purely toroidal flows in axisymmetric plasmas. In addition to the ideal MHD equations, they have made some thermodynamical assumptions: the plasma was taken to be an ideal gas, with larger thermal conductivity along the magnetic field lines than across them, so they have considered the plasma temperature as a surface quantity. Situations where the flow is adiabatic allow the use of plasma entropy as a surface quantity, instead of the temperature.

There are few analytical or semi-analytical solutions for purely rotating equilibria in the literature, for cylindrical [11, 12, 13] and spherical geometries [14]. The latter configuration turns to be interesting to explore astrophysical plasmas, as in magnetic stars. All these solutions share some properties: the magnetic axis (degenerate flux surface where the flux is an extremum) is displaced outwards due to the centrifugal effect, and the isobaric (constant pressure) surfaces no longer coincide with the flux surfaces.

The problem of combined toroidal and poloidal fluxes has been first adressed by Zehrfeld and Green [9], who assumed simple adiabatic flow. Later, Maschke and Perrin [16] have solved such equilbrium equations in the limit of small ratio of poloidal to toroidal magnetic field and small beta. Further developments were due to Kerner and Tokuda [17], Zelazny et al. [19], and Tasso and Thromoulopoulos [18]. A common feature of all these models is the need of supplementing the partial differential equation for the magnetic flux with a Bernoulli-like algebraic equation for the density which contains several other surface quantities.

Zelazny and collaborators [19] have obtained numerical solutions of the differential and algebraic equations (called EGSS, or ''extended Grad-Schlüter-Shafranov'' system), by using the inverse method and Fourier decomposition. They have observed that the differential equation alternates between elliptic and hyperbolic when the Alfvén-Mach number of the poloidal flow (defined as the ratio between the poloidal flow velocity and the Alfvén velocity relative to the poloidal magnetic field) increases. This fact was already analytically investigated by Hameiri [20], and it is of fundamental importance if a well posed boundary value problem is to be solved. If the equation is not elliptic no closed equilibrium configurations are possible.

So far, we have reported studies on isotropic stationary plasmas only. On the other hand, pressure anisotropy due to intense magnetic fields and auxiliary heating methods is an important factor to be taken into account if a general MHD equilibrium equation is to be written down. It is the case, for example, when the auxiliary heating power in toroidal devices is greater than the ohmic heating power. This possibility was already considered as early as in 1961 for static equilibria [21], but only recently analytical solutions for static and rotating anisotropic plasmas were obtained [12, 22]. Clemente has found a method to obtain axisymmetric solutions of the MHD equations with plasma anisotropy in terms of the corresponding solutions of the isotropic set of equations for the static or rotating case.

In this paper we present a more general axisymmetric MHD stationary equilibrium model, which takes into account both poloidal and toroidal plasma flows as well as pressure anisotropy. Besides the usual surface quantities, the poloidal magnetic flux and current functions, other auxiliary surface functions will be defined. The poloidal magnetic flux nonlinear partial differential equation is complemented with a Bernoulli-type algebraic equation which generalizes that introduced for the isotropic case [9, 10, 16]. The system reduces to the previously studied cases when anisotropy tend to vanish. We also discuss the ellipticity of the differential equation and study plasma density bifurcations when the Alfvén-Mach number is varied. This could be related to the problem of L-H transition in tokamaks.

This paper is organized as follows: in the second section we outline the basic MHD and thermodynamical equations to be used and the representations for the magnetic and velocity fields in cylindrical coordinates. Section 3 introduces the auxiliary surface quantities necessary in our model, and in section 4 an energy theorem is derived. The momentum balance equation is derived in section 5, whereas in section 6 we study some particular cases. Section 7 is devoted to an analysis of the ellipticity of the equation, followed by a discussion on plasma density bifurcation. Our conclusions are left to the last section.

II Basic Equations

Let us consider an ideal plasma of electrons and singly charged ions in stationary equilibrium, all partial time derivatives will vanish and nonzero plasma velocity will be allowed. The corresponding ideal anisotropic MHD equations, in standard units, are [23]

is the mass density, n is the particle number density and m_e, m_i are the electronic and ionic masses, respectively. v, E, B, and j are the velocity, electric field, magnetic field and plasma current density, respectively.

If the plasma is subjected to a strong magnetic field and ions have a small Larmor radius with a rapid gyromotion, the pressure tensor may be no longer isotropic. Introducing different plasma pressures: p_^across the magnetic field and p_|| along it, we can introduce the Chew-Goldberger-Low form for the pressure tensor [24]

is a measure of pressure anisotropy, is the identity tensor and B² = B·B.

In order to get to a closed system of equations we must introduce some constitutive thermodynamical assumption so as to obtain an energy equation. In earlier works, where the pressure was supposed isotropic, it has been assumed that either the entropy S or the temperature T are surface quantities. In the former case, an adiabatic equation of state has been used

where g is the specific heat ratio. The specific enthalpy h satisfies the thermodynamical relation

which gives by integration, using (10)

The internal energy e = h/g satisfies the Gibbs' equation

In the case of anisotropic pressure, however, we must use the double adiabatic equations [23]

Using (1), (5) and (6) these equations can be combined, there resulting that

Let us consider cylindrical coordinates (r,f,z), where f is an ignorable quantity, i.e., our configuration exibits symmetry with respect to this coordinate. A magnetic field representation satisfying Eq. (3) can be written in terms of two scalar surface functions Y and I (poloidal flux and current functions, respectively)

We also assume that the plasma can have a rotation around the symmetry axis in both toroidal and poloidal directions. A velocity representation which satisfies Eq. (1) is

where G is a stream function, and F is the corresponding stream function for poloidal flow.

III Auxiliary Surface Quantities

Taking the cross product between (21) and (20), and using axisymmetry we have

From Ohm's law, the azimuthal component of Eq. (6) vanishes. Using (22), it implies that the jacobian of the transformation (r,z)®(F,Y) is identically zero, i.e., F = F(Y) is also a surface function. Hence ÑF = F' ÑY, where the primes denote differentiation with respect to Y.

Now, inserting (22) into Faraday's law, Eq. (5), and using (6), we obtain ÑW × ÑY = 0, where we have defined

and it is also a surface quantity. Physically W(r,z) is a toroidal angular velocity for a rigid surface rotation. This is commonly called ''Ferraro's isorotation law'' [8].

Combining (23) with (21) we have for the plasma velocity the following expression in terms of the magnetic field

Computing the divergence of the pressure tensor, as given in the form (8), and substituting it into the momentum balance equation for equilibria, Eq. (2), we have

where we have used Ampère's law (4). Rearranging terms we have

where we defined the tensor

Dotting (26) with , we find that the divergence of has zero f-component. There results, by using (1) and (2), that (25) can be rewritten as

or, by means of the assumed representations of B and v,

where a new auxiliary surface quantity has been defined as

Let us discuss a special case in order to get a physical insight on the possible motions allowed by the present treatment. In the case of a field-reversed configuration (FRC) we have I = 0, so G = L/F^¢. On the other hand, from (23), we have that G = Wr². However, while L/F^¢ is a surface quantity, it is easy to prove that Wr² is not, hence the equality

is not valid for both L ¹ 0 and F^¢¹ 0.

Therefore, there are only two cases where the equality (31) is compatible with (30) in the a FRC case: i) W ¹ 0 and F^¢ = L = 0, which means only azimuthal rotation (toroidal plasma flow); or ii) W = G = 0 and F^¢¹ 0, i.e., only poloidal flow. Since in FRC some region of vanishing magnetic field, where s_- should vanish, will exist necessarily, eqs.(19) will predict different expressions for the density in the case where perpendicular and parallel plasma pressures tend to be identical. Such expressions will also be incompatible with the corresponding density arising from Bernoulli equation, and we must conclude that both toroidal and poloidal plasma flows are not allowed simultaneously when I = 0. Moreover, poloidal flow will be possible only within a simple adiabatic ideal MHD model, since in such case the density is uniquely determined. We remark that earlier works [10, 11, 20] have dealt with both flows simultaneously also when I=0. When only poloidal flow is present, FRC like equilibria can not be anisotropic in the Chew, Golberger and Low sense. However, they can be anisotropic if only toroidal rotation is present, as it has been recently shown[22].

IV An energy theorem

Dotting Maxwell equations (5) and (4) with B and E, respectively, and subtracting both sides we have Poynting's theorem

Ñ · (E ´ B) = -E · j, (32)

where S = E ×B represents the net energy flux across the surface of a plasma element, and the second term is the ohmic heating term. Using (6) and (2) the latter can be written as

E · j = rv · [( v · Ñ )v] + v · Ñp_^ + v · [B(B · Ñs_)] + s_v · [(B · Ñ)B]. (33)

Using the double adiabatic equations (15,16), this term can be rewritten, after some algebra, in the form of an energy theorem for anisotropic stationary equilbria

Now, using (20) and (6) we have for the divergence of the Poynting vector

and considering that

we have that Equation (34) implies the following relation

and we see that the expression within the brackets is a surface quantity

since B·ÑQ = 0. This is a Bernoulli-type algebraic equation that implicitely defines the plasma density r in this double adiabatic model.

V Momentum balance equation

Consider the equation of motion (2), rewritten by using Ampére's law (4) and opening the divergence of pressure tensor given by (8). It reads

The first term in the left-hand side can be written as

Dotting it with ÑY, and using the representation (24) for the velocity, we have

Similarly, the term involving the pressure tensor would give

The calculation of the last term is analogous to the static, axisymmetric case, and reads

where the Grad-Shafranov operator D^* is given by

Plugging eqs. (40),(41), and (42) into (39), and after some rearrangements, we have the equation

But, using (23), we have that

and substituting it in (44) we have

The term within brackets in the second member of the above equation may be written as

Using the double adiabatic equations (17,18,19), we can show that f and g are actually surface quantities, so that we finally get to the momentum balance equation describing stationary MHD equilibrium of an axissymetric plasma with anisotropy and general flow

Equation (49) together with (38) form a set of coupled equations whose solutions depend on six hypotheses on the functional dependence with Y of the surface quantities W, F^¢, f, g, L and Q.

VI Particular cases

Once we have obtained the general equations for anisotropic equilibria with plasma flows, let us consider some special cases. First, consider the isotropic case, where we have p_^ = p_|| = p, or s_- = 0. In this case the term

or, identifying g = 5/3 and using the adiabatic law (10), as

The equation (49) for isotropic plasmas (simple adiabatic model) assumes the form

Making the following changes of notation: I® Á, W® w, F^¢® u, Q® Q_s, D^*® , Y® y, and L ® Á₀+L₁ we get to the equation presented by Maschke and Perrin [16] for toroidal and poloidal fluxes, to be supplemented by the Bernoulli-type equation (38), which in this limit reduces to

The electric field is E = - ÑF_E, in which the electric potential is also a surface quantity, and it is easy to prove that F_E(Y) = W.

We can further simplify things by switching off the poloidal flux by making F^¢ = 0 in (53), with both G and W nonvanishing, and isotropic pressure. In this case, using (30), we have that L = I, which is the poloidal current flux. The stationary equilibrium equation becomes

is nothing but the plasma enthalpy (12) with a sort of centrifugal correction. Using the thermodynamical relation (14) we rewrite Eq. (55) in the form

which can be identified with the earlier equation of Maschke and Perrin [11] for description of pure toroidal flows with the entropy S as a surface quantity, by making the following changes of notation: I® , W® w, Q® q, D^*® , Y® F, and r® R. We remark that an equation in which the plasma temperature T would be a surface function cannot be derived from an adiabatic model such as used here.

Finally, in the limit of vanishing flow, we have that p^¢ = r(Q^¢-TS^¢), with Q(Y) = h, so the equation is simply the Grad-Schlüter-Shafranov equation for the static axisymmetric case

VII Ellipticity of the equilibrium equation

Let us return to the set of equations characterizing a general stationary MHD equilibrium with anisotropic pressure and poloidal and toroidal flows. It consists of a nonlinear partial differential equation for the poloidal magnetic flux Y, Eq. (49), and an algebraic Bernoulli equation (38), which defines the plasma density r.

An important issue formerly discussed by Hameiri [20] in the case of isotropic pressure and simple adiabatic flow, are the conditions of ellipticity for the partial differential equation. In this case one has to analyse the coefficients of the second order derivatives of Y in Equation (49). This is important because if the equation is elliptic, a boundary value problem may be posed, after defining proper boundary conditions at the plasma boundary and/or external conductors. In such case, it is possible to solve it by using standard numerical (or, hopefully, analytical) techniques, as done by Zelazny and collaborators [19] in the isotropic case. Otherwise, the existence of hyperbolic regimes implies that some kind of shock wave should be generated in the plasma and no closed equilibrium would be possible.

Second-order derivatives (with respect to r and z) appear in (49) both in the term containing the Grad-Shafranov operator D^*Y and in the term formed by the dot product of ÑY with F^¢Ñ(F^¢/r)+Ñs_-. So, we may write the quadratic form

where the coefficients depend on |ÑY|², F^¢2 / r and s_-. The determinant of the symmetric matrix of these coefficients is

(60)

Equation (49) will be elliptic if Det > 0, parabolic if Det = 0, and hyperbolic if Det < 0.

Considering the corresponding terms in (49) there results

where we defined

and, after using (8) and (23), we introduced

and the dots represents differentiation with respect to |ÑY|².

Now using the Bernoulli equation (38), and after a lengthy but straightforward calculation, we find that the determinant is

If s_- = 0, the determinant (65) coincides with that formerly derived by Hameiri for isotropic pressure [20].

For large aspect Tokamak equilibria the poloidal field component | ÑY| /r is small when compared to the total magnetic field B. If the plasma betas (64) (ratios between kinetic and magnetic pressures) are also small enough, we can develop the denominator of Det to the first order in ,b_||, and b_^.

There results, in this case, that Det is positively definite when the following conditions hold

As in the simple adiabatic case we have three regions of ellipticity for the parameter , which can be interpreted as the square of the poloidal Alfvén Mach number. We note that the second region is separated from the first one by a very small gap.

We stress that plasma density, in the present double adiabatic model, is implicitly determined by the Bernoulli equation (38). Substituting there the double adiabatic equations we find, after using the definitions of f and g provided by (48), two auxiliary equations for the density:

where it must be taken into account that B depends also on r . The plus or minus sign distinguishes different branches for r. Their equivalence defines one of the bifurcation points for the density.

Using the expression (38) for Q, the vanishing of the square root in (70) defines a bifurcation point for the density corresponding to

where the Alfvén Mach number was taken as the bifurcation parameter. As it can be seen from the ellipticity conditions, such bifurcation point falls just inside of the second region of ellipticity for the nonlinear partial differential equation for the poloidal magnetic flux. Its connection with an eventual transition similar to the L-H transition in tokamaks (where an abrupt change from low density to high density equilibria is observed when the power of neutral beam injection is increased) is quite suggestive.

VIII Conclusions

Plasma rotation in fusion plasmas generated in Tokamaks and Reversed Field Pinches is a subject of growing interest, since it should be related to the existence or not of equilibrium configurations, as well as the presence of some types of instabilities. Hence, it would be necessary to develop theoretical models to describe the main effects of rotation without resorting to direct numerical 3-D MHD computer simulations. This can be done by using a double adiabatic ideal MHD model with axisymmetry, which includes plasma flow in both poloidal and toroidal directions and pressure anisotropy, which are possible conditions in the present generation of strong power injection and high magnetic field tokamaks.

Our theoretical model, in accordance with previous approaches with less generality, is constituted by a nonlinear partial differential equation for the poloidal magnetic flux, plus an algebraic Bernoulli-type equation which involves thermodynamical assumptions. The system of equation reduces to the previously studied cases in the appropriate limits of vanishing anisotropy and/or plasma flows.

The differential equation involves six surface quantities which depend on the poloidal flux, so in order to obtain a solution, hypotheses should be specified for all these functions. However, before proceeding to a solution of a boundary value problem, it should be investigated whether or not the differential equation is elliptic. This has been examined in this paper by studying the values of the poloidal Alfvén-Mach number for which the equation is elliptic when a tokamak-like equilibrium is assumed. We have found three intervals in which the condition of ellipticity is fulfilled.

We have proved that for a field reversed configuration, where the magnetic field vanishes at some place inside the plasma, and anisotropy is considered, the plasma flow must be purely poloidal or purely toroidal. Hence, the general equations here derived cannot be used in such cases without proper modifications. However, for tokamaks the theory should work.

Another result is concerned with the definition of plasma density, which is implicitely fixed by Bernoulli equation. Our analysis shows the existence of two branches of density, characterizing a density bifurcation when the poloidal Alfvén-Mach number varies continuously. The critical bifurcation parameter is located in one of the three ellipticity regions.

This could be related to the transition between a L-mode to H-mode, observed in present tokamaks with high power neutral beam injection, which is characterized by a sudden increase in the plasma density and improved confinement leading to a better performance. In spite of intense theoretical effort, the L-H transition in Tokamaks is poorly understood, so our results could shed some light on this subject, indicating that the combined action of pressure anisotropy and plasma flows could be the triggering mechanism for the transition to the H-mode regime. However, further work has to be done in order to give more precise answers to this question.

Even for low beta tokamaks our general system of equation seems to be too difficult to allow an analytical solution, even using very simple profiles for the surface functions. However, provided we stay within the ellipticity intervals, a numerical solution could be sought for, after defining proper boundary conditions at the plasma boundary and external conductors. This is beyond the scope of the present paper.

Acknowledgements

We acknowledge useful discussions with Dr. George Morales (UCLA). This work was made possible with partial financial support of the following government agencies: FAPESP (São Paulo), CNPq (Brazil) and NSF (United States, through an international cooperation project).

References

[1] M. G. Bell, Nucl. Fusion 19, 33 (1979).

[2] S. Suckewer, H. P. Eubank, R. J. Goldston, E. Hinnov, E., and N. R. Sauthoff, Phys. Rev. Lett. 43, 207 (1979).

[3] K. Brau, M. Bitter, R. J. Goldston, D. Manos, K. McGuire, and S. Suckewer, Nucl. Fusion 23, 1643 (1983).

[4] R. C. Isler et al., Nucl. Fusion 23, 1017 (1983).

[5] R. Taylor et al. in: Proceedings of 13th IAEA Conf. on Plasma Physics and Controlled Nuclear Fusion Research (Washington, 1990), paper CN-53/A-6-5 (Vienna: IAEA).

[6] L. Carraro, M. E. Puiatti, F. Sattin, P. Scarin, and M. Valisa, Plas. Phys. Control. Fusion 40, 1021 (1998).

[7] R. K. Linford, W. T. Armstrong, D. A. Platts, and E. G. Sherwood, in: Proc. IAEA Conf. on Plasma Physics and Controlled Nuclear Fusion Research (Vienna: IAEA, 1979), Vol. II, p. 447.

[8] B. Chakraborty, Principles of Plasma Mechanics (Wiley, New York, 1990).

[9] H. P. Zehrfeld and B. J. Green, Nucl. Fusion 10, 251 (1970).

[10] A. I. Morozov and L. S. Solov'ev, in: Reviews of Plasma Physics (Ed.: M. A. Leontovich), Vol. 8, Chap. 2 (Consultants Bureau, New York, 1980)

[11] E. K. Maschke and H. J. Perrin, Plasma Phys. 22, 579 (1980).

[12] R. A. Clemente and R. Farengo, Phys. Fluids B 27, 776 (1984).

[13] O. Missiato and J. P. Sudano, in: Proceedings of the First Latin-American Workshop on Plasma Physics and Controlled Nuclear Fusion Research (Sociedade Brasileira de Física, São Paulo, 1982), Vol. 1, p. 264.

[14] R. L. Viana, R. A. Clemente, and S. R. Lopes, Plasma Phys. Control. Fusion 39, 197 (1997).

[15] C. Plumpton and V. C. A. Ferraro, Ap. J. 121, 168 (1955).

[16] E. K. Maschke and H. J. Perrin, Phys. Lett. A 102, 106 (1984).

[17] W. Kerner and S. Tokuda, Z. Naturforsch. 42a, 1154 (1987).

[18] H. Tasso and G. N. Thromoulopoulos, Phys. Plasmas 5, 2378 (1998).

[19] Z. Zelazny, R. Stankiewicz, A. Galkowski, and S. Potempski, Plasma Phys. Control. Fusion, 35, 1215 (1993).

[20] E. Hameiri, Phys. Fluids 26, 230 (1983).

[21] C. Mercier and M. Cotsaftis, Nucl. Fusion 1, 121 (1961).

[22] R. A. Clemente, Nucl. Fusion 33, 963 (1993); R. A. Clemente, Plasma Phys. Control. Fusion 36, 707 (1994).

[23] W. M. Stacey Jr., Fusion Plasma Analysis, (Wiley, New York, 1981).

[24] G. F. Chew, M. L. Goldberger and F. E. Low, Proc. Roy. Soc. (London) A236, 112 (1956).

Publication Dates

• Publication in this collection
  07 Aug 2000
• Date of issue
  Sept 1999

History

• Received
  11 Dec 1998