A Multiscale Approach Using Patches of Finite Elements for Solving Wave Propagation Problems in Microwave Discharge Plasma

1. Introduction

There is an increasing interest of microwave plasma sources at atmospheric pressure for medicine, waste treatment and more specifically for aeronautic and space applications. Energy deposition by microwave discharges has received significant interest in the past several decades as a promising technique for aerodynamic flow control at high speed or plasma assisted combustion [1]. For instance, recent experiences [2] have shown a dramatic reduction of the surface pressure coefficients due to the strong interaction of a microwave discharge with a Mach 3 incoming flow. In the domain of space propulsion, a concept of rocket based on beamed energy propulsion [3] has been proposed for an additional space launch system in the future. Solid rocket boosters are replaced by microwave rocket in which thrust is generated through repetitively pulsed microwave detonation. Properties of microwave plasma are also used in the plasma arrays with 2-D periodicity producing frequency regions of forbidden propagation like band gaps in a photonic crystal [4].

Modeling and numerical simulation are a fundamental step to understand the various mechanisms of the formation of the microwave discharge and how it interacts with the ambient gas. Experiences carried out at MIT [5–8] have shown a microwave breakdown in air at atmospheric pressure using a 110 GHz pulsed gyratron in the microwave power range. Numerical simulations based on a plasma fluid model coupled with Maxwell's equations have been able to reproduce many of the experimental features in 2D [9–11] inlcuding the gas heating effects [12] and in 3D [13]. The spatial and time scales involved in the mutual interaction of the microwave fields and plasma are particularly disparate making the full simulation of such a non-linear dynamics computationally expensive and often impracticable. Indeed, to properly capture the typical plasma diffusion length, the spatial step size needs to be very small with respect to the wavelength (typically 5 or 10 times smaller) leading to time steps on the order of 10⁻¹⁴ s. As such, an effort has to be done to increase the performance of the simulation tools. In microwave plasma modeling, the Maxwell's equations are generally solved by the Finite-Difference-Time-Domain (FDTD) method using explicit time-integration techniques requiring strong time step restrictions for stability. In the past, some efforts haven been accomplished to reduce the computational cost of the FDTD method by implicitly integrating in time the dynamic fields [14]. However, the numerical diffusion of the implicit technique produces a loss of accuracy in the results, which may become significant in simulations over long time scales. A multiscale strategy has to be developed to deal with the different characteristic lengths. The FDTD methods are very efficient from the computational point of view, but present strong geometrical constraints and then, are not directly applicable to sub-gridding. Several FDTD-adaptation techniques have been proposed in the past. On the one hand, we find the space time subgrids where each grid runs at a time step defined by its own CFL stability limit while interpolations are connecting the grids together [15, 16]. Differently, subgrids can run at the time step of a coarse grid and the fine grids are stabilized by means of either an unconditionally stable method such as the Alternating-Direction-Implicit (ADI) FDTD method [17], or a priori removal of unstable grid eigenmodes or filtering of unstable spatial harmonics [18].

This paper presents an adaptation of the “method of finite element patches” to electromagnetic wave propagation problems and an application to simulation of microwave discharge. The patch-based method is an efficient approach to deal with the multiscale behavior of a problem and the coupling of different models without any constraint of conformity between the meshes. It is a Schwarz domain decomposition algorithm and has been first proposed by Glowinski et al. [19] for linear elliptic problems and then extended to derive diffusion [20]. It has similarities with the Chimera method [21], used by the Computational Fluid Dynamics community, or with the Arlequin method [22] for structural design. An extension to non-linear elliptic problems has been also successfully proposed by Brunet et al. [23] and applied to simulation of a plasma around a negatively charged array of solar generator interconnects. In this study, we derive the multiscale resolution algorithm for the variational formulation of the standard Yee scheme [24] and we evaluate its efficiency on standard grid configurations such as nested Cartesian grids. The proposed patch-based approach provides the finite-difference type scheme with the local-grid refinement capacity while preserving its stability. Note that such a method is not limited to the studied scheme, rather it can be inherently extended to any other variational formulation, such as Nedelec's mixed finite element formulations [25], Discontinuous Galerkin formulations [26, 27], etc. This indeed highlights one of the main strengths of the proposed approach in microwave discharge plasma modeling.

The paper is organized as follows. In section 2, the physical model of microwave discharge plasma is presented. In section 3, the patch-based method is derived for the Yee scheme and a theoretical estimate of its computational complexity is given. In section 4, numerical experiments are presented to illustrate the accuracy and the computational efficiency of the proposed method in one dimension, for academic cases, and in two dimensions, simulating a microwave breakdown in air under atmospheric pressure.

2. Physical Model

The physical model describing the dynamics of microwave streamers consists of plasma transport equations coupled with Maxwell's equations via the electron current density. Particularly, microwave oscillations are described by Maxwell's equations:

$\begin{array}{l}{ϵ}_0\frac{\partial \tilde{\text{E}}}{\partial t}=\nabla \times \tilde{\text{H}}-\tilde{\text{J}}\\ {\mu }_0\frac{\partial \tilde{\text{H}}}{\partial t}=-\nabla \times \tilde{\text{E}}& (1)\end{array}$

where $\tilde{\text{E}}(\text{x},t)$ and $\tilde{\text{H}}(\text{x},t)$ are the electric and magnetic fields, ϵ₀ and μ₀ are the free space permittivity and permeability. $\tilde{\text{J}}(\text{x},t)$ is the electron current density through which the coupling with plasma dynamics is established:

$\begin{array}{ll}\tilde{\text{J}}=-e{n}_e{\tilde{\text{v}}}_e.& (2)\end{array}$

In Equation (2), e, n_e(x, t), and ${\tilde{\text{v}}}_e(\text{x},t)$ are the electron charge, density, and mean velocity, respectively. Due to the relatively large ion mass the ion motion is neglected with respect to the electron motion. The instantaneous electron velocity in Equation (2) is determined by the approximate electron momentum transfer equation:

$\begin{array}{ll}\frac{\partial {\tilde{\text{v}}}_e}{\partial t}=-\frac{e\tilde{\text{E}}}{m_e}-{\nu }_m{\tilde{\text{v}}}_e,& (3)\end{array}$

where m_e is the electron mass and ν_m the electron-neutral momentum transfer frequency given by [28]

$\begin{array}{ll}{\nu }_m=\left(5.3\times 10^9\frac{\text{Hz}}{\text{Torr}}\right)p,& (4)\end{array}$

where ν_m is in s⁻¹ and p is the gas pressure at ambient temperature in Torr. It is important to note that Equation (3) does not take the force exerted by the magnetic field into account. As a matter of fact, in our conditions the magnetic term $({\tilde{\text{v}}}_e\times \tilde{\text{B}})$, with $\tilde{\text{B}}={\mu }_0\tilde{\text{H}}$, is negligible with respect to the electric one $\tilde{\text{E}}$ since the mean electron velocity is much smaller than the speed of light c, $\left|{\tilde{\text{v}}}_e\right|\ll c$, and the ratio $\left|\tilde{\text{E}}\right|/\left|\tilde{\text{B}}\right|$ is of the same order of c. The plasma is supposed to be quasi-neutral and its evolution in time and space is obtained from a system of reaction-diffusion equations on the electron density and mean energy averaged over a microwave field period (note that only quantities varying over the wave cycle are denoted with tilde)

$\begin{array}{l}\frac{\partial n_e}{\partial t}-\nabla ·(D_{eff}\nabla n_e)={\nu }_{eff}{n}_e-r{n}_e^2\\ \frac{\partial n_{\epsilon }}{\partial t}-\frac{5}{3}\nabla ·(D_{eff}\nabla n_{\epsilon })=-n_e(e<{\tilde{\text{v}}}_e·\tilde{\text{E}}{>}_{T_M}+\text{Θ}).& (5)\end{array}$

The right-hand side of the first equation represents the charged particle production and losses. ν_eff is an effective electron impact ionization frequency including ionization and attachment and r is the electron-ion recombination coefficient. ν_eff depends on the local electron distribution function and its macroscopic description is obtained in this context by means of the local mean energy approximation. Therefore, the effective ionization coefficient is expressed as a non-linear function of the local mean electron energy ε_e(x, t) which is obtained by the second equation of system (5) through the intermediate variable, the energy density n_ε = n_eε_e. Note that the kinetic energy is neglected with respect to the thermal energy so that ε_e = 3k_BT_e/2 (since the electron distribution is not Maxwellian this equation defines the electron temperature T_e, with k_B being the Boltzmann constant). The right-hand side of this equation represents the energy gain and losses during inelastic and elastic collisions (we consider in this paper only collisions with molecules in the ground state) per unit time by electrons under a varying electric field. Particularly, the term $-e<{\tilde{\text{v}}}_e·\tilde{\text{E}}{>}_{T_M}$ expresses the mean energy gained per unit time by an electron (it is averaged over the entire microwave period T_M = 2π/ω, with ω being the wave angular frequency) and its value depends on the phase shift between the instantaneous electron velocity and the microwave field

$\begin{array}{ll}-e<{\tilde{\text{v}}}_e·\tilde{\text{E}}{>}_{T_M}=-\frac{e}{T_M}{\displaystyle {\int }_t^{t+T_M}}{\tilde{\text{v}}}_e·\tilde{\text{E}}dt.& (6)\end{array}$

