PvdZon 2012 - les09

MagnetohydroSTATICS

Contents

1.

Equations for static equilibria

2.

Cylindrical symmetry

the diffuse linear plasma column (1D tokamak, loops, etc.)

3.

The Grad-Shafranov equation

two-dimensional (cartesian coordinates), derivation

4.

2D Tokamak equilibria

solution strategy, Soloviev solution,. . .

5.

Astrophysical equilibria

solution strategy, coronal loop models,. . .

Motivation

•

MHD equilibrium is an absolute requirement for thermonuclear fusion in a tokamak

⇒

MHD spectroscopy: small perturbations of static equilibrium

•

solar coronal structures evolve on time scales much larger than the dynamic time scale

(waves)

⇒

sunspot seismology: small perturbations of static equilibrium

⇒

AR seismology: small perturbations of static equilibrium

⇒

coronal loop seismology: small perturbations of static equilibrium

Equations for static equilibria

⇒

MHD equations with:

∂

∂t

= 0

and

v = 0

∂ρ

∂t

+ ∇ · (ρv) = 0

,

ρ

Dv

Dt

+

∇p − j × B − ρg = 0

(force balance)

,

∂p

∂t

+ v · ∇p +

5 3

p

∇ · v = 0

,

∂B

∂t

+ ∇ × E = 0

,

j = ∇ × B

(Amp `ere’s law)

,

E + v × B

= 0 ,

whereas

⇒

the dimensionless MHS equations read:

∇p = j × B + ρg

(force balance)

j = ∇ × B

(Amp `ere’s law)

∇ · B = 0

(solenoidal condition)

+ BCs: e.g.

n · B = 0

on an imposed ‘wall’ (tokamak)

Cylindrical symmetry

•

diffuse cylindrical plasma column = one of the most useful models for study of confined plasmas

⇒

probably the most widely studied model in plasma stability theory

•

used as 1st approximation to toroidal systems (tokamaks, loops, flux ropes,. . . ) (2nd dir. of inhomogeneity: much more complicated analysis (PDEs))

Diffuse cylindrical plasma column with a helical magnetic fieldB, surface currents atr = a, and surrounded by a vacuum magnetic field ˆB.

• a r z θ B μ−1

•

consider diffuse plasma in cylinder of radius

a

•

equilibrium equations (neglecting gravity (often OK!)):

j × B = ∇p ,

j = ∇ × B ,

∇ · B = 0

•

in cylindrical coordinates

r, θ, z

(only

r

-dependence, _∂θ∂

=

_∂z∂

= 0

, and 

≡

_drd ):

∇ · B = 0 ⇒

∂rB

r

∂r

= 0 → rB

r

=

const (0 inr=0)

−→

B

_r

= 0,

p



= j

_θ

B

_z

− j

_z

B

_θ

,

j

_θ

= −B

_z

,

j

_z

=

1

r

(rB

θ

)



•

eliminating

j

θ and

j

z :

[ p(r) +

1₂

B

2

(r)]



+

B

2 θ

(r)

r

= 0

(a single ODE!)

⇒

we may choose two of the three profiles

p

(r), B

_θ

(r), B

_z

(r)

arbitrarily

Equilibrium of a z-pinch

•

current in

z

-direction, magnetic field in

θ

-direction, gradients in

r

-direction

⇒

equations reduce to:

dp

dr

= −j

z

B

θ

and j

z

=

1

r

d

dr

(rB

θ

)

⇒

dp

dr

= −

B

θ

r

d

dr

(rB

θ

)

(this is the only restriction on these profiles)

•

choose, e.g., a constant current profile

j

_z = const

⇒ B

θ

=

1₂

rj

z

,

and p

= p

c

(1 − r

2

/a

2

) , p

c

≡

1₄

a

2

j

_z2 a b 0 r j_z a b r 0 B_θ p r a b 0

•

surrounding vacuum:

ˆB

