Elsevier

Journal of Computational Physics

Volume 295, 15 August 2015, Pages 24-45
Journal of Computational Physics

Stable explicit coupling of the Yee scheme with a linear current model in fluctuating magnetized plasmas

https://doi.org/10.1016/j.jcp.2015.03.069Get rights and content

Abstract

This work analyzes the stability of the Yee scheme for non-stationary Maxwell's equations coupled with a linear current model with density fluctuations. We show that the usual procedure may yield unstable scheme for physical situations that correspond to strongly magnetized plasmas in X-mode (TE) polarization. We propose to use first order clustered discretization of the vectorial product that gives back a stable coupling. We validate the schemes on some test cases representative of direct numerical simulations of X-mode in a magnetic fusion plasma including turbulence.

Introduction

The study of the impact of high turbulence level on the wave propagation introduces unexpected behavior on the stability of the numerical scheme, indeed numerical instabilities may arise after a very huge number of time-steps, as it will be shown. At the same time due to the high turbulence level, it becomes irrelevant to use the standard analysis using dispersion relation, that is the reason we propose to use another analysis based on energy to establish the stability criterion of the numerical scheme. A detailed study is provided to explain the way of thinking, and to facilitate the understanding of this work. Up to now, the role of turbulence on the wave propagation was ignored, but it becomes more and more important to explain experiments. The numerical modeling of these observed macroscopic phenomenons resulting from turbulent processes in magnetic plasmas requires huge numerical simulations. Such simulations should capture the statistical properties so that one can realize relevant averaging in time and space. This is particularly true to explain the macroscopic wave behavior which takes place on time scales that diagnostics cannot measure. An example which motivates our work is the recovery of the smooth envelope of a probing beam: to be able to compare with the theory, we need to simulate the non-stationary Maxwell equations with a linear current model, see Eqs. (1), (2) below, in a domain with very fine non-homogeneities (since the underlying plasma is turbulent) and over a simulation time equivalent to a half-million wave periods. This is required to study in detail the role of the wavenumber spectrum on the beam widening [7], [8], [26], and is illustrated with a snapshot of a numerical calculation of a wave in a plasma in Fig. 1. The Gaussian shape is recovered at coarse scale using some time averaging, but at the same time very fine spatial scale details are visible. Since the energy of the wave is conserved during its propagation (we will prove this property for the model used in this work), it seems unavoidable to ask for similar energy preservation properties for the numerical methods in order to capture both the fine and coarse scales in the computations.

The previous example corresponding to beam broadening induced by wave propagation in turbulent plasma is illustrative of diagnostics and wave heating systems used in magnetized fusion plasmas. Hot topics for the ITER design concern the wavenumber resolution of the so-called Doppler reflectometry or on the choice of the probing frequency of the Coherent Thomson Scattering using X-mode as well as the beam widening induced by the turbulence on the electron cyclotron heating to determine whether neo-classical tearing modes inducing big islands can be controlled or not [21]. The same questions arise for the lower-hybrid heating system to evaluate correctly what is the wavenumber spectrum launched into the plasma in view of the efficiency prediction of the non-inductive current driven system and the definition of the energy deposition zone. In the near future, wave polarization changes will be a subject of great interest, as it is becoming with O–X mode heating scenario for the stellarator or tokamak [18], which are subject to anomalous reflection due to turbulence [22], and also for the new diagnostic concepts based on wave polarization changes induced by linear mode conversion or turbulence in inhomogeneous magnetized plasmas of space, astrophysics or fusion.

To our knowledge the state of the art of the computations done for such problems use simplified wave models taking into account only one physical mechanism (refraction in [21], diffraction effects in [31] and recently both in [16], [17]) but not all the wave structure as a Maxwell's equation solver with a linear current model can predict. This strongly questions the numerical stability of the simulations, indeed our numerical tests show a fundamental stability issue on extraordinary mode (also called TE or X-mode) which are more demanding in terms of stability than O-mode computations. Because an ideal scheme should be fast enough to reach the requirements stressed above for the comparison between theory and full-wave simulations, we are strongly interested in explicit FDTD (Finite Difference Time Domain) schemes, indeed they are still the cheapest in terms of CPU. Many FDTD methods were developed to improve the performances and the stability in the past few years, such as EJ collocated FDTD [34], Runge–Kutta exponential time differencing formulation (RKETD) FDTD [20], matrix exponential (ME) FDTD [14], and exponential time differencing (ETD) FDTD [15]. Unfortunately they cannot be easily used to model the general dispersive media. On the other hand, in order to overcome the Courant–Friedrich–Levy (CFL) constraint of the conventional FDTD method, the one-step leapfrog ADI-FDTD method has been developed [4]. It originates from the conventional ADI-FDTD for which we refer to [35]. We finally mention the fact that conservation of energy for a fourth-order FDTD scheme is addressed in [19], and that massive direct simulations of wave propagation in a different physical context can be found in the recent references [12], [10], [23]. A recent variant of FDTD methods with higher computational efficiency is the implicit method [27] or [30], [29], the latter being proved to be stable for varying coefficients. However ITER size simulations for the high frequency cases with small spatial scale turbulence that we consider in Fig. 1 seem not reachable using implicit schemes such as [30]. The reason is the very large number of cells (typically 103d with d=2,3 the dimension of the problem) which induces severe CPU constraints. This is why we study in this work the low cost explicit version of the Yee scheme coupled with a linear current. The stability of the method will be performed with an energy analysis which seems to us very convenient to handle strong gradients of the coefficients.

