这是用户在 2024-6-4 23:00 为 file:///E:/paper-fruit/bao19/zh/bao19.html 保存的双语快照页面，由 沉浸式翻译 提供双语支持。了解如何保存？

Global simulation of ion temperature gradient instabilities in a field-reversed configuration
场反转配置中离子温度梯度不稳定性的全局模拟1

Cite as: Phys. Plasmas 26, 042506 (2019); doi: 10.1063/1.5087079
引用： Phys. Plasmas 26， 042506 （2019）;doi： 10.1063/1.5087079Submitted: 27 December 2018. Accepted: 25 March 2019 .
提交日期：2018年12月27日。录用日期： 2019-03-25 .Published Online: 15 April 2019
网络出版日期： 2019-04-15

1. J. Bao, 1 , 2 , 3 , a ) 1 , 2 , 3 , a ) ^(1,2,3,a)){ }^{1,2,3, a)}1,2,3,a) C. K. Lau, 4 4 ^(4){ }^{4}4 (D) Z. Lin, 7 , b ) 7 , b ) ^(7,b)){ }^{7, b)}7,b) (D) H. Y. Wang, 1 , 5 1 , 5 ^(1,5){ }^{1,5}1,5 D. P. Fulton, 4 4 ^(4){ }^{4}4 S. Dettrick, 4 4 ^(4){ }^{4}4 and T. Tajima 4 4 ^(4){ }^{4}4
1. J. Bao， 1 , 2 , 3 , a ) 1 , 2 , 3 , a ) ^(1,2,3,a)){ }^{1,2,3, a)}1,2,3,a) C. K. Lau， 4 4 ^(4){ }^{4}4 （D） Z. Lin， 7 , b ) 7 , b ) ^(7,b)){ }^{7, b)}7,b) （D） H. Y. Wang， 1 , 5 1 , 5 ^(1,5){ }^{1,5}1,5 D. P. Fulton， 4 4 ^(4){ }^{4}4 S. Dettrick， 4 4 ^(4){ }^{4}4 and T. Tajima 4 4 ^(4){ }^{4}4

2. AFFILIATIONS 2. 隶属关系

'University of California, Irvine, California 92697, USA
'加州大学欧文分校， 加利福尼亚州 92697， 美国
2 2 ^(2){ }^{2}2 Beijing National Laboratory for Condensed Matter Physics and CAS Key Laboratory of Soft Matter Physics, Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China
2 2 ^(2){ }^{2}2 中国科学院物理研究所， 北京凝聚态物理国家实验室， 中国科学院软物质物理重点实验室， 北京 100190
3 3 ^(3){ }^{3}3 University of Chinese Academy of Sciences, Beijing 100049, China
3 3 ^(3){ }^{3}3 中国科学院大学， 北京100049， 中国
4 T A E 4 T A E ^(4)TAE{ }^{4} T A E4TAE Technologies, Inc., 19631 Pauling, Foothill Ranch, California 92610, USA
4 T A E 4 T A E ^(4)TAE{ }^{4} T A E4TAE Technologies， Inc.， 19631 Pauling， Foothill Ranch， California 92610， 美国
5 5 ^(5){ }^{5}5 Fusion Simulation Center, Peking University, Beijing 100871, China
5 5 ^(5){ }^{5}5 北京大学聚变模拟中心，北京100871，中国
a)jianbao0626@gmail.com
b) zhihongl@uci.edu b） zhihongl@uci.edu

3. ABSTRACT 3. 摘要

We investigate the global properties of drift waves in the beam driven field-reversed configuration (FRC), the C2-U device, in which the central FRC and its scrape-off layer (SOL) plasma are connected with the formation sections and divertors. The ion temperature gradient modes are globally connected and unstable across these regions, while they are linearly stable inside the FRC separatrix. The unstable global drift waves in the SOL show an axially varying structure that is less intense near the central FRC region and the mirror throat areas, while being more robust in the bad curvature formation exit areas.
我们研究了束流驱动场反转配置（FRC）中漂移波的全局特性，即C2-U器件，其中中央FRC及其刮擦层（SOL）等离子体与地层部分和分流器相连。离子温度梯度模式在这些区域之间是全局相连的，并且不稳定，而它们在FRC分离层内是线性稳定的。SOL中不稳定的全局漂移波呈现轴向变化结构，在FRC中心区域和镜像喉咙区域附近强度较低，而在不良曲率形成出口区域则更为稳健。
Published under license by AIP Publishing. https://doi.org/10.1063/1.5087079
经AIP Publishing许可出版。https://doi.org/10.1063/1.5087079

4. INTRODUCTION 4. 引言

A field-reversed configuration (FRC) is an elongated prolate compact toroid (CT) with magnetic fields predominantly along the poloidal direction, which consists of a core with closed field lines and a scrapeoff layer (SOL) with open field lines. 1 1 ^(1){ }^{1}1 As a fusion reactor concept, FRC has many advantages: The average beta (the ratio between plasma kinetic pressure and magnetic energy density) is close to unity, which suggests a much cheaper fusion energy than the tokamak. The compact shape and simple geometry of FRC also lead to construction convenience and high magnetic efficiency. The SOL region extends to the device ends and forms natural divertors, which are far away from the core region and allow extraction of fusion energy without restriction. TAE Technologies, Inc. has launched a series of FRC experiments. 2 − 7 2 − 7 ^(2-7){ }^{2-7}2−7 A significant energetic ion population generated from neutral beam injection (NBI) can suppress the macroinstability. 8 , 9 8 , 9 ^(8,9){ }^{8,9}8,9 Meanwhile, the large orbit size effects of the energetic particle would not destabilize the microturbulence and the ion scale turbulent transport is suppressed. 10 10 ^(10){ }^{10}10 These experimental efforts lead to the sustainment of beam-driven hot FRC plasmas for more than 5 m s 5 m s 5ms5 \mathrm{~ms}5 ms in C-2U experiments, which is in the confinement regime limited by turbulent transport. In FRC, the open field line SOL region is connected to the closed field line core region, and the turbulence in these two regions affects each other. In recent C 2 − U C 2 − U C2-U\mathrm{C} 2-\mathrm{U}C2−U experiment, it is found that ion-scale turbulence fluctuation is suppressed in the core, while in the SOL, ionand electron-scale turbulence is observed. 10 10 ^(10){ }^{10}10 Thus, it is important to understand the transport mechanism in FRC for the improvement of plasma confinement. 1D and 2D magnetohydrodynamic (MHD) codes have been built up to model the global FRC transport, which requires turbulence simulation codes to provide the transport coeffcients. 11 − 13 11 − 13 ^(11-13)^{11-13}11−13 First-principles particle-in-cell simulation is a powerful tool to study the fusion plasmas combining with theory, which has successful applications in understanding the anomalous transport in tokamak plasmas. 14 − 24 14 − 24 ^(14-24){ }^{14-24}14−24 For the particle-in-cell study of turbulence in FRC, some pioneer works 25 − 27 25 − 27 ^(25-27){ }^{25-27}25−27 have been performed based on the state-of-the-art fusion plasma simulation code: Gyrokinetic Toroidal Code (GTC). GTC has been successfully applied to simulate microturbulence, 21 , 24 21 , 24 ^(21,24){ }^{21,24}21,24 energetic particle transport, 28 28 ^(28){ }^{28}28 Alfven eigenmodes, 29 , 30 29 , 30 ^(29,30){ }^{29,30}29,30 and MHD instabilities 31 , 32 31 , 32 ^(31,32){ }^{31,32}31,32 in toroidal plasmas. An upgrade to the FRC geometry in the Boozer coordinates of GTC has been carried out by Fulton et al., 25 , 26 25 , 26 ^(25,26){ }^{25,26}25,26 and local gyrokinetic particle simulation study by Lau et al. 10 , 27 10 , 27 ^(10,27){ }^{10,27}10,27 shows that the drift wave is stable in the FRC core due to the large orbit size, magnetic well geometry, and short electron transit length. Meanwhile, a new global particle-in-cell FRC code, ANC, has been developed by incorporating core and SOL regions across the separatrix. ANC simulations show that ion scale turbulence can spread from the SOL to the core. 33 , 34 33 , 34 ^(33,34){ }^{33,34}33,34
磁场反转构型 （FRC） 是一种细长的匯实紧凑型环形线圈 （CT），其磁场主要沿极向方向，由具有闭合磁力线的磁芯和具有开放磁力线的刮擦层 （SOL） 组成。 1 1 ^(1){ }^{1}1 作为一种聚变反应堆概念，FRC具有许多优点：平均β（等离子体动压和磁能密度之间的比率）接近统一，这表明聚变能量比托卡马克便宜得多。FRC的紧凑形状和简单的几何形状也带来了施工便利性和高磁效率。SOL区域延伸到器件末端并形成天然分流器，远离核心区域，允许不受限制地提取聚变能。TAE Technologies， Inc. 推出了一系列 FRC 实验。 2 − 7 2 − 7 ^(2-7){ }^{2-7}2−7 中性束注入（NBI）产生的显著高能离子群可以抑制宏观不稳定性。 8 , 9 8 , 9 ^(8,9){ }^{8,9}8,9 同时，高能粒子的大轨道尺寸效应不会破坏微湍流的稳定性，抑制离子尺度的湍流输运。 10 10 ^(10){ }^{10}10 这些实验工作导致光束驱动的热FRC等离子体的维持时间超过 5 m s 5 m s 5ms5 \mathrm{~ms}5 ms C-2U实验，C-2U实验处于受湍流传输限制的约束状态。在FRC中，开放场线SOL区域与封闭场线核心区域相连，这两个区域的湍流相互影响。在最近的 C 2 − U C 2 − U C2-U\mathrm{C} 2-\mathrm{U}C2−U 实验中，发现核心中离子尺度的湍流波动受到抑制，而在溶胶中观察到离子和电子尺度的湍流。 10 10 ^(10){ }^{10}10 因此，了解FRC中改善血浆限制的转运机制非常重要。 已经建立了一维和二维磁流体动力学 （MHD） 代码来模拟全球 FRC 传输，这需要湍流模拟代码来提供传输系数。 11 − 13 11 − 13 ^(11-13)^{11-13}11−13 第一性原理粒子在细胞模拟是研究聚变等离子体与理论相结合的有力工具，在理解托卡马克等离子体的异常输运方面具有成功的应用。 14 − 24 14 − 24 ^(14-24){ }^{14-24}14−24 对于FRC中湍流的粒子细胞研究，基于最先进的聚变等离子体模拟代码：回动环形码（GTC）进行了一些开创性的工作 25 − 27 25 − 27 ^(25-27){ }^{25-27}25−27 。GTC已成功应用于模拟环形等离子体中的微湍流、 21 , 24 21 , 24 ^(21,24){ }^{21,24}21,24 高能粒子输运、 28 28 ^(28){ }^{28}28 Alfven特征模态 29 , 30 29 , 30 ^(29,30){ }^{29,30}29,30 和MHD不稳定性 31 , 32 31 , 32 ^(31,32){ }^{31,32}31,32 。Fulton等人对GTC的Boozer坐标中的FRC几何进行了升级， 25 , 26 25 , 26 ^(25,26){ }^{25,26}25,26 Lau等人 10 , 27 10 , 27 ^(10,27){ }^{10,27}10,27 的局部陀螺动力学粒子模拟研究表明，由于轨道尺寸大、磁井几何形状大、电子传输长度短，漂移波在FRC核心中是稳定的。同时，通过整合整个分离区的核心和 SOL 区域，开发了一种新的全球细胞颗粒 FRC 代码 ANC。ANC模拟表明，离子级湍流可以从SOL扩散到核心。 33 , 34 33 , 34 ^(33,34){ }^{33,34}33,34
In order to study the turbulent transports globally up to the divertor region for FRC, a new GTC family code, gyrokinetic toroidal code-X (GTC-X), is developed in this work by refactoring the coordinate system and geometry of the original GTC code, i.e., change the Boozer coordinates to cylindrical coordinates and change the geometry from tokamak and stellarator to FRC. GTC-X enables the crossseparatrix simulation with a field aligned mesh covering the whole geometry of the FRC. Compared to the original GTC code, both the particle trajectory and the Poisson solver are newly written as well as the simulation grids in the GTC-X code. The GTC mixed-model OpenMP-MPI parallelization 35 35 ^(35){ }^{35}35 is adopted in GTC-X. This paper mainly presents the numerical developments, code verification, and initial results of ion temperature gradient (ITG) modes in the global FRC geometry. GTC-X global simulations show that ITG is unstable in the SOL and stable in the core, which is consistent with previous local simulations and experimental observations. We find that the ITG mode grows along the field line direction in the SOL and shows an axial variation. The maximum amplitude of the ITG mode is in the formation region with bad curvature, while the mode amplitude is small in the central FRC region. The mode structure in the SOL is sensitive to the parallel domain size, which experiences a transition from even parity to odd parity when increasing the domain size.
为了研究FRC全球至分流器区域的湍流传输，通过重构原始GTC代码的坐标系和几何结构，将Boozer坐标更改为圆柱坐标，将几何形状从托卡马克和恒星器更改为FRC，开发了一种新的GTC系列代码，即陀螺动力学环形码-X（GTC-X）。GTC-X 通过覆盖 FRC 整个几何形状的现场对齐网格实现交叉分离模拟。与原始 GTC 代码相比，粒子轨迹和泊松求解器以及 GTC-X 代码中的仿真网格都是新编写的。GTC-X采用GTC混合模型OpenMP-MPI并行化 35 35 ^(35){ }^{35}35 。本文主要介绍了全球FRC几何中离子温度梯度（ITG）模式的数值发展、代码验证和初步结果。GTC-X全局仿真结果显示，ITG在SOL中不稳定，在核心中稳定，这与之前的局部仿真和实验观测结果一致。我们发现 ITG 模式在 SOL 中沿场线方向增长并显示出轴向变化。ITG模式的最大振幅在曲率不好的地层区域，而模态振幅在FRC中心区域较小。SOL 中的模式结构对并行域大小很敏感，当增加域大小时，并行域大小会经历从偶奇偶校验到奇偶校验的过渡。
This paper is organized as follows. In Sec. II, we introduce the global FRC geometry implementation. The gyrokinetic particle simulation model for FRC is described in Sec. III. The benchmark simulation results are shown in Sec. IV. In Sec. V, the global simulation of ITG modes is described. The conclusion is discussed in Sec. VI.
本文的组织结构如下。在第二部分中，我们介绍了全局FRC几何实现。FRC的陀螺动力学粒子模拟模型在第III节中进行了描述。基准仿真结果显示在第 IV 节中。在第五节中，描述了ITG模式的全局仿真。结论在第六节中讨论。

5. GLOBAL FRC GEOMETRY IMPLEMENTATION
5. 全局FRC几何实现

In order to avoid the singularity of magnetic coordinates at the separatrix, 22 , 36 22 , 36 ^(22,36){ }^{22,36}22,36 we adapt the cylindrical coordinate system for global FRC simulation with ( R , ζ , Z ) ( R , ζ , Z ) (R,zeta,Z)(R, \zeta, Z)(R,ζ,Z), where the 3 independent unit vectors satisfy the right hand rule: R ^ × ζ ^ ⋅ Z ^ = 1 R ^ × ζ ^ ⋅ Z ^ = 1 hat(R)xx hat(zeta)* hat(Z)=1\hat{\boldsymbol{R}} \times \hat{\zeta} \cdot \hat{\boldsymbol{Z}}=1R^×ζ^⋅Z^=1. The poloidal magnetic flux ψ ψ psi\psiψ of FRC equilibrium and cylindrical coordinates used in GTC-X is shown in Fig. 1, which is calculated by an axisymmetric force balance FRC equilibrium solver: LR_eqMI code. 37 37 ^(37){ }^{37}37 The equilibrium box size is normalized by the radial position of magnetic axis: R = R 0 = 26.8 c m R = R 0 = 26.8 c m R=R_(0)=26.8cmR=R_{0}=26.8 \mathrm{~cm}R=R0=26.8 cm, i.e., the distance between the magnetic axis and the cylinder axis. There are several mirror plugs in the SOL region aiming at decreasing the particle end loss, and the expanded divertors are located at the ends of open field lines, where we can apply the edge biasing via the plasma-gun electrodes to improve the confinement. 4 4 ^(4){ }^{4}4 In this section, based on the characteristics of FRC equilibrium, we will introduce the algorithms used in GTC-X for global particle-in-cell modeling of FRC.
为了避免磁坐标在分离点处的奇异性， 22 , 36 22 , 36 ^(22,36){ }^{22,36}22,36 我们适配圆柱坐标系进行 ( R , ζ , Z ) ( R , ζ , Z ) (R,zeta,Z)(R, \zeta, Z)(R,ζ,Z) 全局FRC模拟，其中3个独立的单位向量满足右手定则： R ^ × ζ ^ ⋅ Z ^ = 1 R ^ × ζ ^ ⋅ Z ^ = 1 hat(R)xx hat(zeta)* hat(Z)=1\hat{\boldsymbol{R}} \times \hat{\zeta} \cdot \hat{\boldsymbol{Z}}=1R^×ζ^⋅Z^=1 。GTC-X中使用 ψ ψ psi\psiψ 的FRC平衡和圆柱坐标的极化磁通量如图1所示，该磁通量由轴对称力平衡FRC平衡求解器：LR_eqMI代码计算。 37 37 ^(37){ }^{37}37 平衡箱尺寸由磁轴的径向位置归一化： R = R 0 = 26.8 c m R = R 0 = 26.8 c m R=R_(0)=26.8cmR=R_{0}=26.8 \mathrm{~cm}R=R0=26.8 cm ，即磁轴和圆柱轴之间的距离。SOL区域中有几个镜像塞，旨在减少粒子端损耗，扩展的分流器位于开放场线的末端，在那里我们可以通过等离子枪电极应用边缘偏置来改善约束。 4 4 ^(4){ }^{4}4 在本节中，我们将基于FRC平衡的特性，介绍GTC-X中用于FRC全局粒子单元建模的算法。

6. A. Magnetic field representation in cylindrical coordinates
6. A. 圆柱坐标中的磁场表示

The magnetic field and associated derivatives commonly appear in the particle dynamic equations for the simulation of magnetized plasmas; thus, it is important that the magnetic field satisfies ∇ ⋅ B ∇ ⋅ B grad*B\nabla \cdot \boldsymbol{B}∇⋅B = 0 = 0 =0=0=0 numerically. The magnetic coordinates enable the free divergence representation for the magnetic field as B = ∇ α × ∇ β B = ∇ α × ∇ β B=grad alpha xx grad beta\mathbf{B}=\nabla \alpha \times \nabla \betaB=∇α×∇β, where α α alpha\alphaα and β β beta\betaβ are coordinates which vary along the directions orthogonal to the magnetic field. However, magnetic coordinates fail to address the simulation containing different geometric topologies with a separatrix. 36 36 ^(36){ }^{36}36 Thus, we apply the cylindrical coordinate system as the basic coordinates. In order to guarantee the free divergence property for the magnetic field in cylindrical coordinates, we use the poloidal magnetic flux ψ ψ psi\psiψ to calculate the magnetic field components and their derivatives and, thus, enforce the consistency between each component. In FRC, the equilibrium magnetic field B B B\mathbf{B}B has no toroidal component and can be expressed as B = ∇ ψ × ∇ ζ = B R R ^ + B Z Z ^ B = ∇ ψ × ∇ ζ = B R R ^ + B Z Z ^ B=grad psi xx grad zeta=B_(R) hat(R)+B_(Z) hat(Z)\mathbf{B}=\nabla \psi \times \nabla \zeta=B_{R} \hat{\mathbf{R}}+B_{Z} \hat{\mathbf{Z}}B=∇ψ×∇ζ=BRR^+BZZ^. The magnetic field strength in radial and axial directions can then be derived as
磁场和相关导数通常出现在磁化等离子体模拟的粒子动力学方程中;因此，磁场在 = 0 = 0 =0=0=0 数值上满足 ∇ ⋅ B ∇ ⋅ B grad*B\nabla \cdot \boldsymbol{B}∇⋅B 是很重要的。磁坐标使磁场的自由发散表示为 B = ∇ α × ∇ β B = ∇ α × ∇ β B=grad alpha xx grad beta\mathbf{B}=\nabla \alpha \times \nabla \betaB=∇α×∇β ，其中 α α alpha\alphaα 和 β β beta\betaβ 是沿与磁场正交的方向变化的坐标。然而，磁坐标无法解决包含不同几何拓扑的模拟问题。 36 36 ^(36){ }^{36}36 因此，我们应用圆柱坐标系作为基本坐标。为了保证圆柱坐标中磁场的自由发散性，我们使用极化磁通量 ψ ψ psi\psiψ 来计算磁场分量及其导数，从而加强每个分量之间的一致性。在 FRC 中，平衡磁场 B B B\mathbf{B}B 没有环形分量，可以表示为 B = ∇ ψ × ∇ ζ = B R R ^ + B Z Z ^ B = ∇ ψ × ∇ ζ = B R R ^ + B Z Z ^ B=grad psi xx grad zeta=B_(R) hat(R)+B_(Z) hat(Z)\mathbf{B}=\nabla \psi \times \nabla \zeta=B_{R} \hat{\mathbf{R}}+B_{Z} \hat{\mathbf{Z}}B=∇ψ×∇ζ=BRR^+BZZ^ 。径向和轴向的磁场强度可以推导出为
B R = − 1 R ∂ ψ ∂ Z B R = − 1 R ∂ ψ ∂ Z B_(R)=-(1)/(R)(del psi)/(del Z)B_{R}=-\frac{1}{R} \frac{\partial \psi}{\partial Z}BR=−1R∂ψ∂Z
and
B Z = 1 R ∂ ψ ∂ R B Z = 1 R ∂ ψ ∂ R B_(Z)=(1)/(R)(del psi)/(del R)B_{Z}=\frac{1}{R} \frac{\partial \psi}{\partial R}BZ=1R∂ψ∂R
From equilibrium calculated by the LR_eqMI code, we can get the value of poloidal magnetic flux ψ ψ psi\psiψ over the whole FRC geometry on the equal space grids in the ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) plane as shown in Fig. 1 , and L S R = 150 L S R = 150 LSR=150L S R=150LSR=150 and L S Z = 401 L S Z = 401 LSZ=401L S Z=401LSZ=401 are the equilibrium radial and axial grid numbers, respectively. By using the value on coarse equilibrium grids, we can use the quadratic spline function to calculate ψ ψ psi\psiψ at the arbitrary location ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) inside the equilibrium domain as
根据LR_eqMI码计算的平衡，我们可以得到 ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) 平面内等空间网格上整个FRC几何形状的极化磁通量 ψ ψ psi\psiψ 值，如图1所示， L S R = 150 L S R = 150 LSR=150L S R=150LSR=150 分别 L S Z = 401 L S Z = 401 LSZ=401L S Z=401LSZ=401 是平衡径向和轴向网格数。通过使用粗平衡网格上的值，我们可以使用二次样条函数在平衡域内的任意位置 ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) 计算 ψ ψ psi\psiψ 为
FIG. 1. Contour plot of ψ / | ψ 0 | ψ / ψ 0 psi//|psi_(0)|\psi /\left|\psi_{0}\right|ψ/|ψ0| for global FRC geometry, where ψ 0 ψ 0 psi_(0)\psi_{0}ψ0 is the poloidal magnetic flux value at the magnetic axis as shown by the green star. The black solid lines represent the different field lines (contour line of ψ ψ psi\psiψ ), and the red line represents the separatrix. The arrows denote the directions of cylindrical coordinates.
图 1.全局FRC几何形状的 ψ / | ψ 0 | ψ / ψ 0 psi//|psi_(0)|\psi /\left|\psi_{0}\right|ψ/|ψ0| 等值线图，其中 ψ 0 ψ 0 psi_(0)\psi_{0}ψ0 是磁轴处的极化磁通量值，如绿色星形所示。黑色实线表示不同的场线（等值线）， ψ ψ psi\psiψ 红线表示分离线。箭头表示圆柱坐标的方向。
ψ ( R , Z ) = ψ ( 1 , i , j ) + ψ ( 2 , i , j ) Δ R + ψ ( 3 , i , j ) Δ R 2 + ψ ( 4 , i , j ) Δ Z + ψ ( 5 , i , j ) Δ R Δ Z + ψ ( 6 , i , j ) Δ R 2 Δ Z + ψ ( 7 , i , j ) Δ Z 2 + ψ ( 8 , i , j ) Δ R Δ Z 2 + ψ ( 9 , i , j ) Δ R 2 Δ Z 2 , ψ ( R , Z ) = ψ ( 1 , i , j ) + ψ ( 2 , i , j ) Δ R + ψ ( 3 , i , j ) Δ R 2 + ψ ( 4 , i , j ) Δ Z + ψ ( 5 , i , j ) Δ R Δ Z + ψ ( 6 , i , j ) Δ R 2 Δ Z + ψ ( 7 , i , j ) Δ Z 2 + ψ ( 8 , i , j ) Δ R Δ Z 2 + ψ ( 9 , i , j ) Δ R 2 Δ Z 2 , {:[psi(R","Z)=psi(1","i","j)+psi(2","i","j)Delta R+psi(3","i","j)DeltaR^(2)],[+psi(4","i","j)Delta Z+psi(5","i","j)Delta R Delta Z+psi(6","i","j)DeltaR^(2)Delta Z],[+psi(7","i","j)DeltaZ^(2)+psi(8","i","j)Delta R DeltaZ^(2)+psi(9","i","j)DeltaR^(2)DeltaZ^(2)","]:}\begin{aligned} \psi(R, Z) & =\psi(1, i, j)+\psi(2, i, j) \Delta R+\psi(3, i, j) \Delta R^{2} \\ & +\psi(4, i, j) \Delta Z+\psi(5, i, j) \Delta R \Delta Z+\psi(6, i, j) \Delta R^{2} \Delta Z \\ & +\psi(7, i, j) \Delta Z^{2}+\psi(8, i, j) \Delta R \Delta Z^{2}+\psi(9, i, j) \Delta R^{2} \Delta Z^{2}, \end{aligned}ψ(R,Z)=ψ(1,i,j)+ψ(2,i,j)ΔR+ψ(3,i,j)ΔR2+ψ(4,i,j)ΔZ+ψ(5,i,j)ΔRΔZ+ψ(6,i,j)ΔR2ΔZ+ψ(7,i,j)ΔZ2+ψ(8,i,j)ΔRΔZ2+ψ(9,i,j)ΔR2ΔZ2,
where i ∈ [ 1 , L S R ] i ∈ [ 1 , L S R ] i in[1,LSR]i \in[1, L S R]i∈[1,LSR] and j ∈ [ 1 , L S Z ] j ∈ [ 1 , L S Z ] j in[1,LSZ]j \in[1, L S Z]j∈[1,LSZ] are the radial and axial indexes of equilibrium grids, and Δ R = R − R i , Δ Z = Z − Z j , R i ≤ R < R i + 1 Δ R = R − R i , Δ Z = Z − Z j , R i ≤ R < R i + 1 Delta R=R-R_(i),Delta Z=Z-Z_(j),R_(i) <= R < R_(i+1)\Delta R=R-R_{i}, \Delta Z=Z-Z_{j}, R_{i} \leq R<R_{i+1}ΔR=R−Ri,ΔZ=Z−Zj,Ri≤R<Ri+1 and Z j ≤ Z < Z j + 1 Z j ≤ Z < Z j + 1 Z_(j) <= Z < Z_(j+1)Z_{j} \leq Z<Z_{j+1}Zj≤Z<Zj+1.
其中 i ∈ [ 1 , L S R ] i ∈ [ 1 , L S R ] i in[1,LSR]i \in[1, L S R]i∈[1,LSR] 和 j ∈ [ 1 , L S Z ] j ∈ [ 1 , L S Z ] j in[1,LSZ]j \in[1, L S Z]j∈[1,LSZ] 是平衡网格的径向和轴向索引，和 Δ R = R − R i , Δ Z = Z − Z j , R i ≤ R < R i + 1 Δ R = R − R i , Δ Z = Z − Z j , R i ≤ R < R i + 1 Delta R=R-R_(i),Delta Z=Z-Z_(j),R_(i) <= R < R_(i+1)\Delta R=R-R_{i}, \Delta Z=Z-Z_{j}, R_{i} \leq R<R_{i+1}ΔR=R−Ri,ΔZ=Z−Zj,Ri≤R<Ri+1 Z j ≤ Z < Z j + 1 Z j ≤ Z < Z j + 1 Z_(j) <= Z < Z_(j+1)Z_{j} \leq Z<Z_{j+1}Zj≤Z<Zj+1 。
It is straightforward to show that ∇ ⋅ B = 0 ∇ ⋅ B = 0 grad*B=0\nabla \cdot \boldsymbol{B}=0∇⋅B=0 is guaranteed theoretically and numerically
很明显， ∇ ⋅ B = 0 ∇ ⋅ B = 0 grad*B=0\nabla \cdot \boldsymbol{B}=0∇⋅B=0 这在理论上和数值上都是有保证的
∇ ⋅ B = 1 R ∂ ∂ R ( R B R ) + ∂ B Z ∂ Z = 0 . ∇ ⋅ B = 1 R ∂ ∂ R R B R + ∂ B Z ∂ Z = 0 . grad*B=(1)/(R)(del)/(del R)(RB_(R))+(delB_(Z))/(del Z)=0.\nabla \cdot \boldsymbol{B}=\frac{1}{R} \frac{\partial}{\partial R}\left(R B_{R}\right)+\frac{\partial B_{Z}}{\partial Z}=0 .∇⋅B=1R∂∂R(RBR)+∂BZ∂Z=0.

