Computation of two-fluid flowing equilibrium of spherical torus plasma using multi-grid method

Abstract

Two-fluid (ions and electrons) flowing equilibrium configurations of spherical torus plasma are numerically determined by using the multi-grid method. The axisymmetric equilibrium equations consist of a pair of second-order partial differential equations for the magnetic (electron) and ion stream functions, and Bernoulli equations for the density. It is found from the numerical results that the multi-grid method with the damped Jacobi method in the smoothing step is effective for solving these equations with numerical stability, high accuracy, and high speed.

Keywords

Two-fluid plasma flowing equilibrium spherical torus multi-grid method numerical analysis

1. Introduction

Spherical torus (ST) [1] which is the tokamak with the very small aspect ratio less than 2 is considered as candidates for magnetic confinement fusion reactors because of its high confinement, stability, and beta (ratio of plasma pressure to magnetic pressure). Recently, it is recognized that the effect of plasma flow is important in understanding the confinement of high beta torus plasmas, and that two-fluid effects which consists of ion diamagnetic and inertial effects are significant for smaller scale length, higher beta, and flow velocity approaching the ion diamagnetic drift velocity [2]. The formalism for two-fluid flowing equilibria is developed [3]. This system is described by a pair of second-order partial differential equations for two surface variables, i.e., the magnetic and ion stream functions and Bernoulli equations for the density. In order to solve this system, three arbitrary functions for each species are required to be assumed as functions of the surface variables stated above. The purpose of this study is to solve these equations by the multi-grid method [4] and to investigate the properties of numerical convergence. We focus our attention on the equilibria for geometry and boundary conditions of the rapidly rotating, high-performance National Spherical Torus eXperiment (NSTX) device [1].

Fig. 1.

Computational domain and boundary condition in cylindrical coordinate (x, θ, z), corresponding to the confinement region in the NSTX device.

2. Numerical model

Let us assume axisymmetry about NSTX geometric axis in cylindrical coordinates (r, θ, z), as shown in Fig. 1. An axisymmetric two-fluid flowing equilibrium [5] is described by a pair of second-order partial differential equations for ion surface variable Y (r, z) and electron surface variable 𝜓(r, z) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle \bar{{\psi}}_{i}^{\prime }r^{2}\boldsymbol{\nabla}\cdot \left(\frac{\bar{{\psi}}_{i}^{\prime }}{n}\displaystyle \frac{\boldsymbol{\nabla}Y}{r^{2}}\right)=\displaystyle \frac{r}{{\varepsilon}}(B_{{\theta}}\bar{{\psi}}_{i}^{\prime }-nu_{{\theta}})+nr^{2}(F_{i}^{\prime }-T_{i}^{\prime }\ln n),\end{eqnarray}$ (1) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle r^{2}\boldsymbol{\nabla}\cdot \left(\displaystyle \frac{\boldsymbol{\nabla}{\psi}}{r^{2}}\right)=\frac{r}{{\varepsilon}}(B_{{\theta}}\bar{{\psi}}_{e}^{\prime }-nu_{{\theta}})-nr^{2}(F_{e}^{\prime }-T_{e}^{\prime }\ln n),\end{eqnarray}$ (2) and Bernoulli equations lead to an explicit equation for the density n (r, z) $\begin{eqnarray}\displaystyle n=\exp \left[\displaystyle \frac{F_{i}+F_{e}+u^{2}/2}{T_{i}+T_{e}}-1\right]. & & \displaystyle\end{eqnarray}$ (3) Here u , B _θ, and ϵ are the ion flow velocity, toroidal magnetic field, and two-fluid parameter, respectively. Note that ϵ value is 0.033 in the NSTX experiment, and that 𝜓(r, z) corresponds to the familiar poloidal flux function. For each species (𝛼 = i for the ion fluid and 𝛼 = e for the electron one), the poloidal flow stream function $\bar{{\psi}}_{{\alpha}}$ , generalized Gibbs free energy F _𝛼, and temperature function T _𝛼 are arbitrary surface functions of their respective surface variables. The three arbitrary functions can be assumed so as to reflect the experimental data by choosing appropriate function forms (see Appendix of Ref. [6]). We transform from the coordinates (r, z) to (x, z) by using x ≡ r ²∕2, since Eq. (2) is singular at the position, r = 0. Due to this transformation, Eq. (2) can be rewitten as $\begin{eqnarray}\displaystyle \displaystyle 2x\frac{\partial ^{2}{\psi}}{\partial x^{2}}+\displaystyle \frac{\partial ^{2}{\psi}}{\partial z^{2}}=\displaystyle \frac{\sqrt{2x}}{{\varepsilon}}(B_{{\theta}}\bar{{\psi}}_{e}^{\prime }-nu_{{\theta}})-2nx(F_{e}^{\prime }-T_{e}^{\prime }\ln n). & & \displaystyle\end{eqnarray}$ (4) The computational domain and boundary conditions are shown in Fig. 1.

