Understanding the impact of insulating and conducting endplate boundary conditions on turbulence in CSDX through nonlocal simulations

Edge and Scrape-Off Layer (SOL) turbulence has been a subject of interest in magnetically confined fusion for many decades. It is well-understood that turbulent transport is a major constraint in determining the performance of tokamaks.¹ Many experimental measurements, turbulence modeling, and simulations have tried to understand instabilities and transport in the SOL.^2–5 Nevertheless, the complicated geometry of larger toroidal magnetically confined devices makes detailed code-experiment validation studies challenging. Linear plasma devices such as the Controlled Shear Decorrelation Experiment (CSDX) are useful in studying some of the basic underlying physics due to their simple geometry. CSDX inherently operates at low temperature with plasma densities such that collisions are important, similar to the tokamak SOL. Moreover, these conditions permit the use of experimental techniques such as multi-tip Langmuir probes and fast imaging of turbulence physics, enabling measurement of particle and momentum fluxes, zonal flow formation, and spectral energy transfer across a range of conditions.^6–9 

There have been numerous previous simulation studies on linear plasma devices, including on VINETA,¹⁰ and the Large Plasma Device (LAPD)^11–13 to name a few. The previous simulations of CSDX by Holland et al.¹⁴ used a simple 2D Hasegawa-Wakatani model in the local limit to capture the basic physics of the resistive drift-wave instability. While that work showed self-regulating drift-wave/zonal flow dynamics similar to experiments, this simulation overestimated turbulence levels by an order of magnitude. D'Ippolito et al.¹⁵ simulated an axially averaged domain in slab geometry using the SOLT code, which differentiated insulating and conducting endplates by adding a sheath closure term. In these simulations, the time averaged quantities and frequency domain result of these simulations were in better qualitative agreement with experimental results.

This paper aims to build upon these previous CSDX simulations by adding three new key elements. First, a 3D capability to describe the axial dynamics of turbulence is added, including proper sheath boundary conditions. The choice of proper boundary conditions is known to be important in the edge simulations.¹⁰ For instance, Russell et al.¹⁶ show that sheath dissipation can be a sink for momentum transport in the SOL and can thereby affect turbulence levels. Second, these simulations evolve electron temperature fluctuations to enable synthetic diagnostics studies of probes and camera. Mahdizadeh et al.¹⁷ report experimental measurements from a low temperature plasma, in which the electron temperature fluctuations are small enough not to affect floating potential measurements significantly. However, the simulations by Nold et al.¹⁸ showed that in higher temperature plasmas, electron temperature fluctuations are more significant and the synthetic floating potential is strongly distorted when compared to the actual plasma potential. Before the work presented here, there has not been much focused simulation study on the effect of T_e fluctuations in low temperature plasmas. Third, unlike previous simulations of CSDX, radially varying plasma profiles and parameters are retained, since microturbulence gradient lengthscales are comparable to the device lengthscale.

The model in this paper is used to study the impact of changing endplates on the underlying instabilities in CSDX and will be the basis for validation studies. Here, we will investigate the physics differences of insulating versus conducting endplate on the operative linear instability and its nonlinear saturation. A more detailed code-experiment validation study incorporating quantification of uncertainties and synthetic diagnostics with validation metrics will be presented in a future paper.

The rest of this paper is organized as follows: In Sec. II, we discuss the experimental conditions of the CSDX linear device. In Sec. III, we introduce the turbulence model, calculation of plasma parameters, and code implementation. The physics of the linear instabilities operative in CSDX are discussed in Sec. IV. Sec. V documents differences in the nonlinear saturation of drift-wave turbulence with insulating and conducting endplates. We summarize the results in Sec. VI.

Motivated by a desire for a better understanding of magnetized plasma turbulence dynamics, this study of the CSDX linear plasma device aims to define a simple self-regulating description of drift-wave turbulence/zonal flow dynamics.¹⁹ Previously, CSDX experimental results⁶ showed that the dominant instability in CSDX is the density gradient-driven collisional drift-wave instability. Experimental measurements^8,20–22 showed these modes as saturated due to the nonlinear re-distribution of turbulent energy in CSDX via wave-wave coupling and energy transfer to large scale zonal flow without any apparent momentum source.

Thakur et al.²³ experimentally investigated the effects of changing the endplates from insulating boundary condition (IBC) to the conducting boundary condition (CBC) on the turbulence induced zonal flow. The experiment was conducted by changing the type of endplates while keeping other input parameters such as input gas pressure, magnetic field, and source power constant. Equilibrium density, equilibrium electron temperature, and mean plasma potential profiles were measured through RF compensated swept Langmuir probes. In both cases, similar equilibrium profiles of density and electron temperature profiles were measured within error bars (see Figs. 1(a) and 1(b)). Moreover, the CBC case potential profile is dominated by electron temperature profile when the endplate is grounded. However, in the IBC case, the potential of the endplate floats to whatever value necessary to force the parallel current to vanish, resulting in much smaller mean plasma potential gradients (see Fig. 1(c)). The Ion temperature was measured using Laser Induced Fluorescence (LIF) technique^24,25 and the neutral gas temperature was estimated through high resolution (⁠ $10 m Å$⁠) spectroscopy. The typical values of plasma parameters in the core of CSDX for these experiments are shown in Table I.

