Numerical Magnetohydrodynamics for High-Beta Plasmas

Numerical Magnetohydrodynamics for High-Beta Plasmas JEREMIAH U . BRACKBILL UNIVERSITY OF CALIFORNIA LOS ALAMOS SCIENTIFIC LABORATORY LOS ALAMOS, N E W MEXICO

I. Introduction . . . . . . . . . . . . 1 II. N u m e r i c a l M e t h o d s . . . . . . . . . . . 3 A . A Description of Eulerian, Lagrangian, and Generalized C o m p u t a t i o n Meshes . . . . . . . . . . . . 4 B . A Survey o f N u m e r i c a l M e t h o d s and Applications 8 III. T h e C o m p u t a t i o n o f Convective Transport 10 A . Properties o f Approximations to the Convective Derivative . . . 1 0 Β. Nonlinearly Stable Approximations to Convective Transport . . . 1 3 IV. A Generalized M e s h M e t h o d for Μ H D 17 A . Difference Equations for the Lagrangian Phase of a Generalized M e s h Calculation . . . . . . . . . . . . 1 8 Β. T h e R e z o n e Phase of a Generalized M e s h Calculation . . . . 26 V. Applications 29 A . A Sharp Boundary Calculation of the Theta Pinch . . . . . 29 Β. T h e R o t a t i n g Theta Pinch 31 C. T h e Internal K i n k M o d e Instability 34 VI. Conclusions . . . . . . . . . . . . 38 Appendix. . . . . . . . . . . . . 38 References . . . . . . . . . . . . 39

I. Introduction

MAGNETOHYDRODYNAMICS ( M H D ) is t h e least s o p h i s t i c a t e d m o d e l o f a m a g ­ netically confined p l a s m a which describes the interaction between a magnetic field a n d a p l a s m a self-consistently. T h e M H D m o d e l t r e a t s t h e p l a s m a a s t h o u g h it w e r e a c h a r g e - n e u t r a l fluid i n l o c a l t h e r m o d y n a m i c e q u i l i b r i u m , a n d t h u s n e g l e c t s all b u t a s m a l l p a r t o f t h e p h y s i c s o f p l a s m a s . T h a t s m a l l part, however, describes the transfer of m o m e n t u m a n d energy between the p l a s m a a n d t h e m a g n e t i c field, a n d t h a t is sufficient t o d e s c r i b e t h e effect o f t h e g e o m e t r y o f t h e field o n t h e g r o s s e q u i l i b r i u m a n d s t a b i l i t y o f a h i g h - b e t a plasma. E v e n in t h e i r s i m p l e s t f o r m , t h e M H D e q u a t i o n s c o m p r i s e a c o u p l e d 1

2

JEREMIAH U. BRACKBILL

s y s t e m o f n o n l i n e a r , p a r t i a l differential e q u a t i o n s w h i c h a r e difficult, if n o t i m p o s s i b l e , t o solve a n a l y t i c a l l y . S o m e t i m e s t h e M H D e q u a t i o n s c a n b e r e d u c e d , b e c a u s e o f s y m m e t r i e s o r b e c a u s e o f o r d e r i n g , t o a single o r d i n a r y differential e q u a t i o n w h o s e s o l u t i o n d e s c r i b e s t h e l i n e a r s t a b i l i t y o f a p l a s m a in e q u i l i b r i u m w i t h a m a g n e t i c field. M o r e often, t h e e q u a t i o n s c a n n o t b e r e d u c e d , a n d t h e M H D e q u a t i o n s in t w o , a n d o f t e n t h r e e , d i m e n s i o n s m u s t b e s o l v e d . I n s u c h c a s e s , t h e r e is n o o t h e r w a y t o l e a r n w h a t t h e M H D m o d e l p r e d i c t s a b o u t t h e stability of a p a r t i c u l a r e x p e r i m e n t t h a n b y n u m e r i c a l l y solving the equations. O n e h a s o n l y t o e x a m i n e t h e list o f p r o b l e m s o f c u r r e n t i n t e r e s t c o m p i l e d b y B o d i n (1972) t o see t h a t m a n y o f t h e m will o n l y yield t o n u m e r i c a l cal­ culation. T w o of the p r o b l e m s are concerned with axisymmetric a n d n o n axisymmetric toroidal confinement. Several two-dimensional numerical c o m p u t a t i o n s of the axisymmetric case have already been performed (Lui a n d C h u , 1 9 7 4 ; H o f m a n n , 1974), a n d it is c l e a r t h a t fully t h r e e - d i m e n s i o n a l c o m p u t a t i o n s will a l s o b e n e c e s s a r y b e f o r e instabilities in t h e s e c o n f i n e m e n t geometries can be simulated. Three-dimensional simulations are an ambitious u n d e r t a k i n g w h i c h t h e c u r r e n t i n t e r e s t in c o n t r o l l e d t h e r m o n u c l e a r r e s e a r c h c a n s u p p o r t if t h e r e s u l t s w a r r a n t it. T o e v a l u a t e w h e t h e r n u m e r i c a l c o m p u t a t i o n s c a n p l a y a useful r o l e in t h e s o l u t i o n o f p r o b l e m s in C T R , o n e s h o u l d e x a m i n e n o t o n l y t h e o p p o r t u n i t i e s b u t a l s o t h e c o n s t r a i n t s o n n u m e r i c a l w o r k . O n e c o n s t r a i n t is t h a t n o w a n d for t h e n e x t several y e a r s , t h e c o m p u t a t i o n a l p o w e r a t t h e d i s p o s a l o f a r e s e a r c h e r u s i n g n u m e r i c a l m e t h o d s is b a r e l y sufficient t o s t o r e a n d p r o c e s s e n o u g h i n f o r m a t i o n t o d e s c r i b e t h e c u r r e n t l y i n t e r e s t i n g m a g n e t i c confine­ m e n t g e o m e t r i e s . W h e n o n e c o n s i d e r s t h a t a c o m p a r a b l e a m o u n t o f infor­ m a t i o n h a s b e e n u s e d for n u m e r i c a l c o m p u t a t i o n s i n t w o d i m e n s i o n s , o n e c a n see t h a t m o r e a c c u r a c y is r e q u i r e d o f n u m e r i c a l a p p r o x i m a t i o n s in t h r e e d i m e n s i o n s t h a n in t w o if useful r e s u l t s for t i m e - d e p e n d e n t p r o b l e m s a r e t o be obtained. A review o f t h e a c c u r a c y o f m a n y n u m e r i c a l s o l u t i o n a l g o r i t h m s in c u r r e n t u s e reveals t h a t t h e c o m p u t a t i o n o f c o n v e c t i v e t r a n s p o r t is significantly less a c c u r a t e t h a n t h e c o m p u t a t i o n o f t h e o t h e r t e r m s in t h e e q u a t i o n s d e s c r i b i n g m a g n e t o h y d r o d y n a m i c flow. F u r t h e r m o r e it is c l e a r f r o m t h e l i t e r a t u r e o f n u m e r i c a l m e t h o d s for h y d r o d y n a m i c flows t h a t c o n v e c t i v e t r a n s p o r t h a s long been a recognized problem. W h y the a p p r o x i m a t i o n of the convective d e r i v a t i v e s h o u l d c a u s e p r o b l e m s is t h e s u b j e c t o f a d i s c u s s i o n i n S e c t i o n I I I w h i c h calls t h e r e a d e r ' s a t t e n t i o n t o t h e r e l e v a n t l i t e r a t u r e , a n d r e v i e w s t h e current w o r k o n i m p r o v i n g the accuracy a n d the stability of numerical a p p r o x i m a t i o n s t o this term. Sections IV a n d V present a numerical m e t h o d which avoids, as m u c h as possible, having t o a p p r o x i m a t e the convective derivative.

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

3

T h e a p p r o x i m a t i o n o f c o n v e c t i v e t r a n s p o r t is a m a j o r t h e m e i n t h i s a r t i c l e , b o t h b e c a u s e t h e w o r k d e s c r i b e d in S e c t i o n s I I I - V h a s significantly i m p r o v e d t h e a c c u r a c y o f n u m e r i c a l s o l u t i o n s o f t h e M H D e q u a t i o n s , a n d b e c a u s e it is a f u n d a m e n t a l yet n o t generally recognized p r o b l e m . T h e a r t i c l e b e g i n s w i t h a b r i e f i n t r o d u c t i o n t o n u m e r i c a l m e t h o d s for m a g n e t o h y d r o d y n a m i c s , including a survey of published m e t h o d s a n d applications.

II. Numerical M e t h o d s There has evolved a s t a n d a r d a p p r o a c h to the numerical solution of the p a r t i a l differential e q u a t i o n s d e s c r i b i n g t i m e - d e p e n d e n t m a g n e t o h y d r o d y n a m i c flow. Finite-difference a p p r o x i m a t i o n s t o t h e differential e q u a t i o n s are solved o n a c o m p u t a t i o n m e s h at a sequence of time-steps. O t h e r m e t h o d s , s u c h a s p a r t i c l e d e s c r i p t i o n s o f t h e fluid, h a v e b e e n u s e d t o s i m u l a t e t h e p l a s m a f o c u s e x p e r i m e n t ( B u t l e r et al., 1 9 6 9 ; R o b e r t s a n d P o t t e r , 1968) a n d t h e t h e t a p i n c h e n d loss p r o b l e m ( T u c k , 1968). H o w e v e r , b e c a u s e o f t h e c u r r e n t i n t e r e s t i n t h e g r o s s s t a b i l i t y o f v a r i o u s c o n f i n e m e n t g e o m e t r i e s ( B o d i n , 1972) a n d in t h e t r a n s p o r t p r o c e s s e s i n p l a s m a s , finite-difference a p p r o x i m a t i o n s a r e n o w u s e d a l m o s t exclusively. T h e r e is a l s o a s t a n d a r d a p p r o a c h t o a n o t h e r p r o b l e m . I n p i n c h d i s c h a r g e s , for e x a m p l e , t h e b u l k o f t h e p l a s m a is s w e p t i n w a r d b y a n i n c o m i n g m a g n e t i c field l e a v i n g b e h i n d a l o w d e n s i t y p l a s m a i m m e r s e d in a s t r o n g m a g n e t i c field. I n t h e l o w d e n s i t y p l a s m a , t h e c h a r a c t e r i s t i c s i g n a l s p e e d , t h e A l f v e n s p e e d , is v e r y l a r g e a n d c o n s e q u e n t l y t h e C o u r a n t c o n d i t i o n o n t h e t i m e - s t e p ( R i c h t m y e r a n d M o r t o n , 1967) is very r e s t r i c t i v e . S e v e r a l s o l u t i o n s h a v e b e e n offered for t h i s p r o b l e m . B o r i s (1970) offered t h e s u g g e s t i o n t h a t a relativistic m a s s c o r r e c t i o n b e a p p l i e d in w h i c h a n artificially l o w m a x i m u m s i g n a l s p e e d is s u b s t i t u t e d for t h e s p e e d o f light. By c h o o s i n g t h e a p p r o p r i a t e v a l u e for t h e signal s p e e d , t h e d e s i r e d m i n i m u m t i m e - s t e p c a n b e set. L u i a n d C h u (1974) address the same p r o b l e m in a pinch discharge calculation by setting the m i n i m u m a l l o w a b l e d e n s i t y in t h e m e s h t o 1 5 % o f t h e initial a v e r a g e d e n s i t y . T h e r e is g e n e r a l a g r e e m e n t n o w t h a t a n i m p l i c i t f o r m u l a t i o n o f t h e e q u a t i o n s o f m o t i o n s offers t h e b e s t s o l u t i o n . W i t h a n i m p l i c i t f o r m u l a t i o n , t h e t i m e - s t e p is n o l o n g e r d e t e r m i n e d b y t h e m a x i m u m s i g n a l s p e e d in t h e m e s h b e c a u s e t h e equations are unconditionally stable. Implicit equations require the solution o f s y s t e m s o f s i m u l t a n e o u s e q u a t i o n s for w h i c h v a r i o u s i t e r a t i v e m e t h o d s a r e used, including time-step splitting, alternating direction implicit, a n d even successive o v e r r e l a x a t i o n m e t h o d s ( R i c h t m y e r a n d M o r t o n , 1967). A l l o f these m e t h o d s are evidently economical. F i n a l l y , t h e r e is t h e p r o b l e m o f t h e a c c u r a t e c o m p u t a t i o n o f c o n v e c t i v e

4

JEREMIAH U. BRACKBILL

t r a n s p o r t . A s will b e d i s c u s s e d i n S e c t i o n I I I , t h e r e is s o m e e v i d e n c e t h a t t h e accuracy of m a n y a p p r o x i m a t i o n s t o the convective derivative m u s t be deliber­ a t e l y r e d u c e d t o a v o i d n o n l i n e a r i n s t a b i l i t i e s . S t a b l e s c h e m e s a r e often o b ­ t a i n e d b y a d d i n g c o n s i d e r a b l e diifusion t o t h e difference e q u a t i o n s m a k i n g m e a n i n g f u l t h r e e - d i m e n s i o n a l c a l c u l a t i o n s m o r e difficult t o d o w i t h p r e s e n t day computers. M o r e accurate approximations which are nonlinearly stable a r e d e s c r i b e d i n S e c t i o n I I I . I n S e c t i o n s I V a n d V , a s e c o n d a p p r o a c h is discussed in which the convective transport, a n d therefore the need to a d d s t a b i l i z i n g diifusion, is m u c h r e d u c e d b y u s i n g t h e g e n e r a l i z e d m e s h d e s c r i b e d in t h e f o l l o w i n g s e c t i o n . ( I n a d d i t i o n , t h e e q u a t i o n s f o r E u l e r i a n , L a g r a n g i a n , a n d g e n e r a l i z e d m e s h e s a r e g i v e n , a s well a s a s u r v e y o f M H D c o m p u t a t i o n s using each kind of mesh, in order that the discussion in the remaining sections m a y proceed in a reasonably informed manner.)

A. A DESCRIPTION OF EULERIAN, LAGRANGIAN, AND GENERALIZED COMPUTATION MESHES T h e first s t e p i n f o r m u l a t i n g a n u m e r i c a l a l g o r i t h m for s o l v i n g t h e M H D e q u a t i o n s is t o c h o o s e i n a w a y t o r e p r e s e n t t h e p l a s m a . T y p i c a l l y , t h e p l a s m a is r e p r e s e n t e d b y a finite n u m b e r o f n u m b e r s w h i c h give t h e p o s i t i o n , velocity, d e n s i t y , t e m p e r a t u r e , a n d m a g n e t i c field i n t e n s i t y o v e r t h e p h y s i c a l d o m a i n . T h e s e n u m b e r s a r e s t o r e d a t a n a r r a y o f p o i n t s , called m e s h p o i n t s , c o m p r i s i n g a c o m p u t a t i o n m e s h . E a c h m e s h p o i n t is a s s o c i a t e d w i t h a cell w h i c h is u s e d a s a c o n t r o l v o l u m e i n c o n s t r u c t i n g c o n s e r v a t i v e difference e q u a t i o n s . A c o m ­ p u t a t i o n m e s h m a y b e E u l e r i a n , L a g r a n g i a n , o r g e n e r a l i z e d . O n e k i n d is distinguished from a n o t h e r by the relative m o t i o n of the m e s h points a n d the fluid, a n d b y t h e c o m p l e x i t y o f t h e difference e q u a t i o n s a p p r o x i m a t i n g s p a t i a l derivatives. A n E u l e r i a n m e s h is s t a t i o n a r y w i t h r e s p e c t t o t h e l a b o r a t o r y f r a m e , a n d a L a g r a n g i a n m e s h is s t a t i o n a r y w i t h r e s p e c t t o t h e p l a s m a . A g e n e r a l i z e d mesh m a y be Eulerian, Lagrangian, or neither: the m o t i o n of the grid points w i t h r e s p e c t t o b o t h t h e l a b o r a t o r y f r a m e a n d t h e p l a s m a is a r b i t r a r y . T h e grid m o t i o n determines h o w m u c h convective transport m u s t be c o m p u t e d . I n the following sections, brief descriptions are given of the Eulerian, L a g r a n g i a n , a n d g e n e r a l i z e d m e s h e s w i t h t h e a p p r o p r i a t e e q u a t i o n s for e a c h m e s h . A f t e r t h e d e s c r i p t i o n s , a b r i e f s u r v e y is p r e s e n t e d o f t h e c a l c u l a t i o n s employing each kind of mesh. 1. The Eulerian

Computation

Mesh

E u l e r i a n m e s h e s a r e s i m p l e a n d , in m a n y c a s e s , e c o n o m i c a l t o u s e . Since t h e m e s h is s t a t i o n a r y , t h e z o n i n g m a y b e c h o s e n t o simplify t h e n u m e r i c a l

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

5

a p p r o x i m a t i o n o f s p a t i a l d e r i v a t i v e s . W h e n t h e m e s h is o r t h o g o n a l , a r a p i d s o l u t i o n a l g o r i t h m s u c h a s t i m e - s t e p s p l i t t i n g ( R i c h t m y e r a n d M o r t o n , 1967) c a n be used. M o s t i m p o r t a n t , the c o n s e r v a t i o n f o r m of the ideal M H D e q u a t i o n s is differenced, a s s u r i n g t h e r i g o r o u s c o n s e r v a t i o n o f m a s s , m o m e n ­ t u m , m a g n e t i c flux, a n d e n e r g y . I n c o n s e r v a t i o n f o r m , t h e c o n t i n u i t y , i n d u c ­ tion, m o m e n t u m , a n d energy equations are written dp/dt

+ V · (pu) = 0,

(1)

+ V χ (u χ Β ) = 0,

(2)

+ V · (pun - Q ) = 0,