The mean power loss per electron is identified by

$\begin{array}{ll}\text{Θ}={\nu }_{\epsilon }{\epsilon }_e& (7)\end{array}$

where ν_ε is the energy relaxation frequency and it takes into account the rate coefficients related to elastic and inelastic collisions. Note that in the local mean energy approximation all rate and transport coefficients depend on ε_e as well as ω/N (with N being the gas density). As a result, Θ is a non-linear function of ε_e. The other transport coefficient of Equations (5), D_eff depends on the electron temperature T_e as well as ε_e. D_eff is the effective diffusion coefficient and describes the plasma diffusion. Its heuristically form, derived by [9] and [29], allows to properly describe the transition from ambipolar (in the bulk of plasma) to free diffusion (at its edge)

$\begin{array}{ll}{D}_{eff}=\frac{\alpha D_e+D_a}{1+\alpha }& (8)\end{array}$

with α = ν_iτ_M and τ_M is the Maxwell relaxation time given by τ_M = ϵ₀/(en_e(μ_i + μ_e)). The free electron and ambipolar diffusion coefficients are given by D_e = μ_ek_BT_e/e and D_a = μ_ik_BT_e/e, respectively. The electron mobility μ_e is related to the electron-neutral collision frequency by μ_e = e/(m_eν_m) and we assume that μ_i is kept constant and equal to μ_e/200. Finally regarding the recombination coefficient, r exhibits an electron temperature dependence which is here simply expressed as $r=\alpha \times 10^{-13}{(300/T_e)}^{-1/2}$ m⁻³ s⁻¹, with T_e being in Kelvin and α = 0.1 (according to [10]). The mean electron energy dependence of the rate coefficients, such as the effective ionization and the energy relaxation frequency, and the transport coefficient are obtained from solutions of the Boltzmann equation under a uniform field. The BOLSIG+ solver [30] is used for this purpose by employing the electron-neutral cross-section data set of Biagi in the LXCat database [31]. Figure 1 shows a typical trend of the reduced effective ionization frequency and the reduced energy relaxation frequency used in the model for air as a function of the mean electron energy for a given reduced wave frequency. In this paper, the gas heating effects on the plasma discharge behavior are neglected as we focus on the first stages of the streamer dynamics after microwave breakdown in air at atmospheric pressure (p = 760 Torr).

FIGURE 1

Figure 1. Reduced effective ionization frequency ν_eff/N (A) and reduced energy relaxation frequency ν_ε/N (B) for air as function of the mean electron energy ε_e. They are obtained at ω = 2π110 × 10⁹rad/s and p = 760 Torr from the BOLSIG+ solver using the Biagi input data.

Note that the presented plasma fluid model is the generalization of the simplified one employing the local effective field approximation which has been demonstrated to provide excellent results for microwave breakdown simulations under various conditions [10, 13, 32]. More details on their physical derivation and accuracy in conditions similar to those studied in this paper (high-frequency microwave fields, breakdown voltages, atmospheric pressure) can be found in [10, 11, 29].

3. Numerical Model

The time-domain Maxwell's equations coupled with the plasma momentum transfer Equation (1, 2, 3) are solved with a two-level iterative algorithm based on a finite element approximation and an explicit leapfrog time-stepping technique using a unique constant time step. The plasma fluid equations (5) are solved using a second-order central difference scheme explicitly integrated in time with an Euler scheme as presented in Arcese et al. [11]. In such a discretization, the terms on the right-hand side of Equations (5) are treated using an implicit time-stepping technique as proposed by Hagelaar and Kroesen [33] in order to avoid instabilities in the global resolution. Remark that the time step Δt_plasma related to Equations (5) is set equal to the electromagnetic wave period T_M that is much larger than the time step Δt related to the augmented Maxwell system (1, 2, 3), as we have separated the time scales of the microwave oscillations and the plasma evolution. The microwave-plasma coupling is hence described as follows: the microwave solution obtained during one wave cycle over the global computational domain (through the proposed multiscale method), consisting of two nested Cartesian grids having different sizes and levels of refinement, yields the time-average electron power absorption which is injected in the fluid plasma model to simulate the evolution of the plasma density over the next wave cycle in the patch domain, corresponding to the finest grid. The plasma density is then injected back into the microwave equations to update the power absorption. Note that, in this paper we limit ourselves to the case in which the plasma is only defined in the patch domain because of the standard approximation scheme employed for plasma equations. In the following sections, it is presented the patch-based resolution algorithm for the treatment of the microwave fields. The shared-memory parallelization using OpenMP is realized for the numerical solver.

3.1. Patch-Based Method for Solving the Coupled Maxwell's Equations and the Momentum Transfer Equation

3.1.1. Overview

The patch technique was first introduced by Glowinski et al. [19] to numerically solve elliptic problems with multiscale behavior using multiple levels of not necessarily nested grids. It is a finite element domain decomposition method with complete overlapping and is similar to a Chimera method used in Computational Fluid Dynamics [21]. Being compared to the Fast Adaptive Composite grid method [34–36] (as it is based on a multiplicative Schwarz algorithm) it offers much more flexibility since no conformity between the meshes at the different scales is required. Its iterative relaxed algorithm is straightforward. The problem is successively solved on a coarse mesh covering the whole computational domain and on a single or multiple patches with finer meshes where a better accuracy is needed. Successive corrections to the solution in the patches are thus calculated. A detailed analysis of this method and its convergence properties for linear elliptic problems is carried-out in Wagner [37] and a fast converging variant of the method has been presented [38]. Its extension to a class of non-linear elliptic problems is also studied in Brunet et al. [23]. Hereafter, we briefly introduce such a method before extending it to electromagnetic wave propagation problems.

Consider a linear elliptic problem on a domain Ω ⊂ ℝⁿ which weak formulation can be written

$\begin{array}{ll}\langle {L}(u)|\phi \rangle =\langle f|\phi \rangle \forall \phi \in H_0^1(\text{Ω}),& (9)\end{array}$

where ${L}(·)$ is a continuous, linear, symmetric, strongly elliptic operator (with H¹ being the standard Sobolev space) and f ∈ H⁻¹(Ω) (as boundary conditions we are simply setting u = 0 on ∂Ω). Suppose the solution u of the problem varies rapidly in a small sub-domain Λ ⊂ Ω, called patch (or in multiple sub-domains) and varies slowly in Ω \ Λ. We build a finite dimensional subspace $V_H\subset H_0^1(\Omega )$ associated to a coarse mesh on Ω and a subspace $V_h\subset H_0^1(\Omega )$ associated to a finer mesh on Λ. By setting then V_Hh = V_H + V_h, we approximate the solution $u\in H_0^1(\Omega )$ of the continuous problem (9) with u_Hh ∈ V_Hh which satisfies

$\begin{array}{ll}\langle {L}(u_{Hh})|\phi \rangle =\langle f|\phi \rangle \forall \phi \in V_{Hh}.& (10)\end{array}$

Note that such a problem is not trivial since it is generally not possible to determine a finite element basis of the total space V_Hh and hence to directly determine the solution u_Hh. Therefore, the aim of the patch technique is to evaluate u_Hh in an iterative manner only using the given basis V_H and V_h. The working principle of the algorithm is the following. It starts solving the problem on the coarse mesh in order to have an initial solution. At each iteration, it successively solves the residual of the problem, $R^n=f-{L}(u^n)$, on the patch mesh and the coarse one. The global solution is then updated using a relaxation parameter ω_r. The optimal choice of this parameter depends on the meshes and it allows to optimize the convergence rate of the algorithm [37]. The structure of the algorithm is presented in Algorithm (1).

ALGORITHM 1

Algorithm 1: Patch algorithm when using only a single patch.

Note that the global solution corresponds to the sum of the coarse one and the fine one that is u = u_H + u_h. The algorithm consists in evaluating two successive corrections to the solution in the patch and over the coarse domain. An extensive analysis has shown the strong dependence of the speed convergence of the algorithm (1) with respect to mainly the geometry of meshes and their mutual positions, hence with respect to an abstract angle between their respective finite element spaces, and the relaxation [19]. One can distinguish different types of grid configurations, from a boundary conforming case with structured grids to a boundary non-conforming case with unstructured grids, however for this study we only focus on nested Cartesian grids.

3.1.2. Application to the Electromagnetic Wave Propagation Problem

In order to extend the patch technique to our electromagnetic wave propagation problem we first need to employ a weak formulation for the discretization of Equations (1, 2, 3) and a time-stepping technique for time integration. It has been proved that the standard Yee scheme can be reformulated in terms of finite elements on orthogonal meshes when using a mass-lump technique [39]. Therefore, in this paper we make use of this variational formulation coupled with an explicit leapfrog time-stepping technique using a single constant time step. In the first part of this section we give a detailed presentation of the above-stated formulation, while in the last part we derive the related patch algorithm. For the sake of presentation, we rewrite our system of equations in a concise form (the tilde character is neglected to lighten the notation) assuming a perfect conductor as boundary of Ω

