Two-fluid sub-grid-scale viscosity in nonlinear simulation of ballooning modes in a heliotron device

Here we assess the applicability of the LES approach to the study of mild saturation by revisiting our earlier simulations of the single-fluid MHD equations reported in [7]. The two-fluid term, the gyro-viscosity and the binormal parallel heat conduction are missing from the system of equations (1)–(7). Nevertheless, it is still worth revisiting the numerical data in order to evaluate whether the introduction of the LES approach is appropriate or not.

The initial force-balance equilibrium of ${{\beta}_{0}}=3.6 \%$ and the magnetic axis position ${{R}_{\text{ax}}}=3.6$ m shown in [7] are used again for the LESs in this paper, for the purpose of comparing our new LES results with the earlier studies. This initial force-balance state is computed by the use of the HINT code [19] for an ideal MHD, namely one excluding all the dissipative coefficients in equations (1)–(7). Neither these equilibrium components nor the fluctuation components are separated from each other, but consist of the variables solved in (1)–(7) together. Though the HINT code can construct an equilibrium with magnetic islands, the equilibrium used here does not initially have a magnetic island. It is also noted that the HINT solution can serve as the equilibrium state of the extended MHD equations (1)–(7) in the limit of the small η and ${{\kappa}_{\bot}}$, because ρ is initially uniform in this equilibrium. In figure 1(a), the radial profile of the poloidal (m) and toroidal (n) Fourier coefficient m/n  =  0/0 of the pressure projected on the Boozer coordinate is shown. In the figure, the symbol ψ represents the initial flux function. Figure 1(b) is the rotational transform of the equilibrium. The rational surfaces around $\iota /2\pi \simeq 0.5$ are most unstable, so we expect the growth of unstable modes with m/n  =  2/1, 4/2 and so on. Such a configuration is similar to the experimental one reported by Sakakibara et al [17], in which it is reported that the β-value can be raised after removing the resonant magnetic surface of $\iota /2\pi =0.5$, and that the MHD activity does not bring about the total collapse of the pressure. While a pressure flattening at $\iota /2\pi =0.5$ is observed in the course of the experiment, the flattening disappears and a clear peaked profile is achieved after the removal of the resonant magnetic surface. All the simulations in [7] and in this paper start from this HINT equilibrium, the magnetic configuration of which is similar to that in the experiment in [17]. In [7], many unstable ballooning modes grow simultaneously under the relatively weak dissipation of $\mu =\eta ={{\kappa}_{\bot}}=1\times {{10}^{-6}}$, but with a large parallel heat conductivity ${{\kappa}_{\parallel}}=1\times {{10}^{-2}}$. (Note that ${{\kappa}_{\bigwedge}}=0$ in the simulations.) These unstable modes are resonant to the rational surfaces around $\iota /2\pi \simeq 0.5$, and begin to interact with each other in the nonlinear stage. In the nonlinear course of the time evolution, the plasma motion is quite turbulent, causing it to mix in the core and peripheral regions. As a consequence of the growth, the plasma core pressure collapses, contrasting with the experimental result in [17].

Zoom In Zoom Out Reset image size

Since LES is a numerical method, which is based on the assumption that the flow is turbulent and the nonlinear energy transfer from a large scale to a small scale can be replaced by a Smagorinsky-type (an eddy-viscosity-type) model, we should first verify that the plasma dynamics is turbulent at a nonlinear stage of time evolution. In figure 2, the pressure spectrum at the time of the nonlinear saturation is shown as a function of the poloidal wave number m for a computation of ${{\kappa}_{\parallel}}=1\times {{10}^{-2}}$ and ${{\kappa}_{\bot}}=\mu =\eta =1\times {{10}^{-6}}$. The plot is scattered for the toroidal wavenumber n  <  20 so that we see 20 points for a single value of m. It is seen in figure 2 that the Fourier energies of the pressure are excited over a wide range of the wavenumbers, and show a tendency in which the Fourier energies may be scaled by the power of the wavenumber m, as we often see in the MHD and Hall MHD turbulence spectrum [15, 26]. In the figure, it is difficult to find reminiscent linear eigenfunctions once formed at the linear stage of the simulation. Some n  =  2 coefficients have the largest values of all ns. However, it looks fully turbulent when we see the pressure fluctuation on a poloidal cross-section, being different from those observed at the linear stage (the figure is omitted), due to the strong stretching, entrainment and mixing of the pressure on a poloidal cross-section by the vortices generated originally by the unstable ballooning modes. Thus, the dominance of the n  =  2 Fourier coefficients does not necessarily mean the dominance of the ballooning mode itself, but rather that of a structure which has survived the nonlinear stage.