(3)

+ V · (pen - Q · u) = 0,

(4)

dB/dt dpu/dt dpe/dt

w h e r e ρ is t h e m a s s d e n s i t y o f t h e p l a s m a , u t h e p l a s m a v e l o c i t y , Β t h e m a g n e t i c field i n t e n s i t y , Q t h e stress t e n s o r , a n d e t h e specific t o t a l p l a s m a e n e r g y . T h e stress t e n s o r , Q , a n d p l a s m a e n e r g y , e, a r e g i v e n in m k s u n i t s b y Q = (1/μ)ΒΒ - £[(1/μ)Β · Β + Ρ ] II

(5)

and e = i ( U . u ) + i ( B . B ) / / z p + /,

(6)

w h e r e μ is t h e p e r m e a b i l i t y , / is t h e specific i n t e r n a l e n e r g y o f t h e p l a s m a , a n d Ρ is t h e p l a s m a p r e s s u r e e q u a l t o a f u n c t i o n o f ρ a n d /. W h e n t h e r e a r e l a r g e g r a d i e n t s i n v e l o c i t y , v i s c o s i t y m u s t b e a d d e d t o t h e stress t e n s o r ( Z e P d o v i c h a n d R a i z e r , 1967). A v i s c o s i t y f o r m u l a t i o n is g i v e n b y L a n d a u a n d Lifschitz (1958). O t h e r t r a n s p o r t t e r m s , s u c h a s m a s s diffusion, t h e r m a l c o n d u c t i v i t y , a n d resistive diffusion m u s t s o m e t i m e s b e a d d e d . E x t e n s i v e d i s c u s s i o n s o f t h e n u m e r i c a l p r o b l e m s e n c o u n t e r e d in a d d i n g t h e s e t e r m s a r e c o n t a i n e d in r e v i e w a r t i c l e s b y R o b e r t s a n d P o t t e r (1968) a n d K i l l e e n (1972). 2. Lagrangian

Computation

Meshes

A L a g r a n g i a n m e s h m a y b e b u i l t u p f r o m cells o f a r b i t r a r y s h a p e a n d size. T h e cells i n t w o d i m e n s i o n s , h o w e v e r , a r e u s u a l l y a r b i t r a r y q u a d r i l a t e r a l s b e c a u s e s u c h cells h a v e s i m p l e logical r e l a t i o n s h i p s w i t h o n e a n o t h e r . T r i ­ a n g u l a r cells h a v e s o m e a d v a n t a g e s ; t h e y a r e a l w a y s c o n v e x , a n d i n t e r p o l a t i o n s in t h e i r i n t e r i o r a r e l i n e a r . H o w e v e r , w i t h t r i a n g u l a r z o n e s , a t a b l e o f n e a r e s t n e i g h b o r s m u s t b e k e p t b e c a u s e t h e r e is n o l o n g e r a s i m p l e r e l a t i o n s h i p b e t w e e n t h e p h y s i c a l a n d logical p o s i t i o n o f a cell w i t h r e s p e c t t o its n e i g h b o r s .

6

JEREMIAH U. BRACKBILL

D e r i v i n g difference e q u a t i o n s for a L a g r a n g i a n m e s h is n o t a s a u t o m a t i c a s it is for a r e g u l a r E u l e r i a n m e s h . T h e m e t h o d o f d e r i v a t i o n r e q u i r e s a c o o r d i n a t e t r a n s f o r m a t i o n f r o m t h e n o n r e c t i l i n e a r p h y s i c a l s p a c e t o a recti­ l i n e a r n a t u r a l o r logical s p a c e b e f o r e a n a l o g s t o o r d i n a r y different e q u a t i o n s c a n b e c o n s t r u c t e d . A d e r i v a t i o n o f g e n e r a l i z e d difference e q u a t i o n s for t h r e e d i m e n s i o n a l M H D flow is g i v e n i n S e c t i o n I V . T h e L a g r a n g i a n e q u a t i o n s for i d e a l M H D , c o r r e s p o n d i n g t o t h e E u l e r i a n equations, Eqs. (l)-(4), are written dp/dt

+ p ( V - u ) = 0,

(7)

έ/Β/Λ + Β ( ν · π ) - ( Β · V ) u = 0,

(8)

ρ du/dt

- V · Q = 0,

(9)

ρ dijdt + P ( V · u) = 0.

(10)

and

T h e m o s t o b v i o u s difference b e t w e e n t h e s e e q u a t i o n s a n d t h e E u l e r i a n e q u a t i o n s is t h e a b s e n c e o f c o n v e c t i v e t r a n s p o r t t e r m s . T h i s is c h a r a c t e r i s t i c o f t h e L a g r a n g i a n e q u a t i o n s , in t h a t t h e y d e s c r i b e t h e e v o l u t i o n o f a p a r t i c u l a r e l e m e n t o f t h e fluid, r a t h e r t h a n t h e t i m e v a r i a t i o n o f fluid v a r i a b l e s a t a p a r t i c u l a r p o i n t in s p a c e . 3. The Generalized

Computation

Mesh

I n a generalized mesh, the relative m o t i o n of the p l a s m a a n d the c o m ­ p u t a t i o n m e s h is c o m p l e t e l y a r b i t r a r y . T h a t is, it is p o s s i b l e t o specify t h e velocity of the c o m p u t a t i o n m e s h separately from the velocity of t h e p l a s m a . T h e e q u a t i o n s for a g e n e r a l i z e d m e s h a r e d e r i v e d e i t h e r b y a d d i n g t o a L a g r a n g i a n c a l c u l a t i o n a r e z o n e p h a s e in w h i c h t h e c o n v e c t i v e t r a n s p o r t d u e t o r e l a t i v e m o t i o n b e t w e e n t h e p l a s m a a n d t h e m e s h is c o m p u t e d , o r b y r e w r i t i n g t h e E u l e r i a n e q u a t i o n s , E q s . ( l ) - ( 4 ) , in a m o v i n g f r a m e . There are m a n y advantages to be gained from a generalized mesh. Con­ vective t r a n s p o r t m a y b e r e d u c e d , in m a n y p r o b l e m s o f i n t e r e s t , b y t r a n s ­ forming a w a y the bulk m o t i o n of the p l a s m a as suggested by R o b e r t s a n d P o t t e r (1968). C o n t a c t s u r f a c e s , a s b e t w e e n a p l a s m a a n d a v a c u u m , m a y b e r e s o l v e d easily b y a L a g r a n g i a n i n t e r f a c e w i t h i n a g e n e r a l i z e d m e s h . F o r o t h e r a p p l i c a t i o n s , t h e m e s h p o i n t s m a y b e c o n s t r a i n e d t o lie o n flux surfaces t o allow the accurate calculation of anisotropic thermal conduction (Hertweck a n d S c h n e i d e r , 1970). F i n a l l y , t h e m e s h p o i n t s c a n b e c o n c e n t r a t e d in a r e g i o n o f t h e m e s h w h e r e r e s o l u t i o n is r e q u i r e d , e v e n a s t h a t r e g i o n m o v e s t h r o u g h

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS s p a c e ( H i r t et al, flexibility

7

1974). I n s u m m a r y , a g e n e r a l i z e d m e s h offers e n o r m o u s

in d e s i g n i n g t h e c a l c u l a t i o n t o fit t h e a p p l i c a t i o n .

A s m e n t i o n e d a b o v e , t h e e q u a t i o n s for a g e n e r a l i z e d m e s h c a n b e o b t a i n e d by transforming the Eulerian equations to a m o v i n g coordinate system. In a frame m o v i n g with velocity u', the e q u a t i o n s c o r r e s p o n d i n g t o E q s . ( l ) - ( 4 ) are written dp/dt