7. B. Field line coordinates for the perturbed field calculation
7. B. 扰动场计算的场线坐标

In magnetized plasmas, the gyrocenter drift motion across the magnetic field is much slower than the parallel motion along the field line; thus, the wave pattern is always anisotropic in the parallel and perpendicular directions with k ∥ ≪ k ⊥ ( k ∥ k ∥ ≪ k ⊥ k ∥ k_(||)≪k_(_|_)(k_(||):}k_{\|} \ll k_{\perp}\left(k_{\|}\right.k∥≪k⊥(k∥and k ⊥ k ⊥ k_(_|_)k_{\perp}k⊥ are the parallel and perpendicular wave vectors). In order to improve the numerical efficiency and accuracy, the field aligned mesh is widely adapted for particle-in-cell simulation of magnetized plasmas, i.e., the grids are aligned along the magnetic field direction with only a small number in the parallel direction, which can dramatically suppress the high k ∥ k ∥ k_(||)k_{\|}k∥noise and save computational cost without sacrificing key physics dominated by small k ∥ ⋅ 23 , 38 k ∥ ⋅ 23 , 38 k_(||)*^(23,38)k_{\|} \cdot{ }^{23,38}k∥⋅23,38 In global FRC simulation, we setup a field aligned mesh in both core and scrape-off layer (SOL) regions across the separatrix in cylindrical coordinates. Due to the fact that the magnetic field is not uniform in FRC, the field aligned mesh is not regular in cylindrical coordinates. For solving perturbed fields as well as particle-grid gather-scatter operation, we create the field line coordinates for core ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) and SOL ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) regions on the poloidal plane, separately, where S c S c S_(c)S_{c}Sc and S S S S S_(S)S_{S}SS represent the normalized field line distances along the magnetic field line direction in core and SOL regions, and the mesh is regular in the corresponding field line coordinates. Thus, in GTC-X, we use two different coordinate systems: cylindrical coordinates and field line coordinates to represent the location.
在磁化等离子体中，陀螺中心在磁场上的漂移运动比沿磁场线的平行运动慢得多;因此，波型在平行和垂直方向上 k ∥ ≪ k ⊥ ( k ∥ k ∥ ≪ k ⊥ k ∥ k_(||)≪k_(_|_)(k_(||):}k_{\|} \ll k_{\perp}\left(k_{\|}\right.k∥≪k⊥(k∥ 始终是各向异性的，并且是 k ⊥ k ⊥ k_(_|_)k_{\perp}k⊥ 平行和垂直的波矢量）。为了提高数值效率和精度，磁场对齐网格被广泛应用于磁化等离子体的粒子在单元模拟中，即网格沿磁场方向排列，平行方向上只有少量数字，可以显著抑制高 k ∥ k ∥ k_(||)k_{\|}k∥ 噪声，节省计算成本，同时又不牺牲以小 k ∥ ⋅ 23 , 38 k ∥ ⋅ 23 , 38 k_(||)*^(23,38)k_{\|} \cdot{ }^{23,38}k∥⋅23,38 为主的关键物理场在全局FRC仿真中，我们在圆柱坐标的分离线的核心和刮擦层（SOL）区域设置了一个场对齐的网格。由于 FRC 中的磁场不均匀，因此圆柱坐标中的磁场对齐网格不规则。为了求解扰动场以及粒子网格聚集散射操作，我们分别在极化平面上创建核心 ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) 和 SOL ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) 区域的场线坐标，其中 S c S c S_(c)S_{c}Sc 和 S S S S S_(S)S_{S}SS 表示核心和 SOL 区域中沿磁力线方向的归一化场线距离，网格在相应的场线坐标中是规则的。因此，在 GTC-X 中，我们使用两种不同的坐标系：圆柱坐标和场线坐标来表示位置。
The simulation domain is different from equilibrium shown in Fig. 1. Because drift wave instabilities and associated transports are anisotropic in perpendicular and parallel directions, we choose the simulation domain based on perpendicular coordinates ψ ψ psi\psiψ (we do not need to consider about ζ ζ zeta\zetaζ domain because it is toroidally symmetric from ( 0 , 2 π ) ( 0 , 2 π ) (0,2pi)(0,2 \pi)(0,2π) ), i.e., the inner boundary in the core and the outer boundary in SOL are labeled by poloidal magnetic flux: ψ 0 ψ 0 psi_(0)\psi_{0}ψ0 and ψ 1 ψ 1 psi_(1)\psi_{1}ψ1. Furthermore, left and right boundaries in the SOL region are given by Z 0 Z 0 Z_(0)Z_{0}Z0 and Z 1 Z 1 Z_(1)Z_{1}Z1, where Z 0 = − Z 1 Z 0 = − Z 1 Z_(0)=-Z_(1)Z_{0}=-Z_{1}Z0=−Z1 is symmetric with respect to outer midplane Z = 0 Z = 0 Z=0Z=0Z=0.
仿真域与图1所示的平衡域不同。由于漂移波的不稳定性和相关的输运在垂直和平行方向上是各向异性的，因此我们根据垂直坐标选择仿真域 ψ ψ psi\psiψ （我们不需要考虑 ζ ζ zeta\zetaζ 域，因为它是环形对称 ( 0 , 2 π ) ( 0 , 2 π ) (0,2pi)(0,2 \pi)(0,2π) 的），即磁芯中的内边界和SOL中的外边界由极化磁通量标记： ψ 0 ψ 0 psi_(0)\psi_{0}ψ0 和 ψ 1 ψ 1 psi_(1)\psi_{1}ψ1 .此外，SOL 区域中的左右边界由 Z 0 Z 0 Z_(0)Z_{0}Z0 和 Z 1 Z 1 Z_(1)Z_{1}Z1 给出，其中 Z 0 = − Z 1 Z 0 = − Z 1 Z_(0)=-Z_(1)Z_{0}=-Z_{1}Z0=−Z1 相对于外中平面 Z = 0 Z = 0 Z=0Z=0Z=0 是对称的。
First, we define S c S c S_(c)S_{c}Sc and S S S S S_(S)S_{S}SS by tracing each field line on each ψ ψ psi\psiψ grid in the core and SOL regions as
首先， S S S S S_(S)S_{S}SS 我们将核心区域和 SOL 区域中每个 ψ ψ psi\psiψ 网格上的每条场线定义为 S c S c S_(c)S_{c}Sc
ψ ( i c ) = ψ 0 + δ ψ c × ( i c − 1 ) , ψ i c = ψ 0 + δ ψ c × i c − 1 , psi(i_(c))=psi_(0)+deltapsi_(c)xx(i_(c)-1),\psi\left(i_{c}\right)=\psi_{0}+\delta \psi_{c} \times\left(i_{c}-1\right),ψ(ic)=ψ0+δψc×(ic−1),
where δ ψ c = ( ψ X − ψ 0 ) / ( l s p c − 1 ) , ψ X δ ψ c = ψ X − ψ 0 / ( l s p c − 1 ) , ψ X deltapsi_(c)=(psi_(X)-psi_(0))//(lspc-1),psi_(X)\delta \psi_{c}=\left(\psi_{X}-\psi_{0}\right) /(l s p c-1), \psi_{X}δψc=(ψX−ψ0)/(lspc−1),ψX is the value of ψ ψ psi\psiψ at the separatrix, l s p c l s p c lspcl s p clspc is the spline resolution in the core region, and 1 ≤ i c ≤ l s p c 1 ≤ i c ≤ l s p c 1 <= i_(c) <= lspc1 \leq i_{c} \leq l s p c1≤ic≤lspc. And
其中 δ ψ c = ( ψ X − ψ 0 ) / ( l s p c − 1 ) , ψ X δ ψ c = ψ X − ψ 0 / ( l s p c − 1 ) , ψ X deltapsi_(c)=(psi_(X)-psi_(0))//(lspc-1),psi_(X)\delta \psi_{c}=\left(\psi_{X}-\psi_{0}\right) /(l s p c-1), \psi_{X}δψc=(ψX−ψ0)/(lspc−1),ψX 是 ψ ψ psi\psiψ 的值 在分离点上， l s p c l s p c lspcl s p clspc 是核心区域中的样条分辨率，以及 1 ≤ i c ≤ l s p c 1 ≤ i c ≤ l s p c 1 <= i_(c) <= lspc1 \leq i_{c} \leq l s p c1≤ic≤lspc 。和
ψ ( i S ) = ψ X + δ ψ S × ( i S − 1 ) , ψ i S = ψ X + δ ψ S × i S − 1 , psi(i_(S))=psi_(X)+deltapsi_(S)xx(i_(S)-1),\psi\left(i_{S}\right)=\psi_{X}+\delta \psi_{S} \times\left(i_{S}-1\right),ψ(iS)=ψX+δψS×(iS−1),
where δ ψ S = ( ψ 1 − ψ X ) / ( l s p s − 1 ) δ ψ S = ψ 1 − ψ X / ( l s p s − 1 ) deltapsi_(S)=(psi_(1)-psi_(X))//(lsps-1)\delta \psi_{S}=\left(\psi_{1}-\psi_{X}\right) /(l s p s-1)δψS=(ψ1−ψX)/(lsps−1), lsps is the spline resolution in SOL region, and 1 ≤ i S ≤ l 1 ≤ i S ≤ l 1 <= i_(S) <= l1 \leq i_{S} \leq l1≤iS≤l sps. Both S c S c S_(c)S_{c}Sc and S S S S S_(S)S_{S}SS are normalized by the field line length at each ψ ψ psi\psiψ grid and range from 0 to 1 . In the core region, S c S c S_(c)S_{c}Sc starts at the outer midplane with Z = 0 Z = 0 Z=0Z=0Z=0 and increases along the clockwise direction, and grids at S c = 0 S c = 0 S_(c)=0S_{c}=0Sc=0 and S c = 1 S c = 1 S_(c)=1S_{c}=1Sc=1 are overlapped. In the SOL region, the parallel coordinate S S S S S_(S)S_{S}SS starts on the left boundary Z 0 Z 0 Z_(0)Z_{0}Z0 and ends on the right boundary Z 1 Z 1 Z_(1)Z_{1}Z1.
其中 δ ψ S = ( ψ 1 − ψ X ) / ( l s p s − 1 ) δ ψ S = ψ 1 − ψ X / ( l s p s − 1 ) deltapsi_(S)=(psi_(1)-psi_(X))//(lsps-1)\delta \psi_{S}=\left(\psi_{1}-\psi_{X}\right) /(l s p s-1)δψS=(ψ1−ψX)/(lsps−1) ，lsps 是 SOL 区域中的样条分辨率，sps 1 ≤ i S ≤ l 1 ≤ i S ≤ l 1 <= i_(S) <= l1 \leq i_{S} \leq l1≤iS≤l .两者都 S c S c S_(c)S_{c}Sc S S S S S_(S)S_{S}SS 由每个 ψ ψ psi\psiψ 网格处的场线长度归一化，范围从 0 到 1。在核心区域中， S c S c S_(c)S_{c}Sc 从外中平面开始，沿顺时针方向增加 Z = 0 Z = 0 Z=0Z=0Z=0 ，网格在 S c = 0 S c = 0 S_(c)=0S_{c}=0Sc=0 和 S c = 1 S c = 1 S_(c)=1S_{c}=1Sc=1 处重叠。在 SOL 区域中，平行坐标 S S S S S_(S)S_{S}SS 从左边界 Z 0 Z 0 Z_(0)Z_{0}Z0 开始，到右边界 Z 1 Z 1 Z_(1)Z_{1}Z1 结束。
Next, considering the properties of geometry topology of core and SOL regions in FRC, forward spline functions S c = S c [ ψ ( R , Z ) S c = S c [ ψ ( R , Z ) S_(c)=S_(c)[psi(R,Z)S_{c}=S_{c}[\psi(R, Z)Sc=Sc[ψ(R,Z), θ ( R , Z ) ] θ ( R , Z ) ] theta(R,Z)]\theta(R, Z)]θ(R,Z)] and S S = S S [ ψ ( R , Z ) , Z ] S S = S S [ ψ ( R , Z ) , Z ] S_(S)=S_(S)[psi(R,Z),Z]S_{S}=S_{S}[\psi(R, Z), Z]SS=SS[ψ(R,Z),Z] are created for the transformation from cylindrical coordinates to field line coordinates, where sin ⁡ ( θ ) sin ⁡ ( θ ) sin(theta)\sin (\theta)sin⁡(θ) = Z / Z 2 + ( R − R 0 ) 2 = Z / Z 2 + R − R 0 2 =Z//sqrt(Z^(2)+(R-R_(0))^(2))=Z / \sqrt{Z^{2}+\left(R-R_{0}\right)^{2}}=Z/Z2+(R−R0)2 is the geometric angle with respect to the magnetic axis position ( R = R 0 , Z = 0 ) R = R 0 , Z = 0 (R=R_(0),Z=0)\left(R=R_{0}, Z=0\right)(R=R0,Z=0). The uniform scale of geometric angle θ θ theta\thetaθ for creating spline function S c ( ψ , θ ) S c ( ψ , θ ) S_(c)(psi,theta)S_{c}(\psi, \theta)Sc(ψ,θ) is
接下来，考虑FRC中磁芯和SOL区域的几何拓扑性质，正向样条函数 S c = S c [ ψ ( R , Z ) S c = S c [ ψ ( R , Z ) S_(c)=S_(c)[psi(R,Z)S_{c}=S_{c}[\psi(R, Z)Sc=Sc[ψ(R,Z) ， θ ( R , Z ) ] θ ( R , Z ) ] theta(R,Z)]\theta(R, Z)]θ(R,Z)] 并 S S = S S [ ψ ( R , Z ) , Z ] S S = S S [ ψ ( R , Z ) , Z ] S_(S)=S_(S)[psi(R,Z),Z]S_{S}=S_{S}[\psi(R, Z), Z]SS=SS[ψ(R,Z),Z] 创建了从圆柱坐标到场线坐标的转换，其中 sin ⁡ ( θ ) sin ⁡ ( θ ) sin(theta)\sin (\theta)sin⁡(θ) = Z / Z 2 + ( R − R 0 ) 2 = Z / Z 2 + R − R 0 2 =Z//sqrt(Z^(2)+(R-R_(0))^(2))=Z / \sqrt{Z^{2}+\left(R-R_{0}\right)^{2}}=Z/Z2+(R−R0)2 是相对于磁轴位置 ( R = R 0 , Z = 0 ) R = R 0 , Z = 0 (R=R_(0),Z=0)\left(R=R_{0}, Z=0\right)(R=R0,Z=0) 的几何角度。用于创建样条函数 S c ( ψ , θ ) S c ( ψ , θ ) S_(c)(psi,theta)S_{c}(\psi, \theta)Sc(ψ,θ) 的几何角度 θ θ theta\thetaθ 的均匀比例为
θ ( j c ) = δ θ × ( j c − 1 ) θ j c = δ θ × j c − 1 theta(j_(c))=delta theta xx(j_(c)-1)\theta\left(j_{c}\right)=\delta \theta \times\left(j_{c}-1\right)θ(jc)=δθ×(jc−1)
where δ θ = 2 π / ( l s t c − 1 ) δ θ = 2 π / ( l s t c − 1 ) delta theta=2pi//(lstc-1)\delta \theta=2 \pi /(l s t c-1)δθ=2π/(lstc−1), lstc is the spline resolution in the θ θ theta\thetaθ direction, and 1 ≤ j c ≤ l s t c 1 ≤ j c ≤ l s t c 1 <= j_(c) <= lstc1 \leq j_{c} \leq l s t c1≤jc≤lstc. By tracing each field line and calculating the S c S c S_(c)S_{c}Sc value on the uniform scales of ψ ψ psi\psiψ and θ θ theta\thetaθ, the spline function S c ( ψ , θ ) S c ( ψ , θ ) S_(c)(psi,theta)S_{c}(\psi, \theta)Sc(ψ,θ) can be derived as
其中 δ θ = 2 π / ( l s t c − 1 ) δ θ = 2 π / ( l s t c − 1 ) delta theta=2pi//(lstc-1)\delta \theta=2 \pi /(l s t c-1)δθ=2π/(lstc−1) ， lstc 是方向上的 θ θ theta\thetaθ 样条分辨率，而 1 ≤ j c ≤ l s t c 1 ≤ j c ≤ l s t c 1 <= j_(c) <= lstc1 \leq j_{c} \leq l s t c1≤jc≤lstc .通过跟踪每条场线并计算 ψ ψ psi\psiψ 和 θ θ theta\thetaθ 的均匀刻度上的 S c S c S_(c)S_{c}Sc 值，样条函数 S c ( ψ , θ ) S c ( ψ , θ ) S_(c)(psi,theta)S_{c}(\psi, \theta)Sc(ψ,θ) 可以导出为
S c ( ψ , θ ) = S c ( 1 , i c , j c ) + S c ( 2 , i c , j c ) Δ ψ + S c ( 3 , i c , j c ) Δ ψ 2 + S c ( 4 , i c , j c ) Δ θ + S c ( 5 , i c , j c ) Δ ψ Δ θ + S c ( 6 , i c , j c ) Δ ψ 2 Δ θ + S c ( 7 , i c , j c ) Δ θ 2 + S c ( 8 , i c , j c ) Δ ψ Δ θ 2 + S c ( 9 , i c , j c ) Δ ψ 2 Δ θ 2 . S c ( ψ , θ ) = S c 1 , i c , j c + S c 2 , i c , j c Δ ψ + S c 3 , i c , j c Δ ψ 2 + S c 4 , i c , j c Δ θ + S c 5 , i c , j c Δ ψ Δ θ + S c 6 , i c , j c Δ ψ 2 Δ θ + S c 7 , i c , j c Δ θ 2 + S c 8 , i c , j c Δ ψ Δ θ 2 + S c 9 , i c , j c Δ ψ 2 Δ θ 2 . {:[S_(c)(psi","theta)=S_(c)(1,i_(c),j_(c))+S_(c)(2,i_(c),j_(c))Delta psi+S_(c)(3,i_(c),j_(c))Deltapsi^(2)],[+S_(c)(4,i_(c),j_(c))Delta theta+S_(c)(5,i_(c),j_(c))Delta psi Delta theta+S_(c)(6,i_(c),j_(c))Deltapsi^(2)Delta theta],[+S_(c)(7,i_(c),j_(c))Deltatheta^(2)+S_(c)(8,i_(c),j_(c))Delta psi Deltatheta^(2)+S_(c)(9,i_(c),j_(c))Deltapsi^(2)Deltatheta^(2).]:}\begin{aligned} S_{c}(\psi, \theta)= & S_{c}\left(1, i_{c}, j_{c}\right)+S_{c}\left(2, i_{c}, j_{c}\right) \Delta \psi+S_{c}\left(3, i_{c}, j_{c}\right) \Delta \psi^{2} \\ & +S_{c}\left(4, i_{c}, j_{c}\right) \Delta \theta+S_{c}\left(5, i_{c}, j_{c}\right) \Delta \psi \Delta \theta+S_{c}\left(6, i_{c}, j_{c}\right) \Delta \psi^{2} \Delta \theta \\ & +S_{c}\left(7, i_{c}, j_{c}\right) \Delta \theta^{2}+S_{c}\left(8, i_{c}, j_{c}\right) \Delta \psi \Delta \theta^{2}+S_{c}\left(9, i_{c}, j_{c}\right) \Delta \psi^{2} \Delta \theta^{2} . \end{aligned}Sc(ψ,θ)=Sc(1,ic,jc)+Sc(2,ic,jc)Δψ+Sc(3,ic,jc)Δψ2+Sc(4,ic,jc)Δθ+Sc(5,ic,jc)ΔψΔθ+Sc(6,ic,jc)Δψ2Δθ+Sc(7,ic,jc)Δθ2+Sc(8,ic,jc)ΔψΔθ2+Sc(9,ic,jc)Δψ2Δθ2.
The uniform scale of Z Z ZZZ for creating spline function S S ( ψ , Z ) S S ( ψ , Z ) S_(S)(psi,Z)S_{S}(\psi, Z)SS(ψ,Z) is
用于创建样条函数 S S ( ψ , Z ) S S ( ψ , Z ) S_(S)(psi,Z)S_{S}(\psi, Z)SS(ψ,Z) 的 Z Z ZZZ 均匀比例为
Z ( j S ) = δ Z S × ( j S − 1 ) Z j S = δ Z S × j S − 1 Z(j_(S))=deltaZ_(S)xx(j_(S)-1)Z\left(j_{S}\right)=\delta Z_{S} \times\left(j_{S}-1\right)Z(jS)=δZS×(jS−1)
where δ Z S = ( Z 1 − Z 0 ) / ( l s z s − 1 ) δ Z S = Z 1 − Z 0 / ( l s z s − 1 ) deltaZ_(S)=(Z_(1)-Z_(0))//(lszs-1)\delta Z_{S}=\left(Z_{1}-Z_{0}\right) /(l s z s-1)δZS=(Z1−Z0)/(lszs−1), lszs is the spline resolution in the Z Z ZZZ direction, and 1 ≤ j S ≤ l s z s 1 ≤ j S ≤ l s z s 1 <= j_(S) <= lszs1 \leq j_{S} \leq l s z s1≤jS≤lszs. By tracing each field line and calculating the S S S S S_(S)S_{S}SS value on the uniform scales of ψ ψ psi\psiψ and Z , S S ( ψ , Z ) Z , S S ( ψ , Z ) Z,S_(S)(psi,Z)Z, S_{S}(\psi, Z)Z,SS(ψ,Z) is given as
其中 δ Z S = ( Z 1 − Z 0 ) / ( l s z s − 1 ) δ Z S = Z 1 − Z 0 / ( l s z s − 1 ) deltaZ_(S)=(Z_(1)-Z_(0))//(lszs-1)\delta Z_{S}=\left(Z_{1}-Z_{0}\right) /(l s z s-1)δZS=(Z1−Z0)/(lszs−1) ，lszs 是方向上的 Z Z ZZZ 样条分辨率，而 1 ≤ j S ≤ l s z s 1 ≤ j S ≤ l s z s 1 <= j_(S) <= lszs1 \leq j_{S} \leq l s z s1≤jS≤lszs 。通过跟踪每条场线并计算 ψ ψ psi\psiψ 和 Z , S S ( ψ , Z ) Z , S S ( ψ , Z ) Z,S_(S)(psi,Z)Z, S_{S}(\psi, Z)Z,SS(ψ,Z) 的均匀刻度上的 S S S S S_(S)S_{S}SS 值，给出为
S S ( ψ , Z ) = S S ( 1 , i S , j S ) + S S ( 2 , i S , j S ) Δ ψ + S S ( 3 , i S , j S ) Δ ψ 2 + S S ( 4 , i S , j S ) Δ Z + S S ( 5 , i S , j S ) Δ ψ Δ Z + S S ( 6 , i S , j S ) Δ ψ 2 Δ Z + S S ( 7 , i S , j S ) Δ Z 2 + S S ( 8 , i S , j S ) Δ ψ Δ Z 2 + S S ( 9 , i S , j S ) Δ ψ 2 Δ Z 2 . S S ( ψ , Z ) = S S 1 , i S , j S + S S 2 , i S , j S Δ ψ + S S 3 , i S , j S Δ ψ 2 + S S 4 , i S , j S Δ Z + S S 5 , i S , j S Δ ψ Δ Z + S S 6 , i S , j S Δ ψ 2 Δ Z + S S 7 , i S , j S Δ Z 2 + S S 8 , i S , j S Δ ψ Δ Z 2 + S S 9 , i S , j S Δ ψ 2 Δ Z 2 . {:[S_(S)(psi","Z)=S_(S)(1,i_(S),j_(S))+S_(S)(2,i_(S),j_(S))Delta psi+S_(S)(3,i_(S),j_(S))Deltapsi^(2)],[+S_(S)(4,i_(S),j_(S))Delta Z+S_(S)(5,i_(S),j_(S))Delta psi Delta Z+S_(S)(6,i_(S),j_(S))Deltapsi^(2)Delta Z],[+S_(S)(7,i_(S),j_(S))DeltaZ^(2)+S_(S)(8,i_(S),j_(S))Delta psi DeltaZ^(2)+S_(S)(9,i_(S),j_(S))Deltapsi^(2)DeltaZ^(2).]:}\begin{aligned} S_{S}(\psi, Z)= & S_{S}\left(1, i_{S}, j_{S}\right)+S_{S}\left(2, i_{S}, j_{S}\right) \Delta \psi+S_{S}\left(3, i_{S}, j_{S}\right) \Delta \psi^{2} \\ & +S_{S}\left(4, i_{S}, j_{S}\right) \Delta Z+S_{S}\left(5, i_{S}, j_{S}\right) \Delta \psi \Delta Z+S_{S}\left(6, i_{S}, j_{S}\right) \Delta \psi^{2} \Delta Z \\ & +S_{S}\left(7, i_{S}, j_{S}\right) \Delta Z^{2}+S_{S}\left(8, i_{S}, j_{S}\right) \Delta \psi \Delta Z^{2}+S_{S}\left(9, i_{S}, j_{S}\right) \Delta \psi^{2} \Delta Z^{2} . \end{aligned}SS(ψ,Z)=SS(1,iS,jS)+SS(2,iS,jS)Δψ+SS(3,iS,jS)Δψ2+SS(4,iS,jS)ΔZ+SS(5,iS,jS)ΔψΔZ+SS(6,iS,jS)Δψ2ΔZ+SS(7,iS,jS)ΔZ2+SS(8,iS,jS)ΔψΔZ2+SS(9,iS,jS)Δψ2ΔZ2.
Then, we create the inverse spline functions R c ( ψ , S c ) R c ψ , S c R_(c)(psi,S_(c))R_{c}\left(\psi, S_{c}\right)Rc(ψ,Sc) and Z c ( ψ , S c ) Z c ψ , S c Z_(c)(psi,S_(c))Z_{c}\left(\psi, S_{c}\right)Zc(ψ,Sc) for the core region, and R S ( ψ , S S ) R S ψ , S S R_(S)(psi,S_(S))R_{S}\left(\psi, S_{S}\right)RS(ψ,SS) and Z S ( ψ , S S ) Z S ψ , S S Z_(S)(psi,S_(S))Z_{S}\left(\psi, S_{S}\right)ZS(ψ,SS) for the SOL region. The uniform scales of S c S c S_(c)S_{c}Sc and S S S S S_(S)S_{S}SS for the spline function are
然后，我们为核心区域 R S ( ψ , S S ) R S ψ , S S R_(S)(psi,S_(S))R_{S}\left(\psi, S_{S}\right)RS(ψ,SS) Z S ( ψ , S S ) Z S ψ , S S Z_(S)(psi,S_(S))Z_{S}\left(\psi, S_{S}\right)ZS(ψ,SS) 创建反样条函数 R c ( ψ , S c ) R c ψ , S c R_(c)(psi,S_(c))R_{c}\left(\psi, S_{c}\right)Rc(ψ,Sc) 和 Z c ( ψ , S c ) Z c ψ , S c Z_(c)(psi,S_(c))Z_{c}\left(\psi, S_{c}\right)Zc(ψ,Sc) 以及 SOL 区域。样条函数 S c S c S_(c)S_{c}Sc 的均匀尺度为 S S S S S_(S)S_{S}SS
S c ( j c ) = δ S c × ( j c − 1 ) S c j c = δ S c × j c − 1 S_(c)(j_(c))=deltaS_(c)xx(j_(c)-1)S_{c}\left(j_{c}\right)=\delta S_{c} \times\left(j_{c}-1\right)Sc(jc)=δSc×(jc−1)
where δ S c = 1 / ( l s s c − 1 ) , l s s c δ S c = 1 / ( l s s c − 1 ) , l s s c deltaS_(c)=1//(lssc-1),lssc\delta S_{c}=1 /(l s s c-1), l s s cδSc=1/(lssc−1),lssc is the spline resolution in the parallel direction in the core region, and 1 ≤ j c ≤ l s s c 1 ≤ j c ≤ l s s c 1 <= j_(c) <= lssc1 \leq j_{c} \leq l s s c1≤jc≤lssc
其中 δ S c = 1 / ( l s s c − 1 ) , l s s c δ S c = 1 / ( l s s c − 1 ) , l s s c deltaS_(c)=1//(lssc-1),lssc\delta S_{c}=1 /(l s s c-1), l s s cδSc=1/(lssc−1),lssc 是核心区域平行方向上的样条分辨率，以及 1 ≤ j c ≤ l s s c 1 ≤ j c ≤ l s s c 1 <= j_(c) <= lssc1 \leq j_{c} \leq l s s c1≤jc≤lssc
S S ( j S ) = δ S S × ( j S − 1 ) , S S j S = δ S S × j S − 1 , S_(S)(j_(S))=deltaS_(S)xx(j_(S)-1),S_{S}\left(j_{S}\right)=\delta S_{S} \times\left(j_{S}-1\right),SS(jS)=δSS×(jS−1),
where δ S S = 1 / ( l s s s − 1 ) δ S S = 1 / ( l s s s − 1 ) deltaS_(S)=1//(lsss-1)\delta S_{S}=1 /(l s s s-1)δSS=1/(lsss−1), lsss is the spline resolution in the parallel direction in the SOL region, and 1 ≤ j S ≤ l s s s 1 ≤ j S ≤ l s s s 1 <= j_(S) <= lsss1 \leq j_{S} \leq l s s s1≤jS≤lsss. It is straightforward to get the values R c [ ψ ( i c ) , S c ( j c ) ] , Z c [ ψ ( i c ) , S c ( j c ) ] , R S [ ψ ( i S ) , S S ( j S ) ] R c ψ i c , S c j c , Z c ψ i c , S c j c , R S ψ i S , S S j S R_(c)[psi(i_(c)),S_(c)(j_(c))],Z_(c)[psi(i_(c)),S_(c)(j_(c))],R_(S)[psi(i_(S)),S_(S)(j_(S))]R_{c}\left[\psi\left(i_{c}\right), S_{c}\left(j_{c}\right)\right], Z_{c}\left[\psi\left(i_{c}\right), S_{c}\left(j_{c}\right)\right], R_{S}\left[\psi\left(i_{S}\right), S_{S}\left(j_{S}\right)\right]Rc[ψ(ic),Sc(jc)],Zc[ψ(ic),Sc(jc)],RS[ψ(iS),SS(jS)], and Z S [ ψ ( i S ) , S S ( j S ) ] Z S ψ i S , S S j S Z_(S)[psi(i_(S)),S_(S)(j_(S))]Z_{S}\left[\psi\left(i_{S}\right), S_{S}\left(j_{S}\right)\right]ZS[ψ(iS),SS(jS)] on the uniform scales of ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) and ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS), and the quadratic spline functions R c ( ψ , S c ) , Z c ( ψ , S c ) , R S ( ψ , S S ) R c ψ , S c , Z c ψ , S c , R S ψ , S S R_(c)(psi,S_(c)),Z_(c)(psi,S_(c)),R_(S)(psi,S_(S))R_{c}\left(\psi, S_{c}\right), Z_{c}\left(\psi, S_{c}\right), R_{S}\left(\psi, S_{S}\right)Rc(ψ,Sc),Zc(ψ,Sc),RS(ψ,SS), and Z S ( ψ , S S ) Z S ψ , S S Z_(S)(psi,S_(S))Z_{S}\left(\psi, S_{S}\right)ZS(ψ,SS) can then be created as
其中 δ S S = 1 / ( l s s s − 1 ) δ S S = 1 / ( l s s s − 1 ) deltaS_(S)=1//(lsss-1)\delta S_{S}=1 /(l s s s-1)δSS=1/(lsss−1) ，lsss 是 SOL 区域中平行方向上的样条分辨率，而 1 ≤ j S ≤ l s s s 1 ≤ j S ≤ l s s s 1 <= j_(S) <= lsss1 \leq j_{S} \leq l s s s1≤jS≤lsss 。在 和 的均匀尺度上，得到值 R c [ ψ ( i c ) , S c ( j c ) ] , Z c [ ψ ( i c ) , S c ( j c ) ] , R S [ ψ ( i S ) , S S ( j S ) ] R c ψ i c , S c j c , Z c ψ i c , S c j c , R S ψ i S , S S j S R_(c)[psi(i_(c)),S_(c)(j_(c))],Z_(c)[psi(i_(c)),S_(c)(j_(c))],R_(S)[psi(i_(S)),S_(S)(j_(S))]R_{c}\left[\psi\left(i_{c}\right), S_{c}\left(j_{c}\right)\right], Z_{c}\left[\psi\left(i_{c}\right), S_{c}\left(j_{c}\right)\right], R_{S}\left[\psi\left(i_{S}\right), S_{S}\left(j_{S}\right)\right]Rc[ψ(ic),Sc(jc)],Zc[ψ(ic),Sc(jc)],RS[ψ(iS),SS(jS)] 和 以及二次样条函数 R c ( ψ , S c ) , Z c ( ψ , S c ) , R S ( ψ , S S ) R c ψ , S c , Z c ψ , S c , R S ψ , S S R_(c)(psi,S_(c)),Z_(c)(psi,S_(c)),R_(S)(psi,S_(S))R_{c}\left(\psi, S_{c}\right), Z_{c}\left(\psi, S_{c}\right), R_{S}\left(\psi, S_{S}\right)Rc(ψ,Sc),Zc(ψ,Sc),RS(ψ,SS) 是很简单的，然后 Z S ( ψ , S S ) Z S ψ , S S Z_(S)(psi,S_(S))Z_{S}\left(\psi, S_{S}\right)ZS(ψ,SS) 可以创建为 ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) Z S [ ψ ( i S ) , S S ( j S ) ] Z S ψ i S , S S j S Z_(S)[psi(i_(S)),S_(S)(j_(S))]Z_{S}\left[\psi\left(i_{S}\right), S_{S}\left(j_{S}\right)\right]ZS[ψ(iS),SS(jS)]
R c ( ψ , S c ) = R c ( 1 , i c , j c ) + R c ( 2 , i c , j c ) Δ ψ + R c ( 3 , i c , j c ) Δ ψ 2 + R c ( 4 , i c , j c ) Δ S + R c ( 5 , i c , j c ) Δ ψ Δ S + R c ( 6 , i c , j c ) Δ ψ 2 Δ S + R c ( 7 , i c , j c ) Δ S 2 + R c ( 8 , i c , j c ) Δ ψ Δ S 2 + R c ( 9 , i c , j c ) Δ ψ 2 Δ S 2 , R c ψ , S c = R c 1 , i c , j c + R c 2 , i c , j c Δ ψ + R c 3 , i c , j c Δ ψ 2 + R c 4 , i c , j c Δ S + R c 5 , i c , j c Δ ψ Δ S + R c 6 , i c , j c Δ ψ 2 Δ S + R c 7 , i c , j c Δ S 2 + R c 8 , i c , j c Δ ψ Δ S 2 + R c 9 , i c , j c Δ ψ 2 Δ S 2 , {:[R_(c)(psi,S_(c))=R_(c)(1,i_(c),j_(c))+R_(c)(2,i_(c),j_(c))Delta psi+R_(c)(3,i_(c),j_(c))Deltapsi^(2)],[+R_(c)(4,i_(c),j_(c))Delta S+R_(c)(5,i_(c),j_(c))Delta psi Delta S+R_(c)(6,i_(c),j_(c))Deltapsi^(2)Delta S],[+R_(c)(7,i_(c),j_(c))DeltaS^(2)+R_(c)(8,i_(c),j_(c))Delta psi DeltaS^(2)+R_(c)(9,i_(c),j_(c))Deltapsi^(2)DeltaS^(2)","]:}\begin{aligned} R_{c}\left(\psi, S_{c}\right)= & R_{c}\left(1, i_{c}, j_{c}\right)+R_{c}\left(2, i_{c}, j_{c}\right) \Delta \psi+R_{c}\left(3, i_{c}, j_{c}\right) \Delta \psi^{2} \\ & +R_{c}\left(4, i_{c}, j_{c}\right) \Delta S+R_{c}\left(5, i_{c}, j_{c}\right) \Delta \psi \Delta S+R_{c}\left(6, i_{c}, j_{c}\right) \Delta \psi^{2} \Delta S \\ & +R_{c}\left(7, i_{c}, j_{c}\right) \Delta S^{2}+R_{c}\left(8, i_{c}, j_{c}\right) \Delta \psi \Delta S^{2}+R_{c}\left(9, i_{c}, j_{c}\right) \Delta \psi^{2} \Delta S^{2}, \end{aligned}Rc(ψ,Sc)=Rc(1,ic,jc)+Rc(2,ic,jc)Δψ+Rc(3,ic,jc)Δψ2+Rc(4,ic,jc)ΔS+Rc(5,ic,jc)ΔψΔS+Rc(6,ic,jc)Δψ2ΔS+Rc(7,ic,jc)ΔS2+Rc(8,ic,jc)ΔψΔS2+Rc(9,ic,jc)Δψ2ΔS2,
and
Z c ( ψ , S c ) = Z c ( 1 , i c , j c ) + Z c ( 2 , i c , j c ) Δ ψ + Z c ( 3 , i c , j c ) Δ ψ 2 + Z c ( 4 , i c , j c ) Δ S + Z c ( 5 , i c , j c ) Δ ψ Δ S + Z c ( 6 , i c , j c ) Δ ψ 2 Δ S + Z c ( 7 , i c , j c ) Δ S 2 + Z c ( 8 , i c , j c ) Δ ψ Δ S 2 + Z c ( 9 , i c , j c ) Δ ψ 2 Δ S 2 , Z c ψ , S c = Z c 1 , i c , j c + Z c 2 , i c , j c Δ ψ + Z c 3 , i c , j c Δ ψ 2 + Z c 4 , i c , j c Δ S + Z c 5 , i c , j c Δ ψ Δ S + Z c 6 , i c , j c Δ ψ 2 Δ S + Z c 7 , i c , j c Δ S 2 + Z c 8 , i c , j c Δ ψ Δ S 2 + Z c 9 , i c , j c Δ ψ 2 Δ S 2 , {:[Z_(c)(psi,S_(c))=Z_(c)(1,i_(c),j_(c))+Z_(c)(2,i_(c),j_(c))Delta psi+Z_(c)(3,i_(c),j_(c))Deltapsi^(2)],[+Z_(c)(4,i_(c),j_(c))Delta S+Z_(c)(5,i_(c),j_(c))Delta psi Delta S+Z_(c)(6,i_(c),j_(c))Deltapsi^(2)Delta S],[+Z_(c)(7,i_(c),j_(c))DeltaS^(2)+Z_(c)(8,i_(c),j_(c))Delta psi DeltaS^(2)+Z_(c)(9,i_(c),j_(c))Deltapsi^(2)DeltaS^(2)","]:}\begin{aligned} Z_{c}\left(\psi, S_{c}\right)= & Z_{c}\left(1, i_{c}, j_{c}\right)+Z_{c}\left(2, i_{c}, j_{c}\right) \Delta \psi+Z_{c}\left(3, i_{c}, j_{c}\right) \Delta \psi^{2} \\ & +Z_{c}\left(4, i_{c}, j_{c}\right) \Delta S+Z_{c}\left(5, i_{c}, j_{c}\right) \Delta \psi \Delta S+Z_{c}\left(6, i_{c}, j_{c}\right) \Delta \psi^{2} \Delta S \\ & +Z_{c}\left(7, i_{c}, j_{c}\right) \Delta S^{2}+Z_{c}\left(8, i_{c}, j_{c}\right) \Delta \psi \Delta S^{2}+Z_{c}\left(9, i_{c}, j_{c}\right) \Delta \psi^{2} \Delta S^{2}, \end{aligned}Zc(ψ,Sc)=Zc(1,ic,jc)+Zc(2,ic,jc)Δψ+Zc(3,ic,jc)Δψ2+Zc(4,ic,jc)ΔS+Zc(5,ic,jc)ΔψΔS+Zc(6,ic,jc)Δψ2ΔS+Zc(7,ic,jc)ΔS2+Zc(8,ic,jc)ΔψΔS2+Zc(9,ic,jc)Δψ2ΔS2,
where Δ ψ = ψ − ψ ( i c ) Δ ψ = ψ − ψ i c Delta psi=psi-psi(i_(c))\Delta \psi=\psi-\psi\left(i_{c}\right)Δψ=ψ−ψ(ic) and Δ S = S c − S ( j c ) ⋅ ψ ( i c ) ≤ ψ < ψ ( i c + 1 ) Δ S = S c − S j c ⋅ ψ i c ≤ ψ < ψ i c + 1 Delta S=S_(c)-S(j_(c))*psi(i_(c)) <= psi < psi(i_(c)+1)\Delta S=S_{c}-S\left(j_{c}\right) \cdot \psi\left(i_{c}\right) \leq \psi<\psi\left(i_{c}+1\right)ΔS=Sc−S(jc)⋅ψ(ic)≤ψ<ψ(ic+1), and S ( j c ) ≤ S c < S ( j c + 1 ) S j c ≤ S c < S j c + 1 S(j_(c)) <= S_(c) < S(j_(c)+1)S\left(j_{c}\right) \leq S_{c}<S\left(j_{c}+1\right)S(jc)≤Sc<S(jc+1).
其中 Δ ψ = ψ − ψ ( i c ) Δ ψ = ψ − ψ i c Delta psi=psi-psi(i_(c))\Delta \psi=\psi-\psi\left(i_{c}\right)Δψ=ψ−ψ(ic) 和 Δ S = S c − S ( j c ) ⋅ ψ ( i c ) ≤ ψ < ψ ( i c + 1 ) Δ S = S c − S j c ⋅ ψ i c ≤ ψ < ψ i c + 1 Delta S=S_(c)-S(j_(c))*psi(i_(c)) <= psi < psi(i_(c)+1)\Delta S=S_{c}-S\left(j_{c}\right) \cdot \psi\left(i_{c}\right) \leq \psi<\psi\left(i_{c}+1\right)ΔS=Sc−S(jc)⋅ψ(ic)≤ψ<ψ(ic+1) 和 S ( j c ) ≤ S c < S ( j c + 1 ) S j c ≤ S c < S j c + 1 S(j_(c)) <= S_(c) < S(j_(c)+1)S\left(j_{c}\right) \leq S_{c}<S\left(j_{c}+1\right)S(jc)≤Sc<S(jc+1) 和 。
R S ( ψ , S S ) = R S ( 1 , i S , j S ) + R S ( 2 , i S , j S ) Δ ψ + R S ( 3 , i S , j S ) Δ ψ 2 + R S ( 4 , i S , j S ) Δ S + R S ( 5 , i S , j S ) Δ ψ Δ S + R S ( 6 , i S , j S ) Δ ψ 2 Δ S + R S ( 7 , i S , j S ) Δ S 2 + R S ( 8 , i S , j S ) Δ ψ Δ S 2 + R S ( 9 , i S , j S ) Δ ψ 2 Δ S 2 , R S ψ , S S = R S 1 , i S , j S + R S 2 , i S , j S Δ ψ + R S 3 , i S , j S Δ ψ 2 + R S 4 , i S , j S Δ S + R S 5 , i S , j S Δ ψ Δ S + R S 6 , i S , j S Δ ψ 2 Δ S + R S 7 , i S , j S Δ S 2 + R S 8 , i S , j S Δ ψ Δ S 2 + R S 9 , i S , j S Δ ψ 2 Δ S 2 , {:[R_(S)(psi,S_(S))=R_(S)(1,i_(S),j_(S))+R_(S)(2,i_(S),j_(S))Delta psi+R_(S)(3,i_(S),j_(S))Deltapsi^(2)],[+R_(S)(4,i_(S),j_(S))Delta S+R_(S)(5,i_(S),j_(S))Delta psi Delta S+R_(S)(6,i_(S),j_(S))Deltapsi^(2)Delta S],[+R_(S)(7,i_(S),j_(S))DeltaS^(2)+R_(S)(8,i_(S),j_(S))Delta psi DeltaS^(2)+R_(S)(9,i_(S),j_(S))Deltapsi^(2)DeltaS^(2)","]:}\begin{aligned} R_{S}\left(\psi, S_{S}\right)= & R_{S}\left(1, i_{S}, j_{S}\right)+R_{S}\left(2, i_{S}, j_{S}\right) \Delta \psi+R_{S}\left(3, i_{S}, j_{S}\right) \Delta \psi^{2} \\ & +R_{S}\left(4, i_{S}, j_{S}\right) \Delta S+R_{S}\left(5, i_{S}, j_{S}\right) \Delta \psi \Delta S+R_{S}\left(6, i_{S}, j_{S}\right) \Delta \psi^{2} \Delta S \\ & +R_{S}\left(7, i_{S}, j_{S}\right) \Delta S^{2}+R_{S}\left(8, i_{S}, j_{S}\right) \Delta \psi \Delta S^{2}+R_{S}\left(9, i_{S}, j_{S}\right) \Delta \psi^{2} \Delta S^{2}, \end{aligned}RS(ψ,SS)=RS(1,iS,jS)+RS(2,iS,jS)Δψ+RS(3,iS,jS)Δψ2+RS(4,iS,jS)ΔS+RS(5,iS,jS)ΔψΔS+RS(6,iS,jS)Δψ2ΔS+RS(7,iS,jS)ΔS2+RS(8,iS,jS)ΔψΔS2+RS(9,iS,jS)Δψ2ΔS2,
and
Z S ( ψ , S S ) = Z S ( 1 , i S , j S ) + Z S ( 2 , i S , j S ) Δ ψ + Z S ( 3 , i S , j S ) Δ ψ 2 + Z S ( 4 , i S , j S ) Δ S + Z S ( 5 , i S , j S ) Δ ψ Δ S + Z S ( 6 , i S , j S ) Δ ψ 2 Δ S + Z S ( 7 , i S , j S ) Δ S 2 + Z S ( 8 , i S , j S ) Δ ψ Δ S 2 + Z S ( 9 , i S , j S ) Δ ψ 2 Δ S 2 , Z S ψ , S S = Z S 1 , i S , j S + Z S 2 , i S , j S Δ ψ + Z S 3 , i S , j S Δ ψ 2 + Z S 4 , i S , j S Δ S + Z S 5 , i S , j S Δ ψ Δ S + Z S 6 , i S , j S Δ ψ 2 Δ S + Z S 7 , i S , j S Δ S 2 + Z S 8 , i S , j S Δ ψ Δ S 2 + Z S 9 , i S , j S Δ ψ 2 Δ S 2 , {:[Z_(S)(psi,S_(S))=Z_(S)(1,i_(S),j_(S))+Z_(S)(2,i_(S),j_(S))Delta psi+Z_(S)(3,i_(S),j_(S))Deltapsi^(2)],[+Z_(S)(4,i_(S),j_(S))Delta S+Z_(S)(5,i_(S),j_(S))Delta psi Delta S+Z_(S)(6,i_(S),j_(S))Deltapsi^(2)Delta S],[+Z_(S)(7,i_(S),j_(S))DeltaS^(2)+Z_(S)(8,i_(S),j_(S))Delta psi DeltaS^(2)+Z_(S)(9,i_(S),j_(S))Deltapsi^(2)DeltaS^(2)","]:}\begin{aligned} Z_{S}\left(\psi, S_{S}\right)= & Z_{S}\left(1, i_{S}, j_{S}\right)+Z_{S}\left(2, i_{S}, j_{S}\right) \Delta \psi+Z_{S}\left(3, i_{S}, j_{S}\right) \Delta \psi^{2} \\ & +Z_{S}\left(4, i_{S}, j_{S}\right) \Delta S+Z_{S}\left(5, i_{S}, j_{S}\right) \Delta \psi \Delta S+Z_{S}\left(6, i_{S}, j_{S}\right) \Delta \psi^{2} \Delta S \\ & +Z_{S}\left(7, i_{S}, j_{S}\right) \Delta S^{2}+Z_{S}\left(8, i_{S}, j_{S}\right) \Delta \psi \Delta S^{2}+Z_{S}\left(9, i_{S}, j_{S}\right) \Delta \psi^{2} \Delta S^{2}, \end{aligned}ZS(ψ,SS)=ZS(1,iS,jS)+ZS(2,iS,jS)Δψ+ZS(3,iS,jS)Δψ2+ZS(4,iS,jS)ΔS+ZS(5,iS,jS)ΔψΔS+ZS(6,iS,jS)Δψ2ΔS+ZS(7,iS,jS)ΔS2+ZS(8,iS,jS)ΔψΔS2+ZS(9,iS,jS)Δψ2ΔS2,
where Δ ψ = ψ − ψ ( i S ) Δ ψ = ψ − ψ i S Delta psi=psi-psi(i_(S))\Delta \psi=\psi-\psi\left(i_{S}\right)Δψ=ψ−ψ(iS) and Δ S = S S − S ( j S ) . ψ ( i S ) ≤ ψ < ψ ( i S + 1 ) Δ S = S S − S j S . ψ i S ≤ ψ < ψ i S + 1 Delta S=S_(S)-S(j_(S)).psi(i_(S)) <= psi < psi(i_(S)+1)\Delta S=S_{S}-S\left(j_{S}\right) . \psi\left(i_{S}\right) \leq \psi<\psi\left(i_{S}+1\right)ΔS=SS−S(jS).ψ(iS)≤ψ<ψ(iS+1) and S ( j S ) ≤ S S < S ( j S + 1 ) S j S ≤ S S < S j S + 1 S(j_(S)) <= S_(S) < S(j_(S)+1)S\left(j_{S}\right) \leq S_{S}<S\left(j_{S}+1\right)S(jS)≤SS<S(jS+1). By using Eqs. (7)-(10), we could compute the field aligned mesh in cylindrical coordinates with a given regular mesh in the core region ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) and in the SOL region ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS), respectively. An example of global field aligned mesh is given in Fig. 2(a). The field aligned grids are irregular in cylindrical coordinates. For the overlap part of core and SOL regions at the separatrix ( R > 0 ) ( R > 0 ) (R > 0)(R>0)(R>0), the grid positions are determined by using R c ( ψ , S c ) R c ψ , S c R_(c)(psi,S_(c))R_{c}\left(\psi, S_{c}\right)Rc(ψ,Sc) and Z c ( ψ , S c ) Z c ψ , S c Z_(c)(psi,S_(c))Z_{c}\left(\psi, S_{c}\right)Zc(ψ,Sc), which are shared by both core and SOL regions with the scopes of two coordinates: ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) and ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) as shown by the black stars. In magnetic coordinates, the grids in SOL and core regions are regular inside each domain, respectively. The shared grids at the separatrix are designed regularly with the interior grids in the core region (keep the same parallel coordinate S c S c S_(c)S_{c}Sc value), which are shown in Figs. 2(b) and 2(c).
其中 Δ ψ = ψ − ψ ( i S ) Δ ψ = ψ − ψ i S Delta psi=psi-psi(i_(S))\Delta \psi=\psi-\psi\left(i_{S}\right)Δψ=ψ−ψ(iS) 和 Δ S = S S − S ( j S ) . ψ ( i S ) ≤ ψ < ψ ( i S + 1 ) Δ S = S S − S j S . ψ i S ≤ ψ < ψ i S + 1 Delta S=S_(S)-S(j_(S)).psi(i_(S)) <= psi < psi(i_(S)+1)\Delta S=S_{S}-S\left(j_{S}\right) . \psi\left(i_{S}\right) \leq \psi<\psi\left(i_{S}+1\right)ΔS=SS−S(jS).ψ(iS)≤ψ<ψ(iS+1) 和 S ( j S ) ≤ S S < S ( j S + 1 ) S j S ≤ S S < S j S + 1 S(j_(S)) <= S_(S) < S(j_(S)+1)S\left(j_{S}\right) \leq S_{S}<S\left(j_{S}+1\right)S(jS)≤SS<S(jS+1) .通过使用方程。（7）-（10）中，我们可以分别计算在核心区域 ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) 和SOL区域 ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) 中给定规则网格的圆柱坐标中的场对齐网格。图2（a）给出了全局场对齐网格的示例。场对齐的网格在圆柱坐标中是不规则的。对于分离点 ( R > 0 ) ( R > 0 ) (R > 0)(R>0)(R>0) 处核心区域和 SOL 区域的重叠部分，网格位置由 和 R c ( ψ , S c ) R c ψ , S c R_(c)(psi,S_(c))R_{c}\left(\psi, S_{c}\right)Rc(ψ,Sc) 确定，它们由核心区域和 Z c ( ψ , S c ) Z c ψ , S c Z_(c)(psi,S_(c))Z_{c}\left(\psi, S_{c}\right)Zc(ψ,Sc) SOL 区域共享，范围为两个坐标： ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) 和 ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) ，如黑色星号所示。在磁坐标中，SOL 和核心区域中的网格分别在每个域内是规则的。分离点处的共享网格与核心区域的内部网格有规律地设计（保持相同的平行坐标 S c S c S_(c)S_{c}Sc 值），如图2（b）和图2（c）所示。
In GTC-X, particle dynamic equations are evolved in cylindrical coordinates ( R , ζ , Z ) ( R , ζ , Z ) (R,zeta,Z)(R, \zeta, Z)(R,ζ,Z) to avoid the singularity at the separatrix. However, the field aligned mesh is irregular in ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) space, which is difficult to carry out particle-grid gather-scatter operation for PIC simulation. It is noted that the field aligned grids are labeled by both cylindrical and magnetic coordinates: ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) and ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) in the core region and ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) and ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) in the SOL region. The mesh is regular in magnetic coordinates in the core and SOL regions, respectively, as shown in Fig. 2. Thus, we can transform the particle coordinate from ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) to ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) and ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) by using Eqs. (3), (5), and (6), and then a simple linear interpolation between the particle location and the regular rectangular mesh can be used for particle-grid gather-scatter operation,
在GTC-X中，粒子动力学方程在圆柱坐标中演化， ( R , ζ , Z ) ( R , ζ , Z ) (R,zeta,Z)(R, \zeta, Z)(R,ζ,Z) 以避免分离点处的奇点。然而，场对齐网格在空间上 ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) 是不规则的，难以进行粒子-网格聚集-散射操作进行PIC仿真。需要注意的是，磁场对齐的网格由圆柱坐标和磁坐标标记： ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) 在核心区域和 ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) SOL 区域。网格分别在核心区域和 SOL 区域的磁坐标中呈规则，如图 2 所示。因此，我们可以使用方程将粒子坐标从 ( R , Z ) ( R , Z ) (R,Z)(R, Z)(R,Z) 到 和 ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) 通过 进行转换。（3）、（5）和（6），然后使用粒子位置和规则矩形网格之间的简单线性插值进行粒子网格聚集-散射操作，
FIG. 2. (a) Global field aligned mesh in cylindrical coordinates. Field aligned mesh mapping from cylindrical coordinates to magnetic coordinates: (b) core region grids in ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) coordinates and ( c ) ( c ) (c)(c)(c) SOL region grids in ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) coordinates. The black stars represent the shared grids at the separatrix. The grids shown here are only for illustrating the algorithm, which are much sparser than the ones used in realistic simulation. where we implement 31 D 31 D 31D31 \mathrm{D}31D linear interpolations along ψ ψ psi\psiψ, ζ ζ zeta\zetaζ, and S c / S S S c / S S S_(c)//S_(S)S_{c} / S_{S}Sc/SS for the simulations in 3 dimensional ( ψ , ζ , S c / S S ) ψ , ζ , S c / S S (psi,zeta,S_(c)//S_(S))\left(\psi, \zeta, S_{c} / S_{S}\right)(ψ,ζ,Sc/SS) space.
图 2.（a） 圆柱坐标中的全局场对齐网格。从圆柱坐标到磁坐标的场对齐网格映射：（b） 坐标中 ( ψ , S c ) ψ , S c (psi,S_(c))\left(\psi, S_{c}\right)(ψ,Sc) 的核心区域网格和 ( c ) ( c ) (c)(c)(c) 坐标中的 ( ψ , S S ) ψ , S S (psi,S_(S))\left(\psi, S_{S}\right)(ψ,SS) SOL 区域网格。黑色星星代表分离点的共享网格。此处显示的网格仅用于说明算法，这些网格比实际模拟中使用的网格要稀疏得多。我们实现沿 ψ ψ psi\psiψ 、 ζ ζ zeta\zetaζ 和 S c / S S S c / S S S_(c)//S_(S)S_{c} / S_{S}Sc/SS 的 31 D 31 D 31D31 \mathrm{D}31D 线性插值，用于 3 维 ( ψ , ζ , S c / S S ) ψ , ζ , S c / S S (psi,zeta,S_(c)//S_(S))\left(\psi, \zeta, S_{c} / S_{S}\right)(ψ,ζ,Sc/SS) 空间的模拟。