Measurements on CSDX show that a change of endplate from CBC to IBC can affect the turbulence structure, resulting in broadband turbulence in the IBC case and coherent modes in the CBC case (Fig. 2). The parallel current induced through sheath driven instability influences the perpendicular momentum balance of conducting case. Fig. 3 shows the difference in density fluctuations, floating potential fluctuations, and azimuthal flow between IBC and CBC experiments.

Simulations of insulating versus conducting case in this paper have been done with the same equilibrium fit (black lines in Fig. 1) to equilibrium profiles and plasma parameters (taken from CSDX 1 kG CBC experiment) in order to focus our studies of change in the turbulence due to the changes of the boundary condition. In this paper, we introduce our model and explore the sheath loss effects between the two simulation results.

As discussed in Sec. II, previous experiments suggest that the main instability active in CSDX is the resistive drift-wave. There are four main ratios in CSDX which determine the assumptions we can make in the derivation of a model of this experiment. Since measurements show that the ion temperature is relatively small compared to the electron temperature (⁠ $T i / T e ≈ 0.1$⁠), we will assume cold ions with $T ̃ i = 0$⁠. The second ratio of interest is $ρ s / L n ≈ 0.5$ which shows that radial global effects cannot be neglected; retaining radially varying profiles (nonlocality) is necessary. The third ratio is $λ m f p , e / L z ≈ 0.05$ which indicates collisional plasma behavior; thus fluid model can be used to describe the turbulence. For simplicity, we have neglected axial equilibrium gradients, since $k ∥ / k ⊥ ≈ 0.01$⁠. Here, L_n is the typical lengthscale of the density gradient, L_z is the approximate axial lenghtscale, and $λ m f p , e$ is the mean free path of electrons against momentum-exchange collision with ions and neutrals.

Here, we present the numerical model of drift-reduced nonlocal cold ion fluid equations which evolves fluctuations of density $n ̃$⁠, vorticity $Ω ̃$ and electron temperature $T ̃ e$⁠, in the presence of a fixed density profile $n 0 ( r )$⁠, potential profile $ϕ 0 ( r )$⁠, and electron temperature profile, $T e 0 ( r )$⁠. The effect of the ion parallel momentum is neglected here and deferred to future studies. The set of normalized partial differential equations reads as

where $ϕ ̃$ and $j ̃ ∥$ are the normalized potential and parallel current fluctuations. Momentum inertia is neglected, $Ω ̃ = ∇ ⊥ · ( n 0 ∇ ⊥ ϕ ̃ )$ is the generalized vorticity fluctuations, $Ω 0 = ∂ r ( r n 0 ∂ r ϕ 0 ) / r$⁠, the bracket ${ f , g } = ∂ r f ∂ θ g / r − ∂ θ f ∂ r g / r$ indicates the Poisson bracket operator containing advection nonlinearities where a Bousinesq approximation has been applied to the vorticity bracket, the tilde superscript denotes first order fluctuation quantity, 0 subscript denotes a zeroth order equilibrium profile, and Z subscript denotes the axially and azimuthally averaged zonal component of a fluctuating field or quantity, e.g., $f ̃ Z = ∬ f ̃ d θ d z / ∬ d θ d z$⁠. The mean profiles are then the sum of the equilibrium and time averaged zonal profiles, e.g., $f m e a n = f 0 + ⟨ f ̃ Z ⟩$⁠. Because we do not know the exact deposition profiles of these sources, or other processes such as ionization and recombination rates, we prescribed the equilibrium density and temperature profiles in our simulation, fitted to match experimental equilibrium measurements, and then damp the zonal density and temperature components to prevent quasilinear relaxation or flattening of the profiles. In a fixed-gradient approach, simulations use ad-hoc density and temperature sources (⁠ $− ν s n ̃ Z$ and $− ν s T ̃ e Z$ terms) in order to keep the equilibrium profiles from relaxing away from their experimental shapes. These terms damp zonal average of density and electron temperature at each time step. The CSDX helicon source acts as plasma density and heat source, without acting as momentum source. Hence, the zonal potential (⁠ $ϕ ̃ Z$⁠) is allowed to evolve in the simulation, allowing zonal flows to arise. ν_s is damping rate of zonal component of density (⁠ $n ̃ Z$⁠) and electron temperature (⁠ $T ̃ e Z$⁠) fluctuations.

Moreover, we refer to Eq. (4) in this paper as the adiabatic response of the system, which comes from electron parallel force balance equation between the electron pressure force and the electrostatic force while neglecting electron inertia.

D_n and $D T e$ are the normalized basic electron collisional diffusivity coefficients of density and electron temperature approximated through

is the normalized ion-ion viscosity calculated from a generalized Braginskii formula²⁶ instead of a simple magnetized limit viscosity, in which $X = Ω c i τ i$ is the magnetization factor. The normalized ion-neutral collision frequency is given as $ν i n = 6 × 10 − 19 n n T i / a e 0 m i$⁠, with n_n being the neutrals density, $χ ∥ e = 2.11 m i c s τ e / a m e$ is the normalized parallel electron heat conductivity, and $η ∥ = 0.51 a m e / m i c s τ e$ is the normalized parallel resistivity of plasma. The reference value for the density is $n ¯ = 9.9 × 10 12 cm − 3$ at r = 0, and the reference temperature is $T ¯ e = 3.1 eV$ at r = 0. The perpendicular spatial scales are normalized to gyroradius ρ_s, axial direction is normalized to radius $a = 8 c m$⁠, and timescales are normalized to $a / c s$⁠, where c_s is ion sound velocity at the axis.