$\begin{array}{r}{ϵ}_0\frac{\partial \text{E}}{\partial t}=\nabla \times \text{H}+e{n}_e{\text{v}}_e\text{ in Ω }\times ]0,T[\\ {\mu }_0\frac{\partial \text{H}}{\partial t}=-\nabla \times \text{E}\text{ in Ω }\times ]0,T[\\ \frac{\partial {\text{v}}_e}{\partial t}=-\frac{e\text{E}}{m_e}-{\nu }_m{\text{v}}_e\text{ in Ω }\times ]0,T[\\ \text{E}\times \text{n}=0\text{ on }\partial \text{Ω }\times ]0,T[,& (11)\end{array}$

where t ∈]0, T[ and the initial conditions are such that

$\begin{array}{ll}\text{E}(\text{x},0)={\text{E}}_0(\text{x}),\text{H}(\text{x},0)={\text{H}}_0(\text{x}),{\text{v}}_e(\text{x},0)=0\forall \text{x}\in \text{Ω}.& (12)\end{array}$

Note that, as consequence of the boundary condition on the electric field and the third equation in (11), the boundary condition on the mean velocity is satisfied at any time as well, v_e × n = 0.

3.1.2.1. Mixed finite element formulation

In the first step of our mixed finite element approximation we define a variational formulation of system (11). Of course, several choices are possible, but we exactly opt for a formulation based on the following functional spaces (in 3D):

$\begin{array}{l}(\text{E},{\text{v}}_e)\in H_0(\text{c}\text{u}\text{r}\text{l},\text{ Ω})\times H_0(\text{c}\text{u}\text{r}\text{l},\text{ Ω})\\ \text{H}\in {\left[L^2(\text{Ω})\right]}^3& (13)\end{array}$

where H₀(curl, Ω) = {ϕ ∈ H(curl, Ω) such that ϕ × n = 0 on ∂Ω} is a closed subspace of the Hilbert space

$\begin{array}{l}H(\text{c}\text{u}\text{r}\text{l},\text{ Ω})=\left\{\varphi \in {\left[L^2(\text{Ω})\right]}^3\text{such that }\nabla \times \varphi \in {\left[L^2(\text{Ω})\right]}^3\right\},\end{array}$

with L²(Ω) being the space of functions whose square is integrable over Ω and n being the outward normal unit. Such a particular choice is justified by the fact that the electric field assumes a leading role in microwave plasma problems as it drives the plasma dynamics unlike the magnetic field, therefore more regularity on the E-field is required. Indeed, the above formulation is such that E keeps its physical character (tangential continuity) by belonging to H₀(curl, Ω). Instead, a weakened regularity is provided on H by seeking it in [L²(Ω)]³ which is larger than H₀(curl, Ω). Moreover, as the mean electron velocity is directly linked to the electric field via an ordinary differential equation in system (11), v_e is consistently sought in the same functional space of E.

By multiplying the first and third Equation (11) by ϕ ∈ H₀(curl, Ω) and $\widetilde{\text{ϕ}}\in H_0(\text{curl},\Omega )$, respectively, and the second one by ψ ∈ [L²(Ω)]³ and then integrating by parts the stiffness integral of the first equation (which corresponds to ∇ × H), we obtain

$\begin{array}{l}{ϵ}_0\frac{\partial }{\partial t}{\displaystyle {\int }_{\text{Ω}}}\text{E}·\varphi d\text{x}={\displaystyle {\int }_{\text{Ω}}}\text{H}·(\nabla \times \varphi )d\text{x}+{\displaystyle {\int }_{\partial \text{Ω}}}\left[\text{n}\times (\text{H}\times \text{n})\right]·\\ (\varphi \times \text{n})d\gamma +e{\displaystyle {\int }_{\text{Ω}}}{n}_e{\text{v}}_e·\varphi d\text{x}\\ {\mu }_0\frac{\partial }{\partial t}{\displaystyle {\int }_{\text{Ω}}}\text{H}·\psi d\text{x}=-{\displaystyle {\int }_{\text{Ω}}}(\nabla \times \text{E})·\psi d\text{x}\\ \frac{\partial }{\partial t}{\displaystyle {\int }_{\text{Ω}}}{\text{v}}_e·\tilde{\varphi }d\text{x}=-\frac{e}{m_e}{\displaystyle {\int }_{\text{Ω}}}\text{E}·\tilde{\varphi }d\text{x}-{\displaystyle {\int }_{\text{Ω}}}{\nu }_m{\text{v}}_e·\tilde{\varphi }d\text{x}.& (14)\end{array}$

According to the choice of the functional space H₀(curl, Ω) (which takes into account a medium with a perfectly conducting boundary condition, i.e., E × n = 0 for the E-field) the boundary term in the first Equation (14) vanishes since [n × (H × n)] · (ϕ × n) = H · (ϕ × n). As a result we obtain the following variational problem:

Find $(\text{E}(·,t),\text{H}(·,t),{\text{v}}_e(·,t))\in H_0(\text{c}\text{u}\text{r}\text{l},\Omega )\times {\left[L^2(\Omega )\right]}^3\times H_0(\text{c}\text{u}\text{r}\text{l},\Omega )$, ∀t ∈]0, T[ such that

$\begin{array}{l}{ϵ}_0\frac{\partial }{\partial t}{\displaystyle {\int }_{\text{Ω}}}\text{E}·\varphi d\text{x}={\displaystyle {\int }_{\text{Ω}}}\text{H}·(\nabla \times \varphi )d\text{x}+e{\displaystyle {\int }_{\text{Ω}}}{n}_e{\text{v}}_e·\varphi d\text{x}\\ \forall \varphi \in H_0(\text{c}\text{u}\text{r}\text{l},\text{ Ω}),\\ {\mu }_0\frac{\partial }{\partial t}{\displaystyle {\int }_{\text{Ω}}}\text{H}·\psi d\text{x}=-{\displaystyle {\int }_{\text{Ω}}}(\nabla \times \text{E})·\psi d\text{x}\forall \psi \in {\left[L^2(\text{Ω})\right]}^3,\\ \frac{\partial }{\partial t}{\displaystyle {\int }_{\text{Ω}}}{\text{v}}_e·\tilde{\varphi }d\text{x}=-\frac{e}{m_e}{\displaystyle {\int }_{\text{Ω}}}\text{E}·\tilde{\varphi }d\text{x}-{\displaystyle {\int }_{\text{Ω}}}{\nu }_m{\text{v}}_e·\tilde{\varphi }d\text{x}\\ \forall \tilde{\varphi }\in H_0(\text{curl},\text{ Ω}),& (15)\end{array}$

which are subject to the initial conditions (12). Let us remark that the above formulation is a slight modification (an additional ordinary differential equation on the mean velocity appears as the plasma density is non-zero) of the one analyzed by Monk [40] for the standard Maxwell system where no additional boundary condition on the magnetic field is necessary.

In a more compact form, the problem (15) appears as

Find $\text{u}(·,t)\in \text{V}\equiv H_0(\text{c}\text{u}\text{r}\text{l},\Omega )\times {\left[L^2(\Omega )\right]}^3\times H_0(\text{c}\text{u}\text{r}\text{l},\Omega )$, ∀t ∈]0, T[ such that

$\begin{array}{ll}<{M}\frac{d\text{u}}{dt}|\phi >=<{A}\text{u}|\tilde{{A}}\phi >+<\tilde{{M}}\text{u}|\phi >\forall \phi \in \text{V}& (16)\end{array}$

with the matrices ${M}$, $\tilde{{M}}$, ${A}$, and ${Ã}$ being

$\begin{array}{c}{M}=\left(\begin{array}{ccc}{ϵ}_0& 0& 0\\ 0& {\mu }_0& 0\\ 0& 0& 1\\ \end{array}\right);\tilde{{M}}=\left(\begin{array}{ccc}0& 0& e{n}_e\\ 0& 0& 0\\ -\frac{e}{m_e}& 0& -{\nu }_m\\ \end{array}\right);\\ {A}=\left(\begin{array}{ccc}0& 1& 0\\ -{\nabla }_{\times }& 0& 0\\ 0& 0& 0\\ \end{array}\right);\tilde{{A}}=\left(\begin{array}{ccc}{\nabla }_{\times }& 0& 0\\ 0& 1& 0\\ 0& 0& 0\\ \end{array}\right),\end{array}$

the vectors u and φ being

$\begin{array}{c}\text{u}=\left(\begin{array}{c}\text{E}\\ \text{H}\\ {\text{v}}_e\\ \end{array}\right);\phi =\left(\begin{array}{c}\varphi \\ \psi \\ \tilde{\varphi }\\ \end{array}\right),\end{array}$

and where < ·|·> here denotes the standard scalar product on L²(Ω).

Now, we choose conforming finite elements based on the first family of Nedelec's edge elements [25] in order to build the finite dimensional approximation space V_h such that V_h ⊂ V. Particularly, we focus on mass-lumped spectral elements having a cubic shape in which the basis function space is defined as a tensor product of 1D polynomial spaces. Considering a cubic mesh of Ω, ${{T}}_h$, in which the characteristic dimension of each element is h, the approximation space is thus identified by ${\text{V}}_h={\text{S}}_{h,0}^r\times {\text{D}}_h^r\times {\text{S}}_{h,0}^r$ where, for this specific case

$\begin{array}{l}{\text{S}}_{h,0}^r=\{{\varphi }_h\in H_0(\text{curl},\text{ Ω})\text{ s}.\text{t}.{\varphi }_h{|}_K\in Q_{r-1,r,r}\times Q_{r,r-1,r}\\ \times Q_{r,r,r-1},\forall K\in {{T}}_h\}\\ {\text{D}}_h^r=\{{\psi }_h\in {\left[L^2(\text{Ω})\right]}^3\text{ s}.\text{t}.{\psi }_h{|}_K\in Q_{r,r-1,r-1}\times Q_{r-1,r,r-1}\\ \times Q_{r-1,r-1,r},\forall K\in {{T}}_h\}& (17)\end{array}$

with Q_{l, m, n} being the space of polynomials of maximum degree l in x, m in y and n in z defined on the 3D element $K\in {{T}}_h$. The semi-discrete variational formulation can be then written as

Find u_h(·, t) ∈ V_h, ∀t ∈]0, T[ such that

$\begin{array}{ll}<{M}\frac{d{\text{u}}_h}{dt}|\phi >=<{A}{\text{u}}_h|\tilde{{A}}\phi >+<\tilde{{M}}{\text{u}}_h|\phi >\forall \phi \in {\text{V}}_h& (18)\end{array}$

$\begin{array}{l}{\text{ u}}_h(\text{x},0)={\text{u}}_h^0(\text{x}),\end{array}$

where the approximated solution u_h reads as

$\begin{array}{ll}{\text{u}}_h(\text{x},t)={\displaystyle \sum _{\begin{array}{l}i=1\end{array}}^N}{u}_h^i(t){\phi }^i(\text{x}),& (19)\end{array}$

with φⁱ and $u_h^i$ being the vector-valued basis function and the degree of freedom for u_h, respectively, and N being the dimension of V_h. Remark that we use a general notation in (19), but actually the degrees of freedom for E_h, H_h, and v_{e, h} are located on the mesh elements in a nodal way accounting for the tensor character of the polynomial spaces (17) [39].

By substituting (19) in Equation (18), we can express the problem in a matrix form in which the mass (M_h, ${\tilde{M}}_h$) and the stiffness (A_h) matrices are identified

$\begin{array}{ll}{M}_h\frac{d{U}_h}{dt}=(A_h+{\tilde{M}}_h)U_h,& (20)\end{array}$

and where $U_{h_i}(t)=u_h^i(t)$, M_h, ${\tilde{M}}_h$ and A_h are matrices with coefficients, respectively

$\begin{array}{ll}{M}_{h_{i,j}}=<{M}{\phi }^i|{\phi }^j>;{\tilde{M}}_{h_{i,j}}=<\tilde{{M}}{\phi }^i|{\phi }^j>;A_{h_{i,j}}=<{A}{\phi }^i|\tilde{{A}}{\phi }^j>.& (21)\end{array}$

Previous works have proved that the proposed space approximation exactly provides the Yee scheme on orthogonal meshes with hexahedral elements [39]. This is indeed possible by choosing the degree of polynomials r = 1 and making use of a mass-lumping technique (which means that the above mass matrices reduce to block-diagonal matrices on the specific elements) as well as an explicit leapfrog time-stepping scheme. In what follows, we benefit of such a space approximation in order to derive the two-level resolution algorithm.

3.1.2.2. Approximation in time

We recognize that the semi-discrete problem (18) has the same form of (10) however, the bilinear form $<{L}(·)|\phi >$ consists here of a temporal part and a spatial one, particularly

$\begin{array}{ll}<{L}({\text{u}}_h)|\phi >=<{M}\frac{d{\text{u}}_h}{dt}|\phi >-<{A}{\text{u}}_h|\tilde{{A}}\phi >-<\tilde{{M}}{\text{u}}_h|\phi >.& (22)\end{array}$

In this paper, we choose a second-order accurate explicit leapfrog time-integration technique which is the same provided by the standard Yee scheme. In particular, the electromagnetic and mean electron velocity fields are staggered so that the E-field is updated midway during each time-step between successive updates of the H-field and v_e-field in order to provide no dissipation in the numerical wave propagation. By denoting t_n = nΔt with n ∈ ℕ and Δt the time step chosen for time discretization, the system (18) thus reads in a matrix form as

$\begin{array}{c}\left[\begin{array}{ccc}{M}_h^{{ϵ}_0}& 0& 0\\ -\text{Δ}t{A}_h^E& M_h^{{\mu }_0}& 0\\ -\text{Δ}t{M}_h^E& 0& M_h^{v_e}-\frac{\text{Δ}t}{2}{M}_h^{{\nu }_m}\end{array}\right]\\ U_h^{n+1}=\left[\begin{array}{ccc}{M}_h^{{ϵ}_0}& \text{Δ}t{A}_h^H& \text{Δ}t{M}_h^{v_e/E}\\ 0& M_h^{{\mu }_0}& 0\\ 0& 0& M_h^{v_e}+\frac{\text{Δ}t}{2}{M}_h^{{\nu }_m}\end{array}\right]U_h^n& (23)\end{array}$

where $U_h^n$ and $U_h^{n+1}$ are the successive time approximations of U_h(t). Let us underline that such a time-stepping notation translates into the leapfrog time arrangement of the electromagnetic and the velocity field as follows: at $U_h^n$ it is associated the quantities $E_h^n$, $H_h^{n+1/2}$, and $v_{e,h}^{n+1/2}$ estimated at times t_n, t_n+1/2 and t_n+1/2, respectively, whereas $U_h^{n+1}$ designates the quantities $E_h^{n+1}$, $H_h^{n+3/2}$, and $v_{e,h}^{n+3/2}$ estimated at times t_n+1, t_n+3/2, and t_n+3/2, respectively.

The stiffness and mass (which have all a block-diagonal structure according to the above space approximation) matrices appearing in the scheme are defined by following the same notation of the previous section, hence their coefficients read as

$\begin{array}{l}{M}_{h_{i,j}}^{{ϵ}_0}=<{ϵ}_0{\varphi }^i|{\varphi }^j>;M_{h_{i,j}}^{{\mu }_0}=<{\mu }_0{\psi }^i|{\psi }^j>;M_{h_{i,j}}^{v_e}=<{\tilde{\varphi }}^i|{\tilde{\varphi }}^j>;\\ M_{h_{i,j}}^E=-<\left(\frac{e}{m_e}\right){\varphi }^i|{\tilde{\varphi }}^j>;M_{h_{i,j}}^{v_e/E}=<(e{n}_e){\tilde{\varphi }}^i|{\varphi }^j>;\\ M_{h_{i,j}}^{{\nu }_m}=-<{\nu }_m{\tilde{\varphi }}^i|{\tilde{\varphi }}^j>;\\ A_{h_{i,j}}^E=-<\nabla \times {\varphi }^i|{\psi }^j>;A_{h_{i,j}}^H=<{\psi }^i|\nabla \times {\varphi }^j>,\end{array}$

with $A_h^E=-{\left(A_h^H\right)}^T$. Note that the coupling between E and v_e is integrated explicitly using the same leapfrog arrangement while the source term of the equation on v_e, corresponding to the neutral collision frequency contribution, is treated implicitly using a central discretization in order to alleviate the global stability condition. To figure out how electromagnetic and velocity fields are staggered both in space and time in the proposed scheme on a Cartesian mesh, we rewrite (23) for the simplified 2D transverse-electric (TE) case:

$\begin{array}{l}{ϵ}_0{E}_{x,(i+1/2,j)}^{n+1}={ϵ}_0{E}_{x,(i+1/2,j)}^n+\frac{\text{Δ}t}{\text{Δ}y}{\displaystyle (H_{z,(i+1/2,j+1/2)}^{n+1/2}-H_{z,(i+1/2,j-1/2)}^{n+1/2})}\\ +e\text{Δ}t{n}_{e,(i+1/2,j+1/2)}{v}_{e,x,(i+1/2,j)}^{n+1/2},\\ {ϵ}_0{E}_{y,(i,j+1/2)}^{n+1}={ϵ}_0{E}_{y,(i,j+1/2)}^n-\frac{\text{Δ}t}{\text{Δ}x}{\displaystyle (H_{z,(i+1/2,j+1/2)}^{n+1/2}-H_{z,(i-1/2,j+1/2)}^{n+1/2})}\\ +e\text{Δ}t{n}_{e,(i+1/2,j+1/2)}{v}_{e,y,(i,j+1/2)}^{n+1/2},\\ {\mu }_0{H}_{z,(i+1/2,j+1/2)}^{n+3/2}={\mu }_0{H}_{z,(i+1/2,j+1/2)}^{n+1/2}-\text{Δ}t{\displaystyle [{\displaystyle \left(\frac{E_{y,(i+1,j+1/2)}^{n+1}-E_{y,(i,j+1/2)}^{n+1}}{\text{Δ}x}\right)}}\\ -{\displaystyle \left(\frac{E_{x,(i+1/2,j+1)}^{n+1}-E_{x,(i+1/2,j)}^{n+1}}{\text{Δ}y}\right)]},\\ v_{e,x,(i+1/2,j)}^{n+3/2}(1+\frac{\text{Δ}t}{2}{\nu }_m)=v_{e,x,(i+1/2,j)}^{n+1/2}-\frac{e\text{Δ}t}{m_e}{E}_{x,(i+1/2,j)}^{n+1}\\ -\frac{\text{Δ}t}{2}{\nu }_m{v}_{e,x,(i+1/2,j)}^{n+1/2},\\ v_{e,y,(i,j+1/2)}^{n+3/2}(1+\frac{\text{Δ}t}{2}{\nu }_m)=v_{e,y,(i,j+1/2)}^{n+1/2}-\frac{e\text{Δ}t}{m_e}{E}_{y,(i,j+1/2)}^{n+1}\\ -\frac{\text{Δ}t}{2}{\nu }_m{v}_{e,y,(i,j+1/2)}^{n+1/2},& (24)\end{array}$

where Δx and Δy are the spatial steps chosen for the spatial discretization such that the grid coordinates are x_i = iΔx and y_j = jΔy with (i, j) ∈ ℕ. Such a time-stepping scheme is generally more practical compared to an implicit version, however, it mandates an upper bound on the time-step to ensure numerical stability. Through a discrete energy estimate, based on the analysis of [14], we obtain the stability condition related to our scheme. The constraint on the time step is the following:

$\begin{array}{ll}\text{Δ}t\le \frac{1}{2}min{\displaystyle ({\displaystyle (c\sqrt{\frac{1}{\text{Δ}{x}^2}+\frac{1}{\text{Δ}{y}^2}}{)}^{-1}},{\displaystyle ({\omega }_p{)}^{-1}})}& (25)\end{array}$

with ${\omega }_p=\underset{i,j}{max}(e^2{n}_{e,(i+1/2,j+1/2)})/({ϵ}_0{m}_e){}^{1/2}$ the plasma frequency. One can remark that in our case the value of the plasma frequency (even for plasma densities reaching values of 10²² m⁻³) is such that its contribution on the stability condition is negligible with respect to the first argument of the minimum in (25) which is exactly the condition provided by the Yee scheme for the Maxwell's equations having no electron current terms.

3.1.2.3. Derivation of the patch method for the Yee scheme

At this step, we derive the patch method for the introduced space-time approximation. The system (23) is written in the following concise form

$\begin{array}{ll}(I_h-\text{Δ}t{R}_h^{+}-\text{Δ}t{L}_h^{+})U_h^{n+1}=(I_h+\text{Δ}t{R}_h^{-}+\text{Δ}t{L}_h^{-})U_h^n& (26)\end{array}$

with

$\begin{array}{l}{I}_h=\left[\begin{array}{ccc}{M}_h^{{ϵ}_0}& 0& 0\\ 0& M_h^{{\mu }_0}& 0\\ 0& 0& M_h^{v_e}\end{array}\right];R_h^{+}=\left[\begin{array}{ccc}0& 0& 0\\ A_h^E& 0& 0\\ 0& 0& 0\end{array}\right];\\ R_h^{-}=\left[\begin{array}{ccc}0& A_h^H& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right];\\ L_h^{+}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ M_h^E& 0& \frac{1}{2}{M}_h^{{\nu }_m}\end{array}\right];L_h^{-}=\left[\begin{array}{ccc}0& 0& M_h^{v_e/E}\\ 0& 0& 0\\ 0& 0& \frac{1}{2}{M}_h^{{\nu }_m}\end{array}\right]& (27)\end{array}$

