Gyrokinetic continuum simulations of plasma turbulence in the Texas Helimak

TL;DR

This paper presents the first gyrokinetic continuum simulations of plasma turbulence in the Texas Helimak, capturing key features of the device and comparing results with experimental data to understand turbulence behavior.

Contribution

It introduces the first gyrokinetic simulations of plasma turbulence in the Helimak, including device-specific features like bad-curvature instabilities and sheath boundaries.

Findings

Simulations show qualitative agreement with experimental profiles and fluctuation amplitudes.

Both simulations and experiments exhibit blob-like turbulence statistics.

Quantitative differences highlight areas for further model refinement.

Abstract

The first gyrokinetic simulations of plasma turbulence in the Texas Helimak device, a simple magnetized torus, are presented. The device has features similar to the scrape-off layer region of tokamaks, such as bad-curvature-driven instabilities and sheath boundary conditions on the end plates, which are included in these simulations. Comparisons between simulations and measurements from the experiment show similarities, including equilibrium profiles and fluctuation amplitudes that approach experimental values, but also some important quantitative differences. Both experimental and simulation results exhibit turbulence statistics that are characteristic of blob transport.

Figures

Tables

Table 1. Table 1: Summary of Helimak simulation parameters

Background magnetic field	$B_0=0.1$ T
Background density	$n_0={10}^{16}$ m-3
Background electron temperature	$T_{e0}=10$ eV
Background ion temperature	$T_{i0}=1$ eV
Density gradient scale length	$L_n=0.1$ m
Electron gyroradius	${\rho }_e=1.02\times {10}^{-2}$ m
Electron mean free path	${\lambda }_{ee}$ = 98.4 m
Ion gyroradius	${\rho }_i=6.46\times {10}^{-3}$ m
Ion sound gyroradius	${\rho }_{\mathrm{s0}}=2.04\times {10}^{-2}$ m
Ion sound speed	$c_s=4.90\times {10}^3$ m/s
Ion mean free path	${\lambda }_{ii}$ = 1.39 m
Ion cyclotron frequency	${\Omega }_{ci}=2.40\times {10}^5$ s-1
Ion-ion collision frequency	${\nu }_{ii}=1.11\times {10}^3$ s-1
Ion-neutral collison frequency	${\nu }_{i0}\ge 185$ s-1

Equations

$$\frac{\partial {\mathcal{J}}_s{f}_s}{\partial t}$$

$$\{F,G\}$$

$$H_s=\frac{1}{2}m{v}_{\parallel }^2+\mu B+q_s⟨\varphi {⟩}_{\alpha },$$

$$-\nabla \cdot \left(\frac{n_{i0}^g{q}_i^2{\rho }_{\mathrm{s}0}^2}{T_{e0}}{\nabla }_{\perp }\varphi \right)={\sigma }_g=q_i{n}_i^g(\mathbf{R},t)-e{n}_e(\mathbf{R},t),$$

$$Z$$

$$\phi$$

$$S_s(\mathbf{R},v_{\parallel },\mu )=S_0\exp [-(x-x_{\mathrm{src}}{)}^2/(2{\sigma }_{\mathrm{src}}^2)]$$

$$\times F_M(v_{\parallel },\mu ;T_{s,\mathrm{src}}),$$

$$\frac{d\mathcal{N}}{dt}=\int d\mathbf{R}\left[S_n(x)-\frac{n(x,z)}{{\tau }_{\parallel }}\right]=0,$$

$$n(x,z)=n_p\exp \left[-(x-x_{\mathrm{src}}{)}^2/\left(2{\sigma }_n^2\right)\right]$$

$$\times \frac{1+\sqrt{1-z^2/(L_z/2{)}^2}}{2}.$$

$${\mathcal{E}}_s=\int d\mathbf{z}{H}_s{f}_s(\mathbf{z}).$$

$${\tilde{n}}_w(t_i)=\left[1-\text{cos}\left(\frac{2\pi }{N_t}i\right)\right]\tilde{n}(t_i).$$

$$\mathcal{P}(f)=\frac{⟨|\tilde{n}(f){|}^2⟩}{{\sum }_f\Delta f|\tilde{n}(f){|}^2},$$

Full text

I Introduction

In magnetically confined fusion devices such as tokamaks, the scrape-off layer (SOL) lies outside the last closed magnetic flux surface and consists of open magnetic field lines that intersect material walls. This region is characterized by fast transport along field lines toward the divertor plates or limiters and turbulent cross-field transport across field lines to the wall. Based on current experiments, the SOL heat-flux width is expected to be very narrow (only several gyroradii) and may not scale with machine size.1, 2 Increasing turbulence could mitigate this by diffusing the power over a wider area. Also, turbulent coherent structures called blobs (also referred to as filaments) 3, 4, 5, 6, 7, 8, 9 with higher temperatures and densities than the background plasma, are convected across field lines in the SOL and can collide with chamber walls before reaching the divertor plate, potentially damaging plasma-facing components and injecting impurities into the core plasma.10 Thus, the balance between parallel and cross-field transport in the SOL plays an important role in determining how heat and particles are exhausted in advanced fusion devices, 11, 12 and accurate modeling of this region is essential to predict the performance of future fusion reactors such as ITER.