Zoom In Zoom Out Reset image size

Second, we estimate how much the SGS terms can contribute to the mild saturation of unstable modes by making use of our earlier numerical data. While numerical results in [7] show a large collapse at the plasma core, the growth of unstable modes is mildly saturated for the same initial equilibrium when the viscosity is set higher ($\mu =5\times {{10}^{-4}}$ in [5]). Thus, estimating how a large ${{\mu}_{\text{SGS}}}$ can grow in a simulation in comparison to $\mu =5\times {{10}^{-4}}$ provides an insight for the magnitude of ${{\mu}_{\text{SGS}}}$. The time evolution of ${{\mu}_{\text{SGS}}}$ in equation (27) by using the numerical data in [7] is plotted in figure 3(a), and the isosurfaces of ${{\mu}_{\text{SGS}}}$ at the time of nonlinear saturation are shown in figure 3(b) for ${{C}_{\nu}}={{C}_{\eta}}=0.046$ according to [16] for incompressible channel MHD turbulence. (For Hall MHD turbulence, we consider that ${{C}_{\nu}}={{C}_{\eta}}=0.345$ can work better than ${{C}_{\nu}}={{C}_{\eta}}=0.046$ for reproduction of the spatial structures, although it is too diffusive to reproduce the energy spectrum [18].) With respect to the importance of the parameters in the simulations, $\mu =\kappa =\eta =1\times {{10}^{-6}}$ and the number of grid points $193\times 193\times 640$ in the reference simulation in [7] as well as $\mu =5\times {{10}^{-4}}$ and $\kappa =\eta =1\times {{10}^{-6}}$ have been set to the same resolution as in [5]. Note that the growth of the unstable modes in the latter computation is mildly saturated and that well-formed magnetic surfaces are recovered autonomously after saturation. The plot indicates that the truncated SGS components can play the role of a viscosity that is as high as ${{\nu}_{\text{SGS}}}$ to the GS (low-wavenumber) modes through nonlinear couplings. Since the maximum value of the SGS viscosity in figure 3(a) is higher than $\mu =5\times {{10}^{-4}}$, it is expected that the introduction of the SGS viscosity can lower the saturation level of the instability, although it is also considered that mild saturation of the instability cannot only be achieved by the introduction of the SGS effects, because the unstable low modes can grow linearly before the SGS effects begin to work. Since the estimation of the ${{\nu}_{\text{SGS}}}$ is for the grid points of $193\times 193\times 640$, the role of the SGS components on the GS components can be larger for a smaller number of grid points.

Zoom In Zoom Out Reset image size

Based on an assessment of the LES approach, we carry out the LESs of ballooning modes with a number of grid points $97\times 97\times 640$. The number of grid points is 1/4 of our earlier computations, and the time step width is also coarser. This coarse-graining in both time and space decreases the computational cost significantly. We do not carry out LESs with a smaller number of grid points (larger grid width), because our earlier studies on homogeneous Hall MHD turbulence [15, 27] show that very coarse low-pass filtering can bring about a considerable change in the nonlinear interaction of the Hall term at low wavenumbers. This simply means that solving the equations with a very coarse grid width cannot express the nonlinear dynamics of the low modes, even if we compensate for some aspects of the nonlinear interactions using the SGS models. As far as we have examined, the grid width should be as small as, or smaller than, the ion skin-depth in order to avoid significant change—unless we develop a new SGS model which can correct the change of the nonlinear interaction at low wavenumbers. (This is another way of mentioning the importance of turbulent high-wavenumber components to unstable low modes.) Thus we set the number of grid points $97\times 97\times 640$ in this paper for the Hall parameter ${{\epsilon}_{H}}=0.0125$ (the ion skin-depth is comparable to the grid width) or to ${{\epsilon}_{H}}=0.05$ (the ion skin-depth is larger than the grid width).