+ V · [p(u - a')] + p(V · u') = 0,

(11)

dB/dt

+ V · (u - u ' ) B - V · B u + B ( V · u') = 0,

(12)

dpu/dt

+ V · [pu(u ~ u') - Q ] + pu(V · u') = 0,

(13)

and dpe/dt

+ V · lpe(u

- u ) - Q · u] + pe(W · u') = 0,

(14)

w h e r e , i n t h e i n d u c t i o n e q u a t i o n , t e r m s o f t h e f o r m V · ( A B ) m e a n V is s c a l a r m u l t i p l i e d w i t h A b u t differentiates b o t h A a n d B . W h e n u ' is z e r o , t h e s e equations reduce t o the conservation equations, Eqs. (l)-(4), a n d when u' a n d u are equal, they reduce t o the Lagrangian equations, Eqs. (7)-(10). Alternatively, the equations can be separated into Lagrangian a n d con­ vective t r a n s p o r t p h a s e s a s p r o p o s e d b y H i r t et al. (1974). I n t h i s f o r m u l a t i o n , t h e t i m e v a r i a t i o n i n t h e m o v i n g f r a m e is c o m p u t e d i n t w o s t e p s . I n a L a g r a n g i a n p h a s e , E q s . ( 7 ) - ( 1 0 ) a r e s o l v e d , a n d in a c o n v e c t i o n p h a s e , t h e c o n v e c t i v e t r a n s p o r t is c o m p u t e d f r o m t h e i n t e g r a l e q u a t i o n s

U dv P

d r dt Jv

=ί

[A ·

(u -v O p ]

du — dt

,

MV)

dtJy

pudV

=

r Jv

rp

dV-

d r

— ί and

BdV

=

pedV

=4 ί

ds,

(15)

J [Ä . (u - u ) pu] ds,

(16)

f

(17)

;

S(K)

s ( F[ Ä ) .

(u-u')B]^,

dt Jv

ΙΓ dtjy

dtiy

pedV-

ί

My)

[A ·

(u — 'u)p^] ds,

(18)

w h e r e Κ is a c o n t r o l v o l u m e , a n d s(V) is its b o u n d a r y w i t h o u t w a r d d i r e c t e d u n i t n o r m a l Ä. B o t h t h e c o n s e r v a t i o n f o r m a n d t h e t w o - p h a s e f o r m p e r m i t

8

JEREMIAH U. BRACKBILL

arbitrary relative m o t i o n between the p l a s m a a n d the mesh, a n d b o t h reduce t o t h e L a g r a n g i a n f o r m w h e n t h e velocity o f t h e m e s h a n d t h e velocity o f t h e fluid a r e e q u a l . T h e t w o - p h a s e f o r m u l a t i o n d o e s p e r m i t i m p l i c i t e q u a t i o n s o f m o t i o n a n d explicit c o n v e c t i o n t e r m s . Since t h e m o r e s t r i n g e n t s t a b i l i t y condition m u s t be m e t by the equations of m o t i o n , gains in r u n n i n g speed a r e m a d e b y m a k i n g t h e m i m p l i c i t e v e n w h e n t h e c o n v e c t i o n t e r m s a r e explicit.

B. A SURVEY OF NUMERICAL METHODS AND APPLICATIONS 1. Eulerian

Computations

Since t h e review b y R o b e r t s a n d P o t t e r (1968), t h e r e h a v e b e e n m a n y a p p l i c a t i o n s of E u l e r i a n c o m p u t a t i o n s t o M H D p r o b l e m s . F r e e m a n a n d L a n e (1969) r e p o r t e d a n M H D m e t h o d for t w o - d i m e n s i o n a l a x i s y m m e t r i c flows u s i n g a n explicit, L a x - W e n d r o f f t i m e a d v a n c e m e n t algorithm ( R i c h t m y e r , 1963), w i t h a L a p i d u s diifusion t e r m ( L a p i d u s , 1967) a d d e d t o e a c h o f t h e difference e q u a t i o n s c o r r e s p o n d i n g t o E q s . ( l ) - ( 4 ) . T h i s m e t h o d w a s later applied t o the simulation of the interaction of a p l a s m o i d with a n a x i s y m m e t r i c m a g n e t i c field ( F r e e m a n , 1971). R o b e r t s a n d B o r i s (1969) r e p o r t e d a m e t h o d for t h r e e - d i m e n s i o n a l flow i n w h i c h a u t o m a t i c finite difference e q u a t i o n g e n e r a t o r s for s p a t i a l d e r i v a t i v e s w e r e i n c o r p o r a t e d i n t o a n explicit t i m e a d v a n c e m e n t a l g o r i t h m . L i n d e m u t h a n d K i l l e e n (1973) h a v e p u b l i s h e d a n a l g o r i t h m f o r t w o - d i m e n s i o n a l , a x i s y m m e t r i c flow in w h i c h all t e r m s a p p e a r i n g in t h e e q u a t i o n , i n c l u d i n g t r a n s p o r t t e r m s , a r e differenced implicitly. T h e s o l u t i o n o f t h e implicitly f o r m u l a t e d e q u a t i o n s is p e r f o r m e d w i t h a n a l t e r n a t i n g d i r e c t i o n i m p l i c i t m e t h o d ( P e a c e m a n a n d R a c h f o r d , 1955). A s in t h e e a r l i e r c o d e o f F r e e m a n a n d L a n e , diifusion t e r m s a r e a d d e d t o e a c h o f t h e e q u a t i o n s t o s m o o t h t h e s o l u t i o n . H o w e v e r , o n l y t h e m a s s diffusion coefficient is p u r e l y n u m e r i c a l (I. L i n d e m u t h , p r i v a t e c o m m u n i c a t i o n , 1975). D u c h s (1968) h a s p e r f o r m e d t w o - d i m e n s i o n a l s t u d i e s o f t h e r o t a t i o n i n d u c e d b y t r a n s v e r s e m u l t i p o l e fields i n a p l a s m a c o n f i n e d b y a t h e t a p i n c h field. T h e e q u a t i o n s o f m o t i o n for a two-fluid p l a s m a in a ζ = c o n s t , p l a n e of a n axisymmetric t h e t a pinch are solved, including t h e Hall term, b u t n e g l e c t i n g a n i s o t r o p i c p r e s s u r e effects a n d e l e c t r o n i n e r t i a t e r m s . L u i a n d C h u (1974) h a v e a p p l i e d a t w o - d i m e n s i o n a l n u m e r i c a l m e t h o d t o a t i m e - d e p e n d e n t , b o u n d a r y v a l u e p r o b l e m in a n a x i s y m m e t r i c t o r u s . T h e m o t i o n is c o m p u t e d in a p o l o i d a l p l a n e w i t h e q u a t i o n s w h i c h i n c l u d e t h e r m a l a n d resistive diffusion, b u t n o explicit viscosity. T h e v a l u e o f t h e m a g n e t i c flux o n t h e b o u n d a r y is p r e s c r i b e d a s a f u n c t i o n o f t i m e . A s t h e p l a s m a is p i n c h e d i n w a r d b y a n i n c r e a s i n g , a p p l i e d m a g n e t i c field, a l o w d e n s i t y p l a s m a r e g i o n is left b e h i n d . T o avoid p r o b l e m s with large Alfven speeds a n d the c o n s e q u e n t restrictions on the m a x i m u m time-step consistent with accurate c o m p u t a t i o n , the density

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS in t h e " v a c u u m " r e g i o n is n o t p e r m i t t e d t o fall b e l o w 1 5 % o f t h e i n i t i a l , u n i f o r m d e n s i t y . T h e diffusion e q u a t i o n s a r e a d v a n c e d b y m e a n s o f a t w o - s t e p s c h e m e w h e r e v a l u e s o f t h e diffusion coefficients a r e a d v a n c e d in t i m e o n t h e first s t e p , a n d i m p l i c i t e q u a t i o n s a r e u s e d t o a d v a n c e t h e d e p e n d e n t v a r i a b l e s o n the second step. T h e e q u a t i o n s of m o t i o n are also written implicitly. B o t h t h e diffusion e q u a t i o n s a n d t h e e q u a t i o n s o f m o t i o n a r e s o l v e d b y m e a n s o f a n alternating direction implicit algorithm. A two-dimensional, M H D , multifluid code which includes t r a n s p o r t d u e t o m i c r o t u r b u l e n c e is r e p o r t e d b y W a g n e r a n d M a n h e i m e r (1973). A n i s o t r o p i c h e a t c o n d u c t i o n a n d t h e r m o e l e c t r i c m a g n e t i c field g e n e r a t i o n a r e a l s o i n c l u d e d in the model. T h e code has been applied t o the simulation of the return c u r r e n t f r o m a relativistic e l e c t r o n b e a m flowing t h r o u g h a p l a s m a . Three-dimensional, nonlinear, ideal M H D calculations have been reported b y W o o t e n et al. (1974). C a l c u l a t i o n s o f t h e n o n l i n e a r e v o l u t i o n o f fixed b o u n d a r y k i n k m o d e s w e r e p e r f o r m e d w i t h a n explicit, l e a p f r o g a d v a n c e m e n t a l g o r i t h m . T h e m e t h o d is a n e x t e n s i o n o f a c o m p u t a t i o n a l m e t h o d f o r t h e s o l u t i o n o f t h e l i n e a r i z e d M H D e q u a t i o n s r e p o r t e d b y B a t e m a n et al. (1974). 2 . Lagrangian

Computations

One-dimensional Lagrangian c o m p u t a t i o n s of pinch discharges have been p e r f o r m e d for m a n y y e a r s ( H a i n et al., 1960), b u t t h e r e is o n l y o n e recently reported L a g r a n g i a n calculation in t w o dimensions. T h e calculation o f t h e n o n l i n e a r e v o l u t i o n o f t h e k i n k m o d e in a T o k a m a k b y W h i t e et al. (1974) f o l l o w s t h e m o t i o n o f t h e b o u n d a r y o f a n i n c o m p r e s s i b l e p l a s m a d u e t o the p e r t u r b a t i o n of a n unstable equilibrium. A cylindrical a p p r o x i m a t i o n t o t h e t o r o i d a l g e o m e t r y o f t h e p l a s m a is u s e d , a n d helical s y m m e t r y is a s s u m e d . I n t h i s g e o m e t r y a n d b e c a u s e o f t h e i n c o m p r e s s i b i l i t y o f t h e p l a s m a , t h e field a n d velocity m a y be c o m p u t e d f r o m potential functions. A closed sequence o f s t r a i g h t line s e g m e n t s lying o n t h e p l a s m a b o u n d a r y f o r m t h e L a g r a n g i a n computation mesh. 3. Computations

with a Generalized

Mesh

H e r t w e c k a n d S c h n e i d e r (1970) r e p o r t e d a m e t h o d f o r c a l c u l a t i n g e n d losses in a t h e t a p i n c h u s i n g flux c o o r d i n a t e s . T h e m e s h is s t a t i o n a r y in t h e axial direction, b u t m o v e s radially so t h a t the relative velocity between the flux s u r f a c e s a n d t h e m e s h is z e r o w h e n t h e c o n d u c t i v i t y is infinite. T h e differential e q u a t i o n s a r e t r a n s f o r m e d i n t o t h e n o n o r t h o g o n a l , t i m e - d e p e n d e n t c o o r d i n a t e s a n d t h e n differenced. T h e differential e q u a t i o n s a p p e a r q u i t e c o m p l e x in t h e s e c o o r d i n a t e s , b u t a r e s i m p l y E q s . (11)—(14) w i t h r e f e r e n c e t o t h e m o v i n g f r a m e . A s H e r t w e c k a n d S c h n e i d e r (1970) p o i n t o u t , t h e a c c u r a t e

9

10

JEREMIAH U. BRACKBILL

c a l c u l a t i o n of field a l i g n e d t h e r m a l c o n d u c t i o n is o n l y p o s s i b l e u s i n g flux c o o r d i n a t e s . T h e r e s u l t s o f a c a l c u l a t i o n o f e n d losses i n c l u d i n g t h e r m a l c o n d u c t i o n w e r e l a t e r r e p o r t e d b y S c h n e i d e r (1972). A m e t h o d d e v e l o p e d b y A n d e r s o n (1975) a l s o u s e s flux c o o r d i n a t e s . H o w e v e r , t h e c o m p u t a t i o n m e s h is o r t h o g o n a l . Z o n e lines p e r p e n d i c u l a r t o t h e flux surfaces a r e r e c o n s t r u c t e d e a c h c o m p u t a t i o n t i m e - s t e p m a k i n g p o s ­ sible t h e u s e o f a t i m e - s t e p s p l i t t i n g a l g o r i t h m t o a d v a n c e t h e s o l u t i o n i n t i m e ( R i c h t m y e r a n d M o r t o n , 1967). W i t h t i m e - s t e p s p l i t t i n g a m u l t i d i m e n s i o n a l c a l c u l a t i o n is p e r f o r m e d b y a s e q u e n c e o f o n e - d i m e n s i o n a l o p e r a t i o n s s o t h a t it is c o n v e n i e n t t o a p p l y t h e flux c o r r e c t e d t r a n s p o r t a l g o r i t h m o f B o r i s a n d B o o k (1973) t o r e d u c e c o m p u t a t i o n a l diffusion. A n i m p l i c i t l y f o r m u l a t e d m e t h o d i n w h i c h t h e c o m p u t a t i o n is split i n t o Lagrangian a n d convective transport phases has been applied to M H D c o m p u t a t i o n s b y B r a c k b i l l a n d P r a c h t (1973). T w o - d i m e n s i o n a l c a l c u l a t i o n s o f a z - p i n c h in r - 0 c o o r d i n a t e s a n d a x i s y m m e t r i c c a l c u l a t i o n s o f a h o t d i a m a g n e t i c pellet e x p a n d i n g i n t o a m a g n e t i c field a r e p e r f o r m e d u s i n g a n a l m o s t - L a g r a n g i a n mesh. In a n a l m o s t - L a g r a n g i a n mesh, the relative m o t i o n b e t w e e n t h e m e s h a n d t h e p l a s m a is k e p t t o a m i n i m u m t o r e d u c e c o m p u t a t i o n a l diffusion. A n a x i s y m m e t r i c c a l c u l a t i o n o f a s h a r p b o u n d a r y t h e t a p i n c h ( B r a c k b i l l , 1973) w a s r e p o r t e d i n w h i c h a v a c u u m field s o l u t i o n w a s c o u p l e d t o p l a s m a flow a c r o s s a L a g r a n g i a n i n t e r f a c e . A s i m i l a r m e t h o d for r e s o l v i n g t h e p l a s m a - v a c u u m i n t e r f a c e w a s r e p o r t e d b y L u i (1973). L u i u s e d a n E u l e r i a n m e s h w i t h a r b i t r a r i l y s h a p e d cells t o resolve curved b o u n d a r i e s , a n d L a g r a n g i a n m e s h points to resolve a m o v i n g v a c u u m - p l a s m a i n t e r f a c e . T h e m e t h o d w a s a p p l i e d t o t h e s t u d y o f t h e reflection a n d focusing of shocks by curved walls. A m e t h o d for c o m p u t i n g t h r e e - d i m e n s i o n a l M H D flow w i t h a g e n e r a l i z e d m e s h is p r e s e n t e d i n S e c t i o n I V , a n d t h e r e s u l t s o f its a p p l i c a t i o n t o t h e rotating theta pinch a n d the internal kink m o d e are displayed.

III. T h e Computation of Convective Transport

A . PROPERTIES OF APPROXIMATIONS TO THE CONVECTIVE DERIVATIVE T h e f u n d a m e n t a l difficulty w i t h a c a l c u l a t i o n u s i n g a n E u l e r i a n m e s h is t h e a c c u r a t e c o m p u t a t i o n o f c o n v e c t i v e t r a n s p o r t . W h e t h e r t h e e q u a t i o n s for a n E u l e r i a n m e s h a r e w r i t t e n in c o n s e r v a t i o n f o r m , E q s . ( l ) - ( 4 ) , o r in n o n c o n s e r v a t i o n f o r m , E q s . ( 7 ) - ( 1 0 ) o r E q s . (15)—(18), it is t h e a p p r o x i m a t i o n o f t h e c o n v e c t i v e d e r i v a t i v e , u · V, w h i c h d e t e r m i n e s t h e o v e r a l l a c c u r a c y a n d stability o f t h e s o l u t i o n . O n e m i g h t r e a s o n a b l y a s k : S i m p l e difference s c h e m e s approximating the convective derivative can be m a d e t o have the same formal

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

11

a c c u r a c y a s a n y o t h e r difference e q u a t i o n , s o w h y d o e s t h e c o m p u t a t i o n o f t h e c o n v e c t i v e d e r i v a t i v e i n t r o d u c e e r r o r s ? T h e a n s w e r is t h a t if t h e t e r m s in t h e difference e q u a t i o n s c o r r e s p o n d i n g t o t h e c o n v e c t i v e d e r i v a t i v e a r e w r i t t e n i n a s i m p l e w a y , t h e e q u a t i o n s a r e o f t e n u n s t a b l e . O t h e r , less a c c u r a t e difference e q u a t i o n s w h i c h a r e s t a b l e m u s t b e u s e d i n s t e a d ( R i c h t m y e r , 1963).

1. Nonlinear

Stability

L i n e a r s t a b i l i t y a n a l y s i s is n o t a l w a y s i n f o r m a t i v e a b o u t t h e s t a b i l i t y o f difference a p p r o x i m a t i o n s t o t h e c o n v e c t i v e d e r i v a t i v e . I t is n e c e s s a r y t h a t t h e difference a p p r o x i m a t i o n s b e l i n e a r l y s t a b l e , b u t l i n e a r l y s t a b l e difference equations are

often

nonlinearly

unstable.

For

example,

Richtmyer

and

M o r t o n (1967) d i s c u s s t h e a p p l i c a t i o n o f t h e l e a p f r o g difference e q u a t i o n t o a nonlinear convective transport equation, 2

du/dt + (d/dx)(iu )

= 0,

0 < χ < 1,

t^

0.

(19)

T h e " l e a p f r o g " a p p r o x i m a t i o n t o t h e e q u a t i o n is w r i t t e n »+it

ty

= *-iUj

2

2

+ Ηδί/δχ)Κ\+1)

- (Χ-χ) ],

(20)

a n d is l i n e a r l y s t a b l e ( R i c h t m y e r , 1963). N e v e r t h e l e s s , it is e x p l o s i v e l y u n s t a b l e w h e n n u m e r i c a l l y s o l v e d . T h a t is, t h e s o l u t i o n g r o w s e x p o n e n t i a l l y i n s o m e 1

c a s e s w i t h a g r o w t h r a t e w h i c h g o e s a s (St)' .

A n analysis of the p r o b l e m

i n d i c a t e s t h e l e a p f r o g a p p r o x i m a t i o n is s t a b l e w h e n St is sufficiently s m a l l s o t h a t " t h e t r u n c a t i o n e r r o r i n t h e difference e q u a t i o n s g e n e r a t e s sufficiently s m a l l p e r t u r b a t i o n s . " T h e difference b e t w e e n l i n e a r a n d n o n l i n e a r p r o b l e m s w h i c h c a u s e s t h e difference i n s t a b i l i t y is t h a t " t h e b e h a v i o r d e p e n d s m a r k e d l y o n the relative m a g n i t u d e of the p e r t u r b a t i o n s " ( R i c h t m y e r a n d

Morton,

1967). S i m i l a r c o n c l u s i o n s a r e d r a w n b y F o r n b e r g (1973), i n a n a n a l y s i s o f t h e l e a p f r o g a n d C r a n k - N i c o l s o n a p p r o x i m a t i o n s t o E q . (19). 2. The Stability Equations

of Numerical

with Variable

Approximations

Coefficients:

Kreiss's

to Linear

Partial

Differential

Theorem

T h e analysis of the stability of n u m e r i c a l a p p r o x i m a t i o n s t o linear partial differential e q u a t i o n s w i t h v a r i a b l e coefficients, s u c h a s t h e e q u a t i o n du/dt

= a(x)

du/dx

= 0,

(21)

is t h e n e x t s t e p b e y o n d l i n e a r s t a b i l i t y a n a l y s i s i n d e t e r m i n i n g t h e s t a b i l i t y o f n u m e r i c a l a p p r o x i m a t i o n s t o n o n l i n e a r e q u a t i o n s . A s u m m a r y is g i v e n o f a n

JEREMIAH U. BRACKBILL

12

i m p o r t a n t t h e o r e m b y K r e i s s (1964) w h o s e r e l e v a n c e t o t h e p r e s e n t d i s c u s s i o n is s h o w n in t h e f o l l o w i n g s e c t i o n . K r e i s s ' s t h e o r e m c a n b e s u m m a r i z e d , for o u r p u r p o s e s , in t h e f o l l o w i n g w a y . C o n s i d e r a l i n e a r l y s t a b l e difference e q u a t i o n a p p r o x i m a t i o n t o E q . (21). If t h e t r u e s o l u t i o n t o E q . (21) w e r e s u b s t i t u t e d i n t o t h e a p p r o x i m a t i o n , t h e r i g h t a n d left sides w o u l d b e u n e q u a l m n b y t e r m s o f o r d e r dx a n d 3t , w h e r e δχ is t h e m e s h s p a c i n g a n d St is t h e t i m e s t e p . T h e s e t e r m s a r e t h e t r u n c a t i o n e r r o r s in t h e a p p r o x i m a t i o n . T h e n if a diffusionlike t e r m exists o r is a d d e d t o t h e difference e q u a t i o n w i t h p o s i t i v e diffusion coefficients o f a t least o n e l o w e r o r d e r t h a n t h e t r u n c a t i o n e r r o r s , m 1_ 1 t h a t is, o f o r d e r < 5 x a n d öt"' o r less, t h e e q u a t i o n will b e s t a b l e . A n e q u a t i o n w i t h s u c h a diffusion t e r m is c a l l e d p o s i t i v e l y s t a b l e b y R i c h t m y e r (1963). E v i d e n t l y , t h e r o l e o f t h e diffusion is t o a d d sufficient s m o o t h i n g t o t h e s o l u t i o n s o t h a t p e r t u r b a t i o n s d u e t o t r u n c a t i o n e r r o r s in t h e difference e q u a t i o n a r e k e p t a t a sufficiently l o w level. 3. A Heuristic

Stability

Theory

K r e i s s ' s t h e o r e m s h o w s t h a t t h e p r e s e n c e o f e n o u g h diffusion t o d a m p t r u n c a t i o n - e r r o r - p r o d u c e d p e r t u r b a t i o n s is a n e c e s s a r y a n d sufficient c o n ­ d i t i o n for t h e stability o f difference a p p r o x i m a t i o n s t o l i n e a r differential e q u a t i o n s w i t h v a r i a b l e coefficients. T h a t t h i s c o n d i t i o n is, a t least, a n e c e s s a r y c o n d i t i o n for t h e stability o f difference a p p r o x i m a t i o n s t o n o n l i n e a r e q u a t i o n s is s u g g e s t e d b y t h e w o r k o f H i r t (1968) in d e v e l o p i n g a h e u r i s t i c s t a b i l i t y t h e o r y . I n his s t u d y , finite difference a p p r o x i m a t i o n s a r e r e d u c e d t o differential e q u a t i o n s b y e x p a n d i n g e a c h o f t h e i r t e r m s in a T a y l o r series. T h e l o w e s t o r d e r t e r m s in t h e e x p a n s i o n s a r e t h e differential e q u a t i o n s b e i n g a p p r o x i m a t e d . H i g h e r o r d e r t e r m s in δχ a n d δί a r e t r u n c a t i o n e r r o r s . I t is H i r t ' s o b s e r v a t i o n t h a t s o m e t r u n c a t i o n errors are associated with n o n l i n e a r instabilities of the difference e q u a t i o n s . Specifically, h e e x h i b i t s cases w h e r e t h e coefficient o f t h e s e c o n d s p a t i a l d e r i v a t i v e o f t h e d e n s i t y in t h e differential e q u a t i o n d e r i v e d f r o m t h e difference e q u a t i o n a p p r o x i m a t i n g t h e c o n t i n u i t y e q u a t i o n a p p e a r s t o d e t e r m i n e t h e stability o f a s y s t e m o f c o u p l e d e q u a t i o n s d e s c r i b i n g fluid flow. T h a t is, w h e n t h e coefficient is n e g a t i v e , t h e r e is antidiffusion in t h e continuity e q u a t i o n which destabilizes the entire system of equations. T h e c o n v e r s e o f K r e i s s ' s t h e o r e m for l i n e a r e q u a t i o n s s e e m s t o b e t r u e for n o n ­ l i n e a r e q u a t i o n s a l s o . W i t h o u t t h e a d d i t i o n of p o s i t i v e diffusion o f t h e s a m e o r d e r a s t h e t r u n c a t i o n e r r o r s , t h e difference e q u a t i o n s a r e u n s t a b l e . F o r stability, t h e e q u a t i o n s m u s t b e d i s s i p a t i v e o r p o s i t i v e l y s t a b l e in t h e sense d e s c r i b e d b y K r e i s s . E v i d e n t l y , H i r t ' s h e u r i s t i c stability t h e o r y is p o i n t i n g t o t h e e x i s t e n c e of a s e c o n d n e c e s s a r y c o n d i t i o n for t h e stability o f t h e difference e q u a t i o n s . T h a t is, n o t o n l y m u s t t h e difference e q u a t i o n s b e linearly s t a b l e , b u t t h e y m u s t a l s o satisfy K r e i s s ' s t h e o r e m .

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

13

A t r u n c a t i o n e r r o r a n a l y s i s o f v a r i o u s difference a p p r o x i m a t i o n s t o E q . (15), t h e c o n v e c t i v e d e r i v a t i v e o f t h e d e n s i t y , is g i v e n b y H i r t (1968). A m o n g t h e e q u a t i o n s for w h i c h his r e s u l t s a r e a p p l i c a b l e a r e t h r e e w h i c h a r e l i n e a r l y s t a b l e ( R i c h t m y e r , 1963); u p s t r e a m o r d o n o r cell differencing, l e a p f r o g differencing, a n d i m p l i c i t s p a t i a l l y - c e n t e r e d differencing, T h e s e difference a p p r o x i m a t i o n s a r e listed in T a b l e I, a l o n g w i t h t h e o r d e r o f a c c u r a c y . I n T a b l e I I , t h e c o r r e s p o n d i n g coefficients o f t h e s e c o n d d e r i v a t i v e o f t h e d e n s i t y a r e listed, i n c l u d i n g n o n l i n e a r t e r m s . I n e a c h c a s e , t h e n o n l i n e a r t e r m s a r e p r o p o r t i o n a l t o velocity gradients, a n d w o u l d a p p e a r n o m a t t e r which f o r m o f t h e c o n t i n u i t y e q u a t i o n w e r e a p p r o x i m a t e d . (It is a l s o t r u e t h a t i n e a c h c a s e , t h e n o n l i n e a r t e r m s a r e p r o p o r t i o n a l t o g r a d i e n t s o f u", t h e r e l a t i v e v e l o c i t y b e t w e e n t h e p l a s m a a n d t h e m e s h . T h u s , it is p o s s i b l e t o t r a n s f o r m t o a c o o r d i n a t e f r a m e w h e r e n o n l i n e a r t r u n c a t i o n e r r o r s a r e negligible, a fact which can be exploited with a generalized mesh method.) A s w e c a n see f r o m t h e e n t r i e s in T a b l e I I , all o f t h e m a s s diffusion c o ­ efficients b e c o m e n e g a t i v e w i t h sufficiently l a r g e velocity g r a d i e n t s . T h i s r e s u l t is in a g r e e m e n t w i t h R i c h t m y e r ' s (1963) o b s e r v a t i o n s . I t a l s o h e l p s o n e t o u n d e r s t a n d h o w t h e a d d i t i o n o f viscosity t o t h e m o m e n t u m e q u a t i o n e n h a n c e s t h e s t a b i l i t y o f difference e q u a t i o n s , for t h e effect o f viscosity is t o d e c r e a s e velocity g r a d i e n t s . S i m i l a r e r r o r s a r e f o u n d in t h e difference a p p r o x i m a t i o n s t o E q s . (16)—(18), t h e c o n v e c t i v e d e r i v a t i v e o f t h e m o m e n t u m , m a g n e t i c flux, a n d e n e r g y . E v i d e n t l y , t h e difficulty w i t h m o s t a p p r o x i m a t i o n s t o c o n v e c t i v e t r a n s p o r t is t h a t t h e y d o n o t satisfy t h e s e c o n d n e c e s s a r y c o n d i t i o n for s t a b i l i t y g i v e n by Kreiss's t h e o r e m . A s H i r t ' s analysis shows, the violation of the a s s u m p t i o n s o f K r e i s s ' s t h e o r e m is often a s s o c i a t e d w i t h n o n l i n e a r n u m e r i c a l i n s t a b i l i t i e s .

B . NONLINEARLY STABLE APPROXIMATIONS TO CONVECTIVE TRANSPORT 1. Truncation

Error

Corrections

A m s d e n a n d H i r t (1973) a n d R i v a r d et al (1973) h a v e a p p l i e d H i r t ' s h e u r i s t i c stability t h e o r y t o t h e a p p l i c a t i o n o f t r u n c a t i o n e r r o r c o r r e c t i o n s o f t h e difference e q u a t i o n s for fluid flow. F o r e x a m p l e , t o stabilize t h e l e a p ­ frog differencing o f t h e c o n t i n u i t y e q u a t i o n ( e n t r y b , T a b l e I), a m a s s diffusion t e r m is a d d e d w i t h diffusion coefficient, κ, given b y κ = φ χ

2

du"/dx

(22)

w h e r e α is t y p i c a l l y e q u a l t o 2. W i t h o u t d e c r e a s i n g t h e o r d e r o f a c c u r a c y o f t h e o r i g i n a l difference a p p r o x i m a t i o n , diffusion is a d d e d t o t h e difference e q u a t i o n t o stabilize it. T h e a d d e d diffusion s h o u l d b e sufficient b e c a u s e it

c. Implicit, spatially-centered differencing

b. Leapfrog differencing

a. Upstream differencing

Name

TABLE I

n n Pj = Pj —

+ 1

n

Vj -

1 ("+ (pu)j+ι

+ 12 /

> 0

1 — n + (pu)j-1)

-n(pu)j-1)

uJ

^«^+ι/2>η^+ι/2-<ηΡ^ί/2>Λ^-ι,2)

n (n(pu)j+i

= < lnpj,

-

pJ =

+ 1

pj = ~n 1pJ

Difference equation

DIFFERENCE APPROXIMATIONS TO THE CONTINUITY EQUATION

0(<5* 3,
0(<5JC3, <5r3)

0(δχ\δί2)

Truncation error

14 JEREMIAH U. BRACKBILL

FI«-i

—u 2

2

δχ2 du"

2 dx

~ ~T~dx~

δχ du"

2

~~2~dx

~

δχ2 du"

8

öxAd*u"

"Ιϊ'δχ1

öxAd3u"

~~V2~dx*

öx+d3u"

dx ~lÄ'dx1

2

δχ3 d2u"

Mass diffusion coefficient"

' u" is the relative velocity between the mesh and the plasma (u" = u — u').

c. Implicit, spatially centered differencing

b. Leapfrog differencing

a. Upstream differencing

Name

MASS DIFFUSION COEFFICIENTS DUE TO TRUNCATION ERRORS

TABLE II

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS 15

16

JEREMIAH U. BRACKBILL

c a n c e l s t h e l o w e s t o r d e r t r u n c a t i o n e r r o r , a n d a d d s s o m e p o s i t i v e diffusion t o c a n c e l h i g h e r o r d e r e r r o r s a s well. S o m e c o n f i d e n c e i n t h i s a p p r o a c h is d e r i v e d f r o m its r e d u c i n g t o a p r e s c r i p t i o n w h i c h satisfies K r e i s s ' s t h e o r e m for l i n e a r e q u a t i o n s w i t h v a r i a b l e coefficients. R i v a r d et al. (1973) a p p l y s i m i l a r t r u n c a t i o n e r r o r c o r r e c t i o n s t o all o f t h e fluid e q u a t i o n s . A s o m e w h a t s i m i l a r a p p r o a c h is t a k e n b y L a p i d u s (1967), in t h a t diffusion is a d d e d t o every e q u a t i o n . H o w e v e r , i n h i s a p p r o a c h a n d i n its a p p l i c a t i o n s t o M H D c a l c u l a t i o n s ( F r e e m a n a n d L a n e , 1 9 6 8 ; L i n d e m u t h a n d K i l l e e n , 1973) t h e a d d e d diffusion is n o t c o m p u t e d f r o m a t r u n c a t i o n error analysis. 2. Flux-Corrected

Transport

A m e t h o d a p p r o x i m a t i n g c o n v e c t i v e t r a n s p o r t , w h i c h is d i s c u s s e d i n d e t a i l i n a n o t h e r article i n t h i s v o l u m e , h a s b e e n d e v e l o p e d b y B o r i s a n d B o o k (1973). T h e i r m e t h o d r e q u i r e s t h a t t h e s o l u t i o n o f t h e difference e q u a t i o n for t h e c o n v e c t i v e t r a n s p o r t o f m a s s a n d e n e r g y b e p o s i t i v e e v e r y w h e r e . T h e m e t h o d is c o m p o s e d o f t w o s t e p s . I n t h e first s t e p , a l a r g e v e l o c i t y - i n d e p e n d e n t diffusion is a d d e d t o t h e difference a p p r o x i m a t i o n o f t h e c o n v e c t i v e d e r i v a t i v e . T h e a d d i t i o n o f t h e diffusion a s s u r e s t h a t t h e s o l u t i o n o b t a i n e d i n t h e first s t e p is p o s i t i v e e v e r y w h e r e . I n t h e s e c o n d s t e p t h e diffusion a d d e d i n t h e first s t e p is s u b t r a c t e d e v e r y w h e r e , s u b j e c t t o t h e c o n s t r a i n t t h a t t h e s o l u t i o n obtained be positive everywhere. I n t h e first s t e p , t h e c o n v e c t i v e d e r i v a t i v e i n t h e c o n t i n u i t y e q u a t i o n is a p p r o x i m a t e d in o n e dimension by t h e equation f

n

Pj = "Pj " W öt/Sx)( pj+1

- -py-i)

+ [* + W'ttlöxfWp;^

- 2

n Pj

+ »p,-i),

(23)

w h e r e ρ d e n o t e s t h e s o l u t i o n o b t a i n e d o n t h e first s t e p . T h e difference e q u a t i o n is t h e L a x - W e n d r o f f e q u a t i o n ( R i c h t m y e r , 1963) p l u s a n a d d i t i o n a l velocityi n d e p e n d e n t diffusion. T o first o r d e r i n δχ a n d St, t h e difference e q u a t i o n a p p r o x i m a t e s t h e differential e q u a t i o n 2

dp/dt = - u" dp/dx + i(Sx /St)

2

2

(d p/dx ).

(24) 2

T h e diffusion coefficient f o r t h e v e l o c i t y - i n d e p e n d e n t diffusion is t h u s %öx /St. T h i s diffusion is s u b t r a c t e d in t h e s e c o n d s t e p o f t h e a l g o r i t h m , w h e r e n e g a t i v e diffusion is c o m p u t e d . T h e n e g a t i v e o r antidiffusion is c o m p u t e d subject t o the constraint that t h e second step create n o n e w m a x i m a a n d m i n i m a in t h e s o l u t i o n , a n d t h a t it n o t a c c e n t u a t e a l r e a d y existing e x t r e m a . I n t h e a n t i -

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

17

diffusion s t e p , m a s s a n d e n e r g y fluxes b e t w e e n p a i r s o f cells a r e c o m p u t e d , a s s u r i n g t h a t t h e t o t a l m a s s a n d t o t a l e n e r g y a r e c o n s e r v e d . T h e antidiffusive mass

fluxes,/,

are computed from the equation

/y+i/2 = S G N A

i + j 2/

max{0, m i n ( A , _ 1

A,.+ 3 / S 2 G N A , . +

/2

S G N Δ 7 · + ι / 2, " i " | Δ 7 · + 1 / 2| ,

1 )/ }2 ,

(25)

where A J + i/2 = Pj+i The symbol

~ Pj-

signifies t h a t h i g h e r o r d e r c o r r e c t i o n s a r e m a d e s o t h a t t h e

a n t i d i f f u s i o n i n t h e s e c o n d s t e p c a n c e l s t h e diffusion a d d e d i n t h e first s t e p t o f o u r t h o r d e r i n δχ a n d δί w h e r e v e r p o s s i b l e . T h e n e w d e n s i t y is t h e n c o m ­ puted from the equation n+l

Pj

= Pj-fJ+1/2+fj-1/2.

(26)

Similar e q u a t i o n s are written for the energy. T h e flux e q u a t i o n , E q . (25), a l l o w s n o a d j u s t m e n t s i n t h e s e c o n d s t e p t o l o c a l m a x i m a a n d m i n i m a . W h e r e t h e v a r i a t i o n o f ßj is m o n o t o n i c , t h e a n t i diffusion will c o n t i n u e u n t i l e i t h e r Δ 7 · _ 1 equals either p n+1

s o l u t i o n , p,is

j

+l

o r pj^1.

2/

or Δ 7 +

3 2/

is z e r o , t h a t is, u n t i l ßj

O n e c a n see f r o m t h i s h o w t h e p o s i t i v i t y o f t h e

p r e s e r v e d , for t h e m i n i m u m v a l u e o f pj o v e r t h e e n t i r e d o m a i n

cannot be decreased t o w a r d zero, a n d n o other value can be m a d e smaller than the minimum. Flux-corrected approximations to the convective derivative have been u s e d i n c o m p u t a t i o n s o f t w o - a n d t h r e e - d i m e n s i o n a l flow e m p l o y i n g a t i m e step splitting algorithm t o a d v a n c e the solution in time. I n such c o m p u t a t i o n s , o n e - d i m e n s i o n a l difference e q u a t i o n s i n e a c h o r t h o g o n a l c o o r d i n a t e d i r e c t i o n a r e s o l v e d i n a cyclic f a s h i o n . T h u s , a o n e - d i m e n s i o n a l f o r m u l a t i o n o f

flux-

corrected transport can be used.

I V . A Generalized M e s h M e t h o d for M H D I n t h i s s e c t i o n , a m e t h o d is p r e s e n t e d f o r c o m p u t i n g t h r e e - d i m e n s i o n a l M H D flow o n a g e n e r a l i z e d m e s h . T h e e q u a t i o n s o f m o t i o n a r e split i n t o t w o p h a s e s a s p r o p o s e d b y H i r t et al. (1974) a n d d i s c u s s e d i n S e c t i o n I I . I n t h e first p h a s e , t h e L a g r a n g i a n e q u a t i o n s o f m o t i o n , E q s . ( 7 ) - ( 1 0 ) , a r e s o l v e d e x a c t l y a s o n e w o u l d solve t h e m for a p u r e L a g r a n g i a n c o m p u t a t i o n . I n t h e s e c o n d p h a s e , t h e t r a n s p o r t o f m a s s , m o m e n t u m , flux, a n d e n e r g y d u e t o r e l a t i v e m o t i o n b e t w e e n t h e p l a s m a a n d t h e m e s h is c o m p u t e d f r o m E q s .

18

JEREMIAH U. BRACKBILL

(15)—(18). ( T h e s e c o n d p h a s e is c h a r a c t e r i z e d a s a r e z o n e , w h e r e t h e r e z o n e is w i t h r e s p e c t t o t h e L a g r a n g i a n f r a m e , r a t h e r t h a n t h e fixed, E u l e r i a n f r a m e . ) T h e p r e s c r i p t i o n f o r t h e g r i d velocity, u', is c o m p l e t e l y a r b i t r a r y . By c h o o s i n g u' p r o p e r l y , t h e g e n e r a l i z e d m e s h m a y b e m a d e E u l e r i a n , L a g r a n g i a n , flux c o o r d i n a t e ( H e r t w e c k a n d S c h n e i d e r , 1970), o r a l m o s t - L a g r a n g i a n ( B r a c k b i l l a n d P r a c h t , 1973). T h e a l m o s t - L a g r a n g i a n p r e s c r i p t i o n for u' will be outlined in a brief p a r a g r a p h in this section. T h e b u l k o f t h i s s e c t i o n will b e d e v o t e d t o d e v e l o p i n g g e n e r a l i z e d difference e q u a t i o n s for the L a g r a n g i a n p h a s e of the calculation. These e q u a t i o n s describe the d y n a m i c a l evolution of the system, a n d m u s t be solved n o m a t t e r w h i c h r e z o n i n g p r e s c r i p t i o n is e m p l o y e d in t h e r e z o n e p h a s e . I t will b e s h o w n h o w the e q u a t i o n s are derived, h o w they are solved, h o w b o u n d a r y conditions are applied, a n d w h a t their properties are.

A. DIFFERENCE EQUATIONS FOR THE LAGRANGIAN PHASE OF A GENERALIZED MESH CALCULATION Difference e q u a t i o n s for n o n r e c t i l i n e a r m e s h e s m a y b e d e r i v e d b y t h e finite-element m e t h o d . A s t h e m e t h o d is o u t l i n e d b y Z i e n k i e w i c z (1971), t h e flow field is s u b d i v i d e d i n t o e l e m e n t s , w h e r e e a c h e l e m e n t is a s s o c i a t e d w i t h a n u m b e r o f m e s h p o i n t s a t w h i c h t h e v a l u e s o f t h e flow v a r i a b l e s a r e s t o r e d . A t e v e r y p o i n t in e a c h e l e m e n t , t h e flow v a r i a b l e s a r e d e t e r m i n e d b y i n t e r ­ p o l a t i o n so t h a t they are c o n t i n u o u s between elements. T o derive e q u a t i o n s o f m o t i o n , a f u n c t i o n a l is defined a n d m i n i m i z e d o v e r t h e flow field w i t h t h e v a l u e o f v a r i a b l e s a t m e s h p o i n t s a s p a r a m e t e r s . I n M H D flow, t h e a p p r o p r i a t e f u n c t i o n a l is t h e L a g r a n g i a n . In contrast with this a p p r o a c h the derivation presented here proceeds d i r e c t l y f r o m t h e e q u a t i o n s o f m o t i o n . V a r i a b l e s a r e defined o v e r t h e m e s h a s i n t h e finite e l e m e n t m e t h o d , b u t t h e k i n e m a t i c a l e q u a t i o n s a r e a p p r o x i ­ m a t e d b y a c o o r d i n a t e t r a n s f o r m a t i o n like t h a t u s e d b y S c h u l z (1964). A d e r i v a t i o n for t h e d y n a m i c a l e q u a t i o n after G o a d (1960) u s i n g d ' A l e m b e r t ' s principle and the principle of virtual w o r k completes the formulation of L a g r a n g i a n difference e q u a t i o n s . 1. Kinematical

Equations

T h e k i n e m a t i c a l e v o l u t i o n o f t h e t h e r m o d y n a m i c v a r i a b l e s is c o n t a i n e d in t h e c o n t i n u i t y , i n d u c t i o n , a n d i n t e r n a l e n e r g y e q u a t i o n s , E q s . (7), (8), a n d (10). T h e s e e q u a t i o n s d e s c r i b e t h e e v o l u t i o n o f t h e t h e r m o d y n a m i c v a r i a b l e s ρ, B, a n d i as the p l a s m a dilates a n d shears. W h e n the m o m e n t u m e q u a t i o n , E q . (9), is s o l v e d s i m u l t a n e o u s l y a n d c o n s i s t e n t l y w i t h t h e k i n e m a t i c e q u a t i o n s , t h e full d y n a m i c s o f t h e s y s t e m d e s c r i b e d b y t h e s e e q u a t i o n s a r e o b t a i n e d .

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

19

I t is c o n v e n i e n t , w h e n d e r i v i n g t h e c o m p l e t e s y s t e m o f difference e q u a t i o n s , t o c o n s i d e r t h e k i n e m a t i c a l a n d d y n a m i c a l e q u a t i o n s s e p a r a t e l y . T h a t is, difference e q u a t i o n s a p p r o x i m a t i o n s t o E q s . (7), (8), a n d (10) a r e d e r i v e d s e p a r a t e l y f r o m t h e difference e q u a t i o n for E q . (9). T h e difference e q u a t i o n s a r e s o l v e d o n a c o m p u t a t i o n m e s h c o m p o s e d o f t h e six-surfaced c o n f o r m i n g cells d e s c r i b e d b y P r a c h t (1975). T h e s e cells a r e holomorphic to the unit cube; the coordinate transformation from a point in t h e u n i t c u b e (ξ, η, γ) t o a p o i n t in p h y s i c a l s p a c e χ is t r i l i n e a r , a n d is given b y 2

3

χ = [ x ^ ( l - η) + χ ξη 5

4

+ x ( l - ξ)η + x ( l - 0 ( 1 - ij)] (1 - v) 6

+ [ x £ ( l - η) + χ ξη

7

8

+ x ( l - ξ)η + x ( l - ξ)(I

- q ) ] ν,

(27)

w h e r e χ* is t h e p o s i t i o n o f t h e Ith v e r t e x o f t h e c o m p u t a t i o n cell w i t h i n w h i c h l χ lies. W h e n b o t h x* a n d x* = u a r e s t o r e d , t h e e l e m e n t s o f t h e r a t e o f s t r a i n t e n s o r , Saß, w r i t t e n Sxß

= dujdxß9

oc9ß = 1 , 2 , 3 ,

(28)

a r e d e t e r m i n e d in t h e i n t e r i o r o f a cell. A c o n s i s t e n t k i n e m a t i c a l d e s c r i p t i o n o f t h e flow f o l l o w s w h e n t h e t h e r m o d y n a m i c v a r i a b l e s a r e c o n s t a n t w i t h i n e a c h cell a n d d i s c o n t i n u o u s a t cell b o u n d a r i e s , a n d ΞΛβ is r e p l a c e d b y its a v e r a g e v a l u e o v e r t h e cell in E q s . (7), (8), a n d (10). T h e a v e r a g e v a l u e o f Saß o v e r a cell is c o m p u t e d f r o m t h e e q u a t i o n



W h e n (Saßy

=

±fydVSxß.

(29)

is s u b s t i t u t e d for Saß t h e k i n e m a t i c a l e q u a t i o n s a r e w r i t t e n rfp/A + p < S „ > = 0,

(30)

ρ di/dt + P < S „ > = 0,

(31)

and dBJdt

+ Ba - Bi < 5 t e> = 0,

(32)

w i t h t h e c o n v e n t i o n t h a t t h e r e p e a t e d s u b s c r i p t s i,j,k are to be s u m m e d o v e r . T h e t i m e d e r i v a t i v e s a n d t h e t i m e levels o f o t h e r t e r m s a p p e a r i n g in t h e s e e q u a t i o n s a r e n o t yet specified o r a p p r o x i m a t e d . T h e i n t e g r a l in E q . (29) is e v a l u a t e d b y r e w r i t i n g t h e i n t e g r a n d s o t h a t it

20

JEREMIAH U. BRACKBILL

d e p e n d s explicitly o n ξ9 η, a n d ν,

dV(cuJdxß)

= άξ άη dv d (ua, Xj, xk)/d (£, η, ν) J

(33)

w h e r e ßjk is a cyclic p e r m u t a t i o n of 1 , 2 , 3 . T h e i n t e g r a t i o n v o l u m e in ξ, η, ν s p a c e is t h e u n i t c u b e . I n t e g r a t i n g E q . (29) o v e r this v o l u m e yields a l i n e a r e q u a t i o n for e a c h e l e m e n t o f t h e r a t e of s t r a i n t e n s o r w r i t t e n 8

1

"<*> = κ Σ ^ ' S

(34)

1= 1

l

w h e r e / labels t h e vertices of t h e cell, a n d cß a r e coefficients w h i c h d e p e n d entirely o n t h e g e o m e t r y of t h e cell. l T h e g e o m e t r i c coefficients, cß, o f w h i c h t h e r e a r e t w e n t y - f o u r for e a c h cell, s u m m a r i z e in s t o r a b l e f o r m t h e g e o m e t r y of t h e m e s h . T h e y a p p e a r in all e q u a t i o n s w h i c h refer t o t h i s g e o m e t r y , a s in t h e c a l c u l a t i o n of t h e cell volume.

V=

ΣΛ'

/=i

* = 1,2, or 3,

(35)

a n d t h e c a l c u l a t i o n o f g r a d i e n t s of t h e r m o d y n a m i c v a r i a b l e s a t cell vertices. They are recomputed a n d stored only when the c o m p u t a t i o n mesh moves. T h e g e o m e t r i c coefficients a r e w r i t t e n in a c o m p a c t v e c t o r f o r m in t h e Appendix. 2. The Dynamical

Equations

T h e s y s t e m of difference e q u a t i o n s in s p a c e is c o m p l e t e d b y d e r i v i n g a n a p p r o x i m a t i o n t o E q . (9), t h e m o m e n t u m e q u a t i o n , b y g e n e r a l i z i n g G o a d ' s (1960) d e r i v a t i o n of difference e q u a t i o n s for t w o - d i m e n s i o n a l h y d r o d y n a m i c s . T h i s d e r i v a t i o n e m p l o y s t h e c o m p u t a t i o n of v i r t u a l w o r k a n d d ' A l e m b e r t ' s p r i n c i p l e t o f o r m a r e l a t i o n b e t w e e n c h a n g e s in t h e p o t e n t i a l e n e r g y o f t h e c o m p u t a t i o n cells a n d t h e c o r r e s p o n d i n g a c c e l e r a t i o n of t h e i r vertices. C o n s i d e r a s m a l l e l a b o r a t i o n of t h e d e s c r i p t i o n of t h e flow d e v e l o p e d in l t h e p r e v i o u s s e c t i o n . T o e a c h v e r t e x , t h e r e will n o w b e a s s i g n e d a m a s s , m , w h i c h is c o m p u t e d b y s o m e yet t o b e specified r u l e f r o m t h e d i s t r i b u t e d m a s s e s of t h e cells n e i g h b o r i n g v e r t e x /. I n a d d i t i o n , t h e p o t e n t i a l e n e r g y o f a cell is defined b y t h e e q u a t i o n W =

dVpw,

(36)

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

21

w h e r e w is t h e specific p o t e n t i a l e n e r g y w = i + Β · Β/2μρ.

(37)

A n e v o l u t i o n e q u a t i o n for w, c o r r e s p o n d i n g t o t h e e q u a t i o n for t h e i n t e r n a l e n e r g y , E q . (10), is given b y ρ dw/dt - ( Q · V) · u = 0.

(38)

T h i s e q u a t i o n is n o w u s e d t o c o m p u t e c h a n g e s in t h e cell p o t e n t i a l e n e r g y c o r r e s p o n d i n g t o d i s p l a c e m e n t s o f a v e r t e x l a b e l e d /. T h e c h a n g e i n W o v e r a t i m e i n t e r v a l , δί, is c o m p u t e d b y i n t e g r a t i n g E q . (38) o v e r t h e v o l u m e o f a cell, a n d o v e r δί, δψ

= jvdV

= W(t + δί) - W(t)

j^dt'QijSji.

(39)

T h e stress t e n s o r is c o n s t a n t w i t h i n a cell a n d m a y b e f a c t o r e d f r o m t h e v o l u m e i n t e g r a l . C o m b i n i n g E q s . (34) a n d (39) yields a n e q u a t i o n for t h e r e m a i n i n g v o l u m e i n t e g r a l for δΨ w r i t t e n

Σον)·

+ö

8W = \'t 'dt'[Q,

(40)

W h e n s o m e a v e r a g e v a l u e o v e r t h e t i m e i n t e r v a l δί is s u b s t i t u t e d for Qi3 a n d Cj in E q . (40), a n e q u a t i o n b e t w e e n c h a n g e s in W a n d d i s p l a c e m e n t s of t h e vertices r e s u l t s , a n d is w r i t t e n

w h e r e ds} is given b y ös! = j'^dt'uf.

(42)

I n w o r d s , E q . (41) says t h a t t h e c h a n g e in t h e p o t e n t i a l e n e r g y of a cell is e q u a l t o t h e n e g a t i v e of t h e p r o d u c t of t h e force e x e r t e d b y t h e cell o n e a c h v e r t e x a n d t h e d i s p l a c e m e n t of t h a t v e r t e x , bW = -

F

l

1

F · 5s ,

(43)

1= 1

where F

l

is t h e force e x e r t e d o n t h e v e r t e x / b y t h e cell. T h e force o n e a c h

22

JEREMIAH U. BRACKBILL

v e r t e x is r e l a t e d t o t h e a c c e l e r a t i o n e x p e r i e n c e d b y t h a t v e r t e x by d ' A l e m b e r t ' s equation, 1

l l

(Έ -m u )-Ss

l

= 0.

(44)

T h e r e f o r e , t h e c o m p o n e n t s of t h e a c c e l e r a t i o n o f v e r t e x / d u e t o stresses l a c t i n g i n t e r i o r t o a cell a r e given b y t h e r a t i o of bW t o t h e p r o d u c t of m a n d t h e c o m p o n e n t s of ös\ uj

=

1

-(l/mOftrfC, .

(45)

T h e r e s u l t a n t a c c e l e r a t i o n is c o m p u t e d b y s u m m i n g t h e c o n t r i b u t i o n s f r o m all e i g h t n e i g h b o r i n g cells, a n d is given b y

(46) l

w h e r e c f is t h e i t h c o m p o n e n t of t h e g e o m e t r i c coefficient a s s i g n e d t o t h e / t h v e r t e x of cell k. S o m e of t h e p r o p e r t i e s o f t h e e q u a t i o n of m o t i o n a r e easily s h o w n . F o r e x a m p l e , it c a n b e s h o w n t h a t w h e n t h e g e o m e t r i c coefficients a r e s u m m e d o v e r t h e vertices of a cell, /, t h e r e s u l t is z e r o . T h u s , t h e c o n t r i b u t i o n t o t h e c h a n g e in t h e m o m e n t u m o f p l a s m a f r o m stresses a c t i n g i n t e r i o r t o e a c h cell is z e r o , a n d m o m e n t u m is c o n s e r v e d . S i m i l a r l y , it is easily s h o w n t h a t w h e n t h e g e o m e t r i c coefficients a t a v e r t e x a r e s u m m e d o v e r t h e s u r r o u n d i n g cells, k, t h e r e s u l t is z e r o . T h u s , w h e n t h e stress is c o n s t a n t , a v e r t e x e x p e r i e n c e s n o acceleration.

3. Differencing

in Time and the Conservation

of

Energy

U n c o n d i t i o n a l l y s t a b l e difference e q u a t i o n s in t i m e a r e o b t a i n e d w h e n t h e e q u a t i o n s a r e m a d e implicit. W h e n t h e y a r e c e n t e r e d in t i m e , t h e m a g n e t i c flux a n d t h e t o t a l e n e r g y a r e c o n s e r v e d exactly. Since t h e e v o l u t i o n of i n t e r n a l e n e r g y is c a l c u l a t e d s e p a r a t e l y , t h e r e is n o s p u r i o u s e x c h a n g e of e n e r g y b e t w e e n t h e m a g n e t i c field a n d t h e p l a s m a , a n d l o c a l e n e r g y c o n s e r v a t i o n in t h e sense d e s c r i b e d b y F r o m m (1961) is o b t a i n e d . However, time-centered equations are n o t used. In the one-dimensional s t u d i e s w h i c h p r e c e d e d t h e f o r m u l a t i o n of t h e t w o - a n d t h r e e - d i m e n s i o n a l a l g o r i t h m s , it w a s f o u n d t h a t i m p l i c i t , first o r d e r e q u a t i o n s w e r e less dis­ persive t h a n s e c o n d o r d e r e q u a t i o n s , a n d n o t a p p r e c i a b l y m o r e diffusive. T o o b t a i n e x a c t e n e r g y c o n s e r v a t i o n w i t h t h e first o r d e r e q u a t i o n s , E q . (38)

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

23

for t h e p o t e n t i a l e n e r g y is s o l v e d r a t h e r t h a n E q . (31) for t h e p l a s m a i n t e r n a l e n e r g y . T h e f o r w a r d b i a s e d t i m e a n d s p a c e difference a p p r o x i m a t i o n s for t h e continuity, potential energy, a n d magnetic induction equations are written n + ni

p

» + i p _ »p + n+ip(n+ iw

_

+

\SHy

_ Q^+^\sJty

v)

δί = 0,

(47)

δί = 0,

(48)

and n

+

n

n+1

n+

% - Ba + ( Ba \SHy

- "

+ 1

A l"

+ 1

< S f a» δί = 0,

(49)

w h e r e t h e c o m p o n e n t s of t h e stress t e n s o r a r e given b y μζ)αβ

n 2

n+ll2

n+ n

η +ί

= +V Ba Bß

- (

Ρμ + % B^

öocß.

(50)

T h e t i m e level of e a c h t e r m in t h e s e e q u a t i o n s is d e n o t e d b y t h e l e a d i n g s u p e r s c r i p t , w i t h η d e n o t i n g t i m e level ί a n d η + 1 t i m e level ί + δί. T h e t i m e levels of t h e r a t e of s t r a i n t e n s o r c o r r e s p o n d t o t h e t i m e levels a s s i g n e d t o t h e velocities a p p e a r i n g in t h e m . T i m e - c e n t e r e d t e r m s a r e d e n o t e d b y η + \, a n d a r e e v a l u a t e d b y l i n e a r a v e r a g i n g b e t w e e n t i m e level η a n d n+\. The momen­ t u m e q u a t i o n , E q . (46), is a p p r o x i m a t e d b y +

»V

- V

c

= "Σ

(Q*i ?)

^Im\

(51)

k

w h e r e t h e v e r t e x m a s s , m\ is e q u a l t o t h e l i n e a r a v e r a g e o f t h e m a s s e s o f t h e e i g h t n e i g h b o r i n g cells. T o t a l e n e r g y is e x a c t l y c o n s e r v e d b e c a u s e t h e s u m o v e r cells o f t h e c h a n g e in p o t e n t i a l e n e r g y , δΨ, is e q u a l t o t h e s u m o v e r vertices of t h e c h a n g e in k i n e t i c e n e r g y , δΚΕ. C o n s i d e r t h e c h a n g e in t h e k i n e t i c e n e r g y o f v e r t e x /, given b y 1

δΚΕ

=

n+1/2 l

k

uj YJ(Q jic?)öi.

(52)

k

T h e s u m of t h e a s s o c i a t e d c h a n g e s in t h e p o t e n t i a l e n e r g y in t h e e i g h t n e i g h ­ b o r i n g cells is given b y

k

T h e s u m of δΨ a n d δΚΕ SW+

δΚΕ

=

is given b y

X5$(-" V i* + V +1/

k

c

kn+1/

k

\Sv> )

δ* = 0,

(54)

24

JEREMIAH U. BRACKBILL

w h e r e E q . (34) h a s b e e n u s e d . By c o n s t r u c t i o n , t o t a l e n e r g y is e x a c t l y c o n ­ served over the c o m p u t a t i o n mesh. T h e i n t e r n a l e n e r g y is c o m p u t e d b y s u b t r a c t i n g t h e m a g n e t i c field e n e r g y from the potential energy, i

n +

i

Β)/μρ].

=

(55)

B e c a u s e t h e p o s i t i o n o f a v e r t e x is a d v a n c e d w i t h t h e f o l l o w i n g first o r d e r equation +

» ixJ

=

"χ' + " α

+

ν<5ί,

(56) 2

r a t h e r t h a n a t i m e - c e n t e r e d e q u a t i o n , e r r o r s o f o r d e r dt a r e b e i n g a b s o r b e d by the internal energy. W e r e this t o p r o v e t r o u b l e s o m e , Eqs. (47)-(49) a n d (55) c o u l d b e c e n t e r e d in t i m e t o r e m o v e t h i s e r r o r . T h e t i m e levels in t h e stress t e n s o r , given b y E q . (50), a r e a l r e a d y a p p r o p r i a t e l y c h o s e n . 4. The Solution

of the Implicit

Equations

of

Motion

I n i n c o m p r e s s i b l e flow, t h e fluid v e l o c i t y is divergence-free. W h e n V · u is initially z e r o , a p r e s s u r e field w h i c h a s s u r e s t h a t V · u r e m a i n s z e r o is a s o l u ­ tion of Poisson's equation, ΨΡ'

= R,

(57)

w h e r e P' is e q u a l t o P/p a n d R c o n t a i n s b o d y forces a n d p e r h a p s v i s c o u s stresses. I n a m e t h o d d e v e l o p e d b y H a r l o w a n d W e l c h (1965), t h e e q u a t i o n is solved for P' for e a c h cycle b e f o r e s o l v i n g a n explicit difference a p p r o x i ­ m a t i o n t o t h e m o m e n t u m e q u a t i o n . T h e r e s u l t i n g velocity field is d i v e r g e n c e free. I t w a s l a t e r n o t e d ( H a r l o w a n d A m s d e n , 1971) t h a t a relatively m i n o r m o d i f i c a t i o n o f t h e i n c o m p r e s s i b l e flow e q u a t i o n s r e s u l t s in e q u a t i o n s w h i c h c a n b e u s e d for flow a t all s p e e d s . T h e m e t h o d u s i n g t h e s e e q u a t i o n s t r e a t s t h e d e n s i t y in t h e e q u a t i o n o f s t a t e a n d t h e d e n s i t y a n d velocity in t h e c o n ­ t i n u i t y e q u a t i o n s implicitly. F o r i n c o m p r e s s i b l e flow, t h e r e s u l t i n g e q u a t i o n s r e d u c e t o t h o s e u s e d for i n c o m p r e s s i b l e flow. F o r c o m p r e s s i b l e flow, t h e equations are simply implicit formulations of the usual Eulerian equations. I n a v a r i a t i o n o f t h e o r i g i n a l E u l e r i a n m e t h o d for flow a t all s p e e d s ( H i r t et al., 1974), t h e s o l u t i o n of P o i s s o n ' s e q u a t i o n for t h e p r e s s u r e is o b t a i n e d b y iteratively s o l v i n g t h e c o n t i n u i t y a n d m o m e n t u m e q u a t i o n s s i m u l t a n e o u s l y . T h i s m e t h o d is c o n v e n i e n t l y a p p l i e d t o L a g r a n g i a n difference e q u a t i o n s . I t is a l s o easily e x t e n d e d t o M H D flow.

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

25

I n t h e e x t e n s i o n o f t h e m e t h o d t o M H D , t h e m a g n e t i c field is t r e a t e d i m p l i c i t l y in t h e i n d u c t i o n a n d m o m e n t u m e q u a t i o n s . A c o n s i s t e n t s o l u t i o n o f t h e c o n t i n u i t y , i n d u c t i o n , a n d m o m e n t u m e q u a t i o n s is o b t a i n e d iteratively b y a n a d a p t a t i o n o f t h e successive o v e r r e l a x a t i o n s c h e m e d e v e l o p e d b y H i r t et al. (1974). I n t h i s s c h e m e , t h e i t e r a t i o n v a r i a b l e s a r e t h e d e n s i t y a n d m a g n e t i c field i n t e n s i t y . W i t h i n a cell, e a c h of t h e s e v a r i a b l e s is a d j u s t e d t o r e d u c e t h e r e s i d u a l e r r o r in t h e s o l u t i o n o f t h e c o n t i n u i t y a n d i n d u c t i o n e q u a t i o n s , E q s . (7) a n d (8). S u b s e q u e n t l y , c h a n g e s in t h e velocity a t t h e v e r t i c e s o f t h e cell d u e t o c h a n g e s in t h e stress t e n s o r a r e c o m p u t e d f r o m t h e e q u a t i o n s o f m o t i o n . T h e c h a n g e s in t h e v e l o c i t y a t t h e vertices o f a cell aifect t h e r e s i d u a l e r r o r in n e i g h b o r i n g cells t h r o u g h t h e r a t e o f s t r a i n t e n s o r s o t h a t t h e r e s i d u a l e r r o r p r o p a g a t e s f r o m cell t o cell d u r i n g a n i t e r a t i o n p a s s a s in a G a u s s - S e i d e l iteration. T h e a d j u s t m e n t w i t h i n e a c h cell is p e r f o r m e d b y m e a n s o f a N e w t o n R a p h s o n iteration. T h e derivatives of the residual errors with respect t o the iteration variables are c o m p u t e d from analytical expressions obtained by formally differentiating the e q u a t i o n s of m o t i o n , a n d the complete system of f o u r s i m u l t a n e o u s e q u a t i o n s is s o l v e d . T h e d e r i v a t i v e s i n c l u d e n o t o n l y t h e direct variation of the residuals with the iteration variables, b u t also the i m p l i c i t v a r i a t i o n s t h r o u g h t h e r a t e o f s t r a i n t e n s o r . A s is d i s c u s s e d b y H i r t et al. (1974), t h e i n c l u s i o n o f t h e i m p l i c i t v a r i a t i o n m a k e s t h e i t e r a t i o n s t a b l e for all signal s p e e d s b e c a u s e it c a u s e s t h e a d j u s t m e n t o f t h e i t e r a t i o n v a r i a b l e t o be b o u n d e d . W h e n a n u m b e r of e q u a t i o n s are being solved simultaneously, it is n e c e s s a r y for t h e s a m e r e a s o n t o c o m p u t e all t e r m s in t h e J a c o b i a n m a t r i x . E v e n t h o u g h t h e i m p l i c i t e q u a t i o n s a r e s t a b l e for a l a r g e r t i m e - s t e p , in p r a c t i c e , t h e t i m e - s t e p is l i m i t e d t o o n e w h i c h satisfies t h e i n e q u a l i t y , \u\St
0 < / < l / 4 .

(58)

W h e n t h i s i n e q u a l i t y is satisfied, t h e r e l a t i v e c h a n g e s in t h e t h e r m o d y n a m i c q u a n t i t i e s o v e r a t i m e - s t e p , St, a r e o f t h e o r d e r o f / . W h e n / is 0.25 o r less, a b o u t t e n i t e r a t i o n s a r e all t h a t a r e n e c e s s a r y for t h e r e l a t i v e e r r o r in t h e solu­ 5 t i o n of t h e c o n t i n u i t y a n d i n d u c t i o n e q u a t i o n s t o b e r e d u c e d t o o n e p a r t in 1 0 for l o w s p e e d , c o m p r e s s i b l e flows.

5. Boundary

Conditions

W h e n t h e b o u n d a r y is a c l o s e d c o n d u c t i n g s u r f a c e a n d t h e n o r m a l c o m ­ p o n e n t o f t h e m a g n e t i c field a t t h e b o u n d a r y is z e r o , t h e b o u n d a r y c o n d i t i o n is given b y u

ή < 0,

(59)

26

JEREMIAH U. BRACKBILL

w h e r e ft is t h e o u t w a r d d i r e c t e d n o r m a l t o t h e s u r f a c e a n d u is t h e v e l o c i t y o f a b o u n d a r y v e r t e x . I n g e n e r a l , h o w e v e r , it is n o t sufficient t o s i m p l y r e a d j u s t t h e n o r m a l v e l o c i t y t o satisfy E q . (59), for t w o r e a s o n s . F i r s t , w h e n stresses o n a b o u n d a r y a r e c o m p u t e d f r o m o n e side, fictitious t a n g e n t i a l stresses d u e t o t r u n c a t i o n errors are exerted o n b o u n d a r y vertices. Second, the e q u a t i o n o f m o t i o n is o n l y first o r d e r a c c u r a t e i n δί so t h a t a b o u n d a r y v e r t e x drifts 2

off t h e b o u n d a r y b y a n a m o u n t o f o r d e r δί

at each time-step. A general

b o u n d a r y t r e a t m e n t is d e r i v e d w h i c h is a n e x t e n s i o n o f free-surface b o u n d a r y c o n d i t i o n s ( H i r t ei al., 1970) t o a free-slip b o u n d a r y . T h e b o u n d a r y o f t h e p l a s m a is t r e a t e d a s a free s u r f a c e u p o n w h i c h t h e c o n d u c t i n g b o u n d a r y e x e r t s a force o f c o n s t r a i n t . T h e a p p r o p r i a t e b o u n d a r y c o n d i t i o n for t h e stress a t a free s u r f a c e is given by Q · ft = 0.

(60)

I n o r d e r t o satisfy E q . (60), t h e m e s h is e x t e n d e d a s d e s c r i b e d by H i r t ei al. (1970), a n d t h e stress, Q , is e x t r a p o l a t e d b e y o n d t h e b o u n d a r y . T h e e x t r a ­ p o l a t i o n o f t h e field c o m p o n e n t s is s u c h t h a t Β · ft o n t h e b o u n d a r y is z e r o , a n d t h a t Β χ ft e x t e r i o r t o t h e m e s h i n c l u d e s g e o m e t r i c c o r r e c t i o n s for t h e curvature of the b o u n d a r y . N e x t , t h e c o n d i t i o n o n t h e n o r m a l velocity, E q . (59), is satisfied b y c o m ­ p u t i n g a n o r m a l force o f c o n s t r a i n t . W h e r e t h e e q u a t i o n for t h e b o u n d a r y is g i v e n b y / ( x ) = 0,

(61)

t h e d i s p l a c e m e n t o f a b o u n d a r y v e r t e x , ηλ, d u e t o t h e a c t i o n o f t h e force o f c o n s t r a i n t , is g i v e n b y t h e s o l u t i o n o f t h e e q u a t i o n /(x +

m + λή) =

0,

(62)

w h e r e u i n c l u d e s c o n t r i b u t i o n s f r o m all t h e stresses a c t i n g i n t e r i o r t o t h e b o u n d a r y . F o r e x a m p l e , w h e n t h e b o u n d a r y is a r i g h t c i r c u l a r c y l i n d e r , λ is given by the e q u a t i o n 2

λ

+ 2(ft · ü — Riot)λ

+ [ ü · ü2R(n

· fi)]/<5f = 0,

(63)

w h e r e λ < 0 c o n s i s t e n t w i t h E q . (59). T h e s m a l l e r m a g n i t u d e r o o t is a l w a y s chosen. T h e larger r o o t corresponds to moving the vertex to a diametrically o p p o s e d p o i n t o n the surface.

B . THE REZONE PHASE OF A GENERALIZED MESH CALCULATION T h e t r a n s p o r t of m a s s , m o m e n t u m , m a g n e t i c flux, a n d e n e r g y f r o m cell t o cell is c o m p u t e d f r o m E q s . ( 1 5 ) - ( 1 8 ) . T h e difference e q u a t i o n s for t h e c o n ­ v e c t i o n p h a s e i n t w o d i m e n s i o n s a r e g i v e n b y H i r t et al. (1974) a n d P r a c h t

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

27

(1975). T h e e x t e n s i o n t o i n c l u d e t h e t r a n s p o r t o f m a g n e t i c flux is d e s c r i b e d in B r a c k b i l l a n d P r a c h t (1973). I n a d d i t i o n , a d e s c r i p t i o n o f a t w o - d i m e n s i o n a l c o m p u t e r c o d e is g i v e n in A m s d e n a n d H i r t (1973). W i t h t h i s n u m b e r o f d e t a i l e d references f o r w h a t is essentially a s t a n d a r d a n d w e l l - k n o w n m e t h o d for t h e c o m p u t a t i o n o f c o n v e c t i v e t r a n s p o r t , t h e r e is n o n e e d for a n o t h e r description here. W h e n t h e g r i d velocity, u', is z e r o s o t h a t t h e m e s h is E u l e r i a n , t h e g e n ­ e r a l i z e d m e s h m e t h o d is p r e y t o all t h e ills w h i c h b e s e t a n y o t h e r E u l e r i a n c a l c u l a t i o n . T h a t is, t h e n u m e r i c a l diffusion i n t r o d u c e d b y t h e d o n o r cell e q u a t i o n seriously reduces the overall accuracy of the calculation. A n appli­ c a t i o n o f t h e t r u n c a t i o n e r r o r c o r r e c t i o n p r o c e d u r e d e v e l o p e d b y R i v a r d et al (1973) t o t h e c a l c u l a t i o n o f t h e c o n v e c t i v e t r a n s p o r t s h o u l d significantly i m p r o v e the accuracy of the rezone phase. A p p l y i n g this p r o c e d u r e t o t w o a n d t h r e e - d i m e n s i o n a l M H D c a l c u l a t i o n s is relatively s t r a i g h t f o r w a r d . O u r u s u a l w a y o f d e a l i n g w i t h t h e a c c u r a c y p r o b l e m in r e z o n i n g is t o d o as little o f it a s p o s s i b l e . B y r e d u c i n g t h e r e l a t i v e velocity b e t w e e n t h e p l a s m a a n d t h e m e s h , t h e t r u n c a t i o n e r r o r s in all t h e a p p r o x i m a t i o n s t o t h e c o n v e c t i v e d e r i v a t i v e listed in S e c t i o n I I I a r e r e d u c e d . T h u s , a n a l m o s t - L a g r a n g i a n p r e s c r i p t i o n for t h e g r i d v e l o c i t y , t h a t is, a p r e s c r i p t i o n in w h i c h u — u' is small everywhere, reduces t r a n s p o r t t o the point where the errors introduced b y its a p p r o x i m a t i o n a r e insignificant. A n a l m o s t - L a g r a n g i a n p r e s c r i p t i o n u s i n g t h e s m o o t h i n g a l g o r i t h m is o u t l i n e d in t h e f o l l o w i n g s e c t i o n . 1. The Almost-Lagrangian

Mesh

A n a l m o s t - L a g r a n g i a n calculation exploits the ability of a L a g r a n g i a n m e s h t o b e d i s t o r t e d w i t h o u t affecting t h e a c c u r a c y o f t h e g e n e r a l i z e d differ­ ence e q u a t i o n s . S o m e consideration of the c o o r d i n a t e t r a n s f o r m a t i o n given b y E q . (27) will c o n v i n c e t h e r e a d e r t h a t t h e f o r m a l a c c u r a c y o f t h e difference e q u a t i o n s is p r e s e r v e d a s l o n g a s e v e r y cell in t h e m e s h in c o n v e x . T h e f o r m a l a c c u r a c y is l o s t w h e n a cell b e c o m e s c o n c a v e , for a c o n c a v e cell h a s n o i m a g e in n a t u r a l c o o r d i n a t e s . T h u s , t h e d u t y o f t h e r e z o n e in a n a l m o s t - L a g r a n g i a n c a l c u l a t i o n is t o p r e v e n t c o n c a v e cells. F o r m a n y flows, t h i s r e q u i r e s little relative m o t i o n between the p l a s m a a n d the mesh, provided a n a p p r o p r i a t e m e s h p r e s c r i p t i o n is u s e d . O n e s u c h p r e s c r i p t i o n is t h e s m o o t h i n g a l g o r i t h m , w h i c h is o u t l i n e d in t h e n e x t p a r a g r a p h . 2. The Smoothing

Algorithm

O u r m o s t successful p r e s c r i p t i o n for a n a l m o s t - L a g r a n g i a n c a l c u l a t i o n is t h e s m o o t h i n g a l g o r i t h i m , w h i c h w a s first s u g g e s t e d t o u s b y P . B r o w n e ( u n p u b l i s h e d n o t e s , 1972) a n d is d e s c r i b e d in B r a c k b i l l a n d P r a c h t (1973).

28

JEREMIAH U. BRACKBILL

T h i s a l g o r i t h m is d e r i v e d b y m i n i m i z i n g a n i n t e g r a l / , given b y

/ = I dV'imf

+W

2

+ (Vv) ],