The SOL in tokamaks is generally characterized by steep temperature and density gradients, geometry that includes an X-point, and complex physics, such as plasma-wall interactions, impurity transport, radiation, neutral recycling, and a plasma sheath. Basic plasma physics experiments with simpler magnetic configurations and lower temperature plasmas than tokamaks permit more comprehensive probe diagnostics and wider parameter scans than is generally available in fusion devices and are, thus, useful for validating analytical and numerical models of the SOL. For example, simple magnetized torus (SMT) experiments, such as the Texas Helimak13, 14 and TORPEX,15, 16, 17 are experimental approximations to a sheared slab geometry with bad curvature. They use vertical and toroidal field coils to create open, helical magnetic-field-line configurations with curvature and shear. Though SMT’s do not contain the complicated geometry of a tokamak SOL, such as an X-point, the dimensionless parameters (see table 2.3 in Ref. 18) and helical, open field lines are similar to that of a tokamak SOL (Fig. 1). Thus, experimental data from these devices can be compared with results from analytic and numerical models of SOL turbulence, as has been done with various fluid codes. 19, 20, 17, 21, 22

The main turbulent drives in the Helimak are the interchange and drift wave instabilities.13, 23, 24, 14, 18 Li et al. developed Helimak simulations to study these instabilities using an electrostatic drift-reduced Braginskii model with cold ions.19, 20 These fluid models reproduced some key experimental features, and more complete fluid models may explain SMT and tokamak-SOL plasmas well in some collisional parameter regimes. However, kinetic models are required to capture effects that can be significant in lower-collisionality parameter regimes, such as trapped particles, some non-linear wave–particle interactions, and non-Maxwellian features in the particle distribution functions.25, 26, 27 This paper details the first kinetic simulations of the Helimak. We use the Gkeyll computational framework to solve a five-dimensional (three spatial dimensions and two velocity dimensions) gyrokinetic equation, building on work by Shi et al. using straight and helical open-field-line configurations.28, 29, 30, 31 We maintain similarities with the set-up of the fluid simulations in Refs. 19 and 20 for comparison with those results.

In Refs. 30 and 31, Gkeyll simulations were performed with SMT geometry and NSTX-SOL parameters for a thin radial region ($\Delta R/R\sim 0.1$). These simulations neglected geometrical variation of SOL flux surfaces with poloidal angle, such as flux expansion and magnetic shear. Therefore, detailed comparisons with experimental data were not included, though there were interesting qualitative similarities. Here we use Gkeyll to model the Helimak, where the geometry of the simulation is much closer to that in the experiment, which only has a toroidal magnetic field plus a vertical magnetic field. The experiment has some magnetic shear13, 18 that is ignored at present in the simulation. The Helimak simulation also spans a much broader range of radius, $\Delta R/R\sim 0.5$. We include detailed comparisons with experimental measurements, including the time-averaged density and temperature profiles, turbulence fluctuation amplitudes, and correlation functions.

The ratio of the vertical magnetic field to toroidal magnetic field $B_{Z}/B_{\varphi}$ can be varied experimentally from 0.16 to 0.003 by changing resistance in the vertical field coils.24 The transition between various turbulence regimes in SMT’s as $B_{Z}/B_{\varphi}$ is varied has been studied numerically and experimentally.21, 22, 15, 16 In Ref. 22, Ricci and Rogers used the electrostatic, drift-reduced Braginskii equations32 to analyze the transition between different turbulence regimes, including ideal interchange, resistive-interchange, and drift-interchange, in SMT’s. They predicted only the first and last are present in the Helimak. For the ideal interchange instability, $k_{\parallel}=0$, and the dispersion relation predicts a growth rate of $\gamma_{I}\sim c_{s}/\sqrt{RL_{p}}$, where $c_{s}$ is the ion sound speed, $R$ is the major radius of the SMT, and $L_{p}$ is the pressure gradient scale length. In the high-pitch-angle configurations of the Helimak that we consider in the present paper, the experimentally measured parallel wavelength is typically much larger than the magnetic-field-line connection length, resulting in $k_{\parallel}\simeq 0$.18 Hence, the ideal interchange instability is expected to dominate. In the low-pitch-angle configuration, with finite $k_{\parallel}$, the drift-interchange instability is the main turbulent drive.

The plasma response to sheath physics in the Helimak sets up a radially-varying electric potential, producing a vertical $\bm{E}\times\bm{B}$ flow that varies with radius. A bias voltage can be applied to conducting plates on which magnetic field lines terminate to modify the vertical sheared flow and study the effect of velocity shear on turbulence suppression. In this paper, we only consider the case where all conducting plates are grounded and present the results of limiter-biasing simulations in a future paper.

This paper is organized as follows. In Sec. II, we review the model equations and parameters we use in our simulations. In Sec. III, we describe the results of our non-linear simulation, make comparisons with experimental data, identify the dominant instability, and briefly discuss kinetic effects. In Sec. IV, we summarize our findings and describe future directions of the code development to improve comparisons with experimental data.

II Model Equations

The Gkeyll code uses a nodal discontinuous Galerkin computational method for the spatial discretization and an explicit third-order Runge–Kutta method to discretize in time.33, 34 In our simulations, we solve the full-$f$ gyrokinetic equation in the long-wavelength, zero-Larmor-radius limit using the gyrocenter distribution function $f_{s}(\bm{R},v_{\parallel},\mu,t)$:

[TABLE]

where $s$ refers to the species, $C[f_{s}]$ is a simplified nonlinear Fokker-Planck collision operator, also called the Lenard-Bernstein or Dougherty operator, and $S_{s}$ is a source term, which we describe in more detail below. The Jacobian is $\mathcal{J}=B_{\parallel}^{*}$, where $B_{\parallel}^{*}=\bm{b}\cdot\bm{B}_{\parallel}^{*}$, and $\bm{B}_{\parallel}^{*}=\bm{B}+(Bv_{\parallel}/\Omega_{s})\nabla\times\bm{b}$. We use $\bm{B}_{\parallel}^{*}\simeq\bm{B}$. The phase-space advection velocities $\dot{\bm{R}}=\{\bm{R},H\}$ and $\dot{v}_{\parallel}=\{v_{\parallel},H\}$ are defined in terms of the Poisson bracket