The simulation has been performed in cylindrical coordinates, using an annular domain with infinitesimal inner radius (1% of outer radius) to avoid the numerical singularity at the axis. Dirichlet boundary conditions (⁠ $f ̃ = 0$⁠) have been used for evolving fluctuation properties at the outer radial boundary, while Neumann boundary conditions (⁠ $∂ f ̃ / ∂ r = 0$⁠) have been applied to evolve fluctuations at the inner radial boundary. The axial boundary condition for insulating or conducting endplate is non-trivial. For a perfect insulating endplate, the fluctuations of potential, electron temperature, and the parallel current at the sheath vanish. Moreover, since the insulating endplate can develop arbitrary wall potential value to set $j ∥ = 0$⁠, we have assumed no equilibrium plasma potential (⁠ $ϕ 0 = 0$⁠), and allow the $ϕ ̃ Z$ to be self-generated by the turbulence. In contrast, for a perfect electrical and heat conducting endplate, we used Bohm-like boundary condition,²⁷ where fluctuations are set to have zero axial gradient and parallel current fluctuations go to $j ̃ ∥ = ± ( ϕ ̃ − Λ T ̃ e )$⁠, while equilibrium potential $ϕ 0$ mimics equilibrium electron temperature $T e 0$ as long as the conducting endplates are unbiased such that $j ∥ 0 = 0$⁠. Hence, we are imposing $ϕ 0 = Λ T e 0$ throughout the domain in our model as we considered negligible axial equilibrium gradients and then allow the $ϕ ̃ Z$ to be self-generated. A more detailed derivation of these axial boundary conditions is given in  Appendix A.

In these CSDX experiments, the perpendicular length scales are comparable to the gyroradius. Therefore, it is important to retain the radial variation of profiles. Equilibrium profiles were determined by fitting a general functional form $n 0 = ( n ¯ − n b a s e ) exp ( − ( r / l n ) n α / n α ) + n b a s e$ using the least square method to the measurements taken from experimental conditions.²³ The general function fits to the experimental equilibrium profiles are shown in Figs. 1(a) and 1(b). Plasma parameters, such as the parallel electron heat conductivity and resistivity, then vary radially in a self-consistent manner according to equilibrium density and temperature profile fits, shown in Figs. 4(a) and 4(b). However, the ion-ion viscosity is taken to be radially constant for the simplicity of numerical implementation (see red dashed line in Fig. 4(c)). Our calculations show that the radial standard deviation of viscosity away from the mean is less than 20%. To estimate ion-neutral drag in the experiment, by knowing the injected gas pressure at the wall and assuming ideal gas behavior, we can approximate the neutral gas density from neutral gas temperature. Although in general, the ion-neutral drag can radially be varying, in this paper, we have assumed a radially constant drag $ν i n = 0.17 c s / a$⁠, and the effects of a radially varying drag profile are deferred to future studies.

These equations are solved in the BOUT++ Boundary Turbulence plasma fluid framework.²⁸ The derivatives, except for Poisson bracket terms, in the radial and axial directions are performed with a standard fourth order central finite-difference scheme, and in the azimuthal direction, a Fast Fourier Transform is used. The PVODE solver²⁹ is used for time integration, which combines a Backward Difference Formula with a preconditioned Generalized Minimal Residual (GMRES) method. For the nonlinear advection terms, the Arakawa discretization³⁰ is implemented for the Poisson brackets.

To ensure that the code is solving the model equations correctly, we have verified the code in both the linear limit and nonlinearly saturated regime. First, we note Dudson et al.³¹ have verified many of the differential operators we use in BOUT++ framework through method of manufactured solutions (MMS). Second, in the linear local limit, we have verified the linear eigenmode growth rates and frequencies against the analytical dispersion relation separately for drift-waves, Kelvin-Helmholtz (KH), and sheath-driven temperature gradient instability. Moreover, the energy balance of the code is investigated in Secs. IV–V and shows energy conservation with less than 1% error.

In this section, we discuss the linear physics of conducting versus insulating endplate simulations. Linearly, the IBC case allows only drift-waves to grow and propagate within the plasma. The dominant instability in the IBC case is the density gradient-driven drift-waves.³² Electron temperature gradient-driven drift-waves (thermal drift-waves) are also present, but are less prominent and usually suppressed via parallel conduction. However, in the CBC case, since the equilibrium plasma potential is locked to the equilibrium electron temperature by the endplate sheath physics, sheath-driven linear transverse Kelvin Helmholtz (KH) instability due to the $− a ρ s 1 r ∂ Ω 0 ∂ r ∂ ϕ ̃ ∂ θ + a ρ s 1 r ∂ ϕ 0 ∂ r ∂ Ω ̃ ∂ θ$ term, rotational (centrifugal) modes due to the $a ρ s 1 r ∂ ϕ 0 ∂ r ( ∂ n 0 ∂ r ∂ 2 ϕ ̃ ∂ r ∂ θ − ∂ 2 ϕ 0 ∂ r 2 ∂ n ̃ ∂ θ )$ term, as well as the sheath-driven electron temperature gradient instability due to non-zero parallel current $j ̃ ∥ = ± ( ϕ ̃ − Λ T ̃ e )$ at the sheath boundary, appear in addition to the drift-waves.³³ 