(64)

over the c o m p u t a t i o n mesh. T h e Euler equation corresponding to this mini­ m i z a t i o n p r i n c i p l e c a n b e s o l v e d directly b y i n t e r c h a n g i n g t h e d e p e n d e n t a n d i n d e p e n d e n t v a r i a b l e s ( T h o m p s o n et al, 1974). O u r a p p r o a c h is t o a p p l y t h e finite-element m e t h o d . T h e i n t e g r a l / is e v a l u a t e d b y c o n s t r u c t i n g e l e m e n t s in w h i c h t h e n a t u r a l c o o r d i n a t e s , ξ, η, a n d ν, a r e piecewise l i n e a r f u n c t i o n s o f t h e p h y s i c a l c o o r d i n a t e s , x, y, a n d z. T h e i n t e g r a l is t h e n m i n i m i z e d b y t r e a t i n g t h e p o s i t i o n o f e a c h m e s h p o i n t a s a p a r a m e t e r . T h e r e s u l t o f a single s w e e p o v e r t h e m e s h is a t r i a l v a l u e for t h e d i s p l a c e m e n t o f e a c h v e r t e x , δχ, f r o m t h e p o s i t i o n a t w h i c h / w o u l d b e m i n i m i z e d . T h e essence o f t h e a l m o s t - L a g r a n g i a n m e s h is t h a t t h i s d i s p l a c e ­ m e n t is r e d u c e d b y a f r a c t i o n a l a m o u n t , ω, e a c h cycle, given b y ω = δί/τ,