[TABLE]

The gyrocenter Hamiltonian is

[TABLE]

where $\langle\rangle_{\alpha}$ is the gyro-average. In the long wavelength limit, $\langle\phi\rangle_{\alpha}=\phi$. In the continuous-time limit, the DG advection scheme we use is energy-conserving. For more details of our computational scheme, see Refs. 29, 30, 31.

We solve for the electrostatic potential using the long-wavelength, gyrokinetic Poisson equation with a linear ion polarization density

[TABLE]

where $\rho_{\mathrm{s}0}=c_{\mathrm{s}0}/\Omega_{i}$ is the ion sound gyroradius and $c_{\mathrm{s}0}=\sqrt{T_{e0}/m_{i}}$ is the ion sound speed. (Note that the left hand side is actually independent of electron temperature since $\rho_{s0}^{2}\propto T_{e0}$, but it is convenient to recover conventional normalizations so that $\nabla_{\perp}^{2}$ is normalized by $\rho_{s0}^{2}$, and $e\phi$ is normalized by $T_{e0}$.) In these simulations, we assume $q_{i}=e$, which appears to be a reasonable approximation since spectroscopic measurements indicate no ionzation states greater than $Z=1$ for an argon plasma in the Helimak. The background ion guiding-center density $n_{i0}^{g}$ is taken to be constant in space and time. Nonlinear polarization with a spatially- and time-dependent density is important future work.

We artificially increased the electron mass so that $m_{i}/m_{e}=\nobreak 400$, which is about 180 times smaller than the actual argon ion-to-electron mass ratio but twice the ratio used in previous fluid simulations.20 The electron mean free path, $\lambda_{ee}=v_{t,e}/\nu_{ee}$, is the same as experiment since the thermal speed and collision frequency contain the same mass dependency. The increased electron mass does significantly reduce electron parallel thermal conduction, but it is still fast enough that electron temperature is uniform along field lines.

We use a non-orthogonal field-line-following coordinate system35, 36, 37 to model the helical field lines: $z$ is the distance along the field, $x$ is the radial coordinate, and $y$ is the bi-normal coordinate. Field-line curvature enters through the second term in $\bm{B}^{*}$, which contributes terms proportional to $\nabla\times\bm{b}\cdot\nabla y=-1/x$. The field strength is dominated by the vacuum toroidal field; therefore, we assume $B(x)=B_{0}(R_{0}/x)$. Using the convention for field-line-following coordinate systems35, the direction of increasing $z$ points in the clockwise direction when viewed from above. The Helimak experimental magnetic field points in the counterclockwise, or $-z$, direction. Transformation to a cylindrical coordinate system $(R,\varphi,Z)$ is given by

[TABLE]

where $L_{y}$ is the extent of the domain in $y$ and $L_{z}$ is both the magnetic field line connection length and extent of the domain in $z$. See Fig. 2 for more details. The unit tangent vectors in the figure are defined as $\hat{\bf{e}}_{y}=\partial\bm{R}/\partial y/|\partial\bm{R}/\partial y|$ and $\hat{\bf{e}}_{z}=\partial\bm{R}/\partial z/|\partial\bm{R}/\partial z|$.

Conducting-sheath boundary conditions29, 30, 31 are used in the $z$ direction, where field-lines intersect conducting plates, for the distribution function $f$. These permit local parallel current fluctuations into and out of the wall. Since we assume quasineutrality in our model, we do not resolve the Debye sheath, which contains $n_{i}>n_{e}$ and is on the order of a few Debye lengths. The wall is assumed to be just outside the simulation domain. We solve for the sheath potential $\phi_{\mathrm{sh}}$ at the boundary using the gyrokinetic Poisson equation (4). This sets a cut-off parallel velocity for electrons, $\frac{1}{2}m_{e}v_{c}^{2}=e\phi_{\mathrm{sh}}$. Electrons with velocities greater than $v_{c}$ leave the domain, while those with smaller velocities are reflected. Dirichlet boundary conditions are used in $x$ for the potential ($\phi=0$), and periodic boundary conditions are used in $y$ for $f$ and $\phi$.

The Helimak simulations with Gkeyll evolve the electron and ion distribution functions for an argon discharge. We chose simulation parameters based on typical experimental parameters13, 23 and previous fluid simulations,19, 20 using $B_{0}=0.1$ T and $n_{i0}^{g}=n_{e0}=10^{16}$ m*-3*. Since the electron–ion thermal equilibration time is much longer than a charge-exchange time and particle confinement times for parallel loss, $T_{i0}\ll T_{e0}$ in the experiment. We use $T_{e0}=10$ eV based on experimental measurements. Previous work14, 19 estimates $T_{i0}\sim 0.1$ eV, though this has not been accurately measured. Argon that is ionized at different locations will be accelerated for different lengths of time towards the device ends by the pre-sheath potential, which varies along a field line. This will produce an ion temperature that is a significant fraction of the pre-sheath potential drop, which is of order $T_{e0}$. Therefore, we use $T_{i0}=1$ eV. The equilibrium temperatures define the velocity grid and appear in the the source terms, which we explain in more detail below. The particle temperatures $T_{i}$ and $T_{e}$ vary spatially and are evolved self-consistently in time. A gyrokinetic model is justified for these parameters since the ion-ion collision frequency, $\nu_{ii}=\nobreak 1.112\times\nobreak 10^{3}$ s*-1*, is much less than the ion cyclotron frequency, $\Omega_{ci}=\nobreak 2.40\times 10^{5}$ s*-1*, and the ion gyroradius, $\rho_{i}=\nobreak 6.46\times\nobreak 10^{-3}$ m, is much less than the density gradient scale length, $L_{n}=0.1$ m. See Table 1 for a summary of key simulation parameters.

