MHD simulation of jet

Computer Physics Communications 24 (1981) 343 —354 North-Holland Publishing Company

343

MHD SIMULATION OF JET M.F. TURNER and J.A. WESSON Cuiham Laboratory, Abingdon, Oxon, 0X14 3DB, UK (Euratom/UKAEA Fusion Association)

1. Introduction We are currently developing a code to provide a 3-dimensional MHD description of a toroidal plasma. In particular we aim to simulate the behaviour of the JET experiment, hence the development work is being carried out in JET geometry with its tight aspect ratio and D-shaped cross-section. The objective is to produce a code which, in its final form, can be used to simulate the 3D time dependent evolution of any toroidal, curved boundary plasma and in particular to investigate the development of the various modes of instability, We present here an account of the planned development of the code together with some of the early

Phase I A very simple model in which diffusion of electric and magnetic fields only is simulated with the electrical conductivity specified. In this case no plasma motion is allowed so the plasma is represented as a rigid conductor. The model at this stage is axisymmetric and 2-dimensional. Phase II In this first upgrade of the model we allow motion to take place within the plasma but the density is assumed to remain constant. A simplified energy equation is used at this stage but ohmic heating is included.

results.

Phase III

2. Basic framework

We now further develop the model by introducing the complete energy equation, including additional heating and energy losses, and also the continuity equation. The electrical conductivity is now calculated using the resulting temperature.

We employ an explicit 2-step Lax—Wendroff scheme on a rectangular Eulerian mesh in the minor cross-section to follow the time development of a one fluid MHD model. The toroidal dimension will be treated by means of a Fourier expansion. The calculation is on the resistive time scale and special techniques must be used to avoid fast hydromagnetic waves. Rapid transport along magnetic field lines will be included allowing a description of the loss of confinement resulting from instabilities.

Phase IV

With the model as defined in phase III the 3rd dimension is introduced.

4. Development progress Phase I

At this stage the model was a straightforward axi-

symmetric, 2-dimensional representation of the Max3. Development procedure In order to maintain confidence in the model the development has been divided into several phases as follows:

0010-4655/8 1/0000—0000/$02.75 © 1981 North-Holland

well equations: V XB

=p~j,

VXE= VB=O,

B/at

M.F. Turner, J.A. Wesson / 3D MHD simulation of JET

344 10

1(~

/

0.8

/

0.8

/ 06

0.6-

//

E~

0

/

0.4

04-

/ /

02-

C

0

02-

I

0.2

0.4 IC-

0.6

0

0.1

)