8. Laplacian operator 8. 拉普拉斯算子

In the global simulation of FRC geometry with the field aligned mesh, the perpendicular Laplacian operator is discretized in ( ψ , ζ , S c ) ψ , ζ , S c (psi,zeta,S_(c))\left(\psi, \zeta, S_{c}\right)(ψ,ζ,Sc) space for the core region and in ( ψ , ζ , S S ) ψ , ζ , S S (psi,zeta,S_(S))\left(\psi, \zeta, S_{S}\right)(ψ,ζ,SS) space for the SOL region, respectively, since the irregular mesh in cylindrical coordinates becomes regular in magnetic coordinates.
在具有场对齐网格的 FRC 几何体的全局仿真中，垂直拉普拉斯算子分别在 ( ψ , ζ , S c ) ψ , ζ , S c (psi,zeta,S_(c))\left(\psi, \zeta, S_{c}\right)(ψ,ζ,Sc) 核心区域的空间和 SOL 区域的空间 ( ψ , ζ , S S ) ψ , ζ , S S (psi,zeta,S_(S))\left(\psi, \zeta, S_{S}\right)(ψ,ζ,SS) 中离散化，因为圆柱坐标中的不规则网格在磁坐标中变得规则。
The Laplacian operator can be expanded in a generalized coordinates
拉普拉斯算子可以在广义坐标中展开
∇ 2 f = 1 J ∂ ∂ ξ α ( J ∇ ξ α ⋅ ∇ ξ β ∂ f ∂ ξ β ) ∇ 2 f = 1 J ∂ ∂ ξ α J ∇ ξ α ⋅ ∇ ξ β ∂ f ∂ ξ β grad^(2)f=(1)/(J)(del)/(delxi^(alpha))(J gradxi^(alpha)*gradxi^(beta)(del f)/(delxi^(beta)))\nabla^{2} f=\frac{1}{J} \frac{\partial}{\partial \xi^{\alpha}}\left(J \nabla \xi^{\alpha} \cdot \nabla \xi^{\beta} \frac{\partial f}{\partial \xi^{\beta}}\right)∇2f=1J∂∂ξα(J∇ξα⋅∇ξβ∂f∂ξβ)
where f f fff represents an arbitrary scalar field, ξ α , ξ β ξ α , ξ β xi^(alpha),xi^(beta)\xi^{\alpha}, \xi^{\beta}ξα,ξβ refer to the three dimensional coordinates, and J J JJJ is Jacobian and defined as J − 1 ϵ α β γ J − 1 ϵ α β γ J^(-1)epsilon^(alpha beta gamma)J^{-1} \epsilon^{\alpha \beta \gamma}J−1ϵαβγ = ∇ ξ α × ∇ ξ β ⋅ ∇ ξ γ , ∇ ξ α = ∂ ξ α ∂ R R ^ + ∂ ξ α ∂ ζ 1 R ζ ^ + ∂ ξ α ∂ Z Z ^ = ∇ ξ α × ∇ ξ β ⋅ ∇ ξ γ , ∇ ξ α = ∂ ξ α ∂ R R ^ + ∂ ξ α ∂ ζ 1 R ζ ^ + ∂ ξ α ∂ Z Z ^ =gradxi^(alpha)xx gradxi^(beta)*gradxi^(gamma),gradxi^(alpha)=(delxi^(alpha))/(del R) hat(R)+(delxi^(alpha))/(del zeta)(1)/(R) hat(zeta)+(delxi^(alpha))/(del Z) hat(Z)=\nabla \xi^{\alpha} \times \nabla \xi^{\beta} \cdot \nabla \xi^{\gamma}, \nabla \xi^{\alpha}=\frac{\partial \xi^{\alpha}}{\partial R} \hat{R}+\frac{\partial \xi^{\alpha}}{\partial \zeta} \frac{1}{R} \hat{\zeta}+\frac{\partial \xi^{\alpha}}{\partial Z} \hat{Z}=∇ξα×∇ξβ⋅∇ξγ,∇ξα=∂ξα∂RR^+∂ξα∂ζ1Rζ^+∂ξα∂ZZ^, and ϵ α β γ = 1 ϵ α β γ = 1 epsilon^(alpha beta gamma)=1\epsilon^{\alpha \beta \gamma}=1ϵαβγ=1, and it becomes antisymmetric when the indices change. 39 39 ^(39){ }^{39}39
其中 f f fff 表示任意标量场， ξ α , ξ β ξ α , ξ β xi^(alpha),xi^(beta)\xi^{\alpha}, \xi^{\beta}ξα,ξβ 指三维坐标，是 J J JJJ 雅可比坐标，定义为 J − 1 ϵ α β γ J − 1 ϵ α β γ J^(-1)epsilon^(alpha beta gamma)J^{-1} \epsilon^{\alpha \beta \gamma}J−1ϵαβγ = ∇ ξ α × ∇ ξ β ⋅ ∇ ξ γ , ∇ ξ α = ∂ ξ α ∂ R R ^ + ∂ ξ α ∂ ζ 1 R ζ ^ + ∂ ξ α ∂ Z Z ^ = ∇ ξ α × ∇ ξ β ⋅ ∇ ξ γ , ∇ ξ α = ∂ ξ α ∂ R R ^ + ∂ ξ α ∂ ζ 1 R ζ ^ + ∂ ξ α ∂ Z Z ^ =gradxi^(alpha)xx gradxi^(beta)*gradxi^(gamma),gradxi^(alpha)=(delxi^(alpha))/(del R) hat(R)+(delxi^(alpha))/(del zeta)(1)/(R) hat(zeta)+(delxi^(alpha))/(del Z) hat(Z)=\nabla \xi^{\alpha} \times \nabla \xi^{\beta} \cdot \nabla \xi^{\gamma}, \nabla \xi^{\alpha}=\frac{\partial \xi^{\alpha}}{\partial R} \hat{R}+\frac{\partial \xi^{\alpha}}{\partial \zeta} \frac{1}{R} \hat{\zeta}+\frac{\partial \xi^{\alpha}}{\partial Z} \hat{Z}=∇ξα×∇ξβ⋅∇ξγ,∇ξα=∂ξα∂RR^+∂ξα∂ζ1Rζ^+∂ξα∂ZZ^ 和 ϵ α β γ = 1 ϵ α β γ = 1 epsilon^(alpha beta gamma)=1\epsilon^{\alpha \beta \gamma}=1ϵαβγ=1 ，当索引发生变化时，它变为反对称。 39 39 ^(39){ }^{39}39
In the core region with ( ψ , ζ , S c ) ψ , ζ , S c (psi,zeta,S_(c))\left(\psi, \zeta, S_{c}\right)(ψ,ζ,Sc), considering S c = S c [ ψ ( R , Z ) , θ ( R , Z ) ] S c = S c [ ψ ( R , Z ) , θ ( R , Z ) ] S_(c)=S_(c)[psi(R,Z),theta(R,Z)]S_{c}=S_{c}[\psi(R, Z), \theta(R, Z)]Sc=Sc[ψ(R,Z),θ(R,Z)], the Laplacian is written as
在核心区域中 ( ψ , ζ , S c ) ψ , ζ , S c (psi,zeta,S_(c))\left(\psi, \zeta, S_{c}\right)(ψ,ζ,Sc) ，考虑到 S c = S c [ ψ ( R , Z ) , θ ( R , Z ) ] S c = S c [ ψ ( R , Z ) , θ ( R , Z ) ] S_(c)=S_(c)[psi(R,Z),theta(R,Z)]S_{c}=S_{c}[\psi(R, Z), \theta(R, Z)]Sc=Sc[ψ(R,Z),θ(R,Z)] ，拉普拉斯被写成
∇ c 2 f = g ψ ψ ∂ 2 f ∂ ψ 2 + 2 g c ψ S ∂ 2 f ∂ ψ ∂ S c + g c S S ∂ 2 f ∂ S c 2 + g ζ ζ ∂ 2 f ∂ ζ 2 + 1 J c [ ∂ ( J c g ψ ψ ) ∂ ψ + ∂ ( J c g c ψ S ) ∂ S c ] ∂ f ∂ ψ + 1 J c [ ∂ ( J c g c ψ S ) ∂ ψ + ∂ ( J c g c S S ) ∂ S c ] ∂ f ∂ S c , ∇ c 2 f = g ψ ψ ∂ 2 f ∂ ψ 2 + 2 g c ψ S ∂ 2 f ∂ ψ ∂ S c + g c S S ∂ 2 f ∂ S c 2 + g ζ ζ ∂ 2 f ∂ ζ 2 + 1 J c ∂ J c g ψ ψ ∂ ψ + ∂ J c g c ψ S ∂ S c ∂ f ∂ ψ + 1 J c ∂ J c g c ψ S ∂ ψ + ∂ J c g c S S ∂ S c ∂ f ∂ S c , {:[grad_(c)^(2)f=g^(psi psi)(del^(2)f)/(delpsi^(2))+2g_(c)^(psi S)(del^(2)f)/(del psi delS_(c))+g_(c)^(SS)(del^(2)f)/(delS_(c)^(2))+g^(zeta zeta)(del^(2)f)/(delzeta^(2))],[+(1)/(J_(c))[(del(J_(c)g^(psi psi)))/(del psi)+(del(J_(c)g_(c)psi S))/(delS_(c))](del f)/(del psi)],[+(1)/(J_(c))[(del(J_(c)g_(c)^(psi S)))/(del psi)+(del(J_(c)g_(c)^(SS)))/(delS_(c))](del f)/(delS_(c))","]:}\begin{aligned} \nabla_{c}^{2} f= & g^{\psi \psi} \frac{\partial^{2} f}{\partial \psi^{2}}+2 g_{c}^{\psi S} \frac{\partial^{2} f}{\partial \psi \partial S_{c}}+g_{c}^{S S} \frac{\partial^{2} f}{\partial S_{c}^{2}}+g^{\zeta \zeta} \frac{\partial^{2} f}{\partial \zeta^{2}} \\ & +\frac{1}{J_{c}}\left[\frac{\partial\left(J_{c} g^{\psi \psi}\right)}{\partial \psi}+\frac{\partial\left(J_{c} g_{c} \psi S\right)}{\partial S_{c}}\right] \frac{\partial f}{\partial \psi} \\ & +\frac{1}{J_{c}}\left[\frac{\partial\left(J_{c} g_{c}^{\psi S}\right)}{\partial \psi}+\frac{\partial\left(J_{c} g_{c}^{S S}\right)}{\partial S_{c}}\right] \frac{\partial f}{\partial S_{c}}, \end{aligned}∇c2f=gψψ∂2f∂ψ2+2gcψS∂2f∂ψ∂Sc+gcSS∂2f∂Sc2+gζζ∂2f∂ζ2+1Jc[∂(Jcgψψ)∂ψ+∂(JcgcψS)∂Sc]∂f∂ψ+1Jc[∂(JcgcψS)∂ψ+∂(JcgcSS)∂Sc]∂f∂Sc,
where 哪里
J c − 1 = ∇ ψ × ∇ ζ ⋅ ∇ S c = ( ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ ) × 1 R ζ ^ ⋅ [ ∂ S c ∂ ψ ( ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ ) + ∂ S c ∂ θ ( ∂ θ ∂ R R ^ + ∂ θ ∂ Z Z ^ ) ] = 1 R ∂ ψ ∂ R ∂ S c ∂ θ ∂ θ ∂ Z − 1 R ∂ ψ ∂ Z ∂ S c ∂ θ ∂ θ ∂ R , g ψ ψ = ∇ ψ ⋅ ∇ ψ = ( ∂ ψ ∂ R ) 2 + ( ∂ ψ ∂ Z ) 2 , g c S S = ∇ S c ⋅ ∇ S c = ( ∂ S c ∂ ψ ∂ ψ ∂ R + ∂ S c ∂ θ ∂ θ ∂ R ) 2 + ( ∂ S c ∂ ψ ∂ ψ ∂ Z + ∂ S c ∂ θ ∂ θ ∂ Z ) 2 , g c ψ S = ∇ ψ ⋅ ∇ S c = ∂ S c ∂ ψ [ ( ∂ ψ ∂ R ) 2 + ( ∂ ψ ∂ Z ) 2 ] + ∂ S c ∂ θ [ ∂ ψ ∂ R ∂ θ ∂ R + ∂ ψ ∂ Z ∂ θ ∂ Z ] , g ζ ζ = ∇ ζ ⋅ ∇ ζ = 1 R 2 . J c − 1 = ∇ ψ × ∇ ζ ⋅ ∇ S c = ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ × 1 R ζ ^ ⋅ ∂ S c ∂ ψ ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ + ∂ S c ∂ θ ∂ θ ∂ R R ^ + ∂ θ ∂ Z Z ^ = 1 R ∂ ψ ∂ R ∂ S c ∂ θ ∂ θ ∂ Z − 1 R ∂ ψ ∂ Z ∂ S c ∂ θ ∂ θ ∂ R , g ψ ψ = ∇ ψ ⋅ ∇ ψ = ∂ ψ ∂ R 2 + ∂ ψ ∂ Z 2 , g c S S = ∇ S c ⋅ ∇ S c = ∂ S c ∂ ψ ∂ ψ ∂ R + ∂ S c ∂ θ ∂ θ ∂ R 2 + ∂ S c ∂ ψ ∂ ψ ∂ Z + ∂ S c ∂ θ ∂ θ ∂ Z 2 , g c ψ S = ∇ ψ ⋅ ∇ S c = ∂ S c ∂ ψ ∂ ψ ∂ R 2 + ∂ ψ ∂ Z 2 + ∂ S c ∂ θ ∂ ψ ∂ R ∂ θ ∂ R + ∂ ψ ∂ Z ∂ θ ∂ Z , g ζ ζ = ∇ ζ ⋅ ∇ ζ = 1 R 2 . {:[J_(c)^(-1)=grad psi xx grad zeta*gradS_(c)],[=((del psi)/(del R)( hat(R))+(del psi)/(del Z)( hat(Z)))xx(1)/(R) hat(zeta)*[(delS_(c))/(del psi)((del psi)/(del R)( hat(R))+(del psi)/(del Z)( hat(Z))):}],[{:+(delS_(c))/(del theta)((del theta)/(del R)( hat(R))+(del theta)/(del Z)( hat(Z)))]],[=(1)/(R)(del psi)/(del R)(delS_(c))/(del theta)(del theta)/(del Z)-(1)/(R)(del psi)/(del Z)(delS_(c))/(del theta)(del theta)/(del R)","],[g^(psi psi)=grad psi*grad psi=((del psi)/(del R))^(2)+((del psi)/(del Z))^(2)","],[g_(c)^(SS)=gradS_(c)*gradS_(c)=((delS_(c))/(del psi)(del psi)/(del R)+(delS_(c))/(del theta)(del theta)/(del R))^(2)+((delS_(c))/(del psi)(del psi)/(del Z)+(delS_(c))/(del theta)(del theta)/(del Z))^(2)","],[g_(c)^(psi S)=grad psi*gradS_(c)],[=(delS_(c))/(del psi)[((del psi)/(del R))^(2)+((del psi)/(del Z))^(2)]+(delS_(c))/(del theta)[(del psi)/(del R)(del theta)/(del R)+(del psi)/(del Z)(del theta)/(del Z)]","],[g^(zeta zeta)=grad zeta*grad zeta=(1)/(R^(2)).]:}\begin{aligned} & J_{c}^{-1}= \nabla \psi \times \nabla \zeta \cdot \nabla S_{c} \\ &=\left(\frac{\partial \psi}{\partial R} \hat{R}+\frac{\partial \psi}{\partial Z} \hat{Z}\right) \times \frac{1}{R} \hat{\zeta} \cdot\left[\frac{\partial S_{c}}{\partial \psi}\left(\frac{\partial \psi}{\partial R} \hat{R}+\frac{\partial \psi}{\partial Z} \hat{Z}\right)\right. \\ &\left.+\frac{\partial S_{c}}{\partial \theta}\left(\frac{\partial \theta}{\partial R} \hat{R}+\frac{\partial \theta}{\partial Z} \hat{Z}\right)\right] \\ &= \frac{1}{R} \frac{\partial \psi}{\partial R} \frac{\partial S_{c}}{\partial \theta} \frac{\partial \theta}{\partial Z}-\frac{1}{R} \frac{\partial \psi}{\partial Z} \frac{\partial S_{c}}{\partial \theta} \frac{\partial \theta}{\partial R}, \\ & g^{\psi \psi}=\nabla \psi \cdot \nabla \psi=\left(\frac{\partial \psi}{\partial R}\right)^{2}+\left(\frac{\partial \psi}{\partial Z}\right)^{2}, \\ & g_{c}^{S S}=\nabla S_{c} \cdot \nabla S_{c}=\left(\frac{\partial S_{c}}{\partial \psi} \frac{\partial \psi}{\partial R}+\frac{\partial S_{c}}{\partial \theta} \frac{\partial \theta}{\partial R}\right)^{2}+\left(\frac{\partial S_{c}}{\partial \psi} \frac{\partial \psi}{\partial Z}+\frac{\partial S_{c}}{\partial \theta} \frac{\partial \theta}{\partial Z}\right)^{2}, \\ & g_{c}^{\psi S}= \nabla \psi \cdot \nabla S_{c} \\ &= \frac{\partial S_{c}}{\partial \psi}\left[\left(\frac{\partial \psi}{\partial R}\right)^{2}+\left(\frac{\partial \psi}{\partial Z}\right)^{2}\right]+\frac{\partial S_{c}}{\partial \theta}\left[\frac{\partial \psi}{\partial R} \frac{\partial \theta}{\partial R}+\frac{\partial \psi}{\partial Z} \frac{\partial \theta}{\partial Z}\right], \\ & g^{\zeta \zeta}=\nabla \zeta \cdot \nabla \zeta=\frac{1}{R^{2}} . \end{aligned}Jc−1=∇ψ×∇ζ⋅∇Sc=(∂ψ∂RR^+∂ψ∂ZZ^)×1Rζ^⋅[∂Sc∂ψ(∂ψ∂RR^+∂ψ∂ZZ^)+∂Sc∂θ(∂θ∂RR^+∂θ∂ZZ^)]=1R∂ψ∂R∂Sc∂θ∂θ∂Z−1R∂ψ∂Z∂Sc∂θ∂θ∂R,gψψ=∇ψ⋅∇ψ=(∂ψ∂R)2+(∂ψ∂Z)2,gcSS=∇Sc⋅∇Sc=(∂Sc∂ψ∂ψ∂R+∂Sc∂θ∂θ∂R)2+(∂Sc∂ψ∂ψ∂Z+∂Sc∂θ∂θ∂Z)2,gcψS=∇ψ⋅∇Sc=∂Sc∂ψ[(∂ψ∂R)2+(∂ψ∂Z)2]+∂Sc∂θ[∂ψ∂R∂θ∂R+∂ψ∂Z∂θ∂Z],gζζ=∇ζ⋅∇ζ=1R2.
In the SOL region with ( ψ , ζ , S S ) ψ , ζ , S S (psi,zeta,S_(S))\left(\psi, \zeta, S_{S}\right)(ψ,ζ,SS), considering S S = S S [ ψ ( R S S = S S [ ψ ( R S_(S)=S_(S)[psi(RS_{S}=S_{S}[\psi(RSS=SS[ψ(R, Z ) , Z ] Z ) , Z ] Z),Z]Z), Z]Z),Z], the Laplacian is written as
在 SOL 区域中 ( ψ , ζ , S S ) ψ , ζ , S S (psi,zeta,S_(S))\left(\psi, \zeta, S_{S}\right)(ψ,ζ,SS) ，考虑到 S S = S S [ ψ ( R S S = S S [ ψ ( R S_(S)=S_(S)[psi(RS_{S}=S_{S}[\psi(RSS=SS[ψ(R ， Z ) , Z ] Z ) , Z ] Z),Z]Z), Z]Z),Z] 拉普拉斯被写成
∇ S 2 f = g ψ ψ ∂ 2 f ∂ ψ 2 + 2 g S ψ S ∂ 2 f ∂ ψ ∂ S S + g S S S ∂ 2 f ∂ S S 2 + g ζ ζ ∂ 2 f ∂ ζ 2 + 1 J S [ ∂ ( J S g ψ ψ ) ∂ ψ + ∂ ( J S g S ψ S ) ∂ S S ] ∂ f ∂ ψ + 1 J S [ ∂ ( J S g S ψ S ) ∂ ψ + ∂ ( J S g S S S ) ∂ S S ] ∂ f ∂ S S , ∇ S 2 f = g ψ ψ ∂ 2 f ∂ ψ 2 + 2 g S ψ S ∂ 2 f ∂ ψ ∂ S S + g S S S ∂ 2 f ∂ S S 2 + g ζ ζ ∂ 2 f ∂ ζ 2 + 1 J S ∂ J S g ψ ψ ∂ ψ + ∂ J S g S ψ S ∂ S S ∂ f ∂ ψ + 1 J S ∂ J S g S ψ S ∂ ψ + ∂ J S g S S S ∂ S S ∂ f ∂ S S , {:[grad_(S)^(2)f=g^(psi psi)(del^(2)f)/(delpsi^(2))+2g_(S)^(psi S)(del^(2)f)/(del psi delS_(S))+g_(S)^(SS)(del^(2)f)/(delS_(S)^(2))+g^(zeta zeta)(del^(2)f)/(delzeta^(2))],[+(1)/(J_(S))[(del(J_(S)g^(psi psi)))/(del psi)+(del(J_(S)g_(S)^(psi S)))/(delS_(S))](del f)/(del psi)],[+(1)/(J_(S))[(del(J_(S)g_(S)^(psi S)))/(del psi)+(del(J_(S)g_(S)^(SS)))/(delS_(S))](del f)/(delS_(S))","]:}\begin{aligned} \nabla_{S}^{2} f= & g^{\psi \psi} \frac{\partial^{2} f}{\partial \psi^{2}}+2 g_{S}^{\psi S} \frac{\partial^{2} f}{\partial \psi \partial S_{S}}+g_{S}^{S S} \frac{\partial^{2} f}{\partial S_{S}^{2}}+g^{\zeta \zeta} \frac{\partial^{2} f}{\partial \zeta^{2}} \\ & +\frac{1}{J_{S}}\left[\frac{\partial\left(J_{S} g^{\psi \psi}\right)}{\partial \psi}+\frac{\partial\left(J_{S} g_{S}^{\psi S}\right)}{\partial S_{S}}\right] \frac{\partial f}{\partial \psi} \\ & +\frac{1}{J_{S}}\left[\frac{\partial\left(J_{S} g_{S}^{\psi S}\right)}{\partial \psi}+\frac{\partial\left(J_{S} g_{S}^{S S}\right)}{\partial S_{S}}\right] \frac{\partial f}{\partial S_{S}}, \end{aligned}∇S2f=gψψ∂2f∂ψ2+2gSψS∂2f∂ψ∂SS+gSSS∂2f∂SS2+gζζ∂2f∂ζ2+1JS[∂(JSgψψ)∂ψ+∂(JSgSψS)∂SS]∂f∂ψ+1JS[∂(JSgSψS)∂ψ+∂(JSgSSS)∂SS]∂f∂SS,
where 哪里
J S − 1 = ∇ ψ × ∇ ζ ⋅ ∇ S S = ( ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ ) × 1 R ζ ^ ⋅ [ ∂ S S ∂ ψ ( ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ ) + ∂ S S ∂ Z Z ^ ] = 1 R ∂ ψ ∂ R ∂ S S ∂ Z , g S S S = ∇ S S ⋅ ∇ S S = ( ∂ S S ∂ ψ ∂ ψ ∂ R ) 2 + ( ∂ S S ∂ ψ ∂ ψ ∂ Z + ∂ S S ∂ Z ) 2 , g S ψ S = ∇ ψ ⋅ ∇ S S = ∂ S S ∂ ψ [ ( ∂ ψ ∂ R ) 2 + ( ∂ ψ ∂ Z ) 2 ] + ∂ ψ ∂ Z ∂ S S ∂ Z . J S − 1 = ∇ ψ × ∇ ζ ⋅ ∇ S S = ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ × 1 R ζ ^ ⋅ ∂ S S ∂ ψ ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ + ∂ S S ∂ Z Z ^ = 1 R ∂ ψ ∂ R ∂ S S ∂ Z , g S S S = ∇ S S ⋅ ∇ S S = ∂ S S ∂ ψ ∂ ψ ∂ R 2 + ∂ S S ∂ ψ ∂ ψ ∂ Z + ∂ S S ∂ Z 2 , g S ψ S = ∇ ψ ⋅ ∇ S S = ∂ S S ∂ ψ ∂ ψ ∂ R 2 + ∂ ψ ∂ Z 2 + ∂ ψ ∂ Z ∂ S S ∂ Z . {:[J_(S)^(-1)=grad psi xx grad zeta*gradS_(S)],[=((del psi)/(del R)( hat(R))+(del psi)/(del Z)( hat(Z)))xx(1)/(R) hat(zeta)*[(delS_(S))/(del psi)((del psi)/(del R)( hat(R))+(del psi)/(del Z)( hat(Z)))+(delS_(S))/(del Z)( hat(Z))]],[=(1)/(R)(del psi)/(del R)(delS_(S))/(del Z)","],[g_(S)^(SS)=gradS_(S)*gradS_(S)=((delS_(S))/(del psi)(del psi)/(del R))^(2)+((delS_(S))/(del psi)(del psi)/(del Z)+(delS_(S))/(del Z))^(2)","],[g_(S)^(psi S)=grad psi*gradS_(S)=(delS_(S))/(del psi)[((del psi)/(del R))^(2)+((del psi)/(del Z))^(2)]+(del psi)/(del Z)(delS_(S))/(del Z).]:}\begin{aligned} J_{S}^{-1}= & \nabla \psi \times \nabla \zeta \cdot \nabla S_{S} \\ = & \left(\frac{\partial \psi}{\partial R} \hat{R}+\frac{\partial \psi}{\partial Z} \hat{Z}\right) \times \frac{1}{R} \hat{\zeta} \cdot\left[\frac{\partial S_{S}}{\partial \psi}\left(\frac{\partial \psi}{\partial R} \hat{R}+\frac{\partial \psi}{\partial Z} \hat{Z}\right)+\frac{\partial S_{S}}{\partial Z} \hat{Z}\right] \\ = & \frac{1}{R} \frac{\partial \psi}{\partial R} \frac{\partial S_{S}}{\partial Z}, \\ & g_{S}^{S S}=\nabla S_{S} \cdot \nabla S_{S}=\left(\frac{\partial S_{S}}{\partial \psi} \frac{\partial \psi}{\partial R}\right)^{2}+\left(\frac{\partial S_{S}}{\partial \psi} \frac{\partial \psi}{\partial Z}+\frac{\partial S_{S}}{\partial Z}\right)^{2}, \\ g_{S}^{\psi S} & =\nabla \psi \cdot \nabla S_{S}=\frac{\partial S_{S}}{\partial \psi}\left[\left(\frac{\partial \psi}{\partial R}\right)^{2}+\left(\frac{\partial \psi}{\partial Z}\right)^{2}\right]+\frac{\partial \psi}{\partial Z} \frac{\partial S_{S}}{\partial Z} . \end{aligned}JS−1=∇ψ×∇ζ⋅∇SS=(∂ψ∂RR^+∂ψ∂ZZ^)×1Rζ^⋅[∂SS∂ψ(∂ψ∂RR^+∂ψ∂ZZ^)+∂SS∂ZZ^]=1R∂ψ∂R∂SS∂Z,gSSS=∇SS⋅∇SS=(∂SS∂ψ∂ψ∂R)2+(∂SS∂ψ∂ψ∂Z+∂SS∂Z)2,gSψS=∇ψ⋅∇SS=∂SS∂ψ[(∂ψ∂R)2+(∂ψ∂Z)2]+∂ψ∂Z∂SS∂Z.
For the shared grids between core and SOL regions at the separatrix, the perpendicular Laplacian operator is expanded as
对于分离点处核心区域和 SOL 区域之间的共享网格，垂直拉普拉斯算子扩展为
∇ ⊥ , X 2 f = g ψ ψ ψ ∂ 2 f ∂ ψ 2 | b 0 × ζ + 1 J X ∂ J X g ψ ψ ∂ ψ ∂ f ∂ ψ | b 0 × ζ + g ζ ζ ∂ 2 f ∂ ζ ∇ ⊥ , X 2 f = g ψ ψ ψ ∂ 2 f ∂ ψ 2 b 0 × ζ + 1 J X ∂ J X g ψ ψ ∂ ψ ∂ f ∂ ψ b 0 × ζ + g ζ ζ ∂ 2 f ∂ ζ grad_(_|_,X)^(2)f=g^(psi psi psi)(del^(2)f)/(delpsi^(2))|_(b_(0)xx zeta)+(1)/(J_(X))(delJ_(X)g^(psi psi))/(del psi)(del f)/(del psi)|_(b_(0)xx zeta)+g^(zeta zeta)(del^(2)f)/(del zeta)\nabla_{\perp, X}^{2} f=\left.g^{\psi \psi \psi} \frac{\partial^{2} f}{\partial \psi^{2}}\right|_{b_{0} \times \zeta}+\left.\frac{1}{J_{X}} \frac{\partial J_{X} g^{\psi \psi}}{\partial \psi} \frac{\partial f}{\partial \psi}\right|_{b_{0} \times \zeta}+g^{\zeta \zeta} \frac{\partial^{2} f}{\partial \zeta}∇⊥,X2f=gψψψ∂2f∂ψ2|b0×ζ+1JX∂JXgψψ∂ψ∂f∂ψ|b0×ζ+gζζ∂2f∂ζ
where | b 0 × ζ ζ b 0 × ζ ζ |_(b_(0)xx zeta zeta)\left.\right|_{b_{0} \times \zeta \zeta}|b0×ζζ represents the partial derivative with respect to ψ ψ psi\psiψ among the orthogonal grids along the b 0 × ζ b 0 × ζ b_(0)xx zeta\boldsymbol{b}_{0} \times \zetab0×ζ direction, and
其中 | b 0 × ζ ζ b 0 × ζ ζ |_(b_(0)xx zeta zeta)\left.\right|_{b_{0} \times \zeta \zeta}|b0×ζζ 表示相对于沿 b 0 × ζ b 0 × ζ b_(0)xx zeta\boldsymbol{b}_{0} \times \zetab0×ζ 方向的正交网格 ψ ψ psi\psiψ 之间的偏导数，以及
J X − 1 = ∇ ψ × ∇ ζ ⋅ ∇ S X = ( ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ ) × 1 R ζ ^ ⋅ [ ∂ S X ∂ R R ^ + ∂ S X ∂ Z Z ^ ] = 1 R ∂ ψ ∂ R ∂ S X ∂ Z − 1 R ∂ ψ ∂ Z ∂ S X ∂ R . J X − 1 = ∇ ψ × ∇ ζ ⋅ ∇ S X = ∂ ψ ∂ R R ^ + ∂ ψ ∂ Z Z ^ × 1 R ζ ^ ⋅ ∂ S X ∂ R R ^ + ∂ S X ∂ Z Z ^ = 1 R ∂ ψ ∂ R ∂ S X ∂ Z − 1 R ∂ ψ ∂ Z ∂ S X ∂ R . {:[J_(X)^(-1)=grad psi xx grad zeta*gradS_(X)],[=((del psi)/(del R)( hat(R))+(del psi)/(del Z)( hat(Z)))xx(1)/(R) hat(zeta)*[(delS_(X))/(del R)( hat(R))+(delS_(X))/(del Z)( hat(Z))]],[=(1)/(R)(del psi)/(del R)(delS_(X))/(del Z)-(1)/(R)(del psi)/(del Z)(delS_(X))/(del R).]:}\begin{aligned} J_{X}^{-1} & =\nabla \psi \times \nabla \zeta \cdot \nabla S_{X} \\ & =\left(\frac{\partial \psi}{\partial R} \hat{R}+\frac{\partial \psi}{\partial Z} \hat{Z}\right) \times \frac{1}{R} \hat{\zeta} \cdot\left[\frac{\partial S_{X}}{\partial R} \hat{R}+\frac{\partial S_{X}}{\partial Z} \hat{Z}\right] \\ & =\frac{1}{R} \frac{\partial \psi}{\partial R} \frac{\partial S_{X}}{\partial Z}-\frac{1}{R} \frac{\partial \psi}{\partial Z} \frac{\partial S_{X}}{\partial R} . \end{aligned}JX−1=∇ψ×∇ζ⋅∇SX=(∂ψ∂RR^+∂ψ∂ZZ^)×1Rζ^⋅[∂SX∂RR^+∂SX∂ZZ^]=1R∂ψ∂R∂SX∂Z−1R∂ψ∂Z∂SX∂R.
Taking advantage of the toroidal symmetry of FRC, we could transform the Laplacian operator into Fourier space with respect to the toroidal angle ζ ζ zeta\zetaζ, which can avoid solving the 3 dimensional matrix. In this paper, we can also simplify Eqs. (12)-(14) assuming k ⊥ ≫ k k ⊥ ≫ k k_(_|_)≫kk_{\perp} \gg kk⊥≫k and k ∥ ≪ k ⊥ k ∥ ≪ k ⊥ k_(||)≪k_(_|_)k_{\|} \ll k_{\perp}k∥≪k⊥, where κ = ∇ B 0 / B 0 κ = ∇ B 0 / B 0 kappa=gradB_(0)//B_(0)\kappa=\nabla B_{0} / B_{0}κ=∇B0/B0. Thus, for each toroidal mode n n nnn ( ∂ / ∂ ζ = i n ζ ) ( ∂ / ∂ ζ = i n ζ ) (del//del zeta=in zeta)(\partial / \partial \zeta=i n \zeta)(∂/∂ζ=inζ), the Laplacian operators in the core region, SOL region, and at the separatrix can be written as
利用FRC的环形对称性，我们可以将拉普拉斯算子变换为相对于环形角 ζ ζ zeta\zetaζ 的傅里叶空间，从而避免求解三维矩阵。在本文中，我们还可以简化方程。（12）-（14） 假设 k ⊥ ≫ k k ⊥ ≫ k k_(_|_)≫kk_{\perp} \gg kk⊥≫k 和 k ∥ ≪ k ⊥ k ∥ ≪ k ⊥ k_(||)≪k_(_|_)k_{\|} \ll k_{\perp}k∥≪k⊥ ，其中 κ = ∇ B 0 / B 0 κ = ∇ B 0 / B 0 kappa=gradB_(0)//B_(0)\kappa=\nabla B_{0} / B_{0}κ=∇B0/B0 .因此，对于每个环形模 n n nnn ( ∂ / ∂ ζ = i n ζ ) ( ∂ / ∂ ζ = i n ζ ) (del//del zeta=in zeta)(\partial / \partial \zeta=i n \zeta)(∂/∂ζ=inζ) 态，核心区域、SOL 区域和分离点处的拉普拉斯算子可以写成
∇ ⊥ c 2 f n ≈ g ψ ψ ∂ 2 f n ∂ ψ 2 + 2 g c ψ S ∂ 2 f n ∂ ψ ∂ S c + g c S S ∂ 2 f n ∂ S c 2 − n 2 g ζ ζ f n , ∇ ⊥ S 2 f n ≈ g ψ ψ ∂ 2 f n ∂ ψ 2 + 2 g S ψ S ∂ 2 f n ∂ ψ ∂ S S + g S S S ∂ 2 f n ∂ S S 2 − n 2 g ζ ζ f n , ∇ ⊥ X f n 2 ≈ g ψ ψ ψ ∂ 2 f n ∂ ψ 2 | b 0 × ζ − n 2 g ζ ζ f n , ∇ ⊥ c 2 f n ≈ g ψ ψ ∂ 2 f n ∂ ψ 2 + 2 g c ψ S ∂ 2 f n ∂ ψ ∂ S c + g c S S ∂ 2 f n ∂ S c 2 − n 2 g ζ ζ f n , ∇ ⊥ S 2 f n ≈ g ψ ψ ∂ 2 f n ∂ ψ 2 + 2 g S ψ S ∂ 2 f n ∂ ψ ∂ S S + g S S S ∂ 2 f n ∂ S S 2 − n 2 g ζ ζ f n , ∇ ⊥ X f n 2 ≈ g ψ ψ ψ ∂ 2 f n ∂ ψ 2 b 0 × ζ − n 2 g ζ ζ f n , {:[grad_(_|_ c)^(2)f_(n)~~g^(psi psi)(del^(2)f_(n))/(delpsi^(2))+2g_(c)^(psi S)(del^(2)f_(n))/(del psi delS_(c))+g_(c)^(SS)(del^(2)f_(n))/(delS_(c)^(2))-n^(2)g^(zeta zeta)f_(n)","],[grad_(_|_ S)^(2)f_(n)~~g^(psi psi)(del^(2)f_(n))/(delpsi^(2))+2g_(S)^(psi S)(del^(2)f_(n))/(del psi delS_(S))+g_(S)^(SS)(del^(2)f_(n))/(delS_(S)^(2))-n^(2)g^(zeta zeta)f_(n)","],[grad_(_|_ Xf_(n))^(2)~~g^(psi psi psi)(del^(2)f_(n))/(delpsi^(2))|_(b_(0)xx zeta)-n^(2)g^(zeta zeta)f_(n)","]:}\begin{gathered} \nabla_{\perp c}^{2} f_{n} \approx g^{\psi \psi} \frac{\partial^{2} f_{n}}{\partial \psi^{2}}+2 g_{c}^{\psi S} \frac{\partial^{2} f_{n}}{\partial \psi \partial S_{c}}+g_{c}^{S S} \frac{\partial^{2} f_{n}}{\partial S_{c}^{2}}-n^{2} g^{\zeta \zeta} f_{n}, \\ \nabla_{\perp S}^{2} f_{n} \approx g^{\psi \psi} \frac{\partial^{2} f_{n}}{\partial \psi^{2}}+2 g_{S}^{\psi S} \frac{\partial^{2} f_{n}}{\partial \psi \partial S_{S}}+g_{S}^{S S} \frac{\partial^{2} f_{n}}{\partial S_{S}^{2}}-n^{2} g^{\zeta \zeta} f_{n}, \\ \left.\nabla_{\perp X f_{n}}^{2} \approx g^{\psi \psi \psi} \frac{\partial^{2} f_{n}}{\partial \psi^{2}}\right|_{b_{0} \times \zeta}-n^{2} g^{\zeta \zeta} f_{n}, \end{gathered}∇⊥c2fn≈gψψ∂2fn∂ψ2+2gcψS∂2fn∂ψ∂Sc+gcSS∂2fn∂Sc2−n2gζζfn,∇⊥S2fn≈gψψ∂2fn∂ψ2+2gSψS∂2fn∂ψ∂SS+gSSS∂2fn∂SS2−n2gζζfn,∇⊥Xfn2≈gψψψ∂2fn∂ψ2|b0×ζ−n2gζζfn,
where f n f n f_(n)f_{n}fn is the n n nnn toroidal mode component of f f fff.
其中 f n f n f_(n)f_{n}fn 是 的 n n nnn f f fff 环形模态分量。
The grids shown in Figs. 2 and 3 are only for illustrating the field aligned grid algorithm in cylindrical coordinates clearly, which are much sparser than the realistic simulation.
图 2 和图 3 中所示的网格仅用于清楚地说明圆柱坐标中的场对齐网格算法，比实际模拟要稀疏得多。
FIG. 3. The green + symbols denote the orthogonal grids in the SOL and core regions with respect to the shared grids at the separatrix, and the dotted line in magenta is along the perpendicular direction b 0 × ζ ^ b 0 × ζ ^ b_(0)xx hat(zeta)\mathbf{b}_{0} \times \hat{\zeta}b0×ζ^ from the shared grids as shown by the black stars. The grids shown here are only for illustrating the algorithm, which are much sparser than the ones used in realistic simulation.
图 3.绿色 + 符号表示 SOL 和核心区域中相对于分离点处共享网格的正交网格，洋红色的虚线沿共享网格的垂直方向 b 0 × ζ ^ b 0 × ζ ^ b_(0)xx hat(zeta)\mathbf{b}_{0} \times \hat{\zeta}b0×ζ^ ，如黑色星形所示。此处显示的网格仅用于说明算法，这些网格比实际模拟中使用的网格要稀疏得多。