The toroidal vacuum chamber of the Helimak is rectangular with a radial extent 0.6 m $\leq R\leq 1.6$ m and a vertical height $H=2.0$ m. By changing the current through the vertical field coils, the magnetic-field-line connection length $L_{z}$ can be varied from 12 m to 500 m. In our simulations, we use a slab-like metric with a spatial Jacobian that is independent of $x$. To simulate the entire volume of the Helimak, the configuration space extents are $x\in[0.6,1.6]$ m and $z\in[-L_{z}/2,L_{z}/2]$. $L_{y}$ varies inversely with the connection length: $L_{y}=H/N$, where $N=L_{z}/(2\pi R_{0})$ is the number of field-line turns and $R_{0}=1.1$ m. This gives $y\in[-L_{y}/2,L_{y}/2]$. In the experiment, there is magnetic shear due to the variation in the field with respect to $R$, which is not included in this model. This also means that the connection length $L_{z}$ and binormal length of periodicity $L_{y}$ do not vary in $x$ as they realistically should. The parallel velocity-space extents are given by $v_{\parallel,e}\in[-v_{e,\text{max}},v_{e,\text{max}}]$ and $v_{\parallel,i}\in[-v_{i,\text{max}},v_{i,\text{max}}]$, where $v_{e,\text{max}}=4v_{te}=4\sqrt{T_{e0}/m_{e}}$ and $v_{i,\text{max}}=6c_{\mathrm{s}}=6\sqrt{T_{e0}/m_{i}}$. The perpendicular velocity-space extents for each species are given by $\mu_{s}\in[0,3m_{s}v_{s,\text{max}}^{2}/(16B_{0})]$. The grid resolution is $(N_{x},N_{y},N_{z},N_{v_{\parallel}},N_{\mu})=(48,24,16,10,5)$ with uniform spacing. The numerical solutions for the distribution function and the Hamiltonian are projected onto piecewise-linear basis functions using the discontinuous Galerkin scheme. Because the gyrokinetic system is high in number of dimensions and, hence, computationally intensive, convergence tests were performed only in the $x$ direction for a low-resolution, ($N_{x}$,$N_{y}$,$N_{z}$,$N_{v_{\parallel}}$,$N_{\mu}$) = ($N_{x}$, 12, 8, 10, 5), and a higher-resolution, ($N_{x}$,$N_{y}$,$N_{z}$,$N_{v_{\parallel}}$,$N_{\mu}$) = ($N_{x}$, 24, 16, 10, 5), case. Unpublished resolution studies in other directions have been done in the past for various conditions. Power in the Fourier transform in $x$ of electron density fluctuations are compared for different resolutions in Fig. 3. Fluctuations are calculated as $\tilde{n}=n-\langle n\rangle_{y}$, where the last term denotes an average in the $y$-direction. The power has been averaged in $y$, $z$, and in time from 10 to 16 ms. The decay of the power at high $k_{x}$ indicates that the simulations are not changed significantly with increased resolution. Furthermore, equilibrium profiles and turbulence statistics were similar for the different cases.

In the Helimak, electrons are preferentially heated by electron-cyclotron and upper-hybrid resonances. The latter is believed to dominate and is localized from 1.0–1.1 m radially.18 The upper-hybrid resonance heating is difficult to model accurately due to its density dependence. The source is also challenging to model because the particle source rate and power deposition into the plasma are not accurately known. It is estimated that over 90$\%$ of the 6 kW input power is lost to radiation by electron-impact excitation of neutral argon. This process depends on the neutral density profile. We estimated the total radiative cooling using the predicted radiative cooling rate $L_{Z}$ in Ref. 38. The total power radiated in the device is given by $P_{\text{rad}}\mathcal{V}=L_{Z}n_{e}n_{\text{imp}}\mathcal{V}\simeq 5.5$ kW, assuming $T_{e}=10$ eV, $n_{e}=1\times 10^{16}$ m*-3*, and the total Argon density, including neutrals, is $n_{\text{imp}}=4\times 10^{17}$ m*-3*. $\mathcal{V}$ is the total volume of Helimak vacuum vessel. Neutral interactions and radiation are neglected in our present source model except indirectly through the usage of a net source rate, as described below. We estimate the ion-neutral collision frequency $\nu_{i0}\geq 185$ s*-1* from the electron-neutral frequency24, which differs by a factor of $\sqrt{m_{e}/m_{i}}$, assuming ions travel at the sound speed along field lines. This is much less than the gyrofrequency, though close to the ion transit frequency $L_{z}/(2c_{\mathrm{s}})=245$ s*-1*. Thus, the effect of neutrals is important and a feature we plan to include in future work.

We currently employ the simplified source model used in Refs. 29, 30, 31. For both species, we assume a source of the form

[TABLE]

where $s$ refers to species, $F_{M}$ is a normalized, non-drifting Maxwellian, $x_{\mathrm{src}}=1.0$ m is the source location, and $\sigma_{n}=0.01$ m is the source width. $T_{s,\mathrm{src}}$ is the source temperature. The source is independent of $y$, $z$, and $t$. Previous fluid simulations use $\sigma_{\mathrm{src}}\sim 0.1$ m. We chose a narrower source because the RF resonance is probably very narrow, though the location of the upper-hybrid resonance fluctuates as the density fluctuates. Equilibrium profiles are broadened by turbulence and wider than the assumed source width. Low-resolution scans in $\sigma_{\mathrm{src}}$ from 0.1 to 0.01 m produced little effect on resulting equilibrium profiles.