3. Numerical method

3.1. Computation procedure

Under the above assumptions and boundary conditions, Eqs (1)–(4) are numerically solved by using the following iterative procedure. (1st step) We assume the right-hand side of Eq. (4), and then solve Eq. (4) with the second-order finite difference approximation by using the multi-grid method. (2nd step) We update Y (r, z) and B _θ(r, z) by using $\begin{eqnarray}\displaystyle Y={\psi}+{\varepsilon}ru_{{\theta}}, & & \displaystyle\end{eqnarray}$ (5) and $\begin{eqnarray}\displaystyle B_{{\theta}}=\displaystyle \frac{K({\psi})}{r}+\frac{\bar{{\psi}}_{i}(Y)-\bar{{\psi}}_{i}({\psi})}{{\varepsilon}r}, & & \displaystyle\end{eqnarray}$ (6) where the nearby-fluid model [7], $\begin{eqnarray}\displaystyle \bar{{\psi}}_{i}({\psi})=\bar{{\psi}}_{e}({\psi})+{\varepsilon}K({\psi}), & & \displaystyle\end{eqnarray}$ (7) is employed in the Ampere’s law, $\begin{eqnarray}\displaystyle B_{{\theta}}=\displaystyle \frac{\bar{{\psi}}_{i}(Y)-\bar{{\psi}}_{e}({\psi})}{{\varepsilon}r}. & & \displaystyle\end{eqnarray}$ (8) (3rd step) We update the toroidal ion flow velocity u _θ(r, z) by using Eq. (1) solved for u _θ(r, z), $\begin{eqnarray}\displaystyle u_{{\theta}}=\displaystyle \frac{\bar{{\psi}}_{i}^{\prime }}{n}B_{{\theta}}+{\varepsilon}r\left(F_{i}^{\prime }-T_{i}^{\prime }\ln n\right)-{\varepsilon}\frac{\bar{{\psi}}_{i}^{\prime }}{n}r\boldsymbol{\nabla}\cdot \left(\frac{\bar{{\psi}}_{i}^{\prime }}{n}\frac{\boldsymbol{\nabla}Y}{r^{2}}\right), & & \displaystyle\end{eqnarray}$ (9)through Y (r, z) and B _θ(r, z) obtained in the 2nd step, and then modify Y (r, z) and B _θ(r, z) by using u _θ(r, z) updated above. (4th step) We update the density n (r, z) by using Eq. (3) through the poloidal ion flow velocity, $\begin{eqnarray}\displaystyle \boldsymbol{u}_{p}=\displaystyle \frac{{\nabla}\bar{{\psi}}_{i}}{nr}, & & \displaystyle\end{eqnarray}$ (10) and u _θ(r, z) obtained in the 3rd step. (5th step) We check the convergence of 𝜓(r, z) by using the criteria, $\begin{eqnarray}\displaystyle \displaystyle \frac{\Vert {\psi}^{(k+1)}-{\psi}^{(k)}\Vert }{\Vert {\psi}^{(k+1)}\Vert }\leq {\varepsilon}_{outer}, & & \displaystyle\end{eqnarray}$ (11)and then update the right-hand side of Eq. (4) by using 𝜓(r, z), Y (r, z), B _θ(r, z), u _θ(r, z), and n (r, z). The above steps are iterated until the sufficient convergence has been achieved. Here, k is an iterarion number of the outer loop constructed in this iteration, and ϵ_outer is a value of the convergence criteria. The multi-grid iteration constructs an inner loop.

3.2. Multi-grid method