In figures 4 and 5, a typical example of an LES result is shown by use of the parameters ${{C}_{\nu}}={{C}_{\kappa}}={{C}_{\eta}}=0.046$ and ${{\epsilon}_{H}}=0.0125$ (run LES41). See table 1 for the typical LESs and parameters used for each of the runs. In figure 4(a), the isosurfaces of the pressure and the streamlines at the end of the linear stage of an LES, magnified around a poloidal cross-section, are shown. The streamlines show that flow is generated in the poloidal (binormal) direction. Since the poloidal flow does not appear in the single-fluid MHD simulations and appears as the result of the introduction of the two-fluid effect, it is considered to have been generated by the diamagnetic effect. Because of the flow in the poloidal direction, we can observe that the pressure profile is elongated in the poloidal direction, causing thermal mixing and homogenization in the direction of a fully nonlinear stage (the figure is omitted). In figure 4(b), the isosurfaces of the pressure, magnetic field lines (thick tubes) and streamlines (thin tubes) at a nonlinear stage are shown for the fully torus geometry. The streamlines are spirally elongated along the magnetic field lines in some regions, and drift in the poloidal direction, forming a turbulent mixed drift-vortex structure in some other regions. In this sense, plasma motion is turbulent in some regions while a simple vortex motion along the magnetic field lines can dominate in some others.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. The parameter sets studied in this paper.

Run-name	${{\kappa}_{\parallel}}$	${{\epsilon}_{H}}$	SGS model	${{C}_{\nu}}$	${{C}_{\eta}}$	${{C}_{\kappa}}$	Low-pass filtera
LES41	10−1	0.0125	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.046	0.046	0.046	off
LES42	10−1	0.0125	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	0.46	0.46	off
LES45	10−1	0.0125	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	0.46	0.046	off
LES49	1	0.0125	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.046	0.046	0.046	off
LES50	1	0.0125	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	0.46	0.46	off
LES118	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	0.92	0.46	off
LES119	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	1.38	0.46	off
LES120	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.23	1.38	0.46	off
LES121	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.092	1.38	0.46	off
LES122	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	0.92	0.46	on
LES123	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	1.38	0.46	on
LES124	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.23	1.38	0.46	on
LES125	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.092	1.38	0.46	on
LES126	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	0.92	0.92	off
LES127	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	1.38	0.92	off
LES128	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.23	1.38	0.92	off
LES129	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.092	1.38	0.92	off
LES130	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	0.92	0.92	on
LES131	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.46	1.38	0.92	on
LES132	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.23	1.38	0.92	on
LES133	10−1	0.05	${{\mu}_{\text{SGS}}}=\mu _{\text{SGS}}^{a}$	0.092	1.38	0.92	on

^aBesides the low-pass filter for the LES, another low-pass filter of the fourth order explicit formulation presented in [28] is operated on the fluctuation components of the density, momentum and pressure fields at the end of the fourth stage of the Runge–Kutta marching step in order to remove high-wavenumber components when it is 'on'.

In figure 5, the time evolution of the Fourier energy of each toroidal mode n of (a) parallel and (b) normal components of the velocity obtained by the LES is shown. In figure 5 the Fourier energies grow exponentially at the beginning of the time evolution. The growth rate is largest in n  =  1 due to a large ${{\kappa}_{\parallel}}=1\times {{10}^{-1}}$. Furthermore, the parallel component is much larger than the normal component, being similar to our single-fluid MHD results in [6, 7]. These results show that our LES reproduces a qualitative nature in the ballooning modes of LHD with a higher numerical resolution in [7]. We have also compared the results with a simulation that does not have the SGS terms. Although a simulation which is carried out without the SGS terms and without a sufficient resolution cannot go beyond nonlinear saturation, it can reproduce the beginning of the linear regime of low modes. So far as we have compared the growth rates between the two runs, LES41, and the run without the SGS terms, the growth rates of the low modes differ by only a few percent or less (the figure is omitted). This is quite natural, because the SGS models are essentially nonlinear, so that the SGS terms begin to work in the nonlinear phase, and do not contaminate the linear regime. Note here that the plasma dynamics can be contaminated by the SGS terms when the SGS coefficients are too large, and we must restrict ourselves so that the SGS coefficients are not too large. Since the SGS coefficient is given by equation (27), making ${{C}_{\nu}}$ and ${{C}_{\eta}}$ ten times larger simultaneously can be translated as $\Delta$ to ${{10}^{1/4}}\simeq 1.78$. In other words, the filter length $\Delta$ is effectively doubled by enlarging ${{C}_{\nu}}$ and ${{C}_{\eta}}$ ten times, and the unstable mode resolved in the grid scale is halved. As a consequence, it should appear in an LES as the energy loss from the GS to SGS (lower GS energy in comparison to an LES with smaller Smagorinsky coefficients.) We intend to study the low-mode—$n\leqslant 3$ in the coarsest case and $n\leqslant 6$ for a more preferable case—without the contamination of the SGS terms. In [16], ${{C}_{\nu}}={{C}_{\eta}}=0.046$ is used based on their turbulence theory. Here, we should not forget that the applicability of the LES approach is based on our earlier research in which the forward energy transfer becomes dominant at scales smaller than the ion skin depth. This also restricts the values of ${{C}_{\nu}}$ and ${{C}_{\eta}}$, or $\Delta$ equivalently. This is the origin of the SGS coefficients for LES41, although ${{C}_{\nu}}={{C}_{\eta}}=0.345$ can work better than ${{C}_{\nu}}={{C}_{\eta}}=0.046$ for reproduction of the spatial structures in turbulence with the Hall effects [5]. Supposing the resolution of $n\leqslant 3$ in the coarsest case, we consider ${{C}_{\nu}}={{C}_{\eta}}\leqslant 1$ to be preferable, and ${{C}_{\kappa}}=1.38$ adopted in some runs can be quite marginal, although there is no clear estimation of ${{C}_{\kappa}}$, because this is the first case to introduce the SGS model for our extended MHD system. This is the reason why we have restricted the Smagorinsky coefficients in the range seen in table 1. Since we carry out LESs which aim to solve low modes with the help of SGS effects, we should verify how low modes are solved in the linear stage of a simulation. The pressure and velocity Fourier coefficients in a computation in [7] (the number of grid points $193\times 193\times 640$ and the parameters $\mu =\eta ={{\kappa}_{\parallel}}={{\kappa}_{\bot}}=1\times {{10}^{-6}}$) agree well with those of the linear eigenfunctions at each m and n for $n\leqslant 12$. We make use of the results there as reference data. Since the numerical resolution in the radial and poloidal directions in the LESs here is half that of the reference data, we can expect the resolution of the Fourier coefficients of n  <  6 for the linear phase of our LESs.