Since the actual particle and power source rates are not well known, we estimate the particle source rate $S_{0}$ and source temperature $T_{s,\mathrm{src}}$ using a 1D fluid transport model (see Appendix B in Ref. 30 or Appendix A in Ref. 31) to approximately produce the observed quasi-steady-state density and temperature profiles. We assume a balance of source rates and losses using

[TABLE]

where $\mathcal{N}$ is the total number of particles and $\tau_{\parallel}=\nobreak L_{z}/(2c_{\mathrm{s}})$ is the parallel transit time. We assume a steady-state density profile with a variation in $z$ based on a 1D steady-state single fluid calculation:

[TABLE]

The peak density $n_{p}$ is calculated such that $\int d\bm{R}\,\,n(x,z)/\mathcal{V}$ is equal to the background density. We define the density source for each species as $S_{s,n}(x)\equiv\int d\mathbf{v}\,S_{s}(\bm{R},\mu,v_{\parallel})$ and use Eqs. 8 and II to estimate $S_{0}\approx 9.77\times 10^{19}$ m*-3s-1* for the source amplitude. We estimate $T_{s,\mathrm{src}}=\frac{5}{3}T_{s0}$ for the source temperature in Eq. II based on the temperature at the midpoint in $z$ of a 1D steady-state fluid calculation.30, 31 However, that model assumes only convective outflows and neglects the fast parallel thermal transport of electrons, so it was necessary to double this value for electrons ($T_{e,\mathrm{src}}=\frac{10}{3}T_{e0})$ to approximately match experimental electron temperature profiles.

We expect the quasi-steady-state conditions of the simulation to be independent of initial conditions. Therefore, we chose initial conditions with a narrow density gradient scale length and a non-zero flow velocity as a function of the parallel coordinate, computed from simplified 1D fluid models following Ref. 30 and 31. The density initial conditions are pictured in the left plot of Fig. 4. We also used a density floor of 10$\%$ to minimize unphysical negative values in our full-$f$ distribution function. The code includes a scheme to correct for negative values of the distribution function, detailed in Refs. 29 and 30, which adds a small amount of numerical heating ($\sim 10\%$) to the simulations in the quasi-steady state. A new version of the code is under development and has a more robust method of preventing negative values of the distribution function.

Scans of the ion-to-electron mass ratio were performed in previous helical open-field-line simulations with Gkeyll and resulted in no significant changes to turbulence statistics.30, 31 The difference between the true and reduced mass ratio for argon is much larger in the case of the Helimak. The large parallel ion transit time requires longer simulation times to reach saturation and is computationally expensive. By solving for linear perturbations in the parallel current, $\tilde{j}_{\parallel}$, near the sheath and using $j_{\parallel e}+j_{\parallel i}=0$, we estimated that the reduced ion-to-electron mass ratio would under-predict the response of electrostatic potential to temperature fluctuations by approximately $50\%$. A comparison of low-resolution simulations, $(N_{x},N_{y},N_{z},N_{v_{\parallel}},N_{\mu})=(24,12,8,10,5)$, with the true and reduced mass ratio supports this prediction. A new modal discontinuous Galerkin version of Gkeyll with automatic code generation is under development and is much faster. It should allow routine use of higher ion-to-electron mass ratios in the future.

III Simulation results

We ran the simulation using the Skylake (SKX) compute nodes on the Stampede2 cluster at the Texas Advanced Computing Center. A simulation of a 16 ms argon discharge for $L_{z}=40$ m required approximately 180,000 CPU-hours. This is about 4 times the ion transit time, $\tau_{\parallel}=L_{z}/(2c_{\mathrm{s}})=4.084$ ms, which is the time it takes an ion traveling at its characteristic sound speed, with $T_{e0}=10$ eV, to traverse half the magnetic-field-line connection length. We compared with experimental data from a discharge with a similar connection length at $R_{0}=1.1$ m, corresponding to $B_{Z}\simeq 0.005$ T and $B_{Z}/B_{\varphi}\simeq 0.05$. We selected this configuration as it results in the shortest connection length for which magnetic field lines terminate on conducting plates on the bottom and top of the vacuum chamber.

Figure 5 shows the total particle count and total energy, $\mathcal{E}_{\mathrm{tot}}=\mathcal{E}_{e}+\mathcal{E}_{i}$, of the simulation. Particle energy for a single species is calculated as the integration of the product of the particle Hamiltonian and distribution function over all phase space:

[TABLE]

Particle count is still increasing, though beginning to level-off by the end of the simulation. The energy appears to saturate around 10 ms. The large spike in energy at early time is due to numerical heating added from the scheme to correct for negative distribution function values. The nearly constant energy but increasing particle count indicates a slight decrease in temperature. We calculate equilibrium profiles and turbulence statistics over the interval from 10 to 16 ms.

In the simulation, turbulent structures are visible at 0.15 ms, as seen in the center plot of figure 4, which contains snapshots of the ion guiding-center density in the non-orthogonal field-line-following coordinate system. The right-hand plot shows a non-linear turbulent state at 1 ms. According to the energy plot in Fig. 5, a quasi-steady state is reached at approximately 10 ms. Figure 6 (Multimedia view) depicts electron density, electron temperature, and electrostatic potential at $t=$ 10 ms in the field-line-following coordinate system at the midpoint along $z$. Using the transformation in Eqs. 5 and 6, this data was mapped to the cylindrical coordinate system $(R,\varphi,Z)$ and is shown in Fig. 7 in the $(R,Z)$ plane at a slice in $\varphi$. $N=L_{z}/(2\pi R_{0})=5.8$ is the number of toroidal turns of the field line, and this periodicity is apparent in the figure. For more details of the mapping, see Fig. 2.