In Fig. 5, we show a comparison of the linear growth rate and real frequency versus azimuthal mode number for the local pure drift-wave in the IBC case, the nonlocal pure drift-wave in the IBC case, and nonlocal growth rates of the CBC case. For the nonlocal linear runs, we have solved the Eqs. (1)–(4) without nonlinear advection terms, and the equilibrium profiles are taken from Fig. 1. However, in the local limit, the equilibrium density and temperature are taken to be radially constant and equal to their on-axis values, and the equilibrium density and temperature gradients taken to be independent parameters which are also radially uniform, and solving for Eqs. (1)–(4) in the linear limit. The linear growth rates (frequencies) of the eigenmodes have been obtained by calculating imaginary (real) part of $∂ ln ϕ ̃ m / ∂ t$ after convergence in our initial value solver, where m denotes the azimuthal mode number. We observe that inclusion of the global radial profile effects decreases the growth rate, and shifts the linear frequencies, showing that even in the linear limit retaining the radial profile of input parameters is important. In the CBC case, various types of instabilities (density gradient-driven drift-waves, sheath-driven linear KH and centrifugal force, and sheath-driven electron temperature gradient instability) compete with each other. Consequently, the linear eigenmodes in the CBC have smaller growth rates as compared to the drift-wave growth rates found in the IBC case.

To understand the mechanism behind these different linear behaviors, we have used an energy balance analysis. The energy balance equation can be obtained by multiplying Eqs. (1)–(3) by $T e 0 n ̃ / n 0 , − ϕ ̃$⁠, and $3 n 0 T ̃ e / 2 T e 0$⁠, respectively, so that sum of energy equations conserve the adiabatic response. The total energy balance equation, reads as

The equation is normalized to the instantaneous total energy of system $E t o t = 1 2 ( ⟨ T e 0 n 0 n ̃ 2 ⟩ + ⟨ n 0 ( ∇ ⊥ ϕ ̃ ) 2 ⟩ + 3 2 ⟨ n 0 T e 0 T ̃ e 2 ⟩ )$⁠, where $S n = − a ρ s ⟨ T e 0 ∂ r ln ( n 0 ) n ̃ v ̃ r ⟩ / E t o t$ is the extraction of energy from the equilibrium density gradient, $S T e = − a ρ s ⟨ 3 2 n 0 ∂ r ln ( T e 0 ) T ̃ e v ̃ r ⟩ / E t o t$ is the linear energy extraction from the equilibrium electron temperature gradient, $S K H = a ρ s ⟨ 1 r ∂ r ϕ 0 ϕ ̃ ∂ θ Ω ̃ ⟩ / E t o t$ is the KH energy extraction, and $S r o t = a ρ s ⟨ 1 r ϕ ̃ ∂ r ϕ 0 ( ∂ r n 0 ∂ r θ 2 ϕ ̃ − ∂ r r 2 ϕ 0 ∂ θ n ̃ ) ⟩ / E t o t$ is the energy extraction through rotational (centrifugal) modes. $S j ∥ = − ⟨ η ∥ n 0 j ̃ ∥ 2 ⟩ + ⟨ ∇ ∥ ( T e 0 n 0 n ̃ j ̃ ∥ + 1.71 T ̃ e j ̃ ∥ − ϕ ̃ j ̃ ∥ ) ⟩ / E t o t$ is the parallel current energy sink in which the first term is parallel resistive dissipation and the second term is due to sheath losses, $S ⊥ = − μ ⊥ ⟨ ϕ ̃ ∇ ⊥ 2 Ω ̃ ⟩ / E t o t + D n ⟨ T e 0 n 0 n ̃ ∇ ⊥ 2 n ̃ ⟩ / E t o t + D T e ⟨ 3 n 0 2 T e 0 T ̃ e ∇ ⊥ 2 T ̃ e ⟩ / E t o t$ is the perpendicular energy dissipation rate, $S d r a g = ⟨ ν i n ϕ ̃ Ω ̃ ⟩ / E t o t$ is the ion-neutral drag energy dissipation rate, $S χ ∥ e = ⟨ 3 n 0 2 T e 0 χ ∥ e T ̃ e ∇ ∥ 2 T ̃ e ⟩ / E t o t$ is the parallel conductivity energy dissipation rate, $S s o u r c e = ν s ( ⟨ T e 0 n 0 n ̃ Z 2 ⟩ + ⟨ 3 n 0 2 T e 0 T ̃ e Z 2 ⟩ ) / E t o t$ is the ad-hoc source energy rate, and $S a d v =⟨ T e 0 n 0 n ̃ v ̃ · ∇ ⊥ n ̃ ⟩/ E t o t +⟨− ϕ ̃ v ̃ · ∇ ⊥ Ω ̃ ⟩ / E t o t + ⟨ 3 n 0 2 T e 0 T ̃ e v ̃ · ∇ ⊥ T ̃ e ⟩ / E t o t$ is the energy rate of nonlinear advection terms. The averaging operator $⟨ f ⟩ = ∫ d z ∫ r d r ∑ m f$ denotes a volume average operator, where m is the azimuthal mode number.