In figure 6, the radial profiles of the Fourier coefficients of the pressure for (a) n  =  1 and (b) n  =  3 with some poloidal wavenumber m in the run LES41 are plotted as a function of the initial magnetic flux function $\sqrt{\psi}$. A combination of the wavenumbers m and n is provided for $\iota =\mid m/n\mid \simeq 0.5$ because the pressure gradient is steepest there. The radial profile is plotted at a time when the time evolution shifts from the linear stage (exponential growth of the Fourier energy) to the nonlinear stage. It is seen in figure 6(a) that the Fourier coefficients at $\mid m/n\mid \simeq 0.5$ have the largest amplitude at that time, and the Fourier coefficients peak at $\sqrt{\psi}\simeq 0.4$ where the pressure gradient is steepest. (The vertical axis in the plots is essentially arbitrary in the sense that the Fourier coefficients form a linear eigenfunction of the unstable modes for each n.) In figure 6(b), the Fourier coefficient of m/n  =  −5/3—not of m/n  =  −6/3—has the largest amplitude and its peak appears at $\sqrt{\psi}\simeq 0.43$. The dominance of m/n  =  −5/3 for the n  =  3 mode, in spite of the large pressure gradient at $\iota /2\pi \simeq 0.5$, is a consequence of the coupling with the m/n  =  1/0 coefficient, which is a residue from obtaining the initial equilibrium by using the HINT code [19]. (The influence of m/n  =  1/0 has been already reported in [6].) Since it is the result coming from the initial equilibrium, we consider the linear nature to be well resolved for n  =  3 too. Fourier coefficients in LES41 for $n\geqslant 5$ do not necessarily form a clear radial profile, partially because the high parallel heat conductivity ${{\kappa}_{\parallel}}=1\times {{10}^{-1}}$ tends to suppress the formation of the Fourier coefficients to be like linear eigenfunctions as n becomes larger (a computation with ${{\kappa}_{\parallel}}={{\kappa}_{\bot}}=1\times {{10}^{-6}}$ has been used for studying the radial profile in [7]), and partially because of the influence of the SGS terms.

Zoom In Zoom Out Reset image size

In figure 6(c), the radial profile of the n  =  1 Fourier coefficients in LES42 is plotted at the same time as (a) and (b). The outline of the radial profile is similar to that in figure 6(a), specifically for the main harmonic m/n  =  −2/1. However, we may find that the main harmonic m/n  =  −2/1 is slightly different between figures 6(a) and (c), at $\sqrt{\psi}\simeq 0.6$. The harmonic in figure 6(c) changes more slowly than that in (a). Although figure 6(c) is for the plot of the linear stage in LES42, it might come from the SGS coefficients in LES42, especially from those for the heat conduction ${{C}_{\kappa}}$, which are larger than those in LES41 by ten times or more. In other words, adopting SGS coefficients that are too large can contaminate the linear stage of unstable modes, even though the SGS models in equations (28)–(30) are provided so as not to contaminate the linear stage. Based on the parameter surveys over a wide range of the SGS coefficients ${{C}_{\nu}}$, ${{C}_{\kappa}}$ and ${{C}_{\eta}}$, we restrict our computational works here to those shown in table 1 so that the nonlinearity in the SGS terms does not contaminate the linear evolution too much.