The multi-grid method (MGM) consists of smoothing step and coarse-grid correction step [4]. The high frequency components of the residual of solution 𝜓(r, z) can be eliminated in the smoothing step. We employ three types of classical iterative method, the damped Jacobi, Gauss-Seidel, and SOR methods as the sweeper to check the smoothing effect. Equation (4) with the second-order finite difference approximation can result in the following matrix form, $\begin{eqnarray}\displaystyle L_{l}\boldsymbol{\it\psi}_{l}=\boldsymbol{f}_{l}. & & \displaystyle\end{eqnarray}$ (12) Here, the subscript l is called level number, and the grid number N _l is changed with the relation, N _l−1 = N _l∕2. Application of 𝜈 times classical iterative method to $\boldsymbol{\it\psi}_{l}^{(m)}$ is denoted by $\begin{eqnarray}\displaystyle \bar{\boldsymbol{\it\psi}}_{l}:=S_{l}^{{\nu}}(\boldsymbol{\it\psi}_{l}^{(m)},\boldsymbol{f}_{l}), & & \displaystyle\end{eqnarray}$ (13) where $\bar{\boldsymbol{\it\psi}}_{l}$ is an iterative solution, and m is the mth iteration of MGM. In order to check the smoothing effect, we define the defect norm, $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle \Vert \boldsymbol{d}_{l}^{({\nu})}\Vert \equiv \Vert L_{l}\bar{\boldsymbol{\it\psi}}_{l}^{({\nu})}-\boldsymbol{f}\Vert ,\end{eqnarray}$ (14) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \displaystyle \Vert \boldsymbol{d}_{l}^{({\nu})}\Vert \equiv \max _{0\leq i\leq N_{l,x}\atop 0\leq j\leq N_{l,z}}|\boldsymbol{d}_{l,i,j}^{({\nu})}|.\end{eqnarray}$ (15) Here, N _{l, x} and N _{l, z} are the grid number of x- and z-directions at level l, respectively. We check the smoothing effect by using the criteria, $\begin{eqnarray}\displaystyle \displaystyle \frac{\Vert \boldsymbol{d}_{l}^{({\nu}+1)}\Vert }{\Vert \boldsymbol{d}_{l}^{({\nu})}\Vert }\geq {\varepsilon}_{d}, & & \displaystyle\end{eqnarray}$ (16) where the value of ϵ_d is set by 0 < ϵ_d <1, because it approaches unity as the iteration is repeated.

The next step is the coarse-grid correction which can eliminate the low frequency components of the residual of solution 𝜓(r, z). Substituting $\bar{\boldsymbol{\it\psi}}_{l}$ that the high frequency components of the residual is eliminated in the smoothing step into the equation L _l𝝍_l − f _l = 0, we obtain the defect, $\begin{eqnarray}\displaystyle \boldsymbol{d}_{l}=L_{l}\bar{\boldsymbol{\it\psi}}_{l}-\boldsymbol{f}_{l} & & \displaystyle\end{eqnarray}$ (17) of $\bar{\boldsymbol{\it\psi}}_{l}$ . The residual v _l is expressed by v _l = $\bar{\boldsymbol{\psi}}_l$ −𝝍_l. Here, 𝝍_l is the exact solution. We obtain the equation for v _l, $\begin{eqnarray}\displaystyle L_{l}\boldsymbol{v}_{l}=\boldsymbol{d}_{l}. & & \displaystyle\end{eqnarray}$ (18) Solving Eq. (18) is equivalent to solving the orignal Eq. (12). The residual v _l obtained in the smoothing step gets a smooth grid function, and can be approximated on coarse grids. To approximate Eq. (18) by a coarse-grid equation $\begin{eqnarray}\displaystyle L_{l-1}\boldsymbol{v}_{l-1}=\boldsymbol{d}_{l-1}, & & \displaystyle\end{eqnarray}$ (19) we have to choose appropriate d _l−1 which is produced by d _l so that v _l−1 can provide the approximation to v _l. Hence, a linear mapping r, called a restriction, defines d _l−1 by $\begin{eqnarray}\displaystyle \boldsymbol{d}_{l-1}\equiv r\boldsymbol{d}_{l}. & & \displaystyle\end{eqnarray}$ (20) The solution v _l−1 of Eq. (19) can be obtained by a direct or iterative method. We have to interpolate v _l−1 to obtain approximation to v _l by $\begin{eqnarray}\displaystyle \tilde{\boldsymbol{v}}_{l}\equiv p\boldsymbol{v}_{l-1}, & & \displaystyle\end{eqnarray}$ (21) where p is called a prolongation which describes a coarse-to-fine interpolation. By using the smoothed iterative solution, $\bar{\boldsymbol{\psi}}_l$ and $\tilde{\boldsymbol{v}}_{l}=∼p\boldsymbol{v}_{l-1}$ , the approximate solution of Eq. (12) is updated by $\begin{eqnarray}\displaystyle \boldsymbol{\it\psi}_{l}^{new}\equiv \bar{\boldsymbol{\it\psi}}_{l}-\tilde{\boldsymbol{v}}_{l}. & & \displaystyle\end{eqnarray}$ (22) The step from $\bar{\boldsymbol{\it\psi}}_{l}$ to $\boldsymbol{\it\psi}_{l}^{new}$ by Eq. (22) is called the coarse-grid correction. The combination of smoothing step and coarse-grid correction is called a two-grid method, since two levels l and l– 1 are involved. The coarse grid equation (19) is the same form as the original equation (12) (only l is replaced by l– 1). The two-grid method with levels $\{l-1,l-2\}$ instead of $\{l,l-1\}$ can be used as an iterative solver of Eq. (19). The combination yields a three-grid method involving the levels $\{l,l-1,l-2\}$ . The two-grid method can be extended to the MGM. The convergence criteria of $\boldsymbol{\it\psi}_{l}^{new}$ is similar to Eq. (11). The V-cycle is adopted as the iteration process of the MGM.