Figure 8 compares equilibrium profiles of simulation and experimental data. Experimental Langmuir probes are located on conducting plates at the top and bottom of the vacuum vessel. The probe measurements of the equilibrium profiles have about a 50% uncertainty, and experimental data is generally fit to a smooth curve. There is a top–bottom asymmetry in experimental density profiles. This is not present in the simulation, since our model is symmetric in $z$. Thus, in all comparisons we calculate simulation data at the minimum of the $z$ domain ($z=-L_{z}/2$). One reason for this asymmetry could be the vertical, $Z$ in the lab frame or $z$ in the simulation, $\bm{E}\times\bm{B}$ drift, which we presently neglect in our model. Discussing this issue is complicated because of the non-orthogonal coordinate system that is used, but to be precise, the advection due to a radial electric field, $V_{E\times B}\cdot\nabla=V_{E\times B}\cdot\nabla y\,\partial/\partial y+V_{E\times B}\cdot\nabla z\,\partial/\partial z$ is approximated by $V_{E\times B}\cdot\nabla y\partial/\partial y$, which advects in $y$ at fixed $z$ and, thus, at fixed height $Z$. This is a common assumption for core flux-tube codes assuming $k_{\perp}\gg k_{||}$, but can break down for long wavelength components or small magnetic pitch angle. Based on experimental flow data from Ref. 23, the vertical $\bm{E}\times\bm{B}$ flow is 2–3 times the vertical component of the sonic outflows along the field lines for simulated parameters. This factor increases as the toroidal field component $B_{T}$ increases relative to the vertical component $B_{Z}$ and may influence observed density profiles. Additionally, the microwave waveguide enters the machine at the bottom and may explain the similarity in bottom density profiles for different pitch angles which is not observed for top profiles.18 Experimental temperature profiles are fairly top–bottom symmetric and centered around 1.1 m, irrespective of pitch angle.18 Since electron parallel heat conduction is rapid, the electron-temperature steady-state profiles would be less affected by the vertical $\bm{E}\times\bm{B}$ drift or the location of the microwave guide. The vertical component of $\bm{E}\times\bm{B}$ flow is likely important in transport balances, though not as important for turbulence fluctuations. A full analysis of the effects of mean vertical $\bm{E}\times\bm{B}$ flows on the asymmetry in the experimental profiles involves solving a complex non-linear equation and is not within the scope of this paper.

We model the source for each species using a single Maxwellian with a flat temperature profile and a density profile centered at 1.0 m. We chose the location of the source based on experimental steady-state density profile. Figure 8(a) depicts the narrow electron density source and resulting equilibrium profiles which were broadened by turbulence. The simulation density peak is close in magnitude to experimental profile peaks. However, the simulation fails to capture the steep gradient just to the right of the peak and the flat profile at higher radii that is present in the experiment. This might be improved by implementing a density-weighted Poisson equation, instead of the current one Eq. 4, in which the background ion guiding center density $n_{i0}^{g}$ is constant in space and time. This approximation is similar to the Boussinesq approximation used in many fluid models.19, 20, 39 Some simulations have relaxed this approximation40, 41, 42 but have not been tested with Helimak parameters.

Figure 8(b) shows that the simulation produced a temperature profile centered to the left of the experimental profile. The peak magnitude and the gradient to the right of the peak were both less than experiment. The simulated plasma potential profile is shown in the same plot as $e\phi/\Lambda$, where $\Lambda=\text{ln}\sqrt{m_{i}/(2\pi m_{e}})$ for the reduced ion-to-electron mass ratio and the true mass ratio. An estimate of the simulated floating potential using $\phi_{\text{fl}}\simeq\phi-T_{e}\Lambda/e$ agrees qualitatively with experimental measurements of the same quantity.13, 23, 18

The simulated turbulence profiles come closer to experimentally measured levels than previous fluid simulations.19 Turbulence levels are calculated as the root-mean-square of density fluctuations, $\tilde{n}=n-\langle n\rangle$, normalized to the mean electron density $\langle n\rangle$. In the Helimak, ion saturation current measurements, $I_{\rm{sat}}=n_{e}\sqrt{T_{e}}$, are used as a proxy for density in density fluctuation statistics and have very small associated errors. For comparison, we also calculate simulation density fluctuations using $I_{\rm{sat}}=n_{e}\sqrt{T_{e}}$ but, henceforth, refer to all density fluctuations as $\tilde{n}$, as opposed to $\tilde{I}_{\rm{sat}}$, as is typically done in Helimak experimental papers. In the simulation, angle brackets $\langle\rangle$ denote an average in time, performed from 10 to 16 ms, and an average in the azimuthal angle. Plot 9(a) is normalized to a local mean density, while 9(b) is normalized to a global mean density. The latter shows slightly better agreement, though both differ at large radii. Using a density-weighted Poisson equation may also improve agreement here.