being the matrices associated to the approximation space V_h.

By introducing now two finite dimensional spaces V_H and V_h associated to a coarse mesh on Ω and a fine mesh on Λ, respectively, such as V_H + V_h = V_Hh ⊂ V, as defined in section 3.1.1, the numerical solution U_Hh of problem (26) defined on V_Hh is approximated as U_h + U_H = U_Hh ∈ V_Hh through the iterative Algorithm (2), where n and p are the indexes concerning the time and the iterative advancement, respectively.

ALGORITHM 2

Algorithm 2: Patch algorithm applied to the weak formulation of the Yee scheme when using only a single patch.

The convergence is checked on the L²-norm error of the discrete energy between two successive iterations. Let us stress that the above algorithm is for a two-level configuration where only a single patch is introduced into the global domain as the cases studied in this paper (see numerical results in section 4). Its extension to multiple patches is straightforward [19] and will be the object of future works. Figure 2 illustrates how electromagnetic and velocity fields and plasma variables are spatially located in our computational domain.

FIGURE 2

Figure 2. Simplified plan view of a 2D computational domain containing a single patch. The Cartesian grids are nested and the fine one (which is inside the central coarse cell) has a spatial refinement factor of 2:1. E- and v_e-field are represented by circles, H-field is represented by squares and plasma variables (i.e., n_e, ε_e, ν_eff, etc.) are represented by red crosses. Full and empty circles (or squares) denote the coarse and the fine grid, respectively.