4. Numerical results

Before investigating the properties of numerical convergence of solution 𝜓(r, z), we show contours of 𝜓(r, z) in Fig. 2. The contours indicate the reconstructed equilibrium solution 𝜓(r, z), corresponding to the magnetic flux surfaces in the NSTX device based on the input parameters and surface function coefficients (see Table I of Ref. [6]). The calculation conditions are the MGM (l = 4 → l = 1, V-cycle) with the damped Jacobi method (damped coefficient, ω_d = 0.9 and smoothing criteria, ϵ_d = 0.9), (N _{l, x} N _{l, z}) = (128, 68), the SOR method at l = 1 (acceleration coefficient, ω_s = 1.3 and convergence criteria, ϵ_s = 10⁻⁵), convergence criteria of the MGM, ϵ_MGM = 10⁻⁵, and convergence criteria of the outer loop, ϵ_outer = 10⁻⁵.

Fig. 2.

Typical two-fluid ST equilibrium of NSTX device.

Figure 3 shows dependency of the average value at z = 0 of righ-hand side of Eq. (4) on the grid number. The calculation conditions are the MGM (V-cycle) with the damped Jacobi method (ω_d = 0.9 and ϵ_d = 0.9), the SOR method at l = 1 (ω_s = 1.3 and ϵ_s = 10⁻⁵), ϵ_MGM = 10⁻⁵, and ϵ_outer = 10⁻⁵. This average value converges with inverse square of the grid number. The covergence property of the average value indicates to be numerically stable.

Fig. 3.

Dependency of the average value of righ-hand side of Eq. (4) on the grid number.

Fig. 4.

Dependency of the residual of 𝜓(r, z) on the iteration number of the MGM.

We investigate dependency of the numercal convergence on the sweeper (the damped Jacobi, Gauss-Seidel, and SOR methods) in the smoothing step. Figure 4 shows dependency of the residual of 𝜓(r, z) on the iteration number of the MGM. The calculation conditions are the MGM (l = 4 → l = 1, V-cycle) with the damped Jacobi method (ω_d = 0.9), the MGM (l = 4 → l = 1, V-cycle) with the Gauss-Seidel method, the MGM (l = 4 → l = 1, V-cycle) with the SOR method (ω_SOR = 1.3) (N _{l, x} N _{l, z}) = (128, 168), ϵ_d = 0.9, the SOR method at l = 1 (ω_s = 1.3 and ϵ_s = 10⁻⁵) ϵ_MGM = 10⁻⁵, and ϵ_outer = 10⁻⁵ The residual decreases with the iteration number of the MGM. The residual reduction rate of the damped Jacobi method is larger than that of the Gauss-Seidel and SOR methods. The residual reduction rate of the Gauss-Seidel method is kept at approximately constant, and almost the same as that of the SOR method.

