Stable explicit coupling of the Yee scheme with a linear current model in fluctuating magnetized plasmas
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 with 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 coupled with a linear equation for the electronic current density , The unknowns are the electromagnetic field , with the usual notation that , 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 , and vanishing electric field and current . The coefficients are the mass and charge of the electron, the permittivity and permeability of vacuum. The electronic density 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 . The plasma frequency and cyclotron frequency are defined by At any time, the total energy of the system is naturally the sum of the electromagnetic energy and the kinetic energy of electrons 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 has extremely strong gradients in space induced by the high level of the turbulence 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 can be considered small, even zero in first approximation. One may also consider that 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 and the hybrid-resonance is such that . 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 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 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 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 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 and the hybrid-resonance can show up as the result of the resonant interaction of with the rest of the system. The equations write 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 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 fields that allows to derive stable and explicit schemes when used in conjunction with a proper operator 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 . 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)
- et al.
Hybrid resonance of Maxwell's equations in slab geometry
J. Math. Pures Appl.
(May 2014) - et al.
Simulation of laser propagation in a plasma with a frequency wave equation
J. Comput. Phys.
(1 February 2008) - et al.
Effects of non-Maxwellian species on ion cyclotron waves propagation and absorption in magnetically confined plasmas
Phys. Plasmas
(2005) - et al.
Runge–Kutta exponential time differencing FDTD method for anisotropic magnetized plasma
IEEE Antennas Wirel. Propag. Lett.
(2008) - et al.
Electron-cyclotron wave scattering by edge density fluctuations in ITER
Phys. Plasmas
(2009) Simulation of microwave propagation in a fusion plasma
(January 2011)Kinetic Theory of Plasma Waves – Homogeneous Plasmas
(1998)- et al.
Amplification and absorption of electromagnetic waves in overdense plasmas
Plasma Phys.
(1974) - et al.
A leapfrog formulation of the 3D ADI-FDTD algorithm
Int. J. Numer. Model.
(2009) An analysis of new and existing FDTD methods for isotropic cold plasma and a method for improving their accuracy
IEEE Trans. Antennas Propag.
(1997)
Studies on O-mode reflectometry spectra simulations with velocity shear layer
Nucl. Fusion
A numerical study of forward- and back-scattering signatures on Doppler reflectometry signals
IEEE Trans. Plasma Sci.
Developments on reflectometry simulations for fusion plasmas: applications to ITER position reflectometry
J. Plasma Phys.
3D finite volume simulation of acoustic waves in the earth atmosphere
Comput. Fluids
FDTD modeling for dispersive media using matrix exponential method
IEEE Microw. Wirel. Compon. Lett.
FDTD implementation for magnetoplasma medium using exponential time differencing
IEEE Microw. Wirel. Compon. Lett.
Three-dimensional FDTD simulation of electromagnetic wave transformation in a dynamic inhomogeneous magnetized plasma
IEEE Antennas Wirel. Propag. Lett.
Cited by (20)
An overview of the evolution of the modeling of reflectometry diagnostics in fusion plasmas using finite-difference time-domain codes
2024, Fusion Engineering and DesignSimulation and data processing techniques to design optimized PPR systems on plasma fusion devices
2024, Computer Physics Communications3D full-wave computation of RF modes in magnetised plasmas
2019, Computer Physics CommunicationsDispersion relations in hot magnetized plasmas
2018, Journal of Mathematical Analysis and ApplicationsCitation 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].
Stable perfectly matched layers for a cold plasma in a strong background magnetic field
2017, Journal of Computational PhysicsSimulation as a tool to improve wave heating in fusion plasmas
2015, Journal of Plasma Physics