Note that in the above algorithm some transformation matrices, i.e., the matrices having the double subscript Hh (or hH), appear on the right-hand side of equations of the problems defined both on the coarse and the fine domain. These enable the corrections to the solution defined on a grid by projecting the solution defined and just calculated on the other grid. They mathematically correspond to the mixed term scalar products, wherein the basis functions of both fine and coarse subspaces appear. For instance, looking at the term I_Hh, the coefficients of the related mass matrix $M_{Hh}^{{ϵ}_0}$ read

$\begin{array}{l}{M}_{H{h}_{i,j}}^{{ϵ}_0}=<{ϵ}_0{\varphi }_H^i|{\varphi }_h^j>,\end{array}$

and the same, but in the opposite way, holds for the term I_hH. Remark that all vector-valued basis functions belonging to the fine subspace V_h are by definition prolonged to the coarse subspace V_H hence, they vanish on the boundaries of the patch domain.

The relaxation parameter ω_r provides a particular flexibility to the algorithm with respect to the conformity between the coarse and the fine discretization [37]. This feature makes the patch method rather suitable for adaptive mesh refinements. Indeed, in case of refinement/derefinement (conforming or not) one can increase/decrease the dimension of the related approximation spaces (or change their mutual intersection as well) and hence calculating/deleting the related coupling terms while optimizing the convergence by means of relaxation. However, in this paper we focus on conforming grid constellations setting the relaxation parameter to a constant value, ω_r = 1. As a concluding remark, we highlight that the standard Yee scheme preserves a discrete energy when no electron currents are present in the Maxwell's equations and for perfectly conducting boundary conditions. As such, we expect to hold the same property in the patch version of the scheme (since a weak formulation is employed for discretization) which guarantees the stability (see numerical results in section 4.1.2).

3.2. Computational Complexity Analysis of the Patch Technique

In this section, we analyze the efficiency from a performance point of view of the proposed patch-based scheme for the treatment of the electromagnetic fields. A theoretical estimate of the expected performance gain with respect to the standard Yee scheme is given.

Let us consider a computational domain having a characteristic length L which is discretized into a uniform mesh of N points in order to properly resolve the finest characteristic length of the wave propagation problem. The grid spatial step is identified with h. The computational complexity of the standard Yee scheme in solving on such a grid the D-dimension problem per time step is O(N^D). Now, by focusing on a computational domain containing a local patch having a characteristic length L_p, we define the total number of points in the patch grid as N_f = N/α, with α = L/L_p, in order to ensure the same spatial step h of the previous uniformly refined mesh, and the total number of points in the coarse grid as N_c = N/β, with β = H/h and H being the coarse spatial step. We also define the number of iterations needed for convergence by the patch algorithm per time step as k. According to this notation, the computational cost of the patch-based scheme is O(kN^D(1/α^D + 1/β^D)) and consequently, the performance gain between the two approaches results

$\begin{array}{ll}\eta ={\displaystyle [k{\displaystyle (\frac{1}{{\alpha }^D}+\frac{1}{{\beta }^D})}{]}^{-1}}.& (28)\end{array}$

One can readily remark that the patch version of the Yee scheme can exhibit important gains with respect to the standard one without patches at a given k for large values of the parameters α and β (note that α > 1 and β > 1 generally). Of course, the number of iterations for convergence k depends on the characteristics of the patch and coarse meshes and on the physical problem as well. However, for the 1D and 2D cases analyzed in this paper, k is on the order of few tens and appears to slightly vary with the patch refinement (see the numerical results in section 4). We also recognize that this gain is more important at higher dimensions, i.e., in 2D and even more in 3D, as α and β are to the power of the problem dimension D. Therefore, given a D-dimension problem significant gains of performance can be reached by choosing α with the same order of magnitude of β and much larger than unity. These configurations exactly correspond to problems presenting strong multiscale behaviors as it occurs in our case concerning the dynamics of microwave streamers. As a practical example: in our case the main lengths characterizing the physics are the wavelength of the high-frequency microwave, that is approximately few millimeters (e.g., 2.7 mm for a wave frequency of 110 GHz), and the propagation front length of plasma, which at atmospheric pressure and in conditions analyzed in this paper (see the 2D numerical results in section 4) is approximately 10 μm. By setting L = 2.7 mm and supposing that the plasma front is sufficiently contained in the patch, hence setting for instance L_p = 50 μm, we choose the ratio of the coarse and the fine spatial grid equal to the ratio of these characteristic lengths such as α = β = 54. In the 2D case and for a given k ranging from 5 to 20, we thus expect interesting theoretical performance gains varying from one to two orders of magnitude, namely from 72.9 to 291.6.

4. Numerical Results

In this section, we assess the efficiency of the proposed multiscale method for the wave propagation problem in microwave discharge plasma. Simplified settings are first analyzed in 1D, then 2D simulations of a more realistic case, such as the microwave breakdown in air, are performed. In both cases, the computational domain consists of two nested Cartesian grids: a coarse grid covering the global domain and a fine one covering a small region of that domain wherein the plasma properties are evaluated.

4.1. One-Dimensional Simulations

4.1.1. Interpretation of the Patch Results