In figure 7, the radial profile of the m/n  =  0/0 component of the pressure is plotted at some moments from the initial to the saturation time. Figure 7(a) is for an LES of (a) ${{C}_{\nu}}={{C}_{\kappa}}={{C}_{\eta}}=0.046$, and figure 7(b) is for an LES of (b) ${{C}_{\nu}}$  =  ${{C}_{\kappa}}$  =  ${{C}_{\eta}}=0.46$. Both of the two runs give nonlinear saturation even though the numerical resolution is not very fine. (The saturated state is not perfectly in a steady state because of the absence of the source and sinks. However, the profiles are in a slowly decaying state after the saturation of the energy growth and a drastic change in the profiles is not observed in the phase.) Though the SGS viscosity becomes high in the turbulent regions of the LESs, as we have estimated in the previous section, the plasma core in figure 7(a) finally collapses. The collapse is a consequence of the growth of the low modes without artificial damping by a large viscosity. The growth is sharply different from that of earlier simulations with $\mu =5\times {{10}^{-4}}$ in [5], in which the growth rate of the linear low modes is considerably lowered by the viscosity. This shows that a high and uniform viscosity should not be adopted for the study of instability, and that the SGS viscosity which represents the energy transfer from the GS to the SGS components should be used instead. While figure 7(a) gives a deteriorated profile of P₀₀, the set of SGS coefficients in (b) gives a better result than (a).

Zoom In Zoom Out Reset image size

Recall here that flattening of the total pressure is not observed in experiments [17]. In other words, we consider the baseline as keeping a peaked pressure at the core ($\sqrt{\psi}\simeq 0$). Since figure 7(b) allows total collapse of the pressure to be avoided, meaning a low but finite peak is kept at the center of the core, the LES approach is considered to be able to reasonably compensate for a lack of numerical resolution in a simulation and enable a smaller, quicker numerical simulation. However, we must keep in mind that the SGS coefficients should be calibrated more carefully, so that the SGS models represent the influence of turbulent fluctuation on linear low modes, in accordance with simulations with a higher resolution as well as experiments. Also note that the computational cost is about 1/32 or less than the earlier computations reported in [10], in which the number of grid points is $193\times 193\times 640$ and a strong numerical filter is applied to eliminate very fine (both in time and space) motion.

Now we carry out the LESs of the ballooning modes to see how the SGS coefficients can influence the saturated profile of the pressure, for SGS coefficients that are larger than those in LES41 and LES42. Although appropriate values of the SGS coefficients should be determined through calibration with large computations of very high resolution, without SGS models as well as with many experimental results, the parameter survey over a wide range of SGS coefficients can clarify the possible contributions of the SGS terms to the saturation profile of the pressure and dynamics of the low modes.

In figure 8(a), saturated P₀₀ in runs LES118, LES119 and LES120 are plotted as a function of $\sqrt{\psi}$. While P₀₀ in LES118 and LES120 are hollow at the plasma core, P₀₀ in LES119 peaks at $\sqrt{\psi}\simeq 0$. The pressure deformation associated with advection is considered to be relatively suppressed in LES119 because ${{C}_{\eta}}$ in LES119 is larger than that in LES118 and LES120. This observation means that turbulent fluctuation on the SGS scale can interact with unstable modes and suppress the pressure deformation through the nonlinear interaction. In figure 8(b), the change of P₀₀ in time in the run LES119 is plotted. As we have already seen for other parameter sets, the P₀₀ profile deteriorates as the unstable modes grow. However, at the final stage of the saturation, the plasma core avoids becoming hollow. This suggests that turbulent fluctuations, the roles of which are represented through SGS terms, can work to terminate the final crash by strong nonlinear couplings, although it is difficult for them to suppress the growth of unstable low modes completely, because unstable low modes can grow to sufficiently large amplitudes before turbulent fluctuation terminates the growth of the low modes.

Zoom In Zoom Out Reset image size

Finally, it is mentioned that some parameters have been examined but not presented in this paper. For example, we have examined the influence of an additional low-pass filter on the strong stabilization of the computations. A distinction of the adoption of the additional low-pass filter is indicated by 'on'/'off' in the table. It is noted here that the introduction of the low-pass filter does not affect numerical results significantly.