- Email: [email protected]
Multi-scale gas discharge simulations using asynchronous adaptive mesh refinement
Computer Physics Communications 181 (2010) 247–258
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Multi-scale gas discharge simulations using asynchronous adaptive mesh refinement Thomas Unfer a,b,c,∗ , Jean-Pierre Boeuf a,b , François Rogier c , Frédéric Thivet d a
Université de Toulouse, UPS, INPT, LAPLACE (Laboratoire Plasma et Conversion d’Energie), 118 route de Narbonne, F-31062 Toulouse cedex 9, France CNRS, LAPLACE, F-31062 Toulouse, France c ONERA, 2 avenue Edouard Belin, 31400 Toulouse, France d Université de Toulouse, ISAE, 10 avenue Edouard Belin, 31055 Toulouse cedex 4, France b
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 2 July 2008 Received in revised form 22 September 2009 Accepted 22 September 2009 Available online 3 October 2009 Keywords: Multi-scale Asynchronous Transport problem Local time stepping Gas discharges CFL condition Adaptive mesh refinement
The breakdown of a gas gap at high products of pd (pressure × distance) is a multi-scale phenomenon in both time and space. This is especially true when the plasma is interacting with a gas flow, a problem of considerable recent interest in the context of aerodynamic applications of surface discharges. This paper presents a contribution to the numerical modeling of such discharges. We describe here a new approach for adaptive meshing which is suitable for use with the explicit asynchronous integration scheme described in our previous publication. Rather than relying on a family of nested grids as is commonly done, this technique is based on a single unstructured mesh with possible non-conforming cells at the interface between coarse and fine areas. Substantial computational time saving has been achieved for a surface dielectric barrier discharge configuration of the kind proposed as plasma actuators for flow control. © 2009 Elsevier B.V. All rights reserved.
1. Introduction In the recent years, substantial research efforts have been carried out throughout the world aiming at controlling aerodynamic air flow with plasma actuators [1,2]. Surface dielectric barrier discharges (DBDs) have been proposed as plasma actuators able to modify the boundary layer of a flow around an airfoil. The numerical simulation of these plasma actuators is difficult because of the space and time scales which must be taken into account span many orders of magnitude. One useful approach is adaptive meshing, and Adaptive Mesh Refinement (AMR) techniques have already been successfully applied in the field of gas discharges to investigate the physics of the propagation of fast ionization waves known as streamers [3]. In a previous paper [4], we derived an asynchronous integration procedure for the charged particle transport equations on which fluid models of gas discharge phenomena are based. This technique reduces both numerical diffusion and CPU time and is particularly adapted for meshes with strong variation in cell size. When Berger and Oliger first developed AMR for block-structured meshes in the
*
Corresponding author at: Université de Toulouse, UPS, INPT, LAPLACE (Laboratoire Plasma et Conversion d’Energie), 118 route de Narbonne, F-31062 Toulouse cedex 9, France. E-mail address: [email protected] (T. Unfer). 0010-4655/$ – see front matter doi:10.1016/j.cpc.2009.09.017
©
2009 Elsevier B.V. All rights reserved.
mid 1980’s [5], they introduced a family of nested grids with time synchronization between large and fine grids based on fractional time steps. The time steps on the finer grids were chosen so as to keep the CFL number constant with respect to the larger grids. In this paper the concept of a family of nested grids is not used. In the proposed Asynchronous Adaptive Mesh Refinement (AAMR) method the family of nested grids can be replaced to advantage by a structure introduced by Burgarelli et al. in [6], namely an Autonomous Leaves Graph (ALG). This paper is organized as follows: Section 2 briefly explains the motivation of this work. The model for a plasma discharge in air is presented in Section 3, and the principles of the asynchronous time integration are summarized in Section 4. Section 5 introduces the specific data structures needed for AAMR, and Section 6 explains how the asynchronous scheme is to be used with an unstructured mesh with non-conforming cells. The AAMR algorithm is then introduced in Section 7, and some simulation examples from numerical calculations and comparisons of discharges in the typical surface DBD configuration used for aerodynamic flow control are presented in Section 8. 2. Motivation of this work There is a growing interest in numerical simulations of the generation and evolution of atmospheric pressure air plasma and,
248
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
especially, in the interaction of the plasma with a gas flow. These phenomena are very demanding in terms of CPU time because of the large variation in time and length scales that are involved in the physics of the plasma and plasma/flow interaction. For example, the time scale associated with the coupling between the charged particle transport equations and Poisson’s equation is known as the dielectric relaxation time τd , and a condition for numerical stability in explicit integration schemes is that τd not be exceeded. In atmospheric pressure discharges τd can be as small as 10−12 s. The CFL constraint associated with the electron transport depends on the grid size but, in practice, can be lower than 10−12 s. One concept of actuator for aerodynamic flow control is an asymmetrical surface Dielectric Barrier Discharge (DBD) with a sinusoidal voltage waveform in the kHz frequency range. Momentum transfer from the ions to the neutral gas molecules induces a net gas flow on a time scale 10−4 to 10−3 s. Consequently the time scales involved in this problem extend over more than ten orders of magnitude. The problem is also covering space scales of various orders of magnitude since the Debye sheath in atmospheric plasmas can be as small as a few micrometers whereas the plasma extension of a filamentary discharge above a surface can be several centimeters generating an airflow over a few tens of centimeters. As a result the physics of the discharge is truly multi-scale and it is therefore essential for modeling to develop dedicated numerical tools that are able to tackle these issues. We show in this paper that associating the asynchronous time integration method described previously in [4] with adaptive mesh refinement techniques can lead to excellent synergy and considerable computation time saving. 3. A gas discharge model for air A minimal air plasma model is considered, consisting of only three species, electrons (e), one type of positive ion (p) and one type of negative ion (n). The chemistry accounts for only four reactions: direct ionization by electron impact, electron attachment to form a negative ion, and electron/positive ion and negative ion/positive ion recombination. A detailed description of this model can be found in [7]. While this plasma chemistry is admittedly quite simple compared to a real situation, we expect that the simple model will capture the essential physics of the plasma/flow coupling. The discharge model is based on the following fluid equations:
• Continuity equations
∂ ne Γe = (α − η)Γe − re− p ne n p , + ∇. ∂t ∂np Γp = α Γe − re− p ne n p − rn− p nn n p , + ∇. ∂t ∂ nn Γn = ηΓe − rn− p nn n p , + ∇. ∂t
(1) (2) (3)
ne , n p , nn are respectively the electron, positive ion and negative ion densities. Γe , Γp , Γn are the corresponding particle fluxes. α is the ionization coefficient (first Townsend coefficient), η is the attachment coefficient. Volume recombination is described using re− p and rn− p which are the electron/positive ion and negative ion/positive ion recombination coefficients. • Momentum equation in drift-diffusion form
kB Te ∇ ne , Γe = μe ( E ) −ne E − e kB T p np , ∇ Γp = μ p ( E ) n p E − e
(4) (5)
Fig. 1. Simulation domain for surface dielectric barrier discharge for flow control.
kB Tn nn , ∇ Γn = μn ( E ) −nn E − e
(6)
μe , μ p , μn are respectively the electron, positive ion and negative ion mobilities and are function of the local electric field E only. The diffusion coefficients are expressed as a function of the various mobilities and temperatures using the Einstein relation. In practice the electron temperature T e is taken about 1 eV and the ion temperatures T p , T n are taken equal to the background gas temperature. • Poisson’s equation ε E ) = e (n p − ne − nn ), ∇.(
(7)
ε is the dielectric permittivity of the media (air or dielectric layer), e is the elementary charge. In collisional discharges, and this is particularly true at atmospheric pressure, the electron energy transport is not very important. In fact the energy gained by the charged particles in the electric field tends to be locally balanced by the losses due to collisions with the neutral gas. Under these conditions, the electron energy equation can be replaced by the so-called Local Field Approximation (LFA). The loss of precision is minimal for conditions of interest and considerable computational time can be saved. With the LFA, the transport and rate coefficients (mobilities, reaction coefficients) depend the local value of electric field (or, more precisely, on the reduced electric field, E / N, the ratio of the electric field strength to the neutral density), and their values are obtained by solving the steady state and homogeneous (no spatial gradients) Boltzmann equation. Transport and rate coefficients are calculated using the freeware BOLSIG+ [8,9], and the data tables in tabular form are input to the code. We consider the same numerical experiment as in our previous paper [4]. The geometry is reproduced in Fig. 1. The charging of the dielectric layer by the charged species and the induced potential due to the surface charge is taken into account. 4. The principles of the asynchronous scheme 4.1. Asynchronous integration The asynchronous scheme is a finite volume scheme for unsteady calculations, which is particularly suitable for multi-scale transport problem. The details of the algorithm have been published in [4]. For convenience the principles of the asynchronous time integration are briefly summarized in this section for a onedimensional situation. The extension to multiple dimensions is obvious:
∂ n ∂Γ + = S. ∂t ∂x
(8)
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
249
The asynchronous scheme assumes that each flux or source term is updated independently according to a refresh time tag with respect to a global simulation clock. ni are the density values at the cell centers xi = i x, v i + 1 are the velocity values taken at the cell 2
interfaces xi + 1 = (i + 12 )x. For the scheme to be stable, each local 2 time step is limited by the local value of the CFL condition. Each conservative variable has its own value, value time tag and variation rate. The variation rate represents the instantaneous balance of fluxes and source term of the conservative variable. The idea is to let each component of this balance to evolve independently in time:
Γi + 1 Γi − 1 ∂n S ∂n Γ + ∂n Γ − ∂n 2 2 = Si − + = + + . ∂t i x x ∂t i ∂t i ∂t i Different time variables are needed: t simulation : the current discrete simulation time, t ni : time tag associated with the stored value of density ni , t iΓ : refresh time tag of flux Γi , t iS : refresh time tag of source term S i . During the simulation, t simulation jumps discretely from the most urgent refresh time tag to the next most urgent. Initialization 1. Initialize all density and fluxes. 2. Initialize all refresh time tags to the initial time. Proceed until the simulation time is completed 1. Find the most urgent flux or source term to be refreshed. 2. t simulation becomes this most urgent refresh time tag. 3. Compute the value of the density on both side of the interface needed for flux computation (only a few cells are involved) or of the cell for source term computation: • n j = n j + ∂∂nt (t simulation − t nj ),
• t nj = t simulation with j ∈ [i , i + 1] for first order flux, j ∈ [i − 1, i + 2] for second order flux (e.g. MUSCL) or j = i for source term.
4. Compute the new value of the flux or source term and update its variation rate according to:
•
∂n Γ + ∂t i ∂n Γ − ∂ t i +1
=−
Γi + 1
2 x , Γi + 1
S
or ∂∂nt i = S i . 5. Compute the new refresh time tag of this flux or source term according to the local CFL condition or reaction kinetic t iΓ = t iΓ + t CFL when updating a flux or t iS = t iS + t S when updating a source. 6. If the solution is to be monitored before the next most urgent refresh time tag, build the solution at the output time tag and store it.
•
output ni
=
2
x
∂ n S ∂ n Γ + ∂ n Γ − output = ni + + + t − t ni ∂t i ∂t i ∂t i
for all cells. 4.2. Coupling with Poisson’s equation Poisson’s equation is solved at each cell center with time steps corresponding to the global dielectric relaxation time (minimum of all the local dielectric relaxation times). The electric field is defined on each interface between two cells and is advanced in time between updates of the solution of Poisson’s equation using the following estimator with j c , the conduction current going through
Fig. 2. Flux data type.
the interface:
∂ E i , j jt −t − jtc (t ) = c . ∂t ε0
(9)
Some justification of this coupling can be found in [4]. 5. Specific data structures for AAMR The basis of gas discharge modeling is to integrate in time the transport equations for the charged particles coupled with Poisson’s equation for the electric field. This essentially implies solving repetitively an elliptic equation (Poisson’s equation). For purposes of accuracy when solving Poisson’s equation, it is important that neighboring cells do not differ in size by more than a factor two. This requirement can be fulfilled if, for instance, the mesh has a quadtree structure as defined in [10]. The idea here is to create a simplified ALG (Autonomous Leaves Graph) structure, mimicking a quadtree. Because the full features of ALG are not needed to generate a quadtree-like mesh, a simpler structure can be defined which is expected to speed somewhat the computations. The resulting data structure is suitable for both asynchronous time integration of transport equations and local refinement/merging operations over an unstructured mesh. Moreover, because the asynchronous scheme keeps jumping from one memory area to another while updating interface fluxes, the memory cache misses are much more likelier to occur for asynchronous time integration than for a standard integration scheme. In order to mitigate this effect the concept of collocated data (i.e. keeping all the data corresponding to one single cell in the same memory location) is essential for asynchronous time integration. The useful data can be split in two parts: cell-centered data (densities, for instance) and interface data (flux). Consequently two data types have been defined, one gathering all interface data, termed “Flux data type”, and the other gathering the cell-centered data, or “Cell data type”. 5.1. Flux data type The Flux data type consists of all the data required for computation at the interface between two cells for each transport equation (e.g., fluxes and electric field). It points towards two cells of the Cell data type, one being the upstream cell called Cell_out in Fig. 2 and the other being the downstream cell called Cell_in in the same Fig. 2. The flux of a conservative variable is counted positively if it goes from Cell_out to Cell_in. 5.2. Cell data type The Cell data type consists of all the cell-centered data for each equation (value, variation rate, time of value, etc.) plus additional parameters (boundary, Poisson matrix coefficients, etc.). It also contains additional information about the geometry, position and size of the cell, level of refinement and coordinates in a virtual uniform mesh. In two dimension the Cell data type points towards at least 4 Flux objects (Flux_East_A, Flux_South_A, Flux_West_A, Flux_North_A) when the cell is connected to other cells of the same size or greater. Up to 4 additional Flux objects pointers can be used (Flux_East_B, Flux_South_B, Flux_West_B, Flux_North_B) when the cell is connected to smaller cells (see Fig. 3). The Cell data type also contains additional pointers to the corresponding actions which are used for the asynchronous integra-
250
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
Fig. 3. Cell data type in 2D. Fig. 5. Schematic view of a local refinement.
Fig. 6. Wrong AAMR scenario.
Fig. 4. Schematic view of a regular mesh.
tions (see Section 7.1, these pointers are used to recall an action after an AAMR event so that the fluxes can be refreshed immediately). The Cell data type is additionally connected to the next and previous cell so as to have a chained list of the cells within the mesh (see Section 7.2). 5.3. The mesh: A composite object structure The mesh is built using fluxes to connect cells and is an Autonomous Leaves Graph. By construction, the fluxes are oriented along the positive directions of the mesh. A regular mesh only uses the “_A” type pointers to flux as shown in Fig. 4, whereas the “_B” type pointers to flux are used to connect cells of different sizes around some refined zone (see Fig. 5). Boundary cells can be connected to a single ghost cell. This object-like construction of the mesh satisfies the two following properties:
• Any block refined mesh that can be generated using classical AMR can be reproduced using a single unstructured mesh using this construction. • Any further refinement or merging can be done while modifying only the storage of the local elements in memory and has no effect on the rest of the mesh. 5.4. Adaptive meshing rules Two principles are postulated to develop the algorithm. These are, at each time, the mesh must preserve proper cells connectivity and cells may differ in size from their neighbors by no more than a factor two. To satisfy these principles a few simple rules can be formulated to avoid any “illegal” refinement or merging to occur.
• AAMR Rule #1: In N dimension a cell can be refined in 2 N if its size is larger than or equal to the size of any of its direct neighbors. (Direct neighbor means that the cells are connected by a flux, N is the dimension of the simulation.)
• AAMR Rule #2: 2 N cells of identical size can be merged if their size is smaller than or equal to the size of any of their direct neighbors. Note. These are the minimal rules. However, in practice and for precision purposes, especially when coupling drift-diffusion equations with Poisson’s equation which is our present concern, “direct neighbors” can be extended to “neighbors” which means taking into account the 2 N diagonal neighbors (North–East, South–East, South–West, North–West in 2D). This prevent cells from differing in size of a factor 4 with their diagonal neighbors. This feature is required for accurately solving Poisson’s equation. These two rules are sufficient to insure that the mesh integrity is preserved and that the neighboring size principles are not violated. However it is not sufficient to insure that the adaptive meshing process is reversible, the AAMR process may get trapped in some configuration. Fig. 6 shows an example in 1D which explains why things can go wrong without an additional rule. In fact the third AAMR event in the figure, i.e. “the merging of cell 4 & 5” is problematic because afterward cell 3 or cell 6 can no longer be merged. This behavior shall be forbidden. Each cell shall “remember” a virtual parent cell of greater size in case of merging so as to make the AAMR process reversible. This can be done by assigning coordinates to each cell corresponding to its coordinates in an hypothetical regular mesh of the same size over the whole domain and following Rule #3:
• AAMR Rule #3: Only a cell with odd coordinates is available for merging with cells with the successive even coordinates. 6. Using the asynchronous scheme over a non-conforming unstructured grid 6.1. Mass conservation Consider now the asynchronous scheme in a finite volume formulation with no source term and zero flux at the boundaries. The
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
251
Fig. 8. Time/space diagram of the asynchronous scheme assuming positive drift velocity.
• First case: t p tq t p +1 tq+1 t q +1
ni
t
= ni q −
i =1
N Cell N faces(i ) ∂n Γ t ni + t ∂ t i→ j j =1
i =1
N Cell
nti
=
(10)
i =1
because by definition of the algorithm:
t q +1
ni
Γi → j ∂n Γ ∂n Γ =− =− . ∂ t i→ j x ∂ t j →i
(11)
N Cell is the total number of cells, N faces(i ) is the total number of faces of the cell i. Fig. 7 shows the electron density relative variation in the numerical experiment described in Section 8 with the following hypotheses: both the electron source term and the electron fluxes at the boundaries are taken equal to zero. This test basically evaluates the numerical performance in term of mass conservation under drift diffusion in an electric field with a non-conforming mesh. Density is conserved up to round-off error precision.
The demonstration is given in one dimension for the advection equation for the sake of simplicity. With the notations of Fig. 8, suppose that the drift velocity is smooth and the local time steps are proportional to the local CFL condition, i.e.:
kCFL x v i− 1
kCFL
x
kCFL
x
2
t p =
kCFL x v i+ 1 2
= tq
1−
v
v
∂ v x . 1+ ∂x v
∂ v x kCFL ∈ ]0, 1[, ∂ x 2v
∂ v x 1+ ∂ x 2v
t = ni q v i
(12)
+ ni
t
− ni q ∂(nv ) ∂ n (t p − tq ) + = −v i τ. tq ∂x ∂t tq x
(14)
Hence there is consistency since t p − tq = O (tq x) and
τ = O (tq ).
• Second case: t p tq t p +1 t p +2 tq+1 t q +1
ni
t
1
t
t
t
p p +1 tq ni q v i + 1 − τ ni − v − t p ni − v 1 i − 12 1 i − 12 2 t p +2 , − (tq − τ − t p )ni − v 1 i− 1
= ni q −
x
2
∂ v x ∂n ∂n − v i x + v i (τ + t p ), ∂x 2 ∂x ∂t
v i ∂n ∂(nv ) tq (t p − tq )τ = ni − tq + ∂x tq x ∂ t + t p τ + t p (t p − tq ) . (15)
t p +2 ni v i− 1 2
6.2. Consistency under local CFL condition
tq =
∂ v x , ∂x 2 ∂n ∂n tp tq tq ∂ v x ni v i − 1 = ni v i − ni − v i x − v i (t p − τ ), 2 ∂x 2 ∂x ∂t ∂n ∂n t p +1 tq tq ∂ v x ni v i − 1 = ni v i − ni − v i x + v i τ , 2 ∂x 2 ∂x ∂t
∂ n (t p − tq ) ∂(nv ) t q +1 tq ni = ni − tq + vi τ , ∂x ∂t tq x tq ni v i + 1 2
N Cell
i =1
t
2
scheme is clearly conservative:
nti +t =
t
p tq ni q v i + 1 − τ ni − v 1 i − 12 2 t p +1 − ni − v (tq − τ ) , 1 i− 1
Fig. 7. Mass conservation test.
N Cell
1
x
t q +1
ni
t = ni q v i
t − ni q
In this case, necessary gives consistency.
τ tq − t p so τ = O (tq x), which
• Third case: t p tq tq+1 t p +1
t q +1
ni (13)
For x v ∂∂ vx the following condition is imposed on the local time steps: |t p − tq | tq , this means that three cases have to be considered:
t q +1
ni
t tp tq ni q v i + 1 − tq ni − , v 1 i − 12 2
∂ n t p − τ ∂(nv ) t . = ni q − tq + vi ∂x ∂ t x t
= ni q −
1
x
(16)
In this case, necessary tq τ t p so t p − τ = O (tq x), which gives consistency.
252
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
6.4. CFL condition on a conforming and a non-conforming cell
Table 1 Unstable counterexample in 2D with a 20% desynchronization. Time
0
Density Rate Flux #1 Flux #2
0.4
0.5
0.8
−2
0.2 −1.2
−0.004
1 1
1 0.2
0.08 0.28 0.08 0.2
1
With the notation of Fig. 3, let us consider, for instance that Flux_East_A is positive. Consequently its refresh rate is function of all the fluxes emptying the cell and is equal to the cell time step:
kCFL
t Flux = t Cell = Faces
As a result the first-order asynchronous scheme is convergent. The local CFL condition hypothesis is essential for the convergence of the scheme. 6.3. An alternative formulation of the asynchronous scheme We showed in our previous paper [4] that each local contributor (fluxes and source terms) to the variation rate of the density can be updated asynchronously, and this principle has been used successfully for 2D calculations. However, it is possible to build a small counterexample showing that the stability of the scheme is not assured for arbitrary values of the initial desynchronization between the fluxes. If one considers a unit-sized cell with 20% desynchronization between two fluxes emptying the cell. The local time steps last both 0.5 and the desynchronization is 0.1. Table 1 shows the computation. Furthermore depending on the problem, this highly segregated treatment may not be the optimum in term of computational time or memory load. If any flux or source is independent, it means managing a corresponding action in the scheduler. The algorithm can be optimized using the fact that the CFL condition on a single flux is dependent on the fluxes emptying the same cell in the other directions. That is why the CFL condition can be seen as a condition on the upwind cell. Actually a flux shall be refreshed more often than it takes for the upwind cell to be emptied. Consequently it makes sense to synchronize all the fluxes which are emptying a single cell. For gas discharges simulations with a formulation similar to the model defined in Section 3, all the source terms can also be synchronized with the electron fluxes because the formulation of the chosen chemistry is expressed with respect to the electron flux using Townsend coefficients (see [7]). As a result, a single asynchronous action is required for each specie in each cell. Nevertheless if at a certain point of the simulation, a numerical flux changes direction, it has to be synchronized with the cell that is now being emptied by this flux. The asynchronous action of this cell shall be “recalled”, i.e. this cell shall be refreshed immediately. The alternative formulation makes the scheme faster because less actions are used in the Discrete Time Scheduler leading to less memory usage. Furthermore the stability condition in multiple dimensions can be proved in the case of the first order upwind asynchronous scheme if the leaving fluxes are synchronized. Supposing v x > 0 and v y > 0, the stability condition is equal to the local CFL condition with the notation of Fig. 8:
t q +1 ni , j
t = ni ,q j
1
−
x
t
t
tq ni ,q j v ix+ 1 + tq ni ,q j v 2
t q +1
n ( z )i j v x
−
i − 12
tq t q +1
ni
t
ni q
t q +1
dz −
n ( z )i j v
y i + 12
y
kCFL is the stability margin (kCFL ∈ ]0, 0.5[ for MUSCL scheme). The fluxes contribution are prescribed as follows, f CFL Flux = 0 if the actual numerical flux is directed towards the cell and f CFL Flux = V CFL when the flux is directed outwards. The length of the cell Cell Length is taken equal to the cell size if there is only one flux per cell face or half the size if there are two. Additional stability can be introduced by defining an equivalent CFL condition for the source term with f CFL Source Term being the stability constraint of the unstable part of the source term (negative part). This definition is equivalent to the classical multiple dimensions time step for a synchronous explicit scheme. For instance with v x > 0 and v y > 0, vx vy , f CFL Flux y − = 0, f CFL Flux x+ = and f CFL Flux x− = 0, f CFL Flux x+ = x y x y
the time step is t Cell = kCFL v x y + v y x . If during the simulation the sign of the flux changes then it recombines its refresh rate and CFL contribution with the other cell and this cell is refreshed immediately. For gas discharge modeling the equivalent CFL velocity used is:
v CFL = μ( E ) E +
D
xC 1 − xC 2
+ x
n i μi
ε0
.
(19)
D is the diffusion coefficient, the two factor in front of the x2 ) can be dropped because classical diffusion time step (t = 2D the influence of the diffusion is taken into account for each fluxes which are actually emptying the cell. So for a purely diffusive problem, the local time step used are twice larger except in the cell holding the maximum density. xC 1 − xC 2 is the distance between the two cell centers, x is the cell size. When an AAMR event occurs, the new cells and their direct neighbors have to be refreshed immediately to insure the stability of the scheme, because the scheduled actions for the neighbors were planned assuming the previous configuration of the mesh. 6.5. Numerical diffusion When writing second-order Taylor development of the firstorder asynchronous scheme of the advection equation and isolating the numerical diffusion coefficient, one obtains in the first case:
(t p − τ )2 τ + τ 2 (tq − τ ) 2 2tq x (t p − tq )τ x . + + tq 2v
D async = v 2 −
tq
−v
(20)
The first and last terms also appear when calculating the numerical diffusion coefficient for the synchronous scheme giving x (1 − vxt ). In contrast for the asynchronous rise to D sync = v 2 scheme those two terms cancel each other at least at first order in x when kCFL = 1. The third term is also second order in x, only the second term as to be evaluated knowing that to first order tq t p and (tq − τ )τ 4 q with 0 τ tq . Then the asynchronous diffusion coefficient is lower than:
dz ,
D async
tq
tq x tq y 1− v 1 − v 1 . x i + 2 x i + 2
(18)
,
t 2
i − 12
f Flux + f Source
(17)
v x 8
.
(21)
This means that if the velocity variation over the mesh is greater than 25%, the asynchronous scheme becomes less diffusive
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
253
Fig. 9. Enhanced Discrete Time Scheduler. Fig. 10. Poisson’s equation discretization on non-conforming cell.
than its synchronous equivalent. From a physical point of view, the numerical diffusion in the synchronous scheme can be seen as a consequence of the fact that in the area where the CFL number is lower than one, information travels faster than it should. The asynchronous scheme cancels this effect to first order because information cannot go faster than the fastest local wave corresponding to the local CFL. However the asynchronous scheme introduces diffusion due to the local desynchronization τ (second term in Eq. (20)). Notes.
• Case one and three are corner cases. Nevertheless, to first order, they lead to no numerical diffusion coefficient of zero.
• The estimation of the asynchronous numerical diffusion coeft ficient correspond to the worst case τ = 2 q . Actually τ keeps varying in time so the average value of numerical diffusion should be lower.
6.6. Extended limiter function For the non-conforming cell, the definition of the slope used in the limiter function has to be extended. See the density splitting function in Section 7.3.2 for the details. 7. The asynchronous adaptive mesh refinement algorithm 7.1. Extension of the asynchronous scheme The asynchronous scheme is compatible with adaptive meshing because it is based on interface synchronization. However, the stability condition must be insured for the new cells and their neighboring cells after an AAMR event. Because the neighboring cell refreshment was planned according to the previous mesh configuration, the scheduling has become obsolete and needs to be updated after an AAMR event. The corresponding actions have to be recalled so that the cell can be refreshed immediately. Hence, the Discrete Time Scheduler has to manage the additional functions of task recalling, task creation, task suppression. This is accomplished by adding new pointers into the structure to retrieve the pending actions from the mesh and the task from the pending action. (See Fig. 9.) 7.2. Poisson’s equation Poisson’s equation has to be solved on the unstructured mesh. Using the notation of Fig. 10, the equation is discretized classically for each cell using a finite volume formulation of the electric field = e (n p − ne − nn )). The electric field flux conservation ( εE .dl
values at the interface have to be expressed in function of the potential values in the surrounding cells. For example on the upper east face of the cell O between point A and point B the potential derivative across the vertical face is obtained using:
∂V V East − V O VA − VB =α +β ∂x dEast− O AB with V A =
1 (V O 4
+ V North B + V North East + V East ), V B = V East + V East B ), α = 43 , β = − 13 .
(22) 1 (V O 3
+
When summing up the contributions of each face, one gets the full Poisson matrix. In a homogeneous area the scheme is equivalent to the standard five points scheme, whereas a thirteen points scheme is used for a cell completely surrounded by smaller neighbors. Unfortunately, Poisson’s matrix loses its convenient symmetry, so a conjugate gradient solver cannot be used. This non-symmetry can be seen, for instance, if one looks at the discretization of Poisson’s equation in cell O . The discretized equation in cell O uses the potential value in the North East cell. However, because the discretization in the North East cell uses the standard five points scheme, it does not use the value of the potential in the cell O . So, the matrix is clearly non-symmetrical. We decided to solve Poisson’s equation with the restarted GMRES algorithm [11]. From a practical point of view, the Poisson matrix is stored twice in memory — once in the object-mesh as part of the cell data type and another time in a 13 × X array where X is the cell number with a correspondence look-up table. Any AAMR event will lead to a modification of the Poisson matrix within the objectmesh. As a result the array and look-up table have to be updated before solving the equation. The GMRES algorithm is used to iterate only within the array to speed up calculations because looping within the object-mesh is rather slow. 7.3. The asynchronous adaptive mesh refinement procedure 7.3.1. General algorithm With respect to the asynchronous algorithm (see Fig. 11), the AAMR algorithm consists of adding a testing subroutine which may trigger an AAMR event (see Fig. 12). There are many ways to implement this testing, for instance, some specific asynchronous task dedicated to the testing could be added into the Enhanced Discrete Time Scheduler (see Section 7.1). It makes sense to associate the AAMR testing frequency to the refresh rate of the electron fluxes. The detailed pseudo-code of the “Test for AAMR” function follows: 1. if the source term has been refreshed more than n (user defined) times since last AAMR testing of the current cell: (a) if local refinement criterion is satisfied for the current cell i. Refine the current cell
254
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
Fig. 11. Asynchronous algorithm overview.
Fig. 13. How to split cells in 2D.
(i) update the cell chained list (j) solve Poisson’s equation 2. else (a) for all the neighboring cells within the mesh (Ghost cell excluded) i. if a neighboring cell is greater than the current cell, refine this cell (b) try again to refine the current cell
Fig. 12. AAMR algorithm overview.
(b) else if the current cell is available for merging (AAMR Rule #3) i. if the local merging criterion is satisfied for the current cell A. Merge current cell with its neighbors ii. else refresh the cell (c) else refresh the cell 2. else refresh the cell 7.3.2. The cell refinement procedure The cell refinement procedure consists of replacing a cell with smaller cells when the refinement criterion is met. However, it may happen that a cell cannot by refined because one of its neighbors is of greater size. In such case, the AAMR algorithm shall start by refining the neighboring cell before even if this cell does not need to be refined. The cell refinement procedure is therefore recursive: 1. If the current cell can be refined (a) update the density of the neighboring cells (b) disconnect the current cell from the fluxes coming from its neighbors (c) suppress the cell action task (d) replace the cell with 2 N new cells and fill them with the density splitting function, fill the Poisson matrix, etc. (e) create the new inner fluxes and possible external fluxes (of the _B type for the neighbors) (f) re-connect the fluxes (g) create the new cell action tasks (h) recall the surrounding cells
Note. Step (j) is obviously specific to gas discharge modeling. When Poisson’s equation is solved following an AAMR event, the next update of the electric field can be postponed by one dielectric relaxation time in the asynchronous loop. In case of recursive calls of the refinement procedure, it is only necessary in principle to solve Poisson’s equation once. Density splitting function. This function shall split conservatively and positively the density from the coarse cell to the finer cells. Here is an example in 2D with the notation of Fig. 13. First define the slopes:
sx = minmod
n E − n O n O − nW
s y = minmod
,
dE O
dO W
nN − n O n O − n S
sx B = minmod
s y B = minmod
,
dN O
,
(23)
,
(24)
dO S
n E B − n O n O − nW B
,
dE B O
dO W B
nN B − n O n O − n S B
,
dN B O
dO S B
.
,
(25)
(26)
Then split the density:
n1 = n − s x n2 = n + s x
x 2
x
n3 = n − sxB n4 = n + sxB
2
+ sy
+ syB
x 2
x 2
y
− sy
y 2
y
− syB
(27)
,
2
2
,
(28)
,
(29)
y 2
.
(30)
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
This splitting is conservative because n1 2x
n3
x y 2
2
+ n4
x y 2
2
y 2
+ n2 2x 2y +
= nx y. The positivity in insured by the use
of the limited slopes. The limiter function for MUSCL scheme used in flux computation is extended to the non-conforming cells using the same
Fig. 14. Discharge current vs simulated time.
255
slopes. If a coarse cell is surrounded by finer cells as in Fig. 13, then the limiter uses the slope s x for the upper fluxes (E or W) and sxB for the lower fluxes (E_B or W_B). If for instance there is one single cell on the east side then the average slope sx +2sxB is used to estimate the density on the right side.
Fig. 15. CPU time and cell number ratio vs simulated time.
Fig. 16. AAMR mesh evolution with time, ith level represents the size of the local square cells with a size of
L 2(i −1)
.
256
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
Fig. 17. Comparison between AAMR (left) and asynchronous scheme (right) at t = 45 ns from top to bottom: electron density (log(m−3 )), positive ion density (log(m−3 )), negative ion density (log(m−3 )) and equipotential lines (V).
7.3.3. The cells merging procedure 1. update the cells to be merged (this may already be done while evaluating the merging criterion) 2. disconnect the cells to be merged from the fluxes connected to their neighbors 3. suppress the useless actions 4. suppress the useless fluxes and cells 5. fill the merged cell with the averaged values of the ex-finer cells (density, etc.) 6. reconnect the external fluxes 7. recall the merged cell and the surrounding cells 8. update the cell chained list 9. solve Poisson’s equation
7.3.4. Refinement/merging criterion As for classical AMR [5] the refinement and merging criteria can be either physical or based on Richardson extrapolation. For gas discharge modeling, a simple criterion based on physical considerations is used. The goal of the criterion is to capture the plasma zone along with its sheaths. We showed in our previous paper [4] that the mesh convergence occurs in this type of problem when the Debye length is resolved. So we defined a pseudo Debye length by:
δˆ =
2ε0 k B T e e 2 (ne + nn + n p )
.
(31)
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
257
This length is roughly equal to the Debye length within the plasma where ne + nn n p (quasi-neutrality) and is also quite small in the sheath region where the positive ion density is high and electron density is low. It is mandatory to properly refine the sheath area because all the ionization occurs there (in the Local Field Approximation).
• Refining criterion for the cell:
δˆ Cell size.
(32)
• Merging criterion on the average pseudo Debye length for the cells to merge:
ˆ 3 (Cell size). δ
(33)
Note. The merging criterion has been numerically adjusted to avoid any unnecessary refining/merging toggling of the same cell. It corresponds to a 50% hysteresis. 8. Simulation results
Fig. 18. Discharge current sensitivity to the mesh with the asynchronous scheme.
8.1. Reference test case Simulation parameters:
• • • • • •
Gas: Air Electron temperature: T e = 1 eV Ion temperature: T p = T n = 0.01 eV Applied voltage step: +1200 V Size: h = 150 μm, w = 50 μm, L = 400 μm Grid size: 512 × 256 for reference case and somewhere between 3188 and 512 × 256 for AAMR • Dielectric relative permittivity: εr = 10 • Flux scheme: The drift part of the flux is computed using MUSCL reconstruction with the minmod limiter, the diffusion part is computed using centered differences. • Simulated time: 150 nanoseconds. The upper electrode tip is sharp, and therefore the electric field is enhanced in this region. That is why we used a pre-refined mesh near the electrode tip. This is done by defining a maximum size for each cell which depends on its location with respect to the upper electrode. This “geometrical” mesh can be seen in Fig. 16 at t = 0 ns. In the early stages of the discharge formation, the mesh hardly differs from the geometrical mesh. As the plasma forms and starts to expand along the surface, the cell number grows and the refined area tends to follow the expanding plasma tip. Later in time, the dielectric layer becomes charged, the potential difference in the gas gap decreases, and the plasma starts to decay by recombination and ambipolar diffusion. Meanwhile, the cells in the mesh are starting to merge. Fig. 17 shows the density of each species compared to the reference case computed with the asynchronous scheme. Let us evaluate the performances of the AAMR solution with respect to the asynchronous solution in terms of accuracy and computation time. The AAMR solution captures well the trends in the plasma evolution in low density areas but differs in magnitude from the reference case. However, in the plasma zone the solutions agree. This deviation of the AAMR solution from the reference solution in low density areas is a consequence of the refinement strategy. In effect the refinement criterion chosen does not focus on low density areas where there is no plasma. The refinement/merging criteria and thresholds have been chosen in order to minimize the computational time by looking at the reference solution for macroscopic quantities such as current. Indeed from a macroscopical point of view both discharges are very similar, as can be seen in Fig. 14 where the currents from each case
Fig. 19. Discharge current sensitivity to the mesh with AAMR.
are plotted. As far as CPU time is concerned the gain can be correlated to the number of cells and for this example is more than a factor ten (see Fig. 15). 8.2. Mesh convergence When using the asynchronous scheme with an uniform mesh, it was not possible to reach mesh convergence due to memory limitation. However, Figs. 18 and 19 show that the current of the discharges are very similar for AAMR and the reference case with the same equivalent resolution. Nevertheless additional noise is present in the current calculated with the AAMR scheme. With AAMR, mesh convergence could be reached about an equivalent grid resolution of 1024 × 512. 8.3. Multi-scale convergence With the help of AAMR, a so-called multi-scale convergence test can be constructed. The idea is to keep constant the physically significant parameters of the reference test such as voltage, dielectric layer width and relative permittivity, while progressively increasing the gas volume which is taken into account. The scaling has been
258
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
Fig. 22. CPU time vs simulated time for various scale parameter L.
Fig. 20. Schematic view of the simulated domains with scale parameter L.
for the simulation of discharges at atmospheric pressure where the range of time and length scales is problematic, and is particularly suitable for the modeling of sharp localized phenomena such as streamers. AAMR has been applied here to the modeling of a surface dielectric barrier discharge, and we find a gain of CPU time of more than a factor of 10 with respect to the asynchronous scheme for an equivalent precision. This represents roughly more than a factor of one hundred with respect to a classical explicit integration method with a unique integration time step and an uniform grid. Coupling of transient plasma and Navier–Stokes equations for a complete description of plasma–flow interaction should be feasible using AAMR. Work is in progress in this direction. Furthermore AAMR could be applied to various multi-scale transport problems in other fields. Acknowledgements This work was partially supported by the EADS foundation (one of the authors Thomas Unfer is supported by a fellowship of the EADS foundation). References
Fig. 21. Discharge current vs simulated time for various scale parameter L.
done with the length parameter L. L 0 is the reference case and L i +1 = 2L i with the notations of Fig. 20. For the various cases the size of the smallest cell is kept constant to maintain the same level of resolution, whereas twice larger cells may be used for the low density areas in L i +1 with respect to L i . As far as discharge current is concerned, the results are very good. In fact the differences are lower than the noise level, as shown in Fig. 21. Furthermore, AAMR performs impressively in term of CPU time when increasing the scale parameter (see Fig. 22). This test shows the full power of the method. It is another step forward from the asynchronous scheme towards multi-scale multi-physical modeling of sharp localized phenomena. 9. Conclusion We have shown in this paper that the Adaptive Mesh Refinement (AMR) technique described here, coupled with the asynchronous time integration method described in a previous paper, is a very promising technique for multi-scale gas discharge problems. The resulting Asynchronous Adaptive Mesh Refinement (AAMR) method proposed in the present paper can be extremely efficient
[1] J.R. Roth, Aerodynamic flow acceleration using paraelectric and peristaltic electrohydrodynamic effects of a one atmosphere uniform glow discharge plasma, Phys. Plasmas 10 (2003) 2117–2126. [2] J. Pons, E. Moreau, G. Touchard, Asymmetric surface dielectric barrier discharge in air at atmospheric pressure: Electrical properties and induced airflow characteristics, J. Phys. D: Appl. Phys. 38 (2005) 3635–3642. [3] C. Montijn, W. Hundsdorfer, U. Ebert, An adaptive grid refinement strategy for the simulation of negative streamers, J. Comput. Phys. 219 (2006) 801–835. [4] T. Unfer, J.P. Boeuf, F. Rogier, F. Thivet, An asynchronous scheme with local time stepping for multi-scale transport problems: Application to gas discharges, J. Comput. Phys. 227 (2007) 898. [5] M. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys. 53 (1984) 484–512. [6] D. Burgarelli, M. Kischinhevsky, R.J. Biezuner, A new adaptive mesh refinement strategy for numerically solving evolutionary PDE’s, J. Comput. Appl. Math. 196 (1) (2006) 115–131. [7] Y. Lagmich, T. Callegari, L. C. Pitchford, J.-P. Boeuf, Model description of surface dielectric barrier discharges. [8] Bolsig+, www.laplace.univ-tlse.fr/-BOLSIG-Resolution-de-l-equation-Valeurs (numériques utilisées dans le modèle). [9] G.J.M. Hagelaar, L.C. Pitchford, Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models, Plasma Sources Sci. Technol. 14 (2005) 722–733. [10] H. Samet, The quadtree and related hierarchical data structures, Comput. Surv. 16 (2) (1984). [11] Y. Saad, M. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (3) (1986).
Contents lists available at ScienceDirect
Computer Physics Communications www.elsevier.com/locate/cpc
Multi-scale gas discharge simulations using asynchronous adaptive mesh refinement Thomas Unfer a,b,c,∗ , Jean-Pierre Boeuf a,b , François Rogier c , Frédéric Thivet d a
Université de Toulouse, UPS, INPT, LAPLACE (Laboratoire Plasma et Conversion d’Energie), 118 route de Narbonne, F-31062 Toulouse cedex 9, France CNRS, LAPLACE, F-31062 Toulouse, France c ONERA, 2 avenue Edouard Belin, 31400 Toulouse, France d Université de Toulouse, ISAE, 10 avenue Edouard Belin, 31055 Toulouse cedex 4, France b
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 2 July 2008 Received in revised form 22 September 2009 Accepted 22 September 2009 Available online 3 October 2009 Keywords: Multi-scale Asynchronous Transport problem Local time stepping Gas discharges CFL condition Adaptive mesh refinement
The breakdown of a gas gap at high products of pd (pressure × distance) is a multi-scale phenomenon in both time and space. This is especially true when the plasma is interacting with a gas flow, a problem of considerable recent interest in the context of aerodynamic applications of surface discharges. This paper presents a contribution to the numerical modeling of such discharges. We describe here a new approach for adaptive meshing which is suitable for use with the explicit asynchronous integration scheme described in our previous publication. Rather than relying on a family of nested grids as is commonly done, this technique is based on a single unstructured mesh with possible non-conforming cells at the interface between coarse and fine areas. Substantial computational time saving has been achieved for a surface dielectric barrier discharge configuration of the kind proposed as plasma actuators for flow control. © 2009 Elsevier B.V. All rights reserved.
1. Introduction In the recent years, substantial research efforts have been carried out throughout the world aiming at controlling aerodynamic air flow with plasma actuators [1,2]. Surface dielectric barrier discharges (DBDs) have been proposed as plasma actuators able to modify the boundary layer of a flow around an airfoil. The numerical simulation of these plasma actuators is difficult because of the space and time scales which must be taken into account span many orders of magnitude. One useful approach is adaptive meshing, and Adaptive Mesh Refinement (AMR) techniques have already been successfully applied in the field of gas discharges to investigate the physics of the propagation of fast ionization waves known as streamers [3]. In a previous paper [4], we derived an asynchronous integration procedure for the charged particle transport equations on which fluid models of gas discharge phenomena are based. This technique reduces both numerical diffusion and CPU time and is particularly adapted for meshes with strong variation in cell size. When Berger and Oliger first developed AMR for block-structured meshes in the
*
Corresponding author at: Université de Toulouse, UPS, INPT, LAPLACE (Laboratoire Plasma et Conversion d’Energie), 118 route de Narbonne, F-31062 Toulouse cedex 9, France. E-mail address: [email protected] (T. Unfer). 0010-4655/$ – see front matter doi:10.1016/j.cpc.2009.09.017
©
2009 Elsevier B.V. All rights reserved.
mid 1980’s [5], they introduced a family of nested grids with time synchronization between large and fine grids based on fractional time steps. The time steps on the finer grids were chosen so as to keep the CFL number constant with respect to the larger grids. In this paper the concept of a family of nested grids is not used. In the proposed Asynchronous Adaptive Mesh Refinement (AAMR) method the family of nested grids can be replaced to advantage by a structure introduced by Burgarelli et al. in [6], namely an Autonomous Leaves Graph (ALG). This paper is organized as follows: Section 2 briefly explains the motivation of this work. The model for a plasma discharge in air is presented in Section 3, and the principles of the asynchronous time integration are summarized in Section 4. Section 5 introduces the specific data structures needed for AAMR, and Section 6 explains how the asynchronous scheme is to be used with an unstructured mesh with non-conforming cells. The AAMR algorithm is then introduced in Section 7, and some simulation examples from numerical calculations and comparisons of discharges in the typical surface DBD configuration used for aerodynamic flow control are presented in Section 8. 2. Motivation of this work There is a growing interest in numerical simulations of the generation and evolution of atmospheric pressure air plasma and,
248
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
especially, in the interaction of the plasma with a gas flow. These phenomena are very demanding in terms of CPU time because of the large variation in time and length scales that are involved in the physics of the plasma and plasma/flow interaction. For example, the time scale associated with the coupling between the charged particle transport equations and Poisson’s equation is known as the dielectric relaxation time τd , and a condition for numerical stability in explicit integration schemes is that τd not be exceeded. In atmospheric pressure discharges τd can be as small as 10−12 s. The CFL constraint associated with the electron transport depends on the grid size but, in practice, can be lower than 10−12 s. One concept of actuator for aerodynamic flow control is an asymmetrical surface Dielectric Barrier Discharge (DBD) with a sinusoidal voltage waveform in the kHz frequency range. Momentum transfer from the ions to the neutral gas molecules induces a net gas flow on a time scale 10−4 to 10−3 s. Consequently the time scales involved in this problem extend over more than ten orders of magnitude. The problem is also covering space scales of various orders of magnitude since the Debye sheath in atmospheric plasmas can be as small as a few micrometers whereas the plasma extension of a filamentary discharge above a surface can be several centimeters generating an airflow over a few tens of centimeters. As a result the physics of the discharge is truly multi-scale and it is therefore essential for modeling to develop dedicated numerical tools that are able to tackle these issues. We show in this paper that associating the asynchronous time integration method described previously in [4] with adaptive mesh refinement techniques can lead to excellent synergy and considerable computation time saving. 3. A gas discharge model for air A minimal air plasma model is considered, consisting of only three species, electrons (e), one type of positive ion (p) and one type of negative ion (n). The chemistry accounts for only four reactions: direct ionization by electron impact, electron attachment to form a negative ion, and electron/positive ion and negative ion/positive ion recombination. A detailed description of this model can be found in [7]. While this plasma chemistry is admittedly quite simple compared to a real situation, we expect that the simple model will capture the essential physics of the plasma/flow coupling. The discharge model is based on the following fluid equations:
• Continuity equations
∂ ne Γe = (α − η)Γe − re− p ne n p , + ∇. ∂t ∂np Γp = α Γe − re− p ne n p − rn− p nn n p , + ∇. ∂t ∂ nn Γn = ηΓe − rn− p nn n p , + ∇. ∂t
(1) (2) (3)
ne , n p , nn are respectively the electron, positive ion and negative ion densities. Γe , Γp , Γn are the corresponding particle fluxes. α is the ionization coefficient (first Townsend coefficient), η is the attachment coefficient. Volume recombination is described using re− p and rn− p which are the electron/positive ion and negative ion/positive ion recombination coefficients. • Momentum equation in drift-diffusion form
kB Te ∇ ne , Γe = μe ( E ) −ne E − e kB T p np , ∇ Γp = μ p ( E ) n p E − e
(4) (5)
Fig. 1. Simulation domain for surface dielectric barrier discharge for flow control.
kB Tn nn , ∇ Γn = μn ( E ) −nn E − e
(6)
μe , μ p , μn are respectively the electron, positive ion and negative ion mobilities and are function of the local electric field E only. The diffusion coefficients are expressed as a function of the various mobilities and temperatures using the Einstein relation. In practice the electron temperature T e is taken about 1 eV and the ion temperatures T p , T n are taken equal to the background gas temperature. • Poisson’s equation ε E ) = e (n p − ne − nn ), ∇.(
(7)
ε is the dielectric permittivity of the media (air or dielectric layer), e is the elementary charge. In collisional discharges, and this is particularly true at atmospheric pressure, the electron energy transport is not very important. In fact the energy gained by the charged particles in the electric field tends to be locally balanced by the losses due to collisions with the neutral gas. Under these conditions, the electron energy equation can be replaced by the so-called Local Field Approximation (LFA). The loss of precision is minimal for conditions of interest and considerable computational time can be saved. With the LFA, the transport and rate coefficients (mobilities, reaction coefficients) depend the local value of electric field (or, more precisely, on the reduced electric field, E / N, the ratio of the electric field strength to the neutral density), and their values are obtained by solving the steady state and homogeneous (no spatial gradients) Boltzmann equation. Transport and rate coefficients are calculated using the freeware BOLSIG+ [8,9], and the data tables in tabular form are input to the code. We consider the same numerical experiment as in our previous paper [4]. The geometry is reproduced in Fig. 1. The charging of the dielectric layer by the charged species and the induced potential due to the surface charge is taken into account. 4. The principles of the asynchronous scheme 4.1. Asynchronous integration The asynchronous scheme is a finite volume scheme for unsteady calculations, which is particularly suitable for multi-scale transport problem. The details of the algorithm have been published in [4]. For convenience the principles of the asynchronous time integration are briefly summarized in this section for a onedimensional situation. The extension to multiple dimensions is obvious:
∂ n ∂Γ + = S. ∂t ∂x
(8)
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
249
The asynchronous scheme assumes that each flux or source term is updated independently according to a refresh time tag with respect to a global simulation clock. ni are the density values at the cell centers xi = i x, v i + 1 are the velocity values taken at the cell 2
interfaces xi + 1 = (i + 12 )x. For the scheme to be stable, each local 2 time step is limited by the local value of the CFL condition. Each conservative variable has its own value, value time tag and variation rate. The variation rate represents the instantaneous balance of fluxes and source term of the conservative variable. The idea is to let each component of this balance to evolve independently in time:
Γi + 1 Γi − 1 ∂n S ∂n Γ + ∂n Γ − ∂n 2 2 = Si − + = + + . ∂t i x x ∂t i ∂t i ∂t i Different time variables are needed: t simulation : the current discrete simulation time, t ni : time tag associated with the stored value of density ni , t iΓ : refresh time tag of flux Γi , t iS : refresh time tag of source term S i . During the simulation, t simulation jumps discretely from the most urgent refresh time tag to the next most urgent. Initialization 1. Initialize all density and fluxes. 2. Initialize all refresh time tags to the initial time. Proceed until the simulation time is completed 1. Find the most urgent flux or source term to be refreshed. 2. t simulation becomes this most urgent refresh time tag. 3. Compute the value of the density on both side of the interface needed for flux computation (only a few cells are involved) or of the cell for source term computation: • n j = n j + ∂∂nt (t simulation − t nj ),
• t nj = t simulation with j ∈ [i , i + 1] for first order flux, j ∈ [i − 1, i + 2] for second order flux (e.g. MUSCL) or j = i for source term.
4. Compute the new value of the flux or source term and update its variation rate according to:
•
∂n Γ + ∂t i ∂n Γ − ∂ t i +1
=−
Γi + 1
2 x , Γi + 1
S
or ∂∂nt i = S i . 5. Compute the new refresh time tag of this flux or source term according to the local CFL condition or reaction kinetic t iΓ = t iΓ + t CFL when updating a flux or t iS = t iS + t S when updating a source. 6. If the solution is to be monitored before the next most urgent refresh time tag, build the solution at the output time tag and store it.
•
output ni
=
2
x
∂ n S ∂ n Γ + ∂ n Γ − output = ni + + + t − t ni ∂t i ∂t i ∂t i
for all cells. 4.2. Coupling with Poisson’s equation Poisson’s equation is solved at each cell center with time steps corresponding to the global dielectric relaxation time (minimum of all the local dielectric relaxation times). The electric field is defined on each interface between two cells and is advanced in time between updates of the solution of Poisson’s equation using the following estimator with j c , the conduction current going through
Fig. 2. Flux data type.
the interface:
∂ E i , j jt −t − jtc (t ) = c . ∂t ε0
(9)
Some justification of this coupling can be found in [4]. 5. Specific data structures for AAMR The basis of gas discharge modeling is to integrate in time the transport equations for the charged particles coupled with Poisson’s equation for the electric field. This essentially implies solving repetitively an elliptic equation (Poisson’s equation). For purposes of accuracy when solving Poisson’s equation, it is important that neighboring cells do not differ in size by more than a factor two. This requirement can be fulfilled if, for instance, the mesh has a quadtree structure as defined in [10]. The idea here is to create a simplified ALG (Autonomous Leaves Graph) structure, mimicking a quadtree. Because the full features of ALG are not needed to generate a quadtree-like mesh, a simpler structure can be defined which is expected to speed somewhat the computations. The resulting data structure is suitable for both asynchronous time integration of transport equations and local refinement/merging operations over an unstructured mesh. Moreover, because the asynchronous scheme keeps jumping from one memory area to another while updating interface fluxes, the memory cache misses are much more likelier to occur for asynchronous time integration than for a standard integration scheme. In order to mitigate this effect the concept of collocated data (i.e. keeping all the data corresponding to one single cell in the same memory location) is essential for asynchronous time integration. The useful data can be split in two parts: cell-centered data (densities, for instance) and interface data (flux). Consequently two data types have been defined, one gathering all interface data, termed “Flux data type”, and the other gathering the cell-centered data, or “Cell data type”. 5.1. Flux data type The Flux data type consists of all the data required for computation at the interface between two cells for each transport equation (e.g., fluxes and electric field). It points towards two cells of the Cell data type, one being the upstream cell called Cell_out in Fig. 2 and the other being the downstream cell called Cell_in in the same Fig. 2. The flux of a conservative variable is counted positively if it goes from Cell_out to Cell_in. 5.2. Cell data type The Cell data type consists of all the cell-centered data for each equation (value, variation rate, time of value, etc.) plus additional parameters (boundary, Poisson matrix coefficients, etc.). It also contains additional information about the geometry, position and size of the cell, level of refinement and coordinates in a virtual uniform mesh. In two dimension the Cell data type points towards at least 4 Flux objects (Flux_East_A, Flux_South_A, Flux_West_A, Flux_North_A) when the cell is connected to other cells of the same size or greater. Up to 4 additional Flux objects pointers can be used (Flux_East_B, Flux_South_B, Flux_West_B, Flux_North_B) when the cell is connected to smaller cells (see Fig. 3). The Cell data type also contains additional pointers to the corresponding actions which are used for the asynchronous integra-
250
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
Fig. 3. Cell data type in 2D. Fig. 5. Schematic view of a local refinement.
Fig. 6. Wrong AAMR scenario.
Fig. 4. Schematic view of a regular mesh.
tions (see Section 7.1, these pointers are used to recall an action after an AAMR event so that the fluxes can be refreshed immediately). The Cell data type is additionally connected to the next and previous cell so as to have a chained list of the cells within the mesh (see Section 7.2). 5.3. The mesh: A composite object structure The mesh is built using fluxes to connect cells and is an Autonomous Leaves Graph. By construction, the fluxes are oriented along the positive directions of the mesh. A regular mesh only uses the “_A” type pointers to flux as shown in Fig. 4, whereas the “_B” type pointers to flux are used to connect cells of different sizes around some refined zone (see Fig. 5). Boundary cells can be connected to a single ghost cell. This object-like construction of the mesh satisfies the two following properties:
• Any block refined mesh that can be generated using classical AMR can be reproduced using a single unstructured mesh using this construction. • Any further refinement or merging can be done while modifying only the storage of the local elements in memory and has no effect on the rest of the mesh. 5.4. Adaptive meshing rules Two principles are postulated to develop the algorithm. These are, at each time, the mesh must preserve proper cells connectivity and cells may differ in size from their neighbors by no more than a factor two. To satisfy these principles a few simple rules can be formulated to avoid any “illegal” refinement or merging to occur.
• AAMR Rule #1: In N dimension a cell can be refined in 2 N if its size is larger than or equal to the size of any of its direct neighbors. (Direct neighbor means that the cells are connected by a flux, N is the dimension of the simulation.)
• AAMR Rule #2: 2 N cells of identical size can be merged if their size is smaller than or equal to the size of any of their direct neighbors. Note. These are the minimal rules. However, in practice and for precision purposes, especially when coupling drift-diffusion equations with Poisson’s equation which is our present concern, “direct neighbors” can be extended to “neighbors” which means taking into account the 2 N diagonal neighbors (North–East, South–East, South–West, North–West in 2D). This prevent cells from differing in size of a factor 4 with their diagonal neighbors. This feature is required for accurately solving Poisson’s equation. These two rules are sufficient to insure that the mesh integrity is preserved and that the neighboring size principles are not violated. However it is not sufficient to insure that the adaptive meshing process is reversible, the AAMR process may get trapped in some configuration. Fig. 6 shows an example in 1D which explains why things can go wrong without an additional rule. In fact the third AAMR event in the figure, i.e. “the merging of cell 4 & 5” is problematic because afterward cell 3 or cell 6 can no longer be merged. This behavior shall be forbidden. Each cell shall “remember” a virtual parent cell of greater size in case of merging so as to make the AAMR process reversible. This can be done by assigning coordinates to each cell corresponding to its coordinates in an hypothetical regular mesh of the same size over the whole domain and following Rule #3:
• AAMR Rule #3: Only a cell with odd coordinates is available for merging with cells with the successive even coordinates. 6. Using the asynchronous scheme over a non-conforming unstructured grid 6.1. Mass conservation Consider now the asynchronous scheme in a finite volume formulation with no source term and zero flux at the boundaries. The
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
251
Fig. 8. Time/space diagram of the asynchronous scheme assuming positive drift velocity.
• First case: t p tq t p +1 tq+1 t q +1
ni
t
= ni q −
i =1
N Cell N faces(i ) ∂n Γ t ni + t ∂ t i→ j j =1
i =1
N Cell
nti
=
(10)
i =1
because by definition of the algorithm:
t q +1
ni
Γi → j ∂n Γ ∂n Γ =− =− . ∂ t i→ j x ∂ t j →i
(11)
N Cell is the total number of cells, N faces(i ) is the total number of faces of the cell i. Fig. 7 shows the electron density relative variation in the numerical experiment described in Section 8 with the following hypotheses: both the electron source term and the electron fluxes at the boundaries are taken equal to zero. This test basically evaluates the numerical performance in term of mass conservation under drift diffusion in an electric field with a non-conforming mesh. Density is conserved up to round-off error precision.
The demonstration is given in one dimension for the advection equation for the sake of simplicity. With the notations of Fig. 8, suppose that the drift velocity is smooth and the local time steps are proportional to the local CFL condition, i.e.:
kCFL x v i− 1
kCFL
x
kCFL
x
2
t p =
kCFL x v i+ 1 2
= tq
1−
v
v
∂ v x . 1+ ∂x v
∂ v x kCFL ∈ ]0, 1[, ∂ x 2v
∂ v x 1+ ∂ x 2v
t = ni q v i
(12)
+ ni
t
− ni q ∂(nv ) ∂ n (t p − tq ) + = −v i τ. tq ∂x ∂t tq x
(14)
Hence there is consistency since t p − tq = O (tq x) and
τ = O (tq ).
• Second case: t p tq t p +1 t p +2 tq+1 t q +1
ni
t
1
t
t
t
p p +1 tq ni q v i + 1 − τ ni − v − t p ni − v 1 i − 12 1 i − 12 2 t p +2 , − (tq − τ − t p )ni − v 1 i− 1
= ni q −
x
2
∂ v x ∂n ∂n − v i x + v i (τ + t p ), ∂x 2 ∂x ∂t
v i ∂n ∂(nv ) tq (t p − tq )τ = ni − tq + ∂x tq x ∂ t + t p τ + t p (t p − tq ) . (15)
t p +2 ni v i− 1 2
6.2. Consistency under local CFL condition
tq =
∂ v x , ∂x 2 ∂n ∂n tp tq tq ∂ v x ni v i − 1 = ni v i − ni − v i x − v i (t p − τ ), 2 ∂x 2 ∂x ∂t ∂n ∂n t p +1 tq tq ∂ v x ni v i − 1 = ni v i − ni − v i x + v i τ , 2 ∂x 2 ∂x ∂t
∂ n (t p − tq ) ∂(nv ) t q +1 tq ni = ni − tq + vi τ , ∂x ∂t tq x tq ni v i + 1 2
N Cell
i =1
t
2
scheme is clearly conservative:
nti +t =
t
p tq ni q v i + 1 − τ ni − v 1 i − 12 2 t p +1 − ni − v (tq − τ ) , 1 i− 1
Fig. 7. Mass conservation test.
N Cell
1
x
t q +1
ni
t = ni q v i
t − ni q
In this case, necessary gives consistency.
τ tq − t p so τ = O (tq x), which
• Third case: t p tq tq+1 t p +1
t q +1
ni (13)
For x v ∂∂ vx the following condition is imposed on the local time steps: |t p − tq | tq , this means that three cases have to be considered:
t q +1
ni
t tp tq ni q v i + 1 − tq ni − , v 1 i − 12 2
∂ n t p − τ ∂(nv ) t . = ni q − tq + vi ∂x ∂ t x t
= ni q −
1
x
(16)
In this case, necessary tq τ t p so t p − τ = O (tq x), which gives consistency.
252
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
6.4. CFL condition on a conforming and a non-conforming cell
Table 1 Unstable counterexample in 2D with a 20% desynchronization. Time
0
Density Rate Flux #1 Flux #2
0.4
0.5
0.8
−2
0.2 −1.2
−0.004
1 1
1 0.2
0.08 0.28 0.08 0.2
1
With the notation of Fig. 3, let us consider, for instance that Flux_East_A is positive. Consequently its refresh rate is function of all the fluxes emptying the cell and is equal to the cell time step:
kCFL
t Flux = t Cell = Faces
As a result the first-order asynchronous scheme is convergent. The local CFL condition hypothesis is essential for the convergence of the scheme. 6.3. An alternative formulation of the asynchronous scheme We showed in our previous paper [4] that each local contributor (fluxes and source terms) to the variation rate of the density can be updated asynchronously, and this principle has been used successfully for 2D calculations. However, it is possible to build a small counterexample showing that the stability of the scheme is not assured for arbitrary values of the initial desynchronization between the fluxes. If one considers a unit-sized cell with 20% desynchronization between two fluxes emptying the cell. The local time steps last both 0.5 and the desynchronization is 0.1. Table 1 shows the computation. Furthermore depending on the problem, this highly segregated treatment may not be the optimum in term of computational time or memory load. If any flux or source is independent, it means managing a corresponding action in the scheduler. The algorithm can be optimized using the fact that the CFL condition on a single flux is dependent on the fluxes emptying the same cell in the other directions. That is why the CFL condition can be seen as a condition on the upwind cell. Actually a flux shall be refreshed more often than it takes for the upwind cell to be emptied. Consequently it makes sense to synchronize all the fluxes which are emptying a single cell. For gas discharges simulations with a formulation similar to the model defined in Section 3, all the source terms can also be synchronized with the electron fluxes because the formulation of the chosen chemistry is expressed with respect to the electron flux using Townsend coefficients (see [7]). As a result, a single asynchronous action is required for each specie in each cell. Nevertheless if at a certain point of the simulation, a numerical flux changes direction, it has to be synchronized with the cell that is now being emptied by this flux. The asynchronous action of this cell shall be “recalled”, i.e. this cell shall be refreshed immediately. The alternative formulation makes the scheme faster because less actions are used in the Discrete Time Scheduler leading to less memory usage. Furthermore the stability condition in multiple dimensions can be proved in the case of the first order upwind asynchronous scheme if the leaving fluxes are synchronized. Supposing v x > 0 and v y > 0, the stability condition is equal to the local CFL condition with the notation of Fig. 8:
t q +1 ni , j
t = ni ,q j
1
−
x
t
t
tq ni ,q j v ix+ 1 + tq ni ,q j v 2
t q +1
n ( z )i j v x
−
i − 12
tq t q +1
ni
t
ni q
t q +1
dz −
n ( z )i j v
y i + 12
y
kCFL is the stability margin (kCFL ∈ ]0, 0.5[ for MUSCL scheme). The fluxes contribution are prescribed as follows, f CFL Flux = 0 if the actual numerical flux is directed towards the cell and f CFL Flux = V CFL when the flux is directed outwards. The length of the cell Cell Length is taken equal to the cell size if there is only one flux per cell face or half the size if there are two. Additional stability can be introduced by defining an equivalent CFL condition for the source term with f CFL Source Term being the stability constraint of the unstable part of the source term (negative part). This definition is equivalent to the classical multiple dimensions time step for a synchronous explicit scheme. For instance with v x > 0 and v y > 0, vx vy , f CFL Flux y − = 0, f CFL Flux x+ = and f CFL Flux x− = 0, f CFL Flux x+ = x y x y
the time step is t Cell = kCFL v x y + v y x . If during the simulation the sign of the flux changes then it recombines its refresh rate and CFL contribution with the other cell and this cell is refreshed immediately. For gas discharge modeling the equivalent CFL velocity used is:
v CFL = μ( E ) E +
D
xC 1 − xC 2
+ x
n i μi
ε0
.
(19)
D is the diffusion coefficient, the two factor in front of the x2 ) can be dropped because classical diffusion time step (t = 2D the influence of the diffusion is taken into account for each fluxes which are actually emptying the cell. So for a purely diffusive problem, the local time step used are twice larger except in the cell holding the maximum density. xC 1 − xC 2 is the distance between the two cell centers, x is the cell size. When an AAMR event occurs, the new cells and their direct neighbors have to be refreshed immediately to insure the stability of the scheme, because the scheduled actions for the neighbors were planned assuming the previous configuration of the mesh. 6.5. Numerical diffusion When writing second-order Taylor development of the firstorder asynchronous scheme of the advection equation and isolating the numerical diffusion coefficient, one obtains in the first case:
(t p − τ )2 τ + τ 2 (tq − τ ) 2 2tq x (t p − tq )τ x . + + tq 2v
D async = v 2 −
tq
−v
(20)
The first and last terms also appear when calculating the numerical diffusion coefficient for the synchronous scheme giving x (1 − vxt ). In contrast for the asynchronous rise to D sync = v 2 scheme those two terms cancel each other at least at first order in x when kCFL = 1. The third term is also second order in x, only the second term as to be evaluated knowing that to first order tq t p and (tq − τ )τ 4 q with 0 τ tq . Then the asynchronous diffusion coefficient is lower than:
dz ,
D async
tq
tq x tq y 1− v 1 − v 1 . x i + 2 x i + 2
(18)
,
t 2
i − 12
f Flux + f Source
(17)
v x 8
.
(21)
This means that if the velocity variation over the mesh is greater than 25%, the asynchronous scheme becomes less diffusive
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
253
Fig. 9. Enhanced Discrete Time Scheduler. Fig. 10. Poisson’s equation discretization on non-conforming cell.
than its synchronous equivalent. From a physical point of view, the numerical diffusion in the synchronous scheme can be seen as a consequence of the fact that in the area where the CFL number is lower than one, information travels faster than it should. The asynchronous scheme cancels this effect to first order because information cannot go faster than the fastest local wave corresponding to the local CFL. However the asynchronous scheme introduces diffusion due to the local desynchronization τ (second term in Eq. (20)). Notes.
• Case one and three are corner cases. Nevertheless, to first order, they lead to no numerical diffusion coefficient of zero.
• The estimation of the asynchronous numerical diffusion coeft ficient correspond to the worst case τ = 2 q . Actually τ keeps varying in time so the average value of numerical diffusion should be lower.
6.6. Extended limiter function For the non-conforming cell, the definition of the slope used in the limiter function has to be extended. See the density splitting function in Section 7.3.2 for the details. 7. The asynchronous adaptive mesh refinement algorithm 7.1. Extension of the asynchronous scheme The asynchronous scheme is compatible with adaptive meshing because it is based on interface synchronization. However, the stability condition must be insured for the new cells and their neighboring cells after an AAMR event. Because the neighboring cell refreshment was planned according to the previous mesh configuration, the scheduling has become obsolete and needs to be updated after an AAMR event. The corresponding actions have to be recalled so that the cell can be refreshed immediately. Hence, the Discrete Time Scheduler has to manage the additional functions of task recalling, task creation, task suppression. This is accomplished by adding new pointers into the structure to retrieve the pending actions from the mesh and the task from the pending action. (See Fig. 9.) 7.2. Poisson’s equation Poisson’s equation has to be solved on the unstructured mesh. Using the notation of Fig. 10, the equation is discretized classically for each cell using a finite volume formulation of the electric field = e (n p − ne − nn )). The electric field flux conservation ( εE .dl
values at the interface have to be expressed in function of the potential values in the surrounding cells. For example on the upper east face of the cell O between point A and point B the potential derivative across the vertical face is obtained using:
∂V V East − V O VA − VB =α +β ∂x dEast− O AB with V A =
1 (V O 4
+ V North B + V North East + V East ), V B = V East + V East B ), α = 43 , β = − 13 .
(22) 1 (V O 3
+
When summing up the contributions of each face, one gets the full Poisson matrix. In a homogeneous area the scheme is equivalent to the standard five points scheme, whereas a thirteen points scheme is used for a cell completely surrounded by smaller neighbors. Unfortunately, Poisson’s matrix loses its convenient symmetry, so a conjugate gradient solver cannot be used. This non-symmetry can be seen, for instance, if one looks at the discretization of Poisson’s equation in cell O . The discretized equation in cell O uses the potential value in the North East cell. However, because the discretization in the North East cell uses the standard five points scheme, it does not use the value of the potential in the cell O . So, the matrix is clearly non-symmetrical. We decided to solve Poisson’s equation with the restarted GMRES algorithm [11]. From a practical point of view, the Poisson matrix is stored twice in memory — once in the object-mesh as part of the cell data type and another time in a 13 × X array where X is the cell number with a correspondence look-up table. Any AAMR event will lead to a modification of the Poisson matrix within the objectmesh. As a result the array and look-up table have to be updated before solving the equation. The GMRES algorithm is used to iterate only within the array to speed up calculations because looping within the object-mesh is rather slow. 7.3. The asynchronous adaptive mesh refinement procedure 7.3.1. General algorithm With respect to the asynchronous algorithm (see Fig. 11), the AAMR algorithm consists of adding a testing subroutine which may trigger an AAMR event (see Fig. 12). There are many ways to implement this testing, for instance, some specific asynchronous task dedicated to the testing could be added into the Enhanced Discrete Time Scheduler (see Section 7.1). It makes sense to associate the AAMR testing frequency to the refresh rate of the electron fluxes. The detailed pseudo-code of the “Test for AAMR” function follows: 1. if the source term has been refreshed more than n (user defined) times since last AAMR testing of the current cell: (a) if local refinement criterion is satisfied for the current cell i. Refine the current cell
254
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
Fig. 11. Asynchronous algorithm overview.
Fig. 13. How to split cells in 2D.
(i) update the cell chained list (j) solve Poisson’s equation 2. else (a) for all the neighboring cells within the mesh (Ghost cell excluded) i. if a neighboring cell is greater than the current cell, refine this cell (b) try again to refine the current cell
Fig. 12. AAMR algorithm overview.
(b) else if the current cell is available for merging (AAMR Rule #3) i. if the local merging criterion is satisfied for the current cell A. Merge current cell with its neighbors ii. else refresh the cell (c) else refresh the cell 2. else refresh the cell 7.3.2. The cell refinement procedure The cell refinement procedure consists of replacing a cell with smaller cells when the refinement criterion is met. However, it may happen that a cell cannot by refined because one of its neighbors is of greater size. In such case, the AAMR algorithm shall start by refining the neighboring cell before even if this cell does not need to be refined. The cell refinement procedure is therefore recursive: 1. If the current cell can be refined (a) update the density of the neighboring cells (b) disconnect the current cell from the fluxes coming from its neighbors (c) suppress the cell action task (d) replace the cell with 2 N new cells and fill them with the density splitting function, fill the Poisson matrix, etc. (e) create the new inner fluxes and possible external fluxes (of the _B type for the neighbors) (f) re-connect the fluxes (g) create the new cell action tasks (h) recall the surrounding cells
Note. Step (j) is obviously specific to gas discharge modeling. When Poisson’s equation is solved following an AAMR event, the next update of the electric field can be postponed by one dielectric relaxation time in the asynchronous loop. In case of recursive calls of the refinement procedure, it is only necessary in principle to solve Poisson’s equation once. Density splitting function. This function shall split conservatively and positively the density from the coarse cell to the finer cells. Here is an example in 2D with the notation of Fig. 13. First define the slopes:
sx = minmod
n E − n O n O − nW
s y = minmod
,
dE O
dO W
nN − n O n O − n S
sx B = minmod
s y B = minmod
,
dN O
,
(23)
,
(24)
dO S
n E B − n O n O − nW B
,
dE B O
dO W B
nN B − n O n O − n S B
,
dN B O
dO S B
.
,
(25)
(26)
Then split the density:
n1 = n − s x n2 = n + s x
x 2
x
n3 = n − sxB n4 = n + sxB
2
+ sy
+ syB
x 2
x 2
y
− sy
y 2
y
− syB
(27)
,
2
2
,
(28)
,
(29)
y 2
.
(30)
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
This splitting is conservative because n1 2x
n3
x y 2
2
+ n4
x y 2
2
y 2
+ n2 2x 2y +
= nx y. The positivity in insured by the use
of the limited slopes. The limiter function for MUSCL scheme used in flux computation is extended to the non-conforming cells using the same
Fig. 14. Discharge current vs simulated time.
255
slopes. If a coarse cell is surrounded by finer cells as in Fig. 13, then the limiter uses the slope s x for the upper fluxes (E or W) and sxB for the lower fluxes (E_B or W_B). If for instance there is one single cell on the east side then the average slope sx +2sxB is used to estimate the density on the right side.
Fig. 15. CPU time and cell number ratio vs simulated time.
Fig. 16. AAMR mesh evolution with time, ith level represents the size of the local square cells with a size of
L 2(i −1)
.
256
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
Fig. 17. Comparison between AAMR (left) and asynchronous scheme (right) at t = 45 ns from top to bottom: electron density (log(m−3 )), positive ion density (log(m−3 )), negative ion density (log(m−3 )) and equipotential lines (V).
7.3.3. The cells merging procedure 1. update the cells to be merged (this may already be done while evaluating the merging criterion) 2. disconnect the cells to be merged from the fluxes connected to their neighbors 3. suppress the useless actions 4. suppress the useless fluxes and cells 5. fill the merged cell with the averaged values of the ex-finer cells (density, etc.) 6. reconnect the external fluxes 7. recall the merged cell and the surrounding cells 8. update the cell chained list 9. solve Poisson’s equation
7.3.4. Refinement/merging criterion As for classical AMR [5] the refinement and merging criteria can be either physical or based on Richardson extrapolation. For gas discharge modeling, a simple criterion based on physical considerations is used. The goal of the criterion is to capture the plasma zone along with its sheaths. We showed in our previous paper [4] that the mesh convergence occurs in this type of problem when the Debye length is resolved. So we defined a pseudo Debye length by:
δˆ =
2ε0 k B T e e 2 (ne + nn + n p )
.
(31)
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
257
This length is roughly equal to the Debye length within the plasma where ne + nn n p (quasi-neutrality) and is also quite small in the sheath region where the positive ion density is high and electron density is low. It is mandatory to properly refine the sheath area because all the ionization occurs there (in the Local Field Approximation).
• Refining criterion for the cell:
δˆ Cell size.
(32)
• Merging criterion on the average pseudo Debye length for the cells to merge:
ˆ 3 (Cell size). δ
(33)
Note. The merging criterion has been numerically adjusted to avoid any unnecessary refining/merging toggling of the same cell. It corresponds to a 50% hysteresis. 8. Simulation results
Fig. 18. Discharge current sensitivity to the mesh with the asynchronous scheme.
8.1. Reference test case Simulation parameters:
• • • • • •
Gas: Air Electron temperature: T e = 1 eV Ion temperature: T p = T n = 0.01 eV Applied voltage step: +1200 V Size: h = 150 μm, w = 50 μm, L = 400 μm Grid size: 512 × 256 for reference case and somewhere between 3188 and 512 × 256 for AAMR • Dielectric relative permittivity: εr = 10 • Flux scheme: The drift part of the flux is computed using MUSCL reconstruction with the minmod limiter, the diffusion part is computed using centered differences. • Simulated time: 150 nanoseconds. The upper electrode tip is sharp, and therefore the electric field is enhanced in this region. That is why we used a pre-refined mesh near the electrode tip. This is done by defining a maximum size for each cell which depends on its location with respect to the upper electrode. This “geometrical” mesh can be seen in Fig. 16 at t = 0 ns. In the early stages of the discharge formation, the mesh hardly differs from the geometrical mesh. As the plasma forms and starts to expand along the surface, the cell number grows and the refined area tends to follow the expanding plasma tip. Later in time, the dielectric layer becomes charged, the potential difference in the gas gap decreases, and the plasma starts to decay by recombination and ambipolar diffusion. Meanwhile, the cells in the mesh are starting to merge. Fig. 17 shows the density of each species compared to the reference case computed with the asynchronous scheme. Let us evaluate the performances of the AAMR solution with respect to the asynchronous solution in terms of accuracy and computation time. The AAMR solution captures well the trends in the plasma evolution in low density areas but differs in magnitude from the reference case. However, in the plasma zone the solutions agree. This deviation of the AAMR solution from the reference solution in low density areas is a consequence of the refinement strategy. In effect the refinement criterion chosen does not focus on low density areas where there is no plasma. The refinement/merging criteria and thresholds have been chosen in order to minimize the computational time by looking at the reference solution for macroscopic quantities such as current. Indeed from a macroscopical point of view both discharges are very similar, as can be seen in Fig. 14 where the currents from each case
Fig. 19. Discharge current sensitivity to the mesh with AAMR.
are plotted. As far as CPU time is concerned the gain can be correlated to the number of cells and for this example is more than a factor ten (see Fig. 15). 8.2. Mesh convergence When using the asynchronous scheme with an uniform mesh, it was not possible to reach mesh convergence due to memory limitation. However, Figs. 18 and 19 show that the current of the discharges are very similar for AAMR and the reference case with the same equivalent resolution. Nevertheless additional noise is present in the current calculated with the AAMR scheme. With AAMR, mesh convergence could be reached about an equivalent grid resolution of 1024 × 512. 8.3. Multi-scale convergence With the help of AAMR, a so-called multi-scale convergence test can be constructed. The idea is to keep constant the physically significant parameters of the reference test such as voltage, dielectric layer width and relative permittivity, while progressively increasing the gas volume which is taken into account. The scaling has been
258
T. Unfer et al. / Computer Physics Communications 181 (2010) 247–258
Fig. 22. CPU time vs simulated time for various scale parameter L.
Fig. 20. Schematic view of the simulated domains with scale parameter L.
for the simulation of discharges at atmospheric pressure where the range of time and length scales is problematic, and is particularly suitable for the modeling of sharp localized phenomena such as streamers. AAMR has been applied here to the modeling of a surface dielectric barrier discharge, and we find a gain of CPU time of more than a factor of 10 with respect to the asynchronous scheme for an equivalent precision. This represents roughly more than a factor of one hundred with respect to a classical explicit integration method with a unique integration time step and an uniform grid. Coupling of transient plasma and Navier–Stokes equations for a complete description of plasma–flow interaction should be feasible using AAMR. Work is in progress in this direction. Furthermore AAMR could be applied to various multi-scale transport problems in other fields. Acknowledgements This work was partially supported by the EADS foundation (one of the authors Thomas Unfer is supported by a fellowship of the EADS foundation). References
Fig. 21. Discharge current vs simulated time for various scale parameter L.
done with the length parameter L. L 0 is the reference case and L i +1 = 2L i with the notations of Fig. 20. For the various cases the size of the smallest cell is kept constant to maintain the same level of resolution, whereas twice larger cells may be used for the low density areas in L i +1 with respect to L i . As far as discharge current is concerned, the results are very good. In fact the differences are lower than the noise level, as shown in Fig. 21. Furthermore, AAMR performs impressively in term of CPU time when increasing the scale parameter (see Fig. 22). This test shows the full power of the method. It is another step forward from the asynchronous scheme towards multi-scale multi-physical modeling of sharp localized phenomena. 9. Conclusion We have shown in this paper that the Adaptive Mesh Refinement (AMR) technique described here, coupled with the asynchronous time integration method described in a previous paper, is a very promising technique for multi-scale gas discharge problems. The resulting Asynchronous Adaptive Mesh Refinement (AAMR) method proposed in the present paper can be extremely efficient
[1] J.R. Roth, Aerodynamic flow acceleration using paraelectric and peristaltic electrohydrodynamic effects of a one atmosphere uniform glow discharge plasma, Phys. Plasmas 10 (2003) 2117–2126. [2] J. Pons, E. Moreau, G. Touchard, Asymmetric surface dielectric barrier discharge in air at atmospheric pressure: Electrical properties and induced airflow characteristics, J. Phys. D: Appl. Phys. 38 (2005) 3635–3642. [3] C. Montijn, W. Hundsdorfer, U. Ebert, An adaptive grid refinement strategy for the simulation of negative streamers, J. Comput. Phys. 219 (2006) 801–835. [4] T. Unfer, J.P. Boeuf, F. Rogier, F. Thivet, An asynchronous scheme with local time stepping for multi-scale transport problems: Application to gas discharges, J. Comput. Phys. 227 (2007) 898. [5] M. Berger, J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, J. Comput. Phys. 53 (1984) 484–512. [6] D. Burgarelli, M. Kischinhevsky, R.J. Biezuner, A new adaptive mesh refinement strategy for numerically solving evolutionary PDE’s, J. Comput. Appl. Math. 196 (1) (2006) 115–131. [7] Y. Lagmich, T. Callegari, L. C. Pitchford, J.-P. Boeuf, Model description of surface dielectric barrier discharges. [8] Bolsig+, www.laplace.univ-tlse.fr/-BOLSIG-Resolution-de-l-equation-Valeurs (numériques utilisées dans le modèle). [9] G.J.M. Hagelaar, L.C. Pitchford, Solving the Boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models, Plasma Sources Sci. Technol. 14 (2005) 722–733. [10] H. Samet, The quadtree and related hierarchical data structures, Comput. Surv. 16 (2) (1984). [11] Y. Saad, M. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (3) (1986).