The plasma potential equilibrium profile generates a sheared flow perpendicular to the magnetic field. This flow is calculated as $V_{E}=-(d\phi_{0}/dx)/B_{0}$, with the negative sign resulting from the direction of the experimental magnetic field, and the corresponding shear as $\gamma_{E}=|dV_{E}/dx|$. As mentioned above, we currently neglect the vertical $\bm{E}\times\bm{B}$ in our simulation, and the flow in 10(a) is in the $y$, or toroidal, direction. We approximate the experimental flow and shear from the experimental temperature profiles by assuming $\phi_{0}\approx\Lambda T_{e}/e$, though this calculation overestimates the flow velocities as compared with published experimental measurements.23, 20 We note that this flow is mostly in the vertical direction. Previous research has shown that velocity shear in strongly magnetized plasmas has a stabilizing effect by breaking apart turbulent eddies.43 Therefore, we compare shear to the local linear growth rate of the interchange instability $\gamma_{I}=c_{s}/\sqrt{RL_{p}}$ in Fig. 10(b). The pressure gradient scale length is defined as $1/L_{p}=-(dp/dx)/p$. Turbulence statistics, such as power spectra and auto-correlation functions, are calculated at radial locations where $V_{E}\simeq 0$, at $R=1.09$ m for the simulation and $R=1.17$ m for the experiment, and where $V_{E}$ is approximately maximal, at $R=1.20$ m for the simulation and $R=1.25$ m for the experiment.

The power spectra are calculated from the Fourier transform of the electron density fluctuations $\tilde{n}=n-\langle n\rangle$ with respect to time. Since the density fluctuations are not periodic at the boundaries $t_{\text{min}}$ and $t_{\text{max}}$ of this time series, we multiply the density fluctuations by a sinusoidal window prior to performing the FFT:

[TABLE]

The integer $i$ runs from 0 to $N_{t}$, ensuring that the time series is 0 and thus periodic. The power spectra is then normalized to the total energy

[TABLE]

where $f$ denotes frequency and not the distribution function. The frequency grid spacing $\Delta f$ is 168 Hz. The same window and sampling frequency were used for both simulation and experimental data. The experimental data uses a longer time series and was binned to 6 ms intervals and then averaged.

Figure 11 compares the experimental and simulated power spectra multiplied by the frequency, $f\mathcal{P}(f)$, at radial locations where $\bm{E}\times\bm{B}$ flow is zero and maximal. For the radial values where flow is maximal, the experimental spectrum is weakly peaked near 1.2 kHz, in accordance with previous results,14, 19 while the simulation is a bit noisier and peaks at a higher frequency, in the 2-4 kHz range. The simulation spectra appear to have less relative power in frequencies just below the peak and greater power above the peak when compared to the respective experimental spectra. Furthermore, simulation spectra decay more quickly at high frequencies.

Figure 12 depicts the autocorrelation functions, which are calculated as $C(\tau)=\langle\tilde{n}(t)\tilde{n}(t+\tau)\rangle/\langle\tilde{n}(t)^{2}\rangle$. Brackets denote a time average, and $\tau$ is the time difference. In both sets of data, the curves where $V_{E}$ is maximal exhibit more oscillatory behavior. The correlation time for the simulation data is shorter, as we would expect given the results from the power spectra.

Figure 13 compares the radial correlation length as a function of radius. The experimental and simulation graphs show opposite trends, though the magnitude is similar. Reduced turbulence levels have been shown to correspond to a reduction in the radial correlation length.23, 43 However, the correlation lengths are larger for the simulation at $R>1.1$ m, though the turbulence fluctuations are less than experiment in this region.

Figure 14 compares the probability density functions (PDF’s) at radial locations where $V_{E}$ is zero and maximal. Experimental values were measured at the top of the device and calculated as $\tilde{n}/\tilde{n}_{\text{rms}}$. The positive tails of the simulation PDF’s approach the experimental values until 5$\tilde{n}/\tilde{n}_{rms}$, where they end. Presumably, if the simulations were run for a longer time, the tail in the simulation PDF would extend to higher $\tilde{n}/\tilde{n}_{\rm rms}$. For radii where $V_{E}$ is maximal, negative tails in the simulation are greater than experiment. The simulations have lower density fluctuations at large radii, compared with experiment, which means that the negative fluctuations in density are less affected by the positivity of total density and could possibly explain this difference in the negative tails.

In figure 15(a) experimental and simulation density fluctuations exhibit positive skewness, which is a characteristic of blob transport.44 Skewness is calculated as $E[\tilde{n}^{3}]/\sigma^{3}$, where $E[...]$ is the expectation value and $\sigma$ is the standard deviation of the density fluctuations. There are similarities between the simulation and experimental measurements from the top at radii greater than 1.2 m. However, the simulation does not reproduce the large values observed in the experiment at $R=1.1$ m for the top measurements and at $R=1.0$ for the bottom measurements. We compare excess kurtosis in Fig. 15(b), which is calculated as $E[\tilde{n}^{4}]/\sigma^{4}-3$. Positive excess kurtosis is also a signature of blob transport.44 Though the simulation matches the top experimental curve at radii greater than 1.25 m, it fails to capture the large peaks in excess kurtosis that appear in both experimental data sets. We also note that values of skewness and kurtosis from Gkeyll simulations with NSTX-like parameters are comparable to those calculated here.30, 31

Figure 16 shows a plot of the electron density $(x,z)$ at the midpoint in $y$, and there is little structure along $z$ in the bad curvature region $R>1.0$ m, except in the overall fall-off of density approaching the conducting-sheath boundaries. This suggests interchange-like turbulence, with $k_{z}\approx 0$. Figure 17 depicts the power in the Fourier transform of the density fluctuations, $\tilde{n}=n-\langle n\rangle_{y}$, in $z$, averaged in $x$ and $y$ and in time from 10 to 16 ms. The zeroth mode is more than an order of magnitude greater than the other modes, suggesting that the average $k_{z}\approx 0$ and, hence, interchange-driven turbulence dominates. Limitations in computational resources prevented high-resolution comparisons with lower pitch angles, particularly because the parallel ion transit times are longer and require longer simulation times to saturate. Low-resolution test cases demonstrated greater $k_{z}$ structure as the pitch angle was decreased, possibly indicating a transition to the drift-interchange-driven turbulence regime as predicted in Refs. 16 and 22.