satisfies:

∇ × ˆB = 0

,

∇ · ˆB = 0

⇒

ˆj

z

= 0

and

B

ˆ

θ

(r) = B

θ

(a) a/r

(this radially decaying magnetic field is produced by the total current

I

z flowing within

the plasma interval

0 ≤ r ≤ a

)

•

the configuration may be closed off by putting a perfectly conducting wall at some radius

r

= b

a b 0 r j_z a b r 0 B_θ p r a b 0

Equilibrium of a

θ

-pinch

•

current induced in

θ

-direction, magnetic field in

z

-direction, gradients in

r

-direction

⇒

equations reduce to:

dp

dr

= j

θ

B

z

and j

θ

= −

dB

z

dr

⇒

dp

dr

= −

d

dr



B

_z2

2



⇒ p +

B

z2

2

=

const

–

extreme case: skin-current

θ

-pinch

⇒

fast current induction

⇒

‘magnetic piston’

⇒ p

1

+



B

_z2

2



1

=



B

_z2

2



2 (

p

= nkT

⇒ nT

high, but

τ

short)

j B₁ B 2 j p μ₀ 2 B2

Tokamak:

delicate balance between equilibrium & stability

z - pinch:

very unstable (remains so in a torus)

θ - pinch:

end-losses

(in torus: no equilibrium)

B j

B j

‘Straight tokamak’ limit

•

identifying the ends of a cylinder of finite length

L

= 2πR

₀

⇒

1st approximation of a slender torus with small inverse aspect ratio



≡ a/R

₀

 1

Slender torus with inverse aspect ratio  ≡ a/R₀  1represented as a periodic cylinder with length2πR₀.

R_o a r z R_o 2π a ϕ

q

(r) ≡

rB

z

(r)

R

₀

B

θ

(r)

≡

1

μ

(r)R

₀

inverse pitchμand safety factorqof the mag-netic field lines in a ‘straight tokamak’ periodic

cylinder model of a toroidal plasma.

0

B

2

π

r

2

π

R

0

2

π

μ

=

q

.

2

π

R

0

⇒ q

is ’safety factor’ because

q <

1 ⇒

MHD instabilities

• q(0) > 1

required for stability internal kink modes

• q(1) > 1

required for stability external kink modes

Example

: class of ’1D tokamak’-equilibria (model II)

•

choose:

B

z

= 1

j

z

=

j

0

(1 − r

2

)

ν

ρ

= 1 − (1 −

d

)r

2 with

0 ≤ r ≤ 1

• j

z

=

1

r

d

dr

(rB

θ

) ⇒ rB

θ

= j

0



(1 − r

2

)

ν

r dr

⇒ rB

θ

= −

j

₀

2



(1 − r

2

)

ν

d

(1 − r

2

) = −

j

0

2(ν + 1)

(1 − r

2

)

ν+1

_{+ C}

choose

C

so that

rB

θ

(0) = 0

(

B

θ

(0)

finite)

⇒ B

θ

=

1

2r

j

₀

ν

+ 1