The extra factors n₀ and $T e 0$ are added to the energy terms in order to show energy conservation of adiabatic response in Eq. (4). While the physical energy contains extra factors of $n 0$ and $T e 0$⁠, the physical energy does not preserve the property of conservative nonlinearities in Eqs. (1)–(3), and therefore produces a more complicated analysis.¹² Therefore, since we are using a nonlocal model and $∇ ⊥ · T e 0 n 0$ is non-zero, we have decided to show conservation of the adiabatic response in our formulation for better understanding of parallel energy conservation. We should also note that the terms with $⟨ f ∂ θ f ⟩$ yield zero volume averaged energy due to the azimuthal symmetry, and hence, they are not shown in our calculations.

The balance between the right and left hand sides of Eq. (5) determines the conservation of energy in the system, and evaluation of the difference between these two quantifies the accuracy of the simulation in conserving energy. Fig. 6 shows the time traces of the normalized energy rates in the system for the IBC and CBC linear runs in experimental condition. The dashed lines show the time traces of the energy error which are much smaller than other physical terms in the system (less than 1% of total energy rate of system), demonstrating very good conservation of energy in both simulations. In the linear limit, the normalized energy rates asymptote to a steady value once the dominant modes establish themselves, and effects of the initial conditions dissipate away (or at least become negligible). We should note in the linear limit, S_adv and S_source are zero, and thus are not shown in Fig. 6.

To understand the physics of the instabilities, we performed a density profile scan below and above the experimental condition by changing the exponential exponent parameter $n α$ in our functional form for $n 0 ( r )$⁠, causing a decrease or increase in the density gradient (⁠ $1 / L n = − d ln n 0 / d r$⁠). The results of these scans are shown in Fig. 7. We can observe that in the IBC case (Fig. 7(a)), $S T e$ is small because the maximum electron temperature gradient (⁠ $1 / L T e = − d ln T e 0 / d r$⁠) is located in the region of low density, and also due to the presence of strong parallel damping of electron temperature fluctuations. The drift-waves are primarily driven by the density gradient as seen by the dominant S_n term. Therefore, we can deduce that the electron temperature gradient drift-wave is not a main contributor of instability to the CSDX IBC case, even with very small density gradients profiles.

In the CBC case (Fig. 7(b)), we can observe some notable changes compared to the IBC case. First, the presence of conducting sheath allows additional linear instabilities such as sheath-driven electron temperature gradient instability, plus linear transverse KH and rotational (centrifugal) instability generated by the sheath-driven equilibrium E × B shear. Changes in the density profile then leads to different dominant linear instabilities regimes, due to variations in the electron diamagnetic frequency (⁠ $ω * e = k θ ρ s c s 1 L n$⁠) and parallel electrical conductivity (⁠ $σ ∥ = 1 / η ∥$⁠). Jassby³⁴ demonstrated different regimes of drift-waves and KH instability, depending on electron diamagnetic frequency and parallel electrical conductivity in $k ⊥ ≫ k ∥$ limit. In Fig. 8, we show the contour plots of linear growth rates and the density-potential cross-phase $ζ ( n ̃ , ϕ ̃ ) = arg ( n ̃ m * ϕ ̃ m )$ for the CBC linear density gradient scan study. Through this plot, we can designate which mode is unstable for a given set of parameters. We should note that due to the nonlocal nature of our model, the radial structure of the cross-phase is not radially uniform and different instabilities may be dominant locally at certain radii; however, this simple analysis can give us insight about the most dominant general instability regime via radially averaged cross-phase values.

For the CBC case, at a low density gradient, where $max ( 1 / L T e ) > max ( 1 / L n )$⁠, we can observe in Figs. 7(b) and 8 that sheath imposed linear KH (S_KH) and sheath driven electron temperature gradient modes (⁠ $S T e$⁠) act as main driving sources of instability, while density-gradient drift-waves and centrifugal instabilities are present but less prominent. The presence of $S T e$ at low density gradients is mainly because of sheath driven electron temperature modes, and smaller parallel conduction of electron temperature fluctuations (smaller $χ ∥ e$⁠). Here, most of density-gradient driven drift-waves are strongly damped through parallel current dissipation $S j ∥$⁠. From Fig. 8 when $1 / L n$ is small, we can observe $5 < m < 8$ modes are unstable and have average density-potential cross-phase $π / 2 < ζ ( n ̃ , ϕ ̃ ) < π$⁠, which is an indicator of dominant linear KH instability.³⁴ 

As $1 / L n$ is increased to approach the experimental value, we observe that the S_n becomes the most dominant instability, while S_KH is still present and contributing to the linear drive of the system. The sheath-driven electron temperature gradient modes (⁠ $S T e$⁠) become smaller and eventually become damped due the nonlinear dependence of the parallel heat conductivity $χ ∥ e$ on density. In Fig. 8, we can see the formation of a new set of modes (⁠ $1 < m < 4$⁠) as $1 / L n$ increased, which have $ζ ( n ̃ , ϕ ̃ ) < π / 4$ density-potential cross-phase, indicating the dominance of drift-waves. The CSDX CBC experimental condition lies within a hybrid resistive drift-wave/KH regime, where the flow driven modes are opposed to drift-wave, and may be expected to have different turbulence characteristics compared to IBC case. At even higher density gradients, S_n becomes the main drive of instabilities and S_KH contribution becomes smaller, and the system becomes dominated by drift-waves. The CBC density profile scan within different regimes of instabilities is in qualitative agreement with the picture depicted by Jassby.³⁴ 