One way kinetic effects enter is via the fact that collisions that are not strong enough to fully thermalize electrons in the tail above the sheath potential cutoff velocity during a single bounce of electrons between the sheaths at the two ends of the field line. In some places this may reduce the electron tail and reduce the electron sheath heat-transmission coefficient, but in places where turbulence has carried energetic electrons from a hot region to a region that is colder on average, it may enhance the electron sheath heat-transmission coefficient. The coefficient is given by $\gamma_{e}=q_{\parallel}/(\Gamma_{\parallel,n}T_{e})$, with quantities evaluated at the sheath entrance, $z=L_{z}/2$, and at the midpoint in $y$. We plot the time-evolution of this quantity in figure 18(a) for radial values where $V_{E}$ is zero and maximal. For a Maxwellian distribution function, this coefficient is $2+e\phi_{\rm sh}/T_{e}$, and it is often assumed $e\phi_{\rm sh}/T_{e}\sim 3$.10 In Fig. 18(b), we normalize the electron sheath heat-transmission coefficient to this quantity, using the time-dependent quantities $\phi_{\rm sh}$ and $T_{e}$ rather than assuming $e\phi_{\rm sh}/T_{e}=3$. In the simulation, the sheath heat-transmission coefficient relative to the Maxwellian value fluctuates in time, exceeding the expected Maxwellian coefficient by approximately 10–20%, suggesting relatively small, but non-negligible kinetic effects for these parameters. However, a more detailed convergence study in velocity resolution would be needed to quantify how much of this is due to numerical discretization versus kinetic effects.

Kinetic effects may be much more important in other regimes, such as in some previous PIC simulations of ELM-type events45 on JET, where they observed that $\gamma_{e}$ was 5 to 6 times the nominal Maxwellian value. We leave more detailed analysis of kinetic effects to future simulations of such regimes and also when we extend the code, as planned, to handle the closed-field-line pedestal and SOL regions of tokamaks simultaneously.

IV Conclusion

This work represents the first extensive comparison of kinetic simulations of a toroidal open-field line basic plasma physics device with experimental data. The results demonstrate good progress towards modeling turbulence on helical open-field lines in tokamak SOL-like conditions with gyrokinetic equations. These results indicate which aspects of the model could be improved prior to simulating a more complicated tokamak SOL where kinetic effects are likely more significant.

Simulation density profiles approached Helimak experimental values of magnitude and gradient scale length but failed to caputure the top-bottom asymmetry. The simulation also under-predicted the magnitude and gradient of the temperature profiles. Turbulence fluctuation levels were closer to experiment than previous fluid simulations, particularly near $R=1$ m, but still less than experiment at higher radii. Some qualitative features of turbulence statistics were captured, such as oscillations in the autocorrelation function where $V_{E}$ is maximal, radial correlation lengths of similar magnitude, and positive skewness and excess kurtosis of the distribution function in the bad-curvature region. Notable differences were that the power spectra from the simulation peaked at higher frequencies than experiment, radial correlation values have opposing trends, and the simulation did not capture the largest values of skewness and excess kurtosis at some experimental radial locations. Previous research suggests that interchange turbulence should dominate for this connection length, and this was confirmed with an average $k_{\parallel}\approx 0$ in our simulations. Kinetic effects were observed in fluctuations of the electron sheath heat-transmission coefficient.

It is clear that key physical and geometric features need to be added to our model in the future for better overall agreement with experimental turbulence statistics. A more self-consistent source model that includes radiation and the density-dependent upper-hybrid resonance might improve the agreement of equilibrium profiles. Including a vertical ($z$-directed, in our non-orthogonal coordinates) $\bm{E}\times\bm{B}$ flow may account for some of the top-bottom asymmetry that is visible in the density equilibrium profiles. Since the experimental background density varies radially, a more accurate model would use the density-weighted Poisson equation (instead of 4), which includes the time and spatial variation in the ion guiding-center density on the left-hand side of the equation. Other geometric features to include are magnetic shear and the radial variation in the connection length. These may affect where the transition from interchange to drift-wave turbulence occurs by giving the turbulence shorter parallel correlation length and breaking up the strong oscillatory dynamics present in the simulation.

A new, faster version of the Gkeyll code is being developed and will allow for scans in the ion-to-electron mass ratio and collision frequencies to study their effect on turbulence statistics. Scans in connection length could also be performed to study the transition from interchange to drift-interchange turbulence. The effect of neutral particles is significant in this parameter regime, and future development for the code will likely include an appropriate neutral model. Improvements to the conducting-sheath boundary conditions will also be considered.

Acknowledgements.

T.N.B. would like to thank F. Ebrahimi for reviewing the manuscript and for helpful discussions on turbulence spectra and convergence testing. We also thank N. Mandell for helpful discussions about non-orthogonal coordinate system. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. We thank D.R. Hatch, J.M. Tenbarge, and F.L. Waelbroeck for contributing computational resources through their allocations on the Texas Advanced Computing Center Stampede2 cluster. We also thank J. Juno for helping set up Gkeyll on Stampede. This work was supported by the U.S. Department of Energy SCGSR program under contract DE-SC0014664, by DOE Contract DE-FG02-04ER-54742, through the Institute of Fusion Studies at the University of Texas at Austin, and DOE contract DE-AC02-09CH11466, through the Princeton Plasma Physics Laboratory (PPPL). G.W.H. and A.H. were supported in part by the SciDAC Partnership for Multiscale Gyrokinetics and the SciDAC Partnership for High-Fidelity Boundary Plasma Simulations. A.H. was also supported in part by the PPPL Laboratory Directed Research and Development program.