Fig. 1. Penetration of electric field to the axis of a cylindrical plasma in form of rigid conductor. Solid line shows values from the exact result. E50) [ —i1—2 E~o L

c

s-’ /c~7~t\Jo(a0r/a) U exp—I————2—I n—i \,.~oa ! a0J1 (cs,~)

where the on’s are the zeros of theJ0 Bessel function. Circles on the dashed line give computed results with 24 X 24 mesh. Crosses give results obtained with a 48 X 48 mesh.

2.75 171

~

12.1 6 21

03

I 0.4

I 0.5

IC_lit

Fig. 2. Penetration of electric field into JET shape plasma in form of rigid conductor (20 X 32 mesh). The electric field is taken atR = 3.2 m, z = 0.

Phase II

To allow motion of the plasma we introduce the equation

~=c1(jXB— together with Ohm’s law E = o~j,to describe the field diffusion into a rigid conductor. The equations are written in cylindrical coordinates based on the major axis of the torus. A simple explicit scheme was used on a rectangular mesh. The model accurately described the field penetration into a cylinder [1] as shown in fig. 1. The field penetration into a rigid conductor in the JET geometry shown below was then investigated and the result is presented in fig. 2.

0.2

Vp)+V-c2Vu.

The ultimate objective is to arrange c1 so that we smoothly obtain the displacements required to make / X B = Vp. The value of c1 is chosen to be as small as possible without impeding the resistive flows. The viscosity term is used to damp out any obtrusive oscilla. tions. The pressure is determined from the simplified energy equation including Ohmic heating and the calculated velocities are used in the extended Ohm’s law E + u X B = 77/. To solve this new set of equations we now introduce the Lax—Wendroff scheme which will form the basic numerical method as the code is developed into its final form. The structure of this scheme is such that computations take place on two independent but interlocking meshes, thus one half of the compu-

MF. Turner, J.A. Wesson

/ 3D MHD simulation of JET

Fig. 3. Initialisation of poloidal magnetic field, current and pressure for phase II test case (20 X 32 mesh). Bq~

345

0= 1 T, E40 =

total current (1)

=

Picture

Parameter depicted

2 3

poloidal toroidal current current poloidal magnetic field pressure

4

1

3.7 MA. Maximum value 2 5.6 0.6 KA/m2 MA/rn 0.59 T 0.15 MPa

tations may be regarded as redundant. We prefer the retain computations from both meshes so that at the boundary the full rectangular mesh is available for the calculation of boundary conditions. The results

from the two meshes are then effectively linked together at the boundary. Using this scheme on a 20 X 32 mesh calculations were carried out for 2D axisymmetric JET geometry.

V/rn,

ME. Turner, J.A. Wesson / 3D MHD simulation of JET

346

~ ~

J~\~

\I\\~C\~

\ \\

2 1

\~~

~j~P

I

II

\\

~

lI~~

Fig. 5. Poloidal magnetic field, poloidal current, pressure and velocity in phase II test case after 3 s. Toroidal current at this time is shown in fig. 4. The circles in velocity plot indicate values of toroidal velocity. Complete circles are negative values, split circles are positive values. Picture

Parameter depicted

Maximum value

1 2 3 4

velocity poloidal current poloidal magnetic field pressure

0.56 0.34 0.87 0.15

m/s(POL), 0.05 m/s(TOR) 2 KA/m T MPa

Fig. 4. Time development of toroidal current in phase II

test case

(see page 346).

Picture

Time (s)

Maximum value (MA/rn2)

Total current (MA)

1 2

0.0 0.02

0.63 0.86

3.7 4.1

3 4

0.06 0.40

0.87 0.81

4.2 4.6

5

1.50

0.74

4.6

6

3.0

0.73

4.6

ME. Turner, J.A. Wesson / 3D MHD simulation of JET

348

For a test case the system was initialised with a Solovev [2] equilibrium having an almost flat current profile and zero velocity. This initialisation is depicted in fig. 3 which shows the toroidal and poloidal currents, the pressure distribution and the poloidal magnetic field. The shaped resistivity, 77(R), required to produce this equilibrium was then replaced by a uniform resistivity and the system allowed to move to a new equilibrium state. For this run the coefficient c1

plot the value at the plasma edge but only at the first internal mesh point. It is clear that by three seconds the current distribution has substantially adjusted to the change in the resistivity profile, but the final equilibrium has not been reached. The state of the system after 3 s is shown in fig. 5. Phase III The full energy and continuity equations are now introduced leading to the following complete set of equations for a single fluid:

in the velocity eauation has value 1.3 X 106 and there is no viscosity. The evolution of the current density profile is shown in fig. 4. The apparent irregularity at the boundary in the plots is due to the graphical method which does not

-~—

1

.~

.

I .

—v x E,

.\

_~‘___~•••

/

/

=

I

.~

I

3

Fig. 6a. Evolution of temperature over first second for case 1 of phase III (see table on page 350).

--

4

M. F. Turner, J.A. Wesson / 3D MHD simulation of JET

E

=77/—uXB, =

V

X

B,

at

=c (jXB—Vp)+ V-c2Vu.

~

=_v.(~pu)_pV.U— V-q+Q—L, =

— V~(nu)

Q

= =

heating = 77/2 energy loss (e.g. radiation), at present given 2)2pwhereFisa function of distance in from theformcL(1 —F the edge of the plasma,

n electrondensity. For the first runs with this model we initialised the system with zero current, small initial temperature and uniform electron density. With a background toroidal field (B~cx 1/R) everywhere, a fixed electric field (E 4 cx l/R) is imposed at the boundary and the

‘

=

=—KVT,

L

l~o

at

q

349

2nT

p

I\Afi~i.’

2

/~,

/

1(1 (1

3

Fig. 6b. Evolution of toroidal current over first second for case 1 of phase III (see table on page 350).

L

4

ME. Turner, J.A. Wesson / 3D MHD simulation of JET

350

Co

i:0~’~

~

~~=::~gc~

i~ ~

-

I

~

~

.

~

..

.

..~.e, a

o..-’-~

.

a

a

~ ~

Fig. 6c. Evolution of velocity over first second for case 1 of phase III. Picture (6a, b, c)

Maximum values

Total current (MA)

time (s)

temperature (key)

tor. current density 2) (MAim

pol. velocity (m/s)

tor. velocity

1

0.16

0.20

0.47

0.74

0.02

1.15

2 3 4

0.30 0.5 1.0

0.26 0.43 0.74

0.66 0.99 1.60

1.27 1.46 2.36

0.04 0.10 1.00

2.0 3.0 4.9

(mis)

M.F. Turner, J.A. Wesson / 3D MHD simulation of JET

time evolution of the system studied. The conditions on the boundary of the plasma are:

351

toroidal voltage (E~,)= E~1~0R0 /R

TA, and a uniform pressure consistent with this ternperature and the initial electron density. For the first runs parameter values have been chosen as follows:

toroidal field (B4) = B40R0 /R

B40

temperature (7) given small value (TA) poloidal field (Br) parallel to boundary to simulate

3 initial electron density = 3.0 X 1019 m thermal conductivity (K) = 0.7 X 1020 m~ s~

=

1 T,

TA

40 eV,

equilibrium control. .

.

electrical resistivity = 2.8 X 1 08/(TkeV)3”2 ~2m

The system is initialised with uniform temperature

1

I

/1

I

I I

2

> I

~

.-~1-~~_~>

/1/

I

i~

C

~

I

Fig. 7a. Evolution of temperature over first five seconds for case 2 of phase III (see table on page 353).

ME. Turner, J.A. Wesson / 3D MHD simulation of JET

352

5. Results

energy loss coefficient, CL = 0

Several preliminary runs have been carried out for various values of applied electric field (E40), energy loss coefficient (CL) and velocity coefficients c1 and c2. We present two examples here. Case 1 velocity equation coefficient, c1 viscosity coefficient, c 2 = 2.0 ,

3

=

1.3

x io—

applied electric field atR ~=R0,E40= 0.2 V m’ Evolution of the principal variables over the first second is depicted in fig. 6. We see that the temperature rises at the edge initially before building up in the centre as a result of energy transport. The toroidal current howeverdevelops remainsinasthe a skin current. Circulation of the velocity upper and lower halves of the plasma section and the general level of the velocity rises steadily. In the present run this gradual

‘a

Fig. 7b. Evolution of toroidal current over first five seconds for case 2 of phase 111 (see table on page 353).

M.F. Turner, J.A. Wesson / 3D MHD simulation of JET

OS 0

~ ~

~ ~

~.

353

-

-

-~

LI

Fig. 7c. EvolutionTime of velocity over first five seconds of case 2 of phase III. Picture Maximum values (7a, b, c) (s) temperature tor. current pol. velocity (keV)

density 2) (MA/rn

tor. velocity

(m/s)

(mis)

Total current (MA)

1

0.20

0.07

0.09

2 3 4 5

0.90 1.15 1.35 1.65

0.15 0.18 0.22 0.27

0.17 0.24 0.24 0.34

(1) 0.83 (2) 0.86 (3) 0.83

(1) 0.03 (2) 0.03 (3) 0.04

0.3 0.8 1.1 1.3 1.5

6

5.20

0.69

0.64

(4) 0.76

(4) 0.84

2.8

ME. Turner, J.A. Wesson / 3D MHD simulation of JET

354

rise appears eventually to lead to the development of numerical instabilities, Case 2 Here we reduce the coefficient in the velocity equation to 40% of its case 1 value and we also considerably reduce the viscosity. Some energy loss is introduced and the applied electric field is halved, We now have the following set of coefficients Ci

CL

— —

4

5.0 X 10 5.0 s~

— ,

,

C

2

—0.15

E

=

0.1 V m’

0

The evolution over 5 s is depicted in fig. 7. The aim here is to move towards a stationary equilibrium with a peaked current profile. As the run progresses is gradually reduced in order to encourage central peaking of the current. We see that introduction of energy losses nearer the plasma edge causes central peaking of the temperature to occur more readily. This leads to reduced resistivity at the centre and the skin effect is much less marked than that of case 1. However around t = 0.9 s the current density begins to adopt a form having peaks at the top and bottom of the plasma section and is reduced at the inner edge. This is in spite of the fact that the toroidal electric field is held proportional to 1 /R at the edge of the plasma. Looking at the form of the velocity field at this time we see that the plasma is moving across the poloidal magnetic field from the inner to the outer edge but is approximately parallel to it at the top and bottom. Thus there is a large u X B component in the mid-section which is strongest at the inner edge

because of the higher magnetic field there. This u X B component substantially reduces the effective electric field at the inner edge but has little effect at the top and bottom allowing the toroidal current density to peak there. Later, at about 1.3 s, circulations begin to appear in the velocity, so that now the velocity and magnetic field vectors become more nearly parallel at the inner edge and the current rises there giving a smoother shape to the current density profile. At 5 s the current density profile is still gradually peaking and the velocity flows have not yet fully evolved into the expected form for the final equiibrium.

6. Summary Substantial progress has been made in the development of the computer code for the simulation of the JET experiment. Accurate calculation of magnetic field diffusion has been demonstrated and the development of equilibria has been followed over several seconds. Further investigations will be made into the operation of the code in this form before moving onto the next phase when computations will be extended into the toroidal dimension.

References [1] J.A. Wesson, Plasma Phys. 4 (1962) 175. [2] L.S. Solovev, Nucl. Fusion 25 (1968) 400.