(65)

w h e r e τ is a m e s h r e l a x a t i o n t i m e , w h i c h is t y p i c a l l y t e n t i m e - s t e p s o r m o r e . T h u s , a t e a c h m e s h p o i n t , t h e g r i d v e l o c i t y t o b e s u b s t i t u t e d i n t o E q s . (15)—(18) t o c o m p u t e t r a n s p o r t is given b y u' = δχ/τ.

(66)

T h i s velocity is s m a l l c o m p a r e d w i t h t h e velocity d e t e r m i n i n g t h e t i m e - s t e p w h e n τ is l a r g e c o m p a r e d w i t h δί (cf. S e c t i o n I V , A , 4). T h e effect of t h e s m o o t h i n g a l g o r i t h m o n a d i s t o r t e d m e s h is b e s t s h o w n b y a n e x a m p l e . I n F i g . 1, t h e a p p l i c a t i o n o f t h e a l g o r i t h m t o a d i s t o r t e d ,

(a)

(b)

(c)

FIG. 1. T h e effect o f the s m o o t h i n g algorithm o n a distorted Lagrangian mesh can be seen in these plots o f a computation mesh. T h e Lagrangian m e s h is s h o w n in 1(a), and the mesh after 10 and 50 applications of the s m o o t h i n g algorithm is s h o w n in 1(b) and 1(c). All of the cells are already convex in the mesh s h o w n in 1(b), demonstrating that a Lagrangian mesh is transformed into an acceptable mesh by only a small amount of rezoning.

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