Building upon the linear analysis of Sec. IV, nonlinear simulations of the IBC and CBC cases were performed. The simulations include a growing linear phase followed by the development of a saturated turbulence regime. Fig. 9 shows the time trace of the volume averaged total energy for both insulating and conducting simulation. Both simulations were run with the same fit to equilibrium profiles shown in Figs. 1(a) and 1(b). We observe that the average total energy of the system in the CBC simulation is smaller than in the IBC case and saturates more slowly, consistent with expectations from linear analysis which gave a smaller γ_lin in the CBC case, while the IBC exhibits a more pronounced time-variation in the total energy $E t o t ( t )$⁠.

The damping rate of the zonal density and zonal electron temperature fluctuations was chosen to be large enough that the time averaged zonal fluctuations remain very small relative to the input equilibrium profiles. In  Appendix B, we present the results demonstrating convergence at the simulation saturation statistics with a large enough value of ν_s.

Snapshots of the $r − θ$ distribution of density fluctuations filtered for finite wave numbers in the two cases are shown in Fig. 10. We can statistically observe that $m ≈ 2$ modes dominate both simulations. In the IBC case, density blobs are mostly being formed in the region of the maximum density gradient, propagate in the electron diamagnetic direction, and they exhibit a radial extension that becomes tilted azimuthally as the structure approaches the shear layer. Similar qualitative behavior has been reported experimentally,²³ as well.

Moreover, in Fig. 11, we have shown the snapshots of the r – z distribution of density fluctuations filtered for finite wave numbers for the IBC and CBC simulation. We observe the dominant axial mode is $n ≈ 1$ for both cases, indicating the dominance of drift-waves instability. However, the shape of eigenfunctions is changing from IBC to CBC case, due to the effects of the axial boundary change. We should note that axial equilibrium gradients are not included in the current model and its effects may change the shape of eigenfunctions. Currently, there is no axial measurement of fluctuations in the experiment, and further study on this subject is needed.

An energy balance analysis (analogous to what was introduced in the Sec. IV) was performed for nonlinear simulations as an additional verification test. As before, the energy error obtained is much smaller than other physical terms in the system demonstrating energy conservation in the nonlinearly saturated regime. Details of convergence of the energy rates with grid resolution is shown in  Appendix C. Fig. 12 shows the mean averaged values of energy rates in the saturated turbulent regime. We observe a smaller value of the density-gradient drive (S_n) in the CBC case than the IBC case. Moreover, as shown in Sec. IV, S_KH and S_RIC are present as other mechanisms of turbulence drive in the CBC simulation; however, they are smaller compared to S_n similar to the linear energy analysis consistent with the drift-turbulence dominated system. One notable difference we observe between the linear energy rates and the nonlinear saturation energy rates is the larger energy dissipation rates in the nonlinear simulation, mainly due to the damping of the turbulence-generated m = 0 plasma potential energy by the ion-ion viscosity and ion-neutral drag.

To further understand saturation mechanism and drift-wave/zonal flow paradigm in the nonlinear simulations, we investigate the spectral energy balance of the vorticity equation. The spectral energy balance of the vorticity equation (Eq. (2)) for the azimuthal mode m reads as

where the star superscript denotes the complex conjugate, and the averaging operator is defined as $⟨ f ⟩ = ℜ ( ∫ d z ∫ r d r f ¯ ) t$⁠, and time averaged variation of kinetic energy $∂ t ⟨ E k i n m ⟩$ in the saturation regime goes to zero. In addition, the $⟨ ϕ ̃ m * ∇ ∥ j ̃ ∥ m ⟩$ term can be rewritten as a sum of the injected energy into potential fluctuations from the adiabatic response $⟨ j ̃ ∥ m * ∇ ∥ ϕ ̃ m ⟩$ plus the potential sheath boundary loss term $− ⟨ ∇ ∥ ( ϕ ̃ m j ̃ ∥ m * ) ⟩$⁠. Thus, we can rewrite the generalized spectral kinetic energy balance as

In Fig. 13, we plot the time the averaged spectral energy rates of the IBC and CBC simulation in the saturated regime, normalized to the total kinetic energy. We observe the IBC simulation has more active modes (⁠ $2 ≤ m ≤ 8$⁠) contributing to the energy transfer into the m = 0 zonal flow. Although in the CBC simulation the energy injection from the adiabatic response is smaller than in the IBC case, the CBC simulation exhibits additional KH and centrifugal energy injection into finite m fluctuations. In addition, the potential sheath boundary loss term is an important energy sink in the CBC simulation, resulting in significant damping of kinetic energy. Consequently, due to the larger sheath boundary loss of vorticity energy balance in the CBC case, the nonlinear kinetic energy advection term $T u = − a ρ s ⟨ ϕ ̃ m * { ϕ ̃ , Ω ̃ } m ⟩$ is weaker in transferring energy from the finite m modes to the m = 0 zonal potential, resulting in weaker inverse cascade dynamics for the CBC case. Thus, we observe the zonal potential energy to be smaller in the conducting case. As a result, the IBC simulation exhibits more prominent 2D turbulence dynamics, in which the most linearly unstable modes have a negative T_u, and thus losing kinetic energy, while the higher order modes (m > 6) and the zonal mode (m = 0) have $T u > 0$ and are receiving kinetic energy from the intermediate fluctuations. These observations of the inverse energy transfer dynamics are also qualitatively consistent with the previous experimental measurement of frequency space kinetic energy transfer in various CSDX experiments.^20,22,23