Figure 5 shows dependency of the residual of 𝜓(r, z) on the iteration number of the outer loop. The calculation conditions are the MGM (l = 7 → l = 1, V-cycle) with the damped Jacobi method (ω_d = 0.9), the MGM (l = 7 → l = 1, V-cycle) with the Gauss-Seidel method, the MGM (l = 7 → l = 1, V-cycle) with the SOR method (ω_SOR = 1.3), the SOR method at l = 1 (ω_s = 1.3 and ϵ_s = 10⁻⁵) (N _{l, x}, N _{l, z}) = (1024,1344), ϵ_d = 0.9, ϵ_MGM = 10⁻⁵, and ϵ_outer = 10⁻¹⁰. The residual decreases with the iteration number of the outer loop. The residual reduction rate of the damped Jacobi method is larger than that of the Gauss-Seidel and SOR methods, due to the good covergence of the damped Jacobi method as shown in Fig. 4.

Fig. 5.

Dependency of the residual of 𝜓(r, z) on the iteration number of the outer loop.

Figure 6 shows dependency of the CPU time on the residual of 𝜓(r, z) The calculation conditions are the MGM (l = 6 → l = 1, V-cycle) with the damped Jacobi method (ω_d = 0.9), the MGM (l = 6 → l = 1, V-cycle) with the Gauss-Seidel method, and the MGM (l = 6 → l = 1, V-cycle) with the SOR method (ω_SOR = 1.3), the SOR method at l = 1 (ω_s = 1.3 and ϵ_s = 10⁻⁵), (N _{l, x}, N _{l, z}) = (512, 672)ϵ_d = 0.9, ϵ_MGM = 10⁻⁵, and the SOR method (ω_SOR = 1.5). The CPU time of the MGM with the damped Jacobi method is about 60 times shorter than that of the SOR method at the residual value of 8.7 ×10⁻¹¹. The CPU time of the MGM with the Gauss-Seidel method is almost the same as that of the MGM with the SOR method.

Fig. 6.

Dependency of the CPU time on the residual of 𝜓(r, z).

Figure 7 shows dependency of the CPU time on the grid number. The calculation conditions are the MGM (V-cycle) with the damped Jacobi method (ω_d = 0.9), MGM (V-cycle) with the Gauss-Seidel method, and the MGM (V-cycle) with the SOR method (ω_SOR = 1.3), the SOR method at l = 1 (ω_s = 1.3 and ϵ_s = 10⁻⁵), ϵ_d = 0.9, ϵ_MGM = 10⁻⁵, ϵ_outer = 10⁻⁵, and the SOR method (ω_SOR = 1.5). The CPU time of the MGM with the damped Jacobi method is about 16 times shorter than that of the SOR method at the square of grid number of 1.4 ×10⁶.

Fig. 7.

Dependency of the CPU time on the grid number.

Consequently, the multi-grid method with the damped Jacobi method is effective for computing the two-fluid flowing equilibrium equations with numerical stability, high accuracy, and high speed.

5. Conclusions

We have numerically determined the two-fluid (ions and electrons) flowing equilibrium configurations of spherical torus plasma by using the MGM. The axisymmetric equilibrium equations consist of a pair of second-order partial differential equations for the magnetic (electron) and ion stream functions, and Bernoulli equations for the density. The CPU time of the MGM with the damped Jacobi method is about 60 times shorter than that of the SOR method from the dependency of the CPU time on the residual. The MGM with the damped Jacobi method is effective for solving the two-fluid equilibrium equations with numerical stability, high accuracy, and high speed.

References

Peng

Y.-K.M.

, A component test facility based on the spherical tokamak, Plasma Phys. Controlled Fusion 47 (2005), B263–B283.

Yamada

, Equilibrium analysis of a flowing two-fluid plasma, Phys. Plasmas 9 (2002), 4605–4614.

Steinhauer

L.C.

, Formalism for multi-fluid equilibria with flow, Phys. Plasmas 6 (1999), 2734–2741.

Hackbush

, Multi-Grid Methods and Applications, Springer-Verlag, Berlin, 1985.

Ishida

Steinhauer

L.C.

and Peng

Y.-K.M.

, Two-fluid low-collisionality equilibrium model and application to spherical torus plasmas, Phys. Plasmas 17 (2010), 122507-1–122507-10.

Ishida

and Steinhauer

L.C.

, Spherical torus equilibria reconstructed by a two-fluid low-collisionality model, Phys. Plasmas 19 (2012), 102512-1–102512-10.

Steinhauer

L.C.

and Ishida

, Nearby-fluids equilibria. I. Formalism and transition to single-fluid magnetohydrodynamics, Phys. Plasmas 13 (2006), 052513-1–052513-9.