29

L a g r a n g i a n m e s h is s h o w n . I n F i g . l a , t h e m e s h is clearly p a t h o l o g i c a l . S o m e z o n e s e v e n o v e r l a p . A f t e r 10 i t e r a t i o n s o f t h e a l g o r i t h m w i t h ω = 0 . 2 5 , e v e r y cell in t h e m e s h s h o w n i n F i g . l b is c o n v e x . F i n a l l y , after 50 i t e r a t i o n s , t h e m e s h a p p e a r s in F i g . l c t o h a v e c o n v e r g e d t o a s y m m e t r i c , p o l a r m e s h .

V . Applications T h e results of calculations using the generalized m e s h m e t h o d are pre­ s e n t e d . T h e s e c a l c u l a t i o n s d e m o n s t r a t e h o w t h e flexibility o f t h e m e t h o d c a n b e e x p l o i t e d . F o r e x a m p l e , in t h e c a l c u l a t i o n o f t h e i m p l o s i o n o f a n a x i ­ s y m m e t r i c t h e t a p i n c h , t h e L a g r a n g i a n c h a r a c t e r o f t h e m e s h is u s e d t o r e s o l v e t h e b o u n d a r y b e t w e e n t h e p l a s m a a n d t h e v a c u u m , a n d t h u s t o define t h e b o u n d a r i e s o f t h e r e g i o n o v e r w h i c h t h e v a c u u m field is c o m p u t e d . I n a d d i t i o n , t h e r e s u l t s s h o w t h e a c c u r a c y w h i c h is a c h i e v e d b y p r e s e n t i n g a d i r e c t c o m ­ parison between linear theory a n d a nonlinear calculation of the rotating theta pinch.