As explained above, in the patch technique the global numerical solution is given by the sum of the solution evaluated on the coarse domain and the one evaluated in the patch. Practically, the recovering of the global solution is possible through a simply interpolation of the coarse/fine solution or both on a given grid. One can remark that by definition the fine solution is defined over the global domain (as the coarse solution), but it is zero outside the patch. In order to illustrate how these solutions appear in a classic numerical experiment, we present the case of propagation of a modulated Gaussian pulse through a static plasma behaving as a conductor (the plasma density is overcritical and follows a Gaussian profile with a maximum density of 10²² m⁻³). The patch sufficiently contains the plasma density profile and it is located in the center of the coarse domain which has perfect boundary conditions. In this case, the plasma dynamics is not resolved as we consider a constant density plasma structure. In Figure 3, one can easily check the solution of the electric field (in black solid line) on the global computational domain obtained by the fine and coarse solutions (in red dashed and dash-dotted blue lines, respectively).

FIGURE 3

Figure 3. Time evolution of the electric field concerning the 1D propagation of a modulated Gaussian pulse through a static plasma [particularly, the field at different times from (A)–(D)]. The Gaussian pulse width is large enough to see about 3 cycles of modulation. The plasma density is a Gaussian with maximum 10²² m⁻³ (with standard deviation of 10⁻⁵ m) in the center of patch identified by two vertical lines in the figures. The patch is about 13 times smaller than the coarse domain and the refinement ratio between coarse and fine grids is equal to H/h = 10, with h and H being the spatial fine and coarse steps, respectively. The global solution of the electric field is represented with a black solid line, the solutions calculated on the coarse and fine grids are in blue dash-dotted and red dashed lines, respectively.

4.1.2. Numerical Analysis

It is interesting to investigate the quality of the numerical solution when introducing a well-suited local patch having different refinements. For doing this, we consider a simplified version of the Maxwell problem (15) with a stiff source term. Specifically, the test case consists of an 1D standing wave formation (the wave frequency is f = 110 GHz) in a cavity, with perfectly conducting walls, wherein a static plasma with overcritical density (the plasma density follows a Gaussian profile with maximum $n_0=5\times 10^{22}$ m⁻³ and standard deviation σ = 4 × 10⁻⁵ m) is introduced. For comparisons with an "exact" solution we simplify the problem (15) by considering only the Maxwell's equations with a current density term, within the mean electron velocity is simply set as v_e = μ_eE and no longer governed by the presented ordinary differential equation. An ad hoc source term, consistent with the standing wave solution, is thus introduced on the right-hand side of these simplified equations. In particular, the source term related to the E-field equation (the one related to the H-field equation is set equal to zero) in 1D reads as

$\begin{array}{ll}{f}_E(x,t)=-e{n}_0{\mu }_e{E}_{0}exp\left(-\frac{{(x-x_0)}^2}{2{\sigma }^2}\right)sin(2\pi k_{x}x)cos(2\pi ft),& (29)\end{array}$

where E₀ = 1 V/m, k_x = f/c, parameters e and μ_e are defined as in section (2) and x₀ is the mean of the plasma Gaussian distribution such that the exact solution to the problem is given by E = E₀ sin (2π k_xx) cos (2πft) and H = ϵ₀(f/k_x) E₀ cos (2πk_xx) sin (2πft). The peak of the plasma density profile is located in the first quarter of the computational domain, i.e., on the first wave anti-node, where a local patch is applied in order to get a better precision of the solution. The patch size is kept constant and 10 times smaller than the size of the coarse domain. Figure 4 gives a schematic representation of the computational domain and the 1D standing wave solution. Because of the stiffness of the introduced source term, we expect that, by progressively refining the mesh (locally or uniformly), the numerical solution tends to the exact one.

FIGURE 4

Figure 4. Electric field of the considered 1D standing wave having the frequency equal to 110 GHz (A). Example of the 1D computational domain (with perfect electrical conductor, PEC, boundary conditions) consisting of two nested Cartesian grids (B). The patch is 10 times smaller than the coarse domain having the size equal to the wavelength λ, L = λ = 2.725 × 10⁻³ m.

Figure 5 reports the relative error of the calculated electric field E_Hh (on the global domain) to the exact one E for three different coarse spatial steps H, H = λ/80, λ/160, λ/320 (with λ being the wavelength of the standing wave of Figure 4), and for different levels of patch refinement, H/h = 1, 2, 4, 8, 16, 32 (with h being the spatial fine step). Firstly, we remark that for a given, constant coarse spatial step the global error reduces with patch refinement, in particular, it reduces by approximately two orders of magnitude for an important refinement H/h = 32 (e.g., for H = λ/80, the error trend is represented with the red dashed line in Figure 5). This improvement of the solution accuracy with patch refinement is however bounded for a fixed H, as the error in the domain with coarse discretization dominates and the global error stagnates. In such a case, a better global precision can be obtained through a further refinement of the coarse grid. Figure 5 shows the global improvement of the solution accuracy between the cases with H = λ/80 (red dashed line), H = λ/160 (blue dashed line) and H = λ/320 (black dashed line). Note that the error when using a patch, with a given refinement, is considerably the same compared to the case without a patch, presenting the same refinement everywhere. This can be verified in Figure 5, by, e.g., comparing the error obtained in the case without patch, represented with the black circle (where H = λ/320 and H/h = 1) and the one related to the case with patch, represented with the blue cross (where H = λ/160 and H/h = 2) or with the red diamond (where H = λ/80 and H/h = 4). From these results, it is evident that the local patch shows its usefulness as long as the solution is well approached by the coarse grid. Furthermore, we point out the that the presence of the patch does not affect the consistency of the proposed finite element scheme. The latter indeed preserves the convergence order of the Yee finite-difference type scheme, which is second-order accurate both in space and time, for several patch refinements (the rate of convergence two for the L²-norm related to the electric field can be checked in Figure 5). For all presented cases, we have numerically assessed the algorithm stability through the estimation of the discrete electromagnetic energy.

FIGURE 5

Figure 5. The L²-norm relative error of the global numerical solution of the electric field E_Hh to the exact electric field E, $\left|\right|E_{Hh}-E|{|}_{L^2}/|\left|E\right|{|}_{L^2}$, as function of the fine and the coarse spatial step, h and H, respectively. The refinement level of the coarse grid varies between three values: H = λ/80, λ/160, λ/320. Concerning the patch, the refinement level varies from H/h = 1 to H/h = 32. The computational domain and the patch size are the same of the case in Figure 4.

Concerning the theoretical performance gain, a slightly advantage is expected when using the patch method compared to the standard one (without patch and with the same level of refinement everywhere) for important refinements (see Table 1). This is because the gain is limited by the coarse-fine characteristic lengths ratio α as the analysis in section 3.2 suggests. In higher dimension, the performance gain is more important at lower refinement ratios (see 2D numerical results in section 4.2).

TABLE 1

Table 1. The performance gain between the patch method and the standard one (without patch and with the same level of refinement everywhere) for the case of Figure 5 (black dashed line, H = λ/320) estimated using formula (28), with k being the average number of iterations needed for convergence in the patch algorithm per wave cycle (k_max and k_min are the maximum and the minimum number of iterations, respectively), α = 10, β = 2, 4, 8, 16, 32 and D = 1.

4.2. Two-Dimensional Simulations of Microwave Breakdown in Air

4.2.1. Simulation Conditions

Considering the 2D case for the numerical model presented in section 3, we simulate the formation of a microwave streamer at the antinode of a standing wave created by the intersection of two identical waves having a TE polarization with opposed wave vectors. These simulations correspond to the experiments of [41, 42] in which two-mirror microwave resonator has been used to focus the electromagnetic power in a small volume. A plasma discharge has been initiated because of the larger electric field with respect to the critical one for breakdown, and propagated in the direction of the electric field's polarization. The mechanism responsible of the plasma evolution, which is a combination of diffusion and strong ionization at the streamer tips due to the local enhancement of the electric field, has been demonstrated by detailed numerical experiments [32]. The authors have used an explicit formulation of the FDTD method, i.e., the standard Yee scheme, for the related plasma-EM (electromagnetic waves) model. An implicit formulation (i.e. ADI) for solving the Maxwell's equations in the same model has been proposed as well to reduce the huge computational cost required by the explicit formulation in 2D [14]. The configuration of our simulation is similar to the above-mentioned numerical experiments and can be summarized as follows:

• TE mode y-polarized plane waves with opposite x-directed wave vectors of amplitude E₀ = 2.5 MV/m and frequency f = 110 GHz each, forming a stationary wave having a maximum rms (root mean squared) field E_rms = 3.5 MV/m (which is larger than the critical field for breakdown in air at atmospheric pressure, E_c = 2.5 MV/m),

• atmospheric pressure, p = 760 Torr,

• initial plasma density of a Gaussian profile with maximum of 10¹⁵ m⁻³ in the center of the simulation domain and standard deviation of 6 × 10⁻⁵ m,

• computational domain, having dimension λ × λ, is divided into a scattered field (SF) and a total field (TF) domain in order to inject the waves. A local patch, having a fixed size such that ${\alpha }^2={(L/L_p)}^2=10.07$ (the patch size is kept constant during the simulation), is introduced into the domain and centered on the peak of the initial plasma density profile, which is set to zero in the coarse domain outside the patch. Convolutional perfectly matching layer (CPML) absorbing boundary conditions are thus implemented in order to minimize reflections into the global domain (see Figure 6A),

