DOE/ER/03077-201

MF-102

Courant Institute of

Mathematical Sciences

Magneto-Fluid Dynamics Division

Axisym metric, Non-ideal MHD

States with Steady Flow

W. Kerner and H. Weitzner

(A

O

1

r

_

•

•

•

U.S. Department of Energy Report

Plasma Physics

October 1983

NEW YORK UNIVERSITY

DISCLAIMER

This report v/as prepared as an account of work sponsored

by an agency of the United States Government. Melther

the United States Government nor any aeency thereof, nor

any of their employees, makes any warranty, express or

implied, or assumes "any lep;al liability or responsibility

for the accuracy, completeness, or usefulness of any

information, apparatus, product, or process disclosed, or

represents that its use would not infringe privately

owned rights. Reference herein to any specific commercial

product, process, or service by trade name, trademark,

manufacturer, or otherwise, does nbt necessarily constitute

or imply its endorsement, recommendation, or favoring by

the United States Government or any agency thereof. The

views and opinions of authors expressed herein do not

necessarily state or reflect those of the United States

Government or any agency thereof.

Printed in U.S.A.

Available from

National Technical Information Service

U.S. Department of Commerce

5285 Port Royal Road

Springfield, VA 22l6l

UC 20

UNCLASSIFIED

New York University

Courant Institute of Mathematical Sciences

Magneto-Fluid Dynamics Division

MF-102 DOE/ER/03077-201

AXISYMMETRIC, NON- IDEAL MHD STATES WITH STEADY FLOW

W- Kerner and H. Weitzner

U . S . Department of Energy

Contract No. DE-AC02-76ER03077

November 198 3

UNCLASSIFIED

L.^

Axisymmetric , Non-ideal MHD States with Steady Flow

W. Kerner* and H. Weitzner,

Courant Institute of Mathematical Sciences, New York University,

New York, N.Y, 10012, USA.

Abstract

Toroidal plasma configurations with steady flow are studied

in the framework of non-ideal MHD theory. The properties of

the resulting set of equations are examined.

The numerical solution of the two-dimensional, non-linear

system appears feasible, although the large variation in

the transport coefficients creates considerable numerical

problems .

c

'>

* Permanent address:

Max-Planck-Institut fur Plasmaphysik , EURATOM Association,

D-8046 Garching, Fed. Rep. of Germany.

- 2 -

O. Introduction

Tokamak discharges can be sustained for several seconds with

an energy confinement time approaching 100 milliseconds. These

experiments are, in general, not terminated by instabilities

in the form of a disruption. A suitable description of the

plasma behaviour is given by the macroscopic model (MDH) .

To reproduce the essential features of the experiment, the MHD

equations have to incorporate the non-linear dependence and

to include non-ideal effects. The full time-dependent problem

is, in our opinion, still far too complex for sufficiently

accurate numerical treatment. We therefore treat the long-time

evolution as an equilibrium problem, where the plasma passes

through a sequence of equilibrium states. By imposing axi-

symmetry, which is adequate for tokamak configurations, the

corresponding problem reduces to determining two-dimensional

equilibria. This is feasible with existing numerical techniques

and computing facilities.

The ideal MHD model with scalar pressure has been remarkably

successful for describing a plasma. Owing to the characteristic

Alfven time scale of the order of microseconds, the plasma cannot

be too far from an equilibrium. This feature is built into

H. Grad's 1 -1/2 D transport scheme /I/, where the plasma passes

through two-dimensional equilibria obyeing the equation

^ p = J X B and the profiles evolve as surface quantities.

This model is, however, only adequate if the plasma flow is

small and if the pressure and density are, basically, surface

quantities.

The large amount of additional heating in the form of neutral

beam injection, which is employed in all major tokamak experi-

ments, acts as a source for toroidal flow with a flow velocity

approaching the ion sound speed. The induced poloidal flow, on

the other hand, does not exceed a certain value, its damping

out leading to a poloidal dependence of the pressure and density

on magnetic surfaces. The ideal MHD model can easily be extended

to include flow. A code for computing such equilibria with

flow has been developed by Kerner and Jandl /2/. (Other

- 3 -

contributions and results are referenced in 111 .)

This model, however, neglects dissipative effects and is not

sufficiently general with respect to boundary conditions

and sources.

If resistivity is taken into account, toroidal equilibria

require a flow for their existence, as is pointed out in

Ref. /3/. Pf irsch-Schluter diffusion, Ref. /4/, sustains such

a pressure-driven flow. For ohmic discharges these flows are

usually small. But the additional heating drastically enlarges

the flow. The viscosity is just as important, and so is the

energy flux due to temperature gradients. Taking into account

the resistivity, viscosity and heat conductivity yield a set

of equations which allow a realistic simulation of an experiment.

The hyperbolic character of the continuity equation requires

sources for its solution. The remaining set forms an elliptic

set. The requirement of an elliptic characteristic is an important

point in our analysis since for mixed systems we expect tremen-

dous numerical difficulties. The essential role of the

continuity equation is then apparent. The choice of boundary

data and source terms should select certain equilibria and

forbid others. A close connection with experimental data is

possible. We thus can hope for results explaining the density

limit being observed in the experiments.

For our macroscopic model transport coefficients are required.

In the collision-dominated regime the coefficients are

explicitly known. More appropriate is the neoclassical regime,

where the trapping of particles is taken into account.

The coefficients have different values along and perpendicular

to the magnetic field. The consequences of this anisotropy

with respect to the numerical approximation are discussed.

The basic assumption in this paper is that the additional heating

deposited in the plasma in the form of neutral beam injection

and wave heating makes a two-dimensional treatment of the

equilibrium problem necessary - owing to the poloidal variation

of important quantities such as pressure and density.

- 4 -

The physical model, its mathematical properties and the

consequences for the numerical solution are analyzed.

The single-fluid model containing non-linear and non-ideal

features exhibits such a complexity that a detailed dis-

cussion of its properties - especially with respect to a

numerical approximation - is both useful and necessary.

We are aware of the fundamental role of the transport co-

efficients for our model. Experimental data show that some

coefficients differ from their classical or neoclassical

values and cause anomalous transport. A very interesting

aspect of our model is therefore its potential for determining

the transport coefficients.

The paper is organized as follows:

Section I presents the fluid equations used throughout the

analysis. In Sect. II the constraint of incompressibility is

incorporated into this model, and the equations are derived

and discussed. The non-ideal, compressible fluid is treated

in Sect. III. After the energy equation, the validity of the

incompressibility assumption is discussed. The general pressure

tensor is derived in Sect. IV, with only covariance and symmetry

properties being used. Such a derivation is useful for under-

standing the macroscopic features of this tensor, especially

if one aims to use a simpler form instead of the full tensor.

The transport coefficients are listed in Sect. V. The discussion

and the conclusions are then finally presented in Sect. VI.

- 5 -

I. Fluid Equations

To begin with, we list the MHD equations for a single-fluid

theory. This model relates the density 3 ' velocity u ,

scalar pressure p, pressure tensor P , internal energy e and

the magnetic field B:

(1) Continuity: ^ ^ ^ ^- l%Si I * 9s

(^ ^ ^-"^^ U = - Vp V T xS * \/-P

(2) Momentum: Yv\3 1^-^ * U.-V ) U = - Vp v ^

(3) Energy

S l^ * u-^") ^ - -^ " ^-[^ '^"-^^

where 2 denotes the heat flux and ^ and xi^ mass, reap, energy sources

(4) Ohm's law: ^\ ^ '^uA'^^ ^ ^*^''^'

where Oft is the resistivity tensor, tt\^ the Hall constant,

E the electric field and J the current.

(5) Maxwell: ^ ■- - V* E

By means of the thermodynamic relations the pressure is derived

from the internal energy and entropy S:

(32') V(\^suV V\- \1% + V^x\)^.\?i - \/-}X^£;i,

(33') \ \) (\^.5 0^,^"^ X \)X-^i = '^» ^•j^^'^.o. - eocp. ^p.,>^^'X,j^-^

If there is no flow, then it holds that ^s and w * 0.

From eq. (32) it follows that Y = 'A (^^) and eq. (26) reduces

to the Grad-Shaf ranov equation:

(34) -p,^ - YX^^h^ = V- (^^/y^^ .

The system (30' - 33') decouples into two pairs of equations.

The functions H* and 'K. are determined by eqs. (31') and (33').

Once ^ and % are known, the functions % and to are the

solutions of the linear system 30' and 32'. In the toroidal

system PO - 33) the functions % and Co also occur in eq . (33).

If % and Co are known , the functions ^ and 'X are determined

by eqs. (31) and (33). This defines an iteration scheme for the

numerical solution.

We examine the equations for the poloidal quantities X and ^

with given U> and 'X, for the toroidal case, eqs. (31) and (33).

For non-zero viscosity, ^ -^ O , this is a fourth-order

elliptic equation for \ , together with a second-order elliptic

equation for ^ . For this system only boundary data can be pre-

scribed and no profiles can be imposed, except the dependence

Q 5= 5^V( on the continuity equation. The other two equations

for the toroidal quantities to and % , eqs. (30) and (32) , with

given 'X and ^> also form an elliptic set for Oi and % .

- 13 -

The case with purely toroidal flow, i.e. \»0, is possible

if To'^o - ^c^o ^ ^r^V- iV^/r^^

which determines M^ . Ohm's law for % , eq. (30) , and eq. (32)

form a closed set for tc,^ and H* . Now ^ is no longer a function

of \ , but is arbitrary. The pressure balance, eq. (26) , where

(35) V(p . ^,V ^ W » V[^.)['^i^^^-X^)-0

is then highly degenerate. The two components of this equation

determine the remaining thermodynamic unknowns = , the system has quite different proper-

ties. The momentum equation (33) then reduces to

(36

The solvability condition for this equation is less clear. For

given H^ one has to solve a hyperbolic equation for ^\. This

implies that on closed field lines there is a solvability con-

straint. The free data on ^JX are needed for this constraint.

The discussion of the energy equation in the next section makes

it clear that the assumption of incompressibility is not satis-

fied for equilibria with flow. On the other hand, the system

(30 - 33) is quite complicated. We therefore conclude that the

considerable effort for its numerical solution may not be

worthwhile .

- 14 -

III. Compressible Fluid

In this section the energy equation (3) is included and the

assumption of incompressibility is given up. At first the general

Ohm's law, eq. (4) , is discussed. The resistivity tensor M is

defined as

(37)

^ = /V^^ X - i^i - ^^^^^ y

where 1£ is the unit matrix, b = B/\B| the unit vector in the

direction of the magnetic field and 'VV), and 'V^j^^ the values of

the resistivity along and perpendicular to the magnetic field.

This implies that

(38) V^.^= 'V^il^-^„^ - %4, - /V\i^l ^ 'V^4.

where w = ( ^- b^ b

is the component of the current parallel to the magnetic field

and J^ the perpendicular one.

Using eq. (11) for the current yields

(39) ^, = ( w . w - \^o r) /^'^" ^ = ^u ^

With the introduction of the electric field

E = - Vc^ + toEo ^& ,

where C^ is the scalar potential (see eq. (16)), Ohm's law

takes the form

(40) f^^-^ ^ {y- '^x') l^ -^ ^H 4^ ^ ' -^ ' ^0^0^ ^ ^-^^

The projection onto ^^ yields

- ToEo - M. UH'

■^ '^^ Wx VX-ue.

- 15 -

Next we examine the \}^ - component of the Maxwell equation

which is equivalent to

and, furthermore,

We utilize eq. (41) to eliminate "\q and with

^"^ V|JVM'\^ V /v^„ %^

we eventually get the result

(42) \i. [f^ x)^") + wx ^Cr^V^e =

= M-^ ^ VX A (/V)„-M^V^'"V^ ^\ , r i

'v^^ \^V\^ 4 ^ V^^'OX-VQ ^ x\Jy-^[t]- Y-^rT ~ }^T^>

and the z- component

(54)

Clearly this system is elliptic for non-zero viscosity.

Finally, we examine the energy equation. Introducing the

temperature instead of the internal energy, we obtain the

following energy equation:

(55)

The heat flux is related to the temperature

(56) -^ =r kVT - W, ^J * vcj^^J - v:^ bxT/T

where parallel VCy, and perpendicular K^ and V^^ are thermal

conductivities. The energy flux due to the viscosity is given

by

i? : vu = s ip.^ ^;

- 19 -

Inserting the terms of the pressure tensor yields

-f.u^a^ ii].

With these expressions the energy equation, which determines

the temperature, reads

(58) f^(^V)l + ^T V-li - VUVT^ - /V)i^i * /V|„^„ + f -.Vu

The condition of incompressibility can be derived from the

energy equation. Equivalent to the energy relation (3) or (55)

is the entropy equation

With p = p ( 3 , ^ ) we can easily write an equation for the

pressure as

or

In the case of non-dissipate flow, where the right-hand side

vanishes, we recover the usual conditions for the validity of

incompressible flow, viz. the flow velocity must be small and

the pressure variations must be small compared with the mean

thermodynamic pressure. If dissipation is included, it follows

that, in addition to the above conditions, the dissipation must

not be so large that the right-hand side becomes comparable

with the left. Physically, the dissipation-induced pressure

20 -

variations must be small compared with the other pressure

variations in the system.

In a confiend plasma with substantial pressure variation, the

assumption of incompressibility is not justified but is rather

poor.

At the end of this section we summarize the set (S2) of com-

pressible equations for axisymmetry:

(1) ^ ^-U * ^. ^ « g5j

(41)

(42) "^-[ji "0^) V VH'^ V(7i).^ =

= ^. j ^f MX * -^'^^^ VH> -^

V^^V^ (toEo - juk-W ^ 'viH ^e)x^cv■u'xv■^'■

(52)

(58)

^^ 'v^il^H'l^ + 'v\„ X^

i ( - 1 . \ ^

(53) vw^M*l()A^-^\K-^x^v.Oe,

(53-) Wv^A^-Vv^ ((i - ^],^^ Y^^.-'^-h^-^ * ii^'^^\C^I^^'^-^^^

(54') A»^3 u-^UJ * V^ X x/x-^i -V-^^w,

(58')

- 22 -

If there is no flow, i.e. M = O / it follows from eq. (53)

that X^XC^^ and eq. (51) reduces to the Grad-Shaf ranov

equation .

In the case with purely toroidal flow, i.e. ox •* v ^ o ,

it follows that SX.-V - and ^-xa - 0. The continuity equation

then does not determine the density. Equations (41), (42) and

(53) form a closed set for ^' ,X and co . The two remaining

poloidal components of the momentum equation determine the

thermodynamic quantities «^ and p . As discussed for incom-

pressible flow, in the case with viscosity, the momentum imbalance

may act as a source for poloidal rotation even if no rotation is

imposed at the wall.

For zero viscosity, u. « O , the system has quite different

properties. The momentum equation reduces to

W\

^oA- Vu = - Vp - -^ \J^ - ^ ^'>^ + ve vxx^H^.^6

For given ^ one has to solve a system of mixed type with

sonic transition, which causes difficulties as discussed in

Sect. II.

- 23 -

IV. General Pressure Tensor

In the analysis so far a simplified pressure tensor has

been taken into account. A strong magnetic field, however,

requires a more general tensor than that from ordinary fluid

theory. This tensor is derived from microscopic theory;

see, for example, Braginskii / 5 /.

In our opinion it is useful to formulate the stress-strain

relations in a manner independent of the origin of the viscous

forces, but dependent only on covariance properties and simple

physical bypotheses. In particular, we assume that the stresses

are linear functionals of the strain matrix S :

(59) g.. = ^^ 4 — *

We can generalize this treatment to allow the coefficients

of the linear relation to depend on the invariants of that form,

namely on the traces Tr _S, Tr S_^ and Tr _S ^ . Additionally, we

assume that the stresses can depend on the polar vector B,

but on no other quantities. Again, the coefficients in the

linear stress-strain relation can be generalized to depend on

the invariants found from S and B. The stress matrix P must

be symmetric and invariant under inversion of coordinates,

since B is not. We shall construct the corresponding form and

show that it is essentially equivalent to that of Braginskii

with arbitrary coefficients. We shall adjoin three additional

hypotheses. We postulate that the viscous forces must increase

the entropy of the system and we assume that the resulting

viscous equations must generate an elliptic stress v ■ ^ •

Finally, for convenience, but not essentially, we assume that

similar to an ordinary fluid there is no bulk viscosity, or

that the viscous stress only depends on

.60, ^J.. . ^,. - i 6,1; I =i:^|-:: - |=J.i(^-^

- 24 -

It is convenient to employ the vector B in two distinct forms,

We use B interchangeably as a row or column vector and we

introduce the standard matrix from electromagnetic theory.

(61)

B> '

e>u\

/

which is invariant under inversion of coordinates,

For any vector V it holds that

^ V = B^ '^ y

while

and

The stress tensor must be even in B and of any order in B .

Although we shall drop the bulk viscosity, for completeness

we give the most general stress tensor proportional to (V-x\ ) .

The most general form consistent with the covariance properties

is clearly

where uj* and u^ are the two bulk viscosity coefficients,

which we set equal to zero.

We construct the stress tensor from W^i , B and B.

We systematize the procedure by taking terms of increasing

degree in B and B. Clearly, the only term of degree zero is