[1 − (1 − r

2

₎

ν+1

_]

• j

θ

= −

dB

z

dr

= 0 ⇒

dp

dr

= −j

z

B

θ

= −

1

2r

j

₀2

ν

+ 1

[1 − (1 − r

2

₎

ν+1

_{](1 − r}

2

₎

ν

⇒ p = −

1

2

j

₀2

ν

+ 1



[1 − (1 − r

2

₎

ν+1

](1 − r

2

)

ν

r

2

r

dr

⇒

e.g.

ν

= 1

:

p

=

j02 4

(

(1−r 2₎3 6

+

(1−r 2₎2 4

)

vacuum profiles

:

• ˆj

z

= 0,

ˆp = 0

⇒

choose

B

ˆ

z

= 1

⇒ ˆ

B

θ

=

_2r1 _νj₊₁0

⇒

no surface current

•

a r B B^ z θ

Remark

:

P

= p +

B₂2

≈

Bz2 2 when

β

=

B2p2

 1

– e.g. tokamak:

B

= 3

Tesla

⇒ P = 3.6 × 10

6

N/m

2

= 36

atm

β

= 3% ⇒ p ≈ 1

atm

⇒

pressure effects small, not well described in cylindrical geometry

⇒

2D toroidal geometry required!

The Grad-Shafranov equation





derivation in cartesian coordinates

•

force balance:

∇p = j × B + ρg

with

g = −g e

z and

ρ

= p/T

(dim.less)

• ∇ · B = 0 ⇒ ∃ψ(x, z) : B = ∇ψ(x, z) × e

_

_

_

y poloidal

+ B

_{  }

y

e

y toroidal

⇒ B

x

= −

∂ψ_∂z

,

B

z

=

∂ψ_∂x

∇ · j = 0 ⇒ ∃I(x, z) : j = ∇I(x, z) × e

y

+ j

y

e

y

⇒

j

_x

= −

∂I_∂z

,

j

_z

=

_∂x∂I

•

Amp `ere:

j = ∇ × B ⇒

j

x

= −

∂By ∂z

,

j

z

=

∂By ∂x

⇒

I

≡ B

y

⇒ j

y

=

∂B_∂zx

−

∂B_∂xz

= −[

∂

2

ψ

∂z

2

+

∂

2

ψ

∂x

2







≡Δ∗_ψ

]

•

force balance:

•

force balance:

∂

p

∂y

= j

z

B

x

− j

x

B

z

= −

∂I

∂x

∂ψ

∂z

+

∂I

∂z

∂ψ

∂x

= 0

⇔ I = I(ψ)

∂p

∂x

= j

y

B

z

−

j

z

B

y

= j

y

∂ψ

∂x

−

∂I

∂ψ

∂ψ

∂x

I

= (j

y

− I



_I

₎

∂ψ

∂x

•

force balance:

∂y

= j

z

B

x

− j

x

B

z

= −

∂x

∂z

+

∂z

∂x

= 0

⇔ I = I(ψ)

∂p

∂x

= j

y

B

z

−

j

z

B

y

= j

y

∂ψ

∂x

−

∂I

∂ψ

∂ψ

∂x

I

= (j

y

− I



_I

₎

∂ψ

∂x

∂p

∂z

= j

x

B

y

− j

y

B

x

−

p

T

g

= (j

y

− I



_I

₎

∂ψ

∂z

−

p

T

g

•

force balance:

∂

p

∂y

= j

z

B

x

− j

x

B

z

= −

∂I

∂x

∂ψ

∂z

+

∂I

∂z

∂ψ

∂x

= 0

⇔ I = I(ψ)

∂p

∂x

= j

y

B

z

−

j

z

B

y

= j

y

∂ψ

∂x

−

∂I

∂ψ

∂ψ

∂x

I

= (j

y

− I



_I

₎

∂ψ

∂x

∂p

∂z

= j

x

B

y

− j

y

B

x

−

p

T

g

= (j

y

− I



_I

₎

∂ψ

∂z

−

p

T

g

•

now use

ψ

as a coordinate!

⇒ p(x, z) ≡ p(ψ, z)

and combine: (

x

-comp) ∂ψ_∂z - (

z

-comp) ∂ψ_∂x

⇒

∂p

∂x

∂ψ

∂z

−

∂p

∂z

∂ψ

∂x

=

p

T

g

∂ψ

∂x

•

force balance:

∂y

= j

z

B

x

− j

x

B

z

= −

∂x

∂z

+

∂z

∂x

= 0

⇔ I = I(ψ)

∂p

∂x

= j

y

B

z

−

j

z

B

y

= j

y

∂ψ

∂x

−

∂I

∂ψ

∂ψ

∂x

I

= (j

y

− I



_I

₎

∂ψ

∂x

∂p

∂z

= j

x

B

y

− j

y

B

x

−

p

T

g

= (j

y

− I



_I

₎

∂ψ

∂z

−

p

T

g

•

now use

ψ

as a coordinate!

⇒

p

(x, z) ≡ p(ψ, z)

and combine: (

x

-comp) ∂ψ_∂z - (

z

-comp) ∂ψ_∂x

⇒

∂p

∂x

∂ψ

∂z

−

∂p

∂z

∂ψ

∂x

=

p

T

g

∂ψ

∂x

⇒

∂p

∂ψ

∂ψ

∂x

∂ψ

∂z

−

(

∂p

∂ψ

∂ψ

∂z

+

∂p

∂z

)

∂ψ

∂x

=

p

T

g

∂ψ

∂x

⇒ (

∂p

∂z

+

p

T

g

)

∂ψ

∂x

= 0

∂ψ ∂x =0

=⇒

∂p

∂z

+

p

T

g

= 0

⇒

∂

ln p

∂z

= −

g

T

⇒ p = p

0

(ψ) e

− ₀z _Tg dz

•

substituting this in the force balance makes both the

x

- and

z

-components reduce to:

∂p

∂ψ

∂ψ

∂z

+

∂p

∂z



=−_Tp g

= −(Δ

∗

ψ

+

1

2

∂I

2

∂ψ

)

∂ψ

∂z

−

p

T

g

⇒

Δ

∗

ψ

+

∂

∂ψ

I

2

2

+ p

0

(ψ) e

− z 0 Tg dz



= 0

TOKAMAK equilibria

  

strategy

(here

g

= 0 ⇒ p = p

₀

(ψ)

)

•

choose

I

(ψ)

and

p

₀

(ψ)

(on the basis of data fitting)

•

solve the Grad-Shafranov equation

⇒ ψ = ψ(x, z)

(solution consists of closed magnetic surfaces, centered around a magnetic axis) (this is actually a difficult numerical problem!!)

•

construct a flux coordinate system

(ψ, θ, ϕ)

ϕ

: as in cylindrical coordinates

ψ

: labels flux surfaces

θ

: can be constructed so that magnetic field lines are straight!

•

all geometrical and physical quantities appearing in the operators of the MHD equations can than be expressed in terms of the

(ψ, θ, ϕ)

coordinates

The Grad-Shafranov equation

  

in cylindrical coordinates...

•

in cylindrical coordinates, the 2D equilibrium is described by an elliptic non-linear PDE,

the Grad-Shafranov equation:

[ Δ

∗

Ψ ≡ ] R

∂

∂R

1

R

∂

Ψ

∂R

+

∂

2

Ψ

∂Z

2

= −II



_{− R}

2

_p



_{[ = Rj}

ϕ

]

which has to satisfy the boundary condition

Ψ = const

(on the plasma cross-section)

(the special symbol

Δ

∗ indicates a Laplacian-like operator with the positions of the factors

R

and

1/R

reversed as compared to the usual Laplacian)

2D Tokamak equilibria

i

R

R

₀

a

φ

R

R

_e

Z

Schematic representation of the tokamak geometry with cylindrical coordinates (R, φ, Z). R_e and R_i are the two

values of the major radius where the last closed magnetic surface crosses the equatorial plane, where R_i < R_e; a= (R_e− R_i)/2 and R₀ = (R_e+ R_i)/2.

Soloviev equilibrium

The so-called Soloviev equilibrium is an analytic solution of the Grad-Shafranov equation where the two equilbrium profiles

P

(ψ)

and

II



(ψ)

are independent of

ψ

.

In terms of the ellipticity

E

= b/a

and the triangularity

τ

of the plasma boundary, the general up-down symmetric solution is given by:

ψ

= (x −



2

(1 − x

2

))

2

+ ((1 −



2

4

)(1 + x)

2

+ τx(1 +

x

2

))

y

2

E

2

,

where

x

and

y

are the normalized coordinates:

x

= (R − R

₀

)/a

and

y

= Z/a

used in the poloidal plane and



= a/R

₀

- 1 . 0 0 - . 5 0 . 0 0 . 5 0 1 . 0 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 0 1 . 0 0 X PS I F LUX - 1 . 0 0 - . 5 0 . 0 0 . 5 0 1 . 0 0 . 0 0 . 2 0 . 4 0 . 6 0 . 8 0 1 . 0 0 X P/P 0 PRESSURE - 1 . 0 0 - . 5 0 . 0 0 . 5 0 1 . 0 0 . 0 0 1 . 0 0 2 . 0 0 3 . 0 0 4 . 0 0 5 . 0 0 X J CURRENT DENS I TY - 1 . 0 0 - . 5 0 . 0 0 . 5 0 1 . 0 0 1 . 0 0 2 . 0 0 3 . 0 0 4 . 0 0 5 . 0 0 X q q - PROF I L E

A

D

A

D

The ψ, P(ψ), j_φ(ψ), and q(ψ)-profiles for the equilibria that yield maximum β for four different pressure profiles. From

Astrophysical equilibria

•

same solution method would be possible but would lead to problem: there is no info about the profiles

p

₀

(ψ)

and

I

(ψ)

!

•

moreover:

p

is not a flux function anymore (gravity!)

•

alternative solution strategy used in literature:

–

fix the geometry of the magnetic field, i.e. choose

ψ

(x, z)

and

I

(ψ)

!

–

determine

p

(ψ, z)

, i.e. the pressure required to keep this magnetic field in equilib-rium!

⇒

the Grad-Shafranov equation then reduces to a first-order ODE for

p

₀

(ψ)

!

Δ

∗

ψ

+

∂

∂ψ

⎡

⎢

⎣

I

₂

2

+ p

0

(ψ)e

−



z 0

g

T

dz



⎤

⎥

⎦ = 0

•

writing

Δ

∗

ψ

as

f

(ψ, z)

and integrating the Grad-Shafranov (in

ψ

) then yields:



ψ 0

Δ

∗

_ψdψ



₊

I

2

2

+ p

0

(ψ)e

−



z 0

g

T

dz



= Q(z)

where

Q

(z)

is an arbitrary function of only

z

(as a result of the integration in

ψ

)

•

there are many such solutions published for loops, arcades, etc. e.g. Low (’79):

ψ

=

2(1 + z)

1 + x

2

_{+ (1 + z}

2

₎

and

I

≡ 0

⇒ Δ

∗

_ψ

₌

−16(1+z) [1+x2_+(1+z2_)]3

=

_[1+x2_+(1+z−8 2_)]2

2(1 + z)

1 + x

2

_{+ (1 + z}

2

₎







=ψ

=

_(1+z)−2 2

ψ

3

⇒



ψ 0

Δ

∗

_ψdψ



₌



ψ 0

−2

(1 + z)

2

ψ

3

dψ



=

−2

(1 + z)

2

ψ

4

4

⇒

choose:

Q

(z) =

Q

0

(1 + z)

2 with

Q

0 a constant

⇒

the integrated GS equation then reduces to

Q

₀

(1 + z)

2

+

ψ

4

2(1 + z)

2

= p

0

(ψ)e

−



z 0

g

T

dz

 yielding

p

= p

₀



8(1 + z)

2

[1 + x

2

_{+ (1 + z}

2

_)]

4

+

Q

₀

(1 + z)

2



ρ

= ρ

₀



16(1 + z)

[1 + x

2

_{+ (1 + z}

2

_)]

4

+

2Q

0

(1 + z)

3



T

= T

₀



1 + z

2



• Task

: compare this solution with that for

ψ

=

√

2(

√

2 + z)

x

2

+ (

√

2 + z)

2

,

I

≡ 0