A. A SHARP BOUNDARY CALCULATION OF THE THETA PINCH T h e c o m p r e s s i o n o f a t h e t a p i n c h is d r i v e n b y a n a p p l i e d a z i m u t h a l c u r r e n t o n a c o n d u c t i n g cylindrical wall s u r r o u n d i n g the plasma. T h e current p r o ­ d u c e s a m a g n e t i c field a t t h e w a l l w h i c h d r i v e s t h e p l a s m a t o w a r d t h e axis o f s y m m e t r y . S o m e m a g n e t i c field diffuses i n t o t h e p l a s m a , a n d s o m e p l a s m a is n o t p i c k e d u p b y t h e i n c o m i n g field. T h a t is, t h e b o u n d a r y b e t w e e n t h e p l a s m a r e g i o n a n d t h e field r e g i o n is n o t s h a r p . N e v e r t h e l e s s , t h e i n e r t i a l p r o p e r t i e s o f t h e p l a s m a left b e h i n d a r e n o t i m p o r t a n t , n o r is a n y c u r r e n t c a r r i e d b y t h i s low density plasma. T h u s , the sharp b o u n d a r y model mentioned by R o b e r t s a n d P o t t e r (1968) c a n b e u s e d . I n t h i s m o d e l , t h e r e is a n i n t e r f a c e s e p a r a t i n g t h e p l a s m a a n d t h e v a c u u m . I n t h e v a c u u m , t h e e q u a t i o n for t h e a z i m u t h a l c o m p o n e n t o f t h e v e c t o r p o t e n t i a l ( L e w i s , 1966), is solved. I n t h e p l a s m a , t h e o r d i n a r y M H D e q u a t i o n s , i n c l u d i n g resistive diffusion b u t n o t m a s s diffusion, a r e s o l v e d . A t t h e i n t e r f a c e , t h e c o n t i n u i t y o f stress c o n d i t i o n a c r o s s t h e free s u r f a c e o f t h e p l a s m a is m a i n t a i n e d t o first o r d e r in St, t h e t i m e - s t e p . T h e e q u a t i o n f o r t h e v e c t o r p o t e n t i a l is s o l v e d b y a n a p p l i c a t i o n o f t h e finite-element m e t h o d . T h e i n t e g r a l t o b e m i n i m i z e d is w r i t t e n ( M o r s e a n d F e s h b a c k , 1953) (67) w h e r e r a n d ζ a r e t h e r a d i a l a n d a x i a l c o o r d i n a t e s . T h e i n t e g r a l is m i n i m i z e d with t h e v a l u e o f Αθ a t e a c h m e s h p o i n t a s a p a r a m e t e r . T h e r e s u l t i n g s y s t e m

JEREMIAH U. BRACKBILL

30

o f l i n e a r e q u a t i o n s is s o l v e d b y successive o v e r r e l a x a t i o n ( Y o u n g , 1962), b u t a n alternating direction implicit m e t h o d can also be used. T h e b o u n d a r y c o n d i t i o n s for t h e v a c u u m r e g i o n a r e t h e v a l u e o f Αθ o n t h e v a c u u m i n t e r f a c e , a n d e i t h e r Αθ o r dAe/dn, w h e r e η is t h e n o r m a l c o o r d i n a t e t o t h e w a l l , a t e a c h p o i n t o n t h e wall. T h e v e c t o r p o t e n t i a l a t t h e i n t e r f a c e is d e t e r m i n e d b y i n t e g r a t i n g t h e c h a n g e in Αθ a w a y f r o m t h e axis, w h e r e d(rAe) is g i v e n b y d(rAe)

= r (Bz dr - Br dz).

(68)

T h e field c o m p o n e n t s , Br a n d 2? z, a r e k n o w n in t h e i n t e r i o r of t h e p l a s m a . B e c a u s e t h e v a l u e o f t h e n o r m a l c o m p o n e n t o f t h e m a g n e t i c field o n e i t h e r side o f t h e i n t e r f a c e is d e t e r m i n e d b y t h e v a r i a t i o n o f Ae a l o n g t h e i n t e r f a c e , Β · Α is c o n t i n u o u s a c r o s s t h e i n t e r f a c e . S o m e r e s u l t s o f a s h a r p b o u n d a r y t h e t a p i n c h c a l c u l a t i o n a r e s h o w n in F i g . 2. I n t h e s e c a l c u l a t i o n s t h e c u r r e n t a t t h e w a l l is given b y , dAJdn

2

= (α + β c o s k z ) t + γ,

(69)

w h e r e α, β, a n d γ a r e c o n s t a n t s o f o r d e r u n i t y . O n e - h a l f t h e w a v e l e n g t h c o r r e s p o n d i n g t o k is c o m p u t e d , a n d p e r i o d i c b o u n d a r y c o n d i t i o n s a r e a p p l i e d a t ζ — 0 a n d ζ = π/k. Initially, t h e c o m p u t a t i o n m e s h a p p e a r s a s s h o w n in F i g . 2 a , w i t h l a r g e z o n e s r e s o l v i n g t h e p l a s m a a n d s m a l l z o n e s t h e v a c u u m . T h e axis of t h e p i n c h is o n t h e left a n d t h e c o n d u c t i n g w a l l is o n t h e r i g h t . I n F i g . 2 b , t h e p l a s m a is

(a)

(b)

(c)

(d)

F I G . 2. The results of a sharp boundary calculation of a theta pinch are s h o w n . The computation mesh and magnetic field lines at the initial time in 2(a) and 2(b), and just before m a x i m u m compression in 2(c) and 2(d), show the plasma, depicted by dots in 2(b) and 2(d), being compressed by the magnetic field entering o n the right in each frame o n t o the axis of symmetry o n the left. The leftmost magnetic field line in 2(d) nearly coincides with the z o n e line in 2(c) which delineates the v a c u u m - p l a s m a interface.

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS

31

d e n o t e d b y m a r k e r p a r t i c l e s , a n d flux surfaces b y c o n t o u r s in t h e v a c u u m r e g i o n . A t a l a t e r t i m e , t h e c o m p u t a t i o n m e s h a p p e a r s a s s h o w n in F i g . 2c. T h e i n t e r f a c e , w h i c h c o i n c i d e s w i t h t h e left b o u n d a r y o f t h e fifth z o n e f r o m t h e r i g h t , h a s m o v e d t o t h e left. T h e m e s h h a s in fact f o l l o w e d t h e p l a s m a , s h o w n a t t h e c o r r e s p o n d i n g t i m e in F i g . 2 d , a s it is c o m p r e s s e d b y t h e m a g n e t i c field. ( T h e v a l u e o f m a g n e t i c flux a s s o c i a t e d w i t h e a c h c o n t o u r is c o n s t a n t . ) T h e a d v a n t a g e s o f a g e n e r a l i z e d L a g r a n g i a n m e s h a r e e v i d e n t in t h i s c a l c u l a t i o n . F i r s t , t h e i n t e r f a c e is r e s o l v e d b y t h e m e s h s i m p l y b y c a u s i n g p o i n t s o n t h e i n t e r f a c e t o m o v e w i t h t h e l o c a l fluid velocity. S e c o n d , t h e p l a s m a is r e s o l v e d a t t h e final t i m e e v e n a s it o c c u p i e s a s m a l l e r a n d s m a l l e r f r a c t i o n o f t h e d o m a i n , b e c a u s e it c a r r i e s its z o n e s w i t h it a s it m o v e s . I t w o u l d n o t b e i m p o s s i b l e t o u s e fluid m a r k e r s ( A m s d e n a n d H a r l o w , 1970) t o r e s o l v e t h e i n t e r f a c e , b u t it w o u l d b e difficult t o d o it a s easily o r a s well.

B . THE ROTATING THETA PINCH T h e r o t a t i n g t h e t a p i n c h is initially in u n i f o r m r o t a t i o n w i t h i n a c y l i n d r i c a l , c o n d u c t i n g shell. A o n e - d i m e n s i o n a l e q u i l i b r i u m is g i v e n b y t h e e q u a t i o n 2