The saturated amplitude profiles of the fluctuations for both the cases are shown in Fig. 14. We observe that the plasma potential and density fluctuations have the largest fluctuations amplitudes, and the electron temperature fluctuations are small compared to $n ̃$ and $ϕ ̃$ (see Figs. 14(a)–14(c)). The RMS density peaks near $r ≈ 2.5 cm$⁠, consistent with the experimental measurement of density (see Fig. 3(a)). However, the plasma potential RMS amplitudes peak more toward the edge (⁠ $r ≈ 4 cm$ for CBC simulation and $r ≈ 6 cm$ for IBC simulation), while the experimental floating potential measurements peak at $r ≈ 3 cm$ (see Fig. 1(b)). Additional simulations not shown here indicate that the location of the peak RMS plasma potential fluctuation is sensitive to the details of the equilibrium profile as well as to the radial profile of ion-neutral drag.

Fig. 14(d) shows total mean plasma potential profile (the sum of the equilibrium potential $ϕ 0$ and zonal potential $ϕ ̃ Z$⁠). We observe that the IBC case exhibits a completely different behavior compared to the CBC case. The generated zonal electric field in IBC simulation is radially inward, causing an E × B drift in the electron diamagnetic direction. However in the CBC simulation, the mean potential profile is mostly dominated by sheath-driven equilibrium electron temperature with an outward electric field, while the CBC turbulence-driven component of the zonal potential gives a radially outward electric field in the core and radially inward electric field near the edge, qualitatively consistent with the experimental observations (see Fig. 1(c)).

Figs. 14(e) and 14(f) show smaller Reynolds stress and zonal flow in the CBC case, indicating that the conducting endplate partially short-circuits the Reynolds stress. This phenomenon can be explained with the charge conservation equation $∇ · j ̃ = ∇ ⊥ · j ̃ ⊥ + ∇ ∥ j ̃ ∥ = 0$⁠. Surface integrating current conservation over a flux tube in the plasma, we can write for radial current at radius r as

Since in the CBC simulation, the parallel current is non-zero at the sheath boundary $j ̃ ∥ ( z = L z ) = ϕ ̃ − Λ T ̃ e$⁠, there is finite radial current opposing the polarization drift. However, in the IBC simulation, this radial current vanishes due to the zero parallel current at the sheath boundary. As a result, the turbulence inverse energy cascade is weaker in the case of CBC simulation. This simple explanation is also consistent with our spectral kinetic energy balance analysis shown in Fig. 13. Experimental measurements²³ also observe similar qualitative reduction in Reynolds stress and zonal flow for conducting endplate case.

Unlike D'Ippolito et al.¹⁵ simulation which had a dipolar azimuthal velocity with zero net azimuthal momentum, our IBC simulation (Fig. 14(f)) shows non-symmetric lobes of $v E × B$ due to the presence of ion-neutral drag and nonlocality of vorticity fluctuations in the cylindrical geometry. The CBC simulation shows two weaker lobes of $v E × B$⁠, which is also non-symmetric. Yan et al.³⁵ discusses the experimental evidence of net plasma rotation in the absence of momentum input due to the ion-neutral damping in the outer region of plasma. Here, in the simulations, even a radially uniform ion-neutral drag profile creates residual stress acting on the net plasma rotation. However, the radial location of shear is different from the experimental measurements. The generation of the zonal flow is a highly nonlinear process, and many factors and input profile uncertainties can affect its amplitude and radial structure.