9. PHYSICS MODEL 9. 物理模型

10. A. Formulation 10. A. 配方

We use the electrostatic Vlasov-Poisson system for physics simulation in this paper. Particle dynamics is described by the gyrokinetic equation using gyrocenter position R R R\boldsymbol{R}R, magnetic moment μ μ mu\muμ, and parallel velocity v ∥ v ∥ v_(||)v_{\|}v∥as independent variables in five dimensional phase space
本文使用静电Vlasov-Poisson系统进行物理模拟。粒子动力学由陀螺动力学方程描述，使用陀螺中心位置 R R R\boldsymbol{R}R 、磁矩 μ μ mu\muμ 和平行速度 v ∥ v ∥ v_(||)v_{\|}v∥ 作为五维相空间中的自变量
( ∂ ∂ t + R ˙ ⋅ ∇ + v ˙ ∥ ∂ ∂ v ∥ ) f α ( R , v ∥ , μ , t ) = 0 , R ˙ = v ∥ b + v E + v d , v ∥ ˙ = − 1 m α B ∗ B ∥ ∗ ⋅ ( Z α ∇ ⟨ δ ϕ ⟩ + μ ∇ B ) , ∂ ∂ t + R ˙ ⋅ ∇ + v ˙ ∥ ∂ ∂ v ∥ f α R , v ∥ , μ , t = 0 , R ˙ = v ∥ b + v E + v d , v ∥ ˙ = − 1 m α B ∗ B ∥ ∗ ⋅ Z α ∇ ⟨ δ ϕ ⟩ + μ ∇ B , {:[((del)/(del t)+(R^(˙))*grad+v^(˙)_(||)(del)/(delv_(||)))f_(alpha)(R,v_(||),mu,t)=0","],[R^(˙)=v_(||)b+v_(E)+v_(d)","],[v_(||)^(˙)=-(1)/(m_(alpha))(B^(**))/(B_(||)^(**))*(Z_(alpha)grad(:delta phi:)+mu grad B)","]:}\begin{gathered} \left(\frac{\partial}{\partial t}+\dot{\boldsymbol{R}} \cdot \nabla+\dot{v}_{\|} \frac{\partial}{\partial v_{\|}}\right) f_{\alpha}\left(\boldsymbol{R}, v_{\|}, \mu, t\right)=0, \\ \dot{\boldsymbol{R}}=v_{\|} \boldsymbol{b}+\boldsymbol{v}_{\boldsymbol{E}}+\boldsymbol{v}_{\boldsymbol{d}}, \\ \dot{v_{\|}}=-\frac{1}{m_{\alpha}} \frac{\boldsymbol{B}^{*}}{B_{\|}^{*}} \cdot\left(Z_{\alpha} \nabla\langle\delta \phi\rangle+\mu \nabla B\right), \end{gathered}(∂∂t+R˙⋅∇+v˙∥∂∂v∥)fα(R,v∥,μ,t)=0,R˙=v∥b+vE+vd,v∥˙=−1mαB∗B∥∗⋅(Zα∇⟨δϕ⟩+μ∇B),
where Z α , m α Z α , m α Z_(alpha),m_(alpha)Z_{\alpha}, m_{\alpha}Zα,mα, and f α f α f_(alpha)f_{\alpha}fα are the charge, mass, and distribution function of α α alpha\alphaα species. B B B\boldsymbol{B}B is the equilibrium magnetic field, b = B / B , B ∗ = B b = B / B , B ∗ = B b=B//B,B^(**)=B\boldsymbol{b}=\boldsymbol{B} / B, \boldsymbol{B}^{*}=\boldsymbol{B}b=B/B,B∗=B + ( B v ∥ / Ω c α ) ∇ × b + B v ∥ / Ω c α ∇ × b +(Bv_(||)//Omega_(c alpha))grad xx b+\left(B v_{\|} / \Omega_{c \alpha}\right) \nabla \times \boldsymbol{b}+(Bv∥/Ωcα)∇×b, and B ∥ ∗ = b ⋅ B ∗ . δ ϕ B ∥ ∗ = b ⋅ B ∗ . δ ϕ B_(||)^(**)=b*B^(**).delta phiB_{\|}^{*}=\boldsymbol{b} \cdot \boldsymbol{B}^{*} . \delta \phiB∥∗=b⋅B∗.δϕ is the electrostatic potential. ⟨ ⋯ ⟩ = ( 1 / 2 π ) ∫ d x d ξ ( ⋯ ) δ ( R + ρ α − x ) ⟨ ⋯ ⟩ = ( 1 / 2 π ) ∫ d x d ξ ( ⋯ ) δ R + ρ α − x (:cdots:)=(1//2pi)int dxd xi(cdots)delta(R+rho_(alpha)-x)\langle\cdots\rangle=(1 / 2 \pi) \int \boldsymbol{d} x d \xi(\cdots) \delta\left(\boldsymbol{R}+\boldsymbol{\rho}_{\boldsymbol{\alpha}}-\boldsymbol{x}\right)⟨⋯⟩=(1/2π)∫dxdξ(⋯)δ(R+ρα−x) represents the gyrophase average, ξ ξ xi\xiξ is the gyrophase angle, x x x\boldsymbol{x}x represents the particle position, ρ α = b × v ⊥ / Ω c α ρ α = b × v ⊥ / Ω c α rho_(alpha)=b xxv_(_|_)//Omega_(c alpha)\boldsymbol{\rho}_{\boldsymbol{\alpha}}=\boldsymbol{b} \times \boldsymbol{v}_{\perp} / \Omega_{c \alpha}ρα=b×v⊥/Ωcα is the particle gyroradius, and Ω c α Ω c α Omega_(c alpha)\Omega_{c \alpha}Ωcα is the particle cyclotron frequency. v E v E v_(E)\boldsymbol{v}_{\boldsymbol{E}}vE is the E × B E × B E xx B\boldsymbol{E} \times \boldsymbol{B}E×B velocity, and v d v d v_(d)\boldsymbol{v}_{\boldsymbol{d}}vd is the magnetic drift velocity, which are given as
其中 Z α , m α Z α , m α Z_(alpha),m_(alpha)Z_{\alpha}, m_{\alpha}Zα,mα ， 和 f α f α f_(alpha)f_{\alpha}fα 是物种的 α α alpha\alphaα 电荷、质量和分布函数。 B B B\boldsymbol{B}B 是平衡磁场， b = B / B , B ∗ = B b = B / B , B ∗ = B b=B//B,B^(**)=B\boldsymbol{b}=\boldsymbol{B} / B, \boldsymbol{B}^{*}=\boldsymbol{B}b=B/B,B∗=B + ( B v ∥ / Ω c α ) ∇ × b + B v ∥ / Ω c α ∇ × b +(Bv_(||)//Omega_(c alpha))grad xx b+\left(B v_{\|} / \Omega_{c \alpha}\right) \nabla \times \boldsymbol{b}+(Bv∥/Ωcα)∇×b B ∥ ∗ = b ⋅ B ∗ . δ ϕ B ∥ ∗ = b ⋅ B ∗ . δ ϕ B_(||)^(**)=b*B^(**).delta phiB_{\|}^{*}=\boldsymbol{b} \cdot \boldsymbol{B}^{*} . \delta \phiB∥∗=b⋅B∗.δϕ 是静电势。 ⟨ ⋯ ⟩ = ( 1 / 2 π ) ∫ d x d ξ ( ⋯ ) δ ( R + ρ α − x ) ⟨ ⋯ ⟩ = ( 1 / 2 π ) ∫ d x d ξ ( ⋯ ) δ R + ρ α − x (:cdots:)=(1//2pi)int dxd xi(cdots)delta(R+rho_(alpha)-x)\langle\cdots\rangle=(1 / 2 \pi) \int \boldsymbol{d} x d \xi(\cdots) \delta\left(\boldsymbol{R}+\boldsymbol{\rho}_{\boldsymbol{\alpha}}-\boldsymbol{x}\right)⟨⋯⟩=(1/2π)∫dxdξ(⋯)δ(R+ρα−x) 表示陀螺相位平均值， ξ ξ xi\xiξ 表示陀螺相位角， x x x\boldsymbol{x}x 表示粒子位置， ρ α = b × v ⊥ / Ω c α ρ α = b × v ⊥ / Ω c α rho_(alpha)=b xxv_(_|_)//Omega_(c alpha)\boldsymbol{\rho}_{\boldsymbol{\alpha}}=\boldsymbol{b} \times \boldsymbol{v}_{\perp} / \Omega_{c \alpha}ρα=b×v⊥/Ωcα 表示粒子回转半径， Ω c α Ω c α Omega_(c alpha)\Omega_{c \alpha}Ωcα 表示粒子回旋加速器频率。 v E v E v_(E)\boldsymbol{v}_{\boldsymbol{E}}vE 是 E × B E × B E xx B\boldsymbol{E} \times \boldsymbol{B}E×B 速度， v d v d v_(d)\boldsymbol{v}_{\boldsymbol{d}}vd 是磁漂移速度，其表示为
v E = c b × ∇ ⟨ δ ϕ ⟩ B ∥ ∗ , v d = c m α v ∥ 2 Z α B ∥ ∗ b × ( b ⋅ ∇ b ) + c μ Z α B ∥ ∗ b × ∇ B . v E = c b × ∇ ⟨ δ ϕ ⟩ B ∥ ∗ , v d = c m α v ∥ 2 Z α B ∥ ∗ b × ( b ⋅ ∇ b ) + c μ Z α B ∥ ∗ b × ∇ B . {:[v_(E)=(cb xx grad(:delta phi:))/(B_(||)^(**))","],[v_(d)=(cm_(alpha)v_(||)^(2))/(Z_(alpha)B_(||)^(**))b xx(b*grad b)+(c mu)/(Z_(alpha)B_(||)^(**))b xx grad B.]:}\begin{gathered} \boldsymbol{v}_{\boldsymbol{E}}=\frac{c \boldsymbol{b} \times \nabla\langle\delta \phi\rangle}{B_{\|}^{*}}, \\ \boldsymbol{v}_{\boldsymbol{d}}=\frac{c m_{\alpha} v_{\|}^{2}}{Z_{\alpha} B_{\|}^{*}} \boldsymbol{b} \times(\boldsymbol{b} \cdot \nabla \boldsymbol{b})+\frac{c \mu}{Z_{\alpha} B_{\|}^{*}} \boldsymbol{b} \times \nabla B . \end{gathered}vE=cb×∇⟨δϕ⟩B∥∗,vd=cmαv∥2ZαB∥∗b×(b⋅∇b)+cμZαB∥∗b×∇B.
In order to minimize the particle noise, the perturbative δ f δ f delta f\delta fδf simulation method 16 , 17 16 , 17 ^(16,17){ }^{16,17}16,17 is applied. The particle distribution is decomposed into equilibrium and perturbed parts as f α = f α 0 ( R , μ , v ∥ ) + δ f α ( R , μ , v ∥ , t ) f α = f α 0 R , μ , v ∥ + δ f α R , μ , v ∥ , t f_(alpha)=f_(alpha0)(R,mu,v_(||))+deltaf_(alpha)(R,mu,v_(||),t)f_{\alpha}=f_{\alpha 0}\left(\boldsymbol{R}, \mu, v_{\|}\right)+\delta f_{\alpha}\left(\boldsymbol{R}, \mu, v_{\|}, t\right)fα=fα0(R,μ,v∥)+δfα(R,μ,v∥,t), and the equilibrium f α 0 f α 0 f_(alpha0)f_{\alpha 0}fα0 satisfies the following equation:
为了最小化粒子噪声，采用扰动 δ f δ f delta f\delta fδf 模拟方法 16 , 17 16 , 17 ^(16,17){ }^{16,17}16,17 。粒子分布分解为平衡部分和扰动部分 f α = f α 0 ( R , μ , v ∥ ) + δ f α ( R , μ , v ∥ , t ) f α = f α 0 R , μ , v ∥ + δ f α R , μ , v ∥ , t f_(alpha)=f_(alpha0)(R,mu,v_(||))+deltaf_(alpha)(R,mu,v_(||),t)f_{\alpha}=f_{\alpha 0}\left(\boldsymbol{R}, \mu, v_{\|}\right)+\delta f_{\alpha}\left(\boldsymbol{R}, \mu, v_{\|}, t\right)fα=fα0(R,μ,v∥)+δfα(R,μ,v∥,t) ，均衡 f α 0 f α 0 f_(alpha0)f_{\alpha 0}fα0 满足以下方程：
L 0 f α 0 = 0 L 0 f α 0 = 0 L_(0)f_(alpha0)=0L_{0} f_{\alpha 0}=0L0fα0=0
where L 0 = ∂ / ∂ t + ( v ∥ b + v d ) ⋅ ∇ − ( μ / m α ) B ∗ ⋅ ∇ B / B ∥ ∗ ( ∂ / ∂ v ∥ ) L 0 = ∂ / ∂ t + v ∥ b + v d ⋅ ∇ − μ / m α B ∗ ⋅ ∇ B / B ∥ ∗ ∂ / ∂ v ∥ L_(0)=del//del t+(v_(||)b+v_(d))*grad-(mu//m_(alpha))B^(**)*grad B//B_(||)^(**)(del//delv_(||))L_{0}=\partial / \partial t+\left(v_{\|} \boldsymbol{b}+\boldsymbol{v}_{\boldsymbol{d}}\right) \cdot \nabla-\left(\mu / m_{\alpha}\right) \boldsymbol{B}^{*} \cdot \nabla B / B_{\|}^{*}\left(\partial / \partial v_{\|}\right)L0=∂/∂t+(v∥b+vd)⋅∇−(μ/mα)B∗⋅∇B/B∥∗(∂/∂v∥). Because v d v d v_(d)\boldsymbol{v}_{\boldsymbol{d}}vd is only in the ζ ^ ζ ^ hat(zeta)\hat{\zeta}ζ^ direction, the particle drift orbit width is zero in FRC geometry, then f α 0 = n α 0 ( m α 2 π T α 0 ) 1.5 exp ⁡ [ − m α v ∥ 2 + 2 μ B 2 T α 0 ] f α 0 = n α 0 m α 2 π T α 0 1.5 exp ⁡ − m α v ∥ 2 + 2 μ B 2 T α 0 f_(alpha0)=n_(alpha0)((m_(alpha))/(2piT_(alpha0)))^(1.5)exp[-(m_(alpha)v_(||)^(2)+2mu B)/(2T_(alpha0))]f_{\alpha 0}=n_{\alpha 0}\left(\frac{m_{\alpha}}{2 \pi T_{\alpha 0}}\right)^{1.5} \exp \left[-\frac{m_{\alpha} v_{\|}^{2}+2 \mu B}{2 T_{\alpha 0}}\right]fα0=nα0(mα2πTα0)1.5exp⁡[−mαv∥2+2μB2Tα0] is the exact solution of Eq. (21), and n α 0 ( ψ ) n α 0 ( ψ ) n_(alpha0)(psi)n_{\alpha 0}(\psi)nα0(ψ) and T α 0 ( ψ ) T α 0 ( ψ ) T_(alpha0)(psi)T_{\alpha 0}(\psi)Tα0(ψ) are the 1 D 1 D 1D1 \mathrm{D}1D function of magnetic flux surface. Subtracting Eq. (18) by Eq. (21), we have the equation for perturbed distribution δ f α δ f α deltaf_(alpha)\delta f_{\alpha}δfα
其中 L 0 = ∂ / ∂ t + ( v ∥ b + v d ) ⋅ ∇ − ( μ / m α ) B ∗ ⋅ ∇ B / B ∥ ∗ ( ∂ / ∂ v ∥ ) L 0 = ∂ / ∂ t + v ∥ b + v d ⋅ ∇ − μ / m α B ∗ ⋅ ∇ B / B ∥ ∗ ∂ / ∂ v ∥ L_(0)=del//del t+(v_(||)b+v_(d))*grad-(mu//m_(alpha))B^(**)*grad B//B_(||)^(**)(del//delv_(||))L_{0}=\partial / \partial t+\left(v_{\|} \boldsymbol{b}+\boldsymbol{v}_{\boldsymbol{d}}\right) \cdot \nabla-\left(\mu / m_{\alpha}\right) \boldsymbol{B}^{*} \cdot \nabla B / B_{\|}^{*}\left(\partial / \partial v_{\|}\right)L0=∂/∂t+(v∥b+vd)⋅∇−(μ/mα)B∗⋅∇B/B∥∗(∂/∂v∥) .因为 v d v d v_(d)\boldsymbol{v}_{\boldsymbol{d}}vd 只是在方向上 ζ ^ ζ ^ hat(zeta)\hat{\zeta}ζ^ ，粒子漂移轨道宽度在FRC几何中为零，所以 f α 0 = n α 0 ( m α 2 π T α 0 ) 1.5 exp ⁡ [ − m α v ∥ 2 + 2 μ B 2 T α 0 ] f α 0 = n α 0 m α 2 π T α 0 1.5 exp ⁡ − m α v ∥ 2 + 2 μ B 2 T α 0 f_(alpha0)=n_(alpha0)((m_(alpha))/(2piT_(alpha0)))^(1.5)exp[-(m_(alpha)v_(||)^(2)+2mu B)/(2T_(alpha0))]f_{\alpha 0}=n_{\alpha 0}\left(\frac{m_{\alpha}}{2 \pi T_{\alpha 0}}\right)^{1.5} \exp \left[-\frac{m_{\alpha} v_{\|}^{2}+2 \mu B}{2 T_{\alpha 0}}\right]fα0=nα0(mα2πTα0)1.5exp⁡[−mαv∥2+2μB2Tα0] 是方程（21）的精确解，和 n α 0 ( ψ ) n α 0 ( ψ ) n_(alpha0)(psi)n_{\alpha 0}(\psi)nα0(ψ) 是 T α 0 ( ψ ) T α 0 ( ψ ) T_(alpha0)(psi)T_{\alpha 0}(\psi)Tα0(ψ) 磁通量面的 1 D 1 D 1D1 \mathrm{D}1D 函数。用方程（21）减去方程（18），我们得到扰动分布 δ f α δ f α deltaf_(alpha)\delta f_{\alpha}δfα 方程
L δ f α = − δ L f α 0 , L δ f α = − δ L f α 0 , L deltaf_(alpha)=-delta Lf_(alpha0),L \delta f_{\alpha}=-\delta L f_{\alpha 0},Lδfα=−δLfα0,
where δ L = v E ⋅ ∇ − ( Z α / m α / B ∥ ∗ ) B ∗ ⋅ ∇ ⟨ δ ϕ ⟩ ( ∂ / ∂ v ∥ ) δ L = v E ⋅ ∇ − Z α / m α / B ∥ ∗ B ∗ ⋅ ∇ ⟨ δ ϕ ⟩ ∂ / ∂ v ∥ delta L=v_(E)*grad-(Z_(alpha)//m_(alpha)//B_(||)^(**))B^(**)*grad(:delta phi:)(del//delv_(||))\delta L=\boldsymbol{v}_{\boldsymbol{E}} \cdot \nabla-\left(Z_{\alpha} / m_{\alpha} / B_{\|}^{*}\right) \boldsymbol{B}^{*} \cdot \nabla\langle\delta \phi\rangle\left(\partial / \partial v_{\|}\right)δL=vE⋅∇−(Zα/mα/B∥∗)B∗⋅∇⟨δϕ⟩(∂/∂v∥)and L = L 0 L = L 0 L=L_(0)L=L_{0}L=L0 + δ L + δ L +delta L+\delta L+δL. Defining particle weight as w α = δ f α / f α w α = δ f α / f α w_(alpha)=deltaf_(alpha)//f_(alpha)w_{\alpha}=\delta f_{\alpha} / f_{\alpha}wα=δfα/fα, we can derive the weight equation from Eq. (22) as
其中 δ L = v E ⋅ ∇ − ( Z α / m α / B ∥ ∗ ) B ∗ ⋅ ∇ ⟨ δ ϕ ⟩ ( ∂ / ∂ v ∥ ) δ L = v E ⋅ ∇ − Z α / m α / B ∥ ∗ B ∗ ⋅ ∇ ⟨ δ ϕ ⟩ ∂ / ∂ v ∥ delta L=v_(E)*grad-(Z_(alpha)//m_(alpha)//B_(||)^(**))B^(**)*grad(:delta phi:)(del//delv_(||))\delta L=\boldsymbol{v}_{\boldsymbol{E}} \cdot \nabla-\left(Z_{\alpha} / m_{\alpha} / B_{\|}^{*}\right) \boldsymbol{B}^{*} \cdot \nabla\langle\delta \phi\rangle\left(\partial / \partial v_{\|}\right)δL=vE⋅∇−(Zα/mα/B∥∗)B∗⋅∇⟨δϕ⟩(∂/∂v∥) 和 L = L 0 L = L 0 L=L_(0)L=L_{0}L=L0 + δ L + δ L +delta L+\delta L+δL .将粒子重量定义为 w α = δ f α / f α w α = δ f α / f α w_(alpha)=deltaf_(alpha)//f_(alpha)w_{\alpha}=\delta f_{\alpha} / f_{\alpha}wα=δfα/fα ，我们可以从方程（22）中推导出重量方程为
d w α d t = ( 1 − w α ) [ − v E ⋅ ∇ f α 0 f α 0 | v ⊥ − Z α v ∥ b ⋅ ∇ ⟨ δ ϕ ⟩ T α 0 − Z α T α 0 ( μ b × ∇ B 0 m α Ω c α + v ∥ 2 Ω c α b × ( b ⋅ ∇ b ) ) ⋅ ∇ ⟨ δ ϕ ⟩ ] , d w α d t = 1 − w α − v E ⋅ ∇ f α 0 f α 0 v ⊥ − Z α v ∥ b ⋅ ∇ ⟨ δ ϕ ⟩ T α 0 − Z α T α 0 μ b × ∇ B 0 m α Ω c α + v ∥ 2 Ω c α b × ( b ⋅ ∇ b ) ⋅ ∇ ⟨ δ ϕ ⟩ , {:[(dw_(alpha))/(dt)=(1-w_(alpha))[-v_(E)*(gradf_(alpha0))/(f_(alpha0))|_(v_(_|_))-(Z_(alpha)v_(||)b*grad(:delta phi:))/(T_(alpha0)):}],[{:-(Z_(alpha))/(T_(alpha0))((mu b xx gradB_(0))/(m_(alpha)Omega_(c alpha))+(v_(||)^(2))/(Omega_(c alpha))b xx(b*grad b))*grad(:delta phi:)]","]:}\begin{aligned} \frac{d w_{\alpha}}{d t}= & \left(1-w_{\alpha}\right)\left[-\left.\boldsymbol{v}_{\boldsymbol{E}} \cdot \frac{\nabla f_{\alpha 0}}{f_{\alpha 0}}\right|_{v_{\perp}}-\frac{Z_{\alpha} v_{\|} \boldsymbol{b} \cdot \nabla\langle\delta \phi\rangle}{T_{\alpha 0}}\right. \\ & \left.-\frac{Z_{\alpha}}{T_{\alpha 0}}\left(\frac{\mu \boldsymbol{b} \times \nabla B_{0}}{m_{\alpha} \Omega_{c \alpha}}+\frac{v_{\|}^{2}}{\Omega_{c \alpha}} \boldsymbol{b} \times(\boldsymbol{b} \cdot \nabla \boldsymbol{b})\right) \cdot \nabla\langle\delta \phi\rangle\right], \end{aligned}dwαdt=(1−wα)[−vE⋅∇fα0fα0|v⊥−Zαv∥b⋅∇⟨δϕ⟩Tα0−ZαTα0(μb×∇B0mαΩcα+v∥2Ωcαb×(b⋅∇b))⋅∇⟨δϕ⟩],
where we have used the chain rule ∇ f α 0 | v ⊥ = ∇ f α 0 | μ + ( μ f α 0 ∇ B / T α 0 ) ∇ f α 0 v ⊥ = ∇ f α 0 μ + μ f α 0 ∇ B / T α 0 gradf_(alpha0)|_(v_(_|_))= gradf_(alpha0)|_(mu)+(muf_(alpha0)grad B//T_(alpha0))\left.\nabla f_{\alpha 0}\right|_{v_{\perp}}=\left.\nabla f_{\alpha 0}\right|_{\mu}+\left(\mu f_{\alpha 0} \nabla B / T_{\alpha 0}\right)∇fα0|v⊥=∇fα0|μ+(μfα0∇B/Tα0) in the derivation of Eq. (23) from Eq. (22).
其中，我们在从方程（22）推导方程（23）时使用了链式法则 ∇ f α 0 | v ⊥ = ∇ f α 0 | μ + ( μ f α 0 ∇ B / T α 0 ) ∇ f α 0 v ⊥ = ∇ f α 0 μ + μ f α 0 ∇ B / T α 0 gradf_(alpha0)|_(v_(_|_))= gradf_(alpha0)|_(mu)+(muf_(alpha0)grad B//T_(alpha0))\left.\nabla f_{\alpha 0}\right|_{v_{\perp}}=\left.\nabla f_{\alpha 0}\right|_{\mu}+\left(\mu f_{\alpha 0} \nabla B / T_{\alpha 0}\right)∇fα0|v⊥=∇fα0|μ+(μfα0∇B/Tα0) 。
The gyrokinetic Vlasov equation is used for ion species, and its perturbed density is
陀螺动力学弗拉索夫方程用于离子种类，其扰动密度为
⟨ δ n i ( x , t ) ⟩ = ∫ δ f i ( R , μ , v ∥ , t ) d v d R d ξ δ ( R + ρ i − x ) / ( 2 π ) δ n i ( x , t ) = ∫ δ f i R , μ , v ∥ , t d v d R d ξ δ R + ρ i − x / ( 2 π ) (:deltan_(i)(x,t):)=int deltaf_(i)(R,mu,v_(||),t)dvdRd xi delta(R+rho_(i)-x)//(2pi)\left\langle\delta n_{i}(\boldsymbol{x}, t)\right\rangle=\int \delta f_{i}\left(\boldsymbol{R}, \mu, v_{\|}, t\right) \boldsymbol{d} \boldsymbol{v} \boldsymbol{d} \boldsymbol{R} d \xi \delta\left(\boldsymbol{R}+\boldsymbol{\rho}_{\boldsymbol{i}}-\boldsymbol{x}\right) /(2 \pi)⟨δni(x,t)⟩=∫δfi(R,μ,v∥,t)dvdRdξδ(R+ρi−x)/(2π)
Electron dynamics is assumed as adiabatic for simplicity, and the elec tron perturbed density is
为简单起见，假设电子动力学为绝热，电子加速器扰动密度为
δ n e ( x , t ) = e δ ϕ T e 0 n e 0 . δ n e ( x , t ) = e δ ϕ T e 0 n e 0 . deltan_(e)(x,t)=(e delta phi)/(T_(e0))n_(e0).\delta n_{e}(\boldsymbol{x}, t)=\frac{e \delta \phi}{T_{e 0}} n_{e 0} .δne(x,t)=eδϕTe0ne0.
The gyrokinetic Poisson's equation is 14 14 ^(14){ }^{14}14
陀螺动力学泊松方程为 14 14 ^(14){ }^{14}14
Z i 2 n i 0 T i 0 ( δ ϕ − δ ϕ ~ ) = Z i ⟨ δ n i ⟩ − e δ n e Z i 2 n i 0 T i 0 ( δ ϕ − δ ϕ ~ ) = Z i δ n i − e δ n e (Z_(i)^(2)n_(i0))/(T_(i0))(delta phi- widetilde(delta phi))=Z_(i)(:deltan_(i):)-e deltan_(e)\frac{Z_{i}^{2} n_{i 0}}{T_{i 0}}(\delta \phi-\widetilde{\delta \phi})=Z_{i}\left\langle\delta n_{i}\right\rangle-e \delta n_{e}Zi2ni0Ti0(δϕ−δϕ~)=Zi⟨δni⟩−eδne
where n i 0 n i 0 n_(i0)n_{i 0}ni0 and T i 0 T i 0 T_(i0)T_{i 0}Ti0 are the equilibrium ion density and temperature. δ ϕ ~ δ ϕ ~ widetilde(delta phi)\widetilde{\delta \phi}δϕ~ is the double gyrophase average of electrostatic potential for ion species, which is given as
其中 n i 0 n i 0 n_(i0)n_{i 0}ni0 和 T i 0 T i 0 T_(i0)T_{i 0}Ti0 是平衡离子密度和温度。 δ ϕ ~ δ ϕ ~ widetilde(delta phi)\widetilde{\delta \phi}δϕ~ 是离子种类静电势的双陀螺相平均值，其表示为
δ ϕ ~ ( x , t ) = 1 n i 0 ∫ f i 0 ( R , μ , v ∥ ) ⟨ δ ϕ ⟩ d v d R d ξ δ ( R + ρ i − x ) / ( 2 π ) , δ ϕ ~ ( x , t ) = 1 n i 0 ∫ f i 0 R , μ , v ∥ ⟨ δ ϕ ⟩ d v d R d ξ δ R + ρ i − x / ( 2 π ) , widetilde(delta phi)(x,t)=(1)/(n_(i0))intf_(i0)(R,mu,v_(||))(:delta phi:)dvdRd xi delta(R+rho_(i)-x)//(2pi),\widetilde{\delta \phi}(\boldsymbol{x}, t)=\frac{1}{n_{i 0}} \int f_{i 0}\left(\boldsymbol{R}, \mu, v_{\|}\right)\langle\delta \phi\rangle \boldsymbol{d} \boldsymbol{v} \boldsymbol{d} \boldsymbol{R} d \xi \delta\left(\boldsymbol{R}+\boldsymbol{\rho}_{\boldsymbol{i}}-\boldsymbol{x}\right) /(2 \pi),δϕ~(x,t)=1ni0∫fi0(R,μ,v∥)⟨δϕ⟩dvdRdξδ(R+ρi−x)/(2π),
where ∫ d v = 2 π B m i ∫ d v ∥ d μ ∫ d v = 2 π B m i ∫ d v ∥ d μ int dv=(2pi B)/(m_(i))int dv_(||)d mu\int \boldsymbol{d} v=\frac{2 \pi B}{m_{i}} \int d v_{\|} d \mu∫dv=2πBmi∫dv∥dμ. 其中 ∫ d v = 2 π B m i ∫ d v ∥ d μ ∫ d v = 2 π B m i ∫ d v ∥ d μ int dv=(2pi B)/(m_(i))int dv_(||)d mu\int \boldsymbol{d} v=\frac{2 \pi B}{m_{i}} \int d v_{\|} d \mu∫dv=2πBmi∫dv∥dμ .

11. B. Implementation of dynamic equations in cylindrical coordinates
11. B. 圆柱坐标中动力学方程的实现

The gyrocenter equation of motion is
陀螺中心运动方程为
R ˙ = v ∥ b ⏟ { parallel motion } + v ∥ 2 Ω c α ∇ × b ⏟ { curvaturedrift } + μ b × ∇ B m α Ω c α ⏟ { grad − B drift } + c b × ∇ ⟨ δ ϕ ⟩ B ⏟ { E × B drift } , R ˙ = v ∥ b ⏟ { parallel motion  } + v ∥ 2 Ω c α ∇ × b ⏟ { curvaturedrift  } + μ b × ∇ B m α Ω c α ⏟ { grad  − B  drift  } + c b × ∇ ⟨ δ ϕ ⟩ B ⏟ { E × B  drift  } , R^(˙)=ubrace(v_(||)bubrace)_({"parallel motion "})+ubrace((v_(||)^(2))/(Omega_(c alpha))grad xxbubrace)_({"curvaturedrift "})+ubrace((mubxx grad B)/(m_(alpha)Omega_(c alpha))ubrace)_({"grad "-B" drift "})ubrace(+(cbxx grad(:delta phi:))/(B)ubrace)_({E xx B" drift "}),\dot{\mathbf{R}}=\underbrace{v_{\|} \mathbf{b}}_{\{\text {parallel motion }\}}+\underbrace{\frac{v_{\|}^{2}}{\Omega_{c \alpha}} \nabla \times \mathbf{b}}_{\{\text {curvaturedrift }\}}+\underbrace{\frac{\mu \mathbf{b} \times \nabla B}{m_{\alpha} \Omega_{c \alpha}}}_{\{\text {grad }-B \text { drift }\}} \underbrace{+\frac{c \mathbf{b} \times \nabla\langle\delta \phi\rangle}{B}}_{\{E \times B \text { drift }\}},R˙=v∥b⏟{parallel motion }+v∥2Ωcα∇×b⏟{curvaturedrift }+μb×∇BmαΩcα⏟{grad −B drift }+cb×∇⟨δϕ⟩B⏟{E×B drift },
with 跟
v ∥ b ⏟ { parallel motion } = v ∥ ( b R R ^ + b Z Z ^ ) , v ∥ 2 Ω c α ∇ × b ⏟ { curvature drift } = [ 1 m α Ω c α m α v ∥ 2 B ( b Z ∂ B ∂ R − b R ∂ B ∂ Z ) + v ∥ 2 Ω c α B ( ∂ B R ∂ Z − ∂ B Z ∂ R ) ] ζ ^ , v ∥ b ⏟ { parallel motion  } = v ∥ b R R ^ + b Z Z ^ , v ∥ 2 Ω c α ∇ × b ⏟ { curvature drift  } = 1 m α Ω c α m α v ∥ 2 B b Z ∂ B ∂ R − b R ∂ B ∂ Z + v ∥ 2 Ω c α B ∂ B R ∂ Z − ∂ B Z ∂ R ζ ^ , {:[ubrace(v_(||)bubrace)_({"parallel motion "})=v_(||)(b_(R)( hat(R))+b_(Z)( hat(Z)))","],[ubrace((v_(||)^(2))/(Omega_(c alpha))grad xxbubrace)_({"curvature drift "})=[(1)/(m_(alpha)Omega_(c alpha))(m_(alpha)v_(||)^(2))/(B)(b_(Z)(del B)/(del R)-b_(R)(del B)/(del Z)):}],[{:+(v_(||)^(2))/(Omega_(c alpha)B)((delB_(R))/(del Z)-(delB_(Z))/(del R))] hat(zeta)","]:}\begin{aligned} \underbrace{v_{\|} \mathbf{b}}_{\{\text {parallel motion }\}}= & v_{\|}\left(b_{R} \hat{\mathbf{R}}+b_{Z} \hat{\mathbf{Z}}\right), \\ \underbrace{\frac{v_{\|}^{2}}{\Omega_{c \alpha}} \nabla \times \mathbf{b}}_{\{\text {curvature drift }\}}= & {\left[\frac{1}{m_{\alpha} \Omega_{c \alpha}} \frac{m_{\alpha} v_{\|}^{2}}{B}\left(b_{Z} \frac{\partial B}{\partial R}-b_{R} \frac{\partial B}{\partial Z}\right)\right.} \\ & \left.+\frac{v_{\|}^{2}}{\Omega_{c \alpha} B}\left(\frac{\partial B_{R}}{\partial Z}-\frac{\partial B_{Z}}{\partial R}\right)\right] \hat{\zeta}, \end{aligned}v∥b⏟{parallel motion }=v∥(bRR^+bZZ^),v∥2Ωcα∇×b⏟{curvature drift }=[1mαΩcαmαv∥2B(bZ∂B∂R−bR∂B∂Z)+v∥2ΩcαB(∂BR∂Z−∂BZ∂R)]ζ^,
μ b × ∇ B m α Ω c α ⏟ { grad − B drift } = μ m α Ω c α ( b Z ∂ B ∂ R − b R ∂ B ∂ Z ) ζ ^ , c b × ∇ ⟨ δ ϕ ⟩ B ⏟ E × B drift = c B [ − 1 R b Z ∂ ⟨ δ ϕ ⟩ ∂ ζ R ^ + ( b Z ∂ ⟨ δ ϕ ⟩ ∂ R − b R ∂ ⟨ δ ϕ ⟩ ∂ Z ) ζ ^ + 1 R b R ∂ ⟨ δ ϕ ⟩ ∂ ζ Z ^ ] , μ b × ∇ B m α Ω c α ⏟ { grad  − B  drift  } = μ m α Ω c α b Z ∂ B ∂ R − b R ∂ B ∂ Z ζ ^ , c b × ∇ ⟨ δ ϕ ⟩ B ⏟ E × B  drift  = c B − 1 R b Z ∂ ⟨ δ ϕ ⟩ ∂ ζ R ^ + b Z ∂ ⟨ δ ϕ ⟩ ∂ R − b R ∂ ⟨ δ ϕ ⟩ ∂ Z ζ ^ + 1 R b R ∂ ⟨ δ ϕ ⟩ ∂ ζ Z ^ , {:[ubrace((mubxx grad B)/(m_(alpha)Omega_(c alpha))ubrace)_({"grad "-B" drift "})=(mu)/(m_(alpha)Omega_(c alpha))(b_(Z)(del B)/(del R)-b_(R)(del B)/(del Z)) hat(zeta)","],[ubrace((cbxx grad(:delta phi:))/(B)ubrace)_(E xx B" drift ")=(c)/(B)[-(1)/(R)b_(Z)(del(:delta phi:))/(del zeta)( hat(R))+(b_(Z)(del(:delta phi:))/(del R)-b_(R)(del(:delta phi:))/(del Z))( hat(zeta)):}],[{:+(1)/(R)b_(R)(del(:delta phi:))/(del zeta)( hat(Z))]","]:}\begin{aligned} \underbrace{\frac{\mu \mathbf{b} \times \nabla B}{m_{\alpha} \Omega_{c \alpha}}}_{\{\text {grad }-B \text { drift }\}}= & \frac{\mu}{m_{\alpha} \Omega_{c \alpha}}\left(b_{Z} \frac{\partial B}{\partial R}-b_{R} \frac{\partial B}{\partial Z}\right) \hat{\zeta}, \\ \underbrace{\frac{c \mathbf{b} \times \nabla\langle\delta \phi\rangle}{B}}_{E \times B \text { drift }}= & \frac{c}{B}\left[-\frac{1}{R} b_{Z} \frac{\partial\langle\delta \phi\rangle}{\partial \zeta} \hat{\mathbf{R}}+\left(b_{Z} \frac{\partial\langle\delta \phi\rangle}{\partial R}-b_{R} \frac{\partial\langle\delta \phi\rangle}{\partial Z}\right) \hat{\zeta}\right. \\ & \left.+\frac{1}{R} b_{R} \frac{\partial\langle\delta \phi\rangle}{\partial \zeta} \hat{\mathbf{Z}}\right], \end{aligned}μb×∇BmαΩcα⏟{grad −B drift }=μmαΩcα(bZ∂B∂R−bR∂B∂Z)ζ^,cb×∇⟨δϕ⟩B⏟E×B drift =cB[−1RbZ∂⟨δϕ⟩∂ζR^+(bZ∂⟨δϕ⟩∂R−bR∂⟨δϕ⟩∂Z)ζ^+1RbR∂⟨δϕ⟩∂ζZ^],
and
v ˙ ∥ = − 1 m α 1 B ( B + B v ∥ Ω c α ∇ × b ) ⋅ ( Z α ∇ ⟨ δ ϕ ⟩ + μ ∇ B ) ⏟ { I } = − 1 m α b ⋅ ( Z α ∇ ⟨ δ ϕ ⟩ + μ ∇ B ) ⏟ { I } − 1 m α v ∥ Ω c α ∇ × b ⋅ ( Z α ∇ ⟨ δ ϕ ⟩ + μ ∇ B ) ⏟ { I I } = − Z α m α ∂ ⟨ δ ϕ ⟩ ∂ S − μ m α ( b R ∂ B ∂ R + b Z ∂ B ∂ Z ) ⏟ { I I } − Z α m α v ∥ Ω c α 1 R ∂ ⟨ δ ϕ ⟩ ∂ ζ [ ( b Z ∂ B ∂ R − b R ∂ B ∂ Z ) + ( ∂ B R ∂ Z − ∂ B Z ∂ R ) ] ⏟ . v ˙ ∥ = − 1 m α 1 B B + B v ∥ Ω c α ∇ × b ⋅ Z α ∇ ⟨ δ ϕ ⟩ + μ ∇ B ⏟ { I } = − 1 m α b ⋅ Z α ∇ ⟨ δ ϕ ⟩ + μ ∇ B ⏟ { I } − 1 m α v ∥ Ω c α ∇ × b ⋅ Z α ∇ ⟨ δ ϕ ⟩ + μ ∇ B ⏟ { I I } = − Z α m α ∂ ⟨ δ ϕ ⟩ ∂ S − μ m α b R ∂ B ∂ R + b Z ∂ B ∂ Z ⏟ { I I } − Z α m α v ∥ Ω c α 1 R ∂ ⟨ δ ϕ ⟩ ∂ ζ b Z ∂ B ∂ R − b R ∂ B ∂ Z + ∂ B R ∂ Z − ∂ B Z ∂ R ⏟ . {:[v^(˙)_(||)=ubrace(-(1)/(m_(alpha))(1)/(B)(B+(Bv_(||))/(Omega_(c alpha))grad xxb)*(Z_(alpha)grad(:delta phi:)+mu grad B)ubrace)_({I})],[=ubrace(-(1)/(m_(alpha))b*(Z_(alpha)grad(:delta phi:)+mu grad B)ubrace)_({I})ubrace(-(1)/(m_(alpha))(v_(||))/(Omega_(c alpha))grad xxb*(Z_(alpha)grad(:delta phi:)+mu grad B)ubrace)_({II})],[=ubrace(-(Z_(alpha))/(m_(alpha))(del(:delta phi:))/(del S)-(mu)/(m_(alpha))(b_(R)(del B)/(del R)+b_(Z)(del B)/(del Z))ubrace)_({II})],[ubrace(-(Z_(alpha))/(m_(alpha))(v_(||))/(Omega_(c alpha))(1)/(R)(del(:delta phi:))/(del zeta)[(b_(Z)(del B)/(del R)-b_(R)(del B)/(del Z))+((delB_(R))/(del Z)-(delB_(Z))/(del R))]ubrace).]:}\begin{aligned} \dot{v}_{\|} & =\underbrace{-\frac{1}{m_{\alpha}} \frac{1}{B}\left(\mathbf{B}+\frac{B v_{\|}}{\Omega_{c \alpha}} \nabla \times \mathbf{b}\right) \cdot\left(Z_{\alpha} \nabla\langle\delta \phi\rangle+\mu \nabla B\right)}_{\{I\}} \\ & =\underbrace{-\frac{1}{m_{\alpha}} \mathbf{b} \cdot\left(Z_{\alpha} \nabla\langle\delta \phi\rangle+\mu \nabla B\right)}_{\{I\}} \underbrace{-\frac{1}{m_{\alpha}} \frac{v_{\|}}{\Omega_{c \alpha}} \nabla \times \mathbf{b} \cdot\left(Z_{\alpha} \nabla\langle\delta \phi\rangle+\mu \nabla B\right)}_{\{I I\}} \\ & =\underbrace{-\frac{Z_{\alpha}}{m_{\alpha}} \frac{\partial\langle\delta \phi\rangle}{\partial S}-\frac{\mu}{m_{\alpha}}\left(b_{R} \frac{\partial B}{\partial R}+b_{Z} \frac{\partial B}{\partial Z}\right)}_{\{I I\}} \\ & \underbrace{-\frac{Z_{\alpha}}{m_{\alpha}} \frac{v_{\|}}{\Omega_{c \alpha}} \frac{1}{R} \frac{\partial\langle\delta \phi\rangle}{\partial \zeta}\left[\left(b_{Z} \frac{\partial B}{\partial R}-b_{R} \frac{\partial B}{\partial Z}\right)+\left(\frac{\partial B_{R}}{\partial Z}-\frac{\partial B_{Z}}{\partial R}\right)\right]} . \end{aligned}v˙∥=−1mα1B(B+Bv∥Ωcα∇×b)⋅(Zα∇⟨δϕ⟩+μ∇B)⏟{I}=−1mαb⋅(Zα∇⟨δϕ⟩+μ∇B)⏟{I}−1mαv∥Ωcα∇×b⋅(Zα∇⟨δϕ⟩+μ∇B)⏟{II}=−Zαmα∂⟨δϕ⟩∂S−μmα(bR∂B∂R+bZ∂B∂Z)⏟{II}−Zαmαv∥Ωcα1R∂⟨δϕ⟩∂ζ[(bZ∂B∂R−bR∂B∂Z)+(∂BR∂Z−∂BZ∂R)]⏟.
The weight equation is
权重方程为
where 哪里
− v E ⋅ ∇ f i 0 f i 0 | v ⊥ ⏟ wdrive = c B 0 b Z 1 R ∂ ⟨ δ ϕ ⟩ ∂ ζ [ 1 n i 0 ∂ n i 0 ∂ R + ( ϵ T i 0 − 1.5 ) 1 T i 0 ∂ T i 0 ∂ R ] − c B 0 b R 1 R ∂ ⟨ δ ϕ ⟩ ∂ ζ [ 1 n i 0 ∂ n i 0 ∂ Z + ( ϵ T i 0 − 1.5 ) 1 T i 0 ∂ T i 0 ∂ Z ] , − v E ⋅ ∇ f i 0 f i 0 v ⊥ ⏟ wdrive  = c B 0 b Z 1 R ∂ ⟨ δ ϕ ⟩ ∂ ζ 1 n i 0 ∂ n i 0 ∂ R + ϵ T i 0 − 1.5 1 T i 0 ∂ T i 0 ∂ R − c B 0 b R 1 R ∂ ⟨ δ ϕ ⟩ ∂ ζ 1 n i 0 ∂ n i 0 ∂ Z + ϵ T i 0 − 1.5 1 T i 0 ∂ T i 0 ∂ Z , {:[ubrace(-v_(E)*(gradf_(i0))/(f_(i0))|_(v_(_|_))ubrace)_("wdrive ")=(c)/(B_(0))b_(Z)(1)/(R)(del(:delta phi:))/(del zeta)[(1)/(n_(i0))(deln_(i0))/(del R)+((epsilon)/(T_(i0))-1.5)(1)/(T_(i0))(delT_(i0))/(del R)]],[-(c)/(B_(0))b_(R)(1)/(R)(del(:delta phi:))/(del zeta)[(1)/(n_(i0))(deln_(i0))/(del Z)+((epsilon)/(T_(i0))-1.5)(1)/(T_(i0))(delT_(i0))/(del Z)]","]:}\begin{aligned} \underbrace{-\left.\boldsymbol{v}_{\boldsymbol{E}} \cdot \frac{\nabla f_{i 0}}{f_{i 0}}\right|_{v_{\perp}}}_{\text {wdrive }}= & \frac{c}{B_{0}} b_{Z} \frac{1}{R} \frac{\partial\langle\delta \phi\rangle}{\partial \zeta}\left[\frac{1}{n_{i 0}} \frac{\partial n_{i 0}}{\partial R}+\left(\frac{\epsilon}{T_{i 0}}-1.5\right) \frac{1}{T_{i 0}} \frac{\partial T_{i 0}}{\partial R}\right] \\ & -\frac{c}{B_{0}} b_{R} \frac{1}{R} \frac{\partial\langle\delta \phi\rangle}{\partial \zeta}\left[\frac{1}{n_{i 0}} \frac{\partial n_{i 0}}{\partial Z}+\left(\frac{\epsilon}{T_{i 0}}-1.5\right) \frac{1}{T_{i 0}} \frac{\partial T_{i 0}}{\partial Z}\right], \end{aligned}−vE⋅∇fi0fi0|v⊥⏟wdrive =cB0bZ1R∂⟨δϕ⟩∂ζ[1ni0∂ni0∂R+(ϵTi0−1.5)1Ti0∂Ti0∂R]−cB0bR1R∂⟨δϕ⟩∂ζ[1ni0∂ni0∂Z+(ϵTi0−1.5)1Ti0∂Ti0∂Z],
with ϵ = 0.5 m i v ∥ 2 + μ B ϵ = 0.5 m i v ∥ 2 + μ B epsilon=0.5m_(i)v_(||)^(2)+mu B\epsilon=0.5 m_{i} v_{\|}^{2}+\mu Bϵ=0.5miv∥2+μB, 和 ϵ = 0.5 m i v ∥ 2 + μ B ϵ = 0.5 m i v ∥ 2 + μ B epsilon=0.5m_(i)v_(||)^(2)+mu B\epsilon=0.5 m_{i} v_{\|}^{2}+\mu Bϵ=0.5miv∥2+μB ，
Z i T i 0 v ∥ E ∥ ⏟ wpara = − Z i T i 0 v ∥ b ⋅ ∇ ⟨ δ ϕ ⟩ = − Z i v ∥ T i 0 ∂ ⟨ δ ϕ ⟩ ∂ S , Z i T i 0 v ∥ E ∥ ⏟ wpara  = − Z i T i 0 v ∥ b ⋅ ∇ ⟨ δ ϕ ⟩ = − Z i v ∥ T i 0 ∂ ⟨ δ ϕ ⟩ ∂ S , ubrace((Z_(i))/(T_(i0))v_(||)E_(||)ubrace)_("wpara ")=-(Z_(i))/(T_(i0))v_(||)b*grad(:delta phi:)=-(Z_(i)v_(||))/(T_(i0))(del(:delta phi:))/(del S),\underbrace{\frac{Z_{i}}{T_{i 0}} v_{\|} E_{\|}}_{\text {wpara }}=-\frac{Z_{i}}{T_{i 0}} v_{\|} \boldsymbol{b} \cdot \nabla\langle\delta \phi\rangle=-\frac{Z_{i} v_{\|}}{T_{i 0}} \frac{\partial\langle\delta \phi\rangle}{\partial S},ZiTi0v∥E∥⏟wpara =−ZiTi0v∥b⋅∇⟨δϕ⟩=−Ziv∥Ti0∂⟨δϕ⟩∂S,
with S = S c × L S = S c × L S=S_(c)xx LS=S_{c} \times LS=Sc×L or S = S S × L S = S S × L S=S_(S)xx LS=S_{S} \times LS=SS×L, and L L LLL is the field line length, and
with S = S c × L S = S c × L S=S_(c)xx LS=S_{c} \times LS=Sc×L 或 S = S S × L S = S S × L S=S_(S)xx LS=S_{S} \times LS=SS×L ，是 L L LLL 场线长度，并且
Z i T i 0 v d ⋅ E ⊥ ⏟ wdrift = − Z i T i 0 ( μ b × ∇ B m i Ω c i + v ∥ 2 Ω c α ∇ × b ) ⋅ ∇ ⟨ δ ϕ ⟩ = − Z i T i 0 1 R 1 B ∂ ⟨ δ ϕ ⟩ ∂ ζ [ c Z i ( μ + m i v ∥ 2 B ) ( b Z ∂ B ∂ R − b R ∂ B ∂ Z ) + v ∥ 2 Ω c i ( ∂ B R ∂ Z − ∂ B Z ∂ R ) ] . Z i T i 0 v d ⋅ E ⊥ ⏟ wdrift  = − Z i T i 0 μ b × ∇ B m i Ω c i + v ∥ 2 Ω c α ∇ × b ⋅ ∇ ⟨ δ ϕ ⟩ = − Z i T i 0 1 R 1 B ∂ ⟨ δ ϕ ⟩ ∂ ζ c Z i μ + m i v ∥ 2 B b Z ∂ B ∂ R − b R ∂ B ∂ Z + v ∥ 2 Ω c i ∂ B R ∂ Z − ∂ B Z ∂ R . {:[ubrace((Z_(i))/(T_(i0))v_(d)*E_(_|_)ubrace)_("wdrift ")=-(Z_(i))/(T_(i0))((mu b xx grad B)/(m_(i)Omega_(ci))+(v_(||)^(2))/(Omega_(c alpha))grad xxb)*grad(:delta phi:)],[=-(Z_(i))/(T_(i0))(1)/(R)(1)/(B)(del(:delta phi:))/(del zeta)[(c)/(Z_(i))(mu+(m_(i)v_(||)^(2))/(B))(b_(Z)(del B)/(del R)-b_(R)(del B)/(del Z)):}],[{:+(v_(||)^(2))/(Omega_(ci))((delB_(R))/(del Z)-(delB_(Z))/(del R))].]:}\begin{aligned} \underbrace{\frac{Z_{i}}{T_{i 0}} \boldsymbol{v}_{\boldsymbol{d}} \cdot \boldsymbol{E}_{\perp}}_{\text {wdrift }}= & -\frac{Z_{i}}{T_{i 0}}\left(\frac{\mu \boldsymbol{b} \times \nabla B}{m_{i} \Omega_{c i}}+\frac{v_{\|}^{2}}{\Omega_{c \alpha}} \nabla \times \mathbf{b}\right) \cdot \nabla\langle\delta \phi\rangle \\ = & -\frac{Z_{i}}{T_{i 0}} \frac{1}{R} \frac{1}{B} \frac{\partial\langle\delta \phi\rangle}{\partial \zeta}\left[\frac{c}{Z_{i}}\left(\mu+\frac{m_{i} v_{\|}^{2}}{B}\right)\left(b_{Z} \frac{\partial B}{\partial R}-b_{R} \frac{\partial B}{\partial Z}\right)\right. \\ & \left.+\frac{v_{\|}^{2}}{\Omega_{c i}}\left(\frac{\partial B_{R}}{\partial Z}-\frac{\partial B_{Z}}{\partial R}\right)\right] . \end{aligned}ZiTi0vd⋅E⊥⏟wdrift =−ZiTi0(μb×∇BmiΩci+v∥2Ωcα∇×b)⋅∇⟨δϕ⟩=−ZiTi01R1B∂⟨δϕ⟩∂ζ[cZi(μ+miv∥2B)(bZ∂B∂R−bR∂B∂Z)+v∥2Ωci(∂BR∂Z−∂BZ∂R)].
In GTC-X, the gyro-average is performed analytically by multiplying the fields with Bessel function as ⟨ δ ϕ ⟩ = δ ϕ J 0 ( k ζ ρ i ) ⟨ δ ϕ ⟩ = δ ϕ J 0 k ζ ρ i (:delta phi:)=delta phiJ_(0)(k_(zeta)rho_(i))\langle\delta \phi\rangle=\delta \phi J_{0}\left(k_{\zeta} \rho_{i}\right)⟨δϕ⟩=δϕJ0(kζρi) and ⟨ δ n i ⟩ = δ n i J 0 ( k ζ ρ i ) δ n i = δ n i J 0 k ζ ρ i (:deltan_(i):)=deltan_(i)J_(0)(k_(zeta)rho_(i))\left\langle\delta n_{i}\right\rangle=\delta n_{i} J_{0}\left(k_{\zeta} \rho_{i}\right)⟨δni⟩=δniJ0(kζρi), where ρ i = v t h , i / Ω c i ρ i = v t h , i / Ω c i rho_(i)=v_(th,i)//Omega_(ci)\rho_{i}=v_{t h, i} / \Omega_{c i}ρi=vth,i/Ωci is the ion gyroradius, v t h , i = T i / m i , Ω c i = Z i B / c m i , k ζ = n / R v t h , i = T i / m i , Ω c i = Z i B / c m i , k ζ = n / R v_(th,i)=sqrt(T_(i)//m_(i)),Omega_(ci)=Z_(i)B//cm_(i),k_(zeta)=n//Rv_{t h, i}=\sqrt{T_{i} / m_{i}}, \Omega_{c i}=Z_{i} B / c m_{i}, k_{\zeta}=n / Rvth,i=Ti/mi,Ωci=ZiB/cmi,kζ=n/R is the toroidal wave vector, n n nnn is the toroidal mode number, and R R RRR is the radial position, which is only valid for single toroidal mode simulation. The radial component of perpendicular wave vector k r k r k_(r)k_{r}kr is not considered for the gyroaverage in the current implementation for simplicity. In the future work, a more realistic gyro-average using the 4 point average method 21 21 ^(21){ }^{21}21 will be implemented in GTC-X for multiple mode simulation, which is used in another FRC drift wave code ANC.
在 GTC-X 中，陀螺平均值是通过将贝塞尔函数为 ⟨ δ ϕ ⟩ = δ ϕ J 0 ( k ζ ρ i ) ⟨ δ ϕ ⟩ = δ ϕ J 0 k ζ ρ i (:delta phi:)=delta phiJ_(0)(k_(zeta)rho_(i))\langle\delta \phi\rangle=\delta \phi J_{0}\left(k_{\zeta} \rho_{i}\right)⟨δϕ⟩=δϕJ0(kζρi) 和 ⟨ δ n i ⟩ = δ n i J 0 ( k ζ ρ i ) δ n i = δ n i J 0 k ζ ρ i (:deltan_(i):)=deltan_(i)J_(0)(k_(zeta)rho_(i))\left\langle\delta n_{i}\right\rangle=\delta n_{i} J_{0}\left(k_{\zeta} \rho_{i}\right)⟨δni⟩=δniJ0(kζρi) 的场相乘来执行的，其中 ρ i = v t h , i / Ω c i ρ i = v t h , i / Ω c i rho_(i)=v_(th,i)//Omega_(ci)\rho_{i}=v_{t h, i} / \Omega_{c i}ρi=vth,i/Ωci 是离子陀螺半径， v t h , i = T i / m i , Ω c i = Z i B / c m i , k ζ = n / R v t h , i = T i / m i , Ω c i = Z i B / c m i , k ζ = n / R v_(th,i)=sqrt(T_(i)//m_(i)),Omega_(ci)=Z_(i)B//cm_(i),k_(zeta)=n//Rv_{t h, i}=\sqrt{T_{i} / m_{i}}, \Omega_{c i}=Z_{i} B / c m_{i}, k_{\zeta}=n / Rvth,i=Ti/mi,Ωci=ZiB/cmi,kζ=n/R 是环形波矢量， n n nnn 是环形模态数， R R RRR 是径向位置，仅对单环形模态模拟有效。为简单起见，在当前实现中，陀螺平均值不考虑垂直波矢量 k r k r k_(r)k_{r}kr 的径向分量。在未来的工作中，将在 GTC-X 中使用 4 点平均法 21 21 ^(21){ }^{21}21 实现更真实的陀螺仪平均值，用于多模仿真，该仿真用于另一个 FRC 漂移波代码 ANC。
The comparisons between gyrocenter and fully kinetic particle trajectories with different particle energies and locations are shown in Fig. 4. Both gyrocenter and fully kinetic particle pitch angles at the outer midplane are tan ⁡ θ = v ∥ / v ⊥ = 1 tan ⁡ θ = v ∥ / v ⊥ = 1 tan theta=v_(||)//v_(_|_)=1\tan \theta=v_{\|} / v_{\perp}=1tan⁡θ=v∥/v⊥=1, which are trapped between two mirror throats. It is seen that gyrokinetic description is suitable for SOL region simulation with low temperature and high magnetic field, and its fidelity decreases with increasing the particle energy and moving toward the core region. Recently, a theoretical work has illustrated that the fidelity of the gyrokinetic equation is well achieved in the nonuniform magnetic field with ϵ = ρ / L < 1 , 40 ϵ = ρ / L < 1 , 40 epsilon=rho//L < 1,^(40)\epsilon=\rho / L<1,{ }^{40}ϵ=ρ/L<1,40 where ρ ρ rho\rhoρ is the gyroradius and L L LLL is the magnetic field gradient scale length. With regard to FRC geometry, the gyrokinetic description of thermal ion can capture the gyro-orbit but not betatron and figure- 8 orbit. 41 41 ^(41){ }^{41}41 However, a strong finite Larmor radius effect stabilizes the drift wave instability near or inside the core region, which has been observed in experiment 10 10 ^(10){ }^{10}10 and by previous local simulations. 25 − 27 25 − 27 ^(25-27){ }^{25-27}25−27 Thus, in this paper, we focus on gyrokinetic simulation of ITG instability in FRC as the first step to demonstrate the code capability for GTC-X as well as reveal the global nature of ITG instability in the SOL region.
具有不同粒子能量和位置的陀螺中心和全动力学粒子轨迹之间的比较如图4所示。外中平面的陀螺中心和全动力学粒子俯仰角都被困在 tan ⁡ θ = v ∥ / v ⊥ = 1 tan ⁡ θ = v ∥ / v ⊥ = 1 tan theta=v_(||)//v_(_|_)=1\tan \theta=v_{\|} / v_{\perp}=1tan⁡θ=v∥/v⊥=1 两个镜喉之间。可以看出，陀螺动力学描述适用于低温高磁场的SOL区域模拟，其保真度随着粒子能量的增加和向核心区域移动而降低。最近，一项理论工作表明，陀螺动力学方程的保真度在非均匀磁场中得到了很好的实现， ϵ = ρ / L < 1 , 40 ϵ = ρ / L < 1 , 40 epsilon=rho//L < 1,^(40)\epsilon=\rho / L<1,{ }^{40}ϵ=ρ/L<1,40 其中 ρ ρ rho\rhoρ 陀螺半径和 L L LLL 磁场梯度尺度长度。关于FRC几何形状，热离子的陀螺动力学描述可以捕获陀螺轨道，但不能捕获βtron和8字形轨道。 41 41 ^(41){ }^{41}41 然而，强烈的有限拉莫尔半径效应稳定了核心区域附近或内部的漂移波不稳定性，这在实验 10 10 ^(10){ }^{10}10 和先前的局部模拟中已经观察到。 25 − 27 25 − 27 ^(25-27){ }^{25-27}25−27 因此，本文重点研究了FRC中ITG不稳定性的陀螺动力学仿真，作为展示GTC-X代码能力的第一步，并揭示了SOL区域ITG不稳定性的全局性质。

12. Poisson solver 12. 泊松求解器

Poisson's equation is solved in a semispectral form. Applying Padé approximation, δ ϕ ^ ≈ 1 1 + k ⊥ 2 ρ i 2 δ ϕ δ ϕ ^ ≈ 1 1 + k ⊥ 2 ρ i 2 δ ϕ widehat(delta phi)~~(1)/(1+k_(_|_)^(2)rho_(i)^(2))delta phi\widehat{\delta \phi} \approx \frac{1}{1+k_{\perp}^{2} \rho_{i}^{2}} \delta \phiδϕ^≈11+k⊥2ρi2δϕ, Eq. (26) can be written as
泊松方程以半光谱形式求解。应用帕德近似， δ ϕ ^ ≈ 1 1 + k ⊥ 2 ρ i 2 δ ϕ δ ϕ ^ ≈ 1 1 + k ⊥ 2 ρ i 2 δ ϕ widehat(delta phi)~~(1)/(1+k_(_|_)^(2)rho_(i)^(2))delta phi\widehat{\delta \phi} \approx \frac{1}{1+k_{\perp}^{2} \rho_{i}^{2}} \delta \phiδϕ^≈11+k⊥2ρi2δϕ 方程（26）可以写成
e 2 n i 0 T e 0 [ − Z i 2 e 2 ρ s 2 ( 1 + e 2 Z i 2 n e 0 n i 0 T i T e ) ∇ ⊥ 2 + n e 0 n i 0 ] δ ψ = Z i δ n i 1 1 + e 2 Z i 2 n e 0 n i 0 T i 0 T e 0 , e 2 n i 0 T e 0 − Z i 2 e 2 ρ s 2 1 + e 2 Z i 2 n e 0 n i 0 T i T e ∇ ⊥ 2 + n e 0 n i 0 δ ψ = Z i δ n i 1 1 + e 2 Z i 2 n e 0 n i 0 T i 0 T e 0 , {:[(e^(2)n_(i0))/(T_(e0))[-(Z_(i)^(2))/(e^(2))rho_(s)^(2)(1+(e^(2))/(Z_(i)^(2))(n_(e0))/(n_(i0))(T_(i))/(T_(e)))grad_(_|_)^(2)+(n_(e0))/(n_(i0))]delta psi],[=Z_(i)deltan_(i)(1)/(1+(e^(2))/(Z_(i)^(2))(n_(e0))/(n_(i0))(T_(i0))/(T_(e0)))","]:}\begin{gathered} \frac{e^{2} n_{i 0}}{T_{e 0}}\left[-\frac{Z_{i}^{2}}{e^{2}} \rho_{s}^{2}\left(1+\frac{e^{2}}{Z_{i}^{2}} \frac{n_{e 0}}{n_{i 0}} \frac{T_{i}}{T_{e}}\right) \nabla_{\perp}^{2}+\frac{n_{e 0}}{n_{i 0}}\right] \delta \psi \\ =Z_{i} \delta n_{i} \frac{1}{1+\frac{e^{2}}{Z_{i}^{2}} \frac{n_{e 0}}{n_{i 0}} \frac{T_{i 0}}{T_{e 0}}}, \end{gathered}e2ni0Te0[−Zi2e2ρs2(1+e2Zi2ne0ni0TiTe)∇⊥2+ne0ni0]δψ=Ziδni11+e2Zi2ne0ni0Ti0Te0,
FIG. 4. Comparison of the trajectory between the fully kinetic particle and the gyrocenter for 44.4 e V 44.4 e V 44.4eV44.4 \mathrm{eV}44.4eV and 1022 e V 1022 e V 1022eV1022 \mathrm{eV}1022eV deuterium cases. The black solid line represents the separatrix, and green solid lines are the contour of magnetic flux. where δ ψ = δ ϕ − T i 0 Z i n i 0 δ n i / ( 1 + e 2 Z i 2 n e 0 n i 0 T i 0 T e 0 ) , ρ s = C s / Ω c i δ ψ = δ ϕ − T i 0 Z i n i 0 δ n i / 1 + e 2 Z i 2 n e 0 n i 0 T i 0 T e 0 , ρ s = C s / Ω c i delta psi=delta phi-(T_(i0))/(Z_(i)n_(i0))deltan_(i)//(1+(e^(2))/(Z_(i)^(2))(n_(e0))/(n_(i0))(T_(i0))/(T_(e0))),rho_(s)=C_(s)//Omega_(ci)\delta \psi=\delta \phi-\frac{T_{i 0}}{Z_{i} n_{i 0}} \delta n_{i} /\left(1+\frac{e^{2}}{Z_{i}^{2}} \frac{n_{e 0}}{n_{i 0}} \frac{T_{i 0}}{T_{e 0}}\right), \rho_{s}=C_{s} / \Omega_{c i}δψ=δϕ−Ti0Zini0δni/(1+e2Zi2ne0ni0Ti0Te0),ρs=Cs/Ωci, and C s C s C_(s)C_{s}Cs = T e / m i = T e / m i =sqrt(T_(e)//m_(i))=\sqrt{T_{e} / m_{i}}=Te/mi is the ion sound speed. By applying the Fourier transform in the toroidal direction to Eq. (31), we have
图 4.完全动力学粒子与氘 44.4 e V 44.4 e V 44.4eV44.4 \mathrm{eV}44.4eV 1022 e V 1022 e V 1022eV1022 \mathrm{eV}1022eV 形体的陀螺中心之间的轨迹比较。黑色实线表示分离线，绿色实线表示磁通量的轮廓。其中 δ ψ = δ ϕ − T i 0 Z i n i 0 δ n i / ( 1 + e 2 Z i 2 n e 0 n i 0 T i 0 T e 0 ) , ρ s = C s / Ω c i δ ψ = δ ϕ − T i 0 Z i n i 0 δ n i / 1 + e 2 Z i 2 n e 0 n i 0 T i 0 T e 0 , ρ s = C s / Ω c i delta psi=delta phi-(T_(i0))/(Z_(i)n_(i0))deltan_(i)//(1+(e^(2))/(Z_(i)^(2))(n_(e0))/(n_(i0))(T_(i0))/(T_(e0))),rho_(s)=C_(s)//Omega_(ci)\delta \psi=\delta \phi-\frac{T_{i 0}}{Z_{i} n_{i 0}} \delta n_{i} /\left(1+\frac{e^{2}}{Z_{i}^{2}} \frac{n_{e 0}}{n_{i 0}} \frac{T_{i 0}}{T_{e 0}}\right), \rho_{s}=C_{s} / \Omega_{c i}