• nested Cartesian grids, as shown in Figure 6B, having the coarse spatial step H = λ/100 and the patch steps h = H/4 and h = H/10,

• stopping criterion on the discrete electromagnetic energy (in the iterative patch algorithm) set to 10⁻⁵. The relaxation parameter ω_r is set equal to 1.

FIGURE 6

Figure 6. Schematic representation of the 2D computational domain (A) and illustration of the used grid constellation (B) nested Cartesian grids.

Note that the plasma dynamics is resolved in the patch domain as the plasma density is supposed only defined on the patch grid. In this study, the size of patch is chosen large enough to capture the multiscale features of the plasma and the electromagnetic solution.

4.2.2. Results

Figure 7A exhibits the evolution of plasma density at different times obtained by the patch-based approach for a refinement H/h = 10 and the standard one when no patch is used (H/h = 1). The initial plasma starts to grow due to the overcritical field first uniformly. When its density is no longer negligible with respect to the cut-off density (which in our conditions is approximately 8.6 × 10²⁰ m⁻³ [32]), the plasma starts to interact with the electromagnetic field of the standing wave by behaving as a conductor. Due to polarization effects, the electric field is enhanced at the poles of the plasmoid in the field polarization direction leading to an increase of the ionization in those regions. The plasmoid thus elongates faster in the field direction forming a microwave streamer. The electromagnetic energy is no longer absorbed by the plasma as its density becomes overcritical, as it occurs in the streamer center. The stronger reflection of the waves causes a decrease of the rms field below to the critical field leading to decay of the plasma density in those regions. The results obtained by locally refining the region containing the plasma shows remarkably differences on the plasma density distribution compared to the case without patch. This is mainly due to a better resolution of the density gradients in the plasma front when increasing the refinement, which become extremely sharp when the streamer is formed (the differences are more important in Figure 7A at times 80 ns and 110 ns with respect to previous times where the plasma density gradients are smoother). This discrepancy in resolving the characteristic diffusion length of the plasma front leads to significant differences on the rms electric field distribution as well (see Figure 8). For the coarser grid an erroneous estimation of the plasma density is exhibited (see Figure 9A).

FIGURE 7

Figure 7. Contour plots showing the evolution of plasma density (A) in a microwave streamer from 20 to 110 ns without (top) and with patch H/h = 10 (bottom) in the patched region. The initial plasma density is a Gaussian with maximum 10¹⁵ m⁻³ in the center of the simulation domain and standard deviation 6 × 10⁻⁵ m. The maximum densities at the successive times are (the color scale is normalized at each time for the upper and lower plots), respectively 1.39 × 10¹⁷, 1.05 × 10²¹, 3.29 × 10²¹, and 4.38 × 10²¹ m⁻³. Contour plot showing the plasma density (B) at 110 ns in the reference case, without patch, having a uniform refinement of the whole computational domain, H_ref = H/10 (with H_ref being the step of the uniformly spaced grid), obtained in conditions of (A). The color scale is the same of the contour plot related to 110 ns in (A).

FIGURE 8

Figure 8. Contour plots showing the rms electric field distribution at 110 ns without (Left), with (Middle) patch H/h = 10 and for the reference case (Right) obtained in conditions of Figure 7 in the patched region. The color scale is between the minimum value of 2.38 MV/m and the maximum value of 5.05 MV/m (the scale is normalized for the three cases).

FIGURE 9

Figure 9. Profile of plasma density (A) and rms electric field (B) along the streamer axis, corresponding to conditions of Figure 7 at 110 ns, without (blue solid line) and with (black solid line) patch H/h = 10 and for the reference case (red solid line with circles).

Therefore, we recognize the efficiency of the local patch in improving the accuracy of the solution, both on the plasma density and on the electric field, as the refinement is increased. Note that the simulation results corresponding to the case with patch (h = H/10) are remarkably similar to results obtained in the case, without patch, on a mesh with the same fine refinement everywhere (and with the same CFL number), H_ref = H/10 (with H_ref being the step of the uniformly spaced grid), here refereed to as reference case (see Figure 7B for the plasma profile and Figure 8 for the rms electric field distribution both estimated at 110 ns). A more quantitative comparison is given in Figure 9 where the profiles of plasma density and rms electric field along the streamer axis at 110 ns are obtained in the case with patch and in the reference case. The results are in excellent agreement: the error, between the patch method and the standard FDTD method, concerning the plasma density and the streamer's elongation is less than 1%, whereas the rms electric field error is less than 2% (see Table 2).

TABLE 2

Table 2. Streamer's length (in units of λ), maximum plasma density (in m⁻³) and maximum rms electric field (in MV/m) at 110 ns for the case using the local patch and the reference case corresponding to conditions of Figure 7.

From a performance point of view, the theoretical gain expected when using the local patch with a refinement of H/h = 10 is about 15% compared to the case having a uniform refinement (of the same level) of the whole computational domain (see Table 3) according to the computational complexity analysis of section 3.2, as the average number of iterations to cope with the iterative patch algorithm is about k = 7.92 (here k corresponds to the average of all iterations needed for convergence over all wave periods of the simulation). Note that the configuration here analyzed is not the optimal one, as the patched region is sized to contain the whole developed microwave streamer. In so doing, the patch characteristic length becomes important with respect to the coarse grid step, i.e., L_p ≈ λ/3 >> λ/100 ≡ H, and we cannot expect important performance gains by considerably refining the patch (e.g., by one order of magnitude, as estimated in section 3.2) since the characteristic lengths ratio α dominates the refinement ratio β.

TABLE 3

Table 3. The performance gain between the patch method and the standard one (without patch) for the case of Figure 7 estimated using formula (28) with k being the average number of iterations needed for convergence by the patch algorithm over all wave periods (k_max and k_min are the maximum and the minimum number of iterations, respectively), α = 3.1740, β = 4, 10 and D = 2 (see section 3.2).

In our case, it is clear that the optimal configuration would be to only localize the patch on the plasma front where plasma density gradients are sharp and the electric field rapidly varies (e.g., in reference to condition of Figure 9, by locally refining regions around the streamer tips and leaving a coarser grid elsewhere). Practically, it would mean to add multiple patches having different refinement levels into the domain and adapt their size dynamically in order to follow the plasma evolution. This however is beyond the scope of this paper in which we introduce the patch-based approach and analyze its efficiency on more standard configurations.

5. Conclusion and Perspectives

We have proposed a new method for numerically solving electromagnetic wave propagation problems in plasma having a strong multiscale behavior such as the microwave discharge plasma modeling. The method is derived for a particular finite element formulation of the Maxwell's equations, which are coupled with the plasma transport equations via the time-varying electron current density, advanced in time with an explicit leapfrog time-stepping technique using a single constant time step. The resolution algorithm uses a patch of finite elements to calculate in an iterative manner the solution in the plasma region where a better precision is required. Numerical experiments have been performed demonstrating the usefulness of the patch technique on simplified configurations such as two-level nested Cartesian grid constellations. The application of a local well-suited patch considerably improves the solution precision while preserving the stability and the second-order accuracy of the Yee finite-element type scheme. The proposed method has been validated with 2D simulations of a microwave breakdown in air at atmospheric pressure.

However, a much better computational efficiency of the method can be provided by using patches adequately sized and refined around the propagation front of plasma. An adaptive version of the proposed algorithm is readily derivable as it is:

• Particularly flexible with respect to the conformity of the meshes,

• Extensible by definition to multiple patches having different level of refinements,

• Inherently adaptable to parallel computing due to the independence of problems defined on non-overlapping patches.

Further improvements on the time-stepping technique can be envisaged as well by employing a time integration using local time steps (e.g., by combining explicit and implicit time-stepping schemes). Moreover, a deep understanding of the convergence properties of the iterative algorithm for the electromagnetic wave propagation problem necessitates and they will be studied in future works. We conclude that both the inherent versatility and the particular efficiency of the proposed method can make it suitable to properly simulate the dynamics of microwave plasma discharges in fully three dimensions by providing more insights of their nature and behavior. This, of course, would allow the estimation of their applicability in the context of aerodynamic flow control.

Author Contributions

EA developed the numerical method and led the manuscript preparation. FR proposed the method and contributed to its development and to the article preparation. J-PB proposed the physical model describing the microwave discharge plasma. All three authors EA, FR and J-PB contributed to the analysis and discussion of the results.

Funding

