8 Equations of guiding-center motion—outdated, will be deleted

YJ’s remark: This section is outdated because I have found a more accurate form of the equations of the guiding center motion, as given in Sec. 1. Besides to be more accurate, the new form is compact and suitable for numerical implementation

The equations for the guiding center motion are given by (refer to my note “collisionless_drift_kinetic_equation.tm”)

               (           )
X˙ = v∥b + 1-b×   μ-∇B + v2∥κ ,
          Ω      m
(147)

and

       μ
v˙∥ = −m-b ⋅∇B + v∥κ ⋅ ˙X,
(148)

where μ mv2(2B), b = B∕B. The term vκ˙
X in Eq. (148) can be further simplified by using Eq. (147), which gives

       μ         v∥μ
v˙∥ = −m-b ⋅∇B +  mΩ-κ⋅b × ∇B
(149)

where use has been made of that the curvature κ is perpendicular to b.

[Benchmark: In GTC simulation (refer to H. Zhang’s paper[6]), the time derivative of v is given by

      B + Bv∥∕Ωα∇ × b
˙v∥ = −------mB--------⋅μ∇B
(150)

          μ-        v∥μ-
=⇒  ˙v∥ = − m b⋅∇B − mΩ ∇ × b ⋅∇B,
(151)

Next, I prove Eq. (149) is equivalent to Eq. (151). Using κ b⋅∇b = b×∇×b, the second term on the right-hand side of Eq. (149) is written as

v∥μ-                v∥μ-
m Ω κ⋅b × ∇B   =  − m Ω(b ×∇ × b )⋅(b × ∇B )
               =  − v∥μ∇ × b ⋅∇B + v∥μ-(b ⋅∇B )(∇ ×b ⋅b)
                    m Ω            m Ω
               ≈  − v∥μ∇ × b ⋅∇B + 0,                          (152)
                    m Ω
where in obtaining the last equality, use has been made of that ∇× b b 0 (note that this is correct to the order considered here, I will discuss this later). Using  Eq. (152) in Eq. (149) yields
      μ         v∥μ
˙v∥ = − m-b⋅∇B − m-Ω∇ × b ⋅∇B,
(153)

which is identical with Eq. (151).]

 

8.1 Equilibrium magnetic field in tokamak

The tokamak equilibrium magnetic field can be written

B = ∇Ψ × ∇ ϕ+ g(Ψ)∇ ϕ.
(154)

In my code, the values of the two free functions, Ψ = Ψ(R,Z) and g(Ψ), which specify the magnetic field, is read from the output file “G-eqdsk-file” of EFIT code. Using Eq. (154), the axisymmetric equilibrium magnetic field can be written as

BR = −-1∂Ψ-,
      R ∂Z
(155)

BZ =  1-∂Ψ,
      R ∂R
(156)

     g(Ψ)
Bϕ =  R  .
(157)

The partial derivative of the component of the magnetic field is written as

∂BR-  -1-∂Ψ-   1-∂2Ψ--
∂R  = R2 ∂Z −  R∂Z ∂R
(158)

∂BR-    1-∂2Ψ-
∂Z  = − R ∂Z2
(159)

∂BZ-    -1-∂Ψ-  1-∂2Ψ-
∂R  = − R2 ∂R + R ∂R2
(160)

∂BZ-  -1-∂2Ψ--
∂Z  = R ∂R ∂Z
(161)

∂Bϕ-    1--      1-′   ∂Ψ-
∂R  = − R2g(Ψ) + Rg (Ψ)∂R
(162)

∂Bϕ-= -1g′(Ψ )∂Ψ-
∂Z    R      ∂Z
(163)

         ∘ --------------------
       1   (∂ Ψ)2   (∂ Ψ)2
⇒  B = --   ---   +  ---   + g2
       R    ∂R        ∂Z
(164)

∂B         1      1 1  1 (  ∂Ψ ∂2Ψ    ∂Ψ  ∂2Ψ      ∂g )
∂R-  =  − R2-BR + R-2BR-- 2 ∂R-∂R2-+ 2∂Z-∂Z∂R--+ 2g∂R-
                   (     2         2        )
     =  − B-+ --12  ∂-Ψ∂-Ψ2 + ∂Ψ--∂-Ψ--+g ∂g-                  (165)
          R   BR    ∂R ∂R     ∂Z ∂Z∂R     ∂R
               (                            )
∂B-  =  -11 -1-- 2∂Ψ--∂2Ψ--+ 2∂Ψ-∂2Ψ-+ 2g ∂g
∂Z      R 2 BR    ∂R ∂R∂Z     ∂Z ∂Z2     ∂Z
          1  (∂ Ψ ∂2Ψ    ∂Ψ ∂2Ψ     ∂g)
     =  BR2-  ∂R-∂R-∂Z-+ ∂Z-∂Z2- +g ∂Z-                   (166)

In my numerical code, the numerical data of the poloidal flux function Ψ(R,Z) and toroidal field function g(Ψ) are read in from the output G-eqdsk file of the EFIT code. Then all the partial derivatives are calculated by using central finite-difference. The linear interpolation is used to evaluate the various quantities that are needed at the instantaneous location of guiding-centers to push the orbits.

 

8.1.1 Solovév equilibrium

When I began to write the guiding center orbit code, in order to avoid the numerical interpolating, I use Solovev’s analytic equilibrium. (The latest version of my code constructs magnetic field by reading the output G-eqdsk file of the EFIT code and thus can treat general tokamak magnetic field.) The Solovév equilibrium is an analytic equilibrium in which the poloidal flux function Ψ is given by

       B    [       κ2         ]
Ψ = ---20--  R2Z2 + -0(R2 − R20)2 ,
    2R 0κ0q0         4
(167)

where B0, R0, κ0, q0 are constant parameters. Using Eq. (167), the partial derivatives are written as

∂Ψ      B                         ∂Ψ     B
---= ---20-- [2RZ2 + κ20(R2 − R20)R],---= --2-0-R2Z.
∂R   2R 0κ0q0                     ∂Z   R 0κ0q0
(168)

∂2Ψ      B
---2 =---20-- [2Z2 + κ20(3R2 − R20)]
∂R    2R 0κ0q0
(169)

∂2Ψ     B        ∂2Ψ      2B
---2 = -2-0--R2,------= --2-0-RZ
∂Z     R0κ0q0   ∂R ∂Z   R 0κ0q0
(170)

(171)

The toroidal field function g is a constant function, g = cgR0B0, where cg is a dimensionless constant.