Fig. 15 shows the potential fluctuations autopower spectra in frequency and azimuthal wavenumber space. The IBC case shows more broadband turbulence, where the energy is mostly accumulated in low frequencies consistent with the spectral transfer results, while the CBC case exhibits a more quasi-coherent narrow band spectra. In both cases, the power spectra peak not far from the Doppler shifted linear frequencies (shown as the white curve in Figs. 15(c) and 15(d) calculated from nonlocal linear frequencies (Fig. 5(b)), which shows that the linear instabilities continue to play an important role for setting the turbulence dispersion relation. We can clearly observe coherent features in the m – f contour of the CBC simulation. One of the main differences observed in qualitative comparison with experiment is that the potential fluctuation amplitude is peaked more radially outward near the edge of the simulations. Simulation plasma potential power spectral density contours (Figs. 15(a) and 15(b)) qualitatively look like the experimental floating potential autopower contours (Fig. 2). A direct quantitative simulation-experiment comparison requires a detailed validation study and use of synthetic diagnostics, and will be reported in a separate publication.

In this paper, we have illuminated how changing the parallel boundary conditions have an effect on the development of electrostatic drift turbulence and sheared E × B flows in a linear plasma column. Linearly, our analysis showed that with insulating endplates, the dominant instability is the density gradient-driven collisional drift-wave. With conducting endplates, the plasma lies in a hybrid drift-wave/KH regime, where the density gradient-driven drift-wave is still the dominant instability, but KH instability generated by the sheath-driven equilibrium E × B shear is also an important player and integral to the model. Both sheath-driven electron temperature gradient and rotational (centrifugal) force are feeble and compete against density gradient drift-waves. Presumably, in CBC case, by decreasing the density gradient, the system would go to KH dominated regime.

Moreover, we can conclude that endplates can have significant effect on vorticity equation balance. Conducting endplate allows the current and the heat flux to escape the device. These sheath losses dissipate energy from the system and weakens the Reynolds stress. In general, the endplate change affects the momentum equation, polarization current, vorticity flux, and shear layer generation. These results are qualitatively consistent with the experimental report of Thakur et al.²³ and previous 2D slab simulations of D'Ippolito et al.¹⁵ In addition, we have shown that in the CSDX insulating case, the electron temperature does not significantly alter the nature of instabilities, while in the CSDX conducting case, the electron temperature can be an essential part of the physics by imposing an equilibrium potential profile on the plasma that can counteract the turbulence driven zonal flow. However, in both cases, electron temperature fluctuations may become important in validation studies and should be neglected only with considerable deliberation.

These simulations provide a basis for validation studies of turbulence and transport in CSDX linear device. Experimental studies on CSDX have used saturation current and floating potential measured by Langmuir probes, and light intensity fluctuations measured by fast imaging camera, to interpret the turbulent transport. Many of these measurement techniques are also being used in larger magnetically confined devices, where these measurements might differ from real density and the plasma potential values because of the effects of electron temperature fluctuations. These effects have been neglected in CSDX experimental studies, due to diagnostic limitations. Our model allows us to translate back the real physical quantities into synthetically measured quantities and compare them to experimental measurements. A future paper will provide a detailed validation study and will discuss the effects of profile uncertainty on the degree of agreement between simulation and experiment.

This research was supported by the U.S. Department of Energy Grant No. DE-FG02-06ER54871. We thank B. D. Dudson, B. Friedman, X. Q. Xu, J. Loizu, P. H. Diamond, and A. Ashourvan for the useful discussions.

Here, we present the derivation of axial boundary conditions of the sheath entrance for the conducting and the insulating endplates. When the magnetic field is perpendicular to an absorbing wall, there is no magnetic presheath and the plasma-wall transition will consist of two regions between main bulk plasma and the wall: first, a Debye sheath which is in contact with the wall, and a collisional presheath in between the Debye sheath and unaffected plasma.²⁷ The Debye sheath is where the recombination with the wall occurs, has a size proportional to Debye Length (λ_D), and both quasi-neutrality and the ion-drift approximation break down in this region. The collisional presheath size usually scales with λ_mfp and is quasi-neutral, and in it, the ions are magnetized and accelerated towards the wall, approximately reaching the plasma sound speed c_s at the sheath entrance.

Since our fluid model is quasi-neutral and uses the ion-drift approximation, we need to find the boundary values at the sheath entrance. In order to describe steady-state dynamics of collisional presheath, we will follow the procedure used by Loizu et al.³⁶ but for a simple perpendicular magnetic field to derive a rigorous system of equations, including ion continuity, ion and electron parallel momentum, and electron temperature equations valid in the collisional presheath region.

Specifically, the logic used here in the derivations is analogous to the derivation of conducting magnetic pre-sheath boundary condition done in Loizu et al.,³⁶ but in the limit of the magnetic field exactly perpendicular to the conducting endplate. We start by deriving the boundary condition for a conducting endplate, where the electrical and heat conductivity of sheath is very large.

The steady-state continuity equation for ions is

where $S p i$ contains other ionization/recombination processes. We can rewrite the equation as

The fourth and fifth terms cancel out each other in the above equation. Since equilibrium gradients at the radial boundaries are vanishing in our experiment (see Fig. 1(d)), we can omit the third and sixth terms in the above equation up to first order expansion; therefore, we can further simplify the equation

Now let us consider the parallel component of ion momentum equation

Thus, in the electrostatic limit the equation reduces to

Turning to the parallel component of electron momentum

by assuming negligible electron inertia and high conductivity, the equation reduces to

Inclusion of non-isothermal electrons requires the use of the electron heat equation. By assuming very large electron heat conductivity at the sheath entrance, and neglecting inertia and diffusion, and omitting the transverse terms, the electron heat equation reads as

Thus

We can rewrite the system of equations (Eqs. (A2), (A4), (A5), and (A6)) in vector form $M X = S$

Note that all variables are total values (zeroth order plus fluctuations). Assuming source/sink terms do not affect axial gradients in the presheath region (i.e., assuming $| M i j X j | ≫ | S j |$⁠), we can solve for $det ( M ) = 0$⁠, obtaining the value of ion parallel velocity at sheath entrance as

To obtain the parallel electron velocity at the sheath entrance, we use the Boltzmann relation and ambipolarity of electron flux between sheath entrance and the wall, and we assume $λ D ≫ ρ e$ such that the electrons stay magnetized all the way to the wall. This relation is not strongly true for the CSDX parameters, but for the simplicity we will use this assumption.

Assuming Maxwellian electrons in the sheath, we can write Boltzmann relation for electrons between sheath entrance and wall

where n_e and $ϕ$ are the electron density and potential at the sheath entrance, and $n e w$ and $ϕ w$ are the electron density and potential at the wall, respectively. By assuming grounded wall (⁠ $ϕ w = 0$⁠) and also assuming electrons hit the wall with their thermal velocity, the parallel electron velocity can be obtained through flux ambipolarity between the sheath entrance and wall as

in which $Λ = log m e / 2 π m i$ is the sheath factor, and for Argon gas $Λ A r = 5.1$⁠. Therefore, we can obtain closure terms needed in Eq. (A8) as

Using these closures for Argon gas, the ion parallel velocity (Eq. (A8)) reduces to

Hence, the axial electron temperature gradient and potential gradient at the sheath edge are given as