( ,

(70)

w h e r e Ω is t h e a n g u l a r velocity o f t h e p l a s m a , a n d Β is t h e a x i a l field. F o r 0 ^ \k\ < kc a n d m^O, where k a n d m are the axial a n d a z i m u t h a l waven u m b e r s , r e s p e c t i v e l y , t h e e q u i l i b r i u m is u n s t a b l e t o p e r t u r b a t i o n s in r a d i a l position of the form i(r)

= ^ ^ η

K k

nc

o s ( k z

+ mO).

(71)

W h e n £ ( r ) is t h e r a d i a l e i g e n f u n c t i o n g i v e n b y l i n e a r i z i n g a n d s o l v i n g t h e M H D e q u a t i o n s , t h e v a r i a t i o n in t i m e o f t h e e i g e n f u n c t i o n is given b y

ξ(ί)

= i(0)exp(yO,

w h e r e γ is t h e g r o w t h r a t e o f t h e i n s t a b i l i t y . T h e d o m a i n is p e r i o d i c in Θ a n d z, a n d is b o u n d e d a t r e q u a l t o R c o n d u c t i n g , c y l i n d r i c a l shell. W h e n k e q u a l s z e r o , g r a d i e n t s in t h e a x i a l d i r e c t i o n a r e z e r o a n d t h e is c o n f i n e d t o t h e r-θ p l a n e . A t h r e e - d i m e n s i o n a l c o m p u t a t i o n o f t h e d i m e n s i o n a l flow is p e r f o r m e d b y r e d u c i n g t h e n u m b e r o f z o n e s in t h e direction to one.

(72)

by a flow twoaxial

A detailed c o m p a r i s o n between linear theory a n d the results of the c o m ­ p u t a t i o n is m a d e p o s s i b l e b y t h e u s e o f t r a c e p a r t i c l e s ( P r a c h t , 1975), w h i c h

32

JEREMIAH U. BRACKBILL

m o v e w i t h t h e l o c a l fluid velocity s o t h a t t h e i r t r a j e c t o r y is a l s o t h e t r a j e c t o r y o f t h e fluid e l e m e n t w i t h i n w h i c h t h e y lie. I n stability c a l c u l a t i o n s , t r a c e p a r t i c l e s a r e u s e d t o f o r m a s h a d o w g r i d o f L a g r a n g i a n c o o r d i n a t e s w h i c h is initially i d e n t i c a l w i t h t h e c o m p u t a t i o n m e s h b u t w h i c h is n o l o n g e r i d e n t i c a l w h e n r e l a t i v e m o t i o n b e t w e e n t h e fluid a n d t h e m e s h o c c u r s b e c a u s e o f r e z o n i n g . If R0 is t h e initial r a d i a l c o o r d i n a t e o f a t r a c e p a r t i c l e i n c y l i n d r i c a l c o o r d i n a t e s , its p o s i t i o n a t a n y t i m e t c a n b e w r i t t e n 0, z ) = R0 + Σ ξ κ m e x p [i(kz + mQ +

R(R09

J],

(73)

k, m

w h e r e ξΚηι a n d
("*)

1

/2

(74)

and 1

Φ = tan " [Im (/)/Re (/)],

(75)

w h e r e / is g i v e n b y

7=

i

ΙοΊο

αΘ

dzR R

t ( o,0,z)

- R 0 ] . t cxpl-i(kz

+ m9)l

(76)

T h e ξΚτη c o m p u t e d f r o m E q . (74) c a n b e c o m p a r e d w i t h t h e ξΚηχ o f E q . (72), w h i c h a r e p r e d i c t e d b y l i n e a r t h e o r y . I n F i g . 3 , t h e g r o w t h o f ξ 0 > 3 in t i m e is s h o w n f o r t h r e e different c o m p u t a t i o n s . C u r v e 1, t h e l o w e s t c u r v e , is t h e r e s u l t o f a n E u l e r i a n c a l c u l a t i o n w i t h d o n o r cell differencing o f t h e c o n v e c t i o n t e r m s . D o n o r cell differencing is a c c u r a t e o n l y t o first o r d e r a n d e x t r e m e l y diffusive, a n d is clearly i n a d e q u a t e for stability c a l c u l a t i o n s w h e r e flow is o c c u r r i n g , f o r it r a d i c a l l y r e d u c e s t h e growth rate from that obtained with m o r e accurate calculations. F o r example, the results displayed in curve n u m b e r t w o , t h e middle curve, were obtained by u s i n g i n t e r p o l a t e d d o n o r cell differencing o f t h e c o n v e c t i o n t e r m s ( A m s d e n a n d H i r t , 1973), a s e c o n d o r d e r a c c u r a t e s c h e m e . B e y o n d o n e - h a l f r o t a t i o n p e r i o d , t h e g r o w t h r a t e is n e a r l y c o n s t a n t a n d c o m p a r a b l e t o t h e g r o w t h r a t e m e a s u r e d f r o m t h e u p p e r c u r v e , w h i c h is t h e r e s u l t o f a n a l m o s t - L a g r a n g i a n calculation. In the almost-Lagrangian calculation, the mesh relaxation time, τ, is a p p r o x i m a t e l y 2§bt a n d a d o n o r cell a p p r o x i m a t i o n t o t h e c o n v e c t i o n t e r m s is e m p l o y e d . A l i n e a r g r o w t h p h a s e c a n b e seen b e t w e e n 0.25 a n d 1.2 r o t a t i o n p e r i o d s , w h e r e t h e g r o w t h is e x p o n e n t i a l w i t h a very slowly d e c l i n i n g g r o w t h

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS 3 3

F I G . 3. The growth of a perturbation of an unstable equilibrium in a rotating theta pinch is plotted as c o m p u t e d with Eulerian and almost-Lagrangian computation meshes.

r a t e . S a t u r a t i o n c a n b e seen after 1.2 r o t a t i o n p e r i o d s w h e n t h e g r o w t h r a t e rapidly declines t o zero. T h e r a d i a l e i g e n f u n c t i o n a t o n e r o t a t i o n p e r i o d is s h o w n in F i g . 4. T h e c u r v e f r o m a n a l m o s t - L a g r a n g i a n c a l c u l a t i o n is s o m e w h a t b r o a d e r b u t s i m i l a r in s h a p e t o t h e initial p e r t u r b a t i o n . O n t h e o t h e r h a n d , t h e r a d i a l e i g e n f u n c t i o n for t h e E u l e r i a n c a l c u l a t i o n h a s d e v e l o p e d a d o u b l e - h u m p e d s t r u c t u r e w h i c h suggests t h a t d i s p e r s i o n e r r o r s h a v e s p a t i a l l y s e p a r a t e d m o d e s c o n t a i n e d in t h e initial e i g e n f u n c t i o n . T h e E u l e r i a n c a l c u l a t i o n is t e r m i n a t e d b y a n in­ stability o f t h e t y p e d e s c r i b e d in S e c t i o n I I I , a t a b o u t o n e r o t a t i o n p e r i o d . A n essential difficulty w i t h t h r e e - d i m e n s i o n a l c a l c u l a t i o n s is t h a t t h e a m o u n t o f i n f o r m a t i o n it is p o s s i b l e t o s t o r e is o n l y e n o u g h t o finely resolve p o r t i o n s o f t h e flow. I n t h e c a s e o f t h e m = 3 p e r t u r b a t i o n , t h e r a d i a l eigen­ f u n c t i o n is m o s t efficiently r e s o l v e d b y c l u s t e r i n g t h e p o i n t s r a d i a l l y a b o u t r = 1. W h e n t h e a l m o s t - L a g r a n g i a n m e s h is u s e d , t h e e i g e n f u n c t i o n is r e s o l v e d a t all t i m e s b e c a u s e t h e m e s h p o i n t s m o v e w i t h t h e flow. By c o n t r a s t , t h e E u l e r i a n m e s h p r o v i d e s fine r e s o l u t i o n of a p o r t i o n o f t h e flow o n l y u n t i l it m o v e s o u t o f t h a t r e g i o n o f t h e m e s h w h i c h c a n p r o v i d e it. I n g e n e r a l , o p t i m u m o v e r a l l r e s o l u t i o n w i t h a n E u l e r i a n m e s h is given b y a u n i f o r m d i s t r i b u t i o n

34

JEREMIAH U. BRACKBILL

— Theoretical >— Eulerian — ALE

Radius (Plasma radii) F I G . 4. The computed and linear theoretical perturbation amplitudes for a rotating theta pinch are plotted at a time corresponding to one rotation period of the pinch.

o f p o i n t s in s p a c e e v e n t h o u g h s u c h a d i s t r i b u t i o n is n o t o p t i m u m a t a n y given t i m e d u r i n g t h e c a l c u l a t i o n .

C . THE INTERNAL KINK MODE INSTABILITY T h e i n t e r n a l k i n k m o d e i n s t a b i l i t y in a diffuse screw p i n c h is d i s c u s s e d b y G o e d b l o e d a n d H a g e b e u k (1972). A r e c e n t p a p e r b y R o s e n b l u t h et al. (1973) includes b o t h linear a n d nonlinear results. I n t h e s e r e s u l t s , t h e g r o w t h o f a p e r t u r b a t i o n of a n u n s t a b l e e q u i l i b r i u m is c o m p u t e d . T h e p e r t u r b a t i o n is given b y ξ, = ξ0 e x p [ / ( / : z + m 0 ) ] .

(77)

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS 35 T h e p l a s m a e q u i l i b r i u m is defined b y t h e r a d i a l force b a l a n c e e q u a t i o n dP/dr

=

JZB(

(78)

w h e r e Jz is t h e a x i a l c u r r e n t d e n s i t y in t h e p l a s m a , a n d Be is t h e a z i m u t h a l field. T h e p l a s m a is infinitely c o n d u c t i n g , a n d fills t h e i n t e r i o r o f a c o n ­ d u c t i n g c y l i n d e r . T h e r a d i a l v a r i a t i o n o f t h e c u r r e n t d e n s i t y is w r i t t e n

r > a

(79)

w h e r e
F I G . 5. The surface of a cylindrical plasma in equilibrium is represented by an interior surface of the computation mesh. T h e lines exterior to the plasma surface outline the boundary of the mesh.

36

JEREMIAH U. BRACKBILL

(α)

(b)

F I G . 6. A cross section o f a cylindrical computation mesh is s h o w n for an internal kink m o d e calculation. In 6(a), the initial mesh is shown. In 6(b), the mesh is s h o w n at a time when an m = 1 perturbation has carried the plasma, which occupies the inner nine radial zones, almost t o the wall. 2

ρ = 1 0 " , t o t h e m a x i m u m d e n s i t y , ρ = 1. T h e v e l o c i t y , r e p r e s e n t e d b y v e c t o r s p o i n t i n g f r o m a m e s h p o i n t i n t h e d i r e c t i o n o f flow, a r e n o n z e r o i n i t i a l l y . A s B a t e m a n et al. ( 1 9 7 4 ) n o t e , it is e a s i e r t o g e n e r a t e c o n s i s t e n t i n i t i a l c o n ­ d i t i o n s w h e n a v e l o c i t y p e r t u r b a t i o n is u s e d . T h e p a t t e r n o f t h e v e l o c i t y v e c t o r s is t y p i c a l o f a n m = 1 p e r t u r b a t i o n f o r i n c o m p r e s s i b l e flow, w i t h t h e p l a s m a

(a)

(b)

F I G . 7. T h e isodensity contours corresponding t o the mesh shown in Fig. 6 are plotted. The contours labeled Η and L represent densities o f 1.0 and 0.01 (in arbitrary units), respectively. T h e values o f the density associated with intermediate contours form a g e o ­ metric progression.

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS 37

(α)

(b)

FIG. 8. The velocity vectors corresponding to the mesh in Fig. 6 are shown. Each vector being at a vertex and points in the direction of flow. The maximum velocity represented by a vector in 8(b) is roughly twice the maximum in 8(a). i n t e r i o r t o t h e s i n g u l a r surface ( w h e r e q = 1) m o v i n g b o d i l y t o w a r d t h e w a l l , a n d the p l a s m a exterior t o the singular surface being displaced. I n F i g s . 6 b a n d 8 b , t h e fully d e v e l o p e d i n s t a b i l i t y is s h o w n a t a t i m e c o r r e s p o n d i n g t o o n e Alfven t r a n s i t t i m e . T h e surface of t h e p l a s m a in F i g . 9, c o r r e s p o n d i n g t o F i g . 5, is c i r c u l a r in c r o s s s e c t i o n , b u t h a s suffered a helical

F I G . 9. The plasma surface is represented by an interior surface of the mesh, as in Fig. 5, at a time corresponding to that drawn in Figs. 6b, 7b, and 8b. The helical distortion of the plasma c o l u m n characteristic of an m = 1 perturbation is clearly shown.

38

JEREMIAH U. BRACKBILL

d i s p l a c e m e n t f r o m t h e g e o m e t r i c axis. T h e c o m p u t a t i o n m e s h in F i g . 6 b , c o r r e s p o n d i n g t o Fig. 6a, shows the distortion of the a l m o s t - L a g r a n g i a n m e s h d u e t o t h e p l a s m a flow. T h e c u r r e n t c a r r y i n g p l a s m a , w h i c h o c c u p i e d t h e i n n e r eleven r o w s o f z o n e s in t h e m e s h initially, h a s c a r r i e d z o n e s w i t h it a s it m o v e s t o w a r d t h e w a l l . T h e d e n s i t y g r a d i e n t a t t h e e d g e o f t h e c u r r e n t c a r r y i n g p l a s m a is s h o w n in F i g . 7 b . I t s size, a s i n d i c a t e d b y t h e s p a c i n g o f t h e c o n t o u r s , is n o t d i m i n i s h e d s u b s t a n t i a l l y f r o m its initial v a l u e , w h i c h c a n b e d e d u c e d f r o m F i g . 7a. F i n a l l y , t h e velocity v e c t o r s in F i g . 8 b s h o w s i m i l a r flow p a t t e r n s t o t h o s e s h o w n i n F i g . 8a a t t h e initial t i m e . T h e r e is n o t e n o u g h s p a c e h e r e for a d e t a i l e d a n a l y s i s of t h e s e r e s u l t s , b u t w e c a n c o n c l u d e f r o m t h e m t h a t t h e i n s t a b i l i t y a t l a r g e a m p l i t u d e is n o t q u a l i t a t i v e l y different f r o m t h e i n s t a b i l i t y a t s m a l l a m p l i t u d e . E v e n t u a l l y , flux t r a p p i n g b e t w e e n t h e w a l l a n d t h e p l a s m a will slow t h e m o t i o n o f t h e p l a s m a t o w a r d t h e w a l l , b u t , o t h e r t h a n flux t r a p p i n g , t h e r e s e e m s t o b e n o e v i d e n t s a t u r a t i o n m e c h a n i s m o p e r a t i n g . A m o r e careful a n a l y s i s of t h e in­ t e r n a l k i n k m o d e i n s t a b i l i t y is p l a n n e d w i t h t h e p a r t i c l e d i a g n o s t i c s d e s c r i b e d i n t h e last s e c t i o n .

V I . Conclusions T h e e m p h a s i s in t h i s article h a s b e e n o n a f u n d a m e n t a l p r o b l e m i n t h e s o l u t i o n of t h e e q u a t i o n s for i d e a l m a g n e t o h y d r o d y n a m i c flow. T h i s p r o b l e m , n a m e l y , t h e a c c u r a t e c o m p u t a t i o n o f c o n v e c t i v e t r a n s p o r t , e v i d e n t l y arises b e c a u s e o f n o n l i n e a r instabilities in n u m e r i c a l a p p r o x i m a t i o n s t o t h e c o n ­ vective d e r i v a t i v e . T h e r e is e v i d e n c e t h a t t h e s e instabilities a r e d u e t o c e r t a i n t r u n c a t i o n e r r o r s w h i c h p r o d u c e n e g a t i v e diffusion. T h e s e e r r o r s c a n b e c o m p u t e d , a n d t h e i r effect c a n b e s u p p r e s s e d e i t h e r b y a d d i n g c o m p e n s a t i n g p o s i t i v e diffusion o r b y r e d u c i n g t h e r e l a t i v e m o t i o n b e t w e e n t h e c o m p u t a t i o n mesh a n d the p l a s m a with a generalized mesh. A c o m p a r i s o n a m o n g the r e s u l t s o f l i n e a r stability t h e o r y , a g e n e r a l i z e d m e s h c a l c u l a t i o n , a n d a n u n ­ corrected Eulerian calculation supports the conclusion that more accurate c o m p u t a t i o n o f c o n v e c t i v e t r a n s p o r t l e a d s t o significant i n c r e a s e s in o v e r a l l accuracy. Similar improvements can be expected with corrected Eulerian calculations.

Appendix A c o m p a c t v e c t o r f o r m for t h e g e o m e t r i c coefficients, ca\ d i s c u s s e d in S e c t i o n I V , A , 1, h a s b e e n d e v e l o p e d b y D . C . B a r n e s ( u n p u b l i s h e d n o t e s , 1 1975). W h e r e τ is t h e p o s i t i o n v e c t o r for v e r t e x /, t h e g e o m e t r i c coefficient

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS 3 9 c o r r e s p o n d i n g t o v e r t e x 1 [ i n a cell w h o s e v e r t i c e s a r e l a b e l e d a s i n E q . ( 2 7 ) ] is w r i t t e n c

1

= ^-(r + r

2

5

χ r χ r

5

8

+ r + r

5

8

χ r χ r

4

4

+ r

+ r

4

4

χ r χ r

2

3

+ r

+ r

2

3

χ r

6

+ r

6

χ r

5

2

χ r ).

T h e coefficients f o r v e r t i c e s 2 - 4 a r e o b t a i n e d b y a cyclic p e r m u t a t i o n o f t h e i n d i c e s w i t h i n t h e t w o g r o u p s 1-4 a n d 5 - 8 . S i m i l a r l y , t h e g e o m e t r i c coefficient f o r v e r t e x 5 is w r i t t e n c

5

= ^-(r + r

6

1

χ r χ r

7

6

+ r + r

7

6

χ r χ r

8

8

+ r + r

8

8

χ r χ r

4

1

+ r + r

4

1

χ r

2

+ r

2

χ r

6

1

χ r ).

5

F r o m c , t h e coefficients f o r v e r t i c e s 6 - 8 a r e o b t a i n e d b y a cyclic p e r m u a t t i o n o f t h e i n d i c e s w i t h i n t h e t w o g r o u p s 1-4 a n d 5 - 8 .

ACKNOWLEDGMENTS

The author has received substantial help with this article and with the methods described in it from many people. A m o n g them are W . E . Pracht, w h o s e generously provided three-dimensional hydrodynamics program simplified enormously the task o f writing a c o d e for magnetohydrodynamics; D . C. Barnes, w h o checked and m a d e many improvements and extensions t o the equations described in Section I V ; J. P. Freidberg, w h o suggested many o f the problems t o which the generalized mesh m e t h o d has been applied; and C. W . N i e l s o n , w h o s e critical reading o f the manuscript is responsible for many improvements. T h e author also thanks W . B . G o a d and C. W . Hirt for many helpful conversations and suggestions. This work was supported by U S E R D A , Contract N o . W - 7 4 0 5 - E N G . 36.

REFERENCES

A m s d e n , Α . Α . , and Harlow, F . H . (1970). " T h e S M A C M e t h o d : A Numerical Technique for Calculating Incompressible Fluid F l o w s , " R e p . N o . L A - 4 3 7 0 , L o s A l a m o s Sei. L a b . , Los A l a m o s , N e w M e x i c o . A m s d e n Α . Α . , and Hirt, C. W . (1973). " Y A Q U I : A n Arbitrary Lagrangian-Eulerian Computer Program for Fluid F l o w at A l l Speeds," R e p . N o . L A - 5 1 0 0 , L o s A l a m o s Sei. Lab., L o s A l a m o s , N e w M e x i c o . Anderson, D . V. (1975). / . Comput. Phys. 17, 2 4 6 . Bateman, G . , Schneider, W . , and G r o s s m a n n , W . (1974). Nucl. Fusion 14, 669. B o d i n , Ν . A . B . (1972). Nucl. Fusion 1 2 , 721. Boris, J. P. (1970). " A Physically Motivated Solution o f the Alfven P r o b l e m , " M e m o . R e p . N o . 2167, N a v a l R e s . Lab., W a s h i n g t o n , D . C . Boris, J. P., and B o o k , D . L. (1973). J. Comput. Phys. 1 1 , 38.

40

JEREMIAH U. BRACKBILL

Brackbill, J. U . (1973). Proc. Conf Numer. Simul. Plasmas, 6th, 1973 Lawrence Livermore Lab. Conf. R e p . 730804. Brackbill, J. U . , a n d Pracht, W . E . (1973). / . Comput. Phys. 13, 455. Butler, T. D . , Henins, L, Jahoda, F . C , Marshall, J., and Morse, R. L. (1969). Phys. Fluids 12,1904. D u c h s , D . (1968). Phys. Fluids 1 1 , 2 0 1 0 . Fornberg, Β. (1973). Math Comput. 27, 4 5 . Freeman, J. R. (1971). Nucl. Fusion 1 1 , 4 2 5 . Freeman, J. R., a n d Lane, F . O. (1968). Proc. APS Top. Conf. Numer. Simul. Plasmas, 1968 L o s A l a m o s Sei. Lab. R e p . N o . L A - 3 9 9 0 . F r o m m , J. E. (1961). "Lagrangian Difference Approximations for Fluid D y n a m i c s , " R e p . N o . L A - 2 5 3 5 , L o s A l a m o s Sei. L a b . , L o s A l a m o s , N e w Mexico. G o a d , W . Β. (1960). W A T : A Numerical M e t h o d for T w o - D i m e n s i o n a l Unsteady Fluid F l o w , " R e p . L A M S - 2 3 6 5 , L o s A l a m o s Sei. L a b . , L o s A l a m o s , N e w M e x i c o . G o e d b l o e d , J. P., and Hagebeuk, H . J. L. (1972). Phys. Fluids 1 5 , 1 0 9 0 . H a i n , Κ., Hain, G., Roberts, Κ. V., Roberts, S. J., and Koppendorfer, W. (1960). Z. Naturforsch. A 15, 1039. Harlow, F . H . , a n d A m s d e n , A . A . (1971). / . Comput. Phys. 8, 197. Harlow, F . H . , and Welch, J. E. (1965). Phys. Fluids 8, 842. Hertweck, F . , and Schneider, W . (1970). " A T w o - D i m e n s i o n a l Computer Programme for Solving the M H D Equations for the Theta Pinch in a Time D e p e n d e n t Coordinate System, R e p . N o . I P P 1/110. Inst. Plasma Phys., Garching. Hirt, C. W. (1968). / . Comput. Phys. 2 , 339. Hirt, C. W . , C o o k , J. L., and Butler, T. D . (1970). / . Comput. Phys. 5, 103. Hirt, C. W . , A m s d e n , Α . Α . , and C o o k , J. L. (1974). / . Comput. Phys. 14, 2 2 7 . H o f m a n n , J. (1974). Nucl. Fusion 14, 438. Killeen, J. (1972). In "Information Processing 7 1 " (C. V. Freiman, J. E. Griffith, and J. L. Rosenfeld, eds.), p . 1191. N o r t h - H o l l a n d Publ., Amsterdam. Kreiss, H. O. (1964). Comm. Pure Appl Math. 17, 335. Landau, L. D . , and Lifschitz, Ε. M . (1959). "Fluid Mechanics." Addison-Wesley (Pergamon), Reading, Massachusetts. Lapidus, A . (1967). / . Comput. Phys. 2 , 154. Lewis, H . R. (1966). J. Appl. Phys. 37, 2 5 4 1 . Lindemuth, I., and Killeen, J. (1973). / . Comput. Phys. 13, 181. Lui, H . C. (1973). " F L I C Codes o f Shock Focusing in a Coaxial Electromagnetic Shock T u b e , " Lab. R e p . N o . 60, Columbia University, N e w York. Lui, H . C , and C h u , C. K. (1974). Int. Conf. Numer. Methods Fluid Dyn. 4th, 1974 p. 263. Morse, P. M . , and Feshback, H . (1953). " M e t h o d s of Theoretical Physics," Chapter 3. M c G r a w - H i l l , N e w York. Peaceman, D . W . , a n d Rachford, Η . Η . (1955). / . Soc. Ind. Appl. Math. 3 , 2 8 . Pracht, W . Ε. (1975). / . Comput. Phys. 17, 132. Richtmyer, R. D . (1963). " A Survey of Difference M e t h o d s for Non-steady Fluid D y n a m i c s , " N C A R Tech. N o t e 6 3 - 2 . N a t . Cent. A t m o s . R e s . , Boulder, Colorado. Richtmyer, R. D . , and M o r t o n , K. W . (1967). "Difference M e t h o d s for Initial Value Problems." Wiley (Interscience), N e w York. Rivard, W. C , Farmer, Ο. Α . , Butler, J. D . , and O'Rourke, P. J. (1973). " A M e t h o d for Increased Accuracy in Eulerian Fluid D y n a m i c s Calculations," R e p . N o . L A - 5 4 2 6 - M S , Los A l a m o s Sei. L a b . , L o s A l a m o s , N e w M e x i c o . Roberts, Κ. V., and Boris, J. P. (1969). 3rd Annu. Numer. Plasma Simul. Conf., 1969 Paper n o . 32.

NUMERICAL MAGNETOHYDRODYNAMICS FOR HIGH-BETA PLASMAS 4 1 Roberts. Κ. V., and Potter, D . E . (1968). Methods Compt. Phys. 9, 339. Rosenbluth, Μ . N . , D a g a z i a n , R. Y . , a n d Rutherford, P. H . (1973). Phys. Fluids 16, 1894. Schneider, W . (1972). Z. Phys. 2 5 2 , 147. Schulz, W . D . (1964). Methods Comput. Phys. 3 , 1. T h o m p s o n , J. F . , T h a m e s , F . C , a n d Mastin, C . W . (1974). / . Comput. Phys. 15, 299. Tuck, J. (1968). Eur. Conf. Controlled Fusion Plasma Phys. Res., 2nd, 1968 V o l . 2, p. 595. Wagner, C. E . , a n d Manheimer, W . M . (1973). Proc. Conf. Numer. Simul. Plasmas, 6th 1973 Lawrence Livermore L a b . Conf. R e p . 730804. White, R., M o n t i c e l l o , D . , R o s e n b l u t h , Μ . N . , and Strauss, N . (1974). Plasma Phys. Fusion Res., 1974, V o l . 1 , p. 495. W o o t e n , J., Hicks, H . R., Bateman, G., a n d D o r y , R. A . (1974). "Preliminary Results o f the 3 D N o n l i n e a r Ideal M H D C o d e , " O R N L T M 4784. Hollifield N a t . L a b . , Oak Ridge, Tennessee. Y o u n g , D . M . (1962). In " A Survey o f Numerical A n a l y s i s " (J. T o d d , ed.), Chapter 11. M c G r a w - H i l l , N e w York. ZePdovich, Y a . Β . , a n d Raizer, Y u . P. (1967). "Physics of Shock Waves and H i g h Tempera­ ture H y d r o d y n a m i c P h e n o m e n a , " 2nd e d . , V o l . 2 . A c a d e m i c Press, N e w York. Zienkiewicz, O. C. (1971). " T h e Finite-Element M e t h o d in Engineering Science." M c G r a w Hill, N e w Y o r k .