This research is supported by the French Ministry of Defense (Direction Générale de l'Armement, DGA). EA benefits from an ONERA/DGA fellowship.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

1. Knight DD. Survey of aerodynamic drag reduction at high speed by energy deposition. J Propul Power. (2008) 24:1153–67. doi: 10.2514/1.24595

CrossRef Full Text | Google Scholar

2. Lashkov V, Mashek I, Anisimov Y, Ivanov V, Kolesnichenko Y, Ryvkin M. Gas dynamic ef fect of microwave discharge on supersonic cone-shaped bodies. In: 42nd AIAA Aerospace Sciences Meeting and Exhibit. Reno, NV (2004). p. 671.

Google Scholar

3. Fukunari M, Komurasaki K, Nakamura Y, Oda Y, Sakamoto K. Rocket Propulsion Powered Using a Gyrotron. J Energy Power Eng. (2017) 11:363–371. doi: 10.17265/1934-8975/2017.06.001

CrossRef Full Text | Google Scholar

4. Sakai O, Tachibana K. Plasmas as metamaterials: a review. Plasma Sour Sci Technol. (2012) 21:013001. doi: 10.1088/0963-0252/21/1/013001

CrossRef Full Text | Google Scholar

5. Hidaka Y, Choi E, Mastovsky I, Shapiro M, Sirigiri J, Temkin R. Observation of large arrays of plasma filaments in air breakdown by 1.5-MW 110-GHz gyrotron pulses. Phys Rev Lett. (2008) 100:035003. doi: 10.1103/PhysRevLett.100.035003

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Hidaka Y, Choi E, Mastovsky I, Shapiro M, Sirigiri J, Temkin R, et al. Plasma structures observed in gas breakdown using a 1.5 MW, 110 GHz pulsed gyrotron. Phys Plasmas. (2009) 16:055702. doi: 10.1063/1.3083218

CrossRef Full Text | Google Scholar

7. Cook AM, Hummelt JS, Shapiro MA, Temkin RJ. Millimeter wave scattering and diffraction in 110 GHz air breakdown plasma. Phys Plasmas. (2013) 20:043507. doi: 10.1063/1.4798424

CrossRef Full Text | Google Scholar

8. Schaub S, Hummelt J, Guss W, Shapiro M, Temkin R. Electron density and gas density measurements in a millimeter-wave discharge. Phys Plasmas. (2016) 23:083512. doi: 10.1063/1.4959171

CrossRef Full Text | Google Scholar

9. Boeuf JP, Chaudhury B, Zhu GQ. Theory and modeling of self-organization and propagation of filamentary plasma arrays in microwave breakdown at atmospheric pressure. Phys Rev Lett. (2010) 104:015002. doi: 10.1103/PhysRevLett.104.015002

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Chaudhury B, Boeuf JP, Zhu GQ. Pattern formation and propagation during microwave breakdown. Phys Plasmas. (2010) 17:123505. doi: 10.1063/1.3517177

CrossRef Full Text | Google Scholar

11. Arcese E, Rogier F, Boeuf JP. Plasma fluid modeling of microwave streamers: approximations and accuracy. Phys Plasmas. (2017) 24:113517. doi: 10.1063/1.5006651

CrossRef Full Text | Google Scholar

12. Kourtzanidis K, Rogier F, Boeuf JP. Gas heating effects on the formation and propagation of a microwave streamer in air. J Appl Phys. (2015) 118:103301. doi: 10.1063/1.4930163

CrossRef Full Text | Google Scholar

13. Kourtzanidis K, Boeuf J, Rogier F. Three dimensional simulations of pattern formation during high-pressure, freely localized microwave breakdown in air. Phys Plasmas. (2014) 21:123513. doi: 10.1063/1.4905071

CrossRef Full Text | Google Scholar

14. Kourtzanidis K, Rogier F, Boeuf JP. ADI-FDTD modeling of microwave plasma discharges in air towards fully three-dimensional simulations. Comput Phys Commun. (2015) 195:49–60. doi: 10.1016/j.cpc.2015.04.018

CrossRef Full Text | Google Scholar

15. Xiao K, Pommerenke DJ, Drewniak JL. A three-dimensional FDTD subgridding algorithm with separated temporal and spatial interfaces and related stability analysis. IEEE Trans Antennas Propag. (2007) 55:1981–90. doi: 10.1109/TAP.2007.900180

CrossRef Full Text | Google Scholar

16. Berenger JP. Extension of the FDTD Huygens subgridding algorithm to two dimensions. IEEE Trans Antennas Propag. (2009) 57:3860–7. doi: 10.1109/TAP.2009.2031906

CrossRef Full Text | Google Scholar

17. Ahmed I, Chen ZD. A hybrid ADI-FDTD subgridding scheme for efficient electromagnetic computation. Int J Numer Modell. (2004) 17:237–49. doi: 10.1002/jnm.543

CrossRef Full Text | Google Scholar

18. Li X, Sarris CD, Triverio P. Structure-preserving reduction of finite-difference time-domain equations with controllable stability beyond the CFL limit. IEEE Trans Microwave Theor Techn. (2014) 62:3228–38. doi: 10.1109/TMTT.2014.2366140

CrossRef Full Text | Google Scholar

19. Glowinski R, He J, Lozinski A, Rappaz J, Wagner J. Finite element approximation of multi-scale elliptic problems using patches of elements. Numerische Mathematik. (2005) 101:663–87. doi: 10.1007/s00211-005-0614-5

CrossRef Full Text | Google Scholar

20. Lozinski A, Pironneau O. Numerical zoom for advection diffusion problems with localized multiscales. Numer Methods Partial Differ Equ. (2011) 27:197–207. doi: 10.1002/num.20642

CrossRef Full Text | Google Scholar

21. Steger JL, Benek JA. On the use of composite grid schemes in computational aerodynamics. Comput Methods Appl Mech Eng. (1987) 64:301–20. doi: 10.1016/0045-7825(87)90045-4

CrossRef Full Text | Google Scholar

22. Dhia HB, Rateau G. The Arlequin method as a flexible engineering design tool. Int J Numer Methods Eng. (2005) 62:1442–62. doi: 10.1002/nme.1229

CrossRef Full Text | Google Scholar

23. Brunet A, Sarrailh P, Rogier F, Roussel JF, Payan D. Nonlinear patch method and application. In: European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2016). Crete Island (2016).

Google Scholar

24. Yee K. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans Antennas Propag. (1966) 14:302–7. doi: 10.1109/TAP.1966.1138693

CrossRef Full Text | Google Scholar

25. Nédélec JC. A new family of mixed finite elements in R3. Numerische Mathematik. (1986) 50:57–81. doi: 10.1007/BF01389668

CrossRef Full Text

26. Hesthaven JS, Warburton T. Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications. New York, NY: Springer-Verlag (2007).

Google Scholar

27. Yan S, Lin CP, Arslanbekov RR, Kolobov VI, Jin JM. A discontinuous galerkin time-domain method with dynamically adaptive Cartesian mesh for computational electromagnetics. IEEE Trans Antennas Prop. (2017) 65:3122–33. doi: 10.1109/TAP.2017.2689066

CrossRef Full Text | Google Scholar

28. MacDonald AD. Microwave Breakdown in Gases. 1st Edn. New York, NY: Wiley (1966).

Google Scholar

29. Zhu GQ, Boeuf JP, Chaudhury B. Ionization–diffusion plasma front propagation in a microwave field. Plasma Sour Sci Technol. (2011) 20:035007. doi: 10.1088/0963-0252/20/3/035007

CrossRef Full Text | Google Scholar

30. Hagelaar GJM, Pitchford LC. Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models. Plasma Sour Sci Technol. (2005) 14:722. doi: 10.1088/0963-0252/14/4/011

CrossRef Full Text | Google Scholar

31. Database. Available online at: www.lxcat.net (Accessed September 1, 2016)

32. Chaudhury B, Boeuf JP, Zhu GQ, Pascal O. Physics and modelling of microwave streamers at atmospheric pressure. J Appl Phys. (2011) 110:113306. doi: 10.1063/1.3665202

CrossRef Full Text | Google Scholar

33. Hagelaar G, Kroesen G. Speeding up fluid models for gas discharges by implicit treatment of the electron energy source term. J Comput Phys. (2000) 159:1–12. doi: 10.1006/jcph.2000.6445

CrossRef Full Text | Google Scholar

34. McCormick S. Fast adaptive composite grid (FAC) methods: Theory for the variational case. In: Defect Correction Methods. Vienna: Springer (1984). p. 115–21.

Google Scholar

35. McCormick S, Thomas J. The fast adaptive composite grid (FAC) method for elliptic equations. Math Comput. (1986) 46:439–56. doi: 10.1090/S0025-5718-1986-0829618-X

CrossRef Full Text | Google Scholar

36. McCormick S, Ruge J. Unigrid for multigrid simulation. Math Comput. (1983) 41:43–62. doi: 10.1090/S0025-5718-1983-0701623-0

CrossRef Full Text | Google Scholar

37. Wagner J. Finite Element Methods With Patches and Applications. Lausanne: EPFL (2006).

Google Scholar

38. He J, Lozinski A, Rappaz J. Accelerating the method of finite element patches using approximately harmonic functions. Comptes Rendus Mathematique. (2007) 345:107–12. doi: 10.1016/j.crma.2007.06.006

CrossRef Full Text | Google Scholar

39. Cohen GC. Higher-Order Numerical Methods for Transient Wave Equations. Heidelberg: Springer (2001).

Google Scholar

40. Monk P. An analysis of Nedelec's method for the spatial discretization of Maxwell's equations. J Comput Appl Math. (1993) 47:101–21. doi: 10.1016/0377-0427(93)90093-Q

CrossRef Full Text | Google Scholar

41. Vikharev A, Eremin B. Microwave discharge in a quasioptical resonator. Zh Éksp Teor Fiz. (1975) 68:452–5.

Google Scholar

42. Barashenkov V, Grachev L, Esakov I, Kostenko B, Khodataev K, Yur ev M. Breakdown in air in a rising microwave field. Tech Phys. (2000) 45:1265–1270. doi: 10.1134/1.1318961

CrossRef Full Text | Google Scholar