The general model problem considered in this work is the non-stationary Maxwell system{ε0tE+curlH=Jμ0tH+curlE=0 coupled with a linear equation for the electronic current density J=eNeue,metue=e(E+ueB0). The unknowns are the electromagnetic field (E,H), with the usual notation that H=Bμ0, and the electronic current is J. Eq. (2) is the linearized Newton law for electrons where the Lorentz force on the right hand side is obtained by linearization [2], [3], [13], [24] of the full Newton law in a strongly magnetized medium with bulk magnetic field B00, and vanishing electric field and current E0=J0=0. The coefficients are the mass me and charge e<0 of the electron, the permittivity ε0 and permeability μ0 of vacuum. The electronic density Ne=Ne(t,x) is a prescribed function of the space and time variables. The background magnetic field could as well be dependent of the space and time variables B0=B0(t,x). The plasma frequency and cyclotron frequency are defined byωp=e2Nemeε0 and ωc=|eB0|me. At any time, the total energy of the system is naturally the sum of the electromagnetic energy and the kinetic energy of electronsE=Ω(ε0|E|22+|B|22μ0+me|J|22|e|2Ne)dx. This quantity will play the key role in our analysis.

The specific numerical difficulty addressed here is the stability of such couplings in the regime where the electronic density Ne has extremely strong gradients in space induced by the high level of the turbulence|Ne|>>1. See for example the density profiles used in the numerical section where the density fluctuations are able to remove plasma and the wave behavior is identical to those in vacuum that is the reason why the cut-off frequency goes to zero, no cut-off layer in the vacuum zones. On the other hand |B0| can be considered small, even zero in first approximation. One may also consider that tNe0 since this very simplified situation is representative of some reflectometry problems of fusion plasmas in X-mode polarization. Numerical instabilities that can occur for such problems are usually related with the existence of different cut-offs and resonances [24], [25]: the cyclotron resonance is such that ω=ωc(x) and the hybrid-resonance is such that ω2=ωp(x)2+ωc(x)2. Since for a given ω there may always exist a solution x to these equations, it can yield a singular solution in time harmonic domain [11] and so be the germ of numerical instabilities for time domain simulations. This phenomenon will be evidenced in our numerical section for the hybrid resonance. The cut-off ω=ωp(x) is not a singularity but a mild transition for propagative zone to elliptic (non-propagative). For reasons explained a the beginning of this work, we desire to develop an unconditional stable explicit algorithm that can handle all these situations. This is why we do not rely on a dispersion analysis which is valid mainly for locally frozen coefficients problems, but more on the energy criterion.

To be more specific let us assume now that b=B0|B0|=(0,0,1)t is constant. It is not difficult to see that the initial problem (1), (2) can be decomposed in two different polarizations [24], [25]. The first polarization is transverse magnetic so the electric field together with the electric current is always parallel to b. In this case the vector product bJ vanishes and the equation can be simplified a lot. This situation that does not yield resonances but only cut-off is called O-mode in plasma physics. The other polarization is transverse electric so that bJ is a priori non-zero. This polarization is called X-mode in the language of plasma physics. It is fundamentally the one we are interested in since the cyclotron resonance ω=ωc(x) and the hybrid-resonance ω2=ωp(x)2+ωc(x)2 can show up as the result of the resonant interaction of bJ0 with the rest of the system. The equations write{ε0tEx+yHz=Jx,Jx=eNeux,ε0tEyxHz=Jy,Jy=eNeuy,μ0tHz+xEyyEx=0,metux=eEx+euyBz0,metuy=eEyeuxBz0. This system, as the previous one (1), (2), is linear and will be used for the numerical tests, even if the analysis is done for the full 3D problem.

The main contributions of this work are in two directions: first we develop a theoretical framework that provides stable couplings in the energy norm between the Yee scheme (in any dimension) and a linear current model. Second, we design two simple explicit coupling schemes that are proved to be stable in the energy norm, a property that is generally not properly satisfied by existing FDTD schemes for such problems. The stability is achieved for all electronic densities, including cases with cut-offs or resonances in the computational domain. We will also prove that the CFL condition of one of our method is exactly the same as the one of the Yee scheme for Maxwell's equations in vacuum. Numerical tests performed for (6) with physically based coefficients will illustrate the properties of the new scheme.

The plan is as follows. In Section 2 we recall the basics of the Yee scheme coupled with a linear current model. We shall outline the formal procedure that allows to derive explicit schemes and apply a standard stability analysis to different leap-frog time discretizations. Because in Yee schemes the different components of the vector fields are defined on staggered grids, care must be taken when defining the discrete curl and vector product operators. For this reason we introduce abstract notations for these operators and we identify the conditions under which both the stability analysis and the formal derivation of an explicit scheme can be reproduced at the discrete level. In Section 3 we describe specific operators that satisfy these conditions, a property that is not fulfilled in usual FDTD schemes, see Remark 3.5. Equipped with these operators we then propose two explicit schemes that are proved to be stable under standard CFL conditions, and we give detailed formulas specifying the abstract operators. In Section 4 we finally show numerical experiments that demonstrate the improved stability of this approach.

Section snippets

Yee scheme with linear current model

The standard method to compute a numerical approximation of (1), (2) on a Cartesian grid consists of coupling the Yee scheme with a linear discretization of the current, see [32], [1]. To have better insight in the structure of the coupled problem (1), (2) we eliminate the velocity in terms of the current and discuss a discretization procedure. Here the data Ne is constant in time. Extension to the time-dependent case is postponed to Appendix C and Section 4.2.4.

Thus, the model problem is

Staggered grid operators and explicit schemes

This section is the core of this work. In the previous section we have studied under which conditions the abstract semi-implicit schemes (15) and (21) were stable. We now turn to the problem of deriving fully explicit formulations that are stable. Specifically, we will design a vector product operator acting on edgewise vector fieldsbh:VedgeVedge that allows to derive stable and explicit schemes when used in conjunction with a proper operatorS(u):VedgeVedge for the multiplication of an

Numerical results

Our numerical tests are relevant for reflectometry studies in fusion plasmas or ionospheric probing. The domain of computation is 2D and we consider X-mode configuration, meaning that polarization is transverse electric as in (6). The 2D implementation maintains all the difficulties studied in this work and is simpler to implement and validate. We present three test problems on time-independent plasmas which are the main subject of this work, together with one additional time-dependent plasma

Conclusion

We have detailed the influence of strong gradients in the plasma parameters on the stability in the energy norm of a Yee scheme coupled with a linear current model. In such conditions the WKB approximation is no longer valid and the use of the dispersion relation for stability analysis seems to be inappropriate in the studied cases. The main asset is the identification of an instability phenomenon triggered by strong gradients of Ne. This is related to the use of a centered evaluation of the

Acknowledgements

The authors acknowledge the support of ANR under contract ANR-12-BS01-0006-01. Moreover, this work was carried out within the framework of the European Fusion Development Agreement and the French Research Federation for Fusion Studies. It is supported by the European Communities under the contract of Association between Euratom and CEA, and has received funding from the European Union's Horizon 2020 Research and Innovation Programme under grant agreement number 633053. IST activities also

References (36)

  • F. da Silva et al.

    Studies on O-mode reflectometry spectra simulations with velocity shear layer

    Nucl. Fusion

    (2006)
  • F. da Silva et al.

    A numerical study of forward- and back-scattering signatures on Doppler reflectometry signals

    IEEE Trans. Plasma Sci.

    (2010)
  • F. da Silva et al.

    Developments on reflectometry simulations for fusion plasmas: applications to ITER position reflectometry

    J. Plasma Phys.

    (2006)
  • F. da Silva, S. Heuraux, T. Ribeiro, B. Scott, Development of a 2D full-wave JE-FDTD Maxwell X-mode code for...
  • S. Del Pino et al.

    3D finite volume simulation of acoustic waves in the earth atmosphere

    Comput. Fluids

    (2009)
  • D.Y. Heh et al.

    FDTD modeling for dispersive media using matrix exponential method

    IEEE Microw. Wirel. Compon. Lett.

    (2009)
  • S. Huang et al.

    FDTD implementation for magnetoplasma medium using exponential time differencing

    IEEE Microw. Wirel. Compon. Lett.

    (2005)
  • Joo Hwa Lee et al.

    Three-dimensional FDTD simulation of electromagnetic wave transformation in a dynamic inhomogeneous magnetized plasma

    IEEE Antennas Wirel. Propag. Lett.

    (1999)
  • Cited by (20)

    • Dispersion relations in hot magnetized plasmas

      2018, Journal of Mathematical Analysis and Applications
      Citation Excerpt :

      They furnish the phase velocity, the group velocity and the refractive index. They allow to determine the eikonal equation, and they are crucial in reflectometry [11,21]. They are a prerequisite to understand turbulence phenomena [7].

    View all citing articles on Scopus
    View full text